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 *
11 * The surface is kept small so that the software WebGPU adapter can run it.
12 */
13import { roundedCube } from '../src/surface.ts';
14import { green, greenDn, selftestRhs, selftestPoints, relError, ALPHA, type CVec } from '../src/problem.ts';
15import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
16import { requestSolverDevice } from '../src/gpu/device.ts';
17import { Operator } from '../src/gpu/operator.ts';
19declare global {
20 interface Window {
21 __RESULTS__?: { passed: number; failed: number; fatal?: string; lines: string[] };
22 }
23}
25const logEl = document.getElementById('log')!;
26const lines: string[] = [];
27let passed = 0;
28let failed = 0;
29const log = (s: string, cls = '') => {
30 lines.push(s);
31 const span = document.createElement('span');
32 span.textContent = s + '\n';
33 if (cls) span.className = cls;
34 logEl.append(span);
35 console.log(s);
36};
37const check = (name: string, val: number, tol: number) => {
38 const ok = val < tol;
39 ok ? passed++ : failed++;
40 log(`${ok ? 'ok ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`, ok ? 'ok' : 'fail');
41};
43async function main() {
44 const t0 = performance.now();
45 if (logEl.textContent === 'starting…') logEl.textContent = '';
46 const gpu = await requestSolverDevice();
47 const device = gpu.device;
48 const lim = device.limits;
49 log(`device: maxStorageBufferBindingSize ${(lim.maxStorageBufferBindingSize / 2 ** 20).toFixed(0)} MB, workgroup storage ${lim.maxComputeWorkgroupStorageSize} B`);
51 const params = new URLSearchParams(location.search);
52 const m = Number(params.get('m') ?? 3);
53 const norder = Number(params.get('norder') ?? 6);
54 const k = Number(params.get('k') ?? 4);
55 const a = 1;
56 const p = Number(params.get('p') ?? 8);
57 const surf = roundedCube(a, p, m, norder);
58 const npts = surf.npts;
59 const npols = surf.npols;
60 log(`surface: p = ${p}, ${surf.npatches} patches of order ${norder}, ${npts} nodes; k = ${k}`);
61 const xs = surf.srcvals.subarray(0, 3 * npts);
62 const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
64 const lam = -0.2;
65 const op = new Operator(gpu, surf, k, ALPHA, lam, 0, { progress: () => {} });
66 const info = await op.setup();
67 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`);
69 // 1. near weights against the reference, for three nodes of one patch
70 {
71 const { pairs, W } = await op.readWallWeights();
72 const tp = 4;
73 const worst = new Float64Array(6);
74 for (const loc of [0, Math.floor(npols / 2), npols - 2]) {
75 const i = tp * npols + loc;
76 const xt = [xs[i], xs[npts + i], xs[2 * npts + i]];
77 const nt = [nrm[i], nrm[npts + i], nrm[2 * npts + i]];
78 for (let pp = pairs.rowPtr[i]; pp < pairs.rowPtr[i + 1]; pp++) {
79 const j = pairs.pairP[pp];
80 const self = j === tp;
81 const I = self ? refSelf(surf, j, surf.uvNodes[i], surf.uvNodes[npts + i], xt, nt, k) : refLeaf(surf, j, xt, nt, k);
82 const Wref = momentsToWeights(surf, I);
83 const errs: string[] = [];
84 for (let mm = 0; mm < 6; mm++) {
85 let num = 0;
86 let den = 0;
87 for (let l = 0; l < 2 * npols; l++) {
88 const g = W[2 * ((pp * 6 + mm) * npols) + l];
89 num += (g - Wref[mm][l]) ** 2;
90 den += Wref[mm][l] ** 2;
91 }
92 const e = Math.sqrt(num / den);
93 worst[mm] = Math.max(worst[mm], e);
94 errs.push(e.toExponential(1));
95 }
96 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];
97 log(` node ${i} patch ${j} (${self ? 'self' : 'near'}, dist/rad ${d.toFixed(2)}): rel err per kernel ${errs.join(' ')}`);
98 }
99 }
100 for (let mm = 0; mm < 6; mm++) check(`near weights, kernel ${mm}, worst relative error`, worst[mm], 1e-4);
101 }
103 // 2. Green's identity through the operator
104 {
105 const xe = [1.7, 0.6, -0.5];
106 const u = green(xs, xe, k);
107 const un = greenDn(xs, nrm, xe, k);
108 const S = await op.applySingle(4, un);
109 const D = await op.applySingle(5, u);
110 const res = new Float64Array(2 * npts);
111 for (let i = 0; i < 2 * npts; i++) res[i] = S[i] - D[i] - u[i] / 2;
112 let num = 0;
113 let den = 0;
114 for (let i = 0; i < npts; i++) {
115 num += (res[2 * i] ** 2 + res[2 * i + 1] ** 2) * surf.wts[i];
116 den += (u[2 * i] ** 2 + u[2 * i + 1] ** 2) * surf.wts[i];
117 }
118 check("Green's identity S[du/dn] - D[u] - u/2 on the walls (relative, weighted)", Math.sqrt(num / den), 1e-4);
119 }
121 // 3. the self-test problem
122 {
123 const xe = [1.7, 0.6, -0.5];
124 const rhs = selftestRhs(xs, nrm, xe, k, lam);
125 const iters: string[] = [];
126 const sol = await op.solve(rhs, 1e-7, 100, (it, r) => iters.push(`${it}:${r.toExponential(1)}`));
127 log(` GMRES: ${sol.niter} iterations, residual ${sol.rres.toExponential(2)}, ${sol.seconds.toFixed(2)} s`);
128 const xt = selftestPoints(30, a);
129 const ut = await op.eval(xt, sol.x, sol.phi1);
130 const uex = green(xt, xe, k);
131 check('self-test: relative error at 30 interior points', relError(ut, uex), 2e-4);
132 const ub = await op.evalOnSurface(sol.x, sol.phi1);
133 const ubex = green(xs, xe, k);
134 check('self-test: relative error on the walls', relError(ub, ubex), 2e-4);
135 }
137 op.destroy();
138 log(`${passed} passed, ${failed} failed, ${((performance.now() - t0) / 1000).toFixed(1)} s`, failed ? 'fail' : 'ok');
139 window.__RESULTS__ = { passed, failed, lines };
140}
142main().catch((e) => {
143 const msg = e instanceof Error ? `${e.message}\n${e.stack ?? ''}` : String(e);
144 log(`FATAL ${msg}`, 'fail');
145 window.__RESULTS__ = { passed, failed: failed + 1, fatal: msg, lines };
146});