1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 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";
7test("CPU time loop reproduces the Python run (dent-free, float64)", async () => {
8 const setup = await buildSetup(fixture.config);
9 assert.equal(setup.substeps, fixture.timestep.substeps);
10 assert.ok(Math.abs(setup.dt / fixture.timestep.dt - 1) < 1e-12);
11 const sim = new CpuSimulator(setup);
12 const front: number[] = [], back: number[] = [];
13 const py = fixture.run, hit = py.hits[0], last = py.energy_log[py.energy_log.length - 1];
14 const nSamples = Math.round(last[0] * setup.fs); // Python ran whole blocks, slightly past the duration
15 const chunk = 147;
16 for (let done = 0; done < nSamples; done += chunk) {
17 const r = await sim.runChunk(Math.min(chunk, nSamples - done));
18 front.push(...r.pressureFront); back.push(...r.pressureBack);
19 }
20 const d = await sim.diagnostics();
21 assert.ok(Math.abs(d.time - last[0]) < 1e-9);
22 const rel = (a: number, b: number) => Math.abs(a - b) / Math.abs(b);
23 assert.ok(rel(d.hits[0].peakForce, hit.peak_force) < 1e-3, `peak ${d.hits[0].peakForce} vs ${hit.peak_force}`);
24 assert.ok(Math.abs(d.hits[0].contactTime - hit.contact_time) <= 2 * setup.dt + 1e-12);
25 assert.ok(rel(d.hits[0].velocity, hit.rebound_velocity) < 1e-3, `rebound ${d.hits[0].velocity} vs ${hit.rebound_velocity}`);
26 // the run is strongly nonlinear (max|w|/h ~ 10), so small scheme differences are amplified: 1 percent
27 assert.ok(rel(d.kinetic, last[1]) < 1e-2 && rel(d.membrane, last[2]) < 1e-2 && rel(d.bending, last[3]) < 1e-2,
28 `energies ${[d.kinetic, d.membrane, d.bending]} vs ${last.slice(1, 4)}`);
29 assert.ok(rel(d.exciters, last[4]) < 1e-2 && rel(d.maxWoverH, last[5]) < 1e-2);
30 // the lagged projection leaves a rigid velocity of one step's increment: tiny against the vibration
31 assert.ok(sim.rigidContent() < 1e-3 * Math.sqrt(2 * d.kinetic), `rigid content ${sim.rigidContent()}`);
32 // radiation: a different decimation and delay scheme, so compare loosely (rms and correlation)
33 const pu = front.map((f, i) => f - back[i]), ref: number[] = py.pressure_unbaffled;
34 const n = Math.min(pu.length, ref.length) - 4;
35 const rms = (a: number[]) => Math.sqrt(a.slice(0, n).reduce((s, x) => s + x * x, 0) / n);
36 let best = 0;
37 for (let lag = -3; lag <= 3; lag++) {
38 let dot = 0;
39 for (let i = 4; i < n; i++) dot += pu[i] * ref[i + lag];
40 best = Math.max(best, dot / (rms(pu) * rms(ref) * (n - 4)));
41 }
42 assert.ok(Math.abs(rms(pu) / rms(ref) - 1) < 0.1, `rms ${rms(pu)} vs ${rms(ref)}`);
43 assert.ok(best > 0.95, `correlation ${best}`);
44});
46test("linear run matches Python tightly (no chaotic amplification)", async () => {
47 const cfg = structuredClone(fixture.config);
48 cfg.numerics.nonlinear = false;
49 const setup = await buildSetup(cfg);
50 const sim = new CpuSimulator(setup);
51 const py = fixture.run_linear, hit = py.hits[0], last = py.energy_log[py.energy_log.length - 1];
52 const nSamples = Math.round(last[0] * setup.fs);
53 const front: number[] = [], back: number[] = [];
54 for (let done = 0; done < nSamples; done += 200) {
55 const r = await sim.runChunk(Math.min(200, nSamples - done));
56 front.push(...r.pressureFront); back.push(...r.pressureBack);
57 }
58 const d = await sim.diagnostics();
59 const rel = (a: number, b: number) => Math.abs(a - b) / Math.abs(b);
60 assert.ok(rel(d.hits[0].peakForce, hit.peak_force) < 1e-5, `peak ${d.hits[0].peakForce} vs ${hit.peak_force}`);
61 assert.ok(rel(d.hits[0].velocity, hit.rebound_velocity) < 1e-5);
62 assert.ok(rel(d.kinetic, last[1]) < 1e-4 && rel(d.membrane, last[2]) < 1e-4 && rel(d.bending, last[3]) < 1e-4,
63 `energies ${[d.kinetic, d.membrane, d.bending]} vs ${last.slice(1, 4)}`);
64 assert.ok(rel(d.maxWoverH, last[5]) < 1e-4);
65 const pu = front.map((f, i) => f - back[i]), ref: number[] = py.pressure_unbaffled;
66 const n = Math.min(pu.length, ref.length) - 4;
67 const rms = (a: number[]) => Math.sqrt(a.slice(0, n).reduce((s, x) => s + x * x, 0) / n);
68 assert.ok(Math.abs(rms(pu) / rms(ref) - 1) < 0.1, `rms ${rms(pu)} vs ${rms(ref)}`);
69});