Run the time loop on the GPU in the web app when WebGPU is available - simulated_instruments/cymbal2

simulated_instruments / cymbal2

Run the time loop on the GPU in the web app when WebGPU is available

Jeremy Magland <jmagland@flatironinstitute.org> committed 1 month ago commit 65ecc0383eaa parent 1793cec Browse files

9 changed files+1066−82

README.mdmodified+1−1View file

@@ -47,6 +47,6 @@ Numerically, the solver departs from the paper in one structural way: instead of

4747

4848 ## Defaults

4949

50-The paper's measured instrument: an 18-inch medium crash, 23 cm radius, 4 cm rise, 1 mm bronze at 8860 kg/m^3, on its measured washer stiffness, struck by the paper's mallet (25 mm, 28 g head, on the measured stick) at 0.6 of the radius. The grid is the paper's, 24 radial intervals by 60 azimuthal points, and the damping constants are the paper's fitted sigma0 = 1.34 s^-1 and sigma1 = 1.2e-3, which put T60 near 5 s at the bottom of the spectrum and near 2 s at 5 kHz; the default tail of 6 s carries the ring to the end. The default stroke (5.75 m/s hand, the tip arriving at 8.2 m/s, 311 N peak contact) was chosen against the authors' published demonstration recording: it matches that recording's band balance to within about 2 dB from 100 Hz to 12 kHz, its high-frequency buildup over the first half second, and its decay envelope, the residual being a few dB of extra energy below 700 Hz, where the pickup hears the washer wobble. The same stroke's peak pickup acceleration, 2.4e4 m/s^2, sits at the scale of the paper's Fig. 7 waveforms. A default render integrates about 580,000 steps at 96 kHz and takes a few minutes in MATLAB or under the numbl CLI; `oversample = 1` halves that at some cost in headroom for hard strokes. The mesh resolves the dynamics reliably to roughly 2 kHz by the paper's own estimate; the crash extends far above that, but its top octaves are dispersion-warped and weak, which is the main respect in which this model is duller than a real crash. Raising `n_radial` and `n_theta` together sharpens it at quadratic cost.

50+The paper's measured instrument: an 18-inch medium crash, 23 cm radius, 4 cm rise, 1 mm bronze at 8860 kg/m^3, on its measured washer stiffness, struck by the paper's mallet (25 mm, 28 g head, on the measured stick) at 0.6 of the radius. The grid is the paper's, 24 radial intervals by 60 azimuthal points, and the damping constants are the paper's fitted sigma0 = 1.34 s^-1 and sigma1 = 1.2e-3, which put T60 near 5 s at the bottom of the spectrum and near 2 s at 5 kHz; the default tail of 6 s carries the ring to the end. The default stroke (5.75 m/s hand, the tip arriving at 8.2 m/s, 311 N peak contact) was chosen against the authors' published demonstration recording: it matches that recording's band balance to within about 2 dB from 100 Hz to 12 kHz, its high-frequency buildup over the first half second, and its decay envelope, the residual being a few dB of extra energy below 700 Hz, where the pickup hears the washer wobble. The same stroke's peak pickup acceleration, 2.4e4 m/s^2, sits at the scale of the paper's Fig. 7 waveforms. A default render integrates about 580,000 steps at 96 kHz and takes a few minutes in MATLAB or under the numbl CLI; `oversample = 1` halves that at some cost in headroom for hard strokes. In a browser with WebGPU the web app runs the same integration on the GPU in well under a minute (see `webapp/README.md`). The mesh resolves the dynamics reliably to roughly 2 kHz by the paper's own estimate; the crash extends far above that, but its top octaves are dispersion-warped and weak, which is the main respect in which this model is duller than a real crash. Raising `n_radial` and `n_theta` together sharpens it at quadratic cost.

5151

5252 `strike_velocity` is the one control that matters most: 0.5 m/s is a tap that stays nearly linear and rings like the bell of a gong, 2 m/s begins to wash, and 4 m/s and above crash, with the high-frequency energy arriving over hundreds of milliseconds after the attack. `strike_count` with `strike_rate_hz` plays a repeated figure; each stroke's contact speed varies with `velocity_spread`. The washer parameters, the grip (`grip_clamped`, `grip_spring`), and the stick and head properties are all exposed and all audible in the attack.

webapp/README.mdmodified+3−1View file

@@ -17,7 +17,9 @@ The self-test runs the closed-form and reproducibility checks and the pinned dia

1717

1818 ## Where the time goes

1919

20-Everywhere. Unlike the sibling `cymbal` app, whose expensive stage (the eigensolves) could be cached and whose synthesis was an oscillator bank, here the sound is a nonlinear time-domain integration of the whole shell, several hundred thousand steps of coupled solves, and every parameter feeds it. There is therefore no debounced auto-render: the Render preview button runs a cheap pass (no oversampling, tail clamped to 1.5 s; the azimuthal grid is kept, because hard strokes diverge without the full wavenumber ladder), Render full runs the real thing, and both report elapsed time and can be cancelled. Expect tens of seconds for a preview and minutes for a full render in the browser; the only piece the worker caches is the cheap linear mode table that the spectrum plot overlays.

20+Everywhere. Unlike the sibling `cymbal` app, whose expensive stage (the eigensolves) could be cached and whose synthesis was an oscillator bank, here the sound is a nonlinear time-domain integration of the whole shell, several hundred thousand steps of coupled solves, and every parameter feeds it. There is therefore no debounced auto-render: the Render preview button runs a cheap pass (no oversampling, tail clamped to 1.5 s; the azimuthal grid is kept, because hard strokes diverge without the full wavenumber ladder), Render full runs the real thing, and both report elapsed time and can be cancelled. The only piece the worker caches is the cheap linear mode table that the spectrum plot overlays.

21+

22+When the browser exposes WebGPU, the time loop runs on the GPU instead of through the interpreter: the worker still executes the unmodified setup stage (grid, per-wavenumber inverses, mallet, stroke schedule) and output stage through numbl, but hands the operators to a WGSL port of `cymbal2_sim.m` (`src/worker/gpuSim.ts`, eight compute dispatches per time step; the design and its measurements are documented in the [cymbal2-webgpu](https://vault1.magland.org/proof-of-concept/cymbal2-webgpu) exploration). This is 20x or more faster than interpreting the loop, even on an integrated GPU, and the diagnostics panel says which path a render took. The GPU integrates in single precision; the contact diagnostics match the double-precision interpreter to a fraction of a percent, and the trajectory of a crash is chaotic either way. Without WebGPU the interpreter path runs as before, and expect tens of seconds for a preview and many minutes for a full render. One interaction with the code editor: the WGSL mirrors the shipped `cymbal2_sim.m`, so editing that one file falls back to the interpreter, where the edit is live (a warning in the diagnostics says so); edits to every other `.m` file flow through numbl on both paths. `npm i --no-save webgpu && npx tsx scripts/gpu-selftest.ts` runs the GPU loop from Node on desktop WebGPU and checks the pinned diagnostics.

2123

2224 ## Reading and editing the code in the app

2325

webapp/package-lock.jsonmodified+8−0View file

@@ -23,6 +23,7 @@

2323 "@types/react": "^18.3.5",

2424 "@types/react-dom": "^18.3.0",

2525 "@vitejs/plugin-react": "^4.3.1",

26+ "@webgpu/types": "^0.1.72",

2627 "tsx": "^4.19.0",

2728 "typescript": "^5.5.4",

2829 "vite": "^5.4.3"

@@ -1517,6 +1518,13 @@

15171518 "vite": "^4.2.0 || ^5.0.0 || ^6.0.0 || ^7.0.0"

15181519 }

15191520 },

1521+ "node_modules/@webgpu/types": {

1522+ "version": "0.1.72",

1523+ "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.72.tgz",

1524+ "integrity": "sha512-0cF7RFM2edNoiIS1ODJp0/Gzv4/xSXhwoR0YCza+OWpJWtn4wmo9DvK91aLlH9+uUnwIriP7ZiC3WitmyhuzBw==",

1525+ "dev": true,

1526+ "license": "BSD-3-Clause"

1527+ },

