Fortran 90 (floating point error)

[MelihAltunan]

Hello,

I am dealing with the code below and i am receiving 'runtime error M6104:MATH floating point errorverflow'. i searched the error online and it is written the output values have high range than input values that i represent as input. I am using Microsoft Developer Studio and since my computer 64-bit when i debugged the program shut down. So, i couldn't see which lines create problems. If you examine the code and give me any clue about what is the problem, i will be extremely pleased.

The code;

! computer code for lid-driven cavity
parameter (n=100,m=100)
real f(0:8,0:n,0:m)
real rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
integer i
real tminv(0:8,0:8),sm(0:8),tm(0:8,0:8),stmiv(0:8,0:8)
real ev(0:8,0:8)
open(2,file='uvfield')
open(3,file='uvely')
open(4,file='vvelx')
open(8,file='timeu')
open(10,file='tmat')
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
cx(:)=(/0.,1.,0.,-1.,0.,1.,-1.,-1.,1./)
cy(:)=(/0.,0.,1.,0.,-1.,1.,1.,-1.,-1./)
tm(0,:)=(/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/)
tm(1,:)=(/-4.,-1.,-1.,-1.,-1.,-2.0,-2.0,-2.0,-2.0/)
tm(2,:)=(/-4.0,-2.,-2.,-2.,-2.,1.0,1.0,1.0,1.0/)
tm(3,:)=(/0.0,1.0,0.0,-1.,0.0,1.0,-1.,-1.,1.0/)
tm(4,:)=(/0.0,-2.,0.0,2.0,0.0,1.0,-1.,-1.,1.0/)
tm(5,:)=(/0.0,0.0,1.0,0.0,-1.,1.0,1.0,-1.,-1.0/)
tm(6,:)=(/0.0,0.0,-2.,0.0,2.0,1.0,1.0,-1.,-1./)
tm(7,:)=(/0.0,0.0,-1.,1.0,-1.,0.0,0.0,0.0,0.0/)
tm(8,:)=(/0.0,0.0,0.0,0.0,0.0,1.0,-1.,1.0,-1./)
a1=1./36.
tminv(0,:)=(/4.*a1,-4.*a1,4.*a1,0.,0.,0.,0.,0.,0./)
tminv(1,:)=(/4.*a1,-a1,-2.*a1,6.*a1,-6.*a1,0.,0.,9.*a1,0.0/)
tminv(2,:)=(/4.*a1,-a1,-2.*a1,0.,0.0,6.*a1,-6.*a1,-9.*a1,0.0/)
tminv(3,:)=(/4.*a1,-a1,-2.*a1,-6.*a1,6.*a1,0.,0.,9.*a1,0.0/)
tminv(4,:)=(/4.*a1,-a1,-2.*a1,0.,0.,-6.*a1,6.*a1,-9.*a1,0./)
tminv(5,:)=(/4.*a1,2.*a1,a1,6.*a1,3.*a1,6.*a1,3.*a1,0.,9.*a1/)
tminv(6,:)=(/4.*a1,2.*a1,a1,-6.*a1,-3.*a1,6.*a1,3.*a1,0.,-9.*a1/)
tminv(7,:)=(/4.*a1,2.*a1,a1,-6.*a1,-3.*a1,-6.*a1,-3.*a1,0.,9.*a1/)
tminv(8,:)=(/4.*a1,2.*a1,a1,6.*a1,3.*a1,-6.*a1,3.*a1,0.,-9.*a1/)
do i=0,8
do j=0,8
sumcc=0.0
do l=0,8
sumcc=sumcc+tminv(i,l)*tm(l,j)
end do
ev(i,j)=sumcc
end do
end do
do i=0,8
print*,(ev(i,j),j=0,8)
end do

uo=0.05
rhoo=1.00
dx=1.0
dy=dx
dt=1.0
alpha=0.001
Re=uo*m/alpha
print*,"Re=",Re
omega=1.0/(3.*alpha+0.5)
tau=1./omega

sm(:)=(/1.0,1.4,1.4,1.0,1.2,1.0,1.2,tau,tau/)
do i=0,8
do j=0,8
stmiv(i,j)=tminv(i,j)*sm(j)

end do
end do
do i=0,8
print*,(stmiv(i,j),j=0,8)
end do

mstep=1000
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do i=1,n-1
u(i,m)=uo
v(i,m)=0.0
end do
!main loop
do kk=1,mstep
call collesion(u,v,f,rho,n,m,tm,stmiv)
call streaming (f,n,m)
! ---------------
call sfbound(f,n,m,uo)
call rhouv(f,rho,u,v,cx,cy,n,m)
print*,u(0,m/2),v(0,m/2),rho(0,m/2),u(n,m/2),v(n,m/2),rho(n,m/2)
write(8,*)kk,u(n/2,m/2),v(n/2,m/2)
END DO
!end of main loop
call result(u,v,rho,uo,n,m)
stop
end
!end of main program

