tests/cpuRef.test.ts at main - proof-of-concept/steelsim41

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, uniformPlate, type Fields } from "../src/sim/cpuRef.ts";

6import { type Grid, type Material, type Plate, countInside, insideNode, membraneStiffness, mountPatch, nodeMass, nodeX, nodeY, withEllipseMask, withMask } from "../src/params.ts";

8// Steel with nu = 0.3, the value the reference tables use; 0.79 mm thick.

9const mat: Material = { E: 170e9, nu: 0.3, rho: 7900 };

10const H = 0.79e-3;

12function rectangle(nx: number, ny: number, Lx = 0.3, Ly = 0.2): Grid {

13 return { nx, ny, hx: Lx / (nx - 1), hy: Ly / (ny - 1), n: nx * ny };

14}

16function lcg(seed: number): () => number {

17 return () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);

18}

20function massOf(pl: Plate, k: number): number {

21 return nodeMass(pl, k % pl.grid.nx, Math.floor(k / pl.grid.nx));

22}

24/** A plate whose thickness rises linearly from 0.5 to 1.5 mm across the width. */

25function taperedPlate(g: Grid, w0?: Float64Array, mount: Plate["mount"] = null): Plate {

26 const pl = uniformPlate(g, mat, H, w0);

27 for (let j = 0; j < g.ny; j++) {

28 for (let i = 0; i < g.nx; i++) {

29 if (insideNode(g, i, j)) pl.h[j * g.nx + i] = 1e-3 * (1 + nodeX(g, i) / ((g.nx - 1) * g.hx));

30 }

31 }

32 return { ...pl, mount };

33}

35/** Eigenvalues of a symmetric matrix (row-major, n x n) by cyclic Jacobi. */

36function symEig(A: Float64Array, n: number): number[] {

37 const a = Float64Array.from(A);

38 for (let sweep = 0; sweep < 60; sweep++) {

39 let off = 0;

40 for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] ** 2;

41 if (off < 1e-22) break;

42 for (let p = 0; p < n; p++) {

43 for (let q = p + 1; q < n; q++) {

44 const apq = a[p * n + q];

45 if (Math.abs(apq) < 1e-300) continue;

46 const theta = (a[q * n + q] - a[p * n + p]) / (2 * apq);

47 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));

48 const c = 1 / Math.sqrt(t * t + 1), s = t * c;

49 for (let k = 0; k < n; k++) {

50 const akp = a[k * n + p], akq = a[k * n + q];

51 a[k * n + p] = c * akp - s * akq;

52 a[k * n + q] = s * akp + c * akq;

53 }

54 for (let k = 0; k < n; k++) {

55 const apk = a[p * n + k], aqk = a[q * n + k];

56 a[p * n + k] = c * apk - s * aqk;

57 a[q * n + k] = s * apk + c * aqk;

58 }

59 }

60 }

61 }

62 const ev: number[] = [];

63 for (let k = 0; k < n; k++) ev.push(a[k * n + k]);

64 return ev.sort((x, y) => x - y);

65}

67/** Mass-symmetrized linear stiffness M^{-1/2} K M^{-1/2} of one block ("w" or "uv") over the nodes of the plate. */

68function stiffness(pl: Plate, block: "w" | "uv"): { K: Float64Array; n: number } {

69 const g = pl.grid;

70 const nodes: number[] = [];

71 for (let k = 0; k < g.n; k++) if (massOf(pl, k) > 0) nodes.push(k);

72 const N = nodes.length;

73 const n = block === "w" ? N : 2 * N;

74 const K = new Float64Array(n * n);

75 const x = fields(g.n), y = fields(g.n);

76 const eps = 1e-9;

77 for (let col = 0; col < n; col++) {

78 x.u.fill(0); x.v.fill(0); x.w.fill(0);

79 if (block === "w") x.w[nodes[col]] = eps;

80 else if (col < N) x.u[nodes[col]] = eps;

81 else x.v[nodes[col - N]] = eps;

82 accelerations(pl, x, y);

83 const mc = massOf(pl, nodes[col % N]);

84 for (let row = 0; row < n; row++) {

85 const acc = block === "w" ? y.w[nodes[row]] : row < N ? y.u[nodes[row]] : y.v[nodes[row - N]];

86 K[row * n + col] = (-acc / eps) * Math.sqrt(massOf(pl, nodes[row % N]) / mc);

87 }

88 }

89 return { K, n };

90}

