3D Cavity implementation issues

Ashok · June 24, 2014, 2:33pm

Hello everyone,

I have written a 3D cavity code which uses D3Q19 lattice model. I have some difficulties in debugging it still. Will be very helpful if someone notices the error. Thanks in advance.

[code=“fortran”]
module var
implicit none
save

integer, parameter :: sp = kind(1.0)
integer, parameter :: dp = kind(1.D0)

integer :: xmin, xmax, ymin, ymax, zmin, zmax
integer :: MaxSteps, iStep

real(kind = dp) :: error, tau, Re, Uo, eps, errr, tu

! Distributions weights (D3Q19: W0 = 1/3, W1 = 1/18, W2 = 1/36)
      real(kind = dp), parameter :: W0 = 0.33333333333333333333D0
  real(kind = dp), parameter :: W1 = 0.05555555555555555556D0
  real(kind = dp), parameter :: W2 = 0.02777777777777777778D0

   ! Modified distribution weight (to avoid operations: WiC = Wi*invCs_sq)
 real(kind = dp), parameter :: W0C = 1.00000000000000000000D0
 real(kind = dp), parameter :: W1C = 0.16666666666666666667D0
 real(kind = dp), parameter :: W2C = 0.08333333333333333333D0

real(kind = dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: f, feq, ftemp

real(kind = dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: u, ut
real(kind = dp), ALLOCATABLE, DIMENSION(:,:,:) :: rho

real(kind = dp), dimension(0:18) :: cx, cy, cz

real(kind = dp), parameter :: Cs       = 0.57735026918962576451D0
real(kind = dp), parameter :: Cs_sq    = 0.33333333333333333333D0
real(kind = dp), parameter :: invCs_sq = 3.00000000000000000000D0

   end module var

!------------------------------------------------------------------------------------------------------------------
program cavity
use var
implicit none

	call parameters
	call memalloc
	call initial

!------------------Begining of time steps-------------------
do iStep = 1, MaxSteps

	    call streaming
	    call boundary
	    call hydrodynamics
	    call collision
	    call stats

	    if(errr.lt.eps ) then
	       call output
	       stop
	    end if

	end do
	!	call output

	end program cavity

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

	subroutine parameters
	use var
	implicit none


	xmin = 1
	xmax = 41
	ymin = 1
	ymax = 41
	zmin = 1
	zmax = 41

	MaxSteps = 80001

	Uo = 0.1
	Re = 100.D0

	tu = (3.D0*Uo*dble(xmax - 1))/Re
	tau = tu + 0.5D0

	print*,'tau =',tau

	eps = 1e-06

	end subroutine parameters

!----------------------------------------------------------------------------------

	subroutine MemAlloc
	use var
            implicit none
	
	allocate( rho(xmin:xmax,ymin:ymax,zmin:zmax) )		

	allocate(  u(xmin:xmax,ymin:ymax,zmin:zmax,1:3) )
	allocate( ut(xmin:xmax,ymin:ymax,zmin:zmax,1:3) )

	allocate(     f(xmin:xmax,ymin:ymax,zmin:zmax,0:18) )
	allocate(   feq(xmin:xmax,ymin:ymax,zmin:zmax,0:18) )		
	allocate( ftemp(xmin:xmax,ymin:ymax,zmin:zmax,0:18) )

	end subroutine MemAlloc

!----------------------------------------------------------------------------------------------------------------

	subroutine initial
	use var
	implicit none

! Local variables
integer :: i,j,k,n
real(kind = dp) :: ux, uy, uz, rhon, U_sq
real(kind = dp) :: v

	cx = (/0,1,-1,0,0,0,0,1,1,-1,-1, &
	       1,-1,1,-1,0,0,0,0/)

	cy = (/0,0,0,1,-1,0,0,1,-1,1,-1, &
	       0,0,0,0,1,1,-1,-1/)

	cz = (/0,0,0,0,0,1,-1,0,0,0,0,   &
	       1,1,-1,-1,1,-1,1,-1/)	 

	do i = xmin,xmax
	do j = ymin,ymax
	do k = zmin,zmax
	
	   u(i,j,k,3) = 0.D0
	   u(i,j,k,2) = 0.D0

	 if(k.ne.zmax) then
	   u(i,j,k,1) = 0.D0
	 else
	   u(i,j,k,1) = Uo
	 end if

	   rho(i,j,k) = 1.D0

	   ux   = u(i,j,k,1)
	   uy   = u(i,j,k,2)
	   uz   = u(i,j,k,3)
	   U_sq = 0.5D0*( ux*ux + uy*uy + uz*uz )
	   rhon = rho(i,j,k)		
	 
	   v = ux*cx(0) + uy*cy(0) + uz*cz(0)
	   feq(i,j,k, 0) = rhon*(W0 + W0C*( v + 1.5D0*v*v - U_sq ) )

	   do n = 1,6		
	      v = ux*cx(n) + uy*cy(n) + uz*cz(n)
	      feq(i,j,k,n) = rhon*(W1 + W1C*( v + 1.5D0*v*v - U_sq ) )
	   end do

	   do n = 7,18		
	      v = ux*cx(n) + uy*cy(n) + uz*cz(n)
	      feq(i,j,k,n) = rhon*(W2 + W2C*( v + 1.5D0*v*v - U_sq ) )
	   end do		   
	   
	   do n = 0,18
	      f(i,j,k,n) = feq(i,j,k,n) 
	   end do

	end do
	end do
	end do


	end subroutine initial

!-------------------------------------------------------------------------------------------------------------------------

	subroutine streaming
	use var
	implicit none

! Local variables
integer:: i, j, k, n
!---------------------------------------------------
!------------------streaming------------------------
do i = xmin,xmax
do j = ymin,ymax
do k = zmin,zmax
do n = 0,18
ftemp(i,j,k,n) = f(i,j,k,n)
end do
end do
end do
end do
!---------------------------------------------------
! Direction 0----[0,0,0]
do i = xmin,xmax
do j = ymin,ymax
do k = zmin,zmax-1
f(i,j,k,0) = ftemp(i,j,k,0)
end do
end do
end do

! Direction 1-----[1,0,0]
do i = xmin+1,xmax
do j = ymin,ymax
do k = zmin,zmax-1
f(i,j,k,1) = ftemp(i-1,j,k,1)
end do
end do
end do

! Direction 2-----[-1,0,0]
do i = xmax-1,xmin,-1
do j = ymin,ymax
do k = zmin,zmax-1
f(i,j,k,2) = ftemp(i+1,j,k,2)
end do
end do
end do

! Direction 3-----[0,1,0]
do i = xmin,xmax
do j = ymin+1,ymax
do k = zmin,zmax-1
f(i,j,k,3) = ftemp(i,j-1,k,3)
end do
end do
end do

! Direction 4-----[0,-1,0]
do i = xmin,xmax
do j = ymax-1,ymin,-1
do k = zmin,zmax-1
f(i,j,k,4) = ftemp(i,j+1,k,4)
end do
end do
end do

! Direction 5-----[0,0,1]
do i = xmin,xmax
do j = ymin,ymax
do k = zmin+1,zmax-1
f(i,j,k,5) = f(i,j,k-1,5)
end do
end do
end do

! Direction 6-----[0,0,-1]
do i = xmin,xmax
do j = ymin,ymax
do k = zmax-1,zmin,-1
f(i,j,k,6) = ftemp(i,j,k+1,6)
end do
end do
end do

! Direction 7-----[1,1,0]
do i = xmin+1,xmax
do j = ymin+1,ymax
do k = zmin,zmax-1
f(i,j,k,7) = f(i-1,j-1,k,7)
end do
end do
end do

! Direction 8-----[1,-1,0]
do i = xmin+1,xmax
do j = ymax-1,ymin,-1
do k = zmin,zmax-1
f(i,j,k,8) = ftemp(i-1,j+1,k,8)
end do
end do
end do

! Direction 9-----[-1,1,0]
do i = xmax-1,xmin,-1
do j = ymin+1,ymax
do k = zmin,zmax-1
f(i,j,k,9) = ftemp(i+1,j-1,k,9)
end do
end do
end do

! Direction 10-----[-1,-1,0]
do i = xmax-1,xmin,-1
do j = ymax-1,ymin,-1
do k = zmin,zmax-1
f(i,j,k,10) = ftemp(i+1,j+1,k,10)
end do
end do
end do

! Direction 11-----[1,0,1]
do i = xmin+1,xmax
do j = ymin,ymax
do k = zmin+1,zmax-1
f(i,j,k,11) = ftemp(i-1,j,k-1,11)
end do
end do
end do

! Direction 12-----[-1,0,1]
do i = xmax-1,xmin,-1
do j = ymin,ymax
do k = zmin+1,zmax-1
f(i,j,k,12) = ftemp(i+1,j,k-1,12)
end do
end do
end do

! Direction 13-----[1,0,-1]
do i = xmin+1,xmax
do j = ymin,ymax
do k = zmax-1,zmin,-1
f(i,j,k,13) = ftemp(i-1,j,k+1,13)
end do
end do
end do

! Direction 14-----[-1,0,-1]
do i = xmax-1,xmin,-1
do j = ymin,ymax
do k = zmax-1,zmin,-1
f(i,j,k,14) = ftemp(i+1,j,k+1,14)
end do
end do
end do

! Direction 15-----[0,1,1]
do i = xmin,xmax
do j = ymin+1,ymax
do k = zmin+1,zmax-1
f(i,j,k,15) = ftemp(i,j-1,k-1,15)
end do
end do
end do

! Direction 16-----[0,1,-1]
do i = xmin,xmax
do j = ymin+1,ymax
do k = zmax-1,zmin,-1
f(i,j,k,16) = ftemp(i,j-1,k+1,16)
end do
end do
end do

! Direction 17-----[0,-1,1]
do i = xmin,xmax
do j = ymax-1,ymin,-1
do k = zmin+1,zmax-1
f(i,j,k,17) = ftemp(i,j+1,k-1,17)
end do
end do
end do

! Direction 18-----[0,-1,-1]
do i = xmin,xmax
do j = ymax-1,ymin,-1
do k = zmax-1,zmin,-1
f(i,j,k,18) = ftemp(i,j+1,k+1,18)
end do
end do
end do

	return
	end subroutine streaming

!---------------------------------------------------------------------------------------------------------------------------

	subroutine boundary
	use var
	implicit none

! Locat variables
integer::i,j,k, n

! Velocity boundary condition
do i = xmin,xmax
do j = ymin,ymax
k = zmax
do n = 0,18
f(i,j,k,n) = feq(i,j,k,n)
end do
end do
end do

!--------------------Bounce back----------------------
! On side planes, excluding edges and corner points
!-----------------------------------------------------
! Wall X = xmin
do j = ymin+1, ymax-1
do k = zmin+1, zmax-1
i = xmin
f(i,j,k,1 ) = f(i,j,k,2 )
f(i,j,k,7 ) = f(i,j,k,10)
f(i,j,k,8 ) = f(i,j,k,9 )
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,13) = f(i,j,k,12)
end do
end do

! Wall X = xmax
do j = ymin+1,ymax-1
do k = zmin+1,zmax-1
i = xmax
f(i,j,k,2 ) = f(i,j,k,1 )
f(i,j,k,9 ) = f(i,j,k,8 )
f(i,j,k,10) = f(i,j,k,7 )
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,14) = f(i,j,k,11)
end do
end do

