Tangent 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 runs - proof-of-concept/cymsim

proof-of-concept / cymsim

Tangent 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 runs

Jeremy Magland <jmagland@flatironinstitute.org> committed 3 days ago commit de5d1ce9412e parent 9558a4c Browse files

8 changed files+406−101

docs/format.mdmodified+5−1View file

@@ -129,7 +129,11 @@ where the fractional delay $r_e f_s / c$ in samples must be applied with a band-

129129

130130 The state needed to continue the run from `n_steps_done`. Its content depends on the program that wrote it, named in the group attribute `writer`; readers of the format may ignore the group, and a program should use only a state written by itself. The Python package (`writer = "cymsim-python"`, attributes `n_steps_done`, `dtype`) stores `q` $(N, 5)$ and `v` $(N, 5)$, the nodal displacements $(u, v, w, \theta_x, \theta_y)$ and velocities; `exciter_z`, `exciter_v`, `peak_force`, `contact_steps` per hit; `ring_front`, `ring_back`, the pending retarded-time pressure contributions at the sub-step rate; `decimator_tail_front`, `decimator_tail_back`, the sub-steps still needed by audio samples that are not yet final; and, when recording, `acceleration_ring` $(W, N)$, the partial sums of the not-yet-final acceleration samples; `decimator_tail_probes` $(P, L)$, the decimator tails of the probe signals; and the *overhang*, the samples past `duration` that were already final when the run stopped (the integrator overshoots the duration by $W/2$ audio samples and by the rounding to whole blocks): `overhang_front`, `overhang_unbaffled` (float64), `overhang_probes`, `overhang_acceleration_data`, `overhang_acceleration_scale`, `overhang_displacement_time` and `overhang_displacement_w`, each empty in a partial file. `cymsim run --resume` continues a run from this group, for the same configuration with a longer `duration` (or a shorter one, which cuts the records); the continued run is identical to an uninterrupted one on the same device. A state without the overhang datasets was written before runs could be continued and cannot be resumed. The web companion stores its own float32 GPU state under `writer = "cymsimweb"`.

131131

132-## 12. Reading the file elsewhere

132+## 12. Tangent: `tangent` (optional)

133+

134+Present for a run with `numerics.tangent: true`. The run then also integrates its tangent: the derivative of the state with respect to $\delta$ when every prescribed peak force is multiplied by $(1 + \delta)$, carried exactly by forward-mode differentiation of the time step (`jax.jvp`). It is the limit of $(\text{twin} - \text{run}) / \delta$ for a twin run whose peak forces are larger by the relative amount $\delta$, and does not saturate as a finite twin does. The group holds the tangent's outputs in the layout of the run's own: `audio/baffled`, `audio/unbaffled` (Pa per unit $\delta$), `probes/normal_acceleration` (when there are probes), `displacement/w` and `displacement/time` (when frames are kept; the rate is that of section 7); `norms` $(B, 4)$, one row after every block, with columns (attribute `columns`) `time tangent_norm2 kinetic linear_strain`: the squared energy norm of the tangent state, $T(\dot{q}') + U_\text{lin}(q')$, and the two parts of the same norm of the run's state, $T(\dot{q})$ and $U_\text{lin}(q)$, where $T$ is the kinetic energy and $U_\text{lin}$ the strain energy without the von Kármán terms (a quadratic form); and `state`, the tangent's integrator state in the layout of section 11, from which a tangent run continues. The relative separation of an infinitesimal twin is $\delta \sqrt{\texttt{tangent\_norm2} / (\texttt{kinetic} + \texttt{linear\_strain})}$.

135+

136+## 13. Reading the file elsewhere

133137

134138 Any HDF5 reader works. In Python without cymsim:

135139

src/cymsim/config.pymodified+5−0View file

@@ -139,6 +139,7 @@ class Numerics:

139139 steps_per_block: int = 20000 # sub-steps per compiled block

140140 dtype: str = "float32" # simulation precision: "float32" or "float64"

141141 nonlinear: bool = True # False drops the quadratic strain terms (linear reference run)

142+ tangent: bool = False # also integrate the tangent: d/d delta of the run when every prescribed peak force is x (1 + delta)

142143 n_progress: int = 20 # number of progress/energy reports during the run

143144

144145

@@ -186,6 +187,10 @@ class SimulationConfig:

186187 raise ValueError("numerics.cfl_safety must lie in (0, 1]")

187188 if self.mesh.element_size <= 0 or self.mesh.element_size > p.radius / 2:

188189 raise ValueError("mesh.element_size must be positive and well below the radius")

190+ if n.tangent and self.output.record_acceleration:

191+ raise ValueError("numerics.tangent does not support output.record_acceleration")

192+ if n.tangent and not any(isinstance(h.exciter, ForceExciter) for h in self.hits):

193+ raise ValueError("numerics.tangent differentiates the prescribed peak forces: it needs a force exciter")

189194 if self.output.displacement_rate < 0:

190195 raise ValueError("output.displacement_rate must be non-negative")

191196 for i, h in enumerate(self.hits):

src/cymsim/cymsim_file.pymodified+72−21View file

@@ -23,8 +23,8 @@ import yaml

2323 from .config import SimulationConfig

2424 from .decimation import CUTOFF, DECIMATION_SPAN, KAISER_BETA

2525 from .mesh import Mesh

26-from .simulate import (AccelerationRecord, DisplacementRecord, EnergyRecord, HitSummary,

27- SimulationResult, SimulatorState)

26+from .simulate import (NORM_COLUMNS, AccelerationRecord, DisplacementRecord, EnergyRecord, HitSummary,

27+ SimulationResult, SimulatorState, TangentRecord)

