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*.
12> [!WARNING]
13> **The geometry is rendered, not yet solved on.** The Laplace–Beltrami
14> operator in the models is still the round sphere's — the term that carries
15> the shape is a placeholder that is identically zero. On anything but the
16> sphere you are looking at the sphere's pattern painted onto that surface, not
17> the pattern that surface would grow. Everything the correction needs in order
18> to be dropped in — the embedding, the split of the operator, the iterative
19> solve, the unrolled loop — is built and tested. See
20> [The geometry is not in the operator yet](#the-geometry-is-not-in-the-operator-yet).
22## What a surface is here
24A geometry is an embedding of the sphere into R³: three scalar fields x, y, z
25over the (θ, φ) parametrization, each carried as spherical-harmonic
26coefficients. The unit sphere is the case where all three are pure degree-1
27harmonics.
29You write one down as MATLAB, in [`geometries/`](geometries/):
31```matlab
32function [gx, gy, gz] = shape(theta, phi, waist, stretch)
33 st = sin(theta);
34 r = 1 - waist * (st .^ 2);
35 gx = r .* (st .* cos(phi));
36 gy = r .* (st .* sin(phi));
37 gz = (1 + stretch) * (r .* cos(theta));
38end
39```
41That is ordinary element-wise MATLAB and goes through the same compiler and the
42same WGSL backend the models do. It is evaluated once on the solver's grid, and
43then **analysed into coefficients**, which is the form everything downstream
44uses. Two things follow from going through the coefficients rather than keeping
45the pointwise values:
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.
57Four geometries ship: [sphere](geometries/sphere.m) (the reference case),
58[ellipsoid](geometries/ellipsoid.m), [peanut](geometries/peanut.m) — a dumbbell
59whose waist is a saddle — and [bumpy](geometries/bumpy.m). Each is editable in
60the page, with its own parameters. Changing a shape does not recompile the
61solver and does not disturb the run: the geometry is data whose shape in the
62bindings depends only on the grid, so a swap is six buffer writes and the
63pattern carries straight on.
65A **morph** slider blends the drawn surface back to the unit sphere. The
66parametrization is the sphere's either way, so sweeping it shows which point
67went where.
69## The scheme, and where the geometry enters
71It solves the N-species system
73```
74d(u_k)/dt = D_k*lap_g(u_k) + f_k(t, u_1, ..., u_N), k = 1, ..., N
75```
77where `lap_g` is the Laplace–Beltrami operator of the surface. On the round
78sphere `lap_g` is diagonal in spherical-harmonic space with eigenvalues
79`-l(l+1)`, which is what makes turing-sphere's implicit diffusion a single
80divide. On a general surface it is not diagonal, and not even constant-
81coefficient, so that divide has to become a solve.
83The models split the operator:
85```
86lap_g = lap_s + dlap
87```
89with `lap_s` the round-sphere one. `(I - dt*D*lap_s)` is still exactly
90invertible, so the implicit step
92```
93(I - dt*D*lap_g) Unew = B
94```
96rearranges into a fixed point that keeps the whole geometry on the right-hand
97side,
99```
100Unew = (B + dt*D*dlap(Unew)) ./ (1 + dt*D*lam)
101```
103and the loop iterates it from the round-sphere answer. That is preconditioned
104Richardson, with the operator we can invert exactly as the preconditioner; it
105converges while `dt*D*dlap` stays small against `(I - dt*D*lap_s)`, which is
106what would keep the cost to a few transforms per step rather than a full
107elliptic solve. Written out, the whole of
108[`models/schnakenberg.m`](models/schnakenberg.m)'s step is:
110```matlab
111function [Un, Vn, u, v] = step(U, V, lam, gx, gy, gz, a, b, D1, D2, dt, niter)
112 u = synth(U);
113 v = synth(V);
114 uuv = u .* u .* v;
116 Bu = U + dt * analys(a - u + uuv);
117 Bv = V + dt * analys(b - uuv);
119 Un = Bu ./ (1 + (dt * D1) * lam);
120 Vn = Bv ./ (1 + (dt * D2) * lam);
122 for k = 1:niter
123 dLu = 0 * Un; % <- the placeholder
124 dLv = 0 * Vn;
125 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
126 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
127 end
128end
129```
131Written this way rather than as a residual correction on purpose: with `dlap`
132zero, every iterate is *bit for bit* the first line, with no cancellation to
133round differently. So the sphere case is not "close to" turing-sphere, it is
134the same arithmetic, and the tests assert exactly that — the state after 20
135steps is identical at 0, 1 and 4 iterations.
137### The geometry is not in the operator yet
139What belongs where `dLu` is now is `dlap = lap_g - lap_s` applied to the current
140iterate. Getting it needs two things this repo does not have:
1421. **The induced metric**, `g_ij = ∂_i X · ∂_j X` for `X = (gx, gy, gz)`. The
143 geometry is static and low-degree, so this is a one-off precomputation, not
144 per-step work — but it needs θ- and φ-derivatives of the embedding.
1452. **Surface derivatives of the field**, per iteration. In the round frame this
146 is the spheroidal transform pair — SHTNS's `SHsph_to_spat` and
147 `spat_to_SHsph`, i.e. `grad_s` and `div_s` — which lets the operator be
148 written as `div_s(A grad_s f)` with `A` built from the metric, with no
149 explicit `1/sin θ` to go singular at the poles.
151Both need Legendre *derivative* tables, which the vendored WGSL transforms under
152[`src/sht/`](src/sht/) do not implement — they are scalar synthesis and analysis
153only. That is the missing piece, and it is a substantial addition to the
154transforms rather than a change to the models. Until it lands, the models take
155`gx, gy, gz` (the surface on the grid) and `Gx, Gy, Gz` (the same surface as
156coefficients) as arguments and do not use them, and the app says so.
158### `for` loops, unrolled
160A plan is a fixed list of GPU operations with no branching, which is what makes
161a timestep pure command recording — one submit, no CPU in the loop. A counted
162loop still fits: the planner
163([`src/mgpu/plan.ts`](src/mgpu/plan.ts)) unrolls it, planning the body once per
164iteration.
166Nothing else had to change for that, because numbl gives a variable one cName
167for every assignment to it: the buffer an iteration writes is the buffer the
168next one reads, which is exactly a loop-carried value. The loop variable gets no
169buffer at all — it is bound as a derived scalar to that iteration's literal, so
170a kernel reading `k` folds the number in.
172Two consequences worth stating:
174- **The bounds must be known when the model compiles.** `niter` is supplied as a
175 fixed scalar rather than a tunable one, so changing it recompiles — unlike a
176 parameter, which is a uniform. A runtime bound is refused at compile time with
177 a source position, not silently mis-compiled, and there is a test for that.
178- **Fusion survives.** numbl's inline pass recurses into loop bodies, so a line
179 inside the loop is still one kernel. It runs there with no protected names,
180 though, which means an assignment whose only visible use is later in the same
181 body can be elided — correct for a body-local temp, wrong if something outside
182 the loop wanted it. [`src/mgpu/compile.ts`](src/mgpu/compile.ts) snapshots what
183 each loop body assigns before the pass and refuses the ones that escape, so
184 that case is a compile error rather than a stale read.
186Unrolling is exactly linear in the trip count: 2 GPU ops per species per
187iteration, asserted in the tests.
189## MATLAB, compiled to WebGPU
191Unchanged from turing-sphere, and it now compiles the geometry files too. numbl
192parses and lowers each function for the concrete argument types of the current
193grid; its inline pass folds single-use temps back into their consumer, so one
194line of MATLAB becomes one expression tree; and this repo emits one WGSL compute
195kernel per element-wise statement
196([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). `synth` / `analys` are external
197operations whose type rules numbl learns from a `.mtoc2.js` workspace file, and
198which the backend maps onto the spherical-harmonic pipelines. Anything it cannot
199express is refused at compile time with a source position.
201The Schnakenberg step above compiles to 17 GPU operations at one solve
202iteration: 4 transforms, 11 generated kernels, and 2 buffer copies feeding the
203new state back.
205Two consequences carried over:
207- **The step is synchronous.** WebGPU's encode path is synchronous and every
208 pipeline is built once at compile time, so a timestep is pure command
209 recording; the only `await` in the loop is the single readback per rendered
210 frame.
211- **Parameters are uniforms, not constants.** Moving a slider rewrites a small
212 buffer instead of triggering a recompile. Editing the MATLAB recompiles;
213 changing `dt` does not. `niter` is the deliberate exception, above.
215## Provenance
217- **turing-sphere**, which this is a fork of: the solver, the transforms
218 backend, the compilation path, the benchmarks and the analytic tests.
219- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
220 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
221 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Vendored under
222 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference transform
223 used for testing.
224- **Rendering:** three.js meshes with per-vertex colormaps, adapted from the
225 `SphereEmbedding` view in
226 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
227 extension package ([`src/render/`](src/render/)). That view displays a
228 time-varying embedded geometry with fields on it, which is the same picture
229 this draws — including its sphere/surface morph, which turing-sphere had
230 dropped as having nothing to morph to.
232turing-sphere additionally carries a comparison against a native build of
233upstream SHTNS ([`bench/shtns/`](https://github.com/concept-collection/turing-sphere/tree/main/bench/shtns)).
234That is not duplicated here: the transforms are the same code, and its C-side
235transcription of the model would have to be maintained against a step this
236project intends to change.
238Because the algorithm is compiled to compute shaders, **WebGPU is required** —
239there is no CPU fallback (the f64 CPU transform remains, for tests).
241## Numerics
243- Grid: Gauss–Legendre × equispaced-φ, dealiased for the cubic reactions with
244 the `(pdeg+1)` rule: `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1`
245 (rounded up to a power of two for the GPU FFT path). At the default lmax 63
246 that is a 128×256 grid.
247- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
248 coefficients for m ≥ 0, m-major ordering.
249- fp32 transforms introduce ~1e-6 relative error per step; for pattern formation
250 from 1e-2 seeded noise this is inconsequential. The geometry goes through one
251 analysis/synthesis round trip and picks up the same round-off: the unit sphere
252 comes back with radius 1 to ~2e-5 under Dawn, ~4e-4 under SwiftShader.
253- The shipped geometries are all degree ≤ 5, far below any lmax the app offers,
254 so band-limiting removes nothing from them. A shape you write yourself may not
255 be so lucky — see the note in [`geometries/bumpy.m`](geometries/bumpy.m).
257## Desktop vs browser
259[`scripts/bench.ts`](scripts/bench.ts) runs the same thing the app runs — same
260`.m`, same generated WGSL, same transforms — from Node on desktop WebGPU (Google
261Dawn), and the app prints the command line that reproduces whatever it is
262currently simulating:
264```
265npm run bench -- --preset schnak-spots --geometry ellipsoid --lmax 63 --niter 1 \
266 --steps 2000 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05 \
267 --gax 1.5 --gay 1 --gaz 0.6
268```
270Copy it from under the stats line and compare the `ms/step` it reports with the
271app's. Both sides go through the one shared
272[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
273command, the benchmark parses it back — so there is no second copy of the
274defaults for the two runs to drift apart on. Geometry parameters take a `g`
275prefix (`--gwaist`) so a shape parameter can never collide with a model one.
277The app reports **two** numbers and only the first is comparable to the
278benchmark: `solver` is the batch of steps alone, waited for but not read back;
279`ms/frame` additionally carries a GPU→CPU readback per species, the
280colormapping, and the vertex upload. Those per-frame costs are fixed and do not
281shrink when the GPU gets faster, so on a quick GPU a frame can easily cost ten
282times the steps inside it. That is expected and is not the solver being slower
283in the browser.
285To attribute the gap rather than guess at it:
287```
288node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
289```
291measures the same solver work in both — batched, nothing read back, no rendering
292on either side — and reports each with its CPU-encoding share, the Fourier
293stage, and the adapter. It stops you first if the two are not even the same
294device, which is a common cause of "the browser is much slower". Both sides
295resolve the geometry and the iteration count from the same constants, because
296the iteration count is unrolled into the step and a mismatch would compare two
297different amounts of work.
299The app's **Benchmark** button runs the same measurement in the page, plus the
300**ramp** — the first third of the run against the last. GPUs downclock when
301idle and an animation-paced loop leaves them idle most of every frame, so a
302large ramp means the steady-state number is limited by clocks rather than work.
304### Is it really the same computation?
306```
307node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
308```
310runs one identical spec on the desktop and in a real browser and compares the
311final spectral state. The pipeline is deterministic given (model source,
312geometry, parameters, lmax, niter, seed, steps), so the two should agree to fp32
313round-off — not bit for bit, since GPUs differ in fused-multiply-add and other
314latitude fp32 allows. It also reports which Fourier stage each side chose, since
315FFT and DFT are genuinely different algorithms that round differently.
317Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), an
318optional dependency so that an unsupported platform fails the install of that
319package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
320`--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
321(`;`-separated) passes Dawn options through.
323## Tests
325There is no second implementation of the solver to diff against, so the `.m`
326path is checked against **closed-form answers** and against **exact structural
327properties**. Four modules, run in both environments:
329[`test/analyticChecks.ts`](test/analyticChecks.ts) — cases whose evolution is
330known exactly, run through the whole real pipeline. All three are statements
331about the round sphere, so all three build on the sphere geometry:
333- **A** — a linear reaction leaves every mode independent, growing by exactly
334 `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. Pins the transform round trip, the
335 eigenvalue mapping, the IMEX update and the state feedback at once. ~2e-7 over
336 20 steps.
337- **B** — a nonlinear reaction on a uniform field stays uniform, so each step is
338 exactly the scalar ODE map. 1.5e-8 over 25 steps.
339- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
340 linearized 2×2 IMEX recurrence, and `(l=24, m=7)` is confirmed unstable.
341 Looser (~4e-3) because fp32 keeps about four digits of a perturbation that
342 small.
344[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface and the loop:
346- every geometry compiles and closes; the sphere has radius 1 everywhere and is
347 **exactly degree 1** in the harmonics, which is what makes the reference case
348 exact rather than merely accurate;
349- the peanut matches its own closed-form radial profile at every grid point, and
350 **the same coefficients give the same surface on a 2× grid** — the 2× Gauss
351 latitudes share no point with the 1× ones, so agreeing there is agreeing
352 everywhere, which is what "rendered exactly, not subdivided" means;
353- unrolling is **exactly linear** in the trip count, and the state after 20 steps
354 is **bit-identical** at 0, 1 and 4 iterations;
355- a runtime loop bound is refused at compile time;
356- swapping the surface mid-run leaves the spectral state untouched.
358[`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
359and asserts **how many kernels it compiles to**, split into the base step and
360what one solve iteration adds. That is a fusion guard: if numbl's inline pass
361stops folding, the results stay correct while every operator becomes its own
362dispatch, which is invisible in the numbers.
364[`test/transformChecks.ts`](test/transformChecks.ts) compares the WGSL transforms
365against shtns-webgpu's f64 CPU twin.
367- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
368 `--skip-without-gpu` lets a machine without one say so and move on (which is
369 what CI does, since the browser suite covers the same modules).
370- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
371 Also runs the soak. A few geometry tolerances are set by SwiftShader's fp32,
372 which is about an order of magnitude looser than Dawn's.
374Other commands:
376- `npm run bench -- --help` — the desktop benchmark.
377- `npm run bench:sht -- --help` — the transforms alone, no solver.
378- `npx vite-node scripts/diagnose-sht.ts` — when the transform tests fail on a
379 GPU, say *which* stage is wrong.
380- `npx vite-node scripts/diagnose-leg.ts [--m 0]` — read the Legendre recurrence
381 out of the production shader term by term.
382- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
383 pattern saturates rather than decaying or diverging.
384- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
385 sampling JS heap and catching crashes.
386- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
387 demo after a number of steps.
388- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL.
389- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
90108a6Command for testing against a reference implementationOwen Melia 391### Testing against a reference implementation
393Reference solutions live in the sibling
394[turing-surface-test-data](https://github.com/concept-collection/turing-surface-test-data)
395repo, so an independently-written solver never has to depend on this one.
396`cases/schnakenberg-ellipsoid.md` there specifies the one case this repo
397currently ships a reference for;
398[`docs/ellipsoid-reference-spec.md`](docs/ellipsoid-reference-spec.md)
399restates it in this repo's own terms.
401`npm run ref -- --in <file>` (`scripts/ref.ts`) loads a reference file, runs
402the solver from its exact initial spectral state to the same physical end
403time, and reports the relative-L2 error against its final state (plus a
404geometry sanity check). `--niter` overrides the surface-correction iteration
405count independent of the file, and `--tolerance` turns the check into a
406pass/fail for CI.
410```
411npm install
412npm run dev # local dev server
413npm run build # type-check + production build to dist/
414```
416### The numbl dependency
418numbl is a local `file:../../numbl` dependency, so a sibling checkout of
419[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
420compiler internals — parser, lowerer, IR, inline pass — which its package
421`exports` map does not publish, so they are reached through the `numbl-src` path
422alias in [`vite.config.ts`](vite.config.ts).
424The exact surface we depend on is written down in
425[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
426*that*, not against numbl's sources. This keeps this project's compiler settings
427independent of numbl's, and means a change to one of those shapes upstream
428breaks the build here with a clear diff rather than deep inside numbl's tree.
429The `For` IR node is spelled out there, since the planner now walks it.
431CI clones numbl to the sibling path that the `file:` dependency expects, pinned
432to a commit, with `--ignore-scripts` (npm runs a linked package's `prepare`
433script, and numbl's is husky). numbl's own `node_modules` are not needed: the
434slice we import is self-contained TypeScript.
436The `scripts/*.ts` entry points that touch the compiler go through `vite-node`,
437so they resolve imports exactly as the browser build does. Plain `node` cannot:
438numbl's sources import each other as `./foo.js` while the files are `.ts`.
440Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
442## License
444CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).