1import { describe, expect, it } from "vitest";
2import {
3 FD_METHOD,
4 SCENARIO_PRESETS,
5 buildPlanBundle,
6 cloneSpec,
7 parseMethodPlan,
8} from "../src/index";
9import { maxAbs, runScript, runner, scalar, solve, testSpec, vector } from "./helpers";
11describe("the hanging plate", () => {
12 it("gets the compound-pendulum stiffness from the prestress operator to rounding", async () => {
13 const spec = testSpec();
14 const out = await runScript(
15 spec,
16 [
17 "mdl = plate_model(params);",
18 "d = norm([mdl.mnt.xc mdl.mnt.yc] - mdl.mnt.pivot);",
19 "out.k_swing = mdl.R(:, 2)' * mdl.Kg * mdl.R(:, 2);",
20 "out.mgd = mdl.mnt.mass * mdl.P.g * d;",
21 "out.k_spin = mdl.R(:, 3)' * mdl.Kg * mdl.R(:, 3);",
22 "out.k_heave = mdl.R(:, 1)' * mdl.Kg * mdl.R(:, 1);",
23 "out.mass_check = mdl.ps.mass_check;",
24 "out.kT = mdl.Krr(1, 1);",
25 "out.T_over_L = mdl.mnt.mass * mdl.P.g / params.mount.length;",
26 "out.hz = mdl.rigid_hz(:)';",
27 "out.G = mdl.G(:)'; out.Krr = mdl.Krr(:)';",
28 ].join("\n"),
29 );
30 const kSwing = scalar(out, "k_swing");
31 const mgd = scalar(out, "mgd");
32 expect(Math.abs(kSwing - mgd) / mgd).toBeLessThan(1e-10);
33 expect(Math.abs(scalar(out, "k_spin"))).toBeLessThan(1e-9 * mgd);
34 expect(Math.abs(scalar(out, "k_heave"))).toBeLessThan(1e-9 * mgd);
35 expect(scalar(out, "mass_check")).toBeCloseTo(1, 10);
36 expect(scalar(out, "kT")).toBeCloseTo(scalar(out, "T_over_L"), 8);
37 // the two coupled pendulum modes of thread and plate, from the same 3 x 3 system
38 const G = vector(out, "G");
39 const K = vector(out, "Krr");
40 const a = K[4] * G[0] - K[0] * G[4] * 0; // placeholder to keep the arrays used
41 void a;
42 const g11 = G[0],
43 g12 = G[1],
44 g22 = G[4];
45 const k11 = K[0],
46 k22 = K[4];
47 // det(K - w^2 G) = 0 for the 2 x 2 sway-swing block
48 const A = g11 * g22 - g12 * g12;
49 const B = -(k11 * g22 + k22 * g11);
50 const C = k11 * k22;
51 const w2a = (-B - Math.sqrt(B * B - 4 * A * C)) / (2 * A);
52 const w2b = (-B + Math.sqrt(B * B - 4 * A * C)) / (2 * A);
53 const hz = Array.from(vector(out, "hz")).sort((p, q) => p - q);
54 expect(hz[1]).toBeCloseTo(Math.sqrt(w2a) / (2 * Math.PI), 6);
55 expect(hz[2]).toBeCloseTo(Math.sqrt(w2b) / (2 * Math.PI), 6);
56 // a 300 x 200 mm steel sheet on a 300 mm thread: sway near 0.8 Hz, swing near 2.9 Hz
57 expect(hz[1]).toBeGreaterThan(0.7);
58 expect(hz[1]).toBeLessThan(0.9);
59 expect(hz[2]).toBeGreaterThan(2.6);
60 expect(hz[2]).toBeLessThan(3.2);
61 });
63 it("hangs with its centre of mass below an off-axis hole", async () => {
64 const spec = testSpec();
65 spec.mount = {
66 type: "thread",
67 x: 0.05,
68 y: 0.19,
69 length: 0.3,
70 thread:
71 spec.mount.type === "thread"
72 ? spec.mount.thread
73 : { name: "nylon", diameter: 0.5e-3, G: 1e9 },
74 };
75 const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
76 expect(run.ok).toBe(true);
77 if (!run.ok) return;
78 const plan = parseMethodPlan(run.value);
79 // the hole snaps to the nearest grid node
80 const hx = Math.round(0.05 / plan.dx) * plan.dx;
81 const hy = Math.round(0.19 / plan.dy) * plan.dy;
82 const dx = plan.xc - hx;
83 const dy = plan.yc - hy;
84 const n = Math.hypot(dx, dy);
85 expect(plan.gx).toBeCloseTo(dx / n, 6);
86 expect(plan.gy).toBeCloseTo(dy / n, 6);
87 expect(plan.gz).toBe(0);
88 expect(plan.tension.length).toBe(1);
89 expect(plan.tension[0]).toBeCloseTo(plan.mass * 9.81, 6);
90 expect(plan.sag).toBe(0);
91 });
93 it("is unstable when hung from below its centre of mass, and says so", async () => {
94 const spec = testSpec();
95 spec.mount = {
96 type: "thread",
97 x: 0.15,
98 y: 0.01,
99 length: 0.3,
100 thread: { name: "nylon", diameter: 0.5e-3, G: 1e9 },
101 };
102 // the plate hangs from the hole in any case: plate_mount turns it so the mass is below the hole
103 const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
104 expect(run.ok).toBe(true);
105 if (!run.ok) return;
106 const plan = parseMethodPlan(run.value);
107 expect(plan.gy).toBeCloseTo(1, 6);
108 expect(plan.rigid_hz.every((h) => Number.isFinite(h))).toBe(true);
109 });
111 it("shares the weight between two threads and restrains the yaw", async () => {
112 const spec = testSpec();
113 spec.mount = { type: "threads", x1: 0.06, y1: 0.19, x2: 0.24, y2: 0.19, length: 0.3 };
114 const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
115 expect(run.ok).toBe(true);
116 if (!run.ok) return;
117 const plan = parseMethodPlan(run.value);
118 expect(plan.tension.length).toBe(2);
119 expect(plan.tension[0] + plan.tension[1]).toBeCloseTo(plan.mass * 9.81, 6);
120 expect(plan.tension[0]).toBeCloseTo(plan.tension[1], 6);
121 expect(plan.rigid_hz.every((h) => Number.isFinite(h) && h > 0.3)).toBe(true);
122 expect(plan.rigid_names[2]).toMatch(/yaw/);
123 });
125 it("refuses two threads with the centre of mass outside them", async () => {
126 const spec = testSpec();
127 spec.mount = { type: "threads", x1: 0.2, y1: 0.19, x2: 0.28, y2: 0.19, length: 0.3 };
128 const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
129 expect(run.ok).toBe(false);
130 if (!run.ok) expect(run.error).toMatch(/slack/);
131 });
133 it("releases from a pull of the asked size and swings", async () => {
134 const spec = testSpec();
135 spec.excitation = { type: "release", x: 0.15, y: 0.02, displacement: 0.002, radius: 0.015 };
136 spec.probe = { x: 0.15, y: 0.02 };
137 spec.fd.duration = 0.001;
138 const { result } = await solve(spec);
139 expect(result.info.stable).toBe(true);
140 // the probe node is near the pull point; the total departure there starts near the pull
141 // (the probe sits a row above the pull's centre, where the swing is a little smaller)
142 expect(Math.abs(result.probe.w[0])).toBeGreaterThan(0.0015);
143 expect(Math.abs(result.probe.w[0])).toBeLessThan(0.0022);
144 // the swing coordinate carries most of it and the spin none
145 const c = result.rigid.series;
146 expect(Math.abs(c[1][0])).toBeGreaterThan(0.005);
147 expect(Math.abs(c[2][0])).toBeLessThan(1e-9);
148 expect(result.info.notes.some((n) => /released/.test(n))).toBe(true);
149 });
150});
152describe("plates on supports", () => {
153 it("carries the weight on the pads at rest, sagging into them", async () => {
154 const spec = testSpec("sheet-on-foam");
155 const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
156 expect(run.ok).toBe(true);
157 if (!run.ok) return;
158 const plan = parseMethodPlan(run.value);
159 expect(plan.gz).toBe(-1);
160 expect(plan.sag).toBeGreaterThan(0);
161 const { result } = await solve(spec);
162 const total = result.supports.static.reduce((a, b) => a + b, 0);
163 expect(total).toBeCloseTo(result.static.mass * 9.81, 6);
164 expect(result.supports.static.every((f) => f > 0)).toBe(true);
165 expect(result.static.sag).toBeGreaterThan(1e-4); // soft foam: a fraction of a millimetre
166 expect(result.info.stable).toBe(true);
167 });
169 it("lifts off the pads under a hard strike and lands again", async () => {
170 const spec = testSpec("sheet-on-foam");
171 spec.excitation = {
172 type: "impulse",
173 x: 0.15,
174 y: 0.1,
175 force: 300,
176 duration: 1e-3,
177 radius: 0.02,
178 };
179 spec.fd.duration = 0.12;
180 const { result } = await solve(spec);
181 expect(result.info.stable).toBe(true);
182 const f = result.supports.series[0];
183 const n = f.length;
184 const early = f.subarray(0, Math.floor(n / 3));
185 const late = f.subarray(Math.floor((2 * n) / 3));
186 expect(Math.min(...Array.from(early))).toBeLessThan(1e-6); // off the pad
187 expect(Math.max(...Array.from(late))).toBeGreaterThan(result.supports.static[0] * 0.5); // back on it
188 });
190 it("holds the cymbal between two preloaded felts", async () => {
191 const spec = testSpec("cymbal-on-stand");
192 spec.fd.nx = 21;
193 spec.fd.ny = 21;
194 const { result } = await solve(spec);
195 expect(result.supports.names).toEqual(["lower felt", "upper felt"]);
196 const [lower, upper] = result.supports.static;
197 expect(lower).toBeGreaterThan(0);
198 expect(upper).toBeLessThan(0);
199 expect(lower + upper).toBeCloseTo(result.static.mass * 9.81, 6);
200 expect(result.info.stable).toBe(true);
201 });
203 it("sags like a cantilever beam in the vise", async () => {
204 const spec = testSpec("cantilever-in-vise");
205 spec.fd.nx = 31;
206 spec.fd.ny = 7;
207 spec.excitation = { type: "impulse", x: 0.2, y: 0.025, force: 1, duration: 1e-3, radius: 0.01 };
208 const out = await runScript(
209 spec,
210 [
211 "mdl = plate_model(params);",
212 "w0 = reshape(mdl.st.w0, mdl.ny, mdl.nx);",
213 "out.tip = w0(round((mdl.ny + 1) / 2), end);",
214 "held = mdl.mnt.held; out.xclamp = max(mdl.dom.X(held));",
215 "out.mu = mdl.P.mu; out.D = mdl.P.D; out.E = mdl.P.material.E; out.h = params.plate.h; out.Lx = params.plate.Lx;",
216 ].join("\n"),
217 );
218 const Lfree = scalar(out, "Lx") - scalar(out, "xclamp");
219 const q = scalar(out, "mu") * 9.81; // load per unit area
220 const E = scalar(out, "E");
221 const h = scalar(out, "h");
222 // beam theory for a narrow strip: w = q L^4 / (8 E I / b) with I / b = h^3 / 12
223 const beam = (q * Lfree ** 4) / (8 * E * (h ** 3 / 12));
224 const tip = -scalar(out, "tip");
225 expect(tip).toBeGreaterThan(0.85 * beam);
226 expect(tip).toBeLessThan(1.15 * beam);
227 });
229 it("has no rigid-body coordinates in the vise and three elsewhere", async () => {
230 const vise = await solve(testSpec("cantilever-in-vise"));
231 expect(vise.result.rigid.nr).toBe(0);
232 expect(vise.result.rigid.series.length).toBe(0);
233 const hand = await solve(testSpec("disk-in-hand"));
234 expect(hand.result.rigid.nr).toBe(3);
235 expect(hand.result.rigid.hz.every((h) => Number.isFinite(h) && h >= 0)).toBe(true);
236 expect(hand.result.rigid.hz.filter((h) => h > 0).length).toBeGreaterThanOrEqual(2);
237 });
238});
240describe("the scheme", () => {
241 it("keeps the swing of a struck hanging plate bounded and its split exact", async () => {
242 const spec = testSpec();
243 spec.fd.duration = 0.45;
244 spec.damping.sigma0 = 0;
245 spec.damping.sigma1 = 0;
246 spec.output.nframes = 10;
247 const out = await runScript(
248 spec,
249 [
250 "S = plate_fd_setup(params); mdl = S.mdl;",
251 "Rec = plate_fd_loop(S);",
252 "E = zeros(3, S.nf); for k = 1:S.nf, E(:, k) = mdl.Q' * reshape(Rec.Wf(:, :, k), [], 1); end",
253 "out.leak = max(max(abs(E(1:2, :))));",
254 "out.swing = max(abs(Rec.C(2, :))); out.sway = max(abs(Rec.C(1, :)));",
255 "n = size(Rec.C, 2); out.early = max(abs(Rec.C(2, 1:floor(n / 3)))); out.late = max(abs(Rec.C(2, floor(2 * n / 3):end)));",
256 "% the angular impulse of the strike about the hole, and the swing it can give at most",
257 "J = trapz(S.tt, S.F); rd = (params.excitation.x - mdl.mnt.pivot(1)) * mdl.mnt.ed(1) + (params.excitation.y - mdl.mnt.pivot(2)) * mdl.mnt.ed(2);",
258 "out.bound = abs(J * rd) / mdl.G(2, 2) / (2 * pi * min(mdl.rigid_hz(mdl.rigid_hz > 0.1)));",
259 ].join("\n"),
260 );
261 expect(scalar(out, "leak")).toBeLessThan(1e-12);
262 // the swing stays of the order the impulse allows, and does not grow
263 expect(scalar(out, "swing")).toBeLessThan(3 * scalar(out, "bound"));
264 expect(scalar(out, "late")).toBeLessThan(
265 1.5 * Math.max(scalar(out, "early"), scalar(out, "swing") * 0.5),
266 );
267 });
269 it("factors every operator so the port can apply it as stencils and weights", async () => {
270 const out = await runScript(
271 testSpec(),
272 [
273 "S = plate_fd_setup(params);",
274 "u = S.keep; Pz = spdiags(u, 0, S.N, S.N);",
275 "KBf = Pz * (S.Gs' * S.WB * S.Gs) * Pz; out.eKB = full(max(max(abs(KBf - S.KB)))) / full(max(max(abs(S.KB))));",
276 "KGf = Pz * (S.Gs' * S.WG * S.Gs) * Pz + spdiags(S.ks .* u, 0, S.N, S.N); out.eKG = full(max(max(abs(KGf - S.KG)))) / max(full(max(max(abs(S.KG)))), eps);",
277 "Lf = spdiags(S.inva, 0, S.N, S.N) * (S.Gs' * S.WL * S.Gs) * Pz; out.eL = full(max(max(abs(Lf - S.Lm)))) / full(max(max(abs(S.Lm))));",
278 "[gi, gj, gv] = find(S.Gs); out.eInt = max(abs(gv ./ S.Gscale(gi) - round(gv ./ S.Gscale(gi))));",
279 ].join("\n"),
280 );
281 expect(scalar(out, "eKB")).toBeLessThan(1e-12);
282 expect(scalar(out, "eKG")).toBeLessThan(1e-12);
283 expect(scalar(out, "eL")).toBeLessThan(1e-12);
284 expect(scalar(out, "eInt")).toBeLessThan(1e-9);
285 });
287 it("plans a time step at the stability bound", async () => {
288 const spec = testSpec();
289 const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
290 expect(run.ok).toBe(true);
291 if (!run.ok) return;
292 const plan = parseMethodPlan(run.value);
293 const { Lx, Ly, h } = spec.plate;
294 const { E, nu, rho } = spec.material;
295 const kappa = Math.sqrt((E * h ** 3) / (12 * (1 - nu ** 2)) / (rho * h));
296 const dx = Lx / (spec.fd.nx - 1);
297 const dy = Ly / (spec.fd.ny - 1);
298 const lam = 4 / dx ** 2 + 4 / dy ** 2;
299 const s1 = spec.damping.sigma1;
300 const bending = 2 / (lam * (s1 + Math.sqrt(kappa ** 2 + s1 ** 2)));
301 // the prestress and the thread spring lower the bound by a negligible amount
302 expect(plan.dt_max).toBeLessThanOrEqual(bending);
303 expect(plan.dt_max).toBeGreaterThan(0.999 * bending);
304 expect(plan.dt).toBeCloseTo(spec.fd.cfl * plan.dt_max, 12);
305 expect(plan.cell_steps).toBe(plan.ncells * plan.nsteps);
306 });
308 it("is stable just below the CFL limit and diverges just above it", async () => {
309 const spec = testSpec();
310 spec.fd.cfl = 0.98;
311 const stable = await solve(spec);
312 expect(stable.result.info.stable).toBe(true);
313 expect(stable.result.probe.w.every(Number.isFinite)).toBe(true);
315 spec.fd.cfl = 1.03;
316 const blown = await solve(spec);
317 expect(blown.result.info.stable).toBe(false);
318 expect(blown.result.info.diverged_step).toBeGreaterThan(0);
319 expect(blown.result.probe.w.every(Number.isFinite)).toBe(true);
320 expect(blown.result.frames.w.data.every((v) => Number.isFinite(v) || Number.isNaN(v))).toBe(
321 true,
322 );
323 expect(blown.result.info.notes.some((n) => /diverged/.test(n))).toBe(true);
324 });
326 it("runs every scenario on a coarse grid", async () => {
327 for (const sc of SCENARIO_PRESETS) {
328 const spec = testSpec(sc.id);
329 const { result } = await solve(spec);
330 expect(result.info.stable, sc.id).toBe(true);
331 expect(maxAbs(result.probe.w), sc.id).toBeGreaterThan(0);
332 expect(result.audio.p.every(Number.isFinite), sc.id).toBe(true);
333 expect(result.frames.w.nf, sc.id).toBe(result.frames.t.length);
334 expect(result.supports.series.length, sc.id).toBe(result.supports.names.length);
335 }
336 }, 240_000);
338 it("reports its problem size and a series at the sample rate", async () => {
339 const spec = testSpec();
340 const { result } = await solve(spec);
341 expect(result.info.ncells).toBe(spec.fd.nx * spec.fd.ny);
342 expect(result.info.extra.cell_steps).toBe(result.info.ncells * result.info.nsteps);
343 const expected = Math.floor(result.info.nsteps * result.info.dt * spec.output.fs);
344 expect(result.t.length).toBe(expected);
345 expect(result.audio.p.length).toBe(expected);
346 expect(result.audio.model).toBe("unbaffled");
347 });
349 it("gives the same physics with the baffled listener except for the sound", async () => {
350 const a = testSpec();
351 const b = cloneSpec(a);
352 b.listener.model = "baffled";
353 const ra = await solve(a);
354 const rb = await solve(b);
355 expect(ra.result.probe.w).toEqual(rb.result.probe.w);
356 expect(maxAbs(ra.result.audio.p)).toBeLessThan(maxAbs(rb.result.audio.p) * 1.5);
357 expect(ra.result.audio.p).not.toEqual(rb.result.audio.p);
358 });
359});