Testing

Physics Simulation Testing: Validation Strategies for FEA and CFD Solvers

HelpMeTest

HelpMeTest

— 7 min read
Physics Simulation Testing: Validation Strategies for FEA and CFD Solvers

Physics simulation software is tested differently from business software. There is rarely a "correct answer" in a database to compare against. Instead, verification and validation (V&V) relies on mathematical relationships: convergence rates, conservation laws, symmetry properties, and known analytical solutions. This guide covers the full stack of simulation testing strategies, with concrete examples in Python using FEniCS-style PDE solvers.

The V&V Framework

The simulation community distinguishes two categories of testing:

Verification — does the code correctly implement the mathematical model? Unit tests, convergence tests, and the method of manufactured solutions all verify the implementation.

Validation — does the mathematical model correctly represent the physical reality? This requires comparison against experimental data or high-fidelity reference simulations.

This guide focuses on verification, which is automatable. Validation requires experimental data and domain expertise.

Method of Manufactured Solutions (MMS)

MMS is the most rigorous technique for verifying PDE solvers. The idea:

  1. Choose a smooth, known solution u_exact
  2. Compute the forcing function f that makes u_exact satisfy your PDE
  3. Solve the PDE with that forcing function
  4. Compare the numerical solution against u_exact

This works for any PDE, any boundary condition, any discretization — without needing an analytical solution to a physically meaningful problem.

MMS Example: Poisson Equation

The Poisson equation: -∇²u = f in Ω, with Dirichlet boundary conditions.

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

def manufactured_poisson_1d(n: int) -> tuple[np.ndarray, np.ndarray, float]:
    """
    Solve -u'' = f on [0, 1] with u(0) = u(1) = 0.

    Manufactured solution: u_exact(x) = sin(pi * x)
    Forcing function: f(x) = pi^2 * sin(pi * x)

    Returns: (numerical solution, exact solution, L2 error)
    """
    dx = 1.0 / (n + 1)
    x = np.linspace(dx, 1.0 - dx, n)  # interior points

    # Exact solution and forcing
    u_exact = np.sin(np.pi * x)
    f = np.pi**2 * np.sin(np.pi * x)

    # Assemble FD stiffness matrix: tridiagonal
    diagonals = [
        -np.ones(n - 1),                    # subdiag
        2.0 * np.ones(n) / dx**2,           # main diag
        -np.ones(n - 1),                    # superdiag
    ]
    A = diags(diagonals, [-1, 0, 1], format='csr') / dx**2

    # Solve
    u_numerical = spsolve(A, f)

    # L2 error
    l2_error = np.sqrt(dx * np.sum((u_numerical - u_exact)**2))

    return u_numerical, u_exact, l2_error


def test_mms_poisson_1d_convergence():
    """
    Finite differences on Poisson should converge at O(h^2).
    Test that the convergence rate is between 1.8 and 2.2.
    """
    mesh_sizes = [20, 40, 80, 160, 320]
    errors = []

    for n in mesh_sizes:
        _, _, l2_error = manufactured_poisson_1d(n)
        errors.append(l2_error)

    # Compute convergence rates between consecutive mesh refinements
    rates = []
    for i in range(1, len(errors)):
        rate = np.log(errors[i-1] / errors[i]) / np.log(2.0)
        rates.append(rate)

    avg_rate = np.mean(rates[-3:])  # Use the last 3 refinements (asymptotic behavior)

    assert 1.8 < avg_rate < 2.2, (
        f"Expected O(h^2) convergence, got rate={avg_rate:.3f}. "
        f"Errors: {[f'{e:.2e}' for e in errors]}"
    )

MMS for 2D Heat Equation

