Continue a run from a .cymsim file (cymsim run --resume, simulate(cfg, resume=...)), identical to the uninterrupted run; keep the outputs past the duration in the state and store probe signals in the file - proof-of-concept/cymsim

proof-of-concept / cymsim

Continue a run from a .cymsim file (cymsim run --resume, simulate(cfg, resume=...)), identical to the uninterrupted run; keep the outputs past the duration in the state and store probe signals in the file

Jeremy Magland <jmagland@flatironinstitute.org> committed 4 days ago commit 6c740b5d8b01 parent 42c1bf8 Browse files

7 changed files+307−51

README.mdmodified+3−0View file

@@ -20,6 +20,7 @@ cymsim simulate my_run.yaml # run; writes <prefix>_rayl

2020 cymsim simulate my_run.yaml --linear --duration 0.5 # linear reference run

2121 cymsim run my_run.yaml -o run.cymsim --save-every 60 # also write a .cymsim file, updated every minute of wall time

2222 cymsim run my_run.yaml -o run.cymsim --record-acceleration # include the audio-rate normal acceleration of every node

23+cymsim run my_run.yaml --resume run.cymsim --duration 4 -o run4.cymsim # continue a saved run to 4 s

2324 cymsim info run.cymsim # summarise a .cymsim file

2425 ```

2526

@@ -43,6 +44,8 @@ w = result.displacement.w # transverse displacement at 60 Hz,

4344

4445 from cymsim import load_cymsim

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

47+cfg.numerics.duration = 4.0

48+longer = simulate(cfg, resume=back) # continue it to 4 s; identical to a 4 s run on the same device

4649 ```

4750

4851 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.

docs/format.mdmodified+8−1View file

@@ -104,6 +104,13 @@ Group attributes `fs`, `units = "Pa"`.

104104

105105 Both are the pressure at `listener/position`, sampled at $t_j = j / f_s$, filtered by the same kernel as the acceleration record. In a partial file $T_a$ is the number of samples that were final at the time of writing. The Python package finalises pressure and acceleration samples on the same schedule, so $T_a = T$ in its files, but a reader should not rely on that.

106106

107+`probes` (optional, present when `output.probes` is not empty; attributes `fs`, `units = "m/s^2"`):

108+

109+| dataset | shape | meaning |

110+|---|---|---|

111+| `xy` | $(P, 2)$ float64 | the probe points in the projected plane |

112+| `normal_acceleration` | $(T_p, P)$ float32 | normal acceleration at each point, interpolated linearly in the containing element and filtered by the same kernel as the listener signals |

113+

107114 ## 9. Reproducing the baffled signal from the record

108115

109116 A reader that has only `acceleration`, `mesh` and `listener` recovers the baffled signal as

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

120127

121128 ## 11. Integrator state: `state` (optional, writer-specific)

122129

123-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. The web companion stores its own float32 GPU state under `writer = "cymsimweb"`. Resuming a float32 run reproduces the uninterrupted run only to the precision of the state that was saved.

130+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"`.

124131

125132 ## 12. Reading the file elsewhere

126133

src/cymsim/cli.pymodified+13−1View file

@@ -33,7 +33,16 @@ def cmd_simulate(args) -> int:

3333

3434 cfg = _load(args.config, args)

3535 sim = Simulator(cfg)

36- result = sim.run(save_path=args.output, save_every=args.save_every)

36+ prev = None

37+ if args.resume is not None:

38+ from .cymsim_file import load_cymsim

39+ prev = load_cymsim(args.resume)

40+ try:

41+ sim.check_resumable(prev)

42+ except ValueError as e:

43+ print(f"cannot continue {args.resume}: {e}", file=sys.stderr)

44+ return 2

45+ result = sim.run(save_path=args.output, save_every=args.save_every, resume=prev)

3746 prefix = cfg.output.prefix if cfg.numerics.nonlinear else cfg.output.prefix + "_linear"

