Laplacian loss term, rigid constraint by default, geometry-keyed mode cache, README on both solvers
11 changed files+50−21
.gitignoremodified+1−0View file
@@ -3,3 +3,4 @@ __pycache__/
33 *.progress.json
44 *.meta.json
55 site/
6+cache/
README.mdmodified+4−2View file
@@ -12,12 +12,14 @@ solvers/ source code of the submissions
1212 web/ this repository's site (built by web/build.py)
1313 ```
1414
15-## The submission hermite-fem
15+## The submissions hermite-fem and hermite-modal
1616
17-`solvers/hermite-fem/fem_shell.py` solves the shallow-shell model of MODEL.md directly in the displacements, with no modal reduction. The annulus is approximated by the square cells of a uniform Cartesian grid whose centers lie inside it (a staircase rim and hole). The three fields u, v, w are expanded in bicubic Hermite (Bogner-Fox-Schmit) elements, C1 across cells, four unknowns per node and field. The weak form is integrated with 3 x 3 Gauss points per cell exactly as written: Marguerre membrane strains, bending strains, stress resultants with the Kelvin-Voigt rates, and the assembled virtual work as the internal force, which a test confirms to be the exact gradient of the discrete strain energy; free edges are the natural boundary conditions and need no treatment. The mass is lumped by diagonal scaling of the consistent element mass. The equations are stepped with the explicit leapfrog scheme at 85 percent of the stability limit, which the in-plane block sets; the lumped masses of the in-plane derivative unknowns are scaled by 16, which halves the step count (14 substeps per output sample at 8 mm, about 28 at 4 mm) and moves the in-plane resonances by 0.15 percent while leaving the bending modes unchanged. The rigid constraint of MODEL.md is enforced by applying each stick's reaction as the rigid inertial load that cancels its net force and moments, and by projecting the six mass-orthonormal rigid fields out of the state at every chunk of steps. The sticks are point masses stepped with the same scheme through the Hertz law, each started at the current surface height under it with relative speed v_0. The Rayleigh integral is accumulated at every step from the transverse acceleration at the Gauss points, spread onto the output samples at the retarded times with cubic Lagrange interpolation, so no acceleration history is stored. The step runs on a GPU through PyTorch, the constitutive and output blocks fused by torch.compile and the whole step captured as a CUDA graph: about 5100 steps per second at 8 mm (25488 unknowns, a quarter of an hour for the seven seconds).
17+`solvers/hermite-fem/fem_shell.py` solves the shallow-shell model of MODEL.md directly in the displacements, with no modal reduction. The annulus is approximated by the square cells of a uniform Cartesian grid whose centers lie inside it (a staircase rim and hole). The three fields u, v, w are expanded in bicubic Hermite (Bogner-Fox-Schmit) elements, C1 across cells, four unknowns per node and field. The weak form is integrated with 3 x 3 Gauss points per cell exactly as written: Marguerre membrane strains, bending strains, stress resultants with the Kelvin-Voigt rates, and the assembled virtual work as the internal force, which a test confirms to be the exact gradient of the discrete strain energy; free edges are the natural boundary conditions and need no treatment. The mass is lumped by diagonal scaling of the consistent element mass. The equations are stepped with the explicit leapfrog scheme at 85 percent of the stability limit, which the in-plane block sets through the tiny lumped masses of its derivative unknowns (33 substeps per output sample at 8 mm on the flat annulus, 40 on the crash cymbal, about twice that at 5 mm). Scaling those masses up would halve the step count, and an option does it, but it pulls the spurious derivative-dominated in-plane modes from 40 to 60 kHz down into the audible band, where the cascade pumps them into steady ringing lines; it is off by default. The rigid constraint of MODEL.md is enforced by applying each stick's reaction as the rigid inertial load that cancels its net force and moments, and by projecting the six mass-orthonormal rigid fields out of the state at every chunk of steps. The three loss terms of MODEL.md are all evaluated at the Gauss points: the Kelvin-Voigt rates inside the stress resultants and the Laplacian term as the virtual work of 2 sigma1 rho h grad(w_t). The sticks are point masses stepped with the same scheme through the Hertz law, each started at the current surface height under it with relative speed v_0. The Rayleigh integral is accumulated at every step from the transverse acceleration at the Gauss points, spread onto the output samples at the retarded times with cubic Lagrange interpolation, so no acceleration history is stored. The step runs on a GPU through PyTorch, the constitutive and output blocks fused by torch.compile and the whole step captured as a CUDA graph: about 5100 steps per second at 8 mm (25488 unknowns, a quarter of an hour for the seven seconds).
1818
1919 Checks: the linearized assembly reproduces the free annulus's lowest frequencies to within 1 percent at 8 mm (21.7, 37.6, 50.6, 85.5 Hz against 21.9, 38.0, 50.9, 85.9), the lumped mass shifts them by under 0.3 percent relative to the consistent mass, the internal force matches the energy gradient to 1e-9, and the retarded accumulation matches a direct integral to 4e-5. On the crash geometry the two tilt modes, exact zero modes of the continuum, come out at 2.95 Hz at 8 mm because the bicubic space cannot represent the in-plane displacement that accompanies a rigid tilt of the cosine bell; the constraint removes them anyway. A run of the flat challenge with the shell left free confirmed what MODEL.md says about the constraint: the shell tilted by about a radian in half a second, the partials glided down, the treble drained away, and the explicit run diverged at 0.66 s.
2020
21+`solvers/hermite-fem/modal_shell.py` solves the same discrete system in its eigenmodes below a cutoff (20 kHz by default): the linear part is propagated exactly by the two-step scheme of the v1 modal solver at a few substeps per sample, the nonlinear force and the Laplacian loss are evaluated at the Gauss points by the same routine as the finite element solver and projected on the modes, and the six rigid fields are excluded from the eigenproblem by a penalty (Woodbury-corrected shift-invert slices). The modes are cached under `solvers/hermite-fem/cache/`, keyed by geometry. The two solvers therefore differ only in the modal truncation and the time integration, and their difference measures the effect of the modes above the cutoff, which the finite element solver keeps and the modal one drops.
22+
2123 ```
2224 bash solvers/hermite-fem/run_chain.sh <cymbal2-repo> # on a CUDA GPU: 8 mm then 4 mm for both challenges, hours
2325 python solvers/hermite-fem/check_modes.py <cymbal2-repo>/challenges/cymbal2-flat-v1/challenge.json 0.008 # linear check, CPU
solutions/cymbal2-crash-v1/hermite-fem/fem-g8.h5added+3−0View file
@@ -0,0 +1,3 @@
1+version https://git-lfs.github.com/spec/v1
2+oid sha256:cbbe994a7ee33b34041065d396bd921c9ed3f3f941ea7bf69694b0c782feb952
3+size 36321414
solutions/cymbal2-crash-v1/hermite-fem/fem-g8.kept-partial.wavadded+0−0View file
Binary file not shown.
solutions/cymbal2-crash-v1/hermite-fem/fem-g8.wavadded+0−0View file
Binary file not shown.
solutions/cymbal2-crash-v1/hermite-modal/modal-g8-f20k-m4.h5added+3−0View file
@@ -0,0 +1,3 @@
1+version https://git-lfs.github.com/spec/v1
2+oid sha256:48446e184cad11e5fb9ece5a33a49a028916fe607531048972b5d60316f67768
3+size 39479365
solutions/cymbal2-crash-v1/hermite-modal/modal-g8-f20k-m4.wavadded+0−0View file
Binary file not shown.
solvers/hermite-fem/fem_shell.pymodified+13−6View file
@@ -192,9 +192,11 @@ class ShellFEM:
192192 trans = me[:, 0::4].sum(axis=1) # translational entries
193193 me = me * (mcell / trans)[:, None]
194194 np.add.at(Mdiag, self.cell_dofs[:, 16 * blk : 16 * (blk + 1)], me)
195- # rotational mass scaling of the in-plane derivative unknowns: their lumped masses are tiny and set the explicit
196- # time step through the in-plane block; scaling them by inplane_scale (16 by default) roughly halves the step count
197- # and moves the in-plane resonances by about 0.15 percent, the bending modes not at all.
195+ # optional scaling of the lumped masses of the in-plane derivative unknowns (--inplane-scale). They are tiny and
196+ # set the explicit time step through the in-plane block, and scaling them by 16 halves the step count; but it
197+ # also pulls the spurious derivative-dominated in-plane modes down from 40 to 60 kHz into the audible band, where
198+ # the nonlinear cascade pumps them and they ring as steady lines (seen at 10.7 kHz on the 8 mm grid). The default
199+ # is therefore 1 (no scaling).
198200 nw = 4 * self.n_nodes
199201 for b in (1, 2):
200202 for d in (1, 2, 3):
@@ -301,6 +303,7 @@ class ShellFEM:
301303 Bw = T64(self.B * np.repeat(self.wq, 9)[:, None]) # (81, 48) with quadrature weights folded in
302304 Zx, Zy, Cm, Db = T64(self.Zx), T64(self.Zy), T64(self.Cm), T64(self.Db)
303305 nu, eta = p["nu"], p["eta"]
306+ sig1rh = T64(2.0 * p.get("sigma1", 0.0) * p["rho"] * self.hq) # 2 sigma1 rho h at the Gauss points (Laplacian damping of w)
304307 lin = 0.0 if self.linear else 1.0
305308 n_dof = self.n_dof
306309
@@ -326,7 +329,9 @@ class ShellFEM:
326329 Mxx = Db * (kxx + nu * kyy)
327330 Myy = Db * (kyy + nu * kxx)
328331 Mxy = Db * (1 - nu) * kxy
329- return torch.stack([Nxx * ax + Nxy * ay, Nyy * ay + Nxy * ax, Mxx, Myy, 2 * Mxy, Nxx, Nxy, Nxy, Nyy], dim=-1) # (C, G, 9)
332+ # generalized stresses conjugate to [w_x, w_y, w_xx, w_yy, w_xy, u_x, u_y, v_x, v_y]; the sigma1 term is the
333+ # virtual work 2 sigma1 rho h grad(w_t) . grad(delta w) of the Laplacian loss
334+ return torch.stack([Nxx * ax + Nxy * ay + sig1rh * Dv[..., 0], Nyy * ay + Nxy * ax + sig1rh * Dv[..., 1], Mxx, Myy, 2 * Mxy, Nxx, Nxy, Nxy, Nyy], dim=-1) # (C, G, 9)
330335
331336 constitutive_c = torch.compile(constitutive, dynamic=False) if (compile_ and str(dev).startswith("cuda")) else constitutive
332337
@@ -625,8 +630,10 @@ def main():
625630 tag = take("--tag", str)
626631 c_air = take("--c-air", float)
627632 check = "--check" in args
628- no_drift = "--no-drift" in args
629- inplane_scale = take("--inplane-scale", float, 16.0)
633+ no_drift = "--free" not in args # the rigid constraint of MODEL.md is the default; --free leaves the shell free (diagnostics only)
634+ if "--no-drift" in args:
635+ args.remove("--no-drift")
636+ inplane_scale = take("--inplane-scale", float, 1.0)
630637 compile_ = "--no-compile" not in args
631638 linear = "--linear" in args
632639 fields = "--fields" in args
solvers/hermite-fem/fetch_progress.shmodified+7−6View file
@@ -3,13 +3,14 @@
33 # and print a one-line status per run. Usage: bash fetch_progress.sh
44 here="$(cd "$(dirname "$0")/../.." && pwd)"
55 remote=office:/home/magland/src/vault1/challenges/cymbal2-solutions
6-for c in cymbal2-flat-v1 cymbal2-crash-v1; do
7- mkdir -p "$here/solutions/$c/hermite-fem"
8- scp -q "$remote/solutions/$c/hermite-fem/*.progress.json" "$remote/solutions/$c/hermite-fem/*.partial.wav" "$here/solutions/$c/hermite-fem/" 2>/dev/null
9- for j in "$here"/solutions/$c/hermite-fem/*.progress.json; do
6+for c in cymbal2-flat-v1 cymbal2-crash-v1; do for sub in hermite-fem hermite-modal; do
7+ mkdir -p "$here/solutions/$c/$sub"
8+ scp -q "$remote/solutions/$c/$sub/*.progress.json" "$remote/solutions/$c/$sub/*.partial.wav" "$here/solutions/$c/$sub/" 2>/dev/null
9+ for j in "$here"/solutions/$c/$sub/*.progress.json; do
1010 [ -f "$j" ] && python3 -c "
1111 import json,sys; d=json.load(open('$j'))
12-print(f\"$c {d['tag']}: t = {d['t']:.2f} / {d['duration']:.0f} s, {d['elapsed_s']/60:.0f} min elapsed, {d['eta_s']/60:.0f} min left, peak {d['peak_pa']:.1f} Pa, max|w| {d['max_abs_w_m']*1e3:.1f} mm, sticks: impulse {d['impulse_Ns']} N s, contact {d['contact_ms']} ms, peak force {d['peak_force_N']} N, rebound {d['rebound_mps']} m/s\")"
12+amp = f\"max|w| {d['max_abs_w_m']*1e3:.1f} mm\" if 'max_abs_w_m' in d else f\"{d.get('n_modes','?')} modes\"
13+print(f\"$c {d['tag']}: t = {d['t']:.2f} / {d['duration']:.0f} s, {d['elapsed_s']/60:.0f} min elapsed, {d['eta_s']/60:.0f} min left, peak {d['peak_pa']:.1f} Pa, {amp}, sticks: impulse {d['impulse_Ns']} N s, contact {d['contact_ms']} ms, peak force {d['peak_force_N']} N, rebound {d['rebound_mps']} m/s\")"
1314 done
14-done
15+done; done
1516 exit 0
solvers/hermite-fem/modal_shell.pymodified+18−6View file
@@ -103,7 +103,19 @@ class ModalShell:
103103 K = fem.assemble_linear()
104104 self.log(f" linear stiffness assembled ({K.nnz} nonzeros) in {time.perf_counter() - t0:.1f}s")
105105 t0 = time.perf_counter()
106- self.omega, self.Phi = sliced_modes(K, fem.Mdiag, fem.rigid, f_max, log=log)
106+ cache_dir = Path(__file__).resolve().parent / "cache"
107+ cache_dir.mkdir(exist_ok=True)
108+ import hashlib
109+
110+ geom_key = hashlib.sha1(json.dumps({k: v for k, v in fem.ch["parameters"].items() if k in ("disk_radius", "hole_radius", "bow_height", "bell_radius", "bell_height", "thickness", "thickness_center", "thickness_rim", "E", "nu", "rho")}, sort_keys=True).encode()).hexdigest()[:8]
111+ cache = cache_dir / f"modes-{fem.ch['challenge_id']}-g{fem.a * 1e3:g}-f{f_max / 1e3:g}k-s{fem.inplane_scale:g}-{geom_key}.npz"
112+ if cache.exists():
113+ z = np.load(cache)
114+ self.omega, self.Phi = z["omega"], z["Phi"]
115+ self.log(f" modes loaded from {cache.name}")
116+ else:
117+ self.omega, self.Phi = sliced_modes(K, fem.Mdiag, fem.rigid, f_max, log=log)
118+ np.savez(cache, omega=self.omega, Phi=self.Phi)
107119 self.n_modes = len(self.omega)
108120 self.log(f" {self.n_modes} modes in {time.perf_counter() - t0:.1f}s: f = {np.round(self.omega[:5] / 2 / np.pi, 2)} ... {self.omega[-1] / 2 / np.pi:.0f} Hz")
109121 p = fem.p
@@ -199,7 +211,7 @@ class ModalShell:
199211 Fm = -(PhiT @ f) + w2 * a + eta * w2 * adot - F @ gains
200212 a_next = a1 * a + a2 * a_prev + b * Fm
201213 add = Fm - w2 * a - al * ((a_next - a_prev) / dt)
202- ring[n % L] = add
214+ ring.index_copy_(0, (n % L).view(1), add.view(1, -1))
203215 zd.add_(dt * F / ms)
204216 z.add_(dt * zd)
205217 q_prev.copy_(q)
@@ -212,7 +224,7 @@ class ModalShell:
212224 # folds the table with the ring (delays counted back from that step), then the remaining steps of the sample
213225 step()
214226 idx = (nstep - 1 - ar) % L
215- out[k] = torch.sum(Gt * ring[idx])
227+ out.index_copy_(0, k.view(1), torch.sum(Gt * ring[idx]).view(1))
216228 for _ in range(m - 1):
217229 step()
218230
@@ -252,11 +264,11 @@ class ModalShell:
252264 if fld is not None:
253265 take_frame()
254266 for k in range(nsamp):
267+ k_t.fill_(k)
255268 if graph is not None:
256- k_t.fill_(k)
257269 graph.replay()
258270 else:
259- sample(k)
271+ sample(k_t)
260272 if fld is not None and (k + 1) % fld["every"] == 0 and len(fld["frames"]) < fld["nt"]:
261273 take_frame()
262274 now = time.perf_counter()
@@ -307,7 +319,7 @@ def main():
307319 report = take("--report", float, 60.0)
308320 speed = take("--speed", float)
309321 tag = take("--tag", str)
310- inplane_scale = take("--inplane-scale", float, 16.0)
322+ inplane_scale = take("--inplane-scale", float, 1.0)
311323 fields = "--fields" in args
312324 compile_ = "--no-compile" not in args
313325 use_graph = "--no-graph" not in args
solvers/hermite-fem/run_chain.shmodified+1−1View file
@@ -10,7 +10,7 @@ py="${PYTHON:-/home/magland/src/vault1/challenges/cymbal-env/bin/python}"
1010 run() { # challenge spacing_m
1111 out="$here/../../solutions/$1/hermite-fem"
1212 mkdir -p "$out"
13- "$py" "$here/fem_shell.py" "$spec/challenges/$1/challenge.json" "$out" --spacing "$2" --device cuda --fields --report 60 --no-drift
13+ "$py" "$here/fem_shell.py" "$spec/challenges/$1/challenge.json" "$out" --spacing "$2" --device cuda --fields --report 60
1414 }
1515 run cymbal2-crash-v1 0.008
1616 run cymbal2-crash-v1 0.005