def manufactured_heat_2d(nx: int, ny: int, dt: float, t_final: float):
    """
    Solve u_t = ∇²u + f on [0,1]^2 with Dirichlet BC.

    Manufactured solution: u(x, y, t) = sin(pi*x) * sin(pi*y) * exp(-2*pi^2*t)
    This satisfies u_t = ∇²u exactly (no forcing needed for this choice).
    """
    dx = 1.0 / (nx + 1)
    dy = 1.0 / (ny + 1)
    x = np.linspace(dx, 1.0 - dx, nx)
    y = np.linspace(dy, 1.0 - dy, ny)
    X, Y = np.meshgrid(x, y)

    def u_exact(t):
        return np.sin(np.pi * X) * np.sin(np.pi * Y) * np.exp(-2 * np.pi**2 * t)

    # Initialize
    u = u_exact(0.0)
    n_steps = int(t_final / dt)

    # Forward Euler (first-order in time, second-order in space)
    for _ in range(n_steps):
        laplacian = (
            (np.roll(u, -1, axis=1) - 2*u + np.roll(u, 1, axis=1)) / dx**2 +
            (np.roll(u, -1, axis=0) - 2*u + np.roll(u, 1, axis=0)) / dy**2
        )
        # Zero boundary (roll introduces periodic artifacts — fix boundary)
        laplacian[:, 0] = laplacian[:, -1] = 0
        laplacian[0, :] = laplacian[-1, :] = 0
        u = u + dt * laplacian

    return u, u_exact(t_final)


def test_heat_equation_mms_spatial_convergence():
    """O(h^2) spatial convergence for heat equation (fixed small dt)."""
    dt = 1e-5
    t_final = 0.01
    mesh_refinements = [10, 20, 40, 80]
    errors = []

    for n in mesh_refinements:
        u_num, u_ex = manufactured_heat_2d(n, n, dt, t_final)
        l2_err = np.sqrt(np.mean((u_num - u_ex)**2))
        errors.append(l2_err)

    rates = [
        np.log(errors[i-1] / errors[i]) / np.log(2.0)
        for i in range(1, len(errors))
    ]
    avg_rate = np.mean(rates[-2:])

    assert avg_rate > 1.7, f"Spatial convergence rate too low: {avg_rate:.3f}"

Convergence Order Testing

Every numerical method has a theoretical convergence order. Testing it is non-negotiable:

Method Expected order Convergence test
Finite differences (2nd order) O(h²) in space Halve h, error should quarter
Runge-Kutta 4 O(Δt⁴) in time Halve Δt, error should drop 16×
Finite elements (linear) O(h²) in L2 Richardson extrapolation
Spectral methods Exponential Error vs wavenumber plot

The convergence rate is computed by fitting a line to log(error) vs log(h):

def compute_convergence_rate(mesh_sizes: list[int], errors: list[float]) -> float:
    """
    Fit a line to log(h) vs log(error) to get the convergence rate.
    Returns the slope (negative means error decreases with refinement).
    """
    log_h = np.log([1.0 / n for n in mesh_sizes])
    log_e = np.log(errors)
    # Linear regression
    coeffs = np.polyfit(log_h, log_e, 1)
    return coeffs[0]  # slope = convergence order

The key detail: use the last few refinements for the rate estimate. Coarse meshes are often in the pre-asymptotic regime where the rate is not yet the theoretical value.

Conservation Law Validation

Physical simulations must conserve mass, energy, and momentum. These conservation tests are independent of knowing the exact solution — they verify that the solver does not introduce or destroy physical quantities.

