2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 1import { describe, expect, it } from "vitest";
2import { LIBRARY_FILES, buildCompareBundle, isPlainStruct, parseCompareResult } from "../src/index";
3import { maxAbs, runner, solve, testSpec } from "./helpers";
5describe("implicit finite differences", () => {
6 it("assembles the same operator as the explicit ghost-node stencils, for both edge conditions", async () => {
7 const check = `
8 for bc = {'simply-supported', 'clamped'}
9 nx = 19; ny = 13; dx = 0.3/(nx-1); dy = 0.2/(ny-1);
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 10 dom = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), nx, ny);
11 [L, B] = plate_matrices(dom, bc{1});
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 12 w = zeros(ny, nx); w(2:end-1, 2:end-1) = reshape(sin(1:(ny-2)*(nx-2)), ny-2, nx-2);
13 Lw = plate_laplacian(plate_pad(w, bc{1}), dx, dy);
14 Bw = zeros(ny, nx);
15 Bw(2:end-1,2:end-1) = (Lw(3:end,2:end-1) - 2*Lw(2:end-1,2:end-1) + Lw(1:end-2,2:end-1))/dy^2 ...
16 + (Lw(2:end-1,3:end) - 2*Lw(2:end-1,2:end-1) + Lw(2:end-1,1:end-2))/dx^2;
17 wi = w(2:end-1,2:end-1); Bi = Bw(2:end-1,2:end-1);
18 err.(strrep(bc{1}, '-', '_')) = norm(B*wi(:) - Bi(:)) / norm(Bi(:));
19 sp.(strrep(bc{1}, '-', '_')) = issparse(B) && issparse(L);
20 end
21 out = struct('err', err, 'sp', sp);
22 `;
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 23 const run = await runner.run({
24 bundle: {
25 files: [{ name: "c.m", source: check }, ...LIBRARY_FILES],
26 entry: "c.m",
27 resultVariable: "out",
28 },
29 });
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 30 expect(run.ok).toBe(true);
31 if (!run.ok || !isPlainStruct(run.value)) return;
32 const err = run.value.fields.err;
33 const sp = run.value.fields.sp;
34 if (!isPlainStruct(err) || !isPlainStruct(sp)) throw new Error("unexpected shape");
35 expect(err.fields.simply_supported as number).toBeLessThan(1e-12);
36 expect(err.fields.clamped as number).toBeLessThan(1e-12);
37 expect(sp.fields.simply_supported).toBe(true);
38 expect(sp.fields.clamped).toBe(true);
39 });
41 it("agrees with the explicit scheme on the same grid and with the modal solution", async () => {
42 const spec = testSpec();
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 43 spec.implicit = {
44 nx: spec.fd.nx,
45 ny: spec.fd.ny,
46 duration: spec.fd.duration,
47 steps_per_sample: 2,
48 theta: 0.25,
49 };
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 50 spec.output.duration = spec.fd.duration;
51 spec.output.grid = [spec.fd.nx, spec.fd.ny];
52 const imp = await solve(spec, "implicit");
53 const exp = await solve(spec, "fd");
54 const modal = await solve(spec, "modal");
55 expect(imp.result.method).toBe("implicit");
56 expect(imp.result.info.stable).toBe(true);
57 expect(imp.result.info.nsteps).toBeLessThan(exp.result.info.nsteps);
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 58 const cmpE = await runner.run({
59 bundle: buildCompareBundle(),
60 inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
61 });
62 const cmpM = await runner.run({
63 bundle: buildCompareBundle(),
64 inputs: { A: modal.run.ok ? modal.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
65 });
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 66 expect(cmpE.ok && cmpM.ok).toBe(true);
67 if (!cmpE.ok || !cmpM.ok) return;
68 const e = parseCompareResult(cmpE.value);
69 const m = parseCompareResult(cmpM.value);
70 // same spatial operator, different time discretization: the two FD schemes agree closely
71 expect(e.w_rel_l2).toBeLessThan(0.02);
72 expect(m.w_rel_l2).toBeLessThan(0.05);
73 expect(m.vm_rel_l2).toBeLessThan(0.15);
74 expect(m.frames_rel_l2.length).toBeGreaterThan(0);
75 });
77 it("is unconditionally stable for theta = 1/4 and blows up for theta = 0 at the same step", async () => {
78 const spec = testSpec();
79 spec.implicit = { nx: 31, ny: 21, duration: 0.01, steps_per_sample: 0.25, theta: 0.25 };
80 const coarse = await solve(spec, "implicit");
81 expect(coarse.result.info.stable).toBe(true);
82 expect(maxAbs(coarse.result.probe.w)).toBeLessThan(1e-3);
83 spec.implicit = { nx: 31, ny: 21, duration: 0.002, steps_per_sample: 1, theta: 0 };
84 const explicitLike = await solve(spec, "implicit");
85 expect(explicitLike.result.info.stable).toBe(false);
86 expect(explicitLike.result.info.diverged_step).toBeGreaterThan(0);
87 expect(explicitLike.result.probe.w.every(Number.isFinite)).toBe(true);
88 });
90 it("supports clamped edges", async () => {
91 const spec = testSpec();
92 spec.bc = "clamped";
93 spec.implicit = { nx: 31, ny: 21, duration: 0.003, steps_per_sample: 1, theta: 0.25 };
94 const r = await solve(spec, "implicit");
95 expect(r.result.info.stable).toBe(true);
96 expect(maxAbs(r.result.probe.w)).toBeGreaterThan(0);
97 });
98});