1import test from "node:test";
2import assert from "node:assert/strict";
3import { buildSetup } from "../src/model/setup.ts";
4import { CpuSimulator } from "../src/cpu/simulator.ts";
5import { fixture } from "./helpers.ts";
7/** Relative rms difference over the common length. */
8function relRms(a: ArrayLike<number>, b: ArrayLike<number>): number {
9 let num = 0, den = 0;
10 for (let i = 0; i < Math.min(a.length, b.length); i++) { num += (a[i] - b[i]) ** 2; den += b[i] ** 2; }
11 return Math.sqrt(num / den);
12}
14test("CPU time loop reproduces the Python run (dent-free, float64)", async () => {
15 const setup = await buildSetup(fixture.config);
16 assert.equal(setup.substeps, fixture.timestep.substeps);
17 assert.ok(Math.abs(setup.dt / fixture.timestep.dt - 1) < 1e-12);
18 const sim = new CpuSimulator(setup, { recordAcceleration: true });
19 const front: number[] = [], back: number[] = [], acc: Float32Array[] = [], N = setup.mesh.n;
20 const py = fixture.run, hit = py.hits[0], last = py.energy_log[py.energy_log.length - 1];
21 const nSamples = Math.round(last[0] * setup.fs); // Python ran whole blocks, slightly past the duration
22 const chunk = 147;
23 for (let done = 0; done < nSamples; done += chunk) {
24 const r = await sim.runChunk(Math.min(chunk, nSamples - done));
25 front.push(...r.pressureFront); back.push(...r.pressureBack);
26 for (let i = 0; i < r.pressureFront.length; i++) acc.push(r.acceleration!.slice(i * N, (i + 1) * N));
27 }
28 const d = await sim.diagnostics();
29 assert.ok(Math.abs(d.time - last[0]) < 1e-9);
30 const rel = (a: number, b: number) => Math.abs(a - b) / Math.abs(b);
31 assert.ok(rel(d.hits[0].peakForce, hit.peak_force) < 1e-3, `peak ${d.hits[0].peakForce} vs ${hit.peak_force}`);
32 assert.ok(Math.abs(d.hits[0].contactTime - hit.contact_time) <= 2 * setup.dt + 1e-12);
33 assert.ok(rel(d.hits[0].velocity, hit.rebound_velocity) < 1e-3, `rebound ${d.hits[0].velocity} vs ${hit.rebound_velocity}`);
34 // the run is strongly nonlinear (max|w|/h ~ 10), so small scheme differences are amplified: 1 percent
35 assert.ok(rel(d.kinetic, last[1]) < 1e-2 && rel(d.membrane, last[2]) < 1e-2 && rel(d.bending, last[3]) < 1e-2,
36 `energies ${[d.kinetic, d.membrane, d.bending]} vs ${last.slice(1, 4)}`);
37 assert.ok(rel(d.exciters, last[4]) < 1e-2 && rel(d.maxWoverH, last[5]) < 1e-2);
38 // the lagged projection leaves a rigid velocity of one step's increment: tiny against the vibration
39 assert.ok(sim.rigidContent() < 1e-3 * Math.sqrt(2 * d.kinetic), `rigid content ${sim.rigidContent()}`);
40 // radiation and the acceleration record: the same decimation kernel and normal component as Python;
41 // the delay is applied at sub-step resolution there and by a 16-tap interpolator here
42 const pu = front.map((f, i) => f - back[i]), ref: number[] = py.pressure_unbaffled;
43 assert.ok(relRms(pu, ref) < 1e-3, `unbaffled rel rms error ${relRms(pu, ref)}`);
44 assert.ok(relRms(front, py.pressure_rayleigh) < 1e-3, `baffled rel rms error ${relRms(front, py.pressure_rayleigh)}`);
45 py.acceleration_nodes.forEach((nd: number, k: number) => {
46 const mine = acc.map((row) => row[nd]);
47 assert.ok(relRms(mine, py.acceleration[k]) < 1e-3, `acceleration of node ${nd}: rel rms error ${relRms(mine, py.acceleration[k])}`);
48 });
49});
51test("linear run matches Python tightly (no chaotic amplification)", async () => {
52 const cfg = structuredClone(fixture.config);
53 cfg.numerics.nonlinear = false;
54 const setup = await buildSetup(cfg);
55 const sim = new CpuSimulator(setup);
56 const py = fixture.run_linear, hit = py.hits[0], last = py.energy_log[py.energy_log.length - 1];
57 const nSamples = Math.round(last[0] * setup.fs);
58 const front: number[] = [], back: number[] = [];
59 for (let done = 0; done < nSamples; done += 200) {
60 const r = await sim.runChunk(Math.min(200, nSamples - done));
61 front.push(...r.pressureFront); back.push(...r.pressureBack);
62 }
63 const d = await sim.diagnostics();
64 const rel = (a: number, b: number) => Math.abs(a - b) / Math.abs(b);
65 assert.ok(rel(d.hits[0].peakForce, hit.peak_force) < 1e-5, `peak ${d.hits[0].peakForce} vs ${hit.peak_force}`);
66 assert.ok(rel(d.hits[0].velocity, hit.rebound_velocity) < 1e-5);
67 assert.ok(rel(d.kinetic, last[1]) < 1e-4 && rel(d.membrane, last[2]) < 1e-4 && rel(d.bending, last[3]) < 1e-4,
68 `energies ${[d.kinetic, d.membrane, d.bending]} vs ${last.slice(1, 4)}`);
69 assert.ok(rel(d.maxWoverH, last[5]) < 1e-4);
70 const pu = front.map((f, i) => f - back[i]);
71 assert.ok(relRms(pu, py.pressure_unbaffled) < 1e-3, `unbaffled rel rms error ${relRms(pu, py.pressure_unbaffled)}`);
72});