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