d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% Four ways to end a grid, for waves that ought to leave it.
2%
3% The region of interest is the square of (N+1)-by-(N+1) nodes that u0
4% covers. Around it, a method may add a layer of W extra nodes on every
5% side, which costs work but is never looked at. The interior is Chapter 4's
6% leapfrog scheme; the methods differ only at the edge:
7%
8% 'rigid' a hard wall, du/dn = 0, by mirroring the first interior node
9% into a ghost. Everything comes back, with the same sign.
10% 'mur' Mur's (1981) first-order absorbing condition, a discrete
11% u_t + c u_n = 0 on each edge. It is exact for a wave
12% arriving square-on and reflects a fraction
13% (1 - cos a) / (1 + cos a) of one arriving at angle a.
14% 'sponge' a damping layer: the equation u_tt + 2 s(x) u_t = c^2 Lap u,
15% with s rising from zero to smax across W nodes as the square
16% of the depth, and a zero-displacement edge behind it. A wave
17% is damped on the way in and again on the way back out, and
18% part of it is reflected by the change of s itself.
19% 'pml' a perfectly matched layer (Berenger 1994). The scheme is
20% rewritten in pressure and velocity on a staggered grid,
21% which in the region is the same leapfrog scheme exactly;
22% in the layer the pressure is split into the parts driven by
23% the x and the y velocity, and each is damped only in its own
24% direction. In the continuum this absorbs a wave at any
25% angle and any frequency without reflecting any of it; on
26% the grid, the profile's steps reflect a little.
27%
28% u0 and uold are the region's field at steps 0 and -1. The result U holds
29% the region's field at the nodes `keep` (linear indices into the region)
30% after each of the M steps, one column per step.
32function U = openbox(u0, uold, lam, M, method, W, keep)
33n0 = size(u0, 1);
34if strcmp(method, 'rigid') || strcmp(method, 'mur')
35 W = 0;
36end
37if strcmp(method, 'pml')
38 U = drv_pml(u0, uold, lam, M, W, keep);
39 return;
40end
41n1 = n0 + 2 * W;
42u = zeros(n1, n1);
43uo = zeros(n1, n1);
44u(W + 1:W + n0, W + 1:W + n0) = u0;
45uo(W + 1:W + n0, W + 1:W + n0) = uold;
46un = zeros(n1, n1);
47l2 = lam ^ 2;
48% the damping, per step: s dt at each node, zero in the region itself
49sd = zeros(n1, n1);
50if W > 0
51 % s dt at the outer edge. A stronger sponge damps more but reflects more
52 % from its own gradient; this strength did best of the few tried.
53 smax = 3 * lam / (2 * W) * log(100);
54 for j = 1:n1
55 for i = 1:n1
56 di = max([W + 1 - i, i - (W + n0), 0]);
57 dj = max([W + 1 - j, j - (W + n0), 0]);
58 d = max(di, dj) / W;
59 sd(i, j) = smax * d ^ 2;
60 end
61 end
62end
63rigid = strcmp(method, 'rigid');
64mur = strcmp(method, 'mur');
65km = (lam - 1) / (lam + 1);
66U = zeros(numel(keep), M);
67for n = 1:M
68 for j = 2:n1 - 1
69 for i = 2:n1 - 1
70 a = sd(i, j);
71 un(i, j) = (2 * u(i, j) - (1 - a) * uo(i, j) ...
72 + l2 * (u(i + 1, j) + u(i - 1, j) + u(i, j + 1) + u(i, j - 1) - 4 * u(i, j))) / (1 + a);
73 end
74 end
75 if rigid
76 % the edge node sees its inner neighbour twice
77 for k = 2:n1 - 1
78 un(1, k) = 2 * u(1, k) - uo(1, k) + l2 * (2 * u(2, k) + u(1, k + 1) + u(1, k - 1) - 4 * u(1, k));
79 un(n1, k) = 2 * u(n1, k) - uo(n1, k) + l2 * (2 * u(n1 - 1, k) + u(n1, k + 1) + u(n1, k - 1) - 4 * u(n1, k));
80 un(k, 1) = 2 * u(k, 1) - uo(k, 1) + l2 * (2 * u(k, 2) + u(k + 1, 1) + u(k - 1, 1) - 4 * u(k, 1));
81 un(k, n1) = 2 * u(k, n1) - uo(k, n1) + l2 * (2 * u(k, n1 - 1) + u(k + 1, n1) + u(k - 1, n1) - 4 * u(k, n1));
82 end
83 un(1, 1) = 2 * u(1, 1) - uo(1, 1) + l2 * (2 * u(2, 1) + 2 * u(1, 2) - 4 * u(1, 1));
84 un(1, n1) = 2 * u(1, n1) - uo(1, n1) + l2 * (2 * u(2, n1) + 2 * u(1, n1 - 1) - 4 * u(1, n1));
85 un(n1, 1) = 2 * u(n1, 1) - uo(n1, 1) + l2 * (2 * u(n1 - 1, 1) + 2 * u(n1, 2) - 4 * u(n1, 1));
86 un(n1, n1) = 2 * u(n1, n1) - uo(n1, n1) + l2 * (2 * u(n1 - 1, n1) + 2 * u(n1, n1 - 1) - 4 * u(n1, n1));
87 elseif mur
88 for k = 2:n1 - 1
89 un(1, k) = u(2, k) + km * (un(2, k) - u(1, k));
90 un(n1, k) = u(n1 - 1, k) + km * (un(n1 - 1, k) - u(n1, k));
91 un(k, 1) = u(k, 2) + km * (un(k, 2) - u(k, 1));
92 un(k, n1) = u(k, n1 - 1) + km * (un(k, n1 - 1) - u(k, n1));
93 end
94 un(1, 1) = 0.5 * (un(2, 1) + un(1, 2));
95 un(1, n1) = 0.5 * (un(2, n1) + un(1, n1 - 1));
96 un(n1, 1) = 0.5 * (un(n1 - 1, 1) + un(n1, 2));
97 un(n1, n1) = 0.5 * (un(n1 - 1, n1) + un(n1, n1 - 1));
98 end
99 uo = u;
100 u = un;
101 R = u(W + 1:W + n0, W + 1:W + n0);
102 U(:, n) = R(keep);
103end
104end
106function U = drv_pml(u0, uold, lam, M, W, keep)
107% Leapfrog in pressure p and scaled velocity w = rho c v, staggered: w_x sits
108% between p(i, j) and p(i, j+1), w_y between p(i, j) and p(i+1, j).
109%
110% w <- w - lam grad p, p <- p - lam div w
111%
112% Eliminating w gives Chapter 4's scheme for p, so in the region this is the
113% same computation. In the layer, p = px + py and every field is damped by
114% the profile of its own direction, s(d) = smax (d/W)^2 at depth d, with
115% smax chosen for a nominal reflection of 1e-5 at normal incidence.
116n0 = size(u0, 1);
117n1 = n0 + 2 * W;
118smax = 3 * lam / (2 * W) * log(1e5);
119% damping per step at whole and half nodes, along one axis
120dw = zeros(1, n1);
121dh = zeros(1, n1 - 1);
122for k = 1:n1
123 d = max([W + 1 - k, k - (W + n0), 0]) / W;
124 dw(k) = smax * d ^ 2;
125end
126for k = 1:n1 - 1
127 kk = k + 0.5;
128 d = max([W + 1 - kk, kk - (W + n0), 0]) / W;
129 dh(k) = smax * d ^ 2;
130end
131p0 = zeros(n1, n1);
132p0(W + 1:W + n0, W + 1:W + n0) = u0;
133pm = zeros(n1, n1);
134pm(W + 1:W + n0, W + 1:W + n0) = uold;
135% start at rest: w at step -1/2 is lam/2 times the gradient of p at step 0,
136% which makes p at step 1 equal to what leapfrog would give from p0 and pm
137wx = zeros(n1, n1 - 1);
138wy = zeros(n1 - 1, n1);
139for j = 1:n1 - 1
140 for i = 1:n1
141 wx(i, j) = 0.5 * lam * (p0(i, j + 1) - p0(i, j));
142 end
143end
144for j = 1:n1
145 for i = 1:n1 - 1
146 wy(i, j) = 0.5 * lam * (p0(i + 1, j) - p0(i, j));
147 end
148end
149% the split: all of p0 starts in px, which is harmless because the split
150% only matters inside the layer, where p0 is zero
151px = p0;
152py = zeros(n1, n1);
153pp = p0;
154% the damped update x <- (1 - a) x / (1 + a) - lam (...) / (1 + a), with its
155% two coefficients worked out once per row or column
156fw = (1 - dw / 2) ./ (1 + dw / 2);
157gw = lam ./ (1 + dw / 2);
158fh = (1 - dh / 2) ./ (1 + dh / 2);
159gh = lam ./ (1 + dh / 2);
160U = zeros(numel(keep), M);
161for n = 1:M
162 for j = 1:n1 - 1
163 f = fh(j);
164 g = gh(j);
165 for i = 1:n1
166 wx(i, j) = f * wx(i, j) - g * (pp(i, j + 1) - pp(i, j));
167 end
168 end
169 for j = 1:n1
170 for i = 1:n1 - 1
171 wy(i, j) = fh(i) * wy(i, j) - gh(i) * (pp(i + 1, j) - pp(i, j));
172 end
173 end
174 for j = 2:n1 - 1
175 f = fw(j);
176 g = gw(j);
177 for i = 2:n1 - 1
178 px(i, j) = f * px(i, j) - g * (wx(i, j) - wx(i, j - 1));
179 py(i, j) = fw(i) * py(i, j) - gw(i) * (wy(i, j) - wy(i - 1, j));
180 pp(i, j) = px(i, j) + py(i, j);
181 end
182 end
183 R = pp(W + 1:W + n0, W + 1:W + n0);
184 U(:, n) = R(keep);
185end
186end