15201528 "node_modules/bail": {

15211529 "version": "2.0.2",

15221530 "resolved": "https://registry.npmjs.org/bail/-/bail-2.0.2.tgz",

webapp/package.jsonmodified+1−0View file

@@ -25,6 +25,7 @@

2525 "@types/react": "^18.3.5",

2626 "@types/react-dom": "^18.3.0",

2727 "@vitejs/plugin-react": "^4.3.1",

28+ "@webgpu/types": "^0.1.72",

2829 "tsx": "^4.19.0",

2930 "typescript": "^5.5.4",

3031 "vite": "^5.4.3"

webapp/scripts/gpu-selftest.tsadded+169−0View file

@@ -0,0 +1,169 @@

1+// Runs the webapp's WebGPU time loop (src/worker/gpuSim.ts) from Node against

2+// desktop WebGPU (the optional `webgpu` package, prebuilt Google Dawn) on the

3+// validation suite's reduced configuration, and checks the pinned contact

4+// diagnostics and output rms. This exercises exactly the code the browser

5+// worker runs; a machine without the package or without a GPU reports and

6+// exits 0 so CI can run it opportunistically.

7+// npm i --no-save webgpu && npx tsx scripts/gpu-selftest.ts

8+

9+import { readFileSync, readdirSync } from 'node:fs';

10+import { fileURLToPath } from 'node:url';

11+import path from 'node:path';

12+import { executeCode, RTV, type RuntimeValue, type RuntimeTensor } from 'numbl';

13+

14+const MATLAB_DIR = path.resolve(path.dirname(fileURLToPath(import.meta.url)), '../../matlab');

15+

16+async function installDawn(): Promise<boolean> {

17+ try {

18+ const mod = await import('webgpu');

19+ Object.assign(globalThis, mod.globals);

20+ Object.defineProperty(globalThis, 'navigator', {

21+ value: { gpu: mod.create([]) },

22+ configurable: true,

23+ });

24+ return true;

25+ } catch {

26+ return false;

27+ }

28+}

29+

30+function mfiles() {

31+ return readdirSync(MATLAB_DIR)

32+ .filter((f) => f.endsWith('.m'))

33+ .sort()

34+ .map((name) => ({ name, source: readFileSync(path.join(MATLAB_DIR, name), 'utf8') }));

35+}

36+

37+const scalar = (v: RuntimeValue): number => {

38+ if (typeof v === 'number') return v;

39+ const t = v as RuntimeTensor;

40+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];

41+ throw new Error('expected scalar');

42+};

43+const f32 = (v: RuntimeValue): Float32Array => {

44+ const t = v as RuntimeTensor;

45+ return Float32Array.from(t.data);

46+};

47+

48+async function main() {

49+ if (!(await installDawn())) {

50+ console.log('gpu-selftest: `webgpu` package not installed; skipping (npm i --no-save webgpu)');

51+ return;

52+ }

53+ const { requestGpuDevice, runGpuSim } = await import('../src/worker/gpuSim');

54+

55+ const files = mfiles();

56+ const P = `

57+P = cymbal2_defaults();

58+P.sample_rate = 44100;

59+P.oversample = 1;

60+P.tail_s = 1.5;

61+`;

62+ const r = executeCode(`${P}

63+S = cymbal2_system(P);

64+MB = cymbal2_mallet(P);

65+G = S.G;

66+o_D1 = G.D1; o_D2 = G.D2; o_F = G.F; o_Gi = G.Gi; o_Rd = G.Rd; o_R2d = G.R2d;

67+o_AiB = S.AiB; o_AiC = S.AiC; o_Ai = S.Ai; o_Fi = S.Fi;

68+o_ASiB = MB.ASiB; o_ASiC = MB.ASiC;

69+o_vM = MB.vM; o_vR = MB.vR; o_vK = MB.vK; o_tipcol = MB.tipcol;

70+K = max(1, round(P.strike_count));

71+tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);

72+vstr = zeros(1, K);

73+state = P.seed * 2654435761;

74+for i = 1:K

75+ state = mod(state * 1664525 + 1013904223, 4294967296);

76+ x = state / 4294967296;

77+ vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));

78+end

79+gap = P.strike_gap_m;

80+if gap <= 0

81+ gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;

82+end

83+travel = gap ./ vstr + vstr / (2 * P.strike_accel);

84+tlaunch = tstr - travel;

85+off = -min(0, min(tlaunch));

86+tlaunch = tlaunch + off;

87+duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;

88+nsteps = floor(duration * P.sample_rate * P.oversample);

89+o_tlaunch = tlaunch; o_vstr = vstr;

90+c_nrt = G.nrt; c_Nr = G.Nr; c_nth = G.nth; c_hr = G.hr; c_hth = G.hth;

91+c_ht = S.ht; c_Eh = S.Eh; c_Rs = S.Rs; c_rhoh = S.rhoh;

92+c_fsc = S.forceScale; c_nlr0 = S.eomRow0;

93+c_ls = S.ls; c_js = S.js; c_lp = S.lp; c_jp = S.jp;

94+c_tipi = MB.tipi; c_nfree = MB.nfree; c_ps = MB.ps;

95+c_KM = P.head_stiffness; c_alM = P.head_alpha;

96+c_aH = P.strike_accel; c_aHr = P.retract_accel;

97+c_holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / S.ht));

98+c_missDepth = gap + 0.06;

99+c_gap = gap; c_K = K; c_nsteps = nsteps;

100+`, { displayResults: false }, files, 'stage.m');

101+ const v = r.variableValues;

102+ const s = (n: string) => scalar(v[n]);

103+ const t = (n: string) => f32(v[n]);

104+ const setup = {

105+ ops: {

106+ D1: t('o_D1'), D2: t('o_D2'), F: t('o_F'), Gi: t('o_Gi'),

107+ Rd: t('o_Rd'), R2d: t('o_R2d'),

108+ AiB: t('o_AiB'), AiC: t('o_AiC'), Ai: t('o_Ai'), Fi: t('o_Fi'),

109+ ASiB: t('o_ASiB'), ASiC: t('o_ASiC'),

110+ vM: t('o_vM'), vR: t('o_vR'), vK: t('o_vK'), tipcol: t('o_tipcol'),

111+ },

112+ tlaunch: Array.from((v.o_tlaunch as RuntimeTensor).data),

113+ vstr: Array.from((v.o_vstr as RuntimeTensor).data),

114+ nrt: s('c_nrt'), Nr: s('c_Nr'), nth: s('c_nth'),

115+ hr: s('c_hr'), hth: s('c_hth'),

116+ ht: s('c_ht'), Eh: s('c_Eh'), Rs: s('c_Rs'), rhoh: s('c_rhoh'),

117+ fsc: s('c_fsc'), nlr0: s('c_nlr0'),

118+ ls: s('c_ls'), js: s('c_js'), lp: s('c_lp'), jp: s('c_jp'),

119+ tipi: s('c_tipi'), nfree: s('c_nfree'), ps: s('c_ps'),

120+ KM: s('c_KM'), alM: s('c_alM'), aH: s('c_aH'), aHr: s('c_aHr'),

121+ holdSteps: s('c_holdSteps'), missDepth: s('c_missDepth'),

122+ gap: s('c_gap'), K: s('c_K'), nsteps: s('c_nsteps'),

123+ };

124+

125+ let device;

126+ try {

127+ device = await requestGpuDevice();

128+ } catch (e) {

129+ console.log(`gpu-selftest: no WebGPU device (${e instanceof Error ? e.message : e}); skipping`);

130+ return;

131+ }

132+ const t0 = performance.now();

133+ const sim = await runGpuSim(device, setup, () => {});

134+ device.destroy();

135+ const ms = performance.now() - t0;

136+

137+ const ro = executeCode(`${P}

138+[y, peakAccel] = cymbal2_output(P, acc);

139+`, { displayResults: false, initialVariableValues: { acc: RTV.tensor(Float64Array.from(sim.acc), [1, sim.acc.length]) } }, files, 'stage.m');

140+ const y = f32(ro.variableValues.y);

141+ let ss = 0;

142+ for (let i = 0; i < y.length; i++) ss += y[i] * y[i];

143+ const rms = Math.sqrt(ss / y.length);

144+ let peakF = 0;

145+ for (let i = 0; i < sim.fhist.length; i++) peakF = Math.max(peakF, sim.fhist[i]);

146+ const firstC = sim.fhist.findIndex((x) => x > 0);

147+

148+ console.log(

149+ `gpu-selftest: ${setup.nsteps} steps in ${(ms / 1000).toFixed(1)} s ` +

150+ `(${((1000 * ms) / setup.nsteps).toFixed(0)} us/step); ` +

151+ `peakF ${peakF.toFixed(2)} N (ref 311.14), cvel ${sim.contactVel[0].toFixed(3)} (ref -8.223), ` +

152+ `first contact step ${firstC + 1} (ref 1376), rms ${rms.toFixed(4)} (ref 0.0973)`,

153+ );

154+ const ok =

155+ Math.abs(peakF - 311.14) < 3 &&

156+ Math.abs(sim.contactVel[0] + 8.223) < 0.08 &&

157+ firstC + 1 === 1376 &&

158+ Math.abs(rms - 0.0973) < 0.012;

159+ if (!ok) {

160+ console.error('gpu-selftest: FAILED');

161+ process.exit(1);

162+ }

163+ console.log('gpu-selftest: OK');

164+}

