Adding specification for the Schnakenberg test
2 changed files+140−1
README.mdmodified+11−1View file
@@ -1 +1,11 @@
1-Test data for turing-surface, to be included as a submodule
1+Test data for turing-surface, to be used as a sibling repository.
2+
3+Reference solutions for validating reaction-diffusion solvers on curved
4+surfaces, so independently-written solvers can be checked against each
5+other. Each case below documents the simulation being solved and the data
6+format its reference file(s) use.
7+
8+## Cases
9+
10+- [`cases/schnakenberg-ellipsoid.md`](cases/schnakenberg-ellipsoid.md) —
11+ Schnakenberg reaction-diffusion on a triaxial ellipsoid.
cases/schnakenberg-ellipsoid.mdadded+129−0View file
@@ -0,0 +1,129 @@
1+# Case: Schnakenberg reaction-diffusion on a triaxial ellipsoid
2+
3+A reference solution for validating reaction-diffusion solvers that operate
4+on curved surfaces, specifically the discretization of the surface
5+Laplace-Beltrami operator. Any solver capable of solving reaction-diffusion
6+systems on a closed surface can be checked against this case by starting
7+from the same initial condition, running to the same physical end time, and
8+comparing the resulting fields in spherical-harmonic coefficient space.
9+
10+## Geometry
11+
12+A triaxial ellipsoid, parametrized by colatitude `θ ∈ [0, π]` and azimuth
13+`φ ∈ [0, 2π)`:
14+
15+```
16+x = ax · sin(θ) · cos(φ)
17+y = ay · sin(θ) · sin(φ)
18+z = az · cos(θ)
19+```
20+
21+with semi-axes `ax`, `ay`, `az` (reference defaults: `1.5, 1.0, 0.6`).
22+
23+**The surface is band-limited.** The reference solution was not computed on
24+the exact analytic ellipsoid above, but on its spherical-harmonic
25+reconstruction truncated at a finite degree `lmax`: the analytic `(x, y, z)`
26+was expanded in spherical harmonics and then re-synthesized from that
27+truncated expansion, which very slightly rounds off the surface relative to
28+the ideal ellipsoid. The truncated surface's own coefficients are included
29+in the data file (see below). **A solver being validated against this case
30+should reconstruct its working surface (and induced metric) from those
31+coefficients, not from the analytic formula**, so that any disagreement in
32+the final fields reflects the reaction-diffusion solve and not a difference
33+in which surface was actually used.
34+
35+## Equations
36+
37+The Schnakenberg reaction-diffusion system, with `Δ_g` the Laplace-Beltrami
38+operator of the surface above:
39+
40+```
41+∂u/∂t = D1·Δ_g(u) + a − u + u²v
42+∂v/∂t = D2·Δ_g(v) + b − u²v
43+```
44+
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⁻³ |
51+
52+A solver is free to choose its own numerical method for `Δ_g` and its own
53+timestep — those are not part of what is being tested. What matters is
54+reaching the same physical end time `T`, computed as `steps × dt` from the
55+values stored in the data file's metadata (see below), starting from the
56+supplied initial condition.
57+
58+## Initial condition
59+
60+The initial `u`, `v` fields are supplied as data — spherical-harmonic
61+coefficients at `t = 0` — rather than as a formula to regenerate. A solver
62+should synthesize its starting fields directly from these coefficients. This
63+sidesteps any need to reproduce whatever pseudorandom process originally
64+generated the perturbation: the state at `t = 0` is simply data to load.
65+
66+## Comparison convention
67+
68+To compare against this reference, a solver's own final `u`, `v` fields must
69+be expressed in the same spherical-harmonic convention used throughout this
70+data set:
71+
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:
80+
81+ ```
82+ index(l, m) = m·(lmax + 1) − m·(m − 1)/2 + (l − m)
83+ ```
84+
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), ...]`.
88+
89+A 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
91+reference.
92+
93+## Data file format
94+
95+Each reference run is stored as a single HDF5 file. All coefficient arrays
96+are 32-bit floats (the reference solution itself was computed in single
97+precision, so no comparison should expect agreement finer than that).
98+Scalar metadata is stored as HDF5 attributes, grouped by what they describe:
99+
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+```
120+
121+`lmax` (under `spec/` and `grid/`) is the spherical-harmonic truncation
122+degree that defines `nlm` and the coefficient ordering above; it is the same
123+value everywhere it appears in a given file. `seed` identifies which
124+pseudorandom perturbation produced the initial condition, but is provenance
125+information only — the initial condition itself is fully specified by the
126+`initial/` coefficients, and does not need to be regenerated from the seed.
127+
128+Every dataset in `geometry/`, `initial/`, and `final/` has length `2·nlm`
129+and follows the ordering convention described above.