81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 1import { describe, expect, it } from "vitest";
2import {
3 LIBRARY_FILES,
4 buildCompareBundle,
5 isPlainStruct,
6 methodBlockers,
7 getMethod,
8 parseCompareResult,
9 validateSpec,
10} from "../src/index";
11import { maxAbs, runner, solve, testSpec } from "./helpers";
13/** Run a MATLAB snippet against the library and return the struct it leaves in `out`. */
14async function script(source: string) {
15 const run = await runner.run({
16 bundle: {
17 files: [{ name: "c.m", source }, ...LIBRARY_FILES],
18 entry: "c.m",
19 resultVariable: "out",
20 },
21 });
22 if (!run.ok) throw new Error(`${run.error}\n${run.output.join("")}`);
23 if (!isPlainStruct(run.value)) throw new Error("unexpected shape");
24 return run.value.fields as Record<string, number>;
25}
27function ellipseSpec() {
28 const s = testSpec();
29 s.plate.shape = "ellipse";
30 return s;
31}
33describe("the plate's footprint on the grid", () => {
34 it("is the whole grid for a rectangle and a staircase for an ellipse, with every node carrying area", async () => {
35 const out = await script(`
36 r = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), 19, 13);
37 out.r_inside = nnz(r.inside); out.r_interior = nnz(r.interior); out.r_boundary = nnz(r.boundary);
38 out.r_m_corner = r.m(1, 1); out.r_m_edge = r.m(1, 7); out.r_m_in = r.m(5, 7); out.r_full3 = nnz(r.full3);
39 e = plate_domain(struct('shape', 'ellipse', 'Lx', 0.3, 'Ly', 0.2), 31, 21);
40 out.e_inside = nnz(e.inside); out.e_area = nnz(e.inside) / (31 * 21);
41 out.e_min_m = min(e.m(e.inside)); out.e_max_m = max(e.m(:));
42 out.e_corner = e.inside(1, 1); out.e_centre = e.inside(11, 16);
43 out.e_near_ok = all(e.full3(e.near(e.inside)));
44 out.e_kappa_min = min(e.kappa(e.boundary)); out.e_kappa_max = max(e.kappa(e.boundary));
45 nrm = sqrt(e.normal_x(e.boundary).^2 + e.normal_y(e.boundary).^2);
46 out.e_normal_unit = max(abs(nrm - 1));
47 % a sliver: the disk on a coarse grid still has full stencils
48 d = plate_domain(struct('shape', 'ellipse', 'Lx', 1, 'Ly', 1), 9, 9);
49 out.d_full3 = nnz(d.full3);
50 `);
51 expect(out.r_inside).toBe(19 * 13);
52 expect(out.r_interior).toBe(17 * 11);
53 expect(out.r_boundary).toBe(19 * 13 - 17 * 11);
54 expect(out.r_m_corner).toBe(0.25);
55 expect(out.r_m_edge).toBe(0.5);
56 expect(out.r_m_in).toBe(1);
57 expect(out.r_full3).toBe(17 * 11);
58 // pi/4 of the box, to within the staircase
59 expect(Math.abs(out.e_area - Math.PI / 4)).toBeLessThan(0.05);
60 expect(out.e_min_m).toBeGreaterThanOrEqual(0.25);
61 expect(out.e_max_m).toBe(1);
62 expect(out.e_corner).toBe(false);
63 expect(out.e_centre).toBe(true);
64 expect(out.e_near_ok).toBe(true);
65 // curvature of a 300 x 200 ellipse lies between b/a^2 and a/b^2
66 expect(out.e_kappa_min).toBeGreaterThan(0.1 / 0.15 ** 2 - 1e-9);
67 expect(out.e_kappa_max).toBeLessThan(0.15 / 0.1 ** 2 + 1e-9);
68 expect(out.e_normal_unit).toBeLessThan(1e-12);
69 expect(out.d_full3).toBeGreaterThan(0);
70 });
72 it("keeps the strike and the probe on the plate", () => {
73 const spec = ellipseSpec();
74 expect(validateSpec(spec).some((i) => i.severity === "error")).toBe(false);
75 // inside the box, outside the ellipse
76 spec.excitation.x = 0.01;
77 spec.excitation.y = 0.01;
78 const issues = validateSpec(spec);
79 expect(issues.some((i) => i.path === "excitation" && /ellipse/.test(i.message))).toBe(true);
80 });
81});
83describe("circular plate frequencies", () => {
84 // omega a^2 sqrt(rho h / D) for the disk of radius a, nu = 0.3, from the Bessel-function
85 // frequency equations (clamped: J_n I_{n+1} + I_n J_{n+1} = 0; simply supported:
86 // J_{n+1}/J_n + I_{n+1}/I_n = 2 lambda / (1 - nu)) and, for free edges, the tabulated values
87 // (Leissa, Vibration of Plates, 1969). The staircase edge is first order, so a 25-node grid
88 // is expected within a few percent.
89 it("clamped: 10.216", async () => {
90 const out = await script(`
91 dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
92 [~, K, M] = plate_matrices(dom, 'clamped', 0.3);
93 ev = sort(eig(full(K)));
94 out.f1 = sqrt(ev(1)); out.f2 = sqrt(ev(2)); out.f3 = sqrt(ev(3));
95 out.ratio = ev(end) / (4/dom.dx^2 + 4/dom.dy^2)^2;
96 out.n = size(K, 1);
97 `);
98 expect(Math.abs(out.f1 - 10.216) / 10.216).toBeLessThan(0.02);
99 // the (1, 0) pair, degenerate on the true disk
100 expect(Math.abs(out.f2 - 21.26) / 21.26).toBeLessThan(0.03);
101 expect(Math.abs(out.f3 - 21.26) / 21.26).toBeLessThan(0.03);
102 // the explicit step bound derived for the interior stencil still holds
103 expect(out.ratio).toBeLessThan(1);
104 expect(out.ratio).toBeGreaterThan(0.9);
105 });
107 it("simply supported: 4.935 with the curvature term, the paradox value 5.783 without it", async () => {
108 const out = await script(`
109 dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
110 [L, K] = plate_matrices(dom, 'simply-supported', 0.3);
111 ev = sort(eig(full(K)));
112 out.f1 = sqrt(ev(1)); out.f2 = sqrt(ev(2));
113 out.minev = ev(1) / ev(end);
114 ev0 = sort(eig(full(L * L)));
115 out.g1 = sqrt(ev0(1));
116 out.symm = norm(full(K - K'), 1) / norm(full(K), 1);
117 `);
118 expect(out.symm).toBeLessThan(1e-12);
119 expect(out.minev).toBeGreaterThan(0);
120 expect(Math.abs(out.f1 - 4.935) / 4.935).toBeLessThan(0.03);
121 expect(Math.abs(out.f2 - 13.898) / 13.898).toBeLessThan(0.03);
122 // L*L alone converges to the plate with zero Laplacian on the edge, J_0's first zero squared
123 expect(Math.abs(out.g1 - 5.783) / 5.783).toBeLessThan(0.03);
124 expect(out.g1 / out.f1).toBeGreaterThan(1.12);
125 });
127 it("free: 5.358 and 9.003, with exactly three rigid-body modes", async () => {
128 const out = await script(`
129 dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
130 [L, K, M] = plate_matrices(dom, 'free', 0.3);
131 Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh; S = (S + S') / 2;
132 ev = sort(eig(S));
133 out.zero3 = max(abs(ev(1:3))) / ev(end);
134 out.fourth = ev(4) / ev(end);
135 out.f1 = sqrt(ev(4)); out.f2 = sqrt(ev(5)); out.f3 = sqrt(ev(6));
136 out.ratio = ev(end) / (4/dom.dx^2 + 4/dom.dy^2)^2;
137 X = dom.X(dom.inside); Y = dom.Y(dom.inside);
138 out.nullx = max(abs(K * (X - 1))) * dom.dx^4; out.nully = max(abs(K * (Y - 1))) * dom.dx^4;
139 out.lnull = max(abs(L * ones(size(X))));
140 `);
141 expect(out.zero3).toBeLessThan(1e-9);
142 expect(out.fourth).toBeGreaterThan(1e-6);
143 expect(out.nullx).toBeLessThan(1e-9);
144 expect(out.nully).toBeLessThan(1e-9);
145 expect(out.lnull).toBeLessThan(1e-9);
146 // the (2, 0) pair is split by the staircase's anisotropy; its mean is the disk's value
147 expect(Math.abs((out.f1 + out.f2) / 2 - 5.358) / 5.358).toBeLessThan(0.03);
148 expect(Math.abs(out.f3 - 9.003) / 9.003).toBeLessThan(0.03);
149 expect(out.ratio).toBeLessThan(1);
150 });
151});
153describe("solving on an ellipse", () => {
154 it.each(["simply-supported", "clamped", "free", "pinned-center"] as const)(
155 "runs both finite-difference schemes with %s edges and they agree",
156 async (bc) => {
157 const spec = ellipseSpec();
158 spec.bc = bc;
159 spec.implicit = {
160 nx: spec.fd.nx,
161 ny: spec.fd.ny,
162 duration: spec.fd.duration,
163 steps_per_sample: 2,
164 theta: 0.25,
165 };
166 const exp = await solve(spec, "fd");
167 const imp = await solve(spec, "implicit");
168 for (const { result } of [exp, imp]) {
169 expect(result.info.stable).toBe(true);
170 expect(result.info.extra.shape).toBe("ellipse");
171 expect(result.probe.w.every(Number.isFinite)).toBe(true);
172 expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
173 expect(result.info.notes.some((n) => /staircase/.test(n))).toBe(true);
174 // outside the ellipse the fields are NaN, and the same nodes in every frame and field
175 const { nx, ny, nf, data } = result.frames.w;
176 const nanPer = Array.from(
177 { length: nf },
178 (_, k) =>
179 data.subarray(k * nx * ny, (k + 1) * nx * ny).filter((v) => Number.isNaN(v)).length,
180 );
181 expect(nanPer[0]).toBeGreaterThan(0.1 * nx * ny);
182 expect(nanPer[0]).toBeLessThan(0.3 * nx * ny);
183 expect(new Set(nanPer).size).toBe(1);
184 expect(result.frames.vm.data.filter((v) => Number.isNaN(v)).length).toBe(nanPer[0] * nf);
185 expect(result.envelope.vm.data.filter((v) => Number.isNaN(v)).length).toBe(nanPer[0]);
186 // the corners of the box are off the plate, the centre is on it
187 expect(data[0]).toBeNaN();
188 expect(data[ny - 1 + (nx - 1) * ny]).toBeNaN();
189 expect(Number.isFinite(data[Math.floor(ny / 2) + Math.floor(nx / 2) * ny])).toBe(true);
190 // the probe was kept where it was asked for (0.2, 0.12), on the plate
191 expect(Math.abs(result.probe.x - spec.probe.x)).toBeLessThan(
192 0.5 * (spec.plate.Lx / (nx - 1)),
193 );
194 expect(Math.abs(result.probe.y - spec.probe.y)).toBeLessThan(
195 0.5 * (spec.plate.Ly / (ny - 1)),
196 );
197 }
198 expect(Number(imp.result.info.extra.nunknowns)).toBeLessThan(spec.fd.nx * spec.fd.ny);
199 expect(Number(exp.result.info.extra.nunknowns)).toBe(Number(imp.result.info.extra.nunknowns));
200 const cmp = await runner.run({
201 bundle: buildCompareBundle(),
202 inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
203 });
204 expect(cmp.ok).toBe(true);
205 if (!cmp.ok) return;
206 const c = parseCompareResult(cmp.value);
207 expect(c.w_rel_l2).toBeLessThan(0.03);
208 expect(c.vm_rel_l2).toBeLessThan(0.1);
209 expect(c.frames_rel_l2.length).toBeGreaterThan(0);
210 expect(c.frames_rel_l2.every(Number.isFinite)).toBe(true);
211 },
212 );
214 it("is stiffer than the rectangle of the same box, and diverges above the CFL limit like it", async () => {
215 // A clamped plate on a smaller domain has a higher fundamental (domain monotonicity of the
216 // Dirichlet problem), so the ellipse inscribed in the box must beat the box itself.
217 const out = await script(`
218 r = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), 25, 17);
219 e = plate_domain(struct('shape', 'ellipse', 'Lx', 0.3, 'Ly', 0.2), 25, 17);
220 [~, Kr] = plate_matrices(r, 'clamped', 0.3); [~, Ke] = plate_matrices(e, 'clamped', 0.3);
221 out.rect = min(eig(full(Kr))); out.ell = min(eig(full(Ke)));
222 `);
223 expect(out.ell).toBeGreaterThan(out.rect);
224 // The Rayleigh quotient of the trial function (1 - x^2/a^2 - y^2/b^2)^2 is an upper bound,
225 // omega^2 rho h / D = (40/3) (3/a^4 + 2/(a^2 b^2) + 3/b^4): 10.33 against 10.216 for the disk,
226 // 17.39 for a/b = 1.5. The true value lies about a percent below it, and the staircase on a
227 // 25-node grid a further percent or two below that (the disk gives 10.21 at this resolution).
228 const a = 0.15;
229 const b = 0.1;
230 const bound = Math.sqrt((40 / 3) * (3 / a ** 4 + 2 / (a * b) ** 2 + 3 / b ** 4)) * a ** 2;
231 expect(bound).toBeCloseTo(17.39, 1);
232 const lam = Math.sqrt(out.ell) * a ** 2;
233 expect(lam).toBeLessThan(bound);
234 expect(lam).toBeGreaterThan(0.94 * bound);
235 const ell = ellipseSpec();
236 ell.fd.cfl = 1.04;
237 const blown = await solve(ell, "fd");
238 expect(blown.result.info.stable).toBe(false);
239 expect(blown.result.probe.w.every(Number.isFinite)).toBe(true);
240 });
242 it("is refused by the modal method, in the registry and in the solver", async () => {
243 const spec = ellipseSpec();
244 expect(methodBlockers(getMethod("modal"), spec).some((b) => /ellipse/.test(b))).toBe(true);
245 expect(methodBlockers(getMethod("fd"), spec)).toEqual([]);
246 expect(methodBlockers(getMethod("implicit"), spec)).toEqual([]);
247 await expect(solve(spec, "modal")).rejects.toThrow(/rectangle/);
248 });
250 it("puts the load and the probe on the plate even when the strike is near the edge", async () => {
251 const spec = ellipseSpec();
252 // just inside the ellipse near its right end, footprint spilling over the edge
253 spec.excitation.x = 0.28;
254 spec.excitation.y = 0.1;
255 spec.probe.x = 0.27;
256 spec.probe.y = 0.1;
257 const { result } = await solve(spec, "fd");
258 expect(result.info.stable).toBe(true);
259 expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
260 expect(Number.isFinite(result.probe.x)).toBe(true);
261 });
262});