The .cymsim file format#

A .cymsim file is an HDF5 file that holds one simulation run: the configuration that produced it, the mesh, the listener geometry, the two listener signals, the transverse displacement at a low frame rate for display, optionally the normal acceleration of every node at the audio rate, and, when written by the Python package, the integrator state needed to continue the run. The same format is written and read by the Python package (cymsim.save_cymsim, cymsim.load_cymsim, cymsim run -o) and by the web companion cymsimweb. This document is the specification; the Python implementation is in cymsim_file.py.

The purpose of the acceleration record is to let another program compute the sound of the simulated cymbal for listeners, rooms or radiation models other than the one built into cymsim, for instance by solving the acoustic wave equation with the shell's normal acceleration as a boundary condition. The record is therefore defined so that it needs no further differentiation, and its frequency response is stored with it.

1. Conventions#

All quantities are SI (metres, seconds, kilograms, pascals). The cymbal centre is at the origin, the mean rim plane is $z = 0$ and the bell points to $+z$. The front of the cymbal is its $+z$ face, the face the listener in the default configuration sees. Time series along the first axis start at $t = 0$, the release of the first exciter being at hits[0].time of the configuration. Every dataset that grows with time is resizable along that axis (HDF5 maxshape unlimited), so that a file can be written during a run and extended later. A partial file, written while the run was still in progress, is a valid file whose time series are prefixes of those of the complete file; run.complete says which it is.

Integers are stored as int32, geometry as float64, signals (audio, probes, displacement frames) in the precision of the run, float32 or, for a run with numerics.dtype: float64, float64, and the acceleration record as int16 against a scale per sample (section 5). Strings are variable-length UTF-8. Attributes are used for scalars and short vectors, datasets for arrays.

2. Root#

3. Configuration: config/yaml#

A scalar string dataset holding the run's configuration as YAML text, in the schema of SimulationConfig (the same file that cymsim run and the web app accept). It is the authoritative description of the run; the other groups are derived from it and from the run. Note that output.record_acceleration and output.displacement_rate in this text say which records the writer was asked to keep.

4. Geometry: mesh and listener#

mesh (attributes n_nodes $N$, n_elements $E$):

The nodal normal is the area-weighted mean of the constant slopes of the adjacent elements, $n = (-z_{0,x}, -z_{0,y}, 1) / \sqrt{1 + z_{0,x}^2 + z_{0,y}^2}$. The rest surface of element $e$ is the flat triangle through its three nodes at heights z0.

listener (attributes position $(3,)$, c_air, rho_air, back_diffraction_gain, copied from the configuration):

These are what the built-in radiation model used; a reader that wants the built-in listener signal from the acceleration record needs r_front, area, rho_air and c_air (section 9).

5. Run parameters: run#

Attributes only:

The integrator runs $W/2$ audio samples past duration so that the last samples of the records are final (section 6).

6. Acceleration record: acceleration (optional)#

Present only when the writer was asked to keep it (output.record_acceleration: true). Group attributes: fs, component = "normal", units = "m/s^2", encoding = "int16 * scale[t]", description.

$a_{ij}$ is the acceleration of node $i$ along the nodal normal mesh/normal (front side positive), band-limited to the audio band and sampled at $t_j = j / f_s$. It is not a second difference of stored displacements. In the leapfrog scheme the integrator's velocity increment $(v^{n+1} - v^n) / \Delta t$ is exactly the central second difference of $q$ at sub-step $n$, and this quantity, projected on the normal, is what the time loop filters. The filter is a Kaiser-windowed sinc $h(d)$ over the sub-step offsets $d = -Wm/2, \ldots, Wm/2 - 1$ from the sample centre, normalised to unit DC gain, with $W = 32$, Kaiser $\beta = 8$ and a $-6$ dB point at $0.46 f_s$. Its response is flat to about $0.38 f_s$ (17 kHz at 44.1 kHz), $-0.2$ dB at $0.41 f_s$, $-23$ dB at $f_s / 2$ and below $-80$ dB from $0.54 f_s$; at 44.1 kHz, content of the shell between 22 and 24 kHz therefore folds back into 20 to 22 kHz with 23 to 80 dB of attenuation, and nothing above that survives. A run at a higher numerics.fs moves this edge up. The kernel is linear-phase and its delay is removed, so sample $j$ is centred on $t_j$ exactly. It is stored as taps[w, kk] $= h(kk + (W/2 - 1 - w)\, m)$, the weight with which sub-step $j_0 m + kk$ enters sample $j_0 + w - W/2 + 1$; a reader that wants $h$ itself recovers it as taps[W/2 - 1 - (d div m), d mod m] for $d$ in the range above, or simply reads off the frequency response from the $Wm$ values.

The quantisation is one int16 per node against one float32 scale per sample, $\mathrm{scale}_j = \max_i |a_{ij}| / 32767$ (zero for an all-zero sample). The relative error of a node is at most $1.5 \times 10^{-5}$ of the sample's peak, which is far below the errors of the shell model and adds up incoherently in any radiation sum; in the built-in Rayleigh sum the quantisation changes the pressure by less than $10^{-4}$ relative rms. Storing the peak per sample rather than per file matters because a stroke spans several orders of magnitude within milliseconds.

The dataset is time-major with chunks of 1024 samples by all $N$ nodes (fewer samples per chunk for very large meshes, so that a chunk stays below 16 MB), uncompressed: a time window is one contiguous read, and the quantised data does not compress usefully. A partial file holds the first $T$ samples, all of them final: sample $j$ is final once the run has passed sub-step $(j + W/2)\, m$. The size is $2 N f_s$ bytes per simulated second, about 0.4 GB per second for the default 4657-node mesh at 44.1 kHz.

7. Displacement frames: displacement (optional)#

