Palabos Tutorial 1/3: First steps with Palabos

Introductory comments

The following tutorial provides an overview of the C++ programmer interface of the Palabos library. If instead you are interested in the Python interface, have a look at the Palabos-Python tutorial from DSFD 2010 (PDF). As for the Java interface, no tutorial is available at this moment. The user’s guide provides however installation instructions, and you can then have a look at the provided example applications.

Tutorial 1: First steps with Palabos

Tutorial 1.1: The first code

To get started under Linux or Mac OS X you need to have the C++ frontend of GCC (g++) installed. Under Mac OS X, a convenient way to get GCC is to install xcode. If any of the examples fails to compile under Mac OS X, edit the Makefile, and add the option -DPLB_MAC_OS_X to the entry compileFlags = ....

Note

Windows programmers

Under Windows, you can get started with Palabos in two different ways: either you install the Code::Blocks (codeblocks) programming environment (choose the download that is packaged with the MinGW library), or you run Linux in a virtual machine. The first approach is easier and sufficient for the purposes of this tutorial. The second approach is however more convenient in the long run, because currently, not all features of Palabos (example: the Python and the Java bindings) are available under Windows.

If you choose to work with Code::Blocks, the following command-line instructions do not apply. Instead, use the Code::Blocks IDE and create a project for each code of the tutorial.

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. To use Palabos, you do not need to proceed with an explicit compilation and installation of the library. 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. Finally, it makes it simple to use Palabos on a machine on which you have no administrator rights.

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:

1. /* Code 1.1 in the Palabos tutorial

2. */

3.

4. #include "palabos2D.h"

5. #include "palabos2D.hh"

6. #include <iostream>

7. #include <iomanip>

8.

9. using namespace plb;

10. using namespace std;

11.

12.

13. typedef double T;

14. #define DESCRIPTOR plb::descriptors::D2Q9Descriptor

15.

16.

17. // Initialize the lattice at zero velocity and constant density, except

18. // for a slight density excess on a square sub-domain.

19. void defineInitialDensityAtCenter(MultiBlockLattice2D<T,DESCRIPTOR>& lattice)

20. {

21. // The lattice is of size nx-by-ny

22. const plint nx = lattice.getNx();

23. const plint ny = lattice.getNy();

24.

25. // Create a Box2D which describes the location of cells with a slightly

26. // higher density.

27. plint centralSquareRadius = nx/6;

28. plint centerX = nx/3;

29. plint centerY = ny/4;

30. Box2D centralSquare (

31. centerX - centralSquareRadius, centerX + centralSquareRadius,

32. centerY - centralSquareRadius, centerY + centralSquareRadius );

33.

34. // All cells have initially density rho ...

35. T rho0 = 1.;

36. // .. except for those in the box "centralSquare" which have density

37. // rho+deltaRho

38. T deltaRho = 1.e-4;

39. Array<T,2> u0(0,0);

40.

41. // Initialize constant density everywhere.

42. initializeAtEquilibrium (

43. lattice, lattice.getBoundingBox(), rho0, u0 );

44.

45. // And slightly higher density in the central box.

46. initializeAtEquilibrium (

47. lattice, centralSquare, rho0 + deltaRho, u0 );

48.

49. lattice.initialize();

50. }

51.

52. int main(int argc, char* argv[]) {

53. plbInit(&argc, &argv);

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

55.

56. const plint maxIter = 1000; // Iterate during 1000 steps.

57. const plint nx = 600; // Choice of lattice dimensions.

58. const plint ny = 600;

59. const T omega = 1.; // Choice of the relaxation parameter

60.

61. MultiBlockLattice2D<T, DESCRIPTOR> lattice (

62. nx, ny, new BGKdynamics<T,DESCRIPTOR>(omega) );

63.

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

65.

66. defineInitialDensityAtCenter(lattice);

67.

68. // Main loop over time iterations.

69. for (plint iT=0; iT<maxIter; ++iT) {

70. if (iT%40==0) { // Write an image every 40th time step.

71. pcout << "Writing GIF file at iT=" << iT << endl;

72. // Instantiate an image writer with the color map "leeloo".

73. ImageWriter<T> imageWriter("leeloo");

74. // Write a GIF file with colors rescaled to the range of values

75. // in the matrix

76. imageWriter.writeScaledGif (

77. createFileName("u", iT, 6),

78. *computeVelocityNorm(lattice) );

79. }

80. // Execute lattice Boltzmann iteration.

81. lattice.collideAndStream();

82. }

83. }

