7141e28Add vonkarman-modal submission for plate-strike-nl-v1; viewer removes rigid motionJeremy Magland 1"""Von Karman (Airy stress function) thin-plate solution of plate-strike-nl-v1 by a modal method.
3The transverse displacement is expanded in the linear free-plate modes Phi_k computed with the
4Bogner-Fox-Schmit elements of the linear solver (../kirchhoff-bfs/kl_bfs_plate.py), and the Airy
5stress function in the clamped biharmonic eigenfunctions Psi_p on the same mesh. With
6mass-normalized Phi and L2-normalized Psi, the equations are
8 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,
9 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,
11where 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
12Gautschi-type scheme at a step of 1/(fs m): the damped linear part is propagated exactly (the
13recurrence coefficients are those of the exact damped oscillator), and the strike and the cubic
14coupling force are applied with a static-exact weight, which is second-order accurate and free of
15numerical dispersion in the linear limit. The Rayleigh integral is evaluated with
163 x 3 Gauss points per element; each point's retarded acceleration is a cubic Lagrange
17interpolation of the modal accelerations at neighbouring time steps, and the interpolation weights
18are folded into a per-step table, so the pressure at a step is a short sum over the recent modal
19accelerations. Damping acts on the bending motion only (the in-plane state is quasi-static),
20a further approximation of the challenge model.
22Usage: python vk_modal_plate.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [tag ...]
23 where a tag is n<mesh>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. n24-w150-a150-m2
24"""
25from __future__ import annotations
27import json
28import sys
29import time
30from pathlib import Path
32import numpy as np
33from numpy.polynomial import legendre as leg
34from scipy.io import wavfile
35from scipy.linalg import eigh
37sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "kirchhoff-bfs"))
38from kl_bfs_plate import BFSPlate, load_challenge # noqa: E402
41class VonKarmanModal:
42 def __init__(self, ch: dict, n: int, n_w: int, n_airy: int, force_override: float | None = None):
43 self.ch = ch
44 p = dict(ch["parameters"])
45 if force_override is not None:
46 p["strike_peak_force"] = force_override
47 ch2 = dict(ch)
48 ch2["parameters"] = p
49 self.p = p
50 self.plate = BFSPlate(ch2, n) # linear modes, modal forces, alpha, omega
51 pl = self.plate
52 # transverse modes: the n_w lowest (rigid ones included; they carry no membrane coupling)
53 self.n_w = min(n_w, pl.n_modes)
54 self.omega = pl.omega[: self.n_w]
55 self.alpha = pl.alpha[: self.n_w]
56 self.f = pl.f[: self.n_w]
57 V = pl.V[:, : self.n_w]
58 # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
59 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]])))
60 M1 = pl.M / pl.m # unit-density mass matrix
61 interior = np.setdiff1d(np.arange(pl.ndof), pl.boundary_dofs())
62 z4, psi_int = eigh(Kb[np.ix_(interior, interior)], M1[np.ix_(interior, interior)])
63 self.n_airy = min(n_airy, len(z4))
64 self.zeta4 = z4[: self.n_airy]
65 Psi = np.zeros((pl.ndof, self.n_airy))
66 Psi[interior] = psi_int[:, : self.n_airy]
67 # coupling tensor H[p, i, j] by 3 x 3 Gauss quadrature per element
68 g, wg = leg.leggauss(3)
69 s = (g + 1) / 2
70 ws = wg / 2
71 ex, ey = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
72 XI, ETA = np.meshgrid(s, s, indexing="ij")
73 xq = (ex.ravel()[:, None] + XI.ravel()[None, :]).ravel() * pl.a
74 yq = (ey.ravel()[:, None] + ETA.ravel()[None, :]).ravel() * pl.b
75 WQ = np.tile((ws[:, None] * ws[None, :]).ravel() * pl.a * pl.b, n * n)
76 Bxx, Byy, Bxy = pl.basis_derivs_at(xq, yq)
77 Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
78 PsiQ = pl.basis_at(xq, yq) @ Psi # (Q, n_airy)
79 H = np.zeros((self.n_airy, self.n_w, self.n_w))
80 for pi in range(self.n_airy):
81 wp = (WQ * PsiQ[:, pi])[:, None]
82 H[pi] = (Pxx * wp).T @ Pyy + (Pyy * wp).T @ Pxx - 2.0 * (Pxy * wp).T @ Pxy
83 self.H = H
84 self.Hflat = H.reshape(self.n_airy * self.n_w, self.n_w)
85 self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
86 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
87 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
88 self.tau = r / p["c_air"]
89 self.cq = p["rho_air"] / (2 * np.pi) * WQ / r
90 self.PhiQ = pl.basis_at(xq, yq) @ V # (Q, n_w)
92 def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
93 v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
94 eta = self.coef * (v @ q) # eta_p
95 return v.T @ eta # F_k = sum_p v[p,k] eta_p
97 def run(self, fs: int, nsamp: int, m: int, field_grid=None, field_rate: float | None = None):
98 p = self.p
99 dt = 1.0 / (fs * m)
100 T, Om = p["strike_duration"], 2 * np.pi / p["strike_duration"]
101 nsteps = nsamp * m
102 w2, al, f = self.omega**2, self.alpha, self.f
103 # exact two-step recurrence of the damped linear oscillator: q+ = a1 q + a2 q- + b F
104 rigid = self.omega == 0.0
105 disc = np.sqrt(al**2 - w2 + 0j)
106 E1, E2 = np.exp((-al + disc) * dt), np.exp((-al - disc) * dt)
107 a1 = (E1 + E2).real
108 a2 = -(E1 * E2).real
109 with np.errstate(divide="ignore", invalid="ignore"):
110 b = np.where(rigid, 0.0, (1.0 - a1 - a2) / np.where(rigid, 1.0, w2))
111 if rigid.any():
112 e2 = np.exp(-2.0 * al[rigid] * dt)
113 a1[rigid] = 1.0 + e2
114 a2[rigid] = -e2
115 b[rigid] = dt * (1.0 - e2) / (2.0 * al[rigid])
116 # fractional delays: each quadrature point's retarded acceleration is a 4-point Lagrange
117 # interpolation of the modal accelerations at steps d0-1 .. d0+2; the interpolation weights
118 # are folded into a per-step weight table G[step, mode]
119 tau_steps = self.tau / dt
120 d0 = np.floor(tau_steps).astype(int)
121 fr = tau_steps - d0
122 lag = np.stack([
123 -fr * (fr - 1.0) * (fr - 2.0) / 6.0, # weight of step d0-1 (interpolating at d0+fr)
124 (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0, # d0
125 -(fr + 1.0) * fr * (fr - 2.0) / 2.0, # d0+1
126 (fr + 1.0) * fr * (fr - 1.0) / 6.0, # d0+2
127 ], axis=1)
128 dmax = int(d0.max()) + 2
129 G = np.zeros((dmax + 1, self.n_w))
130 for j, off in enumerate([-1, 0, 1, 2]):
131 dd = d0 + off
132 for u in np.unique(dd):
133 sel = dd == u
134 G[u] += (self.cq[sel] * lag[sel, j]) @ self.PhiQ[sel]
135 ring = np.zeros((dmax + 1, self.n_w)) # ring[n mod (dmax+1)] = qdd at step n
136 out = np.zeros(nsamp)
137 q_prev = np.zeros(self.n_w)
138 q = np.zeros(self.n_w)
139 frames = None
140 if field_grid is not None:
141 fx, fy, ft = field_grid
142 XX, YY = np.meshgrid(fx, fy, indexing="xy")
143 PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.plate.V[:, : self.n_w]
144 frame_every = int(round(fs * m / field_rate))
145 frames = []
146 for nstep in range(nsteps + 1):
147 t = nstep * dt
148 psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
149 Fext = self.nonlinear_force(q) + f * psi # forces other than the linear restoring and damping terms
150 q_next = a1 * q + a2 * q_prev + b * Fext
151 # acceleration from the equation of motion (velocity by centered difference)
152 qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
153 ring[nstep % (dmax + 1)] = qdd
154 if nstep % m == 0: # the ring buffer starts at zero, which is the state before the strike
155 k = nstep // m
156 if k < nsamp:
157 idx = (nstep - np.arange(dmax + 1)) % (dmax + 1)
158 out[k] = np.sum(G * ring[idx])
159 if frames is not None and nstep % frame_every == 0 and len(frames) < len(ft):
160 frames.append((PhiG @ q).astype(np.float32))
161 q_prev, q = q, q_next
162 fields = None
163 if frames is not None:
164 fields = np.stack(frames).reshape(len(frames), len(field_grid[1]), len(field_grid[0]))
165 return out, fields
168def parse_tag(tag: str):
169 parts = dict((kv[0], int(kv[1:])) for kv in tag.split("-"))
170 return parts["n"], parts["w"], parts["a"], parts["m"]
173def main():
174 ch = load_challenge(Path(sys.argv[1]))
175 out_dir = Path(sys.argv[2])
176 args = sys.argv[3:]
177 fields_for = set()
178 force = None
179 if "--fields" in args:
180 i = args.index("--fields")
181 fields_for = set(args[i + 1].split(","))
182 args = args[:i] + args[i + 2:]
183 if "--force" in args:
184 i = args.index("--force")
185 force = float(args[i + 1])
186 args = args[:i] + args[i + 2:]
187 tags = args or ["n12-w50-a50-m1", "n16-w100-a100-m1", "n24-w150-a150-m2", "n24-w200-a200-m2"]
188 out = ch["output"]
189 fs, nsamp = out["sample_rate_hz"], out["num_samples"]
190 meta = []
191 for tag in tags:
192 n, n_w, n_a, m = parse_tag(tag)
193 t0 = time.perf_counter()
194 model = VonKarmanModal(ch, n, n_w, n_a, force)
195 t1 = time.perf_counter()
196 grid = None
197 if tag in fields_for:
198 fg = ch["fields"]
199 nt, ny, nx = fg["shape"]
200 grid = (np.linspace(0, ch["parameters"]["Lx"], nx), np.linspace(0, ch["parameters"]["Ly"], ny), np.arange(nt) / fg["frame_rate_hz"])
201 s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
202 dt_total = time.perf_counter() - t0
203 sid = f"vk-{tag}"
204 out_dir.mkdir(parents=True, exist_ok=True)
205 wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
206 entry = {
207 "solution_id": sid,
208 "output": f"{out_dir.name}/{sid}.wav",
209 "method": "von Karman plate (Airy form), modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
210 "parameters": {"elements_per_side": n, "transverse_modes": model.n_w, "airy_modes": model.n_airy, "oversampling": m, "dt_s": 1.0 / (fs * m), "highest_transverse_mode_hz": round(float(model.omega.max() / (2 * np.pi)), 1)},
211 "timing_seconds": round(dt_total, 3),
212 "peak_si": float(np.max(np.abs(s))),
213 }
214 if fields is not None:
215 import h5py
217 fg = ch["fields"]
218 with h5py.File(out_dir / f"{sid}.h5", "w") as fh:
219 fh.attrs["challenge_id"] = ch["challenge_id"]
220 fh.attrs["solution_id"] = sid
221 fh.attrs["quantity"] = fg["quantity"]
222 fh.attrs["units"] = fg["units"]
223 fh.create_dataset("x", data=grid[0])
224 fh.create_dataset("y", data=grid[1])
225 fh.create_dataset("t", data=grid[2][: fields.shape[0]])
226 fh.create_dataset("u", data=fields, compression="gzip", compression_opts=4, chunks=(8, fields.shape[1], fields.shape[2]))
227 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
228 meta.append(entry)
229 print(f"{tag}: modes w={model.n_w} a={model.n_airy} setup={t1 - t0:.1f}s total={dt_total:.1f}s peak={np.max(np.abs(s)):.4e} Pa", file=sys.stderr)
230 print(json.dumps(meta, indent=1))
233if __name__ == "__main__":
234 main()