Defining boundary conditions

Overview

The library Palabos currently implements a bunch of different boundary conditions for straight, grid-aligned boundaries. For general boundaries, the available options include stair-cased walls through bounce-back nodes, and higher-order accurate curved boundaries.

IMPORTANT: Curved boundaries are available as of version 1.0, and are part of a fully-featured parallel stack, including parallel voxelization and I/O extensions. Curved boundaries are not yet documented. You can however find a full-featured example in showCases/aneurysm.

Grid-aligned boundaries

The class OnLatticeBoundaryConditionXD

Some of the implemented boundary conditions are local (and are therefore implemented as dynamics classes), some are non-local (and are therefore implemented as data processors), and some have both local and non-local components. To offer a uniform interface to the various possibilities, Palabos offers the interface OnLatticeBoundaryConditionXD which is responsible for instantiating dynamics objects, adding data processors, or both, depending on the chosen boundary condition. The following line for example is used to chose a Zou/He boundary condition:

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition =
        createZouHeBoundaryCondition3D<T,DESCRIPTOR>();

The four possibilities currently offered are (see palabos.org-bc for details):

Regularized BC	createLocalBoundaryConditionXD
Skordos BC	createInterpBoundaryConditionXD
Zou/He BC	createZouHeBoundaryConditionXD
Inamuro BC	createInamuroBoundaryConditionXD

With each of these boundary conditions, it is possible to implement velocity Dirichlet boundaries (all velocity components have an imposed value) and pressure boundaries (the pressure is imposed, and the tangential velocity components are zero). Furthermore, various types of Neumann boundary conditions are proposed. In these case, a zero-gradient for a given variable is imposed with first-order accuracy by copying the value of this variable from a neighboring cell.

Boundaries on a regular block domain

The usage of these boundary conditions is somewhat tricky, because Palabos needs to know if a given wall is a straight wall or a corner (2D), or a flat wall, an edge or a corner (3D), and it needs to know the wall orientation. To simplify this, all OnLatticeBoundaryConditionXD objects have a method setVelocityConditionOnBlockBoundaries and a method setPressureConditionOnBlockBoundaries, to set up boundaries which are located on a rectangular domain. If all nodes on the boundary of a given 2D box boundaryBox implement a Dirichlet condition, use the following command:

Box2D boundaryBox(x0,x1, y0,y1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                           lattice, boundaryBox, locationOfBoundaryNodes );

The first argument indicates the boxed shape on which the boundary is located, and the second argument specifies on which node a velocity condition is to be instantiated. If all boundary nodes are located on the exterior bounding box of the lattice, the above command simply becomes:

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);

Finally, if simply all nodes of the exterior bounding box shall implement a velocity condition, this simplifies to:

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);

Now, suppose that the upper and lower walls, including corners, in a nx -by- ny lattice implement a free-slip condition (zero-gradient for tangential velocity components), the left wall a Dirichlet condition, and the right wall and outflow condition (zero-gradient for all velocity components). This is achieved by adding one more parameter to the function call:

Box2D inlet(0, 0, 2, ny-2);
Box2D outlet(nx-1, nx-1, 2, ny-2);
Box2D bottomWall(0, nx-1, 0, 0);
Box2D topWall(0, nx-1, ny-1, ny-1);

boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, inlet );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, outlet, boundary::outflow );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, bottomWall, boundary::freeslip );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, topWall, boundary::freeslip );

Remember that if the boundary is not located on the outer bounds of the lattice, you must use an additional Box3D argument to say where the boundaries are located, additionally to saying which boundary is going to be added. This is necessary, because Palabos needs to know the orientation and the nature (flat wall, corner, etc.) of the boundary nodes.

Zero-gradient conditions

The optional arguments like boundary::outflow are of type boundary::BcType, and can have the following values for velocity boundaries:

