1// The eigen solver: banded factorization and inertia, Lanczos slices against
2// a dense Jacobi solve on small grids (with degenerate pairs), and the
3// free-plate frequencies against the reference values.
4import { test } from "node:test";
5import assert from "node:assert/strict";
6import { makeGrid, withMask, type Plate } from "../src/params.ts";
7import { uniformPlate } from "../src/model/reference.ts";
8import { assemble, stiffnessEntry, type LinearSystem } from "../src/model/assemble.ts";
9import { factorize, solve, inertia } from "../src/modal/banded.ts";
10import { tridiagonalEigen, sliceModes } from "../src/modal/lanczos.ts";
11import { computeModes, localRunner } from "../src/modal/modes.ts";
13const steel = { E: 200e9, nu: 0.3, rho: 7850 };
15/** All eigenvalues of M^{-1/2} K M^{-1/2} by cyclic Jacobi (small systems only). */
16function denseEigenvalues(sys: LinearSystem): number[] {
17 const n = sys.n;
18 const A: number[][] = [];
19 for (let i = 0; i < n; i++) {
20 A.push([]);
21 for (let j = 0; j < n; j++) A[i].push(stiffnessEntry(sys, i, j) / Math.sqrt(sys.mass[i] * sys.mass[j]));
22 }
23 for (let sweep = 0; sweep < 60; sweep++) {
24 let off = 0;
25 for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) off += A[i][j] ** 2;
26 if (off < 1e-24 * n * n) break;
27 for (let p = 0; p < n; p++) {
28 for (let q = p + 1; q < n; q++) {
29 if (Math.abs(A[p][q]) < 1e-300) continue;
30 const theta = (A[q][q] - A[p][p]) / (2 * A[p][q]);
31 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
32 const c = 1 / Math.sqrt(t * t + 1), s = t * c;
33 for (let k = 0; k < n; k++) {
34 const akp = A[k][p], akq = A[k][q];
35 A[k][p] = c * akp - s * akq; A[k][q] = s * akp + c * akq;
36 }
37 for (let k = 0; k < n; k++) {
38 const apk = A[p][k], aqk = A[q][k];
39 A[p][k] = c * apk - s * aqk; A[q][k] = s * apk + c * aqk;
40 }
41 }
42 }
43 }
44 return Array.from({ length: n }, (_, i) => A[i][i]).sort((a, b) => a - b);
45}
47test("the banded LDL^T solves and counts eigenvalues", () => {
48 const sys = assemble(uniformPlate(makeGrid(0.1, 0.08, 0.01), steel, 1e-3));
49 const shift = -((2 * Math.PI * 50) ** 2);
50 const f = factorize(sys.band, sys.mass, sys.n, sys.bw, shift);
51 assert.equal(f.negative, 0, "K + 50 Hz shift is positive definite");
52 const b = new Float64Array(sys.n).map((_, i) => Math.sin(i));
53 const x = new Float64Array(sys.n);
54 solve(f, b, x);
55 // Check (K - shift M) x = b with the banded matvec.
56 const { stiffnessMatvec } = require_assemble();
57 const y = new Float64Array(sys.n);
58 stiffnessMatvec(sys, x, y);
59 let err = 0, nrm = 0;
60 for (let i = 0; i < sys.n; i++) { const r = y[i] - shift * sys.mass[i] * x[i] - b[i]; err += r * r; nrm += b[i] * b[i]; }
61 assert.ok(Math.sqrt(err / nrm) < 1e-9);
62 const dense = denseEigenvalues(sys);
63 for (const f of [30, 100, 300, 1000]) {
64 const lam = (2 * Math.PI * f) ** 2;
65 const expected = dense.filter((v) => v < lam).length;
66 assert.equal(inertia(sys.band, sys.mass, sys.n, sys.bw, lam), expected, `count below ${f} Hz`);
67 }
68});
70import * as assembleModule from "../src/model/assemble.ts";
71function require_assemble(): typeof assembleModule { return assembleModule; }
73test("the tridiagonal eigensolver agrees with a dense solve", () => {
74 const d = [4, 3, 2, 5, 1, 3, 3], e = [1, 0.5, 2, 0.1, 1, 1];
75 const { values, vectors } = tridiagonalEigen(d, e);
76 const n = d.length;
77 for (let k = 0; k < n; k++) {
78 // T v = lambda v
79 for (let i = 0; i < n; i++) {
80 const tv = d[i] * vectors[k][i] + (i > 0 ? e[i - 1] * vectors[k][i - 1] : 0) + (i < n - 1 ? e[i] * vectors[k][i + 1] : 0);
81 assert.ok(Math.abs(tv - values[k] * vectors[k][i]) < 1e-10);
82 }
83 }
84 for (let k = 1; k < n; k++) assert.ok(values[k] >= values[k - 1]);
85});
87test("all modes of a small free square and a small disc match a dense solve, degenerate pairs included", async () => {
88 const cases: { name: string; pl: Plate }[] = [
89 { name: "square", pl: uniformPlate(makeGrid(0.1, 0.1, 0.1 / 9), steel, 1e-3) },
90 { name: "disc", pl: uniformPlate(withMask(makeGrid(0.12, 0.12, 0.01), (x, y) => x * x + y * y <= 0.06 ** 2), steel, 1e-3) },
91 ];
92 for (const c of cases) {
93 const sys = assemble(c.pl);
94 const dense = denseEigenvalues(sys).filter((v) => v > (2 * Math.PI) ** 2);
95 const fMax = Math.sqrt(dense[dense.length - 1]) / (2 * Math.PI) * 1.01;
96 const ms = await computeModes(sys, fMax, localRunner(sys), { sliceSize: 12 });
97 assert.equal(ms.missing, 0, `${c.name}: missing modes`);
98 assert.equal(ms.count, dense.length, `${c.name}: ${ms.count} modes vs ${dense.length} dense`);
99 for (let m = 0; m < ms.count; m++) {
100 const lam = (2 * Math.PI * ms.freqs[m]) ** 2;
101 assert.ok(Math.abs(lam / dense[m] - 1) < 1e-7, `${c.name}: mode ${m} ${lam} vs ${dense[m]}`);
102 }
103 // Mass-normalized shapes: phi' M phi = 1, and K phi = lambda M phi.
104 const { stiffnessMatvec } = require_assemble();
105 const y = new Float64Array(sys.n);
106 for (const m of [0, 1, Math.floor(ms.count / 2), ms.count - 1]) {
107 const phi = Float64Array.from(ms.shapes.subarray(m * sys.n, (m + 1) * sys.n));
108 let mn = 0;
109 for (let i = 0; i < sys.n; i++) mn += sys.mass[i] * phi[i] * phi[i];
110 assert.ok(Math.abs(mn - 1) < 1e-5, `${c.name}: mass norm ${mn}`);
111 stiffnessMatvec(sys, phi, y);
112 const lam = (2 * Math.PI * ms.freqs[m]) ** 2;
113 let r = 0, s = 0;
114 for (let i = 0; i < sys.n; i++) { r += (y[i] - lam * sys.mass[i] * phi[i]) ** 2; s += (lam * sys.mass[i] * phi[i]) ** 2; }
115 // The stored shapes are float32, which the stiffness amplifies by omega_max^2 / omega^2.
116 assert.ok(Math.sqrt(r / s) < 1e-3, `${c.name}: residual ${Math.sqrt(r / s)}`);
117 }
118 // The float64 shapes of one slice satisfy the eigen equation to rounding.
119 const lo = (2 * Math.PI) ** 2, hi = dense[Math.min(9, dense.length - 1)] * 1.0001;
120 const slice = sliceModes(sys.band, sys.mass, sys.n, sys.bw, lo, hi, inertia(sys.band, sys.mass, sys.n, sys.bw, lo), inertia(sys.band, sys.mass, sys.n, sys.bw, hi));
121 assert.equal(slice.lambda.length, slice.expected);
122 for (let m = 0; m < slice.lambda.length; m++) {
123 stiffnessMatvec(sys, slice.shapes[m], y);
124 let r = 0, s = 0;
125 for (let i = 0; i < sys.n; i++) { r += (y[i] - slice.lambda[m] * sys.mass[i] * slice.shapes[m][i]) ** 2; s += (slice.lambda[m] * sys.mass[i] * slice.shapes[m][i]) ** 2; }
126 assert.ok(Math.sqrt(r / s) < 1e-6, `${c.name}: float64 residual ${Math.sqrt(r / s)}`);
127 }
128 // The disc has degenerate pairs: the first two frequencies agree closely.
129 if (c.name === "disc") assert.ok(Math.abs(ms.freqs[1] / ms.freqs[0] - 1) < 0.02, "degenerate pair");
130 }
131});
133test("free square plate frequencies converge to Gorman's values", async () => {
134 // lambda = omega a^2 sqrt(rho h / D) for the completely free square (nu = 0.3).
135 const gorman = [13.468, 19.596, 24.271, 34.801, 34.801, 61.093, 61.093];
136 const a = 0.3, h = 1e-3;
137 const D = (steel.E * h ** 3) / (12 * (1 - steel.nu ** 2));
138 const scale = a * a * Math.sqrt((steel.rho * h) / D);
139 const errors: number[] = [];
140 for (const cells of [10, 20]) {
141 const sys = assemble(uniformPlate(makeGrid(a, a, a / cells), steel, h));
142 const ms = await computeModes(sys, (80 / scale) / (2 * Math.PI), localRunner(sys), { sliceSize: 20 });
143 let worst = 0;
144 for (let m = 0; m < 7; m++) worst = Math.max(worst, Math.abs((2 * Math.PI * ms.freqs[m] * scale) / gorman[m] - 1));
145 errors.push(worst);
146 }
147 assert.ok(errors[0] < 0.1, `10 cells: ${errors[0]}`);
148 assert.ok(errors[1] < 0.025, `20 cells: ${errors[1]}`);
149 assert.ok(errors[1] < errors[0]);
150});
152test("free circular plate frequencies are close to Leissa's values", async () => {
153 // lambda^2 = omega a^2 sqrt(rho h / D), completely free disc, nu = 0.3: (2,0) 5.253, (0,1) 9.084, (3,0) 12.23, (1,1) 20.52, (4,0) 21.6.
154 const leissa = [5.253, 5.253, 9.084, 12.23, 12.23, 20.52, 20.52, 21.6, 21.6];
155 const a = 0.15, h = 1e-3;
156 const D = (steel.E * h ** 3) / (12 * (1 - steel.nu ** 2));
157 const scale = a * a * Math.sqrt((steel.rho * h) / D);
158 const g = withMask(makeGrid(2 * a, 2 * a, (2 * a) / 40), (x, y) => x * x + y * y <= a * a);
159 const sys = assemble(uniformPlate(g, steel, h));
160 const ms = await computeModes(sys, (25 / scale) / (2 * Math.PI), localRunner(sys), { sliceSize: 20 });
161 for (let m = 0; m < leissa.length; m++) {
162 const lam = 2 * Math.PI * ms.freqs[m] * scale;
163 assert.ok(Math.abs(lam / leissa[m] - 1) < 0.03, `mode ${m}: ${lam} vs ${leissa[m]}`);
164 }
165});
167test("a plate on a stand has no rigid modes and bounces at the pad frequency", async () => {
168 const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
169 const pl: Plate = { ...uniformPlate(g, steel, 2e-3), stand: { x: 0, y: 0, pad: { innerRadius: 0.012, outerRadius: 0.03, thickness: 0.01, E: 0.3e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
170 const sys = assemble(pl);
171 const ms = await computeModes(sys, 400, localRunner(sys), { sliceSize: 20 });
172 const { standRates, massProperties } = await import("../src/params.ts");
173 const sr = standRates(pl)!;
174 const { mass } = massProperties(pl);
175 const bounce = Math.sqrt(sr.totalKw / mass) / (2 * Math.PI);
176 // The lowest modes are the two rocking modes and the bounce; the bounce is close to sqrt(K/M).
177 assert.ok(ms.freqs[0] > 1, "no zero-frequency mode");
178 const near = Array.from(ms.freqs).some((f) => Math.abs(f / bounce - 1) < 0.15);
179 assert.ok(near, `bounce ${bounce} Hz among ${Array.from(ms.freqs.slice(0, 5)).map((f) => f.toFixed(1)).join(", ")}`);
180});