2828

2929 FORMAT_VERSION = 1

3030 ENERGY_COLUMNS = "time kinetic membrane bending exciters max_w_over_h"

@@ -129,6 +129,24 @@ class CymsimWriter:

129129 g.create_dataset("normal_acceleration", shape=(0, len(r.probes)), maxshape=(None, len(r.probes)),

130130 dtype=sig, chunks=(4096, len(r.probes)))

131131

132+ if r.tangent is not None:

133+ g = f.create_group("tangent")

134+ g.attrs["description"] = ("derivative of the run with respect to delta when every prescribed peak force "

135+ "is multiplied by (1 + delta): the limit of (twin - run) / delta")

136+ for name in ("baffled", "unbaffled"):

137+ g.create_dataset(f"audio/{name}", shape=(0,), maxshape=(None,), dtype=sig, chunks=(65536,))

138+ if r.probes is not None:

139+ g.create_dataset("probes/normal_acceleration", shape=(0, len(r.probes)), maxshape=(None, len(r.probes)),

140+ dtype=sig, chunks=(4096, len(r.probes)))

141+ if r.tangent.displacement is not None:

142+ g.create_dataset("displacement/w", shape=(0, N), maxshape=(None, N), dtype=sig, chunks=(64, N))

143+ g.create_dataset("displacement/time", shape=(0,), maxshape=(None,), dtype=np.float64, chunks=(4096,))

144+ d = g.create_dataset("norms", shape=(0, 4), maxshape=(None, 4), dtype=np.float64, chunks=(1024, 4))

145+ d.attrs["columns"] = NORM_COLUMNS

146+ d.attrs["description"] = ("after every block: time; kinetic_energy(dv) + linear_energy(dq) of the tangent; "

147+ "kinetic_energy(v) and linear_energy(q) of the run (linear_energy: strain energy "

148+ "without the von Karman terms)")

149+

132150 g = f.create_group("energy")

133151 d = g.create_dataset("table", shape=(0, 6), maxshape=(None, 6), dtype=np.float64, chunks=(1024, 6))

134152 d.attrs["columns"] = ENERGY_COLUMNS

@@ -148,6 +166,23 @@ class CymsimWriter:

148166 ds.resize(len(data), axis=0)

149167 ds[n0:] = data[n0:]

150168

169+ @staticmethod

170+ def _write_state(parent, name: str, s: SimulatorState | None, dtype: str) -> None:

171+ """(Re)write the integrator state ``s`` as the group parent/name."""

172+ if name in parent:

173+ del parent[name]

174+ if s is None:

175+ return

176+ g = parent.create_group(name)

177+ g.attrs["writer"] = STATE_WRITER

178+ g.attrs["n_steps_done"] = int(s.n_steps_done)

179+ g.attrs["dtype"] = dtype

180+ for key in STATE_ARRAYS:

181+ g.create_dataset(key, data=np.asarray(getattr(s, key)))

182+ for key in OPTIONAL_STATE_ARRAYS:

183+ if getattr(s, key) is not None:

184+ g.create_dataset(key, data=np.asarray(getattr(s, key)))

185+

151186 def update(self, result: SimulationResult) -> None:

152187 f, r = self.f, result

153188 self._append(f["audio/baffled"], r.pressure_rayleigh.astype(self.signal_dtype))

@@ -168,19 +203,21 @@ class CymsimWriter:

168203 for name in HIT_FIELDS:

169204 f["hits"][name][...] = np.array([getattr(h, name) for h in r.hits], np.float64)

170205 f["hits/kind"][...] = np.array([h.kind.encode() for h in r.hits], "S16")

171- if "state" in f:

172- del f["state"]

173- if r.state is not None:

174- s = r.state

175- g = f.create_group("state")

176- g.attrs["writer"] = STATE_WRITER

177- g.attrs["n_steps_done"] = int(s.n_steps_done)

178- g.attrs["dtype"] = r.config.numerics.dtype

179- for name in STATE_ARRAYS:

180- g.create_dataset(name, data=np.asarray(getattr(s, name)))

181- for name in OPTIONAL_STATE_ARRAYS:

182- if getattr(s, name) is not None:

183- g.create_dataset(name, data=np.asarray(getattr(s, name)))

206+ self._write_state(f, "state", r.state, r.config.numerics.dtype)

207+ t = r.tangent

208+ if t is not None and "tangent" in f:

209+ g = f["tangent"]

210+ self._append(g["audio/baffled"], t.pressure_rayleigh.astype(self.signal_dtype))

211+ self._append(g["audio/unbaffled"], t.pressure_unbaffled.astype(self.signal_dtype))

212+ if t.probe_acceleration is not None and "probes" in g:

213+ self._append(g["probes/normal_acceleration"], t.probe_acceleration.astype(self.signal_dtype))

214+ if t.displacement is not None and "displacement" in g:

215+ self._append(g["displacement/w"], t.displacement.w.astype(self.signal_dtype))

216+ self._append(g["displacement/time"], t.displacement.time)

217+ g["norms"].resize(len(t.norms), axis=0)

218+ if len(t.norms):

219+ g["norms"][...] = t.norms

220+ self._write_state(g, "state", t.state, r.config.numerics.dtype)

184221 f["run"].attrs["n_steps_done"] = int(r.n_steps_done)

185222 f["run"].attrs["complete"] = np.int32(1 if r.complete else 0) # an integer flag, readable by every HDF5 client

186223 f["run"].attrs["wall_time"] = float(r.wall_time)

