Reference

struct Boundary{S<:Symbol, UR<:UnitRange{<:Integer}}
    name::S         # Boundary patch name
    IDs_range::UR   # range to access boundary info (faces and boundary_cellsID)
end

struct Cell{F<:AbstractFloat, SV3<:SVector{3,F},UR<:UnitRange{<:Integer}}
    centre::SV3     # coordinate of cell centroid
    volume::F       # cell volume
    nodes_range::UR # range to access cell nodes in Mesh3.cell_nodes
    faces_range::UR # range to access cell faces info (faces, neighbours cells, etc.)
end

struct Face3D{
    F<:AbstractFloat, 
    SV2<:SVector{2,<:Integer},
    SV3<:SVector{3,F}, 
    UR<:UnitRange{<:Integer}
    }
    
    nodes_range::UR # range to access face nodes in Mesh3.face_nodes
    ownerCells::SV2 # IDs of face owner cells (always 2)
    centre::SV3     # coordinates of face centre
    normal::SV3     # face normal unit vector
    e::SV3          # unit vector in the direction between owner cells
    area::F         # face area
    delta::F        # distance between owner cells centres (boundary: cell centre to nearest face point)
    weight::F       # linear interpolation weight
end

struct Mesh3{VV, VTF, VR, VO, VI, VS, VB} <: AbstractMesh

3D unstructured mesh. Mesh2 has identical fields. Each cell, face and node property is stored as its own array; mesh.cells, mesh.faces and mesh.nodes return views that index as Cell, Face3D and Node elements. Boundary faces come first in the face arrays.

    cell_centre::VV      # cell centroid coordinates
    cell_volume::VTF     # cell volumes
    cell_nodes_range::VR # range of each cell's nodes in cell_nodes
    cell_faces_range::VR # range of each cell's faces in cell_faces, cell_neighbours, cell_nsign
    cell_nodes::VI       # node IDs of each cell
    cell_faces::VI       # internal face IDs of each cell
    cell_neighbours::VI  # neighbour cell IDs across each cell face
    cell_nsign::VS       # face normal sign per cell face (1 or -1)
    face_nodes_range::VR # range of each face's nodes in face_nodes
    face_ownerCells::VO  # owner cell ID pair of each face (equal on boundary faces)
    face_centre::VV      # face centre coordinates
    face_normal::VV      # face unit normals
    face_e::VV           # unit vectors between owner cell centres
    face_area::VTF       # face areas
    face_delta::VTF      # distance between owner cell centres (boundary: cell centre to nearest face point)
    face_weight::VTF     # linear interpolation weights
    face_nodes::VI       # node IDs of each face
    face_gDiff::VTF      # Laplacian face coefficient (derived, see `_gDiff`)
    boundaries::VB       # boundary patches
    node_coords::VV      # node coordinates
    node_cells_range::VR # range of each node's cells in node_cells
    node_cells::VI       # cell IDs around each node
    get_float::VTF       # one-element array tagging the mesh float type
    get_int::VI          # one-element array tagging the mesh integer type
    boundary_cellsID::VI # owner cell ID of each boundary face
    boundary_cells::VI   # owner cells of boundary faces, each once (derived)
    boundary_cell_faces_range::VR # range of each boundary cell's faces in boundary_cell_faces
    boundary_cell_faces::VI # boundary face IDs grouped by owner cell, ascending within a cell

Kernels that add boundary-face contributions to cells loop over boundary_cells and, for each, over its faces in boundary_cell_faces, so every cell sums its faces in one fixed order: results do not depend on the thread count or the scheduling.

struct Node{SV3<:SVector{3,<:AbstractFloat}, UR<:UnitRange{<:Integer}}
    coords::SV3     # node coordinates
    cells_range::UR # range to access neighbour cells in Mesh3.node_cells
end

pmin, pmax = bounding_box(mesh)

Return the corners of the axis-aligned box enclosing all mesh nodes as SVector{3}s of the mesh float type. On a distributed mesh the box encloses the nodes of all ranks.

reorder_mesh!(mesh::Union{Mesh2,Mesh3}; method=:rcm) -> mesh

Reorder the cells, faces and nodes of a host mesh in place so that neighbouring cells sit close in memory. Cell loops read the values of each cell's neighbours, so meshes whose generator numbers cells without regard to adjacency (typical of unstructured tetrahedral meshes) run faster once reordered. Geometry, face orientation and boundary patches are unchanged; only the numbering differs. The extra memory is a few index vectors of the mesh's integer type and one buffer the size of the largest connectivity list.

• method=:rcm: reverse Cuthill-McKee ordering of the cell graph (bounded index distance between neighbours).
• method=:morton: Z-order curve of the cell centres (compact groups of consecutive cells).

Internal faces follow their lowest-numbered cell, boundary faces their owner cell within each patch, and nodes the order in which the cells first reach them. The mesh is left unchanged unless the new order shortens the mean index distance between neighbouring cells by at least 10%.

Only the mesh in memory changes; mesh files are not touched. Fields built on the reordered mesh, and results written from them, are in the new order: VTK output carries the mesh in every file, and OpenFOAM output replaces a constant/polyMesh numbered differently before writing any field (in a distributed run, distribute and partition_mesh do this when they partition the mesh). Serial meshes only; call before adapt, and before distribute or partition_mesh to reorder a mesh for a distributed run.

Example

mesh = reorder_mesh!(UNV3D_mesh("mesh.unv", scale=0.001))

UNV2D_mesh(meshFile; scale=1, integer_type=Int32, float_type=Float64)

Read and convert 2D UNV mesh file into XCALibre.jl

Input

• meshFile – path to mesh file.

Optional arguments

• scale – used to scale mesh file e.g. scale=0.001 will convert mesh from mm to metres defaults to 1 i.e. no scaling

• integer_type - integer type of the mesh indices; Int64 is needed only when a mesh has more than 2^31 faces, face-node entries or matrix entries, and reading such a mesh as Int32 stops with an error saying so

• float_type - select interger type to use in the mesh (Float32 may be useful on GPU runs)

UNV3D_mesh(unv_mesh::String; scale::Real=1.0, integer_type::Type=Int32, float_type::Type=Float64) -> Mesh3

Constructs a Mesh3 object from a Universal (UNV) file format.

Arguments

• unv_mesh::String: The file path to the UNV mesh file.

Keyword Arguments

• scale::Real=1.0: A scaling factor applied to the nodal coordinates.
• integer_type::Type=Int32: The integer type used for topological indices and connectivity arrays; Int64 is needed only when a mesh has more than 2^31 faces, face-node entries or matrix entries, and reading such a mesh as Int32 stops with an error saying so.
• float_type::Type=Float64: The floating-point type used for coordinates and geometric properties.

Returns

• Mesh3: A fully constructed 3D mesh object containing nodes, cells, faces, boundaries, and populated geometric properties.

FOAM3D_mesh(mesh_file; scale=1, integer_type=Int32, float_type=Float64)

Read and convert 3D OpenFOAM mesh file into XCALibre.jl. Note that, at present, it is not recommended to run 2D cases using meshes imported using this function.

Input

• mesh_file – path to mesh file.

Optional arguments

• scale – used to scale mesh file e.g. scale=0.001 will convert mesh from mm to metres defaults to 1 i.e. no scaling

• integer_type - integer type of the mesh indices; Int64 is needed only when a mesh has more than 2^31 faces, face-node entries or matrix entries, and reading such a mesh as Int32 stops with an error saying so

• float_type - select interger type to use in the mesh (Float32 may be useful on GPU runs)

FirstTouch()

Adapt.jl target that copies every Array of a structure with the backend given to activate_multithread; see first_touch.

activate_multithread(backend::CPU; nthreads=1)

Set the CPU backend that fields, equations and preconditioners are built with, and the number of OpenBLAS threads.

Input arguments

• backend is the only required input which must be CPU() from KernelAbstractions.jl. With CPU(static=true) and pinned threads, arrays are allocated by parallel first touch, so each thread's chunk sits in its NUMA domain; pass the mesh, built before, through first_touch.
• nthreads is the number of BLAS threads (default: 1)

first_touch(mesh)

Copy mesh (or any structure of Arrays) in parallel so that, with CPU(static=true) passed to activate_multithread and pinned threads, each thread first writes, and so places in its NUMA domain, the chunk of cells, faces or nodes it later works on. The mesh is read on one thread, so call it before building the model: mesh = first_touch(mesh). Pays on multi-socket (NUMA) nodes; results are unchanged.

struct ScalarField{VF,M,BC} <: AbstractScalarField
    values::VF  # scalar values at cell centre
    mesh::M     # reference to mesh
    BCs::BC     # store user-provided boundary conditions
end

struct VectorField{S1<:ScalarField,S2,S3,M,BC} <: AbstractVectorField
    x::S1   # x-component is itself a `ScalarField`
    y::S2   # y-component is itself a `ScalarField`
    z::S3   # z-component is itself a `ScalarField`
    mesh::M
    BCs::BC
end

function initialise!(field, value) # dummy function for documentation
    # Assign `value` to field in-place
    nothing
end

This function will set the given field to the value provided in-place. Useful for initialising fields prior to running a simulation.

Input arguments

• field specifies the field to be initialised. The field must be either a AbractScalarField or AbstractVectorField
• value defines the value to be set. This should be a scalar or vector (3 components) depending on the field to be modified e.g. for an AbstractVectorField we can specify as value=[10,0,0]

Note: in most cases the fields to be modified are stored within a physics model i.e. a Physics object. Thus, the argument value must fully qualify the model. For example, if we have created a Physics model named mymodel to set the velocity field, U, we would set the argument field to mymodel.momentum.U. See the example below.

Example

initialise!(mymodel.momentum.U, [2.5, 0, 0])
initialise!(mymodel.momentum.p, 1.25)

initialise! also accepts a functions arguments and initialises a field by evaluating the function at each cell centre. The function should have the signature f(x, y, z) and it must return either a scalar or an SVector type (requirement to work on GPUs)

CrankNicolson(ψ=1)

Crank-Nicolson time scheme in trapezoidal form, with every other term of the equation at tⁿ⁺¹: the Time term is the time derivative at tⁿ⁺¹, ddtⁿ⁺¹ = (1 + ψ)(xⁿ⁺¹ - xⁿ)/Δt - ψ ddtⁿ, where ddtⁿ is kept from the previous step. ψ = 1 is the second-order trapezoidal rule, ψ = 0 is backward Euler, and values in between damp the oscillations the trapezoidal rule leaves undamped. CrankNicolson without arguments is CrankNicolson(1). The first time step of a run is a backward-Euler step.

Dirichlet <: AbstractDirichlet

Dirichlet boundary condition model.

Inputs

• ID Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• value Scalar or Vector value for Dirichlet boundary condition

DirichletFunction(ID, value) <: AbstractDirichlet

Dirichlet boundary condition defined with user-provided function.

Input

• ID Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• value Custom function or struct (<:XCALibreUserFunctor) for Dirichlet boundary condition
• IDs_range Range of indices to access boundary patch faces

Function requirements

The function passed to this boundary condition must have the following signature:

f(coords, time, index) = SVector{3}(ux, uy, uz)

For a vector-value field, or:

f(coords, time, index) = p::Number

for a scalar field. In both cases, coords is the coordinate vector of the face, time is the current time (or iteration counter in steady simulations), and index is the local face index from 1 to N, where N is the number of faces on the boundary. The function must return either a three-component SVector (StaticArrays.jl) or a plain numeric value, depending on whether the boundary condition represents a vector or scalar field.

Empty <: AbstractPhysicalConstraint

Empty boundary condition model (currently only configured for zero gradient)

Fields

• 'ID' – Boundary ID
• value – Scalar or Vector value for Empty boundary condition.

Extrapolated <: AbstractNeumann

This boundary condition extrapolates the face value using the interior cell value. Equivalent to setting a zero gradient boundary condition (semi-implicitly). It applies to both scalar and vector fields.

Example

Extrapolated(:outlet)

FixedHeatFlux <: AbstractNeumann

Fixed wall heat flux boundary condition for a temperature field.

Inputs

• ID Name of the boundary given as a symbol (e.g. :tankWall). Internally it gets replaced with the boundary index ID.
• value Heat flux in W/m^2, positive into the domain.

Example

FixedHeatFlux(:tankWall, 3.5)

Implementation

For a temperature equation containing - Laplacian{scheme}(keff, T), the diffusive term integrated over a boundary cell contributes

-sum_f keff * area * (dT/dn)_out

and a prescribed inward flux q fixes keff * (dT/dn)_out = q on that face (the conductive flux vector is -keff*grad(T), so heat flowing in means the outward normal gradient is positive). The face contribution is therefore the known constant -q*area on the left-hand side, i.e. +q*area on the right.

The coefficient is independent of keff: prescribing the flux prescribes the whole term. That also means q = 0 reduces exactly to Zerogradient, which returns (0, 0) — a useful cross-check on the sign.

FixedTemperature(name::Symbol, model::Enthalpy; T::Number)

Fixed temperature boundary condition model, which allows the user to specify wall temperature that can be translated to the energy specific model, such as sensible enthalpy.

Inputs

• name: Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• T: keyword argument use to define the boundary temperature
• model: defines the underlying Energy model to be used (currently only Enthalpy is available)

Example

FixedTemperature(:inlet, T=300.0, Enthalpy(cp=cp, Tref=288.15))

IEnergy{C,F}

BC converter for internal energy formulation: e = cv*(T - Tref) where cv is the specific heat at constant volume.

Fields

• cv: Specific heat at constant volume.
• Tref: Reference temperature.

Example

FixedTemperature(:inlet, T=300.0, IEnergy(cv=718.0, Tref=288.15))

