1"""The time loop: leapfrog integration of the shell with contact, damping, rigid-mode projection
2and radiation, block-compiled with JAX.
4One sub-step advances the state (q, v) of the shell, the state (zs, vs) of every exciter and the
5two retarded-time pressure rings:
7 f = f_int(q) - contact(q, zs) - 2 sigma1 rho h Lap(v_w)
8 v' = P_rigid [ (v + dt M^-1 f) / (1 + 2 sigma0 dt) ]
9 q' = q + dt v'
10 vs' = vs + dt F / m, zs' = zs + dt vs'
12where P_rigid removes the six rigid-body modes. The sub-step rate is an integer multiple m of the
13audio rate chosen from the linear stability bound. The acceleration that radiates, and that the
14optional record keeps, is the normal component of (v' - v) / dt, which for this scheme is the
15exact central second difference of q. It is decimated to the audio rate inside the loop by the
16polyphase accumulators of ``decimation.py``, so that no second difference is ever taken on
17stored displacements. At the optional probe points (``output.probes``) the same acceleration is
18interpolated, kept at the sub-step rate for one block, and decimated on the host with the same
19kernel as the listener pressure. See docs/model.md, section "Time integration", and docs/format.md.
20"""
21from __future__ import annotations
23import logging
24import math
25import time
26from dataclasses import dataclass, field
28import jax
29import jax.numpy as jnp
30import numpy as np
31from jax import lax
32from scipy.io import wavfile
34from .config import SimulationConfig
35from .decimation import DECIMATION_SPAN, StreamingDecimator, decimation_taps
36from .exciter import HitSchedule, build_patch
37from .geometry import CymbalProfile
38from .mesh import Mesh, disk_mesh
39from .radiation import RadiationModel
40from .shell import N_DOF, ShellModel
42log = logging.getLogger("cymsim")
46class HitSummary:
47 time: float
48 kind: str
49 peak_force: float # [N]
50 contact_time: float # total time with F > 0 [s]
51 rebound_velocity: float # exciter velocity at the end of the run [m/s], positive = upward
52 energy_in: float # kinetic energy at release [J]
53 energy_out: float # kinetic energy left in the exciter [J]
57class EnergyRecord:
58 time: float
59 kinetic: float
60 membrane: float
61 bending: float
62 exciters: float
63 max_w_over_h: float # nonlinearity meter: < 1 effectively linear, 3 to 5 clear cascade
66 def total(self) -> float:
67 return self.kinetic + self.membrane + self.bending + self.exciters
71class AccelerationRecord:
72 """Normal acceleration of every node at the audio rate, int16 against a scale per sample.
74 ``data[j, i]`` * ``scale[j]`` is the acceleration of node i along the rest-surface normal
75 (front side positive) at t = j / fs, low-pass filtered by the kernel ``taps`` (see
76 decimation.py). Sample j is final once the run has passed t = (j + W/2) / fs.
77 """
78 fs: int
79 data: np.ndarray # (T, N) int16
80 scale: np.ndarray # (T,) float32
81 taps: np.ndarray # (W, m) float64 polyphase kernel
84 def n_samples(self) -> int:
85 return len(self.data)
88 def time(self) -> np.ndarray:
89 return np.arange(self.n_samples) / self.fs
91 def as_float(self, start: int = 0, stop: int | None = None) -> np.ndarray:
92 """Decoded accelerations [m/s^2] for samples [start, stop), float32 (T', N)."""
93 return self.data[start:stop].astype(np.float32) * self.scale[start:stop, None]
96 def encode(a: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
97 """(int16 data, float32 scale) of a float block (T', N)."""
98 peak = np.abs(a).max(axis=1)
99 scale = np.where(peak > 0, peak / 32767.0, 1.0).astype(np.float32)
100 return np.round(a / scale[:, None]).astype(np.int16), np.where(peak > 0, scale, 0.0).astype(np.float32)
103 def empty(cls, fs: int, n_nodes: int, taps: np.ndarray) -> "AccelerationRecord":
104 return cls(fs, np.zeros((0, n_nodes), np.int16), np.zeros(0, np.float32), taps)
106 def extend(self, data: np.ndarray, scale: np.ndarray) -> None:
107 self.data = np.concatenate([self.data, data])
108 self.scale = np.concatenate([self.scale, scale])
112class DisplacementRecord:
113 """Transverse displacement w of every node [m] at a low frame rate, for display."""
114 rate: float
115 time: np.ndarray # (Td,)
116 w: np.ndarray # (Td, N) float32
119 def n_frames(self) -> int:
120 return len(self.time)
123 def empty(cls, rate: float, n_nodes: int) -> "DisplacementRecord":
124 return cls(rate, np.zeros(0), np.zeros((0, n_nodes), np.float32))
126 def extend(self, time: np.ndarray, w: np.ndarray) -> None:
127 self.time = np.concatenate([self.time, time])
128 self.w = np.concatenate([self.w, w])
132class SimulatorState:
133 """Everything the Python integrator needs to continue a run from step ``n_steps_done``.
135 Written to the ``state`` group of a .cymsim file. It is specific to this implementation
136 (the web companion keeps its own); readers of the format may ignore it.
137 """
138 n_steps_done: int
139 q: np.ndarray # (N, 5)
140 v: np.ndarray # (N, 5)
141 exciter_z: np.ndarray # (K,)
142 exciter_v: np.ndarray # (K,)
143 peak_force: np.ndarray # (K,)
144 contact_steps: np.ndarray # (K,) int32
145 ring_front: np.ndarray # (L,) pending retarded-time pressure, sub-step rate
146 ring_back: np.ndarray # (L,)
147 decimator_tail_front: np.ndarray # sub-steps still needed by non-final audio samples
148 decimator_tail_back: np.ndarray
149 acceleration_ring: np.ndarray | None = None # (W, N) partial audio samples of the record
153class SimulationResult:
154 config: SimulationConfig
155 fs: int
156 dt: float
157 substeps: int
158 mesh: Mesh
159 z0_node: np.ndarray # (N,) rest height of the nodes [m]
160 h_node: np.ndarray # (N,) thickness at the nodes [m]
161 normal_node: np.ndarray # (N, 3) unit normal of the rest surface, front side
162 r_front: np.ndarray # (E,) element centroid to listener [m]
163 r_back: np.ndarray # (E,) shortest path over the rim to the listener [m]
164 pressure_rayleigh: np.ndarray # baffled Rayleigh integral at the listener [Pa]
165 pressure_unbaffled: np.ndarray # front minus diffracted back term [Pa]
166 hits: list[HitSummary]
167 energy_log: list[EnergyRecord] = field(default_factory=list)
168 acceleration: AccelerationRecord | None = None
169 displacement: DisplacementRecord | None = None
170 complete: bool = True
171 n_steps_done: int = 0
172 state: SimulatorState | None = None
173 wall_time: float = 0.0
174 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]
179 def n_nodes(self) -> int:
180 return self.mesh.n_nodes
183 def n_elements(self) -> int:
184 return self.mesh.n_elements
187 def time(self) -> np.ndarray:
188 return np.arange(len(self.pressure_rayleigh)) / self.fs
191 def duration_done(self) -> float:
192 """Simulated time covered by the run so far [s]."""
193 return self.n_steps_done * self.dt
195 def write_wav(self, prefix: str, normalize: bool = True) -> list[str]:
196 """Write <prefix>_rayleigh.wav and <prefix>_unbaffled.wav as 16-bit PCM. With
197 ``normalize`` the peak is placed at 0.9 full scale; otherwise 1 Pa maps to 0.05 FS."""
198 paths = []
199 for name, p in (("rayleigh", self.pressure_rayleigh), ("unbaffled", self.pressure_unbaffled)):
200 peak = float(np.abs(p).max()) if len(p) else 0.0
201 scale = 0.9 / (peak or 1.0) if normalize else 0.05
202 path = f"{prefix}_{name}.wav"
203 wavfile.write(path, self.fs, (np.clip(p * scale, -1, 1) * 32767).astype(np.int16))
204 paths.append(path)
205 return paths
207 def write_npz(self, prefix: str) -> str:
208 path = f"{prefix}_result.npz"
209 np.savez(path, fs=self.fs, dt=self.dt, pressure_rayleigh=self.pressure_rayleigh,
210 pressure_unbaffled=self.pressure_unbaffled,
211 energy_time=[e.time 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)))
215 return path
217 def save(self, path: str) -> str:
218 """Write a .cymsim file (see docs/format.md)."""
219 from .cymsim_file import save_cymsim
220 return save_cymsim(self, path)
223class Simulator:
224 """Build every part of the model from a SimulationConfig; ``run()`` integrates it."""
226 def __init__(self, config: SimulationConfig):
227 config.validate()
228 self.config = config
229 num, mat, out = config.numerics, config.material, config.output
230 t0 = time.time()
232 self.profile = CymbalProfile(config.profile, config.imperfections)
233 self.mesh = disk_mesh(config.profile.radius, config.mesh.element_size, config.mesh.smoothing)
234 log.info("mesh: %d nodes, %d triangles, %d DOFs (nonlinear = %s)",
235 self.mesh.n_nodes, self.mesh.n_elements, N_DOF * self.mesh.n_nodes, num.nonlinear)
237 self.shell = ShellModel(self.mesh, self.profile, mat, num)
238 sh = self.shell
239 rim = sh.z0_node[self.mesh.boundary]
240 log.info("profile: height %.1f mm, rim warp p-p %.2f mm, %d dents, h in [%.2f, %.2f] mm, "
241 "max |grad z0| = %.3f, mass %.3f kg",
242 sh.z0_node.max() * 1e3, (rim.max() - rim.min()) * 1e3, config.imperfections.n_dents,
243 sh.h_e.min() * 1e3, sh.h_e.max() * 1e3, sh.max_slope, sh.mass)
244 if sh.max_slope > config.profile.max_slope:
245 log.warning("max slope %.2f exceeds the shallow-shell bound %.2f", sh.max_slope, config.profile.max_slope)
247 # ---- time step: integer number of sub-steps per audio sample, inside the stability bound.
248 # The run covers W/2 audio samples beyond the duration so that every sample of the
249 # duration is final after the decimation kernel's look-ahead.
250 self.fs = num.fs
251 self.substeps = math.ceil(1.0 / (self.fs * num.cfl_safety * sh.dt_stable))
252 self.dt = 1.0 / (self.fs * self.substeps)
253 self.n_out = int(round(self.fs * num.duration))
254 self.span = DECIMATION_SPAN
255 self.n_blocks = math.ceil((self.n_out + self.span // 2) * self.substeps / num.steps_per_block)
256 self.n_steps = self.n_blocks * num.steps_per_block
257 xc, yc = self.mesh.centroid[sh.limiting_element]
258 log.info("omega_max: element bound %.3e rad/s (limiting element at r = %.3f m, h = %.2f mm), "
259 "global %.3e rad/s (%d power iterations); %s bound used: dt = %.3e s, "
260 "%d sub-steps per sample, %d steps",
261 sh.omega_max, math.hypot(xc, yc), sh.h_e[sh.limiting_element] * 1e3, sh.omega_global,
262 sh.omega_global_iterations, "global" if sh.omega_step < sh.omega_max else "element",
263 self.dt, self.substeps, self.n_steps)
264 self.taps = decimation_taps(self.substeps)
266 # ---- hits
267 self.patches = [build_patch(h, self.mesh, self.profile) for h in config.hits]
268 self.schedule = HitSchedule(self.patches, self.dt, sh.dtype)
269 for k, p in enumerate(self.patches):
270 tau = p.free_contact_time()
271 if p.prescribed:
272 log.info("hit %d at t = %.3f s: prescribed raised-cosine force at r = %.3f m, peak %.1f N, "
273 "duration %.2f ms", k, p.time, p.hit.r, p.force_peak, p.force_duration * 1e3)
274 continue
275 log.info("hit %d at t = %.3f s: %s, %d-point patch at r = %.3f m, k = %.3g N/m^%.2f, "
276 "m = %.3f kg, v = %.1f m/s, energy %.2f J, free contact ~ %.2f ms "
277 "(direct excitation below ~%.0f Hz)",
278 k, p.time, p.hit.exciter.kind, p.n_points, p.hit.r, p.k_total, p.alpha,
279 p.mass, p.velocity, p.energy, tau * 1e3, 0.5 / tau)
280 if p.time >= num.duration:
281 log.warning("hit %d is scheduled after the end of the run", k)
283 # ---- radiation
284 self.radiation = RadiationModel(self.mesh, self.profile, sh.z0_node, config.listener,
285 self.dt, num.steps_per_block, sh.dtype)
287 # ---- probes: containing element and barycentric weights of each point
288 self.probe_xy = np.array(out.probes, dtype=np.float64).reshape(-1, 2)
289 loc = [self.mesh.locate(x, y) for x, y in self.probe_xy]
290 self.probe_nodes = jnp.asarray(np.array([l[0] for l in loc]).reshape(-1, 3), dtype=jnp.int32)
291 self.probe_w = jnp.asarray(np.array([l[1] for l in loc]).reshape(-1, 3), dtype=sh.dtype)
293 # ---- records
294 self.record = bool(out.record_acceleration)
295 self.disp_rate = float(out.displacement_rate)
296 if self.record:
297 nbytes = self.n_out * (2 * self.mesh.n_nodes + 4)
298 log.info("acceleration record: %d samples x %d nodes, int16, %.2f GB in memory",
299 self.n_out, self.mesh.n_nodes, nbytes / 1e9)
300 if nbytes > 2e9:
301 log.warning("the acceleration record needs %.1f GB; consider a shorter duration or "
302 "a coarser mesh", nbytes / 1e9)
303 S = self.fs * self.substeps # sub-steps per second
304 if self.disp_rate > 0:
305 n_frames = int(math.floor(self.n_out * self.substeps * self.disp_rate / S)) + 1
306 self.frame_steps = np.round(np.arange(n_frames) * S / self.disp_rate).astype(np.int64)
307 else:
308 self.frame_steps = np.zeros(0, np.int64)
309 self._build_step()
310 self.setup_time = time.time() - t0
311 log.info("setup done in %.1f s; backend %s", self.setup_time, jax.devices()[0].platform)
313 # ---- one sub-step and its block loop
314 def _build_step(self):
315 sh, sched, rad = self.shell, self.schedule, self.radiation
316 dt, sigma0 = self.dt, float(self.config.material.sigma0) # a numpy scalar would promote float32
317 m, W, N = self.substeps, self.span, self.mesh.n_nodes
318 spb = self.config.numerics.steps_per_block
319 tri, normal = sh.tri_j, sh.normal_j
320 record, rate = self.record, self.disp_rate
321 taps_j = jnp.asarray(self.taps, sh.dtype)
322 S = float(self.fs * m)
323 self.acc_rows = spb // m + 2 # audio samples emitted per block, at most
324 self.disp_rows = int(math.ceil(spb * rate / S)) + 2 if rate > 0 else 0
326 n_probes = len(self.probe_xy)
327 p_nodes, p_w = self.probe_nodes, self.probe_w
329 def step(n, k_block0, j_block0, st):
330 q, v, zs, vs, bf, bb, Fmax, ncont, ring, acc_out, disp_out, pbuf = st
331 if rate > 0:
332 # frame k sits at sub-step round(k S / rate); emit q at step n if n is one of them
333 k_cand = jnp.round(n * (rate / S))
334 is_frame = jnp.round(k_cand * (S / rate)) == n
335 row = jnp.clip(k_cand.astype(jnp.int32) - k_block0, 0, self.disp_rows - 1)
336 disp_out = disp_out.at[row].add(jnp.where(is_frame, q[:, 2], 0.0))
337 Fk, F, zs, vs = sched.contact(n, q, zs, vs)
338 fo = sh.internal_force(q)
339 fo = sched.apply_to_plate(fo, Fk)
340 fo = fo.at[:, 2].add(-sh.damping_force_w(v[:, 2]))
341 v1 = (v + dt * fo * sh.Minv_j) / (1 + 2 * sigma0 * dt) # sigma0 implicit
342 v1 = sh.project_rigid(v1)
343 q1 = q + dt * v1
344 zs1, vs1 = sched.advance(zs, vs, F, dt)
345 # normal acceleration at step n: central second difference of q, exact for leapfrog
346 a_node = (((v1 - v) / dt)[:, :3] * normal).sum(1)
347 acc_e = a_node[tri].mean(1)
348 bf, bb = rad.scatter(n, acc_e, bf, bb)
349 if n_probes:
350 pbuf = pbuf.at[n % spb].set((a_node[p_nodes] * p_w).sum(1)) # blocks start at multiples of spb
351 if record:
352 j0, kk = n // m, n % m
353 emit = kk == 0
354 j_e = j0 - W // 2 # the sample that just became final
355 slot_e = j_e % W
356 done = ring[slot_e]
357 row = jnp.clip(j_e - j_block0, 0, self.acc_rows - 1)
358 acc_out = acc_out.at[row].add(jnp.where(emit, done, 0.0))
359 ring = ring.at[slot_e].set(jnp.where(emit, 0.0, done))
360 # slot s holds sample j0 + 1 - W/2 + w with s = (j0 + 1 - W/2 + w) mod W
361 r = (j0 + 1 - W // 2) % W
362 tw = jnp.take(taps_j[:, kk], (jnp.arange(W) - r) % W)
363 ring = ring + tw[:, None] * a_node[None, :]
364 return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0), ring, acc_out, disp_out, pbuf)
367 def run_block(n0, k_block0, j_block0, st):
368 return lax.fori_loop(0, spb, lambda k, s: step((n0 + k).astype(jnp.int32), k_block0, j_block0, s), st)
371 def energies(st):
372 q, v, _, vs = st[0], st[1], st[2], st[3]
373 Em, Eb = sh.energy_parts_per_element(q)
374 return (sh.kinetic_energy(v), Em.sum(), Eb.sum(), sched.energy(vs),
375 jnp.abs(q[:, 2] / sh.h_node_j).max())
378 def clear(buf, idx):
379 return buf.at[idx].set(0.0)
381 self._step, self._run_block, self._energies, self._clear = step, run_block, energies, clear
383 def initial_state(self):
384 N, dtype = self.mesh.n_nodes, self.shell.dtype
385 zs, vs = self.schedule.initial_state()
386 bf, bb = self.radiation.initial_buffers()
387 K = self.schedule.n_hits
388 ring = jnp.zeros((self.span, N) if self.record else (0,), dtype)
389 return (jnp.zeros((N, N_DOF), dtype), jnp.zeros((N, N_DOF), dtype), zs, vs, bf, bb,
390 jnp.zeros(K, dtype), jnp.zeros(K, jnp.int32), ring, self._empty_outputs()[0],
391 self._empty_outputs()[1], jnp.zeros((self.config.numerics.steps_per_block, len(self.probe_xy)), dtype))
393 def _empty_outputs(self):
394 N, dtype = self.mesh.n_nodes, self.shell.dtype
395 return (jnp.zeros((self.acc_rows, N) if self.record else (0,), dtype),
396 jnp.zeros((self.disp_rows, N) if self.disp_rate > 0 else (0,), dtype))
398 def _hit_summaries(self, Fmax, ncont, vs) -> list[HitSummary]:
399 hits = []
400 for k, p in enumerate(self.patches):
401 hits.append(HitSummary(p.time, p.hit.exciter.kind, float(Fmax[k]), float(ncont[k]) * self.dt,
402 float(vs[k]), p.energy, 0.5 * p.mass * float(vs[k]) ** 2))
403 return hits
405 def _result(self, state, energy_log, p_front, p_back, acc, disp, n_done, wall, dec_f, dec_b,
406 probe=None) -> SimulationResult:
407 q, v, zs, vs, bf, bb, Fmax, ncont, ring = jax.device_get(state[:9])
408 sim_state = SimulatorState(
409 n_done, np.asarray(q), np.asarray(v), np.asarray(zs), np.asarray(vs), np.asarray(Fmax),
410 np.asarray(ncont), np.asarray(bf), np.asarray(bb), dec_f.tail.copy(), dec_b.tail.copy(),
411 np.asarray(ring) if self.record else None)
412 n = min(len(p_front), len(p_back))
413 p_front, p_back = np.asarray(p_front[:n]), np.asarray(p_back[:n])
414 from . import __version__
415 return SimulationResult(
416 self.config, self.fs, self.dt, self.substeps, self.mesh, self.shell.z0_node, self.shell.h_node,
417 self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front, p_front - p_back,
418 self._hit_summaries(Fmax, ncont, vs), list(energy_log), acc, disp, n_done >= self.n_steps, n_done,
419 sim_state, wall, f"cymsim {__version__}",
420 self.probe_xy if len(self.probe_xy) else None, probe)
422 def run(self, save_path: str | None = None, save_every: float | None = None) -> SimulationResult:
423 """Integrate the run. With ``save_path`` a .cymsim file is written at the end and, if
424 ``save_every`` (seconds of wall time) is given, also periodically during the run, each
425 time as a consistent partial file that can be loaded (see docs/format.md)."""
426 num = self.config.numerics
427 m, spb, N = self.substeps, num.steps_per_block, self.mesh.n_nodes
428 state = self.initial_state()
429 energy_log: list[EnergyRecord] = []
430 dec_f, dec_b = StreamingDecimator(m), StreamingDecimator(m)
431 audio_f, audio_b = [], []
432 n_audio = 0
433 acc = AccelerationRecord.empty(self.fs, N, self.taps) if self.record else None
434 disp = DisplacementRecord.empty(self.disp_rate, N) if self.disp_rate > 0 else None
435 n_probes = len(self.probe_xy)
436 dec_p = [StreamingDecimator(m) for _ in range(n_probes)]
437 audio_p = []
438 probe_now = lambda: np.concatenate(audio_p) if n_probes else None
439 writer = None
440 if save_path is not None:
441 from .cymsim_file import CymsimWriter
442 writer = CymsimWriter(save_path, self._result(state, [], np.zeros(0), np.zeros(0), acc, disp, 0, 0.0,
443 dec_f, dec_b))
444 t_run = time.time()
445 t_saved = t_run
446 report_every = max(1, self.n_blocks // max(1, num.n_progress))
447 log.info("running %d blocks of %d steps (nonlinearity meter: max|w|/h < 1 effectively "
448 "linear; 3 to 5 clear cascade)", self.n_blocks, spb)
449 for blk in range(self.n_blocks):
450 n0, n1 = blk * spb, (blk + 1) * spb
451 k0 = int(np.searchsorted(self.frame_steps, n0))
452 k1 = int(np.searchsorted(self.frame_steps, n1))
453 j_block0 = -(-n0 // m) - self.span // 2
454 acc_out, disp_out = self._empty_outputs()
455 state = state[:9] + (acc_out, disp_out) + state[11:]
456 state = self._run_block(jnp.int32(n0), jnp.int32(k0), jnp.int32(j_block0), state)
458 # ---- pressure: the ring slots of this block's sub-steps are final; read, clear, decimate
459 slots = jnp.asarray(self.radiation.slots(n0, n1))
460 pf, pb = np.asarray(state[4][slots], np.float64), np.asarray(state[5][slots], np.float64)
461 state = state[:4] + (self._clear(state[4], slots), self._clear(state[5], slots)) + state[6:]
462 new_f, new_b = dec_f.push(pf), dec_b.push(pb)
463 keep = max(0, min(len(new_f), self.n_out - n_audio))
464 audio_f.append(new_f[:keep])
465 audio_b.append(new_b[:keep])
466 if n_probes:
467 rows = np.asarray(state[11], np.float64)
468 new_p = np.stack([d.push(rows[:, i]) for i, d in enumerate(dec_p)], 1)
469 audio_p.append(new_p[:keep] if len(new_p) else np.zeros((0, n_probes)))
470 n_audio += keep
472 # ---- records
473 if self.record:
474 n_emit = -(-n1 // m) - -(-n0 // m)
475 rows = np.asarray(state[9][:n_emit], np.float32)
476 j = j_block0 + np.arange(n_emit)
477 sel = (j >= 0) & (j < self.n_out)
478 if sel.any():
479 acc.extend(*AccelerationRecord.encode(rows[sel]))
480 if self.disp_rate > 0 and k1 > k0:
481 frames = np.asarray(state[10][:k1 - k0], np.float32)
482 t_frames = self.frame_steps[k0:k1] * self.dt
483 sel = self.frame_steps[k0:k1] <= self.n_out * m
484 disp.extend(t_frames[sel], frames[sel])
486 if blk % report_every == 0 or blk == self.n_blocks - 1:
487 ke, em, eb, ks, woh = (float(x) for x in self._energies(state))
488 t = n1 * self.dt
489 energy_log.append(EnergyRecord(t, ke, em, eb, ks, woh))
490 log.info(" t = %6.3f s KE %.3e membrane %.3e bending %.3e exciters %.3e J "
491 "max|w|/h = %5.2f [%.0f s wall]", t, ke, em, eb, ks, woh, time.time() - t_run)
492 if not math.isfinite(ke + em + eb):
493 raise FloatingPointError("energy is not finite: the run has blown up "
494 "(lower numerics.cfl_safety)")
495 if writer is not None and save_every is not None and time.time() - t_saved >= save_every \
496 and blk < self.n_blocks - 1:
497 res = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),
498 acc, disp, n1, time.time() - t_run, dec_f, dec_b, probe_now())
499 writer.update(res)
500 t_saved = time.time()
501 log.info(" saved partial run (%.3f s) to %s", res.duration_done, save_path)
503 result = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),
504 acc, disp, self.n_steps, time.time() - t_run, dec_f, dec_b, probe_now())
505 for h in result.hits:
506 log.info("hit (%s): peak force %.1f N, contact %.2f ms, rebound velocity %+.2f m/s "
507 "(energy left in exciter %.2f of %.2f J)", h.kind, h.peak_force,
508 h.contact_time * 1e3, h.rebound_velocity, h.energy_out, h.energy_in)
509 for name, p in (("rayleigh", result.pressure_rayleigh), ("unbaffled", result.pressure_unbaffled)):
510 peak = float(np.abs(p).max())
511 log.info("%-10s: peak %.3f Pa (%.1f dB SPL peak)", name, peak, 20 * np.log10(max(peak, 1e-12) / 2e-5))
512 if writer is not None:
513 writer.update(result)
514 writer.close()
515 log.info("wrote %s", save_path)
516 return result
519def simulate(config: SimulationConfig) -> SimulationResult:
520 return Simulator(config).run()