@@ -229,11 +266,17 @@ def load_cymsim(path: str | Path) -> SimulationResult:

229266 hits = [HitSummary(float(hg["time"][k]), kinds[k], float(hg["peak_force"][k]), float(hg["contact_time"][k]),

230267 float(hg["rebound_velocity"][k]), float(hg["energy_in"][k]), float(hg["energy_out"][k]))

231268 for k in range(len(kinds))]

232- state = None

233- if "state" in f and f["state"].attrs.get("writer", "") == STATE_WRITER:

234- g = f["state"]

235- state = SimulatorState(int(g.attrs["n_steps_done"]), *(g[name][...] for name in STATE_ARRAYS),

236- **{name: g[name][...] for name in OPTIONAL_STATE_ARRAYS if name in g})

269+ state = _read_state(f, "state")

270+ tangent = None

271+ if "tangent" in f:

272+ g = f["tangent"]

273+ tdisp = None

274+ if "displacement" in g:

275+ tdisp = DisplacementRecord(float(f["displacement"].attrs["rate"]), g["displacement/time"][...],

276+ g["displacement/w"][...])

277+ tangent = TangentRecord(g["audio/baffled"][...].astype(np.float64), g["audio/unbaffled"][...].astype(np.float64),

278+ g["probes/normal_acceleration"][...].astype(np.float64) if "probes" in g else None,

279+ tdisp, g["norms"][...], _read_state(g, "state"))

237280 probes = probe_acc = None

238281 if "probes" in f:

239282 probes = f["probes/xy"][...]

@@ -246,7 +289,15 @@ def load_cymsim(path: str | Path) -> SimulationResult:

246289 f["audio/unbaffled"][...].astype(np.float64), hits, energy_log, acc, disp,

247290 bool(run.get("complete", False)), int(run.get("n_steps_done", 0)), state,

248291 float(run.get("wall_time", 0.0)), creator.decode() if isinstance(creator, bytes) else str(creator),

249- probes, probe_acc)

292+ probes, probe_acc, tangent)

293+

294+

295+def _read_state(parent, name: str) -> SimulatorState | None:

296+ if name not in parent or parent[name].attrs.get("writer", "") != STATE_WRITER:

297+ return None

298+ g = parent[name]

299+ return SimulatorState(int(g.attrs["n_steps_done"]), *(g[key][...] for key in STATE_ARRAYS),

300+ **{key: g[key][...] for key in OPTIONAL_STATE_ARRAYS if key in g})

250301

251302

252303 def describe(path: str | Path) -> str:

src/cymsim/exciter.pymodified+5−3View file

@@ -143,12 +143,14 @@ class HitSchedule:

143143 """(face height, velocity) of every exciter before release."""

144144 return (jnp.full(self.n_hits, HOLD_HEIGHT, self.dtype), jnp.zeros(self.n_hits, self.dtype))

145145

146- def contact(self, n, q, zs, vs):

146+ def contact(self, n, q, zs, vs, f_peak=None):

147147 """Contact forces at step n.

148148

149149 Returns the per-point forces Fk (K, P), the total force per hit F (K,), and the exciter

150- state (zs, vs) after applying the release of any hit whose time is step n.

150+ state (zs, vs) after applying the release of any hit whose time is step n. ``f_peak``

151+ replaces the prescribed peak forces (K,), so that they can be differentiated (tangent runs).

151152 """

153+ f_peak = self.f_peak if f_peak is None else f_peak

152154 wp = jnp.einsum("kpi,kpi->kp", self.weights, q[self.nodes, 2]) + self.offsets

153155 release = n == self.step

154156 active = n >= self.step

@@ -159,7 +161,7 @@ class HitSchedule:

159161 # prescribed pulses (k_point is 0 for these, so the contact term above vanishes)

160162 t = (n - self.step).astype(self.dtype) * self.dt

161163 on = self.is_force & (t >= 0) & (t <= self.f_dur)

162- pulse = jnp.where(on, 0.5 * self.f_peak * (1 - jnp.cos(2 * jnp.pi * t / self.f_dur)), 0.0)

164+ pulse = jnp.where(on, 0.5 * f_peak * (1 - jnp.cos(2 * jnp.pi * t / self.f_dur)), 0.0)

163165 Fk = Fk + pulse[:, None] * self.first_point

164166 return Fk, Fk.sum(1), zs, vs

165167

src/cymsim/shell.pymodified+7−0View file

@@ -188,6 +188,13 @@ class ShellModel:

188188 q[self.tri_j], self.Gx_j, self.Gy_j, self.area_j, self.z0x_j, self.z0y_j,

189189 self.Cm_j, self.Ke_j, self.nl)

190190

191+ def linear_energy(self, q):

192+ """Strain energy of q with the quadratic (von Karman) strain terms dropped: a quadratic

193+ form, so that sqrt(kinetic_energy(v) + linear_energy(q)) is a norm on the state."""

194+ Em, Eb = jax.vmap(element_energy_parts, in_axes=(0, 0, 0, 0, 0, 0, 0, 0, None))(

195+ q[self.tri_j], self.Gx_j, self.Gy_j, self.area_j, self.z0x_j, self.z0y_j, self.Cm_j, self.Ke_j, 0.0)

196+ return Em.sum() + Eb.sum()

197+

191198 def total_energy(self, q):

192199 Em, Eb = self.energy_parts_per_element(q)

193200 return Em.sum() + Eb.sum()

src/cymsim/simulate.pymodified+198−75View file

@@ -161,6 +161,24 @@ class SimulatorState:

161161 overhang_displacement_w: np.ndarray | None = None # (Td', N)