92function checkSymmetric(K: Float64Array, n: number): void {

93 let asym = 0, scale = 0;

94 for (let r = 0; r < n; r++) {

95 for (let c = 0; c < n; c++) {

96 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));

97 scale = Math.max(scale, Math.abs(K[r * n + c]));

98 }

99 }

100 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);

101}

102

103/** `gap` is the smallest ratio of the fourth eigenvalue to the largest that counts as non-rigid (a tapered plate on a fine grid has softer low modes). */

104function checkSymmetricWithThreeRigidModes(K: Float64Array, n: number, gap = 1e-4): number[] {

105 checkSymmetric(K, n);

106 const ev = symEig(K, n);

107 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");

108 assert.ok(ev[3] > gap * ev[n - 1], `fourth mode is not rigid (${ev[3]} vs ${ev[n - 1]})`);

109 return ev;

110}

111

112/** Largest relative mismatch between the forces and the centered-difference gradient of the energy. */

113function gradientMismatch(pl: Plate, f: Fields): number {

114 const g = pl.grid;

115 const acc = fields(g.n);

116 accelerations(pl, f, acc);

117 const delta = 1e-7;

118 let worst = 0;

119 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {

120 for (let k = 0; k < g.n; k += 5) {

121 const mass = massOf(pl, k);

122 if (mass === 0) continue;

123 const orig = arr[k];

124 arr[k] = orig + delta;

125 const Ep = potentialEnergy(pl, f);

126 arr[k] = orig - delta;

127 const Em = potentialEnergy(pl, f);

128 arr[k] = orig;

129 const Ffd = -(Ep - Em) / (2 * delta);

130 const Fan = accArr[k] * mass;

131 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));

132 }

133 }

134 return worst;

135}

136

137/** A random state with w of about 25 thicknesses (strongly nonlinear) and a cylindrical 3 mm bow. */

138function randomState(g: Grid, seed: number): { f: Fields; w0: Float64Array } {

139 const rnd = lcg(seed);

140 const f = fields(g.n);

141 const w0 = new Float64Array(g.n);

142 for (let k = 0; k < g.n; k++) {

143 const x = (k % g.nx) / (g.nx - 1) - 0.5;

144 w0[k] = 0.003 * (1 - 4 * x * x);

145 f.w[k] = 0.02 * rnd();

146 f.u[k] = 1e-4 * rnd();

147 f.v[k] = 1e-4 * rnd();

148 }

149 return { f, w0 };

150}

151

152test("bending stiffness is symmetric and has exactly three rigid modes", () => {

153 const { K, n } = stiffness(uniformPlate(rectangle(9, 7), mat, H), "w");

154 checkSymmetricWithThreeRigidModes(K, n);

155});

156

157test("membrane stiffness is symmetric and has exactly three rigid modes", () => {

158 const { K, n } = stiffness(uniformPlate(rectangle(7, 6), mat, H), "uv");

159 checkSymmetricWithThreeRigidModes(K, n);

160});

161

