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);
10 [L, B] = plate_matrices(nx, ny, dx, dy, bc{1});
11 w = zeros(ny, nx); w(2:end-1, 2:end-1) = reshape(sin(1:(ny-2)*(nx-2)), ny-2, nx-2);
12 Lw = plate_laplacian(plate_pad(w, bc{1}), dx, dy);
13 Bw = zeros(ny, nx);
14 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 ...
15 + (Lw(2:end-1,3:end) - 2*Lw(2:end-1,2:end-1) + Lw(2:end-1,1:end-2))/dx^2;
16 wi = w(2:end-1,2:end-1); Bi = Bw(2:end-1,2:end-1);
17 err.(strrep(bc{1}, '-', '_')) = norm(B*wi(:) - Bi(:)) / norm(Bi(:));
18 sp.(strrep(bc{1}, '-', '_')) = issparse(B) && issparse(L);
19 end
20 out = struct('err', err, 'sp', sp);
21 `;
22 const run = await runner.run({ bundle: { files: [{ name: "c.m", source: check }, ...LIBRARY_FILES], entry: "c.m", resultVariable: "out" } });
23 expect(run.ok).toBe(true);
24 if (!run.ok || !isPlainStruct(run.value)) return;
25 const err = run.value.fields.err;
26 const sp = run.value.fields.sp;
27 if (!isPlainStruct(err) || !isPlainStruct(sp)) throw new Error("unexpected shape");
28 expect(err.fields.simply_supported as number).toBeLessThan(1e-12);
29 expect(err.fields.clamped as number).toBeLessThan(1e-12);
30 expect(sp.fields.simply_supported).toBe(true);
31 expect(sp.fields.clamped).toBe(true);
32 });
34 it("agrees with the explicit scheme on the same grid and with the modal solution", async () => {
35 const spec = testSpec();
36 spec.implicit = { nx: spec.fd.nx, ny: spec.fd.ny, duration: spec.fd.duration, steps_per_sample: 2, theta: 0.25 };
37 spec.output.duration = spec.fd.duration;
38 spec.output.grid = [spec.fd.nx, spec.fd.ny];
39 const imp = await solve(spec, "implicit");
40 const exp = await solve(spec, "fd");
41 const modal = await solve(spec, "modal");
42 expect(imp.result.method).toBe("implicit");
43 expect(imp.result.info.stable).toBe(true);
44 expect(imp.result.info.nsteps).toBeLessThan(exp.result.info.nsteps);
45 const cmpE = await runner.run({ bundle: buildCompareBundle(), inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 } });
46 const cmpM = await runner.run({ bundle: buildCompareBundle(), inputs: { A: modal.run.ok ? modal.run.value : 0, B: imp.run.ok ? imp.run.value : 0 } });
47 expect(cmpE.ok && cmpM.ok).toBe(true);
48 if (!cmpE.ok || !cmpM.ok) return;
49 const e = parseCompareResult(cmpE.value);
50 const m = parseCompareResult(cmpM.value);
51 // same spatial operator, different time discretization: the two FD schemes agree closely
52 expect(e.w_rel_l2).toBeLessThan(0.02);
53 expect(m.w_rel_l2).toBeLessThan(0.05);
54 expect(m.vm_rel_l2).toBeLessThan(0.15);
55 expect(m.frames_rel_l2.length).toBeGreaterThan(0);
56 });
58 it("is unconditionally stable for theta = 1/4 and blows up for theta = 0 at the same step", async () => {
59 const spec = testSpec();
60 spec.implicit = { nx: 31, ny: 21, duration: 0.01, steps_per_sample: 0.25, theta: 0.25 };
61 const coarse = await solve(spec, "implicit");
62 expect(coarse.result.info.stable).toBe(true);
63 expect(maxAbs(coarse.result.probe.w)).toBeLessThan(1e-3);
64 spec.implicit = { nx: 31, ny: 21, duration: 0.002, steps_per_sample: 1, theta: 0 };
65 const explicitLike = await solve(spec, "implicit");
66 expect(explicitLike.result.info.stable).toBe(false);
67 expect(explicitLike.result.info.diverged_step).toBeGreaterThan(0);
68 expect(explicitLike.result.probe.w.every(Number.isFinite)).toBe(true);
69 });
71 it("supports clamped edges", async () => {
72 const spec = testSpec();
73 spec.bc = "clamped";
74 spec.implicit = { nx: 31, ny: 21, duration: 0.003, steps_per_sample: 1, theta: 0.25 };
75 const r = await solve(spec, "implicit");
76 expect(r.result.info.stable).toBe(true);
77 expect(maxAbs(r.result.probe.w)).toBeGreaterThan(0);
78 });
79});