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 bound %.3e rad/s (limiting element at r = %.3f m, h = %.2f mm); "
259 "dt = %.3e s, %d sub-steps per sample, %d steps",
260 sh.omega_max, math.hypot(xc, yc), sh.h_e[sh.limiting_element] * 1e3,
261 self.dt, self.substeps, self.n_steps)
262 self.taps = decimation_taps(self.substeps)
264 # ---- hits
265 self.patches = [build_patch(h, self.mesh, self.profile) for h in config.hits]
266 self.schedule = HitSchedule(self.patches, self.dt, sh.dtype)
267 for k, p in enumerate(self.patches):
268 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
273 log.info("hit %d at t = %.3f s: %s, %d-point patch at r = %.3f m, k = %.3g N/m^%.2f, "
274 "m = %.3f kg, v = %.1f m/s, energy %.2f J, free contact ~ %.2f ms "
275 "(direct excitation below ~%.0f Hz)",
276 k, p.time, p.hit.exciter.kind, p.n_points, p.hit.r, p.k_total, p.alpha,
277 p.mass, p.velocity, p.energy, tau * 1e3, 0.5 / tau)
278 if p.time >= num.duration:
279 log.warning("hit %d is scheduled after the end of the run", k)
281 # ---- radiation
282 self.radiation = RadiationModel(self.mesh, self.profile, sh.z0_node, config.listener,
283 self.dt, num.steps_per_block, sh.dtype)
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)
291 # ---- records
292 self.record = bool(out.record_acceleration)
293 self.disp_rate = float(out.displacement_rate)
294 if self.record:
295 nbytes = self.n_out * (2 * self.mesh.n_nodes + 4)
296 log.info("acceleration record: %d samples x %d nodes, int16, %.2f GB in memory",
297 self.n_out, self.mesh.n_nodes, nbytes / 1e9)
298 if nbytes > 2e9:
299 log.warning("the acceleration record needs %.1f GB; consider a shorter duration or "
300 "a coarser mesh", nbytes / 1e9)
301 S = self.fs * self.substeps # sub-steps per second
302 if self.disp_rate > 0:
303 n_frames = int(math.floor(self.n_out * self.substeps * self.disp_rate / S)) + 1
304 self.frame_steps = np.round(np.arange(n_frames) * S / self.disp_rate).astype(np.int64)
305 else:
306 self.frame_steps = np.zeros(0, np.int64)
307 self._build_step()
308 self.setup_time = time.time() - t0
309 log.info("setup done in %.1f s; backend %s", self.setup_time, jax.devices()[0].platform)
311 # ---- one sub-step and its block loop
312 def _build_step(self):
313 sh, sched, rad = self.shell, self.schedule, self.radiation
314 dt, sigma0 = self.dt, self.config.material.sigma0
315 m, W, N = self.substeps, self.span, self.mesh.n_nodes
316 spb = self.config.numerics.steps_per_block
317 tri, normal = sh.tri_j, sh.normal_j
318 record, rate = self.record, self.disp_rate
319 taps_j = jnp.asarray(self.taps, sh.dtype)
320 S = float(self.fs * m)
321 self.acc_rows = spb // m + 2 # audio samples emitted per block, at most
322 self.disp_rows = int(math.ceil(spb * rate / S)) + 2 if rate > 0 else 0
324 n_probes = len(self.probe_xy)
325 p_nodes, p_w = self.probe_nodes, self.probe_w
327 def step(n, k_block0, j_block0, st):
328 q, v, zs, vs, bf, bb, Fmax, ncont, ring, acc_out, disp_out, pbuf = st
329 if rate > 0:
330 # frame k sits at sub-step round(k S / rate); emit q at step n if n is one of them
331 k_cand = jnp.round(n * (rate / S))
332 is_frame = jnp.round(k_cand * (S / rate)) == n
333 row = jnp.clip(k_cand.astype(jnp.int32) - k_block0, 0, self.disp_rows - 1)
334 disp_out = disp_out.at[row].add(jnp.where(is_frame, q[:, 2], 0.0))
335 Fk, F, zs, vs = sched.contact(n, q, zs, vs)
336 fo = sh.internal_force(q)
337 fo = sched.apply_to_plate(fo, Fk)
338 fo = fo.at[:, 2].add(-sh.damping_force_w(v[:, 2]))
339 v1 = (v + dt * fo * sh.Minv_j) / (1 + 2 * sigma0 * dt) # sigma0 implicit
340 v1 = sh.project_rigid(v1)
341 q1 = q + dt * v1
342 zs1, vs1 = sched.advance(zs, vs, F, dt)
343 # normal acceleration at step n: central second difference of q, exact for leapfrog
344 a_node = (((v1 - v) / dt)[:, :3] * normal).sum(1)
345 acc_e = a_node[tri].mean(1)
346 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
349 if record:
350 j0, kk = n // m, n % m
351 emit = kk == 0
352 j_e = j0 - W // 2 # the sample that just became final
353 slot_e = j_e % W
354 done = ring[slot_e]
355 row = jnp.clip(j_e - j_block0, 0, self.acc_rows - 1)
356 acc_out = acc_out.at[row].add(jnp.where(emit, done, 0.0))
357 ring = ring.at[slot_e].set(jnp.where(emit, 0.0, done))
358 # slot s holds sample j0 + 1 - W/2 + w with s = (j0 + 1 - W/2 + w) mod W
359 r = (j0 + 1 - W // 2) % W
360 tw = jnp.take(taps_j[:, kk], (jnp.arange(W) - r) % W)
361 ring = ring + tw[:, None] * a_node[None, :]
362 return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0), ring, acc_out, disp_out, pbuf)
365 def run_block(n0, k_block0, j_block0, st):
366 return lax.fori_loop(0, spb, lambda k, s: step((n0 + k).astype(jnp.int32), k_block0, j_block0, s), st)
369 def energies(st):
370 q, v, _, vs = st[0], st[1], st[2], st[3]
371 Em, Eb = sh.energy_parts_per_element(q)
372 return (sh.kinetic_energy(v), Em.sum(), Eb.sum(), sched.energy(vs),
373 jnp.abs(q[:, 2] / sh.h_node_j).max())
376 def clear(buf, idx):
377 return buf.at[idx].set(0.0)
379 self._step, self._run_block, self._energies, self._clear = step, run_block, energies, clear
381 def initial_state(self):
382 N, dtype = self.mesh.n_nodes, self.shell.dtype
383 zs, vs = self.schedule.initial_state()
384 bf, bb = self.radiation.initial_buffers()
385 K = self.schedule.n_hits
386 ring = jnp.zeros((self.span, N) if self.record else (0,), dtype)
387 return (jnp.zeros((N, N_DOF), dtype), jnp.zeros((N, N_DOF), dtype), zs, vs, bf, bb,
388 jnp.zeros(K, dtype), jnp.zeros(K, jnp.int32), ring, self._empty_outputs()[0],
389 self._empty_outputs()[1], jnp.zeros((self.config.numerics.steps_per_block, len(self.probe_xy)), dtype))
391 def _empty_outputs(self):
392 N, dtype = self.mesh.n_nodes, self.shell.dtype
393 return (jnp.zeros((self.acc_rows, N) if self.record else (0,), dtype),
394 jnp.zeros((self.disp_rows, N) if self.disp_rate > 0 else (0,), dtype))
396 def _hit_summaries(self, Fmax, ncont, vs) -> list[HitSummary]:
397 hits = []
398 for k, p in enumerate(self.patches):
399 hits.append(HitSummary(p.time, p.hit.exciter.kind, float(Fmax[k]), float(ncont[k]) * self.dt,
400 float(vs[k]), p.energy, 0.5 * p.mass * float(vs[k]) ** 2))
401 return hits
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])
406 sim_state = SimulatorState(
407 n_done, np.asarray(q), np.asarray(v), np.asarray(zs), np.asarray(vs), np.asarray(Fmax),
408 np.asarray(ncont), np.asarray(bf), np.asarray(bb), dec_f.tail.copy(), dec_b.tail.copy(),
409 np.asarray(ring) if self.record else None)
410 n = min(len(p_front), len(p_back))
411 p_front, p_back = np.asarray(p_front[:n]), np.asarray(p_back[:n])
412 from . import __version__
413 return SimulationResult(
414 self.config, self.fs, self.dt, self.substeps, self.mesh, self.shell.z0_node, self.shell.h_node,
415 self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front, p_front - p_back,
416 self._hit_summaries(Fmax, ncont, vs), list(energy_log), acc, disp, n_done >= self.n_steps, n_done,
417 sim_state, wall, f"cymsim {__version__}",
418 self.probe_xy if len(self.probe_xy) else None, probe)
420 def run(self, save_path: str | None = None, save_every: float | None = None) -> SimulationResult:
421 """Integrate the run. With ``save_path`` a .cymsim file is written at the end and, if
422 ``save_every`` (seconds of wall time) is given, also periodically during the run, each
423 time as a consistent partial file that can be loaded (see docs/format.md)."""
424 num = self.config.numerics
425 m, spb, N = self.substeps, num.steps_per_block, self.mesh.n_nodes
426 state = self.initial_state()
427 energy_log: list[EnergyRecord] = []
428 dec_f, dec_b = StreamingDecimator(m), StreamingDecimator(m)
429 audio_f, audio_b = [], []
430 n_audio = 0
431 acc = AccelerationRecord.empty(self.fs, N, self.taps) if self.record else None
432 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
437 writer = None
438 if save_path is not None:
439 from .cymsim_file import CymsimWriter
440 writer = CymsimWriter(save_path, self._result(state, [], np.zeros(0), np.zeros(0), acc, disp, 0, 0.0,
441 dec_f, dec_b))
442 t_run = time.time()
443 t_saved = t_run
444 report_every = max(1, self.n_blocks // max(1, num.n_progress))
445 log.info("running %d blocks of %d steps (nonlinearity meter: max|w|/h < 1 effectively "
446 "linear; 3 to 5 clear cascade)", self.n_blocks, spb)
447 for blk in range(self.n_blocks):
448 n0, n1 = blk * spb, (blk + 1) * spb
449 k0 = int(np.searchsorted(self.frame_steps, n0))
450 k1 = int(np.searchsorted(self.frame_steps, n1))
451 j_block0 = -(-n0 // m) - self.span // 2
452 acc_out, disp_out = self._empty_outputs()
453 state = state[:9] + (acc_out, disp_out) + state[11:]
454 state = self._run_block(jnp.int32(n0), jnp.int32(k0), jnp.int32(j_block0), state)
456 # ---- pressure: the ring slots of this block's sub-steps are final; read, clear, decimate
457 slots = jnp.asarray(self.radiation.slots(n0, n1))
458 pf, pb = np.asarray(state[4][slots], np.float64), np.asarray(state[5][slots], np.float64)
459 state = state[:4] + (self._clear(state[4], slots), self._clear(state[5], slots)) + state[6:]
460 new_f, new_b = dec_f.push(pf), dec_b.push(pb)
461 keep = max(0, min(len(new_f), self.n_out - n_audio))
462 audio_f.append(new_f[:keep])
463 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)))
468 n_audio += keep
470 # ---- records
471 if self.record:
472 n_emit = -(-n1 // m) - -(-n0 // m)
473 rows = np.asarray(state[9][:n_emit], np.float32)
474 j = j_block0 + np.arange(n_emit)
475 sel = (j >= 0) & (j < self.n_out)
476 if sel.any():
477 acc.extend(*AccelerationRecord.encode(rows[sel]))
478 if self.disp_rate > 0 and k1 > k0:
479 frames = np.asarray(state[10][:k1 - k0], np.float32)
480 t_frames = self.frame_steps[k0:k1] * self.dt
481 sel = self.frame_steps[k0:k1] <= self.n_out * m
482 disp.extend(t_frames[sel], frames[sel])
484 if blk % report_every == 0 or blk == self.n_blocks - 1:
485 ke, em, eb, ks, woh = (float(x) for x in self._energies(state))
486 t = n1 * self.dt
487 energy_log.append(EnergyRecord(t, ke, em, eb, ks, woh))
488 log.info(" t = %6.3f s KE %.3e membrane %.3e bending %.3e exciters %.3e J "
489 "max|w|/h = %5.2f [%.0f s wall]", t, ke, em, eb, ks, woh, time.time() - t_run)
490 if not math.isfinite(ke + em + eb):
491 raise FloatingPointError("energy is not finite: the run has blown up "
492 "(lower numerics.cfl_safety)")
493 if writer is not None and save_every is not None and time.time() - t_saved >= save_every \
494 and blk < self.n_blocks - 1:
495 res = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),
496 acc, disp, n1, time.time() - t_run, dec_f, dec_b, probe_now())
497 writer.update(res)
498 t_saved = time.time()
499 log.info(" saved partial run (%.3f s) to %s", res.duration_done, save_path)
501 result = self._result(state, energy_log, np.concatenate(audio_f), np.concatenate(audio_b),
502 acc, disp, self.n_steps, time.time() - t_run, dec_f, dec_b, probe_now())
503 for h in result.hits:
504 log.info("hit (%s): peak force %.1f N, contact %.2f ms, rebound velocity %+.2f m/s "
505 "(energy left in exciter %.2f of %.2f J)", h.kind, h.peak_force,
506 h.contact_time * 1e3, h.rebound_velocity, h.energy_out, h.energy_in)
507 for name, p in (("rayleigh", result.pressure_rayleigh), ("unbaffled", result.pressure_unbaffled)):
508 peak = float(np.abs(p).max())
509 log.info("%-10s: peak %.3f Pa (%.1f dB SPL peak)", name, peak, 20 * np.log10(max(peak, 1e-12) / 2e-5))
510 if writer is not None:
511 writer.update(result)
512 writer.close()
513 log.info("wrote %s", save_path)
514 return result
517def simulate(config: SimulationConfig) -> SimulationResult:
518 return Simulator(config).run()