Critical Reading of Simulations

Physical validation criteria, sensitivity analysis, conserved quantities checks, and numerical artifact detection.

1. The Verification and Validation (V&V) Paradigm

In modern computational science and engineering, a simulation is never an empirical proof; it is the numerical solution of an idealized mathematical model. Establishing the scientific integrity and quantitative credibility of computational physics demands strict adherence to the formal framework of Verification and Validation (V&V):

A. Numerical Verification: "Are we solving the equations right?"

Verification is a purely mathematical and computational assessment. It confirms that software algorithms solve the underlying continuum equations with theoretical accuracy, free from programming bugs, spatiotemporal discretization errors, and catastrophic floating-point cancellation. It encompasses:

1. Code Verification: Demonstrating theoretical asymptotic rates of convergence using the Method of Manufactured Solutions (MMS) or known analytical solutions.
2. Solution Verification: Quantifying the numerical discretization uncertainty of a specific simulation run with respect to the continuum mathematical boundary value problem.

B. Physical Validation: "Are we solving the right equations?"

Validation is an empirical physical inquiry. It quantifies how accurately the governing equations and constitutive models (e.g., incompressibility, ideal gas assumption, Hookean elasticity, point charges) represent true experimental reality within a designated operational envelope.

$\text{Physical Reality (Experiment)} \xrightarrow{\text{Validation}} \text{Mathematical Model} \xrightarrow{\text{Verification}} \text{Numerical Simulation}$

2. Systematic Grid Convergence and Roache's GCI Methodology

Proving that a numerical prediction is independent of mesh resolution requires a systematic grid refinement study across at least three successively refined grids with constant grid refinement ratio $r = h_{\text{coarse}} / h_{\text{fine}} > 1.3$.

Let $f_1, f_2, f_3$ be the discrete predictions of a critical physical quantity (e.g., total energy, drag coefficient, peak acoustic pressure) on the fine ($h_1$), medium ($h_2 = r h_1$), and coarse ($h_3 = r^2 h_1$) grids.

Apparent Order of Convergence ($p$):

$p = \frac{\ln\left(\frac{f_3 - f_2}{f_2 - f_1}\right)}{\ln(r)}$

• If calculated $p$ matches the formal Taylor order of the numerical scheme (e.g., $p \approx 2$ for central finite differences), the solution lies within the asymptotic grid convergence regime.
• If $p < 1$ or fluctuates irregularly, the mesh remains too coarse to resolve physical boundary layers or sharp gradients.

Richardson Extrapolation:

Richardson extrapolation estimates the continuous zero-grid-spacing asymptotic limit ($h \to 0$): $f_{\text{extrap}} = f_1 + \frac{f_1 - f_2}{r^p - 1}$

Grid Convergence Index (GCI):

Endorsed by ASME and AIAA standards, the GCI provides a 95% confidence numerical uncertainty estimate with safety factor $F_s = 1.25$ for three-grid studies:

$\text{GCI}_{12} = \frac{F_s \, |e_a|}{r^p - 1} \times 100\%, \quad e_a = \frac{f_1 - f_2}{f_1}$

3. Physical Invariants and Conservation Laws as Integrity Metrics

A simulation may converge smoothly while violating fundamental physical conservation laws due to asymmetric discretization or improper flux balances:

1. Energy Conservation: In unforced conservative systems, relative mechanical energy drift $\Delta E(t) / E_0$ must remain bounded or scale down strictly as $\mathcal{O}(h^p)$.
2. Linear/Angular Momentum and Galilean Invariance: Internal interactions must never accelerate the system center of mass.
3. Quantum Unitarity and Norm Preservation: In quantum simulations, total probability must satisfy $\int |\psi(\vec{x}, t)|^2 d^3x \equiv 1$.
4. Mass and Flux Conservation: Net boundary fluxes must balance volume accumulation rates $\oint \vec{J} \cdot d\vec{A} = -\frac{d}{dt} \int \rho \, dV$.

4. Identifying and Diagnosing Numerical Artifacts

A critical scientist must distinguish physical phenomena from computational artifacts:

5. Uncertainty Quantification (UQ): Aleatory vs Epistemic

Critical assessment separates uncertainty into two distinct categories:

• Aleatory Uncertainty (Irreducible): Inherent physical randomness (thermal fluctuations, turbulent eddies, quantum shot noise). Modeled via probability distributions and Monte Carlo ensembles.
• Epistemic Uncertainty (Reducible): Gaps in physical knowledge (uncertain material conductivity, unmeasured boundary roughness, approximate reaction rates). Mitigated through targeted empirical measurement.

Audit Protocol and Critical Reading of Physical Simulations

Systematic 5-Phase Methodology

To evaluate computational simulation outcomes with scientific rigor, follow this standardized five-phase audit:

