Tutorial 1: First steps with Palabos

Tutorial 1.1: The first code

To get started, download the tarball of the most recent version of Palabos, for example palabos-1.0r0.tgz, and unpack it with the command tar xvfz palabos-1.0r0.tgz. At this point, the code of the library is not explicitly compiled and installed. Instead, the library is compiled on-demand when end-user applications are created. This approach reflects the fact that many Palabos users are also developers, and end-user development is coupled with development of the source core. Furthermore, the approach makes it easy to recompile the code on-demand with or without debug flags, or for sequential/parallel execution.

Change into the directory of the tutorial, palabos/examples/tutorial/section_1, and type make to compile the library and create the executable of the first tutorial code, tutorial_1_1.cpp. Finally, the application is executed by typing ./tutorial_1_1 at the command line.

This application simulates the time evolution of an initial-value problem, without boundary conditions. The boundaries are chosen to be periodic, i.e. flow particles leaving the domain on one boundary re-enter the domain on the opposite boundary. The initial condition has a zero velocity and constant density. On a squared sub-domain of the field, the density is initialized at a slightly higher value, to create a perturbation which then generates a flow pattern.

This setup is chosen mostly to illustrate programming concepts in Palabos, and not to propose an interesting hydrodynamic problem. It should in particular be pointed out that the initial condition does not represent the state of an incompressible fluid, because the zero velocity is incompatible with a non-zero density. After sufficient iterations, the flow automatically converges to a situation which is compatible with the physics of an incompressible, or slightly compressible, flow.

The code tutorial_1_1.cpp is listed below:

/* Code 1.1 in the Palabos tutorial
 */

#include "palabos2D.h"
#include "palabos2D.hh"
#include <iostream>
#include <iomanip>

using namespace plb;
using namespace std;


typedef double T;
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor


// Initialize the lattice at zero velocity and constant density, except
//   for a slight density excess on a square sub-domain.
void defineInitialDensityAtCenter(MultiBlockLattice2D<T,DESCRIPTOR>& lattice)
{
    // The lattice is of size nx-by-ny
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();

    // Create a Box2D which describes the location of cells with a slightly
    //   higher density.
    plint centralSquareRadius = nx/6;
    plint centerX = nx/3;
    plint centerY = ny/4;
    Box2D centralSquare (
            centerX - centralSquareRadius, centerX + centralSquareRadius,
            centerY - centralSquareRadius, centerY + centralSquareRadius );

    // All cells have initially density rho ...
    T rho0 = 1.;
    // .. except for those in the box "centralSquare" which have density
    //    rho+deltaRho
    T deltaRho = 1.e-4;
    Array<T,2> u0(0,0);

    // Initialize constant density everywhere.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), rho0, u0 );

    // And slightly higher density in the central box.
    initializeAtEquilibrium (
           lattice, centralSquare, rho0 + deltaRho, u0 );

    lattice.initialize();
}

int main(int argc, char* argv[]) {
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    const plint maxIter = 1000; // Iterate during 1000 steps.
    const plint nx = 600;       // Choice of lattice dimensions.
    const plint ny = 600;
    const T omega = 1.;        // Choice of the relaxation parameter

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
           nx, ny, new BGKdynamics<T,DESCRIPTOR>(omega) );

    lattice.periodicity().toggleAll(true); // Use periodic boundaries.

    defineInitialDensityAtCenter(lattice);

    // Main loop over time iterations.
    for (plint iT=0; iT<maxIter; ++iT) {
        if (iT%40==0) {  // Write an image every 40th time step.
            pcout << "Writing GIF file at iT=" << iT << endl;
            // Instantiate an image writer with the color map "leeloo".
            ImageWriter<T> imageWriter("leeloo");
            // Write a GIF file with colors rescaled to the range of values
            //   in the matrix
            imageWriter.writeScaledGif (
                    createFileName("u", iT, 6),
                    *computeVelocityNorm(lattice) );
        }
        // Execute lattice Boltzmann iteration.
        lattice.collideAndStream();
    }
}

