2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 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.
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 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.
c47b479Point the demo link back at the original repo's Pages siteJeremy Magland 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,
f168295Movie export: recompute the run from t = 0 into a captioned MP4Jeremy Magland 35colormap selection, and movie download — the run is recomputed from t = 0 and
36encoded to a captioned MP4 in the browser (WebCodecs). The display can
37oversample the solver — the state is spectral, so evaluating it on a finer grid
38for rendering is exact interpolation, not smoothing. This never touches the
39solver or its grid; by default it turns on only when the solver grid is coarse.
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 41Three models are included, one `.m` file each:
43- **[Schnakenberg](models/schnakenberg.m)** — Turing spots (unstable band
44 14 ≤ l ≤ 40, peak l = 24)
45- **[Brusselator](models/brusselator.m)** — stripes and spots from a stiffer reaction
46- **[Allen–Cahn](models/allencahn.m)** — a single species whose interfaces form
47 and coarsen
49## MATLAB, compiled to WebGPU
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 51A model file is ordinary MATLAB defining two functions — `init` builds the initial
52spectral state, `step` advances it one timestep:
54```matlab
55function [Un, Vn, u, v] = step(U, V, lam, a, b, D1, D2, dt)
56 u = synth(U);
57 v = synth(V);
58 uuv = u .* u .* v;
59 Un = (U + dt * analys(a - u + uuv)) ./ (1 + (dt * D1) * lam);
60 Vn = (V + dt * analys(b - uuv)) ./ (1 + (dt * D2) * lam);
61end
62```
64Getting from there to the GPU uses numbl for everything up to the IR, and this
65repo only for the backend:
671. **numbl parses and lowers.** Each function is specialized for the concrete
68 argument types of the current grid, via the same `specializeUserFunction`
69 entry point numbl's own JIT uses. Types and array shapes are fixed at this
70 point, so the backend never has to re-decide what an operation means.
712. **numbl's inline pass fuses.** Lowering emits one statement per *operator*
72 (ANF); `inlinePass` folds single-use temps back into their consumer, so one
73 line of MATLAB becomes one expression tree. `uuv = u .* u .* v` arrives as a
74 single statement, not three.
753. **This repo emits WGSL** ([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). Each
76 element-wise statement becomes one compute kernel that computes one output
77 element per invocation — the WebGPU counterpart of numbl's own C-side fused
78 emitter. Anything it cannot express is refused at compile time with a source
79 position, never silently mis-compiled.
804. **`synth` / `analys` are external operations.** numbl learns their type rules
81 from a `.mtoc2.js` workspace file — its sanctioned extension point for a
82 JS-defined builtin — and the backend maps each call onto the existing
83 spherical-harmonic compute pipelines.
85The Schnakenberg step above compiles to 11 GPU operations: 4 transforms, 5
86generated kernels, and 2 buffer copies feeding the new state back.
88Two consequences worth noting:
90- **The step is synchronous.** WebGPU's encode path (`writeBuffer`, dispatch,
91 `submit`) is all synchronous; only readback and pipeline creation are async, and
92 every pipeline is built once at compile time. So a timestep is pure command
93 recording — the whole batch goes out in one submit, and the only `await` in the
94 loop is the single readback per rendered frame. numbl's own execution being
95 synchronous is therefore not an obstacle: nothing about the algorithm needs to
96 block.
97- **Parameters are uniforms, not constants.** Tunable scalars are deliberately
98 lowered without exact values, so moving a slider rewrites a small buffer
99 instead of triggering a recompile. Editing the MATLAB recompiles; changing `dt`
100 does not.
102## Provenance
104This is the browser port of a MATLAB reference implementation
105(`SphericalReactionDiffusion.m`, "websph"), which defines the solver through a
106four-member porting boundary: `coeffs2vals`, `vals2coeffs`, `grid.lat`,
107`grid.lon`. Profiling of the MATLAB version shows the transforms are ~96% of
108compute, so this port swaps in:
110- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
111 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
112 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Its source is vendored under
113 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 114 transform used for testing. [`bench/shtns/`](bench/shtns/) builds the real
115 SHTNS and measures ours against it — see
116 [Against upstream SHTNS](#against-upstream-shtns).
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 117- **Rendering:** three.js spheres with per-vertex colormaps, adapted from the
118 `SphereEmbedding` view in
119 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
120 extension package ([`src/render/`](src/render/)).
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 121- **Solver:** the MATLAB stayed MATLAB. [`models/`](models/) holds the IMEX loop
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 122 as `.m` files, executed on the GPU by [`src/mgpu/`](src/mgpu/). There is no
123 second implementation: the app, the desktop benchmark and the tests all compile
124 and run the same `.m`.
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 126An earlier version of this repo carried a TypeScript port of the loop alongside
127the `.m`, and used it as the test oracle. That is gone. Two implementations
128agreeing only shows they share assumptions, so the `.m` path is now checked
129against closed-form answers instead — see [Tests](#tests). The one place a second
130implementation is still the right oracle is the transforms themselves, where
131[`src/sht/reference.ts`](src/sht/reference.ts) is shtns-webgpu's own f64
132direct-summation twin.
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 134Because the algorithm is compiled to compute shaders, **WebGPU is required** —
135there is no CPU fallback (the f64 CPU transform remains, for tests).
137## Numerics
139- Grid: Gauss–Legendre × equispaced-phi, dealiased for the cubic reactions with
140 the `(pdeg+1)` rule from the reference implementation:
141 `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1` (rounded up to a power
142 of two for the GPU FFT path). At the default lmax 63 that is a 128×256 grid.
143- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
144 coefficients for m ≥ 0, m-major ordering.
145- fp32 transforms introduce ~1e-6 relative error per step (verified against the
146 f64 CPU path); for pattern formation from 1e-2 seeded noise this is
147 inconsequential.
15a77e2Add a desktop WebGPU benchmark and show its command in the appJeremy Magland 149## Desktop vs browser
151How much does running this in a browser cost? [`scripts/bench.ts`](scripts/bench.ts)
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 152runs the *same* thing — same `.m`, lowered by numbl into the same WGSL kernels,
153over the same transforms — from Node on desktop WebGPU (Google Dawn), and the app
154prints the command line that reproduces whatever it is currently simulating:
156```
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 157npm run bench -- --preset schnak-spots --lmax 63 --steps 2000 \
15a77e2Add a desktop WebGPU benchmark and show its command in the appJeremy Magland 158 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05
159```
161Copy it from under the stats line, run it, and compare the `ms/step` it reports
162with the app's. Both sides go through the one shared
163[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
164command, the benchmark parses it back — so there is no second copy of the
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 165defaults for the two runs to drift apart on. Both then go through the same
166[`ModelSession`](src/mgpu/session.ts), down to the device request in
167`requestShtDevice()` (Dawn is installed under `navigator.gpu` and the WebGPU
168globals, and the rest runs unchanged).
170The benchmark runs under `vite-node`, which is what resolves numbl's compiler
171sources and the `?raw` model imports — plain Node cannot (see
172[The numbl dependency](#the-numbl-dependency)).
174It reports two numbers, because they answer different questions:
176```
177 0.54 ms/step 1857.5 steps/s 92.87 model time/s (batches of 16)
178 one step per submit: 0.74 ms mean · median 0.60 · p05 0.51 · p95 1.29 · min 0.50
179```
181The first is throughput: a batch of steps submitted together and awaited once,
182which is how the app runs and what keeping the state in GPU buffers is for. The
183second is per-step latency, one submit each — comparable to a design that
184synchronises every step, and the only way to get a distribution.
186**What the GPU-resident design is worth.** At lmax 31 on an Intel Xe (Mesa, via
187Dawn) this path runs at **0.25 ms/step**, against **3.01 ms/step** for the
188TypeScript solver this repo used to carry — same machine, same transforms, same
189parameters. A **~12x** difference, and almost all of it is the four per-step
190buffer readbacks that version paid and this one does not. Note that CI, which
191only has a software rasterizer, shows no such gap: there the transforms dominate
192and both designs land within ~10% of each other. The saving is real but it is a
193saving on driver round-trips, so it only appears once the GPU is fast.
e5b7827Fix CI: do not omit optional dependenciesJeremy Magland 195Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), listed
196as an optional dependency so that a platform it has no binaries for fails the
197install of that package alone rather than the whole tree. `npm install` picks it
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 198up; without it there is no desktop GPU to run on and the benchmark says so.
199Those binaries need glibc 2.29+, which rules out older cluster images
200(RHEL/Rocky 8 is 2.28) unless you run inside a container with a newer base. Other
201flags: `--steps`, `--warmup`, `--batch`, `--json`, `--help`;
202`DAWN_FLAGS='backend=vulkan'` (`;`-separated) passes Dawn options through, e.g. to
203pick a backend or to compare against Dawn's own software adapter.
0f3abbdSeparate solver time from frame time, and add a cross-environment checkJeremy Magland 205### Comparing the two honestly
0f3abbdSeparate solver time from frame time, and add a cross-environment checkJeremy Magland 207The app reports **two** numbers, and only the first is comparable to the
208benchmark:
210```
211solver 0.58 ms/step (1724 steps/s) · 12.4 ms/frame incl. readback + render
212```
214`solver` is the batch of steps alone, waited for but not read back — the same
215thing the benchmark's throughput number measures. `ms/frame` additionally carries
216a GPU→CPU readback **per species**, the colormapping, and the vertex upload.
218Those per-frame costs are fixed: they do not shrink when the GPU gets faster. So
219the faster your GPU, the larger the ratio between them — on a quick discrete GPU
220it is easy for a frame to cost ten times the four steps inside it, purely because
221a `mapAsync` round trip in a browser has to drain the queue and cross into the GPU
222process. **That is expected, and it is not the solver being slower in the
223browser.** Compare `solver` with the benchmark's throughput line; comparing
224`ms/frame` against it measures the readback, not the computation.
226Other things the comparison does not control for:
228- the browser's renderer→GPU-process boundary on every submit, where Dawn in Node
229 is in-process; and, for a page that is not cross-origin isolated, coarser
230 `performance.now()`.
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 231- both sides are fp32 throughout, on the same generated kernels, so nothing here
232 is a numerics comparison — only a cost one.
cde22eaStop the per-frame sync from inflating the app's solver numberJeremy Magland 234### Why the browser is slower, and how to find out by how much
236Some gap is real and some is measurement. Four numbers, in increasing order of
237what they include — walk down them and the gap attributes itself:
240node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
241```
243measures the same solver work in both — batched, nothing read back, no rendering
244on either side — and reports each with its CPU-encoding share, the Fourier stage,
245and the adapter. It stops you first if the two are not even the same device: a
246browser quietly falling back to a software adapter is a common cause of "the
247browser is much slower", and then the ratio compares different hardware and means
248nothing.
3221abdAdd a Benchmark button that measures the solver and the GPU's clock rampJeremy Magland 250Or press **Benchmark** in the app: it pauses rendering and runs batches
251continuously for two seconds, reporting the same measurement the terminal makes,
252plus the **ramp** — the first third of the run against the last. GPUs downclock
253when idle and an animation-paced loop leaves them idle most of every frame, so a
254large ramp means the steady-state number is limited by clocks rather than by the
255work.
257By hand, five numbers, in increasing order of what they include:
cde22eaStop the per-frame sync from inflating the app's solver numberJeremy Magland 259| number | includes |
260|---|---|
261| `npm run bench -- --lmax 63` | desktop solver: batched steps, one sync per batch, in-process Dawn |
3221abdAdd a Benchmark button that measures the solver and the GPU's clock rampJeremy Magland 262| the app's **Benchmark** button | browser solver, sustained, no rendering, no pacing |
263| `test.html?soak=2000&lmax=63` → `solver` | the same, without the page around it |
264| the app's `solver` | browser solver, one batch of 32 every two seconds |
cde22eaStop the per-frame sync from inflating the app's solver numberJeremy Magland 265| the app's `ms/frame` | four steps **plus** a readback per species, colormapping and the vertex upload |
267If the soak matches the benchmark, the solver is fine in the browser and
268everything above it is readback and rendering. If the soak is itself slower, the
269remaining suspects are:
271- **the GPU-process boundary.** Every submit and every sync is IPC out of the
272 renderer; Dawn in Node is in-process. This is a fixed per-batch cost, so it hurts
273 most when the GPU is fast. `npm run bench -- --batch 4` makes the desktop pay a
274 sync as often as the app's frame loop does, which shows how much of the gap is
275 just amortization.
17db8f1Add scripts/compare-perf.mjs, and rule out CPU command encodingJeremy Magland 276- **not CPU command encoding**, which is worth ruling out explicitly because it is
277 the obvious suspect: a step is ~47 WebGPU calls, and 32 of them per burst is a
278 lot of JS→GPU traffic. Measured, it goes the other way — 0.009 ms/step in Chrome
279 against 0.062 ms/step under node-webgpu, because Chrome defers commands to the
280 GPU process while node-webgpu validates them inline. Encoding is *cheaper* in
281 the browser. Both `compare-perf.mjs` and the benchmark print it.
cde22eaStop the per-frame sync from inflating the app's solver numberJeremy Magland 282- **competing with the renderer.** The page draws two spheres through WebGL on the
283 same GPU, in its own animation loop. The soak has no renderer, so comparing the
284 soak against the app's `solver` separates contention from everything else.
285- **clocks.** An animation-paced loop leaves the GPU idle for most of each 16 ms
286 frame, so it may never leave its low-power state, while the benchmark hammers it
287 continuously and boosts. On a thermally managed laptop this alone can be worth a
3221abdAdd a Benchmark button that measures the solver and the GPU's clock rampJeremy Magland 288 factor of two, and it is not something the code can fix. The **Benchmark**
289 button's ramp figure measures it directly.
290- **anything else using the GPU.** Another process competing for it changes
291 whichever run overlaps it, which makes a comparison across two separate
292 invocations meaningless. `compare-perf.mjs` runs both sides back to back in one
293 invocation partly for this reason.
294- **not buffer robustness**, another plausible suspect: WebGPU clamps every array
295 access for safety, which could cost real time in the transform kernels' inner
296 loops. Measured with Dawn's `disable_robustness` toggle
297 (`DAWN_FLAGS='enable-dawn-features=disable_robustness' npm run bench`), it makes
298 no difference here at all — 0.59 ms/step either way.
cde22eaStop the per-frame sync from inflating the app's solver numberJeremy Magland 299- **which browser.** WebGPU implementations differ substantially in maturity;
300 Chrome and Safari are not interchangeable for this.
302None of these change *what* is computed — see below for how to confirm that
303independently.
0f3abbdSeparate solver time from frame time, and add a cross-environment checkJeremy Magland 305### Is it really the same computation?
307```
308node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
309```
311runs one identical spec on the desktop and in a real browser and compares the
312final spectral state. The pipeline is deterministic given (model source,
313parameters, lmax, seed, steps) — a seeded PRNG, then fixed arithmetic — so the two
314should agree to fp32 round-off. Both sides build their spec through the same
315`parseArgs`, so neither can quietly use a different default.
317They will *not* agree bit for bit; GPUs differ in fused-multiply-add and other
318latitude fp32 allows. Between Intel Xe (via Dawn) and SwiftShader — about as
319different as two implementations get — 200 steps at lmax 31 agree to a relative
320L2 of **2e-6**.
322It also reports which **Fourier stage** each side chose. `ShtPlan` picks FFT or
323DFT from the device's workgroup-storage and invocation limits, and those are
324genuinely different algorithms that round differently, so a mismatch there
325explains a difference in the values rather than being a symptom of one. The app's
326stats line and the benchmark both print the chosen stage for the same reason.
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 328## Against upstream SHTNS
330The transforms are a WGSL translation of
331[SHTNS](https://nschaeff.bitbucket.io/shtns/), and the tests check them against
332their own f64 CPU twin — which shows they are self-consistent, not how they
333compare with the library they are modeled on. SHTNS itself runs on the CPU with
334hand-tuned SIMD codelets, and on Nvidia GPUs with its own CUDA kernels, including
335a single-precision mode. That is a direct comparison, and
336[`bench/shtns/`](bench/shtns/) makes it:
338```
339cd bench/shtns && ./bootstrap.sh && make # clone SHTns at a pinned commit, build
340node scripts/compare-native.mjs --check # then, from the repo root
341```
343`bootstrap.sh` adds CUDA support when `nvcc` is on `PATH`, so the same tree gives
344you the CPU comparison anywhere and the GPU one on a machine with an Nvidia card.
345`compare-native.mjs` runs every implementation present, back to back in one
346invocation so a second process competing for the GPU affects both sides rather
4166b48Tidy up after the SHTNS comparisonJeremy Magland 347than one, and prints them in one table. On an RTX PRO 6000 Blackwell, at the
348app's default lmax:
350```
351 grid lmax 63 · 128×256 · nlm 2,080 (one synthesis + one analysis per round trip)
4166b48Tidy up after the SHTNS comparisonJeremy Magland 353 webgpu 0.084 ms/round trip 11848/s (baseline) fp32
354 NVIDIA (blackwell), via Dawn · CPU-side launching 0.012 ms/step
355 shtns cuda 0.021 ms/round trip 48009/s 0.25x webgpu fp32
356 NVIDIA RTX PRO 6000 (sm_120, 188 SMs) · CPU-side launching 0.018 ms/step
357 shtns cpu 0.072 ms/round trip 13982/s 0.85x webgpu fp64
359```
4166b48Tidy up after the SHTNS comparisonJeremy Magland 361Read that carefully rather than as "4x". The two GPU rows are limited by different
362things: the WGSL row spends 14% of its time on the CPU and is genuinely GPU-bound,
363while SHTNS spends **86%** — 0.018 ms of 0.021 — queueing its six-or-so kernels, so
364its number is close to what it costs to *submit* a round trip on that host and its
365actual GPU time is below that and unresolved. The 4x is a lower bound on the gap in
366GPU work, not a measurement of it. `compare-native.mjs` flags any row above 50%
367for this reason.
369The other number worth noticing is the third row: one CPU core in fp64 is about
370level with the WGSL transforms on a 188-SM datacentre GPU. At lmax 63 there are
3712,080 coefficients on a 128×256 grid — far too little work to occupy that card, so
372this says more about occupancy than about the shaders. Sweep lmax before drawing
373conclusions, and stop at 511: above that `16*nphi` exceeds the workgroup-storage
374limit, the FFT stage falls back to the O(nphi·mmax) DFT, and the comparison stops
375being about the FFT.
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 377Two things are measured, because they answer different questions:
379- **transforms** (`npm run bench:sht` here, `--mode transform` there) — one
380 spectral → grid → spectral round trip and nothing else. This is the
381 library-against-library number, and since the transforms are ~96% of the
382 solver's compute it is what decides how fast the solver can be.
383- **solver** (`npm run bench` here, `--mode solver` there) — a whole IMEX Euler
384 timestep, which is what the app's `solver` line reports.
386`--check` diffs the final spectral state across implementations, which is what
387makes the timing mean anything: two numbers are only comparable if they are the
388cost of the same computation. That check is possible at all because the spectral
389layout and normalization are SHTNS's own — orthonormal with Condon–Shortley,
390coefficients grouped by `m`, `LM(l,m)` agreeing index for index — so a state can
391be diffed element by element with no reindexing. Over 20 steps, fp32 WGSL against
392fp64 SHTNS agrees to **~1e-6** relative L2, for every model.
394It is also the check on the one second implementation this repo has. The native
395solver cannot run `models/<key>.m` — C has no numbl — so `bench/shtns/spec.h`
396restates the same arithmetic, one line per line of MATLAB. `--check` is what
397keeps that transcription honest, and `compare-native.mjs` refuses to compare two
398runs whose resolved grid or parameters disagree, which is the other way the two
399sides could drift.
401[`bench/shtns/README.md`](bench/shtns/README.md) lists what is *not* identical and
402should be kept in mind when reading the ratio — SHTNS runs its Legendre
403recurrence in fp64 even in fp32 mode for `lmax <= 128` (WebGPU has no fp64 at
404all), the Fourier stages are cuFFT/VkFFT/FFTW against a WGSL FFT, and SHTNS'
405polar optimization is off by default here because we have none.
4166b48Tidy up after the SHTNS comparisonJeremy Magland 407How much the grid size matters is easiest to see on a weak GPU, where there is no
408launch-overhead floor to hide behind. On an Intel Xe iGPU against one core of the
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 409same laptop, one round trip costs:
411| lmax | grid | WGSL (fp32) | SHTNS, 1 CPU core (fp64) |
412|---|---|---|---|
413| 31 | 64×128 | 0.183 ms | 0.017 ms |
414| 63 | 128×256 | 0.250 ms | 0.110 ms |
415| 127 | 256×512 | 0.733 ms | 0.602 ms |
41710x behind at lmax 31, 1.2x at lmax 127 — the same comparison, on the same two
418chips. Whatever a single number says, it is saying it about one grid size.
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 420## Tests
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 422There is no second implementation of the solver to diff against, so the `.m` path
423is checked against **closed-form answers**. Each case is one whose evolution is
424known exactly, run through the whole real pipeline — MATLAB source, numbl
425lowering, generated WGSL, GPU transforms — and compared with arithmetic
426([`test/analyticChecks.ts`](test/analyticChecks.ts)):
428- **A** — a linear reaction `f(u) = c*u` leaves every spherical-harmonic mode
429 independent, growing by exactly `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. This
430 pins the transform round-trip, the eigenvalue mapping, the IMEX update and the
431 state feedback at once, and checks that nothing leaks between modes. Agrees to
432 ~2e-7 over 20 steps.
433- **B** — a nonlinear reaction on a *uniform* field stays uniform and diffusion
434 cannot touch it, so each step is exactly the scalar ODE map. Agrees to 1.5e-8
435 over 25 steps. Checks that a generated kernel evaluates a nonlinear reaction.
436- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
437 linearized 2x2 IMEX recurrence, and the `(l=24, m=7)` mode is confirmed
438 unstable. Looser (~2e-3) because fp32 keeps only about four digits of a
439 perturbation that small.
441Two test models exist only for this: [`test/models/linear.m`](test/models/linear.m)
442and [`test/models/logistic.m`](test/models/logistic.m).
444Alongside those, [`test/modelChecks.ts`](test/modelChecks.ts) compiles every model
445the app offers and asserts **how many kernels it compiles to**. That is a fusion
446guard: numbl's lowering emits one statement per *operator* and its inline pass
447folds them back into per-line expression trees, and if that stops happening the
448results stay correct while every operator becomes its own dispatch. It is
449invisible in the numbers, so it is asserted directly. (It has already caught one
450regression.)
452[`test/transformChecks.ts`](test/transformChecks.ts) is the one remaining
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 453implementation-vs-implementation check inside the suite, comparing the WGSL
454transforms against shtns-webgpu's f64 CPU twin. Comparing them against *upstream*
455SHTNS is a separate, opt-in step, because it needs a native toolchain — see
456[Against upstream SHTNS](#against-upstream-shtns).
458All three modules run in **both** environments, so the two GPU stacks get the same
459guarantees:
461- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
462 pass `--skip-without-gpu` to let a machine without one say so and move on
463 (which is what CI does, since the browser suite covers the same modules).
464- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
465 Also runs the soak.
467Other commands:
469- `npm run bench -- --help` — the desktop benchmark (see
15a77e2Add a desktop WebGPU benchmark and show its command in the appJeremy Magland 470 [Desktop vs browser](#desktop-vs-browser)).
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 471- `npm run bench:sht -- --help` — the transforms alone, with no solver around
472 them, for comparing against upstream SHTNS.
30ed90fAdd scripts/diagnose-sht.ts: which stage of the transform is wrong?Jeremy Magland 473- `npx vite-node scripts/diagnose-sht.ts` — when the transform tests fail on a GPU,
474 say *which* stage is wrong. It reads the intermediate `fm` back out and scores
475 the Legendre and Fourier stages of each direction separately against the f64
476 reference, then breaks the error down by order `m` and by latitude.
4166b48Tidy up after the SHTNS comparisonJeremy Magland 477- `npx vite-node scripts/diagnose-leg.ts [--m 0]` — the follow-up to that: read the
478 Legendre recurrence out of the production shader term by term, by synthesizing a
479 spectrum that is 1 at a single coefficient, and compare each `ỹ_l^m` with the f64
480 reference. The first term that disagrees names the culprit.
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 481- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
482 pattern saturates into O(1)-contrast spots rather than decaying or diverging.
483- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
484 sampling JS heap and catching crashes.
485- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
486 demo after a number of steps.
e7bcd70Add soak, live-check and solver-only soak toolingJeremy Magland 487- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL in a real
488 browser: load, press Run, confirm the solver advances.
0f3abbdSeparate solver time from frame time, and add a cross-environment checkJeremy Magland 489- `node scripts/compare-env.mjs` — run one identical spec on the desktop and in a
490 browser and compare the final state (see
491 [Is it really the same computation?](#is-it-really-the-same-computation)).
17db8f1Add scripts/compare-perf.mjs, and rule out CPU command encodingJeremy Magland 492- `node scripts/compare-perf.mjs` — measure the same solver work in both and split
493 the difference (see
494 [Why the browser is slower](#why-the-browser-is-slower-and-how-to-find-out-by-how-much)).
b689087Benchmark the WGSL transforms against upstream SHTNSJeremy Magland 495- `node scripts/compare-native.mjs` — run one spec through the WGSL transforms and
496 through upstream SHTNS, and line the numbers up (see
497 [Against upstream SHTNS](#against-upstream-shtns)). Needs
498 [`bench/shtns/`](bench/shtns/) built first.
e7bcd70Add soak, live-check and solver-only soak toolingJeremy Magland 499- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
501### A note on canvas resizing
503Early long runs killed the browser after ~700–800 steps. The cause was the
504colorbar's min/max labels changing width as their digit count changed, which
505reflowed the panel, fired the `ResizeObserver`, and called
506`renderer.setSize()` — reallocating the WebGL drawing buffer. Assigning
507`canvas.width` also blanks the canvas even when the value is unchanged, so the
508same bug caused visible flicker. Fixed by giving the colorbar column a fixed
509width and making `SphereScene.resize()` return early on no-op resizes.
4166b48Tidy up after the SHTNS comparisonJeremy Magland 511### A note on the Legendre recurrence on Blackwell
513The first run on an Nvidia GPU — an RTX PRO 6000, driver 590.48, reached through
514Dawn's Vulkan backend — failed 11 of the tests. `synth` was off by 5.5e+3 while
515`analys` was accurate to 7.3e-7, and the solver produced NaN within 40 steps.
517The two diagnostic scripts above were written for it and localized it in two
518steps: `leg_synth` was the only wrong shader, and within it the recurrence was
519right at `l = m` and `l = m+1` and then returned *exactly zero* at `l = m+2`, at
520every latitude. That is not a precision failure. It is
522```wgsl
523let c0 = ab[base + (l + 2u - m)];
524y0 = c0.x * ct * y1 + c0.y * y0; // c0 reads as (0, 0) on the first iteration
525```
527with the `ab` read two lines later working fine. The buffer was not at fault:
528`leg_analys` reads the same array correctly on the same device, and `m = 62, 63`
529— the only orders whose loop breaks before that line — were the only correct
530ones. Nothing about that WGSL is invalid, so it was a miscompiled load.
532Fixed by giving the advance the shape `leg_analys` already used, which that
533driver compiles correctly: both coefficients fetched unconditionally, and the new
534`y0` carried in a temporary rather than assigned and then read back by the `y1`
535update. Two shaders doing the same recurrence should have agreed on form anyway.
537Worth knowing for what it says about the transforms in general: nothing had
538exercised them on Nvidia hardware before, and the existing test caught it
539immediately — it just could not say where. That is what the diagnostics are for.
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 541## Development
543```
544npm install
545npm run dev # local dev server
546npm run build # type-check + production build to dist/
547```
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 549### The numbl dependency
551numbl is a local `file:../../numbl` dependency, so a sibling checkout of
552[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
553compiler internals — parser, lowerer, IR, inline pass — which its package
554`exports` map does not publish, so they are reached through the `numbl-src` path
555alias in [`vite.config.ts`](vite.config.ts).
557The exact surface we depend on is written down in
558[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
559*that*, not against numbl's sources. This keeps this project's compiler settings
560independent of numbl's (its sources do not type-check under the stricter options
561used here), and means a change to one of those shapes upstream breaks the build
562here with a clear diff rather than deep inside numbl's tree.
564The compiler is ~395 kB gzipped and lands in its own chunk. That is the cost of
565compiling MATLAB in the page; a build-time lowering step could remove it at the
566price of no longer being editable live.
568CI clones numbl to the sibling path that `file:` dependency expects, pinned to a
569commit. Two details make that work, both verified by building against a checkout
570that had none of numbl's own dependencies installed:
572- **numbl's `node_modules` are not needed.** The slice we import — parser,
573 lowering, IR, inline pass — is self-contained TypeScript. (Other parts of numbl
574 do import `three`, `react` and `fflate`; we never reach them.)
575- **the install must pass `--ignore-scripts`.** npm runs a linked package's
576 `prepare` script, and numbl's is `husky`, which is not installed in CI.
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 578The `scripts/*.ts` entry points that touch the compiler (the benchmark, the node
579tests, the long run) go through `vite-node`, so they resolve imports exactly as the
580browser build does — the `numbl-src` alias and the `?raw` model imports included.
581Plain `node` cannot: numbl's sources import each other as `./foo.js` while the
582files are `.ts`, which needs a bundler's resolution. Scripts that do not touch the
583compiler (`soak.mjs`, `screenshot.mjs`, `check-live.mjs`, `test-gpu.mjs`) are plain
584`.mjs` and run under `node` directly.
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 586Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
588## License
590CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).