def test_mass_conservation_incompressible_flow():
    """
    Incompressible Navier-Stokes: ∇·u = 0 must hold at every time step.
    Total mass in a closed domain must be constant.
    """
    # Synthetic velocity field: irrotational, divergence-free (e.g., potential flow)
    nx, ny = 50, 50
    dx = dy = 1.0 / nx
    x = np.linspace(0, 1, nx)
    y = np.linspace(0, 1, ny)
    X, Y = np.meshgrid(x, y)

    # Velocity field: u = (cos(pi*x)*sin(pi*y), -sin(pi*x)*cos(pi*y))
    # This is divergence-free analytically
    u = np.cos(np.pi * X) * np.sin(np.pi * Y)
    v = -np.sin(np.pi * X) * np.cos(np.pi * Y)

    # Compute divergence using central differences
    du_dx = (np.roll(u, -1, axis=1) - np.roll(u, 1, axis=1)) / (2 * dx)
    dv_dy = (np.roll(v, -1, axis=0) - np.roll(v, 1, axis=0)) / (2 * dy)
    divergence = du_dx + dv_dy

    # Interior divergence should be near machine epsilon
    interior_div = divergence[2:-2, 2:-2]
    max_div = np.max(np.abs(interior_div))

    assert max_div < 1e-10, (
        f"Velocity field is not divergence-free: max |∇·u| = {max_div:.2e}"
    )


def test_energy_conservation_wave_equation():
    """
    Wave equation u_tt = c^2 u_xx: total energy E = KE + PE must be conserved.
    E(t) / E(0) should stay within 1% for a non-dissipative scheme.
    """
    n = 200
    dx = 1.0 / n
    c = 1.0
    dt = 0.4 * dx / c  # CFL condition

    x = np.linspace(0, 1, n + 1)
    # Gaussian initial condition
    u0 = np.exp(-100 * (x - 0.5)**2)
    u_prev = u0.copy()
    u_curr = u0.copy()  # Zero initial velocity

    def total_energy(u_prev, u_curr):
        # Kinetic energy: (u_t)^2 ≈ ((u_curr - u_prev) / dt)^2
        u_t = (u_curr - u_prev) / dt
        KE = 0.5 * dx * np.sum(u_t**2)
        # Potential energy: c^2 * (u_x)^2
        u_x = (u_curr[1:] - u_curr[:-1]) / dx
        PE = 0.5 * c**2 * dx * np.sum(u_x**2)
        return KE + PE

    E0 = total_energy(u_prev, u_curr)
    n_steps = 100

    for _ in range(n_steps):
        # Leapfrog scheme
        u_next = 2 * u_curr - u_prev + (c * dt / dx)**2 * (
            np.roll(u_curr, -1) - 2 * u_curr + np.roll(u_curr, 1)
        )
        u_next[0] = u_next[-1] = 0  # Dirichlet BC
        u_prev = u_curr
        u_curr = u_next

    E_final = total_energy(u_prev, u_curr)
    energy_ratio = E_final / E0

    assert 0.99 < energy_ratio < 1.01, (
        f"Energy not conserved: E_final/E_0 = {energy_ratio:.6f}"
    )

Benchmark Problem Suites

The NASA Langley Research Center maintains a collection of turbulence model validation cases (turbmodels.larc.nasa.gov). For structural analysis, the NAFEMS benchmark suite provides reference solutions for FEA verification. Using these in your test suite:

import json
from pathlib import Path

class NAFEMSBenchmark:
    """
    NAFEMS LE1: Elliptic membrane under uniform pressure.
    Reference solution: sigma_yy at point D = 92.7 MPa.
    """
    REFERENCE_SIGMA_YY = 92.7e6  # Pa
    TOLERANCE = 0.02  # 2% tolerance on benchmark

    def run_solver(self, mesh_size: float) -> dict:
        """Run the FEA solver and return stress at reference point D."""
        # This would call your actual FEA solver
        # Returning mock data for illustration
        return {"sigma_yy_at_D": 92.7e6 * (1.0 + 0.01 * np.random.randn())}

def test_nafems_le1_elliptic_membrane():
    """NAFEMS LE1 benchmark: linear elastic elliptic membrane."""
    benchmark = NAFEMSBenchmark()
    result = benchmark.run_solver(mesh_size=0.1)

    sigma_computed = result["sigma_yy_at_D"]
    ref = NAFEMSBenchmark.REFERENCE_SIGMA_YY
    rel_error = abs(sigma_computed - ref) / ref

    assert rel_error < NAFEMSBenchmark.TOLERANCE, (
        f"NAFEMS LE1 failed: computed={sigma_computed/1e6:.2f} MPa, "
        f"reference={ref/1e6:.2f} MPa, error={rel_error:.2%}"
    )

