1# acoustic-scattering-3d
3Sound scattering off obstacles in three dimensions, solved live in the
4browser. A second-order leapfrog for the acoustic wave equation runs on a
5cubic grid as WebGPU compute shaders, and the picture is a volume ray march
6that reads the pressure buffer the solver just wrote, with no readback in the
7display path.
9This is the volumetric sibling of
10[acoustic-scattering-2d](https://github.com/concept-collection/acoustic-scattering-2d).
11That project's point is its compiler: the solver and the medium are MATLAB,
12compiled to compute kernels by [numbl](https://numbl.org) live in the page.
13This one starts simpler, with the solver and the scenes written directly in
14WGSL and TypeScript, because the third dimension brings enough new problems of
15its own: memory (a field is n³ floats, so four fields at 192³ are 113 MB),
16a stricter stability limit, and the question of how to look at a wavefield you
17are standing outside of. Bringing the editable-MATLAB arrangement over is the
18obvious next step.
20## The equation
22Pressure only, at constant density:
24```
25p_tt + 2*sig*p_t = c(x)^2 * lap(p) + s(x, t)
26```
28`c` is the sound speed and `sig` the absorption rate, both fields of position
29that the scene defines. Since the density is constant, the impedance ratio
30across an interface is just the speed ratio, so a scatterer much faster than
31its background behaves nearly rigid (sound-hard) and one much slower nearly
32pressure-release (sound-soft). What this formulation cannot do is set
33impedance and speed independently, which needs a variable-density divergence
34form and a second field.
36Centring both the second time derivative and the damping on step *n* gives an
37explicit 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
39the fastest speed anywhere in the medium, so a fast scatterer slows the whole
40run down.
42## How a step runs
44One timestep is one compute dispatch, and the update is in place: a step reads
45the current field `p` at its six neighbours but the previous field `pm` only
46at its own index, so each thread may overwrite `pm[i]` with the new value.
47That leaves two pressure buffers instead of three and no copies per step,
48which matters when a buffer is 28 MB. The two buffers swap roles every step,
49and the renderer is told which one is current.
51A frame's worth of steps (up to 64) is recorded into one command encoder, so
52the source term cannot read anything uploaded between frames. Each step
53instead reads its own slice of a parameter buffer through a dynamic uniform
54offset; the whole batch's parameters, including each step's time, are written
55in one call before the pass. This is the same problem the 2D app solves by
56carrying time as a grid field. The dynamic-offset answer is cheaper (no field,
57no kernel), and is available here because the solver is hand-written rather
58than compiled from MATLAB that only knows about fields.
60## Looking at it
62A two-dimensional field is its own picture; a three-dimensional one is not,
63and every way of drawing it hides something. The app ray-marches the volume:
64each pixel casts a ray through the cube and accumulates colour front to back,
65with the pressure through a diverging colormap about zero and opacity rising
66as a power of |p|, so quiet regions are transparent and wavefronts are what
67you see. The medium is blended in as a grey cloud so the scatterer is visible
68inside the field. Drag to orbit, scroll to zoom.
70Since a volume render of a wavefield is mostly the outside of the wavefield,
71there is a clip plane on x: pull it in and the interior is exposed, which is
72the closest thing here to the 2D picture.
74Two rendering shortcuts are worth knowing. Sampling along the ray is
75nearest-neighbour, because the field lives in a storage buffer rather than a
76filterable 3D texture; with the ray step near the cell size this shows mainly
77as faint stippling on strong fronts. And the compositing is emission only,
78with no lighting, so depth ordering dims what is behind a strong feature but
79nothing casts a shadow.
81## Listening to it
83There is a microphone: the pressure at one grid point, sampled every timestep.
84It is written on the GPU by a one-thread dispatch after each step and appended
85to a trace buffer, because the obvious implementation, reading the field back
86and picking out one number, costs a GPU-to-CPU round trip per step. The whole
87trace comes back once, when there is something to play.
89Everything is SI, so playback is real time: the trace's native rate is 1/dt,
90around 190 kHz at the defaults, which is above what Web Audio will accept, so
91it is resampled to 48 kHz on the way out. The content is band-limited far
92below either rate. The recording is normalized before playback, which
93discards absolute amplitude: that is what the colour scale is for. It restarts
94whenever the run does, and whenever the timestep changes, since a trace is one
95sample per step and two timesteps would be two sample rates in one buffer.
97As in 2D, a pulse is short: a few cycles at an audible frequency is a few
98milliseconds however long the simulation runs. For sustained sound, turn
99`continuous` up and let the run fill some seconds; the trace holds about five
100seconds at the default timestep.
102## What it is honest about
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.
121## Scenes
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. |
130The aperture screen is *slower* than the background rather than faster.
131Reflection at an interface goes as `|c2 - c1|/(c2 + c1)` at constant density,
132so a screen at c = 0.2 reflects about as much as one at c = 5 would, but the
133timestep is set by the fastest speed anywhere on the grid: a slow screen is
134free while a fast one taxes every step of the whole run. What the slow screen
135costs instead is resolution inside itself, which its own absorption swallows.
137## Tests
139```
140npm run smoke # the page itself, in headless Chrome
141```
143The browser check drives the real page: it loads it, waits for the solver to
144take steps and the frame loop to turn, halves the timestep through the slider
145and requires dt to follow, checks the microphone is recording, and swaps the
146scene. It says nothing about whether the picture is right, which is what eyes
147are for.
149## Development
151```
152npm install
153npm run dev
154```
156No dependencies beyond the build tooling: the solver and renderer are plain
157TypeScript and WGSL.
159## License
161Apache-2.0