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

SchemeDescription
SteadyStatesets the time derivative to zero
EulerFirst order implicit Euler scheme
CrankNicolsonSecond order central scheme (not implemented yet)

Laplacian schemes

SchemeDescription
Linear2nd order Gauss gradient scheme with linear interpolation

Divergence schemes

SchemeDescription
Linear2nd order central difference
Upwind1st order upwind scheme
BoundedUpwindBounded version of the Upwind scheme
LUST1st/2nd order mixed scheme (fixed at 75% Linear - 25% Upwind)

Gradient schemes

SchemeDescription
GaussGreen-Gauss uncorrected gradient scheme
MidpointGreen-Gauss skew corrected scheme (2 iterations - hardcoded)

Gradient limiters

SchemeDescription
FaceBasedLimit gradient calculation based on face extrapolation
MFaceBasedLimit gradient calculation based on face extrapolation in multiple dimensions

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

XCALibre.Solve.SchemesType
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)
source

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.

XCALibre.Solve.SolverSetupType
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)
  • 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)
source

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
Note

Only the Jacobi and NormDiagonal preconditioners have GPU ready implementations. At present these have the most robust implementation and they can be used with both CPU and GPU backends. The other preconditioners can only be used on the CPU. Notice that on our tests the LDL preconditioner only works when paired with the Gmres() on the CPU. Also notice that the implementations for DILU (experimental), IC0GPU and ILU0GPU (NVIDIA only), should be considered experimental. Work on improving the offering of preconditioners is ongoing.

Warning

Internally the storage for sparse matrices was moved to the CSR format. Thus, temporarily, we have removed support for LDL and ILU0 preconditioners while we work on CSR compatible implementations.

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