cymsim: nonlinear FEM cymbal simulation package
26 changed files+3059−0
.gitignoreadded+6−0View file
@@ -0,0 +1,6 @@
1+__pycache__/
2+*.egg-info/
3+.pytest_cache/
4+*.wav
5+*.npz
6+*_cascade.png
README.mdadded+61−0View file
@@ -0,0 +1,61 @@
1+# cymsim
2+
3+Explicit nonlinear finite-element simulation of struck cymbals. The cymbal is a Marguerre-von Karman shallow shell with a bell, a bow, a thickness taper and seeded imperfections, discretized on a triangulated disk with constant-strain membrane triangles and Discrete Kirchhoff bending triangles, integrated in time with an explicit leapfrog scheme, struck by one or more rigid exciters through power-law contact springs, and heard at a listener point through a Rayleigh-integral radiation model. The equations, the discretization and the time-stepping method are written out in [docs/model.md](docs/model.md).
4+
5+## Installation
6+
7+```
8+pip install -e ".[dev]" # numpy, scipy, jax (CPU), pyyaml, pytest, matplotlib
9+pip install "jax[cuda12]" # optional: run the time loop on a GPU
10+```
11+
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`.
13+
14+## Usage
15+
16+```
17+cymsim example-config my_run.yaml --exciter mallet # write a default config
18+cymsim check my_run.yaml # mesh, time step, element self-tests, no time stepping
19+cymsim simulate my_run.yaml # run; writes <prefix>_rayleigh.wav, _unbaffled.wav, _result.npz, _cascade.png
20+cymsim simulate my_run.yaml --linear --duration 0.5 # linear reference run
21+```
22+
23+The same from Python:
24+
25+```python
26+from cymsim import SimulationConfig, Hit, StickExciter, simulate
27+
28+cfg = SimulationConfig.from_yaml("examples/crash_mallet.yaml")
29+cfg.hits.append(Hit(time=0.8, r=0.19, phi=2.0, velocity=7.0, exciter=StickExciter()))
30+result = simulate(cfg)
31+result.write_wav("crash")
32+print(result.hits) # peak force, contact time and rebound velocity of each hit
33+print(result.energy_log) # kinetic / membrane / bending / exciter energy and max|w|/h over time
34+```
35+
36+A run is described by one `SimulationConfig` made of `Material`, `Profile`, `Imperfections`, `MeshConfig`, a list of `Hit`s (each with a `StickExciter` or `MalletExciter`), `Listener`, `Numerics` and `Output`. The YAML file mirrors the dataclass fields; omitted fields take the defaults, which reproduce an 18 inch medium-thin B20 crash. Units are SI throughout. See [examples/](examples/) for a mallet stroke, a stick crash, a train of four hits and a quick CPU run.
37+
38+## Layout
39+
40+| module | contents |
41+|---|---|
42+| [config.py](src/cymsim/config.py) | configuration dataclasses, YAML load/save, validation |
43+| [geometry.py](src/cymsim/geometry.py) | mid-surface height z0(x, y) and thickness h(x, y), seeded imperfections |
44+| [mesh.py](src/cymsim/mesh.py) | numpy disk mesher (concentric rings, angular stitching, smoothing), CST operators |
45+| [elements.py](src/cymsim/elements.py) | DKT bending stiffness, Marguerre-von Karman membrane strains, element energies |
46+| [shell.py](src/cymsim/shell.py) | assembled shell: lumped mass, rigid modes, internal force by autodiff, stability bound |
47+| [exciter.py](src/cymsim/exciter.py) | contact patches, power-law contact law, schedule of hits |
48+| [radiation.py](src/cymsim/radiation.py) | Rayleigh integral and the unbaffled (rim-diffraction) approximation |
49+| [simulate.py](src/cymsim/simulate.py) | leapfrog time loop, block compilation, result container |
50+| [diagnostics.py](src/cymsim/diagnostics.py) | band envelopes, onset times, spectral centroid, cascade plot |
51+| [cli.py](src/cymsim/cli.py) | `cymsim simulate`, `cymsim check`, `cymsim example-config` |
52+
53+`prototype.py` is the original single-file prototype that this package reorganizes; it is kept for reference and is not used by the package.
54+
55+## Tests
56+
57+```
58+pytest
59+```
60+
61+The tests check the mesher (manifold, orientation, quality), the DKT element (rigid modes, constant curvature, twist), the membrane strains, the assembled shell (mass, rigid-mode annihilation, force = minus energy gradient, dissipative damping), the configuration round trip, and short coarse simulations with one and two hits.
docs/model.mdadded+243−0View file
@@ -0,0 +1,243 @@
1+# The cymsim model: equations, discretization and method
2+
3+This document states what cymsim solves and how. Section 2 gives the continuous model, section 3 the geometry, section 4 the spatial discretization, section 5 the time integration, section 6 the excitation, section 7 the radiation model and section 8 the diagnostics. Section 9 lists the limitations we are aware of. Each section names the module that implements it. Units are SI throughout.
4+
5+## 1. Overview
6+
7+A crash cymbal is a thin, shallow, curved bronze shell. Struck hard, it does not simply ring at its linear modes: the transverse motion is many times the thickness, membrane stretching couples the modes, and energy injected at low frequency by the stroke climbs to high frequency over tens of milliseconds. This upward cascade is what gives a crash its build-up and its noisy, bright sustain, and it is entirely a nonlinear effect. A model that is meant to reproduce it must therefore keep the geometric nonlinearity of the shell, resolve frequencies up to the audio band, and run long enough for the cascade to develop.
8+
9+cymsim does this with the smallest set of ingredients we found sufficient: the Marguerre-von Karman shallow-shell equations (the von Karman plate with an initial curvature), a flat-facet triangular finite-element discretization, an explicit leapfrog integrator with a stability-based sub-step rate, a rigid exciter coupled through a power-law contact spring, and a Rayleigh-integral listener. Everything from the element energies onward is written with JAX so that the internal force is the exact gradient of the discrete strain energy and the whole time step compiles to one kernel.
10+
11+## 2. Governing equations
12+
13+Implemented in [elements.py](../src/cymsim/elements.py) (strains, energies) and [shell.py](../src/cymsim/shell.py) (damping, mass).
14+
15+### 2.1 Kinematics
16+
17+Let $\Omega$ be the disk of radius $R$ in the $(x, y)$ plane. The undeformed mid-surface is the graph $z = z_0(x, y)$ of a *shallow* surface, meaning $|\nabla z_0| \ll 1$; we warn when $\max |\nabla z_0|$ exceeds 0.4. The displacement of a mid-surface point has in-plane components $u(x, y, t)$, $v(x, y, t)$ along $x$ and $y$ and a transverse component $w(x, y, t)$ along $z$. The Marguerre-von Karman membrane strains are
18+
19+$$
20+\varepsilon_{xx} = u_{,x} + \tfrac12 w_{,x}^2 + z_{0,x}\, w_{,x}, \qquad
21+\varepsilon_{yy} = v_{,y} + \tfrac12 w_{,y}^2 + z_{0,y}\, w_{,y}, \qquad
22+\gamma_{xy} = u_{,y} + v_{,x} + w_{,x} w_{,y} + z_{0,x}\, w_{,y} + z_{0,y}\, w_{,x},
23+$$
24+
25+and the curvatures are those of the plate, $\kappa = (w_{,xx},\, w_{,yy},\, 2 w_{,xy})$. Setting $z_0 = 0$ recovers the von Karman plate. Dropping the quadratic terms $\tfrac12 w_{,x}^2$, $\tfrac12 w_{,y}^2$, $w_{,x} w_{,y}$ gives the *linear shallow shell* with the same geometry; the configuration flag `numerics.nonlinear = false` does exactly this and is the reference against which the cascade is judged. The terms in $z_0$ are the linear coupling between transverse motion and membrane stretch that a curved shell has and a flat plate lacks; they are what make the bell and the bow stiffen the low modes.
26+
27+### 2.2 Constitutive relations and strain energy
28+
29+The material is isotropic with Young's modulus $E$, Poisson ratio $\nu$ and density $\rho$; the thickness $h(x, y)$ varies. Membrane forces and bending moments are $N = C \varepsilon$ and $M = D \kappa$ with
30+
31+$$
32+C = \frac{E h}{1 - \nu^2} \begin{pmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \tfrac{1 - \nu}{2} \end{pmatrix}, \qquad
33+D = \frac{E h^3}{12 (1 - \nu^2)} \begin{pmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \tfrac{1 - \nu}{2} \end{pmatrix}.
34+$$
35+
36+The strain energy is
37+
38+$$
39+U[u, v, w] = \frac12 \int_\Omega \left( \varepsilon^{\mathsf T} C\, \varepsilon + \kappa^{\mathsf T} D\, \kappa \right) dA ,
40+$$
41+
42+the first term being the membrane energy and the second the bending energy. Their ratio during a run is one of the diagnostics of how nonlinear the motion is.
43+
44+### 2.3 Equations of motion and damping
45+
46+With $\rho h$ the mass per unit area and $\delta U / \delta(\cdot)$ the variational derivative,
47+
48+$$
49+\rho h\, \ddot u = -\frac{\delta U}{\delta u} - 2 \sigma_0 \rho h\, \dot u, \qquad
50+\rho h\, \ddot v = -\frac{\delta U}{\delta v} - 2 \sigma_0 \rho h\, \dot v, \qquad
51+\rho h\, \ddot w = -\frac{\delta U}{\delta w} - 2 \sigma_0 \rho h\, \dot w + 2 \sigma_1 \rho h\, \nabla^2 \dot w + f_c ,
52+$$
53+
54+where $f_c$ is the contact pressure of the exciter (section 6). Rotary inertia is neglected in the continuum model (it is reintroduced in modified form by the discretization, section 4.6). The losses follow the two-parameter form used throughout the sound-synthesis literature (Bilbao 2009): $\sigma_0$ gives a frequency-independent decay rate, so that the 60 dB decay time of the low modes is $T_{60} \approx 6.91 / \sigma_0$, and the $\sigma_1$ term gives a wave of wavenumber $k$ the additional decay rate $2 \sigma_1 k^2$. For bending waves $k^2 = \omega \sqrt{\rho h / D}$, so this rate grows linearly with frequency and sets how far up the cascade can reach before losses win. We make no attempt to model the physical loss mechanisms (thermoelastic, viscoelastic, radiation) separately; $\sigma_0$ and $\sigma_1$ are tuning parameters that lump them.
55+
56+### 2.4 Support
57+
58+The shell is free at the rim. A real cymbal rests on a felt washer at the bell, which we do not model; instead the six rigid-body motions of the shell are removed by projection at every step (section 5.3). This is the simplest device that keeps the cymbal from drifting away under the impulse of the stroke.
59+
60+## 3. Geometry
61+
62+Implemented in [geometry.py](../src/cymsim/geometry.py); parameters in `Profile` and `Imperfections`.
63+
64+The profile is written in polar coordinates $(r, \phi)$ on the projected disk. The *shoulder* is the junction between the bell and the bow; its radius $r_{sh}(\phi) = R_{sh} (1 + \epsilon_{oval} \cos(2\phi + \phi_{oval}))$ carries an $m = 2$ ovality. The bow coordinate $s = \mathrm{clip}\big((R - r) / (R - r_{sh}(\phi)),\, 0,\, 1\big)$ runs from 0 at the rim to 1 at the shoulder, and $u = 1 - s$. With $r_b$ the distance from the (offset) bell centre and $b(r_b) = \tfrac12 (1 + \cos \pi t)$, $t = \mathrm{clip}\big((r_b - R_{top}) / (r_{sh} - R_{top}),\, 0,\, 1\big)$, a $C^1$ ramp that is 1 on the bell plateau and 0 beyond the shoulder, the mid-surface height is
65+
66+$$
67+z_0 = H_{bow} \left( 1 - u^{p} \right) + H_{bell}\, b(r_b) + u \sum_m a_m \cos(m \phi + \phi_m) + \sum_i A_i \exp\!\left( -\frac{|x - x_i|^2}{2 s_i^2} \right),
68+$$
69+
70+bow, bell, rim warp and hammer dents in that order. The bow exponent $p = 1$ gives a cone and $p > 1$ a bow that is curved near the rim and flat at the shoulder. The warp amplitudes are $a_m = 2 a_{warp} / m$ for the listed modes with random phases; it vanishes at the shoulder and grows linearly to the rim, so the mean rim height stays $z_0 = 0$. The dents are Gaussian bumps of random sign (a cymbal is hammered from both sides), with centres uniform over the annulus between the shoulder and the rim and with jittered depth and width. The thickness is
71+
72+$$
73+h = h_{rim} + (h_{sh} - h_{rim})\, s + (h_{bell} - h_{sh})\, b(r_b),
74+$$
75+
76+a linear taper along the bow and a smooth step into the bell. All random draws come from `numpy.random.default_rng(seed)` in a fixed order, so a run is reproducible from its configuration. The imperfections matter: a perfectly axisymmetric shell has degenerate mode pairs and symmetric coupling rules, and its cascade is markedly less rich than that of a real cymbal.
77+
78+## 4. Spatial discretization
79+
80+### 4.1 Mesh
81+
82+Implemented in [mesh.py](../src/cymsim/mesh.py). The finite elements live on the flat projected disk; the curvature enters only through $z_0$ in the strains. Nodes are placed on concentric rings of spacing $\Delta r = R / K$ with $K = \mathrm{round}(R / \ell)$ for a target element size $\ell$; ring $k$ carries $n_k = \mathrm{round}(2\pi r_k / \Delta r) \approx 6k$ nodes and alternate rings are rotated by half a step. Consecutive rings are stitched by an angular merge: walking both rings in increasing angle, the next triangle uses whichever ring's next node comes first, which yields $n_{k-1} + n_k$ triangles per annulus and a valid triangulation. A Laplacian smoothing pass (`mesh.smoothing`, default 1) then relaxes the interior nodes with the rim nodes held on the circle. The mean element quality (inradius to circumradius ratio, 1 for an equilateral triangle) is about 0.96 and the minimum about 0.7. More passes improve the mean shape but shrink the smallest elements and thereby the stable time step, by about 10 percent for ten passes; the explicit step is set by the smallest elements in the thick bell (section 5.2), so a mesh graded coarser in the bell would be a worthwhile refinement. The rim is the polygon through the outer ring, so the mesh area is slightly below $\pi R^2$.
83+
84+### 4.2 Degrees of freedom
85+
86+Each node $i$ carries five unknowns $q_i = (u_i, v_i, w_i, \theta_{x,i}, \theta_{y,i})$ with the DKT convention $\theta_x = w_{,y}$ and $\theta_y = -w_{,x}$. The state of the shell is the array $q$ of shape $(N, 5)$.
87+
88+### 4.3 Membrane part: constant-strain triangle with the discrete initial surface
89+
90+On each triangle $u$, $v$ and $w$ are interpolated linearly from their nodal values, so their gradients are constant: $\partial_x \phi = \sum_i G_{x,i} \phi_i$ with $G_x$, $G_y$ the usual CST operators. The initial surface is represented in the same space, that is, $z_{0,x}$ and $z_{0,y}$ on an element are the CST gradients of the nodal heights $z_0(x_i, y_i)$, not the analytic slope at the centroid. This choice is deliberate. A rigid tilt of the shell about the $x$ axis by an angle $a$ displaces a surface point by $v = -a z_0$ and $w = a y$ (section 4.7); its linearized strain is $\varepsilon_{yy} = v_{,y} + z_{0,y} w_{,y} = -a\, (G z_0) + a\, z_{0,y}$, which vanishes exactly only if $z_{0,y}$ is the CST gradient of the nodal heights. With the analytic slope, the mismatch on a dented, coarse element left a rigid tilt with a spurious restoring force of order $10^{-4}$ of a typical stiffness; with the discrete slope the residual is at round-off (about $10^{-16}$ in float64). The membrane energy of element $e$ is $U_{m,e} = \tfrac12 A_e\, \varepsilon^{\mathsf T} C_e\, \varepsilon$ with $C_e$ evaluated at the centroid thickness.
91+
92+### 4.4 Bending part: the Discrete Kirchhoff Triangle
93+
94+Bending uses the DKT element of Batoz, Bathe and Ho (1980). Rotations $\beta_x, \beta_y$ are interpolated quadratically over the triangle from the three nodal $(w, \theta_x, \theta_y)$ triples, with the mid-side values fixed by imposing the Kirchhoff condition (rotation equals slope) discretely along each edge. The resulting curvature operator $B(\xi, \eta)$, with $\kappa = B\, U_e$ for the nine bending DOFs $U_e$ of the element, is linear in the area coordinates, so the element stiffness
95+
96+$$
97+K_e = \int_{A_e} B^{\mathsf T} D_e\, B\; dA
98+$$
99+
100+is integrated exactly by the three-point rule. The bending energy is $U_{b,e} = \tfrac12 U_e^{\mathsf T} K_e U_e$. DKT does not define a transverse displacement field inside the element, only the nodal values, which is why the membrane strains use the CST interpolation of $w$; the pairing of a DKT bending triangle with a CST membrane triangle is the classical flat shell triangle, and the Marguerre terms supply the curvature coupling between them. The element stiffness is computed once in float64 by automatic differentiation of the shape functions and checked at build time: a rigid tilt gives zero force (residual about $10^{-17}$), the state $w = x^2/2$ has energy $\tfrac12 D A_e$ to $10^{-10}$, and the $15 \times 15$ tangent of the full element energy has exactly six zero eigenvalues (three translations and three rotations). `cymsim check` prints these.
101+
102+### 4.5 Internal force by automatic differentiation
103+
104+The discrete strain energy is $U_h(q) = \sum_e \big( U_{m,e}(q_e) + U_{b,e}(q_e) \big)$ with $q_e$ the $(3, 5)$ gather of the element's nodal values. The internal nodal force is
105+
106+$$
107+f_{int}(q) = -\frac{\partial U_h}{\partial q},
108+$$
109+
110+obtained by reverse-mode differentiation (`jax.grad`) of the summed energies; the gather, the per-element evaluation and the scatter are all inside the differentiated function. The force is therefore the exact gradient of the discrete energy, so the semi-discrete system is conservative up to the damping terms, and the energy diagnostic is meaningful. The tangent stiffness needed for the stability bound is the Hessian of the same function at $q = 0$.
111+
112+### 4.6 Lumped mass and rotational inertia
113+
114+The mass is lumped: element $e$ contributes $m_e = \rho h_e A_e / 3$ to the $u$, $v$ and $w$ entries of each of its nodes. The rotational entries receive
115+
116+$$
117+I_e = m_e \left( \frac{h_e^2}{12} + f A_e \right),
118+$$
119+
120+where $h_e^2 / 12$ is the physical rotary inertia per unit mass and $f A_e$ (default $f = 1/64$) is an added inertia. The addition is a numerical device: with the physical value alone, the rotational DOFs of DKT have eigenfrequencies far above any bending wave the mesh can represent and would dictate an unusably small time step. The price is that bending waves whose wavelength is comparable to the element size are slowed slightly, by a relative amount of order $(k h_{eff})^2 / 24$ with $h_{eff}^2 = h_e^2 + 12 f A_e$, which for 6 mm elements is a few percent at 20 kHz and well below one percent at 5 kHz. These are also the waves the mesh resolves least well, so we consider the trade acceptable; $f$ can be reduced at the cost of a smaller time step. A second, optional scaling `inplane_mass_scale` multiplies the in-plane mass; values above 1 slow the fast in-plane waves ($c_L \approx 3.8$ km/s in bronze) and relax the time step further, but they also alter the in-plane dynamics that mediate the nonlinear coupling, so the default is 1.
121+
122+### 4.7 Rigid-body modes of the curved shell
123+
124+The free shell has six rigid-body modes. Three are translations. The in-plane rotation is $(u, v) = (-y, x)$. A tilt about the $x$ axis by an angle $a$ carries the surface point $(x, y, z_0)$ to first order to $(x,\, y - a z_0,\, z_0 + a y)$, so the mode is $v = -z_0$, $w = y$, $\theta_x = 1$; likewise the tilt about $y$ is $u = z_0$, $w = -x$, $\theta_y = 1$. The in-plane components proportional to $z_0$ are what distinguish the curved shell from the plate and are essential: without them a tilt stretches the shell. The six vectors are made $M$-orthonormal by Gram-Schmidt. At build time we verify that the tangent stiffness annihilates all six to round-off. Note that this holds for the linearized operator only: the Marguerre-von Karman strains are invariant under a finite rotation only to first order, so a finite tilt of amplitude $a$ leaves a residual force of order $a^2$. This is a known property of the shallow-shell approximation, not of the discretization, and is harmless here because the projection of section 5.3 keeps the rigid amplitudes at zero.
125+
126+### 4.8 The damping operator
127+
128+The $\sigma_1$ term is discretized in weak form with the CST gradients. The nodal force on $w$ is
129+
130+$$
131+f_{\sigma_1, i} = 2 \sigma_1 \sum_{e \ni i} (\rho h A)_e\, G_{e,i}^{\mathsf T} \big( G_e\, \dot w_e \big),
132+$$
133+
134+which is the assembly of $2 \sigma_1 \rho h \int \nabla \dot w \cdot \nabla \phi_i\, dA$. Its power $\dot w^{\mathsf T} f_{\sigma_1} = 2 \sigma_1 \sum_e (\rho h A)_e |\nabla \dot w|^2 \geq 0$ is non-negative, so the term only ever removes energy, and a uniform velocity is undamped by it. The $\sigma_0$ term is diagonal and acts on all five DOFs.
135+
136+## 5. Time integration
137+
138+Implemented in [simulate.py](../src/cymsim/simulate.py).
139+
140+### 5.1 Leapfrog with implicit $\sigma_0$
141+
142+With $M$ the diagonal lumped mass, one sub-step of size $\Delta t$ is
143+
144+$$
145+v^{n+1} = P \left[ \frac{v^n + \Delta t\, M^{-1} \big( f_{int}(q^n) + f_c(q^n, z_s^n) - f_{\sigma_1}(v^n) \big)}{1 + 2 \sigma_0 \Delta t} \right], \qquad
146+q^{n+1} = q^n + \Delta t\, v^{n+1},
147+$$
148+
149+where $P$ is the rigid-mode projection of section 5.3. This is the symplectic Euler scheme in $(q, v)$, equivalent to the central-difference (Stormer-Verlet) scheme in $q$, second-order accurate and non-dissipative for the conservative part. The $\sigma_0$ term is treated implicitly, which costs nothing because it is diagonal; the $\sigma_1$ term is evaluated explicitly at $v^n$, which is accurate enough because $2 \sigma_1 k^2 \Delta t \ll 1$ for every wavenumber the mesh can carry at the default values. The exciter is advanced by the same symplectic Euler step (section 6.3).
150+
151+### 5.2 Stability and the choice of the time step
152+
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. 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.
154+
155+### 5.3 Rigid-mode projection
156+
157+After each velocity update, $v \leftarrow v - \sum_{k=1}^{6} r_k \big( r_k^{\mathsf T} M v \big)$ with the $M$-orthonormal rigid modes $r_k$ of section 4.7. This removes the linear and angular momentum that the stroke imparts, which is what a stand would do, and it also removes the slow rigid drift that round-off would otherwise accumulate. The projection changes only the rigid part of the velocity and hence removes only rigid kinetic energy; the vibrational energy is untouched. We emphasize that this is a proxy for the support, not a model of it: a felt washer at the bell would also add a local stiffness and some damping.
158+
159+### 5.4 Precision and compilation
160+
161+Set-up (element stiffness, eigenvalue bound, rigid modes) is done in float64; the time loop runs in float32 by default (`numerics.dtype`), which is what a GPU does quickly. At the default parameters the displacement increment $\Delta t\, v$ is about $10^{-3}$ of the displacement itself, so float32 loses a few of its 24 bits per step in $q$ without visible effect on the energy balance; float64 is available for verification. The step is compiled once with `jax.jit` as a `fori_loop` over `steps_per_block` sub-steps, and the Python loop runs over blocks, reporting energies between them.
162+
163+## 6. Excitation
164+
165+Implemented in [exciter.py](../src/cymsim/exciter.py); parameters in `Hit`, `StickExciter`, `MalletExciter`.
166+
167+### 6.1 Contact law
168+
169+Each exciter is a rigid body of mass $m$ that moves only along $z$, with face height $z_s(t)$, coupled to the shell through $P$ springs at patch points $p$. The surface height under point $p$ is $w_p + o_p$, where $w_p = \sum_i \lambda_{p,i} w_i$ is the barycentric interpolation of the nodal transverse displacements in the containing element and $o_p$ is the height of the initial surface above the exciter's face plane, the best-fit plane through the patch points' $z_0$ values. The penetration and the force at point $p$ are
170+
171+$$
172+\delta_p = \max(w_p + o_p - z_s,\, 0), \qquad F_p = \frac{k}{P}\, \delta_p^{\alpha}, \qquad F = \sum_p F_p .
173+$$
174+
175+The plate receives $-F_p$ on the $w$ DOFs of the containing element, distributed with the weights $\lambda_{p,i}$, and the exciter obeys $m \ddot z_s = F$. The exponent $\alpha = 1.5$ is the Hertz law for a hard sphere on a hard surface; $\alpha \approx 2$ to $3$ fits felt. The contact spring is conservative (potential $k \delta^{\alpha+1} / (P(\alpha + 1))$), so the exciter rebounds with whatever energy the shell has not absorbed; there is no contact hysteresis, which slightly overestimates the rebound of a felt mallet. The residual offsets $o_p$ mean that the face touches the curved surface at a few points first, which is what a rigid face on a curved shell does.
176+
177+### 6.2 Contact patches
178+
179+A `StickExciter` contacts along a line of length $\ell$ through the strike point, at an angle `patch_angle` from the radial direction, with $P$ equally spaced springs and a given total stiffness $k$. A `MalletExciter` contacts on a disk: one spring at the centre and six on a hexagon of the head radius, and its stiffness is set from a target compression $\delta_c$ at full stop on a rigid wall, $\tfrac12 m v^2 = k \delta_c^{\alpha + 1} / (\alpha + 1)$. The set-up log reports, for each hit, the contact duration the exciter would have against a rigid wall,
180+
181+$$
182+\tau = \frac{2 c(\alpha)\, \delta_{max}}{v}, \qquad c(\alpha) = \int_0^1 \frac{dx}{\sqrt{1 - x^{\alpha + 1}}},
183+$$
184+
185+which for $\alpha = 1.5$ is Hertz's $2.94\, \delta_{max} / v$, together with the frequency $1 / (2\tau)$ below which the stroke itself injects energy. For the default mallet this is about 120 Hz, so everything the listener hears above a few hundred hertz has come up the cascade. On the shell the actual contact is longer and softer than on a rigid wall because the surface gives way.
186+
187+### 6.3 Schedule of hits
188+
189+A run has a list of hits, each with its own time, position, velocity and exciter. Every exciter exists for the whole run as a row of fixed-shape arrays (hits $\times$ largest patch, padded with zero weights and a very negative offset so that padded points never touch), which lets one compiled step serve any number of hits. Before its release step $n_k = \mathrm{round}(t_k / \Delta t)$ an exciter is inactive: it produces no force and is held far above the shell. At step $n_k$ it is placed at $z_s = \max_p (w_p + o_p) + \text{gap}$, just above the *current* surface, with velocity $-v_k$. This matters for a train of hits, where the surface under a later hit is already moving by several thickness units. After the rebound the exciter simply moves away at constant speed; there is no gravity and no hand. Two consequences follow: an exciter never strikes twice on its own, and if the shell were to overtake a slowly rebounding exciter a second contact would occur (a real player pulls the stick back). The summary of each hit reports the peak force, the total time with $F > 0$ and the rebound velocity.
190+
191+## 7. Radiation
192+
193+Implemented in [radiation.py](../src/cymsim/radiation.py); parameters in `Listener`.
194+
195+The pressure at the listener point $x_L$ from a baffled planar radiator is given by the Rayleigh integral, which we discretize with one point source per element:
196+
197+$$
198+p(x_L, t) = \frac{\rho_{air}}{2\pi} \sum_e \frac{A_e\, \ddot w_e\!\left( t - r_e / c \right)}{r_e},
199+$$
200+
201+where $\ddot w_e$ is the mean of the three nodal accelerations and $r_e$ the distance from the element centroid, placed on the curved surface at height $z_0$, to the listener. This is exact for a rigid piston in an infinite baffle and an approximation here: the surface is shallow rather than flat, the acceleration is taken along $z$ rather than along the local normal, and the air does not load the shell (radiation damping is part of what $\sigma_0$ and $\sigma_1$ lump).
202+
203+A cymbal has no baffle. Its back face radiates with the opposite sign, and that sound reaches the listener by diffracting around the rim; at low frequency the two contributions nearly cancel, which is why an unbaffled plate is a weak low-frequency source. We estimate the back contribution with the same sum but with $r_e$ replaced by the shortest path over the rim, $r_e^{back} = \min_j \big( |x_e - x_j^{rim}| + |x_j^{rim} - x_L| \big)$ over 720 rim points, weighted by `back_diffraction_gain`, and subtract it: $p_{unbaffled} = p_{front} - g\, p_{back}$. This is a Kirchhoff-style edge-diffraction estimate, not a boundary-element solution, and it ignores shadowing of the front by the shell itself. Both signals are returned; the baffled one is a cleaner measure of the shell motion, the unbaffled one sounds more like a cymbal in a room.
204+
205+Retarded times are handled at sub-step resolution: the delay of element $e$ in sub-steps, $d_e = r_e / (c \Delta t)$, is split into an integer part and a fraction, and each element's contribution is added to two adjacent slots of a delay buffer with linear interpolation. At the end the buffers are decimated to the audio rate with a polyphase FIR resampler (`scipy.signal.resample_poly`), which also removes the sub-step-rate content above the audio Nyquist frequency.
206+
207+## 8. Diagnostics
208+
209+Implemented in [simulate.py](../src/cymsim/simulate.py) (energies) and [diagnostics.py](../src/cymsim/diagnostics.py) (spectra). During the run the kinetic, membrane, bending and exciter energies are logged at `n_progress` points together with $\max_i |w_i| / h_i$, the *nonlinearity meter*: below about 1 the run is effectively linear, at 3 to 5 the cascade is clearly developed. The sum of the four energies can only decrease between the release of a hit and the next one, which is a direct check of the integrator. After the run, an STFT of the unbaffled pressure (256-point windows, 32-sample hop, about 0.7 ms) gives per-octave band energy envelopes from 60 Hz to 16 kHz; the *onset* of a band is the first time it reaches 20 dB below its own peak, and the spectral centroid is tracked against time. In a cascade the onsets are staggered with the high bands last and the centroid climbs over tens of milliseconds; in the linear reference run the bands rise together and the centroid is flat.
210+
211+## 9. Limitations and future work
212+
213+The model is a shallow shell, so it is only as good as the assumption $|\nabla z_0| \ll 1$; the default profile reaches a slope of 0.39 at the bell, which is at the edge of what shallow-shell theory is meant for, and a steep bell would call for a proper curved-shell formulation. The elements are flat facets with the curvature entering through $z_0$ alone, and the membrane strain uses the linear interpolation of $w$ while bending uses DKT rotations; this is a standard pairing but not a fully consistent one, and its convergence with element size has not been studied here. The lumped rotational inertia is inflated for the sake of the time step (section 4.6). The support is a rigid-mode projection rather than a model of the stand. The damping has two parameters and no physics behind them. The contact is one-dimensional along $z$, with no friction, no tangential force, no hysteresis and no hand. The radiation is a Rayleigh sum with a heuristic back-side correction.
214+
215+The natural next steps are, roughly in order of expected benefit: a convergence study of the cascade statistics with element size and precision; a bell support with local stiffness and damping in place of the projection; a contact model with hysteretic losses for felt; and a proper unbaffled radiation model. A train of hits already works structurally, but the behaviour of the contact model against a strongly moving surface deserves its own validation.
216+
217+## 10. References
218+
219+- Batoz, J.-L., Bathe, K.-J., Ho, L.-W. (1980). A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15, 1771-1812.
220+- Marguerre, K. (1938). Zur Theorie der gekrummten Platte grosser Formanderung. Proceedings of the 5th International Congress for Applied Mechanics, 93-101.
221+- Bilbao, S. (2009). Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics. Wiley. (Plate and shell schemes; the $\sigma_0$, $\sigma_1$ loss model.)
222+- Chaigne, A., Touze, C., Thomas, O. (2005). Nonlinear vibrations and chaos in gongs and cymbals. Acoustical Science and Technology 26, 403-409.
223+- Touze, C., Bilbao, S., Cadot, O. (2012). Transition scenario to turbulence in thin vibrating plates. Journal of Sound and Vibration 331, 412-433.
224+- Ducceschi, M., Touze, C. (2015). Modal approach for nonlinear vibrations of damped impacted plates: application to sound synthesis of gongs and cymbals. Journal of Sound and Vibration 344, 313-331.
225+- Hughes, T. J. R. (2000). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover. (Element eigenvalue bound for explicit time stepping.)
226+- Hairer, E., Lubich, C., Wanner, G. (2006). Geometric Numerical Integration. Springer. (Symplectic Euler and Stormer-Verlet.)
227+- Pierce, A. D. (1989). Acoustics: An Introduction to Its Physical Principles and Applications. Acoustical Society of America. (Rayleigh integral.)
228+
229+## 11. Symbols and where to find them
230+
231+| symbol | meaning | configuration / code |
232+|---|---|---|
233+| $E, \nu, \rho$ | Young's modulus, Poisson ratio, density | `Material.E`, `.nu`, `.rho` |
234+| $\sigma_0, \sigma_1$ | loss parameters | `Material.sigma0`, `.sigma1` |
235+| $z_0(x, y)$, $h(x, y)$ | mid-surface height, thickness | `CymbalProfile.z0`, `.thickness` |
236+| $C_e, D_e, K_e$ | membrane stiffness, bending rigidity, DKT stiffness of element $e$ | `ShellModel.Cm_e`, `.D_e`, `.Ke` |
237+| $G_x, G_y$ | CST gradient operators | `Mesh.Gx`, `.Gy` |
238+| $M$, $I_e$, $f$ | lumped mass, rotational inertia, its scaling factor | `ShellModel.M`, `.I_e`, `Numerics.rot_inertia_factor` |
239+| $r_k$ | rigid-body modes | `ShellModel.Rm` |
240+| $\omega_{max}$, $\Delta t$, $m$ | stability bound, sub-step, sub-steps per sample | `ShellModel.omega_max`, `Simulator.dt`, `.substeps` |
241+| $k, \alpha, P$ | contact stiffness, exponent, springs per patch | `ContactPatch.k_total`, `.alpha`, `.n_points` |
242+| $z_s, v_s$ | exciter face height and velocity | state entries `zs`, `vs` |
243+| $p_{front}, p_{back}$ | Rayleigh front and diffracted back pressure | `SimulationResult.pressure_rayleigh`, `.pressure_unbaffled` |
examples/crash_mallet.yamladded+72−0View file
@@ -0,0 +1,72 @@
1+# Soft mallet, gong-like stroke on an 18 inch medium-thin crash (the prototype's default run).
2+# The mallet injects only low frequencies; everything above a few hundred Hz comes from the
3+# nonlinear cascade. Units are SI throughout (m, kg, s, Pa). See docs/model.md.
4+material:
5+ E: 110000000000.0
6+ nu: 0.34
7+ rho: 8700.0
8+ sigma0: 0.9
9+ sigma1: 0.0005
10+profile:
11+ radius: 0.2286
12+ bell_height: 0.013
13+ bell_top_radius: 0.025
14+ bell_shoulder_radius: 0.08
15+ bow_height: 0.03
16+ bow_power: 1.4
17+ h_rim: 0.001
18+ h_shoulder: 0.0017
19+ h_bell: 0.0026
20+ max_slope: 0.4
21+imperfections:
22+ seed: 1
23+ n_dents: 400
24+ dent_depth: 0.00025
25+ dent_sigma: 0.007
26+ dent_margin: 0.01
27+ rim_warp: 0.0006
28+ rim_warp_modes:
29+ - 2
30+ - 3
31+ - 4
32+ - 5
33+ shoulder_ovality: 0.03
34+ bell_offset: 0.0015
35+mesh:
36+ element_size: 0.006
37+ smoothing: 1
38+hits:
39+- time: 0.0
40+ r: 0.13716
41+ phi: 0.3
42+ velocity: 8.0
43+ exciter:
44+ kind: mallet
45+ mass: 0.3
46+ alpha: 2.5
47+ compression: 0.012
48+ radius: 0.012
49+ gap: 0.0001
50+listener:
51+ position:
52+ - 0.4
53+ - 0.3
54+ - 0.9
55+ c_air: 343.0
56+ rho_air: 1.2
57+ back_diffraction_gain: 1.0
58+numerics:
59+ fs: 44100
60+ duration: 2.0
61+ cfl_safety: 0.6
62+ rot_inertia_factor: 0.015625
63+ inplane_mass_scale: 1.0
64+ steps_per_block: 20000
65+ dtype: float32
66+ nonlinear: true
67+ n_progress: 20
68+output:
69+ prefix: crash_mallet
70+ wav: true
71+ npz: true
72+ cascade_plot: true
examples/crash_stick.yamladded+72−0View file
@@ -0,0 +1,72 @@
1+# Hard stick shoulder crash near the rim: Hertz contact (alpha = 1.5) along a 3 cm line.
2+material:
3+ E: 110000000000.0
4+ nu: 0.34
5+ rho: 8700.0
6+ sigma0: 0.9
7+ sigma1: 0.0005
8+profile:
9+ radius: 0.2286
10+ bell_height: 0.013
11+ bell_top_radius: 0.025
12+ bell_shoulder_radius: 0.08
13+ bow_height: 0.03
14+ bow_power: 1.4
15+ h_rim: 0.001
16+ h_shoulder: 0.0017
17+ h_bell: 0.0026
18+ max_slope: 0.4
19+imperfections:
20+ seed: 1
21+ n_dents: 400
22+ dent_depth: 0.00025
23+ dent_sigma: 0.007
24+ dent_margin: 0.01
25+ rim_warp: 0.0006
26+ rim_warp_modes:
27+ - 2
28+ - 3
29+ - 4
30+ - 5
31+ shoulder_ovality: 0.03
32+ bell_offset: 0.0015
33+mesh:
34+ element_size: 0.006
35+ smoothing: 1
36+hits:
37+- time: 0.0
38+ r: 0.19430999999999998
39+ phi: 0.3
40+ velocity: 7.0
41+ exciter:
42+ kind: stick
43+ mass: 0.03
44+ k: 500000000.0
45+ alpha: 1.5
46+ patch_length: 0.03
47+ patch_points: 7
48+ patch_angle: 0.6
49+ gap: 0.0001
50+listener:
51+ position:
52+ - 0.4
53+ - 0.3
54+ - 0.9
55+ c_air: 343.0
56+ rho_air: 1.2
57+ back_diffraction_gain: 1.0
58+numerics:
59+ fs: 44100
60+ duration: 2.0
61+ cfl_safety: 0.6
62+ rot_inertia_factor: 0.015625
63+ inplane_mass_scale: 1.0
64+ steps_per_block: 20000
65+ dtype: float32
66+ nonlinear: true
67+ n_progress: 20
68+output:
69+ prefix: crash_stick
70+ wav: true
71+ npz: true
72+ cascade_plot: true
examples/hit_train.yamladded+108−0View file
@@ -0,0 +1,108 @@
1+# A train of four hits: a mallet stroke, two stick crashes on the ringing cymbal, and a second
2+# mallet stroke. Each hit has its own exciter, released at its time just above the moving surface.
3+material:
4+ E: 110000000000.0
5+ nu: 0.34
6+ rho: 8700.0
7+ sigma0: 0.9
8+ sigma1: 0.0005
9+profile:
10+ radius: 0.2286
11+ bell_height: 0.013
12+ bell_top_radius: 0.025
13+ bell_shoulder_radius: 0.08
14+ bow_height: 0.03
15+ bow_power: 1.4
16+ h_rim: 0.001
17+ h_shoulder: 0.0017
18+ h_bell: 0.0026
19+ max_slope: 0.4
20+imperfections:
21+ seed: 1
22+ n_dents: 400
23+ dent_depth: 0.00025
24+ dent_sigma: 0.007
25+ dent_margin: 0.01
26+ rim_warp: 0.0006
27+ rim_warp_modes:
28+ - 2
29+ - 3
30+ - 4
31+ - 5
32+ shoulder_ovality: 0.03
33+ bell_offset: 0.0015
34+mesh:
35+ element_size: 0.006
36+ smoothing: 1
37+hits:
38+- time: 0.0
39+ r: 0.137
40+ phi: 0.3
41+ velocity: 8.0
42+ exciter:
43+ kind: mallet
44+ mass: 0.3
45+ alpha: 2.5
46+ compression: 0.012
47+ radius: 0.012
48+ gap: 0.0001
49+- time: 0.8
50+ r: 0.194
51+ phi: 2.0
52+ velocity: 7.0
53+ exciter:
54+ kind: stick
55+ mass: 0.03
56+ k: 500000000.0
57+ alpha: 1.5
58+ patch_length: 0.03
59+ patch_points: 7
60+ patch_angle: 0.6
61+ gap: 0.0001
62+- time: 1.4
63+ r: 0.194
64+ phi: 2.3
65+ velocity: 5.0
66+ exciter:
67+ kind: stick
68+ mass: 0.03
69+ k: 500000000.0
70+ alpha: 1.5
71+ patch_length: 0.03
72+ patch_points: 7
73+ patch_angle: 0.6
74+ gap: 0.0001
75+- time: 2.0
76+ r: 0.137
77+ phi: -1.0
78+ velocity: 6.0
79+ exciter:
80+ kind: mallet
81+ mass: 0.3
82+ alpha: 2.5
83+ compression: 0.012
84+ radius: 0.012
85+ gap: 0.0001
86+listener:
87+ position:
88+ - 0.4
89+ - 0.3
90+ - 0.9
91+ c_air: 343.0
92+ rho_air: 1.2
93+ back_diffraction_gain: 1.0
94+numerics:
95+ fs: 44100
96+ duration: 3.0
97+ cfl_safety: 0.6
98+ rot_inertia_factor: 0.015625
99+ inplane_mass_scale: 1.0
100+ steps_per_block: 20000
101+ dtype: float32
102+ nonlinear: true
103+ n_progress: 20
104+output:
105+ prefix: train
106+ wav: true
107+ npz: true
108+ cascade_plot: true
examples/quick_cpu.yamladded+71−0View file
@@ -0,0 +1,71 @@
1+# Coarse, short run that finishes in about a minute on a CPU. Element size 20 mm resolves only
2+# the low modes; use it to exercise the pipeline, not for listening.
3+material:
4+ E: 110000000000.0
5+ nu: 0.34
6+ rho: 8700.0
7+ sigma0: 0.9
8+ sigma1: 0.0005
9+profile:
10+ radius: 0.2286
11+ bell_height: 0.013
12+ bell_top_radius: 0.025
13+ bell_shoulder_radius: 0.08
14+ bow_height: 0.03
15+ bow_power: 1.4
16+ h_rim: 0.001
17+ h_shoulder: 0.0017
18+ h_bell: 0.0026
19+ max_slope: 0.4
20+imperfections:
21+ seed: 1
22+ n_dents: 100
23+ dent_depth: 0.00025
24+ dent_sigma: 0.007
25+ dent_margin: 0.01
26+ rim_warp: 0.0006
27+ rim_warp_modes:
28+ - 2
29+ - 3
30+ - 4
31+ - 5
32+ shoulder_ovality: 0.03
33+ bell_offset: 0.0015
34+mesh:
35+ element_size: 0.02
36+ smoothing: 1
37+hits:
38+- time: 0.0
39+ r: 0.13716
40+ phi: 0.3
41+ velocity: 8.0
42+ exciter:
43+ kind: mallet
44+ mass: 0.3
45+ alpha: 2.5
46+ compression: 0.012
47+ radius: 0.012
48+ gap: 0.0001
49+listener:
50+ position:
51+ - 0.4
52+ - 0.3
53+ - 0.9
54+ c_air: 343.0
55+ rho_air: 1.2
56+ back_diffraction_gain: 1.0
57+numerics:
58+ fs: 44100
59+ duration: 0.3
60+ cfl_safety: 0.6
61+ rot_inertia_factor: 0.015625
62+ inplane_mass_scale: 1.0
63+ steps_per_block: 2000
64+ dtype: float32
65+ nonlinear: true
66+ n_progress: 20
67+output:
68+ prefix: quick
69+ wav: true
70+ npz: true
71+ cascade_plot: true
prototype.pyadded+577−0View file
@@ -0,0 +1,577 @@
1+#!/usr/bin/env python3
2+"""
3+crash_cymbal_fem.py — explicit nonlinear FEM simulation of a struck crash cymbal with bell, bow,
4+taper, seeded imperfections, a hard-stick or soft-mallet exciter, a nonlinearity meter, and a
5+spectral-onset diagnostic for the nonlinear upward energy cascade (gong-like "low → high" build-up).
6+
7+Model
8+-----
9+* Marguerre–von Kármán shallow shell: in-plane displacements (u, v) + transverse w on a
10+ mid-surface with initial shape z0(x, y) and spatially varying thickness h(x, y). The FE lives
11+ on the flat projected domain.
12+* Geometry (z0_of / thickness_of): bow + bell + taper, plus
13+ - hammering: N_DENTS random Gaussian dents on the bow, both signs,
14+ - rim warp: low-order azimuthal harmonics growing from the shoulder to the rim,
15+ - shoulder ovality (m = 2 modulation of the bell/bow junction) and a bell-centre offset.
16+ All seeded via SEED, all differentiable, so ∇z0, slope check and rigid modes follow.
17+* Membrane: constant-strain triangle (CST) with the full quadratic Marguerre–von Kármán strain.
18+ NONLINEAR = False removes the quadratic terms (linear shallow shell, same geometry) for
19+ reference runs.
20+* Bending: Discrete Kirchhoff Triangle (DKT, Batoz–Bathe–Ho 1980), nodal DOFs (w, θx, θy),
21+ θx = w,y θy = −w,x.
22+* Internal force = −∇(strain energy) via jax.grad, gather → per-element → scatter.
23+* Time: explicit leapfrog, lumped mass with scaled rotational inertia; dt from a rigorous
24+ element-eigenvalue (Irons) bound. Sub-step rate is an integer multiple of the audio rate.
25+* Losses: σ0 (frequency-independent) + σ1∇² (∝ ω) Rayleigh-type damping.
26+* Support: free–free, 6 rigid-body modes of the curved shell projected out each step.
27+* Exciter: rigid mass contacting through a patch of power-law springs F = k δ^α:
28+ "stick" — Hertz (α = 1.5), hard, ~0.3 ms contact, along a line (shoulder crash);
29+ "mallet" — felt (α ≈ 2.5), soft, ~3 ms contact, on a disk; injects only low frequencies,
30+ so anything above a few hundred Hz must come from the nonlinear cascade.
31+* Radiation at a listener point, retarded time at sub-step resolution:
32+ (a) Rayleigh integral (baffled),
33+ (b) unbaffled approximation: front Rayleigh − back Rayleigh routed via the shortest
34+ path over the rim (Kirchhoff-style edge diffraction; NOT a BEM solve).
35+* Diagnostics: membrane/bending energy split and max|w|/h per block (how nonlinear the run is);
36+ per-octave onset times, spectral centroid vs time, spectrogram + band-envelope PNG.
37+
38+Deps: pip install "jax[cuda12]" gmsh scipy numpy (matplotlib optional, for the plots)
39+"""
40+import math, time
41+import numpy as np
42+import jax
43+jax.config.update("jax_enable_x64", True) # set-up in f64; sim in DTYPE
44+jax.config.update("jax_default_matmul_precision", "highest") # no TF32 in the stiffness matvec
45+import jax.numpy as jnp
46+from jax import lax
47+import gmsh
48+from scipy.signal import resample_poly
49+from scipy.io import wavfile
50+
51+# ============================================================================ parameters
52+DTYPE = jnp.float32 # simulation precision
53+
54+# material — B20 bronze
55+E_MOD = 1.1e11 # Young's modulus [Pa]
56+NU = 0.34
57+RHO = 8700.0 # [kg/m^3]
58+SIGMA0 = 0.9 # [1/s] T60 of lows ≈ 6.9/σ0
59+SIGMA1 = 5e-4 # [m^2/s] ∝ω loss (T60 ≈ 0.6 s @ 5 kHz); sets the cascade ceiling
60+
61+# cymbal profile — 18" medium-thin crash, mean rim height z0 = 0
62+R = 0.2286 # radius [m]
63+BELL_H = 0.013 # bell height above the bow at the shoulder [m]
64+BELL_R_TOP = 0.025 # radius of the (flat) bell plateau [m]
65+BELL_R_SHOULDER = 0.080 # mean bell/bow junction radius [m]
66+BOW_H = 0.030 # bow rise from rim to shoulder [m]
67+BOW_P = 1.4 # bow power: 1 = cone, >1 = curved and flat at the shoulder
68+H_RIM = 1.0e-3 # thickness at the rim [m]
69+H_SHOULDER = 1.7e-3 # thickness at the bell/bow junction [m]
70+H_BELL = 2.6e-3 # thickness in the bell [m]
71+MAX_SLOPE = 0.4 # shallow-shell sanity bound on |∇z0|
72+
73+# symmetry breaking (all seeded)
74+SEED = 1
75+N_DENTS = 400 # hammer marks on the bow
76+DENT_DEPTH = 0.25e-3 # typical dent depth [m], random sign (hammered from both sides)
77+DENT_SIGMA = 7e-3 # dent Gaussian σ [m], ±30 % jitter (≈ 2 elements: keep ≥ H_ELEM)
78+DENT_R_MIN, DENT_R_MAX = BELL_R_SHOULDER + 0.01, R - 0.01
79+RIM_WARP = 0.6e-3 # m = 2 rim-warp amplitude [m]; mode m gets 2·RIM_WARP/m
80+RIM_WARP_MODES = (2, 3, 4, 5)
81+SHOULDER_OVALITY = 0.03 # relative m = 2 modulation of the shoulder radius
82+BELL_OFFSET = 1.5e-3 # bell-centre offset from the geometric centre [m]
83+
84+# mesh
85+H_ELEM = 6e-3 # target element size [m]
86+
87+# excitation
88+EXCITER = "mallet" # "stick": hard shoulder crash; "mallet": soft, gong-like
89+STRIKE_R, STRIKE_PHI = 0.60 * R, 0.3 # contact centre [m], [rad] (0.85 R for the stick)
90+STICK_GAP = 1e-4 # initial clearance above the highest patch point [m]
91+# stick (EXCITER = "stick")
92+V_STRIKE, M_STICK, K_HERTZ = 7.0, 0.030, 5e8 # [m/s], [kg], total Hertz k [N/m^1.5]
93+PATCH_LEN, PATCH_POINTS, PATCH_ANGLE = 0.03, 7, 0.6 # contact line [m], springs, dir. vs radial
94+# mallet (EXCITER = "mallet") — power-law felt contact F = k δ^α, k set from target compression
95+V_MALLET, M_MALLET = 8.0, 0.30 # [m/s], head + effective arm mass [kg] → ~10 J, tam-tam stroke
96+ALPHA_MALLET = 2.5 # felt exponent (1.5 = Hertz, 2–3 = felt)
97+MALLET_COMPRESSION = 12e-3 # felt compression at full stop on a rigid wall [m] → ~3 ms contact
98+MALLET_RADIUS = 0.012 # head contact radius [m], spread over 7 springs
99+
100+NONLINEAR = True # False → linear shallow shell (same geometry); reference run
101+
102+# listener & air (mean rim plane is z = 0, centre at origin; drummer's ear-ish)
103+LISTENER = np.array([0.40, 0.30, 0.90])
104+C_AIR, RHO_AIR = 343.0, 1.2
105+BACK_DIFFRACTION_GAIN = 1.0 # weight of the diffracted back-side term
106+
107+# time / numerics
108+FS = 44100
109+DURATION = 2.0 # [s]
110+CFL_SAFETY = 0.6 # margin for nonlinearity + damping (lower to 0.5 if energy grows)
111+ROT_INERTIA_FACTOR = 1 / 64 # lumped rotational inertia I = m (h²/12 + f·A_e)
112+INPLANE_MASS_SCALE = 1.0 # >1 slows in-plane waves, relaxes dt (see notes)
113+STEPS_PER_BLOCK = 20000
114+OUT_PREFIX = "crash"
115+
116+# ============================================================================ seeded imperfections
117+_rng = np.random.default_rng(SEED)
118+
119+def _uniform_annulus(n, r0, r1):
120+ r = np.sqrt(_rng.uniform(r0 ** 2, r1 ** 2, n))
121+ t = _rng.uniform(0, 2 * np.pi, n)
122+ return r * np.cos(t), r * np.sin(t)
123+
124+_dx, _dy = _uniform_annulus(N_DENTS, DENT_R_MIN, DENT_R_MAX)
125+DENT_X, DENT_Y = jnp.asarray(_dx), jnp.asarray(_dy)
126+DENT_A = jnp.asarray(DENT_DEPTH * _rng.choice([-1.0, 1.0], N_DENTS) * _rng.uniform(0.5, 1.5, N_DENTS))
127+DENT_S = jnp.asarray(DENT_SIGMA * _rng.uniform(0.7, 1.3, N_DENTS))
128+WARP_M = jnp.asarray(RIM_WARP_MODES, dtype=float)
129+WARP_A = jnp.asarray([RIM_WARP * 2.0 / m for m in RIM_WARP_MODES])
130+WARP_PH = jnp.asarray(_rng.uniform(0, 2 * np.pi, len(RIM_WARP_MODES)))
131+OVAL_PH = float(_rng.uniform(0, 2 * np.pi))
132+_bo = _rng.uniform(0, 2 * np.pi)
133+BELL_CX, BELL_CY = BELL_OFFSET * math.cos(_bo), BELL_OFFSET * math.sin(_bo)
134+
135+# ============================================================================ geometry hooks
136+def _polar(x, y):
137+ r = jnp.sqrt(x * x + y * y + 1e-30)
138+ phi = jnp.arctan2(y, jnp.where(r < 1e-6, 1e-6, x)) # finite gradient at the centre
139+ return r, phi
140+
141+def _shoulder_radius(phi):
142+ return BELL_R_SHOULDER * (1.0 + SHOULDER_OVALITY * jnp.cos(2 * phi + OVAL_PH))
143+
144+def _bell_blend(r_bell, r_sh):
145+ """1 on the bell plateau, 0 beyond the shoulder, C¹ cosine ramp in between."""
146+ t = jnp.clip((r_bell - BELL_R_TOP) / (r_sh - BELL_R_TOP), 0.0, 1.0)
147+ return 0.5 * (1.0 + jnp.cos(jnp.pi * t))
148+
149+def z0_of(x, y):
150+ """Mid-surface height z0(x, y) [m]. jnp-differentiable."""
151+ r, phi = _polar(x, y)
152+ r_sh = _shoulder_radius(phi)
153+ r_bell = jnp.sqrt((x - BELL_CX) ** 2 + (y - BELL_CY) ** 2 + 1e-30)
154+ s = jnp.clip((R - r) / (R - r_sh), 0.0, 1.0) # 0 at rim → 1 at shoulder
155+ u = jnp.maximum(1.0 - s, 1e-12) # 1 at rim → 0 at shoulder
156+ bow = BOW_H * (1.0 - u ** BOW_P)
157+ bell = BELL_H * _bell_blend(r_bell, r_sh)
158+ warp = u * jnp.sum(WARP_A * jnp.cos(WARP_M * phi + WARP_PH))
159+ dents = jnp.sum(DENT_A * jnp.exp(-((x - DENT_X) ** 2 + (y - DENT_Y) ** 2) / (2 * DENT_S ** 2)))
160+ return bow + bell + warp + dents
161+
162+def thickness_of(x, y):
163+ """Thickness h(x, y) [m]: linear taper along the bow, smooth step up into the bell."""
164+ r, phi = _polar(x, y)
165+ r_sh = _shoulder_radius(phi)
166+ r_bell = jnp.sqrt((x - BELL_CX) ** 2 + (y - BELL_CY) ** 2 + 1e-30)
167+ s = jnp.clip((R - r) / (R - r_sh), 0.0, 1.0)
168+ bow = H_RIM + (H_SHOULDER - H_RIM) * s
169+ return bow + (H_BELL - H_SHOULDER) * _bell_blend(r_bell, r_sh)
170+
171+_z0_grad = jax.vmap(jax.grad(lambda p: z0_of(p[0], p[1]))) # analytic ∇z0
172+_z0_eval = jax.vmap(z0_of)
173+_h_eval = jax.vmap(thickness_of)
174+
175+# ============================================================================ mesh (gmsh)
176+def make_disk_mesh(radius, h):
177+ gmsh.initialize()
178+ gmsh.option.setNumber("General.Terminal", 0)
179+ gmsh.model.add("cymbal")
180+ gmsh.model.occ.addDisk(0, 0, 0, radius, radius)
181+ gmsh.model.occ.synchronize()
182+ gmsh.option.setNumber("Mesh.MeshSizeMin", 0.8 * h)
183+ gmsh.option.setNumber("Mesh.MeshSizeMax", h)
184+ gmsh.option.setNumber("Mesh.Algorithm", 6) # Frontal-Delaunay
185+ gmsh.option.setNumber("Mesh.Smoothing", 10)
186+ gmsh.model.mesh.generate(2)
187+ tags, xyz, _ = gmsh.model.mesh.getNodes()
188+ _, tri_tags = gmsh.model.mesh.getElementsByType(2) # 3-node triangles
189+ gmsh.finalize()
190+ tags = np.asarray(tags, dtype=np.int64)
191+ xyz = np.asarray(xyz).reshape(-1, 3)
192+ remap = np.full(tags.max() + 1, -1, dtype=np.int64)
193+ remap[tags] = np.arange(len(tags))
194+ tri = remap[np.asarray(tri_tags, dtype=np.int64).reshape(-1, 3)]
195+ used = np.unique(tri)
196+ xy = xyz[used, :2]
197+ inv = np.full(len(tags), -1, dtype=np.int64); inv[used] = np.arange(len(used))
198+ tri = inv[tri]
199+ a = ((xy[tri[:, 1], 0] - xy[tri[:, 0], 0]) * (xy[tri[:, 2], 1] - xy[tri[:, 0], 1])
200+ - (xy[tri[:, 2], 0] - xy[tri[:, 0], 0]) * (xy[tri[:, 1], 1] - xy[tri[:, 0], 1]))
201+ tri[a < 0] = tri[a < 0][:, [0, 2, 1]] # counter-clockwise
202+ return xy, tri
203+
204+# ============================================================================ DKT element
205+def dkt_stiffness(xe, ye, D, nu):
206+ """9x9 DKT bending stiffness, DOFs (w1,θx1,θy1, w2,...), θx = w,y, θy = −w,x."""
207+ i, j = jnp.array([1, 2, 0]), jnp.array([2, 0, 1]) # sides k=4,5,6 ↔ ij = 23, 31, 12
208+ xij, yij = xe[i] - xe[j], ye[i] - ye[j]
209+ l2 = xij ** 2 + yij ** 2
210+ a = -xij / l2
211+ b = 0.75 * xij * yij / l2
212+ c = (0.25 * xij ** 2 - 0.5 * yij ** 2) / l2
213+ d = -yij / l2
214+ e = (0.25 * yij ** 2 - 0.5 * xij ** 2) / l2
215+ a4, a5, a6 = a; b4, b5, b6 = b; c4, c5, c6 = c; d4, d5, d6 = d; e4, e5, e6 = e
216+
217+ def H(s):
218+ xi, eta = s
219+ lam = 1 - xi - eta
220+ N1, N2, N3 = lam * (2 * lam - 1), xi * (2 * xi - 1), eta * (2 * eta - 1)
221+ N4, N5, N6 = 4 * xi * eta, 4 * eta * lam, 4 * xi * lam
222+ Hx = jnp.stack([1.5 * (a6 * N6 - a5 * N5), b5 * N5 + b6 * N6, N1 - c5 * N5 - c6 * N6,
223+ 1.5 * (a4 * N4 - a6 * N6), b6 * N6 + b4 * N4, N2 - c4 * N4 - c6 * N6,
224+ 1.5 * (a5 * N5 - a4 * N4), b4 * N4 + b5 * N5, N3 - c4 * N4 - c5 * N5])
225+ Hy = jnp.stack([1.5 * (d6 * N6 - d5 * N5), -N1 + e5 * N5 + e6 * N6, -b5 * N5 - b6 * N6,
226+ 1.5 * (d4 * N4 - d6 * N6), -N2 + e4 * N4 + e6 * N6, -b6 * N6 - b4 * N4,
227+ 1.5 * (d5 * N5 - d4 * N4), -N3 + e4 * N4 + e5 * N5, -b4 * N4 - b5 * N5])
228+ return jnp.stack([Hx, Hy]) # (2, 9)
229+
230+ x21, x31 = xe[1] - xe[0], xe[2] - xe[0]
231+ y21, y31 = ye[1] - ye[0], ye[2] - ye[0]
232+ det = x21 * y31 - x31 * y21 # = 2A
233+ xi_x, xi_y, eta_x, eta_y = y31 / det, -x31 / det, -y21 / det, x21 / det
234+
235+ def B_at(s):
236+ J = jax.jacfwd(H)(s) # (2, 9, 2)
237+ Hx_xi, Hx_eta, Hy_xi, Hy_eta = J[0, :, 0], J[0, :, 1], J[1, :, 0], J[1, :, 1]
238+ Hx_x = Hx_xi * xi_x + Hx_eta * eta_x; Hx_y = Hx_xi * xi_y + Hx_eta * eta_y
239+ Hy_x = Hy_xi * xi_x + Hy_eta * eta_x; Hy_y = Hy_xi * xi_y + Hy_eta * eta_y
240+ return jnp.stack([Hx_x, Hy_y, Hx_y + Hy_x]) # curvature operator (3, 9)
241+
242+ Db = D * jnp.array([[1, nu, 0], [nu, 1, 0], [0, 0, 0.5 * (1 - nu)]])
243+ gps = jnp.array([[1 / 6, 1 / 6], [2 / 3, 1 / 6], [1 / 6, 2 / 3]]) # exact for quadratic B
244+ Bs = jax.vmap(B_at)(gps) # (3, 3, 9)
245+ return (0.5 * det / 3) * jnp.einsum("gki,kl,glj->ij", Bs, Db, Bs)
246+
247+# ============================================================================ element energy
248+def elem_energy_parts(qe, Gx, Gy, A, z0x, z0y, Cm, Ke):
249+ """qe: (3,5) nodal [u,v,w,θx,θy]. Returns (membrane, bending) energy of one element.
250+ The z0 terms are the *linear* curvature coupling; NONLINEAR gates only the quadratic terms."""
251+ nl = 1.0 if NONLINEAR else 0.0
252+ ux, uy = Gx @ qe[:, 0], Gy @ qe[:, 0]
253+ vx, vy = Gx @ qe[:, 1], Gy @ qe[:, 1]
254+ wx, wy = Gx @ qe[:, 2], Gy @ qe[:, 2]
255+ exx = ux + nl * 0.5 * wx * wx + z0x * wx
256+ eyy = vy + nl * 0.5 * wy * wy + z0y * wy
257+ gxy = uy + vx + nl * wx * wy + z0x * wy + z0y * wx
258+ C11, C12, C33 = Cm
259+ Em = 0.5 * A * (C11 * (exx * exx + eyy * eyy) + 2 * C12 * exx * eyy + C33 * gxy * gxy)
260+ U = qe[:, 2:].reshape(9)
261+ return Em, 0.5 * (U @ Ke @ U)
262+
263+def elem_energy(*args):
264+ Em, Eb = elem_energy_parts(*args)
265+ return Em + Eb
266+
267+# ============================================================================ diagnostics
268+def spectral_onset(p, fs, label, bands=(60, 125, 250, 500, 1000, 2000, 4000, 8000, 16000),
269+ t_max=0.5):
270+ """Per-band energy envelopes at ~6 ms resolution; onset = first time a band reaches −20 dB
271+ of its own peak. Also plots the envelopes (log-log), which shows the cascade far better than
272+ a spectrogram normalised to the low-frequency peak."""
273+ from scipy.signal import stft
274+ f_, t_, Z = stft(p, fs, nperseg=256, noverlap=224)
275+ Pw = np.abs(Z) ** 2
276+ tot_max = Pw.sum(0).max()
277+ envs = []
278+ print(f"\n[{label}] band onset(-20 dB) peak level vs loudest band")
279+ for lo, hi in zip(bands[:-1], bands[1:]):
280+ m = (f_ >= lo) & (f_ < hi)
281+ env = Pw[m].sum(0)
282+ if env.max() <= 0:
283+ continue
284+ env_db = 10 * np.log10(env / env.max() + 1e-30)
285+ onset = t_[np.argmax(env_db > -20)]
286+ print(f" {lo:5d}–{hi:5d} Hz {onset * 1e3:7.1f} ms {t_[np.argmax(env)] * 1e3:7.1f} ms"
287+ f" {10 * np.log10(env.max() / tot_max):6.1f} dB")
288+ envs.append((lo, hi, env))
289+ cen = (Pw * f_[:, None]).sum(0) / (Pw.sum(0) + 1e-30)
290+ for tq in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5):
291+ if tq < t_[-1]:
292+ print(f" spectral centroid at {tq * 1e3:5.0f} ms: {cen[np.searchsorted(t_, tq)]:7.0f} Hz")
293+ try:
294+ import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
295+ fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
296+ sel = t_ <= t_max
297+ ax[0].pcolormesh(t_[sel] * 1e3, f_, 10 * np.log10(Pw[:, sel] / Pw.max() + 1e-12),
298+ vmin=-90, vmax=0, shading="auto")
299+ ax[0].set_yscale("log"); ax[0].set_ylim(50, fs / 2)
300+ ax[0].set_xlabel("t [ms]"); ax[0].set_ylabel("f [Hz]"); ax[0].set_title(label)
301+ for lo, hi, env in envs:
302+ ax[1].plot(t_[sel] * 1e3, 10 * np.log10(env[sel] / tot_max + 1e-30), label=f"{lo}–{hi} Hz")
303+ ax[1].set_ylim(-90, 5); ax[1].set_xscale("log"); ax[1].set_xlim(2, t_max * 1e3)
304+ ax[1].set_xlabel("t [ms]"); ax[1].set_ylabel("band energy [dB re loudest]")
305+ ax[1].legend(fontsize=7, ncol=2); ax[1].grid(alpha=.3)
306+ ax[1].set_title("band envelopes — a cascade = staggered rises, high bands last")
307+ plt.tight_layout(); plt.savefig(f"{label}_cascade.png", dpi=120); plt.close()
308+ print(f" wrote {label}_cascade.png")
309+ except Exception:
310+ pass
311+
312+# ============================================================================ main
313+def main():
314+ t_start = time.time()
315+ tag = OUT_PREFIX if NONLINEAR else OUT_PREFIX + "_linear"
316+ xy, tri = make_disk_mesh(R, H_ELEM)
317+ N, E = len(xy), len(tri)
318+ print(f"mesh: {N} nodes, {E} triangles, {5 * N} DOFs (NONLINEAR = {NONLINEAR})")
319+
320+ # ---- element geometry on the projected domain (numpy, f64)
321+ x, y = xy[tri, 0], xy[tri, 1] # (E, 3)
322+ x21, x31 = x[:, 1] - x[:, 0], x[:, 2] - x[:, 0]
323+ y21, y31 = y[:, 1] - y[:, 0], y[:, 2] - y[:, 0]
324+ area = 0.5 * (x21 * y31 - x31 * y21)
325+ assert (area > 0).all()
326+ Gx = np.stack([y[:, 1] - y[:, 2], y[:, 2] - y[:, 0], y[:, 0] - y[:, 1]], 1) / (2 * area[:, None])
327+ Gy = np.stack([x[:, 2] - x[:, 1], x[:, 0] - x[:, 2], x[:, 1] - x[:, 0]], 1) / (2 * area[:, None])
328+ xc, yc = x.mean(1), y.mean(1)
329+
330+ def locate(px, py):
331+ """Containing element and barycentric weights of a point on the projected domain."""
332+ l1 = ((y[:, 1] - y[:, 2]) * (px - x[:, 2]) + (x[:, 2] - x[:, 1]) * (py - y[:, 2])) / (2 * area)
333+ l2 = ((y[:, 2] - y[:, 0]) * (px - x[:, 2]) + (x[:, 0] - x[:, 2]) * (py - y[:, 2])) / (2 * area)
334+ l3 = 1 - l1 - l2
335+ inside = (l1 >= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9)
336+ assert inside.any(), f"point ({px:.3f}, {py:.3f}) outside mesh"
337+ es = int(np.argmax(inside))
338+ return tri[es], np.array([l1[es], l2[es], l3[es]])
339+
340+ # ---- profile: thickness, initial shape and its analytic gradient
341+ h_e = np.asarray(_h_eval(jnp.asarray(xc), jnp.asarray(yc)))
342+ h_node = np.asarray(_h_eval(jnp.asarray(xy[:, 0]), jnp.asarray(xy[:, 1])))
343+ z0_node = np.asarray(_z0_eval(jnp.asarray(xy[:, 0]), jnp.asarray(xy[:, 1])))
344+ grad_z0 = np.asarray(_z0_grad(jnp.stack([jnp.asarray(xc), jnp.asarray(yc)], 1)))
345+ z0x, z0y = grad_z0[:, 0], grad_z0[:, 1]
346+ slope = np.hypot(z0x, z0y).max()
347+ rim_node = np.hypot(xy[:, 0], xy[:, 1]) > R - 1e-6
348+ print(f"profile: height {z0_node.max() * 1e3:.1f} mm, rim warp p-p "
349+ f"{(z0_node[rim_node].max() - z0_node[rim_node].min()) * 1e3:.2f} mm, "
350+ f"{N_DENTS} dents, h ∈ [{h_e.min() * 1e3:.2f}, {h_e.max() * 1e3:.2f}] mm, "
351+ f"max |∇z0| = {slope:.3f}")
352+ if slope > MAX_SLOPE:
353+ print(f"WARNING: max slope {slope:.2f} exceeds shallow-shell bound {MAX_SLOPE}")
354+
355+ D_e = E_MOD * h_e ** 3 / (12 * (1 - NU ** 2))
356+ Cm_e = (E_MOD * h_e / (1 - NU ** 2))[:, None] * np.array([1.0, NU, 0.5 * (1 - NU)])
357+
358+ # ---- DKT stiffness for every element
359+ Ke64 = jax.vmap(dkt_stiffness, in_axes=(0, 0, 0, None))(jnp.asarray(x), jnp.asarray(y),
360+ jnp.asarray(D_e), NU)
361+ Ke_np = np.asarray(Ke64)
362+
363+ # element self-tests: rigid tilt → zero force; w = ½x² → energy ½·D·A
364+ U_tilt = np.stack([y, np.ones_like(y), np.zeros_like(y)], 2).reshape(E, 9)
365+ res = np.abs(np.einsum("eij,ej->ei", Ke_np, U_tilt)).max() / np.abs(Ke_np).max()
366+ U_cur = np.stack([0.5 * x ** 2, np.zeros_like(x), -x], 2).reshape(E, 9)
367+ Eb = 0.5 * np.einsum("ei,eij,ej->e", U_cur, Ke_np, U_cur)
368+ err = np.abs(Eb / (0.5 * D_e * area) - 1).max()
369+ print(f"DKT checks: rigid-tilt residual {res:.1e}, constant-curvature energy error {err:.1e}")
370+
371+ # ---- lumped mass (per DOF): [u, v, w, θx, θy]
372+ m_e = RHO * h_e * area / 3
373+ I_e = m_e * (h_e ** 2 / 12 + ROT_INERTIA_FACTOR * area)
374+ M = np.zeros((N, 5))
375+ for k in range(3):
376+ np.add.at(M[:, 0], tri[:, k], m_e * INPLANE_MASS_SCALE)
377+ np.add.at(M[:, 1], tri[:, k], m_e * INPLANE_MASS_SCALE)
378+ np.add.at(M[:, 2], tri[:, k], m_e)
379+ np.add.at(M[:, 3], tri[:, k], I_e)
380+ np.add.at(M[:, 4], tri[:, k], I_e)
381+ print(f"plate mass {M[:, 2].sum():.3f} kg")
382+
383+ # ---- stable dt from the element eigenvalue bound (Irons): λ_max(M⁻¹K) ≤ max_e λ_max(M_e⁻¹K_e)
384+ Khess = jax.vmap(jax.hessian(elem_energy))(
385+ jnp.zeros((E, 3, 5)), jnp.asarray(Gx), jnp.asarray(Gy), jnp.asarray(area),
386+ jnp.asarray(z0x), jnp.asarray(z0y), jnp.asarray(Cm_e), Ke64).reshape(E, 15, 15)
387+ Me = np.tile(np.stack([m_e * INPLANE_MASS_SCALE, m_e * INPLANE_MASS_SCALE, m_e, I_e, I_e], 1),
388+ (1, 3)) # (E, 15) node-major
389+ S = np.asarray(Khess) / np.sqrt(Me[:, :, None] * Me[:, None, :])
390+ lam = np.linalg.eigvalsh(S)
391+ n_zero = int((lam[0] < 1e-9 * lam[0].max()).sum())
392+ print(f"element 0 zero-energy modes: {n_zero} (expect 6)")
393+ e_worst = int(np.argmax(lam[:, -1]))
394+ omega_max = math.sqrt(lam[:, -1].max())
395+ dt_bound = 2.0 / omega_max
396+ M_SUB = math.ceil(1.0 / (FS * CFL_SAFETY * dt_bound))
397+ dt = 1.0 / (FS * M_SUB)
398+ n_blocks = math.ceil(DURATION * FS * M_SUB / STEPS_PER_BLOCK)
399+ n_steps = n_blocks * STEPS_PER_BLOCK
400+ print(f"ω_max bound {omega_max:.3e} rad/s (limiting element at r = "
401+ f"{math.hypot(xc[e_worst], yc[e_worst]):.3f} m, h = {h_e[e_worst] * 1e3:.2f} mm)")
402+ print(f"dt = {dt:.3e} s, {M_SUB} sub-steps/sample, {n_steps} steps")
403+
404+ # ---- rigid-body modes of the curved shell, M-orthonormalised
405+ # A rigid tilt drags the in-plane displacements along: rotation about x by φ gives
406+ # v = −φ z0, w = φ y; about y gives u = +φ z0, w = −φ x (first order, shallow-shell).
407+ X, Y = xy[:, 0] - xy[:, 0].mean(), xy[:, 1] - xy[:, 1].mean()
408+ Rm = np.zeros((6, N, 5))
409+ Rm[0, :, 0] = 1; Rm[1, :, 1] = 1 # in-plane translations
410+ Rm[2, :, 0] = -Y; Rm[2, :, 1] = X # in-plane rotation
411+ Rm[3, :, 2] = 1 # z translation
412+ Rm[4, :, 2] = Y; Rm[4, :, 3] = 1; Rm[4, :, 1] = -z0_node # tilt about x (θx = w,y)
413+ Rm[5, :, 2] = -X; Rm[5, :, 4] = 1; Rm[5, :, 0] = z0_node # tilt about y (θy = −w,x)
414+ for k in range(6):
415+ for l in range(k):
416+ Rm[k] -= (Rm[l] * M * Rm[k]).sum() * Rm[l]
417+ Rm[k] /= math.sqrt((Rm[k] * M * Rm[k]).sum())
418+
419+ # ---- contact patch: a set of power-law springs F = k δ^α on a rigid exciter of mass M_EX
420+ xs, ys = STRIKE_R * math.cos(STRIKE_PHI), STRIKE_R * math.sin(STRIKE_PHI)
421+ if EXCITER == "stick":
422+ M_EX, V_EX, ALPHA, K_TOTAL = M_STICK, V_STRIKE, 1.5, K_HERTZ
423+ ang = STRIKE_PHI + PATCH_ANGLE
424+ xi = np.linspace(-0.5, 0.5, PATCH_POINTS) * PATCH_LEN
425+ px, py = xs + xi * math.cos(ang), ys + xi * math.sin(ang)
426+ elif EXCITER == "mallet":
427+ M_EX, V_EX, ALPHA = M_MALLET, V_MALLET, ALPHA_MALLET
428+ # k such that the head stops after MALLET_COMPRESSION on a rigid wall
429+ K_TOTAL = (ALPHA + 1) * 0.5 * M_EX * V_EX ** 2 / MALLET_COMPRESSION ** (ALPHA + 1)
430+ th = np.arange(6) * np.pi / 3
431+ px = np.concatenate([[xs], xs + MALLET_RADIUS * np.cos(th)])
432+ py = np.concatenate([[ys], ys + MALLET_RADIUS * np.sin(th)])
433+ else:
434+ raise ValueError(f"unknown EXCITER {EXCITER!r}")
435+ P = len(px)
436+ loc = [locate(px[k], py[k]) for k in range(P)]
437+ s_idx = np.array([l[0] for l in loc]) # (P, 3)
438+ s_wts = np.array([l[1] for l in loc]) # (P, 3)
439+ z0_p = np.asarray(_z0_eval(jnp.asarray(px), jnp.asarray(py)))
440+ # rigid exciter face: best-fit plane through the surface heights; residual → staggered contact
441+ Apl = np.column_stack([px - xs, py - ys, np.ones(P)])
442+ s_off = z0_p - Apl @ np.linalg.lstsq(Apl, z0_p, rcond=None)[0]
443+ zs0 = s_off.max() + STICK_GAP
444+ tau_est = 2.0 * MALLET_COMPRESSION / V_EX if EXCITER == "mallet" else None
445+ print(f"{EXCITER}: {P}-point patch, k = {K_TOTAL:.3g} N/m^{ALPHA}, m = {M_EX} kg, v = {V_EX} m/s, "
446+ f"energy {0.5 * M_EX * V_EX ** 2:.2f} J"
447+ + (f", free contact ≈ {tau_est * 1e3:.1f} ms → direct excitation below ~{0.5 / tau_est:.0f} Hz"
448+ if tau_est else ""))
449+
450+ # ---- radiation geometry (element centroids on the curved surface, retarded times in sub-steps)
451+ cent = np.stack([xc, yc, z0_node[tri].mean(1)], 1) # (E, 3)
452+ r_front = np.linalg.norm(LISTENER - cent, axis=1)
453+ phi = np.linspace(0, 2 * np.pi, 720, endpoint=False)
454+ rim_xy = np.stack([R * np.cos(phi), R * np.sin(phi)], 1)
455+ rim_z = np.asarray(_z0_eval(jnp.asarray(rim_xy[:, 0]), jnp.asarray(rim_xy[:, 1])))
456+ rim = np.column_stack([rim_xy, rim_z])
457+ d1 = np.linalg.norm(cent[:, None, :] - rim[None, :, :], axis=2)
458+ d2 = np.linalg.norm(rim - LISTENER, axis=1)
459+ r_back = (d1 + d2[None, :]).min(1)
460+ coef_f = RHO_AIR * area / (2 * np.pi * r_front)
461+ coef_b = BACK_DIFFRACTION_GAIN * RHO_AIR * area / (2 * np.pi * r_back)
462+ df, db = r_front / (C_AIR * dt), r_back / (C_AIR * dt)
463+ i0f, frf = np.floor(df).astype(np.int32), df - np.floor(df)
464+ i0b, frb = np.floor(db).astype(np.int32), db - np.floor(db)
465+ buf_len = n_steps + int(max(i0f.max(), i0b.max())) + 2
466+
467+ # ---- device arrays in simulation precision
468+ f = lambda a: jnp.asarray(a, dtype=DTYPE)
469+ tri_j = jnp.asarray(tri, dtype=jnp.int32)
470+ Gx_j, Gy_j, area_j = f(Gx), f(Gy), f(area)
471+ z0x_j, z0y_j, Cm_j, Ke_j = f(z0x), f(z0y), f(Cm_e), f(Ke_np)
472+ rhoh_area_j, hnode_j = f(RHO * h_e * area), f(h_node)
473+ M_j, Minv_j, Rm_j = f(M), f(1.0 / M), f(Rm)
474+ sidx_j, swts_j, soff_j = jnp.asarray(s_idx, dtype=jnp.int32), f(s_wts), f(s_off)
475+ coef_f_j, coef_b_j, frf_j, frb_j = f(coef_f), f(coef_b), f(frf), f(frb)
476+ i0f_j, i0b_j = jnp.asarray(i0f), jnp.asarray(i0b)
477+ k_point = K_TOTAL / P
478+
479+ def total_energy(q):
480+ return jax.vmap(elem_energy)(q[tri_j], Gx_j, Gy_j, area_j, z0x_j, z0y_j, Cm_j, Ke_j).sum()
481+
482+ f_int = jax.grad(lambda q: -total_energy(q)) # internal nodal forces
483+
484+ # shell check: tangent stiffness at rest must annihilate all six rigid modes
485+ q0 = jnp.zeros((N, 5), DTYPE)
486+ K_apply = lambda d: jax.jvp(f_int, (q0,), (d,))[1]
487+ chk = np.random.default_rng(0)
488+ worst = 0.0
489+ for k in range(6):
490+ rnd = chk.standard_normal((N, 5)); rnd *= np.linalg.norm(Rm[k]) / np.linalg.norm(rnd)
491+ worst = max(worst, float(jnp.abs(K_apply(Rm_j[k])).max() / jnp.abs(K_apply(f(rnd))).max()))
492+ print(f"rigid modes: max |K·r| / |K·random| = {worst:.1e} (should be ~1e-6 or smaller)")
493+
494+ def damp_lap(vw):
495+ """−2σ1 ρh ∇²ẇ as a nodal force: 2σ1 · (ρh A_e) · Σ_e G_iᵀ (G ẇ_e)."""
496+ ve = vw[tri_j]
497+ vx, vy = (Gx_j * ve).sum(1), (Gy_j * ve).sum(1)
498+ c = rhoh_area_j[:, None] * (Gx_j * vx[:, None] + Gy_j * vy[:, None])
499+ return 2 * SIGMA1 * jax.ops.segment_sum(c.ravel(), tri_j.ravel(), num_segments=N)
500+
501+ def project_rigid(v):
502+ coef = jnp.einsum("knd,nd->k", Rm_j, M_j * v)
503+ return v - jnp.einsum("k,knd->nd", coef, Rm_j)
504+
505+ def step(n, st):
506+ q, v, zs, vs, bf, bb, Fmax, ncont = st
507+ # exciter / plate contact along the patch
508+ wp = jnp.einsum("pk,pk->p", swts_j, q[sidx_j, 2]) + soff_j # surface vs exciter face
509+ delta = jnp.maximum(wp - zs, 0.0) # (P,)
510+ Fk = k_point * jnp.power(delta, ALPHA)
511+ F = Fk.sum()
512+ # forces
513+ fo = f_int(q)
514+ fo = fo.at[sidx_j, 2].add(-Fk[:, None] * swts_j)
515+ fo = fo.at[:, 2].add(-damp_lap(v[:, 2]))
516+ # leapfrog (σ0 implicit), rigid-mode projection
517+ v1 = (v + dt * fo * Minv_j) / (1 + 2 * SIGMA0 * dt)
518+ v1 = project_rigid(v1)
519+ q1 = q + dt * v1
520+ vs1 = vs + dt * F / M_EX
521+ zs1 = zs + dt * vs1
522+ # radiation: scatter element accelerations into retarded-time buffers
523+ acc = ((v1[:, 2] - v[:, 2]) / dt)[tri_j].mean(1)
524+ sf, sb = coef_f_j * acc, coef_b_j * acc
525+ bf = bf.at[n + i0f_j].add(sf * (1 - frf_j)).at[n + i0f_j + 1].add(sf * frf_j)
526+ bb = bb.at[n + i0b_j].add(sb * (1 - frb_j)).at[n + i0b_j + 1].add(sb * frb_j)
527+ return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0))
528+
529+ @jax.jit
530+ def run_block(n0, st):
531+ return lax.fori_loop(0, STEPS_PER_BLOCK,
532+ lambda k, s: step((n0 + k).astype(jnp.int32), s), st)
533+
534+ @jax.jit
535+ def energies(st):
536+ q, v = st[0], st[1]
537+ Em, Eb = jax.vmap(elem_energy_parts)(q[tri_j], Gx_j, Gy_j, area_j, z0x_j, z0y_j, Cm_j, Ke_j)
538+ return (0.5 * (M_j * v * v).sum(), Em.sum(), Eb.sum(), 0.5 * M_EX * st[3] ** 2,
539+ jnp.abs(q[:, 2] / hnode_j).max())
540+
541+ state = (jnp.zeros((N, 5), DTYPE), jnp.zeros((N, 5), DTYPE),
542+ jnp.asarray(zs0, DTYPE), jnp.asarray(-V_EX, DTYPE),
543+ jnp.zeros(buf_len, DTYPE), jnp.zeros(buf_len, DTYPE),
544+ jnp.asarray(0.0, DTYPE), jnp.asarray(0, jnp.int32))
545+
546+ print(f"setup done in {time.time() - t_start:.1f} s; compiling + running on "
547+ f"{jax.devices()[0].platform} ...")
548+ print(" (nonlinearity meter: max|w|/h ≲ 1 → effectively linear; 3–5 → clear cascade)")
549+ t_run = time.time()
550+ for blk in range(n_blocks):
551+ state = run_block(jnp.asarray(blk * STEPS_PER_BLOCK, jnp.int32), state)
552+ if blk % max(1, n_blocks // 20) == 0 or blk == n_blocks - 1:
553+ ke, em, eb, ks, woh = energies(state)
554+ jax.block_until_ready(ke)
555+ print(f" t = {(blk + 1) * STEPS_PER_BLOCK * dt:6.3f} s KE {float(ke):.3e} "
556+ f"membrane {float(em):.3e} bending {float(eb):.3e} J exciter {float(ks):.3e} J "
557+ f"max|w|/h = {float(woh):5.2f} [{time.time() - t_run:.0f} s wall]")
558+
559+ q, v, zs, vs, bf, bb, Fmax, ncont = jax.device_get(state)
560+ print(f"{EXCITER}: peak force {float(Fmax):.1f} N, contact {float(ncont) * dt * 1e3:.2f} ms, "
561+ f"rebound velocity {float(vs):+.2f} m/s "
562+ f"(energy left in exciter {0.5 * M_EX * float(vs) ** 2:.2f} of {0.5 * M_EX * V_EX ** 2:.2f} J)")
563+
564+ # ---- decimate to audio rate, write, analyse
565+ n_out = int(FS * DURATION)
566+ p_ray = resample_poly(np.asarray(bf[:n_steps], np.float64), 1, M_SUB)[:n_out]
567+ p_back = resample_poly(np.asarray(bb[:n_steps], np.float64), 1, M_SUB)[:n_out]
568+ p_free = p_ray - p_back
569+ for name, p in (("rayleigh", p_ray), ("unbaffled", p_free)):
570+ peak = np.abs(p).max()
571+ print(f"{name:10s}: peak {peak:.3f} Pa ({20 * np.log10(peak / 2e-5):.1f} dB SPL peak)")
572+ wavfile.write(f"{tag}_{name}.wav", FS, (0.9 * p / peak * 32767).astype(np.int16))
573+ print(f"wrote {tag}_rayleigh.wav, {tag}_unbaffled.wav (total {time.time() - t_start:.0f} s)")
574+ spectral_onset(p_free, FS, f"{tag}_unbaffled")
575+
576+if __name__ == "__main__":
577+ main()
pyproject.tomladded+25−0View file
@@ -0,0 +1,25 @@
1+[build-system]
2+requires = ["setuptools>=68"]
3+build-backend = "setuptools.build_meta"
4+
5+[project]
6+name = "cymsim"
7+version = "0.1.0"
8+description = "Explicit nonlinear finite-element simulation of struck cymbals"
9+readme = "README.md"
10+requires-python = ">=3.10"
11+dependencies = ["numpy", "scipy", "jax", "pyyaml"]
12+
13+[project.optional-dependencies]
14+plot = ["matplotlib"]
15+cuda = ["jax[cuda12]"]
16+dev = ["pytest", "matplotlib"]
17+
18+[project.scripts]
19+cymsim = "cymsim.cli:main"
20+
21+[tool.setuptools.packages.find]
22+where = ["src"]
23+
24+[tool.pytest.ini_options]
25+testpaths = ["tests"]
src/cymsim/__init__.pyadded+46−0View file
@@ -0,0 +1,46 @@
1+"""cymsim: explicit nonlinear finite-element simulation of struck cymbals.
2+
3+The model is a Marguerre-von Karman shallow shell on a triangulated disk, discretized with
4+constant-strain membrane triangles and Discrete Kirchhoff bending triangles, integrated in time
5+with an explicit leapfrog scheme, struck by one or more rigid exciters through power-law contact
6+springs, and heard through a Rayleigh-integral radiation model. The equations, the
7+discretization and the time-stepping method are written out in docs/model.md; each module's
8+docstring points to the relevant section.
9+
10+Typical use::
11+
12+ from cymsim import SimulationConfig, simulate
13+ cfg = SimulationConfig.from_yaml("examples/crash_mallet.yaml")
14+ result = simulate(cfg)
15+ result.write_wav("crash")
16+"""
17+import jax
18+
19+# Set-up quantities (element stiffness, eigenvalue bounds, rigid modes) are computed in float64;
20+# the time loop casts everything to the precision chosen in Numerics.dtype.
21+jax.config.update("jax_enable_x64", True)
22+# Never let a matmul silently drop to TF32 inside the stiffness evaluation.
23+jax.config.update("jax_default_matmul_precision", "highest")
24+
25+from .config import ( # noqa: E402
26+ Hit,
27+ Imperfections,
28+ Listener,
29+ MalletExciter,
30+ Material,
31+ MeshConfig,
32+ Numerics,
33+ Output,
34+ Profile,
35+ SimulationConfig,
36+ StickExciter,
37+)
38+from .simulate import SimulationResult, Simulator, simulate # noqa: E402
39+
40+__version__ = "0.1.0"
41+
42+__all__ = [
43+ "Hit", "Imperfections", "Listener", "MalletExciter", "Material", "MeshConfig", "Numerics",
44+ "Output", "Profile", "SimulationConfig", "StickExciter",
45+ "SimulationResult", "Simulator", "simulate",
46+]
src/cymsim/cli.pyadded+106−0View file
@@ -0,0 +1,106 @@
1+"""Command line interface: ``cymsim simulate``, ``cymsim check`` and ``cymsim example-config``."""
2+from __future__ import annotations
3+
4+import argparse
5+import logging
6+import sys
7+
8+from .config import SimulationConfig, default_config
9+
10+
11+def _setup_logging(quiet: bool) -> None:
12+ logging.basicConfig(level=logging.WARNING if quiet else logging.INFO,
13+ format="%(message)s", stream=sys.stdout)
14+
15+
16+def _load(path: str, args) -> SimulationConfig:
17+ cfg = SimulationConfig.from_yaml(path)
18+ if getattr(args, "duration", None) is not None:
19+ cfg.numerics.duration = args.duration
20+ if getattr(args, "out_prefix", None) is not None:
21+ cfg.output.prefix = args.out_prefix
22+ if getattr(args, "linear", False):
23+ cfg.numerics.nonlinear = False
24+ return cfg
25+
26+
27+def cmd_simulate(args) -> int:
28+ from .diagnostics import analyze_cascade, format_cascade_report, plot_cascade
29+ from .simulate import Simulator
30+
31+ cfg = _load(args.config, args)
32+ sim = Simulator(cfg)
33+ result = sim.run()
34+ prefix = cfg.output.prefix if cfg.numerics.nonlinear else cfg.output.prefix + "_linear"
35+ if cfg.output.wav:
36+ print("wrote " + ", ".join(result.write_wav(prefix)))
37+ if cfg.output.npz:
38+ print("wrote " + result.write_npz(prefix))
39+ analysis = analyze_cascade(result.pressure_unbaffled, result.fs)
40+ print(format_cascade_report(analysis, f"{prefix}_unbaffled"))
41+ if cfg.output.cascade_plot:
42+ try:
43+ print("wrote " + plot_cascade(analysis, f"{prefix}_cascade.png", f"{prefix}_unbaffled", result.fs))
44+ except ImportError:
45+ print("matplotlib not available; skipping the cascade plot")
46+ print(f"total wall time {sim.setup_time + result.wall_time:.0f} s")
47+ return 0
48+
49+
50+def cmd_check(args) -> int:
51+ """Build the model and report mesh, time step and element checks without time stepping."""
52+ from .simulate import Simulator
53+
54+ cfg = _load(args.config, args)
55+ sim = Simulator(cfg)
56+ c = sim.shell.element_checks()
57+ print(f"DKT checks: rigid-tilt residual {c.rigid_tilt_residual:.1e}, "
58+ f"constant-curvature energy error {c.constant_curvature_error:.1e}")
59+ print(f"zero-energy modes per element: {c.zero_energy_modes_min} to {c.zero_energy_modes_max} (expect 6)")
60+ print(f"rigid modes: max |K r| / |K random| = {c.rigid_mode_ratio:.1e} (should be ~1e-6 or smaller)")
61+ n = cfg.numerics
62+ print(f"a {n.duration:.2f} s run needs {sim.n_steps} sub-steps of {sim.mesh.n_nodes * 5} DOFs")
63+ return 0
64+
65+
66+def cmd_example_config(args) -> int:
67+ cfg = default_config(args.exciter)
68+ if args.path:
69+ cfg.to_yaml(args.path)
70+ print(f"wrote {args.path}")
71+ else:
72+ import yaml
73+ sys.stdout.write(yaml.safe_dump(cfg.to_dict(), sort_keys=False))
74+ return 0
75+
76+
77+def main(argv=None) -> int:
78+ ap = argparse.ArgumentParser(prog="cymsim", description="Simulate a struck cymbal.")
79+ sub = ap.add_subparsers(dest="command", required=True)
80+
81+ p = sub.add_parser("simulate", help="run a simulation from a YAML config")
82+ p.add_argument("config")
83+ p.add_argument("--duration", type=float, help="override numerics.duration [s]")
84+ p.add_argument("--out-prefix", help="override output.prefix")
85+ p.add_argument("--linear", action="store_true", help="drop the quadratic strain terms (reference run)")
86+ p.add_argument("--quiet", action="store_true")
87+ p.set_defaults(func=cmd_simulate)
88+
89+ p = sub.add_parser("check", help="build the model and print mesh, time-step and element checks")
90+ p.add_argument("config")
91+ p.add_argument("--quiet", action="store_true")
92+ p.set_defaults(func=cmd_check)
93+
94+ p = sub.add_parser("example-config", help="print or write a default config")
95+ p.add_argument("path", nargs="?")
96+ p.add_argument("--exciter", choices=("mallet", "stick"), default="mallet")
97+ p.add_argument("--quiet", action="store_true")
98+ p.set_defaults(func=cmd_example_config)
99+
100+ args = ap.parse_args(argv)
101+ _setup_logging(args.quiet)
102+ return args.func(args)
103+
104+
105+if __name__ == "__main__":
106+ sys.exit(main())
src/cymsim/config.pyadded+263−0View file
@@ -0,0 +1,263 @@
1+"""Configuration dataclasses and their YAML representation.
2+
3+Every physical and numerical parameter of a run lives in one of the dataclasses below, grouped
4+by what it describes: the material, the cymbal profile, the seeded imperfections, the mesh, the
5+list of hits (each with its own exciter), the listener, the numerics and the output. Defaults
6+reproduce the 18-inch medium-thin crash of the original prototype. A SimulationConfig can be
7+built directly in Python or loaded from a YAML file whose keys mirror the field names; a hit's
8+exciter is chosen by its ``kind`` key ("stick" or "mallet").
9+"""
10+from __future__ import annotations
11+
12+import dataclasses
13+import math
14+import types
15+import typing
16+from dataclasses import dataclass, field
17+from pathlib import Path
18+
19+import yaml
20+
21+
22+# ------------------------------------------------------------------------------- material
23+@dataclass
24+class Material:
25+ """Isotropic material of the shell (default: B20 bronze) and its two loss parameters."""
26+ E: float = 1.1e11 # Young's modulus [Pa]
27+ nu: float = 0.34 # Poisson ratio
28+ rho: float = 8700.0 # density [kg/m^3]
29+ sigma0: float = 0.9 # frequency-independent loss [1/s]; T60 of the lows ~ 6.9/sigma0
30+ sigma1: float = 5e-4 # loss proportional to omega [m^2/s]; sets the ceiling of the cascade
31+
32+
33+# -------------------------------------------------------------------------------- profile
34+@dataclass
35+class Profile:
36+ """Axisymmetric part of the cymbal: bow, bell and thickness taper. Mean rim height is z0 = 0."""
37+ radius: float = 0.2286 # outer radius [m] (18 inch)
38+ bell_height: float = 0.013 # bell height above the bow at the shoulder [m]
39+ bell_top_radius: float = 0.025 # radius of the flat bell plateau [m]
40+ bell_shoulder_radius: float = 0.080 # mean radius of the bell/bow junction [m]
41+ bow_height: float = 0.030 # bow rise from the rim to the shoulder [m]
42+ bow_power: float = 1.4 # bow shape exponent: 1 = cone, > 1 = curved, flat at the shoulder
43+ h_rim: float = 1.0e-3 # thickness at the rim [m]
44+ h_shoulder: float = 1.7e-3 # thickness at the shoulder [m]
45+ h_bell: float = 2.6e-3 # thickness in the bell [m]
46+ max_slope: float = 0.4 # warn when max |grad z0| exceeds this (shallow-shell bound)
47+
48+
49+@dataclass
50+class Imperfections:
51+ """Seeded symmetry breaking: hammer dents, rim warp, shoulder ovality and a bell offset."""
52+ seed: int = 1
53+ n_dents: int = 400 # Gaussian hammer marks on the bow, random sign
54+ dent_depth: float = 0.25e-3 # typical dent depth [m] (x 0.5 to 1.5 jitter)
55+ dent_sigma: float = 7e-3 # dent Gaussian sigma [m] (x 0.7 to 1.3 jitter); keep >= element size
56+ dent_margin: float = 0.01 # dents stay this far from the shoulder and the rim [m]
57+ rim_warp: float = 0.6e-3 # m = 2 rim-warp amplitude [m]; mode m gets 2 * rim_warp / m
58+ rim_warp_modes: tuple[int, ...] = (2, 3, 4, 5)
59+ shoulder_ovality: float = 0.03 # relative m = 2 modulation of the shoulder radius
60+ bell_offset: float = 1.5e-3 # bell centre offset from the geometric centre [m]
61+
62+
63+@dataclass
64+class MeshConfig:
65+ element_size: float = 6e-3 # target edge length of the triangles [m]
66+ smoothing: int = 1 # Laplacian smoothing passes; more improves shape, shrinks the smallest elements
67+
68+
69+# ------------------------------------------------------------------------------- exciters
70+@dataclass
71+class StickExciter:
72+ """Hard stick: Hertz contact (alpha = 1.5) along a short line, as in a shoulder crash."""
73+ kind: str = "stick"
74+ mass: float = 0.030 # effective mass [kg]
75+ k: float = 5e8 # total contact stiffness [N/m^alpha]
76+ alpha: float = 1.5 # contact exponent
77+ patch_length: float = 0.03 # length of the contact line [m]
78+ patch_points: int = 7 # springs along the line
79+ patch_angle: float = 0.6 # direction of the line relative to the radial direction [rad]
80+ gap: float = 1e-4 # initial clearance above the surface at release [m]
81+
82+
83+@dataclass
84+class MalletExciter:
85+ """Soft felt mallet: power-law contact on a small disk; k is set from a target compression."""
86+ kind: str = "mallet"
87+ mass: float = 0.30 # head plus effective arm mass [kg]
88+ alpha: float = 2.5 # felt exponent (1.5 = Hertz, 2 to 3 = felt)
89+ compression: float = 12e-3 # felt compression at full stop on a rigid wall [m]
90+ radius: float = 0.012 # head contact radius [m]; 7 springs (centre + hexagon)
91+ gap: float = 1e-4
92+
93+
94+Exciter = typing.Union[StickExciter, MalletExciter]
95+EXCITER_KINDS: dict[str, type] = {"stick": StickExciter, "mallet": MalletExciter}
96+
97+
98+@dataclass
99+class Hit:
100+ """One strike: when, where, how fast, and with what."""
101+ time: float = 0.0 # release time of the exciter [s]
102+ r: float = 0.137 # radial position of the contact centre [m]
103+ phi: float = 0.3 # azimuth of the contact centre [rad]
104+ velocity: float = 8.0 # approach speed [m/s]
105+ exciter: Exciter = field(default_factory=MalletExciter)
106+
107+
108+# ------------------------------------------------------------------------------- listener
109+@dataclass
110+class Listener:
111+ """Listener point and air; the mean rim plane is z = 0 with the cymbal centre at the origin."""
112+ position: tuple[float, float, float] = (0.40, 0.30, 0.90) # [m]
113+ c_air: float = 343.0
114+ rho_air: float = 1.2
115+ back_diffraction_gain: float = 1.0 # weight of the diffracted back-side term (0 = baffled only)
116+
117+
118+# ------------------------------------------------------------------------------- numerics
119+@dataclass
120+class Numerics:
121+ fs: int = 44100 # audio sample rate [Hz]; the sub-step rate is an integer multiple
122+ duration: float = 2.0 # simulated time [s]
123+ cfl_safety: float = 0.6 # fraction of the linear stability limit actually used
124+ rot_inertia_factor: float = 1 / 64 # lumped rotational inertia I = m (h^2/12 + f A_e); see docs
125+ inplane_mass_scale: float = 1.0 # > 1 slows in-plane waves and relaxes the time step; see docs
126+ steps_per_block: int = 20000 # sub-steps per compiled block
127+ dtype: str = "float32" # simulation precision: "float32" or "float64"
128+ nonlinear: bool = True # False drops the quadratic strain terms (linear reference run)
129+ n_progress: int = 20 # number of progress/energy reports during the run
130+
131+
132+@dataclass
133+class Output:
134+ prefix: str = "crash" # output files are <prefix>_<name>.wav / .png / .npz
135+ wav: bool = True
136+ npz: bool = True
137+ cascade_plot: bool = True # needs matplotlib
138+
139+
140+# ----------------------------------------------------------------------------- simulation
141+@dataclass
142+class SimulationConfig:
143+ material: Material = field(default_factory=Material)
144+ profile: Profile = field(default_factory=Profile)
145+ imperfections: Imperfections = field(default_factory=Imperfections)
146+ mesh: MeshConfig = field(default_factory=MeshConfig)
147+ hits: list[Hit] = field(default_factory=lambda: [Hit()])
148+ listener: Listener = field(default_factory=Listener)
149+ numerics: Numerics = field(default_factory=Numerics)
150+ output: Output = field(default_factory=Output)
151+
152+ # ---- validation
153+ def validate(self) -> None:
154+ p, n = self.profile, self.numerics
155+ if not self.hits:
156+ raise ValueError("at least one hit is required")
157+ if n.dtype not in ("float32", "float64"):
158+ raise ValueError(f"numerics.dtype must be 'float32' or 'float64', got {n.dtype!r}")
159+ if not (0 < n.cfl_safety <= 1):
160+ raise ValueError("numerics.cfl_safety must lie in (0, 1]")
161+ if self.mesh.element_size <= 0 or self.mesh.element_size > p.radius / 2:
162+ raise ValueError("mesh.element_size must be positive and well below the radius")
163+ for i, h in enumerate(self.hits):
164+ if h.time < 0:
165+ raise ValueError(f"hit {i}: time must be non-negative")
166+ if h.velocity <= 0:
167+ raise ValueError(f"hit {i}: velocity must be positive (the exciter approaches from above)")
168+ ex = h.exciter
169+ extent = ex.patch_length / 2 if isinstance(ex, StickExciter) else ex.radius
170+ if h.r + extent >= p.radius - self.mesh.element_size:
171+ raise ValueError(f"hit {i}: contact patch reaches the rim (r = {h.r}, extent = {extent})")
172+ if ex.kind not in EXCITER_KINDS:
173+ raise ValueError(f"hit {i}: unknown exciter kind {ex.kind!r}")
174+
175+ # ---- dict / YAML
176+ def to_dict(self) -> dict:
177+ return _plain(dataclasses.asdict(self))
178+
179+ @classmethod
180+ def from_dict(cls, d: dict) -> "SimulationConfig":
181+ return _from_dict(cls, d)
182+
183+ def to_yaml(self, path: str | Path) -> None:
184+ Path(path).write_text(yaml.safe_dump(self.to_dict(), sort_keys=False))
185+
186+ @classmethod
187+ def from_yaml(cls, path: str | Path) -> "SimulationConfig":
188+ with open(path) as f:
189+ d = yaml.safe_load(f) or {}
190+ return cls.from_dict(d)
191+
192+
193+# --------------------------------------------------------------- generic dataclass <-> dict
194+def _plain(obj):
195+ """Recursively turn tuples into lists so that yaml.safe_dump accepts the structure."""
196+ if isinstance(obj, dict):
197+ return {k: _plain(v) for k, v in obj.items()}
198+ if isinstance(obj, (list, tuple)):
199+ return [_plain(v) for v in obj]
200+ return obj
201+
202+
203+def _build(tp, value):
204+ """Construct a value of annotated type ``tp`` from plain YAML data."""
205+ origin = typing.get_origin(tp)
206+ if isinstance(tp, type) and dataclasses.is_dataclass(tp):
207+ return _from_dict(tp, value)
208+ if origin is list:
209+ (item,) = typing.get_args(tp)
210+ return [_build(item, v) for v in value]
211+ if origin is tuple:
212+ return tuple(value)
213+ if origin in (typing.Union, types.UnionType):
214+ # The only union in the schema is the exciter, discriminated by its "kind" key.
215+ if isinstance(value, dict) and "kind" in value:
216+ try:
217+ return _from_dict(EXCITER_KINDS[value["kind"]], value)
218+ except KeyError:
219+ raise ValueError(f"unknown exciter kind {value['kind']!r}; "
220+ f"expected one of {sorted(EXCITER_KINDS)}") from None
221+ if any(isinstance(value, t) for t in typing.get_args(tp)):
222+ return value
223+ raise ValueError("an exciter mapping needs a 'kind' key ('stick' or 'mallet')")
224+ if tp is float and isinstance(value, (int, float)):
225+ return float(value)
226+ return value
227+
228+
229+def _from_dict(cls, d):
230+ if isinstance(d, cls):
231+ return d
232+ if not isinstance(d, dict):
233+ raise ValueError(f"expected a mapping for {cls.__name__}, got {type(d).__name__}")
234+ hints = typing.get_type_hints(cls)
235+ names = {f.name for f in dataclasses.fields(cls)}
236+ unknown = set(d) - names
237+ if unknown:
238+ raise ValueError(f"unknown field(s) for {cls.__name__}: {sorted(unknown)}")
239+ return cls(**{k: _build(hints[k], v) for k, v in d.items()})
240+
241+
242+def default_config(exciter: str = "mallet") -> SimulationConfig:
243+ """The prototype's default run for the given exciter kind."""
244+ if exciter == "mallet":
245+ hit = Hit(time=0.0, r=0.60 * Profile.radius, phi=0.3, velocity=8.0, exciter=MalletExciter())
246+ elif exciter == "stick":
247+ hit = Hit(time=0.0, r=0.85 * Profile.radius, phi=0.3, velocity=7.0, exciter=StickExciter())
248+ else:
249+ raise ValueError(f"unknown exciter kind {exciter!r}")
250+ return SimulationConfig(hits=[hit])
251+
252+
253+def t60_estimate(material: Material, frequency: float) -> float:
254+ """Decay time to -60 dB implied by sigma0 and sigma1 at a given frequency, for a flat plate
255+ with bending wavenumber k^2 = omega / sqrt(D / rho h). Rough, but useful for tuning losses."""
256+ # The sigma1 term acts as 2 sigma1 k^2; for a thin plate k^2 = omega sqrt(rho h / D) which
257+ # depends on thickness; we quote the value for h = 1 mm.
258+ h = 1e-3
259+ D = material.E * h**3 / (12 * (1 - material.nu**2))
260+ omega = 2 * math.pi * frequency
261+ k2 = omega * math.sqrt(material.rho * h / D)
262+ sigma = material.sigma0 + material.sigma1 * k2
263+ return 6.91 / sigma
src/cymsim/diagnostics.pyadded+96−0View file
@@ -0,0 +1,96 @@
1+"""Spectral diagnostics of the radiated pressure: per-octave band envelopes, onset times and the
2+spectral centroid, which together show whether energy climbs from low to high frequency (the
3+nonlinear upward cascade that gives a crash its build-up).
4+"""
5+from __future__ import annotations
6+
7+from dataclasses import dataclass
8+
9+import numpy as np
10+from scipy.signal import stft
11+
12+DEFAULT_BANDS = (60, 125, 250, 500, 1000, 2000, 4000, 8000, 16000)
13+
14+
15+@dataclass
16+class BandEnvelope:
17+ lo: float
18+ hi: float
19+ onset: float # first time the band reaches -20 dB of its own peak [s]
20+ peak_time: float # time of the band's peak [s]
21+ level_db: float # band peak relative to the loudest band [dB]
22+ envelope: np.ndarray # band energy vs time (linear)
23+
24+
25+@dataclass
26+class CascadeAnalysis:
27+ times: np.ndarray
28+ freqs: np.ndarray
29+ power: np.ndarray # STFT power (n_freq, n_time)
30+ bands: list[BandEnvelope]
31+ centroid: np.ndarray # spectral centroid vs time [Hz]
32+
33+ def centroid_at(self, t: float) -> float:
34+ return float(self.centroid[min(np.searchsorted(self.times, t), len(self.times) - 1)])
35+
36+
37+def analyze_cascade(p: np.ndarray, fs: int, bands=DEFAULT_BANDS, nperseg: int = 256,
38+ noverlap: int = 224) -> CascadeAnalysis:
39+ """Band envelopes at roughly 6 ms resolution (for fs = 44.1 kHz and the default hop)."""
40+ f_, t_, Z = stft(p, fs, nperseg=nperseg, noverlap=noverlap)
41+ Pw = np.abs(Z) ** 2
42+ tot_max = Pw.sum(0).max() or 1.0
43+ out = []
44+ for lo, hi in zip(bands[:-1], bands[1:]):
45+ m = (f_ >= lo) & (f_ < hi)
46+ env = Pw[m].sum(0)
47+ if env.max() <= 0:
48+ continue
49+ env_db = 10 * np.log10(env / env.max() + 1e-30)
50+ out.append(BandEnvelope(lo, hi, float(t_[np.argmax(env_db > -20)]), float(t_[np.argmax(env)]),
51+ float(10 * np.log10(env.max() / tot_max)), env))
52+ cen = (Pw * f_[:, None]).sum(0) / (Pw.sum(0) + 1e-30)
53+ return CascadeAnalysis(t_, f_, Pw, out, cen)
54+
55+
56+def format_cascade_report(a: CascadeAnalysis, label: str = "") -> str:
57+ lines = [f"[{label}] band onset(-20 dB) peak level vs loudest band"]
58+ for b in a.bands:
59+ lines.append(f" {b.lo:5.0f}-{b.hi:5.0f} Hz {b.onset * 1e3:7.1f} ms {b.peak_time * 1e3:7.1f} ms"
60+ f" {b.level_db:6.1f} dB")
61+ for tq in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5):
62+ if tq < a.times[-1]:
63+ lines.append(f" spectral centroid at {tq * 1e3:5.0f} ms: {a.centroid_at(tq):7.0f} Hz")
64+ return "\n".join(lines)
65+
66+
67+def plot_cascade(a: CascadeAnalysis, path: str, title: str = "", fs: float = 44100.0, t_max: float = 0.5) -> str:
68+ """Spectrogram and band envelopes (log-log). Requires matplotlib."""
69+ import matplotlib
70+ matplotlib.use("Agg")
71+ import matplotlib.pyplot as plt
72+
73+ tot_max = a.power.sum(0).max() or 1.0
74+ sel = a.times <= t_max
75+ fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
76+ ax[0].pcolormesh(a.times[sel] * 1e3, a.freqs, 10 * np.log10(a.power[:, sel] / a.power.max() + 1e-12),
77+ vmin=-90, vmax=0, shading="auto")
78+ ax[0].set_yscale("log")
79+ ax[0].set_ylim(50, fs / 2)
80+ ax[0].set_xlabel("t [ms]")
81+ ax[0].set_ylabel("f [Hz]")
82+ ax[0].set_title(title)
83+ for b in a.bands:
84+ ax[1].plot(a.times[sel] * 1e3, 10 * np.log10(b.envelope[sel] / tot_max + 1e-30), label=f"{b.lo:.0f}-{b.hi:.0f} Hz")
85+ ax[1].set_ylim(-90, 5)
86+ ax[1].set_xscale("log")
87+ ax[1].set_xlim(2, t_max * 1e3)
88+ ax[1].set_xlabel("t [ms]")
89+ ax[1].set_ylabel("band energy [dB re loudest]")
90+ ax[1].legend(fontsize=7, ncol=2)
91+ ax[1].grid(alpha=0.3)
92+ ax[1].set_title("band envelopes (a cascade shows staggered rises, high bands last)")
93+ plt.tight_layout()
94+ plt.savefig(path, dpi=120)
95+ plt.close(fig)
96+ return path
src/cymsim/elements.pyadded+99−0View file
@@ -0,0 +1,99 @@
1+"""Element formulations: the Discrete Kirchhoff Triangle for bending and the constant-strain
2+triangle for the Marguerre-von Karman membrane strains.
3+
4+Nodal degrees of freedom are (u, v, w, theta_x, theta_y) with theta_x = dw/dy and
5+theta_y = -dw/dx (the DKT convention of Batoz, Bathe and Ho, 1980). Bending enters through a
6+precomputed 9 x 9 element stiffness acting on (w, theta_x, theta_y) of the three nodes; the
7+membrane energy is evaluated from the constant gradients of u, v and w on the element together
8+with the slope of the initial surface. See docs/model.md, section "Spatial discretization".
9+"""
10+from __future__ import annotations
11+
12+import jax
13+import jax.numpy as jnp
14+
15+
16+def dkt_stiffness(xe, ye, D, nu):
17+ """9 x 9 DKT bending stiffness for one triangle.
18+
19+ ``xe``, ``ye``: nodal coordinates (3,), ``D`` the bending rigidity E h^3 / 12 (1 - nu^2).
20+ DOF order is (w1, thx1, thy1, w2, thx2, thy2, w3, thx3, thy3). The curvature operator is
21+ integrated with the three-point rule, which is exact for the quadratic rotation field.
22+ """
23+ i, j = jnp.array([1, 2, 0]), jnp.array([2, 0, 1]) # sides k = 4, 5, 6 <-> ij = 23, 31, 12
24+ xij, yij = xe[i] - xe[j], ye[i] - ye[j]
25+ l2 = xij ** 2 + yij ** 2
26+ a = -xij / l2
27+ b = 0.75 * xij * yij / l2
28+ c = (0.25 * xij ** 2 - 0.5 * yij ** 2) / l2
29+ d = -yij / l2
30+ e = (0.25 * yij ** 2 - 0.5 * xij ** 2) / l2
31+ a4, a5, a6 = a; b4, b5, b6 = b; c4, c5, c6 = c; d4, d5, d6 = d; e4, e5, e6 = e
32+
33+ def H(s):
34+ """Rotation shape functions (beta_x, beta_y) = H(xi, eta) . U, shape (2, 9)."""
35+ xi, eta = s
36+ lam = 1 - xi - eta
37+ N1, N2, N3 = lam * (2 * lam - 1), xi * (2 * xi - 1), eta * (2 * eta - 1)
38+ N4, N5, N6 = 4 * xi * eta, 4 * eta * lam, 4 * xi * lam
39+ Hx = jnp.stack([1.5 * (a6 * N6 - a5 * N5), b5 * N5 + b6 * N6, N1 - c5 * N5 - c6 * N6,
40+ 1.5 * (a4 * N4 - a6 * N6), b6 * N6 + b4 * N4, N2 - c4 * N4 - c6 * N6,
41+ 1.5 * (a5 * N5 - a4 * N4), b4 * N4 + b5 * N5, N3 - c4 * N4 - c5 * N5])
42+ Hy = jnp.stack([1.5 * (d6 * N6 - d5 * N5), -N1 + e5 * N5 + e6 * N6, -b5 * N5 - b6 * N6,
43+ 1.5 * (d4 * N4 - d6 * N6), -N2 + e4 * N4 + e6 * N6, -b6 * N6 - b4 * N4,
44+ 1.5 * (d5 * N5 - d4 * N4), -N3 + e4 * N4 + e5 * N5, -b4 * N4 - b5 * N5])
45+ return jnp.stack([Hx, Hy])
46+
47+ x21, x31 = xe[1] - xe[0], xe[2] - xe[0]
48+ y21, y31 = ye[1] - ye[0], ye[2] - ye[0]
49+ det = x21 * y31 - x31 * y21 # = 2 A
50+ xi_x, xi_y, eta_x, eta_y = y31 / det, -x31 / det, -y21 / det, x21 / det
51+
52+ def B_at(s):
53+ """Curvature operator (kappa_xx, kappa_yy, 2 kappa_xy) = B(s) . U, shape (3, 9)."""
54+ J = jax.jacfwd(H)(s) # (2, 9, 2): d(Hx, Hy)/d(xi, eta)
55+ Hx_xi, Hx_eta, Hy_xi, Hy_eta = J[0, :, 0], J[0, :, 1], J[1, :, 0], J[1, :, 1]
56+ Hx_x = Hx_xi * xi_x + Hx_eta * eta_x
57+ Hx_y = Hx_xi * xi_y + Hx_eta * eta_y
58+ Hy_x = Hy_xi * xi_x + Hy_eta * eta_x
59+ Hy_y = Hy_xi * xi_y + Hy_eta * eta_y
60+ return jnp.stack([Hx_x, Hy_y, Hx_y + Hy_x])
61+
62+ Db = D * jnp.array([[1, nu, 0], [nu, 1, 0], [0, 0, 0.5 * (1 - nu)]])
63+ gps = jnp.array([[1 / 6, 1 / 6], [2 / 3, 1 / 6], [1 / 6, 2 / 3]])
64+ Bs = jax.vmap(B_at)(gps) # (3, 3, 9)
65+ return (0.5 * det / 3) * jnp.einsum("gki,kl,glj->ij", Bs, Db, Bs)
66+
67+
68+def membrane_strains(qe, Gx, Gy, z0x, z0y, nl):
69+ """Marguerre-von Karman membrane strains (e_xx, e_yy, gamma_xy) on one element.
70+
71+ ``qe``: nodal (3, 5) values [u, v, w, thx, thy]; ``Gx``, ``Gy``: CST gradient operators (3,);
72+ ``z0x``, ``z0y``: slope of the initial surface at the element; ``nl``: 1.0 keeps the quadratic
73+ terms, 0.0 gives the linear shallow shell with the same geometry.
74+ """
75+ ux, uy = Gx @ qe[:, 0], Gy @ qe[:, 0]
76+ vx, vy = Gx @ qe[:, 1], Gy @ qe[:, 1]
77+ wx, wy = Gx @ qe[:, 2], Gy @ qe[:, 2]
78+ exx = ux + nl * 0.5 * wx * wx + z0x * wx
79+ eyy = vy + nl * 0.5 * wy * wy + z0y * wy
80+ gxy = uy + vx + nl * wx * wy + z0x * wy + z0y * wx
81+ return exx, eyy, gxy
82+
83+
84+def element_energy_parts(qe, Gx, Gy, A, z0x, z0y, Cm, Ke, nl):
85+ """(membrane, bending) strain energy of one element.
86+
87+ ``Cm`` = (C11, C12, C33) of the membrane stiffness E h / (1 - nu^2) [[1, nu, 0], [nu, 1, 0],
88+ [0, 0, (1 - nu)/2]]; ``Ke`` the 9 x 9 DKT stiffness.
89+ """
90+ exx, eyy, gxy = membrane_strains(qe, Gx, Gy, z0x, z0y, nl)
91+ C11, C12, C33 = Cm
92+ Em = 0.5 * A * (C11 * (exx * exx + eyy * eyy) + 2 * C12 * exx * eyy + C33 * gxy * gxy)
93+ U = qe[:, 2:].reshape(9)
94+ return Em, 0.5 * (U @ Ke @ U)
95+
96+
97+def element_energy(*args):
98+ Em, Eb = element_energy_parts(*args)
99+ return Em + Eb
src/cymsim/exciter.pyadded+152−0View file
@@ -0,0 +1,152 @@
1+"""Exciters, contact patches and the schedule of hits.
2+
3+Each hit is a rigid body of mass m that moves only along z, released at the hit time a small
4+clearance above the surface with downward speed v, and coupled to the shell through P springs
5+with the power law F = k delta^alpha, delta being the local penetration. A stick contacts along
6+a short line, a mallet on a small disk. Because the surface is curved, the rigid face is the
7+best-fit plane through the patch points and the residual heights give a staggered contact. All
8+hits of a run are stacked into fixed-shape arrays (padded to the largest patch) so that the time
9+loop handles any number of hits with the same compiled code. See docs/model.md, section
10+"Excitation".
11+"""
12+from __future__ import annotations
13+
14+import math
15+from dataclasses import dataclass
16+
17+import jax.numpy as jnp
18+import numpy as np
19+from scipy.integrate import quad
20+
21+from .config import Hit, MalletExciter, StickExciter
22+from .geometry import CymbalProfile
23+from .mesh import Mesh
24+
25+HOLD_HEIGHT = 1.0e3 # where an exciter waits before its release [m], far above any motion
26+
27+
28+@dataclass
29+class ContactPatch:
30+ """One hit's exciter, resolved on the mesh."""
31+ hit: Hit
32+ x: np.ndarray # (P,) patch point coordinates on the projected domain
33+ y: np.ndarray
34+ nodes: np.ndarray # (P, 3) containing-element nodes
35+ weights: np.ndarray # (P, 3) barycentric weights
36+ offsets: np.ndarray # (P,) surface height minus the exciter's best-fit face plane
37+ k_total: float # total contact stiffness [N/m^alpha]
38+ alpha: float
39+ mass: float
40+ gap: float
41+ velocity: float
42+ time: float
43+
44+ @property
45+ def n_points(self) -> int:
46+ return len(self.x)
47+
48+ @property
49+ def energy(self) -> float:
50+ return 0.5 * self.mass * self.velocity ** 2
51+
52+ def max_compression(self) -> float:
53+ """Compression at full stop against a rigid wall: m v^2 / 2 = k d^(a+1) / (a+1)."""
54+ a = self.alpha
55+ return ((a + 1) * self.energy / self.k_total) ** (1.0 / (a + 1))
56+
57+ def free_contact_time(self) -> float:
58+ """Contact duration against a rigid wall, 2 (d_max / v) int_0^1 dx / sqrt(1 - x^(a+1)).
59+ For alpha = 1.5 this is the Hertz value 2.94 d_max / v."""
60+ a = self.alpha
61+ c, _ = quad(lambda s: 1.0 / math.sqrt(1.0 - s ** (a + 1)), 0.0, 1.0)
62+ return 2.0 * c * self.max_compression() / self.velocity
63+
64+
65+def build_patch(hit: Hit, mesh: Mesh, profile: CymbalProfile) -> ContactPatch:
66+ ex = hit.exciter
67+ xs, ys = hit.r * math.cos(hit.phi), hit.r * math.sin(hit.phi)
68+ if isinstance(ex, StickExciter):
69+ ang = hit.phi + ex.patch_angle
70+ xi = np.linspace(-0.5, 0.5, ex.patch_points) * ex.patch_length
71+ px, py = xs + xi * math.cos(ang), ys + xi * math.sin(ang)
72+ k_total, alpha = ex.k, ex.alpha
73+ elif isinstance(ex, MalletExciter):
74+ alpha = ex.alpha
75+ # k such that the head stops after ``compression`` on a rigid wall
76+ k_total = (alpha + 1) * 0.5 * ex.mass * hit.velocity ** 2 / ex.compression ** (alpha + 1)
77+ th = np.arange(6) * np.pi / 3
78+ px = np.concatenate([[xs], xs + ex.radius * np.cos(th)])
79+ py = np.concatenate([[ys], ys + ex.radius * np.sin(th)])
80+ else:
81+ raise TypeError(f"unknown exciter type {type(ex).__name__}")
82+ loc = [mesh.locate(px[k], py[k]) for k in range(len(px))]
83+ nodes = np.array([l[0] for l in loc])
84+ weights = np.array([l[1] for l in loc])
85+ z0_p = profile.z0_at(np.stack([px, py], 1))
86+ # rigid exciter face: best-fit plane through the surface heights of the patch points
87+ A = np.column_stack([px - xs, py - ys, np.ones(len(px))])
88+ offsets = z0_p - A @ np.linalg.lstsq(A, z0_p, rcond=None)[0]
89+ return ContactPatch(hit, px, py, nodes, weights, offsets, float(k_total), float(alpha),
90+ float(ex.mass), float(ex.gap), float(hit.velocity), float(hit.time))
91+
92+
93+class HitSchedule:
94+ """All hits of a run as stacked device arrays (K hits, padded to P points each).
95+
96+ Padded points get zero weights and a very negative offset so they never make contact.
97+ """
98+
99+ def __init__(self, patches: list[ContactPatch], dt: float, dtype):
100+ self.patches = patches
101+ K = len(patches)
102+ P = max(p.n_points for p in patches)
103+ nodes = np.zeros((K, P, 3), dtype=np.int32)
104+ weights = np.zeros((K, P, 3))
105+ offsets = np.full((K, P), -1e6)
106+ for k, p in enumerate(patches):
107+ n = p.n_points
108+ nodes[k, :n] = p.nodes
109+ weights[k, :n] = p.weights
110+ offsets[k, :n] = p.offsets
111+ f = lambda a: jnp.asarray(a, dtype=dtype)
112+ self.n_hits, self.n_points = K, P
113+ self.nodes = jnp.asarray(nodes)
114+ self.weights, self.offsets = f(weights), f(offsets)
115+ self.k_point = f([p.k_total / p.n_points for p in patches])
116+ self.alpha = f([p.alpha for p in patches])
117+ self.mass = f([p.mass for p in patches])
118+ self.gap = f([p.gap for p in patches])
119+ self.velocity = f([p.velocity for p in patches])
120+ self.step = jnp.asarray([int(round(p.time / dt)) for p in patches], dtype=jnp.int32)
121+ self.dtype = dtype
122+
123+ def initial_state(self):
124+ """(face height, velocity) of every exciter before release."""
125+ return (jnp.full(self.n_hits, HOLD_HEIGHT, self.dtype), jnp.zeros(self.n_hits, self.dtype))
126+
127+ def contact(self, n, q, zs, vs):
128+ """Contact forces at step n.
129+
130+ Returns the per-point forces Fk (K, P), the total force per hit F (K,), and the exciter
131+ state (zs, vs) after applying the release of any hit whose time is step n.
132+ """
133+ wp = jnp.einsum("kpi,kpi->kp", self.weights, q[self.nodes, 2]) + self.offsets
134+ release = n == self.step
135+ active = n >= self.step
136+ zs = jnp.where(release, wp.max(1) + self.gap, zs)
137+ vs = jnp.where(release, -self.velocity, vs)
138+ delta = jnp.maximum(wp - zs[:, None], 0.0)
139+ Fk = jnp.where(active[:, None], self.k_point[:, None] * jnp.power(delta, self.alpha[:, None]), 0.0)
140+ return Fk, Fk.sum(1), zs, vs
141+
142+ def apply_to_plate(self, fo, Fk):
143+ """Scatter -Fk onto the w DOFs of the patch nodes."""
144+ return fo.at[self.nodes, 2].add(-Fk[:, :, None] * self.weights)
145+
146+ def advance(self, zs, vs, F, dt):
147+ """Symplectic Euler for the exciters (they are integrated with the plate's velocity update)."""
148+ vs1 = vs + dt * F / self.mass
149+ return zs + dt * vs1, vs1
150+
151+ def energy(self, vs):
152+ return 0.5 * (self.mass * vs * vs).sum()
src/cymsim/geometry.pyadded+116−0View file
@@ -0,0 +1,116 @@
1+"""Cymbal mid-surface height z0(x, y) and thickness h(x, y).
2+
3+The shell is described on the flat projected domain, the disk of radius R. The mid-surface height
4+z0 is the sum of a bow, a bell, a low-order rim warp and a set of hammering dents; the thickness
5+tapers linearly along the bow and steps up smoothly into the bell. All pieces are written with
6+jax.numpy so that the slope grad z0, which enters the membrane strains of the shallow shell, is
7+obtained by automatic differentiation rather than by finite differences. Imperfections are drawn
8+from a seeded generator so that a run is reproducible. See docs/model.md, section "Geometry".
9+"""
10+from __future__ import annotations
11+
12+import math
13+
14+import jax
15+import jax.numpy as jnp
16+import numpy as np
17+
18+from .config import Imperfections, Profile
19+
20+
21+class CymbalProfile:
22+ """Callable description of the cymbal surface, including the seeded imperfections.
23+
24+ ``z0(x, y)`` and ``thickness(x, y)`` take scalars and are differentiable; ``z0_at``,
25+ ``grad_z0_at`` and ``thickness_at`` take an (n, 2) array of points.
26+ """
27+
28+ def __init__(self, profile: Profile, imperfections: Imperfections):
29+ self.profile = profile
30+ self.imperfections = imperfections
31+ p, im = profile, imperfections
32+ rng = np.random.default_rng(im.seed)
33+
34+ # hammer dents: uniform over the annulus between the shoulder and the rim, random sign
35+ r0, r1 = p.bell_shoulder_radius + im.dent_margin, p.radius - im.dent_margin
36+ n = int(im.n_dents)
37+ rr = np.sqrt(rng.uniform(r0**2, r1**2, n))
38+ tt = rng.uniform(0.0, 2 * np.pi, n)
39+ self.dent_x = jnp.asarray(rr * np.cos(tt))
40+ self.dent_y = jnp.asarray(rr * np.sin(tt))
41+ self.dent_amp = jnp.asarray(im.dent_depth * rng.choice([-1.0, 1.0], n) * rng.uniform(0.5, 1.5, n))
42+ self.dent_sigma = jnp.asarray(im.dent_sigma * rng.uniform(0.7, 1.3, n))
43+
44+ # rim warp: azimuthal harmonics growing linearly from the shoulder to the rim
45+ modes = tuple(int(m) for m in im.rim_warp_modes)
46+ self.warp_m = jnp.asarray(modes, dtype=float)
47+ self.warp_a = jnp.asarray([im.rim_warp * 2.0 / m for m in modes], dtype=float)
48+ self.warp_phase = jnp.asarray(rng.uniform(0.0, 2 * np.pi, len(modes)))
49+
50+ # shoulder ovality and bell offset
51+ self.oval_phase = float(rng.uniform(0.0, 2 * np.pi))
52+ bo = float(rng.uniform(0.0, 2 * np.pi))
53+ self.bell_cx = im.bell_offset * math.cos(bo)
54+ self.bell_cy = im.bell_offset * math.sin(bo)
55+
56+ self._z0_v = jax.vmap(self.z0)
57+ self._h_v = jax.vmap(self.thickness)
58+ self._g_v = jax.vmap(jax.grad(lambda pt: self.z0(pt[0], pt[1])))
59+
60+ # ---- building blocks
61+ @staticmethod
62+ def _polar(x, y):
63+ r = jnp.sqrt(x * x + y * y + 1e-30)
64+ phi = jnp.arctan2(y, jnp.where(r < 1e-6, 1e-6, x)) # finite gradient at the centre
65+ return r, phi
66+
67+ def _shoulder_radius(self, phi):
68+ p, im = self.profile, self.imperfections
69+ return p.bell_shoulder_radius * (1.0 + im.shoulder_ovality * jnp.cos(2 * phi + self.oval_phase))
70+
71+ def _bell_blend(self, r_bell, r_sh):
72+ """1 on the bell plateau, 0 beyond the shoulder, C1 cosine ramp in between."""
73+ t = jnp.clip((r_bell - self.profile.bell_top_radius) / (r_sh - self.profile.bell_top_radius), 0.0, 1.0)
74+ return 0.5 * (1.0 + jnp.cos(jnp.pi * t))
75+
76+ def _bow_coordinate(self, x, y):
77+ """s in [0, 1]: 0 at the rim, 1 at the shoulder; plus the bell radius and shoulder radius."""
78+ p = self.profile
79+ r, phi = self._polar(x, y)
80+ r_sh = self._shoulder_radius(phi)
81+ r_bell = jnp.sqrt((x - self.bell_cx) ** 2 + (y - self.bell_cy) ** 2 + 1e-30)
82+ s = jnp.clip((p.radius - r) / (p.radius - r_sh), 0.0, 1.0)
83+ return s, phi, r_bell, r_sh
84+
85+ # ---- public scalar functions
86+ def z0(self, x, y):
87+ """Mid-surface height z0(x, y) [m]."""
88+ p = self.profile
89+ s, phi, r_bell, r_sh = self._bow_coordinate(x, y)
90+ u = jnp.maximum(1.0 - s, 1e-12) # 1 at the rim, 0 at the shoulder
91+ bow = p.bow_height * (1.0 - u ** p.bow_power)
92+ bell = p.bell_height * self._bell_blend(r_bell, r_sh)
93+ warp = u * jnp.sum(self.warp_a * jnp.cos(self.warp_m * phi + self.warp_phase))
94+ dents = jnp.sum(self.dent_amp * jnp.exp(
95+ -((x - self.dent_x) ** 2 + (y - self.dent_y) ** 2) / (2 * self.dent_sigma ** 2)))
96+ return bow + bell + warp + dents
97+
98+ def thickness(self, x, y):
99+ """Thickness h(x, y) [m]: linear taper along the bow, smooth step up into the bell."""
100+ p = self.profile
101+ s, phi, r_bell, r_sh = self._bow_coordinate(x, y)
102+ bow = p.h_rim + (p.h_shoulder - p.h_rim) * s
103+ return bow + (p.h_bell - p.h_shoulder) * self._bell_blend(r_bell, r_sh)
104+
105+ # ---- vectorized evaluation on (n, 2) point arrays
106+ def z0_at(self, xy) -> np.ndarray:
107+ xy = jnp.asarray(xy, dtype=jnp.float64)
108+ return np.asarray(self._z0_v(xy[:, 0], xy[:, 1]))
109+
110+ def thickness_at(self, xy) -> np.ndarray:
111+ xy = jnp.asarray(xy, dtype=jnp.float64)
112+ return np.asarray(self._h_v(xy[:, 0], xy[:, 1]))
113+
114+ def grad_z0_at(self, xy) -> np.ndarray:
115+ """Analytic slope (dz0/dx, dz0/dy) at each point, shape (n, 2)."""
116+ return np.asarray(self._g_v(jnp.asarray(xy, dtype=jnp.float64)))
src/cymsim/mesh.pyadded+150−0View file
@@ -0,0 +1,150 @@
1+"""Triangular mesh of the projected disk and the per-element geometry derived from it.
2+
3+The mesher needs only numpy. Nodes are placed on concentric rings with spacing equal to the
4+target element size; ring k carries round(2 pi r_k / dr) nodes, which is close to 6k, so the
5+triangles are nearly equilateral. Alternate rings are rotated by half a step to avoid radially
6+aligned edges. Consecutive rings are stitched by an angular merge: walking both rings with
7+increasing angle, the next triangle uses whichever ring's next node comes first, which yields a
8+valid triangulation of the annulus with (n_inner + n_outer) triangles. The centre is a single node
9+fanned to the first ring. See docs/model.md, section "Spatial discretization".
10+"""
11+from __future__ import annotations
12+
13+from dataclasses import dataclass, field
14+
15+import numpy as np
16+
17+
18+@dataclass
19+class Mesh:
20+ """Nodes ``xy`` (N, 2), counter-clockwise triangles ``tri`` (E, 3) and the rim node indices.
21+
22+ Derived quantities, computed once: element areas, the constant gradient operators Gx, Gy of
23+ the linear (CST) interpolation, and the centroids.
24+ """
25+ xy: np.ndarray
26+ tri: np.ndarray
27+ boundary: np.ndarray
28+ area: np.ndarray = field(init=False, repr=False)
29+ Gx: np.ndarray = field(init=False, repr=False)
30+ Gy: np.ndarray = field(init=False, repr=False)
31+ centroid: np.ndarray = field(init=False, repr=False)
32+
33+ def __post_init__(self):
34+ self.xy = np.asarray(self.xy, dtype=np.float64)
35+ self.tri = np.asarray(self.tri, dtype=np.int64)
36+ x, y = self.xy[self.tri, 0], self.xy[self.tri, 1]
37+ x21, x31 = x[:, 1] - x[:, 0], x[:, 2] - x[:, 0]
38+ y21, y31 = y[:, 1] - y[:, 0], y[:, 2] - y[:, 0]
39+ self.area = 0.5 * (x21 * y31 - x31 * y21)
40+ if not (self.area > 0).all():
41+ raise ValueError("mesh has non-positive element areas (orientation must be counter-clockwise)")
42+ # d(phi)/dx = sum_i Gx[i] phi_i for the linear interpolant on the triangle
43+ self.Gx = np.stack([y[:, 1] - y[:, 2], y[:, 2] - y[:, 0], y[:, 0] - y[:, 1]], 1) / (2 * self.area[:, None])
44+ self.Gy = np.stack([x[:, 2] - x[:, 1], x[:, 0] - x[:, 2], x[:, 1] - x[:, 0]], 1) / (2 * self.area[:, None])
45+ self.centroid = np.stack([x.mean(1), y.mean(1)], 1)
46+
47+ @property
48+ def n_nodes(self) -> int:
49+ return len(self.xy)
50+
51+ @property
52+ def n_elements(self) -> int:
53+ return len(self.tri)
54+
55+ def edge_lengths(self) -> np.ndarray:
56+ """All three edge lengths of every element, shape (E, 3)."""
57+ p = self.xy[self.tri]
58+ return np.linalg.norm(p - np.roll(p, -1, axis=1), axis=2)
59+
60+ def locate(self, px: float, py: float):
61+ """Containing element's nodes (3,) and barycentric weights (3,) of a point of the domain."""
62+ x, y = self.xy[self.tri, 0], self.xy[self.tri, 1]
63+ a2 = 2 * self.area
64+ l1 = ((y[:, 1] - y[:, 2]) * (px - x[:, 2]) + (x[:, 2] - x[:, 1]) * (py - y[:, 2])) / a2
65+ l2 = ((y[:, 2] - y[:, 0]) * (px - x[:, 2]) + (x[:, 0] - x[:, 2]) * (py - y[:, 2])) / a2
66+ l3 = 1 - l1 - l2
67+ inside = (l1 >= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9)
68+ if not inside.any():
69+ raise ValueError(f"point ({px:.4f}, {py:.4f}) lies outside the mesh")
70+ e = int(np.argmax(inside))
71+ return self.tri[e].copy(), np.array([l1[e], l2[e], l3[e]])
72+
73+ def interpolate(self, values: np.ndarray, px: float, py: float) -> float:
74+ nodes, w = self.locate(px, py)
75+ return float(values[nodes] @ w)
76+
77+
78+def _stitch(a_idx, a_ang, b_idx, b_ang):
79+ """Triangulate the annulus between an inner ring (a) and an outer ring (b).
80+
81+ Both rings are given as node indices and angles. Angles are unwrapped so that every outer
82+ angle lies in [a_ang[0], a_ang[0] + 2 pi); then the two rings are merged in increasing angle.
83+ """
84+ na, nb = len(a_idx), len(b_idx)
85+ b_ang = (b_ang - a_ang[0]) % (2 * np.pi) + a_ang[0]
86+ order = np.argsort(b_ang)
87+ b_idx, b_ang = b_idx[order], b_ang[order]
88+ ta = np.append(a_ang, a_ang[0] + 2 * np.pi)
89+ tb = np.append(b_ang, b_ang[0] + 2 * np.pi)
90+ out = []
91+ i = j = 0
92+ while i < na or j < nb:
93+ next_a = ta[i + 1] if i < na else np.inf
94+ next_b = tb[j + 1] if j < nb else np.inf
95+ if next_a <= next_b:
96+ out.append((a_idx[i % na], b_idx[j % nb], a_idx[(i + 1) % na]))
97+ i += 1
98+ else:
99+ out.append((a_idx[i % na], b_idx[j % nb], b_idx[(j + 1) % nb]))
100+ j += 1
101+ return out
102+
103+
104+def _smooth(xy: np.ndarray, tri: np.ndarray, fixed: np.ndarray, n_iter: int) -> np.ndarray:
105+ """Laplacian smoothing: move each free node to the mean of its neighbours."""
106+ N = len(xy)
107+ edges = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]])
108+ src = np.concatenate([edges[:, 0], edges[:, 1]])
109+ dst = np.concatenate([edges[:, 1], edges[:, 0]])
110+ deg = np.bincount(src, minlength=N)
111+ free = np.ones(N, bool)
112+ free[fixed] = False
113+ xy = xy.copy()
114+ for _ in range(n_iter):
115+ acc = np.zeros_like(xy)
116+ np.add.at(acc, src, xy[dst])
117+ mean = acc / deg[:, None]
118+ xy[free] = mean[free]
119+ return xy
120+
121+
122+def disk_mesh(radius: float, element_size: float, smoothing: int = 1) -> Mesh:
123+ """Near-equilateral triangulation of the disk of the given radius. The rim nodes stay on
124+ the circle; interior nodes are relaxed by ``smoothing`` Laplacian passes. One pass removes
125+ the worst stitching triangles; more passes improve the mean shape but shrink the smallest
126+ elements, which lowers the stable time step (see docs/model.md, section 4.1)."""
127+ n_rings = max(2, int(round(radius / element_size)))
128+ dr = radius / n_rings
129+ nodes = [(0.0, 0.0)]
130+ prev_idx, prev_ang = np.array([0]), np.array([0.0])
131+ tris: list[tuple[int, int, int]] = []
132+ for k in range(1, n_rings + 1):
133+ r = k * dr
134+ n_k = max(6, int(round(2 * np.pi * r / dr)))
135+ ang = 2 * np.pi * (np.arange(n_k) + 0.5 * (k % 2)) / n_k
136+ idx = len(nodes) + np.arange(n_k)
137+ nodes.extend(zip(r * np.cos(ang), r * np.sin(ang)))
138+ if k == 1:
139+ tris.extend((0, idx[j], idx[(j + 1) % n_k]) for j in range(n_k))
140+ else:
141+ tris.extend(_stitch(prev_idx, prev_ang, idx, ang))
142+ prev_idx, prev_ang = idx, ang
143+ xy = np.array(nodes, dtype=np.float64)
144+ tri = np.array(tris, dtype=np.int64)
145+ if smoothing > 0:
146+ xy = _smooth(xy, tri, prev_idx, smoothing)
147+ x, y = xy[tri, 0], xy[tri, 1]
148+ signed = (x[:, 1] - x[:, 0]) * (y[:, 2] - y[:, 0]) - (x[:, 2] - x[:, 0]) * (y[:, 1] - y[:, 0])
149+ tri[signed < 0] = tri[signed < 0][:, [0, 2, 1]]
150+ return Mesh(xy=xy, tri=tri, boundary=prev_idx)
src/cymsim/radiation.pyadded+67−0View file
@@ -0,0 +1,67 @@
1+"""Sound pressure at a listener point from the shell's transverse acceleration.
2+
3+The baffled plate is treated with the Rayleigh integral: every element is a point source of
4+volume acceleration A_e w_ddot_e at its centroid on the curved surface, delayed by the travel time
5+to the listener. An unbaffled cymbal also radiates from its back with opposite sign, and that
6+contribution reaches the listener by diffracting around the rim. We approximate it by the same
7+Rayleigh sum with the path length replaced by the shortest route over the rim (a Kirchhoff-type
8+edge-diffraction estimate, not a boundary-element solution) and subtract it from the front term.
9+Retarded times are handled at sub-step resolution with linear interpolation into two delay
10+buffers, which are decimated to the audio rate at the end. See docs/model.md, section
11+"Radiation".
12+"""
13+from __future__ import annotations
14+
15+import jax.numpy as jnp
16+import numpy as np
17+from scipy.signal import resample_poly
18+
19+from .config import Listener
20+from .geometry import CymbalProfile
21+from .mesh import Mesh
22+
23+
24+class RadiationModel:
25+ def __init__(self, mesh: Mesh, profile: CymbalProfile, z0_node: np.ndarray, listener: Listener,
26+ dt: float, n_steps: int, dtype, n_rim: int = 720):
27+ L = np.asarray(listener.position, dtype=float)
28+ self.listener = listener
29+ cent = np.column_stack([mesh.centroid, z0_node[mesh.tri].mean(1)]) # (E, 3) on the surface
30+ r_front = np.linalg.norm(L - cent, axis=1)
31+ # shortest path element -> rim point -> listener
32+ R = float(np.hypot(mesh.xy[:, 0], mesh.xy[:, 1]).max())
33+ phi = np.linspace(0, 2 * np.pi, n_rim, endpoint=False)
34+ rim_xy = np.column_stack([R * np.cos(phi), R * np.sin(phi)])
35+ rim = np.column_stack([rim_xy, profile.z0_at(rim_xy)])
36+ d1 = np.linalg.norm(cent[:, None, :] - rim[None, :, :], axis=2)
37+ d2 = np.linalg.norm(rim - L, axis=1)
38+ r_back = (d1 + d2[None, :]).min(1)
39+ self.r_front, self.r_back = r_front, r_back
40+
41+ c, rho = listener.c_air, listener.rho_air
42+ coef_f = rho * mesh.area / (2 * np.pi * r_front)
43+ coef_b = listener.back_diffraction_gain * rho * mesh.area / (2 * np.pi * r_back)
44+ df, db = r_front / (c * dt), r_back / (c * dt)
45+ i0f, i0b = np.floor(df).astype(np.int32), np.floor(db).astype(np.int32)
46+ self.buffer_length = int(n_steps + max(i0f.max(), i0b.max()) + 2)
47+ f = lambda a: jnp.asarray(a, dtype=dtype)
48+ self.coef_f, self.coef_b = f(coef_f), f(coef_b)
49+ self.frac_f, self.frac_b = f(df - np.floor(df)), f(db - np.floor(db))
50+ self.i0f, self.i0b = jnp.asarray(i0f), jnp.asarray(i0b)
51+ self.dtype = dtype
52+
53+ def initial_buffers(self):
54+ z = jnp.zeros(self.buffer_length, self.dtype)
55+ return z, z
56+
57+ def scatter(self, n, acc_e, bf, bb):
58+ """Add the element accelerations of step n into the retarded-time buffers."""
59+ sf, sb = self.coef_f * acc_e, self.coef_b * acc_e
60+ bf = bf.at[n + self.i0f].add(sf * (1 - self.frac_f)).at[n + self.i0f + 1].add(sf * self.frac_f)
61+ bb = bb.at[n + self.i0b].add(sb * (1 - self.frac_b)).at[n + self.i0b + 1].add(sb * self.frac_b)
62+ return bf, bb
63+
64+ @staticmethod
65+ def to_audio(buffer: np.ndarray, n_steps: int, substeps: int, n_out: int) -> np.ndarray:
66+ """Decimate a sub-step-rate buffer to the audio rate (anti-aliased polyphase resampling)."""
67+ return resample_poly(np.asarray(buffer[:n_steps], dtype=np.float64), 1, substeps)[:n_out]
src/cymsim/shell.pyadded+185−0View file
@@ -0,0 +1,185 @@
1+"""The discrete shell: element data, lumped mass, rigid-body modes, internal force and the
2+stable time-step bound.
3+
4+Given a mesh and a CymbalProfile, this module samples the thickness and the surface slope on the
5+elements, builds the DKT bending stiffness and the membrane stiffness of every element, lumps the
6+mass, and exposes the internal nodal force f_int(q) = -dU/dq obtained by reverse-mode automatic
7+differentiation of the summed element strain energies. It also provides the two operators the
8+integrator needs besides f_int: the sigma1 damping force on w and the projection that removes
9+the six rigid-body modes of the free shell. See docs/model.md, sections "Spatial discretization"
10+and "Time integration".
11+"""
12+from __future__ import annotations
13+
14+import math
15+from dataclasses import dataclass
16+
17+import jax
18+import jax.numpy as jnp
19+import numpy as np
20+
21+from .config import Material, Numerics
22+from .elements import dkt_stiffness, element_energy, element_energy_parts
23+from .geometry import CymbalProfile
24+from .mesh import Mesh
25+
26+N_DOF = 5 # (u, v, w, theta_x, theta_y) per node
27+
28+
29+@dataclass
30+class ElementChecks:
31+ """Residuals of the element self-tests; all should be small (1e-6 or better)."""
32+ rigid_tilt_residual: float # |K_e U_tilt| / |K_e| for a rigid tilt w = y, thx = 1
33+ constant_curvature_error: float # relative error of the energy of w = x^2 / 2
34+ zero_energy_modes_min: int # zero eigenvalues of M_e^-1 K_e per element (expect 6)
35+ zero_energy_modes_max: int
36+ rigid_mode_ratio: float # max |K r| / |K random| over the 6 rigid modes of the shell
37+
38+
39+class ShellModel:
40+ """Assembled shell on a mesh. All device arrays are in the simulation precision ``dtype``."""
41+
42+ def __init__(self, mesh: Mesh, profile: CymbalProfile, material: Material, numerics: Numerics):
43+ self.mesh = mesh
44+ self.material = material
45+ self.numerics = numerics
46+ self.dtype = jnp.float32 if numerics.dtype == "float32" else jnp.float64
47+ self.nl = 1.0 if numerics.nonlinear else 0.0
48+ xy, tri = mesh.xy, mesh.tri
49+ N, E = mesh.n_nodes, mesh.n_elements
50+ x, y = xy[tri, 0], xy[tri, 1]
51+ area = mesh.area
52+
53+ # ---- profile sampled on the mesh: thickness per element and node, height per node.
54+ # The slope that enters the membrane strain is the CST gradient of the nodal
55+ # heights, i.e. the initial surface is represented in the same piecewise-linear
56+ # space as w. This keeps a rigid tilt of the shell exactly force-free; the analytic
57+ # slope at the centroid would not (see docs/model.md).
58+ self.h_e = profile.thickness_at(mesh.centroid)
59+ self.h_node = profile.thickness_at(xy)
60+ self.z0_node = profile.z0_at(xy)
61+ z0_e = self.z0_node[tri]
62+ self.z0x = (mesh.Gx * z0_e).sum(1)
63+ self.z0y = (mesh.Gy * z0_e).sum(1)
64+ self.max_slope = float(np.hypot(self.z0x, self.z0y).max())
65+
66+ # ---- element stiffness: bending rigidity D, membrane stiffness C, DKT matrix
67+ E_, nu, rho = material.E, material.nu, material.rho
68+ self.D_e = E_ * self.h_e ** 3 / (12 * (1 - nu ** 2))
69+ self.Cm_e = (E_ * self.h_e / (1 - nu ** 2))[:, None] * np.array([1.0, nu, 0.5 * (1 - nu)])
70+ self.Ke = np.asarray(jax.vmap(dkt_stiffness, in_axes=(0, 0, 0, None))(
71+ jnp.asarray(x), jnp.asarray(y), jnp.asarray(self.D_e), nu))
72+
73+ # ---- lumped mass per DOF: a third of the element mass to each node; rotational inertia
74+ # I = m (h^2 / 12 + f A_e), the second term being the mass-scaling device discussed
75+ # in the documentation
76+ m_e = rho * self.h_e * area / 3
77+ I_e = m_e * (self.h_e ** 2 / 12 + numerics.rot_inertia_factor * area)
78+ M = np.zeros((N, N_DOF))
79+ for k in range(3):
80+ np.add.at(M[:, 0], tri[:, k], m_e * numerics.inplane_mass_scale)
81+ np.add.at(M[:, 1], tri[:, k], m_e * numerics.inplane_mass_scale)
82+ np.add.at(M[:, 2], tri[:, k], m_e)
83+ np.add.at(M[:, 3], tri[:, k], I_e)
84+ np.add.at(M[:, 4], tri[:, k], I_e)
85+ self.M = M
86+ self.m_e, self.I_e = m_e, I_e
87+ self.mass = float(M[:, 2].sum())
88+
89+ # ---- rigid-body modes of the curved shell, M-orthonormalised (see docs). A tilt drags
90+ # the in-plane displacements along: rotation about x by an angle a gives v = -a z0,
91+ # w = a y; rotation about y gives u = +a z0, w = -a x (first order, shallow shell).
92+ X, Y = xy[:, 0] - xy[:, 0].mean(), xy[:, 1] - xy[:, 1].mean()
93+ Rm = np.zeros((6, N, N_DOF))
94+ Rm[0, :, 0] = 1
95+ Rm[1, :, 1] = 1
96+ Rm[2, :, 0], Rm[2, :, 1] = -Y, X
97+ Rm[3, :, 2] = 1
98+ Rm[4, :, 2], Rm[4, :, 3], Rm[4, :, 1] = Y, 1, -self.z0_node
99+ Rm[5, :, 2], Rm[5, :, 4], Rm[5, :, 0] = -X, 1, self.z0_node
100+ for k in range(6):
101+ for l in range(k):
102+ Rm[k] -= (Rm[l] * M * Rm[k]).sum() * Rm[l]
103+ Rm[k] /= math.sqrt((Rm[k] * M * Rm[k]).sum())
104+ self.Rm = Rm
105+
106+ # ---- stability bound (Irons): lambda_max(M^-1 K) <= max_e lambda_max(M_e^-1 K_e), with
107+ # K_e the tangent stiffness at rest (Hessian of the element energy at q = 0)
108+ Khess = jax.vmap(jax.hessian(element_energy), in_axes=(0, 0, 0, 0, 0, 0, 0, 0, None))(
109+ jnp.zeros((E, 3, N_DOF)), jnp.asarray(mesh.Gx), jnp.asarray(mesh.Gy), jnp.asarray(area),
110+ jnp.asarray(self.z0x), jnp.asarray(self.z0y), jnp.asarray(self.Cm_e), jnp.asarray(self.Ke),
111+ self.nl).reshape(E, 15, 15)
112+ Me = np.tile(np.stack([m_e * numerics.inplane_mass_scale, m_e * numerics.inplane_mass_scale,
113+ m_e, I_e, I_e], 1), (1, 3)) # (E, 15), node-major
114+ S = np.asarray(Khess) / np.sqrt(Me[:, :, None] * Me[:, None, :])
115+ self._lam = np.linalg.eigvalsh(S) # (E, 15) ascending
116+ self.limiting_element = int(np.argmax(self._lam[:, -1]))
117+ self.omega_max = math.sqrt(self._lam[:, -1].max())
118+ self.dt_stable = 2.0 / self.omega_max
119+
120+ # ---- device arrays in simulation precision
121+ f = lambda a: jnp.asarray(a, dtype=self.dtype)
122+ self.tri_j = jnp.asarray(tri, dtype=jnp.int32)
123+ self.Gx_j, self.Gy_j, self.area_j = f(mesh.Gx), f(mesh.Gy), f(area)
124+ self.z0x_j, self.z0y_j, self.Cm_j, self.Ke_j = f(self.z0x), f(self.z0y), f(self.Cm_e), f(self.Ke)
125+ self.rhoh_area_j = f(rho * self.h_e * area)
126+ self.h_node_j = f(self.h_node)
127+ self.M_j, self.Minv_j, self.Rm_j = f(M), f(1.0 / M), f(Rm)
128+ self._sigma1 = float(material.sigma1)
129+ self._N = N
130+ self._f_int = jax.grad(lambda q: -self.total_energy(q))
131+
132+ # ---- energies and forces (plain JAX functions; the integrator jits the whole step)
133+ def energy_parts_per_element(self, q):
134+ """(membrane, bending) energy of every element for nodal state q (N, 5)."""
135+ return jax.vmap(element_energy_parts, in_axes=(0, 0, 0, 0, 0, 0, 0, 0, None))(
136+ q[self.tri_j], self.Gx_j, self.Gy_j, self.area_j, self.z0x_j, self.z0y_j,
137+ self.Cm_j, self.Ke_j, self.nl)
138+
139+ def total_energy(self, q):
140+ Em, Eb = self.energy_parts_per_element(q)
141+ return Em.sum() + Eb.sum()
142+
143+ def internal_force(self, q):
144+ """f_int = -dU/dq, shape (N, 5)."""
145+ return self._f_int(q)
146+
147+ def kinetic_energy(self, v):
148+ return 0.5 * (self.M_j * v * v).sum()
149+
150+ def damping_force_w(self, vw):
151+ """The sigma1 term as a nodal force on w: -2 sigma1 rho h Laplacian(w_dot) integrated
152+ with the CST gradients, i.e. 2 sigma1 sum_e (rho h A)_e G^T (G w_dot_e)."""
153+ ve = vw[self.tri_j]
154+ vx, vy = (self.Gx_j * ve).sum(1), (self.Gy_j * ve).sum(1)
155+ c = self.rhoh_area_j[:, None] * (self.Gx_j * vx[:, None] + self.Gy_j * vy[:, None])
156+ return 2 * self._sigma1 * jax.ops.segment_sum(c.ravel(), self.tri_j.ravel(), num_segments=self._N)
157+
158+ def project_rigid(self, v):
159+ """Remove the M-orthogonal projection onto the six rigid-body modes."""
160+ coef = jnp.einsum("knd,nd->k", self.Rm_j, self.M_j * v)
161+ return v - jnp.einsum("k,knd->nd", coef, self.Rm_j)
162+
163+ # ---- diagnostics
164+ def element_checks(self) -> ElementChecks:
165+ """Element self-tests (float64) and the rigid-mode test of the assembled tangent."""
166+ mesh = self.mesh
167+ x, y = mesh.xy[mesh.tri, 0], mesh.xy[mesh.tri, 1]
168+ E = mesh.n_elements
169+ U_tilt = np.stack([y, np.ones_like(y), np.zeros_like(y)], 2).reshape(E, 9)
170+ res = np.abs(np.einsum("eij,ej->ei", self.Ke, U_tilt)).max() / np.abs(self.Ke).max()
171+ U_cur = np.stack([0.5 * x ** 2, np.zeros_like(x), -x], 2).reshape(E, 9)
172+ Eb = 0.5 * np.einsum("ei,eij,ej->e", U_cur, self.Ke, U_cur)
173+ err = np.abs(Eb / (0.5 * self.D_e * mesh.area) - 1).max()
174+ n_zero = (self._lam < 1e-9 * self._lam[:, -1:]).sum(1)
175+
176+ q0 = jnp.zeros((mesh.n_nodes, N_DOF), self.dtype)
177+ K_apply = lambda d: jax.jvp(self._f_int, (q0,), (d,))[1]
178+ rng = np.random.default_rng(0)
179+ worst = 0.0
180+ for k in range(6):
181+ rnd = rng.standard_normal((mesh.n_nodes, N_DOF))
182+ rnd *= np.linalg.norm(self.Rm[k]) / np.linalg.norm(rnd)
183+ worst = max(worst, float(jnp.abs(K_apply(self.Rm_j[k])).max()
184+ / jnp.abs(K_apply(jnp.asarray(rnd, self.dtype))).max()))
185+ return ElementChecks(float(res), float(err), int(n_zero.min()), int(n_zero.max()), worst)
src/cymsim/simulate.pyadded+243−0View file
@@ -0,0 +1,243 @@
1+"""The time loop: leapfrog integration of the shell with contact, damping, rigid-mode projection
2+and radiation, block-compiled with JAX.
3+
4+One sub-step advances the state (q, v) of the shell, the state (zs, vs) of every exciter and the
5+two retarded-time pressure buffers:
6+
7+ f = f_int(q) - contact(q, zs) - 2 sigma1 rho h Lap(v_w)
8+ v' = P_rigid [ (v + dt M^-1 f) / (1 + 2 sigma0 dt) ]
9+ q' = q + dt v'
10+ vs' = vs + dt F / m, zs' = zs + dt vs'
11+
12+where P_rigid removes the six rigid-body modes. The sub-step rate is an integer multiple of the
13+audio rate chosen from the linear stability bound. See docs/model.md, section "Time integration".
14+"""
15+from __future__ import annotations
16+
17+import logging
18+import math
19+import time
20+from dataclasses import dataclass, field
21+
22+import jax
23+import jax.numpy as jnp
24+import numpy as np
25+from jax import lax
26+from scipy.io import wavfile
27+
28+from .config import SimulationConfig
29+from .exciter import HitSchedule, build_patch
30+from .geometry import CymbalProfile
31+from .mesh import disk_mesh
32+from .radiation import RadiationModel
33+from .shell import N_DOF, ShellModel
34+
35+log = logging.getLogger("cymsim")
36+
37+
38+@dataclass
39+class HitSummary:
40+ time: float
41+ kind: str
42+ peak_force: float # [N]
43+ contact_time: float # total time with F > 0 [s]
44+ rebound_velocity: float # exciter velocity at the end of the run [m/s], positive = upward
45+ energy_in: float # kinetic energy at release [J]
46+ energy_out: float # kinetic energy left in the exciter [J]
47+
48+
49+@dataclass
50+class EnergyRecord:
51+ time: float
52+ kinetic: float
53+ membrane: float
54+ bending: float
55+ exciters: float
56+ max_w_over_h: float # nonlinearity meter: < 1 effectively linear, 3 to 5 clear cascade
57+
58+ @property
59+ def total(self) -> float:
60+ return self.kinetic + self.membrane + self.bending + self.exciters
61+
62+
63+@dataclass
64+class SimulationResult:
65+ config: SimulationConfig
66+ fs: int
67+ dt: float
68+ substeps: int
69+ n_nodes: int
70+ n_elements: int
71+ pressure_rayleigh: np.ndarray # baffled Rayleigh integral at the listener [Pa]
72+ pressure_unbaffled: np.ndarray # front minus diffracted back term [Pa]
73+ hits: list[HitSummary]
74+ energy_log: list[EnergyRecord] = field(default_factory=list)
75+ wall_time: float = 0.0
76+
77+ @property
78+ def time(self) -> np.ndarray:
79+ return np.arange(len(self.pressure_rayleigh)) / self.fs
80+
81+ def write_wav(self, prefix: str, normalize: bool = True) -> list[str]:
82+ """Write <prefix>_rayleigh.wav and <prefix>_unbaffled.wav as 16-bit PCM. With
83+ ``normalize`` the peak is placed at 0.9 full scale; otherwise 1 Pa maps to 0.05 FS."""
84+ paths = []
85+ for name, p in (("rayleigh", self.pressure_rayleigh), ("unbaffled", self.pressure_unbaffled)):
86+ peak = float(np.abs(p).max()) or 1.0
87+ scale = 0.9 / peak if normalize else 0.05
88+ path = f"{prefix}_{name}.wav"
89+ wavfile.write(path, self.fs, (np.clip(p * scale, -1, 1) * 32767).astype(np.int16))
90+ paths.append(path)
91+ return paths
92+
93+ def write_npz(self, prefix: str) -> str:
94+ path = f"{prefix}_result.npz"
95+ np.savez(path, fs=self.fs, dt=self.dt, pressure_rayleigh=self.pressure_rayleigh,
96+ pressure_unbaffled=self.pressure_unbaffled,
97+ energy_time=[e.time for e in self.energy_log],
98+ energy=[[e.kinetic, e.membrane, e.bending, e.exciters, e.max_w_over_h] for e in self.energy_log])
99+ return path
100+
101+
102+class Simulator:
103+ """Build every part of the model from a SimulationConfig; ``run()`` integrates it."""
104+
105+ def __init__(self, config: SimulationConfig):
106+ config.validate()
107+ self.config = config
108+ num, mat = config.numerics, config.material
109+ t0 = time.time()
110+
111+ self.profile = CymbalProfile(config.profile, config.imperfections)
112+ self.mesh = disk_mesh(config.profile.radius, config.mesh.element_size, config.mesh.smoothing)
113+ log.info("mesh: %d nodes, %d triangles, %d DOFs (nonlinear = %s)",
114+ self.mesh.n_nodes, self.mesh.n_elements, N_DOF * self.mesh.n_nodes, num.nonlinear)
115+
116+ self.shell = ShellModel(self.mesh, self.profile, mat, num)
117+ sh = self.shell
118+ rim = sh.z0_node[self.mesh.boundary]
119+ log.info("profile: height %.1f mm, rim warp p-p %.2f mm, %d dents, h in [%.2f, %.2f] mm, "
120+ "max |grad z0| = %.3f, mass %.3f kg",
121+ sh.z0_node.max() * 1e3, (rim.max() - rim.min()) * 1e3, config.imperfections.n_dents,
122+ sh.h_e.min() * 1e3, sh.h_e.max() * 1e3, sh.max_slope, sh.mass)
123+ if sh.max_slope > config.profile.max_slope:
124+ log.warning("max slope %.2f exceeds the shallow-shell bound %.2f", sh.max_slope, config.profile.max_slope)
125+
126+ # ---- time step: integer number of sub-steps per audio sample, inside the stability bound
127+ self.fs = num.fs
128+ self.substeps = math.ceil(1.0 / (self.fs * num.cfl_safety * sh.dt_stable))
129+ self.dt = 1.0 / (self.fs * self.substeps)
130+ self.n_blocks = math.ceil(num.duration * self.fs * self.substeps / num.steps_per_block)
131+ self.n_steps = self.n_blocks * num.steps_per_block
132+ xc, yc = self.mesh.centroid[sh.limiting_element]
133+ log.info("omega_max bound %.3e rad/s (limiting element at r = %.3f m, h = %.2f mm); "
134+ "dt = %.3e s, %d sub-steps per sample, %d steps",
135+ sh.omega_max, math.hypot(xc, yc), sh.h_e[sh.limiting_element] * 1e3,
136+ self.dt, self.substeps, self.n_steps)
137+
138+ # ---- hits
139+ self.patches = [build_patch(h, self.mesh, self.profile) for h in config.hits]
140+ self.schedule = HitSchedule(self.patches, self.dt, sh.dtype)
141+ for k, p in enumerate(self.patches):
142+ tau = p.free_contact_time()
143+ log.info("hit %d at t = %.3f s: %s, %d-point patch at r = %.3f m, k = %.3g N/m^%.2f, "
144+ "m = %.3f kg, v = %.1f m/s, energy %.2f J, free contact ~ %.2f ms "
145+ "(direct excitation below ~%.0f Hz)",
146+ k, p.time, p.hit.exciter.kind, p.n_points, p.hit.r, p.k_total, p.alpha,
147+ p.mass, p.velocity, p.energy, tau * 1e3, 0.5 / tau)
148+ if p.time >= num.duration:
149+ log.warning("hit %d is scheduled after the end of the run", k)
150+
151+ # ---- radiation
152+ self.radiation = RadiationModel(self.mesh, self.profile, sh.z0_node, config.listener,
153+ self.dt, self.n_steps, sh.dtype)
154+ self._build_step()
155+ self.setup_time = time.time() - t0
156+ log.info("setup done in %.1f s; backend %s", self.setup_time, jax.devices()[0].platform)
157+
158+ # ---- one sub-step and its block loop
159+ def _build_step(self):
160+ sh, sched, rad = self.shell, self.schedule, self.radiation
161+ dt, sigma0 = self.dt, self.config.material.sigma0
162+ tri = sh.tri_j
163+
164+ def step(n, st):
165+ q, v, zs, vs, bf, bb, Fmax, ncont = st
166+ Fk, F, zs, vs = sched.contact(n, q, zs, vs)
167+ fo = sh.internal_force(q)
168+ fo = sched.apply_to_plate(fo, Fk)
169+ fo = fo.at[:, 2].add(-sh.damping_force_w(v[:, 2]))
170+ v1 = (v + dt * fo * sh.Minv_j) / (1 + 2 * sigma0 * dt) # sigma0 implicit
171+ v1 = sh.project_rigid(v1)
172+ q1 = q + dt * v1
173+ zs1, vs1 = sched.advance(zs, vs, F, dt)
174+ acc_e = ((v1[:, 2] - v[:, 2]) / dt)[tri].mean(1)
175+ bf, bb = rad.scatter(n, acc_e, bf, bb)
176+ return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0))
177+
178+ steps_per_block = self.config.numerics.steps_per_block
179+
180+ @jax.jit
181+ def run_block(n0, st):
182+ return lax.fori_loop(0, steps_per_block, lambda k, s: step((n0 + k).astype(jnp.int32), s), st)
183+
184+ @jax.jit
185+ def energies(st):
186+ q, v, _, vs = st[0], st[1], st[2], st[3]
187+ Em, Eb = sh.energy_parts_per_element(q)
188+ return (sh.kinetic_energy(v), Em.sum(), Eb.sum(), sched.energy(vs),
189+ jnp.abs(q[:, 2] / sh.h_node_j).max())
190+
191+ self._step, self._run_block, self._energies = step, run_block, energies
192+
193+ def initial_state(self):
194+ N, dtype = self.mesh.n_nodes, self.shell.dtype
195+ zs, vs = self.schedule.initial_state()
196+ bf, bb = self.radiation.initial_buffers()
197+ K = self.schedule.n_hits
198+ return (jnp.zeros((N, N_DOF), dtype), jnp.zeros((N, N_DOF), dtype), zs, vs, bf, bb,
199+ jnp.zeros(K, dtype), jnp.zeros(K, jnp.int32))
200+
201+ def run(self) -> SimulationResult:
202+ num = self.config.numerics
203+ state = self.initial_state()
204+ energy_log: list[EnergyRecord] = []
205+ t_run = time.time()
206+ report_every = max(1, self.n_blocks // max(1, num.n_progress))
207+ log.info("running %d blocks of %d steps (nonlinearity meter: max|w|/h < 1 effectively "
208+ "linear; 3 to 5 clear cascade)", self.n_blocks, num.steps_per_block)
209+ for blk in range(self.n_blocks):
210+ state = self._run_block(jnp.asarray(blk * num.steps_per_block, jnp.int32), state)
211+ if blk % report_every == 0 or blk == self.n_blocks - 1:
212+ ke, em, eb, ks, woh = (float(x) for x in self._energies(state))
213+ t = (blk + 1) * num.steps_per_block * self.dt
214+ energy_log.append(EnergyRecord(t, ke, em, eb, ks, woh))
215+ log.info(" t = %6.3f s KE %.3e membrane %.3e bending %.3e exciters %.3e J "
216+ "max|w|/h = %5.2f [%.0f s wall]", t, ke, em, eb, ks, woh, time.time() - t_run)
217+ if not math.isfinite(ke + em + eb):
218+ raise FloatingPointError("energy is not finite: the run has blown up "
219+ "(lower numerics.cfl_safety)")
220+
221+ q, v, zs, vs, bf, bb, Fmax, ncont = jax.device_get(state)
222+ hits = []
223+ for k, p in enumerate(self.patches):
224+ hits.append(HitSummary(p.time, p.hit.exciter.kind, float(Fmax[k]), float(ncont[k]) * self.dt,
225+ float(vs[k]), p.energy, 0.5 * p.mass * float(vs[k]) ** 2))
226+ log.info("hit %d (%s): peak force %.1f N, contact %.2f ms, rebound velocity %+.2f m/s "
227+ "(energy left in exciter %.2f of %.2f J)", k, p.hit.exciter.kind, hits[-1].peak_force,
228+ hits[-1].contact_time * 1e3, hits[-1].rebound_velocity, hits[-1].energy_out, p.energy)
229+
230+ n_out = int(self.fs * num.duration)
231+ p_front = self.radiation.to_audio(bf, self.n_steps, self.substeps, n_out)
232+ p_back = self.radiation.to_audio(bb, self.n_steps, self.substeps, n_out)
233+ result = SimulationResult(self.config, self.fs, self.dt, self.substeps, self.mesh.n_nodes,
234+ self.mesh.n_elements, p_front, p_front - p_back, hits, energy_log,
235+ time.time() - t_run)
236+ for name, p in (("rayleigh", result.pressure_rayleigh), ("unbaffled", result.pressure_unbaffled)):
237+ peak = float(np.abs(p).max())
238+ log.info("%-10s: peak %.3f Pa (%.1f dB SPL peak)", name, peak, 20 * np.log10(max(peak, 1e-12) / 2e-5))
239+ return result
240+
241+
242+def simulate(config: SimulationConfig) -> SimulationResult:
243+ return Simulator(config).run()
tests/conftest.pyadded+19−0View file
@@ -0,0 +1,19 @@
1+import logging
2+
3+import pytest
4+
5+from cymsim.config import SimulationConfig, default_config
6+
7+logging.getLogger("cymsim").setLevel(logging.WARNING)
8+
9+
10+@pytest.fixture(scope="session")
11+def coarse_config() -> SimulationConfig:
12+ """A 20 mm mesh with few dents: fast enough for CPU tests, still the full model."""
13+ cfg = default_config("mallet")
14+ cfg.mesh.element_size = 0.02
15+ cfg.imperfections.n_dents = 50
16+ cfg.numerics.duration = 0.02
17+ cfg.numerics.steps_per_block = 500
18+ cfg.numerics.n_progress = 2
19+ return cfg
tests/test_config.pyadded+43−0View file
@@ -0,0 +1,43 @@
1+import pytest
2+
3+from cymsim.config import (Hit, MalletExciter, SimulationConfig, StickExciter, default_config)
4+
5+
6+def test_yaml_round_trip(tmp_path):
7+ cfg = default_config("stick")
8+ cfg.hits.append(Hit(time=0.5, r=0.1, phi=1.0, velocity=3.0, exciter=MalletExciter(mass=0.2)))
9+ path = tmp_path / "c.yaml"
10+ cfg.to_yaml(path)
11+ back = SimulationConfig.from_yaml(path)
12+ assert back == cfg
13+ assert isinstance(back.hits[0].exciter, StickExciter)
14+ assert isinstance(back.hits[1].exciter, MalletExciter)
15+ assert back.imperfections.rim_warp_modes == (2, 3, 4, 5)
16+
17+
18+def test_partial_dict_uses_defaults():
19+ cfg = SimulationConfig.from_dict({"hits": [{"time": 0.1, "exciter": {"kind": "stick"}}],
20+ "numerics": {"duration": 0.5}})
21+ assert cfg.hits[0].exciter == StickExciter()
22+ assert cfg.numerics.duration == 0.5 and cfg.material.E == 1.1e11
23+
24+
25+def test_errors():
26+ with pytest.raises(ValueError):
27+ SimulationConfig.from_dict({"material": {"youngs": 1.0}})
28+ with pytest.raises(ValueError):
29+ SimulationConfig.from_dict({"hits": [{"exciter": {"kind": "hammer"}}]})
30+ cfg = default_config()
31+ cfg.hits[0].r = 0.23
32+ with pytest.raises(ValueError):
33+ cfg.validate()
34+ cfg = default_config()
35+ cfg.hits = []
36+ with pytest.raises(ValueError):
37+ cfg.validate()
38+
39+
40+def test_examples_load():
41+ import glob
42+ for path in glob.glob("examples/*.yaml"):
43+ SimulationConfig.from_yaml(path).validate()
tests/test_elements.pyadded+59−0View file
@@ -0,0 +1,59 @@
1+import jax
2+import jax.numpy as jnp
3+import numpy as np
4+
5+from cymsim.elements import dkt_stiffness, element_energy_parts, membrane_strains
6+
7+
8+def _triangle():
9+ return jnp.array([0.0, 0.03, 0.01]), jnp.array([0.0, 0.005, 0.025])
10+
11+
12+def test_dkt_rigid_modes_and_constant_curvature():
13+ xe, ye = _triangle()
14+ D, nu = 10.0, 0.3
15+ Ke = np.asarray(dkt_stiffness(xe, ye, D, nu))
16+ assert np.allclose(Ke, Ke.T)
17+ x, y = np.asarray(xe), np.asarray(ye)
18+ # rigid modes: translation, tilt about x (w = y, thx = 1), tilt about y (w = -x, thy = 1)
19+ for U in (np.tile([1, 0, 0], 3), np.stack([y, np.ones(3), np.zeros(3)], 1).ravel(),
20+ np.stack([-x, np.zeros(3), np.ones(3)], 1).ravel()):
21+ assert np.abs(Ke @ U).max() < 1e-12 * np.abs(Ke).max()
22+ # constant curvature w = x^2 / 2: energy D A / 2
23+ A = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]))
24+ U = np.stack([0.5 * x**2, np.zeros(3), -x], 1).ravel()
25+ assert np.isclose(0.5 * U @ Ke @ U, 0.5 * D * A, rtol=1e-9)
26+ # twist w = x y: energy D (1 - nu) A
27+ U = np.stack([x * y, x, -y], 1).ravel()
28+ assert np.isclose(0.5 * U @ Ke @ U, D * (1 - nu) * A, rtol=1e-9)
29+ eig = np.linalg.eigvalsh(Ke)
30+ assert (eig[:3] < 1e-10 * eig[-1]).all() and eig[3] > 1e-6 * eig[-1]
31+
32+
33+def test_membrane_strains_linear_and_nonlinear():
34+ Gx, Gy = jnp.array([1.0, -1.0, 0.0]), jnp.array([0.0, 1.0, -1.0])
35+ # nodal values arranged so that u_x = 2, v_y = 3, w_x = 0.1, w_y = 0.2 for these operators
36+ qe = jnp.zeros((3, 5)).at[:, 0].set(jnp.array([2.0, 0.0, 0.0])).at[:, 1].set(jnp.array([0.0, 3.0, 0.0])) \
37+ .at[:, 2].set(jnp.array([0.1, 0.0, -0.2]))
38+ exx, eyy, gxy = membrane_strains(qe, Gx, Gy, 0.5, -0.5, 1.0)
39+ assert np.isclose(exx, 2 + 0.5 * 0.01 + 0.5 * 0.1)
40+ assert np.isclose(eyy, 3 + 0.5 * 0.04 - 0.5 * 0.2)
41+ assert np.isclose(gxy, -3 + 0.1 * 0.2 + 0.5 * 0.2 - 0.5 * 0.1) # u_y + v_x = 0 - 3
42+ exx0, eyy0, _ = membrane_strains(qe, Gx, Gy, 0.5, -0.5, 0.0)
43+ assert np.isclose(exx0, 2 + 0.5 * 0.1) and np.isclose(eyy0, 3 - 0.5 * 0.2)
44+
45+
46+def test_energy_hessian_is_symmetric_positive_semidefinite():
47+ xe, ye = _triangle()
48+ Ke = dkt_stiffness(xe, ye, 10.0, 0.3)
49+ x, y = np.asarray(xe), np.asarray(ye)
50+ A = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]))
51+ Gx = jnp.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]]) / (2 * A)
52+ Gy = jnp.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]]) / (2 * A)
53+ Cm = jnp.array([1e6, 3e5, 3.5e5])
54+ f = lambda q: sum(element_energy_parts(q, Gx, Gy, A, 0.1, -0.05, Cm, Ke, 1.0))
55+ H = np.asarray(jax.hessian(f)(jnp.zeros((3, 5)))).reshape(15, 15)
56+ assert np.allclose(H, H.T)
57+ eig = np.linalg.eigvalsh(H)
58+ assert eig[0] > -1e-9 * eig[-1]
59+ assert (eig < 1e-9 * eig[-1]).sum() == 6 # 3 translations + 3 rotations
tests/test_mesh.pyadded+43−0View file
@@ -0,0 +1,43 @@
1+import numpy as np
2+import pytest
3+
4+from cymsim.mesh import disk_mesh
5+
6+
7+@pytest.mark.parametrize("h", [0.03, 0.012])
8+def test_disk_mesh_is_a_valid_manifold(h):
9+ R = 0.2286
10+ m = disk_mesh(R, h)
11+ assert (m.area > 0).all()
12+ # polygon area approaches pi R^2 from below
13+ assert 0.98 * np.pi * R**2 < m.area.sum() <= np.pi * R**2 + 1e-12
14+ edges = np.sort(np.concatenate([m.tri[:, [0, 1]], m.tri[:, [1, 2]], m.tri[:, [2, 0]]]), axis=1)
15+ _, count = np.unique(edges, axis=0, return_counts=True)
16+ assert set(count) == {1, 2} # interior edges shared by two triangles
17+ assert (count == 1).sum() == len(m.boundary) # boundary edges = rim nodes
18+ assert np.allclose(np.hypot(*m.xy[m.boundary].T), R)
19+ # element quality: inradius / circumradius, 1 for equilateral
20+ L = m.edge_lengths()
21+ a, b, c = L.T
22+ s = (a + b + c) / 2
23+ q = 8 * (s - a) * (s - b) * (s - c) / (a * b * c)
24+ assert q.min() > 0.5
25+ assert 0.5 * h < L.mean() < 1.5 * h
26+
27+
28+def test_cst_gradient_is_exact_for_linear_fields():
29+ m = disk_mesh(0.2, 0.03)
30+ f = 2.0 * m.xy[:, 0] - 3.0 * m.xy[:, 1] + 1.0
31+ fe = f[m.tri]
32+ assert np.allclose((m.Gx * fe).sum(1), 2.0)
33+ assert np.allclose((m.Gy * fe).sum(1), -3.0)
34+
35+
36+def test_locate_and_interpolate():
37+ m = disk_mesh(0.2, 0.03)
38+ f = m.xy[:, 0] ** 1 + 0.5 * m.xy[:, 1]
39+ nodes, w = m.locate(0.05, -0.07)
40+ assert np.isclose(w.sum(), 1.0) and (w >= -1e-9).all()
41+ assert np.isclose(m.interpolate(f, 0.05, -0.07), 0.05 - 0.035)
42+ with pytest.raises(ValueError):
43+ m.locate(0.3, 0.0)
tests/test_shell.pyadded+78−0View file
@@ -0,0 +1,78 @@
1+import numpy as np
2+import pytest
3+
4+from cymsim.config import Material, Numerics, Profile, Imperfections
5+from cymsim.geometry import CymbalProfile
6+from cymsim.mesh import disk_mesh
7+from cymsim.shell import ShellModel
8+
9+
10+@pytest.fixture(scope="module")
11+def shell():
12+ prof = CymbalProfile(Profile(), Imperfections(n_dents=30))
13+ mesh = disk_mesh(Profile.radius, 0.025)
14+ return ShellModel(mesh, prof, Material(), Numerics(dtype="float64"))
15+
16+
17+def test_mass_matches_integrated_density(shell):
18+ rho = shell.material.rho
19+ expected = rho * (shell.h_e * shell.mesh.area).sum()
20+ assert np.isclose(shell.mass, expected)
21+ assert (shell.M > 0).all()
22+
23+
24+def test_element_checks(shell):
25+ c = shell.element_checks()
26+ assert c.rigid_tilt_residual < 1e-12
27+ assert c.constant_curvature_error < 1e-8
28+ assert c.zero_energy_modes_min == 6 and c.zero_energy_modes_max == 6
29+ assert c.rigid_mode_ratio < 1e-10
30+
31+
32+def test_rigid_modes_are_M_orthonormal_and_force_free(shell):
33+ import jax
34+ import jax.numpy as jnp
35+ G = np.einsum("knd,nd,lnd->kl", shell.Rm, shell.M, shell.Rm)
36+ assert np.allclose(G, np.eye(6), atol=1e-10)
37+ # the tangent stiffness at rest annihilates every rigid mode
38+ q0 = jnp.zeros(shell.M.shape)
39+ K = lambda d: np.asarray(jax.jvp(shell.internal_force, (q0,), (jnp.asarray(d),))[1])
40+ rng = np.random.default_rng(0)
41+ for k in range(6):
42+ rnd = rng.standard_normal(shell.M.shape)
43+ rnd *= np.linalg.norm(shell.Rm[k]) / np.linalg.norm(rnd)
44+ assert np.abs(K(shell.Rm[k])).max() < 1e-12 * np.abs(K(rnd)).max()
45+ # a finite tilt is force-free only to first order: the residual scales with amplitude squared
46+ r1 = np.abs(np.asarray(shell.internal_force(shell.Rm_j[4] * 1e-3))).max()
47+ r2 = np.abs(np.asarray(shell.internal_force(shell.Rm_j[4] * 2e-3))).max()
48+ assert 3.5 < r2 / r1 < 4.5
49+
50+
51+def test_projection_removes_rigid_motion(shell):
52+ rng = np.random.default_rng(1)
53+ v = rng.standard_normal(shell.M.shape)
54+ p = np.asarray(shell.project_rigid(v))
55+ assert np.allclose(np.einsum("knd,nd->k", shell.Rm, shell.M * p), 0, atol=1e-9)
56+ assert np.allclose(np.asarray(shell.project_rigid(p)), p)
57+
58+
59+def test_internal_force_is_minus_energy_gradient(shell):
60+ rng = np.random.default_rng(2)
61+ q = 1e-4 * rng.standard_normal(shell.M.shape)
62+ d = rng.standard_normal(shell.M.shape)
63+ eps = 1e-6
64+ dE = (float(shell.total_energy(q + eps * d)) - float(shell.total_energy(q - eps * d))) / (2 * eps)
65+ assert np.isclose(-(np.asarray(shell.internal_force(q)) * d).sum(), dE, rtol=1e-5)
66+
67+
68+def test_damping_is_dissipative(shell):
69+ rng = np.random.default_rng(3)
70+ vw = rng.standard_normal(shell.mesh.n_nodes)
71+ fd = np.asarray(shell.damping_force_w(vw))
72+ assert (vw * fd).sum() > 0 # the force -fd does negative work on vw
73+ assert np.isclose(fd.sum(), 0.0, atol=1e-9 * np.abs(fd).sum()) # a uniform velocity is undamped
74+
75+
76+def test_stability_bound_is_finite_and_sane(shell):
77+ assert 1e4 < shell.omega_max < 1e8
78+ assert shell.dt_stable == 2.0 / shell.omega_max
tests/test_simulate.pyadded+59−0View file
@@ -0,0 +1,59 @@
1+import copy
2+
3+import numpy as np
4+import pytest
5+
6+from cymsim.config import Hit, MalletExciter, StickExciter
7+from cymsim.exciter import build_patch
8+from cymsim.simulate import Simulator
9+
10+
11+@pytest.fixture(scope="module")
12+def result(coarse_config):
13+ return Simulator(coarse_config).run()
14+
15+
16+def test_single_hit_run(result, coarse_config):
17+ n_out = int(coarse_config.numerics.fs * coarse_config.numerics.duration)
18+ assert result.pressure_rayleigh.shape == (n_out,)
19+ assert np.isfinite(result.pressure_unbaffled).all()
20+ assert np.abs(result.pressure_unbaffled).max() > 0.1 # a 9.6 J stroke is loud
21+ (h,) = result.hits
22+ assert h.peak_force > 100 and 0.5e-3 < h.contact_time < 20e-3
23+ assert 0 < h.rebound_velocity < coarse_config.hits[0].velocity
24+ assert h.energy_out < h.energy_in
25+ for e in result.energy_log:
26+ assert e.total <= h.energy_in * 1.01 # energy only leaves through damping
27+ assert result.energy_log[-1].max_w_over_h > 1 # this stroke is clearly nonlinear
28+
29+
30+def test_linear_reference_run_is_less_energetic_at_high_frequency(coarse_config):
31+ cfg = copy.deepcopy(coarse_config)
32+ cfg.numerics.nonlinear = False
33+ lin = Simulator(cfg).run()
34+ assert np.isfinite(lin.pressure_unbaffled).all()
35+
36+
37+def test_two_hits(coarse_config):
38+ cfg = copy.deepcopy(coarse_config)
39+ cfg.numerics.duration = 0.03
40+ cfg.hits = [Hit(time=0.0, r=0.137, phi=0.3, velocity=8.0, exciter=MalletExciter()),
41+ Hit(time=0.015, r=0.19, phi=2.0, velocity=7.0, exciter=StickExciter(patch_points=5))]
42+ sim = Simulator(cfg)
43+ assert sim.schedule.n_hits == 2 and sim.schedule.n_points == 7 # padded to the mallet's 7
44+ res = sim.run()
45+ assert len(res.hits) == 2
46+ assert res.hits[1].kind == "stick" and res.hits[1].peak_force > 0
47+ assert 0 < res.hits[1].contact_time < 5e-3
48+ assert res.hits[0].rebound_velocity > 0 and res.hits[1].rebound_velocity > 0
49+
50+
51+def test_contact_patch_geometry(coarse_config):
52+ sim = Simulator(coarse_config)
53+ stick = build_patch(Hit(r=0.19, phi=0.5, velocity=7.0, exciter=StickExciter()), sim.mesh, sim.profile)
54+ assert stick.n_points == 7 and np.isclose(np.hypot(stick.x[3], stick.y[3]), 0.19)
55+ assert np.allclose(stick.weights.sum(1), 1.0)
56+ assert np.isclose(stick.free_contact_time(), 2.94 * stick.max_compression() / 7.0, rtol=1e-2)
57+ mallet = build_patch(Hit(r=0.137, phi=0.3, velocity=8.0, exciter=MalletExciter()), sim.mesh, sim.profile)
58+ assert np.isclose(mallet.max_compression(), MalletExciter.compression)
59+ assert abs(mallet.offsets.max()) < 5e-3