165+

166+main().then(() => process.exit(0)).catch((e) => {

167+ console.error(e);

168+ process.exit(1);

169+});

webapp/src/components/DiagnosticsPanel.tsxmodified+2−1View file

@@ -51,7 +51,8 @@ export default function DiagnosticsPanel({ result }: { result: RenderResponse })

5151 <th>Length</th>

5252 <td>

5353 {d.durationS.toFixed(2)} s · {result.quality} · rendered in{' '}

54- {(d.renderMs / 1000).toFixed(1)} s

54+ {(d.renderMs / 1000).toFixed(1)} s{' '}

55+ {d.backend === 'webgpu' ? 'on the GPU' : 'on the CPU'}

5556 </td>

5657 </tr>

5758 </tbody>

webapp/src/hooks/useRenderer.tsmodified+4−1View file

@@ -40,7 +40,10 @@ export function useRenderer(params: Cymbal2Params, mfiles: MFile[]) {

4040 worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage | ErrorMessage>) => {

4141 const msg = ev.data;

4242 if (msg.type === 'progress') {

43- setState((s) => ({ ...s, stage: msg.stage }));

43+ const pct = msg.stage === 'integrate' && msg.frac > 0.1

44+ ? ` ${Math.round(((msg.frac - 0.1) / 0.85) * 100)}%`

45+ : '';

46+ setState((s) => ({ ...s, stage: msg.stage + pct }));

4447 return;

4548 }

4649 inFlight.current = null;

webapp/src/worker/gpuSim.tsadded+559−0View file

@@ -0,0 +1,559 @@

1+/// <reference types="@webgpu/types" />

2+// WebGPU implementation of the cymbal2 time loop (cymbal2_sim.m), ported from

3+// the prototype in proof-of-concept/cymbal2-webgpu. Everything around the

4+// loop stays in numbl: the caller runs the unmodified MATLAB setup stage

5+// (grid, per-wavenumber inverses, mallet, stroke schedule) and hands the

6+// tensors in through GpuSetup; the acceleration and contact-force histories

7+// come back for the unmodified output stage.

8+//

9+// One time step is eight compute dispatches at full GPU width, the dispatch

10+// boundaries providing global synchronization between stages:

11+// AB hand state machine + Hertz contact + beam update (workgroup 0), and

12+// the shell derivative stencils / Airy right-hand side (all workgroups)

13+// C azimuthal analysis of the Airy right-hand side

14+// D Airy solve, one small radial system per wavenumber

15+// E azimuthal synthesis of the stress function

16+// FG von Karman bracket + forcing (stress-function stencils inline)

17+// H azimuthal analysis of the nonlinear load

18+// I implicit shell update, three radial products per wavenumber

19+// J azimuthal synthesis of the new displacement; pickup acceleration

20+//

21+// Every recorded dispatch is identical: the step index k is a chunk base

22+// (state cell 0, written by the host once per chunk) plus a within-chunk

23+// offset from a dynamic-offset uniform, and all slot rotation is a pure

24+// function of k. Slots rotate every step whether or not the shell is active,

25+// which is exact: before first contact every slot is zero, and after first

26+// contact every step writes. The shell-active flag latches inside AB and is

27+// read by the later stages of the same step one dispatch later; the benign

28+// race inside AB itself can only skip work on a field that is still zero.

29+

30+export interface GpuSetup {

31+ ops: {

32+ D1: Float32Array; D2: Float32Array;

33+ F: Float32Array; Gi: Float32Array;

34+ Rd: Float32Array; R2d: Float32Array;

35+ AiB: Float32Array; AiC: Float32Array; Ai: Float32Array; Fi: Float32Array;

36+ ASiB: Float32Array; ASiC: Float32Array;

37+ vM: Float32Array; vR: Float32Array; vK: Float32Array; tipcol: Float32Array;

38+ };

39+ tlaunch: number[];

40+ vstr: number[];

41+ nrt: number; Nr: number; nth: number;

42+ hr: number; hth: number;

43+ ht: number; Eh: number; Rs: number; rhoh: number;

44+ fsc: number; nlr0: number;

45+ ls: number; js: number; lp: number; jp: number;

46+ tipi: number; nfree: number; ps: number;

47+ KM: number; alM: number; aH: number; aHr: number;

48+ holdSteps: number; missDepth: number; gap: number;

49+ K: number; nsteps: number;

50+}

51+

52+export interface GpuSimResult {

53+ acc: Float32Array;

54+ fhist: Float32Array;

55+ peakXi: number;

56+ peakDisp: number;

57+ contactVel: number[];

58+ strokesPlayed: number;

59+ gpuMs: number;

60+}

61+

62+const fmt = (x: number): string => {

63+ if (!Number.isFinite(x)) throw new Error('non-finite kernel constant');

64+ let s = x.toPrecision(9);

65+ if (!/[.e]/.test(s)) s += '.0';

66+ return s + 'f';

67+};

68+

69+interface Layout {

70+ N: number;

71+ opsData: Float32Array<ArrayBuffer>;

72+ stLen: number;

73+ SC_CVEL: number;

74+}

75+

76+function buildLayoutAndStages(su: GpuSetup, wg: number) {

77+ const nrt = su.nrt | 0;

78+ const nth = su.nth | 0;

79+ const nf = su.nfree | 0;

80+ const N = nrt * nth;

81+ const K = su.K | 0;

82+ const ht = su.ht;

83+

84+ // ---- packed operator buffer ----

85+ const parts: [string, Float32Array, number][] = [

86+ ['D1', su.ops.D1, nrt * nrt],

87+ ['D2', su.ops.D2, nrt * nrt],

88+ ['F', su.ops.F, nth * nth],

89+ ['GI', su.ops.Gi, nth * nth],

90+ ['RD', su.ops.Rd, nrt],

91+ ['R2D', su.ops.R2d, nrt],

92+ ['AIB', su.ops.AiB, nrt * nrt * nth],

93+ ['AIC', su.ops.AiC, nrt * nrt * nth],

94+ ['AI', su.ops.Ai, nrt * nrt * nth],

95+ ['FI', su.ops.Fi, nrt * nrt * nth],

96+ ['ASIB', su.ops.ASiB, nf * nf],

97+ ['ASIC', su.ops.ASiC, nf * nf],

98+ ['VM', su.ops.vM, nf],

99+ ['VR', su.ops.vR, nf],

100+ ['VK', su.ops.vK, nf],

101+ ['TIPCOL', su.ops.tipcol, nf],

102+ ];

103+ // Launch-step thresholds: smallest 0-based kk with (kk+1)*ht >= tlaunch,

104+ // decided in f64 exactly as MATLAB's `k*ht >= tlaunch` decides it.

105+ const kgo = su.tlaunch.map((tl) => {

106+ let g = Math.max(0, Math.ceil(tl / ht) - 1);

107+ while (g > 0 && g * ht >= tl) g -= 1;

108+ while ((g + 1) * ht < tl) g += 1;

109+ return g;

110+ });

111+ parts.push(['KGO', Float32Array.from(kgo), K]);

112+ parts.push(['VSTR', Float32Array.from(su.vstr), K]);

113+ // Transposed beam matrices so each beam-row thread streams its own row.

114+ const tr = (m: Float32Array) => {

115+ const t = new Float32Array(nf * nf);

116+ for (let r = 0; r < nf; r++) for (let p = 0; p < nf; p++) t[p + r * nf] = m[r + p * nf];

117+ return t;

118+ };

119+ parts.push(['ASIBT', tr(su.ops.ASiB), nf * nf]);

120+ parts.push(['ASICT', tr(su.ops.ASiC), nf * nf]);

121+

122+ const off: Record<string, number> = {};

123+ let cursor = 0;

124+ for (const [name, arr, len] of parts) {

125+ if (arr.length !== len) throw new Error(`ops part ${name}: ${arr.length} != ${len}`);

126+ off[name] = cursor;

127+ cursor += len;

128+ }

129+ const opsData = new Float32Array(cursor);

130+ for (const [name, arr] of parts) opsData.set(arr, off[name]);

131+

132+ // ---- state buffer: 16 scalars, cvel, 3 beam slots, 3+3 field slots, scratch ----

133+ const ST_CVEL = 16;

134+ const ST_D = ST_CVEL + K;

135+ const ST_U = ST_D + 3 * nf;

136+ const ST_UH = ST_U + 3 * N;

137+ const scratchNames = ['URR', 'V1', 'V2', 'RHSA', 'RHSAF', 'PH3', 'PHI', 'NL', 'NLF'];

138+ const A: Record<string, number> = {};

139+ let scur = ST_UH + 3 * N;

140+ for (const n of scratchNames) {

141+ A[n] = scur;

142+ scur += N;

143+ }

144+ const stLen = scur;

145+

146+ const c: Record<string, number> = {

147+ HT: ht,

148+ HT2I: 1 / (ht * ht),

149+ GAP: su.gap,

150+ KM: su.KM,

151+ ALM: su.alM,

152+ PS: su.ps,

153+ AH: su.aH,

154+ AHR: su.aHr,

155+ MISSD: su.missDepth,

156+ EHRS: su.Eh / su.Rs,

157+ EH2: su.Eh / 2,

158+ RSI: 1 / su.Rs,

159+ SNL: (ht * ht) / su.rhoh,

160+ FSC: su.fsc,

161+ C1R: 1 / (2 * su.hr),

162+ C2R: 1 / (su.hr * su.hr),

163+ C1T: 1 / (2 * su.hth),

164+ C2T: 1 / (su.hth * su.hth),

165+ };

166+

167+ const header = `

168+const NRT = ${nrt}u;

169+const NTH = ${nth}u;

170+const N = ${N}u;

171+const NF = ${nf}u;

172+const NSTEPS = ${su.nsteps}u;

173+const KS = ${K}u;

174+const NRU = ${su.Nr}u;

175+const LS0 = ${su.ls - 1}u; const JS0 = ${su.js - 1}u;

176+const LP0 = ${su.lp - 1}u; const JP0 = ${su.jp - 1}u;

177+const TIPI = ${su.tipi - 1}u;

178+const NLR0 = ${su.nlr0 - 1}u;

179+const HOLDSTEPS = ${su.holdSteps | 0};

180+${Object.entries(c).map(([k, v]) => `const ${k}: f32 = ${fmt(v)};`).join('\n')}

181+${Object.entries(off).map(([k, v]) => `const OP_${k} = ${v}u;`).join('\n')}

182+const ST_CVEL = ${ST_CVEL}u;

183+const ST_D = ${ST_D}u;

184+const ST_U = ${ST_U}u;

185+const ST_UH = ${ST_UH}u;

186+${scratchNames.map((n) => `const A_${n} = ${A[n]}u;`).join('\n')}

187+

188+@group(0) @binding(0) var<storage, read> ops: array<f32>;

189+@group(0) @binding(1) var<storage, read_write> st: array<f32>;

190+@group(0) @binding(2) var<storage, read_write> outb: array<f32>;

191+@group(0) @binding(3) var<uniform> ctl: vec4<u32>;

192+

193+fn thm(j: u32) -> u32 { return (j + NTH - 1u) % NTH; }

194+fn thp(j: u32) -> u32 { return (j + 1u) % NTH; }

195+`;

196+

197+ const AB = `${header}

198+var<workgroup> wgF: f32;

199+var<workgroup> wgHA: f32;

200+var<workgroup> wgHB: f32;

201+var<workgroup> wgHC: f32;

202+var<workgroup> shD: array<f32, ${nf}>;

203+var<workgroup> shD1: array<f32, ${nf}>;

204+

205+fn u_dr(b: u32, l: u32, j: u32) -> f32 {

206+ if (l < 1u || l > NRT - 2u) { return 0.0; }

207+ return (st[b + (l + 1u) + j * NRT] - st[b + (l - 1u) + j * NRT]) * C1R;

208+}

209+

210+@compute @workgroup_size(${wg})

211+fn main(@builtin(global_invocation_id) gid: vec3<u32>,

212+ @builtin(local_invocation_index) lid: u32,

213+ @builtin(workgroup_id) wgid: vec3<u32>) {

214+ // No early-out on k: the host records exactly nsteps steps, and the

215+ // workgroupBarrier below requires uniform control flow up to it.

216+ let k = u32(st[0u]) + ctl.x;

217+ let uS = ST_U + (k % 3u) * N;

218+ if (wgid.x == 0u) {

219+ if (lid == 0u) {

220+ var shellA = u32(st[1u]); var inC = u32(st[2u]); var phase = u32(st[3u]);

221+ var strokeIdx = u32(st[4u]); var holdLeft = i32(st[5u]); var vtar = st[6u];

222+ var vH = st[7u]; var uHm = st[8u]; var uH0 = st[9u]; var uHp = st[10u];

223+ var peakXi = st[12u]; var peakDisp = st[13u];

224+ let dCur = ST_D + ((k + 2u) % 3u) * NF;

225+ let dOld = ST_D + ((k + 1u) % 3u) * NF;

226+ if (shellA != 0u && k > 0u && (k % 512u) == 0u) {

227+ var pd: f32 = 0.0;

228+ for (var l: u32 = 0u; l < NRT; l = l + 1u) {

229+ pd = max(pd, abs(st[uS + l + JS0 * NRT]));

230+ }

231+ if (pd > peakDisp) { peakDisp = pd; }

232+ }

233+ let xi = st[uS + LS0 + JS0 * NRT] - st[dCur + TIPI] - GAP;

234+ var fM: f32 = 0.0;

235+ if (xi > 0.0) {

236+ fM = KM * pow(xi, ALM);

237+ shellA = 1u;

238+ if (inC == 0u && strokeIdx < KS && st[ST_CVEL + strokeIdx] == 0.0) {

239+ st[ST_CVEL + strokeIdx] = (st[dCur + TIPI] - st[dOld + TIPI]) / HT;

240+ }

241+ inC = 1u;

242+ if (xi > peakXi) { peakXi = xi; }

243+ } else {

244+ inC = 0u;

245+ }

246+ outb[NSTEPS + k] = fM;

247+ wgF = fM;

248+ wgHA = -(uHp - 2.0 * uH0 + uHm);

249+ wgHB = -(HT * 0.5) * (uHp - uHm);

250+ wgHC = -HT * HT * (PS * uH0 + ((1.0 - PS) * 0.5) * (uHp + uHm));

251+ if (phase == 0u) {

252+ if (strokeIdx < KS && f32(k) >= ops[OP_KGO + strokeIdx]) {

253+ phase = 1u; vtar = ops[OP_VSTR + strokeIdx];

254+ }

255+ } else if (phase == 1u) {

256+ vH = max(vH - AH * HT, -vtar);

257+ if (inC != 0u) { phase = 3u; holdLeft = HOLDSTEPS; }

258+ else if (vH <= -vtar) { phase = 2u; }

259+ } else if (phase == 2u) {

260+ if (inC != 0u) { phase = 3u; holdLeft = HOLDSTEPS; }

261+ else if (uHp <= -MISSD) { phase = 4u; }

262+ } else if (phase == 3u) {

263+ holdLeft = holdLeft - 1;

264+ if (holdLeft <= 0) { phase = 4u; }

265+ } else if (phase == 4u) {

266+ vH = min(vH + AHR * HT, 0.5 * vtar);

267+ if (vH > 0.0 && uHp >= -vH * vH / (2.0 * AHR)) { phase = 5u; }

268+ } else {

269+ vH = vH - AHR * HT;

270+ if (vH <= 0.0 || uHp >= 0.0) { vH = 0.0; phase = 0u; strokeIdx = strokeIdx + 1u; }

271+ }

272+ uHm = uH0; uH0 = uHp; uHp = uHp + vH * HT;

273+ st[1u] = f32(shellA); st[2u] = f32(inC); st[3u] = f32(phase);

274+ st[4u] = f32(strokeIdx); st[5u] = f32(holdLeft); st[6u] = vtar;

275+ st[7u] = vH; st[8u] = uHm; st[9u] = uH0; st[10u] = uHp;

276+ st[12u] = peakXi; st[13u] = peakDisp; st[14u] = fM;

277+ }

278+ if (lid < NF) {

279+ shD[lid] = st[ST_D + ((k + 2u) % 3u) * NF + lid];

280+ shD1[lid] = st[ST_D + ((k + 1u) % 3u) * NF + lid];

281+ }

282+ workgroupBarrier();

283+ if (lid < NF) {

284+ let r = lid;

285+ var sum: f32 = ops[OP_VM + r] * wgHA + ops[OP_VR + r] * wgHB + ops[OP_VK + r] * wgHC

286+ + (HT * HT * wgF) * ops[OP_TIPCOL + r];

287+ for (var p: u32 = 0u; p < NF; p = p + 1u) {

288+ sum = sum + ops[OP_ASIBT + p + r * NF] * shD[p]

289+ + ops[OP_ASICT + p + r * NF] * shD1[p];

290+ }

291+ st[ST_D + (k % 3u) * NF + r] = sum;

292+ }

293+ }

294+ let i = gid.x;

295+ if (i >= N) { return; }

296+ if (st[1u] == 0.0) { return; }

297+ let l = i % NRT; let j = i / NRT;

298+ var ur: f32 = 0.0; var urr: f32 = 0.0;

299+ if (l >= 1u && l <= NRT - 2u) {

300+ let um = st[uS + (l - 1u) + j * NRT];

301+ let u0 = st[uS + l + j * NRT];

302+ let up = st[uS + (l + 1u) + j * NRT];

303+ ur = (up - um) * C1R;

304+ urr = (um - 2.0 * u0 + up) * C2R;

305+ }

306+ let ujm = st[uS + l + thm(j) * NRT];

307+ let uj0 = st[uS + l + j * NRT];

308+ let ujp = st[uS + l + thp(j) * NRT];

309+ let v1 = ops[OP_RD + l] * ur + ops[OP_R2D + l] * ((ujm - 2.0 * uj0 + ujp) * C2T);

310+ let urt = (u_dr(uS, l, thp(j)) - u_dr(uS, l, thm(j))) * C1T;

311+ let v2 = ops[OP_RD + l] * urt - ops[OP_R2D + l] * ((ujp - ujm) * C1T);

312+ st[A_URR + i] = urr;

313+ st[A_V1 + i] = v1;

314+ st[A_V2 + i] = v2;

315+ var rv: f32 = 0.0;

316+ if (l >= 3u && l <= NRU - 2u) {

317+ rv = EHRS * (urr + v1) - EH2 * 2.0 * (urr * v1 - v2 * v2);

318+ }

319+ st[A_RHSA + i] = rv;

320+}

321+`;

322+

323+ const gate = `

324+ let k = u32(st[0u]) + ctl.x;

325+ if (k >= NSTEPS) { return; }

326+ let i = gid.x;

327+ if (i >= N) { return; }

328+ if (st[1u] == 0.0) { return; }

329+ let l = i % NRT; let j = i / NRT;`;

330+

331+ const transform = (src: string, dst: string, mat: string) => `${header}

332+@compute @workgroup_size(${wg})

333+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {

334+${gate}

335+ var sum: f32 = 0.0;

336+ for (var q: u32 = 0u; q < NTH; q = q + 1u) {

337+ sum = sum + st[${src} + l + q * NRT] * ops[${mat} + q + j * NTH];

338+ }

339+ st[${dst} + i] = sum;

340+}

341+`;

342+

343+ const D = `${header}

344+@compute @workgroup_size(${wg})

345+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {

346+${gate}

347+ let base = OP_FI + j * NRT * NRT;

348+ var sum: f32 = 0.0;

349+ for (var p: u32 = 0u; p < NRT; p = p + 1u) {

350+ sum = sum + ops[base + l + p * NRT] * st[A_RHSAF + p + j * NRT];

351+ }

352+ st[A_PH3 + i] = sum;

353+}

354+`;

355+

356+ const FG = `${header}

357+fn p_dr(l: u32, j: u32) -> f32 {

358+ if (l < 1u || l > NRT - 2u) { return 0.0; }

359+ return (st[A_PHI + (l + 1u) + j * NRT] - st[A_PHI + (l - 1u) + j * NRT]) * C1R;

360+}

361+@compute @workgroup_size(${wg})

362+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {

363+${gate}

364+ var nl: f32 = 0.0;

365+ if (l >= NLR0 && l < NRU) {

366+ let pm = st[A_PHI + (l - 1u) + j * NRT];

367+ let p0 = st[A_PHI + l + j * NRT];

368+ let pp = st[A_PHI + (l + 1u) + j * NRT];

369+ let pr = (pp - pm) * C1R;

370+ let prr = (pm - 2.0 * p0 + pp) * C2R;

371+ let pjm = st[A_PHI + l + thm(j) * NRT];

372+ let pjp = st[A_PHI + l + thp(j) * NRT];

373+ let w1 = ops[OP_RD + l] * pr + ops[OP_R2D + l] * ((pjm - 2.0 * p0 + pjp) * C2T);

374+ let prt = (p_dr(l, thp(j)) - p_dr(l, thm(j))) * C1T;

375+ let w2 = ops[OP_RD + l] * prt - ops[OP_R2D + l] * ((pjp - pjm) * C1T);

376+ nl = prr * st[A_V1 + i] + st[A_URR + i] * w1 - 2.0 * w2 * st[A_V2 + i]

377+ - (prr + w1) * RSI;

378+ if (l == LS0 && j == JS0) { nl = nl - st[14u] * FSC; }

379+ nl = nl * SNL;

380+ }

381+ st[A_NL + i] = nl;

382+}

383+`;

384+

385+ const I = `${header}

386+@compute @workgroup_size(${wg})

387+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {

388+${gate}

389+ let h3S = ST_UH + (k % 3u) * N;

390+ let h13S = ST_UH + ((k + 2u) % 3u) * N;

391+ let h3nS = ST_UH + ((k + 1u) % 3u) * N;

392+ let base = j * NRT * NRT;

393+ var sum: f32 = 0.0;

394+ for (var p: u32 = 0u; p < NRT; p = p + 1u) {

395+ let pj = p + j * NRT;

396+ sum = sum + ops[OP_AIB + base + l + p * NRT] * st[h3S + pj]

397+ + ops[OP_AIC + base + l + p * NRT] * st[h13S + pj]

398+ + ops[OP_AI + base + l + p * NRT] * st[A_NLF + pj];

399+ }

400+ st[h3nS + i] = sum;

401+}

402+`;

403+

404+ const J = `${header}

405+@compute @workgroup_size(${wg})

406+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {

407+${gate}

408+ let h3nS = ST_UH + ((k + 1u) % 3u) * N;

409+ let uS = ST_U + (k % 3u) * N;

410+ let u1S = ST_U + ((k + 2u) % 3u) * N;

411+ let unS = ST_U + ((k + 1u) % 3u) * N;

412+ var sum: f32 = 0.0;

413+ for (var q: u32 = 0u; q < NTH; q = q + 1u) {

414+ sum = sum + st[h3nS + l + q * NRT] * ops[OP_GI + q + j * NTH];

415+ }

416+ st[unS + i] = sum;

417+ if (l == LP0 && j == JP0) {

418+ outb[k] = (sum - 2.0 * st[uS + i] + st[u1S + i]) * HT2I;

419+ }

420+}

421+`;

422+

423+ const sources = [

424+ AB,

425+ transform('A_RHSA', 'A_RHSAF', 'OP_F'),

426+ D,

427+ transform('A_PH3', 'A_PHI', 'OP_GI'),

428+ FG,

429+ transform('A_NL', 'A_NLF', 'OP_F'),

430+ I,

431+ J,

432+ ];

433+

434+ const layout: Layout = { N, opsData, stLen, SC_CVEL: ST_CVEL };

435+ return { layout, sources, nwg: Math.ceil(N / wg) };

436+}

437+

438+export async function requestGpuDevice(): Promise<GPUDevice> {

439+ const gpu = (navigator as { gpu?: GPU }).gpu;

440+ if (!gpu) throw new Error('WebGPU is not available in this browser');

441+ const adapter = await gpu.requestAdapter();

442+ if (!adapter) throw new Error('no WebGPU adapter');

443+ return adapter.requestDevice();

444+}

445+

446+export async function runGpuSim(

447+ device: GPUDevice,

448+ su: GpuSetup,

449+ onProgress?: (stepsDone: number, nsteps: number) => void,

450+): Promise<GpuSimResult> {

451+ const wg = 64;

452+ const chunk = 1024;

453+ const nsteps = su.nsteps;

454+ const { layout, sources, nwg } = buildLayoutAndStages(su, wg);

455+

456+ const bgl = device.createBindGroupLayout({

457+ entries: [

458+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },

459+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },

460+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },

461+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform', hasDynamicOffset: true, minBindingSize: 16 } },

462+ ],

463+ });

