Time step from the global stability bound (power iteration), with numerics.stability_bound to keep the element bound
7 changed files+72−9
README.mdmodified+1−1View file
@@ -9,7 +9,7 @@ pip install -e ".[dev]" # numpy, scipy, jax (CPU), pyyaml, h5py, pytest
99 pip install "jax[cuda12]" # optional: run the time loop on a GPU
1010 ```
1111
12-Runs at the default resolution (6 mm elements, about 4700 nodes and 40 sub-steps per audio sample) need a GPU to finish in minutes. On a CPU, use a coarser mesh and a shorter duration, as in `examples/quick_cpu.yaml`.
12+Runs at the default resolution (6 mm elements, about 4700 nodes and 38 sub-steps per audio sample) need a GPU to finish in minutes. On a CPU, use a coarser mesh and a shorter duration, as in `examples/quick_cpu.yaml`.
1313
1414 ## Usage
1515
docs/model.mdmodified+2−2View file
@@ -150,7 +150,7 @@ where $P$ is the rigid-mode projection of section 5.3. This is the symplectic Eu
150150
151151 ### 5.2 Stability and the choice of the time step
152152
153-For the linear undamped system the scheme is stable when $\Delta t \leq 2 / \omega_{max}$, with $\omega_{max}^2$ the largest eigenvalue of $M^{-1} K$. Instead of solving that global eigenproblem we use the element bound (Irons; see Hughes 2000): $\lambda_{max}(M^{-1} K) \leq \max_e \lambda_{max}(M_e^{-1} K_e)$, where $K_e$ is the $15 \times 15$ tangent stiffness of the element at rest and $M_e$ its share of the lumped mass. The bound is computed in float64 at build time and the limiting element (the thinnest, near the rim) is reported. The nonlinear membrane stiffening raises the instantaneous frequencies during large motion, so we take only a fraction `cfl_safety` (default 0.6) of the bound. The sub-step count per audio sample is $m = \lceil 1 / (f_s\, c_{safety}\, \Delta t_{bound}) \rceil$ and $\Delta t = 1 / (m f_s)$, so that the radiation buffers decimate exactly to the audio rate. The run is extended by $W/2 = 16$ audio samples beyond the requested duration so that the last samples of the decimated records (section 9) are final. If the energy diagnostic grows during a run the step is too large; the remedy is a smaller `cfl_safety`. The run aborts if the energy becomes non-finite.
153+For the linear undamped system the scheme is stable when $\Delta t \leq 2 / \omega_{max}$, with $\omega_{max}^2$ the largest eigenvalue of $M^{-1} K$. We estimate $\omega_{max}$ from the stiffness assembled from the $15 \times 15$ tangent stiffnesses $K_e$ of the elements at rest, by power iteration on $M^{-1/2} K M^{-1/2}$ (`global_omega_max` in `shell.py`). The start vector $x_i = \sin(i + 1)$ and the stopping rule (the Rayleigh quotient changes by at most $10^{-7}$ relative over 10 iterations) are fixed so that cymsimweb, which implements the same iteration, arrives at the same value; the Rayleigh quotient approaches $\lambda_{max}$ from below, and on the meshes we tried the stopping error in $\omega_{max}$ was below $10^{-5}$. At 6 mm elements the iteration stops after about 70 iterations, on coarse meshes after a few thousand, in under a second either way. If it does not converge in 20000 iterations, the element bound below is used with a warning. The element bound (Irons; see Hughes 2000), $\lambda_{max}(M^{-1} K) \leq \max_e \lambda_{max}(M_e^{-1} K_e)$ with $M_e$ the element's share of the lumped mass, is also computed and reported with its limiting element. It is safe without any iteration but overestimates $\omega_{max}$ by a factor of 1.3 to 1.5 on our meshes (1.45 at 6 mm elements, where the largest global mode is a rotation of about ten nodes near the bell), so a run with `numerics.stability_bound: element`, the only option before October 2026, takes about 1.45 times as many sub-steps. At the default 6 mm and 44.1 kHz the step is 38 sub-steps per sample with the global bound against 54 with the element bound. In a comparison with a float64 run of the mallet example (cymsimweb, first 50 ms), the float32 listener signal at 38 sub-steps was as close to the float64 one as at 54 (relative rms error near $10^{-2}$, octave levels within 0.4 dB), and strong strikes stayed stable down to 25 sub-steps; the leapfrog phase error at these steps is about $10^{-4}$ at 20 kHz, well below the round-off of a float32 run. The nonlinear membrane stiffening raises the instantaneous frequencies during large motion, so we take only a fraction `cfl_safety` (default 0.6) of the bound. The sub-step count per audio sample is $m = \lceil 1 / (f_s\, c_{safety}\, \Delta t_{bound}) \rceil$ and $\Delta t = 1 / (m f_s)$, so that the radiation buffers decimate exactly to the audio rate. The run is extended by $W/2 = 16$ audio samples beyond the requested duration so that the last samples of the decimated records (section 9) are final. If the energy diagnostic grows during a run the step is too large; the remedy is a smaller `cfl_safety`. The run aborts if the energy becomes non-finite.
154154
155155 ### 5.3 Rigid-mode projection
156156
@@ -261,7 +261,7 @@ The natural next steps are, roughly in order of expected benefit: a convergence
261261 | $G_x, G_y$ | CST gradient operators | `Mesh.Gx`, `.Gy` |
262262 | $M$, $I_e$, $f$ | lumped mass, rotational inertia, its scaling factor | `ShellModel.M`, `.I_e`, `Numerics.rot_inertia_factor` |
263263 | $r_k$ | rigid-body modes | `ShellModel.Rm` |
264-| $\omega_{max}$, $\Delta t$, $m$ | stability bound, sub-step, sub-steps per sample | `ShellModel.omega_max`, `Simulator.dt`, `.substeps` |
264+| $\omega_{max}$, $\Delta t$, $m$ | stability bound (global, element), sub-step, sub-steps per sample | `ShellModel.omega_global`, `.omega_max`, `Simulator.dt`, `.substeps` |
265265 | $k, \alpha, P$ | contact stiffness, exponent, springs per patch | `ContactPatch.k_total`, `.alpha`, `.n_points` |
266266 | $z_s, v_s$ | exciter face height and velocity | state entries `zs`, `vs` |
267267 | $p_{front}, p_{back}$ | Rayleigh front and diffracted back pressure | `SimulationResult.pressure_rayleigh`, `.pressure_unbaffled` |
src/cymsim/config.pymodified+3−0View file
@@ -133,6 +133,7 @@ class Numerics:
133133 fs: int = 44100 # audio sample rate [Hz]; the sub-step rate is an integer multiple
134134 duration: float = 2.0 # simulated time [s]
135135 cfl_safety: float = 0.6 # fraction of the linear stability limit actually used
136+ stability_bound: str = "global" # "global": lambda_max(M^-1 K) by power iteration; "element": the Irons element bound (larger, used before 2026-10)
136137 rot_inertia_factor: float = 1 / 64 # lumped rotational inertia I = m (h^2/12 + f A_e); see docs
137138 inplane_mass_scale: float = 1.0 # > 1 slows in-plane waves and relaxes the time step; see docs
138139 steps_per_block: int = 20000 # sub-steps per compiled block
@@ -171,6 +172,8 @@ class SimulationConfig:
171172 raise ValueError("at least one hit is required")
172173 if n.dtype not in ("float32", "float64"):
173174 raise ValueError(f"numerics.dtype must be 'float32' or 'float64', got {n.dtype!r}")
175+ if n.stability_bound not in ("global", "element"):
176+ raise ValueError(f"numerics.stability_bound must be 'global' or 'element', got {n.stability_bound!r}")
174177 if not (0 < n.cfl_safety <= 1):
175178 raise ValueError("numerics.cfl_safety must lie in (0, 1]")
176179 if self.mesh.element_size <= 0 or self.mesh.element_size > p.radius / 2:
src/cymsim/shell.pymodified+39−1View file
@@ -11,6 +11,7 @@ and "Time integration".
1111 """
1212 from __future__ import annotations
1313
14+import logging
1415 import math
1516 from dataclasses import dataclass
1617
@@ -25,6 +26,34 @@ from .mesh import Mesh
2526
2627 N_DOF = 5 # (u, v, w, theta_x, theta_y) per node
2728
29+log = logging.getLogger(__name__)
30+
31+
32+def global_omega_max(K_e: np.ndarray, M: np.ndarray, tri: np.ndarray, tol: float = 1e-7,
33+ max_iter: int = 20000) -> tuple[float, int, bool]:
34+ """sqrt(lambda_max(M^-1 K)) of the stiffness assembled from the element matrices ``K_e``
35+ (E, 15, 15), node-major DOFs, by power iteration on M^-1/2 K M^-1/2 (lumped M, (N, 5)).
36+
37+ The start vector x_i = sin(i + 1) and the stopping rule (the Rayleigh quotient changes by at
38+ most ``tol`` relative over 10 iterations) are fixed, so that cymsimweb, which implements the
39+ same iteration, arrives at the same value. The Rayleigh quotient approaches lambda_max from
40+ below. Returns (omega, iterations, converged)."""
41+ E, L = len(tri), M.size
42+ dof = (N_DOF * tri[:, :, None] + np.arange(N_DOF)).reshape(E, 15)
43+ isq = 1.0 / np.sqrt(M.ravel())
44+ x = np.sin(np.arange(L) + 1.0)
45+ x /= np.linalg.norm(x)
46+ hist = []
47+ for it in range(1, max_iter + 1):
48+ ye = np.einsum("eab,eb->ea", K_e, (x * isq)[dof])
49+ y = np.bincount(dof.ravel(), weights=ye.ravel(), minlength=L) * isq
50+ lam = float(x @ y)
51+ hist.append(lam)
52+ x = y / np.linalg.norm(y)
53+ if it >= 20 and it % 10 == 0 and abs(lam - hist[-11]) <= tol * lam:
54+ return math.sqrt(lam), it, True
55+ return math.sqrt(hist[-1]), max_iter, False
56+
2857
2958 @dataclass
3059 class ElementChecks:
@@ -128,7 +157,16 @@ class ShellModel:
128157 self._lam = np.linalg.eigvalsh(S) # (E, 15) ascending
129158 self.limiting_element = int(np.argmax(self._lam[:, -1]))
130159 self.omega_max = math.sqrt(self._lam[:, -1].max())
131- self.dt_stable = 2.0 / self.omega_max
160+ # ---- the global lambda_max(M^-1 K) of the assembled tangent, which the element bound
161+ # overestimates (by about 1.45 at 6 mm elements); see docs/model.md
162+ self.K_rest = np.asarray(Khess) # (E, 15, 15) tangent at rest
163+ self.omega_global, self.omega_global_iterations, converged = global_omega_max(self.K_rest, M, tri)
164+ if numerics.stability_bound == "global" and not converged:
165+ log.warning("the power iteration for the global stability bound did not converge in %d "
166+ "iterations; using the element bound", self.omega_global_iterations)
167+ use_global = numerics.stability_bound == "global" and converged
168+ self.omega_step = min(self.omega_global, self.omega_max) if use_global else self.omega_max
169+ self.dt_stable = 2.0 / self.omega_step
132170
133171 # ---- device arrays in simulation precision
134172 f = lambda a: jnp.asarray(a, dtype=self.dtype)
src/cymsim/simulate.pymodified+5−3View file
@@ -255,9 +255,11 @@ class Simulator:
255255 self.n_blocks = math.ceil((self.n_out + self.span // 2) * self.substeps / num.steps_per_block)
256256 self.n_steps = self.n_blocks * num.steps_per_block
257257 xc, yc = self.mesh.centroid[sh.limiting_element]
258- log.info("omega_max bound %.3e rad/s (limiting element at r = %.3f m, h = %.2f mm); "
259- "dt = %.3e s, %d sub-steps per sample, %d steps",
260- sh.omega_max, math.hypot(xc, yc), sh.h_e[sh.limiting_element] * 1e3,
258+ log.info("omega_max: element bound %.3e rad/s (limiting element at r = %.3f m, h = %.2f mm), "
259+ "global %.3e rad/s (%d power iterations); %s bound used: dt = %.3e s, "
260+ "%d sub-steps per sample, %d steps",
261+ sh.omega_max, math.hypot(xc, yc), sh.h_e[sh.limiting_element] * 1e3, sh.omega_global,
262+ sh.omega_global_iterations, "global" if sh.omega_step < sh.omega_max else "element",
261263 self.dt, self.substeps, self.n_steps)
262264 self.taps = decimation_taps(self.substeps)
263265
tests/test_shell.pymodified+20−1View file
@@ -75,4 +75,23 @@ def test_damping_is_dissipative(shell):
7575
7676 def test_stability_bound_is_finite_and_sane(shell):
7777 assert 1e4 < shell.omega_max < 1e8
78- assert shell.dt_stable == 2.0 / shell.omega_max
78+ assert shell.omega_global < shell.omega_max # the element bound is an upper bound
79+ assert shell.dt_stable == 2.0 / shell.omega_global # the default numerics.stability_bound is "global"
80+
81+
82+def test_global_bound_matches_the_dense_eigenproblem(shell):
83+ """The power iteration finds the largest eigenvalue of M^-1 K, K assembled from the tangents at rest."""
84+ N, L = shell.mesh.n_nodes, shell.M.size
85+ K = np.zeros((L, L))
86+ dof = (5 * shell.mesh.tri[:, :, None] + np.arange(5)).reshape(-1, 15)
87+ for e in range(len(dof)):
88+ K[np.ix_(dof[e], dof[e])] += shell.K_rest[e]
89+ isq = 1 / np.sqrt(shell.M.ravel())
90+ lam = np.linalg.eigvalsh(isq[:, None] * K * isq[None, :])[-1]
91+ assert abs(shell.omega_global / np.sqrt(lam) - 1) < 1e-6
92+
93+
94+def test_element_bound_on_request(shell):
95+ sh = ShellModel(shell.mesh, CymbalProfile(Profile(), Imperfections(n_dents=30)), Material(),
96+ Numerics(dtype="float64", stability_bound="element"))
97+ assert sh.dt_stable == 2.0 / sh.omega_max
tools/make_web_fixture.pymodified+2−1View file
@@ -54,7 +54,8 @@ def main(path: str) -> None:
5454 lin_cfg = SimulationConfig.from_dict(cfg.to_dict())
5555 lin_cfg.numerics.nonlinear = False
5656 run_linear, _ = run_entry(lin_cfg)
57- fx["timestep"] = {"dt": sim.dt, "substeps": sim.substeps}
57+ fx["timestep"] = {"dt": sim.dt, "substeps": sim.substeps, "omega_global": sim.shell.omega_global,
58+ "omega_global_iterations": sim.shell.omega_global_iterations}
5859 fx["run"], fx["run_linear"] = run, run_linear
5960 taps = sim.taps
6061 fx["decimation"] = {"span": taps.shape[0], "substeps": taps.shape[1],