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/>
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/).
118[`src/solver/`](src/solver/) still holds the TypeScript port of the same loop. The
119app no longer runs it, but it is an *independent* implementation of the scheme,
120which makes it the test oracle: `npm run test:gpu` runs both from the same seeded
121perturbation through the same transforms and compares. It is also where parameter
122metadata (names, defaults, slider ranges) lives, so the two paths cannot be
123configured differently.
125Because the algorithm is now compiled to compute shaders, **WebGPU is required** —
126there is no CPU fallback in the app (the f64 CPU transform remains, for tests).
128## Numerics
130- Grid: Gauss–Legendre × equispaced-phi, dealiased for the cubic reactions with
131 the `(pdeg+1)` rule from the reference implementation:
132 `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1` (rounded up to a power
133 of two for the GPU FFT path). At the default lmax 63 that is a 128×256 grid.
134- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
135 coefficients for m ≥ 0, m-major ordering.
136- fp32 transforms introduce ~1e-6 relative error per step (verified against the
137 f64 CPU path); for pattern formation from 1e-2 seeded noise this is
138 inconsequential.
140## Desktop vs browser
142How much does running this in a browser cost? [`scripts/bench.ts`](scripts/bench.ts)
143runs the reference solver — same WGSL transforms, same parameters — from Node on
144desktop WebGPU (Google Dawn), and the app prints the command line that
145reproduces whatever it is currently simulating:
147```
148node scripts/bench.mjs --preset schnak-spots --lmax 63 --backend webgpu --steps 2000 \
149 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05
150```
152Copy it from under the stats line, run it, and compare the `ms/step` it reports
153with the app's. Both sides go through the one shared
154[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
155command, the benchmark parses it back — so there is no second copy of the
156defaults for the two runs to drift apart on. Node runs the TypeScript sources
157directly, so `src/` is literally the same code in both places, down to the
158device request in `requestShtDevice()` (Dawn is installed under `navigator.gpu`
159and the WebGPU globals, and the rest runs unchanged).
161Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), listed
162as an optional dependency so that a platform it has no binaries for fails the
163install of that package alone rather than the whole tree. `npm install` picks it
164up; without it, only `--backend cpu` runs and the benchmark says so. Those
165binaries need glibc 2.29+, which rules out older cluster images (RHEL/Rocky 8 is
1662.28) unless you run inside a container with a newer base. Other
167flags: `--steps`, `--warmup`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
168(`;`-separated) passes Dawn options through, e.g. to pick a backend or to
169compare against Dawn's own software adapter.
171What the comparison does and does not control for:
173- **it is not the same solver.** The benchmark runs the TypeScript reference; the
174 app runs the `.m` compiled to WGSL. Node cannot load numbl's TypeScript sources
175 (its internal imports are extensionless-`.js`, which needs a bundler's
176 resolution), so the `.m` path is browser-only for now. Same scheme, same
177 transforms, same parameters — but the reaction and the IMEX update happen in f64
178 on the CPU there and in fp32 on the GPU here.
179- the benchmark is **solver only**; the app's `ms/step` includes the per-frame
180 readback amortized over the step batch. For a browser number with no rendering,
181 open `test.html?soak=2000&lmax=63` (that soak also runs the reference solver).
182- the reference pays a buffer readback on *every* transform — four driver
183 round-trips per step — so it measures submit-and-map latency more than
184 arithmetic. The `.m` path keeps everything in GPU buffers and submits once per
185 batch, which is where its advantage should come from. On the software rasterizer
186 in CI the transforms dominate and the two come out within ~10% of each other;
187 the gap on real hardware is untested.
188- the browser adds its own GPU-process boundary and, for a page that is not
189 cross-origin isolated, coarser timers.
191## Tests
193- `npm run bench -- --help` — the desktop benchmark above (see
194 [Desktop vs browser](#desktop-vs-browser)).
195- `npm run test:node` — f64 solver correctness in Node: exact single-mode
196 linear recurrence, exact uniform-state reaction ODE, and the linearized
197 Turing-mode 2×2 IMEX recurrence (all at ~1e-12).
198- `npm run test:gpu` — builds and drives headless Chrome: GPU-vs-CPU transform
199 and solver cross-checks, a 100-step stability run, and then for every `.m`
200 model: that it compiles, that its element-wise lines each fuse into exactly one
201 kernel, and that 10 steps agree with the reference solver from the same seed
202 (they agree to ~1e-7 relative L2 — fp32 round-off).
203- `node scripts/longrun-node.ts` — CPU run to t = 100 confirming pattern
204 saturation.
205- `node scripts/soak.mjs [steps] [lmax] [backend]` — drive the demo for many
206 steps, sampling JS heap and catching crashes. A 900-step run at lmax 63 on
207 software WebGPU (SwiftShader) completes with a flat ~4 MB heap.
208- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot
209 the demo after a number of steps.
210- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL in a real
211 browser: load, press Run, confirm the solver advances.
212- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
214### A note on canvas resizing
216Early long runs killed the browser after ~700–800 steps. The cause was the
217colorbar's min/max labels changing width as their digit count changed, which
218reflowed the panel, fired the `ResizeObserver`, and called
219`renderer.setSize()` — reallocating the WebGL drawing buffer. Assigning
220`canvas.width` also blanks the canvas even when the value is unchanged, so the
221same bug caused visible flicker. Fixed by giving the colorbar column a fixed
222width and making `SphereScene.resize()` return early on no-op resizes.
224## Development
226```
227npm install
228npm run dev # local dev server
229npm run build # type-check + production build to dist/
230```
232### The numbl dependency
234numbl is a local `file:../../numbl` dependency, so a sibling checkout of
235[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
236compiler internals — parser, lowerer, IR, inline pass — which its package
237`exports` map does not publish, so they are reached through the `numbl-src` path
238alias in [`vite.config.ts`](vite.config.ts).
240The exact surface we depend on is written down in
241[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
242*that*, not against numbl's sources. This keeps this project's compiler settings
243independent of numbl's (its sources do not type-check under the stricter options
244used here), and means a change to one of those shapes upstream breaks the build
245here with a clear diff rather than deep inside numbl's tree.
247The compiler is ~395 kB gzipped and lands in its own chunk. That is the cost of
248compiling MATLAB in the page; a build-time lowering step could remove it at the
249price of no longer being editable live.
251CI clones numbl to the sibling path that `file:` dependency expects, pinned to a
252commit. Two details make that work, both verified by building against a checkout
253that had none of numbl's own dependencies installed:
255- **numbl's `node_modules` are not needed.** The slice we import — parser,
256 lowering, IR, inline pass — is self-contained TypeScript. (Other parts of numbl
257 do import `three`, `react` and `fflate`; we never reach them.)
258- **the install must pass `--ignore-scripts`.** npm runs a linked package's
259 `prepare` script, and numbl's is `husky`, which is not installed in CI.
261The `.ts` entry points under `scripts/` are run by Node directly, which strips
262types without being asked only from Node 22.18 / 23.6 / 24 on. Everything here
263works back to 22.6, where stripping exists but is flagged: the npm scripts pass
264`--experimental-strip-types` themselves, and the benchmark — the one command
265that gets copied to other machines — goes through
266[`scripts/bench.mjs`](scripts/bench.mjs), which re-runs itself with the flag
267when it has to. Invoking a `scripts/*.ts` file by hand on 22.6–22.17 needs the
268flag spelled out.
270Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
272## License
274CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).