464+ const plLayout = device.createPipelineLayout({ bindGroupLayouts: [bgl] });

465+

466+ const pipelines: GPUComputePipeline[] = [];

467+ for (const code of sources) {

468+ const module = device.createShaderModule({ code });

469+ const info = await module.getCompilationInfo();

470+ const errs = info.messages.filter((m) => m.type === 'error');

471+ if (errs.length) {

472+ throw new Error('WGSL: ' + errs.map((m) => `${m.lineNum}:${m.linePos} ${m.message}`).join('; '));

473+ }

474+ pipelines.push(device.createComputePipeline({ layout: plLayout, compute: { module, entryPoint: 'main' } }));

475+ }

476+

477+ const mkBuf = (size: number, usage: number) => device.createBuffer({ size, usage });

478+ const opsBuf = mkBuf(layout.opsData.byteLength, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);

479+ device.queue.writeBuffer(opsBuf, 0, layout.opsData);

480+ const stBuf = mkBuf(layout.stLen * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST);

481+ const outBuf = mkBuf(2 * nsteps * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC);

482+ const uniBuf = mkBuf(chunk * 256, GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST);

483+ {

484+ const uni = new Uint32Array(chunk * 64);

485+ for (let s = 0; s < chunk; s++) uni[s * 64] = s;

486+ device.queue.writeBuffer(uniBuf, 0, uni);

487+ }

