90108a6Command for testing against a reference implementationOwen Melia 1# Reference test case: Schnakenberg on a triaxial ellipsoid
3For validating this repo's Laplace-Beltrami surface correction against an
4independently-implemented reference solver, comparing final state in
5spherical-harmonic (SH) coefficient space.
7## Geometry
9A triaxial ellipsoid, `gx = ax·sinθ·cosφ`, `gy = ay·sinθ·sinφ`,
10`gz = az·cosθ` (`geometries/ellipsoid.m`; defaults `ax=1.5, ay=1.0, az=0.6`).
12**Important**: the solver does not run on this analytic surface — it
13band-limits it first, analysing `(gx, gy, gz)` into SH coefficients truncated
14at degree `lmax` and re-synthesizing before use (`src/geom/geometry.ts`). To
15remove geometry-representation error as a confound, a reference file
16supplies these coefficients directly as `/geometry/Gx`, `/geometry/Gy`,
17`/geometry/Gz`. **The reference solver must reconstruct its surface (and
18induced metric) by synthesizing these coefficients, not by evaluating the
19analytic formula above.**
21## Equations
23Schnakenberg reaction-diffusion with the true surface Laplace-Beltrami
24operator `Δ_g` (`models/schnakenberg.m`):
26```
27du/dt = D1·Δ_g(u) + a - u + u²v
28dv/dt = D2·Δ_g(v) + b - u²v
29```
31This repo's internal discretization (`niter` Richardson-iteration count for
32its own `Δ_g` approximation, `dt` for its IMEX-Euler timestep) is not part of
33the equations being tested — the reference solver may use any consistent
34method for `Δ_g` and any timestep. It only needs to reach the same physical
35end time `T = steps · dt`.
37## Initial condition
39Loaded from the reference file's `/initial/U` / `/initial/V` (t=0,
40immediately after seeding, before any step), not regenerated — avoids
41needing to reimplement this repo's PRNG (`src/mgpu/noise.ts`) to get a
42matching initial condition.
44## Output convention (must match exactly, from `src/sht/layout.ts`)
46- Orthonormal spherical harmonics **including Condon-Shortley phase**.
47- Real field ⇒ complex coefficients stored for `m ≥ 0` only:
48 `Q_{l,-m} = (-1)^m · conj(Q_lm)`; `m=0` coefficients have zero imaginary part.
49- **m-major ordering**: for `m = 0..lmax`, for `l = m..lmax`.
50 `index(l,m) = m·(lmax+1) − m·(m−1)/2 + (l−m)`.
51- Flat array, length `2·nlm` with `nlm = (lmax+1)(lmax+2)/2`, `[re,im]`
52 interleaved per coefficient (`qlm[2·index(l,m)]`, `qlm[2·index(l,m)+1]`).
54The reference solver must project its final `u`, `v` onto this same
55convention/truncation and report flat `2·nlm` arrays, to diff directly
56against the reference file's `/final/U` / `/final/V`.
58## HDF5 file layout
60Each reference file is one `.h5` file per run (written with
61[h5wasm](https://github.com/usnistgov/h5wasm); readable from Python with
62`h5py.File(path, "r")`). Coefficient datasets are `float32`, each of length
63`2·nlm` in the convention above. Metadata is stored as attributes, grouped by
64what it describes rather than as a single flat namespace:
66```
67/ (attrs: command, model, species)
68├─ backend/ (attrs: adapter, runtime, precision)
69├─ spec/ (attrs: preset, geometry, lmax, seed, steps, warmup, niter)
70│ ├─ params/ (attrs: the model's own params, e.g. a, b, D1, D2, dt)
71│ └─ geometry_params/ (attrs: the geometry's own params, e.g. ax, ay, az)
72├─ grid/ (attrs: lmax, mmax, nlat, nphi, nlm)
73├─ geometry/
74│ ├─ Gx dataset, float32[2·nlm]
75│ ├─ Gy dataset, float32[2·nlm]
76│ └─ Gz dataset, float32[2·nlm]
77├─ initial/ one dataset per species (e.g. U, V), float32[2·nlm] each
78└─ final/ one dataset per species (e.g. U, V), float32[2·nlm] each
79```
81`species` (root attribute) names which datasets live under `initial/` and
82`final/` — `["U", "V"]` for Schnakenberg. `command` is the equivalent
83`npm run bench --` invocation, for reproducing the run exactly.
85## Parameters
87| name | meaning | default |
88|---|---|---|
89| `a`, `b` | Schnakenberg kinetics | 0.1, 0.9 |
90| `D1`, `D2` | diffusion coefficients | 4e-4, 8e-3 |
91| `ax`, `ay`, `az` | ellipsoid semi-axes | 1.5, 1.0, 0.6 |
92| `lmax` | SH truncation degree | 63 |
93| `T = steps·dt` | physical end time | e.g. 2000·0.05 = 100 |
94| `seed` | provenance only — IC supplied as coefficients | 1 |
96## Checking a run against a reference file
98`npm run ref -- --in <file>` loads a reference file, runs this
99repo's own solver from its exact initial condition to the same physical end
100time, and reports the relative-L2 error of the resulting state against the
101file's final state (and, as a sanity check, of the regenerated geometry
102against the file's own geometry coefficients — this should be ~0 unless
103geometry construction itself has changed). `--niter <n>` overrides the
104solve's own iteration count for the surface correction, independent of what
105the reference file was generated with — useful for seeing how much that
106correction term actually matters for a given run. `--tolerance <n>` turns
107the check into a pass/fail (nonzero exit code on failure), for use in CI.
109## Caveat
111This repo runs fp32 on GPU; expect ~1e-4–1e-6 relative floating-point noise
112on top of any genuine numerical-method disagreement between solvers.