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.
30If the challenge has a shell-of-revolution geometry (cymbal-bow-v1: a paraboloidal bow z = Z(r)), the same
31solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
32system: L(w + w0, Phi) on the right of the bending equation and 2 L(w0, w) in the compatibility equation.
33With G^p_j = int Psi_p L(w0, Phi_j) dA this gives a linear membrane stiffening 2 sum_p |c_p| G^p (G^p)^T,
34which is diagonalized once in the flat-mode basis to give the shell modes (exact in the Gautschi step,
35damped at sigma0 + eta omega^2 / 2 like every mode, as the Kelvin-Voigt law of the challenge implies), and
36quadratic coupling terms F_k = sum_p c_p [ (r^T H^p r + 2 G^p . r) (H^p r)_k + (r^T H^p r) G^p_k ]. The strike
37acts along -z on the curved top face with the normalization of the specification, the listener integral uses
38the distance to the top face, and the field export is the normal displacement n_z w. Shallow-shell theory
39neglects in-plane inertia and terms of the order of the slope squared (0.06 at the rim here).
41With --device cuda the coupling tensor is assembled and the time stepping runs on a GPU through PyTorch
42(the modes and quadrature stay in numpy/scipy); the CPU path is the reference and the two agree to round-off.
44Usage: python vk_modal_disk.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [--duration T] [--seconds S] [--full] [--device cuda] [tag ...]
45 where a tag is g<grid spacing in mm>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. g8-w200-a200-m1
46"""
47from __future__ import annotations
49import json
50import os
51import sys
52import time
53from pathlib import Path
55import numpy as np
56from scipy.io import wavfile
57sys.path.insert(0, str(Path(__file__).resolve().parent))
58from bfs_disk import BFSDisk, load_challenge # noqa: E402
61class Paraboloid:
62 """The bow of cymbal-bow-v1: mid-surface z = Z(r) = H (1 - r^2 / R^2), uniform thickness along the normal."""
64 def __init__(self, R: float, H: float, h: float):
65 self.R, self.H, self.h = R, H, h
68 def from_challenge(ch: dict):
69 g = ch.get("geometry")
70 if not g or g.get("type") != "shell-of-revolution":
71 return None
72 p = ch["parameters"]
73 return Paraboloid(p["disk_radius"], p["bow_height"], p["h"])
75 def z(self, x, y):
76 return self.H * (1.0 - (x * x + y * y) / self.R**2)
78 def grad(self, x, y):
79 c = -2.0 * self.H / self.R**2
80 return c * x, c * y
82 def curvatures(self, x, y):
83 """Second derivatives (Z_xx, Z_yy, Z_xy) of the mid-surface height."""
84 c = -2.0 * self.H / self.R**2
85 return c * np.ones_like(x), c * np.ones_like(x), np.zeros_like(x)
87 def normal(self, x, y):
88 zx, zy = self.grad(x, y)
89 nz = 1.0 / np.sqrt(1.0 + zx * zx + zy * zy)
90 return -zx * nz, -zy * nz, nz
92 def top(self, x, y):
93 """The point of the top face on the normal through the mid-surface point (x, y)."""
94 nx, ny, nz = self.normal(x, y)
95 return x + 0.5 * self.h * nx, y + 0.5 * self.h * ny, self.z(x, y) + 0.5 * self.h * nz
98class VonKarmanDisk:
99 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):
100 self.ch = ch
101 self.device = device
102 self.torch = None
103 if device != "cpu":
104 import torch
106 self.torch = torch
107 p = dict(ch["parameters"])
108 if force_override is not None:
109 p["strike_peak_force"] = force_override
110 if duration_override is not None:
111 p["strike_duration"] = duration_override
112 ch2 = dict(ch)
113 ch2["parameters"] = p
114 self.p = p
115 self.log = log or (lambda s: None)
116 if symmetric_y is None:
117 symmetric_y = p["strike_y"] == 0.0
118 self.symmetric_y = symmetric_y
119 self.profile = Paraboloid.from_challenge(ch)
120 t0 = time.perf_counter()
121 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
122 pl = self.plate
123 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}")
124 self.n_w = pl.n_modes
125 self.omega = pl.omega
126 self.alpha = pl.alpha
127 self.f = pl.f
128 V = pl.V
129 # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
130 t0 = time.perf_counter()
131 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]])))
132 M1 = pl.M / pl.m
133 z4, Psi = pl.eigenmodes(Kb, M1, n_airy, sigma=0.0, dofs=pl.interior_dofs)
134 self.n_airy = len(z4)
135 self.zeta4 = z4
136 self.log(f" airy modes: {self.n_airy} in {time.perf_counter() - t0:.1f}s")
137 # coupling tensor H[p, i, j] by 3 x 3 Gauss quadrature per cell
138 t0 = time.perf_counter()
139 xq, yq, WQ = pl.quadrature(3)
140 Bxx, Byy, Bxy = pl.basis_derivs_at(xq, yq)
141 Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
142 Bn = pl.basis_at(xq, yq)
143 PsiQ = Bn @ Psi # (Q, n_airy)
144 self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
145 if self.torch is None:
146 H = np.zeros((self.n_airy, self.n_w, self.n_w))
147 for pi in range(self.n_airy):
148 wp = (WQ * PsiQ[:, pi])[:, None]
149 A = (Pxx * wp).T @ Pyy
150 H[pi] = A + A.T - 2.0 * (Pxy * wp).T @ Pxy
151 self.H = H
152 self.Hflat = H.reshape(self.n_airy * self.n_w, self.n_w)
153 shape = H.shape
154 else:
155 # the same assembly on the device; the tensor stays there for the time stepping
156 t = self.torch
157 dev = t.device(self.device)
158 Pxx_t, Pyy_t, Pxy_t = (t.as_tensor(a, device=dev) for a in (Pxx, Pyy, Pxy))
159 WPsi = t.as_tensor(WQ[:, None] * PsiQ, device=dev) # (Q, n_airy)
160 H = t.empty((self.n_airy, self.n_w, self.n_w), dtype=t.float64, device=dev)
161 for pi in range(self.n_airy):
162 wp = WPsi[:, pi : pi + 1]
163 A = (Pxx_t * wp).T @ Pyy_t
164 H[pi] = A + A.T - 2.0 * (Pxy_t * wp).T @ Pxy_t
165 t.cuda.synchronize() if self.device.startswith("cuda") else None
166 self.H = None
167 self.Hflat_t = H.reshape(self.n_airy * self.n_w, self.n_w)
168 shape = tuple(H.shape)
169 del Pxx_t, Pyy_t, Pxy_t, WPsi
170 self.log(f" coupling tensor {shape} on {self.device} in {time.perf_counter() - t0:.1f}s")
171 self.G = None # shell coupling (n_airy, n_w), None for the flat disk
172 self.V = V
173 if self.profile is not None:
174 t0 = time.perf_counter()
175 self._apply_profile(xq, yq, WQ, PsiQ, Pxx, Pyy, Pxy)
176 V = self.V
177 self.log(f" shell modes: f_low={np.round(self.omega[self.omega > 0][:4] / (2 * np.pi), 2)} Hz, f_max={self.omega.max() / (2 * np.pi):.0f} Hz in {time.perf_counter() - t0:.1f}s")
178 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
179 _, _, self.tau, self.cq = pl.listener_quadrature(3)
180 self.PhiQ = Bn @ V # (Q, n_w)
182 def _apply_profile(self, xq, yq, WQ, PsiQ, Pxx, Pyy, Pxy):
183 """Marguerre shallow-shell terms for the initial shape: the G tensor, the shell modes (rotation T of the
184 flat modes), the rotated coupling tensors, the strike on the curved top face and the listener distances."""
185 pl, pr, p = self.plate, self.profile, self.p
186 t = self.torch
187 zxx, zyy, zxy = pr.curvatures(xq, yq)
188 # G[p, j] = int Psi_p L(w0, Phi_j) dA, L(w0, Phi) = Z_xx Phi_yy + Z_yy Phi_xx - 2 Z_xy Phi_xy
189 Lw0 = zxx[:, None] * Pyy + zyy[:, None] * Pxx - 2.0 * zxy[:, None] * Pxy # (Q, n_w)
190 G = (WQ[:, None] * PsiQ).T @ Lw0 # (n_airy, n_w)
191 # linear membrane stiffening in the flat-mode basis: K = diag(omega^2) - 2 G^T diag(c) G (c < 0)
192 K = np.diag(self.omega**2) - 2.0 * (G.T * self.coef) @ G
193 K = 0.5 * (K + K.T)
194 w2, T = np.linalg.eigh(K)
195 omega = np.sqrt(np.maximum(w2, 0.0))
196 omega[omega < 1.0] = 0.0 # rigid modes stay rigid (L(w0, harmonic) = 0)
197 # rotate: H'^p = T^T H^p T, G' = G T, f' = T^T f, V' = V T
198 if t is None:
199 Hr = np.einsum("pij,ik,jl->pkl", self.H, T, T, optimize=True)
200 self.H = Hr
201 self.Hflat = Hr.reshape(self.n_airy * self.n_w, self.n_w)
202 else:
203 dev = t.device(self.device)
204 Tt = t.as_tensor(T, device=dev)
205 Ht = self.Hflat_t.reshape(self.n_airy, self.n_w, self.n_w)
206 Ht = t.einsum("pij,ik,jl->pkl", Ht, Tt, Tt)
207 self.Hflat_t = Ht.reshape(self.n_airy * self.n_w, self.n_w)
208 self.G = G @ T
209 self.V = self.V @ T
210 self.omega = omega
211 self.alpha = p["sigma0"] + 0.5 * p["eta"] * omega**2
212 # strike: -int_{S+} p Phi_k dS with p = F/C phi(r'/R_s) psi, r' the horizontal distance of the top-face point,
213 # dS = dA / n_z; quadrature over the projected footprint in polar coordinates about the strike center
214 from numpy.polynomial import legendre as leg
216 Rs, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
217 C = p["strike_normalization_m2"]
218 gr, gw = leg.leggauss(32)
219 rq = (Rs * 1.1) * (gr + 1.0) / 2.0
220 wr = gw * (Rs * 1.1) / 2.0
221 nth = 96
222 th = 2 * np.pi * np.arange(nth) / nth
223 rr, tt = np.meshgrid(rq, th, indexing="ij")
224 wa = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
225 px = (xs + rr * np.cos(tt)).ravel()
226 py = (ys + rr * np.sin(tt)).ravel()
227 tx, ty, _ = pr.top(px, py)
228 rprime = np.sqrt((tx - xs) ** 2 + (ty - ys) ** 2)
229 phi = np.where(rprime <= Rs, 0.5 * (1.0 + np.cos(np.pi * np.minimum(rprime / Rs, 1.0))), 0.0)
230 _, _, nz = pr.normal(px, py)
231 Phi_s = pl.basis_at(px, py) @ self.V
232 self.f = -(p["strike_peak_force"] / C) * ((wa * phi / nz) @ Phi_s)
233 self.total_force_check = float((p["strike_peak_force"] / C) * np.sum(wa * phi / nz)) # should be ~F_max
234 # listener: n . a dS = w_tt dA in the shallow-shell model; distances to the top-face points
235 xq_, yq_, wq_ = pl.quadrature(3)
236 tx, ty, tz = pr.top(xq_, yq_)
237 r = np.sqrt((p["listener_x"] - tx) ** 2 + (p["listener_y"] - ty) ** 2 + (p["listener_z"] - tz) ** 2)
238 pl._listener_q = (xq_, yq_, r / p["c_air"], p["rho_air"] / (2 * np.pi) * wq_ / r)
239 # the plate object doubles as the closed-form linear solver of the shell (for checks): give it the shell modes
240 pl.omega, pl.alpha, pl.f, pl.V = self.omega, self.alpha, self.f, self.V
241 pl.rigid = self.omega == 0.0
242 pl.n_rigid = int(pl.rigid.sum())
243 pl._setup_linear_response()
245 def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
246 v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
247 sq = v @ q # q^T H^p q
248 if self.G is None:
249 return v.T @ (self.coef * sq) # F_k = sum_p c_p (q^T H^p q) v[p,k]
250 g = self.G @ q # G^p . q
251 return v.T @ (self.coef * (sq + 2.0 * g)) + self.G.T @ (self.coef * sq)
253 def _prepare(self, fs: int, nsamp: int, m: int, field_grid, field_rate):
254 p = self.p
255 dt = 1.0 / (fs * m)
256 w2, al = self.omega**2, self.alpha
257 # exact two-step recurrence of the damped linear oscillator: q+ = a1 q + a2 q- + b F
258 rigid = self.omega == 0.0
259 disc = np.sqrt(al**2 - w2 + 0j)
260 E1, E2 = np.exp((-al + disc) * dt), np.exp((-al - disc) * dt)
261 a1 = (E1 + E2).real
262 a2 = -(E1 * E2).real
263 with np.errstate(divide="ignore", invalid="ignore"):
264 b = np.where(rigid, 0.0, (1.0 - a1 - a2) / np.where(rigid, 1.0, w2))
265 if rigid.any():
266 e2 = np.exp(-2.0 * al[rigid] * dt)
267 a1[rigid] = 1.0 + e2
268 a2[rigid] = -e2
269 b[rigid] = dt * (1.0 - e2) / (2.0 * al[rigid])
270 # fractional delays folded into a per-step weight table G[step, mode]
271 tau_steps = self.tau / dt
272 d0 = np.floor(tau_steps).astype(int)
273 fr = tau_steps - d0
274 lag = np.stack([
275 -fr * (fr - 1.0) * (fr - 2.0) / 6.0,
276 (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0,
277 -(fr + 1.0) * fr * (fr - 2.0) / 2.0,
278 (fr + 1.0) * fr * (fr - 1.0) / 6.0,
279 ], axis=1)
280 dmax = int(d0.max()) + 2
281 G = np.zeros((dmax + 1, self.n_w))
282 for j, off in enumerate([-1, 0, 1, 2]):
283 dd = d0 + off
284 for u in np.unique(dd):
285 sel = dd == u
286 G[u] += (self.cq[sel] * lag[sel, j]) @ self.PhiQ[sel]
287 fld = None
288 if field_grid is not None:
289 fx, fy, ft = field_grid
290 XX, YY = np.meshgrid(fx, fy, indexing="xy")
291 inside = self.plate.inside(XX.ravel(), YY.ravel())
292 PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.V
293 if self.profile is not None:
294 PhiG = PhiG * self.profile.normal(XX.ravel(), YY.ravel())[2][:, None] # u . n = n_z w
295 fld = dict(PhiG=PhiG, inside=inside, every=int(round(fs * m / field_rate)), nt=len(ft), ny=len(fy), nx=len(fx))
296 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)
298 def run(self, fs: int, nsamp: int, m: int, field_grid=None, field_rate: float | None = None):
299 pr = self._prepare(fs, nsamp, m, field_grid, field_rate)
300 if self.torch is not None:
301 return self._run_device(pr, m, nsamp)
302 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"))
303 f = self.f
304 ring = np.zeros((dmax + 1, self.n_w))
305 out = np.zeros(nsamp)
306 q_prev = np.zeros(self.n_w)
307 q = np.zeros(self.n_w)
308 frames = [] if fld else None
309 t_start = time.perf_counter()
310 report = max(1, nsteps // 10)
311 for nstep in range(nsteps + 1):
312 t = nstep * dt
313 psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
314 Fext = self.nonlinear_force(q) + f * psi
315 q_next = a1 * q + a2 * q_prev + b * Fext
316 qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
317 ring[nstep % (dmax + 1)] = qdd
318 if nstep % m == 0:
319 k = nstep // m
320 if k < nsamp:
321 idx = (nstep - np.arange(dmax + 1)) % (dmax + 1)
322 out[k] = np.sum(G * ring[idx])
323 if frames is not None and nstep % fld["every"] == 0 and len(frames) < fld["nt"]:
324 fr_ = (fld["PhiG"] @ q).astype(np.float32)
325 fr_[~fld["inside"]] = np.nan
326 frames.append(fr_)
327 q_prev, q = q, q_next
328 if nstep % report == 0 and nstep:
329 if not np.all(np.isfinite(q)):
330 raise RuntimeError(f"solution diverged at step {nstep} (t = {t:.4f} s); increase the oversampling m")
331 el = time.perf_counter() - t_start
332 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}")
333 if not np.all(np.isfinite(out)):
334 raise RuntimeError("output is not finite; increase the oversampling m")
335 fields = np.stack(frames).reshape(len(frames), fld["ny"], fld["nx"]) if frames is not None else None
336 return out, fields
338 def _run_device(self, pr, m: int, nsamp: int):
339 """The same loop with the tensor contraction and the recurrence on the device (PyTorch)."""
340 t = self.torch
341 dev = t.device(self.device)
342 T64 = lambda a: t.as_tensor(np.asarray(a, dtype=np.float64), device=dev)
343 dt, T, Om, nsteps, dmax, fld = pr["dt"], pr["T"], pr["Om"], pr["nsteps"], pr["dmax"], pr["fld"]
344 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))
345 Hf = self.Hflat_t
346 Gt = T64(self.G) if self.G is not None else None
347 L = dmax + 1
348 ar = t.arange(L, device=dev)
349 ring = t.zeros((L, self.n_w), dtype=t.float64, device=dev)
350 out = t.zeros(nsamp, dtype=t.float64, device=dev)
351 q = t.zeros(self.n_w, dtype=t.float64, device=dev)
352 q_prev = t.zeros_like(q)
353 frames = [] if fld else None
354 PhiG = T64(fld["PhiG"]) if fld else None
355 invdt = 1.0 / dt # 2 alpha (q+ - q-) / (2 dt)
356 t_start = time.perf_counter()
357 report = max(1, nsteps // 10)
358 for nstep in range(nsteps + 1):
359 tt = nstep * dt
360 psi = 0.5 * (1.0 - np.cos(Om * tt)) if 0.0 <= tt <= T else 0.0
361 v = (Hf @ q).reshape(self.n_airy, self.n_w)
362 sq = v @ q
363 if Gt is None:
364 Fext = v.T @ (coef * sq)
365 else:
366 Fext = v.T @ (coef * (sq + 2.0 * (Gt @ q))) + Gt.T @ (coef * sq)
367 if psi:
368 Fext = Fext + f * psi
369 q_next = a1 * q + a2 * q_prev + b * Fext
370 qdd = Fext - w2 * q - al * ((q_next - q_prev) * invdt)
371 ring[nstep % L] = qdd
372 if nstep % m == 0:
373 k = nstep // m
374 if k < nsamp:
375 idx = (nstep - ar) % L
376 out[k] = t.sum(G * ring[idx])
377 if frames is not None and nstep % fld["every"] == 0 and len(frames) < fld["nt"]:
378 fr_ = (PhiG @ q).cpu().numpy().astype(np.float32)
379 fr_[~fld["inside"]] = np.nan
380 frames.append(fr_)
381 q_prev, q = q, q_next
382 if nstep % report == 0 and nstep:
383 qc = q.cpu().numpy()
384 if not np.all(np.isfinite(qc)):
385 raise RuntimeError(f"solution diverged at step {nstep} (t = {tt:.4f} s); increase the oversampling m")
386 el = time.perf_counter() - t_start
387 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}")
388 out = out.cpu().numpy()
389 if not np.all(np.isfinite(out)):
390 raise RuntimeError("output is not finite; increase the oversampling m")
391 fields = np.stack(frames).reshape(len(frames), fld["ny"], fld["nx"]) if frames is not None else None
392 return out, fields
395def hardware_string(device: str) -> str:
396 cpu = "unknown CPU"
397 try:
398 for line in open("/proc/cpuinfo"):
399 if line.startswith("model name"):
400 cpu = " ".join(line.split(":", 1)[1].split())
401 break
402 except OSError:
403 pass
404 s = f"{cpu} ({os.cpu_count()} threads), numpy/scipy for the setup"
405 if device != "cpu":
406 import torch
408 s += f"; coupling tensor and time stepping on {torch.cuda.get_device_name(0)} (PyTorch {torch.__version__})"
409 else:
410 s += " and the time stepping"
411 return s
414def parse_tag(tag: str):
415 parts = dict((kv[0], int(kv[1:])) for kv in tag.split("-"))
416 return parts["g"] * 1e-3, parts["w"], parts["a"], parts["m"]
419def main():
420 ch = load_challenge(Path(sys.argv[1]))
421 out_dir = Path(sys.argv[2])
422 args = sys.argv[3:]
423 fields_for = set()
424 force = None
425 seconds = None
427 def take(flag, conv):
428 nonlocal args
429 if flag in args:
430 i = args.index(flag)
431 v = conv(args[i + 1])
432 args = args[:i] + args[i + 2:]
433 return v
434 return None
436 v = take("--fields", lambda s: set(s.split(",")))
437 fields_for = v or set()
438 force = take("--force", float)
439 duration = take("--duration", float)
440 seconds = take("--seconds", float)
441 device = take("--device", str) or "cpu"
442 full = "--full" in args
443 if full:
444 args.remove("--full")
445 tags = args or ["g8-w150-a150-m1", "g8-w300-a300-m1", "g4-w300-a300-m2", "g4-w400-a400-m2"]
446 out = ch["output"]
447 fs = out["sample_rate_hz"]
448 nsamp = out["num_samples"] if seconds is None else int(round(fs * seconds))
449 log = lambda s: print(s, file=sys.stderr, flush=True)
450 meta = []
451 for tag in tags:
452 spacing, n_w, n_a, m = parse_tag(tag)
453 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")
454 t0 = time.perf_counter()
455 model = VonKarmanDisk(ch, spacing, n_w, n_a, force, log=log, symmetric_y=False if full else None, device=device, duration_override=duration)
456 t1 = time.perf_counter()
457 grid = None
458 if tag in fields_for:
459 fg = ch["fields"]
460 nt, ny, nx = fg["shape"]
461 R = ch["parameters"]["disk_radius"]
462 grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
463 s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
464 dt_total = time.perf_counter() - t0
465 sid = f"vk-{tag}"
466 out_dir.mkdir(parents=True, exist_ok=True)
467 scale = out["signal_scale"] if out.get("signal_scale") else 1.0 # a draft challenge may have no scale yet
468 wavfile.write(out_dir / f"{sid}.wav", fs, (s / scale).astype(np.float32))
469 entry = {
470 "solution_id": sid,
471 "output": f"{out_dir.name}/{sid}.wav",
472 "method": ("Marguerre shallow shell (Airy form) on a Cartesian grid masked to the disk, modal method on shell modes built from Bogner-Fox-Schmit plate modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh-type integral with cubic fractional delays" if model.profile is not None else
473 "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"),
474 "parameters": {
475 "grid_spacing_m": spacing,
476 "cells": model.plate.n_cells,
477 "n_dof": model.plate.ndof,
478 "transverse_modes": model.n_w,
479 "basis_even_in_y": model.symmetric_y,
480 "airy_modes": model.n_airy,
481 "oversampling": m,
482 "dt_s": 1.0 / (fs * m),
483 "highest_transverse_mode_hz": round(float(model.omega.max() / (2 * np.pi)), 1),
484 },
485 "timing_seconds": round(dt_total, 3),
486 "hardware": hardware_string(device),
487 "peak_si": float(np.max(np.abs(s))),
488 "setup_seconds": round(t1 - t0, 3),
489 }
490 if fields is not None:
491 import h5py
493 fg = ch["fields"]
494 with h5py.File(out_dir / f"{sid}.h5", "w") as fh:
495 fh.attrs["challenge_id"] = ch["challenge_id"]
496 fh.attrs["solution_id"] = sid
497 fh.attrs["quantity"] = fg["quantity"]
498 fh.attrs["units"] = fg["units"]
499 fh.create_dataset("x", data=grid[0])
500 fh.create_dataset("y", data=grid[1])
501 fh.create_dataset("t", data=grid[2][: fields.shape[0]])
502 fh.create_dataset("u", data=fields, compression="gzip", compression_opts=4, chunks=(8, fields.shape[1], fields.shape[2]))
503 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
504 meta.append(entry)
505 log(f"{tag}: done, setup={t1 - t0:.1f}s total={dt_total:.1f}s peak={np.max(np.abs(s)):.4e} Pa")
506 with open(out_dir / f"{sid}.meta.json", "w") as fh:
507 json.dump(entry, fh, indent=1)
508 print(json.dumps(meta, indent=1))
511if __name__ == "__main__":
512 main()