162test("free square plate frequencies converge to the reference values", () => {

163 // Completely free square plate, nu = 0.3, lambda^2 = omega a^2 sqrt(rho h / D)

164 // (Gorman 1978; Leissa's Ritz values are slightly higher). The scheme is

165 // second order: the errors are about 0.8 percent at 20 cells per side and

166 // 0.3 percent at 30 (the pair at 34.80 is degenerate).

167 const ref = [13.468, 19.596, 24.271, 34.801, 34.801, 61.093, 61.093];

168 const a = 0.3;

169 const D = (mat.E * H ** 3) / (12 * (1 - mat.nu ** 2));

170 const errors: number[][] = [];

171 for (const N of [16, 21]) {

172 const { K, n } = stiffness(uniformPlate(rectangle(N, N, a, a), mat, H), "w");

173 const ev = symEig(K, n);

174 const lam2 = ev.slice(3, 3 + ref.length).map((e) => Math.sqrt(e) * a * a * Math.sqrt((mat.rho * H) / D));

175 errors.push(lam2.map((v, k) => Math.abs(v - ref[k]) / ref[k]));

176 }

177 for (let k = 0; k < ref.length; k++) {

178 assert.ok(errors[1][k] < 0.03, `mode ${k}: error ${(errors[1][k] * 100).toFixed(2)}% at 20 cells`);

179 assert.ok(errors[1][k] < 0.7 * errors[0][k], `mode ${k}: error does not shrink from 15 to 20 cells`);

180 }

181});

182

183test("forces are the negative gradient of the discrete energy, flat and bowed", () => {

184 const g = rectangle(8, 7);

185 const { f, w0 } = randomState(g, 7);

186 assert.ok(gradientMismatch(uniformPlate(g, mat, H), f) < 1e-4, "flat");

187 const bowed = uniformPlate(g, mat, H, w0);

188 assert.ok(gradientMismatch(bowed, f) < 1e-4, "bowed");

189 // The stress-free shape itself carries no energy.

190 assert.ok(potentialEnergy(bowed, fields(g.n)) < 1e-12, "w0 is stress free");

191});

192

193test("loss terms only remove energy", () => {

194 const g = rectangle(9, 8);

195 const pl = { ...uniformPlate(g, mat, H), loss: { sigma0: 1.0, sigma1: 1e-3, sigmaIn: 5.0 } };

196 const rnd = lcg(3);

197 const v = fields(g.n), acc = fields(g.n);

198 for (let k = 0; k < g.n; k++) { v.u[k] = rnd(); v.v[k] = rnd(); v.w[k] = rnd(); }

199 addDamping(pl, v, acc);

200 let power = 0;

201 for (let k = 0; k < g.n; k++) power += massOf(pl, k) * (v.u[k] * acc.u[k] + v.v[k] * acc.v[k] + v.w[k] * acc.w[k]);

202 assert.ok(power < 0, `power ${power} should be negative`);

203});

204

205test("maximum frequency estimate is close to the in-plane CFL estimate", () => {

206 const g: Grid = { nx: 50, ny: 50, hx: 0.004, hy: 0.004, n: 2500 };

207 const wmax = maxAngularFrequency(g, mat, H);

208 const cp = Math.sqrt(membraneStiffness(mat, H) / (mat.rho * H));

209 const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)

210 assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);

211 // The in-plane wave speed does not depend on the thickness, so neither does the estimate.

212 const wmax2 = maxAngularFrequency(g, mat, 2 * H);

213 assert.ok(Math.abs(wmax2 - wmax) < 0.05 * wmax, `wmax ${wmax} vs ${wmax2} at twice the thickness`);

214});

215

216// ---------- Variable thickness ----------

217

218test("tapered plate: both stiffness blocks are symmetric with exactly three rigid modes", () => {

219 const pl = taperedPlate(rectangle(9, 7));

220 const b = stiffness(pl, "w");

221 checkSymmetricWithThreeRigidModes(b.K, b.n);

222 const m = stiffness(taperedPlate(rectangle(7, 6)), "uv");

223 checkSymmetricWithThreeRigidModes(m.K, m.n);

224});

225

226test("tapered plate: forces are the energy gradient, flat and bowed", () => {

227 const g = rectangle(8, 7);

228 const { f, w0 } = randomState(g, 11);

229 assert.ok(gradientMismatch(taperedPlate(g), f) < 1e-4, "flat");

230 assert.ok(gradientMismatch(taperedPlate(g, w0), f) < 1e-4, "bowed");

231});

