1# The .cymsim file format
3A `.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](../src/cymsim/cymsim_file.py).
5The 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.
7## 1. Conventions
9All 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.
11Integers are stored as `int32`, geometry as `float64`, signals as `float32`, 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.
13## 2. Root
15| attribute | type | meaning |
16|---|---|---|
17| `format` | string | always `"cymsim"` |
18| `format_version` | int | `1` for the layout described here |
19| `creator` | string | program and version, e.g. `"cymsim 0.1.0"` or `"cymsimweb 0.1.0"` |
20| `created` | string | ISO 8601 time of the first write, UTC |
22## 3. Configuration: `config/yaml`
24A 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.
26## 4. Geometry: `mesh` and `listener`
28`mesh` (attributes `n_nodes` $N$, `n_elements` $E$):
30| dataset | shape | meaning |
31|---|---|---|
32| `xy` | $(N, 2)$ float64 | node positions in the projected plane |
33| `z0` | $(N,)$ float64 | rest height of the mid-surface at the nodes |
34| `thickness` | $(N,)$ float64 | shell thickness at the nodes |
35| `normal` | $(N, 3)$ float64 | unit normal of the rest surface at the nodes, pointing to the front ($+z$ component positive) |
36| `triangles` | $(E, 3)$ int32 | node indices of each element, counter-clockwise in the plane |
37| `area` | $(E,)$ float64 | projected element areas |
38| `boundary` | $(N_b,)$ int32 | indices of the rim nodes in rim order |
40The 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`.
42`listener` (attributes `position` $(3,)$, `c_air`, `rho_air`, `back_diffraction_gain`, copied from the configuration):
44| dataset | shape | meaning |
45|---|---|---|
46| `r_front` | $(E,)$ float64 | distance from the element centroid on the rest surface to the listener |
47| `r_back` | $(E,)$ float64 | shortest path from the centroid over the rim to the listener |
49These 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).
51## 5. Run parameters: `run`
53Attributes only:
55| attribute | meaning |
56|---|---|
57| `fs` | audio sample rate [Hz]; every audio-rate series is sampled at $t_j = j / f_s$ |
58| `dt` | integrator sub-step [s], $\Delta t = 1 / (m f_s)$ |
59| `substeps` | $m$, sub-steps per audio sample |
60| `duration` | requested simulated time [s] (`numerics.duration`) |
61| `n_samples` | $\mathrm{round}(f_s \cdot \mathrm{duration})$, the length every audio-rate series reaches when complete |
62| `n_steps_done` | sub-steps integrated so far |
63| `complete` | integer flag, 1 when every time series has its full length, 0 for a partial file (an integer rather than an HDF5 boolean so that every HDF5 client reads it) |
64| `decimation_span` | $W$, the number of audio samples spanned by the anti-aliasing kernel (section 6) |
65| `wall_time` | seconds of wall time spent in the time loop (informative) |
67The integrator runs $W/2$ audio samples past `duration` so that the last samples of the records are final (section 6).
69## 6. Acceleration record: `acceleration` (optional)
71Present 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`.
73| dataset | shape | meaning |
74|---|---|---|
75| `normal` | $(T, N)$ int16 | quantised normal acceleration; sample $j$, node $i$ |
76| `scale` | $(T,)$ float32 | scale of each sample: $a_{ij} = \mathrm{normal}_{ij} \cdot \mathrm{scale}_j$ |
77| `taps` | $(W, m)$ float64 | the decimation kernel in polyphase layout (below) |
79$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.
81The 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.
83The 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.
85## 7. Displacement frames: `displacement` (optional)
87Present when `output.displacement_rate` is positive (default 60 Hz). Group attributes: `rate`, `units = "m"`, `description`.
89| dataset | shape | meaning |
90|---|---|---|
91| `w` | $(T_d, N)$ float32 | transverse displacement $w$ of every node at the frame times |
92| `time` | $(T_d,)$ float64 | frame times [s] |
94Frame $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). 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.
96## 8. Listener signals: `audio`
98Group attributes `fs`, `units = "Pa"`.
100| dataset | shape | meaning |
101|---|---|---|
102| `baffled` | $(T_a,)$ float32 | Rayleigh integral of the front face, as if in an infinite baffle |
103| `unbaffled` | $(T_a,)$ float32 | front minus `back_diffraction_gain` times the diffracted back-face term |
105Both 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.
107`probes` (optional, present when `output.probes` is not empty; attributes `fs`, `units = "m/s^2"`):
109| dataset | shape | meaning |
110|---|---|---|
111| `xy` | $(P, 2)$ float64 | the probe points in the projected plane |
112| `normal_acceleration` | $(T_p, P)$ float32 | normal acceleration at each point, interpolated linearly in the containing element and filtered by the same kernel as the listener signals |
114## 9. Reproducing the baffled signal from the record
116A reader that has only `acceleration`, `mesh` and `listener` recovers the baffled signal as
118$$
119p_j = \frac{\rho_{air}}{2 \pi} \sum_e \frac{A_e}{r_e}\, \bar a_e\!\left(t_j - r_e / c\right), \qquad \bar a_e = \tfrac{1}{3} \sum_{i \in e} a_i ,
120$$
122where 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.
124## 10. Diagnostics: `energy` and `hits`
126`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].
128## 11. Integrator state: `state` (optional, writer-specific)
130The 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"`.
132## 12. Reading the file elsewhere
134Any HDF5 reader works. In Python without cymsim:
136```python
137import h5py, numpy as np
138with h5py.File("run.cymsim") as f:
139 fs = f["run"].attrs["fs"]
140 a = f["acceleration/normal"][1000:2000].astype(np.float32) * f["acceleration/scale"][1000:2000, None]
141 xy, z0, n = f["mesh/xy"][...], f["mesh/z0"][...], f["mesh/normal"][...]
142```
144`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.