MeanProfiles(; kwargs...)

Collect profiles of terms, averaged over time and horizontal space. The time average is computed with the trapezoidal rule using samples collected before/after each time step. The averages are saved to HDF5 files, which can be written regularly during a simulation.

Keywords

• profiles = (:vel1, :vel2, :vel3): Terms of which profiles are collected. Currently supported: :vel1, :vel2, :vel3, :adv1, :adv2, :adv3, :sgs11, :sgs12, :sgs13, :sgs22, :sgs23, :sgs33.
• path = "output/profiles": Prefix of the path at which HDF5-files with the profiles are saved. The files are numbered, so the first file would be written to output/profiles-01.h5 for the default path.
• output_frequency: The interval with which the output is written to disk. Note that this should be a multiple of the time step, otherwise some outputs may be skipped.

ProgressMonitor(; kwargs...)

Print information about an ongoing simulation to the standard output at regular intervals of wall time.

Keywords

• frequency = 10: The wall-time interval (in seconds) between each progress report. Reports are always displayed at the end of a time step, if the wall time since the previous report exceeds the specified interval.
• always_update = true: Print a dot after each time step between regular reports to indicate the ongoing progress of the simulation.
• tstep = nothing: Passing the time step $Δt$ to the ProgressMonitor allows computing an advective Courant number $C = Δt / (∂₃^{max} |u|^{max})$. If not set, the advective time scale $1 / (∂₃^{max} |u|^{max})$ is reported instead.

Snapshots(; kwargs...)

Save the instantaneous state $\mathbf{\hat{s}}$ to a collection of Cartesian Binary Data files at regular intervals during the simulation.

Keywords

• frequency: The frequency (in units of simulation time) at which snapshots should be saved. Note that this should be a multiple of the time step $Δt$.
• path = "output/snapshots": Path to a folder inside which the snapshots will be saved.
• centered = true: Save values at the $x₁$- and $x₂$-locations at the center of the intervals obtained by dividing the domain into $N₁×N₂$ segments. If set to false, the locations start at 0 instead and end one grid spacing before the end of the domain.
• precision::DataType: The floating-point data type used in the output, either Float32 or Float64. By default this is set to the data type used for the simulation by the Model.

StepTimer(; kwargs...)

Collect time stamps every $N$ steps to measure the compute time per time step and write it to a JSON file.

Keywords

• path = "output/timestamps.json": The path to the output file.
• frequency = 1: The number of time steps $N$ computed between each collected time stamp.

AB2()

Second-order Adams–Bashforth method for time integration. This corresponds to an explicit linear multi-step methods with parameters $α₀=1$, $β₀=3/2$, and $β₁=−1/2$.

Euler()

First-order explicit (“forward”) Euler method for time integration. This corresponds to an explicit linear multi-step methods with parameters $α₀=1$ and $β₀=1$.

The ODEBuffer structs hold the variables each method requires to compute the next time step without allocating new memory.

Each algorithm should implement a method init_buffer(prob, alg) that takes the ODEMethod as an argument and creates the corresponding buffer, which will then be passed every time perform_step! is called.

The ODEMethod structs are used as arguments to the time integration methods to specify the algorithm used for integration. They are only used to select the correct methods in multiple-dispatch functions.

SSPRK22()

Two-stage second-order strong-stability-preserving Runge–Kutta method for time integration, with the Shu–Osher coefficients $α₁₀=β₁₀=1$ and $α₂₀=α₂₁=β₂₁=1/2$ (see Gottlieb et al. 2009).

SSPRK33()

Three-stage third-order strong-stability-preserving Runge–Kutta method for time integration, with the Shu–Osher coefficients $α₁₀=β₁₀=1$, $α₂₀=3/4$, $α₂₁=β₂₁=1/4$, $α₃₀=1/3$, and $α₃₂=β₃₂=2/3$ (see Gottlieb et al. 2009).

perform_step!(problem, dt, buffer)

Advance the solution problem.u by a time step of dt. The type of buffer determines which integration scheme is used and stores intermediate values required by the scheme.

The method should assume that problem.dt is up-to-date and is allowed to overwrite it. It is not necessary to update problem.dt at the end of the step, as this is done by the caller.

solve!(problem, algorithm, dt, checkpoints = nothing)

Step the time integration problem prob forward in time with the specified algorithm and time step, hitting all the checkpoints exactly during the integration and signaling them to the rate! function of the problem.

