cymsim#

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.

The package is a port of the browser application 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.

Installation#

Python 3.10 or later. The linear model needs only NumPy and SciPy; the full model needs PyTorch, and the plotting helpers need Matplotlib.

pip install -e .            # linear model only
pip install -e ".[full]"    # with PyTorch for the full model
pip install -e ".[all]"     # PyTorch, Matplotlib, IPython

On Linux, pip install torch pulls the CUDA build; for a CPU-only machine pip install torch --index-url https://download.pytorch.org/whl/cpu is much smaller.

Quick start#

import cymsim

plate = cymsim.presets.crash_16().evaluate(spacing=0.008)   # a 16 inch tapered bronze crash
print(plate.summary())

modes = cymsim.compute_modes(plate, fmax=16000)             # about 1900 modes, a few seconds
print(modes.table(10))                                      # frequencies and decay times, by mechanism

strike = cymsim.Strike(x=0.9, y=0.5, velocity=3.0)          # 3 m/s near the edge
listener = cymsim.Listener(distance=0.5, angle=30)
sound = cymsim.synthesize(modes, strike, listener, duration=4.0, movie=True)
sound.to_wav("crash_linear.wav")                            # peak-normalized 16-bit WAV

full = cymsim.simulate(plate, strike, listener, duration=4.0)   # the nonlinear model, on the GPU if present
full.to_wav("crash_full.wav")

In a notebook, sound.audio() returns a playable widget, and cymsim.plot draws the plate, the mode shapes, the spectrogram (with the modal frequencies as ticks), the decay times and the movie frames:

from cymsim import plot
plot.plate(plate, strike=strike, listener=listener)
plot.mode(modes, 5)
plot.spectrogram(full, modes=modes)
plot.decay_times(modes)

The examples/ directory has scripts that run these steps.

Designing a plate#

A Design holds the material, the bounding box, and three fields: the outline (inside), the thickness and the stress-free profile. A field is a number or a function of four arrays, (x, y, r, theta), the point relative to the center of the bounding box, its distance from the center and its angle from the x axis, evaluated with NumPy on the whole grid at once. Lengths are meters.

import numpy as np
from cymsim import Design, BRONZE

inch, mm = 0.0254, 1e-3
R = 8 * inch
ramp = lambda a, b, x: np.clip((x - a) / (b - a), 0, 1)

crash = Design(
    name="16 inch crash with a bow",
    material=BRONZE,
    width=2 * R, height=2 * R,
    inside=lambda x, y, r, th: r <= R,
    thickness=lambda x, y, r, th: 1.2 * mm + (0.7 * mm - 1.2 * mm) * ramp(40 * mm, R, r),
    profile=lambda x, y, r, th: 10 * mm * (1 - (r / R) ** 2),
)
plate = crash.evaluate(spacing=0.008)

The Material carries the elastic constants, the loss factor of the metal's internal friction, and (optionally) the thermal constants from which the thermoelastic damping follows. The decay of the sound is not prescribed: it is derived from these constants and from the radiation of each mode into the air, which the linear model computes from the mode's shape. cymsim.presets has ten designs (flat and bowed steel sheets, an aluminum square, a bronze disc, a 16 inch crash, a 20 inch ride, the crash with a bow and a bell, the same hammered, a steel ellipse, a thunder sheet), each a function returning a Design that can be edited before evaluation.

The profile is treated as a shallow shell, valid while its slopes are small compared with one; plate.profile_stats reports the largest slope, and plate.summary() warns above 0.1. The bell of a real cymbal has slopes of about 0.5, so the curved presets use a flatter bell than a real one.

The two models#

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.

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:

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.

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.

Tests#

pytest                  # about 40 s on a 12-thread CPU
pytest -m "not slow"    # skips the two longer cases

The tests are physics checks rather than comparisons with the browser application: the assembled stiffness has exactly the rigid modes in its null space and equals the linearized force of the full model's kernel to 1e-6 on a bowed, tapered plate with a hole; the kernel's forces are the exact gradient of the discrete energy at a strongly nonlinear state; the loss terms only remove energy; the free square and disc frequencies match Gorman's and Leissa's values as the grid is refined; Zener's loss factor peaks where it should; the per-mode radiation rate reproduces the closed form for a piston and Wallace's efficiencies for simply supported shapes, and the FFT convolution equals the direct double sum; a single mode under the stick follows the closed-form damped sinusoid; the attack and the tail of the synthesized signal join; the compiled step equals the uncompiled one; and a light hit simulated both ways gives spectra that agree.

Limitations#

The linear model has no pitch glide and no wash; a hard hit is only a loud one. The full model is the tool for that, and the two agree for a light hit. The von Karman equations hold for moderate rotations, and the Marguerre treatment of the profile for shallow shells, which strains at the bell of a real cymbal. The plate hangs freely: no stand, no thread, no air loading, and the radiation is that of a plate in a rigid baffle, which overstates the low-frequency level of a plate in free air. The full model carries the losses through a two-parameter fit to an average estimate, so its decay agrees with the modal one on average and not mode by mode. Grid dispersion lowers the high partials (by 3 percent near 1 kHz and 10 percent near 3 kHz for a millimeter of bronze at 8 mm), and the staircase edge places the low frequencies within a few percent. The residual stresses of hammering are absent. The von Karman strains are not invariant under a finite rotation and the equations conserve no angular momentum about the in-plane axes, so after a hard strike (a few joules) the free plate slowly tilts and spins in the model, and a tilt of a few degrees then acts as a shell with that slope and pumps energy into the in-plane grid modes, which comes back as an unphysical swell of the sound a second after the strike; the full model therefore projects the six rigid-body modes of the plate out of the state after every batch (remove_rigid, on by default). The dense eigensolver stores the full matrix, so a three-field system much beyond 10 000 unknowns (a curved 16 inch cymbal at 6 mm) becomes slow on the CPU; spectrum slicing, as the browser application does, would lift that limit.

Layout#

• cymsim/material.py, grid.py, plate.py, design.py, params.py, presets.py: materials, the grid and its mask, the plate on the grid, designs and their evaluation, the strike and the listener, the preset designs.
• cymsim/loss.py: the losses from physical constants (internal friction, Zener's thermoelastic damping, the average radiation efficiency, the outline perimeter) and the two-parameter fit of the full model.
• cymsim/linear/: the sparse assembly of K and M, the eigensolver and the Modes object, the per-mode radiation and thermoelastic rates, and the modal synthesis.
• cymsim/full/: the PyTorch kernel of one time step, the driver with the retarded-time listener sum, and the NumPy energy that the tests check the kernel against.
• cymsim/sound.py, spectrogram.py, wav.py, plot.py: the result object, the STFT, the WAV writer with the level mapping, and the Matplotlib views.