Jorgen
2007-08-31, 06:33
A lot of people that are starting to write user-material subroutines have experience with computer languages such as c, c++, or java, but not much experience with Fortran. This can cause problems since Fortran does not always behave intuitively. One common problem relates to the handling of floating point numbers.
First, however, I want to point out that although most (all?) commercial FE codes are written in Fortran77 and the User's Manuals only talk about Fortran77, it is in my opinion better to use Fortran90. Which is more modern, more pleasant to work with, and backwards compatible with Fortran77.
Now, here's my example about floating point numbers. Suppose that you want calculate the exponential of the floating point number 0.2.
A quick check with Matlab gives:
>> format long
>> exp(-0.2)
ans =
0.818730753077982
One way to do this with Fortran is as follows:
program test
implicit none
real :: z, res
z = -0.2
res = exp(z)
print *, "res = ", res
end program test
Compiling and running this example gives:
% ifort -o test6 test6.f90 ; ./test6
res = 0.8187308
which is a reasonable single precision approximation.
If I want double precision then I use use the following code:
program test
implicit none
double precision :: z, res
z = -0.2D0
res = dexp(z)
print *, "res = ", res
end program test
which when run gives
% ifort -o test6 test6.f90 ; ./test6
res = 0.818730753077982
which gives exactly the same result as Matlab.
Now, what if you forgot to add 'D0' to the floating point variable z:
program test
implicit none
double precision :: z, res
z = -0.2
res = dexp(z)
print *, "res = ", res
end program test
Then you get
% ifort -o test6 test6.f90 ; ./test6
res = 0.818730750637974
which is not the answer you want!
OK, so far so good. But what if you are writing an ABAQUS VUMAT and you want to create both single precision and double precision versions of your code. Do you really need to create two difference source files? That sounds pretty dumb, doesn't it. One way to resolve that problem, and this is the approach that I recommend in general, is to write the code as follows:
program test
implicit none
real :: z, res
z = -0.2
res = exp(z)
print *, "res = ", res
end program test
and then simply tell the compiler if you want single precision or double precision.
That is, single precision:
% ifort -o test6 test6.f90 ; ./test6
res = 0.8187308
and double precision:
% ifort -r8 -o test6 test6.f90 ; ./test6
res = 0.818730753077982
First, however, I want to point out that although most (all?) commercial FE codes are written in Fortran77 and the User's Manuals only talk about Fortran77, it is in my opinion better to use Fortran90. Which is more modern, more pleasant to work with, and backwards compatible with Fortran77.
Now, here's my example about floating point numbers. Suppose that you want calculate the exponential of the floating point number 0.2.
A quick check with Matlab gives:
>> format long
>> exp(-0.2)
ans =
0.818730753077982
One way to do this with Fortran is as follows:
program test
implicit none
real :: z, res
z = -0.2
res = exp(z)
print *, "res = ", res
end program test
Compiling and running this example gives:
% ifort -o test6 test6.f90 ; ./test6
res = 0.8187308
which is a reasonable single precision approximation.
If I want double precision then I use use the following code:
program test
implicit none
double precision :: z, res
z = -0.2D0
res = dexp(z)
print *, "res = ", res
end program test
which when run gives
% ifort -o test6 test6.f90 ; ./test6
res = 0.818730753077982
which gives exactly the same result as Matlab.
Now, what if you forgot to add 'D0' to the floating point variable z:
program test
implicit none
double precision :: z, res
z = -0.2
res = dexp(z)
print *, "res = ", res
end program test
Then you get
% ifort -o test6 test6.f90 ; ./test6
res = 0.818730750637974
which is not the answer you want!
OK, so far so good. But what if you are writing an ABAQUS VUMAT and you want to create both single precision and double precision versions of your code. Do you really need to create two difference source files? That sounds pretty dumb, doesn't it. One way to resolve that problem, and this is the approach that I recommend in general, is to write the code as follows:
program test
implicit none
real :: z, res
z = -0.2
res = exp(z)
print *, "res = ", res
end program test
and then simply tell the compiler if you want single precision or double precision.
That is, single precision:
% ifort -o test6 test6.f90 ; ./test6
res = 0.8187308
and double precision:
% ifort -r8 -o test6 test6.f90 ; ./test6
res = 0.818730753077982