"""The tangent of a run (Simulator._tangent_block) is the derivative of the state with respect to a relative change of the prescribed peak force: it matches central differences of two twin runs.""" import copy import jax.numpy as jnp import numpy as np from cymsim.config import ForceExciter from cymsim.simulate import Simulator DELTA = 1e-6 def _flat_struck(base, peak): cfg = copy.deepcopy(base) cfg.profile.bow_height = cfg.profile.bell_height = 0.0 cfg.profile.h_shoulder = cfg.profile.h_bell = cfg.profile.h_rim cfg.imperfections.n_dents = 0 cfg.hits[0].exciter = ForceExciter(peak=peak, duration=1e-3) cfg.numerics.dtype = "float64" return cfg def _blocks(sim, n_blocks): """q and v after each block of a plain run.""" spb, st, out = sim.config.numerics.steps_per_block, sim.initial_state(), [] for b in range(n_blocks): st = sim._run_block(jnp.int32(b * spb), jnp.int32(0), jnp.int32(0), st) out.append((np.asarray(st[0]), np.asarray(st[1]))) return out def test_tangent_matches_central_differences(coarse_config): F, n_blocks = 1000.0, 6 sim = Simulator(_flat_struck(coarse_config, F)) spb = sim.config.numerics.steps_per_block st = sim.initial_state() dst = Simulator.zero_tangent(st) tangent, primal = [], [] for b in range(n_blocks): st, dst = sim._tangent_block(jnp.int32(b * spb), jnp.int32(0), jnp.int32(0), st, dst) primal.append((np.asarray(st[0]), np.asarray(st[1]))) tangent.append((np.asarray(dst[0]), np.asarray(dst[1]))) plain = _blocks(sim, n_blocks) plus = _blocks(Simulator(_flat_struck(coarse_config, F * (1 + DELTA))), n_blocks) minus = _blocks(Simulator(_flat_struck(coarse_config, F * (1 - DELTA))), n_blocks) for b in range(n_blocks): for i in (0, 1): # q, v assert np.allclose(primal[b][i], plain[b][i], rtol=0, atol=1e-12 * np.abs(plain[b][i]).max()) fd = (plus[b][i] - minus[b][i]) / (2 * DELTA) err = np.linalg.norm(tangent[b][i] - fd) / np.linalg.norm(fd) assert err < 1e-5, (b, i, err) def _tangent_config(base, peak, duration, tangent=True): cfg = _flat_struck(base, peak) cfg.numerics.duration = duration cfg.numerics.tangent = tangent cfg.output.displacement_rate = 500.0 cfg.output.probes = [(0.1, 0.05)] return cfg def test_tangent_run_outputs_match_central_differences(coarse_config): """Sound, probe and displacement frames of the tangent equal (twin+ - twin-) / (2 delta), and the run's own outputs equal those of a plain run.""" F, T = 1000.0, 0.03 res = Simulator(_tangent_config(coarse_config, F, T)).run() plain = Simulator(_tangent_config(coarse_config, F, T, tangent=False)).run() plus = Simulator(_tangent_config(coarse_config, F * (1 + DELTA), T, tangent=False)).run() minus = Simulator(_tangent_config(coarse_config, F * (1 - DELTA), T, tangent=False)).run() tg = res.tangent assert res.complete and tg is not None and len(tg.pressure_unbaffled) == len(res.pressure_unbaffled) pairs = [(res.pressure_unbaffled, plain.pressure_unbaffled, tg.pressure_unbaffled, plus.pressure_unbaffled, minus.pressure_unbaffled), (res.pressure_rayleigh, plain.pressure_rayleigh, tg.pressure_rayleigh, plus.pressure_rayleigh, minus.pressure_rayleigh), (res.probe_acceleration, plain.probe_acceleration, tg.probe_acceleration, plus.probe_acceleration, minus.probe_acceleration), (res.displacement.w, plain.displacement.w, tg.displacement.w, plus.displacement.w, minus.displacement.w)] for own, ref, tan, p, m in pairs: assert np.allclose(own, ref, rtol=0, atol=1e-12 * np.abs(ref).max()) fd = (p - m) / (2 * DELTA) assert np.linalg.norm(tan - fd) / np.linalg.norm(fd) < 1e-5 assert np.array_equal(tg.displacement.time, res.displacement.time) norms = tg.norms from cymsim.simulate import N_NORM_COLUMNS assert norms.shape[1] == N_NORM_COLUMNS and len(norms) == res.n_steps_done // res.config.numerics.steps_per_block assert np.all(norms[:, 1:] >= 0) and norms[-1, 1] > 0 # the parts add up to the norms assert np.allclose(norms[:, 1], norms[:, 4] + norms[:, 5] + norms[:, 6], rtol=1e-12) assert np.allclose(norms[:, 3], norms[:, 7] + norms[:, 8], rtol=1e-12) # the run's kinetic, bending and true membrane energy are its energy log's at the reported times for e in res.energy_log: i = int(np.argmin(np.abs(norms[:, 0] - e.time))) if abs(norms[i, 0] - e.time) < 1e-12: assert np.allclose([norms[i, 2], norms[i, 7], norms[i, 9]], [e.kinetic, e.bending, e.membrane], rtol=1e-9) def test_tangent_run_continues_and_cuts(coarse_config, tmp_path): """A tangent run continued from its .cymsim file equals the uninterrupted one, and a shorter tangent run is cut from a longer one as from a plain run.""" from cymsim.cymsim_file import load_cymsim from cymsim.simulate import cut_run straight = Simulator(_tangent_config(coarse_config, 1000.0, 0.03)).run(save_path=str(tmp_path / "a.cymsim")) Simulator(_tangent_config(coarse_config, 1000.0, 0.015)).run(save_path=str(tmp_path / "b.cymsim")) cont = Simulator(_tangent_config(coarse_config, 1000.0, 0.03)).run(resume=load_cymsim(tmp_path / "b.cymsim")) for a, b in ((straight, cont), (straight.tangent, cont.tangent)): for name in ("pressure_rayleigh", "pressure_unbaffled", "probe_acceleration"): assert np.array_equal(getattr(a, name), getattr(b, name)), name assert np.array_equal(a.displacement.w, b.displacement.w) assert np.array_equal(straight.tangent.norms, cont.tangent.norms) loaded = load_cymsim(tmp_path / "a.cymsim") assert np.array_equal(loaded.tangent.pressure_unbaffled, straight.tangent.pressure_unbaffled) assert np.array_equal(loaded.tangent.norms, straight.tangent.norms) short = cut_run(loaded, _tangent_config(coarse_config, 1000.0, 0.015)) direct = load_cymsim(tmp_path / "b.cymsim") for a, b in ((short, direct), (short.tangent, direct.tangent)): assert np.array_equal(a.pressure_unbaffled, b.pressure_unbaffled) assert np.array_equal(a.displacement.w, b.displacement.w) assert np.array_equal(short.tangent.norms, direct.tangent.norms) assert np.array_equal(short.tangent.state.q, loaded.tangent.state.q)