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