Store benchmark results in a JSON file and compare against them in CI:

def test_regression_solver_output():
    """Regression test: compare against stored reference output."""
    ref_path = Path("tests/reference_data/poisson_n100.npy")
    _, u_exact, _ = manufactured_poisson_1d(n=100)
    u_num, _, _ = manufactured_poisson_1d(n=100)

    if not ref_path.exists():
        # First run: save the reference
        np.save(ref_path, u_num)
        pytest.skip("Reference data created — run again to compare")

    reference = np.load(ref_path)
    np.testing.assert_allclose(
        u_num, reference,
        rtol=1e-12,
        err_msg="Solver output changed from reference. Intentional change? Update reference."
    )

FEniCS Test Harness

FEniCS (now FEniCSx/dolfinx) provides a Python interface to finite element solvers. A minimal test harness:

# Requires: pip install fenics-dolfinx
import pytest

try:
    import dolfinx
    import ufl
    from mpi4py import MPI
    HAS_FENICS = True
except ImportError:
    HAS_FENICS = False

@pytest.mark.skipif(not HAS_FENICS, reason="FEniCSx not installed")
def test_fenics_poisson_mms():
    """MMS verification for FEniCSx Poisson solver."""
    from dolfinx import mesh, fem, default_scalar_type
    from dolfinx.fem.petsc import LinearProblem
    import ufl

    comm = MPI.COMM_WORLD
    domain = mesh.create_unit_square(comm, 32, 32)
    V = fem.functionspace(domain, ("Lagrange", 1))

    # Manufactured solution: u_exact = sin(pi*x)*sin(pi*y)
    x = ufl.SpatialCoordinate(domain)
    u_exact_expr = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
    f_expr = 2 * ufl.pi**2 * u_exact_expr

    # Boundary condition: u = u_exact on boundary
    u_exact_fn = fem.Function(V)
    u_exact_fn.interpolate(
        lambda x: np.sin(np.pi * x[0]) * np.sin(np.pi * x[1])
    )

    facets = mesh.locate_entities_boundary(
        domain, domain.topology.dim - 1,
        lambda x: np.ones(x.shape[1], dtype=bool)
    )
    dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, facets)
    bc = fem.dirichletbc(u_exact_fn, dofs)

    # Variational form
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    f = fem.Function(V)
    f.interpolate(lambda x: 2 * np.pi**2 * np.sin(np.pi * x[0]) * np.sin(np.pi * x[1]))

    a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
    L = f * v * ufl.dx

    problem = LinearProblem(a, L, bcs=[bc])
    uh = problem.solve()

    # Compute L2 error
    error_form = fem.form((uh - u_exact_expr)**2 * ufl.dx)
    l2_error = np.sqrt(comm.allreduce(fem.assemble_scalar(error_form), op=MPI.SUM))

    assert l2_error < 1e-3, f"FEniCSx Poisson L2 error too large: {l2_error:.2e}"

Summary

The testing stack for physics simulation verification:

  • MMS — the canonical technique for verifying any PDE solver; works for any equation, any geometry
  • Convergence order tests — confirm the method achieves its theoretical rate under mesh refinement
  • Conservation law tests — verify mass, energy, and momentum are not artificially created or destroyed
  • Benchmark suites — NAFEMS for FEA, NASA turbulence cases for CFD provide community-validated reference solutions
  • Regression tests — pin solver output to reference arrays; force explicit decisions when results change

The most important insight: in simulation testing, "correct" means mathematically consistent. A solver that conserves energy, achieves the right convergence order, and matches benchmark solutions is very likely correct — even without an exact analytical solution to compare against.

Read more

Start now free