Tutorial - LBM 4: Complex Geometry

A configurable application for LBM simulations with complex geometries DISCLAIMER: for this tutorial, you must have built waLBerla with WALBERLA_BUILD_WITH_OPENMESH enabled!


In this tutorial, the basic LBM setup from Tutorial LBM 1 is extended by the capability of handling complex geometries loaded from mesh files. It has following features:

  • make use of default LBM sweeps from the lbm module
  • can handle some basic boundary conditions: no slip, free slip, pressure, and velocity boundaries
  • initialize geometry of the domain importing a mesh file
  • the boundary and geometry can be configured in a parameter file

The LBM routine mostly follows what has been demonstrated in LBM Tutorial 1. However, as our aim to load arbitrary complex geometry from a mesh file, we cannot specify the geometry via the Boundaries Block in the configuration file. We use waLBerla's mesh module that utilizes OpenMesh for loading, storing and manipulating meshes. Note that you must build waLBerla with OpenMesh capabilities in order to do this tutorial.

Stanford Bunny (loaded from an object file) in a 3D channel

General Setup

In this application, we will again simulate a channel flow around an obstacle. As before, the application will be fully configurable by a configuration file.

Lattice Boltzmann Data Structures

Lattice Model

In this tutorial, we will use the standard D3Q27 stencil with a SRT collision model.

using LatticeModel_T = lbm::D3Q27<lbm::collision_model::SRT>;
using Stencil_T = LatticeModel_T::Stencil;
using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;

As this is only the basic LBM setup, other stencils and collision methods can be included to simulate more accurate, stable and complex flows.

Mesh and Domain

When building waLBerla with the CMAKE configuration WALBERLA_BUILD_WITH_OPENMESH enabled, one can use its capabilities for loading, storing and manipulating meshes. OpenMesh supports many common mesh data formats. We will use an .obj file.

Reading meshes into waLBerla is relatively straightforward, even for large scale simulations. Firstly, we define a block in our parameter file responsible for the mesh and the domain:

meshFile bunny.obj;
domainScaling < 10, 3, 1 >;
dx 3;
cellsPerBlock < 16, 16, 16 >;
periodic < 0, 0, 1 >;

We will comment on the meaning behind these parameters when they are first used in the code snippets.

The name of the file in which our mesh is stored (here: bunny.obj as we want to simulate a flow around the Stanford Bunny) is parsed as usual:

auto domainParameters = walberlaEnv.config()->getOneBlock( "DomainSetup" );
std::string meshFile = domainParameters.getParameter< std::string >( "meshFile" );

Afterwards, the mesh is read in on a single process and broadcasted to all other processes with

auto mesh = make_shared< mesh::TriangleMesh >();
mesh::readAndBroadcast( meshFile, *mesh );

Note that we have requested vertex colors from the mesh for the purpose of coloring the faces of the object according to its vertices in the next step. Therefore we define a helper function vertexToFaceColor that iterates over all faces and colors them in their vertex color. If no uniform coloring of the vertices is given, a default color is taken.

template< typename MeshType >
void vertexToFaceColor( MeshType & mesh, const typename MeshType::Color & defaultColor )
WALBERLA_CHECK( mesh.has_vertex_colors() );
for( auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt )
typename MeshType::Color vertexColor;
bool useVertexColor = true;
auto vertexIt = mesh.fv_iter( *faceIt );
WALBERLA_ASSERT( vertexIt.is_valid() );
vertexColor = mesh.color( *vertexIt );
while( vertexIt.is_valid() && useVertexColor )
if( vertexColor != mesh.color( *vertexIt ) )
useVertexColor = false;
mesh.set_color( *faceIt, useVertexColor ? vertexColor : defaultColor );

After calling this function, we prepare for building the distance octree by precalculating information (e.g. face normals) of the mesh that is required for computing the singed distances from a point to a triangle:

auto triDist = make_shared< mesh::TriangleDistance< mesh::TriangleMesh > >( mesh );

From this information we can finally build the distance octree. It stores information about how close or far boundaries are to each other. Later, this information could be used for e.g. adaptive mesh refinement (note that this will not be covered in this tutorial).

auto distanceOctree = make_shared< mesh::DistanceOctree< mesh::TriangleMesh > >( triDist );