! Wall Y = ymin
do i = xmin+1,xmax-1
do k = zmin+1,zmax-1
j = ymin
f(i,j,k,3 ) = f(i,j,k,4 )
f(i,j,k,7 ) = f(i,j,k,10)
f(i,j,k,9 ) = f(i,j,k,8 )
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,16) = f(i,j,k,17)
end do
end do

! Wall Y = ymax
do i = xmin+1,xmax-1
do k = zmin+1,zmax-1
j = ymax
f(i,j,k,4 ) = f(i,j,k,3 )
f(i,j,k,8 ) = f(i,j,k,9 )
f(i,j,k,10) = f(i,j,k,7 )
f(i,j,k,17) = f(i,j,k,16)
f(i,j,k,18) = f(i,j,k,15)
end do
end do

! Wall Z = zmin
do i = xmin+1,xmax-1
do j = ymin+1,ymax-1
k = zmin
f(i,j,k,5 ) = f(i,j,k,6 )
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,17) = f(i,j,k,16)
end do
end do

!--------------------------------------------
!-------------bounceback on edges------------
!--------------------------------------------
! Line xmin-zmin–[1]
do j = ymin+1,ymax-1
i = xmin
k = zmin
f(i,j,k,1 ) = f(i,j,k,2 )
f(i,j,k,5 ) = f(i,j,k,6 )
f(i,j,k,7 ) = f(i,j,k,10)
f(i,j,k,8 ) = f(i,j,k,9 )
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,12) = feq(i,j,k,12)
f(i,j,k,13) = feq(i,j,k,13)
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,17) = f(i,j,k,16)
end do

