save .cymsim file
16 changed files+1172−79
.gitignoremodified+1−0View file
@@ -1,3 +1,4 @@
1+tmp
12 __pycache__/
23 *.egg-info/
34 .pytest_cache/
README.mdmodified+18−4View file
@@ -5,7 +5,7 @@ Explicit nonlinear finite-element simulation of struck cymbals. The cymbal is a
55 ## Installation
66
77 ```
8-pip install -e ".[dev]" # numpy, scipy, jax (CPU), pyyaml, pytest, matplotlib
8+pip install -e ".[dev]" # numpy, scipy, jax (CPU), pyyaml, h5py, pytest, matplotlib
99 pip install "jax[cuda12]" # optional: run the time loop on a GPU
1010 ```
1111
@@ -18,8 +18,13 @@ cymsim example-config my_run.yaml --exciter mallet # write a default config
1818 cymsim check my_run.yaml # mesh, time step, element self-tests, no time stepping
1919 cymsim simulate my_run.yaml # run; writes <prefix>_rayleigh.wav, _unbaffled.wav, _result.npz, _cascade.png
2020 cymsim simulate my_run.yaml --linear --duration 0.5 # linear reference run
21+cymsim run my_run.yaml -o run.cymsim --save-every 60 # also write a .cymsim file, updated every minute of wall time
22+cymsim run my_run.yaml -o run.cymsim --record-acceleration # include the audio-rate normal acceleration of every node
23+cymsim info run.cymsim # summarise a .cymsim file
2124 ```
2225
26+A `.cymsim` file is an HDF5 file holding the configuration, the mesh, the listener signals, displacement frames at a low rate, optionally the normal acceleration of every node at the audio rate, and the integrator state; the layout is specified in [docs/format.md](docs/format.md). It is meant for taking the shell's motion into other programs (other listeners, rooms, acoustic solvers) and is the file the web companion cymsimweb saves and loads. A partial file written during a run is a consistent prefix of the complete one.
27+
2328 The same from Python:
2429
2530 ```python
@@ -27,10 +32,17 @@ from cymsim import SimulationConfig, Hit, StickExciter, simulate
2732
2833 cfg = SimulationConfig.from_yaml("examples/crash_mallet.yaml")
2934 cfg.hits.append(Hit(time=0.8, r=0.19, phi=2.0, velocity=7.0, exciter=StickExciter()))
35+cfg.output.record_acceleration = True # off by default: about 0.4 GB per simulated second at the default mesh
3036 result = simulate(cfg)
3137 result.write_wav("crash")
38+result.save("crash.cymsim")
3239 print(result.hits) # peak force, contact time and rebound velocity of each hit
3340 print(result.energy_log) # kinetic / membrane / bending / exciter energy and max|w|/h over time
41+a = result.acceleration.as_float(0, 4410) # normal acceleration [m/s^2], first 0.1 s, shape (4410, n_nodes)
42+w = result.displacement.w # transverse displacement at 60 Hz, shape (n_frames, n_nodes)
43+
44+from cymsim import load_cymsim
45+back = load_cymsim("crash.cymsim") # the same SimulationResult, from the file
3446 ```
3547
3648 A run is described by one `SimulationConfig` made of `Material`, `Profile`, `Imperfections`, `MeshConfig`, a list of `Hit`s (each with a `StickExciter` or `MalletExciter`), `Listener`, `Numerics` and `Output`. The YAML file mirrors the dataclass fields; omitted fields take the defaults, which reproduce an 18 inch medium-thin B20 crash. Units are SI throughout. See [examples/](examples/) for a mallet stroke, a stick crash, a train of four hits and a quick CPU run.
@@ -46,9 +58,11 @@ A run is described by one `SimulationConfig` made of `Material`, `Profile`, `Imp
4658 | [shell.py](src/cymsim/shell.py) | assembled shell: lumped mass, rigid modes, internal force by autodiff, stability bound |
4759 | [exciter.py](src/cymsim/exciter.py) | contact patches, power-law contact law, schedule of hits |
4860 | [radiation.py](src/cymsim/radiation.py) | Rayleigh integral and the unbaffled (rim-diffraction) approximation |
49-| [simulate.py](src/cymsim/simulate.py) | leapfrog time loop, block compilation, result container |
61+| [decimation.py](src/cymsim/decimation.py) | the anti-aliasing kernel that takes accelerations and pressure from the sub-step rate to the audio rate |
62+| [simulate.py](src/cymsim/simulate.py) | leapfrog time loop, block compilation, acceleration and displacement records, result container |
63+| [cymsim_file.py](src/cymsim/cymsim_file.py) | the .cymsim HDF5 file: streaming writer, loader, summary |
5064 | [diagnostics.py](src/cymsim/diagnostics.py) | band envelopes, onset times, spectral centroid, cascade plot |
51-| [cli.py](src/cymsim/cli.py) | `cymsim simulate`, `cymsim check`, `cymsim example-config` |
65+| [cli.py](src/cymsim/cli.py) | `cymsim simulate` (`run`), `cymsim check`, `cymsim info`, `cymsim example-config` |
5266
5367 `prototype.py` is the original single-file prototype that this package reorganizes; it is kept for reference and is not used by the package.
5468
@@ -58,4 +72,4 @@ A run is described by one `SimulationConfig` made of `Material`, `Profile`, `Imp
5872 pytest
5973 ```
6074
61-The tests check the mesher (manifold, orientation, quality), the DKT element (rigid modes, constant curvature, twist), the membrane strains, the assembled shell (mass, rigid-mode annihilation, force = minus energy gradient, dissipative damping), the configuration round trip, and short coarse simulations with one and two hits.
75+The tests check the mesher (manifold, orientation, quality), the DKT element (rigid modes, constant curvature, twist), the membrane strains, the assembled shell (mass, rigid-mode annihilation, force = minus energy gradient, dissipative damping), the configuration round trip, the decimation kernel and streaming decimator, short coarse simulations with one and two hits, and the .cymsim file (round trip, partial files, and that the baffled listener signal is recovered from the stored accelerations alone).
docs/format.mdadded+137−0View file
@@ -0,0 +1,137 @@
1+# The .cymsim file format
2+
3+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](../src/cymsim/cymsim_file.py).
4+
5+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.
6+
7+## 1. Conventions
8+
9+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.
10+
11+Integers 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.
12+
13+## 2. Root
14+
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 |
21+
22+## 3. Configuration: `config/yaml`
23+
24+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.
25+
26+## 4. Geometry: `mesh` and `listener`
27+
28+`mesh` (attributes `n_nodes` $N$, `n_elements` $E$):
29+
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 |
39+
40+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`.
41+
42+`listener` (attributes `position` $(3,)$, `c_air`, `rho_air`, `back_diffraction_gain`, copied from the configuration):
43+
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 |
48+
49+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).
50+
51+## 5. Run parameters: `run`
52+
53+Attributes only:
54+
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) |
66+
67+The integrator runs $W/2$ audio samples past `duration` so that the last samples of the records are final (section 6).
68+
69+## 6. Acceleration record: `acceleration` (optional)
70+
71+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`.
72+
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) |
78+
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.
80+
81+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.
82+
83+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.
84+
85+## 7. Displacement frames: `displacement` (optional)
86+
87+Present when `output.displacement_rate` is positive (default 60 Hz). Group attributes: `rate`, `units = "m"`, `description`.
88+
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] |
93+
94+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). 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.
95+
96+## 8. Listener signals: `audio`
97+
98+Group attributes `fs`, `units = "Pa"`.
99+
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 |
104+
105+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.
106+
107+## 9. Reproducing the baffled signal from the record
108+
109+A reader that has only `acceleration`, `mesh` and `listener` recovers the baffled signal as
110+
111+$$
112+p_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 ,
113+$$
114+
115+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.
116+
117+## 10. Diagnostics: `energy` and `hits`
118+
119+`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].
120+
121+## 11. Integrator state: `state` (optional, writer-specific)
122+
123+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. The web companion stores its own float32 GPU state under `writer = "cymsimweb"`. Resuming a float32 run reproduces the uninterrupted run only to the precision of the state that was saved.
124+
125+## 12. Reading the file elsewhere
126+
127+Any HDF5 reader works. In Python without cymsim:
128+
129+```python
130+import h5py, numpy as np
131+with h5py.File("run.cymsim") as f:
132+ fs = f["run"].attrs["fs"]
133+ a = f["acceleration/normal"][1000:2000].astype(np.float32) * f["acceleration/scale"][1000:2000, None]
134+ xy, z0, n = f["mesh/xy"][...], f["mesh/z0"][...], f["mesh/normal"][...]
135+```
136+
137+`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.
docs/model.mdmodified+23−9View file
@@ -1,6 +1,6 @@
11 # The cymsim model: equations, discretization and method
22
3-This document states what cymsim solves and how. Section 2 gives the continuous model, section 3 the geometry, section 4 the spatial discretization, section 5 the time integration, section 6 the excitation, section 7 the radiation model and section 8 the diagnostics. Section 9 lists the limitations we are aware of. Each section names the module that implements it. Units are SI throughout.
3+This document states what cymsim solves and how. Section 2 gives the continuous model, section 3 the geometry, section 4 the spatial discretization, section 5 the time integration, section 6 the excitation, section 7 the radiation model, section 8 the diagnostics and section 9 the records kept for other programs. Section 10 lists the limitations we are aware of. Each section names the module that implements it. Units are SI throughout.
44
55 ## 1. Overview
66
@@ -150,7 +150,7 @@ where $P$ is the rigid-mode projection of section 5.3. This is the symplectic Eu
150150
151151 ### 5.2 Stability and the choice of the time step
152152
153-For the linear undamped system the scheme is stable when $\Delta t \leq 2 / \omega_{max}$, with $\omega_{max}^2$ the largest eigenvalue of $M^{-1} K$. Instead of solving that global eigenproblem we use the element bound (Irons; see Hughes 2000): $\lambda_{max}(M^{-1} K) \leq \max_e \lambda_{max}(M_e^{-1} K_e)$, where $K_e$ is the $15 \times 15$ tangent stiffness of the element at rest and $M_e$ its share of the lumped mass. The bound is computed in float64 at build time and the limiting element (the thinnest, near the rim) is reported. The nonlinear membrane stiffening raises the instantaneous frequencies during large motion, so we take only a fraction `cfl_safety` (default 0.6) of the bound. The sub-step count per audio sample is $m = \lceil 1 / (f_s\, c_{safety}\, \Delta t_{bound}) \rceil$ and $\Delta t = 1 / (m f_s)$, so that the radiation buffers decimate exactly to the audio rate. If the energy diagnostic grows during a run the step is too large; the remedy is a smaller `cfl_safety`. The run aborts if the energy becomes non-finite.
153+For the linear undamped system the scheme is stable when $\Delta t \leq 2 / \omega_{max}$, with $\omega_{max}^2$ the largest eigenvalue of $M^{-1} K$. Instead of solving that global eigenproblem we use the element bound (Irons; see Hughes 2000): $\lambda_{max}(M^{-1} K) \leq \max_e \lambda_{max}(M_e^{-1} K_e)$, where $K_e$ is the $15 \times 15$ tangent stiffness of the element at rest and $M_e$ its share of the lumped mass. The bound is computed in float64 at build time and the limiting element (the thinnest, near the rim) is reported. The nonlinear membrane stiffening raises the instantaneous frequencies during large motion, so we take only a fraction `cfl_safety` (default 0.6) of the bound. The sub-step count per audio sample is $m = \lceil 1 / (f_s\, c_{safety}\, \Delta t_{bound}) \rceil$ and $\Delta t = 1 / (m f_s)$, so that the radiation buffers decimate exactly to the audio rate. The run is extended by $W/2 = 16$ audio samples beyond the requested duration so that the last samples of the decimated records (section 9) are final. If the energy diagnostic grows during a run the step is too large; the remedy is a smaller `cfl_safety`. The run aborts if the energy becomes non-finite.
154154
155155 ### 5.3 Rigid-mode projection
156156
@@ -190,31 +190,43 @@ A run has a list of hits, each with its own time, position, velocity and exciter
190190
191191 ## 7. Radiation
192192
193-Implemented in [radiation.py](../src/cymsim/radiation.py); parameters in `Listener`.
193+Implemented in [radiation.py](../src/cymsim/radiation.py) and [decimation.py](../src/cymsim/decimation.py); parameters in `Listener`.
194194
195195 The pressure at the listener point $x_L$ from a baffled planar radiator is given by the Rayleigh integral, which we discretize with one point source per element:
196196
197197 $$
198-p(x_L, t) = \frac{\rho_{air}}{2\pi} \sum_e \frac{A_e\, \ddot w_e\!\left( t - r_e / c \right)}{r_e},
198+p(x_L, t) = \frac{\rho_{air}}{2\pi} \sum_e \frac{A_e\, \bar a_e\!\left( t - r_e / c \right)}{r_e},
199199 $$
200200
201-where $\ddot w_e$ is the mean of the three nodal accelerations and $r_e$ the distance from the element centroid, placed on the curved surface at height $z_0$, to the listener. This is exact for a rigid piston in an infinite baffle and an approximation here: the surface is shallow rather than flat, the acceleration is taken along $z$ rather than along the local normal, and the air does not load the shell (radiation damping is part of what $\sigma_0$ and $\sigma_1$ lump).
201+where $\bar a_e$ is the mean over the element's three nodes of the *normal acceleration* $a_i = n_i \cdot (\ddot u_i, \ddot v_i, \ddot w_i)$, $n_i$ being the unit normal of the rest surface at node $i$ (the area-weighted mean of the adjacent elements' constant slopes, pointing to the $+z$ side), and $r_e$ the distance from the element centroid, placed on the curved surface at height $z_0$, to the listener. The acceleration is taken from the integrator directly: for the leapfrog scheme $(v^{n+1} - v^n) / \Delta t$ is the exact central second difference of $q$ at step $n$, so no differentiation of stored data is involved. This is exact for a rigid piston in an infinite baffle and an approximation here: the surface is shallow rather than flat, and the air does not load the shell (radiation damping is part of what $\sigma_0$ and $\sigma_1$ lump). At the bell, where the slope reaches 0.39, the normal component differs from $\ddot w$ by several percent, which is why we use it rather than $\ddot w$.
202202
203203 A cymbal has no baffle. Its back face radiates with the opposite sign, and that sound reaches the listener by diffracting around the rim; at low frequency the two contributions nearly cancel, which is why an unbaffled plate is a weak low-frequency source. We estimate the back contribution with the same sum but with $r_e$ replaced by the shortest path over the rim, $r_e^{back} = \min_j \big( |x_e - x_j^{rim}| + |x_j^{rim} - x_L| \big)$ over 720 rim points, weighted by `back_diffraction_gain`, and subtract it: $p_{unbaffled} = p_{front} - g\, p_{back}$. This is a Kirchhoff-style edge-diffraction estimate, not a boundary-element solution, and it ignores shadowing of the front by the shell itself. Both signals are returned; the baffled one is a cleaner measure of the shell motion, the unbaffled one sounds more like a cymbal in a room.
204204
205-Retarded times are handled at sub-step resolution: the delay of element $e$ in sub-steps, $d_e = r_e / (c \Delta t)$, is split into an integer part and a fraction, and each element's contribution is added to two adjacent slots of a delay buffer with linear interpolation. At the end the buffers are decimated to the audio rate with a polyphase FIR resampler (`scipy.signal.resample_poly`), which also removes the sub-step-rate content above the audio Nyquist frequency.
205+Retarded times are handled at sub-step resolution: the delay of element $e$ in sub-steps, $d_e = r_e / (c \Delta t)$, is split into an integer part and a fraction, and each element's contribution is added to two adjacent slots of a ring buffer with linear interpolation. After each block the host reads the slots of the finished sub-steps, clears them and decimates them to the audio rate with the kernel of section 9, the same kernel that filters the acceleration record. Since delay and filtering are both linear and time-invariant, the stored listener signal follows from the stored accelerations up to the linear-interpolation error at the sub-step scale, which is negligible; the test suite rebuilds the baffled signal from the record to $1.5 \times 10^{-4}$ relative rms.
206206
207207 ## 8. Diagnostics
208208
209209 Implemented in [simulate.py](../src/cymsim/simulate.py) (energies) and [diagnostics.py](../src/cymsim/diagnostics.py) (spectra). During the run the kinetic, membrane, bending and exciter energies are logged at `n_progress` points together with $\max_i |w_i| / h_i$, the *nonlinearity meter*: below about 1 the run is effectively linear, at 3 to 5 the cascade is clearly developed. The sum of the four energies can only decrease between the release of a hit and the next one, which is a direct check of the integrator. After the run, an STFT of the unbaffled pressure (256-point windows, 32-sample hop, about 0.7 ms) gives per-octave band energy envelopes from 60 Hz to 16 kHz; the *onset* of a band is the first time it reaches 20 dB below its own peak, and the spectral centroid is tracked against time. In a cascade the onsets are staggered with the high bands last and the centroid climbs over tens of milliseconds; in the linear reference run the bands rise together and the centroid is flat.
210210
211-## 9. Limitations and future work
211+## 9. Records and the .cymsim file
212+
213+Implemented in [decimation.py](../src/cymsim/decimation.py), [simulate.py](../src/cymsim/simulate.py) and [cymsim_file.py](../src/cymsim/cymsim_file.py); parameters in `Output`; file layout in [format.md](format.md).
214+
215+The listener model of section 7 is deliberately simple, and a natural use of the simulation is to feed the shell's motion to a different acoustic model: another listener position, a room, or a solver of the acoustic wave equation with the shell as a vibrating boundary. What such a model needs is the normal acceleration of the surface at the audio rate, and the question is how to obtain it without noise. Differentiating stored displacements twice is the wrong route: a second difference multiplies noise by $\omega^2$, so the top of the band suffers most, it folds whatever the shell does above the audio Nyquist frequency back into the band, and in float32 the three-term difference of nearly equal numbers loses most of its digits. We therefore take the acceleration from the integrator, where it is exact for the scheme (section 7), and the only remaining step is the reduction from the sub-step rate, typically 2.4 MHz, to the audio rate.
216+
217+That reduction has to be a proper low-pass filter, because the cascade puts real energy above 20 kHz and acceleration weights it by $\omega^2$. We use a Kaiser-windowed sinc $h(d)$ spanning $W = 32$ audio samples ($Wm$ sub-steps) with $\beta = 8$ and its $-6$ dB point at $0.46 f_s$: flat to about $0.38 f_s$, $-23$ dB at $f_s / 2$, below $-80$ dB from $0.54 f_s$. It is applied inside the compiled time loop through $W$ accumulator slots per node: at sub-step $j_0 m + kk$ the node adds $a_i\, h(kk + (W/2 - 1 - w)\, m)$ to slot $w$, which belongs to audio sample $j_0 + w - W/2 + 1$; when the sub-step count crosses a sample boundary the oldest slot is complete, is written to the block's output and is reused for the newest sample. The cost is $W$ multiply-adds per node per sub-step, small next to the element force evaluation. The kernel is linear-phase and the read-out removes its delay, so sample $j$ is centred on $t = j / f_s$; the price is that the newest $W/2$ samples are not final, which is why the run is extended by $W/2$ samples (section 5.2) and why a partial save stops $W/2$ samples short of the current time.
218+
219+The record is kept only on request (`output.record_acceleration`, default off), because it is large: one value per node per audio sample, stored as `int16` against one `float32` scale per sample, is $2 N f_s$ bytes per simulated second, about 0.4 GB per second at the default mesh. The per-sample scale keeps the quantisation noise at $1.5 \times 10^{-5}$ of each sample's own peak through the 90 dB or more of a stroke's decay; the noise is incoherent across nodes and changes the built-in Rayleigh sum by less than $10^{-4}$ relative rms. Transverse displacements, in contrast, are kept only at a low frame rate (`output.displacement_rate`, default 60 Hz), in float32 and exactly at the frame times. At 60 Hz they are a strobe of the motion, useful for display and for the slow quasi-static deflection, but not a signal from which anything should be differentiated.
220+
221+Both records, the two listener signals, the configuration, the mesh with its nodal normals, the energy log, the hit summaries and the integrator state are written to an HDF5 file with the extension `.cymsim`, described in [format.md](format.md). The time series are resizable, so `cymsim run -o run.cymsim --save-every 60` writes a consistent partial file every minute and extends it, and a run that is interrupted leaves a loadable file. The web companion writes and reads the same format.
222+
223+## 10. Limitations and future work
212224
213225 The model is a shallow shell, so it is only as good as the assumption $|\nabla z_0| \ll 1$; the default profile reaches a slope of 0.39 at the bell, which is at the edge of what shallow-shell theory is meant for, and a steep bell would call for a proper curved-shell formulation. The elements are flat facets with the curvature entering through $z_0$ alone, and the membrane strain uses the linear interpolation of $w$ while bending uses DKT rotations; this is a standard pairing but not a fully consistent one, and its convergence with element size has not been studied here. The lumped rotational inertia is inflated for the sake of the time step (section 4.6). The support is a rigid-mode projection rather than a model of the stand. The damping has two parameters and no physics behind them. The contact is one-dimensional along $z$, with no friction, no tangential force, no hysteresis and no hand. The radiation is a Rayleigh sum with a heuristic back-side correction.
214226
215227 The natural next steps are, roughly in order of expected benefit: a convergence study of the cascade statistics with element size and precision; a bell support with local stiffness and damping in place of the projection; a contact model with hysteretic losses for felt; and a proper unbaffled radiation model. A train of hits already works structurally, but the behaviour of the contact model against a strongly moving surface deserves its own validation.
216228
217-## 10. References
229+## 11. References
218230
219231 - Batoz, J.-L., Bathe, K.-J., Ho, L.-W. (1980). A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15, 1771-1812.
220232 - Marguerre, K. (1938). Zur Theorie der gekrummten Platte grosser Formanderung. Proceedings of the 5th International Congress for Applied Mechanics, 93-101.
@@ -226,7 +238,7 @@ The natural next steps are, roughly in order of expected benefit: a convergence
226238 - Hairer, E., Lubich, C., Wanner, G. (2006). Geometric Numerical Integration. Springer. (Symplectic Euler and Stormer-Verlet.)
227239 - Pierce, A. D. (1989). Acoustics: An Introduction to Its Physical Principles and Applications. Acoustical Society of America. (Rayleigh integral.)
228240
229-## 11. Symbols and where to find them
241+## 12. Symbols and where to find them
230242
231243 | symbol | meaning | configuration / code |
232244 |---|---|---|
@@ -241,3 +253,5 @@ The natural next steps are, roughly in order of expected benefit: a convergence
241253 | $k, \alpha, P$ | contact stiffness, exponent, springs per patch | `ContactPatch.k_total`, `.alpha`, `.n_points` |
242254 | $z_s, v_s$ | exciter face height and velocity | state entries `zs`, `vs` |
243255 | $p_{front}, p_{back}$ | Rayleigh front and diffracted back pressure | `SimulationResult.pressure_rayleigh`, `.pressure_unbaffled` |
256+| $n_i$, $a_i$ | nodal normal of the rest surface, normal acceleration | `ShellModel.normal_node`, `SimulationResult.acceleration` |
257+| $h(d)$, $W$ | decimation kernel, its span in audio samples | `decimation.decimation_kernel`, `DECIMATION_SPAN` |
pyproject.tomlmodified+1−1View file
@@ -8,7 +8,7 @@ version = "0.1.0"
88 description = "Explicit nonlinear finite-element simulation of struck cymbals"
99 readme = "README.md"
1010 requires-python = ">=3.10"
11-dependencies = ["numpy", "scipy", "jax", "pyyaml"]
11+dependencies = ["numpy", "scipy", "jax", "pyyaml", "h5py"]
1212
1313 [project.optional-dependencies]
1414 plot = ["matplotlib"]
src/cymsim/__init__.pymodified+11−2View file
@@ -13,6 +13,7 @@ Typical use::
1313 cfg = SimulationConfig.from_yaml("examples/crash_mallet.yaml")
1414 result = simulate(cfg)
1515 result.write_wav("crash")
16+ result.save("crash.cymsim") # HDF5 file, see docs/format.md; load_cymsim reads it back
1617 """
1718 import jax
1819
@@ -35,12 +36,20 @@ from .config import ( # noqa: E402
3536 SimulationConfig,
3637 StickExciter,
3738 )
38-from .simulate import SimulationResult, Simulator, simulate # noqa: E402
39+from .cymsim_file import describe, load_cymsim, save_cymsim # noqa: E402
40+from .simulate import ( # noqa: E402
41+ AccelerationRecord,
42+ DisplacementRecord,
43+ SimulationResult,
44+ Simulator,
45+ simulate,
46+)
3947
4048 __version__ = "0.1.0"
4149
4250 __all__ = [
4351 "Hit", "Imperfections", "Listener", "MalletExciter", "Material", "MeshConfig", "Numerics",
4452 "Output", "Profile", "SimulationConfig", "StickExciter",
45- "SimulationResult", "Simulator", "simulate",
53+ "AccelerationRecord", "DisplacementRecord", "SimulationResult", "Simulator", "simulate",
54+ "describe", "load_cymsim", "save_cymsim",
4655 ]
src/cymsim/cli.pymodified+24−3View file
@@ -1,4 +1,5 @@
1-"""Command line interface: ``cymsim simulate``, ``cymsim check`` and ``cymsim example-config``."""
1+"""Command line interface: ``cymsim simulate`` (alias ``run``), ``cymsim check``, ``cymsim info`` and
2+``cymsim example-config``."""
23 from __future__ import annotations
34
45 import argparse
@@ -21,6 +22,8 @@ def _load(path: str, args) -> SimulationConfig:
2122 cfg.output.prefix = args.out_prefix
2223 if getattr(args, "linear", False):
2324 cfg.numerics.nonlinear = False
25+ if getattr(args, "record_acceleration", False):
26+ cfg.output.record_acceleration = True
2427 return cfg
2528
2629
@@ -30,7 +33,7 @@ def cmd_simulate(args) -> int:
3033
3134 cfg = _load(args.config, args)
3235 sim = Simulator(cfg)
33- result = sim.run()
36+ result = sim.run(save_path=args.output, save_every=args.save_every)
3437 prefix = cfg.output.prefix if cfg.numerics.nonlinear else cfg.output.prefix + "_linear"
3538 if cfg.output.wav:
3639 print("wrote " + ", ".join(result.write_wav(prefix)))
@@ -63,6 +66,14 @@ def cmd_check(args) -> int:
6366 return 0
6467
6568
69+def cmd_info(args) -> int:
70+ from .cymsim_file import describe
71+
72+ for path in args.files:
73+ print(describe(path))
74+ return 0
75+
76+
6677 def cmd_example_config(args) -> int:
6778 cfg = default_config(args.exciter)
6879 if args.path:
@@ -78,8 +89,13 @@ def main(argv=None) -> int:
7889 ap = argparse.ArgumentParser(prog="cymsim", description="Simulate a struck cymbal.")
7990 sub = ap.add_subparsers(dest="command", required=True)
8091
81- p = sub.add_parser("simulate", help="run a simulation from a YAML config")
92+ p = sub.add_parser("simulate", aliases=["run"], help="run a simulation from a YAML config")
8293 p.add_argument("config")
94+ p.add_argument("-o", "--output", metavar="FILE.cymsim", help="write the run to a .cymsim file (docs/format.md)")
95+ p.add_argument("--save-every", type=float, metavar="SECONDS",
96+ help="with -o, also save the partial run every SECONDS of wall time")
97+ p.add_argument("--record-acceleration", action="store_true",
98+ help="keep the normal acceleration of every node at the audio rate (output.record_acceleration)")
8399 p.add_argument("--duration", type=float, help="override numerics.duration [s]")
84100 p.add_argument("--out-prefix", help="override output.prefix")
85101 p.add_argument("--linear", action="store_true", help="drop the quadratic strain terms (reference run)")
@@ -91,6 +107,11 @@ def main(argv=None) -> int:
91107 p.add_argument("--quiet", action="store_true")
92108 p.set_defaults(func=cmd_check)
93109
110+ p = sub.add_parser("info", help="summarise .cymsim files")
111+ p.add_argument("files", nargs="+")
112+ p.add_argument("--quiet", action="store_true")
113+ p.set_defaults(func=cmd_info)
114+
94115 p = sub.add_parser("example-config", help="print or write a default config")
95116 p.add_argument("path", nargs="?")
96117 p.add_argument("--exciter", choices=("mallet", "stick"), default="mallet")
src/cymsim/config.pymodified+4−0View file
@@ -135,6 +135,8 @@ class Output:
135135 wav: bool = True
136136 npz: bool = True
137137 cascade_plot: bool = True # needs matplotlib
138+ record_acceleration: bool = False # keep the normal acceleration of every node at the audio rate
139+ displacement_rate: float = 60.0 # frames per second of transverse displacement kept for display; 0 = none
138140
139141
140142 # ----------------------------------------------------------------------------- simulation
@@ -160,6 +162,8 @@ class SimulationConfig:
160162 raise ValueError("numerics.cfl_safety must lie in (0, 1]")
161163 if self.mesh.element_size <= 0 or self.mesh.element_size > p.radius / 2:
162164 raise ValueError("mesh.element_size must be positive and well below the radius")
165+ if self.output.displacement_rate < 0:
166+ raise ValueError("output.displacement_rate must be non-negative")
163167 for i, h in enumerate(self.hits):
164168 if h.time < 0:
165169 raise ValueError(f"hit {i}: time must be non-negative")
src/cymsim/cymsim_file.pyadded+253−0View file
@@ -0,0 +1,253 @@
1+"""Reading and writing .cymsim files (HDF5). The layout is specified in docs/format.md.
2+
3+A .cymsim file holds the configuration as YAML text, the mesh, the listener geometry, the
4+optional record of the normal acceleration of every node at the audio rate, the transverse
5+displacement at a low frame rate, the two listener signals, the energy log, the hit summaries
6+and, for a run written by this package, the integrator state needed to continue it. The time
7+series are resizable along time, so a file can be written while the run is in progress and
8+extended afterwards; a partial file is a consistent prefix of the complete one.
9+
10+``save_cymsim`` writes a finished (or partial) SimulationResult in one go; ``CymsimWriter`` is
11+the streaming form used by ``Simulator.run(save_path=...)``; ``load_cymsim`` reads a file back
12+into a SimulationResult; ``describe`` summarises one without loading the large datasets.
13+"""
14+from __future__ import annotations
15+
16+import datetime as _dt
17+from pathlib import Path
18+
19+import h5py
20+import numpy as np
21+import yaml
22+
23+from .config import SimulationConfig
24+from .decimation import CUTOFF, DECIMATION_SPAN, KAISER_BETA
25+from .mesh import Mesh
26+from .simulate import (AccelerationRecord, DisplacementRecord, EnergyRecord, HitSummary,
27+ SimulationResult, SimulatorState)
28+
29+FORMAT_VERSION = 1
30+ENERGY_COLUMNS = "time kinetic membrane bending exciters max_w_over_h"
31+HIT_FIELDS = ("time", "peak_force", "contact_time", "rebound_velocity", "energy_in", "energy_out")
32+STATE_WRITER = "cymsim-python"
33+
34+
35+class CymsimWriter:
36+ """Create a .cymsim file for a run and extend it as the run proceeds.
37+
38+ ``update(result)`` appends whatever time-series data the result holds beyond what is already
39+ in the file, and overwrites the small state-dependent parts (hits, energy log, integrator
40+ state, progress attributes). The static parts (configuration, mesh, listener, kernel) are
41+ written once at construction.
42+ """
43+
44+ def __init__(self, path: str | Path, result: SimulationResult):
45+ self.path = str(path)
46+ r = result
47+ N, E = r.n_nodes, r.n_elements
48+ f = self.f = h5py.File(self.path, "w")
49+ f.attrs["format"] = "cymsim"
50+ f.attrs["format_version"] = FORMAT_VERSION
51+ f.attrs["creator"] = r.creator or "cymsim"
52+ f.attrs["created"] = _dt.datetime.now(_dt.timezone.utc).isoformat(timespec="seconds")
53+
54+ f.create_dataset("config/yaml", data=yaml.safe_dump(r.config.to_dict(), sort_keys=False),
55+ dtype=h5py.string_dtype())
56+
57+ g = f.create_group("mesh")
58+ g.attrs["n_nodes"], g.attrs["n_elements"] = N, E
59+ g.create_dataset("xy", data=np.asarray(r.mesh.xy, np.float64))
60+ g.create_dataset("z0", data=np.asarray(r.z0_node, np.float64))
61+ g.create_dataset("thickness", data=np.asarray(r.h_node, np.float64))
62+ g.create_dataset("normal", data=np.asarray(r.normal_node, np.float64))
63+ g.create_dataset("triangles", data=np.asarray(r.mesh.tri, np.int32))
64+ g.create_dataset("area", data=np.asarray(r.mesh.area, np.float64))
65+ g.create_dataset("boundary", data=np.asarray(r.mesh.boundary, np.int32))
66+
67+ g = f.create_group("listener")
68+ L = r.config.listener
69+ g.attrs["position"] = np.asarray(L.position, np.float64)
70+ g.attrs["c_air"], g.attrs["rho_air"] = L.c_air, L.rho_air
71+ g.attrs["back_diffraction_gain"] = L.back_diffraction_gain
72+ g.create_dataset("r_front", data=np.asarray(r.r_front, np.float64))
73+ g.create_dataset("r_back", data=np.asarray(r.r_back, np.float64))
74+
75+ g = f.create_group("run")
76+ g.attrs["fs"], g.attrs["dt"], g.attrs["substeps"] = int(r.fs), float(r.dt), int(r.substeps)
77+ g.attrs["duration"] = float(r.config.numerics.duration)
78+ g.attrs["n_samples"] = int(round(r.fs * r.config.numerics.duration))
79+ g.attrs["decimation_span"] = DECIMATION_SPAN
80+ g.attrs["n_steps_done"], g.attrs["complete"] = 0, np.int32(0)
81+
82+ g = f.create_group("audio")
83+ g.attrs["fs"], g.attrs["units"] = int(r.fs), "Pa"
84+ for name in ("baffled", "unbaffled"):
85+ g.create_dataset(name, shape=(0,), maxshape=(None,), dtype=np.float32, chunks=(65536,))
86+
87+ if r.acceleration is not None:
88+ a = r.acceleration
89+ g = f.create_group("acceleration")
90+ g.attrs["fs"] = int(r.fs)
91+ g.attrs["component"] = "normal"
92+ g.attrs["units"] = "m/s^2"
93+ g.attrs["encoding"] = "int16 * scale[t]"
94+ g.attrs["description"] = ("acceleration of every node along the rest-surface normal (front, +z side, "
95+ "positive), low-pass filtered by the kernel in 'taps' and sampled at fs; "
96+ "sample j is centred on t = j / fs")
97+ # time-major chunks of 1024 samples (fewer for very large meshes, keeping a chunk under 16 MB)
98+ rows = int(max(1, min(1024, 16_000_000 // (2 * N))))
99+ g.create_dataset("normal", shape=(0, N), maxshape=(None, N), dtype=np.int16, chunks=(rows, N))
100+ g.create_dataset("scale", shape=(0,), maxshape=(None,), dtype=np.float32, chunks=(65536,))
101+ d = g.create_dataset("taps", data=np.asarray(a.taps, np.float64))
102+ d.attrs["layout"] = "taps[w, kk]: weight of sub-step j0 m + kk in audio sample j0 + w - W/2 + 1"
103+ d.attrs["kaiser_beta"], d.attrs["cutoff"] = KAISER_BETA, CUTOFF
104+
105+ if r.displacement is not None:
106+ g = f.create_group("displacement")
107+ g.attrs["rate"] = float(r.displacement.rate)
108+ g.attrs["units"] = "m"
109+ g.attrs["description"] = "transverse displacement w of every node at the frame times"
110+ g.create_dataset("w", shape=(0, N), maxshape=(None, N), dtype=np.float32, chunks=(64, N))
111+ g.create_dataset("time", shape=(0,), maxshape=(None,), dtype=np.float64, chunks=(4096,))
112+
113+ g = f.create_group("energy")
114+ d = g.create_dataset("table", shape=(0, 6), maxshape=(None, 6), dtype=np.float64, chunks=(1024, 6))
115+ d.attrs["columns"] = ENERGY_COLUMNS
116+
117+ g = f.create_group("hits")
118+ K = len(r.hits)
119+ for name in HIT_FIELDS:
120+ g.create_dataset(name, shape=(K,), dtype=np.float64)
121+ g.create_dataset("kind", shape=(K,), dtype="S16")
122+
123+ @staticmethod
124+ def _append(ds, data):
125+ data = np.asarray(data)
126+ n0 = ds.shape[0]
127+ if len(data) <= n0:
128+ return
129+ ds.resize(len(data), axis=0)
130+ ds[n0:] = data[n0:]
131+
132+ def update(self, result: SimulationResult) -> None:
133+ f, r = self.f, result
134+ self._append(f["audio/baffled"], r.pressure_rayleigh.astype(np.float32))
135+ self._append(f["audio/unbaffled"], r.pressure_unbaffled.astype(np.float32))
136+ if r.acceleration is not None and "acceleration" in f:
137+ self._append(f["acceleration/normal"], r.acceleration.data)
138+ self._append(f["acceleration/scale"], r.acceleration.scale)
139+ if r.displacement is not None and "displacement" in f:
140+ self._append(f["displacement/w"], r.displacement.w)
141+ self._append(f["displacement/time"], r.displacement.time)
142+ table = np.array([[e.time, e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h]
143+ for e in r.energy_log], np.float64).reshape(-1, 6)
144+ f["energy/table"].resize(len(table), axis=0)
145+ if len(table):
146+ f["energy/table"][...] = table
147+ for name in HIT_FIELDS:
148+ f["hits"][name][...] = np.array([getattr(h, name) for h in r.hits], np.float64)
149+ f["hits/kind"][...] = np.array([h.kind.encode() for h in r.hits], "S16")
150+ if "state" in f:
151+ del f["state"]
152+ if r.state is not None:
153+ s = r.state
154+ g = f.create_group("state")
155+ g.attrs["writer"] = STATE_WRITER
156+ g.attrs["n_steps_done"] = int(s.n_steps_done)
157+ g.attrs["dtype"] = r.config.numerics.dtype
158+ for name in ("q", "v", "exciter_z", "exciter_v", "peak_force", "contact_steps", "ring_front",
159+ "ring_back", "decimator_tail_front", "decimator_tail_back"):
160+ g.create_dataset(name, data=np.asarray(getattr(s, name)))
161+ if s.acceleration_ring is not None:
162+ g.create_dataset("acceleration_ring", data=np.asarray(s.acceleration_ring))
163+ f["run"].attrs["n_steps_done"] = int(r.n_steps_done)
164+ f["run"].attrs["complete"] = np.int32(1 if r.complete else 0) # an integer flag, readable by every HDF5 client
165+ f["run"].attrs["wall_time"] = float(r.wall_time)
166+ f.flush()
167+
168+ def close(self) -> None:
169+ self.f.close()
170+
171+
172+def save_cymsim(result: SimulationResult, path: str | Path) -> str:
173+ """Write ``result`` (complete or partial) to a .cymsim file and return the path."""
174+ w = CymsimWriter(path, result)
175+ w.update(result)
176+ w.close()
177+ return str(path)
178+
179+
180+def _read_str(ds) -> str:
181+ v = ds[()]
182+ return v.decode() if isinstance(v, bytes) else str(v)
183+
184+
185+def load_cymsim(path: str | Path) -> SimulationResult:
186+ """Read a .cymsim file into a SimulationResult. Files written by other programs (the web
187+ companion) load too; their ``state`` group, if any, is ignored."""
188+ with h5py.File(str(path), "r") as f:
189+ if f.attrs.get("format", "") not in ("cymsim", b"cymsim"):
190+ raise ValueError(f"{path} is not a .cymsim file")
191+ config = SimulationConfig.from_dict(yaml.safe_load(_read_str(f["config/yaml"])) or {})
192+ mg = f["mesh"]
193+ mesh = Mesh(xy=mg["xy"][...], tri=mg["triangles"][...].astype(np.int64), boundary=mg["boundary"][...].astype(np.int64))
194+ run = f["run"].attrs
195+ fs, dt, substeps = int(run["fs"]), float(run["dt"]), int(run["substeps"])
196+ acc = None
197+ if "acceleration" in f:
198+ g = f["acceleration"]
199+ acc = AccelerationRecord(fs, g["normal"][...], g["scale"][...], g["taps"][...])
200+ disp = None
201+ if "displacement" in f:
202+ g = f["displacement"]
203+ disp = DisplacementRecord(float(g.attrs["rate"]), g["time"][...], g["w"][...])
204+ table = f["energy/table"][...]
205+ energy_log = [EnergyRecord(*row) for row in table.tolist()]
206+ hg = f["hits"]
207+ kinds = [k.decode() if isinstance(k, bytes) else str(k) for k in hg["kind"][...]]
208+ hits = [HitSummary(float(hg["time"][k]), kinds[k], float(hg["peak_force"][k]), float(hg["contact_time"][k]),
209+ float(hg["rebound_velocity"][k]), float(hg["energy_in"][k]), float(hg["energy_out"][k]))
210+ for k in range(len(kinds))]
211+ state = None
212+ if "state" in f and f["state"].attrs.get("writer", "") == STATE_WRITER:
213+ g = f["state"]
214+ state = SimulatorState(
215+ int(g.attrs["n_steps_done"]), g["q"][...], g["v"][...], g["exciter_z"][...], g["exciter_v"][...],
216+ g["peak_force"][...], g["contact_steps"][...], g["ring_front"][...], g["ring_back"][...],
217+ g["decimator_tail_front"][...], g["decimator_tail_back"][...],
218+ g["acceleration_ring"][...] if "acceleration_ring" in g else None)
219+ lg = f["listener"]
220+ creator = f.attrs.get("creator", "")
221+ return SimulationResult(
222+ config, fs, dt, substeps, mesh, mg["z0"][...], mg["thickness"][...], mg["normal"][...],
223+ lg["r_front"][...], lg["r_back"][...], f["audio/baffled"][...].astype(np.float64),
224+ f["audio/unbaffled"][...].astype(np.float64), hits, energy_log, acc, disp,
225+ bool(run.get("complete", False)), int(run.get("n_steps_done", 0)), state,
226+ float(run.get("wall_time", 0.0)), creator.decode() if isinstance(creator, bytes) else str(creator))
227+
228+
229+def describe(path: str | Path) -> str:
230+ """A short human-readable summary of a .cymsim file."""
231+ with h5py.File(str(path), "r") as f:
232+ run, a = f["run"].attrs, f.attrs
233+ N, E = int(f["mesh"].attrs["n_nodes"]), int(f["mesh"].attrs["n_elements"])
234+ fs, dt = int(run["fs"]), float(run["dt"])
235+ done = int(run.get("n_steps_done", 0)) * dt
236+ lines = [f"{path}: cymsim format v{int(a['format_version'])}, written by {a.get('creator', '?')} on {a.get('created', '?')}",
237+ f" {'complete' if run.get('complete', False) else 'PARTIAL'}: {done:.3f} of {float(run['duration']):.3f} s simulated",
238+ f" mesh: {N} nodes, {E} elements; fs = {fs} Hz, {int(run['substeps'])} sub-steps per sample",
239+ f" audio: {f['audio/baffled'].shape[0]} samples ({f['audio/baffled'].shape[0] / fs:.3f} s), baffled and unbaffled"]
240+ if "acceleration" in f:
241+ d = f["acceleration/normal"]
242+ lines.append(f" acceleration (normal, int16): {d.shape[0]} samples x {d.shape[1]} nodes, "
243+ f"{d.shape[0] * d.shape[1] * 2 / 1e6:.1f} MB")
244+ else:
245+ lines.append(" acceleration: not recorded")
246+ if "displacement" in f:
247+ d = f["displacement/w"]
248+ lines.append(f" displacement: {d.shape[0]} frames at {float(f['displacement'].attrs['rate']):g} Hz")
249+ lines.append(f" energy log: {f['energy/table'].shape[0]} entries; hits: {f['hits/time'].shape[0]}")
250+ if "state" in f:
251+ lines.append(f" integrator state present (writer: {f['state'].attrs.get('writer', '?')}, "
252+ f"step {int(f['state'].attrs.get('n_steps_done', 0))})")
253+ return "\n".join(lines)
src/cymsim/decimation.pyadded+88−0View file
@@ -0,0 +1,88 @@
1+"""The anti-aliasing decimation kernel shared by the acceleration record and the listener
2+pressure, and a streaming host-side decimator.
3+
4+The time loop runs at m sub-steps per audio sample. Every quantity that is kept at the audio
5+rate (the normal acceleration of each node, the retarded-time pressure sums) is low-pass
6+filtered by the same finite-impulse-response kernel h(d), d being the offset in sub-steps from
7+the audio sample's centre, and read out once per audio sample. The kernel is a Kaiser-windowed
8+sinc spanning DECIMATION_SPAN audio samples, linear-phase, with its delay removed so that stored
9+sample j is centred exactly on t = j / fs. It is stored in every .cymsim file so that a reader
10+knows the exact frequency response of the record. See docs/format.md and docs/model.md, section
11+"Radiation".
12+
13+``decimation_taps(m)`` returns the kernel in the polyphase layout ``taps[w, kk]`` used by the
14+in-loop accumulators: sub-step s = j0 m + kk contributes ``taps[w, kk]`` times its value to audio
15+sample j = j0 + w - W/2 + 1 (slot w of W). The web companion uses the same layout.
16+"""
17+from __future__ import annotations
18+
19+import numpy as np
20+
21+DECIMATION_SPAN = 32 # audio samples spanned by the kernel (W)
22+KAISER_BETA = 8.0 # about 80 dB of stopband attenuation
23+CUTOFF = 0.46 # -6 dB point in units of fs; passband to about 0.38 fs, stopband from 0.54 fs
24+
25+
26+def decimation_kernel(m: int, span: int = DECIMATION_SPAN, beta: float = KAISER_BETA,
27+ cutoff: float = CUTOFF) -> np.ndarray:
28+ """h(d) for d = -W m / 2 ... W m / 2 - 1, normalised to unit DC gain. Shape (W m,)."""
29+ half = span * m / 2
30+ d = np.arange(-span * m // 2, span * m // 2, dtype=np.float64)
31+ fc = cutoff / m # cycles per sub-step
32+ x = 2 * np.pi * fc * d
33+ sinc = np.where(d == 0, 1.0, np.sin(x) / np.where(d == 0, 1.0, x))
34+ u = d / half
35+ win = np.i0(beta * np.sqrt(np.clip(1 - u * u, 0, None))) / np.i0(beta)
36+ h = sinc * win
37+ return h / h.sum()
38+
39+
40+def decimation_taps(m: int, span: int = DECIMATION_SPAN, beta: float = KAISER_BETA,
41+ cutoff: float = CUTOFF) -> np.ndarray:
42+ """The kernel in polyphase layout (W, m): taps[w, kk] = h(kk + (W/2 - 1 - w) m)."""
43+ h = decimation_kernel(m, span, beta, cutoff)
44+ w = np.arange(span)[:, None]
45+ kk = np.arange(m)[None, :]
46+ d = kk + (span // 2 - 1 - w) * m
47+ return h[d + span * m // 2]
48+
49+
50+class StreamingDecimator:
51+ """Decimate a sub-step-rate stream to the audio rate with ``decimation_kernel``.
52+
53+ ``push(x)`` accepts the next sub-steps (contiguous, in order) and returns every audio sample
54+ that has become final, i.e. sample j once sub-steps up to j m + W m / 2 - 1 are known. The
55+ pre-roll samples with j < 0 are dropped. ``tail`` holds the sub-steps still needed by
56+ samples that are not yet final; it is part of the resumable state.
57+ """
58+
59+ def __init__(self, m: int, span: int = DECIMATION_SPAN, beta: float = KAISER_BETA,
60+ cutoff: float = CUTOFF, tail: np.ndarray | None = None, n_seen: int = 0):
61+ self.m, self.span = m, span
62+ self.h = decimation_kernel(m, span, beta, cutoff)
63+ self.half = span * m // 2
64+ # the sub-steps [n_seen - len(tail), n_seen) are held; the next final sample is j_next.
65+ # A fresh stream starts with W m / 2 zeros: the shell is at rest before t = 0.
66+ self.tail = np.zeros(self.half) if tail is None else np.asarray(tail, dtype=np.float64)
67+ self.n_seen = int(n_seen)
68+ self.j_next = max(0, -((-(self.n_seen - len(self.tail) + self.half)) // m)) # ceil
69+
70+ def push(self, x: np.ndarray) -> np.ndarray:
71+ x = np.asarray(x, dtype=np.float64)
72+ buf = np.concatenate([self.tail, x])
73+ g0 = self.n_seen - len(self.tail) # global index of buf[0]
74+ self.n_seen += len(x)
75+ L = self.span * self.m
76+ # sample j needs buf indices [j m - half - g0, j m + half - g0)
77+ j_last = (self.n_seen - self.half) // self.m
78+ n_j = j_last - self.j_next + 1
79+ if n_j <= 0:
80+ self.tail = buf[max(0, self.j_next * self.m - self.half - g0):]
81+ return np.zeros(0)
82+ i0 = self.j_next * self.m - self.half - g0
83+ assert i0 >= 0
84+ win = np.lib.stride_tricks.sliding_window_view(buf, L)[i0::self.m][:n_j]
85+ out = win @ self.h
86+ self.j_next += n_j
87+ self.tail = buf[self.j_next * self.m - self.half - g0:]
88+ return out
src/cymsim/radiation.pymodified+25−16View file
@@ -1,20 +1,24 @@
1-"""Sound pressure at a listener point from the shell's transverse acceleration.
1+"""Sound pressure at a listener point from the shell's normal acceleration.
22
33 The baffled plate is treated with the Rayleigh integral: every element is a point source of
4-volume acceleration A_e w_ddot_e at its centroid on the curved surface, delayed by the travel time
5-to the listener. An unbaffled cymbal also radiates from its back with opposite sign, and that
4+volume acceleration A_e a_e at its centroid on the curved surface, delayed by the travel time to
5+the listener, where a_e is the mean over the element's nodes of the acceleration along the rest
6+surface normal. An unbaffled cymbal also radiates from its back with opposite sign, and that
67 contribution reaches the listener by diffracting around the rim. We approximate it by the same
78 Rayleigh sum with the path length replaced by the shortest route over the rim (a Kirchhoff-type
89 edge-diffraction estimate, not a boundary-element solution) and subtract it from the front term.
10+
911 Retarded times are handled at sub-step resolution with linear interpolation into two delay
10-buffers, which are decimated to the audio rate at the end. See docs/model.md, section
11-"Radiation".
12+buffers. The buffers are rings that hold a little more than two blocks plus the longest delay;
13+after each block the host reads the finished sub-steps, clears their slots and decimates them to
14+the audio rate with the kernel of ``decimation.py``, the same kernel that filters the recorded
15+node accelerations, so that the stored audio follows from the stored accelerations. See
16+docs/model.md, section "Radiation".
1217 """
1318 from __future__ import annotations
1419
1520 import jax.numpy as jnp
1621 import numpy as np
17-from scipy.signal import resample_poly
1822
1923 from .config import Listener
2024 from .geometry import CymbalProfile
@@ -23,7 +27,7 @@ from .mesh import Mesh
2327
2428 class RadiationModel:
2529 def __init__(self, mesh: Mesh, profile: CymbalProfile, z0_node: np.ndarray, listener: Listener,
26- dt: float, n_steps: int, dtype, n_rim: int = 720):
30+ dt: float, steps_per_block: int, dtype, n_rim: int = 720):
2731 L = np.asarray(listener.position, dtype=float)
2832 self.listener = listener
2933 cent = np.column_stack([mesh.centroid, z0_node[mesh.tri].mean(1)]) # (E, 3) on the surface
@@ -41,9 +45,13 @@ class RadiationModel:
4145 c, rho = listener.c_air, listener.rho_air
4246 coef_f = rho * mesh.area / (2 * np.pi * r_front)
4347 coef_b = listener.back_diffraction_gain * rho * mesh.area / (2 * np.pi * r_back)
48+ self.coef_front_np, self.coef_back_np = coef_f, coef_b
4449 df, db = r_front / (c * dt), r_back / (c * dt)
4550 i0f, i0b = np.floor(df).astype(np.int32), np.floor(db).astype(np.int32)
46- self.buffer_length = int(n_steps + max(i0f.max(), i0b.max()) + 2)
51+ self.max_delay = int(max(i0f.max(), i0b.max())) + 2
52+ # slots of sub-steps [n0, n0 + spb) are read and cleared by the host after the block that
53+ # ends at n0 + spb; the next block writes up to n0 + 2 spb + max_delay, which must not wrap
54+ self.ring_length = int(2 * steps_per_block + self.max_delay + 1)
4755 f = lambda a: jnp.asarray(a, dtype=dtype)
4856 self.coef_f, self.coef_b = f(coef_f), f(coef_b)
4957 self.frac_f, self.frac_b = f(df - np.floor(df)), f(db - np.floor(db))
@@ -51,17 +59,18 @@ class RadiationModel:
5159 self.dtype = dtype
5260
5361 def initial_buffers(self):
54- z = jnp.zeros(self.buffer_length, self.dtype)
62+ z = jnp.zeros(self.ring_length, self.dtype)
5563 return z, z
5664
5765 def scatter(self, n, acc_e, bf, bb):
58- """Add the element accelerations of step n into the retarded-time buffers."""
66+ """Add the element accelerations of step n into the retarded-time rings."""
67+ L = self.ring_length
5968 sf, sb = self.coef_f * acc_e, self.coef_b * acc_e
60- bf = bf.at[n + self.i0f].add(sf * (1 - self.frac_f)).at[n + self.i0f + 1].add(sf * self.frac_f)
61- bb = bb.at[n + self.i0b].add(sb * (1 - self.frac_b)).at[n + self.i0b + 1].add(sb * self.frac_b)
69+ jf, jb = (n + self.i0f) % L, (n + self.i0b) % L
70+ bf = bf.at[jf].add(sf * (1 - self.frac_f)).at[(jf + 1) % L].add(sf * self.frac_f)
71+ bb = bb.at[jb].add(sb * (1 - self.frac_b)).at[(jb + 1) % L].add(sb * self.frac_b)
6272 return bf, bb
6373
64- @staticmethod
65- def to_audio(buffer: np.ndarray, n_steps: int, substeps: int, n_out: int) -> np.ndarray:
66- """Decimate a sub-step-rate buffer to the audio rate (anti-aliased polyphase resampling)."""
67- return resample_poly(np.asarray(buffer[:n_steps], dtype=np.float64), 1, substeps)[:n_out]
74+ def slots(self, n0: int, n1: int) -> np.ndarray:
75+ """Ring slots of sub-steps [n0, n1), all of which are final once step n1 - 1 is done."""
76+ return np.arange(n0, n1) % self.ring_length
src/cymsim/shell.pymodified+14−0View file
@@ -62,6 +62,19 @@ class ShellModel:
6262 self.z0x = (mesh.Gx * z0_e).sum(1)
6363 self.z0y = (mesh.Gy * z0_e).sum(1)
6464 self.max_slope = float(np.hypot(self.z0x, self.z0y).max())
65+ # ---- unit normal of the rest surface at the nodes, from the area-weighted mean of the
66+ # CST slopes of the adjacent elements; +z side is the front (bell up). The normal
67+ # acceleration n . (u_dd, v_dd, w_dd) is what is recorded and what radiates.
68+ gx, gy = np.zeros(N), np.zeros(N)
69+ for k in range(3):
70+ np.add.at(gx, tri[:, k], area * self.z0x)
71+ np.add.at(gy, tri[:, k], area * self.z0y)
72+ a_node = np.zeros(N)
73+ for k in range(3):
74+ np.add.at(a_node, tri[:, k], area)
75+ gx, gy = gx / a_node, gy / a_node
76+ nrm = np.sqrt(1 + gx ** 2 + gy ** 2)
77+ self.normal_node = np.column_stack([-gx / nrm, -gy / nrm, 1 / nrm])
6578
6679 # ---- element stiffness: bending rigidity D, membrane stiffness C, DKT matrix
6780 E_, nu, rho = material.E, material.nu, material.rho
@@ -124,6 +137,7 @@ class ShellModel:
124137 self.z0x_j, self.z0y_j, self.Cm_j, self.Ke_j = f(self.z0x), f(self.z0y), f(self.Cm_e), f(self.Ke)
125138 self.rhoh_area_j = f(rho * self.h_e * area)
126139 self.h_node_j = f(self.h_node)
140+ self.normal_j = f(self.normal_node)
127141 self.M_j, self.Minv_j, self.Rm_j = f(M), f(1.0 / M), f(Rm)
128142 self._sigma1 = float(material.sigma1)
129143 self._N = N
src/cymsim/simulate.pymodified+288−44View file
@@ -2,15 +2,19 @@
22 and radiation, block-compiled with JAX.
33
44 One sub-step advances the state (q, v) of the shell, the state (zs, vs) of every exciter and the
5-two retarded-time pressure buffers:
5+two retarded-time pressure rings:
66
77 f = f_int(q) - contact(q, zs) - 2 sigma1 rho h Lap(v_w)
88 v' = P_rigid [ (v + dt M^-1 f) / (1 + 2 sigma0 dt) ]
99 q' = q + dt v'
1010 vs' = vs + dt F / m, zs' = zs + dt vs'
1111
12-where P_rigid removes the six rigid-body modes. The sub-step rate is an integer multiple of the
13-audio rate chosen from the linear stability bound. See docs/model.md, section "Time integration".
12+where P_rigid removes the six rigid-body modes. The sub-step rate is an integer multiple m of the
13+audio rate chosen from the linear stability bound. The acceleration that radiates, and that the
14+optional record keeps, is the normal component of (v' - v) / dt, which for this scheme is the
15+exact central second difference of q. It is decimated to the audio rate inside the loop by the
16+polyphase accumulators of ``decimation.py``, so that no second difference is ever taken on
17+stored displacements. See docs/model.md, section "Time integration", and docs/format.md.
1418 """
1519 from __future__ import annotations
1620
@@ -26,9 +30,10 @@ from jax import lax
2630 from scipy.io import wavfile
2731
2832 from .config import SimulationConfig
33+from .decimation import DECIMATION_SPAN, StreamingDecimator, decimation_taps
2934 from .exciter import HitSchedule, build_patch
3035 from .geometry import CymbalProfile
31-from .mesh import disk_mesh
36+from .mesh import Mesh, disk_mesh
3237 from .radiation import RadiationModel
3338 from .shell import N_DOF, ShellModel
3439
@@ -60,31 +65,136 @@ class EnergyRecord:
6065 return self.kinetic + self.membrane + self.bending + self.exciters
6166
6267
68+@dataclass
69+class AccelerationRecord:
70+ """Normal acceleration of every node at the audio rate, int16 against a scale per sample.
71+
72+ ``data[j, i]`` * ``scale[j]`` is the acceleration of node i along the rest-surface normal
73+ (front side positive) at t = j / fs, low-pass filtered by the kernel ``taps`` (see
74+ decimation.py). Sample j is final once the run has passed t = (j + W/2) / fs.
75+ """
76+ fs: int
77+ data: np.ndarray # (T, N) int16
78+ scale: np.ndarray # (T,) float32
79+ taps: np.ndarray # (W, m) float64 polyphase kernel
80+
81+ @property
82+ def n_samples(self) -> int:
83+ return len(self.data)
84+
85+ @property
86+ def time(self) -> np.ndarray:
87+ return np.arange(self.n_samples) / self.fs
88+
89+ def as_float(self, start: int = 0, stop: int | None = None) -> np.ndarray:
90+ """Decoded accelerations [m/s^2] for samples [start, stop), float32 (T', N)."""
91+ return self.data[start:stop].astype(np.float32) * self.scale[start:stop, None]
92+
93+ @staticmethod
94+ def encode(a: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
95+ """(int16 data, float32 scale) of a float block (T', N)."""
96+ peak = np.abs(a).max(axis=1)
97+ scale = np.where(peak > 0, peak / 32767.0, 1.0).astype(np.float32)
98+ return np.round(a / scale[:, None]).astype(np.int16), np.where(peak > 0, scale, 0.0).astype(np.float32)
99+
100+ @classmethod
101+ def empty(cls, fs: int, n_nodes: int, taps: np.ndarray) -> "AccelerationRecord":
102+ return cls(fs, np.zeros((0, n_nodes), np.int16), np.zeros(0, np.float32), taps)
103+
104+ def extend(self, data: np.ndarray, scale: np.ndarray) -> None:
105+ self.data = np.concatenate([self.data, data])
106+ self.scale = np.concatenate([self.scale, scale])
107+
108+
109+@dataclass
110+class DisplacementRecord:
111+ """Transverse displacement w of every node [m] at a low frame rate, for display."""
112+ rate: float
113+ time: np.ndarray # (Td,)
114+ w: np.ndarray # (Td, N) float32
115+
116+ @property
117+ def n_frames(self) -> int:
118+ return len(self.time)
119+
120+ @classmethod
121+ def empty(cls, rate: float, n_nodes: int) -> "DisplacementRecord":
122+ return cls(rate, np.zeros(0), np.zeros((0, n_nodes), np.float32))
123+
124+ def extend(self, time: np.ndarray, w: np.ndarray) -> None:
125+ self.time = np.concatenate([self.time, time])
126+ self.w = np.concatenate([self.w, w])
127+
128+
129+@dataclass
130+class SimulatorState:
131+ """Everything the Python integrator needs to continue a run from step ``n_steps_done``.
132+
133+ Written to the ``state`` group of a .cymsim file. It is specific to this implementation
134+ (the web companion keeps its own); readers of the format may ignore it.
135+ """
136+ n_steps_done: int
137+ q: np.ndarray # (N, 5)
138+ v: np.ndarray # (N, 5)
139+ exciter_z: np.ndarray # (K,)
140+ exciter_v: np.ndarray # (K,)
141+ peak_force: np.ndarray # (K,)
142+ contact_steps: np.ndarray # (K,) int32
143+ ring_front: np.ndarray # (L,) pending retarded-time pressure, sub-step rate
144+ ring_back: np.ndarray # (L,)
145+ decimator_tail_front: np.ndarray # sub-steps still needed by non-final audio samples
146+ decimator_tail_back: np.ndarray
147+ acceleration_ring: np.ndarray | None = None # (W, N) partial audio samples of the record
148+
149+
63150 @dataclass
64151 class SimulationResult:
65152 config: SimulationConfig
66153 fs: int
67154 dt: float
68155 substeps: int
69- n_nodes: int
70- n_elements: int
156+ mesh: Mesh
157+ z0_node: np.ndarray # (N,) rest height of the nodes [m]
158+ h_node: np.ndarray # (N,) thickness at the nodes [m]
159+ normal_node: np.ndarray # (N, 3) unit normal of the rest surface, front side
160+ r_front: np.ndarray # (E,) element centroid to listener [m]
161+ r_back: np.ndarray # (E,) shortest path over the rim to the listener [m]
71162 pressure_rayleigh: np.ndarray # baffled Rayleigh integral at the listener [Pa]
72163 pressure_unbaffled: np.ndarray # front minus diffracted back term [Pa]
73164 hits: list[HitSummary]
74165 energy_log: list[EnergyRecord] = field(default_factory=list)
166+ acceleration: AccelerationRecord | None = None
167+ displacement: DisplacementRecord | None = None
168+ complete: bool = True
169+ n_steps_done: int = 0
170+ state: SimulatorState | None = None
75171 wall_time: float = 0.0
172+ creator: str = ""
173+
174+ @property
175+ def n_nodes(self) -> int:
176+ return self.mesh.n_nodes
177+
178+ @property
179+ def n_elements(self) -> int:
180+ return self.mesh.n_elements
76181
77182 @property
78183 def time(self) -> np.ndarray:
79184 return np.arange(len(self.pressure_rayleigh)) / self.fs
80185
186+ @property
187+ def duration_done(self) -> float:
188+ """Simulated time covered by the run so far [s]."""
189+ return self.n_steps_done * self.dt
190+
81191 def write_wav(self, prefix: str, normalize: bool = True) -> list[str]:
82192 """Write <prefix>_rayleigh.wav and <prefix>_unbaffled.wav as 16-bit PCM. With
83193 ``normalize`` the peak is placed at 0.9 full scale; otherwise 1 Pa maps to 0.05 FS."""
84194 paths = []
85195 for name, p in (("rayleigh", self.pressure_rayleigh), ("unbaffled", self.pressure_unbaffled)):
86- peak = float(np.abs(p).max()) or 1.0
87- scale = 0.9 / peak if normalize else 0.05
196+ peak = float(np.abs(p).max()) if len(p) else 0.0
197+ scale = 0.9 / (peak or 1.0) if normalize else 0.05
88198 path = f"{prefix}_{name}.wav"
89199 wavfile.write(path, self.fs, (np.clip(p * scale, -1, 1) * 32767).astype(np.int16))
90200 paths.append(path)
@@ -98,6 +208,11 @@ class SimulationResult:
98208 energy=[[e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h] for e in self.energy_log])
99209 return path
100210
211+ def save(self, path: str) -> str:
212+ """Write a .cymsim file (see docs/format.md)."""
213+ from .cymsim_file import save_cymsim
214+ return save_cymsim(self, path)
215+
101216
102217 class Simulator:
103218 """Build every part of the model from a SimulationConfig; ``run()`` integrates it."""
@@ -105,7 +220,7 @@ class Simulator:
105220 def __init__(self, config: SimulationConfig):
106221 config.validate()
107222 self.config = config
108- num, mat = config.numerics, config.material
223+ num, mat, out = config.numerics, config.material, config.output
109224 t0 = time.time()
110225
111226 self.profile = CymbalProfile(config.profile, config.imperfections)
@@ -123,17 +238,22 @@ class Simulator:
123238 if sh.max_slope > config.profile.max_slope:
124239 log.warning("max slope %.2f exceeds the shallow-shell bound %.2f", sh.max_slope, config.profile.max_slope)
125240
126- # ---- time step: integer number of sub-steps per audio sample, inside the stability bound
241+ # ---- time step: integer number of sub-steps per audio sample, inside the stability bound.
242+ # The run covers W/2 audio samples beyond the duration so that every sample of the
243+ # duration is final after the decimation kernel's look-ahead.
127244 self.fs = num.fs
128245 self.substeps = math.ceil(1.0 / (self.fs * num.cfl_safety * sh.dt_stable))
129246 self.dt = 1.0 / (self.fs * self.substeps)
130- self.n_blocks = math.ceil(num.duration * self.fs * self.substeps / num.steps_per_block)
247+ self.n_out = int(round(self.fs * num.duration))
248+ self.span = DECIMATION_SPAN
249+ self.n_blocks = math.ceil((self.n_out + self.span // 2) * self.substeps / num.steps_per_block)
131250 self.n_steps = self.n_blocks * num.steps_per_block
132251 xc, yc = self.mesh.centroid[sh.limiting_element]
133252 log.info("omega_max bound %.3e rad/s (limiting element at r = %.3f m, h = %.2f mm); "
134253 "dt = %.3e s, %d sub-steps per sample, %d steps",
135254 sh.omega_max, math.hypot(xc, yc), sh.h_e[sh.limiting_element] * 1e3,
136255 self.dt, self.substeps, self.n_steps)
256+ self.taps = decimation_taps(self.substeps)
137257
138258 # ---- hits
139259 self.patches = [build_patch(h, self.mesh, self.profile) for h in config.hits]
@@ -150,7 +270,24 @@ class Simulator:
150270
151271 # ---- radiation
152272 self.radiation = RadiationModel(self.mesh, self.profile, sh.z0_node, config.listener,
153- self.dt, self.n_steps, sh.dtype)
273+ self.dt, num.steps_per_block, sh.dtype)
274+
275+ # ---- records
276+ self.record = bool(out.record_acceleration)
277+ self.disp_rate = float(out.displacement_rate)
278+ if self.record:
279+ nbytes = self.n_out * (2 * self.mesh.n_nodes + 4)
280+ log.info("acceleration record: %d samples x %d nodes, int16, %.2f GB in memory",
281+ self.n_out, self.mesh.n_nodes, nbytes / 1e9)
282+ if nbytes > 2e9:
283+ log.warning("the acceleration record needs %.1f GB; consider a shorter duration or "
284+ "a coarser mesh", nbytes / 1e9)
285+ S = self.fs * self.substeps # sub-steps per second
286+ if self.disp_rate > 0:
287+ n_frames = int(math.floor(self.n_out * self.substeps * self.disp_rate / S)) + 1
288+ self.frame_steps = np.round(np.arange(n_frames) * S / self.disp_rate).astype(np.int64)
289+ else:
290+ self.frame_steps = np.zeros(0, np.int64)
154291 self._build_step()
155292 self.setup_time = time.time() - t0
156293 log.info("setup done in %.1f s; backend %s", self.setup_time, jax.devices()[0].platform)
@@ -159,10 +296,23 @@ class Simulator:
159296 def _build_step(self):
160297 sh, sched, rad = self.shell, self.schedule, self.radiation
161298 dt, sigma0 = self.dt, self.config.material.sigma0
162- tri = sh.tri_j
163-
164- def step(n, st):
165- q, v, zs, vs, bf, bb, Fmax, ncont = st
299+ m, W, N = self.substeps, self.span, self.mesh.n_nodes
300+ spb = self.config.numerics.steps_per_block
301+ tri, normal = sh.tri_j, sh.normal_j
302+ record, rate = self.record, self.disp_rate
303+ taps_j = jnp.asarray(self.taps, sh.dtype)
304+ S = float(self.fs * m)
305+ self.acc_rows = spb // m + 2 # audio samples emitted per block, at most
306+ self.disp_rows = int(math.ceil(spb * rate / S)) + 2 if rate > 0 else 0
307+
308+ def step(n, k_block0, j_block0, st):
309+ q, v, zs, vs, bf, bb, Fmax, ncont, ring, acc_out, disp_out = st
310+ if rate > 0:
311+ # frame k sits at sub-step round(k S / rate); emit q at step n if n is one of them
312+ k_cand = jnp.round(n * (rate / S))
313+ is_frame = jnp.round(k_cand * (S / rate)) == n
314+ row = jnp.clip(k_cand.astype(jnp.int32) - k_block0, 0, self.disp_rows - 1)
315+ disp_out = disp_out.at[row].add(jnp.where(is_frame, q[:, 2], 0.0))
166316 Fk, F, zs, vs = sched.contact(n, q, zs, vs)
167317 fo = sh.internal_force(q)
168318 fo = sched.apply_to_plate(fo, Fk)
@@ -171,15 +321,28 @@ class Simulator:
171321 v1 = sh.project_rigid(v1)
172322 q1 = q + dt * v1
173323 zs1, vs1 = sched.advance(zs, vs, F, dt)
174- acc_e = ((v1[:, 2] - v[:, 2]) / dt)[tri].mean(1)
324+ # normal acceleration at step n: central second difference of q, exact for leapfrog
325+ a_node = (((v1 - v) / dt)[:, :3] * normal).sum(1)
326+ acc_e = a_node[tri].mean(1)
175327 bf, bb = rad.scatter(n, acc_e, bf, bb)
176- return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0))
177-
178- steps_per_block = self.config.numerics.steps_per_block
328+ if record:
329+ j0, kk = n // m, n % m
330+ emit = kk == 0
331+ j_e = j0 - W // 2 # the sample that just became final
332+ slot_e = j_e % W
333+ done = ring[slot_e]
334+ row = jnp.clip(j_e - j_block0, 0, self.acc_rows - 1)
335+ acc_out = acc_out.at[row].add(jnp.where(emit, done, 0.0))
336+ ring = ring.at[slot_e].set(jnp.where(emit, 0.0, done))
337+ # slot s holds sample j0 + 1 - W/2 + w with s = (j0 + 1 - W/2 + w) mod W
338+ r = (j0 + 1 - W // 2) % W
339+ tw = jnp.take(taps_j[:, kk], (jnp.arange(W) - r) % W)
340+ ring = ring + tw[:, None] * a_node[None, :]
341+ return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0), ring, acc_out, disp_out)
179342
180343 @jax.jit
181- def run_block(n0, st):
182- return lax.fori_loop(0, steps_per_block, lambda k, s: step((n0 + k).astype(jnp.int32), s), st)
344+ def run_block(n0, k_block0, j_block0, st):
345+ return lax.fori_loop(0, spb, lambda k, s: step((n0 + k).astype(jnp.int32), k_block0, j_block0, s), st)
183346
184347 @jax.jit
185348 def energies(st):
@@ -188,54 +351,135 @@ class Simulator:
188351 return (sh.kinetic_energy(v), Em.sum(), Eb.sum(), sched.energy(vs),
189352 jnp.abs(q[:, 2] / sh.h_node_j).max())
190353
191- self._step, self._run_block, self._energies = step, run_block, energies
354+ @jax.jit
355+ def clear(buf, idx):
356+ return buf.at[idx].set(0.0)
357+
358+ self._step, self._run_block, self._energies, self._clear = step, run_block, energies, clear
192359
193360 def initial_state(self):
194361 N, dtype = self.mesh.n_nodes, self.shell.dtype
195362 zs, vs = self.schedule.initial_state()
196363 bf, bb = self.radiation.initial_buffers()
197364 K = self.schedule.n_hits
365+ ring = jnp.zeros((self.span, N) if self.record else (0,), dtype)
198366 return (jnp.zeros((N, N_DOF), dtype), jnp.zeros((N, N_DOF), dtype), zs, vs, bf, bb,
199- jnp.zeros(K, dtype), jnp.zeros(K, jnp.int32))
367+ jnp.zeros(K, dtype), jnp.zeros(K, jnp.int32), ring, self._empty_outputs()[0],
368+ self._empty_outputs()[1])
200369
201- def run(self) -> SimulationResult:
370+ def _empty_outputs(self):
371+ N, dtype = self.mesh.n_nodes, self.shell.dtype
372+ return (jnp.zeros((self.acc_rows, N) if self.record else (0,), dtype),
373+ jnp.zeros((self.disp_rows, N) if self.disp_rate > 0 else (0,), dtype))
374+
375+ def _hit_summaries(self, Fmax, ncont, vs) -> list[HitSummary]:
376+ hits = []
377+ for k, p in enumerate(self.patches):
378+ hits.append(HitSummary(p.time, p.hit.exciter.kind, float(Fmax[k]), float(ncont[k]) * self.dt,
379+ float(vs[k]), p.energy, 0.5 * p.mass * float(vs[k]) ** 2))
380+ return hits
381+
382+ def _result(self, state, energy_log, p_front, p_back, acc, disp, n_done, wall, dec_f, dec_b) -> SimulationResult:
383+ q, v, zs, vs, bf, bb, Fmax, ncont, ring, _, _ = jax.device_get(state)
384+ sim_state = SimulatorState(
385+ n_done, np.asarray(q), np.asarray(v), np.asarray(zs), np.asarray(vs), np.asarray(Fmax),
386+ np.asarray(ncont), np.asarray(bf), np.asarray(bb), dec_f.tail.copy(), dec_b.tail.copy(),
387+ np.asarray(ring) if self.record else None)
388+ n = min(len(p_front), len(p_back))
389+ p_front, p_back = np.asarray(p_front[:n]), np.asarray(p_back[:n])
390+ from . import __version__
391+ return SimulationResult(
392+ self.config, self.fs, self.dt, self.substeps, self.mesh, self.shell.z0_node, self.shell.h_node,
393+ self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front, p_front - p_back,
394+ self._hit_summaries(Fmax, ncont, vs), list(energy_log), acc, disp, n_done >= self.n_steps, n_done,
395+ sim_state, wall, f"cymsim {__version__}")
396+
397+ def run(self, save_path: str | None = None, save_every: float | None = None) -> SimulationResult:
398+ """Integrate the run. With ``save_path`` a .cymsim file is written at the end and, if
399+ ``save_every`` (seconds of wall time) is given, also periodically during the run, each
400+ time as a consistent partial file that can be loaded (see docs/format.md)."""
202401 num = self.config.numerics
402+ m, spb, N = self.substeps, num.steps_per_block, self.mesh.n_nodes
203403 state = self.initial_state()
204404 energy_log: list[EnergyRecord] = []
405+ dec_f, dec_b = StreamingDecimator(m), StreamingDecimator(m)
406+ audio_f, audio_b = [], []
407+ n_audio = 0
408+ acc = AccelerationRecord.empty(self.fs, N, self.taps) if self.record else None
409+ disp = DisplacementRecord.empty(self.disp_rate, N) if self.disp_rate > 0 else None
410+ writer = None
411+ if save_path is not None:
412+ from .cymsim_file import CymsimWriter
413+ writer = CymsimWriter(save_path, self._result(state, [], np.zeros(0), np.zeros(0), acc, disp, 0, 0.0,
414+ dec_f, dec_b))
205415 t_run = time.time()
416+ t_saved = t_run
206417 report_every = max(1, self.n_blocks // max(1, num.n_progress))
207418 log.info("running %d blocks of %d steps (nonlinearity meter: max|w|/h < 1 effectively "
208- "linear; 3 to 5 clear cascade)", self.n_blocks, num.steps_per_block)
419+ "linear; 3 to 5 clear cascade)", self.n_blocks, spb)
209420 for blk in range(self.n_blocks):
210- state = self._run_block(jnp.asarray(blk * num.steps_per_block, jnp.int32), state)
421+ n0, n1 = blk * spb, (blk + 1) * spb
422+ k0 = int(np.searchsorted(self.frame_steps, n0))
423+ k1 = int(np.searchsorted(self.frame_steps, n1))
424+ j_block0 = -(-n0 // m) - self.span // 2
425+ acc_out, disp_out = self._empty_outputs()
426+ state = state[:9] + (acc_out, disp_out)
427+ state = self._run_block(jnp.int32(n0), jnp.int32(k0), jnp.int32(j_block0), state)
428+
429+ # ---- pressure: the ring slots of this block's sub-steps are final; read, clear, decimate
430+ slots = jnp.asarray(self.radiation.slots(n0, n1))
431+ pf, pb = np.asarray(state[4][slots], np.float64), np.asarray(state[5][slots], np.float64)
432+ state = state[:4] + (self._clear(state[4], slots), self._clear(state[5], slots)) + state[6:]
433+ new_f, new_b = dec_f.push(pf), dec_b.push(pb)
434+ keep = max(0, min(len(new_f), self.n_out - n_audio))
435+ audio_f.append(new_f[:keep])
436+ audio_b.append(new_b[:keep])
437+ n_audio += keep
438+
439+ # ---- records
440+ if self.record:
441+ n_emit = -(-n1 // m) - -(-n0 // m)
442+ rows = np.asarray(state[9][:n_emit], np.float32)
443+ j = j_block0 + np.arange(n_emit)
444+ sel = (j >= 0) & (j < self.n_out)
445+ if sel.any():
446+ acc.extend(*AccelerationRecord.encode(rows[sel]))
447+ if self.disp_rate > 0 and k1 > k0:
448+ frames = np.asarray(state[10][:k1 - k0], np.float32)
449+ t_frames = self.frame_steps[k0:k1] * self.dt
450+ sel = self.frame_steps[k0:k1] <= self.n_out * m
451+ disp.extend(t_frames[sel], frames[sel])
452+
211453 if blk % report_every == 0 or blk == self.n_blocks - 1:
212454 ke, em, eb, ks, woh = (float(x) for x in self._energies(state))
213- t = (blk + 1) * num.steps_per_block * self.dt
455+ t = n1 * self.dt
214456 energy_log.append(EnergyRecord(t, ke, em, eb, ks, woh))
215457 log.info(" t = %6.3f s KE %.3e membrane %.3e bending %.3e exciters %.3e J "
216458 "max|w|/h = %5.2f [%.0f s wall]", t, ke, em, eb, ks, woh, time.time() - t_run)
217459 if not math.isfinite(ke + em + eb):
218460 raise FloatingPointError("energy is not finite: the run has blown up "
219461 "(lower numerics.cfl_safety)")
220-
221- q, v, zs, vs, bf, bb, Fmax, ncont = jax.device_get(state)
222- hits = []
223- for k, p in enumerate(self.patches):
224- hits.append(HitSummary(p.time, p.hit.exciter.kind, float(Fmax[k]), float(ncont[k]) * self.dt,
225- float(vs[k]), p.energy, 0.5 * p.mass * float(vs[k]) ** 2))
226- log.info("hit %d (%s): peak force %.1f N, contact %.2f ms, rebound velocity %+.2f m/s "
227- "(energy left in exciter %.2f of %.2f J)", k, p.hit.exciter.kind, hits[-1].peak_force,
228- hits[-1].contact_time * 1e3, hits[-1].rebound_velocity, hits[-1].energy_out, p.energy)
229-
230- n_out = int(self.fs * num.duration)
231- p_front = self.radiation.to_audio(bf, self.n_steps, self.substeps, n_out)
232- p_back = self.radiation.to_audio(bb, self.n_steps, self.substeps, n_out)
233- result = SimulationResult(self.config, self.fs, self.dt, self.substeps, self.mesh.n_nodes,
234- self.mesh.n_elements, p_front, p_front - p_back, hits, energy_log,
235- time.time() - t_run)
462+ if writer is not None and save_every is not None and time.time() - t_saved >= save_every \
463+ and blk < self.n_blocks - 1:
464+ res = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),
465+ acc, disp, n1, time.time() - t_run, dec_f, dec_b)
466+ writer.update(res)
467+ t_saved = time.time()
468+ log.info(" saved partial run (%.3f s) to %s", res.duration_done, save_path)
469+
470+ result = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),
471+ acc, disp, self.n_steps, time.time() - t_run, dec_f, dec_b)
472+ for h in result.hits:
473+ log.info("hit (%s): peak force %.1f N, contact %.2f ms, rebound velocity %+.2f m/s "
474+ "(energy left in exciter %.2f of %.2f J)", h.kind, h.peak_force,
475+ h.contact_time * 1e3, h.rebound_velocity, h.energy_out, h.energy_in)
236476 for name, p in (("rayleigh", result.pressure_rayleigh), ("unbaffled", result.pressure_unbaffled)):
237477 peak = float(np.abs(p).max())
238478 log.info("%-10s: peak %.3f Pa (%.1f dB SPL peak)", name, peak, 20 * np.log10(max(peak, 1e-12) / 2e-5))
479+ if writer is not None:
480+ writer.update(result)
481+ writer.close()
482+ log.info("wrote %s", save_path)
239483 return result
240484
241485
tests/test_cymsim_file.pyadded+163−0View file
@@ -0,0 +1,163 @@
1+import copy
2+import dataclasses
3+
4+import h5py
5+import numpy as np
6+import pytest
7+
8+from cymsim.cli import main
9+from cymsim.cymsim_file import CymsimWriter, describe, load_cymsim, save_cymsim
10+from cymsim.decimation import DECIMATION_SPAN
11+from cymsim.simulate import AccelerationRecord, DisplacementRecord, Simulator
12+
13+
14+@pytest.fixture(scope="module")
15+def recorded(coarse_config):
16+ cfg = copy.deepcopy(coarse_config)
17+ cfg.numerics.duration = 0.05
18+ cfg.output.record_acceleration = True
19+ cfg.output.displacement_rate = 200.0
20+ sim = Simulator(cfg)
21+ return sim, sim.run()
22+
23+
24+def test_records_have_the_expected_shapes(recorded):
25+ sim, res = recorded
26+ n_out = int(round(res.fs * res.config.numerics.duration))
27+ assert res.pressure_rayleigh.shape == (n_out,)
28+ assert res.acceleration.data.shape == (n_out, res.n_nodes) and res.acceleration.data.dtype == np.int16
29+ assert res.acceleration.scale.shape == (n_out,)
30+ assert res.acceleration.taps.shape == (DECIMATION_SPAN, res.substeps)
31+ assert res.displacement.w.shape == (11, res.n_nodes) # 0, 5 ms, ..., 50 ms
32+ assert np.allclose(res.displacement.time, np.arange(11) * 0.005, atol=res.dt)
33+ assert res.complete and res.n_steps_done == sim.n_steps
34+ assert res.state is not None and res.state.acceleration_ring.shape == (DECIMATION_SPAN, res.n_nodes)
35+ # the record is quiet before the stroke reaches the plate and loud after it
36+ a = res.acceleration.as_float()
37+ assert np.abs(a[:5]).max() < 1e-3 * np.abs(a).max()
38+ assert np.abs(a).max() > 1e3
39+
40+
41+def test_normal_acceleration_and_displacement_are_consistent(recorded):
42+ # the recorded normal acceleration, doubly integrated, should track the displacement frames
43+ # only loosely (the record is band-limited); instead check the scale: w ~ a / omega^2
44+ sim, res = recorded
45+ a = res.acceleration.as_float()
46+ node = int(np.argmax(np.abs(res.displacement.w[-1])))
47+ w_peak = np.abs(res.displacement.w[:, node]).max()
48+ a_peak = np.abs(a[:, node]).max()
49+ omega = np.sqrt(a_peak / w_peak)
50+ assert 2 * np.pi * 50 < omega < 2 * np.pi * 20000
51+
52+
53+def test_baffled_pressure_follows_from_the_recorded_accelerations(recorded):
54+ """A downstream tool sees only the file: node accelerations at fs, mesh, listener distances.
55+ Rebuilding the Rayleigh sum from those with a band-limited fractional delay must reproduce the
56+ stored baffled signal to a few percent rms. This is the check that the record is usable."""
57+ sim, res = recorded
58+ a = res.acceleration.as_float().astype(np.float64) # (T, N)
59+ tri, area = res.mesh.tri, res.mesh.area
60+ L = res.config.listener
61+ a_e = a[:, tri].mean(2) # (T, E)
62+ coef = L.rho_air * area / (2 * np.pi * res.r_front)
63+ delay = res.r_front / L.c_air * res.fs # audio samples
64+ T, E = a_e.shape
65+ K = 16 # half-width of the interpolator
66+ p = np.zeros(T)
67+ for e in range(E):
68+ i0, fr = int(np.floor(delay[e])), delay[e] - np.floor(delay[e])
69+ k = np.arange(-K, K + 1)
70+ x = k - fr
71+ taps = np.sinc(x) * np.kaiser(2 * K + 1, 8.0) # windowed-sinc fractional delay
72+ taps /= taps.sum()
73+ # p[j] += coef sum_k taps[k] a_e[j - i0 - k]
74+ for kk, tap in zip(k, taps):
75+ s = i0 + kk
76+ if 0 <= s < T:
77+ p[s:] += coef[e] * tap * a_e[:T - s, e]
78+ elif s < 0:
79+ p[:T + s] += coef[e] * tap * a_e[-s:, e]
80+ ref = res.pressure_rayleigh
81+ n = min(len(p), len(ref)) - K
82+ err = np.linalg.norm(p[:n] - ref[:n]) / np.linalg.norm(ref[:n])
83+ assert err < 1e-3, err
84+
85+
86+def _truncate(res, n_samples: int, n_frames: int):
87+ acc = AccelerationRecord(res.fs, res.acceleration.data[:n_samples], res.acceleration.scale[:n_samples],
88+ res.acceleration.taps)
89+ disp = DisplacementRecord(res.displacement.rate, res.displacement.time[:n_frames], res.displacement.w[:n_frames])
90+ return dataclasses.replace(res, pressure_rayleigh=res.pressure_rayleigh[:n_samples],
91+ pressure_unbaffled=res.pressure_unbaffled[:n_samples], acceleration=acc,
92+ displacement=disp, energy_log=res.energy_log[:1], complete=False,
93+ n_steps_done=n_samples * res.substeps)
94+
95+
96+def test_round_trip_and_partial_files(recorded, tmp_path):
97+ sim, res = recorded
98+ path = tmp_path / "run.cymsim"
99+ save_cymsim(res, path)
100+ back = load_cymsim(path)
101+ assert back.config == res.config
102+ assert back.fs == res.fs and back.substeps == res.substeps and np.isclose(back.dt, res.dt)
103+ assert np.array_equal(back.mesh.tri, res.mesh.tri) and np.allclose(back.mesh.xy, res.mesh.xy)
104+ assert np.allclose(back.normal_node, res.normal_node) and np.allclose(back.z0_node, res.z0_node)
105+ assert np.allclose(back.pressure_unbaffled, res.pressure_unbaffled, atol=1e-6 * np.abs(res.pressure_unbaffled).max())
106+ assert np.array_equal(back.acceleration.data, res.acceleration.data)
107+ assert np.array_equal(back.acceleration.scale, res.acceleration.scale)
108+ assert np.array_equal(back.displacement.w, res.displacement.w)
109+ assert back.hits == res.hits and len(back.energy_log) == len(res.energy_log)
110+ assert back.complete and back.state.n_steps_done == res.n_steps_done
111+ assert np.array_equal(back.state.q, res.state.q)
112+ assert back.creator.startswith("cymsim")
113+ assert "complete" in describe(path)
114+
115+ # a partial file is a prefix of the complete one, and extending it gives the complete file
116+ part = tmp_path / "partial.cymsim"
117+ trunc = _truncate(res, 1000, 4)
118+ w = CymsimWriter(part, trunc)
119+ w.update(trunc)
120+ w.close()
121+ p = load_cymsim(part)
122+ assert not p.complete and p.n_steps_done == 1000 * res.substeps
123+ assert np.array_equal(p.acceleration.data, res.acceleration.data[:1000])
124+ assert len(p.pressure_rayleigh) == 1000 and p.displacement.n_frames == 4
125+ assert "PARTIAL" in describe(part)
126+ w = CymsimWriter(part, trunc)
127+ w.update(trunc)
128+ w.update(res)
129+ w.close()
130+ full = load_cymsim(part)
131+ assert full.complete and np.array_equal(full.acceleration.data, res.acceleration.data)
132+ assert np.allclose(full.pressure_rayleigh, res.pressure_rayleigh, atol=1e-6)
133+ assert full.displacement.n_frames == res.displacement.n_frames
134+
135+ with h5py.File(path) as f:
136+ assert f.attrs["format"] == "cymsim"
137+ assert f["acceleration/normal"].chunks[1] == res.n_nodes
138+ assert f["acceleration"].attrs["component"] == "normal"
139+
140+
141+def test_no_records_by_default(coarse_config, tmp_path):
142+ res = Simulator(coarse_config).run()
143+ assert res.acceleration is None and res.displacement is not None # 60 Hz frames by default
144+ assert res.displacement.n_frames == 2 # t = 0 and 1/60 s in a 20 ms run
145+ path = tmp_path / "plain.cymsim"
146+ res.save(str(path))
147+ with h5py.File(path) as f:
148+ assert "acceleration" not in f and "displacement" in f
149+ back = load_cymsim(path)
150+ assert back.acceleration is None and np.allclose(back.pressure_rayleigh, res.pressure_rayleigh, atol=1e-6)
151+
152+
153+def test_cli_run_and_info(coarse_config, tmp_path, capsys):
154+ cfg = copy.deepcopy(coarse_config)
155+ cfg.output.wav = cfg.output.npz = cfg.output.cascade_plot = False
156+ cfg.output.prefix = str(tmp_path / "x")
157+ ypath = tmp_path / "c.yaml"
158+ cfg.to_yaml(ypath)
159+ out = tmp_path / "run.cymsim"
160+ assert main(["run", str(ypath), "-o", str(out), "--save-every", "0", "--record-acceleration", "--quiet"]) == 0
161+ assert main(["info", str(out)]) == 0
162+ text = capsys.readouterr().out
163+ assert "acceleration (normal, int16)" in text and "complete" in text
tests/test_decimation.pyadded+52−0View file
@@ -0,0 +1,52 @@
1+import numpy as np
2+
3+from cymsim.decimation import DECIMATION_SPAN, StreamingDecimator, decimation_kernel, decimation_taps
4+
5+
6+def test_kernel_is_normalised_linear_phase_and_attenuates_the_stopband():
7+ m, fs = 54, 44100
8+ h = decimation_kernel(m)
9+ W = DECIMATION_SPAN
10+ assert len(h) == W * m and np.isclose(h.sum(), 1.0)
11+ d = np.arange(-W * m // 2, W * m // 2)
12+ assert np.allclose(h[1:], h[::-1][:-1]) # symmetric about d = 0 (odd length ignored)
13+ f = np.array([1000.0, 15000.0, 18000.0, 22050.0, 24000.0, 30000.0])
14+ H = np.abs(np.exp(-2j * np.pi * np.outer(f / (fs * m), d)) @ h)
15+ assert H[0] > 0.999 and H[1] > 0.99 and H[2] > 0.95 # passband
16+ assert H[3] < 0.1 and H[4] < 1e-3 and H[5] < 1e-4 # transition and stopband
17+
18+
19+def test_polyphase_layout_matches_the_kernel():
20+ m = 10
21+ h, taps = decimation_kernel(m), decimation_taps(m)
22+ W = DECIMATION_SPAN
23+ assert taps.shape == (W, m)
24+ for w, kk in [(0, 0), (7, 3), (W - 1, m - 1)]:
25+ d = kk + (W // 2 - 1 - w) * m
26+ assert taps[w, kk] == h[d + W * m // 2]
27+ assert np.isclose(taps.sum(), 1.0)
28+
29+
30+def test_streaming_decimator_matches_direct_convolution_for_any_chunking():
31+ m = 7
32+ h = decimation_kernel(m)
33+ W = DECIMATION_SPAN
34+ rng = np.random.default_rng(0)
35+ x = rng.standard_normal(m * 400)
36+ xz = np.concatenate([np.zeros(W * m // 2), x]) # rest before t = 0
37+ n_final = (len(x) - W * m // 2) // m + 1
38+ direct = np.array([xz[j * m:j * m + W * m] @ h for j in range(n_final)])
39+ dec = StreamingDecimator(m)
40+ out, i = [], 0
41+ for L in [1, 3, 50, 999, 7, 300, 2000]:
42+ out.append(dec.push(x[i:i + L]))
43+ i += L
44+ out = np.concatenate(out)
45+ assert len(out) == n_final
46+ assert np.abs(out - direct).max() < 1e-12
47+ # resuming from the saved tail continues the stream exactly
48+ dec2 = StreamingDecimator(m, tail=dec.tail, n_seen=dec.n_seen)
49+ more = rng.standard_normal(m * 50)
50+ a = dec.push(more)
51+ b = dec2.push(more)
52+ assert len(a) == 50 and np.array_equal(a, b)
tools/make_web_fixture.pyadded+70−0View file
@@ -0,0 +1,70 @@
1+"""Refresh the run-dependent parts of cymsimweb's Python reference fixture.
2+
3+Usage: python tools/make_web_fixture.py ../cymsimweb/tests/fixtures/python_reference.json
4+
5+The fixture's configuration (a coarse, dent-free profile) is kept; the ``config``, ``timestep``,
6+``run`` and ``run_linear`` entries are rewritten from a fresh run of this package, including the
7+recorded normal acceleration of a few nodes and a sample of the decimation taps, which the web
8+tests compare against their own implementation. The mesh, shell, DKT and patch entries do not
9+depend on the time loop and are left as they are.
10+"""
11+from __future__ import annotations
12+
13+import dataclasses
14+import json
15+import sys
16+
17+import numpy as np
18+
19+from cymsim.config import SimulationConfig
20+from cymsim.simulate import Simulator
21+
22+ACC_NODES = [0, 50, 120]
23+
24+
25+def run_entry(cfg: SimulationConfig) -> dict:
26+ sim = Simulator(cfg)
27+ res = sim.run()
28+ a = res.acceleration.as_float()
29+ return {
30+ "hits": [dataclasses.asdict(h) for h in res.hits],
31+ "energy_log": [[e.time, e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h] for e in res.energy_log],
32+ "pressure_rayleigh": res.pressure_rayleigh.tolist(),
33+ "pressure_unbaffled": res.pressure_unbaffled.tolist(),
34+ "fs": res.fs,
35+ "acceleration_nodes": ACC_NODES,
36+ "acceleration": [a[:, i].astype(float).tolist() for i in ACC_NODES],
37+ "n_steps": sim.n_steps,
38+ }, sim
39+
40+
41+def main(path: str) -> None:
42+ with open(path) as f:
43+ fx = json.load(f)
44+ cfg = SimulationConfig.from_dict(fx["config"])
45+ cfg.output.record_acceleration = True
46+ cfg.output.displacement_rate = 0.0
47+ sim0 = Simulator(cfg)
48+ # whole blocks must end on an audio sample so that the web loop can stop at the same time
49+ m = sim0.substeps
50+ if cfg.numerics.steps_per_block % m:
51+ cfg.numerics.steps_per_block = m * max(1, cfg.numerics.steps_per_block // m)
52+ fx["config"] = cfg.to_dict()
53+ run, sim = run_entry(cfg)
54+ lin_cfg = SimulationConfig.from_dict(cfg.to_dict())
55+ lin_cfg.numerics.nonlinear = False
56+ run_linear, _ = run_entry(lin_cfg)
57+ fx["timestep"] = {"dt": sim.dt, "substeps": sim.substeps}
58+ fx["run"], fx["run_linear"] = run, run_linear
59+ taps = sim.taps
60+ fx["decimation"] = {"span": taps.shape[0], "substeps": taps.shape[1],
61+ "taps_sample": [[int(w), int(k), float(taps[w, k])] for w, k in
62+ [(0, 0), (15, 0), (15, m - 1), (16, 0), (31, m - 1), (7, m // 2)]],
63+ "normal_first3": sim.shell.normal_node[:3].tolist()}
64+ with open(path, "w") as f:
65+ json.dump(fx, f)
66+ print(f"wrote {path}: {len(run['pressure_unbaffled'])} samples, {sim.n_steps} steps, m = {m}")
67+
68+
69+if __name__ == "__main__":
70+ main(sys.argv[1])