Initial commit: 2D acoustic scattering, live on WebGPU
MATLAB scene and step files, compiled through numbl to WebGPU compute
kernels. Time-domain leapfrog solver with a microphone, real-time
audio playback, and physical (SI) units throughout.
45 changed files+10235−0
.gitignoreadded+3−0View file
@@ -0,0 +1,3 @@
1+node_modules/
2+dist/
3+tmp/
README.mdadded+269−0View file
@@ -0,0 +1,269 @@
1+# acoustic-scattering-2d
2+
3+Sound scattering off obstacles in two dimensions, solved live in the browser.
4+The wave equation and the medium are both written as MATLAB, run by
5+[numbl](https://numbl.org) and compiled to WebGPU compute shaders; the pressure
6+field is drawn straight out of the buffer the solver writes, with no readback in
7+the display path.
8+
9+This is the flat sibling of
10+[turing-surface](https://github.com/concept-collection/turing-surface), which
11+solves reaction-diffusion on closed surfaces with spherical harmonics. The
12+compiler pipeline is the same idea, cut down: there is one host-provided
13+operation instead of six transforms, and the grid is a square.
14+
15+## Two files
16+
17+A **scene** says what the medium is. It is evaluated once, on the CPU, through
18+numbl's interpreter, so it has the whole language available: loops, `if`,
19+indexing, seeded randomness, anything in `tools/`.
20+
21+```matlab
22+function [c, sig] = medium(x, y, L, h, cin, R, absorb)
23+ r = sqrt(x.^2 + y.^2);
24+ inside = 0.5 * (1 - tanh((r - R) / (1.5*h)));
25+
26+ c = 1 + (cin - 1) * inside;
27+ sig = sponge(x, y, L, 0.2*L, 25) + absorb * inside;
28+end
29+```
30+
31+`c` is the sound speed and `sig` the absorption rate. Both are ordinary fields
32+of position, which is what lets one function describe both the scatterer and the
33+open boundary: the absorbing layer around the edge is the statement that the
34+medium swallows sound out there, and a scene is free to put absorption inside
35+the domain too, making a lossy scatterer.
36+
37+A **model** says how the field advances. It is compiled, so it is held to an
38+element-wise subset, plus the Laplacian stencils the host supplies.
39+
40+```matlab
41+function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
42+ ...
43+ sd = sig * dt;
44+ lap = lap2(p);
45+ pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);
46+
47+ pold = p;
48+ tn = t + dt;
49+end
50+```
51+
52+Neither file declares its own parameters. It names the arguments it wants, and
53+the host matches each against what it offers (`src/mgpu/registry.ts` and
54+`src/scene/registry.ts`); a name nobody provides is a compile error with a
55+position in the file, not a silently undefined variable.
56+
57+## The equation
58+
59+Pressure only, at constant density:
60+
61+```
62+p_tt + 2*sig*p_t = c(x,y)^2 * lap(p) + s(x, y, t)
63+```
64+
65+Centring both the second time derivative and the damping on step *n* gives the
66+explicit update above. Since the density is constant, the impedance ratio across
67+an interface is just the speed ratio, so a scatterer much faster than its
68+background behaves nearly rigid (sound-hard) and one much slower nearly
69+pressure-release (sound-soft). What this formulation cannot do is set impedance
70+and speed independently, which needs a variable-density divergence form and a
71+second field.
72+
73+Time is carried as a grid field rather than a number. A batch of timesteps is
74+one replay of a fixed sequence of GPU operations, so nothing the host writes
75+between frames can change inside it; a clock uploaded per frame would stand
76+still for the whole batch. `tn = t + dt` makes the model keep its own clock, at
77+the cost of one buffer and one kernel per step.
78+
79+## How the MATLAB becomes a shader
80+
81+1. **Lowering.** `parseMFile` and `specializeUserFunction` (numbl's own JIT
82+ front end) lower `init` and `step` for the concrete argument types of the
83+ current grid, so every type and shape is fixed at compile time and a step
84+ needs no dynamic dispatch at all.
85+2. **Fusion.** numbl's inline pass folds single-use temporaries into their
86+ consumers. It declines the ones its C back end cannot fuse (`sin`, `exp`,
87+ `tanh`), which WGSL is perfectly happy with, so a second pass
88+ (`src/mgpu/fuse.ts`, adapted from
89+ [math-webgpu-sandbox](https://github.com/concept-collection/math-webgpu-sandbox))
90+ folds those too. The whole source term collapses into the update.
91+3. **Planning.** Each remaining statement becomes one GPU operation: an
92+ element-wise kernel, or a stencil dispatch for `lap2`/`lap4`. The result is a
93+ static op list, so running a timestep is pure command recording, and a whole
94+ frame's worth of steps goes into one submit.
95+
96+On this machine the compiled step is five kernels, one stencil dispatch, and
97+three buffer copies:
98+
99+```
100+kernel u = <element-wise>
101+kernel sd = <element-wise>
102+stencil lap = lap2(p)
103+kernel pn = <element-wise>
104+kernel pold = <element-wise>
105+kernel tn = <element-wise>
106+copy pn -> p
107+copy pold -> pm
108+copy tn -> t
109+```
110+
111+The op list is shown next to the editor in the app, so what a line of MATLAB
112+costs is visible while it is being written.
113+
114+### The one thing that is not element-wise
115+
116+`lap2` and `lap4` are the only places where a grid point reads its neighbours.
117+Keeping them as named operations rather than as array slicing means a model
118+never has to know how the grid is laid out in the buffer, and the host is free
119+to implement each as a single dispatch. `lap2` is the 5-point stencil; `lap4` is
120+the fourth-order 9-point one, which costs twice the reads and buys much less
121+grid dispersion, so a pulse can travel many wavelengths without visibly coming
122+apart.
123+
124+### Kernels that will not fit
125+
126+A fused kernel binds one storage buffer per distinct field its line reads, plus
127+its output and the parameter block, and WebGPU guarantees only eight per compute
128+stage (Chrome gives exactly that on some machines; compatibility mode gives
129+four). The leapfrog update reads nine fields. numbl's inline pass does not know
130+about that limit, and a model has no way to ask it for less, since a temporary
131+used once is exactly what it folds away.
132+
133+So the planner enforces the budget itself: any child subtree that reads more
134+than one field is evaluated into its own buffer and replaced by a reference to
135+it, which leaves the parent reading at most one field per child. The arithmetic
136+is unchanged, in a few more passes over memory, and it only happens on a line
137+that would not otherwise compile. The test suite runs the whole model with the
138+budget squeezed to two fields and checks the answer against the unsplit one.
139+
140+## What it is honest about
141+
142+- **Single precision.** WebGPU has no f64. The stencil differences lose a few
143+ digits to cancellation, which is visible in the accuracy numbers below but far
144+ below the discretization error at any resolution the app runs interactively.
145+- **The absorbing layer is a sponge, not a PML.** An absorbing layer is itself an
146+ impedance mismatch, so it reflects; spreading it over a couple of wavelengths
147+ keeps that small but not zero. Measured, for the shipped profile: about 0.8% of
148+ the incident peak comes back into the interior. A perfectly matched layer would
149+ do much better, particularly at grazing incidence, at the cost of extra fields
150+ and a more involved update.
151+- **The absorbing layer costs domain.** It occupies 20% of the width at each
152+ edge, so the clear region is the middle 60%. Scenes and sources are placed
153+ with that in mind.
154+- **Grid dispersion is real.** A discrete Laplacian gets the wave speed slightly
155+ wrong at wavelengths the grid barely resolves, and the error accumulates with
156+ distance travelled. Raise the frequency in the app until `leapfrog` visibly
157+ trails, then switch to `leapfrog4` and watch it clean up.
158+- **No exact-solution comparison yet.** Scattering by a circular cylinder has a
159+ classical Bessel-Hankel series solution, and comparing against it would put a
160+ number on the total error rather than on the pieces. That is the obvious next
161+ thing to do.
162+
163+## Listening to it
164+
165+There is a microphone: the pressure at one grid point, sampled every timestep.
166+It is written on the GPU by a one-thread dispatch that runs after each step and
167+appends to a buffer, because the obvious implementation, reading the field back
168+and picking out one number, costs a GPU-to-CPU round trip per step and would be
169+slower than the step. The whole trace comes back once, when there is something
170+to play.
171+
172+Turning that into sound needs one choice, because the simulation has no seconds
173+in it: its clock is model time, in which the background speed is 1 and the
174+source frequency `f` is in cycles per model time unit. The choice is made as a
175+pitch. Asking for the source to sound at 440 Hz fixes the playback sample rate
176+at `pitch / (f * dt)`, since the trace holds one sample per timestep.
177+
178+The consequence is worth knowing before it surprises you: **a pulse is short.**
179+The default source is a few dozen cycles, and a few dozen cycles at a musical
180+pitch is a few tens of milliseconds however long you leave the simulation
181+running. What you hear is a click, correctly shaped by whatever the wave did on
182+its way to the microphone, but a click. For sustained sound, turn `continuous`
183+up so the source keeps going, and run long enough to fill some seconds of audio:
184+the speed control goes to 64 timesteps a frame for exactly that. The trace holds
185+262,144 samples, around twelve seconds at a typical playback rate.
186+
187+The recording is normalized before playback, which discards absolute amplitude:
188+that is what the colour scale is for. It restarts whenever the run does, and
189+whenever the timestep changes, since a trace is one sample per step and two
190+timesteps would be two sample rates in one buffer.
191+
192+## Scenes
193+
194+| scene | what it shows |
195+| --- | --- |
196+| Disk | One circular scatterer, the reference case. Slow, fast, or absorbing. |
197+| Room | A square enclosure with a doorway. Close the doorway and it is sealed. |
198+| Two slits | A hard screen with two apertures: diffraction and interference. |
199+| Lens | A smooth slow patch that refracts a plane wave to a focus. |
200+| Random medium | Weak random structure everywhere: multiple scattering, and a coda. |
201+
202+The room is the scene the microphone is for: put a click inside it and what
203+comes back is a direct arrival, then echoes, then reverberation. Its walls are
204+*slower* than the room rather than faster, which is worth explaining. Reflection
205+at an interface goes as `|c2 - c1|/(c2 + c1)` at constant density, so a wall at
206+c = 0.15 reflects about 75% of the pressure amplitude, as a wall at c = 7 would.
207+But the timestep is set by the fastest speed anywhere on the grid, so a rigid
208+wall makes every step of the whole simulation several times smaller while a slow
209+one is free. What a slow wall costs instead is resolution inside itself, which
210+is what the wall's absorption is for: it swallows what gets in, so the badly
211+resolved part never comes back out. Push `wall speed` above 1 for a genuinely
212+rigid room and watch the timestep readout drop to match.
213+
214+A scene may also ask for source settings — the room asks for a point source
215+inside itself, since a plane wave arriving from off-grid would have nothing to
216+do. Those land in the sliders like any other value, and can be moved afterwards.
217+
218+The source blends continuously between a line source (whose far field is a plane
219+wave) and a point source, and between a single Gaussian pulse and a wave that
220+turns on and stays on. Both are parameters rather than modes, because both are
221+one expression in the .m.
222+
223+## Tests
224+
225+```
226+npm run test:node # the solver, on desktop WebGPU (Google Dawn)
227+npm run smoke # the page itself, in headless Chrome
228+```
229+
230+The microphone is checked exactly rather than approximately: after N steps its
231+trace must be N samples long, and its last sample must equal the pressure
232+sitting at the probe's grid point, both being the same f32 written by the same
233+kernel. That pins the probe index, the grid layout, and the fact that the
234+recording dispatch runs once per timestep rather than once per submission.
235+
236+The solver checks are physical wherever they can be. A pulse must travel at the
237+speed the medium says (measured wavefront radius, within 1%); a disk whose speed
238+matches its background must not scatter at all (bit-identical to the uniform
239+medium); a hard disk must scatter strongly; the absorbing layer must leave under
240+2% behind; both schemes must be stable at 95% of the CFL limit they claim. The
241+stencils are compared against a field whose Laplacian is known exactly, which at
242+n = 64 gives 1.1e-2 relative error for `lap2` and 1.9e-4 for `lap4`.
243+
244+The browser check drives the real page: it loads it, waits for the solver to
245+take steps and the frame loop to turn, and swaps the scene. It says nothing
246+about whether the picture is right, which is what eyes are for.
247+
248+`node scripts/smoke.mjs --recompile` goes further and swaps the model, breaks
249+the source deliberately (requiring the failure to be reported rather than
250+thrown), and reverts. That is opt-in because it needs new GPU pipelines built
251+while the page is already drawing, which headless Chrome cannot always do: on
252+some machines every GPU object created after the canvas context has been used
253+fails with "A valid external Instance reference no longer exists", in a
254+thirty-line WebGPU page with none of this project in it. A failure there says
255+as much about the browser as about the app.
256+
257+## Development
258+
259+```
260+npm install # numbl must be checked out as a sibling: ../../numbl
261+npm run dev
262+```
263+
264+The compiler runs from numbl's TypeScript sources directly (through a
265+`numbl-src` vite alias), so no numbl build is needed.
266+
267+## License
268+
269+Apache-2.0
index.htmladded+266−0View file
@@ -0,0 +1,266 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="utf-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1" />
6+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><rect width=%22100%22 height=%22100%22 fill=%22%23f2f2f2%22/><circle cx=%2258%22 cy=%2250%22 r=%2216%22 fill=%22%23888%22/><path d=%22M8 20v60M22 14v72M36 20v60%22 stroke=%22%233b4cc0%22 stroke-width=%228%22 fill=%22none%22/><path d=%22M76 30a24 24 0 0 1 0 40%22 stroke=%22%23b40426%22 stroke-width=%227%22 fill=%22none%22/></svg>" />
7+ <title>acoustic-scattering-2d — 2D scattering, live in the browser</title>
8+ <style>
9+ :root {
10+ --bg: #ffffff;
11+ --ink: #1f2328;
12+ --ink-2: #57606a;
13+ --line: #d0d7de;
14+ --accent: #0969da;
15+ --panel-bg: #f4f6f8;
16+ --tok-com: #6e7781;
17+ --tok-str: #0a3069;
18+ --tok-num: #0550ae;
19+ --tok-kw: #cf222e;
20+ --tok-ext: #8250df;
21+ color-scheme: light dark;
22+ }
23+ @media (prefers-color-scheme: dark) {
24+ :root {
25+ --bg: #14171a;
26+ --ink: #e6e9ec;
27+ --ink-2: #9aa4af;
28+ --line: #333b44;
29+ --accent: #58a6ff;
30+ --panel-bg: #14161c;
31+ --tok-com: #8b949e;
32+ --tok-str: #a5d6ff;
33+ --tok-num: #79c0ff;
34+ --tok-kw: #ff7b72;
35+ --tok-ext: #d2a8ff;
36+ }
37+ }
38+ body {
39+ margin: 0;
40+ background: var(--bg);
41+ color: var(--ink);
42+ font: 15px/1.5 system-ui, -apple-system, sans-serif;
43+ }
44+ main { max-width: 1100px; margin: 0 auto; padding: 20px 16px 48px; }
45+ h1 { font-size: 20px; margin: 0 0 2px; }
46+ .sub { color: var(--ink-2); margin: 0 0 12px; font-size: 13px; }
47+ .sub a { color: var(--accent); }
48+ .controls {
49+ display: flex; flex-wrap: wrap; gap: 8px 14px; align-items: center;
50+ padding: 5px 0;
51+ }
52+ .controls label { color: var(--ink-2); font-size: 13px; white-space: nowrap; }
53+ select, input[type="number"], button {
54+ font: inherit; font-size: 13px;
55+ color: var(--ink); background: var(--bg);
56+ border: 1px solid var(--line); border-radius: 6px;
57+ padding: 4px 8px;
58+ }
59+ input[type="range"] { width: 8em; vertical-align: middle; accent-color: var(--accent); }
60+ button { cursor: pointer; }
61+ button:hover { border-color: var(--accent); }
62+ button.primary { border-color: var(--accent); color: var(--accent); font-weight: 600; min-width: 5.5em; }
63+ .sliders {
64+ display: grid; grid-template-columns: repeat(auto-fill, minmax(230px, 1fr));
65+ gap: 2px 16px; padding: 4px 0;
66+ }
67+ #micbar { gap: 8px 16px; }
68+ #micparams { display: flex; gap: 16px; padding: 0; }
69+ #recinfo { margin-top: 0; }
70+ .slider { display: flex; align-items: center; gap: 8px; font-size: 13px; color: var(--ink-2); }
71+ .slider > span:first-child { flex: 0 0 8.5em; text-align: right; }
72+ .slider > output { flex: 0 0 4em; font-variant-numeric: tabular-nums; color: var(--ink); }
73+ .group-title {
74+ font-size: 12px; text-transform: uppercase; letter-spacing: 0.05em;
75+ color: var(--ink-2); margin: 10px 0 0;
76+ }
77+ #stage {
78+ display: flex; gap: 14px; margin-top: 12px; align-items: stretch;
79+ border: 1px solid var(--line); border-radius: 8px; overflow: hidden;
80+ }
81+ .canvas-box { flex: 1; aspect-ratio: 1 / 1; max-height: 70vh; position: relative; background: var(--panel-bg); }
82+ #view { width: 100%; height: 100%; display: block; cursor: crosshair; touch-action: none; }
83+ .colorbar {
84+ display: flex; flex-direction: column; align-items: center; justify-content: center;
85+ gap: 4px; padding: 8px 4px; background: var(--panel-bg);
86+ width: 60px; flex: none; box-sizing: border-box;
87+ }
88+ .colorbar canvas { border: 1px solid var(--line); border-radius: 2px; }
89+ .colorbar-label { font-size: 11px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
90+ .readout { font-variant-numeric: tabular-nums; color: var(--ink); }
91+ .readout.unstable, .warn { color: #b35900; font-weight: 600; }
92+ .stats { margin-top: 8px; font-size: 13px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
93+ .stats b { color: var(--ink); font-weight: 600; }
94+ #blurb { margin-top: 6px; font-size: 13px; color: var(--ink-2); }
95+ #err {
96+ color: #b35900; white-space: pre-wrap; font-size: 13px;
97+ font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
98+ }
99+ .editor {
100+ margin-top: 12px; border: 1px solid var(--line); border-radius: 8px;
101+ overflow: hidden;
102+ }
103+ .editor-head {
104+ display: flex; gap: 10px; align-items: center; justify-content: space-between;
105+ padding: 6px 10px; font-size: 12px; color: var(--ink-2);
106+ background: var(--panel-bg); border-bottom: 1px solid var(--line);
107+ }
108+ .editor-head button { padding: 2px 10px; font-size: 12px; }
109+ /* The fixed height lives on the row, not on either child: sizing the row
110+ makes both children stretch to one shared pixel height regardless of
111+ either's font size. */
112+ .editor-body { display: flex; align-items: stretch; height: 32em; }
113+ .editor-code { position: relative; flex: 1 1 62%; min-width: 0; }
114+ /* The overlay and the textarea must agree on every metric that affects
115+ where a character lands. Keep these two rules together. */
116+ .editor-code > pre,
117+ .editor-code > textarea {
118+ margin: 0; padding: 10px 12px; border: 0;
119+ box-sizing: border-box; width: 100%; height: 100%;
120+ font: 12.5px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
121+ tab-size: 2;
122+ white-space: pre; overflow-wrap: normal;
123+ }
124+ #highlight {
125+ position: absolute; inset: 0; overflow: hidden;
126+ pointer-events: none; background: var(--bg); color: var(--ink);
127+ }
128+ #source {
129+ position: relative; z-index: 1; display: block;
130+ resize: none; overflow: auto;
131+ background: transparent; color: transparent; caret-color: var(--ink);
132+ }
133+ #source:focus { outline: none; }
134+ /* Transparent text means the selection must be see-through, or selected
135+ code would be invisible. */
136+ #source::selection { background: color-mix(in srgb, var(--accent) 28%, transparent); }
137+ #compiled {
138+ flex: 1 1 38%; min-width: 0; margin: 0; padding: 10px 12px; overflow: auto;
139+ border-left: 1px solid var(--line); background: var(--panel-bg);
140+ font: 12px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
141+ color: var(--ink-2); white-space: pre;
142+ }
143+ @media (max-width: 860px) {
144+ .editor-body { flex-direction: column; height: auto; }
145+ .editor-code { flex: none; height: 26em; }
146+ #compiled { border-left: 0; border-top: 1px solid var(--line); max-height: 12em; }
147+ }
148+ .tok-com { color: var(--tok-com); }
149+ .tok-str { color: var(--tok-str); }
150+ .tok-num { color: var(--tok-num); }
151+ .tok-kw { color: var(--tok-kw); font-weight: 600; }
152+ .tok-ext { color: var(--tok-ext); }
153+ </style>
154+ </head>
155+ <body>
156+ <main>
157+ <h1>acoustic-scattering-2d</h1>
158+ <p class="sub">
159+ Sound scattering off obstacles in two dimensions, solved live on your
160+ GPU. The wave equation and the medium are both the MATLAB below, run in
161+ your browser by <a href="https://numbl.org">numbl</a> and compiled to
162+ WebGPU compute shaders. Edit either and watch it change.
163+ </p>
164+ <p class="sub" id="domaininfo"></p>
165+
166+ <div class="controls">
167+ <label>scene
168+ <select id="scene"></select>
169+ </label>
170+ <label title="The time-stepping scheme. Changing it recompiles.">model
171+ <select id="model"></select>
172+ </label>
173+ <label title="Grid points per side. Changing it recompiles and restarts.">grid
174+ <select id="gridsize">
175+ <option value="128">128²</option>
176+ <option value="256">256²</option>
177+ <option value="384">384²</option>
178+ <option value="512" selected>512²</option>
179+ </select>
180+ </label>
181+ <button id="runpause" class="primary">Run</button>
182+ <button id="restart" title="Back to a silent grid at t = 0">Restart</button>
183+ </div>
184+
185+ <p class="group-title">source</p>
186+ <div class="sliders" id="params"></div>
187+ <p class="group-title" id="scene-title">scene</p>
188+ <div class="sliders" id="sceneparams"></div>
189+
190+ <div class="controls">
191+ <label>colormap
192+ <select id="colormap"></select>
193+ </label>
194+ <label title="Pressure the colour scale saturates at. Auto follows the field.">scale
195+ <select id="scalemode">
196+ <option value="auto" selected>auto</option>
197+ <option value="fixed">hold</option>
198+ </select>
199+ </label>
200+ <label title="Shade the sound-speed field underneath the wave">
201+ <input type="checkbox" id="showmedium" checked /> show medium
202+ </label>
203+ <label id="cfl-label" title="The timestep, as a fraction of the largest one this scheme is stable at on this grid. At 1 and above it diverges, which is worth seeing once.">timestep
204+ <input type="range" id="cfl" min="0.05" max="1.2" step="0.05" value="0.5" />
205+ <output id="dtout" class="readout"></output>
206+ </label>
207+ <label title="Timesteps taken between frames">speed
208+ <select id="spf">
209+ <option value="1">1×</option>
210+ <option value="2">2×</option>
211+ <option value="4" selected>4×</option>
212+ <option value="8">8×</option>
213+ <option value="16">16×</option>
214+ <option value="32">32×</option>
215+ <option value="64">64×</option>
216+ </select>
217+ </label>
218+ </div>
219+
220+ <p class="group-title">microphone — click or drag on the picture to move it</p>
221+ <div class="controls" id="micbar">
222+ <div class="sliders" id="micparams"></div>
223+ <button id="listen" title="Play what the microphone has recorded since the run started, in real time at its real pitch">Listen</button>
224+ <span class="stats" id="recinfo"></span>
225+ </div>
226+
227+ <div id="stage">
228+ <div class="canvas-box"><canvas id="view"></canvas></div>
229+ <div id="colorbar"></div>
230+ </div>
231+ <p class="stats" id="stats"></p>
232+ <p id="blurb"></p>
233+ <p id="err"></p>
234+
235+ <div class="editor">
236+ <div class="editor-head">
237+ <span>
238+ <select id="editor-file" aria-label="file to edit">
239+ <option value="model">model (.m) — the solver</option>
240+ <option value="scene">scene (.m) — the medium</option>
241+ </select>
242+ <span id="editor-title"></span>
243+ </span>
244+ <span>
245+ <button id="recompile" type="button">Run edits</button>
246+ <button id="revert" type="button">Revert</button>
247+ </span>
248+ </div>
249+ <div class="editor-body">
250+ <div class="editor-code">
251+ <pre id="highlight" aria-hidden="true"></pre>
252+ <textarea
253+ id="source"
254+ spellcheck="false"
255+ autocomplete="off"
256+ autocapitalize="off"
257+ aria-label="source (MATLAB)"
258+ ></textarea>
259+ </div>
260+ <pre id="compiled"></pre>
261+ </div>
262+ </div>
263+ </main>
264+ <script type="module" src="/src/main.ts"></script>
265+ </body>
266+</html>
models/leapfrog.madded+64−0View file
@@ -0,0 +1,64 @@
1+% Second-order leapfrog for the 2D acoustic wave equation.
2+%
3+% p_tt + 2*sig*p_t = c^2 * lap(p) + s(x, y, t)
4+%
5+% Pressure only, at constant density, so the medium is the sound speed c
6+% (metres per second) and the absorption sig (inverse seconds) that the scene
7+% defines. Everything here is SI: x, y in metres, t and dt in seconds, f in
8+% hertz — none of it is declared as such anywhere, it simply follows from x,
9+% y and t being metres and seconds, which is the scene's and the app's doing,
10+% not this file's. Centring both the second time
11+% derivative and the damping term on step n,
12+%
13+% (p+ - 2p + p-)/dt^2 + sig*(p+ - p-)/dt = c^2*lap(p) + s
14+%
15+% and solving for p+ gives the update below. It is explicit: the only thing
16+% that couples neighbouring points is `lap2`, the 5-point Laplacian, which the
17+% host supplies as a single GPU dispatch. Everything else here is element-wise
18+% and compiles to one kernel per line.
19+%
20+% Stability wants c*dt/h <= 1/sqrt(2); the app sets dt from the fastest speed
21+% anywhere in the scene, so a fast scatterer slows the whole run down.
22+%
23+% Time is carried as a field rather than a number, because a batch of steps is
24+% one replay of a fixed sequence of GPU operations and nothing the host writes
25+% between frames can change inside it. `tn = t + dt` makes the model keep its
26+% own clock, which is what lets the source term be correct however many steps
27+% are batched into a submit.
28+
29+function [p, pm, t] = init(npts)
30+ p = zeros(npts, 1);
31+ pm = zeros(npts, 1);
32+ t = zeros(npts, 1);
33+end
34+
35+function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
36+ % The source. `cw` blends between a Gaussian pulse (0) and a wave that
37+ % turns on smoothly and stays on (1); `point` blends between a line source
38+ % spanning the grid in y, whose far field is a plane wave, and a point
39+ % source at (x0, y0).
40+ %
41+ % The om^2 is only a choice of units, not a physical amplitude — this is a
42+ % body force of arbitrary strength, not a source with a real acoustic power
43+ % rating. Such a force drives a response that falls off as 1/om^2 — two
44+ % time integrations — so without the factor the field would shrink tenfold
45+ % every time the frequency slider tripled. The equation is linear; scaling
46+ % the source scales the answer and nothing else.
47+ u = (t - t0) / tw;
48+ env = (1 - cw) * exp(-u .* u) + cw * (0.5 * (1 + tanh(u)));
49+ gx = ((x - x0) / w) .^ 2;
50+ gy = point * (((y - y0) / w) .^ 2);
51+ om = 2*pi*f;
52+ s = (om*om) * (env .* sin(om*(t - t0)) .* exp(-(gx + gy)));
53+
54+ % One step. The damping is what the absorbing layer acts through: sig is
55+ % zero over the interior, so there p+ is the plain leapfrog update.
56+ sd = sig * dt;
57+ lap = lap2(p);
58+ pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);
59+
60+ % This step's field becomes the next step's history. A line on its own, so
61+ % it plans as a copy rather than being folded into the update above.
62+ pold = p;
63+ tn = t + dt;
64+end
models/leapfrog4.madded+39−0View file
@@ -0,0 +1,39 @@
1+% Leapfrog with the fourth-order Laplacian, otherwise identical to
2+% models/leapfrog.m. The only change is `lap4` in place of `lap2`.
3+%
4+% What it buys is grid dispersion. A discrete Laplacian gets the wave speed
5+% slightly wrong at wavelengths that the grid barely resolves, and the error
6+% grows with distance travelled: a pulse made of many wavelengths slowly comes
7+% apart, its short components lagging behind, leaving a trailing ripple that
8+% is entirely numerical. The 5-point stencil's phase error goes as (k*h)^2,
9+% the 9-point one's as (k*h)^4, so at the eight or ten points per wavelength
10+% where a run is cheap enough to be interactive, the difference is easy to
11+% see: raise the frequency until models/leapfrog.m visibly trails, then switch
12+% to this one.
13+%
14+% It is not free. The stencil reads twice as many neighbours, and its
15+% stability limit is tighter (c*dt/h <= sqrt(3/8) rather than 1/sqrt(2)), so
16+% the app also takes a slightly smaller timestep. Roughly: a little more work
17+% per step, in exchange for needing a much coarser grid at the same accuracy.
18+
19+function [p, pm, t] = init(npts)
20+ p = zeros(npts, 1);
21+ pm = zeros(npts, 1);
22+ t = zeros(npts, 1);
23+end
24+
25+function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
26+ u = (t - t0) / tw;
27+ env = (1 - cw) * exp(-u .* u) + cw * (0.5 * (1 + tanh(u)));
28+ gx = ((x - x0) / w) .^ 2;
29+ gy = point * (((y - y0) / w) .^ 2);
30+ om = 2*pi*f;
31+ s = (om*om) * (env .* sin(om*(t - t0)) .* exp(-(gx + gy)));
32+
33+ sd = sig * dt;
34+ lap = lap4(p);
35+ pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);
36+
37+ pold = p;
38+ tn = t + dt;
39+end
package-lock.jsonadded+3017−0View file
This diff is 3,022 lines long and is not shown.
package.jsonadded+31−0View file
@@ -0,0 +1,31 @@
1+{
2+ "name": "acoustic-scattering-2d",
3+ "version": "0.1.0",
4+ "description": "2D acoustic scattering solved live in the browser: a MATLAB time-stepper compiled to WebGPU compute kernels",
5+ "type": "module",
6+ "engines": {
7+ "node": ">=22.6"
8+ },
9+ "license": "Apache-2.0",
10+ "scripts": {
11+ "dev": "vite",
12+ "build": "tsc --noEmit && vite build",
13+ "test:node": "vite-node scripts/test-node.ts",
14+ "smoke": "vite build && node scripts/smoke.mjs",
15+ "test": "npm run test:node && npm run smoke"
16+ },
17+ "dependencies": {
18+ "numbl": "file:../../numbl"
19+ },
20+ "optionalDependencies": {
21+ "webgpu": "^0.4.0"
22+ },
23+ "devDependencies": {
24+ "@types/node": "^26.1.1",
25+ "@webgpu/types": "^0.1.44",
26+ "puppeteer-core": "^23.11.1",
27+ "typescript": "^5.5.0",
28+ "vite": "^5.4.0",
29+ "vite-node": "^6.0.0"
30+ }
31+}
scenes/disk.madded+28−0View file
@@ -0,0 +1,28 @@
1+% A single circular scatterer in a uniform background: the reference case.
2+%
3+% In two dimensions this is the classical "scattering by a circular cylinder"
4+% problem — the 2D problem is the 3D one for fields that do not vary along z,
5+% so a disk in the (x, y) plane is the cross-section of an infinite cylinder.
6+% That is the object with an exact series solution in Bessel and Hankel
7+% functions, which makes it the thing to look at first.
8+%
9+% Sound speed only, at constant density, so the impedance ratio across the
10+% interface is just the speed ratio cin: cin much greater than 1 behaves
11+% nearly rigid (sound-hard), cin much less than 1 nearly pressure-release
12+% (sound-soft), and cin = 1 is no scatterer at all. `c0` is the speed of sound
13+% in air (343 m/s), which the app supplies; `cin` is a plain ratio to it, so
14+% the slider means the same thing whatever c0 is.
15+%
16+% `R` and `absorb` are in the app's own units — metres and inverse seconds —
17+% since they are lengths and rates, not ratios to anything.
18+%
19+% The interface is smoothed over about a cell so that the discrete medium is
20+% resolved by the grid. A jump between two neighbouring cells is not: it
21+% scatters the grid's own staircase rather than a circle.
22+function [c, sig] = medium(x, y, L, h, c0, cin, R, absorb)
23+ r = sqrt(x.^2 + y.^2);
24+ inside = 0.5 * (1 - tanh((r - R) / (1.5*h)));
25+
26+ c = c0 * (1 + (cin - 1) * inside);
27+ sig = sponge(x, y, L, 0.2*L, 1700) + absorb * inside;
28+end
scenes/lens.madded+18−0View file
@@ -0,0 +1,18 @@
1+% A gradient-index lens: no interface at all, just a region where sound
2+% travels more slowly.
3+%
4+% The speed dips smoothly to c0*(1 - dn) at the centre of a Gaussian of width
5+% `w` metres, so a plane wave crossing it is retarded most in the middle and
6+% the wavefront curves. The focus lands a little way past the lens and is
7+% much brighter than anything else on the grid, which is worth remembering
8+% when reading the colour scale.
9+%
10+% This is the case where a sharp-interface scatterer has nothing to say: there
11+% is no reflection to speak of, only refraction, and the field is smooth
12+% everywhere. It is also the easiest scene for the solver, since nothing on
13+% the grid is under-resolved.
14+function [c, sig] = medium(x, y, L, h, c0, dn, w, x0)
15+ r2 = ((x - x0).^2 + y.^2) / (w*w);
16+ c = c0 * (1 - dn * exp(-r2));
17+ sig = sponge(x, y, L, 0.2*L, 1700);
18+end
scenes/room.madded+46−0View file
@@ -0,0 +1,46 @@
1+% A room: a square enclosure with a doorway in one wall.
2+%
3+% Close the doorway (gap = 0) and it is a sealed cavity, where a pulse never
4+% leaves and the sound settles into the room's own modes. Open it and the room
5+% becomes a resonator that leaks: sound escapes through the aperture, and the
6+% ringing decays. This is the scene the microphone is for. Put the source and
7+% the microphone inside, listen to a click, and what comes back is a direct
8+% arrival followed by echoes closing up into reverberation.
9+%
10+% `side`, `gap` and `thick` are in metres — a real room, not a fraction of the
11+% domain — so `side` is the half-width and the room is `2*side` across.
12+%
13+% The walls are *slower* than the background rather than faster, and that is a
14+% deliberate trade. Reflection at an interface goes as |c2 - c1|/(c2 + c1)
15+% with the density constant, so a wall at cwall = 0.15 reflects about 75% of
16+% the pressure amplitude, as one at cwall = 7 would. But the timestep is set
17+% by the *fastest* speed anywhere on the grid, so a hard wall would make every
18+% step of the whole simulation several times smaller, while a slow one is
19+% free. What a slow wall costs instead is resolution inside itself: the
20+% wavelength there is shorter by the same factor, and the grid barely
21+% resolves it. That is what the absorption is for. It damps what gets into
22+% the wall so that the badly resolved part never comes back out, which is
23+% also what a real wall does.
24+%
25+% Turn `cwall` above 1 for a genuinely rigid room and watch the timestep in
26+% the readout drop to match.
27+function [c, sig] = medium(x, y, L, h, c0, cwall, side, gap, thick, absorb)
28+ s = 1.5*h;
29+
30+ % The shell between two concentric squares: |x| and |y| both within `side`
31+ % is the room, and out to side + thick is the wall.
32+ r = max(abs(x), abs(y));
33+ outer = 0.5 * (1 - tanh((r - (side + thick)) / s));
34+ inner = 0.5 * (1 - tanh((r - side) / s));
35+ wall = outer - inner;
36+
37+ % The doorway, cut out of the right-hand wall. Written so that gap = 0 shuts
38+ % it exactly: max(0, ...) is identically zero there, where a tanh of the same
39+ % thing would leave a half-open cell behind and quietly ruin the seal.
40+ through = min(1, max(0, gap/2 - abs(y)) / s);
41+ right = 0.5 * (1 + tanh((x - side) / s));
42+ wall = max(0, wall - through .* right);
43+
44+ c = c0 * (1 + (cwall - 1) * wall);
45+ sig = sponge(x, y, L, 0.2*L, 1700) + absorb * wall;
46+end
scenes/slit.madded+21−0View file
@@ -0,0 +1,21 @@
1+% A hard screen with two apertures: Young's two-slit experiment, in sound.
2+%
3+% The screen is a vertical slab of thickness `thick` at x = 0, with two gaps
4+% of width `gap` centred at y = +-sep/2 — all in metres. Its sound speed is
5+% `cwall`, a ratio to the background `c0`; well above 1 it is nearly rigid, so
6+% almost everything that does not go through a slit comes back.
7+%
8+% Drive it with the plane-wave source (the `point source` parameter at 0) and
9+% the two apertures become two sources radiating into the right half-plane.
10+% Where their paths differ by a whole wavelength they add; by half a
11+% wavelength they cancel. The fringe spacing is set by sep/wavelength, so
12+% raising the frequency packs the lobes closer together.
13+function [c, sig] = medium(x, y, L, h, c0, cwall, gap, sep, thick)
14+ s = 1.5*h;
15+ inslab = 0.5 * (1 - tanh((abs(x) - thick/2) / s));
16+ inslit = 0.5 * (1 - tanh((abs(abs(y) - sep/2) - gap/2) / s));
17+ wall = inslab .* (1 - inslit);
18+
19+ c = c0 * (1 + (cwall - 1) * wall);
20+ sig = sponge(x, y, L, 0.2*L, 1700);
21+end
scenes/speckle.madded+34−0View file
@@ -0,0 +1,34 @@
1+% A random smooth medium: many weak scatterers rather than one strong one.
2+%
3+% The speed field is a sum of 60 plane waves with random directions, random
4+% wavenumbers around kc (radians per metre), and random phases — a
5+% band-limited Gaussian random field, built here with an ordinary MATLAB loop
6+% because a scene is evaluated once on the CPU and has the whole language
7+% available, not the element-wise subset a compiled step is held to.
8+%
9+% A pulse crossing this does not scatter once and leave: it scatters weakly
10+% and repeatedly, and what comes out the far side is a spread-out coda rather
11+% than a clean wavefront. Turn `amp` down and the medium becomes transparent;
12+% turn it up and the direct arrival disappears into the noise.
13+%
14+% The field is windowed to the interior so the absorbing layer stays uniform.
15+% Structure inside the sponge would scatter sound back into the domain from
16+% the very place that is supposed to be swallowing it.
17+function [c, sig] = medium(x, y, L, npts, c0, amp, kc, seed)
18+ nm = 60;
19+ rng(seed);
20+ g = zeros(npts, 1);
21+ for k = 1:nm
22+ th = 2*pi*rand();
23+ kk = kc * (0.5 + rand());
24+ ph = 2*pi*rand();
25+ g = g + cos(kk * (cos(th)*x + sin(th)*y) + ph);
26+ end
27+ g = g / sqrt(nm);
28+
29+ r = sqrt(x.^2 + y.^2);
30+ window = 0.5 * (1 - tanh((r - 0.28*L) / (0.03*L)));
31+
32+ c = c0 * (1 + amp * (g .* window));
33+ sig = sponge(x, y, L, 0.2*L, 1700);
34+end
scripts/nodeWebGpu.tsadded+71−0View file
@@ -0,0 +1,71 @@
1+/**
2+ * Desktop WebGPU for the command-line scripts, via the optional `webgpu`
3+ * package (prebuilt Google Dawn).
4+ *
5+ * Installs Dawn under the globals the transform code expects (navigator.gpu,
6+ * GPUBufferUsage, ...) so everything under src/ runs here unchanged —
7+ * including requestShtDevice(), which makes the same device request the
8+ * browser makes.
9+ */
10+
11+export const errMsg = (e: unknown): string =>
12+ e instanceof Error ? e.message : String(e);
13+
14+/**
15+ * Returns a human-readable runtime description. The import specifier is
16+ * indirect so typechecking does not require the optional package.
17+ */
18+export async function installWebGpu(): Promise<string> {
19+ const specifier = 'webgpu';
20+ let mod: {
21+ create: (flags: string[]) => GPU;
22+ globals: Record<string, unknown>;
23+ };
24+ try {
25+ mod = await import(specifier);
26+ } catch (e) {
27+ // Distinguish "not installed" from "installed but the prebuilt Dawn binary
28+ // will not load" — the second is what a machine missing a system library
29+ // looks like, and reporting it as the first sends people in circles.
30+ const detail = errMsg(e);
31+ if (/Cannot find (package|module) '?webgpu'?/.test(detail)) {
32+ throw new Error(
33+ 'desktop WebGPU needs the optional `webgpu` package (prebuilt Google Dawn):\n' +
34+ ' npm install webgpu\n' +
35+ 'It is an optionalDependency, so npm can skip it silently — `npm ls webgpu`\n' +
36+ 'says whether it is there.',
37+ );
38+ }
39+ const glibc = /GLIBC_([0-9.]+)/.exec(detail);
40+ throw new Error(
41+ `the \`webgpu\` package is installed but did not load:\n ${detail}\n` +
42+ (glibc
43+ ? `Dawn's prebuilt binary wants glibc ${glibc[1]} or newer and this host is older\n` +
44+ '(`ldd --version` says how old). No flag bridges that — use a container with a\n' +
45+ 'newer base image, or a newer host.\n'
46+ : 'That is usually the prebuilt Dawn binary missing a system library.\n'),
47+ );
48+ }
49+ Object.assign(globalThis, mod.globals);
50+ // DAWN_FLAGS is ';'-separated because individual Dawn options take
51+ // comma-separated lists, e.g. 'enable-dawn-features=allow_unsafe_apis,...'
52+ const dawnFlags = process.env.DAWN_FLAGS?.split(';').filter(Boolean) ?? [];
53+ Object.defineProperty(globalThis, 'navigator', {
54+ value: { gpu: mod.create(dawnFlags) },
55+ configurable: true,
56+ writable: true,
57+ });
58+ const { version } = await import(`${specifier}/package.json`, {
59+ with: { type: 'json' },
60+ }).then(
61+ (m) => m.default as { version: string },
62+ () => ({ version: '?' }),
63+ );
64+ return `node-webgpu ${version} (Google Dawn)`;
65+}
66+
67+/** The hint to print when Dawn loads but finds no adapter. */
68+export const NO_ADAPTER_HINT =
69+ ' Dawn reaches the GPU through Vulkan on Linux and Windows, Metal on macOS,\n' +
70+ " so a headless box may have no adapter at all. DAWN_FLAGS='backend=vulkan'\n" +
71+ ' makes it explain itself.';
scripts/smoke.mjsadded+258−0View file
@@ -0,0 +1,258 @@
1+/**
2+ * Does the page actually run in a browser?
3+ *
4+ * Serves dist/ and opens it in headless Chrome (hardware WebGPU if there is
5+ * any, SwiftShader otherwise), then waits for the simulation to report that it
6+ * has taken steps and drawn frames. Any console error, page error or failed
7+ * request fails the run.
8+ *
9+ * This checks that the app *works*, not that it looks right — a headless
10+ * browser has no opinion about whether the wavefronts are in the right place.
11+ * Run it after `vite build`: node scripts/smoke.mjs
12+ *
13+ * `--full` additionally swaps the model, edits the source, and plays back the
14+ * microphone. All three need new GPU work started while the page is already
15+ * drawing — a pipeline, a buffer mapping — and headless Chrome cannot always
16+ * do that: on some machines everything of the sort fails with "A valid
17+ * external Instance reference no longer exists", in a thirty-line WebGPU page
18+ * with none of this project in it. So those steps are opt-in, and a failure
19+ * from them says as much about the browser as about the app. What they would
20+ * have covered on the solver side is covered by `npm run test:node`, which
21+ * runs against desktop WebGPU where readback works.
22+ */
23+const withGpuWork = process.argv.includes('--full');
24+import { createServer } from 'node:http';
25+import { readFile } from 'node:fs/promises';
26+import { extname, join } from 'node:path';
27+import puppeteer from 'puppeteer-core';
28+
29+const DIST = new URL('../dist/', import.meta.url).pathname;
30+const CHROME = process.env.CHROME_PATH ?? '/usr/bin/google-chrome';
31+const MIME = {
32+ '.html': 'text/html',
33+ '.js': 'text/javascript',
34+ '.css': 'text/css',
35+ '.json': 'application/json',
36+};
37+
38+const server = createServer(async (req, res) => {
39+ try {
40+ const path = req.url === '/' ? '/index.html' : req.url.split('?')[0];
41+ const data = await readFile(join(DIST, path));
42+ res.writeHead(200, { 'content-type': MIME[extname(path)] ?? 'application/octet-stream' });
43+ res.end(data);
44+ } catch {
45+ res.writeHead(404);
46+ res.end('not found');
47+ }
48+});
49+await new Promise((r) => server.listen(0, '127.0.0.1', r));
50+const port = server.address().port;
51+
52+const flagSets = [
53+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--enable-features=Vulkan'],
54+ [
55+ '--headless=new',
56+ '--no-sandbox',
57+ '--enable-unsafe-webgpu',
58+ '--use-webgpu-adapter=swiftshader',
59+ '--enable-unsafe-swiftshader',
60+ ],
61+];
62+
63+let ok = false;
64+let lastFailure = 'never ran';
65+for (const flags of flagSets) {
66+ const browser = await puppeteer.launch({
67+ executablePath: CHROME,
68+ args: [...flags],
69+ protocolTimeout: 600_000,
70+ });
71+ const problems = [];
72+ try {
73+ const page = await browser.newPage();
74+ page.on('console', (m) => {
75+ if (m.type() === 'error') problems.push(`console: ${m.text()}`);
76+ });
77+ page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`));
78+ page.on('requestfailed', (r) => problems.push(`request failed: ${r.url()}`));
79+
80+ await page.goto(`http://127.0.0.1:${port}/index.html`, { waitUntil: 'load' });
81+
82+ // The page loads paused. Wait for the first compile to have produced an op
83+ // list, then start it.
84+ await page.waitForFunction(
85+ () => (document.getElementById('compiled')?.textContent ?? '').includes('stencil'),
86+ { timeout: 180_000, polling: 500 },
87+ );
88+ // A paused page reports no frame rate, because it takes no steps and draws
89+ // nothing.
90+ await page.waitForFunction(
91+ () => /paused/.test(document.getElementById('stats')?.textContent ?? ''),
92+ { timeout: 60_000, polling: 250 },
93+ );
94+ if (/ms\/frame/.test(await page.$eval('#stats', (n) => n.textContent ?? ''))) {
95+ problems.push('a paused page reported a frame rate');
96+ }
97+ await page.click('#runpause');
98+
99+ /** The stats line reports the step count, so waiting on it says both that
100+ * the solver ran and that the frame loop is turning. */
101+ const running = async (label) => {
102+ await page.waitForFunction(
103+ () => {
104+ const m = /step (\d+)/.exec(document.getElementById('stats')?.textContent ?? '');
105+ return m ? Number(m[1]) > 20 : false;
106+ },
107+ { timeout: 180_000, polling: 500 },
108+ );
109+ console.log(` ${label}: ${(await page.$eval('#stats', (n) => n.textContent)).trim()}`);
110+ };
111+ const compiled = () => page.$eval('#compiled', (n) => n.textContent ?? '');
112+ const errText = () => page.$eval('#err', (n) => n.textContent ?? '');
113+
114+ await running('start');
115+
116+ // The timestep slider must reach the solver: halving the CFL fraction
117+ // halves dt, which the stats line reports.
118+ const dtNow = async () =>
119+ Number(/dt = ([0-9.e+-]+)/.exec(await page.$eval('#stats', (n) => n.textContent ?? ''))?.[1]);
120+ const dtBefore = await dtNow();
121+ await page.evaluate(() => {
122+ const slider = document.getElementById('cfl');
123+ slider.value = '0.25';
124+ slider.dispatchEvent(new Event('input'));
125+ });
126+ await page.waitForFunction(
127+ (before) => {
128+ const m = /dt = ([0-9.e+-]+)/.exec(document.getElementById('stats')?.textContent ?? '');
129+ return m ? Math.abs(Number(m[1]) - before / 2) < before / 20 : false;
130+ },
131+ { timeout: 30_000, polling: 250 },
132+ dtBefore,
133+ );
134+ console.log(` timestep: dt ${dtBefore.toExponential(2)} -> ${(await dtNow()).toExponential(2)} at half the CFL`);
135+
136+ // Clicking the picture must put the microphone there. The canvas shows
137+ // the whole domain, so the centre of it is the origin.
138+ const micReadout = () =>
139+ page.$$eval('#micparams output', (nodes) => nodes.map((n) => Number(n.textContent)));
140+ const box = await page.$eval('#view', (n) => {
141+ const r = n.getBoundingClientRect();
142+ return { x: r.x, y: r.y, w: r.width, h: r.height };
143+ });
144+ await page.mouse.click(box.x + box.w / 2, box.y + box.h / 2);
145+ await new Promise((r) => setTimeout(r, 300));
146+ const [mx, my] = await micReadout();
147+ console.log(` microphone dragged to (${mx}, ${my}) by clicking the centre`);
148+ if (Math.abs(mx) > 0.05 || Math.abs(my) > 0.05) {
149+ problems.push(`clicking the centre put the microphone at (${mx}, ${my})`);
150+ }
151+
152+ // The microphone must be recording, one sample per timestep, and its
153+ // trace must come back off the GPU when asked for.
154+ await page.waitForFunction(
155+ () => /(\d[\d,]*) samples/.test(document.getElementById('recinfo')?.textContent ?? ''),
156+ { timeout: 60_000, polling: 250 },
157+ );
158+ console.log(` microphone: ${(await page.$eval('#recinfo', (n) => n.textContent ?? '')).trim()}`);
159+ const steps = Number(
160+ /step (\d+)/.exec(await page.$eval('#stats', (n) => n.textContent ?? ''))?.[1],
161+ );
162+ const recorded = Number(
163+ /([\d,]+) samples/.exec(await page.$eval('#recinfo', (n) => n.textContent ?? ''))?.[1]
164+ .replace(/,/g, ''),
165+ );
166+ if (!(recorded > 0) || recorded > steps) {
167+ problems.push(`microphone recorded ${recorded} samples in ${steps} steps`);
168+ }
169+ if (withGpuWork) {
170+ await page.click('#listen');
171+ await new Promise((r) => setTimeout(r, 2000));
172+ const listenErr = await errText();
173+ if (listenErr) problems.push(`Listen reported: ${listenErr}`);
174+ }
175+
176+ // Swapping the scene re-evaluates its .m and re-uploads the medium. The
177+ // two-slit screen is much faster than the background, so the reported
178+ // speed range is the evidence.
179+ await page.select('#scene', 'slit');
180+ await page.waitForFunction(
181+ () => /c ∈ \[1, 3\.97\]/.test(document.getElementById('stats')?.textContent ?? ''),
182+ { timeout: 60_000, polling: 250 },
183+ );
184+ await running('after scene swap');
185+
186+ if (!withGpuWork) {
187+ console.log(' (skipping the recompile and playback checks; pass --full to run them)');
188+ } else {
189+ // Swapping the model recompiles, and the fourth-order one must reach for
190+ // the other stencil.
191+ await page.select('#model', 'leapfrog4');
192+ await page.waitForFunction(
193+ () => (document.getElementById('compiled')?.textContent ?? '').includes('lap4'),
194+ { timeout: 120_000, polling: 250 },
195+ );
196+ await running('after model swap');
197+
198+ // A broken edit must be reported rather than thrown, and must not take
199+ // the page down with it.
200+ await page.evaluate(() => {
201+ const ta = document.getElementById('source');
202+ ta.value = ta.value.replace('lap4(p)', 'lap9(p)');
203+ ta.dispatchEvent(new Event('input'));
204+ });
205+ await page.click('#recompile');
206+ await page.waitForFunction(
207+ () => (document.getElementById('err')?.textContent ?? '').length > 0,
208+ { timeout: 120_000, polling: 250 },
209+ );
210+ console.log(` bad edit reported: ${(await errText()).split('\n')[0]}`);
211+
212+ // And reverting must put it back.
213+ await page.click('#revert');
214+ await page.waitForFunction(
215+ () => (document.getElementById('err')?.textContent ?? '').length === 0,
216+ { timeout: 120_000, polling: 250 },
217+ );
218+ await running('after revert');
219+ if ((await compiled()).includes('lap4') === false) {
220+ problems.push('the reverted model did not recompile');
221+ }
222+ }
223+ const report = await page.evaluate(() => ({
224+ stats: document.getElementById('stats')?.textContent ?? '',
225+ err: document.getElementById('err')?.textContent ?? '',
226+ compiled: document.getElementById('compiled')?.textContent ?? '',
227+ // The canvas must have painted something other than the clear colour.
228+ painted: (() => {
229+ const canvas = document.getElementById('view');
230+ return canvas instanceof HTMLCanvasElement && canvas.width > 0;
231+ })(),
232+ }));
233+
234+ console.log(`flags: ${flags.join(' ')}`);
235+ console.log(report.stats.trim());
236+ console.log(report.compiled.split('\n').slice(0, 20).join('\n'));
237+ if (report.err) problems.push(`page error box: ${report.err}`);
238+ if (!report.painted) problems.push('canvas was never sized');
239+ if (problems.length === 0) {
240+ ok = true;
241+ } else {
242+ lastFailure = problems.join('\n');
243+ }
244+ } catch (e) {
245+ lastFailure = [`${e}`, ...problems].join('\n');
246+ } finally {
247+ await browser.close();
248+ }
249+ if (ok) break;
250+}
251+
252+server.close();
253+if (!ok) {
254+ console.error(`smoke: FAILED\n${lastFailure}`);
255+ process.exit(1);
256+}
257+console.log('smoke: the page runs');
258+process.exit(0);
scripts/test-node.tsadded+65−0View file
@@ -0,0 +1,65 @@
1+/**
2+ * The suite on desktop WebGPU (Google Dawn), against the real pipeline:
3+ * MATLAB source -> numbl lowering -> generated WGSL -> GPU.
4+ *
5+ * Run through vite-node, which is what resolves numbl's compiler sources and
6+ * the `?raw` .m imports:
7+ *
8+ * npm run test:node
9+ */
10+import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
11+import { requestAcousticDevice } from '../src/device.ts';
12+import {
13+ stencilChecks,
14+ propagationChecks,
15+ boundaryChecks,
16+ scatteringChecks,
17+ stabilityChecks,
18+ splitChecks,
19+ microphoneChecks,
20+ roomChecks,
21+ planChecks,
22+} from '../test/checks.ts';
23+
24+let failures = 0;
25+const check = (name: string, ok: boolean, detail: string): void => {
26+ console.log(`${ok ? 'PASS' : 'FAIL'} ${name} — ${detail}`);
27+ if (!ok) failures++;
28+};
29+const log = (s: string): void => console.log(s);
30+
31+/**
32+ * These checks compile MATLAB to compute shaders, so they need a GPU.
33+ * `--skip-without-gpu` lets a runner that has none say so and move on; a plain
34+ * local run still fails loudly, so a missing GPU is never mistaken for a pass.
35+ */
36+const skipWithoutGpu = process.argv.includes('--skip-without-gpu');
37+
38+let runtime: string;
39+let device: GPUDevice;
40+try {
41+ runtime = await installWebGpu();
42+ device = await requestAcousticDevice();
43+} catch (e) {
44+ const detail = `${errMsg(e)}\n${NO_ADAPTER_HINT}`;
45+ if (skipWithoutGpu) {
46+ console.log(`SKIP no WebGPU available here, so these checks did not run.\n${detail}`);
47+ process.exit(0);
48+ }
49+ console.error(`test-node: ${detail}`);
50+ process.exit(1);
51+}
52+console.log(`acoustic-scattering-2d tests — ${runtime}\n`);
53+
54+await stencilChecks(device, check, log);
55+await propagationChecks(device, check, log);
56+await boundaryChecks(device, check, log);
57+await scatteringChecks(device, check, log);
58+await stabilityChecks(device, check, log);
59+await splitChecks(device, check, log);
60+await microphoneChecks(device, check, log);
61+await roomChecks(device, check, log);
62+await planChecks(device, check, log);
63+
64+console.log(failures ? `\n${failures} failure(s)` : '\nall checks passed');
65+process.exit(failures ? 1 : 0);
src/audio/play.tsadded+98−0View file
@@ -0,0 +1,98 @@
1+/**
2+ * Turning a recorded trace into something you can hear.
3+ *
4+ * Because the solver is in real SI units, `dt` is a real duration in seconds
5+ * — so the natural playback rate is just 1/dt, sample for sample. Nothing is
6+ * reinterpreted: the trace plays back at the same real-time pace, and the
7+ * same pitch, that a microphone sitting at the probe point would have heard.
8+ * This is more than a convenience. A Courant-limited timestep on a grid fine
9+ * enough to resolve audible frequencies lands, by construction, in the same
10+ * range as an audio sample rate — the app's default settings give dt around
11+ * 20 microseconds, a rate near 48 kHz, which is not a coincidence: both are
12+ * set by "resolve a few centimetres of wave at audio frequency."
13+ *
14+ * `dt` can still land outside what a browser's AudioContext will accept, at
15+ * an unusual grid or CFL setting, so the rate is clamped and the caller is
16+ * told when that happened — the trace then plays sped up or slowed down
17+ * rather than at its real pace, which is worth knowing rather than hiding.
18+ *
19+ * Nothing here is a physical claim about loudness. The trace is normalized so
20+ * that whatever was recorded is audible, which discards exactly the quantity
21+ * (absolute amplitude) that the colour scale already shows.
22+ */
23+
24+/** What a browser will accept as an AudioBuffer sample rate. The spec's range
25+ * is wider than any of this needs; these bounds keep the derived rate inside
26+ * what every implementation supports. */
27+const MIN_RATE = 8000;
28+const MAX_RATE = 192000;
29+
30+export interface PlaybackPlan {
31+ /** Samples per second the trace is played back at. */
32+ rate: number;
33+ /** Seconds of audio. */
34+ duration: number;
35+ /** True if `rate` is the real 1/dt — false if it had to be clamped, in
36+ * which case playback runs faster or slower than the simulation did. */
37+ realTime: boolean;
38+}
39+
40+/** How a recorded trace of `samples` taken at timestep `dt` would be played. */
41+export function planPlayback(samples: number, dt: number): PlaybackPlan {
42+ const wanted = 1 / Math.max(dt, 1e-12);
43+ const rate = Math.min(MAX_RATE, Math.max(MIN_RATE, wanted));
44+ return { rate, duration: samples / rate, realTime: rate === wanted };
45+}
46+
47+let context: AudioContext | null = null;
48+let playing: AudioBufferSourceNode | null = null;
49+
50+/**
51+ * Play a recorded trace. Returns what was actually played, so the caller can
52+ * report it.
53+ *
54+ * Normalized to peak amplitude, and with a few milliseconds of fade at each
55+ * end: a trace that starts or ends away from zero is a step, and a step is a
56+ * click that has nothing to do with the simulation.
57+ */
58+export async function playTrace(
59+ trace: Float32Array,
60+ plan: PlaybackPlan,
61+): Promise<PlaybackPlan> {
62+ if (trace.length === 0) throw new Error('nothing has been recorded yet');
63+ context ??= new AudioContext();
64+ if (context.state === 'suspended') await context.resume();
65+
66+ const buffer = context.createBuffer(1, trace.length, plan.rate);
67+ const channel = buffer.getChannelData(0);
68+ let peak = 0;
69+ for (const v of trace) peak = Math.max(peak, Math.abs(v));
70+ const gain = peak > 0 ? 0.9 / peak : 0;
71+ for (let i = 0; i < trace.length; i++) channel[i] = trace[i] * gain;
72+
73+ const fade = Math.min(Math.round(0.005 * plan.rate), Math.floor(trace.length / 2));
74+ for (let i = 0; i < fade; i++) {
75+ const w = i / fade;
76+ channel[i] *= w;
77+ channel[trace.length - 1 - i] *= w;
78+ }
79+
80+ stop();
81+ const source = context.createBufferSource();
82+ source.buffer = buffer;
83+ source.connect(context.destination);
84+ source.onended = () => {
85+ if (playing === source) playing = null;
86+ };
87+ source.start();
88+ playing = source;
89+ return plan;
90+}
91+
92+export function stop(): void {
93+ if (!playing) return;
94+ playing.stop();
95+ playing = null;
96+}
97+
98+export const isPlaying = (): boolean => playing !== null;
src/audio/recorder.tsadded+226−0View file
@@ -0,0 +1,226 @@
1+/**
2+ * A microphone: the pressure at one grid point, sampled every timestep.
3+ *
4+ * The obvious implementation — read the field back and pick out one number —
5+ * costs a GPU-to-CPU round trip per step, which is more than the step itself.
6+ * So the trace is written on the GPU instead, by a one-thread dispatch that
7+ * runs after each step and appends `p` at the probe point to a buffer. The
8+ * whole trace comes back to the CPU once, when there is something to listen
9+ * to.
10+ *
11+ * One sample per timestep is the natural rate: it is every value the
12+ * simulation has, and nothing is being resampled or interpolated on the way
13+ * in. What that means in seconds is decided at playback (src/audio/play.ts),
14+ * because the simulation has no seconds in it — only model time.
15+ */
16+
17+/** Samples the trace holds: about 12 seconds of audio at a typical playback
18+ * rate, and 1 MB of GPU memory. Recording stops when it is full rather than
19+ * wrapping, so what you hear always starts where the run did. */
20+export const TRACE_CAPACITY = 1 << 18;
21+
22+const SHADER = `
23+struct Probe {
24+ index: u32,
25+ capacity: u32,
26+};
27+
28+@group(0) @binding(0) var<storage, read_write> trace: array<f32>;
29+@group(0) @binding(1) var<storage, read_write> head: array<u32>;
30+@group(0) @binding(2) var<storage, read> field: array<f32>;
31+@group(0) @binding(3) var<uniform> probe: Probe;
32+
33+// One invocation, so the read-modify-write of the head needs no atomic.
34+@compute @workgroup_size(1)
35+fn main() {
36+ let i = head[0];
37+ if (i < probe.capacity) {
38+ trace[i] = field[probe.index];
39+ head[0] = i + 1u;
40+ }
41+}
42+`;
43+
44+export interface RecorderOptions {
45+ device: GPUDevice;
46+ /** The pressure buffer to sample — the host-owned one, which is what both
47+ * `init` and `step` leave their result in. */
48+ field: GPUBuffer;
49+ nx: number;
50+ ny: number;
51+}
52+
53+export class Recorder {
54+ readonly capacity = TRACE_CAPACITY;
55+
56+ #device: GPUDevice;
57+ #pipeline: GPUComputePipeline | null = null;
58+ #layout: GPUBindGroupLayout;
59+ #bindGroup: GPUBindGroup | null = null;
60+ #trace: GPUBuffer;
61+ #head: GPUBuffer;
62+ #probe: GPUBuffer;
63+ #readback: GPUBuffer;
64+ #field: GPUBuffer;
65+ #nx: number;
66+ #ny: number;
67+ /** Samples written since the last clear, as far as the host knows. Counted
68+ * here rather than read back from the GPU: the dispatch runs once per step
69+ * and the host knows exactly how many steps it asked for. */
70+ #count = 0;
71+ #reading = false;
72+
73+ private constructor(init: {
74+ device: GPUDevice;
75+ layout: GPUBindGroupLayout;
76+ trace: GPUBuffer;
77+ head: GPUBuffer;
78+ probe: GPUBuffer;
79+ readback: GPUBuffer;
80+ field: GPUBuffer;
81+ nx: number;
82+ ny: number;
83+ }) {
84+ this.#device = init.device;
85+ this.#layout = init.layout;
86+ this.#trace = init.trace;
87+ this.#head = init.head;
88+ this.#probe = init.probe;
89+ this.#readback = init.readback;
90+ this.#field = init.field;
91+ this.#nx = init.nx;
92+ this.#ny = init.ny;
93+ }
94+
95+ static async create(opts: RecorderOptions): Promise<Recorder> {
96+ const { device } = opts;
97+ const layout = device.createBindGroupLayout({
98+ label: 'recorder',
99+ entries: [
100+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
101+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
102+ {
103+ binding: 2,
104+ visibility: GPUShaderStage.COMPUTE,
105+ buffer: { type: 'read-only-storage' },
106+ },
107+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform' } },
108+ ],
109+ });
110+ const trace = device.createBuffer({
111+ label: 'recorder-trace',
112+ size: 4 * TRACE_CAPACITY,
113+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
114+ });
115+ const head = device.createBuffer({
116+ label: 'recorder-head',
117+ size: 4,
118+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
119+ });
120+ const probe = device.createBuffer({
121+ label: 'recorder-probe',
122+ size: 8,
123+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
124+ });
125+ const readback = device.createBuffer({
126+ label: 'recorder-readback',
127+ size: 4 * TRACE_CAPACITY,
128+ usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
129+ });
130+
131+ const rec = new Recorder({
132+ device, layout, trace, head, probe, readback,
133+ field: opts.field, nx: opts.nx, ny: opts.ny,
134+ });
135+ rec.#pipeline = await device.createComputePipelineAsync({
136+ label: 'recorder',
137+ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }),
138+ compute: {
139+ module: device.createShaderModule({ code: SHADER, label: 'recorder' }),
140+ entryPoint: 'main',
141+ },
142+ });
143+ rec.#bindGroup = device.createBindGroup({
144+ layout,
145+ entries: [
146+ { binding: 0, resource: { buffer: trace } },
147+ { binding: 1, resource: { buffer: head } },
148+ { binding: 2, resource: { buffer: opts.field } },
149+ { binding: 3, resource: { buffer: probe } },
150+ ],
151+ });
152+ rec.clear();
153+ return rec;
154+ }
155+
156+ /** Samples recorded so far. Stops rising once the trace is full. */
157+ get count(): number {
158+ return Math.min(this.#count, this.capacity);
159+ }
160+
161+ get full(): boolean {
162+ return this.#count >= this.capacity;
163+ }
164+
165+ /** Put the microphone at the grid point nearest (x, y). Free: the probe
166+ * index is a uniform, so moving it disturbs neither the run nor the
167+ * recording already made. */
168+ setProbe(ix: number, iy: number): void {
169+ const cx = Math.max(0, Math.min(this.#nx - 1, Math.round(ix)));
170+ const cy = Math.max(0, Math.min(this.#ny - 1, Math.round(iy)));
171+ this.#device.queue.writeBuffer(
172+ this.#probe,
173+ 0,
174+ new Uint32Array([cx + this.#nx * cy, this.capacity]),
175+ );
176+ }
177+
178+ /** Start again from an empty trace. */
179+ clear(): void {
180+ this.#count = 0;
181+ this.#device.queue.writeBuffer(this.#head, 0, new Uint32Array([0]));
182+ }
183+
184+ /** Record one sample. Called once per timestep, inside the step's own
185+ * submission, so no extra work crosses to the host. */
186+ encode(encoder: GPUCommandEncoder): void {
187+ if (!this.#pipeline || !this.#bindGroup || this.full) return;
188+ const pass = encoder.beginComputePass({ label: 'recorder' });
189+ pass.setPipeline(this.#pipeline);
190+ pass.setBindGroup(0, this.#bindGroup);
191+ pass.dispatchWorkgroups(1);
192+ pass.end();
193+ this.#count++;
194+ }
195+
196+ /** The recorded trace. The only readback the microphone ever does. */
197+ async read(): Promise<Float32Array> {
198+ const n = this.count;
199+ if (n === 0) return new Float32Array(0);
200+ if (this.#reading) throw new Error('a trace readback is already in flight');
201+ this.#reading = true;
202+ try {
203+ const enc = this.#device.createCommandEncoder({ label: 'recorder-read' });
204+ enc.copyBufferToBuffer(this.#trace, 0, this.#readback, 0, 4 * n);
205+ this.#device.queue.submit([enc.finish()]);
206+ await this.#readback.mapAsync(GPUMapMode.READ, 0, 4 * n);
207+ const out = new Float32Array(this.#readback.getMappedRange(0, 4 * n).slice(0));
208+ this.#readback.unmap();
209+ return out;
210+ } finally {
211+ this.#reading = false;
212+ }
213+ }
214+
215+ /** The pressure buffer this microphone listens to. */
216+ get field(): GPUBuffer {
217+ return this.#field;
218+ }
219+
220+ destroy(): void {
221+ this.#trace.destroy();
222+ this.#head.destroy();
223+ this.#probe.destroy();
224+ this.#readback.destroy();
225+ }
226+}
src/device.tsadded+45−0View file
@@ -0,0 +1,45 @@
1+/**
2+ * The GPU device, requested the same way everywhere (app, tests, scripts).
3+ *
4+ * One limit matters here. A fused kernel binds one storage buffer per distinct
5+ * grid field its line reads, plus its output and the parameter block, and a
6+ * leapfrog update reads a lot of fields at once — the two pressure histories,
7+ * the clock, the coordinates, the sound speed, the absorption, the Laplacian.
8+ * WebGPU only guarantees 8 storage buffers per compute stage, so we ask for
9+ * whatever the adapter will give up to 16. When that is not enough the planner
10+ * splits the kernel instead of failing (see `fitToBudget` in plan.ts), so this
11+ * is a performance request rather than a requirement.
12+ */
13+export const MAX_STORAGE_BUFFERS = 16;
14+
15+/** The adapter the device came from, kept so its limits and info stay
16+ * available for as long as the device is in use. */
17+let heldAdapter: GPUAdapter | null = null;
18+
19+export async function requestAcousticDevice(): Promise<GPUDevice> {
20+ if (!navigator.gpu) {
21+ throw new Error(
22+ 'this browser has no WebGPU. Chrome and Edge 113+, Safari 26+, and ' +
23+ 'Firefox 141+ on Windows have it; on Linux Firefox and Chrome may need ' +
24+ 'it enabled explicitly.',
25+ );
26+ }
27+ const adapter = await navigator.gpu.requestAdapter({ powerPreference: 'high-performance' });
28+ if (!adapter) throw new Error('WebGPU found no adapter on this machine.');
29+ heldAdapter = adapter;
30+ const want = Math.min(
31+ MAX_STORAGE_BUFFERS,
32+ adapter.limits.maxStorageBuffersPerShaderStage ?? 8,
33+ );
34+ return adapter.requestDevice({
35+ requiredLimits: { maxStorageBuffersPerShaderStage: want },
36+ });
37+}
38+
39+/** The adapter the current device came from, if there is one. */
40+export const currentAdapter = (): GPUAdapter | null => heldAdapter;
41+
42+/** Grid fields one kernel may read, given the device's binding limit: every
43+ * binding but the output and the parameter block. */
44+export const kernelOperandBudget = (device: GPUDevice): number =>
45+ Math.max(2, (device.limits.maxStorageBuffersPerShaderStage ?? 8) - 2);
src/editor/codeEditor.tsadded+92−0View file
@@ -0,0 +1,92 @@
1+/**
2+ * A textarea with syntax highlighting, by overlay.
3+ *
4+ * A textarea cannot colour its own text, so the highlighted source is rendered
5+ * into a <pre> underneath and the textarea sits on top with transparent text and
6+ * a visible caret. The two must agree on every metric that affects layout —
7+ * font, line height, padding, tab size, wrapping — and their scroll offsets are
8+ * kept in sync, or the colours drift away from the characters.
9+ */
10+import { highlightMatlab } from './matlab.ts';
11+
12+export interface CodeEditorOptions {
13+ textarea: HTMLTextAreaElement;
14+ /** The <pre> behind it, holding the highlighted copy. */
15+ overlay: HTMLElement;
16+ /** Names to mark as host-provided operations. */
17+ external?: ReadonlySet<string>;
18+ /** Called on every edit. */
19+ onInput?: (value: string) => void;
20+}
21+
22+export class CodeEditor {
23+ #textarea: HTMLTextAreaElement;
24+ #overlay: HTMLElement;
25+ #external: ReadonlySet<string>;
26+
27+ constructor(opts: CodeEditorOptions) {
28+ this.#textarea = opts.textarea;
29+ this.#overlay = opts.overlay;
30+ this.#external = opts.external ?? new Set();
31+
32+ this.#textarea.addEventListener('input', () => {
33+ this.#repaint();
34+ opts.onInput?.(this.#textarea.value);
35+ });
36+ // Keep the colours under the characters while scrolling.
37+ this.#textarea.addEventListener('scroll', () => this.#syncScroll());
38+ // Tab should indent rather than leave the editor.
39+ this.#textarea.addEventListener('keydown', (e) => this.#onKeyDown(e));
40+ this.#repaint();
41+ }
42+
43+ get value(): string {
44+ return this.#textarea.value;
45+ }
46+
47+ set value(next: string) {
48+ this.#textarea.value = next;
49+ this.#repaint();
50+ }
51+
52+ focus(): void {
53+ this.#textarea.focus();
54+ }
55+
56+ /** Select a character range, scrolling it into view. */
57+ select(start: number, end: number): void {
58+ this.#textarea.focus();
59+ this.#textarea.setSelectionRange(start, end);
60+ // setSelectionRange does not always scroll; nudge the line into view.
61+ const line = this.#textarea.value.slice(0, start).split('\n').length - 1;
62+ const lineHeight = this.#textarea.scrollHeight / Math.max(1, this.#lineCount());
63+ const target = line * lineHeight - this.#textarea.clientHeight / 2;
64+ this.#textarea.scrollTop = Math.max(0, target);
65+ this.#syncScroll();
66+ }
67+
68+ #lineCount(): number {
69+ return this.#textarea.value.split('\n').length + 1; // +1 for the trailing line
70+ }
71+
72+ #onKeyDown(e: KeyboardEvent): void {
73+ if (e.key !== 'Tab' || e.ctrlKey || e.metaKey || e.altKey) return;
74+ e.preventDefault();
75+ const el = this.#textarea;
76+ const { selectionStart: s, selectionEnd: t, value } = el;
77+ el.value = `${value.slice(0, s)} ${value.slice(t)}`;
78+ el.selectionStart = el.selectionEnd = s + 2;
79+ // Let the input listener repaint and notify, as for any other edit.
80+ el.dispatchEvent(new Event('input'));
81+ }
82+
83+ #repaint(): void {
84+ this.#overlay.innerHTML = highlightMatlab(this.#textarea.value, this.#external);
85+ this.#syncScroll();
86+ }
87+
88+ #syncScroll(): void {
89+ this.#overlay.scrollTop = this.#textarea.scrollTop;
90+ this.#overlay.scrollLeft = this.#textarea.scrollLeft;
91+ }
92+}
src/editor/matlab.tsadded+213−0View file
@@ -0,0 +1,213 @@
1+/**
2+ * A small MATLAB tokenizer, for syntax highlighting the model editor.
3+ *
4+ * Only what highlighting needs — comments, literals, numbers, keywords — and
5+ * deliberately not a parser: numbl does the real parsing, and reports errors
6+ * with positions. Tokens preserve the source text exactly, character for
7+ * character, because the highlighted output is overlaid on a textarea and any
8+ * dropped or added character would shift the two out of alignment.
9+ */
10+
11+export type TokenClass = 'com' | 'str' | 'num' | 'kw' | 'ext';
12+
13+export interface Token {
14+ text: string;
15+ cls: TokenClass | null;
16+}
17+
18+const KEYWORDS = new Set([
19+ 'break', 'case', 'catch', 'classdef', 'continue', 'else', 'elseif', 'end',
20+ 'for', 'function', 'global', 'if', 'otherwise', 'parfor', 'persistent',
21+ 'return', 'spmd', 'switch', 'try', 'while',
22+]);
23+
24+const isIdentStart = (c: string): boolean => /[A-Za-z_]/.test(c);
25+const isIdent = (c: string): boolean => /[A-Za-z0-9_]/.test(c);
26+const isDigit = (c: string): boolean => c >= '0' && c <= '9';
27+
28+/**
29+ * In MATLAB `'` is both the transpose operator and the char-literal delimiter.
30+ * It opens a literal unless it directly follows something that can be
31+ * transposed — a value, a closing bracket, or another transpose.
32+ */
33+function quoteIsTranspose(src: string, at: number): boolean {
34+ for (let i = at - 1; i >= 0; i--) {
35+ const c = src[i];
36+ if (c === ' ' || c === '\t') continue;
37+ return isIdent(c) || c === ')' || c === ']' || c === '}' || c === '.' || c === "'";
38+ }
39+ return false;
40+}
41+
42+/**
43+ * Tokenize `src`. `external` names (the operations the host provides, e.g.
44+ * `synth` / `analys`) get their own class so the boundary between the model and
45+ * what it is given is visible in the editor.
46+ */
47+export function tokenizeMatlab(
48+ src: string,
49+ external: ReadonlySet<string> = new Set(),
50+): Token[] {
51+ const out: Token[] = [];
52+ const push = (text: string, cls: TokenClass | null): void => {
53+ if (!text) return;
54+ const last = out[out.length - 1];
55+ if (last && last.cls === cls) last.text += text;
56+ else out.push({ text, cls });
57+ };
58+
59+ let i = 0;
60+ let atLineStart = true;
61+ let inBlockComment = false;
62+
63+ while (i < src.length) {
64+ const c = src[i];
65+
66+ // Block comments: `%{` and `%}` each alone on their line.
67+ if (atLineStart) {
68+ const eol = src.indexOf('\n', i);
69+ const lineEnd = eol === -1 ? src.length : eol;
70+ const line = src.slice(i, lineEnd);
71+ const trimmed = line.trim();
72+ if (!inBlockComment && trimmed === '%{') inBlockComment = true;
73+ else if (inBlockComment && trimmed === '%}') {
74+ push(line, 'com');
75+ i = lineEnd;
76+ inBlockComment = false;
77+ atLineStart = false;
78+ continue;
79+ }
80+ if (inBlockComment) {
81+ push(line, 'com');
82+ i = lineEnd;
83+ atLineStart = false;
84+ continue;
85+ }
86+ }
87+
88+ if (c === '\n') {
89+ push(c, null);
90+ i++;
91+ atLineStart = true;
92+ continue;
93+ }
94+ if (c === ' ' || c === '\t') {
95+ push(c, null);
96+ i++;
97+ continue;
98+ }
99+ atLineStart = false;
100+
101+ // Line comment, including MATLAB's `%%` section markers.
102+ if (c === '%') {
103+ const eol = src.indexOf('\n', i);
104+ const end = eol === -1 ? src.length : eol;
105+ push(src.slice(i, end), 'com');
106+ i = end;
107+ continue;
108+ }
109+
110+ // Line continuation is an operator, but any trailing text is a comment.
111+ if (c === '.' && src.startsWith('...', i)) {
112+ const eol = src.indexOf('\n', i);
113+ const end = eol === -1 ? src.length : eol;
114+ push('...', null);
115+ push(src.slice(i + 3, end), 'com');
116+ i = end;
117+ continue;
118+ }
119+
120+ // Char literal (or transpose).
121+ if (c === "'") {
122+ if (quoteIsTranspose(src, i)) {
123+ push("'", null);
124+ i++;
125+ continue;
126+ }
127+ let j = i + 1;
128+ while (j < src.length && src[j] !== '\n') {
129+ if (src[j] === "'") {
130+ if (src[j + 1] === "'") j += 2; // escaped quote
131+ else {
132+ j++;
133+ break;
134+ }
135+ } else j++;
136+ }
137+ push(src.slice(i, j), 'str');
138+ i = j;
139+ continue;
140+ }
141+
142+ // Double-quoted string.
143+ if (c === '"') {
144+ let j = i + 1;
145+ while (j < src.length && src[j] !== '\n') {
146+ if (src[j] === '"') {
147+ if (src[j + 1] === '"') j += 2;
148+ else {
149+ j++;
150+ break;
151+ }
152+ } else j++;
153+ }
154+ push(src.slice(i, j), 'str');
155+ i = j;
156+ continue;
157+ }
158+
159+ // Number: 12, 1.5, .5, 1e-3, 2i
160+ if (isDigit(c) || (c === '.' && isDigit(src[i + 1]))) {
161+ let j = i;
162+ while (j < src.length && isDigit(src[j])) j++;
163+ if (src[j] === '.') {
164+ j++;
165+ while (j < src.length && isDigit(src[j])) j++;
166+ }
167+ if (src[j] === 'e' || src[j] === 'E') {
168+ let k = j + 1;
169+ if (src[k] === '+' || src[k] === '-') k++;
170+ if (isDigit(src[k])) {
171+ k++;
172+ while (k < src.length && isDigit(src[k])) k++;
173+ j = k;
174+ }
175+ }
176+ if (src[j] === 'i' || src[j] === 'j') j++;
177+ push(src.slice(i, j), 'num');
178+ i = j;
179+ continue;
180+ }
181+
182+ // Identifier / keyword / external operation.
183+ if (isIdentStart(c)) {
184+ let j = i;
185+ while (j < src.length && isIdent(src[j])) j++;
186+ const word = src.slice(i, j);
187+ push(word, KEYWORDS.has(word) ? 'kw' : external.has(word) ? 'ext' : null);
188+ i = j;
189+ continue;
190+ }
191+
192+ push(c, null);
193+ i++;
194+ }
195+
196+ return out;
197+}
198+
199+const escapeHtml = (s: string): string =>
200+ s.replace(/&/g, '&').replace(/</g, '<').replace(/>/g, '>');
201+
202+/** Highlighted HTML for `src`, safe to assign to innerHTML. */
203+export function highlightMatlab(
204+ src: string,
205+ external: ReadonlySet<string> = new Set(),
206+): string {
207+ const html = tokenizeMatlab(src, external)
208+ .map((t) => (t.cls ? `<span class="tok-${t.cls}">${escapeHtml(t.text)}</span>` : escapeHtml(t.text)))
209+ .join('');
210+ // A trailing newline keeps the last line's box height stable, so the overlay
211+ // and the textarea scroll to the same extent.
212+ return `${html}\n`;
213+}
src/grid.tsadded+76−0View file
@@ -0,0 +1,76 @@
1+/**
2+ * The computational grid: a square, cell-centred, uniform in both directions.
3+ *
4+ * Fields are flattened x-fastest — the point (ix, iy) is element `ix + nx*iy`
5+ * — which is the order the stencil shader indexes in and the order the
6+ * renderer reads a row of pixels in. Everything else in the project treats a
7+ * field as an opaque npts x 1 column vector.
8+ *
9+ * Cell-centred rather than node-centred so that no grid point sits exactly on
10+ * the outer boundary: the stencil takes the field outside the domain to be
11+ * zero, and the absorbing layer is meant to have swallowed the wave before it
12+ * gets there.
13+ */
14+export interface Grid {
15+ n: number;
16+ nx: number;
17+ ny: number;
18+ npts: number;
19+ /** Side length of the square domain, in metres, centred on the origin. */
20+ L: number;
21+ /** Grid spacing, L/n, in metres. */
22+ h: number;
23+ /** Coordinates of every point, npts each, x fastest — as the shaders see
24+ * them (f32) and as the scene .m is evaluated at (f64). */
25+ x: Float32Array;
26+ y: Float32Array;
27+ x64: Float64Array;
28+ y64: Float64Array;
29+}
30+
31+export function makeGrid(n: number, L: number): Grid {
32+ const h = L / n;
33+ const npts = n * n;
34+ const x64 = new Float64Array(npts);
35+ const y64 = new Float64Array(npts);
36+ for (let iy = 0; iy < n; iy++) {
37+ const yv = -L / 2 + (iy + 0.5) * h;
38+ for (let ix = 0; ix < n; ix++) {
39+ const k = ix + n * iy;
40+ x64[k] = -L / 2 + (ix + 0.5) * h;
41+ y64[k] = yv;
42+ }
43+ }
44+ return {
45+ n,
46+ nx: n,
47+ ny: n,
48+ npts,
49+ L,
50+ h,
51+ x: new Float32Array(x64),
52+ y: new Float32Array(y64),
53+ x64,
54+ y64,
55+ };
56+}
57+
58+/**
59+ * The timestep the explicit leapfrog is stable at.
60+ *
61+ * All in SI: `h` in metres, `cmax` in metres per second, the result in
62+ * seconds.
63+ *
64+ * Leapfrog on p_tt = c^2 L p is stable while dt^2 c^2 |L|max <= 4, and the
65+ * discrete Laplacian's extreme eigenvalue is what differs between the
66+ * stencils: 8/h^2 for the 5-point one, 32/(3 h^2) for the 9-point
67+ * fourth-order one. That gives c*dt/h <= 1/sqrt(2) = 0.707 and
68+ * c*dt/h <= sqrt(3/8) = 0.612 respectively. `cfl` is the fraction of that
69+ * limit to run at, and `cmax` is the fastest sound speed anywhere in the
70+ * medium — a scatterer faster than the background sets the timestep for the
71+ * whole grid.
72+ */
73+export function stableDt(h: number, cmax: number, order: 2 | 4, cfl = 0.5): number {
74+ const limit = order === 2 ? Math.SQRT1_2 : Math.sqrt(3 / 8);
75+ return (cfl * limit * h) / Math.max(cmax, 1e-12);
76+}
src/main.tsadded+701−0View file
@@ -0,0 +1,701 @@
1+/**
2+ * The app: two MATLAB files, a GPU, and a canvas.
3+ *
4+ * Everything the page does falls into three motions. Changing a *parameter*
5+ * writes a uniform, which is free and does not interrupt the run. Changing the
6+ * *scene* re-evaluates its .m on the CPU and re-uploads two arrays, which is
7+ * cheap and needs no recompile. Changing the *model*, the grid, or either
8+ * file's text recompiles — a fresh session, from source to shaders.
9+ */
10+import { requestAcousticDevice } from './device.ts';
11+import { ModelSession } from './mgpu/session.ts';
12+import { EXTERNAL_OPS } from './mgpu/externals.ts';
13+import { formatFailure, ModelCompileError } from './mgpu/errors.ts';
14+import {
15+ mModels,
16+ mModelByKey,
17+ defaultParams,
18+ type MModel,
19+ type Params,
20+ type ParamSpec,
21+} from './mgpu/registry.ts';
22+import {
23+ mScenes,
24+ mSceneByKey,
25+ defaultSceneParams,
26+ type MScene,
27+} from './scene/registry.ts';
28+import { FieldView } from './render/field.ts';
29+import { Colorbar, fmtValue } from './render/colorbar.ts';
30+import { colormaps, colormapNames } from './render/colormaps.ts';
31+import { CodeEditor } from './editor/codeEditor.ts';
32+import { planPlayback, playTrace, stop as stopAudio } from './audio/play.ts';
33+import { C_AIR, DOMAIN, POOR_RESOLUTION, fmtLength, fmtTime } from './units.ts';
34+
35+const el = <T extends HTMLElement>(id: string): T => {
36+ const node = document.getElementById(id);
37+ if (!node) throw new Error(`missing element #${id}`);
38+ return node as T;
39+};
40+
41+const errBox = el<HTMLParagraphElement>('err');
42+const showError = (e: unknown, source: string): void => {
43+ errBox.textContent = formatFailure(e, source);
44+};
45+const clearError = (): void => {
46+ errBox.textContent = '';
47+};
48+
49+/* ---------------------------------------------------------------- state -- */
50+
51+let model: MModel = mModels[0];
52+let scene: MScene = mScenes[0];
53+// The starting scene may ask for source settings and a microphone position,
54+// the same way it does when picked from the dropdown later.
55+let params: Params = { ...defaultParams(model), ...scene.suggest };
56+let sceneParams: Params = defaultSceneParams(scene);
57+/** The editor's working copies, which may differ from the presets. */
58+const sources = { model: model.source, scene: scene.source };
59+let gridN = 512;
60+let stepsPerFrame = 4;
61+/** Timestep, as a fraction of the largest one the scheme is stable at. */
62+let cfl = 0.5;
63+/** Where the microphone sits, in metres. */
64+let mic = { ...(scene.mic ?? { x: DOMAIN * 0.2, y: 0 }) };
65+/** Paused until asked. The page compiles and draws its silent initial state on
66+ * load, so what is on screen is the medium about to be sounded, and nothing
67+ * moves until Run. */
68+let running = false;
69+let session: ModelSession | null = null;
70+
71+/** Pressure the colormap saturates at, and whether it follows the field. */
72+let scale = 1;
73+let autoScale = true;
74+/** Largest pressure seen since the last restart. The colour scale is not
75+ * allowed to fall far below it, so that once a pulse has left the grid what
76+ * is drawn is an empty grid rather than roundoff at full contrast. */
77+let peakSeen = 0;
78+let colormapName = colormapNames[0];
79+let showMedium = true;
80+/**
81+ * Something on screen would change if we drew now.
82+ *
83+ * A running simulation is dirty every frame by definition. A paused one is
84+ * dirty only when told to be — a new colour scale, a new medium, a resize —
85+ * and otherwise draws nothing at all. Skipping the draw is safe because a
86+ * WebGPU canvas keeps the last frame it presented; without it a paused page
87+ * would still run a full-screen fragment pass sixty times a second, which is
88+ * real GPU work in aid of an unchanging picture.
89+ */
90+let dirty = true;
91+
92+/* ------------------------------------------------------------- the page -- */
93+
94+const device = await requestAcousticDevice().catch((e: unknown) => {
95+ showError(e, '');
96+ return null;
97+});
98+if (!device) throw new Error('no GPU');
99+
100+// A lost device takes everything with it and nothing afterwards will work, so
101+// say so rather than leaving a frozen picture and no explanation.
102+void device.lost.then((info) => {
103+ showError(
104+ new Error(
105+ `the GPU device was lost (${info.reason}): ${info.message}\n` +
106+ 'Reload the page to start again.',
107+ ),
108+ '',
109+ );
110+});
111+
112+const canvas = el<HTMLCanvasElement>('view');
113+const view = new FieldView({ device, canvas, nx: gridN, ny: gridN });
114+view.setColormap(colormaps[colormapName]);
115+const colorbar = new Colorbar(el('colorbar'));
116+colorbar.setColormap(colormaps[colormapName]);
117+
118+// The canvas is sized by CSS; match its backing store when that changes
119+// rather than measuring it every frame.
120+// Resizing reallocates the canvas's backing store, which clears it.
121+new ResizeObserver(() => {
122+ view.resize();
123+ dirty = true;
124+}).observe(canvas);
125+
126+const editor = new CodeEditor({
127+ textarea: el<HTMLTextAreaElement>('source'),
128+ overlay: el('highlight'),
129+ external: EXTERNAL_OPS,
130+ onInput: (value) => {
131+ sources[editorFile.value as 'model' | 'scene'] = value;
132+ el<HTMLButtonElement>('recompile').classList.add('primary');
133+ },
134+});
135+
136+const modelSelect = el<HTMLSelectElement>('model');
137+const sceneSelect = el<HTMLSelectElement>('scene');
138+const editorFile = el<HTMLSelectElement>('editor-file');
139+const colormapSelect = el<HTMLSelectElement>('colormap');
140+
141+for (const m of mModels) {
142+ modelSelect.append(new Option(m.label, m.key));
143+}
144+for (const s of mScenes) {
145+ sceneSelect.append(new Option(s.label, s.key));
146+}
147+for (const name of colormapNames) {
148+ colormapSelect.append(new Option(name, name));
149+}
150+
151+/* ------------------------------------------------------------- controls -- */
152+
153+/**
154+ * One slider per parameter the registry declares. A parameter that is a seed
155+ * gets a button instead: its value picks a draw and means nothing on its own,
156+ * so there is nothing to slide along.
157+ */
158+function buildSliders(
159+ host: HTMLElement,
160+ specs: ParamSpec[],
161+ values: Params,
162+ onChange: (key: string, value: number) => void,
163+): Map<string, (value: number) => void> {
164+ const setters = new Map<string, (value: number) => void>();
165+ host.textContent = '';
166+ for (const spec of specs) {
167+ const row = document.createElement('label');
168+ row.className = 'slider';
169+ if (spec.hint) row.title = spec.hint;
170+ const name = document.createElement('span');
171+ name.textContent = spec.label;
172+ row.append(name);
173+
174+ if (spec.reseed) {
175+ const button = document.createElement('button');
176+ button.type = 'button';
177+ button.textContent = 'new draw';
178+ const out = document.createElement('output');
179+ out.textContent = String(values[spec.key]);
180+ button.addEventListener('click', () => {
181+ const next = 1 + Math.floor(Math.random() * (spec.max - spec.min));
182+ out.textContent = String(next);
183+ onChange(spec.key, next);
184+ });
185+ row.append(button, out);
186+ } else {
187+ const input = document.createElement('input');
188+ input.type = 'range';
189+ input.min = String(spec.min);
190+ input.max = String(spec.max);
191+ input.step = String(spec.step);
192+ input.value = String(values[spec.key]);
193+ const out = document.createElement('output');
194+ out.textContent = fmtValue(values[spec.key]);
195+ input.addEventListener('input', () => {
196+ const v = Number(input.value);
197+ out.textContent = fmtValue(v);
198+ onChange(spec.key, v);
199+ });
200+ row.append(input, out);
201+ setters.set(spec.key, (value) => {
202+ input.value = String(value);
203+ out.textContent = fmtValue(value);
204+ });
205+ }
206+ host.append(row);
207+ }
208+ return setters;
209+}
210+
211+/** Scene edits go through the MATLAB interpreter, which is fast but not free;
212+ * a drag should not queue up one evaluation per pixel. */
213+let sceneTimer = 0;
214+const scheduleSceneUpdate = (): void => {
215+ clearTimeout(sceneTimer);
216+ sceneTimer = window.setTimeout(applyScene, 120);
217+};
218+
219+function applyScene(): void {
220+ if (!session) return;
221+ try {
222+ const before = session.dt;
223+ session.setScene(scene, sceneParams, sources.scene);
224+ clearError();
225+ // A different timestep leaves the leapfrog's two histories half a step
226+ // apart, so that is the one scene change the run cannot survive.
227+ if (Math.abs(session.dt - before) > 1e-12 * before) restart();
228+ showDt();
229+ dirty = true;
230+ view.setSource(
231+ session.gpu.stateBuffer(session.pressureName)!,
232+ session.gpu.stateBuffer('c')!,
233+ session.grid.nx,
234+ session.grid.ny,
235+ );
236+ } catch (e) {
237+ showError(e, sources.scene);
238+ }
239+}
240+
241+/** Sliders for the microphone's position, kept so a drag on the canvas can
242+ * move them too. */
243+let micSliders = new Map<string, (value: number) => void>();
244+
245+/** The microphone's position. Two sliders, in scene coordinates, since that
246+ * is what the scene's own parameters are in — and the canvas itself, which
247+ * is the direct way to place it. */
248+function buildMicControls(): void {
249+ const half = DOMAIN / 2;
250+ const specs: ParamSpec[] = [
251+ { key: 'x', label: 'mic x (m)', value: mic.x, min: -half, max: half, step: 0.05 },
252+ { key: 'y', label: 'mic y (m)', value: mic.y, min: -half, max: half, step: 0.05 },
253+ ];
254+ micSliders = buildSliders(el('micparams'), specs, { x: mic.x, y: mic.y }, (key, value) => {
255+ setMic(key === 'x' ? value : mic.x, key === 'y' ? value : mic.y, false);
256+ });
257+}
258+
259+/** Move the microphone, from wherever the instruction came. */
260+function setMic(x: number, y: number, syncSliders = true): void {
261+ const half = (session?.grid.L ?? DOMAIN) / 2;
262+ mic = {
263+ x: Math.max(-half, Math.min(half, x)),
264+ y: Math.max(-half, Math.min(half, y)),
265+ };
266+ session?.setMic(mic.x, mic.y);
267+ if (syncSliders) {
268+ micSliders.get('x')?.(mic.x);
269+ micSliders.get('y')?.(mic.y);
270+ }
271+ dirty = true;
272+}
273+
274+/**
275+ * Drag the microphone around the picture.
276+ *
277+ * The canvas shows the whole domain, so a pixel maps to a point with nothing
278+ * more than a scale: pointer down puts the microphone where it landed and
279+ * begins a drag, and pointer capture keeps that drag alive if it wanders off
280+ * the canvas.
281+ */
282+function micFromEvent(e: PointerEvent): { x: number; y: number } {
283+ const rect = canvas.getBoundingClientRect();
284+ const L = session?.grid.L ?? DOMAIN;
285+ const u = (e.clientX - rect.left) / rect.width;
286+ // Screen y runs down; the grid's runs up.
287+ const v = (e.clientY - rect.top) / rect.height;
288+ return { x: -L / 2 + u * L, y: L / 2 - v * L };
289+}
290+
291+canvas.addEventListener('pointerdown', (e) => {
292+ try {
293+ canvas.setPointerCapture(e.pointerId);
294+ } catch {
295+ // Best effort: a drag that leaves the canvas just stops tracking.
296+ }
297+ const at = micFromEvent(e);
298+ setMic(at.x, at.y);
299+ e.preventDefault();
300+});
301+canvas.addEventListener('pointermove', (e) => {
302+ if (!canvas.hasPointerCapture(e.pointerId)) return;
303+ const at = micFromEvent(e);
304+ setMic(at.x, at.y);
305+});
306+canvas.addEventListener('pointerup', (e) => {
307+ if (canvas.hasPointerCapture(e.pointerId)) canvas.releasePointerCapture(e.pointerId);
308+});
309+
310+function buildParamControls(): void {
311+ buildSliders(el('params'), model.params, params, (key, value) => {
312+ params = { ...params, [key]: value };
313+ session?.setParams(params);
314+ });
315+ el('scene-title').textContent = `scene — ${scene.blurb}`;
316+ buildSliders(el('sceneparams'), scene.params, sceneParams, (key, value) => {
317+ sceneParams = { ...sceneParams, [key]: value };
318+ scheduleSceneUpdate();
319+ });
320+}
321+
322+/* ------------------------------------------------------------- the loop -- */
323+
324+let frames = 0;
325+let lastFpsAt = performance.now();
326+let msPerFrame = 0;
327+/** A readback in flight; only one at a time, since they share a buffer. */
328+let reading = false;
329+
330+function restart(): void {
331+ if (!session) return;
332+ session.reset();
333+ peakSeen = 0;
334+ scale = autoScale ? 1e-6 : scale;
335+ dirty = true;
336+}
337+
338+/**
339+ * Compile the current sources into a new session and swap it in.
340+ *
341+ * The old session keeps running until the new one exists, and is only torn
342+ * down once the swap has happened. Two reasons: a compile that fails leaves
343+ * something on screen and something to edit rather than a dead page, and
344+ * destroying GPU resources while the next lot of shaders are still compiling
345+ * is exactly the sort of thing a browser is entitled to answer with a lost
346+ * device.
347+ */
348+let building = false;
349+async function rebuild(): Promise<void> {
350+ if (building) return;
351+ building = true;
352+ const old = session;
353+ try {
354+ const next = await ModelSession.create({
355+ device: device!,
356+ model,
357+ params,
358+ source: sources.model,
359+ scene,
360+ sceneParams,
361+ sceneSource: sources.scene,
362+ n: gridN,
363+ L: DOMAIN,
364+ cfl,
365+ });
366+ next.reset();
367+ session = next;
368+ view.setSource(
369+ next.gpu.stateBuffer(next.pressureName)!,
370+ next.gpu.stateBuffer('c')!,
371+ next.grid.nx,
372+ next.grid.ny,
373+ );
374+ next.setMic(mic.x, mic.y);
375+ old?.destroy();
376+ peakSeen = 0;
377+ scale = 1e-6;
378+ dirty = true;
379+ clearError();
380+ showDt();
381+ showRecording();
382+ el<HTMLButtonElement>('recompile').classList.remove('primary');
383+ el('compiled').textContent = describe(next);
384+ } catch (e) {
385+ // Which file the failure belongs to decides which one to show it against.
386+ const which = failingFile(e);
387+ showError(e, sources[which]);
388+ if (e instanceof ModelCompileError && e.start !== undefined) {
389+ editorFile.value = which;
390+ showFile(which);
391+ editor.select(e.start, e.end ?? e.start + 1);
392+ }
393+ } finally {
394+ building = false;
395+ }
396+}
397+
398+/** A scene failure is reported against `medium`, everything else against the
399+ * model's own functions. */
400+const failingFile = (e: unknown): 'model' | 'scene' =>
401+ e instanceof ModelCompileError && e.fn === 'medium' ? 'scene' : 'model';
402+
403+const describe = (s: ModelSession): string => {
404+ const { init, step } = s.describe();
405+ return [
406+ `% init — run once`,
407+ ...init.map((l) => ` ${l}`),
408+ ``,
409+ `% step — run ${stepsPerFrame}x per frame`,
410+ ...step.map((l) => ` ${l}`),
411+ ].join('\n');
412+};
413+
414+/** Where the microphone sits, in grid-index coordinates, for the marker. */
415+function micIndex(s: ModelSession): { ix: number; iy: number } {
416+ const { L, h } = s.grid;
417+ return { ix: (mic.x + L / 2) / h - 0.5, iy: (mic.y + L / 2) / h - 0.5 };
418+}
419+
420+/**
421+ * What the microphone has, and what it would sound like.
422+ *
423+ * The playback rate is just 1/dt — real time, real pitch, no translation —
424+ * because dt is already a real duration. The duration of the clip is still
425+ * worth watching: a pulse is a few dozen cycles, which at an audible pitch is
426+ * a few tens of milliseconds however long the simulation runs. Sustained
427+ * sound needs a sustained source (turn `continuous` up) and a long enough
428+ * run, which is what the speed control is for.
429+ */
430+function showRecording(): void {
431+ const info = el('recinfo');
432+ if (!session) {
433+ info.textContent = '';
434+ return;
435+ }
436+ const n = session.recorder.count;
437+ if (n === 0) {
438+ info.textContent = 'nothing recorded yet — press Run';
439+ return;
440+ }
441+ const plan = planPlayback(n, session.dt);
442+ info.textContent =
443+ `${n.toLocaleString()} samples · ${fmtTime(n * session.dt)} of simulated time · ` +
444+ `${fmtTime(plan.duration)} of audio at ${Math.round(plan.rate).toLocaleString()} Hz` +
445+ (plan.realTime ? '' : ' (sped up — outside the audio range)') +
446+ (session.recorder.full ? ' · buffer full' : '');
447+}
448+
449+/** The timestep the slider currently asks for, as a real duration. */
450+function showDt(): void {
451+ const out = el<HTMLOutputElement>('dtout');
452+ out.textContent = session ? fmtTime(session.dt) : '—';
453+ out.classList.toggle('unstable', cfl >= 1);
454+ el('cfl-label').title =
455+ `The timestep, as a fraction of the largest one this scheme is stable at ` +
456+ `on this grid` +
457+ (session ? ` (dt < ${fmtTime(session.dtLimit)} here)` : '') +
458+ `. At 1 and above it diverges, which is worth seeing once.`;
459+}
460+
461+/**
462+ * How well the grid resolves the current source frequency: the wavelength at
463+ * the background speed, divided by the cell size. Below `POOR_RESOLUTION`
464+ * cells per wavelength, what is on screen is as much grid dispersion as it is
465+ * sound — worth surfacing rather than leaving as a silent limitation.
466+ */
467+function resolutionInfo(s: ModelSession): { text: string; poor: boolean } {
468+ const f = params.f ?? 1;
469+ const wavelength = s.scene.cref / Math.max(f, 1e-9);
470+ const cells = wavelength / s.grid.h;
471+ const poor = cells < POOR_RESOLUTION;
472+ return { text: `λ = ${fmtLength(wavelength)} (${cells.toFixed(1)} cells)`, poor };
473+}
474+
475+function stats(): void {
476+ if (!session) return;
477+ const { grid } = session;
478+ // The rate is reported only when it means something. Paused, the loop takes
479+ // no steps and draws nothing, so any figure here would be measuring the
480+ // display's refresh interval and calling it solver throughput.
481+ const rate =
482+ running && msPerFrame > 0
483+ ? `${msPerFrame.toFixed(1)} ms/frame · ` +
484+ `${((stepsPerFrame * grid.npts) / (msPerFrame * 1e3)).toFixed(0)} Mpoint/s`
485+ : 'paused';
486+ const res = resolutionInfo(session);
487+ el('stats').innerHTML =
488+ `t = <b>${fmtTime(session.t)}</b> · step <b>${session.steps}</b> · ` +
489+ `dt = ${fmtTime(session.dt)} · ${grid.nx}×${grid.ny} over ${fmtLength(grid.L)} ` +
490+ `(${fmtLength(grid.h)} cells) · ` +
491+ `c ∈ [${fmtValue(session.scene.cmin)}, ${fmtValue(session.scene.cmax)}] m/s · ` +
492+ `<span${res.poor ? ' class="warn"' : ''}>${res.text}</span> · ${rate}`;
493+}
494+
495+/**
496+ * Follow the field with the colour scale.
497+ *
498+ * The only readback in the app, and it happens a few times a second rather
499+ * than every frame. Rising fast and falling slowly, so that a pulse arriving
500+ * is not clipped and a pulse leaving does not make the remaining ripples
501+ * flare up to full contrast.
502+ */
503+async function autoscaleStep(): Promise<void> {
504+ if (!session || reading || !autoScale) return;
505+ reading = true;
506+ try {
507+ const p = await session.read(session.pressureName);
508+ let peak = 0;
509+ for (const v of p) peak = Math.max(peak, Math.abs(v));
510+ peakSeen = Math.max(peakSeen, peak);
511+ scale = peak > scale ? peak : 0.97 * scale + 0.03 * peak;
512+ scale = Math.max(scale, 0.02 * peakSeen, 1e-9);
513+ dirty = true;
514+ } catch {
515+ // A rebuild can destroy the buffers mid-read; the next frame recovers.
516+ } finally {
517+ reading = false;
518+ }
519+}
520+
521+let sinceScale = 0;
522+function frame(): void {
523+ if (session && (running || dirty)) {
524+ if (running) session.step(stepsPerFrame);
525+ view.draw({
526+ scale,
527+ cref: session.scene.cref,
528+ cdev: session.scene.cdev,
529+ medium: showMedium ? 0.5 : 0,
530+ mic: micIndex(session),
531+ });
532+ dirty = false;
533+ // Only while running: a paused field cannot have changed since the last
534+ // time its peak was measured, so the readback would be pure waste.
535+ if (running && ++sinceScale >= 6) {
536+ sinceScale = 0;
537+ void autoscaleStep();
538+ }
539+ colorbar.setRange(-scale, scale);
540+ frames++;
541+ }
542+ const now = performance.now();
543+ if (now - lastFpsAt > 400) {
544+ // Frames that drew nothing are not frames; a paused page reports no rate
545+ // rather than the display's refresh interval dressed up as one.
546+ msPerFrame = frames > 0 ? (now - lastFpsAt) / frames : 0;
547+ frames = 0;
548+ lastFpsAt = now;
549+ stats();
550+ showRecording();
551+ }
552+ requestAnimationFrame(frame);
553+}
554+
555+/* -------------------------------------------------------------- wiring --- */
556+
557+function showFile(which: 'model' | 'scene'): void {
558+ editor.value = sources[which];
559+ el('editor-title').textContent =
560+ which === 'model'
561+ ? 'init and step, compiled to WebGPU'
562+ : 'medium(x, y, …) → sound speed and absorption, evaluated once on the CPU';
563+}
564+
565+modelSelect.addEventListener('change', () => {
566+ model = mModelByKey(modelSelect.value) ?? mModels[0];
567+ params = { ...defaultParams(model), ...params };
568+ sources.model = model.source;
569+ if (editorFile.value === 'model') showFile('model');
570+ el('blurb').textContent = `${model.blurb} ${scene.blurb}`;
571+ buildParamControls();
572+ void rebuild();
573+});
574+
575+sceneSelect.addEventListener('change', () => {
576+ scene = mSceneByKey(sceneSelect.value) ?? mScenes[0];
577+ sceneParams = defaultSceneParams(scene);
578+ // A room is no use with the source outside it, so a scene may ask for
579+ // source settings; they land in the sliders like any other value.
580+ if (scene.suggest) params = { ...params, ...scene.suggest };
581+ if (scene.mic) mic = { ...scene.mic };
582+ sources.scene = scene.source;
583+ if (editorFile.value === 'scene') showFile('scene');
584+ el('blurb').textContent = `${model.blurb} ${scene.blurb}`;
585+ buildParamControls();
586+ session?.setParams(params);
587+ applyScene();
588+ restart();
589+});
590+
591+el('gridsize').addEventListener('change', (e) => {
592+ gridN = Number((e.target as HTMLSelectElement).value);
593+ void rebuild();
594+});
595+
596+const listen = el<HTMLButtonElement>('listen');
597+listen.addEventListener('click', () => {
598+ if (!session) return;
599+ const s = session;
600+ const n = s.recorder.count;
601+ if (n === 0) {
602+ showError(
603+ new Error('the microphone has not recorded anything yet — press Run first'),
604+ '',
605+ );
606+ return;
607+ }
608+ listen.disabled = true;
609+ stopAudio();
610+ void s.recorder
611+ .read()
612+ .then((trace) => playTrace(trace, planPlayback(trace.length, s.dt)))
613+ .then(() => clearError())
614+ .catch((e: unknown) => showError(e, ''))
615+ .finally(() => {
616+ listen.disabled = false;
617+ });
618+});
619+
620+el('cfl').addEventListener('input', (e) => {
621+ cfl = Number((e.target as HTMLInputElement).value);
622+ session?.setCfl(cfl);
623+ showDt();
624+});
625+
626+el('spf').addEventListener('change', (e) => {
627+ stepsPerFrame = Number((e.target as HTMLSelectElement).value);
628+ if (session) el('compiled').textContent = describe(session);
629+});
630+
631+colormapSelect.addEventListener('change', () => {
632+ colormapName = colormapSelect.value;
633+ view.setColormap(colormaps[colormapName]);
634+ colorbar.setColormap(colormaps[colormapName]);
635+ dirty = true;
636+});
637+
638+el('scalemode').addEventListener('change', (e) => {
639+ autoScale = (e.target as HTMLSelectElement).value === 'auto';
640+ dirty = true;
641+});
642+
643+el('showmedium').addEventListener('change', (e) => {
644+ showMedium = (e.target as HTMLInputElement).checked;
645+ dirty = true;
646+});
647+
648+
649+const runPause = el<HTMLButtonElement>('runpause');
650+runPause.addEventListener('click', () => {
651+ running = !running;
652+ runPause.textContent = running ? 'Pause' : 'Run';
653+ // Start the frame-rate window fresh, so a resumed run is not averaged
654+ // against the paused frames before it.
655+ frames = 0;
656+ lastFpsAt = performance.now();
657+ dirty = true;
658+});
659+
660+el('restart').addEventListener('click', restart);
661+
662+editorFile.addEventListener('change', () => {
663+ showFile(editorFile.value as 'model' | 'scene');
664+});
665+
666+el('recompile').addEventListener('click', () => {
667+ void rebuild();
668+});
669+
670+el('revert').addEventListener('click', () => {
671+ const which = editorFile.value as 'model' | 'scene';
672+ sources[which] = which === 'model' ? model.source : scene.source;
673+ showFile(which);
674+ void rebuild();
675+});
676+
677+/* ---------------------------------------------------------------- start -- */
678+
679+// The grid dropdown's options are written statically in index.html; annotate
680+// each with the cell size it implies on the actual domain, rather than
681+// hardcoding a number that would drift if DOMAIN ever changed.
682+for (const opt of el<HTMLSelectElement>('gridsize').options) {
683+ const n = Number(opt.value);
684+ opt.textContent = `${n}² (${fmtLength(DOMAIN / n)} cells)`;
685+}
686+el('domaininfo').textContent =
687+ `Domain: ${fmtLength(DOMAIN)} × ${fmtLength(DOMAIN)}, background speed ` +
688+ `${C_AIR} m/s (air).`;
689+
690+modelSelect.value = model.key;
691+sceneSelect.value = scene.key;
692+colormapSelect.value = colormapName;
693+el('blurb').textContent = `${model.blurb} ${scene.blurb}`;
694+buildParamControls();
695+buildMicControls();
696+showFile('model');
697+el<HTMLInputElement>('cfl').value = String(cfl);
698+showDt();
699+showRecording();
700+await rebuild();
701+requestAnimationFrame(frame);
src/mgpu/compile.tsadded+205−0View file
@@ -0,0 +1,205 @@
1+/**
2+ * MATLAB source -> numbl's JIT IR, ready for the WGSL backend.
3+ *
4+ * A model file defines ordinary MATLAB functions; the host specializes the ones
5+ * it needs (`init`, `step`) for the concrete argument types of the current
6+ * grid. This is exactly how numbl drives its own JIT — the caller supplies
7+ * argument types, and lowering fixes every type and shape from there.
8+ *
9+ * Driving it through function signatures rather than injected scope means the
10+ * .m declares what it needs: each parameter name is matched against what the
11+ * host offers, and a name the host does not provide is a compile error rather
12+ * than a silently undefined variable.
13+ *
14+ * Two numbl passes matter here:
15+ * - `specializeUserFunction` lowers one function to IR, one statement per
16+ * operation (ANF), with every node's type fixed.
17+ * - `inlinePass` then folds single-use temps back into their consumer, so a
18+ * source line like `pn = 2*p - pm + cdt2 .* lap` becomes ONE statement whose
19+ * RHS is an expression tree — i.e. one GPU kernel instead of four.
20+ */
21+import { parseMFile } from 'numbl-src/numbl-core/parser/index.ts';
22+import { Workspace, Lowerer, tensorDouble, scalarDouble } from 'numbl-src/numbl-core/jit/index.ts';
23+import { specializeUserFunction } from 'numbl-src/numbl-core/jit/lowering/specialize.ts';
24+import { inlinePass } from 'numbl-src/numbl-core/jit/codegen/inlinePass.ts';
25+import type { IRFunc, IRStmt } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
26+import type { Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
27+import { fuseTemps } from './fuse.ts';
28+import { externalOpFiles, type GridSizes } from './externals.ts';
29+import { ModelCompileError } from './errors.ts';
30+
31+/** What the host can supply for an argument the .m declares. */
32+export type Binding =
33+ /** An array, passed in a GPU buffer. */
34+ | { kind: 'tensor'; shape: number[] }
35+ /** A tunable scalar. Deliberately carries no exact value: an exact scalar
36+ * would be constant-folded into the kernels, so moving a slider would force
37+ * a recompile instead of just rewriting a uniform. */
38+ | { kind: 'param' }
39+ /** A fixed scalar, exact so array constructors reading it keep static
40+ * shapes. */
41+ | { kind: 'const'; value: number };
42+
43+const typeOf = (b: Binding): Type => {
44+ switch (b.kind) {
45+ case 'tensor':
46+ return tensorDouble(b.shape);
47+ case 'param':
48+ return scalarDouble('unknown');
49+ case 'const':
50+ // Carry the sign too: numbl's sign lattice decides, for instance,
51+ // whether sqrt() of a value can go complex.
52+ return scalarDouble(
53+ b.value > 0 ? 'positive' : b.value < 0 ? 'negative' : 'zero',
54+ b.value,
55+ );
56+ }
57+};
58+
59+/** One specialized function, as the planner consumes it. */
60+export interface CompiledFunction {
61+ name: string;
62+ /** Declared arguments, in order, with the cName each lowered to. */
63+ params: { name: string; cName: string; binding: Binding }[];
64+ /** Requested outputs, in order, with the cName holding each result. */
65+ outputs: { name: string; cName: string; ty: Type }[];
66+ /** The lowered body. Read this only after `finish()`: the inline pass
67+ * REPLACES the statement array rather than mutating it, so this is a live
68+ * view of the function rather than a snapshot. */
69+ readonly body: IRStmt[];
70+}
71+
72+/** The shape of a `function` statement in numbl's AST. */
73+interface FunctionDecl {
74+ type: 'Function';
75+ name: string;
76+ params: string[];
77+ outputs: string[];
78+}
79+
80+/**
81+ * A parsed model. Specialize the functions you need, then call `finish()` once
82+ * — the inline pass rewrites every specialization together.
83+ */
84+export class CompiledModel {
85+ #lowerer: Lowerer;
86+ #decls: Map<string, FunctionDecl>;
87+ #bindings: Record<string, Binding>;
88+
89+ constructor(
90+ source: string,
91+ bindings: Record<string, Binding>,
92+ grid: GridSizes,
93+ fileName = 'model.m',
94+ ) {
95+ const ast = parseMFile(source, fileName);
96+ const ws = new Workspace(fileName, []);
97+ ws.addFile({ name: fileName, source, ast });
98+ // lap2 / lap4 become resolvable, with their type rules.
99+ for (const f of externalOpFiles(grid)) ws.addFile(f);
100+ ws.finalize();
101+
102+ this.#bindings = bindings;
103+ this.#lowerer = new Lowerer(ws);
104+ this.#decls = new Map();
105+ for (const stmt of ast.body as { type: string }[]) {
106+ if (stmt.type === 'Function') {
107+ const fn = stmt as unknown as FunctionDecl;
108+ this.#decls.set(fn.name, fn);
109+ }
110+ }
111+ }
112+
113+ /** Names of the functions the file defines. */
114+ functionNames(): string[] {
115+ return [...this.#decls.keys()];
116+ }
117+
118+ /**
119+ * Lower `name` for the current bindings, requesting `nargout` outputs.
120+ * Every declared parameter must name something the host provides.
121+ */
122+ specialize(name: string, nargout: number): CompiledFunction {
123+ const decl = this.#decls.get(name);
124+ if (!decl) {
125+ const defined = this.functionNames();
126+ throw new ModelCompileError(
127+ `the model must define a function named '${name}'` +
128+ (defined.length
129+ ? ` (it defines ${defined.map((n) => `'${n}'`).join(', ')})`
130+ : ' (it defines no functions)'),
131+ );
132+ }
133+ if (decl.outputs.length < nargout) {
134+ throw new ModelCompileError(
135+ `'${name}' must return ${nargout} value${nargout === 1 ? '' : 's'}, ` +
136+ `but declares ${decl.outputs.length}`,
137+ );
138+ }
139+
140+ const bindings = decl.params.map((p) => {
141+ const b = this.#bindings[p];
142+ if (!b) {
143+ const offered = Object.keys(this.#bindings).join(', ');
144+ throw new ModelCompileError(
145+ `'${name}' takes an argument named '${p}', which this app does not ` +
146+ `provide. Available: ${offered}.`,
147+ );
148+ }
149+ return b;
150+ });
151+
152+ const fn: IRFunc = specializeUserFunction.call(
153+ this.#lowerer,
154+ decl,
155+ bindings.map(typeOf),
156+ undefined,
157+ undefined,
158+ undefined,
159+ nargout,
160+ undefined,
161+ );
162+
163+ return {
164+ name,
165+ params: fn.params.map((p, i) => ({
166+ name: p,
167+ cName: fn.cParams[i],
168+ binding: bindings[i],
169+ })),
170+ outputs: fn.outputs.slice(0, nargout).map((o, i) => ({
171+ name: o,
172+ cName: fn.cOutputs[i],
173+ ty: fn.outputTypes[i],
174+ })),
175+ // A getter, not a snapshot: `finish()` runs after every specialization
176+ // and swaps in a rewritten statement array.
177+ get body() {
178+ return fn.body;
179+ },
180+ };
181+ }
182+
183+ /**
184+ * Run the fusion passes over everything specialized so far. They rewrite the
185+ * function bodies in place, so `CompiledFunction`s handed out earlier are
186+ * updated too.
187+ *
188+ * Two of them: numbl's, which folds the temps its own C backend would fuse,
189+ * and this project's (src/mgpu/fuse.ts), which folds the ones it declines —
190+ * `sin`, `exp`, `tanh` and the rest, which WGSL evaluates per element just
191+ * as happily as it does a multiply.
192+ *
193+ * Neither touches a variable the .m names. In particular a bare `pold = p;`
194+ * — the line that turns this step's field into the next step's history —
195+ * survives as its own statement, and plans as the copy it is: numbl's pass
196+ * gives every declared output a protective use count, and ours only folds
197+ * compiler temps.
198+ */
199+ finish(): void {
200+ inlinePass({ topLevelStmts: [], functions: this.#lowerer.specializations });
201+ for (const fn of this.#lowerer.specializations.values()) {
202+ fn.body = fuseTemps(fn.body);
203+ }
204+ }
205+}
src/mgpu/errors.tsadded+105−0View file
@@ -0,0 +1,105 @@
1+/**
2+ * Compile failures, reported in coordinates of the model file the user edits.
3+ *
4+ * Failures arrive from three places, each with its own idea of position:
5+ * numbl's parser (a `position` offset), numbl's lowerer (`UnsupportedConstruct`
6+ * / `JitTypeError`, with a `span`), and this project's WGSL emitter
7+ * (`UnsupportedOnGpu`, carrying the numbl span it was given). All of them are
8+ * offsets into the whole model file — the file is parsed once, and each function
9+ * is specialized from that one AST — so they need only be turned into a line and
10+ * column for the editor.
11+ */
12+
13+/** A compile failure located in the full model source. */
14+export class ModelCompileError extends Error {
15+ /** Offset into the whole .m file, when the failure has a position. */
16+ readonly start?: number;
17+ readonly end?: number;
18+ /** Name of the model function being compiled. */
19+ readonly fn?: string;
20+
21+ constructor(
22+ message: string,
23+ opts: { start?: number; end?: number; fn?: string; cause?: unknown } = {},
24+ ) {
25+ super(message, { cause: opts.cause });
26+ this.name = 'ModelCompileError';
27+ this.start = opts.start;
28+ this.end = opts.end;
29+ this.fn = opts.fn;
30+ }
31+}
32+
33+/** Extract whatever position information an error carries. */
34+function positionOf(e: unknown): { start?: number; end?: number } {
35+ const span = (e as { span?: { start?: unknown; end?: unknown } }).span;
36+ if (span && typeof span.start === 'number') {
37+ return {
38+ start: span.start,
39+ end: typeof span.end === 'number' ? span.end : undefined,
40+ };
41+ }
42+ // numbl's parser SyntaxError reports a bare offset.
43+ const position = (e as { position?: unknown }).position;
44+ if (typeof position === 'number') return { start: position };
45+ return {};
46+}
47+
48+/** Normalize any thrown value into a located `ModelCompileError`. */
49+function asCompileError(e: unknown, fn?: string): ModelCompileError {
50+ if (e instanceof ModelCompileError) return e;
51+ const { start, end } = positionOf(e);
52+ const raw = e instanceof Error ? e.message : String(e);
53+ // numbl's parse errors read as bare token complaints out of context.
54+ const message =
55+ (e as Error)?.name === 'SyntaxError' ? `MATLAB syntax error: ${raw}` : raw;
56+ return new ModelCompileError(message, { fn, start, end, cause: e });
57+}
58+
59+/**
60+ * Run `fn`, locating any compile failure in the model file. Use for whole-file
61+ * phases (parsing) that belong to no single function.
62+ */
63+export function inModel<T>(fn: () => T): T {
64+ try {
65+ return fn();
66+ } catch (e) {
67+ throw asCompileError(e);
68+ }
69+}
70+
71+/** Run `fn`, attributing any compile failure to the model function `name`. */
72+export function inFunction<T>(name: string, fn: () => T): T {
73+ try {
74+ return fn();
75+ } catch (e) {
76+ throw asCompileError(e, name);
77+ }
78+}
79+
80+/** Async form of `inFunction`. */
81+export async function inFunctionAsync<T>(
82+ name: string,
83+ fn: () => Promise<T>,
84+): Promise<T> {
85+ try {
86+ return await fn();
87+ } catch (e) {
88+ throw asCompileError(e, name);
89+ }
90+}
91+
92+/** Render a failure for display: message, section, and 1-based line/column. */
93+export function formatFailure(e: unknown, source: string): string {
94+ const message = e instanceof Error ? e.message : String(e);
95+ if (!(e instanceof ModelCompileError)) return message;
96+ const where: string[] = [];
97+ if (e.start !== undefined && e.start <= source.length) {
98+ const before = source.slice(0, e.start);
99+ const line = before.split('\n').length;
100+ const column = e.start - before.lastIndexOf('\n');
101+ where.push(`line ${line}, column ${column}`);
102+ }
103+ if (e.fn) where.push(`in ${e.fn}()`);
104+ return where.length ? `${message} (${where.join(', ')})` : message;
105+}
src/mgpu/externals.tsadded+86−0View file
@@ -0,0 +1,86 @@
1+/**
2+ * The Laplacian stencils, as external operations the .m can call: `lap2` (the
3+ * 5-point second-order stencil) and `lap4` (the 9-point fourth-order one).
4+ *
5+ * Everything else a model does is element-wise, so these are the only places
6+ * where a grid point reads its neighbours — the one thing the element-wise
7+ * kernel emitter cannot express, since it walks a single linear index across
8+ * every operand. Keeping them as named operations rather than as array slicing
9+ * (`p(2:end-1, :)` and friends) means the .m never has to know how the grid is
10+ * laid out in the buffer, and the host is free to implement the stencil as one
11+ * dispatch (src/mgpu/stencil.ts).
12+ *
13+ * numbl needs only their *type rule* in order to lower a call site. It gets
14+ * that from a `.mtoc2.js` workspace file — numbl's sanctioned extension point
15+ * for a JS-defined builtin (see `mtoc2UserFunctionsByName` in numbl's
16+ * LoweringContext). The file is evaluated in a bare CommonJS sandbox with no
17+ * imports available, so `transfer` builds numbl `Type` objects as plain
18+ * literals, and the grid size is baked in by the generator below (a grid change
19+ * recompiles anyway).
20+ *
21+ * The `emit`/`cBody` exports exist only because the loader's contract requires
22+ * them; we never emit C. The actual implementation is supplied by the WGSL
23+ * backend.
24+ */
25+
26+export interface GridSizes {
27+ /** Grid points, nx*ny. Grid fields are npts x 1 column vectors. */
28+ npts: number;
29+}
30+
31+const numericType = (rows: number, cols: number): string =>
32+ `{ kind: "Numeric", elem: "double", isComplex: false, ` +
33+ `dims: [${dim(rows)}, ${dim(cols)}], shape: [${rows}, ${cols}], sign: "unknown" }`;
34+
35+// numbl's tensorDouble() canonicalizes an extent of 1 to its shared DIM_ONE
36+// singleton; mirror that so types compare equal to host-built ones.
37+const dim = (n: number): string =>
38+ n === 1 ? `{ kind: "exact", value: 1 }` : `{ kind: "exact", value: ${n} }`;
39+
40+/** Source for one stencil's `.mtoc2.js`: grid field in, grid field out. */
41+function stencilSource(name: string, npts: number): string {
42+ return `
43+exports.name = ${JSON.stringify(name)};
44+
45+exports.transfer = function (argTypes, nargout) {
46+ if (argTypes.length !== 1) {
47+ throw new Error("${name} takes exactly one argument, got " + argTypes.length);
48+ }
49+ if (nargout > 1) {
50+ throw new Error("${name} returns one value, but " + nargout + " were requested");
51+ }
52+ var a = argTypes[0];
53+ if (!a || a.kind !== "Numeric" || a.isComplex) {
54+ throw new Error("${name} requires a real numeric array");
55+ }
56+ var s = a.shape;
57+ if (!s || s.length !== 2 || s[0] !== ${npts} || s[1] !== 1) {
58+ throw new Error(
59+ "${name} works on grid fields, so its argument must be ${npts}x1, not " +
60+ (s ? s.join("x") : "unknown shape")
61+ );
62+ }
63+ return [${numericType(npts, 1)}];
64+};
65+
66+// Never called: this project executes the IR on WebGPU and emits no C.
67+exports.emit = function () {
68+ throw new Error("${name}: no C backend (this stencil runs on WebGPU)");
69+};
70+exports.cBody = function () {
71+ return "";
72+};
73+`;
74+}
75+
76+/** Workspace files that make `lap2` / `lap4` resolvable during lowering. */
77+export function externalOpFiles(g: GridSizes): { name: string; source: string }[] {
78+ return [
79+ { name: 'lap2.mtoc2.js', source: stencilSource('lap2', g.npts) },
80+ { name: 'lap4.mtoc2.js', source: stencilSource('lap4', g.npts) },
81+ ];
82+}
83+
84+/** Names the WGSL backend must implement as stencil dispatches rather than
85+ * element-wise kernels. */
86+export const EXTERNAL_OPS = new Set(['lap2', 'lap4']);
src/mgpu/fuse.tsadded+136−0View file
@@ -0,0 +1,136 @@
1+/**
2+ * Fold the single-use ANF temps numbl's inline pass left behind.
3+ *
4+ * numbl's own pass (`inlinePass`) only folds producers its C backend can fuse,
5+ * which leaves out tensor-producing calls — `sin(x)`, `exp(x)`, `tanh(x)`. The
6+ * WGSL emitter fuses all of those happily, so without this pass a source line
7+ * like
8+ *
9+ * s = (amp * env) .* sin(2*pi*f*(t - t0)) .* exp(-(gx + gy));
10+ *
11+ * plans as half a dozen kernels rather than one, each writing a whole grid to
12+ * memory for the next one to read straight back.
13+ *
14+ * Same shape as numbl's pass, deliberately narrower where it matters: only
15+ * compiler temps (`_mtoc2_*`) are folded, so every variable the .m names keeps
16+ * its own buffer and one source line stays one kernel. The producer's RHS must
17+ * be something the emitter can evaluate per element, its result must be used
18+ * exactly once, and nothing between the two statements may write to anything
19+ * it reads.
20+ *
21+ * Adapted from math-webgpu-sandbox's src/mgpu/fuse.ts, trimmed to the
22+ * straight-line bodies this project compiles.
23+ */
24+import type { IRExpr, IRStmt } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
25+import { isGpuFusableExpr } from './wgsl.ts';
26+
27+const isTemp = (cName: string): boolean => cName.startsWith('_mtoc2_');
28+
29+export function fuseTemps(stmts: IRStmt[]): IRStmt[] {
30+ let cur = stmts;
31+ for (let iter = 0; iter < 32; iter++) {
32+ const next = fuseOnePass(cur);
33+ if (next === cur) break;
34+ cur = next;
35+ }
36+ return cur;
37+}
38+
39+/** One sweep, at most one fold. Returns the same array when nothing fired. */
40+function fuseOnePass(stmts: IRStmt[]): IRStmt[] {
41+ const uses = useCounts(stmts);
42+ for (let i = 0; i < stmts.length; i++) {
43+ const p = stmts[i];
44+ if (p.kind !== 'Assign' || !isTemp(p.cName)) continue;
45+ if (uses.get(p.cName) !== 1) continue;
46+ if (!isGpuFusableExpr(p.expr)) continue;
47+
48+ const reads = new Set<string>();
49+ walkVars(p.expr, (c) => reads.add(c));
50+
51+ for (let j = i + 1; j < stmts.length; j++) {
52+ const c = stmts[j];
53+ if (c.kind !== 'Assign') break; // anything else ends the safe window
54+ if (countIn(c.expr, p.cName) > 0) {
55+ // The one use. Fold into it if it is a kernel; if it is a stencil
56+ // call, or anything else that wants its argument in a buffer of its
57+ // own, leave the producer alone.
58+ if (countIn(c.expr, p.cName) === 1 && isGpuFusableExpr(c.expr)) {
59+ c.expr = substitute(c.expr, p.cName, p.expr);
60+ const out = stmts.slice();
61+ out.splice(i, 1);
62+ return out;
63+ }
64+ break;
65+ }
66+ // An intervening write to the temp itself, or to one of its operands,
67+ // invalidates the fold window. Statements that do neither are simply
68+ // skipped over — the stencil dispatch between the source term and the
69+ // update it feeds is exactly that case.
70+ if (c.cName === p.cName || reads.has(c.cName)) break;
71+ }
72+ }
73+ return stmts;
74+}
75+
76+function useCounts(stmts: IRStmt[]): Map<string, number> {
77+ const counts = new Map<string, number>();
78+ const bump = (c: string): void => {
79+ counts.set(c, (counts.get(c) ?? 0) + 1);
80+ };
81+ for (const s of stmts) if (s.kind === 'Assign') walkVars(s.expr, bump);
82+ return counts;
83+}
84+
85+function walkVars(e: IRExpr, visit: (cName: string) => void): void {
86+ const walk = (x: IRExpr): void => {
87+ switch (x.kind) {
88+ case 'Var':
89+ visit(x.cName);
90+ return;
91+ case 'Binary':
92+ walk(x.left);
93+ walk(x.right);
94+ return;
95+ case 'Unary':
96+ walk(x.operand);
97+ return;
98+ case 'Call':
99+ x.args.forEach(walk);
100+ return;
101+ default:
102+ return;
103+ }
104+ };
105+ walk(e);
106+}
107+
108+function countIn(e: IRExpr, cName: string): number {
109+ let n = 0;
110+ walkVars(e, (c) => {
111+ if (c === cName) n++;
112+ });
113+ return n;
114+}
115+
116+/** Replace the (single) `Var` read of `cName` with `replacement`. */
117+function substitute(e: IRExpr, cName: string, replacement: IRExpr): IRExpr {
118+ const sub = (x: IRExpr): IRExpr => {
119+ if (x.kind === 'Var' && x.cName === cName) return replacement;
120+ switch (x.kind) {
121+ case 'Binary':
122+ x.left = sub(x.left);
123+ x.right = sub(x.right);
124+ return x;
125+ case 'Unary':
126+ x.operand = sub(x.operand);
127+ return x;
128+ case 'Call':
129+ for (let i = 0; i < x.args.length; i++) x.args[i] = sub(x.args[i]);
130+ return x;
131+ default:
132+ return x;
133+ }
134+ };
135+ return sub(e);
136+}
src/mgpu/model.tsadded+290−0View file
@@ -0,0 +1,290 @@
1+/**
2+ * A .m model, compiled and running on the GPU.
3+ *
4+ * A model file is ordinary MATLAB: it defines an `init` function that builds
5+ * the initial state and a `step` function that advances it one timestep. Each
6+ * is specialized for the current grid and compiled into a ModelPlan, and both
7+ * operate on the same state buffers (see HostBuffers).
8+ *
9+ * Both functions return the state, in the same order, so their signatures say
10+ * exactly what they produce:
11+ *
12+ * function [p, pm, t] = init(npts)
13+ * function [pn, pold, tn] = step(p, pm, t, c, sig, x, y, dt, f, amp, ...)
14+ *
15+ * The host supplies the things that are setup rather than algorithm: the grid
16+ * coordinates, the medium the scene defines, the timestep, and the parameter
17+ * values. Each argument is matched to the .m's declared parameter name, so the
18+ * file documents its own interface.
19+ *
20+ * `t` is carried as a grid field rather than a scalar, and that is deliberate.
21+ * A batch of timesteps is one replay of a fixed op sequence, so nothing the
22+ * host writes between steps can change inside it — a clock uploaded per frame
23+ * would stand still for the whole batch. Making the model advance its own time
24+ * (`tn = t + dt`) keeps the source term correct however many steps are batched,
25+ * at the cost of one extra buffer and one extra kernel per step, which next to
26+ * the stencil is nothing.
27+ */
28+import { HostBuffers, ModelPlan } from './plan.ts';
29+import type { StencilPlan } from './stencil.ts';
30+import { inFunction, inFunctionAsync, inModel } from './errors.ts';
31+import { CompiledModel, type Binding } from './compile.ts';
32+
33+export interface ModelParams {
34+ [key: string]: number;
35+}
36+
37+/** The grid a model runs on, and the medium it runs in. */
38+export interface GridFields {
39+ nx: number;
40+ ny: number;
41+ /** Grid spacing, the same in x and y. */
42+ h: number;
43+ /** Coordinates of every grid point, npts each, x fastest. */
44+ x: Float32Array;
45+ y: Float32Array;
46+}
47+
48+/** What the scene defines, on the grid. */
49+export interface MediumFields {
50+ /** Sound speed, npts. */
51+ c: Float32Array;
52+ /** Absorption rate (the sponge and any absorbing scatterer), npts. */
53+ sig: Float32Array;
54+}
55+
56+export interface GpuModelOptions {
57+ device: GPUDevice;
58+ stencil: StencilPlan;
59+ grid: GridFields;
60+ medium: MediumFields;
61+ /** Model source (.m text). */
62+ source: string;
63+ /** Parameter names the .m may take as arguments. */
64+ paramNames: string[];
65+ /** State field names, in order (e.g. ['p', 'pm', 't']). */
66+ state: string[];
67+ /** Grid fields one kernel may read, overriding what the device allows.
68+ * Only for tests. */
69+ operandBudget?: number;
70+}
71+
72+/** Names the .m may take for the grid coordinates. */
73+export const GRID_NAMES = ['x', 'y'] as const;
74+/** Names the .m may take for the medium the scene defines. */
75+export const MEDIUM_NAMES = ['c', 'sig'] as const;
76+
77+export class GpuModel {
78+ readonly paramNames: string[];
79+ readonly state: string[];
80+ readonly npts: number;
81+
82+ #device: GPUDevice;
83+ #host: HostBuffers;
84+ #initPlan: ModelPlan;
85+ #stepPlan: ModelPlan;
86+ /** Timestep, host-owned: it follows from the medium and the grid (a CFL
87+ * condition), not from anything the user types, and it is folded into every
88+ * setParams so a .m that takes `dt` is never left with the zero a missing
89+ * parameter would default to. */
90+ #dt = 0;
91+ #readback: GPUBuffer;
92+ /** Which function wrote the state most recently; see `read`. */
93+ #lastRan: 'init' | 'step' = 'init';
94+
95+ private constructor(init: {
96+ device: GPUDevice;
97+ host: HostBuffers;
98+ initPlan: ModelPlan;
99+ stepPlan: ModelPlan;
100+ readback: GPUBuffer;
101+ paramNames: string[];
102+ state: string[];
103+ npts: number;
104+ }) {
105+ this.#device = init.device;
106+ this.#host = init.host;
107+ this.#initPlan = init.initPlan;
108+ this.#stepPlan = init.stepPlan;
109+ this.#readback = init.readback;
110+ this.paramNames = init.paramNames;
111+ this.state = init.state;
112+ this.npts = init.npts;
113+ }
114+
115+ static async create(opts: GpuModelOptions): Promise<GpuModel> {
116+ const { device, stencil, grid, medium, source, paramNames, state } = opts;
117+ const npts = grid.nx * grid.ny;
118+
119+ // What the .m may ask for by name. The grid geometry is exact, so a
120+ // constructor reading it (`zeros(npts, 1)`) keeps a static shape.
121+ const bindings: Record<string, Binding> = {
122+ npts: { kind: 'const', value: npts },
123+ nx: { kind: 'const', value: grid.nx },
124+ ny: { kind: 'const', value: grid.ny },
125+ h: { kind: 'const', value: grid.h },
126+ dt: { kind: 'param' },
127+ };
128+ for (const g of GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
129+ for (const m of MEDIUM_NAMES) bindings[m] = { kind: 'tensor', shape: [npts, 1] };
130+ for (const s of state) bindings[s] = { kind: 'tensor', shape: [npts, 1] };
131+ for (const p of paramNames) bindings[p] = { kind: 'param' };
132+
133+ // Parsing belongs to the file, not to either function.
134+ const compiled = inModel(() => new CompiledModel(source, bindings, { npts }));
135+ const nargout = state.length;
136+ const initFn = inFunction('init', () => compiled.specialize('init', nargout));
137+ const stepFn = inFunction('step', () => compiled.specialize('step', nargout));
138+ compiled.finish();
139+
140+ // Both functions return the state, in order, and both feed it back into
141+ // the shared buffers.
142+ const feedback = [...state];
143+
144+ const host = new HostBuffers(device);
145+ // The host owns the state and the inputs it uploads, whether or not a
146+ // given function happens to take them as arguments — `init` does not read
147+ // `p`, but it writes it, and `step` reads it back.
148+ for (const s of state) host.ensure(s, npts);
149+ for (const g of GRID_NAMES) host.ensure(g, npts);
150+ for (const m of MEDIUM_NAMES) host.ensure(m, npts);
151+
152+ const initPlan = await inFunctionAsync('init', () =>
153+ ModelPlan.create(device, stencil, { fn: initFn, feedback }, host, opts.operandBudget),
154+ );
155+ const stepPlan = await inFunctionAsync('step', () =>
156+ ModelPlan.create(device, stencil, { fn: stepFn, feedback }, host, opts.operandBudget),
157+ );
158+
159+ host.upload('x', grid.x);
160+ host.upload('y', grid.y);
161+ host.upload('c', medium.c);
162+ host.upload('sig', medium.sig);
163+
164+ const readback = device.createBuffer({
165+ label: 'mgpu-readback',
166+ size: 4 * npts,
167+ usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
168+ });
169+
170+ return new GpuModel({
171+ device, host, initPlan, stepPlan, readback, paramNames, state, npts,
172+ });
173+ }
174+
175+ /** The timestep in force. Host-owned; see `#dt`. */
176+ get dt(): number {
177+ return this.#dt;
178+ }
179+
180+ setDt(dt: number): void {
181+ this.#dt = dt;
182+ }
183+
184+ setParams(params: ModelParams): void {
185+ const merged = { dt: this.#dt, ...params };
186+ this.#initPlan.setParams(merged);
187+ this.#stepPlan.setParams(merged);
188+ }
189+
190+ /**
191+ * Swap the medium under a running model. It is data, not code — its shape in
192+ * the bindings depends only on the grid — so changing the scene is two
193+ * buffer writes and needs no recompile.
194+ */
195+ uploadMedium(medium: MediumFields): void {
196+ this.#host.upload('c', medium.c);
197+ this.#host.upload('sig', medium.sig);
198+ }
199+
200+ /** Write a host-owned value directly. Lets a test set up an exact initial
201+ * condition instead of going through `init`. */
202+ upload(name: string, data: Float32Array): void {
203+ this.#host.upload(name, data);
204+ }
205+
206+ /** Run `init`, replacing the state. */
207+ init(): void {
208+ const enc = this.#device.createCommandEncoder({ label: 'mgpu-init' });
209+ this.#initPlan.encodeSteps(enc, 1);
210+ this.#device.queue.submit([enc.finish()]);
211+ this.#lastRan = 'init';
212+ }
213+
214+ /**
215+ * Advance `steps` timesteps. Synchronous — this only records commands and
216+ * submits them; nothing is read back and nothing is awaited.
217+ *
218+ * `after` is recorded once per step, so anything that must see every
219+ * timestep (the microphone) rides along in the same submission.
220+ */
221+ step(steps = 1, after?: (encoder: GPUCommandEncoder) => void): void {
222+ const enc = this.#device.createCommandEncoder({ label: 'mgpu-step' });
223+ this.#stepPlan.encodeSteps(enc, steps, after);
224+ this.#device.queue.submit([enc.finish()]);
225+ this.#lastRan = 'step';
226+ }
227+
228+ /**
229+ * The buffer currently holding a named value. A field the .m computes is
230+ * produced by both functions, into separate buffers (only the state is
231+ * shared), so this resolves to whichever function ran most recently — which
232+ * is what makes the first frame show the initial state rather than an
233+ * unwritten buffer.
234+ */
235+ #locate(name: string): { buffer: GPUBuffer; count: number } | null {
236+ const [first, second] =
237+ this.#lastRan === 'init'
238+ ? [this.#initPlan, this.#stepPlan]
239+ : [this.#stepPlan, this.#initPlan];
240+ const buffer = first.buffer(name) ?? second.buffer(name);
241+ const count = first.elementCount(name) ?? second.elementCount(name);
242+ if (!buffer || count === undefined) return null;
243+ return { buffer, count };
244+ }
245+
246+ /** The GPU buffer a named value would be read from right now. */
247+ valueBuffer(name: string): GPUBuffer | null {
248+ return this.#locate(name)?.buffer ?? null;
249+ }
250+
251+ /**
252+ * The buffer a host-owned field lives in — the state between calls, or an
253+ * input like the sound speed.
254+ *
255+ * This is what the renderer binds, and it must be this rather than
256+ * `valueBuffer`: a bind group is built once and holds a particular buffer,
257+ * while `init` and `step` write their outputs into buffers of their own and
258+ * only agree here, where their feedback copies land. Binding either
259+ * function's private buffer would draw a stale field for half the run.
260+ */
261+ stateBuffer(name: string): GPUBuffer | null {
262+ return this.#host.get(name)?.buffer ?? null;
263+ }
264+
265+ /** Read a named value back to the CPU. The only await in the whole loop. */
266+ async read(name: string): Promise<Float32Array> {
267+ const located = this.#locate(name);
268+ if (!located) throw new Error(`read: the model has no value named '${name}'`);
269+ const { buffer, count } = located;
270+ const enc = this.#device.createCommandEncoder({ label: `mgpu-read-${name}` });
271+ enc.copyBufferToBuffer(buffer, 0, this.#readback, 0, 4 * count);
272+ this.#device.queue.submit([enc.finish()]);
273+ await this.#readback.mapAsync(GPUMapMode.READ, 0, 4 * count);
274+ const out = new Float32Array(this.#readback.getMappedRange(0, 4 * count).slice(0));
275+ this.#readback.unmap();
276+ return out;
277+ }
278+
279+ /** What the .m compiled to, for display. */
280+ describe(): { init: string[]; step: string[] } {
281+ return { init: this.#initPlan.describe(), step: this.#stepPlan.describe() };
282+ }
283+
284+ destroy(): void {
285+ this.#initPlan.destroy();
286+ this.#stepPlan.destroy();
287+ this.#host.destroy();
288+ this.#readback.destroy();
289+ }
290+}
src/mgpu/numbl.d.tsadded+366−0View file
@@ -0,0 +1,366 @@
1+/**
2+ * The numbl compiler surface this project depends on.
3+ *
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+ * scenes — 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.
10+ *
11+ * Declaring the surface here rather than type-checking numbl's sources
12+ * directly keeps this project's compiler settings independent of numbl's, and
13+ * pins the exact contract we rely on. If numbl changes one of these shapes,
14+ * the build breaks here with a clear diff rather than deep inside its tree.
15+ *
16+ * Only the nodes the WGSL backend actually walks are spelled out; every other
17+ * IR kind is collapsed into a catch-all so that unhandled constructs are
18+ * rejected with a message instead of being silently mis-compiled.
19+ */
20+
21+declare module 'numbl-src/numbl-core/jit/lowering/types.ts' {
22+ export type Sign =
23+ | 'positive' | 'nonneg' | 'negative' | 'nonpositive'
24+ | 'zero' | 'nonzero' | 'unknown';
25+
26+ export type DimInfo = { kind: 'exact'; value: number } | { kind: 'unknown' };
27+
28+ export type NumericExact =
29+ | number
30+ | Float64Array
31+ | { re: number; im: number }
32+ | { re: Float64Array; im: Float64Array };
33+
34+ export interface NumericType {
35+ kind: 'Numeric';
36+ elem: 'double' | 'logical' | 'char' | string;
37+ isComplex: boolean;
38+ dims: DimInfo[];
39+ /** Present iff every dim is exact. */
40+ shape?: number[];
41+ sign: Sign;
42+ exact?: NumericExact;
43+ }
44+
45+ /** Everything the WGSL backend rejects. */
46+ export interface NonNumericType {
47+ kind: 'Void' | 'Unknown' | 'String' | 'Handle' | 'Struct' | 'Class' | 'Cell';
48+ }
49+
50+ export type Type = NumericType | NonNumericType;
51+
52+ export function isMultiElement(t: NumericType): boolean;
53+ export function tensorDouble(shape: number[], exact?: Float64Array): NumericType;
54+ export function scalarDouble(sign?: Sign, exact?: number): NumericType;
55+}
56+
57+declare module 'numbl-src/numbl-core/jit/lowering/ir.ts' {
58+ import type { Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
59+
60+ export interface Span {
61+ file: string;
62+ start: number;
63+ end: number;
64+ }
65+
66+ export interface NumLit {
67+ kind: 'NumLit';
68+ value: number;
69+ ty: Type;
70+ span: Span;
71+ }
72+ export interface Var {
73+ kind: 'Var';
74+ name: string;
75+ cName: string;
76+ ty: Type;
77+ span: Span;
78+ }
79+ export interface Binary {
80+ kind: 'Binary';
81+ builtin: string;
82+ left: IRExpr;
83+ right: IRExpr;
84+ ty: Type;
85+ span: Span;
86+ }
87+ export interface Unary {
88+ kind: 'Unary';
89+ builtin: string;
90+ operand: IRExpr;
91+ ty: Type;
92+ span: Span;
93+ }
94+ export interface Call {
95+ kind: 'Call';
96+ cName: string;
97+ name: string;
98+ args: IRExpr[];
99+ ty: Type;
100+ span: Span;
101+ }
102+ /** Any other IR expression kind — rejected by the WGSL emitter. */
103+ export interface OtherExpr {
104+ kind:
105+ | 'ImagLit' | 'StringLit' | 'TensorBuild' | 'TensorConcat' | 'CellLit'
106+ | 'CellEmpty' | 'CellIndexLoad' | 'HandleLit' | 'HandleCaptureLoad'
107+ | 'StructLit' | 'MemberLoad' | 'IndexLoad' | 'IndexSlice' | 'EndRef'
108+ | 'MakeRange';
109+ ty: Type;
110+ span: Span;
111+ }
112+
113+ export type IRExpr = NumLit | Var | Binary | Unary | Call | OtherExpr;
114+
115+ export interface Assign {
116+ kind: 'Assign';
117+ name: string;
118+ cName: string;
119+ ty: Type;
120+ expr: IRExpr;
121+ span: Span;
122+ }
123+ /**
124+ * A counted loop. The planner unrolls it, so only the fields that decide
125+ * the trip count and the loop variable's value are spelled out. `step` is
126+ * already a literal number in the IR — numbl rejects a non-literal step
127+ * during lowering — while `start` and `end` are expressions that must carry
128+ * an exact value for the planner to accept the loop.
129+ */
130+ export interface For {
131+ kind: 'For';
132+ /** Loop variable, as written in the .m. */
133+ varName: string;
134+ /** Loop variable's cName, the key the planner binds its value under. */
135+ cVar: string;
136+ start: IRExpr;
137+ step: number;
138+ end: IRExpr;
139+ body: IRStmt[];
140+ span: Span;
141+ }
142+ /**
143+ * Multi-output call statement: `[a, b] = f(x, y)`. For `isBuiltin: true`
144+ * the builtin's `transfer(argTypes, nargout)` typed the slots during
145+ * lowering; args arrive ANF'd. The planner accepts this only for the
146+ * batched transforms (`synth`/`analys`), where output k is the transform
147+ * of argument k.
148+ */
149+ export interface MultiAssignCall {
150+ kind: 'MultiAssignCall';
151+ cName: string;
152+ name: string;
153+ isBuiltin?: boolean;
154+ args: IRExpr[];
155+ outputs: ReadonlyArray<{
156+ ty: Type;
157+ binding: { name: string; cName: string } | null;
158+ }>;
159+ span: Span;
160+ }
161+
162+ /** Any other IR statement kind — rejected by the planner. */
163+ export interface OtherStmt {
164+ kind:
165+ | 'ExprStmt' | 'If' | 'While' | 'ReturnFromFunction' | 'Break'
166+ | 'Continue' | 'TypeComment' | 'MemberStore'
167+ | 'IndexStore' | 'IndexSliceStore' | 'CellIndexStore';
168+ span: Span;
169+ }
170+
171+ export type IRStmt = Assign | For | MultiAssignCall | OtherStmt;
172+
173+ export interface IRFunc {
174+ name: string;
175+ cName: string;
176+ /** Parameter source names. */
177+ params: string[];
178+ /** Parameter cNames, parallel to `params`. */
179+ cParams: string[];
180+ paramTypes: Type[];
181+ /** Output source names. */
182+ outputs: string[];
183+ /** Output cNames, parallel to `outputs`. */
184+ cOutputs: string[];
185+ outputTypes: Type[];
186+ body: IRStmt[];
187+ span: Span;
188+ }
189+
190+ export interface IRProgram {
191+ topLevelStmts: IRStmt[];
192+ functions: Map<string, IRFunc>;
193+ }
194+}
195+
196+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+
222+ export interface AbstractSyntaxTree {
223+ body: Stmt[];
224+ }
225+ export function parseMFile(input: string, fileName?: string): AbstractSyntaxTree;
226+ export class SyntaxError extends Error {}
227+}
228+
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+
287+declare module 'numbl-src/numbl-core/jit/index.ts' {
288+ import type { AbstractSyntaxTree } from 'numbl-src/numbl-core/parser/index.ts';
289+ import type { IRProgram, IRFunc, Span } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
290+ import type { Type, NumericType, Sign } from 'numbl-src/numbl-core/jit/lowering/types.ts';
291+
292+ export interface WorkspaceFile {
293+ name: string;
294+ source: string;
295+ ast?: AbstractSyntaxTree;
296+ }
297+
298+ export class Workspace {
299+ constructor(mainFile: string, searchPaths?: ReadonlyArray<string>);
300+ addFile(file: WorkspaceFile): void;
301+ finalize(): void;
302+ }
303+
304+ export interface EnvEntry {
305+ cName: string;
306+ ty: Type;
307+ maybeUnassigned?: boolean;
308+ }
309+
310+ export class Lowerer {
311+ constructor(workspace: Workspace);
312+ /** Pre-bindable variable scope: seed host-provided values here. */
313+ env: Map<string, EnvEntry>;
314+ specializations: Map<string, IRFunc>;
315+ lowerProgram(ast: AbstractSyntaxTree): IRProgram;
316+ }
317+
318+ /** Thrown for MATLAB the JIT pipeline cannot lower; carries a source span. */
319+ export class UnsupportedConstruct extends Error {
320+ span?: Span;
321+ }
322+ export class JitTypeError extends Error {
323+ span?: Span;
324+ }
325+
326+ export function tensorDouble(shape: number[], exact?: Float64Array): NumericType;
327+ export function scalarDouble(sign?: Sign, exact?: number): NumericType;
328+ export function isMultiElement(t: NumericType): boolean;
329+}
330+
331+declare module 'numbl-src/numbl-core/jit/lowering/specialize.ts' {
332+ import type { Lowerer } from 'numbl-src/numbl-core/jit/index.ts';
333+ import type { IRFunc, IRExpr, Span } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
334+ import type { Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
335+
336+ /**
337+ * Lower one user function for a concrete argument-type signature. Called with
338+ * a `Lowerer` as `this` (numbl's own JIT does the same), so specializations
339+ * accumulate in `lowerer.specializations`.
340+ */
341+ export function specializeUserFunction(
342+ this: Lowerer,
343+ decl: unknown,
344+ argTypes: Type[],
345+ specSource?: string,
346+ definingFile?: string,
347+ preSeedOutput?: { name: string; ty: Type; initExpr: IRExpr },
348+ nargout?: number,
349+ callSiteSpan?: Span,
350+ ): IRFunc;
351+}
352+
353+declare module 'numbl-src/numbl-core/jit/codegen/inlinePass.ts' {
354+ import type { IRProgram } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
355+ /** Folds single-use ANF temps into their consumer, in place. */
356+ export function inlinePass(prog: IRProgram): void;
357+}
358+
359+declare module 'numbl-src/numbl-core/jit/builtins/index.ts' {
360+ export interface Builtin {
361+ name: string;
362+ /** Safe to evaluate one output element from one input element per slot. */
363+ elementwise?: boolean;
364+ }
365+ export function getBuiltin(name: string): Builtin | undefined;
366+}
src/mgpu/plan.tsadded+742−0View file
@@ -0,0 +1,742 @@
1+/**
2+ * Statement list -> a replayable sequence of GPU operations.
3+ *
4+ * Everything expensive happens once, here: pipeline compilation, buffer
5+ * allocation, bind-group construction. Because numbl fixes every type and
6+ * shape at lowering time, the resulting op sequence is fully static — so
7+ * `encodeSteps` is pure synchronous command recording, with no allocation, no
8+ * pipeline lookup and no readback. That is what lets a whole batch of
9+ * timesteps be encoded into one submit and keeps the CPU out of the loop.
10+ */
11+import { isMultiElement } from 'numbl-src/numbl-core/jit/lowering/types.ts';
12+import type { Assign, IRExpr, IRStmt } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
13+import type { NumericType, Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
14+import type { CompiledFunction } from './compile.ts';
15+import { EXTERNAL_OPS } from './externals.ts';
16+import { StencilPlan, type StencilKind } from './stencil.ts';
17+import { kernelOperandBudget } from '../device.ts';
18+import {
19+ buildKernel,
20+ checkShapes,
21+ UnsupportedOnGpu,
22+ WORKGROUP_SIZE,
23+ type KernelInputs,
24+} from './wgsl.ts';
25+
26+const isNumeric = (t: Type): t is NumericType => t.kind === 'Numeric';
27+const isTensor = (t: Type): boolean => isNumeric(t) && isMultiElement(t);
28+const numel = (t: NumericType): number => (t.shape ?? []).reduce((a, b) => a * b, 1);
29+
30+interface Slot {
31+ buffer: GPUBuffer;
32+ count: number;
33+}
34+
35+const makeBuffer = (device: GPUDevice, label: string, count: number): GPUBuffer =>
36+ device.createBuffer({
37+ label,
38+ size: 4 * count,
39+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST,
40+ });
41+
42+/**
43+ * Buffers for host-bound variables, shared across plans.
44+ *
45+ * A model is two programs — `init` and `step` — compiled separately but
46+ * operating on the same state. `p` in the step must be the very buffer `init`
47+ * wrote, so the buffers for host bindings live here rather than inside either
48+ * plan.
49+ */
50+export class HostBuffers {
51+ #device: GPUDevice;
52+ #slots = new Map<string, Slot>();
53+
54+ constructor(device: GPUDevice) {
55+ this.#device = device;
56+ }
57+
58+ ensure(name: string, count: number): Slot {
59+ const existing = this.#slots.get(name);
60+ if (existing) {
61+ if (existing.count !== count) {
62+ throw new UnsupportedOnGpu(
63+ `'${name}' is ${existing.count} elements in one program and ` +
64+ `${count} in another`,
65+ );
66+ }
67+ return existing;
68+ }
69+ const slot = { buffer: makeBuffer(this.#device, `mgpu-${name}`, count), count };
70+ this.#slots.set(name, slot);
71+ return slot;
72+ }
73+
74+ get(name: string): Slot | undefined {
75+ return this.#slots.get(name);
76+ }
77+
78+ /** Upload initial data for a host binding. */
79+ upload(name: string, data: Float32Array): void {
80+ const slot = this.#slots.get(name);
81+ if (!slot) throw new Error(`upload: no buffer named '${name}'`);
82+ if (data.length !== slot.count) {
83+ throw new Error(
84+ `upload '${name}': expected ${slot.count} elements, got ${data.length}`,
85+ );
86+ }
87+ this.#device.queue.writeBuffer(slot.buffer, 0, data as Float32Array<ArrayBuffer>);
88+ }
89+
90+ destroy(): void {
91+ for (const s of this.#slots.values()) s.buffer.destroy();
92+ this.#slots.clear();
93+ }
94+}
95+
96+type Op =
97+ | {
98+ kind: 'kernel';
99+ pipeline: GPUComputePipeline;
100+ bindGroup: GPUBindGroup;
101+ count: number;
102+ label: string;
103+ /** Set when the kernel had to write to scratch because its output
104+ * aliases one of its inputs; copied back after the dispatch. */
105+ copyBack?: { from: GPUBuffer; to: GPUBuffer; bytes: number };
106+ }
107+ | {
108+ kind: 'stencil';
109+ pipeline: GPUComputePipeline;
110+ bindGroup: GPUBindGroup;
111+ workgroups: number;
112+ label: string;
113+ }
114+ | { kind: 'copy'; from: GPUBuffer; to: GPUBuffer; bytes: number; label: string };
115+
116+export interface PlanSpec {
117+ /** The specialized function this plan executes. */
118+ fn: CompiledFunction;
119+ /** Output index -> host binding name to copy the result into after the run,
120+ * so the next call reads it (the new field feeds the old). */
121+ feedback: (string | null)[];
122+}
123+
124+/**
125+ * Bind group layout for a kernel: the output at 0, `inputs` read-only storage
126+ * buffers after it, then the params buffer.
127+ *
128+ * Declared explicitly rather than with `layout: 'auto'`, because an auto layout
129+ * only contains the bindings the shader actually references — so a kernel that
130+ * happens to use no parameters would drop the params binding and no longer
131+ * match the bind group. An explicit layout may carry bindings the shader
132+ * ignores.
133+ */
134+function kernelLayout(device: GPUDevice, inputs: number): GPUBindGroupLayout {
135+ const readOnly = (binding: number): GPUBindGroupLayoutEntry => ({
136+ binding,
137+ visibility: GPUShaderStage.COMPUTE,
138+ buffer: { type: 'read-only-storage' },
139+ });
140+ return device.createBindGroupLayout({
141+ entries: [
142+ {
143+ binding: 0,
144+ visibility: GPUShaderStage.COMPUTE,
145+ buffer: { type: 'storage' },
146+ },
147+ ...Array.from({ length: inputs }, (_, i) => readOnly(i + 1)),
148+ readOnly(inputs + 1),
149+ ],
150+ });
151+}
152+
153+/**
154+ * Compile one shader into a pipeline.
155+ *
156+ * No validation error scope around it: `createComputePipelineAsync` already
157+ * rejects on a shader that will not compile or a layout that does not match,
158+ * which is the whole of what a scope here would have caught, and the scope
159+ * costs an extra device round trip per pipeline. `getCompilationInfo`, which
160+ * has the line and column within the generated WGSL, is asked for only once
161+ * something has gone wrong, and defensively even then.
162+ *
163+ * In practice the WGSL here is generated, so a shader that fails to compile is
164+ * this project's bug rather than the user's; a mistake in a .m is caught
165+ * earlier, by the emitter, with a position in the MATLAB source.
166+ */
167+async function makePipeline(
168+ device: GPUDevice,
169+ code: string,
170+ label: string,
171+ bindGroupLayout: GPUBindGroupLayout,
172+): Promise<GPUComputePipeline> {
173+ const module = device.createShaderModule({ code, label });
174+ try {
175+ return await device.createComputePipelineAsync({
176+ layout: device.createPipelineLayout({ bindGroupLayouts: [bindGroupLayout] }),
177+ compute: { module, entryPoint: 'main' },
178+ label,
179+ });
180+ } catch (e) {
181+ throw new UnsupportedOnGpu(
182+ `generated WGSL failed to compile for '${label}':\n` +
183+ `${await shaderErrors(module, e)}\n--- shader ---\n${code}`,
184+ );
185+ }
186+}
187+
188+/** Per-line compile errors, if the browser will hand them over. */
189+async function shaderErrors(module: GPUShaderModule, cause: unknown): Promise<string> {
190+ const fallback = cause instanceof Error ? cause.message : String(cause);
191+ try {
192+ const info = await module.getCompilationInfo();
193+ const errors = info.messages.filter((m) => m.type === 'error');
194+ if (!errors.length) return fallback;
195+ return errors.map((e) => ` ${e.lineNum}:${e.linePos} ${e.message}`).join('\n');
196+ } catch {
197+ return fallback;
198+ }
199+}
200+
201+/** A compiled .m function, ready to run on the GPU. */
202+export class ModelPlan {
203+ /** Scalar parameter names, in the order the params buffer expects them. */
204+ readonly paramNames: string[];
205+
206+ #device: GPUDevice;
207+ #ops: Op[];
208+ #owned: GPUBuffer[];
209+ #paramBuf: GPUBuffer;
210+ #paramData: Float32Array;
211+ /** Public name -> buffer, for uploading initial state and reading results. */
212+ #byName: Map<string, Slot>;
213+
214+ private constructor(init: {
215+ device: GPUDevice;
216+ ops: Op[];
217+ byName: Map<string, Slot>;
218+ owned: GPUBuffer[];
219+ paramBuf: GPUBuffer;
220+ paramData: Float32Array;
221+ paramNames: string[];
222+ }) {
223+ this.#device = init.device;
224+ this.#ops = init.ops;
225+ this.#byName = init.byName;
226+ this.#owned = init.owned;
227+ this.#paramBuf = init.paramBuf;
228+ this.#paramData = init.paramData;
229+ this.paramNames = init.paramNames;
230+ }
231+
232+ static async create(
233+ device: GPUDevice,
234+ stencil: StencilPlan,
235+ spec: PlanSpec,
236+ host: HostBuffers,
237+ /** Overrides what the device allows; only tests pass it. */
238+ operandBudget?: number,
239+ ): Promise<ModelPlan> {
240+ const { fn } = spec;
241+
242+ const slots = new Map<string, Slot>();
243+ const byName = new Map<string, Slot>();
244+ const owned: GPUBuffer[] = [];
245+ /** Scalars the .m computes from its parameters, by cName. */
246+ const derivedScalars = new Map<string, { name: string; expr: IRExpr }>();
247+ /** Grid fields one kernel may read on this device (see fitToBudget). */
248+ const budget = operandBudget ?? kernelOperandBudget(device);
249+ /** How many kernels a line has been split into, for naming the pieces. */
250+ let splits = 0;
251+
252+ const alloc = (label: string, count: number): Slot => {
253+ const buffer = makeBuffer(device, label, count);
254+ owned.push(buffer);
255+ return { buffer, count };
256+ };
257+
258+ // Arguments, bound by what the function's signature declares. Array
259+ // arguments come from the shared pool, so a value one function returns is
260+ // the same buffer the next one reads. Scalar parameters share one small
261+ // storage buffer, in signature order.
262+ const paramNames: string[] = [];
263+ const paramSlots = new Map<string, number>();
264+ for (const p of fn.params) {
265+ if (p.binding.kind === 'tensor') {
266+ const count = p.binding.shape.reduce((x, y) => x * y, 1);
267+ const slot = host.ensure(p.name, count);
268+ slots.set(p.cName, slot);
269+ byName.set(p.name, slot);
270+ } else if (p.binding.kind === 'param') {
271+ paramSlots.set(p.cName, paramNames.length);
272+ paramNames.push(p.name);
273+ }
274+ // `const` arguments are exact in the IR and fold into the kernels.
275+ }
276+ const paramData = new Float32Array(Math.max(1, paramNames.length));
277+ const paramBuf = device.createBuffer({
278+ label: 'mgpu-params',
279+ size: 4 * paramData.length,
280+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
281+ });
282+
283+ const ops: Op[] = [];
284+ for (const stmt of fn.body) {
285+ await planStatement(stmt);
286+ }
287+
288+ planFeedback();
289+
290+ return new ModelPlan({ device, ops, byName, owned, paramBuf, paramData, paramNames });
291+
292+ async function planStatement(stmt: IRStmt): Promise<void> {
293+ if (stmt.kind === 'ReturnFromFunction') return; // nothing follows it
294+ if (stmt.kind !== 'Assign') {
295+ throw new UnsupportedOnGpu(
296+ `a model function body may only contain assignments ` +
297+ `(found '${stmt.kind}')`,
298+ stmt.span,
299+ );
300+ }
301+ if (!isNumeric(stmt.ty)) {
302+ throw new UnsupportedOnGpu(`'${stmt.name}' is not a numeric value`, stmt.span);
303+ }
304+ if (!isTensor(stmt.ty)) {
305+ // A scalar the model derives from its parameters (`om = 2*pi*f`). It
306+ // gets no buffer and no dispatch: the kernels that read it bind it as
307+ // a `let` in their prologue.
308+ derivedScalars.set(stmt.cName, { name: stmt.name, expr: stmt.expr });
309+ return;
310+ }
311+ const count = numel(stmt.ty);
312+
313+ // Reuse the destination buffer across steps: the same cName always maps
314+ // to the same buffer, so a step allocates nothing.
315+ let dest = slots.get(stmt.cName);
316+ if (!dest) {
317+ dest = alloc(`mgpu-${stmt.name}`, count);
318+ slots.set(stmt.cName, dest);
319+ } else if (dest.count !== count) {
320+ throw new UnsupportedOnGpu(
321+ `'${stmt.name}' changes size between assignments`,
322+ stmt.span,
323+ );
324+ }
325+ byName.set(stmt.name, dest);
326+
327+ const ext = externalCall(stmt);
328+ if (ext) return planStencil(stmt, ext, dest);
329+
330+ // Element-wise. Checked against the whole line first, so a broadcasting
331+ // mistake is reported against what was written rather than against a
332+ // fragment of it.
333+ checkShapes(stmt.expr, stmt.ty, stmt.name);
334+ const expr = await fitToBudget(stmt.expr, stmt.name, count, stmt.span);
335+ await emitElementwise(stmt.name, expr, stmt.ty, stmt.span, dest, count, stmt.cName);
336+ }
337+
338+ /**
339+ * Emit one element-wise kernel: `expr` evaluated at every index into
340+ * `dest`. `selfCName` is the variable being assigned, if any, so an
341+ * in-place update can be spotted.
342+ */
343+ async function emitElementwise(
344+ name: string,
345+ expr: IRExpr,
346+ ty: NumericType,
347+ span: unknown,
348+ dest: Slot,
349+ count: number,
350+ selfCName?: string,
351+ ): Promise<void> {
352+ // Collect the distinct tensor operands and give them dense binding slots.
353+ const tensors = new Map<string, number>();
354+ collectTensorVars(expr, (cName) => {
355+ if (!tensors.has(cName)) tensors.set(cName, tensors.size);
356+ });
357+
358+ const label = `${name} = <${count} elements, element-wise>`;
359+ const kernel = buildKernel(
360+ { kind: 'Assign', name, ty, expr, span } as unknown as Assign,
361+ {
362+ tensors,
363+ params: paramSlots,
364+ scalars: derivedScalars,
365+ } satisfies KernelInputs,
366+ count,
367+ label,
368+ );
369+
370+ const bindGroupLayout = kernelLayout(device, tensors.size);
371+ const pipeline = await makePipeline(device, kernel.code, label, bindGroupLayout);
372+
373+ // WebGPU forbids aliasing a writable storage binding with another
374+ // binding in the same group, so an in-place update (`p = p + 1`) writes
375+ // to scratch and copies back. Element-wise kernels only ever touch
376+ // their own index, so the copy is the only cost.
377+ const aliased = selfCName !== undefined && tensors.has(selfCName);
378+ const target = aliased ? alloc(`mgpu-${name}-scratch`, count) : dest;
379+
380+ const entries: GPUBindGroupEntry[] = [
381+ { binding: 0, resource: { buffer: target.buffer } },
382+ ];
383+ for (const [cName, i] of tensors) {
384+ const s = slots.get(cName);
385+ if (!s) {
386+ throw new UnsupportedOnGpu(`'${name}' reads a value with no buffer`, span);
387+ }
388+ entries.push({ binding: i + 1, resource: { buffer: s.buffer } });
389+ }
390+ entries.push({ binding: tensors.size + 1, resource: { buffer: paramBuf } });
391+
392+ ops.push({
393+ kind: 'kernel',
394+ pipeline,
395+ bindGroup: device.createBindGroup({ layout: bindGroupLayout, entries }),
396+ count,
397+ label,
398+ copyBack: aliased
399+ ? { from: target.buffer, to: dest.buffer, bytes: 4 * count }
400+ : undefined,
401+ });
402+ }
403+
404+ /**
405+ * Split an expression that reads more grid fields than one kernel may bind.
406+ *
407+ * A kernel binds one storage buffer per distinct field it reads, plus its
408+ * output and the parameter block, and WebGPU guarantees only eight per
409+ * compute stage — fewer in compatibility mode. numbl's inline pass, which
410+ * is what makes one source line become one kernel, does not know about
411+ * that limit, and a model has no way to ask it for less: a temporary used
412+ * once is exactly what it folds away.
413+ *
414+ * So the budget is enforced here instead. Any child subtree that reads
415+ * more than one field is evaluated into its own buffer and replaced by a
416+ * reference to it, which leaves the parent reading at most one field per
417+ * child. The result is the same arithmetic in a few more passes over
418+ * memory, and it only happens on a line that would not otherwise compile.
419+ */
420+ async function fitToBudget(
421+ expr: IRExpr,
422+ hint: string,
423+ count: number,
424+ span: unknown,
425+ ): Promise<IRExpr> {
426+ if (tensorCount(expr) <= budget) return expr;
427+
428+ const fit = async (e: IRExpr): Promise<IRExpr> => {
429+ if (tensorCount(e) <= budget) return e;
430+ const kids = children(e);
431+ if (!kids.length) return e;
432+ const out: IRExpr[] = [];
433+ for (const kid of kids) {
434+ const fitted = await fit(kid);
435+ out.push(tensorCount(fitted) > 1 ? await hoist(fitted) : fitted);
436+ }
437+ return withChildren(e, out);
438+ };
439+
440+ /** Evaluate a subtree into its own buffer and hand back a reference. */
441+ const hoist = async (e: IRExpr): Promise<IRExpr> => {
442+ if (!isNumeric(e.ty) || !isTensor(e.ty)) return e;
443+ const name = `${hint}_part${++splits}`;
444+ const cName = `mgpu_split_${splits}`;
445+ const slot = alloc(`mgpu-${name}`, count);
446+ slots.set(cName, slot);
447+ await emitElementwise(name, e, e.ty, e.span, slot, count);
448+ return { kind: 'Var', name, cName, ty: e.ty, span: e.span } as IRExpr;
449+ };
450+
451+ const fitted = await fit(expr);
452+ if (tensorCount(fitted) > budget) {
453+ throw new UnsupportedOnGpu(
454+ `'${hint}' reads ${tensorCount(fitted)} grid fields at once, and this ` +
455+ `device allows ${budget} per kernel. Compute part of it into a ` +
456+ `named field on a line of its own.`,
457+ span,
458+ );
459+ }
460+ return fitted;
461+ }
462+
463+ /** `lp = lap2(p)`: one stencil dispatch, src and dst distinct. */
464+ async function planStencil(
465+ stmt: Assign,
466+ ext: { name: StencilKind; arg: IRExpr & { kind: 'Var' } },
467+ dest: Slot,
468+ ): Promise<void> {
469+ const argSlot = slots.get(ext.arg.cName);
470+ if (!argSlot) {
471+ throw new UnsupportedOnGpu(
472+ `'${ext.name}' reads '${ext.arg.name}', which has no buffer`,
473+ stmt.span,
474+ );
475+ }
476+ // A stencil reads its neighbours, so unlike an element-wise kernel it
477+ // cannot be routed through scratch and copied back — the neighbours
478+ // would already have been overwritten. WebGPU forbids the aliasing
479+ // outright anyway; refuse rather than silently reroute.
480+ if (argSlot.buffer === dest.buffer) {
481+ throw new UnsupportedOnGpu(
482+ `'${stmt.name} = ${ext.name}(${ext.arg.name})' reads and writes the ` +
483+ `same buffer; assign to a new name instead`,
484+ stmt.span,
485+ );
486+ }
487+ if (dest.count !== stencil.npts) {
488+ throw new UnsupportedOnGpu(
489+ `'${ext.name}' produces a grid field (${stencil.npts} points), but ` +
490+ `'${stmt.name}' holds ${dest.count}`,
491+ stmt.span,
492+ );
493+ }
494+ ops.push({
495+ kind: 'stencil',
496+ pipeline: await stencil.pipeline(ext.name),
497+ bindGroup: stencil.createBinding(argSlot.buffer, dest.buffer),
498+ workgroups: stencil.workgroups,
499+ label: `${stmt.name} = ${ext.name}(${ext.arg.name})`,
500+ });
501+ }
502+
503+ /**
504+ * Feed declared outputs back into the argument buffers they replace, so
505+ * the next call reads what this one produced.
506+ *
507+ * The copies are not independent: a model whose new history field is the
508+ * old current one (`function [pn, pold] = step(p, pm, ...)`) has an output
509+ * whose *source* is another output's *destination*. Doing them in order
510+ * would then copy the new value where the old one was wanted, silently.
511+ * So any source that a previous copy overwrites is staged through scratch
512+ * first — normally none, since a model that writes `pold = p;` gets its
513+ * own buffer from the copy kernel that line plans to.
514+ */
515+ function planFeedback(): void {
516+ const copies: { from: Slot; to: Slot; label: string }[] = [];
517+ fn.outputs.forEach((out, i) => {
518+ const to = spec.feedback[i];
519+ if (!to) return;
520+ const src = slots.get(out.cName);
521+ const dst = host.get(to);
522+ if (!src) {
523+ throw new UnsupportedOnGpu(
524+ `'${fn.name}' declares the output '${out.name}' but never assigns it`,
525+ );
526+ }
527+ if (!dst) throw new UnsupportedOnGpu(`'${to}' is not a host binding`);
528+ if (src.count !== dst.count) {
529+ throw new UnsupportedOnGpu(
530+ `'${out.name}' (${src.count} elements) cannot feed '${to}' (${dst.count})`,
531+ );
532+ }
533+ copies.push({ from: src, to: dst, label: `${out.name} -> ${to}` });
534+ });
535+
536+ const written = new Set<GPUBuffer>();
537+ for (const c of copies) written.add(c.to.buffer);
538+ for (const c of copies) {
539+ // Only a source another copy overwrites needs staging, and only if it
540+ // is not that same copy's own destination (which is a no-op anyway).
541+ if (c.from.buffer !== c.to.buffer && written.has(c.from.buffer)) {
542+ const scratch = alloc(`mgpu-feedback-scratch`, c.from.count);
543+ ops.push({
544+ kind: 'copy',
545+ from: c.from.buffer,
546+ to: scratch.buffer,
547+ bytes: 4 * c.from.count,
548+ label: `${c.label} (staged)`,
549+ });
550+ c.from = scratch;
551+ }
552+ }
553+ for (const c of copies) {
554+ if (c.from.buffer === c.to.buffer) continue; // already in place
555+ ops.push({
556+ kind: 'copy',
557+ from: c.from.buffer,
558+ to: c.to.buffer,
559+ bytes: 4 * c.from.count,
560+ label: c.label,
561+ });
562+ }
563+ }
564+ }
565+
566+ /** Upload parameter values, in `paramNames` order. Cheap — call freely. */
567+ setParams(values: Record<string, number>): void {
568+ this.paramNames.forEach((name, i) => {
569+ const v = values[name];
570+ this.#paramData[i] = Number.isFinite(v) ? v : 0;
571+ });
572+ this.#device.queue.writeBuffer(
573+ this.#paramBuf,
574+ 0,
575+ this.#paramData as Float32Array<ArrayBuffer>,
576+ );
577+ }
578+
579+ /** Buffer holding the named value, or undefined if the .m never binds it. */
580+ buffer(name: string): GPUBuffer | undefined {
581+ return this.#byName.get(name)?.buffer;
582+ }
583+
584+ elementCount(name: string): number | undefined {
585+ return this.#byName.get(name)?.count;
586+ }
587+
588+ /**
589+ * Record `steps` passes of this plan. Synchronous: no awaits, no readback.
590+ * The dispatches share one compute pass, which WebGPU executes in submission
591+ * order with a barrier between them.
592+ *
593+ * `after` runs once per step, inside the same submission — which is what
594+ * lets the microphone sample every timestep rather than every frame.
595+ */
596+ encodeSteps(
597+ encoder: GPUCommandEncoder,
598+ steps: number,
599+ after?: (encoder: GPUCommandEncoder) => void,
600+ ): void {
601+ for (let s = 0; s < steps; s++) {
602+ this.#encodeOps(encoder);
603+ after?.(encoder);
604+ }
605+ }
606+
607+ /** Record one pass over the op sequence into `encoder`. */
608+ #encodeOps(encoder: GPUCommandEncoder): void {
609+ let pass: GPUComputePassEncoder | null = null;
610+ const inPass = (): GPUComputePassEncoder => {
611+ if (!pass) pass = encoder.beginComputePass({ label: 'mgpu-step' });
612+ return pass;
613+ };
614+ const endPass = (): void => {
615+ if (pass) {
616+ pass.end();
617+ pass = null;
618+ }
619+ };
620+ for (const op of this.#ops) {
621+ switch (op.kind) {
622+ case 'kernel': {
623+ const p = inPass();
624+ p.setPipeline(op.pipeline);
625+ p.setBindGroup(0, op.bindGroup);
626+ p.dispatchWorkgroups(Math.ceil(op.count / WORKGROUP_SIZE));
627+ if (op.copyBack) {
628+ endPass();
629+ encoder.copyBufferToBuffer(
630+ op.copyBack.from, 0, op.copyBack.to, 0, op.copyBack.bytes,
631+ );
632+ }
633+ break;
634+ }
635+ case 'stencil': {
636+ const p = inPass();
637+ p.setPipeline(op.pipeline);
638+ p.setBindGroup(0, op.bindGroup);
639+ p.dispatchWorkgroups(op.workgroups);
640+ break;
641+ }
642+ case 'copy':
643+ endPass();
644+ encoder.copyBufferToBuffer(op.from, 0, op.to, 0, op.bytes);
645+ break;
646+ }
647+ }
648+ endPass();
649+ }
650+
651+ /** Human-readable op sequence — what the .m actually compiled to. */
652+ describe(): string[] {
653+ return this.#ops.map((op) => `${op.kind.padEnd(7)} ${op.label}`);
654+ }
655+
656+ destroy(): void {
657+ for (const b of this.#owned) b.destroy();
658+ this.#paramBuf.destroy();
659+ this.#owned.length = 0;
660+ }
661+}
662+
663+/** `lp = lap2(p)` -> the stencil's name and its argument. */
664+function externalCall(
665+ stmt: Assign,
666+): { name: StencilKind; arg: IRExpr & { kind: 'Var' } } | null {
667+ const e = stmt.expr;
668+ if (e.kind !== 'Call' || !EXTERNAL_OPS.has(e.name)) return null;
669+ if (e.args.length !== 1) {
670+ throw new UnsupportedOnGpu(
671+ `'${e.name}' must be applied to a single variable`,
672+ stmt.span,
673+ );
674+ }
675+ const arg = e.args[0];
676+ if (arg.kind !== 'Var') {
677+ throw new UnsupportedOnGpu(
678+ `'${e.name}' must be applied to a variable, not an expression — ` +
679+ `name the field first`,
680+ stmt.span,
681+ );
682+ }
683+ return { name: e.name as StencilKind, arg };
684+}
685+
686+/** Distinct grid fields an expression reads — its storage-buffer cost. */
687+function tensorCount(e: IRExpr): number {
688+ const seen = new Set<string>();
689+ collectTensorVars(e, (c) => seen.add(c));
690+ return seen.size;
691+}
692+
693+/** The subexpressions of a node, in evaluation order. Leaves have none. */
694+function children(e: IRExpr): IRExpr[] {
695+ switch (e.kind) {
696+ case 'Binary':
697+ return [e.left, e.right];
698+ case 'Unary':
699+ return [e.operand];
700+ case 'Call':
701+ return e.args;
702+ default:
703+ return [];
704+ }
705+}
706+
707+/** The same node with its subexpressions replaced. */
708+function withChildren(e: IRExpr, kids: IRExpr[]): IRExpr {
709+ switch (e.kind) {
710+ case 'Binary':
711+ return { ...e, left: kids[0], right: kids[1] };
712+ case 'Unary':
713+ return { ...e, operand: kids[0] };
714+ case 'Call':
715+ return { ...e, args: kids };
716+ default:
717+ return e;
718+ }
719+}
720+
721+function collectTensorVars(e: IRExpr, visit: (cName: string) => void): void {
722+ const walk = (x: IRExpr): void => {
723+ switch (x.kind) {
724+ case 'Var':
725+ if (isTensor(x.ty)) visit(x.cName);
726+ return;
727+ case 'Binary':
728+ walk(x.left);
729+ walk(x.right);
730+ return;
731+ case 'Unary':
732+ walk(x.operand);
733+ return;
734+ case 'Call':
735+ x.args.forEach(walk);
736+ return;
737+ default:
738+ return;
739+ }
740+ };
741+ walk(e);
742+}
src/mgpu/registry.tsadded+153−0View file
@@ -0,0 +1,153 @@
1+/**
2+ * The available models: their MATLAB source, and the metadata the host owns.
3+ *
4+ * A model's *algorithm* lives in its .m file. Everything around it lives here:
5+ * the parameter names the .m may take as arguments, their defaults and slider
6+ * ranges, the state fields it advances, and which of them is the pressure the
7+ * app draws. The .m declares nothing about these — it just names the
8+ * parameters it wants, and `CompiledModel` matches each against this table.
9+ */
10+import leapfrogSource from '../../models/leapfrog.m?raw';
11+import leapfrog4Source from '../../models/leapfrog4.m?raw';
12+import { C_AIR } from '../units.ts';
13+
14+export type Params = Record<string, number>;
15+
16+/** A tunable scalar the .m may take as an argument. */
17+export interface ParamSpec {
18+ key: string;
19+ label: string;
20+ value: number;
21+ min: number;
22+ max: number;
23+ step: number;
24+ /** This parameter is a random seed: its value picks a draw and means
25+ * nothing on its own, so the UI offers a button that jumps to another one
26+ * rather than a slider. */
27+ reseed?: boolean;
28+ /** Shown as a tooltip. */
29+ hint?: string;
30+}
31+
32+export interface MModel {
33+ key: string;
34+ label: string;
35+ blurb: string;
36+ /** State field names the .m advances, in the order its functions return
37+ * them. The first is the pressure, which is what gets drawn. */
38+ state: string[];
39+ params: ParamSpec[];
40+ /** Order of the spatial stencil the .m uses — what the stable timestep
41+ * follows from. */
42+ order: 2 | 4;
43+ /** MATLAB source — the algorithm itself. */
44+ source: string;
45+}
46+
47+/**
48+ * The source parameters, shared by every model here: both differ only in
49+ * their Laplacian.
50+ */
51+const sourceParams: ParamSpec[] = [
52+ {
53+ key: 'f',
54+ label: 'frequency (Hz)',
55+ value: 300,
56+ min: 50,
57+ max: 2000,
58+ step: 10,
59+ hint: `Cycles per second. At the background speed of ${C_AIR} m/s, this is ` +
60+ 'a real audible pitch, and the wavelength (speed / frequency) is what ' +
61+ 'the grid has to resolve — the app warns when it does not.',
62+ },
63+ {
64+ key: 'tw',
65+ label: 'pulse width (s)',
66+ value: 0.01,
67+ min: 0.001,
68+ max: 0.05,
69+ step: 0.001,
70+ hint: 'Duration of the Gaussian envelope, in seconds. Wide means many cycles and a narrow band; narrow means a click.',
71+ },
72+ {
73+ key: 'cw',
74+ label: 'continuous',
75+ value: 0,
76+ min: 0,
77+ max: 1,
78+ step: 0.05,
79+ hint: '0 is a single pulse, 1 a wave that turns on smoothly and stays on. In between is both.',
80+ },
81+ {
82+ key: 'point',
83+ label: 'point source',
84+ value: 0,
85+ min: 0,
86+ max: 1,
87+ step: 0.05,
88+ hint: '0 is a line source spanning the grid, whose far field is a plane wave; 1 is a point at (source x, source y).',
89+ },
90+ {
91+ key: 'x0',
92+ label: 'source x (m)',
93+ value: -3,
94+ min: -4.5,
95+ max: 4.5,
96+ step: 0.1,
97+ },
98+ {
99+ key: 'y0',
100+ label: 'source y (m)',
101+ value: 0,
102+ min: -4.5,
103+ max: 4.5,
104+ step: 0.1,
105+ hint: 'Only meaningful for a point source; a line source spans y.',
106+ },
107+ {
108+ key: 'w',
109+ label: 'source width (m)',
110+ value: 0.05,
111+ min: 0.01,
112+ max: 0.3,
113+ step: 0.01,
114+ hint: 'Physical size of the source. Much smaller than a wavelength and the grid cannot resolve it.',
115+ },
116+ {
117+ key: 't0',
118+ label: 'start time (s)',
119+ value: 0.05,
120+ min: 0.005,
121+ max: 0.2,
122+ step: 0.005,
123+ hint: 'When the pulse is emitted, in seconds. Needs a few pulse widths of head room, or it starts already part-way up its envelope.',
124+ },
125+];
126+
127+const leapfrog: MModel = {
128+ key: 'leapfrog',
129+ label: 'Leapfrog, 5-point Laplacian',
130+ blurb: 'Second order in space and time. The plain scheme.',
131+ state: ['p', 'pm', 't'],
132+ params: sourceParams,
133+ order: 2,
134+ source: leapfrogSource,
135+};
136+
137+const leapfrog4: MModel = {
138+ key: 'leapfrog4',
139+ label: 'Leapfrog, 9-point Laplacian',
140+ blurb: 'Fourth order in space: much less grid dispersion at the same resolution.',
141+ state: ['p', 'pm', 't'],
142+ params: sourceParams,
143+ order: 4,
144+ source: leapfrog4Source,
145+};
146+
147+export const mModels: MModel[] = [leapfrog, leapfrog4];
148+
149+export const mModelByKey = (key: string): MModel | undefined =>
150+ mModels.find((m) => m.key === key);
151+
152+export const defaultParams = (m: { params: ParamSpec[] }): Params =>
153+ Object.fromEntries(m.params.map((p) => [p.key, p.value]));
src/mgpu/session.tsadded+242−0View file
@@ -0,0 +1,242 @@
1+/**
2+ * One running simulation: grid, medium, compiled .m, timestep.
3+ *
4+ * Everything that is not rendering. The app and the tests both go through
5+ * this, so there is one place that decides how a pair of .m files becomes
6+ * something running on the GPU — and nothing about it is browser-specific
7+ * beyond needing a GPUDevice.
8+ */
9+import { makeGrid, stableDt, type Grid } from '../grid.ts';
10+import { StencilPlan } from './stencil.ts';
11+import { GpuModel, type ModelParams } from './model.ts';
12+import { Scene } from '../scene/scene.ts';
13+import { Recorder } from '../audio/recorder.ts';
14+import type { MModel, Params } from './registry.ts';
15+import type { MScene } from '../scene/registry.ts';
16+
17+export interface ModelSessionOptions {
18+ device: GPUDevice;
19+ model: MModel;
20+ params: Params;
21+ /** Override the model source — the editor's working copy. */
22+ source?: string;
23+ scene: MScene;
24+ sceneParams: Params;
25+ /** Override the scene source — the editor's working copy. */
26+ sceneSource?: string;
27+ /** Grid points per side. */
28+ n: number;
29+ /** Side length of the square domain, in metres. Defaults to a small toy
30+ * domain (2 m); the app passes its real DOMAIN constant explicitly. */
31+ L?: number;
32+ /** Fraction of the stability limit to take as the timestep. */
33+ cfl?: number;
34+ /** Grid fields one kernel may read, overriding what the device allows.
35+ * Only for tests, which use it to exercise the planner's kernel splitting
36+ * on a device that would never need it. */
37+ operandBudget?: number;
38+}
39+
40+export class ModelSession {
41+ readonly device: GPUDevice;
42+ readonly model: MModel;
43+ readonly grid: Grid;
44+ readonly gpu: GpuModel;
45+ /** The microphone: the pressure at one point, every timestep. */
46+ readonly recorder: Recorder;
47+
48+ /** Model time and step count since the last reset. */
49+ t = 0;
50+ steps = 0;
51+
52+ #scene: Scene;
53+ #sceneModel: MScene;
54+ #params: Params;
55+ #cfl: number;
56+
57+ private constructor(init: {
58+ device: GPUDevice;
59+ model: MModel;
60+ grid: Grid;
61+ gpu: GpuModel;
62+ recorder: Recorder;
63+ scene: Scene;
64+ sceneModel: MScene;
65+ params: Params;
66+ cfl: number;
67+ }) {
68+ this.device = init.device;
69+ this.model = init.model;
70+ this.grid = init.grid;
71+ this.gpu = init.gpu;
72+ this.recorder = init.recorder;
73+ this.#scene = init.scene;
74+ this.#sceneModel = init.sceneModel;
75+ this.#params = init.params;
76+ this.#cfl = init.cfl;
77+ }
78+
79+ static async create(opts: ModelSessionOptions): Promise<ModelSession> {
80+ const { device, model, params, scene, sceneParams } = opts;
81+ const grid = makeGrid(opts.n, opts.L ?? 2);
82+ const cfl = opts.cfl ?? 0.5;
83+
84+ // The scene first: the model's timestep depends on the fastest speed in
85+ // it, and the medium is an argument the compiled step reads.
86+ const built = Scene.create({
87+ grid,
88+ source: opts.sceneSource ?? scene.source,
89+ paramNames: scene.params.map((p) => p.key),
90+ params: sceneParams,
91+ });
92+
93+ const stencil = StencilPlan.create(device, { nx: grid.nx, ny: grid.ny, h: grid.h });
94+ const gpu = await GpuModel.create({
95+ device,
96+ stencil,
97+ grid,
98+ medium: built,
99+ source: opts.source ?? model.source,
100+ paramNames: model.params.map((p) => p.key),
101+ state: model.state,
102+ operandBudget: opts.operandBudget,
103+ });
104+
105+ // The microphone listens to the host-owned pressure buffer, which is
106+ // where both `init` and `step` leave their result.
107+ const recorder = await Recorder.create({
108+ device,
109+ field: gpu.stateBuffer(model.state[0])!,
110+ nx: grid.nx,
111+ ny: grid.ny,
112+ });
113+
114+ const session = new ModelSession({
115+ device, model, grid, gpu, recorder, scene: built, sceneModel: scene, params, cfl,
116+ });
117+ session.#applyDt();
118+ gpu.setParams(params);
119+ return session;
120+ }
121+
122+ get scene(): Scene {
123+ return this.#scene;
124+ }
125+
126+ get sceneModel(): MScene {
127+ return this.#sceneModel;
128+ }
129+
130+ /** The pressure field's name — the first state field the model advances. */
131+ get pressureName(): string {
132+ return this.model.state[0];
133+ }
134+
135+ get dt(): number {
136+ return this.gpu.dt;
137+ }
138+
139+ /** Fraction of the stability limit the timestep is taken at. */
140+ get cfl(): number {
141+ return this.#cfl;
142+ }
143+
144+ /**
145+ * Change the timestep, as a fraction of what the scheme is stable at.
146+ *
147+ * A fraction rather than a number of seconds, because the stable step
148+ * follows from the grid spacing and the fastest speed in the scene: refine
149+ * the grid or drop in a faster scatterer and a dt that was fine becomes a
150+ * dt that diverges. This way the meaning of the setting survives both.
151+ *
152+ * Applied to the running simulation as it stands. Leapfrog carries two
153+ * fields a timestep apart, so changing dt between steps is inconsistent by
154+ * the size of the change; a small drag is a small transient, and a large
155+ * jump is worth a Restart.
156+ */
157+ setCfl(cfl: number): void {
158+ this.#cfl = cfl;
159+ this.#applyDt();
160+ this.gpu.setParams(this.#params);
161+ // The trace is one sample per timestep, so a recording made at one dt
162+ // cannot be spliced onto one made at another: it would be two different
163+ // sample rates in the same buffer.
164+ this.recorder.clear();
165+ }
166+
167+ /** The largest timestep this scheme is stable at on this grid and medium. */
168+ get dtLimit(): number {
169+ return stableDt(this.grid.h, this.#scene.cmax, this.model.order, 1);
170+ }
171+
172+ #applyDt(): void {
173+ this.gpu.setDt(stableDt(this.grid.h, this.#scene.cmax, this.model.order, this.#cfl));
174+ }
175+
176+ /**
177+ * Swap the medium under a running model. It is data, not code, so this
178+ * needs no recompile — but the timestep follows from it, and changing dt
179+ * mid-run would leave the leapfrog's two histories half a step apart, so
180+ * the caller is expected to reset afterwards.
181+ */
182+ setScene(sceneModel: MScene, params: Params, source?: string): void {
183+ const built = Scene.create({
184+ grid: this.grid,
185+ source: source ?? sceneModel.source,
186+ paramNames: sceneModel.params.map((p) => p.key),
187+ params,
188+ });
189+ this.#scene = built;
190+ this.#sceneModel = sceneModel;
191+ this.gpu.uploadMedium(built);
192+ this.#applyDt();
193+ this.gpu.setParams(this.#params);
194+ this.recorder.clear();
195+ }
196+
197+ setParams(params: ModelParams): void {
198+ this.#params = params;
199+ this.gpu.setParams(params);
200+ }
201+
202+ /** Run `init`: a silent grid at t = 0. */
203+ reset(): void {
204+ this.gpu.init();
205+ this.recorder.clear();
206+ this.t = 0;
207+ this.steps = 0;
208+ }
209+
210+ /** Put the microphone at the grid point nearest the given coordinates. */
211+ setMic(x: number, y: number): void {
212+ const { L, h } = this.grid;
213+ this.recorder.setProbe((x + L / 2) / h - 0.5, (y + L / 2) / h - 0.5);
214+ }
215+
216+ /** Advance `n` steps. Synchronous: records and submits, nothing read back.
217+ * The microphone samples inside the same submission, once per step. */
218+ step(n = 1): void {
219+ this.gpu.step(n, (enc) => this.recorder.encode(enc));
220+ this.t += n * this.dt;
221+ this.steps += n;
222+ }
223+
224+ /** Wait for the submitted steps to finish, without reading anything back. */
225+ sync(): Promise<undefined> {
226+ return this.device.queue.onSubmittedWorkDone();
227+ }
228+
229+ /** Read a named field back to the CPU. */
230+ read(name: string): Promise<Float32Array> {
231+ return this.gpu.read(name);
232+ }
233+
234+ describe(): { init: string[]; step: string[] } {
235+ return this.gpu.describe();
236+ }
237+
238+ destroy(): void {
239+ this.recorder.destroy();
240+ this.gpu.destroy();
241+ }
242+}
src/mgpu/stencil.tsadded+150−0View file
@@ -0,0 +1,150 @@
1+/**
2+ * The Laplacian stencils, as GPU dispatches.
3+ *
4+ * A grid field is an npts x 1 column vector in the buffer, laid out x-fastest:
5+ * the point (ix, iy) is element `ix + nx*iy`. That is the only place in the
6+ * project where the layout matters to anything but the renderer, which is the
7+ * reason the stencils are host-provided operations rather than something a .m
8+ * expresses with array slicing.
9+ *
10+ * Outside the grid the field is taken to be zero, which makes the outer
11+ * boundary sound-soft (a pressure-release wall). Every scene puts an absorbing
12+ * layer in front of it (tools/sponge.m), so in a well set-up run almost nothing
13+ * reaches the wall to be reflected; what does is attenuated on the way out and
14+ * again on the way back. See the README on how good an open boundary that is.
15+ *
16+ * Both stencils divide by h^2, so a .m reads `lap2(p)` as the Laplacian in the
17+ * physical units the scene's coordinates are in. The grid spacing is compiled
18+ * into the shader — the grid is fixed when a model is compiled, and changing it
19+ * recompiles anyway.
20+ */
21+import { UnsupportedOnGpu, WORKGROUP_SIZE } from './wgsl.ts';
22+
23+export type StencilKind = 'lap2' | 'lap4';
24+
25+export interface StencilGrid {
26+ nx: number;
27+ ny: number;
28+ /** Grid spacing, the same in x and y. */
29+ h: number;
30+}
31+
32+/**
33+ * The 5-point second-order Laplacian, and the 9-point fourth-order one.
34+ *
35+ * The fourth-order stencil is the standard (-1/12, 4/3, -5/2, 4/3, -1/12)/h^2
36+ * one-dimensional second derivative applied along each axis, so it is a
37+ * five-point line in x plus a five-point line in y (nine points in all, not a
38+ * 3x3 block). It costs twice the reads of `lap2` and buys two orders of
39+ * accuracy, which for a wave problem shows up as much less grid dispersion:
40+ * a pulse travelling many wavelengths stays a pulse instead of trailing
41+ * numerical ripples. See models/leapfrog4.m.
42+ */
43+function stencilWGSL(kind: StencilKind, grid: StencilGrid, npts: number): string {
44+ const { nx, ny, h } = grid;
45+ const inv = 1 / (h * h);
46+ // A read of a point outside the grid returns zero.
47+ const at = `
48+fn at(ix: i32, iy: i32) -> f32 {
49+ if (ix < 0 || ix >= ${nx} || iy < 0 || iy >= ${ny}) { return 0.0; }
50+ return src[u32(ix + ${nx} * iy)];
51+}
52+`;
53+ const body =
54+ kind === 'lap2'
55+ ? ` let s = at(ix - 1, iy) + at(ix + 1, iy) + at(ix, iy - 1) + at(ix, iy + 1)
56+ - 4.0 * at(ix, iy);`
57+ : ` let s = (-1.0 / 12.0) * (at(ix - 2, iy) + at(ix + 2, iy) + at(ix, iy - 2) + at(ix, iy + 2))
58+ + (4.0 / 3.0) * (at(ix - 1, iy) + at(ix + 1, iy) + at(ix, iy - 1) + at(ix, iy + 1))
59+ - 5.0 * at(ix, iy);`;
60+ return `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
61+@group(0) @binding(1) var<storage, read> src: array<f32>;
62+${at}
63+@compute @workgroup_size(${WORKGROUP_SIZE})
64+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
65+ let i = gid.x;
66+ if (i >= ${npts}u) { return; }
67+ let ix = i32(i % ${nx}u);
68+ let iy = i32(i / ${nx}u);
69+${body}
70+ dst[i] = ${inv}f * s;
71+}
72+`;
73+}
74+
75+/** Compiled stencil pipelines for one grid. Shared by every plan on it. */
76+export class StencilPlan {
77+ readonly grid: StencilGrid;
78+ readonly npts: number;
79+
80+ #device: GPUDevice;
81+ #layout: GPUBindGroupLayout;
82+ #pipelines = new Map<StencilKind, GPUComputePipeline>();
83+
84+ private constructor(device: GPUDevice, grid: StencilGrid, layout: GPUBindGroupLayout) {
85+ this.#device = device;
86+ this.grid = grid;
87+ this.npts = grid.nx * grid.ny;
88+ this.#layout = layout;
89+ }
90+
91+ static create(device: GPUDevice, grid: StencilGrid): StencilPlan {
92+ const layout = device.createBindGroupLayout({
93+ label: 'stencil',
94+ entries: [
95+ {
96+ binding: 0,
97+ visibility: GPUShaderStage.COMPUTE,
98+ buffer: { type: 'storage' },
99+ },
100+ {
101+ binding: 1,
102+ visibility: GPUShaderStage.COMPUTE,
103+ buffer: { type: 'read-only-storage' },
104+ },
105+ ],
106+ });
107+ return new StencilPlan(device, grid, layout);
108+ }
109+
110+ /**
111+ * The pipeline for one stencil, compiled on first use. Both are cheap, but
112+ * a model uses one of them, and compiling only what is asked for keeps the
113+ * op sequence honest about what the .m actually costs.
114+ */
115+ async pipeline(kind: StencilKind): Promise<GPUComputePipeline> {
116+ const existing = this.#pipelines.get(kind);
117+ if (existing) return existing;
118+ const code = stencilWGSL(kind, this.grid, this.npts);
119+ const module = this.#device.createShaderModule({ code, label: kind });
120+ let pipeline: GPUComputePipeline;
121+ try {
122+ pipeline = await this.#device.createComputePipelineAsync({
123+ layout: this.#device.createPipelineLayout({ bindGroupLayouts: [this.#layout] }),
124+ compute: { module, entryPoint: 'main' },
125+ label: kind,
126+ });
127+ } catch (e) {
128+ // No error scope here either; see makePipeline in plan.ts.
129+ throw new UnsupportedOnGpu(
130+ `stencil '${kind}': ${e instanceof Error ? e.message : String(e)}`,
131+ );
132+ }
133+ this.#pipelines.set(kind, pipeline);
134+ return pipeline;
135+ }
136+
137+ createBinding(src: GPUBuffer, dst: GPUBuffer): GPUBindGroup {
138+ return this.#device.createBindGroup({
139+ layout: this.#layout,
140+ entries: [
141+ { binding: 0, resource: { buffer: dst } },
142+ { binding: 1, resource: { buffer: src } },
143+ ],
144+ });
145+ }
146+
147+ get workgroups(): number {
148+ return Math.ceil(this.npts / WORKGROUP_SIZE);
149+ }
150+}
src/mgpu/wgsl.tsadded+420−0View file
@@ -0,0 +1,420 @@
1+/**
2+ * IR expression tree -> one WGSL compute kernel.
3+ *
4+ * This is the WebGPU counterpart of numbl's C-side fused emitter
5+ * (`codegen/emitTensorFused.ts`): for an `Assign` whose right-hand side is
6+ * purely element-wise over operands of the target's shape, emit a single
7+ * kernel that computes one output element per invocation. Because numbl's
8+ * inline pass has already folded the ANF temps back together, one source line
9+ * of MATLAB becomes one kernel.
10+ *
11+ * Everything is f32, matching the existing fp32 WebGPU transform backend.
12+ */
13+import { getBuiltin } from 'numbl-src/numbl-core/jit/builtins/index.ts';
14+import { isMultiElement } from 'numbl-src/numbl-core/jit/lowering/types.ts';
15+import type { IRExpr, Assign } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
16+import type { NumericType, Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
17+
18+/** Raised for a construct the WGSL backend cannot express. Mirrors numbl's
19+ * own decline discipline: fail at compile time with a source span, never
20+ * silently produce something that computes the wrong thing. */
21+export class UnsupportedOnGpu extends Error {
22+ readonly span?: unknown;
23+ constructor(message: string, span?: unknown) {
24+ super(message);
25+ this.name = 'UnsupportedOnGpu';
26+ this.span = span;
27+ }
28+}
29+
30+const isNumeric = (t: Type): t is NumericType => t.kind === 'Numeric';
31+const isTensor = (t: Type): boolean => isNumeric(t) && isMultiElement(t);
32+
33+/** Element-wise binary builtins -> WGSL infix operator. */
34+const BINARY_OPS: Record<string, string> = {
35+ plus: '+',
36+ minus: '-',
37+ times: '*',
38+ rdivide: '/',
39+ // Degenerate to element-wise when at least one side is a scalar; the
40+ // both-tensor (true matrix) case is rejected below.
41+ mtimes: '*',
42+ mrdivide: '/',
43+};
44+
45+/** Element-wise unary builtins -> WGSL prefix operator. */
46+const UNARY_OPS: Record<string, string> = { uminus: '-', uplus: '+' };
47+
48+/** Element-wise builtin calls -> WGSL builtin of the same arity. */
49+const CALL_FNS: Record<string, string> = {
50+ abs: 'abs',
51+ acos: 'acos',
52+ asin: 'asin',
53+ atan: 'atan',
54+ atan2: 'atan2',
55+ ceil: 'ceil',
56+ cos: 'cos',
57+ cosh: 'cosh',
58+ exp: 'exp',
59+ floor: 'floor',
60+ log: 'log',
61+ log2: 'log2',
62+ max: 'max',
63+ min: 'min',
64+ round: 'round',
65+ sign: 'sign',
66+ sin: 'sin',
67+ sinh: 'sinh',
68+ sqrt: 'sqrt',
69+ tan: 'tan',
70+ tanh: 'tanh',
71+};
72+
73+/** Zero-argument builtins that are compile-time constants. numbl lowers `pi`
74+ * as a call rather than folding it, so the backend is where it becomes a
75+ * number. */
76+const CONST_FNS: Record<string, number> = { pi: Math.PI };
77+
78+/** WGSL f32 literal. Must always carry a decimal point or exponent, or WGSL
79+ * infers AbstractInt and rejects the mixed-type arithmetic. */
80+function f32Lit(v: number): string {
81+ if (!Number.isFinite(v)) {
82+ throw new UnsupportedOnGpu(`cannot emit non-finite literal ${v}`);
83+ }
84+ return Number.isInteger(v) && Math.abs(v) < 1e21
85+ ? `${v}.0`
86+ : String(v).includes('e')
87+ ? `${v}f`
88+ : String(v);
89+}
90+
91+/** How a scalar or tensor operand is read inside the kernel. */
92+export interface KernelInputs {
93+ /** cName -> storage binding index, for multi-element tensor operands. */
94+ tensors: Map<string, number>;
95+ /** cName -> slot in the params storage buffer, for runtime scalars. */
96+ params: Map<string, number>;
97+ /** cName -> defining expression, for scalars the .m computes from
98+ * parameters (`us = a + b`). These have no buffer and no param slot; they
99+ * become `let` bindings in the prologue of every kernel that reads them. */
100+ scalars: Map<string, { name: string; expr: IRExpr }>;
101+}
102+
103+/** Mutable state while emitting one kernel. */
104+interface Ctx {
105+ io: KernelInputs;
106+ /** `let` lines to emit before the body, in dependency order. */
107+ prologue: string[];
108+ /** cName -> WGSL identifier, for scalars already bound in the prologue. */
109+ bound: Map<string, string>;
110+}
111+
112+/** WGSL identifier for a derived scalar. Avoids a leading underscore, which
113+ * WGSL reserves. */
114+const scalarIdent = (cName: string): string =>
115+ `s_${cName.replace(/[^A-Za-z0-9_]/g, '_')}`;
116+
117+/**
118+ * Bind a .m-derived scalar in the prologue (once), after whatever it depends
119+ * on, and return its identifier.
120+ */
121+function bindScalar(cName: string, ctx: Ctx): string {
122+ const already = ctx.bound.get(cName);
123+ if (already) return already;
124+ const def = ctx.io.scalars.get(cName)!;
125+ const ident = scalarIdent(cName);
126+ // Claim the name before emitting the RHS so a (malformed) self-reference
127+ // cannot recurse forever.
128+ ctx.bound.set(cName, ident);
129+ const rhs = emitExpr(def.expr, ctx);
130+ ctx.prologue.push(` let ${ident} = ${rhs};`);
131+ return ident;
132+}
133+
134+/**
135+ * Emit the per-element WGSL expression for `e`. `i` is the element index
136+ * variable in scope.
137+ */
138+function emitExpr(e: IRExpr, ctx: Ctx): string {
139+ const io = ctx.io;
140+ switch (e.kind) {
141+ case 'NumLit':
142+ return f32Lit(e.value);
143+
144+ case 'Var': {
145+ if (isTensor(e.ty)) {
146+ const slot = io.tensors.get(e.cName);
147+ if (slot === undefined) {
148+ throw new UnsupportedOnGpu(`no buffer bound for '${e.name}'`, e.span);
149+ }
150+ return `in${slot}[i]`;
151+ }
152+ // Scalar: either an exact compile-time value or a runtime parameter.
153+ if (isNumeric(e.ty) && typeof e.ty.exact === 'number') {
154+ return f32Lit(e.ty.exact);
155+ }
156+ const slot = io.params.get(e.cName);
157+ if (slot !== undefined) return `prm[${slot}]`;
158+ if (io.scalars.has(e.cName)) return bindScalar(e.cName, ctx);
159+ throw new UnsupportedOnGpu(
160+ `scalar '${e.name}' is not a constant, a parameter, or computed in ` +
161+ `this model`,
162+ e.span,
163+ );
164+ }
165+
166+ case 'Binary': {
167+ if ((e.builtin === 'mtimes' || e.builtin === 'mrdivide') &&
168+ isTensor(e.left.ty) && isTensor(e.right.ty)) {
169+ throw new UnsupportedOnGpu(
170+ `matrix '${e.builtin === 'mtimes' ? '*' : '/'}' is not supported; ` +
171+ `use the element-wise form ('.${e.builtin === 'mtimes' ? '*' : '/'}')`,
172+ e.span,
173+ );
174+ }
175+ if (e.builtin === 'power' || e.builtin === 'mpower') {
176+ return emitPower(e.left, e.right, ctx, e.span);
177+ }
178+ const op = BINARY_OPS[e.builtin];
179+ if (!op) {
180+ throw new UnsupportedOnGpu(`operator '${e.builtin}' is not supported`, e.span);
181+ }
182+ return `(${emitExpr(e.left, ctx)} ${op} ${emitExpr(e.right, ctx)})`;
183+ }
184+
185+ case 'Unary': {
186+ const op = UNARY_OPS[e.builtin];
187+ if (!op) {
188+ throw new UnsupportedOnGpu(`unary '${e.builtin}' is not supported`, e.span);
189+ }
190+ return `(${op}${emitExpr(e.operand, ctx)})`;
191+ }
192+
193+ case 'Call': {
194+ // A shape constructor used inside an element-wise expression
195+ // contributes the same constant at every slot, so it needs no buffer.
196+ // (The shape itself is validated against the target by checkShapes.)
197+ if (e.name === 'ones') return '1.0';
198+ if (e.name === 'zeros') return '0.0';
199+
200+ const konst = CONST_FNS[e.name];
201+ if (konst !== undefined && e.args.length === 0) return f32Lit(konst);
202+
203+ const fn = CALL_FNS[e.name];
204+ const b = getBuiltin(e.name);
205+ if (!fn || !b?.elementwise) {
206+ // A call numbl resolved to another function in the file gets a mangled
207+ // specialization name; a builtin keeps its source-level name. Only the
208+ // model's entry points are compiled, so a helper is a distinct failure
209+ // from an unsupported builtin and deserves to say so.
210+ const isUserFunction = e.cName !== e.name;
211+ throw new UnsupportedOnGpu(
212+ isUserFunction
213+ ? `'${e.name}' is a function defined in this model. Only init and ` +
214+ `step are compiled — inline its body into the caller.`
215+ : `'${e.name}' cannot be evaluated element-wise on the GPU`,
216+ e.span,
217+ );
218+ }
219+ return `${fn}(${e.args.map((a) => emitExpr(a, ctx)).join(', ')})`;
220+ }
221+
222+ default:
223+ throw new UnsupportedOnGpu(`'${e.kind}' is not supported on the GPU`, e.span);
224+ }
225+}
226+
227+/**
228+ * `x.^k`. WGSL's `pow` is undefined for a negative base, and these fields go
229+ * negative routinely, so expand small non-negative integer exponents into
230+ * repeated multiplication — which is also what makes `u.^2` free.
231+ */
232+function emitPower(base: IRExpr, exponent: IRExpr, ctx: Ctx, span: unknown): string {
233+ const k =
234+ exponent.kind === 'NumLit'
235+ ? exponent.value
236+ : isNumeric(exponent.ty) && typeof exponent.ty.exact === 'number'
237+ ? exponent.ty.exact
238+ : undefined;
239+ const b = emitExpr(base, ctx);
240+ if (k !== undefined && Number.isInteger(k) && k >= 0 && k <= 8) {
241+ if (k === 0) return '1.0';
242+ // bind once so a compound base expression is not re-evaluated k times
243+ return `pow_i${k}(${b})`;
244+ }
245+ if (k !== undefined && Number.isInteger(k) && k < 0 && k >= -8) {
246+ return `(1.0 / pow_i${-k}(${b}))`;
247+ }
248+ throw new UnsupportedOnGpu(
249+ `'.^' needs a literal integer exponent in [-8, 8] (got ` +
250+ `${k === undefined ? 'a runtime value' : k}); a negative base makes ` +
251+ `WGSL's pow() undefined`,
252+ span,
253+ );
254+}
255+
256+/** Fixed-exponent power helpers, emitted only when used. */
257+function powHelpers(used: Set<number>): string {
258+ const out: string[] = [];
259+ for (const k of [...used].sort((a, b) => a - b)) {
260+ const body =
261+ k === 1 ? 'x' : `x${' * x'.repeat(k - 1)}`;
262+ out.push(`fn pow_i${k}(x: f32) -> f32 { return ${body}; }`);
263+ }
264+ return out.join('\n');
265+}
266+
267+/**
268+ * Can this expression be evaluated inside a fused kernel?
269+ *
270+ * The predicate behind src/mgpu/fuse.ts, and the single statement of what the
271+ * emitter above accepts: everything here is something `emitExpr` can write out
272+ * per element, and everything it declines is something that needs its own
273+ * dispatch. Keep the two in step.
274+ */
275+export function isGpuFusableExpr(e: IRExpr): boolean {
276+ switch (e.kind) {
277+ case 'NumLit':
278+ return true;
279+ case 'Var':
280+ return isNumeric(e.ty);
281+ case 'Binary': {
282+ if (
283+ (e.builtin === 'mtimes' || e.builtin === 'mrdivide') &&
284+ isTensor(e.left.ty) && isTensor(e.right.ty)
285+ ) {
286+ return false;
287+ }
288+ if (e.builtin === 'power' || e.builtin === 'mpower') {
289+ return isGpuFusableExpr(e.left) && isGpuFusableExpr(e.right);
290+ }
291+ return (
292+ e.builtin in BINARY_OPS && isGpuFusableExpr(e.left) && isGpuFusableExpr(e.right)
293+ );
294+ }
295+ case 'Unary':
296+ return e.builtin in UNARY_OPS && isGpuFusableExpr(e.operand);
297+ case 'Call': {
298+ if (e.name === 'zeros' || e.name === 'ones') return true;
299+ if (e.name in CONST_FNS) return e.args.length === 0;
300+ return e.name in CALL_FNS && e.args.every(isGpuFusableExpr);
301+ }
302+ default:
303+ return false;
304+ }
305+}
306+
307+/**
308+ * Reject implicit expansion (broadcasting).
309+ *
310+ * numbl's lowering permits it — `2x4096 .* 1x4096` lowers happily with MATLAB
311+ * expansion semantics — but a kernel that walks one linear index across every
312+ * operand would quietly compute the wrong thing. So every multi-element
313+ * operand must have exactly the target's shape. Scalars are fine: they are
314+ * read from the params buffer or folded in as literals.
315+ */
316+export function checkShapes(e: IRExpr, target: NumericType, name: string): void {
317+ const want = target.shape;
318+ const same = (t: NumericType): boolean => {
319+ const got = t.shape;
320+ return (
321+ !!want && !!got && want.length === got.length &&
322+ want.every((d, i) => d === got[i])
323+ );
324+ };
325+ const walk = (x: IRExpr): void => {
326+ if (isNumeric(x.ty) && isMultiElement(x.ty) && !same(x.ty)) {
327+ const got = x.ty.shape?.join('x') ?? 'dynamic';
328+ throw new UnsupportedOnGpu(
329+ `'${name}' would need implicit expansion: an operand is ${got} but the ` +
330+ `result is ${want?.join('x') ?? 'dynamic'}. Expand it explicitly ` +
331+ `(the GPU kernel walks one index across every operand).`,
332+ x.span,
333+ );
334+ }
335+ switch (x.kind) {
336+ case 'Binary':
337+ walk(x.left);
338+ walk(x.right);
339+ return;
340+ case 'Unary':
341+ walk(x.operand);
342+ return;
343+ case 'Call':
344+ // A shape constructor's own arguments are sizes, not data.
345+ if (x.name !== 'ones' && x.name !== 'zeros') x.args.forEach(walk);
346+ return;
347+ default:
348+ return;
349+ }
350+ };
351+ walk(e);
352+}
353+
354+export const WORKGROUP_SIZE = 64;
355+
356+export interface Kernel {
357+ code: string;
358+ /** Number of output elements. */
359+ count: number;
360+ label: string;
361+}
362+
363+/**
364+ * Build the kernel for one element-wise `Assign`. `io` must already map every
365+ * tensor operand cName to a binding index and every runtime scalar to a
366+ * params slot; the output is binding 0 and the params buffer is the binding
367+ * after the last input.
368+ */
369+export function buildKernel(
370+ stmt: Assign,
371+ io: KernelInputs,
372+ count: number,
373+ label: string,
374+): Kernel {
375+ if (!isNumeric(stmt.ty)) {
376+ throw new UnsupportedOnGpu(`'${stmt.name}' is not a numeric array`, stmt.span);
377+ }
378+ if (stmt.ty.isComplex) {
379+ throw new UnsupportedOnGpu(
380+ `'${stmt.name}' is complex; the GPU backend is real-only (a spectral ` +
381+ `field is carried as a real 2 x nlm array)`,
382+ stmt.span,
383+ );
384+ }
385+
386+ checkShapes(stmt.expr, stmt.ty, stmt.name);
387+ const ctx: Ctx = { io, prologue: [], bound: new Map() };
388+ const body = emitExpr(stmt.expr, ctx);
389+
390+ // pow_iK helpers are discovered during emission; scan the result for them.
391+ const used = new Set<number>();
392+ const emitted = [...ctx.prologue, body].join('\n');
393+ for (const m of emitted.matchAll(/\bpow_i(\d+)\(/g)) used.add(Number(m[1]));
394+
395+ const decls = [`@group(0) @binding(0) var<storage, read_write> out: array<f32>;`];
396+ for (const [, slot] of io.tensors) {
397+ decls.push(
398+ `@group(0) @binding(${slot + 1}) var<storage, read> in${slot}: array<f32>;`,
399+ );
400+ }
401+ // Params live in a read-only storage buffer rather than a uniform block:
402+ // uniform arrays would need 16-byte element stride.
403+ const prmBinding = io.tensors.size + 1;
404+ decls.push(
405+ `@group(0) @binding(${prmBinding}) var<storage, read> prm: array<f32>;`,
406+ );
407+
408+ const code = `${decls.join('\n')}
409+
410+${powHelpers(used)}
411+
412+@compute @workgroup_size(${WORKGROUP_SIZE})
413+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
414+ let i = gid.x;
415+ if (i >= ${count}u) { return; }
416+${ctx.prologue.length ? `${ctx.prologue.join('\n')}\n` : ''} out[i] = ${body};
417+}
418+`;
419+ return { code, count, label };
420+}
src/raw.d.tsadded+15−0View file
@@ -0,0 +1,15 @@
1+/** Vite's `?raw` suffix imports a file's text. Used to load .m model sources. */
2+declare module '*?raw' {
3+ const source: string;
4+ export default source;
5+}
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/render/colorbar.tsadded+47−0View file
@@ -0,0 +1,47 @@
1+import type { ColormapFunc } from './colormaps.ts';
2+
3+/** Compact numeric label: 3 significant digits, trailing zeros trimmed. */
4+export const fmtValue = (v: number): string =>
5+ Number.isFinite(v) ? v.toPrecision(3).replace(/\.?0+$/, '') : '—';
6+
7+/** Vertical colorbar drawn on a small canvas, with min/max labels.
8+ * Adapted from turing-surface's src/render/colorbar.ts. */
9+export class Colorbar {
10+ #canvas: HTMLCanvasElement;
11+ #minLabel: HTMLElement;
12+ #maxLabel: HTMLElement;
13+
14+ constructor(container: HTMLElement) {
15+ container.classList.add('colorbar');
16+ this.#maxLabel = document.createElement('div');
17+ this.#maxLabel.className = 'colorbar-label';
18+ this.#canvas = document.createElement('canvas');
19+ this.#canvas.width = 12;
20+ this.#canvas.height = 160;
21+ this.#minLabel = document.createElement('div');
22+ this.#minLabel.className = 'colorbar-label';
23+ container.append(this.#maxLabel, this.#canvas, this.#minLabel);
24+ }
25+
26+ /** Repaint the gradient. Only when the colormap changes: it is 160 filled
27+ * rows, and the frame loop has better things to do. */
28+ setColormap(cmap: ColormapFunc): void {
29+ const ctx = this.#canvas.getContext('2d');
30+ if (!ctx) return;
31+ const h = this.#canvas.height;
32+ for (let y = 0; y < h; y++) {
33+ const t = 1 - y / (h - 1);
34+ const [r, g, b] = cmap(t);
35+ ctx.fillStyle = `rgb(${r},${g},${b})`;
36+ ctx.fillRect(0, y, this.#canvas.width, 1);
37+ }
38+ }
39+
40+ /** The end labels, which do change as the scale follows the field. */
41+ setRange(vmin: number, vmax: number): void {
42+ const lo = fmtValue(vmin);
43+ const hi = fmtValue(vmax);
44+ if (this.#minLabel.textContent !== lo) this.#minLabel.textContent = lo;
45+ if (this.#maxLabel.textContent !== hi) this.#maxLabel.textContent = hi;
46+ }
47+}
src/render/colormaps.tsadded+113−0View file
@@ -0,0 +1,113 @@
1+/**
2+ * Colormaps: each maps a normalized value in [0, 1] to [r, g, b] in [0, 255].
3+ * Adapted from figpack's SphereEmbedding view (figpack_experimental).
4+ */
5+
6+export type ColormapFunc = (t: number) => [number, number, number];
7+
8+const clamp01 = (t: number) => Math.max(0, Math.min(1, t));
9+
10+// Piecewise-linear interpolation through control points (r, g, b in 0-255)
11+const makeInterpolated = (stops: [number, number, number][]): ColormapFunc => {
12+ const n = stops.length;
13+ return (t: number) => {
14+ t = clamp01(t);
15+ const x = t * (n - 1);
16+ const i = Math.min(n - 2, Math.floor(x));
17+ const f = x - i;
18+ const a = stops[i];
19+ const b = stops[i + 1];
20+ return [
21+ Math.round(a[0] + (b[0] - a[0]) * f),
22+ Math.round(a[1] + (b[1] - a[1]) * f),
23+ Math.round(a[2] + (b[2] - a[2]) * f),
24+ ];
25+ };
26+};
27+
28+// Control points sampled from matplotlib colormaps
29+const viridis = makeInterpolated([
30+ [68, 1, 84],
31+ [72, 40, 120],
32+ [62, 74, 137],
33+ [49, 104, 142],
34+ [38, 130, 142],
35+ [31, 158, 137],
36+ [53, 183, 121],
37+ [109, 205, 89],
38+ [180, 222, 44],
39+ [253, 231, 37],
40+]);
41+
42+const plasma = makeInterpolated([
43+ [13, 8, 135],
44+ [84, 2, 163],
45+ [139, 10, 165],
46+ [185, 50, 137],
47+ [219, 92, 104],
48+ [244, 136, 73],
49+ [254, 188, 43],
50+ [240, 249, 33],
51+]);
52+
53+const inferno = makeInterpolated([
54+ [0, 0, 4],
55+ [40, 11, 84],
56+ [101, 21, 110],
57+ [159, 42, 99],
58+ [212, 72, 66],
59+ [245, 125, 21],
60+ [250, 193, 39],
61+ [252, 255, 164],
62+]);
63+
64+const coolwarm = makeInterpolated([
65+ [59, 76, 192],
66+ [124, 159, 249],
67+ [192, 212, 245],
68+ [242, 242, 242],
69+ [245, 195, 157],
70+ [222, 96, 77],
71+ [180, 4, 38],
72+]);
73+
74+// matplotlib's `seismic`: harder contrast about the middle than coolwarm, and
75+// dark at both ends, which suits a wavefield whose interesting parts are the
76+// extremes.
77+const seismic = makeInterpolated([
78+ [0, 0, 76],
79+ [0, 0, 255],
80+ [255, 255, 255],
81+ [255, 0, 0],
82+ [128, 0, 0],
83+]);
84+
85+const jet = makeInterpolated([
86+ [0, 0, 128],
87+ [0, 0, 255],
88+ [0, 255, 255],
89+ [0, 255, 0],
90+ [255, 255, 0],
91+ [255, 0, 0],
92+ [128, 0, 0],
93+]);
94+
95+const grayscale: ColormapFunc = (t: number) => {
96+ const v = Math.round(clamp01(t) * 255);
97+ return [v, v, v];
98+};
99+
100+/** Diverging maps first: the pressure field is signed and is drawn
101+ * symmetrically about zero, so a map with a distinct middle is what makes
102+ * the wavefronts read. */
103+export const colormaps: Record<string, ColormapFunc> = {
104+ coolwarm,
105+ seismic,
106+ grayscale,
107+ viridis,
108+ plasma,
109+ inferno,
110+ jet,
111+};
112+
113+export const colormapNames = Object.keys(colormaps);
src/render/field.tsadded+288−0View file
@@ -0,0 +1,288 @@
1+/**
2+ * Drawing the pressure field, straight out of the buffer the solver wrote.
3+ *
4+ * There is no readback in the display path: the fragment shader reads the
5+ * solver's storage buffer directly, so a frame costs one draw call and no
6+ * GPU->CPU round trip. (The host does read the field back occasionally, but
7+ * only to decide the colour scale — see the autoscale in main.ts.)
8+ *
9+ * Two things are drawn at once. The pressure goes through a diverging
10+ * colormap about zero, symmetric because the field is signed and a wave is
11+ * not more interesting on one side of zero than the other. The medium is
12+ * blended in underneath as a grey wash proportional to how far the local
13+ * sound speed departs from the background, which shows a hard scatterer as a
14+ * distinct shape and a smooth one as a soft cloud without either needing its
15+ * own kind of drawing.
16+ *
17+ * Sampling is bilinear rather than nearest, so a 256-point grid on a 700-pixel
18+ * canvas looks like a wave rather than like a grid.
19+ */
20+import type { ColormapFunc } from './colormaps.ts';
21+
22+const SHADER = `
23+struct View {
24+ nx: u32,
25+ ny: u32,
26+ scale: f32,
27+ cref: f32,
28+ cdev: f32,
29+ medium: f32,
30+ micX: f32,
31+ micY: f32,
32+ micOn: f32,
33+};
34+
35+@group(0) @binding(0) var<uniform> view: View;
36+@group(0) @binding(1) var<storage, read> field: array<f32>;
37+@group(0) @binding(2) var<storage, read> speed: array<f32>;
38+@group(0) @binding(3) var cmap: texture_2d<f32>;
39+@group(0) @binding(4) var cmapSampler: sampler;
40+
41+struct VSOut {
42+ @builtin(position) pos: vec4f,
43+ @location(0) uv: vec2f,
44+};
45+
46+@vertex
47+fn vs(@builtin(vertex_index) vi: u32) -> VSOut {
48+ // One oversized triangle covering the viewport.
49+ var xy = array<vec2f, 3>(vec2f(-1.0, -3.0), vec2f(-1.0, 1.0), vec2f(3.0, 1.0));
50+ var out: VSOut;
51+ let p = xy[vi];
52+ out.pos = vec4f(p, 0.0, 1.0);
53+ // uv (0,0) at the top-left corner of the domain.
54+ out.uv = vec2f((p.x + 1.0) * 0.5, (1.0 - p.y) * 0.5);
55+ return out;
56+}
57+
58+fn idx(ix: i32, iy: i32) -> u32 {
59+ let cx = clamp(ix, 0, i32(view.nx) - 1);
60+ let cy = clamp(iy, 0, i32(view.ny) - 1);
61+ return u32(cx + i32(view.nx) * cy);
62+}
63+
64+fn bilinear(gx: f32, gy: f32, isField: bool) -> f32 {
65+ let x0 = i32(floor(gx));
66+ let y0 = i32(floor(gy));
67+ let fx = gx - f32(x0);
68+ let fy = gy - f32(y0);
69+ var v00: f32; var v10: f32; var v01: f32; var v11: f32;
70+ if (isField) {
71+ v00 = field[idx(x0, y0)];
72+ v10 = field[idx(x0 + 1, y0)];
73+ v01 = field[idx(x0, y0 + 1)];
74+ v11 = field[idx(x0 + 1, y0 + 1)];
75+ } else {
76+ v00 = speed[idx(x0, y0)];
77+ v10 = speed[idx(x0 + 1, y0)];
78+ v01 = speed[idx(x0, y0 + 1)];
79+ v11 = speed[idx(x0 + 1, y0 + 1)];
80+ }
81+ return mix(mix(v00, v10, fx), mix(v01, v11, fx), fy);
82+}
83+
84+@fragment
85+fn fs(in: VSOut) -> @location(0) vec4f {
86+ // Row 0 of the buffer is the bottom of the picture: y increases upward in
87+ // the grid, downward on the screen.
88+ let gx = in.uv.x * f32(view.nx) - 0.5;
89+ let gy = (1.0 - in.uv.y) * f32(view.ny) - 0.5;
90+
91+ let p = bilinear(gx, gy, true);
92+ let t = clamp(0.5 + 0.5 * p / max(view.scale, 1e-20), 0.0, 1.0);
93+ var rgb = textureSample(cmap, cmapSampler, vec2f(t, 0.5)).rgb;
94+
95+ if (view.medium > 0.0) {
96+ let c = bilinear(gx, gy, false);
97+ let m = clamp(abs(c - view.cref) / max(view.cdev, 1e-20), 0.0, 1.0);
98+ rgb = mix(rgb, vec3f(0.35, 0.35, 0.38), view.medium * m);
99+ }
100+
101+ // The microphone, as a ring sized in grid cells so it stays the same size
102+ // on screen whatever the canvas is scaled to.
103+ if (view.micOn > 0.0) {
104+ let d = length(vec2f(gx - view.micX, gy - view.micY));
105+ let r = 0.022 * f32(view.nx);
106+ let w = 0.005 * f32(view.nx);
107+ if (abs(d - r) < w) {
108+ rgb = mix(rgb, vec3f(1.0, 1.0, 1.0), 0.9);
109+ } else if (d < 0.006 * f32(view.nx)) {
110+ rgb = mix(rgb, vec3f(1.0, 1.0, 1.0), 0.9);
111+ }
112+ }
113+ return vec4f(rgb, 1.0);
114+}
115+`;
116+
117+export interface FieldViewOptions {
118+ device: GPUDevice;
119+ canvas: HTMLCanvasElement;
120+ nx: number;
121+ ny: number;
122+}
123+
124+export class FieldView {
125+ readonly canvas: HTMLCanvasElement;
126+
127+ #device: GPUDevice;
128+ #context: GPUCanvasContext;
129+ #pipeline: GPURenderPipeline;
130+ #layout: GPUBindGroupLayout;
131+ #uniform: GPUBuffer;
132+ #uniformData = new ArrayBuffer(48);
133+ #sampler: GPUSampler;
134+ #cmapTexture: GPUTexture;
135+ #bindGroup: GPUBindGroup | null = null;
136+ #nx: number;
137+ #ny: number;
138+
139+ constructor(opts: FieldViewOptions) {
140+ const { device, canvas } = opts;
141+ this.#device = device;
142+ this.canvas = canvas;
143+ this.#nx = opts.nx;
144+ this.#ny = opts.ny;
145+
146+ const context = canvas.getContext('webgpu');
147+ if (!context) throw new Error('this canvas has no WebGPU context');
148+ this.#context = context;
149+ const format = navigator.gpu.getPreferredCanvasFormat();
150+ context.configure({ device, format, alphaMode: 'opaque' });
151+
152+ this.#layout = device.createBindGroupLayout({
153+ label: 'field-view',
154+ entries: [
155+ { binding: 0, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'uniform' } },
156+ {
157+ binding: 1,
158+ visibility: GPUShaderStage.FRAGMENT,
159+ buffer: { type: 'read-only-storage' },
160+ },
161+ {
162+ binding: 2,
163+ visibility: GPUShaderStage.FRAGMENT,
164+ buffer: { type: 'read-only-storage' },
165+ },
166+ { binding: 3, visibility: GPUShaderStage.FRAGMENT, texture: { sampleType: 'float' } },
167+ { binding: 4, visibility: GPUShaderStage.FRAGMENT, sampler: { type: 'filtering' } },
168+ ],
169+ });
170+
171+ const module = device.createShaderModule({ code: SHADER, label: 'field-view' });
172+ this.#pipeline = device.createRenderPipeline({
173+ label: 'field-view',
174+ layout: device.createPipelineLayout({ bindGroupLayouts: [this.#layout] }),
175+ vertex: { module, entryPoint: 'vs' },
176+ fragment: { module, entryPoint: 'fs', targets: [{ format }] },
177+ primitive: { topology: 'triangle-list' },
178+ });
179+
180+ this.#uniform = device.createBuffer({
181+ label: 'field-view-uniform',
182+ size: this.#uniformData.byteLength,
183+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
184+ });
185+ this.#sampler = device.createSampler({ magFilter: 'linear', minFilter: 'linear' });
186+ this.#cmapTexture = device.createTexture({
187+ label: 'colormap',
188+ size: [256, 1],
189+ format: 'rgba8unorm',
190+ usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST,
191+ });
192+ }
193+
194+ /** Point the view at the buffers of a (new) simulation. */
195+ setSource(pressure: GPUBuffer, speed: GPUBuffer, nx: number, ny: number): void {
196+ this.#nx = nx;
197+ this.#ny = ny;
198+ this.#bindGroup = this.#device.createBindGroup({
199+ layout: this.#layout,
200+ entries: [
201+ { binding: 0, resource: { buffer: this.#uniform } },
202+ { binding: 1, resource: { buffer: pressure } },
203+ { binding: 2, resource: { buffer: speed } },
204+ { binding: 3, resource: this.#cmapTexture.createView() },
205+ { binding: 4, resource: this.#sampler },
206+ ],
207+ });
208+ }
209+
210+ setColormap(cmap: ColormapFunc): void {
211+ const data = new Uint8Array(256 * 4);
212+ for (let i = 0; i < 256; i++) {
213+ const [r, g, b] = cmap(i / 255);
214+ data[4 * i] = r;
215+ data[4 * i + 1] = g;
216+ data[4 * i + 2] = b;
217+ data[4 * i + 3] = 255;
218+ }
219+ this.#device.queue.writeTexture(
220+ { texture: this.#cmapTexture },
221+ data,
222+ { bytesPerRow: 256 * 4 },
223+ { width: 256, height: 1 },
224+ );
225+ }
226+
227+ /**
228+ * Draw one frame. `scale` is the pressure the colormap saturates at, in
229+ * both directions; `cref`/`cdev` describe the medium wash, and `medium` is
230+ * how strongly to apply it (0 turns it off).
231+ */
232+ draw(opts: {
233+ scale: number;
234+ cref: number;
235+ cdev: number;
236+ medium: number;
237+ /** Microphone position, in grid-index coordinates. */
238+ mic?: { ix: number; iy: number } | null;
239+ }): void {
240+ if (!this.#bindGroup) return;
241+ const u32 = new Uint32Array(this.#uniformData);
242+ const f32 = new Float32Array(this.#uniformData);
243+ u32[0] = this.#nx;
244+ u32[1] = this.#ny;
245+ f32[2] = opts.scale;
246+ f32[3] = opts.cref;
247+ f32[4] = opts.cdev;
248+ f32[5] = opts.medium;
249+ f32[6] = opts.mic?.ix ?? 0;
250+ f32[7] = opts.mic?.iy ?? 0;
251+ f32[8] = opts.mic ? 1 : 0;
252+ this.#device.queue.writeBuffer(this.#uniform, 0, this.#uniformData);
253+
254+ const encoder = this.#device.createCommandEncoder({ label: 'field-view' });
255+ const pass = encoder.beginRenderPass({
256+ colorAttachments: [
257+ {
258+ view: this.#context.getCurrentTexture().createView(),
259+ clearValue: { r: 0, g: 0, b: 0, a: 1 },
260+ loadOp: 'clear',
261+ storeOp: 'store',
262+ },
263+ ],
264+ });
265+ pass.setPipeline(this.#pipeline);
266+ pass.setBindGroup(0, this.#bindGroup);
267+ pass.draw(3);
268+ pass.end();
269+ this.#device.queue.submit([encoder.finish()]);
270+ }
271+
272+ /** Match the canvas's backing store to its CSS size. */
273+ resize(): void {
274+ const dpr = Math.min(window.devicePixelRatio || 1, 2);
275+ const rect = this.canvas.getBoundingClientRect();
276+ const w = Math.max(1, Math.round(rect.width * dpr));
277+ const h = Math.max(1, Math.round(rect.height * dpr));
278+ if (this.canvas.width !== w || this.canvas.height !== h) {
279+ this.canvas.width = w;
280+ this.canvas.height = h;
281+ }
282+ }
283+
284+ destroy(): void {
285+ this.#uniform.destroy();
286+ this.#cmapTexture.destroy();
287+ }
288+}
src/scene/registry.tsadded+195−0View file
@@ -0,0 +1,195 @@
1+/**
2+ * The available scenes: their MATLAB source, and the parameters the host
3+ * offers each one.
4+ *
5+ * Same split as the models (src/mgpu/registry.ts): what the medium *is* lives
6+ * in the .m, and the sliders around it live here. A scene also gets `x`, `y`
7+ * (metres), `L`, `h` (metres), `c0` (metres per second — the speed of sound
8+ * in air), `npts`, `nx` and `ny` for free — see src/scene/scene.ts.
9+ *
10+ * Every scene here takes its speed contrasts (`cin`, `cwall`, `dn`, `amp`) as
11+ * plain ratios to `c0` rather than as absolute speeds, so a slider means the
12+ * same thing whatever the background is. Its lengths (`R`, `side`, `gap`,
13+ * `w`, ...) are metres, and its rates (`absorb`) are inverse seconds.
14+ */
15+import type { ParamSpec, Params } from '../mgpu/registry.ts';
16+import diskSource from '../../scenes/disk.m?raw';
17+import roomSource from '../../scenes/room.m?raw';
18+import slitSource from '../../scenes/slit.m?raw';
19+import lensSource from '../../scenes/lens.m?raw';
20+import speckleSource from '../../scenes/speckle.m?raw';
21+
22+export interface MScene {
23+ key: string;
24+ label: string;
25+ blurb: string;
26+ params: ParamSpec[];
27+ /**
28+ * Model parameters this scene wants set when it is chosen — a source
29+ * position, usually. A scene cannot reach into the model's parameters, and
30+ * should not: they belong to different files. But a room is no use with the
31+ * source outside it, so a scene may say what it would like, and the app
32+ * applies it as if the sliders had been moved by hand.
33+ */
34+ suggest?: Params;
35+ /** Where this scene would like the microphone. Inside the room, for a room.
36+ * Applied like `suggest`, and movable afterwards. */
37+ mic?: { x: number; y: number };
38+ source: string;
39+}
40+
41+const disk: MScene = {
42+ key: 'disk',
43+ label: 'Disk',
44+ blurb: 'One circular scatterer — the case with an exact series solution.',
45+ params: [
46+ {
47+ key: 'cin',
48+ label: 'speed inside (×c0)',
49+ value: 0.5,
50+ min: 0.2,
51+ max: 3,
52+ step: 0.05,
53+ hint: 'Ratio to the background speed. Well above 1 is nearly rigid, well below nearly pressure-release, and exactly 1 is no scatterer at all.',
54+ },
55+ { key: 'R', label: 'radius (m)', value: 0.75, min: 0.1, max: 2.5, step: 0.05 },
56+ {
57+ key: 'absorb',
58+ label: 'absorption (1/s)',
59+ value: 0,
60+ min: 0,
61+ max: 1500,
62+ step: 50,
63+ hint: 'Damping inside the disk. Turn it up and the scatterer swallows what enters it instead of ringing.',
64+ },
65+ ],
66+ source: diskSource,
67+};
68+
69+const room: MScene = {
70+ key: 'room',
71+ label: 'Room',
72+ blurb: 'An enclosure with a doorway: echoes, modes, and something to listen to.',
73+ params: [
74+ {
75+ key: 'cwall',
76+ label: 'wall speed (×c0)',
77+ value: 0.15,
78+ min: 0.05,
79+ max: 6,
80+ step: 0.05,
81+ hint: 'Ratio to the background speed. Below 1 the wall is slower than the room and reflects without costing timestep; above 1 it is rigid, and the timestep drops to match. Either way what reflects is |cwall-1|/(cwall+1).',
82+ },
83+ {
84+ key: 'side',
85+ label: 'room half-width (m)',
86+ value: 1.5,
87+ min: 0.5,
88+ max: 2.5,
89+ step: 0.05,
90+ hint: 'The room is 2×side across — a real room, not a fraction of the domain.',
91+ },
92+ {
93+ key: 'gap',
94+ label: 'doorway (m)',
95+ value: 0.9,
96+ min: 0,
97+ max: 2,
98+ step: 0.05,
99+ hint: 'Width of the opening in the right-hand wall — 0.9 m is a typical interior door. At 0 the room is sealed and the sound never leaves.',
100+ },
101+ { key: 'thick', label: 'wall thickness (m)', value: 0.1, min: 0.03, max: 0.3, step: 0.01 },
102+ {
103+ key: 'absorb',
104+ label: 'wall absorption (1/s)',
105+ value: 70,
106+ min: 0,
107+ max: 1500,
108+ step: 10,
109+ hint: 'Damping inside the wall. It swallows whatever gets in rather than letting it rattle around in there, and it is what sets how long the room rings. At 0 a sealed room rings almost forever.',
110+ },
111+ ],
112+ // A click from a point inside the room, off-centre so it does not excite
113+ // only the symmetric modes, heard from the far corner. 700 Hz on a 1.5 m
114+ // half-width room is a wavelength of about half a metre — small enough to
115+ // show interference between the direct sound and the walls, and (at the
116+ // app's default 512-point grid over a 10 m domain) resolved to about 25
117+ // cells per wavelength.
118+ suggest: { point: 1, x0: -0.7, y0: 0.5, w: 0.03, f: 700, tw: 0.003, t0: 0.015, cw: 0 },
119+ mic: { x: 0.9, y: -0.7 },
120+ source: roomSource,
121+};
122+
123+const slit: MScene = {
124+ key: 'slit',
125+ label: 'Two slits',
126+ blurb: 'A hard screen with two apertures: diffraction and interference.',
127+ params: [
128+ {
129+ key: 'cwall',
130+ label: 'wall speed (×c0)',
131+ value: 4,
132+ min: 1.5,
133+ max: 8,
134+ step: 0.1,
135+ hint: 'Ratio to the background speed.',
136+ },
137+ { key: 'gap', label: 'slit width (m)', value: 0.4, min: 0.05, max: 2, step: 0.02 },
138+ { key: 'sep', label: 'slit separation (m)', value: 2, min: 0.5, max: 4, step: 0.05 },
139+ { key: 'thick', label: 'screen thickness (m)', value: 0.3, min: 0.05, max: 1, step: 0.02 },
140+ ],
141+ suggest: { point: 0, f: 400, tw: 0.008, t0: 0.04, cw: 0, x0: -3 },
142+ source: slitSource,
143+};
144+
145+const lens: MScene = {
146+ key: 'lens',
147+ label: 'Lens',
148+ blurb: 'A smooth slow patch that refracts a plane wave to a focus.',
149+ params: [
150+ {
151+ key: 'dn',
152+ label: 'speed dip (fraction of c0)',
153+ value: 0.35,
154+ min: 0,
155+ max: 0.8,
156+ step: 0.01,
157+ hint: 'Fraction of the background speed the centre of the lens is slower by. Deeper means a shorter focal length.',
158+ },
159+ { key: 'w', label: 'lens width (m)', value: 1.5, min: 0.5, max: 3, step: 0.05 },
160+ { key: 'x0', label: 'lens x (m)', value: -1, min: -3, max: 3, step: 0.1 },
161+ ],
162+ suggest: { point: 0, f: 400, tw: 0.008, t0: 0.04, cw: 0, x0: -3 },
163+ source: lensSource,
164+};
165+
166+const speckle: MScene = {
167+ key: 'speckle',
168+ label: 'Random medium',
169+ blurb: 'Weak random structure everywhere: multiple scattering, and a coda.',
170+ params: [
171+ { key: 'amp', label: 'contrast (fraction of c0)', value: 0.15, min: 0, max: 0.4, step: 0.01 },
172+ {
173+ key: 'kc',
174+ label: 'structure scale (rad/m)',
175+ value: 6,
176+ min: 1,
177+ max: 20,
178+ step: 0.5,
179+ hint: 'Centre wavenumber (2π / wavelength) of the random field. Scattering is strongest when this is comparable to the wave’s own.',
180+ },
181+ { key: 'seed', label: 'seed', value: 1, min: 1, max: 9999, step: 1, reseed: true },
182+ ],
183+ suggest: { point: 0, f: 400, tw: 0.008, t0: 0.04, cw: 0, x0: -3 },
184+ source: speckleSource,
185+};
186+
187+// The room first: it is the one that shows the most, and the one the
188+// microphone is for.
189+export const mScenes: MScene[] = [room, disk, slit, lens, speckle];
190+
191+export const mSceneByKey = (key: string): MScene | undefined =>
192+ mScenes.find((s) => s.key === key);
193+
194+export const defaultSceneParams = (s: MScene): Params =>
195+ Object.fromEntries(s.params.map((p) => [p.key, p.value]));
src/scene/scene.tsadded+212−0View file
@@ -0,0 +1,212 @@
1+/**
2+ * The medium: a .m scene file, evaluated once on the grid.
3+ *
4+ * A scene file is ordinary MATLAB defining one function,
5+ *
6+ * function [c, sig] = medium(x, y, <parameters>)
7+ *
8+ * over the solver's grid points. Unlike the models it is *not* compiled to
9+ * WGSL: a model's step runs every frame and must lower to a fixed sequence of
10+ * GPU dispatches, but a scene is evaluated exactly once and survives only as
11+ * two arrays of numbers. So it runs through numbl's CPU interpreter instead,
12+ * which buys the full MATLAB subset — loops, `if`, reductions, indexing,
13+ * seeded randomness via `rng`/`randn`, and anything in tools/ — and f64
14+ * evaluation, where the step dialect is element-wise f32.
15+ *
16+ * `c` is the sound speed in metres per second and `sig` the absorption rate in
17+ * inverse seconds. Both are ordinary
18+ * fields of position, which is what lets one function describe both the
19+ * scatterer and the open boundary: the absorbing layer every scene puts around
20+ * the edge (tools/sponge.m) is just the statement that the medium swallows
21+ * sound out there. A scene is free to put absorption inside the domain too,
22+ * which makes a lossy scatterer.
23+ */
24+import { parseMFile, type FunctionStmt } from 'numbl-src/numbl-core/parser/index.ts';
25+import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
26+import {
27+ RuntimeTensor,
28+ isRuntimeTensor,
29+ type RuntimeValue,
30+} from 'numbl-src/numbl-core/runtime/types.ts';
31+import type { Grid } from '../grid.ts';
32+import { C_AIR } from '../units.ts';
33+import { toolFiles } from '../tools.ts';
34+import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
35+import type { ModelParams } from '../mgpu/model.ts';
36+
37+/** The function a scene file must define. */
38+export const MEDIUM_FN = 'medium';
39+
40+export interface SceneOptions {
41+ grid: Grid;
42+ /** Scene source (.m text). */
43+ source: string;
44+ /** Parameter names the .m may take beyond `x` and `y`. */
45+ paramNames: string[];
46+ params: ModelParams;
47+}
48+
49+export class Scene {
50+ /** Sound speed on the grid, npts. */
51+ readonly c: Float32Array;
52+ /** Absorption rate on the grid, npts. */
53+ readonly sig: Float32Array;
54+ readonly cmin: number;
55+ readonly cmax: number;
56+ /** The background speed, taken to be whatever it is in the corner of the
57+ * domain — which is inside the absorbing layer, where a scene has no
58+ * business putting a scatterer. What the renderer shades departures from. */
59+ readonly cref: number;
60+ /** The largest departure from `cref` anywhere, so the renderer's wash has a
61+ * scale. Zero for a uniform medium, which draws nothing. */
62+ readonly cdev: number;
63+
64+ private constructor(c: Float32Array, sig: Float32Array) {
65+ this.c = c;
66+ this.sig = sig;
67+ let lo = Infinity;
68+ let hi = 0;
69+ for (const v of c) {
70+ if (v < lo) lo = v;
71+ if (v > hi) hi = v;
72+ }
73+ this.cmin = lo;
74+ this.cmax = hi;
75+ this.cref = c[0];
76+ this.cdev = Math.max(Math.abs(hi - this.cref), Math.abs(this.cref - lo));
77+ }
78+
79+ static create(opts: SceneOptions): Scene {
80+ const { grid, source, paramNames, params } = opts;
81+ const [c, sig] = evaluateMedium(source, paramNames, params, grid);
82+ for (let i = 0; i < c.length; i++) {
83+ if (!(c[i] > 0)) {
84+ throw new ModelCompileError(
85+ `the scene's sound speed is ${c[i]} somewhere; it must be positive ` +
86+ `everywhere (the timestep is set by the fastest point, and a zero ` +
87+ `or negative speed has no wave equation)`,
88+ { fn: MEDIUM_FN },
89+ );
90+ }
91+ if (!(sig[i] >= 0)) {
92+ throw new ModelCompileError(
93+ `the scene's absorption is ${sig[i]} somewhere; it must be zero or ` +
94+ `positive (a negative one would amplify rather than absorb)`,
95+ { fn: MEDIUM_FN },
96+ );
97+ }
98+ }
99+ return new Scene(c, sig);
100+ }
101+}
102+
103+/**
104+ * Evaluate the scene file on the grid, through numbl's CPU interpreter.
105+ *
106+ * The .m keeps the same contract a model has: it names the arguments it wants
107+ * — `x`, `y`, and any of the registry's parameters — and the host supplies
108+ * them by name, so their order in the signature is the .m's own business. A
109+ * one-line driver script calls `medium` with exactly the arguments its
110+ * signature declares, with those names pre-bound in the driver's workspace.
111+ */
112+function evaluateMedium(
113+ source: string,
114+ paramNames: string[],
115+ params: ModelParams,
116+ grid: Grid,
117+): [Float32Array, Float32Array] {
118+ const file = `${MEDIUM_FN}.m`;
119+ const ast = inModel(() => parseMFile(source, file));
120+ const fn = ast.body.find(
121+ (s): s is FunctionStmt =>
122+ s.type === 'Function' && (s as FunctionStmt).name === MEDIUM_FN,
123+ );
124+ if (!fn) {
125+ throw new ModelCompileError(`the scene defines no function named '${MEDIUM_FN}'`);
126+ }
127+ if (fn.outputs.length !== 2) {
128+ throw new ModelCompileError(
129+ `'${MEDIUM_FN}' must return two outputs [c, sig] — the sound speed and ` +
130+ `the absorption rate — not ${fn.outputs.length}`,
131+ { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
132+ );
133+ }
134+ // What the grid offers a scene by name, beyond its own parameters: the
135+ // coordinates themselves, in metres, and the numbers that describe them. A
136+ // scene wants `L` to place its absorbing layer relative to the domain and
137+ // `h` to smooth an interface over about a cell, so both are part of the
138+ // contract rather than something every scene has to be told. `c0` is the
139+ // speed of sound in air, which is what "the background" means here.
140+ const vars: Record<string, RuntimeValue> = {
141+ x: new RuntimeTensor(grid.x64, [grid.npts, 1]),
142+ y: new RuntimeTensor(grid.y64, [grid.npts, 1]),
143+ L: grid.L,
144+ h: grid.h,
145+ c0: C_AIR,
146+ npts: grid.npts,
147+ nx: grid.nx,
148+ ny: grid.ny,
149+ };
150+ const known = new Set([...Object.keys(vars), ...paramNames]);
151+ for (const p of fn.params) {
152+ if (!known.has(p)) {
153+ throw new ModelCompileError(
154+ `'${MEDIUM_FN}' takes an argument '${p}' that is neither the grid ` +
155+ `(${Object.keys(vars).join(', ')}) nor one of this scene's parameters` +
156+ (paramNames.length ? ` (${paramNames.join(', ')})` : ''),
157+ { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
158+ );
159+ }
160+ }
161+
162+ for (const name of paramNames) {
163+ const v = params[name];
164+ // Missing parameters read as 0, as ModelPlan.setParams has it.
165+ vars[name] = Number.isFinite(v) ? v : 0;
166+ }
167+
168+ const driver = `[c__, sig__] = ${MEDIUM_FN}(${fn.params.join(', ')});`;
169+ const result = inFunction(MEDIUM_FN, () =>
170+ executeCode(
171+ driver,
172+ { initialVariableValues: vars, displayResults: false, implicitCwdPath: null },
173+ [...toolFiles, { name: file, source }],
174+ 'scene-driver.m',
175+ ),
176+ );
177+
178+ return [
179+ toGridField(result.variableValues['c__'], fn.outputs[0], grid.npts),
180+ toGridField(result.variableValues['sig__'], fn.outputs[1], grid.npts),
181+ ];
182+}
183+
184+/** One returned field -> npts values, rounded to the solver's f32. */
185+function toGridField(
186+ value: RuntimeValue | undefined,
187+ name: string,
188+ npts: number,
189+): Float32Array {
190+ // A uniform field stays scalar in MATLAB; spread it over the grid.
191+ if (typeof value === 'number') return new Float32Array(npts).fill(value);
192+ if (value !== undefined && isRuntimeTensor(value)) {
193+ if (value.imag) {
194+ throw new ModelCompileError(
195+ `the scene's '${name}' is complex; the medium must be real`,
196+ { fn: MEDIUM_FN },
197+ );
198+ }
199+ // A vector of npts values, either orientation. A 2-D reshape is refused
200+ // rather than reordered: the tensor's column-major layout would not match
201+ // the grid's x-fastest rows.
202+ if (value.data.length === npts && value.shape.every((d) => d === 1 || d === npts)) {
203+ return new Float32Array(value.data);
204+ }
205+ throw new ModelCompileError(
206+ `the scene's '${name}' is ${value.shape.join(' x ')}, but the grid wants ` +
207+ `one value per point (${npts} x 1)`,
208+ { fn: MEDIUM_FN },
209+ );
210+ }
211+ throw new ModelCompileError(`the scene's '${name}' is not numeric`, { fn: MEDIUM_FN });
212+}
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 scene is evaluated by numbl's interpreter (see src/scene/scene.ts), which
5+ * resolves a call like `sponge(...)` against the workspace files it is handed.
6+ * Everything in `tools/` is handed to every such run, so any .m can call any
7+ * tool by name — MATLAB's own path semantics, where the file name is the
8+ * 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));
src/units.tsadded+41−0View file
@@ -0,0 +1,41 @@
1+/**
2+ * Everything in SI: metres, seconds, hertz, metres per second.
3+ *
4+ * The app used to be dimensionless — a domain two units across, a background
5+ * speed of one — which is tidy but leaves every number needing a translation
6+ * before it means anything, and leaves the microphone's recording with no
7+ * honest playback rate. In physical units all of that falls out: a wavelength
8+ * is a length you can compare to the room, a timestep is a real duration, and
9+ * the recorded trace plays back in real time at the pitch a microphone there
10+ * would have heard.
11+ *
12+ * It also makes the method's limits visible rather than hidden. A grid solver
13+ * resolves a wavelength with some number of cells, so a fixed grid over a
14+ * fixed domain is a low-frequency method: at 512 points across ten metres,
15+ * 700 Hz is 25 cells per wavelength and 2 kHz is nine. That ratio is now on
16+ * screen, because it is the number that decides whether what you are watching
17+ * is physics or grid dispersion.
18+ */
19+
20+/** Speed of sound in air at about 20 °C, m/s. */
21+export const C_AIR = 343;
22+
23+/** Side of the square domain, metres. Hall-sized: big enough for a room with
24+ * air around it, small enough that a wavefront crosses it in 30 ms. */
25+export const DOMAIN = 10;
26+
27+/** Cells per wavelength below which what is on screen is as much grid
28+ * dispersion as it is sound. */
29+export const POOR_RESOLUTION = 8;
30+
31+/** A length in metres, written the way a person would say it. */
32+export const fmtLength = (m: number): string =>
33+ Math.abs(m) < 1 ? `${(1000 * m).toPrecision(3)} mm` : `${m.toPrecision(3)} m`;
34+
35+/** A duration in seconds, likewise. */
36+export const fmtTime = (s: number): string => {
37+ const a = Math.abs(s);
38+ if (a > 0 && a < 1e-3) return `${(1e6 * s).toPrecision(3)} µs`;
39+ if (a < 1) return `${(1e3 * s).toPrecision(3)} ms`;
40+ return `${s.toPrecision(3)} s`;
41+};
test/checks.tsadded+656−0View file
@@ -0,0 +1,656 @@
1+/**
2+ * What the solver is held to.
3+ *
4+ * These run against the real pipeline — MATLAB source, numbl lowering,
5+ * generated WGSL, GPU — so they check the whole chain rather than any one
6+ * piece of it. The checks are physical wherever they can be: a wave should
7+ * travel at the speed the medium says, a scatterer that matches its
8+ * background should not scatter, an absorbing layer should absorb. What is
9+ * left over is the discretization's error, and that is what the numbers here
10+ * bound.
11+ *
12+ * Everything is in SI, at the app's own domain size (`DOMAIN`, 10 m) and
13+ * background speed (`C_AIR`, 343 m/s) — not a separate "toy" scale — so a
14+ * number that appears here means the same thing it would in the app, and the
15+ * synthetic scenes below (`uniform`, `closedBox`) share the app's own
16+ * absorbing-layer profile rather than inventing their own.
17+ *
18+ * Nothing here checks what the app looks like. That is for a browser.
19+ */
20+import { ModelSession } from '../src/mgpu/session.ts';
21+import type { MModel } from '../src/mgpu/registry.ts';
22+import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';
23+import type { MScene } from '../src/scene/registry.ts';
24+import { mSceneByKey, defaultSceneParams } from '../src/scene/registry.ts';
25+import { C_AIR, DOMAIN } from '../src/units.ts';
26+
27+export type Check = (name: string, ok: boolean, detail: string) => void;
28+export type Log = (s: string) => void;
29+
30+/** A homogeneous medium with the usual absorbing edge — the case every
31+ * scattering result is measured against. Same sponge profile as every real
32+ * scene (src/tools/sponge.m via 0.2*L, 1700), so its absorption behaviour is
33+ * exactly the app's, not a separately tuned stand-in. */
34+const uniform: MScene = {
35+ key: 'uniform',
36+ label: 'Uniform',
37+ blurb: 'No scatterer at all.',
38+ params: [],
39+ source: `
40+function [c, sig] = medium(x, y, L, c0)
41+ c = c0 + 0*x;
42+ sig = sponge(x, y, L, 0.2*L, 1700);
43+end
44+`,
45+};
46+
47+/** The same, with nothing absorbing anywhere: a closed box. */
48+const closedBox: MScene = {
49+ ...uniform,
50+ key: 'closed',
51+ source: `
52+function [c, sig] = medium(x, y, c0)
53+ c = c0 + 0*x;
54+ sig = 0*x;
55+end
56+`,
57+};
58+
59+/** A model whose "step" is one application of a Laplacian stencil, so the
60+ * stencil can be measured on a field we chose. */
61+const stencilProbe = (op: 'lap2' | 'lap4'): MModel => ({
62+ key: `probe-${op}`,
63+ label: `probe ${op}`,
64+ blurb: '',
65+ state: ['p', 'pm', 't'],
66+ params: [],
67+ order: op === 'lap2' ? 2 : 4,
68+ source: `
69+function [p, pm, t] = init(npts)
70+ p = zeros(npts, 1);
71+ pm = zeros(npts, 1);
72+ t = zeros(npts, 1);
73+end
74+
75+function [pn, pold, tn] = step(p, pm, t)
76+ pn = ${op}(p);
77+ pold = pm;
78+ tn = t;
79+end
80+`,
81+});
82+
83+const maxAbs = (a: Float32Array): number => {
84+ let m = 0;
85+ for (const v of a) m = Math.max(m, Math.abs(v));
86+ return m;
87+};
88+
89+/**
90+ * The Laplacian stencils, against a field whose Laplacian is known exactly.
91+ *
92+ * p = sin(kx*x) * sin(ky*y) has lap(p) = -(kx^2 + ky^2) * p. Both stencils
93+ * should reproduce that away from the boundary, the 5-point one to O((k*h)^2)
94+ * and the 9-point one to O((k*h)^4) — which at the resolution used here is
95+ * two orders of magnitude tighter. The margin matters more than either
96+ * number: it is what says the fourth-order stencil is actually fourth order
97+ * and not a mistyped second-order one. Purely a statement about the discrete
98+ * operator, so it does not depend on the medium or the domain size at all.
99+ */
100+export async function stencilChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
101+ const n = 64;
102+ const kx = 4 * Math.PI;
103+ const ky = 3 * Math.PI;
104+ const errs: Record<string, number> = {};
105+
106+ for (const op of ['lap2', 'lap4'] as const) {
107+ const session = await ModelSession.create({
108+ device,
109+ model: stencilProbe(op),
110+ params: {},
111+ scene: uniform,
112+ sceneParams: {},
113+ n,
114+ });
115+ try {
116+ const { grid } = session;
117+ const p = new Float32Array(grid.npts);
118+ for (let i = 0; i < grid.npts; i++) {
119+ p[i] = Math.sin(kx * grid.x64[i]) * Math.sin(ky * grid.y64[i]);
120+ }
121+ session.reset();
122+ session.gpu.upload('p', p);
123+ session.step(1);
124+ const got = await session.read('p');
125+
126+ // Interior only: the stencil takes the field outside the grid to be
127+ // zero, which is not what this analytic field does.
128+ const pad = 3;
129+ const want = -(kx * kx + ky * ky);
130+ let err = 0;
131+ for (let iy = pad; iy < n - pad; iy++) {
132+ for (let ix = pad; ix < n - pad; ix++) {
133+ const i = ix + n * iy;
134+ err = Math.max(err, Math.abs(got[i] - want * p[i]));
135+ }
136+ }
137+ errs[op] = err / Math.abs(want);
138+ log(` ${op}: max relative error ${errs[op].toExponential(2)} at n = ${n}`);
139+ } finally {
140+ session.destroy();
141+ }
142+ }
143+
144+ check(
145+ 'lap2 matches the analytic Laplacian',
146+ errs.lap2 < 0.02,
147+ `relative error ${errs.lap2.toExponential(2)} (expect ~(k*h)^2/12)`,
148+ );
149+ check(
150+ 'lap4 is far more accurate than lap2',
151+ errs.lap4 < errs.lap2 / 10,
152+ `${errs.lap4.toExponential(2)} vs ${errs.lap2.toExponential(2)}`,
153+ );
154+}
155+
156+/**
157+ * A pulse from a point source should be a ring of radius c*(t - t0).
158+ *
159+ * This is the end-to-end statement that the thing solves the wave equation:
160+ * it exercises the source term, the model's own clock, the stencil, and the
161+ * timestep the host computed, and it fails if any of them is wrong by a
162+ * constant factor. 600 Hz on this grid (256 points over 10 m, h = 3.9 cm) is
163+ * about 15 cells per wavelength — comfortably resolved, and short enough that
164+ * a clean, narrow pulse is cheap to run.
165+ */
166+export async function propagationChecks(
167+ device: GPUDevice,
168+ check: Check,
169+ log: Log,
170+): Promise<void> {
171+ const model = mModelByKey('leapfrog')!;
172+ const t0 = 0.005;
173+ const session = await ModelSession.create({
174+ device,
175+ model,
176+ params: {
177+ ...defaultParams(model),
178+ f: 600,
179+ tw: 0.0012,
180+ t0,
181+ cw: 0,
182+ point: 1,
183+ x0: 0,
184+ y0: 0,
185+ w: 0.15,
186+ },
187+ scene: uniform,
188+ sceneParams: {},
189+ n: 512,
190+ L: DOMAIN,
191+ });
192+ try {
193+ session.reset();
194+ // Early enough that the front is still clear of the absorbing layer,
195+ // which starts at |x| = 0.3*DOMAIN and would pull the peak back towards
196+ // the interior.
197+ const until = 0.01;
198+ const steps = Math.round(until / session.dt);
199+ session.step(steps);
200+ const p = await session.read('p');
201+ const t = session.steps * session.dt;
202+
203+ // Where the wavefront is, along +x from the source at the origin. An
204+ // energy centroid rather than the bare peak: at 600 Hz the wavelength is
205+ // 0.57 m, so the tallest individual fringe of an oscillating pulse can
206+ // sit anywhere within half a wavelength of the envelope's true centre
207+ // depending on carrier phase, which swamps the grid's own ~1% resolution.
208+ // Weighting position by p^2 averages over the fringes instead of picking
209+ // whichever one happens to be tallest.
210+ const { grid } = session;
211+ const cutoff = 0.275 * grid.L; // stays clear of the sponge, which starts at 0.3*L
212+ const iy = Math.floor(grid.ny / 2);
213+ let weighted = 0;
214+ let weight = 0;
215+ for (let ix = Math.floor(grid.nx / 2); ix < grid.nx; ix++) {
216+ const i = ix + grid.nx * iy;
217+ if (grid.x64[i] > cutoff) break;
218+ const w2 = p[i] * p[i];
219+ weighted += grid.x64[i] * w2;
220+ weight += w2;
221+ }
222+ const best = weighted / weight;
223+ const want = C_AIR * (t - t0);
224+ const err = Math.abs(best - want) / want;
225+ log(` wavefront at r = ${best.toFixed(4)} m, expected ${want.toFixed(4)} m at t = ${(1000 * t).toFixed(3)} ms`);
226+ check(
227+ 'a pulse travels at the medium speed',
228+ err < 0.03,
229+ `radius off by ${(100 * err).toFixed(1)}%`,
230+ );
231+ } finally {
232+ session.destroy();
233+ }
234+}
235+
236+/**
237+ * The absorbing layer should leave next to nothing behind.
238+ *
239+ * A plane pulse is launched, crosses the grid, and is swallowed. What is
240+ * still in the interior long afterwards is what the sponge reflected, and it
241+ * is the honest measure of how open the open boundary is.
242+ */
243+export async function boundaryChecks(
244+ device: GPUDevice,
245+ check: Check,
246+ log: Log,
247+): Promise<void> {
248+ const model = mModelByKey('leapfrog')!;
249+ const session = await ModelSession.create({
250+ device,
251+ model,
252+ params: {
253+ ...defaultParams(model),
254+ f: 300, tw: 0.006, t0: 0.03, cw: 0, point: 0, x0: -0.25 * DOMAIN,
255+ },
256+ scene: uniform,
257+ sceneParams: {},
258+ n: 256,
259+ L: DOMAIN,
260+ });
261+ try {
262+ session.reset();
263+ const stepsTo = (t: number): number => Math.round(t / session.dt);
264+ session.step(stepsTo(0.07));
265+ const peak = maxAbs(await session.read('p'));
266+
267+ session.step(stepsTo(0.3) - session.steps);
268+ const after = await session.read('p');
269+
270+ // The interior only: the sponge itself is allowed to hold whatever it is
271+ // busy absorbing.
272+ const { grid } = session;
273+ const interior = 0.3 * grid.L; // exactly where the sponge starts
274+ let residual = 0;
275+ for (let i = 0; i < grid.npts; i++) {
276+ if (Math.abs(grid.x64[i]) < interior && Math.abs(grid.y64[i]) < interior) {
277+ residual = Math.max(residual, Math.abs(after[i]));
278+ }
279+ }
280+ const ratio = residual / peak;
281+ log(` peak ${peak.toExponential(2)}, interior residual at t = 0.3 s is ${ratio.toExponential(2)} of it`);
282+ check(
283+ 'the absorbing layer reflects little',
284+ ratio < 0.02,
285+ `residual ${(100 * ratio).toFixed(2)}% of the incident peak`,
286+ );
287+ } finally {
288+ session.destroy();
289+ }
290+}
291+
292+/**
293+ * A scatterer whose speed matches the background is not a scatterer.
294+ *
295+ * Running the disk scene at cin = 1 must reproduce the uniform medium
296+ * exactly, which is a strong statement about the whole scene path: the
297+ * smoothed interface, the coordinates, the upload. And at cin = 3 there must
298+ * be a scattered field worth looking at, or the app would be drawing nothing.
299+ * Uses the real `disk` scene and its real domain, so what is checked is
300+ * exactly what the app runs.
301+ */
302+export async function scatteringChecks(
303+ device: GPUDevice,
304+ check: Check,
305+ log: Log,
306+): Promise<void> {
307+ const model = mModelByKey('leapfrog')!;
308+ const disk = mSceneByKey('disk')!;
309+ const params = {
310+ ...defaultParams(model),
311+ f: 400,
312+ tw: 0.002,
313+ t0: 0.008,
314+ cw: 0,
315+ point: 0,
316+ x0: -0.25 * DOMAIN,
317+ };
318+ const n = 256;
319+ const run = async (scene: MScene, sceneParams: Record<string, number>): Promise<Float32Array> => {
320+ const session = await ModelSession.create({
321+ device, model, params, scene, sceneParams, n, L: DOMAIN,
322+ });
323+ try {
324+ session.reset();
325+ session.step(Math.round(0.022 / session.dt));
326+ return await session.read('p');
327+ } finally {
328+ session.destroy();
329+ }
330+ };
331+
332+ const plain = await run(uniform, {});
333+ const matched = await run(disk, { ...defaultSceneParams(disk), cin: 1, absorb: 0 });
334+ const hard = await run(disk, { ...defaultSceneParams(disk), cin: 3, absorb: 0 });
335+
336+ const peak = maxAbs(plain);
337+ let dMatched = 0;
338+ let dHard = 0;
339+ for (let i = 0; i < plain.length; i++) {
340+ dMatched = Math.max(dMatched, Math.abs(matched[i] - plain[i]));
341+ dHard = Math.max(dHard, Math.abs(hard[i] - plain[i]));
342+ }
343+ log(` matched disk differs by ${(dMatched / peak).toExponential(2)}, hard disk by ${(dHard / peak).toFixed(2)}`);
344+ check(
345+ 'a speed-matched disk does not scatter',
346+ dMatched / peak < 1e-3,
347+ `scattered field ${(dMatched / peak).toExponential(2)} of the incident peak`,
348+ );
349+ check(
350+ 'a hard disk scatters strongly',
351+ dHard / peak > 0.2,
352+ `scattered field ${(100 * dHard / peak).toFixed(0)}% of the incident peak`,
353+ );
354+}
355+
356+/**
357+ * The timestep the host picks should be stable, and near the edge of being
358+ * unstable — a scheme that is merely stable because it is crawling is not
359+ * evidence of anything. Run a closed box (no absorption at all, so nothing
360+ * can hide a slow instability) at 95% of the computed limit and watch it
361+ * bounce around for a long time.
362+ */
363+export async function stabilityChecks(
364+ device: GPUDevice,
365+ check: Check,
366+ log: Log,
367+): Promise<void> {
368+ for (const key of ['leapfrog', 'leapfrog4']) {
369+ const model = mModelByKey(key)!;
370+ const session = await ModelSession.create({
371+ device,
372+ model,
373+ params: {
374+ ...defaultParams(model),
375+ f: 300, tw: 0.003, t0: 0.01, cw: 0, point: 1, x0: 0, y0: 0,
376+ },
377+ scene: closedBox,
378+ sceneParams: {},
379+ n: 128,
380+ L: DOMAIN,
381+ cfl: 0.95,
382+ });
383+ try {
384+ session.reset();
385+ session.step(Math.round(0.03 / session.dt));
386+ const early = maxAbs(await session.read('p'));
387+ session.step(Math.round(0.6 / session.dt));
388+ const late = maxAbs(await session.read('p'));
389+ log(` ${key}: max|p| ${early.toExponential(2)} at t = 30 ms, ${late.toExponential(2)} at t = 630 ms`);
390+ check(
391+ `${key} is stable at 95% of the CFL limit`,
392+ Number.isFinite(late) && late < 5 * early,
393+ `max|p| went from ${early.toExponential(2)} to ${late.toExponential(2)} over 600 ms`,
394+ );
395+ } finally {
396+ session.destroy();
397+ }
398+ }
399+}
400+
401+/**
402+ * The planner's kernel splitting must not change the answer.
403+ *
404+ * On any device worth running this on, the leapfrog update fits in one
405+ * kernel. Squeeze the budget down to two grid fields per kernel — what a
406+ * compatibility-mode device would allow — and the same line has to be
407+ * evaluated in half a dozen pieces through scratch buffers. The arithmetic is
408+ * the same; only the rounding of the intermediates differs, since each piece
409+ * is stored as f32 on the way out.
410+ */
411+export async function splitChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
412+ const model = mModelByKey('leapfrog')!;
413+ const disk = mSceneByKey('disk')!;
414+ const params = {
415+ ...defaultParams(model),
416+ f: 400, tw: 0.002, t0: 0.008, point: 0, x0: -0.25 * DOMAIN,
417+ };
418+ const run = async (operandBudget?: number): Promise<{ p: Float32Array; ops: string[] }> => {
419+ const session = await ModelSession.create({
420+ device,
421+ model,
422+ params,
423+ scene: disk,
424+ sceneParams: defaultSceneParams(disk),
425+ n: 128,
426+ L: DOMAIN,
427+ operandBudget,
428+ });
429+ try {
430+ session.reset();
431+ session.step(Math.round(0.02 / session.dt));
432+ return { p: await session.read('p'), ops: session.describe().step };
433+ } finally {
434+ session.destroy();
435+ }
436+ };
437+
438+ const whole = await run();
439+ const split = await run(2);
440+ const peak = maxAbs(whole.p);
441+ let diff = 0;
442+ for (let i = 0; i < whole.p.length; i++) {
443+ diff = Math.max(diff, Math.abs(whole.p[i] - split.p[i]));
444+ }
445+ log(` ${whole.ops.length} ops whole, ${split.ops.length} split; largest difference ${(diff / peak).toExponential(2)}`);
446+ check(
447+ 'splitting a kernel does not change the answer',
448+ diff / peak < 1e-3,
449+ `fields differ by ${(diff / peak).toExponential(2)} of the peak`,
450+ );
451+ check(
452+ 'a squeezed budget really does split the update',
453+ split.ops.length > whole.ops.length,
454+ `${split.ops.length} ops vs ${whole.ops.length}`,
455+ );
456+}
457+
458+/**
459+ * The microphone records the field, at the point it is pointed at, once per
460+ * timestep.
461+ *
462+ * Checked against the field itself rather than against a description of it:
463+ * after N steps the trace must be N samples long, and its last sample must be
464+ * exactly — not approximately — the pressure sitting at the probe's grid point,
465+ * since both are the same f32 written by the same kernel. That pins the probe
466+ * index, the grid layout, and the fact that the recording dispatch really does
467+ * run once per step rather than once per submission.
468+ */
469+export async function microphoneChecks(
470+ device: GPUDevice,
471+ check: Check,
472+ log: Log,
473+): Promise<void> {
474+ const model = mModelByKey('leapfrog')!;
475+ const session = await ModelSession.create({
476+ device,
477+ model,
478+ params: {
479+ ...defaultParams(model),
480+ f: 600, tw: 0.0012, t0: 0.003, cw: 0, point: 1, x0: 0, y0: 0, w: 0.15,
481+ },
482+ scene: uniform,
483+ sceneParams: {},
484+ n: 128,
485+ L: DOMAIN,
486+ });
487+ try {
488+ const { grid } = session;
489+ const mx = 1.5;
490+ const my = 0.5;
491+ session.setMic(mx, my);
492+ session.reset();
493+ // Long enough for the wave to have reached the microphone and moved on.
494+ const steps = Math.round(0.01 / session.dt);
495+ session.step(steps);
496+
497+ const trace = await session.recorder.read();
498+ const field = await session.read('p');
499+ const ix = Math.round((mx + grid.L / 2) / grid.h - 0.5);
500+ const iy = Math.round((my + grid.L / 2) / grid.h - 0.5);
501+ const at = field[ix + grid.nx * iy];
502+
503+ let peak = 0;
504+ for (const v of trace) peak = Math.max(peak, Math.abs(v));
505+ log(` ${trace.length} samples in ${steps} steps, peak ${peak.toExponential(2)}, last ${trace[trace.length - 1].toExponential(3)} vs field ${at.toExponential(3)}`);
506+
507+ check(
508+ 'the microphone records one sample per timestep',
509+ trace.length === steps,
510+ `${trace.length} samples for ${steps} steps`,
511+ );
512+ check(
513+ 'the microphone records the field at its own grid point',
514+ trace.length > 0 && trace[trace.length - 1] === at,
515+ `last sample ${trace[trace.length - 1]} vs field ${at}`,
516+ );
517+ check(
518+ 'the microphone hears the wave arrive',
519+ peak > 1e-3,
520+ `peak |p| at the microphone was ${peak.toExponential(2)}`,
521+ );
522+
523+ // Moving it must move what it hears: the same run sampled at the origin,
524+ // where a point source is loudest, cannot match a point away from it.
525+ session.setMic(0, 0);
526+ session.reset();
527+ session.step(steps);
528+ const atSource = await session.recorder.read();
529+ let peak2 = 0;
530+ for (const v of atSource) peak2 = Math.max(peak2, Math.abs(v));
531+ log(` peak at the source ${peak2.toExponential(2)}, at r = ${Math.hypot(mx, my).toFixed(2)} m ${peak.toExponential(2)}`);
532+ check(
533+ 'moving the microphone changes what it hears',
534+ peak2 > peak,
535+ `${peak2.toExponential(2)} at the source vs ${peak.toExponential(2)} away from it`,
536+ );
537+ } finally {
538+ session.destroy();
539+ }
540+}
541+
542+/**
543+ * A room rings, and sealing it makes it ring for longer.
544+ *
545+ * Measured as energy in the second half of the microphone's trace against the
546+ * first — a ratio rather than an envelope, because a small room beats between
547+ * its modes and any one window can land in a null. In the open field the same
548+ * pulse passes the microphone once and is gone, which is the contrast that
549+ * makes the number mean something. Uses the real `room` scene's own default
550+ * geometry, at the app's real domain size.
551+ *
552+ * The source is tuned for this test's own grid rather than reused from the
553+ * scene's `suggest`. A wall slower than the background (`cwall` = 0.15) has a
554+ * *shorter* wavelength inside itself than the background does at the same
555+ * frequency, by that same factor — the app's own docs on this scene call this
556+ * out — and an unresolved wall does not behave like a partial reflector, it
557+ * behaves like an absorber: at 256 grid points it swallowed the whole pulse
558+ * in a handful of bounces regardless of the `absorb` parameter, making every
559+ * room in an early version of this check look identically "sealed" no matter
560+ * what. So this uses the app's own 512-point grid, where the wall's own
561+ * wavelength is resolved to about ten cells at 220 Hz rather than four at the
562+ * scene's demo frequency (700 Hz) — and even then the ring is real but not
563+ * long: measured, late/early lands around 0.15, well above the open field's
564+ * 0.004 but short of the naive "rings for a while" threshold a lossless room
565+ * would give. That is the wall's transmission loss actually doing its job,
566+ * not a bug — 74% amplitude reflection per bounce (cwall = 0.15 gives
567+ * |c-1|/(c+1) = 0.74) empties a small room in a few tens of bounces.
568+ */
569+export async function roomChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
570+ const model = mModelByKey('leapfrog')!;
571+ const room = mSceneByKey('room')!;
572+ const params = {
573+ ...defaultParams(model),
574+ point: 1, x0: -0.7, y0: 0.5, w: 0.03, f: 220, tw: 0.002, t0: 0.008, cw: 0,
575+ };
576+ const roomDefaults = defaultSceneParams(room);
577+ const mic = { x: roomDefaults.side * 0.6, y: -roomDefaults.side * 0.5 };
578+
579+ const listen = async (scene: MScene, sceneParams: Record<string, number>) => {
580+ const session = await ModelSession.create({
581+ device, model, params, scene, sceneParams, n: 512, L: DOMAIN,
582+ });
583+ try {
584+ session.setMic(mic.x, mic.y);
585+ session.reset();
586+ session.step(Math.round(0.15 / session.dt));
587+ const trace = await session.recorder.read();
588+ const half = Math.floor(trace.length / 2);
589+ let early = 0;
590+ let late = 0;
591+ for (let i = 0; i < half; i++) early += trace[i] * trace[i];
592+ for (let i = half; i < trace.length; i++) late += trace[i] * trace[i];
593+ return { early, late, ratio: late / Math.max(early, 1e-30) };
594+ } finally {
595+ session.destroy();
596+ }
597+ };
598+
599+ const open = await listen(uniform, {});
600+ const sealed = await listen(room, { ...roomDefaults, gap: 0 });
601+ const wide = await listen(room, { ...roomDefaults, gap: 2 * roomDefaults.side });
602+ log(` late/early energy — open field ${open.ratio.toExponential(3)}, sealed room ${sealed.ratio.toFixed(3)}, wide door ${wide.ratio.toFixed(3)}`);
603+
604+ check(
605+ 'a pulse in the open field does not come back',
606+ open.ratio < 0.02,
607+ `late/early energy ${open.ratio.toExponential(2)}`,
608+ );
609+ check(
610+ 'a room rings after the pulse has passed',
611+ sealed.ratio > 0.05,
612+ `late/early energy ${sealed.ratio.toFixed(3)} inside the room, vs ${open.ratio.toExponential(2)} in the open`,
613+ );
614+ check(
615+ 'sound leaves through the doorway',
616+ sealed.late > 1.2 * wide.late,
617+ `late energy ${sealed.late.toExponential(2)} sealed vs ${wide.late.toExponential(2)} with a wide door`,
618+ );
619+}
620+
621+/**
622+ * What a step compiles to. A guard on the fusion passes: if one of them stops
623+ * firing, the model still gives the right answer, only several times slower,
624+ * and nothing else here would notice.
625+ */
626+export async function planChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
627+ const model = mModelByKey('leapfrog')!;
628+ const disk = mSceneByKey('disk')!;
629+ const session = await ModelSession.create({
630+ device,
631+ model,
632+ params: defaultParams(model),
633+ scene: disk,
634+ sceneParams: defaultSceneParams(disk),
635+ n: 128,
636+ L: DOMAIN,
637+ });
638+ try {
639+ const ops = session.describe().step;
640+ for (const line of ops) log(` ${line}`);
641+ const kernels = ops.filter((o) => o.startsWith('kernel')).length;
642+ const stencils = ops.filter((o) => o.startsWith('stencil')).length;
643+ check(
644+ 'the step uses exactly one stencil dispatch',
645+ stencils === 1,
646+ `${stencils} stencil ops`,
647+ );
648+ check(
649+ 'the source term fuses into the update',
650+ kernels <= 6,
651+ `${kernels} element-wise kernels (5 expected: u, sd, pn, pold, tn)`,
652+ );
653+ } finally {
654+ session.destroy();
655+ }
656+}
tools/sponge.madded+22−0View file
@@ -0,0 +1,22 @@
1+% Absorption profile for an open boundary.
2+%
3+% s = sponge(x, y, L, w, smax)
4+%
5+% All in SI: positions and widths in metres, the result in inverse seconds.
6+% Zero in the interior, ramping up quadratically over a layer of width w
7+% inside each edge of the square [-L/2, L/2]^2 and reaching smax at the wall.
8+% Added to a scene's absorption, this is what makes the finite grid stand in
9+% for an unbounded medium: a wave that leaves the region of interest is
10+% attenuated before it reaches the outer boundary, and whatever reflects off
11+% that boundary is attenuated again on the way back.
12+%
13+% The ramp is gradual on purpose. An absorbing layer is itself an impedance
14+% mismatch, so a sudden one reflects; spreading it over a couple of
15+% wavelengths keeps that reflection small. It is not a perfectly matched
16+% layer, and at grazing incidence it does leak.
17+function s = sponge(x, y, L, w, smax)
18+ dx = max(0, w - (L/2 - abs(x))) / w;
19+ dy = max(0, w - (L/2 - abs(y))) / w;
20+ d = max(dx, dy);
21+ s = smax * d.^2;
22+end
tsconfig.jsonadded+15−0View file
@@ -0,0 +1,15 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "module": "ESNext",
5+ "moduleResolution": "bundler",
6+ "lib": ["ES2022", "DOM", "DOM.Iterable"],
7+ "types": ["@webgpu/types", "node"],
8+ "strict": true,
9+ "noEmit": true,
10+ "allowImportingTsExtensions": true,
11+ "verbatimModuleSyntax": true,
12+ "skipLibCheck": true
13+ },
14+ "include": ["src", "test", "scripts"]
15+}
vite.config.tsadded+28−0View file
@@ -0,0 +1,28 @@
1+import { defineConfig } from 'vite';
2+import { realpathSync } from 'node:fs';
3+import { resolve } from 'node:path';
4+
5+// numbl is a local `file:` dependency, so node_modules/numbl is a symlink to
6+// the sibling checkout. Its package `exports` map only publishes the runtime
7+// entry points, not the compiler internals we need (parser + JIT lowering), so
8+// we reach them through a path alias — the same arrangement turing-surface and
9+// math-webgpu-sandbox use.
10+// Realpath'd through the symlink: dev serves modules under their real ids, so
11+// aliasing the node_modules path would give the same file two identities (one
12+// per spelling) and run its side effects twice — the interpreter's builtin
13+// registry throws on the second.
14+const numblSrc = realpathSync(resolve(import.meta.dirname, 'node_modules/numbl/src'));
15+
16+export default defineConfig({
17+ base: './',
18+ resolve: {
19+ alias: { 'numbl-src': numblSrc },
20+ },
21+ server: {
22+ // the alias resolves outside the project root (through the symlink)
23+ fs: { allow: [import.meta.dirname, numblSrc] },
24+ },
25+ build: {
26+ target: 'es2022',
27+ },
28+});