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 */
36import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
37import { requestAcousticDevice } from '../src/device.ts';
38import { StudioSession, specParams } from '../src/studio/session.ts';
39import { presetByKey, cloneSpec } from '../src/studio/presets.ts';
41const argOf = (name: string, dflt: number): number => {
42 const i = process.argv.indexOf(name);
43 return i >= 0 ? Number(process.argv[i + 1]) : dflt;
44};
45const NX = argOf('--nx', 32);
46const STEPS = argOf('--steps', 1200);
48const 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};
53const 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};
62const 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};
68/* ------------------------------------------------------------- GPU run --- */
70try {
71 await installWebGpu();
72} catch (e) {
73 console.error(errMsg(e), '\n', NO_ADAPTER_HINT);
74 process.exit(1);
75}
76const device = await requestAcousticDevice();
77const preset = presetByKey('drum')!;
78const spec = cloneSpec(preset.spec);
79const { params, sceneParams } = specParams(spec);
80const session = await StudioSession.create({ device, spec, params, sceneParams, nx: NX });
82const air = session.air;
83const g = session.grids['m1'];
84const dt = session.dt;
85const P = session.params;
86console.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);
91session.setMic(preset.mic.x, preset.mic.y, preset.mic.z);
92session.excite();
93session.step(STEPS);
94await session.sync();
95const gpuTrace = await session.recorder.read();
96const gpuU = await session.read('u_m1');
98/* ----------------------------------------------- CPU membrane leapfrog --- */
100const { n1, n2, h1, h2, L1, L2 } = g;
101const mN = g.npts;
102const pin = g.pin;
103const f0 = P.f0_m1;
104const t60 = P.t60_m1;
105const sig1 = P.sig1_m1;
106const gair = P.gair_m1;
107const cm = (2 * f0) / Math.sqrt(1 / (L1 * L1) + 1 / (L2 * L2));
108const dmp = (6.91 / t60) * dt;
110// The strike, exactly as the generated init writes it.
111const um0 = new Float64Array(mN);
112for (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}
119const 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};
135// One membrane step: (u, um) -> un, matching the generated step line.
136const 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};
147// Run the membrane alone, keeping every step's acceleration.
148console.log('\nA/B. membrane: CPU leapfrog and closed-form modes');
149let u = new Float64Array(mN);
150let um = Float64Array.from(um0);
151const 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}
165let okAll = report('CPU leapfrog vs GPU u_m1 (f32)', relErr(u, gpuU), 5e-3);
167/* -------------------------------------------------- modal membrane (B) --- */
169// Interior product-sine modes of the pinned 5-point Laplacian.
170const M1 = n1 - 2;
171const M2 = n2 - 2;
172// Strike projection b[m][n] via two small separable contractions.
173const sinTab1 = new Float64Array(M1 * n1); // sin(m*pi*i/(n1-1))
174const sinTab2 = new Float64Array(M2 * n2);
175for (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));
177for (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));
180const 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};
202const b = projectModes(um0); // alpha[-1] per mode; alpha[0] = 0.
204// Per-mode recursion, exactly the membrane step restricted to one mode.
205const 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};
211// Keep modes the strike actually excites.
212let bMax = 0;
213for (const v of b) bMax = Math.max(bMax, Math.abs(v));
214const kept: { m: number; n: number; lam: number; b: number }[] = [];
215for (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)] });
219console.log(` ${kept.length} of ${M1 * M2} modes carry the strike (rel cutoff 1e-10)`);
221// Modal alpha series and accelerations.
222const 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};
237const modal = kept.map((md) => ({ ...md, ...alphaAcc(md.lam, md.b) }));
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}
253/* ------------------------------------------------------- CPU air replica --- */
255console.log('\nC. coupled CPU replica vs the GPU mic trace');
256const { nx, ny, nz, npts, h } = air;
257const F = session.scene.fields;
258const wall = F.wall;
259const sigF = F.sig;
260const cF = F.c;
261const prof = F.prof_m1;
263// The mic index, as setMic computes it.
264const micIx = Math.max(0, Math.min(nx - 1, Math.round((preset.mic.x + air.Lx / 2) / h - 0.5)));
265const micIy = Math.max(0, Math.min(ny - 1, Math.round((preset.mic.y + air.Ly / 2) / h - 0.5)));
266const micIz = Math.max(0, Math.min(nz - 1, Math.round((preset.mic.z + air.Lz / 2) / h - 0.5)));
267const micIdx = micIx + nx * (micIy + ny * micIz);
269// Bilinear sampling weights of the membrane at each air column (x, y) —
270// exactly what rad_m1 bakes into its WGSL.
271interface RadW { k00: number; k10: number; k01: number; k11: number; w00: number; w10: number; w01: number; w11: number }
272const radW: (RadW | null)[] = [];
273for (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}
293const 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};
299// One air step; drv is the already-assembled source field (or null).
300const 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};
325// The coupled run: modal membrane is already verified ≡ leapfrog, so drive
326// the air from the stored leapfrog accelerations.
327const 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}
342okAll = report('CPU coupled trace vs GPU trace (f32)', relErr(cpuTrace, gpuTrace.subarray(0, STEPS)), 5e-3) && okAll;
344/* ---------------------------------------------------- reciprocity (D/E) --- */
346console.log('\nD. discrete reciprocity on the air');
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.
353const sxA = new Float64Array(M1 * nx);
354const syA = new Float64Array(M2 * ny);
355for (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}
364for (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}
374const 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}
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}
436/* ------------------------------------------------- modal synthesis (E) --- */
438console.log('\nE. modal synthesis vs direct traces');
439const yModal = new Float64Array(STEPS);
440for (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}
450okAll = report('modal synthesis vs CPU coupled trace', relErr(yModal, cpuTrace), 1e-6) && okAll;
451okAll = report('modal synthesis vs GPU trace (f32)', relErr(yModal, gpuTrace.subarray(0, STEPS)), 5e-3) && okAll;
453console.log(okAll ? '\nall checks passed' : '\nSOME CHECKS FAILED');
454session.destroy();
455process.exit(okAll ? 0 : 1);