Phase 1: Conceptual Model and Validity Envelope Audit

1. Identify all foundational theoretical approximations: Is turbulence active where laminar equations were applied? Are small-strain assumptions valid under finite deformations? Are relativistic or quantum corrections negligible?
2. Check characteristic dimensionless parameters (Reynolds $\text{Re}$, Mach $\text{Ma}$, Knudsen $\text{Kn}$, Péclet $\text{Pe}$) to ensure the physical regime matches the domain of applicability of the chosen equations.

Phase 2: Spatiotemporal Mesh Independence Verification

1. Never accept simulation results based on a single grid or a single timestep.
2. Mandate a systematic three-grid refinement study ($h_1 < h_2 < h_3$ with ratio $r \ge 1.3$).
3. Compute the apparent order of convergence $p$ and Roache's GCI:
   • If $\text{GCI}_{12} < 1\% - 2\%$, discretization error is under rigorous control.
   • If $\text{GCI}_{12} > 10\%$, quantitative conclusions are speculative and mathematically unresolved.

Phase 3: Conservation Laws and Invariant Telemetry

1. Plot time series of fundamental invariant drifts: $\Delta E(t) = \frac{|E(t) - E(0)|}{|E(0)|}, \quad \Delta M(t) = \frac{|M(t) - M(0)|}{|M(0)|}$
2. Verify global boundary flux balances: integral influx minus outflux must equal the internal rate of storage change within machine precision.

Phase 4: Active Searching for Numerical Artifacts

1. Symmetry Sanity Checks: If the geometry and boundary conditions are symmetric, does the solution exhibit spurious asymmetry driven by grid alignment or round-off noise?
2. Artificial Boundary Reflections: Check whether outgoing shockwaves or acoustic pulses reflect unphysically from truncation boundaries (indicating missing PML or non-reflecting BCs).
3. Smoothness and Monotonicity Checks: Inspect for high-frequency checkerboard patterns or Gibbs overshoots around sharp interfaces.

Phase 5: Cross-Validation with Benchmarks and Experimental Data