488+ const readOut = mkBuf(2 * nsteps * 4, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);

489+ const readSt = mkBuf(layout.stLen * 4, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);

490+

491+ const bindGroup = device.createBindGroup({

492+ layout: bgl,

493+ entries: [

494+ { binding: 0, resource: { buffer: opsBuf } },

495+ { binding: 1, resource: { buffer: stBuf } },

496+ { binding: 2, resource: { buffer: outBuf } },

497+ { binding: 3, resource: { buffer: uniBuf, offset: 0, size: 16 } },

498+ ],

499+ });

500+

501+ const t0 = performance.now();

502+ const pending: Promise<undefined>[] = [];

503+ let kbase = 0;

504+ while (kbase < nsteps) {

505+ const n = Math.min(chunk, nsteps - kbase);

506+ const kb = new Float32Array(new ArrayBuffer(4));

507+ kb[0] = kbase;

508+ device.queue.writeBuffer(stBuf, 0, kb);

509+ const enc = device.createCommandEncoder();

510+ const pass = enc.beginComputePass();

511+ for (let s = 0; s < n; s++) {

512+ pass.setBindGroup(0, bindGroup, [s * 256]);

513+ for (const p of pipelines) {

514+ pass.setPipeline(p);

515+ pass.dispatchWorkgroups(nwg);

516+ }

517+ }

518+ pass.end();

519+ device.queue.submit([enc.finish()]);

520+ pending.push(device.queue.onSubmittedWorkDone());

521+ if (pending.length >= 2) await pending.shift();

522+ kbase += n;

523+ if (onProgress) onProgress(kbase, nsteps);

524+ }

