1// Float64 reference time loop with exactly the structure of the GPU kernels: an element pass
2// (forces, sigma1 term), a hit pass (contact forces, release), a node pass (gather, contact,
3// lagged rigid-mode projection, leapfrog update, partial sums for the next projection, normal
4// acceleration and its decimation accumulators), the exciter update, and once per audio sample
5// the gather of the listener pressure from the element history. See THEORY.md.
6import { elementForces, strainEnergies } from "../model/forces.ts";
7import { HOLD_HEIGHT } from "../model/exciter.ts";
8import { DECIMATION_SPAN, DELAY_TAPS, delayedValue } from "../model/radiation.ts";
9import type { Setup } from "../model/setup.ts";
10import type { BackendOptions, BackendState, ChunkResult, Diagnostics, SimulatorBackend } from "../model/simulator.ts";
12export class CpuSimulator implements SimulatorBackend {
13 readonly kind = "cpu" as const;
14 private n = 0; // sub-step counter
15 private readonly q: Float64Array; private readonly v: Float64Array;
16 private readonly vPrev: Float64Array; // projected (u, v, w) velocities of the previous step (for the radiated acceleration)
17 private readonly ef: Float64Array; private readonly fc: Float64Array;
18 private readonly zs: Float64Array; private readonly vs: Float64Array; private readonly Fk: Float64Array; private readonly F: Float64Array;
19 private readonly Fmax: Float64Array; private readonly ncont: Uint32Array;
20 private partialsPrev = new Float64Array(12); private partialsNew = new Float64Array(12); // [0,6): Rm.(M v), [6,12): Rm.(M q)
21 private readonly accRing: Float64Array; // (DECIMATION_SPAN, N) partial audio samples of the normal acceleration
22 private readonly elemHist: Float64Array; // (ringLength, E) finished element-mean samples
23 private readonly emitted: Float64Array; // (N,) the node sample finished at the current step
24 private pending: { j: number; front: number; back: number; acc: Float32Array | null }[] = [];
25 samplesDone = 0;
26 readonly maxChunk = 8192;
27 readonly setup: Setup; readonly movieEvery: number; readonly dispEvery: number; readonly recordAcceleration: boolean;
29 constructor(setup: Setup, opts: BackendOptions | number = {}) {
30 if (typeof opts === "number") opts = { movieEvery: opts };
31 this.setup = setup; this.movieEvery = opts.movieEvery ?? 0; this.dispEvery = opts.dispEvery ?? 0; this.recordAcceleration = opts.recordAcceleration ?? false;
32 const N = setup.mesh.n, E = setup.mesh.e, K = setup.schedule.K, P = setup.schedule.P;
33 this.q = new Float64Array(5 * N); this.v = new Float64Array(5 * N); this.vPrev = new Float64Array(3 * N);
34 this.ef = new Float64Array(15 * E); this.fc = new Float64Array(N);
35 this.zs = new Float64Array(K).fill(HOLD_HEIGHT); this.vs = new Float64Array(K);
36 this.Fk = new Float64Array(K * P); this.F = new Float64Array(K); this.Fmax = new Float64Array(K); this.ncont = new Uint32Array(K);
37 this.accRing = new Float64Array(DECIMATION_SPAN * N); this.elemHist = new Float64Array(E * setup.radiation.ringLength);
38 this.emitted = new Float64Array(N);
39 }
41 private step() {
42 const { shell: sh, mesh, schedule: sc, radiation: rad, dt, substeps: m } = this.setup;
43 const N = mesh.n, n = this.n;
44 const sigma0 = this.setup.config.material.sigma0;
45 // ---- element pass
46 elementForces(sh, this.q, this.v, this.ef);
47 // ---- hit pass
48 this.fc.fill(0);
49 for (let k = 0; k < sc.K; k++) {
50 let wmax = -Infinity;
51 const wp = new Float64Array(sc.P);
52 for (let p = 0; p < sc.P; p++) {
53 const b = (k * sc.P + p) * 3;
54 let w = sc.offsets[k * sc.P + p];
55 for (let i = 0; i < 3; i++) w += sc.weights[b + i] * this.q[5 * sc.nodes[b + i] + 2];
56 wp[p] = w; wmax = Math.max(wmax, w);
57 }
58 if (n === sc.releaseStep[k]) { this.zs[k] = wmax + sc.gap[k]; this.vs[k] = -sc.velocity[k]; }
59 const active = n >= sc.releaseStep[k];
60 let Ftot = 0;
61 for (let p = 0; p < sc.P; p++) {
62 const delta = Math.max(wp[p] - this.zs[k], 0);
63 const f = active ? sc.kPoint[k] * Math.pow(delta, sc.alpha[k]) : 0;
64 this.Fk[k * sc.P + p] = f; Ftot += f;
65 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]; }
66 }
67 this.F[k] = Ftot;
68 }
69 // ---- node pass
70 const { M, Minv, Rm, adjOffsets, adjEntries, normalNode } = sh, L = 5 * N, c = this.partialsPrev;
71 this.partialsNew.fill(0);
72 const f = new Float64Array(5);
73 // decimation bookkeeping for the acceleration of the previous step, s = n - 1
74 const s = n - 1, W = DECIMATION_SPAN, taps = rad.taps;
75 const j0 = s >= 0 ? Math.floor(s / m) : 0, kk = s >= 0 ? s - j0 * m : -1;
76 const emit = s >= 0 && kk === 0, jE = j0 - W / 2, slotE = ((jE % W) + W) % W, r = (((j0 + 1 - W / 2) % W) + W) % W;
77 for (let nd = 0; nd < N; nd++) {
78 f.fill(0);
79 for (let a = adjOffsets[nd]; a < adjOffsets[nd + 1]; a++) {
80 const ent = adjEntries[a], e = (ent / 3) | 0, i = ent % 3;
81 for (let d = 0; d < 5; d++) f[d] += this.ef[15 * e + 5 * i + d];
82 }
83 f[2] += this.fc[nd];
84 let acc = 0;
85 for (let d = 0; d < 5; d++) {
86 const idx = 5 * nd + d;
87 // lagged projection: remove the rigid content that the previous step's update introduced
88 let vd = this.v[idx], qd = this.q[idx];
89 for (let k = 0; k < 6; k++) { const r = Rm[k * L + idx]; vd -= c[k] * r; qd -= c[6 + k] * r; }
90 const v1 = (vd + dt * f[d] * Minv[idx]) / (1 + 2 * sigma0 * dt);
91 // radiated acceleration from consecutive projected velocities, so that the rigid increment
92 // that the lagged projection removes one step later never radiates; normal component
93 if (d < 3) { acc += normalNode[3 * nd + d] * (vd - this.vPrev[3 * nd + d]) / dt; this.vPrev[3 * nd + d] = vd; }
94 const q1 = qd + dt * v1;
95 this.q[idx] = q1; this.v[idx] = v1;
96 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; }
97 }
98 if (s >= 0) {
99 if (emit) { this.emitted[nd] = this.accRing[slotE * N + nd]; this.accRing[slotE * N + nd] = 0; }
100 for (let w = 0; w < W; w++) this.accRing[((r + w) % W) * N + nd] += acc * taps[w * m + kk];
101 }
102 }
103 [this.partialsPrev, this.partialsNew] = [this.partialsNew, this.partialsPrev];
104 // ---- exciters
105 for (let k = 0; k < sc.K; k++) {
106 this.vs[k] += (dt * this.F[k]) / sc.mass[k]; this.zs[k] += dt * this.vs[k];
107 this.Fmax[k] = Math.max(this.Fmax[k], this.F[k]); if (this.F[k] > 0) this.ncont[k]++;
108 }
109 this.n++;
110 if (emit && jE >= 0) this.finalizeSample(jE);
111 }
113 /** Audio sample j of the node accelerations is complete: element means into the history, then the listener pressure. */
114 private finalizeSample(j: number) {
115 const { radiation: rad, mesh } = this.setup, E = mesh.e, D = rad.ringLength, { tri } = mesh;
116 const row = (j % D) * E;
117 for (let e = 0; e < E; e++) this.elemHist[row + e] = (this.emitted[tri[3 * e]] + this.emitted[tri[3 * e + 1]] + this.emitted[tri[3 * e + 2]]) / 3;
118 let front = 0, back = 0;
119 for (let e = 0; e < E; e++) {
120 front += rad.coefF[e] * delayedValue(this.elemHist, E, e, D, j, rad.i0F[e], rad.tapsF, e * DELAY_TAPS);
121 back += rad.coefB[e] * delayedValue(this.elemHist, E, e, D, j, rad.i0B[e], rad.tapsB, e * DELAY_TAPS);
122 }
123 this.pending.push({ j, front, back, acc: this.recordAcceleration ? Float32Array.from(this.emitted) : null });
124 }
126 async runChunk(nSamples: number): Promise<ChunkResult> {
127 const m = this.setup.substeps, N = this.setup.mesh.n, frames: Float32Array[] = [], frameSamples: number[] = [];
128 const dispFrames: Float32Array[] = [], dispSamples: number[] = [];
129 const capture = () => { const w = new Float32Array(N); for (let nd = 0; nd < N; nd++) w[nd] = this.q[5 * nd + 2]; return w; };
130 for (let s = 0; s < nSamples * m; s++) {
131 if (this.movieEvery > 0 && this.n % (m * this.movieEvery) === 0) { frames.push(capture()); frameSamples.push(this.n / m); }
132 if (this.dispEvery > 0 && this.n % (m * this.dispEvery) === 0) { dispFrames.push(capture()); dispSamples.push(this.n / m); }
133 this.step();
134 }
135 this.samplesDone += nSamples;
136 const out = this.pending; this.pending = [];
137 const firstSample = out.length ? out[0].j : this.samplesDone;
138 let acceleration: Float32Array | null = null;
139 if (this.recordAcceleration) {
140 acceleration = new Float32Array(out.length * N);
141 out.forEach((o, i) => acceleration!.set(o.acc!, i * N));
142 }
143 return { firstSample, pressureFront: Float32Array.from(out, (o) => o.front), pressureBack: Float32Array.from(out, (o) => o.back),
144 frames, frameSamples, acceleration, dispFrames, dispSamples };
145 }
147 async diagnostics(): Promise<Diagnostics> {
148 const { shell: sh, schedule: sc, dt } = this.setup, N = sh.mesh.n;
149 let ke = 0, woh = 0;
150 for (let i = 0; i < 5 * N; i++) ke += 0.5 * sh.M[i] * this.v[i] * this.v[i];
151 for (let nd = 0; nd < N; nd++) woh = Math.max(woh, Math.abs(this.q[5 * nd + 2] / sh.hNode[nd]));
152 const { membrane, bending } = strainEnergies(sh, this.q);
153 let ex = 0;
154 for (let k = 0; k < sc.K; k++) ex += 0.5 * sc.mass[k] * this.vs[k] * this.vs[k];
155 return { time: this.n * dt, kinetic: ke, membrane, bending, exciters: ex, maxWoverH: woh,
156 hits: Array.from({ length: sc.K }, (_, k) => ({ peakForce: this.Fmax[k], contactTime: this.ncont[k] * dt, velocity: this.vs[k], height: this.zs[k] })) };
157 }
159 /** M-weighted rigid-mode content of the current velocity (the lagged projection leaves one step's increment). */
160 rigidContent(): number {
161 const { shell: sh } = this.setup, L = 5 * sh.mesh.n;
162 let worst = 0;
163 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)); }
164 return worst;
165 }
167 async displacement(): Promise<Float32Array> {
168 const N = this.setup.mesh.n, w = new Float32Array(N);
169 for (let nd = 0; nd < N; nd++) w[nd] = this.q[5 * nd + 2];
170 return w;
171 }
173 async exportState(): Promise<BackendState> {
174 return { kind: "cpu", step: this.n, samplesDone: this.samplesDone, arrays: {
175 q: this.q.slice(), v: this.v.slice(), vPrev: this.vPrev.slice(), partials: this.partialsPrev.slice(),
176 zs: this.zs.slice(), vs: this.vs.slice(), Fmax: this.Fmax.slice(), ncont: this.ncont.slice(),
177 accRing: this.accRing.slice(), elemHist: this.elemHist.slice() } };
178 }
180 async importState(st: BackendState): Promise<void> {
181 if (st.kind !== "cpu") throw new Error("state was written by another backend");
182 const a = st.arrays;
183 for (const [name, dst] of [["q", this.q], ["v", this.v], ["vPrev", this.vPrev], ["partials", this.partialsPrev], ["zs", this.zs], ["vs", this.vs],
184 ["Fmax", this.Fmax], ["ncont", this.ncont], ["accRing", this.accRing], ["elemHist", this.elemHist]] as const) {
185 if (!a[name] || a[name].length !== dst.length) throw new Error(`state array ${name} is missing or has the wrong size`);
186 (dst as Float64Array).set(a[name] as any);
187 }
188 this.n = st.step; this.samplesDone = st.samplesDone; this.pending = [];
189 }
191 dispose() {}
192}