3. Solution verification in practice

Solution verification answers a deceptively simple question: how much numerical error is in my computed answer? It is distinct from code verification (does the solver implement the equations correctly?) and from validation (does the model represent reality?). In practice, solution verification focuses on two dominant error sources: discretization error, arising from representing continuous PDEs on a finite mesh, and iterative error, arising from stopping the linear/nonlinear solver before the algebraic system is fully converged. This chapter shows how to estimate both, how to use Richardson extrapolation to obtain a grid-converged value and observed order of accuracy, and how to set defensible iterative convergence criteria.

3.1 The Error Budget and the Role of Reference Values

For a scalar quantity of interest (QoI) such as drag coefficient, peak stress, or outlet temperature, the computed value on mesh h can be decomposed as:

φ_h = φ_exact + E_discretization(h) + E_iterative + E_roundoff + E_other

Solution verification seeks to bound E_discretization and confirm that E_iterative is negligible relative to it. A useful rule of thumb: iterative error should be at least one to two orders of magnitude smaller than the discretization error you are trying to measure. If the solver is only loosely converged, Richardson extrapolation will fit noise, not physics.

3.2 Iterative Convergence Criteria

Iterative convergence is judged by monitoring residuals and, critically, the QoI itself. Residuals alone are insufficient because they are scaled by the equations and mesh; a "small" residual can still correspond to a drifting QoI.

• Scaled residuals: Require the L2 norm of each equation's residual to fall below a tolerance (e.g., 1e-5 for pressure-based CFD, 1e-8 for structural solvers).
• QoI stabilization: Track the QoI over iterations; require its change over the last N iterations to be below a fraction of the expected discretization error (e.g., < 0.1%).
• Imbalance checks: For conservation laws, verify global mass, momentum, and energy imbalances are below ~0.1% of the flux magnitude.

3.3 Discretization Error Estimation via Grid Refinement

The standard approach is a systematic grid refinement study. Generate at least three meshes with a constant refinement ratio r = h_coarse/h_fine, typically r = √2 (2D) or r = 2^1/3 (3D) for unstructured meshes, or r = 2 for structured meshes. Compute the QoI on each mesh with tight iterative convergence.

Let the QoI values be φ₁ (fine), φ₂ (medium), φ₃ (coarse). The observed order of accuracy p is obtained from:

p = ln[(φ_3 - φ_2)/(φ_2 - φ_1)] / ln(r)

If the solver is second-order accurate and the meshes are in the asymptotic range, p should approach 2. A value far from the formal order (e.g., p < 1 or p > 3) signals that the meshes are too coarse, the refinement is non-uniform, or the solution is not smooth (shocks, singularities).

3.4 Richardson Extrapolation

Richardson extrapolation uses the QoI values and the observed order to estimate the continuum value φ_exact and the discretization error on the finest mesh:

φ_extrapolated = φ_1 + (φ_1 - φ_2)/(r^p - 1)
E_fine ≈ |φ_1 - φ_extrapolated|

The grid convergence index (GCI) provides a conservative uncertainty band:

GCI_fine = Fs * |(φ_1 - φ_2)/φ_1| / (r^p - 1)

where F_s is a safety factor (1.25 when three or more meshes confirm the asymptotic range; 3.0 when only two meshes are available). Report the QoI as φ₁ ± GCI_fine.

3.5 Worked Example: Drag Coefficient on Three Meshes

Consider a 2D airfoil simulation at Re = 6.5 × 10⁶, M = 0.15. Three unstructured meshes are generated with refinement ratio r = √2. The computed drag coefficients are:

Observed order: p = ln[(331.2 − 318.7)/(318.7 − 312.4)] / ln(√2) = ln(12.5/6.3)/0.3466 ≈ 1.98. This is close to the formal second order, confirming the asymptotic range.

Extrapolated value: φ_ext = 312.4 + (312.4 − 318.7)/(2^1.98 − 1) ≈ 312.4 − 6.3/2.96 ≈ 310.3 counts. Discretization error on the fine mesh ≈ 2.1 counts (0.67%). GCI with F_s = 1.25: GCI ≈ 1.25 × 6.3/312.4 / 2.96 ≈ 0.85%. Report C_d = 312.4 ± 2.7 counts.

3.6 Practical Pitfalls and Troubleshooting

• Non-monotonic convergence: If φ₂ is not between φ₁ and φ₃, the asymptotic range is not reached. Refine all meshes or check for mesh quality issues (skewness, aspect ratio).
• Iterative error contamination: If residuals are only reduced to 1e-3, the QoI may still drift by more than the discretization error. Tighten convergence and repeat.
• Refinement ratio too small: r < 1.1 amplifies round-off and discretization noise in the order estimate. Use r ≥ 1.3 where possible.
• QoI selection: Local quantities (peak stress at a singularity) may not converge monotonically. Use integral quantities (total force, mass flow) for order verification.

3.7 Summary

Solution verification is a quantitative discipline, not a checkbox. Iterative convergence must be tight enough that discretization error dominates. A three-mesh refinement study with a constant ratio yields the observed order of accuracy; Richardson extrapolation then provides a continuum estimate and the GCI provides an uncertainty band. When the observed order matches the formal order and the GCI is acceptably small, the solution is verified for the chosen QoI and can proceed to validation against experimental data.