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
240ff79Tangent runs log the energy parts of run and tangent per block (kinetic, bending, linear membrane; the run's true membrane energy), since the linear membrane energy is inflated at large deflection; tangent cache keys carry TANGENT_VERSION = 2Jeremy Magland 87 from cymsim.simulate import N_NORM_COLUMNS
88 assert norms.shape[1] == N_NORM_COLUMNS and len(norms) == res.n_steps_done // res.config.numerics.steps_per_block
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 89 assert np.all(norms[:, 1:] >= 0) and norms[-1, 1] > 0
240ff79Tangent runs log the energy parts of run and tangent per block (kinetic, bending, linear membrane; the run's true membrane energy), since the linear membrane energy is inflated at large deflection; tangent cache keys carry TANGENT_VERSION = 2Jeremy Magland 90 # the parts add up to the norms
91 assert np.allclose(norms[:, 1], norms[:, 4] + norms[:, 5] + norms[:, 6], rtol=1e-12)
92 assert np.allclose(norms[:, 3], norms[:, 7] + norms[:, 8], rtol=1e-12)
93 # the run's kinetic, bending and true membrane energy are its energy log's at the reported times
94 for e in res.energy_log:
95 i = int(np.argmin(np.abs(norms[:, 0] - e.time)))
96 if abs(norms[i, 0] - e.time) < 1e-12:
97 assert np.allclose([norms[i, 2], norms[i, 7], norms[i, 9]], [e.kinetic, e.bending, e.membrane], rtol=1e-9)
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 98
100def test_tangent_run_continues_and_cuts(coarse_config, tmp_path):
101 """A tangent run continued from its .cymsim file equals the uninterrupted one, and a shorter
102 tangent run is cut from a longer one as from a plain run."""
103 from cymsim.cymsim_file import load_cymsim
104 from cymsim.simulate import cut_run
105 straight = Simulator(_tangent_config(coarse_config, 1000.0, 0.03)).run(save_path=str(tmp_path / "a.cymsim"))
106 Simulator(_tangent_config(coarse_config, 1000.0, 0.015)).run(save_path=str(tmp_path / "b.cymsim"))
107 cont = Simulator(_tangent_config(coarse_config, 1000.0, 0.03)).run(resume=load_cymsim(tmp_path / "b.cymsim"))
108 for a, b in ((straight, cont), (straight.tangent, cont.tangent)):
109 for name in ("pressure_rayleigh", "pressure_unbaffled", "probe_acceleration"):
110 assert np.array_equal(getattr(a, name), getattr(b, name)), name
111 assert np.array_equal(a.displacement.w, b.displacement.w)
112 assert np.array_equal(straight.tangent.norms, cont.tangent.norms)
113 loaded = load_cymsim(tmp_path / "a.cymsim")
114 assert np.array_equal(loaded.tangent.pressure_unbaffled, straight.tangent.pressure_unbaffled)
115 assert np.array_equal(loaded.tangent.norms, straight.tangent.norms)
116 short = cut_run(loaded, _tangent_config(coarse_config, 1000.0, 0.015))
117 direct = load_cymsim(tmp_path / "b.cymsim")
118 for a, b in ((short, direct), (short.tangent, direct.tangent)):
119 assert np.array_equal(a.pressure_unbaffled, b.pressure_unbaffled)
120 assert np.array_equal(a.displacement.w, b.displacement.w)
121 assert np.array_equal(short.tangent.norms, direct.tangent.norms)
122 assert np.array_equal(short.tangent.state.q, loaded.tangent.state.q)