d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 1// Tests of the spatial discretization on the float64 reference (no GPU needed).
3import { test } from "node:test";
4import assert from "node:assert/strict";
5import { accelerations, addDamping, fields, maxAngularFrequency, potentialEnergy, type Fields } from "../src/sim/cpuRef.ts";
6import { type Grid, type Material, countInside, nodeWeight, withEllipseMask } from "../src/params.ts";
8// 0.79 mm steel with nu = 0.3, the value the reference tables use.
9const mat: Material = { D: 8.66, C: 1.665e8, G: 5.91e7, nu: 0.3, rhoh: 6.241 };
11function rectangle(nx: number, ny: number, Lx = 0.3, Ly = 0.2): Grid {
12 return { nx, ny, hx: Lx / (nx - 1), hy: Ly / (ny - 1), n: nx * ny };
13}
15function lcg(seed: number): () => number {
16 return () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
17}
19function massOf(g: Grid, k: number): number {
20 return mat.rhoh * g.hx * g.hy * nodeWeight(g, k % g.nx, Math.floor(k / g.nx));
21}
23/** Eigenvalues of a symmetric matrix (row-major, n x n) by cyclic Jacobi. */
24function symEig(A: Float64Array, n: number): number[] {
25 const a = Float64Array.from(A);
26 for (let sweep = 0; sweep < 60; sweep++) {
27 let off = 0;
28 for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] ** 2;
29 if (off < 1e-22) break;
30 for (let p = 0; p < n; p++) {
31 for (let q = p + 1; q < n; q++) {
32 const apq = a[p * n + q];
33 if (Math.abs(apq) < 1e-300) continue;
34 const theta = (a[q * n + q] - a[p * n + p]) / (2 * apq);
35 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
36 const c = 1 / Math.sqrt(t * t + 1), s = t * c;
37 for (let k = 0; k < n; k++) {
38 const akp = a[k * n + p], akq = a[k * n + q];
39 a[k * n + p] = c * akp - s * akq;
40 a[k * n + q] = s * akp + c * akq;
41 }
42 for (let k = 0; k < n; k++) {
43 const apk = a[p * n + k], aqk = a[q * n + k];
44 a[p * n + k] = c * apk - s * aqk;
45 a[q * n + k] = s * apk + c * aqk;
46 }
47 }
48 }
49 }
50 const ev: number[] = [];
51 for (let k = 0; k < n; k++) ev.push(a[k * n + k]);
52 return ev.sort((x, y) => x - y);
53}
55/** Mass-symmetrized linear stiffness M^{-1/2} K M^{-1/2} of one block ("w" or "uv") over the nodes of the sheet. */
56function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
57 const nodes: number[] = [];
58 for (let k = 0; k < g.n; k++) if (massOf(g, k) > 0) nodes.push(k);
59 const N = nodes.length;
60 const n = block === "w" ? N : 2 * N;
61 const K = new Float64Array(n * n);
62 const x = fields(g.n), y = fields(g.n);
63 const eps = 1e-9;
64 for (let col = 0; col < n; col++) {
65 x.u.fill(0); x.v.fill(0); x.w.fill(0);
66 if (block === "w") x.w[nodes[col]] = eps;
67 else if (col < N) x.u[nodes[col]] = eps;
68 else x.v[nodes[col - N]] = eps;
69 accelerations(g, mat, x, y);
70 const mc = massOf(g, nodes[col % N]);
71 for (let row = 0; row < n; row++) {
72 const acc = block === "w" ? y.w[nodes[row]] : row < N ? y.u[nodes[row]] : y.v[nodes[row - N]];
73 K[row * n + col] = (-acc / eps) * Math.sqrt(massOf(g, nodes[row % N]) / mc);
74 }
75 }
76 return { K, n };
77}
79function checkSymmetricWithThreeRigidModes(K: Float64Array, n: number): number[] {
80 let asym = 0, scale = 0;
81 for (let r = 0; r < n; r++) {
82 for (let c = 0; c < n; c++) {
83 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
84 scale = Math.max(scale, Math.abs(K[r * n + c]));
85 }
86 }
87 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
88 const ev = symEig(K, n);
89 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
90 assert.ok(ev[3] > 1e-4 * ev[n - 1], `fourth mode is not rigid (${ev[3]} vs ${ev[n - 1]})`);
91 return ev;
92}
94/** Largest relative mismatch between the forces and the centerd-difference gradient of the energy. */
95function gradientMismatch(g: Grid, f: Fields, w0?: Float64Array): number {
96 const acc = fields(g.n);
97 accelerations(g, mat, f, acc, w0);
98 const delta = 1e-7;
99 let worst = 0;
100 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
101 for (let k = 0; k < g.n; k += 5) {
102 const mass = massOf(g, k);
103 if (mass === 0) continue;
104 const orig = arr[k];
105 arr[k] = orig + delta;
106 const Ep = potentialEnergy(g, mat, f, w0);
107 arr[k] = orig - delta;
108 const Em = potentialEnergy(g, mat, f, w0);
109 arr[k] = orig;
110 const Ffd = -(Ep - Em) / (2 * delta);
111 const Fan = accArr[k] * mass;
112 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
113 }
114 }
115 return worst;
116}
118/** A random state with w of about 25 thicknesses (strongly nonlinear) and a cylindrical 3 mm bow. */
119function randomState(g: Grid, seed: number): { f: Fields; w0: Float64Array } {
120 const rnd = lcg(seed);
121 const f = fields(g.n);
122 const w0 = new Float64Array(g.n);
123 for (let k = 0; k < g.n; k++) {
124 const x = (k % g.nx) / (g.nx - 1) - 0.5;
125 w0[k] = 0.003 * (1 - 4 * x * x);
126 f.w[k] = 0.02 * rnd();
127 f.u[k] = 1e-4 * rnd();
128 f.v[k] = 1e-4 * rnd();
129 }
130 return { f, w0 };
131}
133test("bending stiffness is symmetric and has exactly three rigid modes", () => {
134 const { K, n } = stiffness(rectangle(9, 7), "w");
135 checkSymmetricWithThreeRigidModes(K, n);
136});
138test("membrane stiffness is symmetric and has exactly three rigid modes", () => {
139 const { K, n } = stiffness(rectangle(7, 6), "uv");
140 checkSymmetricWithThreeRigidModes(K, n);
141});
143test("free square plate frequencies converge to the reference values", () => {
144 // Completely free square plate, nu = 0.3, lambda^2 = omega a^2 sqrt(rho h / D)
145 // (Gorman 1978; Leissa's Ritz values are slightly higher). The scheme is
146 // second order: the errors are about 0.8 percent at 20 cells per side and
147 // 0.3 percent at 30 (the pair at 34.80 is degenerate).
148 const ref = [13.468, 19.596, 24.271, 34.801, 34.801, 61.093, 61.093];
149 const a = 0.3;
150 const errors: number[][] = [];
151 for (const N of [16, 21]) {
152 const { K, n } = stiffness(rectangle(N, N, a, a), "w");
153 const ev = symEig(K, n);
154 const lam2 = ev.slice(3, 3 + ref.length).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
155 errors.push(lam2.map((v, k) => Math.abs(v - ref[k]) / ref[k]));
156 }
157 for (let k = 0; k < ref.length; k++) {
158 assert.ok(errors[1][k] < 0.03, `mode ${k}: error ${(errors[1][k] * 100).toFixed(2)}% at 20 cells`);
159 assert.ok(errors[1][k] < 0.7 * errors[0][k], `mode ${k}: error does not shrink from 15 to 20 cells`);
160 }
161});
163test("forces are the negative gradient of the discrete energy, flat and bowed", () => {
164 const g = rectangle(8, 7);
165 const { f, w0 } = randomState(g, 7);
166 assert.ok(gradientMismatch(g, f) < 1e-4, "flat");
167 assert.ok(gradientMismatch(g, f, w0) < 1e-4, "bowed");
168 // The stress-free shape itself carries no energy.
169 assert.ok(potentialEnergy(g, mat, fields(g.n), w0) < 1e-12, "w0 is stress free");
170});
172test("loss terms only remove energy", () => {
173 const g = rectangle(9, 8);
174 const rnd = lcg(3);
175 const v = fields(g.n), acc = fields(g.n);
176 for (let k = 0; k < g.n; k++) { v.u[k] = rnd(); v.v[k] = rnd(); v.w[k] = rnd(); }
177 addDamping(g, 1.0, 1e-3, 5.0, v, acc);
178 let power = 0;
179 for (let k = 0; k < g.n; k++) power += massOf(g, k) * (v.u[k] * acc.u[k] + v.v[k] * acc.v[k] + v.w[k] * acc.w[k]);
180 assert.ok(power < 0, `power ${power} should be negative`);
181});
183test("maximum frequency estimate is close to the in-plane CFL estimate", () => {
184 const g: Grid = { nx: 50, ny: 50, hx: 0.004, hy: 0.004, n: 2500 };
185 const wmax = maxAngularFrequency(g, mat);
186 const cp = Math.sqrt(mat.C / mat.rhoh);
187 const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)
188 assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);
189});
191// ---------- Masked (elliptical) grids ----------
193test("elliptical sheet: both stiffness blocks are symmetric with exactly three rigid modes", () => {
194 const g = withEllipseMask(rectangle(13, 10));
195 assert.ok(countInside(g) < g.n && countInside(g) > g.n / 2, "mask removes some nodes");
196 for (let k = 0; k < g.n; k++) {
197 // Every node of the sheet carries mass (no whiskers left after the mask is trimmed).
198 if (g.mask![k]) assert.ok(massOf(g, k) > 0, `node ${k} has no mass`);
199 }
200 const b = stiffness(g, "w");
201 checkSymmetricWithThreeRigidModes(b.K, b.n);
202 const m = stiffness(g, "uv");
203 checkSymmetricWithThreeRigidModes(m.K, m.n);
204});
206test("elliptical sheet: forces are the energy gradient, flat and bowed", () => {
207 const g = withEllipseMask(rectangle(12, 9));
208 const { f, w0 } = randomState(g, 5);
209 assert.ok(gradientMismatch(g, f) < 1e-4, "flat");
210 assert.ok(gradientMismatch(g, f, w0) < 1e-4, "bowed");
211});
213test("free circular plate frequencies approach the reference values", () => {
214 const nx = 25, a = 0.15; // radius; 24 cells across the diameter
215 const g = withEllipseMask({ nx, ny: nx, hx: (2 * a) / (nx - 1), hy: (2 * a) / (nx - 1), n: nx * nx });
216 const { K, n } = stiffness(g, "w");
217 const ev = symEig(K, n);
218 const lam2 = ev.slice(3, 12).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
219 // Completely free circular plate, nu = 0.3 (Leissa, Table 2.4): modes (2,0) x2, (0,1),
220 // (3,0) x2, (1,1) x2, (4,0) x2 by nodal diameters and circles. The staircase edge is
221 // first-order accurate: about 2 percent here, under 1 percent at 60 cells across.
222 const ref = [5.358, 5.358, 9.003, 12.439, 12.439, 20.475, 20.475, 21.835, 21.835];
223 for (let k = 0; k < ref.length; k++) {
224 const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
225 assert.ok(relErr < 0.04, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
226 }
227});