Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflow
10 changed files+1237−0
.gitattributesadded+1−0View file
@@ -0,0 +1 @@
1+*.h5 filter=lfs diff=lfs merge=lfs -text
.mochi/workflows/results.ymladded+32−0View file
@@ -0,0 +1,32 @@
1+name: Validate and publish
2+on:
3+ push:
4+ workflow_dispatch:
5+concurrency:
6+ group: site
7+ cancel-in-progress: true
8+jobs:
9+ site:
10+ runs-on: ubuntu-latest
11+ steps:
12+ - uses: actions/checkout@v4
13+ - name: Fetch LFS objects (field exports), if git-lfs is available
14+ run: git lfs version && git lfs pull || echo "git-lfs unavailable or pull failed; field files may be pointers"
15+ - name: Get the challenge repository (specifications)
16+ run: git clone --depth 1 --branch main https://vault1.magland.org/challenges/cymbal challenge-repo
17+ - uses: actions/setup-python@v5
18+ with:
19+ python-version: '3.12'
20+ - run: pip install git+https://vault1.magland.org/challenges/challenge-tools
21+ - name: Validate submissions
22+ run: challenge-tools validate --challenges challenge-repo/challenges --solutions solutions
23+ - name: Build this repository's site
24+ run: python web/build.py --out site
25+ - if: github.ref == 'refs/heads/main'
26+ uses: actions/configure-pages@v4
27+ - if: github.ref == 'refs/heads/main'
28+ uses: actions/upload-pages-artifact@v3
29+ with:
30+ path: site
31+ - if: github.ref == 'refs/heads/main'
32+ uses: actions/deploy-pages@v4
README.mdadded+37−0View file
@@ -0,0 +1,37 @@
1+# Solutions to the cymbal challenges
2+
3+This repository holds solutions to the challenges defined in [challenges/cymbal](https://vault1.magland.org/challenges/cymbal). It is one of possibly several solution repositories registered with that challenge repository, which clones each of them, compares every solution against a common reference, and publishes the comparison at [challenges/cymbal/site](https://vault1.magland.org/challenges/cymbal/site/). This repository's own site, at [challenges/cymbal-solutions/site](https://vault1.magland.org/challenges/cymbal-solutions/site/), lists its submissions with audio players and animates the exported displacement fields as a 3D surface.
4+
5+## Layout
6+
7+```
8+solutions/<challenge-id>/<submission-id>.json one submission: metadata and one or more solutions
9+solutions/<challenge-id>/<submission-id>/*.wav the scored output files, referenced from the JSON
10+solutions/<challenge-id>/<submission-id>/*.h5 optional displacement fields for visualization (git LFS)
11+solvers/ source code of the submissions
12+web/ this repository's site (built by web/build.py)
13+```
14+
15+Field exports are about 40 MB each and are stored with git LFS (`.gitattributes` routes `*.h5` there). Cloning without git-lfs leaves pointer files in their place; the validator tolerates that and the site marks such fields as unavailable.
16+
17+## The example submission
18+
19+`vonkarman-modal` solves the von Kármán thin-plate reduction of the challenge model (Airy stress function, in-plane inertia neglected, damping on bending only) by a modal method on a Cartesian grid. The disk is approximated by the square cells of a uniform grid whose centers lie inside it, a staircase boundary, with grid spacings of 8 mm and 4 mm. The linear free-plate modes of that domain and the clamped Airy modes come from Bogner-Fox-Schmit C1 finite elements on the cells and a sparse shift-invert eigensolver; the cubic coupling tensor is integrated on the same cells; the time stepping is a two-step scheme exact in the linear part; and the retarded Rayleigh integral is evaluated with cubic fractional delays. The chain is ordered by grid spacing, number of modes, and time step. It is a thin-plate approximation on an approximate domain, so its converged limit differs from the three-dimensional model that defines the challenge.
20+
21+```
22+bash solvers/vonkarman-modal-disk/run_chain.sh <cymbal-repo>/challenges/cymbal-flat-disk-v1/challenge.json
23+```
24+
25+The chain is `g8-w100-a100-m1`, `g8-w200-a200-m1`, `g4-w200-a200-m2`, `g4-w300-a300-m2` (grid spacing in mm, transverse modes, Airy modes, time steps per output sample). Because the strike is centered on the x axis the solution is even in y, and the modal bases are restricted to that symmetric subspace, which is exact and halves the number of modes needed for a given frequency reach: 200 symmetric modes reach about 6 kHz and 300 about 9 kHz. The 8 mm grid has 2024 cells and 8500 degrees of freedom, the 4 mm grid 8112 cells and 33268. On the 8 mm grid the lowest elastic frequencies are 21.76 and 21.88 Hz for the two-nodal-diameter pair (21.87 Hz analytically), 38.01 Hz for the first axisymmetric mode (38.03) and 87.6 and 90.2 Hz for the four-nodal-diameter pair (89.61); on the 4 mm grid the last pair is 89.03 and 89.19 Hz. The staircase splits degenerate pairs by about a percent at 8 mm and a few tenths of a percent at 4 mm. The linear limit of the nonlinear solver agrees with the closed-form linear modal solution on the same mesh to 5e-4 at one step per sample and 1.3e-4 at two.
26+
27+## Adding a submission
28+
29+Follow the same steps as in [challenges/string-solutions](https://vault1.magland.org/challenges/string-solutions#adding-a-submission-to-this-repository): produce the output exactly as the [specification](https://vault1.magland.org/challenges/cymbal/blob/main/challenges/cymbal-flat-disk-v1/README.md) says, write `solutions/<challenge-id>/<submission-id>.json` following the [schema](https://vault1.magland.org/challenges/challenge-tools/blob/main/schema.json), record known accuracy orderings in `more_accurate_than`, and open a pull request.
30+
31+## Validating locally
32+
33+```
34+pip install git+https://vault1.magland.org/challenges/challenge-tools
35+challenge-tools validate --challenges <path-to-cymbal-repo>/challenges --solutions solutions
36+python web/build.py --out site
37+```
solvers/vonkarman-modal-disk/bfs_disk.pyadded+428−0View file
@@ -0,0 +1,428 @@
1+"""Kirchhoff-Love plate on a Cartesian grid masked to a disk, with Bogner-Fox-Schmit elements.
2+
3+The disk of radius R is approximated by the union of the square cells of a uniform Cartesian grid
4+of spacing a whose centers lie inside the disk (a staircase boundary; the grid is centered on the
5+disk so the approximation has the symmetries of the square lattice). Each cell is a C1-conforming
6+bicubic Hermite rectangle with four degrees of freedom per node (w, w_x, w_y, w_xy). Free edges,
7+including the staircase corners, are natural boundary conditions of the Kirchhoff energy, so
8+nothing is imposed. Element matrices are integrated exactly with 4 x 4 Gauss-Legendre points and
9+assembled into sparse matrices; the lowest modes come from a shift-invert Lanczos eigensolver
10+(ARPACK through scipy) and are mass-normalized. Three of them are rigid-body modes at zero
11+frequency.
12+
13+The class also provides the closed-form linear modal response to the challenge's Hann strike
14+(a particular solution during the pulse and free decay afterwards) and an exact evaluation of the
15+retarded Rayleigh integral, both taken from the square-plate solver of challenges/plate-solutions,
16+so that the linear limit of the nonlinear solver can be checked on the same mesh.
17+
18+Coordinates are those of the challenge: the disk is centered at the origin, w is positive in +z,
19+the strike enters as -p, and w is reported directly as u_z.
20+"""
21+from __future__ import annotations
22+
23+import json
24+from pathlib import Path
25+
26+import numpy as np
27+import scipy.sparse as sp
28+from numpy.polynomial import legendre as leg
29+from scipy.sparse.linalg import eigsh
30+from scipy.spatial import cKDTree
31+
32+RIGID_OMEGA = 1.0 # rad/s; eigenvalues below this are the three rigid-body modes
33+
34+
35+def load_challenge(path: Path) -> dict:
36+ with open(path) as f:
37+ return json.load(f)
38+
39+
40+# ---------------------------------------------------------------- Hermite cubics on [0, 1]
41+def hermite(s: np.ndarray):
42+ """Values and first two derivatives (w.r.t. s) of the four cubic Hermite functions. Shapes (len(s), 4)."""
43+ s = np.asarray(s, dtype=float)
44+ H = np.stack([1 - 3 * s**2 + 2 * s**3, s - 2 * s**2 + s**3, 3 * s**2 - 2 * s**3, -(s**2) + s**3], axis=-1)
45+ dH = np.stack([-6 * s + 6 * s**2, 1 - 4 * s + 3 * s**2, 6 * s - 6 * s**2, -2 * s + 3 * s**2], axis=-1)
46+ d2H = np.stack([-6 + 12 * s, -4 + 6 * s, 6 - 12 * s, -2 + 6 * s], axis=-1)
47+ return H, dH, d2H
48+
49+
50+CORNERS = [(0, 0), (1, 0), (0, 1), (1, 1)]
51+
52+
53+def shape_functions(xi: np.ndarray, eta: np.ndarray, a: float):
54+ """BFS shape functions and second derivatives at local coordinates (xi, eta) for square a x a elements.
55+
56+ Returns N, N_xx, N_yy, N_xy, each of shape (npts, 16). Local dof order: for corner c in
57+ [(0,0), (1,0), (0,1), (1,1)]: w, w_x, w_y, w_xy.
58+ """
59+ Hx, dHx, d2Hx = hermite(xi)
60+ Hy, dHy, d2Hy = hermite(eta)
61+ N = np.zeros((len(xi), 16))
62+ Nxx = np.zeros_like(N)
63+ Nyy = np.zeros_like(N)
64+ Nxy = np.zeros_like(N)
65+ for c, (cx, cy) in enumerate(CORNERS):
66+ ix0, ix1 = (0, 1) if cx == 0 else (2, 3)
67+ iy0, iy1 = (0, 1) if cy == 0 else (2, 3)
68+ combos = [(ix0, iy0, 1.0), (ix1, iy0, a), (ix0, iy1, a), (ix1, iy1, a * a)]
69+ for d, (fx, fy, sc) in enumerate(combos):
70+ k = 4 * c + d
71+ N[:, k] = sc * Hx[:, fx] * Hy[:, fy]
72+ Nxx[:, k] = sc * d2Hx[:, fx] * Hy[:, fy] / a**2
73+ Nyy[:, k] = sc * Hx[:, fx] * d2Hy[:, fy] / a**2
74+ Nxy[:, k] = sc * dHx[:, fx] * dHy[:, fy] / (a * a)
75+ return N, Nxx, Nyy, Nxy
76+
77+
78+class BFSDisk:
79+ def __init__(self, ch: dict, spacing: float, n_modes: int, symmetric_y: bool = False, max_hz: float | None = None):
80+ """symmetric_y restricts the basis to functions even in y (w(x, -y) = w(x, y)), which is exact for a strike
81+ centered on the x axis; the grid is symmetric about that axis by construction. max_hz drops modes above it."""
82+ p = ch["parameters"]
83+ self.p = p
84+ R, h, nu = p["disk_radius"], p["h"], p["nu"]
85+ self.R = R
86+ a = float(spacing)
87+ self.a = a
88+ D = p["E"] * h**3 / (12.0 * (1.0 - nu**2))
89+ m = p["rho"] * h
90+ self.D, self.m = D, m
91+
92+ # ---- the masked grid
93+ half = int(np.ceil(R / a)) + 1
94+ ncell = 2 * half
95+ self.ncell = ncell
96+ self.half = half
97+ self.x0 = -half * a # coordinate of grid line 0
98+ centers = self.x0 + (np.arange(ncell) + 0.5) * a
99+ CX, CY = np.meshgrid(centers, centers, indexing="ij")
100+ active = CX**2 + CY**2 <= R**2
101+ cells = np.argwhere(active) # (nc, 2): ix, iy
102+ self.cells = cells
103+ self.n_cells = len(cells)
104+ self.cell_id = -np.ones((ncell, ncell), dtype=int)
105+ self.cell_id[cells[:, 0], cells[:, 1]] = np.arange(self.n_cells)
106+ self.cell_center = np.stack([centers[cells[:, 0]], centers[cells[:, 1]]], axis=1)
107+ self._tree = cKDTree(self.cell_center)
108+ nn = ncell + 1
109+ gnode = np.zeros((self.n_cells, 4), dtype=int)
110+ for c, (cx, cy) in enumerate(CORNERS):
111+ gnode[:, c] = (cells[:, 0] + cx) * nn + (cells[:, 1] + cy)
112+ touch = np.bincount(gnode.ravel(), minlength=nn * nn)
113+ used = touch > 0
114+ self.node_map = -np.ones(nn * nn, dtype=int)
115+ self.node_map[used] = np.arange(int(used.sum()))
116+ self.n_nodes = int(used.sum())
117+ self.ndof = 4 * self.n_nodes
118+ cnode = self.node_map[gnode] # (nc, 4) compact node ids
119+ gidx = np.flatnonzero(used)
120+ self.node_i, self.node_j = gidx // nn, gidx % nn # grid indices of each compact node
121+ self.cell_dofs = (4 * cnode[:, :, None] + np.arange(4)[None, None, :]).reshape(self.n_cells, 16)
122+ bnodes = self.node_map[used & (touch < 4)]
123+ self.boundary_dofs = (4 * bnodes[:, None] + np.arange(4)[None, :]).ravel()
124+ self.interior_dofs = np.setdiff1d(np.arange(self.ndof), self.boundary_dofs)
125+ self.area = self.n_cells * a * a
126+
127+ # ---- element matrices
128+ g, wg = leg.leggauss(4)
129+ s = (g + 1) / 2
130+ ws = wg / 2
131+ XI, ETA = np.meshgrid(s, s, indexing="ij")
132+ WQ = (ws[:, None] * ws[None, :]).ravel() * a * a
133+ N, Nxx, Nyy, Nxy = shape_functions(XI.ravel(), ETA.ravel(), a)
134+ self._shape = (N, Nxx, Nyy, Nxy, WQ)
135+ C = np.array([[1.0, nu, 0.0], [nu, 1.0, 0.0], [0.0, 0.0, 2.0 * (1.0 - nu)]])
136+ Ke = D * self.element_matrix(C)
137+ Me = np.zeros((16, 16))
138+ for q in range(len(WQ)):
139+ Me += WQ[q] * m * np.outer(N[q], N[q])
140+ self.K = self.assemble(Ke)
141+ self.M = self.assemble(Me)
142+
143+ # ---- modes: the n_modes lowest, mass-normalized (shift-invert Lanczos around a negative shift)
144+ self.symmetric_y = symmetric_y
145+ evals, evecs = self.eigenmodes(self.K, self.M, n_modes, sigma=-1.0e3)
146+ omega = np.sqrt(np.maximum(evals, 0.0))
147+ omega[omega < RIGID_OMEGA] = 0.0
148+ if max_hz is not None:
149+ keep = omega <= 2 * np.pi * max_hz
150+ omega, evecs = omega[keep], evecs[:, keep]
151+ self.omega = omega
152+ self.V = evecs
153+ self.n_modes = len(omega)
154+ self.alpha = p["sigma0"] + 0.5 * p["eta"] * omega**2
155+ self.rigid = omega == 0.0
156+ self.n_rigid = int(self.rigid.sum())
157+
158+ # ---- modal forces: f_n = -int p_hat Phi_n dA over the strike disk (polar Gauss quadrature)
159+ Rs, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
160+ c_phi = 0.5 - 2.0 / np.pi**2
161+ gr, gw = leg.leggauss(24)
162+ r = Rs * (gr + 1.0) / 2.0
163+ wr = gw * Rs / 2.0
164+ nth = 64
165+ th = 2 * np.pi * np.arange(nth) / nth
166+ rr, tt = np.meshgrid(r, th, indexing="ij")
167+ ww = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
168+ px = (xs + rr * np.cos(tt)).ravel()
169+ py = (ys + rr * np.sin(tt)).ravel()
170+ phat = p["strike_peak_force"] / (c_phi * np.pi * Rs**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / Rs))
171+ self.f = -(ww * phat) @ self.modes_at(px, py)
172+
173+ # ---- closed-form linear response constants (see the square-plate solver)
174+ T = p["strike_duration"]
175+ Om = 2 * np.pi / T
176+ self.T, self.Om = T, Om
177+ al, w, f = self.alpha, self.omega, self.f
178+ with np.errstate(divide="ignore", invalid="ignore"):
179+ disc = np.sqrt(al**2 - w**2 + 0j)
180+ self.s1 = -al + disc
181+ self.s2 = -al - disc
182+ H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
183+ self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
184+ self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
185+ rhs1 = -(self.qconst - self.Hc.real)
186+ rhs2 = (1j * Om * self.Hc).real
187+ self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
188+ self.A = rhs1 - self.B
189+ self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
190+ self.vconst = f / (4.0 * al)
191+ self.Ar = -(self.vconst - self.Gc.real)
192+ qT, vT, _ = self._pulse_phase(np.array([T]))
193+ self.qT, self.vT = qT[:, 0], vT[:, 0]
194+ with np.errstate(divide="ignore", invalid="ignore"):
195+ self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
196+ self.A2 = self.qT - self.B2
197+
198+ # ------------------------------------------------------------ symmetry and eigensolves
199+ def symmetric_prolongation(self, dofs: np.ndarray | None = None) -> sp.csc_matrix:
200+ """Sparse matrix whose columns span the functions even in y: each column pairs a node with its mirror image
201+ (w and w_x equal, w_y and w_xy opposite), and nodes on the axis keep only w and w_x. If dofs is given, the
202+ matrix is restricted to those rows (in that order) and to the columns that touch them."""
203+ j_axis = self.half
204+ gi, gj = self.node_i, self.node_j
205+ nn = self.ncell + 1
206+ mirror = self.node_map[gi * nn + (2 * j_axis - gj)]
207+ rows, cols, vals = [], [], []
208+ col = 0
209+ signs = (1.0, 1.0, -1.0, -1.0)
210+ for n in range(self.n_nodes):
211+ if gj[n] > j_axis:
212+ nm = mirror[n]
213+ for d in range(4):
214+ rows += [4 * n + d, 4 * nm + d]
215+ cols += [col, col]
216+ vals += [1.0, signs[d]]
217+ col += 1
218+ elif gj[n] == j_axis:
219+ for d in range(2):
220+ rows.append(4 * n + d)
221+ cols.append(col)
222+ vals.append(1.0)
223+ col += 1
224+ P = sp.csc_matrix((vals, (rows, cols)), shape=(self.ndof, col))
225+ if dofs is not None:
226+ P = P[dofs]
227+ P = P[:, np.unique(P.nonzero()[1])]
228+ return P.tocsc()
229+
230+ def eigenmodes(self, K, M, k: int, sigma: float, dofs: np.ndarray | None = None):
231+ """The k lowest eigenpairs of K v = lambda M v, restricted to the dofs given (default all) and to the
232+ y-symmetric subspace if the mesh was built with symmetric_y. Sorted ascending, M-normalized, and returned
233+ on the full dof set (zero on excluded dofs)."""
234+ Kr, Mr = (K, M) if dofs is None else (K[dofs][:, dofs], M[dofs][:, dofs])
235+ P = self.symmetric_prolongation(dofs) if self.symmetric_y else None
236+ if P is not None:
237+ Kr, Mr = P.T @ Kr @ P, P.T @ Mr @ P
238+ Kr, Mr = sp.csc_matrix(Kr), sp.csc_matrix(Mr)
239+ evals, evecs = eigsh(Kr, k=min(k, Kr.shape[0] - 2), M=Mr, sigma=sigma, which="LM")
240+ order = np.argsort(evals)
241+ evals, evecs = evals[order], evecs[:, order]
242+ if P is not None:
243+ evecs = P @ evecs
244+ Md = M if dofs is None else M[dofs][:, dofs]
245+ evecs = evecs / np.sqrt(np.einsum("ij,ij->j", evecs, Md @ evecs))
246+ if dofs is not None:
247+ full = np.zeros((self.ndof, evecs.shape[1]))
248+ full[dofs] = evecs
249+ evecs = full
250+ return evals, evecs
251+
252+ # ------------------------------------------------------------ assembly
253+ def element_matrix(self, Cmat: np.ndarray) -> np.ndarray:
254+ """Element matrix of int [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA (no stiffness factor)."""
255+ _, Nxx, Nyy, Nxy, WQ = self._shape
256+ Ae = np.zeros((16, 16))
257+ for q in range(len(WQ)):
258+ B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]])
259+ Ae += WQ[q] * (B2.T @ Cmat @ B2)
260+ return Ae
261+
262+ def assemble(self, Ae: np.ndarray) -> sp.csc_matrix:
263+ rows = np.repeat(self.cell_dofs, 16, axis=1).ravel()
264+ cols = np.tile(self.cell_dofs, (1, 16)).ravel()
265+ vals = np.tile(Ae.ravel(), self.n_cells)
266+ A = sp.coo_matrix((vals, (rows, cols)), shape=(self.ndof, self.ndof)).tocsc()
267+ A.sum_duplicates()
268+ return A
269+
270+ # ------------------------------------------------------------ quadrature over the masked domain
271+ def quadrature(self, order: int = 3):
272+ """Gauss-Legendre points and weights over all active cells: xq, yq, wq."""
273+ g, wg = leg.leggauss(order)
274+ s = (g + 1) / 2
275+ ws = wg / 2
276+ XI, ETA = np.meshgrid(s, s, indexing="ij")
277+ a = self.a
278+ xq = (self.x0 + (self.cells[:, 0][:, None] + XI.ravel()[None, :]) * a).ravel()
279+ yq = (self.x0 + (self.cells[:, 1][:, None] + ETA.ravel()[None, :]) * a).ravel()
280+ wq = np.tile((ws[:, None] * ws[None, :]).ravel() * a * a, self.n_cells)
281+ return xq, yq, wq
282+
283+ # ------------------------------------------------------------ evaluation at points
284+ def locate(self, x: np.ndarray, y: np.ndarray, tol: float = 1e-9):
285+ """Cell id and local coordinates of points. Points inside the disk that fall in an inactive cell of the
286+ staircase are assigned the nearest active cell (the element polynomial is extrapolated slightly).
287+ Returns cid (-1 for points outside the disk), xi, eta."""
288+ x = np.asarray(x, dtype=float)
289+ y = np.asarray(y, dtype=float)
290+ ix = np.clip(np.floor((x - self.x0) / self.a).astype(int), 0, self.ncell - 1)
291+ iy = np.clip(np.floor((y - self.x0) / self.a).astype(int), 0, self.ncell - 1)
292+ cid = self.cell_id[ix, iy].copy()
293+ inside = x**2 + y**2 <= self.R**2 * (1 + tol)
294+ fix = inside & (cid < 0)
295+ if fix.any():
296+ _, near = self._tree.query(np.stack([x[fix], y[fix]], axis=1))
297+ cid[fix] = near
298+ cid[~inside] = -1
299+ c = np.maximum(cid, 0)
300+ cx = self.x0 + self.cells[c, 0] * self.a
301+ cy = self.x0 + self.cells[c, 1] * self.a
302+ return cid, (x - cx) / self.a, (y - cy) / self.a
303+
304+ def _sparse_basis(self, vals: np.ndarray, cid: np.ndarray) -> sp.csr_matrix:
305+ npts = len(cid)
306+ ok = cid >= 0
307+ rows = np.repeat(np.arange(npts)[ok], 16)
308+ cols = self.cell_dofs[cid[ok]].ravel()
309+ data = vals[ok].ravel()
310+ return sp.csr_matrix((data, (rows, cols)), shape=(npts, self.ndof))
311+
312+ def basis_at(self, x, y) -> sp.csr_matrix:
313+ """Sparse (npts x ndof) matrix of shape-function values; rows of points outside the disk are zero."""
314+ cid, xi, eta = self.locate(x, y)
315+ N, _, _, _ = shape_functions(xi, eta, self.a)
316+ return self._sparse_basis(N, cid)
317+
318+ def basis_derivs_at(self, x, y):
319+ cid, xi, eta = self.locate(x, y)
320+ _, Nxx, Nyy, Nxy = shape_functions(xi, eta, self.a)
321+ return self._sparse_basis(Nxx, cid), self._sparse_basis(Nyy, cid), self._sparse_basis(Nxy, cid)
322+
323+ def modes_at(self, x, y) -> np.ndarray:
324+ return self.basis_at(x, y) @ self.V
325+
326+ def inside(self, x, y) -> np.ndarray:
327+ return np.asarray(x) ** 2 + np.asarray(y) ** 2 <= self.R**2 * (1 + 1e-9)
328+
329+ # ------------------------------------------------------------ modal time response (linear)
330+ def _pulse_phase(self, t: np.ndarray):
331+ t = t[None, :]
332+ Om = self.Om
333+ e = np.exp(1j * Om * t)
334+ s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
335+ E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
336+ Hc = self.Hc[:, None]
337+ q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
338+ v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
339+ acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
340+ rg = self.rigid
341+ if rg.any():
342+ al = self.alpha[rg][:, None]
343+ Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
344+ Ed = np.exp(-2.0 * al * t)
345+ q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
346+ v[rg] = Ar * Ed + vc - (Gc * e).real
347+ acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
348+ return q, v, acc
349+
350+ def _free_phase(self, tau: np.ndarray):
351+ tau = tau[None, :]
352+ s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
353+ E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
354+ q = (A * E1 + B * E2).real
355+ v = (A * s1 * E1 + B * s2 * E2).real
356+ acc = (A * s1**2 * E1 + B * s2**2 * E2).real
357+ rg = self.rigid
358+ if rg.any():
359+ al = self.alpha[rg][:, None]
360+ qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
361+ Ed = np.exp(-2.0 * al * tau)
362+ q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
363+ v[rg] = vT * Ed
364+ acc[rg] = -2.0 * al * vT * Ed
365+ return q, v, acc
366+
367+ def response(self, t: np.ndarray):
368+ q = np.zeros((self.n_modes, len(t)))
369+ v = np.zeros_like(q)
370+ acc = np.zeros_like(q)
371+ m1 = (t >= 0) & (t <= self.T)
372+ m2 = t > self.T
373+ if m1.any():
374+ q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
375+ if m2.any():
376+ q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
377+ return q, v, acc
378+
379+ # ------------------------------------------------------------ Rayleigh integral (linear, exact retardation)
380+ def listener_quadrature(self, order: int = 3):
381+ p = self.p
382+ xq, yq, wq = self.quadrature(order)
383+ r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
384+ tau = r / p["c_air"]
385+ cq = p["rho_air"] / (2 * np.pi) * wq / r
386+ return xq, yq, tau, cq
387+
388+ def pressure(self, fs: int, nsamp: int) -> np.ndarray:
389+ xq, yq, tau, cq = self.listener_quadrature()
390+ Phi = self.modes_at(xq, yq) # (nq, n_modes)
391+ t = np.arange(nsamp) / fs
392+ out = np.zeros(nsamp)
393+ t_split = self.T + tau.max()
394+ early = (t < t_split) & (t >= tau.min())
395+ late = t >= t_split
396+ ne = int(early.sum())
397+ if ne:
398+ chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
399+ for i0 in range(0, len(xq), chunk):
400+ sl = slice(i0, i0 + chunk)
401+ tt = (t[early][None, :] - tau[sl][:, None]).ravel()
402+ _, _, acc = self.response(tt)
403+ acc = acc.reshape(self.n_modes, -1, ne)
404+ out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
405+ if late.any():
406+ tl = t[late] - self.T
407+ rg = self.rigid
408+ G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
409+ G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
410+ c1 = self.A2 * self.s1**2 * G1
411+ c2 = self.B2 * self.s2**2 * G2
412+ s1, s2 = self.s1.copy(), self.s2.copy()
413+ if rg.any():
414+ al = self.alpha[rg]
415+ Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
416+ c1[rg] = -2.0 * al * self.vT[rg] * Gr
417+ c2[rg] = 0.0
418+ s1[rg] = -2.0 * al
419+ s2[rg] = 0.0
420+ coef = np.concatenate([c1, c2])
421+ rates = np.concatenate([s1, s2])
422+ res = np.zeros(len(tl))
423+ step = max(1, 20_000_000 // len(coef))
424+ for k0 in range(0, len(tl), step):
425+ tk = tl[k0:k0 + step]
426+ res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
427+ out[late] = res
428+ return out
solvers/vonkarman-modal-disk/make_submission.pyadded+58−0View file
@@ -0,0 +1,58 @@
1+"""Assemble the submission JSON from the per-solution meta files written by vk_modal_disk.py.
2+
3+Usage: python make_submission.py <solutions/<challenge-id>/<submission-id>> tag1 tag2 ... (in order of increasing accuracy)
4+"""
5+import json
6+import sys
7+from pathlib import Path
8+
9+HARDWARE = "Intel Core i7-1355U laptop (2 performance + 8 efficiency cores, 40 GB), single process, numpy/scipy with OpenBLAS"
10+
11+
12+def main():
13+ out_dir = Path(sys.argv[1])
14+ tags = sys.argv[2:]
15+ challenge_id = out_dir.parent.name
16+ submission_id = out_dir.name
17+ sols = []
18+ prev = None
19+ for tag in tags:
20+ m = json.load(open(out_dir / f"vk-{tag}.meta.json"))
21+ entry = {
22+ "solution_id": m["solution_id"],
23+ "output": m["output"],
24+ "method": m["method"],
25+ "parameters": m["parameters"],
26+ "timing_seconds": m["timing_seconds"],
27+ "more_accurate_than": [prev] if prev else [],
28+ "hardware": HARDWARE,
29+ }
30+ if "fields" in m:
31+ entry["fields"] = m["fields"]
32+ sols.append(entry)
33+ prev = m["solution_id"]
34+ sub = {
35+ "schema_version": 1,
36+ "submission_id": submission_id,
37+ "challenge_id": challenge_id,
38+ "submitter": {"name": "Jeremy Magland", "email": "jmagland@flatironinstitute.org"},
39+ "date": "2026-09-17",
40+ "description": (
41+ "Von Karman thin-plate reduction (Airy stress function, in-plane inertia neglected, damping on bending only) of the challenge model "
42+ "on a Cartesian grid: the disk is approximated by the square cells whose centers lie inside it (a staircase boundary) at 8 mm and 4 mm spacing. "
43+ "Linear free-plate modes and clamped Airy modes come from Bogner-Fox-Schmit finite elements on the cells and a sparse shift-invert eigensolver; "
44+ "the bases are restricted to functions even in y, which is exact because the strike is centered on the x axis. The cubic coupling tensor is integrated "
45+ "on the cells, the time stepping is a two-step scheme exact in the linear part, and the retarded Rayleigh integral uses cubic fractional delays. "
46+ "The chain is ordered by grid spacing, number of modes, and time step; the mode truncation caps the frequency the cascade can reach."
47+ ),
48+ "code_url": "https://vault1.magland.org/challenges/cymbal-solutions/blob/main/solvers/vonkarman-modal-disk/vk_modal_disk.py",
49+ "viz_url": "https://vault1.magland.org/challenges/cymbal-solutions/site/",
50+ "solutions": sols,
51+ }
52+ with open(out_dir.parent / f"{submission_id}.json", "w") as f:
53+ json.dump(sub, f, indent=1)
54+ print(f"wrote {out_dir.parent / (submission_id + '.json')} with {len(sols)} solutions")
55+
56+
57+if __name__ == "__main__":
58+ main()
solvers/vonkarman-modal-disk/run_chain.shadded+18−0View file
@@ -0,0 +1,18 @@
1+#!/bin/bash
2+# Produce the vonkarman-modal chain for cymbal-flat-disk-v1 (run locally; hours of compute).
3+# Usage: bash run_chain.sh <cymbal-repo>/challenges/cymbal-flat-disk-v1/challenge.json
4+set -e
5+here="$(cd "$(dirname "$0")" && pwd)"
6+ch="$1"
7+out="$here/../../solutions/cymbal-flat-disk-v1/vonkarman-modal"
8+mkdir -p "$out"
9+cd "$here"
10+# the two cheap runs in parallel, then the two 4 mm runs in parallel (memory bandwidth is the limit, so no more than two at once)
11+OMP_NUM_THREADS=4 python3 vk_modal_disk.py "$ch" "$out" g8-w100-a100-m1 > "$out/g8-w100.log" 2>&1 &
12+OMP_NUM_THREADS=4 python3 vk_modal_disk.py "$ch" "$out" g8-w200-a200-m1 > "$out/g8-w200.log" 2>&1 &
13+wait
14+OMP_NUM_THREADS=4 python3 vk_modal_disk.py "$ch" "$out" g4-w200-a200-m2 > "$out/g4-w200.log" 2>&1 &
15+OMP_NUM_THREADS=6 python3 vk_modal_disk.py "$ch" "$out" --fields g4-w300-a300-m2 g4-w300-a300-m2 > "$out/g4-w300.log" 2>&1 &
16+wait
17+python3 make_submission.py "$out" g8-w100-a100-m1 g8-w200-a200-m1 g4-w200-a200-m2 g4-w300-a300-m2
18+rm -f "$out"/*.log
solvers/vonkarman-modal-disk/vk_modal_disk.pyadded+274−0View file
@@ -0,0 +1,274 @@
1+"""Von Karman (Airy stress function) solution of cymbal-flat-disk-v1 by a modal method on a Cartesian grid.
2+
3+The disk is approximated by the square cells of a uniform Cartesian grid whose centers lie inside
4+it (see bfs_disk.py). The transverse displacement is expanded in the linear free-plate modes Phi_k
5+of that staircase domain, computed with Bogner-Fox-Schmit elements, and the Airy stress function
6+in the clamped biharmonic eigenfunctions Psi_p on the same mesh. With mass-normalized Phi and
7+L2-normalized Psi the equations are
8+
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,
11+
12+where 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
13+Gautschi-type scheme at a step of 1/(fs m): the damped linear part is propagated exactly and the
14+strike and the cubic coupling force are applied with a static-exact weight, which is second-order
15+accurate and free of numerical dispersion in the linear limit. The Rayleigh integral is evaluated
16+with 3 x 3 Gauss points per cell; each point's retarded acceleration is a cubic Lagrange
17+interpolation of the modal accelerations at neighbouring time steps, folded into a per-step
18+weight table. Damping acts on the bending motion only (the in-plane state is quasi-static),
19+a further approximation of the challenge model.
20+
21+The 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
23+cascade can reach.
24+
25+Because 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
27+modes needed for a given frequency reach (the coupling cost scales with the cube of the count). --full
28+uses the unrestricted bases instead.
29+
30+Usage: 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+"""
33+from __future__ import annotations
34+
35+import json
36+import sys
37+import time
38+from pathlib import Path
39+
40+import numpy as np
41+from scipy.io import wavfile
42+sys.path.insert(0, str(Path(__file__).resolve().parent))
43+from bfs_disk import BFSDisk, load_challenge # noqa: E402
44+
45+
46+class 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)
95+
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
100+
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
179+
180+
181+def 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"]
184+
185+
186+def 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
193+
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
202+
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
254+
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))
271+
272+
273+if __name__ == "__main__":
274+ main()
web/build.pyadded+63−0View file
@@ -0,0 +1,63 @@
1+"""Assemble this repository's static site: a manifest of the submissions, the solution files, and the viewer page."""
2+from __future__ import annotations
3+
4+import argparse
5+import json
6+import shutil
7+from pathlib import Path
8+
9+ROOT = Path(__file__).resolve().parents[1]
10+LFS_PREFIX = b"version https://git-lfs.github.com/spec/v1"
11+
12+
13+def is_lfs_pointer(p: Path) -> bool:
14+ try:
15+ with open(p, "rb") as f:
16+ return f.read(len(LFS_PREFIX)) == LFS_PREFIX
17+ except OSError:
18+ return False
19+
20+
21+def main():
22+ ap = argparse.ArgumentParser()
23+ ap.add_argument("--out", type=Path, default=ROOT / "site")
24+ args = ap.parse_args()
25+ out: Path = args.out
26+ if out.exists():
27+ shutil.rmtree(out)
28+ out.mkdir(parents=True)
29+ manifest = {"challenges": {}}
30+ for p in sorted((ROOT / "solutions").glob("*/*.json")):
31+ sub = json.loads(p.read_text())
32+ cid = p.parent.name
33+ entry = manifest["challenges"].setdefault(cid, {"submissions": []})
34+ rel_dir = p.parent.relative_to(ROOT)
35+ for sol in sub["solutions"]:
36+ sol["output"] = f"{rel_dir.as_posix()}/{sol['output']}"
37+ if "fields" in sol and not sol["fields"]["path"].startswith("http"):
38+ local = p.parent / sol["fields"]["path"]
39+ sol["fields"]["path"] = f"{rel_dir.as_posix()}/{sol['fields']['path']}"
40+ sol["fields"]["available"] = local.exists() and not is_lfs_pointer(local)
41+ entry["submissions"].append(sub)
42+ shutil.copytree(ROOT / "solutions", out / "solutions")
43+ shutil.copyfile(ROOT / "web" / "index.html", out / "index.html")
44+ shutil.copyfile(ROOT / "web" / "viewer.js", out / "viewer.js")
45+ # the shared Web Audio player, from the installed challenge-tools package
46+ from challenge_tools.build import player_js
47+
48+ (out / "player.js").write_text(player_js())
49+ # suggested playback gains and the challenge parameters (strike point, radius), when the challenge repository is checked out beside us (as in the workflow)
50+ gains = {}
51+ params = {}
52+ for cj in sorted((ROOT / "challenge-repo" / "challenges").glob("*/challenge.json")) if (ROOT / "challenge-repo").exists() else []:
53+ ch = json.loads(cj.read_text())
54+ gains[ch["challenge_id"]] = (ch.get("playback") or {}).get("gain_db")
55+ params[ch["challenge_id"]] = ch.get("parameters", {})
56+ manifest["playback_gain_db"] = gains
57+ manifest["parameters"] = params
58+ (out / "manifest.json").write_text(json.dumps(manifest, indent=1))
59+ print(f"site written to {out}")
60+
61+
62+if __name__ == "__main__":
63+ main()
web/index.htmladded+54−0View file
@@ -0,0 +1,54 @@
1+<!doctype html>
2+<html lang="en">
3+<head>
4+<meta charset="utf-8">
5+<meta name="viewport" content="width=device-width, initial-scale=1">
6+<title>Cymbal solutions</title>
7+<style>
8+:root { --fg: #1b1b1b; --bg: #fff; --muted: #666; --line: #ddd; --accent: #2a6db0; }
9+@media (prefers-color-scheme: dark) { :root { --fg: #e6e6e6; --bg: #151515; --muted: #999; --line: #333; --accent: #7ab3e8; } }
10+body { font-family: system-ui, sans-serif; color: var(--fg); background: var(--bg); margin: 0; padding: 16px; line-height: 1.45; }
11+main { max-width: 1100px; margin: 0 auto; }
12+a { color: var(--accent); }
13+table { border-collapse: collapse; width: 100%; font-size: 0.92em; }
14+th, td { border-bottom: 1px solid var(--line); padding: 6px 8px; text-align: left; vertical-align: top; }
15+td.num { text-align: right; font-variant-numeric: tabular-nums; }
16+.muted { color: var(--muted); }
17+.wrap { overflow-x: auto; }
18+#view { width: 100%; height: 480px; border: 1px solid var(--line); border-radius: 4px; touch-action: none; }
19+.controls { display: flex; flex-wrap: wrap; gap: 12px; align-items: center; margin: 8px 0; }
20+.controls label { display: flex; gap: 6px; align-items: center; }
21+input[type=range] { width: 200px; }
22+button { font: inherit; padding: 4px 10px; }
23+</style>
24+<script type="importmap">
25+{ "imports": {
26+ "three": "https://cdn.jsdelivr.net/npm/three@0.160.0/build/three.module.js",
27+ "three/addons/": "https://cdn.jsdelivr.net/npm/three@0.160.0/examples/jsm/"
28+} }
29+</script>
30+</head>
31+<body><main>
32+<h1>Solutions to the cymbal challenges</h1>
33+<p>This is the site of the solution repository <a href="https://vault1.magland.org/challenges/cymbal-solutions">challenges/cymbal-solutions</a>. The comparison of all registered solutions against the reference is on the <a href="https://vault1.magland.org/challenges/cymbal/site/">challenge repository's results site</a>; this page lists this repository's own submissions and animates the exported displacement fields as a 3D surface.</p>
34+
35+<h2>Displacement animation</h2>
36+<div class="controls">
37+ <label>Solution <select id="fieldSelect"></select></label>
38+ <button id="playBtn">Play</button>
39+ <label>Speed <input id="speed" type="range" min="-3" max="0" step="0.1" value="-2"> <span id="speedLabel"></span></label>
40+ <label>Exaggeration <input id="exag" type="range" min="0" max="3" step="0.1" value="1.5"> <span id="exagLabel"></span></label>
41+ <label><input id="derigid" type="checkbox" checked> remove rigid motion</label>
42+ <label>Time <input id="timeSlider" type="range" min="0" max="1" step="1" value="0"> <span id="timeLabel" class="muted"></span></label>
43+</div>
44+<div id="view"></div>
45+<p id="fieldStatus" class="muted">Loading…</p>
46+
47+<h2>Submissions</h2>
48+<div id="playback"></div>
49+<div id="submissions"></div>
50+</main>
51+<script src="player.js"></script>
52+<script type="module" src="viewer.js"></script>
53+</body>
54+</html>
web/viewer.jsadded+272−0View file
@@ -0,0 +1,272 @@
1+// Lists this repository's submissions and animates the exported HDF5 displacement fields of the disk as a 3D surface.
2+import * as THREE from "three";
3+import { OrbitControls } from "three/addons/controls/OrbitControls.js";
4+import h5wasm from "https://cdn.jsdelivr.net/npm/h5wasm@0.7/dist/esm/hdf5_hl.js";
5+
6+const manifest = await (await fetch("manifest.json")).json();
7+const esc = (s) => String(s).replace(/[&<>"']/g, (c) => ({ "&": "&", "<": "<", ">": ">", '"': """, "'": "'" }[c]));
8+
9+// ---------------------------------------------------------------- submissions table
10+const subsDiv = document.getElementById("submissions");
11+const fieldSolutions = [];
12+for (const [cid, ch] of Object.entries(manifest.challenges)) {
13+ const h3 = document.createElement("h3");
14+ h3.textContent = cid;
15+ subsDiv.appendChild(h3);
16+ for (const sub of ch.submissions) {
17+ const p = document.createElement("p");
18+ p.innerHTML = `<b>${esc(sub.submission_id)}</b> by ${esc(sub.submitter.name)}${sub.date ? " (" + esc(sub.date) + ")" : ""}. ${esc(sub.description || "")}` +
19+ (sub.code_url ? ` <a href="${esc(sub.code_url)}">code</a>` : "");
20+ subsDiv.appendChild(p);
21+ const rows = sub.solutions.map((s) => {
22+ if (s.fields && s.fields.available !== false) fieldSolutions.push({ label: `${sub.submission_id}/${s.solution_id}`, path: s.fields.path, challenge: cid });
23+ const fieldsLink = s.fields ? (s.fields.available === false ? ' <span class="muted">(fields stored in LFS, not published here)</span>' : ` <a href="${esc(s.fields.path)}">fields</a>`) : "";
24+ return `<tr><td>${esc(s.solution_id)}</td><td>${esc(s.method || "")}<br><span class="muted"><code>${esc(JSON.stringify(s.parameters || {}))}</code></span></td>` +
25+ `<td class="num">${s.timing_seconds ?? ""}</td><td class="muted">${esc((s.more_accurate_than || []).join(", "))}</td>` +
26+ `<td><div class="cplayer" data-src="${esc(s.output)}"></div><br><a href="${esc(s.output)}">wav</a>${fieldsLink}</td></tr>`;
27+ }).join("");
28+ const wrap = document.createElement("div");
29+ wrap.className = "wrap";
30+ wrap.innerHTML = `<table><tr><th>Solution</th><th>Method</th><th>Time (s)</th><th>Claimed more accurate than</th><th>Output</th></tr>${rows}</table>`;
31+ subsDiv.appendChild(wrap);
32+ }
33+}
34+
35+const gains = manifest.playback_gain_db || {};
36+const firstGain = Object.values(gains).find((g) => typeof g === "number");
37+ChallengePlayer.controls(document.getElementById("playback"), { gainDb: typeof firstGain === "number" ? firstGain : 20, limiter: true });
38+ChallengePlayer.mount();
39+
40+// ---------------------------------------------------------------- 3D viewer
41+const view = document.getElementById("view");
42+const status = document.getElementById("fieldStatus");
43+const select = document.getElementById("fieldSelect");
44+const playBtn = document.getElementById("playBtn");
45+const speed = document.getElementById("speed");
46+const speedLabel = document.getElementById("speedLabel");
47+const exag = document.getElementById("exag");
48+const exagLabel = document.getElementById("exagLabel");
49+const timeSlider = document.getElementById("timeSlider");
50+const derigid = document.getElementById("derigid");
51+const timeLabel = document.getElementById("timeLabel");
52+
53+const renderer = new THREE.WebGLRenderer({ antialias: true });
54+renderer.setPixelRatio(window.devicePixelRatio || 1);
55+view.appendChild(renderer.domElement);
56+const scene = new THREE.Scene();
57+const dark = matchMedia("(prefers-color-scheme: dark)").matches;
58+scene.background = new THREE.Color(dark ? 0x151515 : 0xffffff);
59+const camera = new THREE.PerspectiveCamera(40, 2, 0.01, 100);
60+camera.position.set(0.45, -0.6, 0.4);
61+camera.up.set(0, 0, 1);
62+const controls = new OrbitControls(camera, renderer.domElement);
63+controls.target.set(0, 0, 0);
64+scene.add(new THREE.AmbientLight(0xffffff, 0.6));
65+const dl = new THREE.DirectionalLight(0xffffff, 0.9);
66+dl.position.set(0.5, -0.8, 1.0);
67+scene.add(dl);
68+
69+let field = null; // {x, y, t, u, nx, ny, nt, umax, emax, plane, vidx, verts, params}
70+let mesh = null;
71+let frame = 0;
72+let playing = false;
73+let lastTs = null;
74+let simTime = 0;
75+
76+function resize() {
77+ const w = view.clientWidth, h = view.clientHeight;
78+ renderer.setSize(w, h, false);
79+ camera.aspect = w / h;
80+ camera.updateProjectionMatrix();
81+}
82+window.addEventListener("resize", () => { resize(); render(); });
83+resize();
84+
85+function buildMesh() {
86+ if (mesh) { scene.remove(mesh); mesh.geometry.dispose(); mesh.material.dispose(); }
87+ const { x, y, nx, ny, vidx, verts, params } = field;
88+ // vertices: the grid points inside the disk; triangles: cells whose four corners are all inside
89+ const pos = new Float32Array(verts.length * 3);
90+ for (let k = 0; k < verts.length; k++) {
91+ const [i, j] = verts[k];
92+ pos[3 * k] = x[i]; pos[3 * k + 1] = y[j]; pos[3 * k + 2] = 0;
93+ }
94+ const idx = [];
95+ for (let j = 0; j < ny - 1; j++) for (let i = 0; i < nx - 1; i++) {
96+ const a = vidx[j * nx + i], b = vidx[j * nx + i + 1], c = vidx[(j + 1) * nx + i], d = vidx[(j + 1) * nx + i + 1];
97+ if (a < 0 || b < 0 || c < 0 || d < 0) continue;
98+ idx.push(a, b, d, a, d, c);
99+ }
100+ const geo = new THREE.BufferGeometry();
101+ geo.setAttribute("position", new THREE.BufferAttribute(pos, 3));
102+ geo.setAttribute("color", new THREE.BufferAttribute(new Float32Array(verts.length * 3), 3));
103+ geo.setIndex(idx);
104+ const mat = new THREE.MeshStandardMaterial({ vertexColors: true, side: THREE.DoubleSide, metalness: 0.1, roughness: 0.7 });
105+ mesh = new THREE.Mesh(geo, mat);
106+ scene.add(mesh);
107+ // the rim at rest, for reference
108+ const R = params.disk_radius || Math.max(...x);
109+ const circ = [];
110+ for (let k = 0; k <= 128; k++) circ.push(new THREE.Vector3(R * Math.cos(2 * Math.PI * k / 128), R * Math.sin(2 * Math.PI * k / 128), 0));
111+ mesh.add(new THREE.Line(new THREE.BufferGeometry().setFromPoints(circ), new THREE.LineBasicMaterial({ color: dark ? 0x555555 : 0xbbbbbb })));
112+ // strike marker
113+ if (typeof params.strike_x === "number") {
114+ const marker = new THREE.Mesh(new THREE.SphereGeometry(0.004, 12, 12), new THREE.MeshBasicMaterial({ color: 0xd04040 }));
115+ marker.position.set(params.strike_x, params.strike_y, 0.004);
116+ mesh.add(marker);
117+ }
118+}
119+
120+function updateSurface() {
121+ if (!field || !mesh) return;
122+ const { u, nx, ny, x, y, plane, verts } = field;
123+ const rigidOff = derigid.checked;
124+ const umax = rigidOff ? field.emax : field.umax;
125+ const pos = mesh.geometry.attributes.position;
126+ const col = mesh.geometry.attributes.color;
127+ const scale = umax > 0 ? (0.03 * Math.pow(10, parseFloat(exag.value)) / umax) * 0.1 : 0; // exaggerated height
128+ const base = frame * nx * ny;
129+ const c0 = plane[3 * frame], gx = plane[3 * frame + 1], gy = plane[3 * frame + 2];
130+ for (let k = 0; k < verts.length; k++) {
131+ const [i, j] = verts[k];
132+ let v = u[base + j * nx + i];
133+ if (rigidOff) v -= c0 + gx * x[i] + gy * y[j];
134+ pos.setZ(k, v * scale);
135+ const s = umax > 0 ? Math.max(-1, Math.min(1, v / umax)) : 0;
136+ col.setXYZ(k, s > 0 ? 1 : 1 + s * 0.8, 1 - Math.abs(s) * 0.7, s < 0 ? 1 : 1 - s * 0.8);
137+ }
138+ pos.needsUpdate = true;
139+ col.needsUpdate = true;
140+ mesh.geometry.computeVertexNormals();
141+ timeSlider.value = String(frame);
142+ timeLabel.textContent = `${(field.t[frame] * 1e3).toFixed(0)} ms`;
143+}
144+
145+function render() { renderer.render(scene, camera); }
146+controls.addEventListener("change", render);
147+
148+function solve3(A, b) {
149+ // Gaussian elimination for the 3 x 3 normal equations of the plane fit
150+ const M = A.map((row, i) => [...row, b[i]]);
151+ for (let c = 0; c < 3; c++) {
152+ let p = c;
153+ for (let r = c + 1; r < 3; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
154+ [M[c], M[p]] = [M[p], M[c]];
155+ if (Math.abs(M[c][c]) < 1e-300) return [0, 0, 0];
156+ for (let r = 0; r < 3; r++) if (r !== c) {
157+ const f = M[r][c] / M[c][c];
158+ for (let k = c; k < 4; k++) M[r][k] -= f * M[c][k];
159+ }
160+ }
161+ return [M[0][3] / M[0][0], M[1][3] / M[1][1], M[2][3] / M[2][2]];
162+}
163+
164+async function loadField(entry) {
165+ const path = entry.path;
166+ status.textContent = `Loading ${path} (this may take a while, the file is large)…`;
167+ playing = false;
168+ playBtn.textContent = "Play";
169+ try {
170+ const { FS } = await h5wasm.ready;
171+ const resp = await fetch(path);
172+ const buf = new Uint8Array(await resp.arrayBuffer());
173+ if (buf.length < 200 && new TextDecoder().decode(buf).startsWith("version https://git-lfs")) {
174+ status.textContent = `${path} is a git LFS pointer; the object was not published with this site.`;
175+ return;
176+ }
177+ const name = "current.h5";
178+ try { FS.unlink(name); } catch (_) {}
179+ FS.writeFile(name, buf);
180+ const file = new h5wasm.File(name, "r");
181+ const x = Float64Array.from(file.get("x").value);
182+ const y = Float64Array.from(file.get("y").value);
183+ const t = Float64Array.from(file.get("t").value);
184+ const u = Float32Array.from(file.get("u").value);
185+ file.close();
186+ const nx = x.length, ny = y.length, nt = t.length;
187+ // grid points inside the disk are those with finite values in the first frame
188+ const vidx = new Int32Array(nx * ny).fill(-1);
189+ const verts = [];
190+ for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
191+ if (Number.isFinite(u[j * nx + i])) { vidx[j * nx + i] = verts.length; verts.push([i, j]); }
192+ }
193+ // best-fit plane per frame over the inside points (rigid translation and rotations of the free disk)
194+ let s1 = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0;
195+ for (const [i, j] of verts) { s1 += 1; sx += x[i]; sy += y[j]; sxx += x[i] * x[i]; sxy += x[i] * y[j]; syy += y[j] * y[j]; }
196+ const A = [[s1, sx, sy], [sx, sxx, sxy], [sy, sxy, syy]];
197+ const plane = new Float64Array(nt * 3);
198+ let umax = 0, emax = 0;
199+ for (let k = 0; k < nt; k++) {
200+ const base = k * nx * ny;
201+ let b0 = 0, bx = 0, by = 0;
202+ for (const [i, j] of verts) { const v = u[base + j * nx + i]; b0 += v; bx += v * x[i]; by += v * y[j]; }
203+ const [c0, gx, gy] = solve3(A, [b0, bx, by]);
204+ plane[3 * k] = c0; plane[3 * k + 1] = gx; plane[3 * k + 2] = gy;
205+ for (const [i, j] of verts) {
206+ const v = u[base + j * nx + i];
207+ const e = v - c0 - gx * x[i] - gy * y[j];
208+ if (Math.abs(v) > umax) umax = Math.abs(v);
209+ if (Math.abs(e) > emax) emax = Math.abs(e);
210+ }
211+ }
212+ const params = (manifest.parameters || {})[entry.challenge] || {};
213+ field = { x, y, t, u, nx, ny, nt, umax, emax, plane, vidx, verts, params };
214+ timeSlider.max = String(field.nt - 1);
215+ frame = 0;
216+ simTime = 0;
217+ buildMesh();
218+ updateSurface();
219+ render();
220+ status.textContent = `${path}: ${verts.length} grid points inside the disk × ${field.nt} frames, max |u| = ${(umax * 1e3).toFixed(2)} mm including rigid motion, ${(emax * 1e3).toFixed(3)} mm elastic`;
221+ } catch (e) {
222+ status.textContent = `Could not load ${path}: ${e}`;
223+ }
224+}
225+
226+function tick(ts) {
227+ if (!playing || !field) return;
228+ if (lastTs !== null) {
229+ const factor = Math.pow(10, parseFloat(speed.value));
230+ simTime += ((ts - lastTs) / 1000) * factor;
231+ const tEnd = field.t[field.nt - 1];
232+ if (simTime > tEnd) simTime = 0;
233+ const dt = field.t[1] - field.t[0];
234+ frame = Math.min(field.nt - 1, Math.round(simTime / dt));
235+ updateSurface();
236+ render();
237+ }
238+ lastTs = ts;
239+ requestAnimationFrame(tick);
240+}
241+
242+playBtn.addEventListener("click", () => {
243+ playing = !playing;
244+ playBtn.textContent = playing ? "Pause" : "Play";
245+ lastTs = null;
246+ if (playing) requestAnimationFrame(tick);
247+});
248+timeSlider.addEventListener("input", () => {
249+ if (!field) return;
250+ frame = parseInt(timeSlider.value, 10);
251+ simTime = field.t[frame];
252+ updateSurface();
253+ render();
254+});
255+const updSpeed = () => { speedLabel.textContent = `${Math.pow(10, parseFloat(speed.value)).toPrecision(2)} × real time`; };
256+const updExag = () => { exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`; updateSurface(); render(); };
257+speed.addEventListener("input", updSpeed);
258+exag.addEventListener("input", updExag);
259+derigid.addEventListener("change", () => { updateSurface(); render(); });
260+updSpeed();
261+exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`;
262+select.addEventListener("change", () => loadField(fieldSolutions[parseInt(select.value, 10)]));
263+
264+fieldSolutions.forEach((f, k) => {
265+ const o = document.createElement("option");
266+ o.value = String(k);
267+ o.textContent = f.label;
268+ select.appendChild(o);
269+});
270+render();
271+if (fieldSolutions.length) loadField(fieldSolutions[0]);
272+else status.textContent = "No solution in this repository exports a displacement field that is published with this site.";