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