Choosing a linear solver

Which linear solver settings suit the implicit schemes, RBF mesh motion and the incompressible solver?

Start from the defaults: they are tuned for each equation, and most cases need no change. The library is set by the hardware. On CPUs zCFD uses Petsc, configured from the control dictionary; on NVIDIA GPUs it uses AMGX, configured from json files read at run time. Change a setting only for a measured reason, such as an outer iteration that stalls, and then only for the one equation concerned. The two backends are tuned independently, so a setting that helps one is not evidence for the other.

A linear system is solved in these cases:

The mean-flow solve of "mf-gmres" is matrix-free GMRES, configured on the scheme itself (see mf-gmres), and takes nothing from these settings. "lu-sgs" and "mf-gmres" also take the assembled turbulence solve wherever the backend forces it: on a CUDA build running on the host only, and on HIP.

Both libraries offer a wide range of solvers and preconditioners, and the settings can make a large difference to the performance and convergence of the solver. For the settings each backend uses, equation by equation, see Side by side: Petsc and AMGX.

The Petsc settings sit under linear_solver_options, in solver_settings:

parameters = {
    "config_version": 2,
    "solver": {
        ...
        "solver_settings": {
            "type": "compressible",
            "linear_solver_options": {
                "flow": {...},
                "turbulence": {...},
                "rbf": {...},
                "double_precision": False,
            },
        },
        ...
    },
    ...
}

Petsc

The options available for the Petsc KSP linear solvers are described in the Petsc KSP manual. zCFD passes them to Petsc from the linear_solver_options dictionary.

Each of flow, turbulence and rbf is a plain dictionary of Petsc command-line flags, with both key and value strings. Where a Petsc option takes no value, give an empty string. Most settings in the Petsc documentation can be passed in this way. The compressible defaults are:

"linear_solver_options": {
    "flow": {
        "-ksp_min_it": "2",
        "-ksp_type": "fgmres",
        "-ksp_rtol": "1.0e-3",
        "-ksp_converged_reason": "",
    },
    "turbulence": {
        "-ksp_min_it": "2",
        "-ksp_type": "fgmres",
        "-ksp_rtol": "1.0e-5",
        "-ksp_converged_reason": "",
    },
    "rbf": {
        "-ksp_type": "fgmres",
        "-ksp_rtol": "1.0e-5",
    },
    "double_precision": False,
}

-ksp_converged_reason prints one line for each linear solve, giving the number of iterations it took and whether it converged. A solve that reports DIVERGED_ITS ran out of iterations rather than diverging. It still returns its best answer, but the tolerance asked for was not met, and a run in which that happens on most cycles is doing work the outer iteration cannot use. -ksp_monitor prints a line for each Krylov iteration. It is useful when investigating a solve, but is not a default, because on a long run it dominates the log; add it to any of these dictionaries when it is needed.

One options namespace per solver

Petsc has a single options database for the whole process. zCFD builds a separate solver for each system it solves, and gives each its own Petsc options prefix, taken from the name of the dictionary it comes from:

So "-ksp_type": "bcgs" in the momentum dictionary reaches Petsc as -momentum_ksp_type bcgs, and applies to that solver alone. The keys in a control file are written without the prefix, and zCFD adds it. The prefix matters only for options given from outside the control file, through the PETSC_OPTIONS environment variable for instance, where the prefix must be written out to reach a single solver.

The Hypre library of algebraic multigrid methods is available through Petsc PCHYPRE, and its settings are given in the same dictionaries. Hypre's BoomerAMG package, as a preconditioner, has performed well with the incompressible solver:

"linear_solver_options": {
    "flow": {
        "-ksp_type": "fgmres",
        "-ksp_rtol": "1.0e-3",
        "-ksp_monitor": "",
        "-ksp_converged_reason": "",
        "-pc_type": "hypre",
        "-pc_hypre_type": "boomeramg",
    },
    ...
}

Further Hypre options are described in the Petsc PCHYPRE documentation.

Incompressible solver

With {"solver_settings": {"type": "incompressible"}}, linear_solver_options gains two further keys, momentum and pressure_correction, and takes its own defaults, modelled on OpenFOAM's simpleFoam.

Why they look the way they do

Under the SIMPLE scheme every equation is solved inside an outer iteration, so its linear solve need not be accurate. It needs to be cheap, because the outer loop revisits it a few thousand times. That is what simpleFoam encodes: across the OpenFOAM tutorials every equation is solved to relTol 0.1, one order of magnitude, and only the kind of solver changes with the equation. The motorBike tutorial is the exception, tightening p to 0.01; the zCFD pressure-correction default follows that tighter figure, and every other equation the 0.1. zCFD follows the same mapping:

The pressure correction is the one equation that needs coarse levels. It is elliptic, and a purely local preconditioner such as ILU cannot carry information across the mesh however strong it is made, so it exhausts its iteration limit without converging. Algebraic multigrid is the cure, and is what OpenFOAM uses for p.

