1# acoustic-scattering-2d
3Sound scattering off obstacles in two dimensions, solved live in the browser.
4The wave equation and the medium are both written as MATLAB, run by
5[numbl](https://numbl.org) and compiled to WebGPU compute shaders; the pressure
6field is drawn straight out of the buffer the solver writes, with no readback in
7the display path.
9This is the flat sibling of
10[turing-surface](https://github.com/concept-collection/turing-surface), which
11solves reaction-diffusion on closed surfaces with spherical harmonics. The
12compiler pipeline is the same idea, cut down: there is one host-provided
13operation instead of six transforms, and the grid is a square.
15## Two files
17A **scene** says what the medium is. It is evaluated once, on the CPU, through
18numbl's interpreter, so it has the whole language available: loops, `if`,
19indexing, seeded randomness, anything in `tools/`.
21```matlab
22function [c, sig] = medium(x, y, L, h, cin, R, absorb)
23 r = sqrt(x.^2 + y.^2);
24 inside = 0.5 * (1 - tanh((r - R) / (1.5*h)));
26 c = 1 + (cin - 1) * inside;
27 sig = sponge(x, y, L, 0.2*L, 25) + absorb * inside;
28end
29```
31`c` is the sound speed and `sig` the absorption rate. Both are ordinary fields
32of position, which is what lets one function describe both the scatterer and the
33open boundary: the absorbing layer around the edge is the statement that the
34medium swallows sound out there, and a scene is free to put absorption inside
35the domain too, making a lossy scatterer.
37A **model** says how the field advances. It is compiled, so it is held to an
38element-wise subset, plus the Laplacian stencils the host supplies.
40```matlab
41function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
42 ...
43 sd = sig * dt;
44 lap = lap2(p);
45 pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);
47 pold = p;
48 tn = t + dt;
49end
50```
52Neither file declares its own parameters. It names the arguments it wants, and
53the host matches each against what it offers (`src/mgpu/registry.ts` and
54`src/scene/registry.ts`); a name nobody provides is a compile error with a
55position in the file, not a silently undefined variable.
57## The equation
59Pressure only, at constant density:
61```
62p_tt + 2*sig*p_t = c(x,y)^2 * lap(p) + s(x, y, t)
63```
65Centring both the second time derivative and the damping on step *n* gives the
66explicit update above. Since the density is constant, the impedance ratio across
67an interface is just the speed ratio, so a scatterer much faster than its
68background behaves nearly rigid (sound-hard) and one much slower nearly
69pressure-release (sound-soft). What this formulation cannot do is set impedance
70and speed independently, which needs a variable-density divergence form and a
71second field.
73Time is carried as a grid field rather than a number. A batch of timesteps is
74one replay of a fixed sequence of GPU operations, so nothing the host writes
75between frames can change inside it; a clock uploaded per frame would stand
76still for the whole batch. `tn = t + dt` makes the model keep its own clock, at
77the cost of one buffer and one kernel per step.
79## How the MATLAB becomes a shader
811. **Lowering.** `parseMFile` and `specializeUserFunction` (numbl's own JIT
82 front end) lower `init` and `step` for the concrete argument types of the
83 current grid, so every type and shape is fixed at compile time and a step
84 needs no dynamic dispatch at all.
852. **Fusion.** numbl's inline pass folds single-use temporaries into their
86 consumers. It declines the ones its C back end cannot fuse (`sin`, `exp`,
87 `tanh`), which WGSL is perfectly happy with, so a second pass
88 (`src/mgpu/fuse.ts`, adapted from
89 [math-webgpu-sandbox](https://github.com/concept-collection/math-webgpu-sandbox))
90 folds those too. The whole source term collapses into the update.
913. **Planning.** Each remaining statement becomes one GPU operation: an
92 element-wise kernel, or a stencil dispatch for `lap2`/`lap4`. The result is a
93 static op list, so running a timestep is pure command recording, and a whole
94 frame's worth of steps goes into one submit.
96On this machine the compiled step is five kernels, one stencil dispatch, and
97three buffer copies:
99```
100kernel u = <element-wise>
101kernel sd = <element-wise>
102stencil lap = lap2(p)
103kernel pn = <element-wise>
104kernel pold = <element-wise>
105kernel tn = <element-wise>
106copy pn -> p
107copy pold -> pm
108copy tn -> t
109```
111The op list is shown next to the editor in the app, so what a line of MATLAB
112costs is visible while it is being written.
114### The one thing that is not element-wise
116`lap2` and `lap4` are the only places where a grid point reads its neighbours.
117Keeping them as named operations rather than as array slicing means a model
118never has to know how the grid is laid out in the buffer, and the host is free
119to implement each as a single dispatch. `lap2` is the 5-point stencil; `lap4` is
120the fourth-order 9-point one, which costs twice the reads and buys much less
121grid dispersion, so a pulse can travel many wavelengths without visibly coming
122apart.
124### Kernels that will not fit
126A fused kernel binds one storage buffer per distinct field its line reads, plus
127its output and the parameter block, and WebGPU guarantees only eight per compute
128stage (Chrome gives exactly that on some machines; compatibility mode gives
129four). The leapfrog update reads nine fields. numbl's inline pass does not know
130about that limit, and a model has no way to ask it for less, since a temporary
131used once is exactly what it folds away.
133So the planner enforces the budget itself: any child subtree that reads more
134than one field is evaluated into its own buffer and replaced by a reference to
135it, which leaves the parent reading at most one field per child. The arithmetic
136is unchanged, in a few more passes over memory, and it only happens on a line
137that would not otherwise compile. The test suite runs the whole model with the
138budget squeezed to two fields and checks the answer against the unsplit one.
140## What it is honest about
142- **Single precision.** WebGPU has no f64. The stencil differences lose a few
143 digits to cancellation, which is visible in the accuracy numbers below but far
144 below the discretization error at any resolution the app runs interactively.
145- **The absorbing layer is a sponge, not a PML.** An absorbing layer is itself an
146 impedance mismatch, so it reflects; spreading it over a couple of wavelengths
147 keeps that small but not zero. Measured, for the shipped profile: about 0.8% of
148 the incident peak comes back into the interior. A perfectly matched layer would
149 do much better, particularly at grazing incidence, at the cost of extra fields
150 and a more involved update.
151- **The absorbing layer costs domain.** It occupies 20% of the width at each
152 edge, so the clear region is the middle 60%. Scenes and sources are placed
153 with that in mind.
154- **Grid dispersion is real.** A discrete Laplacian gets the wave speed slightly
155 wrong at wavelengths the grid barely resolves, and the error accumulates with
156 distance travelled. Raise the frequency in the app until `leapfrog` visibly
157 trails, then switch to `leapfrog4` and watch it clean up.
158- **No exact-solution comparison yet.** Scattering by a circular cylinder has a
159 classical Bessel-Hankel series solution, and comparing against it would put a
160 number on the total error rather than on the pieces. That is the obvious next
161 thing to do.
163## Listening to it
165There is a microphone: the pressure at one grid point, sampled every timestep.
166It is written on the GPU by a one-thread dispatch that runs after each step and
167appends to a buffer, because the obvious implementation, reading the field back
168and picking out one number, costs a GPU-to-CPU round trip per step and would be
169slower than the step. The whole trace comes back once, when there is something
170to play.
172Turning that into sound needs one choice, because the simulation has no seconds
173in it: its clock is model time, in which the background speed is 1 and the
174source frequency `f` is in cycles per model time unit. The choice is made as a
175pitch. Asking for the source to sound at 440 Hz fixes the playback sample rate
176at `pitch / (f * dt)`, since the trace holds one sample per timestep.
178The consequence is worth knowing before it surprises you: **a pulse is short.**
179The default source is a few dozen cycles, and a few dozen cycles at a musical
180pitch is a few tens of milliseconds however long you leave the simulation
181running. What you hear is a click, correctly shaped by whatever the wave did on
182its way to the microphone, but a click. For sustained sound, turn `continuous`
183up so the source keeps going, and run long enough to fill some seconds of audio:
184the speed control goes to 64 timesteps a frame for exactly that. The trace holds
185262,144 samples, around twelve seconds at a typical playback rate.
187The recording is normalized before playback, which discards absolute amplitude:
188that is what the colour scale is for. It restarts whenever the run does, and
189whenever the timestep changes, since a trace is one sample per step and two
190timesteps would be two sample rates in one buffer.
192## Scenes
194| scene | what it shows |
195| --- | --- |
196| Disk | One circular scatterer, the reference case. Slow, fast, or absorbing. |
197| Room | A square enclosure with a doorway. Close the doorway and it is sealed. |
198| Two slits | A hard screen with two apertures: diffraction and interference. |
199| Lens | A smooth slow patch that refracts a plane wave to a focus. |
200| Random medium | Weak random structure everywhere: multiple scattering, and a coda. |
202The room is the scene the microphone is for: put a click inside it and what
203comes back is a direct arrival, then echoes, then reverberation. Its walls are
204*slower* than the room rather than faster, which is worth explaining. Reflection
205at an interface goes as `|c2 - c1|/(c2 + c1)` at constant density, so a wall at
206c = 0.15 reflects about 75% of the pressure amplitude, as a wall at c = 7 would.
207But the timestep is set by the fastest speed anywhere on the grid, so a rigid
208wall makes every step of the whole simulation several times smaller while a slow
209one is free. What a slow wall costs instead is resolution inside itself, which
210is what the wall's absorption is for: it swallows what gets in, so the badly
211resolved part never comes back out. Push `wall speed` above 1 for a genuinely
212rigid room and watch the timestep readout drop to match.
214A scene may also ask for source settings — the room asks for a point source
215inside itself, since a plane wave arriving from off-grid would have nothing to
216do. Those land in the sliders like any other value, and can be moved afterwards.
218The source blends continuously between a line source (whose far field is a plane
219wave) and a point source, and between a single Gaussian pulse and a wave that
220turns on and stays on. Both are parameters rather than modes, because both are
221one expression in the .m.
223## Tests
225```
226npm run test:node # the solver, on desktop WebGPU (Google Dawn)
227npm run smoke # the page itself, in headless Chrome
228```
230The microphone is checked exactly rather than approximately: after N steps its
231trace must be N samples long, and its last sample must equal the pressure
232sitting at the probe's grid point, both being the same f32 written by the same
233kernel. That pins the probe index, the grid layout, and the fact that the
234recording dispatch runs once per timestep rather than once per submission.
236The solver checks are physical wherever they can be. A pulse must travel at the
237speed the medium says (measured wavefront radius, within 1%); a disk whose speed
238matches its background must not scatter at all (bit-identical to the uniform
239medium); a hard disk must scatter strongly; the absorbing layer must leave under
2402% behind; both schemes must be stable at 95% of the CFL limit they claim. The
241stencils are compared against a field whose Laplacian is known exactly, which at
242n = 64 gives 1.1e-2 relative error for `lap2` and 1.9e-4 for `lap4`.
244The browser check drives the real page: it loads it, waits for the solver to
245take steps and the frame loop to turn, and swaps the scene. It says nothing
246about whether the picture is right, which is what eyes are for.
248`node scripts/smoke.mjs --recompile` goes further and swaps the model, breaks
249the source deliberately (requiring the failure to be reported rather than
250thrown), and reverts. That is opt-in because it needs new GPU pipelines built
251while the page is already drawing, which headless Chrome cannot always do: on
252some machines every GPU object created after the canvas context has been used
253fails with "A valid external Instance reference no longer exists", in a
254thirty-line WebGPU page with none of this project in it. A failure there says
255as much about the browser as about the app.
257## Development
259```
260npm install # numbl must be checked out as a sibling: ../../numbl
261npm run dev
262```
264The compiler runs from numbl's TypeScript sources directly (through a
265`numbl-src` vite alias), so no numbl build is needed.
267## License
269Apache-2.0