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, using the twelve coefficients $c_k = r_k^{\mathsf T} M v$ and $d_k = r_k^{\mathsf T} M q$ that the previous node pass summed; performs the leapfrog update $v' = (v + \Delta t M^{-1} f)/(1 + 2 \sigma_0 \Delta t)$, $q' = q + \Delta t v'$; and contributes $r_k^{\mathsf T} M v'$ and $r_k^{\mathsf T} M q'$ to a per-workgroup partial sum (a shared-memory tree reduction) for the next step. 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#

The Python package scatters every element's acceleration into retarded-time buffers at sub-step resolution and decimates the result to the audio rate with a polyphase FIR. 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: each element's acceleration is first decimated to the audio rate, then delayed and summed.

The decimation is a Blackman-windowed sinc with cutoff $0.45 f_s$ spanning eight audio samples, applied through eight accumulator slots per element: at sub-step phase $k$ of audio sample $j_0$ the element adds its acceleration times the tap $h_w(k)$ to slot $w$, which belongs to audio sample $j_0 + w - 3$; when a sample is complete its slot is written to a per-element history ring (of length a little more than the longest delay in audio samples) and the slots shift. A triangular window over two samples was tried first and was found insufficient in the linear reference run, whose sharp contact puts energy above the audio band. The gather is one dispatch of a single workgroup per audio sample: every element's history is read at $j - \lfloor d_e \rfloor$ and $j - \lfloor d_e \rfloor - 1$ with the fractional delay applied by linear interpolation, for the front path and for the path over the rim, and the two sums are reduced in shared memory and written to the chunk's output buffer. Linear interpolation of a fractional delay attenuates the top of the band slightly (up to a few dB at 20 kHz for a half-sample fraction, spread over the elements' different fractions), which we accept for an interactive tool.

The acceleration that radiates is the difference of consecutive projected transverse 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. With the projected differences the two implementations agree to 0.1 percent in rms in both runs.

5. Buffers and precision#

The GPU keeps eleven 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 81 entries of $K_e$, 96 floats per element, with the decimation taps at the tail), element node indices, the state (q, v, the radiated nodal acceleration and the previous projected velocity), the element forces, per-node data (inverse mass, mass and the six rigid modes, 40 floats per node), 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 state (accumulator slots and history) and the radiation coefficients. 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.

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}$, the rms of the unbaffled pressure to 0.1 percent; in the nonlinear run, where the trajectory is sensitive, the same quantities agree to $10^{-3}$ and the pressure to 0.2 percent rms with the peaks at the same sample.

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. 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; 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.