cymsimweb: the method, and how it differs from the Python package#

The equations, the geometry, the finite elements and the time integration are those of the cymsim Python package and are written out in its docs/model.md. This document covers what the web application does differently: the arrangement of one sub-step as two GPU dispatches, the projection of the rigid-body modes with a one-step lag, the radiation evaluated at the audio rate, and the buffer layout and ordering argument of the WebGPU kernels. Section 6 records how the two implementations were checked against each other.

1. What is the same#

The mesh is produced by the same algorithm (concentric rings, angular stitching, Laplacian smoothing) in float64 and agrees with the Python mesher to round-off. The DKT bending stiffness is the same element, with the derivatives of the rotation shape functions written out instead of obtained by automatic differentiation; the membrane strains use the CST gradients of the nodal displacements and of the nodal heights; the lumped mass, the inflated rotational inertia, the six M-orthonormal rigid modes, the Irons element bound on the stable time step and the number of sub-steps per audio sample are computed as in Python. The contact law, the patches, the release just above the moving surface, the two loss terms and the leapfrog update are the same. The imperfections are drawn from a different generator (xoshiro128** rather than numpy's PCG64), so a given seed produces a different dent pattern than in Python; the statistics are the same, and a dent-free geometry is identical in both.

2. The internal force without automatic differentiation#

The Python package obtains the internal force as the gradient of the summed element energies by reverse-mode differentiation. Here the gradient is written out. With the membrane strains

the stress resultants $N = C \varepsilon$ and the CST operators $G_x, G_y$ (so that $\partial_x \phi = \sum_i G_{x,i} \phi_i$), the force on node $i$ of an element of area $A$ is

plus $-K_e U_e$ on the nine bending degrees of freedom, where the linear reference run drops the $w_{,x}, w_{,y}$ terms inside $a_x, a_y$ and inside the strains. The $\sigma_1$ term is added per element as $-2 \sigma_1 (\rho h A)_e G_i^{\mathsf T} (G_e \dot w_e)$. These per-element forces are written to a buffer and summed per node through a precomputed node-to-element adjacency; no atomics are needed and the sum is deterministic. The CPU reference loop (src/cpu/simulator.ts) evaluates exactly these expressions in float64, and the shell tests check that the assembled force is the negative gradient of the energy.

3. One sub-step as two dispatches#

A sub-step needs three global synchronizations in the Python arrangement: after the element forces, after the velocity update (to compute the six rigid-mode coefficients), and after the projection. On the GPU each synchronization is a dispatch boundary, and with two to five million sub-steps per second of sound the dispatch count matters. We use two dispatches by applying the rigid-mode projection with a one-step lag.

Element pass (one thread per element, plus one thread per hit): reads $q$ and $v$, writes the fifteen forces of the element; accumulates the element's radiated acceleration (section 4); the hit threads evaluate the contact forces of their patch points from $q$ and the exciter's face height, applying the release at the hit's step.

Node pass (one thread per node, workgroups of 256): sums the element forces of the node and the contact forces of the patch points it belongs to; removes from $v$ and from $q$ the rigid content that the previous step's update introduced, $v \leftarrow v - R G^{-1} R^{\mathsf T} M v$ and the same for $q$, using the twelve sums $R^{\mathsf T} M v$ and $R^{\mathsf T} M q$ that the previous node pass formed; performs the leapfrog update $v' = (v + \Delta t M^{-1} f)/(1 + 2 \sigma_0 \Delta t)$, $q' = q + \Delta t v'$; and contributes its share of $R^{\mathsf T} M v'$ and $R^{\mathsf T} M q'$ to a per-workgroup partial sum (a shared-memory tree reduction) for the next step. Here $R$ holds the six rigid modes before orthonormalisation (rawRigidModes in src/model/shell.ts: the two in-plane translations, the rotation about $z$, the normal translation and the two tilts of the curved shell), which the thread evaluates from the node's position and rest height, and $G = R^{\mathsf T} M R$ is their $6 \times 6$ Gram matrix, whose inverse is a constant of the shader. $R G^{-1} R^{\mathsf T} M$ is the same projection as $\sum_k r_k r_k^{\mathsf T} M$ with the M-orthonormal modes $r_k$ that the CPU loop uses, but it reads three floats per node instead of thirty. States exported before this arrangement hold the sums for the $r_k$; they are converted on import with the triangular matrix $R^{\mathsf T} M r$, which continues such a state exactly. The first thread also advances the exciters, $v_s' = v_s + \Delta t F / m$, $z_s' = z_s + \Delta t v_s'$.

Because the projection is applied one step late, the state carries at any time the rigid content of exactly one step's increment: a rigid velocity of order $F \Delta t / M_{shell}$ and a rigid displacement of order $\Delta t$ times that, both removed at the next step. Projecting $q$ as well as $v$ matters: without it the rigid displacement accumulates over a contact (about $10^{-5}$ m over a mallet stroke at coarse resolution, which changed the peak contact force by $7 \times 10^{-4}$); with it the linear reference run agrees with Python to $10^{-5}$ in peak force and rebound velocity and $10^{-4}$ in the energies.

The partial sums of two consecutive steps live in two halves of one buffer selected by the parity of the step, so the node pass reads one half while writing the other. The step counter is kept in two slots: the element pass copies slot A to slot B, the node pass reads slot B and writes A + 1 to slot A. Every kernel thus reads a slot that no thread of the same dispatch writes, and no uniform has to be updated per dispatch; a chunk of thousands of sub-steps is encoded as one compute pass with only pipeline switches and dispatches.

4. Radiation at the audio rate, and the acceleration record#

The Python package scatters every element's acceleration into retarded-time buffers at sub-step resolution and decimates the result to the audio rate. A scatter with per-element delays needs atomics or a per-element history far too long to keep at sub-step resolution (about 300 MB at the default resolution), so the web application reverses the order: the acceleration is first decimated to the audio rate, then delayed and summed.

The quantity that is decimated is the same as in Python: the normal acceleration of every node, $a_i = n_i \cdot (\ddot u_i, \ddot v_i, \ddot w_i)$ with $n_i$ the unit normal of the rest surface at the node (section 7 of cymsim's model.md), taken as the difference of consecutive projected velocities divided by $\Delta t$, which for the leapfrog scheme is the exact central second difference of $q$. The decimation kernel is also the same: a Kaiser-windowed sinc spanning $W = 32$ audio samples ($\beta = 8$, $-6$ dB at $0.46 f_s$, below $-80$ dB from $0.54 f_s$; decimationKernel in src/model/radiation.ts reproduces cymsim.decimation.decimation_kernel to $10^{-16}$). It is applied in the node pass through $W$ accumulator slots per node. The node stores $a_i$ at every sub-step phase $k$ of audio sample $j_0$, and at the last phase, $k = M - 1$, it adds $\sum_k h_w(k)\, a_i(k)$ to slot $w$ for every $w$; slot $w$ belongs to audio sample $j_0 + w - W/2 + 1$ and lives at ring position (sample mod $W$). When the sub-step count crosses a sample boundary the oldest slot is complete: its value is written to a per-node history ring (a little longer than the longest delay in audio samples), to the chunk's acceleration output when the record is on, and the slot is cleared for the newest sample. The arithmetic is the same as adding $a_i h_w(k)$ to the $W$ slots at every sub-step, but the slots are read and written once per audio sample instead of once per sub-step, which on an integrated GPU saved about 10% of the time of a sub-step. It costs $W M$ multiply-adds per node per audio sample and $(W + M) N$ floats of state, and it needs no atomics because each node owns its slots. A chunk ends one sub-step before its last phase is processed, so the stored accelerations of the current sample are part of the exported state; a state exported before they were kept continues exactly with them set to zero, since its slots already hold the earlier phases. The finished node samples are what the .cymsim file stores (section 8).

The gather is one dispatch per finished audio sample $j$, with one thread per element. Each thread first forms its element's mean of the node sample from the node history and writes it to a per-element history ring, then reads the element's history at the fractional delays of the front path and of the path over the rim with a 16-tap Kaiser-windowed-sinc interpolator (fractionalDelayTaps), and the two sums are reduced in shared memory to one partial sum per workgroup, which the host adds up when it reads the chunk. (The first version ran the gather in a single workgroup, which took about 12 microseconds per sub-step on an integrated GPU, averaged over the sub-steps of a sample.) The interpolator replaces the linear interpolation of the first version, which attenuated the top of the band by a few dB at a half-sample fraction. Against the Python package, whose delays are applied at sub-step resolution before the same decimation, the listener pressure now agrees to $3 \times 10^{-4}$ relative rms on the test configuration (it was a few percent), and the recorded node accelerations to about $10^{-4}$, the int16 quantisation of the record included.

The acceleration that radiates is the difference of consecutive projected velocities. With the lagged projection, the difference of raw velocities would contain the full rigid acceleration $f_{rigid}/M$ at every step, since each step adds it and the next removes it; a plate that appears to accelerate as a whole radiates like a monopole, which is exactly the component that the front and back terms of an unbaffled plate are meant to cancel. This error was invisible in the nonlinear run and a factor 4.6 in the rms pressure of the linear reference run, where the cancellation is 25:1.

5. Buffers and precision#

The GPU keeps the set-up and the state in storage buffers, at most eight in any one pipeline (the default WebGPU limit): per-element data (CST operators, area, slopes, membrane stiffness, damping factor and the 45 entries of the upper triangle of $K_e$, 58 floats per element), element node indices, the state (q, v and the previous projected $(u, v, w)$ velocities, 13 floats per node), the element forces, per-node data (inverse mass, mass, the position relative to the rigid centre, the rest height and the nodal normal, 16 floats per node, with the decimation taps at the tail), the adjacency in CSR form, the hits (static patch data and the exciter state, one float and one integer buffer), the synchronization buffer (step counters and partial sums), the radiation buffer (the $W N$ accumulator slots, the node history, the current sample's $M N$ sub-step accelerations, and the chunk's acceleration output rows when recording), the element history, the radiation coefficients (delays and interpolator taps, 36 floats per element) and the chunk's pressure output. The per-element and per-node data are stored field-major (field $f$ of element $e$ at $f E + e$), so that the threads of a workgroup read consecutive words; $K_e$, symmetric to round-off, is symmetrised in float64 before packing, and the product $K_e U_e$ is written out term by term in the shader. Together these halved the time of the element pass on an integrated GPU. A uniform holds the chunk's sample and frame bases. Everything on the GPU is float32. The set-up is float64 and is packed once. Diagnostics (energies, exciter state) are read back between chunks; the plate view reads $w$ from the state buffer directly in its vertex shader without a readback. The whole integrator state (state, sync, hits, radiation and element-history buffers, plus the step count) can be read back and written again, which is how a partial .cymsim file is continued.

Float32 is adequate for the leapfrog at these step sizes: the displacement increment per step is about $10^{-3}$ of the displacement, so the update loses a few bits per step in $q$ without a visible effect on the energy balance, and the run is chaotic anyway once max|w|/h exceeds a few units, so no float32 trajectory tracks a float64 one for long. What is checked is the deterministic part (section 6).

6. Verification#

The tests compare three implementations pairwise on a coarse dent-free 18 inch profile (30 mm elements, 227 nodes, ten sub-steps per sample) with a 9.6 J mallet stroke, in the nonlinear and in the linear (reference) setting.

TypeScript set-up against Python: the mesh is identical to round-off; the DKT stiffness agrees with the autodiff element to $10^{-9}$; mass, stability bound, limiting element, element stiffnesses, lumped masses and rigid modes agree to $10^{-8}$ or better; the contact patch (nodes, weights, offsets, stiffness, free contact time) agrees to $10^{-6}$.

CPU loop against Python (float64 both): in the linear run the peak force and rebound velocity agree to $10^{-5}$, the energies and max|w|/h to $10^{-4}$; in the nonlinear run, where the trajectory is sensitive, the same quantities agree to $10^{-3}$. The unbaffled and baffled pressures agree to $2$ to $3 \times 10^{-4}$ relative rms in both runs, and the recorded normal accelerations of individual nodes to about $10^{-4}$ (the test requires $10^{-3}$); the only remaining differences are the delay applied at sub-step resolution in Python against the 16-tap interpolator here, and the int16 quantisation of the record. The decimation taps agree to $10^{-16}$.

GPU kernels against the CPU loop (float32 against float64, under Dawn): peak force to $10^{-4}$, contact time to the step, energies to $10^{-3}$ and the pressure to $10^{-3}$ relative rms in the linear run; $10^{-3}$, $10^{-2}$ and $10^{-2}$ in a short nonlinear run.

7. The displacement movie#

The plate view reads w straight from the simulation's state buffer while a run is in progress, which shows the motion at whatever rate the browser repaints (a few frames per second of wall time, each an arbitrary instant of simulated time). To watch the motion itself, the run also captures frames at a prescribed rate in simulated time: a capture dispatch runs at the start of every sub-step whose audio-sample index is a multiple of movieEvery = round(fs / fps), writing w of every node into a per-chunk buffer that is read back with the chunk's pressure. The CPU loop captures the same instants, and the GPU check compares the two frame by frame (relative rms below 1e-2 in float32).

Frames are stored quantized to 16 bits against each frame's own peak, which costs 2 bytes per node per frame and keeps the relative error near 1e-4 of that frame's amplitude, far below what the colour map and the height exaggeration resolve. The default 2000 fps resolves the low modes of a crash (a few hundred hertz) with several samples per cycle; the contact itself, a few milliseconds long, needs 5000 fps or more. Replay is driven by requestAnimationFrame against wall time divided by the slowdown factor, so the same frames play correctly on any display: the player advances simulated time by the elapsed wall time over the factor and shows the nearest captured frame, rather than stepping one frame per repaint.

8. The .cymsim file#

A run can be saved at any time, complete or not, as a .cymsim file, an HDF5 file whose layout is specified in cymsim's docs/format.md. It holds the configuration as YAML, the mesh with its nodal normals, the listener geometry, the two listener signals, the transverse displacement at a low frame rate (output.displacement_rate, default 60 Hz; a strobe for display, not a signal), optionally the normal acceleration of every node at the audio rate (output.record_acceleration, default off), the energy log, the hit summaries, and this application's integrator state. The file is written with h5wasm (libhdf5 compiled to WebAssembly, loaded on first use) into an in-memory image that is then downloaded, so at export time the memory in use briefly doubles; a running simulation is paused between chunks while the file is built, so that the state and the records are consistent, and continues afterwards.

The acceleration record is the reason for the file: it is what another program needs to compute the sound of the simulated cymbal for listeners, rooms or radiation models other than the built-in one, and it must not be derived from stored displacements (a second difference multiplies noise by $\omega^2$ and folds the shell's ultrasonic content into the band). It is the same audio-rate normal acceleration that feeds the listener sum (section 4), quantised to int16 against one float32 scale per sample, so that a stroke's 90 dB of decay keeps its precision; the relative error of a node is at most $1.5 \times 10^{-5}$ of the sample's peak. The cost is $2 N f_s$ bytes per simulated second, about 0.4 GB per second at the default mesh, kept in memory in blocks as the chunks arrive; the app estimates the size before a run starts, warns above 1 GB and refuses above half the device memory (or 2 GB when the browser does not say). The chunk size is reduced so that a chunk's rows stay under 64 MB on the GPU. The Python test suite rebuilds the baffled listener signal from such a record alone, with a band-limited fractional delay, to $1.5 \times 10^{-4}$ relative rms.

Loading a file restores the editor, the spectrogram and playback, the displacement frames (which also serve as the replay movie, at the file's frame rate), the energy and hit tables and the acceleration record, so that a loaded run can be exported again. A partial file written by this application on the same backend (GPU or CPU) can be continued with Resume, which rebuilds the set-up from the configuration, writes the saved state back into the buffers and runs on, appending to the signals and records; a float32 state reproduces the uninterrupted run only to float32 precision, which in the chaotic regime means the same statistics rather than the same waveform. Files written by the Python package load in the same way; their state group belongs to the Python integrator and is ignored, and their runs cannot be continued here.

9. Limitations#

The limitations of the model are those listed in the Python documentation. Specific to the web application: the imperfections of a seed differ from Python's; the decimation and the fractional delays are approximations described in section 4; the plate view shows $w$ scaled by a user-chosen factor and is not a measurement; the movie is capped at 256 MB, after which capture stops silently apart from a note in the toolbar, and its frames are quantized as described in section 7; and the CPU fallback is a hundred times slower than a discrete GPU. The set-up runs on the main thread and freezes the page for a second or two at full resolution.