232

233test("a uniform plate of twice the thickness has twice the bending frequencies", () => {

234 // Bending frequencies scale with h, membrane frequencies do not.

235 const thin = stiffness(uniformPlate(rectangle(9, 8), mat, H), "w");

236 const thick = stiffness(uniformPlate(rectangle(9, 8), mat, 2 * H), "w");

237 const e1 = symEig(thin.K, thin.n), e2 = symEig(thick.K, thick.n);

238 for (let k = 3; k < 10; k++) assert.ok(Math.abs(e2[k] / e1[k] - 4) < 1e-6, `mode ${k}: ratio ${e2[k] / e1[k]}`);

239 const thinM = stiffness(uniformPlate(rectangle(6, 5), mat, H), "uv");

240 const thickM = stiffness(uniformPlate(rectangle(6, 5), mat, 2 * H), "uv");

241 const m1 = symEig(thinM.K, thinM.n), m2 = symEig(thickM.K, thickM.n);

242 for (let k = 3; k < 10; k++) assert.ok(Math.abs(m2[k] / m1[k] - 1) < 1e-6, `membrane mode ${k}: ratio ${m2[k] / m1[k]}`);

243});

244

245// ---------- Mount ----------

246

247test("mounted plate: the stiffness is symmetric and the rigid modes feel the mount", () => {

248 const g = rectangle(9, 7);

249 const pl = taperedPlate(g, undefined, { x: 0, y: 0, radius: 0.04, stiffness: 1e4, damping: 0 });

250 assert.ok(mountPatch(pl)!.nodes.length > 1, "the mount covers several nodes");

251 for (const block of ["w", "uv"] as const) {

252 const { K, n } = stiffness(pl, block);

253 checkSymmetric(K, n);

254 }

255 // A uniform translation of the plate in any direction feels the whole mount stiffness.

256 const x = fields(g.n), y = fields(g.n);

257 for (const comp of ["u", "v", "w"] as const) {

258 x.u.fill(0); x.v.fill(0); x.w.fill(0);

259 x[comp].fill(1e-6);

260 accelerations(pl, x, y);

261 let F = 0;

262 for (let k = 0; k < g.n; k++) F += massOf(pl, k) * y[comp][k];

263 assert.ok(Math.abs(-F / 1e-6 - 1e4) < 1e-6 * 1e4, `${comp}: mount force ${-F / 1e-6} vs 1e4`);

264 }

265 // A tilt about the x axis feels a restoring moment from the finite mount patch.

266 x.u.fill(0); x.v.fill(0);

267 for (let k = 0; k < g.n; k++) x.w[k] = 1e-6 * nodeY(g, Math.floor(k / g.nx));

268 accelerations(pl, x, y);

269 let M = 0;

270 for (let k = 0; k < g.n; k++) M += massOf(pl, k) * y.w[k] * nodeY(g, Math.floor(k / g.nx));

271 assert.ok(M < 0, `tilt moment ${M} should restore`);

272});

273

274test("mounted plate: forces are the energy gradient and the dashpot removes energy", () => {

275 const g = rectangle(8, 7);

276 const { f, w0 } = randomState(g, 13);

277 const pl = { ...taperedPlate(g, w0, { x: 0.02, y: -0.01, radius: 0.03, stiffness: 5e3, damping: 2 }), loss: { sigma0: 0, sigma1: 0, sigmaIn: 0 } };

278 assert.ok(gradientMismatch(pl, f) < 1e-4, "gradient");

279 const rnd = lcg(5);

280 const v = fields(g.n), acc = fields(g.n);

281 for (let k = 0; k < g.n; k++) { v.u[k] = rnd(); v.v[k] = rnd(); v.w[k] = rnd(); }

282 addDamping(pl, v, acc);

283 let power = 0, touched = 0;

284 for (let k = 0; k < g.n; k++) {

285 power += massOf(pl, k) * (v.u[k] * acc.u[k] + v.v[k] * acc.v[k] + v.w[k] * acc.w[k]);

286 if (acc.w[k] !== 0) touched++;

287 }

288 assert.ok(power < 0, `power ${power} should be negative`);

289 assert.ok(touched > 0 && touched < countInside(g), "the dashpot acts on the mount nodes only");

290});

