Hi Dear,
I tried to write my first ever LBM code to simulate 2D-cavity but it is producing erroneous results infact meaning less, wrong results. I have pasted the whole code and need comments and suggestions how to make it run with good/appropriate results. An urgent help is required. The code:
[code=“fortran”]
program lbe2d
use globalparameters
use intialization
use boundaryconditions
use macroscopicvariables
implicit none
integer :: time, i
character (len=30) :: format
integer, parameter :: iterations = 1000
real*8 :: uu
!> allocate dimensions for dynamic arrays
allocate(f(nx,ny,npop), feq(nx,ny,npop), ftemp(nx,ny,npop), rho(nx,ny), u(nx,ny), v(nx,ny))
!
nu = (2.0d0*tau-1)/6.0d0
uinf = 0.1d0
do i = 1, nx
u(i,ny) = uinf
enddo
print*, u(:,ny)
! Initialization
call initialize
! main loop for time marching
time = 0
timeStep: do
time = time + 1
call streaming
call collision
call bounceback
call macrovar
call compute_feq
write(6,'(i6,a2,3(e16.8,a2))') time,' ', u(nx/2, ny/2)
if (time == iterations) exit
enddo timeStep
deallocate(f, feq, ftemp, rho, u, v)
end program lbe2d
[code="fortran"]
module globalparameters
implicit none
real*8, parameter :: pi = 4.0d0 * atan(1.0d0)
real*8, parameter :: quarter_pi = atan(1.0d0)
real*8, parameter :: cs = 1.0d0/dsqrt(3.0d0)
real*8, parameter :: rho0 = 1.0d0, u0 = 0.0d0, v0 = 0.0d0, tau = 0.6d0
integer, parameter :: nx = 128, ny = 128, npop = 9
real*8, parameter, dimension(npop) :: weights = (/4.d0/9.d0, 1.d0/9.d0, 1.d0/9.d0, 1.d0/9.d0, 1.d0/9.d0, &
1.d0/36.d0, 1.d0/36.d0, 1.d0/36.d0, 1.d0/36.d0/)
real*8, parameter, dimension(npop) :: ex = (/0, 1, 1, 0, -1, -1, -1, 0, 1/)
real*8, parameter, dimension(npop) :: ey = (/0, 0, 1, 1, 1, 0, -1, -1, -1/)
real*8 :: nu, uinf
real*8, allocatable, dimension(:,:) :: u, v, rho, p
real*8, allocatable, dimension(:,:,:) :: f, feq, ftemp
end module globalParameters
[code=“fortran”]
module intialization
!> Initialization subroutine:
use globalparameters
implicit none
contains
subroutine initialize
implicit none
integer:: x, y
character (len=30) :: format
rho = rho0
u = u0
v = v0
!> Calling “compute_feq” subroutine to compute the equilibrium distribution function feq.
call compute_feq
!> Set the initial distribution function to feq i.e. f = feq
f = feq
return
end subroutine initialize
!> This subroutine computes the local equilibrium distribution function.
!> Equilibrium distribution function at all node points based on the
!> current velocity and density fields.
subroutine compute_feq
implicit none
integer :: ie
do ie = 1, npop
feq(:,:,ie) = weights(ie)*rho(:,:) * (1.0d0 + 3.0d0*(ex(ie)*u(:,:) + ey(ie)*v(:,:)) &
+ (9.0d0/2.0d0)*(ex(ie)*u(:,:) + ey(ie)*v(:,:))*(ex(ie)*u(:,:) + ey(ie)*v(:,:)) &
- (3.0d0/2.0d0)*(u(:,:)*u(:,:) + v(:,:)*v(:,:)))
enddo
return
end subroutine compute_feq
!> This subroutine updates the f field at all node points i.e. Streaming
subroutine streaming
implicit none
integer:: ie, x, y, x1, y1
do ie = 1, npop
do y = 1, ny
do x = 1, nx
x1 = 1 + modulo(x - 1 + int(ex(ie)) + nx, nx)
y1 = 1 + modulo(y - 1 + int(ey(ie)) + ny, ny)
ftemp(x1,y1,ie) = f(x,y,ie)
enddo
enddo
enddo
return
end subroutine streaming
!> Collision
subroutine collision
implicit none
integer :: ie
do ie = 1, npop
ftemp(:,:,ie) = f(:,:,ie) - (f(:,:,ie) - feq(:,:,ie))/tau
enddo
return
end subroutine collision
!> End of module initialization
end module intialization
[code="fortran"]
module boundaryconditions
use globalparameters
contains
!> 5 4 3 ^y
!> \|/ | x
!> 6-1-2 --->
!> /|\
!> 7 8 9
subroutine bounceback
! this is for no-slip boundary with Bounce back scheme
implicit none
integer :: ie
! for lower boundary
do ie = 3, 5
ftemp(:,1,ie) = ftemp(:,1,ie+4)
enddo
! for upper boundary
do ie = 7, 9
ftemp(:,ny,ie) = ftemp(:,ny,ie-4)
enddo
! right side boundary
do ie = 5, 7
if (ie == 5) then
ftemp(1,:,ie) = ftemp(1,:,ie+4)
else
ftemp(1,:,ie) = ftemp(1,:,ie-4)
endif
enddo
! left side boundary
do ie = 2, 3
ftemp(nx,:,ie) = ftemp(nx,:,ie+4)
enddo
ftemp(nx,:,9) = ftemp(nx,:,5)
end subroutine bounceback
end module boundaryconditions
[code=“fortran”]
!> This module computes physical/macroscopic variables, i.e. rho, u, and v
module macroscopicvariables
use globalParameters
use intialization
use boundaryconditions
contains
subroutine macrovar
implicit none
integer:: ie
!> Set the distribution function f
f = ftemp
rho = 0.0d0
u = 0.0d0
v = 0.0d0
do ie = 1, npop
rho(:,:) = f(:,:,ie)
u(:,:) = u(:,:) + ex(ie)*f(:,:,ie)
v(:,:) = v(:,:) + ey(ie)*f(:,:,ie)
enddo
p = rho*cs*cs
u = u/rho0
v = v/rho0
return
end subroutine macrovar
end module macroscopicVariables