subroutine collesion(u,v,f,rho,n,m,tm,stmiv)
real f(0:8,0:n,0:m)
real rho(0:n,0:m)
real u(0:n,0:m), v(0:n,0:m)
real tm(0:8,0:8),stmiv(0:8,0:8)
real fmom(0:8,0:n,0:m),fmeq(0:8,0:n,0:m)
!calculate equilibrium moments
do i=0,n
do j=0,m
t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
fmeq(0,i,j)=rho(i,j)
fmeq(1,i,j)=rho(i,j)*(-2.0+3.0*rho(i,j)*t1)
fmeq(2,i,j)=rho(i,j)*(1.0-3.0*rho(i,j)*t1)
fmeq(3,i,j)=rho(i,j)*u(i,j)
fmeq(4,i,j)=-rho(i,j)*u(i,j)
fmeq(5,i,j)=rho(i,j)*v(i,j)
fmeq(6,i,j)=-rho(i,j)*v(i,j)
fmeq(7,i,j)=rho(i,j)*(u(i,j)*u(i,j)-v(i,j)*v(i,j))
fmeq(8,i,j)=rho(i,j)*u(i,j)*v(i,j)
end do
end do
!calculate moments
do i=0,n
do j=0,m
do k=0,8
suma=0.0
do l=0,8
suma=suma+tm(k,l)*f(l,i,j)
end do
fmom(k,i,j)=suma
end do
end do
end do
!collision in the moment space
do i=0,n
do j=0,m
do k=0,8
sumb=0.0
do l=0,8
sumb=sumb+stmiv(k,l)*(fmom(l,i,j)-fmeq(l,i,j))
end do
f(k,i,j)=f(k,i,j)-sumb
end do
end do
end do
return
end
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! streaming
DO j=0,m
DO i=n,1,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=0,n-1 !LEFT TO RIGHT
f(3,i,j)=f(3,i+1,j)
END DO
END DO
DO j=m,1,-1 !TOP TO BOTTOM
DO i=0,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=0,n-1
f(6,i,j)=f(6,i+1,j-1)
END DO
END DO
DO j=0,m-1 !BOTTOM TO TOP
DO i=0,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,1,-1
f(8,i,j)=f(8,i-1,j+1)
END DO
END DO
return
end
subroutine sfbound(f,n,m,uo)
real f(0:8,0:n,0:m)
do j=0,m
! bounce back on west boundary
f(1,0,j)=f(3,0,j)
f(5,0,j)=f(7,0,j)
f(8,0,j)=f(6,0,j)
! bounce back on east boundary
f(3,n,j)=f(1,n,j)
f(7,n,j)=f(5,n,j)
f(6,n,j)=f(8,n,j)
end do
! bounce back on south boundary
do i=0,n
f(2,i,0)=f(4,i,0)
f(5,i,0)=f(7,i,0)
f(6,i,0)=f(8,i,0)
end do
! moving lid, north boundary
do i=1,n-1
rhon=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
f(4,i,m)=f(2,i,m)
f(8,i,m)=f(6,i,m)+rhon*uo/6.0
f(7,i,m)=f(5,i,m)-rhon*uo/6.0
end do
return
end
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m),cx(0:8),cy(0:8)
do j=0,m
do i=0,n
ssum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
end do
end do
do i=1,n
rho(i,m)=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
end do
DO i=1,n
DO j=1,m-1
usum=0.0
vsum=0.0
DO k=0,8
usum=usum+f(k,i,j)*cx(k)
vsum=vsum+f(k,i,j)*cy(k)
END DO
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
END DO
END DO
return
end
subroutine result(u,v,rho,uo,n,m)
real u(0:n,0:m),v(0:n,0:m)
real rho(0:n,0:m),strf(0:n,0:m)
open(5, file='streamf')
! streamfunction calculations
strf(0,0)=0.
do i=0,n
rhoav=0.5*(rho(i-1,0)+rho(i,0))
if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
do j=1,m
rhom=0.5*(rho(i,j)+rho(i,j-1))
strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
end do
end do
!
write(2,*)"VARIABLES =X, Y, U, V, S"
write(2,*)"ZONE ","I=",n+1,"J=",m+1,",","F=BLOCK"
do j=0,m
write(2,*)(i,i=0,n)
end do
do j=0,m
write(2,*)(j,i=0,n)
end do
do j=0,m
write(2,*)(u(i,j),i=0,n)
end do
do j=0,m
write(2,*)(v(i,j),i=0,n)
end do
do j=0,m
write(2,*)(strf(i,j),i=0,n)
end do
do j=0,m
write(3,*)j/float(m),u(n/2,j)/uo,u(3*n/4,j)/uo
end do
do i=0,n
write(4,*) i/float(n),v(i,m/2)/uo
end do
return
end

 

