1// WebGPU version of simulate(): same scheme and result, single precision.
2// See gpu.wgsl for the kernel. The CPU prepares the rest geometry in double
3// precision, chooses the step (prepare.ts), and computes the energies and
4// mode energies from the positions and velocities the GPU returns at each
5// frame.
7import shader from "./gpu.wgsl?raw";
8import { hingeAngle } from "./hinge";
9import { ModeTracker, prepare, staticForces } from "./prepare";
10import { frequencyErrorCents, stepCoefficients } from "./solver";
11import type { SimOptions } from "./solver";
12import type { BeadSystem, SimResult } from "./types";
14const HDR_WORDS = 44;
15const TARGET_MS = 120; // aim for dispatches of about this length (avoids GPU watchdog resets)
17export async function gpuAvailable(): Promise<boolean> {
18 try {
19 return !!(typeof navigator !== "undefined" && navigator.gpu && (await navigator.gpu.requestAdapter()));
20 } catch {
21 return false;
22 }
23}
25// incidence lists: per bead, (slot, coefficient) pairs
26function csr(n: number, entries: [number, number, number][]): { rowptr: Uint32Array; data: Uint32Array } {
27 const count = new Uint32Array(n + 1);
28 for (const [b] of entries) count[b + 1]++;
29 for (let i = 0; i < n; i++) count[i + 1] += count[i];
30 const pos = count.slice(0, n);
31 const data = new Uint32Array(2 * Math.max(1, entries.length));
32 const f = new Float32Array(1);
33 const u = new Uint32Array(f.buffer);
34 for (const [b, slot, coef] of entries) {
35 const q = pos[b]++;
36 data[2 * q] = slot;
37 f[0] = coef;
38 data[2 * q + 1] = u[0];
39 }
40 return { rowptr: count, data };
41}
43export async function simulateGPU(sys: BeadSystem, opts: SimOptions = {}): Promise<SimResult> {
44 const t0 = performance.now();
45 const adapter = await navigator.gpu.requestAdapter();
46 if (!adapter) throw new Error("WebGPU: no adapter");
47 const device = await adapter.requestDevice();
48 try {
49 return await run(device, sys, opts, t0);
50 } finally {
51 device.destroy();
52 }
53}
55async function run(device: GPUDevice, sys: BeadSystem, opts: SimOptions, t0: number): Promise<SimResult> {
56 const { n, pos0, mass, pinned } = sys;
57 const prep = prepare(sys);
58 const warnings = [...prep.warnings];
59 const { substeps, h, nSamples, peRest } = prep;
60 const fs = sys.sampleRate;
61 const Delta = 1 / fs;
62 const cs = stepCoefficients(sys.method);
63 const ns = sys.si.length;
64 const nb = sys.bi.length;
65 const nh = sys.hi.length;
66 const keepD = Array.from(sys.dc.keys()).filter(d => sys.dc[d] !== 0 || sys.dca[d] !== 0);
67 const nd = keepD.length;
69 // ---- parameters (f32) ------------------------------------------------------
70 const pulses = sys.hits.filter(hh => hh.duration > 0);
71 const oInvm = 0;
72 const oMass = n;
73 const oFric = 2 * n;
74 const oS = 3 * n;
75 const oB = oS + 8 * ns;
76 const oH = oB + 8 * nb;
77 const oD = oH + 16 * nh;
78 const oPulse = oD + 8 * nd;
79 const oCs = oPulse + 8 * pulses.length;
80 const P = new Float32Array(oCs + cs.length + 1);
81 for (let i = 0; i < n; i++) {
82 P[oInvm + i] = pinned[i] ? 0 : 1 / mass[i];
83 P[oMass + i] = pinned[i] ? 0 : mass[i];
84 P[oFric + i] = sys.friction[i];
85 }
86 const rel = (a: number, b: number, c: number) => pos0[3 * b + c] - pos0[3 * a + c]; // x_b - x_a at rest
87 for (let s = 0; s < ns; s++) {
88 const i = sys.si[s];
89 const j = sys.sj[s];
90 const r = [rel(i, j, 0), rel(i, j, 1), rel(i, j, 2)];
91 const L = Math.hypot(r[0], r[1], r[2]);
92 P.set([r[0], r[1], r[2], L, L - sys.sL0[s], sys.sk[s], sys.sk3[s]], oS + 8 * s);
93 }
94 for (let t = 0; t < nb; t++) {
95 const [i, j, k] = [sys.bi[t], sys.bj[t], sys.bk[t]];
96 const D = [0, 1, 2].map(c => pos0[3 * i + c] - 2 * pos0[3 * j + c] + pos0[3 * k + c]);
97 P.set([D[0], D[1], D[2], sys.bkb[t], sys.bcb[t]], oB + 8 * t);
98 }
99 for (let q = 0; q < nh; q++) {
100 const [i, j, k, l] = [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]];
101 const R = [i, k, l].flatMap(b => [rel(j, b, 0), rel(j, b, 1), rel(j, b, 2)]);
102 P.set([...R, sys.hkh[q], sys.hch[q]], oH + 16 * q);
103 }
104 keepD.forEach((d, q) => {
105 const [i, j] = [sys.di[d], sys.dj[d]];
106 P.set([rel(i, j, 0), rel(i, j, 1), rel(i, j, 2), sys.dc[d], sys.dca[d]], oD + 8 * q);
107 });
108 pulses.forEach((p, q) => P.set([p.t, p.duration, ...p.force], oPulse + 8 * q));
109 P.set(cs, oCs);
111 // ---- indices (u32) -----------------------------------------------------------
112 const eFb = 3 * ns;
113 const eFh = eFb + 3 * nb;
114 const eDb = 3 * nd;
115 const eDh = eDb + 3 * nb;
116 const forceEntries: [number, number, number][] = [];
117 const dampEntries: [number, number, number][] = [];
118 for (let s = 0; s < ns; s++) {
119 forceEntries.push([sys.si[s], 3 * s, 1], [sys.sj[s], 3 * s, -1]);
120 }
121 for (let t = 0; t < nb; t++) {
122 forceEntries.push([sys.bi[t], eFb + 3 * t, -1], [sys.bj[t], eFb + 3 * t, 2], [sys.bk[t], eFb + 3 * t, -1]);
123 if (sys.bcb[t] !== 0) dampEntries.push([sys.bi[t], eDb + 3 * t, 1], [sys.bj[t], eDb + 3 * t, -2], [sys.bk[t], eDb + 3 * t, 1]);
124 }
125 for (let q = 0; q < nh; q++) {
126 [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]].forEach((b, p) => {
127 forceEntries.push([b, eFh + 12 * q + 3 * p, -1]);
128 if (sys.hch[q] !== 0) dampEntries.push([b, eDh + 12 * q + 3 * p, 1]);
129 });
130 }
131 keepD.forEach((d, q) => dampEntries.push([sys.di[d], 3 * q, 1], [sys.dj[d], 3 * q, -1]));
132 const fcsr = csr(n, forceEntries);
133 const dcsr = csr(n, dampEntries);
134 const pulseInc: number[][] = Array.from({ length: n }, () => []);
135 pulses.forEach((p, q) => p.beads.forEach(b => pulseInc[b].push(q)));
136 const prow = new Uint32Array(n + 1);
137 for (let i = 0; i < n; i++) prow[i + 1] = prow[i] + pulseInc[i].length;
138 const pdat = Uint32Array.from(pulseInc.flat());
139 // impulses, grouped by audio sample
140 const impulses = new Map<number, [number, number, number, number][]>();
141 for (const hit of sys.hits) {
142 if (hit.duration > 0) continue;
143 const s = Math.round(hit.t * fs);
144 if (s >= nSamples) continue;
145 if (!impulses.has(s)) impulses.set(s, []);
146 for (const b of hit.beads) impulses.get(s)!.push([b, ...hit.velocity]);
147 }
148 const impSamples = [...impulses.keys()].sort((a, b) => a - b);
149 const impList: [number, number, number, number][] = [];
150 const impStart = new Map<number, number>();
151 for (const s of impSamples) {
152 impStart.set(s, impList.length);
153 impList.push(...impulses.get(s)!);
154 }
156 const parts: Uint32Array[] = [];
157 let off = 0;
158 const place = (a: Uint32Array) => {
159 const o = off;
160 parts.push(a);
161 off += a.length;
162 return o;
163 };
164 const iS = place(new Uint32Array(Array.from({ length: ns }, (_, s) => [sys.si[s], sys.sj[s]]).flat()));
165 const iB = place(new Uint32Array(Array.from({ length: nb }, (_, t) => [sys.bi[t], sys.bj[t], sys.bk[t]]).flat()));
166 const iH = place(new Uint32Array(Array.from({ length: nh }, (_, q) => [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]]).flat()));
167 const iD = place(new Uint32Array(keepD.flatMap(d => [sys.di[d], sys.dj[d]])));
168 const iFr = place(fcsr.rowptr);
169 const iFe = place(fcsr.data);
170 const iDr = place(dcsr.rowptr);
171 const iDe = place(dcsr.data);
172 const iPr = place(prow);
173 const iPe = place(pdat.length ? pdat : new Uint32Array(1));
174 const iL = place(Uint32Array.from(sys.listenBeads));
175 const impWords = new Uint32Array(4 * Math.max(1, impList.length));
176 {
177 const f = new Float32Array(impWords.buffer);
178 impList.forEach(([b, vx, vy, vz], q) => {
179 impWords[4 * q] = b;
180 f[4 * q + 1] = vx;
181 f[4 * q + 2] = vy;
182 f[4 * q + 3] = vz;
183 });
184 }
185 const iImp = place(impWords);
186 const I = new Uint32Array(off + 1);
187 {
188 let o = 0;
189 for (const a of parts) {
190 I.set(a, o);
191 o += a.length;
192 }
193 }
195 // ---- damping: number of Taylor terms ------------------------------------------
196 let normA = 0;
197 {
198 const row = new Float64Array(n);
199 for (let i = 0; i < n; i++) row[i] = sys.friction[i];
200 keepD.forEach(d => {
201 const c = 2 * (sys.dc[d] + sys.dca[d]);
202 row[sys.di[d]] += c;
203 row[sys.dj[d]] += c;
204 });
205 for (let t = 0; t < nb; t++) {
206 row[sys.bi[t]] += 4 * sys.bcb[t];
207 row[sys.bj[t]] += 8 * sys.bcb[t];
208 row[sys.bk[t]] += 4 * sys.bcb[t];
209 }
210 const g = new Float64Array(12);
211 for (let q = 0; q < nh; q++) {
212 if (sys.hch[q] === 0) continue;
213 hingeAngle(pos0, sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q], g);
214 const total = g.reduce((a, x) => a + Math.abs(x), 0);
215 // 1.5x margin: the gradients change as the shell moves
216 [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]].forEach((b, p) => {
217 row[b] += 1.5 * sys.hch[q] * total * (Math.abs(g[3 * p]) + Math.abs(g[3 * p + 1]) + Math.abs(g[3 * p + 2]));
218 });
219 }
220 for (let i = 0; i < n; i++) if (!pinned[i]) normA = Math.max(normA, row[i] / mass[i]);
221 }
222 const hasDamp = normA > 0;
223 const nSplit = Math.max(1, Math.ceil((Delta * normA) / 0.5));
224 let K = 1;
225 {
226 const z = (Delta / nSplit) * normA;
227 let term = z;
228 while (term > 1e-9 && K < 30) {
229 K++;
230 term *= z / K;
231 }
232 }
234 // ---- buffers ---------------------------------------------------------------------
235 const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
236 const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
237 const frameSize = 7 * n; // u, v, dissipated energy per bead
238 const maxChunk = Math.min(nSamples, 8192);
239 const framesPerChunk = Math.floor(maxChunk / frameStride) + 2;
240 const oFrames = maxChunk;
241 const outWords = oFrames + framesPerChunk * frameSize;
242 const usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC;
243 const mk = (words: number, u = usage) => device.createBuffer({ size: Math.max(16, 4 * words), usage: u });
244 const bufH = device.createBuffer({ size: 4 * HDR_WORDS, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
245 const bufP = mk(P.length);
246 const bufI = mk(I.length);
247 const bufS = mk(17 * n + nh + 1);
248 const bufE = mk(Math.max(3 * ns + 3 * nb + 12 * nh, 3 * nd + 3 * nb + 12 * nh, 1));
249 const bufO = mk(outWords);
250 const bufC = mk(4);
251 const staging = mk(outWords, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
252 device.queue.writeBuffer(bufP, 0, P);
253 device.queue.writeBuffer(bufI, 0, I);
255 // per-phase constants, selected with a dynamic offset
256 const align = device.limits.minUniformBufferOffsetAlignment;
257 const phases: { kick: number; drift: number; tOff: number; gather: number; c: number; saveAv: number; last: number }[] = [];
258 const addPhase = (p: Partial<(typeof phases)[number]>) => {
259 phases.push({ kick: 0, drift: 0, tOff: 0, gather: 0, c: 0, saveAv: 0, last: 0, ...p });
260 return (phases.length - 1) * align;
261 };
262 // During the sub-steps of sample s the sample index on the GPU is already
263 // s + 1 (the record phase advanced it), so t = (s + 1) Delta + tOff.
264 const m = substeps * cs.length;
265 const csAt = (j: number) => cs[j % cs.length];
266 const phStart = addPhase({ kick: 0.5 * csAt(0) * h, drift: csAt(0) * h });
267 const phMid: number[] = [];
268 let tAcc = csAt(0) * h;
269 for (let j = 1; j < m; j++) {
270 phMid.push(addPhase({ gather: 1, kick: 0.5 * (csAt(j - 1) + csAt(j)) * h, drift: csAt(j) * h, tOff: tAcc - Delta }));
271 tAcc += csAt(j) * h;
272 }
273 const phEnd = addPhase({ gather: 1, kick: 0.5 * csAt(m - 1) * h, tOff: tAcc - Delta });
274 const phInit = addPhase({ gather: 1, tOff: -Delta });
275 const phDamp: number[][] = [];
276 for (let p = 0; p < nSplit; p++) {
277 phDamp.push(
278 Array.from({ length: K }, (_, k) =>
279 addPhase({ c: -Delta / nSplit / (k + 1), saveAv: p === 0 && k === 0 ? 1 : 0, last: p === nSplit - 1 && k === K - 1 ? 1 : 0 })
280 )
281 );
282 }
283 const bufQ = device.createBuffer({ size: phases.length * align, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
284 {
285 const q = new ArrayBuffer(phases.length * align);
286 phases.forEach((ph, r) => {
287 const f = new Float32Array(q, r * align, 8);
288 const u = new Uint32Array(q, r * align, 8);
289 f[0] = ph.kick;
290 f[1] = ph.drift;
291 f[2] = ph.tOff;
292 u[3] = ph.gather;
293 f[4] = ph.c;
294 u[5] = ph.saveAv;
295 u[6] = ph.last;
296 });
297 device.queue.writeBuffer(bufQ, 0, q);
298 }
300 const module = device.createShaderModule({ code: shader });
301 const info = await module.getCompilationInfo();
302 const errs = info.messages.filter(msg => msg.type === "error");
303 if (errs.length) throw new Error("WebGPU shader: " + errs.map(msg => `${msg.lineNum}:${msg.linePos} ${msg.message}`).join("; "));
304 const storage = (binding: number, type: GPUBufferBindingType) => ({ binding, visibility: GPUShaderStage.COMPUTE, buffer: { type } });
305 const layout0 = device.createBindGroupLayout({
306 entries: [
307 { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" as const } },
308 storage(1, "read-only-storage"),
309 storage(2, "read-only-storage"),
310 storage(3, "storage"),
311 storage(4, "storage"),
312 storage(5, "storage"),
313 storage(6, "storage"),
314 ],
315 });
316 const layout1 = device.createBindGroupLayout({
317 entries: [{ binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform", hasDynamicOffset: true, minBindingSize: 32 } }],
318 });
319 const plLayout = device.createPipelineLayout({ bindGroupLayouts: [layout0, layout1] });
320 const pipe = (entryPoint: string) => device.createComputePipeline({ layout: plLayout, compute: { module, entryPoint } });
321 const K_ = {
322 elements: pipe("elements"),
323 beads: pipe("beads"),
324 pre: pipe("pre"),
325 dampElements: pipe("dampElements"),
326 dampBeads: pipe("dampBeads"),
327 record: pipe("record"),
328 impulse: pipe("impulse"),
329 initHinges: pipe("initHinges"),
330 resetT: pipe("resetT"),
331 };
332 const bind0 = device.createBindGroup({
333 layout: layout0,
334 entries: [bufH, bufP, bufI, bufS, bufE, bufO, bufC].map((buffer, binding) => ({ binding, resource: { buffer } })),
335 });
336 const bind1 = device.createBindGroup({ layout: layout1, entries: [{ binding: 0, resource: { buffer: bufQ, size: 32 } }] });
338 const hdr = new ArrayBuffer(4 * HDR_WORDS);
339 const hu = new Uint32Array(hdr);
340 const hf = new Float32Array(hdr);
341 const quantity = { displacement: 0, velocity: 1, acceleration: 2 }[sys.listenQuantity];
342 const setHdr = (sample0: number, nImp: number, impOff: number) => {
343 hu.set([n, ns, nb, nh, nd, sys.listenBeads.length, quantity, pulses.length]);
344 hu.set([sample0, frameStride, frameSize, nImp, impOff, oFrames, 0, 0], 8);
345 hu.set([oInvm, oMass, oFric, oS, oB, oH, oD, oPulse], 16);
346 hu.set([iS, iB, iH, iD, iFr, iFe, iDr, iDe, iPr, iPe, iL, iImp], 24);
347 hu.set([eFb, eFh, eDb, eDh], 36);
348 hf.set([Delta, ...sys.listenDir], 40);
349 device.queue.writeBuffer(bufH, 0, hdr);
350 };
351 const groups = (count: number) => Math.max(1, Math.ceil(count / 64));
352 const gBeads = groups(n);
353 const gElems = groups(ns + nb + nh);
354 const gDamp = groups(nd + nb + nh);
355 // only the kernels that read the phase table need it rebound
356 let boundPhase = -1;
357 const go = (pass: GPUComputePassEncoder, p: GPUComputePipeline, nGroups: number, phase = -1) => {
358 pass.setPipeline(p);
359 if (phase >= 0 && phase !== boundPhase) {
360 pass.setBindGroup(1, bind1, [phase]);
361 boundPhase = phase;
362 }
363 pass.dispatchWorkgroups(nGroups);
364 };
365 const encodeSample = (pass: GPUComputePassEncoder) => {
366 go(pass, K_.pre, gBeads);
367 if (quantity !== 2 || !hasDamp) go(pass, K_.record, 1);
368 if (hasDamp) {
369 phDamp.forEach((terms, p) => {
370 if (p > 0) go(pass, K_.resetT, gBeads); // T = v again (pieces of a long damping step)
371 terms.forEach((ph, k) => {
372 go(pass, K_.dampElements, gDamp);
373 go(pass, K_.dampBeads, gBeads, ph);
374 if (quantity === 2 && p === 0 && k === 0) go(pass, K_.record, 1);
375 });
376 });
377 }
378 go(pass, K_.beads, gBeads, phStart);
379 for (const ph of phMid) {
380 go(pass, K_.elements, gElems);
381 go(pass, K_.beads, gBeads, ph);
382 }
383 go(pass, K_.elements, gElems);
384 go(pass, K_.beads, gBeads, phEnd);
385 };
387 // initial state: rest angles of the hinges, then the forces at t = 0
388 setHdr(0, 0, 0);
389 device.queue.writeBuffer(bufC, 0, Uint32Array.of(1));
390 {
391 const enc = device.createCommandEncoder();
392 const pass = enc.beginComputePass();
393 pass.setBindGroup(0, bind0);
394 pass.setBindGroup(1, bind1, [0]);
395 boundPhase = 0;
396 go(pass, K_.initHinges, groups(nh));
397 go(pass, K_.elements, gElems);
398 go(pass, K_.beads, gBeads, phInit);
399 pass.end();
400 device.queue.submit([enc.finish()]);
401 }
403 // ---- outputs --------------------------------------------------------------------
404 const audio = new Float32Array(nSamples);
405 const frameTimes = new Float64Array(nFrames);
406 const energy = new Float64Array(nFrames);
407 const dissipated = new Float64Array(nFrames);
408 const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
409 const animStride = frameStride * Math.max(1, Math.ceil((nFrames * 3 * n * 4) / Math.max(maxAnimBytes, 1)));
410 const animTimes: number[] = [];
411 const animFrames: Float32Array[] = [];
412 const modes = opts.modes ?? null;
413 const nModes = modes ? modes.freqs.length : 0;
414 const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
415 const tracker = modes ? new ModeTracker(modes, mass, pos0) : null;
416 const xf = new Float64Array(3 * n);
417 const vf = new Float64Array(3 * n);
418 const Fs = new Float64Array(3 * n);
419 let frame = 0;
420 let aborted: string | null = null;
422 let sample = 0;
423 let chunk = 16;
424 let lastProgress = performance.now();
425 while (sample < nSamples) {
426 let len = Math.min(chunk, maxChunk, nSamples - sample);
427 const nextImp = impSamples.find(s => s > sample);
428 if (nextImp !== undefined) len = Math.min(len, nextImp - sample);
429 const imps = impulses.get(sample);
430 const tStart = performance.now();
431 setHdr(sample, imps ? imps.length : 0, imps ? impStart.get(sample)! : 0);
432 device.queue.writeBuffer(bufC, 0, Uint32Array.of(sample));
433 const enc = device.createCommandEncoder();
434 const pass = enc.beginComputePass();
435 pass.setBindGroup(0, bind0);
436 pass.setBindGroup(1, bind1, [0]);
437 boundPhase = 0;
438 if (imps) go(pass, K_.impulse, groups(imps.length));
439 for (let q = 0; q < len; q++) encodeSample(pass);
440 pass.end();
441 enc.copyBufferToBuffer(bufO, 0, staging, 0, 4 * outWords);
442 device.queue.submit([enc.finish()]);
443 await staging.mapAsync(GPUMapMode.READ);
444 const out = new Float32Array(staging.getMappedRange().slice(0));
445 staging.unmap();
446 const elapsed = performance.now() - tStart;
448 audio.set(out.subarray(0, len), sample);
449 for (let q = 0; q < len; q++) {
450 if (!Number.isFinite(out[q])) {
451 aborted = `The solution blew up near t = ${((sample + q) / fs).toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
452 break;
453 }
454 }
455 // frames in this chunk
456 const first = Math.ceil(sample / frameStride) * frameStride;
457 for (let fsmp = first; fsmp < sample + len && frame < nFrames && !aborted; fsmp += frameStride) {
458 const o = oFrames + ((fsmp - first) / frameStride) * frameSize;
459 let ke = 0;
460 for (let q = 0; q < 3 * n; q++) {
461 xf[q] = pos0[q] + out[o + q];
462 vf[q] = out[o + 3 * n + q];
463 }
464 for (let i = 0; i < n; i++) ke += 0.5 * (pinned[i] ? 0 : mass[i]) * (vf[3 * i] ** 2 + vf[3 * i + 1] ** 2 + vf[3 * i + 2] ** 2);
465 frameTimes[frame] = fsmp / fs;
466 energy[frame] = ke + staticForces(sys, xf, Fs) - peRest;
467 let diss = 0;
468 for (let i = 0; i < n; i++) diss += out[o + 6 * n + i];
469 dissipated[frame] = diss;
470 if (tracker && modeEnergy) tracker.record(out.subarray(o, o + 3 * n), out.subarray(o + 3 * n, o + 6 * n), modeEnergy, frame * nModes, true);
471 if (fsmp % animStride === 0) {
472 animTimes.push(fsmp / fs);
473 animFrames.push(out.slice(o, o + 3 * n));
474 }
475 opts.onFrame?.(frame, fsmp / fs, xf, vf);
476 frame++;
477 }
478 if (aborted) break;
479 sample += len;
480 // adapt the chunk so each dispatch takes about TARGET_MS
481 chunk = Math.max(16, Math.min(maxChunk, Math.round((len * TARGET_MS) / Math.max(elapsed, 1))));
482 const now = performance.now();
483 if (opts.onProgress && now - lastProgress > 200) {
484 opts.onProgress(sample / fs);
485 lastProgress = now;
486 }
487 }
489 // energy drift after the last excitation
490 const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
491 let energyDrift = 0;
492 let ref = -1;
493 for (let f = 0; f < frame; f++) {
494 if (frameTimes[f] < tExcitationEnd) continue;
495 const total = energy[f] + dissipated[f];
496 if (ref < 0) ref = total;
497 else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
498 }
500 return {
501 audio,
502 sampleRate: fs,
503 dt: h,
504 substeps,
505 omegaMax: prep.omegaMax,
506 frameTimes: frameTimes.slice(0, frame),
507 energy: energy.slice(0, frame),
508 dissipated: dissipated.slice(0, frame),
509 animTimes: Float64Array.from(animTimes),
510 animFrames,
511 modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
512 residualForce: prep.residualForce,
513 maxTension: prep.maxTension,
514 energyDrift,
515 freqErrors: [1000, 2000, 5000, 10000, 20000]
516 .filter(hz => hz < fs / 2)
517 .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
518 elapsedMs: performance.now() - t0,
519 aborted,
520 warnings,
521 };
522}