! Line ymin-zmin—[2]
do i = xmin+1,xmax-1
j = ymin
k = zmin
f(i,j,k,3 ) = f(i,j,k,4 )
f(i,j,k,5 ) = f(i,j,k,6 )
f(i,j,k,7 ) = f(i,j,k,10)
f(i,j,k,9 ) = f(i,j,k,8 )
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,16) = feq(i,j,k,16)
f(i,j,k,17) = feq(i,j,k,17)
end do

! Line xmax-zmin—[3]
do j = ymin+1,ymax-1
i = xmax
k = zmin
f(i,j,k,2 ) = f(i,j,k,1 )
f(i,j,k,5 ) = f(i,j,k,6 )
f(i,j,k,9 ) = f(i,j,k,8 )
f(i,j,k,10) = f(i,j,k,7 )
f(i,j,k,11) = feq(i,j,k,11)
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,14) = feq(i,j,k,14)
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,17) = f(i,j,k,16)
end do

! Line ymax-zmin—[4]
do i = xmin+1,xmax-1
j = ymax
k = zmin
f(i,j,k,4 ) = f(i,j,k,3 )
f(i,j,k,5 ) = f(i,j,k,6 )
f(i,j,k,8 ) = f(i,j,k,9 )
f(i,j,k,10) = f(i,j,k,7 )
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,15) = feq(i,j,k,15)
f(i,j,k,17) = f(i,j,k,16)
f(i,j,k,18) = feq(i,j,k,18)
end do