525+ while (pending.length) await pending.shift();

526+ const gpuMs = performance.now() - t0;

527+

528+ const enc = device.createCommandEncoder();

529+ enc.copyBufferToBuffer(outBuf, 0, readOut, 0, 2 * nsteps * 4);

530+ enc.copyBufferToBuffer(stBuf, 0, readSt, 0, layout.stLen * 4);

531+ device.queue.submit([enc.finish()]);

532+ await readOut.mapAsync(GPUMapMode.READ);

533+ await readSt.mapAsync(GPUMapMode.READ);

534+ const outData = new Float32Array(readOut.getMappedRange().slice(0));

535+ const state = new Float32Array(readSt.getMappedRange().slice(0));

536+ readOut.unmap();

537+ readSt.unmap();

538+ for (const b of [opsBuf, stBuf, outBuf, uniBuf, readOut, readSt]) b.destroy();

539+

540+ const acc = outData.subarray(0, nsteps);

541+ for (let i = 0; i < acc.length; i++) {

542+ if (!Number.isFinite(acc[i])) {

543+ throw new Error(

544+ `cymbal2: simulation diverged at t = ${(i * su.ht).toFixed(4)} s ` +

545+ '(reduce strike_velocity, or raise n_theta / oversample)',

546+ );

547+ }

548+ }

549+

550+ return {

551+ acc,

552+ fhist: outData.subarray(nsteps, 2 * nsteps),

553+ peakXi: state[12],

554+ peakDisp: state[13],

555+ contactVel: Array.from(state.subarray(layout.SC_CVEL, layout.SC_CVEL + su.K)),

556+ strokesPlayed: Math.min(state[4] + (state[3] > 0 ? 1 : 0), su.K),

557+ gpuMs,

558+ };

559+}

webapp/src/worker/renderWorker.tsmodified+319−78View file

@@ -1,18 +1,25 @@

1-// Web Worker that renders the cymbal2 stroke by executing the MATLAB

2-// implementation in ../matlab through numbl. Unlike the sibling cymbal app

3-// there is no cheap stage to cache around: the cost IS the time integration,

4-// one long executeCode call, so requests carry a quality tier

1+// Web Worker that renders the cymbal2 stroke. The MATLAB implementation in

2+// ../matlab remains the source of truth; what varies is who runs the time

3+// loop. When the browser exposes WebGPU, the worker executes the unmodified

4+// setup stage (grid, per-wavenumber inverses, mallet, stroke schedule)

5+// through numbl, hands the operators to a WGSL port of cymbal2_sim.m

6+// (./gpuSim.ts), and feeds the returned acceleration back through the

7+// unmodified cymbal2_output.m; this is 20x or more faster than interpreting

8+// the loop. Without WebGPU, or when cymbal2_sim.m has been edited in the app

9+// (the WGSL mirrors the shipped step, so edits there must fall back), the

10+// original interpreter path runs instead. Requests carry a quality tier

511 // preview -> oversample 1, tail clamped to 1.5 s (the azimuthal grid is

612 // kept: hard strokes need the full wavenumber band to carry

713 // the cascade, and coarsening it destabilizes the integration)

814 // full -> the parameters as given

9-// and the UI shows elapsed time rather than a meaningful progress fraction.

15+// and the UI shows elapsed time, plus integration progress on the GPU path.

1016 // A cheap linear-modes stage (cymbal2_modes per wavenumber) is cached on the

1117 // geometry so the spectrum plot can overlay the mode table immediately.

1218

13-import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';

14-import { hashMFiles, type MFile } from '../mfiles';

19+import { executeCode, RTV, type RuntimeValue, type RuntimeTensor } from 'numbl';

20+import { hashMFiles, DEFAULT_MFILES, type MFile } from '../mfiles';

1521 import { MODES_KEYS, type Cymbal2Params } from '../params';

22+import { requestGpuDevice, runGpuSim, type GpuSetup } from './gpuSim';

1623

1724 export interface RenderRequest {

1825 type: 'render';

@@ -53,6 +60,7 @@ export interface RenderResponse {

5360 renderMs: number;

5461 oversample: number;

5562 nTheta: number;

63+ backend: 'webgpu' | 'cpu';

5664 warnings: string[];

5765 };

5866 }

@@ -138,6 +146,24 @@ function toF32(v: RuntimeValue): Float32Array {

138146 return out;

139147 }

140148

149+// Max-pool the contact force history to about 2000 points for the plot.

150+function poolForce(fhist: Float32Array, simRate: number): { tpool: Float32Array; fpool: Float32Array } {

151+ const nsim = fhist.length;

152+ const np0 = Math.min(2000, nsim);

153+ const blk = Math.ceil(nsim / np0);

154+ const np = Math.ceil(nsim / blk);

155+ const fpool = new Float32Array(np);

156+ const tpool = new Float32Array(np);

157+ for (let b = 0; b < np; b++) {

158+ let m = 0;

159+ const end = Math.min(nsim, (b + 1) * blk);

160+ for (let i = b * blk; i < end; i++) m = Math.max(m, fhist[i]);

161+ fpool[b] = m;

162+ tpool[b] = (b * blk + blk / 2) / simRate;

163+ }

164+ return { tpool, fpool };

165+}

166+

141167 const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;

142168

143169 self.onmessage = (ev: MessageEvent<RenderRequest>) => {

@@ -150,31 +176,247 @@ self.onmessage = (ev: MessageEvent<RenderRequest>) => {

150176 warnings.push(args.map(String).join(' '));

151177 };

152178

179+ render(msg)

180+ .catch((err) => {

181+ const e = err as { message?: string; file?: string; line?: number };

182+ let m: ErrorMessage = {

183+ type: 'error',

184+ id: msg.id,

185+ message: e?.message ? String(e.message) : String(err),

186+ file: typeof e?.file === 'string' ? e.file : undefined,

187+ line: typeof e?.line === 'number' ? e.line : undefined,

188+ };

189+ for (const w of warnings) {

190+ const hit = SKIP_WARNING.exec(w);

191+ if (hit) {

192+ m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };

193+ break;

194+ }

195+ }

196+ (self as unknown as Worker).postMessage(m);

197+ })

198+ .finally(() => {

199+ console.warn = origWarn;

200+ });

201+};

202+

