acoustic-scattering-3d#

Sound scattering off obstacles in three dimensions, solved live in the browser. A second-order leapfrog for the acoustic wave equation runs on a cubic grid as WebGPU compute shaders, and the picture is a volume ray march that reads the pressure buffer the solver just wrote, with no readback in the display path.

This is the volumetric sibling of acoustic-scattering-2d. That project's point is its compiler: the solver and the medium are MATLAB, compiled to compute kernels by numbl live in the page. This one starts simpler, with the solver and the scenes written directly in WGSL and TypeScript, because the third dimension brings enough new problems of its own: memory (a field is n³ floats, so four fields at 192³ are 113 MB), a stricter stability limit, and the question of how to look at a wavefield you are standing outside of. Bringing the editable-MATLAB arrangement over is the obvious next step.

The equation#

Pressure only, at constant density:

p_tt + 2*sig*p_t = c(x)^2 * lap(p) + s(x, t)

c is the sound speed and sig the absorption rate, both fields of position that the scene defines. Since the density is constant, the impedance ratio across an interface is just the speed ratio, so a scatterer much faster than its background behaves nearly rigid (sound-hard) and one much slower nearly pressure-release (sound-soft). What this formulation cannot do is set impedance and speed independently, which needs a variable-density divergence form and a second field.

Centring both the second time derivative and the damping on step n gives an explicit update. The 7-point Laplacian makes the stability condition c*dt/h <= 1/sqrt(3), stricter than the 2D 1/sqrt(2); the app sets dt from the fastest speed anywhere in the medium, so a fast scatterer slows the whole run down.

How a step runs#

One timestep is one compute dispatch, and the update is in place: a step reads the current field p at its six neighbours but the previous field pm only at its own index, so each thread may overwrite pm[i] with the new value. That leaves two pressure buffers instead of three and no copies per step, which matters when a buffer is 28 MB. The two buffers swap roles every step, and the renderer is told which one is current.

A frame's worth of steps (up to 64) is recorded into one command encoder, so the source term cannot read anything uploaded between frames. Each step instead reads its own slice of a parameter buffer through a dynamic uniform offset; the whole batch's parameters, including each step's time, are written in one call before the pass. This is the same problem the 2D app solves by carrying time as a grid field. The dynamic-offset answer is cheaper (no field, no kernel), and is available here because the solver is hand-written rather than compiled from MATLAB that only knows about fields.

Looking at it#

A two-dimensional field is its own picture; a three-dimensional one is not, and every way of drawing it hides something. The app ray-marches the volume: each pixel casts a ray through the cube and accumulates colour front to back, with the pressure through a diverging colormap about zero and opacity rising as a power of |p|, so quiet regions are transparent and wavefronts are what you see. The medium is blended in as a grey cloud so the scatterer is visible inside the field. Drag to orbit, scroll to zoom.

Since a volume render of a wavefield is mostly the outside of the wavefield, there is a clip plane on x: pull it in and the interior is exposed, which is the closest thing here to the 2D picture.

Two rendering shortcuts are worth knowing. Sampling along the ray is nearest-neighbour, because the field lives in a storage buffer rather than a filterable 3D texture; with the ray step near the cell size this shows mainly as faint stippling on strong fronts. And the compositing is emission only, with no lighting, so depth ordering dims what is behind a strong feature but nothing casts a shadow.

Listening to it#

There is a microphone: the pressure at one grid point, sampled every timestep. It is written on the GPU by a one-thread dispatch after each step and appended to a trace buffer, because the obvious implementation, reading the field back and picking out one number, costs a GPU-to-CPU round trip per step. The whole trace comes back once, when there is something to play.

Everything is SI, so playback is real time: the trace's native rate is 1/dt, around 190 kHz at the defaults, which is above what Web Audio will accept, so it is resampled to 48 kHz on the way out. The content is band-limited far below either rate. The recording is normalized before playback, which discards absolute amplitude: that is what the colour scale is for. It restarts whenever the run does, and whenever the timestep changes, since a trace is one sample per step and two timesteps would be two sample rates in one buffer.

As in 2D, a pulse is short: a few cycles at an audible frequency is a few milliseconds however long the simulation runs. For sustained sound, turn continuous up and let the run fill some seconds; the trace holds about five seconds at the default timestep.

What it is honest about#

• Resolution is the whole game. A grid solver resolves a wavelength with some number of cells, and in 3D cells cost their cube. At 128³ over a 2 m domain, 1.5 kHz has about 11 cells per wavelength; the stats line turns the number orange when it drops below 8, at which point what is on screen is as much grid dispersion as sound. There is no lap4 here yet; the 2D project shows what a fourth-order stencil buys.
• The absorbing layer is a sponge, not a PML. Absorption ramps up quadratically over the outer 15% of each face. An absorbing layer is itself an impedance mismatch, so it reflects a little; a perfectly matched layer would do better at the cost of extra fields.
• Single precision. WebGPU has no f64. The stencil differences lose a few digits to cancellation, well below the discretization error at these resolutions.
• No exact-solution comparison yet. Scattering by a sphere has a classical series solution (this is the 3D analogue of the cylinder's Bessel-Hankel series), and comparing against it would put a number on the total error.

Scenes#

The aperture screen is slower than the background rather than faster. Reflection at an interface goes as |c2 - c1|/(c2 + c1) at constant density, so a screen at c = 0.2 reflects about as much as one at c = 5 would, but the timestep is set by the fastest speed anywhere on the grid: a slow screen is free while a fast one taxes every step of the whole run. What the slow screen costs instead is resolution inside itself, which its own absorption swallows.

Tests#

npm run smoke        # the page itself, in headless Chrome

The browser check drives the real page: it loads it, waits for the solver to take steps and the frame loop to turn, halves the timestep through the slider and requires dt to follow, checks the microphone is recording, and swaps the scene. It says nothing about whether the picture is right, which is what eyes are for.

Development#

npm install
npm run dev

No dependencies beyond the build tooling: the solver and renderer are plain TypeScript and WGSL.

License#

Apache-2.0