Global definitions and includes

Line 4-7

The include file palabos2D.h gives access to all declarations in the Palabos project. Additionally, the file palabos2D.hh is included to guarantee access to the full template code, for example to instantiate simulations with different data types than the standard data type double. The distinction between generic and non-generic code is explained in the next tutorial.

Line 11-12

These two declarations guarantee an automatic access to the names in the Palabos code, which is contained in the namespace plb, and to the C++ standard library, which is contained in the namespace std.

Line 14-15

As specified by these two directives, the simulation will be executed with double precision floating point numbers, and on the two-dimensional D2Q9 lattice.

Creating the initial condition

Line 19

The data type at the heart of Palabos is the MultiBlockLatticeXD, where X stands for 2 or 3, for 2D or 3D simulations. As will be shown shortly, the word Multi stands for the fact that in the internal implementation, the lattice is often decomposed into several smaller lattices. On the interface however, and this abstraction mechansim can be confidently used, in a first time at least, to ignore details of the behind-the-scene action.

Line 27-32

The domain is initialized with a slightly exceeding density on a rectangular sub-domain of the full domain. This domain is defined through a geometric object of type Box2D.

Line 42

Initialize all lattice cells at an equilibrium distribution with constant density and zero velocity.

Line 46

Then, redefine the initial values for the cell inside the box centralSquare to have a density exceeding the other ones by deltaRho.

Line 49

The method initialize is called after initialization operations, to prepare the lattice for simulations.

Running the simulation

Line 53

The function plbInit must mandatorily be invoked in the very beginning of the program, to guarantee consistent results between sequential and parallel program runs.

Line 54

Specify the directory where output files will be placed.

Line 61-62

During creation of the multi block, the default dynamics (in this example BGK), also called background dynamics, is specified. This object defines the nature of the collision to be executed by all cells, unless a new individual dynamics is attributed to them in the simulation setup.

Line 64

Periodic boundaries are handled in a special way in Palabos. They can only be imposed on the outer bound of a block-lattice, through a call to the method periodicity().

Line 73

An image writer object provides a convenient way to write GIF images from simulated data, in 2D or in 3D. In the latter case, the image represents for example a slice through the computational domain. Such images are mostly used to monitor the evolution of the simulation, and to get a qualitative idea of the state of the simulation, before proceeding to a detailed evaluation of the data. To check the images, change into the directory tmp during or after the simulation, and visualize them for example with help of the command display. To get an animation from subsequent GIF images, use the command convert (convert -delay 5 u*.gif animation.gif), and visualize the animated gif through a command like animate animation.gif.

Line 81

The method collideAndStream executes a collision and a streaming step on each cell of the lattice. The two steps can also be separated by calling first lattice.collide() and then lattice.stream(). The synchronous execution of collision and streaming is however more efficient, because it requires only a single traversal of the data space to execute a full lattice Boltzmann iteration. The argument true to the method collideAndStream or to the method stream is used to implement periodicity for the domain boundaries. Reversely, the argument false is used for non-periodic boundaries which then implement, by default, a half-way bounce-back condition. The effect of this condition is to simulate no-slip walls located half a cell spacing beyond the outer lattice sites.

This tutorial shows that the central structure in which the data of a simulation is stored, is the MultiBlockLatticeXD. It was chosen this way, because the interface of the MultiBlockLatticeXD mimicks regular data arrays with which researchers and engineers are usually familiar from the implementation of simpler problems in languages like Matlab or Fortran. Another advantage of the “regular-block interface” is the ease with which individual lattice cells are identified and treated, for example for the implementation of boundary conditions.