1. Validate against canonical benchmark problems admitting exact closed-form solutions (Couette/Poiseuille flow, D'Alembert pulses, harmonic oscillator).
2. Perform quantitative comparison against empirical laboratory measurements, superimposing experimental measurement error bars with numerical GCI uncertainty bands.

Worked Examples: Quantitative Simulation Auditing

Example 1: Grid Convergence Index (GCI) Calculation for Lift Coefficient

An aerodynamics team computes the lift coefficient $C_L$ of an airfoil using three systematically refined grids with constant refinement ratio $r = \sqrt{2} \approx 1.414$:

• Grid 3 (coarse, $N_3 = 25000$ cells): $C_{L,3} = 1.0820$
• Grid 2 (medium, $N_2 = 50000$ cells): $C_{L,2} = 1.1240$
• Grid 1 (fine, $N_1 = 100000$ cells): $C_{L,1} = 1.1350$

1. Relative Apparent Discretization Differences:

$e_{23} = \frac{C_{L,2} - C_{L,3}}{C_{L,2}} = \frac{1.1240 - 1.0820}{1.1240} = 0.03737 \, (3.74\%)$ $e_{12} = \frac{C_{L,1} - C_{L,2}}{C_{L,1}} = \frac{1.1350 - 1.1240}{1.1350} = 0.00969 \, (0.97\%)$

2. Apparent Convergence Order ($p$):

$p = \frac{\ln\left(\frac{C_{L,3} - C_{L,2}}{C_{L,2} - C_{L,1}}\right)}{\ln(r)} = \frac{\ln(0.0420 / 0.0110)}{\ln(1.4142)} = \frac{\ln(3.81818)}{0.34657} = \frac{1.33977}{0.34657} \approx 3.866$

Because $p > 0$ and values monotonic ($f_1 > f_2 > f_3$), the sequence converges monotonically toward an asymptotic continuum value.

3. Richardson Continuum Extrapolation ($h \to 0$):

$C_{L,\text{extrap}} = C_{L,1} + \frac{C_{L,1} - C_{L,2}}{r^p - 1} = 1.1350 + \frac{0.0110}{3.818 - 1} = 1.1350 + 0.0039 = 1.1389$

4. Fine Grid GCI with Safety Factor $F_s = 1.25$:

$\text{GCI}_{12} = \frac{F_s \, |e_{12}|}{r^p - 1} = \frac{1.25 \times 0.00969}{2.818} = 0.00430 \implies 0.43\%$

Scientific Audit Verdict: The fine-grid prediction is $C_L = 1.1350 \pm 0.0049$ with numerical discretization uncertainty bounded at $0.43\%$. The spatial mesh is rigorously verified.

Example 2: Distinguishing Physical Friction from Numerical Dissipation

An ideal undamped harmonic oscillator ($m = 1\text{ kg}$, $k = 100\text{ N/m}$, $\omega = 10\text{ rad/s}$, $\gamma = 0$) begins with mechanical energy $E(0) = 50.0\text{ J}$.

After $t = 10\text{ s}$ of simulation (approx 16 natural cycles), telemetry outputs:

• Solver A: $E(10) = 42.1\text{ J}$
• Solver B: $E(10) = 50.000003\text{ J}$
• Solver C: $E(10) = 118.5\text{ J}$

Telemetry Assessment:

1. Solver A (Backward Euler / Excessive Numerical Dissipation): Lost $7.9\text{ J}$ ($15.8\%$ of total energy) with zero physical drag in the model. The algorithm introduces artificial numerical damping that mimics physical viscosity. Audit Verdict: Reject for conservative dynamical modeling.
2. Solver B (Velocity-Verlet or High-Order RK4): Total energy is conserved within a relative drift of $6 \times 10^{-8}$. Audit Verdict: Physically valid and mathematically verified.
3. Solver C (Forward Euler): Total energy grew from $50\text{ J}$ to $118.5\text{ J}$ ($+137\%$). The algorithm continuously pumps unphysical work into the system. Audit Verdict: Spurious instability; invalidate run.

Common Pitfalls in Assessing and Auditing Simulations

1. Conflating Visual Realism with Scientific Truth

• Pitfall: Assuming a 3D simulation with photorealistic graphics and smooth colormaps is physically reliable.
• Consequence: Models hiding severe discretization errors, unphysical energy dissipation, or thermodynamic violations are erroneously accepted.
• Remedy: Mandate rigorous convergence plots, energy invariant tracking, and flux conservation tables before evaluating visual renders.

2. Engaging in Unphysical Parameter Overfitting

• Pitfall: Arbitrarily tweaking unmeasured simulation knobs (e.g., eddy viscosity, contact damping) solely to force agreement with experimental points.
• Consequence: The model loses predictive utility, collapsing immediately under altered operating boundary conditions.
• Remedy: Strictly isolate parameter calibration from true independent experimental validation under unseen configurations.

3. Trusting Single-Grid Simulation Results

• Pitfall: Accepting numerical conclusions produced on a single mesh resolution without grid independence testing.
• Consequence: Inability to determine whether observed structures are real physical features or discrete mesh-induced artifacts.
• Remedy: Always perform systematic grid studies with at least three refinement levels and calculate the formal GCI uncertainty.

4. Conflating Solver Iterative Convergence with Physical Discretization Accuracy

• Pitfall: Believing that because algebraic residuals dropped below $10^{-6}$, the simulation has 6-digit physical accuracy.
• Consequence: Algebraic solver residuals reflect matrix solution quality, not the discretization error between the discrete grid and continuum reality.
• Remedy: Quantify truncation and discretization error through systematic spatial and temporal grid halving.

5. Overlooking Reflections from Artificial Domain Boundaries

• Pitfall: Truncating computational domain boundaries too close to the region of interest without non-reflecting or perfectly matched layers (PML).
• Consequence: Outgoing acoustic, electromagnetic, or hydrodynamic waves bounce off artificial limits, contaminating interior physics.
• Remedy: Place artificial boundaries far from the interest zone and implement non-reflecting absorbing boundary conditions.

Physical Applications of Critical Simulation Auditing

1. Aerospace Aerodynamic Safety and Certification

In aerospace engineering (NASA, ESA, commercial aircraft certification), computational aerodynamic models require formal V&V metrics: $\text{GCI} < 1.0\% \text{ on lift } C_L \text{ and drag } C_D \text{ coefficients}$ Prior to piloted flight testing, hypersonic atmospheric re-entry simulations must demonstrate that peak heat shield temperatures contain no numerical discretization artifacts that could underpredict ablative wear.

2. Global Climate Modeling and Meteorology

In General Circulation Models (GCMs) predicting decadal global climate trends:

• Critical assessment audits climate sensitivity under doubled atmospheric $\text{CO}_2$.
• Methodological V&V isolates intrinsic unforced chaotic climate variability (aleatory uncertainty) from subgrid cloud parameterizations and ice albedo approximations (epistemic uncertainty).

3. Gravitational Wave Astrophysics and Black Hole Mergers

In analyzing gravitational wave chirp detections by LIGO/Virgo/KAGRA, template waveform catalogs for binary black hole inspirals are generated via numerical relativity: $G_{\mu\nu} = \frac{8\pi G}{c^4} T_{\mu\nu}$ Critical auditing verifies that energy and angular momentum carried away by gravitational radiation precisely match the reduction in the system's ADM mass and spin.

4. Nuclear Reactor Thermal-Hydraulics and Safety Licensing

In fission reactors and magnetic fusion experiments (ITER):

• Thermal-hydraulic system codes (e.g., RELAP5, TRACE) undergo strict validation against loss-of-coolant accident (LOCA) test facilities.
• Critical auditing proves that simulated peak cladding temperature limits are not masked or underpredicted by artificial numerical diffusion.

5. Cardiovascular Biomechanics and Cerebral Aneurysm Hemodynamics

In patient-specific computational fluid dynamics for stent deployment and cerebral aneurysm surgery:

• Solvers compute wall shear stress (WSS) along arterial endothelia.
• Grid convergence auditing is paramount: insufficient near-wall boundary layer resolution can underestimate peak wall shear stress by upwards of 50%, generating dangerously inaccurate clinical risk assessments of aneurysm rupture.

Practice Problems: Verification, Validation, and Simulation Auditing

Problem 1: Detecting Pre-Asymptotic Regime via Apparent Convergence Order

A researcher models thermal diffusion with an algorithm possessing a theoretical $\mathcal{O}(\Delta x^2)$ spatial stencil. Running three successively refined meshes ($r = 2$):

• Coarse grid ($h_3 = 0.4\text{ m}$): $T_3 = 320.0\text{ K}$
• Medium grid ($h_2 = 0.2\text{ m}$): $T_2 = 312.0\text{ K}$
• Fine grid ($h_1 = 0.1\text{ m}$): $T_1 = 309.0\text{ K}$

1. Calculate the apparent order of convergence $p$.
2. Contrast calculated $p$ with the theoretical prediction ($p = 2$). What does this reveal regarding the mesh resolution?

Step-by-Step Solution:

1. Inter-grid differences: $f_3 - f_2 = 320.0 - 312.0 = 8.0\text{ K}$ $f_2 - f_1 = 312.0 - 309.0 = 3.0\text{ K}$ With refinement factor $r = 2$: $p = \frac{\ln(8.0 / 3.0)}{\ln(2)} = \frac{\ln(2.6667)}{0.69315} = \frac{0.98083}{0.69315} \approx 1.415$
2. The formal stencil order is $p_{\text{theoretical}} = 2.0$, but the observed empirical order is only $p \approx 1.41$. This discrepancy proves the mesh is still outside the asymptotic convergence window. Local steep thermal gradients or boundary layers remain under-resolved at $h = 0.1\text{ m}$, preventing higher-order Taylor truncation terms from being genuinely negligible.

Problem 2: Mass Conservation Audit in CFD Duct Modeling

In a 2D incompressible fluid simulation ($\rho = 1000\text{ kg/m}^3$), measured mass flow rate at inlet is $\dot{m}_{\text{in}} = 15.42\text{ kg/s}$, while outlet rate is $\dot{m}_{\text{out}} = 14.88\text{ kg/s}$.

1. Determine the relative percentage mass conservation imbalance.
2. If commercial code documentation notes that "residuals dropped below $10^{-4}$," explain why iterative residual convergence does not ensure physical conservation.

Step-by-Step Solution:

1. Percentage mass imbalance: $\Delta \dot{m} = \frac{|15.42 - 14.88|}{15.42} = \frac{0.54}{15.42} \approx 0.0350 \implies 3.50\%$ The numerical solver artificially destroys $3.5\%$ of fluid mass inside the duct.
2. Solver residual convergence measures only algebraic discrepancy of the iterative matrix system ($||A x - b|| < 10^{-4}$), not integral physical flux conservation. If the formulation lacks strict cell-face conservative flux balances, or if the pressure-velocity coupling tolerance is lax, steady mass loss occurs despite nominal algebraic convergence.

Problem 3: Diagnosing Numerical Dispersion in Acoustic Wavefields

An acoustic simulation propagates a Gaussian pulse centered at $f = 1000\text{ Hz}$. After $20\text{ m}$ of travel, high-frequency oscillatory ripples develop behind the main packet that were absent initially.

1. Identify the numerical artifact responsible for these trailing ripples.
2. Outline two concrete corrective measures to eliminate the artifact without altering the physical acoustic properties.

Step-by-Step Solution:

1. The artifact is numerical dispersion. It arises because numerical phase velocity across discrete nodes is frequency-dependent when Courant number $C < 1$ or when resolution drops below 10-15 nodes per wavelength ($\lambda / \Delta x < 15$). High-frequency harmonics lag behind the primary wavepacket, forming trailing ripples.
2. Corrective interventions:
   • Refine the spatial grid to guarantee at least 20 nodes per shortest physical wavelength ($\Delta x \le \lambda_{\min} / 20$).
   • Operate closer to Courant number $C = 1.0$ (where 1D dispersion vanishes), or transition to higher-order compact or spectral difference stencils.