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
14that costs 6 spherical-harmonic transforms per species per iteration** where
15the textbook Cartesian-gradient form needs 12. See
16[The geometry in the operator](#the-geometry-in-the-operator) and
17[docs/reduced-transforms.md](docs/reduced-transforms.md).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 18
19## What a surface is here
21A geometry is an embedding of the sphere into R³: three scalar fields x, y, z
22over the (θ, φ) parametrization, each carried as spherical-harmonic
23coefficients. The unit sphere is the case where all three are pure degree-1
24harmonics.
26You write one down as MATLAB, in [`geometries/`](geometries/):
28```matlab
29function [gx, gy, gz] = shape(theta, phi, waist, stretch)
30 st = sin(theta);
31 r = 1 - waist * (st .^ 2);
32 gx = r .* (st .* cos(phi));
33 gy = r .* (st .* sin(phi));
34 gz = (1 + stretch) * (r .* cos(theta));
35end
36```
38That is ordinary element-wise MATLAB and goes through the same compiler and the
39same WGSL backend the models do. It is evaluated once on the solver's grid, and
40then **analysed into coefficients**, which is the form everything downstream
41uses. Two things follow from going through the coefficients rather than keeping
42the pointwise values:
44- **It is exactly band-limited at lmax.** The surface has as many derivatives as
45 the scheme needs and no aliased content the solver cannot see. What the solver
46 and the renderer both use is the *synthesis* of the coefficients, so for a
47 shape with sharp features the surface being solved on is not quite the one
48 that was written down — which is the honest thing for a spectral method to do.
49- **It can be evaluated on any grid.** The renderer draws the surface on the
50 (possibly finer) display grid by synthesizing the same coefficients there.
51 That is exact interpolation, not subdivision — the same argument that lets the
52 species fields be oversampled, and it is checked directly in the tests.
54Four geometries ship: [sphere](geometries/sphere.m) (the reference case),
55[ellipsoid](geometries/ellipsoid.m), [peanut](geometries/peanut.m) — a dumbbell
56whose waist is a saddle — and [bumpy](geometries/bumpy.m). Each is editable in
57the page, with its own parameters. Changing a shape does not recompile the
58solver 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 59bindings depends only on the grid, so a swap is sixteen buffer writes and the
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 60pattern carries straight on.
62A **morph** slider blends the drawn surface back to the unit sphere. The
63parametrization is the sphere's either way, so sweeping it shows which point
64went where.
66## The scheme, and where the geometry enters
68It solves the N-species system
70```
71d(u_k)/dt = D_k*lap_g(u_k) + f_k(t, u_1, ..., u_N), k = 1, ..., N
72```
74where `lap_g` is the Laplace–Beltrami operator of the surface. On the round
75sphere `lap_g` is diagonal in spherical-harmonic space with eigenvalues
76`-l(l+1)`, which is what makes turing-sphere's implicit diffusion a single
77divide. On a general surface it is not diagonal, and not even constant-
78coefficient, so that divide has to become a solve.
80The models split the operator:
82```
83lap_g = lap_s + dlap
84```
86with `lap_s` the round-sphere one. `(I - dt*D*lap_s)` is still exactly
87invertible, so the implicit step
89```
90(I - dt*D*lap_g) Unew = B
91```
93rearranges into a fixed point that keeps the whole geometry on the right-hand
94side,
96```
97Unew = (B + dt*D*dlap(Unew)) ./ (1 + dt*D*lam)
98```
100and the loop iterates it from the round-sphere answer. That is preconditioned
101Richardson, with the operator we can invert exactly as the preconditioner; it
102converges 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 103what keeps the cost to a few transforms per step rather than a full elliptic
104solve (see [docs/richardson-iteration.md](docs/richardson-iteration.md)). One
105species of [`models/schnakenberg.m`](models/schnakenberg.m)'s solve loop:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 106
107```matlab
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 108for k = 1:niter
109 Fu = Un .* filt; % zero the top 2 degrees before differentiating
110 Ftu = synth(dthetac(Fu)); % sin(theta) * dtheta(u) -- smooth on the sphere
111 Fpu = synth(dphic(Fu)); % dphi(u) -- smooth on the sphere
112 Pu = p1 .* Ftu + p2 .* Fpu; % the two fluxes, also smooth: the precomputed
113 Qu = p2 .* Ftu + q2 .* Fpu; % weights carry every 1/sin(theta) there is
114 Pcu = analys(Pu) .* filt;
115 Qcu = analys(Qu) .* filt;
116 scu = dthetac(Pcu) + dphic(Qcu); % divergence, in coefficient space
117 lapu = r .* synth(scu); % = lap_g(u) on the grid
118 dLu = analys(lapu) + lam .* Un; % dlap = lap_g - lap_s
119 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
121```
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 123On the sphere `dlap` is mathematically zero — `p1 = q2 = 1`, `p2 = 0`,
124`r = 1/sin²θ`, and the composition collapses to `lap_s` — so the sphere case
125reproduces turing-sphere to fp32 round-off, and the tests assert the state
126stays put across 0, 1 and 4 iterations.
128### The geometry in the operator
130`dlap = lap_g - lap_s` is applied to the current iterate at every solve
131iteration, so its transform count is what the whole step's cost scales with.
132Two formulations ship:
1341. **The flux form** (above, all three models): `lap_g u` as the weighted
135 divergence of two weighted fluxes of the sin-scaled derivatives. The
136 weights `p1, p2, q2, r` are grid arrays precomputed once per surface from
137 the embedding's θ/φ tangents
138 ([`src/geom/metric.ts`](src/geom/metric.ts)), chosen so that **every field
139 that gets analysed is a smooth function on the sphere** — the property
140 that makes spherical-harmonic analysis meaningful, and the entire
141 difficulty near the poles. Cost: **6 transforms** per species per
142 iteration (3 syntheses + 3 analyses; `dthetac`/`dphic` are O(nlm)
143 coefficient shuffles, not transforms). The derivation, the smoothness
144 argument and the fp32 error analysis are in
145 [docs/reduced-transforms.md](docs/reduced-transforms.md).
1462. **The Cartesian-gradient form** (Algorithm 4 of `docs/algos.pdf`), kept as
147 a live reference in
148 [`models/schnakenberg_alg4.m`](models/schnakenberg_alg4.m) and selectable
149 in the app: the surface gradient carried as three ambient components
150 through the inverse metric quantities `Vt*/Vp*`. Cost: **12 transforms**
151 per species per iteration. The tests hold both forms to the same answer on
152 a curved surface, and both metric formulations are precomputed and
153 uploaded for every geometry, so either kind of model runs.
155The θ-derivative machinery both forms need — the α± recurrence
156(`sin θ ∂θ Y_l^m = α⁺Y_{l+1}^m + α⁻Y_{l-1}^m`) as a coefficient-space shuffle
157feeding the existing scalar synthesis — lives in
158[`src/sht/deriv.ts`](src/sht/deriv.ts); no Legendre-derivative tables are
159required.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 160
161### `for` loops, unrolled
163A plan is a fixed list of GPU operations with no branching, which is what makes
164a timestep pure command recording — one submit, no CPU in the loop. A counted
165loop still fits: the planner
166([`src/mgpu/plan.ts`](src/mgpu/plan.ts)) unrolls it, planning the body once per
167iteration.
169Nothing else had to change for that, because numbl gives a variable one cName
170for every assignment to it: the buffer an iteration writes is the buffer the
171next one reads, which is exactly a loop-carried value. The loop variable gets no
172buffer at all — it is bound as a derived scalar to that iteration's literal, so
173a kernel reading `k` folds the number in.
175Two consequences worth stating:
177- **The bounds must be known when the model compiles.** `niter` is supplied as a
178 fixed scalar rather than a tunable one, so changing it recompiles — unlike a
179 parameter, which is a uniform. A runtime bound is refused at compile time with
180 a source position, not silently mis-compiled, and there is a test for that.
181- **Fusion survives.** numbl's inline pass recurses into loop bodies, so a line
182 inside the loop is still one kernel. It runs there with no protected names,
183 though, which means an assignment whose only visible use is later in the same
184 body can be elided — correct for a body-local temp, wrong if something outside
185 the loop wanted it. [`src/mgpu/compile.ts`](src/mgpu/compile.ts) snapshots what
186 each loop body assigns before the pass and refuses the ones that escape, so
187 that case is a compile error rather than a stale read.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 189Unrolling is exactly linear in the trip count: 19 GPU ops per species per
190iteration (6 transforms, 4 coefficient shuffles, 9 kernels), asserted in the
191tests.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 192
193## MATLAB, compiled to WebGPU
195Unchanged from turing-sphere, and it now compiles the geometry files too. numbl
196parses and lowers each function for the concrete argument types of the current
197grid; its inline pass folds single-use temps back into their consumer, so one
198line of MATLAB becomes one expression tree; and this repo emits one WGSL compute
199kernel per element-wise statement
200([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). `synth` / `analys` are external
201operations whose type rules numbl learns from a `.mtoc2.js` workspace file, and
202which the backend maps onto the spherical-harmonic pipelines. Anything it cannot
203express is refused at compile time with a source position.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 205The Schnakenberg step compiles to 51 GPU operations at one solve iteration:
20616 transforms, 8 coefficient-space shuffles, 25 generated kernels, and 2
207buffer copies feeding the new state back.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 208
209Two consequences carried over:
211- **The step is synchronous.** WebGPU's encode path is synchronous and every
212 pipeline is built once at compile time, so a timestep is pure command
213 recording; the only `await` in the loop is the single readback per rendered
214 frame.
215- **Parameters are uniforms, not constants.** Moving a slider rewrites a small
216 buffer instead of triggering a recompile. Editing the MATLAB recompiles;
217 changing `dt` does not. `niter` is the deliberate exception, above.
219## Provenance
221- **turing-sphere**, which this is a fork of: the solver, the transforms
222 backend, the compilation path, the benchmarks and the analytic tests.
223- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
224 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
225 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Vendored under
226 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference transform
227 used for testing.
228- **Rendering:** three.js meshes with per-vertex colormaps, adapted from the
229 `SphereEmbedding` view in
230 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
231 extension package ([`src/render/`](src/render/)). That view displays a
232 time-varying embedded geometry with fields on it, which is the same picture
233 this draws — including its sphere/surface morph, which turing-sphere had
234 dropped as having nothing to morph to.
236turing-sphere additionally carries a comparison against a native build of
237upstream SHTNS ([`bench/shtns/`](https://github.com/concept-collection/turing-sphere/tree/main/bench/shtns)).
238That is not duplicated here: the transforms are the same code, and its C-side
239transcription of the model would have to be maintained against a step this
240project intends to change.
242Because the algorithm is compiled to compute shaders, **WebGPU is required** —
243there is no CPU fallback (the f64 CPU transform remains, for tests).
245## Numerics
247- Grid: Gauss–Legendre × equispaced-φ, dealiased for the cubic reactions with
248 the `(pdeg+1)` rule: `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1`
249 (rounded up to a power of two for the GPU FFT path). At the default lmax 63
250 that is a 128×256 grid.
251- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
252 coefficients for m ≥ 0, m-major ordering.
253- fp32 transforms introduce ~1e-6 relative error per step; for pattern formation
254 from 1e-2 seeded noise this is inconsequential. The geometry goes through one
255 analysis/synthesis round trip and picks up the same round-off: the unit sphere
256 comes back with radius 1 to ~2e-5 under Dawn, ~4e-4 under SwiftShader.
257- The shipped geometries are all degree ≤ 5, far below any lmax the app offers,
258 so band-limiting removes nothing from them. A shape you write yourself may not
259 be so lucky — see the note in [`geometries/bumpy.m`](geometries/bumpy.m).
261## Desktop vs browser
263[`scripts/bench.ts`](scripts/bench.ts) runs the same thing the app runs — same
264`.m`, same generated WGSL, same transforms — from Node on desktop WebGPU (Google
265Dawn), and the app prints the command line that reproduces whatever it is
266currently simulating:
268```
269npm run bench -- --preset schnak-spots --geometry ellipsoid --lmax 63 --niter 1 \
270 --steps 2000 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05 \
271 --gax 1.5 --gay 1 --gaz 0.6
272```
274Copy it from under the stats line and compare the `ms/step` it reports with the
275app's. Both sides go through the one shared
276[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
277command, the benchmark parses it back — so there is no second copy of the
278defaults for the two runs to drift apart on. Geometry parameters take a `g`
279prefix (`--gwaist`) so a shape parameter can never collide with a model one.
281The app reports **two** numbers and only the first is comparable to the
282benchmark: `solver` is the batch of steps alone, waited for but not read back;
283`ms/frame` additionally carries a GPU→CPU readback per species, the
284colormapping, and the vertex upload. Those per-frame costs are fixed and do not
285shrink when the GPU gets faster, so on a quick GPU a frame can easily cost ten
286times the steps inside it. That is expected and is not the solver being slower
287in the browser.
289To attribute the gap rather than guess at it:
291```
292node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
293```
295measures the same solver work in both — batched, nothing read back, no rendering
296on either side — and reports each with its CPU-encoding share, the Fourier
297stage, and the adapter. It stops you first if the two are not even the same
298device, which is a common cause of "the browser is much slower". Both sides
299resolve the geometry and the iteration count from the same constants, because
300the iteration count is unrolled into the step and a mismatch would compare two
301different amounts of work.
303The app's **Benchmark** button runs the same measurement in the page, plus the
304**ramp** — the first third of the run against the last. GPUs downclock when
305idle and an animation-paced loop leaves them idle most of every frame, so a
306large ramp means the steady-state number is limited by clocks rather than work.
308### Is it really the same computation?
310```
311node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
312```
314runs one identical spec on the desktop and in a real browser and compares the
315final spectral state. The pipeline is deterministic given (model source,
316geometry, parameters, lmax, niter, seed, steps), so the two should agree to fp32
317round-off — not bit for bit, since GPUs differ in fused-multiply-add and other
318latitude fp32 allows. It also reports which Fourier stage each side chose, since
319FFT and DFT are genuinely different algorithms that round differently.
321Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), an
322optional dependency so that an unsupported platform fails the install of that
323package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
324`--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
325(`;`-separated) passes Dawn options through.
327## Tests
329There is no second implementation of the solver to diff against, so the `.m`
330path is checked against **closed-form answers** and against **exact structural
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 331properties**. Five modules, run in both environments:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 332
333[`test/analyticChecks.ts`](test/analyticChecks.ts) — cases whose evolution is
334known exactly, run through the whole real pipeline. All three are statements
335about the round sphere, so all three build on the sphere geometry:
337- **A** — a linear reaction leaves every mode independent, growing by exactly
338 `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. Pins the transform round trip, the
339 eigenvalue mapping, the IMEX update and the state feedback at once. ~2e-7 over
340 20 steps.
341- **B** — a nonlinear reaction on a uniform field stays uniform, so each step is
342 exactly the scalar ODE map. 1.5e-8 over 25 steps.
343- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
344 linearized 2×2 IMEX recurrence, and `(l=24, m=7)` is confirmed unstable.
345 Looser (~4e-3) because fp32 keeps about four digits of a perturbation that
346 small.
348[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface and the loop:
350- every geometry compiles and closes; the sphere has radius 1 everywhere and is
351 **exactly degree 1** in the harmonics, which is what makes the reference case
352 exact rather than merely accurate;
353- the peanut matches its own closed-form radial profile at every grid point, and
354 **the same coefficients give the same surface on a 2× grid** — the 2× Gauss
355 latitudes share no point with the 1× ones, so agreeing there is agreeing
356 everywhere, which is what "rendered exactly, not subdivided" means;
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 357- unrolling is **exactly linear** in the trip count, and on the sphere — where
358 the geometric correction is mathematically zero — the state after 20 steps
359 stays within fp32 round-off of the 0-iteration one at 1 and 4 iterations;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 360- a runtime loop bound is refused at compile time;
361- swapping the surface mid-run leaves the spectral state untouched.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 363[`test/fluxChecks.ts`](test/fluxChecks.ts) — the six-transform flux-form
364Laplace-Beltrami scheme
365([docs/reduced-transforms.md](docs/reduced-transforms.md)):
367- on the sphere, the precomputed weights match their closed form and the
368 analysed fluxes are **exactly band-limited** (beyond-band tails at f64
369 round-off, ~1e-13), while the deliberately non-smooth control
370 `Q̃/sin θ` keeps a fat tail (~1e-2) — the discrimination the whole scheme
371 rests on;
372- on a non-axisymmetric surface, the flux tails match the Cartesian gradient
373 component's, the doc's §7.1 criterion;
374- the compiled op sequences add **6 transforms per species per iteration
375 against Algorithm 4's 12**, and a real simulation driven by each stays
376 within fp32 accumulation of the other.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 378[`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
379and asserts **how many kernels it compiles to**, split into the base step and
380what one solve iteration adds. That is a fusion guard: if numbl's inline pass
381stops folding, the results stay correct while every operator becomes its own
382dispatch, which is invisible in the numbers.
384[`test/transformChecks.ts`](test/transformChecks.ts) compares the WGSL transforms
385against shtns-webgpu's f64 CPU twin.
387- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
388 `--skip-without-gpu` lets a machine without one say so and move on (which is
389 what CI does, since the browser suite covers the same modules).
390- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
391 Also runs the soak. A few geometry tolerances are set by SwiftShader's fp32,
392 which is about an order of magnitude looser than Dawn's.
394Other commands:
396- `npm run bench -- --help` — the desktop benchmark.
397- `npm run bench:sht -- --help` — the transforms alone, no solver.
398- `npx vite-node scripts/diagnose-sht.ts` — when the transform tests fail on a
399 GPU, say *which* stage is wrong.
400- `npx vite-node scripts/diagnose-leg.ts [--m 0]` — read the Legendre recurrence
401 out of the production shader term by term.
402- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
403 pattern saturates rather than decaying or diverging.
404- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
405 sampling JS heap and catching crashes.
406- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
407 demo after a number of steps.
408- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL.
409- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
411## Development
413```
414npm install
415npm run dev # local dev server
416npm run build # type-check + production build to dist/
417```
419### The numbl dependency
421numbl is a local `file:../../numbl` dependency, so a sibling checkout of
422[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
423compiler internals — parser, lowerer, IR, inline pass — which its package
424`exports` map does not publish, so they are reached through the `numbl-src` path
425alias in [`vite.config.ts`](vite.config.ts).
427The exact surface we depend on is written down in
428[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
429*that*, not against numbl's sources. This keeps this project's compiler settings
430independent of numbl's, and means a change to one of those shapes upstream
431breaks the build here with a clear diff rather than deep inside numbl's tree.
432The `For` IR node is spelled out there, since the planner now walks it.
434CI clones numbl to the sibling path that the `file:` dependency expects, pinned
435to a commit, with `--ignore-scripts` (npm runs a linked package's `prepare`
436script, and numbl's is husky). numbl's own `node_modules` are not needed: the
437slice we import is self-contained TypeScript.
439The `scripts/*.ts` entry points that touch the compiler go through `vite-node`,
440so they resolve imports exactly as the browser build does. Plain `node` cannot:
441numbl's sources import each other as `./foo.js` while the files are `.ts`.
443Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
445## License
447CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).