Choosing a CFD Linear Solver: Comparing Convergence Speeds from Jacobi to BiCGSTAB

In 1952 Lanczos, Hestenes, and Stiefel walked into the same meeting with the same algorithm, written apart. Symmetric positive-definite (SPD) matrices, finished in $N$ iterations. The name that stuck is conjugate gradient (CG). Seventy years later the pressure equation in a CFD code still hunts for its answer in that subspace.

Drag grid size N from 20 to 80 and the four residual curves are drawn again from scratch. Raise the iteration cap and the order in which they hit the floor changes.

Near N=40, CG kinks downward while Jacobi lies almost flat on the log axis. Push N to 80 and the gap widens. The teeth on BiCGSTAB are a residual that will not die smoothly.

Why CFD Solves Linear Systems Iteratively#

Discretize Navier–Stokes and every step drops $A\mathbf{x}=\mathbf{b}$. $A$ is the coefficient matrix, $\mathbf{x}$ the cell-centered unknowns, $\mathbf{b}$ the right-hand side. A million cells makes $A$ a sparse $10^6 \times 10^6$ matrix. LU factorization, the split into lower and upper triangular factors, burns memory and time together at that size. OpenFOAM, SU2, and Gerris iterate.

Start from a guess $\mathbf{x}_0$ and keep editing $\mathbf{x}$ until the residual $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ drops under the tolerance. The definition of "a little" is the solver's name.

The Idea of the Krylov Subspace#

Repeated products of $A$ with the initial residual $\mathbf{r}_0$ fill a space.

$\mathcal{K}_m$ is the Krylov subspace, $m$ its dimension, $\text{span}$ the set of linear combinations. The update stays in $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$, and the leftover $\mathbf{r}_m$ has to meet an optimality condition.

The only required operation is a matrix-vector product, so the sparsity pattern survives. CG minimizes the A-norm error in that space. GMRES minimizes the 2-norm residual. The second one happens in the subspace Arnoldi builds by stacking vectors.

Differences in Solver Characteristics#

The pressure Poisson problem $-\nabla^2 p = f$ is SPD, so CG fits. Momentum with advection is nonsymmetric, so the choice moves to BiCGSTAB or GMRES. OpenFOAM's PCG/PBiCGStab and SU2's BCGSTAB are that split, written into the option names.

Convergence Speed with Python#

The model problem is 1D Poisson, $-u'' = 1$, with $u(0)=u(1)=0$. Finite differences leave a tridiagonal matrix. At $N=60$, Jacobi and CG share the tolerance $10^{-10}$.

import numpy as np
 
def build_poisson_1d(n):
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = 2.0
        if i > 0:     A[i, i - 1] = -1.0
        if i < n - 1: A[i, i + 1] = -1.0
    return A, np.ones(n)
 
def jacobi_sweep(A, b, tol=1e-10, max_iter=2000):
    D  = np.diag(A)
    LU = A - np.diag(D)
    x  = np.zeros_like(b)
    b0 = np.linalg.norm(b)
    history = []
    for _ in range(max_iter):
        x = (b - LU @ x) / D
        res = np.linalg.norm(b - A @ x) / b0
        history.append(res)
        if res < tol:
            break
    return x, history
 
def cg_sweep(A, b, tol=1e-10, max_iter=2000):
    x  = np.zeros_like(b)
    r  = b - A @ x
    p  = r.copy()
    rr = r @ r
    b0 = np.linalg.norm(b)
    history = [np.sqrt(rr) / b0]
    for _ in range(max_iter):
        Ap    = A @ p
        alpha = rr / (p @ Ap)
        x    += alpha * p
        r    -= alpha * Ap
        rr_new = r @ r
        history.append(np.sqrt(rr_new) / b0)
        if history[-1] < tol:
            break
        p  = r + (rr_new / rr) * p
        rr = rr_new
    return x, history
 
A, b = build_poisson_1d(n=60)
_, hist_j  = jacobi_sweep(A, b)
_, hist_cg = cg_sweep(A, b)
print(f"Jacobi: {len(hist_j):4d} iterations, final {hist_j[-1]:.2e}")
print(f"CG:     {len(hist_cg):4d} iterations, final {hist_cg[-1]:.2e}")

The run prints this.

Jacobi: 2000 iterations, final 2.11e-02
CG:       60 iterations, final 1.87e-13

Raise the log-condition slider and the four curves stall in a different place inside one iteration budget. The iteration cap decides where the plot is cut.

At the default $\kappa=10^3$ and about 80 iterations, Jacobi (red) is nearly horizontal and CG (cyan) is steeper. Push $\kappa$ toward $10^6$ and Jacobi still has no floor at 200 iterations. The printed $N=60$ case is one point on that picture: Jacobi at 2.11e-2 after 2000 sweeps, CG at 1.87e-13 after 60.

Turn the printed residual into a per-step factor. From 1 down to 2.11e-2 in 2000 sweeps is a contraction of 0.998 per sweep, the 2000th root of 0.0211. That is the scale of the spectral radius $1-O(1/N^2)$ in the table. At $N=60$, $1/N^2$ is 1/3600, about 0.00028, and 0.998 is that term with a small constant in front.

Keep the factor and $10^{-10}$ costs about 12,000 iterations. From 2.11e-2, about 10,000 sweeps are still missing. Double $N$ and $1/N^2$ falls by four, so the same residual costs four times as many iterations.

CG's 60 iterations equal $N$, so the ceiling of at most $N$ steps is used up. 1.87e-13 is about $1.1 imes 10^{11}$ times smaller than 2.11e-2. Double precision adds about $10^{-16}$ to 1, and two or three digits wore off. A million cells at 8 bytes is 8 MB a vector: 32 MB for four CG vectors, 56 MB for seven BiCGSTAB vectors. Storage is cheap beside the matrix. The tail is not.

Bare iterations, preconditioner off#

Those curves are bare iterations. Turn a preconditioner on and one tick on the horizontal axis means something else. A vertical spike in BiCGSTAB is breakdown: an inner product in a denominator has stuck at zero. More iterations do not walk that back. Change the preconditioner, or move to GMRES. A restart period $m$ near 30 is the usual compromise between memory and that spike.

Practical Solver Selection Guide#

Pressure Poisson, and heat transfer that is only diffusion: CG with ILU(0). If symmetry is intact, another name just stores more vectors.

Nonsymmetric momentum starts on BiCGSTAB. When the residual wobbles, switch to GMRES($m$). Here $m$ is the restart period.

Gauss-Seidel is sequential, so it goes limp on a GPU. The pattern of the last decade is a Jacobi smoother plus algebraic multigrid (AMG).

SU2 also carries LU-SGS (Lower-Upper Symmetric Gauss-Seidel), a stationary iteration beside Krylov. Tied to implicit time stepping on compressible flow, its memory footprint stays smaller. It is not as sturdy as Krylov.

The preconditioner moves the iteration count more than the solver name does. PCG plus AMG is often 100 times faster than bare CG. That factor is whether the ten-thousand-sweep tail gets cut.

Calling 2.11e-2 a converged residual#

Hand the next time step a Poisson residual that is still 2.11e-2 after 2000 sweeps, and that pressure error walks into the next nonlinear iteration. While the matrix stays SPD, the slot belongs to CG. Once advection tips the matrix nonsymmetric, CG is no longer eligible and the handoff is BiCGSTAB. Picking the name first and measuring the condition number later is the failure.

The open question is whether to stack another ten thousand iterations. A V-cycle that clears the error on a coarse grid may cut this $1-O(1/N^2)$ tail, or it may call Jacobi again as the fine-grid smoother. This residual does not answer that.