# 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](https://github.com/usnistgov/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

| name | meaning | default |
|---|---|---|
| `a`, `b` | Schnakenberg kinetics | 0.1, 0.9 |
| `D1`, `D2` | diffusion coefficients | 4e-4, 8e-3 |
| `ax`, `ay`, `az` | ellipsoid semi-axes | 1.5, 1.0, 0.6 |
| `lmax` | SH truncation degree | 63 |
| `T = steps·dt` | physical end time | e.g. 2000·0.05 = 100 |
| `seed` | provenance only — IC supplied as coefficients | 1 |

## 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.

## 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.