162162

163163

164+NORM_COLUMNS = "time tangent_norm2 kinetic linear_strain"

165+

166+

167+@dataclass

168+class TangentRecord:

169+ """The tangent of a run (numerics.tangent): the derivative of its outputs and state with

170+ respect to delta when every prescribed peak force is multiplied by (1 + delta), the limit of

171+ (twin - run) / delta for twin runs. Signals as in SimulationResult; ``norms`` holds, after

172+ every block, the time, the squared energy norm of the tangent state, kinetic_energy(dv) +

173+ linear_energy(dq), and the two parts of the same norm of the run's own state."""

174+ pressure_rayleigh: np.ndarray

175+ pressure_unbaffled: np.ndarray

176+ probe_acceleration: np.ndarray | None

177+ displacement: DisplacementRecord | None

178+ norms: np.ndarray # (B, 4), columns NORM_COLUMNS

179+ state: SimulatorState | None = None

180+

181+

164182 @dataclass

165183 class SimulationResult:

166184 config: SimulationConfig

@@ -186,6 +204,7 @@ class SimulationResult:

186204 creator: str = ""

187205 probes: np.ndarray | None = None # (P, 2) probe points [m]

188206 probe_acceleration: np.ndarray | None = None # (T, P) normal acceleration at the probes [m/s^2]

207+ tangent: TangentRecord | None = None # with numerics.tangent

189208

190209 @property

191210 def n_nodes(self) -> int:

@@ -243,6 +262,24 @@ def _log_until(energy_log: list[EnergyRecord], t_end: float) -> list[EnergyRecor

243262 return [e for e in energy_log if e.time <= t_end * (1 + 1e-12)]

244263

245264

265+def _cut_signals(front, unb, probe, disp, s: SimulatorState, n_out: int, dt: float, substeps: int):

266+ """Signals cut to n_out samples (the state's overhang appended first), and the overhang

267+ fields of the state past n_out."""

268+ cat = lambda a, b: np.concatenate([a, b])

269+ front, unb = cat(front, s.overhang_front), cat(unb, s.overhang_unbaffled)

270+ over = dict(overhang_front=front[n_out:], overhang_unbaffled=unb[n_out:])

271+ if probe is not None:

272+ probe = cat(probe, s.overhang_probes)

273+ over["overhang_probes"] = probe[n_out:]

274+ probe = probe[:n_out]

275+ if disp is not None:

276+ t, w = cat(disp.time, s.overhang_displacement_time), cat(disp.w, s.overhang_displacement_w)

277+ k_out = int(np.count_nonzero(np.round(t / dt) <= n_out * substeps)) # frame times are steps * dt

278+ over.update(overhang_displacement_time=t[k_out:], overhang_displacement_w=w[k_out:])

279+ disp = DisplacementRecord(disp.rate, t[:k_out], w[:k_out])

280+ return front[:n_out], unb[:n_out], probe, disp, over

281+

282+

246283 def cut_run(prev: SimulationResult, config: SimulationConfig) -> SimulationResult | None:

247284 """The run of ``config`` taken from the longer run ``prev`` (whose configuration differs at

248285 most in CONTINUATION_FIELDS), or None if ``prev`` does not reach the end of ``config``'s run.

@@ -255,29 +292,98 @@ def cut_run(prev: SimulationResult, config: SimulationConfig) -> SimulationResul

255292 n_steps = steps_for(n_out, prev.substeps, num.steps_per_block)

256293 if s is None or s.overhang_front is None or s.n_steps_done < n_steps:

257294 return None

258- cat = lambda a, b: np.concatenate([a, b])

259- front, unb = cat(prev.pressure_rayleigh, s.overhang_front), cat(prev.pressure_unbaffled, s.overhang_unbaffled)

260- over = dict(overhang_front=front[n_out:], overhang_unbaffled=unb[n_out:])

261- probe = None

262- if prev.probe_acceleration is not None:

263- probe = cat(prev.probe_acceleration, s.overhang_probes)

264- over["overhang_probes"] = probe[n_out:]

295+ front, unb, probe, disp, over = _cut_signals(prev.pressure_rayleigh, prev.pressure_unbaffled,

296+ prev.probe_acceleration, prev.displacement, s, n_out,

297+ prev.dt, prev.substeps)

265298 acc = prev.acceleration

266299 if acc is not None:

300+ cat = lambda a, b: np.concatenate([a, b])

267301 data, scale = cat(acc.data, s.overhang_acceleration_data), cat(acc.scale, s.overhang_acceleration_scale)

268302 over.update(overhang_acceleration_data=data[n_out:], overhang_acceleration_scale=scale[n_out:])

269303 acc = AccelerationRecord(acc.fs, data[:n_out], scale[:n_out], acc.taps)

270- disp = prev.displacement

271- if disp is not None:

272- t, w = cat(disp.time, s.overhang_displacement_time), cat(disp.w, s.overhang_displacement_w)

273- k_out = int(np.count_nonzero(np.round(t / prev.dt) <= n_out * prev.substeps)) # frame times are steps * dt

274- over.update(overhang_displacement_time=t[k_out:], overhang_displacement_w=w[k_out:])

275- disp = DisplacementRecord(disp.rate, t[:k_out], w[:k_out])

304+ tangent = prev.tangent

305+ if tangent is not None:

306+ ts = tangent.state

307+ tf, tu, tp, td, tover = _cut_signals(tangent.pressure_rayleigh, tangent.pressure_unbaffled,

308+ tangent.probe_acceleration, tangent.displacement, ts, n_out,

309+ prev.dt, prev.substeps)

310+ norms = tangent.norms[tangent.norms[:, 0] <= n_steps * prev.dt * (1 + 1e-12)]

311+ tangent = TangentRecord(tf, tu, tp, td, norms, dataclasses.replace(ts, **tover))

276312 return dataclasses.replace(

277- prev, config=config, pressure_rayleigh=front[:n_out], pressure_unbaffled=unb[:n_out], acceleration=acc,

313+ prev, config=config, pressure_rayleigh=front, pressure_unbaffled=unb, acceleration=acc,

278314 displacement=disp, complete=True, state=dataclasses.replace(s, **over),

279- energy_log=_log_until(prev.energy_log, n_steps * prev.dt),

280- probe_acceleration=None if probe is None else probe[:n_out])

315+ energy_log=_log_until(prev.energy_log, n_steps * prev.dt), probe_acceleration=probe, tangent=tangent)

316+

317+

318+class _Stream:

319+ """What a state tuple streams out block by block (the run's, or its tangent's): the listener

320+ pressure and the probe signals through their decimators, and the displacement frames."""

321+

322+ def __init__(self, sim: "Simulator", prev=None):

323+ """``prev``: (pressure_rayleigh, pressure_unbaffled, probe_acceleration, displacement,

324+ state) of the run being continued, or None for a run from rest."""

325+ m, P = sim.substeps, len(sim.probe_xy)

326+ self.sim, self.P = sim, P

327+ self.disp = (DisplacementRecord.empty(sim.disp_rate, sim.mesh.n_nodes, np.dtype(sim.shell.dtype))

328+ if sim.disp_rate > 0 else None)

329+ if prev is None:

330+ self.f, self.u, self.p = [np.zeros(0)], [np.zeros(0)], [np.zeros((0, P))]

331+ self.dec_f, self.dec_b = StreamingDecimator(m), StreamingDecimator(m)

332+ self.dec_p = [StreamingDecimator(m) for _ in range(P)]

333+ return

334+ front, unb, probe, disp, s = prev

335+ n = s.n_steps_done

336+ self.f, self.u = [front, s.overhang_front], [unb, s.overhang_unbaffled]

337+ self.p = [probe, s.overhang_probes] if P else [np.zeros((0, 0))]

338+ self.dec_f = StreamingDecimator(m, tail=s.decimator_tail_front, n_seen=n)

339+ self.dec_b = StreamingDecimator(m, tail=s.decimator_tail_back, n_seen=n)

340+ self.dec_p = [StreamingDecimator(m, tail=t, n_seen=n) for t in s.decimator_tail_probes[:P]] if P else []

341+ if self.disp is not None:

342+ self.disp.extend(disp.time, disp.w)

343+ self.disp.extend(s.overhang_displacement_time, s.overhang_displacement_w)

344+

345+ def take(self, state, n0: int, n1: int, k0: int, k1: int):

346+ """Collect the outputs of the block [n0, n1) from ``state``; returns the state with the

347+ pressure ring slots it read cleared."""

348+ sim = self.sim

349+ # the ring slots of this block's sub-steps are final; read, clear, decimate

350+ slots = jnp.asarray(sim.radiation.slots(n0, n1))

351+ pf, pb = np.asarray(state[4][slots], np.float64), np.asarray(state[5][slots], np.float64)

352+ state = state[:4] + (sim._clear(state[4], slots), sim._clear(state[5], slots)) + state[6:]

353+ # (samples past the duration are kept too; _result moves them to the state's overhang)

354+ new_f, new_b = self.dec_f.push(pf), self.dec_b.push(pb)

355+ self.f.append(new_f)

356+ self.u.append(new_f - new_b)

357+ if self.P:

358+ rows = np.asarray(state[11], np.float64)

359+ self.p.append(np.stack([d.push(rows[:, i]) for i, d in enumerate(self.dec_p)], 1))

360+ if self.disp is not None and k1 > k0:

361+ frames = np.asarray(state[10][:k1 - k0], np.dtype(sim.shell.dtype))

362+ self.disp.extend(sim.frame_steps[k0:k1] * sim.dt, frames)

363+ return state

364+

365+ def split(self, n_done: int, state, record_ring=None):

366+ """(front, unbaffled, probe, displacement) up to the duration and the SimulatorState

367+ after n_done steps, which keeps everything past the duration as its overhang."""

368+ sim = self.sim

369+ q, v, zs, vs, bf, bb, Fmax, ncont = jax.device_get(state[:8])

370+ front, unb = np.concatenate(self.f), np.concatenate(self.u)

371+ n = min(len(front), len(unb))

372+ probe = np.concatenate(self.p) if self.P else None

373+ empty = SimulatorState(0, *(None,) * 10, overhang_front=np.zeros(0), overhang_unbaffled=np.zeros(0),

374+ overhang_probes=np.zeros((0, self.P)) if self.P else None,

375+ overhang_displacement_time=np.zeros(0) if self.disp is not None else None,

376+ overhang_displacement_w=(np.zeros((0, sim.mesh.n_nodes), self.disp.w.dtype)

377+ if self.disp is not None else None))

378+ front, unb, probe, disp, over = _cut_signals(front[:n], unb[:n], probe, self.disp, empty, sim.n_out,

379+ sim.dt, sim.substeps)

380+ ncont = np.asarray(ncont) if np.issubdtype(np.asarray(ncont).dtype, np.integer) \

381+ else np.zeros(np.shape(ncont), np.int32) # a tangent's float0 counts

382+ st = SimulatorState(

383+ n_done, np.asarray(q), np.asarray(v), np.asarray(zs), np.asarray(vs), np.asarray(Fmax), ncont,

384+ np.asarray(bf), np.asarray(bb), self.dec_f.tail.copy(), self.dec_b.tail.copy(),

385+ record_ring, np.stack([d.tail for d in self.dec_p]) if self.dec_p else None, **over)

386+ return front, unb, probe, disp, st

281387

282388

283389 class Simulator:

@@ -387,7 +493,7 @@ class Simulator:

387493 n_probes = len(self.probe_xy)

388494 p_nodes, p_w = self.probe_nodes, self.probe_w

389495

390- def step(n, k_block0, j_block0, st):

496+ def step(n, k_block0, j_block0, st, f_peak=None):

391497 q, v, zs, vs, bf, bb, Fmax, ncont, ring, acc_out, disp_out, pbuf = st

392498 if rate > 0:

393499 # frame k sits at sub-step round(k S / rate); emit q at step n if n is one of them

@@ -395,7 +501,7 @@ class Simulator:

395501 is_frame = jnp.round(k_cand * (S / rate)) == n

396502 row = jnp.clip(k_cand.astype(jnp.int32) - k_block0, 0, self.disp_rows - 1)

397503 disp_out = disp_out.at[row].add(jnp.where(is_frame, q[:, 2], 0.0))

398- Fk, F, zs, vs = sched.contact(n, q, zs, vs)

504+ Fk, F, zs, vs = sched.contact(n, q, zs, vs, f_peak)

399505 fo = sh.internal_force(q)

400506 fo = sched.apply_to_plate(fo, Fk)

401507 fo = fo.at[:, 2].add(-sh.damping_force_w(v[:, 2]))

@@ -428,6 +534,18 @@ class Simulator:

428534 def run_block(n0, k_block0, j_block0, st):

429535 return lax.fori_loop(0, spb, lambda k, s: step((n0 + k).astype(jnp.int32), k_block0, j_block0, s), st)

430536

537+ @jax.jit

538+ def tangent_block(n0, k_block0, j_block0, st, dst):

539+ """One block of the run together with its tangent dst: the derivative of the state

540+ with respect to delta when every prescribed peak force is multiplied by (1 + delta),

541+ carried exactly by forward-mode differentiation of the step."""

542+ fp = sched.f_peak

543+

544+ def body(k, carry):

545+ n = (n0 + k).astype(jnp.int32)

546+ return jax.jvp(lambda s, f: step(n, k_block0, j_block0, s, f), (carry[0], fp), (carry[1], fp))

547+ return lax.fori_loop(0, spb, body, (st, dst))

548+

431549 @jax.jit

432550 def energies(st):

433551 q, v, _, vs = st[0], st[1], st[2], st[3]

@@ -439,7 +557,21 @@ class Simulator:

439557 def clear(buf, idx):

440558 return buf.at[idx].set(0.0)

441559

560+ @jax.jit

561+ def norms(st, dst):

562+ """Squared energy norm of the tangent state, and the two parts of the run's own."""

563+ return (sh.kinetic_energy(dst[1]) + sh.linear_energy(dst[0]),

564+ sh.kinetic_energy(st[1]), sh.linear_energy(st[0]))

565+

442566 self._step, self._run_block, self._energies, self._clear = step, run_block, energies, clear

567+ self._tangent_block, self._norms = tangent_block, norms

568+

569+ @staticmethod

570+ def zero_tangent(state):

571+ """A zero tangent for a device state tuple (float0 for its integer arrays)."""

572+ return jax.tree_util.tree_map(

573+ lambda a: np.zeros(a.shape, jax.dtypes.float0) if jnp.issubdtype(a.dtype, jnp.integer)

574+ else jnp.zeros_like(a), state)

443575

444576 def initial_state(self):

445577 N, dtype = self.mesh.n_nodes, self.shell.dtype

@@ -461,6 +593,11 @@ class Simulator:

461593 f(s.peak_force), jnp.asarray(s.contact_steps, jnp.int32), ring, acc_out, disp_out,

462594 jnp.zeros((self.config.numerics.steps_per_block, len(self.probe_xy)), dtype))

