Seed runs from smooth random fields, and add the blob geometry
Both halves of this need randn, seeded reproducibly, over the whole
grid -- which the compiled WGSL dialect does not have and, for a
build-time quantity, has no reason to. So the shared move is to run the
non-step MATLAB through numbl's *interpreter* instead, in f64, with the
full subset available: loops, arrays, min/max, legendre, rng/randn.
1. Geometry files no longer compile to WGSL. A model's step runs every
frame and must lower to a fixed sequence of dispatches; a shape is
evaluated exactly once at build time and survives only as
coefficients, so nothing was buying. Geometry.create drops its device
argument and calls executeCode on a one-line driver that invokes
shape() with exactly the arguments its signature declares -- the same
name-matching contract the compiled path had. Everything downstream
(analys, the metric, the flux weights) is untouched.
2. tools/ is a sibling of models/ and geometries/: shared MATLAB every
interpreter run can call by file name, as on MATLAB's path. It holds
randnfunsphere and randnfun3, ported from chebfun and keeping their
upstream signatures. Not available to the models' step, which
compiles to WGSL where none of this exists.
3. geometries/blob.m is surfacefun's blob -- the sphere warped by
randnfunsphere, rescaled to [-1, 1]. amp sets how far it departs from
the sphere, lambda how fine its lobes are, and the seed parameter is
rendered as a Re-seed shape button rather than a number box, since
its value picks a draw and means nothing on its own.
4. init() is seeded from randnfun3 rather than white noise: chebfun's
randnfun3 on the surface's bounding box, restricted to the surface by
evaluating it at the grid points, the way surfacefun seeds a run.
The models say so themselves -- init(lam3, gx, gy, gz, ...) -- and
the seed wavelength is a control in the page.
A band-limited seed is fully resolved by the grid where white noise
was whatever the grid happened to alias: its energy above degree 20
is 4e-14 of the total, and the flux-form and Algorithm-4 operators
now track each other to 2.7e-6 through a run instead of 3.6e-4.
randnfun3 splits across the CPU/GPU line, and the split is forced.
Drawing the modes needs randn and a sqrt(nnz) normalization, so it is
interpreter MATLAB (a few thousand coefficients, ~5 ms); evaluating is
npts x nmodes (~6e7 terms at the default lambda), so it is a WGSL kernel
reached as an external operation the way synth is, with the coefficient
table filled in behind the call. lambda is not hidden: the plan records
which parameter the .m asked with, and the host draws from that value.
Nothing caps lambda but memory and patience -- the mode table grows to
whatever is asked for, and the only refusal is a table that could not be
built at all, reported with the mode count it wanted. Being slow is the
caller's business; freezing the browser is not, and at 0.01 (54M modes,
26 s) neither half could stay where it was. The draw moves to a worker
(13 s of synchronous interpreter time would stop the page painting), and
the GPU sum is split across 16 dispatches with the submission ended at
each one, since a browser's GPU process is shared with compositing and
one long dispatch risks the watchdog killing the device. Measured during
a seed: 731 animation frames, no stalled sample. Slices past the end of
a small table exit immediately, so a coarse lambda pays nothing.
seed() is consequently async, and the device now asks for the adapter's
full storage-buffer limit at creation so a browser's 128 MB default is
not what decides how fine lambda can be.
Two things worth knowing about lambda: it is an *absolute* length in the
surface's own units, as in chebfun, not a fraction of the surface's
size. And it is useful only down to about 2*pi/lmax (0.1 at the default
lmax 63) -- past that init's own analys discards what the grid cannot
hold, and the seed gets weaker rather than finer while costing eight
times as much per halving. Raising lmax moves that floor down. Both are
tabulated in the README.
Also fixes the numbl alias in vite.config.ts to resolve through the
symlink's realpath. dev serves modules under their real ids, so aliasing
the node_modules path gave the same file two identities and ran its side
effects twice -- harmless until something imported the interpreter,
whose builtin registry throws on the second registration. The
production build resolved it once either way, which is why the tests
never caught it.
The seed field's WGSL sum is checked against the same modes summed in
f64 on the CPU (1.98e-6 over ~1,400 terms). A kernel misreading the
packed mode table would still produce a smooth random-looking field,
which no "looks patterned" check would catch.
32 changed files+1680−228
README.mdmodified+132−15View file
@@ -36,10 +36,12 @@ function [gx, gy, gz] = shape(theta, phi, waist, stretch)
3636 end
3737 ```
3838
39-That is ordinary element-wise MATLAB and goes through the same compiler and the
40-same WGSL backend the models do. It is evaluated once on the solver's grid, and
41-then **analysed into coefficients**, which is the form everything downstream
42-uses. Two things follow from going through the coefficients rather than keeping
39+That is ordinary MATLAB. Unlike the models it is not compiled to WGSL: a shape
40+is evaluated exactly once at build time, so it runs through numbl's CPU
41+interpreter instead, in f64, with the full MATLAB subset available — loops,
42+arrays, `min`/`max`, `legendre`, seeded randomness via `rng`/`randn`. The
43+result is then **analysed into coefficients**, which is the form everything
44+downstream uses. Two things follow from going through the coefficients rather than keeping
4345 the pointwise values:
4446
4547 - **It is exactly band-limited at lmax.** The surface has as many derivatives as
@@ -52,10 +54,16 @@ the pointwise values:
5254 That is exact interpolation, not subdivision — the same argument that lets the
5355 species fields be oversampled, and it is checked directly in the tests.
5456
55-Four geometries ship: [sphere](geometries/sphere.m) (the reference case),
57+Five geometries ship: [sphere](geometries/sphere.m) (the reference case),
5658 [ellipsoid](geometries/ellipsoid.m), [peanut](geometries/peanut.m) — a dumbbell
57-whose waist is a saddle — and [bumpy](geometries/bumpy.m). Each is editable in
58-the page, with its own parameters. Changing a shape does not recompile the
59+whose waist is a saddle — [bumpy](geometries/bumpy.m), and one random one:
60+[blob](geometries/blob.m), surfacefun's blob — the sphere warped by a smooth
61+random function built from chebfun's `randnfunsphere` construction (random
62+spherical-harmonic coefficients up to degree ⌊2π/λ⌋, rescaled to [−1, 1]).
63+It is seeded, so the same seed always gives the same shape; `amp` sets how far
64+it departs from the sphere, `λ` how fine its lobes are, and **Re-seed shape**
65+draws another one. Each geometry is
66+editable in the page, with its own parameters. Changing a shape does not recompile the
5967 solver and does not disturb the run: the geometry is data whose shape in the
6068 bindings depends only on the grid, so a swap is sixteen buffer writes and the
6169 pattern carries straight on.
@@ -64,6 +72,104 @@ A **morph** slider blends the drawn surface back to the unit sphere. The
6472 parametrization is the sphere's either way, so sweeping it shows which point
6573 went where.
6674
75+### Seeding, and `tools/`
76+
77+A run starts from the uniform steady state plus a small perturbation, and that
78+perturbation is a *smooth* random field rather than white noise: chebfun's
79+[`randnfun3`](tools/randnfun3.m) on the surface's bounding box, restricted to
80+the surface by evaluating it at the grid points — the way surfacefun seeds a
81+run. Each model's `init` says so itself:
82+
83+```matlab
84+function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)
85+ f = randnfun3(lam3, gx, gy, gz);
86+ ...
87+```
88+
89+A band-limited seed is fully resolved by the grid, where white noise is
90+whatever the grid happened to alias: the tests measure its energy above degree
91+20 at 5e-14 of the total, and the flux-form and Algorithm-4 operators now
92+track each other to 3e-6 through a run instead of 4e-4. The **seed λ** control
93+sets the field's wavelength; smaller means finer features to grow from. It is
94+an *absolute* length in the surface's own units, as in chebfun — not a
95+fraction of the surface's size — so a larger surface draws more modes at the
96+same λ.
97+
98+**λ is useful down to about 2π/lmax, and no further.** A field of wavelength λ
99+on a unit-radius surface carries angular content up to degree ≈ 2π/λ, so at
100+the default lmax 63 the grid holds everything down to λ ≈ 0.1. Past that,
101+`init`'s own `analys` discards what the grid cannot represent, and the seed
102+gets *weaker* rather than finer while costing eight times as much per halving:
103+
104+| λ | 2π/λ | rms of the resolved seed | energy above l=55 | peak degree |
105+|---|---|---|---|---|
106+| 0.5 | 13 | 2.6e-2 | 1e-8 | 8 |
107+| 0.2 | 31 | 2.6e-2 | 1e-8 | 10 |
108+| 0.1 | 63 | 2.4e-2 | 0.10 | 44 |
109+| 0.05 | 126 | 1.5e-2 | 0.20 | 48 |
110+| 0.03 | 209 | 9.6e-3 | 0.27 | 63 |
111+
112+Raising lmax moves that floor down, and the seed really does get finer: at
113+lmax 127 the same λ=0.05 keeps its full amplitude (2.5e-2 against 1.5e-2 at
114+lmax 63) with its peak at degree 79 instead of pinned to the band edge, and
115+λ=0.1 becomes *fully* resolved (2e-8 of its energy in the top decile, against
116+1e-1 at lmax 63 — so even 0.1 is slightly under-resolved on the default grid).
117+
118+Note that lmax cuts both ways: it quadruples npts, so every λ also costs four
119+times as much to sum.
120+
121+**Nothing caps λ but memory and patience.** The mode table grows to whatever
122+is asked for and the only refusal is a table that could not be built at all,
123+reported with the mode count it wanted rather than silently truncated. On a
124+128×256 grid:
125+
126+| λ | modes | seed time |
127+|---|---|---|
128+| 0.05 | 480,431 | 0.26 s |
129+| 0.03 | 2,094,657 | 0.98 s |
130+| 0.02 | 6,882,185 | 3.1 s |
131+| 0.015 | 16,092,829 | 7.4 s |
132+| 0.01 | 53,574,764 | 25.6 s |
133+
134+Being slow is the caller's business; **freezing the browser is not**, and at
135+these times neither half of the work can be left where it was:
136+
137+- The draw is synchronous interpreter time — 13 s at λ=0.01 — which on the
138+ main thread stops the page painting and gets it offered up for killing. It
139+ runs on a worker instead
140+ ([`randnfun3.worker.ts`](src/mgpu/randnfun3.worker.ts)); it touches no GPU
141+ and no DOM, so nothing about it needed that thread. Measured during a seed:
142+ 731 animation frames, no stalled sample.
143+- The GPU sum is split across a fixed 16 dispatches (`randnfun3Chunks`)
144+ accumulating into the same output, and `submitYielding` ends the submission
145+ at each one. A browser's GPU process is shared with compositing, so a single
146+ submission running tens of seconds stops *every* tab painting, and one
147+ dispatch that long risks the watchdog killing the device outright. Slices
148+ past the end of a small table exit immediately, so a coarse λ pays nothing.
149+
150+The device is also asked for the adapter's full storage-buffer limit at
151+creation ([`src/sht/sht.ts`](src/sht/sht.ts)), so a browser's 128 MB default
152+is not what decides how fine λ can be. `seed()` is consequently async.
153+
154+`randnfun3` splits across the CPU/GPU line, and the split is forced rather
155+than chosen. Drawing the modes needs `randn` and a `sqrt(nnz)` normalization,
156+neither of which exists in the compiled WGSL dialect, so the draw is MATLAB in
157+[`tools/randnfun3.m`](tools/randnfun3.m) run by the interpreter — a few
158+thousand coefficients, ~5 ms. Evaluating is `npts × nmodes` (~6e7 terms at the
159+default λ), so that is a WGSL kernel
160+([`src/mgpu/randnfun3.ts`](src/mgpu/randnfun3.ts)) reached as an external
161+operation, the way `synth` is. The coefficient table is filled in behind the
162+call, as `synth` hides its Legendre matrices; λ is not hidden, and the plan
163+records which parameter the `.m` asked with so the host draws from that value.
164+
165+[`tools/`](tools/) is the shared MATLAB every interpreter run can call, by file
166+name, as on MATLAB's path — currently `randnfun3` and
167+[`randnfunsphere`](tools/randnfunsphere.m), which `blob.m` is written on. Both
168+keep their upstream signatures, including options nothing shipped uses yet
169+(`randnfunsphere`'s `'monochromatic'`), because the point of a tool is that a
170+geometry you write next can reach for it. Tools are not available to the
171+models' *step*, which compiles to WGSL where none of this exists.
172+
67173 ## The scheme, and where the geometry enters
68174
69175 It solves the N-species system
@@ -228,7 +334,8 @@ tests.
228334
229335 ## MATLAB, compiled to WebGPU
230336
231-Unchanged from turing-sphere, and it now compiles the geometry files too. numbl
337+Unchanged from turing-sphere. This is the models' path — the geometry files
338+instead run once through numbl's CPU interpreter, as above. numbl
232339 parses and lowers each function for the concrete argument types of the current
233340 grid; its inline pass folds single-use temps back into their consumer, so one
234341 line of MATLAB becomes one expression tree; and this repo emits one WGSL compute
@@ -309,12 +416,14 @@ there is no CPU fallback (the f64 CPU transform remains, for tests).
309416 - Spectral layout: SHTNS conventions — orthonormal + Condon–Shortley, complex
310417 coefficients for m ≥ 0, m-major ordering.
311418 - 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
419+ from a 1e-2 seeded perturbation this is inconsequential. The geometry goes through one
313420 analysis/synthesis round trip and picks up the same round-off: the unit sphere
314421 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).
422+- The shipped analytic geometries are all degree ≤ 5, and the random ones stay
423+ near degree 13 at their finest slider settings — far below any lmax the app
424+ offers, so band-limiting removes little to nothing from them. A shape you
425+ write yourself may not be so lucky — see the note in
426+ [`geometries/bumpy.m`](geometries/bumpy.m).
318427
319428 ## Desktop vs browser
320429
@@ -403,9 +512,10 @@ about the round sphere, so all three build on the sphere geometry:
403512 Looser (~4e-3) because fp32 keeps about four digits of a perturbation that
404513 small.
405514
406-[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface and the loop:
515+[`test/geometryChecks.ts`](test/geometryChecks.ts) — the surface, the loop, and
516+the seed:
407517
408-- every geometry compiles and closes; the sphere has radius 1 everywhere and is
518+- every geometry evaluates and closes; the sphere has radius 1 everywhere and is
409519 **exactly degree 1** in the harmonics, which is what makes the reference case
410520 exact rather than merely accurate;
411521 - the peanut matches its own closed-form radial profile at every grid point, and
@@ -416,7 +526,14 @@ about the round sphere, so all three build on the sphere geometry:
416526 the geometric correction is mathematically zero — the state after 20 steps
417527 stays within fp32 round-off of the 0-iteration one at 1 and 4 iterations;
418528 - a runtime loop bound is refused at compile time;
419-- swapping the surface mid-run leaves the spectral state untouched.
529+- swapping the surface mid-run leaves the spectral state untouched;
530+- the seed field's **WGSL sum matches the same modes summed in f64 on the
531+ CPU** (1.8e-6 over ~1,400 terms) — a kernel misreading the packed mode table
532+ would still produce a smooth random-looking field, which no "looks patterned"
533+ check would catch; the same seed redraws the same field and a different one
534+ does not; the field is band-limited (5e-14 of its energy above degree 20)
535+ with λ setting the scale; and a λ finer than the mode table holds is refused
536+ rather than silently truncated.
420537
421538 [`test/fluxChecks.ts`](test/fluxChecks.ts) — the six-transform flux-form
422539 Laplace-Beltrami scheme
geometries/blob.madded+19−0View file
@@ -0,0 +1,19 @@
1+% A random blob: the sphere, radius-modulated by a smooth random function
2+% on the sphere — surfacefun's blob, built on chebfun's randnfunsphere
3+% (tools/randnfunsphere.m).
4+%
5+% `seed` picks the draw; the same seed always gives the same blob. `scale`
6+% is the random function's wavelength, so smaller means finer lobes.
7+
8+function [gx, gy, gz] = shape(theta, phi, amp, scale, seed)
9+ rng(seed);
10+ f = randnfunsphere(scale, theta, phi);
11+ % blob.m's normalization: shift nonnegative, rescale to [-1, 1].
12+ f = f + abs(min(f));
13+ f = 2*(f/max(f)) - 1;
14+ r = 1 + amp*f;
15+ st = sin(theta);
16+ gx = r .* (st .* cos(phi));
17+ gy = r .* (st .* sin(phi));
18+ gz = r .* cos(theta);
19+end
geometries/sphere.mmodified+5−3View file
@@ -1,9 +1,11 @@
11 % The unit sphere — the reference case.
22 %
33 % A geometry file defines shape(theta, phi, ...) -> gx, gy, gz: the surface
4-% over the solver's grid (all npts x 1), compiled to WebGPU like the models.
5-% The host analyses the result into spherical-harmonic coefficients,
6-% band-limited at lmax.
4+% over the solver's grid (all npts x 1). Unlike the models it runs once, on
5+% the CPU through numbl's interpreter, so the full MATLAB subset is
6+% available — loops, arrays, min/max, legendre, seeded randomness via
7+% rng/randn. The host analyses the result into spherical-harmonic
8+% coefficients, band-limited at lmax.
79
810 function [gx, gy, gz] = shape(theta, phi)
911 st = sin(theta);
index.htmlmodified+3−0View file
@@ -239,6 +239,9 @@
239239 </label>
240240 <button id="runpause" class="primary">Run</button>
241241 <button id="benchmark">Benchmark</button>
242+ <label title="Wavelength of the smooth random field the initial condition is seeded from (chebfun's randnfun3, restricted to the surface). An absolute length in the surface's own units, as in chebfun — not a fraction of its size — so smaller means finer features to grow from. Nothing caps it but memory and patience — 0.02 takes about 3 s to seed, 0.01 about 25 s — but it is only *useful* down to about 2*pi/lmax (0.1 at lmax 63): below that the field carries more detail than the grid holds, init's analys discards it, and the seed gets weaker rather than finer. Raise lmax to go finer.">seed λ
243+ <input id="lam3" type="number" min="0" step="0.05" value="0.5">
244+ </label>
242245 <button id="reseed">Re-seed</button>
243246 <button id="resetview">Reset view</button>
244247 <button id="movietoggle" title="Export the run as an MP4 movie">Export movie</button>
models/allencahn.mmodified+3−2View file
@@ -4,8 +4,9 @@
44 %
55 % Same scheme as models/schnakenberg.m.
66
7-function [U, u] = init(noise)
8- U = analys(noise);
7+% Seeded from a smooth random field -- see models/schnakenberg.m.
8+function [U, u] = init(lam3, gx, gy, gz)
9+ U = analys(0.01 * randnfun3(lam3, gx, gy, gz));
910 u = synth(U);
1011 end
1112
models/brusselator.mmodified+4−2View file
@@ -6,8 +6,10 @@
66 % Same scheme as models/schnakenberg.m, including the grouped transforms:
77 % [a, b] = synth(x, y) runs the group as batched Legendre dispatches.
88
9-function [U, V, u, v] = init(noise, A, B)
10- [U, V] = analys(A + noise, (B / A) * ones(numel(noise), 1));
9+% Seeded from a smooth random field -- see models/schnakenberg.m.
10+function [U, V, u, v] = init(lam3, gx, gy, gz, A, B)
11+ f = randnfun3(lam3, gx, gy, gz);
12+ [U, V] = analys(A + 0.01*f, (B / A) * ones(numel(f), 1));
1113 [u, v] = synth(U, V);
1214 end
1315
models/schnakenberg.mmodified+8−2View file
@@ -15,10 +15,16 @@
1515 % docs/reduced-transforms.md, and models/schnakenberg_alg4.m
1616 % for the original form kept as a live reference.
1717
18-function [U, V, u, v] = init(noise, a, b)
18+% The uniform steady state, perturbed by a smooth random field: chebfun's
19+% randnfun3 on the surface's bounding box, restricted to the surface by
20+% evaluating it at the grid points -- the way surfacefun seeds a run. lam3
21+% is its wavelength; the draw is seeded on the host, the sum over its
22+% Fourier modes runs on the GPU (src/mgpu/randnfun3.ts).
23+function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)
24+ f = randnfun3(lam3, gx, gy, gz);
1925 us = a + b;
2026 vs = b / (us * us);
21- [U, V] = analys(us + noise, vs * ones(numel(noise), 1));
27+ [U, V] = analys(us + 0.01*f, vs * ones(numel(f), 1));
2228 [u, v] = synth(U, V);
2329 end
2430
models/schnakenberg_alg4.mmodified+5−3View file
@@ -17,11 +17,13 @@
1717 % (docs/reduced-transforms.md); this variant is kept live
1818 % for A/B comparison, in the app and in the tests.
1919
20-function [U, V, u, v] = init(noise, a, b)
20+% Seeded from a smooth random field -- see models/schnakenberg.m.
21+function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)
22+ f = randnfun3(lam3, gx, gy, gz);
2123 us = a + b;
2224 vs = b / (us * us);
23- U = analys(us + noise);
24- V = analys(vs * ones(numel(noise), 1));
25+ U = analys(us + 0.01*f);
26+ V = analys(vs * ones(numel(f), 1));
2527 u = synth(U);
2628 v = synth(V);
2729 end
scripts/bench.tsmodified+2−2View file
@@ -176,7 +176,7 @@ try {
176176 geometryParams: spec.geometryParams,
177177 niter: spec.niter,
178178 });
179- session.seed(spec.seed);
179+ await session.seed(spec.seed);
180180
181181 const plan = session.describe();
182182 const kernels = plan.step.filter((l) => l.startsWith('kernel')).length;
@@ -259,7 +259,7 @@ try {
259259 let digest = null;
260260 let state: Float32Array | null = null;
261261 if (wantDigest) {
262- session.seed(spec.seed);
262+ await session.seed(spec.seed);
263263 session.step(spec.steps);
264264 await done();
265265 state = await session.read(model.state[0]);
scripts/longrun-node.tsmodified+1−1View file
@@ -19,7 +19,7 @@ const device = await requestShtDevice().catch((e: unknown) => {
1919 throw new Error(`${errMsg(e)}\n${NO_ADAPTER_HINT}`);
2020 });
2121 const session = await ModelSession.create({ device, model, params, lmax });
22-session.seed(1);
22+await session.seed(1);
2323 console.log(`longrun — models/${model.key}.m at lmax ${lmax}, ${runtime}\n`);
2424
2525 const nsteps = Math.round(100 / params.dt);
src/geom/geometry.tsmodified+196−143View file
@@ -6,9 +6,13 @@
66 *
77 * function [gx, gy, gz] = shape(theta, phi, <parameters>)
88 *
9- * over the solver's (theta, phi) grid — the same element-wise MATLAB the models
10- * are written in, compiled by the same backend into the same kind of WGSL
11- * kernel. It is evaluated once, on the CPU's behalf, and then *analysed*: the
9+ * over the solver's (theta, phi) grid. Unlike the models it is *not* compiled
10+ * to WGSL: a model's step runs every frame and must lower to a fixed sequence
11+ * of GPU dispatches, but a shape is evaluated exactly once at build time and
12+ * survives only as coefficients. So it runs through numbl's CPU interpreter
13+ * instead, which buys the full MATLAB subset — loops, arrays, reductions,
14+ * `legendre`, seeded randomness via `rng`/`randn` — and f64 evaluation, where
15+ * the step dialect is element-wise f32. The result is then *analysed*: the
1216 * canonical geometry this project carries is the three sets of coefficients
1317 * `X`, `Y`, `Z`, one per Cartesian component of the embedding.
1418 *
@@ -28,20 +32,25 @@
2832 * The unit sphere is the case where `x`, `y`, `z` are pure degree-1 harmonics
2933 * and everything downstream reduces to turing-sphere.
3034 */
35+import { parseMFile, type FunctionStmt } from 'numbl-src/numbl-core/parser/index.ts';
36+import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
37+import {
38+ RuntimeTensor,
39+ isRuntimeTensor,
40+ type RuntimeValue,
41+} from 'numbl-src/numbl-core/runtime/types.ts';
3142 import { ShtPlan } from '../sht/sht.ts';
3243 import type { ShtConfig } from '../sht/layout.ts';
3344 import type { DerivPlan } from '../sht/deriv.ts';
3445 import { computeMetric, computeFluxMetric } from './metric.ts';
35-import { HostBuffers, ModelPlan } from '../mgpu/plan.ts';
36-import { CompiledModel, type Binding } from '../mgpu/compile.ts';
37-import { inFunction, inFunctionAsync, inModel } from '../mgpu/errors.ts';
46+import { toolFiles } from '../tools.ts';
47+import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
3848 import type { ModelParams } from '../mgpu/model.ts';
3949
4050 /** The function a geometry file must define. */
4151 export const SHAPE_FN = 'shape';
4252
4353 export interface GeometryOptions {
44- device: GPUDevice;
4554 /** The solver's transform plan — the grid the shape is evaluated on. */
4655 sht: ShtPlan;
4756 cfg: ShtConfig;
@@ -146,125 +155,90 @@ export class Geometry {
146155 }
147156
148157 /**
149- * Compile the shape file, evaluate it once on the solver grid, and reduce it
150- * to coefficients. Everything here happens at build time — a geometry never
151- * takes part in the timestep — so it reads back through the CPU freely.
158+ * Evaluate the shape file once on the solver grid and reduce it to
159+ * coefficients. Everything here happens at build time — a geometry never
160+ * takes part in the timestep — so the .m runs on the CPU (see
161+ * `evaluateShape`) and only the analysis onward touches the GPU.
152162 */
153163 static async create(opts: GeometryOptions): Promise<Geometry> {
154- const { device, sht, cfg, source, paramNames, params, deriv } = opts;
164+ const { sht, cfg, source, paramNames, params, deriv } = opts;
155165 const npts = cfg.nlat * cfg.nphi;
156- const nlm = sht.nlm;
157166
158- const bindings: Record<string, Binding> = {
159- theta: { kind: 'tensor', shape: [npts, 1] },
160- phi: { kind: 'tensor', shape: [npts, 1] },
161- npts: { kind: 'const', value: npts },
162- };
163- for (const p of paramNames) bindings[p] = { kind: 'param' };
167+ const { theta, phi } = gridAngles(sht, cfg);
168+ const raw = evaluateShape(source, paramNames, params, theta, phi, npts);
164169
165- const compiled = inModel(() => new CompiledModel(source, bindings, { npts, nlm }));
166- const fn = inFunction(SHAPE_FN, () => compiled.specialize(SHAPE_FN, 3));
167- compiled.finish();
168-
169- const host = new HostBuffers(device);
170- host.ensure('theta', npts);
171- host.ensure('phi', npts);
172-
173- const plan = await inFunctionAsync(SHAPE_FN, () =>
174- // Nothing feeds back: the three outputs are read once and the plan is
175- // thrown away.
176- ModelPlan.create(device, sht, { fn, feedback: [null, null, null] }, host),
177- );
178-
179- try {
180- const { theta, phi } = gridAngles(sht, cfg);
181- host.upload('theta', theta);
182- host.upload('phi', phi);
183- plan.setParams(params);
184-
185- const enc = device.createCommandEncoder({ label: 'geometry-shape' });
186- plan.encodeSteps(enc, 1);
187- device.queue.submit([enc.finish()]);
188-
189- const raw = await Promise.all(
190- fn.outputs.map((out) => readBuffer(device, plan, out.name, npts)),
191- );
192- // Coefficients first, then back to the grid: what the solver and the
193- // renderer both see is the band-limited surface, not the raw .m output.
194- const [X, Y, Z] = [
195- await sht.analys(raw[0]),
196- await sht.analys(raw[1]),
197- await sht.analys(raw[2]),
198- ];
199- const [x, y, z] = [
200- await sht.synth(X),
201- await sht.synth(Y),
202- await sht.synth(Z),
203- ];
170+ // Coefficients first, then back to the grid: what the solver and the
171+ // renderer both see is the band-limited surface, not the raw .m output.
172+ const [X, Y, Z] = [
173+ await sht.analys(raw[0]),
174+ await sht.analys(raw[1]),
175+ await sht.analys(raw[2]),
176+ ];
177+ const [x, y, z] = [
178+ await sht.synth(X),
179+ await sht.synth(Y),
180+ await sht.synth(Z),
181+ ];
204182
205- // Inverse metric quantities (algos.tex Algorithm 2): theta/phi
206- // derivatives of the embedding's coefficients, contracted through the
207- // inverse first fundamental form. Depends only on the geometry, so
208- // this is a one-off alongside x,y,z above, not per-step work.
209- const Xt = await deriv.dtheta(X);
210- const Xp = await deriv.dphi(X);
211- const Yt = await deriv.dtheta(Y);
212- const Yp = await deriv.dphi(Y);
213- const Zt = await deriv.dtheta(Z);
214- const Zp = await deriv.dphi(Z);
215- const { Vtx, Vty, Vtz, Vpx, Vpy, Vpz } = computeMetric(npts, Xt, Xp, Yt, Yp, Zt, Zp);
183+ // Inverse metric quantities (algos.tex Algorithm 2): theta/phi
184+ // derivatives of the embedding's coefficients, contracted through the
185+ // inverse first fundamental form. Depends only on the geometry, so
186+ // this is a one-off alongside x,y,z above, not per-step work.
187+ const Xt = await deriv.dtheta(X);
188+ const Xp = await deriv.dphi(X);
189+ const Yt = await deriv.dtheta(Y);
190+ const Yp = await deriv.dphi(Y);
191+ const Zt = await deriv.dtheta(Z);
192+ const Zp = await deriv.dphi(Z);
193+ const { Vtx, Vty, Vtz, Vpx, Vpy, Vpz } = computeMetric(npts, Xt, Xp, Yt, Yp, Zt, Zp);
216194
217- // Flux-form metric weights for the six-transform scheme, built from the
218- // *undivided* theta tangents sin(theta)*X_theta (smooth on the sphere,
219- // unlike X_theta itself) and the same X_phi as above. Also a one-off;
220- // the f64 combination happens on the CPU, rounded to f32 for upload.
221- const sXtx = await deriv.sinDtheta(X);
222- const sXty = await deriv.sinDtheta(Y);
223- const sXtz = await deriv.sinDtheta(Z);
224- const flux = computeFluxMetric(npts, sXtx, sXty, sXtz, Xp, Yp, Zp);
195+ // Flux-form metric weights for the six-transform scheme, built from the
196+ // *undivided* theta tangents sin(theta)*X_theta (smooth on the sphere,
197+ // unlike X_theta itself) and the same X_phi as above. Also a one-off;
198+ // the f64 combination happens on the CPU, rounded to f32 for upload.
199+ const sXtx = await deriv.sinDtheta(X);
200+ const sXty = await deriv.sinDtheta(Y);
201+ const sXtz = await deriv.sinDtheta(Z);
202+ const flux = computeFluxMetric(npts, sXtx, sXty, sXtz, Xp, Yp, Zp);
225203
226- // The preconditioner scale — see the Jhat field comment. The symbol
227- // matrix in the orthonormal frame is S = (1/J)[[p1,p2],[p2,q2]] with
228- // 1/J = r sin^2(theta); its entries are the bounded quantities
229- // g^tt, sin g^tp, sin^2 g^pp, so the eigenvalue extremes are clean to
230- // take over the grid. det S = 1/J^2, so the area factor comes along
231- // for free. f64 throughout.
232- let muMin = Infinity;
233- let muMax = 0;
234- let Jmin = Infinity;
235- let Jmax = 0;
236- for (let i = 0; i < cfg.nlat; i++) {
237- const ct = sht.cosTheta[i];
238- const st2 = Math.max(0, 1 - ct * ct);
239- for (let j = 0; j < cfg.nphi; j++) {
240- const k = i * cfg.nphi + j;
241- const invJ = flux.r[k] * st2;
242- const s11 = flux.p1[k] * invJ;
243- const s12 = flux.p2[k] * invJ;
244- const s22 = flux.q2[k] * invJ;
245- const mean = (s11 + s22) / 2;
246- const disc = Math.sqrt(((s11 - s22) / 2) ** 2 + s12 * s12);
247- if (mean - disc < muMin) muMin = mean - disc;
248- if (mean + disc > muMax) muMax = mean + disc;
249- const J = 1 / invJ;
250- if (J < Jmin) Jmin = J;
251- if (J > Jmax) Jmax = J;
252- }
204+ // The preconditioner scale — see the Jhat field comment. The symbol
205+ // matrix in the orthonormal frame is S = (1/J)[[p1,p2],[p2,q2]] with
206+ // 1/J = r sin^2(theta); its entries are the bounded quantities
207+ // g^tt, sin g^tp, sin^2 g^pp, so the eigenvalue extremes are clean to
208+ // take over the grid. det S = 1/J^2, so the area factor comes along
209+ // for free. f64 throughout.
210+ let muMin = Infinity;
211+ let muMax = 0;
212+ let Jmin = Infinity;
213+ let Jmax = 0;
214+ for (let i = 0; i < cfg.nlat; i++) {
215+ const ct = sht.cosTheta[i];
216+ const st2 = Math.max(0, 1 - ct * ct);
217+ for (let j = 0; j < cfg.nphi; j++) {
218+ const k = i * cfg.nphi + j;
219+ const invJ = flux.r[k] * st2;
220+ const s11 = flux.p1[k] * invJ;
221+ const s12 = flux.p2[k] * invJ;
222+ const s22 = flux.q2[k] * invJ;
223+ const mean = (s11 + s22) / 2;
224+ const disc = Math.sqrt(((s11 - s22) / 2) ** 2 + s12 * s12);
225+ if (mean - disc < muMin) muMin = mean - disc;
226+ if (mean + disc > muMax) muMax = mean + disc;
227+ const J = 1 / invJ;
228+ if (J < Jmin) Jmin = J;
229+ if (J > Jmax) Jmax = J;
253230 }
254- const Jhat = 2 / (muMin + muMax);
255-
256- return new Geometry({
257- x, y, z, X, Y, Z, Vtx, Vty, Vtz, Vpx, Vpy, Vpz,
258- p1: new Float32Array(flux.p1),
259- p2: new Float32Array(flux.p2),
260- q2: new Float32Array(flux.q2),
261- r: new Float32Array(flux.r),
262- Jhat, muMin, muMax, Jmin, Jmax,
263- });
264- } finally {
265- plan.destroy();
266- host.destroy();
267231 }
232+ const Jhat = 2 / (muMin + muMax);
233+
234+ return new Geometry({
235+ x, y, z, X, Y, Z, Vtx, Vty, Vtz, Vpx, Vpy, Vpz,
236+ p1: new Float32Array(flux.p1),
237+ p2: new Float32Array(flux.p2),
238+ q2: new Float32Array(flux.q2),
239+ r: new Float32Array(flux.r),
240+ Jhat, muMin, muMax, Jmin, Jmax,
241+ });
268242 }
269243
270244 /**
@@ -300,14 +274,15 @@ export class Geometry {
300274 }
301275 }
302276
303-/** The (theta, phi) of every grid point, flattened phi-fastest as the fields are. */
277+/** The (theta, phi) of every grid point, flattened phi-fastest as the fields
278+ * are — in f64, the precision the shape is evaluated at. */
304279 function gridAngles(
305280 sht: ShtPlan,
306281 cfg: ShtConfig,
307-): { theta: Float32Array; phi: Float32Array } {
282+): { theta: Float64Array; phi: Float64Array } {
308283 const { nlat, nphi } = cfg;
309- const theta = new Float32Array(nlat * nphi);
310- const phi = new Float32Array(nlat * nphi);
284+ const theta = new Float64Array(nlat * nphi);
285+ const phi = new Float64Array(nlat * nphi);
311286 for (let i = 0; i < nlat; i++) {
312287 const th = Math.acos(Math.max(-1, Math.min(1, sht.cosTheta[i])));
313288 for (let j = 0; j < nphi; j++) {
@@ -318,30 +293,108 @@ function gridAngles(
318293 return { theta, phi };
319294 }
320295
321-async function readBuffer(
322- device: GPUDevice,
323- plan: ModelPlan,
296+/**
297+ * Evaluate the shape file on the grid, through numbl's CPU interpreter.
298+ *
299+ * The .m keeps the same contract it had as a compiled model: it names the
300+ * arguments it wants — `theta`, `phi`, and any of the registry's parameters —
301+ * and the host supplies them by name, so their order in the signature is the
302+ * .m's own business. A one-line driver script calls `shape` with exactly the
303+ * arguments its signature declares, with those names pre-bound in the
304+ * driver's workspace.
305+ */
306+function evaluateShape(
307+ source: string,
308+ paramNames: string[],
309+ params: ModelParams,
310+ theta: Float64Array,
311+ phi: Float64Array,
312+ npts: number,
313+): [Float32Array, Float32Array, Float32Array] {
314+ const file = `${SHAPE_FN}.m`;
315+ const ast = inModel(() => parseMFile(source, file));
316+ const fn = ast.body.find(
317+ (s): s is FunctionStmt =>
318+ s.type === 'Function' && (s as FunctionStmt).name === SHAPE_FN,
319+ );
320+ if (!fn) {
321+ throw new ModelCompileError(
322+ `the geometry defines no function named '${SHAPE_FN}'`,
323+ );
324+ }
325+ if (fn.outputs.length !== 3) {
326+ throw new ModelCompileError(
327+ `'${SHAPE_FN}' must return three outputs [gx, gy, gz], not ${fn.outputs.length}`,
328+ { fn: SHAPE_FN, start: fn.span.start, end: fn.span.end },
329+ );
330+ }
331+ const known = new Set(['theta', 'phi', ...paramNames]);
332+ for (const p of fn.params) {
333+ if (!known.has(p)) {
334+ throw new ModelCompileError(
335+ `'${SHAPE_FN}' takes an argument '${p}' that is neither the grid ` +
336+ `(theta, phi) nor one of this geometry's parameters` +
337+ (paramNames.length ? ` (${paramNames.join(', ')})` : ''),
338+ { fn: SHAPE_FN, start: fn.span.start, end: fn.span.end },
339+ );
340+ }
341+ }
342+
343+ const vars: Record<string, RuntimeValue> = {
344+ theta: new RuntimeTensor(theta, [npts, 1]),
345+ phi: new RuntimeTensor(phi, [npts, 1]),
346+ };
347+ for (const name of paramNames) {
348+ const v = params[name];
349+ // Missing parameters read as 0, as ModelPlan.setParams has it.
350+ vars[name] = Number.isFinite(v) ? v : 0;
351+ }
352+
353+ const driver = `[gx__, gy__, gz__] = ${SHAPE_FN}(${fn.params.join(', ')});`;
354+ const result = inFunction(SHAPE_FN, () =>
355+ executeCode(
356+ driver,
357+ { initialVariableValues: vars, displayResults: false, implicitCwdPath: null },
358+ [...toolFiles, { name: file, source }],
359+ 'geometry-driver.m',
360+ ),
361+ );
362+
363+ return [
364+ toGridField(result.variableValues['gx__'], fn.outputs[0], npts),
365+ toGridField(result.variableValues['gy__'], fn.outputs[1], npts),
366+ toGridField(result.variableValues['gz__'], fn.outputs[2], npts),
367+ ];
368+}
369+
370+/** One returned coordinate → npts values, rounded to the transforms' f32. */
371+function toGridField(
372+ value: RuntimeValue | undefined,
324373 name: string,
325- count: number,
326-): Promise<Float32Array> {
327- const buffer = plan.buffer(name);
328- if (!buffer) {
329- throw new Error(`the geometry never assigns '${name}'`);
374+ npts: number,
375+): Float32Array {
376+ // A constant coordinate stays scalar in MATLAB; spread it over the grid.
377+ if (typeof value === 'number') return new Float32Array(npts).fill(value);
378+ if (value !== undefined && isRuntimeTensor(value)) {
379+ if (value.imag) {
380+ throw new ModelCompileError(
381+ `the geometry's '${name}' is complex; coordinates must be real`,
382+ { fn: SHAPE_FN },
383+ );
384+ }
385+ // A vector of npts values, either orientation. A 2-D reshape is refused
386+ // rather than reordered: the tensor's column-major layout would not match
387+ // the grid's phi-fastest rows.
388+ if (value.data.length === npts && value.shape.every((d) => d === 1 || d === npts)) {
389+ return new Float32Array(value.data);
390+ }
391+ throw new ModelCompileError(
392+ `the geometry's '${name}' is ${value.shape.join(' x ')}, but the grid ` +
393+ `wants one value per point (${npts} x 1)`,
394+ { fn: SHAPE_FN },
395+ );
330396 }
331- const staging = device.createBuffer({
332- label: `geometry-read-${name}`,
333- size: 4 * count,
334- usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
397+ throw new ModelCompileError(`the geometry's '${name}' is not numeric`, {
398+ fn: SHAPE_FN,
335399 });
336- try {
337- const enc = device.createCommandEncoder({ label: `geometry-read-${name}` });
338- enc.copyBufferToBuffer(buffer, 0, staging, 0, 4 * count);
339- device.queue.submit([enc.finish()]);
340- await staging.mapAsync(GPUMapMode.READ);
341- const out = new Float32Array(staging.getMappedRange().slice(0));
342- staging.unmap();
343- return out;
344- } finally {
345- staging.destroy();
346- }
347400 }
src/geom/registry.tsmodified+14−1View file
@@ -13,6 +13,7 @@ import sphereSource from '../../geometries/sphere.m?raw';
1313 import ellipsoidSource from '../../geometries/ellipsoid.m?raw';
1414 import peanutSource from '../../geometries/peanut.m?raw';
1515 import bumpySource from '../../geometries/bumpy.m?raw';
16+import blobSource from '../../geometries/blob.m?raw';
1617 import type { ParamSpec, Params } from '../mgpu/registry.ts';
1718
1819 export interface MGeometry {
@@ -67,7 +68,19 @@ const bumpy: MGeometry = {
6768 source: bumpySource,
6869 };
6970
70-export const mGeometries: MGeometry[] = [sphere, ellipsoid, peanut, bumpy];
71+const blob: MGeometry = {
72+ key: 'blob',
73+ label: 'Blob',
74+ blurb: 'The sphere warped by a smooth random function — a fresh shape per seed.',
75+ params: [
76+ { key: 'amp', label: 'amp', value: 0.5, min: 0, max: 0.8, step: 0.05 },
77+ { key: 'scale', label: 'λ', value: 1, min: 0.5, max: 3, step: 0.1 },
78+ { key: 'seed', label: 'seed', value: 1, min: 0, max: 9999, step: 1, reseed: true },
79+ ],
80+ source: blobSource,
81+};
82+
83+export const mGeometries: MGeometry[] = [sphere, ellipsoid, peanut, bumpy, blob];
7184
7285 export const mGeometryByKey = (key: string): MGeometry | undefined =>
7386 mGeometries.find((g) => g.key === key);
src/main.tsmodified+56−3View file
@@ -45,6 +45,7 @@ const elColormap = $<HTMLSelectElement>('colormap');
4545 const elRunPause = $<HTMLButtonElement>('runpause');
4646 const elBenchmark = $<HTMLButtonElement>('benchmark');
4747 const elReseed = $<HTMLButtonElement>('reseed');
48+const elLam3 = $<HTMLInputElement>('lam3');
4849 const elResetView = $<HTMLButtonElement>('resetview');
4950 const elMovieToggle = $<HTMLButtonElement>('movietoggle');
5051 const elMovieBar = $('moviebar');
@@ -275,6 +276,28 @@ function buildGeomParamInputs(): void {
275276 tag.textContent = `${geometry.key}.m`;
276277 elGeomParams.append(tag);
277278 for (const spec of geometry.params) {
279+ // A random seed picks a draw and means nothing on its own, so it gets a
280+ // button to the next one rather than a box to type a number into. The
281+ // shape changes; the simulation running on it does not restart.
282+ if (spec.reseed) {
283+ const button = document.createElement('button');
284+ button.textContent = 'Re-seed shape';
285+ button.title =
286+ `Draw another ${geometry.label.toLowerCase()} — a new random surface, ` +
287+ `leaving the pattern running on it alone.`;
288+ button.addEventListener('click', () => {
289+ const span = spec.max - spec.min;
290+ let next = geomParams[spec.key];
291+ // Never hand back the shape that is already on screen.
292+ while (next === geomParams[spec.key]) {
293+ next = spec.min + Math.floor(Math.random() * (span + 1));
294+ }
295+ geomParams[spec.key] = next;
296+ viewChange = viewChange.then(() => applyGeometry());
297+ });
298+ elGeomParams.append(button);
299+ continue;
300+ }
278301 const label = document.createElement('label');
279302 label.textContent = `${spec.label} `;
280303 const input = document.createElement('input');
@@ -372,6 +395,35 @@ elGeometry.addEventListener('change', () => {
372395 applyGeometryChoice(elGeometry.value);
373396 viewChange = viewChange.then(() => applyGeometry());
374397 });
398+// The seed field's wavelength: a uniform plus a host-side redraw, so it
399+// reseeds the run in place rather than recompiling it. Too small a value asks
400+// for more Fourier modes than the table holds, which `drawModes` refuses —
401+// report that like any other failure instead of leaving the run half-seeded.
402+elLam3.addEventListener('change', () => {
403+ const v = Number(elLam3.value);
404+ if (!Number.isFinite(v) || v <= 0) return;
405+ // Changing the wavelength redraws the field, which restarts the run — so
406+ // pause first, exactly as the Re-seed button does. Without it the reseed's
407+ // readback races the pump's own, and the two collide on the staging buffer.
408+ setRunning(false);
409+ viewChange = viewChange.then(async () => {
410+ if (!session) return;
411+ const previous = session.lam3;
412+ try {
413+ session.setLam3(v);
414+ await reseed();
415+ elErr.textContent = '';
416+ } catch (e) {
417+ // Too fine a wavelength asks for more Fourier modes than the table
418+ // holds. Put the working value back rather than leaving the run seeded
419+ // from a field that was never drawn.
420+ elErr.textContent = e instanceof Error ? e.message : String(e);
421+ session.setLam3(previous);
422+ elLam3.value = String(previous);
423+ await reseed();
424+ }
425+ });
426+});
375427 // Morph is pure rendering: no readback, no GPU work, just the vertex buffer.
376428 elMorph.addEventListener('input', () => {
377429 morph = Number(elMorph.value);
@@ -635,6 +687,7 @@ async function rebuild(): Promise<void> {
635687 geometryParams: geomParams,
636688 geometrySource: geomSource(),
637689 niter: Number(elNiter.value),
690+ lam3: Number(elLam3.value),
638691 });
639692 } catch (e) {
640693 reportCompileError(e);
@@ -642,7 +695,7 @@ async function rebuild(): Promise<void> {
642695 }
643696 if (gen !== generation) return;
644697
645- session.seed(seed);
698+ await session.seed(seed);
646699
647700 const plan = session.describe();
648701 elCompiled.textContent =
@@ -672,7 +725,7 @@ async function rebuild(): Promise<void> {
672725 async function reseed(): Promise<void> {
673726 if (!session) return;
674727 const gen = generation;
675- session.seed(seed);
728+ await session.seed(seed);
676729 if (gen !== generation) return;
677730 for (const r of ranges) {
678731 r.lo = NaN;
@@ -957,7 +1010,7 @@ async function recordMovie(): Promise<void> {
9571010 try {
9581011 // Reset the color-range smoothing as a re-seed does, so the shading
9591012 // evolves in the movie the way it did live.
960- session.seed(seed);
1013+ await session.seed(seed);
9611014 seeded = true;
9621015 for (const r of ranges) {
9631016 r.lo = NaN;
src/mgpu/externals.tsmodified+61−1View file
@@ -153,10 +153,70 @@ export function externalOpFiles(g: GridSizes): { name: string; source: string }[
153153 name: 'dphig.mtoc2.js',
154154 source: transformSource('dphig', g.npts, 1, g.npts, 1),
155155 },
156+ {
157+ // The seeded random field a model's `init` starts from
158+ // (src/mgpu/randnfun3.ts): a wavelength and the three surface
159+ // coordinates in, one value per grid point out.
160+ name: 'randnfun3.mtoc2.js',
161+ source: randnfun3Source(g),
162+ },
156163 ];
157164 }
158165
166+/** Source for `randnfun3`'s `.mtoc2.js`: `f = randnfun3(lambda, x, y, z)`. */
167+function randnfun3Source(g: GridSizes): string {
168+ return `
169+exports.name = "randnfun3";
170+
171+exports.transfer = function (argTypes, nargout) {
172+ if (argTypes.length !== 4) {
173+ throw new Error(
174+ "randnfun3 takes a wavelength and the three surface coordinates -- " +
175+ "randnfun3(lambda, gx, gy, gz) -- got " + argTypes.length + " argument(s)"
176+ );
177+ }
178+ if (nargout > 1) {
179+ throw new Error("randnfun3 returns one value, but " + nargout + " were requested");
180+ }
181+ var lam = argTypes[0];
182+ if (!lam || lam.kind !== "Numeric" || lam.isComplex) {
183+ throw new Error("randnfun3's wavelength must be a real number");
184+ }
185+ var ls = lam.shape;
186+ if (!ls || ls.length !== 2 || ls[0] !== 1 || ls[1] !== 1) {
187+ throw new Error(
188+ "randnfun3's wavelength must be a single number, not a " +
189+ (ls ? ls.join("x") : "unknown shape") + " array"
190+ );
191+ }
192+ var names = ["gx", "gy", "gz"];
193+ for (var i = 1; i < 4; i++) {
194+ var a = argTypes[i];
195+ if (!a || a.kind !== "Numeric" || a.isComplex) {
196+ throw new Error("randnfun3 requires real numeric arrays (" + names[i - 1] + ")");
197+ }
198+ var s = a.shape;
199+ if (!s || s.length !== 2 || s[0] !== ${g.npts} || s[1] !== 1) {
200+ throw new Error(
201+ "randnfun3 evaluates on the grid, so " + names[i - 1] +
202+ " must be ${g.npts}x1, not " + (s ? s.join("x") : "unknown shape")
203+ );
204+ }
205+ }
206+ return [${numericType(g.npts, 1)}];
207+};
208+
209+// Never called: this project executes the IR on WebGPU and emits no C.
210+exports.emit = function () {
211+ throw new Error("randnfun3: no C backend (this runs on WebGPU)");
212+};
213+exports.cBody = function () {
214+ return "";
215+};
216+`;
217+}
218+
159219 /** Names the WGSL backend must implement as GPU encodes rather than kernels. */
160220 export const EXTERNAL_OPS = new Set([
161- 'synth', 'analys', 'dtheta', 'dphi', 'dthetac', 'dphic', 'dphig',
221+ 'synth', 'analys', 'dtheta', 'dphi', 'dthetac', 'dphic', 'dphig', 'randnfun3',
162222 ]);
src/mgpu/model.tsmodified+28−7View file
@@ -20,7 +20,8 @@
2020 import { ShtPlan } from '../sht/sht.ts';
2121 import type { DerivPlan } from '../sht/deriv.ts';
2222 import { lmIndex, type ShtConfig } from '../sht/layout.ts';
23-import { HostBuffers, ModelPlan } from './plan.ts';
23+import { HostBuffers, ModelPlan, type Randnfun3Lambda } from './plan.ts';
24+import { MODE_BUFFER } from './randnfun3.ts';
2425 import { inFunction, inFunctionAsync, inModel } from './errors.ts';
2526 import { CompiledModel, type Binding } from './compile.ts';
2627
@@ -208,6 +209,10 @@ export class GpuModel {
208209 // so swapping the surface updates it with no recompile. The session
209210 // folds the current geometry's value into every setParams call.
210211 bindings['jhat'] = { kind: 'param' };
212+ // The wavelength of the seeded random field (src/mgpu/randnfun3.ts).
213+ // A uniform like jhat, not a const: changing it redraws the field
214+ // without recompiling the step.
215+ bindings['lam3'] = { kind: 'param' };
211216 }
212217 for (const s of state) bindings[s] = { kind: 'tensor', shape: [2, nlm] };
213218 for (const p of paramNames) bindings[p] = { kind: 'param' };
@@ -324,12 +329,28 @@ export class GpuModel {
324329 this.#jhat = geometry.Jhat;
325330 }
326331
327- /** Upload the seeded perturbation and run `init`. */
328- init(noise: Float32Array): void {
329- this.#host.upload('noise', noise);
330- const enc = this.#device.createCommandEncoder({ label: 'mgpu-init' });
331- this.#initPlan.encodeSteps(enc, 1);
332- this.#device.queue.submit([enc.finish()]);
332+ /** The wavelength this model's `init` asked `randnfun3` for, or null if it
333+ * seeds some other way. The session resolves it and draws the modes. */
334+ get randnfun3Lambda(): Randnfun3Lambda | null {
335+ return this.#initPlan.randnfun3Lambda;
336+ }
337+
338+ /**
339+ * Upload the seeded initial data and run `init`.
340+ *
341+ * Both inputs are optional in the sense that a .m uses one or the other:
342+ * `modes` is the random field's coefficient table for a model that calls
343+ * `randnfun3`, `noise` the plain grid field for one that takes `noise`
344+ * directly (the analytic test models inject exact initial conditions that
345+ * way). Only what the plan actually bound is uploaded.
346+ */
347+ async init(noise: Float32Array, modes: Float32Array | null): Promise<void> {
348+ if (this.#host.get('noise')) this.#host.upload('noise', noise);
349+ // Sized to the wavelength, so this may reallocate and rebind.
350+ if (modes) this.#initPlan.uploadRandnfun3Table(this.#host, modes);
351+ // Submitted in pieces: a fine seed wavelength makes the mode sum long
352+ // enough that one submission would stall the browser's compositor.
353+ await this.#initPlan.submitYielding('mgpu-init');
333354 this.#lastRan = 'init';
334355 }
335356
src/mgpu/numbl.d.tsmodified+90−5View file
@@ -1,10 +1,12 @@
11 /**
22 * The numbl compiler surface this project depends on.
33 *
4- * We reach past numbl's published entry points into its JIT internals (parser,
5- * lowerer, IR, inline pass), which its package `exports` map does not expose.
6- * Those imports resolve through the `numbl-src` alias in vite.config.ts; these
7- * declarations are what TypeScript checks against.
4+ * We reach past numbl's published entry points into its internals — the JIT
5+ * side (parser, lowerer, IR, inline pass) that compiles the models, and the
6+ * interpreter side (executeCode, runtime values) that evaluates the
7+ * geometries — which its package `exports` map does not expose. Those imports
8+ * resolve through the `numbl-src` alias in vite.config.ts; these declarations
9+ * are what TypeScript checks against.
810 *
911 * Declaring the surface here rather than type-checking numbl's sources
1012 * directly keeps this project's compiler settings independent of numbl's, and
@@ -192,13 +194,96 @@ declare module 'numbl-src/numbl-core/jit/lowering/ir.ts' {
192194 }
193195
194196 declare module 'numbl-src/numbl-core/parser/index.ts' {
197+ export interface ParseSpan {
198+ start: number;
199+ end: number;
200+ }
201+
202+ /** The one parse-tree node this project inspects (src/geom/geometry.ts,
203+ * finding `shape` and its argument names). */
204+ export interface FunctionStmt {
205+ type: 'Function';
206+ name: string;
207+ params: string[];
208+ outputs: string[];
209+ span: ParseSpan;
210+ }
211+
212+ /** Any other statement in a file's body — opaque to this project. Its
213+ * `type` is some other literal; narrowing to FunctionStmt goes through an
214+ * explicit type guard rather than the discriminant. */
215+ export interface OtherParseStmt {
216+ type: string;
217+ span: ParseSpan;
218+ }
219+
220+ export type Stmt = FunctionStmt | OtherParseStmt;
221+
195222 export interface AbstractSyntaxTree {
196- body: unknown[];
223+ body: Stmt[];
197224 }
198225 export function parseMFile(input: string, fileName?: string): AbstractSyntaxTree;
199226 export class SyntaxError extends Error {}
200227 }
201228
229+declare module 'numbl-src/numbl-core/runtime/types.ts' {
230+ /** A numeric array: f64 data in column-major order, with its shape. */
231+ export class RuntimeTensor {
232+ readonly kind: 'tensor';
233+ data: Float64Array;
234+ /** Present iff the value is complex. */
235+ imag: Float64Array | undefined;
236+ shape: number[];
237+ constructor(data: Float64Array, shape: number[], imag?: Float64Array);
238+ }
239+
240+ /** Every other value kind the interpreter can hold, collapsed. */
241+ export interface OtherRuntimeValue {
242+ readonly kind: string;
243+ }
244+
245+ export type RuntimeValue =
246+ | number
247+ | boolean
248+ | string
249+ | RuntimeTensor
250+ | OtherRuntimeValue;
251+
252+ export function isRuntimeTensor(value: RuntimeValue): value is RuntimeTensor;
253+}
254+
255+declare module 'numbl-src/numbl-core/executeCode.ts' {
256+ import type { RuntimeValue } from 'numbl-src/numbl-core/runtime/types.ts';
257+
258+ export interface ExecOptions {
259+ /** Variables pre-bound in the script's workspace before it runs. */
260+ initialVariableValues?: Record<string, RuntimeValue>;
261+ displayResults?: boolean;
262+ onOutput?: (text: string) => void;
263+ /** null opts out of scanning a working directory for .m files. */
264+ implicitCwdPath?: string | null;
265+ }
266+
267+ export interface ExecWorkspaceFile {
268+ name: string;
269+ source: string;
270+ }
271+
272+ export interface ExecResult {
273+ output: string[];
274+ /** The script's workspace after it ran. */
275+ variableValues: Record<string, RuntimeValue>;
276+ }
277+
278+ /** Run a script through numbl's interpreter (with its JS-JIT), CPU-side. */
279+ export function executeCode(
280+ source: string,
281+ options?: ExecOptions,
282+ workspaceFiles?: ExecWorkspaceFile[],
283+ mainFileName?: string,
284+ ): ExecResult;
285+}
286+
202287 declare module 'numbl-src/numbl-core/jit/index.ts' {
203288 import type { AbstractSyntaxTree } from 'numbl-src/numbl-core/parser/index.ts';
204289 import type { IRProgram, IRFunc, Span } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
src/mgpu/plan.tsmodified+291−12View file
@@ -21,6 +21,13 @@ import { ShtPlan, type ShtBinding, type ShtBatchBinding, type ShtDphigBinding }
2121 import { DerivPlan, type DerivBinding } from '../sht/deriv.ts';
2222 import type { CompiledFunction } from './compile.ts';
2323 import { EXTERNAL_OPS } from './externals.ts';
24+import {
25+ MODE_BUFFER,
26+ INITIAL_MODES,
27+ modeTableLength,
28+ randnfun3Chunks,
29+ randnfun3WGSL,
30+} from './randnfun3.ts';
2431 import {
2532 buildKernel,
2633 UnsupportedOnGpu,
@@ -96,6 +103,35 @@ export class HostBuffers {
96103 return this.#slots.get(name);
97104 }
98105
106+ /**
107+ * Replace a slot's buffer with a larger one. Only for buffers whose size is
108+ * not fixed by the grid — the randnfun3 mode table, which grows with the
109+ * wavelength asked for. The caller must rebuild any bind group holding the
110+ * old buffer; it is destroyed here.
111+ */
112+ resize(name: string, count: number): Slot {
113+ const existing = this.#slots.get(name);
114+ if (!existing) throw new Error(`resize: no buffer named '${name}'`);
115+ if (count <= existing.count) return existing;
116+ existing.buffer.destroy();
117+ const slot = { buffer: makeBuffer(this.#device, `mgpu-${name}`, count), count };
118+ this.#slots.set(name, slot);
119+ return slot;
120+ }
121+
122+ /** Upload into the front of a slot, leaving any tail as it was. For a
123+ * variable-length payload in a buffer sized to its high-water mark. */
124+ uploadInto(name: string, data: Float32Array): void {
125+ const slot = this.#slots.get(name);
126+ if (!slot) throw new Error(`uploadInto: no buffer named '${name}'`);
127+ if (data.length > slot.count) {
128+ throw new Error(
129+ `uploadInto '${name}': ${data.length} elements into a ${slot.count}-element buffer`,
130+ );
131+ }
132+ this.#device.queue.writeBuffer(slot.buffer, 0, data as Float32Array<ArrayBuffer>);
133+ }
134+
99135 /** Upload initial data for a host binding. */
100136 upload(name: string, data: Float32Array): void {
101137 const slot = this.#slots.get(name);
@@ -124,6 +160,9 @@ type Op =
124160 /** Set when the kernel had to write to scratch because its output
125161 * aliases one of its inputs; copied back after the dispatch. */
126162 copyBack?: { from: GPUBuffer; to: GPUBuffer; bytes: number };
163+ /** End the submission here when run through `submitYielding`, so the
164+ * GPU is handed back between chunks of a long seed. */
165+ yieldAfter?: boolean;
127166 }
128167 | { kind: 'synth' | 'analys'; binding: ShtBinding; label: string }
129168 | { kind: 'synth-batch' | 'analys-batch'; binding: ShtBatchBinding; labels: string[] }
@@ -295,6 +334,9 @@ async function makePipeline(
295334 export class ModelPlan {
296335 /** Scalar parameter names, in the order the params buffer expects them. */
297336 readonly paramNames: string[];
337+ /** The wavelength this plan's `randnfun3` call asked for, or null if it
338+ * makes none. The host draws the coefficient table from it. */
339+ readonly randnfun3Lambda: Randnfun3Lambda | null;
298340
299341 #device: GPUDevice;
300342 #sht: ShtPlan;
@@ -303,6 +345,7 @@ export class ModelPlan {
303345 #owned: GPUBuffer[];
304346 #paramBuf: GPUBuffer;
305347 #paramData: Float32Array;
348+ #rebindRandnfun3: ((table: GPUBuffer) => void) | null;
306349 /** Public name -> buffer, for uploading initial state and reading results. */
307350 #byName: Map<string, Slot>;
308351
@@ -316,6 +359,8 @@ export class ModelPlan {
316359 paramBuf: GPUBuffer;
317360 paramData: Float32Array;
318361 paramNames: string[];
362+ randnfun3Lambda: Randnfun3Lambda | null;
363+ rebindRandnfun3: ((table: GPUBuffer) => void) | null;
319364 }) {
320365 this.#device = init.device;
321366 this.#sht = init.sht;
@@ -326,6 +371,30 @@ export class ModelPlan {
326371 this.#paramBuf = init.paramBuf;
327372 this.#paramData = init.paramData;
328373 this.paramNames = init.paramNames;
374+ this.randnfun3Lambda = init.randnfun3Lambda;
375+ this.#rebindRandnfun3 = init.rebindRandnfun3;
376+ }
377+
378+ /**
379+ * Point the randnfun3 dispatch at a mode table big enough for `data`,
380+ * growing the buffer if this wavelength needs more modes than the last one,
381+ * and upload it.
382+ */
383+ uploadRandnfun3Table(host: HostBuffers, data: Float32Array): void {
384+ const slot = host.get(MODE_BUFFER);
385+ if (!slot || !this.#rebindRandnfun3) return;
386+ if (data.length > slot.count) {
387+ const max = this.#device.limits.maxStorageBufferBindingSize;
388+ if (4 * data.length > max) {
389+ throw new Error(
390+ `randnfun3: this wavelength needs a ${(4 * data.length / 1e6).toFixed(0)} MB ` +
391+ `mode table, past this device's ${(max / 1e6).toFixed(0)} MB limit ` +
392+ `on a single buffer. Use a larger lambda.`,
393+ );
394+ }
395+ this.#rebindRandnfun3(host.resize(MODE_BUFFER, data.length).buffer);
396+ }
397+ host.uploadInto(MODE_BUFFER, data);
329398 }
330399
331400 static async create(
@@ -375,6 +444,14 @@ export class ModelPlan {
375444 usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
376445 });
377446
447+ /** Set when the .m calls `randnfun3`: which wavelength it asked for, so
448+ * the host draws the coefficient table the kernel reads from exactly
449+ * that value (src/mgpu/randnfun3.ts). */
450+ let randnfun3Lambda: Randnfun3Lambda | null = null;
451+ /** Rebuilds the randnfun3 dispatch's bind group after the mode table is
452+ * reallocated for a finer wavelength. */
453+ let rebindRandnfun3: ((table: GPUBuffer) => void) | null = null;
454+
378455 const planned: Planned[] = [];
379456 for (const stmt of fn.body) {
380457 await planStatement(stmt);
@@ -413,6 +490,7 @@ export class ModelPlan {
413490
414491 return new ModelPlan({
415492 device, sht, deriv, ops, byName, owned, paramBuf, paramData, paramNames,
493+ randnfun3Lambda, rebindRandnfun3,
416494 });
417495
418496 async function planStatement(stmt: IRStmt): Promise<void> {
@@ -457,14 +535,19 @@ export class ModelPlan {
457535
458536 const ext = externalCall(stmt);
459537 if (ext) {
460- const argSlot = slots.get(ext.argCName);
538+ if (ext.name === 'randnfun3') {
539+ await planRandnfun3(stmt, ext.args, dest);
540+ return;
541+ }
542+ const arg = ext.args[0] as IRExpr & { kind: 'Var' };
543+ const argSlot = slots.get(arg.cName);
461544 if (!argSlot) {
462545 throw new UnsupportedOnGpu(
463- `'${ext.name}' reads '${ext.argName}', which has no buffer`,
546+ `'${ext.name}' reads '${arg.name}', which has no buffer`,
464547 stmt.span,
465548 );
466549 }
467- const label = `${stmt.name} = ${ext.name}(${ext.argName})`;
550+ const label = `${stmt.name} = ${ext.name}(${arg.name})`;
468551 if (ext.name === 'dphig') {
469552 // Grid -> grid, staged through the plan's fm scratch; safe even
470553 // in place, so no aliasing guard is needed.
@@ -504,7 +587,7 @@ export class ModelPlan {
504587 // silently reroute.
505588 if (argSlot.buffer === dest.buffer) {
506589 throw new UnsupportedOnGpu(
507- `'${stmt.name} = ${ext.name}(${ext.argName})' reads and ` +
590+ `'${stmt.name} = ${ext.name}(${arg.name})' reads and ` +
508591 `writes the same buffer; assign to a new name instead`,
509592 stmt.span,
510593 );
@@ -658,6 +741,115 @@ export class ModelPlan {
658741 }
659742 }
660743
744+ /**
745+ * `f = randnfun3(lambda, gx, gy, gz)`: the seeded random field, summed
746+ * over its Fourier modes at every surface point.
747+ *
748+ * One dispatch, one thread per point. The coefficient table is not an
749+ * argument — it is a host buffer this plan binds and the host refills per
750+ * seed, the way `synth` reads Legendre matrices the .m never names. What
751+ * the .m *does* choose is the wavelength, which is recorded here so the
752+ * host draws the table for exactly that value.
753+ */
754+ async function planRandnfun3(
755+ stmt: Assign,
756+ args: IRExpr[],
757+ dest: Slot,
758+ ): Promise<void> {
759+ const lam = args[0];
760+ const lambda: Randnfun3Lambda | null =
761+ lam.kind === 'NumLit'
762+ ? { kind: 'const', value: lam.value }
763+ : lam.kind === 'Var' && paramSlots.has(lam.cName)
764+ ? { kind: 'param', name: lam.name }
765+ : null;
766+ if (!lambda) {
767+ throw new UnsupportedOnGpu(
768+ `randnfun3's wavelength is drawn on the host before the step runs, ` +
769+ `so it must be a number or a model parameter — not a value ` +
770+ `computed on the GPU`,
771+ stmt.span,
772+ );
773+ }
774+ if (randnfun3Lambda && !sameLambda(randnfun3Lambda, lambda)) {
775+ throw new UnsupportedOnGpu(
776+ `this function calls randnfun3 with two different wavelengths; ` +
777+ `one coefficient table is drawn per plan, so only one is supported`,
778+ stmt.span,
779+ );
780+ }
781+ randnfun3Lambda = lambda;
782+
783+ const points = args.slice(1).map((a) => {
784+ const v = a as IRExpr & { kind: 'Var' };
785+ const slot = slots.get(v.cName);
786+ if (!slot) {
787+ throw new UnsupportedOnGpu(
788+ `randnfun3 reads '${v.name}', which has no buffer`,
789+ stmt.span,
790+ );
791+ }
792+ return { slot, name: v.name };
793+ });
794+
795+ const modes = host.ensure(MODE_BUFFER, modeTableLength(INITIAL_MODES));
796+ const label =
797+ `${stmt.name} = randnfun3(${
798+ lambda.kind === 'const' ? lambda.value : lambda.name
799+ }, ${points.map((p) => p.name).join(', ')})`;
800+
801+ const bindGroupLayout = device.createBindGroupLayout({
802+ label: 'mgpu-randnfun3',
803+ entries: [0, 1, 2, 3, 4].map((binding) => ({
804+ binding,
805+ visibility: GPUShaderStage.COMPUTE,
806+ buffer: { type: binding === 0 ? ('storage' as const) : ('read-only-storage' as const) },
807+ })),
808+ });
809+ // The table is sized to whatever wavelength is actually asked for, so a
810+ // finer one reallocates it — and with it these bind groups, which are
811+ // the only things holding the old buffer.
812+ const bind = (table: GPUBuffer): GPUBindGroup =>
813+ device.createBindGroup({
814+ layout: bindGroupLayout,
815+ entries: [
816+ { binding: 0, resource: { buffer: dest.buffer } },
817+ ...points.map((p, i) => ({
818+ binding: i + 1,
819+ resource: { buffer: p.slot.buffer },
820+ })),
821+ { binding: 4, resource: { buffer: table } },
822+ ],
823+ });
824+
825+ // One dispatch per slice of the mode table — see randnfun3Chunks. Each
826+ // reads the same table and accumulates into the same output, so they
827+ // share a bind group and differ only in their compiled slice index.
828+ const ops: (Op & { kind: 'kernel' })[] = [];
829+ for (let chunk = 0; chunk < randnfun3Chunks; chunk++) {
830+ const chunkLabel = `${label} [${chunk + 1}/${randnfun3Chunks}]`;
831+ const op = {
832+ kind: 'kernel' as const,
833+ pipeline: await makePipeline(
834+ device,
835+ randnfun3WGSL(dest.count, chunk),
836+ chunkLabel,
837+ bindGroupLayout,
838+ ),
839+ bindGroup: bind(modes.buffer),
840+ count: dest.count,
841+ label: chunkLabel,
842+ yieldAfter: true,
843+ };
844+ ops.push(op);
845+ planned.push(op);
846+ }
847+ rebindRandnfun3 = (table: GPUBuffer): void => {
848+ const group = bind(table);
849+ for (const op of ops) op.bindGroup = group;
850+ };
851+ }
852+
661853 /**
662854 * Unroll a counted loop into the op sequence.
663855 *
@@ -739,13 +931,64 @@ export class ModelPlan {
739931 return this.#byName.get(name)?.count;
740932 }
741933
934+ /**
935+ * Run one pass of this plan, submitting in pieces so the GPU is not held for
936+ * the whole of it.
937+ *
938+ * For `init` only, and only because the seed field's mode sum can be huge:
939+ * at a fine wavelength the dispatches add up to tens of seconds, and a
940+ * browser's GPU process is shared with compositing, so one submission that
941+ * long stops the whole browser painting — the user's tabs included. Ops
942+ * marked `yieldAfter` (the randnfun3 chunks) end their submission and give
943+ * the queue back before the next one is recorded, which turns a freeze into
944+ * a wait. Everything else is recorded exactly as `encodeSteps` would.
945+ */
946+ async submitYielding(label: string): Promise<void> {
947+ let encoder = this.#device.createCommandEncoder({ label });
948+ let any = false;
949+ for (const group of this.#yieldGroups()) {
950+ if (any) {
951+ // Let the queue drain, then hand the event loop back, so compositing
952+ // and input get a turn between chunks.
953+ await this.#device.queue.onSubmittedWorkDone();
954+ await new Promise((r) => setTimeout(r, 0));
955+ encoder = this.#device.createCommandEncoder({ label });
956+ }
957+ this.#encodeOps(encoder, group);
958+ this.#device.queue.submit([encoder.finish()]);
959+ any = true;
960+ }
961+ if (!any) {
962+ this.#encodeOps(encoder, []);
963+ this.#device.queue.submit([encoder.finish()]);
964+ }
965+ }
966+
967+ /** The op list split at every `yieldAfter` boundary. */
968+ *#yieldGroups(): Generator<Op[]> {
969+ let group: Op[] = [];
970+ for (const op of this.#ops) {
971+ group.push(op);
972+ if (op.kind === 'kernel' && op.yieldAfter) {
973+ yield group;
974+ group = [];
975+ }
976+ }
977+ if (group.length) yield group;
978+ }
979+
742980 /**
743981 * Record `steps` timesteps. Synchronous: no awaits, no readback. All of the
744982 * ops share one compute pass, which WebGPU executes in submission order
745983 * with a barrier between dispatches.
746984 */
747985 encodeSteps(encoder: GPUCommandEncoder, steps: number): void {
748- for (let s = 0; s < steps; s++) {
986+ for (let s = 0; s < steps; s++) this.#encodeOps(encoder, this.#ops);
987+ }
988+
989+ /** Record one pass over `ops` into `encoder`. */
990+ #encodeOps(encoder: GPUCommandEncoder, ops: Op[]): void {
991+ {
749992 let pass: GPUComputePassEncoder | null = null;
750993 const inPass = (): GPUComputePassEncoder => {
751994 if (!pass) pass = encoder.beginComputePass({ label: 'mgpu-step' });
@@ -757,7 +1000,7 @@ export class ModelPlan {
7571000 pass = null;
7581001 }
7591002 };
760- for (const op of this.#ops) {
1003+ for (const op of ops) {
7611004 switch (op.kind) {
7621005 case 'kernel': {
7631006 const p = inPass();
@@ -847,22 +1090,58 @@ export class ModelPlan {
8471090 }
8481091 }
8491092
850-/** `x = synth(y)` / `x = analys(y)` -> the call's name and argument. */
1093+/**
1094+ * `x = synth(y)` / `x = randnfun3(lam, gx, gy, gz)` -> the call's name and
1095+ * arguments.
1096+ *
1097+ * Every external op but `randnfun3` takes exactly one array; `randnfun3`
1098+ * takes a wavelength and the three surface coordinates. Its wavelength may
1099+ * be a literal, so arguments are returned as expressions and the caller
1100+ * decides which it needs as a buffer.
1101+ */
8511102 function externalCall(
8521103 stmt: Assign,
853-): { name: string; argCName: string; argName: string } | null {
1104+): { name: string; args: IRExpr[] } | null {
8541105 const e = stmt.expr;
8551106 if (e.kind !== 'Call' || !EXTERNAL_OPS.has(e.name)) return null;
856- if (e.args.length !== 1 || e.args[0].kind !== 'Var') {
1107+ const arity = e.name === 'randnfun3' ? 4 : 1;
1108+ if (e.args.length !== arity) {
8571109 throw new UnsupportedOnGpu(
858- `'${e.name}' must be applied to a single variable`,
1110+ arity === 1
1111+ ? `'${e.name}' must be applied to a single variable`
1112+ : `'${e.name}' takes ${arity} arguments, got ${e.args.length}`,
8591113 stmt.span,
8601114 );
8611115 }
862- const arg = e.args[0];
863- return { name: e.name, argCName: arg.cName, argName: arg.name };
1116+ // Only the wavelength may be something other than a plain variable.
1117+ for (let i = e.name === 'randnfun3' ? 1 : 0; i < e.args.length; i++) {
1118+ if (e.args[i].kind !== 'Var') {
1119+ throw new UnsupportedOnGpu(
1120+ `'${e.name}' must be applied to variables, not expressions`,
1121+ stmt.span,
1122+ );
1123+ }
1124+ }
1125+ return { name: e.name, args: e.args };
8641126 }
8651127
1128+/** A `randnfun3` wavelength argument: a literal, or the parameter to read it
1129+ * from when the host fills the coefficient table. */
1130+export type Randnfun3Lambda =
1131+ | { kind: 'const'; value: number }
1132+ | { kind: 'param'; name: string };
1133+
1134+const sameLambda = (a: Randnfun3Lambda, b: Randnfun3Lambda): boolean =>
1135+ a.kind === 'const' && b.kind === 'const'
1136+ ? a.value === b.value
1137+ : a.kind === 'param' && b.kind === 'param' && a.name === b.name;
1138+
1139+/** The wavelength value a plan's `randnfun3` call resolves to. */
1140+export const resolveLambda = (
1141+ lambda: Randnfun3Lambda,
1142+ params: Record<string, number>,
1143+): number => (lambda.kind === 'const' ? lambda.value : params[lambda.name]);
1144+
8661145 function collectTensorVars(e: IRExpr, visit: (cName: string) => void): void {
8671146 const walk = (x: IRExpr): void => {
8681147 switch (x.kind) {
src/mgpu/randnfun3.tsadded+294−0View file
@@ -0,0 +1,294 @@
1+/**
2+ * `randnfun3` — a smooth random function in 3D, evaluated at the surface.
3+ *
4+ * chebfun's randnfun3 is a random trig series on a box: a few thousand
5+ * Fourier modes with independent normal coefficients, confined to a ball for
6+ * isotropy and normalized to unit variance. Restricting it to a surface is
7+ * just evaluating it at the surface's points, which is what a model's `init`
8+ * wants for a seeded initial condition (surfacefun seeds exactly this way).
9+ *
10+ * The work splits in two, and the split is forced rather than chosen:
11+ *
12+ * - **Drawing the modes needs `randn`**, which the compiled WGSL dialect has
13+ * no counterpart for, and `sqrt(nnz)` normalization, which is a reduction.
14+ * Both are a few lines of MATLAB, so the draw lives in
15+ * `tools/randnfun3.m` and runs in numbl's interpreter — a few thousand
16+ * numbers, ~5 ms.
17+ * - **Evaluating is npts x nmodes**, ~6e7 terms at the default lambda. That
18+ * is the whole cost, and it is what this file's kernel does on the GPU.
19+ *
20+ * So the .m calls `f = randnfun3(lambda, gx, gy, gz)` — chebfun's signature,
21+ * lambda in and values out — and the coefficient table is filled in behind it
22+ * by the host, the way `synth` hides its Legendre matrices. lambda is not
23+ * decorative: the plan records which parameter the .m passed, and the host
24+ * draws the table from *that* parameter's value (src/mgpu/plan.ts,
25+ * `randnfun3Lambda`), so changing it in the .m changes the field.
26+ */
27+import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
28+import { isRuntimeTensor } from 'numbl-src/numbl-core/runtime/types.ts';
29+import { toolFiles } from '../tools.ts';
30+
31+/**
32+ * Dispatches the mode sum is split across.
33+ *
34+ * lambda is an absolute length and the mode count goes as its inverse cube,
35+ * so halving lambda costs eight times the work — there is no natural ceiling
36+ * to put on that, and nothing in the method breaks as it grows. It just gets
37+ * slower, which is the caller's business. What is *not* the caller's business
38+ * is a browser's GPU-process watchdog, which kills the device outright when a
39+ * single dispatch runs too long; a fine wavelength would otherwise turn "this
40+ * takes a while" into "device lost".
41+ *
42+ * So the sum is split into a fixed number of dispatches, each covering its own
43+ * slice of the table and accumulating into the same output. The count is fixed
44+ * at plan time (the op sequence has no runtime branching) and the slice bounds
45+ * come from the table's header, so one plan serves any wavelength. Slices that
46+ * fall past the end of a small table exit immediately, which is why a coarse
47+ * wavelength pays nothing for the split.
48+ */
49+const CHUNKS = 16;
50+
51+/** Floats the table needs for `nmodes` modes. */
52+export const modeTableLength = (nmodes: number): number =>
53+ HEADER + STRIDE * nmodes;
54+
55+/** Modes a table holds, from its header. */
56+export const modeCount = (table: Float32Array): number => table[0];
57+
58+/**
59+ * Largest table this will try to build, in f32. Not a policy about how fine a
60+ * wavelength is sensible — that is the caller's call, and a fine one is
61+ * merely slow — but the point past which the draw would fail anyway: the
62+ * host-side Float32Array alone would be 8 GB. The device's own
63+ * storage-buffer limit is checked separately, when the buffer is allocated.
64+ */
65+const MAX_TABLE_FLOATS = 2 ** 31;
66+
67+/** What the table starts at, before any seed has been drawn. Big enough for
68+ * the default wavelength on the shipped surfaces, so the common case never
69+ * reallocates. */
70+export const INITIAL_MODES = 4096;
71+
72+/** Wavelength of the seeded field when the app names none. Fine enough to
73+ * give a Turing pattern plenty to grow from, coarse enough that the draw is
74+ * ~1,400 modes rather than the ~11,500 of the slider's finest setting. */
75+export const DEFAULT_LAMBDA = 0.5;
76+
77+/** Floats before the first mode: `[nmodes, 0, 0, 0]`. The count travels in
78+ * the buffer rather than a second binding, so the kernel needs one storage
79+ * buffer and the host one write. */
80+const HEADER = 4;
81+/** Floats per mode: kx, ky, kz, real, imag. */
82+const STRIDE = 5;
83+
84+/** The name the coefficient buffer takes in the plan's HostBuffers. */
85+export const MODE_BUFFER = 'randnfun3_modes';
86+
87+/** How many dispatches `randnfun3WGSL` must be planned as. */
88+export const randnfun3Chunks = CHUNKS;
89+
90+/**
91+ * One thread per surface point, summing this chunk's slice of the modes.
92+ *
93+ * The inner loop is a dot product, a cos, a sin and two multiply-adds, over a
94+ * table small enough (~1,400 modes at the default lambda) to sit in cache for
95+ * every thread. Chunk 0 initializes the output and the rest accumulate onto
96+ * it; dispatches within one compute pass are ordered, so the reads see the
97+ * previous chunk's writes. Nothing here is per-step work: `init` runs once a
98+ * seed.
99+ */
100+export function randnfun3WGSL(npts: number, chunk: number): string {
101+ return `
102+@group(0) @binding(0) var<storage, read_write> outf: array<f32>;
103+@group(0) @binding(1) var<storage, read> px: array<f32>;
104+@group(0) @binding(2) var<storage, read> py: array<f32>;
105+@group(0) @binding(3) var<storage, read> pz: array<f32>;
106+// [nmodes, _, _, _], then kx, ky, kz, re, im per mode.
107+@group(0) @binding(4) var<storage, read> modes: array<f32>;
108+
109+@compute @workgroup_size(64)
110+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
111+ let i = gid.x;
112+ if (i >= ${npts}u) { return; }
113+ let n = u32(modes[0]);
114+ // This chunk's slice. Ceiling division, so the last slices are the short
115+ // ones and an empty slice costs a single comparison.
116+ let per = (n + ${CHUNKS}u - 1u) / ${CHUNKS}u;
117+ let lo = min(${chunk}u * per, n);
118+ let hi = min(lo + per, n);
119+ var acc = 0.0;
120+ if (lo < hi) {
121+ let x = px[i];
122+ let y = py[i];
123+ let z = pz[i];
124+ for (var m = lo; m < hi; m = m + 1u) {
125+ let b = ${HEADER}u + m * ${STRIDE}u;
126+ let t = modes[b] * x + modes[b + 1u] * y + modes[b + 2u] * z;
127+ acc = acc + modes[b + 3u] * cos(t) - modes[b + 4u] * sin(t);
128+ }
129+ }
130+${chunk === 0 ? ' outf[i] = acc;' : ' outf[i] = outf[i] + acc;'}
131+}
132+`;
133+}
134+
135+/**
136+ * Modes a wavelength will draw on a box, without drawing them: chebfun's
137+ * cube size, times the fraction its isotropy ball keeps (pi/6 of a cube,
138+ * approached from below at small m). Used to price a wavelength up front.
139+ */
140+function plannedModes(lambda: number, box: BoundingBox): number {
141+ const side = (w: number): number => 2 * Math.round((1.2 * w) / lambda + 2) + 1;
142+ const cube =
143+ side(box.x1 - box.x0) * side(box.y1 - box.y0) * side(box.z1 - box.z0);
144+ return Math.ceil((Math.PI / 6) * cube);
145+}
146+
147+/** The box a random field is drawn over: the surface's own bounding box. */
148+export interface BoundingBox {
149+ x0: number; x1: number;
150+ y0: number; y1: number;
151+ z0: number; z1: number;
152+}
153+
154+/** The bounding box of a surface, as `Geometry` holds its coordinates. */
155+export function boundingBox(
156+ x: Float32Array,
157+ y: Float32Array,
158+ z: Float32Array,
159+): BoundingBox {
160+ const box = {
161+ x0: Infinity, x1: -Infinity,
162+ y0: Infinity, y1: -Infinity,
163+ z0: Infinity, z1: -Infinity,
164+ };
165+ for (let i = 0; i < x.length; i++) {
166+ if (x[i] < box.x0) box.x0 = x[i];
167+ if (x[i] > box.x1) box.x1 = x[i];
168+ if (y[i] < box.y0) box.y0 = y[i];
169+ if (y[i] > box.y1) box.y1 = y[i];
170+ if (z[i] < box.z0) box.z0 = z[i];
171+ if (z[i] > box.z1) box.z1 = z[i];
172+ }
173+ return box;
174+}
175+
176+/**
177+ * Draw a field's modes and pack them for the GPU: `tools/randnfun3.m` run
178+ * through the interpreter, seeded, then interleaved into the buffer layout
179+ * above. Column-major out of MATLAB, interleaved on the way in.
180+ */
181+export function drawModes(
182+ lambda: number,
183+ box: BoundingBox,
184+ seed: number,
185+ /** Points the field will be summed at, for the cost budget. */
186+ npts: number,
187+): Float32Array {
188+ if (!(lambda > 0) || !Number.isFinite(lambda)) {
189+ throw new Error(`randnfun3: lambda must be a positive number, got ${lambda}`);
190+ }
191+ // The mode count follows from lambda and the box alone, so a table that
192+ // cannot be built is refused before anything is drawn. The only ceiling is
193+ // what fits: how slow a fine wavelength is, is the caller's to decide.
194+ const planned = plannedModes(lambda, box);
195+ if (modeTableLength(planned) > MAX_TABLE_FLOATS) {
196+ throw new Error(
197+ `randnfun3: lambda ${lambda} needs about ` +
198+ `${planned.toLocaleString()} Fourier modes on this surface, a ` +
199+ `${((4 * modeTableLength(planned)) / 1e9).toFixed(1)} GB table. ` +
200+ `lambda is an absolute length, so a larger surface needs more modes ` +
201+ `for the same value, and halving it costs eight times as many.`,
202+ );
203+ }
204+ const result = executeCode(
205+ 'rng(seed); [k, c] = randnfun3(lambda, [x0 x1 y0 y1 z0 z1]);',
206+ {
207+ initialVariableValues: { lambda, seed, ...box },
208+ displayResults: false,
209+ implicitCwdPath: null,
210+ },
211+ toolFiles,
212+ 'randnfun3-driver.m',
213+ );
214+ const k = result.variableValues['k'];
215+ const c = result.variableValues['c'];
216+ if (!k || !c || !isRuntimeTensor(k) || !isRuntimeTensor(c)) {
217+ throw new Error("randnfun3: tools/randnfun3.m did not return [k, c] arrays");
218+ }
219+ const nmodes = k.shape[0];
220+ const out = new Float32Array(modeTableLength(nmodes));
221+ out[0] = nmodes;
222+ for (let i = 0; i < nmodes; i++) {
223+ const b = HEADER + STRIDE * i;
224+ out[b] = k.data[i]; // kx
225+ out[b + 1] = k.data[nmodes + i]; // ky
226+ out[b + 2] = k.data[2 * nmodes + i]; // kz
227+ out[b + 3] = c.data[i]; // real
228+ out[b + 4] = c.data[nmodes + i]; // imag
229+ }
230+ return out;
231+}
232+
233+/**
234+ * `drawModes` on a worker thread, so a fine wavelength does not freeze the
235+ * page (src/mgpu/randnfun3.worker.ts).
236+ *
237+ * Falls back to drawing in place where there is no `Worker` — the node test
238+ * runner and the desktop benchmark, neither of which has an event loop it
239+ * would matter to. Failures surface as a rejection either way, so a caller
240+ * never has to know which path ran.
241+ */
242+export function drawModesAsync(
243+ lambda: number,
244+ box: BoundingBox,
245+ seed: number,
246+ npts: number,
247+): Promise<Float32Array> {
248+ if (typeof Worker === 'undefined') {
249+ try {
250+ return Promise.resolve(drawModes(lambda, box, seed, npts));
251+ } catch (e) {
252+ return Promise.reject(e instanceof Error ? e : new Error(String(e)));
253+ }
254+ }
255+ const w = drawWorker();
256+ const id = nextDrawId++;
257+ return new Promise((resolve, reject) => {
258+ pendingDraws.set(id, { resolve, reject });
259+ w.postMessage({ id, lambda, box, seed, npts });
260+ });
261+}
262+
263+let worker: Worker | null = null;
264+let nextDrawId = 1;
265+const pendingDraws = new Map<
266+ number,
267+ { resolve: (t: Float32Array) => void; reject: (e: Error) => void }
268+>();
269+
270+/** The draw worker, started on first use and kept for the session — starting
271+ * one re-parses numbl, which costs more than a coarse draw does. */
272+function drawWorker(): Worker {
273+ if (worker) return worker;
274+ worker = new Worker(new URL('./randnfun3.worker.ts', import.meta.url), {
275+ type: 'module',
276+ });
277+ worker.onmessage = (e: MessageEvent<{ id: number; table?: Float32Array; error?: string }>): void => {
278+ const waiting = pendingDraws.get(e.data.id);
279+ if (!waiting) return;
280+ pendingDraws.delete(e.data.id);
281+ if (e.data.error !== undefined) waiting.reject(new Error(e.data.error));
282+ else waiting.resolve(e.data.table!);
283+ };
284+ worker.onerror = (e: ErrorEvent): void => {
285+ // A worker that died takes every outstanding draw with it.
286+ for (const [, waiting] of pendingDraws) {
287+ waiting.reject(new Error(`randnfun3 draw worker failed: ${e.message}`));
288+ }
289+ pendingDraws.clear();
290+ worker?.terminate();
291+ worker = null;
292+ };
293+ return worker;
294+}
src/mgpu/randnfun3.worker.tsadded+40−0View file
@@ -0,0 +1,40 @@
1+/**
2+ * Drawing a seed field's Fourier modes, off the main thread.
3+ *
4+ * The draw is `tools/randnfun3.m` in numbl's interpreter, and its cost goes as
5+ * the inverse cube of the wavelength: milliseconds at the default lambda, but
6+ * ~13 s at lambda 0.01. Synchronous JS that long does not merely feel slow —
7+ * it blocks the event loop outright, so the page stops painting and the
8+ * browser offers to kill it. Nothing about the draw needs the main thread
9+ * (it touches no GPU and no DOM), so it runs here and the result is
10+ * transferred back.
11+ */
12+import { drawModes, type BoundingBox } from './randnfun3.ts';
13+
14+export interface DrawRequest {
15+ id: number;
16+ lambda: number;
17+ box: BoundingBox;
18+ seed: number;
19+ npts: number;
20+}
21+
22+export type DrawReply =
23+ | { id: number; table: Float32Array; error?: undefined }
24+ | { id: number; table?: undefined; error: string };
25+
26+self.onmessage = (e: MessageEvent<DrawRequest>): void => {
27+ const { id, lambda, box, seed, npts } = e.data;
28+ let reply: DrawReply;
29+ let transfer: Transferable[] = [];
30+ try {
31+ const table = drawModes(lambda, box, seed, npts);
32+ reply = { id, table };
33+ transfer = [table.buffer];
34+ } catch (err) {
35+ reply = { id, error: err instanceof Error ? err.message : String(err) };
36+ }
37+ (self as unknown as {
38+ postMessage: (m: DrawReply, t: Transferable[]) => void;
39+ }).postMessage(reply, transfer);
40+};
src/mgpu/registry.tsmodified+7−0View file
@@ -26,6 +26,13 @@ export interface ParamSpec {
2626 min: number;
2727 max: number;
2828 step: number;
29+ /**
30+ * This parameter is a random seed: its value picks a draw and means nothing
31+ * on its own, so the UI offers a button that jumps to another one rather
32+ * than a box to type a number into. `min`/`max` still bound what the button
33+ * picks.
34+ */
35+ reseed?: boolean;
2936 }
3037
3138 export interface MModel {
src/mgpu/session.tsmodified+58−8View file
@@ -11,6 +11,8 @@ import { DerivPlan } from '../sht/deriv.ts';
1111 import { gridForLmax, type ShtConfig } from '../sht/layout.ts';
1212 import { GpuModel, type ModelParams } from './model.ts';
1313 import { seededNoise } from './noise.ts';
14+import { boundingBox, drawModesAsync, DEFAULT_LAMBDA } from './randnfun3.ts';
15+import { resolveLambda } from './plan.ts';
1416 import type { MModel } from './registry.ts';
1517 import { Geometry } from '../geom/geometry.ts';
1618 import { mGeometryByKey, defaultGeometryParams, SPHERE_KEY, type MGeometry } from '../geom/registry.ts';
@@ -36,6 +38,9 @@ export interface ModelSessionOptions {
3638 * is unrolled into the op sequence, so a change recompiles.
3739 */
3840 niter?: number;
41+ /** Wavelength of the seeded random field a model's `init` draws
42+ * (src/mgpu/randnfun3.ts). Redrawn on the next seed, never recompiled. */
43+ lam3?: number;
3944 }
4045
4146 export class ModelSession {
@@ -62,6 +67,10 @@ export class ModelSession {
6267 /** Display-only transforms on the oversampled grid; null at 1x. */
6368 #displaySht: ShtPlan | null;
6469 #oversample: number;
70+ /** Wavelength of the seeded random field, and the seed it was drawn from —
71+ * kept so changing one can redraw with the other unchanged. */
72+ #lam3: number;
73+ #seed = 1;
6574
6675 private constructor(init: {
6776 device: GPUDevice;
@@ -76,6 +85,7 @@ export class ModelSession {
7685 geometryModel: MGeometry;
7786 deriv: DerivPlan;
7887 niter: number;
88+ lam3: number;
7989 }) {
8090 this.device = init.device;
8191 this.model = init.model;
@@ -90,6 +100,7 @@ export class ModelSession {
90100 this.#geometryModel = init.geometryModel;
91101 this.#deriv = init.deriv;
92102 this.niter = init.niter;
103+ this.#lam3 = init.lam3;
93104 }
94105
95106 get geometry(): Geometry {
@@ -138,7 +149,6 @@ export class ModelSession {
138149 // buffers of numbers (the embedding, and both metric formulations built
139150 // on it: the inverse metric quantities and the flux-form weights).
140151 const geometry = await Geometry.create({
141- device,
142152 sht,
143153 cfg,
144154 source: opts.geometrySource ?? geometryModel.source,
@@ -158,10 +168,11 @@ export class ModelSession {
158168 deriv,
159169 niter,
160170 });
161- gpu.setParams(params);
171+ const lam3 = opts.lam3 ?? DEFAULT_LAMBDA;
172+ gpu.setParams({ lam3, ...params });
162173 return new ModelSession({
163174 device, model, cfg, sht, displaySht, gpu, params, oversample,
164- geometry, geometryModel, deriv, niter,
175+ geometry, geometryModel, deriv, niter, lam3,
165176 });
166177 } catch (e) {
167178 // The transform plans own GPU buffers; do not leak them on a compile error.
@@ -194,7 +205,6 @@ export class ModelSession {
194205 source?: string,
195206 ): Promise<void> {
196207 const next = await Geometry.create({
197- device: this.device,
198208 sht: this.sht,
199209 cfg: this.cfg,
200210 source: source ?? geometryModel.source,
@@ -240,16 +250,56 @@ export class ModelSession {
240250 old?.destroy();
241251 }
242252
243- /** Run `init` from a seeded perturbation, resetting model time. */
244- seed(seed: number): void {
245- this.gpu.init(seededNoise(this.npts, this.model.seedAmp, seed));
253+ /**
254+ * Run `init` from a seeded perturbation, resetting model time.
255+ *
256+ * A model that calls `randnfun3` gets its coefficient table drawn here:
257+ * over the current surface's bounding box, at the wavelength its own .m
258+ * asked for. The draw is host-side MATLAB (a few ms); the evaluation at
259+ * every grid point is the GPU kernel inside `init`.
260+ */
261+ async seed(seed: number): Promise<void> {
262+ const lambda = this.gpu.randnfun3Lambda;
263+ // Drawn on a worker: at a fine wavelength this is seconds of interpreter
264+ // time, and it must not be seconds of frozen page.
265+ const modes = lambda
266+ ? await drawModesAsync(
267+ resolveLambda(lambda, this.#mergedParams()),
268+ boundingBox(this.#geometry.x, this.#geometry.y, this.#geometry.z),
269+ seed,
270+ this.npts,
271+ )
272+ : null;
273+ await this.gpu.init(seededNoise(this.npts, this.model.seedAmp, seed), modes);
274+ this.#seed = seed;
246275 this.t = 0;
247276 this.steps = 0;
248277 }
249278
279+ /** The model's parameters plus the ones the host owns. */
280+ #mergedParams(): ModelParams {
281+ return { lam3: this.#lam3, ...this.#params };
282+ }
283+
250284 setParams(params: ModelParams): void {
251285 this.#params = params;
252- this.gpu.setParams(params);
286+ this.gpu.setParams(this.#mergedParams());
287+ }
288+
289+ /** Wavelength of the seeded random field. Redraws on the next seed. */
290+ get lam3(): number {
291+ return this.#lam3;
292+ }
293+
294+ /**
295+ * Change the random field's wavelength. Nothing recompiles — lam3 is a
296+ * uniform and the coefficient table is host-drawn — but the field itself
297+ * only changes on the next `seed`, which is where it is drawn.
298+ */
299+ setLam3(lambda: number): void {
300+ if (lambda === this.#lam3) return;
301+ this.#lam3 = lambda;
302+ this.gpu.setParams(this.#mergedParams());
253303 }
254304
255305 /** Advance `n` steps. Synchronous: records and submits, nothing read back. */
src/raw.d.tsmodified+10−0View file
@@ -3,3 +3,13 @@ declare module '*?raw' {
33 const source: string;
44 export default source;
55 }
6+
7+/** Vite's `import.meta.glob`, used to load every .m in tools/ at once
8+ * (src/tools.ts). Only the eager + `?raw` form this project uses is
9+ * declared — it returns each match's text, keyed by path. */
10+interface ImportMeta {
11+ glob(
12+ pattern: string,
13+ options: { query: '?raw'; eager: true; import: 'default' },
14+ ): Record<string, string>;
15+}
src/sht/sht.tsmodified+13−1View file
@@ -776,8 +776,20 @@ export async function requestShtDevice(): Promise<GPUDevice> {
776776 // timestamp-query is only used by the profiling scripts, but it has to be
777777 // requested at device creation, and asking costs nothing when unused.
778778 if (adapter.features.has('timestamp-query')) features.push('timestamp-query');
779+ // The seed field's mode table is the one buffer whose size is not fixed by
780+ // the grid — it grows with how fine a wavelength is asked for
781+ // (src/mgpu/randnfun3.ts), and a browser's default 128 MB storage-buffer
782+ // limit is well below what the adapter will actually give. Ask for the
783+ // adapter's own maximum so the wavelength is limited by the hardware rather
784+ // than by a default.
785+ const maxStorage = adapter.limits.maxStorageBufferBindingSize;
786+ const maxBuffer = adapter.limits.maxBufferSize;
779787 return adapter.requestDevice({
780788 requiredFeatures: features,
781- requiredLimits: { maxComputeWorkgroupStorageSize: wgStorage },
789+ requiredLimits: {
790+ maxComputeWorkgroupStorageSize: wgStorage,
791+ maxStorageBufferBindingSize: maxStorage,
792+ maxBufferSize: maxBuffer,
793+ },
782794 });
783795 }
src/tools.tsadded+27−0View file
@@ -0,0 +1,27 @@
1+/**
2+ * The shared MATLAB utilities in `tools/`, as interpreter workspace files.
3+ *
4+ * A geometry or a seeding draw is evaluated by numbl's interpreter (see
5+ * src/geom/geometry.ts), which resolves a call like `randnfunsphere(...)`
6+ * against the workspace files it is handed. Everything in `tools/` is handed
7+ * to every such run, so any .m can call any tool by name — MATLAB's own path
8+ * semantics, where the file name is the function name.
9+ *
10+ * These are *not* available to the models: a model's step compiles to WGSL,
11+ * where none of this exists.
12+ */
13+const sources = import.meta.glob('../tools/*.m', {
14+ query: '?raw',
15+ eager: true,
16+ import: 'default',
17+}) as Record<string, string>;
18+
19+export interface ToolFile {
20+ name: string;
21+ source: string;
22+}
23+
24+/** Every tool, named as MATLAB wants it (`randnfunsphere.m`). */
25+export const toolFiles: ToolFile[] = Object.entries(sources)
26+ .map(([path, source]) => ({ name: path.slice(path.lastIndexOf('/') + 1), source }))
27+ .sort((a, b) => a.name.localeCompare(b.name));
test/analyticChecks.tsmodified+2−3View file
@@ -65,7 +65,6 @@ async function makeModel(
6565 const sht = await ShtPlan.create(device, cfg);
6666 const deriv = await DerivPlan.create(device, sht);
6767 const geometry = await Geometry.create({
68- device,
6968 sht,
7069 cfg,
7170 source: mGeometryByKey(SPHERE_KEY)!.source,
@@ -163,7 +162,7 @@ export async function analyticChecks(
163162
164163 // Uniform initial field: stays uniform, and diffusion cannot touch it.
165164 const field = new Float32Array(npts).fill(u0);
166- gpu.init(field);
165+ await gpu.init(field, null);
167166 const Ustart = await gpu.read('U');
168167 gpu.step(nsteps);
169168 const Uend = await gpu.read('U');
@@ -220,7 +219,7 @@ export async function analyticChecks(
220219
221220 // Seed the exact homogeneous fixed point by handing init a zero
222221 // perturbation, then add a small single-mode bump to u only.
223- gpu.init(new Float32Array(npts));
222+ await gpu.init(new Float32Array(npts), null);
224223 const l = 24;
225224 const m = 7;
226225 const idx = lmIndex(lmax, l, m);
test/fluxChecks.tsmodified+2−2View file
@@ -239,7 +239,7 @@ export async function fluxChecks(
239239 const sht = await ShtPlan.create(device, cfg);
240240 const deriv = await DerivPlan.create(device, sht);
241241 const geometry = await Geometry.create({
242- device, sht, cfg,
242+ sht, cfg,
243243 source: g.source,
244244 paramNames: g.params.map((p) => p.key),
245245 params: defaultGeometryParams(g),
@@ -359,7 +359,7 @@ export async function fluxChecks(
359359 ).length,
360360 );
361361 if (niter === 1) {
362- session.seed(1);
362+ await session.seed(1);
363363 session.step(STEPS);
364364 states.push(await session.read('U'));
365365 }
test/geometryChecks.tsmodified+173−6View file
@@ -37,6 +37,7 @@ import {
3737 SPHERE_KEY,
3838 } from '../src/geom/registry.ts';
3939 import { ModelCompileError } from '../src/mgpu/errors.ts';
40+import { boundingBox, drawModes, DEFAULT_LAMBDA } from '../src/mgpu/randnfun3.ts';
4041 import type { Check, Log } from './analyticChecks.ts';
4142
4243 const LMAX = 31;
@@ -56,7 +57,6 @@ async function buildGeometry(device: GPUDevice, key: string) {
5657 const sht = await ShtPlan.create(device, cfg);
5758 const deriv = await DerivPlan.create(device, sht);
5859 const geometry = await Geometry.create({
59- device,
6060 sht,
6161 cfg,
6262 source: g.source,
@@ -295,7 +295,7 @@ export async function geometryChecks(
295295 device, model, params, lmax: LMAX, niter,
296296 });
297297 ops.push(session.describe().step.length);
298- session.seed(1);
298+ await session.seed(1);
299299 session.step(STEPS);
300300 states.push(await session.read('U'));
301301 session.destroy();
@@ -361,7 +361,7 @@ export async function geometryChecks(
361361 device, model, params, lmax: SWEEP_LMAX,
362362 geometry: peanut, geometryParams: peanutParams, niter,
363363 });
364- session.seed(1);
364+ await session.seed(1);
365365 session.step(STEPS);
366366 states.push(await session.read('U'));
367367 session.destroy();
@@ -416,7 +416,7 @@ export async function geometryChecks(
416416 geometry: geomSpec, geometryParams: defaultGeometryParams(geomSpec),
417417 niter,
418418 });
419- session.seed(1);
419+ await session.seed(1);
420420 session.step(STEPS);
421421 const values = await session.read('u');
422422 const finite = values.every((v) => Number.isFinite(v));
@@ -469,7 +469,7 @@ export async function geometryChecks(
469469 `rate ${((g.muMax - g.muMin) / (g.muMax + g.muMin)).toFixed(3)} ` +
470470 `vs plain ${(g.muMax - 1).toFixed(2)}`;
471471 }
472- session.seed(1);
472+ await session.seed(1);
473473 session.step(STEPS);
474474 const values = await session.read('u');
475475 outcomes.push(values.every((v) => Number.isFinite(v)));
@@ -512,7 +512,7 @@ export async function geometryChecks(
512512 const session = await ModelSession.create({
513513 device, model, params: defaultParams(model), lmax: LMAX,
514514 });
515- session.seed(1);
515+ await session.seed(1);
516516 session.step(STEPS);
517517 const before = await session.read('U');
518518
@@ -534,6 +534,173 @@ export async function geometryChecks(
534534 );
535535 session.destroy();
536536 }
537+
538+ await randnfun3Checks(device, check, log);
539+}
540+
541+/**
542+ * The seeded initial condition: chebfun's randnfun3, drawn on the host and
543+ * summed on the GPU (src/mgpu/randnfun3.ts).
544+ *
545+ * The split is the thing worth testing. The draw is MATLAB whose distribution
546+ * is checked directly, and the sum is a WGSL kernel checked against the same
547+ * modes evaluated in f64 on the CPU — if the kernel's indexing into the packed
548+ * mode table were wrong it would still produce a smooth random-looking field,
549+ * which is exactly the kind of wrong no "looks patterned" check would catch.
550+ */
551+async function randnfun3Checks(
552+ device: GPUDevice,
553+ check: Check,
554+ log: Log,
555+): Promise<void> {
556+ const model = mModelByKey('schnakenberg')!;
557+ const params = defaultParams(model);
558+ const make = (lam3: number): Promise<ModelSession> =>
559+ ModelSession.create({ device, model, params, lmax: LMAX, lam3 });
560+
561+ // ---- the GPU sum matches the same modes evaluated on the CPU -----------
562+ {
563+ const session = await make(DEFAULT_LAMBDA);
564+ await session.seed(3);
565+ // `u` after init is the steady state plus 0.01*f, so the field is
566+ // recovered by removing the model's own uniform offset.
567+ const u = await session.read('u');
568+ const g = session.geometry;
569+ const modes = drawModes(
570+ DEFAULT_LAMBDA,
571+ boundingBox(g.x, g.y, g.z),
572+ 3,
573+ g.x.length,
574+ );
575+ const nmodes = modes[0];
576+
577+ // The same sum in f64, straight from the packed table the GPU read.
578+ let maxErr = 0;
579+ let amp = 0;
580+ const us = params.a + params.b;
581+ for (let i = 0; i < g.x.length; i++) {
582+ let f = 0;
583+ for (let j = 0; j < nmodes; j++) {
584+ const b = 4 + 5 * j;
585+ const t = modes[b] * g.x[i] + modes[b + 1] * g.y[i] + modes[b + 2] * g.z[i];
586+ f += modes[b + 3] * Math.cos(t) - modes[b + 4] * Math.sin(t);
587+ }
588+ const want = us + 0.01 * f;
589+ maxErr = Math.max(maxErr, Math.abs(u[i] - want));
590+ amp = Math.max(amp, Math.abs(0.01 * f));
591+ }
592+ log(` randnfun3: ${nmodes} modes at lambda ${DEFAULT_LAMBDA}, |perturbation| up to ${amp.toExponential(2)}`);
593+ check(
594+ 'randnfun3: the GPU sum matches the same modes summed on the CPU',
595+ // fp32 over ~1400 terms against f64, on a field of amplitude ~1e-2.
596+ maxErr < 2e-6 && amp > 1e-3,
597+ `max |GPU - CPU| = ${maxErr.toExponential(2)}, perturbation amplitude ${amp.toExponential(2)}`,
598+ );
599+ session.destroy();
600+ }
601+
602+ // ---- a seed reproduces, a different seed does not ----------------------
603+ {
604+ const a = await make(DEFAULT_LAMBDA);
605+ await a.seed(11);
606+ const first = await a.read('u');
607+ await a.seed(11);
608+ const again = await a.read('u');
609+ await a.seed(12);
610+ const other = await a.read('u');
611+ let same = true;
612+ let differs = false;
613+ for (let i = 0; i < first.length; i++) {
614+ if (first[i] !== again[i]) same = false;
615+ if (first[i] !== other[i]) differs = true;
616+ }
617+ check(
618+ 'randnfun3: the same seed redraws the same field, a different one does not',
619+ same && differs,
620+ same ? (differs ? 'reproducible and seed-dependent' : 'seed 12 gave seed 11 back') : 'not reproducible',
621+ );
622+ a.destroy();
623+ }
624+
625+ // ---- the field is smooth, and lambda sets how smooth -------------------
626+ //
627+ // This is what randnfun3 buys over the white noise it replaced: the seed is
628+ // band-limited, so it is fully resolved by the grid instead of being
629+ // whatever the grid happened to alias. Measured as the share of spectral
630+ // energy above degree 20 — near zero for a smooth field, and larger for a
631+ // shorter wavelength, which is the direction lambda is supposed to move it.
632+ {
633+ const tail = async (lam3: number): Promise<number> => {
634+ const session = await make(lam3);
635+ await session.seed(5);
636+ const U = await session.read('U');
637+ let lo = 0;
638+ let hi = 0;
639+ for (let m = 0; m <= LMAX; m++) {
640+ for (let l = m; l <= LMAX; l++) {
641+ const i = lmIndex(LMAX, l, m);
642+ const e = U[2 * i] ** 2 + U[2 * i + 1] ** 2;
643+ if (l > 20) hi += e;
644+ else lo += e;
645+ }
646+ }
647+ session.destroy();
648+ return hi / (lo + hi);
649+ };
650+ const coarse = await tail(1);
651+ const fine = await tail(0.4);
652+ log(` randnfun3: energy above l=20 is ${coarse.toExponential(2)} at lambda 1, ${fine.toExponential(2)} at lambda 0.4`);
653+ check(
654+ 'randnfun3: the seed is band-limited, and lambda sets its scale',
655+ coarse < 1e-3 && fine > coarse,
656+ `tail ${coarse.toExponential(2)} (lambda 1) < ${fine.toExponential(2)} (lambda 0.4)`,
657+ );
658+ }
659+
660+ // ---- a finer wavelength grows the table rather than being capped -------
661+ //
662+ // The mode table is sized to the wavelength asked for, so going finer
663+ // reallocates it and rebinds the dispatch. Getting that wrong would leave
664+ // the kernel reading a destroyed buffer or a stale one, so check that a
665+ // fine field is actually there and actually different.
666+ {
667+ const session = await make(DEFAULT_LAMBDA);
668+ await session.seed(21);
669+ const coarse = await session.read('u');
670+ session.setLam3(0.12);
671+ await session.seed(21);
672+ const fine = await session.read('u');
673+ let differs = false;
674+ let finite = true;
675+ for (let i = 0; i < fine.length; i++) {
676+ if (!Number.isFinite(fine[i])) finite = false;
677+ if (fine[i] !== coarse[i]) differs = true;
678+ }
679+ check(
680+ 'randnfun3: a finer wavelength grows the mode table and rebinds',
681+ finite && differs,
682+ finite ? 'redrew finer, buffer rebound' : 'field went non-finite after resize',
683+ );
684+ session.destroy();
685+ }
686+
687+ // ---- a wavelength past the cost budget is refused, not truncated -------
688+ {
689+ const session = await make(DEFAULT_LAMBDA);
690+ let message = '';
691+ try {
692+ session.setLam3(1e-4);
693+ await session.seed(1);
694+ } catch (e) {
695+ message = e instanceof Error ? e.message : String(e);
696+ }
697+ check(
698+ 'randnfun3: a wavelength whose table could not be built is refused',
699+ message.includes('Fourier modes on this surface'),
700+ message ? `refused: ${message.slice(0, 62)}…` : 'drew it anyway',
701+ );
702+ session.destroy();
703+ }
537704 }
538705
539706 /** Index of the entry minimizing `score`, over the first `n` entries. */
test/modelChecks.tsmodified+3−3View file
@@ -119,7 +119,7 @@ export async function modelChecks(
119119 `${kernels} kernels (expected ${expected})`,
120120 );
121121
122- session.seed(1);
122+ await session.seed(1);
123123 session.step(STEPS);
124124
125125 // Every rendered field must be finite and have developed some contrast.
@@ -168,7 +168,7 @@ export async function modelChecks(
168168 .describe()
169169 .step.filter((l) => l.includes('[batch lane')).length;
170170 }
171- session.seed(1);
171+ await session.seed(1);
172172 session.step(STEPS);
173173 states.push(await session.read('U'));
174174 session.destroy();
@@ -252,7 +252,7 @@ export async function modelChecks(
252252 lmax: LMAX,
253253 oversample: 2,
254254 });
255- session.seed(1);
255+ await session.seed(1);
256256 session.step(STEPS);
257257
258258 const fine = await session.readSpecies(0);
test/test-page.tsmodified+2−2View file
@@ -87,7 +87,7 @@ async function soak(steps: number, lmax: number): Promise<void> {
8787 geometryParams: defaultGeometryParams(geometry),
8888 niter: DEFAULT_NITER,
8989 });
90- session.seed(5);
90+ await session.seed(5);
9191 log(
9292 `soak: ${steps} steps at lmax ${lmax} ` +
9393 `(grid ${session.cfg.nlat}x${session.cfg.nphi}, ${geometry.key}, ` +
@@ -191,7 +191,7 @@ async function dumpState(q: URLSearchParams): Promise<void> {
191191 geometryParams: spec.geometryParams,
192192 niter: spec.niter,
193193 });
194- session.seed(spec.seed);
194+ await session.seed(spec.seed);
195195 session.step(spec.steps);
196196 await session.sync();
197197 const state = await session.read(model.state[0]);
tools/randnfun3.madded+66−0View file
@@ -0,0 +1,66 @@
1+% Smooth random function in 3D — chebfun's randnfun3, as the Fourier modes
2+% it is built from rather than as a chebfun3.
3+%
4+% [K, C] = randnfun3(LAMBDA, DOM) draws a random trig series on the box
5+% DOM = [x0 x1 y0 y1 z0 z1] with maximum frequency about 2*pi/LAMBDA in
6+% each direction and standard normal distribution N(0,1) at each point.
7+% K is nmodes x 3 (angular wavenumbers) and C is nmodes x 2 (real and
8+% imaginary parts), defining
9+%
10+% f(x,y,z) = sum_j C(j,1)*cos(K(j,:)*[x;y;z]) - C(j,2)*sin(K(j,:)*[x;y;z])
11+%
12+% Seed the draw with rng(...) before calling.
13+%
14+% chebfun returns a chebfun3 and evaluates it later; this project has no
15+% such object, and the sum above is what the GPU evaluates at the surface
16+% points (src/mgpu/randnfun3.ts). Splitting it here is also what keeps the
17+% draw in MATLAB: randn has no counterpart in the compiled WGSL dialect.
18+
19+function [k, c] = randnfun3(lambda, dom)
20+ % chebfun's nonperiodic path builds a periodic function on a domain about
21+ % 20% larger and restricts it. Restriction is free when evaluating at
22+ % points, so we keep the enlarged period and never form the smaller one.
23+ m = round(1.2*(dom(2)-dom(1))/lambda + 2);
24+ n = round(1.2*(dom(4)-dom(3))/lambda + 2);
25+ p = round(1.2*(dom(6)-dom(5))/lambda + 2);
26+ m2 = 2*m+1;
27+ n2 = 2*n+1;
28+ p2 = 2*p+1;
29+ N = m2*n2*p2;
30+
31+ % chebfun draws the whole cube (column-major) before masking; drawing in
32+ % that same order keeps a seed meaning the same thing here as there.
33+ cr = randn(N, 1);
34+ ci = randn(N, 1);
35+
36+ % The cube's integer wavenumbers, -m:m x -n:n x -p:p in column-major order.
37+ i = (0:N-1).';
38+ jx = mod(i, m2) - m;
39+ jy = mod(floor(i/m2), n2) - n;
40+ jz = floor(i/(m2*n2)) - p;
41+
42+ % Confine to a ball for isotropy.
43+ keep = ((jx/m).^2 + (jy/n).^2 + (jz/p).^2) <= 1;
44+ jx = jx(keep);
45+ jy = jy(keep);
46+ jz = jz(keep);
47+ cr = cr(keep);
48+ ci = ci(keep);
49+
50+ % Normalize so the variance is 1 at each point.
51+ s = 1/sqrt(numel(cr));
52+ cr = s*cr;
53+ ci = s*ci;
54+
55+ % Angular wavenumbers on the enlarged period, which is a whole number of
56+ % wavelengths on each side.
57+ kx = 2*pi*jx/(m*lambda);
58+ ky = 2*pi*jy/(n*lambda);
59+ kz = 2*pi*jz/(p*lambda);
60+
61+ % Fold the box's origin into the phase, so evaluating is a plain sum over
62+ % cos(k.x) and sin(k.x) with no offset left to carry.
63+ ph = -(kx*dom(1) + ky*dom(3) + kz*dom(5));
64+ k = [kx, ky, kz];
65+ c = [cr.*cos(ph) - ci.*sin(ph), cr.*sin(ph) + ci.*cos(ph)];
66+end
tools/randnfunsphere.madded+59−0View file
@@ -0,0 +1,59 @@
1+% Smooth random function on the unit sphere — chebfun's randnfunsphere,
2+% evaluated at the given (theta, phi) instead of returned as a spherefun.
3+%
4+% F = randnfunsphere(LAMBDA, THETA, PHI) is a combination of all spherical
5+% harmonics up to degree floor(2*pi/LAMBDA) with independent N(0,1)
6+% coefficients, normalized so the variance is 1 at each point.
7+%
8+% randnfunsphere(LAMBDA, THETA, PHI, 'monochromatic') uses only the
9+% harmonics of that one degree, so every component has the same wave
10+% number — chebfun's 'monochrome' option.
11+%
12+% Seed the draw with rng(...) before calling. This project has no chebfun
13+% objects: what would be a spherefun there is returned here as values on the
14+% grid the caller passes in.
15+
16+function f = randnfunsphere(lambda, theta, phi, type)
17+ if ( nargin < 4 )
18+ type = 'white';
19+ end
20+ % The unit sphere has circumference 2*pi, matching randnfun's deg = L/lambda.
21+ deg = floor(2*pi/lambda);
22+ if ( strncmpi(type, 'm', 1) )
23+ c = randn(2*deg+1, 1);
24+ c = sqrt(4*pi/numel(c)) * c; % normalize so the variance is 1
25+ f = sphHarmSumFixedDeg(theta, phi, deg, c);
26+ else
27+ c = randn((deg+1)^2, 1);
28+ c = sqrt(4*pi/numel(c)) * c; % normalize so the variance is 1
29+ f = sphHarmSum(theta, phi, deg, c);
30+ end
31+end
32+
33+% All spherical harmonics up to degree deg, with coefficients ordered by
34+% degree and order (0, -1,0,1, -2,-1,0,1,2, ...). Order +m carries
35+% cos(m*phi), order -m carries sin(m*phi).
36+function f = sphHarmSum(theta, phi, deg, c)
37+ f = 1/sqrt(4*pi) * c(1) * ones(size(theta));
38+ k = 1; % coefficients consumed so far
39+ for l = 1:deg
40+ cl = c(k+1 : k+2*l+1); % this degree's orders, -l..l
41+ k = k + 2*l + 1;
42+ f = f + sphHarmSumFixedDeg(theta, phi, l, cl);
43+ end
44+end
45+
46+% All spherical harmonics of the single degree l.
47+function f = sphHarmSumFixedDeg(theta, phi, l, c)
48+ m = (0:l).';
49+ a = (-1).^m ./ sqrt((1 + double(m==0)) * pi);
50+ costh = cos(theta(:)).'; % legendre wants cos(theta), in a row
51+ G = legendre(l, costh, 'norm'); % (l+1) x npts
52+ f = 0 * theta;
53+ for mm = 0:l
54+ f = f + a(mm+1) * c(l+1+mm) * (G(mm+1,:).' .* cos(mm*phi));
55+ if mm > 0
56+ f = f + a(mm+1) * c(l+1-mm) * (G(mm+1,:).' .* sin(mm*phi));
57+ end
58+ end
59+end
vite.config.tsmodified+6−1View file
@@ -1,4 +1,5 @@
11 import { defineConfig } from 'vite';
2+import { realpathSync } from 'node:fs';
23 import { resolve } from 'node:path';
34
45 // numbl is a local `file:` dependency, so node_modules/numbl is a symlink to
@@ -8,7 +9,11 @@ import { resolve } from 'node:path';
89 // express this — Node rejects node_modules targets — and plain Node could not
910 // resolve numbl's internal `.js`->`.ts` imports anyway, which is why the GPU
1011 // tests run in the browser harness rather than under `node`.)
11-const numblSrc = resolve(import.meta.dirname, 'node_modules/numbl/src');
12+// Realpath'd through the symlink: dev serves modules under their real ids, so
13+// aliasing the node_modules path would give the same file two identities (one
14+// per spelling) and run its side effects twice — the interpreter's builtin
15+// registry throws on the second.
16+const numblSrc = realpathSync(resolve(import.meta.dirname, 'node_modules/numbl/src'));
1217
1318 export default defineConfig({
1419 base: './',