[gsal]

How about you learn to do some basic debugging of your own?

Here is one idea...you can go incremental starting from the beginning or you can take a bisection approach.

Insert a stop statement after a few lines of code; make sure you do not put the stop statement inside if-then-else or do loops...put it after a few relevant lines of code. If the program does not fail but simply stops...move the stop statement a few more relevant lines down; when you get to a subroutine call, move the stop statement beyond just one sub call at a time until the program fails...then, you will know which subroutine the problem is in...then move into the sub and start inserting the stop statement there...

 

[AlephZero]

Insert a stop statement after a few lines of code; ...

You could do all that. Or, you could find out how to display the "stack trace" window in MS Developer, and just read the answer from the screen.

You should be able to look at the values of variables etc after the code has crashed, and work backwards to find what caused the error. (You might have set a variable to the wrong value somewhere before it actually caused the crash).

If you can see that some variables have the wrong values but you can't figure out why the values are wrong, you can set a "watch point" breakpoint to stop the program when a variable is set to that value, and then work backwards from there till you find what caused the error.

 

[MelihAltunan]

As following what you've said, i think i found which lines create problem.
real stmiv(0:8,0:8)
real fmom(0:8,0:n,0:m),fmeq(0:8,0:n,0:m)
do i=0,n
do j=0,m
do k=0,8
sumb=0.0
do l=0,8
sumb=sumb+stmiv(k,l)*(fmom(l,i,j)-fmeq(l,i,j))
end do
f(k,i,j)=f(k,i,j)-sumb
end do
end do
end do

do you have any idea why it is happening?

 

[MelihAltunan]

i did bisectional approach, writing stop statement. Unfournately i am not able to look at the variables after the code has crushed because the program shut down itself(that's why i am looking for compiler workin 64 bit system). Now, i am researching "stack trace", thank you for all your help.

 

[gsal]

If you know how to use the development environment and step through the code and debug...go for it...it is just the I have never had a need to such programs and I feel that, more often than not, they get in the way of my working style.

So, I would continue my advice without the debugger.

First of all, if it is possible, I would reduce the size of the problem from n=100, m=100 to something smaller that still shows the problem.

Right before the nested loops, I would write out the matrices fmom and fmeq and make sure they have reasonable values...if they do...

...then, now that you know which set of do loops the problem is in, I would insert a "write" statement and write out the indices i,j,k and the value of f() and see it progress towards overflow until crash time.

In order to see the output from the write statement, you really need to flush the output queue, otherwise the program will crash without emptying the queue and you will not see what you want...I would do something like:

do i=0,n
    do j=0,m
        do k=0,8
            sumb=0.0
            do l=0,8
                sumb=sumb+stmiv(k,l)*(fmom(l,i,j)-fmeq(l,i,j))
            end do
            f(k,i,j)=f(k,i,j)-sumb
            write(0,*) i, j, k, f(i,j,k)
            flush(0)
        end do
    end do
end do

I am thinking the problem is in f() because it is being manipulated through a lot more loops than sumb, but the error message says "overflow", so maybe the problem is in sumb...then again I am not a computer science kind of guy but maybe a very large negative value also causes "overflow"? or is it then called "underflow"...I don't know...

 

[AlephZero]

If you know how to use the development environment and step through the code and debug...go for it...it is just the I have never had a need to such programs and I feel that, more often than not, they get in the way of my working style.

I'm no great fan of debuggers either. IMO the best way to debug code is by thinking, not playing with a user interface!

But if the OP is already using MS Developer Studio (I assume that's part of the same product line as the free MS Visual Studio) it does no harm to learn how to use it.

Looking at a stack trace "for free" on the screen certainly beats constructing one by hand from 500 printed pages of a memory dump in hexadecimal (which is the way I learned how to debug code!)

 

[MelihAltunan]

I found silverfrost fortran which can work properly on 64 bit system. With using this program debugger, i am able to trace steps. Actually, it is amazing watch the code step by step, you can understand where exactly your error is. Anyway, thank you so much for your comments, i feel like i have lots of great things to learn from you.

 

[Mark44]

I'm no great fan of debuggers either. IMO the best way to debug code is by thinking, not playing with a user interface!

I am a great fan of debuggers. I agree that there is no substitute for thinking about the algorithm before coding it. Although some of the early debuggers were difficult to use and required a fairly steep learning curve, the ones I've used for the past 20 years (Borland and Microsoft) aren't that difficult to get started with, and are very powerful.

I did some x86 assembly work back in the early 90s, and being able to look at the CPU registers and memory was crucial to what I was doing.

But if the OP is already using MS Developer Studio (I assume that's part of the same product line as the free MS Visual Studio) it does no harm to learn how to use it.

Looking at a stack trace "for free" on the screen certainly beats constructing one by hand from 500 printed pages of a memory dump in hexadecimal (which is the way I learned how to debug code!)