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

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

各位大佬，最近刚学lbm，利用Fortran进行编程，但是代码输出数据全部是零，不知道该怎么改了，特来求教，请解答，以下附属代码

[Fortran] 纯文本查看 复制代码

001parameter (n=1000,m=40)
002real f(0:8,0:n,0:m)
003real feq(0:8,0:n,0:m),rho(0:n,0:m)
004real w(0:8), cx(0:8),cy(0:8)
005real u(0:n,0:m), v(0:n,0:m)
006real ssum, usum, vsum
007integer i
008open(2,file='uvfield')
009open(3,file='uvely')
010open(4,file='vvelx')
011open(8,file='timeu')
012uo=0.20
013sumvelo=0.0
014rhoo=5.00 ! 虽然不知道为什么在后边密度的计算没有加上rhoo
015dx=1.0
016dy=dx
017dt=1.0
018alpha=0.02
019Re=uo*m/alpha
020print *, "Re=", Re
021omega=1.0/(3.*alpha+0.5)
022mstep=40000      ! 时间步长的确定需要考虑
023cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
024cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
025w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
026! 初始条件的确定  对于初始条件的设定，是流体没有进入空间之前的设定
027do j=21,m
028do i=0,40
029rho(i,j)=rhoo
030u(i,j)=0.0
031v(i,j)=0.0
032end do
033    end do
034do j=0,m
035    do i=41,n
036        rho(i,j)=rhoo
037        u(i,j)=0.0
038        v(i,j)=0.0
039    end do
040end do
041! 边界条件设定   这个是流体没有进入前，也就是流体状态的设定
042do j=21,m
043    u(0,j)=uo
044    v(0,j)=0.0
045    end do
046! main loop
0471 do kk=1,mstep
048call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
049call streaming(f,n,m)
050! ——————————–
051call sfbound(f,n,m,uo)
052call rhouv(f,rho,u,v,cx,cy,n,m)
053!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)
054write(8,*) kk,u(n/2,m/2),v(n/2,m/2)  
055    end do
056! end of the main loop
057call result(u,v,rho,uo,n,m)
058stop
059end 
060subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
061real f(0:8,0:n,0:m)
062real feq(0:8,0:n,0:m),rho(0:n,0:m)
063real w(0:8), cx(0:8),cy(0:8)
064real u(0:n,0:m), v(0:n,0:m)
065do j=21,m    ! j=20是代表障碍物，21代表了流体
066    do i=0,40
067        t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
068DO k=0,8
069t2=u(i,j)*cx(k)+v(i,j)*cy(k)
070feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
071f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
072END DO
073END DO
074END DO
075do j=0,m
076    do i=41,n
077        t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
078DO k=0,8
079t2=u(i,j)*cx(k)+v(i,j)*cy(k)
080feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
081f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
082END DO
083END DO
084END DO
085return
086    end
087subroutine streaming(f,n,m)
088real f(0:8,0:n,0:m)
089! 没有障碍物的流情况 这一部分存在着命名的重叠，先写后改 这一部分的节点是有问题的
090! 存在障碍物的部分流
091do j=22,m
092    do i=40,1,-1
093        f(1,i,j)=f(1,i-1,j)
094    end do
095    do i=0,39
096        f(3,i,j)=f(3,i+1,j)
097    end do
098end do
099do j=40,22,-1
100    do i=0,40
101        f(2,i,j)=f(2,i,j-1)
102    end do
103    do i=40,1,-1
104        f(5,i,j)=f(5,i-1,j-1)
105    end do
106    do i=0,39
107        f(6,i,j)=f(6,i+1,j-1)
108    end do
109end do
110do j=21,39
111    do i=0,40
112        f(4,i,j)=f(4,i,j+1)
113    end do
114    do i=0,39
115        f(7,i,j)=f(7,i+1,j+1)
116    end do
117    do i=40,1,-1
118        f(8,i,j)=f(8,i-1,j+1)
119    end do
120end do
121DO j=0,m
122DO i=n,42,-1 !RIGHT TO LEFT
123f(1,i,j)=f(1,i-1,j)
124END DO
125DO i=41,n-1 !LEFT TO RIGHT
126f(3,i,j)=f(3,i+1,j)
127END DO
128END DO
129DO j=m,1,-1 !TOP TO BOTTOM
130DO i=41,n
131f(2,i,j)=f(2,i,j-1)
132END DO
133DO i=n,42,-1
134f(5,i,j)=f(5,i-1,j-1)
135END DO
136DO i=41,n-1
137f(6,i,j)=f(6,i+1,j-1)
138END DO
139END DO
140DO j=0,m-1 !BOTTOM TO TOP
141DO i=41,n
142f(4,i,j)=f(4,i,j+1)
143END DO
144DO i=41,n-1
145f(7,i,j)=f(7,i+1,j+1)
146END DO
147DO i=n,42,-1
148f(8,i,j)=f(8,i-1,j+1)
149END DO
150END DO
151return
152    end
153subroutine sfbound(f,n,m,uo)
154real f(0:8,0:n,0:m)
155real uo,rhow
156! 有障碍物和没有障碍物的底边条件    这一部分没验证输出数据
157do i=0,n
158    f(2,i,0)=f(4,i,0)
159    f(5,i,0)=f(7,i,0)
160    f(6,i,0)=f(8,i,0)
161    end do
162! 存在障碍物北边界条件
163do i=0,40
164    f(2,i,20)=f(4,i,20)
165    f(5,i,20)=f(7,i,20)
166    f(6,i,20)=f(8,i,20)
167end do
168! 存在障碍物右边界条件
169do j=0,20
170    f(1,40,j)=f(3,40,j)
171    f(5,40,j)=f(7,40,j)
172    f(8,40,j)=f(6,40,j)
173end do
174! 左边没有障碍物的已知速度边界条件
175do j=21,m
176    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)
177    f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.
178    f(5,0,j)=f(7,0,j)+rhow*uo/6.
179    f(8,0,j)=f(6,0,j)+rhow*uo/6.
180end do
181! 上方边界
182do i=0,n
183    f(4,i,m)=f(2,i,m)
184    f(7,i,m)=f(5,i,m)
185    f(8,i,m)=f(6,i,m)
186end do
187! 没有障碍物的右边界
188do j=1,m
189    f(1,n,j)=2.*f(1,n-1,j)-f(1,n-2,j)
190    f(5,n,j)=2.*f(5,n-1,j)-f(5,n-2,j)
191    f(8,n,j)=2.*f(8,n-1,j)-f(8,n-2,j)
192end do
193return
194end
195subroutine rhouv(f,rho,u,v,cx,cy,n,m)
196real 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)
197integer :: i, j, k
198real :: ssum, usum, vsum
199! 有障碍物的一部分
200do j=21,m
201    do i=0,40
202        ssum=0.0
203            do k=0,8
204            ssum=ssum+f(k,i,j)
205            end do
206            rho(i,j)=ssum    
207            usum=0.0
208            vsum=0.0
209            do k=0,8
210            usum=usum+f(k,i,j)*cx(k)
211            vsum=vsum+f(k,i,j)*cy(k) 
212            !print*,usum! 为什么？？？？？？
213            end do
214            u(i,j)=usum/rho(i,j)
215            v(i,j)=vsum/rho(i,j)
216            !print*,u(i,j),v(i,j)
217        end do
218end do
219! 这一部分错了
220do j=0,m
221    do i=41,n
222        ssum=0.0
223        usum=0.0
224        vsum=0.0
225        do k=0,8
226            ssum=ssum+f(k,i,j)
227        end do
228        rho(i,j)=ssum  
229        do k=0,8
230            usum=usum+f(k,i,j)*cx(k)
231            vsum=vsum+f(k,i,j)*cy(k)
232        end do
233        u(i,j)=usum/rho(i,j)
234        v(i,j)=vsum/rho(i,j)
235    end do
236end do
237do j=1,m
238    v(n,j)=0.0
239end do
240do j=0,20
241    do i=0,40
242        u(i,j)=0.0
243        v(i,j)=0.0
244    end do
245end do
246return
247    end
248subroutine result(u,v,rho,uo,n,m) 
249! 流函数这一部分我太会
250   real u(0:n,0:m),v(0:n,0:m),rho(0:n,0:m),strf(0:n,0:m)
251    integer :: i, j
252    real :: rhoav, rhom
253    open(5, file='streamf')
254    strf(0,0)=0.
255    do i=1,n
256        rhoav=0.5*(rho(i-1,0)+rho(i,0))
257        strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
258    end do
259    do i=0,40
260        do j=21,m
261            rhom=0.5*(rho(i,j)+rho(i,j-1))
262            strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
263        end do
264    end do
265    do i=41,n
266        do j=1,m
267            rhom=0.5*(rho(i,j)+rho(i,j-1))
268            strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
269        end do
270    end do
271    ! 写入数据到输出文件
272    write(2,*)"VARIABLES =X, Y, U, V, S"
273    write(2,*)"ZONE ","I=",n+1,"J=",m+1,",","F=BLOCK"
274    do j=0,m
275        do i=0,n
276            write(2,*)i,j,u(i,j),v(i,j),strf(i,j)
277        end do
278    end do
279    do j=0,m
280        write(3,*)j/float(m),u(n/2,j)/uo,u(n/4,j)/uo,u(3*n/4,j)/uo
281    end do
282    do i=0,n
283        write(4,*) i/float(n),v(i,m/2)/uo
284    end do
285end 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，如果是逻辑问题只能你自己找原因。
为什么不找个开源库， ...

