Preprocessor

The preprocessor preproc is responsible for splitting the computational domain (described either analytically or through STL files) into cubes, and for elaborating the geometrical data needed by solver. These tasks are dealt with in four separate parts of the source, that communicate with each other through a small set of variables. The modular structure of preproc enables it to deal with multiple configurations. Parallel execution is supported on CPUs only: since the computational demand of preproc is limited, parallelism is only required to access large quantities of memory.

In the case directory, preproc writes preproc-log.dat, which contains useful information on the various stages of pre-processing, some warnings when the number of points is modified to make it even, and important information that is needed later, as the values of nDecompCubes and nTotCubes. Additional useful information is contained in the shared file general-info.dat.

The parts of preproc concern:

1. identification of the bounding box (BB);
2. decomposition of the BB into cubes;
3. calculation of the IB coefficients and of the geometrical features of the boundary conditions;
4. pruning of the empty cubes, and final write to disk.

Different types of BCs are dealt with differently:

• no-slip BC (type w) are imposed via the IB coefficients;
• periodic BC (type p) are hard-coded in the communications between cubes (this is why periodic BC must be applied on whole faces of the BB);
• the other types are dealt with in two steps; preproc only manages their geometrical description, that is read later by solver, which links it to the corresponding BC type and executes the corresponding routine.

The main program

The file preproc.cpl contains all the general-purpose operations, including:

• initialization and identification of the bounding box (BB);
• mesh generation;
• decomposition of the BB into cubes;
• call to possible module-specific routines required for customization;
• calculation of the IB coefficients and of the geometrical features of the boundary conditions;
• pruning empty cubes and writing data to memory.

Bounding box

The grid discretizes the region of space contained in the bounding box. The BB is a cuboid defined by its lower and upper boundary in the three directions; the coordinates of the BB boundaries are stored in the arrays bbMin and bbMax.

For analytically-described geometries, instead, it is up to the user to define the computational domain and its BB.

For STL geometries, preproc automatically retrieves the boundaries of the BB by scanning all the triangles contained in the STL files, to identify minima and maxima of the coordinates across all triangles, stored in bbMin and bbMax. In a second step, and only for the non-periodic directions, the BB is slightly shifted to avoid limit cases where a triangle happens to have one vertex exactly on the BB: this would cause finite-precision issues in the ray-casting step afterwards. Each non-periodic planar BB face, where the grid starts from, is thus shifted by $8 \times \epsilon$ outwards towards the solid; $\epsilon$ is the machine epsilon of float variables (since the STL triangles are by definition saved in single precision).

Mesh generation

The IBM adopts a staggered grid, with pressure defined at the center of a cubic cell and velocities defined at its faces. The code supports both uniform and non-uniform rectilinear grids, with arbitrary (one-dimensional) variations of the grid step size in every direction. The choice between uniform and non-uniform grid is made at compile time, by setting/unsetting the corresponding macro GRID_IS_NONUNIFORM in the file headers.cpl.

The process of mesh generation in the $x$ direction is described here with some detail (see the figure below for a sketch); the same goes for the other directions. The program generates Nx collocation points for the u velocity component; they fall within the interval [bbMin(0), bbMax(0)] defined by the extrema of the BB in the $x$ direction. TheseNx points are where the discrete Navier--Stokes equations might be solved; whether or not in a specific point this actually takes place (because the point is internal, i.e. in the fluid) depends on the geometrical complexity of the domain. For a simple channel flow, for example, all points inside the bounding box are internal points, since the bounding box and the computational domain coincide.

Figure Example of uniform grid in the $x$ direction with Nx=10. Blue dots represent pressure collocation points (automatically defined by preproc). Green arrows are the u velocity collocation points assigned by the user. Red arrows are the ghost nodes (automatically defined by preproc).

Ghost nodes are required to compute $x$ derivatives in the two extreme points at i=1 and i=Nx. The ghost nodes, at i=0 and i=Nx+1, are positioned outside the BB: in this way, pressure collocation points, automatically placed by preproc half-way between two velocity collocation points, end up exactly on the side of the BB. Therefore, given bbMin(0), bbMax(0) and the internal velocity collocation points x(1), x(2), ... x(Nx) computed in grids.cpl according to the specific grid distribution strategy, preproc computes the ghost points such that:

${ \left\lbrace \begin{aligned} &x^u_0 = x_{min} - 0.5 (x^u_1 - x_{min}) \\ &x^u_{Nx+1} = x_{max} + 0.5 (x_{max} - x^u_{Nx}) \\ &x^p_i = \frac{x^u_{i-1} + x^u_i}{2} \qquad \text{for $i=1:N_x+1$} \\ &x^p_0 = x^p_1 - (x^u_1 - x^u_0) \end{aligned} \right.}$

The obvious consequence is that the first and last points must be inside the BB for the grid to be valid, i.e.:

$x^u_1 > x_{min}, \quad x^u_{Nx} < x_{max} .$

This requirement is not checked, and must be fulfilled when creating a new grid type.

Also, the iterative algorithm used to enforce incompressibility implies a further restriction when periodic boundary conditions are used. Along periodic directions, the number of points given as input by the user must be even. Therefore, if an odd number is given, preproc automatically adds one point to make it even, and issues a warning in preproc-log.dat.

Non-uniform grids

The type of grid step variation is selected, independently for each direction, in the section NONUNIFORM_GRID of the settings.in file, by assigning the variables xGridType, yGridType and zGridType. The available mesh types can be found in modules-preproc/grids.cpl. Uniform (type u), hyperbolic tangent (type h) and cosine-spaced (type c) grids are implemented. The framework is such that new grid types can be easily added: one simply writes in grids.cpl the line of code for the desired variation of the grid step.

When using a non-uniform grid along periodic directions with STL geometry, additional requirements exist. They are addressed later.

Analytic-geometry mode

To have preproc work in this mode, the macro #define GEOMETRY_IS_ANALYTIC must be defined (which is mutually exclusive with the macro GEOMETRY_IS_STL) in headers.cpl before compilation. An error message is issued if none or both macros are defined.

All the analytic-geometry actions are defined through functions contained in modules-preproc/geometry-analytic.cpl. For generality, functions are defined as <function name><Geometry> in a geometry-specific file placed in the geometry-analytic directory; they are then linked to the global function it refers to via the line <function name>==<function name><Geometry>. This ensures that only the geometry-specific function has to be specified, while the modules-preproc/geometry-analytic module remains unchanged.

Setting up preproc in analytic-geometry mode requires four steps:

1. Fill in the section ANALYTIC_GEOMETRY of settings.in, by providing the dimensions of the bounding box (Lx, Ly, Lz), the number of computational boundaries (nFaces), and a list of nFaces types of boundary conditions, eventually completed by specific parameters.
2. Define the number of planes in which BC other than w and p are applied. This is done by the subroutine computenOtherBCFiles<Geometry>. An example is computenOtherBCFilesChannel in modules-common/specifics-channel.cpl.
3. Define the boolean function inBody<geometry>(x,y,z), which contains the analytical description of the geometry, and returns YES only if the point of coordinates (x,y,z) belongs to the solid. An example can be found in the modules-common/specifics-channel.cpl module, where inBodyChannel is defined.
4. Define the buildBoundingSquare<Geometry> subroutine. This is actually enabled twice for each non-periodic and non-wall direction only (once for the lower and once for the upper boundary). For each boundary, it fills the structure BSquareDelimiter by specifying the edges of the plane containing the boundary condition prescribed in the three Cartesian directions, and the direction of the wall-normal vector. An example for the plane channel flow is available in modules-common/specifics-channel.cpl: in this case the boundaries of the computational domain coincide with those where the boundary conditions have to be applied; however, other flows, like the flow around a solid body, can be studied with this infrastructure.

Once these functions are defined, preproc looks for cubes intersected by the boundary conditions chosen, flags every relevant point, and stores this information in the cube-specific file.

STL-geometry mode

To have preproc work in this mode, the macro #define GEOMETRY_IS_STL must be defined (which is mutually exclusive with the macro GEOMETRY_IS_ANALYTIC) in headers.cpl before compilation. An error message is issued if none or both macros are defined.

Setting up preproc in STL-geometry mode requires to provide a list of STL files describing the input geometry. The surfaces delimiting the computational domain can be wall surfaces, periodic planes and surfaces where other BC are applied. The names of the STL files are listed in, and read from the input file settings.in. All the STL-specific actions are defined through functions contained in modules-preproc/geometry-stl.cpl.

The set of STL files must satisfy the following constraints:

• the computational domain (i.e. the fluid region) must be fully enclosed by the surface defined by the set of STL files;
• the normals of the triangles of the STL files must point towards the fluid domain;
• apart from those STL files where the no-slip BC w is used, the STL files must describe a planar surface, orthogonal to one of the coordinate axes. preproc checks planarity of these surfaces, and returns an error if the check fails (within a tolerance that depends on the grid spacing);
• periodic BC are applied to whole planar faces of the BB; fluid regions extending beyond a periodic plane will simply not be considered as fluid regions, and the behavior of preproc is in this case undefined;
• STL files for which periodic conditions are selected must be a pair: one STL for the lower boundary, and one for the upper boundary;
• the $y$ direction cannot be periodic, unless the calculation is tri-periodic, as for homogeneous isotropic turbulence;
• for an indefinite wall delimited by periodic planes, the IB coefficients near the boundary are computed assuming that the wall is symmetric with respect to the periodic plane. When using grid refinement in periodic directions, this condition is particularly strict as it requires the grid to be symmetric near the boundary of the domain. Referring to the grid figure above, this implies:

${\left \lbrace \begin{aligned} &x^u_1 - x_{min} = x_{max} - x^u_{Nx} \\ &x^u_2 - x_{min} = x_{max} - x^u_{Nx-1} \end{aligned} \right.}$

• the surface described by one STL file can take only one BC; to be described by the same STL file, two distinct surfaces must have the same orientation (and the same normal) with respect to the fluid domain, the same type of BC, and the same boundary value imposed at the boundary;
• multiple wall files are allowed; they are treated as one.

A set of STL files satisfying these constraints is processed correctly by preproc and solver. Not all constraints are (and some cannot be) checked at run time.

Identification of the BCs

The boundary conditions are elaborated by preproc to extract the geometrical information needed by solver.

preproc checks first if periodic BC are present, and makes sure that all the non-wall STL files meet the requirement of being planar and orthogonal to one coordinate axis. Then, for each STL file, preproc finds the cubes intersected by a certain BC plane. For each such cube, a list is compiled containing all the cube grid points belonging to the BC plane: a boolean is attached to each point of the BC plane, set to YES when the point is inside the triangles of the STL file of the BC. The check is carried out at the pressure points; the corresponding velocity points are also flagged. The lists and their length are eventually saved by preproc in each cube-specific file. This information enables solver to apply the BC only at the flagged points, which are actually inside the STL file.

Besides identifying points where BC will be applied, the STL files also define the interface between the computational domain and the external (solid) volume. As described below, a ray-casting algorithm is used to discriminate between internal (fluid) and external (solid) regions. The ray-casting algorithm is run along the $y$ direction, and spans the whole domain. To correctly enumerate the fluid/solid interfaces, STL lists which do not intersect the cube but are aligned with the cube in the ray-casting direction are needed too. These are only used when computing the IB coefficients; therefore, they are not saved together with the ones which intersect the cube, as they do not represent a BC to be applied inside the cube.

Calculation of the IB coefficients

The IB coefficients are computed by the subroutine calcImbCoeffs. Due to the large memory requirements, this calculation is carried out independently for each $y-z$ slab of each cube. The memory requirement ensues from the ray-casting algorithm (currently implemented in the $y$ direction), which must allocate enough points to cover the whole domain in that direction. The IB coefficients are computed by summing all the 6 contributions coming from the stencil of the 3D Laplacian operator. These are stored, for each slab, in a data structure called star: the subroutine scanSTL computes the value of star for each point.

The routine scanSTL considers a $y-z$ slab, and performs the following operations:

• it updates the value of star in the points where the intersection with a triangle of the wall STL file is found;
• it executes the ray-casting algorithm to discriminate solid and fluid regions, based on the number of solid-fluid interfaces, including the BC surfaces.

At the higher level, for each $y-z$ slab, calcImbCoeffs performs the following actions:

• it initializes star to internal;
• it calls scanSTL, which fills star;
• it computes the IB coefficients for the slab, only in the region of interest for the cube;
• it sets to SOLID the IB coefficients of the points where BCs other than periodic and wall are imposed;
• it writes the IB coefficients to disk.

Ray-casting in scanSTL and usage of star

The information on which point is internal or external, as well as on the distance of each grid point along the three Cartesian directions from the STL surface is stored in the data structure star, later used to compute the IBM coefficients in the subroutine calcImbCoeffs.

star is an array of size $N_y \times N_z \times 3 \times 2$. The first two indices span the $y-z$ plane, the third index signals which direction is being considered, and the last index stores the sign of the distance from the nearest wall.

To compute the relative position of each point from the nearest wall, the function intersect is used for each STL triangle, each grid point, and the three directions. The triangle is projected over the plane perpendicular to each Cartesian axis, the baricentric coordinates are computed to establish whether the grid point falls within the projected triangle; when it does, the (signed) distance from the grid point to the triangle along the considered direction is computed. If this distance is lower than the distance between grid points, and of the previous distance stored in star, it is then saved in star.

Once star is preliminary filled in this way, a ray-casting algorithm is used along the $y$ direction, taking into account the number of interfaces between solid and fluid from solid walls. This information is then stored in star, which is set to zero for the external points. There is also a second step, where boundary conditions other than w cutting the $y-z$ plane are considered; star is modified accordingly, by removing domain regions that are excluded by their presence.

Pruning the empty cubes

Before saving cubes data to disk, empty cubes are pruned, i.e. the cubes whose number of fluid points is zero are removed. Once computed, the cube-specific number of fluid points is broadcast and saved in a data structure within the cube file. This allows in principle to implement a static load balancing for solver, where the number of internal cube points can be used as weight without loading all the IB coefficients to memory. Once each MPI rank has received the number of internal points of the cubes residing in other ranks, the connectivity of the cubes is re-established. The connectivity of an empty cube ID is set to -1 (cube IDs start from 1). Only at this point, data files for non-pruned cubes are written to disk.