import copy import dataclasses import h5py import numpy as np import pytest from cymsim.cli import main from cymsim.cymsim_file import CymsimWriter, describe, load_cymsim, save_cymsim from cymsim.decimation import DECIMATION_SPAN from cymsim.simulate import AccelerationRecord, DisplacementRecord, Simulator @pytest.fixture(scope="module") def recorded(coarse_config): cfg = copy.deepcopy(coarse_config) cfg.numerics.duration = 0.05 cfg.output.record_acceleration = True cfg.output.displacement_rate = 200.0 sim = Simulator(cfg) return sim, sim.run() def test_records_have_the_expected_shapes(recorded): sim, res = recorded n_out = int(round(res.fs * res.config.numerics.duration)) assert res.pressure_rayleigh.shape == (n_out,) assert res.acceleration.data.shape == (n_out, res.n_nodes) and res.acceleration.data.dtype == np.int16 assert res.acceleration.scale.shape == (n_out,) assert res.acceleration.taps.shape == (DECIMATION_SPAN, res.substeps) assert res.displacement.w.shape == (11, res.n_nodes) # 0, 5 ms, ..., 50 ms assert np.allclose(res.displacement.time, np.arange(11) * 0.005, atol=res.dt) assert res.complete and res.n_steps_done == sim.n_steps assert res.state is not None and res.state.acceleration_ring.shape == (DECIMATION_SPAN, res.n_nodes) # the record is quiet before the stroke reaches the plate and loud after it a = res.acceleration.as_float() assert np.abs(a[:5]).max() < 1e-3 * np.abs(a).max() assert np.abs(a).max() > 1e3 def test_normal_acceleration_and_displacement_are_consistent(recorded): # the recorded normal acceleration, doubly integrated, should track the displacement frames # only loosely (the record is band-limited); instead check the scale: w ~ a / omega^2 sim, res = recorded a = res.acceleration.as_float() node = int(np.argmax(np.abs(res.displacement.w[-1]))) w_peak = np.abs(res.displacement.w[:, node]).max() a_peak = np.abs(a[:, node]).max() omega = np.sqrt(a_peak / w_peak) assert 2 * np.pi * 50 < omega < 2 * np.pi * 20000 def test_baffled_pressure_follows_from_the_recorded_accelerations(recorded): """A downstream tool sees only the file: node accelerations at fs, mesh, listener distances. Rebuilding the Rayleigh sum from those with a band-limited fractional delay must reproduce the stored baffled signal to a few percent rms. This is the check that the record is usable.""" sim, res = recorded a = res.acceleration.as_float().astype(np.float64) # (T, N) tri, area = res.mesh.tri, res.mesh.area L = res.config.listener a_e = a[:, tri].mean(2) # (T, E) coef = L.rho_air * area / (2 * np.pi * res.r_front) delay = res.r_front / L.c_air * res.fs # audio samples T, E = a_e.shape K = 16 # half-width of the interpolator p = np.zeros(T) for e in range(E): i0, fr = int(np.floor(delay[e])), delay[e] - np.floor(delay[e]) k = np.arange(-K, K + 1) x = k - fr taps = np.sinc(x) * np.kaiser(2 * K + 1, 8.0) # windowed-sinc fractional delay taps /= taps.sum() # p[j] += coef sum_k taps[k] a_e[j - i0 - k] for kk, tap in zip(k, taps): s = i0 + kk if 0 <= s < T: p[s:] += coef[e] * tap * a_e[:T - s, e] elif s < 0: p[:T + s] += coef[e] * tap * a_e[-s:, e] ref = res.pressure_rayleigh n = min(len(p), len(ref)) - K err = np.linalg.norm(p[:n] - ref[:n]) / np.linalg.norm(ref[:n]) assert err < 1e-3, err def _truncate(res, n_samples: int, n_frames: int): acc = AccelerationRecord(res.fs, res.acceleration.data[:n_samples], res.acceleration.scale[:n_samples], res.acceleration.taps) disp = DisplacementRecord(res.displacement.rate, res.displacement.time[:n_frames], res.displacement.w[:n_frames]) return dataclasses.replace(res, pressure_rayleigh=res.pressure_rayleigh[:n_samples], pressure_unbaffled=res.pressure_unbaffled[:n_samples], acceleration=acc, displacement=disp, energy_log=res.energy_log[:1], complete=False, n_steps_done=n_samples * res.substeps) def test_round_trip_and_partial_files(recorded, tmp_path): sim, res = recorded path = tmp_path / "run.cymsim" save_cymsim(res, path) back = load_cymsim(path) assert back.config == res.config assert back.fs == res.fs and back.substeps == res.substeps and np.isclose(back.dt, res.dt) assert np.array_equal(back.mesh.tri, res.mesh.tri) and np.allclose(back.mesh.xy, res.mesh.xy) assert np.allclose(back.normal_node, res.normal_node) and np.allclose(back.z0_node, res.z0_node) assert np.allclose(back.pressure_unbaffled, res.pressure_unbaffled, atol=1e-6 * np.abs(res.pressure_unbaffled).max()) assert np.array_equal(back.acceleration.data, res.acceleration.data) assert np.array_equal(back.acceleration.scale, res.acceleration.scale) assert np.array_equal(back.displacement.w, res.displacement.w) assert back.hits == res.hits and len(back.energy_log) == len(res.energy_log) assert back.complete and back.state.n_steps_done == res.n_steps_done assert np.array_equal(back.state.q, res.state.q) assert back.creator.startswith("cymsim") assert "complete" in describe(path) # a partial file is a prefix of the complete one, and extending it gives the complete file part = tmp_path / "partial.cymsim" trunc = _truncate(res, 1000, 4) w = CymsimWriter(part, trunc) w.update(trunc) w.close() p = load_cymsim(part) assert not p.complete and p.n_steps_done == 1000 * res.substeps assert np.array_equal(p.acceleration.data, res.acceleration.data[:1000]) assert len(p.pressure_rayleigh) == 1000 and p.displacement.n_frames == 4 assert "PARTIAL" in describe(part) w = CymsimWriter(part, trunc) w.update(trunc) w.update(res) w.close() full = load_cymsim(part) assert full.complete and np.array_equal(full.acceleration.data, res.acceleration.data) assert np.allclose(full.pressure_rayleigh, res.pressure_rayleigh, atol=1e-6) assert full.displacement.n_frames == res.displacement.n_frames with h5py.File(path) as f: assert f.attrs["format"] == "cymsim" assert f["acceleration/normal"].chunks[1] == res.n_nodes assert f["acceleration"].attrs["component"] == "normal" def test_no_records_by_default(coarse_config, tmp_path): res = Simulator(coarse_config).run() assert res.acceleration is None and res.displacement is not None # 60 Hz frames by default assert res.displacement.n_frames == 2 # t = 0 and 1/60 s in a 20 ms run path = tmp_path / "plain.cymsim" res.save(str(path)) with h5py.File(path) as f: assert "acceleration" not in f and "displacement" in f back = load_cymsim(path) assert back.acceleration is None and np.allclose(back.pressure_rayleigh, res.pressure_rayleigh, atol=1e-6) def test_cli_run_and_info(coarse_config, tmp_path, capsys): cfg = copy.deepcopy(coarse_config) cfg.output.wav = cfg.output.npz = cfg.output.cascade_plot = False cfg.output.prefix = str(tmp_path / "x") ypath = tmp_path / "c.yaml" cfg.to_yaml(ypath) out = tmp_path / "run.cymsim" assert main(["run", str(ypath), "-o", str(out), "--save-every", "0", "--record-acceleration", "--quiet"]) == 0 assert main(["info", str(out)]) == 0 text = capsys.readouterr().out assert "acceleration (normal, int16)" in text and "complete" in text