f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 1import { test } from "node:test";
2import assert from "node:assert/strict";
3import { accelerations, potentialEnergy, maxAngularFrequency, type Fields } from "../src/sim/cpuRef.ts";
4import { type Grid, type Material, trapWeight } from "../src/params.ts";
6const mat: Material = { D: 8.66, C: 1.665e8, G: 5.91e7, nu: 0.3, rhoh: 6.241 };
8function fields(n: number): Fields {
9 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
10}
12/** Jacobi eigenvalues of a symmetric matrix (row-major, n x n). */
13function symEig(A: Float64Array, n: number): number[] {
14 const a = Float64Array.from(A);
15 for (let sweep = 0; sweep < 60; sweep++) {
16 let off = 0;
17 for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] ** 2;
18 if (off < 1e-22) break;
19 for (let p = 0; p < n; p++) {
20 for (let q = p + 1; q < n; q++) {
21 const apq = a[p * n + q];
22 if (Math.abs(apq) < 1e-300) continue;
23 const theta = (a[q * n + q] - a[p * n + p]) / (2 * apq);
24 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
25 const c = 1 / Math.sqrt(t * t + 1), s = t * c;
26 for (let k = 0; k < n; k++) {
27 const akp = a[k * n + p], akq = a[k * n + q];
28 a[k * n + p] = c * akp - s * akq;
29 a[k * n + q] = s * akp + c * akq;
30 }
31 for (let k = 0; k < n; k++) {
32 const apk = a[p * n + k], aqk = a[q * n + k];
33 a[p * n + k] = c * apk - s * aqk;
34 a[q * n + k] = s * apk + c * aqk;
35 }
36 }
37 }
38 }
39 const ev: number[] = [];
40 for (let k = 0; k < n; k++) ev.push(a[k * n + k]);
41 return ev.sort((x, y) => x - y);
42}
44/** Mass-symmetrized linear stiffness of one block ("w" or "uv"). */
45function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
46 const N = g.n;
47 const n = block === "w" ? N : 2 * N;
48 const K = new Float64Array(n * n);
49 const x = fields(N), y = fields(N);
50 const eps = 1e-9;
51 const mass = (k: number) => mat.rhoh * g.hx * g.hy * trapWeight(g, k % g.nx, Math.floor(k / g.nx));
52 for (let col = 0; col < n; col++) {
53 x.u.fill(0); x.v.fill(0); x.w.fill(0);
54 if (block === "w") x.w[col] = eps;
55 else if (col < N) x.u[col] = eps;
56 else x.v[col - N] = eps;
57 accelerations(g, mat, x, y);
58 const mc = mass(col % N);
59 for (let row = 0; row < n; row++) {
60 const acc = block === "w" ? y.w[row] : row < N ? y.u[row] : y.v[row - N];
61 const mr = mass(row % N);
62 // K_sym = M^{1/2} (M^{-1} K) M^{-1/2}
63 K[row * n + col] = (-acc / eps) * Math.sqrt(mr / mc);
64 }
65 }
66 return { K, n };
67}
69test("bending stiffness is symmetric and has exactly three rigid modes", () => {
70 const nx = 9, ny = 7;
71 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
72 const { K, n } = stiffness(g, "w");
73 let asym = 0, scale = 0;
74 for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
75 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
76 scale = Math.max(scale, Math.abs(K[r * n + c]));
77 }
78 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
79 const ev = symEig(K, n);
80 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
81 assert.ok(ev[3] > 1e-4 * ev[n - 1], "fourth mode is not rigid");
82});
84test("membrane stiffness is symmetric and has exactly three rigid modes", () => {
85 const nx = 7, ny = 6;
86 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
87 const { K, n } = stiffness(g, "uv");
88 let asym = 0, scale = 0;
89 for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
90 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
91 scale = Math.max(scale, Math.abs(K[r * n + c]));
92 }
93 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym}`);
94 const ev = symEig(K, n);
95 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
96 assert.ok(ev[3] > 1e-4 * ev[n - 1], "fourth mode is not rigid");
97});
99test("free square plate frequencies match the Leissa values", () => {
100 const nx = 21, ny = 21, a = 0.3;
101 const g: Grid = { nx, ny, hx: a / (nx - 1), hy: a / (ny - 1), n: nx * ny };
102 const { K, n } = stiffness(g, "w");
103 const ev = symEig(K, n);
104 const lam2 = ev.slice(3, 8).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
105 // Completely free square plate, nu = 0.3 (Leissa / Gorman): first five nonzero modes.
106 const ref = [13.468, 19.596, 24.27, 34.801, 34.801];
107 for (let k = 0; k < ref.length; k++) {
108 const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
109 assert.ok(relErr < 0.03, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
110 }
111});
113test("forces are the negative gradient of the discrete energy (nonlinear state)", () => {
114 const nx = 8, ny = 7;
115 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
116 const f = fields(g.n), acc = fields(g.n);
117 let seed = 7;
118 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
119 for (let k = 0; k < g.n; k++) {
120 f.w[k] = 0.02 * rnd(); // ~ 25 thicknesses, strongly nonlinear
121 f.u[k] = 1e-4 * rnd();
122 f.v[k] = 1e-4 * rnd();
123 }
124 accelerations(g, mat, f, acc);
125 const delta = 1e-7;
126 let worst = 0;
127 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
128 for (let k = 0; k < g.n; k += 5) {
129 const mass = mat.rhoh * g.hx * g.hy * trapWeight(g, k % nx, Math.floor(k / nx));
130 const orig = arr[k];
131 arr[k] = orig + delta;
132 const Ep = potentialEnergy(g, mat, f);
133 arr[k] = orig - delta;
134 const Em = potentialEnergy(g, mat, f);
135 arr[k] = orig;
136 const Ffd = -(Ep - Em) / (2 * delta);
137 const Fan = accArr[k] * mass;
138 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
139 }
140 }
141 assert.ok(worst < 1e-4, `gradient mismatch ${worst}`);
142});
144test("maximum frequency estimate is close to the in-plane CFL estimate", () => {
145 const g: Grid = { nx: 50, ny: 50, hx: 0.004, hy: 0.004, n: 2500 };
146 const wmax = maxAngularFrequency(g, mat);
147 const cp = Math.sqrt(mat.C / mat.rhoh);
148 const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)
149 assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);
150});
152test("forces are the energy gradient with an initial-shape imperfection", () => {
153 const nx = 8, ny = 7;
154 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
155 const f = fields(g.n), acc = fields(g.n);
156 const w0 = new Float64Array(g.n);
157 let seed = 11;
158 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
159 for (let k = 0; k < g.n; k++) {
160 const x = (k % nx) / (nx - 1) - 0.5;
161 w0[k] = 0.003 * (1 - 4 * x * x); // cylindrical bow, 3 mm
162 f.w[k] = 0.01 * rnd();
163 f.u[k] = 1e-4 * rnd();
164 f.v[k] = 1e-4 * rnd();
165 }
166 accelerations(g, mat, f, acc, w0);
167 const delta = 1e-7;
168 let worst = 0;
169 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
170 for (let k = 0; k < g.n; k += 5) {
171 const mass = mat.rhoh * g.hx * g.hy * trapWeight(g, k % nx, Math.floor(k / nx));
172 const orig = arr[k];
173 arr[k] = orig + delta;
174 const Ep = potentialEnergy(g, mat, f, w0);
175 arr[k] = orig - delta;
176 const Em = potentialEnergy(g, mat, f, w0);
177 arr[k] = orig;
178 const Ffd = -(Ep - Em) / (2 * delta);
179 const Fan = accArr[k] * mass;
180 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
181 }
182 }
183 assert.ok(worst < 1e-4, `gradient mismatch ${worst}`);
184 // The stress-free shape itself carries no membrane energy.
185 const rest = fields(g.n);
186 assert.ok(potentialEnergy(g, mat, rest, w0) < 1e-12, "w0 is stress free");
187});