Store signals (audio, probes, displacement frames) in the run's precision, float64 for a float64 run; run the tests on the CPU, where runs are reproducible bit for bit
5 changed files+52−18
docs/format.mdmodified+6−6View file
@@ -8,7 +8,7 @@ The purpose of the acceleration record is to let another program compute the sou
88
99 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.
1010
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.
11+Integers are stored as `int32`, geometry as `float64`, signals (audio, probes, displacement frames) in the precision of the run, `float32` or, for a run with `numerics.dtype: float64`, `float64`, and the acceleration record as `int16` against a scale per sample (section 5). Strings are variable-length UTF-8. Attributes are used for scalars and short vectors, datasets for arrays.
1212
1313 ## 2. Root
1414
@@ -88,10 +88,10 @@ Present when `output.displacement_rate` is positive (default 60 Hz). Group attri
8888
8989 | dataset | shape | meaning |
9090 |---|---|---|
91-| `w` | $(T_d, N)$ float32 | transverse displacement $w$ of every node at the frame times |
91+| `w` | $(T_d, N)$ float32/64 | transverse displacement $w$ of every node at the frame times |
9292 | `time` | $(T_d,)$ float64 | frame times [s] |
9393
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.
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; twice that for a float64 run). The web app's slow-motion replay uses its own capture at a much higher frame rate and does not write it to the file.
9595
9696 ## 8. Listener signals: `audio`
9797
@@ -99,8 +99,8 @@ Group attributes `fs`, `units = "Pa"`.
9999
100100 | dataset | shape | meaning |
101101 |---|---|---|
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 |
102+| `baffled` | $(T_a,)$ float32/64 | Rayleigh integral of the front face, as if in an infinite baffle |
103+| `unbaffled` | $(T_a,)$ float32/64 | front minus `back_diffraction_gain` times the diffracted back-face term |
104104
105105 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.
106106
@@ -109,7 +109,7 @@ Both are the pressure at `listener/position`, sampled at $t_j = j / f_s$, filter
109109 | dataset | shape | meaning |
110110 |---|---|---|
111111 | `xy` | $(P, 2)$ float64 | the probe points in the projected plane |
112-| `normal_acceleration` | $(T_p, P)$ float32 | normal acceleration at each point, interpolated linearly in the containing element and filtered by the same kernel as the listener signals |
112+| `normal_acceleration` | $(T_p, P)$ float32/64 | normal acceleration at each point, interpolated linearly in the containing element and filtered by the same kernel as the listener signals |
113113
114114 ## 9. Reproducing the baffled signal from the record
115115
src/cymsim/cymsim_file.pymodified+14−7View file
@@ -37,6 +37,12 @@ OPTIONAL_STATE_ARRAYS = ("acceleration_ring", "decimator_tail_probes", "overhang
3737 "overhang_displacement_time", "overhang_displacement_w")
3838
3939
40+def signal_dtype(config: SimulationConfig):
41+ """Signals (audio, probes, displacement frames) are stored in the run's precision, so that the
42+ small differences between runs of a float64 simulation survive."""
43+ return np.float64 if config.numerics.dtype == "float64" else np.float32
44+
45+
4046 class CymsimWriter:
4147 """Create a .cymsim file for a run and extend it as the run proceeds.
4248
@@ -50,6 +56,7 @@ class CymsimWriter:
5056 self.path = str(path)
5157 r = result
5258 N, E = r.n_nodes, r.n_elements
59+ sig = self.signal_dtype = signal_dtype(r.config)
5360 f = self.f = h5py.File(self.path, "w")
5461 f.attrs["format"] = "cymsim"
5562 f.attrs["format_version"] = FORMAT_VERSION
@@ -87,7 +94,7 @@ class CymsimWriter:
8794 g = f.create_group("audio")
8895 g.attrs["fs"], g.attrs["units"] = int(r.fs), "Pa"
8996 for name in ("baffled", "unbaffled"):
90- g.create_dataset(name, shape=(0,), maxshape=(None,), dtype=np.float32, chunks=(65536,))
97+ g.create_dataset(name, shape=(0,), maxshape=(None,), dtype=sig, chunks=(65536,))
9198
9299 if r.acceleration is not None:
93100 a = r.acceleration
@@ -112,7 +119,7 @@ class CymsimWriter:
112119 g.attrs["rate"] = float(r.displacement.rate)
113120 g.attrs["units"] = "m"
114121 g.attrs["description"] = "transverse displacement w of every node at the frame times"
115- g.create_dataset("w", shape=(0, N), maxshape=(None, N), dtype=np.float32, chunks=(64, N))
122+ g.create_dataset("w", shape=(0, N), maxshape=(None, N), dtype=sig, chunks=(64, N))
116123 g.create_dataset("time", shape=(0,), maxshape=(None,), dtype=np.float64, chunks=(4096,))
117124
118125 if r.probes is not None:
@@ -120,7 +127,7 @@ class CymsimWriter:
120127 g.attrs["fs"], g.attrs["units"] = int(r.fs), "m/s^2"
121128 g.create_dataset("xy", data=np.asarray(r.probes, np.float64))
122129 g.create_dataset("normal_acceleration", shape=(0, len(r.probes)), maxshape=(None, len(r.probes)),
123- dtype=np.float32, chunks=(4096, len(r.probes)))
130+ dtype=sig, chunks=(4096, len(r.probes)))
124131
125132 g = f.create_group("energy")
126133 d = g.create_dataset("table", shape=(0, 6), maxshape=(None, 6), dtype=np.float64, chunks=(1024, 6))
@@ -143,16 +150,16 @@ class CymsimWriter:
143150
144151 def update(self, result: SimulationResult) -> None:
145152 f, r = self.f, result
146- self._append(f["audio/baffled"], r.pressure_rayleigh.astype(np.float32))
147- self._append(f["audio/unbaffled"], r.pressure_unbaffled.astype(np.float32))
153+ self._append(f["audio/baffled"], r.pressure_rayleigh.astype(self.signal_dtype))
154+ self._append(f["audio/unbaffled"], r.pressure_unbaffled.astype(self.signal_dtype))
148155 if r.acceleration is not None and "acceleration" in f:
149156 self._append(f["acceleration/normal"], r.acceleration.data)
150157 self._append(f["acceleration/scale"], r.acceleration.scale)
151158 if r.displacement is not None and "displacement" in f:
152- self._append(f["displacement/w"], r.displacement.w)
159+ self._append(f["displacement/w"], r.displacement.w.astype(self.signal_dtype))
153160 self._append(f["displacement/time"], r.displacement.time)
154161 if r.probe_acceleration is not None and "probes" in f:
155- self._append(f["probes/normal_acceleration"], r.probe_acceleration.astype(np.float32))
162+ self._append(f["probes/normal_acceleration"], r.probe_acceleration.astype(self.signal_dtype))
156163 table = np.array([[e.time, e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h]
157164 for e in r.energy_log], np.float64).reshape(-1, 6)
158165 f["energy/table"].resize(len(table), axis=0)
src/cymsim/simulate.pymodified+5−5View file
@@ -114,15 +114,15 @@ class DisplacementRecord:
114114 """Transverse displacement w of every node [m] at a low frame rate, for display."""
115115 rate: float
116116 time: np.ndarray # (Td,)
117- w: np.ndarray # (Td, N) float32
117+ w: np.ndarray # (Td, N) in the run's precision
118118
119119 @property
120120 def n_frames(self) -> int:
121121 return len(self.time)
122122
123123 @classmethod
124- def empty(cls, rate: float, n_nodes: int) -> "DisplacementRecord":
125- return cls(rate, np.zeros(0), np.zeros((0, n_nodes), np.float32))
124+ def empty(cls, rate: float, n_nodes: int, dtype=np.float32) -> "DisplacementRecord":
125+ return cls(rate, np.zeros(0), np.zeros((0, n_nodes), dtype))
126126
127127 def extend(self, time: np.ndarray, w: np.ndarray) -> None:
128128 self.time = np.concatenate([self.time, time])
@@ -540,7 +540,7 @@ class Simulator:
540540 m, spb, N = self.substeps, num.steps_per_block, self.mesh.n_nodes
541541 n_probes = len(self.probe_xy)
542542 acc = AccelerationRecord.empty(self.fs, N, self.taps) if self.record else None
543- disp = DisplacementRecord.empty(self.disp_rate, N) if self.disp_rate > 0 else None
543+ disp = DisplacementRecord.empty(self.disp_rate, N, np.dtype(self.shell.dtype)) if self.disp_rate > 0 else None
544544 if resume is None:
545545 state, n_done, wall0 = self.initial_state(), 0, 0.0
546546 energy_log: list[EnergyRecord] = []
@@ -612,7 +612,7 @@ class Simulator:
612612 if sel.any():
613613 acc.extend(*AccelerationRecord.encode(rows[sel]))
614614 if self.disp_rate > 0 and k1 > k0:
615- frames = np.asarray(state[10][:k1 - k0], np.float32)
615+ frames = np.asarray(state[10][:k1 - k0], np.dtype(self.shell.dtype))
616616 disp.extend(self.frame_steps[k0:k1] * self.dt, frames)
617617
618618 if blk % report_every == 0 or blk == self.n_blocks - 1:
tests/conftest.pymodified+6−0View file
@@ -1,7 +1,13 @@
11 import logging
2+import os
23
34 import pytest
45
6+# The tests run on the CPU, where runs are reproducible bit for bit (on a GPU, scatter-adds by
7+# atomics make two identical runs differ slightly). JAX_PLATFORMS=cuda runs them on a GPU, where
8+# the tests comparing separate runs exactly are expected to fail.
9+os.environ.setdefault("JAX_PLATFORMS", "cpu")
10+
511 from cymsim.config import SimulationConfig, default_config
612
713 logging.getLogger("cymsim").setLevel(logging.WARNING)
tests/test_resume.pymodified+21−0View file
@@ -115,3 +115,24 @@ def test_cli_resume(coarse_config, straight, tmp_path, capsys):
115115 cfg.to_yaml(ypath)
116116 assert main(["run", str(ypath), "--resume", str(straight), "-o", str(out), "--quiet"]) == 2
117117 assert "hits" in capsys.readouterr().err
118+
119+
120+def test_float64_runs_store_float64_signals(coarse_config, tmp_path):
121+ """A float64 run keeps its signals in float64, and continues exactly."""
122+ def cfg(duration):
123+ c = _config(coarse_config, duration)
124+ c.numerics.dtype = "float64"
125+ return c
126+
127+ straight = tmp_path / "straight.cymsim"
128+ Simulator(cfg(0.035)).run(save_path=str(straight))
129+ with h5py.File(straight) as f:
130+ for name in ("audio/baffled", "audio/unbaffled", "probes/normal_acceleration", "displacement/w"):
131+ assert f[name].dtype == np.float64, name
132+ res = load_cymsim(straight)
133+ assert res.displacement.w.dtype == np.float64
134+ first = tmp_path / "first.cymsim"
135+ Simulator(cfg(0.02)).run(save_path=str(first))
136+ out = tmp_path / "continued.cymsim"
137+ Simulator(cfg(0.035)).run(save_path=str(out), resume=load_cymsim(first))
138+ _assert_same_files(out, straight)