291

292// ---------- Masked grids ----------

293

294test("elliptical plate: both stiffness blocks are symmetric with exactly three rigid modes", () => {

295 const g = withEllipseMask(rectangle(13, 10));

296 const pl = uniformPlate(g, mat, H);

297 assert.ok(countInside(g) < g.n && countInside(g) > g.n / 2, "mask removes some nodes");

298 for (let k = 0; k < g.n; k++) {

299 // Every node of the plate carries mass (no whiskers left after the mask is trimmed).

300 if (g.mask![k]) assert.ok(massOf(pl, k) > 0, `node ${k} has no mass`);

301 }

302 const b = stiffness(pl, "w");

303 checkSymmetricWithThreeRigidModes(b.K, b.n);

304 const m = stiffness(pl, "uv");

305 checkSymmetricWithThreeRigidModes(m.K, m.n);

306});

307

308test("disc with a center hole: symmetric stiffness with three rigid modes, forces are the energy gradient", () => {

309 const nx = 15, a = 0.1;

310 const g0: Grid = { nx, ny: nx, hx: (2 * a) / (nx - 1), hy: (2 * a) / (nx - 1), n: nx * nx };

311 const g = withMask(g0, (x, y) => Math.hypot(x, y) <= a && Math.hypot(x, y) >= 0.02);

312 assert.ok(!insideNode(g, 7, 7), "the center node lies in the hole");

313 const pl = taperedPlate(g);

314 const b = stiffness(pl, "w");

315 checkSymmetricWithThreeRigidModes(b.K, b.n, 1e-5);

316 const { f } = randomState(g, 9);

317 assert.ok(gradientMismatch(pl, f) < 1e-4, "gradient");

318});

319

320test("free circular plate frequencies approach the reference values", () => {

321 const nx = 25, a = 0.15; // radius; 24 cells across the diameter

322 const g = withEllipseMask({ nx, ny: nx, hx: (2 * a) / (nx - 1), hy: (2 * a) / (nx - 1), n: nx * nx });

323 const { K, n } = stiffness(uniformPlate(g, mat, H), "w");

324 const ev = symEig(K, n);

325 const D = (mat.E * H ** 3) / (12 * (1 - mat.nu ** 2));

326 const lam2 = ev.slice(3, 12).map((e) => Math.sqrt(e) * a * a * Math.sqrt((mat.rho * H) / D));

327 // Completely free circular plate, nu = 0.3 (Leissa, Table 2.4): modes (2,0) x2, (0,1),

328 // (3,0) x2, (1,1) x2, (4,0) x2 by nodal diameters and circles. The staircase edge is

329 // first-order accurate: about 2 percent here, under 1 percent at 60 cells across.

330 const ref = [5.358, 5.358, 9.003, 12.439, 12.439, 20.475, 20.475, 21.835, 21.835];

331 for (let k = 0; k < ref.length; k++) {

332 const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];

333 assert.ok(relErr < 0.04, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);

334 }

335});

336

337test("node coordinates are centered on the bounding box", () => {

338 const g = rectangle(9, 7, 0.4, 0.3);

339 assert.ok(Math.abs(nodeX(g, 0) + 0.2) < 1e-12 && Math.abs(nodeX(g, 8) - 0.2) < 1e-12);

340 assert.ok(Math.abs(nodeY(g, 0) + 0.15) < 1e-12 && Math.abs(nodeY(g, 6) - 0.15) < 1e-12);

341});