LUSTDC <: AbstractScheme

LUST in deferred-correction form. The matrix holds the upwind coefficients, so it stays diagonally dominant, and the face correction $\frac{3}{4}(\phi_{lin} - \phi_U)$ is added explicitly from the current field, with $\phi_{lin}$ the linearly interpolated and $\phi_U$ the upwind face value. A converged steady solution has the face values of LUST, $\phi_f = \frac{3}{4}\phi_{lin} + \frac{1}{4}\phi_U$. Boundary faces use the LUST boundary coefficients (periodic patches keep its implicit coupling). No gradient is needed. Intended for steady runs: in transient runs the correction lags one time step, so results differ from LUST by $O(\Delta t)$ and the stable time step is smaller; use LUST for transient and LES runs. This is not OpenFOAM's LUST, which keeps the linear part implicit and corrects the remaining quarter explicitly with the upwind-cell gradient, as LinearUpwind does.

LUSTLU <: AbstractScheme

LUST plus an explicit second-order correction of its upwind share, for transient and LES runs. The matrix holds the LUST coefficients (75 % linear, 25 % upwind), and the upwind quarter of the face value is corrected explicitly with the gradient of the upwind cell, as in LinearUpwind:

$\phi_f = \frac{3}{4}\phi_{lin} + \frac{1}{4}\left[\phi_U + (\nabla\phi)_U \cdot (\mathbf{x}_f - \mathbf{x}_U)\right]$.

The face value is second order on any mesh and LUST's first-order dissipation is removed. The gradient is supplied through the equation's operator, Divergence{LUSTLU}(flux, phi, grad); the solvers pass the gradient of U (limited with Schemes(limiter=...)). Boundary faces take the LUST coefficients; periodic patches also get the correction. Without a gradient the scheme reduces to LUST (as in the film and multiphase solvers). It is not available for he.

LeastSquares <: AbstractScheme

Gradient scheme, selected with Schemes(gradient=LeastSquares). The gradient of a cell is the weighted least-squares fit $\min \sum_f w_f (\phi_f - \phi_P - \nabla\phi\cdot\mathbf{d}_f)^2$ over its faces, with $\mathbf{d}_f$ the vector from the cell centre to the neighbour centre (internal faces) or to the face centre (boundary faces, which use the boundary value) and $w_f = 1/|\mathbf{d}_f|^2$, so near neighbours dominate on stretched meshes. It is exact for linear fields on any mesh, skewed or stretched, where Gauss is not, and is computed cell by cell from a 3 x 3 normal matrix inverted once per Grad.

LimitedLinearUpwind <: AbstractScheme

Bounded form of LinearUpwind. The face correction $\delta$ is limited per face with the downwind value, $0 \le \delta/\Delta \le \min(r, 1)$, where $\Delta = \phi_D - \phi_U$ and $r = 2 (\nabla\phi)_U \cdot (\mathbf{x}_D - \mathbf{x}_U)/\Delta - 1$. The face value then lies between the upwind and downwind values and local extrema are convected with upwind values. Steady solutions stay bounded on regular meshes, and over- and undershoots are much reduced on irregular ones. The limiter is inactive for linear fields on uniform meshes ($r = 1$) and reduces the scheme to first order at extrema and discontinuities.

LinearUpwind <: AbstractScheme

Second-order upwind divergence scheme. The matrix holds the upwind coefficients and the face value is corrected explicitly (deferred correction) with the gradient of the upwind cell, $\phi_f = \phi_U + (\nabla\phi)_U \cdot (\mathbf{x}_f - \mathbf{x}_U)$. The gradient is supplied through the equation's operator, Divergence{LinearUpwind}(flux, phi, grad); the solvers pass the gradient of U (limited with Schemes(limiter=...)). Boundary faces use upwind values. Without a gradient the scheme reduces to Upwind.

Neumann <: AbstractNeumann

Neumann boundary condition model to set the face normal gradient at the boundary explicitly

Inputs

• ID Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• value Scalar providing face normal gradient

Example

Neumann(:outlet, 0)

NeumannFunction(ID, value) <: AbstractNeumann

Neumann boundary condition defined with user-provided function.

Input

• ID Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• value Custom function defining the desired Neumann boundary condition.

Function requirements

The function passed to this boundary condition has not yet been implemented. However, users can pass a custom struct to specialise the internal implementations of many functions. By default, at present, this function will assign a zero gradient boundary condition.

Outlet <: AbstractNeumann

Outlet boundary condition with backflow prevention.

For the density-based solver, this behaves as a zero-gradient (extrapolated) outlet when flow is leaving the domain. When backflow is detected (interior normal velocity directed into the domain), an impermeable wall flux is applied to prevent spurious mass ingestion. This is equivalent to a one-way valve at the boundary.

For pressure-based solvers (SIMPLE/PISO), this behaves identically to Zerogradient and applies to both scalar and vector fields.

Example

Outlet(:outlet)

struct Periodic{I,V,R<:UnitRange} <: AbstractPhysicalConstraint
    ID::I
    value::V
end

Implicit implementation of Periodic boundary condition. Note that to apply this condition two periodic patch pairs need to be constructed using the function construct_periodic. The implementation currently requires conformant patch pairs.

Fields

• ID is the name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• value tuple containing information needed to apply this boundary

Robin <: AbstractBoundary

Robin (mixed) boundary condition a·φ + b·∇φ·n = value for scalar fields. a=1, b=0 recovers Dirichlet and a=0, b=1 recovers Neumann. Not supported by the density-based (Godunov) solver (run! throws an ArgumentError).

Inputs

• ID Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• a coefficient of the boundary value (default 1)
• b coefficient of the face normal gradient (default 0)
• value right-hand side of the constraint (default 0)

Example

Robin(:wall, a=1.0, b=0.5, value=10.0)

RotatingWall <: AbstractPhysicalConstraint

RotatingWall boundary condition model for no-slip for rotating RotatingWalls (rotating around an axis). It has been strictly implemented to work with vectors only.

Inputs

• rpm provides the rotating speed (internally converted to radian/second)
• centre vector indicating the location of the rotation centre e.g. [0,0,0]
• axis unit vector indicating the rotation axis e.g. [1,0,0]

Example

RotatingWall(:inner_wall, rpm=100, centre=[0,0,0], axis=[0,0,1])

Slip <: AbstractPhysicalConstraint

Slip boundary condition for vector and scalar fields. Vectors keep both tangential components and have the face-normal component removed, so the patch is impermeable but exerts no tangential shear. Scalars use an explicit zero-normal-gradient condition. Use Wall instead when the patch should apply no-slip.

Input

• ID is the boundary name (for example, :plate). It is replaced by the boundary index during boundary assignment.

Example

Slip(:plate)

Symmetry <: AbstractPhysicalConstraint

Symmetry boundary condition for vector and scalar fields. Scalars use an explicit zero-normal-gradient condition. Vectors remove the face-normal component while retaining both tangential components.

Input

• ID is the boundary name (for example, :freestream). It is replaced by the boundary index during boundary assignment.

Example

Symmetry(:freestream)

TimeRate(mesh)
rate!(dst, rate, term, q, q0, time)

Rate dq/dt of an explicitly evaluated quantity (such as dp/dt in an energy equation), with the scheme of the equation's Time term term: (q - qⁿ)/Δt, or under Crank-Nicolson (1 + c)(q - qⁿ)/Δt - c (dq/dt)ⁿ, the rate of the step before being kept by rate.

Wall <: AbstractDirichlet

Wall boundary condition model for no-slip or moving walls (linear motion). It should be applied to the velocity vector, and in most cases, its scalar variant should be applied to scalars.

Inputs

• ID represents the name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID
• value should be given as a vector for the velocity e.g. [10,0,0]. For scalar fields such as the pressure the value entry can be omitted or set to zero explicitly.

Examples

Wall(:plate, [0, 0, 0]) # no-slip wall condition for velocity
Wall(:plate) # corresponding definition for scalars, e.g. pressure
Wall(:plate, 0) # alternative definition for scalars, e.g. pressure

Zerogradient <: AbstractNeumann

Zerogradient boundary condition model (explicitly applied to the boundary)

Input

• ID Name of the boundary given as a symbol (e.g. :inlet). Internally it gets replaced with the boundary index ID

Example

Zerogradient(:inlet)

construct_periodic(mesh, backend, patch1::Symbol, patch2::Symbol; tol=1e-8)

Function for construction of periodic boundary conditions.

Input

• mesh – Mesh.
• backend – Backend configuraton.
• patch1 – Primary periodic patch ID.
• patch2 – Neighbour periodic patch ID.
• tol – keyword argument providing the tolerance used to find matching face pairs

Output

• periodic::Tuple - tuple containing boundary definitions for patch1 and patch2 i.e. (periodic1, periodic2). The fields of periodic1 and periodic2 are

  • ID – Index to access boundary information in mesh object
  • value represents a PeriodicValue struct with the following fields:
    • patchID – boundary/patch ID
    • transform – stores information to apply the patch pair matching e.g. LinearTransform.
    • face_map – vector providing indices to faces of match patch
    • ismaster – flat to identify one of the patch pairs as the main patch

Example

- `periodic = construct_periodic(mesh, CPU(), :top, :bottom)` - Example using CPU 
backend with periodic boundaries named `top` and `bottom`.

time_level!(term, prev, rho_prev, time, config)

Returns the field the Time term reads as its old level when the equation is discretised at time. A Crank-Nicolson term advances its stored levels on the first call of a new time and returns its source field; other schemes return prev.

macro define_boundary(boundary, operator, definition)
    quote
        @inline (bc::$boundary)(
            term::Operator{F,P,I,$operator}, cellID, zcellID, cells, faces, fID, i, component, time
            ) where {F,P,I} = $definition
    end |> esc
end

Macro to reduce boilerplate code when defining boundary conditions (implemented as functors) and provides access to key fields needed in the implementation of boundary conditions, such as the boundary cell and face objects (more details below)

Input arguments

• boundary specifies the boundary type being defined
• operator specifies the operator to which the condition applies e.g. Laplacian
• definition provides the implementation details

Available fields

• term reference to operator on which the boundary applies (gives access to the field and mesh)
• cellID ID of the corresponding boundary cell
• zcellID sparse matrix linear index for the cell
• cells the mesh cells; read a boundary cell property as cells.volume[cellID]
• faces the mesh faces; read a boundary face property as faces.area[fID]
• fID ID of the boundary face
• cell, face whole boundary cell and face elements, built only when the body names them (slower to compile than column reads)
• i local index of the boundary faces within a kernel or loop
• component for vectors this specifies the components being evaluated (access as component.value). For scalars component = nothing
• time provides the current simulation time. This only applies to time dependent boundary implementation defined as functions or neural networks.

Example

Below the use of this macro is illustrated for the implementation of a Dirichlet boundary condition acting on the Laplacian using the Linear scheme:

@define_boundary Dirichlet Laplacian{Linear} begin
    J = term.flux[fID]      # extract operator flux
    area, delta = faces.area[fID], faces.delta[fID] # boundary face information
    flux = J*area/delta     # calculate the face flux
    ap = term.sign*(-flux)  # diagonal (cell) matrix coefficient
    ap, ap*bc.value         # return `ap` and `an`
end

When called, this functor will return two values ap and an, where ap is the cell contribution for approximating the boundary face value, and an is the explicit part of the face value approximation i.e. ap contributes to the diagonal of the sparse matrix (left-hand side) and an is the explicit contribution assigned to the solution vector b on the right-hand of the linear system of equations $Ax = b$

AMG(; coarsening=SmoothAggregation(), smoother=AMGJacobi(), coarse_solve=OnDevice(), kwargs...) <: PreconditionerType

XCALibre's native algebraic multigrid, used as SolverSetup(solver=Cg(), preconditioner=AMG(...)) (symmetric systems such as pressure), with solver=Bicgstab() (non-symmetric systems) or with solver=Richardson() (repeated V-cycles, no Krylov acceleration). On a distributed mesh PETSc's GAMG() is used instead, with a one-time warning. Mixed precision comes from coarse_storage=Float32, not from precision. See the AMG section of the user guide for the keywords.

AMGGaussSeidel(; sweep=AMGSymmetricSweep(), iterations=1)

Gauss-Seidel AMG smoother (CPU only). sweep selects the traversal order: AMGSymmetricSweep() (forward + backward), AMGForwardSweep(), or AMGBackwardSweep(). iterations sets the number of sweeps per smoother call. omega is fixed at 1 (no relaxation); use AMGSOR for a tunable direct relaxation factor.

AMGJacobi(; omega=4/3)

Weighted-Jacobi AMG smoother. omega is a lambda_max-scaled coefficient: the effective relaxation factor applied each sweep is ω_eff = omega / lambda_max, where lambda_max is the spectral radius of D⁻¹A at that level. The default omega=4/3 gives ω_eff ≈ 2/3 for Poisson-like operators (where lambda_max ≈ 2), which damps the upper half of the spectrum optimally. Valid scaled range: (0, 2) — values ≥ 2 are clamped for SPD stability. Note: unlike AMGSOR, omega here is NOT a direct relaxation factor.

AMGSOR(; omega=1.0, sweep=AMGSymmetricSweep(), iterations=1)

Successive Over-Relaxation AMG smoother (CPU only). omega is the direct relaxation factor applied each sweep: x_new = (1-omega)*x_old + omega*(D\(b - L*x)). Stable range (0, 2); omega=1 recovers Gauss-Seidel. Unlike AMGJacobi, omega is NOT scaled by lambda_max — it is used as-is.

AdaptiveTimeStepping(; 
    # keyword arguments

    maxCo=0.75,
    maxAlphaCo=0.5,
    minShrink=0.1,
    maxGrow=1.2
)

