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*.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 12The geometry is in the operator: the models evaluate the surface
13Laplace–Beltrami operator `lap_g` inside the implicit solve, in a **flux form
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 14that costs 5 spherical-harmonic transforms per species per iteration** (plus
15one Legendre-free FFT derivative) where the textbook Cartesian-gradient form
16needs 12. See
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 17[The geometry in the operator](#the-geometry-in-the-operator) and
18[docs/reduced-transforms.md](docs/reduced-transforms.md).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 19
20## What a surface is here
22A geometry is an embedding of the sphere into R³: three scalar fields x, y, z
23over the (θ, φ) parametrization, each carried as spherical-harmonic
24coefficients. The unit sphere is the case where all three are pure degree-1
25harmonics.
27You write one down as MATLAB, in [`geometries/`](geometries/):
29```matlab
30function [gx, gy, gz] = shape(theta, phi, waist, stretch)
31 st = sin(theta);
32 r = 1 - waist * (st .^ 2);
33 gx = r .* (st .* cos(phi));
34 gy = r .* (st .* sin(phi));
35 gz = (1 + stretch) * (r .* cos(theta));
36end
37```
39That is ordinary element-wise MATLAB and goes through the same compiler and the
40same WGSL backend the models do. It is evaluated once on the solver's grid, and
41then **analysed into coefficients**, which is the form everything downstream
42uses. Two things follow from going through the coefficients rather than keeping
43the pointwise values:
45- **It is exactly band-limited at lmax.** The surface has as many derivatives as
46 the scheme needs and no aliased content the solver cannot see. What the solver
47 and the renderer both use is the *synthesis* of the coefficients, so for a
48 shape with sharp features the surface being solved on is not quite the one
49 that was written down — which is the honest thing for a spectral method to do.
50- **It can be evaluated on any grid.** The renderer draws the surface on the
51 (possibly finer) display grid by synthesizing the same coefficients there.
52 That is exact interpolation, not subdivision — the same argument that lets the
53 species fields be oversampled, and it is checked directly in the tests.
55Four geometries ship: [sphere](geometries/sphere.m) (the reference case),
56[ellipsoid](geometries/ellipsoid.m), [peanut](geometries/peanut.m) — a dumbbell
57whose waist is a saddle — and [bumpy](geometries/bumpy.m). Each is editable in
58the page, with its own parameters. Changing a shape does not recompile the
59solver and does not disturb the run: the geometry is data whose shape in the
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 60bindings depends only on the grid, so a swap is sixteen buffer writes and the
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 61pattern carries straight on.
63A **morph** slider blends the drawn surface back to the unit sphere. The
64parametrization is the sphere's either way, so sweeping it shows which point
65went where.
67## The scheme, and where the geometry enters
69It solves the N-species system
71```
72d(u_k)/dt = D_k*lap_g(u_k) + f_k(t, u_1, ..., u_N), k = 1, ..., N
73```
75where `lap_g` is the Laplace–Beltrami operator of the surface. On the round
76sphere `lap_g` is diagonal in spherical-harmonic space with eigenvalues
77`-l(l+1)`, which is what makes turing-sphere's implicit diffusion a single
78divide. On a general surface it is not diagonal, and not even constant-
79coefficient, so that divide has to become a solve.
81The models split the operator:
83```
84lap_g = lap_s + dlap
85```
87with `lap_s` the round-sphere one. `(I - dt*D*lap_s)` is still exactly
88invertible, so the implicit step
90```
91(I - dt*D*lap_g) Unew = B
92```
94rearranges into a fixed point that keeps the whole geometry on the right-hand
95side,
97```
98Unew = (B + dt*D*dlap(Unew)) ./ (1 + dt*D*lam)
99```
101and the loop iterates it from the round-sphere answer. That is preconditioned
102Richardson, with the operator we can invert exactly as the preconditioner; it
103converges while `dt*D*dlap` stays small against `(I - dt*D*lap_s)`, which is
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 104what keeps the cost to a few transforms per step rather than a full elliptic
105solve (see [docs/richardson-iteration.md](docs/richardson-iteration.md)). One
106species of [`models/schnakenberg.m`](models/schnakenberg.m)'s solve loop:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 107
108```matlab
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 109lamJ = lam ./ jhat; % mean-J preconditioner eigenvalues (below)
110...
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 111for k = 1:niter
112 Fu = Un .* filt; % zero the top 2 degrees before differentiating
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 113 vtu = dthetac(Fu);
114 vpu = dphic(Fu);
115 [Ftu, Fpu] = synth(vtu, vpu); % sin(theta)*dtheta(u), dphi(u) -- smooth on
116 % the sphere, one batched dispatch
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 117 Pu = p1 .* Ftu + p2 .* Fpu; % the two fluxes, also smooth: the precomputed
118 Qu = p2 .* Ftu + q2 .* Fpu; % weights carry every 1/sin(theta) there is
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 120 Pcu = PAu .* filt;
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 121 scu = dthetac(Pcu); % theta part of the divergence, coefficients
122 Lu = synth(scu); % sin(theta) * dtheta(P) on the grid
123 dQu = dphig(Qu); % d/dphi is diagonal in the Fourier index:
124 % two FFT stages, no Legendre work at all
125 lapu = r .* (Lu + dQu); % = lap_g(u) on the grid
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 126 dLu = (analys(lapu) + lamJ .* Un) .* filt; % dlap, projected onto the band
127 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
129```
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 131On the sphere `dlap` is mathematically zero — `p1 = q2 = 1`, `p2 = 0`,
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 132`r = 1/sin²θ`, `jhat = 1`, and the composition collapses to `lap_s` — so the
133sphere case reproduces turing-sphere to fp32 round-off, and the tests assert
134the state stays put across 0, 1 and 4 iterations.
136**The preconditioner folds in the symbol of the operator.** `jhat` is the
137host's minimax scale `2/(μmin + μmax)` over the eigenvalues `μ(x)` of the
138operator's principal symbol — the inverse squared principal stretches of
139the embedding, direction included, read straight off the flux-metric
140arrays (`S = (1/J)·[[p1,p2],[p2,q2]]`). Preconditioning with `lam/jhat`
141then contracts every mode *and every direction* at rate
142`(μmax − μmin)/(μmax + μmin) < 1` on any surface, where the plain `lam`
143diverges wherever `μ > 2` — peanut reaches `μ = 6.2`. A det-based mean of
144the area factor (μ's geometric mean, exact only for conformal surfaces) is
145not enough: it under-corrects anisotropic stretching and leaves directional
146high-degree bands with amplification > 1, which surfaced as patterns going
147high-frequency and diverging as `niter` or `lmax` grew. The answer never
148depends on `jhat` — the `lamJ` term added inside `dLu` is the term divided
149back out — only the convergence rate does.
151**The correction is projected onto the band** (`.* filt` on `dLu`,
152matching algos.tex Algorithm 5's zeroing of the top coefficients). Without
153it the top two degrees iterate toward the *undiffused* `Bu` — each solve
154iteration strips a bit more of their implicit diffusion, at species-
155dependent rates, which manufactures a spurious Turing band at the band
156edge: visible on the round sphere as top-degree energy growing ~3%/step at
157`lmax 127, niter 8`. With both fixes the whole niter × geometry sweep
158converges, the spectral centroid of the pattern is resolution-independent
159(l ≈ 26 at lmax 63 and 127 alike), and `jhat: 1` is kept as the divergent
160control in the tests.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 161
162### The geometry in the operator
164`dlap = lap_g - lap_s` is applied to the current iterate at every solve
165iteration, so its transform count is what the whole step's cost scales with.
166Two formulations ship:
1681. **The flux form** (above, all three models): `lap_g u` as the weighted
169 divergence of two weighted fluxes of the sin-scaled derivatives. The
170 weights `p1, p2, q2, r` are grid arrays precomputed once per surface from
171 the embedding's θ/φ tangents
172 ([`src/geom/metric.ts`](src/geom/metric.ts)), chosen so that **every field
173 that gets analysed is a smooth function on the sphere** — the property
174 that makes spherical-harmonic analysis meaningful, and the entire
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 175 difficulty near the poles. Cost: **5 Legendre transforms** per species
176 per iteration (3 syntheses + 2 analyses; the phi flux never needs the
177 Legendre basis — `dphig` differentiates it on the grid with two FFT
178 stages, masking m past the top-degree filter — and `dthetac`/`dphic`
179 are O(nlm) coefficient shuffles). The derivation, the smoothness
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 180 argument and the fp32 error analysis are in
181 [docs/reduced-transforms.md](docs/reduced-transforms.md).
1822. **The Cartesian-gradient form** (Algorithm 4 of `docs/algos.pdf`), kept as
183 a live reference in
184 [`models/schnakenberg_alg4.m`](models/schnakenberg_alg4.m) and selectable
185 in the app: the surface gradient carried as three ambient components
186 through the inverse metric quantities `Vt*/Vp*`. Cost: **12 transforms**
187 per species per iteration. The tests hold both forms to the same answer on
188 a curved surface, and both metric formulations are precomputed and
189 uploaded for every geometry, so either kind of model runs.
191The θ-derivative machinery both forms need — the α± recurrence
192(`sin θ ∂θ Y_l^m = α⁺Y_{l+1}^m + α⁻Y_{l-1}^m`) as a coefficient-space shuffle
193feeding the existing scalar synthesis — lives in
194[`src/sht/deriv.ts`](src/sht/deriv.ts); no Legendre-derivative tables are
195required.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 196
197### `for` loops, unrolled
199A plan is a fixed list of GPU operations with no branching, which is what makes
200a timestep pure command recording — one submit, no CPU in the loop. A counted
201loop still fits: the planner
202([`src/mgpu/plan.ts`](src/mgpu/plan.ts)) unrolls it, planning the body once per
203iteration.
205Nothing else had to change for that, because numbl gives a variable one cName
206for every assignment to it: the buffer an iteration writes is the buffer the
207next one reads, which is exactly a loop-carried value. The loop variable gets no
208buffer at all — it is bound as a derived scalar to that iteration's literal, so
209a kernel reading `k` folds the number in.
211Two consequences worth stating:
213- **The bounds must be known when the model compiles.** `niter` is supplied as a
214 fixed scalar rather than a tunable one, so changing it recompiles — unlike a
215 parameter, which is a uniform. A runtime bound is refused at compile time with
216 a source position, not silently mis-compiled, and there is a test for that.
217- **Fusion survives.** numbl's inline pass recurses into loop bodies, so a line
218 inside the loop is still one kernel. It runs there with no protected names,
219 though, which means an assignment whose only visible use is later in the same
220 body can be elided — correct for a body-local temp, wrong if something outside
221 the loop wanted it. [`src/mgpu/compile.ts`](src/mgpu/compile.ts) snapshots what
222 each loop body assigns before the pass and refuses the ones that escape, so
223 that case is a compile error rather than a stale read.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 225Unrolling is exactly linear in the trip count: 19 GPU ops per species per
226iteration (6 transforms, 4 coefficient shuffles, 9 kernels), asserted in the
227tests.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 228
229## MATLAB, compiled to WebGPU
231Unchanged from turing-sphere, and it now compiles the geometry files too. numbl
232parses and lowers each function for the concrete argument types of the current
233grid; its inline pass folds single-use temps back into their consumer, so one
234line of MATLAB becomes one expression tree; and this repo emits one WGSL compute
235kernel per element-wise statement
236([`src/mgpu/wgsl.ts`](src/mgpu/wgsl.ts)). `synth` / `analys` are external
237operations whose type rules numbl learns from a `.mtoc2.js` workspace file, and
238which the backend maps onto the spherical-harmonic pipelines. Anything it cannot
239express is refused at compile time with a source position.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 241The Schnakenberg step compiles to 51 GPU operations at one solve iteration:
24216 transforms, 8 coefficient-space shuffles, 25 generated kernels, and 2
243buffer copies feeding the new state back.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 244
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 245**Transforms batch.** The expensive part of every Legendre stage is
246generating the associated Legendre values on the fly by recurrence — work
247that depends only on the grid, not on the field. `synth`/`analys` therefore
248take multiple fields, and a grouped call runs as one batched dispatch: one
249walk of the recurrence, one accumulator lane per field —
251```matlab
252[Ftu, Fpu, Ftv, Fpv] = synth(vtu, vpu, vtv, vpv); % one Legendre dispatch
253```
255The grouping is a promise of independence, never of a lane width: the
256planner ([`src/mgpu/plan.ts`](src/mgpu/plan.ts), `materializeTransforms`)
257chunks each group into whatever the device supports — one ×4 batch under the
258default WebGPU limits, or scalar dispatches with `SHT_BATCH=0` for A/B — so
259the same source runs anywhere. Ungrouped transforms that happen to sit on
260consecutive independent lines are batched the same way. Per-lane arithmetic
261is identical to the scalar kernels', so batched and scalar plans produce
262bit-identical states, asserted in the tests along with compile-time refusal
263of a group that drops one of its outputs. All 16 transforms of the step
264above land in batches, worth ~25% of the whole step (0.88 vs 1.14 ms/step at
265lmax 127, 2 iterations, on bumpy).
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 267Two consequences carried over:
269- **The step is synchronous.** WebGPU's encode path is synchronous and every
270 pipeline is built once at compile time, so a timestep is pure command
271 recording; the only `await` in the loop is the single readback per rendered
272 frame.
273- **Parameters are uniforms, not constants.** Moving a slider rewrites a small
274 buffer instead of triggering a recompile. Editing the MATLAB recompiles;
275 changing `dt` does not. `niter` is the deliberate exception, above.
277## Provenance
279- **turing-sphere**, which this is a fork of: the solver, the transforms
280 backend, the compilation path, the benchmarks and the analytic tests.
281- **Transforms:** [shtns-webgpu](https://github.com/concept-collection/shtns-webgpu) —
282 fp32 spherical harmonic transforms in WGSL compute shaders, modeled on
283 [SHTNS](https://nschaeff.bitbucket.io/shtns/). Vendored under
284 [`src/sht/`](src/sht/) (CECILL-2.1), including the f64 CPU reference transform
285 used for testing.
286- **Rendering:** three.js meshes with per-vertex colormaps, adapted from the
287 `SphereEmbedding` view in
288 [figpack](https://github.com/flatironinstitute/figpack)'s experimental
289 extension package ([`src/render/`](src/render/)). That view displays a
290 time-varying embedded geometry with fields on it, which is the same picture
291 this draws — including its sphere/surface morph, which turing-sphere had
292 dropped as having nothing to morph to.
294turing-sphere additionally carries a comparison against a native build of
295upstream SHTNS ([`bench/shtns/`](https://github.com/concept-collection/turing-sphere/tree/main/bench/shtns)).
296That is not duplicated here: the transforms are the same code, and its C-side
297transcription of the model would have to be maintained against a step this
298project intends to change.
300Because the algorithm is compiled to compute shaders, **WebGPU is required** —
301there is no CPU fallback (the f64 CPU transform remains, for tests).
303## Numerics
305- Grid: Gauss–Legendre × equispaced-φ, dealiased for the cubic reactions with
306 the `(pdeg+1)` rule: `nlat ≥ ((pdeg+1)·lmax+1)/2`, `nphi ≥ (pdeg+1)·lmax+1`
307 (rounded up to a power of two for the GPU FFT path). At the default lmax 63
308 that is a 128×256 grid.
309- Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
310 coefficients for m ≥ 0, m-major ordering.
311- fp32 transforms introduce ~1e-6 relative error per step; for pattern formation
312 from 1e-2 seeded noise this is inconsequential. The geometry goes through one
313 analysis/synthesis round trip and picks up the same round-off: the unit sphere
314 comes back with radius 1 to ~2e-5 under Dawn, ~4e-4 under SwiftShader.
315- The shipped geometries are all degree ≤ 5, far below any lmax the app offers,
316 so band-limiting removes nothing from them. A shape you write yourself may not
317 be so lucky — see the note in [`geometries/bumpy.m`](geometries/bumpy.m).
319## Desktop vs browser
321[`scripts/bench.ts`](scripts/bench.ts) runs the same thing the app runs — same
322`.m`, same generated WGSL, same transforms — from Node on desktop WebGPU (Google
323Dawn), and the app prints the command line that reproduces whatever it is
324currently simulating:
326```
327npm run bench -- --preset schnak-spots --geometry ellipsoid --lmax 63 --niter 1 \
328 --steps 2000 --seed 1 --a 0.1 --b 0.9 --D1 0.0004 --D2 0.008 --dt 0.05 \
329 --gax 1.5 --gay 1 --gaz 0.6
330```
332Copy it from under the stats line and compare the `ms/step` it reports with the
333app's. Both sides go through the one shared
334[`src/bench/runSpec.ts`](src/bench/runSpec.ts) — the app formats a run into that
335command, the benchmark parses it back — so there is no second copy of the
336defaults for the two runs to drift apart on. Geometry parameters take a `g`
337prefix (`--gwaist`) so a shape parameter can never collide with a model one.
339The app reports **two** numbers and only the first is comparable to the
340benchmark: `solver` is the batch of steps alone, waited for but not read back;
341`ms/frame` additionally carries a GPU→CPU readback per species, the
342colormapping, and the vertex upload. Those per-frame costs are fixed and do not
343shrink when the GPU gets faster, so on a quick GPU a frame can easily cost ten
344times the steps inside it. That is expected and is not the solver being slower
345in the browser.
347To attribute the gap rather than guess at it:
349```
350node scripts/compare-perf.mjs [--lmax 63] [--steps 300]
351```
353measures the same solver work in both — batched, nothing read back, no rendering
354on either side — and reports each with its CPU-encoding share, the Fourier
355stage, and the adapter. It stops you first if the two are not even the same
356device, which is a common cause of "the browser is much slower". Both sides
357resolve the geometry and the iteration count from the same constants, because
358the iteration count is unrolled into the step and a mismatch would compare two
359different amounts of work.
361The app's **Benchmark** button runs the same measurement in the page, plus the
362**ramp** — the first third of the run against the last. GPUs downclock when
363idle and an animation-paced loop leaves them idle most of every frame, so a
364large ramp means the steady-state number is limited by clocks rather than work.
366### Is it really the same computation?
368```
369node scripts/compare-env.mjs [--lmax 31] [--steps 200] [--preset schnak-spots]
370```
372runs one identical spec on the desktop and in a real browser and compares the
373final spectral state. The pipeline is deterministic given (model source,
374geometry, parameters, lmax, niter, seed, steps), so the two should agree to fp32
375round-off — not bit for bit, since GPUs differ in fused-multiply-add and other
376latitude fp32 allows. It also reports which Fourier stage each side chose, since
377FFT and DFT are genuinely different algorithms that round differently.
379Desktop WebGPU comes from the `webgpu` package (prebuilt Dawn, ~70 MB), an
380optional dependency so that an unsupported platform fails the install of that
381package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
382`--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
383(`;`-separated) passes Dawn options through.
385## Tests
387There is no second implementation of the solver to diff against, so the `.m`
388path is checked against **closed-form answers** and against **exact structural
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 389properties**. Five modules, run in both environments:
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 390
391[`test/analyticChecks.ts`](test/analyticChecks.ts) — cases whose evolution is
392known exactly, run through the whole real pipeline. All three are statements
393about the round sphere, so all three build on the sphere geometry:
395- **A** — a linear reaction leaves every mode independent, growing by exactly
396 `(1 + dt*c) / (1 + dt*D*l(l+1))` per step. Pins the transform round trip, the
397 eigenvalue mapping, the IMEX update and the state feedback at once. ~2e-7 over
398 20 steps.
399- **B** — a nonlinear reaction on a uniform field stays uniform, so each step is
400 exactly the scalar ODE map. 1.5e-8 over 25 steps.
401- **C** — a 1e-6 perturbation of the Schnakenberg fixed point follows the
402 linearized 2×2 IMEX recurrence, and `(l=24, m=7)` is confirmed unstable.
403 Looser (~4e-3) because fp32 keeps about four digits of a perturbation that
404 small.
406[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface and the loop:
408- every geometry compiles and closes; the sphere has radius 1 everywhere and is
409 **exactly degree 1** in the harmonics, which is what makes the reference case
410 exact rather than merely accurate;
411- the peanut matches its own closed-form radial profile at every grid point, and
412 **the same coefficients give the same surface on a 2× grid** — the 2× Gauss
413 latitudes share no point with the 1× ones, so agreeing there is agreeing
414 everywhere, which is what "rendered exactly, not subdivided" means;
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 415- unrolling is **exactly linear** in the trip count, and on the sphere — where
416 the geometric correction is mathematically zero — the state after 20 steps
417 stays within fp32 round-off of the 0-iteration one at 1 and 4 iterations;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 418- a runtime loop bound is refused at compile time;
419- swapping the surface mid-run leaves the spectral state untouched.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 421[`test/fluxChecks.ts`](test/fluxChecks.ts) — the six-transform flux-form
422Laplace-Beltrami scheme
423([docs/reduced-transforms.md](docs/reduced-transforms.md)):
425- on the sphere, the precomputed weights match their closed form and the
426 analysed fluxes are **exactly band-limited** (beyond-band tails at f64
427 round-off, ~1e-13), while the deliberately non-smooth control
428 `Q̃/sin θ` keeps a fat tail (~1e-2) — the discrimination the whole scheme
429 rests on;
430- on a non-axisymmetric surface, the flux tails match the Cartesian gradient
431 component's, the doc's §7.1 criterion;
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 432- the compiled op sequences add **5 Legendre transforms per species per
433 iteration against Algorithm 4's 12**, and a real simulation driven by
434 each stays within fp32 accumulation of the other.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 435
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 436[`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
437and asserts **how many kernels it compiles to**, split into the base step and
438what one solve iteration adds. That is a fusion guard: if numbl's inline pass
439stops folding, the results stay correct while every operator becomes its own
440dispatch, which is invisible in the numbers.
442[`test/transformChecks.ts`](test/transformChecks.ts) compares the WGSL transforms
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 443against shtns-webgpu's f64 CPU twin, and holds every compiled batch width to
444the scalar transforms lane by lane; a model run with `SHT_BATCH=0` must
445reproduce the batched run's state exactly.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 446
447- `npm run test:node` — under Dawn on the desktop, via `vite-node`. Needs a GPU;
448 `--skip-without-gpu` lets a machine without one say so and move on (which is
449 what CI does, since the browser suite covers the same modules).
450- `npm run test:gpu` — builds and drives headless Chrome, on SwiftShader in CI.
451 Also runs the soak. A few geometry tolerances are set by SwiftShader's fp32,
452 which is about an order of magnitude looser than Dawn's.
454Other commands:
456- `npm run bench -- --help` — the desktop benchmark.
457- `npm run bench:sht -- --help` — the transforms alone, no solver.
458- `npx vite-node scripts/diagnose-sht.ts` — when the transform tests fail on a
459 GPU, say *which* stage is wrong.
460- `npx vite-node scripts/diagnose-leg.ts [--m 0]` — read the Legendre recurrence
461 out of the production shader term by term.
462- `npx vite-node scripts/longrun-node.ts [lmax]` — run to t = 100 and confirm the
463 pattern saturates rather than decaying or diverging.
464- `node scripts/soak.mjs [steps] [lmax]` — drive the demo for many steps,
465 sampling JS heap and catching crashes.
466- `node scripts/screenshot.mjs out.png [light|dark] [minSteps]` — screenshot the
467 demo after a number of steps.
468- `node scripts/check-live.mjs [url]` — smoke-check a deployed URL.
469- `test.html?soak=<steps>&lmax=<n>` — solver-only soak with no rendering.
471## Development
473```
474npm install
475npm run dev # local dev server
476npm run build # type-check + production build to dist/
477```
479### The numbl dependency
481numbl is a local `file:../../numbl` dependency, so a sibling checkout of
482[numbl](https://github.com/flatironinstitute/numbl) is required. We use its
483compiler internals — parser, lowerer, IR, inline pass — which its package
484`exports` map does not publish, so they are reached through the `numbl-src` path
485alias in [`vite.config.ts`](vite.config.ts).
487The exact surface we depend on is written down in
488[`src/mgpu/numbl.d.ts`](src/mgpu/numbl.d.ts) and TypeScript checks against
489*that*, not against numbl's sources. This keeps this project's compiler settings
490independent of numbl's, and means a change to one of those shapes upstream
491breaks the build here with a clear diff rather than deep inside numbl's tree.
492The `For` IR node is spelled out there, since the planner now walks it.
494CI clones numbl to the sibling path that the `file:` dependency expects, pinned
495to a commit, with `--ignore-scripts` (npm runs a linked package's `prepare`
496script, and numbl's is husky). numbl's own `node_modules` are not needed: the
497slice we import is self-contained TypeScript.
499The `scripts/*.ts` entry points that touch the compiler go through `vite-node`,
500so they resolve imports exactly as the browser build does. Plain `node` cannot:
501numbl's sources import each other as `./foo.js` while the files are `.ts`.
503Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to `main`.
505## License
507CECILL-2.1 (inherited from SHTNS via shtns-webgpu, whose sources are vendored).