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.
30With --device cuda the coupling tensor is assembled and the time stepping runs on a GPU through PyTorch
31(the modes and quadrature stay in numpy/scipy); the CPU path is the reference and the two agree to round-off.
33Usage: python vk_modal_disk.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [--seconds S] [--full] [--device cuda] [tag ...]
34 where a tag is g<grid spacing in mm>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. g8-w200-a200-m1
35"""
36from __future__ import annotations
38import json
39import os
40import sys
41import time
42from pathlib import Path
44import numpy as np
45from scipy.io import wavfile
46sys.path.insert(0, str(Path(__file__).resolve().parent))
47from bfs_disk import BFSDisk, load_challenge # noqa: E402
50class VonKarmanDisk:
51 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"):
52 self.ch = ch
53 self.device = device
54 self.torch = None
55 if device != "cpu":
56 import torch
58 self.torch = torch
59 p = dict(ch["parameters"])
60 if force_override is not None:
61 p["strike_peak_force"] = force_override
62 ch2 = dict(ch)
63 ch2["parameters"] = p
64 self.p = p
65 self.log = log or (lambda s: None)
66 if symmetric_y is None:
67 symmetric_y = p["strike_y"] == 0.0
68 self.symmetric_y = symmetric_y
69 t0 = time.perf_counter()
70 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
71 pl = self.plate
72 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}")
73 self.n_w = pl.n_modes
74 self.omega = pl.omega
75 self.alpha = pl.alpha
76 self.f = pl.f
77 V = pl.V
78 # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
79 t0 = time.perf_counter()
80 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]])))
81 M1 = pl.M / pl.m
82 z4, Psi = pl.eigenmodes(Kb, M1, n_airy, sigma=0.0, dofs=pl.interior_dofs)
83 self.n_airy = len(z4)
84 self.zeta4 = z4
85 self.log(f" airy modes: {self.n_airy} in {time.perf_counter() - t0:.1f}s")
86 # coupling tensor H[p, i, j] by 3 x 3 Gauss quadrature per cell
87 t0 = time.perf_counter()
88 xq, yq, WQ = pl.quadrature(3)
89 Bxx, Byy, Bxy = pl.basis_derivs_at(xq, yq)
90 Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
91 Bn = pl.basis_at(xq, yq)
92 PsiQ = Bn @ Psi # (Q, n_airy)
93 self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
94 if self.torch is None:
95 H = np.zeros((self.n_airy, self.n_w, self.n_w))
96 for pi in range(self.n_airy):
97 wp = (WQ * PsiQ[:, pi])[:, None]
98 A = (Pxx * wp).T @ Pyy
99 H[pi] = A + A.T - 2.0 * (Pxy * wp).T @ Pxy
100 self.H = H
101 self.Hflat = H.reshape(self.n_airy * self.n_w, self.n_w)
102 shape = H.shape
103 else:
104 # the same assembly on the device; the tensor stays there for the time stepping
105 t = self.torch
106 dev = t.device(self.device)
107 Pxx_t, Pyy_t, Pxy_t = (t.as_tensor(a, device=dev) for a in (Pxx, Pyy, Pxy))
108 WPsi = t.as_tensor(WQ[:, None] * PsiQ, device=dev) # (Q, n_airy)
109 H = t.empty((self.n_airy, self.n_w, self.n_w), dtype=t.float64, device=dev)
110 for pi in range(self.n_airy):
111 wp = WPsi[:, pi : pi + 1]
112 A = (Pxx_t * wp).T @ Pyy_t
113 H[pi] = A + A.T - 2.0 * (Pxy_t * wp).T @ Pxy_t
114 t.cuda.synchronize() if self.device.startswith("cuda") else None
115 self.H = None
116 self.Hflat_t = H.reshape(self.n_airy * self.n_w, self.n_w)
117 shape = tuple(H.shape)
118 del Pxx_t, Pyy_t, Pxy_t, WPsi
119 self.log(f" coupling tensor {shape} on {self.device} in {time.perf_counter() - t0:.1f}s")
120 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
121 _, _, self.tau, self.cq = pl.listener_quadrature(3)
122 self.PhiQ = Bn @ V # (Q, n_w)
124 def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
125 v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
126 eta = self.coef * (v @ q) # eta_p
127 return v.T @ eta # F_k = sum_p v[p,k] eta_p
129 def _prepare(self, fs: int, nsamp: int, m: int, field_grid, field_rate):
130 p = self.p
131 dt = 1.0 / (fs * m)
132 w2, al = self.omega**2, self.alpha
133 # exact two-step recurrence of the damped linear oscillator: q+ = a1 q + a2 q- + b F
134 rigid = self.omega == 0.0
135 disc = np.sqrt(al**2 - w2 + 0j)
136 E1, E2 = np.exp((-al + disc) * dt), np.exp((-al - disc) * dt)
137 a1 = (E1 + E2).real
138 a2 = -(E1 * E2).real
139 with np.errstate(divide="ignore", invalid="ignore"):
140 b = np.where(rigid, 0.0, (1.0 - a1 - a2) / np.where(rigid, 1.0, w2))
141 if rigid.any():
142 e2 = np.exp(-2.0 * al[rigid] * dt)
143 a1[rigid] = 1.0 + e2
144 a2[rigid] = -e2
145 b[rigid] = dt * (1.0 - e2) / (2.0 * al[rigid])
146 # fractional delays folded into a per-step weight table G[step, mode]
147 tau_steps = self.tau / dt
148 d0 = np.floor(tau_steps).astype(int)
149 fr = tau_steps - d0
150 lag = np.stack([
151 -fr * (fr - 1.0) * (fr - 2.0) / 6.0,
152 (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0,
153 -(fr + 1.0) * fr * (fr - 2.0) / 2.0,
154 (fr + 1.0) * fr * (fr - 1.0) / 6.0,
155 ], axis=1)
156 dmax = int(d0.max()) + 2
157 G = np.zeros((dmax + 1, self.n_w))
158 for j, off in enumerate([-1, 0, 1, 2]):
159 dd = d0 + off
160 for u in np.unique(dd):
161 sel = dd == u
162 G[u] += (self.cq[sel] * lag[sel, j]) @ self.PhiQ[sel]
163 fld = None
164 if field_grid is not None:
165 fx, fy, ft = field_grid
166 XX, YY = np.meshgrid(fx, fy, indexing="xy")
167 inside = self.plate.inside(XX.ravel(), YY.ravel())
168 PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.plate.V
169 fld = dict(PhiG=PhiG, inside=inside, every=int(round(fs * m / field_rate)), nt=len(ft), ny=len(fy), nx=len(fx))
170 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)
172 def run(self, fs: int, nsamp: int, m: int, field_grid=None, field_rate: float | None = None):
173 pr = self._prepare(fs, nsamp, m, field_grid, field_rate)
174 if self.torch is not None:
175 return self._run_device(pr, m, nsamp)
176 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"))
177 f = self.f
178 ring = np.zeros((dmax + 1, self.n_w))
179 out = np.zeros(nsamp)
180 q_prev = np.zeros(self.n_w)
181 q = np.zeros(self.n_w)
182 frames = [] if fld else None
183 t_start = time.perf_counter()
184 report = max(1, nsteps // 10)
185 for nstep in range(nsteps + 1):
186 t = nstep * dt
187 psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
188 Fext = self.nonlinear_force(q) + f * psi
189 q_next = a1 * q + a2 * q_prev + b * Fext
190 qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
191 ring[nstep % (dmax + 1)] = qdd
192 if nstep % m == 0:
193 k = nstep // m
194 if k < nsamp:
195 idx = (nstep - np.arange(dmax + 1)) % (dmax + 1)
196 out[k] = np.sum(G * ring[idx])
197 if frames is not None and nstep % fld["every"] == 0 and len(frames) < fld["nt"]:
198 fr_ = (fld["PhiG"] @ q).astype(np.float32)
199 fr_[~fld["inside"]] = np.nan
200 frames.append(fr_)
201 q_prev, q = q, q_next
202 if nstep % report == 0 and nstep:
203 if not np.all(np.isfinite(q)):
204 raise RuntimeError(f"solution diverged at step {nstep} (t = {t:.4f} s); increase the oversampling m")
205 el = time.perf_counter() - t_start
206 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}")
207 if not np.all(np.isfinite(out)):
208 raise RuntimeError("output is not finite; increase the oversampling m")
209 fields = np.stack(frames).reshape(len(frames), fld["ny"], fld["nx"]) if frames is not None else None
210 return out, fields
212 def _run_device(self, pr, m: int, nsamp: int):
213 """The same loop with the tensor contraction and the recurrence on the device (PyTorch)."""
214 t = self.torch
215 dev = t.device(self.device)
216 T64 = lambda a: t.as_tensor(np.asarray(a, dtype=np.float64), device=dev)
217 dt, T, Om, nsteps, dmax, fld = pr["dt"], pr["T"], pr["Om"], pr["nsteps"], pr["dmax"], pr["fld"]
218 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))
219 Hf = self.Hflat_t
220 L = dmax + 1
221 ar = t.arange(L, device=dev)
222 ring = t.zeros((L, self.n_w), dtype=t.float64, device=dev)
223 out = t.zeros(nsamp, dtype=t.float64, device=dev)
224 q = t.zeros(self.n_w, dtype=t.float64, device=dev)
225 q_prev = t.zeros_like(q)
226 frames = [] if fld else None
227 PhiG = T64(fld["PhiG"]) if fld else None
228 invdt = 1.0 / dt # 2 alpha (q+ - q-) / (2 dt)
229 t_start = time.perf_counter()
230 report = max(1, nsteps // 10)
231 for nstep in range(nsteps + 1):
232 tt = nstep * dt
233 psi = 0.5 * (1.0 - np.cos(Om * tt)) if 0.0 <= tt <= T else 0.0
234 v = (Hf @ q).reshape(self.n_airy, self.n_w)
235 eta = coef * (v @ q)
236 Fext = v.T @ eta
237 if psi:
238 Fext = Fext + f * psi
239 q_next = a1 * q + a2 * q_prev + b * Fext
240 qdd = Fext - w2 * q - al * ((q_next - q_prev) * invdt)
241 ring[nstep % L] = qdd
242 if nstep % m == 0:
243 k = nstep // m
244 if k < nsamp:
245 idx = (nstep - ar) % L
246 out[k] = t.sum(G * ring[idx])
247 if frames is not None and nstep % fld["every"] == 0 and len(frames) < fld["nt"]:
248 fr_ = (PhiG @ q).cpu().numpy().astype(np.float32)
249 fr_[~fld["inside"]] = np.nan
250 frames.append(fr_)
251 q_prev, q = q, q_next
252 if nstep % report == 0 and nstep:
253 qc = q.cpu().numpy()
254 if not np.all(np.isfinite(qc)):
255 raise RuntimeError(f"solution diverged at step {nstep} (t = {tt:.4f} s); increase the oversampling m")
256 el = time.perf_counter() - t_start
257 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}")
258 out = out.cpu().numpy()
259 if not np.all(np.isfinite(out)):
260 raise RuntimeError("output is not finite; increase the oversampling m")
261 fields = np.stack(frames).reshape(len(frames), fld["ny"], fld["nx"]) if frames is not None else None
262 return out, fields
265def hardware_string(device: str) -> str:
266 cpu = "unknown CPU"
267 try:
268 for line in open("/proc/cpuinfo"):
269 if line.startswith("model name"):
270 cpu = " ".join(line.split(":", 1)[1].split())
271 break
272 except OSError:
273 pass
274 s = f"{cpu} ({os.cpu_count()} threads), numpy/scipy for the setup"
275 if device != "cpu":
276 import torch
278 s += f"; coupling tensor and time stepping on {torch.cuda.get_device_name(0)} (PyTorch {torch.__version__})"
279 else:
280 s += " and the time stepping"
281 return s
284def parse_tag(tag: str):
285 parts = dict((kv[0], int(kv[1:])) for kv in tag.split("-"))
286 return parts["g"] * 1e-3, parts["w"], parts["a"], parts["m"]
289def main():
290 ch = load_challenge(Path(sys.argv[1]))
291 out_dir = Path(sys.argv[2])
292 args = sys.argv[3:]
293 fields_for = set()
294 force = None
295 seconds = None
297 def take(flag, conv):
298 nonlocal args
299 if flag in args:
300 i = args.index(flag)
301 v = conv(args[i + 1])
302 args = args[:i] + args[i + 2:]
303 return v
304 return None
306 v = take("--fields", lambda s: set(s.split(",")))
307 fields_for = v or set()
308 force = take("--force", float)
309 seconds = take("--seconds", float)
310 device = take("--device", str) or "cpu"
311 full = "--full" in args
312 if full:
313 args.remove("--full")
314 tags = args or ["g8-w150-a150-m1", "g8-w300-a300-m1", "g4-w300-a300-m2", "g4-w400-a400-m2"]
315 out = ch["output"]
316 fs = out["sample_rate_hz"]
317 nsamp = out["num_samples"] if seconds is None else int(round(fs * seconds))
318 log = lambda s: print(s, file=sys.stderr, flush=True)
319 meta = []
320 for tag in tags:
321 spacing, n_w, n_a, m = parse_tag(tag)
322 log(f"{tag}: force={force if force is not None else ch['parameters']['strike_peak_force']} N")
323 t0 = time.perf_counter()
324 model = VonKarmanDisk(ch, spacing, n_w, n_a, force, log=log, symmetric_y=False if full else None, device=device)
325 t1 = time.perf_counter()
326 grid = None
327 if tag in fields_for:
328 fg = ch["fields"]
329 nt, ny, nx = fg["shape"]
330 R = ch["parameters"]["disk_radius"]
331 grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
332 s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
333 dt_total = time.perf_counter() - t0
334 sid = f"vk-{tag}"
335 out_dir.mkdir(parents=True, exist_ok=True)
336 wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
337 entry = {
338 "solution_id": sid,
339 "output": f"{out_dir.name}/{sid}.wav",
340 "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",
341 "parameters": {
342 "grid_spacing_m": spacing,
343 "cells": model.plate.n_cells,
344 "n_dof": model.plate.ndof,
345 "transverse_modes": model.n_w,
346 "basis_even_in_y": model.symmetric_y,
347 "airy_modes": model.n_airy,
348 "oversampling": m,
349 "dt_s": 1.0 / (fs * m),
350 "highest_transverse_mode_hz": round(float(model.omega.max() / (2 * np.pi)), 1),
351 },
352 "timing_seconds": round(dt_total, 3),
353 "hardware": hardware_string(device),
354 "peak_si": float(np.max(np.abs(s))),
355 "setup_seconds": round(t1 - t0, 3),
356 }
357 if fields is not None:
358 import h5py
360 fg = ch["fields"]
361 with h5py.File(out_dir / f"{sid}.h5", "w") as fh:
362 fh.attrs["challenge_id"] = ch["challenge_id"]
363 fh.attrs["solution_id"] = sid
364 fh.attrs["quantity"] = fg["quantity"]
365 fh.attrs["units"] = fg["units"]
366 fh.create_dataset("x", data=grid[0])
367 fh.create_dataset("y", data=grid[1])
368 fh.create_dataset("t", data=grid[2][: fields.shape[0]])
369 fh.create_dataset("u", data=fields, compression="gzip", compression_opts=4, chunks=(8, fields.shape[1], fields.shape[2]))
370 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
371 meta.append(entry)
372 log(f"{tag}: done, setup={t1 - t0:.1f}s total={dt_total:.1f}s peak={np.max(np.abs(s)):.4e} Pa")
373 with open(out_dir / f"{sid}.meta.json", "w") as fh:
374 json.dump(entry, fh, indent=1)
375 print(json.dumps(meta, indent=1))
378if __name__ == "__main__":
379 main()