Computational Model

The first step of a simulation consists of setting up a Model. This struct holds the current state of the (discretized) quantities $qᵢ$ as well as all the data that is required to compute their current rate of change $∂qᵢ/∂t$.

Once a Model has been set up, its dynamics can be simulated as discussed in Evolution in Time.

Model(resolution, domain, processes)

A Model provides a discretized representation of the dynamics of the specified processes within the specified domain.

It contains discretized state variables such as the velocity field that are initially set to zero, as well as all the information needed to efficiently compute the rate of change of these state variables. The required state variables are automatically determined from the specified processes.


  • resolution::NTuple{3, Integer}: The number of modes along the $x_1$ and $x_2$ direction, as well as the number of grid points along the $x_3$ direction. The first two values should be odd and will be reduced by one if they are even to ensure that the resolved wavenumbers are symmetric around zero.
  • domain::Domain: The simulated universum, represented with the Domain type.
  • processes: A list of the physical processes that govern the dynamical behavior of the state variables. For convenience, the list can be any iterable collection and can also be nested.


  • comm::MPI.Comm: The MPI communicator that the model will make use of. If not specified, MPI.COMM_WORLD is used. Note that specifying a different communicator is poorly tested currently and may lead to unexpected behavior.


The Domain describes the mathematical domain in which the governing equations are assumed to hold and its relation to the Cartesian domain in which the computations are performed.

Along the first two coordinate directions $x₁$ and $x₂$ (usually streamwise and cross-stream), the domain is assumed to be periodic with a repeating unit of size $L₁×L₂$, i.e.

\[q(x₁, x₂, x₃) = q(x₁+L₁, x₂, x₃) = q(x₁, x₂+L₂, x₃) \,.\]

Along the third coordinate direction $x₃$ (usually wall-normal/vertical), the domain is assumed to span an range of size $L₃$ with boundary conditions specified at the lower and upper limits. The $x₃$-range is represented as a function of $ζ ∈ [0, 1]$ such as $x₃(ζ) = L₃ζ$. By specifying a custom, non-linear function, the computational domain can be deformed such that a discretization that is equidistant in $ζ$ results in grid stretching/refinement in the $x₃$-domain.

The mathematical domain is therefore defined as

\[(x₁, x₂, x₃) ∈ (-∞, +∞) × (-∞, +∞) × \left[x₃(ζ=0), x₃(ζ=1)\right] \,,\]

where boundary conditions have to be specified at $x_3(\zeta=0)$ and $x_3(\zeta=1)$.

Domain([T=Float64,] dimensions, lower_boundary, upper_boundary, [mapping])

Three-dimensional Cartesian domain of size dimensions, periodic along the first two coordinate directions and delimited by lower_boundary and upper_boundary along the third coordinate direction.


  • T::DataType = Float64: Type that is used for coordinates and fields inside domain.
  • dimensions::Tuple: Size of the domain. The third dimension can be specified as a single value, in which case the domain is assumed to start at $x_3=0$ or as a tuple with the minimum and maximum $x_3$ values. If it is omitted, the default of $x_3 ∈ [0,1]$ is assumed, unless a custom mapping is provided.
  • lower_boundary, upper_boundary: Boundary definitions.
  • mapping: A non-linear mapping from the interval $[0,1]$ to the range of $x_3$ values, instead of the default linear mapping. The mapping can either be specified as a tuple of two functions representing the mapping and its derivative, or using the predefined SinusoidalMapping. In the former case, the third element of dimensions is ignored if specified; in the latter case the mapping is adjusted to the domain size.

Domain Boundaries

While the domain is assumed to be periodic along the first two coordinate directions, the behaviors at the minimum and maximum $x₃$-value need to be prescribed. While it is possible to specify the mathematical boundary conditions of each state variable using a CustomBoundary, the boundary definitions SmoothWall, RoughWall, and FreeSlipBoundary provide a convenient shorthand for common types of boundaries and also allow processes to adapt their behavior to the physical properties of those boundaries.

RoughWall(roughness, [von_karman_constant = 0.4])

A wall that is aerodynamically rough with a roughness length scale $z₀$. The mean streamwise velocity near the wall can be assumed to follow a logarithmic profile with the specified von Kármán constant.

CustomBoundary(; boundary_behaviors...)

A boundary that explicitly specifies the mathematical boundary conditions for all state variables. Boundary conditions can be specified as keyword arguments, where the key is the name of the field and the value is either one of :dirichlet or :neumann for homogeneous boundary conditions, or a Pair that includes the non-zero value, such as :dirichlet => 1.


CustomBoundary(vel1 = :dirichlet => 1, vel2 = :dirichlet, vel3 = :neumann)

Domain Mappings

