Reference test case: Schnakenberg on a triaxial ellipsoid#

For validating this repo's Laplace-Beltrami surface correction against an independently-implemented reference solver, comparing final state in spherical-harmonic (SH) coefficient space.

Geometry#

A triaxial ellipsoid, gx = ax·sinθ·cosφ, gy = ay·sinθ·sinφ, gz = az·cosθ (geometries/ellipsoid.m; defaults ax=1.5, ay=1.0, az=0.6).

Important: the solver does not run on this analytic surface — it band-limits it first, analysing (gx, gy, gz) into SH coefficients truncated at degree lmax and re-synthesizing before use (src/geom/geometry.ts). To remove geometry-representation error as a confound, a reference file supplies these coefficients directly as /geometry/Gx, /geometry/Gy, /geometry/Gz. The reference solver must reconstruct its surface (and induced metric) by synthesizing these coefficients, not by evaluating the analytic formula above.

Equations#

Schnakenberg reaction-diffusion with the true surface Laplace-Beltrami operator Δ_g (models/schnakenberg.m):

du/dt = D1·Δ_g(u) + a - u + u²v
dv/dt = D2·Δ_g(v) + b     - u²v

This repo's internal discretization (niter Richardson-iteration count for its own Δ_g approximation, dt for its IMEX-Euler timestep) is not part of the equations being tested — the reference solver may use any consistent method for Δ_g and any timestep. It only needs to reach the same physical end time T = steps · dt.

Initial condition#

Loaded from the reference file's /initial/U / /initial/V (t=0, immediately after seeding, before any step), not regenerated — avoids needing to reimplement this repo's PRNG (src/mgpu/noise.ts) to get a matching initial condition.

Output convention (must match exactly, from src/sht/layout.ts)#

• Orthonormal spherical harmonics including Condon-Shortley phase.
• Real field ⇒ complex coefficients stored for m ≥ 0 only: Q_{l,-m} = (-1)^m · conj(Q_lm); m=0 coefficients have zero imaginary part.
• m-major ordering: for m = 0..lmax, for l = m..lmax. index(l,m) = m·(lmax+1) − m·(m−1)/2 + (l−m).
• Flat array, length 2·nlm with nlm = (lmax+1)(lmax+2)/2, [re,im] interleaved per coefficient (qlm[2·index(l,m)], qlm[2·index(l,m)+1]).

The reference solver must project its final u, v onto this same convention/truncation and report flat 2·nlm arrays, to diff directly against the reference file's /final/U / /final/V.

HDF5 file layout#

Each reference file is one .h5 file per run (written with h5wasm; readable from Python with h5py.File(path, "r")). Coefficient datasets are float32, each of length 2·nlm in the convention above. Metadata is stored as attributes, grouped by what it describes rather than as a single flat namespace:

/                             (attrs: command, model, species)
├─ backend/                   (attrs: adapter, runtime, precision)
├─ spec/                      (attrs: preset, geometry, lmax, seed, steps, warmup, niter)
│  ├─ params/                 (attrs: the model's own params, e.g. a, b, D1, D2, dt)
│  └─ geometry_params/        (attrs: the geometry's own params, e.g. ax, ay, az)
├─ grid/                      (attrs: lmax, mmax, nlat, nphi, nlm)
├─ geometry/
│  ├─ Gx                      dataset, float32[2·nlm]
│  ├─ Gy                      dataset, float32[2·nlm]
│  └─ Gz                      dataset, float32[2·nlm]
├─ initial/                   one dataset per species (e.g. U, V), float32[2·nlm] each
└─ final/                     one dataset per species (e.g. U, V), float32[2·nlm] each

species (root attribute) names which datasets live under initial/ and final/ — ["U", "V"] for Schnakenberg. command is the equivalent npm run bench -- invocation, for reproducing the run exactly.

Parameters#

Checking a run against a reference file#

npm run ref -- --in <file> loads a reference file, runs this repo's own solver from its exact initial condition to the same physical end time, and reports the relative-L2 and relative-L-infinity (max-norm) error of the resulting state against the file's final state (and, as a sanity check, of the regenerated geometry against the file's own geometry coefficients — this should be ~0 unless geometry construction itself has changed). --niter <n> overrides the solve's own iteration count for the surface correction, independent of what the reference file was generated with — useful for seeing how much that correction term actually matters for a given run. --tolerance <n> and --tolerance-linf <n> each independently turn their metric into a pass/fail (nonzero exit code on failure), for use in CI.

The browser demo runs the same check visually: Compare → Reference file… loads a reference file into the convergence study, seeds every variant from its exact initial state, runs them side by side to its end time, and shows its final state as one extra static row — on the file's own surface, with each variant's relative-L2 distance to it updating live. Both readers share one parser (src/compare/referenceCase.ts), so the layout above is interpreted identically on the CLI and in the page.

Note the files record only the two endpoint states (initial/, final/) — no intermediate snapshots — so the comparison is meaningful at the end time; the live Δ before that reads as "distance still to the final state".

Caveat#

This repo runs fp32 on GPU; expect ~1e-4–1e-6 relative floating-point noise on top of any genuine numerical-method disagreement between solvers.