3847 if cfg.output.wav:

3948 print("wrote " + ", ".join(result.write_wav(prefix)))

@@ -96,6 +105,9 @@ def main(argv=None) -> int:

96105 help="with -o, also save the partial run every SECONDS of wall time")

97106 p.add_argument("--record-acceleration", action="store_true",

98107 help="keep the normal acceleration of every node at the audio rate (output.record_acceleration)")

108+ p.add_argument("--resume", metavar="FILE.cymsim",

109+ help="continue the run saved in FILE.cymsim (same configuration up to the duration and the "

110+ "output naming) instead of starting from rest")

99111 p.add_argument("--duration", type=float, help="override numerics.duration [s]")

100112 p.add_argument("--out-prefix", help="override output.prefix")

101113 p.add_argument("--linear", action="store_true", help="drop the quadratic strain terms (reference run)")

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

@@ -153,6 +153,13 @@ class Output:

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

154154

155155

156+# The fields in which a continued run may differ from the run it continues: the length of the run

157+# and the naming of the output files. Two configurations that agree on everything else describe

158+# runs that are prefixes of one another (Simulator.run(resume=...)).

159+CONTINUATION_FIELDS = (("numerics", "duration"), ("output", "prefix"), ("output", "wav"), ("output", "npz"),

160+ ("output", "cascade_plot"))

161+

162+

156163 # ----------------------------------------------------------------------------- simulation

157164 @dataclass

158165 class SimulationConfig:

@@ -202,6 +209,13 @@ class SimulationConfig:

202209 def to_dict(self) -> dict:

203210 return _plain(dataclasses.asdict(self))

204211

212+ def identity_dict(self) -> dict:

213+ """``to_dict()`` without CONTINUATION_FIELDS: what a run and its continuation share."""

214+ d = self.to_dict()

215+ for section, name in CONTINUATION_FIELDS:

216+ d[section].pop(name, None)

217+ return d

218+

205219 @classmethod

206220 def from_dict(cls, d: dict) -> "SimulationConfig":

207221 return _from_dict(cls, d)

@@ -217,6 +231,14 @@ class SimulationConfig:

217231

218232

219233 # --------------------------------------------------------------- generic dataclass <-> dict

234+def dict_differences(a, b, path: str = "") -> list[str]:

235+ """Dotted paths at which two plain nested structures differ."""

236+ if isinstance(a, dict) and isinstance(b, dict):

237+ return [d for k in sorted(set(a) | set(b), key=str)

238+ for d in dict_differences(a.get(k), b.get(k), f"{path}.{k}" if path else str(k))]

239+ return [] if a == b else [path or "(root)"]

240+

241+

220242 def _plain(obj):

221243 """Recursively turn tuples into lists so that yaml.safe_dump accepts the structure."""

222244 if isinstance(obj, dict):

src/cymsim/cymsim_file.pymodified+29−10View file

@@ -30,6 +30,11 @@ FORMAT_VERSION = 1

3030 ENERGY_COLUMNS = "time kinetic membrane bending exciters max_w_over_h"

3131 HIT_FIELDS = ("time", "peak_force", "contact_time", "rebound_velocity", "energy_in", "energy_out")

3232 STATE_WRITER = "cymsim-python"

33+STATE_ARRAYS = ("q", "v", "exciter_z", "exciter_v", "peak_force", "contact_steps", "ring_front", "ring_back",

34+ "decimator_tail_front", "decimator_tail_back")

35+OPTIONAL_STATE_ARRAYS = ("acceleration_ring", "decimator_tail_probes", "overhang_front", "overhang_unbaffled",

36+ "overhang_probes", "overhang_acceleration_data", "overhang_acceleration_scale",

37+ "overhang_displacement_time", "overhang_displacement_w")

3338

3439

3540 class CymsimWriter:

@@ -110,6 +115,13 @@ class CymsimWriter:

