2 * The app: problem controls, the solve pipeline (near quadrature, optional
3 * self-test, room problem, energy balance, field on three planes), and the
4 * views. The operator is kept between solves when only the wall admittance
5 * or the source changes, since its near weights depend only on k and the
6 * surface.
7 */
8import { roundedCube, type Surface } from './surface.ts';
9import {
10 ALPHA,
11 energyBalance,
12 green,
13 inside,
14 maskedPoints,
15 relError,
16 roomRhs,
17 selftestPoints,
18 selftestRhs,
19 slicePoints,
20 type CVec,
21} from './problem.ts';
22import { adapterDescription, Gpu } from './gpu/device.ts';
23import { Operator } from './gpu/operator.ts';
24import { colorLimits, DB_MAX, DB_MIN, type Field, type Mode, type SliceField } from './view/field.ts';
25import { SliceView } from './view/slices.ts';
26import { RoomView } from './view/room3d.ts';
27import { drawColorbar, DIVERGING, SEQUENTIAL } from './view/colormaps.ts';
28import { drawConvergence } from './view/convergence.ts';
30const $ = <T extends HTMLElement>(id: string): T => {
31 const el = document.getElementById(id);
32 if (!el) throw new Error(`no #${id}`);
33 return el as T;
34};
36const A = 1.0;
37const EXTERIOR_SOURCE = [1.7, 0.6, -0.5];
39interface Params {
40 k: number;
41 beta: number;
42 m: number;
43 norder: number;
44 p: number;
45 x0: [number, number, number];
46 eps: number;
47 ngrid: number;
48 selftest: boolean;
49}
51const PRESETS: Record<string, Partial<Params>> = {
52 quick: { k: 6, beta: 0.2, m: 2, norder: 6 },
53 small: { k: 10, beta: 0.2, m: 4, norder: 6 },
54 medium: { k: 10, beta: 0.2, m: 6, norder: 8 },
55 large: { k: 40, beta: 0.05, m: 10, norder: 8 },
56};
57const DEFAULTS: Params = { k: 10, beta: 0.2, m: 4, norder: 6, p: 8, x0: [-0.35, 0.2, 0.1], eps: 1e-5, ngrid: 101, selftest: false };
59const ui = {
60 preset: $<HTMLSelectElement>('preset'),
61 k: $<HTMLInputElement>('k'),
62 beta: $<HTMLInputElement>('beta'),
63 m: $<HTMLInputElement>('m'),
64 norder: $<HTMLSelectElement>('norder'),
65 p: $<HTMLSelectElement>('p'),
66 x0: [$<HTMLInputElement>('x0'), $<HTMLInputElement>('y0'), $<HTMLInputElement>('z0')],
67 eps: $<HTMLSelectElement>('eps'),
68 ngrid: $<HTMLSelectElement>('ngrid'),
69 selftest: $<HTMLInputElement>('selftest'),
70 solve: $<HTMLButtonElement>('solve'),
71 cancel: $<HTMLButtonElement>('cancel'),
72 derived: $('derived'),
73 progress: $<HTMLProgressElement>('progress'),
74 status: $('status'),
75 err: $('err'),
76 mode: $<HTMLSelectElement>('mode'),
77 animate: $<HTMLInputElement>('animate'),
78 showWalls: $<HTMLInputElement>('showWalls'),
79 showPlanes: [$<HTMLInputElement>('showXY'), $<HTMLInputElement>('showXZ'), $<HTMLInputElement>('showYZ')],
80 results: $<HTMLTableElement>('results'),
81 conv: $<HTMLCanvasElement>('conv'),
82 empty: $('empty'),
83 cbar: $<HTMLCanvasElement>('cbar'),
84 cbarHi: $('cbarHi'),
85 cbarLo: $('cbarLo'),
86 hover: $('hover'),
87 device: $('device'),
88};
90const fmtBytes = (b: number) => (b >= 2 ** 30 ? `${(b / 2 ** 30).toFixed(2)} GB` : `${(b / 2 ** 20).toFixed(0)} MB`);
91const fmtS = (s: number) => (s < 10 ? `${s.toFixed(2)} s` : `${s.toFixed(1)} s`);
93// ------------------------------------------------------------------ params
95function readParams(): Params {
96 return {
97 k: Number(ui.k.value),
98 beta: Number(ui.beta.value),
99 m: Math.round(Number(ui.m.value)),
100 norder: Number(ui.norder.value),
101 p: Number(ui.p.value),
102 x0: ui.x0.map((e) => Number(e.value)) as [number, number, number],
103 eps: Number(ui.eps.value),
104 ngrid: Number(ui.ngrid.value),
105 selftest: ui.selftest.checked,
106 };
107}
109function writeParams(p: Params) {
110 ui.k.value = String(p.k);
111 ui.beta.value = String(p.beta);
112 ui.m.value = String(p.m);
113 ui.norder.value = String(p.norder);
114 ui.p.value = String(p.p);
115 p.x0.forEach((v, i) => (ui.x0[i].value = String(+v.toFixed(3))));
116 ui.eps.value = String(p.eps);
117 if (![...ui.eps.options].some((o) => o.value === ui.eps.value)) ui.eps.value = '1e-5';
118 ui.ngrid.value = String(p.ngrid);
119 ui.selftest.checked = p.selftest;
120 matchPreset();
121}
123function matchPreset() {
124 const p = readParams();
125 const name = Object.entries(PRESETS).find(([, v]) => v.k === p.k && v.beta === p.beta && v.m === p.m && v.norder === p.norder)?.[0];
126 ui.preset.value = name ?? 'custom';
127}
129function validate(p: Params): string | null {
130 if (!(p.k > 0 && p.k <= 80)) return 'the wavenumber must be between 0 and 80';
131 if (!(p.beta > 0)) return 'the wall admittance must be positive (absorbing walls)';
132 if (!(p.m >= 1 && p.m <= 16)) return 'patches per face edge must be between 1 and 16';
133 if (!inside(p.x0, A, p.p, 0.04)) return 'the source must be inside the room, a little away from the walls';
134 return null;
135}
137// ------------------------------------------------------------------ derived
139let deviceLimit = Infinity;
141function describe(p: Params): { html: string; blocking: boolean } {
142 const npts = 6 * p.m * p.m * (p.norder + 1) ** 2;
143 const wl = (2 * Math.PI) / p.k;
144 const perWl = ((p.norder + 1) * p.m * wl) / (2 * A);
145 const rfac = p.norder <= 6 ? 2.0 : 1.25;
146 // near patches per node, measured on the named cases: about 9 at rfac 2, 5.5 at rfac 1.25, capped by the patch count
147 const perNode = Math.min(6 * p.m * p.m, rfac >= 2 ? 9.5 : 5.8);
148 const wbytes = npts * perNode * 6 * (p.norder + 1) ** 2 * 8;
149 const parts = [
150 `<b>${npts.toLocaleString()}</b> points`,
151 `room <b>${((2 * A) / wl).toFixed(1)}</b> wavelengths wide`,
152 `<span class="${perWl < 5 ? 'warn' : ''}" title="nodes per wavelength along a wall, on average">${perWl.toFixed(1)} nodes per wavelength</span>`,
153 `<span class="${wbytes > deviceLimit ? 'bad' : wbytes > 0.5 * deviceLimit ? 'warn' : ''}" title="storage for the near weights, estimated; this GPU allows ${fmtBytes(deviceLimit)} per buffer">near weights ≈ ${fmtBytes(wbytes)}</span>`,
154 ];
155 let blocking = false;
156 if (wbytes > 1.3 * deviceLimit) {
157 parts.push(`<span class="bad">probably more than this GPU allows (${fmtBytes(deviceLimit)})</span>`);
158 blocking = false; // an estimate only; let the operator decide
159 }
160 if (perWl < 5) parts.push('<span class="warn">under-resolved: expect large errors</span>');
161 const err = validate(p);
162 if (err) {
163 parts.push(`<span class="bad">${err}</span>`);
164 blocking = true;
165 }
166 return { html: parts.join(' · '), blocking };
167}
169function updateDerived() {
170 const d = describe(readParams());
171 ui.derived.innerHTML = d.html;
172 if (!running) ui.solve.disabled = d.blocking;
173}
175// ------------------------------------------------------------------ status
177let running = false;
178let cancelled = false;
179class Cancelled extends Error {}
181function setStatus(text: string, fraction?: number) {
182 ui.status.textContent = text;
183 if (fraction === undefined) ui.progress.removeAttribute('value');
184 else ui.progress.value = fraction;
185}
187function checkCancel() {
188 if (cancelled) throw new Cancelled('cancelled');
189}
191// ------------------------------------------------------------------ results
193type Row = [string, string, string?];
195function showResults(rows: Row[]) {
196 const tb = document.createElement('tbody');
197 for (const [a, b, cls] of rows) {
198 const tr = document.createElement('tr');
199 const td1 = document.createElement('td');
200 td1.textContent = a;
201 const td2 = document.createElement('td');
202 td2.textContent = b;
203 if (cls) td2.className = cls;
204 tr.append(td1, td2);
205 tb.append(tr);
206 }
207 ui.results.replaceChildren(tb);
208}
210const quality = (e: number, tol: number) => (e < 10 * tol ? '' : e < 1e-2 ? 'warn' : 'bad');
212// ------------------------------------------------------------------ main
214async function start() {
215 if (!navigator.gpu) throw new Error('this browser has no WebGPU. Chrome and Edge 113+, Safari 26+, and Firefox 141+ on Windows have it; on Linux, Chrome may need chrome://flags/#enable-unsafe-webgpu.');
216 const adapter = await navigator.gpu.requestAdapter({ powerPreference: 'high-performance' });
217 if (!adapter) throw new Error('WebGPU found no adapter on this machine.');
218 const lim = adapter.limits;
219 const device = await adapter.requestDevice({
220 requiredLimits: { maxStorageBufferBindingSize: lim.maxStorageBufferBindingSize, maxBufferSize: lim.maxBufferSize },
221 });
222 device.lost.then((info) => {
223 if (info.reason !== 'destroyed') {
224 ui.err.textContent = `The GPU device was lost (${info.message}). This can happen when one dispatch runs too long; reload the page to continue.`;
225 }
226 });
227 deviceLimit = Math.min(device.limits.maxStorageBufferBindingSize, device.limits.maxBufferSize);
228 const software = /swiftshader/i.test(adapterDescription(adapter));
229 ui.device.textContent = `WebGPU adapter: ${adapterDescription(adapter)}; largest buffer ${fmtBytes(deviceLimit)}.` + (software ? ' This is a software adapter, so only the quick case is practical.' : '');
230 const gpu = new Gpu(device, adapter);
232 const slices = new SliceView(
233 $('slices'),
234 (x) => {
235 const p = readParams();
236 if (!inside(x, A, p.p, 0.04)) return;
237 x.forEach((v, i) => (ui.x0[i].value = String(+v.toFixed(3))));
238 updateDerived();
239 if (!running) void solve();
240 },
241 (t) => (ui.hover.textContent = t),
242 );
243 const room = new RoomView(device, $<HTMLCanvasElement>('room'));
245 let op: Operator | null = null;
246 let opKey = '';
247 let surf: Surface | null = null;
248 let field: Field | null = null;
249 let conv: { label: string; res: number[]; color: string }[] = [];
250 let lastEps = 1e-5;
252 const redrawConv = () => drawConvergence(ui.conv, conv, lastEps);
254 const updateColorbar = () => {
255 const mode = ui.mode.value as Mode;
256 drawColorbar(ui.cbar, mode === 'pressure' ? DIVERGING : SEQUENTIAL);
257 if (!field) {
258 ui.cbarHi.textContent = ui.cbarLo.textContent = '';
259 } else if (mode === 'pressure') {
260 ui.cbarHi.textContent = `+${field.umax.toPrecision(2)}`;
261 ui.cbarLo.textContent = `−${field.umax.toPrecision(2)}`;
262 } else {
263 ui.cbarHi.textContent = `+${DB_MAX} dB`;
264 ui.cbarLo.textContent = `${DB_MIN} dB`;
265 }
266 };
268 async function solve() {
269 const p = readParams();
270 const bad = validate(p);
271 if (bad) {
272 ui.err.textContent = bad;
273 return;
274 }
275 running = true;
276 cancelled = false;
277 ui.solve.disabled = true;
278 ui.cancel.disabled = false;
279 ui.err.textContent = '';
280 lastEps = p.eps;
281 const rows: Row[] = [];
282 // the previous solution no longer matches the controls, so take it down
283 field = null;
284 slices.clear();
285 room.clear();
286 ui.empty.textContent = 'Solving…';
287 ui.empty.style.display = '';
288 showResults([['solving…', '']]);
289 updateColorbar();
290 conv = [];
291 redrawConv();
292 try {
293 const t0 = performance.now();
294 const key = `${p.k}|${p.m}|${p.norder}|${p.p}`;
295 const lam = -p.beta;
296 let setupS = 0;
297 if (!op || key !== opKey) {
298 op?.destroy();
299 op = null;
300 opKey = '';
301 setStatus('building the surface…');
302 await new Promise((r) => setTimeout(r, 0));
303 surf = roundedCube(A, p.p, p.m, p.norder);
304 const o = new Operator(gpu, surf, p.k, ALPHA, lam, 0, {
305 progress: (stage, f) => {
306 checkCancel();
307 setStatus(`${stage}: ${(100 * f).toFixed(0)}%`, f);
308 },
309 });
310 setStatus('near quadrature…', 0);
311 const info = await o.setup();
312 op = o;
313 opKey = key;
314 setupS = info.setupSeconds;
315 }
316 const o = op!;
317 const s = surf!;
318 o.lam = lam;
319 const info = o.info;
320 rows.push(['points', s.npts.toLocaleString()]);
321 rows.push(['near pairs · quadtree leaves', `${info.npairs.toLocaleString()} · ${info.nleaves.toLocaleString()}`]);
322 rows.push(['near weights', fmtBytes(info.weightsBytes)]);
323 rows.push(['near quadrature (setup)', setupS ? fmtS(setupS) : `reused (${fmtS(info.setupSeconds)} earlier)`]);
324 const npts = s.npts;
325 const xs = s.srcvals.subarray(0, 3 * npts);
326 const nrm = s.srcvals.subarray(9 * npts, 12 * npts);
327 conv = [];
329 if (p.selftest) {
330 const series = { label: 'self-test', res: [] as number[], color: '#8b949e' };
331 conv.push(series);
332 const sol = await o.solve(selftestRhs(xs, nrm, EXTERIOR_SOURCE, p.k, lam), p.eps, 400, (it, r) => {
333 checkCancel();
334 series.res.push(r);
335 redrawConv();
336 setStatus(`self-test GMRES: iteration ${it}, residual ${r.toExponential(1)}`);
337 });
338 setStatus('self-test: evaluating…');
339 const xt = selftestPoints(50, A);
340 const ut = await o.eval(xt, sol.x, sol.phi1, 'self-test field');
341 const ub = await o.evalOnSurface(sol.x, sol.phi1);
342 const e1 = relError(ut, green(xt, EXTERIOR_SOURCE, p.k));
343 const e2 = relError(ub, green(xs, EXTERIOR_SOURCE, p.k));
344 rows.push(['self-test error, interior · walls', `${e1.toExponential(1)} · ${e2.toExponential(1)}`, quality(Math.max(e1, e2), p.eps)]);
345 }
347 const series = { label: 'room problem', res: [] as number[], color: '#2f81f7' };
348 conv.push(series);
349 const sol = await o.solve(roomRhs(xs, nrm, p.x0, p.k, lam), p.eps, 400, (it, r) => {
350 checkCancel();
351 series.res.push(r);
352 redrawConv();
353 setStatus(`GMRES: iteration ${it}, residual ${r.toExponential(1)}`);
354 });
355 rows.push(['GMRES iterations · residual', `${sol.niter} · ${sol.rres.toExponential(1)}`]);
356 rows.push(['solve', `${fmtS(sol.seconds)} (${((1000 * sol.seconds) / sol.niter).toFixed(0)} ms per iteration)`]);
358 // energy balance
359 setStatus('energy balance…');
360 const v0 = await o.eval(Float64Array.from(p.x0), sol.x, sol.phi1, 'field at the source');
361 const vw = await o.evalOnSurface(sol.x, sol.phi1);
362 const gw = green(xs, p.x0, p.k);
363 const uWall = new Float64Array(2 * npts);
364 for (let i = 0; i < 2 * npts; i++) uWall[i] = gw[i] + vw[i];
365 const eb = energyBalance(p.k, lam, [v0[0], v0[1]], s.wts, uWall);
366 rows.push(['power: source · walls', `${eb.pSrc.toFixed(6)} · ${eb.pWall.toFixed(6)}`]);
367 rows.push(['energy mismatch', eb.mismatch.toExponential(1), quality(eb.mismatch, p.eps)]);
369 // the field on the three planes through the source
370 const tf = performance.now();
371 const { grid, slices: sl } = slicePoints(p.ngrid, p.x0, A, p.p);
372 const sfs: SliceField[] = [];
373 for (const [si, slice] of sl.entries()) {
374 const { pts, index } = maskedPoints(slice.pts, slice.mask);
375 const v = await o.eval(pts, sol.x, sol.phi1, `field on the ${slice.name}-plane (${si + 1} of 3)`);
376 const g = green(pts, p.x0, p.k);
377 const u = new Float64Array(2 * slice.mask.length).fill(NaN);
378 for (let q = 0; q < index.length; q++) {
379 u[2 * index[q]] = g[2 * q] + v[2 * q];
380 u[2 * index[q] + 1] = g[2 * q + 1] + v[2 * q + 1];
381 }
382 sfs.push({ slice, u });
383 }
384 const nfield = sl.reduce((acc, x) => acc + x.mask.reduce((a, b) => a + b, 0), 0);
385 rows.push(['field on the planes', `${fmtS((performance.now() - tf) / 1000)} (${nfield.toLocaleString()} points)`]);
386 rows.push(['total', fmtS((performance.now() - t0) / 1000)]);
387 const lims = colorLimits(sfs, p.x0, p.k);
388 field = { k: p.k, beta: p.beta, x0: [...p.x0], a: A, p: p.p, ngrid: p.ngrid, grid, slices: sfs, uWall, ...lims };
389 slices.setField(field);
390 room.setField(s, field);
391 ui.empty.style.display = 'none';
392 updateColorbar();
393 showResults(rows);
394 redrawConv();
395 setStatus(`done in ${fmtS((performance.now() - t0) / 1000)}`, 1);
396 } catch (e) {
397 ui.empty.textContent = 'No solution: press Solve to compute the field.';
398 if (e instanceof Cancelled) {
399 setStatus('cancelled', 0);
400 // a cancelled setup leaves no usable operator
401 if (op && !op.info) {
402 op.destroy();
403 op = null;
404 opKey = '';
405 }
406 } else {
407 console.error(e instanceof Error ? `${e.message}\n${e.stack}` : String(e));
408 ui.err.textContent = e instanceof Error ? e.message : String(e);
409 setStatus('failed', 0);
410 op?.destroy();
411 op = null;
412 opKey = '';
413 }
414 if (rows.length) showResults(rows);
415 } finally {
416 running = false;
417 ui.cancel.disabled = true;
418 updateDerived();
419 }
420 }
422 ui.solve.addEventListener('click', () => void solve());
423 ui.cancel.addEventListener('click', () => {
424 cancelled = true;
425 setStatus('cancelling…');
426 });
427 ui.preset.addEventListener('change', () => {
428 const pr = PRESETS[ui.preset.value];
429 if (pr) writeParams({ ...readParams(), ...pr });
430 updateDerived();
431 });
432 for (const el of [ui.k, ui.beta, ui.m, ui.norder, ui.p, ...ui.x0, ui.eps, ui.ngrid, ui.selftest]) {
433 el.addEventListener('input', () => {
434 matchPreset();
435 updateDerived();
436 });
437 }
438 ui.mode.addEventListener('change', () => {
439 updateColorbar();
440 room.invalidate();
441 });
442 for (const el of [ui.showWalls, ...ui.showPlanes]) el.addEventListener('change', () => room.invalidate());
443 ui.x0.forEach((el) => el.addEventListener('keydown', (e) => e.key === 'Enter' && !running && void solve()));
445 // the frame loop: one period of the pressure every 2 s
446 let phase = 0;
447 let last = performance.now();
448 const frame = (t: number) => {
449 const dt = (t - last) / 1000;
450 last = t;
451 const animating = ui.animate.checked && ui.mode.value === 'pressure' && field !== null;
452 if (animating) phase = (phase + dt * Math.PI) % (2 * Math.PI);
453 const mode = ui.mode.value as Mode;
454 const ph = mode === 'pressure' ? phase : 0;
455 if (field && (animating || slicesDirty)) {
456 slices.draw(mode, ph);
457 slicesDirty = false;
458 }
459 // ?noroom skips the 3D view: in headless Chrome, presenting to a WebGPU canvas breaks later async calls
460 if (!noRoom) room.draw(mode, ph, { walls: ui.showWalls.checked, planes: ui.showPlanes.map((c) => c.checked) as [boolean, boolean, boolean] }, animating);
461 requestAnimationFrame(frame);
462 };
463 let slicesDirty = true;
464 const noRoom = location.search.includes('noroom');
465 ui.mode.addEventListener('change', () => (slicesDirty = true));
466 ui.animate.addEventListener('change', () => (slicesDirty = true));
467 requestAnimationFrame(frame);
469 writeParams({ ...DEFAULTS, x0: [...DEFAULTS.x0] });
470 updateDerived();
471 updateColorbar();
472 redrawConv();
473 setStatus(software ? 'ready (software adapter: choose the quick case)' : 'ready', 0);
474 if (software) {
475 writeParams({ ...readParams(), ...PRESETS.quick });
476 updateDerived();
477 }
478 (window as unknown as { __app__: unknown }).__app__ = { solve, get field() { return field; } };
479}
481start().catch((e) => {
482 console.error(e);
483 ui.err.textContent = e instanceof Error ? e.message : String(e);
484 setStatus('unavailable', 0);
485});
487export type { CVec };