1import { describe, expect, it } from "vitest";
2import {
3 audibleFreqs,
4 isolatedAudible,
5 maxAbs,
6 reference,
7 referenceWav,
8 runScript,
9 scalar,
10 solve,
11 testSpec,
12} from "./helpers";
14/**
15 * The finite-element method against the diagnostics of
16 * aluminum_plate_fem_sim.py (sim-docs) at the meshes the reference was run
17 * at: the same element, assembly, thread spring, damping, coupled contact
18 * and radiation model. The eigensolver differs (subspace iteration here,
19 * ARPACK there) and the contact integrator differs (fixed-step Runge-Kutta
20 * here, DOP853 there), so the agreement is to the accuracy of both, which
21 * is far tighter than the physics.
22 */
23function checkAgainst(
24 name: string,
25 result: Awaited<ReturnType<typeof solve>>["result"],
26 wavTol: number,
27) {
28 const ref = reference(name);
29 const f = audibleFreqs(result);
30 expect(f.length).toBe(ref.modes.length);
31 ref.modes.forEach((m, k) => {
32 expect(Math.abs(f[k] - m.frequency_hz) / m.frequency_hz).toBeLessThan(2e-7);
33 });
34 expect(result.info.n_rigid).toBe(ref.derived.low_frequency_momentum_modes);
35 expect(result.info.ndof).toBe(ref.derived.fem_degrees_of_freedom);
36 for (const k of isolatedAudible(result)) {
37 const m = ref.modes[k - result.info.n_rigid];
38 expect(result.modes.q[k]).toBeCloseTo(m.q_total, 2);
39 expect(Math.abs(result.modes.hit[k])).toBeCloseTo(Math.abs(m.hit_participation), 3);
40 }
41 const c = result.contact;
42 const d = ref.derived;
43 const rel = (x: number, y: number) => Math.abs(x - y) / Math.abs(y);
44 expect(rel(c.tc, d.contact_duration_s)).toBeLessThan(1e-5);
45 expect(rel(c.impulse, d.impact_impulse_Ns)).toBeLessThan(1e-4);
46 expect(rel(c.peak_force, d.peak_contact_force_N)).toBeLessThan(1e-5);
47 expect(rel(c.hammer_exit_v, d.hammer_exit_velocity_m_s_positive_into_plate)).toBeLessThan(1e-5);
48 expect(
49 rel(c.plate_exit_v, d.plate_contact_point_exit_velocity_m_s_positive_into_plate),
50 ).toBeLessThan(1e-5);
51 expect(c.effective_restitution).toBeCloseTo(d.effective_relative_restitution, 5);
52 const wav = referenceWav(name);
53 const a = result.audio.p;
54 const g = result.audio.gain;
55 let worst = 0;
56 for (let i = 0; i < wav.samples.length; i++)
57 worst = Math.max(worst, Math.abs(a[i] * g - wav.samples[i]));
58 expect(worst).toBeLessThan(wavTol);
59}
61describe("finite elements against the reference simulator", () => {
62 it("reproduces the centre hit on a 16 x 16 mesh", async () => {
63 const { result } = await solve(testSpec("center-hit"), "fem");
64 // degenerate pairs mix differently between the two eigensolvers, which moves the sound by a little
65 checkAgainst("fem_center_16", result, 5e-4);
66 expect(result.contact.hammer_v.length).toBe(result.contact.t.length);
67 expect(result.contact.tc).toBeGreaterThan(result.contact.hertz_duration);
68 expect(result.contact.peak_force).toBeLessThan(result.contact.hertz_peak_force / 4);
69 });
71 it("reproduces the centre hit on a 24 x 24 mesh", async () => {
72 const spec = testSpec("center-hit");
73 spec.fem.mesh_x = 24;
74 spec.fem.mesh_y = 24;
75 const { result } = await solve(spec, "fem");
76 checkAgainst("fem_center_24", result, 5e-4);
77 });
79 it("reproduces the off-centre hit on a 32 x 32 mesh", async () => {
80 const spec = testSpec("off-center-hit");
81 spec.fem.mesh_x = 32;
82 spec.fem.mesh_y = 32;
83 const { result } = await solve(spec, "fem");
84 checkAgainst("fem_off_32", result, 5e-4);
85 expect(Math.abs(result.modes.hit[result.info.n_rigid])).toBeGreaterThan(1);
86 });
88 it("converges the main partial from above as the mesh is refined", async () => {
89 const at = async (n: number) => {
90 const spec = testSpec();
91 spec.fem.mesh_x = n;
92 spec.fem.mesh_y = n;
93 spec.render.duration = 0.01;
94 const { result } = await solve(spec, "fem");
95 return audibleFreqs(result)[2];
96 };
97 const f16 = await at(16);
98 const f24 = await at(24);
99 expect(f16).toBeGreaterThan(f24);
100 expect(f24).toBeGreaterThan(reference("rr_center").modes[2].frequency_hz);
101 });
103 it("finds the rigid-body modes and mass-normalizes every mode", async () => {
104 const spec = testSpec();
105 const out = await runScript(
106 spec,
107 [
108 "P = plate_setup(params);",
109 "A = plate_fem_assemble(params, P);",
110 "[lam, V, info] = plate_fem_eigs(A.K, A.M, 12, 30, 60, 1e-10);",
111 "out.converged = double(info.converged);",
112 "out.residual = info.residual;",
113 "out.iterations = info.iterations;",
114 "out.lam1 = lam(1); out.lam3 = lam(3); out.lam4 = lam(4);",
115 "out.norm4 = V(:, 4)' * (A.M * V(:, 4));",
116 "out.orth = abs(V(:, 4)' * (A.M * V(:, 5)));",
117 "out.pendulum_hz = sqrt(max(lam(3), 0)) / (2 * pi);",
118 ].join("\n"),
119 );
120 expect(scalar(out, "converged")).toBe(1);
121 // eigenvalues converge quadratically in the residual: a 1e-10 eigenvalue tolerance goes with residuals near 1e-5
122 expect(scalar(out, "residual")).toBeLessThan(1e-4);
123 expect(scalar(out, "iterations")).toBeLessThan(40);
124 // two free rotations at zero, then the pendulum, then the first bending mode
125 expect(Math.abs(scalar(out, "lam1"))).toBeLessThan(1);
126 expect(scalar(out, "pendulum_hz")).toBeGreaterThan(1);
127 expect(scalar(out, "pendulum_hz")).toBeLessThan(5);
128 expect(Math.sqrt(scalar(out, "lam4")) / (2 * Math.PI)).toBeCloseTo(3636.156, 2);
129 expect(scalar(out, "norm4")).toBeCloseTo(1, 8);
130 expect(scalar(out, "orth")).toBeLessThan(1e-8);
131 });
133 it("carries the mode shapes, curvatures and frames the views draw", async () => {
134 const spec = testSpec();
135 spec.render.duration = 0.02;
136 const { result } = await solve(spec, "fem");
137 expect(result.grid.x.length).toBe(spec.fem.mesh_x + 1);
138 expect(result.modes.shapes.ny).toBe(spec.fem.mesh_y + 1);
139 expect(result.modes.shapes.nf).toBe(result.modes.n);
140 expect(result.frames.w_total.nf).toBe(spec.render.nframes);
141 expect(maxAbs(result.frames.vm.data)).toBeGreaterThan(0);
142 // the rigid modes carry no curvature
143 const n = result.modes.shapes.nx * result.modes.shapes.ny;
144 let kmax = 0;
145 for (let i = 0; i < n; i++) kmax = Math.max(kmax, Math.abs(result.frames.w_total.data[i]));
146 expect(result.info.n_rigid).toBe(3);
147 expect(kmax).toBeGreaterThanOrEqual(0);
148 });
149});