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