203+// Everything the response assembly needs, produced by either backend.

204+interface SimOutcome {

205+ y32: Float32Array;

206+ peakAccel: number;

207+ nDropped: number;

208+ tpool: Float32Array;

209+ fpool: Float32Array;

210+ backend: 'webgpu' | 'cpu';

211+ d: {

212+ strokesPlayed: number;

213+ meanV: number;

214+ contactVel: number;

215+ peakForceN: number;

216+ contactMs: number;

217+ contactEvents: number;

218+ peakXiMm: number;

219+ peakDispM: number;

220+ peakDispThickness: number;

221+ nsteps: number;

222+ simRate: number;

223+ strikeRadiusM: number;

224+ pickupRadiusM: number;

225+ RsM: number;

226+ };

227+}

228+

229+// The stroke-schedule block from cymbal2_sim.m, verbatim, so the GPU render

230+// integrates exactly the schedule the interpreter would.

231+const SCHEDULE_SRC = `

232+K = max(1, round(P.strike_count));

233+tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);

234+vstr = zeros(1, K);

235+state = P.seed * 2654435761;

236+for i = 1:K

237+ state = mod(state * 1664525 + 1013904223, 4294967296);

238+ x = state / 4294967296;

239+ vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));

240+end

241+gap = P.strike_gap_m;

242+if gap <= 0

243+ gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;

244+end

245+travel = gap ./ vstr + vstr / (2 * P.strike_accel);

246+tlaunch = tstr - travel;

247+off = -min(0, min(tlaunch));

248+tlaunch = tlaunch + off;

249+duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;

250+nsteps = floor(duration * P.sample_rate * P.oversample);

251+`;

252+

253+async function renderOnGpu(

254+ mfiles: MFile[],

255+ params: Cymbal2Params,

256+ onProgress: (frac: number) => void,

257+): Promise<SimOutcome> {

258+ // A WebGPU device first: if there is none, fail before spending seconds on

259+ // the setup stage.

260+ const device = await requestGpuDevice();

153261 try {

154- render(msg);

155- } catch (err) {

156- const e = err as { message?: string; file?: string; line?: number };

157- let m: ErrorMessage = {

158- type: 'error',

159- id: msg.id,

160- message: e?.message ? String(e.message) : String(err),

161- file: typeof e?.file === 'string' ? e.file : undefined,

162- line: typeof e?.line === 'number' ? e.line : undefined,

262+ const r = runStage(mfiles, `

263+${pLiteral(params)}

264+S = cymbal2_system(P);

265+MB = cymbal2_mallet(P);

266+G = S.G;

267+o_D1 = G.D1; o_D2 = G.D2;

268+o_F = G.F; o_Gi = G.Gi;

269+o_Rd = G.Rd; o_R2d = G.R2d;

270+o_AiB = S.AiB; o_AiC = S.AiC; o_Ai = S.Ai; o_Fi = S.Fi;

271+o_ASiB = MB.ASiB; o_ASiC = MB.ASiC;

272+o_vM = MB.vM; o_vR = MB.vR; o_vK = MB.vK; o_tipcol = MB.tipcol;

273+${SCHEDULE_SRC}

274+o_tlaunch = tlaunch; o_vstr = vstr;

275+c_nrt = G.nrt; c_Nr = G.Nr; c_nth = G.nth; c_hr = G.hr; c_hth = G.hth;

276+c_ht = S.ht; c_Eh = S.Eh; c_Rs = S.Rs; c_rhoh = S.rhoh;

277+c_fsc = S.forceScale; c_nlr0 = S.eomRow0; c_nDropped = S.nDropped;

278+c_ls = S.ls; c_js = S.js; c_lp = S.lp; c_jp = S.jp;

279+c_strikeR = S.strikeRadiusM; c_pickupR = S.pickupRadiusM;

280+c_tipi = MB.tipi; c_nfree = MB.nfree; c_ps = MB.ps;

281+c_holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / S.ht));

282+c_missDepth = gap + 0.06;

283+c_gap = gap; c_K = K; c_nsteps = nsteps;

284+c_simRate = P.sample_rate * P.oversample;

285+`);

286+ const v = r.variableValues;

287+ const s = (n: string) => asScalar(v[n]);

288+ const t = (n: string) => toF32(v[n]);

289+ const setup: GpuSetup = {

290+ ops: {

291+ D1: t('o_D1'), D2: t('o_D2'),

292+ F: t('o_F'), Gi: t('o_Gi'),

293+ Rd: t('o_Rd'), R2d: t('o_R2d'),

294+ AiB: t('o_AiB'), AiC: t('o_AiC'), Ai: t('o_Ai'), Fi: t('o_Fi'),

295+ ASiB: t('o_ASiB'), ASiC: t('o_ASiC'),

296+ vM: t('o_vM'), vR: t('o_vR'), vK: t('o_vK'), tipcol: t('o_tipcol'),

297+ },

298+ tlaunch: Array.from((v.o_tlaunch as RuntimeTensor).data),

299+ vstr: Array.from((v.o_vstr as RuntimeTensor).data),

300+ nrt: s('c_nrt'), Nr: s('c_Nr'), nth: s('c_nth'),

301+ hr: s('c_hr'), hth: s('c_hth'),

302+ ht: s('c_ht'), Eh: s('c_Eh'), Rs: s('c_Rs'), rhoh: s('c_rhoh'),

303+ fsc: s('c_fsc'), nlr0: s('c_nlr0'),

304+ ls: s('c_ls'), js: s('c_js'), lp: s('c_lp'), jp: s('c_jp'),

305+ tipi: s('c_tipi'), nfree: s('c_nfree'), ps: s('c_ps'),

306+ KM: params.head_stiffness, alM: params.head_alpha,

307+ aH: params.strike_accel, aHr: params.retract_accel,

308+ holdSteps: s('c_holdSteps'), missDepth: s('c_missDepth'),

309+ gap: s('c_gap'), K: s('c_K'), nsteps: s('c_nsteps'),

163310 };

164- for (const w of warnings) {

165- const hit = SKIP_WARNING.exec(w);

166- if (hit) {

167- m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };

168- break;

169- }

311+ const simRate = s('c_simRate');

312+

313+ const sim = await runGpuSim(device, setup, (done, n) => onProgress(done / n));

314+

315+ // The output stage through the unmodified cymbal2_output.m.

316+ const accT = RTV.tensor(Float64Array.from(sim.acc), [1, sim.acc.length]);

317+ const ro = runStage(mfiles, `

318+${pLiteral(params)}

319+[y, peakAccel] = cymbal2_output(P, acc);

320+`, { acc: accT });

321+

322+ // Contact bookkeeping from the force history, as cymbal2_sim does it.

323+ let peakF = 0, contactSteps = 0, events = 0, prev = 0;

324+ for (let i = 0; i < sim.fhist.length; i++) {

325+ const f = sim.fhist[i];

326+ if (f > peakF) peakF = f;

327+ if (f > 0) { contactSteps++; if (prev === 0) events++; }

328+ prev = f;

170329 }

171- (self as unknown as Worker).postMessage(m);

330+ const { tpool, fpool } = poolForce(sim.fhist, simRate);

331+ const meanV = setup.vstr.reduce((a, b) => a + b, 0) / setup.vstr.length;

332+

333+ return {

334+ y32: toF32(ro.variableValues.y),

335+ peakAccel: asScalar(ro.variableValues.peakAccel),

336+ nDropped: s('c_nDropped'),

337+ tpool,

338+ fpool,

339+ backend: 'webgpu',

340+ d: {

341+ strokesPlayed: sim.strokesPlayed,

342+ meanV,

343+ contactVel: sim.contactVel[0] ?? 0,

344+ peakForceN: peakF,

345+ contactMs: (1000 * contactSteps) / simRate,

346+ contactEvents: events,

347+ peakXiMm: sim.peakXi * 1000,

348+ peakDispM: sim.peakDisp,

349+ peakDispThickness: sim.peakDisp / params.thickness_m,

350+ nsteps: setup.nsteps,

351+ simRate,

352+ strikeRadiusM: s('c_strikeR'),

353+ pickupRadiusM: s('c_pickupR'),

354+ RsM: setup.Rs,

355+ },

356+ };

172357 } finally {

173- console.warn = origWarn;

358+ device.destroy();

174359 }

175-};

360+}

176361

