Verification and validation
How we know each layer is right. Tests live next to the code (unit), in crates/*/tests (integration), and in experiments/ (studies whose outputs are checked in CI at small size). Policy details follow the global testing skill.
1. Test pyramid
flowchart BT
U[Unit tests<br/>axes, constants, formulas,<br/>single recipes] --> P[Property tests<br/>conservation per pair,<br/>symmetry, positivity]
P --> O[Oracle and parity tests<br/>thesis Python, coag_py,<br/>scipy Radau, analytic solutions]
O --> S[Statistical tests<br/>estimator bias and variance<br/>over replicas]
S --> E[Experiment regression<br/>small-size versions of E1 to E16<br/>with stored reference outputs]
2. Oracles
| Oracle | Used for | How |
|---|---|---|
| Analytic solutions | Constant, sum, product kernels (pure coagulation); aggregation-fragmentation steady state if a closed form is available (Sekiyama et al. 2017) | Closed forms from Ohtsuki 1990, Wattis 2006, coag_py/mc_sampling/mc_run.py (with its Bessel-function asymptotics). |
Thesis Python (msc-thesis.py) | Disk, kinematics, outcomes, kernels, thesis-estimator runs | Fixture generator script runs the thesis code on fixed configs and stores arrays as .npy under tests/fixtures/thesis/ with the source commit hash. Generated once, committed. |
coag_py (Dullemond) | Kovetz and Olund kernel, cancellation handling, backward Euler, MC prototype | Same fixture approach, tests/fixtures/coag_py/. |
scipy solve_ivp(Radau) | Our Radau IIA, full trajectories | Python bridge (PyO3) or fixture files; equal rtol/atol, analytic Jacobian. |
| DustPy (Stammler and Birnstiel 2022) | Steady state of growth and fragmentation at one radius with transport off; the DustPy-like robustness setup | Project-local uv environment; matched settings recorded; steady-state distributions compared. |
| Conservation laws | Every kernel, every step | Built-in assertions in debug builds; property tests in CI. |
Fixtures are generated data from deterministic configs, not real observational data.
3. Invariants checked automatically
| Invariant | Tolerance | Where |
|---|---|---|
| Per-pair mass defect |Σ_k m_k X_kij| (mass gained equals mass lost) | ≤ 10 ε_mach (m_i + m_j) with cancellation handling | Property test over all pairs, every outcome model. |
| Per-pair number change for hit-and-stick | exactly -1 per collision | Property test. |
| Symmetry K_kij = K_kji (symmetrized form) | exact | Property test. |
| Global mass error δ_M(t) | ≤ 10⁻¹² (coag), documented bound (frag), at round-off growth rate √(steps) ε | Every integration test. |
| Positivity of N_k | no negative values beyond -10⁻³⁰ relative; negativity events counted and reported | Every integration test (replaces silent clipping, finding F7). |
| Sampled step mass conservation | same as full kernel, for every sample | Every sampled integration test. |
| Estimator unbiasedness | |mean − truth| < 3 SE over R replicas (R = 10³ in CI, 10⁵ in studies) | Statistical tests with fixed seeds. |
| Precision independence | f64 and double-double agree within the f64 error budget; the mass error floor drops with double-double | Extended-precision tests on small grids (ADR 0016). |
| Determinism | bit-identical output for identical config and seed, across thread counts | CI runs with 1 and N threads. |
Statistical tests use fixed seeds, so they are deterministic; their thresholds are chosen so a correct estimator fails with probability below 10⁻⁶ per test, documented next to each test.
4. Error metrics (replacing the thesis definitions)
Let ρ_k = m_k N_k and ρ = Σ ρ_k.
- Mass error: δ_M(t) = (ρ(t) - ρ(0)) / ρ(0). (Thesis chapter 3 definition.)
- Accuracy error vs reference (full kernel, same integrator, tight tolerances): mass-weighted L1 distance of the normalized distributions, e_acc(t) = ½ Σ_k | ρ_k(t)/ρ(t) − ρ^ref_k(t)/ρ^ref(t) |, which lies in [0, 1], is defined when bins are empty, and equals the fraction of mass “in the wrong bins”. Report also the log-space L2 error on bins above a mass floor for the shape at the small-mass tail.
- Time-lag diagnostic: mass-weighted mean log-mass ⟨ln m⟩(t) vs reference; separates a systematic lag (bias, F1) from noise.
- Steady-state noise floor: time-RMS of ‖dN/dt‖ (mass weighted) over the last decade in time, and the spread of ρ_k across replicas.
- Estimator diagnostics on fixed N: bias B_k, variance V_k, and the mass-weighted mean squared error MSE = Σ_k m_k² (B_k² + V_k) / ρ².
- Cost: recipe entries evaluated, kernel entries touched, Jacobian non-zeros, LU/Krylov flops, wall-clock (single thread, release build, fixed machine noted in the manifest).
- Work-normalized efficiency: η = 1 / (MSE × cost); the comparison metric for estimators. For full runs, work-precision diagrams (e_acc at t_end vs cost).
- Accuracy targets: e_acc = 1 % (headline), 10 % (quick look), 0.1 % where affordable; each at transient times and at steady state.
5. Integrator verification
| Check | Method |
|---|---|
| Order of convergence | Constant kernel with analytic solution; fixed step sequences; slope fit. Euler 1, DOPRI 5, backward Euler 1, linearized implicit Euler 1, Rosenbrock-W 2 or 3 (per chosen scheme), Radau IIA 5 (3 for stiff order reduction in some components, documented). |
| Stability on stiff problems | Full thesis kernel at large Δt: implicit methods must stay positive and bounded; explicit methods must fail beyond their CFL estimate. |
| Step-size control | Accepted/rejected step statistics match scipy within a factor of 2 on identical problems. |
| Jacobian | Analytic vs finite difference on random states, relative error < 10⁻⁷. |
| Steady state | Full model: RMS dN/dt decays to round-off level, matching thesis Fig. 103 behavior. |
6. Physics validation
| Check | Reference |
|---|---|
| Disk profiles | Thesis figures (temperature, sound speed, gas density 2×2); diskmodel.py. |
| Δv maps, St, C_ij | Thesis Figs. relative velocity (4 components, total), cross section, rate coefficient; Ormel and Cuzzi (2007) limiting cases. |
| Fragmentation probability | Thesis Fig. collision outcome probabilities; Stammler and Birnstiel (2022). |
| Steady state with fragmentation | Qualitative comparison with Birnstiel 2010 / DustPy at similar settings (peak size set by v_frag); optional direct DustPy run at one radius. |
| Bouncing | Dominik and Dullemond (2024) qualitative results (narrow distribution, small-grain removal). |
7. Continuous integration budget
just check must finish in under 5 minutes on the M1 Pro: unit, property, parity, and statistical tests at small sizes (N_m ≤ 50, R ≤ 10³). Full studies run through heavy and are not part of check.