This code consists of three files:
Here is gridsize.inc:
PARAMETER( nx = 128, nz = 128 )
Next, we have the global variables file - common.inc:
C Variables and common blocks used in lbm.f
DOUBLE PRECISION x(0:nx+1),z(0:nz+1)
DOUBLE PRECISION u(0:nx+1,0:nz+1),w(0:nx+1,0:nz+1)
DOUBLE PRECISION f(0:nx+1,0:nz+1,0:8)
DOUBLE PRECISION fprev(0:nx+1,0:nz+1,0:8)
DOUBLE PRECISION feq(0:nx+1,0:nz+1,0:8),omega(0:nx+1,0:nz+1,0:8)
DOUBLE PRECISION feq_top(0:nx+1,0:8),feq_init(0:nx+1,0:nz+1,0:8)
DOUBLE PRECISION ex(0:8),ez(0:8),alpha(0:8)
DOUBLE PRECISION den(0:nx+1,0:nz+1),p(0:nx+1,0:nz+1)
COMMON /BLOCK_01/ xl,zl,x,z,dx,dz
COMMON /BLOCK_03/ zero,one,rinf,TINY,PI
COMMON /BLOCK_04/ epsnss
COMMON /BLOCK_05/ den,p
COMMON /BLOCK_06/ u,w
COMMON /BLOCK_07/ f,fprev,feq,feq_top,feq_init,omega
COMMON /BLOCK_08/ ex,ez,alpha
COMMON /BLOCK_14/ dt
COMMON /BLOCK_20/ Re,UINF,densf,tau
And finally, the LBM implementation file - lbm2d.f:
C23456789012345678901234567890123456789012345678901234567890123456789012
C ------------------------------------------------------------------
C LATTICE BOLTZMANN METHOD FOR SIMULATING DRIVEN CAVITY FLOW
C ------------------------------------------------------------------
C
C * 2-Dimensional geometry - Uniform grid along both X and Z
C * Bhatnagar-Gross-Krook (BGK) FINITE DIFFERENCE discretization
C * 9-velocity model { 0 1 2 3 4 5 6 7 8 }
C * Unsteady, incompressible (almost !), viscous flow
C * Bounce-back boundary conditions on stationary walls
C * Postprocessing using MATLAB
C
C BENCHMARK PROBLEMS
C
C 1. The 2-Dimensional driven cavity
C
C Written By : Abhijit S. Joshi
C Last modified : Thu 25 Oct 2001 at 12:30 hrs
C Remarks : Studying the effect of Mach number on accuracy
C References : [1] Simulation of Cavity Flow by the Lattice
C Boltzmann method - Hou et.al.(1995)
C Journal of Computational Physics,118,329-347
C [2] Accuracy and efficiency study of LBM for
C steady state flow simulations - Lai et.al.(2001)
C Numerical Heat Transfer(B),39,21-43
C ------------------------------------------------------------------
C Main Program
C ------------------------------------------------------------------
INCLUDE 'gridsize.inc' ! contains { nx , nz }
INCLUDE 'common.inc' ! contains declarations & common
Re = 100.0D0 ! Reynolds number
CALL grid ! set up grid for the problem
CALL initial ! specify initial conditions
CALL time_integration ! Evolve the "f" distributions
CALL write_results ! create output files
END
C ------------------------------------------------------------------
C * SETTING UP THE GRID *
C
C nx = no. of nodes along the X direction (1,2,3,....nx)
C nz = no. of nodes along the Z direction (1,2,3,....nz)
C ( x(i), z(k) ) = X-Z coordinates of node (i,k)
C
C ------------------------------------------------------------------
SUBROUTINE grid
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
WRITE(6,*)
WRITE(6,*)'Grid generation........'
WRITE(6,*)
c Calculating grid spacing along X and Z.
xl = 1.0D0 ! domain length along X in lattice units
zl = 1.0D0 ! domain length along Z in lattice units
dx = xl/FLOAT(nx-1) ! lattice spacing along X
dz = zl/FLOAT(nz-1) ! lattice spacing along Z
dt = dx ! time step
c Calculating the X and Z coordinates of the node points
x(1) = 0.0D0
x(nx) = xl
DO i=2,nx-1
x(i) = x(i-1) + dx
ENDDO
z(1) = 0.0D0
z(nz) = zl
DO k=2,nz-1
z(k) = z(k-1) + dz
ENDDO
c calculating the "direction vectors" for LBM. The 9-velocity model
c is used here with the numbers running from {0,1,2,3,4,5,6,7,8}
ex(0) = 0.0D0
ez(0) = 0.0D0
ex(1) = 1.0D0
ez(1) = 0.0D0
ex(2) = 1.0D0
ez(2) = 1.0D0
ex(3) = 0.0D0
ez(3) = 1.0D0
ex(4) = -1.0D0
ez(4) = 1.0D0
ex(5) = -1.0D0
ez(5) = 0.0D0
ex(6) = -1.0D0
ez(6) = -1.0D0
ex(7) = 0.0D0
ez(7) = -1.0D0
ex(8) = 1.0D0
ez(8) = -1.0D0
END
C ------------------------------------------------------------------
C This is the initialization for the 2-dimensional Lid driven cavity
C flow. The velocity is set to zero everywhere inside the cavity and
C density is set equal to unity. Later, velocity at the lid is set
C to UINF, which is designed to match the desired Reynolds number.
C As far as the LBM is concerned, the equilibrium distribution
C functions are calculated here based on the initial known velocity
C field both on the boundaries and interior nodes.
C ------------------------------------------------------------------
SUBROUTINE initial
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
WRITE(6,*)
WRITE(6,*)'Initialization........'
WRITE(6,*)
c Initially the fluid inside the cavity is at rest, and at t = 0
c the top layer starts moving at speed "U" to the right. The values
c of the properties are adjusted to obtain the desired Reynolds no.
densf = 1.0D0 ! reference density
tau = 0.60D0 ! relaxation time
viskin = (2.0D0*tau - 1.0D0)*dx**2
1 / (6.0D0*dt) ! kinematic viscosity
UINF = viskin * Re / xl ! suitable lid velocity
cs = dx/(dt*SQRT(3.0D0)) ! speed of sound
WRITE(6,*)'Parameter values used in the simulation:'
WRITE(6,*)
WRITE(6,*)'Reynolds number (Re) = ',Re
WRITE(6,*)'Domain size (xl zl) = ',xl,zl
WRITE(6,*)'Grid size (nx nz) = ',nx,nz
WRITE(6,*)'grid spacing dx = dz = ',dx
WRITE(6,*)'Time step (dt) = ',dt
WRITE(6,*)'Reference density = ',densf
WRITE(6,*)'Kinematic viscosity = ',viskin
WRITE(6,*)'Relaxation time = ',tau
WRITE(6,*)'Lid velocity (UINF) = ',UINF
WRITE(6,*)'Speed of sound (cs) = ',cs
WRITE(6,*)'Mach number = ',UINF/cs
WRITE(6,*)
DO k=1,nz
DO i=1,nx
den(i,k) = densf ! set fluid density
u(i,k) = 0.0D0 ! set fluid velocity = 0
w(i,k) = 0.0D0
ENDDO
ENDDO
DO i=1,nx ! top lid moves to the right
u(i,nz) = UINF
ENDDO
c weight functions used in the equilibrium distribution
alpha(0) = 4.0D0/9.0D0
DO id = 1,7,2
alpha(id) = 1.0D0/9.0D0
ENDDO
DO id = 2,8,2
alpha(id) = 1.0D0/36.0D0
ENDDO
c The equilibrium particle distributions are computed based on
c the initial velocity distributions.
CALL equilibrium_distribution
c initial equilibrium distributions are stored here
DO id=0,8
DO i=1,nx
DO k=1,nz
feq_init(i,k,id) = feq(i,k,id)
ENDDO
ENDDO
ENDDO
c particle distributions are assigned equilibrium values
DO id=0,8
DO k=1,nz
DO i=1,nx
f(i,k,id) = feq(i,k,id)
fprev(i,k,id) = f(i,k,id)
ENDDO
ENDDO
ENDDO
c Does the equilibrium distribution give back the initial density
c and velocity everywhere in the domain ?
WRITE(6,*)'Checking density & velocity after initialization'
CALL density ! recompute the density den(i,k)
CALL velocity ! recompute the velocity u(i,k), w(i,k)
d_den = 0.0D0
d_vex = 0.0D0
d_vez = 0.0D0
DO i=1,nx
d_den = MAX( DABS(den(i,nz) - densf), d_den)
d_vex = MAX( DABS(u(i,nz) - UINF), d_vex)
d_vez = MAX( DABS(w(i,nz) - 0.0D0), d_vez)
ENDDO
DO i=1,nx
DO k=1,nz-1
d_den = MAX( DABS(den(i,k) - densf), d_den)
d_vex = MAX( DABS(u(i,k) - 0.0D0), d_vex)
d_vez = MAX( DABS(w(i,k) - 0.0D0), d_vez)
ENDDO
ENDDO
WRITE(6,*)
WRITE(6,*)'Max.change in density = ',d_den
WRITE(6,*)' x-velocity = ',d_vex
WRITE(6,*)' z-velocity = ',d_vez
WRITE(6,*)
END
C ------------------------------------------------------------------
C Calculate the Equilibrium particle distribution functions at all
C node points based on the current velocity & density fields.
C ------------------------------------------------------------------
SUBROUTINE equilibrium_distribution
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
DO id = 0,8
DO k=1,nz
DO i=1,nx
feq(i,k,id) = den(i,k) * alpha(id) * ( 1.0
1 + 3.0 * (ex(id)*u(i,k) + ez(id)*w(i,k))
2 + 4.5 * (ex(id)*u(i,k) + ez(id)*w(i,k))**2
3 - 1.5 * (u(i,k)**2 + w(i,k)**2) )
ENDDO
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C TIME INTEGRATION ROUTINE
C
C The "f" distributions {0,1,2,3,4,5,6,7,8} are advanced in time
C based on the discrete LGBK model. (Explicit method; No iteration)
C ------------------------------------------------------------------
SUBROUTINE time_integration
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
time = 0.0D0
WRITE(6,*)
WRITE(6,*)'TIME INTEGRATION BEGINS'
WRITE(6,*)
c epsnss = 1.0D-6 ! criterion for steady state
DO itime=1,100000 ! time-integration begins
time = time + dt ! advance time to the next level
CALL collision ! evaluate the R.H.S.
CALL streaming ! update the node points
CALL boundary_conditions ! bounce back scheme implemented
CALL density ! construct the density field
CALL velocity ! construct the velocity field
CALL check_steady_state(dfdt)
IF(MOD(itime,100).EQ.0) WRITE(6,100)INT(itime),time,dfdt
IF(dfdt.LT.epsnss) GOTO 200
ENDDO ! time-integration ends
WRITE(6,*)'WARNING:- Terminating before steady state'
RETURN
200 WRITE(6,*)'Steady state reached. Terminating computation'
WRITE(6,*)
100 FORMAT(1X,'time step = ',I8,3X,'t = ',E12.5,3X,'dfdt = ',E12.5)
END
C ------------------------------------------------------------------
C Calculate the R.H.S. of the discrete LBGK equation
C ------------------------------------------------------------------
SUBROUTINE collision
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
CALL equilibrium_distribution
DO id = 0,8
DO k=1,nz
DO i=1,nx
omega(i,k,id) = - ( fprev(i,k,id) - feq(i,k,id) ) / tau
ENDDO
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C Update the "f" fields at all node points
C ------------------------------------------------------------------
SUBROUTINE streaming
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
DO id=0,8
DO k=1,nz
DO i=1,nx
idelx = INT( ex(id) )
idelz = INT( ez(id) )
f(i+idelx,k+idelz,id) = fprev(i,k,id) + omega(i,k,id)
ENDDO
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C Bounce Back Scheme
C ------------------------------------------------------------------
SUBROUTINE boundary_conditions
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
c top distribution functions set to their initial equilibrium values
DO id=0,8
DO i=1,nx
f(i,nz,id) = feq_init(i,nz,id)
ENDDO
ENDDO
c Bounce back scheme is used to determine the distribution functions
c on walls with the "no-slip" velocity boundary conditions
c Left wall ( f8 f1 f2 )
DO k=2,nz-1
f(1,k,8) = f(1,k,4)
f(1,k,1) = f(1,k,5)
f(1,k,2) = f(1,k,6)
ENDDO
c Right wall ( f4 f5 f6 )
DO k=2,nz-1
f(nx,k,4) = f(nx,k,8)
f(nx,k,5) = f(nx,k,1)
f(nx,k,6) = f(nx,k,2)
ENDDO
c Bottom wall ( f2 f3 f4 )
DO i=2,nx-1
f(i,1,2) = f(i,2,6)
f(i,1,3) = f(i,2,7)
f(i,1,4) = f(i,2,8)
ENDDO
c Bottom left corner ( f8 f1 f2 f3 f4 )
f(1,1,1) = f(1,1,5)
f(1,1,2) = f(1,1,6)
f(1,1,3) = f(1,1,7)
f(1,1,4) = ( den(1,1) - (f(1,1,0) + f(1,1,1) + f(1,1,2) + f(1,1,3)
1 + f(1,1,5) + f(1,1,6) + f(1,1,7)) ) * 0.5
f(1,1,8) = f(1,1,4)
c Bottom right corner ( f2 f3 f4 f5 f6 )
f(nx,1,3) = f(nx,1,7)
f(nx,1,4) = f(nx,1,8)
f(nx,1,5) = f(nx,1,1)
f(nx,1,2) = ( den(nx,1) - (f(nx,1,0) + f(nx,1,1) + f(nx,1,3)
1 + f(nx,1,4) + f(nx,1,5) + f(nx,1,7)
2 + f(nx,1,8)) ) * 0.5
f(nx,1,6) = f(nx,1,2)
END
C ------------------------------------------------------------------
C Compute the density at all node points
C ------------------------------------------------------------------
SUBROUTINE density
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
DO k=1,nz
DO i=1,nx
sum = 0.0D0
DO id = 0,8
sum = sum + f(i,k,id)
ENDDO
den(i,k) = sum
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C Compute the velocity at all node points
C ------------------------------------------------------------------
SUBROUTINE velocity
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
DO k=1,nz
DO i=1,nx
sumx = 0.0D0
sumz = 0.0D0
DO id = 0,8
sumx = sumx + f(i,k,id)*ex(id)
sumz = sumz + f(i,k,id)*ez(id)
ENDDO
u(i,k) = sumx / den(i,k)
w(i,k) = sumz / den(i,k)
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C Compute the pressure at all node points
C ------------------------------------------------------------------
SUBROUTINE pressure
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
c In the LBM, the density is not really a constant even for
c incompressible flows, and the pressure is proportional to the
c density. But in this problem, the pressure is not required in the
c calculation.
C_s = 1.0D0/DSQRT(3.0D0) ! speed of sound in the fluid
DO i=1,nx
DO k=1,nz
p(i,k) = C_s**2 * den(i,k)
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C Checking if the system has reached steady state
C ------------------------------------------------------------------
SUBROUTINE check_steady_state(dfdt)
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
dfdt = 0.0D0
DO k=1,nz
DO i=1,nx
DO id = 0,8
chf = DABS( f(i,k,id) - fprev(i,k,id) )
dfdt = MAX( dfdt, chf )
fprev(i,k,id) = f(i,k,id)
ENDDO
ENDDO
ENDDO
END
C ------------------------------------------------------------------
C Post-Processing
C ------------------------------------------------------------------
SUBROUTINE write_results
INCLUDE 'gridsize.inc'
INCLUDE 'common.inc'
c Velocity field information for MATLAB
OPEN(11,FILE = 'DAT/X.dat')
OPEN(12,FILE = 'DAT/Z.dat')
OPEN(13,FILE = 'DAT/U.dat')
OPEN(14,FILE = 'DAT/W.dat')
REWIND(11)
REWIND(12)
REWIND(13)
REWIND(14)
DO i=1,nx
DO k=1,nz
WRITE(11,*)x(i)
ENDDO
ENDDO
DO i=1,nx
DO k=1,nz
WRITE(12,*)z(k)
ENDDO
ENDDO
DO i=1,nx
DO k=1,nz
WRITE(13,*)u(i,k)
ENDDO
ENDDO
DO i=1,nx
DO k=1,nz
WRITE(14,*)w(i,k)
ENDDO
ENDDO
c Centerline velocity profiles for the driven cavity in
c dimensionless form
OPEN(15,FILE = 'DAT/uz.dat')
REWIND(15)
DO k=1,nz
WRITE(15,*)u((nx/2)+1,k)/UINF,z(k)/zl
ENDDO
OPEN(16,FILE = 'DAT/wx.dat')
REWIND(16)
DO i=1,nx
WRITE(16,*)x(i)/xl,w(i,(nz/2)+1)/UINF
ENDDO
c Scalar fields (contour plotting / 3-D surface plotting)
OPEN(21,FILE='DAT/x.m')
OPEN(22,FILE='DAT/z.m')
OPEN(23,FILE='DAT/u.m')
OPEN(24,FILE='DAT/w.m')
OPEN(25,FILE='DAT/den.m')
OPEN(40,FILE='DAT/f0.m')
OPEN(41,FILE='DAT/f1.m')
OPEN(42,FILE='DAT/f2.m')
OPEN(43,FILE='DAT/f3.m')
OPEN(44,FILE='DAT/f4.m')
OPEN(45,FILE='DAT/f5.m')
OPEN(46,FILE='DAT/f6.m')
OPEN(47,FILE='DAT/f7.m')
OPEN(48,FILE='DAT/f8.m')
REWIND(21)
REWIND(22)
REWIND(23)
REWIND(24)
REWIND(25)
REWIND(40)
REWIND(41)
REWIND(42)
REWIND(43)
REWIND(44)
REWIND(45)
REWIND(46)
REWIND(47)
REWIND(48)
C Write data for postprocessing using MATLAB (contour plotting)
WRITE(21,*)'x = [...'
WRITE(22,*)'z = [...'
WRITE(23,*)'u = [...'
WRITE(24,*)'w = [...'
WRITE(25,*)'den = [...'
WRITE(40,*)'f0 = [...'
WRITE(41,*)'f1 = [...'
WRITE(42,*)'f2 = [...'
WRITE(43,*)'f3 = [...'
WRITE(44,*)'f4 = [...'
WRITE(45,*)'f5 = [...'
WRITE(46,*)'f6 = [...'
WRITE(47,*)'f7 = [...'
WRITE(48,*)'f8 = [...'
DO k=1,nz
WRITE(21,100)(x(i),i=1,nx)
WRITE(22,100)(z(k),i=1,nx)
WRITE(23,100)(u(i,k),i=1,nx)
WRITE(24,100)(w(i,k),i=1,nx)
WRITE(25,100)(den(i,k),i=1,nx)
WRITE(40,100)(f(i,k,0),i=1,nx)
WRITE(41,100)(f(i,k,1),i=1,nx)
WRITE(42,100)(f(i,k,2),i=1,nx)
WRITE(43,100)(f(i,k,3),i=1,nx)
WRITE(44,100)(f(i,k,4),i=1,nx)
WRITE(45,100)(f(i,k,5),i=1,nx)
WRITE(46,100)(f(i,k,6),i=1,nx)
WRITE(47,100)(f(i,k,7),i=1,nx)
WRITE(48,100)(f(i,k,8),i=1,nx)
ENDDO
WRITE(21,*)'] ;'
WRITE(22,*)'] ;'
WRITE(23,*)'] ;'
WRITE(24,*)'] ;'
WRITE(25,*)'] ;'
WRITE(40,*)'] ;'
WRITE(41,*)'] ;'
WRITE(42,*)'] ;'
WRITE(43,*)'] ;'
WRITE(44,*)'] ;'
WRITE(45,*)'] ;'
WRITE(46,*)'] ;'
WRITE(47,*)'] ;'
WRITE(48,*)'] ;'
100 FORMAT(1X,256(E12.6,2X))
END
C ------------------------------------------------------------------
C23456789012345678901234567890123456789012345678901234567890123456789012
Note that this code is provided for historical reasons only. Perhaps it may be useful to understand some of the implementation details about the LBM that are hard to grasp just by reading papers.
If I were to write a similar code today in FORTRAN, there are several things I would do differently:
Global variables (the common blocks) make it very hard to debug and maintain code. A better approach is to use local variables inside each subroutine and function files and pass required variables as parameters.
It is always a good idea to declare each and every variable explicitly. No exceptions.
In the olden days of FORTRAN-77, we used implicit variables, where anything beginning with {I, J, K, L, M, N} were automatically integers. This is an absolute no-no now.
Always use:
PROGRAM lbm
IMPLICIT NONE
.
.
.
END
I was a big fan of using Matlab for post-processing results from my LBM codes.
The biggest problem with this approach is that you need Matlab, which is not free.
A much better option is to use open-source visualization tools like ParaView and Visit.
I now prefer that code write output data using XDMF + HDF5. This can be easily scaled across thousands of cores and can be read in parallel by ParaView and Visit.
You will note that there are no comments or instructions about how to build and run the code.
In this case, the build process was simple:
gfortran lbm2d.f -o lbm.x
But, in general, it is always good practice to use a Makefile.
It seems really shocking to me now, but I did not use any kind of version control tools at that time.
I would never dream of writing anything these days without using Git.