1function R = plate_fd_loop(S)
2% PLATE_FD_LOOP The explicit time stepping of plate_fd.
3%
4% R = plate_fd_loop(S), with S from plate_fd_setup, advances
5%
6% wn = c_new * (2 w - c_old wm - k2dt2 B + s1dt2 (Lw - Lw_m) + fdt2 F(n) s)
7%
8% from rest for S.nsteps steps, where Lw is the five-point Laplacian of the
9% current field, B the biharmonic (the Laplacian applied twice, the 13-point
10% stencil) and Lw_m the Laplacian of the previous field. On a rectangle with
11% supported or clamped edges the operators are the dense stencils with one
12% layer of ghost nodes (plate_pad) and the edge nodes are held at zero.
13% Otherwise (S.sparse) they are the sparse operators S.Lm, S.Km of
14% plate_matrices on every node, divided by the node weights S.mvec; the
15% nodes that are not unknowns (off the plate, on a supported ring, the pin)
16% are zeroed with S.keep, and the S.nr rigid-body motions of a free plate
17% are projected out of the new field with S.Rb, S.Qb.
18%
19% What the loop records, and returns in R:
20% w_p, wxx_p, wyy_p, wxy_p displacement and curvatures at the probe node,
21% 1 x (nsteps + 1); entry n + 1 is step n (t = n dt)
22% Wf the field at the frame steps, ny x nx x nf
23% Pbins radiated pressure per listener delay bin and
24% audio sample, nbins x na: the acceleration of
25% every node at every step, weighted by
26% S.wgt and box-averaged onto audio samples
27% diverged_step 0, or the step at which max |w| reached
28% S.blowup (the run stops there)
29% notes remarks about the run
30% backend 'numbl': where the loop ran
31% loop_s wall time of the loop by tic/toc
32%
33% The lab has a WebGPU port of this function (packages/core/src/gpu/fdLoop.ts)
34% that computes the same record from the same S, in single precision. It
35% carries the state as (w, w - wm) rather than (w, wm), which is the same
36% arithmetic rearranged so that the rounding of w does not feed back into the
37% velocity; keep the two implementations in step.
39t_loop = tic;
40nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; nsteps = S.nsteps; fs = S.fs;
41bc = S.P.bc; use_sparse = S.sparse; nr = S.nr;
42s = S.s; F = S.F; ip = S.ip; jp = S.jp;
43frame_steps = S.frame_steps; nf = S.nf;
44bins = S.bins; wgt = S.wgt; nbins = S.Lis.nbins; na = S.na;
45c_old = S.c_old; c_new = S.c_new; k2dt2 = S.k2dt2; s1dt2 = S.s1dt2; fdt2 = S.fdt2; blowup = S.blowup;
46if use_sparse
47 Lm = S.Lm; Km = S.Km; mvec = S.mvec; keep = S.keep; Rb = S.Rb; Qb = S.Qb;
48end
50% State
51w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
52w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
53Wf = zeros(ny, nx, nf);
54Pbins = zeros(nbins, na);
55diverged_step = 0;
56notes = {};
57kf = 1;
59for nstep = 1:nsteps
60 if use_sparse
61 Lw = reshape((Lm * w(:)) ./ mvec, ny, nx);
62 B = reshape((Km * w(:)) ./ mvec, ny, nx);
63 else
64 Lw = plate_laplacian(plate_pad(w, bc), dx, dy);
65 B = zeros(ny, nx);
66 B(2:end - 1, 2:end - 1) = (Lw(3:end, 2:end - 1) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(1:end - 2, 2:end - 1)) / dy^2 ...
67 + (Lw(2:end - 1, 3:end) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(2:end - 1, 1:end - 2)) / dx^2;
68 end
69 wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
70 if use_sparse
71 wn = wn .* keep;
72 if nr > 0
73 wn = wn - reshape(Rb * (Qb' * wn(:)), ny, nx);
74 end
75 else
76 wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;
77 end
79 % Probe displacement and curvatures from the 3 x 3 neighbourhood
80 w_p(nstep + 1) = wn(ip, jp);
81 wxx_p(nstep + 1) = (wn(ip, jp + 1) - 2 * wn(ip, jp) + wn(ip, jp - 1)) / dx^2;
82 wyy_p(nstep + 1) = (wn(ip + 1, jp) - 2 * wn(ip, jp) + wn(ip - 1, jp)) / dy^2;
83 wxy_p(nstep + 1) = (wn(ip + 1, jp + 1) - wn(ip + 1, jp - 1) - wn(ip - 1, jp + 1) + wn(ip - 1, jp - 1)) / (4 * dx * dy);
85 % Radiated pressure: acceleration at this step, weighted and delayed per node
86 acc = (wn - 2 * w + wm) / dt^2;
87 ia = floor(nstep * dt * fs) + 1;
88 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [nbins 1]) * (dt * fs);
90 if kf <= nf && nstep == frame_steps(kf)
91 Wf(:, :, kf) = wn;
92 kf = kf + 1;
93 end
95 peak = max(abs(wn(:)));
96 if ~(peak < blowup) % also catches NaN and Inf
97 diverged_step = nstep;
98 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
99 break
100 end
101 wm = w; w = wn; Lw_m = Lw;
102end
104R.w_p = w_p; R.wxx_p = wxx_p; R.wyy_p = wyy_p; R.wxy_p = wxy_p;
105R.Wf = Wf;
106R.Pbins = Pbins;
107R.diverged_step = diverged_step;
108R.notes = notes;
109R.backend = 'numbl';
110R.loop_s = toc(t_loop);
111end