Hi every one
I have written a program for solving Lid Driven Cavity for Re=400 in Fortran. But unfortunately I can get a correct result. please help me to correct all errors in this code.
% c6 c2 c5
% \ | /
% c3 -c0 - c1
% / |
% c7 c4 c8
Module D2Q9
double precision ::w(0:8) =(/4.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0 &
,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0/)
integer:: cx(0:8)=(/0, 1, 0, -1, 0, 1, -1, -1, 1/)
integer:: cy(0:8)=(/0, 0, 1, 0, -1, 1, 1, -1, -1/)
End Module D2Q9
Module simParam
integer,parameter :: n=250
integer,parameter :: m=250
integer,parameter :: Nt=500000
double precision,parameter :: u0=0.2d0
double precision,parameter :: TS=0.d0
double precision,parameter :: TN=1.d0
double precision,parameter :: TW=0.d0
double precision,parameter :: Pr=0.71d0
double precision,parameter :: Re=400.d0
double precision,parameter :: visco=u0*dble(m+1)/Re
double precision,parameter :: alpha=visco/pr
double precision,parameter :: rho0=1.d0
double precision,parameter :: Th0=0.d0
double precision,parameter :: dx=1.d0
double precision,parameter :: dy=1.d0
double precision,parameter :: dt=1.d0
End Module SimParam
program LBMPD2Q9
use D2Q9
use SimParam
implicit none
double precision,dimension(:,:,:),allocatable :: f,feq,g,geq
double precision,dimension(:,:),allocatable:: rho,th,u,v,up,thp,uSqr
double precision:: s1,s2,s3,s4,err1,omega,omegat,snul,snur,rnur,rnul
integer:: i,j,time,k
allocate( u(0:n,0:m) )
allocate( v(0:n,0:m) )
allocate( uSqr(0:n,0:m) )
allocate( Thp(0:n,0:m) )
allocate( up(0:n,0:m) )
allocate( Th(0:n,0:m) )
allocate( rho(0:n,0:m) )
allocate( f(0:8,0:n,0:m) )
allocate( feq(0:8,0:n,0:m) )
allocate( g(0:8,0:n,0:m) )
allocate( geq(0:8,0:n,0:m) )
open(2,file=‘uvfield.data’)
open(3,file=‘uvely1.data’)
open(7,file=‘uvely2.data’)
open(6,file=‘uvely3.data’)
open(8,file=‘uvely4.data’)
open(9,file=‘timeuv.data’)
open(4,file=‘vvelx.data’)
call computeomega(omega,omegat,alpha)
call initialization(rho,u,v,th,f,feq,g,geq,uSqr)
call writeInput(omega,omegat)
! print *, th(0,m/2),rho(0,m/2),u(0,m-1)
! stop
do time=1,Nt
!compute velocity
call computeFeq(feq,rho,u,v,uSqr)
call collisionF(f,feq,omega)
call streaming(f)
call boundF(f,u,rho)
call computeU(rho,f,u,v,uSqr)
!compute temperature
call computeGeq(geq,th,u,v,uSqr)
call collisionG(g,geq,omegat)
call streaming(g)
call boundG(g)
call computeT(th,g)
! print , th(n/2,m/2),rho(0,m/2),u(n,m-1)
write(9,) time,th(n/2,m/2),rho(n/2,m/2),u(n,m/2),v(n,m/2)
s1=0.d0
s2=0.d0
! s3=0.d0
! s4=0.d0
do j=0,m
do i=0,N
s1=s1+abs(th(i,j)-thp(i,j))
s2=s2+abs(thp(i,j))
! s3=s3+abs(v(i,j)-vp(i,j))
! s4=s4+abs(vp(i,j))
enddo
enddo
err1=s1/s2
! err2=s3/s4
if (err1<=1.e-8) then
exit
else
do j=0,m
do i=0,N
thp(i,j)=th(i,j)
! vp(i,j)=v(i,j)
enddo
enddo
endif
if(mod(time,5)==0)then
print*,time,err1 !,err2
endif
!Nusselt number Calculation
snul=0.0
snur=0.0
do i=0,n
rnul=(th(i,0)-th(i,1))*float(n)
rnur=(th(i,m-1)-th(i,m))float(n)
snul=snul+rnul
snur=snur+rnur
!print,i/float(n),snul,snur
end do
enddo
call output(th,u,v,rho)
deallocate( u )
deallocate( v )
deallocate( uSqr )
deallocate( Thp )
deallocate( up )
deallocate( Th )
deallocate( rho )
deallocate( f )
deallocate( feq )
deallocate( g )
deallocate( geq )
end program LBMPD2Q9
! ========================================================
! Print out simulation parameters to screen
! ========================================================
SUBROUTINE writeInput(omega,omegat)
USE simParam
implicit none
double precision, INTENT(IN):: omega,omegat
print*, “Re =”, Re
print*, “visco =”, visco
print*, “alpha =”, alpha
print*, “omega =”, omega
print*, “omegat =”, omegat
return
END SUBROUTINE writeInput
! ========================================================
! compute the relaxation time
! ========================================================
subroutine computeomega(omega,omegat)
use SimParam
implicit none
double precision::omega,omegat
omega =1.0/(3.0d0visco+0.5d0)
omegat=1.0/(3.0d0alpha+0.5d0)
return
end subroutine computeomega
! ========================================================
! initialization
! ========================================================
subroutine initialization(rho,u,v,th,f,feq,g,geq,uSqr)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT):: rho(0:n,0:m) ,u(0:n,0:m) ,v(0:n,0:m) ,f(0:8,0:n,0:m),th(0:n,0:m)
double precision, INTENT(INOUT)::uSqr(0:n,0:m) ,g(0:8,0:n,0:m) ,feq(0:8,0:n,0:m),geq(0:8,0:n,0:m)
double precision::t2,t1
integer::i,j,k
do i=0,n
do j=0,m
rho(i,j)=rho0
u(i,j)=0.d0
v(i,j)=0.d0
th(i,j)=Th0
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)=feq(k,i,j)
geq(k,i,j)=th(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
g(k,i,j)=geq(k,i,j)
enddo
end do
end do
do i=1,n-1
u(i,m)=u0
v(i,m)=0.0
u(i,0)=0.d0
v(i,0)=0.d0
Th(i,m)=TN
th(i,0)=TS
end do
do j=0,m,1
u(0,j)=0.d0
u(N,j)=0.d0
v(0,j)=0.d0
v(N,j)=0.d0
enddo
return
end subroutine initialization
! ========================================================
! Compute equilibrium distribution Feq
! ========================================================
SUBROUTINE computeFeq(feq,rho,u,v,uSqr)
use D2Q9
use SimParam
implicit none
double precision, INTENT(IN)::uSqr(0:n,0:m),rho(0:n,0:m) ,u(0:n,0:m) ,v(0:n,0:m)
double precision, INTENT(INOUT) ::feq(0:8,0:n,0:m)
double precision:: t2,t1
integer::i,j,k
DO i=0,n
DO j=0,m
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)
END DO
END DO
END DO
return
end SUBROUTINE computeFeq
! ========================================================
! Compute equilibrium distribution Geq
! ========================================================
SUBROUTINE computeGeq(geq,th,u,v,uSqr)
use D2Q9
use SimParam
implicit none
double precision, INTENT(IN) ::uSqr(0:n,0:m),th(0:n,0:m) ,u(0:n,0:m) ,v(0:n,0:m)
double precision, INTENT(INOUT) ::geq(0:8,0:n,0:m)
double precision:: t2,t1
integer::i,j,k
DO i=0,n
DO j=0,m
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)
geq(k,i,j)=th(i,j)*w(k)*(1.0 +3.0*t2+4.50*t2*t2-1.50*t1)
END DO
END DO
END DO
return
end SUBROUTINE computeGeq
! ========================================================
! compute collision of the velocity
! ========================================================
subroutine collisionF(f,feq,omega)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT):: f(0:8,0:n,0:m)
double precision, INTENT(IN) :: omega, feq(0:8,0:n,0:m)
integer::i,j,k
DO i=0,n
DO j=0,m
DO k=0,8
f(k,i,j)=f(k,i,j)-omega*(f(k,i,j)-feq(k,i,j))
END DO
END DO
END DO
return
end subroutine collisionF
! ========================================================
! compute collision of the temperature
! ========================================================
subroutine collisionG(g,geq,omegat)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT):: g(0:8,0:n,0:m)
double precision, INTENT(IN) :: geq(0:8,0:n,0:m),omegat
integer::i,j,k
DO i=0,n
DO j=0,m
DO k=0,8
g(k,i,j)=g(k,i,j)-omegat*(g(k,i,j)-geq(k,i,j))
END DO
END DO
END DO
return
end subroutine collisionG
! ========================================================
! compute streaming step
! ========================================================
subroutine streaming(f)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT):: f(0:8,0:n,0:m)
integer::i,j
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 streaming
! ========================================================
! compute boundaries conditions of the velocity
! ========================================================
subroutine boundF(f,u,rho)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT) :: u(0:n,0:m),rho(0:n,0:m)
double precision, INTENT(INOUT) :: f(0:8,0:n,0:m)
double precision:: rhow,Rhoe,rhon
integer::i,j,k
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=1,n-1
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)+rhonu0/6.0
f(7,i,m)=f(5,i,m)-rhonu0/6.0
!f(7,i,m)=f(5,i,m)+6rho(i,m)u0cx(7)w(7)
!f(8,i,m)=f(6,i,m)+6rho(i,m)u0cx(8)*w(8)
end do
return
end subroutine boundF
! ========================================================
! compute boundaries conditions of the temperature
! ========================================================
subroutine boundG(g)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT):: g(0:8,0:n,0:m)
integer::i,j,k
! Boundary conditions
! West boundary condition, T=0.
do j=0,m
g(1,0,j)=-g(3,0,j)
g(5,0,j)=-g(7,0,j)
g(8,0,j)=-g(6,0,j)
end do
! East boundary condition, T=0.
do j=0,m
g(6,n,j)=-g(8,n,j)
g(3,n,j)=-g(1,n,j)
g(7,n,j)=-g(5,n,j)
end do
! Top boundary conditions, T=tw=1.0
do i=0,n
g(2,i,0)=TS*(w(2)+w(4))-g(4,i,0)
g(6,i,0)=TS*(w(6)+w(8))-g(8,i,0)
g(5,i,0)=TS*(w(5)+w(7))-g(7,i,0)
end do
do i=0,n
g(8,i,m)=TN*(w(8)+w(6))-g(6,i,m)
g(7,i,m)=TN*(w(7)+w(5))-g(5,i,m)
g(4,i,m)=TN*(w(4)+w(2))-g(2,i,m)
end do
return
end subroutine boundG
! ========================================================
! compute velocity
! ========================================================
subroutine computeU(rho,f,u,v,uSqr)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT)::rho(0:n,0:m) ,u(0:n,0:m) ,v(0:n,0:m),uSqr(0:n,0:m)
double precision, INTENT(IN) ::f(0:8,0:n,0:m)
double precision::sum,usum,vsum
integer::i,j,k
do i=0,n
do j=0,m
! rho(i,j)=f(0,i,j)+f(1,i,j)+f(2,i,j)+f(3,i,j)+f(4,i,j)+f(5,i,j)+f(6,i,j)+f(7,i,j)+f(8,i,j)
! u(i,j) = (f(1,i,j)+ f(5,i,j)+f(8,i,j) - f(3,i,j)- f(6,i,j) -f(7,i,j)) / rho(i,j)
! v(i,j) = (f(2,i,j)+ f(5,i,j)+f(6,i,j) - f(4,i,j) - f(7,i,j)-f(8,i,j)) / rho(i,j)
! uSqr(i,j)=u(i,j)*u(i,j)+v(i,j)*v(i,j)
end do
end do
do j=1,m-1
do i=1,n-1
sum=0.0
do k=0,8
sum=sum+f(k,i,j)
end do
rho(i,j)=sum
end do
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 computeU
! ========================================================
! compute temperature
! ========================================================
subroutine computeT(th,g)
use D2Q9
use SimParam
implicit none
double precision, INTENT(INOUT):: th(0:n,0:m)
double precision, INTENT(IN) ::g(0:8,0:n,0:m)
double precision:: sumt
integer::i,j,k
! do i=0,n
! do j=0,m
! th(i,j)=g(0,i,j)+g(1,i,j)+g(2,i,j)+g(3,i,j)+g(4,i,j)+g(5,i,j)+g(6,i,j)+g(7,i,j)+g(8,i,j)
! end do
! end do
do j=1,m-1
do i=1,n-1
sumt=0.0
do k=0,8
sumt=sumt+g(k,i,j)
end do
th(i,j)=sumt
end do
end do
return
end subroutine computeT
! ========================================================
! compute the result
! ========================================================
subroutine output(th,u,v,rho)
use D2Q9
use SimParam
implicit none
double precision:: th(0:n,0:m) ,u(0:n,0:m) ,v(0:n,0:m),rho(0:n,0:m)
integer::i,j
write(2,)“VARIABLES =X, Y, U, V, T”
write(2,)“ZONE “,“I=”,n+1,“J=”,m+1,”,”,“F=Block”
do j=0,m
write(2,)(i/float(n),i=0,n)
end do
do j=0,m
write(2,)(j/float(n),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,)(th(i,j),i=0,n)
end do
do j=0,m
write(3,)j/float(m),u(n-30,j)/u0
write(7,)j/float(m),u(n-20,j)/u0
write(6,)j/float(m),u(n-10,j)/u0
write(8,)j/float(m),u(n/2,j)/u0
end do
do i=0,n
write(4,) i/float(n),v(i,m/2)/u0
end do
return
end subroutine output