463595

596+ def _tangent_from(self, s: SimulatorState):

597+ """The device tangent tuple rebuilt from a saved tangent SimulatorState."""

598+ st = self._state_from(s)

599+ return st[:7] + (np.zeros(st[7].shape, jax.dtypes.float0),) + st[8:]

600+

464601 def check_resumable(self, prev: SimulationResult) -> None:

465602 """Raise ValueError unless this simulator can continue the run ``prev``."""

466603 s = prev.state

@@ -478,6 +615,8 @@ class Simulator:

478615 raise ValueError(f"the earlier run's state is {s.q.dtype}, this run is {np.dtype(self.shell.dtype)}")

479616 if s.n_steps_done % self.config.numerics.steps_per_block:

480617 raise ValueError("the earlier run stopped inside a block")

618+ if self.config.numerics.tangent and (prev.tangent is None or prev.tangent.state is None):

619+ raise ValueError("the earlier run has no tangent state")

481620

482621 def _empty_outputs(self):

483622 N, dtype = self.mesh.n_nodes, self.shell.dtype

@@ -491,36 +630,29 @@ class Simulator:

491630 float(vs[k]), p.energy, 0.5 * p.mass * float(vs[k]) ** 2))

492631 return hits

493632

494- def _result(self, state, energy_log, p_front, p_unb, acc, disp, n_done, wall, dec_f, dec_b, dec_p,

495- probe) -> SimulationResult:

633+ def _result(self, state, energy_log, stream: _Stream, acc, n_done, wall, tangent=None) -> SimulationResult:

496634 """The result so far. Every series holds all final samples; those past the duration are

497- moved to the state's overhang."""

498- q, v, zs, vs, bf, bb, Fmax, ncont, ring = jax.device_get(state[:9])

635+ moved to the state's overhang. ``tangent``: (tangent state, its _Stream, norm rows)."""

499636 n_out = self.n_out

500- n = min(len(p_front), len(p_unb))

501- p_front, p_unb = np.asarray(p_front[:n]), np.asarray(p_unb[:n])

502- acc_over = (None, None)

637+ Fmax, ncont, vs = jax.device_get((state[6], state[7], state[3]))

638+ acc_over = {}

503639 if acc is not None:

504- acc_over = (acc.data[n_out:], acc.scale[n_out:])

640+ acc_over = dict(overhang_acceleration_data=acc.data[n_out:], overhang_acceleration_scale=acc.scale[n_out:])

505641 acc = AccelerationRecord(acc.fs, acc.data[:n_out], acc.scale[:n_out], acc.taps)

506- disp_over = (None, None)

507- if disp is not None:

508- k_out = int(np.searchsorted(self.frame_steps, n_out * self.substeps, side="right"))

509- disp_over = (disp.time[k_out:], disp.w[k_out:])

510- disp = DisplacementRecord(disp.rate, disp.time[:k_out], disp.w[:k_out])

511- sim_state = SimulatorState(

512- n_done, np.asarray(q), np.asarray(v), np.asarray(zs), np.asarray(vs), np.asarray(Fmax),

513- np.asarray(ncont), np.asarray(bf), np.asarray(bb), dec_f.tail.copy(), dec_b.tail.copy(),

514- np.asarray(ring) if self.record else None,

515- np.stack([d.tail for d in dec_p]) if dec_p else None,

516- p_front[n_out:], p_unb[n_out:], None if probe is None else probe[n_out:], *acc_over, *disp_over)

642+ ring = np.asarray(jax.device_get(state[8])) if self.record else None

643+ p_front, p_unb, probe, disp, sim_state = stream.split(n_done, state, ring)

644+ sim_state = dataclasses.replace(sim_state, **acc_over)

645+ trec = None

646+ if tangent is not None:

647+ dstate, tstream, norms = tangent

648+ tf, tu, tp, td, tst = tstream.split(n_done, dstate)

649+ trec = TangentRecord(tf, tu, tp, td, np.array(norms, np.float64).reshape(-1, 4), tst)

517650 from . import __version__

518651 return SimulationResult(

519652 self.config, self.fs, self.dt, self.substeps, self.mesh, self.shell.z0_node, self.shell.h_node,

520- self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front[:n_out], p_unb[:n_out],

653+ self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front, p_unb,

521654 self._hit_summaries(Fmax, ncont, vs), list(energy_log), acc, disp, n_done >= self.n_steps, n_done,

522- sim_state, wall, f"cymsim {__version__}",

523- self.probe_xy if len(self.probe_xy) else None, None if probe is None else probe[:n_out])

655+ sim_state, wall, f"cymsim {__version__}", self.probe_xy if len(self.probe_xy) else None, probe, trec)

524656