110115 g.create_dataset("w", shape=(0, N), maxshape=(None, N), dtype=np.float32, chunks=(64, N))

111116 g.create_dataset("time", shape=(0,), maxshape=(None,), dtype=np.float64, chunks=(4096,))

112117

118+ if r.probes is not None:

119+ g = f.create_group("probes")

120+ g.attrs["fs"], g.attrs["units"] = int(r.fs), "m/s^2"

121+ g.create_dataset("xy", data=np.asarray(r.probes, np.float64))

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

123+ dtype=np.float32, chunks=(4096, len(r.probes)))

124+

113125 g = f.create_group("energy")

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

115127 d.attrs["columns"] = ENERGY_COLUMNS

@@ -139,6 +151,8 @@ class CymsimWriter:

139151 if r.displacement is not None and "displacement" in f:

140152 self._append(f["displacement/w"], r.displacement.w)

141153 self._append(f["displacement/time"], r.displacement.time)

154+ if r.probe_acceleration is not None and "probes" in f:

155+ self._append(f["probes/normal_acceleration"], r.probe_acceleration.astype(np.float32))

142156 table = np.array([[e.time, e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h]

143157 for e in r.energy_log], np.float64).reshape(-1, 6)

144158 f["energy/table"].resize(len(table), axis=0)

@@ -155,11 +169,11 @@ class CymsimWriter:

155169 g.attrs["writer"] = STATE_WRITER

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

157171 g.attrs["dtype"] = r.config.numerics.dtype

158- for name in ("q", "v", "exciter_z", "exciter_v", "peak_force", "contact_steps", "ring_front",

159- "ring_back", "decimator_tail_front", "decimator_tail_back"):

172+ for name in STATE_ARRAYS:

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

161- if s.acceleration_ring is not None:

162- g.create_dataset("acceleration_ring", data=np.asarray(s.acceleration_ring))

174+ for name in OPTIONAL_STATE_ARRAYS:

175+ if getattr(s, name) is not None:

176+ g.create_dataset(name, data=np.asarray(getattr(s, name)))

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

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

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

@@ -211,11 +225,12 @@ def load_cymsim(path: str | Path) -> SimulationResult:

211225 state = None

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

213227 g = f["state"]

214- state = SimulatorState(

215- int(g.attrs["n_steps_done"]), g["q"][...], g["v"][...], g["exciter_z"][...], g["exciter_v"][...],

216- g["peak_force"][...], g["contact_steps"][...], g["ring_front"][...], g["ring_back"][...],

217- g["decimator_tail_front"][...], g["decimator_tail_back"][...],

218- g["acceleration_ring"][...] if "acceleration_ring" in g else None)

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

229+ **{name: g[name][...] for name in OPTIONAL_STATE_ARRAYS if name in g})

230+ probes = probe_acc = None

231+ if "probes" in f:

232+ probes = f["probes/xy"][...]

233+ probe_acc = f["probes/normal_acceleration"][...].astype(np.float64)

219234 lg = f["listener"]

220235 creator = f.attrs.get("creator", "")

221236 return SimulationResult(

@@ -223,7 +238,8 @@ def load_cymsim(path: str | Path) -> SimulationResult:

223238 lg["r_front"][...], lg["r_back"][...], f["audio/baffled"][...].astype(np.float64),

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

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

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

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

242+ probes, probe_acc)

227243

228244

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

@@ -246,6 +262,9 @@ def describe(path: str | Path) -> str:

246262 if "displacement" in f:

247263 d = f["displacement/w"]

248264 lines.append(f" displacement: {d.shape[0]} frames at {float(f['displacement'].attrs['rate']):g} Hz")

265+ if "probes" in f:

266+ d = f["probes/normal_acceleration"]

267+ lines.append(f" probes: {d.shape[1]} points, {d.shape[0]} samples")

249268 lines.append(f" energy log: {f['energy/table'].shape[0]} entries; hits: {f['hits/time'].shape[0]}")

