docs: measured GPU timings (RTX 6000 Ada) for the eigensolver and the full model
2 changed files+15−5
README.mdmodified+13−3View file
@@ -1,6 +1,6 @@
11 # cymsim
22
3-A Python package for simulating a drumstick strike on a thin metal plate or a cymbal. A plate is described in Python (material with its loss constants, bounding box, outline, thickness field, an optional stress-free profile) and sampled on a regular grid. The sound is computed in one of two ways. The **linear model** computes the modes of the plate once per design and synthesizes any strike from them in about a second, on the CPU. The **full model** solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly in time with PyTorch, on a GPU when one is available; it takes about a minute per second of sound on a GPU (several minutes on a CPU) and captures the pitch glide and the wash of a hard hit that the linear model cannot. Both use the same grid, so the linear sound is the small-amplitude limit of the full one, and a light hit sounds the same either way.
3+A Python package for simulating a drumstick strike on a thin metal plate or a cymbal. A plate is described in Python (material with its loss constants, bounding box, outline, thickness field, an optional stress-free profile) and sampled on a regular grid. The sound is computed in one of two ways. The **linear model** computes the modes of the plate once per design and synthesizes any strike from them in about a second, on the CPU. The **full model** solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly in time with PyTorch, on a GPU when one is available; it takes a few seconds per second of sound on a GPU (about 90 s on a CPU) and captures the pitch glide and the wash of a hard hit that the linear model cannot. Both use the same grid, so the linear sound is the small-amplitude limit of the full one, and a light hit sounds the same either way.
44
55 The package is a port of the browser application [cymbalsim](../cymbalsim) to a library that can be scripted, so that designs, strikes and the numerics can be explored and compared programmatically. The model and the numerical method are described in [THEORY.md](THEORY.md).
66
@@ -75,9 +75,19 @@ The profile is treated as a shallow shell, valid while its slopes are small comp
7575
7676 ## The two models
7777
78-`compute_modes(plate, fmax)` assembles the linearized stiffness as a sparse matrix and computes all modes below `fmax` with a dense LAPACK solve restricted to that band (SciPy), or on a CUDA GPU with `torch.linalg.eigh` when one is available (`solver="torch"`). A flat plate has one unknown per node and its 16 inch crash at 8 mm takes about 4 s; a plate with a profile couples the transverse to the in-plane motion, has three unknowns per node, and takes about 40 s on the CPU. The result, a `Modes` object, holds the frequencies, the mass-normalized shapes on the grid, and the three loss rates of each mode; `modes.save` and `Modes.load` store it. `synthesize(modes, strike, listener, duration)` then produces a `Sound` in about a second; the strike and the listener can be changed at no further cost.
78+`compute_modes(plate, fmax)` assembles the linearized stiffness as a sparse matrix and computes all modes below `fmax` with a dense LAPACK solve restricted to that band (SciPy), or on a CUDA GPU with `torch.linalg.eigh` when one is available (`solver="torch"`). A flat plate has one unknown per node and its 16 inch crash at 8 mm takes about 3 s on a 12-thread CPU (0.3 s on an RTX 6000 Ada); a plate with a profile couples the transverse to the in-plane motion, has three unknowns per node, and takes about 25 to 40 s on the CPU (1.6 s on the GPU). The result, a `Modes` object, holds the frequencies, the mass-normalized shapes on the grid, and the three loss rates of each mode; `modes.save` and `Modes.load` store it. `synthesize(modes, strike, listener, duration)` then produces a `Sound` in about a second; the strike and the listener can be changed at no further cost.
7979
80-`simulate(plate, strike, listener, duration, device=None, dtype="float32")` runs the full model. The time step is set by the in-plane wave speed (about 1.9 microseconds for bronze at 8 mm, so 4 s of sound is about 2 million steps). The step is written as a pure function of PyTorch tensors and compiled with `torch.compile`, which fuses its many small operations; on a CPU this gives about 90 s per second of sound for the 16 inch crash, and a GPU is expected to be several times faster (the code path is the same, but we have not timed it on a GPU). `Simulation` exposes the run batch by batch, with access to the state, for experiments. Single precision is the default, with the rounding safeguards described in THEORY.md; `dtype="float64"` is available.
80+`simulate(plate, strike, listener, duration, device=None, dtype="float32")` runs the full model. The time step is set by the in-plane wave speed (about 1.9 microseconds for bronze at 8 mm, so 4 s of sound is about 2 million steps). The step is written as a pure function of PyTorch tensors and compiled with `torch.compile`, which fuses its many small operations; on CUDA the steps of one audio sample are further captured into a CUDA graph and replayed, so that a sample costs one launch. Measured times per second of sound, in single precision:
81+
82+| Case | Nodes | 12-thread CPU | RTX 6000 Ada |
83+| --- | --- | --- | --- |
84+| 16 inch crash, 8 mm grid | 2 156 | 90 s | 3.9 s |
85+| 16 inch crash with bell, 8 mm | 2 156 | | 3.9 s |
86+| 16 inch crash, 4 mm grid | 8 389 | | 8.9 s |
87+| Thunder sheet, 8 mm | 8 855 | | 7.2 s |
88+| Thunder sheet, 4 mm | 35 190 | | 18.5 s |
89+
90+The GPU is far from saturated at these sizes (about 7 microseconds per step whatever the grid, up to some ten thousand nodes), so a finer grid costs little more than the shorter time step it requires. Double precision (`dtype="float64"`) is about 7 times slower on that GPU. `Simulation` exposes the run batch by batch, with access to the state, for experiments. Single precision is the default, with the rounding safeguards described in THEORY.md.
8191
8292 Both return a `Sound` with the listener pressure in pascals, the optional movie of the transverse displacement, and what the run reported (peak force, contact time). `Sound.to_wav` writes a 16-bit file at the peak, at a loudness target, or at an absolute level; `Sound.spectrogram` computes an STFT.
8393
THEORY.mdmodified+2−2View file
@@ -65,7 +65,7 @@ Where $D$ is constant the interior stencil is the standard 13-point biharmonic.
6565
6666 **Modes.** Collecting the linear part on the grid gives $M \ddot q + K q = f$, with $q$ the vector of unknowns (the transverse displacement of every plate node for a flat plate; $u$, $v$, $w$ interleaved for a plate with a profile), $M$ the diagonal lumped mass, and $K$ the symmetric stiffness, the Hessian of $U$ at the rest state; it is assembled as a sparse matrix from the stencil operators of Section 2, and the tests check that $Kq$ equals the linearized force of the full model to rounding. The modes are the solutions of $K\phi = \omega^2 M \phi$, normalized so that $\phi^T M \phi = 1$; the rigid modes at $\omega = 0$ are dropped, since a hanging plate's rigid motion is taken up by the thread and does not radiate in the model below.
6767
68-Here cymsim departs from the browser application, which slices the spectrum and solves each slice by shift-invert Lanczos in a pool of workers. With LAPACK at hand we form $A = M^{-1/2} K M^{-1/2}$ as a dense matrix and call the driver that computes the eigenpairs in a given interval, $(2\pi \cdot 1\,\mathrm{Hz})^2 < \lambda \le (2\pi f_{\max})^2$. A 16 inch cymbal at 8 mm has about 2150 nodes and 1900 bending modes below 16 kHz, computed in about 4 s; with a profile the same cymbal has 6500 unknowns and takes about 40 s on a 12-thread CPU, or a few seconds on a GPU with `torch.linalg.eigh`. The dense approach stores $n^2$ numbers and its cost grows as $n^3$, so it is the right tool up to roughly 10 000 unknowns; beyond that, spectrum slicing would be the remedy.
68+Here cymsim departs from the browser application, which slices the spectrum and solves each slice by shift-invert Lanczos in a pool of workers. With LAPACK at hand we form $A = M^{-1/2} K M^{-1/2}$ as a dense matrix and call the driver that computes the eigenpairs in a given interval, $(2\pi \cdot 1\,\mathrm{Hz})^2 < \lambda \le (2\pi f_{\max})^2$. A 16 inch cymbal at 8 mm has about 2150 nodes and 1900 bending modes below 16 kHz, computed in about 3 s on a 12-thread CPU or 0.3 s on an RTX 6000 Ada with `torch.linalg.eigh`; with a profile the same cymbal has 6500 unknowns and takes 25 to 40 s on the CPU, or 1.6 s on the GPU. The dense approach stores $n^2$ numbers and its cost grows as $n^3$, so it is the right tool up to roughly 10 000 unknowns; beyond that, spectrum slicing would be the remedy.
6969
7070 **Modal equations.** Expanding $q = \sum_m \phi_m q_m(t)$ turns the plate into independent oscillators $\ddot q_m + 2\sigma_m \dot q_m + \omega_m^2 q_m = \phi_m^T f(t)$, with $\sigma_m$ the sum of the three loss rates of the mode. The internal friction gives $\eta\omega_m/2$; the thermoelastic loss gives $\sum_k \eta_{te}(\omega_m, h_k)\,U_{b,k}/\omega_m$, the local loss factor weighted by the bending energy $U_{b,k}$ of the mode at each node, which reduces to $\eta_{te}\omega_m/2$ on a flat uniform plate. The radiation rate is the power the mode radiates over twice its energy. A baffled plate whose nodes move with the velocities $v_k$ radiates the power $v^T R v$ from one face, with the radiation resistance matrix $R_{jk} = (\rho_{\mathrm{air}}\omega^2/4\pi c)\,\Delta A_j \Delta A_k\, \mathrm{sinc}(k_a r_{jk})$ of Elliott and Johnson, $k_a = \omega/c$; both faces radiate, so
7171
@@ -104,7 +104,7 @@ The full model keeps the nonlinear strains of Section 1 and steps the semi-discr
104104
105105 **Single precision.** Three places in the kernel would otherwise turn rounding into audible noise. The first is the membrane strain of a curved design: written directly, $\varepsilon_{xx} = u_x + \tfrac12[(w_0+w)_x^2 - w_{0,x}^2]$ subtracts two nearly equal numbers once $w$ is small compared with the profile, and the rounding error of the two squares is proportional to $w_{0,x}^2$, a constant that does not fall as the vibration decays and acts as a broadband forcing. The strains are therefore written so that every term carries the deviation, $\varepsilon_{xx} = u_x + w_x(w_{0,x} + \tfrac12 w_x)$ with $w$ the deviation alone. The second is the update of $w$: with $w$ of a few millimeters and $\Delta t \approx 2\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, so the update is compensated, with the rounding error of each addition carried in a separate array and added back at the next step. The third is the listener sum, which forms the acceleration as the difference of successive delayed velocities at the audio rate rather than sampling the instantaneous acceleration. What remains is a noise floor proportional to the signal, roughly 75 dB below the peak in the top octave for a strongly curved design and 100 dB for a flat one; `dtype="float64"` removes it at a cost in speed on consumer GPUs.
106106
107-**The kernel.** Here cymsim departs from the browser application's hand-written WGSL kernel. The step is a pure function of PyTorch tensors over the $(n_y, n_x)$ grid: the second differences and the cell quantities are slices and pads of the whole arrays, the membrane forces of the four Gauss points are accumulated per cell corner and scattered to the nodes by four padded additions, the stick is a handful of scalar tensor operations, and the transverse velocity at the sample instants is written into a ring buffer from which the retarded-time listener sum is formed once per batch of one movie frame. As written the step is some two hundred small tensor operations, and on the CPU or a GPU alike it would be bound by their dispatch; `torch.compile` fuses them into a few kernels, which makes the step about six times faster on the CPU (0.17 ms for the 16 inch crash, about 90 s per second of sound) and is what makes a GPU run at all. The same code runs on CUDA, on Apple's MPS and on the CPU.
107+**The kernel.** Here cymsim departs from the browser application's hand-written WGSL kernel. The step is a pure function of PyTorch tensors over the $(n_y, n_x)$ grid: the second differences and the cell quantities are slices and pads of the whole arrays, the membrane forces of the four Gauss points are accumulated per cell corner and scattered to the nodes by four padded additions, the stick is a handful of scalar tensor operations, and the transverse velocity at the sample instants is written into a ring buffer from which the retarded-time listener sum is formed once per batch of one movie frame. As written the step is some two hundred small tensor operations, and on the CPU or a GPU alike it would be bound by their dispatch; `torch.compile` fuses them into a few kernels, which makes the step about six times faster on the CPU (0.17 ms for the 16 inch crash, about 90 s per second of sound). On a GPU the fused step is still bound by the Python call and the kernel launches (about 0.1 ms), so the steps of one audio sample are captured into a CUDA graph on static buffers and replayed, one launch per sample; the step then costs about 7 microseconds on an RTX 6000 Ada, 3.9 s per second of sound for the 16 inch crash at 8 mm, and hardly more at four times the nodes. The same code runs on CUDA, on Apple's MPS and on the CPU.
108108
109109 ## 5. Validation
110110