Runtime and solvers

Final steps before launching a simulation

Runtime setup

At this stage in the setup workflow, the user defines the runtime information to configure the runtime behaviour of the flow solver, including the time step to use (only meaningful for transient solutions), as well as information about how often to write results to disk. XCALibre.jl provides a the Runtime to perform this operation.

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

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

Configuration object

Once all the simulation configuration information has been defined, from discretisation scheme to runtime information, all settings must be wrapped in a Configuration object. The definition, including expected input arguments, for the Configuration object are detailed below.

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

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

Initialising fields

The last (optional) step before running the simulation is to provide an initial guess for all the fields being solved. Although this step is optional, in most cases the flow solvers will perform better when initialised. To set an initial value for a field, the initialise! function is provided, which assigns a starting value to a given field.

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

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

A uniform initial guess leaves a velocity field that does not satisfy continuity, which the first few solver iterations must work off before making physical progress. For external and internal flows around bodies, a potential-flow field is a much better starting point: it already satisfies continuity and carries the shape of the geometry.

potential_flow! solves a velocity-potential equation on the current mesh and sets the velocity to the resulting divergence-free field. By default it starts from rest: the internal velocity is set to zero first, so the flow is driven by the velocity boundary conditions alone (pass from_rest=false to project the current field instead). Boundary conditions for the potential are inferred from those already assigned to pressure, so no extra setup is needed.

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

Call it after initialise! and before run!:

initialise!(model.momentum.U, velocity)
initialise!(model.momentum.p, 0.0)

potential_flow!(model, config; ncorrectors=5)

residuals = run!(model, config)

On meshes with appreciable non-orthogonality, pass ncorrectors to run non-orthogonal correctors on the potential equation. Set it to zero for orthogonal meshes.

AMG preconditioner

XCALibre's native algebraic multigrid is a preconditioner: solver picks the outer iteration and preconditioner = AMG(...) the multigrid. Use solver = Cg() for symmetric systems such as pressure equations, solver = Bicgstab() for non-symmetric ones, and solver = Richardson() for a standalone multigrid solve (repeated V-cycles with no Krylov acceleration).

SolverSetup(
    solver = Cg(),
    preconditioner = AMG(
        coarsening = SmoothAggregation(),
        smoother = AMGJacobi()
    ),
    convergence = 1e-7,
    relax = 1.0,
    rtol = 0.0,
    atol = 1e-5
)

The same settings run on a distributed mesh, where PETSc's GAMG() takes the place of AMG() with a one-time warning (the AMG keywords are not carried over). The other way round, GAMG() or BoomerAMG() on a serial mesh runs the native AMG(). The rules for every preconditioner are listed in the preconditioners section of the numerical setup page.

For non-symmetric systems use solver = Bicgstab(), which gives AMG-preconditioned stabilised biconjugate gradient. A typical case is the pressure equation in compressible flow, where the implicit pressure convection is upwinded and therefore non-symmetric. Cg() rejects a non-symmetric matrix. examples/2D_cylinder_transonic_RANS_AMG_BICGStab.jl uses it for the pressure equation of a transonic RANS case.

Choice of coarsening for non-symmetric operators: Bicgstab() has been tested with the default SmoothAggregation(). RugeStuben() and Geometric() are accepted but have not been validated on non-symmetric operators - RugeStuben() builds its strength of connection from each row's entries only, not from the transpose - so prefer the default unless you have checked the alternative on your case.

If Bicgstab() breaks down - the shadow residual going orthogonal, or the stabilising step collapsing - the half-step iterate is kept and the solve restarts with a fresh shadow residual, at most twice.

scale_correction (on by default) is supported. It makes the V-cycle preconditioner depend on its input, but the solver is right-preconditioned, so the solution and residual it carries stay consistent, and convergence is confirmed against the true residual b - A*x.

SolverSetup(
    solver = Bicgstab(),
    preconditioner = AMG(smoother = AMGGaussSeidel(sweep = AMGForwardSweep())),
    convergence = 1e-7,
    relax = 1.0,
    rtol = 0.0,
    atol = 1e-5     
)

Solving the coarsest level

The coarsest level of the hierarchy is small but is visited on every V-cycle, so on a GPU the default host round trip can dominate. coarse_solve selects how it is solved:

• OnDevice(; max_rows = 512) (default) — device-resident factorisation, no host round trip. A coarsest level with more rows than max_rows falls back to the host, since the factor is rebuilt on every coarse refresh.
• OnHost() — always copy the coarsest right-hand side to the host and solve there. Also written coarse_solve = CPU().
• OnDeviceJacobi(; omega, iterations) and OnDeviceChebyshev(; degree, eig_ratio, lambda_scale) — apply a fixed number of smoother sweeps on the device instead of solving. Both are constant linear operators, so they are valid with solver = Cg() as well as solver = Richardson().
• OnDeviceKrylov(; solver, rtol, atol, itmax) — run an inner Krylov solve on the device. An inner Krylov solve is nonlinear in its right-hand side, which breaks the fixed-operator assumption of preconditioned CG, so this one requires solver = Richardson() and is rejected with solver = Cg(). Passing a GPU backend, as in coarse_solve = CUDABackend(), selects it.