The coordinates along the third dimension can be specified as an arbitrary function $x₃(ζ)$, where $ζ∈[0,1]$. This coordinate transform can be used to reduce the discretization errors of values that are computed on a grid that is uniform in $ζ$ (i.e. grid stretching). To specify a custom mapping, both the function $x₃(ζ)$ and its derivative $\mathrm{d}x₃/\mathrm{d}ζ$ have to be provided so $x₃$-derivatives can be computed with the chain rule. The SinusoidalMapping provides a predefined coordinate transform that automatically adjusts to the specified domain size and boundaries.

SinusoidalMapping(η, variant = :auto)

Define a transformed vertical coordinate with a mapping in the form of

$x_3 = \frac{\sin(ζ η π/2)}{\sin(η π/2)}$

appropriately rescaled such that it maps the range $0≤ζ≤1$ to the $x_3$-range of the Domain.

The parameter $0<η<1$ controls the strength of the grid stretching, where values close to $0$ result in a more equidistant spacing and values close to $1$ result in a higher density of grid points close to the wall(s). The value $η=1$ is not allowed since it produces a vanishing derivative of the mapping function at the wall.

The variant defines at which of the boundaries the coordinates are refined and can be set to :below, :above, :symmetric (both boundaries refined), or :auto (refined for boundaries that produce a boundary-layer).



Within a semiperiodic Domain, a quantity $q$ can be represented as the Fourier series

\[q(x₁, x₂, x₃) = \sum_{\substack{-∞ < κ₁ < ∞ \\ -∞ < κ₂ < ∞}} \hat{q}^{κ₁κ₂}\!\left(x₃(ζ)\right) \,\mathrm{e}^{i 2 π κ₁ x₁ / L₁} \,\mathrm{e}^{i 2 π κ₂ x₂ / L₂} \,.\]

These continuous fields are discretized by resolving $N₁ × N₂$ Fourier modes at $N₃$ grid points.

Along the first two coordinate directions ($i=1,2$), the normalized (integer) wavenumbers $|κᵢ| < Nᵢ / 2$ are resolved. Since this range is symmetric around $κᵢ=0$, there are always an odd number of resolved wavenumbers. When specifying $N₁$ and $N₂$, an even number is therefore equivalent to the next-smaller odd number.

Along the third coordinate direction, the $ζ$-range $[0, 1]$ is split into $N₃$ sections. Numerical solutions are then computed either at the center $ζ_C$ of each section or at the interfaces $ζ_I$ between the sections, with

\[ζ_C ∈ \left\{ \frac{1/2}{N₃}, \frac{3/2}{N₃}, …, \frac{N₃-1/2}{N₃} \right\} \quad \text{and} \quad ζ_I ∈ \left\{ \frac{1}{N₃}, \frac{2}{N₃}, …, \frac{N₃-1}{N₃} \right\} \,.\]

This staggered arrangement can reduce discretization errors and avoid instabilities due to odd–even decoupling. To simplify the notation, we also introduce the definitions

\[Δζ ≡ \frac{1}{N₃} \quad \text{and} \quad ζ¯ ≡ ζ - \frac{Δζ}{2} \quad \text{and} \quad ζ⁺ ≡ ζ + \frac{Δζ}{2}\]

to refer to the grid spacing and the (staggered) neighbors of a $ζ$-location.

A quantity $q(x₁, x₂, x₃)$ is therefore represented by a vector $\bm{\hat{q}}$ of the $N₁×N₂×N₃$ modal/nodal values $\hat{q}^{κ₁κ₂ζ_C}$, or the $N₁×N₂×(N₃-1)$ values $\hat{q}^{κ₁κ₂ζ_I}$ for quantities that are discretized on $ζ_I$ nodes. The elements of $\bm{\hat{q}}$ are complex numbers but since $q$ is real-valued, the values for the wavenumbers $(-κ₁, -κ₂)$ are the complex conjugate of the values for $(κ₁, κ₂)$.

Physical Processes

The Model represents the dynamics of a number of quantities $q_i$ that are each represented by a vector $\bm{\hat{q}}_i$ as described in Discretization. Together they make up the state $\bm{\hat{s}} = \begin{pmatrix} \bm{\hat{q}₁}, \bm{\hat{q}₂}, … \end{pmatrix}^\intercal$. The state usually consists of the three components of the velocity field plus perhaps a number of scalar fields.

The evolution of this state is governed by an equation in the form of

\[\frac{∂ \bm{\hat{s}}}{∂ t} = \sum_i f_i(\bm{\hat{s}}, t) \,,\]

where $f_i$ are functions that describe the contribution of each physical process (molecular diffusion, advection, etc.) to the rate of change of the state. The Model is configured with the desired subset of the Available Processes, which then determine the dynamical behavior of a simulation.

While processes are free to choose how they compute their contribution $f_i$, they generally rely on a shared approach for discrete derivatives and non-linear expressions. Some processes are implemented as a Projections, i.e. they directly modify the state $\bm{\hat{s}}$ instead of contributing to its rate of change. These can be used to enforce invariants of the state, such as enforcing incompressibility with the Pressure process.