If you are unfamiliar with object-oriented programming, it is likely that at this point you start feeling nervous about this form of data encapsulation and the apparent loss of control over the details of the program workflow. If you stay focused, you will however progressively understand that the behind-the-scene action is simple enough and by no means hidden from the understanding of the programmer. Instead, object-oriented mechanisms make it easier to keep control over a structured and well understood course of action. The principles and mechanisms of the MultiBlockLatticeXD are presented and trained in the Tutorial 2. In particular, Tutorial 2.1 shows that the MultiBlockLatticeXD could in practice be replaced by an AtomicBlockLatticeXD, which stands for a simple regular array, and you should feel free to go ahead and do so. Be aware though that there are in practice only disadvantages in making this substitution, as one loses for example the possibility to parallelize the program. On the other hand, it has however been observed that the use of the simple AtomicBlockLatticeXD can help overcome psychological barriers by providing a sense of control in a first stage of gaining familiarity with Palabos.

Tutorial 1.2: Initializing the lattice

In the previous lesson, the initial condition was created by assigning a constant value of density and velocity to domains of rectangular shape. A more detailed initialization of lattice cells is obtained by writing functions in which a different value of density and velocity is attributed to each cell. This is illustrated in the tutorial code 1.2, tutorial_1_2.cpp. To compile this code, edit the file Makefile, and replace tutorial_1_1 by tutorial_1_2 in the corresponding line.

The relevant parts of the code tutorial_1_2.cpp are listed below:

const plint maxIter = 1000; // Iterate during 1000 steps.
const plint nx = 600;       // Choice of lattice dimensions.
const plint ny = 600;
const T omega = 1.;        // Choice of the relaxation parameter

T rho0 = 1.; // All cells have initially density rho ...
// .. except for those inside the disk which have density
//    rho+deltaRho
T deltaRho = 1.e-4;
Array<T,2> u0(0,0);

void initializeConstRho(plint iX, plint iY, T& rho, Array<T,2>& u) {
    u = u0;
    rho = rho0 + deltaRho;
}

void initializeRhoOnDisk(plint iX, plint iY, T& rho, Array<T,2>& u) {
    plint radius = nx/6;
    plint centerX = nx/3;
    plint centerY = ny/4;
    u = u0;
    if( (iX-centerX)*(iX-centerX) + (iY-centerY)*(iY-centerY) < radius*radius) {
        rho = rho0 + deltaRho;
    }
    else {
        rho = rho0;
    }
}

// Initialize the lattice at zero velocity and constant density, except
//   for a slight density excess on a circular sub-domain.
void defineInitialDensityAtCenter(MultiBlockLattice2D<T,DESCRIPTOR>& lattice)
{
    // Initialize constant density everywhere.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), rho0, u0 );

    // And slightly higher density in the central box.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), initializeRhoOnDisk );

    lattice.initialize();
}

Line 12

The function initializeConstRho can be used to obtain exactly the same effect as in the previous lesson: assign a constant density to each cell of a sub-domain.

Line 17

The function initializeRhoOnDisk on the other hand attributes a different density only to cells contained inside a disk. The initial condition is therefore more natural, having a circular instead of a rectangular shape.

Line 39

This time, the function initializeAtEquilibrium is called with the function initializeRhoOnDisk as an argument, instead of a constant value of density and velocity. More generally, the function initializeAtEquilibrium behaves like an algorithm of the C++ standard template library. It’s argument can be a functional in the general sense, i.e. either a classical function or an object which behaves like a function by overloading the function call operator. The advantage of using an object instead of a function is that objects can store internal data and can therefore have a customized behavior. In such a case, the function initializeRhoOnDisk could for example be customized with respect to the center and radius of the disk. The usage of such a function object is illustrated in the tutorial code 1.4.

Tutorial 1.3: Compilation options

So far, we have compiled the programs using the default compiler options. Edit the file Makefile to see other available options

Precompiled mode

Obsolete: the precompiled mode is not available any more in Version 0.7 or newer.

Debug mode

If your program crashes, you may want to recompile in debug mode, with the flag debug=true. In this mode, the program executes additional error checking, which can help in locating the source of the error. Furthermore, debug information is put into the object file, after which you can use a debugger such as gdb.

Parallel mode

