af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 1"""Von Karman (Airy stress function) solution of cymbal-flat-disk-v1 by a modal method on a Cartesian grid.
3The disk is approximated by the square cells of a uniform Cartesian grid whose centers lie inside
4it (see bfs_disk.py). The transverse displacement is expanded in the linear free-plate modes Phi_k
5of that staircase domain, computed with Bogner-Fox-Schmit elements, and the Airy stress function
6in the clamped biharmonic eigenfunctions Psi_p on the same mesh. With mass-normalized Phi and
7L2-normalized Psi the equations are
9 q_k'' + 2 alpha_k q_k' + omega_k^2 q_k = f_k psi(t/T_s) + sum_{j,p} H^p_{kj} q_j eta_p,
10 eta_p = -(E h / (2 zeta_p^4)) sum_{i,j} H^p_{ij} q_i q_j, H^p_{ij} = int Psi_p L(Phi_i, Phi_j) dA,
12where L(f, g) = f_xx g_yy + f_yy g_xx - 2 f_xy g_xy (Touze et al.). Time stepping uses a two-step
13Gautschi-type scheme at a step of 1/(fs m): the damped linear part is propagated exactly and the
14strike and the cubic coupling force are applied with a static-exact weight, which is second-order
15accurate and free of numerical dispersion in the linear limit. The Rayleigh integral is evaluated
16with 3 x 3 Gauss points per cell; each point's retarded acceleration is a cubic Lagrange
17interpolation of the modal accelerations at neighbouring time steps, folded into a per-step
18weight table. Damping acts on the bending motion only (the in-plane state is quasi-static),
19a further approximation of the challenge model.
21The convergence chain has three axes: the grid spacing (8 mm, 4 mm), the number of modes kept
22(transverse and Airy), and the time step. The mode truncation caps the frequency the nonlinear
23cascade can reach.
25Because the strike is centered on the x axis, the solution is even in y, and by default the modal bases
26(transverse and Airy) are restricted to that symmetric subspace, which is exact and halves the number of
27modes needed for a given frequency reach (the coupling cost scales with the cube of the count). --full
28uses the unrestricted bases instead.
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 30With --device cuda the coupling tensor is assembled and the time stepping runs on a GPU through PyTorch
31(the modes and quadrature stay in numpy/scipy); the CPU path is the reference and the two agree to round-off.
f5791abStrike duration override for previewsJeremy Magland 33Usage: python vk_modal_disk.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [--duration T] [--seconds S] [--full] [--device cuda] [tag ...]
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 34 where a tag is g<grid spacing in mm>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. g8-w200-a200-m1
35"""
36from __future__ import annotations
38import json
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 39import os
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 40import sys
41import time
42from pathlib import Path
44import numpy as np
45from scipy.io import wavfile
46sys.path.insert(0, str(Path(__file__).resolve().parent))
47from bfs_disk import BFSDisk, load_challenge # noqa: E402
50class VonKarmanDisk:
f5791abStrike duration override for previewsJeremy Magland 51 def __init__(self, ch: dict, spacing: float, n_w: int, n_airy: int, force_override: float | None = None, log=None, symmetric_y: bool | None = None, device: str = "cpu", duration_override: float | None = None):
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 52 self.ch = ch
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 53 self.device = device
54 self.torch = None
55 if device != "cpu":
56 import torch
58 self.torch = torch
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 59 p = dict(ch["parameters"])
60 if force_override is not None:
61 p["strike_peak_force"] = force_override
63 p["strike_duration"] = duration_override
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 64 ch2 = dict(ch)
65 ch2["parameters"] = p
66 self.p = p
67 self.log = log or (lambda s: None)
68 if symmetric_y is None:
69 symmetric_y = p["strike_y"] == 0.0
70 self.symmetric_y = symmetric_y
71 t0 = time.perf_counter()
72 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
73 pl = self.plate
74 self.log(f" linear modes: {pl.n_modes}{' (even in y)' if symmetric_y else ''} on {pl.n_cells} cells ({pl.ndof} dof) in {time.perf_counter() - t0:.1f}s, f_max={pl.omega[-1] / (2 * np.pi):.0f} Hz, rigid={pl.n_rigid}")
75 self.n_w = pl.n_modes
76 self.omega = pl.omega
77 self.alpha = pl.alpha
78 self.f = pl.f
79 V = pl.V
80 # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
81 t0 = time.perf_counter()
82 Kb = pl.assemble(pl.element_matrix(np.array([[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]])))
83 M1 = pl.M / pl.m
84 z4, Psi = pl.eigenmodes(Kb, M1, n_airy, sigma=0.0, dofs=pl.interior_dofs)
85 self.n_airy = len(z4)
86 self.zeta4 = z4
87 self.log(f" airy modes: {self.n_airy} in {time.perf_counter() - t0:.1f}s")
88 # coupling tensor H[p, i, j] by 3 x 3 Gauss quadrature per cell
89 t0 = time.perf_counter()
90 xq, yq, WQ = pl.quadrature(3)
91 Bxx, Byy, Bxy = pl.basis_derivs_at(xq, yq)
92 Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
93 Bn = pl.basis_at(xq, yq)
94 PsiQ = Bn @ Psi # (Q, n_airy)
95 self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 96 if self.torch is None:
97 H = np.zeros((self.n_airy, self.n_w, self.n_w))
98 for pi in range(self.n_airy):
99 wp = (WQ * PsiQ[:, pi])[:, None]
100 A = (Pxx * wp).T @ Pyy
101 H[pi] = A + A.T - 2.0 * (Pxy * wp).T @ Pxy
102 self.H = H
103 self.Hflat = H.reshape(self.n_airy * self.n_w, self.n_w)
104 shape = H.shape
105 else:
106 # the same assembly on the device; the tensor stays there for the time stepping
107 t = self.torch
108 dev = t.device(self.device)
109 Pxx_t, Pyy_t, Pxy_t = (t.as_tensor(a, device=dev) for a in (Pxx, Pyy, Pxy))
110 WPsi = t.as_tensor(WQ[:, None] * PsiQ, device=dev) # (Q, n_airy)
111 H = t.empty((self.n_airy, self.n_w, self.n_w), dtype=t.float64, device=dev)
112 for pi in range(self.n_airy):
113 wp = WPsi[:, pi : pi + 1]
114 A = (Pxx_t * wp).T @ Pyy_t
115 H[pi] = A + A.T - 2.0 * (Pxy_t * wp).T @ Pxy_t
116 t.cuda.synchronize() if self.device.startswith("cuda") else None
117 self.H = None
118 self.Hflat_t = H.reshape(self.n_airy * self.n_w, self.n_w)
119 shape = tuple(H.shape)
120 del Pxx_t, Pyy_t, Pxy_t, WPsi
121 self.log(f" coupling tensor {shape} on {self.device} in {time.perf_counter() - t0:.1f}s")
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 122 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
123 _, _, self.tau, self.cq = pl.listener_quadrature(3)
124 self.PhiQ = Bn @ V # (Q, n_w)
126 def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
127 v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
128 eta = self.coef * (v @ q) # eta_p
129 return v.T @ eta # F_k = sum_p v[p,k] eta_p
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 131 def _prepare(self, fs: int, nsamp: int, m: int, field_grid, field_rate):
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 132 p = self.p
133 dt = 1.0 / (fs * m)
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 134 w2, al = self.omega**2, self.alpha
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 135 # exact two-step recurrence of the damped linear oscillator: q+ = a1 q + a2 q- + b F
136 rigid = self.omega == 0.0
137 disc = np.sqrt(al**2 - w2 + 0j)
138 E1, E2 = np.exp((-al + disc) * dt), np.exp((-al - disc) * dt)
139 a1 = (E1 + E2).real
140 a2 = -(E1 * E2).real
141 with np.errstate(divide="ignore", invalid="ignore"):
142 b = np.where(rigid, 0.0, (1.0 - a1 - a2) / np.where(rigid, 1.0, w2))
143 if rigid.any():
144 e2 = np.exp(-2.0 * al[rigid] * dt)
145 a1[rigid] = 1.0 + e2
146 a2[rigid] = -e2
147 b[rigid] = dt * (1.0 - e2) / (2.0 * al[rigid])
148 # fractional delays folded into a per-step weight table G[step, mode]
149 tau_steps = self.tau / dt
150 d0 = np.floor(tau_steps).astype(int)
151 fr = tau_steps - d0
152 lag = np.stack([
153 -fr * (fr - 1.0) * (fr - 2.0) / 6.0,
154 (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0,
155 -(fr + 1.0) * fr * (fr - 2.0) / 2.0,
156 (fr + 1.0) * fr * (fr - 1.0) / 6.0,
157 ], axis=1)
158 dmax = int(d0.max()) + 2
159 G = np.zeros((dmax + 1, self.n_w))
160 for j, off in enumerate([-1, 0, 1, 2]):
161 dd = d0 + off
162 for u in np.unique(dd):
163 sel = dd == u
164 G[u] += (self.cq[sel] * lag[sel, j]) @ self.PhiQ[sel]
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 165 fld = None
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 166 if field_grid is not None:
167 fx, fy, ft = field_grid
168 XX, YY = np.meshgrid(fx, fy, indexing="xy")
169 inside = self.plate.inside(XX.ravel(), YY.ravel())
170 PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.plate.V
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 171 fld = dict(PhiG=PhiG, inside=inside, every=int(round(fs * m / field_rate)), nt=len(ft), ny=len(fy), nx=len(fx))
172 return dict(dt=dt, T=p["strike_duration"], Om=2 * np.pi / p["strike_duration"], nsteps=nsamp * m, a1=a1, a2=a2, b=b, w2=w2, al=al, G=G, dmax=dmax, fld=fld)
174 def run(self, fs: int, nsamp: int, m: int, field_grid=None, field_rate: float | None = None):
175 pr = self._prepare(fs, nsamp, m, field_grid, field_rate)
176 if self.torch is not None:
177 return self._run_device(pr, m, nsamp)
178 dt, T, Om, nsteps, a1, a2, b, w2, al, G, dmax, fld = (pr[k] for k in ("dt", "T", "Om", "nsteps", "a1", "a2", "b", "w2", "al", "G", "dmax", "fld"))
179 f = self.f
180 ring = np.zeros((dmax + 1, self.n_w))
181 out = np.zeros(nsamp)
182 q_prev = np.zeros(self.n_w)
183 q = np.zeros(self.n_w)
184 frames = [] if fld else None
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 185 t_start = time.perf_counter()
186 report = max(1, nsteps // 10)
187 for nstep in range(nsteps + 1):
188 t = nstep * dt
189 psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
190 Fext = self.nonlinear_force(q) + f * psi
191 q_next = a1 * q + a2 * q_prev + b * Fext
192 qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
193 ring[nstep % (dmax + 1)] = qdd
194 if nstep % m == 0:
195 k = nstep // m
196 if k < nsamp:
197 idx = (nstep - np.arange(dmax + 1)) % (dmax + 1)
198 out[k] = np.sum(G * ring[idx])
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 199 if frames is not None and nstep % fld["every"] == 0 and len(frames) < fld["nt"]:
200 fr_ = (fld["PhiG"] @ q).astype(np.float32)
201 fr_[~fld["inside"]] = np.nan
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 202 frames.append(fr_)
203 q_prev, q = q, q_next
204 if nstep % report == 0 and nstep:
205 if not np.all(np.isfinite(q)):
206 raise RuntimeError(f"solution diverged at step {nstep} (t = {t:.4f} s); increase the oversampling m")
207 el = time.perf_counter() - t_start
208 self.log(f" step {nstep}/{nsteps} t={t:.2f}s elapsed={el:.0f}s eta={el * (nsteps - nstep) / nstep:.0f}s max|q|={np.abs(q).max():.3e}")
209 if not np.all(np.isfinite(out)):
210 raise RuntimeError("output is not finite; increase the oversampling m")
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 211 fields = np.stack(frames).reshape(len(frames), fld["ny"], fld["nx"]) if frames is not None else None
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 212 return out, fields
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 214 def _run_device(self, pr, m: int, nsamp: int):
215 """The same loop with the tensor contraction and the recurrence on the device (PyTorch)."""
216 t = self.torch
217 dev = t.device(self.device)
218 T64 = lambda a: t.as_tensor(np.asarray(a, dtype=np.float64), device=dev)
219 dt, T, Om, nsteps, dmax, fld = pr["dt"], pr["T"], pr["Om"], pr["nsteps"], pr["dmax"], pr["fld"]
220 a1, a2, b, w2, al, G, f, coef = (T64(x) for x in (pr["a1"], pr["a2"], pr["b"], pr["w2"], pr["al"], pr["G"], self.f, self.coef))
221 Hf = self.Hflat_t
222 L = dmax + 1
223 ar = t.arange(L, device=dev)
224 ring = t.zeros((L, self.n_w), dtype=t.float64, device=dev)
225 out = t.zeros(nsamp, dtype=t.float64, device=dev)
226 q = t.zeros(self.n_w, dtype=t.float64, device=dev)
227 q_prev = t.zeros_like(q)
228 frames = [] if fld else None
229 PhiG = T64(fld["PhiG"]) if fld else None
5b29a8cFix the damping term in the device loopJeremy Magland 230 invdt = 1.0 / dt # 2 alpha (q+ - q-) / (2 dt)
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 231 t_start = time.perf_counter()
232 report = max(1, nsteps // 10)
233 for nstep in range(nsteps + 1):
234 tt = nstep * dt
235 psi = 0.5 * (1.0 - np.cos(Om * tt)) if 0.0 <= tt <= T else 0.0
236 v = (Hf @ q).reshape(self.n_airy, self.n_w)
237 eta = coef * (v @ q)
238 Fext = v.T @ eta
239 if psi:
240 Fext = Fext + f * psi
241 q_next = a1 * q + a2 * q_prev + b * Fext
5b29a8cFix the damping term in the device loopJeremy Magland 242 qdd = Fext - w2 * q - al * ((q_next - q_prev) * invdt)
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 243 ring[nstep % L] = qdd
244 if nstep % m == 0:
245 k = nstep // m
246 if k < nsamp:
247 idx = (nstep - ar) % L
248 out[k] = t.sum(G * ring[idx])
249 if frames is not None and nstep % fld["every"] == 0 and len(frames) < fld["nt"]:
250 fr_ = (PhiG @ q).cpu().numpy().astype(np.float32)
251 fr_[~fld["inside"]] = np.nan
252 frames.append(fr_)
253 q_prev, q = q, q_next
254 if nstep % report == 0 and nstep:
255 qc = q.cpu().numpy()
256 if not np.all(np.isfinite(qc)):
257 raise RuntimeError(f"solution diverged at step {nstep} (t = {tt:.4f} s); increase the oversampling m")
258 el = time.perf_counter() - t_start
259 self.log(f" step {nstep}/{nsteps} t={tt:.2f}s elapsed={el:.0f}s eta={el * (nsteps - nstep) / nstep:.0f}s max|q|={np.abs(qc).max():.3e}")
260 out = out.cpu().numpy()
261 if not np.all(np.isfinite(out)):
262 raise RuntimeError("output is not finite; increase the oversampling m")
263 fields = np.stack(frames).reshape(len(frames), fld["ny"], fld["nx"]) if frames is not None else None
264 return out, fields
267def hardware_string(device: str) -> str:
268 cpu = "unknown CPU"
269 try:
270 for line in open("/proc/cpuinfo"):
271 if line.startswith("model name"):
272 cpu = " ".join(line.split(":", 1)[1].split())
273 break
274 except OSError:
275 pass
276 s = f"{cpu} ({os.cpu_count()} threads), numpy/scipy for the setup"
277 if device != "cpu":
278 import torch
280 s += f"; coupling tensor and time stepping on {torch.cuda.get_device_name(0)} (PyTorch {torch.__version__})"
281 else:
282 s += " and the time stepping"
283 return s
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 285
286def parse_tag(tag: str):
287 parts = dict((kv[0], int(kv[1:])) for kv in tag.split("-"))
288 return parts["g"] * 1e-3, parts["w"], parts["a"], parts["m"]
291def main():
292 ch = load_challenge(Path(sys.argv[1]))
293 out_dir = Path(sys.argv[2])
294 args = sys.argv[3:]
295 fields_for = set()
296 force = None
297 seconds = None
299 def take(flag, conv):
300 nonlocal args
301 if flag in args:
302 i = args.index(flag)
303 v = conv(args[i + 1])
304 args = args[:i] + args[i + 2:]
305 return v
306 return None
308 v = take("--fields", lambda s: set(s.split(",")))
309 fields_for = v or set()
310 force = take("--force", float)
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 312 seconds = take("--seconds", float)
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 313 device = take("--device", str) or "cpu"
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 314 full = "--full" in args
315 if full:
316 args.remove("--full")
317 tags = args or ["g8-w150-a150-m1", "g8-w300-a300-m1", "g4-w300-a300-m2", "g4-w400-a400-m2"]
318 out = ch["output"]
319 fs = out["sample_rate_hz"]
320 nsamp = out["num_samples"] if seconds is None else int(round(fs * seconds))
321 log = lambda s: print(s, file=sys.stderr, flush=True)
322 meta = []
323 for tag in tags:
324 spacing, n_w, n_a, m = parse_tag(tag)
f5791abStrike duration override for previewsJeremy Magland 325 log(f"{tag}: force={force if force is not None else ch['parameters']['strike_peak_force']} N, duration={duration if duration is not None else ch['parameters']['strike_duration']} s")
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 326 t0 = time.perf_counter()
f5791abStrike duration override for previewsJeremy Magland 327 model = VonKarmanDisk(ch, spacing, n_w, n_a, force, log=log, symmetric_y=False if full else None, device=device, duration_override=duration)
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 328 t1 = time.perf_counter()
329 grid = None
330 if tag in fields_for:
331 fg = ch["fields"]
332 nt, ny, nx = fg["shape"]
333 R = ch["parameters"]["disk_radius"]
334 grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
335 s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
336 dt_total = time.perf_counter() - t0
337 sid = f"vk-{tag}"
338 out_dir.mkdir(parents=True, exist_ok=True)
339 wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
340 entry = {
341 "solution_id": sid,
342 "output": f"{out_dir.name}/{sid}.wav",
343 "method": "von Karman plate (Airy form) on a Cartesian grid masked to the disk, modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
344 "parameters": {
345 "grid_spacing_m": spacing,
346 "cells": model.plate.n_cells,
347 "n_dof": model.plate.ndof,
348 "transverse_modes": model.n_w,
349 "basis_even_in_y": model.symmetric_y,
350 "airy_modes": model.n_airy,
351 "oversampling": m,
352 "dt_s": 1.0 / (fs * m),
353 "highest_transverse_mode_hz": round(float(model.omega.max() / (2 * np.pi)), 1),
354 },
355 "timing_seconds": round(dt_total, 3),
f827f8eGPU path for the tensor assembly and time stepping (PyTorch); record the hardware per solutionJeremy Magland 356 "hardware": hardware_string(device),
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 357 "peak_si": float(np.max(np.abs(s))),
358 "setup_seconds": round(t1 - t0, 3),
359 }
360 if fields is not None:
361 import h5py
363 fg = ch["fields"]
364 with h5py.File(out_dir / f"{sid}.h5", "w") as fh:
365 fh.attrs["challenge_id"] = ch["challenge_id"]
366 fh.attrs["solution_id"] = sid
367 fh.attrs["quantity"] = fg["quantity"]
368 fh.attrs["units"] = fg["units"]
369 fh.create_dataset("x", data=grid[0])
370 fh.create_dataset("y", data=grid[1])
371 fh.create_dataset("t", data=grid[2][: fields.shape[0]])
372 fh.create_dataset("u", data=fields, compression="gzip", compression_opts=4, chunks=(8, fields.shape[1], fields.shape[2]))
373 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
374 meta.append(entry)
375 log(f"{tag}: done, setup={t1 - t0:.1f}s total={dt_total:.1f}s peak={np.max(np.abs(s)):.4e} Pa")
376 with open(out_dir / f"{sid}.meta.json", "w") as fh:
377 json.dump(entry, fh, indent=1)
378 print(json.dumps(meta, indent=1))
381if __name__ == "__main__":
382 main()