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
Scheme | Description |
---|---|
SteadyState | sets the time derivative to zero |
Euler | First order implicit Euler scheme |
CrankNicolson | Second order central scheme (not implemented yet) |
Laplacian schemes
Scheme | Description |
---|---|
Linear | 2nd order Gauss gradient scheme with linear interpolation |
Divergence schemes
Scheme | Description |
---|---|
Linear | 2nd order central difference |
Upwind | 1st order upwind scheme |
BoundedUpwind | Bounded version of the Upwind scheme |
LUST | 1st/2nd order mixed scheme (fixed at 75% Linear - 25% Upwind) |
Gradient schemes
Scheme | Description |
---|---|
Gauss | Green-Gauss uncorrected gradient scheme |
Midpoint | Green-Gauss skew corrected scheme (2 iterations - hardcoded) |
Gradient limiters
Scheme | Description |
---|---|
FaceBased | Limit gradient calculation based on face extrapolation |
MFaceBased | Limit 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.Schemes
— TypeSchemes(;
# 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 isSteadyState
)divergence
: is used to set the divergence scheme (default isLinear
)laplacian
: is used to set the laplacian scheme (default isLinear
)gradient
: is used to set the gradient scheme (default isGauss
)limiter
: is used to specify if gradient limiters should be used, currently supported limiters includeFaceBased
andMFaceBased
(default isnothing
)
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. forU
.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.SolverSetup
— TypeSolverSetup(;
# 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 ofBicgstab()
,Cg()
,Gmres()
which are re-exported in XCALibre.jlpreconditioner
: instance of preconditioner to be used e.g. Jacobi()convergence
sets the stopping criteria of this fieldrelax
: specifies the relaxation factor to be used e.g. set to 1 for no relaxationsmoother
: specifies smoothing method to be applied before discretisation.JacobiSmoother
: is currently the only choice (defaults tonothing
)limit
: used in some solvers to bound the solution within these limits e.g. (min, max). It defaults tonothing
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 toFloat64
)
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
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.
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
)
)