OpenFOAM Linear Solvers: GAMG, PCG and smoothSolver Guide

The linear solver settings in fvSolution directly control how fast each pressure or velocity equation is solved per iteration — and whether it converges at all. Choosing the wrong solver or misconfiguring tolerances is a common cause of slow convergence or unexplained divergence. This guide explains when to use each solver and how to tune the key parameters.

By CFDpilot · Updated April 2026

1. The four main OpenFOAM linear solvers

GAMG — Geometric-Algebraic Multi-Grid

GAMG is the best choice for the pressure equation (p or p_rgh) on any mesh with more than ~50 000 cells. It builds a hierarchy of coarser representations of the system and solves on all levels, making it very efficient for elliptic equations like the pressure Poisson equation. It is not well-suited for strongly asymmetric matrices (e.g. highly convection-dominated transport equations).

GAMG works by progressively coarsening the matrix system: it groups cells together (agglomeration) to create a hierarchy of systems from fine to coarse, solves approximately on each level using a smoother, then corrects the fine-level solution using the coarse-level error estimate. This V-cycle or W-cycle approach means it converges in far fewer matrix-vector products than a simple iterative solver.

PCG — Preconditioned Conjugate Gradient

PCG works for symmetric positive-definite systems. It is a good alternative to GAMG for pressure on smaller meshes or when GAMG is causing instability. The most common preconditioner is DIC (Diagonal Incomplete-Cholesky).

PBiCGStab — Preconditioned Bi-Conjugate Gradient Stabilised

PBiCGStab is the standard solver for asymmetric systems such as the momentum equation (U) and turbulence transport equations (k, epsilon, omega). The preconditioner is typically DILU (Diagonal Incomplete-LU).

smoothSolver

smoothSolver applies a simple iterative smoother (GaussSeidel or symGaussSeidel) without a coarse-grid correction. It converges slowly but is very robust. Use it as a fallback for turbulence equations when PBiCGStab diverges, or inside GAMG as the smoother on coarse levels.

2. Typical fvSolution for incompressible steady-state (simpleFoam)

// system/fvSolution
solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        tolerance       1e-6;
        relTol          0.1;
    }
    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }
    k
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }
    epsilon
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }
}

The SIMPLE algorithm controls the outer pressure-velocity coupling loop. After the linear solvers converge for each equation, the SIMPLE corrector updates the pressure and velocity fields using relaxation factors. The linear solver tolerance determines inner loop accuracy; the SIMPLE residual controls accepts these approximate solutions as sufficient for the outer iteration.

3. GAMG key settings explained

nCellsInCoarsestLevel: The minimum number of cells in the coarsest grid level. A value of 4–20 is typical. Too small and the coarse solve becomes trivial and slow; too large and GAMG stops building levels early, losing its advantage.

agglomerator: Controls how cells are grouped to form coarser levels. faceAreaPair is the default and works well in most cases. algebraicPair can be faster on some meshes.

smoother: Applied within each GAMG level. GaussSeidel is standard. For symmetric systems, symGaussSeidel converges faster per sweep. DICGaussSeidel is more expensive per sweep but needs fewer iterations.

mergeLevels: Number of coarsening steps merged together. 1 (default) is safe; 2 can speed things up but may reduce robustness.

The number of GAMG levels created depends on the mesh size. OpenFOAM prints the number of levels to stdout when you run with the -v verbose flag. For a 1-million cell mesh you typically get 5–7 levels, which is near-optimal. Fewer than 4 levels on a large mesh indicates that nCellsInCoarsestLevel is too large.

// Additional GAMG options for difficult meshes
nPreSweeps      0;    // smoother sweeps before coarse solve
nPostSweeps     2;    // smoother sweeps after coarse solve
nFinestSweeps   2;    // sweeps on finest level
cacheAgglomeration true;  // cache coarsening for speed

4. tolerance vs relTol: understanding the difference

