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.

Always run at least three meshes. Two meshes give you a ratio but no confidence interval; three or more let you confirm monotonic convergence and estimate the asymptotic range.

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:

  1. Select three geometrically similar meshes with refinement ratio r ≥ 1.3.
  2. Compute the observed order p from the fine/medium/coarse triplet.
  3. Calculate the relative error e = |(f_medium − f_fine) / f_fine|.
  4. 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:

  1. Pick smooth manufactured fields for all primitives (velocity, pressure, temperature, species).
  2. Symbolically differentiate and form the residual S = L(u_m) − R(u_m), where L is the differential operator and R is the source.
  3. Implement S as a source term in the solver, either analytically or via a symbolic library.
  4. Apply manufactured Dirichlet or Neumann boundary conditions derived from u_m.
  5. Run on a sequence of refined meshes and compute the L2 or L∞ error norm against u_m.
  6. Fit log(error) versus log(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.

MMS verifies the equation, not the physics. A manufactured solution that satisfies the discrete equations perfectly still tells you nothing about whether the turbulence model is appropriate for your application. MMS is a code-verification tool, not a validation tool.

Choosing the Right Tool for the Job

TechniqueWhat It VerifiesTypical OutputWhen to Use
Order-of-accuracy testFormal order of the discretizationObserved p vs. formal pAfter any scheme change or refactor
GCIDiscretization uncertainty on a production meshPercentage band on QoIReporting results for a specific case
MMSCorrect implementation of all PDE termsError norm slope, term-by-term coverageNew 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.

Code Verification Workflow: MMS → Order Test → GCI Method of Manufactured Solutions 1. Choose u_m(x,y,t) 2. Compute residual S = L(u_m) − R(u_m) 3. Add S as source term 4. Apply BCs from u_m 5. Refine mesh, compute error norm → Exact solution known by construction Order-of-Accuracy Test p = ln[(f_c − f_m)/(f_m − f_f)] / ln(r) Run 3+ meshes, r ≥ 1.3 Fit log(error) vs log(h) Compare p_obs to p_formal → Slope should match scheme order Grid Convergence Index GCI = F_s · e / (r^p − 1) F_s = 1.25 (3 meshes) e = |(f_m − f_f)/f_f| Report as % uncertainty band → Discretization uncertainty on QoI Convergence Plot (log-log) log(h) — mesh size log(error) slope = 2 (ideal) observed p ≈ 1.95 Failure Symptoms & Diagnosis • p collapses to 0.5 → mesh too coarse or solution non-smooth • Error norm plateaus → inconsistent BC or missing source term • p oscillates → non-asymptotic range, refine further • GCI > tolerance → refine mesh or improve quality • p matches formal order → code verified for that term MMS verifies the equation; GCI quantifies discretization uncertainty; order test confirms the scheme.