acoustic-scattering-2d#

Sound scattering off obstacles in two dimensions, solved live in the browser. The wave equation and the medium are both written as MATLAB, run by numbl and compiled to WebGPU compute shaders; the pressure field is drawn straight out of the buffer the solver writes, with no readback in the display path.

This is the flat sibling of turing-surface, which solves reaction-diffusion on closed surfaces with spherical harmonics. The compiler pipeline is the same idea, cut down: there is one host-provided operation instead of six transforms, and the grid is a square.

Two files#

A scene says what the medium is. It is evaluated once, on the CPU, through numbl's interpreter, so it has the whole language available: loops, if, indexing, seeded randomness, anything in tools/.

function [c, sig] = medium(x, y, L, h, cin, R, absorb)
  r = sqrt(x.^2 + y.^2);
  inside = 0.5 * (1 - tanh((r - R) / (1.5*h)));

  c = 1 + (cin - 1) * inside;
  sig = sponge(x, y, L, 0.2*L, 25) + absorb * inside;
end

c is the sound speed and sig the absorption rate. Both are ordinary fields of position, which is what lets one function describe both the scatterer and the open boundary: the absorbing layer around the edge is the statement that the medium swallows sound out there, and a scene is free to put absorption inside the domain too, making a lossy scatterer.

A model says how the field advances. It is compiled, so it is held to an element-wise subset, plus the Laplacian stencils the host supplies.

function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
  ...
  sd = sig * dt;
  lap = lap2(p);
  pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);

  pold = p;
  tn = t + dt;
end

Neither file declares its own parameters. It names the arguments it wants, and the host matches each against what it offers (src/mgpu/registry.ts and src/scene/registry.ts); a name nobody provides is a compile error with a position in the file, not a silently undefined variable.

The equation#

Pressure only, at constant density:

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

Centring both the second time derivative and the damping on step n gives the explicit update above. 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.

Time is carried as a grid field rather than a number. A batch of timesteps is one replay of a fixed sequence of GPU operations, so nothing the host writes between frames can change inside it; a clock uploaded per frame would stand still for the whole batch. tn = t + dt makes the model keep its own clock, at the cost of one buffer and one kernel per step.

How the MATLAB becomes a shader#

