Add a prescribed-force exciter and point probes of the normal acceleration; accept Python 3.9 in the YAML loader - proof-of-concept/cymsim

proof-of-concept / cymsim

Add a prescribed-force exciter and point probes of the normal acceleration; accept Python 3.9 in the YAML loader

Jeremy Magland <jmagland@flatironinstitute.org> committed 6 days ago commit e19a54051204 parent e30c940 Browse files

6 changed files+134−25

README.mdmodified+1−1View file

@@ -45,7 +45,7 @@ from cymsim import load_cymsim

4545 back = load_cymsim("crash.cymsim") # the same SimulationResult, from the file

4646 ```

4747

48-A run is described by one `SimulationConfig` made of `Material`, `Profile`, `Imperfections`, `MeshConfig`, a list of `Hit`s (each with a `StickExciter` or `MalletExciter`), `Listener`, `Numerics` and `Output`. The YAML file mirrors the dataclass fields; omitted fields take the defaults, which reproduce an 18 inch medium-thin B20 crash. Units are SI throughout. See [examples/](examples/) for a mallet stroke, a stick crash, a train of four hits and a quick CPU run.

48+A run is described by one `SimulationConfig` made of `Material`, `Profile`, `Imperfections`, `MeshConfig`, a list of `Hit`s (each with a `StickExciter`, a `MalletExciter`, or a `ForceExciter` for a prescribed force pulse), `Listener`, `Numerics` and `Output`. The YAML file mirrors the dataclass fields; omitted fields take the defaults, which reproduce an 18 inch medium-thin B20 crash. Units are SI throughout. See [examples/](examples/) for a mallet stroke, a stick crash, a train of four hits and a quick CPU run.

4949

5050 ## Layout

5151

docs/model.mdmodified+13−1View file

@@ -162,7 +162,7 @@ Set-up (element stiffness, eigenvalue bound, rigid modes) is done in float64; th

162162

163163 ## 6. Excitation

164164

165-Implemented in [exciter.py](../src/cymsim/exciter.py); parameters in `Hit`, `StickExciter`, `MalletExciter`.

165+Implemented in [exciter.py](../src/cymsim/exciter.py); parameters in `Hit`, `StickExciter`, `MalletExciter`, `ForceExciter`.

166166

167167 ### 6.1 Contact law

168168

@@ -188,6 +188,16 @@ which for $\alpha = 1.5$ is Hertz's $2.94\, \delta_{max} / v$, together with the

188188

189189 A run has a list of hits, each with its own time, position, velocity and exciter. Every exciter exists for the whole run as a row of fixed-shape arrays (hits $\times$ largest patch, padded with zero weights and a very negative offset so that padded points never touch), which lets one compiled step serve any number of hits. Before its release step $n_k = \mathrm{round}(t_k / \Delta t)$ an exciter is inactive: it produces no force and is held far above the shell. At step $n_k$ it is placed at $z_s = \max_p (w_p + o_p) + \text{gap}$, just above the *current* surface, with velocity $-v_k$. This matters for a train of hits, where the surface under a later hit is already moving by several thickness units. After the rebound the exciter simply moves away at constant speed; there is no gravity and no hand. Two consequences follow: an exciter never strikes twice on its own, and if the shell were to overtake a slowly rebounding exciter a second contact would occur (a real player pulls the stick back). The summary of each hit reports the peak force, the total time with $F > 0$ and the rebound velocity.

190190

191+### 6.4 Prescribed force

192+

193+A `ForceExciter` replaces the contact by a force given in advance, a raised cosine of peak $F_0$ and length $T$ applied at one point from the hit time $t_k$:

194+

195+$$

196+F(t) = \tfrac12 F_0 \left(1 - \cos \frac{2\pi (t - t_k)}{T}\right), \qquad t_k \le t \le t_k + T,

197+$$

198+

199+distributed to the containing element's $w$ DOFs with the barycentric weights and pushing the surface down. The force does not depend on the motion, so there is no rebound, and the hit's velocity is ignored. This is not a model of a stick or a mallet. It is the excitation used by modal plate models (Ducceschi and Touzé 2015, among others), and it exists so that cymsim can be run on exactly their problems.

200+

191201 ## 7. Radiation

192202

193203 Implemented in [radiation.py](../src/cymsim/radiation.py) and [decimation.py](../src/cymsim/decimation.py); parameters in `Listener`.

@@ -218,6 +228,8 @@ That reduction has to be a proper low-pass filter, because the cascade puts real

218228

219229 The record is kept only on request (`output.record_acceleration`, default off), because it is large: one value per node per audio sample, stored as `int16` against one `float32` scale per sample, is $2 N f_s$ bytes per simulated second, about 0.4 GB per second at the default mesh. The per-sample scale keeps the quantisation noise at $1.5 \times 10^{-5}$ of each sample's own peak through the 90 dB or more of a stroke's decay; the noise is incoherent across nodes and changes the built-in Rayleigh sum by less than $10^{-4}$ relative rms. Transverse displacements, in contrast, are kept only at a low frame rate (`output.displacement_rate`, default 60 Hz), in float32 and exactly at the frame times. At 60 Hz they are a strobe of the motion, useful for display and for the slow quasi-static deflection, but not a signal from which anything should be differentiated.

220230

231+For comparisons at a single point, `output.probes` lists points $(x, y)$ on the projected domain where the normal acceleration is kept at the audio rate in float64: the nodal values are interpolated with the barycentric weights of the containing element, held for one block at the sub-step rate, and decimated on the host with the same kernel. A probe placed on a node reproduces that node's column of the acceleration record. Probe signals are returned in the result and its `.npz` file but are not yet written to the `.cymsim` file.

232+

221233 Both records, the two listener signals, the configuration, the mesh with its nodal normals, the energy log, the hit summaries and the integrator state are written to an HDF5 file with the extension `.cymsim`, described in [format.md](format.md). The time series are resizable, so `cymsim run -o run.cymsim --save-every 60` writes a consistent partial file every minute and extends it, and a run that is interrupted leaves a loadable file. The web companion writes and reads the same format.

222234

223235 ## 10. Limitations and future work

src/cymsim/config.pymodified+25−6View file

@@ -5,7 +5,7 @@ by what it describes: the material, the cymbal profile, the seeded imperfections

55 list of hits (each with its own exciter), the listener, the numerics and the output. Defaults

66 reproduce the 18-inch medium-thin crash of the original prototype. A SimulationConfig can be

77 built directly in Python or loaded from a YAML file whose keys mirror the field names; a hit's

8-exciter is chosen by its ``kind`` key ("stick" or "mallet").

8+exciter is chosen by its ``kind`` key ("stick", "mallet" or "force").

99 """

1010 from __future__ import annotations

1111

@@ -91,8 +91,20 @@ class MalletExciter:

9191 gap: float = 1e-4

9292

9393

94-Exciter = typing.Union[StickExciter, MalletExciter]

95-EXCITER_KINDS: dict[str, type] = {"stick": StickExciter, "mallet": MalletExciter}

94+@dataclass

95+class ForceExciter:

96+ """Prescribed force at one point, independent of the motion: a raised cosine

97+ F(t) = peak (1 - cos(2 pi t / duration)) / 2 for 0 <= t <= duration after the hit time,

98+ pushing the surface down. There is no contact and no rebound. This is the excitation of

99+ modal plate models such as Ducceschi and Touze (2015), and is meant for comparison with them;

100+ ``Hit.velocity`` is ignored."""

101+ kind: str = "force"

102+ peak: float = 100.0 # peak force [N]

103+ duration: float = 2e-3 # total length of the pulse [s]

104+

105+

106+Exciter = typing.Union[StickExciter, MalletExciter, ForceExciter]

107+EXCITER_KINDS: dict[str, type] = {"stick": StickExciter, "mallet": MalletExciter, "force": ForceExciter}

96108

97109

98110 @dataclass

@@ -137,6 +149,7 @@ class Output:

137149 cascade_plot: bool = True # needs matplotlib

138150 record_acceleration: bool = False # keep the normal acceleration of every node at the audio rate

139151 displacement_rate: float = 60.0 # frames per second of transverse displacement kept for display; 0 = none

152+ probes: list[tuple[float, float]] = field(default_factory=list) # points (x, y) [m] where the normal acceleration is kept at the audio rate

140153

141154

142155 # ----------------------------------------------------------------------------- simulation

@@ -170,11 +183,17 @@ class SimulationConfig:

170183 if h.velocity <= 0:

171184 raise ValueError(f"hit {i}: velocity must be positive (the exciter approaches from above)")

172185 ex = h.exciter

173- extent = ex.patch_length / 2 if isinstance(ex, StickExciter) else ex.radius

186+ extent = (ex.patch_length / 2 if isinstance(ex, StickExciter)

187+ else ex.radius if isinstance(ex, MalletExciter) else 0.0)

174188 if h.r + extent >= p.radius - self.mesh.element_size:

175189 raise ValueError(f"hit {i}: contact patch reaches the rim (r = {h.r}, extent = {extent})")

176190 if ex.kind not in EXCITER_KINDS:

177191 raise ValueError(f"hit {i}: unknown exciter kind {ex.kind!r}")

192+ if isinstance(ex, ForceExciter) and (ex.peak < 0 or ex.duration <= 0):

193+ raise ValueError(f"hit {i}: a force pulse needs peak >= 0 and duration > 0")

194+ for i, (x, y) in enumerate(self.output.probes):

195+ if math.hypot(x, y) >= p.radius:

196+ raise ValueError(f"output.probes[{i}] = ({x}, {y}) lies outside the cymbal")

178197

179198 # ---- dict / YAML

180199 def to_dict(self) -> dict:

@@ -214,7 +233,7 @@ def _build(tp, value):

214233 return [_build(item, v) for v in value]

215234 if origin is tuple:

216235 return tuple(value)

217- if origin in (typing.Union, types.UnionType):

236+ if origin in (typing.Union, getattr(types, "UnionType", typing.Union)): # UnionType: Python >= 3.10

218237 # The only union in the schema is the exciter, discriminated by its "kind" key.

219238 if isinstance(value, dict) and "kind" in value:

220239 try:

@@ -224,7 +243,7 @@ def _build(tp, value):

224243 f"expected one of {sorted(EXCITER_KINDS)}") from None

225244 if any(isinstance(value, t) for t in typing.get_args(tp)):

226245 return value

227- raise ValueError("an exciter mapping needs a 'kind' key ('stick' or 'mallet')")

246+ raise ValueError("an exciter mapping needs a 'kind' key ('stick', 'mallet' or 'force')")

228247 if tp is float and isinstance(value, (int, float)):

229248 return float(value)

230249 return value

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

@@ -4,7 +4,9 @@ Each hit is a rigid body of mass m that moves only along z, released at the hit

44 clearance above the surface with downward speed v, and coupled to the shell through P springs

55 with the power law F = k delta^alpha, delta being the local penetration. A stick contacts along

66 a short line, a mallet on a small disk. Because the surface is curved, the rigid face is the

7-best-fit plane through the patch points and the residual heights give a staggered contact. All

7+best-fit plane through the patch points and the residual heights give a staggered contact. A

8+"force" exciter instead applies a prescribed raised-cosine force at one point, for comparison

9+with models driven by a given force; it has no contact and no motion of its own. All

810 hits of a run are stacked into fixed-shape arrays (padded to the largest patch) so that the time

911 loop handles any number of hits with the same compiled code. See docs/model.md, section

1012 "Excitation".

@@ -18,7 +20,7 @@ import jax.numpy as jnp

1820 import numpy as np

1921 from scipy.integrate import quad

2022

21-from .config import Hit, MalletExciter, StickExciter

23+from .config import ForceExciter, Hit, MalletExciter, StickExciter

2224 from .geometry import CymbalProfile

2325 from .mesh import Mesh

2426

@@ -40,6 +42,12 @@ class ContactPatch:

4042 gap: float

4143 velocity: float

4244 time: float

45+ force_peak: float = 0.0 # prescribed pulse (force exciter only) [N]

46+ force_duration: float = 0.0

47+

48+ @property

49+ def prescribed(self) -> bool:

50+ return self.force_duration > 0

4351

4452 @property

4553 def n_points(self) -> int:

@@ -56,7 +64,9 @@ class ContactPatch:

5664

5765 def free_contact_time(self) -> float:

5866 """Contact duration against a rigid wall, 2 (d_max / v) int_0^1 dx / sqrt(1 - x^(a+1)).

59- For alpha = 1.5 this is the Hertz value 2.94 d_max / v."""

67+ For alpha = 1.5 this is the Hertz value 2.94 d_max / v. For a prescribed force, the pulse length."""

68+ if self.prescribed:

69+ return self.force_duration

6070 a = self.alpha

6171 c, _ = quad(lambda s: 1.0 / math.sqrt(1.0 - s ** (a + 1)), 0.0, 1.0)

6272 return 2.0 * c * self.max_compression() / self.velocity

@@ -77,6 +87,10 @@ def build_patch(hit: Hit, mesh: Mesh, profile: CymbalProfile) -> ContactPatch:

7787 th = np.arange(6) * np.pi / 3

7888 px = np.concatenate([[xs], xs + ex.radius * np.cos(th)])

7989 py = np.concatenate([[ys], ys + ex.radius * np.sin(th)])

90+ elif isinstance(ex, ForceExciter):

91+ nodes, w = mesh.locate(xs, ys)

92+ return ContactPatch(hit, np.array([xs]), np.array([ys]), nodes[None], w[None], np.zeros(1), 0.0, 1.0,

93+ 1.0, 0.0, 0.0, float(hit.time), float(ex.peak), float(ex.duration))

8094 else:

8195 raise TypeError(f"unknown exciter type {type(ex).__name__}")

8296 loc = [mesh.locate(px[k], py[k]) for k in range(len(px))]

@@ -118,6 +132,11 @@ class HitSchedule:

118132 self.gap = f([p.gap for p in patches])

119133 self.velocity = f([p.velocity for p in patches])

120134 self.step = jnp.asarray([int(round(p.time / dt)) for p in patches], dtype=jnp.int32)

135+ self.is_force = jnp.asarray([p.prescribed for p in patches])

136+ self.f_peak = f([p.force_peak for p in patches])

137+ self.f_dur = f([p.force_duration if p.prescribed else 1.0 for p in patches])

138+ self.first_point = f(np.arange(P)[None, :] == 0)

139+ self.dt = dt

121140 self.dtype = dtype

122141

123142 def initial_state(self):

@@ -137,6 +156,11 @@ class HitSchedule:

137156 vs = jnp.where(release, -self.velocity, vs)

138157 delta = jnp.maximum(wp - zs[:, None], 0.0)

139158 Fk = jnp.where(active[:, None], self.k_point[:, None] * jnp.power(delta, self.alpha[:, None]), 0.0)

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

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

161+ 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)

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

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

141165

142166 def apply_to_plate(self, fo, Fk):

@@ -144,8 +168,9 @@ class HitSchedule:

144168 return fo.at[self.nodes, 2].add(-Fk[:, :, None] * self.weights)

145169

146170 def advance(self, zs, vs, F, dt):

147- """Symplectic Euler for the exciters (they are integrated with the plate's velocity update)."""

148- vs1 = vs + dt * F / self.mass

171+ """Symplectic Euler for the exciters (they are integrated with the plate's velocity update).

172+ A prescribed force does not move its exciter."""

173+ vs1 = vs + dt * jnp.where(self.is_force, 0.0, F) / self.mass

149174 return zs + dt * vs1, vs1

150175

151176 def energy(self, vs):

src/cymsim/simulate.pymodified+42−11View file

@@ -14,7 +14,9 @@ audio rate chosen from the linear stability bound. The acceleration that radiate

1414 optional record keeps, is the normal component of (v' - v) / dt, which for this scheme is the

1515 exact central second difference of q. It is decimated to the audio rate inside the loop by the

1616 polyphase accumulators of ``decimation.py``, so that no second difference is ever taken on

17-stored displacements. See docs/model.md, section "Time integration", and docs/format.md.

17+stored displacements. At the optional probe points (``output.probes``) the same acceleration is

18+interpolated, kept at the sub-step rate for one block, and decimated on the host with the same

19+kernel as the listener pressure. See docs/model.md, section "Time integration", and docs/format.md.

1820 """

1921 from __future__ import annotations

2022

@@ -170,6 +172,8 @@ class SimulationResult:

170172 state: SimulatorState | None = None

171173 wall_time: float = 0.0

172174 creator: str = ""

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

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

173177

174178 @property

175179 def n_nodes(self) -> int:

@@ -205,7 +209,9 @@ class SimulationResult:

205209 np.savez(path, fs=self.fs, dt=self.dt, pressure_rayleigh=self.pressure_rayleigh,

206210 pressure_unbaffled=self.pressure_unbaffled,

207211 energy_time=[e.time for e in self.energy_log],

208- energy=[[e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h] for e in self.energy_log])

212+ energy=[[e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h] for e in self.energy_log],

213+ **({} if self.probe_acceleration is None else

214+ dict(probes=self.probes, probe_acceleration=self.probe_acceleration)))

209215 return path

210216

211217 def save(self, path: str) -> str:

@@ -260,6 +266,10 @@ class Simulator:

260266 self.schedule = HitSchedule(self.patches, self.dt, sh.dtype)

261267 for k, p in enumerate(self.patches):

262268 tau = p.free_contact_time()

269+ if p.prescribed:

270+ log.info("hit %d at t = %.3f s: prescribed raised-cosine force at r = %.3f m, peak %.1f N, "

271+ "duration %.2f ms", k, p.time, p.hit.r, p.force_peak, p.force_duration * 1e3)

272+ continue

263273 log.info("hit %d at t = %.3f s: %s, %d-point patch at r = %.3f m, k = %.3g N/m^%.2f, "

264274 "m = %.3f kg, v = %.1f m/s, energy %.2f J, free contact ~ %.2f ms "

265275 "(direct excitation below ~%.0f Hz)",

@@ -272,6 +282,12 @@ class Simulator:

272282 self.radiation = RadiationModel(self.mesh, self.profile, sh.z0_node, config.listener,

273283 self.dt, num.steps_per_block, sh.dtype)

274284

285+ # ---- probes: containing element and barycentric weights of each point

286+ self.probe_xy = np.array(out.probes, dtype=np.float64).reshape(-1, 2)

287+ loc = [self.mesh.locate(x, y) for x, y in self.probe_xy]

288+ self.probe_nodes = jnp.asarray(np.array([l[0] for l in loc]).reshape(-1, 3), dtype=jnp.int32)

289+ self.probe_w = jnp.asarray(np.array([l[1] for l in loc]).reshape(-1, 3), dtype=sh.dtype)

290+

275291 # ---- records

276292 self.record = bool(out.record_acceleration)

277293 self.disp_rate = float(out.displacement_rate)

@@ -305,8 +321,11 @@ class Simulator:

305321 self.acc_rows = spb // m + 2 # audio samples emitted per block, at most

306322 self.disp_rows = int(math.ceil(spb * rate / S)) + 2 if rate > 0 else 0

307323

324+ n_probes = len(self.probe_xy)

325+ p_nodes, p_w = self.probe_nodes, self.probe_w

326+

308327 def step(n, k_block0, j_block0, st):

309- q, v, zs, vs, bf, bb, Fmax, ncont, ring, acc_out, disp_out = st

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

310329 if rate > 0:

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

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

@@ -325,6 +344,8 @@ class Simulator:

325344 a_node = (((v1 - v) / dt)[:, :3] * normal).sum(1)

326345 acc_e = a_node[tri].mean(1)

327346 bf, bb = rad.scatter(n, acc_e, bf, bb)

347+ if n_probes:

348+ pbuf = pbuf.at[n % spb].set((a_node[p_nodes] * p_w).sum(1)) # blocks start at multiples of spb

328349 if record:

329350 j0, kk = n // m, n % m

330351 emit = kk == 0

@@ -338,7 +359,7 @@ class Simulator:

338359 r = (j0 + 1 - W // 2) % W

339360 tw = jnp.take(taps_j[:, kk], (jnp.arange(W) - r) % W)

340361 ring = ring + tw[:, None] * a_node[None, :]

341- return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0), ring, acc_out, disp_out)

362+ return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0), ring, acc_out, disp_out, pbuf)

342363

343364 @jax.jit

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

@@ -365,7 +386,7 @@ class Simulator:

365386 ring = jnp.zeros((self.span, N) if self.record else (0,), dtype)

366387 return (jnp.zeros((N, N_DOF), dtype), jnp.zeros((N, N_DOF), dtype), zs, vs, bf, bb,

367388 jnp.zeros(K, dtype), jnp.zeros(K, jnp.int32), ring, self._empty_outputs()[0],

368- self._empty_outputs()[1])

389+ self._empty_outputs()[1], jnp.zeros((self.config.numerics.steps_per_block, len(self.probe_xy)), dtype))

369390

370391 def _empty_outputs(self):

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

@@ -379,8 +400,9 @@ class Simulator:

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

380401 return hits

381402

382- def _result(self, state, energy_log, p_front, p_back, acc, disp, n_done, wall, dec_f, dec_b) -> SimulationResult:

383- q, v, zs, vs, bf, bb, Fmax, ncont, ring, _, _ = jax.device_get(state)

403+ def _result(self, state, energy_log, p_front, p_back, acc, disp, n_done, wall, dec_f, dec_b,

404+ probe=None) -> SimulationResult:

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

384406 sim_state = SimulatorState(

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

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

@@ -392,7 +414,8 @@ class Simulator:

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

393415 self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front, p_front - p_back,

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

395- sim_state, wall, f"cymsim {__version__}")

417+ sim_state, wall, f"cymsim {__version__}",

418+ self.probe_xy if len(self.probe_xy) else None, probe)

396419

397420 def run(self, save_path: str | None = None, save_every: float | None = None) -> SimulationResult:

398421 """Integrate the run. With ``save_path`` a .cymsim file is written at the end and, if

@@ -407,6 +430,10 @@ class Simulator:

407430 n_audio = 0

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

409432 disp = DisplacementRecord.empty(self.disp_rate, N) if self.disp_rate > 0 else None

433+ n_probes = len(self.probe_xy)

434+ dec_p = [StreamingDecimator(m) for _ in range(n_probes)]

435+ audio_p = []

436+ probe_now = lambda: np.concatenate(audio_p) if n_probes else None

410437 writer = None

411438 if save_path is not None:

412439 from .cymsim_file import CymsimWriter

@@ -423,7 +450,7 @@ class Simulator:

423450 k1 = int(np.searchsorted(self.frame_steps, n1))

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

425452 acc_out, disp_out = self._empty_outputs()

426- state = state[:9] + (acc_out, disp_out)

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

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

428455

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

@@ -434,6 +461,10 @@ class Simulator:

434461 keep = max(0, min(len(new_f), self.n_out - n_audio))

435462 audio_f.append(new_f[:keep])

436463 audio_b.append(new_b[:keep])

464+ if n_probes:

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

466+ new_p = np.stack([d.push(rows[:, i]) for i, d in enumerate(dec_p)], 1)

467+ audio_p.append(new_p[:keep] if len(new_p) else np.zeros((0, n_probes)))

437468 n_audio += keep

438469

439470 # ---- records

@@ -462,13 +493,13 @@ class Simulator:

462493 if writer is not None and save_every is not None and time.time() - t_saved >= save_every \

463494 and blk < self.n_blocks - 1:

464495 res = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),

465- acc, disp, n1, time.time() - t_run, dec_f, dec_b)

496+ acc, disp, n1, time.time() - t_run, dec_f, dec_b, probe_now())

466497 writer.update(res)

467498 t_saved = time.time()

468499 log.info(" saved partial run (%.3f s) to %s", res.duration_done, save_path)

469500

470501 result = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),

471- acc, disp, self.n_steps, time.time() - t_run, dec_f, dec_b)

502+ acc, disp, self.n_steps, time.time() - t_run, dec_f, dec_b, probe_now())

472503 for h in result.hits:

473504 log.info("hit (%s): peak force %.1f N, contact %.2f ms, rebound velocity %+.2f m/s "

474505 "(energy left in exciter %.2f of %.2f J)", h.kind, h.peak_force,

tests/test_simulate.pymodified+23−1View file

@@ -3,7 +3,7 @@ import copy

33 import numpy as np

44 import pytest

55

6-from cymsim.config import Hit, MalletExciter, StickExciter

6+from cymsim.config import ForceExciter, Hit, MalletExciter, StickExciter

77 from cymsim.exciter import build_patch

88 from cymsim.simulate import Simulator

99

@@ -57,3 +57,25 @@ def test_contact_patch_geometry(coarse_config):

5757 mallet = build_patch(Hit(r=0.137, phi=0.3, velocity=8.0, exciter=MalletExciter()), sim.mesh, sim.profile)

5858 assert np.isclose(mallet.max_compression(), MalletExciter.compression)

5959 assert abs(mallet.offsets.max()) < 5e-3

60+

61+

62+def test_prescribed_force_and_probes(coarse_config):

63+ cfg = copy.deepcopy(coarse_config)

64+ cfg.hits = [Hit(time=0.002, r=0.15, phi=0.7, exciter=ForceExciter(peak=50.0, duration=2e-3))]

65+ cfg.output.record_acceleration = True

66+ sim = Simulator(cfg)

67+ node = int(np.argmin(np.hypot(sim.mesh.xy[:, 0] - 0.1, sim.mesh.xy[:, 1] - 0.05)))

68+ cfg.output.probes = [tuple(sim.mesh.xy[node]), (0.05, -0.12)]

69+ res = Simulator(cfg).run()

70+ (h,) = res.hits

71+ assert h.kind == "force"

72+ assert np.isclose(h.peak_force, 50.0, rtol=1e-3)

73+ assert np.isclose(h.contact_time, 2e-3, rtol=1e-2)

74+ assert h.energy_in == 0 and h.rebound_velocity == 0

75+ a = res.probe_acceleration

76+ n_out = int(cfg.numerics.fs * cfg.numerics.duration)

77+ assert a.shape == (n_out, 2) and np.isfinite(a).all()

78+ assert np.abs(a[: int(0.002 * cfg.numerics.fs) - 16]).max() < 1e-6 * np.abs(a).max() # nothing before the hit

79+ # a probe on a node sees that node's acceleration (same decimation kernel; int16 record)

80+ rec = res.acceleration.as_float()[:, node]

81+ assert np.abs(a[:, 0] - rec).max() < 1e-3 * np.abs(rec).max()