Global definitions and includes

Line 4-7The 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-12These 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-15As 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 19The 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-32The 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 42Initialize all lattice cells at an equilibrium distribution with constant density and zero velocity.Line 46Then, redefine the initial values for the cell inside the box centralSquare to have a density exceeding the other ones by deltaRho.Line 49The method initialize is called after initialization operations, to prepare the lattice for simulations.

Running the simulation

Line 53The function plbInit must mandatorily be invoked in the very beginning of the program, to guarantee consistent results between sequential and parallel program runs.Line 54Specify the directory where output files will be placed.Line 61-62During 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 64Periodic 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 73An 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 81The 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 MultiBlockLatticeXDcould 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:

1. const plint maxIter = 1000; // Iterate during 1000 steps.

2. const plint nx = 600; // Choice of lattice dimensions.

3. const plint ny = 600;

4. const T omega = 1.; // Choice of the relaxation parameter

5.

6. T rho0 = 1.; // All cells have initially density rho ...

7. // .. except for those inside the disk which have density

8. // rho+deltaRho

9. T deltaRho = 1.e-4;

10. Array<T,2> u0(0,0);

11.

12. void initializeConstRho(plint iX, plint iY, T& rho, Array<T,2>& u) {

13. u = u0;

14. rho = rho0 + deltaRho;

15. }

16.

17. void initializeRhoOnDisk(plint iX, plint iY, T& rho, Array<T,2>& u) {

18. plint radius = nx/6;

19. plint centerX = nx/3;

20. plint centerY = ny/4;

21. u = u0;

22. if( (iX-centerX)*(iX-centerX) + (iY-centerY)*(iY-centerY) < radius*radius) {

23. rho = rho0 + deltaRho;

24. }

25. else {

26. rho = rho0;

27. }

28. }

29.

30. // Initialize the lattice at zero velocity and constant density, except

31. // for a slight density excess on a circular sub-domain.

32. void defineInitialDensityAtCenter(MultiBlockLattice2D<T,DESCRIPTOR>& lattice)

33. {

34. // Initialize constant density everywhere.

35. initializeAtEquilibrium (

36. lattice, lattice.getBoundingBox(), rho0, u0 );

37.

38. // And slightly higher density in the central box.

39. initializeAtEquilibrium (

40. lattice, lattice.getBoundingBox(), initializeRhoOnDisk );

41.

42. lattice.initialize();

43. }

Line 12The 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 17The 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 39This 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

Debug mode

Debug mode is activated in Palabos by setting 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 like gdb.

Debug mode is activated by default in all Palabos examples. As a general rule, we recommend that you keep this flag activated even in production mode. It decreases the overall performance of Palabos by a few percents only, and provides helpful insights whenever your program crashes.

Parallel mode

All codes presented in this tutorial work in serial and parallel. Be aware, though, that most examples perform lots of output operations, such as, the generation of GIF images. They are therefore unlikely to run significantly faster in parallel, unless you comment out the output operations. Compilation for parallel execution is achieved by selecting the parallel compiler in the line parallelCXX=... of the Makefile, and by setting the flag MPIparallel=true. The program can the be executed through a call to mpirun, mpiexec, or something similar.

Parallel compilation is activated by default in all Palabos examples: most computers nowadays have multi-core CPUs which you can exploit by running the Palabos applications in parallel.

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:

1. int main(int argc, char* argv[]) {

2. plbInit(&argc, &argv);

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

4.

5. MultiBlockLattice2D<T, DESCRIPTOR> lattice (

6. nx, ny, new BGKdynamics<T,DESCRIPTOR>(omega) );

7.

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

9.

10. defineInitialDensityAtCenter(lattice);

11.

12. // First part: loop over time iterations.

13. for (plint iT=0; iT<maxIter; ++iT) {

14. lattice.collideAndStream();

15. }

16.

17. // Second part: Data analysis.

18. Array<T,2> velocity;

19. lattice.get(nx/2, ny/2).computeVelocity(velocity);

20. pcout << "Velocity in the middle of the lattice: ("

21. << velocity[0] << "," << velocity[1] << ")" << endl;

22.

23. pcout << "Velocity norm along a horizontal line: " << endl;

24. Box2D line(0, 100, ny/2, ny/2);

25. pcout << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

26.

27. plb_ofstream ofile("profile.dat");

28. ofile << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

29.

30. pcout << "Average density in the domain: " << computeAverageDensity(lattice) << endl;

31. pcout << "Average energy in the domain: " << computeAverageEnergy(lattice) << endl;

32. }