Constructs an AdaptiveTimeStepping object used to control automatic time-step adjustment based on the Courant number.

This struct is passed optionally to Runtime and enables adaptive time stepping in transient simulations. If not provided, a fixed time step is used.

Input arguments

• maxCo::AbstractFloat: target maximum Courant number. The time step will be adjusted such that the computed Courant number approaches this value.
• maxAlphaCo::AbstractFloat: target maximum Alpha Courant number. The time step will be adjusted such that the computed Courant number approaches this value.
• minShrink::AbstractFloat: lower bound on the multiplicative factor applied to the current time step. Prevents excessively large reductions in a single update.
• maxGrow::AbstractFloat: upper bound on the multiplicative factor applied to the current time step. Prevents excessive time-step growth.

BoomerAMG(; kwargs...) <: PreconditionerType

HYPRE BoomerAMG via PETSc (-pc_type hypre). On a serial mesh XCALibre's native AMG runs instead, with a one-time warning. Needs a PETSc with hypre, which the stock Float64 libraries include (Float32 ones do not). SPD systems only (e.g. pressure). Use it where GAMG converges poorly or the strongest reduction per solve matters; GAMG is the default recommendation, since it needs no special build and is the better AMG on the GPU. Like any AMG, results vary slightly with the rank count, and it needs more memory than Jacobi.

Defaults are tuned for 3D (BOOMERAMG_3D_DEFAULTS): PETSc's own BoomerAMG defaults build a high-complexity hierarchy in 3D that degrades as the mesh grows. On a GPU with a CUDA-built hypre, PETSc switches to PMIS coarsening and l1-Jacobi smoothing, so GPU and CPU residuals differ.

freeze=N holds the whole preconditioner fixed for N solves and rebuilds it from the current matrix on the Nth (default 10); in between, Krylov iterates on the current matrix with a hierarchy built from an older one. In SIMPLE the pressure matrix changes VALUES only, so a stale hierarchy stays a good preconditioner, and a full hypre rebuild every solve otherwise dominates runtime. freeze=1 rebuilds every solve. Longer freezes cut setup cost further but let the stale hierarchy degrade the pressure residual.

On GPU fields the solver refuses BoomerAMG() unless PETSc's hypre reports device execution, since a host-only hypre crashes there; device=true skips that check for a build the query misreads.

Each other keyword k=v overrides a default and is passed to PETSc as -pc_hypre_boomeramg_<k> v. Common knobs: strong_threshold (0.5-0.7 in 3D), coarsen_type ("HMIS"/"PMIS"/"Falgout"), interp_type ("ext+i"/"classical"), P_max (interpolation entries per row), agg_nl (aggressive-coarsening levels), relax_type_all (smoother, e.g. "SOR/Jacobi"/"Chebyshev"), grid_sweeps_all (smoother sweeps). See the PETSc -pc_hypre_boomeramg_* options for the full list; anything else can go through petsc_options.

FullPrecision()

Default SolverSetup precision: the linear system is solved in the mesh's float type.

GAMG(; kwargs...) <: PreconditionerType

PETSc native aggregation AMG (-pc_type gamg) on distributed meshes; on a serial mesh XCALibre's native AMG runs instead, with a one-time warning (GAMG's keywords are not carried over). No hypre build required. SPD systems only (e.g. pressure). The recommended AMG for distributed pressure solves on large meshes, on CPU and GPU. It needs more memory than Jacobi and pays off only when each rank holds a large partition. The hierarchy is built per partition, so residuals vary slightly with the rank count; use Jacobi where results must match across rank counts exactly.

Defaults set reuse_interpolation=true: because the mesh never refines, the pressure matrix keeps a FIXED sparsity pattern, so GAMG builds the aggregation + prolongation P once and recomputes only the coarse operators (RAP) and smoothers each solve — the hierarchy stays numerically current at a fraction of a full setup. freeze=N additionally holds the whole preconditioner fixed for N solves, skipping even that update, and rebuilds it on the Nth (default 25; freeze=1 updates every solve). Because the coefficient drift between rebuilds is small, the freeze leaves the pressure residual essentially unchanged while removing most of the setup cost.

Each keyword k=v overrides a default and is passed as -pc_gamg_<k> v, e.g. GAMG(threshold=0.02, square_graph=1). See the PETSc -pc_gamg_* options.

Geometric(; merge_levels=1)

OpenFOAM-style geometric agglomeration (faceAreaPair): cells are merged pairwise into coarse clusters by greatest face weight, using |a_ij| as the face weight (proportional to face-area/delta in FVM). merge_levels sets the number of pairwise passes per coarse level, giving cluster size ~2^merge_levels. Prolongation is unsmoothed piecewise-constant injection (additive correction); coarse operators are formed algebraically by Galerkin RAP.

struct JacobiSmoother{L,F,V} <: AbstractSmoother
    loops::L
    omega::F
    x_temp::V
end

Structure to hold information for using the weighted Jacobi smoother.

Fields

• loops represents the number of smoothing iterations.
• omega represents the relaxation weight, 1 corresponds to no weighting. Typically a weight of 2/3 is used to filter high frequencies in the residual field.
• x_temp is a vector used internally to store intermediate solutions.

Constructors

JacobiSmoother(mesh::AbstractMesh)

This constructor takes a mesh object as input to create the smoother. The number of cells in the mesh is used to allocate x_temp. The number of loops and smoother factor are set to 5 and 1, respectively.

JacobiSmoother(; domain, loops, omega=2/3)

Alternative constructor for JacobiSmoother.

keyword arguments

• domain represents a mesh object of type AbstractMesh.
• loops is the number of iterations to be used
• omega represents the weighting factor, 1 does not relax the system, 2/3 is found to work well for smoothing high frequencies in the residual field

LaggedTerms(eqn, fields; evaluate=nothing, levels=())
extrapolate_lagged!(lagged::LaggedTerms, time, config)

Explicitly treated (lagged) fields of a transient solver, such as the convecting face flux, the effective viscosity and the velocity gradient, which the solver recomputes later in each step. A lag => field pair is a lagged copy of a transported field. Called once at the start of each step, extrapolate_lagged! sets them to the time level the time scheme of eqn integrates at: Euler keeps their values from the end of the previous step; Crank-Nicolson extrapolates them to tⁿ⁺¹ from the two previous steps, with weights from the time-step history, after evaluate() has set them from the initial state on the first step. levels lists (eqn, field, rho) of equations whose Time term takes its old level at the start of the step, before the solver changes rho (such as a pressure equation with ψ updated by the energy equation).

LocalTimeStepping(; maxCo=1.0, startCo=maxCo, ramp=0, turbulence_dt=:min)

Local (per-cell) pseudo-time stepping for the explicit density-based solver (Fluid{SupersonicFlow}, godunov!). Pass it as adaptive to Runtime to converge steady problems faster: every cell advances with its own time step dt_i = CFL / (0.5 Σ_f (|U_i·n_f| + a_i) A_f / V_i + ν_eff Σ_f A_f² / V_i²), so the result of a converged run does not depend on dt and the iterations are not a time history. It combines with FEuler(), RK2() and the multistage schemes RK3(), RK4(), RK5().

Input arguments

• maxCo::AbstractFloat: Courant number of the cells (final value of the ramp).
• startCo::AbstractFloat: Courant number of the first iteration.
• ramp::Integer: number of iterations over which the Courant number grows linearly from startCo to maxCo (0 = no ramp). Useful to start from a uniform field.
• turbulence_dt::Symbol: the implicit turbulence equations take one time step, the :min (default), :mean or :max of the cell time steps of the iteration.

MixedF16()

As MixedF32 with the off-diagonal matrix entries stored in Float16, which halves the matrix traffic again; the Krylov iterations run in Float32. The system is scaled symmetrically to unit diagonal to fit Float16's range, the diagonal is kept in Float32 and compensated so the rounded matrix keeps the exact row sums of the scaled system, and each correction is applied with a minimal-residual step length computed in full precision. A single correction reduces the residual by about 1e-3 at most, so tighter rtol takes several corrections.

Serial meshes and the Jacobi preconditioner only: the scaled matrix has a unit diagonal, so Jacobi is applied as part of the scaling at no extra cost, while factorisation preconditioners would need their own higher-precision copy of the matrix, removing the saving. Useful on GPUs; CPUs without native Float16 arithmetic run it slower than MixedF32.

MixedF32()

SolverSetup precision that solves each linear system for its correction in Float32. The residual r = b - Ax and the update x += d are formed in the mesh's float type, so the solution keeps full accuracy while the Krylov iterations, matrix and preconditioner use Float32 storage, which halves the memory traffic of each iteration. A single Float32 solve cannot reduce the residual by much more than 1e-5; when rtol asks for more, the correction is repeated (iterative refinement) until the full-precision residual meets rtol, so tight tolerances, e.g. in transient runs, are met as in full precision.

Supported with the Krylov solvers (Cg, Cgs, Bicgstab, Gmres) and every serial preconditioner, which is built on the Float32 matrix. On distributed meshes PETSc solves the correction with its Float32 library (PETSc_jll loads it beside the Float64 one), so PETSc preconditioners (e.g. GAMG) are available too.

OnDevice(; max_rows=512)

Solve the coarsest level on the device (no host round-trip): a device-resident factorization (Cholesky/LU on CUDA; generic backends use a device dense inverse + GEMV fallback). Eliminates a per-cycle fine→coarse host sync/transfer (the win at MPI/multi-GPU scale). The factor is rebuilt each coarse refresh (~n^3), so coarsest levels larger than max_rows fall back to OnHost — max_rows=512 is a conservative no-regression default (well-coarsened 3D reaches a tiny coarsest). On a CPU backend this is equivalent to OnHost.

OnDeviceChebyshev(; degree=10, eig_ratio=10.0, lambda_scale=1.1)