All codes developed in this tutorial work in serial and parallel. Be aware, though, that most codes execute massively output operations, like the generation of GIF images. They are therefore unlikely to run significantly faster in parallel, unless you comment out the output operations. To compile for parallel execution, choose your C++ compilation command which adds MPI functionality, and define it in the line parallelCXX=.... Then, set the flag MPIparallel=true. Finally, execute your program with the command mpirun or mpiexec or similar.

Tutorial 1.4: Data analysis

In this lesson, a few approaches to analyzing the results of a simulation are reviewed. The corresponding code is found in the file tutorial_1_4.cpp

Let us look at the function main in this code:

int main(int argc, char* argv[]) {
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
           nx, ny, new BGKdynamics<T,DESCRIPTOR>(omega) );

    lattice.periodicity().toggleAll(true); // Set periodic boundaries.

    defineInitialDensityAtCenter(lattice);

    // First part: loop over time iterations.
    for (plint iT=0; iT<maxIter; ++iT) {
        lattice.collideAndStream();
    }

    // Second part: Data analysis.
    Array<T,2> velocity;
    lattice.get(nx/2, ny/2).computeVelocity(velocity);
    pcout << "Velocity in the middle of the lattice: ("
          << velocity[0] << "," << velocity[1] << ")" << endl;

    pcout << "Velocity norm along a horizontal line: " << endl;
    Box2D line(0, 100, ny/2, ny/2);
    pcout << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

    plb_ofstream ofile("profile.dat");
    ofile << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

    pcout << "Average density in the domain: " << computeAverageDensity(lattice) << endl;
    pcout << "Average energy in the domain: " << computeAverageEnergy(lattice) << endl;
}

Line 19

The method get delivers access to a lattice cell. A cell is an object with which one can not only access the actual variables defined on the cell, but also perform useful computations, such as, evaluate macroscopic variables. In this example, the velocity is computed on a cell in the middle of the domain.

Line 25

With the method extractSubDomain, a rectangular sub domain of the full lattice is extracted. In the present example, the velocity norm is computed on the extracted domain and printed to the terminal.

Line 27

In the same way, C++ streams are used to write the data into a file instead of the terminal. Make sure to use the data type plb_ofstream instead of ofstream to guarantee working conditions in parallel programs.

Line 30-31

Many functions are predefined for computing average values, like the average density or the average kinetic energy in the present example.

The data written to the file profile.dat is conveniently post-processed with a mathematics processing tool. For example, type the two following lines at the command prompt of the program Octave to plot a curve of the velocity profile:

load profile.dat
plot(profile)

Tutorial 1.5: Boundary conditions

For the first time, we implement a simulation which represents a well-known physical situation: a Poiseuille flow. This flow evolves in a channel with no-slip walls, and the flow velocity is everywhere parallel to the channel walls. The analytical solution of this flow is represented by a parabolic profile for the velocity component parallel to the channel walls. In the simulation, the analytical solution is used to implement Dirichlet boundary condition for the velocity on the domain inlet and outlet. As an initial condition, a zero velocity field is used, and the simulation is started to converge towards the Poiseuille solution in each point of the domain.

Note that the simulation has initially a discontinuity between the zero velocity in the initial condition and the finite velocity in the boundary condition. If the simulation becomes numerically unstable when you play with the parameters, try either increasing the lattice resolution or decreasing the velocity in lattice units.

This time, the full code is reprinted in the tutorial:

#include "palabos2D.h"
#include "palabos2D.hh"

#include <vector>
#include <iostream>
#include <iomanip>

/* Code 1.5 in the Palabos tutorial
 */

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor


/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const& parameters) {
    T y = (T)iY / parameters.getResolution();
    return 4.*parameters.getLatticeU() * (y-y*y);
}

