Numerical setup

Available discretisation schemes, linear solvers and preconditioners

Discretisation schemes

Part of the methodology to solve the various model equations using the Finite Volume Method (FVM) is to discretise each equation term, essentially, this process linearises the partial differential equation so that it can be represented as a system of linear equations, which can be solved using linear algebra along with iterative solvers. This section presents the discretisation schemes currently available in XCALibre.jl.

Discretisation schemes in XCALibre.jl are organised under the abstract type AbstractScheme. As shown previously, a list of available schemes can be found using the subtypes function:

AbstractScheme
├─ BoundedUpwind
├─ CrankNicolson
├─ Euler
├─ Gauss
├─ LUST
├─ Linear
├─ Midpoint
├─ SteadyState
└─ Upwind

Schemes available

Time schemes

Laplacian schemes

Divergence schemes

Gradient schemes

Gradient limiters

The construction of all gradient limiters currently requires users to pass the mesh object that will be used for the simulations e.g. MFaceBased(mesh). Notice that the mesh object must be moved to the same device where the simulation will be performed. To prevent mistakes and avoid errors, it is recommended to construct the limiters using the simulation Physics object directly since it provide a reference to the mesh (domain) e.g. if a Physics object named model has been created, construct the limiter as MFaceBased(model.domain)

Specifying schemes

XCALibre.jl flow solvers offer considerable flexibility to users for defining discretisation schemes. However, this means that discretisation schemes must be specified for every term of every equation solved. The schemes must be provided as a NamedTuple where each keyword corresponds to the fields being solved, e.g. (U = ..., p = ..., k = ..., <field> = ...). To facilitate this process, the Schemes function is provided. Used without any inputs Schemes uses the default values provided (see details below).

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

For example, below we set the schemes for the U and p fields. Notice that in the first case the schemes will take their default values (entry for p). In the case of U, we are only changing the setting for the divergence scheme to Upwind.

using XCALibre

# Note: this example assumes a Physics object named `model` already exists

schemes = (
    p = Schemes(), # no input provided (will use defaults)
    U = Schemes(divergence = Upwind, gradient=Midpoint, limiter=MFaceBased(model.domain)),
)

Linear solvers

Linear solvers in XCALibre.jl are provided by Krylov.jl. The following solvers types are re-exported in XCALibre.jl

• Bicgstab() is a general purpose linear solver. Works well with non-symmetric matrices e.g. for U.
• Cg() is particular strong with symmetric matrices e.g to solve the pressure equation.
• Gmres() is a general solver. We have found it works best on the CPU backend.

For more information on these solvers you can review the excellent documentation provided by the Krylov.jl team.

XCALibre.jl provides the SolverSetup convenience function for setting solvers. See details below.

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

Preconditioners

XCALibre.jl offers a range of preconditioners which are subtypes of the abstract type PreconditionerType, exploring its subtypes we can find a list of the currently available preconditioners:

PreconditionerType
├─ LDIVPreconditioner
│  └─ DILU
└─ MULPreconditioner
   ├─ IC0GPU
   ├─ ILU0GPU
   ├─ Jacobi
   └─ NormDiagonal

Below an example is provided in context. Here, we are setting solvers for both the velocity field U and the pressure field p and packing them into a NamedTuple "solvers". The Jacobi preconditioner is used in both solvers. Notice that preconditioners are specified with an instance of their type i.e. Jacobi(). Internally, the preconditioner instance is used for dispatch. This tupple will then be passed on to create the final Configuration object.

using XCALibre

# Note: this example assumes a Physics object named `model` already exists

solvers = (
    U = SolverSetup(
        solver      = Bicgstab(), # Gmres()
        preconditioner = Jacobi(),
        convergence = 1e-7,
        relax       = 0.7,
        rtol = 1e-4,
        atol = 1e-10
    ),
    p = SolverSetup(
        solver      = Cg(), # Bicgstab(), Gmres()
        preconditioner = Jacobi(),
        convergence = 1e-7,
        relax       = 0.7,
        rtol = 1e-4,
        atol = 1e-10
    )
)