525657 def run(self, save_path: str | None = None, save_every: float | None = None,

526658 resume: SimulationResult | None = None, checkpoint=None,

@@ -538,39 +670,38 @@ class Simulator:

538670 covers the duration, nothing is integrated and its records are cut to the duration."""

539671 num = self.config.numerics

540672 m, spb, N = self.substeps, num.steps_per_block, self.mesh.n_nodes

541- n_probes = len(self.probe_xy)

542673 acc = AccelerationRecord.empty(self.fs, N, self.taps) if self.record else None

543- disp = DisplacementRecord.empty(self.disp_rate, N, np.dtype(self.shell.dtype)) if self.disp_rate > 0 else None

674+ tangent = bool(num.tangent)

675+ dstate = tstream = norms = None

544676 if resume is None:

545677 state, n_done, wall0 = self.initial_state(), 0, 0.0

546678 energy_log: list[EnergyRecord] = []

547- audio_f, audio_u, audio_p = [np.zeros(0)], [np.zeros(0)], [np.zeros((0, n_probes))]

548- dec_f, dec_b = StreamingDecimator(m), StreamingDecimator(m)

549- dec_p = [StreamingDecimator(m) for _ in range(n_probes)]

679+ stream = _Stream(self)

680+ if tangent:

681+ dstate, tstream, norms = self.zero_tangent(state), _Stream(self), []

550682 else:

551683 self.check_resumable(resume)

552684 s = resume.state

553685 state, n_done, wall0 = self._state_from(s), s.n_steps_done, resume.wall_time

554686 energy_log = list(resume.energy_log)

555- audio_f = [resume.pressure_rayleigh, s.overhang_front]

556- audio_u = [resume.pressure_unbaffled, s.overhang_unbaffled]

557- audio_p = [resume.probe_acceleration, s.overhang_probes] if n_probes else [np.zeros((0, 0))]

558- dec_f = StreamingDecimator(m, tail=s.decimator_tail_front, n_seen=n_done)

559- dec_b = StreamingDecimator(m, tail=s.decimator_tail_back, n_seen=n_done)

560- dec_p = [StreamingDecimator(m, tail=t, n_seen=n_done) for t in s.decimator_tail_probes[:n_probes]] \

561- if n_probes else []

687+ stream = _Stream(self, (resume.pressure_rayleigh, resume.pressure_unbaffled, resume.probe_acceleration,

688+ resume.displacement, s))

689+ if tangent:

690+ t = resume.tangent

691+ dstate = self._tangent_from(t.state)

692+ tstream = _Stream(self, (t.pressure_rayleigh, t.pressure_unbaffled, t.probe_acceleration,

693+ t.displacement, t.state))

694+ norms = [tuple(r) for r in t.norms]

562695 if acc is not None:

563696 acc.extend(resume.acceleration.data, resume.acceleration.scale)

564697 acc.extend(s.overhang_acceleration_data, s.overhang_acceleration_scale)

565- if disp is not None:

566- disp.extend(resume.displacement.time, resume.displacement.w)

567- disp.extend(s.overhang_displacement_time, s.overhang_displacement_w)

568698 if n_done >= self.n_steps:

569699 energy_log = _log_until(energy_log, self.n_steps * self.dt)

700+ if tangent:

701+ norms = [r for r in norms if r[0] <= self.n_steps * self.dt * (1 + 1e-12)]

570702 log.info("continuing the run from t = %.3f s", n_done * self.dt)

571- probe_now = lambda: np.concatenate(audio_p) if n_probes else None

572- result_now = lambda n, wall: self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_u),

573- acc, disp, n, wall, dec_f, dec_b, dec_p, probe_now())

703+ result_now = lambda n, wall: self._result(state, energy_log, stream, acc, n, wall,

704+ (dstate, tstream, norms) if tangent else None)

574705 writer = None

575706 if save_path is not None:

576707 from .cymsim_file import CymsimWriter

@@ -588,22 +719,17 @@ class Simulator:

588719 j_block0 = -(-n0 // m) - self.span // 2

589720 acc_out, disp_out = self._empty_outputs()

590721 state = state[:9] + (acc_out, disp_out) + state[11:]

591- state = self._run_block(jnp.int32(n0), jnp.int32(k0), jnp.int32(j_block0), state)

592-

593- # ---- pressure: the ring slots of this block's sub-steps are final; read, clear, decimate

594- slots = jnp.asarray(self.radiation.slots(n0, n1))

595- pf, pb = np.asarray(state[4][slots], np.float64), np.asarray(state[5][slots], np.float64)

596- state = state[:4] + (self._clear(state[4], slots), self._clear(state[5], slots)) + state[6:]

597- # (samples past the duration are kept too; _result moves them to the state's overhang)

598- new_f, new_b = dec_f.push(pf), dec_b.push(pb)

599- audio_f.append(new_f)

600- audio_u.append(new_f - new_b)

601- if n_probes:

602- rows = np.asarray(state[11], np.float64)

603- audio_p.append(np.stack([d.push(rows[:, i]) for i, d in enumerate(dec_p)], 1))

722+ if tangent:

723+ dstate = dstate[:9] + self._empty_outputs() + dstate[11:]

724+ state, dstate = self._tangent_block(jnp.int32(n0), jnp.int32(k0), jnp.int32(j_block0), state, dstate)

725+ dstate = tstream.take(dstate, n0, n1, k0, k1)

726+ norms.append((n1 * self.dt, *(float(x) for x in self._norms(state, dstate))))

727+ else:

728+ state = self._run_block(jnp.int32(n0), jnp.int32(k0), jnp.int32(j_block0), state)

729+ state = stream.take(state, n0, n1, k0, k1)

604730 n_done = n1

605731

606- # ---- records

732+ # ---- the acceleration record

607733 if self.record:

608734 n_emit = -(-n1 // m) - -(-n0 // m)

609735 rows = np.asarray(state[9][:n_emit], np.float32)

@@ -611,9 +737,6 @@ class Simulator:

611737 sel = j >= 0

612738 if sel.any():

613739 acc.extend(*AccelerationRecord.encode(rows[sel]))

614- if self.disp_rate > 0 and k1 > k0:

615- frames = np.asarray(state[10][:k1 - k0], np.dtype(self.shell.dtype))

616- disp.extend(self.frame_steps[k0:k1] * self.dt, frames)

617740

618741 if blk % report_every == 0 or blk == self.n_blocks - 1:

619742 ke, em, eb, ks, woh = (float(x) for x in self._energies(state))

tests/test_cache.pymodified+1−1View file

@@ -14,7 +14,7 @@ from cymsim.config import ForceExciter, Hit, MalletExciter, SimulationConfig, St

1414 # sha256 of the full default configuration. When this test fails, a default has changed or a field

1515 # was added. If the results of some existing configuration change, bump MODEL_VERSION in cache.py;

1616 # either way, update the fingerprint.

17-DEFAULTS_FINGERPRINT = "4b9d969c01b47ddabc6ee40ac7315d0c1ed1ac998f4a95b120ed4e668268672c"

17+DEFAULTS_FINGERPRINT = "55d1085684bbe86effb8ac29dd02132d69c19611831cbc03dd59ea7645ec134e"

1818

1919

2020 def test_key_format():

tests/test_tangent.pyadded+113−0View file

@@ -0,0 +1,113 @@

1+"""The tangent of a run (Simulator._tangent_block) is the derivative of the state with respect to a

2+relative change of the prescribed peak force: it matches central differences of two twin runs."""

3+import copy

4+

5+import jax.numpy as jnp

6+import numpy as np

7+

8+from cymsim.config import ForceExciter

9+from cymsim.simulate import Simulator

10+

11+DELTA = 1e-6

12+

13+

14+def _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

22+

23+

24+def _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

31+

32+

33+def 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)

53+

54+

55+def _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

62+

63+

64+def 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

89+

90+

91+def 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)