GPU coarse solve: apply a FIXED-degree Chebyshev polynomial smoother (x init 0) to the coarsest level. Like OnDeviceJacobi it is a constant linear operator (b-independent), so valid with solver=Cg() as well as solver=Richardson(); fully on device. Chebyshev targets the spectrum [lambda_max/eig_ratio, lambda_max] (lambdamax estimated per level; `lambdascaleinflates the upper bound), giving faster coarse error reduction per matvec than Jacobi at equal cost. On a CPU backend this is equivalent toOnHost`.

OnDeviceJacobi(; omega=4/3, iterations=10)

GPU coarse solve: apply a FIXED number of weighted-Jacobi sweeps (x init 0) to the truncated, large coarsest level instead of a direct factor. A fixed-sweep Jacobi is a constant linear operator p(A_c)·b, so — unlike OnDeviceKrylov — it is valid as an outer-CG preconditioner: allowed with solver=Cg() as well as solver=Richardson(). Fully on device (no host round-trip), removing the coarse host sync point at MPI/multi-GPU scale. Cheaper per cycle than a direct solve but only approximate, so outer iterations may rise; best when the coarsest is large and well-conditioned (huge coarsening ratio). On a CPU backend this is equivalent to OnHost.

OnDeviceKrylov(; solver=nothing, rtol=1e-2, atol=0.0, itmax=50)

GPU strategy: truncate the hierarchy at a LARGE coarsest level (set max_coarse_rows high, e.g. 4000-8000) so every level has enough work for good GPU occupancy, then solve that coarsest level with an in-place Krylov solver (Cg/Bicgstab) + Jacobi preconditioner — the configuration GPUs parallelise best. Reuses XCALibre's Krylov.jl infrastructure (no host round-trip, no extra alloc). solver=nothing auto-picks Cg (symmetric coarsest) or Bicgstab. Requires SolverSetup(solver=Richardson()): an inner Krylov solve is nonlinear in b, which breaks PCG's fixed-SPD-preconditioner assumption. On a CPU backend this is equivalent to OnHost.

OnHost()

Solve the coarsest level on the host (sparse LU/QR). On a GPU backend the per-cycle coarse solve copies the rhs to host and the solution back — fastest single-device for a largish coarsest.

Richardson()

Stationary iteration x += M⁻¹(b - Ax) with no Krylov acceleration. Use it with preconditioner = AMG(...) for a standalone multigrid solve (repeated V-cycles); on a distributed mesh it maps to PETSc's richardson.

Runtime(; 
        # keyword arguments

        iterations::I, 
        write_interval::I, 
        time_step::N,
        adaptive::A
    ) where {I<:Integer,N<:Number} = begin
    
    # returned Runtime struct
    Runtime{I<:Integer,F<:AbstractFloat}
        (
            iterations=iterations, 
            dt=time_step, 
            write_interval=write_interval,
            adaptive=adaptive
        )
end

This is a convenience function to set the top-level runtime information. The inputs are all keyword arguments and provide basic information to flow solvers just before running a simulation.

Input arguments

• iterations::Integer: specifies the number of iterations in a simulation run.
• write_interval::Integer: defines how often simulation results are written to file (on the current working directory). The interval is currently based on number of iterations. Set to -1 to run without writing results to file.
• time_step::AbstractFloat: the time step to use in the simulation. Notice that for steady solvers this is simply a counter and it is recommended to simply use 1.
• adaptive::Union{Nothing, AdaptiveTimeStepping}: optionally enables adaptive time stepping. Pass an AdaptiveTimeStepping object to automatically adjust dt based on the Courant number during transient simulations. Defaults to nothing, meaning a fixed time step is used.

Example

runtime = Runtime(iterations=2000, time_step=1, write_interval=2000)

Schemes(;
    # keyword arguments and their default values
    time=SteadyState,
    divergence=Linear, 
    laplacian=Linear, 
    gradient=Gauss,
    limiter=nothing) = begin

    # Returns Schemes struct used to configure discretisation
    
    Schemes(
        time=time,
        divergence=divergence,
        laplacian=laplacian,
        gradient=gradient,
        limiter=limiter
    )   
end

The Schemes struct is used at the top-level API to help users define discretisation schemes for every field solved.

Inputs

• time: is used to set the time schemes (default is SteadyState)
• divergence: is used to set the divergence scheme (default is Linear)
• laplacian: is used to set the laplacian scheme (default is Linear)
• gradient: is used to set the gradient scheme (default is Gauss)
• limiter: is used to specify if gradient limiters should be used, currently supported limiters include CellBased, Venkatakrishnan, FaceBased and MFaceBased (default is nothing)

SolverSetup(; 
        # required keyword arguments 

        solver::S, 
        preconditioner::PT, 
        convergence, 
        relax,

        # optional keyword arguments

        float_type=Float64,
        smoother=nothing,
        limit=nothing,
        itmax::Integer=1000, 
        atol=(eps(float_type))^0.9,
        rtol=float_type(1e-1),
        precision=FullPrecision()

    ) where {S,PT<:PreconditionerType} = begin

        return SolverSetup(kwargs...)  
end

This function is used to provide solver settings that will be used internally in XCALibre.jl. It returns a SolverSetup object with solver settings that are used internally by the flow solvers.

Input arguments

• solver: the outer iteration, one of Cg(), Cgs(), Bicgstab(), Gmres() (Krylov.jl, re-exported in XCALibre.jl) or Richardson() (with AMG only). Gmres() restarts every 30 iterations.
• preconditioner: the method, e.g. Jacobi(), DILU(), AMG(), GAMG(). Where a preconditioner is not available (e.g. AMG on a distributed mesh) a named substitute runs with a one-time warning; see the preconditioners section of the numerical setup page
• convergence: residual target for this field that stops the outer (e.g. SIMPLE) iteration; it does not control the linear solver, which stops on atol, rtol and itmax. The residual, Σ|b - Ax| / Σ|Aᵢᵢxᵢ|, is taken on the field before the linear solve, so it does not depend on the solver settings.
• relax: specifies the relaxation factor to be used e.g. set to 1 for no relaxation
• smoother: specifies smoothing method to be applied before discretisation. JacobiSmoother: is currently the only choice (defaults to nothing)
• limit: used in some solvers to bound the solution within these limits e.g. (min, max). It defaults to nothing
• itmax: maximum number of iterations in a single solver pass (defaults to 1000). Also applies to PETSc solves.
• atol: absolute tolerance of the linear solve, in the field's units (defaults to eps(float_type)^0.9).
• rtol: relative tolerance of the linear solve (defaults to 1e-1). Every solver, preconditioner, precision and backend (including AMG and PETSc) stops when ‖D⁻¹(b - Ax)‖₂ ≤ atol + rtol·‖D⁻¹(b - Ax₀)‖₂, D = |diag(A)|: the size of a Jacobi correction, with rtol relative to that of the initial guess x₀ (the current field values), not to ‖b‖.
• float_type: specifies the floating point type to be used by the solver. It is also used to estimate the absolute tolerance for the solver (defaults to Float64)
• precision: precision of the linear solve, FullPrecision() (default), MixedF32() or MixedF16().

MeshWave <: AbstractWallDistance

Geometric wall distance (default). Each cell carries the location of its nearest wall point and the Euclidean distance to it. Cells touching a wall node start from the exact closest point on the wall faces through their nodes; sweeps over the internal faces then let a cell take a neighbour's wall point whenever that point is closer, until no cell changes. Each sweep only reads face neighbours, so it runs unchanged on distributed meshes (one halo exchange per sweep) and the result does not depend on the partitioning.

Example

• RANS{KOmegaSST}(walls=(:wall,), wall_distance=MeshWave())

Poisson(; iterations=1000) <: AbstractWallDistance

Wall distance from the solution φ of ∇²φ = -1 (φ = 0 on walls), y = -|∇φ| + √(|∇φ|² + 2φ). It needs solvers.y and schemes.y in the configuration. Exact for a single plane wall, it loses accuracy away from walls, e.g. upstream of a leading edge.

Example

• RANS{KOmegaSST}(walls=(:wall,), wall_distance=Poisson(iterations=2000))

Venkatakrishnan(K=5.0)

Cell-based gradient limiter with the smooth limiter function of Venkatakrishnan (1995), selected with Schemes(limiter=Venkatakrishnan(K)). Like CellBased it scales the cell gradient by one factor in [0, 1] so that the values extrapolated to the faces stay within the minimum and maximum of the cell and its face neighbours, but the factor is a smooth function of the ratio of allowed to extrapolated change instead of min(1, ratio). Limiting therefore starts gradually (the factor is 1 only from a ratio of 2) and does not switch on and off from one iteration to the next.

K ≥ 0 controls the threshold ε² = (K h)³ below which differences are not limited, with h the cell length scale (cube root of the volume, square root of the area on 2D meshes). ε² is compared with squared differences of the limited field, so its effect depends on the units of the field: for velocities and pressures in SI units and cell sizes well below 1 m it is negligible unless K is large. K = 0 limits all differences.

div!(phi::VectorField, Γf, tensor, BCs, config; Γc=nothing)

Cell values of ∇·(Γ tensor): phi[i] = (1/V_i) Σ_f Γf[f] (tensor_f ⋅ S_f), where tensor is any cell tensor field indexed as tensor[i] (including lazy forms such as Dev2(T(gradU))) and S_f is the outward face area vector. At internal faces tensor_f is the linear interpolation of the two cell values with the mesh weights; at boundary faces it is the owner-cell value. With a cell coefficient Γc (any cell-indexable value, e.g. ν + ν_t) internal faces take the interpolated product, w Γc[P] tensor[P] + (1 - w) Γc[N] tensor[N], i.e. the face value of the cell flux tensor Γ·tensor; boundary faces keep Γf[f] times the owner value, and periodic faces take the internal-face value between the two cells of the pair. Interpolating Γ and the tensor separately instead lets a large neighbour viscosity multiply the steep gradient of a thin near-wall cell (a face weight close to 1 on the thin side), which makes an explicit source built this way unstable. Faces of Slip, Symmetry and Empty boundaries carry no flux. BCs are the boundary conditions of the field tensor derives from (e.g. boundaries.U). Ghost cells of a distributed mesh must hold current values of tensor.

interpolate!(phif::FaceScalarField, phi::ScalarField, mdotf::FaceScalarField, scheme::Divergence, config)

Internal face values of phi convected by the face flux mdotf with the face interpolation of the divergence scheme scheme, e.g. Divergence{LUST}(), so that mdotf*phif is the face flux the implicit operator Divergence{LUST}(mdotf, phi) assembles. Upwind and BoundedUpwind take the upwind cell value, Linear the linear interpolation, and LUST and LUSTDC 75% linear and 25% upwind. Boundary face values are left unchanged. Any other scheme throws an ArgumentError; the energy equation reconstructs K for LinearUpwind and LimitedLinearUpwind itself, with the gradient of K.

total_volume(mesh, config) → Tf

Sum of all cell volumes in the mesh domain. Runs on the mesh backend.

volume_average(phi::ScalarField, config) → Tf
volume_average(phi::VectorField, config) → Vector{3}

Volume-averaged mean: (∫ phi dV) / (∫ dV).

volume_integral(phi::ScalarField, config) → Tf
volume_integral(phi::VectorField, config) → Vector{3}

Volume integral of a field: ∫ phi dV (component-wise for vector fields). Runs on the field's backend (CPU or GPU).

wall_distance!(model, walls, config; method=MeshWave())

Compute the wall distance model.turbulence.y to the patches named in walls with method (MeshWave() or Poisson()). Returns a new Configuration whose boundaries include the conditions of y.

weighted_volume_integral(phi::ScalarField, weight_func, config) → Tf
weighted_volume_integral(phi::VectorField, weight_func, config) → Vector{3}

∫ phi(x) * w(x, y, z) dV where w = weight_func(x, y, z) is evaluated at each cell centroid (component-wise for vector fields). Runs on the field's backend. On GPU backends weight_func is compiled into the kernel, so it must not allocate, and it must not read non-const globals or capture arrays. Capture parameters with let (e.g. w = let x0 = 0.5; (x, y, z) -> x - x0 end) or use const globals; a named function needs no @inline.

struct Momentum{V,S,SS} <: AbstractMomentumModel
    U::V 
    p::S 
    sources::SS
end

Momentum model containing key momentum fields.

Fields

• 'U' – Velocity VectorField.
• 'p' – Pressure ScalarField.
• 'sources' – Momentum model sources.

Examples

• `Momentum(mesh::AbstractMesh)

Compressible <: AbstractCompressible

Compressible fluid model containing fluid field parameters for compressible flows with constant parameters - ideal gas with constant viscosity.

Fields

• 'nu' – Fluid kinematic viscosity.
• 'cp' – Fluid specific heat capacity.
• gamma – Ratio of specific heats.
• Pr – Fluid Prantl number.

Examples

• Fluid{Compressible}(; nu=1E-5, cp=1005.0, gamma=1.4, Pr=0.7) - Constructur with default values.

Conduction <: AbstractEnergyModel

Type that represents energy conduction model for solids.

Fields

• T Temperature scalar field on cells (ScalarField).
• Tf Temperature scalar field on faces (FaceScalarField).

ConstEos <: AbstractEosModel

Constant density equation of state model.

Fields

• 'rho' – Constant density value.

ConstMu <: AbstractViscosityModel

Constant dynamic viscosity model.

Fields

• 'mu' – Dynamic viscosity value [Pa⋅s].

Fluid <: AbstractFluid

Abstract fluid model type for constructing new fluid models.

Fields

• 'args' – Model arguments.

Stores all fluid-specific constants required for the Helmholtz Equation of State, ancillary equations, and associated property calculations.

Incompressible <: AbstractIncompressible

Incompressible fluid model containing fluid field parameters for incompressible flows.

Fields

• 'nu' – Fluid kinematic viscosity.
• 'rho' – Fluid density.

Examples

• Fluid{Incompressible}(nu=0.001, rho=1.0) - Constructor with default values.

Incompressible_MRF <: AbstractIncompressible

Incompressible fluid model containing fluid field parameters for incompressible flows that utilise multiple reference frames (MRF).

Fields

• 'nu' – Fluid kinematic viscosity.
• 'rho' – Fluid density.
• 'frames' – Reference frames information.

Examples

• Fluid{Incompressible}(nu=0.001, rho=1.0) - Constructor with default values.

InternalEnergy <: AbstractEnergyModel

Type that represents the internal energy transport model, coefficients and respective fields. The solved variable is e = cv*(T - Tref) where cv = cp/gamma.

Fields

• he: Internal energy ScalarField (e = cv*(T - Tref)).
• T: Temperature ScalarField.
• hef: Internal energy FaceScalarField.
• Tf: Temperature FaceScalarField.
• K: Specific kinetic energy ScalarField.
• S_he: Energy source term ScalarField (-p*div(U) for internal energy).
• coeffs: A tuple of model coefficients.

KEquation <: AbstractTurbulenceModel

KEquation LES model containing all Smagorinksy field parameters.

Fields

• nut – Eddy viscosity ScalarField.
• nutf – Eddy viscosity FaceScalarField.
• coeffs – Model coefficients.

KOmega <: AbstractTurbulenceModel

kOmega model containing all kOmega field parameters.

Fields

• k – Turbulent kinetic energy ScalarField.
• omega – Specific dissipation rate ScalarField.
• nut – Eddy viscosity ScalarField.
• nutf – Eddy viscosity FaceScalarField.
• coeffs – Model coefficients.

KOmegaLKE <: AbstractTurbulenceModel

kOmega model containing all kOmega field parameters.

Fields

• k – Turbulent kinetic energy ScalarField.
• omega – Specific dissipation rate ScalarField.
• kl – ScalarField.
• nut – Eddy viscosity ScalarField.
• kf – Turbulent kinetic energy FaceScalarField.
• omegaf – Specific dissipation rate FaceScalarField.
• klf – FaceScalarField.
• nutf – Eddy viscosity FaceScalarField.
• coeffs – Model coefficients.
• Tu – Freestream turbulence intensity for model.
• y – Near-wall distance for model.

KOmega <: AbstractTurbulenceModel

kOmega model containing all kOmega field parameters.

Fields

• 'k' – Turbulent kinetic energy ScalarField.
• 'omega' – Specific dissipation rate ScalarField.
• 'nut' – Eddy viscosity ScalarField.
• 'kf' – Turbulent kinetic energy FaceScalarField.
• 'omegaf' – Specific dissipation rate FaceScalarField.
• 'nutf' – Eddy viscosity FaceScalarField.
• 'coeffs' – Model coefficients.

LES <: AbstractLESModel

Abstract LES model type for constructing LES models.

Fields

• 'args' – Model arguments.

Laminar <: AbstractTurbulenceModel

Laminar model definition for physics API.

Mixture(; diameter=1.0e-3) <: AbstractMultiphaseModel

Manninen drift-flux mixture-model settings.

Fields

• diameter – Dispersed-phase particle/bubble diameter [m].

Multiphase <: AbstractMultiphase

Multiphase fluid model containing multiple phases and their interaction properties.

Fields

• 'model' – Multiphase model selecting the solver pathway (VOF or Mixture).
• 'phases' – Tuple of PhaseState structures.
• 'physics_properties' – NamedTuple of physical models (drag, surface tension, etc.).
• 'volume_fraction' – Index of the phase tracked by the volume fraction field.
• 'alpha' – Volume fraction ScalarField.
• 'alphaf' – Volume fraction FaceScalarField.
• 'rho' – Mixture density ScalarField.
• 'rhof' – Mixture density FaceScalarField.
• 'nu' – Mixture kinematic viscosity ScalarField.
• 'nuf' – Mixture kinematic viscosity FaceScalarField.
• 'p_rgh' – Dynamic pressure ScalarField.
• 'p_rghf' – Dynamic pressure FaceScalarField.

Phase <: AbstractPhase

Configuration structure for a single fluid phase.

Fields

• rho – Density model (Equation of State) for the phase.
• mu – Viscosity model for the phase.

struct Physics{T,F,SO,M,Tu,E,D,BI}
    time::T
    fluid::F
    solid::SO
    momentum::M
    turbulence::Tu
    energy::E
    domain::D
    boundary_info::BI
    wall_info
end

XCALibre's parametric Physics type for user-level API. Also used to dispatch flow solvers.

Fields

• time::Union{Steady, Transient} – User-provided time model.
• fluid::AbstractFluid – User-provided Fluid` model.
• solid::AbstractSolid – User-provided Solid` model.
• momentum – Momentum model. Currently this is auto-generated by the Physics constructor
• turbulence::AbstractTurbulenceModel – User-provided Turbulence` model.
• energy:AbstractEnergyModel – User-provided Energy model.
• domain::AbstractMesh – User-provided Mesh. Must be adapted to target device before constructing a Physics object.
• boundary_info::boundary_info – Mesh boundary information. Auto-generated by the Physics constructor
• wall_info – Wall boundary names and wall distance method ((walls, method)) for turbulence models requiring wall distance. Auto-generated by the Physics constructor

Physics(; time, fluid, solid, turbulence, energy, domain)::Physics{T,F,SO,M,Tu,E,D,BI}

Physics constructor part of the top-level API. It can be used to define the Physics and models relevant to a simulation. This constructor uses keyword arguments to allow users to fine-tune their simulations, whilst some fields are auto-generated behind the scenes for convenience (Momentum and boundary_info). Where:

• time - specified the time model (Steady or Transient)
• fluid - specifies the type of fluid (Incompressible, etc.)
• turbulence - specified the Turbulence model
• energy - specifies the Energy treatment
• domain - provides the mesh to used (must be adapted to the target backend device)

RANS <: AbstractRANSModel

Abstract RANS model type for consturcting RANS models.

Fields

• 'args' – Model arguments.

SensibleEnthalpy <: AbstractEnergyModel

Type that represents energy model, coefficients and respective fields.

Fields

• he: Sensible enthalpy ScalarField.
• T: Temperature ScalarField.
• hef: Sensible enthalpy FaceScalarField.
• Tf: Temperature FaceScalarField.
• K: Specific kinetic energy ScalarField.
• S_he: Energy source term ScalarField (dp/dt for enthalpy).
• coeffs: A tuple of model coefficients.

Smagorinsky <: AbstractTurbulenceModel

Smagorinsky LES model containing all Smagorinksy field parameters.

Fields

• nut – Eddy viscosity ScalarField.
• nutf – Eddy viscosity FaceScalarField.
• coeffs – Model coefficients.

Solid <: AbstractSolid

Abstract solid model type for constructing new solid models.

Fields

• 'args' – Model arguments.

Steady

Steady model for Physics model API.

Examples

• Steady()

SupersonicFlow <: AbstractCompressible

Fluid model for density-based (explicit) supersonic flow solver. Uses Rusanov (Local Lax-Friedrichs) flux with Forward Euler time integration.

Fields

• nu – Kinematic viscosity (ConstantScalar).
• rho – Density field (ScalarField, updated each iteration).
• nuf – Face kinematic viscosity.
• rhof – Face density field.
• cp – Specific heat at constant pressure (ConstantScalar).
• gamma – Ratio of specific heats (ConstantScalar).
• Pr – Prandtl number (ConstantScalar).
• R – Specific gas constant cp*(1 - 1/gamma) (ConstantScalar).

Examples

• Fluid{SupersonicFlow}(nu=1E-5, cp=1005.0, gamma=1.4, Pr=0.7) - Constructor with default values.

Transient

Transient model for Physics model API.

Examples

• Transient()

VOF(; sigma=0.0, cAlpha=1.0) <: AbstractMultiphaseModel

Volume-of-Fluid interface-capturing settings.

Fields

• sigma – Surface tension coefficient [N/m].
• cAlpha – Interface compression coefficient (MULES), default is 1.0.

WALE <: AbstractTurbulenceModel

Wall-Adapting Local Eddy-viscosity (WALE) LES model (Nicoud & Ducros, 1999) containing all WALE field parameters. The eddy viscosity is computed algebraically as

nut = (C*Δ)^2 * (Sd:Sd)^(3/2) / ((S:S)^(5/2) + (Sd:Sd)^(5/4))

where S is the strain rate tensor and Sd is the traceless symmetric part of the squared velocity gradient tensor. The default model coefficient is C = 0.325.

Fields

• nut – Eddy viscosity ScalarField.
• nutf – Eddy viscosity FaceScalarField.
• coeffs – Model coefficients.

WeaklyCompressible <: AbstractCompressible

Weakly compressible fluid model containing fluid field parameters for weakly compressible flows with constant parameters - ideal gas with constant viscosity.

Fields

• 'nu' – Fluid kinematic viscosity.
• 'cp' – Fluid specific heat capacity.
• gamma – Ratio of specific heats.
• Pr – Fluid Prandtl number.

Examples

• Fluid{WeaklyCompressible}(; nu=1E-5, cp=1005.0, gamma=1.4, Pr=0.7) - Constructor with

default values.

Main high-level wrapper. Determines the fluid phase (liquid, vapor, two-phase, or supercritical) and computes a comprehensive set of thermodynamic properties for the given state (T, P).

Computes the isobaric expansion coefficient in 1/K.

Computes the isobaric (constant pressure) heat capacity in J/(mol*K).

Computes the isochoric (constant volume) heat capacity in J/(mol*K).

change(model::Physics, property, value) => updatedModel::Physics

A convenience function to change properties of an exisitng Physics model.

Input arguments

• model::Physics a Physics model to modify
• property is a symbol specifying the property to change
• value is the new setting for the specified property

Output

This function return a new Physics object

Example

To change a model to run a transient simulation e.g. after converging in steady state

modelTransient = change(model, :time, Transient())

Computes pressure (in Pa) for a given temperature (K) and molar density (mol/m^3).

Computes the ideal-gas isochoric heat capacity (Cv0) in J/(mol*K) at a given temperature.

Calculates the derivative of saturation pressure with respect to temperature (dPsat/dT) using the ancillary equation.

eddy_viscosity!(turbulenceModel, model, S, time, config)

Sets ν_t (cell and face) from the current fields without advancing any transport equation.

energy!(energy::EnergyEquationModel, model::Physics{T1,F,SO,M,Tu,E,D,BI}, prevP, prevRhoK, mdotf, gradP, gradU, rho, mueff, time, config
) where {T1,F,SO,M,Tu,E,D,BI}

Run energy transport equations (sensible enthalpy or internal energy).

Input

• energy: EnergyEquationModel.
• model: Physics model defined by user.
• prevP: Previous pressure cell values.
• prevRhoK: Previous rho*K values.
• mdotf: Face mass flow.
• gradP: Pressure gradient.
• gradU: Velocity gradient.
• rho: Density ScalarField.
• mueff: Effective viscosity FaceScalarField.
• time: Simulation runtime.
• config: Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Computes the total enthalpy (ideal + residual) in J/mol.

Computes the total entropy (ideal + residual) in J/(mol*K).

Calculates saturation properties (pressure, liquid density, vapor density) at a given temperature. Uses a secant solver to find the pressure where the liquid and vapor Gibbs free energies are equal.

Iteratively solves for the saturation temperature (in K) for a given pressure (Pa) using the ancillary equation.

Computes the Gibbs free energy in J/mol for a given state.

Computes the ideal-gas enthalpy (h0) in J/mol at a given temperature.

initialise(turbulence::KEquation, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
) where {T,F,SO,M,Tu,E,D,BI}

Initialisation of turbulent transport equations.

Input

• turbulence – turbulence model.
• model – Physics model defined by user.
• mdtof – Face mass flow.
• peqn – Pressure equation.
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

• KEquationModel( turbulence, Δ, magS, ModelState((), false) ) – Turbulence model structure.

initialise(turbulence::KOmega, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
) where {T,F,SO,M,Tu,E,D,BI}

Initialisation of turbulent transport equations.

Input

• turbulence – turbulence model.
• model – Physics model defined by user.
• mdtof – Face mass flow.
• peqn – Pressure equation.
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

• KOmegaModel( turbulence, k_eqn, ω_eqn, state ) – Turbulence model structure.

initialise(turbulence::KOmegaLKE, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
) where {T,F,SO,M,Tu,E,D,BI}

Initialisation of turbulent transport equations.

Input

• turbulence – turbulence model.
• model – Physics model defined by user.
• mdtof – Face mass flow.
• peqn – Pressure equation.
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

• KOmegaLKEModel( turbulence, k_eqn, ω_eqn, kl_eqn, nueffkLS, nueffkS, nueffωS, nuL, nuts, Ω, γ, fv, normU, Reυ, divU, S2, ReLambda, ∇k, ∇ω, state ) – Turbulence model structure.

initialise(turbulence::KOmega, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
) where {T,F,SO,M,Tu,E,D,BI}

Initialisation of turbulent transport equations.

Input

• turbulence – turbulence model.
• model – Physics model defined by user.
• mdtof – Face mass flow.
• peqn – Pressure equation.
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

• KOmegaModel(k_eqn, ω_eqn) – Turbulence model structure.

function initialise(
    turbulence::Laminar, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
    ) where {T,F,SO,M,Tu,E,D,BI}
return LaminarModel(), config

end

Initialisation of turbulent transport equations.

Input

• turbulence – turbulence model.
• model – Physics model defined by user.
• mdtof – Face mass flow.
• peqn – Pressure equation.
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

• LaminarModel() – Turbulence model structure.

initialise(turbulence::Smagorinsky, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
) where {T,F,SO,M,Tu,E,D,BI}

Initialisation of turbulent transport equations.

Input

• turbulence: turbulence model.
• model: Physics model defined by user.
• mdtof: Face mass flow.
• peqn: Pressure equation.
• config: Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

Returns a structure holding the fields and data needed for this model

SmagorinskyModel(
    turbulence, 
    Δ, 
    ModelState((), false)
)

initialise(turbulence::WALE, model::Physics{T,F,SO,M,Tu,E,D,BI}, mdotf, peqn, config
) where {T,F,SO,M,Tu,E,D,BI}

Initialisation of turbulent transport equations.

Input

• turbulence: turbulence model.
• model: Physics model defined by user.
• mdtof: Face mass flow.
• peqn: Pressure equation.
• config: Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Output

Returns a structure holding the fields and data needed for this model

WALEModel(
    turbulence,
    Δ,
    ModelState((), false),
    wall_scratch
)

Computes the total internal energy (ideal + residual) in J/mol.

Computes the isothermal compressibility in 1/Pa (required for thermal conductivity computation).

Computes a set of key thermodynamic properties from molar density and temperature, and converts them from molar to mass-specific units.

turbulence!(rans::KOmegaLKEModel, model::Physics{T,F,SO,M,Turb,E,D,BI}, S, prev, time, config ) where {T,F,SO,M,Turb<:AbstractTurbulenceModel,E,D,BI}

Run turbulence model transport equations.

Input

• rans::KOmegaLKEModel – KOmega turbulence model.
• model – Physics model defined by user.
• S – Strain rate tensor.
• prev – Previous field.
• time –
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

turbulence!(les::KEquationModel, model::Physics{T,F,SO,M,Tu,E,D,BI}, S, S2, prev, time, config
) where {T,F,SO,M,Tu<:AbstractTurbulenceModel,E,D,BI}

Run turbulence model transport equations.

Input

• les::KEquationModel – KEquation LES turbulence model.
• model – Physics model defined by user.
• S – Strain rate tensor.
• S2 – Square of the strain rate magnitude.
• prev – Previous field.
• time –
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

turbulence!(rans::KOmegaModel, model::Physics{T,F,SO,M,Tu,E,D,BI}, S, prev, time, config
) where {T,F,SO,M,Tu<:AbstractTurbulenceModel,E,D,BI}

Run turbulence model transport equations.

Input

• rans::KOmegaModel – KOmega turbulence model.
• model – Physics model defined by user.
• S – Strain rate tensor.
• prev – Previous field.
• time –
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

turbulence!(rans::LaminarModel, model::Physics{T,F,SO,M,Tu,E,D,BI}, S, prev, time, config
) where {T,F,SO,M,Tu<:Laminar,E,D,BI}

Run turbulence model transport equations.

Input

• rans::LaminarModel – Laminar turbulence model.
• model – Physics model defined by user.
• S – Strain rate tensor.
• prev – Previous field.
• time –
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

turbulence!(les::SmagorinskyModel, model::Physics{T,F,SO,M,Tu,E,D,BI}, S, S2, prev, time, config
) where {T,F,SO,M,Tu<:AbstractTurbulenceModel,E,D,BI}

Run turbulence model transport equations.

Input

• les::SmagorinskyModel: Smagorinsky LES turbulence model.
• model: Physics model defined by user.
• S: Strain rate tensor.
• S2: Square of the strain rate magnitude.
• prev: Previous field.
• time: current simulation time
• config: Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

turbulence!(les::WALEModel, model::Physics{T,F,SO,M,Tu,E,D,BI}, S, prev, time, config
) where {T,F,SO,M,Tu<:AbstractTurbulenceModel,E,D,BI}

Update the WALE eddy viscosity (algebraic model, no transport equations are solved).

Input

• les::WALEModel: WALE LES turbulence model.
• model: Physics model defined by user.
• S: Strain rate tensor.
• prev: Previous field.
• time: current simulation time
• config: Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

turbulence!(rans::KOmegaModel{E1,E2,S1}, model::Physics{T,F,SO,M,Tu,E,D,BI}, S, prev, time, config
) where {T,F,SO,M,Tu<:KOmega,E,D,BI,E1,E2,S1}

Run turbulence model transport equations.

Input

• rans::KOmegaModel{E1,E2,S1} – KOmega turbulence model.
• model – Physics model defined by user.
• S – Strain rate tensor.
• prev – Previous field.
• time –
• config – Configuration structure defined by user with solvers, schemes, runtime and hardware structures set.

Computes the ideal-gas internal energy (u0) in J/mol at a given temperature.

Calculates the saturation pressure (in Pa) for a given temperature using a simplified ancillary equation.

struct Configuration{SC,SL,RT,HW,PP}
    schemes::SC
    solvers::SL
    runtime::RT
    hardware::HW
    postprocess::PP
end

The Configuration type is passed to all flow solvers and provides all the relevant information to run a simulation.

Inputs

• schemes::NamedTuple this keyword argument is used to pass distretisation scheme information to flow solvers. See Numerical setup for details.
• solvers::NamedTuple this keyword argument is used to pass the configurations for the linear solvers for each field information to flow solvers. See Runtime and solvers for details.
• runtime::NamedTuple this keyword argument is used to pass runtime information to the flow solvers. See Runtime and solvers for details.
• hardware::NamedTuple this keyword argument is used to pass the hardware configuration and backend settings to the flow solvers. See Pre-processing for details.

Optional keywords

• postprocess this keyword argument is used to pass any fields that need to be post-processed. See Post-processing for details.

Example

config = Configuration(
    solvers=solvers, schemes=schemes, runtime=runtime, hardware=hardware, boundaries=BCs)

hardware = Hardware(backend, workgroup)

Struct used to configure the backend.

Inputs

• backend: used to specify the backend e.g. CPU(), CUDABackend() or other backends supported by KernelAbstraction.jl
• workgroup::Int this is an integer specifying the number of workers that cooperate in a parallel run. For GPUs this could be set to the size of the device's warp e.g. workgroup = 32. On CPUs, the default value in KernelAbstractions.jl is currently workgroup = 1024.

Output

This function returns a Hardware object with the fields backend and workgroup which are accessed by internally in XCALibre.jl to execute a given kernel in the target backend.

Forward Euler (1st order) explicit time stepping.

HLLC (Harten-Lax-van Leer Contact) flux scheme.

MUSCL{L}()

2nd-order MUSCL reconstruction with slope limiter L (VanLeer, MinMod, or Superbee). Set via reconstruction = MUSCL{VanLeer}() in the schemes NamedTuple.

MinMod limiter (most diffusive TVD): ψ(r) = max(0, min(1, r))

SSP-RK2 / Heun's method (2nd order) explicit time stepping. W^(1) = W^n - (dt/V)R(W^n) W^(2) = W^(1) - (dt/V)R(W^(1)) W^{n+1} = 0.5*(W^n + W^(2))

Multistage Runge-Kutta schemes of Jameson (low storage, optimised for steady-state convergence): W^(0) = W^n, W^(k) = W^n - α_k(dt/V)R(W^(k-1)), W^{n+1} = W^(N) with RK3() α = (0.1481, 0.4, 1), RK4() α = (1/4, 1/3, 1/2, 1) and RK5() α = (1/4, 1/6, 3/8, 1/2, 1). They extend the stability limit (CFL up to ≈ N-1) and are meant for steady problems together with LocalTimeStepping; they are not time accurate. Viscous fluxes and turbulence are evaluated with the fields at the start of the step in every stage.

Rusanov (Local Lax-Friedrichs) flux scheme.

Superbee limiter (least diffusive TVD): ψ(r) = max(0, min(2r, 1), min(r, 2))

Van Leer smooth limiter: ψ(r) = (r + |r|) / (1 + |r|)

blend_mixture_nu!(nu_field, alpha_field, rho_field, mu_0, mu_1)

Mixture kinematic viscosity built from dynamic viscosity field:

nu = (mu_0 * alpha + mu_1 * (1 - alpha)) / rho_blend

blend_properties!(property_field, alpha_field, property_0, property_1)

Linearly blends a per-phase scalar property using the volume-fraction field:

property_field = property_0 * alpha + property_1 * (1 - alpha)

cell_grad_magnitude!(mag_field, grad, config)

Cell-centred |∇alpha| from a gradient field, used as the |∇alpha|-weighted face interpolation weight for kappa_f so that near wall cells with small |∇alpha| don't ruin the surface-tension force.

cell_mueff(rho, nu, turbulence)

Cell μeff = ρ(ν + νt) for transpose_stress! in the compressible and multiphase solvers, the cell counterpart of their face mueff = rhof*(nuf + nutf). rho and nu must be cell values (constant or cell fields); nothing otherwise.

cell_nueff(nu, turbulence)

Cell νeff = ν + νt for transpose_stress! in the incompressible solvers (ν alone for laminar flow). nu must be a cell value (constant or cell field); nothing otherwise, and the term then uses the face viscosity at every face.

compression_flux!(phirf, ∇alphaf, mdotf, C_alpha, config)

Anti-diffusive compression face flux:

phir_f · Sf = min(Cα · |phi_f|/|Sf|, Phi_max) · (nnhatf · Sf)

Return fixed dt from runtime.dt[1] (no kernel launch).

compute_gh!(gh, g, config)

Computes g . x at cell centres. Used for hydrostatic pressure reconstruction.

compute_ghf!(ghf, g, config)

Computes g . x at face centres.

correct_velocity_rgh!(U, Hv, ∇p, rD, config)

Slightly modified version of correct_velocity! for convenience.

cpiso!(model, config;
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true)

Compressible and transient variant of the PISO algorithm with a sensible enthalpy transport equation for the energy.

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)
• outer_loops passes over momentum, energy and pressure per time step, so that the density of the time terms matches the state equation at the new time level (default = 1)
• transpose_stress include the explicit viscous stress ∇·(μ_eff dev2((∇U)ᵀ)) in the momentum equation (default = true)

