1# cymbal2-webgpu
3A limitation of the [cymbal2](https://vault1.magland.org/simulated_instruments/cymbal2) instrument is render time: the sound is a nonlinear time-domain integration of a von Karman shell, several hundred thousand strictly sequential steps, and the browser app executes those steps through numbl's interpreter at roughly 1.9 ms per step, which puts a full default render near twenty minutes. Here we explore moving the time loop, and only the time loop, onto the GPU through WebGPU, keeping numbl as the runtime for everything around it. On the modest integrated GPU of a laptop (Intel Iris Xe, running at its 300 MHz idle clock throughout, see below), the same render takes 83 microseconds per step: 48 seconds for the full six-second crash, about 4.4x faster than MATLAB R2026b, 7.5x faster than the numbl CLI with its native addon, and 23x faster than the numbl-in-browser path the webapp actually uses. The audio is single precision throughout; the contact diagnostics match the f64 reference to 0.03%, first contact lands on exactly the same step, and the acoustic rms agrees within 3%, which is within the tolerance the cymbal2 validation suite already grants the nonlinear trajectory.
5## Division of labor
7The design principle is that the MATLAB code in `simulated_instruments/cymbal2/matlab` remains the single source of truth, and the GPU replaces only `cymbal2_sim.m`, the step loop. numbl executes the unmodified setup stage (`cymbal2_system`, `cymbal2_mallet`, and the stroke schedule block copied verbatim from `cymbal2_sim`) and hands the resulting tensors to the host, which packs them into one read-only storage buffer; numbl likewise executes the unmodified output stage (`cymbal2_output`) on the acceleration signal the GPU returns. The WGSL mirrors the step structure of `cymbal2_sim.m` line for line: the same implicit per-wavenumber update, the same boundary rows (already baked into the stored inverses), the same Hertz contact, hand state machine, and beam update, with the sparse radial and azimuthal difference matrices applied as stencils and the dense transforms, per-wavenumber solves, and beam matrices applied as stored matrices. Scalar parameters, node indices, and the launch-step thresholds of the stroke schedule are baked into the generated shader source, and the thresholds are computed on the host in f64 with exactly the comparison `cymbal2_sim` uses, so the gesture is step-for-step aligned with the CPU reference.
9```
10numbl (setup: grid, per-wavenumber inverses, mallet, schedule)
11 -> WebGPU (the time loop: 8 compute dispatches per step)
12 -> numbl (output: DC blocker, decimation, normalize)
13```
15Usage, from this directory (`npm install` first; the `webgpu` package is prebuilt Google Dawn):
17```
18node src/run.mjs # full default render -> cymbal2-gpu.wav
19node src/run.mjs --pinned # the validation suite's reduced configuration
20node src/run.mjs --compare # GPU vs the unmodified cymbal2_sim through numbl
21node src/run.mjs --mode mega # the single-workgroup variant, for comparison
22node src/run.mjs --velocity 4 --count 3 --rate 2 # parameter overrides
23```
25## The step as eight dispatches
27One time step is eight compute dispatches at full GPU width, with the dispatch boundaries providing the global synchronization the stages need:
29| stage | work |
30|---|---|
31| AB | hand state machine, Hertz contact, and beam update (workgroup 0); shell derivative stencils and the Airy right-hand side (all workgroups) |
32| C | azimuthal analysis of the Airy right-hand side |
33| D | Airy solve, one small radial system per wavenumber |
34| E | azimuthal synthesis of the stress function |
35| FG | the von Karman bracket and forcing, stress-function stencils recomputed inline |
36| H | azimuthal analysis of the nonlinear load |
37| I | the implicit shell update, three small radial products per wavenumber |
38| J | azimuthal synthesis of the new displacement; pickup acceleration |
40Because every recorded dispatch must be identical, the step index k is the sum of a chunk base written by the host into the state buffer once per 1024 steps and a within-chunk offset read from a dynamic-offset uniform, and all slot rotation is a pure function of k: the three-generation displacement and wavenumber-state slots rotate every step whether or not the shell is active, which is harmless before first contact because every slot is still zero, and exact afterwards because every active step writes. The shell-active flag latches inside AB; the stages of the same step read it one dispatch later, so they always observe the latch, and the one benign race inside AB itself can only skip work on a field that is still identically zero. The nonlinear trajectory of a cymbal is chaotic, so f32 cannot track the f64 reference for long; what the comparison pins down is the deterministic part (contact step, peak force, contact speed, the first few hundred steps of the acceleration) and the statistics of the rest.
42## What we measured, including what failed
44The exploration tried three architectures, and the measurements on this device (Intel Iris Xe integrated, Mesa/Dawn, one workgroup occupying at most one of six subslices) are worth recording plainly.
46**A single-workgroup megakernel** (`--mode mega`), which loops over thousands of steps inside one dispatch with workgroup barriers between stages, was the first correct implementation and the natural design, since it puts nothing on the CPU and nothing in the dispatch queue. It runs at 396 us/step at workgroup size 256 and 151 us/step at 1024, and profiling showed why: barriers cost almost nothing (about 0.03 us each), but a single workgroup carries too few hardware threads to hide storage-load latency, so the dense 60-point reductions run at about 6 GFLOP/s no matter how they are written. We tried register-accumulator chains, fully unrolled scalar accumulators, vec4 loads with padded layouts, and shared-memory staging; none moved the number, because the bottleneck is latency times too-few-threads, not instruction count. Notably the megakernel keeps the GPU at its full 1300 MHz clock, being continuously busy, and is still slower than the multi-dispatch backend running at 300 MHz.
48**The multi-dispatch backend** (the default) is the opposite trade: full width, but a dispatch boundary between stages. Per-dispatch overhead is about 2 us, and the eight stages sum to 83 us/step. The costs are now spread across the whole GPU, and the same reduction that was latency-bound in one workgroup is unremarkable across twenty-five.
50**A persistent kernel** (`--mode persistent`), thirteen workgroups launched once and synchronized between stages by a global spin barrier on storage atomics, would remove the dispatch boundaries entirely and was the most interesting candidate. It runs without deadlock on this device, and it is wrong: outputs are nondeterministic and diverge shortly after contact. The reason is fundamental rather than a bug in the barrier: WGSL provides no device-scope memory fence, its atomics are relaxed, and `storageBarrier()` is workgroup-scoped, so non-atomic data written by one workgroup is not guaranteed visible to another inside a dispatch, and on this driver it demonstrably is not. Passing every array through atomics would make it correct and slow. The mode is kept as a documented negative result.
52**The frequency governor matters more than any of this.** Throughout the multi-dispatch runs the GPU sits at its 300 MHz minimum clock; the stream of small dispatches separated by hazard barriers never looks busy enough to the i915 governor to ramp, while the megakernel's long dispatches do (1300 MHz), and pinning the clock requires root. Attempts to keep the clock up from inside the application, a concurrent second-device load and a spinner workgroup inside each dispatch, both made things slower: the GPU time-slices between contexts too coarsely, and the spinner extends every dispatch since the barriers wait for it. The practical conclusion is that the 83 us/step figure is a floor set by the governor, not by the work: the arithmetic itself would take roughly a quarter as long at full clock, and a discrete GPU, or a browser whose compositor keeps clocks up, should land well below this number without any code change.
54Timings for the full default render (578k steps at 96 kHz, 18-inch crash, 5.75 m/s stroke), all on the same laptop:
56| path | per step | full render |
57|---|---|---|
58| numbl in Node/browser (the webapp today) | ~1900 us | ~18 min |
59| numbl CLI with native addon | ~620 us | ~6 min |
60| MATLAB R2026b | ~370 us | ~3.6 min |
61| WebGPU multi-dispatch (this, at 300 MHz) | 83 us | 48 s |
63The `src/probe*.mjs` scripts are the microbenchmarks behind these statements. One caution learned the hard way: short-lived timing harnesses on this stack mislead, both from startup transients and from what appears to be early resolution of completion promises under deep queues, and the numbers above come from the instrumented production runner on full-length renders, not from the probes.
65## Correctness
67`--compare` runs the unmodified `cymbal2_sim` through numbl on the reduced pinned configuration (44.1 kHz, no oversampling, 1.5 s tail) and checks the GPU trace against it: first contact on the identical step (1375), peak force 311.23 N against 311.14, tip contact speed -8.224 m/s against -8.223, peak indentation 5.644 mm against 5.643, four contact events in both, and relative L2 of the acceleration over the first 256 steps after contact of about 3e-2, which is f32-versus-f64 on a trajectory with positive Lyapunov exponent, not an indexing error (an indexing error produces order-one differences immediately). The output-stage rms lands at 0.1001 against the suite's pinned 0.0973 with its percent-level tolerance. The megakernel and multi-dispatch backends agree with each other to the same order, differing only in floating-point evaluation order.
69## What this means for the webapp
71Nothing here is Node-specific: the host code speaks standard WebGPU, numbl already runs in the render worker, and the tensor marshalling in `src/setup.mjs` is the same `executeCode` pattern `renderWorker.ts` uses. Integrating this into the cymbal2 webapp would mean adding the kernel generator and the multi-dispatch runner to the worker, using it when `navigator.gpu` yields a device, and keeping the pure-numbl path as the fallback; the preview/full tier distinction would survive unchanged, though at these speeds a preview is nearly interactive (the pinned-size render takes 5.6 s here) and a full render drops from twenty minutes to under a minute. Two caveats carry over from the findings above: the browser path needs the same watchdog-conscious chunking this runner already does (1024 steps per submit, two in flight), and edited MATLAB code in the app's editor would still render through the interpreter, since the WGSL mirrors `cymbal2_sim.m` rather than compiling it; a numbl-to-WGSL compilation of the step, the way turing-surface compiles its models, is the natural next step but a substantially larger one, since this loop is matmul- and branch-shaped rather than elementwise.