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