0

我正在处理下面的代码,并且收到“运行时错误 M6104:MATH 浮点错误:溢出”。我在网上搜索了错误,并写道输出值的范围比我表示为输入的输入值高。我正在使用 Microsoft Developer Studio,因为当我调试程序时我的计算机是 64 位的,所以我关闭了。所以,我看不出哪些线路会产生问题。如果您检查代码并向我提供有关问题所在的任何线索,我将非常高兴。

编码;

! 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
4

1 回答 1

0

你是如何运行它的?如果您在 Visual Studio 中使用 F5 运行,它应该在代码中断时停止。

小事: a) 尝试为输出通道使用更高的数字。大多数编译器使用 1、2、5、6 之类的低数字作为控制台的输入/来自控制台 b)您似乎没有关闭任何通道 c)在子程序结果中,您打开了流 5,但根本不使用它。d) 如果你缩进你的代码而不是从第 1 列开始,它会更容易阅读。例如

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

e) 通常在 VS 中,浮点异常是关闭的。启动调试器,在第一个可执行行上设置断点,转到 Debug/Exceptions/Win32 exceptions。打开所有浮点异常。f) Visual Studio 中的条件中断可能非常慢,尤其是对于 9x101x101 矩阵。你可以做类似的事情

sumb = sumb + ...
if (sumb .gt. 10000.0) then
   print *,'break here'
end if

并设置打印语句所在的断点。

于 2013-04-13T16:55:45.700 回答