1% Driver for Chapter 7: letting sound out of a grid.
2%
3% The editable script supplies
4%
5% U = openbox(u0, uold, lam, M, method, W, keep)
6%
7% which runs leapfrog on a square region with some treatment of its edges,
8% and returns the region's field at the nodes `keep` after every step.
9%
10% The reference is exact without being a formula. An explicit scheme moves
11% information by at most one node per step, so the same leapfrog scheme on a
12% grid that extends M nodes further in every direction cannot feel its own
13% edges anywhere in the region during M steps: inside the region, it is the
14% scheme's free-space solution, to rounding error. The driver runs that
15% reference with its own code. Whatever differs from it in the region is
16% what the edge treatment sent back, and nothing else.
17%
18% mode 'snap' one method: snapshots, and the reflection it produced.
19% mode 'sweep' every method at several layer widths: reflection against
20% cost.
21% mode 'hear' a longer run with a click, heard at one point, for a rigid
22% box and for two absorbing edges.
23p = wv_params();
24out = struct();
25N = p.N;
26n0 = N + 1;
27lam = p.lam;
28h = p.L / N;
29c = 343;
30dt = lam * h / c;
31x = (0:N)' * h;
32[X, Y] = meshgrid(x, x);
33sg = p.sigma;
34u0 = exp(-((X - p.xs) .^ 2 + (Y - p.ys) .^ 2) / (2 * sg ^ 2));
35% released from rest: uold = u0 + (lam^2/2) Lap_h u0, the same for every run
36lap = zeros(n0, n0);
37lap(2:N, 2:N) = u0(3:n0, 2:N) + u0(1:N - 1, 2:N) + u0(2:N, 3:n0) + u0(2:N, 1:N - 1) - 4 * u0(2:N, 2:N);
38um = u0 + 0.5 * lam ^ 2 * lap;
40if strcmp(p.mode, 'hear')
41 M = round(p.T / dt);
42 ir = round(p.xr / h) + 1;
43 jr = round(p.yr / h) + 1;
44 meths = p.methods;
45 Yh = zeros(M, numel(meths));
46 for k = 1:numel(meths)
47 U = openbox(u0, um, lam, M, meths{k}, p.W, (jr - 1) * n0 + ir);
48 Yh(:, k) = U(1, :)';
49 end
50 s = max(abs(Yh(:)));
51 fs = 1 / dt;
52 clips = cell(1, numel(meths));
53 for k = 1:numel(meths)
54 clips{k} = wv_audio(Yh(:, k), fs, s);
55 end
56 out.clips = clips;
57 out.methods = meths;
58 out.fs = fs;
59 out.t = (1:M) * dt * 1000;
60 nkeep = min(M, round(0.08 / dt));
61 out.tt = out.t(1:nkeep);
62 out.traces = wv_rows(Yh(1:nkeep, :)' / s);
63else
64 M = round(p.T / dt);
65 Uref = drv_reference(u0, um, lam, M);
66 rn = sum(Uref(:) .^ 2);
67 if strcmp(p.mode, 'snap')
68 t0 = tic;
69 U = openbox(u0, um, lam, M, p.method, p.W, (1:n0 * n0)');
70 out.ms = toc(t0) * 1000;
71 E = U - Uref;
72 out.refl = 10 * log10(sum(E(:) .^ 2) / rn);
73 ks = min(M, round(p.tsnap / dt));
74 out.tsnap = ks * dt * 1000;
75 out.ref = wv_rows(reshape(Uref(:, ks), n0, n0));
76 out.num = wv_rows(reshape(U(:, ks), n0, n0));
77 out.err = wv_rows(reshape(E(:, ks), n0, n0));
78 out.emax = max(abs(E(:)));
79 out.umax = max(abs(Uref(:, ks)));
80 out.urefmax = max(abs(Uref(:)));
81 % the reflected energy in the region, step by step
82 out.t = (1:M) * dt * 1000;
83 out.eng = 10 * log10(max(sum(E .^ 2, 1) / max(sum(Uref .^ 2, 1)), 1e-30));
84 out.L = p.L;
85 else
86 meths = p.methods;
87 Ws = p.Ws(:)';
88 R = NaN(numel(meths), numel(Ws));
89 TM = NaN(numel(meths), numel(Ws));
90 for k = 1:numel(meths)
91 for j = 1:numel(Ws)
92 if (strcmp(meths{k}, 'rigid') || strcmp(meths{k}, 'mur')) && j > 1
93 continue;
94 end
95 t0 = tic;
96 U = openbox(u0, um, lam, M, meths{k}, Ws(j), (1:n0 * n0)');
97 TM(k, j) = toc(t0) * 1000;
98 E = U - Uref;
99 R(k, j) = 10 * log10(sum(E(:) .^ 2) / rn);
100 end
101 end
102 out.methods = meths;
103 out.Ws = Ws;
104 out.refl = wv_rows(R);
105 out.ms = wv_rows(TM);
106 end
107end
109wv_write(out);
111function U = drv_reference(u0, uold, lam, M)
112% The same scheme on a grid M + 2 nodes larger on every side, with its far
113% edges held at zero: nothing from those edges reaches the region in M steps.
114n0 = size(u0, 1);
115e = M + 2;
116n1 = n0 + 2 * e;
117u = zeros(n1, n1);
118uo = zeros(n1, n1);
119u(e + 1:e + n0, e + 1:e + n0) = u0;
120uo(e + 1:e + n0, e + 1:e + n0) = uold;
121un = zeros(n1, n1);
122l2 = lam ^ 2;
123U = zeros(n0 * n0, M);
124for n = 1:M
125 % only the nodes that can have become nonzero need updating
126 lo = max(2, e + 1 - n);
127 hi = min(n1 - 1, e + n0 + n);
128 for j = lo:hi
129 for i = lo:hi
130 un(i, j) = 2 * u(i, j) - uo(i, j) + l2 * (u(i + 1, j) + u(i - 1, j) + u(i, j + 1) + u(i, j - 1) - 4 * u(i, j));
131 end
132 end
133 uo = u;
134 u = un;
135 U(:, n) = reshape(u(e + 1:e + n0, e + 1:e + n0), n0 * n0, 1);
136end
137end