! Line xmin-ymin—[5]
do k = zmin+1,zmax-1
i = xmin
j = ymin
f(i,j,k,1 ) = f(i,j,k,2 )
f(i,j,k,3 ) = f(i,j,k,4 )
f(i,j,k,7 ) = f(i,j,k,10)
f(i,j,k,8 ) = feq(i,j,k,8)
f(i,j,k,9 ) = feq(i,j,k,9)
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,13) = f(i,j,k,12)
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,16) = f(i,j,k,17)
end do

! Line xmax-ymin—[6]
do k = zmin+1,zmax-1
i = xmax
j = ymin
f(i,j,k,2 ) = f(i,j,k,1 )
f(i,j,k,3 ) = f(i,j,k,4 )
f(i,j,k,7 ) = feq(i,j,k,7)
f(i,j,k,9 ) = f(i,j,k,8 )
f(i,j,k,10) = feq(i,j,k,10)
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,14) = f(i,j,k,11)
f(i,j,k,15) = f(i,j,k,18)
f(i,j,k,16) = f(i,j,k,17)
end do

! Line xmax-ymax—[7]
do k = zmin+1,zmax-1
i = xmax
j = ymax
f(i,j,k,2 ) = f(i,j,k,1 )
f(i,j,k,4 ) = f(i,j,k,3 )
f(i,j,k,8 ) = feq(i,j,k,8)
f(i,j,k,9 ) = feq(i,j,k,9)
f(i,j,k,10) = f(i,j,k,7 )
f(i,j,k,12) = f(i,j,k,13)
f(i,j,k,14) = f(i,j,k,11)
f(i,j,k,17) = f(i,j,k,16)
f(i,j,k,18) = f(i,j,k,15)
end do

