WebGPU port of the fmm3dbie impedance solver as an interactive web app
27 changed files+6430−0
.github/workflows/deploy.ymladded+38−0View file
@@ -0,0 +1,38 @@
1+name: deploy
2+on:
3+ push:
4+ branches: [main]
5+ workflow_dispatch:
6+
7+permissions:
8+ contents: read
9+ pages: write
10+ id-token: write
11+
12+concurrency:
13+ group: pages
14+ cancel-in-progress: true
15+
16+jobs:
17+ build-deploy:
18+ runs-on: ubuntu-latest
19+ environment:
20+ name: github-pages
21+ url: ${{ steps.deployment.outputs.page_url }}
22+ steps:
23+ - uses: actions/checkout@v4
24+ - uses: actions/setup-node@v4
25+ with:
26+ node-version: 24
27+ cache: npm
28+ - run: npm ci
29+ # `build` runs tsc --noEmit first, so this is the type check too. The
30+ # GPU tests are not run here: they need a WebGPU-capable headless Chrome.
31+ # Run `npm run test:gpu` and `node scripts/smoke.mjs` locally before pushing.
32+ - run: npm run build
33+ - uses: actions/configure-pages@v5
34+ - uses: actions/upload-pages-artifact@v3
35+ with:
36+ path: dist
37+ - id: deployment
38+ uses: actions/deploy-pages@v4
.gitignoreadded+3−0View file
@@ -0,0 +1,3 @@
1+node_modules/
2+dist/
3+tmp/
README.mdadded+43−0View file
@@ -0,0 +1,43 @@
1+# fmm3dbie-webgpu-poc
2+
3+A limitation of the [CUDA proof of concept](https://vault1.magland.org/proof-of-concept/fmm3dbie-gpu-poc) for fmm3dbie's interior Helmholtz impedance solver is that it needs an NVIDIA GPU, CUDA and a Python environment, which makes it hard to try and hard to share. Here we present a port of the same solver to WebGPU, as a web page. It runs on whatever GPU the browser can reach, including integrated graphics, and solves the same showcase problem interactively.
4+
5+**Live page:** https://fmm3dbie-webgpu-poc--proof-of-concept.vault1-sites.magland.org/
6+
7+The problem is a point source at a fixed frequency in a closed room with partially absorbing walls. The room is the superellipsoid |x|^p + |y|^p + |z|^p = 1, and with time dependence e^{-iωt} the pressure amplitude satisfies Δu + k²u = −δ(x − x₀) inside and ∂u/∂n − ikβu = 0 on the walls, where β > 0 is the wall admittance. The page lets the user set k, β, the discretization, the room shape and the source position, and shows the field on three planes through the source and on the walls, in 3D, animated over one period.
8+
9+## Method
10+
11+The equation, representation and discretization are those of the CUDA version: the correction v = u − G is represented as S_k[σ] + iD_k S_{i|k|}[σ], on 6m² curved quadrilateral patches with (order + 1)² Gauss–Legendre nodes each. The near part of each layer potential is computed once with adaptive quadrature (polar coordinates with a sinh substitution on the target's own patch, a quadtree on neighbouring patches) and stored as weights on the GPU. The far part is a direct O(N²) sum over oversampled nodes, recomputed at every GMRES iteration. GMRES runs on the CPU in double precision; only the operator is applied on the GPU.
12+
13+Unlike the CUDA version, which uses double precision for the near quadrature, everything on the GPU is single precision, since WebGPU has no double precision type. Three changes make this sufficient:
14+
15+1. For a target on the wall, the vector from the target to a nearby quadrature node is formed in the coordinates of the superellipsoid's analytic parametrization rather than as a difference of two rounded positions. It then keeps full relative precision however short it is, which the double-layer kernels need, since their numerators (y − x)·n are O(r²).
16+2. The weakly singular difference D′_k − D′_{i|k|} is evaluated from its power series when kr < 1, where the direct formula subtracts two nearly equal O(1/r³) terms.
17+3. sin and cos are evaluated with our own range reduction and polynomials. WGSL only requires them to be accurate to 2⁻¹¹ in absolute terms, and some implementations (SwiftShader among them) are no better than that, which limited the near weights to about 1e-4.
18+
19+## Accuracy
20+
21+`npm run test:gpu` compares the near weights with a brute-force double precision reference, checks Green's identity through the full operator, and solves a problem with a known solution (the field of a source outside the room). On 54 patches of order 6 at k = 4 with the SwiftShader software adapter, the near weights agree with the reference to 7e-6 or better, the Green's identity residual is 3.4e-6 (3.6e-6 in double precision on the CPU for the same discretization), and the self-test error is 1.6e-6 in the interior and 9.0e-6 on the walls. For comparison, the CUDA version reports 2.2e-6 on its small case.
22+
23+For sharper rooms the discretization, not the precision, limits the accuracy. At p = 30 the same surface gives a Green's identity residual of 1.1e-4 on the GPU and 1.1e-4 in double precision; refining to 96 patches of order 8 brings the self-test error down to 1.7e-5.
24+
25+## Running it
26+
27+ npm install
28+ npm run dev # the app at http://localhost:5173
29+ npm run test:node # CPU checks of the surface and the reference quadrature
30+ npm run test:gpu # the GPU validation suite in headless Chrome
31+ npm run smoke # build, then solve the quick case in the app, headless
32+
33+The headless tests use a hardware WebGPU adapter when Chrome finds one and SwiftShader otherwise. Note that in headless Chrome, presenting a frame to a WebGPU canvas makes later asynchronous WebGPU calls fail ("a valid external Instance reference no longer exists"), so the smoke test opens the page with `?noroom`, which skips the 3D view. We have not seen this in a browser with a window.
34+
35+## Limitations
36+
37+The near weights take 6 × (order + 1)² × 8 bytes per pair of a node and a near patch. At order 8 this is several hundred MB for the medium case (17,496 points) and more than most browsers allow in one buffer for the large case (48,600 points). Splitting the weights across buffers, or recomputing them, would remove this limit.
38+
39+The far field is a direct sum, so each operator application costs O(N²); a GPU FMM would replace it for larger problems without changing the near weights.
40+
41+Only the rounded-cube family of rooms is available, because the single precision near quadrature relies on its exact parametrization. A general surface would need the interpolated geometry and a different way of forming short displacements accurately.
42+
43+Timings on the software adapter are not representative, and we have not yet measured the named cases on hardware GPUs.
index.htmladded+180−0View file
@@ -0,0 +1,180 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="utf-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1" />
6+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><rect width=%22100%22 height=%22100%22 rx=%2224%22 fill=%22%230b0e12%22/><rect x=%2218%22 y=%2218%22 width=%2264%22 height=%2264%22 rx=%2218%22 fill=%22none%22 stroke=%22%239aa4af%22 stroke-width=%225%22/><circle cx=%2242%22 cy=%2256%22 r=%228%22 fill=%22%23b2182b%22/><circle cx=%2242%22 cy=%2256%22 r=%2219%22 fill=%22none%22 stroke=%22%232166ac%22 stroke-width=%225%22/></svg>" />
7+ <title>fmm3dbie-webgpu-poc: a room with impedance walls, solved on your GPU</title>
8+ <style>
9+ :root {
10+ --bg: #ffffff;
11+ --ink: #1f2328;
12+ --ink-2: #57606a;
13+ --line: #d0d7de;
14+ --accent: #0969da;
15+ --panel-bg: #f6f8fa;
16+ --warn: #9a5b00;
17+ --bad: #cf222e;
18+ color-scheme: light dark;
19+ }
20+ @media (prefers-color-scheme: dark) {
21+ :root {
22+ --bg: #14171a;
23+ --ink: #e6e9ec;
24+ --ink-2: #9aa4af;
25+ --line: #333b44;
26+ --accent: #58a6ff;
27+ --panel-bg: #1a1e23;
28+ --warn: #e3a341;
29+ --bad: #ff7b72;
30+ }
31+ }
32+ * { box-sizing: border-box; }
33+ body { margin: 0; background: var(--bg); color: var(--ink); font: 15px/1.5 system-ui, -apple-system, sans-serif; }
34+ main { max-width: 1180px; margin: 0 auto; padding: 20px 16px 56px; }
35+ h1 { font-size: 20px; margin: 0 0 2px; }
36+ h2 { font-size: 15px; margin: 22px 0 6px; }
37+ .sub { color: var(--ink-2); margin: 0 0 12px; font-size: 13.5px; max-width: 78ch; }
38+ a { color: var(--accent); }
39+ .panel { background: var(--panel-bg); border: 1px solid var(--line); border-radius: 8px; padding: 10px 14px; }
40+ .grid { display: grid; grid-template-columns: repeat(auto-fill, minmax(150px, 1fr)); gap: 6px 16px; }
41+ .field { display: flex; flex-direction: column; gap: 2px; font-size: 12.5px; color: var(--ink-2); }
42+ .field > span { white-space: nowrap; }
43+ input[type='number'], select, button {
44+ font: inherit; font-size: 13.5px; color: var(--ink); background: var(--bg);
45+ border: 1px solid var(--line); border-radius: 6px; padding: 4px 8px; min-width: 0;
46+ }
47+ input[type='number'] { width: 100%; font-variant-numeric: tabular-nums; }
48+ .xyz { display: grid; grid-template-columns: 1fr 1fr 1fr; gap: 4px; }
49+ button { cursor: pointer; }
50+ button:hover:not(:disabled) { border-color: var(--accent); }
51+ button:disabled { opacity: 0.5; cursor: default; }
52+ button.primary { border-color: var(--accent); color: var(--accent); font-weight: 600; min-width: 6em; }
53+ .row { display: flex; flex-wrap: wrap; align-items: center; gap: 8px 14px; margin-top: 10px; }
54+ .row label { font-size: 13px; color: var(--ink-2); white-space: nowrap; }
55+ .derived { font-size: 13px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
56+ .derived b { color: var(--ink); font-weight: 600; }
57+ .warn { color: var(--warn); }
58+ .bad { color: var(--bad); }
59+ progress { width: 220px; height: 10px; vertical-align: middle; accent-color: var(--accent); }
60+ #status { font-size: 13px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
61+ #err { color: var(--bad); white-space: pre-wrap; font-size: 13px; font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace; margin: 8px 0 0; }
62+ #err:empty { display: none; }
63+ .view-bar { display: flex; flex-wrap: wrap; gap: 6px 16px; align-items: center; margin: 14px 0 6px; font-size: 13px; color: var(--ink-2); }
64+ .stage { display: grid; grid-template-columns: minmax(0, 1fr) 360px; gap: 14px; align-items: start; }
65+ @media (max-width: 900px) { .stage { grid-template-columns: minmax(0, 1fr); } }
66+ .room { position: relative; border: 1px solid var(--line); border-radius: 8px; overflow: hidden; background: #0b0e12; aspect-ratio: 4 / 3; }
67+ #room { width: 100%; height: 100%; display: block; cursor: grab; touch-action: none; }
68+ #room:active { cursor: grabbing; }
69+ .room .hint { position: absolute; left: 10px; bottom: 8px; font-size: 11.5px; color: #9aa4af; pointer-events: none; }
70+ .room .empty { position: absolute; inset: 0; display: flex; align-items: center; justify-content: center; color: #9aa4af; font-size: 14px; text-align: center; padding: 20px; }
71+ .cbar { position: absolute; right: 10px; top: 10px; bottom: 10px; display: flex; flex-direction: column; align-items: center; gap: 4px; color: #c9d1d9; font-size: 11px; font-variant-numeric: tabular-nums; pointer-events: none; }
72+ .cbar canvas { flex: 1; width: 14px; min-height: 60px; border: 1px solid #333b44; border-radius: 2px; }
73+ .side { display: flex; flex-direction: column; gap: 12px; }
74+ table { border-collapse: collapse; width: 100%; font-size: 13px; font-variant-numeric: tabular-nums; }
75+ td { padding: 3px 4px; border-bottom: 1px solid var(--line); vertical-align: top; }
76+ td:first-child { color: var(--ink-2); }
77+ td:last-child { text-align: right; }
78+ #conv { width: 100%; height: 150px; display: block; }
79+ .slices { display: grid; grid-template-columns: repeat(3, minmax(0, 1fr)); gap: 12px; margin-top: 14px; }
80+ @media (max-width: 700px) { .slices { grid-template-columns: minmax(0, 1fr); } }
81+ figure.slice { margin: 0; }
82+ figure.slice canvas { width: 100%; aspect-ratio: 1; display: block; border: 1px solid var(--line); border-radius: 6px; background: repeating-conic-gradient(var(--panel-bg) 0 25%, var(--bg) 0 50%) 0 0 / 16px 16px; cursor: crosshair; }
83+ figure.slice figcaption { font-size: 12px; color: var(--ink-2); margin-top: 3px; }
84+ #hover { min-height: 1.5em; font-size: 12.5px; color: var(--ink-2); font-variant-numeric: tabular-nums; margin-top: 4px; }
85+ .notes { font-size: 13.5px; color: var(--ink-2); max-width: 82ch; }
86+ .notes p { margin: 6px 0; }
87+ .notes code { font-size: 12.5px; }
88+ #device { font-size: 12px; color: var(--ink-2); margin-top: 18px; }
89+ </style>
90+ </head>
91+ <body>
92+ <main>
93+ <h1>fmm3dbie-webgpu-poc</h1>
94+ <p class="sub">
95+ A point source at a fixed frequency in a closed room with partially absorbing walls, solved with a high-order boundary integral equation on your GPU through WebGPU, in single precision. This is a browser port of a <a href="https://vault1.magland.org/proof-of-concept/fmm3dbie-gpu-poc">CUDA proof of concept</a> for
96+ <a href="https://github.com/fastalgorithms/fmm3dbie">fmm3dbie</a>'s interior Helmholtz impedance solver. Choose the problem, press Solve, then turn the room with the mouse and click a slice to move the source.
97+ </p>
98+
99+ <div class="panel">
100+ <div class="grid">
101+ <label class="field" title="Named cases from the CUDA proof of concept, plus a quick one for slow GPUs"><span>case</span>
102+ <select autocomplete="off" id="preset">
103+ <option value="quick">quick (1,176 points)</option>
104+ <option value="small" selected>small (4,704 points)</option>
105+ <option value="medium">medium (17,496 points)</option>
106+ <option value="large">large (48,600 points)</option>
107+ <option value="custom">custom</option>
108+ </select>
109+ </label>
110+ <label class="field" title="Wavenumber k = 2π / wavelength, for a room of half-width 1"><span>wavenumber k</span><input autocomplete="off" id="k" type="number" min="0.5" max="60" step="0.5" /></label>
111+ <label class="field" title="Normalized wall admittance β > 0; smaller β means more reflective walls"><span>wall admittance β</span><input autocomplete="off" id="beta" type="number" min="0.001" max="5" step="0.01" /></label>
112+ <label class="field" title="Patches per edge of each cube face (6 m² patches in all)"><span>patches per face edge m</span><input autocomplete="off" id="m" type="number" min="1" max="16" step="1" /></label>
113+ <label class="field" title="Polynomial order per patch; (order + 1)² nodes each"><span>order per patch</span>
114+ <select autocomplete="off" id="norder"><option>4</option><option>5</option><option>6</option><option>7</option><option>8</option><option>9</option><option>10</option></select>
115+ </label>
116+ <label class="field" title="The walls are |x|^p + |y|^p + |z|^p = 1: p = 2 is a sphere, larger p a sharper cube, which needs more patches or a higher order for the same accuracy"><span>room shape exponent p</span>
117+ <select autocomplete="off" id="p"><option>2</option><option>4</option><option>6</option><option selected>8</option><option>10</option><option>12</option><option>16</option><option>20</option><option>30</option></select>
118+ </label>
119+ <label class="field" title="Source position; must be inside the room. You can also click a slice."><span>source (x, y, z)</span>
120+ <span class="xyz"><input autocomplete="off" id="x0" type="number" step="0.05" /><input autocomplete="off" id="y0" type="number" step="0.05" /><input autocomplete="off" id="z0" type="number" step="0.05" /></span>
121+ </label>
122+ <label class="field" title="Relative residual at which GMRES stops"><span>GMRES tolerance</span>
123+ <select autocomplete="off" id="eps"><option value="1e-4">1e-4</option><option value="1e-5" selected>1e-5</option><option value="1e-6">1e-6</option></select>
124+ </label>
125+ <label class="field" title="Grid points per side of the three slices"><span>slice resolution</span>
126+ <select autocomplete="off" id="ngrid"><option>61</option><option selected>101</option><option>151</option></select>
127+ </label>
128+ </div>
129+ <div class="row">
130+ <button id="solve" class="primary">Solve</button>
131+ <button id="cancel" disabled>Cancel</button>
132+ <label title="Also solve a problem with a known solution (the field of a source outside the room) and report the error"><input autocomplete="off" id="selftest" type="checkbox" /> self-test</label>
133+ <span class="derived" id="derived"></span>
134+ </div>
135+ <div class="row" style="margin-top: 6px"><progress id="progress" max="1" value="0"></progress><span id="status">loading…</span></div>
136+ <pre id="err"></pre>
137+ </div>
138+
139+ <div class="view-bar">
140+ <label>show
141+ <select autocomplete="off" id="mode">
142+ <option value="pressure">pressure Re(u e^−iωt)</option>
143+ <option value="level">level (dB re median)</option>
144+ </select>
145+ </label>
146+ <label><input autocomplete="off" id="animate" type="checkbox" checked /> animate over one period</label>
147+ <label><input autocomplete="off" id="showWalls" type="checkbox" checked /> walls</label>
148+ <label><input autocomplete="off" id="showXY" type="checkbox" checked /> xy</label>
149+ <label><input autocomplete="off" id="showXZ" type="checkbox" checked /> xz</label>
150+ <label><input autocomplete="off" id="showYZ" type="checkbox" /> yz</label>
151+ </div>
152+
153+ <div class="stage">
154+ <div class="room">
155+ <canvas id="room"></canvas>
156+ <div class="empty" id="empty">Press Solve to compute the field.</div>
157+ <div class="cbar"><span id="cbarHi"></span><canvas id="cbar" width="14" height="256"></canvas><span id="cbarLo"></span></div>
158+ <div class="hint">drag to turn · scroll to zoom · the nearer walls are cut away</div>
159+ </div>
160+ <div class="side">
161+ <div class="panel">
162+ <table id="results"><tbody><tr><td>no solve yet</td><td></td></tr></tbody></table>
163+ </div>
164+ <div class="panel"><canvas id="conv"></canvas></div>
165+ </div>
166+ </div>
167+
168+ <div class="slices" id="slices"></div>
169+ <div id="hover"></div>
170+
171+ <h2>What is computed</h2>
172+ <div class="notes">
173+ <p>The room is the superellipsoid |x|<sup>p</sup> + |y|<sup>p</sup> + |z|<sup>p</sup> = 1, with a point source inside and impedance walls (∂u/∂n − ikβu = 0). We solve fmm3dbie's combined field integral equation on 6m² curved patches: near interactions use precomputed quadrature weights, far interactions a direct sum, and GMRES runs on the CPU. The GPU work is in single precision.</p>
174+ <p>The self-test solves a problem with a known solution and reports the error. The energy mismatch compares the power from the source with the power absorbed by the walls. Large cases may need more GPU memory than the browser allows.</p>
175+ </div>
176+ <div id="device"></div>
177+ </main>
178+ <script type="module" src="/src/main.ts"></script>
179+ </body>
180+</html>
package-lock.jsonadded+2142−0View file
This diff is 2,147 lines long and is not shown.
package.jsonadded+25−0View file
@@ -0,0 +1,25 @@
1+{
2+ "name": "fmm3dbie-webgpu-poc",
3+ "version": "0.1.0",
4+ "description": "A boundary integral equation solver for the Helmholtz equation in a room with impedance walls, running on WebGPU in the browser",
5+ "type": "module",
6+ "engines": {
7+ "node": ">=22.6"
8+ },
9+ "license": "Apache-2.0",
10+ "scripts": {
11+ "dev": "vite",
12+ "build": "tsc --noEmit && vite build",
13+ "test:node": "node --experimental-strip-types test/test-node.ts",
14+ "test:gpu": "vite build && node scripts/test-gpu.mjs",
15+ "smoke": "vite build && node scripts/smoke.mjs",
16+ "test": "npm run test:node && npm run test:gpu"
17+ },
18+ "devDependencies": {
19+ "@types/node": "^26.1.1",
20+ "@webgpu/types": "^0.1.44",
21+ "puppeteer-core": "^23.11.1",
22+ "typescript": "^5.5.0",
23+ "vite": "^5.4.0"
24+ }
25+}
scripts/smoke.mjsadded+94−0View file
@@ -0,0 +1,94 @@
1+/**
2+ * Does the app work in a browser? Serves dist/, opens index.html in headless
3+ * Chrome (hardware WebGPU if there is any, SwiftShader otherwise), solves
4+ * the quick case, and waits for the status line to report
5+ * success. Any console error, page error or failed request fails the run. A
6+ * screenshot goes to tmp/smoke.png.
7+ *
8+ * Run after `vite build`: node scripts/smoke.mjs [case] (default quick)
9+ */
10+import { createServer } from 'node:http';
11+import { mkdir, readFile } from 'node:fs/promises';
12+import { extname, join } from 'node:path';
13+import puppeteer from 'puppeteer-core';
14+
15+const DIST = new URL('../dist/', import.meta.url).pathname;
16+const TMP = new URL('../tmp/', import.meta.url).pathname;
17+const CHROME = process.env.CHROME_PATH ?? '/usr/bin/google-chrome';
18+const CASE = process.argv[2] ?? 'quick';
19+const MIME = { '.html': 'text/html', '.js': 'text/javascript', '.css': 'text/css', '.json': 'application/json', '.svg': 'image/svg+xml' };
20+
21+const server = createServer(async (req, res) => {
22+ try {
23+ const path = req.url === '/' ? '/index.html' : req.url.split('?')[0];
24+ const data = await readFile(join(DIST, path));
25+ res.writeHead(200, { 'content-type': MIME[extname(path)] ?? 'application/octet-stream' });
26+ res.end(data);
27+ } catch {
28+ res.writeHead(404);
29+ res.end('not found');
30+ }
31+});
32+await new Promise((r) => server.listen(0, '127.0.0.1', r));
33+const port = server.address().port;
34+await mkdir(TMP, { recursive: true });
35+
36+const flagSets = [
37+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--enable-features=Vulkan'],
38+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--use-webgpu-adapter=swiftshader', '--enable-unsafe-swiftshader'],
39+];
40+if (process.env.SWIFTSHADER_ONLY) flagSets.shift();
41+
42+let ok = false;
43+for (const flags of flagSets) {
44+ const browser = await puppeteer.launch({ executablePath: CHROME, args: [...flags, '--window-size=1280,1400'], protocolTimeout: 1_800_000 });
45+ const problems = [];
46+ try {
47+ const page = await browser.newPage();
48+ await page.setViewport({ width: 1280, height: 1400 });
49+ page.on('console', (m) => {
50+ if (process.env.VERBOSE) console.log(` [${m.type()}] ${m.text()}`);
51+ if (m.type() === 'error') problems.push(`console: ${m.text()}`);
52+ });
53+ page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`));
54+ page.on('requestfailed', (r) => problems.push(`request failed: ${r.url()}`));
55+ await page.goto(`http://127.0.0.1:${port}/index.html?noroom`, { waitUntil: 'load' });
56+ await page.waitForFunction(() => /^ready/.test(document.getElementById('status')?.textContent ?? ''), { timeout: 60_000, polling: 250 });
57+ console.log(await page.$eval('#device', (e) => e.textContent));
58+ await page.select('#preset', CASE);
59+ await page.select('#ngrid', '61');
60+ const t0 = Date.now();
61+ await page.click('#solve');
62+ let lastStatus = '';
63+ const poll = setInterval(async () => {
64+ try {
65+ const s = await page.$eval('#status', (e) => e.textContent);
66+ if (s !== lastStatus && !/iteration/.test(s)) console.log(` ${s}`);
67+ lastStatus = s;
68+ } catch {}
69+ }, 2000);
70+ await page.waitForFunction(() => /^(done|failed)/.test(document.getElementById('status')?.textContent ?? ''), { timeout: 1_800_000, polling: 500 });
71+ clearInterval(poll);
72+ const status = await page.$eval('#status', (e) => e.textContent);
73+ const err = await page.$eval('#err', (e) => e.textContent);
74+ const rows = await page.$$eval('#results tr', (trs) => trs.map((tr) => [...tr.children].map((td) => td.textContent).join(': ')));
75+ console.log(`status: ${status} (${((Date.now() - t0) / 1000).toFixed(1)} s)`);
76+ for (const r of rows) console.log(` ${r}`);
77+ if (err) problems.push(`error line: ${err}`);
78+ await new Promise((r) => setTimeout(r, 1000));
79+ await page.screenshot({ path: join(TMP, 'smoke.png'), fullPage: true });
80+ console.log(`screenshot: ${join(TMP, 'smoke.png')}`);
81+ if (!/^done/.test(status)) problems.push(`status: ${status}`);
82+ ok = problems.length === 0;
83+ } catch (e) {
84+ problems.push(`harness: ${e.message}`);
85+ } finally {
86+ await browser.close();
87+ }
88+ for (const p of problems) console.log(`problem: ${p}`);
89+ if (ok) break;
90+ console.log('retrying with next flag set…');
91+}
92+server.close();
93+console.log(ok ? 'smoke test passed' : 'smoke test FAILED');
94+process.exit(ok ? 0 : 1);
scripts/test-gpu.mjsadded+72−0View file
@@ -0,0 +1,72 @@
1+/**
2+ * Headless GPU test runner: serves dist/, opens test.html in headless Chrome
3+ * (hardware WebGPU if there is any, the SwiftShader software adapter
4+ * otherwise), and reports the suite's results.
5+ *
6+ * Run after `vite build`: node scripts/test-gpu.mjs [?m=2&norder=4&k=4]
7+ */
8+import { createServer } from 'node:http';
9+import { readFile } from 'node:fs/promises';
10+import { extname, join } from 'node:path';
11+import puppeteer from 'puppeteer-core';
12+
13+const DIST = new URL('../dist/', import.meta.url).pathname;
14+const CHROME = process.env.CHROME_PATH ?? '/usr/bin/google-chrome';
15+const QUERY = process.argv[2] ?? '';
16+const MIME = {
17+ '.html': 'text/html',
18+ '.js': 'text/javascript',
19+ '.css': 'text/css',
20+ '.json': 'application/json',
21+};
22+
23+const server = createServer(async (req, res) => {
24+ try {
25+ const path = req.url === '/' ? '/index.html' : req.url.split('?')[0];
26+ const data = await readFile(join(DIST, path));
27+ res.writeHead(200, { 'content-type': MIME[extname(path)] ?? 'application/octet-stream' });
28+ res.end(data);
29+ } catch {
30+ res.writeHead(404);
31+ res.end('not found');
32+ }
33+});
34+await new Promise((r) => server.listen(0, '127.0.0.1', r));
35+const port = server.address().port;
36+
37+const flagSets = [
38+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--enable-features=Vulkan'],
39+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--use-webgpu-adapter=swiftshader', '--enable-unsafe-swiftshader'],
40+];
41+if (process.env.SWIFTSHADER_ONLY) flagSets.shift();
42+
43+let final = null;
44+for (const flags of flagSets) {
45+ const browser = await puppeteer.launch({ executablePath: CHROME, args: flags, protocolTimeout: 1_800_000 });
46+ try {
47+ const page = await browser.newPage();
48+ page.on('console', (msg) => {
49+ const t = msg.text();
50+ if (!/^(ok|FAIL|surface|device|operator|FATAL|\s)/.test(t)) console.log(` [page] ${t}`);
51+ else console.log(t);
52+ });
53+ page.on('pageerror', (err) => console.log(` [pageerror] ${err.message}`));
54+ await page.goto(`http://127.0.0.1:${port}/test.html${QUERY}`, { waitUntil: 'load' });
55+ const results = await page.waitForFunction(() => window.__RESULTS__, { timeout: 1_800_000, polling: 500 });
56+ final = await results.jsonValue();
57+ } catch (e) {
58+ console.error(`run with flags [${flags.join(' ')}] failed: ${e.message}`);
59+ } finally {
60+ await browser.close();
61+ }
62+ if (final && !final.fatal) break;
63+ if (final?.fatal) console.log(`fatal: ${final.fatal}`);
64+ console.log('retrying with next flag set…');
65+}
66+server.close();
67+if (!final) {
68+ console.error('no results');
69+ process.exit(1);
70+}
71+console.log(`\n${final.passed} passed, ${final.failed} failed`);
72+process.exit(final.failed ? 1 : 0);
src/gauss.tsadded+41−0View file
@@ -0,0 +1,41 @@
1+/** Gauss-Legendre nodes (ascending) and weights on [-1, 1], by Newton iteration. */
2+export function leggauss(n: number): { x: Float64Array; w: Float64Array } {
3+ const x = new Float64Array(n);
4+ const w = new Float64Array(n);
5+ const m = (n + 1) >> 1;
6+ for (let i = 0; i < m; i++) {
7+ let z = Math.cos((Math.PI * (i + 0.75)) / (n + 0.5));
8+ let pp = 0;
9+ for (let iter = 0; iter < 100; iter++) {
10+ let p1 = 1;
11+ let p2 = 0;
12+ for (let j = 1; j <= n; j++) {
13+ const p3 = p2;
14+ p2 = p1;
15+ p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
16+ }
17+ pp = (n * (z * p1 - p2)) / (z * z - 1);
18+ const z1 = z;
19+ z = z1 - p1 / pp;
20+ if (Math.abs(z - z1) < 1e-15) break;
21+ }
22+ x[i] = -z;
23+ x[n - 1 - i] = z;
24+ w[i] = w[n - 1 - i] = 2 / ((1 - z * z) * pp * pp);
25+ }
26+ return { x, w };
27+}
28+
29+/** P_0(x), ..., P_norder(x) for each x: array of shape (norder+1, n), index k*n + i. */
30+export function legendreVander(x: ArrayLike<number>, norder: number): Float64Array {
31+ const n = x.length;
32+ const P = new Float64Array((norder + 1) * n);
33+ for (let i = 0; i < n; i++) {
34+ P[i] = 1;
35+ if (norder >= 1) P[n + i] = x[i];
36+ for (let k = 1; k < norder; k++) {
37+ P[(k + 1) * n + i] = ((2 * k + 1) * x[i] * P[k * n + i] - k * P[(k - 1) * n + i]) / (k + 1);
38+ }
39+ }
40+ return P;
41+}
src/gmres.tsadded+178−0View file
@@ -0,0 +1,178 @@
1+/**
2+ * Unrestarted GMRES for (zid I + K) x = rhs with the weighted inner product
3+ * <u, v> = sum conj(u) v w, following fmm3dbie's zgmres_guru. The Krylov
4+ * vectors and the Hessenberg system are kept in double precision on the CPU;
5+ * only the operator K is applied on the GPU.
6+ */
7+import type { CVec } from './problem.ts';
8+
9+export type Apply = (x: CVec) => Promise<CVec>;
10+
11+export interface GmresResult {
12+ x: CVec;
13+ niter: number;
14+ rres: number;
15+ residuals: number[];
16+}
17+
18+export async function gmres(
19+ applyK: Apply,
20+ rhs: CVec,
21+ wts: Float64Array,
22+ zid: [number, number],
23+ eps: number,
24+ numit: number,
25+ onIter?: (it: number, rres: number) => void,
26+): Promise<GmresResult> {
27+ const n = wts.length;
28+ const wnorm = (v: CVec) => {
29+ let s = 0;
30+ for (let i = 0; i < n; i++) s += (v[2 * i] ** 2 + v[2 * i + 1] ** 2) * wts[i];
31+ return Math.sqrt(s);
32+ };
33+ const rb = wnorm(rhs);
34+ const V: CVec[] = [];
35+ const m1 = numit + 1;
36+ const Hr = new Float64Array(m1 * numit);
37+ const Hi = new Float64Array(m1 * numit);
38+ const csr = new Float64Array(numit);
39+ const csi = new Float64Array(numit);
40+ const snr = new Float64Array(numit);
41+ const sni = new Float64Array(numit);
42+ const sr = new Float64Array(m1);
43+ const si = new Float64Array(m1);
44+ sr[0] = rb;
45+ const v0 = new Float64Array(2 * n);
46+ for (let i = 0; i < 2 * n; i++) v0[i] = rhs[i] / rb;
47+ V.push(v0);
48+ const residuals: number[] = [];
49+ let itDone = 0;
50+ for (let it = 0; it < numit; it++) {
51+ let w = await applyK(V[it]);
52+ w = new Float64Array(w); // do not alias the operator's output
53+ // classical Gram-Schmidt against all previous vectors, repeated once
54+ for (let pass = 0; pass < 2; pass++) {
55+ for (let k = 0; k <= it; k++) {
56+ const vk = V[k];
57+ let hr = 0;
58+ let hi = 0;
59+ for (let i = 0; i < n; i++) {
60+ const ar = vk[2 * i];
61+ const ai = vk[2 * i + 1];
62+ const br = w[2 * i] * wts[i];
63+ const bi = w[2 * i + 1] * wts[i];
64+ hr += ar * br + ai * bi; // conj(a) b
65+ hi += ar * bi - ai * br;
66+ }
67+ Hr[k * numit + it] += hr;
68+ Hi[k * numit + it] += hi;
69+ for (let i = 0; i < n; i++) {
70+ const vr = vk[2 * i];
71+ const vi = vk[2 * i + 1];
72+ w[2 * i] -= hr * vr - hi * vi;
73+ w[2 * i + 1] -= hr * vi + hi * vr;
74+ }
75+ }
76+ }
77+ Hr[it * numit + it] += zid[0];
78+ Hi[it * numit + it] += zid[1];
79+ const hn = wnorm(w);
80+ Hr[(it + 1) * numit + it] = hn;
81+ Hi[(it + 1) * numit + it] = 0;
82+ const vn = new Float64Array(2 * n);
83+ if (hn > 0) for (let i = 0; i < 2 * n; i++) vn[i] = w[i] / hn;
84+ V.push(vn);
85+ // apply the previous Givens rotations to the new column
86+ for (let k = 0; k < it; k++) {
87+ const ar = Hr[k * numit + it];
88+ const ai = Hi[k * numit + it];
89+ const br = Hr[(k + 1) * numit + it];
90+ const bi = Hi[(k + 1) * numit + it];
91+ // t1 = cs a + sn b ; t2 = -conj(sn) a + cs b
92+ const t1r = csr[k] * ar - csi[k] * ai + snr[k] * br - sni[k] * bi;
93+ const t1i = csr[k] * ai + csi[k] * ar + snr[k] * bi + sni[k] * br;
94+ const t2r = -(snr[k] * ar + sni[k] * ai) + csr[k] * br - csi[k] * bi;
95+ const t2i = -(snr[k] * ai - sni[k] * ar) + csr[k] * bi + csi[k] * br;
96+ Hr[k * numit + it] = t1r;
97+ Hi[k * numit + it] = t1i;
98+ Hr[(k + 1) * numit + it] = t2r;
99+ Hi[(k + 1) * numit + it] = t2i;
100+ }
101+ const ar = Hr[it * numit + it];
102+ const ai = Hi[it * numit + it];
103+ const br = Hr[(it + 1) * numit + it];
104+ const bi = Hi[(it + 1) * numit + it];
105+ const aa = Math.hypot(ar, ai);
106+ const r = Math.hypot(aa, Math.hypot(br, bi));
107+ if (r === 0) {
108+ csr[it] = 1;
109+ csi[it] = 0;
110+ snr[it] = 0;
111+ sni[it] = 0;
112+ } else {
113+ csr[it] = aa / r;
114+ csi[it] = 0;
115+ // sn = (a/|a|) conj(b) / r
116+ let pr = 1;
117+ let pi = 0;
118+ if (aa > 0) {
119+ pr = ar / aa;
120+ pi = ai / aa;
121+ }
122+ snr[it] = (pr * br + pi * bi) / r;
123+ sni[it] = (pi * br - pr * bi) / r;
124+ }
125+ // H[it, it] = cs a + sn b ; H[it+1, it] = 0
126+ Hr[it * numit + it] = csr[it] * ar - csi[it] * ai + snr[it] * br - sni[it] * bi;
127+ Hi[it * numit + it] = csr[it] * ai + csi[it] * ar + snr[it] * bi + sni[it] * br;
128+ Hr[(it + 1) * numit + it] = 0;
129+ Hi[(it + 1) * numit + it] = 0;
130+ // s[it+1] = -conj(sn) s[it] ; s[it] = cs s[it]
131+ const s0r = sr[it];
132+ const s0i = si[it];
133+ sr[it + 1] = -(snr[it] * s0r + sni[it] * s0i);
134+ si[it + 1] = -(snr[it] * s0i - sni[it] * s0r);
135+ sr[it] = csr[it] * s0r - csi[it] * s0i;
136+ si[it] = csr[it] * s0i + csi[it] * s0r;
137+ itDone = it + 1;
138+ const rres = Math.hypot(sr[it + 1], si[it + 1]) / rb;
139+ residuals.push(rres);
140+ onIter?.(itDone, rres);
141+ if (rres < eps) break;
142+ }
143+ // back substitution on the upper triangular system
144+ const yr = new Float64Array(itDone);
145+ const yi = new Float64Array(itDone);
146+ for (let i = itDone - 1; i >= 0; i--) {
147+ let tr = sr[i];
148+ let ti = si[i];
149+ for (let j = i + 1; j < itDone; j++) {
150+ const hr = Hr[i * numit + j];
151+ const hi = Hi[i * numit + j];
152+ tr -= hr * yr[j] - hi * yi[j];
153+ ti -= hr * yi[j] + hi * yr[j];
154+ }
155+ const hr = Hr[i * numit + i];
156+ const hi = Hi[i * numit + i];
157+ const d = hr * hr + hi * hi;
158+ yr[i] = (tr * hr + ti * hi) / d;
159+ yi[i] = (ti * hr - tr * hi) / d;
160+ }
161+ const x = new Float64Array(2 * n);
162+ for (let j = 0; j < itDone; j++) {
163+ const vj = V[j];
164+ for (let i = 0; i < n; i++) {
165+ const vr = vj[2 * i];
166+ const vi = vj[2 * i + 1];
167+ x[2 * i] += yr[j] * vr - yi[j] * vi;
168+ x[2 * i + 1] += yr[j] * vi + yi[j] * vr;
169+ }
170+ }
171+ const Kx = await applyK(x);
172+ const res = new Float64Array(2 * n);
173+ for (let i = 0; i < n; i++) {
174+ res[2 * i] = zid[0] * x[2 * i] - zid[1] * x[2 * i + 1] + Kx[2 * i] - rhs[2 * i];
175+ res[2 * i + 1] = zid[0] * x[2 * i + 1] + zid[1] * x[2 * i] + Kx[2 * i + 1] - rhs[2 * i + 1];
176+ }
177+ return { x, niter: itDone, rres: wnorm(res) / rb, residuals };
178+}
src/gpu/device.tsadded+180−0View file
@@ -0,0 +1,180 @@
1+/** The GPU device, plus small helpers for buffers, uniforms, readback and paced dispatching. */
2+
3+export async function requestSolverDevice(): Promise<Gpu> {
4+ if (!navigator.gpu) {
5+ throw new Error(
6+ 'this browser has no WebGPU. Chrome and Edge 113+, Safari 26+, and Firefox 141+ on Windows have it; ' +
7+ 'on Linux, Chrome may need chrome://flags/#enable-unsafe-webgpu and Firefox dom.webgpu.enabled.',
8+ );
9+ }
10+ const adapter = await navigator.gpu.requestAdapter({ powerPreference: 'high-performance' });
11+ if (!adapter) throw new Error('WebGPU found no adapter on this machine.');
12+ const lim = adapter.limits;
13+ // The near weights of a medium-sized problem are a few hundred MB, so ask
14+ // for the largest buffers the adapter will give rather than the 128 MB default.
15+ const device = await adapter.requestDevice({
16+ requiredLimits: {
17+ maxStorageBufferBindingSize: lim.maxStorageBufferBindingSize,
18+ maxBufferSize: lim.maxBufferSize,
19+ maxComputeWorkgroupStorageSize: Math.max(16384, Math.min(lim.maxComputeWorkgroupStorageSize, 32768)),
20+ },
21+ });
22+ device.lost.then((info) => {
23+ if (info.reason !== 'destroyed') console.error(`WebGPU device lost: ${info.message}`);
24+ });
25+ return new Gpu(device, adapter);
26+}
27+
28+export function adapterDescription(adapter: GPUAdapter | null): string {
29+ if (!adapter) return 'no adapter';
30+ const info = (adapter as GPUAdapter & { info?: GPUAdapterInfo }).info;
31+ if (!info) return 'unknown adapter';
32+ return [info.description, info.device, info.architecture, info.vendor].filter((s) => s).join(', ') || 'unknown adapter';
33+}
34+
35+/** A device wrapper with the conveniences the solver needs. */
36+export class Gpu {
37+ private pipelines = new Map<string, GPUComputePipeline>();
38+
39+ /**
40+ * The adapter is held for as long as the device is used: in Chrome, async
41+ * device calls fail ("a valid external Instance reference no longer exists")
42+ * once the adapter has been garbage collected.
43+ */
44+ constructor(
45+ readonly device: GPUDevice,
46+ readonly adapter?: GPUAdapter,
47+ ) {}
48+
49+ storage(data: ArrayBufferView | number, extraUsage: GPUBufferUsageFlags = 0): GPUBuffer {
50+ const size = typeof data === 'number' ? data : data.byteLength;
51+ const buf = this.device.createBuffer({
52+ size: Math.max(16, Math.ceil(size / 16) * 16),
53+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC | extraUsage,
54+ });
55+ if (typeof data !== 'number') this.device.queue.writeBuffer(buf, 0, data.buffer, data.byteOffset, data.byteLength);
56+ return buf;
57+ }
58+
59+ uniform(bytes: number): GPUBuffer {
60+ return this.device.createBuffer({ size: Math.max(16, Math.ceil(bytes / 16) * 16), usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
61+ }
62+
63+ write(buf: GPUBuffer, data: ArrayBufferView, offset = 0) {
64+ this.device.queue.writeBuffer(buf, offset, data.buffer, data.byteOffset, data.byteLength);
65+ }
66+
67+ /** Compile a pipeline once; throws with the WGSL diagnostics if the shader does not compile. */
68+ async pipeline(key: string, code: string, entry: string): Promise<GPUComputePipeline> {
69+ let p = this.pipelines.get(key);
70+ if (!p) {
71+ const module = this.device.createShaderModule({ code, label: key });
72+ try {
73+ p = await this.device.createComputePipelineAsync({ layout: 'auto', compute: { module, entryPoint: entry }, label: key });
74+ } catch (e) {
75+ throw new Error(`could not compile the ${key} kernel: ${e instanceof Error ? e.message : String(e)}`);
76+ }
77+ this.pipelines.set(key, p);
78+ }
79+ return p;
80+ }
81+
82+ bind(pipeline: GPUComputePipeline, buffers: GPUBuffer[]): GPUBindGroup {
83+ return this.device.createBindGroup({
84+ layout: pipeline.getBindGroupLayout(0),
85+ entries: buffers.map((b, i) => ({ binding: i, resource: { buffer: b } })),
86+ });
87+ }
88+
89+ dispatch(pipeline: GPUComputePipeline, group: GPUBindGroup, nGroups: number) {
90+ const enc = this.device.createCommandEncoder();
91+ const pass = enc.beginComputePass();
92+ pass.setPipeline(pipeline);
93+ pass.setBindGroup(0, group);
94+ pass.dispatchWorkgroups(nGroups);
95+ pass.end();
96+ this.device.queue.submit([enc.finish()]);
97+ }
98+
99+ clear(buf: GPUBuffer) {
100+ const enc = this.device.createCommandEncoder();
101+ enc.clearBuffer(buf);
102+ this.device.queue.submit([enc.finish()]);
103+ }
104+
105+ async read(buf: GPUBuffer, bytes = buf.size): Promise<ArrayBuffer> {
106+ const staging = this.device.createBuffer({ size: bytes, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
107+ const enc = this.device.createCommandEncoder();
108+ enc.copyBufferToBuffer(buf, 0, staging, 0, bytes);
109+ this.device.queue.submit([enc.finish()]);
110+ await staging.mapAsync(GPUMapMode.READ);
111+ const copy = staging.getMappedRange().slice(0);
112+ staging.unmap();
113+ staging.destroy();
114+ return copy;
115+ }
116+
117+ async sync() {
118+ await this.device.queue.onSubmittedWorkDone();
119+ }
120+}
121+
122+/** Little-endian writer for uniform structs of 32-bit scalars. */
123+export class UniformWriter {
124+ private view: DataView;
125+ private off = 0;
126+ readonly bytes: ArrayBuffer;
127+
128+ constructor(nWords: number) {
129+ this.bytes = new ArrayBuffer(Math.ceil((4 * nWords) / 16) * 16);
130+ this.view = new DataView(this.bytes);
131+ }
132+ u32(v: number) {
133+ this.view.setUint32(this.off, v >>> 0, true);
134+ this.off += 4;
135+ return this;
136+ }
137+ i32(v: number) {
138+ this.view.setInt32(this.off, v | 0, true);
139+ this.off += 4;
140+ return this;
141+ }
142+ f32(v: number) {
143+ this.view.setFloat32(this.off, v, true);
144+ this.off += 4;
145+ return this;
146+ }
147+ /** Skip to the next 16-byte boundary (vec4 alignment). */
148+ align16() {
149+ this.off = Math.ceil(this.off / 16) * 16;
150+ return this;
151+ }
152+}
153+
154+/**
155+ * Chooses how much work to submit per dispatch so that each one takes about
156+ * `target` seconds: long enough to keep the GPU busy, short enough that the
157+ * page stays responsive, progress can be reported, and no driver watchdog is
158+ * tripped. Work is measured in arbitrary units chosen by the caller.
159+ */
160+export class Pacer {
161+ private rate = 0; // units per second, 0 = unknown
162+
163+ constructor(
164+ private target = 0.15,
165+ private initial = 1,
166+ ) {}
167+
168+ /** Units to submit next, given a hard cap. */
169+ next(cap: number): number {
170+ if (this.rate <= 0) return Math.min(cap, this.initial);
171+ return Math.max(1, Math.min(cap, Math.floor(this.rate * this.target)));
172+ }
173+
174+ /** Record that `units` took `seconds`. */
175+ record(units: number, seconds: number) {
176+ if (seconds <= 0) return;
177+ const r = units / Math.max(seconds, 1e-3);
178+ this.rate = this.rate > 0 ? 0.5 * (this.rate + r) : r;
179+ }
180+}
src/gpu/operator.tsadded+635−0View file
@@ -0,0 +1,635 @@
1+/**
2+ * The discretized operator of fmm3dbie's right-preconditioned combined field
3+ * impedance equation (helm_rpcomb_imp), on WebGPU. A port of
4+ * fmm3dbie_gpu/solver.py from the CUDA proof of concept.
5+ *
6+ * u = S_k[sigma] + i alpha D_k S_kappa[sigma], kappa = i|k|,
7+ *
8+ * with the boundary condition du/dn + i k lambda u = f. The discretized
9+ * equation is
10+ *
11+ * z sigma + S'_k sigma + i alpha S'_kappa phi2 + i alpha (D'_k - D'_kappa) phi1
12+ * + i k lambda [S_k sigma + i alpha D_k phi1 + i alpha w phi1] = f,
13+ *
14+ * phi1 = S_kappa sigma, phi2 = S'_kappa sigma,
15+ * z = (-1)^ifinout / 2 - i alpha / 4, w = 1/2 (exterior), -1/2 (interior).
16+ *
17+ * Every layer potential is the sum of a far part (smooth rule on the
18+ * oversampled nodes of the patches that are not near the target, the dense
19+ * kernel) and a near part (precomputed weights for the near patches, the
20+ * sparse kernel). The near weights are computed once, in chunks of pairs, and
21+ * kept on the GPU. GMRES runs on the CPU in double precision and calls
22+ * applyK for each Krylov vector.
23+ */
24+import { leggauss } from '../gauss.ts';
25+import { oversampleMatrix, type Surface } from '../surface.ts';
26+import type { CVec } from '../problem.ts';
27+import { gmres, type GmresResult } from '../gmres.ts';
28+import { Gpu, Pacer, UniformWriter } from './device.ts';
29+import { BUILD_LEAVES, DENSE_APPLY, DENSE_WG, NEAR_WEIGHTS, NT, SPARSE_APPLY } from './shaders.ts';
30+
31+export interface OperatorOptions {
32+ /** near radius factor: a patch is near a target within rfac * (patch radius) of its centroid */
33+ rfac?: number;
34+ rfacQ?: number;
35+ nq?: number;
36+ nqr?: number;
37+ nqx?: number;
38+ subLen?: number;
39+ oscLen?: number;
40+ kRMax?: number;
41+ nlmax?: number;
42+ nover?: number;
43+ /** seconds of GPU work per dispatch */
44+ pace?: number;
45+ progress?: (stage: string, fraction: number) => void;
46+}
47+
48+export interface NearPairs {
49+ pairT: Int32Array;
50+ pairP: Int32Array;
51+ rowPtr: Uint32Array;
52+ /** per pair: (u_t, v_t, kind, 0); kind 0 off-surface target, 1 on-surface, 2 self patch */
53+ pairF: Float32Array;
54+ npairs: number;
55+}
56+
57+/** A set of targets with everything needed to apply the operator at them. */
58+interface TargetSet {
59+ n: number;
60+ targ: GPUBuffer; // 2 vec4 per target: [x, 0], [n, 0]
61+ pairs: NearPairs;
62+ pairsI: GPUBuffer;
63+ pairsF: GPUBuffer;
64+ rowPtr: GPUBuffer;
65+ pairP: GPUBuffer;
66+ W: GPUBuffer;
67+ nleaves: number;
68+}
69+
70+export interface SetupInfo {
71+ npts: number;
72+ noversampled: number;
73+ npairs: number;
74+ nleaves: number;
75+ weightsBytes: number;
76+ setupSeconds: number;
77+}
78+
79+const now = () => performance.now() / 1000;
80+
81+export class Operator {
82+ readonly k: number;
83+ readonly alpha: number;
84+ /** wall impedance; the near weights do not depend on it, so it may be changed between solves */
85+ lam: number;
86+ readonly ifinout: number;
87+ readonly norder: number;
88+ readonly npols: number;
89+ readonly npts: number;
90+ readonly nover: number;
91+ readonly npolsO: number;
92+ readonly rfac: number;
93+ readonly quad: { nq: number; nqr: number; nqx: number; subLen: number; oscLen: number };
94+ readonly leaf: { rfacQ: number; nlmax: number; kRMax: number };
95+ readonly progress: (stage: string, fraction: number) => void;
96+ readonly pace: number;
97+ info!: SetupInfo;
98+
99+ private gpu: Gpu;
100+ private surf: Surface;
101+ private r2: Float64Array;
102+ private overMat: Float64Array;
103+ private gpar!: GPUBuffer;
104+ private tables!: GPUBuffer;
105+ private srcgeo!: GPUBuffer;
106+ private srcq!: GPUBuffer;
107+ private qbuf!: GPUBuffer;
108+ private out!: GPUBuffer;
109+ private wall!: TargetSet;
110+ private nself = 0;
111+ private selfItems: Float32Array[] = []; // per surface node, its polar segments (8 floats each)
112+ private uniformLeaves: GPUBuffer;
113+ private uniformWeights: GPUBuffer;
114+ private uniformDense: GPUBuffer;
115+ private uniformSparse: GPUBuffer;
116+ private destroyed = false;
117+
118+ constructor(gpu: Gpu, surf: Surface, k: number, alpha: number, lam: number, ifinout = 0, opts: OperatorOptions = {}) {
119+ this.gpu = gpu;
120+ this.surf = surf;
121+ this.k = k;
122+ this.alpha = alpha;
123+ this.lam = lam;
124+ this.ifinout = ifinout;
125+ this.norder = surf.norder;
126+ this.npols = surf.npols;
127+ this.npts = surf.npts;
128+ if (this.norder > 10) throw new Error('order at most 10');
129+ this.rfac = opts.rfac ?? (this.norder <= 6 ? 2.0 : 1.25); // fmm3dbie's get_rfacs
130+ this.nover = opts.nover ?? this.norder + 2;
131+ this.npolsO = (this.nover + 1) ** 2;
132+ let radMax = 0;
133+ for (let j = 0; j < surf.npatches; j++) radMax = Math.max(radMax, surf.rads[j]);
134+ // the radial rule must resolve the basis and exp(ikr) along a ray
135+ const nqr = opts.nqr ?? Math.floor(Math.min(Math.max(12, this.norder + 2 * Math.abs(k) * radMax + 6), 32));
136+ this.quad = { nq: opts.nq ?? 8, nqr, nqx: opts.nqx ?? 8, subLen: opts.subLen ?? 1.0, oscLen: opts.oscLen ?? 2.0 };
137+ this.leaf = { rfacQ: opts.rfacQ ?? 2.0, nlmax: opts.nlmax ?? 14, kRMax: opts.kRMax ?? 3.0 };
138+ this.progress = opts.progress ?? (() => {});
139+ this.pace = opts.pace ?? 0.15;
140+ this.r2 = new Float64Array(surf.npatches);
141+ for (let j = 0; j < surf.npatches; j++) this.r2[j] = (this.rfac * surf.rads[j]) ** 2;
142+ this.overMat = oversampleMatrix(this.norder, this.nover);
143+ this.uniformLeaves = gpu.uniform(64);
144+ this.uniformWeights = gpu.uniform(64);
145+ this.uniformDense = gpu.uniform(64);
146+ this.uniformSparse = gpu.uniform(64);
147+ }
148+
149+ /** Upload the geometry and compute the near weights for the surface nodes. */
150+ async setup(): Promise<SetupInfo> {
151+ const t0 = now();
152+ const { gpu, surf, npols, npts } = this;
153+ // patch charts: 4 vec4 per patch, [n, s0], [e1, s1], [e2, t0], [0, t1]
154+ const gp = new Float32Array(16 * surf.npatches);
155+ for (let j = 0; j < surf.npatches; j++) {
156+ const par = surf.exact.params.subarray(13 * j, 13 * j + 13);
157+ gp.set([par[0], par[1], par[2], par[9], par[3], par[4], par[5], par[10], par[6], par[7], par[8], par[11], 0, 0, 0, par[12]], 16 * j);
158+ }
159+ this.gpar = gpu.storage(gp);
160+ // umat and the quadrature rules in one table
161+ const { nq, nqr, nqx } = this.quad;
162+ const gl = leggauss(nq);
163+ const glr = leggauss(nqr);
164+ const glx = leggauss(nqx);
165+ const tab = new Float32Array(npols * npols + 2 * nq + 2 * nqr + 2 * nqx);
166+ tab.set(surf.umat, 0);
167+ let o = npols * npols;
168+ tab.set(gl.x, o);
169+ tab.set(gl.w, o + nq);
170+ o += 2 * nq;
171+ for (let i = 0; i < nqr; i++) {
172+ tab[o + i] = 0.5 * (glr.x[i] + 1);
173+ tab[o + nqr + i] = 0.5 * glr.w[i];
174+ }
175+ o += 2 * nqr;
176+ tab.set(glx.x, o);
177+ tab.set(glx.w, o + nqx);
178+ this.tables = gpu.storage(tab);
179+ // oversampled sources for the far field
180+ const ov = surf.oversampled(this.nover);
181+ const N = ov.w.length;
182+ const sg = new Float32Array(8 * N);
183+ for (let i = 0; i < N; i++) {
184+ sg[8 * i] = ov.x[i];
185+ sg[8 * i + 1] = ov.x[N + i];
186+ sg[8 * i + 2] = ov.x[2 * N + i];
187+ sg[8 * i + 3] = ov.w[i];
188+ sg[8 * i + 4] = ov.n[i];
189+ sg[8 * i + 5] = ov.n[N + i];
190+ sg[8 * i + 6] = ov.n[2 * N + i];
191+ }
192+ this.srcgeo = gpu.storage(sg);
193+ this.srcq = gpu.storage(16 * N);
194+ this.qbuf = gpu.storage(2 * 8 * npts);
195+ this.out = gpu.storage(6 * 8 * npts);
196+
197+ // self-patch polar segments for every node, from the first fundamental form at the node
198+ this.selfItems = new Array(npts);
199+ for (let i = 0; i < npts; i++) this.selfItems[i] = this.selfSegments(i);
200+
201+ const xs = surf.srcvals.subarray(0, 3 * npts);
202+ const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
203+ this.wall = await this.makeTargets(xs, nrm, true, 'near quadrature');
204+ await gpu.sync();
205+ const setupSeconds = now() - t0;
206+ this.info = {
207+ npts,
208+ noversampled: N,
209+ npairs: this.wall.pairs.npairs,
210+ nleaves: this.wall.nleaves,
211+ weightsBytes: this.wall.W.size,
212+ setupSeconds,
213+ };
214+ return this.info;
215+ }
216+
217+ destroy() {
218+ this.destroyed = true;
219+ for (const b of [this.gpar, this.tables, this.srcgeo, this.srcq, this.qbuf, this.out, this.uniformLeaves, this.uniformWeights, this.uniformDense, this.uniformSparse]) b?.destroy();
220+ if (this.wall) this.freeTargets(this.wall);
221+ }
222+
223+ /** Bytes of near weights for npairs pairs. */
224+ weightsBytes(npairs: number): number {
225+ return npairs * 6 * this.npols * 8;
226+ }
227+
228+ /** Estimated bytes of near weights for the surface nodes, before setup. */
229+ estimateWeightsBytes(): number {
230+ const { npts } = this;
231+ const xs = this.surf.srcvals.subarray(0, 3 * npts);
232+ let cnt = 0;
233+ const { npatches } = this.surf;
234+ const cen = this.surf.centroids;
235+ // sample every 16th node
236+ for (let i = 0; i < npts; i += 16) {
237+ for (let j = 0; j < npatches; j++) {
238+ const d2 = (xs[i] - cen[j]) ** 2 + (xs[npts + i] - cen[npatches + j]) ** 2 + (xs[2 * npts + i] - cen[2 * npatches + j]) ** 2;
239+ if (d2 <= this.r2[j]) cnt++;
240+ }
241+ }
242+ return this.weightsBytes(Math.ceil((cnt * npts) / Math.ceil(npts / 16)));
243+ }
244+
245+ // ---------------------------------------------------------------- near pairs
246+
247+ /** Pairs (target, patch) with |x_t - c_j|^2 <= r2_j, grouped by target, patches ascending. */
248+ nearPairs(x: Float64Array, onSurface: boolean): NearPairs {
249+ const n = x.length / 3;
250+ const { npatches, centroids: cen, patchId, uvNodes } = this.surf;
251+ const r2 = this.r2;
252+ const tList: number[] = [];
253+ const pList: number[] = [];
254+ const fList: number[] = [];
255+ const rowPtr = new Uint32Array(n + 1);
256+ const params = this.surf.exact.params;
257+ for (let i = 0; i < n; i++) {
258+ const xi = x[i];
259+ const yi = x[n + i];
260+ const zi = x[2 * n + i];
261+ for (let j = 0; j < npatches; j++) {
262+ const d2 = (xi - cen[j]) ** 2 + (yi - cen[npatches + j]) ** 2 + (zi - cen[2 * npatches + j]) ** 2;
263+ if (d2 > r2[j]) continue;
264+ tList.push(i);
265+ pList.push(j);
266+ if (!onSurface) {
267+ fList.push(0, 0, 0, 0);
268+ } else if (patchId[i] === j) {
269+ fList.push(uvNodes[i], uvNodes[n + i], 2, 0);
270+ } else {
271+ // the target mapped into this patch's face chart: c_t = x_t / (x_t . n)
272+ const par = params.subarray(13 * j, 13 * j + 13);
273+ const xn = xi * par[0] + yi * par[1] + zi * par[2];
274+ if (xn < 0.2) {
275+ fList.push(0, 0, 0, 0);
276+ } else {
277+ const cx = xi / xn;
278+ const cy = yi / xn;
279+ const cz = zi / xn;
280+ const st = cx * par[3] + cy * par[4] + cz * par[5];
281+ const tt = cx * par[6] + cy * par[7] + cz * par[8];
282+ fList.push((2 * (st - par[9])) / (par[10] - par[9]) - 1, (2 * (tt - par[11])) / (par[12] - par[11]) - 1, 1, 0);
283+ }
284+ }
285+ }
286+ rowPtr[i + 1] = tList.length;
287+ }
288+ return { pairT: Int32Array.from(tList), pairP: Int32Array.from(pList), rowPtr, pairF: Float32Array.from(fList), npairs: tList.length };
289+ }
290+
291+ /** The sinh-substituted segments along the four edges for the self patch of node i: 8 floats each. */
292+ private selfSegments(i: number): Float32Array {
293+ const { npts } = this;
294+ const sv = this.surf.srcvals;
295+ const u0 = this.surf.uvNodes[i];
296+ const v0 = this.surf.uvNodes[npts + i];
297+ const ru = [sv[3 * npts + i], sv[4 * npts + i], sv[5 * npts + i]];
298+ const rv = [sv[6 * npts + i], sv[7 * npts + i], sv[8 * npts + i]];
299+ const E = ru[0] ** 2 + ru[1] ** 2 + ru[2] ** 2;
300+ const F = ru[0] * rv[0] + ru[1] * rv[1] + ru[2] * rv[2];
301+ const G = rv[0] ** 2 + rv[1] ** 2 + rv[2] ** 2;
302+ const det = Math.sqrt(Math.max(E * G - F * F, 1e-300));
303+ const k = Math.abs(this.k);
304+ const { subLen, oscLen } = this.quad;
305+ const items: number[] = [];
306+ for (let e = 0; e < 4; e++) {
307+ const h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e];
308+ if (Math.abs(h) < 1e-12) continue;
309+ let x1: number, x2: number, xc: number, heff: number, arc: number;
310+ if (e < 2) {
311+ x1 = -1 - v0;
312+ x2 = 1 - v0;
313+ xc = (-h * F) / G;
314+ heff = (Math.abs(h) * det) / G;
315+ arc = Math.sqrt(G);
316+ } else {
317+ x1 = -1 - u0;
318+ x2 = 1 - u0;
319+ xc = (-h * F) / E;
320+ heff = (Math.abs(h) * det) / E;
321+ arc = Math.sqrt(E);
322+ }
323+ const xi0 = Math.asinh((x1 - xc) / heff);
324+ const xi1 = Math.asinh((x2 - xc) / heff);
325+ // sub-intervals of length <= subLen in xi whose extent along the edge is <= dxmax
326+ const dxmax = k * arc > 0 ? oscLen / (k * arc) : Infinity;
327+ let xi = xi0;
328+ let cnt = 0;
329+ while (xi < xi1 && cnt < 256) {
330+ const xh = xc + heff * Math.sinh(xi);
331+ let xn = Math.min(xi + subLen, Math.asinh((xh + dxmax - xc) / heff));
332+ if (xn > xi1 || cnt === 255) xn = xi1;
333+ items.push(1, e, xi, xn, h, xc, heff, 0);
334+ xi = xn;
335+ cnt++;
336+ }
337+ }
338+ return Float32Array.from(items);
339+ }
340+
341+ // ------------------------------------------------------------ near weights
342+
343+ /** Build the near pairs, quadtree leaves and near weights for a set of targets. */
344+ private async makeTargets(x: Float64Array, nrm: Float64Array | null, onSurface: boolean, stage: string): Promise<TargetSet> {
345+ const { gpu } = this;
346+ const n = x.length / 3;
347+ const tg = new Float32Array(8 * n);
348+ for (let i = 0; i < n; i++) {
349+ tg[8 * i] = x[i];
350+ tg[8 * i + 1] = x[n + i];
351+ tg[8 * i + 2] = x[2 * n + i];
352+ if (nrm) {
353+ tg[8 * i + 4] = nrm[i];
354+ tg[8 * i + 5] = nrm[n + i];
355+ tg[8 * i + 6] = nrm[2 * n + i];
356+ }
357+ }
358+ const pairs = this.nearPairs(x, onSurface);
359+ const wbytes = this.weightsBytes(pairs.npairs);
360+ const lim = gpu.device.limits;
361+ if (wbytes > lim.maxStorageBufferBindingSize || wbytes > lim.maxBufferSize) {
362+ throw new Error(
363+ `the near weights need ${(wbytes / 2 ** 20).toFixed(0)} MB but this GPU allows buffers of at most ` +
364+ `${(Math.min(lim.maxStorageBufferBindingSize, lim.maxBufferSize) / 2 ** 20).toFixed(0)} MB; use fewer patches or a lower order`,
365+ );
366+ }
367+ const pi = new Int32Array(2 * Math.max(1, pairs.npairs));
368+ for (let p = 0; p < pairs.npairs; p++) {
369+ pi[2 * p] = pairs.pairT[p];
370+ pi[2 * p + 1] = pairs.pairP[p];
371+ }
372+ const set: TargetSet = {
373+ n,
374+ targ: gpu.storage(tg),
375+ pairs,
376+ pairsI: gpu.storage(pi),
377+ pairsF: gpu.storage(pairs.npairs ? pairs.pairF : new Float32Array(4)),
378+ rowPtr: gpu.storage(pairs.rowPtr),
379+ pairP: gpu.storage(pairs.npairs ? new Uint32Array(pairs.pairP.buffer, pairs.pairP.byteOffset, pairs.npairs) : new Uint32Array(1)),
380+ W: gpu.storage(Math.max(16, wbytes)),
381+ nleaves: 0,
382+ };
383+ if (pairs.npairs > 0) set.nleaves = await this.computeWeights(set, stage);
384+ return set;
385+ }
386+
387+ private freeTargets(s: TargetSet) {
388+ for (const b of [s.targ, s.pairsI, s.pairsF, s.rowPtr, s.pairP, s.W]) b.destroy();
389+ }
390+
391+ /** Quadtree leaves and near weights for all pairs of a target set, in paced chunks. Returns the leaf count. */
392+ private async computeWeights(set: TargetSet, stage: string): Promise<number> {
393+ const { gpu, quad, leaf } = this;
394+ const { npairs, pairF } = set.pairs;
395+ const leavesPipe = await gpu.pipeline('leaves', BUILD_LEAVES, 'build_leaves');
396+ const weightsPipe = await gpu.pipeline('weights', NEAR_WEIGHTS, 'near_weights');
397+ const pacer = new Pacer(this.pace, 512);
398+ let nleaves = 0;
399+ let p0 = 0;
400+ let itemsPerPair = 40;
401+ const maxItemsBytes = 48 * 2 ** 20;
402+ while (p0 < npairs) {
403+ if (this.destroyed) throw new Error('operator destroyed');
404+ const cap = Math.min(65535, npairs - p0, Math.max(256, Math.floor(maxItemsBytes / (32 * itemsPerPair))));
405+ const n = pacer.next(cap);
406+ const t0 = now();
407+ // pass 1: leaf counts
408+ const counts = gpu.storage(4 * n);
409+ const uw = new UniformWriter(9);
410+ uw.u32(n).u32(p0).u32(0).i32(leaf.nlmax).f32(leaf.rfacQ).f32(Math.abs(this.k)).f32(leaf.kRMax).f32(this.surf.exact.a).i32(this.surf.exact.p);
411+ gpu.write(this.uniformLeaves, new Uint8Array(uw.bytes));
412+ const dummyItems = gpu.storage(16);
413+ const dummyOffsets = gpu.storage(16);
414+ gpu.dispatch(leavesPipe, gpu.bind(leavesPipe, [this.uniformLeaves, set.pairsI, set.pairsF, set.targ, this.gpar, counts, dummyOffsets, dummyItems]), Math.ceil(n / 64));
415+ const cnt = new Uint32Array(await gpu.read(counts, 4 * n));
416+ // self pairs: the CPU segments
417+ const ptr = new Uint32Array(n + 1);
418+ for (let q = 0; q < n; q++) {
419+ const p = p0 + q;
420+ const c = pairF[4 * p + 2] === 2 ? this.selfItems[set.pairs.pairT[p]].length / 8 : cnt[q];
421+ if (pairF[4 * p + 2] !== 2) nleaves += cnt[q];
422+ ptr[q + 1] = ptr[q] + c;
423+ }
424+ const nitems = ptr[n];
425+ const items = new Float32Array(8 * Math.max(1, nitems));
426+ for (let q = 0; q < n; q++) {
427+ const p = p0 + q;
428+ if (pairF[4 * p + 2] === 2) items.set(this.selfItems[set.pairs.pairT[p]], 8 * ptr[q]);
429+ }
430+ const itemsBuf = gpu.storage(items);
431+ const ptrBuf = gpu.storage(ptr);
432+ // pass 2: fill the leaves
433+ const uw2 = new UniformWriter(9);
434+ uw2.u32(n).u32(p0).u32(1).i32(leaf.nlmax).f32(leaf.rfacQ).f32(Math.abs(this.k)).f32(leaf.kRMax).f32(this.surf.exact.a).i32(this.surf.exact.p);
435+ gpu.write(this.uniformLeaves, new Uint8Array(uw2.bytes));
436+ gpu.dispatch(leavesPipe, gpu.bind(leavesPipe, [this.uniformLeaves, set.pairsI, set.pairsF, set.targ, this.gpar, counts, ptrBuf, itemsBuf]), Math.ceil(n / 64));
437+ // the weights
438+ const uw3 = new UniformWriter(10);
439+ uw3.u32(n).u32(p0).i32(this.norder).i32(this.npols).i32(quad.nq).i32(quad.nqr).i32(quad.nqx).f32(this.k).f32(this.surf.exact.a).i32(this.surf.exact.p);
440+ gpu.write(this.uniformWeights, new Uint8Array(uw3.bytes));
441+ gpu.dispatch(weightsPipe, gpu.bind(weightsPipe, [this.uniformWeights, set.pairsI, set.pairsF, set.targ, this.gpar, ptrBuf, itemsBuf, this.tables, set.W]), n);
442+ await gpu.sync();
443+ const dt = now() - t0;
444+ pacer.record(n, dt);
445+ itemsPerPair = Math.max(1, nitems / n);
446+ for (const b of [counts, dummyItems, dummyOffsets, itemsBuf, ptrBuf]) b.destroy();
447+ p0 += n;
448+ this.progress(stage, p0 / npairs);
449+ }
450+ return nleaves;
451+ }
452+
453+ // ------------------------------------------------------------- application
454+
455+ /** Values of q at the oversampled nodes, per patch: q_o = M q. */
456+ private oversample(q: CVec, out: Float32Array, offset: number) {
457+ const { npols, npolsO, overMat } = this;
458+ const np = this.surf.npatches;
459+ for (let j = 0; j < np; j++) {
460+ for (let kk = 0; kk < npolsO; kk++) {
461+ let sr = 0;
462+ let si = 0;
463+ for (let m = 0; m < npols; m++) {
464+ const M = overMat[kk * npols + m];
465+ sr += M * q[2 * (j * npols + m)];
466+ si += M * q[2 * (j * npols + m) + 1];
467+ }
468+ const idx = 4 * (j * npolsO + kk) + offset;
469+ out[idx] = sr;
470+ out[idx + 1] = si;
471+ }
472+ }
473+ }
474+
475+ /** Far and near parts of the selected kernels at a target set. sel: 6 entries in {0, 1, 2}; reqs: [kernel, density] pairs. */
476+ private async applyAt(set: TargetSet, qa: CVec, qb: CVec | null, sel: number[], reqs: [number, number][]): Promise<Float32Array> {
477+ const { gpu, npts, npolsO } = this;
478+ const N = npolsO * this.surf.npatches;
479+ const sq = new Float32Array(4 * N);
480+ this.oversample(qa, sq, 0);
481+ if (qb) this.oversample(qb, sq, 2);
482+ gpu.write(this.srcq, sq);
483+ const qd = new Float32Array(2 * 2 * npts);
484+ for (let i = 0; i < 2 * npts; i++) qd[i] = qa[i];
485+ if (qb) for (let i = 0; i < 2 * npts; i++) qd[2 * npts + i] = qb[i];
486+ gpu.write(this.qbuf, qd);
487+ const out = set.n === this.wall.n ? this.out : gpu.storage(6 * 8 * set.n);
488+ gpu.clear(out);
489+ const selmask = sel.reduce((acc, s, m) => acc | (s << (4 * m)), 0);
490+ const densePipe = await gpu.pipeline('dense', DENSE_APPLY, 'dense_apply');
491+ const group = gpu.bind(densePipe, [this.uniformDense, set.targ, this.srcgeo, this.srcq, set.rowPtr, set.pairP, out]);
492+ const pacer = new Pacer(this.pace, Math.max(1, Math.floor(4e8 / (set.n * npolsO * 60))));
493+ const np = this.surf.npatches;
494+ let j0 = 0;
495+ while (j0 < np) {
496+ const n = pacer.next(np - j0);
497+ const t0 = now();
498+ const uw = new UniformWriter(6);
499+ uw.u32(set.n).u32(j0).u32(j0 + n).u32(npolsO).f32(this.k).u32(selmask);
500+ gpu.write(this.uniformDense, new Uint8Array(uw.bytes));
501+ gpu.dispatch(densePipe, group, Math.ceil(set.n / DENSE_WG));
502+ await gpu.sync();
503+ pacer.record(n, now() - t0);
504+ j0 += n;
505+ }
506+ if (set.pairs.npairs > 0 && reqs.length > 0) {
507+ const sparsePipe = await gpu.pipeline('sparse', SPARSE_APPLY, 'sparse_apply');
508+ const uw = new UniformWriter(8);
509+ uw.u32(set.n).u32(this.npols).u32(npts).u32(reqs.length);
510+ for (let r = 0; r < 4; r++) uw.u32(r < reqs.length ? reqs[r][0] | (reqs[r][1] << 8) : 0);
511+ gpu.write(this.uniformSparse, new Uint8Array(uw.bytes));
512+ gpu.dispatch(sparsePipe, gpu.bind(sparsePipe, [this.uniformSparse, set.rowPtr, set.pairP, set.W, this.qbuf, out]), Math.ceil(set.n / 128));
513+ }
514+ const res = new Float32Array(await gpu.read(out, 6 * 8 * set.n));
515+ if (out !== this.out) out.destroy();
516+ return res;
517+ }
518+
519+ zid(): [number, number] {
520+ return [(-1) ** this.ifinout / 2, -this.alpha / 4];
521+ }
522+
523+ /** The operator without the identity term z. */
524+ async applyK(sigma: CVec): Promise<CVec> {
525+ const { npts, k, alpha: a, lam } = this;
526+ const F = await this.applyAt(this.wall, sigma, null, [1, 1, 1, 0, 1, 0], [[0, 0], [1, 0], [2, 0], [4, 0]]);
527+ const phi1 = new Float64Array(2 * npts);
528+ const phi2 = new Float64Array(2 * npts);
529+ const res = new Float64Array(2 * npts);
530+ for (let i = 0; i < npts; i++) {
531+ phi1[2 * i] = F[2 * (npts + i)];
532+ phi1[2 * i + 1] = F[2 * (npts + i) + 1];
533+ phi2[2 * i] = F[2 * (2 * npts + i)];
534+ phi2[2 * i + 1] = F[2 * (2 * npts + i) + 1];
535+ const s4r = F[2 * (4 * npts + i)];
536+ const s4i = F[2 * (4 * npts + i) + 1];
537+ // res = F0 + i k lam F4
538+ res[2 * i] = F[2 * i] - k * lam * s4i;
539+ res[2 * i + 1] = F[2 * i + 1] + k * lam * s4r;
540+ }
541+ const G = await this.applyAt(this.wall, phi1, phi2, [0, 0, 2, 1, 0, 1], [[2, 1], [3, 0], [5, 0]]);
542+ const sgn = (-1) ** this.ifinout;
543+ for (let i = 0; i < npts; i++) {
544+ const g2r = G[2 * (2 * npts + i)];
545+ const g2i = G[2 * (2 * npts + i) + 1];
546+ const g3r = G[2 * (3 * npts + i)];
547+ const g3i = G[2 * (3 * npts + i) + 1];
548+ const g5r = G[2 * (5 * npts + i)];
549+ const g5i = G[2 * (5 * npts + i) + 1];
550+ // res += i a (G2 + G3) - k a lam G5 + sgn k a lam phi1 / 2
551+ res[2 * i] += -a * (g2i + g3i) - k * a * lam * g5r + (sgn * k * a * lam * phi1[2 * i]) / 2;
552+ res[2 * i + 1] += a * (g2r + g3r) - k * a * lam * g5i + (sgn * k * a * lam * phi1[2 * i + 1]) / 2;
553+ }
554+ return res;
555+ }
556+
557+ /** phi1 = S_kappa[sigma] on the surface. */
558+ async sKappa(sigma: CVec): Promise<CVec> {
559+ const F = await this.applyAt(this.wall, sigma, null, [0, 1, 0, 0, 0, 0], [[1, 0]]);
560+ const { npts } = this;
561+ const phi1 = new Float64Array(2 * npts);
562+ for (let i = 0; i < 2 * npts; i++) phi1[i] = F[2 * npts + i];
563+ return phi1;
564+ }
565+
566+ /** One of the six operators applied to q on the surface (principal value), for tests. */
567+ async applySingle(m: number, q: CVec): Promise<CVec> {
568+ const sel = [0, 0, 0, 0, 0, 0];
569+ sel[m] = 1;
570+ const F = await this.applyAt(this.wall, q, null, sel, [[m, 0]]);
571+ const { npts } = this;
572+ const out = new Float64Array(2 * npts);
573+ for (let i = 0; i < 2 * npts; i++) out[i] = F[2 * m * npts + i];
574+ return out;
575+ }
576+
577+ async solve(rhs: CVec, epsGmres = 1e-5, numit = 200, onIter?: (it: number, rres: number) => void): Promise<GmresResult & { phi1: CVec; seconds: number }> {
578+ const t0 = now();
579+ const r = await gmres((x) => this.applyK(x), rhs, this.surf.wts, this.zid(), epsGmres, numit, onIter);
580+ const phi1 = await this.sKappa(r.x);
581+ return { ...r, phi1, seconds: now() - t0 };
582+ }
583+
584+ /** u = S_k sigma + i alpha D_k phi1 at the surface nodes, as the limit from the side selected by ifinout. */
585+ async evalOnSurface(sigma: CVec, phi1: CVec): Promise<CVec> {
586+ const F = await this.applyAt(this.wall, sigma, phi1, [0, 0, 0, 0, 1, 2], [[4, 0], [5, 1]]);
587+ return this.combine(F, this.npts, phi1);
588+ }
589+
590+ /** u = S_k sigma + i alpha D_k phi1 at off-surface targets x (3, n), in chunks. */
591+ async eval(x: Float64Array, sigma: CVec, phi1: CVec, stage = 'field'): Promise<CVec> {
592+ const n = x.length / 3;
593+ const u = new Float64Array(2 * n);
594+ const chunk = 4096;
595+ for (let i0 = 0; i0 < n; i0 += chunk) {
596+ const m = Math.min(chunk, n - i0);
597+ const xs = new Float64Array(3 * m);
598+ for (let d = 0; d < 3; d++) xs.set(x.subarray(d * n + i0, d * n + i0 + m), d * m);
599+ const set = await this.makeTargets(xs, null, false, stage);
600+ const F = await this.applyAt(set, sigma, phi1, [0, 0, 0, 0, 1, 2], [[4, 0], [5, 1]]);
601+ this.freeTargets(set);
602+ u.set(this.combine(F, m, null), 2 * i0);
603+ this.progress(stage, (i0 + m) / n);
604+ }
605+ return u;
606+ }
607+
608+ private combine(F: Float32Array, n: number, phi1OnSurface: CVec | null): CVec {
609+ const a = this.alpha;
610+ const u = new Float64Array(2 * n);
611+ const jump = phi1OnSurface ? (a * (-1) ** (this.ifinout + 1)) / 2 : 0;
612+ for (let i = 0; i < n; i++) {
613+ const s4r = F[2 * (4 * n + i)];
614+ const s4i = F[2 * (4 * n + i) + 1];
615+ const d5r = F[2 * (5 * n + i)];
616+ const d5i = F[2 * (5 * n + i) + 1];
617+ // u = S + i a D (+ i a (-1)^(ifinout+1) phi1 / 2 on the surface)
618+ u[2 * i] = s4r - a * d5i;
619+ u[2 * i + 1] = s4i + a * d5r;
620+ if (phi1OnSurface) {
621+ u[2 * i] += -jump * phi1OnSurface[2 * i + 1];
622+ u[2 * i + 1] += jump * phi1OnSurface[2 * i];
623+ }
624+ }
625+ return u;
626+ }
627+
628+ /** The near weights of the surface nodes, for tests: (npairs, 6, npols) complex interleaved. */
629+ async readWallWeights(): Promise<{ pairs: NearPairs; W: Float32Array }> {
630+ const W = new Float32Array(await this.gpu.read(this.wall.W, this.weightsBytes(this.wall.pairs.npairs)));
631+ return { pairs: this.wall.pairs, W };
632+ }
633+}
634+
635+export { NT };
src/gpu/shaders.tsadded+555−0View file
@@ -0,0 +1,555 @@
1+/**
2+ * WGSL compute kernels for the Helmholtz layer potentials on quadrilateral
3+ * patch surfaces, in single precision. A port of fmm3dbie_gpu/cuda_ops.py
4+ * from the CUDA proof of concept, with two changes that single precision
5+ * requires:
6+ *
7+ * 1. For a target on the surface, the displacement y - x_t from the target to
8+ * a quadrature node is formed in the coordinates of the superellipsoid face
9+ * chart rather than as a difference of two rounded positions, so that it
10+ * carries full relative precision however small it is. Without this the
11+ * double-layer kernels, whose numerators (y - x_t).n are O(r^2), lose most
12+ * of their digits near the target.
13+ * 2. The weakly singular difference D'_k - D'_kappa is evaluated from a short
14+ * power series when k r < 1, where the direct formula subtracts two nearly
15+ * equal O(1/r^3) terms.
16+ *
17+ * Six kernels are handled together, in the order used by fmm3dbie's
18+ * helm_rpcomb_imp solver:
19+ *
20+ * 0: S'_k 1: S_kappa 2: S'_kappa 3: D'_k - D'_kappa 4: S_k 5: D_k
21+ *
22+ * with kappa = i|k| and k real, G_k = exp(ikr)/(4 pi r), S' the normal
23+ * derivative at the target, D the normal derivative at the source.
24+ */
25+
26+export const NT = 64; // threads per workgroup in near_weights
27+export const NORD1_MAX = 11; // norder <= 10
28+export const NPOLS_MAX = 121;
29+export const DENSE_WG = 128;
30+
31+const COMMON = /* wgsl */ `
32+const OVER4PI: f32 = 0.07957747154594767;
33+
34+fn cmul(a: vec2f, b: vec2f) -> vec2f {
35+ return vec2f(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
36+}
37+
38+// (cos x, sin x) to a few ulp. WGSL only asks sin and cos for an absolute error
39+// of 2^-11 on [-pi, pi], and some implementations (SwiftShader among them) use
40+// approximations that coarse, so we reduce by pi/2 in three parts (Cody-Waite)
41+// and use the minimax polynomials of Cephes on [-pi/4, pi/4].
42+fn sincos(x: f32) -> vec2f {
43+ let q = round(x * 0.63661977236758134);
44+ var r = x - q * 1.5703125;
45+ r = r - q * 4.837512969970703125e-4;
46+ r = r - q * 7.54978995489188216e-8;
47+ let z = r * r;
48+ let s = r + r * z * (-1.6666654611e-1 + z * (8.3321608736e-3 + z * (-1.9515295891e-4)));
49+ let c = 1.0 - 0.5 * z + z * z * (4.166664568298827e-2 + z * (-1.388731625493765e-3 + z * 2.443315711809948e-5));
50+ let iq = i32(q) & 3;
51+ switch (iq) {
52+ case 0: { return vec2f(c, s); }
53+ case 1: { return vec2f(-s, c); }
54+ case 2: { return vec2f(-c, -s); }
55+ default: { return vec2f(s, -c); }
56+ }
57+}
58+
59+fn ipow(x: f32, n: i32) -> f32 {
60+ var r = 1.0;
61+ for (var i = 0; i < n; i++) { r *= x; }
62+ return r;
63+}
64+
65+// sum_{j=0}^{m} c^j d^{m-j}, so that c^{m+1} - d^{m+1} = (c - d) * geomsum(c, d, m)
66+fn geomsum(c: f32, d: f32, m: i32) -> f32 {
67+ var acc = 1.0;
68+ var dm = 1.0;
69+ for (var j = 1; j <= m; j++) { dm *= d; acc = c * acc + dm; }
70+ return acc;
71+}
72+
73+// The six kernels at one (target, source) pair, times a real weight w.
74+// dx = x_t - y. K[3] is only computed when need3 is set.
75+fn kernels6(dx: vec3f, nt: vec3f, ns: vec3f, w: f32, k: f32, need3: bool) -> array<vec2f, 6> {
76+ let r2 = dot(dx, dx);
77+ let r = sqrt(r2);
78+ let rinv = 1.0 / r;
79+ let r2inv = rinv * rinv;
80+ let dnt = dot(dx, nt);
81+ let dns = dot(dx, ns);
82+ let kr = k * r;
83+ let E0 = OVER4PI * rinv * w;
84+ let ek = sincos(kr);
85+ let ec = exp(-kr);
86+ let Ek = E0 * ek;
87+ let Ec = E0 * ec;
88+ let fk = vec2f(1.0, -kr);
89+ let fc = 1.0 + kr;
90+ let fkEk = cmul(fk, Ek);
91+ var K: array<vec2f, 6>;
92+ K[4] = Ek;
93+ K[1] = vec2f(Ec, 0.0);
94+ K[0] = -dnt * r2inv * fkEk;
95+ K[2] = vec2f(-dnt * fc * Ec * r2inv, 0.0);
96+ K[5] = dns * r2inv * fkEk;
97+ K[3] = vec2f(0.0, 0.0);
98+ if (need3) {
99+ let nn = dot(ns, nt);
100+ var A: vec2f; // ((1 - ikr) e^{ikr} - (1 + kr) e^{-kr})
101+ var B: vec2f; // ((3 - 3ikr - k^2r^2) e^{ikr} - (3 + 3kr + k^2r^2) e^{-kr})
102+ if (kr < 1.0) {
103+ // both are sum_{n>=2} b_n (kr)^n (i^n - (-1)^n), with b_n = (1-n)/n! and (n-1)(n-3)/n!
104+ A = vec2f(0.0, 0.0);
105+ B = vec2f(0.0, 0.0);
106+ var t = 0.5 * kr * kr;
107+ for (var n = 2; n <= 14; n++) {
108+ var cn: vec2f;
109+ switch (n % 4) {
110+ case 0: { cn = vec2f(0.0, 0.0); }
111+ case 1: { cn = vec2f(1.0, 1.0); }
112+ case 2: { cn = vec2f(-2.0, 0.0); }
113+ default: { cn = vec2f(1.0, -1.0); }
114+ }
115+ let nf = f32(n);
116+ A += ((1.0 - nf) * t) * cn;
117+ B += ((nf - 1.0) * (nf - 3.0) * t) * cn;
118+ t *= kr / (nf + 1.0);
119+ }
120+ } else {
121+ A = cmul(fk, ek) - vec2f(fc * ec, 0.0);
122+ let gk = vec2f(3.0 - kr * kr, -3.0 * kr);
123+ let gc = 3.0 + 3.0 * kr + kr * kr;
124+ B = cmul(gk, ek) - vec2f(gc * ec, 0.0);
125+ }
126+ K[3] = E0 * (nn * r2inv * A - dns * dnt * r2inv * r2inv * B);
127+ }
128+ return K;
129+}
130+
131+fn legendre(x: f32, norder: i32, P: ptr<function, array<f32, ${NORD1_MAX}>>) {
132+ (*P)[0] = 1.0;
133+ if (norder >= 1) { (*P)[1] = x; }
134+ for (var n = 1; n < norder; n++) {
135+ (*P)[n + 1] = ((2.0 * f32(n) + 1.0) * x * (*P)[n] - f32(n) * (*P)[n - 1]) / f32(n + 1);
136+ }
137+}
138+`;
139+
140+/** The superellipsoid face chart; needs a gpar binding. */
141+const CHART = /* wgsl */ `
142+struct Geo { x: vec3f, ru: vec3f, rv: vec3f, d: vec3f }
143+
144+// xyz, r_u, r_v of the superellipsoid face chart of patch ip at (u, v), and the
145+// displacement d = y - x_t. With rel set, x_t is the chart point at (ut, vt) and d
146+// is formed in chart coordinates; otherwise d = y - xt for a given xt.
147+// gpar: 4 vec4 per patch, [n, s0], [e1, s1], [e2, t0], [0, t1].
148+fn chart_rel(ip: i32, u: f32, v: f32, ut: f32, vt: f32, rel: bool, xt: vec3f, ga: f32, gpe: i32) -> Geo {
149+ let g0 = gpar[4 * ip]; let g1 = gpar[4 * ip + 1]; let g2 = gpar[4 * ip + 2]; let g3 = gpar[4 * ip + 3];
150+ let nf = g0.xyz; let e1 = g1.xyz; let e2 = g2.xyz;
151+ let hs = 0.5 * (g1.w - g0.w); let ht = 0.5 * (g3.w - g2.w);
152+ let s = g0.w + (u + 1.0) * hs; let t = g2.w + (v + 1.0) * ht;
153+ let c = nf + s * e1 + t * e2;
154+ let S = ipow(c.x, gpe) + ipow(c.y, gpe) + ipow(c.z, gpe);
155+ let N = pow(S, 1.0 / f32(gpe));
156+ let gN = vec3f(ipow(c.x / N, gpe - 1), ipow(c.y / N, gpe - 1), ipow(c.z / N, gpe - 1));
157+ let dNu = dot(gN, e1); let dNv = dot(gN, e2);
158+ var out: Geo;
159+ out.x = ga * c / N;
160+ out.ru = ga * hs * (e1 / N - c * dNu / (N * N));
161+ out.rv = ga * ht * (e2 / N - c * dNv / (N * N));
162+ if (rel) {
163+ let st = g0.w + (ut + 1.0) * hs; let tt = g2.w + (vt + 1.0) * ht;
164+ let ct = nf + st * e1 + tt * e2;
165+ let dc = ((u - ut) * hs) * e1 + ((v - vt) * ht) * e2;
166+ let St = ipow(ct.x, gpe) + ipow(ct.y, gpe) + ipow(ct.z, gpe);
167+ let Nt = pow(St, 1.0 / f32(gpe));
168+ // N^p - Nt^p, then N - Nt, both without cancellation
169+ let dP = dc.x * geomsum(c.x, ct.x, gpe - 1) + dc.y * geomsum(c.y, ct.y, gpe - 1) + dc.z * geomsum(c.z, ct.z, gpe - 1);
170+ let dN = dP / geomsum(N, Nt, gpe - 1);
171+ out.d = ga * (dc * Nt - ct * dN) / (N * Nt);
172+ } else {
173+ out.d = out.x - xt;
174+ }
175+ return out;
176+}
177+
178+`;
179+
180+/**
181+ * Quadtree leaves on the reference square of a near (non-self) patch, refined
182+ * while the target is within rfac_q times the radius of a leaf or the leaf is
183+ * large compared with the wavelength. Two passes: counts only (fill = 0), then
184+ * fill using offsets. Self pairs are skipped; their items come from the CPU.
185+ * Items are 2 vec4: [kind, cu, cv, h], [0, 0, 0, 0] with kind 0 for a leaf.
186+ */
187+export const BUILD_LEAVES = /* wgsl */ `
188+struct Params {
189+ n: u32,
190+ pair_offset: u32,
191+ fill: u32,
192+ nlmax: i32,
193+ rfac_q: f32,
194+ zk: f32,
195+ kR_max: f32,
196+ ga: f32,
197+ gpe: i32,
198+}
199+@group(0) @binding(0) var<uniform> U: Params;
200+@group(0) @binding(1) var<storage, read> pairs_i: array<vec2i>;
201+@group(0) @binding(2) var<storage, read> pairs_f: array<vec4f>;
202+@group(0) @binding(3) var<storage, read> targ: array<vec4f>;
203+@group(0) @binding(4) var<storage, read> gpar: array<vec4f>;
204+@group(0) @binding(5) var<storage, read_write> counts: array<u32>;
205+@group(0) @binding(6) var<storage, read> offsets: array<u32>;
206+@group(0) @binding(7) var<storage, read_write> items: array<vec4f>;
207+${COMMON}
208+${CHART}
209+
210+@compute @workgroup_size(64)
211+fn build_leaves(@builtin(global_invocation_id) gid: vec3u) {
212+ let pl = gid.x;
213+ if (pl >= U.n) { return; }
214+ let p = pl + U.pair_offset;
215+ let ip = pairs_i[p].y;
216+ let it = pairs_i[p].x;
217+ let kind = i32(pairs_f[p].z);
218+ if (kind == 2) {
219+ if (U.fill == 0u) { counts[pl] = 0u; }
220+ return;
221+ }
222+ let tx = targ[2 * it].xyz;
223+ var su: array<f32, 64>;
224+ var sv: array<f32, 64>;
225+ var sh: array<f32, 64>;
226+ var sd: array<i32, 64>;
227+ var sp = 1;
228+ su[0] = 0.0; sv[0] = 0.0; sh[0] = 1.0; sd[0] = 0;
229+ var n = 0u;
230+ let rq2 = U.rfac_q * U.rfac_q;
231+ let kmax2 = U.kR_max * U.kR_max;
232+ let k2 = U.zk * U.zk;
233+ while (sp > 0) {
234+ sp--;
235+ let cu = su[sp]; let cv = sv[sp]; let h = sh[sp]; let d = sd[sp];
236+ var refine = false;
237+ if (d < U.nlmax) {
238+ // centroid and radius from a 3x3 sample grid
239+ var xc: array<vec3f, 9>;
240+ var cen = vec3f(0.0);
241+ for (var q = 0; q < 9; q++) {
242+ let g = chart_rel(ip, cu + h * f32(q % 3 - 1), cv + h * f32(q / 3 - 1), 0.0, 0.0, false, vec3f(0.0), U.ga, U.gpe);
243+ xc[q] = g.x;
244+ cen += g.x / 9.0;
245+ }
246+ var R2 = 0.0;
247+ for (var q = 0; q < 9; q++) { let e = xc[q] - cen; R2 = max(R2, dot(e, e)); }
248+ let e = tx - cen;
249+ let d2 = dot(e, e);
250+ refine = (d2 < rq2 * R2) || (k2 * R2 > kmax2);
251+ }
252+ if (refine && sp <= 60) {
253+ let hh = 0.5 * h;
254+ su[sp] = cu - hh; sv[sp] = cv - hh; sh[sp] = hh; sd[sp] = d + 1; sp++;
255+ su[sp] = cu + hh; sv[sp] = cv - hh; sh[sp] = hh; sd[sp] = d + 1; sp++;
256+ su[sp] = cu + hh; sv[sp] = cv + hh; sh[sp] = hh; sd[sp] = d + 1; sp++;
257+ su[sp] = cu - hh; sv[sp] = cv + hh; sh[sp] = hh; sd[sp] = d + 1; sp++;
258+ } else {
259+ if (U.fill != 0u) {
260+ let o = 2u * (offsets[pl] + n);
261+ items[o] = vec4f(0.0, cu, cv, h);
262+ items[o + 1u] = vec4f(0.0);
263+ }
264+ n++;
265+ }
266+ }
267+ if (U.fill == 0u) { counts[pl] = n; }
268+}
269+`;
270+
271+/**
272+ * Near-field weights, one workgroup per (target, patch) pair. The nodes of the
273+ * pair's items (self-patch polar segments or quadtree leaves) are evaluated
274+ * ${NT} at a time; each thread then accumulates the Legendre moments
275+ * I_n = sum K P_n |J| w for its basis indices n, and finally the moments are
276+ * turned into node weights w_j = sum_n U[n, j] I_n.
277+ *
278+ * Self items are 2 vec4: [1, edge, xi0, xi1], [h, xc, heff, 0]: on edge e the
279+ * polar rule uses rho in [0, 1] along the ray and the sinh substitution
280+ * x = xc + heff sinh(xi), xi in [xi0, xi1], along the edge.
281+ * tables: umat (npols*npols) followed by the GL rules gl_x, gl_w (nq), glr_x,
282+ * glr_w (nqr, on [0, 1]) and glx_x, glx_w (nqx).
283+ */
284+export const NEAR_WEIGHTS = /* wgsl */ `
285+struct Params {
286+ n: u32,
287+ pair_offset: u32,
288+ norder: i32,
289+ npols: i32,
290+ nq: i32,
291+ nqr: i32,
292+ nqx: i32,
293+ zk: f32,
294+ ga: f32,
295+ gpe: i32,
296+}
297+@group(0) @binding(0) var<uniform> U: Params;
298+@group(0) @binding(1) var<storage, read> pairs_i: array<vec2i>;
299+@group(0) @binding(2) var<storage, read> pairs_f: array<vec4f>;
300+@group(0) @binding(3) var<storage, read> targ: array<vec4f>;
301+@group(0) @binding(4) var<storage, read> gpar: array<vec4f>;
302+@group(0) @binding(5) var<storage, read> item_ptr: array<u32>;
303+@group(0) @binding(6) var<storage, read> items: array<vec4f>;
304+@group(0) @binding(7) var<storage, read> tables: array<f32>;
305+@group(0) @binding(8) var<storage, read_write> W: array<vec2f>;
306+${COMMON}
307+${CHART}
308+
309+const NT = ${NT};
310+var<workgroup> sK: array<vec2f, ${6 * NT}>;
311+var<workgroup> sLu: array<f32, ${NORD1_MAX * NT}>;
312+var<workgroup> sLv: array<f32, ${NORD1_MAX * NT}>;
313+var<workgroup> sM: array<vec2f, ${6 * NPOLS_MAX}>;
314+var<workgroup> sNtot: i32;
315+
316+@compute @workgroup_size(${NT})
317+fn near_weights(@builtin(workgroup_id) wid: vec3u, @builtin(local_invocation_id) lid: vec3u) {
318+ let pl = wid.x;
319+ let p = pl + U.pair_offset;
320+ let t = i32(lid.x);
321+ let it = pairs_i[p].x;
322+ let ip = pairs_i[p].y;
323+ let ut = pairs_f[p].x;
324+ let vt = pairs_f[p].y;
325+ let kind = i32(pairs_f[p].z);
326+ let is_self = kind == 2;
327+ let onsurf = kind >= 1;
328+ let xt = targ[2 * it].xyz;
329+ let nt = targ[2 * it + 1].xyz;
330+ let n1 = U.norder + 1;
331+ let npols = U.npols;
332+ let i0 = item_ptr[pl];
333+ let nodes_per_item = select(U.nq * U.nq, U.nqr * U.nqx, is_self);
334+ if (t == 0) { sNtot = i32(item_ptr[pl + 1u] - i0) * nodes_per_item; }
335+ let ntot = workgroupUniformLoad(&sNtot);
336+ let oq = npols * npols; // gl tables start after umat
337+ let or = oq + 2 * U.nq; // glr_x
338+ let ox = or + 2 * U.nqr; // glx_x
339+
340+ var Ir: array<vec2f, 12>;
341+ for (var a = 0; a < 12; a++) { Ir[a] = vec2f(0.0, 0.0); }
342+
343+ for (var base = 0; base < ntot; base += NT) {
344+ let g = base + t;
345+ var K: array<vec2f, 6>;
346+ var Lu: array<f32, ${NORD1_MAX}>;
347+ var Lv: array<f32, ${NORD1_MAX}>;
348+ let live = g < ntot;
349+ if (live) {
350+ let item = g / nodes_per_item;
351+ let loc = g - item * nodes_per_item;
352+ let it0 = items[2u * (i0 + u32(item))];
353+ let it1 = items[2u * (i0 + u32(item)) + 1u];
354+ var u: f32; var v: f32; var wq: f32;
355+ if (is_self) {
356+ let e = i32(it0.y);
357+ let ir = loc / U.nqx;
358+ let ix = loc - ir * U.nqx;
359+ let h = it1.x; let xc = it1.y; let heff = it1.z;
360+ let rho = tables[or + ir];
361+ let wr = tables[or + U.nqr + ir];
362+ let xi = it0.z + 0.5 * (it0.w - it0.z) * (tables[ox + ix] + 1.0);
363+ let wxi = 0.5 * (it0.w - it0.z) * tables[ox + U.nqx + ix];
364+ let x = xc + heff * sinh(xi);
365+ let dx = heff * cosh(xi) * wxi;
366+ if (e < 2) { u = ut + rho * h; v = vt + rho * x; } else { u = ut + rho * x; v = vt + rho * h; }
367+ wq = abs(h) * rho * wr * dx;
368+ } else {
369+ let a = loc / U.nq;
370+ let b = loc - a * U.nq;
371+ u = it0.y + it0.w * tables[oq + a];
372+ v = it0.z + it0.w * tables[oq + b];
373+ wq = it0.w * it0.w * tables[oq + U.nq + a] * tables[oq + U.nq + b];
374+ }
375+ let G = chart_rel(ip, u, v, ut, vt, onsurf, xt, U.ga, U.gpe);
376+ var ns = cross(G.ru, G.rv);
377+ let jac = length(ns);
378+ ns = ns / jac;
379+ K = kernels6(-G.d, nt, ns, wq * jac, U.zk, true);
380+ legendre(u, U.norder, &Lu);
381+ legendre(v, U.norder, &Lv);
382+ }
383+ for (var m = 0; m < 6; m++) { sK[m * NT + t] = select(vec2f(0.0, 0.0), K[m], live); }
384+ for (var a = 0; a < n1; a++) {
385+ sLu[a * NT + t] = select(0.0, Lu[a], live);
386+ sLv[a * NT + t] = select(0.0, Lv[a], live);
387+ }
388+ workgroupBarrier();
389+ for (var acc = 0; acc < 2; acc++) {
390+ let n = t + acc * NT;
391+ if (n < npols) {
392+ let a = n % n1;
393+ let b = n / n1;
394+ for (var i = 0; i < NT; i++) {
395+ let P = sLu[a * NT + i] * sLv[b * NT + i];
396+ for (var m = 0; m < 6; m++) { Ir[acc * 6 + m] += P * sK[m * NT + i]; }
397+ }
398+ }
399+ }
400+ workgroupBarrier();
401+ }
402+
403+ // moments -> node weights, w_j = sum_n U[n, j] I_n
404+ for (var acc = 0; acc < 2; acc++) {
405+ let n = t + acc * NT;
406+ if (n < npols) { for (var m = 0; m < 6; m++) { sM[m * ${NPOLS_MAX} + n] = Ir[acc * 6 + m]; } }
407+ }
408+ workgroupBarrier();
409+ for (var j = t; j < npols; j += NT) {
410+ for (var m = 0; m < 6; m++) {
411+ var w = vec2f(0.0, 0.0);
412+ for (var n = 0; n < npols; n++) { w += tables[n * npols + j] * sM[m * ${NPOLS_MAX} + n]; }
413+ W[(u32(p) * 6u + u32(m)) * u32(npols) + u32(j)] = w;
414+ }
415+ }
416+}
417+`;
418+
419+/**
420+ * Far-field pass over the patches j0 <= j < j1: out[m][i] += sum over patches j
421+ * not near target i of sum_l K_m(x_i, y_l) w_l q_l over the oversampled nodes
422+ * of patch j, for the kernels selected in selmask (4 bits per kernel: 0 skip,
423+ * 1 density qa, 2 density qb). A patch is near target i when it appears in
424+ * the target's near list (row_ptr, pair_p, ascending), so the split between
425+ * this sum and the near weights is exact.
426+ * srcgeo: 2 vec4 per node, [x, w], [n, 0]; srcq: vec4 per node (qa, qb).
427+ */
428+export const DENSE_APPLY = /* wgsl */ `
429+struct Params {
430+ ntarg: u32,
431+ j0: u32,
432+ j1: u32,
433+ npols_o: u32,
434+ zk: f32,
435+ selmask: u32,
436+}
437+@group(0) @binding(0) var<uniform> U: Params;
438+@group(0) @binding(1) var<storage, read> targ: array<vec4f>;
439+@group(0) @binding(2) var<storage, read> srcgeo: array<vec4f>;
440+@group(0) @binding(3) var<storage, read> srcq: array<vec4f>;
441+@group(0) @binding(4) var<storage, read> row_ptr: array<u32>;
442+@group(0) @binding(5) var<storage, read> pair_p: array<u32>;
443+@group(0) @binding(6) var<storage, read_write> out: array<vec2f>;
444+${COMMON}
445+
446+const WG = ${DENSE_WG}u;
447+var<workgroup> tile: array<vec4f, ${3 * DENSE_WG}>;
448+
449+@compute @workgroup_size(${DENSE_WG})
450+fn dense_apply(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id) lid: vec3u) {
451+ let i = gid.x;
452+ let valid = i < U.ntarg;
453+ var xt = vec3f(0.0);
454+ var nt = vec3f(0.0);
455+ var cur = 0u;
456+ var end = 0u;
457+ if (valid) {
458+ xt = targ[2u * i].xyz;
459+ nt = targ[2u * i + 1u].xyz;
460+ cur = row_ptr[i];
461+ end = row_ptr[i + 1u];
462+ while (cur < end && pair_p[cur] < U.j0) { cur++; }
463+ }
464+ var sel: array<u32, 6>;
465+ for (var m = 0u; m < 6u; m++) { sel[m] = (U.selmask >> (4u * m)) & 3u; }
466+ let need3 = sel[3] != 0u;
467+ var acc: array<vec2f, 6>;
468+ for (var m = 0; m < 6; m++) { acc[m] = vec2f(0.0, 0.0); }
469+ for (var j = U.j0; j < U.j1; j++) {
470+ var far = valid;
471+ if (valid && cur < end && pair_p[cur] == j) { far = false; cur++; }
472+ for (var base = 0u; base < U.npols_o; base += WG) {
473+ let l = base + lid.x;
474+ if (l < U.npols_o) {
475+ let s = j * U.npols_o + l;
476+ tile[3u * lid.x] = srcgeo[2u * s];
477+ tile[3u * lid.x + 1u] = srcgeo[2u * s + 1u];
478+ tile[3u * lid.x + 2u] = srcq[s];
479+ }
480+ workgroupBarrier();
481+ if (far) {
482+ let nl = min(WG, U.npols_o - base);
483+ var part: array<vec2f, 6>;
484+ for (var m = 0; m < 6; m++) { part[m] = vec2f(0.0, 0.0); }
485+ for (var s = 0u; s < nl; s++) {
486+ let g0 = tile[3u * s];
487+ let g1 = tile[3u * s + 1u];
488+ let q = tile[3u * s + 2u];
489+ let K = kernels6(xt - g0.xyz, nt, g1.xyz, g0.w, U.zk, need3);
490+ for (var m = 0u; m < 6u; m++) {
491+ if (sel[m] == 1u) { part[m] += cmul(K[m], q.xy); }
492+ else if (sel[m] == 2u) { part[m] += cmul(K[m], q.zw); }
493+ }
494+ }
495+ for (var m = 0; m < 6; m++) { acc[m] += part[m]; }
496+ }
497+ workgroupBarrier();
498+ }
499+ }
500+ if (valid) {
501+ for (var m = 0u; m < 6u; m++) {
502+ if (sel[m] != 0u) { out[m * U.ntarg + i] += acc[m]; }
503+ }
504+ }
505+}
506+`;
507+
508+/**
509+ * Near contributions: out[m][i] += sum over the near pairs p of target i of
510+ * sum_l W[p][m][l] q[d][patch nodes l], for up to four requests (m, d) packed
511+ * as m | d << 8 in req.
512+ */
513+export const SPARSE_APPLY = /* wgsl */ `
514+struct Params {
515+ ntarg: u32,
516+ npols: u32,
517+ npts: u32,
518+ nreq: u32,
519+ req: vec4u,
520+}
521+@group(0) @binding(0) var<uniform> U: Params;
522+@group(0) @binding(1) var<storage, read> row_ptr: array<u32>;
523+@group(0) @binding(2) var<storage, read> pair_p: array<u32>;
524+@group(0) @binding(3) var<storage, read> W: array<vec2f>;
525+@group(0) @binding(4) var<storage, read> q: array<vec2f>;
526+@group(0) @binding(5) var<storage, read_write> out: array<vec2f>;
527+
528+fn cmul(a: vec2f, b: vec2f) -> vec2f {
529+ return vec2f(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
530+}
531+
532+@compute @workgroup_size(128)
533+fn sparse_apply(@builtin(global_invocation_id) gid: vec3u) {
534+ let i = gid.x;
535+ if (i >= U.ntarg) { return; }
536+ var acc: array<vec2f, 4>;
537+ for (var r = 0u; r < 4u; r++) { acc[r] = vec2f(0.0, 0.0); }
538+ for (var p = row_ptr[i]; p < row_ptr[i + 1u]; p++) {
539+ let j = pair_p[p];
540+ for (var r = 0u; r < U.nreq; r++) {
541+ let m = U.req[r] & 0xffu;
542+ let d = U.req[r] >> 8u;
543+ let wbase = (p * 6u + m) * U.npols;
544+ let qbase = d * U.npts + j * U.npols;
545+ var s = vec2f(0.0, 0.0);
546+ for (var l = 0u; l < U.npols; l++) { s += cmul(W[wbase + l], q[qbase + l]); }
547+ acc[r] += s;
548+ }
549+ }
550+ for (var r = 0u; r < U.nreq; r++) {
551+ let m = U.req[r] & 0xffu;
552+ out[m * U.ntarg + i] += acc[r];
553+ }
554+}
555+`;
src/main.tsadded+487−0View file
@@ -0,0 +1,487 @@
1+/**
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+ */
8+import { roundedCube, type Surface } from './surface.ts';
9+import {
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';
22+import { adapterDescription, Gpu } from './gpu/device.ts';
23+import { Operator } from './gpu/operator.ts';
24+import { colorLimits, DB_MAX, DB_MIN, type Field, type Mode, type SliceField } from './view/field.ts';
25+import { SliceView } from './view/slices.ts';
26+import { RoomView } from './view/room3d.ts';
27+import { drawColorbar, DIVERGING, SEQUENTIAL } from './view/colormaps.ts';
28+import { drawConvergence } from './view/convergence.ts';
29+
30+const $ = <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+};
35+
36+const A = 1.0;
37+const EXTERIOR_SOURCE = [1.7, 0.6, -0.5];
38+
39+interface 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+}
50+
51+const 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+};
57+const 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 };
58+
59+const 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+};
89+
90+const fmtBytes = (b: number) => (b >= 2 ** 30 ? `${(b / 2 ** 30).toFixed(2)} GB` : `${(b / 2 ** 20).toFixed(0)} MB`);
91+const fmtS = (s: number) => (s < 10 ? `${s.toFixed(2)} s` : `${s.toFixed(1)} s`);
92+
93+// ------------------------------------------------------------------ params
94+
95+function 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+}
108+
109+function 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+}
122+
123+function 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+}
128+
129+function 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+}
136+
137+// ------------------------------------------------------------------ derived
138+
139+let deviceLimit = Infinity;
140+
141+function 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+}
168+
169+function updateDerived() {
170+ const d = describe(readParams());
171+ ui.derived.innerHTML = d.html;
172+ if (!running) ui.solve.disabled = d.blocking;
173+}
174+
175+// ------------------------------------------------------------------ status
176+
177+let running = false;
178+let cancelled = false;
179+class Cancelled extends Error {}
180+
181+function 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+}
186+
187+function checkCancel() {
188+ if (cancelled) throw new Cancelled('cancelled');
189+}
190+
191+// ------------------------------------------------------------------ results
192+
193+type Row = [string, string, string?];
194+
195+function 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+}
209+
210+const quality = (e: number, tol: number) => (e < 10 * tol ? '' : e < 1e-2 ? 'warn' : 'bad');
211+
212+// ------------------------------------------------------------------ main
213+
214+async 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);
231+
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'));
244+
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;
251+
252+ const redrawConv = () => drawConvergence(ui.conv, conv, lastEps);
253+
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+ };
267+
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 = [];
328+
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+ }
346+
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)`]);
357+
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)]);
368+
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+ }
421+
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()));
444+
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);
468+
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+}
480+
481+start().catch((e) => {
482+ console.error(e);
483+ ui.err.textContent = e instanceof Error ? e.message : String(e);
484+ setStatus('unavailable', 0);
485+});
486+
487+export type { CVec };
src/problem.tsadded+168−0View file
@@ -0,0 +1,168 @@
1+/**
2+ * The showcase problem: a point source at x0 inside a rounded cube (the
3+ * superellipsoid |x|^p + |y|^p + |z|^p = a^p) with impedance walls. With time
4+ * dependence exp(-i omega t) the pressure amplitude u satisfies
5+ *
6+ * (Delta + k^2) u = -delta(x - x0) in the room,
7+ * du/dn + i k lambda u = 0 on the walls,
8+ *
9+ * with n the outward normal and lambda = -beta, where beta > 0 is the wall
10+ * admittance. We write u = G + v with G the free-space field of the source and
11+ * find v from fmm3dbie's right-preconditioned combined field equation.
12+ *
13+ * Complex vectors are interleaved Float64Arrays: [re0, im0, re1, im1, ...].
14+ */
15+
16+export type CVec = Float64Array;
17+
18+export const ALPHA = 1.0;
19+
20+/** Free-space Green's function exp(ikr)/(4 pi r) at points x (3, n) for a source x0. */
21+export function green(x: Float64Array, x0: ArrayLike<number>, k: number): CVec {
22+ const n = x.length / 3;
23+ const out = new Float64Array(2 * n);
24+ for (let i = 0; i < n; i++) {
25+ const r = Math.hypot(x[i] - x0[0], x[n + i] - x0[1], x[2 * n + i] - x0[2]);
26+ const s = 1 / (4 * Math.PI * r);
27+ out[2 * i] = s * Math.cos(k * r);
28+ out[2 * i + 1] = s * Math.sin(k * r);
29+ }
30+ return out;
31+}
32+
33+/** Normal derivative of the Green's function at x with unit normals nrm (3, n). */
34+export function greenDn(x: Float64Array, nrm: Float64Array, x0: ArrayLike<number>, k: number): CVec {
35+ const n = x.length / 3;
36+ const out = new Float64Array(2 * n);
37+ for (let i = 0; i < n; i++) {
38+ const dx = x[i] - x0[0];
39+ const dy = x[n + i] - x0[1];
40+ const dz = x[2 * n + i] - x0[2];
41+ const r = Math.hypot(dx, dy, dz);
42+ const s = 1 / (4 * Math.PI * r);
43+ const gr = s * Math.cos(k * r);
44+ const gi = s * Math.sin(k * r);
45+ const dn = (dx * nrm[i] + dy * nrm[n + i] + dz * nrm[2 * n + i]) / r;
46+ // g * (ik - 1/r) * dn
47+ const fr = -1 / r;
48+ const fi = k;
49+ out[2 * i] = (gr * fr - gi * fi) * dn;
50+ out[2 * i + 1] = (gr * fi + gi * fr) * dn;
51+ }
52+ return out;
53+}
54+
55+export function inside(x: ArrayLike<number>, a: number, p: number, margin = 0): boolean {
56+ return Math.pow(Math.abs(x[0] / a), p) + Math.pow(Math.abs(x[1] / a), p) + Math.pow(Math.abs(x[2] / a), p) < (1 - margin) ** p;
57+}
58+
59+/** Impedance data for v = u - G: -(dG/dn + i k lambda G) at the wall nodes. */
60+export function roomRhs(xs: Float64Array, nrm: Float64Array, x0: ArrayLike<number>, k: number, lam: number): CVec {
61+ const g = green(xs, x0, k);
62+ const gn = greenDn(xs, nrm, x0, k);
63+ const n = g.length / 2;
64+ const out = new Float64Array(2 * n);
65+ for (let i = 0; i < n; i++) {
66+ // i k lam g = (-k lam g_im, k lam g_re)
67+ out[2 * i] = -(gn[2 * i] - k * lam * g[2 * i + 1]);
68+ out[2 * i + 1] = -(gn[2 * i + 1] + k * lam * g[2 * i]);
69+ }
70+ return out;
71+}
72+
73+/** Impedance data of the field of a source outside the room (the self-test). */
74+export function selftestRhs(xs: Float64Array, nrm: Float64Array, xe: ArrayLike<number>, k: number, lam: number): CVec {
75+ const r = roomRhs(xs, nrm, xe, k, lam);
76+ for (let i = 0; i < r.length; i++) r[i] = -r[i];
77+ return r;
78+}
79+
80+/** n deterministic pseudo-random points in the cube of half-width 0.8 a. */
81+export function selftestPoints(n: number, a: number): Float64Array {
82+ const out = new Float64Array(3 * n);
83+ let s = 12345;
84+ const rnd = () => {
85+ s = (Math.imul(s, 1664525) + 1013904223) >>> 0;
86+ return s / 4294967296;
87+ };
88+ for (let i = 0; i < 3 * n; i++) out[i] = 0.8 * a * (2 * rnd() - 1);
89+ return out;
90+}
91+
92+export interface Slice {
93+ name: 'xy' | 'xz' | 'yz';
94+ axes: [number, number];
95+ fixed: number;
96+ pts: Float64Array; // (3, ngrid^2), row-major over the grid with the first axis varying fastest
97+ mask: Uint8Array;
98+}
99+
100+/** Points on the three planes through x0; each grid runs from -a to a. */
101+export function slicePoints(ngrid: number, x0: ArrayLike<number>, a: number, p: number): { grid: Float64Array; slices: Slice[] } {
102+ const grid = new Float64Array(ngrid);
103+ for (let i = 0; i < ngrid; i++) grid[i] = -a + (2 * a * i) / (ngrid - 1);
104+ const defs: [Slice['name'], [number, number], number][] = [
105+ ['xy', [0, 1], 2],
106+ ['xz', [0, 2], 1],
107+ ['yz', [1, 2], 0],
108+ ];
109+ const slices: Slice[] = [];
110+ for (const [name, axes, fixed] of defs) {
111+ const N = ngrid * ngrid;
112+ const pts = new Float64Array(3 * N);
113+ const mask = new Uint8Array(N);
114+ for (let row = 0; row < ngrid; row++) {
115+ for (let col = 0; col < ngrid; col++) {
116+ const q = row * ngrid + col;
117+ const x = [0, 0, 0];
118+ x[axes[0]] = grid[col];
119+ x[axes[1]] = grid[row];
120+ x[fixed] = x0[fixed];
121+ pts[q] = x[0];
122+ pts[N + q] = x[1];
123+ pts[2 * N + q] = x[2];
124+ mask[q] = inside(x, a, p, 1e-3) ? 1 : 0;
125+ }
126+ }
127+ slices.push({ name, axes, fixed, pts, mask });
128+ }
129+ return { grid, slices };
130+}
131+
132+/** Power radiated by the source, k/(4 pi) + Im v(x0), and absorbed by the walls, -k Re(lambda) int |u|^2 dS. */
133+export function energyBalance(k: number, lam: number, vAtSource: [number, number], wts: Float64Array, uWall: CVec) {
134+ const pSrc = k / (4 * Math.PI) + vAtSource[1];
135+ let s = 0;
136+ for (let i = 0; i < wts.length; i++) s += wts[i] * (uWall[2 * i] ** 2 + uWall[2 * i + 1] ** 2);
137+ const pWall = -k * lam * s;
138+ return { pSrc, pWall, mismatch: Math.abs(pSrc - pWall) / Math.abs(pSrc) };
139+}
140+
141+/** Relative 2-norm difference of two complex vectors. */
142+export function relError(u: CVec, uex: CVec): number {
143+ let num = 0;
144+ let den = 0;
145+ for (let i = 0; i < u.length; i++) {
146+ num += (u[i] - uex[i]) ** 2;
147+ den += uex[i] ** 2;
148+ }
149+ return Math.sqrt(num / den);
150+}
151+
152+/** Select the columns of x (3, n) where mask is set. */
153+export function maskedPoints(x: Float64Array, mask: Uint8Array): { pts: Float64Array; index: Int32Array } {
154+ const n = mask.length;
155+ let cnt = 0;
156+ for (let i = 0; i < n; i++) cnt += mask[i];
157+ const pts = new Float64Array(3 * cnt);
158+ const index = new Int32Array(cnt);
159+ let q = 0;
160+ for (let i = 0; i < n; i++) {
161+ if (!mask[i]) continue;
162+ pts[q] = x[i];
163+ pts[cnt + q] = x[n + i];
164+ pts[2 * cnt + q] = x[2 * n + i];
165+ index[q++] = i;
166+ }
167+ return { pts, index };
168+}
src/reference.tsadded+186−0View file
@@ -0,0 +1,186 @@
1+/**
2+ * Double precision CPU reference for the tests: the six kernels, and brute
3+ * force near weights for one (target, patch) pair (a port of
4+ * tests/test_quadrature.py from the CUDA proof of concept).
5+ *
6+ * The six kernels, in the order used by fmm3dbie's helm_rpcomb_imp solver:
7+ * 0: S'_k 1: S_kappa 2: S'_kappa 3: D'_k - D'_kappa 4: S_k 5: D_k
8+ * with kappa = i|k|, S' the normal derivative at the target and D the normal
9+ * derivative at the source.
10+ */
11+import { leggauss } from './gauss.ts';
12+import { basis2d, type Surface } from './surface.ts';
13+
14+const OVER4PI = 1 / (4 * Math.PI);
15+
16+/** K (6, n) complex, interleaved per kernel: out[m][2q], out[m][2q+1]. */
17+export function kernels6(xt: ArrayLike<number>, nt: ArrayLike<number>, y: Float64Array, ns: Float64Array, k: number): Float64Array[] {
18+ const n = y.length / 3;
19+ const K = Array.from({ length: 6 }, () => new Float64Array(2 * n));
20+ for (let q = 0; q < n; q++) {
21+ const dx = xt[0] - y[q];
22+ const dy = xt[1] - y[n + q];
23+ const dz = xt[2] - y[2 * n + q];
24+ const r2 = dx * dx + dy * dy + dz * dz;
25+ const r = Math.sqrt(r2);
26+ const dnt = dx * nt[0] + dy * nt[1] + dz * nt[2];
27+ const dns = dx * ns[q] + dy * ns[n + q] + dz * ns[2 * n + q];
28+ const nn = ns[q] * nt[0] + ns[n + q] * nt[1] + ns[2 * n + q] * nt[2];
29+ const e0 = OVER4PI / r;
30+ const Ekr = e0 * Math.cos(k * r);
31+ const Eki = e0 * Math.sin(k * r);
32+ const Ec = e0 * Math.exp(-k * r);
33+ // fk = 1 - i k r, fc = 1 + k r
34+ const fkr = 1;
35+ const fki = -k * r;
36+ const fc = 1 + k * r;
37+ // gk = 3 - 3 i k r - k^2 r^2, gc = 3 + 3 k r + k^2 r^2
38+ const gkr = 3 - k * k * r2;
39+ const gki = -3 * k * r;
40+ const gc = 3 + 3 * k * r + k * k * r2;
41+ const fkEkr = fkr * Ekr - fki * Eki;
42+ const fkEki = fkr * Eki + fki * Ekr;
43+ const gkEkr = gkr * Ekr - gki * Eki;
44+ const gkEki = gkr * Eki + gki * Ekr;
45+ K[4][2 * q] = Ekr;
46+ K[4][2 * q + 1] = Eki;
47+ K[1][2 * q] = Ec;
48+ K[0][2 * q] = (-dnt * fkEkr) / r2;
49+ K[0][2 * q + 1] = (-dnt * fkEki) / r2;
50+ K[2][2 * q] = (-dnt * fc * Ec) / r2;
51+ K[5][2 * q] = (dns * fkEkr) / r2;
52+ K[5][2 * q + 1] = (dns * fkEki) / r2;
53+ const c1 = nn / r2;
54+ const c2 = (dns * dnt) / (r2 * r2);
55+ K[3][2 * q] = c1 * fkEkr - c2 * gkEkr - (c1 * fc * Ec - c2 * gc * Ec);
56+ K[3][2 * q + 1] = c1 * fkEki - c2 * gkEki;
57+ }
58+ return K;
59+}
60+
61+/** Moments sum_nodes K_m P_n |J| w -> (6, npols) complex. */
62+export function momentsAt(surf: Surface, patch: number, uv: Float64Array, w: Float64Array, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number): Float64Array[] {
63+ const n = w.length;
64+ const g = surf.evalGeometry(patch, uv);
65+ const y = g.subarray(0, 3 * n);
66+ const ns = new Float64Array(3 * n);
67+ const jac = new Float64Array(n);
68+ for (let q = 0; q < n; q++) {
69+ const ru = [g[3 * n + q], g[4 * n + q], g[5 * n + q]];
70+ const rv = [g[6 * n + q], g[7 * n + q], g[8 * n + q]];
71+ const cx = ru[1] * rv[2] - ru[2] * rv[1];
72+ const cy = ru[2] * rv[0] - ru[0] * rv[2];
73+ const cz = ru[0] * rv[1] - ru[1] * rv[0];
74+ jac[q] = Math.hypot(cx, cy, cz);
75+ ns[q] = cx / jac[q];
76+ ns[n + q] = cy / jac[q];
77+ ns[2 * n + q] = cz / jac[q];
78+ }
79+ const K = kernels6(xt, nt, y, ns, k);
80+ const P = basis2d(uv, surf.norder);
81+ const npols = surf.npols;
82+ const I = Array.from({ length: 6 }, () => new Float64Array(2 * npols));
83+ for (let m = 0; m < 6; m++) {
84+ for (let nn = 0; nn < npols; nn++) {
85+ let sr = 0;
86+ let si = 0;
87+ for (let q = 0; q < n; q++) {
88+ const f = P[nn * n + q] * jac[q] * w[q];
89+ sr += K[m][2 * q] * f;
90+ si += K[m][2 * q + 1] * f;
91+ }
92+ I[m][2 * nn] = sr;
93+ I[m][2 * nn + 1] = si;
94+ }
95+ }
96+ return I;
97+}
98+
99+function addInto(I: Float64Array[], J: Float64Array[]) {
100+ for (let m = 0; m < 6; m++) for (let i = 0; i < I[m].length; i++) I[m][i] += J[m][i];
101+}
102+
103+/** Self-patch moments: polar rules on the triangles (target, edge), split at the foot of the perpendicular, many nodes. */
104+export function refSelf(surf: Surface, patch: number, u0: number, v0: number, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number, nx = 300, nr = 60): Float64Array[] {
105+ const gr = leggauss(nr);
106+ const gx = leggauss(nx);
107+ const I = Array.from({ length: 6 }, () => new Float64Array(2 * surf.npols));
108+ for (let e = 0; e < 4; e++) {
109+ const h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e];
110+ if (Math.abs(h) < 1e-14) continue;
111+ const [x1, x2] = e < 2 ? [-1 - v0, 1 - v0] : [-1 - u0, 1 - u0];
112+ for (const xend of [x1, x2]) {
113+ if (Math.abs(xend) < 1e-14) continue;
114+ const xi1 = Math.asinh(Math.abs(xend) / Math.abs(h));
115+ const n = nr * nx;
116+ const uv = new Float64Array(2 * n);
117+ const w = new Float64Array(n);
118+ for (let a = 0; a < nr; a++) {
119+ const rho = 0.5 * (gr.x[a] + 1);
120+ const wr = 0.5 * gr.w[a];
121+ for (let b = 0; b < nx; b++) {
122+ const xi = 0.5 * xi1 * (gx.x[b] + 1);
123+ const x = Math.sign(xend) * Math.abs(h) * Math.sinh(xi);
124+ const dx = Math.abs(h) * Math.cosh(xi) * 0.5 * xi1 * gx.w[b];
125+ const q = a * nx + b;
126+ if (e < 2) {
127+ uv[q] = u0 + rho * h;
128+ uv[n + q] = v0 + rho * x;
129+ } else {
130+ uv[q] = u0 + rho * x;
131+ uv[n + q] = v0 + rho * h;
132+ }
133+ w[q] = wr * dx * Math.abs(h) * rho;
134+ }
135+ }
136+ addInto(I, momentsAt(surf, patch, uv, w, xt, nt, k));
137+ }
138+ }
139+ return I;
140+}
141+
142+/** Near-patch moments with a uniform quadtree of the given depth and an nq x nq rule per leaf. */
143+export function refLeaf(surf: Surface, patch: number, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number, depth = 5, nq = 8): Float64Array[] {
144+ const g = leggauss(nq);
145+ const nl = 2 ** depth;
146+ const h = 1 / nl;
147+ const n = nl * nl * nq * nq;
148+ const uv = new Float64Array(2 * n);
149+ const w = new Float64Array(n);
150+ let q = 0;
151+ for (let i = 0; i < nl; i++) {
152+ for (let j = 0; j < nl; j++) {
153+ const cu = -1 + h + 2 * h * i;
154+ const cv = -1 + h + 2 * h * j;
155+ for (let a = 0; a < nq; a++) {
156+ for (let b = 0; b < nq; b++) {
157+ uv[q] = cu + h * g.x[a];
158+ uv[n + q] = cv + h * g.x[b];
159+ w[q] = h * h * g.w[a] * g.w[b];
160+ q++;
161+ }
162+ }
163+ }
164+ }
165+ return momentsAt(surf, patch, uv, w, xt, nt, k);
166+}
167+
168+/** Node weights w_j = sum_n U[n, j] I_n from the moments: (6, npols) complex. */
169+export function momentsToWeights(surf: Surface, I: Float64Array[]): Float64Array[] {
170+ const npols = surf.npols;
171+ const W = Array.from({ length: 6 }, () => new Float64Array(2 * npols));
172+ for (let m = 0; m < 6; m++) {
173+ for (let j = 0; j < npols; j++) {
174+ let sr = 0;
175+ let si = 0;
176+ for (let n = 0; n < npols; n++) {
177+ const u = surf.umat[n * npols + j];
178+ sr += u * I[m][2 * n];
179+ si += u * I[m][2 * n + 1];
180+ }
181+ W[m][2 * j] = sr;
182+ W[m][2 * j + 1] = si;
183+ }
184+ }
185+ return W;
186+}
src/surface.tsadded+325−0View file
@@ -0,0 +1,325 @@
1+/**
2+ * Quadrilateral patch surfaces in the fmm3dbie convention (iptype 11), and the
3+ * rounded cube (superellipsoid) built from them. A direct port of
4+ * fmm3dbie_gpu/surface.py from the CUDA proof of concept.
5+ *
6+ * Each patch is the image of the reference square [-1, 1]^2 and carries
7+ * (norder+1)^2 tensor Gauss-Legendre nodes, with the u index varying fastest,
8+ * and an unnormalized Legendre basis P_j(u) P_i(v) in the same order. Arrays
9+ * with a shape (r, c) are stored row-major in one Float64Array, index r*c + i.
10+ */
11+import { leggauss, legendreVander } from './gauss.ts';
12+
13+/** Reference nodes (2, npols) with u in the first row, and the tensor weights. */
14+export function discNodes(norder: number): { uv: Float64Array; w: Float64Array } {
15+ const n1 = norder + 1;
16+ const { x, w } = leggauss(n1);
17+ const npols = n1 * n1;
18+ const uv = new Float64Array(2 * npols);
19+ const ww = new Float64Array(npols);
20+ for (let i = 0; i < n1; i++) {
21+ for (let j = 0; j < n1; j++) {
22+ const n = i * n1 + j;
23+ uv[n] = x[j];
24+ uv[npols + n] = x[i];
25+ ww[n] = w[i] * w[j];
26+ }
27+ }
28+ return { uv, w: ww };
29+}
30+
31+/** Basis values (npols, npts): index n = i*(norder+1) + j is P_j(u) P_i(v). */
32+export function basis2d(uv: Float64Array, norder: number): Float64Array {
33+ const npts = uv.length / 2;
34+ const n1 = norder + 1;
35+ const Pu = legendreVander(uv.subarray(0, npts), norder);
36+ const Pv = legendreVander(uv.subarray(npts), norder);
37+ const B = new Float64Array(n1 * n1 * npts);
38+ for (let i = 0; i < n1; i++) {
39+ for (let j = 0; j < n1; j++) {
40+ const n = i * n1 + j;
41+ for (let q = 0; q < npts; q++) B[n * npts + q] = Pv[i * npts + q] * Pu[j * npts + q];
42+ }
43+ }
44+ return B;
45+}
46+
47+/** U (npols, npols) with coefs = U @ vals for values at the tensor GL nodes. */
48+export function valsToCoefsMatrix(norder: number): Float64Array {
49+ const n1 = norder + 1;
50+ const { x, w } = leggauss(n1);
51+ const P = legendreVander(x, norder);
52+ const U1 = new Float64Array(n1 * n1);
53+ for (let k = 0; k < n1; k++) for (let j = 0; j < n1; j++) U1[k * n1 + j] = ((2 * k + 1) / 2) * P[k * n1 + j] * w[j];
54+ const npols = n1 * n1;
55+ const U = new Float64Array(npols * npols);
56+ for (let i = 0; i < n1; i++) {
57+ for (let ii = 0; ii < n1; ii++) {
58+ for (let j = 0; j < n1; j++) {
59+ for (let jj = 0; jj < n1; jj++) U[(i * n1 + ii) * npols + (j * n1 + jj)] = U1[i * n1 + j] * U1[ii * n1 + jj];
60+ }
61+ }
62+ }
63+ return U;
64+}
65+
66+/** M (npols_o, npols) with values at the order-nover nodes = M @ values at the order-norder nodes. */
67+export function oversampleMatrix(norder: number, nover: number): Float64Array {
68+ const { uv } = discNodes(nover);
69+ const npolsO = (nover + 1) ** 2;
70+ const npols = (norder + 1) ** 2;
71+ const B = basis2d(uv, norder); // (npols, npols_o)
72+ const U = valsToCoefsMatrix(norder);
73+ const M = new Float64Array(npolsO * npols);
74+ for (let k = 0; k < npolsO; k++) {
75+ for (let m = 0; m < npols; m++) {
76+ let s = 0;
77+ for (let n = 0; n < npols; n++) s += B[n * npolsO + k] * U[n * npols + m];
78+ M[k * npols + m] = s;
79+ }
80+ }
81+ return M;
82+}
83+
84+/** (n, e1, e2) with e1 x e2 = n for the six faces of the cube. */
85+export const FACES: [number[], number[], number[]][] = [
86+ [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
87+ [[-1, 0, 0], [0, 0, 1], [0, 1, 0]],
88+ [[0, 1, 0], [0, 0, 1], [1, 0, 0]],
89+ [[0, -1, 0], [1, 0, 0], [0, 0, 1]],
90+ [[0, 0, 1], [1, 0, 0], [0, 1, 0]],
91+ [[0, 0, -1], [0, 1, 0], [1, 0, 0]],
92+];
93+
94+/** m+1 breakpoints on [-1, 1], blending uniform and sine-clustered spacing. */
95+export function gradedBreaks(m: number, grading: number): Float64Array {
96+ const br = new Float64Array(m + 1);
97+ for (let i = 0; i <= m; i++) {
98+ const t = -1 + (2 * i) / m;
99+ br[i] = (1 - grading) * t + grading * Math.sin(0.5 * Math.PI * t);
100+ }
101+ return br;
102+}
103+
104+/**
105+ * xyz, r_u, r_v (9, n) of the face chart par = (n, e1, e2, s0, s1, t0, t1) of the
106+ * superellipsoid |x|^p = a^p at reference points uv (2, n).
107+ */
108+export function superellipsoidChart(par: ArrayLike<number>, a: number, p: number, uv: Float64Array): Float64Array {
109+ const n = uv.length / 2;
110+ const out = new Float64Array(9 * n);
111+ const hs = 0.5 * (par[10] - par[9]);
112+ const ht = 0.5 * (par[12] - par[11]);
113+ const c = [0, 0, 0];
114+ const gN = [0, 0, 0];
115+ for (let q = 0; q < n; q++) {
116+ const s = par[9] + (uv[q] + 1) * hs;
117+ const t = par[11] + (uv[n + q] + 1) * ht;
118+ let S = 0;
119+ for (let d = 0; d < 3; d++) {
120+ c[d] = par[d] + s * par[3 + d] + t * par[6 + d];
121+ S += Math.pow(c[d], p);
122+ }
123+ const N = Math.pow(S, 1 / p);
124+ let dNu = 0;
125+ let dNv = 0;
126+ for (let d = 0; d < 3; d++) {
127+ gN[d] = Math.pow(c[d] / N, p - 1);
128+ dNu += gN[d] * par[3 + d];
129+ dNv += gN[d] * par[6 + d];
130+ }
131+ for (let d = 0; d < 3; d++) {
132+ out[d * n + q] = (a * c[d]) / N;
133+ out[(3 + d) * n + q] = a * hs * (par[3 + d] / N - (c[d] * dNu) / (N * N));
134+ out[(6 + d) * n + q] = a * ht * (par[6 + d] / N - (c[d] * dNv) / (N * N));
135+ }
136+ }
137+ return out;
138+}
139+
140+export interface ExactGeometry {
141+ params: Float64Array; // (npatches, 13)
142+ a: number;
143+ p: number;
144+}
145+
146+export class Surface {
147+ readonly npols: number;
148+ readonly npts: number;
149+ readonly npatches: number;
150+ readonly uvs: Float64Array; // (2, npols) reference nodes
151+ readonly wref: Float64Array;
152+ readonly umat: Float64Array; // (npols, npols)
153+ readonly srccoefs: Float64Array; // (9, npts)
154+ readonly jac: Float64Array;
155+ readonly wts: Float64Array;
156+ readonly centroids: Float64Array; // (3, npatches)
157+ readonly rads: Float64Array;
158+ readonly patchId: Int32Array;
159+ readonly uvNodes: Float64Array; // (2, npts)
160+ readonly norder: number;
161+ readonly srcvals: Float64Array; // (12, npts): xyz, r_u, r_v, unit normal
162+ readonly exact: ExactGeometry;
163+
164+ constructor(norder: number, srcvals: Float64Array, exact: ExactGeometry) {
165+ this.norder = norder;
166+ this.srcvals = srcvals;
167+ this.exact = exact;
168+ this.npols = (norder + 1) ** 2;
169+ this.npts = srcvals.length / 12;
170+ if (this.npts % this.npols !== 0) throw new Error('srcvals is not a whole number of patches');
171+ this.npatches = this.npts / this.npols;
172+ const { uv, w } = discNodes(norder);
173+ this.uvs = uv;
174+ this.wref = w;
175+ this.umat = valsToCoefsMatrix(norder);
176+ const { npols, npts, npatches } = this;
177+ this.srccoefs = new Float64Array(9 * npts);
178+ for (let c = 0; c < 9; c++) {
179+ for (let pch = 0; pch < npatches; pch++) {
180+ const base = pch * npols;
181+ for (let n = 0; n < npols; n++) {
182+ let s = 0;
183+ for (let j = 0; j < npols; j++) s += this.umat[n * npols + j] * srcvals[c * npts + base + j];
184+ this.srccoefs[c * npts + base + n] = s;
185+ }
186+ }
187+ }
188+ this.jac = new Float64Array(npts);
189+ this.wts = new Float64Array(npts);
190+ for (let i = 0; i < npts; i++) {
191+ const ru = [srcvals[3 * npts + i], srcvals[4 * npts + i], srcvals[5 * npts + i]];
192+ const rv = [srcvals[6 * npts + i], srcvals[7 * npts + i], srcvals[8 * npts + i]];
193+ const cx = ru[1] * rv[2] - ru[2] * rv[1];
194+ const cy = ru[2] * rv[0] - ru[0] * rv[2];
195+ const cz = ru[0] * rv[1] - ru[1] * rv[0];
196+ this.jac[i] = Math.hypot(cx, cy, cz);
197+ this.wts[i] = this.jac[i] * w[i % npols];
198+ }
199+ // centroid and radius of each patch from a 5x5 sample grid of the exact chart
200+ const g5 = new Float64Array(50);
201+ for (let i = 0; i < 5; i++) {
202+ for (let j = 0; j < 5; j++) {
203+ g5[i * 5 + j] = -1 + 0.5 * i;
204+ g5[25 + i * 5 + j] = -1 + 0.5 * j;
205+ }
206+ }
207+ this.centroids = new Float64Array(3 * npatches);
208+ this.rads = new Float64Array(npatches);
209+ for (let pch = 0; pch < npatches; pch++) {
210+ const g = this.evalGeometry(pch, g5);
211+ const cen = [0, 0, 0];
212+ for (let k = 0; k < 25; k++) for (let d = 0; d < 3; d++) cen[d] += g[d * 25 + k] / 25;
213+ let r2 = 0;
214+ for (let k = 0; k < 25; k++) {
215+ let s = 0;
216+ for (let d = 0; d < 3; d++) s += (g[d * 25 + k] - cen[d]) ** 2;
217+ r2 = Math.max(r2, s);
218+ }
219+ for (let d = 0; d < 3; d++) this.centroids[d * npatches + pch] = cen[d];
220+ this.rads[pch] = Math.sqrt(r2);
221+ }
222+ this.patchId = new Int32Array(npts);
223+ this.uvNodes = new Float64Array(2 * npts);
224+ for (let i = 0; i < npts; i++) {
225+ this.patchId[i] = Math.floor(i / npols);
226+ this.uvNodes[i] = uv[i % npols];
227+ this.uvNodes[npts + i] = uv[npols + (i % npols)];
228+ }
229+ }
230+
231+ /** xyz, r_u, r_v (9, n) of the exact chart of one patch at reference points uv (2, n). */
232+ evalGeometry(patch: number, uv: Float64Array): Float64Array {
233+ const e = this.exact;
234+ return superellipsoidChart(e.params.subarray(13 * patch, 13 * patch + 13), e.a, e.p, uv);
235+ }
236+
237+ /** Nodes of the order-nover rule on every patch: xyz (3, N), unit normal (3, N), weights (N). */
238+ oversampled(nover: number): { x: Float64Array; n: Float64Array; w: Float64Array } {
239+ const { uv, w: wo } = discNodes(nover);
240+ const npolsO = (nover + 1) ** 2;
241+ const N = npolsO * this.npatches;
242+ const x = new Float64Array(3 * N);
243+ const nrm = new Float64Array(3 * N);
244+ const w = new Float64Array(N);
245+ for (let pch = 0; pch < this.npatches; pch++) {
246+ const g = this.evalGeometry(pch, uv);
247+ for (let k = 0; k < npolsO; k++) {
248+ const i = pch * npolsO + k;
249+ const ru = [g[3 * npolsO + k], g[4 * npolsO + k], g[5 * npolsO + k]];
250+ const rv = [g[6 * npolsO + k], g[7 * npolsO + k], g[8 * npolsO + k]];
251+ const cx = ru[1] * rv[2] - ru[2] * rv[1];
252+ const cy = ru[2] * rv[0] - ru[0] * rv[2];
253+ const cz = ru[0] * rv[1] - ru[1] * rv[0];
254+ const jac = Math.hypot(cx, cy, cz);
255+ x[i] = g[k];
256+ x[N + i] = g[npolsO + k];
257+ x[2 * N + i] = g[2 * npolsO + k];
258+ nrm[i] = cx / jac;
259+ nrm[N + i] = cy / jac;
260+ nrm[2 * N + i] = cz / jac;
261+ w[i] = jac * wo[k];
262+ }
263+ }
264+ return { x, n: nrm, w };
265+ }
266+}
267+
268+/**
269+ * The superellipsoid |x|^p + |y|^p + |z|^p = a^p with m x m patches per cube
270+ * face, graded toward the face edges where the curvature is concentrated.
271+ */
272+export function roundedCube(a: number, p: number, m: number, norder: number, grading = 0.5): Surface {
273+ const { uv } = discNodes(norder);
274+ const npols = (norder + 1) ** 2;
275+ const br = gradedBreaks(m, grading);
276+ const npatches = 6 * m * m;
277+ const npts = npatches * npols;
278+ const params = new Float64Array(13 * npatches);
279+ const srcvals = new Float64Array(12 * npts);
280+ let pch = 0;
281+ for (const [n, e1, e2] of FACES) {
282+ for (let i = 0; i < m; i++) {
283+ for (let j = 0; j < m; j++) {
284+ const par = [...n, ...e1, ...e2, br[i], br[i + 1], br[j], br[j + 1]];
285+ params.set(par, 13 * pch);
286+ const g = superellipsoidChart(par, a, p, uv);
287+ for (let k = 0; k < npols; k++) {
288+ const idx = pch * npols + k;
289+ for (let c = 0; c < 9; c++) srcvals[c * npts + idx] = g[c * npols + k];
290+ const ru = [g[3 * npols + k], g[4 * npols + k], g[5 * npols + k]];
291+ const rv = [g[6 * npols + k], g[7 * npols + k], g[8 * npols + k]];
292+ const cx = ru[1] * rv[2] - ru[2] * rv[1];
293+ const cy = ru[2] * rv[0] - ru[0] * rv[2];
294+ const cz = ru[0] * rv[1] - ru[1] * rv[0];
295+ const len = Math.hypot(cx, cy, cz);
296+ srcvals[9 * npts + idx] = cx / len;
297+ srcvals[10 * npts + idx] = cy / len;
298+ srcvals[11 * npts + idx] = cz / len;
299+ }
300+ pch++;
301+ }
302+ }
303+ }
304+ return new Surface(norder, srcvals, { params, a, p });
305+}
306+
307+/** Exact volume of |x|^p + |y|^p + |z|^p <= a^p. */
308+export function superellipsoidVolume(a: number, p: number): number {
309+ return (8 * a ** 3 * gamma(1 + 1 / p) ** 3) / gamma(1 + 3 / p);
310+}
311+
312+/** Lanczos approximation of the Gamma function (used only for the volume check). */
313+function gamma(z: number): number {
314+ const g = 7;
315+ const c = [
316+ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059,
317+ 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
318+ ];
319+ if (z < 0.5) return Math.PI / (Math.sin(Math.PI * z) * gamma(1 - z));
320+ z -= 1;
321+ let x = c[0];
322+ for (let i = 1; i < g + 2; i++) x += c[i] / (z + i);
323+ const t = z + g + 0.5;
324+ return Math.sqrt(2 * Math.PI) * Math.pow(t, z + 0.5) * Math.exp(-t) * x;
325+}
src/view/colormaps.tsadded+34−0View file
@@ -0,0 +1,34 @@
1+/** Colormaps as 256-entry RGBA lookup tables, built from control points. */
2+
3+function lut(stops: string[]): Uint8Array<ArrayBuffer> {
4+ const rgb = stops.map((h) => [1, 3, 5].map((i) => parseInt(h.slice(i, i + 2), 16)));
5+ const out = new Uint8Array(256 * 4);
6+ for (let i = 0; i < 256; i++) {
7+ const t = (i / 255) * (rgb.length - 1);
8+ const j = Math.min(rgb.length - 2, Math.floor(t));
9+ const f = t - j;
10+ for (let c = 0; c < 3; c++) out[4 * i + c] = Math.round(rgb[j][c] * (1 - f) + rgb[j + 1][c] * f);
11+ out[4 * i + 3] = 255;
12+ }
13+ return out;
14+}
15+
16+/** ColorBrewer RdBu, reversed: blue for negative pressure, red for positive. */
17+export const DIVERGING = lut(['#053061', '#2166ac', '#4393c3', '#92c5de', '#d1e5f0', '#f7f7f7', '#fddbc7', '#f4a582', '#d6604d', '#b2182b', '#67001f']);
18+
19+/** Viridis, for the level in dB. */
20+export const SEQUENTIAL = lut(['#440154', '#482878', '#3e4989', '#31688e', '#26828e', '#1f9e89', '#35b779', '#6ece58', '#b5de2b', '#fde725']);
21+
22+/** A vertical colorbar into a canvas. */
23+export function drawColorbar(canvas: HTMLCanvasElement, map: Uint8Array) {
24+ const ctx = canvas.getContext('2d')!;
25+ const img = ctx.createImageData(1, 256);
26+ for (let i = 0; i < 256; i++) img.data.set(map.subarray(4 * (255 - i), 4 * (255 - i) + 4), 4 * i);
27+ const tmp = document.createElement('canvas');
28+ tmp.width = 1;
29+ tmp.height = 256;
30+ tmp.getContext('2d')!.putImageData(img, 0, 0);
31+ ctx.imageSmoothingEnabled = true;
32+ ctx.clearRect(0, 0, canvas.width, canvas.height);
33+ ctx.drawImage(tmp, 0, 0, canvas.width, canvas.height);
34+}
src/view/convergence.tsadded+68−0View file
@@ -0,0 +1,68 @@
1+/** GMRES relative residual against iteration, on a log scale. */
2+export function drawConvergence(canvas: HTMLCanvasElement, series: { label: string; res: number[]; color: string }[], eps: number) {
3+ const dpr = window.devicePixelRatio || 1;
4+ const W = canvas.clientWidth || 300;
5+ const H = canvas.clientHeight || 140;
6+ canvas.width = Math.round(W * dpr);
7+ canvas.height = Math.round(H * dpr);
8+ const ctx = canvas.getContext('2d')!;
9+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
10+ ctx.clearRect(0, 0, W, H);
11+ const style = getComputedStyle(canvas);
12+ const ink = style.getPropertyValue('--ink-2').trim() || '#888';
13+ const line = style.getPropertyValue('--line').trim() || '#ccc';
14+ const L = 44;
15+ const R = 8;
16+ const T = 8;
17+ const B = 22;
18+ const nmax = Math.max(10, ...series.map((s) => s.res.length));
19+ const lo = Math.floor(Math.log10(Math.min(eps, ...series.flatMap((s) => s.res)) * 0.5));
20+ const hi = 0;
21+ const X = (i: number) => L + ((W - L - R) * i) / nmax;
22+ const Y = (r: number) => T + ((H - T - B) * (hi - Math.log10(r))) / (hi - lo);
23+ ctx.font = '11px system-ui, sans-serif';
24+ ctx.fillStyle = ink;
25+ ctx.strokeStyle = line;
26+ ctx.lineWidth = 1;
27+ for (let e = lo; e <= hi; e++) {
28+ if ((e - lo) % Math.max(1, Math.ceil((hi - lo) / 5)) !== 0) continue;
29+ const y = Y(10 ** e);
30+ ctx.beginPath();
31+ ctx.moveTo(L, y);
32+ ctx.lineTo(W - R, y);
33+ ctx.stroke();
34+ ctx.textAlign = 'right';
35+ ctx.textBaseline = 'middle';
36+ ctx.fillText(`1e${e}`, L - 4, y);
37+ }
38+ ctx.textAlign = 'center';
39+ ctx.textBaseline = 'top';
40+ ctx.fillText('GMRES iteration', (L + W - R) / 2, H - B + 6);
41+ // tolerance
42+ ctx.setLineDash([4, 3]);
43+ ctx.beginPath();
44+ ctx.moveTo(L, Y(eps));
45+ ctx.lineTo(W - R, Y(eps));
46+ ctx.stroke();
47+ ctx.setLineDash([]);
48+ for (const s of series) {
49+ if (!s.res.length) continue;
50+ ctx.strokeStyle = s.color;
51+ ctx.lineWidth = 2;
52+ ctx.beginPath();
53+ ctx.moveTo(X(0), Y(1));
54+ s.res.forEach((r, i) => ctx.lineTo(X(i + 1), Y(Math.max(r, 10 ** lo))));
55+ ctx.stroke();
56+ }
57+ let ly = T + 4;
58+ for (const s of series) {
59+ if (!s.res.length) continue;
60+ ctx.fillStyle = s.color;
61+ ctx.fillRect(W - R - 110, ly + 3, 10, 3);
62+ ctx.fillStyle = ink;
63+ ctx.textAlign = 'left';
64+ ctx.textBaseline = 'top';
65+ ctx.fillText(s.label, W - R - 96, ly - 2);
66+ ly += 14;
67+ }
68+}
src/view/field.tsadded+60−0View file
@@ -0,0 +1,60 @@
1+/**
2+ * The solved field as the views need it, and how a complex value becomes a
3+ * color: the pressure at phase t, Re(u exp(-i omega t)) = Re u cos wt + Im u sin wt,
4+ * on a symmetric diverging scale, or the level 20 log10(|u| / ref) in dB.
5+ */
6+import type { Slice } from '../problem.ts';
7+
8+export type Mode = 'pressure' | 'level';
9+
10+export interface SliceField {
11+ slice: Slice;
12+ /** interleaved complex, NaN outside the room */
13+ u: Float64Array;
14+}
15+
16+export interface Field {
17+ k: number;
18+ beta: number;
19+ x0: [number, number, number];
20+ a: number;
21+ p: number;
22+ ngrid: number;
23+ grid: Float64Array;
24+ slices: SliceField[];
25+ /** wall field at the nodes, interleaved complex */
26+ uWall: Float64Array;
27+ /** color limit for the pressure: 99th percentile of |u| more than half a wavelength from the source */
28+ umax: number;
29+ /** reference for the level: the median of |u| on the same points */
30+ ref: number;
31+}
32+
33+export const DB_MIN = -30;
34+export const DB_MAX = 15;
35+
36+/** Color limits from the slices, as plot.py computes them. */
37+export function colorLimits(slices: SliceField[], x0: ArrayLike<number>, k: number): { umax: number; ref: number } {
38+ const vals: number[] = [];
39+ const rmin = Math.PI / k;
40+ for (const { slice, u } of slices) {
41+ const n = slice.mask.length;
42+ for (let q = 0; q < n; q++) {
43+ if (!slice.mask[q]) continue;
44+ const r = Math.hypot(slice.pts[q] - x0[0], slice.pts[n + q] - x0[1], slice.pts[2 * n + q] - x0[2]);
45+ if (r <= rmin) continue;
46+ const v = Math.hypot(u[2 * q], u[2 * q + 1]);
47+ if (Number.isFinite(v)) vals.push(v);
48+ }
49+ }
50+ if (!vals.length) return { umax: 1, ref: 1 };
51+ vals.sort((x, y) => x - y);
52+ return { umax: vals[Math.floor(0.99 * (vals.length - 1))], ref: vals[Math.floor(0.5 * (vals.length - 1))] };
53+}
54+
55+/** Position in [0, 1] on the colormap of a complex value, or NaN. */
56+export function colorParam(re: number, im: number, mode: Mode, c: number, s: number, umax: number, ref: number): number {
57+ if (mode === 'pressure') return 0.5 + (0.5 * (re * c + im * s)) / umax;
58+ const db = 20 * Math.log10(Math.hypot(re, im) / ref);
59+ return (db - DB_MIN) / (DB_MAX - DB_MIN);
60+}
src/view/room3d.tsadded+486−0View file
@@ -0,0 +1,486 @@
1+/**
2+ * A WebGPU view of the room: the walls colored by the wall field, the three
3+ * planes through the source, and the source itself. Walls whose outward
4+ * normal faces the camera are cut away in the fragment shader, so the view
5+ * looks into the room through the nearer walls. Drag to orbit, scroll or
6+ * pinch to zoom.
7+ */
8+import { basis2d, type Surface } from '../surface.ts';
9+import { DIVERGING, SEQUENTIAL } from './colormaps.ts';
10+import { DB_MAX, DB_MIN, type Field, type Mode } from './field.ts';
11+
12+/** Multisampling; 1 when the query string has ?nomsaa (SwiftShader's GPU process does not survive the resolve in headless Chrome). */
13+const SAMPLES = typeof location !== 'undefined' && location.search.includes('nomsaa') ? 1 : 4;
14+
15+const SHARED = /* wgsl */ `
16+struct View {
17+ viewProj: mat4x4f,
18+ eye: vec4f,
19+ p1: vec4f, // cos wt, sin wt, umax, ref
20+ p2: vec4f, // mode (0 pressure, 1 level), db min, db max, 0
21+}
22+@group(0) @binding(0) var<uniform> V: View;
23+@group(0) @binding(1) var cmapD: texture_2d<f32>;
24+@group(0) @binding(2) var cmapS: texture_2d<f32>;
25+@group(0) @binding(3) var samp: sampler;
26+
27+fn colorOf(re: f32, im: f32) -> vec3f {
28+ if (V.p2.x < 0.5) {
29+ let t = clamp(0.5 + 0.5 * (re * V.p1.x + im * V.p1.y) / V.p1.z, 0.0, 1.0);
30+ return textureSampleLevel(cmapD, samp, vec2f(t, 0.5), 0.0).rgb;
31+ }
32+ let db = 20.0 * log2(max(length(vec2f(re, im)) / V.p1.w, 1e-30)) * 0.30102999566;
33+ let t = clamp((db - V.p2.y) / (V.p2.z - V.p2.y), 0.0, 1.0);
34+ return textureSampleLevel(cmapS, samp, vec2f(t, 0.5), 0.0).rgb;
35+}
36+`;
37+
38+const WALL = /* wgsl */ `
39+${SHARED}
40+struct VOut {
41+ @builtin(position) pos: vec4f,
42+ @location(0) world: vec3f,
43+ @location(1) nrm: vec3f,
44+ @location(2) u: vec2f,
45+}
46+@vertex fn vs(@location(0) p: vec3f, @location(1) n: vec3f, @location(2) u: vec2f) -> VOut {
47+ var o: VOut;
48+ o.pos = V.viewProj * vec4f(p, 1.0);
49+ o.world = p;
50+ o.nrm = n;
51+ o.u = u;
52+ return o;
53+}
54+@fragment fn fs(i: VOut) -> @location(0) vec4f {
55+ let toEye = normalize(V.eye.xyz - i.world);
56+ let n = normalize(i.nrm);
57+ if (dot(n, toEye) > 0.0) { discard; }
58+ let shade = 0.6 + 0.4 * abs(dot(n, toEye));
59+ return vec4f(colorOf(i.u.x, i.u.y) * shade, 1.0);
60+}
61+`;
62+
63+const PLANE = /* wgsl */ `
64+${SHARED}
65+@group(0) @binding(4) var field: texture_2d<f32>;
66+struct VOut {
67+ @builtin(position) pos: vec4f,
68+ @location(0) uv: vec2f,
69+}
70+@vertex fn vs(@location(0) p: vec3f, @location(1) uv: vec2f) -> VOut {
71+ var o: VOut;
72+ o.pos = V.viewProj * vec4f(p, 1.0);
73+ o.uv = uv;
74+ return o;
75+}
76+@fragment fn fs(i: VOut) -> @location(0) vec4f {
77+ // bilinear by hand: r32float textures are not filterable everywhere
78+ let dim = vec2f(textureDimensions(field));
79+ let x = clamp(i.uv * (dim - 1.0), vec2f(0.0), dim - 1.0);
80+ let x0 = vec2i(floor(x));
81+ let x1 = min(x0 + 1, vec2i(dim) - 1);
82+ let f = x - vec2f(x0);
83+ let a = textureLoad(field, x0, 0).xy;
84+ let b = textureLoad(field, vec2i(x1.x, x0.y), 0).xy;
85+ let c = textureLoad(field, vec2i(x0.x, x1.y), 0).xy;
86+ let d = textureLoad(field, x1, 0).xy;
87+ // NaN marks points outside the room
88+ if (any(a != a) || any(b != b) || any(c != c) || any(d != d)) { discard; }
89+ let u = mix(mix(a, b, f.x), mix(c, d, f.x), f.y);
90+ return vec4f(colorOf(u.x, u.y), 1.0);
91+}
92+`;
93+
94+const MARKER = /* wgsl */ `
95+${SHARED}
96+@vertex fn vs(@location(0) p: vec3f) -> @builtin(position) vec4f {
97+ return V.viewProj * vec4f(p, 1.0);
98+}
99+@fragment fn fs() -> @location(0) vec4f {
100+ return vec4f(0.05, 0.05, 0.05, 1.0);
101+}
102+`;
103+
104+type Mat4 = Float32Array;
105+const BG = { r: 0.043, g: 0.055, b: 0.071, a: 1 };
106+
107+function perspective(fovy: number, aspect: number, near: number, far: number): Mat4 {
108+ const f = 1 / Math.tan(fovy / 2);
109+ const m = new Float32Array(16);
110+ m[0] = f / aspect;
111+ m[5] = f;
112+ m[10] = far / (near - far);
113+ m[11] = -1;
114+ m[14] = (near * far) / (near - far);
115+ return m;
116+}
117+
118+function lookAt(eye: number[], target: number[], up: number[]): Mat4 {
119+ const z = norm(sub(eye, target));
120+ const x = norm(cross(up, z));
121+ const y = cross(z, x);
122+ const m = new Float32Array(16);
123+ m.set([x[0], y[0], z[0], 0, x[1], y[1], z[1], 0, x[2], y[2], z[2], 0, -dot(x, eye), -dot(y, eye), -dot(z, eye), 1]);
124+ return m;
125+}
126+
127+function mul(a: Mat4, b: Mat4): Mat4 {
128+ const m = new Float32Array(16);
129+ for (let c = 0; c < 4; c++) for (let r = 0; r < 4; r++) {
130+ let s = 0;
131+ for (let k = 0; k < 4; k++) s += a[k * 4 + r] * b[c * 4 + k];
132+ m[c * 4 + r] = s;
133+ }
134+ return m;
135+}
136+
137+const sub = (a: number[], b: number[]) => [a[0] - b[0], a[1] - b[1], a[2] - b[2]];
138+const dot = (a: number[], b: number[]) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
139+const cross = (a: number[], b: number[]) => [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
140+const norm = (a: number[]) => {
141+ const l = Math.hypot(a[0], a[1], a[2]);
142+ return [a[0] / l, a[1] / l, a[2] / l];
143+};
144+
145+export interface Show {
146+ walls: boolean;
147+ planes: [boolean, boolean, boolean];
148+}
149+
150+export class RoomView {
151+ private ctx: GPUCanvasContext;
152+ private format: GPUTextureFormat;
153+ private wallPipe: GPURenderPipeline;
154+ private planePipe: GPURenderPipeline;
155+ private markerPipe: GPURenderPipeline;
156+ private ubuf: GPUBuffer;
157+ private cmapD: GPUTexture;
158+ private cmapS: GPUTexture;
159+ private sampler: GPUSampler;
160+ private baseGroup: GPUBindGroup;
161+ private wall: { vbuf: GPUBuffer; ibuf: GPUBuffer; count: number } | null = null;
162+ private planes: { vbuf: GPUBuffer; tex: GPUTexture; group: GPUBindGroup }[] = [];
163+ private marker: { vbuf: GPUBuffer; count: number } | null = null;
164+ private msaa: GPUTexture | null = null;
165+ private depth: GPUTexture | null = null;
166+ private theta = 0.75;
167+ private phi = 0.95;
168+ private dist = 5.2;
169+ private field: Field | null = null;
170+ private dirty = true;
171+
172+ constructor(
173+ private device: GPUDevice,
174+ private canvas: HTMLCanvasElement,
175+ ) {
176+ this.ctx = canvas.getContext('webgpu') as GPUCanvasContext;
177+ this.format = navigator.gpu.getPreferredCanvasFormat();
178+ this.ctx.configure({ device, format: this.format, alphaMode: 'opaque' });
179+ this.ubuf = device.createBuffer({ size: 112, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
180+ const mkmap = (lut: Uint8Array<ArrayBuffer>) => {
181+ const t = device.createTexture({ size: [256, 1], format: 'rgba8unorm', usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST });
182+ device.queue.writeTexture({ texture: t }, lut, { bytesPerRow: 1024 }, [256, 1]);
183+ return t;
184+ };
185+ this.cmapD = mkmap(DIVERGING);
186+ this.cmapS = mkmap(SEQUENTIAL);
187+ this.sampler = device.createSampler({ magFilter: 'linear', minFilter: 'linear', addressModeU: 'clamp-to-edge' });
188+
189+ const baseEntries: GPUBindGroupLayoutEntry[] = [
190+ { binding: 0, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: 'uniform' } },
191+ { binding: 1, visibility: GPUShaderStage.FRAGMENT, texture: { sampleType: 'float' } },
192+ { binding: 2, visibility: GPUShaderStage.FRAGMENT, texture: { sampleType: 'float' } },
193+ { binding: 3, visibility: GPUShaderStage.FRAGMENT, sampler: { type: 'filtering' } },
194+ ];
195+ const baseLayout = device.createBindGroupLayout({ entries: baseEntries });
196+ const planeLayout = device.createBindGroupLayout({
197+ entries: [...baseEntries, { binding: 4, visibility: GPUShaderStage.FRAGMENT, texture: { sampleType: 'unfilterable-float' } }],
198+ });
199+ this.baseGroup = device.createBindGroup({ layout: baseLayout, entries: this.baseEntries() });
200+ const common = (code: string, layout: GPUBindGroupLayout, buffers: GPUVertexBufferLayout[], topology: GPUPrimitiveTopology = 'triangle-list'): GPURenderPipeline => {
201+ const module = device.createShaderModule({ code });
202+ return device.createRenderPipeline({
203+ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }),
204+ vertex: { module, entryPoint: 'vs', buffers },
205+ fragment: { module, entryPoint: 'fs', targets: [{ format: this.format }] },
206+ primitive: { topology, cullMode: 'none' },
207+ depthStencil: { format: 'depth24plus', depthWriteEnabled: true, depthCompare: 'less' },
208+ multisample: { count: SAMPLES },
209+ });
210+ };
211+ this.wallPipe = common(WALL, baseLayout, [
212+ {
213+ arrayStride: 32,
214+ attributes: [
215+ { shaderLocation: 0, offset: 0, format: 'float32x3' },
216+ { shaderLocation: 1, offset: 12, format: 'float32x3' },
217+ { shaderLocation: 2, offset: 24, format: 'float32x2' },
218+ ],
219+ },
220+ ]);
221+ this.planePipe = common(PLANE, planeLayout, [
222+ {
223+ arrayStride: 20,
224+ attributes: [
225+ { shaderLocation: 0, offset: 0, format: 'float32x3' },
226+ { shaderLocation: 1, offset: 12, format: 'float32x2' },
227+ ],
228+ },
229+ ]);
230+ this.planeLayout = planeLayout;
231+ this.markerPipe = common(MARKER, baseLayout, [{ arrayStride: 12, attributes: [{ shaderLocation: 0, offset: 0, format: 'float32x3' }] }]);
232+ this.attachControls();
233+ new ResizeObserver(() => {
234+ this.dirty = true;
235+ }).observe(canvas);
236+ }
237+
238+ private planeLayout: GPUBindGroupLayout;
239+
240+ private baseEntries(): GPUBindGroupEntry[] {
241+ return [
242+ { binding: 0, resource: { buffer: this.ubuf } },
243+ { binding: 1, resource: this.cmapD.createView() },
244+ { binding: 2, resource: this.cmapS.createView() },
245+ { binding: 3, resource: this.sampler },
246+ ];
247+ }
248+
249+ private attachControls() {
250+ const c = this.canvas;
251+ const pointers = new Map<number, { x: number; y: number }>();
252+ let pinch = 0;
253+ c.addEventListener('pointerdown', (e) => {
254+ c.setPointerCapture(e.pointerId);
255+ pointers.set(e.pointerId, { x: e.clientX, y: e.clientY });
256+ if (pointers.size === 2) {
257+ const [a, b] = [...pointers.values()];
258+ pinch = Math.hypot(a.x - b.x, a.y - b.y);
259+ }
260+ });
261+ c.addEventListener('pointermove', (e) => {
262+ const p = pointers.get(e.pointerId);
263+ if (!p) return;
264+ if (pointers.size === 1) {
265+ this.theta -= (e.clientX - p.x) * 0.008;
266+ this.phi = Math.min(Math.PI - 0.05, Math.max(0.05, this.phi - (e.clientY - p.y) * 0.008));
267+ }
268+ p.x = e.clientX;
269+ p.y = e.clientY;
270+ if (pointers.size === 2) {
271+ const [a, b] = [...pointers.values()];
272+ const d = Math.hypot(a.x - b.x, a.y - b.y);
273+ if (pinch > 0) this.zoom(pinch / d);
274+ pinch = d;
275+ }
276+ this.dirty = true;
277+ });
278+ const up = (e: PointerEvent) => {
279+ pointers.delete(e.pointerId);
280+ pinch = 0;
281+ };
282+ c.addEventListener('pointerup', up);
283+ c.addEventListener('pointercancel', up);
284+ c.addEventListener(
285+ 'wheel',
286+ (e) => {
287+ e.preventDefault();
288+ this.zoom(Math.exp(e.deltaY * 0.001));
289+ },
290+ { passive: false },
291+ );
292+ }
293+
294+ private zoom(f: number) {
295+ this.dist = Math.min(12, Math.max(2.2, this.dist * f));
296+ this.dirty = true;
297+ }
298+
299+ /** Drop the displayed field, for while a new one is computed. */
300+ clear() {
301+ this.wall?.vbuf.destroy();
302+ this.wall?.ibuf.destroy();
303+ this.wall = null;
304+ for (const p of this.planes) {
305+ p.vbuf.destroy();
306+ p.tex.destroy();
307+ }
308+ this.planes = [];
309+ this.marker?.vbuf.destroy();
310+ this.marker = null;
311+ this.field = null;
312+ this.dirty = true;
313+ }
314+
315+ /** The wall mesh: every patch sampled on an nv x nv grid, with the wall field interpolated by the patch's Legendre expansion. */
316+ setField(surf: Surface, f: Field, nv = 9) {
317+ const { device } = this;
318+ this.field = f;
319+ const npols = surf.npols;
320+ const npatches = surf.npatches;
321+ const npts = surf.npts;
322+ const uv = new Float64Array(2 * nv * nv);
323+ for (let a = 0; a < nv; a++) for (let b = 0; b < nv; b++) {
324+ uv[a * nv + b] = -1 + (2 * b) / (nv - 1);
325+ uv[nv * nv + a * nv + b] = -1 + (2 * a) / (nv - 1);
326+ }
327+ const B = basis2d(uv, surf.norder); // (npols, nv^2)
328+ const ng = nv * nv;
329+ const verts = new Float32Array(8 * ng * npatches);
330+ const idx = new Uint32Array(6 * (nv - 1) * (nv - 1) * npatches);
331+ const cr = new Float64Array(npols);
332+ const ci = new Float64Array(npols);
333+ let ii = 0;
334+ for (let j = 0; j < npatches; j++) {
335+ const g = surf.evalGeometry(j, uv);
336+ for (let n = 0; n < npols; n++) {
337+ let sr = 0;
338+ let si = 0;
339+ for (let l = 0; l < npols; l++) {
340+ const U = surf.umat[n * npols + l];
341+ sr += U * f.uWall[2 * (j * npols + l)];
342+ si += U * f.uWall[2 * (j * npols + l) + 1];
343+ }
344+ cr[n] = sr;
345+ ci[n] = si;
346+ }
347+ for (let q = 0; q < ng; q++) {
348+ const o = 8 * (j * ng + q);
349+ const ru = [g[3 * ng + q], g[4 * ng + q], g[5 * ng + q]];
350+ const rv = [g[6 * ng + q], g[7 * ng + q], g[8 * ng + q]];
351+ const nn = norm(cross(ru, rv));
352+ let ur = 0;
353+ let ui = 0;
354+ for (let n = 0; n < npols; n++) {
355+ ur += cr[n] * B[n * ng + q];
356+ ui += ci[n] * B[n * ng + q];
357+ }
358+ verts.set([g[q], g[ng + q], g[2 * ng + q], nn[0], nn[1], nn[2], ur, ui], o);
359+ }
360+ for (let a = 0; a < nv - 1; a++) for (let b = 0; b < nv - 1; b++) {
361+ const v0 = j * ng + a * nv + b;
362+ idx.set([v0, v0 + 1, v0 + nv + 1, v0, v0 + nv + 1, v0 + nv], ii);
363+ ii += 6;
364+ }
365+ }
366+ void npts;
367+ this.wall?.vbuf.destroy();
368+ this.wall?.ibuf.destroy();
369+ const vbuf = device.createBuffer({ size: verts.byteLength, usage: GPUBufferUsage.VERTEX | GPUBufferUsage.COPY_DST });
370+ device.queue.writeBuffer(vbuf, 0, verts);
371+ const ibuf = device.createBuffer({ size: idx.byteLength, usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
372+ device.queue.writeBuffer(ibuf, 0, idx);
373+ this.wall = { vbuf, ibuf, count: idx.length };
374+
375+ // the three planes, as textures of (Re u, Im u) with NaN outside the room
376+ for (const p of this.planes) {
377+ p.vbuf.destroy();
378+ p.tex.destroy();
379+ }
380+ this.planes = f.slices.map(({ slice, u }) => {
381+ const n = f.ngrid;
382+ const tex = device.createTexture({ size: [n, n], format: 'rg32float', usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST });
383+ const data = new Float32Array(2 * n * n);
384+ for (let q = 0; q < n * n; q++) {
385+ data[2 * q] = slice.mask[q] ? u[2 * q] : NaN;
386+ data[2 * q + 1] = slice.mask[q] ? u[2 * q + 1] : NaN;
387+ }
388+ device.queue.writeTexture({ texture: tex }, data, { bytesPerRow: 8 * n }, [n, n]);
389+ const [i, jx] = slice.axes;
390+ const corner = (s: number, t: number) => {
391+ const x = [0, 0, 0];
392+ x[i] = -f.a + 2 * f.a * s;
393+ x[jx] = -f.a + 2 * f.a * t;
394+ x[slice.fixed] = f.x0[slice.fixed];
395+ return [...x, s, t];
396+ };
397+ const quad = new Float32Array([...corner(0, 0), ...corner(1, 0), ...corner(1, 1), ...corner(0, 0), ...corner(1, 1), ...corner(0, 1)]);
398+ const vb = device.createBuffer({ size: quad.byteLength, usage: GPUBufferUsage.VERTEX | GPUBufferUsage.COPY_DST });
399+ device.queue.writeBuffer(vb, 0, quad);
400+ const group = device.createBindGroup({ layout: this.planeLayout, entries: [...this.baseEntries(), { binding: 4, resource: tex.createView() }] });
401+ return { vbuf: vb, tex, group };
402+ });
403+
404+ // the source: a small octahedron
405+ const r = 0.035;
406+ const [x, y, z] = f.x0;
407+ const P = [
408+ [x + r, y, z], [x - r, y, z], [x, y + r, z], [x, y - r, z], [x, y, z + r], [x, y, z - r],
409+ ];
410+ const faces = [[0, 2, 4], [2, 1, 4], [1, 3, 4], [3, 0, 4], [2, 0, 5], [1, 2, 5], [3, 1, 5], [0, 3, 5]];
411+ const mv = new Float32Array(faces.flat().flatMap((k) => P[k]));
412+ this.marker?.vbuf.destroy();
413+ const mb = device.createBuffer({ size: mv.byteLength, usage: GPUBufferUsage.VERTEX | GPUBufferUsage.COPY_DST });
414+ device.queue.writeBuffer(mb, 0, mv);
415+ this.marker = { vbuf: mb, count: faces.length * 3 };
416+ this.dirty = true;
417+ }
418+
419+ private targets() {
420+ const dpr = Math.min(2, window.devicePixelRatio || 1);
421+ const w = Math.max(1, Math.round(this.canvas.clientWidth * dpr));
422+ const h = Math.max(1, Math.round(this.canvas.clientHeight * dpr));
423+ if (this.canvas.width !== w || this.canvas.height !== h || !this.depth) {
424+ this.canvas.width = w;
425+ this.canvas.height = h;
426+ this.msaa?.destroy();
427+ this.depth?.destroy();
428+ this.msaa = this.device.createTexture({ size: [w, h], format: this.format, sampleCount: SAMPLES, usage: GPUTextureUsage.RENDER_ATTACHMENT });
429+ this.depth = this.device.createTexture({ size: [w, h], format: 'depth24plus', sampleCount: SAMPLES, usage: GPUTextureUsage.RENDER_ATTACHMENT });
430+ }
431+ return { w, h };
432+ }
433+
434+ draw(mode: Mode, phase: number, show: Show, animating: boolean) {
435+ if (!animating && !this.dirty) return;
436+ this.dirty = false;
437+ const { device } = this;
438+ const { w, h } = this.targets();
439+ const eye = [this.dist * Math.sin(this.phi) * Math.cos(this.theta), this.dist * Math.sin(this.phi) * Math.sin(this.theta), this.dist * Math.cos(this.phi)];
440+ const vp = mul(perspective((35 * Math.PI) / 180, w / h, 0.1, 50), lookAt(eye, [0, 0, 0], [0, 0, 1]));
441+ const f = this.field;
442+ const u = new Float32Array(28);
443+ u.set(vp, 0);
444+ u.set([eye[0], eye[1], eye[2], 1], 16);
445+ u.set([Math.cos(phase), Math.sin(phase), f?.umax ?? 1, f?.ref ?? 1], 20);
446+ u.set([mode === 'pressure' ? 0 : 1, DB_MIN, DB_MAX, 0], 24);
447+ device.queue.writeBuffer(this.ubuf, 0, u);
448+ const enc = device.createCommandEncoder();
449+ const pass = enc.beginRenderPass({
450+ colorAttachments: [
451+ SAMPLES > 1
452+ ? { view: this.msaa!.createView(), resolveTarget: this.ctx.getCurrentTexture().createView(), clearValue: BG, loadOp: 'clear', storeOp: 'discard' }
453+ : { view: this.ctx.getCurrentTexture().createView(), clearValue: BG, loadOp: 'clear', storeOp: 'store' },
454+ ],
455+ depthStencilAttachment: { view: this.depth!.createView(), depthClearValue: 1, depthLoadOp: 'clear', depthStoreOp: 'discard' },
456+ });
457+ if (this.wall && show.walls) {
458+ pass.setPipeline(this.wallPipe);
459+ pass.setBindGroup(0, this.baseGroup);
460+ pass.setVertexBuffer(0, this.wall.vbuf);
461+ pass.setIndexBuffer(this.wall.ibuf, 'uint32');
462+ pass.drawIndexed(this.wall.count);
463+ }
464+ if (this.planes.length) {
465+ pass.setPipeline(this.planePipe);
466+ this.planes.forEach((p, s) => {
467+ if (!show.planes[s]) return;
468+ pass.setBindGroup(0, p.group);
469+ pass.setVertexBuffer(0, p.vbuf);
470+ pass.draw(6);
471+ });
472+ }
473+ if (this.marker) {
474+ pass.setPipeline(this.markerPipe);
475+ pass.setBindGroup(0, this.baseGroup);
476+ pass.setVertexBuffer(0, this.marker.vbuf);
477+ pass.draw(this.marker.count);
478+ }
479+ pass.end();
480+ device.queue.submit([enc.finish()]);
481+ }
482+
483+ invalidate() {
484+ this.dirty = true;
485+ }
486+}
src/view/slices.tsadded+146−0View file
@@ -0,0 +1,146 @@
1+/**
2+ * The three planes through the source as 2D images, redrawn every frame so
3+ * that the pressure can be animated over one period. Clicking a plane moves
4+ * the source to that point (its third coordinate is the plane's).
5+ */
6+import { DIVERGING, SEQUENTIAL } from './colormaps.ts';
7+import { colorParam, type Field, type Mode } from './field.ts';
8+
9+const AXIS = ['x', 'y', 'z'];
10+
11+export class SliceView {
12+ private canvases: HTMLCanvasElement[] = [];
13+ private titles: HTMLElement[] = [];
14+ private img: ImageData[] = [];
15+ private field: Field | null = null;
16+
17+ constructor(
18+ container: HTMLElement,
19+ private onPick: (x: [number, number, number]) => void,
20+ private onHover: (text: string) => void,
21+ ) {
22+ for (let s = 0; s < 3; s++) {
23+ const fig = document.createElement('figure');
24+ fig.className = 'slice';
25+ const cap = document.createElement('figcaption');
26+ const cv = document.createElement('canvas');
27+ cv.width = 300;
28+ cv.height = 300;
29+ fig.append(cv, cap);
30+ container.append(fig);
31+ this.canvases.push(cv);
32+ this.titles.push(cap);
33+ cv.addEventListener('click', (ev) => this.pick(s, ev));
34+ cv.addEventListener('pointermove', (ev) => this.hover(s, ev));
35+ cv.addEventListener('pointerleave', () => this.onHover(''));
36+ }
37+ }
38+
39+ /** Blank the planes, for while a new field is computed. */
40+ clear() {
41+ this.field = null;
42+ this.img = [];
43+ for (const [s, cv] of this.canvases.entries()) {
44+ cv.getContext('2d')!.clearRect(0, 0, cv.width, cv.height);
45+ this.titles[s].textContent = '';
46+ }
47+ this.onHover('');
48+ }
49+
50+ setField(f: Field) {
51+ this.field = f;
52+ this.img = f.slices.map(() => new ImageData(f.ngrid, f.ngrid));
53+ f.slices.forEach(({ slice }, s) => {
54+ const [i, j] = slice.axes;
55+ this.titles[s].textContent = `${slice.name}-plane at ${AXIS[slice.fixed]} = ${f.x0[slice.fixed].toFixed(2)} (${AXIS[i]} across, ${AXIS[j]} up)`;
56+ });
57+ }
58+
59+ /** Coordinates of a pointer event in the room, for slice s. */
60+ private at(s: number, ev: PointerEvent | MouseEvent): { x: [number, number, number]; q: number } | null {
61+ const f = this.field;
62+ if (!f) return null;
63+ const cv = this.canvases[s];
64+ const r = cv.getBoundingClientRect();
65+ const fx = (ev.clientX - r.left) / r.width;
66+ const fy = 1 - (ev.clientY - r.top) / r.height;
67+ const { slice } = f.slices[s];
68+ const x: [number, number, number] = [...f.x0] as [number, number, number];
69+ x[slice.axes[0]] = -f.a + 2 * f.a * fx;
70+ x[slice.axes[1]] = -f.a + 2 * f.a * fy;
71+ const col = Math.round(fx * (f.ngrid - 1));
72+ const row = Math.round(fy * (f.ngrid - 1));
73+ if (col < 0 || row < 0 || col >= f.ngrid || row >= f.ngrid) return null;
74+ return { x, q: row * f.ngrid + col };
75+ }
76+
77+ private pick(s: number, ev: MouseEvent) {
78+ const p = this.at(s, ev);
79+ if (p) this.onPick(p.x);
80+ }
81+
82+ private hover(s: number, ev: PointerEvent) {
83+ const p = this.at(s, ev);
84+ const f = this.field;
85+ if (!p || !f) return;
86+ const { u, slice } = f.slices[s];
87+ if (!slice.mask[p.q]) {
88+ this.onHover(`(${p.x.map((v) => v.toFixed(2)).join(', ')}): outside the room`);
89+ return;
90+ }
91+ const re = u[2 * p.q];
92+ const im = u[2 * p.q + 1];
93+ const db = 20 * Math.log10(Math.hypot(re, im) / f.ref);
94+ this.onHover(`(${p.x.map((v) => v.toFixed(2)).join(', ')}): u = ${re.toExponential(3)} ${im < 0 ? '−' : '+'} ${Math.abs(im).toExponential(3)}i, level ${db.toFixed(1)} dB · click to move the source here`);
95+ }
96+
97+ draw(mode: Mode, phase: number) {
98+ const f = this.field;
99+ if (!f) return;
100+ const c = Math.cos(phase);
101+ const sn = Math.sin(phase);
102+ const map = mode === 'pressure' ? DIVERGING : SEQUENTIAL;
103+ f.slices.forEach(({ slice, u }, s) => {
104+ const img = this.img[s];
105+ const n = f.ngrid;
106+ for (let row = 0; row < n; row++) {
107+ for (let col = 0; col < n; col++) {
108+ const q = row * n + col;
109+ const o = 4 * ((n - 1 - row) * n + col);
110+ if (!slice.mask[q]) {
111+ img.data[o + 3] = 0;
112+ continue;
113+ }
114+ const t = colorParam(u[2 * q], u[2 * q + 1], mode, c, sn, f.umax, f.ref);
115+ const i = Math.max(0, Math.min(255, Math.round(t * 255))) * 4;
116+ img.data[o] = map[i];
117+ img.data[o + 1] = map[i + 1];
118+ img.data[o + 2] = map[i + 2];
119+ img.data[o + 3] = 255;
120+ }
121+ }
122+ const cv = this.canvases[s];
123+ const ctx = cv.getContext('2d')!;
124+ const tmp = new OffscreenCanvas(n, n);
125+ tmp.getContext('2d')!.putImageData(img, 0, 0);
126+ ctx.clearRect(0, 0, cv.width, cv.height);
127+ ctx.imageSmoothingEnabled = true;
128+ ctx.imageSmoothingQuality = 'high';
129+ ctx.drawImage(tmp, 0, 0, cv.width, cv.height);
130+ // the source
131+ const [i, j] = slice.axes;
132+ const px = ((f.x0[i] + f.a) / (2 * f.a)) * cv.width;
133+ const py = (1 - (f.x0[j] + f.a) / (2 * f.a)) * cv.height;
134+ ctx.beginPath();
135+ ctx.arc(px, py, 5, 0, 2 * Math.PI);
136+ ctx.lineWidth = 1.5;
137+ ctx.strokeStyle = '#111';
138+ ctx.stroke();
139+ ctx.beginPath();
140+ ctx.arc(px, py, 6.5, 0, 2 * Math.PI);
141+ ctx.strokeStyle = '#fff';
142+ ctx.lineWidth = 1;
143+ ctx.stroke();
144+ });
145+ }
146+}
test.htmladded+17−0View file
@@ -0,0 +1,17 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="utf-8" />
5+ <title>fmm3dbie-webgpu-poc validation</title>
6+ <style>
7+ body { font: 14px/1.4 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace; margin: 16px; }
8+ .ok { color: #1a7f37; }
9+ .fail { color: #cf222e; font-weight: 600; }
10+ </style>
11+ </head>
12+ <body>
13+ <h1>fmm3dbie-webgpu-poc validation suite</h1>
14+ <pre id="log">starting…</pre>
15+ <script type="module" src="/test/test-page.ts"></script>
16+ </body>
17+</html>
test/test-node.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+/**
2+ * CPU-only checks, no GPU: the surface discretization (area and volume of the
3+ * superellipsoid against the exact values) and Green's identity
4+ * S_k[du/dn] - D_k[u] = u/2 at wall nodes using the brute-force reference
5+ * weights. Run with: node --experimental-strip-types test/test-node.ts
6+ */
7+import { roundedCube, superellipsoidVolume } from '../src/surface.ts';
8+import { green, greenDn } from '../src/problem.ts';
9+import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
10+import { leggauss } from '../src/gauss.ts';
11+
12+let failed = 0;
13+const check = (name: string, val: number, tol: number) => {
14+ const ok = val < tol;
15+ if (!ok) failed++;
16+ console.log(`${ok ? 'ok ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`);
17+};
18+
19+// Gauss-Legendre against known values
20+{
21+ const { x, w } = leggauss(5);
22+ check('leggauss(5) node', Math.abs(x[4] - 0.906179845938664), 1e-14);
23+ check('leggauss(5) weight', Math.abs(w[2] - 0.568888888888889), 1e-14);
24+}
25+
26+const a = 1;
27+const p = 8;
28+const surf = roundedCube(a, p, 3, 6);
29+console.log(`surface: ${surf.npatches} patches, ${surf.npts} nodes`);
30+
31+// volume by the divergence theorem: V = (1/3) int x . n dS
32+{
33+ let vol = 0;
34+ let area = 0;
35+ const n = surf.npts;
36+ for (let i = 0; i < n; i++) {
37+ const xn = surf.srcvals[i] * surf.srcvals[9 * n + i] + surf.srcvals[n + i] * surf.srcvals[10 * n + i] + surf.srcvals[2 * n + i] * surf.srcvals[11 * n + i];
38+ vol += (xn * surf.wts[i]) / 3;
39+ area += surf.wts[i];
40+ }
41+ const vex = superellipsoidVolume(a, p);
42+ check('volume by divergence theorem', Math.abs(vol - vex) / vex, 1e-6);
43+ console.log(` area ${area.toFixed(6)}, volume ${vol.toFixed(8)} (exact ${vex.toFixed(8)})`);
44+}
45+
46+// Green's identity at a few nodes of one patch, brute force over all patches
47+{
48+ const k = 6;
49+ const xe = [1.7, 0.6, -0.5];
50+ const n = surf.npts;
51+ const npols = surf.npols;
52+ const xs = surf.srcvals.subarray(0, 3 * n);
53+ const nrm = surf.srcvals.subarray(9 * n, 12 * n);
54+ const u = green(xs, xe, k);
55+ const un = greenDn(xs, nrm, xe, k);
56+ const P = 13;
57+ for (const loc of [0, 24, 47]) {
58+ const i = P * npols + loc;
59+ const xt = [xs[i], xs[n + i], xs[2 * n + i]];
60+ const nt = [nrm[i], nrm[n + i], nrm[2 * n + i]];
61+ let Sr = 0;
62+ let Si = 0;
63+ let Dr = 0;
64+ let Di = 0;
65+ for (let j = 0; j < surf.npatches; j++) {
66+ let I;
67+ if (j === P) I = refSelf(surf, j, surf.uvNodes[i], surf.uvNodes[n + i], xt, nt, k, 200, 40);
68+ else {
69+ 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];
70+ I = refLeaf(surf, j, xt, nt, k, d < 3 ? 5 : 2);
71+ }
72+ const W = momentsToWeights(surf, I);
73+ for (let l = 0; l < npols; l++) {
74+ const q = j * npols + l;
75+ Sr += W[4][2 * l] * un[2 * q] - W[4][2 * l + 1] * un[2 * q + 1];
76+ Si += W[4][2 * l] * un[2 * q + 1] + W[4][2 * l + 1] * un[2 * q];
77+ Dr += W[5][2 * l] * u[2 * q] - W[5][2 * l + 1] * u[2 * q + 1];
78+ Di += W[5][2 * l] * u[2 * q + 1] + W[5][2 * l + 1] * u[2 * q];
79+ }
80+ }
81+ const rr = Sr - Dr - u[2 * i] / 2;
82+ const ri = Si - Di - u[2 * i + 1] / 2;
83+ check(`Green's identity at node ${loc} of patch ${P}`, Math.hypot(rr, ri) / Math.hypot(u[2 * i], u[2 * i + 1]), 1e-5);
84+ }
85+}
86+
87+if (failed) {
88+ console.log(`${failed} check(s) failed`);
89+ process.exit(1);
90+}
91+console.log('all checks passed');
test/test-page.tsadded+146−0View file
@@ -0,0 +1,146 @@
1+/**
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+ */
13+import { roundedCube } from '../src/surface.ts';
14+import { green, greenDn, selftestRhs, selftestPoints, relError, ALPHA, type CVec } from '../src/problem.ts';
15+import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
16+import { requestSolverDevice } from '../src/gpu/device.ts';
17+import { Operator } from '../src/gpu/operator.ts';
18+
19+declare global {
20+ interface Window {
21+ __RESULTS__?: { passed: number; failed: number; fatal?: string; lines: string[] };
22+ }
23+}
24+
25+const logEl = document.getElementById('log')!;
26+const lines: string[] = [];
27+let passed = 0;
28+let failed = 0;
29+const 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+};
37+const 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+};
42+
43+async 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`);
50+
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);
63+
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`);
68+
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+ }
102+
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+ }
120+
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+ }
136+
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+}
141+
142+main().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+});
tsconfig.jsonadded+15−0View file
@@ -0,0 +1,15 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "module": "ESNext",
5+ "moduleResolution": "bundler",
6+ "lib": ["ES2022", "DOM", "DOM.Iterable"],
7+ "types": ["@webgpu/types", "node"],
8+ "strict": true,
9+ "noEmit": true,
10+ "allowImportingTsExtensions": true,
11+ "verbatimModuleSyntax": true,
12+ "skipLibCheck": true
13+ },
14+ "include": ["src", "test"]
15+}
vite.config.tsadded+15−0View file
@@ -0,0 +1,15 @@
1+import { defineConfig } from 'vite';
2+import { resolve } from 'node:path';
3+
4+export default defineConfig({
5+ base: './',
6+ build: {
7+ target: 'es2022',
8+ rollupOptions: {
9+ input: {
10+ main: resolve(__dirname, 'index.html'),
11+ test: resolve(__dirname, 'test.html'),
12+ },
13+ },
14+ },
15+});