Available Processes

Most processes have an argument field::Symbol that specifies the name of the quantity they act on. By convention, the names used for the velocity components are :vel1, :vel2, and :vel3.


In principle, each process is free to define its own way of handling spatial derivatives, as long as they can be applied to the discretized form of the state. However, the available processes all rely on exact frequency-domain derivatives along the $x₁$- and $x₂$-directions and second-order finite differences along the $x₃$-direction.

To simplify the notation for discretized processes, the definitions

\[∂₁(κ₁) ≡ \frac{i 2 π κ₁}{L₁} \quad \text{and} \quad ∂₂(κ₂) ≡ \frac{i 2 π κ₂}{L₂} \quad \text{and} \quad ∂₃(ζ) ≡ \left( Δζ \frac{\mathrm{d}x₃}{\mathrm{d}ζ} \bigg|_{ζ} \right)^{-1}\]

are introduced for the “derivation factors” that appear in expressions of discretized spatial derivatives.

This gives the following expressions for the discrete computation of derivatives:

\[\widehat{\frac{∂^n q}{∂x₁^n}}^{κ₁κ₂ζ} = ∂₁(κ₁)^n \, \hat{q}^{κ₁κ₂ζ} \quad \text{and} \quad \widehat{\frac{∂^n q}{∂x₂^n}}^{κ₁κ₂ζ} = ∂₂(κ₂)^n \, \hat{q}^{κ₁κ₂ζ}\]

\[\widehat{\frac{∂q}{∂x₃}}^{κ₁κ₂ζ} = ∂₃(ζ) \left( \hat{q}^{κ₁κ₂ζ⁺} − \hat{q}^{κ₁κ₂ζ¯} \right) + \mathcal{O}(Δζ²)\]

\[\widehat{\frac{∂²q}{∂x₃²}}^{κ₁κ₂ζ} = ∂₃(ζ) \left( ∂₃(ζ⁺) \left( \hat{q}^{κ₁κ₂(ζ+Δζ)} − \hat{q}^{κ₁κ₂ζ} \right) − ∂₃(ζ¯) \left( \hat{q}^{κ₁κ₂ζ} − \hat{q}^{κ₁κ₂(ζ−Δζ)} \right) \right) + \mathcal{O}(Δζ²)\]

Near the boundary, one-sided finite differences with the same order of accuracy are used where necessary. To compute second derivatives at the first $ζ_C$-location, $q$ is extrapolated to $ζ=\frac{−Δζ}{2}$ with

\[\hat{q}^{κ₁κ₂(−Δζ/2)} = \frac{16}{5}\, \hat{q}^{κ₁κ₂0} − 3\, \hat{q}^{κ₁κ₂(Δζ/2)} + 1\, \hat{q}^{κ₁κ₂(3Δζ/2)} − \frac{1}{5}\, \hat{q}^{κ₁κ₂(5Δζ/2)}\]

where $\hat{q}^{κ₁κ₂0}$ is a Dirichlet boundary condition. The upper boundary at $ζ=1$ is treated in an analogous manner.

Non-linear Processes

Performing non-linear operations on frequency-domain representations of a quantity can be challenging. Non-linear processes therefore rely on a physical-domain representation computed at $N₁^\mathrm{PD} × N₂^\mathrm{PD}$ points. These values can be obtained with an inverse FFT of the appropriate size. After computing non-linear expressions at these points, the result is transformed back to a frequency-domain representation with a forward FFT. Note that unless the physical-domain resolution is fine enough to represent all wavenumbers produced by the non-linear operations, such a computation will produce aliasing errors in addition to the truncation errors that are inherent to the frequency-domain discretization.


Projection processes are defined as a linear function $f_i(\bm{\hat{s}}, t) = F_i \, λ(\bm{\hat{s}}, t)$ of some quantity $λ$ that enforces an affine constraint of the form $C_i \bm{\hat{s}} = \bm{c}_i$ on the state $\bm{\hat{s}}$. The main role of projections is to accommodate the Pressure process for incompressible flows, but they can also be used for other functionality such as the ConstantMean source term.

In keeping with the discussion of time-integration methods, the explicit form of $f_i$ can be written as

\[f_i(\bm{\hat{s}}, t) = − F_i (C_i F_i)^{−1} C_i \sum_{j ≠ i} f_j(\bm{\hat{s}}, t) \,,\]

but the implementation relies on a projection function

\[P_i(\bm{\hat{s}}) = \bm{\hat{s}} − F_i (C_i F_i)^{−1} \left(C_i \bm{\hat{s}} − \bm{c}_i\right)\]



There is no special functionality to simultaneously enforce multiple constraints. If more than one projections are enabled, the $P_i(\bm{\hat{s}})$ are applied consecutively. This is valid as long as all projections are independent with $C_i F_j = 0$ for all combinations $(P_i, P_j)$ of enabled projection.