Stencils

Defines a D-dimensional stencil info with Q directions.

The stencil module consists of two parts:

• Definition of all directions and their properties: Directions.h
• Automatically generated stencil classes. To generate new stencils, look at generate.py, to adapt how the stencil classes look like, change Stencil.in.h (which is used as input by the generate script). This generation process cannot be done using C++ template mechanisms, so here the python script is used as a simple form of code generation.

Directions

In Directions.h, all possible directions are defined as an enumeration. Each direction has several properties. For example: the inverse direction, the offset in x,y,z direction, etc. These properties are stored as constant arrays which can be indexed using the direction enum value:

cx[W]          // = 1
dirToString[W] // = "W"
inverseDir[W]  // = E
dirToString[inverseDir[W]] // = "E"

The only context in which direction enums should be used as integers, are above property arrays. Do not use them as indices to the fourth coordinate of fields!

An alternative implementation would have been to store directions as classes which hold these properties (cx, inverseDir, ...) as members. This was not done since this prevents some compile-time optimizations of the compiler. For best possible performance, these compile-time optimizations are, however, essential! Above code snippets can all be evaluated at compile time, since all arrays are declared constant.

Stencil Classes

Generation Mechanism

With all possible directions and their properties defined, a stencil class now only has to pick a subset of these directions. The generation of the stencil classes are automated by the python script generate.py. The script uses the file Stencil.in.h as template and substitutes strings according to the stencil definition. Which stencils are generated is specified at the beginning of the script:

directions = ['C','N','S','W','E','T','B',
              'NW','NE','SW','SE','TN','TS',
              'TW','TE','BN','BS','BW','BE',
              'TNE','TNW','TSE','TSW','BNE','BNW','BSE','BSW'];
 
# edit this to add new stencils
stencils = [
            { 'name' : 'D2Q4', 'dim' : 2, 'dirs' : directions[1: 5] },
            { 'name' : 'D2Q5', 'dim' : 2, 'dirs' : directions[ : 5] },
            { 'name' : 'D3Q7', 'dim' : 3, 'dirs' : directions[ : 7] },
            { 'name' : 'D3Q19','dim' : 3, 'dirs' : directions[ :19] },
            { 'name' : 'D3Q27','dim' : 3, 'dirs' : directions[ :27] }
            #also possible: pick directions as you need them in arbitrary order:
            #{ 'name' : 'CustomStencil','dim' : 3, 'dirs' : ['S','NW','N'] }
           ]

Stencil Class

A stencil is defined by picking a subset of all possible directions. These directions are listed in the dir array. So when iterating a stencil, one iterates the dir array which contains the directions. To get the properties for these directions, the global arrays of Directions.h can be used:

using namespace stencil;
for( uint_t i = 0; i < $name::Size; ++i )
{
   int cx = cx[ $name::dir[i] ];
   Direction inverse = inverseDir[ $name::dir[i] ];
}

Since all these arrays are constant, the lookup can be resolved at compile time and no performance overhead should be generated.

For more convenient iteration, use the provided iterator. The following code is equivalent to the code example above:

using namespace stencil;
for( auto dir = $name::begin(); dir != $name::end(); ++dir )
{
  int cx = dir.cx();
  Direction inverse = dir.inverseDir();
}

UseCases

Stencils and Fields

In LBM simulation, one typically stores a PDF for every direction contained in a stencil. To write code that is independent of the stencil, the stencil type is usually passed as template parameter. The size of the f coordinate (= fourth dimension of the field) equals the number of directions contained in the stencil and can be obtained by Stencil::Size (or simply Stencil::Q).

One can iterate all directions of the stencil using the Stencil::dir[] array or the Stencil::iterator obtained via Stencil::begin(), Stencil::beginNoCenter(), and Stencil::end().

The reverse mapping is also possible: When one has a stencil::Direction enum at hand and needs to determine the index of this direction, the Stencil::idx[] array has to be used. Do not use the enum directly as index: field( x, y, z, W ) is not correct in the general case! Typically you want field( x, y, z, Stencil::idx[W] )!

There are special cases (for example D3Q19 or D3Q27) where both variants are coincidentally the same.

field( x, y, z, W )               // most probably WRONG !
field( x, y, z, Stencil::idx[W] ) // typically what you want and therefore CORRECT !