The coupled scheme is the exception. There the block-4 velocity-pressure solve is the step itself rather than an inner iteration, so flow keeps a tight tolerance. simpleFoam has no counterpart to it.

The defaults:

"linear_solver_options": {
    # Coupled scheme only: this solve is the step, so it is solved tightly.
    "flow": {
        "-ksp_min_it": "2",
        "-ksp_type": "fgmres",
        "-ksp_rtol": "1.0e-7",
        "-ksp_atol": "1.0e-16",
        "-ksp_max_it": "50",
        "-pc_type": "bjacobi",
        "-sub_pc_type": "ilu",
        "-ksp_converged_reason": "",
    },
    # simpleFoam's k/omega/epsilon/nuTilda entry. Also used for the energy equation.
    "turbulence": {
        "-ksp_min_it": "1",
        "-ksp_type": "bcgs",
        "-ksp_rtol": "0.1",
        "-ksp_atol": "1.0e-16",
        "-ksp_max_it": "50",
        "-pc_type": "bjacobi",
        "-sub_pc_type": "ilu",
        "-ksp_converged_reason": "",
    },
    # SIMPLE only: simpleFoam's U entry.
    "momentum": {
        "-ksp_min_it": "1",
        "-ksp_type": "bcgs",
        "-ksp_rtol": "0.1",
        "-ksp_atol": "1.0e-12",
        "-ksp_max_it": "50",
        "-pc_type": "bjacobi",
        "-sub_pc_type": "ilu",
        "-ksp_converged_reason": "",
    },
    # SIMPLE only: simpleFoam's p entry.
    "pressure_correction": {
        "-ksp_min_it": "1",
        "-ksp_type": "cg",
        "-ksp_rtol": "0.01",
        "-ksp_atol": "1.0e-12",
        "-ksp_max_it": "100",
        "-pc_type": "hypre",
        "-pc_hypre_type": "boomeramg",
        "-zcfd_constant_nullspace": "1",
        "-zcfd_pc_reuse_cycles": "20",
        "-ksp_converged_reason": "",
    },
    "rbf": {
        "-ksp_type": "fgmres",
        "-ksp_rtol": "1.0e-5",
    },
    "double_precision": True,
}

-zcfd_constant_nullspace is not a Petsc option: zCFD consumes it and never passes it on. It attaches a constant null space to the pressure-correction operator, which keeps the solve well posed when the domain is enclosed and the pressure has no datum. zCFD detaches it again if the boundary conditions include a fixed-pressure outlet, since that outlet then sets the datum.

-zcfd_pc_reuse_cycles is not a Petsc option either; zCFD consumes it too. It sets how often the preconditioner is rebuilt. The matrix is reassembled at every outer iteration, so by default Petsc rebuilds the multigrid hierarchy before every solve, and that is the dominant cost. Profiled with -log_view on a 52k-cell plate SIMPLE case, PCSetUp took 5.01 s of the 8.50 s spent in KSPSolve, 59% of all linear-solver time, for a matrix whose sparsity never changes. Building it once and freezing it is no answer either: on the plate case that takes the pressure solve from 1.26 iterations per cycle over the first fifty cycles to 8.3 over the last fifty, as the start-up transient moves the operator away from it. Rebuilding every N solves cuts the setup cost by a factor of N, while keeping the preconditioner at most N cycles stale; it is the trade OpenFOAM's GAMG makes with cacheAgglomeration. Set it to 1 to rebuild for every solve. Whatever the value, it cannot change the answer: the Krylov solver applies the true operator, so a stale preconditioner costs iterations and nothing else.

pressure_correction and momentum apply only when the incompressible solver uses the SIMPLE coupling scheme, {"convergence_control": {"scheme": {"name": "simple"}}}. The coupled scheme builds neither solver, and setting either key under it is rejected.

AMGX

The AMGX settings for the mean-flow, turbulence and RBF linear systems are read from json files, not from the control dictionary: ZCFD_HOME/amgx.json, ZCFD_HOME/turbamgx.json and ZCFD_HOME/RBF_amgx.json. The files are in AMGX's own configuration format, and are unrelated to the config_version key of the zCFD control file; the "config_version": 2 below is the version of AMGX's file format, not zCFD's. The mean-flow settings are:

{
  "config_version": 2,
  "determinism_flag": 1,
  "solver": {
    "preconditioner": {
      "error_scaling": 0,
      "print_grid_stats": 0,
      "max_uncolored_percentage": 0.05,
      "algorithm": "AGGREGATION",
      "solver": "AMG",
      "smoother": "MULTICOLOR_GS",
      "presweeps": 0,
      "selector": "SIZE_8",
      "coarse_solver": "NOSOLVER",
      "max_iters": 1,
      "postsweeps": 3,
      "min_coarse_rows": 32,
      "relaxation_factor": 0.75,
      "scope": "amg",
      "max_levels": 40,
      "matrix_coloring_scheme": "PARALLEL_GREEDY",
      "cycle": "V"
    },
    "use_scalar_norm": 1,
    "solver": "FGMRES",
    "print_solve_stats": 1,
    "obtain_timings": 1,
    "max_iters": 10,
    "monitor_residual": 1,
    "gmres_n_restart": 5,
    "convergence": "RELATIVE_INI_CORE",
    "scope": "main",
    "tolerance": 1e-3,
    "norm": "L2"
  }
}

