Solvers & Convergence

After discretizing the governing equations, we have a massive system of linear equations. For a 3D mesh with 10 million cells solving for 5 variables, that's 50 million unknowns. Direct solution methods like Gaussian elimination have $O(n^3)$ complexity — utterly impractical. Instead, CFD relies on iterative solvers that progressively refine an initial guess.

The Linear System

Matrix Form

The discretized equations form:

$$\mathbf{A}\mathbf{x} = \mathbf{b}$$

Sparsity Pattern

For a 5-point stencil in 2D (or 7-point in 3D), each equation involves only a few neighbors:

$$a_P \phi_P = a_E \phi_E + a_W \phi_W + a_N \phi_N + a_S \phi_S + b_P$$

Basic Iterative Methods

Jacobi Method

The simplest iterative approach: solve for each unknown using values from the previous iteration.

Algorithm:

$$\phi_P^{n+1} = \frac{1}{a_P}\left( a_E \phi_E^n + a_W \phi_W^n + a_N \phi_N^n + a_S \phi_S^n + b_P \right)$$

Gauss-Seidel Method

Use the latest available values immediately:

$$\phi_P^{n+1} = \frac{1}{a_P}\left( a_E \phi_E^n + a_W \phi_W^{n+1} + a_N \phi_N^n + a_S \phi_S^{n+1} + b_P \right)$$

Properties:

Faster convergence than Jacobi (about 2x)
Sequential updates reduce parallelism
Sweep direction matters (influences which neighbors are "new")

Successive Over-Relaxation (SOR)

Accelerate Gauss-Seidel with an over-relaxation factor $\omega$:

$$\phi_P^{n+1} = (1-\omega)\phi_P^n + \omega \phi_P^{GS}$$

Where $\phi_P^{GS}$ is the Gauss-Seidel update.

Optimal $\omega$:

$\omega = 1$: Standard Gauss-Seidel
$1 < \omega < 2$: Over-relaxation (faster)
$\omega \approx 1.5 - 1.9$: Typical for CFD
$\omega \geq 2$: Unstable

Residual Monitoring

Observe how residuals drop during iteration and identify convergence patterns.

What Is a Residual?

The residual measures how well the current solution satisfies the equations:

$$r_P = b_P - a_P \phi_P + \sum a_{nb} \phi_{nb}$$

For a perfect solution, $r_P = 0$ everywhere.

Scaled Residuals

Raw residuals depend on variable magnitudes. Scaled residuals are normalized:

$$R_\phi = \frac{\sum_{cells} |r_P|}{\sum_{cells} |a_P \phi_P|}$$

Alternatively, normalize by inlet flux or other reference values.

Typical Convergence Criteria

Variable	Typical Target
Continuity	$10^{-4}$ to $10^{-6}$
Momentum	$10^{-4}$ to $10^{-6}$
Energy	$10^{-6}$ to $10^{-8}$
Turbulence ($k$, $\epsilon$)	$10^{-4}$ to $10^{-5}$

Monitoring Beyond Residuals

Residuals alone don't guarantee accuracy. Also monitor:

Mass balance — Net mass flow should be zero (closed domain) or match inlet (open)
Key quantities — Drag coefficient, pressure drop, exit temperature
Point monitors — Values at specific locations should stabilize

Convergence Patterns

Healthy Convergence

Residuals drop monotonically
3-4 orders of magnitude reduction
Key quantities stabilize

Warning Signs

Pattern	Possible Cause
Stalled residuals	Under-relaxation too low, mesh issues
Oscillating residuals	Under-relaxation too high, transient physics
Divergence	Bad mesh, wrong BC, too aggressive settings
Slow convergence	Poor initial guess, highly coupled physics

Multigrid Methods

See how multigrid accelerates convergence by solving on multiple grid levels.

The Problem with Smooth Errors

Basic iterative methods eliminate high-frequency errors (cell-to-cell oscillations) quickly but struggle with low-frequency errors (smooth, gradual variations). On fine meshes, smooth errors look like constants locally.

The Multigrid Idea

Low-frequency errors on a fine grid become high-frequency on a coarser grid. The multigrid strategy:

Smooth on fine grid (removes high-frequency errors)
Restrict residual to coarser grid
Solve (approximately) on coarse grid
Prolongate correction back to fine grid
Smooth again on fine grid

V-Cycle Algorithm

