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
66ae13eUpdated command now tracks L_infty error too.Owen Melia 100time, and reports the relative-L2 and relative-L-infinity (max-norm) error of
101the resulting state against the file's final state (and, as a sanity check,
102of the regenerated geometry against the file's own geometry coefficients —
103this should be ~0 unless geometry construction itself has changed). `--niter
104<n>` overrides the solve's own iteration count for the surface correction,
105independent of what the reference file was generated with — useful for seeing
106how much that correction term actually matters for a given run. `--tolerance
107<n>` and `--tolerance-linf <n>` each independently turn their metric into a
108pass/fail (nonzero exit code on failure), for use in CI.
c90d0e2Check reference files in the browser's compare modeJeremy Magland 110The browser demo runs the same check visually: **Compare → Reference file…**
111loads a reference file into the convergence study, seeds every variant from
112its exact initial state, runs them side by side to its end time, and shows
113its final state as one extra static row — on the file's own surface, with
114each variant's relative-L2 distance to it updating live. Both readers share
115one parser (`src/compare/referenceCase.ts`), so the layout above is
116interpreted identically on the CLI and in the page.
118Note the files record only the two endpoint states (`initial/`, `final/`) —
119no intermediate snapshots — so the comparison is meaningful at the end time;
120the live Δ before that reads as "distance still to the final state".
124This repo runs fp32 on GPU; expect ~1e-4–1e-6 relative floating-point noise
125on top of any genuine numerical-method disagreement between solvers.