boundary::dirichlet (default value)	Dirichlet condition: imposed velocity value.
boundary::outflow or boundary::neumann	Zero-gradient for all velocity components.
boundary::freeslip	Zero-gradient for tangential velocity components, zero value for the others.
boundary::normalOutflow	Zero-gradient for normal velocity components, zero value for the others.

For pressure boundaries, you have the following choice:

boundary::dirichlet (default value)	Dirichlet condition: imposed pressure value, zero value for the tangential velocity components.
boundary::neumann	Zero-gradient for the pressure, zero value for the tangential velocity components.

Boundary conditions cannot be overriden

It is not possible to override the type of a boundary. Once a boundary node is a Dirichlet velocity node, it stays a Dirichlet velocity node. The result of re-defining it as, say, a pressure boundary node, is undefined. Therefore, boundary conditions must be defined carefully, and piece-wise, as shown in the example above. On the other hand, the value imposed on the boundary (i.e. the velocity value on a Dirichlet velocity boundary) can be changed as often as needed.

Setting the velocity or pressure value on Dirichlet boundaries

A constant velocity value is imposed with the function setBoundaryVelocity. The following line sets the velocity values to zero on all nodes of a 2D domain which have previously been defined to be Dirichlet velocity nodes:

setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );

On all other nodes, this command has no effect. To set a constant pressure value (or equivalently, density value), use the command setBoundaryDensity:

setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. );

A non-constant, space dependent velocity resp. pressure profile can be defined by replacing the velocity resp. pressure argument by either a function or by a function object (an instance of a class which overrides the function call operator) which yields a value of the velocity resp. density as a function of space position. An example is provided in the file examples/showCases/poiseuille/poiseuille.cpp.

Interior and exterior boundaries

Sometimes, the boundary is not located on the surface of a regular, block-shaped domain. In this case, the functions like setVelocityConditionOnBlockBoundaries are useless, and the boundary must be manually constructed. As an example, here’s how to construct manually a free-slip condition along the top-wall, including the corner nodes:

// The automatic creation of a top-wall ...
Box2D topWall(0, nx-1, ny-1, ny-1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               topWall, lattice, boundary::freeslip );

// ... can be replaced by a manual construction as follows:
Box2D topLid(1, nx-2, ny-1, ny-1);
boundaryCondition->addVelocityBoundary1P( topLid, lattice, boundary::freeslip );
boundaryCondition->addExternalVelocityCornerNP( 0, ny-1, lattice, boundary::freeslip );
boundaryCondition->addExternalVelocityCornerPP( nx-1, ny-1, lattice, boundary::freeslip );

A distinction is made between external and internal corners, depending on whether the corner is convex or concave. On the following example geometry, you’ll find five external and one internal corner:

The extensions like 1P`, ``NP, and PP at the end of the methods of the boundary-condition object are used to indicate the orientation of the wall normal, pointing outside the fluid domain, as shown on the figure above. On a straight wall, the code 1P means: “the wall normal points into positive y-direction”. Likewise, the inlet would be labeled with the code 0N as in “negative x-direction”. On a corner, the code NP means “negative x-direction and positive y-direction”. It is important to mention that this wall normal is purely geometrical and does not depend on whether a given wall has the function of an inlet or an outlet. In both cases, the wall normal points away from the fluid, into the non-fluid area.

In 3D, you’ll use the following function for plane walls:

addVelocityBoundaryDO

where the direction D can be 0 for x, 1 for y, or 2 for z, and the orientation O has the value P for positive and N for negative. Edges can be internal or external. For example:

addInternalVelocityEdge0NP
addExternalVelocityEdge1PN

where the first digit of the code indicates the axis to which the edge is parallel, and the two subsequent digits indicate the orientation of the edge inside the plane normal to the edge, in the same way as the corner nodes in 2D. The axes are counted periodically: 0NP means x-plane, negative y-value, and positive z-value, whereas 1PN means y-plane, positive z-value, and negative x-value. Finally, 3D corners are constructed in the same way as in 2D, through a function call like the following:

addInternalVelocityCornerPNP
addExternalVelocityCornerNNN

Periodic boundaries

The concept of periodicity in Palabos is different from the approach chosen for the other types of boundary conditions. Periodicity can only be implemented on opposite, outer boundaries of an atomic-block or a multi-block. This behavior works also if the multi-block has a sparse memory implementation. In this case, all fluid nodes which are in contact with an outer boundary of the multi-block can be periodic, if the corresponding node on the opposite wall is a fluid node.

In the case of an atomic-block (remember that this case is uninteresting, because you should always work with multi-blocks anyway), periodicity is simply a property of the streaming operator. It has no effect for scalar-fields and tensor-fields. In the case of a multi-block however, periodicity can also affect scalar-fields and tensor-fields. Being periodic in this case means that if you define a non-local data processors which accesses nearest neighbor nodes, the value of these neighbor nodes is determined by periodicity along an outer boundary of the multi-block.

By default, all blocks are non-periodic: as mentioned previously, the behavior of outer boundary nodes defaults to one of the versions of bounce-back algorithm. To turn on periodicity along, say, the x-direction, write a command like

lattice.periodicity.toggle(0, true);  // Periodic block-lattice.
scalarField.periodicity.toggle(0, true);  // Periodic scalar-field.

You can also define all boundaries to be periodic through a single command:

lattice.periodicity.toggleAll(true);

Please note that the periodicity or non-periodicity of a block should be defined as soon as possible after the creation of the block, because Palabos needs to know if boundaries are periodic or not when adding or executing data processors.

As a last remark, you should be aware that boundaries of Neumann type (outflow, free-slip, adiabatic thermal wall, etc.) should never be defined to be periodic. This wouldn’t make sense anyway, and it leads to an undefined behavior (aka the program crashes).

Bounce-back

In Palabos, all outer boundaries of a lattice which are not periodic automatically implement a version of the bounce-back boundary algorithm, according to the following algorithm. In pre-collision state, all unknown populations are assigned the post-collision value on the same node at the end of the previous iteration, from the opposite location. We will call this algorithm Bounce-Back 1.

A bounce-back condition can also be used elsewhere in the domain, by explicitly assigning a dynamics object of type BounceBack<T,Descriptor> to chosen cells. On these newly assigned nodes, the collision step is replaced by a bounce-back procedure, during which each population is replaced by the population corresponding to the opposite direction. Such a node can not be considered any more as a fluid node. Computing velocity-moments of the populations on a bounce-back node yields for example arbitrary results, because some of the populations (those which are not incoming from the fluid) have arbitrary values. Consequently, you cannot compute macroscopic quantities like density and velocity on these nodes. Instead, these nodes should be considered as pure algorithmic entities which have the purpose to re-inject into the fluid all populations that leave the domain. Consider fluid nodes adjacent to this type of bounce-back cells. If you think about it, you’ll realize that these nodes behave exactly as if they were cells of type Bounce-Back 1, with a small difference: unknown populations at time t get the post-collision value from opposite populations from time t-2 (and not time t-1 as in Bounce-Back 1. We will refer to this version of bounce-back as Bounce-Back 2.

The effect of both versions of bounce-back is to produce a no-slip wall located half-way between nodes. In the case of Bounce-Back 1, this is half a cell-width beyond the bounce-back fluid cell, whereas in the case of Bounce-Back 2 it is half-way between the last fluid cell and the non-fluid bounce-back node. Obviously, Bounce-Back 2 wastes numerical precision over Bounce-Back 1 by leaping over a time step, and therefore leads to a lower-order accuracy in time of the numerical method. Furthermore, both versions of bounce-back represent the wall location through a staircase approximation, which is low-order accurate (“Order-h”) in space. On the other hand, Bounce-Back 2 offers a huge advantage: its implementation is independent of the wall orientation. If you decide that a node is a bounce-back node, you simply say so; no need to know if it is a corner, a left wall, or a right wall. Constructing a geometry is then as easy as filling areas with bounce-back nodes. In many cases, this ease of use makes up, to a large extent, for the low-order accuracy of the wall representation. When setting up a new simulation, it is always good to try Bounce-Back 2 as a first attempt, as the quality of the obtained results is often surprising. In highly complex geometry such as porous media, bounce-back is often even considered superior to other approaches.

Bounce-back domain from analytical description

The Rule 0 in Palabos can be stated often enough: don’t write loops over space in end-user code. This also true, of course, when assigning bounce-back dynamics to lattice cells. To guarantee reasonable performance, this is task is performed inside a data processor, or even better, by using the wrapper function defineDynamics. This process is illustrated in the example program examples/codesByTopic/bounceBack/instantiateCylinder.cpp. First, a class is written which defines the location of the bounce-back nodes as a function of the cell coordinates iX and iY. In the present case, the nodes are located inside a 2D circular domain:

template<typename T>
class CylinderShapeDomain2D : public DomainFunctional2D {
public:
    CylinderShapeDomain2D(plint cx_, plint cy_, plint radius)
        : cx(cx_),
          cy(cy_),
          radiusSqr(util::sqr(radius))
    { }
    // The function-call operator is overridden to specify the location
    //   of bounce-back nodes.
    virtual bool operator() (plint iX, plint iY) const {
        return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiusSqr;
    }
    virtual CylinderShapeDomain2D<T>* clone() const {
        return new CylinderShapeDomain2D<T>(*this);
    }
private:
    plint cx, cy;
    plint radiusSqr;
};

An instance of this class is then provided as an argument to defineDynamics:

defineDynamics( lattice, lattice.getBoundingBox(),
                new CylinderShapeDomain2D<T>(cx, cy, radius),
                new BounceBack<T,DESCRIPTOR> );

The second argument is of type Box2D, and it can be used to improve the efficiency of this function call: the bounce-back nodes are attributed on cells on the intersection between this box and the domain specified by the domain-functional.

If your domain is specified as the union of a collection of simpler domains, you can use defineDynamics iteratively on each of the simpler domains. If the domain is the intersection or any other geometric operation of simpler domains, you’ll need to play with boolean operators inside the definition of your domain-functional.

Bounce-back domain from boolean mask

In this section, it is assumed that the geometry of the domain is prescribed from an external source. This is for example the case when you simulate a porous media and possess data from a scan of the porous material in question. This easiest approach consists in storing this data as an ASCII file which distinguishes fluid nodes from solid nodes by zeros and ones, separated with spaces, tabs, or newline characters. The file represents the content of a regular array and stores only raw data, no information on the matrix dimensions. The data layout must be the same as in Palabos: an increment of the z-coordinate represents a continuous progress in the file (in 3D), followed by the y-coordinate, and then the x-coordinate. This data is simply flushed from the file into a scalar-field. The following code is taken from the example examples/codesByTopic/io/loadGeometry.cpp:

MultiScalarField2D<T> boolMask(nx, ny);
plb_ifstream ifile("geometry.dat");
ifile >> boolMask;

Note that currently, no error-checking is implemented for such I/O operations. It is your responsibility to ensure that the dimensions of the scalar-field correspond to the size of the data in the geometry file. As a next step, the bounce-back nodes ar instantiated using one of the versions of the function defineDynamics:

defineDynamics(lattice, boolMask, new BounceBack<T,DESCRIPTOR>, true);

This function is currently only defined for multi-blocks using the same data type (T in this case), which explains why the scalar-field boolMask is of type T and not bool, as one could have expected. The scalar-field is not bool, as one could have expected, but the real type T. This is a

Domain creation from an STL file

A full-featured development framework in computational fluid dynamics is expected to be able to read general geometry descriptions (for example the output of a CAD program), and use them to create corresponding boundary conditions in the simulation. In such a general form, this is currently not possible with Palabos. A framework is however offered to read STL files, and to model the corresponding domain by means of bounce-back boundary nodes. The STL file format offers a simple way to describe triangulated surfaces without any additional attributes. It is convenient, because it is widely used, and because most CAD file formats can to some extent be converted to STL. At the same time, the range of application of this format is restricted because of the inability to store additional information in the file, such as the type of the boundary (no-slip, Dirichlet, Neumann, etc.), the value of the velocity in the Dirichlet case, and more. In Palabos, the automatism for the setup of simulations from STL files is therefore restricted to the creation of no-slip walls with bounce-back nodes. Other types of boundaries must still be created manually, as explained below.

To set up a simulation from an STL file, the domain must first be voxelized, i.e. converted from a surface description to a volume description. This is done with help of the open-source software CVMLCPP by Fokko Beekhof. As a first step, you need to install CVMLCPP, which is fortunately an easy thing to do. Simply download the source code of CVMLCPP, and put it into the directory $(PALABOS)/externalLibraries/. In the end, the directory in which the source code of CVMLCPP is located must have the name $(PALABOS)/externalLibraries/cvmlcpp/; rename the directory correspondingly if needed. Don’t compile or install anything: the code of CVMLCPP consists of generic template files which are simply included in Palabos end-user programs. For CVMLCPP to work, you need also to install the open-source library boost, either from the repository of any major Linux distribution, or from the website.

The following explanations are developed around the example of an STL file representing a part of a blood vessel, with a pathological deformation known under the name of aneurysm. It will be shown how to simulate a fluid flow in a domain with the shape prescribe by the domain. Please note that the STL file was designed for the sole purpose of serving as an example in Palabos and does not originate from real-life data. You can use this file and the presented simulations to improve your understanding of Palabos, but in their present form they are inadequate for drawing scientific or medical conclusions.

Full domain creation

The CVMLCPP library can be used to voxelize an STL file with a regular full-domain representation. In this case, the domain is surrounded by a rectangular box, and all spaces inside this box, but outside the fluid domain, are filled with bounce-back node. While this is a simple approach, it presents obvious efficiency issues. The proportion of fluid cells inside this box can be well below 50 percent, and the rest consists of wasted memory. Furthermore, parsing these nodes in the LB code is time consuming, even if their treatment is cheaper than the treatment of normal fluid cells. Finally, the voxelization is not parallelized, and you can easily run into memory issues as a full domain representation is held in the memory of a single processor. For these reasons, a more efficient sparse domain approach is presented in the next section.

An example of a full-domain voxelization is presented in the file examples/codesByTopic/readSTL/fullDomainFromSTL.cpp. The relevant lines of the code are:

ScalarField3D<T>* boolMask = boolMaskFromSTL<T> (
        "aneurysm.stl", parameters.getResolution() );

MultiBlockLattice3D<T,DESCRIPTOR> lattice (
        boolMask->getNx(), boolMask->getNy(), boolMask->getNz(),
        new BGKdynamics<T,DESCRIPTOR>( parameters.getOmega() ) );

// Set up the inlet and outlets here

dynamicsFromBoolMask(*boolMask, lattice, new BounceBack<T,DESCRIPTOR>(1.));

The command boolMaskFromSTL invokes CVMLCPP’s voxelizer to create a scalar field filled with 1’s and 0’s, depending on whether a cell corresponds to the fluid domain or to the solid domain. The data type of this field is T instead of bool, because this is required for instantiating bounce-back nodes from the information in the bool-mask (it is currently not possible to create couplings between blocks based on different data types in Palabos). Subsequently, the function dynamicsFromBoolMask simply invokes a data processor to parse the bool-mask and attribute bounce-back nodes to the lattice correspondingly.

Now comes the less elegant part of the whole process: the creation of inlets and outlets. These must necessarily be defined on a plane surface which is orthogonal to one of the major grid axes, because curved boundaries are not yet implemented in Palabos. The location of the inlet and outlet must be determined manually, for example by inspecting the STL file with a tool like Paraview. It is sufficient to know in which plane the inlet and outlet is located, and have a vague idea of its location in the plane. Then, define inlets and outlets on a rectangular domain which exceeds the fluid domain, to be sure that all of the fluid is covered, and finally erase excess boundary nodes by overwriting them with bounce-back nodes with the function dynamicsFromBoolMask. Similarly, remaining fluid domain behind the inlet, or beyond the outlet, can be erased manually by overriding it with bounce-back nodes. You can imagine bounce-back nodes a bit like the eraser tool in an image manipulation program with which you accurately adjust lines from a first sloppy drawing. This procedure is illustrated on the following figure:

Here, the inlet and the two outlets have been defined with a slightly exceeding domain (gray-colored platelets), and bounce-back nodes have been written in the red-colored boxed to erase before-inlet and after-outlet areas. After this, the function dynamicsFromBoolMask is called to trim the excess in the inlet/outlet regions.

Sparse domain creation

As promised, en efficient sparse domain version of the voxelization process is now presented. It is based on CVMLCPP’s ability to create an octree structure the leaves of which consist of blocks covering the fluid domain. A lot of memory is saved already during the voxelization process (this is important because the voxelization is non-parallel), because only the content of the blocks in the octree is voxelized. Furthermore, no memory is allocated for blocks consisting only of fluid cells. The block structure created in CVMLCPP is then reused in Palabos to create the internal sparse structure of a multi-block. This leads not only to a memory- and CPU-efficient implementation of the simulation, but also to an easy way of parallelizing the program. The blocks created in CVMLCPP being all of equal size, they are simply attributed to the available processors in equal numbers. This is one of these rare situations in CFD in which a very powerful solver is obtained with only simple ingredients, by putting together lattice Boltzmann, bounce-back boundary treatment, CVMLCPP’s octree domain discretization, and Palabos’ parallel multi-block domain representation.

Caution

As of the 0.6 release of Palabos, the CVMLCPP library is using slightly different algorithms for full-domain and sparse-domain voxelization. The two domains are therefore not exactly the same. While it is possible to compare the overall result of the full and the sparse versions, the numerical results cannot be found to be exactly identical.

The following figure shows how the fluid domain is covered by regular blocks in CVMLCPP to obtain a sparse domain representation. The memory requirements are progressively decreased by increasing the depth of the octree decomposition, indicated by the parameter level in the figure:

An example of a sparse-domain voxelization is presented in the file examples/codesByTopic/readSTL/sparseDomainFromSTL.cpp. The relevant lines of the code are:

int limitDepth = 2;
TreeReader3D<int>* treeReader =
    createUniformTreeReader3D("aneurysm.stl", parameters.getNx(), limitDepth);

MultiBlockLattice3D<T,DESCRIPTOR>* lattice = generateSparseMultiBlock (
        *treeReader,
        new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

// Set up the inlet and outlets here

bool wallFlag = false;
dynamicsFromOctree(*lattice, *treeReader, new BounceBack<T,DESCRIPTOR>(1.), wallFlag);

Choosing the level of the octree decomposition has an important impact on the efficiency of the simulation. A small level leads to a decomposition in large blocks, which is penalizing by the presence of many unnecessary bounce-back nodes, while a large level leads to an abundance of small blocks, which ultimately impacts the efficiency because of the communication overhead between blocks. On a parallel machine, the level must also be chosen large enough to have enough blocks to assign to each processor.

The treatment of the inlet and the outlets is exactly identical in the sparse-domain and in the full-domain case. Actually, the two codes fullDomainFromSTL and sparseDomainFromSTL are identical except for the creation of the multi-block and the instantiation of the bounce-back wall nodes from the STL file.