250269 if "state" in f:

251270 lines.append(f" integrator state present (writer: {f['state'].attrs.get('writer', '?')}, "

src/cymsim/simulate.pymodified+115−39View file

@@ -31,7 +31,7 @@ import numpy as np

3131 from jax import lax

3232 from scipy.io import wavfile

3333

34-from .config import SimulationConfig

34+from .config import CONTINUATION_FIELDS, SimulationConfig, dict_differences

3535 from .decimation import DECIMATION_SPAN, StreamingDecimator, decimation_taps

3636 from .exciter import HitSchedule, build_patch

3737 from .geometry import CymbalProfile

@@ -147,6 +147,17 @@ class SimulatorState:

147147 decimator_tail_front: np.ndarray # sub-steps still needed by non-final audio samples

148148 decimator_tail_back: np.ndarray

149149 acceleration_ring: np.ndarray | None = None # (W, N) partial audio samples of the record

150+ decimator_tail_probes: np.ndarray | None = None # (P, L) the same for the probe signals

151+ # Final outputs past the end of the duration: the run overshoots it by the kernel's half-span

152+ # and by the rounding to whole blocks. A continued run begins its records with them. None in

153+ # a state written before runs could be continued.

154+ overhang_front: np.ndarray | None = None # (T',) baffled pressure [Pa]

155+ overhang_unbaffled: np.ndarray | None = None # (T',)

156+ overhang_probes: np.ndarray | None = None # (T', P)

157+ overhang_acceleration_data: np.ndarray | None = None # (T', N) int16

158+ overhang_acceleration_scale: np.ndarray | None = None

159+ overhang_displacement_time: np.ndarray | None = None # (Td',)

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

150161

151162

152163 @dataclass

@@ -302,7 +313,8 @@ class Simulator:

302313 "a coarser mesh", nbytes / 1e9)

303314 S = self.fs * self.substeps # sub-steps per second

304315 if self.disp_rate > 0:

305- n_frames = int(math.floor(self.n_out * self.substeps * self.disp_rate / S)) + 1

316+ # every frame the run passes, including those past the duration (the state's overhang)

317+ n_frames = int(math.floor(self.n_steps * self.disp_rate / S)) + 1

306318 self.frame_steps = np.round(np.arange(n_frames) * S / self.disp_rate).astype(np.int64)

307319 else:

308320 self.frame_steps = np.zeros(0, np.int64)

@@ -390,6 +402,34 @@ class Simulator:

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

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

392404

405+ def _state_from(self, s: SimulatorState):

406+ """The device state tuple of ``initial_state`` rebuilt from a saved SimulatorState."""

407+ dtype = self.shell.dtype

408+ f = lambda a: jnp.asarray(a, dtype)

409+ ring = f(s.acceleration_ring) if self.record else jnp.zeros((0,), dtype)

410+ acc_out, disp_out = self._empty_outputs()