! Line xmin-ymax—[8]
do k = zmin+1,zmax-1
i = xmin
j = ymax
f(i,j,k,1 ) = f(i,j,k,2 )
f(i,j,k,4 ) = f(i,j,k,3 )
f(i,j,k,7 ) = feq(i,j,k,7)
f(i,j,k,8 ) = f(i,j,k,9 )
f(i,j,k,10) = feq(i,j,k,10)
f(i,j,k,11) = f(i,j,k,14)
f(i,j,k,13) = f(i,j,k,12)
f(i,j,k,17) = f(i,j,k,16)
f(i,j,k,18) = f(i,j,k,15)
end do
!---------------------------------------------------
!----------------Corner points----------------------
!---------------------------------------------------
! xmin-ymin-zmin—[1]
f(xmin,ymin,xmin,1 ) = f(xmin,ymin,zmin,2 )
f(xmin,ymin,xmin,3 ) = f(xmin,ymin,zmin,4 )
f(xmin,ymin,xmin,5 ) = f(xmin,ymin,zmin,6 )
f(xmin,ymin,xmin,7 ) = f(xmin,ymin,zmin,10)
f(xmin,ymin,zmin,8 ) = feq(xmin,ymin,zmin,8 )
f(xmin,ymin,zmin,9 ) = feq(xmin,ymin,zmin,9 )
f(xmin,ymin,zmin,11) = f(xmin,ymin,zmin,14)
f(xmin,ymin,zmin,12) = feq(xmin,ymin,zmin,12)
f(xmin,ymin,zmin,13) = feq(xmin,ymin,zmin,13)
f(xmin,ymin,zmin,15) = f(xmin,ymin,zmin,18)
f(xmin,ymin,zmin,16) = feq(xmin,ymin,zmin,16)
f(xmin,ymin,zmin,17) = feq(xmin,ymin,zmin,17)

! xmax-ymin-zmin—[2]
f(xmax,ymin,zmin,2 ) = f(xmax,ymin,zmin,1 )
f(xmax,ymin,zmin,3 ) = f(xmax,ymin,zmin,4 )
f(xmax,ymin,zmin,5 ) = f(xmax,ymin,zmin,6 )
f(xmax,ymin,zmin,7 ) = feq(xmax,ymin,zmin,7 )
f(xmax,ymin,zmin,9 ) = f(xmax,ymin,zmin,8 )
f(xmax,ymin,zmin,10) = feq(xmax,ymin,zmin,10)
f(xmax,ymin,zmin,11) = feq(xmax,ymin,zmin,11)
f(xmax,ymin,zmin,12) = f(xmax,ymin,zmin,13)
f(xmax,ymin,zmin,14) = feq(xmax,ymin,zmin,14)
f(xmax,ymin,zmin,15) = f(xmax,ymin,zmin,18)
f(xmax,ymin,zmin,16) = feq(xmax,ymin,zmin,16)
f(xmax,ymin,zmin,17) = feq(xmax,ymin,zmin,17)

! xmax-ymax-zmin—[3]
f(xmax,ymax,zmin,2 ) = f(xmax,ymax,zmin,1 )
f(xmax,ymax,zmin,4 ) = f(xmax,ymax,zmin,3 )
f(xmax,ymax,zmin,5 ) = f(xmax,ymax,zmin,6 )
f(xmax,ymax,zmin,8 ) = feq(xmax,ymax,zmin,8 )
f(xmax,ymax,zmin,9 ) = feq(xmax,ymax,zmin,9 )
f(xmax,ymax,zmin,10) = f(xmax,ymax,zmin,7 )
f(xmax,ymax,zmin,11) = feq(xmax,ymax,zmin,11)
f(xmax,ymax,zmin,12) = f(xmax,ymax,zmin,13)
f(xmax,ymax,zmin,14) = feq(xmax,ymax,zmin,14)
f(xmax,ymax,zmin,15) = feq(xmax,ymax,zmin,15)
f(xmax,ymax,zmin,17) = f(xmax,ymax,zmin,16)
f(xmax,ymax,zmin,18) = feq(xmax,ymax,zmin,18)