Output

• Ux Vector of x-velocity residuals for each iteration.
• Uy Vector of y-velocity residuals for each iteration.
• Uz Vector of y-velocity residuals for each iteration.
• p Vector of pressure residuals for each iteration.

csimple!(
    model_in, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true,
    transpose_stress=true, consistent=false
)

Compressible variant of the SIMPLE algorithm with a sensible enthalpy transport equation for the energy.

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)
• transpose_stress include the explicit viscous stress ∇·(μ_eff dev2((∇U)ᵀ)) in the momentum equation (default = true)
• consistent use the consistent (SIMPLEC) pressure–velocity coupling (default = false). The pressure equation and the velocity and face-flux corrections then use rDt = V/(aP + Σ aN), the inverse row sum of the under-relaxed momentum matrix (aN ≤ 0 are its off-diagonal coefficients), instead of rD = V/aP. Unlike rD, rDt does not shrink with the velocity relaxation factor, so the pressure Laplacian keeps its weight against the compressibility term ∇·(ψf ((H/aP)·S) p) of Compressible fluids, which can keep steady transonic runs stable where the default coupling amplifies round-off until the run diverges. The converged solution is unchanged; only the path to it differs. Rows with ΣN aN + aP ≤ 10⁻³ aP keep rD

Output

• Ux Vector of x-velocity residuals for each iteration.
• Uy Vector of y-velocity residuals for each iteration.
• Uz Vector of y-velocity residuals for each iteration.
• p Vector of pressure residuals for each iteration.
• e Vector of energy residuals for each iteration.

