A plucked dulcimer string and its box, as two coupled wave equations on WebGPU
44 changed files+10687−0
.github/workflows/deploy.ymladded+66−0View file
@@ -0,0 +1,66 @@
1+name: deploy
2+on:
3+ push:
4+ branches: [main]
5+ workflow_dispatch:
6+
7+permissions:
8+ contents: read
9+ pages: write
10+ id-token: write
11+
12+concurrency:
13+ group: pages
14+ cancel-in-progress: true
15+
16+jobs:
17+ build-deploy:
18+ runs-on: ubuntu-latest
19+ environment:
20+ name: github-pages
21+ url: ${{ steps.deployment.outputs.page_url }}
22+ steps:
23+ - uses: actions/checkout@v4
24+ - uses: actions/setup-node@v4
25+ with:
26+ node-version: 24
27+ cache: npm
28+ # numbl is a `file:../../numbl` dependency: we use its compiler internals
29+ # (parser, lowerer, IR, inline pass) and its interpreter, which its
30+ # published package `exports` do not expose. Clone it where that relative
31+ # path expects it. Pinned to the exact commit this project was built and
32+ # tested against, so CI builds exactly what was verified locally rather
33+ # than whatever numbl's main happens to be on the day it runs.
34+ #
35+ # numbl's own dependencies are NOT needed: the slice we import is
36+ # self-contained TypeScript, verified by building against a checkout with
37+ # no node_modules. One file in it is generated rather than committed,
38+ # though — the interpreter's stdlib bundle, which numbl gitignores — so a
39+ # bare checkout is missing it and `executeCode.ts` (used to evaluate scene
40+ # .m files on the CPU) fails to resolve it. Its generator only reads .m
41+ # files off disk, so plain `node` (type-stripping, unflagged since 22.18)
42+ # runs it without installing anything.
43+ - name: Check out numbl (sibling dependency)
44+ env:
45+ NUMBL_REF: 3901dbcbfa127240ee57a158184b2aa235478754
46+ run: |
47+ git clone --filter=blob:none --no-checkout \
48+ https://github.com/flatironinstitute/numbl.git "$GITHUB_WORKSPACE/../../numbl"
49+ git -C "$GITHUB_WORKSPACE/../../numbl" checkout --quiet "$NUMBL_REF"
50+ node "$GITHUB_WORKSPACE/../../numbl/scripts/bundle-stdlib.ts"
51+ # --ignore-scripts: npm runs a linked package's `prepare` script, and
52+ # numbl's is husky, which is not installed here.
53+ - run: npm ci --ignore-scripts
54+ # The checks compile MATLAB to compute shaders, so they need a GPU;
55+ # --skip-without-gpu lets this runner (which has none) say so and move on.
56+ - run: npm run test:node -- --skip-without-gpu
57+ - run: npm run build
58+ # Pages must already be enabled with "GitHub Actions" as the source; the
59+ # workflow token cannot create the site itself (`enablement: true` fails
60+ # with "Resource not accessible by integration").
61+ - uses: actions/configure-pages@v5
62+ - uses: actions/upload-pages-artifact@v3
63+ with:
64+ path: dist
65+ - id: deployment
66+ uses: actions/deploy-pages@v4
.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
README.mdadded+119−0View file
@@ -0,0 +1,119 @@
1+# dulcimer
2+
3+One plucked string and the box it sounds over, simulated as two coupled wave
4+equations and solved live in the browser. Both solvers are written as MATLAB,
5+run by [numbl](https://numbl.org) and compiled to WebGPU compute shaders; the
6+pressure field is drawn straight out of the buffer the solver writes, a
7+microphone records the air at one point every timestep, and the recording
8+plays back at the pitch a microphone there would have heard.
9+
10+This is the instrument-shaped sibling of
11+[acoustic-scattering-2d](https://github.com/concept-collection/acoustic-scattering-2d),
12+whose MATLAB-to-WGSL compiler pipeline it reuses, and of
13+[acoustic-scattering-3d](https://github.com/concept-collection/acoustic-scattering-3d),
14+whose volume renderer it adapts. What is new here is that two equations run
15+together, on two different grids, with the coupling between them part of the
16+editable model.
17+
18+## The equations
19+
20+The string carries transverse displacement u (metres) on a line of nodes:
21+
22+```
23+u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
24+```
25+
26+cs = 2 Ls f0 is the wave speed a string of length Ls needs to sound the
27+fundamental f0. The fourth-derivative term is bending stiffness, expressed
28+through the inharmonicity coefficient B (partial n lands near
29+n f0 sqrt(1 + B n^2)); sig0 is plain decay, quoted as a 60 dB time; sig1
30+damps high frequencies faster than low, which is why a plucked note starts
31+bright and mellows. The pluck is an initial condition: a triangle drawn to
32+the pluck point, released from rest. This is the standard stiff-string
33+formulation (Bilbao, Numerical Sound Synthesis, ch. 7), with the standard
34+explicit scheme.
35+
36+The air carries pressure p on a rectangular grid around the instrument:
37+
38+```
39+p_tt + 2 sig p_t = c^2 lapw(p, wall) + s
40+```
41+
42+`lapw` is a wall-masked Laplacian, and it is how the body is rigid. Each face
43+of the 7-point stencil is scaled by the mask at the neighbour it reads, so a
44+face into the shell carries no flux: the discrete Neumann (sound-hard)
45+condition, in the divergence form div(w grad p). The flat siblings make
46+walls out of fast material instead, which is unavailable here: wood at about
47+4000 m/s would cut the global timestep twelve-fold, where a mask costs
48+nothing. What a mask cannot do is absorb, so the scene paints a thin lossy
49+skin over the shell's surfaces, and that parameter is what tames the
50+cavity's ring.
51+
52+The coupling from string to air is one-way and comes in two routes, each with
53+its own gain. `string radiation` injects the string's acceleration along its
54+own line of cells: the direct route, idealized, since a thin string is in
55+reality a very poor radiator. `bridge drive` takes the string's pull at its
56+bridge end (the tension times the arriving slope) and drives a patch of air
57+just above the top plate: the real instrument's main route, minus the
58+plate's own resonances, which this model does not carry. Either gain at zero
59+switches that route off. The air does not push back on the string; at these
60+amplitudes the back-reaction is far below everything else the model already
61+neglects.
62+
63+## Two files, two grids
64+
65+As in the 2d sibling, the *model* (init and step, both grids together) is
66+compiled: numbl lowers it to typed IR and this project's backend emits one
67+WGSL kernel per source line, with five host-provided operations — `dxx`,
68+`dxxxx` on the string, `lapw` in the air, and `spread`/`bridge` carrying the
69+string into the air — as the only places anything reads a neighbour. The
70+*scene* (the body: walls, absorption, and the two coupling profiles) is
71+interpreted on the CPU once per edit, so it has the whole MATLAB subset
72+available and costs no recompile.
73+
74+The two grids are not independent. The air's CFL condition fixes the
75+timestep, and the string then chooses the finest node spacing that is stable
76+at that dt for anything the sliders can reach, so moving a parameter never
77+forces a recompile. In numbers, at the default 128-cell grid: 7.8 mm cells,
78+dt = 6.6 µs (a 152 kHz sample rate), and a string of 59 nodes.
79+
80+## Listening
81+
82+The microphone is a one-thread GPU dispatch riding in the same submission as
83+each timestep, so recording costs nothing and no readback happens until you
84+ask to listen. Because everything is in SI units, playback needs no
85+translation: one sample per timestep at 1/dt per second is real time at real
86+pitch. A .wav download (48 kHz, resampled) is one button over.
87+
88+Two ways to run. *Watching* takes a few timesteps per display frame, the
89+wave crawling in slow motion. *Render note* runs the solver flat out with no
90+display until the requested seconds of audio exist, then plays them. On a
91+discrete GPU a second of audio takes very roughly half a minute to render at
92+the default grid; on an integrated one it can take several minutes, which is
93+what the draft 64 grid is for (honest only to about 2.7 kHz, but several
94+times faster).
95+
96+## Honest limitations
97+
98+The grid resolves sound to about 5.5 kHz at the default size, so the top of
99+the timbre is simply absent, and the note is duller than a real instrument.
100+The soundboard has no modes of its own: the bridge drives the air directly,
101+so the body colours the sound only through its cavity and hole. The wall
102+mask is perfectly rigid apart from its absorption skin, and the box walls
103+are drawn at grid resolution, so plates thinner than about two cells leak.
104+The coupling gains are physically arbitrary; the equations are linear, so
105+they set relative balance, not absolute loudness.
106+
107+## Running it
108+
109+```
110+npm install
111+npm run dev # local dev server
112+npm test # solver checks against desktop WebGPU (Google Dawn)
113+npm run smoke # headless-browser check of the built page (environment permitting)
114+```
115+
116+`npm test` compiles the actual .m files to actual shaders and checks physics:
117+the measured pitch of the plucked string, the 60 dB decay time, that a sealed
118+box keeps sound out and a sound hole lets the cavity speak, and that the
119+microphone's spectrum sits on the string's partial comb.
index.htmladded+322−0View file
@@ -0,0 +1,322 @@
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/><rect x=%2210%22 y=%2258%22 width=%2280%22 height=%2224%22 rx=%224%22 fill=%22%23a07040%22/><circle cx=%2250%22 cy=%2270%22 r=%227%22 fill=%22%23f2f2f2%22/><path d=%22M10 48 Q 50 28 90 48%22 stroke=%22%23b40426%22 stroke-width=%224%22 fill=%22none%22/><path d=%22M10 48 L90 48%22 stroke=%22%233b4cc0%22 stroke-width=%222%22/></svg>" />
7+ <title>dulcimer — a plucked string and its box, in wave equations</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:disabled { cursor: default; opacity: 0.5; }
63+ button.primary { border-color: var(--accent); color: var(--accent); font-weight: 600; min-width: 5.5em; }
64+ .sliders {
65+ display: grid; grid-template-columns: repeat(auto-fill, minmax(240px, 1fr));
66+ gap: 2px 16px; padding: 4px 0;
67+ }
68+ #micbar { gap: 8px 16px; }
69+ #micparams { display: flex; flex-wrap: wrap; gap: 0 16px; padding: 0; }
70+ #recinfo { margin-top: 0; }
71+ .slider { display: flex; align-items: center; gap: 8px; font-size: 13px; color: var(--ink-2); }
72+ .slider > span:first-child { flex: 0 0 9.5em; text-align: right; }
73+ .slider > output { flex: 0 0 4.5em; font-variant-numeric: tabular-nums; color: var(--ink); }
74+ .group-title {
75+ font-size: 12px; text-transform: uppercase; letter-spacing: 0.05em;
76+ color: var(--ink-2); margin: 10px 0 0;
77+ }
78+ progress { width: 10em; height: 0.9em; }
79+ #stringbox {
80+ margin-top: 12px; border: 1px solid var(--line); border-radius: 8px;
81+ background: var(--panel-bg); position: relative;
82+ }
83+ #stringplot { width: 100%; height: 92px; display: block; }
84+ #stringlabel {
85+ position: absolute; top: 4px; left: 8px; font-size: 11px;
86+ color: var(--ink-2); font-variant-numeric: tabular-nums;
87+ }
88+ #stage {
89+ display: flex; gap: 14px; margin-top: 12px; align-items: stretch;
90+ border: 1px solid var(--line); border-radius: 8px; overflow: hidden;
91+ }
92+ .canvas-box {
93+ flex: 1; aspect-ratio: 16 / 10; max-height: 70vh; position: relative;
94+ background: #0b0e12;
95+ }
96+ #view { width: 100%; height: 100%; display: block; cursor: grab; touch-action: none; }
97+ #overlay { position: absolute; inset: 0; width: 100%; height: 100%; pointer-events: none; }
98+ .colorbar {
99+ display: flex; flex-direction: column; align-items: center; justify-content: center;
100+ gap: 4px; padding: 8px 4px; background: var(--panel-bg);
101+ width: 60px; flex: none; box-sizing: border-box;
102+ }
103+ .colorbar canvas { border: 1px solid var(--line); border-radius: 2px; }
104+ .colorbar-label { font-size: 11px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
105+ .readout { font-variant-numeric: tabular-nums; color: var(--ink); }
106+ .warn { color: #b35900; font-weight: 600; }
107+ .stats { margin-top: 8px; font-size: 13px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
108+ .stats b { color: var(--ink); font-weight: 600; }
109+ #err {
110+ color: #b35900; white-space: pre-wrap; font-size: 13px;
111+ font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
112+ }
113+ .editor {
114+ margin-top: 12px; border: 1px solid var(--line); border-radius: 8px;
115+ overflow: hidden;
116+ }
117+ .editor-head {
118+ display: flex; gap: 10px; align-items: center; justify-content: space-between;
119+ padding: 6px 10px; font-size: 12px; color: var(--ink-2);
120+ background: var(--panel-bg); border-bottom: 1px solid var(--line);
121+ }
122+ .editor-head button { padding: 2px 10px; font-size: 12px; }
123+ /* The fixed height lives on the row, not on either child: sizing the row
124+ makes both children stretch to one shared pixel height regardless of
125+ either's font size. */
126+ .editor-body { display: flex; align-items: stretch; height: 32em; }
127+ .editor-code { position: relative; flex: 1 1 62%; min-width: 0; }
128+ /* The overlay and the textarea must agree on every metric that affects
129+ where a character lands. Keep these two rules together. */
130+ .editor-code > pre,
131+ .editor-code > textarea {
132+ margin: 0; padding: 10px 12px; border: 0;
133+ box-sizing: border-box; width: 100%; height: 100%;
134+ font: 12.5px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
135+ tab-size: 2;
136+ white-space: pre; overflow-wrap: normal;
137+ }
138+ #highlight {
139+ position: absolute; inset: 0; overflow: hidden;
140+ pointer-events: none; background: var(--bg); color: var(--ink);
141+ }
142+ #source {
143+ position: relative; z-index: 1; display: block;
144+ resize: none; overflow: auto;
145+ background: transparent; color: transparent; caret-color: var(--ink);
146+ }
147+ #source:focus { outline: none; }
148+ /* Transparent text means the selection must be see-through, or selected
149+ code would be invisible. */
150+ #source::selection { background: color-mix(in srgb, var(--accent) 28%, transparent); }
151+ #compiled {
152+ flex: 1 1 38%; min-width: 0; margin: 0; padding: 10px 12px; overflow: auto;
153+ border-left: 1px solid var(--line); background: var(--panel-bg);
154+ font: 12px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
155+ color: var(--ink-2); white-space: pre;
156+ }
157+ @media (max-width: 860px) {
158+ .editor-body { flex-direction: column; height: auto; }
159+ .editor-code { flex: none; height: 26em; }
160+ #compiled { border-left: 0; border-top: 1px solid var(--line); max-height: 12em; }
161+ }
162+ .tok-com { color: var(--tok-com); }
163+ .tok-str { color: var(--tok-str); }
164+ .tok-num { color: var(--tok-num); }
165+ .tok-kw { color: var(--tok-kw); font-weight: 600; }
166+ .tok-ext { color: var(--tok-ext); }
167+ </style>
168+ </head>
169+ <body>
170+ <main>
171+ <h1>dulcimer</h1>
172+ <p class="sub">
173+ One plucked string and the box it sounds over, simulated as two coupled
174+ wave equations — a stiff damped string, and 3D acoustics around a rigid
175+ body — both written as the MATLAB below, run in your browser by
176+ <a href="https://numbl.org">numbl</a> and compiled to WebGPU compute
177+ shaders. Pluck it, watch the sound, put a microphone in the air, and
178+ listen.
179+ </p>
180+ <p class="sub" id="domaininfo"></p>
181+
182+ <div class="controls">
183+ <label title="Air grid points along the string; across and up get half each. Changing it recompiles.">grid
184+ <select id="gridsize">
185+ <option value="64">64 (draft)</option>
186+ <option value="96">96</option>
187+ <option value="128" selected>128</option>
188+ <option value="160">160</option>
189+ <option value="192">192</option>
190+ </select>
191+ </label>
192+ <label title="String length. Changing it recompiles.">string length
193+ <select id="stringlen">
194+ <option value="0.45">0.45 m</option>
195+ <option value="0.6" selected>0.6 m</option>
196+ <option value="0.72">0.72 m</option>
197+ </select>
198+ </label>
199+ <button id="runpause" class="primary">Run</button>
200+ <button id="pluck" title="Draw the string back into its pluck and release it, clearing the recording">Pluck</button>
201+ <label title="Timesteps taken between display frames while watching live. Below 1×, one step is taken every few frames.">speed
202+ <select id="spf">
203+ <option value="0.25">¼×</option>
204+ <option value="0.5">½×</option>
205+ <option value="1">1×</option>
206+ <option value="4">4×</option>
207+ <option value="16" selected>16×</option>
208+ <option value="64">64×</option>
209+ <option value="256">256×</option>
210+ </select>
211+ </label>
212+ <span style="flex:1"></span>
213+ <label>render
214+ <select id="renderdur">
215+ <option value="1">1 s</option>
216+ <option value="2" selected>2 s</option>
217+ <option value="3">3 s</option>
218+ <option value="5">5 s</option>
219+ </select>
220+ </label>
221+ <button id="rendernote" title="Pluck, then run the solver flat out with no display until this much of the note is recorded, and play it">Render note</button>
222+ <progress id="renderprog" value="0" max="1" hidden></progress>
223+ </div>
224+
225+ <p class="group-title">string</p>
226+ <div class="sliders" id="params"></div>
227+ <p class="group-title" id="scene-title">body</p>
228+ <div class="sliders" id="sceneparams"></div>
229+
230+ <p class="group-title">microphone</p>
231+ <div class="controls" id="micbar">
232+ <div class="sliders" id="micparams"></div>
233+ <button id="listen" title="Play what the microphone has recorded since the pluck, in real time at its real pitch">Listen</button>
234+ <button id="download" title="Save the recording as a 48 kHz .wav">Save .wav</button>
235+ <span class="stats" id="recinfo"></span>
236+ </div>
237+
238+ <div class="controls">
239+ <label>colormap
240+ <select id="colormap"></select>
241+ </label>
242+ <label title="Pressure the colour scale saturates at. Auto follows the field.">scale
243+ <select id="scalemode">
244+ <option value="auto" selected>auto</option>
245+ <option value="fixed">hold</option>
246+ </select>
247+ </label>
248+ <label>opacity
249+ <input type="range" id="opacity" min="0.2" max="6" step="0.1" value="2" />
250+ </label>
251+ <label>contrast
252+ <input type="range" id="contrast" min="0.5" max="3" step="0.1" value="1.6" />
253+ </label>
254+ <label title="Ray-march samples per pixel ray">quality
255+ <select id="quality">
256+ <option value="96">96</option>
257+ <option value="192" selected>192</option>
258+ <option value="320">320</option>
259+ </select>
260+ </label>
261+ <label title="Cut the picture open along the string: everything on the far side of this y-plane is hidden">clip
262+ <input type="range" id="clipy" min="-0.5" max="0.5" step="0.01" value="0.5" />
263+ </label>
264+ <label title="How many times larger than life the string's displacement is drawn in the 3D view">string ×
265+ <input type="range" id="exaggerate" min="1" max="100" step="1" value="25" />
266+ </label>
267+ </div>
268+ <div class="controls">
269+ <label><input type="checkbox" id="showfield" checked /> sound</label>
270+ <label><input type="checkbox" id="showbody" checked /> body</label>
271+ <label><input type="checkbox" id="showstring" checked /> string</label>
272+ <label><input type="checkbox" id="showwire" checked /> outline</label>
273+ <label><input type="checkbox" id="showmic" checked /> mic</label>
274+ </div>
275+
276+ <div id="stringbox">
277+ <canvas id="stringplot"></canvas>
278+ <span id="stringlabel"></span>
279+ </div>
280+
281+ <div id="stage">
282+ <div class="canvas-box">
283+ <canvas id="view"></canvas>
284+ <canvas id="overlay"></canvas>
285+ </div>
286+ <div id="colorbar"></div>
287+ </div>
288+ <p class="stats" id="stats"></p>
289+ <p id="err"></p>
290+
291+ <div class="editor">
292+ <div class="editor-head">
293+ <span>
294+ <select id="editor-file" aria-label="file to edit">
295+ <option value="model">model (.m) — the two solvers</option>
296+ <option value="scene">scene (.m) — the body</option>
297+ </select>
298+ <span id="editor-title"></span>
299+ </span>
300+ <span>
301+ <button id="recompile" type="button">Run edits</button>
302+ <button id="revert" type="button">Revert</button>
303+ </span>
304+ </div>
305+ <div class="editor-body">
306+ <div class="editor-code">
307+ <pre id="highlight" aria-hidden="true"></pre>
308+ <textarea
309+ id="source"
310+ spellcheck="false"
311+ autocomplete="off"
312+ autocapitalize="off"
313+ aria-label="source (MATLAB)"
314+ ></textarea>
315+ </div>
316+ <pre id="compiled"></pre>
317+ </div>
318+ </div>
319+ </main>
320+ <script type="module" src="/src/main.ts"></script>
321+ </body>
322+</html>
models/dulcimer.madded+80−0View file
@@ -0,0 +1,80 @@
1+% A plucked dulcimer string, and the air it sounds in. Two wave equations,
2+% advanced together, both compiled to WebGPU kernels.
3+%
4+% The string carries transverse displacement u (metres) on ns nodes along x:
5+%
6+% u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
7+%
8+% cs = 2 Ls f0 is the wave speed a string of length Ls needs to sound the
9+% fundamental f0. The u_xxxx term is bending stiffness, which sharpens the
10+% upper partials (inharmonicity B: partial n lands near n*f0*sqrt(1 + B n^2)).
11+% sig0 is plain decay — 6.91/t60 makes the amplitude fall 60 dB in t60
12+% seconds — and sig1 damps high frequencies faster than low, which is why a
13+% plucked note starts bright and mellows as it rings. The pluck is an initial
14+% condition: init below shapes the string into a triangle and releases it
15+% from rest. Both ends are pinned by the `pin` mask.
16+%
17+% The air carries pressure p on the box grid:
18+%
19+% p_tt + 2 sig p_t = c^2 lapw(p, wall) + s
20+%
21+% lapw is the wall-masked Laplacian: the dulcimer body's shell enters as the
22+% mask `wall` (0 in the shell, 1 in air), which makes the shell rigid — the
23+% Neumann condition, not a fast material, so it costs no timestep. The
24+% source s is the string, coupled two ways with an independent gain on each:
25+%
26+% gline the string radiates directly: its acceleration, spread along its
27+% line of cells (a thin string is a poor radiator in reality, so
28+% this is the idealized version of that);
29+% gbridge the string's pull on the bridge drives a patch of air above the
30+% top plate — the instrument's actual mechanism, minus the plate's
31+% own resonances, which this model does not carry.
32+%
33+% Either gain at zero switches that route off entirely. The absolute source
34+% strength is arbitrary (the equation is linear, and playback is normalized);
35+% the bridge branch carries a cs^2/Ls factor so that equal gains are of
36+% comparable loudness rather than one branch drowning the other.
37+
38+function [u, um, p, pm] = init(xs, npts, Ls, pluckpos, amp)
39+ % A triangle peaked at the pluck point, `amp` metres high, released from
40+ % rest (um = u). Written via abs because min(a, b) = (a + b - |a - b|)/2.
41+ a = xs / (pluckpos * Ls);
42+ b = (Ls - xs) / (Ls - pluckpos * Ls);
43+ u = amp * 0.5 * (a + b - abs(a - b));
44+ um = u;
45+ p = zeros(npts, 1);
46+ pm = zeros(npts, 1);
47+end
48+
49+function [un, uold, pn, pold] = step(u, um, p, pm, c, sig, wall, lineprof, boardprof, pin, dt, Ls, f0, B, t60, sig1, gline, gbridge)
50+ % -- the string ---------------------------------------------------------
51+ cs = 2 * Ls * f0; % tension, expressed as a wave speed
52+ kap2 = B * (cs * Ls / pi)^2; % stiffness, expressed as inharmonicity
53+ sig0 = 6.91 / t60;
54+
55+ lu = dxx(u);
56+ lum = dxx(um);
57+ l4 = dxxxx(u);
58+ d = sig0 * dt;
59+ un = pin .* ((2*u - (1 - d)*um + (dt*dt)*((cs*cs)*lu - kap2*l4) + (2*sig1*dt)*(lu - lum)) ./ (1 + d));
60+
61+ % -- the coupling -------------------------------------------------------
62+ % The string's acceleration, sampled at each air cell's own x; and the
63+ % slope at the bridge (the pull of the tension on it), broadcast. The
64+ % scene says where each acts, through lineprof and boardprof. cs^2/Ls
65+ % puts the bridge branch on the same footing as the acceleration branch.
66+ acc = (un - 2*u + um) / (dt*dt);
67+ sl = spread(acc);
68+ fb = bridge(un);
69+ s = gline * (lineprof .* sl) - (gbridge * (cs*cs) / Ls) * (boardprof .* fb);
70+
71+ % -- the air ------------------------------------------------------------
72+ lp = lapw(p, wall);
73+ sd = sig * dt;
74+ pn = wall .* ((2*p - (1 - sd) .* pm + (c*dt).^2 .* lp + (dt*dt) * s) ./ (1 + sd));
75+
76+ % This step's fields become the next step's history. Lines of their own,
77+ % so each plans as the copy it is.
78+ uold = u;
79+ pold = p;
80+end
package-lock.jsonadded+3475−0View file
This diff is 3,480 lines long and is not shown.
package.jsonadded+31−0View file
@@ -0,0 +1,31 @@
1+{
2+ "name": "dulcimer",
3+ "version": "0.1.0",
4+ "description": "A plucked dulcimer string and its resonating body, simulated by wave equations in the browser: MATLAB solvers 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"
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/box.madded+72−0View file
@@ -0,0 +1,72 @@
1+% The dulcimer's body: a rigid box under the string, with a round sound hole
2+% in its top plate.
3+%
4+% The top plate lies at z = 0. The box hangs below it (down to z = -boxd) and
5+% the string runs along x just above it, at height `gap`. Everything here is
6+% metres and inverse seconds.
7+%
8+% The shell is a *mask*, not a material: `wall` is 0 inside the solid and 1
9+% in air, and the solver's masked Laplacian drops any flux into a masked
10+% cell, which is the rigid (Neumann) condition. So the walls reflect
11+% everything and cost nothing — no timestep penalty, unlike a fast material.
12+% What a rigid mask cannot do is absorb, so `absorb` paints a thin layer of
13+% ordinary volume absorption over the shell's surfaces (the smoothed mask's
14+% transition skin, 4*wall*(1-wall), which peaks exactly at the surface); that
15+% is what lets the cavity's ring be tamed like real wood tames it.
16+%
17+% The masks are drawn with a sharp profile (about half a cell) on purpose. A
18+% box is grid-aligned, so sharpness costs no staircase artefacts, and a plate
19+% only two cells thick has to reach mask ≈ 0 at its centre or it leaks.
20+function [c, sig, wall, lineprof, boardprof] = medium(x, y, z, h, Lx, Ly, Lz, c0, Ls, body, boxl, boxw, boxd, thick, holer, holex, gap, patchr, absorb)
21+ s = 0.45 * h;
22+
23+ % The solid: the outer box minus its interior cavity.
24+ ox = 0.5 * (1 - tanh((abs(x) - boxl/2) / s));
25+ oy = 0.5 * (1 - tanh((abs(y) - boxw/2) / s));
26+ oz = 0.5 * (1 - tanh((abs(z + boxd/2) - boxd/2) / s));
27+ ix = 0.5 * (1 - tanh((abs(x) - (boxl/2 - thick)) / s));
28+ iy = 0.5 * (1 - tanh((abs(y) - (boxw/2 - thick)) / s));
29+ iz = 0.5 * (1 - tanh((abs(z + boxd/2) - (boxd/2 - thick)) / s));
30+ shell = max(0, ox .* oy .* oz - ix .* iy .* iz);
31+
32+ % The sound hole: a cylinder of radius holer at (holex, 0), cut out of the
33+ % top plate only.
34+ r = sqrt((x - holex).^2 + y.^2);
35+ hole = 0.5 * (1 - tanh((r - holer) / s));
36+ top = 0.5 * (1 - tanh((abs(z + thick/2) - thick/2) / s));
37+ shell = max(0, shell - hole .* top);
38+
39+ % `body` at 0 takes the box away entirely: a bare string in open air. The
40+ % bridge patch stays where the plate would have been, so the bridge route
41+ % keeps working as an abstract source even with nothing to push against.
42+ shell = body * shell;
43+
44+ wall = 1 - shell;
45+
46+ % Uniform air. The walls being a mask is what keeps this uniform — and the
47+ % timestep as large as air allows.
48+ c = c0;
49+
50+ % Open boundary at the domain edge, plus the wall's surface absorption.
51+ sig = sponge3(x, y, z, Lx, Ly, Lz, 0.15*Ly, 9000) + absorb * (4 * shell .* (1 - shell));
52+
53+ % Where the string radiates directly: a Gaussian tube two cells wide
54+ % around the string's line (y = 0, z = gap), fading out at its ends.
55+ % Normalized so its cross-section integrates to one — the tube is a stand-in
56+ % for a line, and without the 1/(pi a^2) its strength would depend on the
57+ % grid that draws it.
58+ zs = gap;
59+ a = 2*h;
60+ tube = exp(-(y.^2 + (z - zs).^2) / a^2) / (pi * a^2);
61+ span = 0.5 * (1 + tanh((Ls/2 - abs(x)) / (1.5*h)));
62+ lineprof = tube .* span .* wall;
63+
64+ % Where the bridge force drives the air: a patch of the air just above the
65+ % top plate, around the bridge end of the string (x = +Ls/2). This stands
66+ % in for the top plate moving; the plate's own modes are not modelled.
67+ % Normalized to integrate to one over its volume, like the tube.
68+ w = 1.5*h;
69+ pr = exp(-((x - Ls/2).^2 + y.^2) / patchr^2) / (pi * patchr^2);
70+ skin = exp(-(z - w).^2 / w^2) / (sqrt(pi) * w);
71+ boardprof = pr .* skin .* wall;
72+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+205−0View file
@@ -0,0 +1,205 @@
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 compile to finish, starts
6+ * the run, and checks that steps are being taken and the microphone is
7+ * filling. Any console error, page error or failed 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 string is in the right place.
11+ * Run it after `vite build`: node scripts/smoke.mjs
12+ *
13+ * `--full` additionally exercises a broken edit, the revert, and playback.
14+ * Those need new GPU work started while the page is already drawing, which
15+ * headless Chrome cannot always do (see the 2d sibling's smoke script for
16+ * the history); a failure from them says as much about the browser as about
17+ * the app. What they would have covered on the solver side is covered by
18+ * `npm run test:node`, which runs against desktop WebGPU.
19+ */
20+const withGpuWork = process.argv.includes('--full');
21+import { createServer } from 'node:http';
22+import { readFile } from 'node:fs/promises';
23+import { extname, join } from 'node:path';
24+import puppeteer from 'puppeteer-core';
25+
26+const DIST = new URL('../dist/', import.meta.url).pathname;
27+const CHROME = process.env.CHROME_PATH ?? '/usr/bin/google-chrome';
28+const MIME = {
29+ '.html': 'text/html',
30+ '.js': 'text/javascript',
31+ '.css': 'text/css',
32+ '.json': 'application/json',
33+};
34+
35+const server = createServer(async (req, res) => {
36+ try {
37+ const path = req.url === '/' ? '/index.html' : req.url.split('?')[0];
38+ const data = await readFile(join(DIST, path));
39+ res.writeHead(200, { 'content-type': MIME[extname(path)] ?? 'application/octet-stream' });
40+ res.end(data);
41+ } catch {
42+ res.writeHead(404);
43+ res.end('not found');
44+ }
45+});
46+await new Promise((r) => server.listen(0, '127.0.0.1', r));
47+const port = server.address().port;
48+
49+const flagSets = [
50+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--enable-features=Vulkan'],
51+ [
52+ '--headless=new',
53+ '--no-sandbox',
54+ '--enable-unsafe-webgpu',
55+ '--use-webgpu-adapter=swiftshader',
56+ '--enable-unsafe-swiftshader',
57+ ],
58+];
59+
60+let ok = false;
61+let lastFailure = 'never ran';
62+for (const flags of flagSets) {
63+ const browser = await puppeteer.launch({
64+ executablePath: CHROME,
65+ args: [...flags],
66+ protocolTimeout: 600_000,
67+ });
68+ const problems = [];
69+ try {
70+ const page = await browser.newPage();
71+ page.on('console', (m) => {
72+ if (m.type() === 'error') problems.push(`console: ${m.text()}`);
73+ });
74+ page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`));
75+ page.on('requestfailed', (r) => problems.push(`request failed: ${r.url()}`));
76+
77+ await page.goto(`http://127.0.0.1:${port}/index.html`, { waitUntil: 'load' });
78+
79+ // The page loads paused, at the draft grid to keep SwiftShader honest.
80+ // Wait for the first compile to have produced an op list.
81+ await page.waitForFunction(
82+ () => (document.getElementById('compiled')?.textContent ?? '').includes('external'),
83+ { timeout: 180_000, polling: 500 },
84+ );
85+ const errText = () => page.$eval('#err', (n) => n.textContent ?? '');
86+ if (await errText()) problems.push(`compile reported: ${await errText()}`);
87+
88+ // The compiled plan must contain every coupling op.
89+ const compiled = await page.$eval('#compiled', (n) => n.textContent ?? '');
90+ for (const op of ['lapw(p, wall)', 'spread(acc)', 'bridge(un)', 'dxxxx(u)']) {
91+ if (!compiled.includes(op)) problems.push(`compiled plan is missing ${op}`);
92+ }
93+
94+ // Drop to the draft grid so a software rasterizer can take steps at all.
95+ await page.select('#gridsize', '64');
96+ await page.waitForFunction(
97+ () => /64×32×32/.test(document.getElementById('domaininfo')?.textContent ?? ''),
98+ { timeout: 180_000, polling: 250 },
99+ );
100+
101+ // A paused page reports no frame rate.
102+ if (/ms\/frame/.test(await page.$eval('#stats', (n) => n.textContent ?? ''))) {
103+ problems.push('a paused page reported a frame rate');
104+ }
105+ await page.click('#runpause');
106+
107+ /** The stats line reports the step count, so waiting on it says both that
108+ * the solver ran and that the frame loop is turning. */
109+ const running = async (label) => {
110+ await page.waitForFunction(
111+ () => {
112+ const m = /step ([\d,]+)/.exec(document.getElementById('stats')?.textContent ?? '');
113+ return m ? Number(m[1].replace(/,/g, '')) > 50 : false;
114+ },
115+ { timeout: 180_000, polling: 500 },
116+ );
117+ console.log(` ${label}: ${(await page.$eval('#stats', (n) => n.textContent)).trim()}`);
118+ };
119+ await running('start');
120+
121+ // The microphone must be filling, one sample per timestep.
122+ await page.waitForFunction(
123+ () => /recorded/.test(document.getElementById('recinfo')?.textContent ?? ''),
124+ { timeout: 60_000, polling: 250 },
125+ );
126+ console.log(` microphone: ${(await page.$eval('#recinfo', (n) => n.textContent ?? '')).trim()}`);
127+
128+ // A body edit re-evaluates the scene .m without recompiling: shrinking
129+ // the sound hole to nothing must not error and must keep stepping.
130+ await page.evaluate(() => {
131+ for (const row of document.querySelectorAll('#sceneparams label.slider')) {
132+ if (!row.querySelector('span')?.textContent?.startsWith('sound hole radius')) continue;
133+ const input = row.querySelector('input');
134+ input.value = '0';
135+ input.dispatchEvent(new Event('input'));
136+ }
137+ });
138+ await new Promise((r) => setTimeout(r, 1500));
139+ if (await errText()) problems.push(`scene edit reported: ${await errText()}`);
140+ await running('after body edit');
141+
142+ if (!withGpuWork) {
143+ console.log(' (skipping the broken-edit and playback checks; pass --full to run them)');
144+ } else {
145+ // A broken edit must be reported rather than thrown, and must not take
146+ // the page down with it.
147+ await page.evaluate(() => {
148+ const ta = document.getElementById('source');
149+ ta.value = ta.value.replace('lapw(p, wall)', 'lapnope(p, wall)');
150+ ta.dispatchEvent(new Event('input'));
151+ });
152+ await page.click('#recompile');
153+ await page.waitForFunction(
154+ () => (document.getElementById('err')?.textContent ?? '').length > 0,
155+ { timeout: 120_000, polling: 250 },
156+ );
157+ console.log(` bad edit reported: ${(await errText()).split('\n')[0]}`);
158+
159+ // And reverting must put it back.
160+ await page.click('#revert');
161+ await page.waitForFunction(
162+ () => (document.getElementById('err')?.textContent ?? '').length === 0,
163+ { timeout: 120_000, polling: 250 },
164+ );
165+ await running('after revert');
166+
167+ await page.click('#listen');
168+ await new Promise((r) => setTimeout(r, 2000));
169+ const listenErr = await errText();
170+ if (listenErr) problems.push(`Listen reported: ${listenErr}`);
171+ }
172+
173+ const report = await page.evaluate(() => ({
174+ stats: document.getElementById('stats')?.textContent ?? '',
175+ err: document.getElementById('err')?.textContent ?? '',
176+ painted: (() => {
177+ const canvas = document.getElementById('view');
178+ return canvas instanceof HTMLCanvasElement && canvas.width > 0;
179+ })(),
180+ }));
181+
182+ console.log(`flags: ${flags.join(' ')}`);
183+ console.log(report.stats.trim());
184+ if (report.err) problems.push(`page error box: ${report.err}`);
185+ if (!report.painted) problems.push('canvas was never sized');
186+ if (problems.length === 0) {
187+ ok = true;
188+ } else {
189+ lastFailure = problems.join('\n');
190+ }
191+ } catch (e) {
192+ lastFailure = [`${e}`, ...problems].join('\n');
193+ } finally {
194+ await browser.close();
195+ }
196+ if (ok) break;
197+}
198+
199+server.close();
200+if (!ok) {
201+ console.error(`smoke: FAILED\n${lastFailure}`);
202+ process.exit(1);
203+}
204+console.log('smoke: the page runs');
205+process.exit(0);
scripts/test-node.tsadded+63−0View file
@@ -0,0 +1,63 @@
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+ planChecks,
14+ pitchChecks,
15+ decayChecks,
16+ wallChecks,
17+ holeChecks,
18+ spectrumChecks,
19+ bareChecks,
20+ splitChecks,
21+} from '../test/checks.ts';
22+
23+let failures = 0;
24+const check = (name: string, ok: boolean, detail: string): void => {
25+ console.log(`${ok ? 'PASS' : 'FAIL'} ${name} — ${detail}`);
26+ if (!ok) failures++;
27+};
28+const log = (s: string): void => console.log(s);
29+
30+/**
31+ * These checks compile MATLAB to compute shaders, so they need a GPU.
32+ * `--skip-without-gpu` lets a runner that has none say so and move on; a plain
33+ * local run still fails loudly, so a missing GPU is never mistaken for a pass.
34+ */
35+const skipWithoutGpu = process.argv.includes('--skip-without-gpu');
36+
37+let runtime: string;
38+let device: GPUDevice;
39+try {
40+ runtime = await installWebGpu();
41+ device = await requestAcousticDevice();
42+} catch (e) {
43+ const detail = `${errMsg(e)}\n${NO_ADAPTER_HINT}`;
44+ if (skipWithoutGpu) {
45+ console.log(`SKIP no WebGPU available here, so these checks did not run.\n${detail}`);
46+ process.exit(0);
47+ }
48+ console.error(`test-node: ${detail}`);
49+ process.exit(1);
50+}
51+console.log(`dulcimer tests — ${runtime}\n`);
52+
53+await planChecks(device, check, log);
54+await pitchChecks(device, check, log);
55+await decayChecks(device, check, log);
56+await wallChecks(device, check, log);
57+await holeChecks(device, check, log);
58+await spectrumChecks(device, check, log);
59+await bareChecks(device, check, log);
60+await splitChecks(device, check, log);
61+
62+console.log(failures ? `\n${failures} failure(s)` : '\nall checks passed');
63+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+231−0View file
@@ -0,0 +1,231 @@
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 7 seconds of audio at the default
18+ * timestep, and 4 MB of GPU memory. Recording stops when it is full rather
19+ * than wrapping, so what you hear always starts where the pluck did. */
20+export const TRACE_CAPACITY = 1 << 20;
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+ nz: number;
52+}
53+
54+export class Recorder {
55+ readonly capacity = TRACE_CAPACITY;
56+
57+ #device: GPUDevice;
58+ #pipeline: GPUComputePipeline | null = null;
59+ #layout: GPUBindGroupLayout;
60+ #bindGroup: GPUBindGroup | null = null;
61+ #trace: GPUBuffer;
62+ #head: GPUBuffer;
63+ #probe: GPUBuffer;
64+ #readback: GPUBuffer;
65+ #field: GPUBuffer;
66+ #nx: number;
67+ #ny: number;
68+ #nz: number;
69+ /** Samples written since the last clear, as far as the host knows. Counted
70+ * here rather than read back from the GPU: the dispatch runs once per step
71+ * and the host knows exactly how many steps it asked for. */
72+ #count = 0;
73+ #reading = false;
74+
75+ private constructor(init: {
76+ device: GPUDevice;
77+ layout: GPUBindGroupLayout;
78+ trace: GPUBuffer;
79+ head: GPUBuffer;
80+ probe: GPUBuffer;
81+ readback: GPUBuffer;
82+ field: GPUBuffer;
83+ nx: number;
84+ ny: number;
85+ nz: number;
86+ }) {
87+ this.#device = init.device;
88+ this.#layout = init.layout;
89+ this.#trace = init.trace;
90+ this.#head = init.head;
91+ this.#probe = init.probe;
92+ this.#readback = init.readback;
93+ this.#field = init.field;
94+ this.#nx = init.nx;
95+ this.#ny = init.ny;
96+ this.#nz = init.nz;
97+ }
98+
99+ static async create(opts: RecorderOptions): Promise<Recorder> {
100+ const { device } = opts;
101+ const layout = device.createBindGroupLayout({
102+ label: 'recorder',
103+ entries: [
104+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
105+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
106+ {
107+ binding: 2,
108+ visibility: GPUShaderStage.COMPUTE,
109+ buffer: { type: 'read-only-storage' },
110+ },
111+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform' } },
112+ ],
113+ });
114+ const trace = device.createBuffer({
115+ label: 'recorder-trace',
116+ size: 4 * TRACE_CAPACITY,
117+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
118+ });
119+ const head = device.createBuffer({
120+ label: 'recorder-head',
121+ size: 4,
122+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
123+ });
124+ const probe = device.createBuffer({
125+ label: 'recorder-probe',
126+ size: 8,
127+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
128+ });
129+ const readback = device.createBuffer({
130+ label: 'recorder-readback',
131+ size: 4 * TRACE_CAPACITY,
132+ usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
133+ });
134+
135+ const rec = new Recorder({
136+ device, layout, trace, head, probe, readback,
137+ field: opts.field, nx: opts.nx, ny: opts.ny, nz: opts.nz,
138+ });
139+ rec.#pipeline = await device.createComputePipelineAsync({
140+ label: 'recorder',
141+ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }),
142+ compute: {
143+ module: device.createShaderModule({ code: SHADER, label: 'recorder' }),
144+ entryPoint: 'main',
145+ },
146+ });
147+ rec.#bindGroup = device.createBindGroup({
148+ layout,
149+ entries: [
150+ { binding: 0, resource: { buffer: trace } },
151+ { binding: 1, resource: { buffer: head } },
152+ { binding: 2, resource: { buffer: opts.field } },
153+ { binding: 3, resource: { buffer: probe } },
154+ ],
155+ });
156+ rec.clear();
157+ return rec;
158+ }
159+
160+ /** Samples recorded so far. Stops rising once the trace is full. */
161+ get count(): number {
162+ return Math.min(this.#count, this.capacity);
163+ }
164+
165+ get full(): boolean {
166+ return this.#count >= this.capacity;
167+ }
168+
169+ /** Put the microphone at the grid point nearest (ix, iy, iz). Free: the
170+ * probe index is a uniform, so moving it disturbs neither the run nor the
171+ * recording already made. */
172+ setProbe(ix: number, iy: number, iz: number): void {
173+ const cx = Math.max(0, Math.min(this.#nx - 1, Math.round(ix)));
174+ const cy = Math.max(0, Math.min(this.#ny - 1, Math.round(iy)));
175+ const cz = Math.max(0, Math.min(this.#nz - 1, Math.round(iz)));
176+ this.#device.queue.writeBuffer(
177+ this.#probe,
178+ 0,
179+ new Uint32Array([cx + this.#nx * (cy + this.#ny * cz), this.capacity]),
180+ );
181+ }
182+
183+ /** Start again from an empty trace. */
184+ clear(): void {
185+ this.#count = 0;
186+ this.#device.queue.writeBuffer(this.#head, 0, new Uint32Array([0]));
187+ }
188+
189+ /** Record one sample. Called once per timestep, inside the step's own
190+ * submission, so no extra work crosses to the host. */
191+ encode(encoder: GPUCommandEncoder): void {
192+ if (!this.#pipeline || !this.#bindGroup || this.full) return;
193+ const pass = encoder.beginComputePass({ label: 'recorder' });
194+ pass.setPipeline(this.#pipeline);
195+ pass.setBindGroup(0, this.#bindGroup);
196+ pass.dispatchWorkgroups(1);
197+ pass.end();
198+ this.#count++;
199+ }
200+
201+ /** The recorded trace. The only readback the microphone ever does. */
202+ async read(): Promise<Float32Array> {
203+ const n = this.count;
204+ if (n === 0) return new Float32Array(0);
205+ if (this.#reading) throw new Error('a trace readback is already in flight');
206+ this.#reading = true;
207+ try {
208+ const enc = this.#device.createCommandEncoder({ label: 'recorder-read' });
209+ enc.copyBufferToBuffer(this.#trace, 0, this.#readback, 0, 4 * n);
210+ this.#device.queue.submit([enc.finish()]);
211+ await this.#readback.mapAsync(GPUMapMode.READ, 0, 4 * n);
212+ const out = new Float32Array(this.#readback.getMappedRange(0, 4 * n).slice(0));
213+ this.#readback.unmap();
214+ return out;
215+ } finally {
216+ this.#reading = false;
217+ }
218+ }
219+
220+ /** The pressure buffer this microphone listens to. */
221+ get field(): GPUBuffer {
222+ return this.#field;
223+ }
224+
225+ destroy(): void {
226+ this.#trace.destroy();
227+ this.#head.destroy();
228+ this.#probe.destroy();
229+ this.#readback.destroy();
230+ }
231+}
src/audio/wav.tsadded+59−0View file
@@ -0,0 +1,59 @@
1+/**
2+ * A recorded trace as a .wav file.
3+ *
4+ * The trace holds one sample per timestep, 1/dt of them a second — around
5+ * 150 kHz at the default grid, which is a legal but eccentric rate for a
6+ * .wav. It is resampled to 48 kHz by linear interpolation, which is harmless
7+ * here: the content is band-limited far below either rate (the grid stops
8+ * resolving sound around 5 kHz). Normalized to peak, like playback, and
9+ * faded a few milliseconds at each end so the file does not open and close
10+ * with a click.
11+ */
12+
13+export const WAV_RATE = 48000;
14+
15+export function traceToWav(trace: Float32Array, dt: number): Blob {
16+ const n = Math.max(1, Math.round(trace.length * dt * WAV_RATE));
17+ let peak = 0;
18+ for (const v of trace) peak = Math.max(peak, Math.abs(v));
19+ const gain = peak > 0 ? 0.9 / peak : 0;
20+
21+ const samples = new Float32Array(n);
22+ for (let i = 0; i < n; i++) {
23+ const s = i / (WAV_RATE * dt);
24+ const k = Math.min(trace.length - 2, Math.floor(s));
25+ const f = Math.min(1, s - k);
26+ samples[i] = gain * ((1 - f) * trace[k] + f * trace[k + 1]);
27+ }
28+ const fade = Math.min(Math.round(0.005 * WAV_RATE), Math.floor(n / 2));
29+ for (let i = 0; i < fade; i++) {
30+ const w = i / fade;
31+ samples[i] *= w;
32+ samples[n - 1 - i] *= w;
33+ }
34+
35+ const bytes = 44 + 2 * n;
36+ const buf = new ArrayBuffer(bytes);
37+ const view = new DataView(buf);
38+ const str = (off: number, s: string): void => {
39+ for (let i = 0; i < s.length; i++) view.setUint8(off + i, s.charCodeAt(i));
40+ };
41+ str(0, 'RIFF');
42+ view.setUint32(4, bytes - 8, true);
43+ str(8, 'WAVE');
44+ str(12, 'fmt ');
45+ view.setUint32(16, 16, true); // PCM chunk size
46+ view.setUint16(20, 1, true); // PCM
47+ view.setUint16(22, 1, true); // mono
48+ view.setUint32(24, WAV_RATE, true);
49+ view.setUint32(28, 2 * WAV_RATE, true); // byte rate
50+ view.setUint16(32, 2, true); // block align
51+ view.setUint16(34, 16, true); // bits per sample
52+ str(36, 'data');
53+ view.setUint32(40, 2 * n, true);
54+ for (let i = 0; i < n; i++) {
55+ const v = Math.max(-1, Math.min(1, samples[i]));
56+ view.setInt16(44 + 2 * i, Math.round(v * 32767), true);
57+ }
58+ return new Blob([buf], { type: 'audio/wav' });
59+}
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+ * field its line reads, plus its output and the parameter block, and the air
6+ * update reads a lot of fields at once — the two pressure histories, the
7+ * medium, the wall mask, the Laplacian, the coupling source. WebGPU only
8+ * guarantees 8 storage buffers per compute stage, so we ask for whatever the
9+ * adapter will give up to 16. When that is not enough the planner splits the
10+ * kernel instead of failing (see `fitToBudget` in plan.ts), so this is a
11+ * 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+145−0View file
@@ -0,0 +1,145 @@
1+/**
2+ * The two computational grids: a rectangular box of air, and a line of
3+ * string.
4+ *
5+ * The air grid is cell-centred and uniform, with the same spacing h along
6+ * every axis; the box is twice as long in x (the string direction) as in y
7+ * and z, so nx = 2*ny = 2*nz and a field is nx*ny*nz f32 values, flattened
8+ * x-fastest: the point (ix, iy, iz) is element `ix + nx*(iy + ny*iz)`.
9+ * Cell-centred so that no point sits exactly on the outer boundary: the
10+ * stencil takes the field outside the domain to be zero, and the absorbing
11+ * layer is meant to have swallowed the wave before it gets there.
12+ *
13+ * The string grid is node-centred: ns points from x = 0 to x = Ls inclusive,
14+ * spacing hs = Ls/(ns-1), so its endpoints land exactly on the terminations
15+ * that pin them. In the air's coordinates the string runs from -Ls/2 to
16+ * +Ls/2 along x, at y = 0 and a height above the top plate that the scene
17+ * decides.
18+ */
19+import { CFL } from './units.ts';
20+
21+export interface AirGrid {
22+ nx: number;
23+ ny: number;
24+ nz: number;
25+ npts: number;
26+ /** Domain extents in metres, centred on the origin. */
27+ Lx: number;
28+ Ly: number;
29+ Lz: number;
30+ /** Grid spacing, the same in all three directions. */
31+ h: number;
32+ /** Coordinates of every grid point, npts each, x fastest — as the shaders
33+ * see them (f32) and as the scene .m is evaluated at (f64). */
34+ x: Float32Array;
35+ y: Float32Array;
36+ z: Float32Array;
37+ x64: Float64Array;
38+ y64: Float64Array;
39+ z64: Float64Array;
40+}
41+
42+export function makeAirGrid(nx: number, Lx: number): AirGrid {
43+ const ny = nx / 2;
44+ const nz = nx / 2;
45+ if (!Number.isInteger(ny)) throw new Error(`grid nx = ${nx} must be even`);
46+ const h = Lx / nx;
47+ const Ly = ny * h;
48+ const Lz = nz * h;
49+ const npts = nx * ny * nz;
50+ const x64 = new Float64Array(npts);
51+ const y64 = new Float64Array(npts);
52+ const z64 = new Float64Array(npts);
53+ for (let iz = 0; iz < nz; iz++) {
54+ const zv = -Lz / 2 + (iz + 0.5) * h;
55+ for (let iy = 0; iy < ny; iy++) {
56+ const yv = -Ly / 2 + (iy + 0.5) * h;
57+ const base = nx * (iy + ny * iz);
58+ for (let ix = 0; ix < nx; ix++) {
59+ const k = base + ix;
60+ x64[k] = -Lx / 2 + (ix + 0.5) * h;
61+ y64[k] = yv;
62+ z64[k] = zv;
63+ }
64+ }
65+ }
66+ return {
67+ nx, ny, nz, npts, Lx, Ly, Lz, h,
68+ x: new Float32Array(x64),
69+ y: new Float32Array(y64),
70+ z: new Float32Array(z64),
71+ x64, y64, z64,
72+ };
73+}
74+
75+export interface StringGrid {
76+ ns: number;
77+ /** String length, metres. */
78+ Ls: number;
79+ /** Node spacing, Ls/(ns-1). */
80+ hs: number;
81+ /** Node positions, 0 to Ls inclusive. */
82+ xs: Float32Array;
83+ xs64: Float64Array;
84+ /** 1 at interior nodes, 0 at the two pinned ends. The compiled step
85+ * multiplies its update by this, which is what terminates the string. */
86+ pin: Float32Array;
87+}
88+
89+/**
90+ * The timestep the explicit leapfrog on the air grid is stable at, seconds.
91+ *
92+ * Leapfrog on p_tt = c^2 lap(p) is stable while dt^2 c^2 |lap|max <= 4, and
93+ * the 7-point Laplacian's extreme eigenvalue is 12/h^2, so the condition is
94+ * c*dt/h <= 1/sqrt(3). `cmax` is the fastest sound speed anywhere in the
95+ * medium; since the body's walls are a mask rather than a fast material (see
96+ * src/mgpu/ops.ts), cmax is normally just the speed of air.
97+ */
98+export function stableDt(h: number, cmax: number, cfl = CFL): number {
99+ return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
100+}
101+
102+/**
103+ * The coarsest string spacing the stiff-string scheme demands at timestep dt,
104+ * in metres.
105+ *
106+ * The scheme (models/dulcimer.m) is the standard explicit leapfrog for
107+ * u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
108+ * whose stability condition (Bilbao, Numerical Sound Synthesis, ch. 7) is
109+ * hs^2 >= (a + sqrt(a^2 + 16 kap^2 dt^2)) / 2, a = cs^2 dt^2 + 4 sig1 dt.
110+ */
111+export function stableHs(dt: number, cs: number, kap2: number, sig1: number): number {
112+ const a = cs * cs * dt * dt + 4 * sig1 * dt;
113+ return Math.sqrt((a + Math.sqrt(a * a + 16 * kap2 * dt * dt)) / 2);
114+}
115+
116+/**
117+ * Build the string grid for length Ls at timestep dt, sized so that the
118+ * scheme is stable for every value the parameter sliders can reach.
119+ *
120+ * The air's CFL condition fixes dt, so the string has no say in the timestep;
121+ * what it gets to choose is its own spacing, and the finest stable spacing is
122+ * what maximizes the string's bandwidth. Sizing for the sliders' worst case
123+ * (highest fundamental, most stiffness, most frequency-dependent damping)
124+ * rather than their current values means moving a slider never forces a
125+ * recompile: the grid stays valid, merely a little coarser than that setting
126+ * alone would need.
127+ */
128+export function makeStringGrid(
129+ Ls: number,
130+ dt: number,
131+ worst: { f0: number; B: number; sig1: number },
132+): StringGrid {
133+ const cs = 2 * Ls * worst.f0;
134+ const kap2 = (worst.B * cs * cs * Ls * Ls) / (Math.PI * Math.PI);
135+ // 0.95: a little margin under the exact limit, as CFL is for the air.
136+ const hsMin = stableHs(dt, cs, kap2, worst.sig1) / 0.95;
137+ const ns = Math.max(8, Math.min(256, Math.floor(Ls / hsMin) + 1));
138+ const hs = Ls / (ns - 1);
139+ const xs64 = new Float64Array(ns);
140+ for (let i = 0; i < ns; i++) xs64[i] = i * hs;
141+ const pin = new Float32Array(ns).fill(1);
142+ pin[0] = 0;
143+ pin[ns - 1] = 0;
144+ return { ns, Ls, hs, xs: new Float32Array(xs64), xs64, pin };
145+}
src/main.tsadded+693−0View file
@@ -0,0 +1,693 @@
1+/**
2+ * The app: two MATLAB files, a GPU, a string, a box, and a microphone.
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
6+ * the *body* re-evaluates the scene .m on the CPU and re-uploads five arrays,
7+ * which is cheap and needs no recompile. Changing the grid, the string
8+ * length, or either file's text recompiles — a fresh session, from source to
9+ * shaders.
10+ *
11+ * There are two ways to run. *Watching*: a few timesteps per frame, the wave
12+ * crawling in slow motion. *Rendering a note*: the solver flat out with no
13+ * display until a chosen duration of audio exists, then playback. Both feed
14+ * the same GPU-side microphone.
15+ */
16+import { requestAcousticDevice } from './device.ts';
17+import { ModelSession } from './mgpu/session.ts';
18+import { EXTERNAL_OPS } from './mgpu/externals.ts';
19+import { formatFailure, ModelCompileError } from './mgpu/errors.ts';
20+import {
21+ dulcimerModel,
22+ defaultParams,
23+ type Params,
24+ type ParamSpec,
25+} from './mgpu/registry.ts';
26+import { boxScene, defaultSceneParams, type MScene } from './scene/registry.ts';
27+import { VolumeView, cameraFrame, type Camera } from './render/volume.ts';
28+import { Overlay } from './render/overlay.ts';
29+import { StringPlot } from './render/stringplot.ts';
30+import { Colorbar, fmtValue } from './render/colorbar.ts';
31+import { colormaps, colormapNames } from './render/colormaps.ts';
32+import { CodeEditor } from './editor/codeEditor.ts';
33+import { planPlayback, playTrace, stop as stopAudio } from './audio/play.ts';
34+import { traceToWav } from './audio/wav.ts';
35+import { C_AIR, POOR_RESOLUTION, fmtLength, fmtTime } from './units.ts';
36+
37+const el = <T extends HTMLElement>(id: string): T => {
38+ const node = document.getElementById(id);
39+ if (!node) throw new Error(`missing element #${id}`);
40+ return node as T;
41+};
42+
43+const errBox = el<HTMLParagraphElement>('err');
44+const showError = (e: unknown, source: string): void => {
45+ errBox.textContent = formatFailure(e, source);
46+};
47+const clearError = (): void => {
48+ errBox.textContent = '';
49+};
50+
51+/* ---------------------------------------------------------------- state -- */
52+
53+const model = dulcimerModel;
54+const scene: MScene = boxScene;
55+let params: Params = defaultParams(model);
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 gridNx = 128;
60+let Ls = 0.6;
61+/** Timesteps per display frame. May be fractional — ¼× means one step every
62+ * fourth frame — which is what `stepDebt` accumulates toward. */
63+let stepsPerFrame = 16;
64+let stepDebt = 0;
65+/** Where the microphone sits, in metres. Off to the side and above, out in
66+ * the air the instrument radiates into. */
67+const mic = { x: 0.12, y: 0.08, z: 0.1 };
68+/** Paused until asked. The page compiles and draws its silent initial state
69+ * on load — the string drawn back into its pluck — and nothing moves until
70+ * Run (or Render note). */
71+let running = false;
72+let session: ModelSession | null = null;
73+
74+const camera: Camera = { az: -2.2, el: 0.45, dist: 1.35 };
75+const viewState = { opacity: 2, contrast: 1.6, quality: 192, clipFrac: 0.5, exaggerate: 25 };
76+const show = { field: true, body: true, string: true, wire: true, mic: true };
77+
78+/** Pressure the colormap saturates at, and whether it follows the field. */
79+let scale = 1e-6;
80+let autoScale = true;
81+/** Largest pressure seen since the last pluck. The colour scale is not
82+ * allowed to fall far below it, so that once the note has faded what is
83+ * drawn is quiet air rather than roundoff at full contrast. */
84+let peakSeen = 0;
85+let colormapName = colormapNames[0];
86+/** The last string displacement read back, for the plot and the overlay. */
87+let lastString: Float32Array | null = null;
88+/** Something on screen would change if we drew now. A running simulation is
89+ * dirty every frame by definition; a paused one only when told. */
90+let dirty = true;
91+/** A note render in progress (the flag doubles as its cancel switch). */
92+let renderingNote = false;
93+
94+/* ------------------------------------------------------------- the page -- */
95+
96+const device = await requestAcousticDevice().catch((e: unknown) => {
97+ showError(e, '');
98+ return null;
99+});
100+if (!device) throw new Error('no GPU');
101+
102+// A lost device takes everything with it and nothing afterwards will work, so
103+// say so rather than leaving a frozen picture and no explanation.
104+void device.lost.then((info) => {
105+ showError(
106+ new Error(
107+ `the GPU device was lost (${info.reason}): ${info.message}\n` +
108+ 'Reload the page to start again.',
109+ ),
110+ '',
111+ );
112+});
113+
114+const canvas = el<HTMLCanvasElement>('view');
115+const view = new VolumeView(device, canvas);
116+view.setColormap(colormaps[colormapName]);
117+const overlay = new Overlay(el<HTMLCanvasElement>('overlay'));
118+const plot = new StringPlot(el<HTMLCanvasElement>('stringplot'));
119+const colorbar = new Colorbar(el('colorbar'));
120+colorbar.setColormap(colormaps[colormapName]);
121+const darkMedia = matchMedia('(prefers-color-scheme: dark)');
122+
123+new ResizeObserver(() => {
124+ view.resize();
125+ overlay.resize();
126+ plot.resize();
127+ dirty = true;
128+}).observe(canvas);
129+
130+const editor = new CodeEditor({
131+ textarea: el<HTMLTextAreaElement>('source'),
132+ overlay: el('highlight'),
133+ external: new Set([...EXTERNAL_OPS.keys(), 'sponge3']),
134+ onInput: (value) => {
135+ sources[editorFile.value as 'model' | 'scene'] = value;
136+ el<HTMLButtonElement>('recompile').classList.add('primary');
137+ },
138+});
139+const editorFile = el<HTMLSelectElement>('editor-file');
140+
141+const colormapSelect = el<HTMLSelectElement>('colormap');
142+for (const name of colormapNames) colormapSelect.append(new Option(name, name));
143+colormapSelect.value = colormapName;
144+
145+/* ------------------------------------------------------------- controls -- */
146+
147+/** One slider per parameter the registry declares. */
148+function buildSliders(
149+ host: HTMLElement,
150+ specs: ParamSpec[],
151+ values: Params,
152+ onChange: (key: string, value: number) => void,
153+): Map<string, (value: number) => void> {
154+ const setters = new Map<string, (value: number) => void>();
155+ host.textContent = '';
156+ for (const spec of specs) {
157+ const row = document.createElement('label');
158+ row.className = 'slider';
159+ if (spec.hint) row.title = spec.hint;
160+ const name = document.createElement('span');
161+ name.textContent = spec.label;
162+ const input = document.createElement('input');
163+ input.type = 'range';
164+ input.min = String(spec.min);
165+ input.max = String(spec.max);
166+ input.step = String(spec.step);
167+ input.value = String(values[spec.key]);
168+ const out = document.createElement('output');
169+ out.textContent = fmtValue(values[spec.key]);
170+ input.addEventListener('input', () => {
171+ const v = Number(input.value);
172+ out.textContent = fmtValue(v);
173+ onChange(spec.key, v);
174+ });
175+ row.append(name, input, out);
176+ host.append(row);
177+ setters.set(spec.key, (value) => {
178+ input.value = String(value);
179+ out.textContent = fmtValue(value);
180+ });
181+ }
182+ return setters;
183+}
184+
185+/** Scene edits go through the MATLAB interpreter, which is fast but not free;
186+ * a drag should not queue up one evaluation per pixel. */
187+let sceneTimer = 0;
188+const scheduleSceneUpdate = (): void => {
189+ clearTimeout(sceneTimer);
190+ sceneTimer = window.setTimeout(applyScene, 150);
191+};
192+
193+function applyScene(): void {
194+ if (!session) return;
195+ try {
196+ session.setScene(scene, sceneParams, sources.scene);
197+ clearError();
198+ // A scene edit that raises the fastest speed invalidates the timestep
199+ // (and with it the string grid), which only a rebuild can fix.
200+ if (Math.abs(session.dtWanted - session.dt) > 1e-12 * session.dt) {
201+ void rebuild();
202+ return;
203+ }
204+ dirty = true;
205+ } catch (e) {
206+ showError(e, sources.scene);
207+ }
208+}
209+
210+let micSliders = new Map<string, (value: number) => void>();
211+
212+function buildMicControls(): void {
213+ const specs: ParamSpec[] = [
214+ { key: 'x', label: 'mic x (m)', value: mic.x, min: -0.45, max: 0.45, step: 0.01 },
215+ { key: 'y', label: 'mic y (m)', value: mic.y, min: -0.22, max: 0.22, step: 0.01 },
216+ { key: 'z', label: 'mic z (m)', value: mic.z, min: -0.22, max: 0.22, step: 0.01 },
217+ ];
218+ micSliders = buildSliders(
219+ el('micparams'),
220+ specs,
221+ { x: mic.x, y: mic.y, z: mic.z },
222+ (key, value) => {
223+ mic[key as 'x' | 'y' | 'z'] = value;
224+ session?.setMic(mic.x, mic.y, mic.z);
225+ dirty = true;
226+ },
227+ );
228+}
229+
230+function buildParamControls(): void {
231+ buildSliders(el('params'), model.params, params, (key, value) => {
232+ params = { ...params, [key]: value };
233+ session?.setParams(params);
234+ });
235+ el('scene-title').textContent = `body — ${scene.blurb}`;
236+ buildSliders(el('sceneparams'), scene.params, sceneParams, (key, value) => {
237+ sceneParams = { ...sceneParams, [key]: value };
238+ scheduleSceneUpdate();
239+ });
240+}
241+
242+/* ------------------------------------------------------------- the loop -- */
243+
244+let frames = 0;
245+let lastFpsAt = performance.now();
246+let msPerFrame = 0;
247+/** A readback in flight; only one at a time, since they share a buffer. */
248+let reading = false;
249+
250+function pluck(): void {
251+ if (!session) return;
252+ session.pluck();
253+ peakSeen = 0;
254+ if (autoScale) scale = 1e-6;
255+ dirty = true;
256+ showRecording();
257+}
258+
259+/**
260+ * Compile the current sources into a new session and swap it in.
261+ *
262+ * The old session keeps running until the new one exists, and is only torn
263+ * down once the swap has happened: a compile that fails leaves something on
264+ * screen and something to edit rather than a dead page, and destroying GPU
265+ * resources while the next lot of shaders are still compiling is exactly the
266+ * sort of thing a browser is entitled to answer with a lost device.
267+ */
268+let building = false;
269+async function rebuild(): Promise<void> {
270+ if (building) return;
271+ building = true;
272+ renderingNote = false;
273+ const old = session;
274+ try {
275+ const next = await ModelSession.create({
276+ device: device!,
277+ model,
278+ params,
279+ source: sources.model,
280+ scene,
281+ sceneParams,
282+ sceneSource: sources.scene,
283+ nx: gridNx,
284+ Ls,
285+ });
286+ next.pluck();
287+ session = next;
288+ view.setSource(
289+ next.gpu.stateBuffer(next.pressureName)!,
290+ next.gpu.stateBuffer('wall')!,
291+ {
292+ nx: next.air.nx, ny: next.air.ny, nz: next.air.nz,
293+ Lx: next.air.Lx, Ly: next.air.Ly, Lz: next.air.Lz, h: next.air.h,
294+ },
295+ );
296+ next.setMic(mic.x, mic.y, mic.z);
297+ old?.destroy();
298+ lastString = null;
299+ peakSeen = 0;
300+ scale = 1e-6;
301+ dirty = true;
302+ clearError();
303+ showRecording();
304+ showStatics();
305+ el<HTMLButtonElement>('recompile').classList.remove('primary');
306+ el('compiled').textContent = describe(next);
307+ } catch (e) {
308+ // Which file the failure belongs to decides which one to show it against.
309+ const which = failingFile(e);
310+ showError(e, sources[which]);
311+ if (e instanceof ModelCompileError && e.start !== undefined) {
312+ editorFile.value = which;
313+ showFile(which);
314+ editor.select(e.start, e.end ?? e.start + 1);
315+ }
316+ } finally {
317+ building = false;
318+ }
319+}
320+
321+/** A scene failure is reported against `medium`, everything else against the
322+ * model's own functions. */
323+const failingFile = (e: unknown): 'model' | 'scene' =>
324+ e instanceof ModelCompileError && e.fn === 'medium' ? 'scene' : 'model';
325+
326+const describe = (s: ModelSession): string => {
327+ const { init, step } = s.describe();
328+ return [
329+ `% init — one pluck`,
330+ ...init.map((l) => ` ${l}`),
331+ ``,
332+ `% step — every timestep`,
333+ ...step.map((l) => ` ${l}`),
334+ ].join('\n');
335+};
336+
337+/** What the microphone has, and what it would sound like. */
338+function showRecording(): void {
339+ const info = el('recinfo');
340+ if (!session) {
341+ info.textContent = '';
342+ return;
343+ }
344+ const n = session.recorder.count;
345+ if (n === 0) {
346+ info.textContent = 'nothing recorded yet — Run, or Render note';
347+ return;
348+ }
349+ const plan = planPlayback(n, session.dt);
350+ info.textContent =
351+ `${fmtTime(n * session.dt)} recorded · plays at ` +
352+ `${Math.round(plan.rate).toLocaleString()} Hz` +
353+ (plan.realTime ? '' : ' (rate clamped)') +
354+ (session.recorder.full ? ' · buffer full' : '');
355+}
356+
357+/** The facts that only change on a rebuild. */
358+function showStatics(): void {
359+ if (!session) return;
360+ const { air, string } = session;
361+ const fmax = C_AIR / (POOR_RESOLUTION * air.h);
362+ el('domaininfo').textContent =
363+ `The air is a ${fmtLength(air.Lx)} × ${fmtLength(air.Ly)} × ${fmtLength(air.Lz)} box ` +
364+ `at ${air.nx}×${air.ny}×${air.nz} cells of ${fmtLength(air.h)} — honest to about ` +
365+ `${(fmax / 1000).toFixed(1)} kHz. The string has ${string.ns} nodes; ` +
366+ `dt = ${fmtTime(session.dt)} (${Math.round(1 / session.dt / 1000)} kHz).`;
367+}
368+
369+function stats(): void {
370+ if (!session) return;
371+ const rate =
372+ running && msPerFrame > 0 ? `${msPerFrame.toFixed(1)} ms/frame` : 'paused';
373+ const f = params.f0 ?? 0;
374+ const partials = Math.max(1, Math.floor(C_AIR / (POOR_RESOLUTION * session.air.h) / Math.max(f, 1)));
375+ el('stats').innerHTML =
376+ `t = <b>${fmtTime(session.t)}</b> · step <b>${session.steps.toLocaleString()}</b> · ` +
377+ `${fmtValue(f)} Hz fundamental — the air carries its first ~${partials} partials · ${rate}`;
378+}
379+
380+/**
381+ * Follow the field with the colour scale, and keep the string plot fed.
382+ *
383+ * The only readbacks in the app, a few times a second rather than every
384+ * frame, and strictly one at a time — they share a staging buffer.
385+ */
386+async function pollFields(): Promise<void> {
387+ if (!session || reading) return;
388+ reading = true;
389+ try {
390+ const u = await session.read(session.displacementName);
391+ lastString = u;
392+ if (autoScale) {
393+ const p = await session.read(session.pressureName);
394+ let peak = 0;
395+ for (const v of p) peak = Math.max(peak, Math.abs(v));
396+ peakSeen = Math.max(peakSeen, peak);
397+ scale = peak > scale ? peak : 0.97 * scale + 0.03 * peak;
398+ scale = Math.max(scale, 0.05 * peakSeen, 1e-12);
399+ }
400+ dirty = true;
401+ } catch {
402+ // A rebuild can destroy the buffers mid-read; the next poll recovers.
403+ } finally {
404+ reading = false;
405+ }
406+}
407+
408+let lastPollAt = 0;
409+function frame(): void {
410+ if (session && (running || dirty)) {
411+ if (running && !renderingNote) {
412+ // Fractional speeds accumulate a debt and step when it reaches a whole
413+ // timestep, so ¼× is one step every fourth frame rather than nothing.
414+ stepDebt += stepsPerFrame;
415+ const n = Math.floor(stepDebt);
416+ if (n > 0) {
417+ session.step(n);
418+ stepDebt -= n;
419+ }
420+ }
421+ const cf = cameraFrame(camera, canvas.width / Math.max(canvas.height, 1));
422+ view.draw({
423+ frame: cf,
424+ scale: Math.max(scale, 1e-20),
425+ opacity: show.field ? viewState.opacity : 0,
426+ contrast: viewState.contrast,
427+ steps: viewState.quality,
428+ clipY: viewState.clipFrac * session.air.Ly,
429+ body: show.body ? 0.6 : 0,
430+ mic: show.mic ? mic : null,
431+ });
432+ overlay.draw({
433+ frame: cf,
434+ geometry: {
435+ Ls: session.string.Ls,
436+ stringZ: scene.stringZ(sceneParams),
437+ boxl: sceneParams.boxl ?? 0,
438+ boxw: sceneParams.boxw ?? 0,
439+ boxd: sceneParams.boxd ?? 0,
440+ holer: sceneParams.holer ?? 0,
441+ holex: sceneParams.holex ?? 0,
442+ },
443+ string: show.string ? lastString : null,
444+ exaggerate: viewState.exaggerate,
445+ wireframe: show.wire,
446+ });
447+ plot.draw(lastString, 1.1 * (params.amp ?? 0.002), darkMedia.matches);
448+ el('stringlabel').textContent =
449+ `string displacement, full scale ±${fmtValue(1100 * (params.amp ?? 0.002))} mm`;
450+ dirty = false;
451+ colorbar.setRange(-scale, scale);
452+ frames++;
453+
454+ const now = performance.now();
455+ if ((running || renderingNote) && now - lastPollAt > 250) {
456+ lastPollAt = now;
457+ void pollFields();
458+ }
459+ }
460+ const now = performance.now();
461+ if (now - lastFpsAt > 400) {
462+ msPerFrame = frames > 0 ? (now - lastFpsAt) / frames : 0;
463+ frames = 0;
464+ lastFpsAt = now;
465+ stats();
466+ showRecording();
467+ }
468+ requestAnimationFrame(frame);
469+}
470+
471+/* ------------------------------------------------------- rendering a note -- */
472+
473+/**
474+ * Pluck, then run the solver as fast as the GPU will go — no display frames
475+ * in the way, just batches of steps with a sync between them so the queue
476+ * never runs unboundedly ahead — until the requested duration of audio is in
477+ * the trace. Then play it.
478+ */
479+async function renderNote(): Promise<void> {
480+ if (!session || building) return;
481+ const btn = el<HTMLButtonElement>('rendernote');
482+ if (renderingNote) {
483+ renderingNote = false; // cancel: the loop below notices and stops
484+ return;
485+ }
486+ renderingNote = true;
487+ setRunning(false);
488+ btn.textContent = 'Cancel';
489+ const prog = el<HTMLProgressElement>('renderprog');
490+ prog.hidden = false;
491+ prog.value = 0;
492+ const seconds = Number(el<HTMLSelectElement>('renderdur').value);
493+ try {
494+ pluck();
495+ const s = session;
496+ const total = Math.min(
497+ Math.ceil(seconds / s.dt),
498+ s.recorder.capacity,
499+ );
500+ const batch = 512;
501+ let done = 0;
502+ while (done < total && renderingNote && session === s) {
503+ const n = Math.min(batch, total - done);
504+ s.step(n);
505+ done += n;
506+ await s.sync();
507+ prog.value = done / total;
508+ }
509+ if (renderingNote && session === s) {
510+ const trace = await s.recorder.read();
511+ stopAudio();
512+ await playTrace(trace, planPlayback(trace.length, s.dt));
513+ clearError();
514+ }
515+ } catch (e) {
516+ showError(e, '');
517+ } finally {
518+ renderingNote = false;
519+ prog.hidden = true;
520+ btn.textContent = 'Render note';
521+ dirty = true;
522+ showRecording();
523+ }
524+}
525+
526+/* -------------------------------------------------------------- wiring --- */
527+
528+function showFile(which: 'model' | 'scene'): void {
529+ editor.value = sources[which];
530+ el('editor-title').textContent =
531+ which === 'model'
532+ ? 'init and step — string and air together, compiled to WebGPU'
533+ : 'medium(x, y, z, …) → walls and coupling, evaluated once on the CPU';
534+}
535+
536+el('gridsize').addEventListener('change', (e) => {
537+ gridNx = Number((e.target as HTMLSelectElement).value);
538+ void rebuild();
539+});
540+
541+el('stringlen').addEventListener('change', (e) => {
542+ Ls = Number((e.target as HTMLSelectElement).value);
543+ void rebuild();
544+});
545+
546+const runPause = el<HTMLButtonElement>('runpause');
547+const setRunning = (r: boolean): void => {
548+ running = r;
549+ runPause.textContent = r ? 'Pause' : 'Run';
550+ frames = 0;
551+ lastFpsAt = performance.now();
552+ dirty = true;
553+};
554+runPause.addEventListener('click', () => setRunning(!running && !renderingNote));
555+
556+el('pluck').addEventListener('click', pluck);
557+
558+el('spf').addEventListener('change', (e) => {
559+ stepsPerFrame = Number((e.target as HTMLSelectElement).value);
560+});
561+
562+el('rendernote').addEventListener('click', () => void renderNote());
563+
564+const listen = el<HTMLButtonElement>('listen');
565+listen.addEventListener('click', () => {
566+ if (!session) return;
567+ const s = session;
568+ if (s.recorder.count === 0) {
569+ showError(new Error('the microphone has not recorded anything yet — Run, or Render note'), '');
570+ return;
571+ }
572+ listen.disabled = true;
573+ stopAudio();
574+ void s.recorder
575+ .read()
576+ .then((trace) => playTrace(trace, planPlayback(trace.length, s.dt)))
577+ .then(() => clearError())
578+ .catch((e: unknown) => showError(e, ''))
579+ .finally(() => {
580+ listen.disabled = false;
581+ });
582+});
583+
584+el('download').addEventListener('click', () => {
585+ if (!session) return;
586+ const s = session;
587+ if (s.recorder.count === 0) {
588+ showError(new Error('the microphone has not recorded anything yet — Run, or Render note'), '');
589+ return;
590+ }
591+ void s.recorder
592+ .read()
593+ .then((trace) => {
594+ const url = URL.createObjectURL(traceToWav(trace, s.dt));
595+ const a = document.createElement('a');
596+ a.href = url;
597+ a.download = 'dulcimer.wav';
598+ a.click();
599+ setTimeout(() => URL.revokeObjectURL(url), 10_000);
600+ })
601+ .catch((e: unknown) => showError(e, ''));
602+});
603+
604+colormapSelect.addEventListener('change', () => {
605+ colormapName = colormapSelect.value;
606+ view.setColormap(colormaps[colormapName]);
607+ colorbar.setColormap(colormaps[colormapName]);
608+ dirty = true;
609+});
610+
611+el('scalemode').addEventListener('change', (e) => {
612+ autoScale = (e.target as HTMLSelectElement).value === 'auto';
613+});
614+
615+const bindRange = (id: string, apply: (v: number) => void): void => {
616+ el<HTMLInputElement>(id).addEventListener('input', (e) => {
617+ apply(Number((e.target as HTMLInputElement).value));
618+ dirty = true;
619+ });
620+};
621+bindRange('opacity', (v) => (viewState.opacity = v));
622+bindRange('contrast', (v) => (viewState.contrast = v));
623+bindRange('clipy', (v) => (viewState.clipFrac = v));
624+bindRange('exaggerate', (v) => (viewState.exaggerate = v));
625+el('quality').addEventListener('change', (e) => {
626+ viewState.quality = Number((e.target as HTMLSelectElement).value);
627+ dirty = true;
628+});
629+
630+const bindCheck = (id: string, apply: (v: boolean) => void): void => {
631+ el<HTMLInputElement>(id).addEventListener('change', (e) => {
632+ apply((e.target as HTMLInputElement).checked);
633+ dirty = true;
634+ });
635+};
636+bindCheck('showfield', (v) => (show.field = v));
637+bindCheck('showbody', (v) => (show.body = v));
638+bindCheck('showstring', (v) => (show.string = v));
639+bindCheck('showwire', (v) => (show.wire = v));
640+bindCheck('showmic', (v) => (show.mic = v));
641+
642+// --- orbit ---------------------------------------------------------------
643+let dragging = false;
644+let last = [0, 0];
645+canvas.addEventListener('pointerdown', (e) => {
646+ dragging = true;
647+ last = [e.clientX, e.clientY];
648+ canvas.setPointerCapture(e.pointerId);
649+});
650+canvas.addEventListener('pointermove', (e) => {
651+ if (!dragging) return;
652+ camera.az -= (e.clientX - last[0]) * 0.008;
653+ camera.el = Math.min(1.5, Math.max(-1.5, camera.el + (e.clientY - last[1]) * 0.008));
654+ last = [e.clientX, e.clientY];
655+ dirty = true;
656+});
657+canvas.addEventListener('pointerup', () => (dragging = false));
658+canvas.addEventListener('pointercancel', () => (dragging = false));
659+canvas.addEventListener(
660+ 'wheel',
661+ (e) => {
662+ e.preventDefault();
663+ camera.dist = Math.min(5, Math.max(0.7, camera.dist * Math.exp(e.deltaY * 0.001)));
664+ dirty = true;
665+ },
666+ { passive: false },
667+);
668+
669+editorFile.addEventListener('change', () => {
670+ showFile(editorFile.value as 'model' | 'scene');
671+});
672+
673+el('recompile').addEventListener('click', () => {
674+ void rebuild();
675+});
676+
677+el('revert').addEventListener('click', () => {
678+ const which = editorFile.value as 'model' | 'scene';
679+ sources[which] = which === 'model' ? model.source : scene.source;
680+ showFile(which);
681+ void rebuild();
682+});
683+
684+/* ---------------------------------------------------------------- start -- */
685+
686+buildParamControls();
687+buildMicControls();
688+showFile('model');
689+darkMedia.addEventListener('change', () => (dirty = true));
690+await rebuild();
691+// The initial pluck shape, so the page opens showing the string drawn back.
692+void pollFields();
693+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+110−0View file
@@ -0,0 +1,110 @@
1+/**
2+ * The host-provided operations a model .m can call — the only places where a
3+ * point reads anything but itself.
4+ *
5+ * Everything else a model does is element-wise, so these five are where the
6+ * physics that couples neighbours (and couples the two grids) lives:
7+ *
8+ * dxx(u) second difference along the string (string -> string)
9+ * dxxxx(u) fourth difference along the string (string -> string)
10+ * lapw(p, w) wall-masked 7-point Laplacian in the air (air, air -> air)
11+ * spread(a) sample a string field at each air point's x (string -> air)
12+ * bridge(u) du/dx at the string's bridge end, broadcast (string -> air)
13+ *
14+ * numbl needs only their *type rule* in order to lower a call site. It gets
15+ * that from a `.mtoc2.js` workspace file — numbl's sanctioned extension point
16+ * for a JS-defined builtin (see `mtoc2UserFunctionsByName` in numbl's
17+ * LoweringContext). The file is evaluated in a bare CommonJS sandbox with no
18+ * imports available, so `transfer` builds numbl `Type` objects as plain
19+ * literals, and the grid sizes are baked in by the generator below (a grid
20+ * change recompiles anyway).
21+ *
22+ * The `emit`/`cBody` exports exist only because the loader's contract requires
23+ * them; we never emit C. The actual implementation is supplied by the WGSL
24+ * backend (src/mgpu/ops.ts).
25+ */
26+
27+export interface GridSizes {
28+ /** Air grid points, nx*ny*nz. Air fields are npts x 1 column vectors. */
29+ npts: number;
30+ /** String nodes. String fields are ns x 1 column vectors. */
31+ ns: number;
32+}
33+
34+/** One external operation's shape contract. */
35+export interface ExternalOpSpec {
36+ name: string;
37+ /** Element counts of the arguments, in order. */
38+ args: ('air' | 'string')[];
39+ out: 'air' | 'string';
40+}
41+
42+export const EXTERNAL_OP_SPECS: ExternalOpSpec[] = [
43+ { name: 'dxx', args: ['string'], out: 'string' },
44+ { name: 'dxxxx', args: ['string'], out: 'string' },
45+ { name: 'lapw', args: ['air', 'air'], out: 'air' },
46+ { name: 'spread', args: ['string'], out: 'air' },
47+ { name: 'bridge', args: ['string'], out: 'air' },
48+];
49+
50+/** Names the WGSL backend must implement as dispatches rather than
51+ * element-wise kernels, with the argument count of each. */
52+export const EXTERNAL_OPS = new Map(EXTERNAL_OP_SPECS.map((s) => [s.name, s]));
53+
54+const numericType = (rows: number): string =>
55+ `{ kind: "Numeric", elem: "double", isComplex: false, ` +
56+ `dims: [{ kind: "exact", value: ${rows} }, { kind: "exact", value: 1 }], ` +
57+ `shape: [${rows}, 1], sign: "unknown" }`;
58+
59+/** Source for one op's `.mtoc2.js`: fields in, field out, shapes checked. */
60+function opSource(spec: ExternalOpSpec, g: GridSizes): string {
61+ const size = (k: 'air' | 'string'): number => (k === 'air' ? g.npts : g.ns);
62+ const checks = spec.args
63+ .map((kind, i) => {
64+ const n = size(kind);
65+ const what = kind === 'air' ? 'an air field' : 'a string field';
66+ return `
67+ var a${i} = argTypes[${i}];
68+ if (!a${i} || a${i}.kind !== "Numeric" || a${i}.isComplex) {
69+ throw new Error("${spec.name}: argument ${i + 1} must be a real numeric array");
70+ }
71+ var s${i} = a${i}.shape;
72+ if (!s${i} || s${i}.length !== 2 || s${i}[0] !== ${n} || s${i}[1] !== 1) {
73+ throw new Error(
74+ "${spec.name}: argument ${i + 1} must be ${what} (${n}x1), not " +
75+ (s${i} ? s${i}.join("x") : "unknown shape")
76+ );
77+ }`;
78+ })
79+ .join('\n');
80+ return `
81+exports.name = ${JSON.stringify(spec.name)};
82+
83+exports.transfer = function (argTypes, nargout) {
84+ if (argTypes.length !== ${spec.args.length}) {
85+ throw new Error("${spec.name} takes ${spec.args.length} argument(s), got " + argTypes.length);
86+ }
87+ if (nargout > 1) {
88+ throw new Error("${spec.name} returns one value, but " + nargout + " were requested");
89+ }
90+${checks}
91+ return [${numericType(size(spec.out))}];
92+};
93+
94+// Never called: this project executes the IR on WebGPU and emits no C.
95+exports.emit = function () {
96+ throw new Error("${spec.name}: no C backend (this op runs on WebGPU)");
97+};
98+exports.cBody = function () {
99+ return "";
100+};
101+`;
102+}
103+
104+/** Workspace files that make the ops resolvable during lowering. */
105+export function externalOpFiles(g: GridSizes): { name: string; source: string }[] {
106+ return EXTERNAL_OP_SPECS.map((spec) => ({
107+ name: `${spec.name}.mtoc2.js`,
108+ source: opSource(spec, g),
109+ }));
110+}
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+309−0View file
@@ -0,0 +1,309 @@
1+/**
2+ * The .m model, compiled and running on the GPU.
3+ *
4+ * The model file is ordinary MATLAB: it defines an `init` function that builds
5+ * the initial state — the plucked string and a silent air field — and a `step`
6+ * function that advances both one timestep. Each is specialized for the
7+ * current grids and compiled into a ModelPlan, and both operate on the same
8+ * state buffers (see HostBuffers).
9+ *
10+ * Both functions return the state, in the same order, so their signatures say
11+ * exactly what they produce:
12+ *
13+ * function [u, um, p, pm] = init(xs, npts, Ls, pluckpos, amp)
14+ * function [un, uold, pn, pold] = step(u, um, p, pm, ...)
15+ *
16+ * The state spans two grids — `u`, `um` are string fields (ns nodes), `p`,
17+ * `pm` air fields (npts cells) — and the compiled statements mix freely:
18+ * each line is one kernel over its own field's size, and the external ops
19+ * (src/mgpu/ops.ts) are where a value crosses from one grid to the other.
20+ *
21+ * The host supplies the things that are setup rather than algorithm: the grid
22+ * coordinates, the medium and coupling profiles the scene defines, the
23+ * timestep, and the parameter values. Each argument is matched to the .m's
24+ * declared parameter name, so the file documents its own interface.
25+ *
26+ * There is no clock here, and that is not an oversight: the pluck is an
27+ * initial condition, not a driven source, so nothing in the model needs to
28+ * know what time it is.
29+ */
30+import { HostBuffers, ModelPlan } from './plan.ts';
31+import type { OpPlan } from './ops.ts';
32+import { inFunction, inFunctionAsync, inModel } from './errors.ts';
33+import { CompiledModel, type Binding } from './compile.ts';
34+import type { AirGrid, StringGrid } from '../grid.ts';
35+
36+export interface ModelParams {
37+ [key: string]: number;
38+}
39+
40+/** What the scene defines, on the air grid: npts values each. */
41+export interface MediumFields {
42+ /** Sound speed, m/s. */
43+ c: Float32Array;
44+ /** Absorption rate, 1/s (the sponge and any wall absorption). */
45+ sig: Float32Array;
46+ /** 1 in air, 0 in the body's solid shell. What `lapw` masks by. */
47+ wall: Float32Array;
48+ /** Where the string radiates directly: a tube around the string line. */
49+ lineprof: Float32Array;
50+ /** Where the bridge force drives the air: a patch above the top plate. */
51+ boardprof: Float32Array;
52+}
53+
54+/** One state field the .m advances, and which grid it lives on. */
55+export interface StateField {
56+ name: string;
57+ grid: 'air' | 'string';
58+}
59+
60+export interface GpuModelOptions {
61+ device: GPUDevice;
62+ ops: OpPlan;
63+ air: AirGrid;
64+ string: StringGrid;
65+ medium: MediumFields;
66+ /** Model source (.m text). */
67+ source: string;
68+ /** Parameter names the .m may take as arguments. */
69+ paramNames: string[];
70+ /** State fields the .m advances, in the order its functions return them. */
71+ state: StateField[];
72+ /** Grid fields one kernel may read, overriding what the device allows.
73+ * Only for tests. */
74+ operandBudget?: number;
75+}
76+
77+/** Names the .m may take for the air grid coordinates. */
78+export const AIR_GRID_NAMES = ['x', 'y', 'z'] as const;
79+/** Names the .m may take for the string grid: node positions and the pin
80+ * mask that terminates the ends. */
81+export const STRING_GRID_NAMES = ['xs', 'pin'] as const;
82+/** Names the .m may take for what the scene defines. */
83+export const MEDIUM_NAMES = ['c', 'sig', 'wall', 'lineprof', 'boardprof'] as const;
84+
85+export class GpuModel {
86+ readonly paramNames: string[];
87+ readonly state: StateField[];
88+ readonly npts: number;
89+ readonly ns: number;
90+
91+ #device: GPUDevice;
92+ #host: HostBuffers;
93+ #initPlan: ModelPlan;
94+ #stepPlan: ModelPlan;
95+ /** Timestep, host-owned: it follows from the grid and the medium (a CFL
96+ * condition), not from anything the user types, and it is folded into every
97+ * setParams so the .m's `dt` is never left with the zero a missing
98+ * parameter would default to. */
99+ #dt = 0;
100+ #readback: GPUBuffer;
101+ /** Which function wrote the state most recently; see `read`. */
102+ #lastRan: 'init' | 'step' = 'init';
103+
104+ private constructor(init: {
105+ device: GPUDevice;
106+ host: HostBuffers;
107+ initPlan: ModelPlan;
108+ stepPlan: ModelPlan;
109+ readback: GPUBuffer;
110+ paramNames: string[];
111+ state: StateField[];
112+ npts: number;
113+ ns: number;
114+ }) {
115+ this.#device = init.device;
116+ this.#host = init.host;
117+ this.#initPlan = init.initPlan;
118+ this.#stepPlan = init.stepPlan;
119+ this.#readback = init.readback;
120+ this.paramNames = init.paramNames;
121+ this.state = init.state;
122+ this.npts = init.npts;
123+ this.ns = init.ns;
124+ }
125+
126+ static async create(opts: GpuModelOptions): Promise<GpuModel> {
127+ const { device, ops, air, string, medium, source, paramNames, state } = opts;
128+ const npts = air.npts;
129+ const ns = string.ns;
130+ const sizeOf = (grid: 'air' | 'string'): number => (grid === 'air' ? npts : ns);
131+
132+ // What the .m may ask for by name. The grid geometry is exact, so a
133+ // constructor reading it (`zeros(npts, 1)`) keeps a static shape.
134+ const bindings: Record<string, Binding> = {
135+ npts: { kind: 'const', value: npts },
136+ nx: { kind: 'const', value: air.nx },
137+ ny: { kind: 'const', value: air.ny },
138+ nz: { kind: 'const', value: air.nz },
139+ h: { kind: 'const', value: air.h },
140+ ns: { kind: 'const', value: ns },
141+ hs: { kind: 'const', value: string.hs },
142+ Ls: { kind: 'const', value: string.Ls },
143+ dt: { kind: 'param' },
144+ };
145+ for (const g of AIR_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
146+ for (const g of STRING_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [ns, 1] };
147+ for (const m of MEDIUM_NAMES) bindings[m] = { kind: 'tensor', shape: [npts, 1] };
148+ for (const s of state) bindings[s.name] = { kind: 'tensor', shape: [sizeOf(s.grid), 1] };
149+ for (const p of paramNames) bindings[p] = { kind: 'param' };
150+
151+ // Parsing belongs to the file, not to either function.
152+ const compiled = inModel(() => new CompiledModel(source, bindings, { npts, ns }));
153+ const nargout = state.length;
154+ const initFn = inFunction('init', () => compiled.specialize('init', nargout));
155+ const stepFn = inFunction('step', () => compiled.specialize('step', nargout));
156+ compiled.finish();
157+
158+ // Both functions return the state, in order, and both feed it back into
159+ // the shared buffers.
160+ const feedback = state.map((s) => s.name);
161+
162+ const host = new HostBuffers(device);
163+ // The host owns the state and the inputs it uploads, whether or not a
164+ // given function happens to take them as arguments.
165+ for (const s of state) host.ensure(s.name, sizeOf(s.grid));
166+ for (const g of AIR_GRID_NAMES) host.ensure(g, npts);
167+ for (const g of STRING_GRID_NAMES) host.ensure(g, ns);
168+ for (const m of MEDIUM_NAMES) host.ensure(m, npts);
169+
170+ const initPlan = await inFunctionAsync('init', () =>
171+ ModelPlan.create(device, ops, { fn: initFn, feedback }, host, opts.operandBudget),
172+ );
173+ const stepPlan = await inFunctionAsync('step', () =>
174+ ModelPlan.create(device, ops, { fn: stepFn, feedback }, host, opts.operandBudget),
175+ );
176+
177+ host.upload('x', air.x);
178+ host.upload('y', air.y);
179+ host.upload('z', air.z);
180+ host.upload('xs', string.xs);
181+ host.upload('pin', string.pin);
182+ for (const m of MEDIUM_NAMES) host.upload(m, medium[m]);
183+
184+ const readback = device.createBuffer({
185+ label: 'mgpu-readback',
186+ size: 4 * Math.max(npts, ns),
187+ usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
188+ });
189+
190+ return new GpuModel({
191+ device, host, initPlan, stepPlan, readback, paramNames, state, npts, ns,
192+ });
193+ }
194+
195+ /** The timestep in force. Host-owned; see `#dt`. */
196+ get dt(): number {
197+ return this.#dt;
198+ }
199+
200+ setDt(dt: number): void {
201+ this.#dt = dt;
202+ }
203+
204+ setParams(params: ModelParams): void {
205+ const merged = { dt: this.#dt, ...params };
206+ this.#initPlan.setParams(merged);
207+ this.#stepPlan.setParams(merged);
208+ }
209+
210+ /**
211+ * Swap the medium under a running model. It is data, not code — its shape in
212+ * the bindings depends only on the grid — so changing the body's geometry is
213+ * five buffer writes and needs no recompile.
214+ */
215+ uploadMedium(medium: MediumFields): void {
216+ for (const m of MEDIUM_NAMES) this.#host.upload(m, medium[m]);
217+ }
218+
219+ /** Write a host-owned value directly. Lets a test set up an exact initial
220+ * condition instead of going through `init`. */
221+ upload(name: string, data: Float32Array): void {
222+ this.#host.upload(name, data);
223+ }
224+
225+ /** Run `init`, replacing the state. */
226+ init(): void {
227+ const enc = this.#device.createCommandEncoder({ label: 'mgpu-init' });
228+ this.#initPlan.encodeSteps(enc, 1);
229+ this.#device.queue.submit([enc.finish()]);
230+ this.#lastRan = 'init';
231+ }
232+
233+ /**
234+ * Advance `steps` timesteps. Synchronous — this only records commands and
235+ * submits them; nothing is read back and nothing is awaited.
236+ *
237+ * `after` is recorded once per step, so anything that must see every
238+ * timestep (the microphone) rides along in the same submission.
239+ */
240+ step(steps = 1, after?: (encoder: GPUCommandEncoder) => void): void {
241+ const enc = this.#device.createCommandEncoder({ label: 'mgpu-step' });
242+ this.#stepPlan.encodeSteps(enc, steps, after);
243+ this.#device.queue.submit([enc.finish()]);
244+ this.#lastRan = 'step';
245+ }
246+
247+ /**
248+ * The buffer currently holding a named value. A field the .m computes is
249+ * produced by both functions, into separate buffers (only the state is
250+ * shared), so this resolves to whichever function ran most recently — which
251+ * is what makes the first frame show the initial state rather than an
252+ * unwritten buffer.
253+ */
254+ #locate(name: string): { buffer: GPUBuffer; count: number } | null {
255+ const [first, second] =
256+ this.#lastRan === 'init'
257+ ? [this.#initPlan, this.#stepPlan]
258+ : [this.#stepPlan, this.#initPlan];
259+ const buffer = first.buffer(name) ?? second.buffer(name);
260+ const count = first.elementCount(name) ?? second.elementCount(name);
261+ if (!buffer || count === undefined) return null;
262+ return { buffer, count };
263+ }
264+
265+ /** The GPU buffer a named value would be read from right now. */
266+ valueBuffer(name: string): GPUBuffer | null {
267+ return this.#locate(name)?.buffer ?? null;
268+ }
269+
270+ /**
271+ * The buffer a host-owned field lives in — the state between calls, or an
272+ * input like the wall mask.
273+ *
274+ * This is what the renderer binds, and it must be this rather than
275+ * `valueBuffer`: a bind group is built once and holds a particular buffer,
276+ * while `init` and `step` write their outputs into buffers of their own and
277+ * only agree here, where their feedback copies land. Binding either
278+ * function's private buffer would draw a stale field for half the run.
279+ */
280+ stateBuffer(name: string): GPUBuffer | null {
281+ return this.#host.get(name)?.buffer ?? null;
282+ }
283+
284+ /** Read a named value back to the CPU. The only await in the whole loop. */
285+ async read(name: string): Promise<Float32Array> {
286+ const located = this.#locate(name);
287+ if (!located) throw new Error(`read: the model has no value named '${name}'`);
288+ const { buffer, count } = located;
289+ const enc = this.#device.createCommandEncoder({ label: `mgpu-read-${name}` });
290+ enc.copyBufferToBuffer(buffer, 0, this.#readback, 0, 4 * count);
291+ this.#device.queue.submit([enc.finish()]);
292+ await this.#readback.mapAsync(GPUMapMode.READ, 0, 4 * count);
293+ const out = new Float32Array(this.#readback.getMappedRange(0, 4 * count).slice(0));
294+ this.#readback.unmap();
295+ return out;
296+ }
297+
298+ /** What the .m compiled to, for display. */
299+ describe(): { init: string[]; step: string[] } {
300+ return { init: this.#initPlan.describe(), step: this.#stepPlan.describe() };
301+ }
302+
303+ destroy(): void {
304+ this.#initPlan.destroy();
305+ this.#stepPlan.destroy();
306+ this.#host.destroy();
307+ this.#readback.destroy();
308+ }
309+}
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/ops.tsadded+250−0View file
@@ -0,0 +1,250 @@
1+/**
2+ * The external operations, as GPU dispatches.
3+ *
4+ * An air field is an npts x 1 column vector laid out x-fastest — the point
5+ * (ix, iy, iz) is element `ix + nx*(iy + ny*iz)` — and a string field is ns
6+ * nodes from x = 0 to x = Ls. Those layouts, and the placement of the string
7+ * inside the air's coordinates, matter only here and to the renderer; a .m
8+ * never sees them, which is the reason these are host-provided operations
9+ * rather than something a model expresses with array slicing.
10+ *
11+ * `lapw` is the wall-masked Laplacian, and it is what makes the body's shell
12+ * rigid. Each of the six face terms is scaled by the mask at the neighbour it
13+ * reads, so a face into the wall contributes nothing — which is exactly the
14+ * discrete Neumann (zero normal velocity, sound-hard) condition, in the
15+ * divergence form div(w grad p). The alternative the flat siblings use, a
16+ * wall as a fast material, is unusable here: wood at ~4000 m/s would cut the
17+ * global timestep twelve-fold, where a mask costs nothing. Outside the domain
18+ * the field is taken to be zero (a soft outer boundary); the scene's
19+ * absorbing layer is meant to have swallowed the wave before it matters.
20+ *
21+ * `spread` and `bridge` are the two couplings from the string into the air.
22+ * `spread` gives each air point the string value at its own x (linearly
23+ * interpolated, zero beyond the ends): multiplied by a scene-built line
24+ * profile, that is the string radiating directly. `bridge` broadcasts the
25+ * string's slope at its bridge end to every air point: multiplied by a
26+ * scene-built patch profile, that is the bridge force driving the top plate.
27+ * Both read tiny buffers and write big ones, so each is one cheap dispatch.
28+ */
29+import { UnsupportedOnGpu, WORKGROUP_SIZE } from './wgsl.ts';
30+import { EXTERNAL_OPS } from './externals.ts';
31+
32+export type OpKind = 'dxx' | 'dxxxx' | 'lapw' | 'spread' | 'bridge';
33+
34+export interface OpGeometry {
35+ /** Air grid. */
36+ nx: number;
37+ ny: number;
38+ nz: number;
39+ h: number;
40+ /** String grid. */
41+ ns: number;
42+ hs: number;
43+ /** World x of the string's first node (the nut); the bridge is the last. */
44+ xs0: number;
45+ /** World x extent of the air domain, for mapping voxel index to metres. */
46+ Lx: number;
47+}
48+
49+function opWGSL(kind: OpKind, g: OpGeometry): { code: string; outCount: number; args: number } {
50+ const { nx, ny, nz, h, ns, hs } = g;
51+ const npts = nx * ny * nz;
52+ const inv2 = 1 / (hs * hs);
53+ const inv4 = 1 / (hs * hs * hs * hs);
54+
55+ const stringAt = `
56+fn at(j: i32) -> f32 {
57+ if (j < 0 || j >= ${ns}) { return 0.0; }
58+ return src[u32(j)];
59+}
60+`;
61+
62+ switch (kind) {
63+ case 'dxx':
64+ return {
65+ args: 1,
66+ outCount: ns,
67+ code: `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
68+@group(0) @binding(1) var<storage, read> src: array<f32>;
69+${stringAt}
70+@compute @workgroup_size(${WORKGROUP_SIZE})
71+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
72+ let i = gid.x;
73+ if (i >= ${ns}u) { return; }
74+ let j = i32(i);
75+ dst[i] = ${inv2}f * (at(j - 1) - 2.0 * at(j) + at(j + 1));
76+}
77+`,
78+ };
79+
80+ case 'dxxxx':
81+ return {
82+ args: 1,
83+ outCount: ns,
84+ code: `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
85+@group(0) @binding(1) var<storage, read> src: array<f32>;
86+${stringAt}
87+@compute @workgroup_size(${WORKGROUP_SIZE})
88+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
89+ let i = gid.x;
90+ if (i >= ${ns}u) { return; }
91+ let j = i32(i);
92+ dst[i] = ${inv4}f * (at(j - 2) - 4.0 * at(j - 1) + 6.0 * at(j) - 4.0 * at(j + 1) + at(j + 2));
93+}
94+`,
95+ };
96+
97+ case 'lapw':
98+ return {
99+ args: 2,
100+ outCount: npts,
101+ code: `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
102+@group(0) @binding(1) var<storage, read> src: array<f32>;
103+@group(0) @binding(2) var<storage, read> wall: array<f32>;
104+
105+fn idx(ix: i32, iy: i32, iz: i32) -> i32 {
106+ return ix + ${nx} * (iy + ${ny} * iz);
107+}
108+fn inside(ix: i32, iy: i32, iz: i32) -> bool {
109+ return ix >= 0 && ix < ${nx} && iy >= 0 && iy < ${ny} && iz >= 0 && iz < ${nz};
110+}
111+// One face's contribution: masked by the wall at the neighbour, so a face
112+// into the wall carries no flux — the Neumann (rigid) condition.
113+fn face(ix: i32, iy: i32, iz: i32, pc: f32) -> f32 {
114+ if (!inside(ix, iy, iz)) { return -pc; } // outside the domain: p = 0, open
115+ let k = u32(idx(ix, iy, iz));
116+ return wall[k] * (src[k] - pc);
117+}
118+
119+@compute @workgroup_size(${WORKGROUP_SIZE})
120+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
121+ let i = gid.x;
122+ if (i >= ${npts}u) { return; }
123+ let ix = i32(i % ${nx}u);
124+ let iy = i32((i / ${nx}u) % ${ny}u);
125+ let iz = i32(i / ${nx * ny}u);
126+ let pc = src[i];
127+ let s = face(ix - 1, iy, iz, pc) + face(ix + 1, iy, iz, pc)
128+ + face(ix, iy - 1, iz, pc) + face(ix, iy + 1, iz, pc)
129+ + face(ix, iy, iz - 1, pc) + face(ix, iy, iz + 1, pc);
130+ dst[i] = ${1 / (h * h)}f * wall[i] * s;
131+}
132+`,
133+ };
134+
135+ case 'spread':
136+ return {
137+ args: 1,
138+ outCount: npts,
139+ code: `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
140+@group(0) @binding(1) var<storage, read> src: array<f32>;
141+
142+@compute @workgroup_size(${WORKGROUP_SIZE})
143+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
144+ let i = gid.x;
145+ if (i >= ${npts}u) { return; }
146+ let ix = i % ${nx}u;
147+ let xw = ${-g.Lx / 2}f + (f32(ix) + 0.5) * ${h}f;
148+ let u = (xw - ${g.xs0}f) / ${hs}f;
149+ let j = i32(floor(u));
150+ if (j < 0 || j >= ${ns - 1}) { dst[i] = 0.0; return; }
151+ let f = u - f32(j);
152+ dst[i] = mix(src[u32(j)], src[u32(j) + 1u], f);
153+}
154+`,
155+ };
156+
157+ case 'bridge':
158+ return {
159+ args: 1,
160+ outCount: npts,
161+ code: `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
162+@group(0) @binding(1) var<storage, read> src: array<f32>;
163+
164+@compute @workgroup_size(${WORKGROUP_SIZE})
165+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
166+ let i = gid.x;
167+ if (i >= ${npts}u) { return; }
168+ // du/dx at the last node. The end is pinned, so this is the string's
169+ // arriving slope — what the tension pulls the bridge with.
170+ dst[i] = (src[${ns - 1}u] - src[${ns - 2}u]) * ${1 / hs}f;
171+}
172+`,
173+ };
174+ }
175+}
176+
177+/** Compiled op pipelines for one pair of grids. Shared by every plan. */
178+export class OpPlan {
179+ readonly geometry: OpGeometry;
180+
181+ #device: GPUDevice;
182+ #layouts = new Map<number, GPUBindGroupLayout>();
183+ #pipelines = new Map<OpKind, { pipeline: GPUComputePipeline; outCount: number; args: number }>();
184+
185+ constructor(device: GPUDevice, geometry: OpGeometry) {
186+ this.#device = device;
187+ this.geometry = geometry;
188+ }
189+
190+ /** The op's shape contract, for the planner's checks. */
191+ spec(kind: OpKind): { argCounts: number[]; outCount: number } {
192+ const s = EXTERNAL_OPS.get(kind);
193+ if (!s) throw new UnsupportedOnGpu(`unknown external op '${kind}'`);
194+ const { nx, ny, nz, ns } = this.geometry;
195+ const size = (k: 'air' | 'string'): number => (k === 'air' ? nx * ny * nz : ns);
196+ return { argCounts: s.args.map(size), outCount: size(s.out) };
197+ }
198+
199+ #layout(args: number): GPUBindGroupLayout {
200+ const existing = this.#layouts.get(args);
201+ if (existing) return existing;
202+ const entries: GPUBindGroupLayoutEntry[] = [
203+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
204+ ];
205+ for (let i = 0; i < args; i++) {
206+ entries.push({
207+ binding: i + 1,
208+ visibility: GPUShaderStage.COMPUTE,
209+ buffer: { type: 'read-only-storage' },
210+ });
211+ }
212+ const layout = this.#device.createBindGroupLayout({ label: `op-${args}`, entries });
213+ this.#layouts.set(args, layout);
214+ return layout;
215+ }
216+
217+ /** The pipeline for one op, compiled on first use. */
218+ async pipeline(kind: OpKind): Promise<{ pipeline: GPUComputePipeline; outCount: number; args: number }> {
219+ const existing = this.#pipelines.get(kind);
220+ if (existing) return existing;
221+ const { code, outCount, args } = opWGSL(kind, this.geometry);
222+ const module = this.#device.createShaderModule({ code, label: kind });
223+ let pipeline: GPUComputePipeline;
224+ try {
225+ pipeline = await this.#device.createComputePipelineAsync({
226+ layout: this.#device.createPipelineLayout({ bindGroupLayouts: [this.#layout(args)] }),
227+ compute: { module, entryPoint: 'main' },
228+ label: kind,
229+ });
230+ } catch (e) {
231+ // No error scope here either; see makePipeline in plan.ts.
232+ throw new UnsupportedOnGpu(
233+ `op '${kind}': ${e instanceof Error ? e.message : String(e)}`,
234+ );
235+ }
236+ const built = { pipeline, outCount, args };
237+ this.#pipelines.set(kind, built);
238+ return built;
239+ }
240+
241+ createBinding(kind: OpKind, srcs: GPUBuffer[], dst: GPUBuffer): GPUBindGroup {
242+ const entries: GPUBindGroupEntry[] = [{ binding: 0, resource: { buffer: dst } }];
243+ srcs.forEach((s, i) => entries.push({ binding: i + 1, resource: { buffer: s } }));
244+ return this.#device.createBindGroup({ layout: this.#layout(srcs.length), entries });
245+ }
246+
247+ workgroups(outCount: number): number {
248+ return Math.ceil(outCount / WORKGROUP_SIZE);
249+ }
250+}
src/mgpu/plan.tsadded+761−0View file
@@ -0,0 +1,761 @@
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 { OpPlan, type OpKind } from './ops.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: 'external';
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+ external: OpPlan,
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 planExternal(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 = lapw(p, wall)`: one dispatch, sources and dst distinct. */
464+ async function planExternal(
465+ stmt: Assign,
466+ ext: { name: OpKind; args: (IRExpr & { kind: 'Var' })[] },
467+ dest: Slot,
468+ ): Promise<void> {
469+ const contract = external.spec(ext.name);
470+ const argSlots = ext.args.map((arg, i) => {
471+ const s = slots.get(arg.cName);
472+ if (!s) {
473+ throw new UnsupportedOnGpu(
474+ `'${ext.name}' reads '${arg.name}', which has no buffer`,
475+ stmt.span,
476+ );
477+ }
478+ if (s.count !== contract.argCounts[i]) {
479+ throw new UnsupportedOnGpu(
480+ `'${ext.name}' wants ${contract.argCounts[i]} elements for its ` +
481+ `argument ${i + 1}, but '${arg.name}' holds ${s.count}`,
482+ stmt.span,
483+ );
484+ }
485+ // An op reads its neighbours, so unlike an element-wise kernel it
486+ // cannot be routed through scratch and copied back — the neighbours
487+ // would already have been overwritten. WebGPU forbids the aliasing
488+ // outright anyway; refuse rather than silently reroute.
489+ if (s.buffer === dest.buffer) {
490+ throw new UnsupportedOnGpu(
491+ `'${stmt.name} = ${ext.name}(...)' reads and writes the same ` +
492+ `buffer through '${arg.name}'; assign to a new name instead`,
493+ stmt.span,
494+ );
495+ }
496+ return s;
497+ });
498+ if (dest.count !== contract.outCount) {
499+ throw new UnsupportedOnGpu(
500+ `'${ext.name}' produces a ${contract.outCount}-point field, but ` +
501+ `'${stmt.name}' holds ${dest.count}`,
502+ stmt.span,
503+ );
504+ }
505+ const built = await external.pipeline(ext.name);
506+ ops.push({
507+ kind: 'external',
508+ pipeline: built.pipeline,
509+ bindGroup: external.createBinding(
510+ ext.name,
511+ argSlots.map((s) => s.buffer),
512+ dest.buffer,
513+ ),
514+ workgroups: external.workgroups(contract.outCount),
515+ label: `${stmt.name} = ${ext.name}(${ext.args.map((a) => a.name).join(', ')})`,
516+ });
517+ }
518+
519+ /**
520+ * Feed declared outputs back into the argument buffers they replace, so
521+ * the next call reads what this one produced.
522+ *
523+ * The copies are not independent: a model whose new history field is the
524+ * old current one (`function [pn, pold] = step(p, pm, ...)`) has an output
525+ * whose *source* is another output's *destination*. Doing them in order
526+ * would then copy the new value where the old one was wanted, silently.
527+ * So any source that a previous copy overwrites is staged through scratch
528+ * first — normally none, since a model that writes `pold = p;` gets its
529+ * own buffer from the copy kernel that line plans to.
530+ */
531+ function planFeedback(): void {
532+ const copies: { from: Slot; to: Slot; label: string }[] = [];
533+ fn.outputs.forEach((out, i) => {
534+ const to = spec.feedback[i];
535+ if (!to) return;
536+ const src = slots.get(out.cName);
537+ const dst = host.get(to);
538+ if (!src) {
539+ throw new UnsupportedOnGpu(
540+ `'${fn.name}' declares the output '${out.name}' but never assigns it`,
541+ );
542+ }
543+ if (!dst) throw new UnsupportedOnGpu(`'${to}' is not a host binding`);
544+ if (src.count !== dst.count) {
545+ throw new UnsupportedOnGpu(
546+ `'${out.name}' (${src.count} elements) cannot feed '${to}' (${dst.count})`,
547+ );
548+ }
549+ copies.push({ from: src, to: dst, label: `${out.name} -> ${to}` });
550+ });
551+
552+ const written = new Set<GPUBuffer>();
553+ for (const c of copies) written.add(c.to.buffer);
554+ for (const c of copies) {
555+ // Only a source another copy overwrites needs staging, and only if it
556+ // is not that same copy's own destination (which is a no-op anyway).
557+ if (c.from.buffer !== c.to.buffer && written.has(c.from.buffer)) {
558+ const scratch = alloc(`mgpu-feedback-scratch`, c.from.count);
559+ ops.push({
560+ kind: 'copy',
561+ from: c.from.buffer,
562+ to: scratch.buffer,
563+ bytes: 4 * c.from.count,
564+ label: `${c.label} (staged)`,
565+ });
566+ c.from = scratch;
567+ }
568+ }
569+ for (const c of copies) {
570+ if (c.from.buffer === c.to.buffer) continue; // already in place
571+ ops.push({
572+ kind: 'copy',
573+ from: c.from.buffer,
574+ to: c.to.buffer,
575+ bytes: 4 * c.from.count,
576+ label: c.label,
577+ });
578+ }
579+ }
580+ }
581+
582+ /** Upload parameter values, in `paramNames` order. Cheap — call freely. */
583+ setParams(values: Record<string, number>): void {
584+ this.paramNames.forEach((name, i) => {
585+ const v = values[name];
586+ this.#paramData[i] = Number.isFinite(v) ? v : 0;
587+ });
588+ this.#device.queue.writeBuffer(
589+ this.#paramBuf,
590+ 0,
591+ this.#paramData as Float32Array<ArrayBuffer>,
592+ );
593+ }
594+
595+ /** Buffer holding the named value, or undefined if the .m never binds it. */
596+ buffer(name: string): GPUBuffer | undefined {
597+ return this.#byName.get(name)?.buffer;
598+ }
599+
600+ elementCount(name: string): number | undefined {
601+ return this.#byName.get(name)?.count;
602+ }
603+
604+ /**
605+ * Record `steps` passes of this plan. Synchronous: no awaits, no readback.
606+ * The dispatches share one compute pass, which WebGPU executes in submission
607+ * order with a barrier between them.
608+ *
609+ * `after` runs once per step, inside the same submission — which is what
610+ * lets the microphone sample every timestep rather than every frame.
611+ */
612+ encodeSteps(
613+ encoder: GPUCommandEncoder,
614+ steps: number,
615+ after?: (encoder: GPUCommandEncoder) => void,
616+ ): void {
617+ for (let s = 0; s < steps; s++) {
618+ this.#encodeOps(encoder);
619+ after?.(encoder);
620+ }
621+ }
622+
623+ /** Record one pass over the op sequence into `encoder`. */
624+ #encodeOps(encoder: GPUCommandEncoder): void {
625+ let pass: GPUComputePassEncoder | null = null;
626+ const inPass = (): GPUComputePassEncoder => {
627+ if (!pass) pass = encoder.beginComputePass({ label: 'mgpu-step' });
628+ return pass;
629+ };
630+ const endPass = (): void => {
631+ if (pass) {
632+ pass.end();
633+ pass = null;
634+ }
635+ };
636+ for (const op of this.#ops) {
637+ switch (op.kind) {
638+ case 'kernel': {
639+ const p = inPass();
640+ p.setPipeline(op.pipeline);
641+ p.setBindGroup(0, op.bindGroup);
642+ p.dispatchWorkgroups(Math.ceil(op.count / WORKGROUP_SIZE));
643+ if (op.copyBack) {
644+ endPass();
645+ encoder.copyBufferToBuffer(
646+ op.copyBack.from, 0, op.copyBack.to, 0, op.copyBack.bytes,
647+ );
648+ }
649+ break;
650+ }
651+ case 'external': {
652+ const p = inPass();
653+ p.setPipeline(op.pipeline);
654+ p.setBindGroup(0, op.bindGroup);
655+ p.dispatchWorkgroups(op.workgroups);
656+ break;
657+ }
658+ case 'copy':
659+ endPass();
660+ encoder.copyBufferToBuffer(op.from, 0, op.to, 0, op.bytes);
661+ break;
662+ }
663+ }
664+ endPass();
665+ }
666+
667+ /** Human-readable op sequence — what the .m actually compiled to. */
668+ describe(): string[] {
669+ return this.#ops.map((op) => `${op.kind.padEnd(7)} ${op.label}`);
670+ }
671+
672+ destroy(): void {
673+ for (const b of this.#owned) b.destroy();
674+ this.#paramBuf.destroy();
675+ this.#owned.length = 0;
676+ }
677+}
678+
679+/** `lp = lapw(p, wall)` -> the op's name and its arguments. */
680+function externalCall(
681+ stmt: Assign,
682+): { name: OpKind; args: (IRExpr & { kind: 'Var' })[] } | null {
683+ const e = stmt.expr;
684+ if (e.kind !== 'Call' || !EXTERNAL_OPS.has(e.name)) return null;
685+ const spec = EXTERNAL_OPS.get(e.name)!;
686+ if (e.args.length !== spec.args.length) {
687+ throw new UnsupportedOnGpu(
688+ `'${e.name}' takes ${spec.args.length} argument(s), got ${e.args.length}`,
689+ stmt.span,
690+ );
691+ }
692+ const args = e.args.map((arg) => {
693+ if (arg.kind !== 'Var') {
694+ throw new UnsupportedOnGpu(
695+ `'${e.name}' must be applied to variables, not expressions — ` +
696+ `name the field first`,
697+ stmt.span,
698+ );
699+ }
700+ return arg;
701+ });
702+ return { name: e.name as OpKind, args };
703+}
704+
705+/** Distinct grid fields an expression reads — its storage-buffer cost. */
706+function tensorCount(e: IRExpr): number {
707+ const seen = new Set<string>();
708+ collectTensorVars(e, (c) => seen.add(c));
709+ return seen.size;
710+}
711+
712+/** The subexpressions of a node, in evaluation order. Leaves have none. */
713+function children(e: IRExpr): IRExpr[] {
714+ switch (e.kind) {
715+ case 'Binary':
716+ return [e.left, e.right];
717+ case 'Unary':
718+ return [e.operand];
719+ case 'Call':
720+ return e.args;
721+ default:
722+ return [];
723+ }
724+}
725+
726+/** The same node with its subexpressions replaced. */
727+function withChildren(e: IRExpr, kids: IRExpr[]): IRExpr {
728+ switch (e.kind) {
729+ case 'Binary':
730+ return { ...e, left: kids[0], right: kids[1] };
731+ case 'Unary':
732+ return { ...e, operand: kids[0] };
733+ case 'Call':
734+ return { ...e, args: kids };
735+ default:
736+ return e;
737+ }
738+}
739+
740+function collectTensorVars(e: IRExpr, visit: (cName: string) => void): void {
741+ const walk = (x: IRExpr): void => {
742+ switch (x.kind) {
743+ case 'Var':
744+ if (isTensor(x.ty)) visit(x.cName);
745+ return;
746+ case 'Binary':
747+ walk(x.left);
748+ walk(x.right);
749+ return;
750+ case 'Unary':
751+ walk(x.operand);
752+ return;
753+ case 'Call':
754+ x.args.forEach(walk);
755+ return;
756+ default:
757+ return;
758+ }
759+ };
760+ walk(e);
761+}
src/mgpu/registry.tsadded+143−0View file
@@ -0,0 +1,143 @@
1+/**
2+ * The model: its MATLAB source, and the metadata the host owns.
3+ *
4+ * The model's *algorithm* lives in models/dulcimer.m. Everything around it
5+ * lives here: the parameter names the .m may take as arguments, their
6+ * defaults and slider ranges, the state fields it advances, and which of
7+ * them is the pressure the app draws and the microphone records. The .m
8+ * declares nothing about these — it just names the parameters it wants, and
9+ * `GpuModel` matches each against this table.
10+ */
11+import dulcimerSource from '../../models/dulcimer.m?raw';
12+import type { StateField } from './model.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+ /** Shown as a tooltip. */
25+ hint?: string;
26+}
27+
28+export interface MModel {
29+ key: string;
30+ label: string;
31+ blurb: string;
32+ /** State fields the .m advances, in the order its functions return them. */
33+ state: StateField[];
34+ /** The air pressure field — what gets drawn and recorded. */
35+ pressure: string;
36+ /** The string displacement field — what the string plot shows. */
37+ displacement: string;
38+ params: ParamSpec[];
39+ /** MATLAB source — the algorithm itself. */
40+ source: string;
41+}
42+
43+export const dulcimerModel: MModel = {
44+ key: 'dulcimer',
45+ label: 'Dulcimer',
46+ blurb:
47+ 'A stiff, damped string released from a triangular pluck, driving the ' +
48+ 'acoustic wave equation around a rigid box.',
49+ state: [
50+ { name: 'u', grid: 'string' },
51+ { name: 'um', grid: 'string' },
52+ { name: 'p', grid: 'air' },
53+ { name: 'pm', grid: 'air' },
54+ ],
55+ pressure: 'p',
56+ displacement: 'u',
57+ params: [
58+ {
59+ key: 'f0',
60+ label: 'fundamental (Hz)',
61+ value: 294,
62+ min: 100,
63+ max: 600,
64+ step: 1,
65+ hint: 'The pitch the string is tuned to. 294 Hz is the D above middle C, a common dulcimer melody string. The wave speed on the string follows: cs = 2·Ls·f0.',
66+ },
67+ {
68+ key: 'B',
69+ label: 'inharmonicity',
70+ value: 0.0001,
71+ min: 0,
72+ max: 0.002,
73+ step: 0.00002,
74+ hint: 'Bending stiffness, as the inharmonicity coefficient B: partial n sounds near n·f0·sqrt(1 + B·n²). Zero is an ideal string; ~1e-4 is a light steel string; higher starts to sound bell-like.',
75+ },
76+ {
77+ key: 't60',
78+ label: 'decay t60 (s)',
79+ value: 4,
80+ min: 0.2,
81+ max: 8,
82+ step: 0.1,
83+ hint: 'Seconds for the string to decay 60 dB. What the note’s overall ring is.',
84+ },
85+ {
86+ key: 'sig1',
87+ label: 'brightness decay (m²/s)',
88+ value: 0.005,
89+ min: 0,
90+ max: 0.05,
91+ step: 0.001,
92+ hint: 'Frequency-dependent damping. High partials die faster than low ones, so the note starts bright and mellows; zero keeps it buzzing to the end.',
93+ },
94+ {
95+ key: 'pluckpos',
96+ label: 'pluck position',
97+ value: 0.22,
98+ min: 0.08,
99+ max: 0.92,
100+ step: 0.01,
101+ hint: 'Where along the string it is plucked, as a fraction of its length. Near the middle favours the odd partials (hollower); near the end excites them all (brighter). Takes effect on the next pluck.',
102+ },
103+ {
104+ key: 'amp',
105+ label: 'pluck height (m)',
106+ value: 0.002,
107+ min: 0.0002,
108+ max: 0.005,
109+ step: 0.0002,
110+ hint: 'How far the string is pulled before release. The equations are linear, so this scales everything and changes nothing else. Takes effect on the next pluck.',
111+ },
112+ {
113+ key: 'gline',
114+ label: 'string radiation',
115+ value: 0.3,
116+ min: 0,
117+ max: 2,
118+ step: 0.05,
119+ hint: 'Gain on the string radiating directly along its length. A thin string barely does this in reality, so it is the idealized route; 0 switches it off.',
120+ },
121+ {
122+ key: 'gbridge',
123+ label: 'bridge drive',
124+ value: 1,
125+ min: 0,
126+ max: 2,
127+ step: 0.05,
128+ hint: 'Gain on the string’s pull at the bridge driving the top-plate patch — the real instrument’s main route into the air. 0 switches it off.',
129+ },
130+ ],
131+ source: dulcimerSource,
132+};
133+
134+/** The slider maxima the string grid must stay stable for, whatever the user
135+ * drags to (see makeStringGrid). */
136+export function worstCase(model: MModel): { f0: number; B: number; sig1: number } {
137+ const max = (key: string): number =>
138+ model.params.find((p) => p.key === key)?.max ?? 0;
139+ return { f0: max('f0'), B: max('B'), sig1: max('sig1') };
140+}
141+
142+export const defaultParams = (m: MModel): Params =>
143+ Object.fromEntries(m.params.map((p) => [p.key, p.value]));
src/mgpu/session.tsadded+252−0View file
@@ -0,0 +1,252 @@
1+/**
2+ * One running simulation: the two grids, the body, the compiled .m, the
3+ * timestep, the microphone.
4+ *
5+ * Everything that is not rendering. The app and the tests both go through
6+ * this, so there is one place that decides how a pair of .m files becomes
7+ * something running on the GPU — and nothing about it is browser-specific
8+ * beyond needing a GPUDevice.
9+ *
10+ * The order of construction matters and is worth spelling out. The air grid
11+ * fixes h; the scene fixes the fastest speed; together they fix dt (the CFL
12+ * condition). Only then can the string grid be sized, because the string has
13+ * no say in the timestep — the air dictates dt, and the string chooses the
14+ * finest spacing that is stable at that dt for anything the sliders can ask
15+ * (see makeStringGrid). Then the ops and the model compile against both
16+ * grids.
17+ */
18+import {
19+ makeAirGrid,
20+ makeStringGrid,
21+ stableDt,
22+ type AirGrid,
23+ type StringGrid,
24+} from '../grid.ts';
25+import { OpPlan } from './ops.ts';
26+import { GpuModel, type ModelParams } from './model.ts';
27+import { Scene } from '../scene/scene.ts';
28+import { Recorder } from '../audio/recorder.ts';
29+import { worstCase, type MModel, type Params } from './registry.ts';
30+import type { MScene } from '../scene/registry.ts';
31+import { C_AIR, CFL, DOMAIN_X } from '../units.ts';
32+
33+export interface ModelSessionOptions {
34+ device: GPUDevice;
35+ model: MModel;
36+ params: Params;
37+ /** Override the model source — the editor's working copy. */
38+ source?: string;
39+ scene: MScene;
40+ sceneParams: Params;
41+ /** Override the scene source — the editor's working copy. */
42+ sceneSource?: string;
43+ /** Air grid points along x; y and z get half each. */
44+ nx: number;
45+ /** Domain length along x, metres. */
46+ Lx?: number;
47+ /** String length, metres. */
48+ Ls: number;
49+ /** Grid fields one kernel may read, overriding what the device allows.
50+ * Only for tests, which use it to exercise the planner's kernel splitting
51+ * on a device that would never need it. */
52+ operandBudget?: number;
53+}
54+
55+export class ModelSession {
56+ readonly device: GPUDevice;
57+ readonly model: MModel;
58+ readonly air: AirGrid;
59+ readonly string: StringGrid;
60+ readonly gpu: GpuModel;
61+ /** The microphone: the pressure at one point, every timestep. */
62+ readonly recorder: Recorder;
63+
64+ /** Model time and step count since the last pluck. */
65+ t = 0;
66+ steps = 0;
67+
68+ #scene: Scene;
69+ #sceneModel: MScene;
70+ #params: Params;
71+ #dt: number;
72+
73+ private constructor(init: {
74+ device: GPUDevice;
75+ model: MModel;
76+ air: AirGrid;
77+ string: StringGrid;
78+ gpu: GpuModel;
79+ recorder: Recorder;
80+ scene: Scene;
81+ sceneModel: MScene;
82+ params: Params;
83+ dt: number;
84+ }) {
85+ this.device = init.device;
86+ this.model = init.model;
87+ this.air = init.air;
88+ this.string = init.string;
89+ this.gpu = init.gpu;
90+ this.recorder = init.recorder;
91+ this.#scene = init.scene;
92+ this.#sceneModel = init.sceneModel;
93+ this.#params = init.params;
94+ this.#dt = init.dt;
95+ }
96+
97+ static async create(opts: ModelSessionOptions): Promise<ModelSession> {
98+ const { device, model, params, scene, sceneParams } = opts;
99+ const air = makeAirGrid(opts.nx, opts.Lx ?? DOMAIN_X);
100+
101+ // The scene first: the timestep depends on the fastest speed in it.
102+ const built = Scene.create({
103+ air,
104+ Ls: opts.Ls,
105+ source: opts.sceneSource ?? scene.source,
106+ paramNames: scene.params.map((p) => p.key),
107+ params: sceneParams,
108+ });
109+
110+ const dt = stableDt(air.h, Math.max(built.cmax, C_AIR), CFL);
111+ const string = makeStringGrid(opts.Ls, dt, worstCase(model));
112+
113+ const ops = new OpPlan(device, {
114+ nx: air.nx,
115+ ny: air.ny,
116+ nz: air.nz,
117+ h: air.h,
118+ ns: string.ns,
119+ hs: string.hs,
120+ xs0: -string.Ls / 2,
121+ Lx: air.Lx,
122+ });
123+
124+ const gpu = await GpuModel.create({
125+ device,
126+ ops,
127+ air,
128+ string,
129+ medium: built,
130+ source: opts.source ?? model.source,
131+ paramNames: model.params.map((p) => p.key),
132+ state: model.state,
133+ operandBudget: opts.operandBudget,
134+ });
135+ gpu.setDt(dt);
136+
137+ // The microphone listens to the host-owned pressure buffer, which is
138+ // where both `init` and `step` leave their result.
139+ const recorder = await Recorder.create({
140+ device,
141+ field: gpu.stateBuffer(model.pressure)!,
142+ nx: air.nx,
143+ ny: air.ny,
144+ nz: air.nz,
145+ });
146+
147+ const session = new ModelSession({
148+ device, model, air, string, gpu, recorder,
149+ scene: built, sceneModel: scene, params, dt,
150+ });
151+ gpu.setParams(params);
152+ return session;
153+ }
154+
155+ get scene(): Scene {
156+ return this.#scene;
157+ }
158+
159+ get sceneModel(): MScene {
160+ return this.#sceneModel;
161+ }
162+
163+ /** The air pressure field's name — what gets drawn and recorded. */
164+ get pressureName(): string {
165+ return this.model.pressure;
166+ }
167+
168+ /** The string displacement field's name — what the string plot shows. */
169+ get displacementName(): string {
170+ return this.model.displacement;
171+ }
172+
173+ get dt(): number {
174+ return this.#dt;
175+ }
176+
177+ /**
178+ * Swap the body under a running model. It is data, not code, so this needs
179+ * no recompile. The timestep is deliberately left alone: it was set from
180+ * max(cmax, c_air) at build time, and a scene edit that *raises* the
181+ * fastest speed above that needs a rebuild anyway (the caller compares
182+ * `dtWanted` and rebuilds when they disagree).
183+ */
184+ setScene(sceneModel: MScene, params: Params, source?: string): void {
185+ const built = Scene.create({
186+ air: this.air,
187+ Ls: this.string.Ls,
188+ source: source ?? sceneModel.source,
189+ paramNames: sceneModel.params.map((p) => p.key),
190+ params,
191+ });
192+ this.#scene = built;
193+ this.#sceneModel = sceneModel;
194+ this.gpu.uploadMedium(built);
195+ }
196+
197+ /** The timestep the current scene would ask for. Differs from `dt` only
198+ * when a scene edit changed the fastest speed, which calls for a rebuild. */
199+ get dtWanted(): number {
200+ return stableDt(this.air.h, Math.max(this.#scene.cmax, C_AIR), CFL);
201+ }
202+
203+ setParams(params: ModelParams): void {
204+ this.#params = params;
205+ this.gpu.setParams(params);
206+ }
207+
208+ /** Run `init`: the string drawn into its pluck, silent air, t = 0. */
209+ pluck(): void {
210+ this.gpu.init();
211+ this.recorder.clear();
212+ this.t = 0;
213+ this.steps = 0;
214+ }
215+
216+ /** Put the microphone at the grid point nearest (x, y, z), metres. */
217+ setMic(x: number, y: number, z: number): void {
218+ const { Lx, Ly, Lz, h } = this.air;
219+ this.recorder.setProbe(
220+ (x + Lx / 2) / h - 0.5,
221+ (y + Ly / 2) / h - 0.5,
222+ (z + Lz / 2) / h - 0.5,
223+ );
224+ }
225+
226+ /** Advance `n` steps. Synchronous: records and submits, nothing read back.
227+ * The microphone samples inside the same submission, once per step. */
228+ step(n = 1): void {
229+ this.gpu.step(n, (enc) => this.recorder.encode(enc));
230+ this.t += n * this.#dt;
231+ this.steps += n;
232+ }
233+
234+ /** Wait for the submitted steps to finish, without reading anything back. */
235+ sync(): Promise<undefined> {
236+ return this.device.queue.onSubmittedWorkDone();
237+ }
238+
239+ /** Read a named field back to the CPU. */
240+ read(name: string): Promise<Float32Array> {
241+ return this.gpu.read(name);
242+ }
243+
244+ describe(): { init: string[]; step: string[] } {
245+ return this.gpu.describe();
246+ }
247+
248+ destroy(): void {
249+ this.recorder.destroy();
250+ this.gpu.destroy();
251+ }
252+}
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/overlay.tsadded+154−0View file
@@ -0,0 +1,154 @@
1+/**
2+ * The line work over the volume: the string, the body's wireframe, the sound
3+ * hole, the bridge and nut.
4+ *
5+ * A 2D canvas stacked on the WebGPU one, redrawn each frame from the string
6+ * readback and the same camera the ray march uses (cameraFrame), so its lines
7+ * land on the volume's pixels. Lines composited over a volume have no correct
8+ * depth — the overlay is always on top — which is the price of keeping the
9+ * march simple; these are markers and guides, not things in the scene.
10+ *
11+ * The string's actual displacement is millimetres on a half-metre picture, so
12+ * it is drawn exaggerated (the same factor the string plot states), and it is
13+ * the one honest exaggeration in the app: everything else is to scale.
14+ */
15+import type { CameraFrame } from './volume.ts';
16+
17+export interface OverlayGeometry {
18+ /** String extent (centred on x) and height above the top plate. */
19+ Ls: number;
20+ stringZ: number;
21+ /** Body box: length, width, depth below z = 0. */
22+ boxl: number;
23+ boxw: number;
24+ boxd: number;
25+ /** Sound hole. */
26+ holer: number;
27+ holex: number;
28+}
29+
30+export interface OverlayOptions {
31+ frame: CameraFrame;
32+ geometry: OverlayGeometry;
33+ /** String displacement, metres, at the ns nodes; null hides the string. */
34+ string: Float32Array | null;
35+ /** Displacement multiplier for visibility. */
36+ exaggerate: number;
37+ wireframe: boolean;
38+}
39+
40+type P3 = [number, number, number];
41+
42+export class Overlay {
43+ readonly canvas: HTMLCanvasElement;
44+ #ctx: CanvasRenderingContext2D;
45+
46+ constructor(canvas: HTMLCanvasElement) {
47+ this.canvas = canvas;
48+ const ctx = canvas.getContext('2d');
49+ if (!ctx) throw new Error('overlay canvas has no 2d context');
50+ this.#ctx = ctx;
51+ }
52+
53+ resize(): void {
54+ const rect = this.canvas.getBoundingClientRect();
55+ const w = Math.max(1, Math.round(rect.width));
56+ const h = Math.max(1, Math.round(rect.height));
57+ if (this.canvas.width !== w || this.canvas.height !== h) {
58+ this.canvas.width = w;
59+ this.canvas.height = h;
60+ }
61+ }
62+
63+ clear(): void {
64+ this.#ctx.clearRect(0, 0, this.canvas.width, this.canvas.height);
65+ }
66+
67+ draw(o: OverlayOptions): void {
68+ const ctx = this.#ctx;
69+ const { width, height } = this.canvas;
70+ ctx.clearRect(0, 0, width, height);
71+
72+ const f = o.frame;
73+ /** World point -> pixel, or null when behind the camera. */
74+ const project = (q: P3): [number, number] | null => {
75+ const dx = q[0] - f.eye[0];
76+ const dy = q[1] - f.eye[1];
77+ const dz = q[2] - f.eye[2];
78+ const zc = dx * f.fwd[0] + dy * f.fwd[1] + dz * f.fwd[2];
79+ if (zc < 1e-6) return null;
80+ const xc = dx * f.right[0] + dy * f.right[1] + dz * f.right[2];
81+ const yc = dx * f.up[0] + dy * f.up[1] + dz * f.up[2];
82+ const nx = xc / (zc * f.tanHalfFov * f.aspect);
83+ const ny = yc / (zc * f.tanHalfFov);
84+ return [((nx + 1) / 2) * width, ((1 - ny) / 2) * height];
85+ };
86+
87+ const polyline = (pts: P3[], close = false): void => {
88+ ctx.beginPath();
89+ let started = false;
90+ for (const q of pts) {
91+ const s = project(q);
92+ if (!s) {
93+ started = false;
94+ continue;
95+ }
96+ if (started) ctx.lineTo(s[0], s[1]);
97+ else {
98+ ctx.moveTo(s[0], s[1]);
99+ started = true;
100+ }
101+ }
102+ if (close) ctx.closePath();
103+ ctx.stroke();
104+ };
105+
106+ const g = o.geometry;
107+
108+ if (o.wireframe) {
109+ // The body: a box from z = -boxd to 0.
110+ const hx = g.boxl / 2;
111+ const hy = g.boxw / 2;
112+ ctx.lineWidth = 1;
113+ ctx.strokeStyle = 'rgba(160, 170, 185, 0.65)';
114+ const corners = (zv: number): P3[] => [
115+ [-hx, -hy, zv], [hx, -hy, zv], [hx, hy, zv], [-hx, hy, zv],
116+ ];
117+ polyline(corners(0), true);
118+ polyline(corners(-g.boxd), true);
119+ for (let i = 0; i < 4; i++) {
120+ const top = corners(0)[i];
121+ const bot = corners(-g.boxd)[i];
122+ polyline([top, bot]);
123+ }
124+ // The sound hole, on the top plate.
125+ if (g.holer > 0) {
126+ const circle: P3[] = [];
127+ for (let i = 0; i <= 48; i++) {
128+ const a = (2 * Math.PI * i) / 48;
129+ circle.push([g.holex + g.holer * Math.cos(a), g.holer * Math.sin(a), 0]);
130+ }
131+ polyline(circle);
132+ }
133+ }
134+
135+ if (o.string) {
136+ const u = o.string;
137+ const ns = u.length;
138+ const hs = g.Ls / (ns - 1);
139+ // Displacement is in y (the pluck direction), exaggerated to be seen.
140+ const pts: P3[] = [];
141+ for (let i = 0; i < ns; i++) {
142+ pts.push([-g.Ls / 2 + i * hs, o.exaggerate * u[i], g.stringZ]);
143+ }
144+ ctx.lineWidth = 1.6;
145+ ctx.strokeStyle = 'rgba(255, 214, 130, 0.95)';
146+ polyline(pts);
147+ // Nut and bridge: short posts down to the plate.
148+ ctx.lineWidth = 2;
149+ ctx.strokeStyle = 'rgba(220, 226, 235, 0.8)';
150+ polyline([[-g.Ls / 2, 0, g.stringZ], [-g.Ls / 2, 0, 0]]);
151+ polyline([[g.Ls / 2, 0, g.stringZ], [g.Ls / 2, 0, 0]]);
152+ }
153+ }
154+}
src/render/stringplot.tsadded+67−0View file
@@ -0,0 +1,67 @@
1+/**
2+ * The string on its own: u(x), drawn flat.
3+ *
4+ * The 3D view shows the string in context but a millimetre of displacement
5+ * needs exaggeration there; this plot gives the displacement an honest axis.
6+ * The vertical scale is fixed to the pluck height rather than following the
7+ * field, so the decay of the note is visible as the curve settling, not as an
8+ * axis chasing it.
9+ */
10+export class StringPlot {
11+ #canvas: HTMLCanvasElement;
12+ #ctx: CanvasRenderingContext2D;
13+
14+ constructor(canvas: HTMLCanvasElement) {
15+ this.#canvas = canvas;
16+ const ctx = canvas.getContext('2d');
17+ if (!ctx) throw new Error('string plot canvas has no 2d context');
18+ this.#ctx = ctx;
19+ }
20+
21+ resize(): void {
22+ const rect = this.#canvas.getBoundingClientRect();
23+ const w = Math.max(1, Math.round(rect.width * devicePixelRatio));
24+ const h = Math.max(1, Math.round(rect.height * devicePixelRatio));
25+ if (this.#canvas.width !== w || this.#canvas.height !== h) {
26+ this.#canvas.width = w;
27+ this.#canvas.height = h;
28+ }
29+ }
30+
31+ /** Draw the displacement `u` (metres) against a fixed range ±`umax`. */
32+ draw(u: Float32Array | null, umax: number, dark: boolean): void {
33+ const ctx = this.#ctx;
34+ const { width: w, height: h } = this.#canvas;
35+ ctx.clearRect(0, 0, w, h);
36+
37+ const ink = dark ? 'rgba(230, 233, 236, 0.9)' : 'rgba(31, 35, 40, 0.9)';
38+ const faint = dark ? 'rgba(154, 164, 175, 0.35)' : 'rgba(87, 96, 106, 0.35)';
39+ const pad = 6 * devicePixelRatio;
40+
41+ // Rest line and end posts.
42+ ctx.lineWidth = devicePixelRatio;
43+ ctx.strokeStyle = faint;
44+ ctx.beginPath();
45+ ctx.moveTo(pad, h / 2);
46+ ctx.lineTo(w - pad, h / 2);
47+ ctx.stroke();
48+ ctx.beginPath();
49+ ctx.moveTo(pad, pad);
50+ ctx.lineTo(pad, h - pad);
51+ ctx.moveTo(w - pad, pad);
52+ ctx.lineTo(w - pad, h - pad);
53+ ctx.stroke();
54+
55+ if (!u || u.length < 2 || !(umax > 0)) return;
56+ ctx.lineWidth = 1.6 * devicePixelRatio;
57+ ctx.strokeStyle = ink;
58+ ctx.beginPath();
59+ for (let i = 0; i < u.length; i++) {
60+ const x = pad + ((w - 2 * pad) * i) / (u.length - 1);
61+ const y = h / 2 - (h / 2 - pad) * Math.max(-1, Math.min(1, u[i] / umax));
62+ if (i === 0) ctx.moveTo(x, y);
63+ else ctx.lineTo(x, y);
64+ }
65+ ctx.stroke();
66+ }
67+}
src/render/volume.tsadded+379−0View file
@@ -0,0 +1,379 @@
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-to-CPU round trip. Adapted from acoustic-scattering-3d's renderer; the
7+ * differences are a rectangular domain instead of a cube, and the body drawn
8+ * from the wall mask instead of from a speed contrast.
9+ *
10+ * The picture is a ray march. Each pixel casts one ray, intersects it with
11+ * the domain box, and steps along it accumulating emission front to back:
12+ * the pressure goes through a diverging colormap about zero, and the opacity
13+ * goes as a power of |p|, so quiet regions are transparent and the
14+ * wavefronts are what you see. The body is added as a grey emission where
15+ * the wall mask says solid, which shows the box and its sound hole without
16+ * needing its own kind of drawing.
17+ *
18+ * Two honest limitations. Sampling is nearest-neighbour, not trilinear: the
19+ * field lives in a storage buffer rather than a filterable 3D texture, so
20+ * trilinear would be eight fetches per sample and the march takes tens of
21+ * millions of samples a frame. With the ray step set near the cell size the
22+ * difference is visible mainly as a faint stippling on strong wavefronts.
23+ * And the compositing is emission only, with no lighting and no shadowing,
24+ * so what is behind a strong feature is dimmed but never occluded correctly.
25+ *
26+ * A clip plane on y is provided because a volume render of a wavefield is
27+ * mostly the outside of a wavefield. Pulling the clip in cuts the picture
28+ * lengthwise through the string and the cavity, which is the view that
29+ * actually shows the instrument working.
30+ */
31+import type { ColormapFunc } from './colormaps.ts';
32+
33+const SHADER = `
34+struct View {
35+ eye: vec4f, // xyz: eye position, w: tan(fov/2)
36+ right: vec4f, // xyz: camera right, w: aspect ratio
37+ up: vec4f, // xyz: camera up, w: pressure the colormap saturates at
38+ fwd: vec4f, // xyz: camera forward, w: contrast exponent
39+ bg: vec4f, // rgb: background, w: body strength
40+ dims: vec4f, // Lx, Ly, Lz, h
41+ grid: vec4f, // nx, ny, nz, ray steps
42+ misc: vec4f, // opacity, clip y, (unused), (unused)
43+ m: vec4f, // xyz: microphone position, w: whether to draw it
44+};
45+
46+@group(0) @binding(0) var<uniform> V: View;
47+@group(0) @binding(1) var<storage, read> p: array<f32>;
48+@group(0) @binding(2) var<storage, read> wall: array<f32>;
49+@group(0) @binding(3) var cmap: texture_2d<f32>;
50+@group(0) @binding(4) var samp: sampler;
51+
52+struct VSOut {
53+ @builtin(position) pos: vec4f,
54+ @location(0) ndc: vec2f,
55+};
56+
57+@vertex
58+fn vs(@builtin(vertex_index) vi: u32) -> VSOut {
59+ // One oversized triangle covering the viewport.
60+ var xy = array<vec2f, 3>(vec2f(-1.0, -3.0), vec2f(-1.0, 1.0), vec2f(3.0, 1.0));
61+ var out: VSOut;
62+ let q = xy[vi];
63+ out.pos = vec4f(q, 0.0, 1.0);
64+ out.ndc = q;
65+ return out;
66+}
67+
68+fn voxel(q: vec3f) -> u32 {
69+ let n = vec3i(i32(V.grid.x), i32(V.grid.y), i32(V.grid.z));
70+ let h = V.dims.w;
71+ let g = clamp(vec3i(floor((q + 0.5 * V.dims.xyz) / h)), vec3i(0), n - vec3i(1));
72+ return u32(g.x + n.x * (g.y + n.y * g.z));
73+}
74+
75+// A point on the box's surface is on an edge when two of its three distances
76+// to the bounding planes vanish, so the test is on the median of the three.
77+fn edge(q: vec3f) -> f32 {
78+ let d = abs(abs(q) - 0.5 * V.dims.xyz);
79+ let lo = min(d.x, min(d.y, d.z));
80+ let hi = max(d.x, max(d.y, d.z));
81+ let mid = d.x + d.y + d.z - lo - hi;
82+ let w = 0.004 * V.dims.x;
83+ return 1.0 - smoothstep(w, 2.0 * w, mid);
84+}
85+
86+@fragment
87+fn fs(in: VSOut) -> @location(0) vec4f {
88+ let half = 0.5 * V.dims.xyz;
89+ let eye = V.eye.xyz;
90+ let dir = normalize(V.fwd.xyz
91+ + in.ndc.x * V.right.w * V.eye.w * V.right.xyz
92+ + in.ndc.y * V.eye.w * V.up.xyz);
93+
94+ let inv = 1.0 / dir;
95+ let ta = (-half - eye) * inv;
96+ let tb = (half - eye) * inv;
97+ let lo = min(ta, tb);
98+ let hi = max(ta, tb);
99+ let t1 = min(min(hi.x, hi.y), hi.z);
100+ var t0 = max(max(lo.x, lo.y), lo.z);
101+
102+ var col = V.bg.rgb;
103+ if (t1 <= max(t0, 0.0)) { return vec4f(col, 1.0); }
104+ t0 = max(t0, 0.0);
105+
106+ let lineCol = vec3f(0.42, 0.47, 0.55);
107+ col = mix(col, lineCol, 0.5 * edge(eye + dir * t1));
108+
109+ let steps = i32(V.grid.w);
110+ let dl = (t1 - t0) / f32(steps);
111+ // Opacity is quoted per cell, so a longer ray step is proportionally more
112+ // opaque and the picture does not change brightness with the quality knob.
113+ let unit = dl / V.dims.w;
114+ let scale = max(V.up.w, 1e-20);
115+ let grey = vec3f(0.55, 0.58, 0.63);
116+
117+ var acc = vec3f(0.0);
118+ var alpha = 0.0;
119+ for (var k = 0; k < steps; k = k + 1) {
120+ if (alpha > 0.995) { break; }
121+ let q = eye + dir * (t0 + (f32(k) + 0.5) * dl);
122+ if (q.y > V.misc.y) { continue; }
123+ let i = voxel(q);
124+ let v = clamp(p[i] / scale, -1.0, 1.0);
125+ var a = pow(abs(v), V.fwd.w) * V.misc.x * unit;
126+ var rgb = textureSampleLevel(cmap, samp, vec2f(0.5 + 0.5 * v, 0.5), 0.0).rgb;
127+ if (V.bg.w > 0.0) {
128+ let m = clamp(1.0 - wall[i], 0.0, 1.0);
129+ let am = m * V.bg.w * unit;
130+ let tot = a + am;
131+ if (tot > 1e-12) { rgb = (a * rgb + am * grey) / tot; }
132+ a = tot;
133+ }
134+ a = clamp(a, 0.0, 1.0);
135+ acc = acc + (1.0 - alpha) * a * rgb;
136+ alpha = alpha + (1.0 - alpha) * a;
137+ }
138+ col = acc + (1.0 - alpha) * col;
139+ col = mix(col, lineCol, 0.75 * edge(eye + dir * t0));
140+
141+ // The microphone, a dot in a ring at the ray's closest approach to it.
142+ // Drawn on top rather than composited into the march: a marker, not a
143+ // thing in the scene.
144+ if (V.m.w > 0.0) {
145+ let toM = V.m.xyz - eye;
146+ let tm = dot(toM, dir);
147+ if (tm > 0.0) {
148+ let d = length(toM - tm * dir) / V.dims.x;
149+ let dotm = 1.0 - smoothstep(0.004, 0.007, d);
150+ let ring = 1.0 - smoothstep(0.0015, 0.004, abs(d - 0.016));
151+ col = mix(col, vec3f(1.0), 0.9 * max(dotm, ring));
152+ }
153+ }
154+ return vec4f(col, 1.0);
155+}
156+`;
157+
158+/** Where the camera is looking from. Angles in radians, distance in metres. */
159+export interface Camera {
160+ az: number;
161+ el: number;
162+ dist: number;
163+}
164+
165+/** The camera's orthonormal frame, shared with the overlay so its lines land
166+ * on the volume's pixels. */
167+export interface CameraFrame {
168+ eye: [number, number, number];
169+ right: [number, number, number];
170+ up: [number, number, number];
171+ fwd: [number, number, number];
172+ tanHalfFov: number;
173+ aspect: number;
174+}
175+
176+export const TAN_HALF_FOV = Math.tan((32 * Math.PI) / 360);
177+
178+const cross = (a: number[], b: number[]): [number, number, number] => [
179+ a[1] * b[2] - a[2] * b[1],
180+ a[2] * b[0] - a[0] * b[2],
181+ a[0] * b[1] - a[1] * b[0],
182+];
183+const norm = (a: [number, number, number]): [number, number, number] => {
184+ const m = Math.hypot(a[0], a[1], a[2]) || 1;
185+ return [a[0] / m, a[1] / m, a[2] / m];
186+};
187+
188+export function cameraFrame(camera: Camera, aspect: number): CameraFrame {
189+ const { az, el, dist } = camera;
190+ const ce = Math.cos(el);
191+ const toEye: [number, number, number] = [
192+ ce * Math.cos(az),
193+ ce * Math.sin(az),
194+ Math.sin(el),
195+ ];
196+ const eye: [number, number, number] = [toEye[0] * dist, toEye[1] * dist, toEye[2] * dist];
197+ const fwd: [number, number, number] = [-toEye[0], -toEye[1], -toEye[2]];
198+ const right = norm(cross([0, 0, 1], fwd));
199+ const up = cross(fwd, right);
200+ return { eye, right, up, fwd, tanHalfFov: TAN_HALF_FOV, aspect };
201+}
202+
203+export interface DrawOptions {
204+ frame: CameraFrame;
205+ /** Pressure the colormap saturates at, in both directions. */
206+ scale: number;
207+ opacity: number;
208+ contrast: number;
209+ /** Samples along each ray. */
210+ steps: number;
211+ /** Everything with y above this is not drawn, in metres. */
212+ clipY: number;
213+ /** Strength of the body's grey emission, 0 to turn it off. */
214+ body: number;
215+ /** Microphone position in metres, or null to not draw it. */
216+ mic?: { x: number; y: number; z: number } | null;
217+}
218+
219+export interface VolumeDims {
220+ nx: number;
221+ ny: number;
222+ nz: number;
223+ Lx: number;
224+ Ly: number;
225+ Lz: number;
226+ h: number;
227+}
228+
229+export class VolumeView {
230+ readonly canvas: HTMLCanvasElement;
231+
232+ #device: GPUDevice;
233+ #context: GPUCanvasContext;
234+ #pipeline: GPURenderPipeline;
235+ #layout: GPUBindGroupLayout;
236+ #uniform: GPUBuffer;
237+ #host = new ArrayBuffer(9 * 16);
238+ #sampler: GPUSampler;
239+ #cmapTexture: GPUTexture;
240+ #bindGroup: GPUBindGroup | null = null;
241+ #dims: VolumeDims | null = null;
242+ /** Background, matched to the CSS so the box sits on the page. */
243+ bg: [number, number, number] = [0.043, 0.055, 0.071];
244+
245+ constructor(device: GPUDevice, canvas: HTMLCanvasElement) {
246+ this.#device = device;
247+ this.canvas = canvas;
248+
249+ const context = canvas.getContext('webgpu');
250+ if (!context) throw new Error('this canvas has no WebGPU context');
251+ this.#context = context;
252+ const format = navigator.gpu.getPreferredCanvasFormat();
253+ context.configure({ device, format, alphaMode: 'opaque' });
254+
255+ this.#layout = device.createBindGroupLayout({
256+ label: 'volume-view',
257+ entries: [
258+ { binding: 0, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'uniform' } },
259+ { binding: 1, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'read-only-storage' } },
260+ { binding: 2, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'read-only-storage' } },
261+ { binding: 3, visibility: GPUShaderStage.FRAGMENT, texture: { sampleType: 'float' } },
262+ { binding: 4, visibility: GPUShaderStage.FRAGMENT, sampler: { type: 'filtering' } },
263+ ],
264+ });
265+
266+ const module = device.createShaderModule({ code: SHADER, label: 'volume-view' });
267+ this.#pipeline = device.createRenderPipeline({
268+ label: 'volume-view',
269+ layout: device.createPipelineLayout({ bindGroupLayouts: [this.#layout] }),
270+ vertex: { module, entryPoint: 'vs' },
271+ fragment: { module, entryPoint: 'fs', targets: [{ format }] },
272+ primitive: { topology: 'triangle-list' },
273+ });
274+
275+ this.#uniform = device.createBuffer({
276+ label: 'volume-view-uniform',
277+ size: this.#host.byteLength,
278+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
279+ });
280+ this.#sampler = device.createSampler({ magFilter: 'linear', minFilter: 'linear' });
281+ this.#cmapTexture = device.createTexture({
282+ label: 'colormap',
283+ size: [256, 1],
284+ format: 'rgba8unorm',
285+ usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST,
286+ });
287+ }
288+
289+ /** Point the view at the buffers of a (new) simulation: the host-owned
290+ * pressure state and the wall mask, which is what the body is drawn from. */
291+ setSource(pressure: GPUBuffer, wall: GPUBuffer, dims: VolumeDims): void {
292+ this.#dims = dims;
293+ this.#bindGroup = this.#device.createBindGroup({
294+ layout: this.#layout,
295+ entries: [
296+ { binding: 0, resource: { buffer: this.#uniform } },
297+ { binding: 1, resource: { buffer: pressure } },
298+ { binding: 2, resource: { buffer: wall } },
299+ { binding: 3, resource: this.#cmapTexture.createView() },
300+ { binding: 4, resource: this.#sampler },
301+ ],
302+ });
303+ }
304+
305+ setColormap(cmap: ColormapFunc): void {
306+ const data = new Uint8Array(256 * 4);
307+ for (let i = 0; i < 256; i++) {
308+ const [r, g, b] = cmap(i / 255);
309+ data[4 * i] = r;
310+ data[4 * i + 1] = g;
311+ data[4 * i + 2] = b;
312+ data[4 * i + 3] = 255;
313+ }
314+ this.#device.queue.writeTexture(
315+ { texture: this.#cmapTexture },
316+ data,
317+ { bytesPerRow: 256 * 4 },
318+ { width: 256, height: 1 },
319+ );
320+ }
321+
322+ /** Draw one frame. */
323+ draw(o: DrawOptions): void {
324+ const bg = this.#bindGroup;
325+ const dims = this.#dims;
326+ if (!bg || !dims) return;
327+
328+ const { eye, right, up, fwd } = o.frame;
329+ const f = new Float32Array(this.#host);
330+ f.set([eye[0], eye[1], eye[2], o.frame.tanHalfFov], 0);
331+ f.set([right[0], right[1], right[2], o.frame.aspect], 4);
332+ f.set([up[0], up[1], up[2], o.scale], 8);
333+ f.set([fwd[0], fwd[1], fwd[2], o.contrast], 12);
334+ f.set([this.bg[0], this.bg[1], this.bg[2], o.body], 16);
335+ f.set([dims.Lx, dims.Ly, dims.Lz, dims.h], 20);
336+ f.set([dims.nx, dims.ny, dims.nz, o.steps], 24);
337+ f.set([o.opacity, o.clipY, 0, 0], 28);
338+ f.set([o.mic?.x ?? 0, o.mic?.y ?? 0, o.mic?.z ?? 0, o.mic ? 1 : 0], 32);
339+ this.#device.queue.writeBuffer(this.#uniform, 0, this.#host);
340+
341+ const enc = this.#device.createCommandEncoder({ label: 'volume-view' });
342+ const pass = enc.beginRenderPass({
343+ colorAttachments: [
344+ {
345+ view: this.#context.getCurrentTexture().createView(),
346+ clearValue: { r: this.bg[0], g: this.bg[1], b: this.bg[2], a: 1 },
347+ loadOp: 'clear',
348+ storeOp: 'store',
349+ },
350+ ],
351+ });
352+ pass.setPipeline(this.#pipeline);
353+ pass.setBindGroup(0, bg);
354+ pass.draw(3);
355+ pass.end();
356+ this.#device.queue.submit([enc.finish()]);
357+ }
358+
359+ /**
360+ * Match the canvas's backing store to its CSS size. Device pixel ratio is
361+ * ignored: the march is the whole cost of a frame and it scales with
362+ * pixels, so a retina display would pay four times over for a picture that
363+ * is already smooth.
364+ */
365+ resize(): void {
366+ const rect = this.canvas.getBoundingClientRect();
367+ const w = Math.max(1, Math.round(rect.width));
368+ const h = Math.max(1, Math.round(rect.height));
369+ if (this.canvas.width !== w || this.canvas.height !== h) {
370+ this.canvas.width = w;
371+ this.canvas.height = h;
372+ }
373+ }
374+
375+ destroy(): void {
376+ this.#uniform.destroy();
377+ this.#cmapTexture.destroy();
378+ }
379+}
src/scene/registry.tsadded+120−0View file
@@ -0,0 +1,120 @@
1+/**
2+ * The scene: its MATLAB source, and the parameters the host offers it.
3+ *
4+ * Same split as the model (src/mgpu/registry.ts): what the body *is* lives in
5+ * the .m, and the sliders around it live here. A scene also gets `x`, `y`,
6+ * `z` (metres), `Lx`, `Ly`, `Lz`, `h` (metres), `c0` (m/s), `Ls` (the
7+ * string's length, metres) and the grid's counts for free — see
8+ * src/scene/scene.ts.
9+ */
10+import type { ParamSpec, Params } from '../mgpu/registry.ts';
11+import boxSource from '../../scenes/box.m?raw';
12+
13+export interface MScene {
14+ key: string;
15+ label: string;
16+ blurb: string;
17+ params: ParamSpec[];
18+ /** The string's height above the top plate, from this scene's parameters —
19+ * the host needs it to draw the string and place defaults, and the scene
20+ * needs it to build the line profile, so the scene's `gap` parameter is
21+ * the one shared source. */
22+ stringZ: (params: Params) => number;
23+ source: string;
24+}
25+
26+export const boxScene: MScene = {
27+ key: 'box',
28+ label: 'Box',
29+ blurb: 'A rigid box with a round sound hole in its top plate, under the string.',
30+ params: [
31+ {
32+ key: 'body',
33+ label: 'body (0 = removed)',
34+ value: 1,
35+ min: 0,
36+ max: 1,
37+ step: 1,
38+ hint: 'Set to 0 to take the box away entirely and hear the bare string in open air. The bridge patch stays, so the bridge-drive route still sounds.',
39+ },
40+ {
41+ key: 'boxl',
42+ label: 'body length (m)',
43+ value: 0.66,
44+ min: 0.3,
45+ max: 0.9,
46+ step: 0.01,
47+ hint: 'Along the string. The default just outreaches the string, so the bridge end sits over the body.',
48+ },
49+ { key: 'boxw', label: 'body width (m)', value: 0.16, min: 0.08, max: 0.3, step: 0.005 },
50+ {
51+ key: 'boxd',
52+ label: 'body depth (m)',
53+ value: 0.08,
54+ min: 0.03,
55+ max: 0.2,
56+ step: 0.005,
57+ hint: 'How far the box hangs below its top plate. Deeper means a bigger cavity and a lower air resonance.',
58+ },
59+ {
60+ key: 'thick',
61+ label: 'wall thickness (m)',
62+ value: 0.016,
63+ min: 0.008,
64+ max: 0.04,
65+ step: 0.002,
66+ hint: 'Below about two cells the plates start to leak — the grid cannot hold a thinner wall closed.',
67+ },
68+ {
69+ key: 'holer',
70+ label: 'sound hole radius (m)',
71+ value: 0.028,
72+ min: 0,
73+ max: 0.08,
74+ step: 0.002,
75+ hint: 'The opening in the top plate, toward the string. At 0 the box is sealed; the hole’s size against the cavity’s volume is what sets the air resonance (the Helmholtz pitch).',
76+ },
77+ {
78+ key: 'holex',
79+ label: 'sound hole x (m)',
80+ value: 0,
81+ min: -0.3,
82+ max: 0.3,
83+ step: 0.01,
84+ },
85+ {
86+ key: 'gap',
87+ label: 'string height (m)',
88+ value: 0.02,
89+ min: 0.01,
90+ max: 0.06,
91+ step: 0.002,
92+ hint: 'How far above the top plate the string runs.',
93+ },
94+ {
95+ key: 'patchr',
96+ label: 'bridge patch radius (m)',
97+ value: 0.05,
98+ min: 0.02,
99+ max: 0.12,
100+ step: 0.005,
101+ hint: 'Size of the top-plate region the bridge force drives (the stand-in for the moving plate).',
102+ },
103+ {
104+ key: 'absorb',
105+ label: 'wall absorption (1/s)',
106+ value: 800,
107+ min: 0,
108+ max: 2500,
109+ step: 25,
110+ hint: 'A thin lossy skin on the shell’s surfaces. At 0 the rigid cavity rings almost forever; turn it up until the box stops sounding metallic.',
111+ },
112+ ],
113+ stringZ: (params) => params.gap ?? 0.02,
114+ source: boxSource,
115+};
116+
117+export const mScenes: MScene[] = [boxScene];
118+
119+export const defaultSceneParams = (s: MScene): Params =>
120+ Object.fromEntries(s.params.map((p) => [p.key, p.value]));
src/scene/scene.tsadded+228−0View file
@@ -0,0 +1,228 @@
1+/**
2+ * The instrument's body and the coupling geometry: a .m scene file, evaluated
3+ * once on the air grid.
4+ *
5+ * A scene file is ordinary MATLAB defining one function,
6+ *
7+ * function [c, sig, wall, lineprof, boardprof] = medium(x, y, z, <...>)
8+ *
9+ * over the solver's grid points. Unlike the model it is *not* compiled to
10+ * WGSL: the step runs every timestep and must lower to a fixed sequence of
11+ * GPU dispatches, but a scene is evaluated exactly once and survives only as
12+ * five arrays of numbers. So it runs through numbl's CPU interpreter instead,
13+ * which buys the full MATLAB subset — loops, `if`, indexing, anything in
14+ * tools/ — and f64 evaluation, where the step dialect is element-wise f32.
15+ *
16+ * The five outputs:
17+ * c sound speed, m/s (normally just c0 everywhere — the walls are
18+ * a mask, not a material, so they cost no timestep)
19+ * sig absorption rate, 1/s (the sponge at the domain edge, plus any
20+ * absorption the body's surfaces are given)
21+ * wall 1 in air, 0 in the body's solid shell; what lapw masks by
22+ * lineprof where the string radiates directly (a tube around its line)
23+ * boardprof where the bridge force drives the air (a patch on the top plate)
24+ */
25+import { parseMFile, type FunctionStmt } from 'numbl-src/numbl-core/parser/index.ts';
26+import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
27+import {
28+ RuntimeTensor,
29+ isRuntimeTensor,
30+ type RuntimeValue,
31+} from 'numbl-src/numbl-core/runtime/types.ts';
32+import type { AirGrid } from '../grid.ts';
33+import { C_AIR } from '../units.ts';
34+import { toolFiles } from '../tools.ts';
35+import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
36+import type { MediumFields, ModelParams } from '../mgpu/model.ts';
37+
38+/** The function a scene file must define. */
39+export const MEDIUM_FN = 'medium';
40+
41+const OUTPUTS = ['c', 'sig', 'wall', 'lineprof', 'boardprof'] as const;
42+
43+export interface SceneOptions {
44+ air: AirGrid;
45+ /** String length, metres — the scene builds the coupling profiles around
46+ * the string, so it needs to know where the string is. */
47+ Ls: number;
48+ /** Scene source (.m text). */
49+ source: string;
50+ /** Parameter names the .m may take beyond the grid's own. */
51+ paramNames: string[];
52+ params: ModelParams;
53+}
54+
55+export class Scene implements MediumFields {
56+ readonly c: Float32Array;
57+ readonly sig: Float32Array;
58+ readonly wall: Float32Array;
59+ readonly lineprof: Float32Array;
60+ readonly boardprof: Float32Array;
61+ readonly cmin: number;
62+ readonly cmax: number;
63+ /** The background speed, taken from a corner of the domain — inside the
64+ * absorbing layer, where a scene has no business putting anything. */
65+ readonly cref: number;
66+
67+ private constructor(fields: Record<(typeof OUTPUTS)[number], Float32Array>) {
68+ this.c = fields.c;
69+ this.sig = fields.sig;
70+ this.wall = fields.wall;
71+ this.lineprof = fields.lineprof;
72+ this.boardprof = fields.boardprof;
73+ let lo = Infinity;
74+ let hi = 0;
75+ for (const v of this.c) {
76+ if (v < lo) lo = v;
77+ if (v > hi) hi = v;
78+ }
79+ this.cmin = lo;
80+ this.cmax = hi;
81+ this.cref = this.c[0];
82+ }
83+
84+ static create(opts: SceneOptions): Scene {
85+ const fields = evaluateScene(opts);
86+ const { c, sig, wall } = fields;
87+ for (let i = 0; i < c.length; i++) {
88+ if (!(c[i] > 0)) {
89+ throw new ModelCompileError(
90+ `the scene's sound speed is ${c[i]} somewhere; it must be positive ` +
91+ `everywhere (the timestep is set by the fastest point, and a zero ` +
92+ `or negative speed has no wave equation)`,
93+ { fn: MEDIUM_FN },
94+ );
95+ }
96+ if (!(sig[i] >= 0)) {
97+ throw new ModelCompileError(
98+ `the scene's absorption is ${sig[i]} somewhere; it must be zero or ` +
99+ `positive (a negative one would amplify rather than absorb)`,
100+ { fn: MEDIUM_FN },
101+ );
102+ }
103+ // The mask multiplies Laplacian fluxes; outside [0, 1] it would add
104+ // energy or invert a face. Clamp rather than refuse: a smoothed
105+ // difference of indicators dips a hair below zero in f64 routinely.
106+ if (wall[i] < 0) wall[i] = 0;
107+ else if (wall[i] > 1) wall[i] = 1;
108+ }
109+ return new Scene(fields);
110+ }
111+}
112+
113+/**
114+ * Evaluate the scene file on the grid, through numbl's CPU interpreter.
115+ *
116+ * The .m keeps the same contract the model has: it names the arguments it
117+ * wants — the coordinates, the domain's numbers, and any of the registry's
118+ * parameters — and the host supplies them by name, so their order in the
119+ * signature is the .m's own business.
120+ */
121+function evaluateScene(
122+ opts: SceneOptions,
123+): Record<(typeof OUTPUTS)[number], Float32Array> {
124+ const { air, Ls, source, paramNames, params } = opts;
125+ const file = `${MEDIUM_FN}.m`;
126+ const ast = inModel(() => parseMFile(source, file));
127+ const fn = ast.body.find(
128+ (s): s is FunctionStmt =>
129+ s.type === 'Function' && (s as FunctionStmt).name === MEDIUM_FN,
130+ );
131+ if (!fn) {
132+ throw new ModelCompileError(`the scene defines no function named '${MEDIUM_FN}'`);
133+ }
134+ if (fn.outputs.length !== OUTPUTS.length) {
135+ throw new ModelCompileError(
136+ `'${MEDIUM_FN}' must return ${OUTPUTS.length} outputs ` +
137+ `[${OUTPUTS.join(', ')}], not ${fn.outputs.length}`,
138+ { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
139+ );
140+ }
141+ // What the grid offers a scene by name, beyond its own parameters: the
142+ // coordinates in metres, the numbers that describe the domain, the speed
143+ // of sound in air, and the string's length — the scene builds the coupling
144+ // profiles around the string, so it needs to know where the string is.
145+ const vars: Record<string, RuntimeValue> = {
146+ x: new RuntimeTensor(air.x64, [air.npts, 1]),
147+ y: new RuntimeTensor(air.y64, [air.npts, 1]),
148+ z: new RuntimeTensor(air.z64, [air.npts, 1]),
149+ Lx: air.Lx,
150+ Ly: air.Ly,
151+ Lz: air.Lz,
152+ h: air.h,
153+ c0: C_AIR,
154+ npts: air.npts,
155+ nx: air.nx,
156+ ny: air.ny,
157+ nz: air.nz,
158+ Ls,
159+ };
160+ const known = new Set([...Object.keys(vars), ...paramNames]);
161+ for (const p of fn.params) {
162+ if (!known.has(p)) {
163+ throw new ModelCompileError(
164+ `'${MEDIUM_FN}' takes an argument '${p}' that is neither the grid ` +
165+ `(${Object.keys(vars).join(', ')}) nor one of this scene's parameters` +
166+ (paramNames.length ? ` (${paramNames.join(', ')})` : ''),
167+ { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
168+ );
169+ }
170+ }
171+
172+ for (const name of paramNames) {
173+ const v = params[name];
174+ // Missing parameters read as 0, as ModelPlan.setParams has it.
175+ vars[name] = Number.isFinite(v) ? v : 0;
176+ }
177+
178+ const outNames = OUTPUTS.map((o) => `${o}__`);
179+ const driver = `[${outNames.join(', ')}] = ${MEDIUM_FN}(${fn.params.join(', ')});`;
180+ const result = inFunction(MEDIUM_FN, () =>
181+ executeCode(
182+ driver,
183+ { initialVariableValues: vars, displayResults: false, implicitCwdPath: null },
184+ [...toolFiles, { name: file, source }],
185+ 'scene-driver.m',
186+ ),
187+ );
188+
189+ const fields = {} as Record<(typeof OUTPUTS)[number], Float32Array>;
190+ OUTPUTS.forEach((name, i) => {
191+ fields[name] = toGridField(
192+ result.variableValues[outNames[i]],
193+ fn.outputs[i],
194+ air.npts,
195+ );
196+ });
197+ return fields;
198+}
199+
200+/** One returned field -> npts values, rounded to the solver's f32. */
201+function toGridField(
202+ value: RuntimeValue | undefined,
203+ name: string,
204+ npts: number,
205+): Float32Array {
206+ // A uniform field stays scalar in MATLAB; spread it over the grid.
207+ if (typeof value === 'number') return new Float32Array(npts).fill(value);
208+ if (value !== undefined && isRuntimeTensor(value)) {
209+ if (value.imag) {
210+ throw new ModelCompileError(
211+ `the scene's '${name}' is complex; the medium must be real`,
212+ { fn: MEDIUM_FN },
213+ );
214+ }
215+ // A vector of npts values, either orientation. A reshape is refused
216+ // rather than reordered: the tensor's column-major layout would not match
217+ // the grid's x-fastest order.
218+ if (value.data.length === npts && value.shape.every((d) => d === 1 || d === npts)) {
219+ return new Float32Array(value.data);
220+ }
221+ throw new ModelCompileError(
222+ `the scene's '${name}' is ${value.shape.join(' x ')}, but the grid wants ` +
223+ `one value per point (${npts} x 1)`,
224+ { fn: MEDIUM_FN },
225+ );
226+ }
227+ throw new ModelCompileError(`the scene's '${name}' is not numeric`, { fn: MEDIUM_FN });
228+}
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+43−0View file
@@ -0,0 +1,43 @@
1+/**
2+ * Everything in SI: metres, seconds, hertz, metres per second.
3+ *
4+ * Physical units are what make the numbers here mean something. The string's
5+ * fundamental is a real pitch, the box is a real size you could build, the
6+ * timestep is a real duration, and the microphone's trace plays back at the
7+ * pitch a microphone there would have heard. They also keep the method's
8+ * limits visible: a grid solver resolves a wavelength with some number of
9+ * cells, so a fixed grid is a low-frequency method, and how low is a number
10+ * the app can show rather than hide.
11+ */
12+
13+/** Speed of sound in air at about 20 °C, m/s. */
14+export const C_AIR = 343;
15+
16+/** The domain, metres: a box around the instrument, twice as long along the
17+ * string (x) as across (y) and up (z), centred on the origin. The top plate
18+ * of the dulcimer body sits at z = 0 and the string runs along x just above
19+ * it. */
20+export const DOMAIN_X = 1.0;
21+export const DOMAIN_YZ = 0.5;
22+
23+/** Fraction of the stability limit the timestep is taken at. Fixed rather
24+ * than a control: everything downstream (the string grid, the recording's
25+ * sample rate) follows from dt, and one good value is worth more here than a
26+ * slider. */
27+export const CFL = 0.5;
28+
29+/** Cells per wavelength below which what is on screen (and in the recording)
30+ * is as much grid dispersion as it is sound. */
31+export const POOR_RESOLUTION = 8;
32+
33+/** A length in metres, written the way a person would say it. */
34+export const fmtLength = (m: number): string =>
35+ Math.abs(m) < 1 ? `${(1000 * m).toPrecision(3)} mm` : `${m.toPrecision(3)} m`;
36+
37+/** A duration in seconds, likewise. */
38+export const fmtTime = (s: number): string => {
39+ const a = Math.abs(s);
40+ if (a > 0 && a < 1e-3) return `${(1e6 * s).toPrecision(3)} µs`;
41+ if (a < 1) return `${(1e3 * s).toPrecision(3)} ms`;
42+ return `${s.toPrecision(3)} s`;
43+};
test/checks.tsadded+289−0View file
@@ -0,0 +1,289 @@
1+/**
2+ * The checks, against the real pipeline: MATLAB source -> numbl lowering ->
3+ * generated WGSL -> GPU. Physics first — does the string sound its pitch,
4+ * does it decay on schedule, does the sealed box actually seal — and the
5+ * planner's mechanics after.
6+ *
7+ * Everything runs on deliberately small grids: these are correctness checks,
8+ * and a 32-point air grid already carries every code path the 128-point one
9+ * does.
10+ */
11+import { ModelSession } from '../src/mgpu/session.ts';
12+import { dulcimerModel, defaultParams, type Params } from '../src/mgpu/registry.ts';
13+import { boxScene, defaultSceneParams } from '../src/scene/registry.ts';
14+import { planPlayback } from '../src/audio/play.ts';
15+
16+type Check = (name: string, ok: boolean, detail: string) => void;
17+type Log = (s: string) => void;
18+
19+const peakOf = (a: Float32Array): number => {
20+ let m = 0;
21+ for (const v of a) m = Math.max(m, Math.abs(v));
22+ return m;
23+};
24+
25+async function makeSession(
26+ device: GPUDevice,
27+ nx: number,
28+ params: Params = {},
29+ sceneParams: Params = {},
30+ operandBudget?: number,
31+): Promise<ModelSession> {
32+ return ModelSession.create({
33+ device,
34+ model: dulcimerModel,
35+ params: { ...defaultParams(dulcimerModel), ...params },
36+ scene: boxScene,
37+ sceneParams: { ...defaultSceneParams(boxScene), ...sceneParams },
38+ nx,
39+ Ls: 0.6,
40+ operandBudget,
41+ });
42+}
43+
44+/** The compiled plan: external ops present, pluck in place, mic riding along. */
45+export async function planChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
46+ const s = await makeSession(device, 32);
47+ const { step } = s.describe();
48+ for (const l of step) log(` ${l}`);
49+ const externals = step.filter((l) => l.startsWith('external'));
50+ const wanted = ['dxx(u)', 'dxx(um)', 'dxxxx(u)', 'spread(acc)', 'bridge(un)', 'lapw(p, wall)'];
51+ for (const w of wanted) {
52+ check(
53+ `the step uses ${w}`,
54+ externals.some((l) => l.includes(w)),
55+ externals.length ? externals.join('; ') : 'no external ops planned',
56+ );
57+ }
58+
59+ s.pluck();
60+ const u = await s.read('u');
61+ const amp = defaultParams(dulcimerModel).amp;
62+ const peak = peakOf(u);
63+ check(
64+ 'the pluck draws the string to its set height',
65+ Math.abs(peak - amp) < 0.15 * amp,
66+ `peak |u| = ${peak.toExponential(3)}, amp = ${amp}`,
67+ );
68+ check('the ends are pinned', u[0] === 0 && u[u.length - 1] === 0, `u[0]=${u[0]}, u[end]=${u[u.length - 1]}`);
69+
70+ s.step(100);
71+ check(
72+ 'the microphone samples once per timestep',
73+ s.recorder.count === 100,
74+ `${s.recorder.count} samples after 100 steps`,
75+ );
76+ const plan = planPlayback(s.recorder.count, s.dt);
77+ check(
78+ 'playback is real time at the solver rate',
79+ plan.realTime && Math.abs(plan.rate * s.dt - 1) < 1e-9,
80+ `rate ${plan.rate.toFixed(0)} Hz, dt ${s.dt.toExponential(3)}`,
81+ );
82+ s.destroy();
83+}
84+
85+/** The string sounds the pitch it is tuned to. */
86+export async function pitchChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
87+ // An ideal string (no stiffness, no damping to speak of), watched at a
88+ // point for a dozen periods; the zero crossings say the frequency.
89+ const f0 = 294;
90+ const s = await makeSession(device, 32, { f0, B: 0, sig1: 0, t60: 8 });
91+ s.pluck();
92+ const node = Math.round(s.string.ns * 0.4);
93+ const periods = 12;
94+ const stepsTotal = Math.round(periods / f0 / s.dt);
95+ const chunk = 10;
96+ const series: number[] = [];
97+ for (let done = 0; done < stepsTotal; done += chunk) {
98+ s.step(chunk);
99+ const u = await s.read('u');
100+ series.push(u[node]);
101+ }
102+ let crossings = 0;
103+ for (let i = 1; i < series.length; i++) {
104+ if ((series[i - 1] < 0 && series[i] >= 0) || (series[i - 1] >= 0 && series[i] < 0)) crossings++;
105+ }
106+ const measured = crossings / 2 / (stepsTotal * s.dt);
107+ log(` ${crossings} zero crossings over ${(stepsTotal * s.dt * 1000).toFixed(1)} ms -> ${measured.toFixed(1)} Hz`);
108+ check(
109+ `the string sounds its fundamental (${f0} Hz)`,
110+ Math.abs(measured - f0) < 0.04 * f0,
111+ `measured ${measured.toFixed(1)} Hz`,
112+ );
113+
114+ // d'Alembert: at half a period the string is the mirror of its pluck, so
115+ // the displacement at the pluck point flips sign and shrinks to the
116+ // triangle's value at the mirrored position.
117+ s.pluck();
118+ const u0 = await s.read('u');
119+ s.step(Math.round(1 / f0 / 2 / s.dt));
120+ const u1 = await s.read('u');
121+ const at = Math.round(0.22 * (s.string.ns - 1));
122+ check(
123+ 'half a period later the pluck point has swung through zero',
124+ u0[at] > 0 && u1[at] < 0,
125+ `u ${u0[at].toExponential(2)} -> ${u1[at].toExponential(2)}`,
126+ );
127+ s.destroy();
128+}
129+
130+/** The decay knob means what it says. */
131+export async function decayChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
132+ const t60 = 0.5;
133+ const s = await makeSession(device, 32, { t60, B: 0, sig1: 0 });
134+ s.pluck();
135+ const before = peakOf(await s.read('u'));
136+ // Half of t60: amplitude should be down 30 dB, i.e. to about 3.2%.
137+ s.step(Math.round(t60 / 2 / s.dt));
138+ const after = peakOf(await s.read('u'));
139+ const db = 20 * Math.log10(after / before);
140+ log(` |u| ${before.toExponential(2)} -> ${after.toExponential(2)} in ${t60 / 2} s (${db.toFixed(1)} dB)`);
141+ check(
142+ 't60 decays the string on schedule',
143+ Math.abs(db + 30) < 4,
144+ `${db.toFixed(1)} dB over t60/2, want -30`,
145+ );
146+ s.destroy();
147+}
148+
149+/** The wall mask is a wall: a sealed box keeps the sound out. */
150+export async function wallChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
151+ // Bridge drive only (the direct route off), so the source sits above the
152+ // plate; a microphone inside a sealed box should hear far less than one
153+ // outside beside it.
154+ const seconds = 0.06;
155+ const run = async (mic: [number, number, number]): Promise<number> => {
156+ const s = await makeSession(
157+ device,
158+ 64,
159+ { gline: 0, gbridge: 1 },
160+ { holer: 0, absorb: 0, boxd: 0.12, thick: 0.02 },
161+ );
162+ s.setMic(...mic);
163+ s.pluck();
164+ s.step(Math.round(seconds / s.dt));
165+ const trace = await s.recorder.read();
166+ s.destroy();
167+ return peakOf(trace);
168+ };
169+ const inside = await run([0, 0, -0.06]);
170+ const outside = await run([0, 0, 0.08]);
171+ log(` |p| inside the sealed box ${inside.toExponential(2)}, outside ${outside.toExponential(2)}`);
172+ check(
173+ 'a sealed box keeps the sound out',
174+ inside < 0.05 * outside,
175+ `inside/outside = ${(inside / outside).toExponential(2)}`,
176+ );
177+}
178+
179+/** The sound hole lets the cavity speak: opening it raises what gets in. */
180+export async function holeChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
181+ const seconds = 0.06;
182+ const run = async (holer: number): Promise<number> => {
183+ const s = await makeSession(
184+ device,
185+ 64,
186+ { gline: 0, gbridge: 1 },
187+ { holer, absorb: 0, boxd: 0.12, thick: 0.02 },
188+ );
189+ s.setMic(0, 0, -0.06);
190+ s.pluck();
191+ s.step(Math.round(seconds / s.dt));
192+ const trace = await s.recorder.read();
193+ s.destroy();
194+ return peakOf(trace);
195+ };
196+ const sealed = await run(0);
197+ const open = await run(0.05);
198+ log(` |p| in the cavity: sealed ${sealed.toExponential(2)}, open ${open.toExponential(2)}`);
199+ check(
200+ 'the sound hole lets the cavity speak',
201+ open > 5 * sealed,
202+ `open/sealed = ${(open / sealed).toExponential(2)}`,
203+ );
204+}
205+
206+/** Air sound sits on the string's partial comb, not somewhere else. */
207+export async function spectrumChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
208+ const f0 = 294;
209+ const s = await makeSession(device, 64, { f0 });
210+ s.setMic(0.12, 0.08, 0.1);
211+ s.pluck();
212+ const seconds = 0.12;
213+ const total = Math.round(seconds / s.dt);
214+ for (let done = 0; done < total; done += 2048) {
215+ s.step(Math.min(2048, total - done));
216+ await s.sync();
217+ }
218+ const trace = await s.recorder.read();
219+ const dft = (f: number): number => {
220+ let re = 0;
221+ let im = 0;
222+ for (let i = 0; i < trace.length; i++) {
223+ const w = 2 * Math.PI * f * i * s.dt;
224+ re += trace[i] * Math.cos(w);
225+ im -= trace[i] * Math.sin(w);
226+ }
227+ return Math.hypot(re, im) / trace.length;
228+ };
229+ const comb = (dft(f0) + dft(2 * f0) + dft(3 * f0)) / 3;
230+ const off = (dft(1.41 * f0) + dft(2.53 * f0)) / 2;
231+ log(` comb ${comb.toExponential(2)}, off-comb ${off.toExponential(2)}`);
232+ check(
233+ 'the microphone hears the string’s partials',
234+ comb > 8 * off,
235+ `comb/off = ${(comb / off).toFixed(1)}`,
236+ );
237+ s.destroy();
238+}
239+
240+/** Taking the body away leaves open air that still carries the string. */
241+export async function bareChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
242+ const s = await makeSession(device, 32, { gline: 1, gbridge: 0 }, { body: 0 });
243+ let minWall = 1;
244+ for (const v of s.scene.wall) minWall = Math.min(minWall, v);
245+ check(
246+ 'removing the body leaves pure air',
247+ minWall > 0.999,
248+ `min wall mask = ${minWall}`,
249+ );
250+ s.setMic(0.1, 0.05, 0.05);
251+ s.pluck();
252+ s.step(Math.round(0.03 / s.dt));
253+ const trace = await s.recorder.read();
254+ const peak = peakOf(trace);
255+ log(` bare-string trace peak ${peak.toExponential(2)}`);
256+ check('the bare string still sounds', peak > 0, `trace peak ${peak.toExponential(2)}`);
257+ s.destroy();
258+}
259+
260+/** A starved operand budget splits kernels instead of failing. */
261+export async function splitChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
262+ const s = await makeSession(device, 32, {}, {}, 2);
263+ const { step } = s.describe();
264+ const parts = step.filter((l) => l.includes('_part')).length;
265+ log(` ${parts} split kernels at budget 2`);
266+ check('a tight binding budget splits the update', parts > 0, `${parts} split kernels`);
267+ s.pluck();
268+ s.step(50);
269+ const p = await s.read('p');
270+ let bad = 0;
271+ for (const v of p) if (!Number.isFinite(v)) bad++;
272+ check('the split model still runs', bad === 0, `${bad} non-finite values`);
273+
274+ // The split result must agree with the unsplit one.
275+ const s2 = await makeSession(device, 32);
276+ s2.pluck();
277+ s2.step(50);
278+ const p2 = await s2.read('p');
279+ let worst = 0;
280+ const scale = Math.max(peakOf(p), 1e-30);
281+ for (let i = 0; i < p.length; i++) worst = Math.max(worst, Math.abs(p[i] - p2[i]));
282+ check(
283+ 'split and unsplit agree',
284+ worst < 1e-4 * scale,
285+ `worst |diff| = ${worst.toExponential(2)} of peak ${scale.toExponential(2)}`,
286+ );
287+ s.destroy();
288+ s2.destroy();
289+}
tools/sponge3.madded+23−0View file
@@ -0,0 +1,23 @@
1+% Absorption profile for an open boundary, in three dimensions.
2+%
3+% s = sponge3(x, y, z, Lx, Ly, Lz, 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 face of the box [-Lx/2, Lx/2] x [-Ly/2, Ly/2] x [-Lz/2, Lz/2]
8+% and reaching smax at the wall. Added to a scene's absorption, this is what
9+% makes the finite grid stand in for open air: a wave that leaves the
10+% instrument is attenuated before it reaches the outer boundary, and whatever
11+% reflects off 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 = sponge3(x, y, z, Lx, Ly, Lz, w, smax)
18+ dx = max(0, w - (Lx/2 - abs(x))) / w;
19+ dy = max(0, w - (Ly/2 - abs(y))) / w;
20+ dz = max(0, w - (Lz/2 - abs(z))) / w;
21+ d = max(max(dx, dy), dz);
22+ s = smax * d.^2;
23+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
9+// acoustic-scattering-2d and 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+});