411+ return (f(s.q), f(s.v), f(s.exciter_z), f(s.exciter_v), f(s.ring_front), f(s.ring_back),

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

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

414+

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

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

417+ s = prev.state

418+ if s is None:

419+ raise ValueError("the earlier run has no integrator state of this program")

420+ if s.overhang_front is None:

421+ raise ValueError("the earlier run was saved before cymsim could continue runs; run it again")

422+ diffs = dict_differences(prev.config.identity_dict(), self.config.identity_dict())

423+ if diffs:

424+ raise ValueError("the earlier run has a different configuration (" + ", ".join(diffs) + "); only "

425+ + ", ".join(".".join(f) for f in CONTINUATION_FIELDS) + " may change")

426+ if (prev.fs, prev.substeps, prev.dt, prev.n_nodes) != (self.fs, self.substeps, self.dt, self.mesh.n_nodes):

427+ raise ValueError("the earlier run has a different mesh or time step; was cymsim changed since?")

428+ if np.dtype(s.q.dtype) != np.dtype(self.shell.dtype):

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

430+ if s.n_steps_done % self.config.numerics.steps_per_block:

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

432+

393433 def _empty_outputs(self):

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

395435 return (jnp.zeros((self.acc_rows, N) if self.record else (0,), dtype),

@@ -402,51 +442,92 @@ class Simulator:

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

403443 return hits

404444

405- def _result(self, state, energy_log, p_front, p_back, acc, disp, n_done, wall, dec_f, dec_b,

406- probe=None) -> SimulationResult:

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

446+ probe) -> SimulationResult:

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

448+ moved to the state's overhang."""

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

450+ n_out = self.n_out

451+ n = min(len(p_front), len(p_unb))

452+ p_front, p_unb = np.asarray(p_front[:n]), np.asarray(p_unb[:n])

453+ acc_over = (None, None)

454+ if acc is not None:

455+ acc_over = (acc.data[n_out:], acc.scale[n_out:])

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

457+ disp_over = (None, None)

458+ if disp is not None:

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

460+ disp_over = (disp.time[k_out:], disp.w[k_out:])

461+ disp = DisplacementRecord(disp.rate, disp.time[:k_out], disp.w[:k_out])

408462 sim_state = SimulatorState(

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

410464 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])

465+ np.asarray(ring) if self.record else None,

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

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

414468 from . import __version__

415469 return SimulationResult(

416470 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,

471+ self.shell.normal_node, self.radiation.r_front, self.radiation.r_back, p_front[:n_out], p_unb[:n_out],

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

419473 sim_state, wall, f"cymsim {__version__}",

420- self.probe_xy if len(self.probe_xy) else None, probe)

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

421475

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

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

477+ resume: SimulationResult | None = None) -> SimulationResult:

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