godunov!(model, config; output=VTK())

Explicit Godunov-type density-based solver for compressible flows, including supersonic flows with shockwaves. Solves the Euler/Navier-Stokes equations in conservative form [ρ, ρU, ρE] using finite-volume flux reconstruction. Dispatched from run! when model.fluid isa SupersonicFlow.

Input arguments

• model reference to a Physics model defined by the user with fluid = Fluid{SupersonicFlow}(...).
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())

Scheme options (set via schemes in Configuration)

• flux — Riemann solver: Rusanov() (default) or HLLC()
• reconstruction — spatial reconstruction: Upwind() (1st order, default) or MUSCL{L}() (2nd order) where L is a slope limiter: VanLeer, MinMod, or Superbee
• time_stepping — explicit time integration: FEuler() (1st order, default) or RK2() (2nd order SSP) or the steady multistage schemes RK3(), RK4(), RK5()

Runtime(adaptive=LocalTimeStepping(maxCo=..., ramp=...)) enables per-cell pseudo-time stepping (steady problems).

forces=(patches=[:wall], interval=10) records pressure and viscous force histories in residuals.forces.

On a distributed mesh (see distribute) each rank advances its own cells; residuals, time steps and forces are reduced over all ranks. Write results with output=OpenFOAM() (3D meshes) or set write_interval=-1.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.rho) with the following entries:

• rho Vector of density residuals for each time step.

interpolate_weighted!(phif, phi, weight_field, config)

Weighted face interpolation: phif[f] = (phi_c1·w_c1 + phi_c2·w_c2) / (w_c1 + w_c2). With w = |∇alpha| it preserves kappa_f at interfaces and zeros it in bulk where kappa is noisy.

laplace!(model_in, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true)

Top-level entry point for solving the Laplace (heat conduction) equation on model.domain. Optionally runs in steady or transient mode and can call Conduction energy model for high fidelity (k and cp are recomputed at each iteration).

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• T Vector of T residuals for each iteration.

mules_limit!(mp_model::AbstractMultiphaseModel,
             phiAf, alpha_prev, phiLf,
             Pplus, Pminus, Qplus, Qminus, Rplus, Rminus,
             alphaMaxLocal, alphaMinLocal,
             dt, mesh, config)

Unified MULES (Zalesak FCT) flux-limiter for explicit alpha transport, shared by both VOF and Mixture sub models. Limits the anti-diffusive face flux phiAf = phiHf - phiLf so the explicit update α^{n+1} = α^n - dt/V · div(phiLf + phiAf) stays bounded.

Logic:

1. mules_set_bounds!: sets per-cell alphaMax, alphaMin Dispatched:
   • VOF local neighbour extrema of alpha_prev between [0,1]; prevents isolated alpha=1 cells; preserves sharp interface.
   • Mixture hard [0,1] limit; must tolerate uniform alpha without freezing the anti-diffusive drift flux.
2. Build low-order α*, P+-, Q+-
3. Ratios: R+- = clamp(Q+-/P+-, 0, 1).
4. Per-face apply lambda_f (scales phiAf).

Boundary faces are left unlimited (lambda = 1).

nhat_prep!(nhatf_prep, alpha, ∇alphaf, config)

Builds a face-normal unit vector field from the face-interpolated alpha grad.

phi_gf!(phi_gf, rho, ghf, rDf, model, config)

Builds the gravity (buoyancy) contribution to the face mass flux. On each face computes

phi_gf[f] = -ghf[f] · area · snGrad(rho) · rDf[f]

It is summed into mdotf before the pressure solve and reconstructed to a cell vector (phi_g) for the velocity correction. A face with no density jump (single-phase region) contributes zero.

cpiso!(model, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true)

Incompressible and transient variant of the SIMPLE algorithm to solving coupled momentum and mass conservation equations.

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)
• transpose_stress include the explicit viscous stress ∇·(ν_eff dev2((∇U)ᵀ)) in the momentum equation (default = true)

Output

• Ux Vector of x-velocity residuals for each iteration.
• Uy Vector of y-velocity residuals for each iteration.
• Uz Vector of y-velocity residuals for each iteration.
• p Vector of pressure residuals for each iteration.

potential_flow!(model, config; ncorrectors=0, pref=nothing, time=0, petsc_options="", from_rest=true)

Set the velocity to a divergence-free potential-flow field. With from_rest=true (default) the internal velocity is first set to zero, so the face flux to be corrected comes from the velocity boundary conditions only; with from_rest=false the current velocity field is projected. A nonzero initial field that is not a discrete gradient (a uniform field on a non-orthogonal mesh, for instance) keeps a solenoidal part after the projection, which shows up as spurious circulation, notably in enclosed or near-stagnant regions. Velocity-potential boundary conditions are inferred from pressure: fixed pressure becomes fixed zero potential, periodic patches remain periodic, and other patches use zero normal gradient.

The corrected face-volume flux is returned with the residual of the last potential solve, taken before that solve.

pressure_grad!(p_rgh, ∇p_rghf_deconstructed, phi_gf, rDf, config)

Builds the face-normal pressure-gradient field used to correct the cell velocity (Rhie-Chow consistent). On each face computes

∇p_rghf_deconstructed[f] = (phi_gf[f] - snGrad(p_rgh)·area·rDf[f]) / (rDf[f])

The result is a face scalar later reconstructed to a cell vector (reconstruct!) and fed to correct_velocity_rgh!.

prime_face_nut!(model, boundaries, config)

Set the face eddy viscosity from the current cell values at the start of every run, as turbulence! does: interpolated to the faces, boundary values from the conditions, and the wall-function values on wall faces. turbulence! updates the face values only after the first momentum solve, so without this the first iterations would see face values that do not match the cell νt used by the explicit transpose stress: zero after initialisation (the implicit Laplacian then holds ν alone while the explicit term carries the initial νt, typically 10³-10⁴ ν), or values left by an earlier run on fields since re-initialised. Either makes the explicit part dominate the first iterations, which can trigger a divergence.

