1/**
2 * What the solver is held to.
3 *
4 * These run against the real pipeline — MATLAB source, numbl lowering,
5 * generated WGSL, GPU — so they check the whole chain rather than any one
6 * piece of it. The checks are physical wherever they can be: a wave should
7 * travel at the speed the medium says, a scatterer that matches its
8 * background should not scatter, an absorbing layer should absorb. What is
9 * left over is the discretization's error, and that is what the numbers here
10 * bound.
11 *
12 * Everything is in SI, at the app's own domain size (`DOMAIN`, 10 m) and
13 * background speed (`C_AIR`, 343 m/s) — not a separate "toy" scale — so a
14 * number that appears here means the same thing it would in the app, and the
15 * synthetic scenes below (`uniform`, `closedBox`) share the app's own
16 * absorbing-layer profile rather than inventing their own.
17 *
18 * Nothing here checks what the app looks like. That is for a browser.
19 */
20import { ModelSession } from '../src/mgpu/session.ts';
21import type { MModel } from '../src/mgpu/registry.ts';
22import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';
23import type { MScene } from '../src/scene/registry.ts';
24import { mSceneByKey, defaultSceneParams } from '../src/scene/registry.ts';
25import { C_AIR, DOMAIN } from '../src/units.ts';
27export type Check = (name: string, ok: boolean, detail: string) => void;
28export type Log = (s: string) => void;
30/** A homogeneous medium with the usual absorbing edge — the case every
31 * scattering result is measured against. Same sponge profile as every real
32 * scene (src/tools/sponge.m via 0.2*L, 1700), so its absorption behaviour is
33 * exactly the app's, not a separately tuned stand-in. */
34const uniform: MScene = {
35 key: 'uniform',
36 label: 'Uniform',
37 blurb: 'No scatterer at all.',
38 params: [],
39 source: `
40function [c, sig] = medium(x, y, L, c0)
41 c = c0 + 0*x;
42 sig = sponge(x, y, L, 0.2*L, 1700);
43end
44`,
45};
47/** The same, with nothing absorbing anywhere: a closed box. */
48const closedBox: MScene = {
49 ...uniform,
50 key: 'closed',
51 source: `
52function [c, sig] = medium(x, y, c0)
53 c = c0 + 0*x;
54 sig = 0*x;
55end
56`,
57};
59/** A model whose "step" is one application of a Laplacian stencil, so the
60 * stencil can be measured on a field we chose. */
61const stencilProbe = (op: 'lap2' | 'lap4'): MModel => ({
62 key: `probe-${op}`,
63 label: `probe ${op}`,
64 blurb: '',
65 state: ['p', 'pm', 't'],
66 params: [],
67 order: op === 'lap2' ? 2 : 4,
68 source: `
69function [p, pm, t] = init(npts)
70 p = zeros(npts, 1);
71 pm = zeros(npts, 1);
72 t = zeros(npts, 1);
73end
75function [pn, pold, tn] = step(p, pm, t)
76 pn = ${op}(p);
77 pold = pm;
78 tn = t;
79end
80`,
81});
83const maxAbs = (a: Float32Array): number => {
84 let m = 0;
85 for (const v of a) m = Math.max(m, Math.abs(v));
86 return m;
87};
89/**
90 * The Laplacian stencils, against a field whose Laplacian is known exactly.
91 *
92 * p = sin(kx*x) * sin(ky*y) has lap(p) = -(kx^2 + ky^2) * p. Both stencils
93 * should reproduce that away from the boundary, the 5-point one to O((k*h)^2)
94 * and the 9-point one to O((k*h)^4) — which at the resolution used here is
95 * two orders of magnitude tighter. The margin matters more than either
96 * number: it is what says the fourth-order stencil is actually fourth order
97 * and not a mistyped second-order one. Purely a statement about the discrete
98 * operator, so it does not depend on the medium or the domain size at all.
99 */
100export async function stencilChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
101 const n = 64;
102 const kx = 4 * Math.PI;
103 const ky = 3 * Math.PI;
104 const errs: Record<string, number> = {};
106 for (const op of ['lap2', 'lap4'] as const) {
107 const session = await ModelSession.create({
108 device,
109 model: stencilProbe(op),
110 params: {},
111 scene: uniform,
112 sceneParams: {},
113 n,
114 });
115 try {
116 const { grid } = session;
117 const p = new Float32Array(grid.npts);
118 for (let i = 0; i < grid.npts; i++) {
119 p[i] = Math.sin(kx * grid.x64[i]) * Math.sin(ky * grid.y64[i]);
120 }
121 session.reset();
122 session.gpu.upload('p', p);
123 session.step(1);
124 const got = await session.read('p');
126 // Interior only: the stencil takes the field outside the grid to be
127 // zero, which is not what this analytic field does.
128 const pad = 3;
129 const want = -(kx * kx + ky * ky);
130 let err = 0;
131 for (let iy = pad; iy < n - pad; iy++) {
132 for (let ix = pad; ix < n - pad; ix++) {
133 const i = ix + n * iy;
134 err = Math.max(err, Math.abs(got[i] - want * p[i]));
135 }
136 }
137 errs[op] = err / Math.abs(want);
138 log(` ${op}: max relative error ${errs[op].toExponential(2)} at n = ${n}`);
139 } finally {
140 session.destroy();
141 }
142 }
144 check(
145 'lap2 matches the analytic Laplacian',
146 errs.lap2 < 0.02,
147 `relative error ${errs.lap2.toExponential(2)} (expect ~(k*h)^2/12)`,
148 );
149 check(
150 'lap4 is far more accurate than lap2',
151 errs.lap4 < errs.lap2 / 10,
152 `${errs.lap4.toExponential(2)} vs ${errs.lap2.toExponential(2)}`,
153 );
154}
156/**
157 * A pulse from a point source should be a ring of radius c*(t - t0).
158 *
159 * This is the end-to-end statement that the thing solves the wave equation:
160 * it exercises the source term, the model's own clock, the stencil, and the
161 * timestep the host computed, and it fails if any of them is wrong by a
162 * constant factor. 600 Hz on this grid (256 points over 10 m, h = 3.9 cm) is
163 * about 15 cells per wavelength — comfortably resolved, and short enough that
164 * a clean, narrow pulse is cheap to run.
165 */
166export async function propagationChecks(
167 device: GPUDevice,
168 check: Check,
169 log: Log,
170): Promise<void> {
171 const model = mModelByKey('leapfrog')!;
172 const t0 = 0.005;
173 const session = await ModelSession.create({
174 device,
175 model,
176 params: {
177 ...defaultParams(model),
178 f: 600,
179 tw: 0.0012,
180 t0,
181 cw: 0,
182 point: 1,
183 x0: 0,
184 y0: 0,
185 w: 0.15,
186 },
187 scene: uniform,
188 sceneParams: {},
189 n: 512,
190 L: DOMAIN,
191 });
192 try {
193 session.reset();
194 // Early enough that the front is still clear of the absorbing layer,
195 // which starts at |x| = 0.3*DOMAIN and would pull the peak back towards
196 // the interior.
197 const until = 0.01;
198 const steps = Math.round(until / session.dt);
199 session.step(steps);
200 const p = await session.read('p');
201 const t = session.steps * session.dt;
203 // Where the wavefront is, along +x from the source at the origin. An
204 // energy centroid rather than the bare peak: at 600 Hz the wavelength is
205 // 0.57 m, so the tallest individual fringe of an oscillating pulse can
206 // sit anywhere within half a wavelength of the envelope's true centre
207 // depending on carrier phase, which swamps the grid's own ~1% resolution.
208 // Weighting position by p^2 averages over the fringes instead of picking
209 // whichever one happens to be tallest.
210 const { grid } = session;
211 const cutoff = 0.275 * grid.L; // stays clear of the sponge, which starts at 0.3*L
212 const iy = Math.floor(grid.ny / 2);
213 let weighted = 0;
214 let weight = 0;
215 for (let ix = Math.floor(grid.nx / 2); ix < grid.nx; ix++) {
216 const i = ix + grid.nx * iy;
217 if (grid.x64[i] > cutoff) break;
218 const w2 = p[i] * p[i];
219 weighted += grid.x64[i] * w2;
220 weight += w2;
221 }
222 const best = weighted / weight;
223 const want = C_AIR * (t - t0);
224 const err = Math.abs(best - want) / want;
225 log(` wavefront at r = ${best.toFixed(4)} m, expected ${want.toFixed(4)} m at t = ${(1000 * t).toFixed(3)} ms`);
226 check(
227 'a pulse travels at the medium speed',
228 err < 0.03,
229 `radius off by ${(100 * err).toFixed(1)}%`,
230 );
231 } finally {
232 session.destroy();
233 }
234}
236/**
237 * The absorbing layer should leave next to nothing behind.
238 *
239 * A plane pulse is launched, crosses the grid, and is swallowed. What is
240 * still in the interior long afterwards is what the sponge reflected, and it
241 * is the honest measure of how open the open boundary is.
242 */
243export async function boundaryChecks(
244 device: GPUDevice,
245 check: Check,
246 log: Log,
247): Promise<void> {
248 const model = mModelByKey('leapfrog')!;
249 const session = await ModelSession.create({
250 device,
251 model,
252 params: {
253 ...defaultParams(model),
254 f: 300, tw: 0.006, t0: 0.03, cw: 0, point: 0, x0: -0.25 * DOMAIN,
255 },
256 scene: uniform,
257 sceneParams: {},
258 n: 256,
259 L: DOMAIN,
260 });
261 try {
262 session.reset();
263 const stepsTo = (t: number): number => Math.round(t / session.dt);
264 session.step(stepsTo(0.07));
265 const peak = maxAbs(await session.read('p'));
267 session.step(stepsTo(0.3) - session.steps);
268 const after = await session.read('p');
270 // The interior only: the sponge itself is allowed to hold whatever it is
271 // busy absorbing.
272 const { grid } = session;
273 const interior = 0.3 * grid.L; // exactly where the sponge starts
274 let residual = 0;
275 for (let i = 0; i < grid.npts; i++) {
276 if (Math.abs(grid.x64[i]) < interior && Math.abs(grid.y64[i]) < interior) {
277 residual = Math.max(residual, Math.abs(after[i]));
278 }
279 }
280 const ratio = residual / peak;
281 log(` peak ${peak.toExponential(2)}, interior residual at t = 0.3 s is ${ratio.toExponential(2)} of it`);
282 check(
283 'the absorbing layer reflects little',
284 ratio < 0.02,
285 `residual ${(100 * ratio).toFixed(2)}% of the incident peak`,
286 );
287 } finally {
288 session.destroy();
289 }
290}
292/**
293 * A scatterer whose speed matches the background is not a scatterer.
294 *
295 * Running the disk scene at cin = 1 must reproduce the uniform medium
296 * exactly, which is a strong statement about the whole scene path: the
297 * smoothed interface, the coordinates, the upload. And at cin = 3 there must
298 * be a scattered field worth looking at, or the app would be drawing nothing.
299 * Uses the real `disk` scene and its real domain, so what is checked is
300 * exactly what the app runs.
301 */
302export async function scatteringChecks(
303 device: GPUDevice,
304 check: Check,
305 log: Log,
306): Promise<void> {
307 const model = mModelByKey('leapfrog')!;
308 const disk = mSceneByKey('disk')!;
309 const params = {
310 ...defaultParams(model),
311 f: 400,
312 tw: 0.002,
313 t0: 0.008,
314 cw: 0,
315 point: 0,
316 x0: -0.25 * DOMAIN,
317 };
318 const n = 256;
319 const run = async (scene: MScene, sceneParams: Record<string, number>): Promise<Float32Array> => {
320 const session = await ModelSession.create({
321 device, model, params, scene, sceneParams, n, L: DOMAIN,
322 });
323 try {
324 session.reset();
325 session.step(Math.round(0.022 / session.dt));
326 return await session.read('p');
327 } finally {
328 session.destroy();
329 }
330 };
332 const plain = await run(uniform, {});
333 const matched = await run(disk, { ...defaultSceneParams(disk), cin: 1, absorb: 0 });
334 const hard = await run(disk, { ...defaultSceneParams(disk), cin: 3, absorb: 0 });
336 const peak = maxAbs(plain);
337 let dMatched = 0;
338 let dHard = 0;
339 for (let i = 0; i < plain.length; i++) {
340 dMatched = Math.max(dMatched, Math.abs(matched[i] - plain[i]));
341 dHard = Math.max(dHard, Math.abs(hard[i] - plain[i]));
342 }
343 log(` matched disk differs by ${(dMatched / peak).toExponential(2)}, hard disk by ${(dHard / peak).toFixed(2)}`);
344 check(
345 'a speed-matched disk does not scatter',
346 dMatched / peak < 1e-3,
347 `scattered field ${(dMatched / peak).toExponential(2)} of the incident peak`,
348 );
349 check(
350 'a hard disk scatters strongly',
351 dHard / peak > 0.2,
352 `scattered field ${(100 * dHard / peak).toFixed(0)}% of the incident peak`,
353 );
354}
356/**
357 * The timestep the host picks should be stable, and near the edge of being
358 * unstable — a scheme that is merely stable because it is crawling is not
359 * evidence of anything. Run a closed box (no absorption at all, so nothing
360 * can hide a slow instability) at 95% of the computed limit and watch it
361 * bounce around for a long time.
362 */
363export async function stabilityChecks(
364 device: GPUDevice,
365 check: Check,
366 log: Log,
367): Promise<void> {
368 for (const key of ['leapfrog', 'leapfrog4']) {
369 const model = mModelByKey(key)!;
370 const session = await ModelSession.create({
371 device,
372 model,
373 params: {
374 ...defaultParams(model),
375 f: 300, tw: 0.003, t0: 0.01, cw: 0, point: 1, x0: 0, y0: 0,
376 },
377 scene: closedBox,
378 sceneParams: {},
379 n: 128,
380 L: DOMAIN,
381 cfl: 0.95,
382 });
383 try {
384 session.reset();
385 session.step(Math.round(0.03 / session.dt));
386 const early = maxAbs(await session.read('p'));
387 session.step(Math.round(0.6 / session.dt));
388 const late = maxAbs(await session.read('p'));
389 log(` ${key}: max|p| ${early.toExponential(2)} at t = 30 ms, ${late.toExponential(2)} at t = 630 ms`);
390 check(
391 `${key} is stable at 95% of the CFL limit`,
392 Number.isFinite(late) && late < 5 * early,
393 `max|p| went from ${early.toExponential(2)} to ${late.toExponential(2)} over 600 ms`,
394 );
395 } finally {
396 session.destroy();
397 }
398 }
399}
401/**
402 * The planner's kernel splitting must not change the answer.
403 *
404 * On any device worth running this on, the leapfrog update fits in one
405 * kernel. Squeeze the budget down to two grid fields per kernel — what a
406 * compatibility-mode device would allow — and the same line has to be
407 * evaluated in half a dozen pieces through scratch buffers. The arithmetic is
408 * the same; only the rounding of the intermediates differs, since each piece
409 * is stored as f32 on the way out.
410 */
411export async function splitChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
412 const model = mModelByKey('leapfrog')!;
413 const disk = mSceneByKey('disk')!;
414 const params = {
415 ...defaultParams(model),
416 f: 400, tw: 0.002, t0: 0.008, point: 0, x0: -0.25 * DOMAIN,
417 };
418 const run = async (operandBudget?: number): Promise<{ p: Float32Array; ops: string[] }> => {
419 const session = await ModelSession.create({
420 device,
421 model,
422 params,
423 scene: disk,
424 sceneParams: defaultSceneParams(disk),
425 n: 128,
426 L: DOMAIN,
427 operandBudget,
428 });
429 try {
430 session.reset();
431 session.step(Math.round(0.02 / session.dt));
432 return { p: await session.read('p'), ops: session.describe().step };
433 } finally {
434 session.destroy();
435 }
436 };
438 const whole = await run();
439 const split = await run(2);
440 const peak = maxAbs(whole.p);
441 let diff = 0;
442 for (let i = 0; i < whole.p.length; i++) {
443 diff = Math.max(diff, Math.abs(whole.p[i] - split.p[i]));
444 }
445 log(` ${whole.ops.length} ops whole, ${split.ops.length} split; largest difference ${(diff / peak).toExponential(2)}`);
446 check(
447 'splitting a kernel does not change the answer',
448 diff / peak < 1e-3,
449 `fields differ by ${(diff / peak).toExponential(2)} of the peak`,
450 );
451 check(
452 'a squeezed budget really does split the update',
453 split.ops.length > whole.ops.length,
454 `${split.ops.length} ops vs ${whole.ops.length}`,
455 );
456}
458/**
459 * The microphone records the field, at the point it is pointed at, once per
460 * timestep.
461 *
462 * Checked against the field itself rather than against a description of it:
463 * after N steps the trace must be N samples long, and its last sample must be
464 * exactly — not approximately — the pressure sitting at the probe's grid point,
465 * since both are the same f32 written by the same kernel. That pins the probe
466 * index, the grid layout, and the fact that the recording dispatch really does
467 * run once per step rather than once per submission.
468 */
469export async function microphoneChecks(
470 device: GPUDevice,
471 check: Check,
472 log: Log,
473): Promise<void> {
474 const model = mModelByKey('leapfrog')!;
475 const session = await ModelSession.create({
476 device,
477 model,
478 params: {
479 ...defaultParams(model),
480 f: 600, tw: 0.0012, t0: 0.003, cw: 0, point: 1, x0: 0, y0: 0, w: 0.15,
481 },
482 scene: uniform,
483 sceneParams: {},
484 n: 128,
485 L: DOMAIN,
486 });
487 try {
488 const { grid } = session;
489 const mx = 1.5;
490 const my = 0.5;
491 session.setMic(mx, my);
492 session.reset();
493 // Long enough for the wave to have reached the microphone and moved on.
494 const steps = Math.round(0.01 / session.dt);
495 session.step(steps);
497 const trace = await session.recorder.read();
498 const field = await session.read('p');
499 const ix = Math.round((mx + grid.L / 2) / grid.h - 0.5);
500 const iy = Math.round((my + grid.L / 2) / grid.h - 0.5);
501 const at = field[ix + grid.nx * iy];
503 let peak = 0;
504 for (const v of trace) peak = Math.max(peak, Math.abs(v));
505 log(` ${trace.length} samples in ${steps} steps, peak ${peak.toExponential(2)}, last ${trace[trace.length - 1].toExponential(3)} vs field ${at.toExponential(3)}`);
507 check(
508 'the microphone records one sample per timestep',
509 trace.length === steps,
510 `${trace.length} samples for ${steps} steps`,
511 );
512 check(
513 'the microphone records the field at its own grid point',
514 trace.length > 0 && trace[trace.length - 1] === at,
515 `last sample ${trace[trace.length - 1]} vs field ${at}`,
516 );
517 check(
518 'the microphone hears the wave arrive',
519 peak > 1e-3,
520 `peak |p| at the microphone was ${peak.toExponential(2)}`,
521 );
523 // Moving it must move what it hears: the same run sampled at the origin,
524 // where a point source is loudest, cannot match a point away from it.
525 session.setMic(0, 0);
526 session.reset();
527 session.step(steps);
528 const atSource = await session.recorder.read();
529 let peak2 = 0;
530 for (const v of atSource) peak2 = Math.max(peak2, Math.abs(v));
531 log(` peak at the source ${peak2.toExponential(2)}, at r = ${Math.hypot(mx, my).toFixed(2)} m ${peak.toExponential(2)}`);
532 check(
533 'moving the microphone changes what it hears',
534 peak2 > peak,
535 `${peak2.toExponential(2)} at the source vs ${peak.toExponential(2)} away from it`,
536 );
537 } finally {
538 session.destroy();
539 }
540}
542/**
543 * A room rings, and sealing it makes it ring for longer.
544 *
545 * Measured as energy in the second half of the microphone's trace against the
546 * first — a ratio rather than an envelope, because a small room beats between
547 * its modes and any one window can land in a null. In the open field the same
548 * pulse passes the microphone once and is gone, which is the contrast that
549 * makes the number mean something. Uses the real `room` scene's own default
550 * geometry, at the app's real domain size.
551 *
552 * The source is tuned for this test's own grid rather than reused from the
553 * scene's `suggest`. A wall slower than the background (`cwall` = 0.15) has a
554 * *shorter* wavelength inside itself than the background does at the same
555 * frequency, by that same factor — the app's own docs on this scene call this
556 * out — and an unresolved wall does not behave like a partial reflector, it
557 * behaves like an absorber: at 256 grid points it swallowed the whole pulse
558 * in a handful of bounces regardless of the `absorb` parameter, making every
559 * room in an early version of this check look identically "sealed" no matter
560 * what. So this uses the app's own 512-point grid, where the wall's own
561 * wavelength is resolved to about ten cells at 220 Hz rather than four at the
562 * scene's demo frequency (700 Hz) — and even then the ring is real but not
563 * long: measured, late/early lands around 0.15, well above the open field's
564 * 0.004 but short of the naive "rings for a while" threshold a lossless room
565 * would give. That is the wall's transmission loss actually doing its job,
566 * not a bug — 74% amplitude reflection per bounce (cwall = 0.15 gives
567 * |c-1|/(c+1) = 0.74) empties a small room in a few tens of bounces.
568 */
569export async function roomChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
570 const model = mModelByKey('leapfrog')!;
571 const room = mSceneByKey('room')!;
572 const params = {
573 ...defaultParams(model),
574 point: 1, x0: -0.7, y0: 0.5, w: 0.03, f: 220, tw: 0.002, t0: 0.008, cw: 0,
575 };
576 const roomDefaults = defaultSceneParams(room);
577 const mic = { x: roomDefaults.side * 0.6, y: -roomDefaults.side * 0.5 };
579 const listen = async (scene: MScene, sceneParams: Record<string, number>) => {
580 const session = await ModelSession.create({
581 device, model, params, scene, sceneParams, n: 512, L: DOMAIN,
582 });
583 try {
584 session.setMic(mic.x, mic.y);
585 session.reset();
586 session.step(Math.round(0.15 / session.dt));
587 const trace = await session.recorder.read();
588 const half = Math.floor(trace.length / 2);
589 let early = 0;
590 let late = 0;
591 for (let i = 0; i < half; i++) early += trace[i] * trace[i];
592 for (let i = half; i < trace.length; i++) late += trace[i] * trace[i];
593 return { early, late, ratio: late / Math.max(early, 1e-30) };
594 } finally {
595 session.destroy();
596 }
597 };
599 const open = await listen(uniform, {});
600 const sealed = await listen(room, { ...roomDefaults, gap: 0 });
601 const wide = await listen(room, { ...roomDefaults, gap: 2 * roomDefaults.side });
602 log(` late/early energy — open field ${open.ratio.toExponential(3)}, sealed room ${sealed.ratio.toFixed(3)}, wide door ${wide.ratio.toFixed(3)}`);
604 check(
605 'a pulse in the open field does not come back',
606 open.ratio < 0.02,
607 `late/early energy ${open.ratio.toExponential(2)}`,
608 );
609 check(
610 'a room rings after the pulse has passed',
611 sealed.ratio > 0.05,
612 `late/early energy ${sealed.ratio.toFixed(3)} inside the room, vs ${open.ratio.toExponential(2)} in the open`,
613 );
614 check(
615 'sound leaves through the doorway',
616 sealed.late > 1.2 * wide.late,
617 `late energy ${sealed.late.toExponential(2)} sealed vs ${wide.late.toExponential(2)} with a wide door`,
618 );
619}
621/**
622 * What a step compiles to. A guard on the fusion passes: if one of them stops
623 * firing, the model still gives the right answer, only several times slower,
624 * and nothing else here would notice.
625 */
626export async function planChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
627 const model = mModelByKey('leapfrog')!;
628 const disk = mSceneByKey('disk')!;
629 const session = await ModelSession.create({
630 device,
631 model,
632 params: defaultParams(model),
633 scene: disk,
634 sceneParams: defaultSceneParams(disk),
635 n: 128,
636 L: DOMAIN,
637 });
638 try {
639 const ops = session.describe().step;
640 for (const line of ops) log(` ${line}`);
641 const kernels = ops.filter((o) => o.startsWith('kernel')).length;
642 const stencils = ops.filter((o) => o.startsWith('stencil')).length;
643 check(
644 'the step uses exactly one stencil dispatch',
645 stencils === 1,
646 `${stencils} stencil ops`,
647 );
648 check(
649 'the source term fuses into the update',
650 kernels <= 6,
651 `${kernels} element-wise kernels (5 expected: u, sd, pn, pold, tn)`,
652 );
653 } finally {
654 session.destroy();
655 }
656}