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.
30Usage: python vk_modal_disk.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [--seconds S] [--full] [tag ...]
31 where a tag is g<grid spacing in mm>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. g8-w200-a200-m1
32"""
33from __future__ import annotations
35import json
36import sys
37import time
38from pathlib import Path
40import numpy as np
41from scipy.io import wavfile
42sys.path.insert(0, str(Path(__file__).resolve().parent))
43from bfs_disk import BFSDisk, load_challenge # noqa: E402
46class VonKarmanDisk:
47 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):
48 self.ch = ch
49 p = dict(ch["parameters"])
50 if force_override is not None:
51 p["strike_peak_force"] = force_override
52 ch2 = dict(ch)
53 ch2["parameters"] = p
54 self.p = p
55 self.log = log or (lambda s: None)
56 if symmetric_y is None:
57 symmetric_y = p["strike_y"] == 0.0
58 self.symmetric_y = symmetric_y
59 t0 = time.perf_counter()
60 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
61 pl = self.plate
62 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}")
63 self.n_w = pl.n_modes
64 self.omega = pl.omega
65 self.alpha = pl.alpha
66 self.f = pl.f
67 V = pl.V
68 # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
69 t0 = time.perf_counter()
70 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]])))
71 M1 = pl.M / pl.m
72 z4, Psi = pl.eigenmodes(Kb, M1, n_airy, sigma=0.0, dofs=pl.interior_dofs)
73 self.n_airy = len(z4)
74 self.zeta4 = z4
75 self.log(f" airy modes: {self.n_airy} in {time.perf_counter() - t0:.1f}s")
76 # coupling tensor H[p, i, j] by 3 x 3 Gauss quadrature per cell
77 t0 = time.perf_counter()
78 xq, yq, WQ = pl.quadrature(3)
79 Bxx, Byy, Bxy = pl.basis_derivs_at(xq, yq)
80 Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
81 Bn = pl.basis_at(xq, yq)
82 PsiQ = Bn @ Psi # (Q, n_airy)
83 H = np.zeros((self.n_airy, self.n_w, self.n_w))
84 for pi in range(self.n_airy):
85 wp = (WQ * PsiQ[:, pi])[:, None]
86 A = (Pxx * wp).T @ Pyy
87 H[pi] = A + A.T - 2.0 * (Pxy * wp).T @ Pxy
88 self.H = H
89 self.Hflat = H.reshape(self.n_airy * self.n_w, self.n_w)
90 self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
91 self.log(f" coupling tensor {H.shape} in {time.perf_counter() - t0:.1f}s")
92 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
93 _, _, self.tau, self.cq = pl.listener_quadrature(3)
94 self.PhiQ = Bn @ V # (Q, n_w)
96 def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
97 v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
98 eta = self.coef * (v @ q) # eta_p
99 return v.T @ eta # F_k = sum_p v[p,k] eta_p
101 def run(self, fs: int, nsamp: int, m: int, field_grid=None, field_rate: float | None = None):
102 p = self.p
103 dt = 1.0 / (fs * m)
104 T, Om = p["strike_duration"], 2 * np.pi / p["strike_duration"]
105 nsteps = nsamp * m
106 w2, al, f = self.omega**2, self.alpha, self.f
107 # exact two-step recurrence of the damped linear oscillator: q+ = a1 q + a2 q- + b F
108 rigid = self.omega == 0.0
109 disc = np.sqrt(al**2 - w2 + 0j)
110 E1, E2 = np.exp((-al + disc) * dt), np.exp((-al - disc) * dt)
111 a1 = (E1 + E2).real
112 a2 = -(E1 * E2).real
113 with np.errstate(divide="ignore", invalid="ignore"):
114 b = np.where(rigid, 0.0, (1.0 - a1 - a2) / np.where(rigid, 1.0, w2))
115 if rigid.any():
116 e2 = np.exp(-2.0 * al[rigid] * dt)
117 a1[rigid] = 1.0 + e2
118 a2[rigid] = -e2
119 b[rigid] = dt * (1.0 - e2) / (2.0 * al[rigid])
120 # fractional delays folded into a per-step weight table G[step, mode]
121 tau_steps = self.tau / dt
122 d0 = np.floor(tau_steps).astype(int)
123 fr = tau_steps - d0
124 lag = np.stack([
125 -fr * (fr - 1.0) * (fr - 2.0) / 6.0,
126 (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0,
127 -(fr + 1.0) * fr * (fr - 2.0) / 2.0,
128 (fr + 1.0) * fr * (fr - 1.0) / 6.0,
129 ], axis=1)
130 dmax = int(d0.max()) + 2
131 G = np.zeros((dmax + 1, self.n_w))
132 for j, off in enumerate([-1, 0, 1, 2]):
133 dd = d0 + off
134 for u in np.unique(dd):
135 sel = dd == u
136 G[u] += (self.cq[sel] * lag[sel, j]) @ self.PhiQ[sel]
137 ring = np.zeros((dmax + 1, self.n_w))
138 out = np.zeros(nsamp)
139 q_prev = np.zeros(self.n_w)
140 q = np.zeros(self.n_w)
141 frames = None
142 if field_grid is not None:
143 fx, fy, ft = field_grid
144 XX, YY = np.meshgrid(fx, fy, indexing="xy")
145 inside = self.plate.inside(XX.ravel(), YY.ravel())
146 PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.plate.V
147 frame_every = int(round(fs * m / field_rate))
148 frames = []
149 t_start = time.perf_counter()
150 report = max(1, nsteps // 10)
151 for nstep in range(nsteps + 1):
152 t = nstep * dt
153 psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
154 Fext = self.nonlinear_force(q) + f * psi
155 q_next = a1 * q + a2 * q_prev + b * Fext
156 qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
157 ring[nstep % (dmax + 1)] = qdd
158 if nstep % m == 0:
159 k = nstep // m
160 if k < nsamp:
161 idx = (nstep - np.arange(dmax + 1)) % (dmax + 1)
162 out[k] = np.sum(G * ring[idx])
163 if frames is not None and nstep % frame_every == 0 and len(frames) < len(ft):
164 fr_ = (PhiG @ q).astype(np.float32)
165 fr_[~inside] = np.nan
166 frames.append(fr_)
167 q_prev, q = q, q_next
168 if nstep % report == 0 and nstep:
169 if not np.all(np.isfinite(q)):
170 raise RuntimeError(f"solution diverged at step {nstep} (t = {t:.4f} s); increase the oversampling m")
171 el = time.perf_counter() - t_start
172 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}")
173 if not np.all(np.isfinite(out)):
174 raise RuntimeError("output is not finite; increase the oversampling m")
175 fields = None
176 if frames is not None:
177 fields = np.stack(frames).reshape(len(frames), len(field_grid[1]), len(field_grid[0]))
178 return out, fields
181def parse_tag(tag: str):
182 parts = dict((kv[0], int(kv[1:])) for kv in tag.split("-"))
183 return parts["g"] * 1e-3, parts["w"], parts["a"], parts["m"]
186def main():
187 ch = load_challenge(Path(sys.argv[1]))
188 out_dir = Path(sys.argv[2])
189 args = sys.argv[3:]
190 fields_for = set()
191 force = None
192 seconds = None
194 def take(flag, conv):
195 nonlocal args
196 if flag in args:
197 i = args.index(flag)
198 v = conv(args[i + 1])
199 args = args[:i] + args[i + 2:]
200 return v
201 return None
203 v = take("--fields", lambda s: set(s.split(",")))
204 fields_for = v or set()
205 force = take("--force", float)
206 seconds = take("--seconds", float)
207 full = "--full" in args
208 if full:
209 args.remove("--full")
210 tags = args or ["g8-w150-a150-m1", "g8-w300-a300-m1", "g4-w300-a300-m2", "g4-w400-a400-m2"]
211 out = ch["output"]
212 fs = out["sample_rate_hz"]
213 nsamp = out["num_samples"] if seconds is None else int(round(fs * seconds))
214 log = lambda s: print(s, file=sys.stderr, flush=True)
215 meta = []
216 for tag in tags:
217 spacing, n_w, n_a, m = parse_tag(tag)
218 log(f"{tag}: force={force if force is not None else ch['parameters']['strike_peak_force']} N")
219 t0 = time.perf_counter()
220 model = VonKarmanDisk(ch, spacing, n_w, n_a, force, log=log, symmetric_y=False if full else None)
221 t1 = time.perf_counter()
222 grid = None
223 if tag in fields_for:
224 fg = ch["fields"]
225 nt, ny, nx = fg["shape"]
226 R = ch["parameters"]["disk_radius"]
227 grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
228 s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
229 dt_total = time.perf_counter() - t0
230 sid = f"vk-{tag}"
231 out_dir.mkdir(parents=True, exist_ok=True)
232 wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
233 entry = {
234 "solution_id": sid,
235 "output": f"{out_dir.name}/{sid}.wav",
236 "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",
237 "parameters": {
238 "grid_spacing_m": spacing,
239 "cells": model.plate.n_cells,
240 "n_dof": model.plate.ndof,
241 "transverse_modes": model.n_w,
242 "basis_even_in_y": model.symmetric_y,
243 "airy_modes": model.n_airy,
244 "oversampling": m,
245 "dt_s": 1.0 / (fs * m),
246 "highest_transverse_mode_hz": round(float(model.omega.max() / (2 * np.pi)), 1),
247 },
248 "timing_seconds": round(dt_total, 3),
249 "peak_si": float(np.max(np.abs(s))),
250 "setup_seconds": round(t1 - t0, 3),
251 }
252 if fields is not None:
253 import h5py
255 fg = ch["fields"]
256 with h5py.File(out_dir / f"{sid}.h5", "w") as fh:
257 fh.attrs["challenge_id"] = ch["challenge_id"]
258 fh.attrs["solution_id"] = sid
259 fh.attrs["quantity"] = fg["quantity"]
260 fh.attrs["units"] = fg["units"]
261 fh.create_dataset("x", data=grid[0])
262 fh.create_dataset("y", data=grid[1])
263 fh.create_dataset("t", data=grid[2][: fields.shape[0]])
264 fh.create_dataset("u", data=fields, compression="gzip", compression_opts=4, chunks=(8, fields.shape[1], fields.shape[2]))
265 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
266 meta.append(entry)
267 log(f"{tag}: done, setup={t1 - t0:.1f}s total={dt_total:.1f}s peak={np.max(np.abs(s)):.4e} Pa")
268 with open(out_dir / f"{sid}.meta.json", "w") as fh:
269 json.dump(entry, fh, indent=1)
270 print(json.dumps(meta, indent=1))
273if __name__ == "__main__":
274 main()