reconstruct!(phi::VectorField, psif::FaceScalarField, moments, config)

Least-squares reconstruction of a cell-centred vector field phi from a face-normal scalar field psif. Both internal and boundary faces contribute to the least-squares system. moments is a preallocated (ncells, 9) scratch array. Required for discretisation consistency.

1st-order reconstruction — returns cell-centred values unchanged.

2nd-order MUSCL reconstruction with slope limiter L. Slope ratio (Darwish & Moukalled): r = 2*(∇φ·dLR)/Δ − 1. Velocity components reconstructed via gradU * dLR (matrix–vector product). ρ and p clamped to positive floor to prevent NaN in sqrt(γp/ρ).

function run!(
    model::Physics, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    )

    # here an internal function is used for solver dispatch
    return residuals
end

This is the top level API function to initiate a simulation. It uses the user-provided model defined as a Physics object to dispatch to the appropriate solver.

Dispatched flow solvers

• Steady incompressible (SIMPLE algorithm for coupling)
• Transient incompressible (PISO algorithm for coupling)
• Steady weakly compressible (SIMPLE algorithm for coupling)
• Transient weakly compressible (PISO algorithm for coupling)

On a distributed mesh only the incompressible SIMPLE and PISO solvers, the Laplace solver and the density-based solver (SupersonicFlow) are supported; any other solver or model errors rather than solving each rank's block on its own.

Calling run! again on the same model continues from its current fields, including the face mass flux of the last iteration, and retraces the uninterrupted run. Re-initialising the fields (e.g. with initialise!) between calls starts afresh.

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)
• progress show a progress bar with the iteration and residuals (default = true). Use progress=false for large-scale runs that are not on a local PC, e.g. cluster batch jobs. Distributed runs never show it
• transpose_stress include the explicit part of the viscous stress, ∇·(μ_eff dev2((∇U)ᵀ)), in the momentum equation of the incompressible, compressible and multiphase solvers (default = true). It is zero for a constant viscosity in incompressible flow and needed wherever the effective viscosity varies, e.g. with a turbulence model; false drops it
• consistent steady solvers (simple! and csimple!) only: use the consistent (SIMPLEC) pressure–velocity coupling, which needs little or no pressure under-relaxation and keeps steady transonic runs stable for common relaxation factors; see simple! and csimple! (default = false)
• restart resume an incompressible run on a distributed mesh from results it wrote with output=OpenFOAM(): the iteration (steady) or time (transient) of a written time directory, or that directory's name. The resumed run continues from that point up to iterations and retraces the uninterrupted run (default = nothing)

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux). The fields available within the returned residuals tuple depend on the solver used. For example, for an incompressible solver, a x-momentum equation residual can be retrieved accessing the Ux field i.e. residuals.Ux. Look at reference guide for each dispatch method to find out which fields are available.

Example

residuals = run!(model, config) 

# to access the pressure residual

residuals.p 

run!(
    model::Physics{T,F,M,Tu,E,D,BI}, config;
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T,F<:SupersonicFlow,M,Tu,E,D,BI}

Calls the explicit density-based solver using the Rusanov (Local Lax-Friedrichs) flux for supersonic/hypersonic compressible flows.

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user.
• output select the format used for simulation results (default = VTK())

Output

Returns a NamedTuple with the field:

• rho - Vector of density L2 residuals for each iteration.

run!(
    model::Physics{T,F,M,Tu,E,D,BI}, config;
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T<:Steady,F<:Incompressible,M,Tu,E,D,BI} = 
begin
    residuals = simple!(model, config, pref=pref)
    return residuals
end

Calls the incompressible steady solver using the SIMPLE algorithm.

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux - Vector of x-velocity residuals for each iteration.
• Uy - Vector of y-velocity residuals for each iteration.
• Uz - Vector of y-velocity residuals for each iteration.
• p - Vector of pressure residuals for each iteration.

run!(
    model::Physics{T,F,M,Tu,E,D,BI}, config;
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T<:Steady,F<:Incompressible,M,Tu,E,D,BI} = 
begin
    residuals = simple!(model, config, pref=pref)
    return residuals
end

Calls the incompressible steady solver using the SIMPLE-MRF algorithm.

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux - Vector of x-velocity residuals for each iteration.
• Uy - Vector of y-velocity residuals for each iteration.
• Uz - Vector of y-velocity residuals for each iteration.
• p - Vector of pressure residuals for each iteration.

run!(
    model::Physics{T,F,M,Tu,E,D,BI}, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T<:Steady,F<:WeaklyCompressible,M,Tu,E,D,BI} = 
begin
    residuals = csimple!(model, config, pref=pref); #, pref=0.0)
    return residuals
end

Calls the compressible steady solver using the SIMPLE algorithm for weakly compressible fluids.

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux - Vector of x-velocity residuals for each iteration.
• Uy - Vector of y-velocity residuals for each iteration.
• Uz - Vector of y-velocity residuals for each iteration.
• p - Vector of pressure residuals for each iteration.
• e - Vector of energy residuals for each iteration.

run!(
    model::Physics{T,F,M,Tu,E,D,BI}; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T<:Transient,F<:WeaklyCompressible,M,Tu,E,D,BI} = 
begin
    residuals = cpiso!(model, config)
    return residuals
end

Calls the compressible transient solver using the PISO algorithm for weakly compressible fluids.

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux - Vector of x-velocity residuals for each iteration.
• Uy - Vector of y-velocity residuals for each iteration.
• Uz - Vector of y-velocity residuals for each iteration.
• p - Vector of pressure residuals for each iteration.
• e - Vector of energy residuals for each iteration.

run!(
    model::Physics{T,F,SO,M,Tu,E,D,BI}, config;
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T,F,SO<:Uniform,M,Tu,E,D,BI} = 
begin
    residuals = laplace!(model, config, pref=pref)
    return residuals
end

Top-level entry point for solving the Laplace (heat conduction) equation on model.domain. Optionally runs in steady or transient mode and can call CryogenicConduction energy model for high fidelity (k and cp are recomputed at each iteration).

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• T - (Vector?) of temperature residuals for each iteration.

run!(
    model::Physics{T,F,M,Tu,E,D,BI}, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true
    ) where{T<:Transient,F<:Incompressible,M,Tu,E,D,BI} = 
begin
    residuals = piso!(model, config, pref=pref); #, pref=0.0)
    return residuals
end

Calls the incompressible transient solver using the PISO algorithm.

Input

• model represents the Physics model defined by user.
• config Configuration structure defined by user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM() (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only.

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux - Vector of x-velocity residuals for each iteration.
• Uy - Vector of y-velocity residuals for each iteration.
• Uz - Vector of y-velocity residuals for each iteration.
• p - Vector of pressure residuals for each iteration.

simple!(model_in, config;
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true,
    transpose_stress=true, consistent=false)

Incompressible variant of the SIMPLE algorithm to solving coupled momentum and mass conservation equations.

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)
• transpose_stress include the explicit viscous stress ∇·(ν_eff dev2((∇U)ᵀ)) in the momentum equation (default = true)
• consistent use the consistent (SIMPLEC) pressure–velocity coupling (default = false). The pressure equation and the velocity and face-flux corrections then use rDt = V/(aP + Σ aN), the inverse row sum of the under-relaxed momentum matrix, instead of rD = V/aP, with the terms (rDt - rD)∇p added to the predicted velocity and face fluxes. rDt does not shrink with the velocity relaxation factor, so the pressure correction is not over-damped and pressure needs little or no under-relaxation; the converged solution is unchanged. Rows with ΣN aN + aP ≤ 10⁻³ a_P keep rD

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux Vector of x-velocity residuals for each iteration.
• Uy Vector of y-velocity residuals for each iteration.
• Uz Vector of y-velocity residuals for each iteration.
• p Vector of pressure residuals for each iteration.

simple!(model_in, config; 
    output=VTK(), pref=nothing, ncorrectors=0, inner_loops=0, progress=true)

Incompressible variant of the SIMPLE algorithm to solving coupled momentum and mass conservation equations.

Input arguments

• model reference to a Physics model defined by the user.
• config Configuration structure defined by the user with solvers, schemes, runtime and hardware structures configuration details.
• output select the format used for simulation results from VTK() or OpenFOAM (default = VTK())
• pref Reference pressure value for cases that do not have a pressure defining BC. Incompressible solvers only (default = nothing)
• ncorrectors number of non-orthogonality correction loops (default = 0)
• inner_loops number to inner loops used in transient solver based on PISO algorithm (default = 0)
• transpose_stress include the explicit viscous stress ∇·(ν_eff dev2((∇U)ᵀ)) in the momentum equation (default = true)

Output

This function returns a NamedTuple for accessing the residuals (e.g. residuals.Ux) with the following entries:

• Ux Vector of x-velocity residuals for each iteration.
• Uy Vector of y-velocity residuals for each iteration.
• Uz Vector of y-velocity residuals for each iteration.
• p Vector of pressure residuals for each iteration.

Forward-Euler step: W^{n+1} = W^n - (dt/V)*R(W^n). Call recover_primitives! after.

SSP-RK2 step: W^(1)=W^n-(dt/V)R(W^n), W^(2)=W^(1)-(dt/V)R(W^(1)), W^{n+1}=0.5*(W^n+W^(2)). R(W^n) must be pre-computed. Call recover_primitives! after.

surface_tension_flux!(rDf, sigma, kappaf, alpha, phi_gf, config)

CSF surface-tension contribution to the face flux: phi_gf -= σ · κf · (∇αf · Sf) · rDf.

transpose_stress!(source, mueff, gradU, U_BCs, time, config; cell_mueff=nothing)

