2 * GPU validation, run in a browser (scripts/test-gpu.mjs drives it headlessly):
3 *
4 * 1. near weights against the brute-force double precision reference, for the
5 * self pair and the near pairs of a few nodes;
6 * 2. Green's identity S_k[du/dn] - D_k[u] = u/2 at all wall nodes through the
7 * full operator (far field plus near weights);
8 * 3. the self-test problem: the field of a source outside the room, solved
9 * with the impedance equation and compared with the exact field inside.
10 *
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 11 * The same three checks then run on a coarse cymbal (178 patches), which uses
12 * the interpolated geometry and the exterior equation: the near weights of a face
13 * patch and a rim patch, Green's identity for a source outside the cymbal, and
14 * the exterior self-test with sources spread inside the cymbal (?cymbal=0 skips it,
15 * ?room=0 skips the rounded cube).
16 *
17 * ?bench=small|medium|large instead runs a named cymbal case of the CUDA proof of
18 * concept (same surface, same unknowns) and reports its self-test error, the
19 * scattered field at the listener, the reciprocity check and the timings, for a
20 * direct comparison with run.py --problem cymbal.
21 *
22 * The surfaces are kept small so that the software WebGPU adapter can run them.
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 24import { roundedCube, type Surface } from '../src/surface.ts';
25import { cymbal, CymbalShape, CYMBAL_DEFAULTS } from '../src/cymbal.ts';
26import { selftestSources, selftestPoints as cymbalPoints, sourceSum, LISTENER, C_SOUND } from '../src/cymbalProblem.ts';
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 27import { green, greenDn, selftestRhs, selftestPoints, relError, ALPHA, type CVec } from '../src/problem.ts';
28import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 29import { adapterDescription, requestSolverDevice } from '../src/gpu/device.ts';
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 30import { Operator } from '../src/gpu/operator.ts';
32declare global {
33 interface Window {
34 __RESULTS__?: { passed: number; failed: number; fatal?: string; lines: string[] };
35 }
36}
38const logEl = document.getElementById('log')!;
39const lines: string[] = [];
40let passed = 0;
41let failed = 0;
42const log = (s: string, cls = '') => {
43 lines.push(s);
44 const span = document.createElement('span');
45 span.textContent = s + '\n';
46 if (cls) span.className = cls;
47 logEl.append(span);
48 console.log(s);
49};
50const check = (name: string, val: number, tol: number) => {
51 const ok = val < tol;
52 ok ? passed++ : failed++;
53 log(`${ok ? 'ok ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`, ok ? 'ok' : 'fail');
54};
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 56/** Near weights against the brute-force reference for three nodes of each patch in tps. */
57async function checkWeights(op: Operator, surf: Surface, tps: number[], k: number, label: string) {
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 58 const npts = surf.npts;
59 const npols = surf.npols;
60 const xs = surf.srcvals.subarray(0, 3 * npts);
61 const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 62 const { pairs, W } = await op.readWallWeights();
63 const worst = new Float64Array(6);
64 for (const tp of tps) {
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 65 for (const loc of [0, Math.floor(npols / 2), npols - 2]) {
66 const i = tp * npols + loc;
67 const xt = [xs[i], xs[npts + i], xs[2 * npts + i]];
68 const nt = [nrm[i], nrm[npts + i], nrm[2 * npts + i]];
69 for (let pp = pairs.rowPtr[i]; pp < pairs.rowPtr[i + 1]; pp++) {
70 const j = pairs.pairP[pp];
71 const self = j === tp;
72 const I = self ? refSelf(surf, j, surf.uvNodes[i], surf.uvNodes[npts + i], xt, nt, k) : refLeaf(surf, j, xt, nt, k);
73 const Wref = momentsToWeights(surf, I);
74 const errs: string[] = [];
75 for (let mm = 0; mm < 6; mm++) {
76 let num = 0;
77 let den = 0;
78 for (let l = 0; l < 2 * npols; l++) {
79 const g = W[2 * ((pp * 6 + mm) * npols) + l];
80 num += (g - Wref[mm][l]) ** 2;
81 den += Wref[mm][l] ** 2;
82 }
83 const e = Math.sqrt(num / den);
84 worst[mm] = Math.max(worst[mm], e);
85 errs.push(e.toExponential(1));
86 }
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 87 const kind = ['plain', 'preimage', 'self'][pairs.pairF[8 * pp + 2]];
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 88 const d = Math.hypot(xt[0] - surf.centroids[j], xt[1] - surf.centroids[surf.npatches + j], xt[2] - surf.centroids[2 * surf.npatches + j]) / surf.rads[j];
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 89 log(` node ${i} patch ${j} (${kind}, dist/rad ${d.toFixed(2)}): rel err per kernel ${errs.join(' ')}`);
91 }
92 }
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 93 for (let mm = 0; mm < 6; mm++) check(`${label}near weights, kernel ${mm}, worst relative error`, worst[mm], 1e-4);
94}
96async function main() {
97 const t0 = performance.now();
98 if (logEl.textContent === 'starting…') logEl.textContent = '';
99 const gpu = await requestSolverDevice();
100 const device = gpu.device;
101 const lim = device.limits;
102 log(`device: ${adapterDescription(gpu.adapter ?? null)}; maxStorageBufferBindingSize ${(lim.maxStorageBufferBindingSize / 2 ** 20).toFixed(0)} MB, workgroup storage ${lim.maxComputeWorkgroupStorageSize} B`);
104 const params = new URLSearchParams(location.search);
105 const bench = params.get('bench');
106 if (bench) {
107 await cymbalBench(gpu, bench);
108 log(`done, ${((performance.now() - t0) / 1000).toFixed(1)} s`);
109 window.__RESULTS__ = { passed, failed, lines };
110 return;
111 }
112 if (params.get('room') !== '0') await roomTests(gpu, params);
113 if (params.get('cymbal') !== '0') await cymbalTests(gpu);
114 log(`${passed} passed, ${failed} failed, ${((performance.now() - t0) / 1000).toFixed(1)} s`, failed ? 'fail' : 'ok');
115 window.__RESULTS__ = { passed, failed, lines };
116}
118/** The rounded cube: exact geometry, interior equation. */
119async function roomTests(gpu: Awaited<ReturnType<typeof requestSolverDevice>>, params: URLSearchParams) {
120 const m = Number(params.get('m') ?? 3);
121 const norder = Number(params.get('norder') ?? 6);
122 const k = Number(params.get('k') ?? 4);
123 const a = 1;
124 const p = Number(params.get('p') ?? 8);
125 const surf = roundedCube(a, p, m, norder);
126 const npts = surf.npts;
127 const npols = surf.npols;
128 log(`surface: p = ${p}, ${surf.npatches} patches of order ${norder}, ${npts} nodes; k = ${k}`);
129 const xs = surf.srcvals.subarray(0, 3 * npts);
130 const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
132 const lam = -0.2;
133 const op = new Operator(gpu, surf, k, ALPHA, lam, 0, { progress: () => {} });
134 const info = await op.setup();
135 log(`operator: ${info.npairs} near pairs, ${info.nleaves} leaves, ${info.noversampled} oversampled nodes, weights ${(info.weightsBytes / 2 ** 20).toFixed(1)} MB, setup ${info.setupSeconds.toFixed(2)} s`);
137 // 1. near weights against the reference, for three nodes of one patch
138 await checkWeights(op, surf, [4], k, '');
140 // 2. Green's identity through the operator
141 {
142 const xe = [1.7, 0.6, -0.5];
143 const u = green(xs, xe, k);
144 const un = greenDn(xs, nrm, xe, k);
145 const S = await op.applySingle(4, un);
146 const D = await op.applySingle(5, u);
147 const res = new Float64Array(2 * npts);
148 for (let i = 0; i < 2 * npts; i++) res[i] = S[i] - D[i] - u[i] / 2;
149 let num = 0;
150 let den = 0;
151 for (let i = 0; i < npts; i++) {
152 num += (res[2 * i] ** 2 + res[2 * i + 1] ** 2) * surf.wts[i];
153 den += (u[2 * i] ** 2 + u[2 * i + 1] ** 2) * surf.wts[i];
154 }
155 check("Green's identity S[du/dn] - D[u] - u/2 on the walls (relative, weighted)", Math.sqrt(num / den), 1e-4);
156 }
158 // 3. the self-test problem
159 {
160 const xe = [1.7, 0.6, -0.5];
161 const rhs = selftestRhs(xs, nrm, xe, k, lam);
162 const iters: string[] = [];
163 const sol = await op.solve(rhs, 1e-7, 100, (it, r) => iters.push(`${it}:${r.toExponential(1)}`));
164 log(` GMRES: ${sol.niter} iterations, residual ${sol.rres.toExponential(2)}, ${sol.seconds.toFixed(2)} s`);
165 const xt = selftestPoints(30, a);
166 const ut = await op.eval(xt, sol.x, sol.phi1);
167 const uex = green(xt, xe, k);
168 check('self-test: relative error at 30 interior points', relError(ut, uex), 2e-4);
169 const ub = await op.evalOnSurface(sol.x, sol.phi1);
170 const ubex = green(xs, xe, k);
171 check('self-test: relative error on the walls', relError(ub, ubex), 2e-4);
172 }
174 op.destroy();
175}
177main().catch((e) => {
178 const msg = e instanceof Error ? `${e.message}\n${e.stack ?? ''}` : String(e);
179 log(`FATAL ${msg}`, 'fail');
180 window.__RESULTS__ = { passed, failed: failed + 1, fatal: msg, lines };
181});
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 182
183/** The same checks on a coarse cymbal: interpolated geometry, exterior equation, lambda = 0. */
184async function cymbalTests(gpu: Awaited<ReturnType<typeof requestSolverDevice>>) {
185 const shape = new CymbalShape({ ...CYMBAL_DEFAULTS, tRim: 0.02 });
186 const surf = cymbal(shape, 0.1, 5);
187 const npts = surf.npts;
188 const k = (2 * Math.PI * 1000) / 343;
189 log(`cymbal surface: rim ${shape.par.tRim * 1000} mm, ${surf.npatches} patches of order ${surf.norder}, ${npts} nodes; k = ${k.toFixed(2)}`);
190 const xs = surf.srcvals.subarray(0, 3 * npts);
191 const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
192 const op = new Operator(gpu, surf, k, ALPHA, 0, 1, { progress: () => {} });
193 const info = await op.setup();
194 log(`operator: ${info.npairs} near pairs, ${info.nleaves} leaves, weights ${(info.weightsBytes / 2 ** 20).toFixed(1)} MB, setup ${info.setupSeconds.toFixed(2)} s`);
196 // a patch on the top face near the bell, and the first rim patch
197 const nface = (surf.npatches - surf.layout.nrim) / 2;
198 await checkWeights(op, surf, [surf.layout.mC ** 2 + 2, 2 * nface], k, 'cymbal: ');
200 // Green's identity for a source outside the body: S[du/dn] - D[u] = u/2 inside it
201 {
202 const xe = [0.05, -0.1, 0.3];
203 const u = green(xs, xe, k);
204 const un = greenDn(xs, nrm, xe, k);
205 const S = await op.applySingle(4, un);
206 const D = await op.applySingle(5, u);
207 let num = 0;
208 let den = 0;
209 for (let i = 0; i < npts; i++) {
210 const rr = S[2 * i] - D[2 * i] - u[2 * i] / 2;
211 const ri = S[2 * i + 1] - D[2 * i + 1] - u[2 * i + 1] / 2;
212 num += (rr * rr + ri * ri) * surf.wts[i];
213 den += (u[2 * i] ** 2 + u[2 * i + 1] ** 2) * surf.wts[i];
214 }
215 check("cymbal: Green's identity on the surface (relative, weighted)", Math.sqrt(num / den), 1e-3);
216 }
218 // the exterior self-test: sources spread inside the cymbal
219 {
220 const { x0s, ws } = selftestSources(shape);
221 const rhs = sourceSum(xs, x0s, ws, k, nrm);
222 const sol = await op.solve(rhs, 1e-7, 100);
223 log(` GMRES: ${sol.niter} iterations, residual ${sol.rres.toExponential(2)}, ${sol.seconds.toFixed(2)} s`);
224 const xt = cymbalPoints(shape, 30);
225 const ut = await op.eval(xt, sol.x, sol.phi1);
226 check('cymbal: self-test, relative error at 30 points around the cymbal', relError(ut, sourceSum(xt, x0s, ws, k)), 2e-3);
227 const ub = await op.evalOnSurface(sol.x, sol.phi1);
228 check('cymbal: self-test, relative error on the surface', relError(ub, sourceSum(xs, x0s, ws, k)), 5e-3);
229 }
230 op.destroy();
231}
233/** A named cymbal case of the CUDA proof of concept, for comparison with run.py. */
234async function cymbalBench(gpu: Awaited<ReturnType<typeof requestSolverDevice>>, name: string) {
235 const cases: Record<string, { freq: number; h: number; norder: number }> = {
236 small: { freq: 1000, h: 0.06, norder: 6 },
237 medium: { freq: 3000, h: 0.04, norder: 8 },
238 large: { freq: 8000, h: 0.03, norder: 8 },
239 };
240 const c = cases[name];
241 if (!c) throw new Error(`no case ${name}`);
242 const shape = new CymbalShape();
243 const surf = cymbal(shape, c.h, c.norder);
244 const npts = surf.npts;
245 const k = (2 * Math.PI * c.freq) / C_SOUND;
246 log(`cymbal ${name}: ${surf.npatches} patches of order ${c.norder}, ${npts} points, k = ${k.toFixed(2)}`);
247 const xs = surf.srcvals.subarray(0, 3 * npts);
248 const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
249 const op = new Operator(gpu, surf, k, ALPHA, 0, 1, { progress: () => {} });
250 const info = await op.setup();
251 log(`setup (near quadrature) ${info.setupSeconds.toFixed(2)} s: ${info.npairs} near pairs, ${info.nleaves} leaves, weights ${(info.weightsBytes / 2 ** 20).toFixed(0)} MB`);
252 const eps = 1e-6;
253 {
254 const { x0s, ws } = selftestSources(shape);
255 const sol = await op.solve(sourceSum(xs, x0s, ws, k, nrm), eps, 400);
256 const xt = cymbalPoints(shape, 50);
257 const e1 = relError(await op.eval(xt, sol.x, sol.phi1), sourceSum(xt, x0s, ws, k));
258 const e2 = relError(await op.evalOnSurface(sol.x, sol.phi1), sourceSum(xs, x0s, ws, k));
259 log(`self-test: ${sol.niter} iterations, ${sol.seconds.toFixed(2)} s; relative error ${e1.toExponential(2)} off the surface, ${e2.toExponential(2)} on it`);
260 }
261 const B = [-0.15, 0.05, -0.12];
262 const rhs = (x0: number[]) => {
263 const g = greenDn(xs, nrm, x0, k);
264 for (let i = 0; i < g.length; i++) g[i] = -g[i];
265 return g;
266 };
267 const solL = await op.solve(rhs(LISTENER), eps, 400);
268 log(`main problem: ${solL.niter} iterations, residual ${solL.rres.toExponential(2)}, solve ${solL.seconds.toFixed(2)} s`);
269 const vL = await op.eval(Float64Array.from([LISTENER[0], B[0], LISTENER[1], B[1], LISTENER[2], B[2]]), solL.x, solL.phi1);
270 log(`scattered field at the listener ${vL[0].toExponential(6)} ${vL[1] < 0 ? '-' : '+'} ${Math.abs(vL[1]).toExponential(6)}i`);
271 const solB = await op.solve(rhs(B), eps, 400);
272 const vB = await op.eval(Float64Array.from(LISTENER), solB.x, solB.phi1);
273 const rec = Math.hypot(vL[2] - vB[0], vL[3] - vB[1]) / Math.hypot(vL[2], vL[3]);
274 log(`reciprocity: at B due to the listener ${vL[2].toExponential(6)} ${vL[3] < 0 ? '-' : '+'} ${Math.abs(vL[3]).toExponential(6)}i, relative difference ${rec.toExponential(2)}`);
275 op.destroy();
276}