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*.
13> **The geometry is in the operator.** The models solve with the surface's
14> Laplace–Beltrami operator, iterated by a fixed-count preconditioned
15> Richardson solve ([`solvers/richardson.m`](solvers/richardson.m), applying
16> [`lib/dlap.m`](lib/dlap.m)). The iteration count is fixed at compile time
17> with no residual check, so a shape/timestep/diffusivity combination outside
18> its convergence radius diverges over many steps rather than being caught —
19> the tests pin the known cases. See
20> [Where the geometry enters the operator](#where-the-geometry-enters-the-operator).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 21
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
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 106what keeps the cost to a few transforms per step rather than a full elliptic
107solve.
109The pieces of that sentence are separate files, because they are separate
110ideas. The **operator** — `dlap` applied to a spectral field — is
111[`lib/dlap.m`](lib/dlap.m). The **solver** — the fixed point above, iterated
112`niter` times — is [`solvers/richardson.m`](solvers/richardson.m):
114```matlab
115function X = richardson(B, dtD, lam, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter)
116 X = B ./ (1 + dtD * lam);
117 for k = 1:niter
118 dL = dlap(X, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
119 X = (B + dtD * dL) ./ (1 + dtD * lam);
120 end
121end
122```
124And a **model** is a reaction plus one solve per species — the whole of
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 125[`models/schnakenberg.m`](models/schnakenberg.m)'s step is:
127```matlab
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 128function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, niter)
130 v = synth(V);
131 uuv = u .* u .* v;
133 Bu = U + dt * analys(a - u + uuv);
134 Bv = V + dt * analys(b - uuv);
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 136 Un = richardson(Bu, dt * D1, lam, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter);
137 Vn = richardson(Bv, dt * D2, lam, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter);
139```
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 141Trying a different solver against the same operator is a change to those two
142call lines: every solver composes from `dlap` (the matvec is
143`(1 + dtD.*lam).*x - dtD.*dlap(x)`, the preconditioner the elementwise
144divide), and which one a model calls is part of what compiles — swapping
145recompiles, like changing `niter` already does. The solver is written as a
146full re-evaluation rather than an accumulated correction on purpose: where
147`dlap` computes to zero there is no correction to mis-round, and the divide is
148turing-sphere's arithmetic unchanged.
150Three solvers ship. [`solvers/bicgstab.m`](solvers/bicgstab.m) solves the
151same system by preconditioned BiCGSTAB — same `dlap`, same preconditioner, a
152Krylov recurrence instead of a stationary one, at two `dlap` evaluations per
153iteration instead of one. Its scalars (`rho`, `alpha`, `omega`) never touch
154the CPU: `dot` is a GPU reduction into a 1-element buffer, the recurrences on
155its results compile to 1-element kernels, and a single-element value
156broadcasts into the vector updates. Inner products carry the half-spectrum
157weight `wlm` (m > 0 counts twice), making them the real L2 inner products on
158the sphere. With no residual test, every ratio `a/b` is written in the
159guarded form `a*b/(b*b + 1e-30)`, so a converged (or broken-down) iteration
160goes stationary instead of dividing noise by noise. The difference is not
161academic: at the app's default lmax, Schnakenberg on the peanut sits outside
162the Richardson iteration's convergence radius for `niter ≥ 2` and diverges,
163while BiCGSTAB on the identical operator converges monotonically — the tests
164pin both behaviors, side by side.
166[`solvers/gmres.m`](solvers/gmres.m) is right-preconditioned GMRES(niter) —
167one Arnoldi sweep, no restart — with the residual minimized over the whole
168Krylov space. Its bookkeeping is what the other solvers never need: a basis
169of niter+1 spectral fields, a Hessenberg matrix, Givens rotations, a
170triangular back-substitution. The basis lives in a *bank* (`getslab` /
171`setslab`: the k-th 2 × nlm field of a wider array), the small matrices are
172element-addressed (`getat` / `setat`), and both are functional updates the
173planner compiles to static-offset buffer copies — MATLAB's own `H(i,j) = h`
174cannot lower, because numbl must prove an indexed write in bounds before the
175loop unrolls, and a loop variable has no value yet at that point. Written as
176calls, the index resolves at *planning*, where unrolling has made it a
177literal. The same resolution lets an inner loop bound depend on the outer
178loop's variable, which is what makes the `for i = 1:j` orthogonalization
179sweep compile.
181### Where the geometry enters the operator
183[`lib/dlap.m`](lib/dlap.m) is Algorithm 3 of the evolving-surface notes: the
184field's θ/φ derivatives (the `dtheta`/`dphi` transforms,
185[`src/sht/deriv.ts`](src/sht/deriv.ts)) are contracted through the inverse
186metric quantities into a tangential gradient; each Cartesian component is
187re-analysed and differentiated again; the results recombine into the surface
188divergence, and `lam .* F` adds back what the round-sphere part already
189carries. The metric quantities `Vt*`/`Vp*`
190([`src/geom/metric.ts`](src/geom/metric.ts)) are built once from the
191embedding's derivatives when the geometry is (re)built — the geometry is
192static, so per step they are just six more buffers the kernels read. `filt`
193zeroes the top two spectral degrees wherever the operator re-differentiates,
194because the derivative recurrences cannot exactly represent a derivative
195there.
197On the sphere `dlap` computes to (numerical) zero, so any `niter` lands
198within transform round-off of the exact round-sphere answer — asserted in the
199tests. Off the sphere the correction genuinely moves the answer, and
200convergence is a real constraint: the fixed-count loop has no residual check,
201so the tests also pin which shape/niter combinations are known to sit outside
202the convergence radius and diverge.
204### Subroutines
206A model file is not limited to `init` and `step`: it can define further
207functions and call them, and every model compiles against the shared library
208files — [`lib/`](lib/) for operators, [`solvers/`](solvers/) for solvers —
209with MATLAB's visibility rules (a file's namesake function is public; a
210model-local function of the same name shadows it). numbl specializes each
211callee for the argument types at its call sites, and the host then splices
212the lowered body into the caller, one clone per call site
213([`src/mgpu/inlineCalls.ts`](src/mgpu/inlineCalls.ts)): arguments bind by
214renaming rather than copying, and assignments to a callee output become
215assignments to the caller's variable, which is what lets a solver iterate its
216result in place. Expansion runs before the fusion pass, so a call fuses
217exactly as the same code written inline would — the boundary costs nothing,
218and `describe()`'s op listing names the expanded internals
219(`richardson#1.X`). Recursion cannot unroll into a fixed op sequence and is
220refused at compile time, like a runtime loop bound.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 221
222### `for` loops, unrolled
224A plan is a fixed list of GPU operations with no branching, which is what makes
225a timestep pure command recording — one submit, no CPU in the loop. A counted
226loop still fits: the planner
227([`src/mgpu/plan.ts`](src/mgpu/plan.ts)) unrolls it, planning the body once per
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 228iteration. The loop that matters is `solvers/richardson.m`'s `for k = 1:niter`,
229expanded into each model's step at every solve call site.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 230
231Nothing else had to change for that, because numbl gives a variable one cName
232for every assignment to it: the buffer an iteration writes is the buffer the
233next one reads, which is exactly a loop-carried value. The loop variable gets no
234buffer at all — it is bound as a derived scalar to that iteration's literal, so
235a kernel reading `k` folds the number in.
237Two consequences worth stating:
239- **The bounds must be known when the model compiles.** `niter` is supplied as a
240 fixed scalar rather than a tunable one, so changing it recompiles — unlike a
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 241 parameter, which is a uniform. A bound may also be an enclosing unrolled
242 loop's variable (`for i = 1:j` — each unrolled `j` plans its own inner trip
243 count, which is how GMRES's triangular sweeps compile). A genuinely runtime
244 bound is refused at compile time with a source position, not silently
245 mis-compiled, and there is a test for that.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 246- **Fusion survives.** numbl's inline pass recurses into loop bodies, so a line
247 inside the loop is still one kernel. It runs there with no protected names,
248 though, which means an assignment whose only visible use is later in the same
249 body can be elided — correct for a body-local temp, wrong if something outside
250 the loop wanted it. [`src/mgpu/compile.ts`](src/mgpu/compile.ts) snapshots what
251 each loop body assigns before the pass and refuses the ones that escape, so
252 that case is a compile error rather than a stale read.
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 254Unrolling is exactly linear in the trip count: 26 GPU ops per species per
255iteration (the operator's twelve transforms and the solve's kernels), asserted
256in the tests.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 257
258## MATLAB, compiled to WebGPU
260Unchanged from turing-sphere, and it now compiles the geometry files too. numbl
261parses and lowers each function for the concrete argument types of the current
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 262grid; user-function calls are expanded into the caller, one clone per call
263site; the inline pass folds single-use temps back into their consumer, so one
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 264line of MATLAB becomes one expression tree; and this repo emits one WGSL compute
265kernel per element-wise statement
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 266([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). `synth` / `analys` (and the
267derivative pair `dtheta` / `dphi`) are external operations whose type rules
268numbl learns from a `.mtoc2.js` workspace file, and which the backend maps onto
269the spherical-harmonic pipelines; `dot` is one more, mapped onto a
270single-dispatch reduction ([`src/mgpu/reduce.ts`](src/mgpu/reduce.ts)) whose
2711-element result stays on the GPU — scalars computed from it become 1-element
272kernels, and reading one inside a vector expression broadcasts it. The
273indexed-access ops (`getslab`/`setslab`, `getat`/`setat`) compile to
274static-offset buffer copies, their indices evaluated at planning time where
275the unrolled loop's variable is a literal. Anything it cannot express is
276refused at compile time with a source position.
278The Schnakenberg step above compiles to 65 GPU operations at one solve
279iteration: 28 transforms, 35 generated kernels, and 2 buffer copies feeding the
282Two consequences carried over:
284- **The step is synchronous.** WebGPU's encode path is synchronous and every
285 pipeline is built once at compile time, so a timestep is pure command
286 recording; the only `await` in the loop is the single readback per rendered
287 frame.
288- **Parameters are uniforms, not constants.** Moving a slider rewrites a small
289 buffer instead of triggering a recompile. Editing the MATLAB recompiles;
290 changing `dt` does not. `niter` is the deliberate exception, above.
292## Provenance
294- **turing-sphere**, which this is a fork of: the solver, the transforms
295 backend, the compilation path, the benchmarks and the analytic tests.
296- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
297 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
298 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Vendored under
299 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference transform
300 used for testing.
301- **Rendering:** three.js meshes with per-vertex colormaps, adapted from the
302 `SphereEmbedding` view in
303 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
304 extension package ([`src/render/`](src/render/)). That view displays a
305 time-varying embedded geometry with fields on it, which is the same picture
306 this draws — including its sphere/surface morph, which turing-sphere had
307 dropped as having nothing to morph to.
309turing-sphere additionally carries a comparison against a native build of
310upstream SHTNS ([`bench/shtns/`](https://github.com/concept-collection/turing-sphere/tree/main/bench/shtns)).
311That is not duplicated here: the transforms are the same code, and its C-side
312transcription of the model would have to be maintained against a step this
313project intends to change.
315Because the algorithm is compiled to compute shaders, **WebGPU is required** —
316there is no CPU fallback (the f64 CPU transform remains, for tests).
318## Numerics
320- Grid: Gauss–Legendre × equispaced-φ, dealiased for the cubic reactions with
321 the `(pdeg+1)` rule: `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1`
322 (rounded up to a power of two for the GPU FFT path). At the default lmax 63
323 that is a 128×256 grid.
324- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
325 coefficients for m ≥ 0, m-major ordering.
326- fp32 transforms introduce ~1e-6 relative error per step; for pattern formation
327 from 1e-2 seeded noise this is inconsequential. The geometry goes through one
328 analysis/synthesis round trip and picks up the same round-off: the unit sphere
329 comes back with radius 1 to ~2e-5 under Dawn, ~4e-4 under SwiftShader.
330- The shipped geometries are all degree ≤ 5, far below any lmax the app offers,
331 so band-limiting removes nothing from them. A shape you write yourself may not
332 be so lucky — see the note in [`geometries/bumpy.m`](geometries/bumpy.m).
334## Desktop vs browser
336[`scripts/bench.ts`](scripts/bench.ts) runs the same thing the app runs — same
337`.m`, same generated WGSL, same transforms — from Node on desktop WebGPU (Google
338Dawn), and the app prints the command line that reproduces whatever it is
339currently simulating:
341```
342npm run bench -- --preset schnak-spots --geometry ellipsoid --lmax 63 --niter 1 \
343 --steps 2000 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05 \
344 --gax 1.5 --gay 1 --gaz 0.6
345```
347Copy it from under the stats line and compare the `ms/step` it reports with the
348app's. Both sides go through the one shared
349[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
350command, the benchmark parses it back — so there is no second copy of the
351defaults for the two runs to drift apart on. Geometry parameters take a `g`
352prefix (`--gwaist`) so a shape parameter can never collide with a model one.
354The app reports **two** numbers and only the first is comparable to the
355benchmark: `solver` is the batch of steps alone, waited for but not read back;
356`ms/frame` additionally carries a GPU→CPU readback per species, the
357colormapping, and the vertex upload. Those per-frame costs are fixed and do not
358shrink when the GPU gets faster, so on a quick GPU a frame can easily cost ten
359times the steps inside it. That is expected and is not the solver being slower
360in the browser.
362To attribute the gap rather than guess at it:
364```
365node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
366```
368measures the same solver work in both — batched, nothing read back, no rendering
369on either side — and reports each with its CPU-encoding share, the Fourier
370stage, and the adapter. It stops you first if the two are not even the same
371device, which is a common cause of "the browser is much slower". Both sides
372resolve the geometry and the iteration count from the same constants, because
373the iteration count is unrolled into the step and a mismatch would compare two
374different amounts of work.
376The app's **Benchmark** button runs the same measurement in the page, plus the
377**ramp** — the first third of the run against the last. GPUs downclock when
378idle and an animation-paced loop leaves them idle most of every frame, so a
379large ramp means the steady-state number is limited by clocks rather than work.
381### Is it really the same computation?
383```
384node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
385```
387runs one identical spec on the desktop and in a real browser and compares the
388final spectral state. The pipeline is deterministic given (model source,
389geometry, parameters, lmax, niter, seed, steps), so the two should agree to fp32
390round-off — not bit for bit, since GPUs differ in fused-multiply-add and other
391latitude fp32 allows. It also reports which Fourier stage each side chose, since
392FFT and DFT are genuinely different algorithms that round differently.
394Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), an
395optional dependency so that an unsupported platform fails the install of that
396package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
397`--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
398(`;`-separated) passes Dawn options through.
400## Tests
402There is no second implementation of the solver to diff against, so the `.m`
403path is checked against **closed-form answers** and against **exact structural
404properties**. Four modules, run in both environments:
406[`test/analyticChecks.ts`](test/analyticChecks.ts) — cases whose evolution is
407known exactly, run through the whole real pipeline. All three are statements
408about the round sphere, so all three build on the sphere geometry:
410- **A** — a linear reaction leaves every mode independent, growing by exactly
411 `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. Pins the transform round trip, the
412 eigenvalue mapping, the IMEX update and the state feedback at once. ~2e-7 over
413 20 steps.
414- **B** — a nonlinear reaction on a uniform field stays uniform, so each step is
415 exactly the scalar ODE map. 1.5e-8 over 25 steps.
416- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
417 linearized 2×2 IMEX recurrence, and `(l=24, m=7)` is confirmed unstable.
418 Looser (~4e-3) because fp32 keeps about four digits of a perturbation that
419 small.
421[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface and the loop:
423- every geometry compiles and closes; the sphere has radius 1 everywhere and is
424 **exactly degree 1** in the harmonics, which is what makes the reference case
425 exact rather than merely accurate;
426- the peanut matches its own closed-form radial profile at every grid point, and
427 **the same coefficients give the same surface on a 2× grid** — the 2× Gauss
428 latitudes share no point with the 1× ones, so agreeing there is agreeing
429 everywhere, which is what "rendered exactly, not subdivided" means;
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 430- unrolling is **exactly linear** in the trip count; on the sphere the
431 geometric correction computes to (numerical) zero, so 0, 1 and 4 iterations
432 agree to transform round-off; on the peanut it **measurably moves the
433 answer**;
434- a niter × geometry sweep stays finite except the combinations **known to sit
435 outside the Richardson convergence radius**, which are pinned as diverging —
436 and on exactly those combinations **bicgstab and gmres keep converging**,
437 also pinned;
438- at equal niter, **bicgstab and gmres land far closer to the converged
439 answer** than richardson on the same operator;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 440- a runtime loop bound is refused at compile time;
441- swapping the surface mid-run leaves the spectral state untouched.
443[`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
444and asserts **how many kernels it compiles to**, split into the base step and
445what one solve iteration adds. That is a fusion guard: if numbl's inline pass
446stops folding, the results stay correct while every operator becomes its own
59f3e22Factor the solver out of the models; add BiCGSTAB and GMRESJeremy Magland 447dispatch, which is invisible in the numbers. It also compiles a model built of
448**user-defined subroutines** — multi-output, scalar-returning, a solver-like
449local with its own loop — asserts recursion is refused, and checks the
450**reduction and indexing primitives** directly against exact expected values:
451the `dot` sum and the `wlm`-weighted inner product against the CPU, scalar
452arithmetic on a GPU-resident result bit for bit, the 1-element broadcast, and
453element/slab round-trips through `setat`/`getat` and `setslab`/`getslab`.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 454
455[`test/transformChecks.ts`](test/transformChecks.ts) compares the WGSL transforms
456against shtns-webgpu's f64 CPU twin.
458- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
459 `--skip-without-gpu` lets a machine without one say so and move on (which is
460 what CI does, since the browser suite covers the same modules).
461- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
462 Also runs the soak. A few geometry tolerances are set by SwiftShader's fp32,
463 which is about an order of magnitude looser than Dawn's.
465Other commands:
467- `npm run bench -- --help` — the desktop benchmark.
468- `npm run bench:sht -- --help` — the transforms alone, no solver.
469- `npx vite-node scripts/diagnose-sht.ts` — when the transform tests fail on a
470 GPU, say *which* stage is wrong.
471- `npx vite-node scripts/diagnose-leg.ts [--m 0]` — read the Legendre recurrence
472 out of the production shader term by term.
473- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
474 pattern saturates rather than decaying or diverging.
475- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
476 sampling JS heap and catching crashes.
477- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
478 demo after a number of steps.
479- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL.
480- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
482## Development
484```
485npm install
486npm run dev # local dev server
487npm run build # type-check + production build to dist/
488```
490### The numbl dependency
492numbl is a local `file:../../numbl` dependency, so a sibling checkout of
493[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
494compiler internals — parser, lowerer, IR, inline pass — which its package
495`exports` map does not publish, so they are reached through the `numbl-src` path
496alias in [`vite.config.ts`](vite.config.ts).
498The exact surface we depend on is written down in
499[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
500*that*, not against numbl's sources. This keeps this project's compiler settings
501independent of numbl's, and means a change to one of those shapes upstream
502breaks the build here with a clear diff rather than deep inside numbl's tree.
503The `For` IR node is spelled out there, since the planner now walks it.
505CI clones numbl to the sibling path that the `file:` dependency expects, pinned
506to a commit, with `--ignore-scripts` (npm runs a linked package's `prepare`
507script, and numbl's is husky). numbl's own `node_modules` are not needed: the
508slice we import is self-contained TypeScript.
510The `scripts/*.ts` entry points that touch the compiler go through `vite-node`,
511so they resolve imports exactly as the browser build does. Plain `node` cannot:
512numbl's sources import each other as `./foo.js` while the files are `.ts`.
514Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
516## License
518CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).