8354482cymmodes: linear modes of a cymsim cymbal driven in the full nonlinear model, with the energy of every mode over timeJeremy Magland 1// Float64 reference time loop with exactly the structure of the GPU kernels: an element pass
2// (forces, sigma1 term, radiation accumulation), a hit pass (contact forces, release), a node
3// pass (gather, contact, lagged rigid-mode projection, leapfrog update, partial sums for the
4// next projection), and the exciter update. Relative to cymsimweb, this copy adds a prescribed
5// initial state, a distributed harmonic drive on the transverse DOFs, and snapshots of the state
6// projected onto a set of transverse modes.
7import { elementForces, strainEnergies } from "../model/forces.ts";
8import { HOLD_HEIGHT } from "../model/exciter.ts";
9import { DECIMATION_SPAN } from "../model/radiation.ts";
10import type { Setup } from "../model/setup.ts";
11import type { BackendOptions, ChunkResult, Diagnostics, SimulatorBackend } from "../model/simulator.ts";
13export class CpuSimulator implements SimulatorBackend {
14 readonly kind = "cpu" as const;
15 private n = 0; // sub-step counter
16 private readonly q: Float64Array; private readonly v: Float64Array; private readonly nodeAcc: Float64Array;
17 private readonly vwPrev: Float64Array; // projected transverse velocity of the previous step (for the radiated acceleration)
18 private readonly ef: Float64Array; private readonly fc: Float64Array;
19 private readonly zs: Float64Array; private readonly vs: Float64Array; private readonly Fk: Float64Array; private readonly F: Float64Array;
20 private readonly Fmax: Float64Array; private readonly ncont: Uint32Array;
21 private partialsPrev = new Float64Array(12); private partialsNew = new Float64Array(12); // [0,6): Rm.(M v), [6,12): Rm.(M q)
22 private readonly acc: Float64Array; // (E, DECIMATION_SPAN) decimation accumulators
23 private readonly history: Float64Array;
24 private pending: { j: number; front: number; back: number }[] = [];
25 samplesDone = 0;
26 readonly maxChunk = 8192;
27 readonly setup: Setup; readonly opts: BackendOptions;
29 constructor(setup: Setup, opts: BackendOptions = {}) {
30 this.setup = setup; this.opts = opts;
31 const N = setup.mesh.n, E = setup.mesh.e, K = setup.schedule.K, P = setup.schedule.P;
32 this.q = new Float64Array(5 * N); this.v = new Float64Array(5 * N); this.nodeAcc = new Float64Array(N); this.vwPrev = new Float64Array(N);
33 this.ef = new Float64Array(15 * E); this.fc = new Float64Array(N);
34 this.zs = new Float64Array(K).fill(HOLD_HEIGHT); this.vs = new Float64Array(K);
35 this.Fk = new Float64Array(K * P); this.F = new Float64Array(K); this.Fmax = new Float64Array(K); this.ncont = new Uint32Array(K);
36 this.acc = new Float64Array(E * DECIMATION_SPAN); this.history = new Float64Array(E * setup.radiation.ringLength);
37 }
39 setState(q: Float64Array, v: Float64Array) {
40 this.q.set(q); this.v.set(v);
41 const N = this.setup.mesh.n;
42 for (let nd = 0; nd < N; nd++) this.vwPrev[nd] = v[5 * nd + 2];
43 }
45 private step() {
46 const { shell: sh, mesh, schedule: sc, radiation: rad, dt, substeps: m } = this.setup;
47 const { tri } = mesh, N = mesh.n, E = mesh.e, n = this.n;
48 const sigma0 = this.setup.config.material.sigma0;
49 // ---- element pass
50 elementForces(sh, this.q, this.v, this.ef);
51 if (n >= 1) {
52 const s = n - 1, j0 = Math.floor(s / m), kk = s - j0 * m, W = DECIMATION_SPAN, taps = rad.taps;
53 if (kk === 0 && j0 > 0) this.finalizeSample(j0 - W / 2); // slot 0 = sample j0 - (W/2 - 1) of the previous j0
54 for (let e = 0; e < E; e++) {
55 const a = (this.nodeAcc[tri[3 * e]] + this.nodeAcc[tri[3 * e + 1]] + this.nodeAcc[tri[3 * e + 2]]) / 3;
56 for (let w = 0; w < W; w++) this.acc[e * W + w] += a * taps[w * m + kk];
57 }
58 }
59 // ---- hit pass
60 this.fc.fill(0);
61 for (let k = 0; k < sc.K; k++) {
62 let wmax = -Infinity;
63 const wp = new Float64Array(sc.P);
64 for (let p = 0; p < sc.P; p++) {
65 const b = (k * sc.P + p) * 3;
66 let w = sc.offsets[k * sc.P + p];
67 for (let i = 0; i < 3; i++) w += sc.weights[b + i] * this.q[5 * sc.nodes[b + i] + 2];
68 wp[p] = w; wmax = Math.max(wmax, w);
69 }
70 if (n === sc.releaseStep[k]) { this.zs[k] = wmax + sc.gap[k]; this.vs[k] = -sc.velocity[k]; }
71 const active = n >= sc.releaseStep[k];
72 let Ftot = 0;
73 for (let p = 0; p < sc.P; p++) {
74 const delta = Math.max(wp[p] - this.zs[k], 0);
75 const f = active ? sc.kPoint[k] * Math.pow(delta, sc.alpha[k]) : 0;
76 this.Fk[k * sc.P + p] = f; Ftot += f;
77 if (f > 0) { const b = (k * sc.P + p) * 3; for (let i = 0; i < 3; i++) this.fc[sc.nodes[b + i]] -= f * sc.weights[b + i]; }
78 }
79 this.F[k] = Ftot;
80 }
81 // ---- drive: distributed transverse force fNode * env * sin(omega t) added with the contacts
82 const drv = this.opts.drive;
83 if (drv && n < drv.stopStep) {
84 const env = Math.min(n / drv.rampSteps, 1) * Math.sin(drv.omega * dt * n);
85 for (let nd = 0; nd < N; nd++) this.fc[nd] += drv.fNode[nd] * env;
86 }
87 // ---- node pass
88 const { M, Minv, Rm, adjOffsets, adjEntries } = sh, L = 5 * N, c = this.partialsPrev;
89 this.partialsNew.fill(0);
90 const f = new Float64Array(5);
91 for (let nd = 0; nd < N; nd++) {
92 f.fill(0);
93 for (let a = adjOffsets[nd]; a < adjOffsets[nd + 1]; a++) {
94 const ent = adjEntries[a], e = (ent / 3) | 0, i = ent % 3;
95 for (let d = 0; d < 5; d++) f[d] += this.ef[15 * e + 5 * i + d];
96 }
97 f[2] += this.fc[nd];
98 for (let d = 0; d < 5; d++) {
99 const idx = 5 * nd + d;
100 // lagged projection: remove the rigid content that the previous step's update introduced
101 let vd = this.v[idx], qd = this.q[idx];
102 for (let k = 0; k < 6; k++) { const r = Rm[k * L + idx]; vd -= c[k] * r; qd -= c[6 + k] * r; }
103 const v1 = (vd + dt * f[d] * Minv[idx]) / (1 + 2 * sigma0 * dt);
104 // radiated acceleration from consecutive projected velocities, so that the rigid increment
105 // that the lagged projection removes one step later never radiates
106 if (d === 2) { this.nodeAcc[nd] = (vd - this.vwPrev[nd]) / dt; this.vwPrev[nd] = vd; }
107 const q1 = qd + dt * v1;
108 this.q[idx] = q1; this.v[idx] = v1;
109 for (let k = 0; k < 6; k++) { const r = Rm[k * L + idx] * M[idx]; this.partialsNew[k] += r * v1; this.partialsNew[6 + k] += r * q1; }
110 }
111 }
112 [this.partialsPrev, this.partialsNew] = [this.partialsNew, this.partialsPrev];
113 // ---- exciters
114 for (let k = 0; k < sc.K; k++) {
115 this.vs[k] += (dt * this.F[k]) / sc.mass[k]; this.zs[k] += dt * this.vs[k];
116 this.Fmax[k] = Math.max(this.Fmax[k], this.F[k]); if (this.F[k] > 0) this.ncont[k]++;
117 }
118 this.n++;
119 }
121 /** Slot 0 (audio sample j) of the accumulators is complete: store it, shift the slots, evaluate the listener pressure. */
122 private finalizeSample(j: number) {
123 const { radiation: rad, mesh } = this.setup, E = mesh.e, D = rad.ringLength, W = DECIMATION_SPAN;
124 for (let e = 0; e < E; e++) {
125 if (j >= 0) this.history[e * D + (j % D)] = this.acc[e * W];
126 for (let w = 0; w < W - 1; w++) this.acc[e * W + w] = this.acc[e * W + w + 1];
127 this.acc[e * W + W - 1] = 0;
128 }
129 if (j < 0) return;
130 let front = 0, back = 0;
131 for (let e = 0; e < E; e++) {
132 const jf = j - rad.i0F[e], jb = j - rad.i0B[e];
133 if (jf >= 1) front += rad.coefF[e] * (this.history[e * D + (jf % D)] * (1 - rad.frF[e]) + this.history[e * D + ((jf - 1) % D)] * rad.frF[e]);
134 if (jb >= 1) back += rad.coefB[e] * (this.history[e * D + (jb % D)] * (1 - rad.frB[e]) + this.history[e * D + ((jb - 1) % D)] * rad.frB[e]);
135 }
136 this.pending.push({ j, front, back });
137 }
139 /** Project the current state onto the modes: 2 Nm values, q_m then qdot_m. */
140 private projectNow(): Float32Array {
141 const pr = this.opts.projection!;
142 const N = this.setup.mesh.n, Nm = pr.nModes, out = new Float32Array(2 * Nm);
143 for (let m = 0; m < Nm; m++) {
144 let sq = 0, sv = 0;
145 const row = m * N;
146 for (let nd = 0; nd < N; nd++) { const p = pr.phiM[row + nd]; sq += p * this.q[5 * nd + 2]; sv += p * this.v[5 * nd + 2]; }
147 out[m] = sq; out[Nm + m] = sv;
148 }
149 return out;
150 }
152 async runChunk(nSamples: number): Promise<ChunkResult> {
153 const m = this.setup.substeps, N = this.setup.mesh.n, frames: Float32Array[] = [], frameSamples: number[] = [];
154 const proj: Float32Array[] = [], projSamples: number[] = [];
155 const movieEvery = this.opts.movieEvery ?? 0, pr = this.opts.projection;
156 const tm = this.opts.traceMode, tracing = pr && tm !== undefined;
157 const traceBase = this.samplesDone;
158 const trace = new Float32Array(tracing ? nSamples : 0);
159 for (let s = 0; s < nSamples * m; s++) {
160 if (movieEvery > 0 && this.n % (m * movieEvery) === 0) {
161 const w = new Float32Array(N);
162 for (let nd = 0; nd < N; nd++) w[nd] = this.q[5 * nd + 2];
163 frames.push(w); frameSamples.push(this.n / m);
164 }
165 if (pr && this.n % (m * pr.every) === 0) { proj.push(this.projectNow()); projSamples.push(this.n / m); }
166 if (tracing && this.n % m === 0) {
167 let sq = 0;
168 const row = tm! * N;
169 for (let nd = 0; nd < N; nd++) sq += pr!.phiM[row + nd] * this.q[5 * nd + 2];
170 trace[this.n / m - traceBase] = sq;
171 }
172 this.step();
173 }
174 this.samplesDone += nSamples;
175 const out = this.pending; this.pending = [];
176 const firstSample = out.length ? out[0].j : this.samplesDone;
177 return { firstSample, pressureFront: Float32Array.from(out, (o) => o.front), pressureBack: Float32Array.from(out, (o) => o.back),
178 frames, frameSamples, proj, projSamples, trace, traceBase };
179 }
181 async diagnostics(): Promise<Diagnostics> {
182 const { shell: sh, schedule: sc, dt } = this.setup, N = sh.mesh.n;
183 let ke = 0, woh = 0;
184 for (let i = 0; i < 5 * N; i++) ke += 0.5 * sh.M[i] * this.v[i] * this.v[i];
185 for (let nd = 0; nd < N; nd++) woh = Math.max(woh, Math.abs(this.q[5 * nd + 2] / sh.hNode[nd]));
186 const { membrane, bending } = strainEnergies(sh, this.q);
187 let ex = 0;
188 for (let k = 0; k < sc.K; k++) ex += 0.5 * sc.mass[k] * this.vs[k] * this.vs[k];
189 return { time: this.n * dt, kinetic: ke, membrane, bending, exciters: ex, maxWoverH: woh,
190 hits: Array.from({ length: sc.K }, (_, k) => ({ peakForce: this.Fmax[k], contactTime: this.ncont[k] * dt, velocity: this.vs[k], height: this.zs[k] })) };
191 }
193 /** M-weighted rigid-mode content of the current velocity (the lagged projection leaves one step's increment). */
194 rigidContent(): number {
195 const { shell: sh } = this.setup, L = 5 * sh.mesh.n;
196 let worst = 0;
197 for (let k = 0; k < 6; k++) { let c = 0; for (let i = 0; i < L; i++) c += sh.Rm[k * L + i] * sh.M[i] * this.v[i]; worst = Math.max(worst, Math.abs(c)); }
198 return worst;
199 }
201 async displacement(): Promise<Float32Array> {
202 const N = this.setup.mesh.n, w = new Float32Array(N);
203 for (let nd = 0; nd < N; nd++) w[nd] = this.q[5 * nd + 2];
204 return w;
205 }
207 dispose() {}
208}