! xmin-ymax-zmin—[4]
f(xmin,ymax,zmin,1 ) = f(xmin,ymax,zmin,2 )
f(xmin,ymax,zmin,4 ) = f(xmin,ymax,zmin,3 )
f(xmin,ymax,zmin,5 ) = f(xmin,ymax,zmin,6 )
f(xmin,ymax,zmin,7 ) = feq(xmin,ymax,zmin,7 )
f(xmin,ymax,zmin,8 ) = f(xmin,ymax,zmin,9 )
f(xmin,ymax,zmin,10) = feq(xmin,ymax,zmin,10)
f(xmin,ymax,zmin,11) = f(xmin,ymax,zmin,14)
f(xmin,ymax,zmin,12) = feq(xmin,ymax,zmin,12)
f(xmin,ymax,zmin,13) = feq(xmin,ymax,zmin,13)
f(xmin,ymax,zmin,15) = feq(xmin,ymax,zmin,15)
f(xmin,ymax,zmin,17) = f(xmin,ymax,zmin,16)
f(xmin,ymax,zmin,18) = feq(xmin,ymax,zmin,18)

	return
	end subroutine boundary

!-------------------------------------------------------------------------------------------------------------------------------

	subroutine hydrodynamics
	use var
	implicit none

! Local variables
integer::i,j,k
integer::i1,i2,j1,j2,k1,k2
real(kind = dp):: invRho, ux, uy, uz

	do i = xmin,xmax
	do j = ymin,ymax
	do k = zmin,zmax-1

	   ut(i,j,k,1) = u(i,j,k,1)
	   ut(i,j,k,2) = u(i,j,k,2)
	   ut(i,j,k,3) = u(i,j,k,3)

	   rho(i,j,k) = f(i,j,k, 0)               &
		      + f(i,j,k, 1) + f(i,j,k, 2) &
                          + f(i,j,k, 3) + f(i,j,k, 4) &
	              + f(i,j,k, 5) + f(i,j,k, 6) &
		      + f(i,j,k, 7) + f(i,j,k, 8) &
                          + f(i,j,k, 9) + f(i,j,k,10) &
	              + f(i,j,k,11) + f(i,j,k,12) &
		      + f(i,j,k,13) + f(i,j,k,14) &
		      + f(i,j,k,15) + f(i,k,j,16) &
		      + f(i,j,k,17) + f(i,j,k,18) 
	              		   

	   u(i,j,k,1) = f(i,j,k, 1)*cx(1)  + f(i,j,k, 2)*cx(2)  &
                          + f(i,j,k, 3)*cx(3)  + f(i,j,k, 4)*cx(4)  &
	              + f(i,j,k, 5)*cx(5)  + f(i,j,k, 6)*cx(6)  &
		      + f(i,j,k, 7)*cx(7)  + f(i,j,k, 8)*cx(8)  &
                          + f(i,j,k, 9)*cx(9)  + f(i,j,k,10)*cx(10) &
	              + f(i,j,k,11)*cx(11) + f(i,j,k,12)*cx(12) &
		      + f(i,j,k,13)*cx(13) + f(i,j,k,14)*cx(14) &
		      + f(i,j,k,15)*cx(15) + f(i,k,j,16)*cx(16) &
		      + f(i,j,k,17)*cx(17) + f(i,j,k,18)*cx(18) 

	   u(i,j,k,2) = f(i,j,k, 1)*cy(1)  + f(i,j,k, 2)*cy(2)  &
                          + f(i,j,k, 3)*cy(3)  + f(i,j,k, 4)*cy(4)  &
	              + f(i,j,k, 5)*cy(5)  + f(i,j,k, 6)*cy(6)  &
		      + f(i,j,k, 7)*cy(7)  + f(i,j,k, 8)*cy(8)  &
                          + f(i,j,k, 9)*cy(9)  + f(i,j,k,10)*cy(10) &
	              + f(i,j,k,11)*cy(11) + f(i,j,k,12)*cy(12) &
		      + f(i,j,k,13)*cy(13) + f(i,j,k,14)*cy(14) &
		      + f(i,j,k,15)*cy(15) + f(i,k,j,16)*cy(16) &
		      + f(i,j,k,17)*cy(17) + f(i,j,k,18)*cy(18) 

	   u(i,j,k,3) = f(i,j,k, 1)*cz(1)  + f(i,j,k, 2)*cz(2)  &
                          + f(i,j,k, 3)*cz(3)  + f(i,j,k, 4)*cz(4)  &
	              + f(i,j,k, 5)*cz(5)  + f(i,j,k, 6)*cz(6)  &
		      + f(i,j,k, 7)*cz(7)  + f(i,j,k, 8)*cz(8)  &
                          + f(i,j,k, 9)*cz(9)  + f(i,j,k,10)*cz(10) &
	              + f(i,j,k,11)*cz(11) + f(i,j,k,12)*cz(12) &
		      + f(i,j,k,13)*cz(13) + f(i,j,k,14)*cz(14) &
		      + f(i,j,k,15)*cz(15) + f(i,k,j,16)*cz(16) &
		      + f(i,j,k,17)*cz(17) + f(i,j,k,18)*cz(18) 

	   !invRho = 1.D0/rho(i,j,k)
	
	   u(i,j,k,1) = u(i,j,k,1)/rho(i,j,k)
	   u(i,j,k,2) = u(i,j,k,2)/rho(i,j,k)
	   u(i,j,k,3) = u(i,j,k,3)/rho(i,j,k)

	end do
	end do
	end do

