3Reaction–diffusion systems (Turing patterns) on **closed surfaces given by
4spherical-harmonic embeddings**, solved live in the browser with a spectral
5method whose transforms run on the GPU via WebGPU.
7This is the sibling of
8[turing-sphere](https://github.com/concept-collection/turing-sphere), which
9solves the same systems on the round sphere. Everything there is here; what is
10added is a *surface*.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 12The geometry is in the operator: the models evaluate the surface
13Laplace–Beltrami operator `lap_g` inside the implicit solve, in a **flux form
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 14that costs 5 spherical-harmonic transforms per species per iteration** (plus
15one Legendre-free FFT derivative) where the textbook Cartesian-gradient form
16needs 12. See
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 17[The geometry in the operator](#the-geometry-in-the-operator) and
18[docs/reduced-transforms.md](docs/reduced-transforms.md).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 19
20## What a surface is here
22A geometry is an embedding of the sphere into R³: three scalar fields x, y, z
23over the (θ, φ) parametrization, each carried as spherical-harmonic
24coefficients. The unit sphere is the case where all three are pure degree-1
25harmonics.
27You write one down as MATLAB, in [`geometries/`](geometries/):
29```matlab
30function [gx, gy, gz] = shape(theta, phi, waist, stretch)
31 st = sin(theta);
32 r = 1 - waist * (st .^ 2);
33 gx = r .* (st .* cos(phi));
34 gy = r .* (st .* sin(phi));
35 gz = (1 + stretch) * (r .* cos(theta));
36end
37```
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 39That is ordinary MATLAB. Unlike the models it is not compiled to WGSL: a shape
40is evaluated exactly once at build time, so it runs through numbl's CPU
41interpreter instead, in f64, with the full MATLAB subset available — loops,
42arrays, `min`/`max`, `legendre`, seeded randomness via `rng`/`randn`. The
43result is then **analysed into coefficients**, which is the form everything
44downstream uses. Two things follow from going through the coefficients rather than keeping
47- **It is exactly band-limited at lmax.** The surface has as many derivatives as
48 the scheme needs and no aliased content the solver cannot see. What the solver
49 and the renderer both use is the *synthesis* of the coefficients, so for a
50 shape with sharp features the surface being solved on is not quite the one
51 that was written down — which is the honest thing for a spectral method to do.
52- **It can be evaluated on any grid.** The renderer draws the surface on the
53 (possibly finer) display grid by synthesizing the same coefficients there.
54 That is exact interpolation, not subdivision — the same argument that lets the
55 species fields be oversampled, and it is checked directly in the tests.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 57Five geometries ship: [sphere](geometries/sphere.m) (the reference case),
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 58[ellipsoid](geometries/ellipsoid.m), [peanut](geometries/peanut.m) — a dumbbell
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 59whose waist is a saddle — [bumpy](geometries/bumpy.m), and one random one:
60[blob](geometries/blob.m), surfacefun's blob — the sphere warped by a smooth
61random function built from chebfun's `randnfunsphere` construction (random
62spherical-harmonic coefficients up to degree ⌊2π/λ⌋, rescaled to [−1, 1]).
63It is seeded, so the same seed always gives the same shape; `amp` sets how far
64it departs from the sphere, `λ` how fine its lobes are, and **Re-seed shape**
65draws another one. Each geometry is
66editable in the page, with its own parameters. Changing a shape does not recompile the
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 67solver and does not disturb the run: the geometry is data whose shape in the
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 68bindings depends only on the grid, so a swap is sixteen buffer writes and the
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 69pattern carries straight on.
71A **morph** slider blends the drawn surface back to the unit sphere. The
72parametrization is the sphere's either way, so sweeping it shows which point
73went where.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 75### Seeding, and `tools/`
77A run starts from the uniform steady state plus a small perturbation, and that
78perturbation is a *smooth* random field rather than white noise: chebfun's
79[`randnfun3`](tools/randnfun3.m) on the surface's bounding box, restricted to
80the surface by evaluating it at the grid points — the way surfacefun seeds a
81run. Each model's `init` says so itself:
83```matlab
84function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)
85 f = randnfun3(lam3, gx, gy, gz);
86 ...
87```
89A band-limited seed is fully resolved by the grid, where white noise is
90whatever the grid happened to alias: the tests measure its energy above degree
9120 at 5e-14 of the total, and the flux-form and Algorithm-4 operators now
92track each other to 3e-6 through a run instead of 4e-4. The **seed λ** control
93sets the field's wavelength; smaller means finer features to grow from. It is
94an *absolute* length in the surface's own units, as in chebfun — not a
95fraction of the surface's size — so a larger surface draws more modes at the
96same λ.
98**λ is useful down to about 2π/lmax, and no further.** A field of wavelength λ
99on a unit-radius surface carries angular content up to degree ≈ 2π/λ, so at
100the default lmax 63 the grid holds everything down to λ ≈ 0.1. Past that,
101`init`'s own `analys` discards what the grid cannot represent, and the seed
102gets *weaker* rather than finer while costing eight times as much per halving:
104| λ | 2π/λ | rms of the resolved seed | energy above l=55 | peak degree |
105|---|---|---|---|---|
106| 0.5 | 13 | 2.6e-2 | 1e-8 | 8 |
107| 0.2 | 31 | 2.6e-2 | 1e-8 | 10 |
108| 0.1 | 63 | 2.4e-2 | 0.10 | 44 |
109| 0.05 | 126 | 1.5e-2 | 0.20 | 48 |
110| 0.03 | 209 | 9.6e-3 | 0.27 | 63 |
112Raising lmax moves that floor down, and the seed really does get finer: at
113lmax 127 the same λ=0.05 keeps its full amplitude (2.5e-2 against 1.5e-2 at
114lmax 63) with its peak at degree 79 instead of pinned to the band edge, and
115λ=0.1 becomes *fully* resolved (2e-8 of its energy in the top decile, against
1161e-1 at lmax 63 — so even 0.1 is slightly under-resolved on the default grid).
118Note that lmax cuts both ways: it quadruples npts, so every λ also costs four
119times as much to sum.
121**Nothing caps λ but memory and patience.** The mode table grows to whatever
122is asked for and the only refusal is a table that could not be built at all,
123reported with the mode count it wanted rather than silently truncated. On a
124128×256 grid:
126| λ | modes | seed time |
127|---|---|---|
128| 0.05 | 480,431 | 0.26 s |
129| 0.03 | 2,094,657 | 0.98 s |
130| 0.02 | 6,882,185 | 3.1 s |
131| 0.015 | 16,092,829 | 7.4 s |
132| 0.01 | 53,574,764 | 25.6 s |
134Being slow is the caller's business; **freezing the browser is not**, and at
135these times neither half of the work can be left where it was:
137- The draw is synchronous interpreter time — 13 s at λ=0.01 — which on the
138 main thread stops the page painting and gets it offered up for killing. It
139 runs on a worker instead
140 ([`randnfun3.worker.ts`](src/mgpu/randnfun3.worker.ts)); it touches no GPU
141 and no DOM, so nothing about it needed that thread. Measured during a seed:
142 731 animation frames, no stalled sample.
143- The GPU sum is split across a fixed 16 dispatches (`randnfun3Chunks`)
144 accumulating into the same output, and `submitYielding` ends the submission
145 at each one. A browser's GPU process is shared with compositing, so a single
146 submission running tens of seconds stops *every* tab painting, and one
147 dispatch that long risks the watchdog killing the device outright. Slices
148 past the end of a small table exit immediately, so a coarse λ pays nothing.
150The device is also asked for the adapter's full storage-buffer limit at
151creation ([`src/sht/sht.ts`](src/sht/sht.ts)), so a browser's 128 MB default
152is not what decides how fine λ can be. `seed()` is consequently async.
154`randnfun3` splits across the CPU/GPU line, and the split is forced rather
155than chosen. Drawing the modes needs `randn` and a `sqrt(nnz)` normalization,
156neither of which exists in the compiled WGSL dialect, so the draw is MATLAB in
157[`tools/randnfun3.m`](tools/randnfun3.m) run by the interpreter — a few
158thousand coefficients, ~5 ms. Evaluating is `npts × nmodes` (~6e7 terms at the
159default λ), so that is a WGSL kernel
160([`src/mgpu/randnfun3.ts`](src/mgpu/randnfun3.ts)) reached as an external
161operation, the way `synth` is. The coefficient table is filled in behind the
162call, as `synth` hides its Legendre matrices; λ is not hidden, and the plan
163records which parameter the `.m` asked with so the host draws from that value.
165[`tools/`](tools/) is the shared MATLAB every interpreter run can call, by file
166name, as on MATLAB's path — currently `randnfun3` and
167[`randnfunsphere`](tools/randnfunsphere.m), which `blob.m` is written on. Both
168keep their upstream signatures, including options nothing shipped uses yet
169(`randnfunsphere`'s `'monochromatic'`), because the point of a tool is that a
170geometry you write next can reach for it. Tools are not available to the
171models' *step*, which compiles to WGSL where none of this exists.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 173## The scheme, and where the geometry enters
175It solves the N-species system
177```
178d(u_k)/dt = D_k*lap_g(u_k) + f_k(t, u_1, ..., u_N), k = 1, ..., N
179```
181where `lap_g` is the Laplace–Beltrami operator of the surface. On the round
182sphere `lap_g` is diagonal in spherical-harmonic space with eigenvalues
183`-l(l+1)`, which is what makes turing-sphere's implicit diffusion a single
184divide. On a general surface it is not diagonal, and not even constant-
185coefficient, so that divide has to become a solve.
187The models split the operator:
189```
190lap_g = lap_s + dlap
191```
193with `lap_s` the round-sphere one. `(I - dt*D*lap_s)` is still exactly
194invertible, so the implicit step
196```
197(I - dt*D*lap_g) Unew = B
198```
200rearranges into a fixed point that keeps the whole geometry on the right-hand
201side,
203```
204Unew = (B + dt*D*dlap(Unew)) ./ (1 + dt*D*lam)
205```
207and the loop iterates it from the round-sphere answer. That is preconditioned
208Richardson, with the operator we can invert exactly as the preconditioner; it
209converges while `dt*D*dlap` stays small against `(I - dt*D*lap_s)`, which is
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 210what keeps the cost to a few transforms per step rather than a full elliptic
211solve (see [docs/richardson-iteration.md](docs/richardson-iteration.md)). One
212species of [`models/schnakenberg.m`](models/schnakenberg.m)'s solve loop:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 213
214```matlab
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 215lamJ = lam ./ jhat; % mean-J preconditioner eigenvalues (below)
216...
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 217for k = 1:niter
218 Fu = Un .* filt; % zero the top 2 degrees before differentiating
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 219 vtu = dthetac(Fu);
220 vpu = dphic(Fu);
221 [Ftu, Fpu] = synth(vtu, vpu); % sin(theta)*dtheta(u), dphi(u) -- smooth on
222 % the sphere, one batched dispatch
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 223 Pu = p1 .* Ftu + p2 .* Fpu; % the two fluxes, also smooth: the precomputed
224 Qu = p2 .* Ftu + q2 .* Fpu; % weights carry every 1/sin(theta) there is
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 226 Pcu = PAu .* filt;
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 227 scu = dthetac(Pcu); % theta part of the divergence, coefficients
228 Lu = synth(scu); % sin(theta) * dtheta(P) on the grid
229 dQu = dphig(Qu); % d/dphi is diagonal in the Fourier index:
230 % two FFT stages, no Legendre work at all
231 lapu = r .* (Lu + dQu); % = lap_g(u) on the grid
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 232 dLu = (analys(lapu) + lamJ .* Un) .* filt; % dlap, projected onto the band
233 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
235```
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 237On the sphere `dlap` is mathematically zero — `p1 = q2 = 1`, `p2 = 0`,
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 238`r = 1/sin²θ`, `jhat = 1`, and the composition collapses to `lap_s` — so the
239sphere case reproduces turing-sphere to fp32 round-off, and the tests assert
240the state stays put across 0, 1 and 4 iterations.
242**The preconditioner folds in the symbol of the operator.** `jhat` is the
243host's minimax scale `2/(μmin + μmax)` over the eigenvalues `μ(x)` of the
244operator's principal symbol — the inverse squared principal stretches of
245the embedding, direction included, read straight off the flux-metric
246arrays (`S = (1/J)·[[p1,p2],[p2,q2]]`). Preconditioning with `lam/jhat`
247then contracts every mode *and every direction* at rate
248`(μmax − μmin)/(μmax + μmin) < 1` on any surface, where the plain `lam`
249diverges wherever `μ > 2` — peanut reaches `μ = 6.2`. A det-based mean of
250the area factor (μ's geometric mean, exact only for conformal surfaces) is
251not enough: it under-corrects anisotropic stretching and leaves directional
252high-degree bands with amplification > 1, which surfaced as patterns going
253high-frequency and diverging as `niter` or `lmax` grew. The answer never
254depends on `jhat` — the `lamJ` term added inside `dLu` is the term divided
255back out — only the convergence rate does.
257**The correction is projected onto the band** (`.* filt` on `dLu`,
258matching algos.tex Algorithm 5's zeroing of the top coefficients). Without
259it the top two degrees iterate toward the *undiffused* `Bu` — each solve
260iteration strips a bit more of their implicit diffusion, at species-
261dependent rates, which manufactures a spurious Turing band at the band
262edge: visible on the round sphere as top-degree energy growing ~3%/step at
263`lmax 127, niter 8`. With both fixes the whole niter × geometry sweep
264converges, the spectral centroid of the pattern is resolution-independent
265(l ≈ 26 at lmax 63 and 127 alike), and `jhat: 1` is kept as the divergent
266control in the tests.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 267
268### The geometry in the operator
270`dlap = lap_g - lap_s` is applied to the current iterate at every solve
271iteration, so its transform count is what the whole step's cost scales with.
272Two formulations ship:
2741. **The flux form** (above, all three models): `lap_g u` as the weighted
275 divergence of two weighted fluxes of the sin-scaled derivatives. The
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 276 weights `p2, r, dp1, dq2, jinv` are grid arrays precomputed once per
277 surface from the embedding's θ/φ tangents
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 278 ([`src/geom/metric.ts`](src/geom/metric.ts)), chosen so that **every field
279 that gets analysed is a smooth function on the sphere** — the property
280 that makes spherical-harmonic analysis meaningful, and the entire
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 281 difficulty near the poles. The divergence is split against the round
282 sphere: the sphere's share of it is `-jinv .* lap_s(u)`, exact in
283 spectral space, so `r ~ 1/sin²θ` multiplies only the geometry deviation.
284 Without that split `r` amplifies the polar round-off of the whole flux
285 into a static forcing that nucleates a spot at the pole on every seed.
286 Cost: **6 Legendre transforms** per species per iteration (4 syntheses —
287 two gradient, one divergence, one for the sphere's `-lam .* u`, which
288 rides in the gradient's batch — plus 2 analyses; the phi flux never needs
289 the Legendre basis, `dphig` differentiates it on the grid with two FFT
290 stages, masking m past the top-degree filter, and `dthetac`/`dphic`
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 291 are O(nlm) coefficient shuffles). The derivation, the smoothness
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 292 argument and the fp32 error analysis are in
293 [docs/reduced-transforms.md](docs/reduced-transforms.md).
2942. **The Cartesian-gradient form** (Algorithm 4 of `docs/algos.pdf`), kept as
295 a live reference in
296 [`models/schnakenberg_alg4.m`](models/schnakenberg_alg4.m) and selectable
297 in the app: the surface gradient carried as three ambient components
298 through the inverse metric quantities `Vt*/Vp*`. Cost: **12 transforms**
299 per species per iteration. The tests hold both forms to the same answer on
300 a curved surface, and both metric formulations are precomputed and
301 uploaded for every geometry, so either kind of model runs.
303The θ-derivative machinery both forms need — the α± recurrence
304(`sin θ ∂θ Y_l^m = α⁺Y_{l+1}^m + α⁻Y_{l-1}^m`) as a coefficient-space shuffle
305feeding the existing scalar synthesis — lives in
306[`src/sht/deriv.ts`](src/sht/deriv.ts); no Legendre-derivative tables are
307required.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 308
309### `for` loops, unrolled
311A plan is a fixed list of GPU operations with no branching, which is what makes
312a timestep pure command recording — one submit, no CPU in the loop. A counted
313loop still fits: the planner
314([`src/mgpu/plan.ts`](src/mgpu/plan.ts)) unrolls it, planning the body once per
315iteration.
317Nothing else had to change for that, because numbl gives a variable one cName
318for every assignment to it: the buffer an iteration writes is the buffer the
319next one reads, which is exactly a loop-carried value. The loop variable gets no
320buffer at all — it is bound as a derived scalar to that iteration's literal, so
321a kernel reading `k` folds the number in.
323Two consequences worth stating:
325- **The bounds must be known when the model compiles.** `niter` is supplied as a
326 fixed scalar rather than a tunable one, so changing it recompiles — unlike a
327 parameter, which is a uniform. A runtime bound is refused at compile time with
328 a source position, not silently mis-compiled, and there is a test for that.
329- **Fusion survives.** numbl's inline pass recurses into loop bodies, so a line
330 inside the loop is still one kernel. It runs there with no protected names,
331 though, which means an assignment whose only visible use is later in the same
332 body can be elided — correct for a body-local temp, wrong if something outside
333 the loop wanted it. [`src/mgpu/compile.ts`](src/mgpu/compile.ts) snapshots what
334 each loop body assigns before the pass and refuses the ones that escape, so
335 that case is a compile error rather than a stale read.
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 337Unrolling is exactly linear in the trip count: 18 GPU ops per species per
338iteration (7 transforms, 3 coefficient shuffles, 8 kernels), asserted in the
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 339tests.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 340
341## MATLAB, compiled to WebGPU
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 343Unchanged from turing-sphere. This is the models' path — the geometry files
344instead run once through numbl's CPU interpreter, as above. numbl
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 345parses and lowers each function for the concrete argument types of the current
346grid; its inline pass folds single-use temps back into their consumer, so one
347line of MATLAB becomes one expression tree; and this repo emits one WGSL compute
348kernel per element-wise statement
349([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). `synth` / `analys` are external
350operations whose type rules numbl learns from a `.mtoc2.js` workspace file, and
351which the backend maps onto the spherical-harmonic pipelines. Anything it cannot
352express is refused at compile time with a source position.
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 354The Schnakenberg step compiles to 50 GPU operations at one solve iteration:
35518 transforms, 6 coefficient-space shuffles, 24 generated kernels, and 2
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 356buffer copies feeding the new state back.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 357
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 358**Transforms batch.** The expensive part of every Legendre stage is
359generating the associated Legendre values on the fly by recurrence — work
360that depends only on the grid, not on the field. `synth`/`analys` therefore
361take multiple fields, and a grouped call runs as one batched dispatch: one
362walk of the recurrence, one accumulator lane per field —
364```matlab
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 365[Ftu, Fpu, Ftv, Fpv, Su, Sv] = synth(vtu, vpu, vtv, vpv, lam .* Fu, lam .* Fv);
368The grouping is a promise of independence, never of a lane width: the
369planner ([`src/mgpu/plan.ts`](src/mgpu/plan.ts), `materializeTransforms`)
370chunks each group into whatever the device supports — one ×4 batch under the
371default WebGPU limits, or scalar dispatches with `SHT_BATCH=0` for A/B — so
372the same source runs anywhere. Ungrouped transforms that happen to sit on
373consecutive independent lines are batched the same way. Per-lane arithmetic
374is identical to the scalar kernels', so batched and scalar plans produce
375bit-identical states, asserted in the tests along with compile-time refusal
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 376of a group that drops one of its outputs. All 16 Legendre transforms of the step
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 377above land in batches, worth ~25% of the whole step (0.88 vs 1.14 ms/step at
378lmax 127, 2 iterations, on bumpy).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 379
380Two consequences carried over:
382- **The step is synchronous.** WebGPU's encode path is synchronous and every
383 pipeline is built once at compile time, so a timestep is pure command
384 recording; the only `await` in the loop is the single readback per rendered
385 frame.
386- **Parameters are uniforms, not constants.** Moving a slider rewrites a small
387 buffer instead of triggering a recompile. Editing the MATLAB recompiles;
388 changing `dt` does not. `niter` is the deliberate exception, above.
390## Provenance
392- **turing-sphere**, which this is a fork of: the solver, the transforms
393 backend, the compilation path, the benchmarks and the analytic tests.
394- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
395 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
396 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Vendored under
397 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference transform
398 used for testing.
399- **Rendering:** three.js meshes with per-vertex colormaps, adapted from the
400 `SphereEmbedding` view in
401 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
402 extension package ([`src/render/`](src/render/)). That view displays a
403 time-varying embedded geometry with fields on it, which is the same picture
404 this draws — including its sphere/surface morph, which turing-sphere had
405 dropped as having nothing to morph to.
407turing-sphere additionally carries a comparison against a native build of
408upstream SHTNS ([`bench/shtns/`](https://github.com/concept-collection/turing-sphere/tree/main/bench/shtns)).
409That is not duplicated here: the transforms are the same code, and its C-side
410transcription of the model would have to be maintained against a step this
411project intends to change.
413Because the algorithm is compiled to compute shaders, **WebGPU is required** —
414there is no CPU fallback (the f64 CPU transform remains, for tests).
416## Numerics
418- Grid: Gauss–Legendre × equispaced-φ, dealiased for the cubic reactions with
419 the `(pdeg+1)` rule: `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1`
420 (rounded up to a power of two for the GPU FFT path). At the default lmax 63
421 that is a 128×256 grid.
422- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
423 coefficients for m ≥ 0, m-major ordering.
424- fp32 transforms introduce ~1e-6 relative error per step; for pattern formation
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 425 from a 1e-2 seeded perturbation this is inconsequential. The geometry goes through one
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 426 analysis/synthesis round trip and picks up the same round-off: the unit sphere
427 comes back with radius 1 to ~2e-5 under Dawn, ~4e-4 under SwiftShader.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 428- The shipped analytic geometries are all degree ≤ 5, and the random ones stay
429 near degree 13 at their finest slider settings — far below any lmax the app
430 offers, so band-limiting removes little to nothing from them. A shape you
431 write yourself may not be so lucky — see the note in
432 [`geometries/bumpy.m`](geometries/bumpy.m).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 433
434## Desktop vs browser
436[`scripts/bench.ts`](scripts/bench.ts) runs the same thing the app runs — same
437`.m`, same generated WGSL, same transforms — from Node on desktop WebGPU (Google
438Dawn), and the app prints the command line that reproduces whatever it is
439currently simulating:
441```
442npm run bench -- --preset schnak-spots --geometry ellipsoid --lmax 63 --niter 1 \
443 --steps 2000 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05 \
444 --gax 1.5 --gay 1 --gaz 0.6
445```
447Copy it from under the stats line and compare the `ms/step` it reports with the
448app's. Both sides go through the one shared
449[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
450command, the benchmark parses it back — so there is no second copy of the
451defaults for the two runs to drift apart on. Geometry parameters take a `g`
452prefix (`--gwaist`) so a shape parameter can never collide with a model one.
454The app reports **two** numbers and only the first is comparable to the
455benchmark: `solver` is the batch of steps alone, waited for but not read back;
456`ms/frame` additionally carries a GPU→CPU readback per species, the
457colormapping, and the vertex upload. Those per-frame costs are fixed and do not
458shrink when the GPU gets faster, so on a quick GPU a frame can easily cost ten
459times the steps inside it. That is expected and is not the solver being slower
460in the browser.
462To attribute the gap rather than guess at it:
464```
465node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
466```
468measures the same solver work in both — batched, nothing read back, no rendering
469on either side — and reports each with its CPU-encoding share, the Fourier
470stage, and the adapter. It stops you first if the two are not even the same
471device, which is a common cause of "the browser is much slower". Both sides
472resolve the geometry and the iteration count from the same constants, because
473the iteration count is unrolled into the step and a mismatch would compare two
474different amounts of work.
476The app's **Benchmark** button runs the same measurement in the page, plus the
477**ramp** — the first third of the run against the last. GPUs downclock when
478idle and an animation-paced loop leaves them idle most of every frame, so a
479large ramp means the steady-state number is limited by clocks rather than work.
481### Is it really the same computation?
483```
484node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
485```
487runs one identical spec on the desktop and in a real browser and compares the
488final spectral state. The pipeline is deterministic given (model source,
489geometry, parameters, lmax, niter, seed, steps), so the two should agree to fp32
490round-off — not bit for bit, since GPUs differ in fused-multiply-add and other
491latitude fp32 allows. It also reports which Fourier stage each side chose, since
492FFT and DFT are genuinely different algorithms that round differently.
494Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), an
495optional dependency so that an unsupported platform fails the install of that
496package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
497`--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
498(`;`-separated) passes Dawn options through.
500## Tests
502There is no second implementation of the solver to diff against, so the `.m`
503path is checked against **closed-form answers** and against **exact structural
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 504properties**. Five modules, run in both environments:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 505
506[`test/analyticChecks.ts`](test/analyticChecks.ts) — cases whose evolution is
507known exactly, run through the whole real pipeline. All three are statements
508about the round sphere, so all three build on the sphere geometry:
510- **A** — a linear reaction leaves every mode independent, growing by exactly
511 `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. Pins the transform round trip, the
512 eigenvalue mapping, the IMEX update and the state feedback at once. ~2e-7 over
513 20 steps.
514- **B** — a nonlinear reaction on a uniform field stays uniform, so each step is
515 exactly the scalar ODE map. 1.5e-8 over 25 steps.
516- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
517 linearized 2×2 IMEX recurrence, and `(l=24, m=7)` is confirmed unstable.
518 Looser (~4e-3) because fp32 keeps about four digits of a perturbation that
519 small.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 521[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface, the loop, and
522the seed:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 523
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 524- every geometry evaluates and closes; the sphere has radius 1 everywhere and is
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 525 **exactly degree 1** in the harmonics, which is what makes the reference case
526 exact rather than merely accurate;
527- the peanut matches its own closed-form radial profile at every grid point, and
528 **the same coefficients give the same surface on a 2× grid** — the 2× Gauss
529 latitudes share no point with the 1× ones, so agreeing there is agreeing
530 everywhere, which is what "rendered exactly, not subdivided" means;
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 531- unrolling is **exactly linear** in the trip count, and on the sphere — where
532 the geometric correction is mathematically zero — the state after 20 steps
533 stays within fp32 round-off of the 0-iteration one at 1 and 4 iterations;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 534- a runtime loop bound is refused at compile time;
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 535- swapping the surface mid-run leaves the spectral state untouched;
536- the seed field's **WGSL sum matches the same modes summed in f64 on the
537 CPU** (1.8e-6 over ~1,400 terms) — a kernel misreading the packed mode table
538 would still produce a smooth random-looking field, which no "looks patterned"
539 check would catch; the same seed redraws the same field and a different one
540 does not; the field is band-limited (5e-14 of its energy above degree 20)
541 with λ setting the scale; and a λ finer than the mode table holds is refused
542 rather than silently truncated.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 543
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 544[`test/fluxChecks.ts`](test/fluxChecks.ts) — the six-transform flux-form
545Laplace-Beltrami scheme
546([docs/reduced-transforms.md](docs/reduced-transforms.md)):
548- on the sphere, the precomputed weights match their closed form and the
549 analysed fluxes are **exactly band-limited** (beyond-band tails at f64
550 round-off, ~1e-13), while the deliberately non-smooth control
551 `Q̃/sin θ` keeps a fat tail (~1e-2) — the discrimination the whole scheme
552 rests on;
553- on a non-axisymmetric surface, the flux tails match the Cartesian gradient
554 component's, the doc's §7.1 criterion;
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 555- the compiled op sequences add **5 Legendre transforms per species per
556 iteration against Algorithm 4's 12**, and a real simulation driven by
557 each stays within fp32 accumulation of the other.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 558
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 559[`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
560and asserts **how many kernels it compiles to**, split into the base step and
561what one solve iteration adds. That is a fusion guard: if numbl's inline pass
562stops folding, the results stay correct while every operator becomes its own
563dispatch, which is invisible in the numbers.
565[`test/transformChecks.ts`](test/transformChecks.ts) compares the WGSL transforms
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 566against shtns-webgpu's f64 CPU twin, and holds every compiled batch width to
567the scalar transforms lane by lane; a model run with `SHT_BATCH=0` must
568reproduce the batched run's state exactly.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 569
570- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
571 `--skip-without-gpu` lets a machine without one say so and move on (which is
572 what CI does, since the browser suite covers the same modules).
573- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
574 Also runs the soak. A few geometry tolerances are set by SwiftShader's fp32,
575 which is about an order of magnitude looser than Dawn's.
577Other commands:
579- `npm run bench -- --help` — the desktop benchmark.
580- `npm run bench:sht -- --help` — the transforms alone, no solver.
581- `npx vite-node scripts/diagnose-sht.ts` — when the transform tests fail on a
582 GPU, say *which* stage is wrong.
583- `npx vite-node scripts/diagnose-leg.ts [--m 0]` — read the Legendre recurrence
584 out of the production shader term by term.
585- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
586 pattern saturates rather than decaying or diverging.
587- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
588 sampling JS heap and catching crashes.
589- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
590 demo after a number of steps.
591- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL.
592- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
90108a6Command for testing against a reference implementationOwen Melia 594### Testing against a reference implementation
596Reference solutions live in the sibling
597[turing-surface-test-data](https://github.com/concept-collection/turing-surface-test-data)
598repo, so an independently-written solver never has to depend on this one.
599`cases/schnakenberg-ellipsoid.md` there specifies the one case this repo
600currently ships a reference for;
601[`docs/ellipsoid-reference-spec.md`](docs/ellipsoid-reference-spec.md)
602restates it in this repo's own terms.
604`npm run ref -- --in <file>` (`scripts/ref.ts`) loads a reference file, runs
605the solver from its exact initial spectral state to the same physical end
66ae13eUpdated command now tracks L_infty error too.Owen Melia 606time, and reports the relative-L2 and relative-L-infinity (max-norm) error
607against its final state (plus a geometry sanity check). `--niter` overrides
608the surface-correction iteration count independent of the file, and
609`--tolerance`/`--tolerance-linf` each independently turn their metric into a
3210e7dOpen the reference comparison in one clickJeremy Magland 612The same check runs in the page: **Compare to reference…** picks a `.h5` and
613opens the comparison in one step — the file's own settings (its recorded
614niter, its band, its dt) as the single variant, paused at the file's exact
615initial state, ready to Run. The file defines the whole problem — model,
616parameters, geometry, initial state — and the run stops at the file's end
617time, measured against one extra static row showing its final state on its
618own surface. Watching *where* a variant leaves the reference (rather than
619just reading one number per run) is the point. To widen the study, stop
620comparing, pick more chips, and press Compare — the file stays loaded, with
621the lmax choices floored at its band, since a narrower one could not hold
622its initial state. Reading the file uses
623[h5wasm](https://github.com/usnistgov/h5wasm)'s wasm build, loaded lazily on
624the first file opened.
628```
629npm install
630npm run dev # local dev server
631npm run build # type-check + production build to dist/
632```
634### The numbl dependency
636numbl is a local `file:../../numbl` dependency, so a sibling checkout of
637[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
638compiler internals — parser, lowerer, IR, inline pass — which its package
639`exports` map does not publish, so they are reached through the `numbl-src` path
640alias in [`vite.config.ts`](vite.config.ts).
642The exact surface we depend on is written down in
643[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
644*that*, not against numbl's sources. This keeps this project's compiler settings
645independent of numbl's, and means a change to one of those shapes upstream
646breaks the build here with a clear diff rather than deep inside numbl's tree.
647The `For` IR node is spelled out there, since the planner now walks it.
649CI clones numbl to the sibling path that the `file:` dependency expects, pinned
650to a commit, with `--ignore-scripts` (npm runs a linked package's `prepare`
651script, and numbl's is husky). numbl's own `node_modules` are not needed: the
652slice we import is self-contained TypeScript.
654The `scripts/*.ts` entry points that touch the compiler go through `vite-node`,
655so they resolve imports exactly as the browser build does. Plain `node` cannot:
656numbl's sources import each other as `./foo.js` while the files are `.ts`.
658Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
660## License
662CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).