Preprocessor

The preprocessor preproc is responsible for splitting the computational domain 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. 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);
• other types of BC 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.

preproc can elaborate geometries described either analytically or through STL files. While the latter supports every type of (implemented) BC, the former currently only supports periodic and wall-type BC.

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 STL geometries, preproc automatically retrieves the boundaries of the BB by scanning all the triangles contained in the STL files. For analytically-described geometries, instead, it is up to the user to define the computational domain and its BB.

Mesh generation

XXX ABBIAMO UNA FIGURA DEL CUBO? XXX Figure. CAPTION TEXT

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 or unsetting the corresponding macro GRID_NON_UNIFORM in the file headers.cpl. The type of grid step variation is selected, independently for each direction, in the section GRID_NON_UNIFORM of the settings.in file, by assigning the variables xGridType, yGridType and zGridType. The available mesh types can be found in modules-preproc/grid-refinement.cpl. Currently, only uniform (type u) and hyperbolic tangent (type h) grids are implemented. The framework is easily extended: one simply adds a new grid type and writes in grid-refinement.cpl the line of code for the desired variation of the grid step.

Non-uniform grid

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 points where the Navier--Stokes equations might be solved; whether or not in a specific point this actually takes place (because the point is internal) depends on the geometrical complexity of the domain. For a simple channel flow, for example, all points inside the bounding box are internal (fluid) 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 the program). Green arrows are the u velocity collocation point assigned by the user. Red arrows are the ghost nodes (automatically defined by the program).

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 grid-refinement.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 grid points generated in grid-refinement.cpl must be inside the BB for the grid to be valid, i.e.:

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

This requirement is not checked: when creating a new type of grid refinement, it is mandatory to verify it.

Also, the red-black iterative algorithm used to enforce incompressibility requires 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 the preprocessing log file.

Please note that, when using non-uniform grid generation in periodic directions with STL geometry, additional requirements exist. They are addressed later XXX LINK XXX.

Analytic-geometry mode

In this mode, only periodic and wall-type (no-slip and no-penetration) BC are presently supported. To have preproc work in this mode, the macro #define GEOMETRY_ANALYTIC must be defined (and the macro GEOMETRY_STL not defined) in headers.cpl before compilation. Note that at least one macro is required, and that the two macros are mutually exclusive. An error message is issued if these conditions are not met.

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

• in the section GEOMETRY_ANALYTIC of the settings.in file, the user sets the dimensions of the BB, i.e. Lx, Ly, Lz, and assigns periodicity to any combination of directions (e.g. periodic=1 0 1 means periodic BC will be used in x and z);
• in the file geometry-analytic.cpl, the boolean function inBody(x,y,z) contains the analytical description of the geometry, and returns YES when the point of coordinates (x,y,z) belongs to the solid.

STL-geometry mode

In this mode, all the available BC are supported. To have preproc work in this mode, the macro GEOMETRY_STL must be defined (and the macro GEOMETRY_ANALYTIC not defined) in headers.cpl before compilation.

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 identify 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 a no-slip BC 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;
• the STL files associated with each periodic direction 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 in 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 Figure \ref{fig:x-grid-example} 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 a unique file.

A set of STL files satisfying these constraints is processed correctly by both preproc and solver. Note that not all constraints are or can be checked at runtime.

Detection of the bounding box

The bounding box is the smallest parallelepiped that encloses all the triangles of the STL files. preproc reads all the STL files, identifies minima and maxima of the coordinates across all triangles, and stores them in the variables 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 raytracing step afterwards. Each non-periodic planar BB face, where the grid starts from, is thus shifted by $8 \times\epsilon$ outwards towards the solid, where $\epsilon$ is the machine epsilon of float variables (since the STL triangles are by definition saved in single precision).

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 collocation points; the corresponding velocity collocation points are also flagged. The lists and their length are eventually saved by preproc in each cube-specific file. This information will enable solver to apply the BC only on the labelled 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 its exterior, treated as solid. As described below, when computing the IB coefficients a ray-tracing algorithm is used to discriminate between internal (fluid) and external (solid) regions. The ray-tracing 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-tracing 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 calcImbc. 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-tracing 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 performs the ray-tracing algorithm to identify the solid points which are further away from the triangles.

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

• it initializes star to internal;
• it calls scanSTL;
• 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-tracing

XXX MISSING TEXT XXX

Pruning the empty cubes

Before saving cubes data to disk, empty cubes are pruned, i.e. 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.