1"""The tangent of a run (Simulator._tangent_block) is the derivative of the state with respect to a
2relative change of the prescribed peak force: it matches central differences of two twin runs."""
3import copy
5import jax.numpy as jnp
6import numpy as np
8from cymsim.config import ForceExciter
9from cymsim.simulate import Simulator
11DELTA = 1e-6
14def _flat_struck(base, peak):
15 cfg = copy.deepcopy(base)
16 cfg.profile.bow_height = cfg.profile.bell_height = 0.0
17 cfg.profile.h_shoulder = cfg.profile.h_bell = cfg.profile.h_rim
18 cfg.imperfections.n_dents = 0
19 cfg.hits[0].exciter = ForceExciter(peak=peak, duration=1e-3)
20 cfg.numerics.dtype = "float64"
21 return cfg
24def _blocks(sim, n_blocks):
25 """q and v after each block of a plain run."""
26 spb, st, out = sim.config.numerics.steps_per_block, sim.initial_state(), []
27 for b in range(n_blocks):
28 st = sim._run_block(jnp.int32(b * spb), jnp.int32(0), jnp.int32(0), st)
29 out.append((np.asarray(st[0]), np.asarray(st[1])))
30 return out
33def test_tangent_matches_central_differences(coarse_config):
34 F, n_blocks = 1000.0, 6
35 sim = Simulator(_flat_struck(coarse_config, F))
36 spb = sim.config.numerics.steps_per_block
37 st = sim.initial_state()
38 dst = Simulator.zero_tangent(st)
39 tangent, primal = [], []
40 for b in range(n_blocks):
41 st, dst = sim._tangent_block(jnp.int32(b * spb), jnp.int32(0), jnp.int32(0), st, dst)
42 primal.append((np.asarray(st[0]), np.asarray(st[1])))
43 tangent.append((np.asarray(dst[0]), np.asarray(dst[1])))
44 plain = _blocks(sim, n_blocks)
45 plus = _blocks(Simulator(_flat_struck(coarse_config, F * (1 + DELTA))), n_blocks)
46 minus = _blocks(Simulator(_flat_struck(coarse_config, F * (1 - DELTA))), n_blocks)
47 for b in range(n_blocks):
48 for i in (0, 1): # q, v
49 assert np.allclose(primal[b][i], plain[b][i], rtol=0, atol=1e-12 * np.abs(plain[b][i]).max())
50 fd = (plus[b][i] - minus[b][i]) / (2 * DELTA)
51 err = np.linalg.norm(tangent[b][i] - fd) / np.linalg.norm(fd)
52 assert err < 1e-5, (b, i, err)
55def _tangent_config(base, peak, duration, tangent=True):
56 cfg = _flat_struck(base, peak)
57 cfg.numerics.duration = duration
58 cfg.numerics.tangent = tangent
59 cfg.output.displacement_rate = 500.0
60 cfg.output.probes = [(0.1, 0.05)]
61 return cfg
64def test_tangent_run_outputs_match_central_differences(coarse_config):
65 """Sound, probe and displacement frames of the tangent equal (twin+ - twin-) / (2 delta), and the
66 run's own outputs equal those of a plain run."""
67 F, T = 1000.0, 0.03
68 res = Simulator(_tangent_config(coarse_config, F, T)).run()
69 plain = Simulator(_tangent_config(coarse_config, F, T, tangent=False)).run()
70 plus = Simulator(_tangent_config(coarse_config, F * (1 + DELTA), T, tangent=False)).run()
71 minus = Simulator(_tangent_config(coarse_config, F * (1 - DELTA), T, tangent=False)).run()
72 tg = res.tangent
73 assert res.complete and tg is not None and len(tg.pressure_unbaffled) == len(res.pressure_unbaffled)
74 pairs = [(res.pressure_unbaffled, plain.pressure_unbaffled, tg.pressure_unbaffled,
75 plus.pressure_unbaffled, minus.pressure_unbaffled),
76 (res.pressure_rayleigh, plain.pressure_rayleigh, tg.pressure_rayleigh,
77 plus.pressure_rayleigh, minus.pressure_rayleigh),
78 (res.probe_acceleration, plain.probe_acceleration, tg.probe_acceleration,
79 plus.probe_acceleration, minus.probe_acceleration),
80 (res.displacement.w, plain.displacement.w, tg.displacement.w, plus.displacement.w, minus.displacement.w)]
81 for own, ref, tan, p, m in pairs:
82 assert np.allclose(own, ref, rtol=0, atol=1e-12 * np.abs(ref).max())
83 fd = (p - m) / (2 * DELTA)
84 assert np.linalg.norm(tan - fd) / np.linalg.norm(fd) < 1e-5
85 assert np.array_equal(tg.displacement.time, res.displacement.time)
86 norms = tg.norms
87 assert norms.shape[1] == 4 and len(norms) == res.n_steps_done // res.config.numerics.steps_per_block
88 assert np.all(norms[:, 1:] >= 0) and norms[-1, 1] > 0
91def test_tangent_run_continues_and_cuts(coarse_config, tmp_path):
92 """A tangent run continued from its .cymsim file equals the uninterrupted one, and a shorter
93 tangent run is cut from a longer one as from a plain run."""
94 from cymsim.cymsim_file import load_cymsim
95 from cymsim.simulate import cut_run
96 straight = Simulator(_tangent_config(coarse_config, 1000.0, 0.03)).run(save_path=str(tmp_path / "a.cymsim"))
97 Simulator(_tangent_config(coarse_config, 1000.0, 0.015)).run(save_path=str(tmp_path / "b.cymsim"))
98 cont = Simulator(_tangent_config(coarse_config, 1000.0, 0.03)).run(resume=load_cymsim(tmp_path / "b.cymsim"))
99 for a, b in ((straight, cont), (straight.tangent, cont.tangent)):
100 for name in ("pressure_rayleigh", "pressure_unbaffled", "probe_acceleration"):
101 assert np.array_equal(getattr(a, name), getattr(b, name)), name
102 assert np.array_equal(a.displacement.w, b.displacement.w)
103 assert np.array_equal(straight.tangent.norms, cont.tangent.norms)
104 loaded = load_cymsim(tmp_path / "a.cymsim")
105 assert np.array_equal(loaded.tangent.pressure_unbaffled, straight.tangent.pressure_unbaffled)
106 assert np.array_equal(loaded.tangent.norms, straight.tangent.norms)
107 short = cut_run(loaded, _tangent_config(coarse_config, 1000.0, 0.015))
108 direct = load_cymsim(tmp_path / "b.cymsim")
109 for a, b in ((short, direct), (short.tangent, direct.tangent)):
110 assert np.array_equal(a.pressure_unbaffled, b.pressure_unbaffled)
111 assert np.array_equal(a.displacement.w, b.displacement.w)
112 assert np.array_equal(short.tangent.norms, direct.tangent.norms)
113 assert np.array_equal(short.tangent.state.q, loaded.tangent.state.q)