2. Code verification fundamentals
Code verification answers a single, uncompromising question: does the software correctly solve the mathematical model it claims to solve? This is distinct from validation, which asks whether the model represents physical reality. A solver can be perfectly verified and still be physically wrong, or physically insightful yet numerically broken. Code verification isolates the numerical implementation from the physics, and it does so using three interlocking techniques: order-of-accuracy testing, the Grid Convergence Index (GCI), and the Method of Manufactured Solutions (MMS).
Why Discretization Error Must Be Measured, Not Assumed
Every finite-volume, finite-element, or finite-difference solver replaces continuous derivatives with discrete approximations. A second-order scheme carries a truncation error proportional to h², where h is a representative cell size or element length. The observed order of accuracy p is extracted from a refinement study by comparing solutions on systematically refined meshes:
p = ln( (f_coarse - f_medium) / (f_medium - f_fine) ) / ln(r)
Here r is the uniform refinement ratio (commonly 2 or √2), and f is a scalar quantity of interest — drag coefficient, peak temperature, centerline velocity, or an integral norm such as L2. If the code is implemented correctly and the asymptotic regime is reached, p should match the formal order of the scheme. A second-order upwind convection scheme should return p ≈ 2.0; a first-order scheme should return p ≈ 1.0. When p collapses to 0.5 or oscillates wildly, the mesh is too coarse, the solution is not smooth, or there is a coding error.
The Grid Convergence Index
Roache's Grid Convergence Index converts an observed order into a defensible uncertainty band. It is the standard reporting mechanism in ASME V&V 20 and AIAA guidance. The procedure is:
- Select three geometrically similar meshes with refinement ratio
r ≥ 1.3. - Compute the observed order
pfrom the fine/medium/coarse triplet. - Calculate the relative error
e = |(f_medium − f_fine) / f_fine|. - Apply the GCI formula with a safety factor
F_s(1.25 for three-mesh studies, 3.0 for two-mesh studies):
GCI_fine = F_s * e / (r^p - 1)
The result is a percentage. A GCI of 2.3% on lift coefficient means the fine-mesh value is within 2.3% of the asymptotic grid-independent solution at 95% confidence. If the GCI exceeds your engineering tolerance, refine further or improve mesh quality. Note that GCI quantifies discretization uncertainty only — it says nothing about turbulence model error or boundary condition error.
The Method of Manufactured Solutions
Order-of-accuracy testing needs an exact solution to compare against. For simple cases like Poiseuille flow you have one, but for a compressible reacting solver you do not. MMS solves this by manufacturing a solution. You choose an analytic function for every dependent variable — say u(x,y) = sin(πx)cos(πy) — substitute it into the governing PDE, and compute the residual analytically. That residual becomes a source term added to the code. The code then solves the modified PDE, and the manufactured function is the exact solution by construction.
The workflow is concrete:
- Pick smooth manufactured fields for all primitives (velocity, pressure, temperature, species).
- Symbolically differentiate and form the residual
S = L(u_m) − R(u_m), whereLis the differential operator andRis the source. - Implement
Sas a source term in the solver, either analytically or via a symbolic library. - Apply manufactured Dirichlet or Neumann boundary conditions derived from
u_m. - Run on a sequence of refined meshes and compute the
L2orL∞error norm againstu_m. - Fit
log(error)versuslog(h); the slope is the observed order.
MMS is powerful because it exercises every term in the discretization — convection, diffusion, source, reaction, time integration — without requiring a physically realizable flow. It catches sign errors, missing terms, wrong coefficients, and boundary condition bugs that a simple analytic case would miss.
Choosing the Right Tool for the Job
| Technique | What It Verifies | Typical Output | When to Use |
|---|---|---|---|
| Order-of-accuracy test | Formal order of the discretization | Observed p vs. formal p | After any scheme change or refactor |
| GCI | Discretization uncertainty on a production mesh | Percentage band on QoI | Reporting results for a specific case |
| MMS | Correct implementation of all PDE terms | Error norm slope, term-by-term coverage | New solver, new physics module, regression suite |
Practical Pitfalls
Three failure modes dominate real verification studies. First, non-asymptotic behavior: if the coarse mesh is too coarse, the observed order is meaningless. Always confirm that p stabilizes as h → 0. Second, non-smooth manufactured solutions: MMS requires derivatives of sufficient order to exist; a |x| term will destroy the convergence rate. Third, inconsistent boundary conditions: if the manufactured BC is applied at the wrong location or with the wrong sign, the error norm will plateau instead of converging.
A robust verification suite runs MMS on every build, checks that the observed order matches the formal order within a tolerance (typically ±0.1 for second-order schemes), and archives the GCI for each production mesh. This turns code verification from a one-time exercise into a continuous, auditable process — the foundation on which simulation credibility is built.