2import dataclasses
4import h5py
5import numpy as np
6import pytest
8from cymsim.cli import main
9from cymsim.cymsim_file import CymsimWriter, describe, load_cymsim, save_cymsim
10from cymsim.decimation import DECIMATION_SPAN
11from cymsim.simulate import AccelerationRecord, DisplacementRecord, Simulator
15def 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()
24def 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
41def 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
53def 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
86def _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)
96def 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)
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
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"
141def 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)
153def 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