# 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);
1. decomposition of the BB into cubes;
1. calculation of the IB coefficients and of the geometrical features of the boundary conditions;
1. 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.

<!--An improved domain decomposition strategy as the *building-cube method* \cite{nakahasi-kim-2004} may be implemented, further reducing the operations directly implemented in `preproc.cpl`.-->

## 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

<!-- 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/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.
These`Nx` 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. 

<a name="x-grid-example"></a>

![x-grid-example](assets/x-grid-example.png)
**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](preproc.md#stl-geometry-mode).


## 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.
1. 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`.
1.  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. 
1. 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.