424479 ``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)."""

480+ time as a consistent partial file that can be loaded (see docs/format.md).

481+

482+ With ``resume``, an earlier run (usually loaded from a .cymsim file) of a configuration

483+ that differs from this one at most in CONTINUATION_FIELDS, the run continues from that

484+ run's integrator state, and its records begin with the earlier ones. On the same device

485+ the result is identical to that of an uninterrupted run. If the earlier run already

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

426487 num = self.config.numerics

427488 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

489+ n_probes = len(self.probe_xy)

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

434491 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 = []

492+ if resume is None:

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

494+ energy_log: list[EnergyRecord] = []

495+ audio_f, audio_u, audio_p = [np.zeros(0)], [np.zeros(0)], [np.zeros((0, n_probes))]

496+ dec_f, dec_b = StreamingDecimator(m), StreamingDecimator(m)

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

498+ else:

499+ self.check_resumable(resume)

500+ s = resume.state

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

502+ energy_log = list(resume.energy_log)

503+ audio_f = [resume.pressure_rayleigh, s.overhang_front]

504+ audio_u = [resume.pressure_unbaffled, s.overhang_unbaffled]

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

506+ dec_f = StreamingDecimator(m, tail=s.decimator_tail_front, n_seen=n_done)

507+ dec_b = StreamingDecimator(m, tail=s.decimator_tail_back, n_seen=n_done)

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

509+ if n_probes else []

510+ if acc is not None:

511+ acc.extend(resume.acceleration.data, resume.acceleration.scale)

512+ acc.extend(s.overhang_acceleration_data, s.overhang_acceleration_scale)

513+ if disp is not None:

514+ disp.extend(resume.displacement.time, resume.displacement.w)

515+ disp.extend(s.overhang_displacement_time, s.overhang_displacement_w)

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

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

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

519+ acc, disp, n, wall, dec_f, dec_b, dec_p, probe_now())

439520 writer = None

440521 if save_path is not None:

441522 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))

523+ writer = CymsimWriter(save_path, result_now(n_done, wall0))

444524 t_run = time.time()

445525 t_saved = t_run

446526 report_every = max(1, self.n_blocks // max(1, num.n_progress))

527+ blk0 = n_done // spb

447528 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):

529+ "linear; 3 to 5 clear cascade)", max(0, self.n_blocks - blk0), spb)

530+ for blk in range(blk0, self.n_blocks):

450531 n0, n1 = blk * spb, (blk + 1) * spb

451532 k0 = int(np.searchsorted(self.frame_steps, n0))

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

@@ -459,29 +540,26 @@ class Simulator:

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

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

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

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

462544 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])

545+ audio_f.append(new_f)

546+ audio_u.append(new_f - new_b)

466547 if n_probes:

467548 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

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

550+ n_done = n1

471551

472552 # ---- records

473553 if self.record:

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

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

476556 j = j_block0 + np.arange(n_emit)

477- sel = (j >= 0) & (j < self.n_out)

557+ sel = j >= 0

478558 if sel.any():

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

480560 if self.disp_rate > 0 and k1 > k0:

481561 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])

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

485563

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

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

@@ -494,14 +572,12 @@ class Simulator:

494572 "(lower numerics.cfl_safety)")

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

496574 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())

575+ res = result_now(n_done, wall0 + time.time() - t_run)

499576 writer.update(res)

500577 t_saved = time.time()

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

502579

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())

580+ result = result_now(n_done, wall0 + time.time() - t_run)

505581 for h in result.hits:

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

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

@@ -516,5 +592,5 @@ class Simulator:

516592 return result

517593

518594

519-def simulate(config: SimulationConfig) -> SimulationResult:

520- return Simulator(config).run()

595+def simulate(config: SimulationConfig, resume: SimulationResult | None = None) -> SimulationResult:

596+ return Simulator(config).run(resume=resume)

tests/test_resume.pyadded+117−0View file

@@ -0,0 +1,117 @@

1+"""Continuing a run from a .cymsim file reproduces the uninterrupted run exactly."""

2+import copy

3+

4+import h5py

5+import numpy as np

6+import pytest

7+

8+from cymsim.cli import main

9+from cymsim.cymsim_file import CymsimWriter, load_cymsim

10+from cymsim.simulate import Simulator

11+

12+

13+def _config(base, duration):

14+ cfg = copy.deepcopy(base)

15+ cfg.numerics.duration = duration

16+ cfg.output.record_acceleration = True

17+ cfg.output.displacement_rate = 200.0

18+ cfg.output.probes = [(0.1, 0.05), (0.05, -0.12)]

19+ return cfg

20+

21+

22+def _datasets(path) -> dict:

23+ """Every dataset of a .cymsim file except the energy log, whose report times depend on how