/// A functional, used to initialize the velocity for the boundary conditions
template<typename T>
class PoiseuilleVelocity {
public:
    PoiseuilleVelocity(IncomprFlowParam<T> parameters_)
        : parameters(parameters_)
    { }
    /// This version of the operator returns the velocity only,
    ///    to instantiate the boundary condition.
    void operator()(plint iX, plint iY, Array<T,2>& u) const {
        u[0] = poiseuilleVelocity(iY, parameters);
        u[1] = T();
    }
    /// This version of the operator returns also a constant value for
    ///    the density, to create the initial condition.
    void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
        u[0] = poiseuilleVelocity(iY, parameters);
        u[1] = T();
        rho  = (T)1;
    }
private:
    IncomprFlowParam<T> parameters;
};

void channelSetup (
        MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
        IncomprFlowParam<T> const& parameters,
        OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
    // Create Velocity boundary conditions.
    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

    // Specify the boundary velocity.
    setBoundaryVelocity (
            lattice, lattice.getBoundingBox(),
            PoiseuilleVelocity<T>(parameters) );

    // Create the initial condition.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );

    lattice.initialize();
}

void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
    const plint imSize = 600;
    ImageWriter<T> imageWriter("leeloo");
    imageWriter.writeScaledGif(createFileName("u", iter, 6),
                               *computeVelocityNorm(lattice),
                               imSize, imSize );
}

int main(int argc, char* argv[]) {
    plbInit(&argc, &argv);

    global::directories().setOutputDir("./tmp/");

    // Use the class IncomprFlowParam to convert from
    //   dimensionless variables to lattice units, in the
    //   context of incompressible flows.
    IncomprFlowParam<T> parameters(
            (T) 1e-2,  // Reference velocity (the maximum velocity
                       //   in the Poiseuille profile) in lattice units.
            (T) 100.,  // Reynolds number
            100,       // Resolution of the reference length (channel height).
            2.,        // Channel length in dimensionless variables
            1.         // Channel height in dimensionless variables
    );
    const T imSave   = (T)0.1;  // Time intervals at which to save GIF
                                //   images, in dimensionless time units.
    const T maxT     = (T)3.1;  // Total simulation time, in dimensionless
                                //   time units.

    writeLogFile(parameters, "Poiseuille flow");

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
             parameters.getNx(), parameters.getNy(),
             new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

    OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
        boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

    channelSetup(lattice, parameters, *boundaryCondition);

    // Main loop over time iterations.
    for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
        if (iT%parameters.nStep(imSave)==0 && iT>0) {
            pcout << "Saving Gif at time step " << iT << endl;
            writeGifs(lattice, iT);
        }
        // Execute lattice Boltzmann iteration.
        lattice.collideAndStream();
    }

    delete boundaryCondition;
}

Poiseuille profile

Line 20

This function computes the parabolic Poiseuille profile, in a channel of a given height, and with a given maximum velocity in the middle of the channel. The object IncomprFlowParam is discussed below: it stores various parameters of the simulation.

Line 26-47

This is an example for a so-called function object, or functional. It is a class which overloads the function call operator, operator(). Instances of this class behave like normal functions, but the object is more intelligent than a pure function. It can for example accept parameters at the constructor and possess a configurable behavior. The present function object encapsulates the function poiseuilleVelocity and is configured with an object of type IncomprFlowParam. This is a useful trick to avoid the need for declaring simulation parameters as global variables and making them accessible to everyone, as we did in the code 1.2.

Line 55

Specify that all exterior boundaries of the domain implement a Dirichlet boundary condition for the velocity. At this point, the value of the velocity on the boundaries is yet undefined.

Line 58

Define the value of the velocity on boundary nodes from the analytical Poiseuille profile. Note that although this function is applied to the entire domain, it has an effect only of nodes which have been previously defined as being Dirichlet boundary nodes.

Simulation parameters and boundary condition

Line 86

An object of type IncomprFlowParam stores the parameters of the simulation (Reynolds number, lattice resolution, reference velocity in lattice units, etc.), and performs unit conversions. It is for example used to compute automatically the relaxation parameter omega from the Reynolds number.

Line 105

