Add WebGPU solver with a toggle; faster CPU solver
- Time step from power iteration on the linearized stiffness (the Gershgorin bound overestimated omega_max by 1.4-1.5x on the cymbal)
- CPU: recompute only velocity-dependent forces at each audio sample
- WebGPU solver (src/sim/gpu.ts): one small kernel per phase, many samples per command buffer, f32 with displacement-based state; falls back to the CPU on failure
- WebGPU toggle next to Run (on for drum and cymbal, off for string)
- bench.html compares both solvers on every preset
14 changed files+1529−139
README.mdmodified+15−5View file
@@ -2,7 +2,7 @@
22
33 Live: https://magland.github.io/toy-instruments/ (deployed from `main` by `.github/workflows/pages.yml`).
44
5-A browser gallery of toy physical-model instruments. Pick an instrument, set its geometry, properties and excitation, run the simulation, and get the pickup waveform, a spectrogram, a playable/downloadable `.wav`, and the bead displacements over time. Static site (Vite + TypeScript); the solver runs in a Web Worker.
5+A browser gallery of toy physical-model instruments. Pick an instrument, set its geometry, properties and excitation, run the simulation, and get the pickup waveform, a spectrogram, a playable/downloadable `.wav`, and the bead displacements over time. Static site (Vite + TypeScript); the solver runs in a Web Worker, on the CPU or with WebGPU (see Solvers).
66
77 ```
88 npm install
@@ -18,7 +18,7 @@ npm run build # static site in dist/
1818 - Excitation: pluck (triangle between the nearest supports, released at rest), strike (Gaussian initial velocity), hammer (raised-cosine force pulse). Direction is an angle in the y-z plane.
1919 - Damping from two decay times: `sigma(omega) = sigma0 + sigma1 omega^2`, fitted to T60 at two frequencies for transverse modes. `sigma0` is mass-proportional; `sigma1` is a dashpot on each spring's relative velocity.
2020 - Audio: pickup acceleration (forces / mass, at the full-step velocity), from one bead or summed over all beads, along y, z or x. The `.wav` is normalized to -1 dBFS.
21-- Solver: velocity Verlet (explicit, symplectic without damping). `dt = cfl * 2 / omega_max` (Gershgorin bound with axial `k` and transverse `T/l`), rounded down to a whole number of steps per audio sample.
21+- Solver: velocity Verlet (explicit, symplectic without damping). `dt = cfl * 2 / omega_max`, rounded down to a whole number of steps per audio sample; `omega_max` is the highest frequency of the system linearized at rest (power iteration, axial `k` and transverse `T/l` per spring, bending, plus a mallet's contact stiffness). See Solvers.
2222 - 3D model view (three.js): Setup shows the beads at rest with supports (wall = clamped, translucent plane + ring = sliding), pinned beads (cubes), excitation (pluck shape, strike velocity arrows, hammer force arrows) and pickup (bead + axis arrow). Replay animates the saved frames ("Replay frames" per second, default 2000) with adjustable exaggeration, on a timeline shared with the 2D side/top views and space-time image.
2323 - Checks shown per run: energy balance (energy + dissipated - work in - initial) relative to peak energy, largest strain. Linear transverse and longitudinal modes (Sturm bisection on the tridiagonal stiffness) are marked on the spectrum.
2424
@@ -31,7 +31,7 @@ String tests (`test/string.test.ts`): modes vs the discrete-string formula; smal
3131 - Excitation: felt mallet (point mass, contact force K delta^p on a Gaussian patch, integrated with the beads; contact time and rebounds come out of the simulation), strike, force pulse, press (initial dent).
3232 - Damping as for the string (dashpots calibrated to the springs' transverse stiffness f / a).
3333 - Extra views: top-down displacement map and cross-section along the diameter through the strike point; mallet contact force.
34-- Presets: a fast tom with in-plane stiffness softened ~70x (about 8x slower than real time for 571 beads in Chromium) and a realistic Mylar one (~6x more steps).
34+- Presets: a fast tom with in-plane stiffness softened ~70x (1.5 s of sound in 9 s on the CPU and 5.6 s with WebGPU, 571 beads, Chromium on an Intel Iris Xe laptop) and a realistic Mylar one (~5x more steps).
3535 - Tests (`test/drum.test.ts`): every bead has six springs; a small centre press rings within 1.3% of the ideal membrane fundamental 2.405 c / (2 pi R); the mallet bounces off and energy balances.
3636
3737 ## Cymbal
@@ -43,11 +43,21 @@ String tests (`test/string.test.ts`): modes vs the discrete-string formula; smal
4343 - Excitation: drumstick (Hertz, K delta^1.5) or felt mallet as a 3D point mass with contact along the surface normal, regularized Coulomb friction and a contact reach; hit angle tilts the stroke toward the rim. Also strike and force pulse.
4444 - Sound: far-field microphone via the Rayleigh integral, p = rho / (2 pi) sum A_i a_n,i(t - R_i / c) / R_i (baffled-plate assumption: no diffraction around the rim).
4545 - Tests (`test/cymbal.test.ts`, `test/hinge.test.ts`): hinge gradient vs finite differences; bending energy matches the plate for cylinder, sphere, saddle and twist to 0.1%; a free flat plate's (2,0) mode within 3.4% of Leissa's value on a 9-ring lattice; energy balance with a glancing stick, friction, felts and damping.
46-- Limits: about 45x slower than real time for ~370 beads in Chromium; bending is resolved only to ~1.3 kHz at the default resolution, so it sounds like a small gong, not a crash. A single-bead stick contact chatters on this lattice, so presets spread the contact over ~1.5 bead spacings.
46+- Limits: about 30x slower than real time on the CPU and 14x with WebGPU for ~370 beads (Chromium, Intel Iris Xe laptop); bending is resolved only to ~1.3 kHz at the default resolution, so it sounds like a small gong, not a crash. A single-bead stick contact chatters on this lattice, so presets spread the contact over ~1.5 bead spacings.
47+
48+## Solvers
49+
50+Both solvers integrate the same scheme (`src/sim/solver.ts` on the CPU in double precision, `src/sim/gpu.ts` with WebGPU in single precision); the WebGPU toggle next to Run picks one, and a failed GPU run falls back to the CPU with the reason in the run summary. The toggle is on by default for the drum and cymbal and off for the string, whose few dozen beads run faster on the CPU.
51+
52+- Time step: `omega_max` comes from power iteration on the linearized stiffness rather than a Gershgorin bound, which overestimated it by 1.4 to 1.5x on the cymbal (so steps were that much shorter than needed).
53+- CPU: at each audio sample the solver recomputes only the velocity-dependent forces (damping, friction) at the full-step velocity; spring and hinge elasticity depend on positions, which have not changed. Together with the time step, this made the cymbal presets 1.4 to 1.6x faster.
54+- WebGPU: each phase of a time step (kick and drift, hinge angles and contact, hinge moments, bead forces, kick with rigid-motion removal) is a small kernel; a command buffer holds the dispatches for many samples, and step and sample counters live on the GPU. Forces are gathered per bead from fixed-width adjacency lists, so there are no atomics. To keep small vibrations from vanishing in f32, the state is the displacement from rest, and spring extensions and hinge angles are computed from precomputed rest differences. Sums over beads (rigid motion, energy, audio) go through per-workgroup partials. On the first 50 ms of every preset, the audio agrees with the CPU to within 0.5% RMS (0.1% for most), and the energy-balance check and largest strain match.
55+- Speed, WebGPU vs CPU on an Intel Iris Xe laptop (Chromium): the default cymbal (366 beads) about 2x, a 20-ring cymbal (1444 beads) 4x, 30 rings (3188 beads) 6x; strings (3 to 50 beads) 4 to 80x slower in 50 ms test runs, where the GPU's start-up (about half a second) dominates. A single persistent workgroup running whole time steps per dispatch was tried first; on this GPU its large kernel made every barrier cost several microseconds, and it used only one of the GPU's subslices.
56+- `bench.html` (dev server only, `npm run dev` then `/bench.html?dur=0.05&only=cymbal&rings=20`) runs every preset on both solvers and prints the differences and timings.
4757
4858 ## Layout
4959
50-- `src/sim/` generic 3D bead-spring system, solver, worker.
60+- `src/sim/` generic 3D bead-spring system, CPU solver (`solver.ts`), WebGPU solver (`gpu.ts`), worker.
5161 - `src/instruments/<id>/` model (parameters to system), presets, page. Register new instruments in `src/instruments/registry.ts`.
5262 - `src/ui/instrumentPage.ts` the shared page (controls, 3D model, setup, run, audio, plots); an instrument supplies setup rows and extra cards.
5363 - `src/ui/` reusable controls form, line plot, spectrogram, chain displacement viewer, 3D scene, shared playback timeline, FFT, WAV.
bench.htmladded+8−0View file
@@ -0,0 +1,8 @@
1+<!doctype html>
2+<html lang="en">
3+<head><meta charset="utf-8"><title>Solver bench</title></head>
4+<body>
5+<pre id="out">running…</pre>
6+<script type="module" src="/bench/gpu.ts"></script>
7+</body>
8+</html>
bench/gpu.tsadded+90−0View file
@@ -0,0 +1,90 @@
1+// CPU vs WebGPU solver on every preset (dev only: open /bench.html).
2+// Query: ?dur=0.05&only=Flat&cpu=0&rings=20 (rings: cymbal mesh override)
3+import { simulate } from '../src/sim/solver';
4+import { simulateGPU } from '../src/sim/gpu';
5+import type { BeadSystem, Excitation, SimOptions, SimResult } from '../src/sim/types';
6+import { CYMBAL_PRESETS } from '../src/instruments/bead-cymbal/presets';
7+import { buildCymbal, buildCymbalExcitation, buildMic, deriveCymbal } from '../src/instruments/bead-cymbal/model';
8+import { DRUM_PRESETS } from '../src/instruments/bead-drum/presets';
9+import { buildDrum, buildDrumExcitation, buildDrumPickup, deriveDrum } from '../src/instruments/bead-drum/model';
10+import { PRESETS as STRING_PRESETS } from '../src/instruments/bead-string/presets';
11+import { buildExcitation, buildPickup, buildSystem, derive } from '../src/instruments/bead-string/model';
12+
13+const qs = new URLSearchParams(location.search);
14+const dur = Number(qs.get('dur') ?? 0.05);
15+const only = qs.get('only');
16+const runCpu = qs.get('cpu') !== '0';
17+const rings = qs.get('rings') ? Number(qs.get('rings')) : undefined;
18+const out = document.getElementById('out')!;
19+const lines: string[] = [];
20+const log = (s: string) => {
21+ lines.push(s);
22+ out.textContent = lines.join('\n');
23+};
24+
25+interface Case { name: string; sys: BeadSystem; exc: Excitation; opt: SimOptions }
26+const cases: Case[] = [];
27+for (const pr of CYMBAL_PRESETS) {
28+ const p = rings ? { ...pr.params, rings } : pr.params;
29+ const d = deriveCymbal(p);
30+ const mic = buildMic(p, d);
31+ cases.push({ name: `cymbal: ${pr.name}`, sys: buildCymbal(d), exc: buildCymbalExcitation(p, d), opt: { duration: dur, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay } });
32+}
33+for (const pr of DRUM_PRESETS) {
34+ const p = pr.params;
35+ const d = deriveDrum(p);
36+ cases.push({ name: `drum: ${pr.name}`, sys: buildDrum(d), exc: buildDrumExcitation(p, d), opt: { duration: dur, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: buildDrumPickup(p, d) } });
37+}
38+for (const pr of STRING_PRESETS) {
39+ const p = pr.params;
40+ const d = derive(p);
41+ cases.push({ name: `string: ${pr.name}`, sys: buildSystem(p, d), exc: buildExcitation(p), opt: { duration: dur, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: buildPickup(p) } });
42+}
43+
44+const maxAbs = (a: ArrayLike<number>) => {
45+ let m = 0;
46+ for (let i = 0; i < a.length; i++) m = Math.max(m, Math.abs(a[i]));
47+ return m;
48+};
49+// relative RMS difference of b against a over [0, len)
50+function relDiff(a: ArrayLike<number>, b: ArrayLike<number>, len: number) {
51+ let d = 0, s = 0;
52+ for (let i = 0; i < len; i++) {
53+ d += (a[i] - b[i]) ** 2;
54+ s += a[i] ** 2;
55+ }
56+ return Math.sqrt(d / (s || 1));
57+}
58+const balErr = (r: SimResult) => maxAbs(r.balance) / (maxAbs(r.energy) || 1);
59+
60+(async () => {
61+ const results: unknown[] = [];
62+ for (const c of cases) {
63+ if (only && !c.name.includes(only)) continue;
64+ try {
65+ const g = await simulateGPU(c.sys, c.exc, c.opt);
66+ const row: Record<string, unknown> = { name: c.name, n: c.sys.n, spp: g.stepsPerSample, gpuWall: g.wallSeconds, gpuLoop: g.loopSeconds, usPerStep: (g.loopSeconds * 1e6) / (g.audio.length * g.stepsPerSample), gpuSlow: g.loopSeconds / dur, gpuBal: balErr(g), device: g.device };
67+ if (runCpu) {
68+ const r = simulate(c.sys, c.exc, c.opt);
69+ const L = r.audio.length;
70+ Object.assign(row, {
71+ cpuWall: r.wallSeconds,
72+ speedup: r.wallSeconds / g.wallSeconds,
73+ cpuBal: balErr(r),
74+ audio10ms: relDiff(r.audio, g.audio, Math.min(L, Math.round(0.01 * r.fs))),
75+ audioAll: relDiff(r.audio, g.audio, L),
76+ frames: relDiff(r.frames, g.frames, r.frames.length),
77+ energyEnd: g.energy[g.energy.length - 1] / (r.energy[r.energy.length - 1] || 1),
78+ stretch: [r.maxStretch, g.maxStretch],
79+ contact: r.contactForce ? relDiff(r.contactForce, g.contactForce!, L) : undefined,
80+ });
81+ }
82+ results.push(row);
83+ log(JSON.stringify(row, (_, v) => (typeof v === 'number' ? Number(v.toPrecision(4)) : v)));
84+ } catch (e) {
85+ log(`${c.name}: ERROR ${(e as Error).message}`);
86+ }
87+ }
88+ log('DONE');
89+ (window as unknown as { benchResults: unknown }).benchResults = results;
90+})();
package-lock.jsonmodified+1−0View file
@@ -12,6 +12,7 @@
1212 },
1313 "devDependencies": {
1414 "@types/three": "^0.180.0",
15+ "@webgpu/types": "^0.1.74",
1516 "typescript": "^5.9.0",
1617 "vite": "^7.1.0",
1718 "vitest": "^3.2.0"
package.jsonmodified+1−0View file
@@ -12,6 +12,7 @@
1212 },
1313 "devDependencies": {
1414 "@types/three": "^0.180.0",
15+ "@webgpu/types": "^0.1.74",
1516 "typescript": "^5.9.0",
1617 "vite": "^7.1.0",
1718 "vitest": "^3.2.0"
src/instruments/bead-string/page.tsmodified+1−0View file
@@ -110,6 +110,7 @@ export function mountBeadString(root: HTMLElement) {
110110 mountInstrumentPage<StringParams>(root, {
111111 title: 'String of beads',
112112 wavName: 'bead-string.wav',
113+ webgpuDefault: false, // a few dozen beads: the CPU is faster
113114 presets: PRESETS,
114115 sections: SECTIONS,
115116 setup(p) {
src/sim/gpu.tsadded+1098−0View file
@@ -0,0 +1,1098 @@
1+import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
2+import { timeStep, type Progress } from './solver';
3+
4+// The same scheme as simulate() in solver.ts (velocity Verlet, same time step,
5+// rigid-motion removal, contact, delayed pickup), run on the GPU in single
6+// precision. Each phase of a time step is its own small kernel, and a command
7+// buffer holds the dispatches for many samples; time and sample counters live
8+// on the GPU, so nothing returns to JavaScript until a chunk is done. (One
9+// persistent workgroup running whole time steps was tried first: its large
10+// kernel made every barrier cost several microseconds on an Intel iGPU, and it
11+// could use only one of the GPU's subslices.) Forces are gathered per bead from
12+// fixed-width adjacency lists, so no atomics are needed. Sums over beads go
13+// through per-workgroup partials, which each workgroup of the next kernel adds
14+// up itself (cheaper than a dispatch of its own). Three buffers only: the
15+// per-dispatch cost grows with the number bound. To keep small vibrations from vanishing in
16+// f32, the state is the displacement from rest, and spring extensions and
17+// hinge angles are built from precomputed rest differences.
18+
19+const WGB = 64; // invocations per workgroup
20+
21+interface Packed {
22+ wgsl: string;
23+ ro: Uint32Array; // parameters (f32 bits) and indices
24+ rw: Float32Array; // state, shared scalars, per-workgroup partial sums
25+ outLen: number; // audio, contact force, frame scalars, mallet positions, then frames
26+ framesAt: number;
27+ nSamples: number;
28+ nFrames: number;
29+ samplesPerFrame: number;
30+ dt: number;
31+ stepsPerSample: number;
32+ omegaMax: number;
33+ nBead: number; // workgroups over beads
34+ nHinge: number; // workgroups over hinges
35+ hasMallet: boolean;
36+ freeRigid: boolean;
37+ velForces: boolean;
38+}
39+
40+function f(x: number): string {
41+ if (!Number.isFinite(x)) throw new Error(`non-finite constant ${x}`);
42+ return x.toExponential(9);
43+}
44+
45+class Pack<T extends Float32Array | Uint32Array> {
46+ parts: ArrayLike<number>[] = [];
47+ len = 0;
48+ constructor(private Ctor: new (n: number) => T) {}
49+ add(a: ArrayLike<number>): number {
50+ const off = this.len;
51+ this.parts.push(a);
52+ this.len += a.length;
53+ return off;
54+ }
55+ build(): T {
56+ const out = new this.Ctor(Math.max(4, this.len));
57+ let o = 0;
58+ for (const p of this.parts) {
59+ out.set(p, o);
60+ o += p.length;
61+ }
62+ return out;
63+ }
64+}
65+
66+// rows of variable length -> slot-major table [slot][row] padded with `pad`
67+function padded<T>(rows: T[][], pad: T): { width: number; table: T[] } {
68+ const width = Math.max(1, ...rows.map((r) => r.length));
69+ const table: T[] = [];
70+ for (let q = 0; q < width; q++) for (const r of rows) table.push(q < r.length ? r[q] : pad);
71+ return { width, table };
72+}
73+
74+function pack(sys: BeadSystem, exc: Excitation, opt: SimOptions, sg: boolean): Packed {
75+ const { n, si, sj, sk, sL0, sc, mass, pinned, anchors, rest } = sys;
76+ const mallet = exc.mallet;
77+ const { omegaMax, dt, stepsPerSample } = timeStep(sys, mallet, opt.fs, opt.cfl);
78+ const nSamples = Math.round(opt.duration * opt.fs);
79+ const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
80+ const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
81+
82+ const prm = new Pack(Float32Array);
83+ const idx = new Pack(Uint32Array);
84+
85+ // beads
86+ const invMass = Float64Array.from(mass, (m, i) => (pinned[i] ? 0 : 1 / m));
87+ let M = 0;
88+ const cm = [0, 0, 0];
89+ for (let i = 0; i < n; i++) {
90+ M += mass[i];
91+ for (let c = 0; c < 3; c++) cm[c] += mass[i] * rest[3 * i + c];
92+ }
93+ const r0 = new Float64Array(3 * n);
94+ for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) r0[3 * i + c] = rest[3 * i + c] - cm[c] / M;
95+ const i0 = new Uint32Array(n);
96+ const fr = new Float64Array(n);
97+ if (opt.pickupDelay)
98+ for (let i = 0; i < n; i++) {
99+ i0[i] = Math.floor(opt.pickupDelay[i]);
100+ fr[i] = opt.pickupDelay[i] - i0[i];
101+ }
102+ const H = Math.max(...i0) + 2;
103+ const pw = mallet ? mallet.weights : new Float64Array(n);
104+ const patch: number[] = [];
105+ for (let i = 0; i < n; i++) if (pw[i] > 0) patch.push(i);
106+ let mEff = 1;
107+ if (mallet) {
108+ let inv = 0;
109+ for (const i of patch) inv += (pw[i] * pw[i]) / mass[i];
110+ mEff = 1 / (1 / mallet.mass + inv);
111+ }
112+
113+ const P_MASS = prm.add(mass);
114+ const P_IM = prm.add(invMass);
115+ const P_R0 = prm.add(r0);
116+ const P_W = prm.add(opt.pickup);
117+ const P_FR = prm.add(fr);
118+ const P_PW = prm.add(pw);
119+ const P_FS = prm.add(exc.force ? exc.force.shape : new Float64Array(3 * n));
120+
121+ // springs, gathered per bead: [k, c, rest difference (3), rest length, e0, weight]
122+ // e0 = rest length - natural length (pre-tension); the weight counts a spring
123+ // between two moving beads half at each end. Padding: a zero spring to a fixed point.
124+ const FIXED = 0xffffffff;
125+ const SLIDE = 0x80000000;
126+ const code = (node: number) => (node < n ? node : anchors[node - n].slideBead < 0 ? FIXED : SLIDE + anchors[node - n].slideBead);
127+ const nodePos = (node: number, c: number) => {
128+ if (node < n) return rest[3 * node + c];
129+ const a = anchors[node - n];
130+ return a.slideBead < 0 || c === 0 ? a.pos[c] : rest[3 * a.slideBead + c];
131+ };
132+ type SpringEntry = { other: number; vals: number[] };
133+ const byBead: SpringEntry[][] = Array.from({ length: n }, () => []);
134+ for (let s = 0; s < si.length; s++) {
135+ for (const [me, ot] of [[si[s], sj[s]], [sj[s], si[s]]]) {
136+ if (me >= n || pinned[me]) continue;
137+ const rd = [0, 1, 2].map((c) => nodePos(ot, c) - nodePos(me, c));
138+ const lr = Math.hypot(rd[0], rd[1], rd[2]);
139+ const moving = ot < n && !pinned[ot];
140+ byBead[me].push({ other: code(ot), vals: [sk[s], sc[s], rd[0], rd[1], rd[2], lr, lr - sL0[s], moving ? 0.5 : 1] });
141+ }
142+ }
143+ const sp = padded(byBead, { other: FIXED, vals: [0, 0, 1, 0, 0, 1, 0, 0] });
144+ const P_SF = prm.add(sp.table.flatMap((e) => e.vals));
145+ const I_SO = idx.add(sp.table.map((e) => e.other));
146+
147+ // bending: hinge nh is a dummy with zero angle, rate, moment and gradient
148+ const b = sys.bending;
149+ const nh = b ? b.w0.length : 0;
150+ let P_HF = 0, P_HC = 0, I_HN = 0, I_HCC = 0, I_BHE = 0;
151+ let MAXHC = 1, MAXBH = 1;
152+ if (b) {
153+ const hf: number[] = [];
154+ const hn: number[] = [];
155+ const d = (p: number, q: number) => [0, 1, 2].map((c) => rest[3 * p + c] - rest[3 * q + c]);
156+ for (let e = 0; e < nh; e++) {
157+ const [w0, w1, e0, e1] = [b.w0[e], b.w1[e], b.e0[e], b.e1[e]];
158+ if (Math.max(w0, w1, e0, e1) >= n) throw new Error('hinge on an anchor');
159+ hn.push(w0, w1, e0, e1);
160+ hf.push(b.sign[e], b.theta0[e], ...d(e1, e0), ...d(w0, e0), ...d(w1, e1));
161+ }
162+ // moment of hinge a = sum over its triangles of M[a][c] phi_c
163+ const rows: Map<number, number>[] = Array.from({ length: nh }, () => new Map());
164+ for (let t = 0; t < b.tri.length / 3; t++)
165+ for (let a = 0; a < 3; a++) {
166+ const ea = b.tri[3 * t + a];
167+ if (ea < 0) continue;
168+ for (let c = 0; c < 3; c++) {
169+ const ec = b.tri[3 * t + c];
170+ if (ec >= 0) rows[ea].set(ec, (rows[ea].get(ec) ?? 0) + b.M[9 * t + 3 * a + c]);
171+ }
172+ }
173+ const hc = padded(rows.map((r) => [...r]), [nh, 0]);
174+ MAXHC = hc.width;
175+ // per bead: the hinges it belongs to (hinge * 4 + slot)
176+ const bh: number[][] = Array.from({ length: n }, () => []);
177+ for (let e = 0; e < nh; e++) for (let k = 0; k < 4; k++) if (!pinned[hn[4 * e + k]]) bh[hn[4 * e + k]].push(4 * e + k);
178+ const bhp = padded(bh, 4 * nh);
179+ MAXBH = bhp.width;
180+ P_HF = prm.add(hf);
181+ P_HC = prm.add(hc.table.map(([, v]) => v));
182+ I_HN = idx.add(hn);
183+ I_HCC = idx.add(hc.table.map(([c]) => c));
184+ I_BHE = idx.add(bhp.table);
185+ }
186+ const I_PATCH = idx.add(patch);
187+ const I_I0 = idx.add(i0);
188+
189+ // state
190+ const S_U = 0;
191+ const S_V = S_U + 3 * n;
192+ const S_A = S_V + 3 * n;
193+ const S_G = S_A + 3 * n;
194+ const S_PHI = S_G + 12 * (nh + 1);
195+ const S_RATE = S_PHI + nh + 1;
196+ const S_MOM = S_RATE + nh + 1;
197+ const S_MEL = S_MOM + nh + 1;
198+ const S_HIST = S_MEL + nh + 1;
199+ const state = new Float32Array(S_HIST + n * H);
200+ for (let i = 0; i < n; i++)
201+ if (!pinned[i])
202+ for (let c = 0; c < 3; c++) {
203+ state[S_U + 3 * i + c] = exc.u0[3 * i + c];
204+ state[S_V + 3 * i + c] = exc.v0[3 * i + c];
205+ }
206+ const nBead = Math.ceil(n / WGB);
207+ const nHinge = Math.ceil(nh / WGB);
208+ // read-write buffer: the state, then shared scalars, then per-workgroup partial sums
209+ const SC = state.length;
210+ const Q_BF = SC + 64; // 17 per bead workgroup
211+ const Q_HD = Q_BF + 17 * nBead; // 1 per hinge workgroup (hinge damping power)
212+ const Q_KC = Q_HD + Math.max(1, nHinge); // 6 per bead workgroup
213+ const Q_REC = Q_KC + 6 * nBead; // 3 per bead or hinge workgroup: audio, energy, largest strain
214+ const Q_IRV = Q_REC + 3 * (nBead + nHinge); // 15 per bead workgroup
215+ const rw = new Float32Array(Q_IRV + 15 * nBead);
216+ rw.set(state);
217+ if (mallet) rw.set(mallet.vel0, SC + 3);
218+ // read-only buffer: parameters (as f32 bits), then indices
219+ const prmArr = prm.build();
220+ const idxArr = idx.build();
221+ const ro = new Uint32Array(prmArr.length + idxArr.length);
222+ ro.set(new Uint32Array(prmArr.buffer));
223+ ro.set(idxArr, prmArr.length);
224+ const IO = prmArr.length;
225+
226+ const O_AUD = 0;
227+ const O_CF = nSamples;
228+ const O_FSC = 2 * nSamples;
229+ const O_MP = O_FSC + 6 * nFrames;
230+ const OFR = O_MP + 3 * nFrames;
231+
232+ const force = exc.force;
233+ const sigma1 = b?.sigma1 ?? 0;
234+ const anyDash = sc.some((c) => c !== 0);
235+ const velForces = sys.sigma0 !== 0 || anyDash || sigma1 !== 0 || (mallet?.mu ?? 0) > 0;
236+ const u = (x: number) => `${x}u`;
237+ const consts = `
238+const N = ${u(n)};
239+const NH = ${u(nh)};
240+const NP = ${u(patch.length)};
241+const NBW = ${u(nBead)};
242+const NHW = ${u(nHinge)};
243+const MAXS = ${u(sp.width)};
244+const MAXBH = ${u(MAXBH)};
245+const MAXHC = ${u(MAXHC)};
246+const H = ${u(H)};
247+const SPF = ${u(samplesPerFrame)};
248+const NF = ${u(nFrames)};
249+const DT = ${f(dt)};
250+const S0 = ${f(sys.sigma0)};
251+const S1 = ${f(sigma1)};
252+const MTOT = ${f(M)};
253+const HAS_MALLET = ${!!mallet};
254+const FREE_RIGID = ${!!sys.freeRigid};
255+const HAS_PULSE = ${!!force};
256+const ANY_DASH = ${anyDash};
257+const MM = ${f(mallet?.mass ?? 1)};
258+const MK = ${f(mallet?.K ?? 0)};
259+const MP = ${f(mallet?.p ?? 1)};
260+const MU = ${f(mallet?.mu ?? 0)};
261+const REACH2 = ${f(mallet?.reach !== undefined ? mallet.reach * mallet.reach : -1)};
262+const MEFF = ${f(mEff)};
263+const NRM = vec3f(${(mallet?.normal ?? [0, 0, 0]).map(f).join(', ')});
264+const F0 = ${f(force?.F0 ?? 0)};
265+const TAU = ${f(force?.tau ?? 1)};
266+const P_MASS = ${u(P_MASS)}; const P_IM = ${u(P_IM)}; const P_R0 = ${u(P_R0)}; const P_W = ${u(P_W)};
267+const P_FR = ${u(P_FR)}; const P_PW = ${u(P_PW)}; const P_FS = ${u(P_FS)}; const P_SF = ${u(P_SF)};
268+const P_HF = ${u(P_HF)}; const P_HC = ${u(P_HC)};
269+const I_SO = ${u(I_SO + IO)}; const I_HN = ${u(I_HN + IO)}; const I_HCC = ${u(I_HCC + IO)}; const I_BHE = ${u(I_BHE + IO)};
270+const I_PATCH = ${u(I_PATCH + IO)}; const I_I0 = ${u(I_I0 + IO)};
271+const S_U = ${u(S_U)}; const S_V = ${u(S_V)}; const S_A = ${u(S_A)}; const S_G = ${u(S_G)};
272+const S_PHI = ${u(S_PHI)}; const S_RATE = ${u(S_RATE)}; const S_MOM = ${u(S_MOM)}; const S_MEL = ${u(S_MEL)};
273+const S_HIST = ${u(S_HIST)}; const SC = ${u(SC)};
274+const Q_BF = ${u(Q_BF)}; const Q_HD = ${u(Q_HD)}; const Q_KC = ${u(Q_KC)}; const Q_REC = ${u(Q_REC)}; const Q_IRV = ${u(Q_IRV)};
275+const O_AUD = ${u(O_AUD)}; const O_CF = ${u(O_CF)}; const O_FSC = ${u(O_FSC)}; const O_MP = ${u(O_MP)}; const OFR = ${u(OFR)};
276+`;
277+ return {
278+ wgsl: (sg ? 'enable subgroups;\n' : '') + consts + reducers(sg) + KERNELS.replaceAll('/*SG*/', sg ? ', @builtin(subgroup_invocation_id) sgiIn: u32, @builtin(subgroup_size) sgsIn: u32' : '').replaceAll('/*SGI*/', sg ? 'sgi = sgiIn; sgs = sgsIn;' : ''),
279+ ro,
280+ rw,
281+ outLen: OFR + 3 * n * nFrames,
282+ framesAt: OFR,
283+ nSamples,
284+ nFrames,
285+ samplesPerFrame,
286+ dt,
287+ stepsPerSample,
288+ omegaMax,
289+ nBead,
290+ nHinge,
291+ hasMallet: !!mallet,
292+ freeRigid: !!sys.freeRigid,
293+ velForces,
294+ };
295+}
296+
297+// Workgroup sums: each invocation puts its values, then wsum(K) leaves the
298+// sums in wres[0..K] (all invocations must call it). wmax: the largest value
299+// put at 0. With subgroups: a subgroup sum per value, then one invocation per
300+// value adds the subgroups' sums (subgroups taken to be consecutive
301+// invocations); else a tree in workgroup memory.
302+function reducers(sg: boolean): string {
303+ const common = `
304+var<workgroup> wred: array<f32, 17u * 64u>;
305+var<workgroup> wres: array<f32, 17>;
306+var<private> sgi: u32;
307+var<private> sgs: u32;
308+fn put(l: u32, k: u32, v: f32) { wred[k * 64u + l] = v; }
309+fn put3(l: u32, k: u32, v: vec3f) { put(l, k, v.x); put(l, k + 1u, v.y); put(l, k + 2u, v.z); }
310+`;
311+ if (sg)
312+ return common + `
313+var<workgroup> wsg: array<f32, 17u * 16u>;
314+fn wsum(l: u32, K: u32) {
315+ workgroupBarrier();
316+ let g = l / sgs;
317+ for (var k = 0u; k < K; k++) {
318+ let s = subgroupAdd(wred[k * 64u + l]);
319+ if (sgi == 0u) { wsg[k * 16u + g] = s; }
320+ }
321+ workgroupBarrier();
322+ if (l < K) {
323+ var t = 0.0;
324+ for (var j = 0u; j < 64u / sgs; j++) { t += wsg[l * 16u + j]; }
325+ wres[l] = t;
326+ }
327+ workgroupBarrier();
328+}
329+fn wmax(l: u32) -> f32 {
330+ workgroupBarrier();
331+ let s = subgroupMax(wred[l]);
332+ if (sgi == 0u) { wsg[l / sgs] = s; }
333+ workgroupBarrier();
334+ var m = 0.0;
335+ for (var j = 0u; j < 64u / sgs; j++) { m = max(m, wsg[j]); }
336+ workgroupBarrier();
337+ return m;
338+}
339+`;
340+ return common + `
341+fn wsum(l: u32, K: u32) {
342+ workgroupBarrier();
343+ // each level adds pairs (k, j) and (k, j + s), spread over all invocations
344+ for (var s = 32u; s > 0u; s = s / 2u) {
345+ for (var x = l; x < K * s; x += 64u) {
346+ let k = x / s;
347+ let j = x - k * s;
348+ wred[k * 64u + j] += wred[k * 64u + j + s];
349+ }
350+ workgroupBarrier();
351+ }
352+ if (l < K) { wres[l] = wred[l * 64u]; }
353+ workgroupBarrier();
354+}
355+fn wmax(l: u32) -> f32 {
356+ workgroupBarrier();
357+ for (var s = 32u; s > 0u; s = s / 2u) {
358+ if (l < s) { wred[l] = max(wred[l], wred[l + s]); }
359+ workgroupBarrier();
360+ }
361+ let m = wred[0];
362+ workgroupBarrier();
363+ return m;
364+}
365+`;
366+}
367+
368+// Shared scalars (at SC), written by invocation 0 of workgroup 0:
369+// 0-2 mallet position q, 3-5 velocity, 6-8 acceleration, 9-11 force on the
370+// mallet, 12 contact force, 13 contact energy, 14 friction power,
371+// 17-18 dissipated since the last frame (Kahan sum), 19-20 work in since the
372+// last frame, 44 time steps taken and 45 next sample to record (u32 bits).
373+// Each workgroup keeps its own copy of the rigid-motion correction (wcor).
374+const KERNELS = /* wgsl */ `
375+@group(0) @binding(0) var<storage, read> ro: array<u32>;
376+@group(0) @binding(1) var<storage, read_write> rw: array<f32>;
377+@group(0) @binding(2) var<storage, read_write> ob: array<f32>;
378+
379+fn P(k: u32) -> f32 { return bitcast<f32>(ro[k]); }
380+fn pv(o: u32) -> vec3f { return vec3f(P(o), P(o + 1u), P(o + 2u)); }
381+fn ld(o: u32) -> vec3f { return vec3f(rw[o], rw[o + 1u], rw[o + 2u]); }
382+fn sv(o: u32, x: vec3f) { rw[o] = x.x; rw[o + 1u] = x.y; rw[o + 2u] = x.z; }
383+fn s3(k: u32) -> vec3f { return ld(SC + k); }
384+fn ss3(k: u32, x: vec3f) { sv(SC + k, x); }
385+fn count(k: u32) -> u32 { return bitcast<u32>(rw[SC + k]); }
386+fn gU(i: u32) -> vec3f { return ld(S_U + 3u * i); }
387+fn gV(i: u32) -> vec3f { return ld(S_V + 3u * i); }
388+fn gA(i: u32) -> vec3f { return ld(S_A + 3u * i); }
389+fn sU(i: u32, x: vec3f) { sv(S_U + 3u * i, x); }
390+fn sV(i: u32, x: vec3f) { sv(S_V + 3u * i, x); }
391+fn sA(i: u32, x: vec3f) { sv(S_A + 3u * i, x); }
392+// position of bead i relative to the rest centre of mass
393+fn pos(i: u32) -> vec3f { return pv(P_R0 + 3u * i) + gU(i); }
394+
395+// displacement or velocity of a spring's far node: a bead, a fixed anchor (0)
396+// or a sliding anchor that follows bead b across its rod
397+fn otherU(o: u32) -> vec3f {
398+ if (o == 0xffffffffu) { return vec3f(0.0); }
399+ if (o >= 0x80000000u) { let x = gU(o - 0x80000000u); return vec3f(0.0, x.y, x.z); }
400+ return gU(o);
401+}
402+fn otherV(o: u32) -> vec3f {
403+ if (o == 0xffffffffu) { return vec3f(0.0); }
404+ if (o >= 0x80000000u) { let x = gV(o - 0x80000000u); return vec3f(0.0, x.y, x.z); }
405+ return gV(o);
406+}
407+
408+fn kahan(k: u32, x: f32) {
409+ let o = SC + k;
410+ let y = x - rw[o + 1u];
411+ let t = rw[o] + y;
412+ rw[o + 1u] = (t - rw[o]) - y;
413+ rw[o] = t;
414+}
415+
416+// sums of per-workgroup partials rw[q + w * stride + k] over w < nw into wres[0..K]
417+// (the loads are spread over the invocations: global memory latency, not
418+// arithmetic, is what costs here)
419+fn psum(l: u32, q: u32, K: u32, stride: u32, nw: u32) {
420+ for (var k = 0u; k < K; k++) { put(l, k, 0.0); }
421+ for (var x = l; x < nw * K; x += 64u) {
422+ let w = x / K;
423+ let k = x - w * K;
424+ wred[k * 64u + l] += rw[q + w * stride + k];
425+ }
426+ wsum(l, K);
427+}
428+
429+// sums for the rigid part of a per-bead field x: mass-weighted geometry
430+// (m u, m s s^T as diagonal and off-diagonal) and momentum (m x, m s x x),
431+// s = position relative to the rest centre of mass
432+struct Rig { mu: vec3f, sd: vec3f, so: vec3f, p: vec3f, l: vec3f };
433+fn rigAdd(r: Rig, m: f32, u: vec3f, s: vec3f, x: vec3f) -> Rig {
434+ return Rig(r.mu + m * u, r.sd + m * s * s, r.so + m * vec3f(s.x * s.y, s.x * s.z, s.y * s.z), r.p + m * x, r.l + m * cross(s, x));
435+}
436+fn rigZero() -> Rig { return Rig(vec3f(0.0), vec3f(0.0), vec3f(0.0), vec3f(0.0), vec3f(0.0)); }
437+fn rigPut(l: u32, r: Rig) { put3(l, 0u, r.mu); put3(l, 3u, r.sd); put3(l, 6u, r.so); put3(l, 9u, r.p); put3(l, 12u, r.l); }
438+
439+// This workgroup's rigid-motion correction: 0-2 centre of mass, 3-8 inverse
440+// inertia, 9-11 rigid translation rate, 12-14 rotation rate, 15 kinetic
441+// energy of the rigid part, 16 damping power, 17 external power (from the
442+// bead-force sums).
443+var<workgroup> wcor: array<f32, 18>;
444+fn w3(k: u32) -> vec3f { return vec3f(wcor[k], wcor[k + 1u], wcor[k + 2u]); }
445+fn sw3(k: u32, x: vec3f) { wcor[k] = x.x; wcor[k + 1u] = x.y; wcor[k + 2u] = x.z; }
446+// invocation 0: centre of mass and inverse inertia from geometry sums in wres[0..9]
447+fn rigidGeometry() {
448+ let c = vec3f(wres[0], wres[1], wres[2]) / MTOT;
449+ let sxx = wres[3] - MTOT * c.x * c.x;
450+ let syy = wres[4] - MTOT * c.y * c.y;
451+ let szz = wres[5] - MTOT * c.z * c.z;
452+ let sxy = wres[6] - MTOT * c.x * c.y;
453+ let sxz = wres[7] - MTOT * c.x * c.z;
454+ let syz = wres[8] - MTOT * c.y * c.z;
455+ let a = syy + szz; let b = -sxy; let cc = -sxz;
456+ let d = sxx + szz; let e = -syz; let ff = sxx + syy;
457+ let A00 = d * ff - e * e; let A01 = cc * e - b * ff; let A02 = b * e - cc * d;
458+ let A11 = a * ff - cc * cc; let A12 = b * cc - a * e; let A22 = a * d - b * b;
459+ let det = a * A00 + b * A01 + cc * A02;
460+ let k = select(0.0, 1.0 / det, abs(det) > 0.0);
461+ wcor[3] = A00 * k; wcor[4] = A01 * k; wcor[5] = A02 * k;
462+ wcor[6] = A11 * k; wcor[7] = A12 * k; wcor[8] = A22 * k;
463+ sw3(0u, c);
464+}
465+// invocation 0: rigid translation and rotation rate of a field from its
466+// momentum sums in wres[o..o+6], and its kinetic energy
467+fn rigidPart(o: u32) {
468+ let P = vec3f(wres[o], wres[o + 1u], wres[o + 2u]);
469+ let L = vec3f(wres[o + 3u], wres[o + 4u], wres[o + 5u]) - cross(w3(0u), P);
470+ let w = vec3f(wcor[3] * L.x + wcor[4] * L.y + wcor[5] * L.z, wcor[4] * L.x + wcor[6] * L.y + wcor[7] * L.z, wcor[5] * L.x + wcor[7] * L.y + wcor[8] * L.z);
471+ let vc = P / MTOT;
472+ sw3(9u, vc);
473+ sw3(12u, w);
474+ wcor[15] = 0.5 * dot(P, vc) + 0.5 * dot(w, L);
475+}
476+fn rigidAt(s: vec3f) -> vec3f { return w3(9u) + cross(w3(12u), s - w3(0u)); }
477+
478+// all invocations: the geometry and the rigid part of the accelerations from
479+// the bead-force partials, and the damping and external power
480+fn accelRigid(l: u32) {
481+ psum(l, Q_BF, 17u, 17u, NBW);
482+ if (l == 0u) {
483+ rigidGeometry();
484+ rigidPart(9u);
485+ wcor[16] = wres[15];
486+ wcor[17] = wres[16];
487+ }
488+ workgroupBarrier();
489+}
490+// all invocations: the geometry from the bead-force partials and the rigid
491+// part of the velocities from the kick partials
492+fn velRigid(l: u32) {
493+ psum(l, Q_BF, 9u, 17u, NBW);
494+ if (l == 0u) { rigidGeometry(); }
495+ workgroupBarrier();
496+ psum(l, Q_KC, 6u, 6u, NBW);
497+ if (l == 0u) { rigidPart(0u); }
498+ workgroupBarrier();
499+}
500+
501+// dihedral angle change and its gradient (Bridson et al. 2003; see hinge.ts)
502+fn hinge(e: u32) {
503+ let w0 = ro[I_HN + 4u * e];
504+ let w1 = ro[I_HN + 4u * e + 1u];
505+ let e0 = ro[I_HN + 4u * e + 2u];
506+ let e1 = ro[I_HN + 4u * e + 3u];
507+ let hb = P_HF + 11u * e;
508+ let sg = P(hb);
509+ let ue0 = gU(e0);
510+ let ue1 = gU(e1);
511+ let E = pv(hb + 2u) + (ue1 - ue0);
512+ let P1 = pv(hb + 5u) + (gU(w0) - ue0);
513+ let P2 = pv(hb + 8u) + (gU(w1) - ue1);
514+ let Q1 = P1 - E;
515+ let Q2 = P2 + E;
516+ let N1 = cross(P1, Q1);
517+ let N2 = cross(P2, Q2);
518+ let n1sq = dot(N1, N1);
519+ let n2sq = dot(N2, N2);
520+ let eLen = length(E);
521+ let nn = sqrt(n1sq) * sqrt(n2sq);
522+ let theta = atan2(dot(cross(N1, N2), E) / (nn * eLen), dot(N1, N2) / nn);
523+ let gw0 = (-eLen / n1sq) * N1;
524+ let gw1 = (-eLen / n2sq) * N2;
525+ let t1 = dot(Q1, E) / eLen / n1sq;
526+ let t2 = dot(P2, E) / eLen / n2sq;
527+ let r1 = dot(P1, E) / eLen / n1sq;
528+ let r2 = dot(Q2, E) / eLen / n2sq;
529+ let ge0 = -(t1 * N1 + t2 * N2);
530+ let ge1 = r1 * N1 + r2 * N2;
531+ let g = S_G + 12u * e;
532+ sv(g, sg * gw0);
533+ sv(g + 3u, sg * gw1);
534+ sv(g + 6u, sg * ge0);
535+ sv(g + 9u, sg * ge1);
536+ rw[S_PHI + e] = sg * (theta - P(hb + 1u));
537+ if (S1 != 0.0) {
538+ rw[S_RATE + e] = sg * (dot(gw0, gV(w0)) + dot(gw1, gV(w1)) + dot(ge0, gV(e0)) + dot(ge1, gV(e1)));
539+ }
540+}
541+
542+// mallet contact force from the patch-weighted displacement U and velocity Vp
543+fn contact(U: vec3f, Vp: vec3f) {
544+ let q = s3(0u);
545+ let qv = s3(3u);
546+ var fm = vec3f(0.0);
547+ var fc = 0.0;
548+ var ec = 0.0;
549+ var pf = 0.0;
550+ let d = dot(q - U, NRM);
551+ if (d > 0.0) {
552+ let o = q - U - d * NRM;
553+ if (REACH2 < 0.0 || dot(o, o) <= REACH2) {
554+ let fnorm = MK * pow(d, MP);
555+ fc = fnorm;
556+ ec = fnorm * d / (MP + 1.0);
557+ fm = -fnorm * NRM;
558+ if (MU > 0.0) {
559+ var r = qv - Vp;
560+ r -= dot(r, NRM) * NRM;
561+ let slip = length(r);
562+ let eps = max(1e-3, 4.0 * MU * fnorm * DT / MEFF);
563+ let g = MU * fnorm / sqrt(slip * slip + eps * eps);
564+ fm -= g * r;
565+ pf = g * slip * slip;
566+ }
567+ }
568+ }
569+ ss3(6u, fm / MM);
570+ ss3(9u, fm);
571+ rw[SC + 12u] = fc;
572+ rw[SC + 13u] = ec;
573+ rw[SC + 14u] = pf;
574+}
575+
576+// invocation 0 of workgroup 0 (after the mallet's position and velocity
577+// are final for the step): the mallet's half kick and drift, count the step
578+fn malletDrift() {
579+ if (HAS_MALLET) {
580+ let qv = s3(3u) + 0.5 * DT * s3(6u);
581+ ss3(3u, qv);
582+ ss3(0u, s3(0u) + DT * qv);
583+ }
584+ rw[SC + 44u] = bitcast<f32>(count(44u) + 1u);
585+}
586+
587+// workgroup 0: sum the record partials and write the sample (and frame),
588+// then move to the next sample
589+fn recordSums(l: u32) {
590+ var mx = 0.0;
591+ for (var w = l; w < NBW + NHW; w += 64u) { mx = max(mx, rw[Q_REC + 3u * w + 2u]); }
592+ put(l, 0u, mx);
593+ let largest = wmax(l);
594+ psum(l, Q_REC, 3u, 3u, NBW + NHW);
595+ if (l == 0u) {
596+ let smp = count(45u);
597+ ob[O_AUD + smp] = wres[0];
598+ ob[O_CF + smp] = rw[SC + 12u];
599+ if (smp % SPF == 0u && smp / SPF < NF) {
600+ let fi = smp / SPF;
601+ let b = O_FSC + 6u * fi;
602+ let qv = s3(3u);
603+ ob[b] = wres[1];
604+ ob[b + 1u] = select(0.0, 0.5 * MM * dot(qv, qv) + rw[SC + 13u], HAS_MALLET);
605+ ob[b + 2u] = rw[SC + 17u];
606+ ob[b + 3u] = rw[SC + 19u];
607+ ob[b + 4u] = largest;
608+ rw[SC + 17u] = 0.0; rw[SC + 18u] = 0.0; rw[SC + 19u] = 0.0; rw[SC + 20u] = 0.0;
609+ let q = s3(0u);
610+ ob[O_MP + 3u * fi] = q.x;
611+ ob[O_MP + 3u * fi + 1u] = q.y;
612+ ob[O_MP + 3u * fi + 2u] = q.z;
613+ }
614+ rw[SC + 45u] = bitcast<f32>(smp + 1u);
615+ }
616+}
617+
618+// first kick and drift of a sample's steps; workgroup 0 also writes the
619+// sample recorded just before
620+@compute @workgroup_size(64)
621+fn kickDriftRec(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
622+ /*SGI*/
623+ if (wg.x == 0u) {
624+ recordSums(l);
625+ if (l == 0u) { malletDrift(); }
626+ }
627+ let i = wg.x * 64u + l;
628+ if (i >= N || P(P_IM + i) == 0.0) { return; }
629+ let v = gV(i) + 0.5 * DT * gA(i);
630+ sV(i, v);
631+ sU(i, gU(i) + DT * v);
632+}
633+
634+// later kicks and drifts: first remove the rigid part of the velocities left
635+// by the previous step
636+@compute @workgroup_size(64)
637+fn kickDrift(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
638+ /*SGI*/
639+ if (FREE_RIGID) { velRigid(l); }
640+ if (wg.x == 0u && l == 0u) {
641+ if (FREE_RIGID) { kahan(17u, wcor[15]); }
642+ malletDrift();
643+ }
644+ let i = wg.x * 64u + l;
645+ if (i >= N || P(P_IM + i) == 0.0) { return; }
646+ var v = gV(i);
647+ if (FREE_RIGID) { v -= rigidAt(pos(i)); }
648+ v += 0.5 * DT * gA(i);
649+ sV(i, v);
650+ sU(i, gU(i) + DT * v);
651+}
652+
653+// after the last step of a sample: remove the rigid part of the velocities
654+@compute @workgroup_size(64)
655+fn velCorrect(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
656+ /*SGI*/
657+ velRigid(l);
658+ if (wg.x == 0u && l == 0u) { kahan(17u, wcor[15]); }
659+ let i = wg.x * 64u + l;
660+ if (i >= N) { return; }
661+ sV(i, gV(i) - rigidAt(pos(i)));
662+}
663+
664+// hinge angles, gradients and rates; the last workgroup does the mallet contact
665+@compute @workgroup_size(64)
666+fn hinges(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
667+ /*SGI*/
668+ if (wg.x < NHW) {
669+ let e = wg.x * 64u + l;
670+ if (e < NH) { hinge(e); }
671+ return;
672+ }
673+ var su = vec3f(0.0);
674+ var sv = vec3f(0.0);
675+ for (var j = l; j < NP; j += 64u) {
676+ let i = ro[I_PATCH + j];
677+ let w = P(P_PW + i);
678+ su += w * gU(i);
679+ sv += w * gV(i);
680+ }
681+ put3(l, 0u, su);
682+ put3(l, 3u, sv);
683+ wsum(l, 6u);
684+ if (l == 0u) { contact(vec3f(wres[0], wres[1], wres[2]), vec3f(wres[3], wres[4], wres[5])); }
685+}
686+
687+// hinge moments, elastic and with damping; per-workgroup hinge damping power
688+@compute @workgroup_size(64)
689+fn moments(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
690+ /*SGI*/
691+ let e = wg.x * 64u + l;
692+ var pd = 0.0;
693+ if (e < NH) {
694+ var m = 0.0;
695+ var r = 0.0;
696+ for (var q = 0u; q < MAXHC; q++) {
697+ let j = q * NH + e;
698+ let c = ro[I_HCC + j];
699+ let k = P(P_HC + j);
700+ m += k * rw[S_PHI + c];
701+ if (S1 != 0.0) { r += k * rw[S_RATE + c]; }
702+ }
703+ rw[S_MEL + e] = m;
704+ rw[S_MOM + e] = m + 2.0 * S1 * r;
705+ if (S1 != 0.0) { pd = 2.0 * S1 * rw[S_RATE + e] * r; }
706+ }
707+ put(l, 0u, pd);
708+ wsum(l, 1u);
709+ if (l == 0u) { rw[Q_HD + wg.x] = wres[0]; }
710+}
711+
712+// accelerations (before rigid removal); per-workgroup damping and external
713+// power and sums for the rigid part of the accelerations
714+@compute @workgroup_size(64)
715+fn beadForces(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
716+ /*SGI*/
717+ let i = wg.x * 64u + l;
718+ let t = f32(count(44u)) * DT;
719+ var pd = 0.0;
720+ var pe = 0.0;
721+ var rig = rigZero();
722+ if (i < N) {
723+ let im = P(P_IM + i);
724+ let ui = gU(i);
725+ var ai = vec3f(0.0);
726+ if (im != 0.0) {
727+ let vi = gV(i);
728+ var fi = vec3f(0.0);
729+ for (var q = 0u; q < MAXS; q++) {
730+ let j = q * N + i;
731+ let o = ro[I_SO + j];
732+ let b = P_SF + 8u * j;
733+ let rd = pv(b + 2u);
734+ let du = otherU(o) - ui;
735+ let d = rd + du;
736+ let len = length(d);
737+ let dl = (2.0 * dot(rd, du) + dot(du, du)) / (len + P(b + 5u));
738+ fi += (P(b) * (dl + P(b + 6u)) / len) * d;
739+ if (ANY_DASH) {
740+ let c = P(b + 1u);
741+ let vr = otherV(o) - vi;
742+ fi += c * vr;
743+ pd += P(b + 7u) * c * dot(vr, vr);
744+ }
745+ }
746+ if (NH > 0u) {
747+ for (var q = 0u; q < MAXBH; q++) {
748+ let ent = ro[I_BHE + q * N + i];
749+ let e = ent >> 2u;
750+ fi -= rw[S_MOM + e] * ld(S_G + 12u * e + 3u * (ent & 3u));
751+ }
752+ }
753+ if (HAS_MALLET) { fi -= P(P_PW + i) * s3(9u); }
754+ if (S0 != 0.0) {
755+ let m2 = 2.0 * S0 * P(P_MASS + i);
756+ fi -= m2 * vi;
757+ pd += m2 * dot(vi, vi);
758+ }
759+ if (HAS_PULSE && t >= 0.0 && t <= TAU) {
760+ let fe = F0 * 0.5 * (1.0 - cos(6.283185307 * t / TAU)) * pv(P_FS + 3u * i);
761+ fi += fe;
762+ pe += dot(fe, vi);
763+ }
764+ ai = fi * im;
765+ }
766+ sA(i, ai);
767+ if (FREE_RIGID) { rig = rigAdd(rig, P(P_MASS + i), ui, pv(P_R0 + 3u * i) + ui, ai); }
768+ }
769+ rigPut(l, rig);
770+ put(l, 15u, pd);
771+ put(l, 16u, pe);
772+ wsum(l, 17u);
773+ if (l < 17u) { rw[Q_BF + wg.x * 17u + l] = wres[l]; }
774+}
775+
776+// remove the rigid part of the accelerations, second half kick; per-workgroup
777+// sums for the rigid part of the velocities. Workgroup 0: the mallet's second
778+// half kick, dissipated energy and work in.
779+@compute @workgroup_size(64)
780+fn kick(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
781+ /*SGI*/
782+ if (FREE_RIGID || wg.x == 0u) { accelRigid(l); }
783+ if (wg.x == 0u && l == 0u) {
784+ var hd = 0.0;
785+ for (var w = 0u; w < NHW; w++) { hd += rw[Q_HD + w]; }
786+ if (HAS_MALLET) { ss3(3u, s3(3u) + 0.5 * DT * s3(6u)); }
787+ kahan(17u, DT * (wcor[16] + hd + rw[SC + 14u]));
788+ kahan(19u, DT * wcor[17]);
789+ }
790+ let i = wg.x * 64u + l;
791+ var p = vec3f(0.0);
792+ var lm = vec3f(0.0);
793+ if (i < N) {
794+ var a = gA(i);
795+ var v = gV(i);
796+ if (FREE_RIGID) {
797+ let s = pos(i);
798+ a -= rigidAt(s);
799+ sA(i, a);
800+ if (P(P_IM + i) != 0.0) { v += 0.5 * DT * a; }
801+ let m = P(P_MASS + i);
802+ p = m * v;
803+ lm = m * cross(s, v);
804+ } else if (P(P_IM + i) != 0.0) {
805+ v += 0.5 * DT * a;
806+ }
807+ sV(i, v);
808+ }
809+ if (FREE_RIGID) {
810+ put3(l, 0u, p);
811+ put3(l, 3u, lm);
812+ wsum(l, 6u);
813+ if (l < 6u) { rw[Q_KC + wg.x * 6u + l] = wres[l]; }
814+ }
815+}
816+
817+// per-workgroup audio and (at frames) energy and largest strain for the next
818+// sample: bead workgroups first, then hinge workgroups (hinge energy). With
819+// correct = true the accelerations still hold their rigid part (forces were
820+// just recomputed) and it is removed first.
821+fn record(wg: u32, l: u32, correct: bool) {
822+ let smp = count(45u);
823+ let isFrame = smp % SPF == 0u && smp / SPF < NF;
824+ var ap = 0.0;
825+ var ep = 0.0;
826+ var mx = 0.0;
827+ if (wg < NBW) {
828+ if (correct) { accelRigid(l); }
829+ let i = wg * 64u + l;
830+ if (i < N) {
831+ var a = gA(i);
832+ if (correct) {
833+ a -= rigidAt(pos(i));
834+ sA(i, a);
835+ }
836+ let hb = S_HIST + i * H;
837+ rw[hb + smp % H] = dot(pv(P_W + 3u * i), a);
838+ let i0 = ro[I_I0 + i];
839+ let fr = P(P_FR + i);
840+ if (smp >= i0) { ap += (1.0 - fr) * rw[hb + (smp - i0) % H]; }
841+ if (smp >= i0 + 1u) { ap += fr * rw[hb + (smp - i0 - 1u) % H]; }
842+ if (isFrame) {
843+ let ui = gU(i);
844+ let fo = OFR + (smp / SPF) * 3u * N + 3u * i;
845+ ob[fo] = ui.x;
846+ ob[fo + 1u] = ui.y;
847+ ob[fo + 2u] = ui.z;
848+ let vi = gV(i);
849+ ep += 0.5 * P(P_MASS + i) * dot(vi, vi);
850+ if (P(P_IM + i) != 0.0) {
851+ for (var q = 0u; q < MAXS; q++) {
852+ let j = q * N + i;
853+ let o = ro[I_SO + j];
854+ let b = P_SF + 8u * j;
855+ let rd = pv(b + 2u);
856+ let lr = P(b + 5u);
857+ let e0 = P(b + 6u);
858+ let du = otherU(o) - ui;
859+ let len = length(rd + du);
860+ let dl = (2.0 * dot(rd, du) + dot(du, du)) / (len + lr);
861+ // energy above the pre-tensioned equilibrium: k/2 ((dl + e0)^2 - e0^2)
862+ ep += P(b + 7u) * 0.5 * P(b) * dl * (dl + 2.0 * e0);
863+ if (o < N) { mx = max(mx, abs(dl + e0) / (lr - e0)); }
864+ }
865+ }
866+ }
867+ }
868+ } else if (isFrame) {
869+ let e = (wg - NBW) * 64u + l;
870+ if (e < NH) { ep = 0.5 * rw[S_PHI + e] * rw[S_MEL + e]; }
871+ }
872+ put(l, 0u, mx);
873+ let largest = wmax(l);
874+ put(l, 0u, ap);
875+ put(l, 1u, ep);
876+ wsum(l, 2u);
877+ if (l == 0u) {
878+ rw[Q_REC + wg * 3u] = wres[0];
879+ rw[Q_REC + wg * 3u + 1u] = wres[1];
880+ rw[Q_REC + wg * 3u + 2u] = largest;
881+ }
882+}
883+@compute @workgroup_size(64)
884+fn recordSample(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
885+ /*SGI*/ record(wg.x, l, false); }
886+@compute @workgroup_size(64)
887+fn recordSampleCorrect(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
888+ /*SGI*/ record(wg.x, l, FREE_RIGID); }
889+// writes the last sample
890+@compute @workgroup_size(64)
891+fn recordLast(@builtin(local_invocation_index) l: u32 /*SG*/) {
892+ /*SGI*/ recordSums(l); }
893+
894+// at the start: per-workgroup sums for the rigid part of the initial velocities
895+@compute @workgroup_size(64)
896+fn initRigid(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
897+ /*SGI*/
898+ let i = wg.x * 64u + l;
899+ var rig = rigZero();
900+ if (i < N) {
901+ let u = gU(i);
902+ rig = rigAdd(rig, P(P_MASS + i), u, pv(P_R0 + 3u * i) + u, gV(i));
903+ }
904+ rigPut(l, rig);
905+ wsum(l, 15u);
906+ if (l < 15u) { rw[Q_IRV + wg.x * 15u + l] = wres[l]; }
907+}
908+// and its removal
909+@compute @workgroup_size(64)
910+fn initVel(@builtin(workgroup_id) wg: vec3u, @builtin(local_invocation_index) l: u32 /*SG*/) {
911+ /*SGI*/
912+ psum(l, Q_IRV, 15u, 15u, NBW);
913+ if (l == 0u) {
914+ rigidGeometry();
915+ rigidPart(9u);
916+ }
917+ workgroupBarrier();
918+ let i = wg.x * 64u + l;
919+ if (i >= N) { return; }
920+ sV(i, gV(i) - rigidAt(pos(i)));
921+}
922+`;
923+
924+const ENTRY = ['kickDriftRec', 'kickDrift', 'velCorrect', 'hinges', 'moments', 'beadForces', 'kick', 'recordSample', 'recordSampleCorrect', 'recordLast', 'initRigid', 'initVel'] as const;
925+type Entry = (typeof ENTRY)[number];
926+
927+export async function simulateGPU(
928+ sys: BeadSystem,
929+ exc: Excitation,
930+ opt: SimOptions,
931+ progress?: Progress,
932+): Promise<SimResult & { device: string; loopSeconds: number }> {
933+ const t0 = performance.now();
934+ const gpu = navigator.gpu;
935+ if (!gpu) throw new Error('WebGPU is not available in this browser');
936+ const adapter = await gpu.requestAdapter({ powerPreference: 'high-performance' });
937+ if (!adapter) throw new Error('no WebGPU adapter');
938+ const lim = adapter.limits;
939+ const sg = adapter.features.has('subgroups');
940+ const pk = pack(sys, exc, opt, sg);
941+ const biggest = Math.max(pk.ro.byteLength, pk.rw.byteLength, pk.outLen * 4);
942+ if (biggest > lim.maxStorageBufferBindingSize || biggest > lim.maxBufferSize)
943+ throw new Error(`a GPU buffer would need ${(biggest / 1e6).toFixed(0)} MB (limit ${(Math.min(lim.maxStorageBufferBindingSize, lim.maxBufferSize) / 1e6).toFixed(0)} MB); lower the replay frame rate or the duration`);
944+ const device = await adapter.requestDevice({
945+ requiredFeatures: sg ? ['subgroups' as GPUFeatureName] : [],
946+ requiredLimits: {
947+ maxStorageBufferBindingSize: Math.max(biggest, 1 << 20),
948+ maxBufferSize: Math.max(biggest, 1 << 20),
949+ },
950+ });
951+ let lost = '';
952+ device.lost.then((info) => (lost = info.message || 'GPU device lost'));
953+ try {
954+ device.pushErrorScope('validation');
955+ const module = device.createShaderModule({ code: pk.wgsl });
956+ const info = await module.getCompilationInfo();
957+ const errs = info.messages.filter((m) => m.type === 'error');
958+ if (errs.length) throw new Error(`shader: ${errs.map((m) => `${m.lineNum}:${m.linePos} ${m.message}`).join('; ')}`);
959+ const bgl = device.createBindGroupLayout({
960+ entries: (['read-only-storage', 'storage', 'storage'] as const).map((type, binding) => ({ binding, visibility: GPUShaderStage.COMPUTE, buffer: { type } })),
961+ });
962+ const layout = device.createPipelineLayout({ bindGroupLayouts: [bgl] });
963+ const pipes = {} as Record<Entry, GPUComputePipeline>;
964+ await Promise.all(ENTRY.map(async (e) => (pipes[e] = await device.createComputePipelineAsync({ layout, compute: { module, entryPoint: e } }))));
965+ const mk = (data: Float32Array | Uint32Array | null, size: number) => {
966+ const buf = device.createBuffer({ size: Math.max(16, size), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC });
967+ if (data) device.queue.writeBuffer(buf, 0, data as Float32Array<ArrayBuffer>);
968+ return buf;
969+ };
970+ const outB = mk(null, pk.outLen * 4);
971+ const bind = device.createBindGroup({
972+ layout: bgl,
973+ entries: [mk(pk.ro, pk.ro.byteLength), mk(pk.rw, pk.rw.byteLength), outB].map((buffer, binding) => ({ binding, resource: { buffer } })),
974+ });
975+
976+ const { nBead, nHinge, hasMallet, freeRigid, velForces, stepsPerSample } = pk;
977+ const nh = sys.bending ? sys.bending.w0.length : 0;
978+ type D = [Entry, number];
979+ const forces: D[] = [];
980+ if (nHinge + (hasMallet ? 1 : 0) > 0) forces.push(['hinges', nHinge + (hasMallet ? 1 : 0)]);
981+ if (nh > 0) forces.push(['moments', nHinge]);
982+ forces.push(['beadForces', nBead]);
983+ // one sample: record it, then the steps to the next one; the accelerations
984+ // hold their rigid part after a force evaluation until a kick removes it
985+ const sample = (first: boolean, last: boolean): D[] => {
986+ const raw = first || velForces;
987+ const list: D[] = [[raw ? 'recordSampleCorrect' : 'recordSample', nBead + nHinge]];
988+ if (last) return [...list, ['recordLast', 1]];
989+ for (let j = 0; j < stepsPerSample; j++) list.push([j === 0 ? 'kickDriftRec' : 'kickDrift', nBead], ...forces, ['kick', nBead]);
990+ if (freeRigid) list.push(['velCorrect', nBead]);
991+ if (velForces) list.push(...forces);
992+ return list;
993+ };
994+ const start: D[] = [...(freeRigid ? ([['initRigid', nBead], ['initVel', nBead]] as D[]) : []), ...forces];
995+ const encode = (pass: GPUComputePassEncoder, list: D[]) => {
996+ for (const [e, w] of list) {
997+ pass.setPipeline(pipes[e]);
998+ pass.dispatchWorkgroups(w);
999+ }
1000+ };
1001+ const submit = (fill: (pass: GPUComputePassEncoder) => void) => {
1002+ const enc = device.createCommandEncoder();
1003+ const pass = enc.beginComputePass();
1004+ pass.setBindGroup(0, bind);
1005+ fill(pass);
1006+ pass.end();
1007+ device.queue.submit([enc.finish()]);
1008+ return device.queue.onSubmittedWorkDone();
1009+ };
1010+ await submit((pass) => encode(pass, start));
1011+ const err = await device.popErrorScope();
1012+ if (err) throw new Error(err.message);
1013+
1014+ // chunks of samples, about 50 ms of GPU work each, two in flight so the
1015+ // GPU does not wait while the next one is encoded
1016+ const tLoop = performance.now();
1017+ const mid = sample(false, false);
1018+ let chunk = 16;
1019+ let pending: Promise<void> | null = null;
1020+ let pendingEnd = 0;
1021+ let tPrev = performance.now();
1022+ for (let s = 0; s < pk.nSamples; ) {
1023+ const e = Math.min(pk.nSamples, s + chunk);
1024+ const done = submit((pass) => {
1025+ for (let k = s; k < e; k++) encode(pass, k === 0 || k === pk.nSamples - 1 ? sample(k === 0, k === pk.nSamples - 1) : mid);
1026+ });
1027+ if (pending) {
1028+ await pending;
1029+ if (lost) throw new Error(lost);
1030+ const now = performance.now();
1031+ chunk = Math.max(1, Math.min(Math.round(chunk * Math.min(2, 50 / Math.max(now - tPrev, 1))), 4096));
1032+ tPrev = now;
1033+ progress?.(pendingEnd / pk.nSamples);
1034+ }
1035+ pending = done;
1036+ pendingEnd = e;
1037+ s = e;
1038+ }
1039+ await pending;
1040+ if (lost) throw new Error(lost);
1041+ progress?.(1);
1042+ const loopSeconds = (performance.now() - tLoop) / 1000;
1043+
1044+ const stage = device.createBuffer({ size: Math.max(16, pk.outLen * 4), usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
1045+ const enc = device.createCommandEncoder();
1046+ enc.copyBufferToBuffer(outB, 0, stage, 0, Math.max(16, pk.outLen * 4));
1047+ device.queue.submit([enc.finish()]);
1048+ await stage.mapAsync(GPUMapMode.READ);
1049+ const out = new Float32Array(stage.getMappedRange().slice(0, pk.outLen * 4));
1050+ stage.destroy();
1051+ if (lost) throw new Error(lost);
1052+
1053+ const { nSamples, nFrames, samplesPerFrame, dt } = pk;
1054+ const audio = out.slice(0, nSamples);
1055+ for (let k = 0; k < audio.length; k++) if (!Number.isFinite(audio[k])) throw new Error('the GPU run blew up (non-finite output)');
1056+ const contactForce = exc.mallet ? out.slice(nSamples, 2 * nSamples) : undefined;
1057+ const O_FSC = 2 * nSamples;
1058+ const O_MP = O_FSC + 6 * nFrames;
1059+ const frameTimes = new Float64Array(nFrames);
1060+ const energy = new Float64Array(nFrames);
1061+ const balance = new Float64Array(nFrames);
1062+ let maxStretch = 0;
1063+ let diss = 0;
1064+ let work = 0;
1065+ let E0 = 0;
1066+ for (let k = 0; k < nFrames; k++) {
1067+ const b = O_FSC + 6 * k;
1068+ frameTimes[k] = k * samplesPerFrame * stepsPerSample * dt;
1069+ energy[k] = out[b];
1070+ const total = out[b] + out[b + 1];
1071+ diss += out[b + 2];
1072+ work += out[b + 3];
1073+ if (k === 0) E0 = total;
1074+ balance[k] = total + diss - work - E0;
1075+ maxStretch = Math.max(maxStretch, out[b + 4]);
1076+ }
1077+ return {
1078+ audio,
1079+ fs: opt.fs,
1080+ frames: out.slice(pk.framesAt),
1081+ frameTimes,
1082+ energy,
1083+ balance,
1084+ dt,
1085+ stepsPerSample,
1086+ omegaMax: pk.omegaMax,
1087+ wallSeconds: (performance.now() - t0) / 1000,
1088+ maxStretch,
1089+ contactForce,
1090+ malletPos: exc.mallet ? out.slice(O_MP, O_MP + 3 * nFrames) : undefined,
1091+ loopSeconds,
1092+ device: [adapter.info?.vendor, adapter.info?.architecture].filter(Boolean).join(' ') || 'GPU',
1093+ };
1094+ } finally {
1095+ device.destroy();
1096+ }
1097+}
1098+
src/sim/run.tsmodified+2−2View file
@@ -5,7 +5,7 @@ export interface Running {
55 cancel: () => void;
66 }
77
8-export function runInWorker(sys: BeadSystem, exc: Excitation, opt: SimOptions, onProgress: (f: number) => void): Running {
8+export function runInWorker(sys: BeadSystem, exc: Excitation, opt: SimOptions, gpu: boolean, onProgress: (f: number) => void): Running {
99 const worker = new Worker(new URL('./worker.ts', import.meta.url), { type: 'module' });
1010 let reject!: (e: Error) => void;
1111 const promise = new Promise<SimResult>((res, rej) => {
@@ -26,7 +26,7 @@ export function runInWorker(sys: BeadSystem, exc: Excitation, opt: SimOptions, o
2626 rej(new Error(e.message));
2727 };
2828 });
29- worker.postMessage({ sys, exc, opt });
29+ worker.postMessage({ sys, exc, opt, gpu });
3030 return {
3131 promise,
3232 cancel: () => {
src/sim/solver.tsmodified+260−118View file
@@ -17,59 +17,136 @@ function nodePositions(sys: BeadSystem): Float64Array {
1717 return X;
1818 }
1919
20-// Upper bound on the highest angular frequency of the linearized system
21-// (Gershgorin on M^-1 K at rest). A spring under tension T has stiffness k along
22-// its axis and T/l across it, so each spring contributes 2 max(k, |T|/l) / m.
23-// A hinge contributes k g g^T (g = gradient of its angle); a mallet a rank-one
24-// contact stiffness between it and the patch.
25-export function omegaMaxBound(sys: BeadSystem, mallet?: Mallet): number {
20+// Linearized stiffness at rest, y = K x (x, y: 3n, beads only; anchors held).
21+// A spring under tension T has stiffness k along its axis and T/l across it;
22+// bending contributes G^T M G per triangle (G = hinge-angle gradients).
23+function restStiffness(sys: BeadSystem): (x: Float64Array, y: Float64Array) => void {
2624 const { n, si, sj, sk, sL0 } = sys;
27- const rowSum = new Float64Array(n);
2825 const X = nodePositions(sys);
29- for (let s = 0; s < si.length; s++) {
30- const a = si[s];
31- const b = sj[s];
32- const l = Math.hypot(X[3 * b] - X[3 * a], X[3 * b + 1] - X[3 * a + 1], X[3 * b + 2] - X[3 * a + 2]);
33- const T = sk[s] * (l - sL0[s]);
34- const kmax = Math.max(sk[s], Math.abs(T) / l);
35- if (a < n) rowSum[a] += 2 * kmax;
36- if (b < n) rowSum[b] += 2 * kmax;
26+ const ns = si.length;
27+ const sd = new Float64Array(3 * ns); // unit axis
28+ const ka = new Float64Array(ns); // axial stiffness minus transverse
29+ const kt = new Float64Array(ns); // transverse T / l
30+ for (let s = 0; s < ns; s++) {
31+ const a = 3 * si[s], b = 3 * sj[s];
32+ const dx = X[b] - X[a], dy = X[b + 1] - X[a + 1], dz = X[b + 2] - X[a + 2];
33+ const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
34+ sd[3 * s] = dx / l;
35+ sd[3 * s + 1] = dy / l;
36+ sd[3 * s + 2] = dz / l;
37+ kt[s] = (sk[s] * (l - sL0[s])) / l;
38+ ka[s] = sk[s] - kt[s];
3739 }
38- const h = sys.bending;
39- if (h) {
40- const kh = hingeStiffness(h);
40+ const b = sys.bending;
41+ const nh = b ? b.w0.length : 0;
42+ const G = new Float64Array(12 * nh);
43+ const nodes = new Int32Array(4 * nh);
44+ if (b) {
4145 const g = new Float64Array(13);
42- for (let e = 0; e < kh.length; e++) {
43- const nodes = [h.w0[e], h.w1[e], h.e0[e], h.e1[e]];
44- hingeAngleGrad(X, nodes[0], nodes[1], nodes[2], nodes[3], g);
45- let l1 = 0;
46- for (let q = 1; q < 13; q++) l1 += Math.abs(g[q]);
47- nodes.forEach((node, k) => {
48- const inf = Math.max(Math.abs(g[1 + 3 * k]), Math.abs(g[2 + 3 * k]), Math.abs(g[3 + 3 * k]));
49- if (node < n) rowSum[node] += kh[e] * inf * l1;
50- });
46+ for (let e = 0; e < nh; e++) {
47+ hingeAngleGrad(X, b.w0[e], b.w1[e], b.e0[e], b.e1[e], g);
48+ for (let q = 0; q < 12; q++) G[12 * e + q] = b.sign[e] * g[1 + q];
49+ nodes.set([b.w0[e], b.w1[e], b.e0[e], b.e1[e]], 4 * e);
5150 }
5251 }
53- const kc = mallet ? contactStiffness(mallet) : 0;
54- let w2 = mallet ? kc / mallet.mass : 0;
55- for (let i = 0; i < n; i++) {
56- if (sys.pinned[i]) continue;
57- const extra = mallet ? kc * mallet.weights[i] : 0;
58- w2 = Math.max(w2, (rowSum[i] + extra) / sys.mass[i]);
52+ const phi = new Float64Array(nh);
53+ const mom = new Float64Array(nh);
54+ const get = (x: Float64Array, node: number, c: number) => (node < n ? x[3 * node + c] : 0);
55+ return (x, y) => {
56+ y.fill(0);
57+ for (let s = 0; s < ns; s++) {
58+ const a = si[s], c = sj[s];
59+ const rx = get(x, c, 0) - get(x, a, 0), ry = get(x, c, 1) - get(x, a, 1), rz = get(x, c, 2) - get(x, a, 2);
60+ const ex = sd[3 * s], ey = sd[3 * s + 1], ez = sd[3 * s + 2];
61+ const pa = ka[s] * (ex * rx + ey * ry + ez * rz);
62+ const fx = pa * ex + kt[s] * rx, fy = pa * ey + kt[s] * ry, fz = pa * ez + kt[s] * rz;
63+ if (a < n) {
64+ y[3 * a] -= fx;
65+ y[3 * a + 1] -= fy;
66+ y[3 * a + 2] -= fz;
67+ }
68+ if (c < n) {
69+ y[3 * c] += fx;
70+ y[3 * c + 1] += fy;
71+ y[3 * c + 2] += fz;
72+ }
73+ }
74+ if (!b) return;
75+ for (let e = 0; e < nh; e++) {
76+ let v = 0;
77+ for (let k = 0; k < 4; k++) {
78+ const nd = nodes[4 * e + k];
79+ if (nd < n) v += G[12 * e + 3 * k] * x[3 * nd] + G[12 * e + 3 * k + 1] * x[3 * nd + 1] + G[12 * e + 3 * k + 2] * x[3 * nd + 2];
80+ }
81+ phi[e] = v;
82+ }
83+ mom.fill(0);
84+ for (let t = 0; t < b.tri.length / 3; t++)
85+ for (let a = 0; a < 3; a++) {
86+ const ea = b.tri[3 * t + a];
87+ if (ea < 0) continue;
88+ for (let c = 0; c < 3; c++) {
89+ const ec = b.tri[3 * t + c];
90+ if (ec >= 0) mom[ea] += b.M[9 * t + 3 * a + c] * phi[ec];
91+ }
92+ }
93+ for (let e = 0; e < nh; e++)
94+ for (let k = 0; k < 4; k++) {
95+ const nd = nodes[4 * e + k];
96+ if (nd < n) for (let c = 0; c < 3; c++) y[3 * nd + c] += mom[e] * G[12 * e + 3 * k + c];
97+ }
98+ };
99+}
100+
101+// Highest angular frequency of the linearized system at rest: power iteration
102+// on M^-1/2 K M^-1/2 (seeded, so the same system always gives the same value).
103+// The Rayleigh quotient approaches the top eigenvalue from below; the contact
104+// stiffness of a mallet is added on top (Weyl: lambda(A + B) <= lambda(A) + lambda(B)).
105+export function omegaMaxEstimate(sys: BeadSystem, mallet?: Mallet, iterations = 400): number {
106+ const { n, pinned } = sys;
107+ const K = restStiffness(sys);
108+ const rs = Float64Array.from(sys.mass, (m, i) => (pinned[i] ? 0 : 1 / Math.sqrt(m)));
109+ let seed = 12345;
110+ const rand = () => ((seed = (seed * 1103515245 + 12345) % 2147483648) / 2147483648) - 0.5;
111+ const v = new Float64Array(3 * n);
112+ const z = new Float64Array(3 * n);
113+ const y = new Float64Array(3 * n);
114+ for (let k = 0; k < 3 * n; k++) v[k] = rs[Math.floor(k / 3)] > 0 ? rand() : 0;
115+ let lam = 0;
116+ for (let it = 0; it < iterations; it++) {
117+ let nv = 0;
118+ for (let k = 0; k < 3 * n; k++) nv += v[k] * v[k];
119+ nv = Math.sqrt(nv);
120+ if (!(nv > 0)) break;
121+ for (let k = 0; k < 3 * n; k++) z[k] = (v[k] / nv) * rs[Math.floor(k / 3)];
122+ K(z, y);
123+ lam = 0;
124+ for (let k = 0; k < 3 * n; k++) {
125+ y[k] *= rs[Math.floor(k / 3)];
126+ lam += (y[k] * v[k]) / nv;
127+ }
128+ v.set(y);
129+ }
130+ let w2 = Math.max(lam, 0);
131+ if (mallet) {
132+ let inv = 1 / mallet.mass;
133+ for (let i = 0; i < n; i++) if (!pinned[i] && mallet.weights[i] > 0) inv += (mallet.weights[i] * mallet.weights[i]) / sys.mass[i];
134+ w2 += contactStiffness(mallet) * inv;
59135 }
60136 return Math.sqrt(w2);
61137 }
62138
63-// per hinge: sum of |M| entries of its rows over the triangles that hold it
64-function hingeStiffness(b: Bending): Float64Array {
65- const k = new Float64Array(b.w0.length);
66- for (let t = 0; t < b.tri.length / 3; t++)
67- for (let a = 0; a < 3; a++) {
68- const e = b.tri[3 * t + a];
69- if (e < 0) continue;
70- for (let c = 0; c < 3; c++) k[e] += Math.abs(b.M[9 * t + 3 * a + c]);
71- }
72- return k;
139+export interface TimeStep {
140+ omegaMax: number;
141+ dt: number;
142+ stepsPerSample: number;
143+}
144+
145+// dt = cfl * 2 / omega_max, rounded down to a whole number of steps per audio sample
146+export function timeStep(sys: BeadSystem, mallet: Mallet | undefined, fs: number, cfl: number): TimeStep {
147+ const omegaMax = omegaMaxEstimate(sys, mallet);
148+ const stepsPerSample = Math.max(1, Math.ceil(omegaMax / (fs * cfl * 2)));
149+ return { omegaMax, dt: 1 / (fs * stepsPerSample), stepsPerSample };
73150 }
74151
75152 // bending energy at positions X, optionally only over triangles with keep(t)
@@ -107,12 +184,10 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
107184 const hPhi = new Float64Array(nh);
108185 const hRate = new Float64Array(nh);
109186 const hMom = new Float64Array(nh);
187+ const hMomEl = new Float64Array(nh); // elastic part of the moments
110188
111189 const mallet = exc.mallet;
112- const omegaMax = omegaMaxBound(sys, mallet);
113- const dtMax = (opt.cfl * 2) / omegaMax;
114- const stepsPerSample = Math.max(1, Math.ceil(1 / (opt.fs * dtMax)));
115- const dt = 1 / (opt.fs * stepsPerSample);
190+ const { omegaMax, dt, stepsPerSample } = timeStep(sys, mallet, opt.fs, opt.cfl);
116191 const nSamples = Math.round(opt.duration * opt.fs);
117192 const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
118193 const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
@@ -120,6 +195,7 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
120195 const X = new Float64Array(3 * nn);
121196 const V = new Float64Array(3 * nn);
122197 const F = new Float64Array(3 * n); // force on beads
198+ const Fel = new Float64Array(3 * n); // elastic (position-dependent) part
123199 const A = new Float64Array(3 * n);
124200 for (let i = 0; i < n; i++) {
125201 for (let c = 0; c < 3; c++) {
@@ -226,14 +302,11 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
226302
227303 const hOut = new Float64Array(13);
228304 let eHinge = 0;
229- // Fills F with total force on beads; sets damping power (W, >= 0) and
230- // external power (W), both at velocity V.
231- let pDamp = 0;
232- let pExt = 0;
233- const computeForces = (t: number) => {
234- F.fill(0);
235- contact();
236- pDamp = pFriction;
305+ const anyDashpot = sc.some((c) => c !== 0);
306+ // Position-dependent forces: spring and bending elasticity into Fel; hinge
307+ // angles, gradients and elastic moments; rigid-body geometry.
308+ const elasticForces = () => {
309+ Fel.fill(0);
237310 for (let s = 0; s < ns; s++) {
238311 const a = 3 * si[s];
239312 const b = 3 * sj[s];
@@ -242,49 +315,27 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
242315 const dz = X[b + 2] - X[a + 2];
243316 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
244317 const g = (sk[s] * (l - sL0[s])) / l;
245- let fx = g * dx;
246- let fy = g * dy;
247- let fz = g * dz;
248- const c = sc[s];
249- if (c !== 0) {
250- const rx = V[b] - V[a];
251- const ry = V[b + 1] - V[a + 1];
252- const rz = V[b + 2] - V[a + 2];
253- fx += c * rx;
254- fy += c * ry;
255- fz += c * rz;
256- pDamp += c * (rx * rx + ry * ry + rz * rz);
257- }
258318 if (si[s] < n) {
259- F[a] += fx;
260- F[a + 1] += fy;
261- F[a + 2] += fz;
319+ Fel[a] += g * dx;
320+ Fel[a + 1] += g * dy;
321+ Fel[a + 2] += g * dz;
262322 }
263323 if (sj[s] < n) {
264- F[b] -= fx;
265- F[b + 1] -= fy;
266- F[b + 2] -= fz;
324+ Fel[b] -= g * dx;
325+ Fel[b + 1] -= g * dy;
326+ Fel[b + 2] -= g * dz;
267327 }
268328 }
269329 eHinge = 0;
270330 if (bend) {
271- const s1 = bend.sigma1;
272331 for (let e = 0; e < nh; e++) {
273332 const th = hingeAngleGrad(X, bend.w0[e], bend.w1[e], bend.e0[e], bend.e1[e], hOut);
274333 const sg = bend.sign[e];
275334 hPhi[e] = sg * (th - bend.theta0[e]);
276335 const o = 12 * e;
277336 for (let q = 0; q < 12; q++) hG[o + q] = sg * hOut[1 + q];
278- if (s1 !== 0) {
279- const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
280- hRate[e] =
281- hG[o] * V[w0] + hG[o + 1] * V[w0 + 1] + hG[o + 2] * V[w0 + 2] +
282- hG[o + 3] * V[w1] + hG[o + 4] * V[w1 + 1] + hG[o + 5] * V[w1 + 2] +
283- hG[o + 6] * V[e0] + hG[o + 7] * V[e0 + 1] + hG[o + 8] * V[e0 + 2] +
284- hG[o + 9] * V[e1] + hG[o + 10] * V[e1 + 1] + hG[o + 11] * V[e1 + 2];
285- }
286337 }
287- hMom.fill(0);
338+ hMomEl.fill(0);
288339 const tri = bend.tri;
289340 const Mt = bend.M;
290341 for (let t = 0; t < nt; t++) {
@@ -292,20 +343,76 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
292343 const ea = tri[3 * t + a];
293344 if (ea < 0) continue;
294345 let m = 0;
295- let r = 0;
296346 for (let c = 0; c < 3; c++) {
297347 const ec = tri[3 * t + c];
298- if (ec < 0) continue;
299- const M = Mt[9 * t + 3 * a + c];
300- m += M * hPhi[ec];
301- r += M * hRate[ec];
348+ if (ec >= 0) m += Mt[9 * t + 3 * a + c] * hPhi[ec];
302349 }
303350 eHinge += 0.5 * hPhi[ea] * m;
304- if (s1 !== 0) {
351+ hMomEl[ea] += m;
352+ }
353+ }
354+ }
355+ if (sys.freeRigid) rigidGeometry();
356+ };
357+
358+ // Fills F with total force on beads and A with accelerations, from the
359+ // elastic part (elasticForces must have run at this position) and the
360+ // velocity-dependent part at V; sets damping power (W, >= 0) and external
361+ // power (W), both at V.
362+ let pDamp = 0;
363+ let pExt = 0;
364+ const dynamicForces = (t: number) => {
365+ F.set(Fel);
366+ contact();
367+ pDamp = pFriction;
368+ if (anyDashpot)
369+ for (let s = 0; s < ns; s++) {
370+ const c = sc[s];
371+ if (c === 0) continue;
372+ const a = 3 * si[s];
373+ const b = 3 * sj[s];
374+ const rx = V[b] - V[a];
375+ const ry = V[b + 1] - V[a + 1];
376+ const rz = V[b + 2] - V[a + 2];
377+ pDamp += c * (rx * rx + ry * ry + rz * rz);
378+ if (si[s] < n) {
379+ F[a] += c * rx;
380+ F[a + 1] += c * ry;
381+ F[a + 2] += c * rz;
382+ }
383+ if (sj[s] < n) {
384+ F[b] -= c * rx;
385+ F[b + 1] -= c * ry;
386+ F[b + 2] -= c * rz;
387+ }
388+ }
389+ if (bend) {
390+ const s1 = bend.sigma1;
391+ hMom.set(hMomEl);
392+ if (s1 !== 0) {
393+ for (let e = 0; e < nh; e++) {
394+ const o = 12 * e;
395+ const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
396+ hRate[e] =
397+ hG[o] * V[w0] + hG[o + 1] * V[w0 + 1] + hG[o + 2] * V[w0 + 2] +
398+ hG[o + 3] * V[w1] + hG[o + 4] * V[w1 + 1] + hG[o + 5] * V[w1 + 2] +
399+ hG[o + 6] * V[e0] + hG[o + 7] * V[e0 + 1] + hG[o + 8] * V[e0 + 2] +
400+ hG[o + 9] * V[e1] + hG[o + 10] * V[e1 + 1] + hG[o + 11] * V[e1 + 2];
401+ }
402+ const tri = bend.tri;
403+ const Mt = bend.M;
404+ for (let t = 0; t < nt; t++) {
405+ for (let a = 0; a < 3; a++) {
406+ const ea = tri[3 * t + a];
407+ if (ea < 0) continue;
408+ let r = 0;
409+ for (let c = 0; c < 3; c++) {
410+ const ec = tri[3 * t + c];
411+ if (ec >= 0) r += Mt[9 * t + 3 * a + c] * hRate[ec];
412+ }
305413 pDamp += 2 * s1 * hRate[ea] * r;
306- m += 2 * s1 * r;
414+ hMom[ea] += 2 * s1 * r;
307415 }
308- hMom[ea] += m;
309416 }
310417 }
311418 for (let e = 0; e < nh; e++) {
@@ -359,22 +466,49 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
359466 }
360467 if (sys.freeRigid) removeRigid(A);
361468 };
469+ const computeForces = (t: number) => {
470+ elasticForces();
471+ dynamicForces(t);
472+ };
362473
363- // subtract the rigid-body part (mass-weighted) of a per-bead field u;
364- // returns the kinetic energy of that part when u is a velocity
365- const removeRigid = (u: Float64Array): number => {
366- let M = 0, cx = 0, cy = 0, cz = 0;
474+ // centre of mass and inverse inertia tensor at the current positions
475+ let totalMass = 0;
476+ for (let i = 0; i < n; i++) totalMass += mass[i];
477+ const cm = new Float64Array(3);
478+ const Iinv = new Float64Array(6);
479+ const rigidGeometry = () => {
480+ let cx = 0, cy = 0, cz = 0;
367481 for (let i = 0; i < n; i++) {
368- M += mass[i];
369482 cx += mass[i] * X[3 * i];
370483 cy += mass[i] * X[3 * i + 1];
371484 cz += mass[i] * X[3 * i + 2];
372485 }
373- cx /= M;
374- cy /= M;
375- cz /= M;
376- let px = 0, py = 0, pz = 0, lx = 0, ly = 0, lz = 0;
486+ cx /= totalMass;
487+ cy /= totalMass;
488+ cz /= totalMass;
377489 let ixx = 0, iyy = 0, izz = 0, ixy = 0, ixz = 0, iyz = 0;
490+ for (let i = 0; i < n; i++) {
491+ const m = mass[i];
492+ const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
493+ ixx += m * (ry * ry + rz * rz);
494+ iyy += m * (rx * rx + rz * rz);
495+ izz += m * (rx * rx + ry * ry);
496+ ixy -= m * rx * ry;
497+ ixz -= m * rx * rz;
498+ iyz -= m * ry * rz;
499+ }
500+ cm[0] = cx;
501+ cm[1] = cy;
502+ cm[2] = cz;
503+ invert3(ixx, ixy, ixz, iyy, iyz, izz, Iinv);
504+ };
505+
506+ // subtract the rigid-body part (mass-weighted) of a per-bead field u, about
507+ // the geometry from rigidGeometry; returns the kinetic energy of that part
508+ // when u is a velocity
509+ const removeRigid = (u: Float64Array): number => {
510+ const [cx, cy, cz] = cm;
511+ let px = 0, py = 0, pz = 0, lx = 0, ly = 0, lz = 0;
378512 for (let i = 0; i < n; i++) {
379513 const m = mass[i];
380514 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
@@ -385,15 +519,12 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
385519 lx += m * (ry * uz - rz * uy);
386520 ly += m * (rz * ux - rx * uz);
387521 lz += m * (rx * uy - ry * ux);
388- ixx += m * (ry * ry + rz * rz);
389- iyy += m * (rx * rx + rz * rz);
390- izz += m * (rx * rx + ry * ry);
391- ixy -= m * rx * ry;
392- ixz -= m * rx * rz;
393- iyz -= m * ry * rz;
394522 }
395- const [wx, wy, wz] = solve3(ixx, ixy, ixz, iyy, iyz, izz, lx, ly, lz);
396- const vx = px / M, vy = py / M, vz = pz / M;
523+ const [a, b, c, d, e, f] = Iinv;
524+ const wx = a * lx + b * ly + c * lz;
525+ const wy = b * lx + d * ly + e * lz;
526+ const wz = c * lx + e * ly + f * lz;
527+ const vx = px / totalMass, vy = py / totalMass, vz = pz / totalMass;
397528 for (let i = 0; i < n; i++) {
398529 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
399530 u[3 * i] -= vx + (wy * rz - wz * ry);
@@ -402,10 +533,12 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
402533 }
403534 return 0.5 * (px * vx + py * vy + pz * vz) + 0.5 * (wx * lx + wy * ly + wz * lz);
404535 };
405- if (sys.freeRigid) removeRigid(V); // an initial strike carries no rigid motion either
536+ if (sys.freeRigid) {
537+ rigidGeometry();
538+ removeRigid(V); // an initial strike carries no rigid motion either
539+ }
406540
407- const velocityForces =
408- sys.sigma0 !== 0 || sc.some((c) => c !== 0) || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
541+ const velocityForces = sys.sigma0 !== 0 || anyDashpot || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
409542 let maxStretch = 0;
410543 const malletEnergy = () => (mallet ? 0.5 * mallet.mass * (qv[0] ** 2 + qv[1] ** 2 + qv[2] ** 2) + eContact : 0);
411544 // total mechanical energy at the current state; computeForces must have run
@@ -507,8 +640,8 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
507640 if (sys.freeRigid) dissipated += removeRigid(V);
508641 syncAnchors();
509642 }
510- // at a sample time: recompute with the full-step velocity
511- if (velocityForces) computeForces(t);
643+ // at a sample time: recompute with the full-step velocity (positions unchanged)
644+ if (velocityForces) dynamicForces(t);
512645 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
513646 }
514647
@@ -529,11 +662,20 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog
529662 };
530663 }
531664
532-// solve the symmetric 3x3 system [a b c; b d e; c e f] x = (r0, r1, r2)
533-function solve3(a: number, b: number, c: number, d: number, e: number, f: number, r0: number, r1: number, r2: number): [number, number, number] {
665+// inverse of the symmetric 3x3 matrix [a b c; b d e; c e f], as its six
666+// distinct entries (zero if singular)
667+function invert3(a: number, b: number, c: number, d: number, e: number, f: number, out: Float64Array) {
534668 const A00 = d * f - e * e, A01 = c * e - b * f, A02 = b * e - c * d;
535669 const A11 = a * f - c * c, A12 = b * c - a * e, A22 = a * d - b * b;
536670 const det = a * A00 + b * A01 + c * A02;
537- if (!(Math.abs(det) > 0)) return [0, 0, 0];
538- return [(A00 * r0 + A01 * r1 + A02 * r2) / det, (A01 * r0 + A11 * r1 + A12 * r2) / det, (A02 * r0 + A12 * r1 + A22 * r2) / det];
671+ if (!(Math.abs(det) > 0)) {
672+ out.fill(0);
673+ return;
674+ }
675+ out[0] = A00 / det;
676+ out[1] = A01 / det;
677+ out[2] = A02 / det;
678+ out[3] = A11 / det;
679+ out[4] = A12 / det;
680+ out[5] = A22 / det;
539681 }
src/sim/types.tsmodified+2−0View file
@@ -96,4 +96,6 @@ export interface SimResult {
9696 maxStretch: number; // largest |l - L0| / L0 of bead-to-bead springs, at frame times
9797 contactForce?: Float32Array; // mallet contact force at audio rate (N)
9898 malletPos?: Float32Array; // mallet position q (3 per frame) at frame times (m)
99+ backend?: string; // which solver ran, e.g. 'CPU' or 'WebGPU (intel gen-12lp)'
100+ fallback?: string; // why WebGPU was asked for but not used
99101 }
src/sim/worker.tsmodified+19−7View file
@@ -1,19 +1,31 @@
11 import { simulate } from './solver';
2-import type { BeadSystem, Excitation, SimOptions } from './types';
2+import { simulateGPU } from './gpu';
3+import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
34
45 export interface WorkerRequest {
56 sys: BeadSystem;
67 exc: Excitation;
78 opt: SimOptions;
9+ gpu: boolean;
810 }
911
10-self.onmessage = (ev: MessageEvent<WorkerRequest>) => {
11- const { sys, exc, opt } = ev.data;
12+self.onmessage = async (ev: MessageEvent<WorkerRequest>) => {
13+ const { sys, exc, opt, gpu } = ev.data;
14+ const progress = (f: number) => self.postMessage({ type: 'progress', fraction: f });
1215 try {
13- const result = simulate(sys, exc, opt, (f) => self.postMessage({ type: 'progress', fraction: f }));
14- self.postMessage({ type: 'done', result }, {
15- transfer: [result.audio.buffer, result.frames.buffer, result.frameTimes.buffer, result.energy.buffer, result.balance.buffer],
16- });
16+ let result: SimResult | undefined;
17+ let fallback: string | undefined;
18+ if (gpu) {
19+ try {
20+ const r = await simulateGPU(sys, exc, opt, progress);
21+ result = { ...r, backend: `WebGPU (${r.device})` };
22+ } catch (e) {
23+ fallback = (e as Error).message ?? String(e);
24+ }
25+ }
26+ result ??= { ...simulate(sys, exc, opt, progress), backend: 'CPU', fallback };
27+ const transfer = [result.audio.buffer, result.frames.buffer, result.frameTimes.buffer, result.energy.buffer, result.balance.buffer];
28+ self.postMessage({ type: 'done', result }, { transfer: transfer as ArrayBuffer[] });
1729 } catch (e) {
1830 self.postMessage({ type: 'error', message: String(e) });
1931 }
src/style.cssmodified+1−0View file
@@ -116,6 +116,7 @@ input:focus, select:focus { outline: 2px solid var(--accent); outline-offset: -1
116116 .btn:disabled { opacity: 0.5; cursor: default; }
117117 .btn-small { background: transparent; color: var(--text-primary); border: 1px solid var(--border); padding: 4px 10px; font-size: 13px; }
118118 .run-row { display: flex; gap: 8px; align-items: center; margin-top: 14px; }
119+.gpu-toggle { margin-left: auto; }
119120 progress { width: 100%; margin-top: 8px; accent-color: var(--accent); }
120121
121122 .results { display: flex; flex-direction: column; gap: 16px; min-width: 0; }
src/ui/instrumentPage.tsmodified+30−6View file
@@ -5,7 +5,7 @@ import { spectrum } from './fft';
55 import { encodeWav } from './wav';
66 import { BeadScene, type SceneSpec } from './scene3d';
77 import { Timeline } from './timeline';
8-import { omegaMaxBound } from '../sim/solver';
8+import { timeStep } from '../sim/solver';
99 import { runInWorker, type Running } from '../sim/run';
1010 import type { BeadSystem, Excitation, SimResult } from '../sim/types';
1111
@@ -44,6 +44,8 @@ export interface InstrumentConfig<P extends SimParams> {
4444 audioQuantity?: { label: string; unit: string }; // default: acceleration, m/s²
4545 presets: Preset<P>[];
4646 sections: Section<P>[];
47+ // WebGPU on by default (default true); off where the systems are too small to gain from it
48+ webgpuDefault?: boolean;
4749 // instrument-specific summary of the setup; built = null when the parameters are invalid
4850 setup(p: P): { rows: [string, string][]; warnings: string[]; built: Built | null };
4951 summaryRows?(p: P, built: Built, r: SimResult): [string, string][];
@@ -66,6 +68,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
6668 <div class="run-row">
6769 <button class="btn" data-r="run">Run simulation</button>
6870 <button class="btn-small" data-r="cancel" hidden>Cancel</button>
71+ <label class="inline gpu-toggle" data-r="gpu-label"><input type="checkbox" data-r="gpu"> WebGPU</label>
6972 </div>
7073 <progress data-r="progress" max="1" value="0" hidden></progress>
7174 <p class="error" data-r="error" hidden></p>
@@ -113,6 +116,27 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
113116 const cancelBtn = q<HTMLButtonElement>('cancel');
114117 const progress = q<HTMLProgressElement>('progress');
115118 const errorBox = q<HTMLParagraphElement>('error');
119+ const gpuBox = q<HTMLInputElement>('gpu');
120+ const hasGpu = 'gpu' in navigator;
121+ const gpuKey = `toy-instruments.webgpu.${cfg.wavName}`;
122+ let gpuPref: string | null = null;
123+ try {
124+ gpuPref = localStorage.getItem(gpuKey);
125+ } catch {
126+ /* storage unavailable */
127+ }
128+ gpuBox.checked = hasGpu && (gpuPref ? gpuPref === 'on' : cfg.webgpuDefault !== false);
129+ gpuBox.disabled = !hasGpu;
130+ q('gpu-label').title = hasGpu
131+ ? 'Run the solver on the GPU (single precision). Faster for meshes of a few hundred beads or more; slower for small systems. Falls back to the CPU if the GPU cannot run it.'
132+ : 'This browser does not offer WebGPU; the solver runs on the CPU.';
133+ gpuBox.addEventListener('change', () => {
134+ try {
135+ localStorage.setItem(gpuKey, gpuBox.checked ? 'on' : 'off');
136+ } catch {
137+ /* storage unavailable */
138+ }
139+ });
116140
117141 const wave = new LinePlot(q('wave'));
118142 const spec = new Spectrogram(q('spec'));
@@ -141,9 +165,8 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
141165 const { rows, warnings, built } = cfg.setup(p);
142166 if (built) {
143167 scene.showSetup(built.scene);
144- const wmax = omegaMaxBound(built.sys, built.exc.mallet);
145- const spp = Math.max(1, Math.ceil(1 / (p.fs * ((p.cfl * 2) / wmax))));
146- rows.push(['Highest frequency bound', `${fmtNum(wmax / (2 * Math.PI))} Hz`]);
168+ const { omegaMax: wmax, stepsPerSample: spp } = timeStep(built.sys, built.exc.mallet, p.fs, p.cfl);
169+ rows.push(['Highest frequency', `${fmtNum(wmax / (2 * Math.PI))} Hz`]);
147170 rows.push(['Time step', `${fmtNum(1 / (p.fs * spp))} s (${spp} per sample)`]);
148171 // a bending hinge costs about four springs
149172 const elements = built.sys.si.length + 4 * (built.sys.bending?.w0.length ?? 0);
@@ -151,7 +174,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
151174 rows.push(['Element updates', fmtNum(work)]);
152175 const frameMB = (Math.round(p.duration * p.frameRate) * 3 * built.sys.n * 4) / 1e6;
153176 if (frameMB > 400) warnings.push(`Replay frames would take ${fmtNum(frameMB)} MB; lower the replay frame rate or the duration.`);
154- if (work > 3e9) warnings.push(`This run is large (${fmtNum(work)} element updates): expect roughly ${fmtNum(work / 2.2e7 / 60)} min on a desktop.`);
177+ if (work > 3e9) warnings.push(`This run is large (${fmtNum(work)} element updates): expect roughly ${fmtNum(work / 2.2e7 / 60)} min on a desktop CPU (WebGPU is usually much faster).`);
155178 }
156179 q('derived').innerHTML = statRows(rows);
157180 q('warnings').innerHTML = warnings.map((w) => `<li>${w}</li>`).join('');
@@ -213,6 +236,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
213236 for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
214237 const rt = p.duration / r.wallSeconds;
215238 q('summary').innerHTML = statRows([
239+ ['Solver', (r.backend ?? 'CPU') + (r.fallback ? ` (WebGPU failed: ${r.fallback})` : '')],
216240 ['Time step', `${fmtNum(r.dt)} s, ${r.stepsPerSample} per sample`],
217241 ['Wall time', `${fmtNum(r.wallSeconds)} s (${rt >= 1 ? fmtNum(rt) + '× faster than' : fmtNum(1 / rt) + '× slower than'} real time)`],
218242 ['Energy balance error', emax > 0 ? `${fmtNum((bal / emax) * 100)}% of peak energy` : '–'],
@@ -271,7 +295,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
271295 cancelBtn.hidden = false;
272296 progress.hidden = false;
273297 progress.value = 0;
274- running = runInWorker(built.sys, built.exc, opt, (f) => (progress.value = f));
298+ running = runInWorker(built.sys, built.exc, opt, gpuBox.checked, (f) => (progress.value = f));
275299 try {
276300 show(p, built, await running.promise);
277301 } catch (e) {
tsconfig.jsonmodified+1−1View file
@@ -8,7 +8,7 @@
88 "noUnusedLocals": true,
99 "noUnusedParameters": true,
1010 "skipLibCheck": true,
11- "types": ["vite/client"]
11+ "types": ["vite/client", "@webgpu/types"]
1212 },
1313 "include": ["src", "test"]
1414 }