!---------------------------------------------------------------
do i = xmin,xmax
do j = ymin,ymax
k1 = zmin
k2 = zmax

	   u(i,j,k1,1) = 0.D0
	   u(i,j,k1,2) = 0.D0
	   u(i,j,k1,3) = 0.D0

	   u(i,j,k2,1) = Uo
	   u(i,j,k2,2) = 0.D0
	   u(i,j,k2,3) = 0.D0

	   rho(i,j,k2) = 1.D0

	end do
	end do

	do j = ymin,ymax
	do k = zmin,zmax-1
	   i1 = xmin
	   i2 = xmax

	   u(i1,j,k,1) = 0.D0
	   u(i1,j,k,2) = 0.D0
	   u(i1,j,k,3) = 0.D0

	   u(i2,j,k,1) = 0.D0
	   u(i2,j,k,2) = 0.D0
	   u(i2,j,k,3) = 0.D0

	end do
	end do

	do i = xmin,xmax
	do k = zmin,zmax-1
	   j1 = ymin
	   j2 = ymax

	   u(i,j1,k,1) = 0.D0
	   u(i,j1,k,2) = 0.D0
	   u(i,j1,k,3) = 0.D0

	   u(i,j2,k,1) = 0.D0
	   u(i,j2,k,2) = 0.D0
	   u(i,j2,k,3) = 0.D0

	end do
	end do		

	end subroutine hydrodynamics

!------------------------------------------------------------------------------------------------------------------

	subroutine collision
	use var
	implicit none

! Local variables
integer:: i,j,k,n
real(kind = dp):: ux, uy, uz, rhon
real(kind = dp):: v, w, U_sq