1. Lowering. parseMFile and specializeUserFunction (numbl's own JIT front end) lower init and step for the concrete argument types of the current grid, so every type and shape is fixed at compile time and a step needs no dynamic dispatch at all.
2. Fusion. numbl's inline pass folds single-use temporaries into their consumers. It declines the ones its C back end cannot fuse (sin, exp, tanh), which WGSL is perfectly happy with, so a second pass (src/mgpu/fuse.ts, adapted from math-webgpu-sandbox) folds those too. The whole source term collapses into the update.
3. Planning. Each remaining statement becomes one GPU operation: an element-wise kernel, or a stencil dispatch for lap2/lap4. The result is a static op list, so running a timestep is pure command recording, and a whole frame's worth of steps goes into one submit.

On this machine the compiled step is five kernels, one stencil dispatch, and three buffer copies:

kernel  u = <element-wise>
kernel  sd = <element-wise>
stencil lap = lap2(p)
kernel  pn = <element-wise>
kernel  pold = <element-wise>
kernel  tn = <element-wise>
copy    pn -> p
copy    pold -> pm
copy    tn -> t

The op list is shown next to the editor in the app, so what a line of MATLAB costs is visible while it is being written.

The one thing that is not element-wise#

lap2 and lap4 are the only places where a grid point reads its neighbours. Keeping them as named operations rather than as array slicing means a model never has to know how the grid is laid out in the buffer, and the host is free to implement each as a single dispatch. lap2 is the 5-point stencil; lap4 is the fourth-order 9-point one, which costs twice the reads and buys much less grid dispersion, so a pulse can travel many wavelengths without visibly coming apart.

Kernels that will not fit#

A fused kernel binds one storage buffer per distinct field its line reads, plus its output and the parameter block, and WebGPU guarantees only eight per compute stage (Chrome gives exactly that on some machines; compatibility mode gives four). The leapfrog update reads nine fields. numbl's inline pass does not know about that limit, and a model has no way to ask it for less, since a temporary used once is exactly what it folds away.

So the planner enforces the budget itself: any child subtree that reads more than one field is evaluated into its own buffer and replaced by a reference to it, which leaves the parent reading at most one field per child. The arithmetic is unchanged, in a few more passes over memory, and it only happens on a line that would not otherwise compile. The test suite runs the whole model with the budget squeezed to two fields and checks the answer against the unsplit one.

What it is honest about#

• Single precision. WebGPU has no f64. The stencil differences lose a few digits to cancellation, which is visible in the accuracy numbers below but far below the discretization error at any resolution the app runs interactively.
• The absorbing layer is a sponge, not a PML. An absorbing layer is itself an impedance mismatch, so it reflects; spreading it over a couple of wavelengths keeps that small but not zero. Measured, for the shipped profile: about 0.8% of the incident peak comes back into the interior. A perfectly matched layer would do much better, particularly at grazing incidence, at the cost of extra fields and a more involved update.
• The absorbing layer costs domain. It occupies 20% of the width at each edge, so the clear region is the middle 60%. Scenes and sources are placed with that in mind.
• Grid dispersion is real. A discrete Laplacian gets the wave speed slightly wrong at wavelengths the grid barely resolves, and the error accumulates with distance travelled. Raise the frequency in the app until leapfrog visibly trails, then switch to leapfrog4 and watch it clean up.
• No exact-solution comparison yet. Scattering by a circular cylinder has a classical Bessel-Hankel series solution, and comparing against it would put a number on the total error rather than on the pieces. That is the obvious next thing to do.

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 that runs after each step and appends to a buffer, because the obvious implementation, reading the field back and picking out one number, costs a GPU-to-CPU round trip per step and would be slower than the step. The whole trace comes back once, when there is something to play.

Turning that into sound needs one choice, because the simulation has no seconds in it: its clock is model time, in which the background speed is 1 and the source frequency f is in cycles per model time unit. The choice is made as a pitch. Asking for the source to sound at 440 Hz fixes the playback sample rate at pitch / (f * dt), since the trace holds one sample per timestep.

The consequence is worth knowing before it surprises you: a pulse is short. The default source is a few dozen cycles, and a few dozen cycles at a musical pitch is a few tens of milliseconds however long you leave the simulation running. What you hear is a click, correctly shaped by whatever the wave did on its way to the microphone, but a click. For sustained sound, turn continuous up so the source keeps going, and run long enough to fill some seconds of audio: the speed control goes to 64 timesteps a frame for exactly that. The trace holds 262,144 samples, around twelve seconds at a typical playback 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.

Scenes#

The room is the scene the microphone is for: put a click inside it and what comes back is a direct arrival, then echoes, then reverberation. Its walls are slower than the room rather than faster, which is worth explaining. Reflection at an interface goes as |c2 - c1|/(c2 + c1) at constant density, so a wall at c = 0.15 reflects about 75% of the pressure amplitude, as a wall at c = 7 would. But the timestep is set by the fastest speed anywhere on the grid, so a rigid wall makes every step of the whole simulation several times smaller while a slow one is free. What a slow wall costs instead is resolution inside itself, which is what the wall's absorption is for: it swallows what gets in, so the badly resolved part never comes back out. Push wall speed above 1 for a genuinely rigid room and watch the timestep readout drop to match.

A scene may also ask for source settings — the room asks for a point source inside itself, since a plane wave arriving from off-grid would have nothing to do. Those land in the sliders like any other value, and can be moved afterwards.

The source blends continuously between a line source (whose far field is a plane wave) and a point source, and between a single Gaussian pulse and a wave that turns on and stays on. Both are parameters rather than modes, because both are one expression in the .m.

Tests#

npm run test:node    # the solver, on desktop WebGPU (Google Dawn)
npm run smoke        # the page itself, in headless Chrome

The microphone is checked exactly rather than approximately: after N steps its trace must be N samples long, and its last sample must equal the pressure sitting at the probe's grid point, both being the same f32 written by the same kernel. That pins the probe index, the grid layout, and the fact that the recording dispatch runs once per timestep rather than once per submission.

The solver checks are physical wherever they can be. A pulse must travel at the speed the medium says (measured wavefront radius, within 1%); a disk whose speed matches its background must not scatter at all (bit-identical to the uniform medium); a hard disk must scatter strongly; the absorbing layer must leave under 2% behind; both schemes must be stable at 95% of the CFL limit they claim. The stencils are compared against a field whose Laplacian is known exactly, which at n = 64 gives 1.1e-2 relative error for lap2 and 1.9e-4 for lap4.

The browser check drives the real page: it loads it, waits for the solver to take steps and the frame loop to turn, and swaps the scene. It says nothing about whether the picture is right, which is what eyes are for.

node scripts/smoke.mjs --recompile goes further and swaps the model, breaks the source deliberately (requiring the failure to be reported rather than thrown), and reverts. That is opt-in because it needs new GPU pipelines built while the page is already drawing, which headless Chrome cannot always do: on some machines every GPU object created after the canvas context has been used fails with "A valid external Instance reference no longer exists", in a thirty-line WebGPU page with none of this project in it. A failure there says as much about the browser as about the app.

Development#

npm install        # numbl must be checked out as a sibling: ../../numbl
npm run dev

The compiler runs from numbl's TypeScript sources directly (through a numbl-src vite alias), so no numbl build is needed.

License#

Apache-2.0