cbf6d88Initial commit: 2D acoustic scattering, live on WebGPUJeremy Magland 1% Second-order leapfrog for the 2D acoustic wave equation.
2%
3% p_tt + 2*sig*p_t = c^2 * lap(p) + s(x, y, t)
4%
5% Pressure only, at constant density, so the medium is the sound speed c
6% (metres per second) and the absorption sig (inverse seconds) that the scene
7% defines. Everything here is SI: x, y in metres, t and dt in seconds, f in
8% hertz — none of it is declared as such anywhere, it simply follows from x,
9% y and t being metres and seconds, which is the scene's and the app's doing,
10% not this file's. Centring both the second time
11% derivative and the damping term on step n,
12%
13% (p+ - 2p + p-)/dt^2 + sig*(p+ - p-)/dt = c^2*lap(p) + s
14%
15% and solving for p+ gives the update below. It is explicit: the only thing
16% that couples neighbouring points is `lap2`, the 5-point Laplacian, which the
17% host supplies as a single GPU dispatch. Everything else here is element-wise
18% and compiles to one kernel per line.
19%
20% Stability wants c*dt/h <= 1/sqrt(2); the app sets dt from the fastest speed
21% anywhere in the scene, so a fast scatterer slows the whole run down.
22%
23% Time is carried as a field rather than a number, because a batch of steps is
24% one replay of a fixed sequence of GPU operations and nothing the host writes
25% between frames can change inside it. `tn = t + dt` makes the model keep its
26% own clock, which is what lets the source term be correct however many steps
27% are batched into a submit.
29function [p, pm, t] = init(npts)
30 p = zeros(npts, 1);
31 pm = zeros(npts, 1);
32 t = zeros(npts, 1);
33end
35function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
36 % The source. `cw` blends between a Gaussian pulse (0) and a wave that
37 % turns on smoothly and stays on (1); `point` blends between a line source
38 % spanning the grid in y, whose far field is a plane wave, and a point
39 % source at (x0, y0).
40 %
41 % The om^2 is only a choice of units, not a physical amplitude — this is a
42 % body force of arbitrary strength, not a source with a real acoustic power
43 % rating. Such a force drives a response that falls off as 1/om^2 — two
44 % time integrations — so without the factor the field would shrink tenfold
45 % every time the frequency slider tripled. The equation is linear; scaling
46 % the source scales the answer and nothing else.
47 u = (t - t0) / tw;
48 env = (1 - cw) * exp(-u .* u) + cw * (0.5 * (1 + tanh(u)));
49 gx = ((x - x0) / w) .^ 2;
50 gy = point * (((y - y0) / w) .^ 2);
51 om = 2*pi*f;
52 s = (om*om) * (env .* sin(om*(t - t0)) .* exp(-(gx + gy)));
54 % One step. The damping is what the absorbing layer acts through: sig is
55 % zero over the interior, so there p+ is the plain leapfrog update.
56 sd = sig * dt;
57 lap = lap2(p);
58 pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);
60 % This step's field becomes the next step's history. A line on its own, so
61 % it plans as a copy rather than being folded into the update above.
62 pold = p;
63 tn = t + dt;
64end