# Case: Allen–Cahn 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 field 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 field reflects the reaction-diffusion solve and not a difference
in which surface was actually used.

## Equations

The Allen–Cahn equation, with `Δ_g` the Laplace-Beltrami operator of the
surface above. Unlike the Schnakenberg cases in this data set, Allen–Cahn is
a single-species phase-field equation: an interface forms between `u ≈ +1`
and `u ≈ −1` regions and then coarsens over time.

```
∂u/∂t = eps2·Δ_g(u) + u − u³
```

| parameter | meaning | reference value |
|---|---|---|
| `eps2` | interface-width parameter (`ε²`) | 0.01 |

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` field is supplied as data — spherical-harmonic coefficients
at `t = 0` — rather than as a formula to regenerate. A solver should
synthesize its starting field 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` field must be
expressed in the same spherical-harmonic convention used throughout this
data set:

- **Orthonormal** spherical harmonics, **including the Condon-Shortley
  phase**.
- Because `u` is a real-valued field, 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
                            name of the species stored under initial/ and
                            final/, e.g. "U")
├─ 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: eps2, 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 (U), float32[2·nlm] — the initial
│                          condition at t = 0
└─ final/                  one dataset (U), 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.
