3D acoustic scattering solved live in the browser on WebGPU
A leapfrog for the acoustic wave equation on a cubic grid, one compute
dispatch per timestep, drawn by ray-marching the pressure buffer the
solver just wrote. Scenes (sphere, two spheres, aperture, random medium)
are built on the CPU; the open boundary is a sponge layer. A microphone
records one sample per timestep on the GPU and plays back through Web
Audio, resampled to 48 kHz. Sibling of acoustic-scattering-2d, without
the numbl/MATLAB layer for now.
18 changed files+4586−0
.github/workflows/deploy.ymladded+41−0View file
@@ -0,0 +1,41 @@
1+name: deploy
2+on:
3+ push:
4+ branches: [main]
5+ workflow_dispatch:
6+
7+permissions:
8+ contents: read
9+ pages: write
10+ id-token: write
11+
12+concurrency:
13+ group: pages
14+ cancel-in-progress: true
15+
16+jobs:
17+ build-deploy:
18+ runs-on: ubuntu-latest
19+ environment:
20+ name: github-pages
21+ url: ${{ steps.deployment.outputs.page_url }}
22+ steps:
23+ - uses: actions/checkout@v4
24+ - uses: actions/setup-node@v4
25+ with:
26+ node-version: 24
27+ cache: npm
28+ - run: npm ci
29+ # `build` runs tsc --noEmit first, so this is the type check too. The
30+ # smoke test is not run here: it needs a WebGPU-capable headless Chrome,
31+ # which this runner cannot reliably provide. Run `npm run smoke` locally
32+ # before pushing.
33+ - run: npm run build
34+ # Pages must already be enabled with "GitHub Actions" as the source; the
35+ # workflow token cannot create the site itself.
36+ - uses: actions/configure-pages@v5
37+ - uses: actions/upload-pages-artifact@v3
38+ with:
39+ path: dist
40+ - id: deployment
41+ uses: actions/deploy-pages@v4
.gitignoreadded+3−0View file
@@ -0,0 +1,3 @@
1+node_modules/
2+dist/
3+tmp/
README.mdadded+161−0View file
@@ -0,0 +1,161 @@
1+# acoustic-scattering-3d
2+
3+Sound scattering off obstacles in three dimensions, solved live in the
4+browser. A second-order leapfrog for the acoustic wave equation runs on a
5+cubic grid as WebGPU compute shaders, and the picture is a volume ray march
6+that reads the pressure buffer the solver just wrote, with no readback in the
7+display path.
8+
9+This is the volumetric sibling of
10+[acoustic-scattering-2d](https://github.com/concept-collection/acoustic-scattering-2d).
11+That project's point is its compiler: the solver and the medium are MATLAB,
12+compiled to compute kernels by [numbl](https://numbl.org) live in the page.
13+This one starts simpler, with the solver and the scenes written directly in
14+WGSL and TypeScript, because the third dimension brings enough new problems of
15+its own: memory (a field is n³ floats, so four fields at 192³ are 113 MB),
16+a stricter stability limit, and the question of how to look at a wavefield you
17+are standing outside of. Bringing the editable-MATLAB arrangement over is the
18+obvious next step.
19+
20+## The equation
21+
22+Pressure only, at constant density:
23+
24+```
25+p_tt + 2*sig*p_t = c(x)^2 * lap(p) + s(x, t)
26+```
27+
28+`c` is the sound speed and `sig` the absorption rate, both fields of position
29+that the scene defines. Since the density is constant, the impedance ratio
30+across an interface is just the speed ratio, so a scatterer much faster than
31+its background behaves nearly rigid (sound-hard) and one much slower nearly
32+pressure-release (sound-soft). What this formulation cannot do is set
33+impedance and speed independently, which needs a variable-density divergence
34+form and a second field.
35+
36+Centring both the second time derivative and the damping on step *n* gives an
37+explicit update. The 7-point Laplacian makes the stability condition
38+`c*dt/h <= 1/sqrt(3)`, stricter than the 2D `1/sqrt(2)`; the app sets dt from
39+the fastest speed anywhere in the medium, so a fast scatterer slows the whole
40+run down.
41+
42+## How a step runs
43+
44+One timestep is one compute dispatch, and the update is in place: a step reads
45+the current field `p` at its six neighbours but the previous field `pm` only
46+at its own index, so each thread may overwrite `pm[i]` with the new value.
47+That leaves two pressure buffers instead of three and no copies per step,
48+which matters when a buffer is 28 MB. The two buffers swap roles every step,
49+and the renderer is told which one is current.
50+
51+A frame's worth of steps (up to 64) is recorded into one command encoder, so
52+the source term cannot read anything uploaded between frames. Each step
53+instead reads its own slice of a parameter buffer through a dynamic uniform
54+offset; the whole batch's parameters, including each step's time, are written
55+in one call before the pass. This is the same problem the 2D app solves by
56+carrying time as a grid field. The dynamic-offset answer is cheaper (no field,
57+no kernel), and is available here because the solver is hand-written rather
58+than compiled from MATLAB that only knows about fields.
59+
60+## Looking at it
61+
62+A two-dimensional field is its own picture; a three-dimensional one is not,
63+and every way of drawing it hides something. The app ray-marches the volume:
64+each pixel casts a ray through the cube and accumulates colour front to back,
65+with the pressure through a diverging colormap about zero and opacity rising
66+as a power of |p|, so quiet regions are transparent and wavefronts are what
67+you see. The medium is blended in as a grey cloud so the scatterer is visible
68+inside the field. Drag to orbit, scroll to zoom.
69+
70+Since a volume render of a wavefield is mostly the outside of the wavefield,
71+there is a clip plane on x: pull it in and the interior is exposed, which is
72+the closest thing here to the 2D picture.
73+
74+Two rendering shortcuts are worth knowing. Sampling along the ray is
75+nearest-neighbour, because the field lives in a storage buffer rather than a
76+filterable 3D texture; with the ray step near the cell size this shows mainly
77+as faint stippling on strong fronts. And the compositing is emission only,
78+with no lighting, so depth ordering dims what is behind a strong feature but
79+nothing casts a shadow.
80+
81+## Listening to it
82+
83+There is a microphone: the pressure at one grid point, sampled every timestep.
84+It is written on the GPU by a one-thread dispatch after each step and appended
85+to a trace buffer, because the obvious implementation, reading the field back
86+and picking out one number, costs a GPU-to-CPU round trip per step. The whole
87+trace comes back once, when there is something to play.
88+
89+Everything is SI, so playback is real time: the trace's native rate is 1/dt,
90+around 190 kHz at the defaults, which is above what Web Audio will accept, so
91+it is resampled to 48 kHz on the way out. The content is band-limited far
92+below either rate. The recording is normalized before playback, which
93+discards absolute amplitude: that is what the colour scale is for. It restarts
94+whenever the run does, and whenever the timestep changes, since a trace is one
95+sample per step and two timesteps would be two sample rates in one buffer.
96+
97+As in 2D, a pulse is short: a few cycles at an audible frequency is a few
98+milliseconds however long the simulation runs. For sustained sound, turn
99+`continuous` up and let the run fill some seconds; the trace holds about five
100+seconds at the default timestep.
101+
102+## What it is honest about
103+
104+- **Resolution is the whole game.** A grid solver resolves a wavelength with
105+ some number of cells, and in 3D cells cost their cube. At 128³ over a 2 m
106+ domain, 1.5 kHz has about 11 cells per wavelength; the stats line turns the
107+ number orange when it drops below 8, at which point what is on screen is as
108+ much grid dispersion as sound. There is no `lap4` here yet; the 2D project
109+ shows what a fourth-order stencil buys.
110+- **The absorbing layer is a sponge, not a PML.** Absorption ramps up
111+ quadratically over the outer 15% of each face. An absorbing layer is itself
112+ an impedance mismatch, so it reflects a little; a perfectly matched layer
113+ would do better at the cost of extra fields.
114+- **Single precision.** WebGPU has no f64. The stencil differences lose a few
115+ digits to cancellation, well below the discretization error at these
116+ resolutions.
117+- **No exact-solution comparison yet.** Scattering by a sphere has a classical
118+ series solution (this is the 3D analogue of the cylinder's Bessel-Hankel
119+ series), and comparing against it would put a number on the total error.
120+
121+## Scenes
122+
123+| scene | what it shows |
124+| --- | --- |
125+| Sphere | One spherical scatterer, the reference case. Slow, fast, or absorbing. |
126+| Two spheres | Multiple scattering between a pair: the pattern is not the sum of two singles. |
127+| Aperture | A screen with a circular hole: 3D diffraction, which no 2D slit can show. |
128+| Random medium | Weak random structure everywhere: multiple scattering, and a coda. |
129+
130+The aperture screen is *slower* than the background rather than faster.
131+Reflection at an interface goes as `|c2 - c1|/(c2 + c1)` at constant density,
132+so a screen at c = 0.2 reflects about as much as one at c = 5 would, but the
133+timestep is set by the fastest speed anywhere on the grid: a slow screen is
134+free while a fast one taxes every step of the whole run. What the slow screen
135+costs instead is resolution inside itself, which its own absorption swallows.
136+
137+## Tests
138+
139+```
140+npm run smoke # the page itself, in headless Chrome
141+```
142+
143+The browser check drives the real page: it loads it, waits for the solver to
144+take steps and the frame loop to turn, halves the timestep through the slider
145+and requires dt to follow, checks the microphone is recording, and swaps the
146+scene. It says nothing about whether the picture is right, which is what eyes
147+are for.
148+
149+## Development
150+
151+```
152+npm install
153+npm run dev
154+```
155+
156+No dependencies beyond the build tooling: the solver and renderer are plain
157+TypeScript and WGSL.
158+
159+## License
160+
161+Apache-2.0
index.htmladded+180−0View file
@@ -0,0 +1,180 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="utf-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1" />
6+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><rect width=%22100%22 height=%22100%22 fill=%22%230b0e12%22/><path d=%22M18 26h44l20 14v34H38L18 60z%22 fill=%22none%22 stroke=%22%236b7684%22 stroke-width=%224%22/><circle cx=%2248%22 cy=%2252%22 r=%2213%22 fill=%22%23b40426%22/><circle cx=%2248%22 cy=%2252%22 r=%2222%22 fill=%22none%22 stroke=%22%233b4cc0%22 stroke-width=%225%22/></svg>" />
7+ <title>acoustic-scattering-3d — 3D scattering, live in the browser</title>
8+ <style>
9+ :root {
10+ --bg: #ffffff;
11+ --ink: #1f2328;
12+ --ink-2: #57606a;
13+ --line: #d0d7de;
14+ --accent: #0969da;
15+ --panel-bg: #f4f6f8;
16+ color-scheme: light dark;
17+ }
18+ @media (prefers-color-scheme: dark) {
19+ :root {
20+ --bg: #14171a;
21+ --ink: #e6e9ec;
22+ --ink-2: #9aa4af;
23+ --line: #333b44;
24+ --accent: #58a6ff;
25+ --panel-bg: #14161c;
26+ }
27+ }
28+ body {
29+ margin: 0;
30+ background: var(--bg);
31+ color: var(--ink);
32+ font: 15px/1.5 system-ui, -apple-system, sans-serif;
33+ }
34+ main { max-width: 1100px; margin: 0 auto; padding: 20px 16px 48px; }
35+ h1 { font-size: 20px; margin: 0 0 2px; }
36+ .sub { color: var(--ink-2); margin: 0 0 12px; font-size: 13px; }
37+ .sub a { color: var(--accent); }
38+ .controls {
39+ display: flex; flex-wrap: wrap; gap: 8px 14px; align-items: center;
40+ padding: 5px 0;
41+ }
42+ .controls label { color: var(--ink-2); font-size: 13px; white-space: nowrap; }
43+ select, input[type="number"], button {
44+ font: inherit; font-size: 13px;
45+ color: var(--ink); background: var(--bg);
46+ border: 1px solid var(--line); border-radius: 6px;
47+ padding: 4px 8px;
48+ }
49+ input[type="range"] { width: 8em; vertical-align: middle; accent-color: var(--accent); }
50+ button { cursor: pointer; }
51+ button:hover { border-color: var(--accent); }
52+ button.primary { border-color: var(--accent); color: var(--accent); font-weight: 600; min-width: 5.5em; }
53+ .sliders {
54+ display: grid; grid-template-columns: repeat(auto-fill, minmax(240px, 1fr));
55+ gap: 2px 16px; padding: 4px 0;
56+ }
57+ .slider { display: flex; align-items: center; gap: 8px; font-size: 13px; color: var(--ink-2); }
58+ .slider > span:first-child { flex: 0 0 8.5em; text-align: right; }
59+ .slider > output { flex: 0 0 4.5em; font-variant-numeric: tabular-nums; color: var(--ink); }
60+ .group-title {
61+ font-size: 12px; text-transform: uppercase; letter-spacing: 0.05em;
62+ color: var(--ink-2); margin: 10px 0 0;
63+ }
64+ #micbar { gap: 8px 16px; }
65+ #micparams { display: flex; flex-wrap: wrap; gap: 2px 16px; padding: 0; }
66+ #recinfo { margin-top: 0; }
67+ #stage {
68+ display: flex; gap: 0; margin-top: 12px; align-items: stretch;
69+ border: 1px solid var(--line); border-radius: 8px; overflow: hidden;
70+ background: #0b0e12;
71+ }
72+ .canvas-box { flex: 1; aspect-ratio: 1 / 1; max-height: 70vh; position: relative; }
73+ #view { width: 100%; height: 100%; display: block; cursor: grab; touch-action: none; }
74+ #view:active { cursor: grabbing; }
75+ .colorbar {
76+ display: flex; flex-direction: column; align-items: center; justify-content: center;
77+ gap: 4px; padding: 8px 6px; width: 64px; flex: none; box-sizing: border-box;
78+ color: #9aa4af;
79+ }
80+ .colorbar canvas { border: 1px solid #333b44; border-radius: 2px; }
81+ .colorbar-label { font-size: 11px; font-variant-numeric: tabular-nums; }
82+ .readout { font-variant-numeric: tabular-nums; color: var(--ink); }
83+ .readout.unstable, .warn { color: #b35900; font-weight: 600; }
84+ .stats { margin-top: 8px; font-size: 13px; color: var(--ink-2); font-variant-numeric: tabular-nums; }
85+ .stats b { color: var(--ink); font-weight: 600; }
86+ #blurb { margin-top: 6px; font-size: 13px; color: var(--ink-2); }
87+ #err {
88+ color: #b35900; white-space: pre-wrap; font-size: 13px;
89+ font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
90+ }
91+ </style>
92+ </head>
93+ <body>
94+ <main>
95+ <h1>acoustic-scattering-3d</h1>
96+ <p class="sub">
97+ Sound scattering off obstacles in three dimensions, solved live on your
98+ GPU: a leapfrog on a cubic grid, drawn by marching a ray through the
99+ pressure buffer the solver just wrote. Drag to turn the cube, scroll to
100+ zoom, and pull the clip plane in to see inside it. The flat sibling is
101+ <a href="https://github.com/concept-collection/acoustic-scattering-2d">acoustic-scattering-2d</a>.
102+ </p>
103+ <p class="sub" id="domaininfo"></p>
104+
105+ <div class="controls">
106+ <label>scene
107+ <select id="scene"></select>
108+ </label>
109+ <label title="Grid points per side. Changing it restarts the run.">grid
110+ <select id="gridsize"></select>
111+ </label>
112+ <button id="runpause" class="primary">Run</button>
113+ <button id="restart" title="Back to a silent grid at t = 0">Restart</button>
114+ </div>
115+
116+ <p class="group-title">source</p>
117+ <div class="sliders" id="params"></div>
118+ <p class="group-title" id="scene-title">scene</p>
119+ <div class="sliders" id="sceneparams"></div>
120+
121+ <div class="controls">
122+ <label>colormap
123+ <select id="colormap"></select>
124+ </label>
125+ <label title="Pressure the colour scale saturates at. Auto follows the field; hold keeps the value it had when you switched.">scale
126+ <select id="scalemode">
127+ <option value="auto">auto</option>
128+ <option value="fixed" selected>hold</option>
129+ </select>
130+ </label>
131+ <label title="Shade the sound-speed field alongside the wave">
132+ <input type="checkbox" id="showmedium" checked /> show medium
133+ </label>
134+ <label title="Samples along each ray. Fewer is faster and grainier.">quality
135+ <select id="quality">
136+ <option value="96">low</option>
137+ <option value="192" selected>medium</option>
138+ <option value="320">high</option>
139+ </select>
140+ </label>
141+ <label id="cfl-label" title="The timestep, as a fraction of the largest one this scheme is stable at on this grid. At 1 and above it diverges, which is worth seeing once.">timestep
142+ <input type="range" id="cfl" min="0.05" max="1.2" step="0.05" value="0.5" />
143+ <output id="dtout" class="readout"></output>
144+ </label>
145+ <label title="Timesteps taken between frames">speed
146+ <select id="spf">
147+ <option value="1">1×</option>
148+ <option value="2">2×</option>
149+ <option value="4" selected>4×</option>
150+ <option value="8">8×</option>
151+ <option value="16">16×</option>
152+ <option value="32">32×</option>
153+ <option value="64">64×</option>
154+ </select>
155+ </label>
156+ </div>
157+ <div class="sliders" id="viewparams"></div>
158+
159+ <p class="group-title">microphone</p>
160+ <div class="controls" id="micbar">
161+ <div class="sliders" id="micparams"></div>
162+ <button id="listen" title="Play what the microphone has recorded since the run started, in real time">Listen</button>
163+ <span class="stats" id="recinfo"></span>
164+ </div>
165+
166+ <div id="stage">
167+ <div class="canvas-box"><canvas id="view"></canvas></div>
168+ <div class="colorbar" id="colorbar">
169+ <span class="colorbar-label" id="cbhi"></span>
170+ <canvas id="cbar" width="14" height="220"></canvas>
171+ <span class="colorbar-label" id="cblo"></span>
172+ </div>
173+ </div>
174+ <p class="stats" id="stats"></p>
175+ <p id="blurb"></p>
176+ <p id="err"></p>
177+ </main>
178+ <script type="module" src="/src/main.ts"></script>
179+ </body>
180+</html>
package-lock.jsonadded+2142−0View file
This diff is 2,147 lines long and is not shown.
package.jsonadded+23−0View file
@@ -0,0 +1,23 @@
1+{
2+ "name": "acoustic-scattering-3d",
3+ "version": "0.1.0",
4+ "description": "3D acoustic scattering solved live in the browser on WebGPU",
5+ "type": "module",
6+ "engines": {
7+ "node": ">=22.6"
8+ },
9+ "license": "Apache-2.0",
10+ "scripts": {
11+ "dev": "vite",
12+ "build": "tsc --noEmit && vite build",
13+ "smoke": "vite build && node scripts/smoke.mjs",
14+ "test": "npm run smoke"
15+ },
16+ "devDependencies": {
17+ "@types/node": "^26.1.1",
18+ "@webgpu/types": "^0.1.44",
19+ "puppeteer-core": "^23.11.1",
20+ "typescript": "^5.5.0",
21+ "vite": "^5.4.0"
22+ }
23+}
scripts/smoke.mjsadded+176−0View file
@@ -0,0 +1,176 @@
1+/**
2+ * Does the page actually run in a browser?
3+ *
4+ * Serves dist/ and opens it in headless Chrome (hardware WebGPU if there is
5+ * any, SwiftShader otherwise), then waits for the simulation to report that it
6+ * has taken steps and the frame loop to turn. Any console error, page error or
7+ * failed request fails the run.
8+ *
9+ * This checks that the app *works*, not that it looks right — a headless
10+ * browser has no opinion about whether the wavefronts are in the right place.
11+ * Run it after `vite build`: node scripts/smoke.mjs
12+ */
13+import { createServer } from 'node:http';
14+import { readFile } from 'node:fs/promises';
15+import { extname, join } from 'node:path';
16+import puppeteer from 'puppeteer-core';
17+
18+const DIST = new URL('../dist/', import.meta.url).pathname;
19+const CHROME = process.env.CHROME_PATH ?? '/usr/bin/google-chrome';
20+const MIME = {
21+ '.html': 'text/html',
22+ '.js': 'text/javascript',
23+ '.css': 'text/css',
24+ '.json': 'application/json',
25+};
26+
27+const server = createServer(async (req, res) => {
28+ try {
29+ const path = req.url === '/' ? '/index.html' : req.url.split('?')[0];
30+ const data = await readFile(join(DIST, path));
31+ res.writeHead(200, { 'content-type': MIME[extname(path)] ?? 'application/octet-stream' });
32+ res.end(data);
33+ } catch {
34+ res.writeHead(404);
35+ res.end('not found');
36+ }
37+});
38+await new Promise((r) => server.listen(0, '127.0.0.1', r));
39+const port = server.address().port;
40+
41+const flagSets = [
42+ ['--headless=new', '--no-sandbox', '--enable-unsafe-webgpu', '--enable-features=Vulkan'],
43+ [
44+ '--headless=new',
45+ '--no-sandbox',
46+ '--enable-unsafe-webgpu',
47+ '--use-webgpu-adapter=swiftshader',
48+ '--enable-unsafe-swiftshader',
49+ ],
50+];
51+
52+let ok = false;
53+let lastFailure = 'never ran';
54+for (const flags of flagSets) {
55+ const browser = await puppeteer.launch({
56+ executablePath: CHROME,
57+ args: [...flags],
58+ protocolTimeout: 600_000,
59+ });
60+ const problems = [];
61+ try {
62+ const page = await browser.newPage();
63+ page.on('console', (m) => {
64+ if (m.type() === 'error') problems.push(`console: ${m.text()}`);
65+ });
66+ page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`));
67+ page.on('requestfailed', (r) => problems.push(`request failed: ${r.url()}`));
68+
69+ await page.goto(`http://127.0.0.1:${port}/index.html`, { waitUntil: 'load' });
70+
71+ // The page loads paused; the stats line says so once the loop is turning.
72+ await page.waitForFunction(
73+ () => /paused/.test(document.getElementById('stats')?.textContent ?? ''),
74+ { timeout: 120_000, polling: 250 },
75+ );
76+ await page.click('#runpause');
77+
78+ /** The stats line reports the step count, so waiting on it says both that
79+ * the solver ran and that the frame loop is turning. */
80+ const running = async (label) => {
81+ await page.waitForFunction(
82+ () => {
83+ const m = /step (\d+)/.exec(document.getElementById('stats')?.textContent ?? '');
84+ return m ? Number(m[1]) > 20 : false;
85+ },
86+ { timeout: 180_000, polling: 500 },
87+ );
88+ console.log(` ${label}: ${(await page.$eval('#stats', (n) => n.textContent)).trim()}`);
89+ };
90+ await running('start');
91+
92+ // The timestep slider must reach the solver: halving the CFL fraction
93+ // halves dt, which the stats line reports.
94+ const dtNow = async () =>
95+ Number(/dt = ([0-9.e+-]+)/.exec(await page.$eval('#stats', (n) => n.textContent ?? ''))?.[1]);
96+ const dtBefore = await dtNow();
97+ await page.evaluate(() => {
98+ const slider = document.getElementById('cfl');
99+ slider.value = '0.25';
100+ slider.dispatchEvent(new Event('input'));
101+ });
102+ await page.waitForFunction(
103+ (before) => {
104+ const m = /dt = ([0-9.e+-]+)/.exec(document.getElementById('stats')?.textContent ?? '');
105+ return m ? Math.abs(Number(m[1]) - before / 2) < before / 20 : false;
106+ },
107+ { timeout: 30_000, polling: 250 },
108+ dtBefore,
109+ );
110+ console.log(
111+ ` timestep: dt ${dtBefore.toExponential(2)} -> ${(await dtNow()).toExponential(2)} at half the CFL`,
112+ );
113+
114+ // The microphone must be recording, one sample per timestep. (Both
115+ // numbers are host-side counters, so this checks the wiring, not the GPU
116+ // trace — that would need a readback the page only does on Listen.)
117+ await page.waitForFunction(
118+ () => /[\d,]+ samples/.test(document.getElementById('recinfo')?.textContent ?? ''),
119+ { timeout: 60_000, polling: 250 },
120+ );
121+ const recorded = Number(
122+ /([\d,]+) samples/
123+ .exec(await page.$eval('#recinfo', (n) => n.textContent ?? ''))?.[1]
124+ .replace(/,/g, ''),
125+ );
126+ const stepsSoFar = Number(
127+ /step (\d+)/.exec(await page.$eval('#stats', (n) => n.textContent ?? ''))?.[1],
128+ );
129+ console.log(` microphone: ${(await page.$eval('#recinfo', (n) => n.textContent ?? '')).trim()}`);
130+ if (!(recorded > 0) || recorded > stepsSoFar) {
131+ problems.push(`microphone recorded ${recorded} samples in ${stepsSoFar} steps`);
132+ }
133+
134+ // Swapping the scene rebuilds the medium. The aperture screen is slower
135+ // than the background (the smoothing keeps the grid minimum above the
136+ // nominal 0.2), so a reported cmin below 1 is the evidence.
137+ await page.select('#scene', 'aperture');
138+ await page.waitForFunction(
139+ () => /c ∈ \[0\./.test(document.getElementById('stats')?.textContent ?? ''),
140+ { timeout: 60_000, polling: 250 },
141+ );
142+ await running('after scene swap');
143+
144+ const report = await page.evaluate(() => ({
145+ stats: document.getElementById('stats')?.textContent ?? '',
146+ err: document.getElementById('err')?.textContent ?? '',
147+ painted: (() => {
148+ const canvas = document.getElementById('view');
149+ return canvas instanceof HTMLCanvasElement && canvas.width > 0;
150+ })(),
151+ }));
152+
153+ console.log(`flags: ${flags.join(' ')}`);
154+ console.log(report.stats.trim());
155+ if (report.err) problems.push(`page error box: ${report.err}`);
156+ if (!report.painted) problems.push('canvas was never sized');
157+ if (problems.length === 0) {
158+ ok = true;
159+ } else {
160+ lastFailure = problems.join('\n');
161+ }
162+ } catch (e) {
163+ lastFailure = [`${e}`, ...problems].join('\n');
164+ } finally {
165+ await browser.close();
166+ }
167+ if (ok) break;
168+}
169+
170+server.close();
171+if (!ok) {
172+ console.error(`smoke: FAILED\n${lastFailure}`);
173+ process.exit(1);
174+}
175+console.log('smoke: the page runs');
176+process.exit(0);
src/audio.tsadded+69−0View file
@@ -0,0 +1,69 @@
1+/**
2+ * Playing the microphone trace.
3+ *
4+ * The trace holds one sample per timestep, dt seconds apart, so its native
5+ * sample rate is 1/dt — around 190 kHz at the default settings, which is
6+ * beyond what an AudioBuffer must accept. The trace is resampled to 48 kHz by
7+ * linear interpolation, which is harmless here: the content is band-limited
8+ * far below either rate (the source tops out at 4 kHz).
9+ *
10+ * The AudioContext is created once and resumed on every play. Browsers start
11+ * a context suspended unless it is created directly inside a user gesture,
12+ * and the Listen click's gesture is over by the time the trace has come back
13+ * from the GPU — resuming is what makes playback actually sound. This is the
14+ * bug the first version of this file had: it played, silently, into a
15+ * suspended context.
16+ *
17+ * The recording is normalized before playback, which discards absolute
18+ * amplitude: that is what the colour scale is for. A few milliseconds of fade
19+ * are applied at each end — a trace that starts or ends away from zero is a
20+ * step, and a step is a click that has nothing to do with the simulation.
21+ */
22+
23+let context: AudioContext | null = null;
24+let playing: AudioBufferSourceNode | null = null;
25+
26+export interface Played {
27+ /** Seconds of audio actually played. */
28+ duration: number;
29+ /** Peak |p| of the trace before normalization. Zero means silence went by. */
30+ peak: number;
31+}
32+
33+export async function playTrace(samples: Float32Array, dt: number): Promise<Played> {
34+ const rate = 48000;
35+ const n = Math.max(1, Math.round(samples.length * dt * rate));
36+
37+ let peak = 0;
38+ for (let i = 0; i < samples.length; i++) peak = Math.max(peak, Math.abs(samples[i]));
39+ const g = peak > 0 ? 0.9 / peak : 0;
40+
41+ context ??= new AudioContext();
42+ if (context.state === 'suspended') await context.resume();
43+
44+ const buf = context.createBuffer(1, n, rate);
45+ const ch = buf.getChannelData(0);
46+ for (let i = 0; i < n; i++) {
47+ const s = i / (rate * dt);
48+ const k = Math.min(samples.length - 2, Math.floor(s));
49+ const f = Math.min(1, s - k);
50+ ch[i] = g * ((1 - f) * samples[k] + f * samples[k + 1]);
51+ }
52+ const fade = Math.min(Math.round(0.005 * rate), Math.floor(n / 2));
53+ for (let i = 0; i < fade; i++) {
54+ const w = i / fade;
55+ ch[i] *= w;
56+ ch[n - 1 - i] *= w;
57+ }
58+
59+ if (playing) playing.stop();
60+ const src = context.createBufferSource();
61+ src.buffer = buf;
62+ src.connect(context.destination);
63+ src.onended = () => {
64+ if (playing === src) playing = null;
65+ };
66+ src.start();
67+ playing = src;
68+ return { duration: n / rate, peak };
69+}
src/device.tsadded+32−0View file
@@ -0,0 +1,32 @@
1+/** The GPU device, requested the same way everywhere (app and scripts). */
2+export async function requestAcousticDevice(): Promise<GPUDevice> {
3+ if (!navigator.gpu) {
4+ throw new Error(
5+ 'this browser has no WebGPU. Chrome and Edge 113+, Safari 26+, and ' +
6+ 'Firefox 141+ on Windows have it; on Linux Firefox and Chrome may need ' +
7+ 'it enabled explicitly.',
8+ );
9+ }
10+ const adapter = await navigator.gpu.requestAdapter({ powerPreference: 'high-performance' });
11+ if (!adapter) throw new Error('WebGPU found no adapter on this machine.');
12+
13+ // A field of 192^3 f32 is 28 MB and there are four of them, so the default
14+ // 128 MB buffer limit and 256 MB total are both too small at the top grid
15+ // size. Ask for what the adapter will give.
16+ const lim = adapter.limits;
17+ return adapter.requestDevice({
18+ requiredLimits: {
19+ maxStorageBufferBindingSize: lim.maxStorageBufferBindingSize,
20+ maxBufferSize: lim.maxBufferSize,
21+ },
22+ });
23+}
24+
25+/** Largest cube side this device can hold four fields of. */
26+export function maxGridSide(device: GPUDevice): number {
27+ const perField = device.limits.maxStorageBufferBindingSize;
28+ for (const n of [192, 160, 128, 96, 64]) {
29+ if (4 * n * n * n <= perField) return n;
30+ }
31+ return 64;
32+}
src/grid.tsadded+48−0View file
@@ -0,0 +1,48 @@
1+/**
2+ * The computational grid: a cube, cell-centred, uniform in all three
3+ * directions.
4+ *
5+ * Fields are flattened x-fastest then y — the point (ix, iy, iz) is element
6+ * `ix + n*iy + n*n*iz` — which is the order the stencil shader indexes in and
7+ * the order the renderer reads a voxel in. Everything else treats a field as
8+ * an opaque npts-long array of f32.
9+ *
10+ * Cell-centred rather than node-centred so that no grid point sits exactly on
11+ * the outer boundary: the stencil takes the field outside the domain to be
12+ * zero, and the absorbing layer is meant to have swallowed the wave before it
13+ * gets there.
14+ *
15+ * `n` is required to be a multiple of 16 because the update shader dispatches
16+ * 8x8x4 workgroups and does not want a partial one at the far corner. That is
17+ * a convenience, not a physical constraint.
18+ */
19+export interface Grid {
20+ n: number;
21+ npts: number;
22+ /** Side length of the cubic domain, in metres, centred on the origin. */
23+ L: number;
24+ /** Grid spacing, L/n, in metres. */
25+ h: number;
26+}
27+
28+export function makeGrid(n: number, L: number): Grid {
29+ return { n, npts: n * n * n, L, h: L / n };
30+}
31+
32+/** The coordinate of cell index `i` along one axis, in metres. */
33+export const coord = (i: number, g: Grid): number => -g.L / 2 + (i + 0.5) * g.h;
34+
35+/**
36+ * The timestep the explicit leapfrog is stable at, in seconds.
37+ *
38+ * Leapfrog on p_tt = c^2 lap(p) is stable while dt^2 c^2 |lap|max <= 4, and
39+ * the 7-point Laplacian's extreme eigenvalue is 12/h^2, so the condition is
40+ * c*dt/h <= 1/sqrt(3) = 0.577. That is stricter than the 0.707 of the flat
41+ * case, which is the usual price of a dimension. `cfl` is the fraction of the
42+ * limit to run at, and `cmax` is the fastest sound speed anywhere in the
43+ * medium: a scatterer faster than the background sets the timestep for the
44+ * whole grid.
45+ */
46+export function stableDt(h: number, cmax: number, cfl = 0.5): number {
47+ return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
48+}
src/main.tsadded+403−0View file
@@ -0,0 +1,403 @@
1+/**
2+ * The app: wiring between the solver (sim.ts), the scenes (scenes.ts), the
3+ * volume renderer (render/volume.ts) and the page.
4+ */
5+import { requestAcousticDevice, maxGridSide } from './device.ts';
6+import { makeGrid, stableDt, type Grid } from './grid.ts';
7+import { Sim, MAX_STEPS_PER_FRAME, MIC_CAPACITY, type SourceParams } from './sim.ts';
8+import { playTrace } from './audio.ts';
9+import { scenes, sceneKeys, buildMedium, type Medium, type SceneParam } from './scenes.ts';
10+import { VolumeView, type Camera } from './render/volume.ts';
11+import { colormaps, colormapNames } from './render/colormaps.ts';
12+import { C_AIR, DOMAIN, POOR_RESOLUTION, fmtLength, fmtTime } from './units.ts';
13+
14+const $ = <T extends HTMLElement>(id: string): T => {
15+ const el = document.getElementById(id);
16+ if (!el) throw new Error(`no #${id}`);
17+ return el as T;
18+};
19+
20+const err = $('err');
21+const report = (e: unknown) => {
22+ err.textContent = e instanceof Error ? e.message : String(e);
23+};
24+
25+/** A labelled slider with a live readout, in one of the .sliders grids. */
26+function slider(
27+ parent: HTMLElement,
28+ p: SceneParam,
29+ onInput: (v: number) => void,
30+): { set: (v: number) => void } {
31+ const row = document.createElement('label');
32+ row.className = 'slider';
33+ const name = document.createElement('span');
34+ name.textContent = p.label;
35+ const input = document.createElement('input');
36+ input.type = 'range';
37+ input.min = String(p.min);
38+ input.max = String(p.max);
39+ input.step = String(p.step);
40+ input.value = String(p.value);
41+ const out = document.createElement('output');
42+ const fmt = (v: number) => `${+v.toPrecision(3)}${p.unit ? ' ' + p.unit : ''}`;
43+ out.textContent = fmt(p.value);
44+ input.addEventListener('input', () => {
45+ const v = Number(input.value);
46+ out.textContent = fmt(v);
47+ onInput(v);
48+ });
49+ row.append(name, input, out);
50+ parent.append(row);
51+ return {
52+ set: (v) => {
53+ input.value = String(v);
54+ out.textContent = fmt(v);
55+ },
56+ };
57+}
58+
59+async function start() {
60+ const device = await requestAcousticDevice();
61+ device.addEventListener('uncapturederror', (e) => report((e as GPUUncapturedErrorEvent).error.message));
62+
63+ // --- state ---------------------------------------------------------------
64+ const gridSel = $<HTMLSelectElement>('gridsize');
65+ for (const n of [64, 96, 128, 160, 192]) {
66+ if (n > maxGridSide(device)) break;
67+ const opt = document.createElement('option');
68+ opt.value = String(n);
69+ opt.textContent = `${n}³`;
70+ if (n === 128) opt.selected = true;
71+ gridSel.append(opt);
72+ }
73+
74+ let sceneKey = 'sphere';
75+ let grid: Grid = makeGrid(Number(gridSel.value), DOMAIN);
76+ let cfl = 0.5;
77+ let running = false;
78+ let stepsPerFrame = 4;
79+
80+ const source: SourceParams = {
81+ f: 1200,
82+ cycles: 2,
83+ cw: 0,
84+ point: 0,
85+ x0: -0.32 * DOMAIN,
86+ y0: 0,
87+ z0: 0,
88+ w: 0.02 * DOMAIN,
89+ };
90+ const sceneValues: Record<string, Record<string, number>> = {};
91+ for (const k of sceneKeys) {
92+ sceneValues[k] = Object.fromEntries(scenes[k].params.map((p) => [p.key, p.value]));
93+ }
94+
95+ const camera: Camera = { az: -2.35, el: 0.5, dist: 2.6 * DOMAIN };
96+ const viewState = { opacity: 1.6, contrast: 1.6, clip: 0.5 };
97+ let quality = 192;
98+ let colormap = 'coolwarm';
99+ let scaleMode: 'auto' | 'fixed' = 'fixed';
100+ let scale = 1;
101+ let showMedium = true;
102+ /** Microphone position, metres. In the forward-scattering shadow of the
103+ * default sphere, which is where there is something to hear. */
104+ const mic = { x: 0.25 * DOMAIN, y: 0, z: 0 };
105+
106+ let medium: Medium = buildMedium(scenes[sceneKey], grid, sceneValues[sceneKey]);
107+ let sim = new Sim(device, grid, medium, source, stableDt(grid.h, medium.cmax, cfl));
108+ const view = new VolumeView(device, $<HTMLCanvasElement>('view'));
109+ view.setSource(sim.pressures, sim.speed, grid.n, grid.L);
110+ view.setColormap(colormaps[colormap]);
111+
112+ const applyProbe = () => {
113+ const toIdx = (v: number) => Math.floor((v + grid.L / 2) / grid.h);
114+ sim.setProbe(toIdx(mic.x), toIdx(mic.y), toIdx(mic.z));
115+ };
116+ applyProbe();
117+
118+ // --- rebuilds ------------------------------------------------------------
119+ const applyDt = () => {
120+ const dt = stableDt(grid.h, medium.cmax, cfl);
121+ if (dt !== sim.dt) {
122+ sim.dt = dt;
123+ // One sample per step: two timesteps would be two sample rates in one
124+ // trace.
125+ sim.resetTrace();
126+ }
127+ };
128+
129+ /** New medium on the existing grid (a scene parameter moved): rebuild the
130+ * sim but keep the run going conceptually — the field restarts, which is
131+ * honest, since the old field belongs to the old medium. */
132+ const rebuildSim = () => {
133+ sim.destroy();
134+ medium = buildMedium(scenes[sceneKey], grid, sceneValues[sceneKey]);
135+ sim = new Sim(device, grid, medium, source, stableDt(grid.h, medium.cmax, cfl));
136+ view.setSource(sim.pressures, sim.speed, grid.n, grid.L);
137+ applyProbe();
138+ updateStatics();
139+ };
140+
141+ // --- controls ------------------------------------------------------------
142+ const sceneSel = $<HTMLSelectElement>('scene');
143+ for (const k of sceneKeys) {
144+ const opt = document.createElement('option');
145+ opt.value = k;
146+ opt.textContent = scenes[k].label;
147+ sceneSel.append(opt);
148+ }
149+
150+ const sceneParamsBox = $('sceneparams');
151+ const buildSceneSliders = () => {
152+ sceneParamsBox.replaceChildren();
153+ $('scene-title').textContent = `scene — ${scenes[sceneKey].label.toLowerCase()}`;
154+ $('blurb').textContent = scenes[sceneKey].blurb;
155+ for (const p of scenes[sceneKey].params) {
156+ slider(sceneParamsBox, { ...p, value: sceneValues[sceneKey][p.key] }, (v) => {
157+ sceneValues[sceneKey][p.key] = v;
158+ rebuildSim();
159+ });
160+ }
161+ };
162+
163+ sceneSel.addEventListener('change', () => {
164+ sceneKey = sceneSel.value;
165+ const want = scenes[sceneKey].source;
166+ if (want) {
167+ Object.assign(source, want);
168+ const s = source as unknown as Record<string, number>;
169+ for (const [k, set] of Object.entries(sourceSliders)) set.set(s[k]);
170+ }
171+ buildSceneSliders();
172+ rebuildSim();
173+ });
174+
175+ gridSel.addEventListener('change', () => {
176+ grid = makeGrid(Number(gridSel.value), DOMAIN);
177+ rebuildSim();
178+ });
179+
180+ const runBtn = $<HTMLButtonElement>('runpause');
181+ const setRunning = (r: boolean) => {
182+ running = r;
183+ runBtn.textContent = r ? 'Pause' : 'Run';
184+ };
185+ runBtn.addEventListener('click', () => setRunning(!running));
186+ $('restart').addEventListener('click', () => {
187+ sim.restart();
188+ });
189+
190+ // Source sliders. Moving one does not restart: the source term reads these
191+ // every step, so the change simply takes effect.
192+ const sourceSliders: Record<string, { set: (v: number) => void }> = {};
193+ const params = $('params');
194+ const sourceDefs: (SceneParam & { key: keyof SourceParams })[] = [
195+ { key: 'f', label: 'frequency', min: 200, max: 4000, step: 20, value: source.f, unit: 'Hz' },
196+ { key: 'cycles', label: 'pulse length', min: 0.5, max: 12, step: 0.5, value: source.cycles },
197+ { key: 'cw', label: 'continuous', min: 0, max: 1, step: 0.05, value: source.cw },
198+ { key: 'point', label: 'plane ↔ point', min: 0, max: 1, step: 0.05, value: source.point },
199+ { key: 'x0', label: 'x position', min: -0.4 * DOMAIN, max: 0.4 * DOMAIN, step: 0.01, value: source.x0, unit: 'm' },
200+ { key: 'y0', label: 'y position', min: -0.4 * DOMAIN, max: 0.4 * DOMAIN, step: 0.01, value: source.y0, unit: 'm' },
201+ { key: 'z0', label: 'z position', min: -0.4 * DOMAIN, max: 0.4 * DOMAIN, step: 0.01, value: source.z0, unit: 'm' },
202+ ];
203+ for (const p of sourceDefs) {
204+ sourceSliders[p.key] = slider(params, p, (v) => {
205+ source[p.key] = v;
206+ });
207+ }
208+
209+ // The microphone: three position sliders, a Listen button, a readout.
210+ const micparams = $('micparams');
211+ for (const axis of ['x', 'y', 'z'] as const) {
212+ slider(
213+ micparams,
214+ { key: axis, label: `mic ${axis}`, min: -0.45 * DOMAIN, max: 0.45 * DOMAIN, step: 0.01, value: mic[axis], unit: 'm' },
215+ (v) => {
216+ mic[axis] = v;
217+ applyProbe();
218+ },
219+ );
220+ }
221+ let playNote = '';
222+ $('listen').addEventListener('click', () => {
223+ sim
224+ .readTrace()
225+ .then(async (tr) => {
226+ if (!tr || tr.length < 2) return;
227+ const played = await playTrace(tr, sim.dt);
228+ playNote =
229+ played.peak > 0
230+ ? ` — played ${fmtTime(played.duration)}`
231+ : ' — the microphone heard only silence';
232+ })
233+ .catch(report);
234+ });
235+
236+ const viewparams = $('viewparams');
237+ slider(viewparams, { key: 'op', label: 'opacity', min: 0.2, max: 6, step: 0.1, value: viewState.opacity }, (v) => {
238+ viewState.opacity = v;
239+ });
240+ slider(viewparams, { key: 'ct', label: 'contrast', min: 0.5, max: 3, step: 0.1, value: viewState.contrast }, (v) => {
241+ viewState.contrast = v;
242+ });
243+ slider(viewparams, { key: 'clip', label: 'clip x', min: -0.5, max: 0.5, step: 0.01, value: viewState.clip }, (v) => {
244+ viewState.clip = v;
245+ });
246+
247+ const cmapSel = $<HTMLSelectElement>('colormap');
248+ for (const name of colormapNames) {
249+ const opt = document.createElement('option');
250+ opt.value = name;
251+ opt.textContent = name;
252+ cmapSel.append(opt);
253+ }
254+ cmapSel.value = colormap;
255+ cmapSel.addEventListener('change', () => {
256+ colormap = cmapSel.value;
257+ view.setColormap(colormaps[colormap]);
258+ drawColorbar();
259+ });
260+
261+ $<HTMLSelectElement>('scalemode').addEventListener('change', (e) => {
262+ scaleMode = (e.target as HTMLSelectElement).value as typeof scaleMode;
263+ });
264+ $<HTMLInputElement>('showmedium').addEventListener('change', (e) => {
265+ showMedium = (e.target as HTMLInputElement).checked;
266+ });
267+ $<HTMLSelectElement>('quality').addEventListener('change', (e) => {
268+ quality = Number((e.target as HTMLSelectElement).value);
269+ });
270+ $<HTMLSelectElement>('spf').addEventListener('change', (e) => {
271+ stepsPerFrame = Math.min(Number((e.target as HTMLSelectElement).value), MAX_STEPS_PER_FRAME);
272+ });
273+
274+ const cflSlider = $<HTMLInputElement>('cfl');
275+ const dtout = $('dtout');
276+ cflSlider.addEventListener('input', () => {
277+ cfl = Number(cflSlider.value);
278+ applyDt();
279+ updateStatics();
280+ });
281+
282+ // --- orbit ---------------------------------------------------------------
283+ const canvas = $<HTMLCanvasElement>('view');
284+ let dragging = false;
285+ let last = [0, 0];
286+ canvas.addEventListener('pointerdown', (e) => {
287+ dragging = true;
288+ last = [e.clientX, e.clientY];
289+ canvas.setPointerCapture(e.pointerId);
290+ });
291+ canvas.addEventListener('pointermove', (e) => {
292+ if (!dragging) return;
293+ camera.az -= (e.clientX - last[0]) * 0.008;
294+ camera.el = Math.min(1.5, Math.max(-1.5, camera.el + (e.clientY - last[1]) * 0.008));
295+ last = [e.clientX, e.clientY];
296+ });
297+ canvas.addEventListener('pointerup', () => (dragging = false));
298+ canvas.addEventListener('pointercancel', () => (dragging = false));
299+ canvas.addEventListener(
300+ 'wheel',
301+ (e) => {
302+ e.preventDefault();
303+ camera.dist = Math.min(8 * DOMAIN, Math.max(1.1 * DOMAIN, camera.dist * Math.exp(e.deltaY * 0.001)));
304+ },
305+ { passive: false },
306+ );
307+
308+ // --- readouts ------------------------------------------------------------
309+ const stats = $('stats');
310+ const domaininfo = $('domaininfo');
311+
312+ const drawColorbar = () => {
313+ const cb = $<HTMLCanvasElement>('cbar');
314+ const ctx = cb.getContext('2d')!;
315+ const f = colormaps[colormap];
316+ for (let i = 0; i < cb.height; i++) {
317+ const [r, g, b] = f(1 - i / (cb.height - 1));
318+ ctx.fillStyle = `rgb(${r},${g},${b})`;
319+ ctx.fillRect(0, i, cb.width, 1);
320+ }
321+ };
322+ drawColorbar();
323+
324+ const updateStatics = () => {
325+ applyDt();
326+ const lam = C_AIR / source.f;
327+ const cells = lam / grid.h;
328+ domaininfo.textContent =
329+ `A ${fmtLength(grid.L)} cube of air at ${grid.n}³ points (${fmtLength(grid.h)} cells); ` +
330+ `at ${source.f} Hz a wavelength is ${fmtLength(lam)}, ${cells.toFixed(1)} cells.`;
331+ dtout.textContent = fmtTime(sim.dt);
332+ dtout.classList.toggle('unstable', cfl >= 1);
333+ };
334+ updateStatics();
335+
336+ // --- frame loop ----------------------------------------------------------
337+ let frames = 0;
338+ let lastFps = performance.now();
339+ let msPerFrame = 0;
340+
341+ const frame = () => {
342+ view.resize();
343+ if (running) {
344+ sim.run(stepsPerFrame);
345+ if (scaleMode === 'auto') {
346+ // A rebuild can destroy the sim while a readback is in flight; a
347+ // rejection then is about the old sim and not worth reporting.
348+ sim
349+ .peak()
350+ .then((m) => {
351+ if (m !== null && m > 0) scale = Math.max(scale * 0.98, m * 0.85);
352+ })
353+ .catch(() => {});
354+ }
355+ }
356+ view.draw(sim.pressureIndex, {
357+ camera,
358+ scale: Math.max(scale, 1e-20),
359+ opacity: viewState.opacity,
360+ contrast: viewState.contrast,
361+ steps: quality,
362+ clipX: viewState.clip * grid.L,
363+ medium: showMedium ? 0.35 : 0,
364+ cref: C_AIR,
365+ cdev: Math.max(medium.cmax - C_AIR, C_AIR - medium.cmin, 1e-6),
366+ mic,
367+ });
368+
369+ frames++;
370+ const now = performance.now();
371+ if (now - lastFps > 500) {
372+ msPerFrame = (now - lastFps) / frames;
373+ frames = 0;
374+ lastFps = now;
375+
376+ const cellsPerLam = C_AIR / source.f / grid.h;
377+ const resolution =
378+ cellsPerLam < POOR_RESOLUTION
379+ ? ` — <span class="warn">${cellsPerLam.toFixed(1)} cells per wavelength: mostly grid dispersion</span>`
380+ : '';
381+ const rate = running ? ` — ${msPerFrame.toFixed(1)} ms/frame` : ' — paused';
382+ stats.innerHTML =
383+ `step <b>${sim.steps}</b> — t = ${fmtTime(sim.t)} — dt = ${sim.dt.toExponential(3)} s` +
384+ ` — c ∈ [${(medium.cmin / C_AIR).toPrecision(3)}, ${(medium.cmax / C_AIR).toPrecision(3)}]·c₀` +
385+ rate +
386+ resolution;
387+ $('recinfo').textContent =
388+ sim.recorded === 0
389+ ? 'nothing recorded yet'
390+ : `${sim.recorded.toLocaleString('en-US')} samples — ${fmtTime(sim.recorded * sim.dt)}` +
391+ (sim.recorded >= MIC_CAPACITY ? ' (full)' : '') +
392+ playNote;
393+ }
394+ $('cbhi').textContent = `+${scale.toPrecision(2)}`;
395+ $('cblo').textContent = `−${scale.toPrecision(2)}`;
396+ requestAnimationFrame(frame);
397+ };
398+
399+ buildSceneSliders();
400+ requestAnimationFrame(frame);
401+}
402+
403+start().catch(report);
src/render/colormaps.tsadded+113−0View file
@@ -0,0 +1,113 @@
1+/**
2+ * Colormaps: each maps a normalized value in [0, 1] to [r, g, b] in [0, 255].
3+ * Adapted from figpack's SphereEmbedding view (figpack_experimental).
4+ */
5+
6+export type ColormapFunc = (t: number) => [number, number, number];
7+
8+const clamp01 = (t: number) => Math.max(0, Math.min(1, t));
9+
10+// Piecewise-linear interpolation through control points (r, g, b in 0-255)
11+const makeInterpolated = (stops: [number, number, number][]): ColormapFunc => {
12+ const n = stops.length;
13+ return (t: number) => {
14+ t = clamp01(t);
15+ const x = t * (n - 1);
16+ const i = Math.min(n - 2, Math.floor(x));
17+ const f = x - i;
18+ const a = stops[i];
19+ const b = stops[i + 1];
20+ return [
21+ Math.round(a[0] + (b[0] - a[0]) * f),
22+ Math.round(a[1] + (b[1] - a[1]) * f),
23+ Math.round(a[2] + (b[2] - a[2]) * f),
24+ ];
25+ };
26+};
27+
28+// Control points sampled from matplotlib colormaps
29+const viridis = makeInterpolated([
30+ [68, 1, 84],
31+ [72, 40, 120],
32+ [62, 74, 137],
33+ [49, 104, 142],
34+ [38, 130, 142],
35+ [31, 158, 137],
36+ [53, 183, 121],
37+ [109, 205, 89],
38+ [180, 222, 44],
39+ [253, 231, 37],
40+]);
41+
42+const plasma = makeInterpolated([
43+ [13, 8, 135],
44+ [84, 2, 163],
45+ [139, 10, 165],
46+ [185, 50, 137],
47+ [219, 92, 104],
48+ [244, 136, 73],
49+ [254, 188, 43],
50+ [240, 249, 33],
51+]);
52+
53+const inferno = makeInterpolated([
54+ [0, 0, 4],
55+ [40, 11, 84],
56+ [101, 21, 110],
57+ [159, 42, 99],
58+ [212, 72, 66],
59+ [245, 125, 21],
60+ [250, 193, 39],
61+ [252, 255, 164],
62+]);
63+
64+const coolwarm = makeInterpolated([
65+ [59, 76, 192],
66+ [124, 159, 249],
67+ [192, 212, 245],
68+ [242, 242, 242],
69+ [245, 195, 157],
70+ [222, 96, 77],
71+ [180, 4, 38],
72+]);
73+
74+// matplotlib's `seismic`: harder contrast about the middle than coolwarm, and
75+// dark at both ends, which suits a wavefield whose interesting parts are the
76+// extremes.
77+const seismic = makeInterpolated([
78+ [0, 0, 76],
79+ [0, 0, 255],
80+ [255, 255, 255],
81+ [255, 0, 0],
82+ [128, 0, 0],
83+]);
84+
85+const jet = makeInterpolated([
86+ [0, 0, 128],
87+ [0, 0, 255],
88+ [0, 255, 255],
89+ [0, 255, 0],
90+ [255, 255, 0],
91+ [255, 0, 0],
92+ [128, 0, 0],
93+]);
94+
95+const grayscale: ColormapFunc = (t: number) => {
96+ const v = Math.round(clamp01(t) * 255);
97+ return [v, v, v];
98+};
99+
100+/** Diverging maps first: the pressure field is signed and is drawn
101+ * symmetrically about zero, so a map with a distinct middle is what makes
102+ * the wavefronts read. */
103+export const colormaps: Record<string, ColormapFunc> = {
104+ coolwarm,
105+ seismic,
106+ grayscale,
107+ viridis,
108+ plasma,
109+ inferno,
110+ jet,
111+};
112+
113+export const colormapNames = Object.keys(colormaps);
src/render/volume.tsadded+357−0View file
@@ -0,0 +1,357 @@
1+/**
2+ * Drawing the pressure field, straight out of the buffer the solver wrote.
3+ *
4+ * There is no readback in the display path: the fragment shader reads the
5+ * solver's storage buffer directly, so a frame costs one draw call and no
6+ * GPU-to-CPU round trip. (The host does read a reduction back occasionally, to
7+ * decide the colour scale.)
8+ *
9+ * The picture is a ray march. Each pixel casts one ray, intersects it with the
10+ * cube, and steps along it accumulating emission front to back: the pressure
11+ * goes through a diverging colormap about zero, and the opacity goes as a
12+ * power of |p|, so quiet regions are transparent and the wavefronts are what
13+ * you see. The medium is added as a grey emission proportional to how far the
14+ * local sound speed departs from the background, which shows a hard scatterer
15+ * as a distinct shape and a smooth one as a soft cloud without either needing
16+ * its own kind of drawing.
17+ *
18+ * Two honest limitations. Sampling is nearest-neighbour, not trilinear: the
19+ * field lives in a storage buffer rather than a filterable 3D texture, so
20+ * trilinear would be eight fetches per sample and the march takes tens of
21+ * millions of samples a frame. With the ray step set near the cell size the
22+ * difference is visible mainly as a faint stippling on strong wavefronts. And
23+ * the compositing is emission only, with no lighting and no shadowing, so what
24+ * is behind a strong feature is dimmed but never occluded correctly.
25+ *
26+ * A clip plane on x is provided because a volume render of a wavefield is
27+ * mostly the outside of a wavefield. Pulling the clip in is how you see the
28+ * interior, and it is the closest thing here to the flat sibling's picture.
29+ */
30+import type { ColormapFunc } from './colormaps.ts';
31+
32+const SHADER = `
33+struct View {
34+ eye: vec4f, // xyz: eye position, w: tan(fov/2)
35+ right: vec4f, // xyz: camera right, w: aspect ratio
36+ up: vec4f, // xyz: camera up, w: domain side L
37+ fwd: vec4f, // xyz: camera forward, w: pressure the colormap saturates at
38+ bg: vec4f, // rgb: background
39+ a: vec4f, // opacity, medium strength, reference speed, speed spread
40+ b: vec4f, // ray steps, clip plane x, grid side n, contrast exponent
41+ m: vec4f, // xyz: microphone position, w: whether to draw it
42+};
43+
44+@group(0) @binding(0) var<uniform> V: View;
45+@group(0) @binding(1) var<storage, read> p: array<f32>;
46+@group(0) @binding(2) var<storage, read> cs: array<f32>;
47+@group(0) @binding(3) var cmap: texture_2d<f32>;
48+@group(0) @binding(4) var samp: sampler;
49+
50+struct VSOut {
51+ @builtin(position) pos: vec4f,
52+ @location(0) ndc: vec2f,
53+};
54+
55+@vertex
56+fn vs(@builtin(vertex_index) vi: u32) -> VSOut {
57+ // One oversized triangle covering the viewport.
58+ var xy = array<vec2f, 3>(vec2f(-1.0, -3.0), vec2f(-1.0, 1.0), vec2f(3.0, 1.0));
59+ var out: VSOut;
60+ let q = xy[vi];
61+ out.pos = vec4f(q, 0.0, 1.0);
62+ out.ndc = q;
63+ return out;
64+}
65+
66+fn voxel(q: vec3f) -> u32 {
67+ let n = i32(V.b.z);
68+ let h = V.up.w / f32(n);
69+ let g = clamp(vec3i(floor((q + vec3f(0.5 * V.up.w)) / h)), vec3i(0), vec3i(n - 1));
70+ return u32(g.x + n * g.y + n * n * g.z);
71+}
72+
73+// A point on the cube's surface is on an edge when two of its three distances
74+// to the bounding planes vanish, so the test is on the median of the three.
75+fn edge(q: vec3f, half: f32, L: f32) -> f32 {
76+ let d = abs(abs(q) - vec3f(half));
77+ let lo = min(d.x, min(d.y, d.z));
78+ let hi = max(d.x, max(d.y, d.z));
79+ let mid = d.x + d.y + d.z - lo - hi;
80+ let w = 0.004 * L;
81+ return 1.0 - smoothstep(w, 2.0 * w, mid);
82+}
83+
84+@fragment
85+fn fs(in: VSOut) -> @location(0) vec4f {
86+ let L = V.up.w;
87+ let half = 0.5 * L;
88+ let eye = V.eye.xyz;
89+ let dir = normalize(V.fwd.xyz
90+ + in.ndc.x * V.right.w * V.eye.w * V.right.xyz
91+ + in.ndc.y * V.eye.w * V.up.xyz);
92+
93+ let inv = 1.0 / dir;
94+ let ta = (vec3f(-half) - eye) * inv;
95+ let tb = (vec3f(half) - eye) * inv;
96+ let lo = min(ta, tb);
97+ let hi = max(ta, tb);
98+ let t1 = min(min(hi.x, hi.y), hi.z);
99+ var t0 = max(max(lo.x, lo.y), lo.z);
100+
101+ var col = V.bg.rgb;
102+ if (t1 <= max(t0, 0.0)) { return vec4f(col, 1.0); }
103+ t0 = max(t0, 0.0);
104+
105+ let lineCol = vec3f(0.42, 0.47, 0.55);
106+ col = mix(col, lineCol, 0.5 * edge(eye + dir * t1, half, L));
107+
108+ let steps = i32(V.b.x);
109+ let dl = (t1 - t0) / f32(steps);
110+ // Opacity is quoted per cell, so a longer ray step is proportionally more
111+ // opaque and the picture does not change brightness with the quality knob.
112+ let unit = dl / (L / V.b.z);
113+ let scale = max(V.fwd.w, 1e-20);
114+ let grey = vec3f(0.55, 0.58, 0.63);
115+
116+ var acc = vec3f(0.0);
117+ var alpha = 0.0;
118+ for (var k = 0; k < steps; k = k + 1) {
119+ if (alpha > 0.995) { break; }
120+ let q = eye + dir * (t0 + (f32(k) + 0.5) * dl);
121+ if (q.x > V.b.y) { continue; }
122+ let i = voxel(q);
123+ let v = clamp(p[i] / scale, -1.0, 1.0);
124+ var a = pow(abs(v), V.b.w) * V.a.x * unit;
125+ var rgb = textureSampleLevel(cmap, samp, vec2f(0.5 + 0.5 * v, 0.5), 0.0).rgb;
126+ if (V.a.y > 0.0) {
127+ let m = clamp(abs(cs[i] - V.a.z) / max(V.a.w, 1e-20), 0.0, 1.0);
128+ let am = m * V.a.y * unit;
129+ let tot = a + am;
130+ if (tot > 1e-12) { rgb = (a * rgb + am * grey) / tot; }
131+ a = tot;
132+ }
133+ a = clamp(a, 0.0, 1.0);
134+ acc = acc + (1.0 - alpha) * a * rgb;
135+ alpha = alpha + (1.0 - alpha) * a;
136+ }
137+ col = acc + (1.0 - alpha) * col;
138+ col = mix(col, lineCol, 0.75 * edge(eye + dir * t0, half, L));
139+
140+ // The microphone, a dot in a ring at the ray's closest approach to it.
141+ // Drawn on top rather than composited into the march, like the 2D app's
142+ // white ring: a marker, not a thing in the scene.
143+ if (V.m.w > 0.0) {
144+ let toM = V.m.xyz - eye;
145+ let tm = dot(toM, dir);
146+ if (tm > 0.0) {
147+ let d = length(toM - tm * dir) / L;
148+ let dotm = 1.0 - smoothstep(0.004, 0.007, d);
149+ let ring = 1.0 - smoothstep(0.0015, 0.004, abs(d - 0.016));
150+ col = mix(col, vec3f(1.0), 0.9 * max(dotm, ring));
151+ }
152+ }
153+ return vec4f(col, 1.0);
154+}
155+`;
156+
157+/** Where the camera is looking from. Angles in radians, distance in metres. */
158+export interface Camera {
159+ az: number;
160+ el: number;
161+ dist: number;
162+}
163+
164+export interface DrawOptions {
165+ camera: Camera;
166+ /** Pressure the colormap saturates at, in both directions. */
167+ scale: number;
168+ opacity: number;
169+ contrast: number;
170+ /** Samples along each ray. */
171+ steps: number;
172+ /** Everything with x above this is not drawn, in metres. */
173+ clipX: number;
174+ /** Strength of the medium wash, 0 to turn it off. */
175+ medium: number;
176+ cref: number;
177+ cdev: number;
178+ /** Microphone position in metres, or null to not draw it. */
179+ mic?: { x: number; y: number; z: number } | null;
180+}
181+
182+const TAN_HALF_FOV = Math.tan((32 * Math.PI) / 360);
183+
184+const cross = (a: number[], b: number[]) => [
185+ a[1] * b[2] - a[2] * b[1],
186+ a[2] * b[0] - a[0] * b[2],
187+ a[0] * b[1] - a[1] * b[0],
188+];
189+const norm = (a: number[]) => {
190+ const m = Math.hypot(a[0], a[1], a[2]) || 1;
191+ return [a[0] / m, a[1] / m, a[2] / m];
192+};
193+
194+export class VolumeView {
195+ readonly canvas: HTMLCanvasElement;
196+
197+ #device: GPUDevice;
198+ #context: GPUCanvasContext;
199+ #pipeline: GPURenderPipeline;
200+ #layout: GPUBindGroupLayout;
201+ #uniform: GPUBuffer;
202+ #host = new ArrayBuffer(8 * 16);
203+ #sampler: GPUSampler;
204+ #cmapTexture: GPUTexture;
205+ #bindGroups: GPUBindGroup[] = [];
206+ #n = 0;
207+ #L = 1;
208+ /** Background, matched to the CSS so the cube sits on the page. */
209+ bg: [number, number, number] = [0.043, 0.055, 0.071];
210+
211+ constructor(device: GPUDevice, canvas: HTMLCanvasElement) {
212+ this.#device = device;
213+ this.canvas = canvas;
214+
215+ const context = canvas.getContext('webgpu');
216+ if (!context) throw new Error('this canvas has no WebGPU context');
217+ this.#context = context;
218+ const format = navigator.gpu.getPreferredCanvasFormat();
219+ context.configure({ device, format, alphaMode: 'opaque' });
220+
221+ this.#layout = device.createBindGroupLayout({
222+ label: 'volume-view',
223+ entries: [
224+ { binding: 0, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'uniform' } },
225+ { binding: 1, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'read-only-storage' } },
226+ { binding: 2, visibility: GPUShaderStage.FRAGMENT, buffer: { type: 'read-only-storage' } },
227+ { binding: 3, visibility: GPUShaderStage.FRAGMENT, texture: { sampleType: 'float' } },
228+ { binding: 4, visibility: GPUShaderStage.FRAGMENT, sampler: { type: 'filtering' } },
229+ ],
230+ });
231+
232+ const module = device.createShaderModule({ code: SHADER, label: 'volume-view' });
233+ this.#pipeline = device.createRenderPipeline({
234+ label: 'volume-view',
235+ layout: device.createPipelineLayout({ bindGroupLayouts: [this.#layout] }),
236+ vertex: { module, entryPoint: 'vs' },
237+ fragment: { module, entryPoint: 'fs', targets: [{ format }] },
238+ primitive: { topology: 'triangle-list' },
239+ });
240+
241+ this.#uniform = device.createBuffer({
242+ label: 'volume-view-uniform',
243+ size: this.#host.byteLength,
244+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
245+ });
246+ this.#sampler = device.createSampler({ magFilter: 'linear', minFilter: 'linear' });
247+ this.#cmapTexture = device.createTexture({
248+ label: 'colormap',
249+ size: [256, 1],
250+ format: 'rgba8unorm',
251+ usage: GPUTextureUsage.TEXTURE_BINDING | GPUTextureUsage.COPY_DST,
252+ });
253+ }
254+
255+ /**
256+ * Point the view at the buffers of a (new) simulation. Two pressure buffers
257+ * are given because the solver alternates between them; `draw` is told which
258+ * one is current.
259+ */
260+ setSource(pressures: [GPUBuffer, GPUBuffer], speed: GPUBuffer, n: number, L: number): void {
261+ this.#n = n;
262+ this.#L = L;
263+ this.#bindGroups = pressures.map((p) =>
264+ this.#device.createBindGroup({
265+ layout: this.#layout,
266+ entries: [
267+ { binding: 0, resource: { buffer: this.#uniform } },
268+ { binding: 1, resource: { buffer: p } },
269+ { binding: 2, resource: { buffer: speed } },
270+ { binding: 3, resource: this.#cmapTexture.createView() },
271+ { binding: 4, resource: this.#sampler },
272+ ],
273+ }),
274+ );
275+ }
276+
277+ setColormap(cmap: ColormapFunc): void {
278+ const data = new Uint8Array(256 * 4);
279+ for (let i = 0; i < 256; i++) {
280+ const [r, g, b] = cmap(i / 255);
281+ data[4 * i] = r;
282+ data[4 * i + 1] = g;
283+ data[4 * i + 2] = b;
284+ data[4 * i + 3] = 255;
285+ }
286+ this.#device.queue.writeTexture(
287+ { texture: this.#cmapTexture },
288+ data,
289+ { bytesPerRow: 256 * 4 },
290+ { width: 256, height: 1 },
291+ );
292+ }
293+
294+ /** Draw one frame from the pressure buffer with the given index. */
295+ draw(which: number, o: DrawOptions): void {
296+ const bg = this.#bindGroups[which];
297+ if (!bg) return;
298+
299+ const { az, el, dist } = o.camera;
300+ const ce = Math.cos(el);
301+ const toEye = [ce * Math.cos(az), ce * Math.sin(az), Math.sin(el)];
302+ const eye = toEye.map((v) => v * dist);
303+ const fwd = toEye.map((v) => -v);
304+ const right = norm(cross([0, 0, 1], fwd));
305+ const up = cross(fwd, right);
306+
307+ const rect = this.canvas.width / Math.max(this.canvas.height, 1);
308+ const f = new Float32Array(this.#host);
309+ f.set([eye[0], eye[1], eye[2], TAN_HALF_FOV], 0);
310+ f.set([right[0], right[1], right[2], rect], 4);
311+ f.set([up[0], up[1], up[2], this.#L], 8);
312+ f.set([fwd[0], fwd[1], fwd[2], o.scale], 12);
313+ f.set([this.bg[0], this.bg[1], this.bg[2], 1], 16);
314+ f.set([o.opacity, o.medium, o.cref, o.cdev], 20);
315+ f.set([o.steps, o.clipX, this.#n, o.contrast], 24);
316+ f.set([o.mic?.x ?? 0, o.mic?.y ?? 0, o.mic?.z ?? 0, o.mic ? 1 : 0], 28);
317+ this.#device.queue.writeBuffer(this.#uniform, 0, this.#host);
318+
319+ const enc = this.#device.createCommandEncoder({ label: 'volume-view' });
320+ const pass = enc.beginRenderPass({
321+ colorAttachments: [
322+ {
323+ view: this.#context.getCurrentTexture().createView(),
324+ clearValue: { r: this.bg[0], g: this.bg[1], b: this.bg[2], a: 1 },
325+ loadOp: 'clear',
326+ storeOp: 'store',
327+ },
328+ ],
329+ });
330+ pass.setPipeline(this.#pipeline);
331+ pass.setBindGroup(0, bg);
332+ pass.draw(3);
333+ pass.end();
334+ this.#device.queue.submit([enc.finish()]);
335+ }
336+
337+ /**
338+ * Match the canvas's backing store to its CSS size. Device pixel ratio is
339+ * ignored: the march is the whole cost of a frame and it scales with pixels,
340+ * so a retina display would pay four times over for a picture that is
341+ * already smooth.
342+ */
343+ resize(): void {
344+ const rect = this.canvas.getBoundingClientRect();
345+ const w = Math.max(1, Math.round(rect.width));
346+ const h = Math.max(1, Math.round(rect.height));
347+ if (this.canvas.width !== w || this.canvas.height !== h) {
348+ this.canvas.width = w;
349+ this.canvas.height = h;
350+ }
351+ }
352+
353+ destroy(): void {
354+ this.#uniform.destroy();
355+ this.#cmapTexture.destroy();
356+ }
357+}
src/scenes.tsadded+242−0View file
@@ -0,0 +1,242 @@
1+/**
2+ * Scenes: what the medium is.
3+ *
4+ * A scene returns two fields of position, the sound speed c (m/s) and the
5+ * absorption sig (1/s). Both are ordinary fields, which is what lets one
6+ * function describe both the scatterer and the open boundary: the absorbing
7+ * layer around the outside is the statement that the medium swallows sound out
8+ * there, and a scene is free to put absorption inside the domain too, making a
9+ * lossy scatterer.
10+ *
11+ * These are built on the CPU, once, whenever a parameter changes. At 128^3
12+ * that is two million evaluations and takes well under a second; at 192^3 it
13+ * is seven million and is noticeable. The flat sibling writes its scenes as
14+ * MATLAB and evaluates them through numbl's interpreter, which is the more
15+ * interesting arrangement and the obvious thing to bring over here later.
16+ *
17+ * Interfaces are smoothed over about a cell. A jump between two neighbouring
18+ * cells is not resolved by the grid: it scatters the grid's own staircase
19+ * rather than the shape that was asked for.
20+ */
21+import { coord, type Grid } from './grid.ts';
22+import { C_AIR, SPONGE_FRAC, SPONGE_MAX } from './units.ts';
23+
24+export interface SceneParam {
25+ key: string;
26+ label: string;
27+ min: number;
28+ max: number;
29+ step: number;
30+ value: number;
31+ unit?: string;
32+}
33+
34+export interface Medium {
35+ c: Float32Array<ArrayBuffer>;
36+ sig: Float32Array<ArrayBuffer>;
37+ cmin: number;
38+ cmax: number;
39+}
40+
41+export interface Scene {
42+ label: string;
43+ blurb: string;
44+ params: SceneParam[];
45+ /** Source settings this scene wants, applied when it is selected. */
46+ source?: Record<string, number>;
47+ /** Sound speed at (x, y, z), in m/s, and any absorption of its own. */
48+ medium(x: number, y: number, z: number, h: number, v: Record<string, number>): [number, number];
49+}
50+
51+/**
52+ * Absorption profile for an open boundary: zero in the interior, ramping up
53+ * quadratically over a layer of width w inside each face of the cube and
54+ * reaching SPONGE_MAX at the wall. This is what makes the finite grid stand in
55+ * for an unbounded medium.
56+ *
57+ * The ramp is gradual on purpose. An absorbing layer is itself an impedance
58+ * mismatch, so a sudden one reflects; spreading it over a couple of
59+ * wavelengths keeps that small. It is not a perfectly matched layer, and at
60+ * grazing incidence it does leak.
61+ */
62+function sponge(x: number, y: number, z: number, L: number): number {
63+ const w = SPONGE_FRAC * L;
64+ const dx = Math.max(0, w - (L / 2 - Math.abs(x))) / w;
65+ const dy = Math.max(0, w - (L / 2 - Math.abs(y))) / w;
66+ const dz = Math.max(0, w - (L / 2 - Math.abs(z))) / w;
67+ const d = Math.max(dx, Math.max(dy, dz));
68+ return SPONGE_MAX * d * d;
69+}
70+
71+/** A smoothed indicator: 1 well inside the surface, 0 well outside. */
72+const inside = (signedDistance: number, h: number): number =>
73+ 0.5 * (1 - Math.tanh(signedDistance / (1.5 * h)));
74+
75+/** Deterministic value noise on a coarse lattice, trilinearly interpolated. */
76+function makeNoise(m: number, seed: number): (u: number, v: number, w: number) => number {
77+ let s = (seed * 0x9e3779b9) >>> 0;
78+ const rand = () => {
79+ s = (s + 0x6d2b79f5) >>> 0;
80+ let t = s;
81+ t = Math.imul(t ^ (t >>> 15), t | 1);
82+ t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
83+ return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
84+ };
85+ const g = new Float64Array(m * m * m);
86+ for (let i = 0; i < g.length; i++) g[i] = 2 * rand() - 1;
87+ const at = (i: number, j: number, k: number) => {
88+ const w = (a: number) => ((a % m) + m) % m;
89+ return g[w(i) + m * w(j) + m * m * w(k)];
90+ };
91+ // u, v, w in [0, 1) over the whole domain.
92+ return (u, v, w) => {
93+ const fx = u * m;
94+ const fy = v * m;
95+ const fz = w * m;
96+ const i = Math.floor(fx);
97+ const j = Math.floor(fy);
98+ const k = Math.floor(fz);
99+ const sx = fx - i;
100+ const sy = fy - j;
101+ const sz = fz - k;
102+ // Smoothstep on each axis so the field is continuous in its derivative,
103+ // which keeps the medium resolvable by the grid.
104+ const ex = sx * sx * (3 - 2 * sx);
105+ const ey = sy * sy * (3 - 2 * sy);
106+ const ez = sz * sz * (3 - 2 * sz);
107+ const lerp = (a: number, b: number, t: number) => a + (b - a) * t;
108+ const c00 = lerp(at(i, j, k), at(i + 1, j, k), ex);
109+ const c10 = lerp(at(i, j + 1, k), at(i + 1, j + 1, k), ex);
110+ const c01 = lerp(at(i, j, k + 1), at(i + 1, j, k + 1), ex);
111+ const c11 = lerp(at(i, j + 1, k + 1), at(i + 1, j + 1, k + 1), ex);
112+ return lerp(lerp(c00, c10, ey), lerp(c01, c11, ey), ez);
113+ };
114+}
115+
116+/** Rebuilt whenever the random medium's seed or scale changes. */
117+let noiseCache: { key: string; f: (u: number, v: number, w: number) => number } | null = null;
118+
119+export const scenes: Record<string, Scene> = {
120+ sphere: {
121+ label: 'Sphere',
122+ blurb:
123+ 'One spherical scatterer in a uniform background: the reference case, and the ' +
124+ 'three-dimensional problem with a classical series solution. Speeds above 1 behave ' +
125+ 'nearly rigid, below 1 nearly pressure-release, and exactly 1 is no scatterer at all.',
126+ params: [
127+ { key: 'R', label: 'radius', min: 0.05, max: 0.5, step: 0.01, value: 0.25, unit: 'm' },
128+ { key: 'cin', label: 'speed ratio', min: 0.1, max: 4, step: 0.05, value: 2.5 },
129+ { key: 'absorb', label: 'absorption', min: 0, max: 6000, step: 100, value: 0, unit: '1/s' },
130+ ],
131+ medium(x, y, z, h, v) {
132+ const r = Math.sqrt(x * x + y * y + z * z);
133+ const q = inside(r - v.R, h);
134+ return [C_AIR * (1 + (v.cin - 1) * q), v.absorb * q];
135+ },
136+ },
137+
138+ pair: {
139+ label: 'Two spheres',
140+ blurb:
141+ 'Two identical scatterers side by side. Each one reradiates what it receives, ' +
142+ 'including what it receives from the other, so the pattern behind them is not the sum ' +
143+ 'of two single-sphere patterns.',
144+ params: [
145+ { key: 'R', label: 'radius', min: 0.05, max: 0.35, step: 0.01, value: 0.16, unit: 'm' },
146+ { key: 'sep', label: 'separation', min: 0.15, max: 1.0, step: 0.01, value: 0.5, unit: 'm' },
147+ { key: 'cin', label: 'speed ratio', min: 0.1, max: 4, step: 0.05, value: 3 },
148+ ],
149+ medium(x, y, z, h, v) {
150+ const d = v.sep / 2;
151+ const r1 = Math.sqrt(x * x + (y - d) * (y - d) + z * z);
152+ const r2 = Math.sqrt(x * x + (y + d) * (y + d) + z * z);
153+ const q = Math.max(inside(r1 - v.R, h), inside(r2 - v.R, h));
154+ return [C_AIR * (1 + (v.cin - 1) * q), 0];
155+ },
156+ },
157+
158+ aperture: {
159+ label: 'Aperture',
160+ blurb:
161+ 'A screen with a circular hole in it. Behind the hole is the three-dimensional ' +
162+ 'diffraction pattern, which a plane cut through the flat problem cannot show: the ' +
163+ 'Airy-like rings are a property of the circle, not of the slit.',
164+ params: [
165+ { key: 'a', label: 'hole radius', min: 0.04, max: 0.5, step: 0.01, value: 0.15, unit: 'm' },
166+ { key: 'th', label: 'screen thickness', min: 0.02, max: 0.2, step: 0.01, value: 0.06, unit: 'm' },
167+ { key: 'cw', label: 'screen speed', min: 0.05, max: 1, step: 0.05, value: 0.2 },
168+ ],
169+ medium(x, y, z, h, v) {
170+ // The screen is slow rather than fast. Reflection at an interface goes
171+ // as |c2 - c1|/(c2 + c1) at constant density, so c = 0.2 reflects about
172+ // as much as c = 5 would, but the timestep is set by the fastest speed
173+ // anywhere on the grid, so a slow screen is free and a fast one taxes
174+ // every step of the whole run. What a slow screen costs instead is
175+ // resolution inside itself, which its own absorption then swallows.
176+ const slab = inside(Math.abs(x) - v.th / 2, h);
177+ const rho = Math.sqrt(y * y + z * z);
178+ const hole = inside(rho - v.a, h);
179+ const q = slab * (1 - hole);
180+ return [C_AIR * (1 + (v.cw - 1) * q), 4000 * q];
181+ },
182+ },
183+
184+ random: {
185+ label: 'Random medium',
186+ blurb:
187+ 'Weak random structure everywhere, with no scatterer in particular. A pulse through ' +
188+ 'it arrives on time and then keeps arriving: multiple scattering turns the tail into a ' +
189+ 'coda, which is what a real inhomogeneous medium does.',
190+ params: [
191+ { key: 'amp', label: 'contrast', min: 0, max: 0.6, step: 0.01, value: 0.25 },
192+ { key: 'scale', label: 'blob size', min: 0.05, max: 0.6, step: 0.01, value: 0.18, unit: 'm' },
193+ { key: 'seed', label: 'seed', min: 1, max: 20, step: 1, value: 1 },
194+ ],
195+ medium(x, y, z, _h, v) {
196+ // `buildMedium` has rebuilt the lattice for this scale and seed before
197+ // calling us, so it is always here.
198+ const f = noiseCache!.f;
199+ const L = DOM.L;
200+ return [C_AIR * (1 + v.amp * f((x + L / 2) / L, (y + L / 2) / L, (z + L / 2) / L)), 0];
201+ },
202+ },
203+};
204+
205+/** The domain the scenes are currently being evaluated over. Set by `build`;
206+ * the random medium needs it to place its lattice. */
207+const DOM = { L: 2 };
208+
209+export const sceneKeys = Object.keys(scenes);
210+
211+/** Evaluate a scene over the grid, adding the absorbing layer. */
212+export function buildMedium(scene: Scene, grid: Grid, v: Record<string, number>): Medium {
213+ DOM.L = grid.L;
214+ if (scene === scenes.random) {
215+ const m = Math.max(2, Math.round(grid.L / Math.max(v.scale, 1e-3)));
216+ const key = `${m}:${v.seed}`;
217+ if (!noiseCache || noiseCache.key !== key) noiseCache = { key, f: makeNoise(m, v.seed) };
218+ }
219+
220+ const { n, npts, L, h } = grid;
221+ const c = new Float32Array(npts);
222+ const sig = new Float32Array(npts);
223+ let cmin = Infinity;
224+ let cmax = -Infinity;
225+ for (let iz = 0; iz < n; iz++) {
226+ const z = coord(iz, grid);
227+ for (let iy = 0; iy < n; iy++) {
228+ const y = coord(iy, grid);
229+ const row = n * iy + n * n * iz;
230+ for (let ix = 0; ix < n; ix++) {
231+ const x = coord(ix, grid);
232+ const [cv, sv] = scene.medium(x, y, z, h, v);
233+ const k = ix + row;
234+ c[k] = cv;
235+ sig[k] = sv + sponge(x, y, z, L);
236+ if (cv < cmin) cmin = cv;
237+ if (cv > cmax) cmax = cv;
238+ }
239+ }
240+ }
241+ return { c, sig, cmin, cmax };
242+}
src/sim.tsadded+533−0View file
@@ -0,0 +1,533 @@
1+/**
2+ * The solver: a second-order leapfrog for the acoustic wave equation on a
3+ * cubic grid,
4+ *
5+ * p_tt + 2*sig*p_t = c(x)^2 * lap(p) + s(x, t)
6+ *
7+ * Pressure only, at constant density, so the medium is two fields: the sound
8+ * speed c and the absorption sig, both supplied by the scene. Centring the
9+ * second time derivative and the damping on step n gives an explicit update,
10+ * which is one compute dispatch per timestep.
11+ *
12+ * Two things are worth explaining.
13+ *
14+ * **The update is in place.** A step reads p at its six neighbours but reads
15+ * the previous field pm only at its own index, so a thread may overwrite
16+ * pm[i] with the new value: no other thread will read it. That leaves two
17+ * pressure buffers instead of three, which matters at 192^3 where each one is
18+ * 28 MB, and it removes the buffer-to-buffer copy per step. The roles swap
19+ * every step, so which buffer holds the current field depends on the parity of
20+ * the step count; `pressure` reports it and the renderer follows.
21+ *
22+ * **Time is a uniform read at a dynamic offset.** A frame's worth of steps is
23+ * recorded into one command encoder, so nothing written between submits can
24+ * change inside it, and a clock uploaded per frame would stand still for the
25+ * whole batch. Rather than carrying time as a grid field (the flat sibling's
26+ * answer, forced there by its compiler), each step reads its own 64-byte slice
27+ * of one parameter buffer through a dynamic offset. The whole batch's
28+ * parameters are written in one call before the pass.
29+ */
30+import type { Grid } from './grid.ts';
31+
32+/** Bytes of parameters per step. Padded to the uniform dynamic-offset
33+ * alignment, which WebGPU guarantees to be at most 256. */
34+const PARAM_STRIDE = 256;
35+
36+/** Steps that fit in the parameter buffer, and so in one submit. */
37+export const MAX_STEPS_PER_FRAME = 64;
38+
39+/** Workgroups are reduced into this many partial maxima before readback. */
40+const REDUCE_GROUPS = 64;
41+
42+/** Samples the microphone trace can hold: 4 MB of f32, a few seconds of
43+ * audio at the rates the default timestep implies. */
44+export const MIC_CAPACITY = 1 << 20;
45+
46+const STEP_SHADER = `
47+struct Params {
48+ n: u32,
49+ h: f32,
50+ dt: f32,
51+ t: f32,
52+ L: f32,
53+ f: f32,
54+ t0: f32,
55+ tw: f32,
56+ cw: f32,
57+ x0: f32,
58+ y0: f32,
59+ z0: f32,
60+ w: f32,
61+ point: f32,
62+};
63+
64+@group(0) @binding(0) var<uniform> P: Params;
65+@group(1) @binding(0) var<storage, read> p: array<f32>;
66+@group(1) @binding(1) var<storage, read_write> pm: array<f32>;
67+@group(1) @binding(2) var<storage, read> cs: array<f32>;
68+@group(1) @binding(3) var<storage, read> sg: array<f32>;
69+
70+// Outside the domain the field is zero. The absorbing layer is meant to have
71+// swallowed the wave long before it reaches here.
72+fn at(ix: i32, iy: i32, iz: i32) -> f32 {
73+ let n = i32(P.n);
74+ if (ix < 0 || iy < 0 || iz < 0 || ix >= n || iy >= n || iz >= n) { return 0.0; }
75+ return p[u32(ix + n * iy + n * n * iz)];
76+}
77+
78+@compute @workgroup_size(8, 8, 4)
79+fn main(@builtin(global_invocation_id) gid: vec3u) {
80+ let n = P.n;
81+ if (gid.x >= n || gid.y >= n || gid.z >= n) { return; }
82+ let i = gid.x + n * gid.y + n * n * gid.z;
83+ let ix = i32(gid.x);
84+ let iy = i32(gid.y);
85+ let iz = i32(gid.z);
86+
87+ let pc = p[i];
88+ let lap = (at(ix + 1, iy, iz) + at(ix - 1, iy, iz)
89+ + at(ix, iy + 1, iz) + at(ix, iy - 1, iz)
90+ + at(ix, iy, iz + 1) + at(ix, iy, iz - 1)
91+ - 6.0 * pc) / (P.h * P.h);
92+
93+ let x = -0.5 * P.L + (f32(ix) + 0.5) * P.h;
94+ let y = -0.5 * P.L + (f32(iy) + 0.5) * P.h;
95+ let z = -0.5 * P.L + (f32(iz) + 0.5) * P.h;
96+
97+ // The source. \`cw\` blends between a Gaussian pulse (0) and a wave that
98+ // turns on smoothly and stays on (1); \`point\` blends between a planar
99+ // source spanning the grid in y and z, whose far field is a plane wave, and
100+ // a point source at (x0, y0, z0).
101+ //
102+ // The om^2 is a choice of units, not a physical amplitude: a body force of
103+ // fixed strength drives a response falling off as 1/om^2, so without it the
104+ // field would shrink every time the frequency slider went up.
105+ let u = (P.t - P.t0) / P.tw;
106+ let env = (1.0 - P.cw) * exp(-u * u) + P.cw * 0.5 * (1.0 + tanh(u));
107+ let gx = (x - P.x0) / P.w;
108+ let gy = P.point * (y - P.y0) / P.w;
109+ let gz = P.point * (z - P.z0) / P.w;
110+ let om = 6.283185307179586 * P.f;
111+ let s = om * om * env * sin(om * (P.t - P.t0)) * exp(-(gx * gx + gy * gy + gz * gz));
112+
113+ // One step. The damping is what the absorbing layer acts through: sig is
114+ // zero over the interior, so there this is the plain leapfrog update.
115+ let sd = sg[i] * P.dt;
116+ let cd = cs[i] * P.dt;
117+ pm[i] = (2.0 * pc - (1.0 - sd) * pm[i] + cd * cd * lap + P.dt * P.dt * s) / (1.0 + sd);
118+}
119+`;
120+
121+// The microphone: the pressure at one grid point, appended to a trace by a
122+// one-thread dispatch that runs after each step inside the same command
123+// stream. The obvious implementation — reading the field back and picking out
124+// one number — costs a GPU-to-CPU round trip per step and would be slower
125+// than the step; this way the whole trace comes back once, when there is
126+// something to play. One thread and in-order dispatches mean the plain
127+// read-modify-write on the counter is safe.
128+const MIC_SHADER = `
129+struct MicInfo { probe: u32, cap: u32 };
130+@group(0) @binding(0) var<uniform> M: MicInfo;
131+@group(0) @binding(1) var<storage, read> p: array<f32>;
132+@group(0) @binding(2) var<storage, read_write> trace: array<f32>;
133+@group(0) @binding(3) var<storage, read_write> count: array<u32>;
134+
135+@compute @workgroup_size(1)
136+fn main() {
137+ let k = count[0];
138+ if (k < M.cap) {
139+ trace[k] = p[M.probe];
140+ count[0] = k + 1u;
141+ }
142+}
143+`;
144+
145+const REDUCE_SHADER = `
146+@group(0) @binding(0) var<uniform> npts: u32;
147+@group(0) @binding(1) var<storage, read> p: array<f32>;
148+@group(0) @binding(2) var<storage, read_write> out: array<f32>;
149+
150+var<workgroup> sh: array<f32, 256>;
151+
152+// Both loops are given uniform trip counts on purpose: a workgroupBarrier may
153+// only be reached in uniform control flow, and a loop whose exit depends on
154+// the thread index taints everything after it.
155+@compute @workgroup_size(256)
156+fn main(@builtin(global_invocation_id) gid: vec3u,
157+ @builtin(local_invocation_id) lid: vec3u,
158+ @builtin(workgroup_id) wid: vec3u) {
159+ let stride = 256u * ${REDUCE_GROUPS}u;
160+ let per = (npts + stride - 1u) / stride;
161+ var m = 0.0;
162+ for (var k = 0u; k < per; k = k + 1u) {
163+ let i = gid.x + k * stride;
164+ if (i < npts) { m = max(m, abs(p[i])); }
165+ }
166+ sh[lid.x] = m;
167+ workgroupBarrier();
168+ for (var s = 128u; s > 0u; s = s >> 1u) {
169+ if (lid.x < s) { sh[lid.x] = max(sh[lid.x], sh[lid.x + s]); }
170+ workgroupBarrier();
171+ }
172+ if (lid.x == 0u) { out[wid.x] = sh[0]; }
173+}
174+`;
175+
176+/** Everything the source term needs, in SI. */
177+export interface SourceParams {
178+ f: number;
179+ cycles: number;
180+ cw: number;
181+ point: number;
182+ x0: number;
183+ y0: number;
184+ z0: number;
185+ w: number;
186+}
187+
188+export class Sim {
189+ readonly grid: Grid;
190+ /** Model time, seconds. */
191+ t = 0;
192+ steps = 0;
193+ dt: number;
194+ source: SourceParams;
195+
196+ #device: GPUDevice;
197+ #pa: GPUBuffer;
198+ #pb: GPUBuffer;
199+ #c: GPUBuffer;
200+ #sig: GPUBuffer;
201+ #params: GPUBuffer;
202+ #paramsHost = new ArrayBuffer(MAX_STEPS_PER_FRAME * PARAM_STRIDE);
203+ #paramsBG: GPUBindGroup;
204+ #fieldBG: [GPUBindGroup, GPUBindGroup];
205+ #pipeline: GPUComputePipeline;
206+ /** Index of the field bind group to use for the next step. It is also which
207+ * buffer currently holds the field: 0 means `pa`. */
208+ #next = 0;
209+
210+ #reducePipeline: GPUComputePipeline;
211+ #reduceBG: [GPUBindGroup, GPUBindGroup];
212+ #partials: GPUBuffer;
213+ #readback: GPUBuffer;
214+ #reading = false;
215+
216+ /** Samples recorded so far. A host-side mirror of the GPU counter: both
217+ * add one per step until the capacity, so they agree exactly. */
218+ recorded = 0;
219+ #micPipeline: GPUComputePipeline;
220+ #micBG: [GPUBindGroup, GPUBindGroup];
221+ #micUniform: GPUBuffer;
222+ #trace: GPUBuffer;
223+ #micCount: GPUBuffer;
224+ #micRB: GPUBuffer;
225+ #micReading = false;
226+
227+ constructor(
228+ device: GPUDevice,
229+ grid: Grid,
230+ medium: { c: Float32Array<ArrayBuffer>; sig: Float32Array<ArrayBuffer> },
231+ source: SourceParams,
232+ dt: number,
233+ ) {
234+ this.#device = device;
235+ this.grid = grid;
236+ this.source = source;
237+ this.dt = dt;
238+
239+ const bytes = grid.npts * 4;
240+ const field = (label: string) =>
241+ device.createBuffer({
242+ label,
243+ size: bytes,
244+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
245+ });
246+ this.#pa = field('p-a');
247+ this.#pb = field('p-b');
248+ this.#c = field('speed');
249+ this.#sig = field('absorption');
250+ device.queue.writeBuffer(this.#c, 0, medium.c);
251+ device.queue.writeBuffer(this.#sig, 0, medium.sig);
252+
253+ this.#params = device.createBuffer({
254+ label: 'step-params',
255+ size: this.#paramsHost.byteLength,
256+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
257+ });
258+
259+ const module = device.createShaderModule({ code: STEP_SHADER, label: 'leapfrog3d' });
260+ const paramsLayout = device.createBindGroupLayout({
261+ entries: [
262+ {
263+ binding: 0,
264+ visibility: GPUShaderStage.COMPUTE,
265+ buffer: { type: 'uniform', hasDynamicOffset: true, minBindingSize: 64 },
266+ },
267+ ],
268+ });
269+ const fieldLayout = device.createBindGroupLayout({
270+ entries: [
271+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
272+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
273+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
274+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
275+ ],
276+ });
277+ this.#pipeline = device.createComputePipeline({
278+ label: 'leapfrog3d',
279+ layout: device.createPipelineLayout({ bindGroupLayouts: [paramsLayout, fieldLayout] }),
280+ compute: { module, entryPoint: 'main' },
281+ });
282+ this.#paramsBG = device.createBindGroup({
283+ layout: paramsLayout,
284+ entries: [{ binding: 0, resource: { buffer: this.#params, size: 64 } }],
285+ });
286+ const pair = (read: GPUBuffer, write: GPUBuffer) =>
287+ device.createBindGroup({
288+ layout: fieldLayout,
289+ entries: [
290+ { binding: 0, resource: { buffer: read } },
291+ { binding: 1, resource: { buffer: write } },
292+ { binding: 2, resource: { buffer: this.#c } },
293+ { binding: 3, resource: { buffer: this.#sig } },
294+ ],
295+ });
296+ this.#fieldBG = [pair(this.#pa, this.#pb), pair(this.#pb, this.#pa)];
297+
298+ // The colour scale wants max |p| over the whole field, which at 128^3 is
299+ // 8 MB to read back. A workgroup reduction turns it into 64 floats first.
300+ const nptsBuf = device.createBuffer({
301+ size: 4,
302+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
303+ });
304+ device.queue.writeBuffer(nptsBuf, 0, new Uint32Array([grid.npts]));
305+ this.#partials = device.createBuffer({
306+ size: REDUCE_GROUPS * 4,
307+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
308+ });
309+ this.#readback = device.createBuffer({
310+ size: REDUCE_GROUPS * 4,
311+ usage: GPUBufferUsage.COPY_DST | GPUBufferUsage.MAP_READ,
312+ });
313+ const reduceLayout = device.createBindGroupLayout({
314+ entries: [
315+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform' } },
316+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
317+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
318+ ],
319+ });
320+ this.#reducePipeline = device.createComputePipeline({
321+ label: 'maxabs',
322+ layout: device.createPipelineLayout({ bindGroupLayouts: [reduceLayout] }),
323+ compute: {
324+ module: device.createShaderModule({ code: REDUCE_SHADER, label: 'maxabs' }),
325+ entryPoint: 'main',
326+ },
327+ });
328+ const red = (p: GPUBuffer) =>
329+ device.createBindGroup({
330+ layout: reduceLayout,
331+ entries: [
332+ { binding: 0, resource: { buffer: nptsBuf } },
333+ { binding: 1, resource: { buffer: p } },
334+ { binding: 2, resource: { buffer: this.#partials } },
335+ ],
336+ });
337+ this.#reduceBG = [red(this.#pa), red(this.#pb)];
338+
339+ // The microphone. All buffers start zeroed, so the counter needs no init.
340+ this.#micUniform = device.createBuffer({
341+ label: 'mic-info',
342+ size: 8,
343+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
344+ });
345+ this.#trace = device.createBuffer({
346+ label: 'mic-trace',
347+ size: MIC_CAPACITY * 4,
348+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
349+ });
350+ this.#micCount = device.createBuffer({
351+ label: 'mic-count',
352+ size: 4,
353+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
354+ });
355+ this.#micRB = device.createBuffer({
356+ label: 'mic-readback',
357+ size: MIC_CAPACITY * 4,
358+ usage: GPUBufferUsage.COPY_DST | GPUBufferUsage.MAP_READ,
359+ });
360+ const micLayout = device.createBindGroupLayout({
361+ entries: [
362+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform' } },
363+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
364+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
365+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
366+ ],
367+ });
368+ this.#micPipeline = device.createComputePipeline({
369+ label: 'microphone',
370+ layout: device.createPipelineLayout({ bindGroupLayouts: [micLayout] }),
371+ compute: {
372+ module: device.createShaderModule({ code: MIC_SHADER, label: 'microphone' }),
373+ entryPoint: 'main',
374+ },
375+ });
376+ const micOn = (p: GPUBuffer) =>
377+ device.createBindGroup({
378+ layout: micLayout,
379+ entries: [
380+ { binding: 0, resource: { buffer: this.#micUniform } },
381+ { binding: 1, resource: { buffer: p } },
382+ { binding: 2, resource: { buffer: this.#trace } },
383+ { binding: 3, resource: { buffer: this.#micCount } },
384+ ],
385+ });
386+ this.#micBG = [micOn(this.#pa), micOn(this.#pb)];
387+ this.setProbe(0, 0, 0);
388+ }
389+
390+ /** Both pressure buffers; `pressureIndex` says which holds the field now. */
391+ get pressures(): [GPUBuffer, GPUBuffer] {
392+ return [this.#pa, this.#pb];
393+ }
394+ get pressureIndex(): number {
395+ return this.#next;
396+ }
397+ get speed(): GPUBuffer {
398+ return this.#c;
399+ }
400+
401+ /** Point the microphone at a grid cell. Does not restart the trace: moving
402+ * the microphone during a run is a microphone that moved. */
403+ setProbe(ix: number, iy: number, iz: number): void {
404+ const n = this.grid.n;
405+ const cl = (i: number) => Math.min(n - 1, Math.max(0, i));
406+ const probe = cl(ix) + n * cl(iy) + n * n * cl(iz);
407+ this.#device.queue.writeBuffer(this.#micUniform, 0, new Uint32Array([probe, MIC_CAPACITY]));
408+ }
409+
410+ /** Start the recording over. Called on restart, and whenever dt changes:
411+ * the trace is one sample per step, and two timesteps would be two sample
412+ * rates in one buffer. */
413+ resetTrace(): void {
414+ this.#device.queue.writeBuffer(this.#micCount, 0, new Uint32Array([0]));
415+ this.recorded = 0;
416+ }
417+
418+ /** The recorded trace, back from the GPU. Null if a read is in flight. */
419+ async readTrace(): Promise<Float32Array | null> {
420+ if (this.#micReading) return null;
421+ if (this.recorded === 0) return new Float32Array(0);
422+ this.#micReading = true;
423+ try {
424+ const bytes = this.recorded * 4;
425+ const enc = this.#device.createCommandEncoder({ label: 'mic-read' });
426+ enc.copyBufferToBuffer(this.#trace, 0, this.#micRB, 0, bytes);
427+ this.#device.queue.submit([enc.finish()]);
428+ await this.#micRB.mapAsync(GPUMapMode.READ, 0, bytes);
429+ const out = new Float32Array(this.#micRB.getMappedRange(0, bytes).slice(0));
430+ this.#micRB.unmap();
431+ return out;
432+ } finally {
433+ this.#micReading = false;
434+ }
435+ }
436+
437+ /** Back to a silent grid at t = 0. */
438+ restart(): void {
439+ const enc = this.#device.createCommandEncoder();
440+ enc.clearBuffer(this.#pa);
441+ enc.clearBuffer(this.#pb);
442+ this.#device.queue.submit([enc.finish()]);
443+ this.t = 0;
444+ this.steps = 0;
445+ this.#next = 0;
446+ this.resetTrace();
447+ }
448+
449+ /** Take `n` timesteps, all in one submit. */
450+ run(n: number): void {
451+ const count = Math.min(n, MAX_STEPS_PER_FRAME);
452+ const g = this.grid;
453+ const s = this.source;
454+ // A pulse `cycles` long, delayed enough that it starts near zero.
455+ const tw = s.cycles / Math.max(s.f, 1e-6);
456+ const t0 = 2.5 * tw;
457+ for (let k = 0; k < count; k++) {
458+ const off = k * PARAM_STRIDE;
459+ const u32 = new Uint32Array(this.#paramsHost, off, 1);
460+ const f32 = new Float32Array(this.#paramsHost, off, 16);
461+ u32[0] = g.n;
462+ f32[1] = g.h;
463+ f32[2] = this.dt;
464+ f32[3] = this.t + k * this.dt;
465+ f32[4] = g.L;
466+ f32[5] = s.f;
467+ f32[6] = t0;
468+ f32[7] = tw;
469+ f32[8] = s.cw;
470+ f32[9] = s.x0;
471+ f32[10] = s.y0;
472+ f32[11] = s.z0;
473+ f32[12] = s.w;
474+ f32[13] = s.point;
475+ }
476+ this.#device.queue.writeBuffer(this.#params, 0, this.#paramsHost, 0, count * PARAM_STRIDE);
477+
478+ const wg = [g.n / 8, g.n / 8, g.n / 4] as const;
479+ const enc = this.#device.createCommandEncoder({ label: 'steps' });
480+ const pass = enc.beginComputePass();
481+ for (let k = 0; k < count; k++) {
482+ pass.setPipeline(this.#pipeline);
483+ pass.setBindGroup(0, this.#paramsBG, [k * PARAM_STRIDE]);
484+ pass.setBindGroup(1, this.#fieldBG[this.#next]);
485+ pass.dispatchWorkgroups(wg[0], wg[1], wg[2]);
486+ this.#next ^= 1;
487+ // Record the field this step just wrote, which the toggle now points at.
488+ pass.setPipeline(this.#micPipeline);
489+ pass.setBindGroup(0, this.#micBG[this.#next]);
490+ pass.dispatchWorkgroups(1);
491+ }
492+ pass.end();
493+ this.#device.queue.submit([enc.finish()]);
494+ this.t += count * this.dt;
495+ this.steps += count;
496+ this.recorded = Math.min(this.recorded + count, MIC_CAPACITY);
497+ }
498+
499+ /**
500+ * Largest |p| anywhere, for the colour scale. Asynchronous and
501+ * self-throttling: while one request is in flight further ones return null
502+ * rather than queueing up.
503+ */
504+ async peak(): Promise<number | null> {
505+ if (this.#reading) return null;
506+ this.#reading = true;
507+ try {
508+ const enc = this.#device.createCommandEncoder({ label: 'peak' });
509+ const pass = enc.beginComputePass();
510+ pass.setPipeline(this.#reducePipeline);
511+ pass.setBindGroup(0, this.#reduceBG[this.pressureIndex]);
512+ pass.dispatchWorkgroups(REDUCE_GROUPS);
513+ pass.end();
514+ enc.copyBufferToBuffer(this.#partials, 0, this.#readback, 0, REDUCE_GROUPS * 4);
515+ this.#device.queue.submit([enc.finish()]);
516+ await this.#readback.mapAsync(GPUMapMode.READ);
517+ const v = new Float32Array(this.#readback.getMappedRange().slice(0));
518+ this.#readback.unmap();
519+ let m = 0;
520+ for (const x of v) m = Math.max(m, x);
521+ return m;
522+ } finally {
523+ this.#reading = false;
524+ }
525+ }
526+
527+ destroy(): void {
528+ const own = [this.#pa, this.#pb, this.#c, this.#sig, this.#params, this.#partials];
529+ for (const b of [...own, this.#micUniform, this.#trace, this.#micCount]) b.destroy();
530+ // Mapping is asynchronous; destroying a buffer with a pending map is an
531+ // error, so leave the readback buffers to be collected.
532+ }
533+}
src/units.tsadded+40−0View file
@@ -0,0 +1,40 @@
1+/**
2+ * Everything in SI: metres, seconds, hertz, metres per second.
3+ *
4+ * The domain is deliberately small. A three-dimensional grid costs the cube of
5+ * its resolution, so where the flat sibling can afford 512 points across ten
6+ * metres this one gets 128 across two, and the number that decides whether
7+ * what is on screen is physics or grid dispersion is the same either way: how
8+ * many cells fit in a wavelength. At 128 points across 2 m a 1.5 kHz tone is
9+ * about fifteen cells per wavelength, which is comfortable; raise the
10+ * frequency and the app says when it stops being so.
11+ */
12+
13+/** Speed of sound in air at about 20 °C, m/s. */
14+export const C_AIR = 343;
15+
16+/** Side of the cubic domain, metres. Room-corner sized rather than
17+ * hall sized, for the reason above. */
18+export const DOMAIN = 2;
19+
20+/** Fraction of the domain given to the absorbing layer at each face. */
21+export const SPONGE_FRAC = 0.15;
22+
23+/** Absorption rate the sponge reaches at the wall, inverse seconds. */
24+export const SPONGE_MAX = 9000;
25+
26+/** Cells per wavelength below which what is on screen is as much grid
27+ * dispersion as it is sound. */
28+export const POOR_RESOLUTION = 8;
29+
30+/** A length in metres, written the way a person would say it. */
31+export const fmtLength = (m: number): string =>
32+ Math.abs(m) < 1 ? `${(1000 * m).toPrecision(3)} mm` : `${m.toPrecision(3)} m`;
33+
34+/** A duration in seconds, likewise. */
35+export const fmtTime = (s: number): string => {
36+ const a = Math.abs(s);
37+ if (a > 0 && a < 1e-3) return `${(1e6 * s).toPrecision(3)} µs`;
38+ if (a < 1) return `${(1e3 * s).toPrecision(3)} ms`;
39+ return `${s.toPrecision(3)} s`;
40+};
tsconfig.jsonadded+15−0View file
@@ -0,0 +1,15 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "module": "ESNext",
5+ "moduleResolution": "bundler",
6+ "lib": ["ES2022", "DOM", "DOM.Iterable"],
7+ "types": ["@webgpu/types", "node"],
8+ "strict": true,
9+ "noEmit": true,
10+ "allowImportingTsExtensions": true,
11+ "verbatimModuleSyntax": true,
12+ "skipLibCheck": true
13+ },
14+ "include": ["src"]
15+}
vite.config.tsadded+8−0View file
@@ -0,0 +1,8 @@
1+import { defineConfig } from 'vite';
2+
3+export default defineConfig({
4+ base: './',
5+ build: {
6+ target: 'es2022',
7+ },
8+});