1function Rec = plate_fd_loop(S)
2% PLATE_FD_LOOP The explicit time stepping of plate_fd.
3%
4% Rec = plate_fd_loop(S), with S from plate_fd_setup, advances the departure
5% from the static equilibrium, w = R c + phi, for S.nsteps steps: the rigid-
6% body coordinates c as nr scalar equations and the flexural field phi on
7% the grid,
8%
9% cacc = Ginv (R' f - Krr c - AR' phi)
10% cn = c_new (2 c - c_old cm + dt^2 cacc)
11% acc = (f - KB phi - KG phi - AR c) / M - R cacc
12% phin = c_new (2 phi - c_old phim + dt^2 acc)
13%
14% where f is the nodal force at this step: the strike F(n) sa, the contact
15% and dashpot forces about rest, and the sigma1 loss 2 sigma1 M Lm (w - wm) / dt
16% on the backward-difference velocity of the total field w = R c + phi.
17% c_old = 1 - sigma0 dt and c_new = 1 / (1 + sigma0 dt) carry the
18% frequency-independent loss. Every force enters both parts through f, so
19% the sum R cn + phin is exactly the leapfrog step of the whole system and
20% the losses stay dissipative; the term - R cacc keeps phi free of
21% rigid-body content, so the two parts stay what they are called. (A
22% sigma1 loss applied to phi alone is not symmetric on the whole system and
23% was found to pump the swing of a hanging plate slowly.)
24% The scheme is second order in space and time (Bilbao, Numerical Sound
25% Synthesis, 2009, ch. 12).
26%
27% Contacts. A bottom contact with compression gb at rest pushes up by
28% kb (gb - w) while gb - w > 0, a top contact with compression gt pushes
29% down by kt (gt + w) while gt + w > 0; the static parts are subtracted so
30% f is zero at rest, and the dashpots act on the backward-difference
31% velocity while the contact is engaged.
32%
33% What the loop records, and returns in Rec. The series are averaged onto
34% the output samples (S.na of them at S.fs): each step adds its value,
35% weighted by S.wsamp, to the sample that contains its time, so a sample is
36% the mean over its steps.
37% w_p, wxx_p, wyy_p, wxy_p displacement (total departure) and curvatures
38% (flexural) at the probe node, 1 x na
39% C rigid-body coordinates, nr x na
40% Fsup force of each support, nsup x na [N]: the
41% sum over its nodes of the force of its kind
42% (spring, bottom contact or top contact),
43% static part included
44% Pbins radiated pressure per listener delay bin and
45% sample, nbins x na
46% Wf, Cf phi at the frame steps, ny x nx x nf, and c there, nr x nf
47% diverged_step 0, or the step at which the field passed
48% S.blowup (the run stops there)
49% notes, backend, loop_s
50%
51% The lab has a WebGPU port of this function (packages/core/src/gpu/fdLoop.ts)
52% that computes the same record from the same S, in single precision. It
53% carries phi as (phi, phi - phim) and c likewise, the same arithmetic
54% rearranged so that the rounding of the field does not feed back into the
55% velocity; keep the two implementations in step.
57t_loop = tic;
58nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; nsteps = S.nsteps; fs = S.fs; na = S.na;
59N = S.N; nr = S.nr;
60KB = S.KB; KG = S.KG; Lm = S.Lm; Mvec = S.Mvec; keep = S.keep;
61R = S.R; AR = S.AR; Krr = S.Krr; Ginv = S.Ginv;
62kb = S.kb; cb = S.cb; gb = S.gb; kt = S.kt; ct = S.ct; gt = S.gt; cs = S.cs; ks = S.ks;
63fstat_s = S.fstat_s; fstat_b = S.fstat_b; fstat_t = S.fstat_t;
64sa = S.sa; F = S.F; ip = S.ip; jp = S.jp; kp = ip + (jp - 1) * ny;
65frame_steps = S.frame_steps; nf = S.nf;
66ln = S.lis_node; lb = S.lis_bin; lw = S.lis_wgt; nbins = S.nbins;
67sup_ptr = S.sup_ptr; sup_nodes = S.sup_nodes; nsup = S.nsup; sup_kind = S.sup_kind;
68c_old = S.c_old; c_new = S.c_new; s1f = S.s1f; LR = S.LR; dt2 = S.dt2; blowup = S.blowup;
69wsamp = S.wsamp;
70gb0 = max(gb, 0); gt0 = max(gt, 0);
71hasb = any(kb > 0); hast = any(kt > 0); hascs = any(cs > 0);
72rmax = max(abs(R), [], 1);
74% State
75phi = S.phi0; phim = phi; c = S.c0; cm = c;
76Lphi_m = Lm * phi;
77w_p = zeros(1, na); wxx_p = zeros(1, na); wyy_p = zeros(1, na); wxy_p = zeros(1, na);
78C = zeros(nr, na); Fsup = zeros(nsup, na);
79Pbins = zeros(nbins, na);
80Wf = zeros(ny, nx, nf); Cf = zeros(nr, nf);
81diverged_step = 0;
82notes = {};
83kf = 1;
85for nstep = 1:nsteps
86 ia = floor((nstep - 1) * dt * fs) + 1;
87 ws = wsamp(ia);
88 Lphi = Lm * phi;
89 if nr > 0
90 w = R * c + phi;
91 v = ((phi - phim) + R * (c - cm)) / dt;
92 else
93 w = phi;
94 v = (phi - phim) / dt;
95 end
97 % The sigma1 loss on the total velocity, as a nodal force: 2 sigma1 M Lm (w - wm) / dt
98 dLw = Lphi - Lphi_m;
99 if nr > 0
100 dLw = dLw + LR * (c - cm);
101 end
102 f = F(nstep) * sa + s1f .* dLw;
104 % Contact and dashpot forces about rest
105 fsS = fstat_s - ks .* w; % springs: their static force plus the departure
106 fsB = fstat_b; fsT = fstat_t;
107 if hasb
108 compb = gb - w;
109 actb = kb > 0 & compb > 0;
110 fb = kb .* max(compb, 0) - kb .* gb0 - (cb .* v) .* actb;
111 f = f + fb; fsB = fsB + fb;
112 end
113 if hast
114 compt = gt + w;
115 actt = kt > 0 & compt > 0;
116 ft = -kt .* max(compt, 0) + kt .* gt0 - (ct .* v) .* actt;
117 f = f + ft; fsT = fsT + ft;
118 end
119 if hascs
120 fl = -cs .* v;
121 f = f + fl; fsS = fsS + fl;
122 end
124 % Series of the current state, box-averaged onto the output samples
125 w_p(ia) = w_p(ia) + w(kp) * ws;
126 wxx_p(ia) = wxx_p(ia) + (phi(kp + ny) - 2 * phi(kp) + phi(kp - ny)) / dx^2 * ws;
127 wyy_p(ia) = wyy_p(ia) + (phi(kp + 1) - 2 * phi(kp) + phi(kp - 1)) / dy^2 * ws;
128 wxy_p(ia) = wxy_p(ia) + (phi(kp + ny + 1) - phi(kp + ny - 1) - phi(kp - ny + 1) + phi(kp - ny - 1)) / (4 * dx * dy) * ws;
129 if nr > 0
130 C(:, ia) = C(:, ia) + c * ws;
131 end
132 for q = 1:nsup
133 nd = sup_nodes(sup_ptr(q):sup_ptr(q + 1) - 1);
134 if sup_kind(q) == 0
135 fq = fsS(nd);
136 elseif sup_kind(q) == 1
137 fq = fsB(nd);
138 else
139 fq = fsT(nd);
140 end
141 Fsup(q, ia) = Fsup(q, ia) + sum(fq) * ws;
142 end
144 % Rigid-body coordinates
145 if nr > 0
146 cacc = Ginv * (R' * f - Krr * c - AR' * phi);
147 cn = c_new * (2 * c - c_old * cm + dt2 * cacc);
148 acc = (f - KB * phi - KG * phi - AR * c) ./ Mvec - R * cacc;
149 else
150 cn = c;
151 acc = (f - KB * phi - KG * phi) ./ Mvec;
152 end
153 phin = c_new * (2 * phi - c_old * phim + dt2 * acc);
154 phin = phin .* keep;
156 % Radiated pressure: acceleration of the total field at this step, weighted and delayed per entry
157 accw = ((phin - phi) - (phi - phim)) / dt2;
158 if nr > 0
159 accw = accw + R * (((cn - c) - (c - cm)) / dt2);
160 end
161 Pbins(:, ia) = Pbins(:, ia) + accumarray(lb, lw .* accw(ln), [nbins 1]) * ws;
163 if kf <= nf && nstep == frame_steps(kf)
164 Wf(:, :, kf) = reshape(phin, ny, nx);
165 if nr > 0
166 Cf(:, kf) = cn;
167 end
168 kf = kf + 1;
169 end
171 peak = max(abs(phin));
172 if nr > 0
173 peak = peak + sum(abs(cn(:)') .* rmax);
174 end
175 if ~(peak < blowup) % also catches NaN and Inf
176 diverged_step = nstep;
177 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
178 break
179 end
180 phim = phi; phi = phin; cm = c; c = cn; Lphi_m = Lphi;
181end
183Rec.w_p = w_p; Rec.wxx_p = wxx_p; Rec.wyy_p = wyy_p; Rec.wxy_p = wxy_p;
184Rec.C = C; Rec.Fsup = Fsup;
185Rec.Wf = Wf; Rec.Cf = Cf;
186Rec.Pbins = Pbins;
187Rec.diverged_step = diverged_step;
188Rec.notes = notes;
189Rec.backend = 'numbl';
190Rec.loop_s = toc(t_loop);
191end