177-function render(msg: RenderRequest): void {

362+function renderOnCpu(mfiles: MFile[], params: Cymbal2Params): SimOutcome {

363+ const r = runStage(mfiles, `

364+${pLiteral(params)}

365+S = cymbal2_system(P);

366+MB = cymbal2_mallet(P);

367+[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);

368+[y, peakAccel] = cymbal2_output(P, acc);

369+% max-pool the contact force history to about 2000 points for the plot

370+nsim = length(fhist);

371+np = min(2000, nsim);

372+blk = ceil(nsim / np);

373+np = ceil(nsim / blk);

374+fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);

375+tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;

376+d_strokes = dinfo.strokesPlayed;

377+d_meanV = dinfo.meanV;

378+d_cvel = dinfo.contactVel(1);

379+d_peakF = dinfo.peakForceN;

380+d_contactMs = dinfo.contactMs;

381+d_events = dinfo.contactEvents;

382+d_xi = dinfo.peakXiMm;

383+d_disp = dinfo.peakDispM;

384+d_dispH = dinfo.peakDispThickness;

385+d_nsteps = dinfo.nsteps;

386+d_simRate = dinfo.simRate;

387+d_dropped = S.nDropped;

388+d_strikeR = S.strikeRadiusM;

389+d_pickupR = S.pickupRadiusM;

390+d_Rs = S.Rs;

391+`);

392+ const v = r.variableValues;

393+ return {

394+ y32: toF32(v.y),

395+ peakAccel: asScalar(v.peakAccel),

396+ nDropped: asScalar(v.d_dropped),

397+ tpool: toF32(v.tpool),

398+ fpool: toF32(v.fpool),

399+ backend: 'cpu',

400+ d: {

401+ strokesPlayed: asScalar(v.d_strokes),

402+ meanV: asScalar(v.d_meanV),

403+ contactVel: asScalar(v.d_cvel),

404+ peakForceN: asScalar(v.d_peakF),

405+ contactMs: asScalar(v.d_contactMs),

406+ contactEvents: asScalar(v.d_events),

407+ peakXiMm: asScalar(v.d_xi),

408+ peakDispM: asScalar(v.d_disp),

409+ peakDispThickness: asScalar(v.d_dispH),

410+ nsteps: asScalar(v.d_nsteps),

411+ simRate: asScalar(v.d_simRate),

412+ strikeRadiusM: asScalar(v.d_strikeR),

413+ pickupRadiusM: asScalar(v.d_pickupR),

414+ RsM: asScalar(v.d_Rs),

415+ },

416+ };

417+}

418+

419+async function render(msg: RenderRequest): Promise<void> {

178420 const { id, mfiles, quality } = msg;

179421 const params = { ...msg.params };

180422 if (quality === 'preview') {

@@ -209,46 +451,44 @@ end

209451 lruSet(modesCache, modesKey, modesEntry, MODES_CACHE_MAX);

210452 }

211453

212- // ---- the time integration, one long call ----

213- post({ type: 'progress', id, stage: 'integrate', frac: 0.2 });

454+ // ---- the time integration ----

455+ post({ type: 'progress', id, stage: 'integrate', frac: 0.1 });

214456 const t0 = performance.now();

215- const r = runStage(mfiles, `

216-${pLiteral(params)}

217-S = cymbal2_system(P);

218-MB = cymbal2_mallet(P);

219-[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);

220-[y, peakAccel] = cymbal2_output(P, acc);

221-% max-pool the contact force history to about 2000 points for the plot

222-nsim = length(fhist);

223-np = min(2000, nsim);

224-blk = ceil(nsim / np);

225-np = ceil(nsim / blk);

226-fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);

227-tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;

228-d_strokes = dinfo.strokesPlayed;

229-d_meanV = dinfo.meanV;

230-d_cvel = dinfo.contactVel(1);

231-d_peakF = dinfo.peakForceN;

232-d_contactMs = dinfo.contactMs;

233-d_events = dinfo.contactEvents;

234-d_xi = dinfo.peakXiMm;

235-d_disp = dinfo.peakDispM;

236-d_dispH = dinfo.peakDispThickness;

237-d_nsteps = dinfo.nsteps;

238-d_simRate = dinfo.simRate;

239-d_dropped = S.nDropped;

240-d_strikeR = S.strikeRadiusM;

241-d_pickupR = S.pickupRadiusM;

242-d_Rs = S.Rs;

243-`);

457+ const warnings: string[] = [];

458+

459+ // The WGSL mirrors the shipped cymbal2_sim.m, so an edited step falls back

460+ // to the interpreter, where the edit is live. Edits to the setup and output

461+ // stages flow through numbl on both paths.

462+ const simSource = mfiles.find((f) => f.name === 'cymbal2_sim.m')?.source;

463+ const simShipped = DEFAULT_MFILES.find((f) => f.name === 'cymbal2_sim.m')?.source;

464+ const simEdited = simSource !== simShipped;

465+

466+ let out: SimOutcome | null = null;

467+ if (!simEdited) {

468+ try {

469+ out = await renderOnGpu(mfiles, params, (frac) =>

470+ post({ type: 'progress', id, stage: 'integrate', frac: 0.1 + 0.85 * frac }),

471+ );

472+ } catch (err) {

473+ const detail = err instanceof Error ? err.message : String(err);

474+ if (/diverged/.test(detail)) throw err;

475+ // console.warn is intercepted during a render; use info so the reason

476+ // is visible in the devtools console.

477+ console.info(`cymbal2: WebGPU unavailable (${detail}); rendering on the CPU`);

478+ }

479+ }

480+ if (!out) {

481+ out = renderOnCpu(mfiles, params);

482+ if (simEdited) {

483+ // Not a fault: just say why the render took the slow path.

484+ warnings.push('cymbal2_sim.m differs from the shipped source, so the render used the interpreter (the GPU path mirrors the shipped step)');

485+ }

486+ }

244487 const renderMs = performance.now() - t0;

245488

246- post({ type: 'progress', id, stage: 'output', frac: 0.95 });

247- const y32 = toF32(r.variableValues.y);

248- const peakAccel = asScalar(r.variableValues.peakAccel);

249- const nDropped = asScalar(r.variableValues.d_dropped);

489+ post({ type: 'progress', id, stage: 'output', frac: 0.97 });

490+ const { y32, peakAccel, nDropped } = out;

250491

251- const warnings: string[] = [];

252492 if (nDropped > 0) {

253493 warnings.push(

254494 `${nDropped} azimuthal columns were truncated: the radial grid cannot resolve their inner boundary layer; raise n_radial`,

@@ -270,31 +510,32 @@ d_Rs = S.Rs;

270510 quality,

271511 y: y32,

272512 sampleRate: params.sample_rate,

273- forceT: toF32(r.variableValues.tpool),

274- forceN: toF32(r.variableValues.fpool),

513+ forceT: out.tpool,

514+ forceN: out.fpool,

275515 modeF: modesEntry.modeF.slice(),

276516 modeM: modesEntry.modeM.slice(),

277517 diag: {

278- strokesPlayed: asScalar(r.variableValues.d_strokes),

279- meanV: asScalar(r.variableValues.d_meanV),

280- contactVel: asScalar(r.variableValues.d_cvel),

281- peakForceN: asScalar(r.variableValues.d_peakF),

282- contactMs: asScalar(r.variableValues.d_contactMs),

283- contactEvents: asScalar(r.variableValues.d_events),

284- peakXiMm: asScalar(r.variableValues.d_xi),

285- peakDispMm: asScalar(r.variableValues.d_disp) * 1000,

286- peakDispThickness: asScalar(r.variableValues.d_dispH),

518+ strokesPlayed: out.d.strokesPlayed,

519+ meanV: out.d.meanV,

520+ contactVel: out.d.contactVel,

521+ peakForceN: out.d.peakForceN,

522+ contactMs: out.d.contactMs,

523+ contactEvents: out.d.contactEvents,

524+ peakXiMm: out.d.peakXiMm,

525+ peakDispMm: out.d.peakDispM * 1000,

526+ peakDispThickness: out.d.peakDispThickness,

287527 peakAccel,

288- strikeRadiusM: asScalar(r.variableValues.d_strikeR),

289- pickupRadiusM: asScalar(r.variableValues.d_pickupR),

290- RsM: asScalar(r.variableValues.d_Rs),

291- nsteps: asScalar(r.variableValues.d_nsteps),

292- simRate: asScalar(r.variableValues.d_simRate),

528+ strikeRadiusM: out.d.strikeRadiusM,

529+ pickupRadiusM: out.d.pickupRadiusM,

530+ RsM: out.d.RsM,

531+ nsteps: out.d.nsteps,

532+ simRate: out.d.simRate,

293533 nDropped,

294534 durationS: y32.length / params.sample_rate,

295535 renderMs,

296536 oversample: params.oversample,

297537 nTheta: params.n_theta,

538+ backend: out.backend,

298539 warnings,

299540 },

300541 };