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
    weight::F       # linear interpolation weight
end

struct Mesh3{VC, VI, VF<:AbstractArray{<:Face3D}, VB, VN, SV3, UR} <: AbstractMesh
    cells::VC           # vector of cells
    cell_nodes::VI      # vector of indices to access cell nodes
    cell_faces::VI      # vector of indices to access cell faces
    cell_neighbours::VI # vector of indices to access cell neighbours
    cell_nsign::VI      # vector of indices to with face normal correction (1 or -1 )
    faces::VF           # vector of faces
    face_nodes::VI      # vector of indices to access face nodes
    boundaries::VB      # vector of boundaries
    nodes::VN           # vector of nodes
    node_cells::VI      # vector of indices to access node cells
    get_float::SV3      # store mesh float type
    get_int::UR         # store mesh integer type
    boundary_cellsID::VI # vector of indices of boundary cell IDs
end

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

UNV2D_mesh(meshFile; scale=1, integer_type=Int64, 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 - select interger type to use in the mesh (Int32 may be useful on GPU runs)

• 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=Int64, 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=Int64: The integer type used for topological indices and connectivity arrays.
• 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=Int64, 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 - select interger type to use in the mesh (Int32 may be useful on GPU runs)

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

activate_multithread(backend::CPU; nthreads=1) = BLAS.set_num_threads(nthreads)

Convenience function to set number of BLAS threads.

Input arguments

• backend is the only required input which must be CPU() from KernelAbstractions.jl
• nthreads can be used to set the number of BLAS cores (default nthreads=1)

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<:AbstractMesh,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)

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)

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

Neumann <: AbstractNeumann

Neumann boundary condition model to set the gradient at the boundary explicitly (currently only configured for zero gradient)

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

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 <: AbstractDirichlet

Slip 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

Slip(:plate) # slip wall condition

Symmetry <: AbstractBoundary

Symmetry boundary condition vector and scalar fields. Notice that for scalar fields, this boundary condition applies an explicit zero gradient condition. In some rare cases, the use of an Extrapolated condition for scalars may be beneficial (to assign a semi-implicit zero gradient condition)

Input

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

Example

Symmetry(:freestream)

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

macro define_boundary(boundary, operator, definition)
    quote
        @inline (bc::$boundary)(
            term::Operator{F,P,I,$operator}, cellID, zcellID, cell, face, 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
• cell gives access to boundary cell object and corresponding information
• face gives access to boundary face object and corresponding information
• fID ID of the boundary face (to index Mesh2.faces vector)
• 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) = face  # extract 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$

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.

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

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 in BOTH mode=Cg() and mode=AMGSolver(); 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 in BOTH mode=Cg() and mode=AMGSolver(). 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. GATED to mode=AMGSolver(): 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.

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 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(_get_float(region)))^0.9,
        rtol=_get_float(region)(1e-1)

    ) 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: solver object from Krylov.jl and it could be one of Bicgstab(), Cg(), Gmres() which are re-exported in XCALibre.jl
• preconditioner: instance of preconditioner to be used e.g. Jacobi()
• convergence sets the stopping criteria of this field
• 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, or 200 for AMG)
• atol: absolute tolerance for the solver (default to eps(FloatType)^0.9)
• rtol: set relative tolerance for the solver (defaults to 1e-1)
• 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)

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.
• kf – Turbulent kinetic energy FaceScalarField.
• omegaf – Specific dissipation rate FaceScalarField.
• 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 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.

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.

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

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!(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))

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.

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)

Compute CFL-limited dt via per-cell kernel; update runtime.dt in place.

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)

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)

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
)

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)

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)

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)

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)

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)

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.

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

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

Least-squares reconstruction of a cell-centred vector field phi from a face-normal scalar field psif. 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
    )

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

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

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)

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)

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)

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.

well_balanced_pressure_grad!(grad_field, face_buf, p_rgh, rho, ghf,
                              mesh, 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

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.

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)

stress_tensor(U::VectorField, ν, νt, config)

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

viscous_force(patch::Symbol, model, config)

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)

wall_shear_stress(patch::Symbol, model,config)

Function to calculate the wall shear stress acting on a given patch/boundary.

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

setField_Box!(; mesh, field, value::F, min_corner::V, max_corner::V) where {F <: AbstractFloat, V <: AbstractVector}

Sets field values equal to value argument only for cells located within a box defined by min and max corners

Input arguments

• mesh reference to domain inside a Physics model defined by the user.
• field field to be modified.
• value new values for the chosen field.
• min_corner minimum coordinate values of the box in the format of [x, y, z].
• max_corner maximum coordinate values of the box in the format of [x, y, z].

setField_Circle2D!(mesh, field, value::F, centre::V, radius::F) where {F <: AbstractFloat, V <: AbstractVector}

Sets field values equal to value argument only for cells located within a 2D circle defined by its centre and radius

Input arguments

• mesh reference to domain inside a Physics model defined by the user.
• field field to be modified.
• value new values for the chosen field.
• centre coordinate of centre of the circle in the format of [x, y].
• radius length of the radius of the circle.

setField_Expression!(; mesh, field, condition::Function, value_true::F,
                      value_false::Union{F,Nothing}=nothing
                    ) where {F <: AbstractFloat}

Assigns value_true to every cell whose centre (x, y, z) satisfies the user-supplied condition. If value_false is provided, all other cells get value_false; otherwise they are left unchanged.

The condition can capture any variables from the calling scope via Julia's closure syntax - useful for parametric initial conditions like a sinusoidal interface for Rayleigh-Taylor or any expression that doesn't fit a box/sphere/circle.

Input arguments

• mesh reference to domain inside a Physics model defined by the user.
• field field to be modified.
• condition function (x, y, z) -> Bool evaluated at each cell centre.
• value_true value assigned to cells where the predicate returns true.
• value_false (optional) value assigned where the predicate returns false. If omitted, those cells are left unchanged.

Example — Rayleigh-Taylor interface

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

setField_Sphere3D!(mesh, field, value::F, centre::V, radius::F) where {F <: AbstractFloat, V <: AbstractVector}

Sets field values equal to value argument only for cells located within a 3D sphere defined by its centre and radius

Input arguments

• mesh reference to domain inside a Physics model defined by the user.
• field field to be modified.
• value new values for the chosen field.
• centre coordinate of centre of the sphere in the format of [x, y, z].
• radius length of the radius of the sphere.