Simulation Validation

Validating a simulation you intend to show people

If you are visualising physics for an audience, correctness is the product. This is the pattern that makes that claim checkable instead of aspirational: every solver ships checks against an exact result, all of them report into one registry, and nothing renders until they pass.

Structure

physics/harness.py          CHECKS list + check_scalar / check_vectors / record
physics/checks_<area>.py    one module per solver, importing the harness
physics/validate.py         imports every check module, prints, exits non-zero

The harness lives in its own module so check modules and the runner can both import it without a cycle. Each check appends (name, ok, note) — the note is what makes a passing suite useful, because it tells you the margin, not just the verdict.

Compare vector fields by norm, never component-wise

err = max(norm(got - want, axis=-1) / maximum(norm(want, axis=-1), 1e-30))

A component that is analytically zero makes a component-wise relative error meaningless: it divides noise by noise and reports a relative error of 1e8 on a field that is perfectly correct. This wastes a lot of time if you meet it first as a mysterious failure.

Test limits by convergence, not by a point value

Some quantities are first-order results. They are not supposed to match at any finite parameter value, and a test claiming otherwise is lying about what was verified. Assert that the ratio approaches 1 as the small parameter shrinks, and put the sequence in the note:

PASS  perihelion advance -> 6piGM/(c^2 a(1-e^2))   measured/Einstein 1.0081 -> 1.0020 -> 1.0005 as c doubles
PASS  capacitor E -> sigma/eps0 as the gap closes  0.8623 -> 0.9229 -> 0.9555 at gap/size 0.20, 0.10, 0.05
PASS  ring sum -> closed-form drag                 k_sum/k_closed 0.9386 -> 1.0000 -> 1.0000

Each of those started as a failing point-comparison. Each time, the simulation was right and the test was wrong about what it should equal. Before loosening a tolerance, work out whether the quantity is a limit.

Errors falling by 4x per halving of the small parameter is the signature of a first-order expansion — good confirmation you understand the discrepancy rather than papering over it.

Choose the integrator for the invariant you care about

Measured: Boris held |v| to 1.8e-15 over 1500 steps; Verlet's orbital energy drifted 1.5e-14 over 60 orbits. RK4 on the same orbit loses energy monotonically and spirals in — at 30 fps that looks like physics and is numerics.

When you must use RK4, check a conserved quantity it does not know about: the Jacobi constant for the restricted three-body problem held to 3e-14.

Stiff ODEs: RK4 explodes, it does not wobble

RK4 on dv/dt = g - v/tau is unstable beyond about dt = 2.78 tau. Past that it does not degrade gracefully — it returns 1e245 and keeps going. A renderer will happily build a scene from that and produce an empty frame.

This is easy to hit when the settle time is milliseconds and "one frame" is tens of them. Make the integrator refuse:

tau = mass / k
if dt > 0.5 * tau:
    raise ValueError(f"dt={dt:.4g} exceeds half the settle time tau={tau:.4g}; "
                     "RK4 is unstable here. Integrate finely and resample.")

Integrate at tau/20 and np.interp onto your keyframe times.

Finite differences need a step scan

Checking something like div B = 0 numerically has no single correct h: too large and truncation error dominates, too small and cancellation does. Scan h over several decades and take the best, reporting it as such:

PASS  div B = 0 (dipole)   best |div B| / max|dBi/dxi| = 2.50e-08 over h scan

Normalise by the largest individual derivative, or you are testing the field's magnitude rather than its divergence.

Fit, do not average

Recovering an orbit radius by taking the mean of a partial orbit gives the wrong centre. Use a least-squares circle fit, which is robust to partial arcs:

A = column_stack([2*x, 2*y, ones(len(x))])
cx, cy, c = lstsq(A, x**2 + y**2)[0]
radius = sqrt(c + cx**2 + cy**2)

Also sample over whole periods where you can, and compare a chord to the field at the chord's midpoint rather than its start — comparing at the start introduces an O(h) bias that looks like a real error.

Exclude singularities from tracing

Every field line of a point dipole converges on its singularity, where a fixed-step integrator produces garbage — and the point-dipole model is invalid inside the magnet's body anyway. Pass the source bodies as (centre, radius) exclusion blobs and stop a line when it enters one.

Check the vectorised path against the naive one

Keep a slow, obviously-correct reference implementation (a per-source Python loop) and assert the fast einsum path matches it exactly on random points:

PASS  vectorised field == reference sum   rel err 0.00e+00 over 300 random points

This is the cheapest check in the suite and the one most likely to catch a silent indexing mistake during optimisation.

When a check fails, suspect the check

In this build, failing checks broke down roughly as:

• the test was wrong about the expected value (a limit, or an unphysical parameter choice) — most common
• a genuine sign error — a theta-hat built as z - cos(th) r instead of cos(th) r - z, which is -theta-hat. Watch for a relative error of exactly 2.0: that is got = -want, and it is a sign flip every time
• a real numerical problem — the stiff-ODE blow-up above

A relative error of exactly 2.0, or a ratio of exactly -1, should send you looking for a sign before anything else.

Parameters are part of the physics

A solver can be perfect and the scenario nonsense. Two from this build:

• A magnetic bottle with a gyroradius 51% of the coil radius is not adiabatic; the magnetic moment drifted 380% and the particle was not really trapped. At 6% it conserved to 3.3%.
• A copper-tube eddy brake with an implausible bore and a weak magnet gave a terminal velocity of 4.9 m/s. Real parameters — a half-inch pipe and a stack of three 10 mm magnets — gave 5.9 cm/s, matching the lecture demo.

Sanity-check derived quantities against something you can look up, and put the number in the check's note so it stays visible.

Derive rather than quote, then close the loop

Where a closed form is needed, derive it in the module docstring and then verify it against an independent numerical route that shares none of the algebra. The eddy-drag coefficient here was derived by integrating the ring kernel, then confirmed by a discrete ring sum converging on it to 1e-4 — and the ring kernel itself was checked against the already-validated Biot-Savart core.

A quoted constant you half-remember is a liability. A derivation with a convergence test behind it is not.