WebGPU toggle: single-precision GPU integration (displacement-based, multi-dispatch phases, frames read back), linear modes behind a toggle (off by default), shared prepare()
13 changed files+1094−73
README.mdmodified+8−2View file
@@ -9,7 +9,7 @@ The aim is to reproduce the phenomena of drums and cymbals (resonant modes, nonl
99 ```bash
1010 npm install
1111 npm run dev # the app, at the printed localhost URL
12-npm test # solver tests (Node)
12+npm test # solver tests (Node; the WebGPU path is checked in the browser with src/devcompare.ts)
1313 ```
1414
1515 ## The model
@@ -43,9 +43,15 @@ Each run reports:
4343 - the energy drift after the last excitation, with the energy removed by damping added back (it is counted exactly);
4444 - optionally, the relative difference of the listener signal from a run at dt/2 over [0, t] for several t. In a chaotic run this difference grows with t no matter how small dt is, so the short windows are the meaningful ones.
4545
46+## WebGPU
47+
48+The "WebGPU" toggle runs the time integration on the GPU (`src/sim/gpu.ts`, `gpu.wgsl`): the same scheme as on the CPU, in single precision, since WebGPU has no double precision. To keep single precision accurate, the GPU stores each bead's displacement from rest rather than its position, every element works from its rest geometry (computed in double precision on the CPU) plus displacement differences, and a spring's stretch is computed as (2 r0 . du + |du|^2)/(L + |r0|), without subtracting nearly equal lengths. Each phase of a step (element forces, bead update, damping terms, recording) is a dispatch over beads or elements; the phases for a chunk of audio samples are encoded into one command buffer, and the CPU computes energies and mode energies from the positions and velocities read back at each frame.
49+
50+Compared with the CPU (`src/devcompare.ts`, in Chrome on a laptop with an Intel Iris Xe GPU), the listener signals agree to single precision: relative differences of 2e-5 (drum) to 3e-3 (plate, dome) over the first 10 to 30 ms, with matching energy drift and dissipated energy; in a chaotic run the difference then grows like any other perturbation. The GPU has a floor of about 11 to 15x real time, the cost of about 15 dispatches per audio sample, so it is slower than the CPU for small nets (the 122-bead string, 14x against 3x) and faster for large ones: 367-bead dome 17x against 39x, 1015-bead plate 20x against 50x, a 6481-bead drum 97x against 304x. The perturbed twin needs double precision and is skipped when WebGPU is on.
51+
4652 ## Linear modes
4753
48-The stiffness at rest is assembled exactly from the elements. For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, so its modes need only the z coordinates; the app computes them with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example). For a curved net (a dome) the modes need all three coordinates of every free bead (up to 1200 coordinates, about 12 s). The app records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
54+Linear modes are computed only when the "linear modes" toggle is on (off by default), or when the script sets `sys.modes = true`; the mode-energy and decay-time plots need them. The stiffness at rest is assembled exactly from the elements. For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, so its modes need only the z coordinates; the app computes them with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example). For a curved net (a dome) the modes need all three coordinates of every free bead (up to 1200 coordinates, about 12 s). The app records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
4955
5056 The "small-strike reference" option runs the same system with every hit scaled down by 1000 and the output scaled back up, which is what the system would do if it were linear. The app plays it for comparison and reports the difference from the actual run.
5157
package-lock.jsonmodified+8−0View file
@@ -22,6 +22,7 @@
2222 "@types/react": "^19.3.0",
2323 "@types/react-dom": "^19.3.0",
2424 "@vitejs/plugin-react": "^6.1.2",
25+ "@webgpu/types": "^0.1.74",
2526 "tsx": "^4.23.15",
2627 "typescript": "^7.0.2",
2728 "vite": "^8.3.3",
@@ -1438,6 +1439,13 @@
14381439 "url": "https://opencollective.com/vitest"
14391440 }
14401441 },
1442+ "node_modules/@webgpu/types": {
1443+ "version": "0.1.74",
1444+ "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.74.tgz",
1445+ "integrity": "sha512-lgiI4hbuLcI9unnm2cL/tvCaQU45dp0xcLWJh5uB/9MBGvIA4F8XIk7nSiBeg+K4xWGYwgZKn80LmERI5CxoTA==",
1446+ "dev": true,
1447+ "license": "BSD-3-Clause"
1448+ },
14411449 "node_modules/assertion-error": {
14421450 "version": "2.0.1",
14431451 "resolved": "https://registry.npmjs.org/assertion-error/-/assertion-error-2.0.1.tgz",
package.jsonmodified+1−0View file
@@ -25,6 +25,7 @@
2525 "@types/react": "^19.3.0",
2626 "@types/react-dom": "^19.3.0",
2727 "@vitejs/plugin-react": "^6.1.2",
28+ "@webgpu/types": "^0.1.74",
2829 "tsx": "^4.23.15",
2930 "typescript": "^7.0.2",
3031 "vite": "^8.3.3",
src/App.tsxmodified+42−4View file
@@ -13,6 +13,28 @@ const examples = Object.entries(exampleSources)
1313 const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
1414
1515 const STORE_KEY = "bead-percussion:script";
16+
17+// a boolean remembered in this browser (falls back to the default if storage is unavailable)
18+function usePersistent(name: string, initial: boolean): [boolean, (v: boolean) => void] {
19+ const key = `bead-percussion:${name}`;
20+ const [value, setValue] = useState(() => {
21+ try {
22+ const s = localStorage.getItem(key);
23+ return s === null ? initial : s === "1";
24+ } catch {
25+ return initial;
26+ }
27+ });
28+ const set = (v: boolean) => {
29+ setValue(v);
30+ try {
31+ localStorage.setItem(key, v ? "1" : "0");
32+ } catch {
33+ /* storage unavailable */
34+ }
35+ };
36+ return [value, set];
37+}
1638 const load = () => {
1739 try {
1840 return localStorage.getItem(STORE_KEY);
@@ -42,6 +64,9 @@ export default function App() {
4264 const [dtCheck, setDtCheck] = useState(false);
4365 const [linearRef, setLinearRef] = useState(false);
4466 const [twin, setTwin] = useState(false);
67+ const [modesOn, setModesOn] = usePersistent("modes", false);
68+ const [gpuOn, setGpuOn] = usePersistent("gpu", false);
69+ const gpuAvailable = typeof navigator !== "undefined" && "gpu" in navigator;
4570 const [playing, setPlaying] = useState<"main" | "linear" | "twin" | null>(null);
4671 const [runId, setRunId] = useState(0); // remounts panels with per-run state
4772 const worker = useRef<Worker | null>(null);
@@ -94,10 +119,10 @@ export default function App() {
94119 setError(e.message || "Worker error");
95120 setRunning(false);
96121 };
97- const req: RunRequest = { type: "run", script, dtCheck, linear: linearRef, twin };
122+ const req: RunRequest = { type: "run", script, dtCheck, linear: linearRef, twin, modes: modesOn, gpu: gpuOn && gpuAvailable };
98123 w.postMessage(req);
99124 },
100- [script, stopAudio, dtCheck, linearRef, twin]
125+ [script, stopAudio, dtCheck, linearRef, twin, modesOn, gpuOn, gpuAvailable]
101126 );
102127
103128 const stop = () => {
@@ -196,6 +221,18 @@ export default function App() {
196221 )}
197222 </div>
198223 <div className="controls options">
224+ <label title="Compute the linear modes (eigenvalues) and track the energy in each; can take seconds for large or curved nets">
225+ <input type="checkbox" checked={modesOn} onChange={e => setModesOn(e.target.checked)} /> linear modes
226+ </label>
227+ <label
228+ title={
229+ gpuAvailable
230+ ? "Integrate on the GPU with WebGPU (single precision). Faster for large nets; the perturbed twin still runs on the CPU."
231+ : "WebGPU is not available in this browser"
232+ }
233+ >
234+ <input type="checkbox" disabled={!gpuAvailable} checked={gpuOn && gpuAvailable} onChange={e => setGpuOn(e.target.checked)} /> WebGPU
235+ </label>
199236 <label title="Run again at dt/2 and compare the listener signals">
200237 <input type="checkbox" checked={dtCheck} onChange={e => setDtCheck(e.target.checked)} /> dt/2 check
201238 </label>
@@ -304,7 +341,7 @@ function Report({ p }: { p: RunPayload }) {
304341 const rows: [string, string][] = [
305342 ["beads", `${p.sys.n} (${nPinned} pinned)`],
306343 ["springs, bends, hinges, dashpots", `${p.sys.si.length}, ${p.sys.nBends}, ${p.sys.nHinges}, ${p.sys.nDashpots}`],
307- ["method", p.sys.method],
344+ ["method", `${p.sys.method} (${p.sys.engine})`],
308345 ["dt", `${r.dt.toExponential(3)} s (${r.substeps} per audio sample)`],
309346 ["highest linear frequency", `${(r.omegaMax / (2 * Math.PI)).toFixed(0)} Hz`],
310347 ["frequency error", r.freqErrors.map(e => `${e.cents.toFixed(e.cents < 0.1 ? 3 : 2)}¢ at ${e.hz / 1000} kHz`).join(", ")],
@@ -413,5 +450,6 @@ Solver (optional)
413450 sys.dt fixed step (rounded to divide the audio sample interval)
414451 sys.sampleRate audio rate (default 48000)
415452 sys.frameRate energy/mode/animation records per second (default 500)
416- sys.modes compute linear transverse modes (default true; flat systems)
453+ sys.modes true/false: compute linear modes, overriding the page's
454+ "linear modes" toggle (off by default)
417455 Bead indices are 1-based.`;
src/devcompare.tsadded+35−0View file
@@ -0,0 +1,35 @@
1+// Development check: run an example on the CPU and on the GPU and compare.
2+// From the browser console: (await import('/src/devcompare.ts')).compare('drum', 0.1)
3+
4+import { relDiff } from "./sim/analysis";
5+import { simulateGPU } from "./sim/gpu";
6+import { runScript } from "./sim/runScript";
7+import { simulate } from "./sim/solver";
8+
9+const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
10+const exampleSources = import.meta.glob("../matlab/examples/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
11+const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
12+
13+export async function compare(name: string, duration: number, edit: (s: string) => string = s => s) {
14+ const src = edit(exampleSources[`../matlab/examples/${name}.m`]);
15+ const { sys } = runScript(src, helpers, () => {});
16+ sys.duration = duration;
17+ const cpu = simulate(sys);
18+ const gpu = await simulateGPU(sys);
19+ const fs = cpu.sampleRate;
20+ const peak = (a: Float32Array) => a.reduce((m, x) => Math.max(m, Math.abs(x)), 0);
21+ return {
22+ n: sys.n,
23+ substeps: cpu.substeps,
24+ cpuRealTime: cpu.elapsedMs / 1000 / duration,
25+ gpuRealTime: gpu.elapsedMs / 1000 / duration,
26+ aborted: gpu.aborted,
27+ diff: [0.001, 0.003, 0.01, 0.03, 0.1, duration].filter(t => t <= duration).map(t => [t, relDiff(cpu.audio, gpu.audio, Math.round(t * fs))]),
28+ peakCpu: peak(cpu.audio),
29+ peakGpu: peak(gpu.audio),
30+ driftCpu: cpu.energyDrift,
31+ driftGpu: gpu.energyDrift,
32+ energyEnd: [cpu.energy[cpu.energy.length - 1], gpu.energy[gpu.energy.length - 1]],
33+ dissEnd: [cpu.dissipated[cpu.dissipated.length - 1], gpu.dissipated[gpu.dissipated.length - 1]],
34+ };
35+}
src/sim/gpu.tsadded+522−0View file
@@ -0,0 +1,522 @@
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.
6+
7+import shader from "./gpu.wgsl?raw";
8+import { hingeAngle } from "./hinge";
9+import { ModeTracker, prepare, staticForces } from "./prepare";
10+import { frequencyErrorCents, stepCoefficients } from "./solver";
11+import type { SimOptions } from "./solver";
12+import type { BeadSystem, SimResult } from "./types";
13+
14+const HDR_WORDS = 44;
15+const TARGET_MS = 120; // aim for dispatches of about this length (avoids GPU watchdog resets)
16+
17+export 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+}
24+
25+// incidence lists: per bead, (slot, coefficient) pairs
26+function 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+}
42+
43+export 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+}
54+
55+async 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;
68+
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);
110+
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+ }
155+
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+ }
194+
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+ }
233+
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);
254+
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+ }
299+
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 } }] });
337+
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+ };
386+
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+ }
402+
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;
421+
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;
447+
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+ }
488+
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+ }
499+
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+}
src/sim/gpu.wgsladded+291−0View file
@@ -0,0 +1,291 @@
1+// Bead solver on the GPU: the same scheme as solver.ts (velocity Verlet or
2+// its Yoshida composition, with exact damping applied once per audio sample),
3+// in single precision.
4+//
5+// Each phase of a step is its own dispatch over beads or elements; the host
6+// encodes the phases for a chunk of audio samples into one command buffer.
7+// Per-phase constants (kick and drift lengths, Taylor coefficients) come from
8+// a table selected with a dynamic uniform offset, and the current sample
9+// index lives in a GPU buffer, which the record phase advances.
10+//
11+// To keep single precision accurate, the state is the displacement u from
12+// the rest positions, never the absolute position, and every element works
13+// from its rest geometry (computed in double precision on the CPU) plus
14+// displacement differences. Element forces are written to per-element slots
15+// and gathered per bead through an incidence list, since there are no atomic
16+// float additions.
17+
18+const WG: u32 = 64u;
19+
20+struct Hdr {
21+ n: u32, ns: u32, nb: u32, nh: u32,
22+ nd: u32, nl: u32, quantity: u32, npulse: u32,
23+ sample0: u32, frameStride: u32, frameSize: u32, nImp: u32,
24+ impOff: u32, oFrames: u32, pad0: u32, pad1: u32,
25+ // offsets into P (f32 parameters)
26+ oInvm: u32, oMass: u32, oFric: u32, oS: u32,
27+ oB: u32, oH: u32, oD: u32, oPulse: u32,
28+ // offsets into I (u32 indices)
29+ iS: u32, iB: u32, iH: u32, iD: u32,
30+ iFr: u32, iFe: u32, iDr: u32, iDe: u32,
31+ iPr: u32, iPe: u32, iL: u32, iImp: u32,
32+ // offsets into E for the bending and hinge slots (forces / damping)
33+ eFb: u32, eFh: u32, eDb: u32, eDh: u32,
34+ Delta: f32, dirx: f32, diry: f32, dirz: f32,
35+}
36+
37+// per-phase constants
38+struct Ph {
39+ kick: f32, drift: f32, tOff: f32, gather: u32,
40+ c: f32, saveAv: u32, last: u32, pad: u32,
41+}
42+
43+@group(0) @binding(0) var<uniform> H: Hdr;
44+@group(0) @binding(1) var<storage, read> P: array<f32>;
45+@group(0) @binding(2) var<storage, read> I: array<u32>;
46+@group(0) @binding(3) var<storage, read_write> S: array<f32>; // u | v | F | T | Av | KE0 | diss | theta0
47+@group(0) @binding(4) var<storage, read_write> E: array<f32>; // element slots
48+@group(0) @binding(5) var<storage, read_write> O: array<f32>; // audio | frames
49+@group(0) @binding(6) var<storage, read_write> C: array<u32>; // [current sample]
50+@group(1) @binding(0) var<uniform> Q: Ph;
51+
52+fn offV() -> u32 { return 3u * H.n; }
53+fn offF() -> u32 { return 6u * H.n; }
54+fn offT() -> u32 { return 9u * H.n; }
55+fn offAv() -> u32 { return 12u * H.n; }
56+fn offKE() -> u32 { return 15u * H.n; }
57+fn offDiss() -> u32 { return 16u * H.n; }
58+fn offTh0() -> u32 { return 17u * H.n; }
59+
60+fn S3(o: u32) -> vec3<f32> { return vec3<f32>(S[o], S[o + 1u], S[o + 2u]); }
61+fn setS3(o: u32, w: vec3<f32>) { S[o] = w.x; S[o + 1u] = w.y; S[o + 2u] = w.z; }
62+fn setE3(o: u32, w: vec3<f32>) { E[o] = w.x; E[o + 1u] = w.y; E[o + 2u] = w.z; }
63+fn E3(o: u32) -> vec3<f32> { return vec3<f32>(E[o], E[o + 1u], E[o + 2u]); }
64+fn P3(o: u32) -> vec3<f32> { return vec3<f32>(P[o], P[o + 1u], P[o + 2u]); }
65+fn U(b: u32) -> vec3<f32> { return S3(3u * b); }
66+fn V(b: u32) -> vec3<f32> { return S3(offV() + 3u * b); }
67+
68+// Hinge h: angle and gradient from the rest vectors relative to bead j plus
69+// displacement differences.
70+struct HingeG { theta: f32, g: array<vec3<f32>, 4> }
71+
72+fn hingeAt(h: u32) -> HingeG {
73+ let bi = I[H.iH + 4u * h];
74+ let bj = I[H.iH + 4u * h + 1u];
75+ let bk = I[H.iH + 4u * h + 2u];
76+ let bl = I[H.iH + 4u * h + 3u];
77+ let o = H.oH + 16u * h;
78+ let uj = U(bj);
79+ let pa = P3(o) + (U(bi) - uj); // bead i
80+ let pk = P3(o + 3u) + (U(bk) - uj); // bead k
81+ let pl = P3(o + 6u) + (U(bl) - uj); // bead l; bead j is the origin
82+ let e = pk;
83+ let eLen = length(e);
84+ let b = pa - pk;
85+ let c = pl - pk;
86+ let n1 = cross(pa, b);
87+ let n2 = cross(c, pl);
88+ let n1sq = dot(n1, n1);
89+ let n2sq = dot(n2, n2);
90+ var r: HingeG;
91+ r.theta = -atan2(dot(cross(n1, n2), e) / eLen, dot(n1, n2));
92+ let p1 = dot(pa, e) / eLen;
93+ let q1 = dot(b, e) / eLen;
94+ let p2 = dot(pl, e) / eLen;
95+ let q2 = dot(c, e) / eLen;
96+ let m1 = 1.0 / n1sq;
97+ let m2 = 1.0 / n2sq;
98+ r.g[0] = (eLen * m1) * n1;
99+ r.g[1] = q1 * m1 * n1 + q2 * m2 * n2;
100+ r.g[2] = -p1 * m1 * n1 - p2 * m2 * n2;
101+ r.g[3] = (eLen * m2) * n2;
102+ return r;
103+}
104+
105+fn wrapAngle(d: f32) -> f32 {
106+ var x = d;
107+ if (x > 3.14159265) { x -= 6.28318531; }
108+ if (x <= -3.14159265) { x += 6.28318531; }
109+ return x;
110+}
111+
112+// ---- forces ---------------------------------------------------------------------
113+
114+// one invocation per element: springs, then bending triplets, then hinges
115+@compute @workgroup_size(WG)
116+fn elements(@builtin(global_invocation_id) gid: vec3<u32>) {
117+ var q = gid.x;
118+ if (q < H.ns) {
119+ let o = H.oS + 8u * q;
120+ let r0 = P3(o);
121+ let du = U(I[H.iS + 2u * q + 1u]) - U(I[H.iS + 2u * q]);
122+ let r = r0 + du;
123+ let L = length(r);
124+ // L - |r0| without cancellation, plus the rest stretch |r0| - L0
125+ let e = P[o + 4u] + (2.0 * dot(r0, du) + dot(du, du)) / (L + P[o + 3u]);
126+ let f = P[o + 5u] * e + P[o + 6u] * e * e * e;
127+ setE3(3u * q, (f / L) * r);
128+ return;
129+ }
130+ q -= H.ns;
131+ if (q < H.nb) {
132+ let o = H.oB + 8u * q;
133+ let D = P3(o) + U(I[H.iB + 3u * q]) - 2.0 * U(I[H.iB + 3u * q + 1u]) + U(I[H.iB + 3u * q + 2u]);
134+ setE3(H.eFb + 3u * q, P[o + 3u] * D);
135+ return;
136+ }
137+ q -= H.nb;
138+ if (q < H.nh) {
139+ let hg = hingeAt(q);
140+ let s = P[H.oH + 16u * q + 9u] * wrapAngle(hg.theta - S[offTh0() + q]);
141+ for (var p = 0u; p < 4u; p++) { setE3(H.eFh + 12u * q + 3u * p, s * hg.g[p]); }
142+ }
143+}
144+
145+fn pulseForce(b: u32, t: f32) -> vec3<f32> {
146+ var f = vec3<f32>(0.0);
147+ for (var q = I[H.iPr + b]; q < I[H.iPr + b + 1u]; q++) {
148+ let o = H.oPulse + 8u * I[H.iPe + q];
149+ let x = (t - P[o]) / P[o + 1u];
150+ if (x > 0.0 && x < 1.0) { f += (0.5 * (1.0 - cos(6.28318531 * x))) * P3(o + 2u); }
151+ }
152+ return f;
153+}
154+
155+// gather (optional), kick and drift; time t = (sample index) * Delta + tOff
156+@compute @workgroup_size(WG)
157+fn beads(@builtin(global_invocation_id) gid: vec3<u32>) {
158+ let b = gid.x;
159+ if (b >= H.n) { return; }
160+ var F = S3(offF() + 3u * b);
161+ if (Q.gather == 1u) {
162+ F = vec3<f32>(0.0);
163+ for (var q = I[H.iFr + b]; q < I[H.iFr + b + 1u]; q++) {
164+ F += bitcast<f32>(I[H.iFe + 2u * q + 1u]) * E3(I[H.iFe + 2u * q]);
165+ }
166+ if (H.npulse > 0u) { F += pulseForce(b, f32(C[0]) * H.Delta + Q.tOff); }
167+ setS3(offF() + 3u * b, F);
168+ }
169+ let v = V(b) + (Q.kick * P[H.oInvm + b]) * F;
170+ setS3(offV() + 3u * b, v);
171+ setS3(3u * b, U(b) + Q.drift * v);
172+}
173+
174+// ---- sample boundary ----------------------------------------------------------------
175+
176+// frame copy (if this sample is a frame), T = v, kinetic energy before damping
177+@compute @workgroup_size(WG)
178+fn pre(@builtin(global_invocation_id) gid: vec3<u32>) {
179+ let b = gid.x;
180+ if (b >= H.n) { return; }
181+ let sample = C[0];
182+ let v = V(b);
183+ if (sample % H.frameStride == 0u) {
184+ let first = ((H.sample0 + H.frameStride - 1u) / H.frameStride) * H.frameStride;
185+ let o = H.oFrames + ((sample - first) / H.frameStride) * H.frameSize;
186+ let u = U(b);
187+ O[o + 3u * b] = u.x;
188+ O[o + 3u * b + 1u] = u.y;
189+ O[o + 3u * b + 2u] = u.z;
190+ O[o + 3u * H.n + 3u * b] = v.x;
191+ O[o + 3u * H.n + 3u * b + 1u] = v.y;
192+ O[o + 3u * H.n + 3u * b + 2u] = v.z;
193+ O[o + 6u * H.n + b] = S[offDiss() + b];
194+ }
195+ setS3(offT() + 3u * b, v);
196+ S[offKE() + b] = 0.5 * P[H.oMass + b] * dot(v, v);
197+}
198+
199+// T = v, for each further piece of a damping step split into pieces
200+@compute @workgroup_size(WG)
201+fn resetT(@builtin(global_invocation_id) gid: vec3<u32>) {
202+ if (gid.x < H.n) { setS3(offT() + 3u * gid.x, V(gid.x)); }
203+}
204+
205+// one Taylor term of exp(-tau A): the damping forces of the elements from T
206+@compute @workgroup_size(WG)
207+fn dampElements(@builtin(global_invocation_id) gid: vec3<u32>) {
208+ var q = gid.x;
209+ if (q < H.nd) {
210+ let bi = I[H.iD + 2u * q];
211+ let bj = I[H.iD + 2u * q + 1u];
212+ let o = H.oD + 8u * q;
213+ let rw = S3(offT() + 3u * bi) - S3(offT() + 3u * bj);
214+ var f = P[o + 3u] * rw;
215+ if (P[o + 4u] != 0.0) {
216+ let dir = normalize(P3(o) + U(bj) - U(bi));
217+ f += (P[o + 4u] * dot(dir, rw)) * dir;
218+ }
219+ setE3(3u * q, f);
220+ return;
221+ }
222+ q -= H.nd;
223+ if (q < H.nb) {
224+ let cb = P[H.oB + 8u * q + 4u];
225+ let w = S3(offT() + 3u * I[H.iB + 3u * q]) - 2.0 * S3(offT() + 3u * I[H.iB + 3u * q + 1u]) + S3(offT() + 3u * I[H.iB + 3u * q + 2u]);
226+ setE3(H.eDb + 3u * q, cb * w);
227+ return;
228+ }
229+ q -= H.nb;
230+ if (q < H.nh) {
231+ let ch = P[H.oH + 16u * q + 10u];
232+ if (ch == 0.0) { return; } // no damping entries reference these slots
233+ let hg = hingeAt(q);
234+ var s = 0.0;
235+ for (var p = 0u; p < 4u; p++) { s += dot(hg.g[p], S3(offT() + 3u * I[H.iH + 4u * q + p])); }
236+ for (var p = 0u; p < 4u; p++) { setE3(H.eDh + 12u * q + 3u * p, (ch * s) * hg.g[p]); }
237+ }
238+}
239+
240+// A T per bead, then T <- c A T and v <- v + T (c = -tau/k)
241+@compute @workgroup_size(WG)
242+fn dampBeads(@builtin(global_invocation_id) gid: vec3<u32>) {
243+ let b = gid.x;
244+ if (b >= H.n) { return; }
245+ var f = P[H.oFric + b] * S3(offT() + 3u * b);
246+ for (var q = I[H.iDr + b]; q < I[H.iDr + b + 1u]; q++) {
247+ f += bitcast<f32>(I[H.iDe + 2u * q + 1u]) * E3(I[H.iDe + 2u * q]);
248+ }
249+ let a = P[H.oInvm + b] * f;
250+ if (Q.saveAv == 1u) { setS3(offAv() + 3u * b, a); }
251+ let t = Q.c * a;
252+ setS3(offT() + 3u * b, t);
253+ let v = V(b) + t;
254+ setS3(offV() + 3u * b, v);
255+ if (Q.last == 1u) {
256+ S[offDiss() + b] += S[offKE() + b] - 0.5 * P[H.oMass + b] * dot(v, v);
257+ }
258+}
259+
260+// the listener's sample; advances the sample index
261+@compute @workgroup_size(1)
262+fn record() {
263+ let dir = vec3<f32>(H.dirx, H.diry, H.dirz);
264+ var val = 0.0;
265+ for (var q = 0u; q < H.nl; q++) {
266+ let b = I[H.iL + q];
267+ if (H.quantity == 0u) { val += dot(U(b), dir); }
268+ else if (H.quantity == 1u) { val += dot(S3(offT() + 3u * b), dir); } // T holds v at the boundary
269+ else { val += dot(P[H.oInvm + b] * S3(offF() + 3u * b) - S3(offAv() + 3u * b), dir); }
270+ }
271+ O[C[0] - H.sample0] = val / f32(H.nl);
272+ C[0] = C[0] + 1u;
273+}
274+
275+// velocity impulses at the start of a chunk: [bead, vx, vy, vz] per entry
276+@compute @workgroup_size(WG)
277+fn impulse(@builtin(global_invocation_id) gid: vec3<u32>) {
278+ let q = gid.x;
279+ if (q >= H.nImp) { return; }
280+ let o = H.iImp + 4u * (H.impOff + q);
281+ let b = I[o];
282+ if (P[H.oInvm + b] > 0.0) {
283+ setS3(offV() + 3u * b, V(b) + vec3<f32>(bitcast<f32>(I[o + 1u]), bitcast<f32>(I[o + 2u]), bitcast<f32>(I[o + 3u])));
284+ }
285+}
286+
287+// rest angles of the hinges (u = 0)
288+@compute @workgroup_size(WG)
289+fn initHinges(@builtin(global_invocation_id) gid: vec3<u32>) {
290+ if (gid.x < H.nh) { S[offTh0() + gid.x] = hingeAt(gid.x).theta; }
291+}
src/sim/parse.tsmodified+1−1View file
@@ -277,7 +277,7 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
277277 safety: scalar(f.get("safety"), "sys.safety", 0.5),
278278 method,
279279 frameRate: scalar(f.get("frameRate"), "sys.frameRate", 500),
280- modes: scalar(f.get("modes"), "sys.modes", 1) !== 0,
280+ modes: f.has("modes") ? scalar(f.get("modes"), "sys.modes", 1) !== 0 : null,
281281 };
282282 if (!(sys.duration > 0)) throw new SysError("sys.duration must be positive");
283283 if (!(sys.safety > 0 && sys.safety <= 1)) throw new SysError("sys.safety must be in (0, 1]");
src/sim/prepare.tsadded+160−0View file
@@ -0,0 +1,160 @@
1+// Pieces shared by the CPU and GPU solvers: the static forces and energy in
2+// double precision, the checks of the rest state, the choice of time step,
3+// and the projection of the motion onto linear modes.
4+
5+import { angleDiff, hingeAngle } from "./hinge";
6+import { estimateOmegaMax, stabilityLimit } from "./solver";
7+import type { BeadSystem, Modes } from "./types";
8+
9+// Spring, bending and hinge forces at positions x (added into F, which is
10+// cleared first); returns the potential energy.
11+export function staticForces(sys: BeadSystem, x: Float64Array, F: Float64Array): number {
12+ F.fill(0);
13+ let pe = 0;
14+ for (let s = 0; s < sys.si.length; s++) {
15+ const i3 = 3 * sys.si[s];
16+ const j3 = 3 * sys.sj[s];
17+ const dx = x[j3] - x[i3];
18+ const dy = x[j3 + 1] - x[i3 + 1];
19+ const dz = x[j3 + 2] - x[i3 + 2];
20+ const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
21+ const e = L - sys.sL0[s];
22+ const e2 = e * e;
23+ pe += 0.5 * sys.sk[s] * e2 + 0.25 * sys.sk3[s] * e2 * e2;
24+ const g = (sys.sk[s] * e + sys.sk3[s] * e2 * e) / L;
25+ F[i3] += g * dx;
26+ F[i3 + 1] += g * dy;
27+ F[i3 + 2] += g * dz;
28+ F[j3] -= g * dx;
29+ F[j3 + 1] -= g * dy;
30+ F[j3 + 2] -= g * dz;
31+ }
32+ for (let b = 0; b < sys.bi.length; b++) {
33+ const i3 = 3 * sys.bi[b];
34+ const j3 = 3 * sys.bj[b];
35+ const k3 = 3 * sys.bk[b];
36+ for (let c = 0; c < 3; c++) {
37+ const d = x[i3 + c] - 2 * x[j3 + c] + x[k3 + c];
38+ pe += 0.5 * sys.bkb[b] * d * d;
39+ F[i3 + c] -= sys.bkb[b] * d;
40+ F[j3 + c] += 2 * sys.bkb[b] * d;
41+ F[k3 + c] -= sys.bkb[b] * d;
42+ }
43+ }
44+ const g = new Float64Array(12);
45+ for (let h = 0; h < sys.hi.length; h++) {
46+ const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
47+ const d = angleDiff(hingeAngle(x, beads[0], beads[1], beads[2], beads[3], g), sys.hth0[h]);
48+ pe += 0.5 * sys.hkh[h] * d * d;
49+ for (let p = 0; p < 4; p++) for (let c = 0; c < 3; c++) F[3 * beads[p] + c] -= sys.hkh[h] * d * g[3 * p + c];
50+ }
51+ return pe;
52+}
53+
54+export interface Prepared {
55+ peRest: number;
56+ residualForce: number;
57+ maxTension: number;
58+ omegaMax: number;
59+ substeps: number;
60+ h: number;
61+ nSamples: number;
62+ warnings: string[];
63+}
64+
65+export function prepare(sys: BeadSystem): Prepared {
66+ const { n, pinned, pos0 } = sys;
67+ const warnings: string[] = [];
68+ const F = new Float64Array(3 * n);
69+ const peRest = staticForces(sys, pos0, F);
70+ let residualForce = 0;
71+ for (let i = 0; i < n; i++) {
72+ if (!pinned[i]) residualForce = Math.max(residualForce, Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]));
73+ }
74+ let maxTension = 0;
75+ let forceScale = 0; // the larger of the rest tension and the force of a 1e-3 strain on any element
76+ for (let s = 0; s < sys.si.length; s++) {
77+ const i3 = 3 * sys.si[s];
78+ const j3 = 3 * sys.sj[s];
79+ const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
80+ const e = L - sys.sL0[s];
81+ maxTension = Math.max(maxTension, Math.abs(sys.sk[s] * e + sys.sk3[s] * e * e * e));
82+ forceScale = Math.max(forceScale, 1e-3 * sys.sk[s] * sys.sL0[s]);
83+ }
84+ forceScale = Math.max(forceScale, maxTension);
85+ for (let b = 0; b < sys.bi.length; b++) {
86+ const i3 = 3 * sys.bi[b];
87+ const j3 = 3 * sys.bj[b];
88+ const len = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
89+ forceScale = Math.max(forceScale, 1e-3 * sys.bkb[b] * len);
90+ }
91+ const g = new Float64Array(12);
92+ for (let h = 0; h < sys.hi.length; h++) {
93+ hingeAngle(pos0, sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h], g);
94+ forceScale = Math.max(forceScale, 1e-3 * sys.hkh[h] * Math.max(...g.map(Math.abs)));
95+ }
96+ if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
97+ warnings.push(
98+ `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
99+ `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
100+ );
101+ }
102+
103+ const omegaMax = estimateOmegaMax(sys);
104+ const fs = sys.sampleRate;
105+ const hLimit = stabilityLimit(sys.method) / omegaMax;
106+ let substeps: number;
107+ if (sys.dt > 0) {
108+ substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
109+ if (1 / (fs * substeps) > hLimit) {
110+ warnings.push(`dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`);
111+ }
112+ } else {
113+ substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
114+ }
115+ return {
116+ peRest,
117+ residualForce,
118+ maxTension,
119+ omegaMax,
120+ substeps,
121+ h: 1 / (fs * substeps),
122+ nSamples: Math.max(1, Math.round(sys.duration * fs)),
123+ warnings,
124+ };
125+}
126+
127+// Energy of each linear mode, E_k = (qdot_k^2 + omega_k^2 q_k^2)/2, from the
128+// displacements and velocities.
129+export class ModeTracker {
130+ private bz: Float64Array;
131+ private bv: Float64Array;
132+ private om2: Float64Array;
133+ constructor(private modes: Modes, private mass: Float64Array, private pos0: Float64Array) {
134+ this.bz = new Float64Array(modes.dof.length);
135+ this.bv = new Float64Array(modes.dof.length);
136+ this.om2 = Float64Array.from(modes.freqs, f => (2 * Math.PI * f) ** 2);
137+ }
138+ // x: positions (or displacements if relative is true); writes nModes energies into out at offset
139+ record(x: ArrayLike<number>, v: ArrayLike<number>, out: Float32Array, offset: number, relative = false) {
140+ const { dof, psi } = this.modes;
141+ const nf = dof.length;
142+ for (let r = 0; r < nf; r++) {
143+ const q = dof[r];
144+ const sm = Math.sqrt(this.mass[(q / 3) | 0]);
145+ this.bz[r] = sm * (relative ? x[q] : x[q] - this.pos0[q]);
146+ this.bv[r] = sm * v[q];
147+ }
148+ const nModes = this.om2.length;
149+ for (let k = 0; k < nModes; k++) {
150+ let q = 0;
151+ let qd = 0;
152+ const off = k * nf;
153+ for (let r = 0; r < nf; r++) {
154+ q += psi[off + r] * this.bz[r];
155+ qd += psi[off + r] * this.bv[r];
156+ }
157+ out[offset + k] = 0.5 * (qd * qd + this.om2[k] * q * q);
158+ }
159+ }
160+}
src/sim/solver.tsmodified+5−56View file
@@ -10,6 +10,7 @@
1010 // part is second order (Verlet) or fourth order (Yoshida) in the step h.
1111
1212 import { angleDiff, hingeAngle } from "./hinge";
13+import { prepare } from "./prepare";
1314 import type { BeadSystem, Method, Modes, SimResult } from "./types";
1415
1516 const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
@@ -168,64 +169,12 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
168169 }
169170 };
170171
171- // ---- rest state checks ---------------------------------------------------
172+ // ---- rest state checks and time step (shared with the GPU solver) ----------
172173 let pe = internalForces();
173- const peRest = pe; // energy stored by pre-tension; energies are reported relative to rest
174- let residualForce = 0;
175- for (let i = 0; i < n; i++) {
176- if (pinned[i]) continue;
177- const f = Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]);
178- if (f > residualForce) residualForce = f;
179- }
180- let maxTension = 0;
181- for (let s = 0; s < ns; s++) {
182- const i3 = 3 * si[s];
183- const j3 = 3 * sj[s];
184- const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
185- const e = L - sL0[s];
186- maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
187- }
188- // force scale: the larger of the rest tension and the force of a 1e-3 strain
189- // on any spring or bending triplet
190- let forceScale = maxTension;
191- for (let s = 0; s < ns; s++) forceScale = Math.max(forceScale, 1e-3 * sk[s] * sL0[s]);
192- for (let b = 0; b < nb; b++) {
193- const i3 = 3 * bi[b];
194- const j3 = 3 * bj[b];
195- const len = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
196- forceScale = Math.max(forceScale, 1e-3 * bkb[b] * len);
197- }
198- {
199- const g = new Float64Array(12);
200- for (let h = 0; h < nhg; h++) {
201- hingeAngle(x, hi[h], hj[h], hk[h], hl[h], g);
202- forceScale = Math.max(forceScale, 1e-3 * hkh[h] * Math.max(...g.map(Math.abs)));
203- }
204- }
205- if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
206- warnings.push(
207- `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
208- `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
209- );
210- }
211-
212- // ---- time step -----------------------------------------------------------
213- const omegaMax = estimateOmegaMax(sys);
174+ const prep = prepare(sys);
175+ const { peRest, residualForce, maxTension, omegaMax, substeps, h, nSamples } = prep;
176+ warnings.push(...prep.warnings);
214177 const fs = sys.sampleRate;
215- const hLimit = stabilityLimit(sys.method) / omegaMax;
216- let substeps: number;
217- if (sys.dt > 0) {
218- substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
219- if (1 / (fs * substeps) > hLimit) {
220- warnings.push(
221- `dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`
222- );
223- }
224- } else {
225- substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
226- }
227- const h = 1 / (fs * substeps);
228- const nSamples = Math.max(1, Math.round(sys.duration * fs));
229178 const nSteps = nSamples * substeps;
230179 const cs = stepCoefficients(sys.method);
231180
src/sim/types.tsmodified+1−1View file
@@ -60,7 +60,7 @@ export interface BeadSystem {
6060 safety: number; // fraction of the stability limit used when dt is automatic
6161 method: Method;
6262 frameRate: number; // Hz, snapshots for animation and mode energies
63- modes: boolean; // compute linear transverse modes and track their energies
63+ modes: boolean | null; // compute linear modes and track their energies; null: as the page's toggle says
6464 }
6565
6666 export interface SimProgress {
src/worker.tsmodified+19−8View file
@@ -3,7 +3,8 @@
33 import { growthRate, relDiff } from "./sim/analysis";
44 import { computeModes, isPlanar, modalDamping } from "./sim/modes";
55 import { numblEig, runScript } from "./sim/runScript";
6-import { simulate } from "./sim/solver";
6+import { simulateGPU } from "./sim/gpu";
7+import { simulate, type SimOptions } from "./sim/solver";
78 import { runTwin } from "./sim/twin";
89 import type { BeadSystem, SimResult } from "./sim/types";
910
@@ -25,6 +26,7 @@ export interface RunPayload {
2526 nBends: number;
2627 nHinges: number;
2728 listenQuantity: string;
29+ engine: string;
2830 listenBeads: Int32Array;
2931 hitBeads: Int32Array;
3032 method: string;
@@ -49,6 +51,8 @@ export interface RunRequest {
4951 dtCheck: boolean;
5052 linear: boolean;
5153 twin: boolean;
54+ modes: boolean; // the page's "linear modes" toggle (sys.modes in the script overrides it)
55+ gpu: boolean;
5256 }
5357
5458 export type FromWorker =
@@ -63,8 +67,10 @@ const post = (m: FromWorker, transfer: Transferable[] = []) => (self as unknown
6367 const windowsOf = (duration: number) =>
6468 [0.01, 0.03, 0.1, 0.3, 1, duration].filter((t, i, a) => t <= duration && a.indexOf(t) === i);
6569
66-self.onmessage = (ev: MessageEvent<RunRequest>) => {
67- const { script, dtCheck, linear, twin } = ev.data;
70+self.onmessage = async (ev: MessageEvent<RunRequest>) => {
71+ const { script, dtCheck, linear, twin: twinWanted, modes: modesToggle, gpu } = ev.data;
72+ // the GPU path for the main run and the checks; the CPU otherwise
73+ const sim = (s: BeadSystem, o: SimOptions) => (gpu ? simulateGPU(s, o) : Promise.resolve(simulate(s, o)));
6874 try {
6975 post({ type: "stage", text: "Running script" });
7076 const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text }));
@@ -72,7 +78,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
7278
7379 let modes = null;
7480 const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
75- if (sys.modes) {
81+ if (sys.modes ?? modesToggle) {
7682 const planar = isPlanar(sys);
7783 if (planar && nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS} for a flat net).`);
7884 else if (!planar && 3 * nFree > MAX_FULL_DOF) {
@@ -96,11 +102,15 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
96102 // positions at every frame, kept for the perturbed twin
97103 const saved: Float64Array[] = [];
98104 const savedBytes = Math.round(sys.duration * sys.frameRate + 1) * 3 * sys.n * 8;
105+ if (twinWanted && gpu) {
106+ notes.push("Perturbed twin skipped: a perturbation of one part in 1e9 is below single precision; turn off WebGPU to run it.");
107+ }
108+ const twin = twinWanted && !gpu;
99109 const keep = twin && savedBytes <= MAX_TWIN_BYTES;
100110 if (twin && !keep) notes.push(`Perturbed twin skipped: it would store ${(savedBytes / 1e6).toFixed(0)} MB of positions (limit ${MAX_TWIN_BYTES / 1e6} MB). Lower sys.frameRate or sys.duration.`);
101111
102- post({ type: "stage", text: "Simulating" });
103- const result = simulate(sys, {
112+ post({ type: "stage", text: gpu ? "Simulating (WebGPU)" : "Simulating" });
113+ const result = await sim(sys, {
104114 modes,
105115 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating" }),
106116 onFrame: keep ? (_f, _t, x) => saved.push(x.slice()) : undefined,
@@ -110,7 +120,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
110120 if (dtCheck && !result.aborted) {
111121 post({ type: "stage", text: "Simulating again with dt/2" });
112122 const half: BeadSystem = { ...sys, dt: result.dt / 2 };
113- const r2 = simulate(half, {
123+ const r2 = await sim(half, {
114124 maxAnimBytes: 0,
115125 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating with dt/2" }),
116126 });
@@ -132,7 +142,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
132142 force: [h.force[0] * s, h.force[1] * s, h.force[2] * s],
133143 })),
134144 };
135- const r3 = simulate(small, {
145+ const r3 = await sim(small, {
136146 maxAnimBytes: 0,
137147 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Small-strike reference" }),
138148 });
@@ -165,6 +175,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
165175 nBends: sys.bi.length,
166176 nHinges: sys.hi.length,
167177 listenQuantity: sys.listenQuantity,
178+ engine: gpu ? "WebGPU, single precision" : "CPU, double precision",
168179 listenBeads: sys.listenBeads,
169180 hitBeads,
170181 method: sys.method,
tsconfig.jsonmodified+1−1View file
@@ -9,7 +9,7 @@
99 "noUnusedLocals": true,
1010 "noUnusedParameters": true,
1111 "skipLibCheck": true,
12- "types": ["vite/client", "node"],
12+ "types": ["vite/client", "node", "@webgpu/types"],
1313 "noEmit": true,
1414 "allowImportingTsExtensions": true
1515 },