1// The modal drive and projection against the float64 CPU loop, on the small test cymbal.
2//
3// The linear runs are the calibration: released from a mode shape (with the condensation
4// complement) the linearised model must keep the energy in that mode and decay it at exp(-2
5// sigma0 t); driven at resonance with the mass-weighted shape it must respond in that mode only.
6// The linear basis and the simulator are not identical (the simulator carries in-plane and
7// rotational inertia that the condensed eigenproblem drops), so "keeps the energy" is a
8// tolerance, not an identity, and the nonlinear run is compared against the linear one rather
9// than against zero.
10import assert from "node:assert/strict";
11import { test } from "node:test";
12import { cloneConfig } from "../src/config/yaml.ts";
13import { CpuSimulator } from "../src/cpu/simulator.ts";
14import { buildDrive, buildProjection, modalAmplitude, modalInitialState, modeShape, type DriveSpec } from "../src/drive.ts";
15import { ModalTracker } from "../src/modal/tracker.ts";
16import { internalForce } from "../src/model/forces.ts";
17import { buildSetup } from "../src/model/setup.ts";
18import type { BackendOptions } from "../src/model/simulator.ts";
19import { buildCondensationOperator } from "../src/precompute/assemble.ts";
20import { computeModes } from "../src/precompute/pipeline.ts";
21import { miniConfig, noProgress, OPTIONS } from "./helpers.ts";
23const base = miniConfig();
24const modes = computeModes(base, OPTIONS, noProgress);
25const MODE = 3;
27function transverseMass(setup: Awaited<ReturnType<typeof buildSetup>>): Float64Array {
28 const N = setup.mesh.n, Mw = new Float64Array(N);
29 for (let n = 0; n < N; n++) Mw[n] = setup.shell.M[5 * n + 2];
30 return Mw;
31}
33interface RunResult { tracker: ModalTracker; freq: Float64Array }
35async function runDecay(amplitude: number, nonlinear: boolean, o: Float64Array): Promise<RunResult> {
36 const cfg = cloneConfig(base);
37 cfg.numerics.nonlinear = nonlinear;
38 const setup = await buildSetup(cfg);
39 assert.equal(setup.mesh.n, modes.nNodes);
40 const shape = modeShape(modes, MODE, 1);
41 const initial = modalInitialState(shape, o, modalAmplitude(modes, shape, amplitude));
42 const Mw = transverseMass(setup);
43 const opts: BackendOptions = { projection: buildProjection(modes, Mw, 1, Math.max(1, Math.round(setup.fs / 500))) };
44 const sim = new CpuSimulator(setup, opts);
45 sim.setState(initial.q, initial.v);
46 const tracker = new ModalTracker(modes.freq, MODE);
47 let left = setup.nSamples;
48 while (left > 0) { const n = Math.min(2048, left); tracker.addChunk(await sim.runChunk(n), setup.fs); left -= n; }
49 return { tracker, freq: modes.freq };
50}
52function drivenFraction(t: ModalTracker, s: number): number {
53 const total = t.totals[s];
54 return total > 0 ? t.energies[s][t.drivenMode] / total : 1;
55}
57const baseSetup = await buildSetup(base);
58const cond = buildCondensationOperator(baseSetup.mesh, baseSetup.shell, noProgress);
59const complement = cond.complement(modeShape(modes, MODE, 1));
61test("the complement makes the mode an equilibrium of the in-plane and rotation DOFs", async () => {
62 const cfg = cloneConfig(base);
63 cfg.numerics.nonlinear = false;
64 const setup = await buildSetup(cfg);
65 const shape = modeShape(modes, MODE, 1);
66 const a = modalAmplitude(modes, shape, 0.05);
67 const { q } = modalInitialState(shape, complement, a);
68 const f = internalForce(setup.shell, q);
69 const N = setup.mesh.n;
70 let fw = 0, fo = 0;
71 for (let n = 0; n < N; n++) {
72 fw += f[5 * n + 2] ** 2;
73 fo += f[5 * n] ** 2 + f[5 * n + 1] ** 2 + f[5 * n + 3] ** 2 + f[5 * n + 4] ** 2;
74 }
75 assert.ok(Math.sqrt(fo / fw) < 0.05, `off-transverse force ${Math.sqrt(fo / fw)} of the transverse force`);
76 // and the transverse force is -omega^2 Mw w of the eigenproblem
77 const omega = 2 * Math.PI * modes.freq[MODE];
78 let num = 0, den = 0;
79 for (let n = 0; n < N; n++) {
80 const expected = -omega * omega * setup.shell.M[5 * n + 2] * q[5 * n + 2];
81 num += (f[5 * n + 2] - expected) ** 2; den += expected ** 2;
82 }
83 assert.ok(Math.sqrt(num / den) < 0.05, `transverse force off the eigenvalue by ${Math.sqrt(num / den)}`);
84});
86test("projection of the initial state is the driven mode alone", async () => {
87 const cfg = cloneConfig(base);
88 const setup = await buildSetup(cfg);
89 const shape = modeShape(modes, MODE, 1);
90 const a = modalAmplitude(modes, shape, 0.05);
91 const { q } = modalInitialState(shape, complement, a);
92 const proj = buildProjection(modes, transverseMass(setup), 1, 1);
93 const N = setup.mesh.n;
94 let worst = 0;
95 for (let m = 0; m < modes.nModes; m++) {
96 let s = 0;
97 for (let n = 0; n < N; n++) s += proj.phiM[m * N + n] * q[5 * n + 2];
98 if (m === MODE) assert.ok(Math.abs(s / a - 1) < 0.01, `driven coordinate ${s} vs ${a}`);
99 else worst = Math.max(worst, Math.abs(s / a));
100 }
101 assert.ok(worst < 0.02, `foreign coordinate up to ${worst} of the driven one`);
102});
104const linear = await runDecay(0.05, false, complement);
105const linearEnd = linear.tracker.length - 1;
107test("a small-amplitude linear release stays in its mode and decays at exp(-2 sigma0 t)", () => {
108 const t = linear.tracker;
109 assert.ok(t.length > 20);
110 assert.ok(drivenFraction(t, 0) > 0.99, `initial fraction ${drivenFraction(t, 0)}`);
111 assert.ok(drivenFraction(t, linearEnd) > 0.85, `final fraction ${drivenFraction(t, linearEnd)}`);
112 // total modal energy against the exact initial value and the expected decay
113 const shape = modeShape(modes, MODE, 1);
114 const a = modalAmplitude(modes, shape, 0.05);
115 const omega = 2 * Math.PI * modes.freq[MODE];
116 const E0 = 0.5 * omega * omega * a * a;
117 assert.ok(Math.abs(t.totals[0] / E0 - 1) < 0.02, `E(0) = ${t.totals[0]} vs ${E0}`);
118 const span = t.times[linearEnd] - t.times[0];
119 const expected = Math.exp(-2 * base.material.sigma0 * span);
120 const measured = t.totals[linearEnd] / t.totals[0];
121 assert.ok(Math.abs(measured / expected - 1) < 0.15, `decay ${measured} vs ${expected}`);
122});
124test("a strong nonlinear release moves energy into other modes; the linear run does not", async () => {
125 const nl = await runDecay(4, true, complement);
126 const s = Math.min(linearEnd, nl.tracker.length - 1);
127 const linOther = 1 - drivenFraction(linear.tracker, s);
128 const nlOther = 1 - drivenFraction(nl.tracker, s);
129 assert.ok(nlOther > 0.05, `nonlinear transfer only ${nlOther}`);
130 assert.ok(nlOther > 5 * linOther, `nonlinear ${nlOther} vs linear ${linOther}`);
131});
133test("harmonic forcing at resonance excites the driven mode alone (linear)", async () => {
134 const cfg = cloneConfig(base);
135 cfg.numerics.nonlinear = false;
136 cfg.numerics.duration = 0.05;
137 const setup = await buildSetup(cfg);
138 const shape = modeShape(modes, MODE, 1);
139 const spec: DriveSpec = { mode: MODE, mechanism: "forced", amplitude: 0.2, detuneCents: 0, rampTime: 0.005, driveTime: 1 };
140 const sigma = cfg.material.sigma0;
141 const drive = buildDrive(spec, modes, shape, transverseMass(setup), sigma, 2 * Math.PI * modes.freq[MODE], setup.dt);
142 const sim = new CpuSimulator(setup, { drive, traceMode: MODE,
143 projection: buildProjection(modes, transverseMass(setup), 1, Math.max(1, Math.round(setup.fs / 500))) });
144 const tracker = new ModalTracker(modes.freq, MODE);
145 let left = setup.nSamples, traceChecked = 0;
146 while (left > 0) {
147 const n = Math.min(2048, left);
148 const res = await sim.runChunk(n);
149 // the audio-rate trace agrees with the snapshot projection where they coincide
150 assert.equal(res.trace.length, n);
151 for (let s = 0; s < res.proj.length; s++) {
152 const q = res.proj[s][MODE], t = res.trace[res.projSamples[s] - res.traceBase];
153 assert.ok(Math.abs(t - q) <= 1e-6 * Math.max(Math.abs(q), 1e-12), `trace ${t} vs projection ${q}`);
154 traceChecked++;
155 }
156 tracker.addChunk(res, setup.fs);
157 left -= n;
158 }
159 assert.ok(traceChecked > 10, `only ${traceChecked} trace samples checked`);
160 const end = tracker.length - 1;
161 assert.ok(tracker.energies[end][MODE] > 0, "no response to the drive");
162 assert.ok(drivenFraction(tracker, end) > 0.95, `driven fraction ${drivenFraction(tracker, end)}`);
163 // the response grows towards the steady state (far from it after 0.05 s, but it must be rising)
164 const mid = tracker.indexAtTime(0.02);
165 assert.ok(tracker.energies[end][MODE] > 2 * tracker.energies[mid][MODE], "the driven mode is not growing");
166});