Render any scene from precomputation, and split the Play knobs by cost
The solver's cost is the air's: ~114k timesteps per second of audio at the
default grid, whatever is playing. Almost all of that stepping is avoidable.
Every object's discrete operator has exact sine eigenmodes, the coupling is
one-way, and everything is linear, so each modal amplitude follows a scalar
recursion with a closed-form solution and the air sees a sum of fixed spatial
patterns times known scalar series. The mic trace is then those series
convolved with per-pattern impulse responses of the room, which depend only on
geometry.
Discrete reciprocity collapses the precomputation to a single run: the air
update's transfer matrix is symmetric (the masked Laplacian in divergence form
has w_i*w_j/h^2 off-diagonal entries, everything else diagonal), so the mic's
response to a pattern equals that pattern's projection onto the field one
impulse at the mic produces.
- src/studio/modal.ts: every object in closed form. Mode sets per type, with
the pluck and strike excitations as recursion initial conditions.
- src/studio/room.ts: one reciprocal air run, through the same compiler and
lapw op as the solver, projecting onto every pattern group each step.
- src/studio/fast.ts: orchestration, per-object routes (a string's bridge patch
included), and re-synthesis from a kept room.
- scripts/render.ts: --fast, --ir-seconds, --modal-eps, --fmax.
- A Fast render button in the transport bar, with the room cached on geometry
and the mic, and re-recorded when either changes.
Matches the timestepped solver to its own f32 rounding noise: 2.5e-3 for the
drum, 7.8e-5 for the dulcimer, 4.2e-5 for the kit, ~2e-4 for the rest. Twenty
seconds of drum takes 25 s rather than a projected 16 minutes, and every knob
afterwards re-renders in about a second. scripts/modal-poc.ts checks the
derivation itself to machine precision.
The Play tab now separates tuning from shape, since the geometry knobs cost
what tuning never does: a size recompiles, a body dimension re-evaluates the
scene, and either re-records the room.
12 changed files+2049−15
README.mdmodified+86−4View file
@@ -24,8 +24,13 @@ the controls they came for.
2424
2525 **Play** is presets and musical knobs. Pick a scene (dulcimer, drum, chime,
2626 a duet of strings, a drum-and-chime pair), excite it, render a clip, listen,
27-save a .wav. The sliders here are the handful a player would reach for:
28-pitch, decay, brightness, where the pluck lands.
27+save a .wav. The sliders come in two groups, split by what they cost. The
28+tuning group is the handful a player would reach for — pitch, decay,
29+brightness, where the pluck lands — and is free to move, even while the
30+scene plays. The shape group is the instrument itself: each object's size
31+and the body's main dimensions. Those change geometry, so a size recompiles
32+and either kind re-records a Fast render room; the group is set apart so
33+the two costs are never mixed.
2934
3035 **Build** is the studio. Add strings, bars and membranes; move each one in
3136 the air; choose and shape the body; place the microphone; set the air grid.
@@ -156,6 +161,81 @@ memory-bandwidth-bound, which is exactly what a discrete card is good at).
156161 Memory on the host grows as nx³ during scene evaluation; beyond nx = 384
157162 expect several gigabytes.
158163
164+## The fast path: instruments from precomputation
165+
166+A limitation of the renderer above is that its cost is the air's: about
167+114,000 timesteps per second of audio at the default grid, each a handful of
168+passes over the whole box, whatever is playing. Almost all of that stepping
169+is avoidable. Three facts about the models as compiled make it so.
170+Everything is linear and each excitation is an initial condition; the
171+coupling is one-way, so the air never pushes back on the objects; and every
172+object's discrete operator has exact sine eigenmodes (the same fact the
173+spectrum overlay already draws on, and for a bar the reason the
174+fourth-difference operator carries antisymmetric ghosts), so each modal
175+amplitude follows a scalar two-term recursion with a closed-form solution.
176+The source the air sees is then a sum of fixed spatial patterns times known
177+scalar series — one pattern per mode for an object radiating along its line
178+or plane, and for a string one more, the bridge patch its modes' arriving
179+slopes drive — and the microphone's trace is those series convolved with
180+per-pattern impulse responses of the room. The responses depend only on
181+geometry: the body, every object's size and placement, the grid, and the
182+microphone. Pitches, decays, brightness, the excitations and the gains
183+never touch them.
184+
185+One further fact collapses the precomputation to a single run. The air
186+update's transfer matrix is symmetric (the masked Laplacian in divergence
187+form has w_i w_j / h^2 off-diagonal entries, and everything else is
188+diagonal), so discrete reciprocity holds exactly: the mic's response to any
189+source pattern equals the projection of that pattern onto the field an
190+impulse at the mic produces. So
191+
192+```
193+npm run render -- --preset drum --seconds 20 --fast
194+```
195+
196+runs the air once, kicked at the microphone, through the same compiler and
197+the same lapw op as everything else, and projects the field onto every
198+pattern after each step (src/studio/room.ts, src/studio/modal.ts). The
199+projections are separable and ride the step's own submission. The trace is
200+then synthesized on the CPU, one complex one-pole per (mode, kernel) pair,
201+in about a second for any duration and any setting of the musical
202+parameters. Strings and bars have a few dozen modes each and keep all of
203+them; a membrane's set is bounded by `--modal-eps`, which with the kernel
204+length is the fast path's whole approximation. On the machine this was developed on, two seconds of drum at
205+the default grid took 96 s the ordinary way and 26 s with `--fast`; twenty
206+seconds took 25 s against a projected 16 minutes, since only the room
207+recording touches the air. Re-synthesizing from a kept room (see
208+`resynthesize` in src/studio/fast.ts) costs about a second, which is what
209+makes the head effectively playable.
210+
211+In the page, the transport bar's **Fast render** button is the same
212+machinery. The first press records the room, with a progress bar and a
213+cancel; every later press reuses it, so pitches, decays, brightness and the
214+excitations re-render a clip in about a second. The recording is kept until
215+something it depends on changes (the body, an object's size or placement,
216+the grid, the microphone), and a membrane strike moved far enough to want
217+modes outside the recording re-records by itself. Listen, Save .wav and the
218+spectrum all take the synthesized clip; the volume view stays with the live
219+solver, since the fast path produces no field to draw.
220+
221+In the app's cost language this is a fourth tier: model parameters stay
222+free, and it is geometry or microphone edits that re-record the room. The
223+approximations are stated, small, and reported: modes weaker than
224+`--modal-eps` times the strongest are dropped (default 1e-6; tightening it
225+to 1e-10 changes the output by parts in 1e6), and the room response is
226+truncated at `--ir-seconds` (default 0.5 s; the render prints how many dB
227+down the kernels are at that point). Within those, the fast path reproduces
228+the timestepped solver to its own f32 rounding noise — about 2e-3 relative
229+over two seconds of drum, rather less for the string and bar scenes — and
230+the test suite holds it there for the drum, the dulcimer and the kit,
231+retuned and restruck from a single recording. What it does not produce is
232+the pressure field itself, so the volume view still needs the live solver;
233+it needs the generated model .m (a hand-edited scene .m is fine, since a
234+scene is data either way); and it needs the uniform sound speed every
235+generated scene has, since a spatially varying c breaks the symmetry that
236+reciprocity rests on. The derivation is checked to machine precision in
237+scripts/modal-poc.ts.
238+
159239 numbl is a `file:../../numbl` dependency, expected as a sibling checkout two
160240 directories up (the layout of the concept-collection vault).
161241
@@ -164,5 +244,7 @@ the string's measured pitch and decay time, the bar landing on its discrete
164244 dispersion relation, the membrane's (1,1) mode, that a sealed box keeps
165245 sound out and a sound hole lets the cavity speak, that both strings of a
166246 duet reach the microphone, that moving an object moves its coupling profile
167-without a recompile, and that the planner's kernel splitting changes nothing
168-but the dispatch count.
247+without a recompile, that the planner's kernel splitting changes nothing
248+but the dispatch count, and that the fast path's modal synthesis matches the
249+timestepped drum, dulcimer and kit from one room recording per scene,
250+before and after retuning.
index.htmlmodified+11−1View file
@@ -248,9 +248,12 @@
248248 <option value="2" selected>2 s</option>
249249 <option value="3">3 s</option>
250250 <option value="5">5 s</option>
251+ <option value="10">10 s</option>
252+ <option value="20">20 s</option>
251253 </select>
252254 </label>
253255 <button id="renderclip" title="Excite, then run the solver flat out with no display until this much audio is recorded, then play it">Render clip</button>
256+ <button id="fastrender" title="Record the room's impulse responses once (about the cost of half a second of ordinary rendering), then synthesize the clip from every object's closed-form modes. Same numbers as the solver; after the first recording, every knob re-renders in about a second. Body, object size or placement, grid and mic edits re-record the room. Needs the generated model .m.">Fast render</button>
254257 <progress id="renderprog" value="0" max="1" hidden></progress>
255258 <button id="listen" title="Play what the microphone has recorded since the last excite, in real time at its real pitch">Listen</button>
256259 <button id="download" title="Save the recording as a 48 kHz .wav">Save .wav</button>
@@ -269,8 +272,15 @@
269272
270273 <section id="panel-play" class="tabpanel active">
271274 <div class="presets" id="presets"></div>
272- <p class="group-title">the knobs a player would reach for</p>
275+ <p class="group-title">tuning — free to move, even while it plays</p>
273276 <div class="sliders" id="musicalparams"></div>
277+ <p class="group-title">shape — the instrument itself</p>
278+ <p class="sub" style="margin-top:0" title="Sizes recompile the solver (a brief pause); body shape re-evaluates the scene in place. Either way the geometry changed, so the next Fast render re-records its room.">
279+ Sizes recompile; body shape applies in place. Both change the geometry,
280+ so the next Fast render re-records the room. The full set of shape and
281+ placement controls is in the Build tab.
282+ </p>
283+ <div class="sliders" id="shapeparams"></div>
274284 </section>
275285
276286 <section id="panel-build" class="tabpanel">
scripts/modal-poc.tsadded+455−0View file
@@ -0,0 +1,455 @@
1+/**
2+ * Numerical proof of concept for the modal fast path (drum).
3+ *
4+ * The claim to verify, before building anything: the drum's mic trace can be
5+ * reproduced, to floating-point accuracy, without timestepping the membrane
6+ * or (after a one-time precompute) the air. It rests on three facts about
7+ * the discrete scheme as compiled today:
8+ *
9+ * 1. The membrane's 5-point Laplacian with pinned rim has exact product-
10+ * sine eigenmodes, so each modal amplitude obeys a scalar two-term
11+ * recursion — no PDE stepping needed.
12+ * 2. The coupling is one-way and everything is linear, so the source the
13+ * air sees is a sum of fixed spatial patterns (prof .* rad(mode)) times
14+ * those scalar series.
15+ * 3. The air update's transfer matrix is symmetric (the masked Laplacian
16+ * in divergence form has w_i*w_j/h^2 off-diagonal entries; everything
17+ * else is diagonal), so discrete reciprocity holds exactly: one run
18+ * with an impulse at the mic yields every pattern-to-mic impulse
19+ * response at once.
20+ *
21+ * This script builds the real drum session at a small grid, runs the real
22+ * GPU solver, and then checks each link in the chain against f64 CPU
23+ * references:
24+ *
25+ * A. a CPU replica of the membrane leapfrog vs the GPU membrane state;
26+ * B. the modal (closed-form) membrane vs the CPU replica — this is the
27+ * "modes are exact" claim, and should agree to ~1e-12;
28+ * C. a CPU replica of the whole coupled step vs the GPU mic trace;
29+ * D. reciprocity on the air alone: direct pattern-source run vs the
30+ * projected reciprocal run — should agree to f64 roundoff;
31+ * E. end to end: mic trace synthesized from modal amplitudes convolved
32+ * with reciprocal kernels, vs the direct CPU trace and the GPU trace.
33+ *
34+ * Usage: npx vite-node scripts/modal-poc.ts [--nx 32] [--steps 1200]
35+ */
36+import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
37+import { requestAcousticDevice } from '../src/device.ts';
38+import { StudioSession, specParams } from '../src/studio/session.ts';
39+import { presetByKey, cloneSpec } from '../src/studio/presets.ts';
40+
41+const argOf = (name: string, dflt: number): number => {
42+ const i = process.argv.indexOf(name);
43+ return i >= 0 ? Number(process.argv[i + 1]) : dflt;
44+};
45+const NX = argOf('--nx', 32);
46+const STEPS = argOf('--steps', 1200);
47+
48+const rms = (a: Float64Array | Float32Array): number => {
49+ let s = 0;
50+ for (let i = 0; i < a.length; i++) s += a[i] * a[i];
51+ return Math.sqrt(s / a.length);
52+};
53+const relErr = (a: ArrayLike<number>, b: ArrayLike<number>): number => {
54+ let d = 0;
55+ let r = 0;
56+ for (let i = 0; i < a.length; i++) {
57+ d += (a[i] - b[i]) * (a[i] - b[i]);
58+ r += b[i] * b[i];
59+ }
60+ return Math.sqrt(d / Math.max(r, 1e-300));
61+};
62+const report = (label: string, err: number, tol: number): boolean => {
63+ const ok = err < tol;
64+ console.log(` ${ok ? 'ok ' : 'FAIL'} ${label}: rel err ${err.toExponential(2)} (tol ${tol})`);
65+ return ok;
66+};
67+
68+/* ------------------------------------------------------------- GPU run --- */
69+
70+try {
71+ await installWebGpu();
72+} catch (e) {
73+ console.error(errMsg(e), '\n', NO_ADAPTER_HINT);
74+ process.exit(1);
75+}
76+const device = await requestAcousticDevice();
77+const preset = presetByKey('drum')!;
78+const spec = cloneSpec(preset.spec);
79+const { params, sceneParams } = specParams(spec);
80+const session = await StudioSession.create({ device, spec, params, sceneParams, nx: NX });
81+
82+const air = session.air;
83+const g = session.grids['m1'];
84+const dt = session.dt;
85+const P = session.params;
86+console.log(
87+ `drum at nx = ${NX}: air ${air.nx}x${air.ny}x${air.nz}, membrane ${g.n1}x${g.n2}, ` +
88+ `dt = ${(1e6 * dt).toFixed(2)} µs, ${STEPS} steps`,
89+);
90+
91+session.setMic(preset.mic.x, preset.mic.y, preset.mic.z);
92+session.excite();
93+session.step(STEPS);
94+await session.sync();
95+const gpuTrace = await session.recorder.read();
96+const gpuU = await session.read('u_m1');
97+
98+/* ----------------------------------------------- CPU membrane leapfrog --- */
99+
100+const { n1, n2, h1, h2, L1, L2 } = g;
101+const mN = g.npts;
102+const pin = g.pin;
103+const f0 = P.f0_m1;
104+const t60 = P.t60_m1;
105+const sig1 = P.sig1_m1;
106+const gair = P.gair_m1;
107+const cm = (2 * f0) / Math.sqrt(1 / (L1 * L1) + 1 / (L2 * L2));
108+const dmp = (6.91 / t60) * dt;
109+
110+// The strike, exactly as the generated init writes it.
111+const um0 = new Float64Array(mN);
112+for (let k = 0; k < mN; k++) {
113+ const dx = g.xs64[k] - P.strikex_m1 * L1;
114+ const dy = g.ys64[k] - P.strikey_m1 * L2;
115+ const rr = (dx * dx + dy * dy) / (P.striker_m1 * P.striker_m1);
116+ um0[k] = -(dt * P.vel_m1) * pin[k] * Math.exp(-rr);
117+}
118+
119+const lap2 = (u: Float64Array, out: Float64Array): void => {
120+ const ih1 = 1 / (h1 * h1);
121+ const ih2 = 1 / (h2 * h2);
122+ for (let j = 0; j < n2; j++) {
123+ for (let i = 0; i < n1; i++) {
124+ const k = i + n1 * j;
125+ const c = u[k];
126+ const xm = i > 0 ? u[k - 1] : 0;
127+ const xp = i < n1 - 1 ? u[k + 1] : 0;
128+ const ym = j > 0 ? u[k - n1] : 0;
129+ const yp = j < n2 - 1 ? u[k + n1] : 0;
130+ out[k] = ih1 * (xm - 2 * c + xp) + ih2 * (ym - 2 * c + yp);
131+ }
132+ }
133+};
134+
135+// One membrane step: (u, um) -> un, matching the generated step line.
136+const memStep = (u: Float64Array, um: Float64Array, un: Float64Array, lu: Float64Array, lum: Float64Array): void => {
137+ lap2(u, lu);
138+ lap2(um, lum);
139+ for (let k = 0; k < mN; k++) {
140+ un[k] =
141+ pin[k] *
142+ ((2 * u[k] - (1 - dmp) * um[k] + dt * dt * cm * cm * lu[k] + 2 * sig1 * dt * (lu[k] - lum[k])) /
143+ (1 + dmp));
144+ }
145+};
146+
147+// Run the membrane alone, keeping every step's acceleration.
148+console.log('\nA/B. membrane: CPU leapfrog and closed-form modes');
149+let u = new Float64Array(mN);
150+let um = Float64Array.from(um0);
151+const accCpu: Float64Array[] = [];
152+{
153+ const lu = new Float64Array(mN);
154+ const lum = new Float64Array(mN);
155+ for (let k = 0; k < STEPS; k++) {
156+ const un = new Float64Array(mN);
157+ memStep(u, um, un, lu, lum);
158+ const acc = new Float64Array(mN);
159+ for (let i = 0; i < mN; i++) acc[i] = (un[i] - 2 * u[i] + um[i]) / (dt * dt);
160+ accCpu.push(acc);
161+ um = u;
162+ u = un;
163+ }
164+}
165+let okAll = report('CPU leapfrog vs GPU u_m1 (f32)', relErr(u, gpuU), 5e-3);
166+
167+/* -------------------------------------------------- modal membrane (B) --- */
168+
169+// Interior product-sine modes of the pinned 5-point Laplacian.
170+const M1 = n1 - 2;
171+const M2 = n2 - 2;
172+// Strike projection b[m][n] via two small separable contractions.
173+const sinTab1 = new Float64Array(M1 * n1); // sin(m*pi*i/(n1-1))
174+const sinTab2 = new Float64Array(M2 * n2);
175+for (let m = 1; m <= M1; m++)
176+ for (let i = 0; i < n1; i++) sinTab1[(m - 1) * n1 + i] = Math.sin((m * Math.PI * i) / (n1 - 1));
177+for (let n = 1; n <= M2; n++)
178+ for (let j = 0; j < n2; j++) sinTab2[(n - 1) * n2 + j] = Math.sin((n * Math.PI * j) / (n2 - 1));
179+
180+const projectModes = (v: Float64Array): Float64Array => {
181+ // v (n1 x n2) -> coefficients (M1 x M2), normalized.
182+ const tmp = new Float64Array(M1 * n2);
183+ for (let m = 0; m < M1; m++) {
184+ for (let j = 0; j < n2; j++) {
185+ let s = 0;
186+ for (let i = 0; i < n1; i++) s += sinTab1[m * n1 + i] * v[i + n1 * j];
187+ tmp[m + M1 * j] = s;
188+ }
189+ }
190+ const out = new Float64Array(M1 * M2);
191+ const norm = ((n1 - 1) / 2) * ((n2 - 1) / 2);
192+ for (let m = 0; m < M1; m++) {
193+ for (let n = 0; n < M2; n++) {
194+ let s = 0;
195+ for (let j = 0; j < n2; j++) s += sinTab2[n * n2 + j] * tmp[m + M1 * j];
196+ out[m + M1 * n] = s / norm;
197+ }
198+ }
199+ return out;
200+};
201+
202+const b = projectModes(um0); // alpha[-1] per mode; alpha[0] = 0.
203+
204+// Per-mode recursion, exactly the membrane step restricted to one mode.
205+const lamOf = (m: number, n: number): number => {
206+ const s1 = Math.sin((m * Math.PI) / (2 * (n1 - 1)));
207+ const s2 = Math.sin((n * Math.PI) / (2 * (n2 - 1)));
208+ return (4 / (h1 * h1)) * s1 * s1 + (4 / (h2 * h2)) * s2 * s2;
209+};
210+
211+// Keep modes the strike actually excites.
212+let bMax = 0;
213+for (const v of b) bMax = Math.max(bMax, Math.abs(v));
214+const kept: { m: number; n: number; lam: number; b: number }[] = [];
215+for (let m = 1; m <= M1; m++)
216+ for (let n = 1; n <= M2; n++)
217+ if (Math.abs(b[m - 1 + M1 * (n - 1)]) > 1e-10 * bMax)
218+ kept.push({ m, n, lam: lamOf(m, n), b: b[m - 1 + M1 * (n - 1)] });
219+console.log(` ${kept.length} of ${M1 * M2} modes carry the strike (rel cutoff 1e-10)`);
220+
221+// Modal alpha series and accelerations.
222+const alphaAcc = (lam: number, bm: number): { acc: Float64Array; last: number } => {
223+ const acc = new Float64Array(STEPS);
224+ let prev = bm; // alpha[k-1]
225+ let cur = 0; // alpha[k]
226+ const cA = (2 - dt * dt * cm * cm * lam - 2 * sig1 * dt * lam) / (1 + dmp);
227+ const cB = (-(1 - dmp) + 2 * sig1 * dt * lam) / (1 + dmp);
228+ for (let k = 0; k < STEPS; k++) {
229+ const next = cA * cur + cB * prev;
230+ acc[k] = (next - 2 * cur + prev) / (dt * dt);
231+ prev = cur;
232+ cur = next;
233+ }
234+ return { acc, last: cur };
235+};
236+
237+const modal = kept.map((md) => ({ ...md, ...alphaAcc(md.lam, md.b) }));
238+
239+// Reconstruct u after STEPS from modal amplitudes and compare to the CPU run.
240+{
241+ const uModal = new Float64Array(mN);
242+ for (const md of modal) {
243+ for (let j = 0; j < n2; j++) {
244+ const sj = sinTab2[(md.n - 1) * n2 + j];
245+ for (let i = 0; i < n1; i++) {
246+ uModal[i + n1 * j] += md.last * sinTab1[(md.m - 1) * n1 + i] * sj;
247+ }
248+ }
249+ }
250+ okAll = report('modal membrane vs CPU leapfrog', relErr(uModal, u), 1e-9) && okAll;
251+}
252+
253+/* ------------------------------------------------------- CPU air replica --- */
254+
255+console.log('\nC. coupled CPU replica vs the GPU mic trace');
256+const { nx, ny, nz, npts, h } = air;
257+const F = session.scene.fields;
258+const wall = F.wall;
259+const sigF = F.sig;
260+const cF = F.c;
261+const prof = F.prof_m1;
262+
263+// The mic index, as setMic computes it.
264+const micIx = Math.max(0, Math.min(nx - 1, Math.round((preset.mic.x + air.Lx / 2) / h - 0.5)));
265+const micIy = Math.max(0, Math.min(ny - 1, Math.round((preset.mic.y + air.Ly / 2) / h - 0.5)));
266+const micIz = Math.max(0, Math.min(nz - 1, Math.round((preset.mic.z + air.Lz / 2) / h - 0.5)));
267+const micIdx = micIx + nx * (micIy + ny * micIz);
268+
269+// Bilinear sampling weights of the membrane at each air column (x, y) —
270+// exactly what rad_m1 bakes into its WGSL.
271+interface RadW { k00: number; k10: number; k01: number; k11: number; w00: number; w10: number; w01: number; w11: number }
272+const radW: (RadW | null)[] = [];
273+for (let iy = 0; iy < ny; iy++) {
274+ for (let ix = 0; ix < nx; ix++) {
275+ const xw = -air.Lx / 2 + (ix + 0.5) * h;
276+ const yw = -air.Ly / 2 + (iy + 0.5) * h;
277+ const uu = (xw + L1 / 2) / h1;
278+ const vv = (yw + L2 / 2) / h2;
279+ const ju = Math.floor(uu);
280+ const jv = Math.floor(vv);
281+ if (ju < 0 || ju >= n1 - 1 || jv < 0 || jv >= n2 - 1) {
282+ radW[ix + nx * iy] = null;
283+ continue;
284+ }
285+ const fu = uu - ju;
286+ const fv = vv - jv;
287+ radW[ix + nx * iy] = {
288+ k00: ju + n1 * jv, k10: ju + 1 + n1 * jv, k01: ju + n1 * (jv + 1), k11: ju + 1 + n1 * (jv + 1),
289+ w00: (1 - fu) * (1 - fv), w10: fu * (1 - fv), w01: (1 - fu) * fv, w11: fu * fv,
290+ };
291+ }
292+}
293+const radSample = (v: Float64Array, col: number): number => {
294+ const w = radW[col];
295+ if (!w) return 0;
296+ return w.w00 * v[w.k00] + w.w10 * v[w.k10] + w.w01 * v[w.k01] + w.w11 * v[w.k11];
297+};
298+
299+// One air step; drv is the already-assembled source field (or null).
300+const airStep = (p: Float64Array, pm: Float64Array, pn: Float64Array, drv: ((i: number) => number) | null): void => {
301+ const ih2 = 1 / (h * h);
302+ for (let iz = 0; iz < nz; iz++) {
303+ for (let iy = 0; iy < ny; iy++) {
304+ for (let ix = 0; ix < nx; ix++) {
305+ const i = ix + nx * (iy + ny * iz);
306+ const pc = p[i];
307+ let s = 0;
308+ s += ix > 0 ? wall[i - 1] * (p[i - 1] - pc) : -pc;
309+ s += ix < nx - 1 ? wall[i + 1] * (p[i + 1] - pc) : -pc;
310+ s += iy > 0 ? wall[i - nx] * (p[i - nx] - pc) : -pc;
311+ s += iy < ny - 1 ? wall[i + nx] * (p[i + nx] - pc) : -pc;
312+ s += iz > 0 ? wall[i - nx * ny] * (p[i - nx * ny] - pc) : -pc;
313+ s += iz < nz - 1 ? wall[i + nx * ny] * (p[i + nx * ny] - pc) : -pc;
314+ const lp = ih2 * wall[i] * s;
315+ const sd = sigF[i] * dt;
316+ const src = drv ? drv(i) : 0;
317+ pn[i] =
318+ wall[i] *
319+ ((2 * pc - (1 - sd) * pm[i] + cF[i] * dt * cF[i] * dt * lp + dt * dt * src) / (1 + sd));
320+ }
321+ }
322+ }
323+};
324+
325+// The coupled run: modal membrane is already verified ≡ leapfrog, so drive
326+// the air from the stored leapfrog accelerations.
327+const cpuTrace = new Float64Array(STEPS);
328+{
329+ let p = new Float64Array(npts);
330+ let pm = new Float64Array(npts);
331+ let pn = new Float64Array(npts);
332+ for (let k = 0; k < STEPS; k++) {
333+ const acc = accCpu[k];
334+ airStep(p, pm, pn, (i) => gair * prof[i] * radSample(acc, i % (nx * ny)));
335+ cpuTrace[k] = pn[micIdx];
336+ const t = pm;
337+ pm = p;
338+ p = pn;
339+ pn = t;
340+ }
341+}
342+okAll = report('CPU coupled trace vs GPU trace (f32)', relErr(cpuTrace, gpuTrace.subarray(0, STEPS)), 5e-3) && okAll;
343+
344+/* ---------------------------------------------------- reciprocity (D/E) --- */
345+
346+console.log('\nD. discrete reciprocity on the air');
347+
348+// Reciprocal run: unit source at the mic at step 0, projected each step onto
349+// every kept mode's pattern prof .* rad(phi). The projection is separable:
350+// contract z first (plane = sum_z prof .* p), then the two sine factors.
351+// rad(phi_mn)(ix,iy) = sxA[m][ix] * syA[n][iy], because bilinear
352+// interpolation of a product of sines factorizes per axis.
353+const sxA = new Float64Array(M1 * nx);
354+const syA = new Float64Array(M2 * ny);
355+for (let ix = 0; ix < nx; ix++) {
356+ const xw = -air.Lx / 2 + (ix + 0.5) * h;
357+ const uu = (xw + L1 / 2) / h1;
358+ const ju = Math.floor(uu);
359+ if (ju < 0 || ju >= n1 - 1) continue;
360+ const fu = uu - ju;
361+ for (let m = 1; m <= M1; m++)
362+ sxA[(m - 1) * nx + ix] = (1 - fu) * sinTab1[(m - 1) * n1 + ju] + fu * sinTab1[(m - 1) * n1 + ju + 1];
363+}
364+for (let iy = 0; iy < ny; iy++) {
365+ const yw = -air.Ly / 2 + (iy + 0.5) * h;
366+ const vv = (yw + L2 / 2) / h2;
367+ const jv = Math.floor(vv);
368+ if (jv < 0 || jv >= n2 - 1) continue;
369+ const fv = vv - jv;
370+ for (let n = 1; n <= M2; n++)
371+ syA[(n - 1) * ny + iy] = (1 - fv) * sinTab2[(n - 1) * n2 + jv] + fv * sinTab2[(n - 1) * n2 + jv + 1];
372+}
373+
374+const kernels = modal.map(() => new Float64Array(STEPS));
375+{
376+ let p = new Float64Array(npts);
377+ let pm = new Float64Array(npts);
378+ let pn = new Float64Array(npts);
379+ const plane = new Float64Array(nx * ny);
380+ const tmp = new Float64Array(M1 * ny);
381+ for (let k = 0; k < STEPS; k++) {
382+ airStep(p, pm, pn, k === 0 ? (i) => (i === micIdx ? 1 : 0) : null);
383+ // plane(ix,iy) = sum_z prof * pn
384+ plane.fill(0);
385+ for (let iz = 0; iz < nz; iz++) {
386+ const base = nx * ny * iz;
387+ for (let c2 = 0; c2 < nx * ny; c2++) plane[c2] += prof[base + c2] * pn[base + c2];
388+ }
389+ // tmp(m,iy) = sum_ix sxA[m][ix] * plane(ix,iy)
390+ tmp.fill(0);
391+ for (let m = 0; m < M1; m++) {
392+ for (let iy = 0; iy < ny; iy++) {
393+ let s = 0;
394+ for (let ix = 0; ix < nx; ix++) s += sxA[m * nx + ix] * plane[ix + nx * iy];
395+ tmp[m + M1 * iy] = s;
396+ }
397+ }
398+ for (let q = 0; q < modal.length; q++) {
399+ const { m, n } = modal[q];
400+ let s = 0;
401+ for (let iy = 0; iy < ny; iy++) s += syA[(n - 1) * ny + iy] * tmp[m - 1 + M1 * iy];
402+ kernels[q][k] = s;
403+ }
404+ const t = pm;
405+ pm = p;
406+ p = pn;
407+ pn = t;
408+ }
409+}
410+
411+// Direct check for one mode: fire its pattern as an impulse, listen at the
412+// mic; must equal the reciprocal kernel to f64 roundoff.
413+{
414+ const q = modal.findIndex((md) => md.m === 2 && md.n === 3);
415+ const md = modal[q >= 0 ? q : 0];
416+ const pat = new Float64Array(npts);
417+ for (let i = 0; i < npts; i++) {
418+ const col = i % (nx * ny);
419+ pat[i] = prof[i] * sxA[(md.m - 1) * nx + (col % nx)] * syA[(md.n - 1) * ny + ((col / nx) | 0)];
420+ }
421+ const direct = new Float64Array(STEPS);
422+ let p = new Float64Array(npts);
423+ let pm = new Float64Array(npts);
424+ let pn = new Float64Array(npts);
425+ for (let k = 0; k < STEPS; k++) {
426+ airStep(p, pm, pn, k === 0 ? (i) => pat[i] : null);
427+ direct[k] = pn[micIdx];
428+ const t = pm;
429+ pm = p;
430+ p = pn;
431+ pn = t;
432+ }
433+ okAll = report(`reciprocity, mode (${md.m},${md.n})`, relErr(kernels[q >= 0 ? q : 0], direct), 1e-10) && okAll;
434+}
435+
436+/* ------------------------------------------------- modal synthesis (E) --- */
437+
438+console.log('\nE. modal synthesis vs direct traces');
439+const yModal = new Float64Array(STEPS);
440+for (let q = 0; q < modal.length; q++) {
441+ const hk = kernels[q];
442+ const a = modal[q].acc;
443+ // y[t] += gair * sum_k hk[t-k] * a[k]
444+ for (let t = 0; t < STEPS; t++) {
445+ let s = 0;
446+ for (let k = 0; k <= t; k++) s += hk[t - k] * a[k];
447+ yModal[t] += gair * s;
448+ }
449+}
450+okAll = report('modal synthesis vs CPU coupled trace', relErr(yModal, cpuTrace), 1e-6) && okAll;
451+okAll = report('modal synthesis vs GPU trace (f32)', relErr(yModal, gpuTrace.subarray(0, STEPS)), 5e-3) && okAll;
452+
453+console.log(okAll ? '\nall checks passed' : '\nSOME CHECKS FAILED');
454+session.destroy();
455+process.exit(okAll ? 0 : 1);
scripts/render.tsmodified+101−0View file
@@ -14,6 +14,7 @@
1414 *
1515 * npm run render -- scene.json --seconds 10 --out clip.wav
1616 * npm run render -- --preset drum --seconds 5 --nx 192
17+ * npm run render -- --preset drum --seconds 20 --fast
1718 * npm run render -- --list-presets
1819 *
1920 * Options:
@@ -21,6 +22,20 @@
2122 * --nx N override the air grid (even; memory grows as nx^3)
2223 * --mic x,y,z override the microphone position, metres
2324 * --out FILE output path (default: the scene's or preset's name.wav)
25+ * --fast the modal path: record the room's impulse responses
26+ * once by reciprocity, synthesize the trace from every
27+ * object's closed-form modes. Same numbers as the full
28+ * solver up to the options below; the win grows with
29+ * --seconds, since only the room recording touches the
30+ * air grid. See src/studio/fast.ts. Needs the generated
31+ * model .m (a hand-edited scene .m is fine).
32+ * --ir-seconds N fast: kernel (room response) length (default 0.5, or the
33+ * full duration if shorter). Raise it if the decay figure
34+ * the render prints looks high.
35+ * --modal-eps X fast: drop membrane modes weaker than X times the
36+ * strongest (default 1e-6; 0 keeps every mode the strike
37+ * touches; strings and bars always keep every mode)
38+ * --fmax F fast: drop membrane modes above F Hz (default: keep all)
2439 *
2540 * Runs through vite-node (which resolves numbl's compiler sources) on
2641 * desktop WebGPU via the optional `webgpu` package (Google Dawn). On a
@@ -36,6 +51,7 @@ import { parseSceneFile, makeSceneFile, type SceneFile } from '../src/studio/sce
3651 import { PRESETS, cloneSpec, presetByKey } from '../src/studio/presets.ts';
3752 import { specParams } from '../src/studio/session.ts';
3853 import { stableDt } from '../src/grids.ts';
54+import { fastIneligible, fastRender } from '../src/studio/fast.ts';
3955 import { traceToWav, WAV_RATE } from '../src/audio/wav.ts';
4056 import { C_AIR, CFL, DOMAIN_X, POOR_RESOLUTION, fmtTime } from '../src/units.ts';
4157
@@ -57,6 +73,10 @@ let seconds = 5;
5773 let nxOverride: number | null = null;
5874 let micOverride: { x: number; y: number; z: number } | null = null;
5975 let outPath: string | null = null;
76+let fast = false;
77+let irSeconds: number | null = null;
78+let modalEps = 1e-6;
79+let modalFmax = Infinity;
6080
6181 for (let i = 0; i < args.length; i++) {
6282 const a = args[i];
@@ -72,6 +92,10 @@ for (let i = 0; i < args.length; i++) {
7292 else if (a === '--seconds') seconds = Number(next());
7393 else if (a === '--nx') nxOverride = Number(next());
7494 else if (a === '--out') outPath = next();
95+ else if (a === '--fast') fast = true;
96+ else if (a === '--ir-seconds') irSeconds = Number(next());
97+ else if (a === '--modal-eps') modalEps = Number(next());
98+ else if (a === '--fmax') modalFmax = Number(next());
7599 else if (a === '--mic') {
76100 const parts = next().split(',').map(Number);
77101 if (parts.length !== 3 || parts.some((v) => !Number.isFinite(v))) usage();
@@ -135,6 +159,83 @@ try {
135159 process.exit(1);
136160 }
137161
162+/* ---------------------------------------------------- the fast path ----- */
163+
164+if (fast) {
165+ const bad = fastIneligible(file.spec, file.sources?.model);
166+ if (bad) {
167+ console.error(`--fast: ${bad}. Rendering without --fast works for every scene.`);
168+ process.exit(1);
169+ }
170+ const defaults = specParams(file.spec);
171+ const params = { ...defaults.params, ...file.params };
172+ const sceneParams = { ...defaults.sceneParams, ...file.sceneParams };
173+ const dt = stableDt(DOMAIN_X / nx, C_AIR, CFL);
174+ const steps = Math.ceil(seconds / dt);
175+ const ir = irSeconds ?? Math.min(seconds, 0.5);
176+
177+ const t0 = performance.now();
178+ let lastPct = -1;
179+ const result = await fastRender({
180+ device,
181+ spec: file.spec,
182+ params,
183+ sceneParams,
184+ sceneSource: file.sources?.scene,
185+ nx,
186+ mic,
187+ steps,
188+ irSeconds: ir,
189+ modeOpts: { eps: modalEps, fmax: modalFmax },
190+ onProgress: (done, total) => {
191+ const pct = Math.floor((10 * done) / total);
192+ if (pct > lastPct) {
193+ lastPct = pct;
194+ process.stderr.write(` room ${String(10 * pct).padStart(3)}%\n`);
195+ }
196+ },
197+ }).catch((e: unknown) => {
198+ console.error(`--fast failed: ${errMsg(e)}`);
199+ process.exit(1);
200+ });
201+
202+ const { air } = result;
203+ const fmaxHonest = C_AIR / (POOR_RESOLUTION * air.h);
204+ const decayDb = result.kernelDecay > 0 ? 20 * Math.log10(result.kernelDecay) : -Infinity;
205+ const objects = result.objects
206+ .map((or) => {
207+ const g = or.grid;
208+ const size = g.kind === 'line' ? `${g.n1} nodes` : `${g.n1}×${g.n2} nodes`;
209+ return `${or.id} (${size}, ${or.modes.length} modes${or.bridgeRow !== null ? ' + bridge' : ''})`;
210+ })
211+ .join(', ');
212+ console.log(
213+ `air ${air.nx}×${air.ny}×${air.nz} (honest to ~${(fmaxHonest / 1000).toFixed(1)} kHz) · ` +
214+ `in closed form: ${objects}\n` +
215+ `room ${fmtTime(result.room.steps * result.dt)} by reciprocity ` +
216+ `(${result.room.rows} kernels), ` +
217+ `${decayDb === -Infinity ? 'silent' : `${decayDb.toFixed(0)} dB down`} at its end` +
218+ (decayDb > -60 && result.room.steps < steps
219+ ? ' — still ringing; consider a longer --ir-seconds'
220+ : '') +
221+ (result.missing > 0 ? ` · ${result.missing} wanted modes unrecorded` : ''),
222+ );
223+ let peak = 0;
224+ for (const v of result.trace) peak = Math.max(peak, Math.abs(v));
225+ if (!(peak > 0)) {
226+ console.error('warning: the trace is silence — is the gain zero, or the mic inside a wall?');
227+ }
228+ const wav = traceToWav(result.trace, result.dt);
229+ await writeFile(out, Buffer.from(await wav.arrayBuffer()));
230+ const wall = (performance.now() - t0) / 1000;
231+ console.log(
232+ `${fmtTime(steps * result.dt)} of audio in ${wall.toFixed(1)} s ` +
233+ `(${((steps * result.dt) / wall).toFixed(2)}× real time) -> ` +
234+ `${out} (${WAV_RATE / 1000} kHz, ${(steps * result.dt).toFixed(2)} s)`,
235+ );
236+ process.exit(0);
237+}
238+
138239 // Size the trace for the request before the session exists. The estimate
139240 // assumes the medium tops out at the speed of air, which holds for every
140241 // generated scene; a hand-edited faster medium shrinks dt, and the shortfall
scripts/test-node.tsmodified+2−0View file
@@ -21,6 +21,7 @@ import {
2121 placementChecks,
2222 spectrumChecks,
2323 splitChecks,
24+ modalChecks,
2425 } from '../test/checks.ts';
2526
2627 let failures = 0;
@@ -64,6 +65,7 @@ await duetChecks(device, check, log);
6465 await placementChecks(device, check, log);
6566 await spectrumChecks(device, check, log);
6667 await splitChecks(device, check, log);
68+await modalChecks(device, check, log);
6769
6870 console.log(failures ? `\n${failures} failure(s)` : '\nall checks passed');
6971 process.exit(failures ? 1 : 0);
src/main.tsmodified+208−2View file
@@ -28,6 +28,9 @@ import {
2828 } from './studio/types.ts';
2929 import { suffixedSpecs } from './studio/codegen.ts';
3030 import { makeSceneFile, parseSceneFile, serializeSceneFile } from './studio/scenefile.ts';
31+import { fastIneligible, prepareRoom, roomContributions, type FastRoom } from './studio/fast.ts';
32+import { synthesize } from './studio/modal.ts';
33+import { RoomCancelled } from './studio/room.ts';
3134 import {
3235 drawSpectrum, objectColor, predictedPartials, type PartialSet,
3336 } from './studio/spectrum.ts';
@@ -96,6 +99,15 @@ let plotId: string | null = null;
9699 let lastTrace: { trace: Float32Array; dt: number } | null = null;
97100 let dirty = true;
98101 let renderingClip = false;
102+/** The fast path (src/studio/fast.ts): the cached room recording, keyed on
103+ * everything that shapes it — the fourth cost tier. Model parameters reuse
104+ * it; geometry, grid or mic edits key differently and re-record. */
105+let fastRoom: { key: string; fr: FastRoom } | null = null;
106+let fastBusy = false;
107+let fastCancel = false;
108+/** The last fast-synthesized clip. Live recording replaces it: whatever was
109+ * recorded or synthesized most recently is what Listen and Save play. */
110+let fastTrace: { trace: Float32Array; dt: number; note: string } | null = null;
99111
100112 /* ------------------------------------------------------------- the page -- */
101113
@@ -306,6 +318,38 @@ function buildMusicalSliders(): void {
306318 }
307319 }
308320
321+/** The Play tab's geometry knobs — deliberately a separate group from the
322+ * tuning sliders, because they cost what tuning never does: an object size
323+ * recompiles, a body shape re-evaluates the scene, and either way the next
324+ * Fast render re-records its room. Only the specs marked musical appear;
325+ * the full set lives in the Build tab. */
326+function buildShapeSliders(): void {
327+ const host = el('shapeparams');
328+ host.textContent = '';
329+ for (const o of spec.objects) {
330+ const t = templateOf(o);
331+ const sizeSpecs = t.sizeSpecs
332+ .filter((s) => s.musical)
333+ .map((s) => ({ ...s, label: `${o.id} · ${s.label}` }));
334+ addSliders(
335+ host,
336+ sizeSpecs,
337+ Object.fromEntries(t.sizeSpecs.map((s) => [s.key, (o.size as Params)[s.key] ?? s.value])),
338+ (key, value) => {
339+ writeBack();
340+ (o.size as Params)[key] = value;
341+ presetKey = null;
342+ scheduleRebuild();
343+ },
344+ );
345+ }
346+ const body = BODY_TEMPLATES[spec.body.type];
347+ const bodySpecs = body.paramSpecs
348+ .filter((s) => s.musical)
349+ .map((s) => ({ ...s, label: `body · ${s.label}` }));
350+ addSliders(host, bodySpecs, sceneParams, onSceneParam);
351+}
352+
309353 function buildObjectCards(): void {
310354 const host = el('objects');
311355 host.textContent = '';
@@ -401,6 +445,7 @@ function refreshControls(): void {
401445 el<HTMLButtonElement>('excite').textContent = verb();
402446 buildPresetCards();
403447 buildMusicalSliders();
448+ buildShapeSliders();
404449 buildObjectCards();
405450 buildBodyControls();
406451 buildMicControls();
@@ -418,6 +463,7 @@ let reading = false;
418463 function excite(): void {
419464 if (!session) return;
420465 session.excite();
466+ fastTrace = null;
421467 peakSeen = 0;
422468 if (autoScale) scale = 1e-6;
423469 plotPeak.clear();
@@ -470,6 +516,7 @@ async function rebuild(): Promise<void> {
470516 peakSeen = 0;
471517 scale = 1e-6;
472518 lastTrace = null;
519+ fastTrace = null;
473520 dirty = true;
474521 clearError();
475522 showRecording();
@@ -524,6 +571,15 @@ const describe = (s: StudioSession): string => {
524571 /** What the microphone has, and what it would sound like. */
525572 function showRecording(): void {
526573 const info = el('recinfo');
574+ if (fastTrace) {
575+ const { trace, dt, note } = fastTrace;
576+ const plan = planPlayback(trace.length, dt);
577+ info.textContent =
578+ `${fmtTime(trace.length * dt)} synthesized ${note} · plays at ` +
579+ `${Math.round(plan.rate).toLocaleString()} Hz` +
580+ (plan.realTime ? '' : ' (rate clamped)');
581+ return;
582+ }
527583 if (!session) {
528584 info.textContent = '';
529585 return;
@@ -756,7 +812,7 @@ async function analyzeTrace(): Promise<Float32Array | null> {
756812 * the trace. Then play it.
757813 */
758814 async function renderClip(): Promise<void> {
759- if (!session || building) return;
815+ if (!session || building || fastBusy) return;
760816 const btn = el<HTMLButtonElement>('renderclip');
761817 if (renderingClip) {
762818 renderingClip = false; // cancel: the loop below notices and stops
@@ -801,11 +857,132 @@ async function renderClip(): Promise<void> {
801857 }
802858 }
803859
860+/* ------------------------------------------------- the fast render path -- */
861+
862+/** Everything the room recording depends on. A matching key means the cached
863+ * kernels are still the truth; model parameters are deliberately absent. */
864+const fastKey = (): string =>
865+ JSON.stringify({
866+ objects: spec.objects.map((o) => ({ id: o.id, type: o.type, size: o.size })),
867+ body: spec.body.type,
868+ sceneParams,
869+ nx: gridNx,
870+ mic: { x: mic.x, y: mic.y, z: mic.z },
871+ scene: sources.scene,
872+ });
873+
874+/**
875+ * Render the clip from precomputation: record the room once (the only slow
876+ * step — about half a second of ordinary rendering), then synthesize the mic
877+ * trace from the membrane's closed-form modes against the cached kernels.
878+ * Any later change to pitch, decay, brightness, the strike or the gains
879+ * reuses the room and lands here in about a second; geometry, grid and mic
880+ * edits change the key and re-record. The volume view is untouched: the fast
881+ * path produces the trace and nothing else, so the field on screen stays
882+ * whatever the live solver last drew.
883+ */
884+async function fastRenderClip(): Promise<void> {
885+ const btn = el<HTMLButtonElement>('fastrender');
886+ if (fastBusy) {
887+ fastCancel = true; // cancel: the recording and synthesis loops notice
888+ return;
889+ }
890+ if (building || renderingClip) return;
891+ const bad = fastIneligible(spec, sources.model ?? undefined);
892+ if (bad) {
893+ showError(new Error(`Fast render: ${bad}.`), '');
894+ return;
895+ }
896+ fastBusy = true;
897+ fastCancel = false;
898+ setRunning(false);
899+ btn.textContent = 'Cancel';
900+ const prog = el<HTMLProgressElement>('renderprog');
901+ prog.hidden = false;
902+ prog.value = 0;
903+ const seconds = Number(el<HTMLSelectElement>('renderdur').value);
904+ try {
905+ syncParams();
906+ const P = { ...params };
907+ const key = fastKey();
908+
909+ // Two attempts at most: the second only if a new strike reaches well
910+ // outside the cached recording's membrane mode set, which discards it.
911+ for (let attempt = 0; attempt < 2; attempt++) {
912+ let synthBase = 0;
913+ if (!fastRoom || fastRoom.key !== key) {
914+ fastRoom = null;
915+ synthBase = 0.9;
916+ const fr = await prepareRoom({
917+ device: device!,
918+ spec,
919+ params: P,
920+ sceneParams: { ...sceneParams },
921+ sceneSource: sources.scene ?? undefined,
922+ nx: gridNx,
923+ mic: { ...mic },
924+ irSeconds: 0.5,
925+ // Membranes recorded a decade deeper than synthesis keeps, so
926+ // later strikes rarely want modes the recording lacks (strings
927+ // and bars record every mode and never can).
928+ modeOpts: { eps: 1e-7 },
929+ onProgress: (done, total) => {
930+ prog.value = 0.9 * (done / total);
931+ },
932+ cancelled: () => fastCancel,
933+ });
934+ fastRoom = { key, fr };
935+ }
936+
937+ const { fr } = fastRoom;
938+ const total = Math.ceil(seconds / fr.dt);
939+ const { contribs, used, wanted, missing } = roomContributions(fr, spec, P, { eps: 1e-6 });
940+ if (missing > 0.02 * wanted && attempt === 0) {
941+ fastRoom = null; // the strike moved too far; record afresh
942+ continue;
943+ }
944+ const master = new Float64Array(total);
945+ const batch = 64;
946+ for (let q0 = 0; q0 < contribs.length; q0 += batch) {
947+ if (fastCancel) throw new RoomCancelled();
948+ const part = synthesize(contribs.slice(q0, q0 + batch), fr.room.steps, total);
949+ for (let t = 0; t < total; t++) master[t] += part[t];
950+ prog.value = synthBase + (1 - synthBase) * Math.min(1, (q0 + batch) / contribs.length);
951+ await new Promise((r) => setTimeout(r, 0));
952+ }
953+
954+ const note =
955+ `(fast, ${used} modes)` +
956+ (fr.kernelDecay > 1e-3
957+ ? ' · the room was still ringing where its recording ends; the clip truncates it'
958+ : '');
959+ fastTrace = { trace: Float32Array.from(master), dt: fr.dt, note };
960+ lastTrace = { trace: fastTrace.trace, dt: fr.dt };
961+ drawSpec();
962+ clearError();
963+ showRecording();
964+ stopAudio();
965+ await playTrace(fastTrace.trace, planPlayback(fastTrace.trace.length, fr.dt));
966+ break;
967+ }
968+ } catch (e) {
969+ if (!(e instanceof RoomCancelled)) showError(e, '');
970+ } finally {
971+ fastBusy = false;
972+ prog.hidden = true;
973+ btn.textContent = 'Fast render';
974+ showRecording();
975+ }
976+}
977+
804978 /* -------------------------------------------------------------- wiring --- */
805979
806980 const TABS = ['play', 'build', 'lab'] as const;
807981 for (const tab of TABS) {
808982 el(`tab-${tab}`).addEventListener('click', () => {
983+ // The Play shape knobs and the Build cards are two views of the same
984+ // values; rebuilding the controls on a tab switch keeps them in step.
985+ refreshControls();
809986 for (const t of TABS) {
810987 el(`tab-${t}`).classList.toggle('active', t === tab);
811988 el(`panel-${t}`).classList.toggle('active', t === tab);
@@ -928,13 +1105,19 @@ el('add-membrane').addEventListener('click', () => addObject('membrane'));
9281105
9291106 const runPause = el<HTMLButtonElement>('runpause');
9301107 const setRunning = (r: boolean): void => {
1108+ if (r && fastTrace) {
1109+ // Live stepping records over the microphone again; the synthesized clip
1110+ // stops being "the latest recording" the moment the solver moves.
1111+ fastTrace = null;
1112+ showRecording();
1113+ }
9311114 running = r;
9321115 runPause.textContent = r ? 'Pause' : 'Run';
9331116 frames = 0;
9341117 lastFpsAt = performance.now();
9351118 dirty = true;
9361119 };
937-runPause.addEventListener('click', () => setRunning(!running && !renderingClip));
1120+runPause.addEventListener('click', () => setRunning(!running && !renderingClip && !fastBusy));
9381121
9391122 el('excite').addEventListener('click', excite);
9401123
@@ -944,8 +1127,22 @@ el('spf').addEventListener('change', (e) => {
9441127
9451128 el('renderclip').addEventListener('click', () => void renderClip());
9461129
1130+el('fastrender').addEventListener('click', () => void fastRenderClip());
1131+
9471132 const listen = el<HTMLButtonElement>('listen');
9481133 listen.addEventListener('click', () => {
1134+ if (fastTrace) {
1135+ const { trace, dt } = fastTrace;
1136+ listen.disabled = true;
1137+ stopAudio();
1138+ void playTrace(trace, planPlayback(trace.length, dt))
1139+ .then(() => clearError())
1140+ .catch((e: unknown) => showError(e, ''))
1141+ .finally(() => {
1142+ listen.disabled = false;
1143+ });
1144+ return;
1145+ }
9491146 if (!session) return;
9501147 const s = session;
9511148 if (s.recorder.count === 0) {
@@ -964,6 +1161,15 @@ listen.addEventListener('click', () => {
9641161 });
9651162
9661163 el('download').addEventListener('click', () => {
1164+ if (fastTrace) {
1165+ const url = URL.createObjectURL(traceToWav(fastTrace.trace, fastTrace.dt));
1166+ const a = document.createElement('a');
1167+ a.href = url;
1168+ a.download = 'soundbox.wav';
1169+ a.click();
1170+ setTimeout(() => URL.revokeObjectURL(url), 10_000);
1171+ return;
1172+ }
9671173 if (!session) return;
9681174 const s = session;
9691175 if (s.recorder.count === 0) {
src/studio/bodies.tsmodified+4−4View file
@@ -45,11 +45,11 @@ const boxBody: BodyTemplate = {
4545 { key: 'boxl', label: 'length (m)', value: 0.66, min: 0.3, max: 0.9, step: 0.01,
4646 hint: 'Along x. The default just outreaches the default string.' },
4747 { key: 'boxw', label: 'width (m)', value: 0.16, min: 0.08, max: 0.4, step: 0.005 },
48- { key: 'boxd', label: 'depth (m)', value: 0.08, min: 0.03, max: 0.25, step: 0.005,
48+ { key: 'boxd', label: 'depth (m)', value: 0.08, min: 0.03, max: 0.25, step: 0.005, musical: true,
4949 hint: 'How far the box hangs below its top plate (which sits at z = 0). Deeper means a bigger cavity and a lower air resonance.' },
5050 { key: 'thick', label: 'wall thickness (m)', value: 0.016, min: 0.008, max: 0.05, step: 0.002,
5151 hint: 'Plates thinner than about two grid cells leak; the app warns rather than stops you.' },
52- { key: 'holer', label: 'sound hole radius (m)', value: 0.035, min: 0, max: 0.12, step: 0.005,
52+ { key: 'holer', label: 'sound hole radius (m)', value: 0.035, min: 0, max: 0.12, step: 0.005, musical: true,
5353 hint: '0 seals the box. The hole is what lets the cavity speak.' },
5454 { key: 'holex', label: 'sound hole x (m)', value: -0.05, min: -0.3, max: 0.3, step: 0.01 },
5555 { key: 'absorb', label: 'surface absorption (1/s)', value: 800, min: 0, max: 4000, step: 50,
@@ -78,9 +78,9 @@ const shellBody: BodyTemplate = {
7878 label: 'Drum shell',
7979 blurb: 'A rigid cylinder hanging below z = 0, open at the top — put a membrane over it and it is a drum.',
8080 paramSpecs: [
81- { key: 'shellr', label: 'radius (m)', value: 0.17, min: 0.08, max: 0.3, step: 0.005,
81+ { key: 'shellr', label: 'radius (m)', value: 0.17, min: 0.08, max: 0.3, step: 0.005, musical: true,
8282 hint: 'A hair larger than the head that sits over it.' },
83- { key: 'shelld', label: 'depth (m)', value: 0.3, min: 0.05, max: 0.4, step: 0.01,
83+ { key: 'shelld', label: 'depth (m)', value: 0.3, min: 0.05, max: 0.4, step: 0.01, musical: true,
8484 hint: 'The cavity under the head. Deeper means a lower air resonance.' },
8585 { key: 'thick', label: 'wall thickness (m)', value: 0.016, min: 0.008, max: 0.05, step: 0.002 },
8686 { key: 'closed', label: 'closed bottom (0/1)', value: 1, min: 0, max: 1, step: 1,
src/studio/fast.tsadded+290−0View file
@@ -0,0 +1,290 @@
1+/**
2+ * The fast render path: any generated scene rendered from precomputation
3+ * instead of timestepping, at the same numbers the full solver produces.
4+ *
5+ * The pieces are src/studio/modal.ts (every object in closed form) and
6+ * src/studio/room.ts (every pattern-to-mic impulse response from one
7+ * reciprocal air run); this module strings them together the way
8+ * StudioSession strings the full pipeline: evaluate the scene, record the
9+ * room, synthesize the mic trace. Each object contributes its coupling
10+ * routes as pattern groups — one kernel per mode for the rad route, plus,
11+ * for a string, the single bridge-patch kernel its modes' arriving slopes
12+ * drive. What the recording costs is one air-only run of `irSeconds`, once
13+ * per *geometry* — the body, every object's size and placement, the grid,
14+ * and the microphone. Every model parameter (pitches, decays, brightness,
15+ * excitations, gains) stays free afterwards: re-synthesizing from the same
16+ * room is CPU work measured in seconds. That is a fourth cost tier under
17+ * the app's existing three — params free, scene cheap, structure
18+ * recompiles, and now: geometry or mic edits re-record the room.
19+ *
20+ * The trade, stated plainly: the fast path produces the microphone's trace
21+ * and nothing else — no pressure field to draw — and it is exact only up to
22+ * the membrane mode set kept (eps, fmax; strings and bars keep every mode)
23+ * and the kernel length (irSeconds); a room response still ringing at the
24+ * kernel's end is truncated. `kernelDecay` reports how far down the kernels
25+ * are by then, so the truncation is a number rather than a hope.
26+ */
27+import { makeAirGrid, stableDt, type AirGrid, type ObjectGrid } from '../grids.ts';
28+import { Scene } from '../scene/scene.ts';
29+import { C_AIR, CFL, DOMAIN_X } from '../units.ts';
30+import { generate } from './codegen.ts';
31+import {
32+ accPoles, nextPoles, objectModeSet, radSineTables, synthesize,
33+ type Contribution, type ModeSetOptions,
34+} from './modal.ts';
35+import { recordRoom, type Room, type RoomPatternGroup } from './room.ts';
36+import { suffixed, type ObjectType, type Params, type SceneSpec } from './types.ts';
37+
38+/** Why a scene cannot take the fast path, or null if it can. */
39+export function fastIneligible(spec: SceneSpec, modelSource?: string): string | null {
40+ if (modelSource !== undefined) {
41+ return 'the scene carries a hand-edited model .m; the fast path assumes the generated physics';
42+ }
43+ if (spec.objects.length === 0) {
44+ return 'a scene needs at least one object';
45+ }
46+ return null;
47+}
48+
49+export interface FastRoomOptions {
50+ device: GPUDevice;
51+ spec: SceneSpec;
52+ /** Fully merged parameter values (defaults filled), as specParams gives. */
53+ params: Params;
54+ sceneParams: Params;
55+ /** A hand-edited scene .m, if any; evaluated exactly as the session would. */
56+ sceneSource?: string;
57+ nx: number;
58+ Lx?: number;
59+ mic: { x: number; y: number; z: number };
60+ /** Kernel (room response) length, seconds. */
61+ irSeconds: number;
62+ /** Optional cap on the kernel length, steps (the renderer clips it to the
63+ * trace it is about to make). */
64+ maxKernelSteps?: number;
65+ /** Membrane recording set (line objects always record every mode). */
66+ modeOpts?: ModeSetOptions;
67+ onProgress?: (done: number, total: number) => void;
68+ /** Polled during the recording; returning true abandons it. */
69+ cancelled?: () => boolean;
70+}
71+
72+/** One object's slice of the recorded room. */
73+export interface ObjectRoom {
74+ id: string;
75+ type: ObjectType;
76+ grid: ObjectGrid;
77+ /** The modes with recorded rad kernels, in row order from rowStart
78+ * (n = 0 for a line object's modes). */
79+ modes: { m: number; n: number }[];
80+ rowStart: number;
81+ /** The bridge patch's kernel row, strings only. */
82+ bridgeRow: number | null;
83+}
84+
85+export interface FastRoom {
86+ room: Room;
87+ dt: number;
88+ air: AirGrid;
89+ objects: ObjectRoom[];
90+ /** Peak |kernel| in the recording's last tenth over the overall peak — how
91+ * decayed the room response is where the kernel truncates it. */
92+ kernelDecay: number;
93+}
94+
95+export interface FastRenderOptions extends FastRoomOptions {
96+ /** Trace length, steps. */
97+ steps: number;
98+}
99+
100+export interface FastResult extends FastRoom {
101+ trace: Float32Array;
102+ /** Patterns synthesized / wanted / wanted-but-unrecorded, across objects. */
103+ used: number;
104+ wanted: number;
105+ missing: number;
106+}
107+
108+/** Evaluate the scene and record the room: everything that depends only on
109+ * geometry and the microphone, once. */
110+export async function prepareRoom(opts: FastRoomOptions): Promise<FastRoom> {
111+ const { device, spec, params, sceneParams } = opts;
112+ const bad = fastIneligible(spec);
113+ if (bad) throw new Error(bad);
114+
115+ const air = makeAirGrid(opts.nx, opts.Lx ?? DOMAIN_X);
116+ let dt = stableDt(air.h, C_AIR, CFL);
117+ let generated = generate(spec, air, dt);
118+ const buildScene = (): Scene =>
119+ Scene.create({
120+ air,
121+ source: opts.sceneSource ?? generated.sceneSource,
122+ outputs: generated.sceneOutputs,
123+ consts: generated.sceneConsts,
124+ paramNames: generated.sceneParamNames,
125+ params: sceneParams,
126+ });
127+ const scene = buildScene();
128+ // Reciprocity (see room.ts) rests on the air update's symmetry, which a
129+ // spatially varying sound speed breaks; every generated scene is uniform.
130+ if (scene.cmax - scene.cmin > 1e-6 * scene.cmax) {
131+ throw new Error(
132+ 'the scene has a spatially varying sound speed; the fast path needs the uniform air every generated scene has',
133+ );
134+ }
135+ const dtFinal = stableDt(air.h, Math.max(scene.cmax, C_AIR), CFL);
136+ if (dtFinal < dt) {
137+ dt = dtFinal;
138+ generated = generate(spec, air, dt);
139+ }
140+
141+ const micIndex =
142+ Math.max(0, Math.min(air.nx - 1, Math.round((opts.mic.x + air.Lx / 2) / air.h - 0.5))) +
143+ air.nx *
144+ (Math.max(0, Math.min(air.ny - 1, Math.round((opts.mic.y + air.Ly / 2) / air.h - 0.5))) +
145+ air.ny *
146+ Math.max(0, Math.min(air.nz - 1, Math.round((opts.mic.z + air.Lz / 2) / air.h - 0.5))));
147+
148+ // One rad group per object (a membrane records the mode set the options
149+ // keep; a line object records every interior mode, so no later excitation
150+ // can want one the room lacks), plus a string's bridge-patch group.
151+ const groups: RoomPatternGroup[] = [];
152+ const objects: ObjectRoom[] = [];
153+ let row = 0;
154+ for (const o of spec.objects) {
155+ const grid = generated.grids[o.id];
156+ const { sx, sy, M1, M2 } = radSineTables(grid, air);
157+ const modes =
158+ o.type === 'membrane'
159+ ? objectModeSet(o, grid, dt, params, opts.modeOpts).map((md) => ({ m: md.m, n: md.n }))
160+ : Array.from({ length: grid.n1 - 2 }, (_, i) => ({ m: i + 1, n: 0 }));
161+ groups.push({
162+ label: `rad-${o.id}`,
163+ prof: scene.fields[`prof_${o.id}`],
164+ sx, sxRows: M1, sy, syRows: M2,
165+ pairs: modes.map((md) => ({ xi: md.m - 1, yi: Math.max(0, md.n - 1) })),
166+ });
167+ const or: ObjectRoom = {
168+ id: o.id, type: o.type, grid, modes, rowStart: row, bridgeRow: null,
169+ };
170+ row += modes.length;
171+ if (o.type === 'string') {
172+ groups.push({
173+ label: `brg-${o.id}`,
174+ prof: scene.fields[`bprof_${o.id}`],
175+ sx: new Float32Array(air.nx).fill(1), sxRows: 1,
176+ sy: new Float32Array(air.ny).fill(1), syRows: 1,
177+ pairs: [{ xi: 0, yi: 0 }],
178+ });
179+ or.bridgeRow = row;
180+ row += 1;
181+ }
182+ objects.push(or);
183+ }
184+
185+ const irSteps = Math.max(
186+ 1,
187+ Math.min(opts.maxKernelSteps ?? Infinity, Math.round(opts.irSeconds / dt)),
188+ );
189+ const room = await recordRoom({
190+ device,
191+ air,
192+ fields: { c: scene.fields.c, sig: scene.fields.sig, wall: scene.fields.wall },
193+ dt,
194+ micIndex,
195+ groups,
196+ steps: irSteps,
197+ onProgress: opts.onProgress,
198+ cancelled: opts.cancelled,
199+ });
200+
201+ let peak = 0;
202+ let tailPeak = 0;
203+ const tailFrom = Math.floor(0.9 * irSteps);
204+ for (let q = 0; q < room.rows; q++) {
205+ for (let t = 0; t < irSteps; t++) {
206+ const v = Math.abs(room.kernels[q * irSteps + t]);
207+ if (v > peak) peak = v;
208+ if (t >= tailFrom && v > tailPeak) tailPeak = v;
209+ }
210+ }
211+
212+ return { room, dt, air, objects, kernelDecay: peak > 0 ? tailPeak / peak : 0 };
213+}
214+
215+/**
216+ * Everything one synthesis needs from a room under the current parameters —
217+ * the fourth tier's payoff: pitches, decays, brightness, excitations and
218+ * gains all live here, none of them touch the air, so building this and
219+ * running `synthesize` over it is pure CPU work.
220+ *
221+ * A membrane's room only holds kernels for the modes it was recorded with; a
222+ * new strike can reach outside that set, and those modes are dropped and
223+ * counted in `missing` (a recording made with a small eps keeps that rare —
224+ * line objects record every mode and can never miss).
225+ */
226+export function roomContributions(
227+ fr: FastRoom,
228+ spec: SceneSpec,
229+ params: Params,
230+ modeOpts?: ModeSetOptions,
231+): { contribs: Contribution[]; used: number; wanted: number; missing: number } {
232+ const contribs: Contribution[] = [];
233+ let used = 0;
234+ let wanted = 0;
235+ let missing = 0;
236+ const K = fr.room.steps;
237+ const kernelRow = (r: number): Float32Array => fr.room.kernels.subarray(r * K, (r + 1) * K);
238+
239+ for (const o of spec.objects) {
240+ const or = fr.objects.find((x) => x.id === o.id);
241+ if (!or || or.type !== o.type) {
242+ throw new Error(`the room was recorded for a different scene (object '${o.id}')`);
243+ }
244+ const modes = objectModeSet(o, or.grid, fr.dt, params, modeOpts);
245+ const rowOf = new Map<number, number>();
246+ or.modes.forEach((md, q) => rowOf.set((md.m << 16) | md.n, or.rowStart + q));
247+ const gair = params[suffixed('gair', o.id)] ?? 0;
248+ wanted += modes.length;
249+
250+ const present = modes.filter((md) => {
251+ const has = rowOf.has((md.m << 16) | md.n);
252+ if (!has) missing++;
253+ return has;
254+ });
255+ used += present.length;
256+ for (const md of present) {
257+ contribs.push({
258+ poles: accPoles(md, fr.dt, gair),
259+ kernel: kernelRow(rowOf.get((md.m << 16) | md.n)!),
260+ });
261+ }
262+
263+ if (o.type === 'string' && or.bridgeRow !== null) {
264+ // The bridge: the string's tension pulls the patch with the slope its
265+ // modes arrive at the pinned end with. The generated step's term is
266+ // −(gbridge·cs²/L)·bprof·fb, fb = (un[n−1] − un[n−2])/h,
267+ // and un[n−1] = 0, so mode m contributes +gbridge·cs²/(L·h) times
268+ // sin(mπ(n−2)/(n−1)) times its alpha[k+1] series.
269+ const g = or.grid;
270+ const cs = 2 * g.L1 * (params[suffixed('f0', o.id)] ?? 0);
271+ const gbridge = params[suffixed('gbridge', o.id)] ?? 0;
272+ const gain = (gbridge * cs * cs) / (g.L1 * g.h1);
273+ const poles = present.flatMap((md) =>
274+ nextPoles(md, gain * Math.sin((md.m * Math.PI * (g.n1 - 2)) / (g.n1 - 1))),
275+ );
276+ contribs.push({ poles, kernel: kernelRow(or.bridgeRow) });
277+ }
278+ }
279+ return { contribs, used, wanted, missing };
280+}
281+
282+/** The whole fast path in one call: record the room, synthesize the trace. */
283+export async function fastRender(opts: FastRenderOptions): Promise<FastResult> {
284+ const fr = await prepareRoom({ ...opts, maxKernelSteps: opts.steps });
285+ const { contribs, used, wanted, missing } = roomContributions(
286+ fr, opts.spec, opts.params, opts.modeOpts,
287+ );
288+ const trace = synthesize(contribs, fr.room.steps, opts.steps);
289+ return { ...fr, trace, used, wanted, missing };
290+}
src/studio/modal.tsadded+457−0View file
@@ -0,0 +1,457 @@
1+/**
2+ * Every object in closed form: exact modes of the discrete schemes, and the
3+ * mic trace synthesized from them against precomputed room kernels.
4+ *
5+ * This is the mathematics behind the fast render path (src/studio/room.ts
6+ * records the kernels it consumes; src/studio/fast.ts strings the two
7+ * together). The facts it rests on, verified numerically in
8+ * scripts/modal-poc.ts and test/checks.ts:
9+ *
10+ * - Every object's discrete operator has exact sine eigenmodes. The pinned
11+ * 5-point Laplacian's are product sines; the second difference's are
12+ * sines with lam_m = (4/h²) sin²(mπ/2(n−1)); and the fourth difference
13+ * with antisymmetric ghosts is exactly the square of the pinned second
14+ * difference (the README's own point about it), so it shares those sines
15+ * with eigenvalue lam². The damping terms are built from the same
16+ * operators. Each modal amplitude therefore obeys a scalar two-term
17+ * recursion — the leapfrog restricted to one mode — whose solution is a
18+ * closed-form geometric series. No object is ever timestepped.
19+ * - The coupling into the air is linear and one-way, so every source the
20+ * air sees is a sum of *fixed spatial patterns* times scalar series in
21+ * those modal amplitudes: an object's rad route is one pattern per mode
22+ * (prof .* rad(phi_m)), and a string's bridge route is a single pattern
23+ * (bprof) times the slope its modes arrive at the far end with.
24+ * - The whole system is LTI, so the mic trace is the sum over (pattern,
25+ * series) pairs of the series convolved with that pattern's impulse
26+ * response of the room — the kernel, which depends only on geometry.
27+ *
28+ * Each series being exactly Re(c·z^k) per characteristic root, its
29+ * convolution reduces to a complex one-pole recursion over the kernel: the
30+ * whole synthesis is O(poles × samples) scalar work on the CPU, with the
31+ * closed-form exponential tail falling out for free once the kernel ends.
32+ *
33+ * Everything here matches the compiled f32 schemes to their own rounding
34+ * noise; the approximations are the membrane mode set kept (the caller's
35+ * eps/fmax — line objects have few modes and always keep all) and the
36+ * kernel length.
37+ */
38+import type { AirGrid, ObjectGrid } from '../grids.ts';
39+import type { ObjectSpec, Params } from './types.ts';
40+import { suffixed } from './types.ts';
41+
42+/** One root of a mode's recursion, carrying the displacement series:
43+ * alpha[k] += Re(a · z^k). A real root is the zIm = 0 case. */
44+export interface BasePole {
45+ aRe: number;
46+ aIm: number;
47+ zRe: number;
48+ zIm: number;
49+}
50+
51+/** One term of a source series, ready to filter a kernel: Re(c · z^k). */
52+export interface SeriesPole {
53+ cRe: number;
54+ cIm: number;
55+ zRe: number;
56+ zIm: number;
57+}
58+
59+/** One mode of one object: which sine it is (n = 0 on a line grid), its
60+ * discrete frequency and decay under the current parameters, and its
61+ * displacement series as poles. */
62+export interface ModalMode {
63+ /** Sine indices, 1-based; a line mode has n = 0. */
64+ m: number;
65+ n: number;
66+ /** Eigenvalue of minus the discrete second difference / Laplacian. */
67+ lam: number;
68+ /** Discrete frequency (Hz) and 1/e decay time (s) under these params. */
69+ freq: number;
70+ tau: number;
71+ /** Peak source strength |acc coefficient|, for truncation decisions. */
72+ amp: number;
73+ poles: BasePole[];
74+}
75+
76+/** sin(mπi/(n−1)) for m = 1..n−2, i = 0..n−1, row-major by m. */
77+export function sineTable(n: number): Float64Array {
78+ const t = new Float64Array((n - 2) * n);
79+ for (let m = 1; m <= n - 2; m++)
80+ for (let i = 0; i < n; i++) t[(m - 1) * n + i] = Math.sin((m * Math.PI * i) / (n - 1));
81+ return t;
82+}
83+
84+/* ----------------------------------------------------- the recursions --- */
85+
86+/**
87+ * Solve one mode's recursion in closed form.
88+ *
89+ * Restricted to an eigenvalue, every object's step is
90+ * a[k+1] = cA·a[k] + cB·a[k−1],
91+ * and the excitation is the pair (a[0], a[−1]): a strike is a velocity
92+ * (a[0] = 0, a[−1] = b), a pluck a held displacement (a[0] = a[−1] = b).
93+ * The characteristic roots give a[k] as a geometric series; underdamped
94+ * modes (all the audible ones) as one conjugate pair written Re(a·z^k),
95+ * modes pushed past critical damping near the grid's Nyquist as two real
96+ * roots.
97+ */
98+export function solveBasePoles(
99+ cA: number,
100+ cB: number,
101+ a0: number,
102+ am1: number,
103+): BasePole[] {
104+ const disc = cA * cA + 4 * cB;
105+ if (disc < 0) {
106+ const zRe = cA / 2;
107+ const zIm = Math.sqrt(-disc) / 2;
108+ const zAbs2 = zRe * zRe + zIm * zIm;
109+ // alpha[0] = Re(a), alpha[-1] = Re(a·conj(z))/|z|².
110+ const aRe = a0;
111+ const aIm = zIm !== 0 ? (am1 * zAbs2 - aRe * zRe) / zIm : 0;
112+ return [{ aRe, aIm, zRe, zIm }];
113+ }
114+ const sq = Math.sqrt(Math.max(disc, 1e-300));
115+ const r1 = (cA + sq) / 2;
116+ const r2 = (cA - sq) / 2;
117+ if (r1 === r2 || r1 === 0 || r2 === 0) return [];
118+ // a1 + a2 = a0; a1/r1 + a2/r2 = am1.
119+ const det = 1 / r1 - 1 / r2;
120+ const a1 = (am1 - a0 / r2) / det;
121+ const a2 = a0 - a1;
122+ return [
123+ { aRe: a1, aIm: 0, zRe: r1, zIm: 0 },
124+ { aRe: a2, aIm: 0, zRe: r2, zIm: 0 },
125+ ];
126+}
127+
128+const freqOf = (poles: BasePole[], dt: number): { freq: number; tau: number } => {
129+ if (poles.length === 0) return { freq: 0, tau: 0 };
130+ const p = poles[0];
131+ if (p.zIm !== 0) {
132+ const zAbs = Math.hypot(p.zRe, p.zIm);
133+ return {
134+ freq: Math.abs(Math.atan2(p.zIm, p.zRe)) / (2 * Math.PI * dt),
135+ tau: zAbs < 1 ? -dt / Math.log(zAbs) : Infinity,
136+ };
137+ }
138+ const rMax = Math.max(...poles.map((q) => Math.abs(q.zRe)));
139+ return {
140+ freq: poles.some((q) => q.zRe < 0) ? 1 / (2 * dt) : 0,
141+ tau: rMax > 0 && rMax < 1 ? -dt / Math.log(rMax) : Infinity,
142+ };
143+};
144+
145+/** Complex-multiply a base pole's coefficient by f(z). */
146+const scalePole = (p: BasePole, fRe: number, fIm: number): SeriesPole => ({
147+ cRe: p.aRe * fRe - p.aIm * fIm,
148+ cIm: p.aRe * fIm + p.aIm * fRe,
149+ zRe: p.zRe,
150+ zIm: p.zIm,
151+});
152+
153+/** The acceleration series (un − 2u + um)/dt² of a mode, as poles:
154+ * coefficient scaled by (z − 2 + 1/z)/dt² per root. */
155+export function accPoles(mode: ModalMode, dt: number, gain = 1): SeriesPole[] {
156+ return mode.poles.map((p) => {
157+ const zAbs2 = p.zRe * p.zRe + p.zIm * p.zIm;
158+ const fRe = (gain * (p.zRe - 2 + p.zRe / zAbs2)) / (dt * dt);
159+ const fIm = (gain * (p.zIm - p.zIm / zAbs2)) / (dt * dt);
160+ return scalePole(p, fRe, fIm);
161+ });
162+}
163+
164+/** The series alpha[k+1] (the object's *new* displacement at step k, which
165+ * is what a string's bridge slope is read from): coefficient scaled by z. */
166+export function nextPoles(mode: ModalMode, gain = 1): SeriesPole[] {
167+ return mode.poles.map((p) => scalePole(p, gain * p.zRe, gain * p.zIm));
168+}
169+
170+/* ------------------------------------------------------ the mode sets --- */
171+
172+export interface ModeSetOptions {
173+ /** Membranes only: keep modes whose source strength is at least eps × the
174+ * strongest. 0 keeps everything the strike touches. Line objects have
175+ * few modes and always keep all. */
176+ eps?: number;
177+ /** Membranes only: keep modes below this discrete frequency, Hz. */
178+ fmax?: number;
179+}
180+
181+const ampOf = (poles: BasePole[], dt: number): number => {
182+ let amp = 0;
183+ for (const p of accPoles({ poles } as ModalMode, dt)) amp += Math.hypot(p.cRe, p.cIm);
184+ return amp;
185+};
186+
187+/** Eigenvalue of minus the pinned second difference for sine mode m. */
188+export const lineLam = (g: ObjectGrid, m: number): number => {
189+ const s = Math.sin((m * Math.PI) / (2 * (g.n1 - 1)));
190+ return (4 / (g.h1 * g.h1)) * s * s;
191+};
192+
193+/** Eigenvalue of minus the pinned 5-point Laplacian for sine mode (m, n). */
194+export const membraneLam = (g: ObjectGrid, m: number, n: number): number => {
195+ const s1 = Math.sin((m * Math.PI) / (2 * (g.n1 - 1)));
196+ const s2 = Math.sin((n * Math.PI) / (2 * (g.n2 - 1)));
197+ return (4 / (g.h1 * g.h1)) * s1 * s1 + (4 / (g.h2 * g.h2)) * s2 * s2;
198+};
199+
200+/** Project a line field onto every interior sine mode, normalized. */
201+export function projectLine(g: ObjectGrid, v: Float64Array): Float64Array {
202+ const n = g.n1;
203+ const t = sineTable(n);
204+ const out = new Float64Array(n - 2);
205+ for (let m = 0; m < n - 2; m++) {
206+ let s = 0;
207+ for (let i = 0; i < n; i++) s += t[m * n + i] * v[i];
208+ out[m] = s / ((n - 1) / 2);
209+ }
210+ return out;
211+}
212+
213+/** Project a membrane field onto every interior sine mode: two separable
214+ * contractions, out[(m−1) + (n1−2)(n−1)], normalized. */
215+export function projectModes(g: ObjectGrid, v: Float64Array): Float64Array {
216+ const { n1, n2 } = g;
217+ const M1 = n1 - 2;
218+ const M2 = n2 - 2;
219+ const s1 = sineTable(n1);
220+ const s2 = sineTable(n2);
221+ const tmp = new Float64Array(M1 * n2);
222+ for (let m = 0; m < M1; m++) {
223+ for (let j = 0; j < n2; j++) {
224+ let s = 0;
225+ for (let i = 0; i < n1; i++) s += s1[m * n1 + i] * v[i + n1 * j];
226+ tmp[m + M1 * j] = s;
227+ }
228+ }
229+ const out = new Float64Array(M1 * M2);
230+ const norm = ((n1 - 1) / 2) * ((n2 - 1) / 2);
231+ for (let m = 0; m < M1; m++) {
232+ for (let n = 0; n < M2; n++) {
233+ let s = 0;
234+ for (let j = 0; j < n2; j++) s += s2[n * n2 + j] * tmp[m + M1 * j];
235+ out[m + M1 * n] = s / norm;
236+ }
237+ }
238+ return out;
239+}
240+
241+/** One object's parameter, unsuffixed, defaulting hard on absence. */
242+const par = (params: Params, id: string, key: string): number => {
243+ const v = params[suffixed(key, id)];
244+ if (v === undefined) throw new Error(`fast path: parameter '${suffixed(key, id)}' missing`);
245+ return v;
246+};
247+
248+/**
249+ * The mode set for one object under its current parameters, sorted
250+ * strongest first: everything the synthesis needs. For strings and bars
251+ * every interior mode is kept (there are a few dozen); for membranes
252+ * `opts` bounds the set, which with the kernel length is the fast path's
253+ * only approximation.
254+ *
255+ * The excitations are exactly the generated init's: a string's triangular
256+ * pluck held then released (a[0] = a[−1]), a bar's and membrane's
257+ * mallet-shaped initial velocity (a[0] = 0).
258+ */
259+export function objectModeSet(
260+ obj: ObjectSpec,
261+ g: ObjectGrid,
262+ dt: number,
263+ params: Params,
264+ opts: ModeSetOptions = {},
265+): ModalMode[] {
266+ const id = obj.id;
267+ const t60 = par(params, id, 't60');
268+ const sig1 = par(params, id, 'sig1');
269+ const d = (6.91 / t60) * dt;
270+ /** cA/cB from the stiffness's squared angular frequency at eigenvalue lam. */
271+ const coeffs = (w2: number, lam: number): { cA: number; cB: number } => ({
272+ cA: (2 - dt * dt * w2 - 2 * sig1 * dt * lam) / (1 + d),
273+ cB: (-(1 - d) + 2 * sig1 * dt * lam) / (1 + d),
274+ });
275+
276+ if (obj.type === 'membrane') {
277+ const eps = opts.eps ?? 1e-6;
278+ const fmax = opts.fmax ?? Infinity;
279+ const f0 = par(params, id, 'f0');
280+ const cm = (2 * f0) / Math.sqrt(1 / (g.L1 * g.L1) + 1 / (g.L2 * g.L2));
281+ const um0 = new Float64Array(g.npts);
282+ const sx0 = par(params, id, 'strikex') * g.L1;
283+ const sy0 = par(params, id, 'strikey') * g.L2;
284+ const sr = par(params, id, 'striker');
285+ const vel = par(params, id, 'vel');
286+ for (let k = 0; k < g.npts; k++) {
287+ const dx = g.xs64[k] - sx0;
288+ const dy = g.ys64[k] - sy0;
289+ um0[k] = -(dt * vel) * g.pin[k] * Math.exp(-(dx * dx + dy * dy) / (sr * sr));
290+ }
291+ const b = projectModes(g, um0);
292+ const M1 = g.n1 - 2;
293+ const all: ModalMode[] = [];
294+ for (let n = 1; n <= g.n2 - 2; n++) {
295+ for (let m = 1; m <= M1; m++) {
296+ const bmn = b[m - 1 + M1 * (n - 1)];
297+ if (bmn === 0) continue;
298+ const lam = membraneLam(g, m, n);
299+ const { cA, cB } = coeffs(cm * cm * lam, lam);
300+ const poles = solveBasePoles(cA, cB, 0, bmn);
301+ const { freq, tau } = freqOf(poles, dt);
302+ if (poles.length === 0 || freq > fmax) continue;
303+ all.push({ m, n, lam, freq, tau, amp: ampOf(poles, dt), poles });
304+ }
305+ }
306+ all.sort((a, bb) => bb.amp - a.amp);
307+ const ampMax = all.length ? all[0].amp : 0;
308+ return all.filter((md) => md.amp >= eps * ampMax);
309+ }
310+
311+ // A line object: pluck (string) or strike (bar), then the same recursion
312+ // with the type's own stiffness.
313+ const n = g.n1;
314+ const u0 = new Float64Array(n);
315+ let w2: (lam: number) => number;
316+ if (obj.type === 'string') {
317+ const f0 = par(params, id, 'f0');
318+ const B = par(params, id, 'B');
319+ const cs = 2 * g.L1 * f0;
320+ const kap2 = (B * (cs * g.L1) ** 2) / (Math.PI * Math.PI);
321+ w2 = (lam) => cs * cs * lam + kap2 * lam * lam;
322+ const pp = par(params, id, 'pluckpos') * g.L1;
323+ const amp = par(params, id, 'amp');
324+ for (let i = 0; i < n; i++) {
325+ const a = g.xs64[i] / pp;
326+ const bb = (g.L1 - g.xs64[i]) / (g.L1 - pp);
327+ u0[i] = amp * 0.5 * (a + bb - Math.abs(a - bb));
328+ }
329+ } else {
330+ const f1 = par(params, id, 'f1');
331+ const kap = (2 * f1 * g.L1 * g.L1) / Math.PI;
332+ w2 = (lam) => kap * kap * lam * lam;
333+ const sp = par(params, id, 'strikepos') * g.L1;
334+ const sr = par(params, id, 'striker');
335+ const vel = par(params, id, 'vel');
336+ for (let i = 0; i < n; i++) {
337+ const r = (g.xs64[i] - sp) / sr;
338+ u0[i] = -(dt * vel) * g.pin[i] * Math.exp(-r * r);
339+ }
340+ }
341+ const b = projectLine(g, u0);
342+ const out: ModalMode[] = [];
343+ for (let m = 1; m <= n - 2; m++) {
344+ if (b[m - 1] === 0) continue;
345+ const lam = lineLam(g, m);
346+ const { cA, cB } = coeffs(w2(lam), lam);
347+ // A pluck starts held (a[0] = a[−1]); a strike starts moving (a[0] = 0).
348+ const poles =
349+ obj.type === 'string'
350+ ? solveBasePoles(cA, cB, b[m - 1], b[m - 1])
351+ : solveBasePoles(cA, cB, 0, b[m - 1]);
352+ if (poles.length === 0) continue;
353+ const { freq, tau } = freqOf(poles, dt);
354+ out.push({ m, n: 0, lam, freq, tau, amp: ampOf(poles, dt), poles });
355+ }
356+ return out;
357+}
358+
359+/* -------------------------------------------------------- air tables ---- */
360+
361+/** The bilinearly interpolated sine tables the coupling ops effectively
362+ * apply, per air column: rad of a product of sines factorizes per axis,
363+ * and a line object's rad is the x factor alone. */
364+export function radSineTables(
365+ g: ObjectGrid,
366+ air: AirGrid,
367+): { sx: Float32Array; sy: Float32Array; M1: number; M2: number } {
368+ const { n1, n2, h1, h2, L1, L2 } = g;
369+ const M1 = n1 - 2;
370+ const s1 = sineTable(n1);
371+ const sx = new Float32Array(M1 * air.nx);
372+ for (let ix = 0; ix < air.nx; ix++) {
373+ const xw = -air.Lx / 2 + (ix + 0.5) * air.h;
374+ const u = (xw + L1 / 2) / h1;
375+ const j = Math.floor(u);
376+ if (j < 0 || j >= n1 - 1) continue;
377+ const f = u - j;
378+ for (let m = 0; m < M1; m++)
379+ sx[m * air.nx + ix] = (1 - f) * s1[m * n1 + j] + f * s1[m * n1 + j + 1];
380+ }
381+ if (g.kind === 'line') {
382+ // No second axis: the y factor is 1 everywhere (the tube profile in the
383+ // scene's prof carries the y and z structure).
384+ return { sx, sy: new Float32Array(air.ny).fill(1), M1, M2: 1 };
385+ }
386+ const M2 = n2 - 2;
387+ const s2 = sineTable(n2);
388+ const sy = new Float32Array(M2 * air.ny);
389+ for (let iy = 0; iy < air.ny; iy++) {
390+ const yw = -air.Ly / 2 + (iy + 0.5) * air.h;
391+ const v = (yw + L2 / 2) / h2;
392+ const j = Math.floor(v);
393+ if (j < 0 || j >= n2 - 1) continue;
394+ const f = v - j;
395+ for (let n = 0; n < M2; n++)
396+ sy[n * air.ny + iy] = (1 - f) * s2[n * n2 + j] + f * s2[n * n2 + j + 1];
397+ }
398+ return { sx, sy, M1, M2 };
399+}
400+
401+/* --------------------------------------------------------- synthesis ---- */
402+
403+/** One (series, kernel) pair: the poles' series convolved with the kernel. */
404+export interface Contribution {
405+ poles: SeriesPole[];
406+ /** One kernel row, `kernelSteps` samples. */
407+ kernel: Float32Array;
408+}
409+
410+/**
411+ * Synthesize the mic trace: every contribution's kernel filtered by its
412+ * poles, summed. The trace may be longer than the kernels: a pole keeps
413+ * ringing the room's response down after the kernel ends, which is exactly
414+ * the closed-form tail.
415+ */
416+export function synthesize(
417+ contribs: Contribution[],
418+ kernelSteps: number,
419+ totalSteps: number,
420+): Float32Array {
421+ const y = new Float64Array(totalSteps);
422+ for (const c of contribs) {
423+ const h = c.kernel;
424+ for (const pl of c.poles) {
425+ const { cRe, cIm, zRe, zIm } = pl;
426+ let wRe = 0;
427+ let wIm = 0;
428+ let wPk = 0;
429+ const lim = Math.min(kernelSteps, totalSteps);
430+ for (let t = 0; t < lim; t++) {
431+ const nRe = zRe * wRe - zIm * wIm + h[t];
432+ const nIm = zRe * wIm + zIm * wRe;
433+ wRe = nRe;
434+ wIm = nIm;
435+ const m2 = wRe * wRe + wIm * wIm;
436+ if (m2 > wPk) wPk = m2;
437+ y[t] += cRe * wRe - cIm * wIm;
438+ }
439+ // The closed-form tail: the pole rings the kernel's response down with
440+ // no further input. Cut it 240 dB under its own peak — far below both
441+ // f32 and the .wav — which also keeps the recursion out of subnormal
442+ // floats, where every multiply is slow.
443+ const cut = Math.max(wPk * 1e-24, 1e-280);
444+ for (let t = lim; t < totalSteps; t++) {
445+ const nRe = zRe * wRe - zIm * wIm;
446+ const nIm = zRe * wIm + zIm * wRe;
447+ wRe = nRe;
448+ wIm = nIm;
449+ if (wRe * wRe + wIm * wIm < cut) break;
450+ y[t] += cRe * wRe - cIm * wIm;
451+ }
452+ }
453+ }
454+ const out = new Float32Array(totalSteps);
455+ for (let t = 0; t < totalSteps; t++) out[t] = y[t];
456+ return out;
457+}
src/studio/objects.tsmodified+4−4View file
@@ -57,7 +57,7 @@ const stringTemplate: ObjectTemplate = {
5757 blurb: 'A stiff, damped string, plucked; radiates along its line and drives a bridge patch at its far end.',
5858 sizeSpecs: [
5959 {
60- key: 'L1', label: 'length (m)', value: 0.6, min: 0.3, max: 0.9, step: 0.05,
60+ key: 'L1', label: 'length (m)', value: 0.6, min: 0.3, max: 0.9, step: 0.05, musical: true,
6161 hint: 'The speaking length. Changing it recompiles (the grid is sized from it).',
6262 },
6363 ],
@@ -122,7 +122,7 @@ const barTemplate: ObjectTemplate = {
122122 blurb: 'A stiff bar with no tension — a chime; struck with a mallet-shaped velocity, partials at n²·f1.',
123123 sizeSpecs: [
124124 {
125- key: 'L1', label: 'length (m)', value: 0.3, min: 0.15, max: 0.5, step: 0.05,
125+ key: 'L1', label: 'length (m)', value: 0.3, min: 0.15, max: 0.5, step: 0.05, musical: true,
126126 hint: 'Changing it recompiles (the grid is sized from it).',
127127 },
128128 ],
@@ -179,11 +179,11 @@ const membraneTemplate: ObjectTemplate = {
179179 blurb: 'A 2D membrane pinned at its rim — a drum head; struck, with the inharmonic mode set that makes a drum thud rather than chord.',
180180 sizeSpecs: [
181181 {
182- key: 'L1', label: 'width (m)', value: 0.32, min: 0.15, max: 0.5, step: 0.01,
182+ key: 'L1', label: 'width (m)', value: 0.32, min: 0.15, max: 0.5, step: 0.01, musical: true,
183183 hint: 'Extent along x. Changing it recompiles.',
184184 },
185185 {
186- key: 'L2', label: 'depth (m)', value: 0.32, min: 0.15, max: 0.5, step: 0.01,
186+ key: 'L2', label: 'depth (m)', value: 0.32, min: 0.15, max: 0.5, step: 0.01, musical: true,
187187 hint: 'Extent along y. Changing it recompiles.',
188188 },
189189 ],
src/studio/room.tsadded+354−0View file
@@ -0,0 +1,354 @@
1+/**
2+ * The room recording: every pattern-to-mic impulse response, from one run.
3+ *
4+ * The fast render path (src/studio/modal.ts) needs, per coupling pattern —
5+ * a membrane mode's prof .* rad(phi_mn), a string or bar mode's
6+ * prof .* rad(phi_m), a string's bridge patch bprof — the air's response at
7+ * the microphone to an impulse of that pattern. Recording them one run per
8+ * pattern would cost thousands of runs; discrete reciprocity collapses them
9+ * to one. The air update's transfer matrix is symmetric — the masked
10+ * Laplacian in divergence form has w_i·w_j/h² off-diagonal entries and
11+ * everything else is diagonal — so the response at the mic to pattern g
12+ * equals the projection onto g of the field one impulse *at the mic*
13+ * produces. (This needs the uniform sound speed every generated scene has;
14+ * scripts/modal-poc.ts checks the symmetry to machine precision.) So: run
15+ * the air alone, kicked at the microphone, and after every step project the
16+ * field onto every pattern at once.
17+ *
18+ * The air runs through the same pipeline as everything else — an air-only .m
19+ * compiled by the same compiler against the same lapw op, so the fast path's
20+ * numerics are the real solver's, not a reimplementation. The impulse is the
21+ * init: kicking the update with source s at step 0 produces
22+ * p[1] = dt²·wall·s/(1 + sig·dt), which init writes directly.
23+ *
24+ * The projections ride each step's submission, the way the Recorder's mic
25+ * tap does, organized as *pattern groups*: each group contracts the field
26+ * against one 3D profile over z (the sheet, the tube, the bridge patch),
27+ * then against per-axis row tables sx and sy — rad of a product of sines
28+ * factorizes per axis, a line object's y factor is a row of ones, and a
29+ * fixed pattern like the bridge is a single (ones, ones) pair. Kernels
30+ * accumulate on the GPU in chunks and come back to the CPU between chunks.
31+ */
32+import type { AirGrid } from '../grids.ts';
33+import { GpuModel } from '../mgpu/model.ts';
34+import { makeLapw, OpPlan } from '../mgpu/ops.ts';
35+import { WORKGROUP_SIZE } from '../mgpu/wgsl.ts';
36+
37+/** The air-only model the recording runs: the generated step's air lines,
38+ * with the reciprocal impulse as the initial condition. */
39+const roomModel = (): string => `% Generated by soundbox for the room recording: the air alone, kicked at
40+% the microphone. The step is the generated model's air update, verbatim.
41+function [p, pm] = init(npts, wall, sig, imp, dt)
42+ p = (dt*dt) * ((wall .* imp) ./ (1 + sig*dt));
43+ pm = zeros(npts, 1);
44+end
45+
46+function [pn, po] = step(p, pm, c, sig, wall, dt)
47+ lp = lapw(p, wall);
48+ sd = sig * dt;
49+ pn = wall .* ((2*p - (1 - sd) .* pm + (c*dt).^2 .* lp) ./ (1 + sd));
50+ po = p;
51+end
52+`;
53+
54+/** Thrown when a recording is cancelled through `RoomOptions.cancelled`. */
55+export class RoomCancelled extends Error {
56+ constructor() {
57+ super('room recording cancelled');
58+ }
59+}
60+
61+/** One family of patterns sharing a 3D profile: each pair (xi, yi) is the
62+ * pattern prof(x,y,z) · sx[xi](ix) · sy[yi](iy), and produces one kernel
63+ * row. */
64+export interface RoomPatternGroup {
65+ label: string;
66+ /** The 3D profile, air-sized; contracted over z first. */
67+ prof: Float32Array;
68+ /** Row tables: sx is sxRows × nx, sy is syRows × ny. */
69+ sx: Float32Array;
70+ sxRows: number;
71+ sy: Float32Array;
72+ syRows: number;
73+ /** The (sx row, sy row) of each kernel row this group records. */
74+ pairs: { xi: number; yi: number }[];
75+}
76+
77+export interface RoomOptions {
78+ device: GPUDevice;
79+ air: AirGrid;
80+ /** The scene's evaluated medium. */
81+ fields: { c: Float32Array; sig: Float32Array; wall: Float32Array };
82+ dt: number;
83+ /** Flat air index of the microphone cell. */
84+ micIndex: number;
85+ groups: RoomPatternGroup[];
86+ /** Kernel length, steps. */
87+ steps: number;
88+ onProgress?: (done: number, total: number) => void;
89+ /** Polled between batches; returning true abandons the recording (the
90+ * buffers are torn down and RoomCancelled is thrown). */
91+ cancelled?: () => boolean;
92+}
93+
94+export interface Room {
95+ /** One row of `steps` samples per pattern, row-major: the groups'
96+ * `pairs`, concatenated in group order. */
97+ kernels: Float32Array;
98+ steps: number;
99+ rows: number;
100+}
101+
102+const planeShader = (air: AirGrid): string => `
103+@group(0) @binding(0) var<storage, read_write> plane: array<f32>;
104+@group(0) @binding(1) var<storage, read> field: array<f32>;
105+@group(0) @binding(2) var<storage, read> prof: array<f32>;
106+@compute @workgroup_size(${WORKGROUP_SIZE})
107+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
108+ let c = gid.x;
109+ if (c >= ${air.nx * air.ny}u) { return; }
110+ var s = 0.0;
111+ for (var iz = 0u; iz < ${air.nz}u; iz++) {
112+ let i = c + ${air.nx * air.ny}u * iz;
113+ s += prof[i] * field[i];
114+ }
115+ plane[c] = s;
116+}
117+`;
118+
119+const xcShader = (air: AirGrid, rows: number): string => `
120+@group(0) @binding(0) var<storage, read_write> tmp: array<f32>;
121+@group(0) @binding(1) var<storage, read> plane: array<f32>;
122+@group(0) @binding(2) var<storage, read> sx: array<f32>;
123+@compute @workgroup_size(${WORKGROUP_SIZE})
124+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
125+ let k = gid.x;
126+ if (k >= ${rows * air.ny}u) { return; }
127+ let m = k % ${rows}u;
128+ let iy = k / ${rows}u;
129+ var s = 0.0;
130+ for (var ix = 0u; ix < ${air.nx}u; ix++) {
131+ s += sx[m * ${air.nx}u + ix] * plane[ix + ${air.nx}u * iy];
132+ }
133+ tmp[k] = s;
134+}
135+`;
136+
137+const projShader = (
138+ air: AirGrid,
139+ sxRows: number,
140+ nPairs: number,
141+ rowOffset: number,
142+ totalRows: number,
143+ chunk: number,
144+): string => `
145+@group(0) @binding(0) var<storage, read_write> chunkBuf: array<f32>;
146+@group(0) @binding(1) var<storage, read> tmp: array<f32>;
147+@group(0) @binding(2) var<storage, read> sy: array<f32>;
148+@group(0) @binding(3) var<storage, read> pairs: array<u32>; // xi | yi<<16
149+@group(0) @binding(4) var<storage, read> head: array<u32>;
150+@compute @workgroup_size(${WORKGROUP_SIZE})
151+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
152+ let q = gid.x;
153+ if (q >= ${nPairs}u) { return; }
154+ let mn = pairs[q];
155+ let xi = mn & 0xffffu;
156+ let yi = mn >> 16u;
157+ var s = 0.0;
158+ for (var iy = 0u; iy < ${air.ny}u; iy++) {
159+ s += sy[yi * ${air.ny}u + iy] * tmp[xi + ${sxRows}u * iy];
160+ }
161+ let i = head[0];
162+ if (i < ${chunk}u) {
163+ chunkBuf[i * ${totalRows}u + ${rowOffset}u + q] = s;
164+ }
165+}
166+`;
167+
168+const bumpShader = `
169+@group(0) @binding(0) var<storage, read_write> head: array<u32>;
170+@compute @workgroup_size(1)
171+fn main() {
172+ head[0] = head[0] + 1u;
173+}
174+`;
175+
176+/** Run the reciprocal recording and return every pattern's kernel. */
177+export async function recordRoom(opts: RoomOptions): Promise<Room> {
178+ const { device, air, fields, dt, micIndex, groups, steps } = opts;
179+ const rows = groups.reduce((a, g) => a + g.pairs.length, 0);
180+ if (rows === 0) return { kernels: new Float32Array(0), steps, rows };
181+ for (const g of groups) {
182+ for (const p of g.pairs) {
183+ if (p.xi < 0 || p.xi >= g.sxRows || p.yi < 0 || p.yi >= g.syRows) {
184+ throw new Error(`room group '${g.label}': pair (${p.xi},${p.yi}) outside its tables`);
185+ }
186+ }
187+ }
188+
189+ // The air-only model, through the ordinary pipeline.
190+ const imp = new Float32Array(air.npts);
191+ imp[micIndex] = 1;
192+ const gpu = await GpuModel.create({
193+ device,
194+ ops: new OpPlan(device, [makeLapw('lapw', air)]),
195+ consts: { npts: air.npts, nx: air.nx, ny: air.ny, nz: air.nz },
196+ fields: { c: fields.c, sig: fields.sig, wall: fields.wall, imp },
197+ source: roomModel(),
198+ paramNames: [],
199+ state: [
200+ { name: 'p', count: air.npts },
201+ { name: 'pm', count: air.npts },
202+ ],
203+ });
204+ gpu.setDt(dt);
205+ gpu.setParams({});
206+
207+ // Kernels come back in chunks of whole samples; each sample is one f32 per
208+ // pattern row. Cap a chunk at 64 MB.
209+ const chunk = Math.max(1, Math.min(steps, Math.floor((64 << 20) / (4 * rows))));
210+ const owned: GPUBuffer[] = [];
211+ const mk = (label: string, size: number): GPUBuffer => {
212+ const b = device.createBuffer({
213+ label: `room-${label}`,
214+ size,
215+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST,
216+ });
217+ owned.push(b);
218+ return b;
219+ };
220+ try {
221+ const headBuf = mk('head', 4);
222+ const chunkBuf = mk('chunk', 4 * rows * chunk);
223+ const readback = device.createBuffer({
224+ label: 'room-readback',
225+ size: 4 * rows * chunk,
226+ usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
227+ });
228+ owned.push(readback);
229+ device.queue.writeBuffer(headBuf, 0, new Uint32Array([0]));
230+
231+ const pipe = async (code: string, label: string): Promise<GPUComputePipeline> =>
232+ device.createComputePipelineAsync({
233+ label: `room-${label}`,
234+ layout: 'auto',
235+ compute: {
236+ module: device.createShaderModule({ code, label: `room-${label}` }),
237+ entryPoint: 'main',
238+ },
239+ });
240+ const bind = (p: GPUComputePipeline, buffers: GPUBuffer[]): GPUBindGroup =>
241+ device.createBindGroup({
242+ layout: p.getBindGroupLayout(0),
243+ entries: buffers.map((buffer, i) => ({ binding: i, resource: { buffer } })),
244+ });
245+ const upload = (label: string, data: Float32Array | Uint32Array): GPUBuffer => {
246+ const b = mk(label, data.byteLength);
247+ device.queue.writeBuffer(b, 0, data as Float32Array<ArrayBuffer>);
248+ return b;
249+ };
250+
251+ const pField = gpu.stateBuffer('p')!;
252+ const planePipe = await pipe(planeShader(air), 'plane');
253+ const bumpPipe = await pipe(bumpShader, 'bump');
254+ const bumpBind = bind(bumpPipe, [headBuf]);
255+
256+ interface GroupRun {
257+ xcPipe: GPUComputePipeline;
258+ projPipe: GPUComputePipeline;
259+ planeBind: GPUBindGroup;
260+ xcBind: GPUBindGroup;
261+ projBind: GPUBindGroup;
262+ xcThreads: number;
263+ nPairs: number;
264+ }
265+ const runs: GroupRun[] = [];
266+ let rowOffset = 0;
267+ for (const g of groups) {
268+ const planeBuf = mk(`plane-${g.label}`, 4 * air.nx * air.ny);
269+ const tmpBuf = mk(`tmp-${g.label}`, 4 * g.sxRows * air.ny);
270+ const profBuf = upload(`prof-${g.label}`, g.prof);
271+ const sxBuf = upload(`sx-${g.label}`, g.sx);
272+ const syBuf = upload(`sy-${g.label}`, g.sy);
273+ const pairsBuf = upload(
274+ `pairs-${g.label}`,
275+ new Uint32Array(g.pairs.map((p) => p.xi | (p.yi << 16))),
276+ );
277+ const xcPipe = await pipe(xcShader(air, g.sxRows), `xc-${g.label}`);
278+ const projPipe = await pipe(
279+ projShader(air, g.sxRows, g.pairs.length, rowOffset, rows, chunk),
280+ `proj-${g.label}`,
281+ );
282+ runs.push({
283+ xcPipe,
284+ projPipe,
285+ planeBind: bind(planePipe, [planeBuf, pField, profBuf]),
286+ xcBind: bind(xcPipe, [tmpBuf, planeBuf, sxBuf]),
287+ projBind: bind(projPipe, [chunkBuf, tmpBuf, syBuf, pairsBuf, headBuf]),
288+ xcThreads: g.sxRows * air.ny,
289+ nPairs: g.pairs.length,
290+ });
291+ rowOffset += g.pairs.length;
292+ }
293+
294+ const encodeProjection = (enc: GPUCommandEncoder): void => {
295+ const pass = enc.beginComputePass({ label: 'room-projection' });
296+ const go = (p: GPUComputePipeline, b: GPUBindGroup, threads: number): void => {
297+ pass.setPipeline(p);
298+ pass.setBindGroup(0, b);
299+ pass.dispatchWorkgroups(Math.ceil(threads / WORKGROUP_SIZE));
300+ };
301+ for (const r of runs) {
302+ go(planePipe, r.planeBind, air.nx * air.ny);
303+ go(r.xcPipe, r.xcBind, r.xcThreads);
304+ go(r.projPipe, r.projBind, r.nPairs);
305+ }
306+ go(bumpPipe, bumpBind, 1);
307+ pass.end();
308+ };
309+
310+ // Sample 0 is the projection of the post-init field (the kicked air, one
311+ // step in); every later sample rides its step's own submission.
312+ const kernels = new Float32Array(rows * steps);
313+ gpu.init();
314+ {
315+ const enc = device.createCommandEncoder({ label: 'room-first-sample' });
316+ encodeProjection(enc);
317+ device.queue.submit([enc.finish()]);
318+ }
319+ let got = 0; // samples already on the CPU
320+ let pending = 1; // samples sitting in the chunk buffer
321+ const flush = async (): Promise<void> => {
322+ const enc = device.createCommandEncoder({ label: 'room-flush' });
323+ enc.copyBufferToBuffer(chunkBuf, 0, readback, 0, 4 * rows * pending);
324+ device.queue.submit([enc.finish()]);
325+ await readback.mapAsync(GPUMapMode.READ, 0, 4 * rows * pending);
326+ const data = new Float32Array(readback.getMappedRange(0, 4 * rows * pending));
327+ // Chunk samples are row-minor; kernels are row-major.
328+ for (let s = 0; s < pending; s++)
329+ for (let q = 0; q < rows; q++) kernels[q * steps + got + s] = data[s * rows + q];
330+ readback.unmap();
331+ device.queue.writeBuffer(headBuf, 0, new Uint32Array([0]));
332+ got += pending;
333+ pending = 0;
334+ };
335+ // Steps per submission, as the plain renderer batches; the sync between
336+ // batches keeps the queue from running unboundedly ahead, and is what
337+ // paces the progress reports and the cancellation checks.
338+ const batch = 2048;
339+ while (got + pending < steps) {
340+ if (opts.cancelled?.()) throw new RoomCancelled();
341+ const n = Math.min(batch, chunk - pending, steps - got - pending);
342+ gpu.step(n, encodeProjection);
343+ pending += n;
344+ await device.queue.onSubmittedWorkDone();
345+ opts.onProgress?.(got + pending, steps);
346+ if (pending === chunk) await flush();
347+ }
348+ if (pending > 0) await flush();
349+ return { kernels, steps, rows };
350+ } finally {
351+ for (const b of owned) b.destroy();
352+ gpu.destroy();
353+ }
354+}
test/checks.tsmodified+77−0View file
@@ -13,6 +13,8 @@ import { StudioSession, specParams } from '../src/studio/session.ts';
1313 import { cloneSpec, presetByKey } from '../src/studio/presets.ts';
1414 import { makeSceneFile, parseSceneFile, serializeSceneFile } from '../src/studio/scenefile.ts';
1515 import { generate } from '../src/studio/codegen.ts';
16+import { fastRender, prepareRoom, roomContributions } from '../src/studio/fast.ts';
17+import { synthesize } from '../src/studio/modal.ts';
1618 import { makeAirGrid, stableDt } from '../src/grids.ts';
1719 import { C_AIR, CFL } from '../src/units.ts';
1820 import type { Params, SceneSpec } from '../src/studio/types.ts';
@@ -432,3 +434,78 @@ export async function splitChecks(device: GPUDevice, check: Check, log: Log): Pr
432434 s.destroy();
433435 s2.destroy();
434436 }
437+
438+/** The modal fast path (src/studio/fast.ts): every preset's mic trace,
439+ * synthesized from closed-form modes and one reciprocal room recording per
440+ * scene, must match the timestepped solver to its own f32 rounding noise —
441+ * across all three object types, both string routes, and every body. And
442+ * it must keep matching from the *same* room when pitch, decay and strike
443+ * change, since the room depends only on geometry. */
444+export async function modalChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
445+ const relErr = (a: Float32Array, b: Float32Array, n: number): number => {
446+ let d = 0;
447+ let r = 0;
448+ for (let i = 0; i < n; i++) {
449+ d += (a[i] - b[i]) * (a[i] - b[i]);
450+ r += b[i] * b[i];
451+ }
452+ return Math.sqrt(d / Math.max(r, 1e-300));
453+ };
454+ const steps = 6000;
455+
456+ const run = async (sceneSpec: SceneSpec, mic: { x: number; y: number; z: number }, params: Params): Promise<Float32Array> => {
457+ const s = await makeSession(device, sceneSpec, 32, params);
458+ s.setMic(mic.x, mic.y, mic.z);
459+ s.excite();
460+ s.step(steps);
461+ await s.sync();
462+ const trace = await s.recorder.read();
463+ s.destroy();
464+ return trace;
465+ };
466+
467+ // One scene per object type: the drum (membrane over a shell), the
468+ // dulcimer (string, both routes, over the box), and the kit (membrane and
469+ // bar together — two objects in one room).
470+ for (const key of ['drum', 'dulcimer', 'kit'] as const) {
471+ const preset = presetByKey(key)!;
472+ const sceneSpec = cloneSpec(preset.spec);
473+ const defaults = specParams(sceneSpec);
474+ const fast = await fastRender({
475+ device, spec: sceneSpec, params: defaults.params, sceneParams: defaults.sceneParams,
476+ nx: 32, mic: preset.mic, steps, irSeconds: 1e9, modeOpts: { eps: 1e-9 },
477+ });
478+ const gpu = await run(sceneSpec, preset.mic, {});
479+ const err = relErr(fast.trace, gpu, steps);
480+ log(` ${key}: ${fast.room.rows} kernels, ${fast.used} modes; rel err ${err.toExponential(2)}`);
481+ check(`the modal path reproduces the ${key}`, err < 5e-3, `rel err ${err.toExponential(2)}`);
482+ }
483+
484+ // Retune and restrike from one kept room: pitch, decay, brightness and
485+ // excitation move, the room does not. Membranes are recorded a decade
486+ // deeper than synthesis keeps, so a moved strike stays covered.
487+ {
488+ const preset = presetByKey('drum')!;
489+ const sceneSpec = cloneSpec(preset.spec);
490+ const defaults = specParams(sceneSpec);
491+ const fr = await prepareRoom({
492+ device, spec: sceneSpec, params: defaults.params, sceneParams: defaults.sceneParams,
493+ nx: 32, mic: preset.mic, irSeconds: 1e9, maxKernelSteps: steps, modeOpts: { eps: 1e-12 },
494+ });
495+ const edits: Params = { f0_m1: 150, t60_m1: 0.7, sig1_m1: 0.002, strikex_m1: 0.62, striker_m1: 0.03 };
496+ const paramsB = { ...defaults.params, ...edits };
497+ const { contribs, used, wanted, missing } = roomContributions(fr, sceneSpec, paramsB, { eps: 1e-8 });
498+ const traceB = synthesize(contribs, fr.room.steps, steps);
499+ const gpuB = await run(sceneSpec, preset.mic, edits);
500+ const errB = relErr(traceB, gpuB, steps);
501+ log(` retuned drum: ${used}/${wanted} modes (${missing} unrecorded); rel err ${errB.toExponential(2)}`);
502+ // A stray borderline mode may sit outside the recorded set (a mode near
503+ // a node of the first strike); the error bound is the gate, the count
504+ // just guards against wholesale mismatch.
505+ check(
506+ 'the same room serves a retuned, restruck drum',
507+ errB < 5e-3 && missing <= 0.02 * wanted,
508+ `rel err ${errB.toExponential(2)}, ${missing}/${wanted} modes missing from the room`,
509+ );
510+ }
511+}