de5d1ceTangent runs (numerics.tangent): integrate the derivative of the run with respect to a relative change of the prescribed peak forces by forward-mode differentiation of the step; its sound, probes, frames and an energy norm per block are stored in the .cymsim file (group tangent), and tangent runs continue and are cut like plain runsJeremy Magland 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)