These settings perform well for most cases. If the linear system is hard to converge, reducing "tolerance" may help. Otherwise, changing the aggregation selector from "selector": "SIZE_8" to "SIZE_4" or "SIZE_2" builds smaller aggregates, and so more coarse levels, improving convergence at the cost of memory. Increasing "gmres_n_restart" may also improve convergence, again at the cost of memory.

AMGX has too many options to describe here, but the AMGX example configurations cover the different solvers.

Incompressible AMGX configurations

The incompressible solver reads its own set of files, one for each equation, so that each system has a configuration suited to it. Each file is looked for first in the working directory and then in ZCFD_HOME, so a copy placed beside the control file overrides it for that case.

"structure_reuse_levels": -1 in incomp_amgx.json and incomp_pcorr_amgx.json keeps the multigrid hierarchy from cycle to cycle. The matrix coefficients change at every outer iteration but the sparsity pattern never does, so re-aggregating every cycle would be wasted work. Leave "coarse_solver": "NOSOLVER" in incomp_pcorr_amgx.json as it is: with no fixed-pressure outlet the pressure-correction operator is singular, and a direct coarse solve meets a zero pivot on it. That is also why OpenFOAM leaves directSolveCoarsest off.

On NVIDIA GPUs, with AMGX, the one linear-solver setting that still comes from the zCFD control dictionary is double_precision, on linear_solver_options:

"linear_solver_options": {"double_precision": False}

double_precision sets whether the mean-flow linear system is solved in single or double precision. Single precision reduces memory use and speeds up the solver; more demanding cases may need double precision.

Side by side: Petsc and AMGX

The same equation is solved by a different library on CPU and on GPU, and the two are configured in different places: Petsc through the control dictionary, AMGX through its json files. This section sets them side by side, equation by equation, so that the settings a case is running can be read in one place.

The two backends are tuned independently, and are not expected to match. Where they differ, the difference is usually deliberate and measured; the note under Incompressible AMGX configurations gives the two cases where a Petsc setting was tried on AMGX and rejected.

How the settings correspond

The tolerance rows depend on the convergence mode, which is why the same tolerance key means different things in different files:

• COMBINED_REL_INI_ABS, used by the incompressible files, stops when |r| < tolerance or |r|/|r0| <= alt_rel_tolerance. tolerance is then the absolute floor and alt_rel_tolerance the relative one: the two Petsc tolerances, in the opposite order to the one the names suggest.
• RELATIVE_INI_CORE, used by the compressible files, stops when |r|/|r0| <= tolerance. Here tolerance is the relative tolerance, and there is no absolute floor at all.

Incompressible, SIMPLE scheme

Momentum predictor — momentum / prefix momentum_ / incomp_simple_mom_amgx.json:

Pressure correction — pressure_correction / prefix pcorr_ / incomp_pcorr_amgx.json:

This is the closest-matched pair: multigrid at relTol 0.01 on both sides.

Turbulence and energy — turbulence / prefixes turb_ and energy_ / incomp_turb_amgx.json:

The relative tolerances differ by five orders of magnitude, by design; see the note under Incompressible AMGX configurations. The energy equation shares the turbulence dictionary and the same json file on both backends; only its Petsc options prefix is its own.

Incompressible, coupled scheme

Block-4 velocity-pressure — flow / prefix flow_ / incomp_amgx.json. Turbulence and energy are as above; momentum and pressure_correction are not built at all under this scheme.

Compressible

Mean flow: flow / prefix flow_ / amgx.json. Turbulence: turbulence / prefix turb_ / turbamgx.json, or turbamgx_sa.json when the turbulence model is sa-neg.

The tolerances match exactly here. The iteration limits do not: a compressible Petsc solve is bounded only by its tolerance, whereas the AMGX one stops after ten iterations, converged or not.

RBF mesh motion

rbf / prefix rbf_ / RBF_amgx.json. This is the least closely matched pair in zCFD: the two backends use a different Krylov method, a different cycle and different kinds of tolerance.

Summary

Reading the settings a case actually used

On the Petsc path the solver logs every option it sets, under the prefix it sets it under, so the log is the authoritative record of what a run used, including anything overridden from the control dictionary. On the AMGX path, set "print_solve_stats": 1 in the relevant json file to have each solve report its iteration count and residuals.

Reference entries

• linear_solver_options: flow, turbulence, momentum, pressure_correction, rbf and double_precision
• scheme and turbulence_solver