Multigrid Method for CFD Convergence Acceleration

Simple iterative methods like Gauss-Seidel are very good at damping high-frequency error components — those that vary rapidly from cell to cell — but terrible at eliminating low-frequency (smooth, long-wavelength) errors. On a large mesh, errors with wavelengths comparable to the domain size converge so slowly that the solver effectively stalls. Multigrid exploits a key insight: smooth errors on a fine grid look like high-frequency errors on a coarser grid, where an iterative smoother can eliminate them efficiently. By cycling between fine and coarse grids, multigrid achieves $O(N)$ solution complexity — optimal convergence that cannot be beaten asymptotically.

Exactly. The fine-grid residual is "restricted" to a coarser grid, solved (or smoothed) there, and the correction is "prolongated" back to correct the fine-grid solution. Repeating this with multiple coarsening levels — fine → coarse → coarser → coarsest → back up — is what gives the V-cycle its characteristic shape. The coarsest level is solved exactly (it's small enough for a direct solver). The dramatic convergence rate improvement — from thousands of iterations with Gauss-Seidel to just 5–20 V-cycles for the same accuracy — makes multigrid indispensable for large-scale CFD.

For the $N$-point Poisson equation $A\mathbf{u} = \mathbf{f}$ (the pressure equation in incompressible CFD), the spectral radius of the Gauss-Seidel iteration matrix is approximately:

For $N = 100$ unknowns in each direction, $\rho \approx 0.9999$ — each iteration reduces the error by only 0.01%. To reduce the error by a factor of 10 requires $\log(10)/\log(1/\rho) \approx N^2/\pi^2 \cdot \ln(10) \approx 23,000$ iterations. Doubling the grid resolution makes convergence 4× slower. For a 1M cell CFD mesh ($N \approx 100^3 = 10^6$), purely iterative solvers would require astronomical iteration counts. Multigrid breaks this $O(N^2)$-per-digit-of-accuracy curse and achieves $O(N)$.

A single V-cycle for solving $A_h u_h = f_h$ on the fine grid proceeds as follows:

1. Pre-smoothing: Apply $\nu_1$ iterations of the smoother (Gauss-Seidel) on the fine grid to eliminate high-frequency errors.
2. Compute residual: $r_h = f_h - A_h u_h^{(current)}$
3. Restrict residual: $r_{2h} = I_h^{2h} r_h$ — project the residual to the coarse grid.
4. Coarse-grid solve: Solve $A_{2h} e_{2h} = r_{2h}$ (either exactly if the coarse grid is small enough, or recursively with another V-cycle).
5. Prolongate correction: $u_h \leftarrow u_h + I_{2h}^h e_{2h}$ — interpolate the coarse-grid correction back to the fine grid.
6. Post-smoothing: Apply $\nu_2$ more smoother iterations to remove high-frequency errors introduced by the prolongation.

Typical values: $\nu_1 = 1$, $\nu_2 = 2$ for standard CFD problems. The V-cycle costs $O(N)$ operations and reduces the error by a constant factor (typically 0.1–0.3) per cycle, independent of $N$. This gives true $O(N)$ solution complexity.

Two operators are needed:

• Restriction $I_h^{2h}$: Projects the fine-grid residual to the coarse grid. The simplest choice is injection (copy the value at every other grid point). Better accuracy comes from Full Weighting, which takes a weighted average of neighboring fine-grid points:

• Prolongation $I_{2h}^h$: Interpolates the coarse-grid correction back to the fine grid. Linear (bilinear/trilinear) interpolation is standard and is the adjoint of the Full Weighting restriction: $I_{2h}^h = (I_h^{2h})^T$.

For the Galerkin coarse-grid operator condition, the coarse-grid operator should satisfy $A_{2h} = I_h^{2h} A_h I_{2h}^h$. This automatically generates consistent coarse-grid problems and is the standard in algebraic multigrid.

The V-cycle does a single fine→coarse→fine pass; the W-cycle recurses twice at each coarser level, doing more work per cycle but converging faster per cycle:

For most CFD applications (incompressible SIMPLE loops), V-cycle is optimal. W-cycle shows advantage for highly stretched grids (aspect ratios >100:1) or strongly anisotropic problems. FMG is increasingly used to initialize transient simulations — the steady-state FMG solution provides a much better starting point than a uniform initial condition, reducing the unsteady "spin-up" phase.

That's where Algebraic Multigrid (AMG) comes in. AMG constructs coarse-level equations purely from the matrix sparsity pattern and coupling strengths — no geometric mesh hierarchy required. This makes it applicable to any FVM or FEM system on unstructured, polyhedral, or even meshless grids.

• Ruge-Stüben (Classical) AMG: Classifies unknowns into C-points (coarsened/kept) and F-points (removed/interpolated) based on strong coupling. Most common for CFD pressure equations. Implemented in Fluent's AMG solver and OpenFOAM's GAMG.
• Smoothed Aggregation AMG: Clusters neighboring unknowns into aggregates at coarser levels. Effective for structural mechanics and scalar transport problems. Used in TRILINOS ML and ParaSails libraries.
• AMGe (element-based): Uses element stiffness matrices for coarsening. Specialized for FEM systems — gives superior performance for structural dynamics problems.

The price of AMG versus geometric multigrid is higher setup cost — the coarsening and coarse-grid construction typically takes 20–30% of the total solver time, but amortizes well for multi-timestep transient problems.

In OpenFOAM, multigrid is configured in the fvSolution file. Here's a typical GAMG setup for the pressure equation with explanations:

p
{
    solver          GAMG;          // Geometric-Algebraic Multi-Grid
    smoother        GaussSeidel;   // Smoother: GaussSeidel, DILU, DIC
    nPreSweeps      0;             // Smoothing passes before coarsening
    nPostSweeps     2;             // Smoothing passes after correction
    cacheAgglomeration true;       // Cache coarse grid structure (transient)
    nCellsInCoarsestLevel 10;      // Min cells at coarsest level
    agglomerator    faceAreaPair;  // Coarsening algorithm
    mergeLevels     1;             // Merge levels per coarsening
    tolerance       1e-6;          // Absolute residual tolerance
    relTol          0.01;          // Relative tolerance (relax for inner iter)
}

GAMG (Geometric-Algebraic Multi-Grid) is OpenFOAM's primary AMG solver. For the pressure Poisson equation it typically runs 5–10× faster than ICCG (incomplete Cholesky conjugate gradient), making it the default choice for large cases. For velocity and scalar transport equations, PBICG or PBiCGStab with DILU preconditioning is often faster than GAMG due to the non-symmetric nature of those matrices.

Oscillating pressure residuals in Fluent's AMG solver usually indicate one of three things: (1) Aggressive coarsening with poor smoother matching — the coarse-grid corrections overcorrect the fine-grid solution. Try increasing nPostSweeps or switching the smoother from ILU to Gauss-Seidel. (2) Highly anisotropic mesh — cells with very high aspect ratios (long and thin boundary layer cells) degrade the quality of coarse-grid operators. Reduce the relaxation factor for pressure (from 0.3 to 0.2–0.25). (3) Incorrect AMG cycle type — Fluent's F-cycle (a variant of W-cycle) sometimes oscillates on problems with many near-zero eigenvalues. Switching to V-cycle or reducing the F-cycle frequency can help.

Here's a comparison of multigrid implementations across the main platforms:

Several exciting directions are active:

• GPU-AMG: NVIDIA's AmgX library brings AMG to GPU clusters. For pressure equations in large LES/DDES simulations (10M–100M cells), GPU-AMG delivers 10–50× speedup compared to CPU-AMG. Ansys Fluent 2024 includes direct AmgX integration.
• ML-based smoothers: Replacing Gauss-Seidel with a neural network smoother. Trained on representative problem instances, these learned smoothers adapt to local anisotropy and achieve better convergence rates for heterogeneous problems (e.g., multimaterial domains in CHT).
• p-multigrid: Instead of coarsening the mesh, coarsen the polynomial degree in high-order DG (discontinuous Galerkin) or spectral element methods. Especially effective for large-scale aeroelastic and turbomachinery simulations.
• Quantum multigrid (research): Combining quantum computing's matrix arithmetic with multigrid's hierarchical structure is a theoretical direction for exponentially large eigenvalue problems.

Systematic diagnostic steps:

1. Incorrect coarsening algorithm: For highly stretched meshes (boundary layers), standard face-area-pair agglomeration may not create effective coarse grids. Try anisotropic or paired coarsening. In Fluent, switch the AMG coarsening from IBJA to PILUT for highly anisotropic grids.
2. Poor smoother choice: For non-symmetric systems (momentum equations with convection), Gauss-Seidel underperforms. Use DILU (diagonal incomplete LU) for velocity equations.
3. Too few coarsening levels: If the coarsest grid still has many unknowns, low-frequency errors are not fully eliminated. Verify that GAMG actually coarsens to fewer than 10 cells at the coarsest level (check OpenFOAM log output).
4. Anisotropic grid: Highly stretched cells (BL mesh with aspect ratio >100) degrade multigrid performance significantly. Consider using ILU-preconditioned BICGSTAB for the viscous layers rather than AMG.