1# turing-sphere
3Reaction–diffusion systems (Turing patterns) solved **live in the browser on the
4surface of a sphere**, using a spectral spherical-harmonic method with the
5transforms running on the GPU via WebGPU.
7The solver itself is **MATLAB**. The `.m` files under [`models/`](models/) are the
8algorithm — [numbl](https://numbl.org) parses and lowers them in the browser, and
9each element-wise line becomes a WebGPU compute kernel. You can edit the MATLAB
10on the page and watch the pattern change.
12**Live demo:** <https://concept-collection.github.io/turing-sphere-2/>
14## What it does
16It solves the N-species system
18```
19d(u_k)/dt = D_k*lap_s(u_k) + f_k(t, x, y, z, u_1, ..., u_N), k = 1, ..., N
20```
22on the unit sphere, where `lap_s` is the Laplace–Beltrami operator. Diffusion is
23treated implicitly in spherical-harmonic coefficient space, where `lap_s` is
24diagonal with eigenvalues `-l(l+1)`; reaction is treated explicitly on the grid.
25The two are combined with a first-order IMEX Euler step — the entire time loop is
27```
28V_k = synth(U_k) # spectral -> grid
29R_k = analys(f_k(t, x, y, z, V_1..V_N)) # reaction on grid -> spectral
30U_k = (U_k + dt*R_k) / (1 + dt*D_k*l(l+1))
31```
33You watch the patterns emerge in real time on orbitable 3D spheres (one per
34species, cameras synced), with pause/resume, re-seeding, live parameter editing,
35and colormap selection.
37Three models are included, one `.m` file each:
39- **[Schnakenberg](models/schnakenberg.m)** — Turing spots (unstable band
40 14 ≤ l ≤ 40, peak l = 24)
41- **[Brusselator](models/brusselator.m)** — stripes and spots from a stiffer reaction
42- **[Allen–Cahn](models/allencahn.m)** — a single species whose interfaces form
43 and coarsen
45## MATLAB, compiled to WebGPU
47A model file is ordinary MATLAB defining two functions — `init` builds the initial
48spectral state, `step` advances it one timestep:
50```matlab
51function [Un, Vn, u, v] = step(U, V, lam, a, b, D1, D2, dt)
52 u = synth(U);
53 v = synth(V);
54 uuv = u .* u .* v;
55 Un = (U + dt * analys(a - u + uuv)) ./ (1 + (dt * D1) * lam);
56 Vn = (V + dt * analys(b - uuv)) ./ (1 + (dt * D2) * lam);
57end
58```
60Getting from there to the GPU uses numbl for everything up to the IR, and this
61repo only for the backend:
631. **numbl parses and lowers.** Each function is specialized for the concrete
64 argument types of the current grid, via the same `specializeUserFunction`
65 entry point numbl's own JIT uses. Types and array shapes are fixed at this
66 point, so the backend never has to re-decide what an operation means.
672. **numbl's inline pass fuses.** Lowering emits one statement per *operator*
68 (ANF); `inlinePass` folds single-use temps back into their consumer, so one
69 line of MATLAB becomes one expression tree. `uuv = u .* u .* v` arrives as a
70 single statement, not three.
713. **This repo emits WGSL** ([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). Each
72 element-wise statement becomes one compute kernel that computes one output
73 element per invocation — the WebGPU counterpart of numbl's own C-side fused
74 emitter. Anything it cannot express is refused at compile time with a source
75 position, never silently mis-compiled.
764. **`synth` / `analys` are external operations.** numbl learns their type rules
77 from a `.mtoc2.js` workspace file — its sanctioned extension point for a
78 JS-defined builtin — and the backend maps each call onto the existing
79 spherical-harmonic compute pipelines.
81The Schnakenberg step above compiles to 11 GPU operations: 4 transforms, 5
82generated kernels, and 2 buffer copies feeding the new state back.
84Two consequences worth noting:
86- **The step is synchronous.** WebGPU's encode path (`writeBuffer`, dispatch,
87 `submit`) is all synchronous; only readback and pipeline creation are async, and
88 every pipeline is built once at compile time. So a timestep is pure command
89 recording — the whole batch goes out in one submit, and the only `await` in the
90 loop is the single readback per rendered frame. numbl's own execution being
91 synchronous is therefore not an obstacle: nothing about the algorithm needs to
92 block.
93- **Parameters are uniforms, not constants.** Tunable scalars are deliberately
94 lowered without exact values, so moving a slider rewrites a small buffer
95 instead of triggering a recompile. Editing the MATLAB recompiles; changing `dt`
96 does not.
98## Provenance
100This is the browser port of a MATLAB reference implementation
101(`SphericalReactionDiffusion.m`, "websph"), which defines the solver through a
102four-member porting boundary: `coeffs2vals`, `vals2coeffs`, `grid.lat`,
103`grid.lon`. Profiling of the MATLAB version shows the transforms are ~96% of
104compute, so this port swaps in:
106- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
107 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
108 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Its source is vendored under
109 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference
110 transform used for testing.
111- **Rendering:** three.js spheres with per-vertex colormaps, adapted from the
112 `SphereEmbedding` view in
113 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
114 extension package ([`src/render/`](src/render/)).
115- **Solver:** the MATLAB stayed MATLAB. [`models/`](models/) holds the IMEX loop
116 as `.m` files, executed on the GPU by [`src/mgpu/`](src/mgpu/). There is no
117 second implementation: the app, the desktop benchmark and the tests all compile
118 and run the same `.m`.
120An earlier version of this repo carried a TypeScript port of the loop alongside
121the `.m`, and used it as the test oracle. That is gone. Two implementations
122agreeing only shows they share assumptions, so the `.m` path is now checked
123against closed-form answers instead — see [Tests](#tests). The one place a second
124implementation is still the right oracle is the transforms themselves, where
125[`src/sht/reference.ts`](src/sht/reference.ts) is shtns-webgpu's own f64
126direct-summation twin.
128Because the algorithm is compiled to compute shaders, **WebGPU is required** —
129there is no CPU fallback (the f64 CPU transform remains, for tests).
131## Numerics
133- Grid: Gauss–Legendre × equispaced-phi, dealiased for the cubic reactions with
134 the `(pdeg+1)` rule from the reference implementation:
135 `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1` (rounded up to a power
136 of two for the GPU FFT path). At the default lmax 63 that is a 128×256 grid.
137- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
138 coefficients for m ≥ 0, m-major ordering.
139- fp32 transforms introduce ~1e-6 relative error per step (verified against the
140 f64 CPU path); for pattern formation from 1e-2 seeded noise this is
141 inconsequential.
143## Desktop vs browser
145How much does running this in a browser cost? [`scripts/bench.ts`](scripts/bench.ts)
146runs the *same* thing — same `.m`, lowered by numbl into the same WGSL kernels,
147over the same transforms — from Node on desktop WebGPU (Google Dawn), and the app
148prints the command line that reproduces whatever it is currently simulating:
150```
151npm run bench -- --preset schnak-spots --lmax 63 --steps 2000 \
152 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05
153```
155Copy it from under the stats line, run it, and compare the `ms/step` it reports
156with the app's. Both sides go through the one shared
157[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
158command, the benchmark parses it back — so there is no second copy of the
159defaults for the two runs to drift apart on. Both then go through the same
160[`ModelSession`](src/mgpu/session.ts), down to the device request in
161`requestShtDevice()` (Dawn is installed under `navigator.gpu` and the WebGPU
162globals, and the rest runs unchanged).
164The benchmark runs under `vite-node`, which is what resolves numbl's compiler
165sources and the `?raw` model imports — plain Node cannot (see
166[The numbl dependency](#the-numbl-dependency)).
168It reports two numbers, because they answer different questions:
170```
171 0.54 ms/step 1857.5 steps/s 92.87 model time/s (batches of 16)
172 one step per submit: 0.74 ms mean · median 0.60 · p05 0.51 · p95 1.29 · min 0.50
173```
175The first is throughput: a batch of steps submitted together and awaited once,
176which is how the app runs and what keeping the state in GPU buffers is for. The
177second is per-step latency, one submit each — comparable to a design that
178synchronises every step, and the only way to get a distribution.
180**What the GPU-resident design is worth.** At lmax 31 on an Intel Xe (Mesa, via
181Dawn) this path runs at **0.25 ms/step**, against **3.01 ms/step** for the
182TypeScript solver this repo used to carry — same machine, same transforms, same
183parameters. A **~12x** difference, and almost all of it is the four per-step
184buffer readbacks that version paid and this one does not. Note that CI, which
185only has a software rasterizer, shows no such gap: there the transforms dominate
186and both designs land within ~10% of each other. The saving is real but it is a
187saving on driver round-trips, so it only appears once the GPU is fast.
189Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), listed
190as an optional dependency so that a platform it has no binaries for fails the
191install of that package alone rather than the whole tree. `npm install` picks it
192up; without it there is no desktop GPU to run on and the benchmark says so.
193Those binaries need glibc 2.29+, which rules out older cluster images
194(RHEL/Rocky 8 is 2.28) unless you run inside a container with a newer base. Other
195flags: `--steps`, `--warmup`, `--batch`, `--json`, `--help`;
196`DAWN_FLAGS='backend=vulkan'` (`;`-separated) passes Dawn options through, e.g. to
197pick a backend or to compare against Dawn's own software adapter.
199### Comparing the two honestly
201The app reports **two** numbers, and only the first is comparable to the
202benchmark:
204```
205solver 0.58 ms/step (1724 steps/s) · 12.4 ms/frame incl. readback + render
206```
208`solver` is the batch of steps alone, waited for but not read back — the same
209thing the benchmark's throughput number measures. `ms/frame` additionally carries
210a GPU→CPU readback **per species**, the colormapping, and the vertex upload.
212Those per-frame costs are fixed: they do not shrink when the GPU gets faster. So
213the faster your GPU, the larger the ratio between them — on a quick discrete GPU
214it is easy for a frame to cost ten times the four steps inside it, purely because
215a `mapAsync` round trip in a browser has to drain the queue and cross into the GPU
216process. **That is expected, and it is not the solver being slower in the
217browser.** Compare `solver` with the benchmark's throughput line; comparing
218`ms/frame` against it measures the readback, not the computation.
220Other things the comparison does not control for:
222- the browser's renderer→GPU-process boundary on every submit, where Dawn in Node
223 is in-process; and, for a page that is not cross-origin isolated, coarser
224 `performance.now()`.
225- both sides are fp32 throughout, on the same generated kernels, so nothing here
226 is a numerics comparison — only a cost one.
228### Why the browser is slower, and how to find out by how much
230Some gap is real and some is measurement. Four numbers, in increasing order of
231what they include — walk down them and the gap attributes itself:
233```
234node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
235```
237measures the same solver work in both — batched, nothing read back, no rendering
238on either side — and reports each with its CPU-encoding share, the Fourier stage,
239and the adapter. It stops you first if the two are not even the same device: a
240browser quietly falling back to a software adapter is a common cause of "the
241browser is much slower", and then the ratio compares different hardware and means
242nothing.
244Or press **Benchmark** in the app: it pauses rendering and runs batches
245continuously for two seconds, reporting the same measurement the terminal makes,
246plus the **ramp** — the first third of the run against the last. GPUs downclock
247when idle and an animation-paced loop leaves them idle most of every frame, so a
248large ramp means the steady-state number is limited by clocks rather than by the
249work.
251By hand, five numbers, in increasing order of what they include:
253| number | includes |
254|---|---|
255| `npm run bench -- --lmax 63` | desktop solver: batched steps, one sync per batch, in-process Dawn |
256| the app's **Benchmark** button | browser solver, sustained, no rendering, no pacing |
257| `test.html?soak=2000&lmax=63` → `solver` | the same, without the page around it |
258| the app's `solver` | browser solver, one batch of 32 every two seconds |
259| the app's `ms/frame` | four steps **plus** a readback per species, colormapping and the vertex upload |
261If the soak matches the benchmark, the solver is fine in the browser and
262everything above it is readback and rendering. If the soak is itself slower, the
263remaining suspects are:
265- **the GPU-process boundary.** Every submit and every sync is IPC out of the
266 renderer; Dawn in Node is in-process. This is a fixed per-batch cost, so it hurts
267 most when the GPU is fast. `npm run bench -- --batch 4` makes the desktop pay a
268 sync as often as the app's frame loop does, which shows how much of the gap is
269 just amortization.
270- **not CPU command encoding**, which is worth ruling out explicitly because it is
271 the obvious suspect: a step is ~47 WebGPU calls, and 32 of them per burst is a
272 lot of JS→GPU traffic. Measured, it goes the other way — 0.009 ms/step in Chrome
273 against 0.062 ms/step under node-webgpu, because Chrome defers commands to the
274 GPU process while node-webgpu validates them inline. Encoding is *cheaper* in
275 the browser. Both `compare-perf.mjs` and the benchmark print it.
276- **competing with the renderer.** The page draws two spheres through WebGL on the
277 same GPU, in its own animation loop. The soak has no renderer, so comparing the
278 soak against the app's `solver` separates contention from everything else.
279- **clocks.** An animation-paced loop leaves the GPU idle for most of each 16 ms
280 frame, so it may never leave its low-power state, while the benchmark hammers it
281 continuously and boosts. On a thermally managed laptop this alone can be worth a
282 factor of two, and it is not something the code can fix. The **Benchmark**
283 button's ramp figure measures it directly.
284- **anything else using the GPU.** Another process competing for it changes
285 whichever run overlaps it, which makes a comparison across two separate
286 invocations meaningless. `compare-perf.mjs` runs both sides back to back in one
287 invocation partly for this reason.
288- **not buffer robustness**, another plausible suspect: WebGPU clamps every array
289 access for safety, which could cost real time in the transform kernels' inner
290 loops. Measured with Dawn's `disable_robustness` toggle
291 (`DAWN_FLAGS='enable-dawn-features=disable_robustness' npm run bench`), it makes
292 no difference here at all — 0.59 ms/step either way.
293- **which browser.** WebGPU implementations differ substantially in maturity;
294 Chrome and Safari are not interchangeable for this.
296None of these change *what* is computed — see below for how to confirm that
297independently.
299### Is it really the same computation?
301```
302node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
303```
305runs one identical spec on the desktop and in a real browser and compares the
306final spectral state. The pipeline is deterministic given (model source,
307parameters, lmax, seed, steps) — a seeded PRNG, then fixed arithmetic — so the two
308should agree to fp32 round-off. Both sides build their spec through the same
309`parseArgs`, so neither can quietly use a different default.
311They will *not* agree bit for bit; GPUs differ in fused-multiply-add and other
312latitude fp32 allows. Between Intel Xe (via Dawn) and SwiftShader — about as
313different as two implementations get — 200 steps at lmax 31 agree to a relative
314L2 of **2e-6**.
316It also reports which **Fourier stage** each side chose. `ShtPlan` picks FFT or
317DFT from the device's workgroup-storage and invocation limits, and those are
318genuinely different algorithms that round differently, so a mismatch there
319explains a difference in the values rather than being a symptom of one. The app's
320stats line and the benchmark both print the chosen stage for the same reason.
322## Tests
324There is no second implementation of the solver to diff against, so the `.m` path
325is checked against **closed-form answers**. Each case is one whose evolution is
326known exactly, run through the whole real pipeline — MATLAB source, numbl
327lowering, generated WGSL, GPU transforms — and compared with arithmetic
328([`test/analyticChecks.ts`](test/analyticChecks.ts)):
330- **A** — a linear reaction `f(u) = c*u` leaves every spherical-harmonic mode
331 independent, growing by exactly `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. This
332 pins the transform round-trip, the eigenvalue mapping, the IMEX update and the
333 state feedback at once, and checks that nothing leaks between modes. Agrees to
334 ~2e-7 over 20 steps.
335- **B** — a nonlinear reaction on a *uniform* field stays uniform and diffusion
336 cannot touch it, so each step is exactly the scalar ODE map. Agrees to 1.5e-8
337 over 25 steps. Checks that a generated kernel evaluates a nonlinear reaction.
338- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
339 linearized 2x2 IMEX recurrence, and the `(l=24, m=7)` mode is confirmed
340 unstable. Looser (~2e-3) because fp32 keeps only about four digits of a
341 perturbation that small.
343Two test models exist only for this: [`test/models/linear.m`](test/models/linear.m)
344and [`test/models/logistic.m`](test/models/logistic.m).
346Alongside those, [`test/modelChecks.ts`](test/modelChecks.ts) compiles every model
347the app offers and asserts **how many kernels it compiles to**. That is a fusion
348guard: numbl's lowering emits one statement per *operator* and its inline pass
349folds them back into per-line expression trees, and if that stops happening the
350results stay correct while every operator becomes its own dispatch. It is
351invisible in the numbers, so it is asserted directly. (It has already caught one
352regression.)
354[`test/transformChecks.ts`](test/transformChecks.ts) is the one remaining
355implementation-vs-implementation check, comparing the WGSL transforms against
356shtns-webgpu's f64 CPU twin.
358All three modules run in **both** environments, so the two GPU stacks get the same
359guarantees:
361- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
362 pass `--skip-without-gpu` to let a machine without one say so and move on
363 (which is what CI does, since the browser suite covers the same modules).
364- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
365 Also runs the soak.
367Other commands:
369- `npm run bench -- --help` — the desktop benchmark (see
370 [Desktop vs browser](#desktop-vs-browser)).
371- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
372 pattern saturates into O(1)-contrast spots rather than decaying or diverging.
373- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
374 sampling JS heap and catching crashes.
375- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
376 demo after a number of steps.
377- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL in a real
378 browser: load, press Run, confirm the solver advances.
379- `node scripts/compare-env.mjs` — run one identical spec on the desktop and in a
380 browser and compare the final state (see
381 [Is it really the same computation?](#is-it-really-the-same-computation)).
382- `node scripts/compare-perf.mjs` — measure the same solver work in both and split
383 the difference (see
384 [Why the browser is slower](#why-the-browser-is-slower-and-how-to-find-out-by-how-much)).
385- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
387### A note on canvas resizing
389Early long runs killed the browser after ~700–800 steps. The cause was the
390colorbar's min/max labels changing width as their digit count changed, which
391reflowed the panel, fired the `ResizeObserver`, and called
392`renderer.setSize()` — reallocating the WebGL drawing buffer. Assigning
393`canvas.width` also blanks the canvas even when the value is unchanged, so the
394same bug caused visible flicker. Fixed by giving the colorbar column a fixed
395width and making `SphereScene.resize()` return early on no-op resizes.
397## Development
399```
400npm install
401npm run dev # local dev server
402npm run build # type-check + production build to dist/
403```
405### The numbl dependency
407numbl is a local `file:../../numbl` dependency, so a sibling checkout of
408[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
409compiler internals — parser, lowerer, IR, inline pass — which its package
410`exports` map does not publish, so they are reached through the `numbl-src` path
411alias in [`vite.config.ts`](vite.config.ts).
413The exact surface we depend on is written down in
414[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
415*that*, not against numbl's sources. This keeps this project's compiler settings
416independent of numbl's (its sources do not type-check under the stricter options
417used here), and means a change to one of those shapes upstream breaks the build
418here with a clear diff rather than deep inside numbl's tree.
420The compiler is ~395 kB gzipped and lands in its own chunk. That is the cost of
421compiling MATLAB in the page; a build-time lowering step could remove it at the
422price of no longer being editable live.
424CI clones numbl to the sibling path that `file:` dependency expects, pinned to a
425commit. Two details make that work, both verified by building against a checkout
426that had none of numbl's own dependencies installed:
428- **numbl's `node_modules` are not needed.** The slice we import — parser,
429 lowering, IR, inline pass — is self-contained TypeScript. (Other parts of numbl
430 do import `three`, `react` and `fflate`; we never reach them.)
431- **the install must pass `--ignore-scripts`.** npm runs a linked package's
432 `prepare` script, and numbl's is `husky`, which is not installed in CI.
434The `scripts/*.ts` entry points that touch the compiler (the benchmark, the node
435tests, the long run) go through `vite-node`, so they resolve imports exactly as the
436browser build does — the `numbl-src` alias and the `?raw` model imports included.
437Plain `node` cannot: numbl's sources import each other as `./foo.js` while the
438files are `.ts`, which needs a bundler's resolution. Scripts that do not touch the
439compiler (`soak.mjs`, `screenshot.mjs`, `check-live.mjs`, `test-gpu.mjs`) are plain
440`.mjs` and run under `node` directly.
442Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
444## License
446CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).