multigrid(level):
    if level == coarsest:
        solve_directly()
    else:
        smooth(level)                    // Pre-smoothing
        residual = compute_residual()
        coarse_residual = restrict(residual)
        coarse_correction = multigrid(level+1)
        correction = prolongate(coarse_correction)
        apply_correction(correction)
        smooth(level)                    // Post-smoothing

Why Multigrid Is Optimal

For a grid with $N$ cells:

Basic iterative methods: $O(N^2)$ operations to converge
Multigrid: $O(N)$ operations — optimal scaling

This is why multigrid is essential for large-scale CFD.

Algebraic Multigrid (AMG)

Geometric multigrid requires explicit coarse grids. Algebraic multigrid constructs coarse levels directly from matrix coefficients:

Works on unstructured meshes
Automatic coarsening
Used in most commercial solvers

Pressure Equation Challenges

The pressure correction equation (Poisson-type) is often the hardest to solve:

Elliptic nature — Information propagates everywhere
No diagonal dominance — Near-singular for some BCs
Slow convergence — Low-frequency modes dominate

Solution: Always use multigrid or AMG for pressure. SIMPLE iterations per outer loop should be minimized.

Solver Selection

Equation Type Matters

Equation Type	Recommended Solver
Pressure (Poisson)	AMG, Multigrid
Momentum	ILU + BiCGSTAB, GMRES
Scalars (T, k, ε)	Gauss-Seidel, ILU

Krylov Subspace Methods

Modern CFD often uses Krylov methods:

Conjugate Gradient (CG) — Symmetric positive definite only
BiCGSTAB — General non-symmetric systems
GMRES — Robust, needs memory for basis vectors

Combined with preconditioners (ILU, AMG), these are state-of-the-art.

Under-Relaxation Revisited

Why It Helps Convergence

The SIMPLE algorithm over-estimates corrections. Under-relaxation stabilizes:

$$\phi^{new} = \phi^{old} + \alpha (\phi^{calc} - \phi^{old})$$

Effect on Convergence

Under-relaxation	Effect
Too low ($\alpha = 0.1$)	Stable but extremely slow
Optimal	Fast, stable convergence
Too high ($\alpha = 0.9$)	Oscillations or divergence

Adaptive Under-Relaxation

Some solvers adjust $\alpha$ based on residual behavior:

Increase if converging smoothly
Decrease if oscillating

Practical Convergence Strategy

Starting a Difficult Case

Initialize well — Use potential flow or steady boundary values
Start conservative — Low under-relaxation, first-order schemes
Ramp up gradually — Increase $\alpha$, switch to higher-order
Monitor everything — Not just residuals

Convergence Checklist

[ ] Residuals dropped 3+ orders of magnitude
[ ] Mass imbalance < $10^{-6}$ of inlet flux
[ ] Key quantities stabilized (< 0.1% change)
[ ] Physical plausibility verified
[ ] Results independent of further iterations

Steady vs. Transient

Steady-State Iterations

Goal: Find time-independent solution
Use SIMPLE, under-relaxation
Many outer iterations (hundreds to thousands)

Transient Time-Stepping

Goal: Resolve time evolution
Use PISO or dual time-stepping
Few iterations per time step (3-10)
Residuals must converge each step

Common Convergence Issues

Issue: Residuals Stuck at $10^{-3}$

Causes:

Mesh quality issues (high skewness)
Inconsistent boundary conditions
Physics not reaching steady state

Fixes:

Improve mesh quality
Check boundary condition values
Consider transient solution

Issue: One Variable Won't Converge

Often pressure or continuity:

Check outlet BC (needs pressure reference)
Verify mass balance at boundaries
Increase pressure under-relaxation iterations

Issue: Periodic Oscillations

Physical cause: Vortex shedding, unsteady separation Solution: Switch to transient simulation — the flow is genuinely unsteady!

Key Takeaways

Iterative solvers are essential for large CFD systems — direct methods don't scale
Gauss-Seidel converges faster than Jacobi but is less parallelizable
Multigrid achieves optimal $O(N)$ convergence by eliminating smooth errors on coarse grids
Residuals measure equation imbalance — target 3-4 orders reduction
Monitor key quantities — residuals alone don't guarantee accuracy
Under-relaxation stabilizes coupled solvers like SIMPLE
AMG enables multigrid on unstructured meshes automatically

What's Next

With the numerical machinery in place, we turn to one of CFD's greatest challenges: Turbulence Modeling. The next lesson covers Reynolds averaging, RANS equations, and popular turbulence models like k-ε and k-ω SST.