Choose the algorithm which is used to implement the boundary condition. The type LocalBoundaryCondition implements the regularized boundary condition, which is entirely local (it doesn’t need to acces neighbor nodes). The type InterpBoundaryCondition uses finite difference schemes to compute the strain rate tensor on walls, and is therefore non-local.

Tutorial 1.6: 3D Channel flow

To keep everything as simple as possible, all tutorials so far were based on 2D simulations. Turning to 3D is pretty easy, though. Let’s reconsider the Poiseuille flow from the previous tutorial and implement in 3D. The resulting code, shown below, can also be found in the file lesson_1_6.cpp:

#include "palabos3D.h"
#include "palabos3D.hh"
#include <vector>
#include <iostream>
#include <iomanip>


/* Code 1.6 in the Palabos tutorial
 */

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D3Q19Descriptor


/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, plint iZ, IncomprFlowParam<T> const& parameters) {
    T y = (T)iY / parameters.getResolution();
    T z = (T)iZ / parameters.getResolution();
    return 4.*parameters.getLatticeU() * (y-y*y) * (z-z*z);
}

/// A functional, used to initialize the velocity for the boundary conditions
template<typename T>
class PoiseuilleVelocity {
public:
    PoiseuilleVelocity(IncomprFlowParam<T> parameters_)
        : parameters(parameters_)
    { }
    void operator()(plint iX, plint iY, plint iZ, Array<T,3>& u) const {
        u[0] = poiseuilleVelocity(iY, iZ, parameters);
        u[1] = T();
        u[2] = T();
    }
private:
    IncomprFlowParam<T> parameters;
};

void channelSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
                   IncomprFlowParam<T> const& parameters,
                   OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
    // Create Velocity boundary conditions
    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

    setBoundaryVelocity (
            lattice, lattice.getBoundingBox(),
            PoiseuilleVelocity<T>(parameters) );

    lattice.initialize();
}

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
    const plint imSize = 600;
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();
    const plint nz = lattice.getNz();
    Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);
    ImageWriter<T> imageWriter("leeloo");
    imageWriter.writeScaledGif(createFileName("u", iter, 6),
                               *computeVelocityNorm(lattice, slice),
                               imSize, imSize );
}

int main(int argc, char* argv[]) {
    plbInit(&argc, &argv);

    global::directories().setOutputDir("./tmp/");

    // Use the class IncomprFlowParam to convert from
    //   dimensionless variables to lattice units, in the
    //   context of incompressible flows.
    IncomprFlowParam<T> parameters(
            (T) 1e-2,  // Reference velocity (the maximum velocity
                       //   in the Poiseuille profile) in lattice units.
            (T) 100.,  // Reynolds number
            30,        // Resolution of the reference length (channel height).
            3.,        // Channel length in dimensionless variables
            1.,        // Channel height in dimensionless variables
            1.         // Channel depth in dimensionless variables
    );
    const T imSave   = (T)0.02; // Time intervals at which to save GIF
                                //   images, in dimensionless time units.
    const T maxT     = (T)2.5;  // Total simulation time, in dimensionless
                                //   time units.

    writeLogFile(parameters, "3D Poiseuille flow");

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

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
        //boundaryCondition = createInterpBoundaryCondition3D<T,DESCRIPTOR>();
        boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

    channelSetup(lattice, parameters, *boundaryCondition);

    // Main loop over time iterations.
    for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
        if (iT%parameters.nStep(imSave)==0) {
            pcout << "Saving Gif at time step " << iT << endl;
            writeGifs(lattice, iT);
        }
        // Execute lattice Boltzmann iteration.
        lattice.collideAndStream();
    }

    delete boundaryCondition;
}

The few lines which have changed from the 2D to the 3D code are the following:

Line 1 and 3

The header files palabos2D.h and palabos2D.hh are replaced by their 3D counterpart.

Line 16

On this line, the type of lattice is chosen. While in 2D the only available nearest-neighbor lattice is D2Q9 (D2Q7 is currently not supported), there are 3 choices in 3D: D3Q15, D3Q19, and D3Q27.

