`BoundaryLayerDynamics.Logging.MeanProfiles`

— Type`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.

`BoundaryLayerDynamics.Logging.ProgressMonitor`

— Type`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.

`BoundaryLayerDynamics.Logging.Snapshots`

— Type`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`

.

`BoundaryLayerDynamics.Logging.StepTimer`

— Type`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.

`BoundaryLayerDynamics.ODEMethods.AB2`

— Type`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$.

`BoundaryLayerDynamics.ODEMethods.Euler`

— Type`Euler()`

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

`BoundaryLayerDynamics.ODEMethods.ODEBuffer`

— TypeThe `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.

`BoundaryLayerDynamics.ODEMethods.ODEMethod`

— TypeThe `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.

`BoundaryLayerDynamics.ODEMethods.SSPRK22`

— Type`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).

`BoundaryLayerDynamics.ODEMethods.SSPRK33`

— Type`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).

`BoundaryLayerDynamics.ODEMethods.perform_step!`

— Function`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.

`BoundaryLayerDynamics.ODEMethods.solve!`

— Method`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.

`BoundaryLayerDynamics.PhysicalSpace.get_field!`

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

`BoundaryLayerDynamics.Processes.ConstantMean`

— Type`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.

`BoundaryLayerDynamics.Processes.ConstantSource`

— Type`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$.

`BoundaryLayerDynamics.Processes.DistributedBatchLDLt`

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

`BoundaryLayerDynamics.Processes.MolecularDiffusion`

— Type`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$.

`BoundaryLayerDynamics.Processes.MomentumAdvection`

— Type`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.

`BoundaryLayerDynamics.Processes.Pressure`

— Type`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.

`BoundaryLayerDynamics.Processes.StaticSmagorinskyModel`

— Type`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.

`BoundaryLayerDynamics.Processes.add_gradient!`

— FunctionCompute 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.

`BoundaryLayerDynamics.Processes.add_laplacian!`

— FunctionCompute 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.

`BoundaryLayerDynamics.Processes.adjust_d3g3_diagonal`

— MethodFor 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.

`BoundaryLayerDynamics.Processes.solve_pressure!`

— MethodThe 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.

`BoundaryLayerDynamics.Model`

— Type`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.

`BoundaryLayerDynamics.evolve!`

— Method`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`Number`

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

`BoundaryLayerDynamics.BoundaryConditions.boundary_layer`

— Function`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.

`BoundaryLayerDynamics.BoundaryConditions.layers_c2i`

— MethodThis 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.

`BoundaryLayerDynamics.BoundaryConditions.layers_i2c`

— MethodThis 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.

`BoundaryLayerDynamics.Domains.CustomBoundary`

— Type`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)`

`BoundaryLayerDynamics.Domains.Domain`

— Type`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.

`BoundaryLayerDynamics.Domains.FreeSlipBoundary`

— Type`FreeSlipBoundary()`

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

`BoundaryLayerDynamics.Domains.RoughWall`

— Type`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.

`BoundaryLayerDynamics.Domains.SinusoidalMapping`

— Type`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).

`BoundaryLayerDynamics.Domains.SmoothWall`

— Type`SmoothWall()`

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

`BoundaryLayerDynamics.Domains.scalefactor`

— MethodComputes 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.