Case: Schnakenberg reaction-diffusion on a triaxial ellipsoid#

A reference solution for validating reaction-diffusion solvers that operate on curved surfaces, specifically the discretization of the surface Laplace-Beltrami operator. Any solver capable of solving reaction-diffusion systems on a closed surface can be checked against this case by starting from the same initial condition, running to the same physical end time, and comparing the resulting fields in spherical-harmonic coefficient space.

Geometry#

A triaxial ellipsoid, parametrized by colatitude θ ∈ [0, π] and azimuth φ ∈ [0, 2π):

x = ax · sin(θ) · cos(φ)
y = ay · sin(θ) · sin(φ)
z = az · cos(θ)

with semi-axes ax, ay, az (reference defaults: 1.5, 1.0, 0.6).

The surface is band-limited. The reference solution was not computed on the exact analytic ellipsoid above, but on its spherical-harmonic reconstruction truncated at a finite degree lmax: the analytic (x, y, z) was expanded in spherical harmonics and then re-synthesized from that truncated expansion, which very slightly rounds off the surface relative to the ideal ellipsoid. The truncated surface's own coefficients are included in the data file (see below). A solver being validated against this case should reconstruct its working surface (and induced metric) from those coefficients, not from the analytic formula, so that any disagreement in the final fields reflects the reaction-diffusion solve and not a difference in which surface was actually used.

Equations#

The Schnakenberg reaction-diffusion system, with Δ_g the Laplace-Beltrami operator of the surface above:

∂u/∂t = D1·Δ_g(u) + a − u + u²v
∂v/∂t = D2·Δ_g(v) + b     − u²v

A solver is free to choose its own numerical method for Δ_g and its own timestep — those are not part of what is being tested. What matters is reaching the same physical end time T, computed as steps × dt from the values stored in the data file's metadata (see below), starting from the supplied initial condition.

Initial condition#

The initial u, v fields are supplied as data — spherical-harmonic coefficients at t = 0 — rather than as a formula to regenerate. A solver should synthesize its starting fields directly from these coefficients. This sidesteps any need to reproduce whatever pseudorandom process originally generated the perturbation: the state at t = 0 is simply data to load.

Comparison convention#

To compare against this reference, a solver's own final u, v fields must be expressed in the same spherical-harmonic convention used throughout this data set:

• Orthonormal spherical harmonics, including the Condon-Shortley phase.

• Because u and v are real-valued fields, only coefficients for m ≥ 0 are stored, as complex numbers. The m < 0 coefficients are implied by Q(l, −m) = (−1)^m · conj(Q(l, m)), and every m = 0 coefficient has a zero imaginary part.

• Ordering: "m-major" — for m = 0 up to lmax, and within each m, for l = m up to lmax. The flat index of coefficient (l, m) is:

  index(l, m) = m·(lmax + 1) − m·(m − 1)/2 + (l − m)
  

• The total coefficient count is nlm = (lmax + 1)(lmax + 2)/2. Coefficients are stored as a flat array of length 2·nlm, with the real and imaginary part of each coefficient adjacent: [re(0), im(0), re(1), im(1), ...].

A solver's comparison output should be a flat array of the same length, 2·nlm, in the same ordering, for direct numerical comparison against the reference.

Data file format#

Each reference run is stored as a single HDF5 file. All coefficient arrays are 32-bit floats (the reference solution itself was computed in single precision, so no comparison should expect agreement finer than that). Scalar metadata is stored as HDF5 attributes, grouped by what they describe:

/                          (attrs: a short description of the run, and the
                            names of the species stored under initial/ and
                            final/, e.g. "U", "V")
├─ spec/                   run parameters as attributes: lmax, seed, steps,
│                          warmup, and the iteration count used internally
│                          by the reference solver's own Δ_g approximation
│  ├─ params/              the model's own parameters: a, b, D1, D2, dt
│  └─ geometry_params/     the geometry's own parameters: ax, ay, az
├─ grid/                   the reference solver's own working grid, as
│                          attributes: lmax, mmax, nlat, nphi, nlm
├─ geometry/
│  ├─ Gx                   dataset, float32[2·nlm] — the band-limited
│  ├─ Gy                   surface's own spherical-harmonic coefficients,
│  └─ Gz                   one per Cartesian component
├─ initial/                one dataset per species (e.g. U, V), each
│                          float32[2·nlm] — the initial condition at t = 0
└─ final/                  one dataset per species, each float32[2·nlm] —
                           the solution after evolving to t = steps·dt

lmax (under spec/ and grid/) is the spherical-harmonic truncation degree that defines nlm and the coefficient ordering above; it is the same value everywhere it appears in a given file. seed identifies which pseudorandom perturbation produced the initial condition, but is provenance information only — the initial condition itself is fully specified by the initial/ coefficients, and does not need to be regenerated from the seed.

Every dataset in geometry/, initial/, and final/ has length 2·nlm and follows the ordering convention described above.