SolverSetup(solver = Richardson(), preconditioner = AMG(coarsening = SmoothAggregation(),
    smoother = AMGJacobi(), coarse_solve = OnDeviceChebyshev(degree = 8)), convergence = 1e-7, relax = 1.0)

On a CPU backend every OnDevice* option behaves as OnHost().

Launching flow solvers

In XCALibre.jl the run! function is used to start a simulation, which will dispatch to the appropriate flow solver for execution. Once the simulation is complete a NamedTuple containing residual information is returned for users to explore the convergence history of the simulation.

The residual of each field is the normalised residual Σᵢ|bᵢ - (Ax)ᵢ| / Σᵢ|Aᵢᵢxᵢ| of its assembled system, evaluated on the current field x just before the linear solve, so it measures how far the previous iteration's fields are from satisfying this iteration's equations. It does not depend on the linear solver, its preconditioner, precision or tolerances, nor on whether the run is serial, distributed or on a GPU. When the field is zero everywhere (a common initial pressure), the denominator is zero and the unnormalised Σᵢ|bᵢ - (Ax)ᵢ| is reported for that iteration.

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

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

residuals = run!(model, config) 

# to access the pressure residual

residuals.p 

run! shows a progress bar with the current iteration and residuals. Pass progress=false to turn it off:

residuals = run!(model, config; progress=false)

We recommend progress=false for large-scale runs that are not on a local PC, such as batch jobs on a cluster, where the progress bar only fills the job's log. The residual history returned by run! is the same either way.

The momentum equation of every flow solver carries the full viscous stress, ∇·τ with τ = μeff (∇U + (∇U)ᵀ - ⅔ (∇·U) I). The part ∇·(μeff ∇U) is discretised implicitly as a Laplacian, and the rest, ∇·(μeff dev2((∇U)ᵀ)) with dev2(A) = A - ⅔ tr(A) I, is added as an explicit source evaluated from the current velocity gradient. For a constant viscosity in incompressible flow this source is zero, since ∇·(∇U)ᵀ = ∇(∇·U); it matters wherever the effective viscosity varies, as it does with any turbulence model. At boundary faces the source uses the owner-cell gradient with its normal derivative taken from the velocity condition: (Ub − Uc)/δ for fixed values (Wall, Dirichlet, RotatingWall, and DirichletFunction evaluated at the face and the current time), zero for Zerogradient, Extrapolated, Outlet and NeumannFunction, and the mirrored value for Slip and Symmetry; periodic faces are treated as internal faces between the two cells of a pair. Pass `transposestress=false` to leave it out, e.g. to compare with results obtained before it was included:

residuals = run!(model, config; transpose_stress=false)

The steady solvers (simple! for Incompressible fluids and csimple! for Compressible and WeaklyCompressible fluids, used by run! for Steady() models) offer a consistent (SIMPLEC) pressure–velocity coupling. By default, the pressure equation and the velocity correction use rD = V/aP, the inverse diagonal coefficient of the under-relaxed momentum matrix. Because aP grows as the velocity relaxation factor αU decreases, rD scales roughly with αU, so the pressure correction is over-damped and has to be compensated by pressure under-relaxation. With consistent=true the pressure Laplacian and the velocity and face-flux corrections use rDt = V/(aP + Σ aN), the inverse row sum of the momentum matrix, which does not shrink with αU; the extra terms (rDt - rD)∇p are added to the predicted velocity and face fluxes, so the converged solution is the same. Incompressible runs then need little or no pressure under-relaxation. In compressible flow the weak Laplacian of the default coupling also loses weight against the compressibility term ∇·(ψ ((H/aP)·S) p) of a Compressible fluid; in transonic flow this can amplify round-off from one iteration to the next until the run diverges, even for a uniform free stream on some meshes. The consistent coupling keeps the Laplacian's weight and is worth trying when a steady transonic run diverges. The option works on distributed meshes with simple!:

residuals = run!(model, config; consistent=true)

It is off by default, which leaves earlier results unchanged. The transient solvers do not take this option.

Restarting simulations

It should be noted that when running a simulation with run!, the solution fields in the Physics model are mutated. Thus, running the simulation from the previous solution is simply a matter of reissuing the run! function. The continued run also picks up the face mass flux of the last iteration, so a run split into several run! calls (for example, to monitor forces every few hundred iterations) retraces the uninterrupted run. Re-initialising the fields with initialise! between calls starts afresh. At present, this has the side effect of overwriting any existing solution files (.vtk or .vtu). Users must be aware of this behaviour.

In some cases, it may be desirable to solve a problem with a steady solver and use the solution to run transient simulations. This is possible using the change function.

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

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