Add the cymbal problem and GIF export
A listener near a rigid cymbal in free space, solved with the exterior mode of
the impedance equation. Surfaces without an analytic chart use the Legendre
interpolant of the node data, with short displacements formed from divided
differences of the Legendre polynomials so that they keep full relative
precision in single precision. The page gets a problem selector, box slices
with the cymbal outline, levels relative to the free field, and an Export GIF
button that writes the three planes as one GIF. The GPU test suite checks the
cymbal too, and test.html?bench= compares with the CUDA version's cases.
16 changed files+1803−247
README.mdmodified+22−5View file
@@ -4,7 +4,9 @@ A limitation of the [CUDA proof of concept](https://vault1.magland.org/proof-of-
44
55 **Live page:** https://fmm3dbie-webgpu-poc--proof-of-concept.vault1-sites.magland.org/
66
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.
7+The page solves either of the CUDA version's two showcase problems. The first 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+The second is a listener near a stationary, rigid cymbal in free space. The source is placed at the listener, so that by reciprocity the field u(x) is the pressure at the listener due to a unit source at x. The cymbal is 40 cm across, with a raised bell; since fmm3dbie's equations need a closed surface, it is a thin closed body, 20 mm thick at the center and 10 mm (1/40 of the diameter) at its rounded rim. The equation is the exterior mode of the same impedance equation, with λ = 0. The page lets the user set the frequency, the patch size and order, the rim thickness and the listener's position, and shows the field on three planes through a 1 m box around the cymbal, with the level relative to the free field. `?problem=cymbal` opens the page on this problem.
810
911 ## Method
1012
@@ -12,7 +14,7 @@ The equation, representation and discretization are those of the CUDA version: t
1214
1315 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:
1416
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²).
17+1. For a target on the wall, the vector from the target to a nearby quadrature node is formed in the coordinates of the patch's 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²). For the room this is the superellipsoid's analytic parametrization. For the cymbal, and for any other surface, it is the Legendre interpolant of the node data (the geometry fmm3dbie uses), and the vector comes from divided differences of the Legendre polynomials, P_n(u) − P_n(u_t) = (u − u_t) δP_n, with δP_n from its own three-term recurrence. For a target on a neighboring patch, the CPU finds its preimage in the patch's (extended) interpolant by Gauss–Newton, and passes the small remaining gap along with it.
1618 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.
1719 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.
1820
@@ -20,6 +22,18 @@ Unlike the CUDA version, which uses double precision for the near quadrature, ev
2022
2123 `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.
2224
25+The suite then runs the same checks on a coarse cymbal (192 patches of order 5, with a 20 mm rim). On an RTX 6000 Ada the near weights agree with the double precision reference to 1.7e-5 or better, and the exterior self-test error, with sources spread inside the cymbal, is 2.0e-4 around the cymbal and 6.9e-4 on it, which is the discretization error of that coarse surface.
26+
27+`test.html?bench=small|medium|large` runs the named cymbal cases of the CUDA version, on the same surface and unknowns. On the same RTX 6000 Ada (Chrome under Xvfb, Vulkan), with GMRES tolerance 1e-6:
28+
29+| case | points | WebGPU setup + solve | CUDA setup + solve | self-test error, WebGPU / CUDA | reciprocity, WebGPU / CUDA | scattered field at the listener, WebGPU vs CUDA |
30+|---|---|---|---|---|---|---|
31+| small (1 kHz) | 19,012 | 2.6 + 2.5 s | 0.8 + 1.0 s | 1.5e-5 / 1.5e-5 | 7.9e-5 / 7.7e-5 | 1.2e-6 |
32+| medium (3 kHz) | 35,964 | 3.4 + 5.7 s | 1.3 + 1.9 s | 1.2e-6 / 1.1e-6 | 6.0e-6 / 5.2e-6 | 3.3e-7 |
33+| large (8 kHz) | 42,444 | 3.5 + 5.0 s | 1.5 + 1.8 s | 3.2e-6 / 4.6e-6 | 1.6e-6 / 7.0e-7 | 1.2e-6 |
34+
35+The two codes agree to about 1e-6, so on these cases single precision costs no accuracy that the discretization would otherwise give. The WebGPU code is two to three times slower than CUDA on the same GPU. fmm3dbie on the workstation's 32 CPU threads takes 27 to 81 s for these solves.
36+
2337 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.
2438
2539 ## Running it
@@ -29,15 +43,18 @@ For sharper rooms the discretization, not the precision, limits the accuracy. At
2943 npm run test:node # CPU checks of the surface and the reference quadrature
3044 npm run test:gpu # the GPU validation suite in headless Chrome
3145 npm run smoke # build, then solve the quick case in the app, headless
46+ node scripts/smoke.mjs quick cymbal # the same for the cymbal, after a build
47+
48+The page accepts `?problem=cymbal`, `?case=<name>` and `?autosolve`, which solves at once. After a solve, Export GIF writes the three planes side by side as one GIF, animated over one period for the pressure (24 frames) or a single image for the level.
3249
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.
50+The headless tests use a hardware WebGPU adapter when Chrome finds one and SwiftShader otherwise. On a Linux machine without a display, headless Chrome may only find SwiftShader even when there is a GPU; a windowed Chrome under Xvfb (`xvfb-run -a google-chrome --ignore-gpu-blocklist --enable-unsafe-webgpu --enable-features=Vulkan ...`) found the NVIDIA adapter in our tests. 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.
3451
3552 ## Limitations
3653
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.
54+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 room's medium case (17,496 points) and more than most browsers allow in one buffer for its large case (48,600 points). The thin cymbal has more near pairs per node, and its small, medium and large cases need 640, 730 and 780 MB. Splitting the weights across buffers, or recomputing them, would remove this limit.
3855
3956 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.
4057
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.
58+The interpolated geometry with divided differences works for any surface of quadrilateral patches, but only the rounded cube and the cymbal are wired into the page. The far-field energy check for the cymbal runs on the CPU and is skipped for the largest cases.
4259
4360 Timings on the software adapter are not representative, and we have not yet measured the named cases on hardware GPUs.
index.htmlmodified+26−20View file
@@ -4,7 +4,7 @@
44 <meta charset="utf-8" />
55 <meta name="viewport" content="width=device-width, initial-scale=1" />
66 <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>
7+ <title>fmm3dbie-webgpu-poc: a room with impedance walls, or a listener near a cymbal, solved on your GPU</title>
88 <style>
99 :root {
1010 --bg: #ffffff;
@@ -86,37 +86,41 @@
8686 .notes p { margin: 6px 0; }
8787 .notes code { font-size: 12.5px; }
8888 #device { font-size: 12px; color: var(--ink-2); margin-top: 18px; }
89+ .hidden { display: none !important; }
8990 </style>
9091 </head>
9192 <body>
9293 <main>
9394 <h1>fmm3dbie-webgpu-poc</h1>
9495 <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.
96+ Acoustic fields at a fixed frequency, solved with a high-order boundary integral equation on your GPU through WebGPU, in single precision: a point source in a closed room with partially absorbing walls, or a listener near a rigid cymbal in free space. 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
97+ <a href="https://github.com/fastalgorithms/fmm3dbie">fmm3dbie</a>'s Helmholtz impedance solver. Choose the problem, press Solve, then turn the view with the mouse and click a slice to move the source.
9798 </p>
9899
99100 <div class="panel">
100101 <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>
102+ <label class="field" title="The room: a point source in a closed room with impedance walls. The cymbal: a listener near a rigid cymbal in free space, solved with the source at the listener (reciprocity)."><span>problem</span>
103+ <select autocomplete="off" id="problem">
104+ <option value="room" selected>room</option>
105+ <option value="cymbal">listener near a cymbal</option>
108106 </select>
109107 </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>
108+ <label class="field" title="Named cases from the CUDA proof of concept, plus a quick one for slow GPUs"><span>case</span>
109+ <select autocomplete="off" id="preset"></select>
110+ </label>
111+ <label class="field room-only" 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>
112+ <label class="field room-only" 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>
113+ <label class="field room-only" 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>
114+ <label class="field cymbal-only" title="Frequency in Hz; the speed of sound is 343 m/s and the cymbal is 40 cm across"><span>frequency (Hz)</span><input autocomplete="off" id="freq" type="number" min="100" max="16000" step="100" /></label>
115+ <label class="field cymbal-only" title="Target patch size in cm; the mesh is also refined toward the rim and across the bell, whatever this is"><span>patch size (cm)</span><input autocomplete="off" id="h" type="number" min="2" max="20" step="0.5" /></label>
116+ <label class="field cymbal-only" title="Thickness at the rim in mm; the cymbal is twice as thick at the center. Real cymbals are about 1 mm thick; thinner needs more patches."><span>rim thickness (mm)</span><input autocomplete="off" id="tRim" type="number" min="4" max="40" step="1" /></label>
113117 <label class="field" title="Polynomial order per patch; (order + 1)² nodes each"><span>order per patch</span>
114118 <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>
115119 </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>
120+ <label class="field room-only" 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>
117121 <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>
118122 </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>
123+ <label class="field" id="x0label" title="Source position; must be inside the room. You can also click a slice."><span id="x0name">source (x, y, z)</span>
120124 <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>
121125 </label>
122126 <label class="field" title="Relative residual at which GMRES stops"><span>GMRES tolerance</span>
@@ -129,7 +133,7 @@
129133 <div class="row">
130134 <button id="solve" class="primary">Solve</button>
131135 <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>
136+ <label title="Also solve a problem with a known solution and report the error: the field of a source outside the room, or of sources spread inside the cymbal"><input autocomplete="off" id="selftest" type="checkbox" /> self-test</label>
133137 <span class="derived" id="derived"></span>
134138 </div>
135139 <div class="row" style="margin-top: 6px"><progress id="progress" max="1" value="0"></progress><span id="status">loading…</span></div>
@@ -140,14 +144,15 @@
140144 <label>show
141145 <select autocomplete="off" id="mode">
142146 <option value="pressure">pressure Re(u e^−iωt)</option>
143- <option value="level">level (dB re median)</option>
147+ <option value="level" id="levelOption">level (dB re median)</option>
144148 </select>
145149 </label>
146150 <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>
151+ <label><input autocomplete="off" id="showWalls" type="checkbox" checked /> <span id="wallsName">walls</span></label>
148152 <label><input autocomplete="off" id="showXY" type="checkbox" checked /> xy</label>
149153 <label><input autocomplete="off" id="showXZ" type="checkbox" checked /> xz</label>
150154 <label><input autocomplete="off" id="showYZ" type="checkbox" /> yz</label>
155+ <button id="exportGif" disabled title="The three planes side by side as one GIF: animated over one period for the pressure, a single image for the level">Export GIF</button>
151156 </div>
152157
153158 <div class="stage">
@@ -155,7 +160,7 @@
155160 <canvas id="room"></canvas>
156161 <div class="empty" id="empty">Press Solve to compute the field.</div>
157162 <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>
163+ <div class="hint" id="hint">drag to turn · scroll to zoom · the nearer walls are cut away</div>
159164 </div>
160165 <div class="side">
161166 <div class="panel">
@@ -171,7 +176,8 @@
171176 <h2>What is computed</h2>
172177 <div class="notes">
173178 <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>
179+ <p>The cymbal is 40 cm across, stationary and sound-hard, in free space, with a listener nearby. We put the source at the listener: by reciprocity the field u(x) is then the pressure at the listener due to a unit source at x, so one solve gives the listener's response to a source anywhere. The planes cut a 1 m box around the cymbal, which only selects where the field is shown; the integral equation needs no truncation. Since fmm3dbie's equations need a closed surface, the cymbal is a thin closed body with a raised bell, 20 mm thick at the center and thinner toward the rounded rim. The level is shown relative to the free field, so that 0 dB means no effect of the cymbal.</p>
180+ <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, or, for the cymbal, with the power radiated to infinity. Large cases may need more GPU memory than the browser allows.</p>
175181 </div>
176182 <div id="device"></div>
177183 </main>
scripts/smoke.mjsmodified+6−5View file
@@ -3,9 +3,9 @@
33 * Chrome (hardware WebGPU if there is any, SwiftShader otherwise), solves
44 * the quick case, and waits for the status line to report
55 * success. Any console error, page error or failed request fails the run. A
6- * screenshot goes to tmp/smoke.png.
6+ * screenshot goes to tmp/smoke-<problem>.png.
77 *
8- * Run after `vite build`: node scripts/smoke.mjs [case] (default quick)
8+ * Run after `vite build`: node scripts/smoke.mjs [case] [problem] (default quick room; problem room or cymbal)
99 */
1010 import { createServer } from 'node:http';
1111 import { mkdir, readFile } from 'node:fs/promises';
@@ -16,6 +16,7 @@ const DIST = new URL('../dist/', import.meta.url).pathname;
1616 const TMP = new URL('../tmp/', import.meta.url).pathname;
1717 const CHROME = process.env.CHROME_PATH ?? '/usr/bin/google-chrome';
1818 const CASE = process.argv[2] ?? 'quick';
19+const PROBLEM = process.argv[3] ?? 'room';
1920 const MIME = { '.html': 'text/html', '.js': 'text/javascript', '.css': 'text/css', '.json': 'application/json', '.svg': 'image/svg+xml' };
2021
2122 const server = createServer(async (req, res) => {
@@ -52,7 +53,7 @@ for (const flags of flagSets) {
5253 });
5354 page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`));
5455 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.goto(`http://127.0.0.1:${port}/index.html?noroom&problem=${PROBLEM}`, { waitUntil: 'load' });
5657 await page.waitForFunction(() => /^ready/.test(document.getElementById('status')?.textContent ?? ''), { timeout: 60_000, polling: 250 });
5758 console.log(await page.$eval('#device', (e) => e.textContent));
5859 await page.select('#preset', CASE);
@@ -76,8 +77,8 @@ for (const flags of flagSets) {
7677 for (const r of rows) console.log(` ${r}`);
7778 if (err) problems.push(`error line: ${err}`);
7879 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')}`);
80+ await page.screenshot({ path: join(TMP, `smoke-${PROBLEM}.png`), fullPage: true });
81+ console.log(`screenshot: ${join(TMP, `smoke-${PROBLEM}.png`)}`);
8182 if (!/^done/.test(status)) problems.push(`status: ${status}`);
8283 ok = problems.length === 0;
8384 } catch (e) {
src/cymbal.tsadded+395−0View file
@@ -0,0 +1,395 @@
1+/**
2+ * A cymbal as a thin closed body built from quadrilateral patches. A port of
3+ * fmm3dbie_gpu/cymbal.py from the CUDA proof of concept; lengths in meters.
4+ *
5+ * fmm3dbie's integral equations need a closed surface, so the cymbal is not a
6+ * zero-thickness shell but a thin body around an axisymmetric midsurface
7+ * z = F(rho), rho = x^2 + y^2, with thickness T(rho). The midsurface has a raised
8+ * bell at the center and a shallow bow outside it; the thickness decreases from
9+ * the center to the rim, with zero slope at the rim. The top and bottom faces are
10+ * the offset surfaces M +- (T/2) N, with M the midsurface and N its upward unit
11+ * normal, and a half-torus of diameter T closes the body at the rim. Every patch
12+ * lies entirely on one face or on the rim, so every patch is analytic.
13+ *
14+ * The faces are meshed with an O-grid in the (x, y) plane: an equiangular central
15+ * square of m_c x m_c patches, then rings of patches, graded toward the rim and
16+ * refined across the bell. Neighboring rings need not have the same number of
17+ * patches around (the mesh is not conforming). The node derivatives dx/du, dx/dv
18+ * are computed with forward-mode dual numbers, exactly to rounding. The surface
19+ * carries no analytic chart for the solver: the integrals use the Legendre
20+ * interpolant of the node data, as in fmm3dbie.
21+ */
22+import { discNodes, Surface } from './surface.ts';
23+
24+// ------------------------------------------------------------ dual numbers
25+
26+/** a + b eps with eps^2 = 0: a value and its derivative along one direction. */
27+class D {
28+ readonly v: number;
29+ readonly d: number;
30+ constructor(v: number, d: number) {
31+ this.v = v;
32+ this.d = d;
33+ }
34+}
35+type DN = D | number;
36+const lift = (a: DN): D => (typeof a === 'number' ? new D(a, 0) : a);
37+const add = (a: DN, b: DN) => {
38+ const x = lift(a);
39+ const y = lift(b);
40+ return new D(x.v + y.v, x.d + y.d);
41+};
42+const sub = (a: DN, b: DN) => {
43+ const x = lift(a);
44+ const y = lift(b);
45+ return new D(x.v - y.v, x.d - y.d);
46+};
47+const mul = (a: DN, b: DN) => {
48+ const x = lift(a);
49+ const y = lift(b);
50+ return new D(x.v * y.v, x.d * y.v + x.v * y.d);
51+};
52+const div = (a: DN, b: DN) => {
53+ const x = lift(a);
54+ const y = lift(b);
55+ return new D(x.v / y.v, (x.d * y.v - x.v * y.d) / (y.v * y.v));
56+};
57+const dsqrt = (a: D) => {
58+ const s = Math.sqrt(a.v);
59+ return new D(s, a.d / (2 * s));
60+};
61+const dexp = (a: D) => {
62+ const e = Math.exp(a.v);
63+ return new D(e, e * a.d);
64+};
65+const dcos = (a: D) => new D(Math.cos(a.v), -Math.sin(a.v) * a.d);
66+const dsin = (a: D) => new D(Math.sin(a.v), Math.cos(a.v) * a.d);
67+const dtan = (a: D) => {
68+ const t = Math.tan(a.v);
69+ return new D(t, (1 + t * t) * a.d);
70+};
71+
72+// ------------------------------------------------------------------ shape
73+
74+export interface CymbalParams {
75+ /** outer radius, including the rim */
76+ radius: number;
77+ /** thickness at the rim */
78+ tRim: number;
79+ /** thickness at the center / tRim */
80+ taper: number;
81+ bellRadius: number;
82+ bellHeight: number;
83+ bellWidth: number;
84+ /** fall of the midsurface from the center to the rim, apart from the bell */
85+ bowHeight: number;
86+}
87+
88+export const CYMBAL_DEFAULTS: CymbalParams = {
89+ radius: 0.2,
90+ tRim: 0.01,
91+ taper: 2.0,
92+ bellRadius: 0.06,
93+ bellHeight: 0.025,
94+ bellWidth: 0.01,
95+ bowHeight: 0.015,
96+};
97+
98+export class CymbalShape {
99+ readonly par: CymbalParams;
100+ readonly rMid: number;
101+
102+ constructor(par: CymbalParams = CYMBAL_DEFAULTS) {
103+ this.par = par;
104+ this.rMid = par.radius - par.tRim / 2;
105+ }
106+
107+ private bellArg(rho: DN): D {
108+ const { bellRadius: rb, bellWidth: w } = this.par;
109+ return div(sub(rho, rb * rb), 2 * rb * w);
110+ }
111+
112+ /** midsurface height as a function of rho */
113+ F(rho: DN): D {
114+ const { bowHeight, bellHeight } = this.par;
115+ const R2 = this.rMid ** 2;
116+ return add(mul(bowHeight, sub(1, div(rho, R2))), div(bellHeight, add(1, dexp(this.bellArg(rho)))));
117+ }
118+
119+ /** dF/drho */
120+ dF(rho: DN): D {
121+ const { bowHeight, bellHeight, bellRadius: rb, bellWidth: w } = this.par;
122+ const e = dexp(this.bellArg(rho));
123+ const ope = add(1, e);
124+ return sub(-bowHeight / this.rMid ** 2, div(mul(bellHeight / (2 * rb * w), e), mul(ope, ope)));
125+ }
126+
127+ /** thickness as a function of rho */
128+ T(rho: DN): D {
129+ const { tRim, taper } = this.par;
130+ const s = sub(1, div(rho, this.rMid ** 2));
131+ return mul(tRim, add(1, mul(taper - 1, mul(s, s))));
132+ }
133+
134+ /** point of the top (side = 1) or bottom (side = -1) face over (x, y) of the midsurface */
135+ facePoint(x: DN, y: DN, side: number): [D, D, D] {
136+ const rho = add(mul(x, x), mul(y, y));
137+ const g = mul(2, this.dF(rho));
138+ const q = dsqrt(add(1, mul(mul(g, g), rho)));
139+ const s = div(mul(side / 2, this.T(rho)), q);
140+ const sg = mul(s, g);
141+ return [sub(x, mul(sg, x)), sub(y, mul(sg, y)), add(this.F(rho), s)];
142+ }
143+
144+ /** point of the rim at azimuth theta; phi from -pi/2 (bottom face) to pi/2 (top face) */
145+ rimPoint(theta: DN, phi: DN): [D, D, D] {
146+ const R = this.rMid;
147+ const rho = R * R;
148+ const fp = 2 * R * this.dF(rho).v;
149+ const q = Math.sqrt(1 + fp * fp);
150+ const h = this.T(rho).v / 2;
151+ const a = mul(h / q, dcos(lift(phi)));
152+ const b = mul(h / q, dsin(lift(phi)));
153+ const rr = sub(add(R, a), mul(b, fp));
154+ const th = lift(theta);
155+ return [mul(rr, dcos(th)), mul(rr, dsin(th)), add(add(this.F(rho).v, mul(a, fp)), b)];
156+ }
157+
158+ /** closed outline of the body in the (r, z) half plane, counterclockwise: [r[], z[]] */
159+ crossSection(n = 400): [Float64Array, Float64Array] {
160+ const r: number[] = [];
161+ const z: number[] = [];
162+ for (let i = 0; i < n; i++) {
163+ const p = this.facePoint((this.rMid * i) / (n - 1), 0, -1);
164+ r.push(p[0].v);
165+ z.push(p[2].v);
166+ }
167+ const nr = Math.floor(n / 4);
168+ for (let i = 0; i < nr; i++) {
169+ const p = this.rimPoint(0, -Math.PI / 2 + (Math.PI * i) / (nr - 1));
170+ r.push(p[0].v);
171+ z.push(p[2].v);
172+ }
173+ for (let i = n - 1; i >= 0; i--) {
174+ const p = this.facePoint((this.rMid * i) / (n - 1), 0, 1);
175+ r.push(p[0].v);
176+ z.push(p[2].v);
177+ }
178+ return [Float64Array.from(r), Float64Array.from(z)];
179+ }
180+
181+ /** volume by Pappus' theorem on the cross section (a check) */
182+ volume(n = 20000): number {
183+ const [r, z] = this.crossSection(n);
184+ let s = 0;
185+ for (let i = 0; i < r.length; i++) {
186+ const j = (i + 1) % r.length;
187+ s += ((r[i] * r[i] + r[i] * r[j] + r[j] * r[j]) / 6) * (z[j] - z[i]);
188+ }
189+ return 2 * Math.PI * s;
190+ }
191+}
192+
193+/**
194+ * A test for points inside the body, or within margin of its surface. Every
195+ * vertical line meets the cross section in one interval, so the outline, offset
196+ * outward by margin, splits at its outermost point into a lower and an upper
197+ * chain, and a point is inside when z lies between them at its radius.
198+ */
199+export function insideTest(shape: CymbalShape, margin = 0): (x: ArrayLike<number>) => boolean {
200+ const [r0, z0] = shape.crossSection(4000);
201+ const n = r0.length;
202+ const r = new Float64Array(n);
203+ const z = new Float64Array(n);
204+ for (let i = 0; i < n; i++) {
205+ const a = (i + n - 1) % n;
206+ const b = (i + 1) % n;
207+ let tr = r0[b] - r0[a];
208+ let tz = z0[b] - z0[a];
209+ const l = Math.hypot(tr, tz);
210+ tr /= l;
211+ tz /= l;
212+ r[i] = Math.max(0, r0[i] + margin * tz);
213+ z[i] = z0[i] - margin * tr;
214+ }
215+ let im = 0;
216+ for (let i = 1; i < n; i++) if (r[i] > r[im]) im = i;
217+ const chain = (idx: number[]) => {
218+ const cr: number[] = [];
219+ const cz: number[] = [];
220+ let m = -Infinity;
221+ for (const i of idx) {
222+ m = Math.max(m, r[i]);
223+ cr.push(m);
224+ cz.push(z[i]);
225+ }
226+ return { cr, cz };
227+ };
228+ const lo = chain(Array.from({ length: im + 1 }, (_, i) => i));
229+ const up = chain(Array.from({ length: n - im }, (_, i) => n - 1 - i));
230+ const interp = (c: { cr: number[]; cz: number[] }, x: number) => {
231+ const { cr, cz } = c;
232+ if (x <= cr[0]) return cz[0];
233+ let a = 0;
234+ let b = cr.length - 1;
235+ if (x >= cr[b]) return cz[b];
236+ while (b - a > 1) {
237+ const m = (a + b) >> 1;
238+ if (cr[m] <= x) a = m;
239+ else b = m;
240+ }
241+ const t = cr[b] > cr[a] ? (x - cr[a]) / (cr[b] - cr[a]) : 0;
242+ return cz[a] + t * (cz[b] - cz[a]);
243+ };
244+ const rmax = r[im];
245+ return (x) => {
246+ const rr = Math.hypot(x[0], x[1]);
247+ return rr <= rmax && x[2] > interp(lo, rr) && x[2] < interp(up, rr);
248+ };
249+}
250+
251+// ----------------------------------------------------------------- the mesh
252+
253+export interface RingLayout {
254+ mC: number;
255+ c: number;
256+ radii: number[];
257+ counts: number[];
258+ nrim: number;
259+}
260+
261+/**
262+ * Radii and patch counts of the O-grid for a target patch size h: the central
263+ * square has half-width c and mC x mC patches (mC >= 2); ring i lies between
264+ * radii[i] and radii[i + 1] (radii[0] = c means the square) and has counts[i]
265+ * patches around. The ring around the square reaches r = 2c, so that it is 0.59c
266+ * wide at the square's corners.
267+ */
268+export function ringLayout(shape: CymbalShape, h: number, grade = 1.5, cFrac = 0.12): RingLayout {
269+ const R = shape.rMid;
270+ const { tRim, bellRadius, bellWidth } = shape.par;
271+ const c = cFrac * R;
272+ const mC = Math.max(2, Math.ceil((2 * c) / h));
273+ const first = 2 * c;
274+ const b0 = bellRadius - 2.5 * bellWidth;
275+ const b1 = bellRadius + 2.5 * bellWidth;
276+ const wmax = (edge: number) => {
277+ let w = Math.min(h, tRim + (grade - 1) * (R - edge));
278+ if (edge > b0 && edge - w < b1) w = Math.min(w, 2 * bellWidth);
279+ return w;
280+ };
281+ const out = [R];
282+ for (;;) {
283+ const last = out[out.length - 1];
284+ const w = wmax(last);
285+ if (last - w < first + 0.5 * w) break;
286+ out.push(last - w);
287+ }
288+ const radii = [c, first, ...out.reverse()];
289+ const counts: number[] = [];
290+ for (let i = 0; i < radii.length - 1; i++) {
291+ const width = radii[i + 1] - radii[i];
292+ const arc = Math.min(h, 2.5 * width);
293+ counts.push(4 * Math.max(mC, Math.ceil((2 * Math.PI * radii[i + 1]) / (4 * arc))));
294+ }
295+ counts[0] = 4 * mC;
296+ for (let i = 1; i < counts.length; i++) counts[i] = Math.max(counts[i], counts[i - 1]);
297+ return { mC, c, radii, counts, nrim: counts[counts.length - 1] };
298+}
299+
300+type Chart = (u: D, v: D) => [D, D];
301+type PatchMap = (u: D, v: D) => [D, D, D];
302+
303+/** planar charts (u, v) in [-1, 1]^2 -> (x, y) of the face patches, oriented so that d/du x d/dv points up */
304+function faceCharts(L: RingLayout): Chart[] {
305+ const charts: Chart[] = [];
306+ const q = Math.PI / 4;
307+ const { mC, c, radii, counts } = L;
308+ for (let i = 0; i < mC; i++) {
309+ for (let j = 0; j < mC; j++) {
310+ const a0 = -1 + (2 * i) / mC;
311+ const a1 = -1 + (2 * (i + 1)) / mC;
312+ const b0 = -1 + (2 * j) / mC;
313+ const b1 = -1 + (2 * (j + 1)) / mC;
314+ charts.push((u, v) => {
315+ const s = add(a0, mul(0.5 * (a1 - a0), add(u, 1)));
316+ const t = add(b0, mul(0.5 * (b1 - b0), add(v, 1)));
317+ return [mul(c, dtan(mul(q, s))), mul(c, dtan(mul(q, t)))];
318+ });
319+ }
320+ }
321+ for (let ir = 0; ir < radii.length - 1; ir++) {
322+ const r0 = radii[ir];
323+ const r1 = radii[ir + 1];
324+ const n = counts[ir];
325+ const dth = (2 * Math.PI) / n;
326+ for (let k = 0; k < n; k++) {
327+ const th0 = -q + k * dth;
328+ if (ir === 0) {
329+ charts.push((u, v) => {
330+ const th = add(th0, mul(0.5 * dth, add(v, 1)));
331+ const sec = Math.floor((th.v + q) / (2 * q));
332+ const rin = div(c, dcos(sub(th, sec * 2 * q)));
333+ const a = mul(0.5, add(u, 1));
334+ const rad = add(mul(sub(1, a), rin), mul(a, r1));
335+ return [mul(rad, dcos(th)), mul(rad, dsin(th))];
336+ });
337+ } else {
338+ charts.push((u, v) => {
339+ const th = add(th0, mul(0.5 * dth, add(v, 1)));
340+ const rad = add(r0, mul(0.5 * (r1 - r0), add(u, 1)));
341+ return [mul(rad, dcos(th)), mul(rad, dsin(th))];
342+ });
343+ }
344+ }
345+ }
346+ return charts;
347+}
348+
349+/** all patch maps (u, v) -> xyz with outward orientation: top faces, bottom faces, rim */
350+function patchMaps(shape: CymbalShape, h: number, grade: number): { maps: PatchMap[]; layout: RingLayout } {
351+ const layout = ringLayout(shape, h, grade);
352+ const charts = faceCharts(layout);
353+ const maps: PatchMap[] = [];
354+ for (const ch of charts) maps.push((u, v) => shape.facePoint(...ch(u, v), 1));
355+ for (const ch of charts) maps.push((u, v) => shape.facePoint(...ch(v, u), -1)); // swapping u and v flips the normal
356+ const { nrim } = layout;
357+ for (let k = 0; k < nrim; k++) {
358+ const th0 = -Math.PI / 4 + (2 * Math.PI * k) / nrim;
359+ maps.push((u, v) => shape.rimPoint(add(th0, mul(Math.PI / nrim, add(u, 1))), mul(0.5 * Math.PI, v)));
360+ }
361+ return { maps, layout };
362+}
363+
364+/** Node data (12, npts) in the fmm3dbie iptype 11 convention, and the layout. */
365+export function cymbalSrcvals(shape: CymbalShape, h: number, norder: number, grade = 1.5): { srcvals: Float64Array; layout: RingLayout } {
366+ const { uv } = discNodes(norder);
367+ const npols = (norder + 1) ** 2;
368+ const { maps, layout } = patchMaps(shape, h, grade);
369+ const npts = maps.length * npols;
370+ const sv = new Float64Array(12 * npts);
371+ for (const [j, f] of maps.entries()) {
372+ for (let k = 0; k < npols; k++) {
373+ const i = j * npols + k;
374+ const pu = f(new D(uv[k], 1), new D(uv[npols + k], 0));
375+ const pv = f(new D(uv[k], 0), new D(uv[npols + k], 1));
376+ const ru = [pu[0].d, pu[1].d, pu[2].d];
377+ const rv = [pv[0].d, pv[1].d, pv[2].d];
378+ const nx = ru[1] * rv[2] - ru[2] * rv[1];
379+ const ny = ru[2] * rv[0] - ru[0] * rv[2];
380+ const nz = ru[0] * rv[1] - ru[1] * rv[0];
381+ const l = Math.hypot(nx, ny, nz);
382+ const vals = [pu[0].v, pu[1].v, pu[2].v, ...ru, ...rv, nx / l, ny / l, nz / l];
383+ for (let c = 0; c < 12; c++) sv[c * npts + i] = vals[c];
384+ }
385+ }
386+ return { srcvals: sv, layout };
387+}
388+
389+/** The cymbal as a Surface with the interpolated geometry, for a target patch size h (m). */
390+export function cymbal(shape: CymbalShape, h: number, norder: number, grade = 1.5): Surface & { layout: RingLayout } {
391+ const { srcvals, layout } = cymbalSrcvals(shape, h, norder, grade);
392+ const s = new Surface(norder, srcvals, null) as Surface & { layout: RingLayout };
393+ s.layout = layout;
394+ return s;
395+}
src/cymbalProblem.tsadded+202−0View file
@@ -0,0 +1,202 @@
1+/**
2+ * The cymbal problem: a reciprocity solve for a listener near a rigid cymbal
3+ * (problems/cymbal.py in the CUDA proof of concept). A point source at the
4+ * listener's position L radiates in free space, and the only obstacle is a
5+ * stationary, sound-hard cymbal:
6+ *
7+ * (Delta + k^2) u = -delta(x - L) outside the cymbal,
8+ * du/dn = 0 on the cymbal,
9+ *
10+ * with the radiation condition. By reciprocity u(x) is also the pressure at the
11+ * listener due to a unit source at x. We solve the exterior mode of the same
12+ * impedance equation with lambda = 0, and show the field on three planes through
13+ * a 1 m box around the cymbal. Lengths in meters, sound speed 343 m/s.
14+ */
15+import { CymbalShape, insideTest } from './cymbal.ts';
16+import type { CVec, Slice } from './problem.ts';
17+
18+export const C_SOUND = 343;
19+export const RADIUS = 0.2;
20+export const LISTENER: [number, number, number] = [0.25, 0, 0.45];
21+/** the evaluation box: [lo, hi] in x, y, z */
22+export const BOX: [number, number][] = [
23+ [-0.5, 0.5],
24+ [-0.5, 0.5],
25+ [-0.35, 0.65],
26+];
27+/** height of the horizontal slice, below the cymbal */
28+export const SLICE_BELOW = -0.1;
29+
30+/**
31+ * Point sources inside the cymbal whose field is the exact solution outside, for
32+ * the self-test: positions (3, n) and strengths. A single source would be at most
33+ * T/2 = 1 cm from the surface, and its data would vary on that scale, which the
34+ * main problem's data does not. Instead the sources are spread 3 mm apart over
35+ * the midsurface out to r = rMax, with strengths tapering smoothly to zero.
36+ */
37+export function selftestSources(shape: CymbalShape, rMax = 0.12, spacing = 0.003): { x0s: Float64Array; ws: Float64Array } {
38+ const g: number[] = [];
39+ for (let t = -rMax; t <= rMax + spacing / 2; t += spacing) g.push(t);
40+ const pts: number[][] = [];
41+ const ws: number[] = [];
42+ for (const x of g) {
43+ for (const y of g) {
44+ const rho = x * x + y * y;
45+ if (rho >= rMax * rMax) continue;
46+ pts.push([x, y, shape.F(rho).v]);
47+ ws.push(((1 - rho / rMax ** 2) ** 3 * spacing ** 2) / ((Math.PI * rMax ** 2) / 4));
48+ }
49+ }
50+ const n = pts.length;
51+ const x0s = new Float64Array(3 * n);
52+ for (let j = 0; j < n; j++) for (let d = 0; d < 3; d++) x0s[d * n + j] = pts[j][d];
53+ return { x0s, ws: Float64Array.from(ws) };
54+}
55+
56+/** sum_j ws[j] G(x, x0s_j) at points x (3, n), or with nrm given its normal derivative. */
57+export function sourceSum(x: Float64Array, x0s: Float64Array, ws: Float64Array, k: number, nrm?: Float64Array): CVec {
58+ const n = x.length / 3;
59+ const m = ws.length;
60+ const out = new Float64Array(2 * n);
61+ for (let i = 0; i < n; i++) {
62+ let sr = 0;
63+ let si = 0;
64+ for (let j = 0; j < m; j++) {
65+ const dx = x[i] - x0s[j];
66+ const dy = x[n + i] - x0s[m + j];
67+ const dz = x[2 * n + i] - x0s[2 * m + j];
68+ const r = Math.hypot(dx, dy, dz);
69+ const s = ws[j] / (4 * Math.PI * r);
70+ let gr = s * Math.cos(k * r);
71+ let gi = s * Math.sin(k * r);
72+ if (nrm) {
73+ const dn = (dx * nrm[i] + dy * nrm[n + i] + dz * nrm[2 * n + i]) / r;
74+ // g (ik - 1/r) dn
75+ const tr = (-gr / r - gi * k) * dn;
76+ const ti = (gr * k - gi / r) * dn;
77+ gr = tr;
78+ gi = ti;
79+ }
80+ sr += gr;
81+ si += gi;
82+ }
83+ out[2 * i] = sr;
84+ out[2 * i + 1] = si;
85+ }
86+ return out;
87+}
88+
89+/** n points around the cymbal, at least 5 mm from its surface (deterministic). */
90+export function selftestPoints(shape: CymbalShape, n: number): Float64Array {
91+ const near = insideTest(shape, 0.005);
92+ const pts: number[][] = [];
93+ let s = 12345;
94+ const rnd = () => {
95+ s = (Math.imul(s, 1664525) + 1013904223) >>> 0;
96+ return s / 4294967296;
97+ };
98+ while (pts.length < n) {
99+ const x = [-0.3 + 0.6 * rnd(), -0.3 + 0.6 * rnd(), -0.15 + 0.35 * rnd()];
100+ if (!near(x)) pts.push(x);
101+ }
102+ const out = new Float64Array(3 * n);
103+ for (let j = 0; j < n; j++) for (let d = 0; d < 3; d++) out[d * n + j] = pts[j][d];
104+ return out;
105+}
106+
107+/**
108+ * Power radiated to infinity, k int |u_inf|^2 dOmega, for u = G(., x0) + S_k sigma +
109+ * i alpha D_k phi1, with u ~ u_inf exp(ikr)/r. The far field pattern is computed
110+ * with the smooth rule at the nodes, and the sphere integral with Gauss-Legendre
111+ * in cos(theta) and the trapezoid rule in phi, exact for the band limit. For a
112+ * rigid cymbal it equals the power emitted by the source, k/(4 pi) + Im v(x0).
113+ */
114+export function farFieldPower(
115+ k: number,
116+ x0: ArrayLike<number>,
117+ xs: Float64Array,
118+ nrm: Float64Array,
119+ wts: Float64Array,
120+ sigma: CVec,
121+ phi1: CVec,
122+ alpha: number,
123+ leggauss: (n: number) => { x: Float64Array; w: Float64Array },
124+): number {
125+ const n = wts.length;
126+ let rmax = Math.hypot(x0[0], x0[1], x0[2]);
127+ for (let i = 0; i < n; i++) rmax = Math.max(rmax, Math.hypot(xs[i], xs[n + i], xs[2 * n + i]));
128+ const L = Math.ceil(k * rmax) + 20;
129+ const { x: ct, w: wct } = leggauss(L + 1);
130+ const nphi = 2 * L + 2;
131+ // qs = w sigma, qd = i alpha w phi1
132+ const qs = new Float64Array(2 * n);
133+ const qd = new Float64Array(2 * n);
134+ for (let i = 0; i < n; i++) {
135+ qs[2 * i] = wts[i] * sigma[2 * i];
136+ qs[2 * i + 1] = wts[i] * sigma[2 * i + 1];
137+ qd[2 * i] = -alpha * wts[i] * phi1[2 * i + 1];
138+ qd[2 * i + 1] = alpha * wts[i] * phi1[2 * i];
139+ }
140+ let total = 0;
141+ for (let a = 0; a <= L; a++) {
142+ const st = Math.sqrt(1 - ct[a] ** 2);
143+ for (let b = 0; b < nphi; b++) {
144+ const ph = (2 * Math.PI * b) / nphi;
145+ const d = [st * Math.cos(ph), st * Math.sin(ph), ct[a]];
146+ const p0 = -k * (d[0] * x0[0] + d[1] * x0[1] + d[2] * x0[2]);
147+ let ur = Math.cos(p0);
148+ let ui = Math.sin(p0);
149+ for (let i = 0; i < n; i++) {
150+ const ph1 = -k * (d[0] * xs[i] + d[1] * xs[n + i] + d[2] * xs[2 * n + i]);
151+ const er = Math.cos(ph1);
152+ const ei = Math.sin(ph1);
153+ const dn = d[0] * nrm[i] + d[1] * nrm[n + i] + d[2] * nrm[2 * n + i];
154+ // e (qs + (-i k dn) qd)
155+ const cr = qs[2 * i] + k * dn * qd[2 * i + 1];
156+ const ci = qs[2 * i + 1] - k * dn * qd[2 * i];
157+ ur += er * cr - ei * ci;
158+ ui += er * ci + ei * cr;
159+ }
160+ total += ((wct[a] * 2 * Math.PI) / nphi) * (ur * ur + ui * ui);
161+ }
162+ }
163+ return (k * total) / (16 * Math.PI * Math.PI);
164+}
165+
166+/**
167+ * Three planes through the box: xz through the listener, yz through the cymbal's
168+ * axis, and xy below the cymbal. The mask is 1 outside the cymbal and away from
169+ * the source.
170+ */
171+export function boxSlices(ngrid: number, x0: ArrayLike<number>, shape: CymbalShape): Slice[] {
172+ const inside = insideTest(shape, 1e-3);
173+ const defs: [Slice['name'], [number, number], number, number][] = [
174+ ['xz', [0, 2], 1, x0[1]],
175+ ['yz', [1, 2], 0, 0],
176+ ['xy', [0, 1], 2, SLICE_BELOW],
177+ ];
178+ const out: Slice[] = [];
179+ for (const [name, axes, fixed, value] of defs) {
180+ const N = ngrid * ngrid;
181+ const pts = new Float64Array(3 * N);
182+ const mask = new Uint8Array(N);
183+ const [ea, eb] = [BOX[axes[0]], BOX[axes[1]]];
184+ const h = (ea[1] - ea[0]) / (ngrid - 1);
185+ for (let row = 0; row < ngrid; row++) {
186+ for (let col = 0; col < ngrid; col++) {
187+ const q = row * ngrid + col;
188+ const x = [0, 0, 0];
189+ x[axes[0]] = ea[0] + ((ea[1] - ea[0]) * col) / (ngrid - 1);
190+ x[axes[1]] = eb[0] + ((eb[1] - eb[0]) * row) / (ngrid - 1);
191+ x[fixed] = value;
192+ pts[q] = x[0];
193+ pts[N + q] = x[1];
194+ pts[2 * N + q] = x[2];
195+ const atSource = Math.hypot(x[0] - x0[0], x[1] - x0[1], x[2] - x0[2]) < 1e-3 * h;
196+ mask[q] = inside(x) || atSource ? 0 : 1;
197+ }
198+ }
199+ out.push({ name, axes, fixed, value, extent: [[ea[0], ea[1]], [eb[0], eb[1]]], pts, mask });
200+ }
201+ return out;
202+}
src/gpu/operator.tsmodified+126−23View file
@@ -49,7 +49,12 @@ export interface NearPairs {
4949 pairT: Int32Array;
5050 pairP: Int32Array;
5151 rowPtr: Uint32Array;
52- /** per pair: (u_t, v_t, kind, 0); kind 0 off-surface target, 1 on-surface, 2 self patch */
52+ /**
53+ * per pair, 8 floats: (u_t, v_t, kind, 0, e_x, e_y, e_z, 0). kind 0: the displacement to
54+ * the target is a plain difference (off-surface targets, or on-surface targets that do
55+ * not map into the patch's chart); 1: an on-surface target with preimage (u_t, v_t) in
56+ * the patch's chart and gap e = chart(u_t, v_t) - x_t; 2: the target's own patch.
57+ */
5358 pairF: Float32Array;
5459 npairs: number;
5560 }
@@ -150,11 +155,22 @@ export class Operator {
150155 async setup(): Promise<SetupInfo> {
151156 const t0 = now();
152157 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+ let gp: Float32Array;
159+ if (surf.exact) {
160+ // superellipsoid face charts: 4 vec4 per patch, [n, s0], [e1, s1], [e2, t0], [0, t1]
161+ gp = new Float32Array(16 * surf.npatches);
162+ for (let j = 0; j < surf.npatches; j++) {
163+ const par = surf.exact.params.subarray(13 * j, 13 * j + 13);
164+ 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);
165+ }
166+ } else {
167+ // Legendre coefficients of xyz, r_u, r_v: 9 * npols floats per patch, padded to vec4
168+ gp = new Float32Array(4 * Math.ceil((9 * npts) / 4));
169+ for (let j = 0; j < surf.npatches; j++) {
170+ for (let c = 0; c < 9; c++) {
171+ for (let n = 0; n < npols; n++) gp[j * 9 * npols + c * npols + n] = surf.srccoefs[c * npts + j * npols + n];
172+ }
173+ }
158174 }
159175 this.gpar = gpu.storage(gp);
160176 // umat and the quadrature rules in one table
@@ -253,7 +269,7 @@ export class Operator {
253269 const pList: number[] = [];
254270 const fList: number[] = [];
255271 const rowPtr = new Uint32Array(n + 1);
256- const params = this.surf.exact.params;
272+ const exact = this.surf.exact;
257273 for (let i = 0; i < n; i++) {
258274 const xi = x[i];
259275 const yi = x[n + i];
@@ -264,23 +280,27 @@ export class Operator {
264280 tList.push(i);
265281 pList.push(j);
266282 if (!onSurface) {
267- fList.push(0, 0, 0, 0);
283+ fList.push(0, 0, 0, 0, 0, 0, 0, 0);
268284 } else if (patchId[i] === j) {
269- fList.push(uvNodes[i], uvNodes[n + i], 2, 0);
270- } else {
285+ fList.push(uvNodes[i], uvNodes[n + i], 2, 0, 0, 0, 0, 0);
286+ } else if (exact) {
271287 // 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);
288+ const par = exact.params.subarray(13 * j, 13 * j + 13);
273289 const xn = xi * par[0] + yi * par[1] + zi * par[2];
274290 if (xn < 0.2) {
275- fList.push(0, 0, 0, 0);
291+ fList.push(0, 0, 0, 0, 0, 0, 0, 0);
276292 } else {
277293 const cx = xi / xn;
278294 const cy = yi / xn;
279295 const cz = zi / xn;
280296 const st = cx * par[3] + cy * par[4] + cz * par[5];
281297 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);
298+ fList.push((2 * (st - par[9])) / (par[10] - par[9]) - 1, (2 * (tt - par[11])) / (par[12] - par[11]) - 1, 1, 0, 0, 0, 0, 0);
283299 }
300+ } else {
301+ const pre = this.preimage(j, xi, yi, zi);
302+ if (pre) fList.push(pre[0], pre[1], 1, 0, pre[2], pre[3], pre[4], 0);
303+ else fList.push(0, 0, 0, 0, 0, 0, 0, 0);
284304 }
285305 }
286306 rowPtr[i + 1] = tList.length;
@@ -288,6 +308,85 @@ export class Operator {
288308 return { pairT: Int32Array.from(tList), pairP: Int32Array.from(pList), rowPtr, pairF: Float32Array.from(fList), npairs: tList.length };
289309 }
290310
311+ /**
312+ * For the interpolated geometry: the preimage (u, v) of the point x in the Legendre
313+ * interpolant of patch j, extended beyond [-1, 1]^2, by Gauss-Newton from the nearest
314+ * node, and the gap e = x(u, v) - x. Returns [u, v, e_x, e_y, e_z], or null when the
315+ * point does not map close to the patch (then a plain difference is used).
316+ */
317+ private preimage(j: number, x: number, y: number, z: number): number[] | null {
318+ const { npols, npts, norder } = this;
319+ const n1 = norder + 1;
320+ const sc = this.surf.srccoefs;
321+ const sv = this.surf.srcvals;
322+ const base = j * npols;
323+ let best = 0;
324+ let bd = Infinity;
325+ for (let l = 0; l < npols; l++) {
326+ const d = (sv[base + l] - x) ** 2 + (sv[npts + base + l] - y) ** 2 + (sv[2 * npts + base + l] - z) ** 2;
327+ if (d < bd) {
328+ bd = d;
329+ best = l;
330+ }
331+ }
332+ let u = this.surf.uvNodes[base + best];
333+ let v = this.surf.uvNodes[npts + base + best];
334+ const Pu = new Float64Array(n1);
335+ const Pv = new Float64Array(n1);
336+ const Du = new Float64Array(n1);
337+ const Dv = new Float64Array(n1);
338+ const leg = (t: number, P: Float64Array, Dd: Float64Array) => {
339+ P[0] = 1;
340+ Dd[0] = 0;
341+ if (norder >= 1) {
342+ P[1] = t;
343+ Dd[1] = 1;
344+ }
345+ for (let m = 1; m < norder; m++) {
346+ P[m + 1] = ((2 * m + 1) * t * P[m] - m * P[m - 1]) / (m + 1);
347+ Dd[m + 1] = Dd[m - 1] + (2 * m + 1) * P[m];
348+ }
349+ };
350+ const g = new Float64Array(9); // x, x_u, x_v of the xyz interpolant
351+ const evalAt = () => {
352+ leg(u, Pu, Du);
353+ leg(v, Pv, Dv);
354+ g.fill(0);
355+ for (let c = 0; c < 3; c++) {
356+ for (let b = 0; b < n1; b++) {
357+ for (let a = 0; a < n1; a++) {
358+ const cf = sc[c * npts + base + b * n1 + a];
359+ g[c] += cf * Pu[a] * Pv[b];
360+ g[3 + c] += cf * Du[a] * Pv[b];
361+ g[6 + c] += cf * Pu[a] * Dv[b];
362+ }
363+ }
364+ }
365+ };
366+ for (let it = 0; it < 30; it++) {
367+ evalAt();
368+ const r = [x - g[0], y - g[1], z - g[2]];
369+ const a11 = g[3] ** 2 + g[4] ** 2 + g[5] ** 2;
370+ const a12 = g[3] * g[6] + g[4] * g[7] + g[5] * g[8];
371+ const a22 = g[6] ** 2 + g[7] ** 2 + g[8] ** 2;
372+ const b1 = g[3] * r[0] + g[4] * r[1] + g[5] * r[2];
373+ const b2 = g[6] * r[0] + g[7] * r[1] + g[8] * r[2];
374+ const det = a11 * a22 - a12 * a12;
375+ if (!(det > 0)) return null;
376+ const du = (a22 * b1 - a12 * b2) / det;
377+ const dv = (a11 * b2 - a12 * b1) / det;
378+ u += du;
379+ v += dv;
380+ if (Math.abs(u) > 3 || Math.abs(v) > 3) return null;
381+ if (Math.abs(du) + Math.abs(dv) < 1e-13) break;
382+ }
383+ if (Math.abs(u) > 2 || Math.abs(v) > 2) return null;
384+ evalAt();
385+ const e = [g[0] - x, g[1] - y, g[2] - z];
386+ if (Math.hypot(e[0], e[1], e[2]) > 0.1 * this.surf.rads[j]) return null;
387+ return [u, v, ...e];
388+ }
389+
291390 /** The sinh-substituted segments along the four edges for the self patch of node i: 8 floats each. */
292391 private selfSegments(i: number): Float32Array {
293392 const { npts } = this;
@@ -374,7 +473,7 @@ export class Operator {
374473 targ: gpu.storage(tg),
375474 pairs,
376475 pairsI: gpu.storage(pi),
377- pairsF: gpu.storage(pairs.npairs ? pairs.pairF : new Float32Array(4)),
476+ pairsF: gpu.storage(pairs.npairs ? pairs.pairF : new Float32Array(8)),
378477 rowPtr: gpu.storage(pairs.rowPtr),
379478 pairP: gpu.storage(pairs.npairs ? new Uint32Array(pairs.pairP.buffer, pairs.pairP.byteOffset, pairs.npairs) : new Uint32Array(1)),
380479 W: gpu.storage(Math.max(16, wbytes)),
@@ -406,8 +505,12 @@ export class Operator {
406505 const t0 = now();
407506 // pass 1: leaf counts
408507 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);
508+ const ex = this.surf.exact;
509+ const ga = ex ? ex.a : 0;
510+ const gpe = ex ? ex.p : 0;
511+ const gmode = ex ? 0 : 1;
512+ const uw = new UniformWriter(11);
513+ uw.u32(n).u32(p0).u32(0).i32(leaf.nlmax).f32(leaf.rfacQ).f32(Math.abs(this.k)).f32(leaf.kRMax).f32(ga).i32(gpe).i32(gmode).i32(this.norder);
411514 gpu.write(this.uniformLeaves, new Uint8Array(uw.bytes));
412515 const dummyItems = gpu.storage(16);
413516 const dummyOffsets = gpu.storage(16);
@@ -417,26 +520,26 @@ export class Operator {
417520 const ptr = new Uint32Array(n + 1);
418521 for (let q = 0; q < n; q++) {
419522 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];
523+ const c = pairF[8 * p + 2] === 2 ? this.selfItems[set.pairs.pairT[p]].length / 8 : cnt[q];
524+ if (pairF[8 * p + 2] !== 2) nleaves += cnt[q];
422525 ptr[q + 1] = ptr[q] + c;
423526 }
424527 const nitems = ptr[n];
425528 const items = new Float32Array(8 * Math.max(1, nitems));
426529 for (let q = 0; q < n; q++) {
427530 const p = p0 + q;
428- if (pairF[4 * p + 2] === 2) items.set(this.selfItems[set.pairs.pairT[p]], 8 * ptr[q]);
531+ if (pairF[8 * p + 2] === 2) items.set(this.selfItems[set.pairs.pairT[p]], 8 * ptr[q]);
429532 }
430533 const itemsBuf = gpu.storage(items);
431534 const ptrBuf = gpu.storage(ptr);
432535 // 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);
536+ const uw2 = new UniformWriter(11);
537+ uw2.u32(n).u32(p0).u32(1).i32(leaf.nlmax).f32(leaf.rfacQ).f32(Math.abs(this.k)).f32(leaf.kRMax).f32(ga).i32(gpe).i32(gmode).i32(this.norder);
435538 gpu.write(this.uniformLeaves, new Uint8Array(uw2.bytes));
436539 gpu.dispatch(leavesPipe, gpu.bind(leavesPipe, [this.uniformLeaves, set.pairsI, set.pairsF, set.targ, this.gpar, counts, ptrBuf, itemsBuf]), Math.ceil(n / 64));
437540 // 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);
541+ const uw3 = new UniformWriter(11);
542+ uw3.u32(n).u32(p0).i32(this.norder).i32(this.npols).i32(quad.nq).i32(quad.nqr).i32(quad.nqx).f32(this.k).f32(ga).i32(gpe).i32(gmode);
440543 gpu.write(this.uniformWeights, new Uint8Array(uw3.bytes));
441544 gpu.dispatch(weightsPipe, gpu.bind(weightsPipe, [this.uniformWeights, set.pairsI, set.pairsF, set.targ, this.gpar, ptrBuf, itemsBuf, this.tables, set.W]), n);
442545 await gpu.sync();
src/gpu/shaders.tsmodified+92−12View file
@@ -5,11 +5,14 @@
55 * requires:
66 *
77 * 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.
8+ * a quadrature node is formed in the coordinates of the patch chart rather
9+ * than as a difference of two rounded positions, so that it carries full
10+ * relative precision however small it is. Without this the double-layer
11+ * kernels, whose numerators (y - x_t).n are O(r^2), lose most of their
12+ * digits near the target. There are two charts: the superellipsoid face
13+ * chart of the rounded cube (gmode 0), and for any other surface the
14+ * Legendre interpolant of the node data (gmode 1), where the displacement
15+ * comes from divided differences of the Legendre polynomials.
1316 * 2. The weakly singular difference D'_k - D'_kappa is evaluated from a short
1417 * power series when k r < 1, where the direct formula subtracts two nearly
1518 * equal O(1/r^3) terms.
@@ -137,7 +140,7 @@ fn legendre(x: f32, norder: i32, P: ptr<function, array<f32, ${NORD1_MAX}>>) {
137140 }
138141 `;
139142
140-/** The superellipsoid face chart; needs a gpar binding. */
143+/** The patch charts; they need a gpar binding. */
141144 const CHART = /* wgsl */ `
142145 struct Geo { x: vec3f, ru: vec3f, rv: vec3f, d: vec3f }
143146
@@ -175,6 +178,73 @@ fn chart_rel(ip: i32, u: f32, v: f32, ut: f32, vt: f32, rel: bool, xt: vec3f, ga
175178 return out;
176179 }
177180
181+// P_n(x) - P_n(xt) = (x - xt) * dd[n], from the three-term recurrence:
182+// dd[n+1] = ((2n+1) (P_n(xt) + x dd[n]) - n dd[n-1]) / (n+1).
183+fn legendre_dd(x: f32, norder: i32, Pt: ptr<function, array<f32, ${NORD1_MAX}>>, dd: ptr<function, array<f32, ${NORD1_MAX}>>) {
184+ (*dd)[0] = 0.0;
185+ if (norder >= 1) { (*dd)[1] = 1.0; }
186+ for (var n = 1; n < norder; n++) {
187+ (*dd)[n + 1] = ((2.0 * f32(n) + 1.0) * ((*Pt)[n] + x * (*dd)[n]) - f32(n) * (*dd)[n - 1]) / f32(n + 1);
188+ }
189+}
190+
191+fn gcoef(i: u32) -> f32 {
192+ return gpar[i >> 2u][i & 3u];
193+}
194+
195+// xyz, r_u, r_v of the Legendre interpolant of patch ip at (u, v), and d = y - x_t.
196+// gpar holds the coefficients of every patch, 9 * npols floats each (xyz, r_u, r_v).
197+// With rel set, d = x(u, v) - x(ut, vt) + e is formed with divided differences, where
198+// e = x(ut, vt) - x_t is the (small) gap between the target and the interpolant at
199+// its preimage, computed on the CPU; otherwise d = y - xt.
200+fn chart_poly(ip: i32, u: f32, v: f32, ut: f32, vt: f32, rel: bool, xt: vec3f, e: vec3f, norder: i32) -> Geo {
201+ let n1 = norder + 1;
202+ let npols = n1 * n1;
203+ let base = u32(ip) * u32(9 * npols);
204+ var Pu: array<f32, ${NORD1_MAX}>;
205+ var Pv: array<f32, ${NORD1_MAX}>;
206+ legendre(u, norder, &Pu);
207+ legendre(v, norder, &Pv);
208+ var out: Geo;
209+ out.x = vec3f(0.0); out.ru = vec3f(0.0); out.rv = vec3f(0.0);
210+ for (var b = 0; b < n1; b++) {
211+ for (var a = 0; a < n1; a++) {
212+ let P = Pu[a] * Pv[b];
213+ let n = u32(b * n1 + a);
214+ let np = u32(npols);
215+ out.x += P * vec3f(gcoef(base + n), gcoef(base + np + n), gcoef(base + 2u * np + n));
216+ out.ru += P * vec3f(gcoef(base + 3u * np + n), gcoef(base + 4u * np + n), gcoef(base + 5u * np + n));
217+ out.rv += P * vec3f(gcoef(base + 6u * np + n), gcoef(base + 7u * np + n), gcoef(base + 8u * np + n));
218+ }
219+ }
220+ if (rel) {
221+ var Put: array<f32, ${NORD1_MAX}>;
222+ var Pvt: array<f32, ${NORD1_MAX}>;
223+ var du: array<f32, ${NORD1_MAX}>;
224+ var dv: array<f32, ${NORD1_MAX}>;
225+ legendre(ut, norder, &Put);
226+ legendre(vt, norder, &Pvt);
227+ legendre_dd(u, norder, &Put, &du);
228+ legendre_dd(v, norder, &Pvt, &dv);
229+ // P_a(u) P_b(v) - P_a(ut) P_b(vt) = (u - ut) du_a P_b(v) + (v - vt) P_a(ut) dv_b
230+ let hu = u - ut;
231+ let hv = v - vt;
232+ var d = e;
233+ for (var b = 0; b < n1; b++) {
234+ for (var a = 0; a < n1; a++) {
235+ let w = hu * du[a] * Pv[b] + hv * Put[a] * dv[b];
236+ let n = u32(b * n1 + a);
237+ let np = u32(npols);
238+ d += w * vec3f(gcoef(base + n), gcoef(base + np + n), gcoef(base + 2u * np + n));
239+ }
240+ }
241+ out.d = d;
242+ } else {
243+ out.d = out.x - xt;
244+ }
245+ return out;
246+}
247+
178248 `;
179249
180250 /**
@@ -195,6 +265,8 @@ struct Params {
195265 kR_max: f32,
196266 ga: f32,
197267 gpe: i32,
268+ gmode: i32,
269+ norder: i32,
198270 }
199271 @group(0) @binding(0) var<uniform> U: Params;
200272 @group(0) @binding(1) var<storage, read> pairs_i: array<vec2i>;
@@ -214,7 +286,7 @@ fn build_leaves(@builtin(global_invocation_id) gid: vec3u) {
214286 let p = pl + U.pair_offset;
215287 let ip = pairs_i[p].y;
216288 let it = pairs_i[p].x;
217- let kind = i32(pairs_f[p].z);
289+ let kind = i32(pairs_f[2u * p].z);
218290 if (kind == 2) {
219291 if (U.fill == 0u) { counts[pl] = 0u; }
220292 return;
@@ -239,7 +311,11 @@ fn build_leaves(@builtin(global_invocation_id) gid: vec3u) {
239311 var xc: array<vec3f, 9>;
240312 var cen = vec3f(0.0);
241313 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);
314+ let uq = cu + h * f32(q % 3 - 1);
315+ let vq = cv + h * f32(q / 3 - 1);
316+ var g: Geo;
317+ if (U.gmode == 1) { g = chart_poly(ip, uq, vq, 0.0, 0.0, false, vec3f(0.0), vec3f(0.0), U.norder); }
318+ else { g = chart_rel(ip, uq, vq, 0.0, 0.0, false, vec3f(0.0), U.ga, U.gpe); }
243319 xc[q] = g.x;
244320 cen += g.x / 9.0;
245321 }
@@ -293,6 +369,7 @@ struct Params {
293369 zk: f32,
294370 ga: f32,
295371 gpe: i32,
372+ gmode: i32,
296373 }
297374 @group(0) @binding(0) var<uniform> U: Params;
298375 @group(0) @binding(1) var<storage, read> pairs_i: array<vec2i>;
@@ -320,9 +397,10 @@ fn near_weights(@builtin(workgroup_id) wid: vec3u, @builtin(local_invocation_id)
320397 let t = i32(lid.x);
321398 let it = pairs_i[p].x;
322399 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);
400+ let ut = pairs_f[2u * p].x;
401+ let vt = pairs_f[2u * p].y;
402+ let kind = i32(pairs_f[2u * p].z);
403+ let egap = pairs_f[2u * p + 1u].xyz;
326404 let is_self = kind == 2;
327405 let onsurf = kind >= 1;
328406 let xt = targ[2 * it].xyz;
@@ -372,7 +450,9 @@ fn near_weights(@builtin(workgroup_id) wid: vec3u, @builtin(local_invocation_id)
372450 v = it0.z + it0.w * tables[oq + b];
373451 wq = it0.w * it0.w * tables[oq + U.nq + a] * tables[oq + U.nq + b];
374452 }
375- let G = chart_rel(ip, u, v, ut, vt, onsurf, xt, U.ga, U.gpe);
453+ var G: Geo;
454+ if (U.gmode == 1) { G = chart_poly(ip, u, v, ut, vt, onsurf, xt, egap, U.norder); }
455+ else { G = chart_rel(ip, u, v, ut, vt, onsurf, xt, U.ga, U.gpe); }
376456 var ns = cross(G.ru, G.rv);
377457 let jac = length(ns);
378458 ns = ns / jac;
src/main.tsmodified+280−69View file
@@ -1,11 +1,15 @@
11 /**
22 * 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.
3+ * self-test, main problem, energy balance, field on three planes), and the
4+ * views. There are two problems: the room (a point source inside a rounded
5+ * cube with impedance walls, the interior equation) and the cymbal (a listener
6+ * near a rigid cymbal in free space, the exterior equation with lambda = 0,
7+ * solved with the source at the listener). The operator is kept between solves
8+ * when only the wall admittance or the source changes, since its near weights
9+ * depend only on k and the surface.
710 */
811 import { roundedCube, type Surface } from './surface.ts';
12+import { cymbal, CymbalShape, CYMBAL_DEFAULTS, insideTest, ringLayout } from './cymbal.ts';
913 import {
1014 ALPHA,
1115 energyBalance,
@@ -18,14 +22,18 @@ import {
1822 selftestRhs,
1923 slicePoints,
2024 type CVec,
25+ type Slice,
2126 } from './problem.ts';
27+import { BOX, boxSlices, C_SOUND, farFieldPower, LISTENER, RADIUS, selftestPoints as cymbalPoints, selftestSources, sourceSum } from './cymbalProblem.ts';
28+import { leggauss } from './gauss.ts';
2229 import { adapterDescription, Gpu } from './gpu/device.ts';
2330 import { Operator } from './gpu/operator.ts';
24-import { colorLimits, DB_MAX, DB_MIN, type Field, type Mode, type SliceField } from './view/field.ts';
31+import { colorLimits, type Field, type Mode, type SliceField } from './view/field.ts';
2532 import { SliceView } from './view/slices.ts';
2633 import { RoomView } from './view/room3d.ts';
2734 import { drawColorbar, DIVERGING, SEQUENTIAL } from './view/colormaps.ts';
2835 import { drawConvergence } from './view/convergence.ts';
36+import { fieldGif } from './view/exportGif.ts';
2937
3038 const $ = <T extends HTMLElement>(id: string): T => {
3139 const el = document.getElementById(id);
@@ -36,34 +44,66 @@ const $ = <T extends HTMLElement>(id: string): T => {
3644 const A = 1.0;
3745 const EXTERIOR_SOURCE = [1.7, 0.6, -0.5];
3846
47+type Problem = 'room' | 'cymbal';
48+
3949 interface Params {
50+ problem: Problem;
51+ // the room
4052 k: number;
4153 beta: number;
4254 m: number;
43- norder: number;
4455 p: number;
56+ // the cymbal: frequency (Hz), target patch size (cm), rim thickness (mm)
57+ freq: number;
58+ h: number;
59+ tRim: number;
60+ // both
61+ norder: number;
4562 x0: [number, number, number];
4663 eps: number;
4764 ngrid: number;
4865 selftest: boolean;
4966 }
5067
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 },
68+interface Preset {
69+ label: string;
70+ values: Partial<Params>;
71+}
72+
73+const PRESETS: Record<Problem, Record<string, Preset>> = {
74+ room: {
75+ quick: { label: 'quick (1,176 points)', values: { k: 6, beta: 0.2, m: 2, norder: 6 } },
76+ small: { label: 'small (4,704 points)', values: { k: 10, beta: 0.2, m: 4, norder: 6 } },
77+ medium: { label: 'medium (17,496 points)', values: { k: 10, beta: 0.2, m: 6, norder: 8 } },
78+ large: { label: 'large (48,600 points)', values: { k: 40, beta: 0.05, m: 10, norder: 8 } },
79+ },
80+ cymbal: {
81+ quick: { label: 'quick (4,800 points, 20 mm rim)', values: { freq: 1000, h: 10, tRim: 20, norder: 4 } },
82+ small: { label: 'small (19,012 points)', values: { freq: 1000, h: 6, tRim: 10, norder: 6 } },
83+ medium: { label: 'medium (35,964 points)', values: { freq: 3000, h: 4, tRim: 10, norder: 8 } },
84+ large: { label: 'large (42,444 points)', values: { freq: 8000, h: 3, tRim: 10, norder: 8 } },
85+ },
86+};
87+const COMMON: Omit<Params, 'problem' | 'x0'> = { k: 10, beta: 0.2, m: 4, p: 8, freq: 1000, h: 6, tRim: 10, norder: 6, eps: 1e-5, ngrid: 101, selftest: false };
88+const DEFAULTS: Record<Problem, Params> = {
89+ room: { ...COMMON, problem: 'room', x0: [-0.35, 0.2, 0.1] },
90+ cymbal: { ...COMMON, ...PRESETS.cymbal.small.values, problem: 'cymbal', x0: [...LISTENER] },
5691 };
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 };
5892
5993 const ui = {
94+ problem: $<HTMLSelectElement>('problem'),
6095 preset: $<HTMLSelectElement>('preset'),
6196 k: $<HTMLInputElement>('k'),
6297 beta: $<HTMLInputElement>('beta'),
6398 m: $<HTMLInputElement>('m'),
99+ freq: $<HTMLInputElement>('freq'),
100+ h: $<HTMLInputElement>('h'),
101+ tRim: $<HTMLInputElement>('tRim'),
64102 norder: $<HTMLSelectElement>('norder'),
65103 p: $<HTMLSelectElement>('p'),
66104 x0: [$<HTMLInputElement>('x0'), $<HTMLInputElement>('y0'), $<HTMLInputElement>('z0')],
105+ x0label: $('x0label'),
106+ x0name: $('x0name'),
67107 eps: $<HTMLSelectElement>('eps'),
68108 ngrid: $<HTMLSelectElement>('ngrid'),
69109 selftest: $<HTMLInputElement>('selftest'),
@@ -74,9 +114,13 @@ const ui = {
74114 status: $('status'),
75115 err: $('err'),
76116 mode: $<HTMLSelectElement>('mode'),
117+ levelOption: $<HTMLOptionElement>('levelOption'),
77118 animate: $<HTMLInputElement>('animate'),
78119 showWalls: $<HTMLInputElement>('showWalls'),
120+ wallsName: $('wallsName'),
79121 showPlanes: [$<HTMLInputElement>('showXY'), $<HTMLInputElement>('showXZ'), $<HTMLInputElement>('showYZ')],
122+ exportGif: $<HTMLButtonElement>('exportGif'),
123+ hint: $('hint'),
80124 results: $<HTMLTableElement>('results'),
81125 conv: $<HTMLCanvasElement>('conv'),
82126 empty: $('empty'),
@@ -90,15 +134,22 @@ const ui = {
90134 const fmtBytes = (b: number) => (b >= 2 ** 30 ? `${(b / 2 ** 30).toFixed(2)} GB` : `${(b / 2 ** 20).toFixed(0)} MB`);
91135 const fmtS = (s: number) => (s < 10 ? `${s.toFixed(2)} s` : `${s.toFixed(1)} s`);
92136
137+const cymbalShape = (p: Params) => new CymbalShape({ ...CYMBAL_DEFAULTS, tRim: p.tRim / 1000 });
138+const wavenumber = (p: Params) => (p.problem === 'room' ? p.k : (2 * Math.PI * p.freq) / C_SOUND);
139+
93140 // ------------------------------------------------------------------ params
94141
95142 function readParams(): Params {
96143 return {
144+ problem: ui.problem.value as Problem,
97145 k: Number(ui.k.value),
98146 beta: Number(ui.beta.value),
99147 m: Math.round(Number(ui.m.value)),
100- norder: Number(ui.norder.value),
101148 p: Number(ui.p.value),
149+ freq: Number(ui.freq.value),
150+ h: Number(ui.h.value),
151+ tRim: Number(ui.tRim.value),
152+ norder: Number(ui.norder.value),
102153 x0: ui.x0.map((e) => Number(e.value)) as [number, number, number],
103154 eps: Number(ui.eps.value),
104155 ngrid: Number(ui.ngrid.value),
@@ -107,11 +158,18 @@ function readParams(): Params {
107158 }
108159
109160 function writeParams(p: Params) {
161+ if (ui.problem.value !== p.problem) {
162+ ui.problem.value = p.problem;
163+ }
164+ applyProblem(p.problem);
110165 ui.k.value = String(p.k);
111166 ui.beta.value = String(p.beta);
112167 ui.m.value = String(p.m);
113- ui.norder.value = String(p.norder);
114168 ui.p.value = String(p.p);
169+ ui.freq.value = String(p.freq);
170+ ui.h.value = String(p.h);
171+ ui.tRim.value = String(p.tRim);
172+ ui.norder.value = String(p.norder);
115173 p.x0.forEach((v, i) => (ui.x0[i].value = String(+v.toFixed(3))));
116174 ui.eps.value = String(p.eps);
117175 if (![...ui.eps.options].some((o) => o.value === ui.eps.value)) ui.eps.value = '1e-5';
@@ -120,17 +178,50 @@ function writeParams(p: Params) {
120178 matchPreset();
121179 }
122180
181+/** Show the controls and labels of one problem, and its named cases. */
182+function applyProblem(problem: Problem) {
183+ const room = problem === 'room';
184+ document.querySelectorAll('.room-only').forEach((el) => el.classList.toggle('hidden', !room));
185+ document.querySelectorAll('.cymbal-only').forEach((el) => el.classList.toggle('hidden', room));
186+ if (ui.preset.dataset.problem !== problem) {
187+ ui.preset.replaceChildren(
188+ ...Object.entries(PRESETS[problem]).map(([v, pr]) => new Option(pr.label, v)),
189+ new Option('custom', 'custom'),
190+ );
191+ ui.preset.dataset.problem = problem;
192+ }
193+ ui.x0name.textContent = room ? 'source (x, y, z)' : 'listener (x, y, z) in m';
194+ ui.x0label.title = room
195+ ? 'Source position; must be inside the room. You can also click a slice.'
196+ : 'Listener position, where the source is placed (reciprocity); must be outside the cymbal. You can also click a slice.';
197+ ui.wallsName.textContent = room ? 'walls' : 'cymbal';
198+ ui.levelOption.textContent = room ? 'level (dB re median)' : 'level (dB re free field)';
199+ ui.hint.textContent = room ? 'drag to turn · scroll to zoom · the nearer walls are cut away' : 'drag to turn · scroll to zoom';
200+}
201+
202+function presetMatches(p: Params, v: Partial<Params>): boolean {
203+ return Object.entries(v).every(([key, val]) => p[key as keyof Params] === val);
204+}
205+
123206 function matchPreset() {
124207 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];
208+ const name = Object.entries(PRESETS[p.problem]).find(([, pr]) => presetMatches(p, pr.values))?.[0];
126209 ui.preset.value = name ?? 'custom';
127210 }
128211
129212 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';
213+ if (p.problem === 'room') {
214+ if (!(p.k > 0 && p.k <= 80)) return 'the wavenumber must be between 0 and 80';
215+ if (!(p.beta > 0)) return 'the wall admittance must be positive (absorbing walls)';
216+ if (!(p.m >= 1 && p.m <= 16)) return 'patches per face edge must be between 1 and 16';
217+ if (!inside(p.x0, A, p.p, 0.04)) return 'the source must be inside the room, a little away from the walls';
218+ return null;
219+ }
220+ if (!(p.freq >= 100 && p.freq <= 16000)) return 'the frequency must be between 100 Hz and 16 kHz';
221+ if (!(p.h >= 2 && p.h <= 20)) return 'the patch size must be between 2 and 20 cm';
222+ if (!(p.tRim >= 4 && p.tRim <= 40)) return 'the rim thickness must be between 4 and 40 mm';
223+ if (!p.x0.every((v) => Math.abs(v) <= 3)) return 'the listener must be within 3 m of the cymbal';
224+ if (insideTest(cymbalShape(p), 0.01)(p.x0)) return 'the listener must be outside the cymbal, at least 1 cm from it';
134225 return null;
135226 }
136227
@@ -139,25 +230,39 @@ function validate(p: Params): string | null {
139230 let deviceLimit = Infinity;
140231
141232 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);
233+ const npols = (p.norder + 1) ** 2;
145234 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;
235+ let npts: number;
236+ let perWl: number;
237+ let perNode: number;
238+ let size: string;
239+ if (p.problem === 'room') {
240+ npts = 6 * p.m * p.m * npols;
241+ const wl = (2 * Math.PI) / p.k;
242+ perWl = ((p.norder + 1) * p.m * wl) / (2 * A);
243+ // 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
244+ perNode = Math.min(6 * p.m * p.m, rfac >= 2 ? 9.5 : 5.8);
245+ size = `room <b>${((2 * A) / wl).toFixed(1)}</b> wavelengths wide`;
246+ } else {
247+ const L = ringLayout(cymbalShape(p), p.h / 100);
248+ const npatches = 2 * (L.mC ** 2 + L.counts.reduce((a, b) => a + b, 0)) + L.nrim;
249+ npts = npatches * npols;
250+ const wl = C_SOUND / p.freq;
251+ perWl = ((p.norder + 1) * wl) / (p.h / 100);
252+ // measured on the named cases: about 15 near patches per node at rfac 2 and 5 at rfac 1.25 (the thin body)
253+ perNode = rfac >= 2 ? 15 : 5.5;
254+ size = `cymbal <b>${((2 * RADIUS) / wl).toFixed(1)}</b> wavelengths across`;
255+ }
256+ const wbytes = npts * perNode * 6 * npols * 8;
149257 const parts = [
150258 `<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>`,
259+ size,
260+ `<span class="${perWl < 5 ? 'warn' : ''}" title="nodes per wavelength along the surface, on average">${perWl.toFixed(1)} nodes per wavelength</span>`,
153261 `<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>`,
154262 ];
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- }
263+ if (wbytes > 1.3 * deviceLimit) parts.push(`<span class="bad">probably more than this GPU allows (${fmtBytes(deviceLimit)})</span>`);
160264 if (perWl < 5) parts.push('<span class="warn">under-resolved: expect large errors</span>');
265+ let blocking = false;
161266 const err = validate(p);
162267 if (err) {
163268 parts.push(`<span class="bad">${err}</span>`);
@@ -209,6 +314,9 @@ function showResults(rows: Row[]) {
209314
210315 const quality = (e: number, tol: number) => (e < 10 * tol ? '' : e < 1e-2 ? 'warn' : 'bad');
211316
317+/** Pause so that the page can repaint the status line before a long CPU step. */
318+const yieldToPage = () => new Promise((r) => setTimeout(r, 0));
319+
212320 // ------------------------------------------------------------------ main
213321
214322 async function start() {
@@ -226,14 +334,14 @@ async function start() {
226334 });
227335 deviceLimit = Math.min(device.limits.maxStorageBufferBindingSize, device.limits.maxBufferSize);
228336 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.' : '');
337+ ui.device.textContent = `WebGPU adapter: ${adapterDescription(adapter)}; largest buffer ${fmtBytes(deviceLimit)}.` + (software ? ' This is a software adapter, so only the quick cases are practical.' : '');
230338 const gpu = new Gpu(device, adapter);
231339
232340 const slices = new SliceView(
233341 $('slices'),
234342 (x) => {
235343 const p = readParams();
236- if (!inside(x, A, p.p, 0.04)) return;
344+ if (validate({ ...p, x0: x })) return;
237345 x.forEach((v, i) => (ui.x0[i].value = String(+v.toFixed(3))));
238346 updateDerived();
239347 if (!running) void solve();
@@ -248,20 +356,25 @@ async function start() {
248356 let field: Field | null = null;
249357 let conv: { label: string; res: number[]; color: string }[] = [];
250358 let lastEps = 1e-5;
359+ let fieldTitle = '';
360+
361+ /** The three planes as a GIF, in the current display mode. */
362+ const buildGif = () => (field ? fieldGif(field, ui.mode.value as Mode, { title: fieldTitle }) : null);
251363
252364 const redrawConv = () => drawConvergence(ui.conv, conv, lastEps);
253365
254366 const updateColorbar = () => {
255367 const mode = ui.mode.value as Mode;
256- drawColorbar(ui.cbar, mode === 'pressure' ? DIVERGING : SEQUENTIAL);
368+ const diverging = mode === 'pressure' || (field?.level.diverging ?? false);
369+ drawColorbar(ui.cbar, diverging ? DIVERGING : SEQUENTIAL);
257370 if (!field) {
258371 ui.cbarHi.textContent = ui.cbarLo.textContent = '';
259372 } else if (mode === 'pressure') {
260373 ui.cbarHi.textContent = `+${field.umax.toPrecision(2)}`;
261374 ui.cbarLo.textContent = `−${field.umax.toPrecision(2)}`;
262375 } else {
263- ui.cbarHi.textContent = `+${DB_MAX} dB`;
264- ui.cbarLo.textContent = `${DB_MIN} dB`;
376+ ui.cbarHi.textContent = `+${field.level.max} dB`;
377+ ui.cbarLo.textContent = `${field.level.min} dB`;
265378 }
266379 };
267380
@@ -281,6 +394,7 @@ async function start() {
281394 const rows: Row[] = [];
282395 // the previous solution no longer matches the controls, so take it down
283396 field = null;
397+ ui.exportGif.disabled = true;
284398 slices.clear();
285399 room.clear();
286400 ui.empty.textContent = 'Solving…';
@@ -291,17 +405,20 @@ async function start() {
291405 redrawConv();
292406 try {
293407 const t0 = performance.now();
294- const key = `${p.k}|${p.m}|${p.norder}|${p.p}`;
295- const lam = -p.beta;
408+ const isRoom = p.problem === 'room';
409+ const k = wavenumber(p);
410+ const key = isRoom ? `room|${p.k}|${p.m}|${p.norder}|${p.p}` : `cymbal|${p.freq}|${p.h}|${p.norder}|${p.tRim}`;
411+ const lam = isRoom ? -p.beta : 0;
412+ const shape = isRoom ? null : cymbalShape(p);
296413 let setupS = 0;
297414 if (!op || key !== opKey) {
298415 op?.destroy();
299416 op = null;
300417 opKey = '';
301418 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, {
419+ await yieldToPage();
420+ surf = isRoom ? roundedCube(A, p.p, p.m, p.norder) : cymbal(shape!, p.h / 100, p.norder);
421+ const o = new Operator(gpu, surf, k, ALPHA, lam, isRoom ? 0 : 1, {
305422 progress: (stage, f) => {
306423 checkCancel();
307424 setStatus(`${stage}: ${(100 * f).toFixed(0)}%`, f);
@@ -329,24 +446,38 @@ async function start() {
329446 if (p.selftest) {
330447 const series = { label: 'self-test', res: [] as number[], color: '#8b949e' };
331448 conv.push(series);
332- const sol = await o.solve(selftestRhs(xs, nrm, EXTERIOR_SOURCE, p.k, lam), p.eps, 400, (it, r) => {
449+ let rhs: CVec;
450+ let exact: (x: Float64Array) => CVec;
451+ let xt: Float64Array;
452+ if (isRoom) {
453+ rhs = selftestRhs(xs, nrm, EXTERIOR_SOURCE, k, lam);
454+ exact = (x) => green(x, EXTERIOR_SOURCE, k);
455+ xt = selftestPoints(50, A);
456+ } else {
457+ setStatus('self-test: the data of the interior sources…');
458+ await yieldToPage();
459+ const { x0s, ws } = selftestSources(shape!);
460+ rhs = sourceSum(xs, x0s, ws, k, nrm);
461+ exact = (x) => sourceSum(x, x0s, ws, k);
462+ xt = cymbalPoints(shape!, 50);
463+ }
464+ const sol = await o.solve(rhs, p.eps, 400, (it, r) => {
333465 checkCancel();
334466 series.res.push(r);
335467 redrawConv();
336468 setStatus(`self-test GMRES: iteration ${it}, residual ${r.toExponential(1)}`);
337469 });
338470 setStatus('self-test: evaluating…');
339- const xt = selftestPoints(50, A);
340471 const ut = await o.eval(xt, sol.x, sol.phi1, 'self-test field');
341472 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)]);
473+ const e1 = relError(ut, exact(xt));
474+ const e2 = relError(ub, exact(xs));
475+ rows.push([isRoom ? 'self-test error, interior · walls' : 'self-test error, around · on the cymbal', `${e1.toExponential(1)} · ${e2.toExponential(1)}`, quality(Math.max(e1, e2), p.eps)]);
345476 }
346477
347- const series = { label: 'room problem', res: [] as number[], color: '#2f81f7' };
478+ const series = { label: isRoom ? 'room problem' : 'cymbal problem', res: [] as number[], color: '#2f81f7' };
348479 conv.push(series);
349- const sol = await o.solve(roomRhs(xs, nrm, p.x0, p.k, lam), p.eps, 400, (it, r) => {
480+ const sol = await o.solve(roomRhs(xs, nrm, p.x0, k, lam), p.eps, 400, (it, r) => {
350481 checkCancel();
351482 series.res.push(r);
352483 redrawConv();
@@ -357,37 +488,73 @@ async function start() {
357488
358489 // energy balance
359490 setStatus('energy balance…');
491+ await yieldToPage();
360492 const v0 = await o.eval(Float64Array.from(p.x0), sol.x, sol.phi1, 'field at the source');
361493 const vw = await o.evalOnSurface(sol.x, sol.phi1);
362- const gw = green(xs, p.x0, p.k);
494+ const gw = green(xs, p.x0, k);
363495 const uWall = new Float64Array(2 * npts);
364496 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)]);
497+ if (isRoom) {
498+ const eb = energyBalance(k, lam, [v0[0], v0[1]], s.wts, uWall);
499+ rows.push(['power: source · walls', `${eb.pSrc.toFixed(6)} · ${eb.pWall.toFixed(6)}`]);
500+ rows.push(['energy mismatch', eb.mismatch.toExponential(1), quality(eb.mismatch, p.eps)]);
501+ } else {
502+ const pSrc = k / (4 * Math.PI) + v0[1];
503+ // the far field pattern costs (number of directions) x (number of nodes) on the CPU
504+ const L = Math.ceil(k * 0.6) + 20;
505+ if ((L + 1) * (2 * L + 2) * npts < 2e8) {
506+ const pFar = farFieldPower(k, p.x0, xs, nrm, s.wts, sol.x, sol.phi1, ALPHA, leggauss);
507+ const mismatch = Math.abs(pSrc - pFar) / pSrc;
508+ rows.push(['power: source · to infinity', `${pSrc.toFixed(6)} · ${pFar.toFixed(6)}`]);
509+ rows.push(['energy mismatch', mismatch.toExponential(1), quality(mismatch, p.eps)]);
510+ } else {
511+ rows.push(['power from the source', `${pSrc.toFixed(6)} (the far field check is skipped at this size)`]);
512+ }
513+ }
368514
369- // the field on the three planes through the source
515+ // the field on the three planes
370516 const tf = performance.now();
371- const { grid, slices: sl } = slicePoints(p.ngrid, p.x0, A, p.p);
517+ const sl: Slice[] = isRoom ? slicePoints(p.ngrid, p.x0, A, p.p).slices : boxSlices(p.ngrid, p.x0, shape!);
372518 const sfs: SliceField[] = [];
373519 for (const [si, slice] of sl.entries()) {
374520 const { pts, index } = maskedPoints(slice.pts, slice.mask);
375521 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);
522+ const g = green(pts, p.x0, k);
377523 const u = new Float64Array(2 * slice.mask.length).fill(NaN);
524+ const ref = isRoom ? undefined : new Float64Array(slice.mask.length).fill(NaN);
378525 for (let q = 0; q < index.length; q++) {
379526 u[2 * index[q]] = g[2 * q] + v[2 * q];
380527 u[2 * index[q] + 1] = g[2 * q + 1] + v[2 * q + 1];
528+ if (ref) ref[index[q]] = Math.hypot(g[2 * q], g[2 * q + 1]);
381529 }
382- sfs.push({ slice, u });
530+ sfs.push({ slice, u, ref });
383531 }
384532 const nfield = sl.reduce((acc, x) => acc + x.mask.reduce((a, b) => a + b, 0), 0);
385533 rows.push(['field on the planes', `${fmtS((performance.now() - tf) / 1000)} (${nfield.toLocaleString()} points)`]);
386534 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 };
535+ const lims = colorLimits(sfs, p.x0, k);
536+ if (isRoom) {
537+ field = {
538+ k, x0: [...p.x0], ngrid: p.ngrid, slices: sfs, uWall, freeFieldLevel: false, ...lims,
539+ level: { min: -30, max: 15, diverging: false, label: 're median' },
540+ scene: { center: [0, 0, 0], radius: 1, cutaway: true },
541+ };
542+ } else {
543+ const [r, z] = shape!.crossSection();
544+ const c = BOX.map(([lo, hi]) => (lo + hi) / 2) as [number, number, number];
545+ field = {
546+ k, x0: [...p.x0], ngrid: p.ngrid, slices: sfs, uWall, freeFieldLevel: true, ...lims,
547+ level: { min: -15, max: 15, diverging: true, label: 're free field' },
548+ scene: { center: c, radius: 0.45, cutaway: false },
549+ outline: { r, z, radius: RADIUS },
550+ };
551+ }
552+ fieldTitle = isRoom
553+ ? `Point source at (${p.x0.join(', ')}) in a room, k = ${p.k}, wall admittance ${p.beta}`
554+ : `Listener at (${p.x0.join(', ')}) m near a rigid cymbal, ${p.freq} Hz; u(x) is the pressure at the listener due to a source at x`;
389555 slices.setField(field);
390556 room.setField(s, field);
557+ ui.exportGif.disabled = false;
391558 ui.empty.style.display = 'none';
392559 updateColorbar();
393560 showResults(rows);
@@ -420,16 +587,38 @@ async function start() {
420587 }
421588
422589 ui.solve.addEventListener('click', () => void solve());
590+ ui.exportGif.addEventListener('click', () => {
591+ const p = readParams();
592+ setStatus('writing the GIF…');
593+ setTimeout(() => {
594+ const gif = buildGif();
595+ if (!gif) return;
596+ const url = URL.createObjectURL(new Blob([gif as Uint8Array<ArrayBuffer>], { type: 'image/gif' }));
597+ const a = document.createElement('a');
598+ a.href = url;
599+ const what = p.problem === 'room' ? `room-k${p.k}` : `cymbal-${p.freq}Hz`;
600+ a.download = `fmm3dbie-${what}-${ui.mode.value}.gif`;
601+ a.click();
602+ setTimeout(() => URL.revokeObjectURL(url), 10_000);
603+ setStatus(`wrote ${a.download} (${(gif.length / 2 ** 20).toFixed(1)} MB)`, 1);
604+ }, 0);
605+ });
423606 ui.cancel.addEventListener('click', () => {
424607 cancelled = true;
425608 setStatus('cancelling…');
426609 });
610+ ui.problem.addEventListener('change', () => {
611+ const problem = ui.problem.value as Problem;
612+ const d = DEFAULTS[problem];
613+ writeParams({ ...d, x0: [...d.x0], ...(software ? PRESETS[problem].quick.values : {}) });
614+ updateDerived();
615+ });
427616 ui.preset.addEventListener('change', () => {
428- const pr = PRESETS[ui.preset.value];
429- if (pr) writeParams({ ...readParams(), ...pr });
617+ const pr = PRESETS[readParams().problem][ui.preset.value];
618+ if (pr) writeParams({ ...readParams(), ...pr.values });
430619 updateDerived();
431620 });
432- for (const el of [ui.k, ui.beta, ui.m, ui.norder, ui.p, ...ui.x0, ui.eps, ui.ngrid, ui.selftest]) {
621+ for (const el of [ui.k, ui.beta, ui.m, ui.freq, ui.h, ui.tRim, ui.norder, ui.p, ...ui.x0, ui.eps, ui.ngrid, ui.selftest]) {
433622 el.addEventListener('input', () => {
434623 matchPreset();
435624 updateDerived();
@@ -457,25 +646,47 @@ async function start() {
457646 slicesDirty = false;
458647 }
459648 // ?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);
649+ if (!noRoom) {
650+ // the checkboxes by plane name, in the order of the field's planes
651+ const byName = { xy: ui.showPlanes[0].checked, xz: ui.showPlanes[1].checked, yz: ui.showPlanes[2].checked };
652+ const planes = (field?.slices.map((sf) => byName[sf.slice.name]) ?? [false, false, false]) as [boolean, boolean, boolean];
653+ room.draw(mode, ph, { walls: ui.showWalls.checked, planes }, animating);
654+ }
461655 requestAnimationFrame(frame);
462656 };
463657 let slicesDirty = true;
464- const noRoom = location.search.includes('noroom');
658+ // (the test hook __app__.showRoom() turns the view on after a solve, for screenshots)
659+ let noRoom = location.search.includes('noroom');
465660 ui.mode.addEventListener('change', () => (slicesDirty = true));
466661 ui.animate.addEventListener('change', () => (slicesDirty = true));
467662 requestAnimationFrame(frame);
468663
469- writeParams({ ...DEFAULTS, x0: [...DEFAULTS.x0] });
664+ // ?problem=cymbal selects the cymbal at load, ?case=<name> a named case, ?autosolve solves at once
665+ const query = new URLSearchParams(location.search);
666+ const initial: Problem = query.get('problem') === 'cymbal' ? 'cymbal' : 'room';
667+ const d0 = DEFAULTS[initial];
668+ writeParams({ ...d0, x0: [...d0.x0], ...(software ? PRESETS[initial].quick.values : {}) });
470669 updateDerived();
471670 updateColorbar();
472671 redrawConv();
473- setStatus(software ? 'ready (software adapter: choose the quick case)' : 'ready', 0);
474- if (software) {
475- writeParams({ ...readParams(), ...PRESETS.quick });
672+ const named = PRESETS[initial][query.get('case') ?? ''];
673+ if (named) {
674+ writeParams({ ...readParams(), ...named.values });
476675 updateDerived();
477676 }
478- (window as unknown as { __app__: unknown }).__app__ = { solve, get field() { return field; } };
677+ setStatus(software ? 'ready (software adapter: choose the quick case)' : 'ready', 0);
678+ if (query.has('autosolve')) void solve();
679+ (window as unknown as { __app__: unknown }).__app__ = {
680+ solve,
681+ get field() {
682+ return field;
683+ },
684+ showRoom() {
685+ noRoom = false;
686+ room.invalidate();
687+ },
688+ buildGif,
689+ };
479690 }
480691
481692 start().catch((e) => {
src/problem.tsmodified+6−1View file
@@ -93,7 +93,12 @@ export interface Slice {
9393 name: 'xy' | 'xz' | 'yz';
9494 axes: [number, number];
9595 fixed: number;
96+ /** the value of the fixed coordinate */
97+ value: number;
98+ /** [lo, hi] along the two axes */
99+ extent: [[number, number], [number, number]];
96100 pts: Float64Array; // (3, ngrid^2), row-major over the grid with the first axis varying fastest
101+ /** 1 where the field is computed (inside the room, or outside the cymbal) */
97102 mask: Uint8Array;
98103 }
99104
@@ -124,7 +129,7 @@ export function slicePoints(ngrid: number, x0: ArrayLike<number>, a: number, p:
124129 mask[q] = inside(x, a, p, 1e-3) ? 1 : 0;
125130 }
126131 }
127- slices.push({ name, axes, fixed, pts, mask });
132+ slices.push({ name, axes, fixed, value: x0[fixed], extent: [[-a, a], [-a, a]], pts, mask });
128133 }
129134 return { grid, slices };
130135 }
src/surface.tsmodified+22−5View file
@@ -159,9 +159,10 @@ export class Surface {
159159 readonly uvNodes: Float64Array; // (2, npts)
160160 readonly norder: number;
161161 readonly srcvals: Float64Array; // (12, npts): xyz, r_u, r_v, unit normal
162- readonly exact: ExactGeometry;
162+ /** the analytic chart of every patch, or null for the interpolated geometry */
163+ readonly exact: ExactGeometry | null;
163164
164- constructor(norder: number, srcvals: Float64Array, exact: ExactGeometry) {
165+ constructor(norder: number, srcvals: Float64Array, exact: ExactGeometry | null = null) {
165166 this.norder = norder;
166167 this.srcvals = srcvals;
167168 this.exact = exact;
@@ -196,7 +197,7 @@ export class Surface {
196197 this.jac[i] = Math.hypot(cx, cy, cz);
197198 this.wts[i] = this.jac[i] * w[i % npols];
198199 }
199- // centroid and radius of each patch from a 5x5 sample grid of the exact chart
200+ // centroid and radius of each patch from a 5x5 sample grid of its chart
200201 const g5 = new Float64Array(50);
201202 for (let i = 0; i < 5; i++) {
202203 for (let j = 0; j < 5; j++) {
@@ -228,10 +229,26 @@ export class Surface {
228229 }
229230 }
230231
231- /** xyz, r_u, r_v (9, n) of the exact chart of one patch at reference points uv (2, n). */
232+ /**
233+ * xyz, r_u, r_v (9, n) of one patch at reference points uv (2, n): the exact chart
234+ * when there is one, otherwise the Legendre interpolant of the node data, as in
235+ * fmm3dbie.
236+ */
232237 evalGeometry(patch: number, uv: Float64Array): Float64Array {
233238 const e = this.exact;
234- return superellipsoidChart(e.params.subarray(13 * patch, 13 * patch + 13), e.a, e.p, uv);
239+ if (e) return superellipsoidChart(e.params.subarray(13 * patch, 13 * patch + 13), e.a, e.p, uv);
240+ const n = uv.length / 2;
241+ const { npols, npts } = this;
242+ const B = basis2d(uv, this.norder);
243+ const out = new Float64Array(9 * n);
244+ const base = patch * npols;
245+ for (let c = 0; c < 9; c++) {
246+ for (let m = 0; m < npols; m++) {
247+ const cm = this.srccoefs[c * npts + base + m];
248+ for (let q = 0; q < n; q++) out[c * n + q] += cm * B[m * n + q];
249+ }
250+ }
251+ return out;
235252 }
236253
237254 /** Nodes of the order-nover rule on every patch: xyz (3, N), unit normal (3, N), weights (N). */
src/view/exportGif.tsadded+196−0View file
@@ -0,0 +1,196 @@
1+/**
2+ * The three planes side by side as one GIF: animated over one period for the
3+ * pressure (24 frames, looping), a single frame for the level. The planes are
4+ * sampled per pixel from the complex field (bilinear, as in the 3D view), so
5+ * each frame maps directly to colormap indices. Captions, the source marker,
6+ * the cymbal's outline and the colorbar labels are drawn once on a canvas and
7+ * overlaid.
8+ */
9+import { DIVERGING, SEQUENTIAL } from './colormaps.ts';
10+import { encodeGif } from './gif.ts';
11+import type { Field, Mode } from './field.ts';
12+import { drawOutline } from './slices.ts';
13+
14+const AXIS = ['x', 'y', 'z'];
15+/** colormap levels in the palette; the rest are fixed colors */
16+const NCMAP = 240;
17+const FIXED: [number, number, number][] = [
18+ [255, 255, 255], // 240 background
19+ [17, 17, 17], // 241 black
20+ [40, 40, 40], // 242 cymbal
21+ [226, 229, 233], // 243 no field (outside the room, inside the cymbal)
22+ [87, 96, 106], // 244 text
23+ [160, 166, 173], // 245 frame
24+];
25+const BG = NCMAP;
26+const NOFIELD = NCMAP + 3;
27+
28+export interface GifOptions {
29+ /** panel size in pixels */
30+ size?: number;
31+ frames?: number;
32+ title?: string;
33+}
34+
35+export function fieldGif(f: Field, mode: Mode, opts: GifOptions = {}): Uint8Array {
36+ const S = opts.size ?? 300;
37+ const nframes = mode === 'pressure' ? (opts.frames ?? 24) : 1;
38+ const pad = 12;
39+ const titleH = opts.title ? 24 : 0;
40+ const capH = 20;
41+ const cbarH = 36;
42+ const np = f.slices.length;
43+ const W = pad + np * (S + pad);
44+ const H = pad + titleH + S + capH + cbarH + pad;
45+ const diverging = mode === 'pressure' || f.level.diverging;
46+ const lut = diverging ? DIVERGING : SEQUENTIAL;
47+
48+ // the palette: colormap levels, then the fixed colors
49+ const palette = new Uint8Array(768);
50+ for (let i = 0; i < NCMAP; i++) {
51+ const j = Math.round((i * 255) / (NCMAP - 1)) * 4;
52+ palette.set([lut[j], lut[j + 1], lut[j + 2]], 3 * i);
53+ }
54+ FIXED.forEach((c, i) => palette.set(c, 3 * (NCMAP + i)));
55+
56+ // per pixel of each panel: the interpolated field (re, im) and level reference, NaN where there is none
57+ const x0p = pad;
58+ const y0p = pad + titleH;
59+ const re = new Float32Array(np * S * S);
60+ const im = new Float32Array(np * S * S);
61+ const lv = new Float32Array(np * S * S);
62+ f.slices.forEach(({ u, ref, slice }, s) => {
63+ const n = f.ngrid;
64+ for (let py = 0; py < S; py++) {
65+ const gy = (1 - py / (S - 1)) * (n - 1);
66+ const r0 = Math.min(n - 2, Math.floor(gy));
67+ const fy = gy - r0;
68+ for (let px = 0; px < S; px++) {
69+ const gx = (px / (S - 1)) * (n - 1);
70+ const c0 = Math.min(n - 2, Math.floor(gx));
71+ const fx = gx - c0;
72+ const q = [r0 * n + c0, r0 * n + c0 + 1, (r0 + 1) * n + c0, (r0 + 1) * n + c0 + 1];
73+ const w = [(1 - fx) * (1 - fy), fx * (1 - fy), (1 - fx) * fy, fx * fy];
74+ let a = 0;
75+ let b = 0;
76+ let rf = 0;
77+ let ok = true;
78+ for (let t = 0; t < 4; t++) {
79+ if (!slice.mask[q[t]]) ok = false;
80+ a += w[t] * u[2 * q[t]];
81+ b += w[t] * u[2 * q[t] + 1];
82+ rf += w[t] * (ref ? ref[q[t]] : f.ref);
83+ }
84+ const o = s * S * S + py * S + px;
85+ re[o] = ok ? a : NaN;
86+ im[o] = ok ? b : NaN;
87+ lv[o] = ok ? (20 * Math.log10(Math.hypot(a, b) / rf) - f.level.min) / (f.level.max - f.level.min) : NaN;
88+ }
89+ }
90+ });
91+
92+ // the static overlay
93+ const cv = new OffscreenCanvas(W, H);
94+ const ctx = cv.getContext('2d')!;
95+ ctx.font = '12px system-ui, sans-serif';
96+ ctx.textBaseline = 'middle';
97+ if (opts.title) {
98+ ctx.fillStyle = 'rgb(17,17,17)';
99+ ctx.font = '13px system-ui, sans-serif';
100+ ctx.fillText(opts.title, pad, pad + titleH / 2 - 2);
101+ ctx.font = '12px system-ui, sans-serif';
102+ }
103+ f.slices.forEach(({ slice }, s) => {
104+ const left = x0p + s * (S + pad);
105+ const [ea, eb] = slice.extent;
106+ const X = (a: number) => left + ((a - ea[0]) / (ea[1] - ea[0])) * (S - 1);
107+ const Y = (b: number) => y0p + (1 - (b - eb[0]) / (eb[1] - eb[0])) * (S - 1);
108+ ctx.save();
109+ ctx.beginPath();
110+ ctx.rect(left, y0p, S, S);
111+ ctx.clip();
112+ if (f.outline) drawOutline(ctx as unknown as CanvasRenderingContext2D, slice, f.outline, X, Y);
113+ if (Math.abs(f.x0[slice.fixed] - slice.value) < 1e-9) {
114+ const [i, j] = slice.axes;
115+ ctx.beginPath();
116+ ctx.arc(X(f.x0[i]), Y(f.x0[j]), 5, 0, 2 * Math.PI);
117+ ctx.lineWidth = 1.5;
118+ ctx.strokeStyle = 'rgb(17,17,17)';
119+ ctx.stroke();
120+ ctx.beginPath();
121+ ctx.arc(X(f.x0[i]), Y(f.x0[j]), 6.5, 0, 2 * Math.PI);
122+ ctx.lineWidth = 1;
123+ ctx.strokeStyle = 'rgb(255,255,255)';
124+ ctx.stroke();
125+ }
126+ ctx.restore();
127+ ctx.strokeStyle = 'rgb(160,166,173)';
128+ ctx.lineWidth = 1;
129+ ctx.strokeRect(left - 0.5, y0p - 0.5, S + 1, S + 1);
130+ ctx.fillStyle = 'rgb(87,96,106)';
131+ const [i, j] = slice.axes;
132+ ctx.fillText(`${slice.name}-plane at ${AXIS[slice.fixed]} = ${slice.value.toFixed(2)} (${AXIS[i]} across, ${AXIS[j]} up)`, left, y0p + S + capH / 2 + 2);
133+ });
134+ // colorbar: a strip of colormap levels, with its limits
135+ const cbW = Math.min(360, W - 2 * pad - 180);
136+ const cbX = Math.round((W - cbW) / 2);
137+ const cbY = y0p + S + capH + 6;
138+ const cbH = 10;
139+ ctx.fillStyle = 'rgb(87,96,106)';
140+ const lo = mode === 'pressure' ? `−${f.umax.toPrecision(2)}` : `${f.level.min} dB`;
141+ const hi = mode === 'pressure' ? `+${f.umax.toPrecision(2)}` : `+${f.level.max} dB`;
142+ const label = mode === 'pressure' ? 'pressure Re(u e^−iωt)' : `level (dB ${f.level.label})`;
143+ ctx.textAlign = 'right';
144+ ctx.fillText(lo, cbX - 6, cbY + cbH / 2);
145+ ctx.textAlign = 'left';
146+ ctx.fillText(hi, cbX + cbW + 6, cbY + cbH / 2);
147+ ctx.textAlign = 'center';
148+ ctx.fillText(label, cbX + cbW / 2, cbY + cbH + 11);
149+ const ov = ctx.getImageData(0, 0, W, H).data;
150+ // overlay pixels as the nearest fixed color
151+ const overlay = new Int16Array(W * H).fill(-1);
152+ for (let p = 0; p < W * H; p++) {
153+ if (ov[4 * p + 3] < 110) continue;
154+ let best = 0;
155+ let bd = Infinity;
156+ FIXED.forEach((c, i) => {
157+ const d = (c[0] - ov[4 * p]) ** 2 + (c[1] - ov[4 * p + 1]) ** 2 + (c[2] - ov[4 * p + 2]) ** 2;
158+ if (d < bd) {
159+ bd = d;
160+ best = i;
161+ }
162+ });
163+ overlay[p] = NCMAP + best;
164+ }
165+
166+ const frames: Uint8Array[] = [];
167+ const level = (t: number) => Math.max(0, Math.min(NCMAP - 1, Math.round(t * (NCMAP - 1))));
168+ for (let fr = 0; fr < nframes; fr++) {
169+ const ph = (2 * Math.PI * fr) / nframes;
170+ const c = Math.cos(ph);
171+ const sn = Math.sin(ph);
172+ const img = new Uint8Array(W * H).fill(BG);
173+ for (let s = 0; s < np; s++) {
174+ const left = x0p + s * (S + pad);
175+ for (let py = 0; py < S; py++) {
176+ for (let px = 0; px < S; px++) {
177+ const o = s * S * S + py * S + px;
178+ const a = re[o];
179+ let idx: number;
180+ if (Number.isNaN(a)) idx = NOFIELD;
181+ else if (mode === 'pressure') idx = level(0.5 + (0.5 * (a * c + im[o] * sn)) / f.umax);
182+ else idx = level(lv[o]);
183+ img[(y0p + py) * W + left + px] = idx;
184+ }
185+ }
186+ }
187+ for (let x = 0; x < cbW; x++) {
188+ const idx = level(x / (cbW - 1));
189+ for (let y = 0; y < cbH; y++) img[(cbY + y) * W + cbX + x] = idx;
190+ }
191+ for (let p = 0; p < W * H; p++) if (overlay[p] >= 0) img[p] = overlay[p];
192+ frames.push(img);
193+ }
194+ // one period in 2 s, as in the page
195+ return encodeGif(W, H, palette, frames, Math.max(2, Math.round(200 / nframes)));
196+}
src/view/field.tsmodified+20−13View file
@@ -1,7 +1,10 @@
11 /**
22 * The solved field as the views need it, and how a complex value becomes a
33 * 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.
4+ * on a symmetric diverging scale, or the level 20 log10(|u| / ref) in dB. The
5+ * reference of the level is one number for the room (the median of |u|) and the
6+ * free field |G| at each point for the cymbal, so that 0 dB means no effect of
7+ * the cymbal.
58 */
69 import type { Slice } from '../problem.ts';
710
@@ -9,30 +12,34 @@ export type Mode = 'pressure' | 'level';
912
1013 export interface SliceField {
1114 slice: Slice;
12- /** interleaved complex, NaN outside the room */
15+ /** interleaved complex, NaN where the field is not computed */
1316 u: Float64Array;
17+ /** reference of the level at each point; when absent, Field.ref */
18+ ref?: Float64Array;
1419 }
1520
1621 export interface Field {
1722 k: number;
18- beta: number;
23+ /** the source (the listener, for the cymbal) */
1924 x0: [number, number, number];
20- a: number;
21- p: number;
2225 ngrid: number;
23- grid: Float64Array;
2426 slices: SliceField[];
25- /** wall field at the nodes, interleaved complex */
27+ /** surface field at the nodes, interleaved complex */
2628 uWall: Float64Array;
29+ /** the level's reference on the surface is the free field |G| = 1 / (4 pi |x - x0|) rather than ref */
30+ freeFieldLevel: boolean;
2731 /** color limit for the pressure: 99th percentile of |u| more than half a wavelength from the source */
2832 umax: number;
29- /** reference for the level: the median of |u| on the same points */
33+ /** reference for the level when there is no per-point reference: the median of |u| */
3034 ref: number;
35+ /** the level scale in dB, and whether it is centered (diverging colormap) */
36+ level: { min: number; max: number; diverging: boolean; label: string };
37+ /** the 3D view: what to look at, how far, and whether surfaces facing the camera are cut away */
38+ scene: { center: [number, number, number]; radius: number; cutaway: boolean };
39+ /** the cymbal's cross section (r, z) and radius, drawn on the slices */
40+ outline?: { r: Float64Array; z: Float64Array; radius: number };
3141 }
3242
33-export const DB_MIN = -30;
34-export const DB_MAX = 15;
35-
3643 /** Color limits from the slices, as plot.py computes them. */
3744 export function colorLimits(slices: SliceField[], x0: ArrayLike<number>, k: number): { umax: number; ref: number } {
3845 const vals: number[] = [];
@@ -53,8 +60,8 @@ export function colorLimits(slices: SliceField[], x0: ArrayLike<number>, k: numb
5360 }
5461
5562 /** 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 {
63+export function colorParam(re: number, im: number, mode: Mode, c: number, s: number, umax: number, ref: number, level: Field['level']): number {
5764 if (mode === 'pressure') return 0.5 + (0.5 * (re * c + im * s)) / umax;
5865 const db = 20 * Math.log10(Math.hypot(re, im) / ref);
59- return (db - DB_MIN) / (DB_MAX - DB_MIN);
66+ return (db - level.min) / (level.max - level.min);
6067 }
src/view/gif.tsadded+122−0View file
@@ -0,0 +1,122 @@
1+/**
2+ * A minimal animated GIF (GIF89a) encoder: one global palette of 256 colors,
3+ * frames of palette indices, LZW compression, and an infinite loop.
4+ */
5+
6+/** Encode frames (width x height palette indices each) with palette (256 x RGB) and a delay per frame in 1/100 s. */
7+export function encodeGif(width: number, height: number, palette: Uint8Array, frames: Uint8Array[], delayCs: number): Uint8Array {
8+ const out = new ByteWriter();
9+ out.str('GIF89a');
10+ out.u16(width);
11+ out.u16(height);
12+ out.byte(0xf7); // global color table, 8 bits per color, 256 entries
13+ out.byte(0); // background color index
14+ out.byte(0); // pixel aspect ratio
15+ for (let i = 0; i < 768; i++) out.byte(palette[i] ?? 0);
16+ if (frames.length > 1) {
17+ // NETSCAPE2.0 application extension: loop forever
18+ out.bytes([0x21, 0xff, 0x0b]);
19+ out.str('NETSCAPE2.0');
20+ out.bytes([0x03, 0x01, 0x00, 0x00, 0x00]);
21+ }
22+ for (const f of frames) {
23+ out.bytes([0x21, 0xf9, 0x04, 0x00]); // graphic control extension, no transparency
24+ out.u16(delayCs);
25+ out.bytes([0x00, 0x00]);
26+ out.byte(0x2c); // image descriptor
27+ out.u16(0);
28+ out.u16(0);
29+ out.u16(width);
30+ out.u16(height);
31+ out.byte(0); // no local color table, not interlaced
32+ lzw(f, 8, out);
33+ }
34+ out.byte(0x3b);
35+ return out.result();
36+}
37+
38+class ByteWriter {
39+ private buf = new Uint8Array(1 << 16);
40+ private n = 0;
41+
42+ byte(b: number) {
43+ if (this.n === this.buf.length) {
44+ const nb = new Uint8Array(2 * this.buf.length);
45+ nb.set(this.buf);
46+ this.buf = nb;
47+ }
48+ this.buf[this.n++] = b;
49+ }
50+ bytes(bs: number[]) {
51+ for (const b of bs) this.byte(b);
52+ }
53+ u16(v: number) {
54+ this.byte(v & 0xff);
55+ this.byte((v >> 8) & 0xff);
56+ }
57+ str(s: string) {
58+ for (let i = 0; i < s.length; i++) this.byte(s.charCodeAt(i));
59+ }
60+ result(): Uint8Array {
61+ return this.buf.slice(0, this.n);
62+ }
63+}
64+
65+/** LZW-compress indices into GIF image data sub-blocks. */
66+function lzw(indices: Uint8Array, minCodeSize: number, out: ByteWriter) {
67+ out.byte(minCodeSize);
68+ const clear = 1 << minCodeSize;
69+ const eoi = clear + 1;
70+ let codeSize = minCodeSize + 1;
71+ let next = eoi + 1;
72+ // dictionary: (prefix code << 8 | byte) -> code
73+ const dict = new Map<number, number>();
74+ const block: number[] = [];
75+ let acc = 0;
76+ let nbits = 0;
77+ const emit = (code: number) => {
78+ acc |= code << nbits;
79+ nbits += codeSize;
80+ while (nbits >= 8) {
81+ block.push(acc & 0xff);
82+ acc >>>= 8;
83+ nbits -= 8;
84+ if (block.length === 255) {
85+ out.byte(255);
86+ out.bytes(block);
87+ block.length = 0;
88+ }
89+ }
90+ };
91+ emit(clear);
92+ let prefix = indices[0];
93+ for (let i = 1; i < indices.length; i++) {
94+ const c = indices[i];
95+ const key = (prefix << 8) | c;
96+ const code = dict.get(key);
97+ if (code !== undefined) {
98+ prefix = code;
99+ continue;
100+ }
101+ emit(prefix);
102+ if (next < 4096) {
103+ dict.set(key, next++);
104+ // the decoder widens its codes one code later than the encoder adds the entry
105+ if (next > 1 << codeSize && codeSize < 12) codeSize++;
106+ } else {
107+ emit(clear);
108+ dict.clear();
109+ codeSize = minCodeSize + 1;
110+ next = eoi + 1;
111+ }
112+ prefix = c;
113+ }
114+ emit(prefix);
115+ emit(eoi);
116+ if (nbits > 0) block.push(acc & 0xff);
117+ if (block.length) {
118+ out.byte(block.length);
119+ out.bytes(block);
120+ }
121+ out.byte(0);
122+}
src/view/room3d.tsmodified+67−42View file
@@ -1,13 +1,13 @@
11 /**
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.
2+ * A WebGPU view of the scene: the surface colored by the surface field, the
3+ * three planes, and the source. For the room, walls whose outward normal faces
4+ * the camera are cut away in the fragment shader, so the view looks into the
5+ * room through the nearer walls; the cymbal is drawn whole. Drag to orbit,
6+ * scroll or pinch to zoom.
77 */
88 import { basis2d, type Surface } from '../surface.ts';
99 import { DIVERGING, SEQUENTIAL } from './colormaps.ts';
10-import { DB_MAX, DB_MIN, type Field, type Mode } from './field.ts';
10+import type { Field, Mode } from './field.ts';
1111
1212 /** Multisampling; 1 when the query string has ?nomsaa (SwiftShader's GPU process does not survive the resolve in headless Chrome). */
1313 const SAMPLES = typeof location !== 'undefined' && location.search.includes('nomsaa') ? 1 : 4;
@@ -16,21 +16,24 @@ const SHARED = /* wgsl */ `
1616 struct View {
1717 viewProj: mat4x4f,
1818 eye: vec4f,
19- p1: vec4f, // cos wt, sin wt, umax, ref
20- p2: vec4f, // mode (0 pressure, 1 level), db min, db max, 0
19+ p1: vec4f, // cos wt, sin wt, umax, 0
20+ p2: vec4f, // mode (0 pressure, 1 level), db min, db max, level colormap (0 sequential, 1 diverging)
21+ p3: vec4f, // cut away surfaces facing the camera (0, 1), 0, 0, 0
2122 }
2223 @group(0) @binding(0) var<uniform> V: View;
2324 @group(0) @binding(1) var cmapD: texture_2d<f32>;
2425 @group(0) @binding(2) var cmapS: texture_2d<f32>;
2526 @group(0) @binding(3) var samp: sampler;
2627
27-fn colorOf(re: f32, im: f32) -> vec3f {
28+// lref: the reference of the level at this point ('ref' is reserved in WGSL)
29+fn colorOf(re: f32, im: f32, lref: f32) -> vec3f {
2830 if (V.p2.x < 0.5) {
2931 let t = clamp(0.5 + 0.5 * (re * V.p1.x + im * V.p1.y) / V.p1.z, 0.0, 1.0);
3032 return textureSampleLevel(cmapD, samp, vec2f(t, 0.5), 0.0).rgb;
3133 }
32- let db = 20.0 * log2(max(length(vec2f(re, im)) / V.p1.w, 1e-30)) * 0.30102999566;
34+ let db = 20.0 * log2(max(length(vec2f(re, im)) / lref, 1e-30)) * 0.30102999566;
3335 let t = clamp((db - V.p2.y) / (V.p2.z - V.p2.y), 0.0, 1.0);
36+ if (V.p2.w > 0.5) { return textureSampleLevel(cmapD, samp, vec2f(t, 0.5), 0.0).rgb; }
3437 return textureSampleLevel(cmapS, samp, vec2f(t, 0.5), 0.0).rgb;
3538 }
3639 `;
@@ -42,21 +45,23 @@ struct VOut {
4245 @location(0) world: vec3f,
4346 @location(1) nrm: vec3f,
4447 @location(2) u: vec2f,
48+ @location(3) lref: f32,
4549 }
46-@vertex fn vs(@location(0) p: vec3f, @location(1) n: vec3f, @location(2) u: vec2f) -> VOut {
50+@vertex fn vs(@location(0) p: vec3f, @location(1) n: vec3f, @location(2) u: vec2f, @location(3) lref: f32) -> VOut {
4751 var o: VOut;
4852 o.pos = V.viewProj * vec4f(p, 1.0);
4953 o.world = p;
5054 o.nrm = n;
5155 o.u = u;
56+ o.lref = lref;
5257 return o;
5358 }
5459 @fragment fn fs(i: VOut) -> @location(0) vec4f {
5560 let toEye = normalize(V.eye.xyz - i.world);
5661 let n = normalize(i.nrm);
57- if (dot(n, toEye) > 0.0) { discard; }
62+ if (V.p3.x > 0.5 && dot(n, toEye) > 0.0) { discard; }
5863 let shade = 0.6 + 0.4 * abs(dot(n, toEye));
59- return vec4f(colorOf(i.u.x, i.u.y) * shade, 1.0);
64+ return vec4f(colorOf(i.u.x, i.u.y, i.lref) * shade, 1.0);
6065 }
6166 `;
6267
@@ -80,14 +85,15 @@ struct VOut {
8085 let x0 = vec2i(floor(x));
8186 let x1 = min(x0 + 1, vec2i(dim) - 1);
8287 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+ // (Re u, Im u, level reference)
89+ let a = textureLoad(field, x0, 0).xyz;
90+ let b = textureLoad(field, vec2i(x1.x, x0.y), 0).xyz;
91+ let c = textureLoad(field, vec2i(x0.x, x1.y), 0).xyz;
92+ let d = textureLoad(field, x1, 0).xyz;
93+ // NaN marks points where the field is not computed
8894 if (any(a != a) || any(b != b) || any(c != c) || any(d != d)) { discard; }
8995 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);
96+ return vec4f(colorOf(u.x, u.y, u.z), 1.0);
9197 }
9298 `;
9399
@@ -176,7 +182,7 @@ export class RoomView {
176182 this.ctx = canvas.getContext('webgpu') as GPUCanvasContext;
177183 this.format = navigator.gpu.getPreferredCanvasFormat();
178184 this.ctx.configure({ device, format: this.format, alphaMode: 'opaque' });
179- this.ubuf = device.createBuffer({ size: 112, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
185+ this.ubuf = device.createBuffer({ size: 128, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
180186 const mkmap = (lut: Uint8Array<ArrayBuffer>) => {
181187 const t = device.createTexture({ size: [256, 1], format: 'rgba8unorm', usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST });
182188 device.queue.writeTexture({ texture: t }, lut, { bytesPerRow: 1024 }, [256, 1]);
@@ -210,11 +216,12 @@ export class RoomView {
210216 };
211217 this.wallPipe = common(WALL, baseLayout, [
212218 {
213- arrayStride: 32,
219+ arrayStride: 36,
214220 attributes: [
215221 { shaderLocation: 0, offset: 0, format: 'float32x3' },
216222 { shaderLocation: 1, offset: 12, format: 'float32x3' },
217223 { shaderLocation: 2, offset: 24, format: 'float32x2' },
224+ { shaderLocation: 3, offset: 32, format: 'float32' },
218225 ],
219226 },
220227 ]);
@@ -291,8 +298,12 @@ export class RoomView {
291298 );
292299 }
293300
301+ /** the scene's radius: 1 for the room */
302+ private scale = 1;
303+ private center: [number, number, number] = [0, 0, 0];
304+
294305 private zoom(f: number) {
295- this.dist = Math.min(12, Math.max(2.2, this.dist * f));
306+ this.dist = Math.min(12 * this.scale, Math.max(2.2 * this.scale, this.dist * f));
296307 this.dirty = true;
297308 }
298309
@@ -312,10 +323,15 @@ export class RoomView {
312323 this.dirty = true;
313324 }
314325
315- /** The wall mesh: every patch sampled on an nv x nv grid, with the wall field interpolated by the patch's Legendre expansion. */
326+ /** The surface mesh: every patch sampled on an nv x nv grid, with the surface field interpolated by the patch's Legendre expansion. */
316327 setField(surf: Surface, f: Field, nv = 9) {
317328 const { device } = this;
318329 this.field = f;
330+ if (f.scene.radius !== this.scale) {
331+ this.dist *= f.scene.radius / this.scale;
332+ this.scale = f.scene.radius;
333+ }
334+ this.center = f.scene.center;
319335 const npols = surf.npols;
320336 const npatches = surf.npatches;
321337 const npts = surf.npts;
@@ -326,7 +342,7 @@ export class RoomView {
326342 }
327343 const B = basis2d(uv, surf.norder); // (npols, nv^2)
328344 const ng = nv * nv;
329- const verts = new Float32Array(8 * ng * npatches);
345+ const verts = new Float32Array(9 * ng * npatches);
330346 const idx = new Uint32Array(6 * (nv - 1) * (nv - 1) * npatches);
331347 const cr = new Float64Array(npols);
332348 const ci = new Float64Array(npols);
@@ -345,7 +361,7 @@ export class RoomView {
345361 ci[n] = si;
346362 }
347363 for (let q = 0; q < ng; q++) {
348- const o = 8 * (j * ng + q);
364+ const o = 9 * (j * ng + q);
349365 const ru = [g[3 * ng + q], g[4 * ng + q], g[5 * ng + q]];
350366 const rv = [g[6 * ng + q], g[7 * ng + q], g[8 * ng + q]];
351367 const nn = norm(cross(ru, rv));
@@ -355,7 +371,8 @@ export class RoomView {
355371 ur += cr[n] * B[n * ng + q];
356372 ui += ci[n] * B[n * ng + q];
357373 }
358- verts.set([g[q], g[ng + q], g[2 * ng + q], nn[0], nn[1], nn[2], ur, ui], o);
374+ const ref = f.freeFieldLevel ? 1 / (4 * Math.PI * Math.hypot(g[q] - f.x0[0], g[ng + q] - f.x0[1], g[2 * ng + q] - f.x0[2])) : f.ref;
375+ verts.set([g[q], g[ng + q], g[2 * ng + q], nn[0], nn[1], nn[2], ur, ui, ref], o);
359376 }
360377 for (let a = 0; a < nv - 1; a++) for (let b = 0; b < nv - 1; b++) {
361378 const v0 = j * ng + a * nv + b;
@@ -372,26 +389,28 @@ export class RoomView {
372389 device.queue.writeBuffer(ibuf, 0, idx);
373390 this.wall = { vbuf, ibuf, count: idx.length };
374391
375- // the three planes, as textures of (Re u, Im u) with NaN outside the room
392+ // the three planes, as textures of (Re u, Im u, level reference) with NaN where there is no field
376393 for (const p of this.planes) {
377394 p.vbuf.destroy();
378395 p.tex.destroy();
379396 }
380- this.planes = f.slices.map(({ slice, u }) => {
397+ this.planes = f.slices.map(({ slice, u, ref }) => {
381398 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);
399+ const tex = device.createTexture({ size: [n, n], format: 'rgba32float', usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST });
400+ const data = new Float32Array(4 * n * n);
384401 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;
402+ data[4 * q] = slice.mask[q] ? u[2 * q] : NaN;
403+ data[4 * q + 1] = slice.mask[q] ? u[2 * q + 1] : NaN;
404+ data[4 * q + 2] = ref ? ref[q] : f.ref;
387405 }
388- device.queue.writeTexture({ texture: tex }, data, { bytesPerRow: 8 * n }, [n, n]);
406+ device.queue.writeTexture({ texture: tex }, data, { bytesPerRow: 16 * n }, [n, n]);
389407 const [i, jx] = slice.axes;
408+ const [ea, eb] = slice.extent;
390409 const corner = (s: number, t: number) => {
391410 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];
411+ x[i] = ea[0] + (ea[1] - ea[0]) * s;
412+ x[jx] = eb[0] + (eb[1] - eb[0]) * t;
413+ x[slice.fixed] = slice.value;
395414 return [...x, s, t];
396415 };
397416 const quad = new Float32Array([...corner(0, 0), ...corner(1, 0), ...corner(1, 1), ...corner(0, 0), ...corner(1, 1), ...corner(0, 1)]);
@@ -402,7 +421,7 @@ export class RoomView {
402421 });
403422
404423 // the source: a small octahedron
405- const r = 0.035;
424+ const r = 0.035 * this.scale;
406425 const [x, y, z] = f.x0;
407426 const P = [
408427 [x + r, y, z], [x - r, y, z], [x, y + r, z], [x, y - r, z], [x, y, z + r], [x, y, z - r],
@@ -436,14 +455,20 @@ export class RoomView {
436455 this.dirty = false;
437456 const { device } = this;
438457 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]));
458+ const c = this.center;
459+ const eye = [
460+ c[0] + this.dist * Math.sin(this.phi) * Math.cos(this.theta),
461+ c[1] + this.dist * Math.sin(this.phi) * Math.sin(this.theta),
462+ c[2] + this.dist * Math.cos(this.phi),
463+ ];
464+ const vp = mul(perspective((35 * Math.PI) / 180, w / h, 0.01 * this.scale, 50 * this.scale), lookAt(eye, c, [0, 0, 1]));
441465 const f = this.field;
442- const u = new Float32Array(28);
466+ const u = new Float32Array(32);
443467 u.set(vp, 0);
444468 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);
469+ u.set([Math.cos(phase), Math.sin(phase), f?.umax ?? 1, 0], 20);
470+ u.set([mode === 'pressure' ? 0 : 1, f?.level.min ?? -30, f?.level.max ?? 15, f?.level.diverging ? 1 : 0], 24);
471+ u.set([f?.scene.cutaway === false ? 0 : 1, 0, 0, 0], 28);
447472 device.queue.writeBuffer(this.ubuf, 0, u);
448473 const enc = device.createCommandEncoder();
449474 const pass = enc.beginRenderPass({
src/view/slices.tsmodified+57−18View file
@@ -1,10 +1,11 @@
11 /**
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).
2+ * The three planes as 2D images, redrawn every frame so that the pressure can
3+ * be animated over one period. Clicking a plane moves the source (the listener,
4+ * for the cymbal) to that point, whose third coordinate is the plane's. For the
5+ * cymbal, its outline is drawn over the planes.
56 */
67 import { DIVERGING, SEQUENTIAL } from './colormaps.ts';
7-import { colorParam, type Field, type Mode } from './field.ts';
8+import { colorParam, type Field, type Mode, type SliceField } from './field.ts';
89
910 const AXIS = ['x', 'y', 'z'];
1011
@@ -52,7 +53,7 @@ export class SliceView {
5253 this.img = f.slices.map(() => new ImageData(f.ngrid, f.ngrid));
5354 f.slices.forEach(({ slice }, s) => {
5455 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+ this.titles[s].textContent = `${slice.name}-plane at ${AXIS[slice.fixed]} = ${slice.value.toFixed(2)} (${AXIS[i]} across, ${AXIS[j]} up)`;
5657 });
5758 }
5859
@@ -65,9 +66,11 @@ export class SliceView {
6566 const fx = (ev.clientX - r.left) / r.width;
6667 const fy = 1 - (ev.clientY - r.top) / r.height;
6768 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;
69+ const [ea, eb] = slice.extent;
70+ const x: [number, number, number] = [0, 0, 0];
71+ x[slice.fixed] = slice.value;
72+ x[slice.axes[0]] = ea[0] + (ea[1] - ea[0]) * fx;
73+ x[slice.axes[1]] = eb[0] + (eb[1] - eb[0]) * fy;
7174 const col = Math.round(fx * (f.ngrid - 1));
7275 const row = Math.round(fy * (f.ngrid - 1));
7376 if (col < 0 || row < 0 || col >= f.ngrid || row >= f.ngrid) return null;
@@ -83,15 +86,16 @@ export class SliceView {
8386 const p = this.at(s, ev);
8487 const f = this.field;
8588 if (!p || !f) return;
86- const { u, slice } = f.slices[s];
89+ const { u, slice, ref } = f.slices[s];
90+ const at = `(${p.x.map((v) => v.toFixed(2)).join(', ')})`;
8791 if (!slice.mask[p.q]) {
88- this.onHover(`(${p.x.map((v) => v.toFixed(2)).join(', ')}): outside the room`);
92+ this.onHover(`${at}: ${f.outline ? 'inside the cymbal' : 'outside the room'}`);
8993 return;
9094 }
9195 const re = u[2 * p.q];
9296 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`);
97+ const db = 20 * Math.log10(Math.hypot(re, im) / (ref ? ref[p.q] : f.ref));
98+ this.onHover(`${at}: u = ${re.toExponential(3)} ${im < 0 ? '−' : '+'} ${Math.abs(im).toExponential(3)}i, level ${db.toFixed(1)} dB ${f.level.label} · click to move the ${f.outline ? 'listener' : 'source'} here`);
9599 }
96100
97101 draw(mode: Mode, phase: number) {
@@ -99,8 +103,8 @@ export class SliceView {
99103 if (!f) return;
100104 const c = Math.cos(phase);
101105 const sn = Math.sin(phase);
102- const map = mode === 'pressure' ? DIVERGING : SEQUENTIAL;
103- f.slices.forEach(({ slice, u }, s) => {
106+ const map = mode === 'pressure' || f.level.diverging ? DIVERGING : SEQUENTIAL;
107+ f.slices.forEach(({ slice, u, ref }, s) => {
104108 const img = this.img[s];
105109 const n = f.ngrid;
106110 for (let row = 0; row < n; row++) {
@@ -111,7 +115,7 @@ export class SliceView {
111115 img.data[o + 3] = 0;
112116 continue;
113117 }
114- const t = colorParam(u[2 * q], u[2 * q + 1], mode, c, sn, f.umax, f.ref);
118+ const t = colorParam(u[2 * q], u[2 * q + 1], mode, c, sn, f.umax, ref ? ref[q] : f.ref, f.level);
115119 const i = Math.max(0, Math.min(255, Math.round(t * 255))) * 4;
116120 img.data[o] = map[i];
117121 img.data[o + 1] = map[i + 1];
@@ -127,10 +131,15 @@ export class SliceView {
127131 ctx.imageSmoothingEnabled = true;
128132 ctx.imageSmoothingQuality = 'high';
129133 ctx.drawImage(tmp, 0, 0, cv.width, cv.height);
130- // the source
131134 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;
135+ const [ea, eb] = slice.extent;
136+ const X = (a: number) => ((a - ea[0]) / (ea[1] - ea[0])) * cv.width;
137+ const Y = (b: number) => (1 - (b - eb[0]) / (eb[1] - eb[0])) * cv.height;
138+ if (f.outline) drawOutline(ctx, slice, f.outline, X, Y);
139+ // the source, when it lies on this plane
140+ if (Math.abs(f.x0[slice.fixed] - slice.value) > 1e-9) return;
141+ const px = X(f.x0[i]);
142+ const py = Y(f.x0[j]);
134143 ctx.beginPath();
135144 ctx.arc(px, py, 5, 0, 2 * Math.PI);
136145 ctx.lineWidth = 1.5;
@@ -144,3 +153,33 @@ export class SliceView {
144153 });
145154 }
146155 }
156+
157+/** The cymbal: its cross section on a vertical plane, the rim (dashed) on a horizontal one. */
158+export function drawOutline(ctx: CanvasRenderingContext2D, slice: SliceField['slice'], o: NonNullable<Field['outline']>, X: (a: number) => number, Y: (b: number) => number) {
159+ ctx.save();
160+ ctx.strokeStyle = '#222';
161+ ctx.fillStyle = 'rgba(40, 40, 40, 0.85)';
162+ ctx.lineWidth = 1;
163+ if (slice.fixed === 2) {
164+ ctx.setLineDash([4, 3]);
165+ ctx.beginPath();
166+ ctx.arc(X(0), Y(0), Math.abs(X(o.radius) - X(0)), 0, 2 * Math.PI);
167+ ctx.stroke();
168+ } else if (Math.abs(slice.value) < o.radius) {
169+ // the plane is at distance |value| from the axis: horizontal coordinate +-sqrt(r^2 - value^2)
170+ for (const sign of [1, -1]) {
171+ ctx.beginPath();
172+ let started = false;
173+ for (let q = 0; q < o.r.length; q++) {
174+ if (o.r[q] < Math.abs(slice.value)) continue;
175+ const a = sign * Math.sqrt(o.r[q] ** 2 - slice.value ** 2);
176+ if (started) ctx.lineTo(X(a), Y(o.z[q]));
177+ else ctx.moveTo(X(a), Y(o.z[q]));
178+ started = true;
179+ }
180+ ctx.closePath();
181+ ctx.fill();
182+ }
183+ }
184+ ctx.restore();
185+}
test/test-page.tsmodified+164−34View file
@@ -8,12 +8,25 @@
88 * 3. the self-test problem: the field of a source outside the room, solved
99 * with the impedance equation and compared with the exact field inside.
1010 *
11- * The surface is kept small so that the software WebGPU adapter can run it.
11+ * The same three checks then run on a coarse cymbal (178 patches), which uses
12+ * the interpolated geometry and the exterior equation: the near weights of a face
13+ * patch and a rim patch, Green's identity for a source outside the cymbal, and
14+ * the exterior self-test with sources spread inside the cymbal (?cymbal=0 skips it,
15+ * ?room=0 skips the rounded cube).
16+ *
17+ * ?bench=small|medium|large instead runs a named cymbal case of the CUDA proof of
18+ * concept (same surface, same unknowns) and reports its self-test error, the
19+ * scattered field at the listener, the reciprocity check and the timings, for a
20+ * direct comparison with run.py --problem cymbal.
21+ *
22+ * The surfaces are kept small so that the software WebGPU adapter can run them.
1223 */
13-import { roundedCube } from '../src/surface.ts';
24+import { roundedCube, type Surface } from '../src/surface.ts';
25+import { cymbal, CymbalShape, CYMBAL_DEFAULTS } from '../src/cymbal.ts';
26+import { selftestSources, selftestPoints as cymbalPoints, sourceSum, LISTENER, C_SOUND } from '../src/cymbalProblem.ts';
1427 import { green, greenDn, selftestRhs, selftestPoints, relError, ALPHA, type CVec } from '../src/problem.ts';
1528 import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
16-import { requestSolverDevice } from '../src/gpu/device.ts';
29+import { adapterDescription, requestSolverDevice } from '../src/gpu/device.ts';
1730 import { Operator } from '../src/gpu/operator.ts';
1831
1932 declare global {
@@ -40,37 +53,15 @@ const check = (name: string, val: number, tol: number) => {
4053 log(`${ok ? 'ok ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`, ok ? 'ok' : 'fail');
4154 };
4255
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);
56+/** Near weights against the brute-force reference for three nodes of each patch in tps. */
57+async function checkWeights(op: Operator, surf: Surface, tps: number[], k: number, label: string) {
5858 const npts = surf.npts;
5959 const npols = surf.npols;
60- log(`surface: p = ${p}, ${surf.npatches} patches of order ${norder}, ${npts} nodes; k = ${k}`);
6160 const xs = surf.srcvals.subarray(0, 3 * npts);
6261 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);
62+ const { pairs, W } = await op.readWallWeights();
63+ const worst = new Float64Array(6);
64+ for (const tp of tps) {
7465 for (const loc of [0, Math.floor(npols / 2), npols - 2]) {
7566 const i = tp * npols + loc;
7667 const xt = [xs[i], xs[npts + i], xs[2 * npts + i]];
@@ -93,12 +84,58 @@ async function main() {
9384 worst[mm] = Math.max(worst[mm], e);
9485 errs.push(e.toExponential(1));
9586 }
87+ const kind = ['plain', 'preimage', 'self'][pairs.pairF[8 * pp + 2]];
9688 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(' ')}`);
89+ log(` node ${i} patch ${j} (${kind}, dist/rad ${d.toFixed(2)}): rel err per kernel ${errs.join(' ')}`);
9890 }
9991 }
100- for (let mm = 0; mm < 6; mm++) check(`near weights, kernel ${mm}, worst relative error`, worst[mm], 1e-4);
10192 }
93+ for (let mm = 0; mm < 6; mm++) check(`${label}near weights, kernel ${mm}, worst relative error`, worst[mm], 1e-4);
94+}
95+
96+async function main() {
97+ const t0 = performance.now();
98+ if (logEl.textContent === 'starting…') logEl.textContent = '';
99+ const gpu = await requestSolverDevice();
100+ const device = gpu.device;
101+ const lim = device.limits;
102+ log(`device: ${adapterDescription(gpu.adapter ?? null)}; maxStorageBufferBindingSize ${(lim.maxStorageBufferBindingSize / 2 ** 20).toFixed(0)} MB, workgroup storage ${lim.maxComputeWorkgroupStorageSize} B`);
103+
104+ const params = new URLSearchParams(location.search);
105+ const bench = params.get('bench');
106+ if (bench) {
107+ await cymbalBench(gpu, bench);
108+ log(`done, ${((performance.now() - t0) / 1000).toFixed(1)} s`);
109+ window.__RESULTS__ = { passed, failed, lines };
110+ return;
111+ }
112+ if (params.get('room') !== '0') await roomTests(gpu, params);
113+ if (params.get('cymbal') !== '0') await cymbalTests(gpu);
114+ log(`${passed} passed, ${failed} failed, ${((performance.now() - t0) / 1000).toFixed(1)} s`, failed ? 'fail' : 'ok');
115+ window.__RESULTS__ = { passed, failed, lines };
116+}
117+
118+/** The rounded cube: exact geometry, interior equation. */
119+async function roomTests(gpu: Awaited<ReturnType<typeof requestSolverDevice>>, params: URLSearchParams) {
120+ const m = Number(params.get('m') ?? 3);
121+ const norder = Number(params.get('norder') ?? 6);
122+ const k = Number(params.get('k') ?? 4);
123+ const a = 1;
124+ const p = Number(params.get('p') ?? 8);
125+ const surf = roundedCube(a, p, m, norder);
126+ const npts = surf.npts;
127+ const npols = surf.npols;
128+ log(`surface: p = ${p}, ${surf.npatches} patches of order ${norder}, ${npts} nodes; k = ${k}`);
129+ const xs = surf.srcvals.subarray(0, 3 * npts);
130+ const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
131+
132+ const lam = -0.2;
133+ const op = new Operator(gpu, surf, k, ALPHA, lam, 0, { progress: () => {} });
134+ const info = await op.setup();
135+ log(`operator: ${info.npairs} near pairs, ${info.nleaves} leaves, ${info.noversampled} oversampled nodes, weights ${(info.weightsBytes / 2 ** 20).toFixed(1)} MB, setup ${info.setupSeconds.toFixed(2)} s`);
136+
137+ // 1. near weights against the reference, for three nodes of one patch
138+ await checkWeights(op, surf, [4], k, '');
102139
103140 // 2. Green's identity through the operator
104141 {
@@ -135,8 +172,6 @@ async function main() {
135172 }
136173
137174 op.destroy();
138- log(`${passed} passed, ${failed} failed, ${((performance.now() - t0) / 1000).toFixed(1)} s`, failed ? 'fail' : 'ok');
139- window.__RESULTS__ = { passed, failed, lines };
140175 }
141176
142177 main().catch((e) => {
@@ -144,3 +179,98 @@ main().catch((e) => {
144179 log(`FATAL ${msg}`, 'fail');
145180 window.__RESULTS__ = { passed, failed: failed + 1, fatal: msg, lines };
146181 });
182+
183+/** The same checks on a coarse cymbal: interpolated geometry, exterior equation, lambda = 0. */
184+async function cymbalTests(gpu: Awaited<ReturnType<typeof requestSolverDevice>>) {
185+ const shape = new CymbalShape({ ...CYMBAL_DEFAULTS, tRim: 0.02 });
186+ const surf = cymbal(shape, 0.1, 5);
187+ const npts = surf.npts;
188+ const k = (2 * Math.PI * 1000) / 343;
189+ log(`cymbal surface: rim ${shape.par.tRim * 1000} mm, ${surf.npatches} patches of order ${surf.norder}, ${npts} nodes; k = ${k.toFixed(2)}`);
190+ const xs = surf.srcvals.subarray(0, 3 * npts);
191+ const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
192+ const op = new Operator(gpu, surf, k, ALPHA, 0, 1, { progress: () => {} });
193+ const info = await op.setup();
194+ log(`operator: ${info.npairs} near pairs, ${info.nleaves} leaves, weights ${(info.weightsBytes / 2 ** 20).toFixed(1)} MB, setup ${info.setupSeconds.toFixed(2)} s`);
195+
196+ // a patch on the top face near the bell, and the first rim patch
197+ const nface = (surf.npatches - surf.layout.nrim) / 2;
198+ await checkWeights(op, surf, [surf.layout.mC ** 2 + 2, 2 * nface], k, 'cymbal: ');
199+
200+ // Green's identity for a source outside the body: S[du/dn] - D[u] = u/2 inside it
201+ {
202+ const xe = [0.05, -0.1, 0.3];
203+ const u = green(xs, xe, k);
204+ const un = greenDn(xs, nrm, xe, k);
205+ const S = await op.applySingle(4, un);
206+ const D = await op.applySingle(5, u);
207+ let num = 0;
208+ let den = 0;
209+ for (let i = 0; i < npts; i++) {
210+ const rr = S[2 * i] - D[2 * i] - u[2 * i] / 2;
211+ const ri = S[2 * i + 1] - D[2 * i + 1] - u[2 * i + 1] / 2;
212+ num += (rr * rr + ri * ri) * surf.wts[i];
213+ den += (u[2 * i] ** 2 + u[2 * i + 1] ** 2) * surf.wts[i];
214+ }
215+ check("cymbal: Green's identity on the surface (relative, weighted)", Math.sqrt(num / den), 1e-3);
216+ }
217+
218+ // the exterior self-test: sources spread inside the cymbal
219+ {
220+ const { x0s, ws } = selftestSources(shape);
221+ const rhs = sourceSum(xs, x0s, ws, k, nrm);
222+ const sol = await op.solve(rhs, 1e-7, 100);
223+ log(` GMRES: ${sol.niter} iterations, residual ${sol.rres.toExponential(2)}, ${sol.seconds.toFixed(2)} s`);
224+ const xt = cymbalPoints(shape, 30);
225+ const ut = await op.eval(xt, sol.x, sol.phi1);
226+ check('cymbal: self-test, relative error at 30 points around the cymbal', relError(ut, sourceSum(xt, x0s, ws, k)), 2e-3);
227+ const ub = await op.evalOnSurface(sol.x, sol.phi1);
228+ check('cymbal: self-test, relative error on the surface', relError(ub, sourceSum(xs, x0s, ws, k)), 5e-3);
229+ }
230+ op.destroy();
231+}
232+
233+/** A named cymbal case of the CUDA proof of concept, for comparison with run.py. */
234+async function cymbalBench(gpu: Awaited<ReturnType<typeof requestSolverDevice>>, name: string) {
235+ const cases: Record<string, { freq: number; h: number; norder: number }> = {
236+ small: { freq: 1000, h: 0.06, norder: 6 },
237+ medium: { freq: 3000, h: 0.04, norder: 8 },
238+ large: { freq: 8000, h: 0.03, norder: 8 },
239+ };
240+ const c = cases[name];
241+ if (!c) throw new Error(`no case ${name}`);
242+ const shape = new CymbalShape();
243+ const surf = cymbal(shape, c.h, c.norder);
244+ const npts = surf.npts;
245+ const k = (2 * Math.PI * c.freq) / C_SOUND;
246+ log(`cymbal ${name}: ${surf.npatches} patches of order ${c.norder}, ${npts} points, k = ${k.toFixed(2)}`);
247+ const xs = surf.srcvals.subarray(0, 3 * npts);
248+ const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
249+ const op = new Operator(gpu, surf, k, ALPHA, 0, 1, { progress: () => {} });
250+ const info = await op.setup();
251+ log(`setup (near quadrature) ${info.setupSeconds.toFixed(2)} s: ${info.npairs} near pairs, ${info.nleaves} leaves, weights ${(info.weightsBytes / 2 ** 20).toFixed(0)} MB`);
252+ const eps = 1e-6;
253+ {
254+ const { x0s, ws } = selftestSources(shape);
255+ const sol = await op.solve(sourceSum(xs, x0s, ws, k, nrm), eps, 400);
256+ const xt = cymbalPoints(shape, 50);
257+ const e1 = relError(await op.eval(xt, sol.x, sol.phi1), sourceSum(xt, x0s, ws, k));
258+ const e2 = relError(await op.evalOnSurface(sol.x, sol.phi1), sourceSum(xs, x0s, ws, k));
259+ log(`self-test: ${sol.niter} iterations, ${sol.seconds.toFixed(2)} s; relative error ${e1.toExponential(2)} off the surface, ${e2.toExponential(2)} on it`);
260+ }
261+ const B = [-0.15, 0.05, -0.12];
262+ const rhs = (x0: number[]) => {
263+ const g = greenDn(xs, nrm, x0, k);
264+ for (let i = 0; i < g.length; i++) g[i] = -g[i];
265+ return g;
266+ };
267+ const solL = await op.solve(rhs(LISTENER), eps, 400);
268+ log(`main problem: ${solL.niter} iterations, residual ${solL.rres.toExponential(2)}, solve ${solL.seconds.toFixed(2)} s`);
269+ const vL = await op.eval(Float64Array.from([LISTENER[0], B[0], LISTENER[1], B[1], LISTENER[2], B[2]]), solL.x, solL.phi1);
270+ log(`scattered field at the listener ${vL[0].toExponential(6)} ${vL[1] < 0 ? '-' : '+'} ${Math.abs(vL[1]).toExponential(6)}i`);
271+ const solB = await op.solve(rhs(B), eps, 400);
272+ const vB = await op.eval(Float64Array.from(LISTENER), solB.x, solB.phi1);
273+ const rec = Math.hypot(vL[2] - vB[0], vL[3] - vB[1]) / Math.hypot(vL[2], vL[3]);
274+ log(`reciprocity: at B due to the listener ${vL[2].toExponential(6)} ${vL[3] < 0 ? '-' : '+'} ${Math.abs(vL[3]).toExponential(6)}i, relative difference ${rec.toExponential(2)}`);
275+ op.destroy();
276+}