tolerance is the absolute residual threshold. The linear solver stops when the residual drops below this value. This is the final convergence criterion used when relTol is satisfied or when the run is near the final solution.

relTol is the relative reduction factor. The solver stops early if the residual is reduced by this fraction of its initial value. In early iterations, relTol = 0.1 means the solver does just enough work per outer iteration — avoiding wasted inner iterations when the outer residuals are still large.

// Recommended: relTol non-zero during solving, set to 0 for final tight convergence
tolerance   1e-7;   // absolute — final answer quality
relTol      0.1;   // relative — controls work per outer iteration

For final solutions or steady-state convergence studies, set relTol 0 to force the solver to always reach tolerance.

Typical tolerance values for production simulations:

• Pressure (p): tolerance 1e-6 to 1e-8, relTol 0.01 to 0.1 during the run, relTol 0 for the final corrector
• Velocity (U): tolerance 1e-6 to 1e-8, relTol 0.1
• Turbulence (k, epsilon, omega): tolerance 1e-6, relTol 0.1 — turbulence equations are less sensitive to tight linear convergence

5. Typical fvSolution for transient cases (pimpleFoam)

solvers
{
    p
    {
        solver          GAMG;
        smoother        DICGaussSeidel;
        tolerance       1e-6;
        relTol          0.05;
    }
    pFinal
    {
        $p;
        relTol          0;  // tight solve on final PIMPLE corrector
    }
    U
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }
    UFinal
    {
        $U;
        relTol          0;
    }
}

The pFinal and UFinal entries are used by pimpleFoam on the last PIMPLE corrector loop, ensuring tight convergence at each time step while saving cost on intermediate correctors.

In pimpleFoam, the PIMPLE dictionary in fvSolution controls the outer loop:

// system/fvSolution — PIMPLE block
PIMPLE
{
    nOuterCorrectors    2;   // outer PIMPLE loops per time step
    nCorrectors         2;   // inner pressure correctors
    nNonOrthogonalCorrectors 1;
}

With nOuterCorrectors 1, pimpleFoam degenerates to the PISO algorithm. With nOuterCorrectors > 1 and a large Courant number limit, pimpleFoam can take larger time steps than pure PISO while remaining stable.

6. fvSolution for compressible solvers (rhoPimpleFoam)

Compressible solvers add the enthalpy or internal energy equation. The solver configuration must include entries for the energy field e or h as well as the density field rho:

solvers
{
    rho
    {
        solver          diagonal;  // rho is computed explicitly
    }
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.05;
    }
    pFinal  { $p; relTol 0; }
    U
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }
    UFinal  { $U; relTol 0; }
    e
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }
    eFinal  { $e; relTol 0; }
}

The rho solver uses diagonal because density is computed directly from the equation of state (p = rho * R * T) rather than solved iteratively. Only the pressure, momentum, and energy equations require iterative solvers.

7. Performance impact of solver choice

On a 1-million cell RANS case, the solver choice for pressure typically accounts for 60–80% of total wall-clock time per iteration. Switching from PCG+DIC to GAMG+GaussSeidel can reduce solve time by a factor of 3–5. On highly non-orthogonal meshes, GAMG with DICGaussSeidel smoother is more robust than plain GaussSeidel.

If GAMG is slow or diverging, check:

• Is nCellsInCoarsestLevel appropriate? Try 4 or 20.
• Is the mesh highly non-orthogonal? Use nNonOrthogonalCorrectors 2 in SIMPLE or PIMPLE dict.
• For very small meshes (<10 000 cells), PCG+DIC may outperform GAMG.

You can measure solver performance by adding the solverInfo function object to system/controlDict:

functions
{
    solverInfo
    {
        type            solverInfo;
        libs            (utilityFunctionObjects);
        fields          (p U k epsilon);
        writeResiduals  true;
    }
}

This writes per-field residual histories to postProcessing/solverInfo/ as CSV files, which you can plot to identify which equation is limiting convergence and whether your tolerance settings are appropriate.