Present when output.displacement_rate is positive (default 60 Hz). Group attributes: rate, units = "m", description.

Frame $k$ is taken at the sub-step nearest to $k / \mathrm{rate}$; time gives the exact sub-step time. These frames are a strobe of the motion for display, not a signal: at 60 Hz they undersample every mode of the cymbal. They do carry what the acceleration record cannot, the slow quasi-static deflection, and they are cheap ($4 N \cdot \mathrm{rate}$ bytes per second, about 1.1 MB per second at the default mesh; twice that for a float64 run). The web app's slow-motion replay uses its own capture at a much higher frame rate and does not write it to the file.

8. Listener signals: audio#

Group attributes fs, units = "Pa".

Both are the pressure at listener/position, sampled at $t_j = j / f_s$, filtered by the same kernel as the acceleration record. In a partial file $T_a$ is the number of samples that were final at the time of writing. The Python package finalises pressure and acceleration samples on the same schedule, so $T_a = T$ in its files, but a reader should not rely on that.

probes (optional, present when output.probes is not empty; attributes fs, units = "m/s^2"):

9. Reproducing the baffled signal from the record#

A reader that has only acceleration, mesh and listener recovers the baffled signal as

where the fractional delay $r_e f_s / c$ in samples must be applied with a band-limited interpolator (a windowed sinc of a few dozen taps; linear interpolation attenuates the top octave). With a 33-tap Kaiser-windowed sinc the Python test suite reproduces the stored baffled signal to $1.5 \times 10^{-4}$ relative rms on a coarse run. The remaining difference comes from the built-in model applying the delay at sub-step resolution before filtering, and from the interpolator's own response. The back-face term is the same sum with r_back in place of r_front and the sign reversed: the normal acceleration of the back face is the negative of the front's.

10. Diagnostics: energy and hits#

energy/table is a $(N_e, 6)$ float64 dataset with the attribute columns = "time kinetic membrane bending exciters max_w_over_h" (energies in joules, max_w_over_h the nonlinearity meter $\max_i |w_i| / h_i$), one row per progress report. hits holds one entry per configured hit in the datasets time, kind (fixed-length string, "stick" or "mallet"), peak_force [N], contact_time [s], rebound_velocity [m/s, positive upward, as of the last step written], energy_in and energy_out [J].

11. Integrator state: state (optional, writer-specific)#

The state needed to continue the run from n_steps_done. Its content depends on the program that wrote it, named in the group attribute writer; readers of the format may ignore the group, and a program should use only a state written by itself. The Python package (writer = "cymsim-python", attributes n_steps_done, dtype) stores q $(N, 5)$ and v $(N, 5)$, the nodal displacements $(u, v, w, \theta_x, \theta_y)$ and velocities; exciter_z, exciter_v, peak_force, contact_steps per hit; ring_front, ring_back, the pending retarded-time pressure contributions at the sub-step rate; decimator_tail_front, decimator_tail_back, the sub-steps still needed by audio samples that are not yet final; and, when recording, acceleration_ring $(W, N)$, the partial sums of the not-yet-final acceleration samples; decimator_tail_probes $(P, L)$, the decimator tails of the probe signals; and the overhang, the samples past duration that were already final when the run stopped (the integrator overshoots the duration by $W/2$ audio samples and by the rounding to whole blocks): overhang_front, overhang_unbaffled (float64), overhang_probes, overhang_acceleration_data, overhang_acceleration_scale, overhang_displacement_time and overhang_displacement_w, each empty in a partial file. cymsim run --resume continues a run from this group, for the same configuration with a longer duration (or a shorter one, which cuts the records); the continued run is identical to an uninterrupted one on the same device. A state without the overhang datasets was written before runs could be continued and cannot be resumed. The web companion stores its own float32 GPU state under writer = "cymsimweb".

12. Tangent: tangent (optional)#

Present for a run with numerics.tangent: true. The run then also integrates its tangent: the derivative of the state with respect to $\delta$ when every prescribed peak force is multiplied by $(1 + \delta)$, carried exactly by forward-mode differentiation of the time step (jax.jvp). It is the limit of $(\text{twin} - \text{run}) / \delta$ for a twin run whose peak forces are larger by the relative amount $\delta$, and does not saturate as a finite twin does. The group holds the tangent's outputs in the layout of the run's own: audio/baffled, audio/unbaffled (Pa per unit $\delta$), probes/normal_acceleration (when there are probes), displacement/w and displacement/time (when frames are kept; the rate is that of section 7); norms $(B, 4)$, one row after every block, with columns (attribute columns) time tangent_norm2 kinetic linear_strain: the squared energy norm of the tangent state, $T(\dot{q}') + U_\text{lin}(q')$, and the two parts of the same norm of the run's state, $T(\dot{q})$ and $U_\text{lin}(q)$, where $T$ is the kinetic energy and $U_\text{lin}$ the strain energy without the von Kármán terms (a quadratic form); and state, the tangent's integrator state in the layout of section 11, from which a tangent run continues. The relative separation of an infinitesimal twin is $\delta \sqrt{\texttt{tangent\_norm2} / (\texttt{kinetic} + \texttt{linear\_strain})}$.

13. Reading the file elsewhere#

Any HDF5 reader works. In Python without cymsim:

import h5py, numpy as np
with h5py.File("run.cymsim") as f:
    fs = f["run"].attrs["fs"]
    a = f["acceleration/normal"][1000:2000].astype(np.float32) * f["acceleration/scale"][1000:2000, None]
    xy, z0, n = f["mesh/xy"][...], f["mesh/z0"][...], f["mesh/normal"][...]

cymsim info run.cymsim prints a summary, cymsim.load_cymsim returns a SimulationResult whose acceleration.as_float(start, stop) decodes a window of the record.