waLBerla Tutorial 02: LBM Setup

In this tutorial we set up a basic lattice Boltzmann simulation of a periodic channel.

First we create a simple block structure, same as in the last "Game of Life" tutorial.



In [ ]:

    
from waLBerla import *
from material.matplotlib_setup import *

blocks = createUniformBlockGrid( cells=(500,200,1), periodic=(1,0,1) )

For LBM we have the choice of different lattice models. A lattice model in waLBerla consists of:

a stencil that defines the cell neighborhoods. In this tutorial we use a two dimensional D2Q9 stencil.
a force model to use volume forces in the domain. For our setup we need a constant volume force to drive the channel and make use of the "SimpleConstant" model.
a collision model defining the LBM collision operator. Her we use SRT (BGK), however the two relaxation time (TRT) model is also supported and for the D3Q19 model a MRT implementation is available as well

With the lattice model a pdf field is created and added to all blocks. Additionally we create optional adaptors for velocity and density. These adaptors behave similar to fields, however no memory is required for them, the values are computed on-the-fly from the pdfs.



In [ ]:

    
omega = 1.3
collisionModel =lbm.collisionModels.SRT(omega)
forceModel = lbm.forceModels.SimpleConstant( (1e-5,0,0) )
latticeModel = lbm.makeLatticeModel("D2Q9", collisionModel, forceModel)
lbm.addPdfFieldToStorage(blocks, "pdfs", latticeModel, 
                         velocityAdaptor="vel", densityAdaptor="rho", initialDensity=1.0)

To store geometry information a flag field is created and used as input to the LBM boundary handling:



In [ ]:

    
field.addFlagFieldToStorage( blocks, 'flags')
lbm.addBoundaryHandlingToStorage(blocks, 'boundary', 'pdfs', 'flags')

Next we define the setup of our channel. At the northern and southern domain borders no-slip boundaries are required. To set up the boundaries the distributed nature of the domain has to be taken into account. When iterating over the blocks we only get the locally stored blocks. So this loop is automatically parallelized when using multiple MPI processes. For each block we first check if it is located at the border of the global domain, and if so, we set the border slice to no-slip. When defining the slice, the additional 'g' marker is used, to place the boundary in the ghost layer.



In [ ]:

    
def setBoundariesChannel( blocks, boundaryHandlingID ):
    for block in blocks:
        b = block[ boundaryHandlingID ]
        if block.atDomainMinBorder[1]:
            b.forceBoundary('NoSlip',   makeSlice[ :, 0, :, 'g'])
        if block.atDomainMaxBorder[1]:
            b.forceBoundary('NoSlip',   makeSlice[ :,-1, :, 'g'])
        b.fillWithDomain()
setBoundariesChannel( blocks, 'boundary' )

Finally a LBM sweep functor is created, making use of the already created pdf and flag field.



In [ ]:

    
sweep = lbm.makeCellwiseSweep(blocks, "pdfs", flagFieldID='flags', flagList=['fluid'])

Since we have set our domain to be periodic in x-direction,the east and west boundary are taken care of by the ghost layer exchange. Here the communication is set up as explained in the previous tutorial.



In [ ]:

    
communication = createUniformBufferedScheme( blocks, 'D2Q9')
communication.addDataToCommunicate( field.createPackInfo( blocks, 'pdfs') )

One timestep consists of ghost layer communication, boundary handling and an LBM stream-collide step



In [ ]:

    
def run( timesteps ):
    for t in range(timesteps):
        communication()
        for block in blocks: block['boundary']()
        for block in blocks: sweep.streamCollide( block )

To get a more interesting result we place an airfoil shaped object inside our channel. This geometry is loaded from an image placed in the 'material' subfolder. You can upload different images and simulate the flow around these objects:



In [ ]:

    
from waLBerla.geometry_setup import *
setBoundaryFromBlackAndWhiteImage(blocks, "boundary", makeSlice[0.25:0.75, 0.3:0.6 ,0.5], 
                                  "material/wing.png", "NoSlip")

Next we create a video of the density distribution in the channel ( which is proportional to pressure in LBM). Before creating a frame for the video we let the simulation run for 10 timesteps, then we set the density inside the object to NaN to get a better looking plot



In [ ]:

    
def setByMask( blocks, targetField, targetValue, flagField, flagNames ):
    for b in blocks:
        mask = 0
        for flagName in flagNames:
            mask |= b[flagField].flag(flagName)

        targetArr = field.toArray( b[targetField], True )
        flagArr   = field.toArray( b[flagField]  , True )[:,:,:,0]
        targetArr[ np.bitwise_and( flagArr, mask ) > 0, : ] = targetValue


    
def visualizationStep():
    run(10)
    # Set pdfs to NaN in NoSlip cells to better see the obstacle in the visualization:
    setFieldUsingFlagMask(blocks, 'pdfs', np.NaN, 'flags', ['NoSlip'] )



In [ ]:

    
import waLBerla.plot as wplt

run(700) # run the simulation until flow field has stabilized - then no colormap rescaling is required
ani = wplt.scalarFieldAnimation( blocks, 'rho', makeSlice[:,:,0.5], visualizationStep, frames=300)
displayAsHtmlImage(ani)