Even though we have successfully loaded the complex geometry and set up the corresponding distance octree, we have not defined our computational LB domain yet. In this tutorial, the LB domain is defined relatively to the loaded geometry. Henceforth, we calculate the axis-aligned bounding box of the geometry and scale it to our needs. Here, we chose our channel to be 10x3x1 times the size of the Stanford Bunny. This scaling is defined in the parameter file (parameter: domainScaling). As the bunny will be placed in the center of the bounding box, we shift the center to the left such that the bunny will be nearer to the inflow.

auto aabb = computeAABB( *mesh );
aabb.scale( domainScaling );
aabb.setCenter( aabb.center() + 0.2 * Vector3<real_t>(aabb.xSize(), 0, 0) );

Finally, we use a convenient built-in data structure that will be responsible for the creation of the structured block forest.

mesh::ComplexGeometryStructuredBlockforestCreator bfc( aabb, Vector3< real_t >( dx ), mesh::makeExcludeMeshInterior( distanceOctree, dx ) );

The ComplexGeometryStructuredBlockforestCreator takes as arguments the axis-aligned bounding box of the domain, the cell sizes dx and an exclusion function. In this tutorial, we want to exclude the interior of the Stanford Bunny with a maximum error of dx. After setting the periodicity obtained from the parameter file

bfc.setPeriodicity( periodicity );

we create the structured block forest on which the sweeps will be performed

auto blocks = bfc.createStructuredBlockForest( cellsPerBlock );

Boundary Handling

The rest of the tutorial mainly follows Tutorial 01 and should therefore be self-explanatory. The only adaption that has to be made is the treatment of the boundaries. Whereas the boundary conditions of the basic domain (inflow, outflow, no-slip) can be again treated by the lbm::DefaultBoundaryHandlingFactory, the no-slip boundary conditions of the bunny need special attention. After the standard routine of the lbm::DefaultBoundaryHandlingFactory

const FlagUID fluidFlagUID( "Fluid" );
auto boundariesConfig = walberlaEnv.config()->getOneBlock( "Boundaries" );
typedef lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory;
BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage( blocks, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID,
boundariesConfig.getParameter< Vector3<real_t> >( "velocity0", Vector3<real_t>() ),
boundariesConfig.getParameter< Vector3<real_t> >( "velocity1", Vector3<real_t>() ),
boundariesConfig.getParameter< real_t > ( "pressure0", real_c( 1.0 ) ),
boundariesConfig.getParameter< real_t > ( "pressure1", real_c( 1.0 ) ) );
geometry::initBoundaryHandling<BHFactory::BoundaryHandling>( *blocks, boundaryHandlingId, boundariesConfig );
geometry::setNonBoundaryCellsToDomain<BHFactory::BoundaryHandling> ( *blocks, boundaryHandlingId );

we address the complex geometry. In a first step, the boundaries that we colored in vertexToFaceColor are equipped with the no-slip boundary UID and the boundary information is added to the mesh.

mesh::ColorToBoundaryMapper<mesh::TriangleMesh> colorToBoundaryMapper(( mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) ));
colorToBoundaryMapper.set( mesh::TriangleMesh::Color(255,255,255), mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) );
auto boundaryLocations = colorToBoundaryMapper.addBoundaryInfoToMesh( *mesh );

In order to have the mesh information available at postprocessing, we write this information to a VTK file as

mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter( mesh, "meshBoundaries", 1 );
meshWriter.addDataSource( make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >( boundaryLocations ) );
meshWriter.addDataSource( make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >() );
meshWriter.addDataSource( make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >() );

Lastly, the boundary information of the complex geometry is added to the structured block forest and the fluid cells and the boundaries are set accordingly.

mesh::BoundarySetup boundarySetup( blocks, makeMeshDistanceFunction( distanceOctree ), numGhostLayers );
boundarySetup.setDomainCells< BHFactory::BoundaryHandling > ( boundaryHandlingId, mesh::BoundarySetup::OUTSIDE );
boundarySetup.setBoundaries<BHFactory::BoundaryHandling>( boundaryHandlingId, makeBoundaryLocationFunction( distanceOctree, boundaryLocations ), mesh::BoundarySetup::INSIDE );


Now that you have seen how to load and process complex geometries from a mesh file, you can play around with this tutorial and adapt it to your needs.

Suggestions for further improvement are the already mentioned adaptive refinement of the loaded geometry. Load balancing of the resulting blocks is then recommended. Thin can be done by using functionalities of mesh::ComplexGeometryStructuredBlockforestCreator and mesh::MeshWorkloadMemory, respectively.