!-----Calculating equillibrium distribution function-------------------------
do i = xmin,xmax
do j = ymin,ymax
do k = zmin,zmax

	   ux   = u(i,j,k,1)
	   uy   = u(i,j,k,2)
	   uz   = u(i,j,k,3)
	   U_sq = 0.5D0*( ux*ux + uy*uy + uz*uz )
	   rhon = rho(i,j,k)		
	 
	   v = ux*cx(0) + uy*cy(0) + uz*cz(0)
	   feq(i,j,k,0) = rhon*(W0 + W0C*( v + 1.5D0*v*v - U_sq ) )
	
	   do n = 1,6		
	      v = ux*cx(n) + uy*cy(n) + uz*cz(n)
	      feq(i,j,k,n) = rhon*(W1 + W1C*( v + 1.5D0*v*v - U_sq ) )
	   end do

	   do n = 7,18		
	      v = ux*cx(n) + uy*cy(n) + uz*cz(n)
	      feq(i,j,k,n) = rhon*(W2 + W2C*( v + 1.5D0*v*v - U_sq ) )
	   end do

	end do
	end do
	end do

!-------------------------Collision-------------------------------
do i = xmin,xmax
do j = ymin,ymax
do k = zmin,zmax-1
do n = 0,18
f(i,j,k,n) = f(i,j,k,n) - (f(i,j,k,n) - feq(i,j,k,n))/tau
end do
end do
end do
end do
! Velocity boundary condition
do i = xmin,xmax
do j = ymin,ymax
k = zmax
do n = 0,18
f(i,j,k,n) = feq(i,j,k,n)
end do
end do
end do

	end subroutine collision

!--------------------------------------------------------------------------------------------------------
subroutine stats
use var
implicit none

! Local variables
integer :: i,j,k
real(kind = dp) :: du,dv,dw
real(kind = dp) :: err1, err2
real(kind = dp) :: ux, uy, uz

       	err1 = 0.D0
	err2 = 0.D0

          	do i = xmin,xmax
    	do j = ymin,ymax
	do k = zmin,zmax

	   ux = u(i,j,k,1)
	   uy = u(i,j,k,2)
	   uz = u(i,j,k,3)

	   du = u(i,j,k,1) - ut(i,j,k,1)
	   dv = u(i,j,k,2) - ut(i,j,k,2)
	   dw = u(i,j,k,3) - ut(i,j,k,3)

               err1 = err1 + dsqrt( du*du + dv*dv + dw*dw )            
	   err2 = err2 + dsqrt( ux*ux + uy*uy + uz*uz )
                 
     	end do
     	end do
	end do

	errr = err1/err2

!
print*,iStep,errr

	open(unit = 10, file = "error.out", status = "UNKNOWN", position = "APPEND")
	   write(10,'(I9,6ES19.9)')istep,errr
         close(unit = 10)

            end subroutine stats

!----------------------------------------------------------------------------------------------------------------------------------------------

Tanu · June 25, 2014, 2:37am

====================
du = u(i,j,k,1) - ut(i,j,k,1)
dv = u(i,j,k,2) - ut(i,j,k,2)
dw = u(i,j,k,3) - ut(i,j,k,3)

where you define the variable ut? I am not sure may be there are no linked of du, dv and dw with ut(i,j,k,1),ut(i,j,k,2),ut(i,j,k,3) .Please check it once again.

Tanu

Ashok · June 25, 2014, 6:02am

Hi,

ut(i,j,k) are defined in var module. The are used in hydrodynamics.f90 subroutine to store the values of velocity of previous time step.

Tanu · June 25, 2014, 7:46am

Hi Ashok,
If you would like to send your FORTRAN code, then I will try to check it.

With all the best
Tanu

Ashok · June 26, 2014, 7:37am

HI,

Give me your mail ID, I will mail you the code.
Thanks.

arash81069 · June 29, 2014, 7:36pm

where do u find this code ?
the code that u enter here have not syntax error except :

	    if(errr.lt.eps ) then
	       call output
	       stop
	    end if

there is not subroutine witth the name output!!!
if u want give me your fortran code i thinck your code may have some convergence mistake
resad1391@gmail.com
best wishes

3D Cavity implementation issues

==================== du = u(i,j,k,1) - ut(i,j,k,1) dv = u(i,j,k,2) - ut(i,j,k,2) dw = u(i,j,k,3) - ut(i,j,k,3)

====================
du = u(i,j,k,1) - ut(i,j,k,1)
dv = u(i,j,k,2) - ut(i,j,k,2)
dw = u(i,j,k,3) - ut(i,j,k,3)