Transform a field from the frequency domain to an extended set of nodes in the physical domain by adding extra frequencies set to zero.

ConstantMean(field, mean_value = 1)

Source term for a scalar quantity $q$ with a source strength that is constant in space but dynamically adjusted in time to maintain a constant mean value $Q$ for $q$.

Arguments

• field::Symbol: The name of the quantity $q$.
• mean_value::Real: The mean value $Q$ that is maintained.

ConstantSource(field, strength = 1)

Source term for a scalar quantity $q$ with source strength $S$ that is constant in space and time.

Arguments

• field::Symbol: The name of the quantity $q$.
• strength::Real: The source strength $S$.

The DistributedBatchLDLt{P,Tm,Tv,B} holds the data of a decomposition A=L·D·Lᵀ for a batch of symmetric tridiagonal matrices A. The purpose is to solve a system A·x=b efficiently for several A, x, and b. The vectors and matrices are distributed across MPI processes, with P designating the role of the current MPI process. Tm is the element type of the matrices A while Tv is the element type of the vectors x and b.

Since the algorithm for solving the linear system is sequential by nature, there is a tradeoff between solving many systems at once, maximizing the throughput of each individual solution, and solving smaller batches of systems, minimizing the time MPI processes have to wait for each other. This tradeoff is represented by the value of B, which is the number of systems solved in batch at a time. Ideally the optimal value should be found automatically by measuring the computation time for different values of B and choosing the optimal value.

The type is constructed from an iterable type (such as a generator) holding the diagonal vectors dvs and another one holding the off-diagonals evs. The diagonals are distributed amongst MPI processes like variables on H-nodes while the off-diagonals are distributed like variables on V-nodes. An example batch of b-vectors (or x-vectors) is also passed to the constructor, since those might have a different element type than the matrices.

Internally, we store a set of vectors γ and another set of vectors β, representing the L and D matrices of the L·D·Lᵀ decomposition. D is the diagonal matrix with γ as its diagonal vector while L is the matrix with γ¯¹ as the diagonal and β as the lower off-diagonal. Another vector is allocated as buffer for the values exchanged between MPI processes when solving the system A·x=b.

The algorithm used for the decomposition and the solution of the linear systems can be found in section 3.7.2 of the book “Numerical Mathematics” by Quarteroni, Sacco & Saleri (2007).

MolecularDiffusion(field, diffusivity)

Diffusive transport of a scalar quantity $q$ with a diffusion coefficient $D$ that is constant in space and time.

Arguments

• field::Symbol: The name of the quantity $q$.
• diffusivity::Real: The diffusion coefficient $D$.

MomentumAdvection(; dealiasing = :quadratic)

Advective transport of momentum ($u_i$, i.e. normalized by density) in rotation form. Only the velocity–vorticity products are computed, while the gradient of kinetic energy is assumed to be handled as part of a modified pressure term.

Keywords

• dealiasing: The level of dealiasing, determining the physical-domain resolution at which multiplications are performed. The default :quadratic setting relies on the “3/2-rule” to avoid aliasing errors for the products that are computed as part of the momentum advection term. Alternatively, dealiasing can be set to nothing for a physical-domain resolution that matches the frequency-domain resolution (rounded up to an even value), or to a tuple of two integers to set the resolution manually.

Pressure(batch_size = 64)

Transport of momentum by a pressure-like variable that enforces a divergence-free velocity field.

Arguments

• batch_size::Int: The number of wavenumber pairs that are included in each batch of the tri-diagonal solver. The batching serves to stagger the computation such that different MPI ranks can work on different batches at the same time.

StaticSmagorinskyModel(; kwargs...)

Subgrid-scale advective transport of resolved momentum, approximated with the static Smagorinsky model.

Keywords

• Cs = 0.1: The Smagorinsky coefficient $C_\mathrm{S}$.
• dealiasing: The level of dealiasing, determining the physical-domain resolution at nonlinear operations are performed. The default nothing gives a physical-domain resolution that matches the frequency-domain resolution (rounded up to an even value). Alternatively, dealiasing can be set to :quadratic for a resolution based on the “3/2-rule” or to a tuple of two integers to set the resolution manually. Aliasing errors should be reduced for larger resolutions but they can be expected at any resolution.
• wall_damping = true: Adjust the Smagorinsky lenghtscale towards the expected value of $κ x₃$ towards the wall.
• wall_damping_exponent = 2: Exponent of the near-wall adjustment.

