1/**
2 * The checks, against the real pipeline: MATLAB source -> numbl lowering ->
3 * generated WGSL -> GPU. Physics first — does the string sound its pitch,
4 * does it decay on schedule, does the sealed box actually seal — and the
5 * planner's mechanics after.
6 *
7 * Everything runs on deliberately small grids: these are correctness checks,
8 * and a 32-point air grid already carries every code path the 128-point one
9 * does.
10 */
11import { ModelSession } from '../src/mgpu/session.ts';
12import { dulcimerModel, defaultParams, type Params } from '../src/mgpu/registry.ts';
13import { boxScene, defaultSceneParams } from '../src/scene/registry.ts';
14import { planPlayback } from '../src/audio/play.ts';
16type Check = (name: string, ok: boolean, detail: string) => void;
17type Log = (s: string) => void;
19const peakOf = (a: Float32Array): number => {
20 let m = 0;
21 for (const v of a) m = Math.max(m, Math.abs(v));
22 return m;
23};
25async function makeSession(
26 device: GPUDevice,
27 nx: number,
28 params: Params = {},
29 sceneParams: Params = {},
30 operandBudget?: number,
31): Promise<ModelSession> {
32 return ModelSession.create({
33 device,
34 model: dulcimerModel,
35 params: { ...defaultParams(dulcimerModel), ...params },
36 scene: boxScene,
37 sceneParams: { ...defaultSceneParams(boxScene), ...sceneParams },
38 nx,
39 Ls: 0.6,
40 operandBudget,
41 });
42}
44/** The compiled plan: external ops present, pluck in place, mic riding along. */
45export async function planChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
46 const s = await makeSession(device, 32);
47 const { step } = s.describe();
48 for (const l of step) log(` ${l}`);
49 const externals = step.filter((l) => l.startsWith('external'));
50 const wanted = ['dxx(u)', 'dxx(um)', 'dxxxx(u)', 'spread(acc)', 'bridge(un)', 'lapw(p, wall)'];
51 for (const w of wanted) {
52 check(
53 `the step uses ${w}`,
54 externals.some((l) => l.includes(w)),
55 externals.length ? externals.join('; ') : 'no external ops planned',
56 );
57 }
59 s.pluck();
60 const u = await s.read('u');
61 const amp = defaultParams(dulcimerModel).amp;
62 const peak = peakOf(u);
63 check(
64 'the pluck draws the string to its set height',
65 Math.abs(peak - amp) < 0.15 * amp,
66 `peak |u| = ${peak.toExponential(3)}, amp = ${amp}`,
67 );
68 check('the ends are pinned', u[0] === 0 && u[u.length - 1] === 0, `u[0]=${u[0]}, u[end]=${u[u.length - 1]}`);
70 s.step(100);
71 check(
72 'the microphone samples once per timestep',
73 s.recorder.count === 100,
74 `${s.recorder.count} samples after 100 steps`,
75 );
76 const plan = planPlayback(s.recorder.count, s.dt);
77 check(
78 'playback is real time at the solver rate',
79 plan.realTime && Math.abs(plan.rate * s.dt - 1) < 1e-9,
80 `rate ${plan.rate.toFixed(0)} Hz, dt ${s.dt.toExponential(3)}`,
81 );
82 s.destroy();
83}
85/** The string sounds the pitch it is tuned to. */
86export async function pitchChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
87 // An ideal string (no stiffness, no damping to speak of), watched at a
88 // point for a dozen periods; the zero crossings say the frequency.
89 const f0 = 294;
90 const s = await makeSession(device, 32, { f0, B: 0, sig1: 0, t60: 8 });
91 s.pluck();
92 const node = Math.round(s.string.ns * 0.4);
93 const periods = 12;
94 const stepsTotal = Math.round(periods / f0 / s.dt);
95 const chunk = 10;
96 const series: number[] = [];
97 for (let done = 0; done < stepsTotal; done += chunk) {
98 s.step(chunk);
99 const u = await s.read('u');
100 series.push(u[node]);
101 }
102 let crossings = 0;
103 for (let i = 1; i < series.length; i++) {
104 if ((series[i - 1] < 0 && series[i] >= 0) || (series[i - 1] >= 0 && series[i] < 0)) crossings++;
105 }
106 const measured = crossings / 2 / (stepsTotal * s.dt);
107 log(` ${crossings} zero crossings over ${(stepsTotal * s.dt * 1000).toFixed(1)} ms -> ${measured.toFixed(1)} Hz`);
108 check(
109 `the string sounds its fundamental (${f0} Hz)`,
110 Math.abs(measured - f0) < 0.04 * f0,
111 `measured ${measured.toFixed(1)} Hz`,
112 );
114 // d'Alembert: at half a period the string is the mirror of its pluck, so
115 // the displacement at the pluck point flips sign and shrinks to the
116 // triangle's value at the mirrored position.
117 s.pluck();
118 const u0 = await s.read('u');
119 s.step(Math.round(1 / f0 / 2 / s.dt));
120 const u1 = await s.read('u');
121 const at = Math.round(0.22 * (s.string.ns - 1));
122 check(
123 'half a period later the pluck point has swung through zero',
124 u0[at] > 0 && u1[at] < 0,
125 `u ${u0[at].toExponential(2)} -> ${u1[at].toExponential(2)}`,
126 );
127 s.destroy();
128}
130/** The decay knob means what it says. */
131export async function decayChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
132 const t60 = 0.5;
133 const s = await makeSession(device, 32, { t60, B: 0, sig1: 0 });
134 s.pluck();
135 const before = peakOf(await s.read('u'));
136 // Half of t60: amplitude should be down 30 dB, i.e. to about 3.2%.
137 s.step(Math.round(t60 / 2 / s.dt));
138 const after = peakOf(await s.read('u'));
139 const db = 20 * Math.log10(after / before);
140 log(` |u| ${before.toExponential(2)} -> ${after.toExponential(2)} in ${t60 / 2} s (${db.toFixed(1)} dB)`);
141 check(
142 't60 decays the string on schedule',
143 Math.abs(db + 30) < 4,
144 `${db.toFixed(1)} dB over t60/2, want -30`,
145 );
146 s.destroy();
147}
149/** The wall mask is a wall: a sealed box keeps the sound out. */
150export async function wallChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
151 // Bridge drive only (the direct route off), so the source sits above the
152 // plate; a microphone inside a sealed box should hear far less than one
153 // outside beside it.
154 const seconds = 0.06;
155 const run = async (mic: [number, number, number]): Promise<number> => {
156 const s = await makeSession(
157 device,
158 64,
159 { gline: 0, gbridge: 1 },
160 { holer: 0, absorb: 0, boxd: 0.12, thick: 0.02 },
161 );
162 s.setMic(...mic);
163 s.pluck();
164 s.step(Math.round(seconds / s.dt));
165 const trace = await s.recorder.read();
166 s.destroy();
167 return peakOf(trace);
168 };
169 const inside = await run([0, 0, -0.06]);
170 const outside = await run([0, 0, 0.08]);
171 log(` |p| inside the sealed box ${inside.toExponential(2)}, outside ${outside.toExponential(2)}`);
172 check(
173 'a sealed box keeps the sound out',
174 inside < 0.05 * outside,
175 `inside/outside = ${(inside / outside).toExponential(2)}`,
176 );
177}
179/** The sound hole lets the cavity speak: opening it raises what gets in. */
180export async function holeChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
181 const seconds = 0.06;
182 const run = async (holer: number): Promise<number> => {
183 const s = await makeSession(
184 device,
185 64,
186 { gline: 0, gbridge: 1 },
187 { holer, absorb: 0, boxd: 0.12, thick: 0.02 },
188 );
189 s.setMic(0, 0, -0.06);
190 s.pluck();
191 s.step(Math.round(seconds / s.dt));
192 const trace = await s.recorder.read();
193 s.destroy();
194 return peakOf(trace);
195 };
196 const sealed = await run(0);
197 const open = await run(0.05);
198 log(` |p| in the cavity: sealed ${sealed.toExponential(2)}, open ${open.toExponential(2)}`);
199 check(
200 'the sound hole lets the cavity speak',
201 open > 5 * sealed,
202 `open/sealed = ${(open / sealed).toExponential(2)}`,
203 );
204}
206/** Air sound sits on the string's partial comb, not somewhere else. */
207export async function spectrumChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
208 const f0 = 294;
209 const s = await makeSession(device, 64, { f0 });
210 s.setMic(0.12, 0.08, 0.1);
211 s.pluck();
212 const seconds = 0.12;
213 const total = Math.round(seconds / s.dt);
214 for (let done = 0; done < total; done += 2048) {
215 s.step(Math.min(2048, total - done));
216 await s.sync();
217 }
218 const trace = await s.recorder.read();
219 const dft = (f: number): number => {
220 let re = 0;
221 let im = 0;
222 for (let i = 0; i < trace.length; i++) {
223 const w = 2 * Math.PI * f * i * s.dt;
224 re += trace[i] * Math.cos(w);
225 im -= trace[i] * Math.sin(w);
226 }
227 return Math.hypot(re, im) / trace.length;
228 };
229 const comb = (dft(f0) + dft(2 * f0) + dft(3 * f0)) / 3;
230 const off = (dft(1.41 * f0) + dft(2.53 * f0)) / 2;
231 log(` comb ${comb.toExponential(2)}, off-comb ${off.toExponential(2)}`);
232 check(
233 'the microphone hears the string’s partials',
234 comb > 8 * off,
235 `comb/off = ${(comb / off).toFixed(1)}`,
236 );
237 s.destroy();
238}
240/** Taking the body away leaves open air that still carries the string. */
241export async function bareChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
242 const s = await makeSession(device, 32, { gline: 1, gbridge: 0 }, { body: 0 });
243 let minWall = 1;
244 for (const v of s.scene.wall) minWall = Math.min(minWall, v);
245 check(
246 'removing the body leaves pure air',
247 minWall > 0.999,
248 `min wall mask = ${minWall}`,
249 );
250 s.setMic(0.1, 0.05, 0.05);
251 s.pluck();
252 s.step(Math.round(0.03 / s.dt));
253 const trace = await s.recorder.read();
254 const peak = peakOf(trace);
255 log(` bare-string trace peak ${peak.toExponential(2)}`);
256 check('the bare string still sounds', peak > 0, `trace peak ${peak.toExponential(2)}`);
257 s.destroy();
258}
260/** A starved operand budget splits kernels instead of failing. */
261export async function splitChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
262 const s = await makeSession(device, 32, {}, {}, 2);
263 const { step } = s.describe();
264 const parts = step.filter((l) => l.includes('_part')).length;
265 log(` ${parts} split kernels at budget 2`);
266 check('a tight binding budget splits the update', parts > 0, `${parts} split kernels`);
267 s.pluck();
268 s.step(50);
269 const p = await s.read('p');
270 let bad = 0;
271 for (const v of p) if (!Number.isFinite(v)) bad++;
272 check('the split model still runs', bad === 0, `${bad} non-finite values`);
274 // The split result must agree with the unsplit one.
275 const s2 = await makeSession(device, 32);
276 s2.pluck();
277 s2.step(50);
278 const p2 = await s2.read('p');
279 let worst = 0;
280 const scale = Math.max(peakOf(p), 1e-30);
281 for (let i = 0; i < p.length; i++) worst = Math.max(worst, Math.abs(p[i] - p2[i]));
282 check(
283 'split and unsplit agree',
284 worst < 1e-4 * scale,
285 `worst |diff| = ${worst.toExponential(2)} of peak ${scale.toExponential(2)}`,
286 );
287 s.destroy();
288 s2.destroy();
289}