1# Case: Schnakenberg reaction-diffusion on a triaxial ellipsoid
3A reference solution for validating reaction-diffusion solvers that operate
4on curved surfaces, specifically the discretization of the surface
5Laplace-Beltrami operator. Any solver capable of solving reaction-diffusion
6systems on a closed surface can be checked against this case by starting
7from the same initial condition, running to the same physical end time, and
8comparing the resulting fields in spherical-harmonic coefficient space.
10## Geometry
12A triaxial ellipsoid, parametrized by colatitude `θ ∈ [0, π]` and azimuth
13`φ ∈ [0, 2π)`:
15```
16x = ax · sin(θ) · cos(φ)
17y = ay · sin(θ) · sin(φ)
18z = az · cos(θ)
19```
21with semi-axes `ax`, `ay`, `az` (reference defaults: `1.5, 1.0, 0.6`).
23**The surface is band-limited.** The reference solution was not computed on
24the exact analytic ellipsoid above, but on its spherical-harmonic
25reconstruction truncated at a finite degree `lmax`: the analytic `(x, y, z)`
26was expanded in spherical harmonics and then re-synthesized from that
27truncated expansion, which very slightly rounds off the surface relative to
28the ideal ellipsoid. The truncated surface's own coefficients are included
29in the data file (see below). **A solver being validated against this case
30should reconstruct its working surface (and induced metric) from those
31coefficients, not from the analytic formula**, so that any disagreement in
32the final fields reflects the reaction-diffusion solve and not a difference
33in which surface was actually used.
35## Equations
37The Schnakenberg reaction-diffusion system, with `Δ_g` the Laplace-Beltrami
38operator of the surface above:
40```
41∂u/∂t = D1·Δ_g(u) + a − u + u²v
42∂v/∂t = D2·Δ_g(v) + b − u²v
43```
45| parameter | meaning | reference value |
46|---|---|---|
47| `a` | production rate | 0.1 |
48| `b` | production rate | 0.9 |
49| `D1` | diffusion coefficient of `u` | 4×10⁻⁴ |
50| `D2` | diffusion coefficient of `v` | 8×10⁻³ |
52A solver is free to choose its own numerical method for `Δ_g` and its own
53timestep — those are not part of what is being tested. What matters is
54reaching the same physical end time `T`, computed as `steps × dt` from the
55values stored in the data file's metadata (see below), starting from the
56supplied initial condition.
58## Initial condition
60The initial `u`, `v` fields are supplied as data — spherical-harmonic
61coefficients at `t = 0` — rather than as a formula to regenerate. A solver
62should synthesize its starting fields directly from these coefficients. This
63sidesteps any need to reproduce whatever pseudorandom process originally
64generated the perturbation: the state at `t = 0` is simply data to load.
66## Comparison convention
68To compare against this reference, a solver's own final `u`, `v` fields must
69be expressed in the same spherical-harmonic convention used throughout this
70data set:
72- **Orthonormal** spherical harmonics, **including the Condon-Shortley
73 phase**.
74- Because `u` and `v` are real-valued fields, only coefficients for `m ≥ 0`
75 are stored, as complex numbers. The `m < 0` coefficients are implied by
76 `Q(l, −m) = (−1)^m · conj(Q(l, m))`, and every `m = 0` coefficient has a
77 zero imaginary part.
78- **Ordering**: "m-major" — for `m = 0` up to `lmax`, and within each `m`,
79 for `l = m` up to `lmax`. The flat index of coefficient `(l, m)` is:
81 ```
82 index(l, m) = m·(lmax + 1) − m·(m − 1)/2 + (l − m)
83 ```
85- The total coefficient count is `nlm = (lmax + 1)(lmax + 2)/2`. Coefficients
86 are stored as a flat array of length `2·nlm`, with the real and imaginary
87 part of each coefficient adjacent: `[re(0), im(0), re(1), im(1), ...]`.
89A solver's comparison output should be a flat array of the same length,
90`2·nlm`, in the same ordering, for direct numerical comparison against the
91reference.
93## Data file format
95Each reference run is stored as a single HDF5 file. All coefficient arrays
96are 32-bit floats (the reference solution itself was computed in single
97precision, so no comparison should expect agreement finer than that).
98Scalar metadata is stored as HDF5 attributes, grouped by what they describe:
100```
101/ (attrs: a short description of the run, and the
102 names of the species stored under initial/ and
103 final/, e.g. "U", "V")
104├─ spec/ run parameters as attributes: lmax, seed, steps,
105│ warmup, and the iteration count used internally
106│ by the reference solver's own Δ_g approximation
107│ ├─ params/ the model's own parameters: a, b, D1, D2, dt
108│ └─ geometry_params/ the geometry's own parameters: ax, ay, az
109├─ grid/ the reference solver's own working grid, as
110│ attributes: lmax, mmax, nlat, nphi, nlm
111├─ geometry/
112│ ├─ Gx dataset, float32[2·nlm] — the band-limited
113│ ├─ Gy surface's own spherical-harmonic coefficients,
114│ └─ Gz one per Cartesian component
115├─ initial/ one dataset per species (e.g. U, V), each
116│ float32[2·nlm] — the initial condition at t = 0
117└─ final/ one dataset per species, each float32[2·nlm] —
118 the solution after evolving to t = steps·dt
119```
121`lmax` (under `spec/` and `grid/`) is the spherical-harmonic truncation
122degree that defines `nlm` and the coefficient ordering above; it is the same
123value everywhere it appears in a given file. `seed` identifies which
124pseudorandom perturbation produced the initial condition, but is provenance
125information only — the initial condition itself is fully specified by the
126`initial/` coefficients, and does not need to be regenerated from the seed.
128Every dataset in `geometry/`, `initial/`, and `final/` has length `2·nlm`
129and follows the ordering convention described above.