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 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> turns the check into a pass/fail (nonzero exit code on failure), for use in CI.

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.