Line 19The 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 25With 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 27In 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_ofstreaminstead of ofstream to guarantee working conditions in parallel programs.Line 30-31Many 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:

1. #include "palabos2D.h"

2. #include "palabos2D.hh"

3.

4. #include <vector>

5. #include <iostream>

6. #include <iomanip>

7.

8. /* Code 1.5 in the Palabos tutorial

9. */

10.

11. using namespace plb;

12. using namespace std;

13.

14. typedef double T;

15. #define DESCRIPTOR plb::descriptors::D2Q9Descriptor

16.

17.

18. /// Velocity on the parabolic Poiseuille profile

19. T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const& parameters) {

20. T y = (T)iY / parameters.getResolution();

21. return 4.*parameters.getLatticeU() * (y-y*y);

22. }

23.

24. /// A functional, used to initialize the velocity for the boundary conditions

25. template<typename T>

26. class PoiseuilleVelocity {

27. public:

28. PoiseuilleVelocity(IncomprFlowParam<T> parameters_)

29. : parameters(parameters_)

30. { }

31. /// This version of the operator returns the velocity only,

32. /// to instantiate the boundary condition.

33. void operator()(plint iX, plint iY, Array<T,2>& u) const {

34. u[0] = poiseuilleVelocity(iY, parameters);

35. u[1] = T();

36. }

37. /// This version of the operator returns also a constant value for

38. /// the density, to create the initial condition.

39. void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {

40. u[0] = poiseuilleVelocity(iY, parameters);

41. u[1] = T();

42. rho = (T)1;

43. }

44. private:

45. IncomprFlowParam<T> parameters;

46. };

47.

48. void channelSetup (

49. MultiBlockLattice2D<T,DESCRIPTOR>& lattice,

50. IncomprFlowParam<T> const& parameters,

51. OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )

52. {

53. // Create Velocity boundary conditions.

54. boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

55.

56. // Specify the boundary velocity.

57. setBoundaryVelocity (

58. lattice, lattice.getBoundingBox(),

59. PoiseuilleVelocity<T>(parameters) );

60.

61. // Create the initial condition.

62. initializeAtEquilibrium (

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

64.

65. lattice.initialize();

66. }

67.

68. void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)

69. {

70. const plint imSize = 600;

71. ImageWriter<T> imageWriter("leeloo");

72. imageWriter.writeScaledGif(createFileName("u", iter, 6),

73. *computeVelocityNorm(lattice),

74. imSize, imSize );

75. }

76.

77. int main(int argc, char* argv[]) {

78. plbInit(&argc, &argv);

79.

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

81.

82. // Use the class IncomprFlowParam to convert from

83. // dimensionless variables to lattice units, in the

84. // context of incompressible flows.

85. IncomprFlowParam<T> parameters(

86. (T) 1e-2, // Reference velocity (the maximum velocity

87. // in the Poiseuille profile) in lattice units.

88. (T) 100., // Reynolds number

89. 100, // Resolution of the reference length (channel height).

90. 2., // Channel length in dimensionless variables

91. 1. // Channel height in dimensionless variables

92. );

93. const T imSave = (T)0.1; // Time intervals at which to save GIF

94. // images, in dimensionless time units.

95. const T maxT = (T)3.1; // Total simulation time, in dimensionless

96. // time units.

97.

98. writeLogFile(parameters, "Poiseuille flow");

99.

100. MultiBlockLattice2D<T, DESCRIPTOR> lattice (

101. parameters.getNx(), parameters.getNy(),

102. new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

103.

104. OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*

105. boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

106.

107. channelSetup(lattice, parameters, *boundaryCondition);

108.

109. // Main loop over time iterations.

110. for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {

111. if (iT%parameters.nStep(imSave)==0 && iT>0) {

112. pcout << "Saving Gif at time step " << iT << endl;

113. writeGifs(lattice, iT);

114. }

115. // Execute lattice Boltzmann iteration.

116. lattice.collideAndStream();

117. }

118.

119. delete boundaryCondition;

120. }

Poiseuille profile

Line 20This 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-47This 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 55Specify 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 58Define 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 86An 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 105Choose 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:

1. #include "palabos3D.h"

2. #include "palabos3D.hh"

3. #include <vector>

4. #include <iostream>

5. #include <iomanip>

6.

7.

8. /* Code 1.6 in the Palabos tutorial

9. */

10.

11. using namespace plb;

12. using namespace std;

13.

14. typedef double T;

15. #define DESCRIPTOR plb::descriptors::D3Q19Descriptor

16.

17.

18. /// Velocity on the parabolic Poiseuille profile

19. T poiseuilleVelocity(plint iY, plint iZ, IncomprFlowParam<T> const& parameters) {

20. T y = (T)iY / parameters.getResolution();

21. T z = (T)iZ / parameters.getResolution();

22. return 4.*parameters.getLatticeU() * (y-y*y) * (z-z*z);

23. }

24.

25. /// A functional, used to initialize the velocity for the boundary conditions

26. template<typename T>

27. class PoiseuilleVelocity {

28. public:

29. PoiseuilleVelocity(IncomprFlowParam<T> parameters_)

30. : parameters(parameters_)

31. { }

32. void operator()(plint iX, plint iY, plint iZ, Array<T,3>& u) const {

33. u[0] = poiseuilleVelocity(iY, iZ, parameters);

34. u[1] = T();

35. u[2] = T();

36. }

37. private:

38. IncomprFlowParam<T> parameters;

39. };

40.

41. void channelSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,

42. IncomprFlowParam<T> const& parameters,

43. OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )

44. {

45. // Create Velocity boundary conditions

46. boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

47.

48. setBoundaryVelocity (

49. lattice, lattice.getBoundingBox(),

50. PoiseuilleVelocity<T>(parameters) );

51.

52. lattice.initialize();

53. }

54.

55. void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)

56. {

57. const plint imSize = 600;

58. const plint nx = lattice.getNx();

59. const plint ny = lattice.getNy();

60. const plint nz = lattice.getNz();

61. Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);

62. ImageWriter<T> imageWriter("leeloo");

63. imageWriter.writeScaledGif(createFileName("u", iter, 6),

64. *computeVelocityNorm(lattice, slice),

65. imSize, imSize );

66. }

67.

68. int main(int argc, char* argv[]) {

69. plbInit(&argc, &argv);

70.

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

72.

73. // Use the class IncomprFlowParam to convert from

74. // dimensionless variables to lattice units, in the

75. // context of incompressible flows.

76. IncomprFlowParam<T> parameters(

77. (T) 1e-2, // Reference velocity (the maximum velocity

78. // in the Poiseuille profile) in lattice units.

79. (T) 100., // Reynolds number

80. 30, // Resolution of the reference length (channel height).

81. 3., // Channel length in dimensionless variables

82. 1., // Channel height in dimensionless variables

83. 1. // Channel depth in dimensionless variables

84. );

85. const T imSave = (T)0.02; // Time intervals at which to save GIF

86. // images, in dimensionless time units.

87. const T maxT = (T)2.5; // Total simulation time, in dimensionless

88. // time units.

89.

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

91.

92. MultiBlockLattice3D<T, DESCRIPTOR> lattice (

93. parameters.getNx(), parameters.getNy(), parameters.getNz(),

94. new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

95.

96. OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*

97. //boundaryCondition = createInterpBoundaryCondition3D<T,DESCRIPTOR>();

98. boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

99.

100. channelSetup(lattice, parameters, *boundaryCondition);

101.

102. // Main loop over time iterations.

103. for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {

104. if (iT%parameters.nStep(imSave)==0) {

105. pcout << "Saving Gif at time step " << iT << endl;

106. writeGifs(lattice, iT);

107. }

108. // Execute lattice Boltzmann iteration.

109. lattice.collideAndStream();

110. }

111.

112. delete boundaryCondition;

113. }

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

Line 1 and 3The header files palabos2D.h and palabos2D.hh are replaced by their 3D counterpart.Line 16On 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 19To 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-39The 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 41Initialization 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 55It 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 76It 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 92The 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:

1. // Attribute a value to nx, ny, and nz.

2. Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);

3. plb_ofstream ofile("slice.dat");

4. 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:

1. % Matlab/Octave script

2. % Assign a value to nx and ny

3. load slice.dat

4. slice = reshape(slice, nx, ny);

5. 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:

1. void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,

2. IncomprFlowParam<T> const& parameters, plint iter)

3. {

4. T dx = parameters.getDeltaX();

5. T dt = parameters.getDeltaT();

6. VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);

7. vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);

8. vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt);

9. }

The data in the VTK file is represented in dimensionless variables. For this, the variable dx is given to the constructor of the VtkImageOutput3Dobject 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.