Marguerre shallow-shell terms for the bowed disk (cymbal-bow-v1)
4 changed files+144−6
solvers/vonkarman-modal-disk/__pycache__/bfs_disk.cpython-313.pycmodified+0−0View file
Binary file not shown.
solvers/vonkarman-modal-disk/__pycache__/vk_modal_disk.cpython-313.pycmodified+0−0View file
Binary file not shown.
solvers/vonkarman-modal-disk/bfs_disk.pymodified+9−0View file
@@ -170,6 +170,13 @@ class BFSDisk:
170170 phat = p["strike_peak_force"] / (c_phi * np.pi * Rs**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / Rs))
171171 self.f = -(ww * phat) @ self.modes_at(px, py)
172172
173+ self._listener_q = None # optional override of the listener quadrature (used by the shell solver)
174+ self._setup_linear_response()
175+
176+ def _setup_linear_response(self):
177+ """Closed-form constants of the damped modal response to the Hann strike (see the square-plate solver).
178+ Depends on omega, alpha, f and rigid; call again after changing them."""
179+ p = self.p
173180 # ---- closed-form linear response constants (see the square-plate solver)
174181 T = p["strike_duration"]
175182 Om = 2 * np.pi / T
@@ -378,6 +385,8 @@ class BFSDisk:
378385
379386 # ------------------------------------------------------------ Rayleigh integral (linear, exact retardation)
380387 def listener_quadrature(self, order: int = 3):
388+ if self._listener_q is not None:
389+ return self._listener_q
381390 p = self.p
382391 xq, yq, wq = self.quadrature(order)
383392 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
solvers/vonkarman-modal-disk/vk_modal_disk.pymodified+135−6View file
@@ -27,6 +27,17 @@ Because the strike is centered on the x axis, the solution is even in y, and by
2727 modes needed for a given frequency reach (the coupling cost scales with the cube of the count). --full
2828 uses the unrestricted bases instead.
2929
30+If the challenge has a shell-of-revolution geometry (cymbal-bow-v1: a paraboloidal bow z = Z(r)), the same
31+solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
32+system: L(w + w0, Phi) on the right of the bending equation and 2 L(w0, w) in the compatibility equation.
33+With 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,
34+which is diagonalized once in the flat-mode basis to give the shell modes (exact in the Gautschi step,
35+damped at sigma0 + eta omega^2 / 2 like every mode, as the Kelvin-Voigt law of the challenge implies), and
36+quadratic 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
37+acts along -z on the curved top face with the normalization of the specification, the listener integral uses
38+the distance to the top face, and the field export is the normal displacement n_z w. Shallow-shell theory
39+neglects in-plane inertia and terms of the order of the slope squared (0.06 at the rim here).
40+
3041 With --device cuda the coupling tensor is assembled and the time stepping runs on a GPU through PyTorch
3142 (the modes and quadrature stay in numpy/scipy); the CPU path is the reference and the two agree to round-off.
3243
@@ -47,6 +58,43 @@ sys.path.insert(0, str(Path(__file__).resolve().parent))
4758 from bfs_disk import BFSDisk, load_challenge # noqa: E402
4859
4960
61+class Paraboloid:
62+ """The bow of cymbal-bow-v1: mid-surface z = Z(r) = H (1 - r^2 / R^2), uniform thickness along the normal."""
63+
64+ def __init__(self, R: float, H: float, h: float):
65+ self.R, self.H, self.h = R, H, h
66+
67+ @staticmethod
68+ def from_challenge(ch: dict):
69+ g = ch.get("geometry")
70+ if not g or g.get("type") != "shell-of-revolution":
71+ return None
72+ p = ch["parameters"]
73+ return Paraboloid(p["disk_radius"], p["bow_height"], p["h"])
74+
75+ def z(self, x, y):
76+ return self.H * (1.0 - (x * x + y * y) / self.R**2)
77+
78+ def grad(self, x, y):
79+ c = -2.0 * self.H / self.R**2
80+ return c * x, c * y
81+
82+ def curvatures(self, x, y):
83+ """Second derivatives (Z_xx, Z_yy, Z_xy) of the mid-surface height."""
84+ c = -2.0 * self.H / self.R**2
85+ return c * np.ones_like(x), c * np.ones_like(x), np.zeros_like(x)
86+
87+ def normal(self, x, y):
88+ zx, zy = self.grad(x, y)
89+ nz = 1.0 / np.sqrt(1.0 + zx * zx + zy * zy)
90+ return -zx * nz, -zy * nz, nz
91+
92+ def top(self, x, y):
93+ """The point of the top face on the normal through the mid-surface point (x, y)."""
94+ nx, ny, nz = self.normal(x, y)
95+ return x + 0.5 * self.h * nx, y + 0.5 * self.h * ny, self.z(x, y) + 0.5 * self.h * nz
96+
97+
5098 class VonKarmanDisk:
5199 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):
52100 self.ch = ch
@@ -68,6 +116,7 @@ class VonKarmanDisk:
68116 if symmetric_y is None:
69117 symmetric_y = p["strike_y"] == 0.0
70118 self.symmetric_y = symmetric_y
119+ self.profile = Paraboloid.from_challenge(ch)
71120 t0 = time.perf_counter()
72121 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
73122 pl = self.plate
@@ -119,14 +168,87 @@ class VonKarmanDisk:
119168 shape = tuple(H.shape)
120169 del Pxx_t, Pyy_t, Pxy_t, WPsi
121170 self.log(f" coupling tensor {shape} on {self.device} in {time.perf_counter() - t0:.1f}s")
171+ self.G = None # shell coupling (n_airy, n_w), None for the flat disk
172+ self.V = V
173+ if self.profile is not None:
174+ t0 = time.perf_counter()
175+ self._apply_profile(xq, yq, WQ, PsiQ, Pxx, Pyy, Pxy)
176+ V = self.V
177+ self.log(f" shell modes: f_low={np.round(self.omega[self.omega > 0][:4] / (2 * np.pi), 2)} Hz, f_max={self.omega.max() / (2 * np.pi):.0f} Hz in {time.perf_counter() - t0:.1f}s")
122178 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
123179 _, _, self.tau, self.cq = pl.listener_quadrature(3)
124180 self.PhiQ = Bn @ V # (Q, n_w)
125181
182+ def _apply_profile(self, xq, yq, WQ, PsiQ, Pxx, Pyy, Pxy):
183+ """Marguerre shallow-shell terms for the initial shape: the G tensor, the shell modes (rotation T of the
184+ flat modes), the rotated coupling tensors, the strike on the curved top face and the listener distances."""
185+ pl, pr, p = self.plate, self.profile, self.p
186+ t = self.torch
187+ zxx, zyy, zxy = pr.curvatures(xq, yq)
188+ # G[p, j] = int Psi_p L(w0, Phi_j) dA, L(w0, Phi) = Z_xx Phi_yy + Z_yy Phi_xx - 2 Z_xy Phi_xy
189+ Lw0 = zxx[:, None] * Pyy + zyy[:, None] * Pxx - 2.0 * zxy[:, None] * Pxy # (Q, n_w)
190+ G = (WQ[:, None] * PsiQ).T @ Lw0 # (n_airy, n_w)
191+ # linear membrane stiffening in the flat-mode basis: K = diag(omega^2) - 2 G^T diag(c) G (c < 0)
192+ K = np.diag(self.omega**2) - 2.0 * (G.T * self.coef) @ G
193+ K = 0.5 * (K + K.T)
194+ w2, T = np.linalg.eigh(K)
195+ omega = np.sqrt(np.maximum(w2, 0.0))
196+ omega[omega < 1.0] = 0.0 # rigid modes stay rigid (L(w0, harmonic) = 0)
197+ # rotate: H'^p = T^T H^p T, G' = G T, f' = T^T f, V' = V T
198+ if t is None:
199+ Hr = np.einsum("pij,ik,jl->pkl", self.H, T, T, optimize=True)
200+ self.H = Hr
201+ self.Hflat = Hr.reshape(self.n_airy * self.n_w, self.n_w)
202+ else:
203+ dev = t.device(self.device)
204+ Tt = t.as_tensor(T, device=dev)
205+ Ht = self.Hflat_t.reshape(self.n_airy, self.n_w, self.n_w)
206+ Ht = t.einsum("pij,ik,jl->pkl", Ht, Tt, Tt)
207+ self.Hflat_t = Ht.reshape(self.n_airy * self.n_w, self.n_w)
208+ self.G = G @ T
209+ self.V = self.V @ T
210+ self.omega = omega
211+ self.alpha = p["sigma0"] + 0.5 * p["eta"] * omega**2
212+ # strike: -int_{S+} p Phi_k dS with p = F/C phi(r'/R_s) psi, r' the horizontal distance of the top-face point,
213+ # dS = dA / n_z; quadrature over the projected footprint in polar coordinates about the strike center
214+ from numpy.polynomial import legendre as leg
215+
216+ Rs, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
217+ C = p["strike_normalization_m2"]
218+ gr, gw = leg.leggauss(32)
219+ rq = (Rs * 1.1) * (gr + 1.0) / 2.0
220+ wr = gw * (Rs * 1.1) / 2.0
221+ nth = 96
222+ th = 2 * np.pi * np.arange(nth) / nth
223+ rr, tt = np.meshgrid(rq, th, indexing="ij")
224+ wa = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
225+ px = (xs + rr * np.cos(tt)).ravel()
226+ py = (ys + rr * np.sin(tt)).ravel()
227+ tx, ty, _ = pr.top(px, py)
228+ rprime = np.sqrt((tx - xs) ** 2 + (ty - ys) ** 2)
229+ phi = np.where(rprime <= Rs, 0.5 * (1.0 + np.cos(np.pi * np.minimum(rprime / Rs, 1.0))), 0.0)
230+ _, _, nz = pr.normal(px, py)
231+ Phi_s = pl.basis_at(px, py) @ self.V
232+ self.f = -(p["strike_peak_force"] / C) * ((wa * phi / nz) @ Phi_s)
233+ self.total_force_check = float((p["strike_peak_force"] / C) * np.sum(wa * phi / nz)) # should be ~F_max
234+ # listener: n . a dS = w_tt dA in the shallow-shell model; distances to the top-face points
235+ xq_, yq_, wq_ = pl.quadrature(3)
236+ tx, ty, tz = pr.top(xq_, yq_)
237+ r = np.sqrt((p["listener_x"] - tx) ** 2 + (p["listener_y"] - ty) ** 2 + (p["listener_z"] - tz) ** 2)
238+ pl._listener_q = (xq_, yq_, r / p["c_air"], p["rho_air"] / (2 * np.pi) * wq_ / r)
239+ # the plate object doubles as the closed-form linear solver of the shell (for checks): give it the shell modes
240+ pl.omega, pl.alpha, pl.f, pl.V = self.omega, self.alpha, self.f, self.V
241+ pl.rigid = self.omega == 0.0
242+ pl.n_rigid = int(pl.rigid.sum())
243+ pl._setup_linear_response()
244+
126245 def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
127246 v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
128- eta = self.coef * (v @ q) # eta_p
129- return v.T @ eta # F_k = sum_p v[p,k] eta_p
247+ sq = v @ q # q^T H^p q
248+ if self.G is None:
249+ return v.T @ (self.coef * sq) # F_k = sum_p c_p (q^T H^p q) v[p,k]
250+ g = self.G @ q # G^p . q
251+ return v.T @ (self.coef * (sq + 2.0 * g)) + self.G.T @ (self.coef * sq)
130252
131253 def _prepare(self, fs: int, nsamp: int, m: int, field_grid, field_rate):
132254 p = self.p
@@ -167,7 +289,9 @@ class VonKarmanDisk:
167289 fx, fy, ft = field_grid
168290 XX, YY = np.meshgrid(fx, fy, indexing="xy")
169291 inside = self.plate.inside(XX.ravel(), YY.ravel())
170- PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.plate.V
292+ PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.V
293+ if self.profile is not None:
294+ PhiG = PhiG * self.profile.normal(XX.ravel(), YY.ravel())[2][:, None] # u . n = n_z w
171295 fld = dict(PhiG=PhiG, inside=inside, every=int(round(fs * m / field_rate)), nt=len(ft), ny=len(fy), nx=len(fx))
172296 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)
173297
@@ -219,6 +343,7 @@ class VonKarmanDisk:
219343 dt, T, Om, nsteps, dmax, fld = pr["dt"], pr["T"], pr["Om"], pr["nsteps"], pr["dmax"], pr["fld"]
220344 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))
221345 Hf = self.Hflat_t
346+ Gt = T64(self.G) if self.G is not None else None
222347 L = dmax + 1
223348 ar = t.arange(L, device=dev)
224349 ring = t.zeros((L, self.n_w), dtype=t.float64, device=dev)
@@ -234,8 +359,11 @@ class VonKarmanDisk:
234359 tt = nstep * dt
235360 psi = 0.5 * (1.0 - np.cos(Om * tt)) if 0.0 <= tt <= T else 0.0
236361 v = (Hf @ q).reshape(self.n_airy, self.n_w)
237- eta = coef * (v @ q)
238- Fext = v.T @ eta
362+ sq = v @ q
363+ if Gt is None:
364+ Fext = v.T @ (coef * sq)
365+ else:
366+ Fext = v.T @ (coef * (sq + 2.0 * (Gt @ q))) + Gt.T @ (coef * sq)
239367 if psi:
240368 Fext = Fext + f * psi
241369 q_next = a1 * q + a2 * q_prev + b * Fext
@@ -340,7 +468,8 @@ def main():
340468 entry = {
341469 "solution_id": sid,
342470 "output": f"{out_dir.name}/{sid}.wav",
343- "method": "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",
471+ "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
472+ "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"),
344473 "parameters": {
345474 "grid_spacing_m": spacing,
346475 "cells": model.plate.n_cells,