Compute the gradient of a scalar field and add it to a vector field. The scalar field and the horizontal components of the vector field are defined on C-nodes while the vertical component of the vector field are defined on I-nodes. An optional prefactor can be used to rescale the gradient before it is added.

Compute the Laplacian of a scalar field and add it to a different scalar field. Both fields have to be defined on the same set of nodes. An optional prefactor can be used to rescale the Laplacian before it is added.

For k1=k2=0, we want to replace the last equation of the linear system with the equation 3 p(N₃-1/2) - p(N₃-3/2) = 0, which (approximately) sets the mean pressure at the top of the domain to zero. To retain the symmetry of the matrix, this function computes the last value of the off-diagonal and and sets the last value of the diagonal to minus three times that value. For this to work correctly, the last entry of the vector for which the system is solved should also be set to zero. Otherwise, the solver will still work correctly but the pressure will differ by some additive factor.

The pressure solver computes a pressure field p such that the provided velocity field becomes divergence-free when the pressure gradient is subtracted, i.e. ∇·(u − ∇p) = 0. The solver only requires boundary conditions for the vertical velocity component due to the way the staggered grid is set up.

Note that the formulation of this solver does not correspond exactly to the way pressure appears in the Navier-Stokes equations and the problem has to be rescaled to apply the solver. It is important that the provided boundary conditions are compatible with the rescaled problem.

In practice, the solver is employed in two different ways. One option is to apply it to a velocity field to obtain its projection into a divergence-free space. This corresponds more or less to the formulation given above, but the pressure does not really have a physical meaning in this case. The other option is to apply it to one or all terms of the RHS to obtain the pressure field they induce. In this case, the u of the solver is really a ∂u/∂t and the boundary conditions provided to the pressure solver should be those of the time derivative of the velocity, i.e. constant values becomes zero.

The second option is useful for obtaining pressure fields in post-processing, but for time stepping it is best to rely on the first formulation and apply the the solver to u + Δt RHS rather than directly to RHS. Otherwise, errors can accumulate in u and the solution may drift away from a divergence-free flow field over 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.

Arguments

• 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.

Keywords

• 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.

evolve!(model, tspan; dt, method, output)

Simulate a model’s evolution over time, optionally collecting data along the way.

Arguments

• model::Model: The Model containing the current state as well as the discretized processes that describe the rate of change of the state.
• tspan: The time span over which the evolution of the model should be simulated. Can be a single Number with the total duration or a Tuple of Numbers with the start and end times.

Keywords

• dt: The (constant) size of the time step (required). Note that tspan has to be divisible by dt.
• method = SSPRK33(): The time-integration method used to solve the semi-discretized initial-value problem as an ordinary differential equation.
• output = []: A list of output modules that collect data during the simulation.

boundary_layer(layers, ind, bc, extra_layers=())

Pass the ind-th layer away from the boundary to the process responsible for the boundary, avoiding communication if the data is already on the right process. Other processes should not use the return value.

NOTE: This only has been verified to work for C-nodes and probably needs adjustments for I-nodes.

Arguments

• layers: The local layers of each process.
• ind: The index of the desired layer. Positive values are counted from the lower boundary, negative values from the upper boundary.
• bc: The BoundaryCondition for the boundary that the data should be passed to. If communication is required, the buffer of this boundary condition will be used.
• extra_layers: Additional layers that are already available at the boundary process, e.g. from earlier communication. It is assumed that every process has the same number of extra layers and that they are ordered away from the boundary.

This function takes a field defined on C-nodes, converts it into layers, and expands them through communication such that the list includes the layers just above and below all I-nodes. This means passing data down throughout the domain. The boundary condition is only used for passing data, the actual values are unused.

This function takes a field defined on I-nodes, converts it into layers, and expands them through communication such that the list includes the layers just above and below all C-nodes. This means passing data up throughout the domain.

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.

Example

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

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.

Arguments

• 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.

FreeSlipBoundary()

A no-penetration boundary with vanishing gradients for tangential velocity components.

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.

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).

SmoothWall()

A wall that is aerodynamically smooth and has no-slip, no-penetration boundary conditions for the velocity field.

Computes the constant coordinate scaling along dimension dim, i.e. 1/L where L is the domain size. This can only be computed for coordinates that are at most linearly transformed and will return an error otherwise.