Line 19

To set up a Dirichlet condition for the inlet and the outlet, an analytic velocity profile is proposed. There exists an analytical solution for 3D channel flow, but we do not use it here to keep the code simple. Instead, the velocity profile is formulated as a tensor product of the 2D Poiseuille flow, once in each of the two axes parallel to the inlet and outlet.

Line 26-39

The function object which is used to instantiate the boundary conditions has the same form in 3D as in 2D. This time, the function call operator takes three arguments, the three space dimensions, and returns a 3D velocity vector.

Line 41

Initialization of the geometry is exactly identical as in the previous 2D Poiseuille flow, except that the extension 2D is replaced by 3D in all keywords.

Line 55

It remains useful to write regularly GIF snapshots of the flow to monitor the flow evolution in a simple way. A 2D slice needs to be extraced from the full domain to produce an image.

Line 76

It is possible to provide the class IncomprFlowParam with an additional argument to indicate the extent of the domain in z-direction in the 3D case.

Line 92

The central object of the simulation, the MultiBlockLattice2D, is replaced by a MultiBlockLattice3D.

Tutorial 1.7: Post-processing with Paraview

In previous tutorials, two approaches to analyzing the computed numerical data have been used. The first consists of using the class ImageWriter to produce GIF snapshots of the velocity norm, or other scalar flow variables, during the simulation. This approach is particularly useful for obtaining a qualitative impression of the evolution of the flow, and for being able to interrupt the simulation when errors appear. Another, more quantitative approach, consists of writing the data into a file, using a plain ASCII format. In this case, the data stream is linearized, which means that structural information, such as the size of the domain in each space dimension, is lost. In the following code, based on the previous tutorial, the velocity is for example computed on a slice perpendicular to the z-axis, and written into an ASCII file:

// Attribute a value to nx, ny, and nz.
Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);
plb_ofstream ofile("slice.dat");
ofile << setprecision(4) << *computeVelocityNorm(lattice, slice) << endl;

As you can see, the C++ stream operators have been overloaded for this task. This means that the output can be formated using I/O manipulators like setprecision(n) to chose the numerical precision, setw(n) to format each number in a cell of fixed width, or fixed respectively scientific to chose the representation of floating point variables. The flow field stored in the file slice.dat can for example be analyzed with help of the open-source program Octave. For this, you need to know the values of nx, ny and nz, because they are not stored in the file:

% Matlab/Octave script
% Assign a value to nx and ny
load slice.dat
slice = reshape(slice, nx, ny);
imagesc(slice)

In the present tutorial, an additional approach is shown in which the data is written into a file in a so-called VTK format, after which it can be analyzed with a general post-processing tool. The following function is used in the tutorial code tutorial_1_7 to write VTK data of the 3D Poiseuille flow presented in tutorial 1.6:

void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
              IncomprFlowParam<T> const& parameters, plint iter)
{
    T dx = parameters.getDeltaX();
    T dt = parameters.getDeltaT();
    VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);
    vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);
    vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt);
}

The data in the VTK file is represented in dimensionless variables. For this, the variable dx is given to the constructor of the VtkImageOutput3D object in order to indicate the size of a cell spacing. Furthermore, every variable which is written needs to be rescaled according to its units, namely 1. for the density and dx/dt for the velocity. A few comments on the VTK output are in order:

• The VTK files can hold an arbitrary number of fields, which are written one after another with the command writeData. In case of Vector fields, the number of components of the vector needs to be indicated as a template argument.
• The data is written in binary format, using the binary-encoded ASCII format Base64 for compatibility. Therefore, unlike the plain ASCII output used previously, the VTK format does not loose any digit of precision.
• In the above example, the data is converted from double to float to save space on the disk.

The open-source program Paraview provides a convenient interactive interface for analyzing the data. To end this tutorial, you are encouraged to run the program tutorial_1_7, install Paraview on your system and analyze the output. Among others, Paraview can easily produce an animation from the written image series.