24+ the run was split."""

25+ out = {}

26+ with h5py.File(path) as f:

27+ f.visititems(lambda name, obj: out.__setitem__(name, obj[()]) if isinstance(obj, h5py.Dataset) else None)

28+ out.pop("energy/table")

29+ return out

30+

31+

32+def _assert_same_files(a, b):

33+ da, db = _datasets(a), _datasets(b)

34+ assert sorted(da) == sorted(db)

35+ for name in da:

36+ assert np.array_equal(da[name], db[name]), name

37+ with h5py.File(a) as fa, h5py.File(b) as fb:

38+ for key in ("n_steps_done", "complete", "duration", "n_samples"):

39+ assert fa["run"].attrs[key] == fb["run"].attrs[key], key

40+

41+

42+@pytest.fixture(scope="module")

43+def straight(coarse_config, tmp_path_factory):

44+ path = tmp_path_factory.mktemp("straight") / "straight.cymsim"

45+ Simulator(_config(coarse_config, 0.05)).run(save_path=str(path))

46+ return path

47+

48+

49+def test_continued_run_equals_the_uninterrupted_one(coarse_config, straight, tmp_path):

50+ prev = None

51+ for k, duration in enumerate((0.02, 0.035, 0.05)):

52+ path = tmp_path / f"seg{k}.cymsim"

53+ Simulator(_config(coarse_config, duration)).run(save_path=str(path), resume=prev)

54+ prev = load_cymsim(path)

55+ assert prev.complete and len(prev.pressure_rayleigh) == round(44100 * duration)

56+ assert len(prev.state.overhang_front) > 0 # the overshoot past the duration is kept

57+ _assert_same_files(path, straight)

58+ assert len(prev.energy_log) > len(load_cymsim(straight).energy_log)

59+

60+

61+def test_continue_from_a_killed_run(coarse_config, straight, tmp_path, monkeypatch):

62+ """A run killed right after a partial save continues from the partial file."""

63+ path = tmp_path / "killed.cymsim"

64+ update, calls = CymsimWriter.update, []

65+

66+ def update_then_die(self, result):

67+ update(self, result)

68+ calls.append(1)

69+ if len(calls) == 3:

70+ self.close()

71+ raise KeyboardInterrupt

72+

73+ monkeypatch.setattr(CymsimWriter, "update", update_then_die)

74+ with pytest.raises(KeyboardInterrupt):

75+ Simulator(_config(coarse_config, 0.05)).run(save_path=str(path), save_every=0)

76+ monkeypatch.undo()

77+ part = load_cymsim(path)

78+ assert not part.complete and 0 < part.n_steps_done

79+ out = tmp_path / "continued.cymsim"

80+ Simulator(_config(coarse_config, 0.05)).run(save_path=str(out), resume=part)

81+ _assert_same_files(out, straight)

82+

83+

84+def test_shorter_duration_cuts_the_records(coarse_config, straight, tmp_path):

85+ short = Simulator(_config(coarse_config, 0.03)).run()

86+ cut = Simulator(_config(coarse_config, 0.03)).run(resume=load_cymsim(straight))

87+ n = round(44100 * 0.03)

88+ assert len(cut.pressure_rayleigh) == n and cut.complete

89+ assert np.array_equal(cut.pressure_rayleigh.astype(np.float32), short.pressure_rayleigh.astype(np.float32))

90+ assert np.array_equal(cut.acceleration.data, short.acceleration.data)

91+ assert np.array_equal(cut.displacement.w, short.displacement.w)

92+ assert np.array_equal(cut.probe_acceleration.astype(np.float32), short.probe_acceleration.astype(np.float32))

93+

94+

95+def test_a_different_configuration_is_refused(coarse_config, straight):

96+ cfg = _config(coarse_config, 0.08)

97+ cfg.hits[0].velocity *= 1.001

98+ with pytest.raises(ValueError, match=r"hits"):

99+ Simulator(cfg).run(resume=load_cymsim(straight))

100+ cfg = _config(coarse_config, 0.08)

101+ cfg.output.prefix = "elsewhere" # output naming may change

102+ Simulator(cfg).check_resumable(load_cymsim(straight))

103+

104+

105+def test_cli_resume(coarse_config, straight, tmp_path, capsys):

106+ cfg = _config(coarse_config, 0.05)

107+ cfg.output.wav = cfg.output.npz = cfg.output.cascade_plot = False

108+ cfg.output.prefix = str(tmp_path / "x")

109+ ypath = tmp_path / "c.yaml"

110+ cfg.to_yaml(ypath)

111+ out = tmp_path / "longer.cymsim"

112+ assert main(["run", str(ypath), "--resume", str(straight), "--duration", "0.06", "-o", str(out), "--quiet"]) == 0

113+ assert load_cymsim(out).complete and len(load_cymsim(out).pressure_rayleigh) == round(44100 * 0.06)

114+ cfg.hits[0].r += 0.01

115+ cfg.to_yaml(ypath)

116+ assert main(["run", str(ypath), "--resume", str(straight), "-o", str(out), "--quiet"]) == 2

117+ assert "hits" in capsys.readouterr().err