Set source = ∇·(μeff dev2((∇U)ᵀ)) from the cell gradient gradU and the face viscosity mueff (νeff for the incompressible solvers), using the current values of both. With the cell viscosity cell_mueff (see cell_nueff and cell_mueff) internal and periodic faces interpolate the product μeff dev2((∇U)ᵀ) of the two cells, and boundary faces use mueff times the boundary-face gradient (see `transposestressboundary!), whose fixed-value conditions, includingDirichletFunction, take their face value attime`.

well_balanced_pressure_grad!(grad_field, face_buf, p_rgh, rho, ghf,
                              mesh, moments, config; sigma=0, kappaf=nothing,
                              alpha=nothing)

Overrides grad_field.values with a face-snGrad reconstruction of the predictor body force terms:

face_buf[f] = area_f · ( snGrad(p_rgh) + ghf · snGrad(rho) - sigma·kappaf · snGrad(alpha) )

*Important for stability.

FieldAverage(
#required arguments
field;
name::String,

#optional keyword arguments
start::Union{Real,Nothing},
stop::Union{Real,Nothing},
update_interval::Union{Real,Nothing})

Constructor to allocate memory to store the time averaged field. Once created, should be passed to the Configuration object as an argument with keyword postprocess

Input arguments

• field the VectorField or ScalarField to be averaged, e.g , model.momentum.U.
• name::String the name of the field to be averaged, e.g "U_mean", this is used only when exporting to .vtk format

Optional arguments

• start::Union{Real,Nothing} optional keyword which specifies the start time/iteration of the averaging window, for steady simulations, this is in iterations, for transient simulations it is in flow time.
• stop::Union{Real,Nothing} optional keyword which specifies the end iteration/time of the averaging window. Default value is the last iteration/timestep.
• update_interval::Union{Real,Nothing} optional keyword which specifies how often the time average of the field is updated and stored (default value is 1 i.e average updates every timestep/iteration). Note that the frequency of writing the post-processed fields is specified by the write_interval in Configuration.

FieldRMS(
#required arguments
field;
name::String,

#optional keyword arguments
start::Union{Real,Nothing},
stop::Union{Real,Nothing},
update_interval::Union{Real,Nothing})

Constructor to allocate memory to store the root mean square of the fluctuations of a field over the averaging window. Once created, should be passed to the Configuration object as an argument with keyword postprocess

Input arguments

• field the VectorField or ScalarField, e.g , model.momentum.U.
• name::String the name/label of the field, e.g "U_rms", this is used only when exporting to .vtk format

Optional arguments

• start::Union{Real,Nothing} optional keyword which specifies the start of the RMS calculation window, for steady simulations, this is in iterations, for transient simulations it is in flow time.
• stop::Union{Real,Nothing} optional keyword which specifies the end iteration/time of the RMS calculation window. Default value is the last iteration/timestep.
• update_interval::Union{Real,Nothing} optional keyword which specifies how often the RMS of the field is updated and stored (default value is 1 i.e RMS updates every timestep/iteration). Note that the frequency of writing the post-processed fields is specified by the write_interval in Configuration.

Probe(
#required arguments
field,
mesh_cpu; 

#required keyword arguments
location::AbstractVector
name::AbstractString

#optional keyword arguments
output::ProbeOutput 
start::Union{Real,Nothing},
stop::Union{Real,Nothing},
update_interval::Union{Real,Nothing})

Constructor to probe an internal field during runtime and write it out to a .csv file. Once created, should be passed to the Configuration object as an argument with keyword postprocess

Input arguments

• field, the field you want to probe, such as model.momentum.U
• mesh_cpu, the mesh object stored on the CPU needs to be passed, this is to find the nearest cell centre to the location you want to probe. On a distributed mesh pass the distributed mesh (on the CPU) on every rank: the rank owning the nearest cell samples it and rank 0 writes the file

Input Keyword arguments

• location, e.g location = [1.0, 2.0, 0.0] and must be a vector ∈ ℝ³ even for 2D simulations
• name, the name of the output file e.g "Velocity probe at [1.0, 2.0, 3.0]"

Optional arguments

• output::ProbeOutput the format to output the probe, defaults to CSV, can only pass either TXT() or CSV()
• start::Union{Real,Nothing} optional keyword which specifies the start of the probe output window, for steady simulations, this is in iterations, for transient simulations it is in flow time.
• stop::Union{Real,Nothing} optional keyword which specifies the end iteration/time of the probe output window. Default value is the last iteration/timestep.
• update_interval::Union{Real,Nothing} optional keyword which specifies how often the probe's field is output to file (default value is 1 i.e probe writes field value every iteration/timestep).

Qcriterion(
#required arguments
model.momentum.U; 

#optional keyword arguments
start::Union{Real,Nothing},
stop::Union{Real,Nothing},
update_interval::Union{Real,Nothing})

Constructor to allocate memory to store the Q-criterion over the calculation window. Once created, should be passed to the Configuration object as an argument with keyword postprocess

Input arguments

• field, must be model.momentum.U

Optional arguments

• start::Union{Real,Nothing} optional keyword which specifies the start of the Q-criterion calculation window, for steady simulations, this is in iterations, for transient simulations it is in flow time.
• stop::Union{Real,Nothing} optional keyword which specifies the end iteration/time of the Q-criterion calculation window. Default value is the last iteration/timestep.
• update_interval::Union{Real,Nothing} optional keyword which specifies how often the Q-criterion is updated and stored (default value is 1 i.e Q-criterion updates every timestep/iteration). Note that the frequency of writing the post-processed fields is specified by the write_interval in Configuration.

Q is only computed inside the window. Outside it the stored field keeps its last computed value (zero before start), and that value is what the post-processed output files contain.

ReynoldsStress(
#required arguments
model.momentum.U; 

#optional keyword arguments
start::Union{Real,Nothing},
stop::Union{Real,Nothing},
update_interval::Union{Real,Nothing})

Constructor to allocate memory to store the Reynolds Stress Tensor over the calculation window. Once created, should be passed to the Configuration object as an argument with keyword postprocess

Input arguments

• field, must be model.momentum.U

Optional arguments

• start::Union{Real,Nothing} optional keyword which specifies the start of the Reynolds Stress Tensor calculation window, for steady simulations, this is in iterations, for transient simulations it is in flow time.
• stop::Union{Real,Nothing} optional keyword which specifies the end iteration/time of the Reynolds Stress Tensor calculation window. Default value is the last iteration/timestep.
• update_interval::Union{Real,Nothing} optional keyword which specifies how often the Reynolds Stress Tensor is updated and stored (default value is 1 i.e Reynolds Stress Tensor updates every timestep/iteration). Note that the frequency of writing the post-processed fields is specified by the write_interval in Configuration.

function boundary_average(patch::Symbol, field, config; time=0)
    # Extract mesh object
    mesh = field.mesh

    # Determine ID (index) of the boundary patch 
    ID = boundary_index(mesh.boundaries, patch)
    @info "calculating average on patch: $patch at index $ID"
    boundary = mesh.boundaries[ID]
    (; IDs_range) = boundary

    # Create face field of same type provided by user (scalar or vector)
    sum = nothing
    if typeof(field) <: VectorField 
        faceField = FaceVectorField(mesh)
        sum = zeros(_get_float(mesh), 3) # create zero vector
    else
        faceField = FaceScalarField(mesh)
        sum = zero(_get_float(mesh)) # create zero
    end

    # Interpolate CFD results to boundary
    interpolate!(faceField, field, config)
    correct_boundaries!(faceField, field, field.BCs, time, config)

    # Calculate the average
    for fID ∈ IDs_range
        sum += faceField[fID]
    end
    ave = sum/length(IDs_range)

    # return average
    return ave
end

pressure_force(patch::Symbol, model, config)

Calculate the pressure force vector [Fx, Fy, Fz] acting on a given patch/boundary.

Input arguments

• patch::Symbol name of the boundary of interest (as a Symbol)
• model Physics object defining the simulation
• config Configuration object (provides the hardware backend)

On a distributed mesh the call is collective: every rank must make it, and all ranks return the force on the whole patch.

stress_tensor(U::VectorField, ν, νt, config; time=0)

Function to calculate the stress tensor.

Input arguments

• U::VectorField velocity field
• ν laminar viscosity of the fluid
• νt eddy viscosity from turbulence models. Pass ConstantScalar(0) for laminar flows
• config need to pass Configuration object as this contains the boundary conditions

Optional arguments

• time the time at which the velocity boundary conditions are evaluated (default 0), as in viscous_force

viscous_force(patch::Symbol, model, config; time=0)

Calculate the viscous force vector [Fx, Fy, Fz] acting on a given patch/boundary.

Input arguments

• patch::Symbol name of the boundary of interest (as a Symbol)
• model Physics object defining the simulation
• config Configuration object (provides boundary conditions and hardware backend)

Optional arguments

• time the time at which the velocity condition of the patch is evaluated (default 0). It matters only for a DirichletFunction whose function depends on time: pass the solution time (the iteration number for steady runs), as the solvers do when they evaluate the condition.

The force of each face is the viscous stress the momentum equation applies there, -ρ ν_eff (∇U_b + ∇U_bᵀ - ⅔ tr(∇U_b) I)⋅n A. ν_eff = ν + ν_t uses the wall-function value of ν_t on the face, and ∇U_b is the owner-cell gradient (the gradient scheme of U in config.schemes, Gauss if none) with its normal component replaced by the normal gradient of the velocity condition, (U_b - U_c)/δ. This includes the transposed (explicit) stress term of the momentum equation, so the result is the force the solver applies. U_b is the face value of the velocity boundary condition: the fixed vector of Wall and Dirichlet, the function of DirichletFunction evaluated at the face centre and time, or the rotation velocity of RotatingWall. Zero-gradient conditions have a zero normal gradient, so only the tangential derivatives of the cell gradient act there. Patches with a Slip or Symmetry velocity boundary condition exert no shear, so their viscous force is zero.

On a distributed mesh the call is collective: every rank must make it, and all ranks return the force on the whole patch.

wall_shear_stress(patch::Symbol, model, config; time=0)

Function to calculate the wall shear stress acting on a given patch/boundary. Returns the kinematic viscous stress acting on each face of the patch, -ν_eff (∇U_b + ∇U_bᵀ - ⅔ tr(∇U_b) I)⋅n (the face values of viscous_force per unit area and density; zero on Slip and Symmetry patches), as a FaceVectorField of the patch length, and the face centres as a vector.

Input arguments

• patch::Symbol name of the boundary of interest (as a Symbol)
• model instance of Physics object needs to be passed
• config need to pass Configuration object as this contains the boundary conditions

Optional arguments

• time the time at which the velocity condition of the patch is evaluated (default 0), as in viscous_force, which uses the same face values.

setField_Box!(; mesh, field, value::F, min_corner::V, max_corner::V, config)

Sets field values to value for all cells whose centre lies within the axis-aligned box defined by min_corner and max_corner. Runs on the backend and workgroup set in config.hardware.

Warning: if a cell's outer boundary extends outside the box but its centre lies within it, that cell is still counted.

Returns the number of cells set.

setField_Circle2D!(; mesh, field, value::F, centre::V, radius::F, config)

Sets field values to value for cells whose centre is within radius of centre (given as [x, y], assumed to lie in the X-Y plane at z=0). The comparison is made against each cell's full 3-D centre, so a cell whose centre has a non-zero z-offset is measured with that offset intact — this is only exact for meshes lying in the X-Y plane. Runs on the backend and workgroup set in config.hardware.

Returns the number of cells set.

setField_Expression!(; mesh, field, condition, value_true::F, value_false=nothing, config)

Assigns value_true to every cell whose centre (x, y, z) satisfies condition(x, y, z). If value_false is provided all other cells receive value_false; otherwise they are unchanged. Runs on the backend and workgroup set in config.hardware (CPU or GPU). On GPU backends condition is compiled into the kernel, so it must not allocate, and it must not read non-const globals or capture arrays. Capture parameters with let (as below) or use const globals; a named function needs no @inline.

Returns the number of cells where condition returned true.

Example — Rayleigh-Taylor interface

interface = let Ly = 4.0, A = 0.05, λ = 1.0
    (x, y, z) -> y > Ly/2 + A * cos(2π * x / λ)
end
setField_Expression!(
    mesh      = mesh,
    field     = model.fluid.alpha,
    condition = interface,
    value_true  = 1.0,
    value_false = 0.0,
    config      = config,
)

setField_Sphere3D!(; mesh, field, value::F, centre::V, radius::F, config)

Sets field values to value for cells whose centre lies within radius of centre. Runs on the backend and workgroup set in config.hardware.

Returns the number of cells set.

DistributedMesh <: AbstractMesh

Wraps a rank-local mesh (owned + one ghost layer) with partition and communication metadata. All non-metadata properties forward to the wrapped mesh, so fields, Physics and kernels treat it as a normal mesh.

FOAMCase(dir; scale=1, integer_type=Int32, float_type=Float64)

An OpenFOAM case already decomposed into processor<rank>/ directories (by decomposePar or by a distributed XCALibre run with output=OpenFOAM()). distribute(FOAMCase(dir)) gives each rank its own part without any rank reading the global mesh; the keywords are those of FOAM3D_mesh.

bind_device!(backend; comm=MPI.COMM_WORLD)
bind_device!(backend, rank)

Bind this MPI rank to a GPU. Without rank, the rank's position among the ranks on its own node picks the device, so ranks on one node take its devices in turn whatever the global rank order; when a node has more ranks than devices they share, and a warning says so. bind_device!(backend, rank) binds to device rank % ndevices explicitly. No-op on CPU. Call before adapt(backend, dmesh) or building fields/HaloExchange on a GPU backend.

distribute(dir::AbstractString; comm=MPI.COMM_WORLD)

Load an offline decomposition written by partition_mesh: each rank reads only its own rank_<rank>.xdm from dir (no rank-0 memory bottleneck).

distribute(mesh; comm=MPI.COMM_WORLD, periodic_patches=())

Online mesh distribution: rank 0 partitions mesh (Metis k-way) and scatters one DistributedMesh per rank; other ranks may pass nothing as mesh.

periodic_patches takes patch-name pairs, e.g. [(:top, :bottom)], for meshes with periodic boundaries: matched owner cells are contracted in the partition graph so each periodic pair lands on one rank, and construct_periodic on the DistributedMesh then works per rank exactly as in serial.

A constant/polyMesh in the working directory holding mesh numbered differently (e.g. after reorder_mesh!) is rewritten in the order of mesh, so decomposed output maps onto it.

distribute(reader::Function; dir=nothing, key=nothing, comm=MPI.COMM_WORLD, periodic_patches=())

Read and distribute a mesh from a call every rank makes identically. Rank 0 runs reader() to build the global mesh and the other ranks skip it, so no rank == 0 guard appears in the script and no value can differ in type between ranks:

mesh = distribute() do
    UNV3D_mesh(path, scale=0.001)
end

With dir, the decomposition is written there once and every rank then loads only its own part, which removes rank 0's global-mesh memory ceiling from later runs. A decomposition already in dir for the same number of ranks and the same key is reused and reader() is never called; otherwise it is replaced. Without a key only the rank count is compared, so pass one that names the mesh (any hashable value, e.g. the reader's arguments) whenever dir may hold another mesh's parts:

mesh = distribute(dir="parts", key=(path, 0.001)) do
    UNV3D_mesh(path, scale=0.001)
end

MPI is initialised if it is not already.

distribute(case::FOAMCase; comm=MPI.COMM_WORLD)

Load a decomposed OpenFOAM case: each rank reads processor<rank>/constant/polyMesh, and ghost cells are built by one exchange with each neighbour across the processor patches. Run under mpiexec -n <number of processor directories>. Original cell ids come from cellProcAddressing when the case has it, so gather returns fields in the undecomposed order; otherwise in rank-block order.

gather(field, dmesh; comm=dmesh.comm, root=0)

Gather a distributed field to root in ORIGINAL (pre-partition) cell order. Returns a Vector (ScalarField) or (x, y, z) NamedTuple of Vectors (VectorField) on root, nothing on other ranks.

is_root()
is_root(comm)

True on the rank that owns terminal output and result writing, and true in a serial run where MPI was never initialised, so one guard serves both:

is_root() && println("final residual ", residuals.p[end])

mesh_info(path)

Header of a binary mesh file written by partition_mesh, as a NamedTuple: kind (:serial or :partitioned), nranks and rank (1 and 0 for a serial mesh), mesh (Mesh2 or Mesh3), the integer and float types TI and TF, n_owned and n_ghost cells, the format number and the XCALibre and Julia versions that wrote it, and the length of every stored array.

partition_mesh(mesh, nparts; dir, periodic_patches=(), key=nothing)

Offline decomposition: partition mesh into nparts rank-local meshes and write one rank_<r>.xdm per rank into dir, recording key for distribute(reader; dir, key). Load with distribute(dir; comm) under mpiexec -n nparts. Each file is binary with a header (mesh_info) naming its format, kind and rank count; a part of another format or rank count is refused at load with the call that fixes it. A constant/polyMesh in the working directory is brought to the order of mesh as in distribute.

repartition(dm::DistributedMesh; method=:ptscotch, petsc_options="")

Rebalance a distributed mesh in parallel: PETSc's MatPartitioning (method :ptscotch or :parmetis; PT-Scotch balances cell counts more tightly) partitions the distributed cell graph and cells migrate to their new ranks, so no rank ever holds the global mesh. Use it after distribute of a decomposition that is poorly balanced, such as decomposePar -method simple. Needs using PETSc with a PETSc built with the chosen package. PETSc is initialised here if it is not already, so start-up entries of petsc_options must be given here rather than to the solver.