谢谢您给出的建议，在代码中

[Fortran] 纯文本查看 复制代码

01do j=21,m
02    do i=0,40
03        ssum=0.0
04            do k=0,8
05            ssum=ssum+f(k,i,j)
06            end do
07            rho(i,j)=ssum 
08            usum=0.0
09            vsum=0.0
10            do k=0,8
11            usum=usum+f(k,i,j)*cx(k)
12            vsum=vsum+f(k,i,j)*cy(k) 
13            !print*,usum! 为什么？？？？？？
14            end do
15            !print*,usum !这个可以正常输出数据
16            u(i,j)=usum/rho(i,j)
17            v(i,j)=vsum/rho(i,j)
18            !print*,u(i,j),v(i,j)
19        end do
20end do
21! 这一部分错了
22do j=0,m
23    do i=41,n
24        ssum=0.0
25        do k=0,8
26            ssum=ssum+f(k,i,j)
27        end do  
28        rho(i,j)=ssum  
29        usum=0.0
30        vsum=0.0
31        do k=0,8
32            usum=usum+f(k,i,j)*cx(k)
33            vsum=vsum+f(k,i,j)*cy(k)   
34            !print*,usum
35        end do
36        !print*,usum
37        u(i,j)=usum/rho(i,j)
38        v(i,j)=vsum/rho(i,j)
39    end do
40end do

后一部分循环usum结果为零，不知道是什么原因。
代码本身是根据不同条件改写别人已有的代码，再次感谢您的建议

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

已解决

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

规矩勋章