各位大佬，最近刚学lbm，利用Fortran进行编程，但是代码输...-Fortran语法讨论-专业Fortran论坛

小贺发表于 2023-12-3 22:14:06

各位大佬，最近刚学lbm，利用Fortran进行编程，但是代码输...

各位大佬，最近刚学lbm，利用Fortran进行编程，但是代码输出数据全部是零，不知道该怎么改了，特来求教，请解答，以下附属代码
parameter (n=1000,m=40)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),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)
real ssum, usum, vsum
integer i
open(2,file='uvfield')
open(3,file='uvely')
open(4,file='vvelx')
open(8,file='timeu')
uo=0.20
sumvelo=0.0
rhoo=5.00 ! 虽然不知道为什么在后边密度的计算没有加上rhoo
dx=1.0
dy=dx
dt=1.0
alpha=0.02
Re=uo*m/alpha
print *, "Re=", Re
omega=1.0/(3.*alpha+0.5)
mstep=40000    ! 时间步长的确定需要考虑
cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
! 初始条件的确定对于初始条件的设定，是流体没有进入空间之前的设定
do j=21,m
do i=0,40
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do j=0,m
do i=41,n
   rho(i,j)=rhoo
   u(i,j)=0.0
   v(i,j)=0.0
end do
end do
! 边界条件设定这个是流体没有进入前，也就是流体状态的设定
do j=21,m
u(0,j)=uo
v(0,j)=0.0
end do
! main loop
1 do kk=1,mstep
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
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 the main loop
call result(u,v,rho,uo,n,m)
stop
end
subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),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)
do j=21,m ! j=20是代表障碍物，21代表了流体
do i=0,40
   t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO
END DO
END DO
do j=0,m
do i=41,n
   t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO
END DO
END DO
return
end
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! 没有障碍物的流情况这一部分存在着命名的重叠，先写后改这一部分的节点是有问题的
! 存在障碍物的部分流
do j=22,m
do i=40,1,-1
   f(1,i,j)=f(1,i-1,j)
end do
do i=0,39
   f(3,i,j)=f(3,i+1,j)
end do
end do
do j=40,22,-1
do i=0,40
   f(2,i,j)=f(2,i,j-1)
end do
do i=40,1,-1
   f(5,i,j)=f(5,i-1,j-1)
end do
do i=0,39
   f(6,i,j)=f(6,i+1,j-1)
end do
end do
do j=21,39
do i=0,40
   f(4,i,j)=f(4,i,j+1)
end do
do i=0,39
   f(7,i,j)=f(7,i+1,j+1)
end do
do i=40,1,-1
   f(8,i,j)=f(8,i-1,j+1)
end do
end do
DO j=0,m
DO i=n,42,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=41,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=41,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,42,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=41,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=41,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=41,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,42,-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)
real uo,rhow
! 有障碍物和没有障碍物的底边条件这一部分没验证输出数据
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
! 存在障碍物北边界条件
do i=0,40
f(2,i,20)=f(4,i,20)
f(5,i,20)=f(7,i,20)
f(6,i,20)=f(8,i,20)
end do
! 存在障碍物右边界条件
do j=0,20
f(1,40,j)=f(3,40,j)
f(5,40,j)=f(7,40,j)
f(8,40,j)=f(6,40,j)
end do
! 左边没有障碍物的已知速度边界条件
do j=21,m
rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1.-uo)
f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.
f(5,0,j)=f(7,0,j)+rhow*uo/6.
f(8,0,j)=f(6,0,j)+rhow*uo/6.
end do
! 上方边界
do i=0,n
f(4,i,m)=f(2,i,m)
f(7,i,m)=f(5,i,m)
f(8,i,m)=f(6,i,m)
end do
! 没有障碍物的右边界
do j=1,m
f(1,n,j)=2.*f(1,n-1,j)-f(1,n-2,j)
f(5,n,j)=2.*f(5,n-1,j)-f(5,n-2,j)
f(8,n,j)=2.*f(8,n-1,j)-f(8,n-2,j)
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)
integer :: i, j, k
real :: ssum, usum, vsum
! 有障碍物的一部分
do j=21,m
do i=0,40
   ssum=0.0
         do k=0,8
         ssum=ssum+f(k,i,j)
         end do
         rho(i,j)=ssum
         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)
         !print*,usum! 为什么？？？？？？
         end do
         u(i,j)=usum/rho(i,j)
         v(i,j)=vsum/rho(i,j)
         !print*,u(i,j),v(i,j)
   end do
end do
! 这一部分错了
do j=0,m
do i=41,n
   ssum=0.0
   usum=0.0
   vsum=0.0
   do k=0,8
         ssum=ssum+f(k,i,j)
   end do
   rho(i,j)=ssum
   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
do j=1,m
v(n,j)=0.0
end do
do j=0,20
do i=0,40
   u(i,j)=0.0
   v(i,j)=0.0
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),rho(0:n,0:m),strf(0:n,0:m)
integer :: i, j
real :: rhoav, rhom
open(5, file='streamf')
strf(0,0)=0.
do i=1,n
   rhoav=0.5*(rho(i-1,0)+rho(i,0))
   strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
end do
do i=0,40
   do j=21,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
do i=41,n
   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
   do i=0,n
         write(2,*)i,j,u(i,j),v(i,j),strf(i,j)
   end do
end do
do j=0,m
   write(3,*)j/float(m),u(n/2,j)/uo,u(n/4,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
end subroutine result

necrohan 发表于 2023-12-4 13:22:07

我按你的计算3个时间步，在(0,21)那里数据不是0，如果是逻辑问题只能你自己找原因。
为什么不找个开源库，在那个基础上做。

小贺发表于 2023-12-4 15:06:42

necrohan 发表于 2023-12-4 13:22
我按你的计算3个时间步，在(0,21)那里数据不是0，如果是逻辑问题只能你自己找原因。
为什么不找个开源库， ...
谢谢您给出的建议，在代码中
do j=21,m
do i=0,40
   ssum=0.0
         do k=0,8
         ssum=ssum+f(k,i,j)
         end do
         rho(i,j)=ssum
         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)
         !print*,usum! 为什么？？？？？？
         end do
         !print*,usum !这个可以正常输出数据
         u(i,j)=usum/rho(i,j)
         v(i,j)=vsum/rho(i,j)
         !print*,u(i,j),v(i,j)
   end do
end do
! 这一部分错了
do j=0,m
do i=41,n
   ssum=0.0
   do k=0,8
         ssum=ssum+f(k,i,j)
   end do
   rho(i,j)=ssum
   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)
         !print*,usum
   end do
   !print*,usum
   u(i,j)=usum/rho(i,j)
   v(i,j)=vsum/rho(i,j)
end do
end do
后一部分循环usum结果为零，不知道是什么原因。
代码本身是根据不同条件改写别人已有的代码，再次感谢您的建议

小贺发表于 2023-12-6 19:28:15

已解决

页: [1]

Fortran Coder's Archiver

各位大佬，最近刚学lbm，利用Fortran进行编程，但是代码输...