import { describe, expect, it } from "vitest"; import { audibleFreqs, isolatedAudible, maxAbs, reference, referenceWav, runScript, scalar, solve, testSpec, } from "./helpers"; /** * The finite-element method against the diagnostics of * aluminum_plate_fem_sim.py (sim-docs) at the meshes the reference was run * at: the same element, assembly, thread spring, damping, coupled contact * and radiation model. The eigensolver differs (subspace iteration here, * ARPACK there) and the contact integrator differs (fixed-step Runge-Kutta * here, DOP853 there), so the agreement is to the accuracy of both, which * is far tighter than the physics. */ function checkAgainst( name: string, result: Awaited>["result"], wavTol: number, ) { const ref = reference(name); const f = audibleFreqs(result); expect(f.length).toBe(ref.modes.length); ref.modes.forEach((m, k) => { expect(Math.abs(f[k] - m.frequency_hz) / m.frequency_hz).toBeLessThan(2e-7); }); expect(result.info.n_rigid).toBe(ref.derived.low_frequency_momentum_modes); expect(result.info.ndof).toBe(ref.derived.fem_degrees_of_freedom); for (const k of isolatedAudible(result)) { const m = ref.modes[k - result.info.n_rigid]; expect(result.modes.q[k]).toBeCloseTo(m.q_total, 2); expect(Math.abs(result.modes.hit[k])).toBeCloseTo(Math.abs(m.hit_participation), 3); } const c = result.contact; const d = ref.derived; const rel = (x: number, y: number) => Math.abs(x - y) / Math.abs(y); expect(rel(c.tc, d.contact_duration_s)).toBeLessThan(1e-5); expect(rel(c.impulse, d.impact_impulse_Ns)).toBeLessThan(1e-4); expect(rel(c.peak_force, d.peak_contact_force_N)).toBeLessThan(1e-5); expect(rel(c.hammer_exit_v, d.hammer_exit_velocity_m_s_positive_into_plate)).toBeLessThan(1e-5); expect( rel(c.plate_exit_v, d.plate_contact_point_exit_velocity_m_s_positive_into_plate), ).toBeLessThan(1e-5); expect(c.effective_restitution).toBeCloseTo(d.effective_relative_restitution, 5); const wav = referenceWav(name); const a = result.audio.p; const g = result.audio.gain; let worst = 0; for (let i = 0; i < wav.samples.length; i++) worst = Math.max(worst, Math.abs(a[i] * g - wav.samples[i])); expect(worst).toBeLessThan(wavTol); } describe("finite elements against the reference simulator", () => { it("reproduces the centre hit on a 16 x 16 mesh", async () => { const { result } = await solve(testSpec("center-hit"), "fem"); // degenerate pairs mix differently between the two eigensolvers, which moves the sound by a little checkAgainst("fem_center_16", result, 5e-4); expect(result.contact.hammer_v.length).toBe(result.contact.t.length); expect(result.contact.tc).toBeGreaterThan(result.contact.hertz_duration); expect(result.contact.peak_force).toBeLessThan(result.contact.hertz_peak_force / 4); }); it("reproduces the centre hit on a 24 x 24 mesh", async () => { const spec = testSpec("center-hit"); spec.fem.mesh_x = 24; spec.fem.mesh_y = 24; const { result } = await solve(spec, "fem"); checkAgainst("fem_center_24", result, 5e-4); }); it("reproduces the off-centre hit on a 32 x 32 mesh", async () => { const spec = testSpec("off-center-hit"); spec.fem.mesh_x = 32; spec.fem.mesh_y = 32; const { result } = await solve(spec, "fem"); checkAgainst("fem_off_32", result, 5e-4); expect(Math.abs(result.modes.hit[result.info.n_rigid])).toBeGreaterThan(1); }); it("converges the main partial from above as the mesh is refined", async () => { const at = async (n: number) => { const spec = testSpec(); spec.fem.mesh_x = n; spec.fem.mesh_y = n; spec.render.duration = 0.01; const { result } = await solve(spec, "fem"); return audibleFreqs(result)[2]; }; const f16 = await at(16); const f24 = await at(24); expect(f16).toBeGreaterThan(f24); expect(f24).toBeGreaterThan(reference("rr_center").modes[2].frequency_hz); }); it("finds the rigid-body modes and mass-normalizes every mode", async () => { const spec = testSpec(); const out = await runScript( spec, [ "P = plate_setup(params);", "A = plate_fem_assemble(params, P);", "[lam, V, info] = plate_fem_eigs(A.K, A.M, 12, 30, 60, 1e-10);", "out.converged = double(info.converged);", "out.residual = info.residual;", "out.iterations = info.iterations;", "out.lam1 = lam(1); out.lam3 = lam(3); out.lam4 = lam(4);", "out.norm4 = V(:, 4)' * (A.M * V(:, 4));", "out.orth = abs(V(:, 4)' * (A.M * V(:, 5)));", "out.pendulum_hz = sqrt(max(lam(3), 0)) / (2 * pi);", ].join("\n"), ); expect(scalar(out, "converged")).toBe(1); // eigenvalues converge quadratically in the residual: a 1e-10 eigenvalue tolerance goes with residuals near 1e-5 expect(scalar(out, "residual")).toBeLessThan(1e-4); expect(scalar(out, "iterations")).toBeLessThan(40); // two free rotations at zero, then the pendulum, then the first bending mode expect(Math.abs(scalar(out, "lam1"))).toBeLessThan(1); expect(scalar(out, "pendulum_hz")).toBeGreaterThan(1); expect(scalar(out, "pendulum_hz")).toBeLessThan(5); expect(Math.sqrt(scalar(out, "lam4")) / (2 * Math.PI)).toBeCloseTo(3636.156, 2); expect(scalar(out, "norm4")).toBeCloseTo(1, 8); expect(scalar(out, "orth")).toBeLessThan(1e-8); }); it("carries the mode shapes, curvatures and frames the views draw", async () => { const spec = testSpec(); spec.render.duration = 0.02; const { result } = await solve(spec, "fem"); expect(result.grid.x.length).toBe(spec.fem.mesh_x + 1); expect(result.modes.shapes.ny).toBe(spec.fem.mesh_y + 1); expect(result.modes.shapes.nf).toBe(result.modes.n); expect(result.frames.w_total.nf).toBe(spec.render.nframes); expect(maxAbs(result.frames.vm.data)).toBeGreaterThan(0); // the rigid modes carry no curvature const n = result.modes.shapes.nx * result.modes.shapes.ny; let kmax = 0; for (let i = 0; i < n; i++) kmax = Math.max(kmax, Math.abs(result.frames.w_total.data[i])); expect(result.info.n_rigid).toBe(3); expect(kmax).toBeGreaterThanOrEqual(0); }); });