cymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman model
37 changed files+4038−0
.gitignoreadded+10−0View file
@@ -0,0 +1,10 @@
1+__pycache__/
2+*.pyc
3+.venv/
4+dist/
5+build/
6+*.egg-info/
7+.pytest_cache/
8+*.wav
9+*.npz
10+.ipynb_checkpoints/
README.mdadded+101−0View file
@@ -0,0 +1,101 @@
1+# cymsim
2+
3+A Python package for simulating a drumstick strike on a thin metal plate or a cymbal. A plate is described in Python (material with its loss constants, bounding box, outline, thickness field, an optional stress-free profile) and sampled on a regular grid. The sound is computed in one of two ways. The **linear model** computes the modes of the plate once per design and synthesizes any strike from them in about a second, on the CPU. The **full model** solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly in time with PyTorch, on a GPU when one is available; it takes about a minute per second of sound on a GPU (several minutes on a CPU) and captures the pitch glide and the wash of a hard hit that the linear model cannot. Both use the same grid, so the linear sound is the small-amplitude limit of the full one, and a light hit sounds the same either way.
4+
5+The package is a port of the browser application [cymbalsim](../cymbalsim) to a library that can be scripted, so that designs, strikes and the numerics can be explored and compared programmatically. The model and the numerical method are described in [THEORY.md](THEORY.md).
6+
7+## Installation
8+
9+Python 3.10 or later. The linear model needs only NumPy and SciPy; the full model needs PyTorch, and the plotting helpers need Matplotlib.
10+
11+ pip install -e . # linear model only
12+ pip install -e ".[full]" # with PyTorch for the full model
13+ pip install -e ".[all]" # PyTorch, Matplotlib, IPython
14+
15+On Linux, `pip install torch` pulls the CUDA build; for a CPU-only machine `pip install torch --index-url https://download.pytorch.org/whl/cpu` is much smaller.
16+
17+## Quick start
18+
19+```python
20+import cymsim
21+
22+plate = cymsim.presets.crash_16().evaluate(spacing=0.008) # a 16 inch tapered bronze crash
23+print(plate.summary())
24+
25+modes = cymsim.compute_modes(plate, fmax=16000) # about 1900 modes, a few seconds
26+print(modes.table(10)) # frequencies and decay times, by mechanism
27+
28+strike = cymsim.Strike(x=0.9, y=0.5, velocity=3.0) # 3 m/s near the edge
29+listener = cymsim.Listener(distance=0.5, angle=30)
30+sound = cymsim.synthesize(modes, strike, listener, duration=4.0, movie=True)
31+sound.to_wav("crash_linear.wav") # peak-normalized 16-bit WAV
32+
33+full = cymsim.simulate(plate, strike, listener, duration=4.0) # the nonlinear model, on the GPU if present
34+full.to_wav("crash_full.wav")
35+```
36+
37+In a notebook, `sound.audio()` returns a playable widget, and `cymsim.plot` draws the plate, the mode shapes, the spectrogram (with the modal frequencies as ticks), the decay times and the movie frames:
38+
39+```python
40+from cymsim import plot
41+plot.plate(plate, strike=strike, listener=listener)
42+plot.mode(modes, 5)
43+plot.spectrogram(full, modes=modes)
44+plot.decay_times(modes)
45+```
46+
47+The `examples/` directory has scripts that run these steps.
48+
49+## Designing a plate
50+
51+A `Design` holds the material, the bounding box, and three fields: the outline (`inside`), the `thickness` and the stress-free `profile`. A field is a number or a function of four arrays, `(x, y, r, theta)`, the point relative to the center of the bounding box, its distance from the center and its angle from the x axis, evaluated with NumPy on the whole grid at once. Lengths are meters.
52+
53+```python
54+import numpy as np
55+from cymsim import Design, BRONZE
56+
57+inch, mm = 0.0254, 1e-3
58+R = 8 * inch
59+ramp = lambda a, b, x: np.clip((x - a) / (b - a), 0, 1)
60+
61+crash = Design(
62+ name="16 inch crash with a bow",
63+ material=BRONZE,
64+ width=2 * R, height=2 * R,
65+ inside=lambda x, y, r, th: r <= R,
66+ thickness=lambda x, y, r, th: 1.2 * mm + (0.7 * mm - 1.2 * mm) * ramp(40 * mm, R, r),
67+ profile=lambda x, y, r, th: 10 * mm * (1 - (r / R) ** 2),
68+)
69+plate = crash.evaluate(spacing=0.008)
70+```
71+
72+The `Material` carries the elastic constants, the loss factor of the metal's internal friction, and (optionally) the thermal constants from which the thermoelastic damping follows. The decay of the sound is not prescribed: it is derived from these constants and from the radiation of each mode into the air, which the linear model computes from the mode's shape. `cymsim.presets` has ten designs (flat and bowed steel sheets, an aluminum square, a bronze disc, a 16 inch crash, a 20 inch ride, the crash with a bow and a bell, the same hammered, a steel ellipse, a thunder sheet), each a function returning a `Design` that can be edited before evaluation.
73+
74+The profile is treated as a shallow shell, valid while its slopes are small compared with one; `plate.profile_stats` reports the largest slope, and `plate.summary()` warns above 0.1. The bell of a real cymbal has slopes of about 0.5, so the curved presets use a flatter bell than a real one.
75+
76+## The two models
77+
78+`compute_modes(plate, fmax)` assembles the linearized stiffness as a sparse matrix and computes all modes below `fmax` with a dense LAPACK solve restricted to that band (SciPy), or on a CUDA GPU with `torch.linalg.eigh` when one is available (`solver="torch"`). A flat plate has one unknown per node and its 16 inch crash at 8 mm takes about 4 s; a plate with a profile couples the transverse to the in-plane motion, has three unknowns per node, and takes about 40 s on the CPU. The result, a `Modes` object, holds the frequencies, the mass-normalized shapes on the grid, and the three loss rates of each mode; `modes.save` and `Modes.load` store it. `synthesize(modes, strike, listener, duration)` then produces a `Sound` in about a second; the strike and the listener can be changed at no further cost.
79+
80+`simulate(plate, strike, listener, duration, device=None, dtype="float32")` runs the full model. The time step is set by the in-plane wave speed (about 1.9 microseconds for bronze at 8 mm, so 4 s of sound is about 2 million steps). The step is written as a pure function of PyTorch tensors and compiled with `torch.compile`, which fuses its many small operations; on a CPU this gives about 90 s per second of sound for the 16 inch crash, and a GPU is expected to be several times faster (the code path is the same, but we have not timed it on a GPU). `Simulation` exposes the run batch by batch, with access to the state, for experiments. Single precision is the default, with the rounding safeguards described in THEORY.md; `dtype="float64"` is available.
81+
82+Both return a `Sound` with the listener pressure in pascals, the optional movie of the transverse displacement, and what the run reported (peak force, contact time). `Sound.to_wav` writes a 16-bit file at the peak, at a loudness target, or at an absolute level; `Sound.spectrogram` computes an STFT.
83+
84+## Tests
85+
86+ pytest # about 40 s on a 12-thread CPU
87+ pytest -m "not slow" # skips the two longer cases
88+
89+The tests are physics checks rather than comparisons with the browser application: the assembled stiffness has exactly the rigid modes in its null space and equals the linearized force of the full model's kernel to 1e-6 on a bowed, tapered plate with a hole; the kernel's forces are the exact gradient of the discrete energy at a strongly nonlinear state; the loss terms only remove energy; the free square and disc frequencies match Gorman's and Leissa's values as the grid is refined; Zener's loss factor peaks where it should; the per-mode radiation rate reproduces the closed form for a piston and Wallace's efficiencies for simply supported shapes, and the FFT convolution equals the direct double sum; a single mode under the stick follows the closed-form damped sinusoid; the attack and the tail of the synthesized signal join; the compiled step equals the uncompiled one; and a light hit simulated both ways gives spectra that agree.
90+
91+## Limitations
92+
93+The linear model has no pitch glide and no wash; a hard hit is only a loud one. The full model is the tool for that, and the two agree for a light hit. The von Karman equations hold for moderate rotations, and the Marguerre treatment of the profile for shallow shells, which strains at the bell of a real cymbal. The plate hangs freely: no stand, no thread, no air loading, and the radiation is that of a plate in a rigid baffle, which overstates the low-frequency level of a plate in free air. The full model carries the losses through a two-parameter fit to an average estimate, so its decay agrees with the modal one on average and not mode by mode. Grid dispersion lowers the high partials (by 3 percent near 1 kHz and 10 percent near 3 kHz for a millimeter of bronze at 8 mm), and the staircase edge places the low frequencies within a few percent. The residual stresses of hammering are absent. The dense eigensolver stores the full matrix, so a three-field system much beyond 10 000 unknowns (a curved 16 inch cymbal at 6 mm) becomes slow on the CPU; spectrum slicing, as the browser application does, would lift that limit.
94+
95+## Layout
96+
97+- `cymsim/material.py`, `grid.py`, `plate.py`, `design.py`, `params.py`, `presets.py`: materials, the grid and its mask, the plate on the grid, designs and their evaluation, the strike and the listener, the preset designs.
98+- `cymsim/loss.py`: the losses from physical constants (internal friction, Zener's thermoelastic damping, the average radiation efficiency, the outline perimeter) and the two-parameter fit of the full model.
99+- `cymsim/linear/`: the sparse assembly of K and M, the eigensolver and the `Modes` object, the per-mode radiation and thermoelastic rates, and the modal synthesis.
100+- `cymsim/full/`: the PyTorch kernel of one time step, the driver with the retarded-time listener sum, and the NumPy energy that the tests check the kernel against.
101+- `cymsim/sound.py`, `spectrogram.py`, `wav.py`, `plot.py`: the result object, the STFT, the WAV writer with the level mapping, and the Matplotlib views.
THEORY.mdadded+115−0View file
@@ -0,0 +1,115 @@
1+# Model and numerical method
2+
3+This note describes the model behind cymsim: a thin metal plate of arbitrary outline and spatially varying thickness, flat or with a stress-free curved profile, hanging freely, struck once with a drumstick. It follows the note of the browser application cymbalsim, of which cymsim is a port; the differences of implementation are noted where they occur (Sections 3 and 4). The plate is described by a `Design`; the strike and the listener are parameters of the setup. The sound is computed in one of two ways. The *linear model* takes the equations of the plate linearized about its rest shape, computes their modes once per design, and synthesizes the strike from the modes. The *full model* steps the geometrically nonlinear plate equations of Foppl and von Karman explicitly in time; it captures what the linear model cannot, the pitch glide and the transfer of energy toward high frequencies that turn the initial tone of a hard hit into a wash. Both use the same spatial discretization, so the linear model is exactly the small-amplitude limit of the full one.
4+
5+## 1. Plate model
6+
7+The plate occupies a region $\Omega$ of the plane, given by the design as a bounding box and an outline function. It has a thickness $h(x,y)$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$; the material constants are uniform, the thickness need not be. Its motion is described by the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$, all measured from a stress-free shape $w_0(x,y)$, the profile of the design, which is zero for a flat plate. The membrane strains in the Marguerre form are
8+
9+$$
10+\varepsilon_{xx} = u_x + \tfrac12\left[(w_0 + w)_x^2 - w_{0,x}^2\right], \qquad
11+\varepsilon_{yy} = v_y + \tfrac12\left[(w_0 + w)_y^2 - w_{0,y}^2\right], \qquad
12+\gamma_{xy} = u_y + v_x + (w_0 + w)_x (w_0 + w)_y - w_{0,x} w_{0,y},
13+$$
14+
15+so that $u = v = w = 0$ carries no stress. The membrane forces per unit length are $N_{xx} = C(\varepsilon_{xx} + \nu\varepsilon_{yy})$, $N_{yy} = C(\varepsilon_{yy} + \nu\varepsilon_{xx})$, and $N_{xy} = G\gamma_{xy}$, with $C = Eh/(1-\nu^2)$ and $G = Eh/(2(1+\nu))$. Writing $W = w_0 + w$ for the total shape and $D = Eh^3/(12(1-\nu^2))$ for the local bending stiffness, the equations of motion are
16+
17+$$
18+\rho h\, u_{tt} = \partial_x N_{xx} + \partial_y N_{xy}, \qquad
19+\rho h\, v_{tt} = \partial_x N_{xy} + \partial_y N_{yy},
20+$$
21+
22+$$
23+\rho h\, w_{tt} = -\mathcal{B}[w] + \partial_x(N_{xx} W_x + N_{xy} W_y) + \partial_y(N_{xy} W_x + N_{yy} W_y) + \text{(losses)} ,
24+$$
25+
26+where $\mathcal{B}[w]$ is the variable-stiffness bending operator, the Euler-Lagrange derivative of the bending energy
27+
28+$$
29+U_b = \frac12 \int_\Omega D\left[w_{xx}^2 + w_{yy}^2 + 2\nu\, w_{xx}w_{yy} + 2(1-\nu)\, w_{xy}^2\right] dA ,
30+$$
31+
32+which reduces to $D\nabla^4 w$ where $D$ is constant. For $w_0 = 0$ and constant $h$ these are the von Karman equations.
33+
34+**The linear part.** Dropping the terms of second order in the displacements, the strains become $\varepsilon_{xx} = u_x + w_{0,x} w_x$, $\varepsilon_{yy} = v_y + w_{0,y} w_y$, $\gamma_{xy} = u_y + v_x + w_{0,x} w_y + w_{0,y} w_x$. For a flat plate this decouples completely: the transverse motion is the Kirchhoff plate $\rho h w_{tt} = -\mathcal{B}[w]$, and the in-plane motion is neither excited by a transverse strike nor heard. A profile couples the two through its slopes: a curved plate stretches when it bends, which raises the low modes.
35+
36+**The nonlinear part.** The quadratic terms in the strains couple bending to stretching at finite amplitude. A flat plate has a symmetric, hardening nonlinearity; a curved plate has an asymmetric one, may snap between curved states, and softens as it passes through the flat configuration. These terms produce the pitch glide of a strongly struck cymbal and the cascade of energy from the low modes toward high frequencies. They are kept by the full model and dropped by the linear one.
37+
38+**Thickness field.** The thickness enters in three places: the mass per unit area $\rho h$, the membrane stiffnesses $C$ and $G$, proportional to $h$, and the bending stiffness $D$, proportional to $h^3$. A plate of twice the thickness has twice the bending frequencies and the same in-plane frequencies.
39+
40+**Profile.** Only the slopes of the stress-free shape enter the strains, so its mean is immaterial, but its tilt is not: the von Karman strains are not invariant under rotation, and a tilted stress-free shape would add a spurious membrane coupling linear in the tilt. The mass-weighted least-squares plane is therefore subtracted from every profile. The treatment is that of a shallow shell, valid while the slopes are small compared with one; `Plate.profile_stats` reports the largest slope.
41+
42+**Losses.** A partial of angular frequency $\omega$ decays as $e^{-\sigma t}$, and its decay time $T_{60} = \ln(1000)/\sigma$ is the time in which it falls by 60 dB; the rate $\sigma$ is the sum of three mechanisms. *Internal friction* of the metal is a hysteretic loss: each cycle dissipates the fraction $2\pi\eta$ of the stored energy, so a mode decays at $\sigma = \eta\omega/2$, and the loss factor $\eta$ is a material constant of order $10^{-4}$. *Thermoelastic damping* comes from the heat that flows across the thickness between the compressed and the stretched face of a bent plate; Zener's analysis gives the loss factor $\eta_{te}(\omega, h) = \Delta\,\omega\tau/(1 + \omega^2\tau^2)$ with the relaxation strength $\Delta = [E\alpha^2 T/(\rho c_p)]\,(1+\nu)/(1-\nu)$ and the relaxation time $\tau = h^2\rho c_p/(\pi^2 k)$, from the thermal expansion $\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$. *Radiation* into the air is the energy the plate loses as sound, which the linear model computes exactly for each mode from its shape (Section 3) and the full model takes from an analytic average (Section 4). A plate with free edges radiates poorly below its coincidence frequency (about 17 kHz for 1 mm bronze): its modes have no net volume velocity and no dipole moment, and at a free edge the bending moment and the shear vanish. On the presets the radiation overtakes the internal friction only above about 5 kHz.
43+
44+## 2. Spatial discretization
45+
46+The bounding box is covered by a regular grid of $n_x \times n_y$ nodes with spacings $h_x$ and $h_y$ (8 mm by default), and the outline selects the part of the grid that is the plate. All three displacements live at the nodes, and so does the thickness. The forces are the negative gradient of a discrete potential energy $U = U_b + U_m$, so the semi-discrete system is Hamiltonian, its linear part has a symmetric stiffness matrix, and the free edge needs no special formulas.
47+
48+**Outline.** A cell belongs to the plate when its center lies inside the outline, and a node when it is a corner of such a cell; the edge is a staircase of whole cells, and a hole is a region whose cells are absent. A second difference exists where the node and both of its neighbors belong to the plate, a cell contributes membrane energy and the mixed bending term when its four corners do, and the weight $\omega_k$ of a node is the fraction of its four cells that belong to the plate. Deciding membership by cell centers keeps the area unbiased.
49+
50+**Bending.** We evaluate $w_{xx}$ and $w_{yy}$ by three-point second differences at the nodes and $w_{xy}$ by the four-point difference at cell centers, and integrate the bending energy with the node weights, the bending stiffness $D_k$ at the nodes and $D_c$ at the cells (from the mean thickness of the four corners):
51+
52+$$
53+U_b = \frac{h_x h_y}{2}\left[\sum_k \omega_k D_k \left(a_k^2 + b_k^2 + 2\nu a_k b_k\right) + 2(1-\nu)\sum_c D_c\, c_c^2\right].
54+$$
55+
56+Where $D$ is constant the interior stencil is the standard 13-point biharmonic. The three rigid modes $w = 1, x, y$ lie in the null space for any thickness field and any outline.
57+
58+**Membrane.** The cells are bilinear (Q4) finite elements with 2 by 2 Gauss integration, using the strains of Section 1 with $w_0$ and $w$ interpolated bilinearly and the stiffnesses of the cell's mean thickness. One-point integration would be cheaper but admits hourglass modes, which the nonlinear coupling would excite.
59+
60+**Mass.** The mass matrix is lumped: the mass of node $k$ is $m_k = \rho h_k h_x h_y \omega_k$.
61+
62+**Accuracy.** For the rectangle the scheme converges at second order: the first seven frequencies of a completely free square plate are within 8 percent of Gorman's reference values at 10 cells per side and within 2 percent at 20. The staircase edge of a curved outline is only first-order accurate; the first nine frequencies of a free circular plate are within 3 percent of Leissa's values at 40 cells across, and the degenerate pairs, which the staircase splits, differ by about 1 percent. Grid dispersion affects every partial regardless of the outline: a plane bending wave of wavenumber $k$ along a grid axis has the discrete frequency ratio $(2 - 2\cos kh)/(kh)^2$, so a partial comes out low by 3 percent at $kh \approx 0.60$ and by 10 percent at $kh \approx 1.12$; for 1 mm bronze on the 8 mm grid these points are near 1.1 kHz and 3.6 kHz. `Plate.dispersion_frequency` reports them. The high partials are therefore audibly flattened at 8 mm, in both models alike, and a finer grid is a matter of patience rather than of the method.
63+
64+## 3. The linear model
65+
66+**Modes.** Collecting the linear part on the grid gives $M \ddot q + K q = f$, with $q$ the vector of unknowns (the transverse displacement of every plate node for a flat plate; $u$, $v$, $w$ interleaved for a plate with a profile), $M$ the diagonal lumped mass, and $K$ the symmetric stiffness, the Hessian of $U$ at the rest state; it is assembled as a sparse matrix from the stencil operators of Section 2, and the tests check that $Kq$ equals the linearized force of the full model to rounding. The modes are the solutions of $K\phi = \omega^2 M \phi$, normalized so that $\phi^T M \phi = 1$; the rigid modes at $\omega = 0$ are dropped, since a hanging plate's rigid motion is taken up by the thread and does not radiate in the model below.
67+
68+Here cymsim departs from the browser application, which slices the spectrum and solves each slice by shift-invert Lanczos in a pool of workers. With LAPACK at hand we form $A = M^{-1/2} K M^{-1/2}$ as a dense matrix and call the driver that computes the eigenpairs in a given interval, $(2\pi \cdot 1\,\mathrm{Hz})^2 < \lambda \le (2\pi f_{\max})^2$. A 16 inch cymbal at 8 mm has about 2150 nodes and 1900 bending modes below 16 kHz, computed in about 4 s; with a profile the same cymbal has 6500 unknowns and takes about 40 s on a 12-thread CPU, or a few seconds on a GPU with `torch.linalg.eigh`. The dense approach stores $n^2$ numbers and its cost grows as $n^3$, so it is the right tool up to roughly 10 000 unknowns; beyond that, spectrum slicing would be the remedy.
69+
70+**Modal equations.** Expanding $q = \sum_m \phi_m q_m(t)$ turns the plate into independent oscillators $\ddot q_m + 2\sigma_m \dot q_m + \omega_m^2 q_m = \phi_m^T f(t)$, with $\sigma_m$ the sum of the three loss rates of the mode. The internal friction gives $\eta\omega_m/2$; the thermoelastic loss gives $\sum_k \eta_{te}(\omega_m, h_k)\,U_{b,k}/\omega_m$, the local loss factor weighted by the bending energy $U_{b,k}$ of the mode at each node, which reduces to $\eta_{te}\omega_m/2$ on a flat uniform plate. The radiation rate is the power the mode radiates over twice its energy. A baffled plate whose nodes move with the velocities $v_k$ radiates the power $v^T R v$ from one face, with the radiation resistance matrix $R_{jk} = (\rho_{\mathrm{air}}\omega^2/4\pi c)\,\Delta A_j \Delta A_k\, \mathrm{sinc}(k_a r_{jk})$ of Elliott and Johnson, $k_a = \omega/c$; both faces radiate, so
71+
72+$$
73+\sigma_{\mathrm{rad},m} = \frac{\rho_{\mathrm{air}}\,\omega_m^2}{2\pi c} \sum_{j,k} \Delta A_j\,\Delta A_k\,\phi_m(j)\,\phi_m(k)\,\frac{\sin(k_a r_{jk})}{k_a r_{jk}},
74+$$
75+
76+which for a small piston at low frequency reduces to the textbook resistance $\rho_{\mathrm{air}}\omega^2 S^2/(2\pi c)$ divided by the mass. On a regular grid the kernel depends only on the offset $j - k$, so the double sum is a convolution, and by Parseval it is a sum over the spectrum: two 2-D FFTs of the zero-padded grid per mode, batched over the modes, rather than $N^2$ terms.
77+
78+**Strike.** The stick tip is a free mass $M_s$ (30 g by default) arriving with speed $v_0$ at a point of the plate. The contact patch consists of the plate nodes within the tip radius (5 mm) plus the nearest one; its contact coordinate is the mass-weighted mean transverse displacement $w_c = \sum_k (m_k/M_p) w_k$ of these nodes, the penetration is $\delta = w_c - z$ with $z$ the tip position, and the force is the Hertz law $F = K\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ and $K = 10^7\ \mathrm{N/m^{1.5}}$ by default, a softer contact than the Hertz value for a hard wooden tip on metal, which stands for the compliance of a real tip. The force is spread over the patch in proportion to the node masses, which makes the pair conjugate: the modal force is $-F g_m$ and the contact coordinate is $\sum_m g_m q_m$ with the same gains $g_m = \sum_k (m_k/M_p)\, \phi_m(k)$. The stick obeys $M_s \ddot z = F$. The contact is stepped at 2 microseconds; over each step the force is held constant, and every oscillator is advanced by its exact solution under a constant force, so the step is unconditionally stable for the plate and only the stick is integrated explicitly. Once the contact has been open for 5 ms the stick is withdrawn. After that instant $t_1$ every mode decays freely, and its state is one complex amplitude, $q_m(t) = \mathrm{Re}\left[A_m e^{s_m (t - t_1)}\right]$ with $s_m = -\sigma_m + i\sqrt{\omega_m^2 - \sigma_m^2}$.
79+
80+**Radiation.** The listener model is the Rayleigh integral for a baffled plate,
81+
82+$$
83+p(\mathbf r, t) = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \ddot w_k(t - \tau_k), \qquad \tau_k = r_k / c,
84+$$
85+
86+with $c = 343$ m/s; the common delay is discarded and only the spread of the delays across the plate is kept. For a mode in free decay the delay is a phase, so the whole sum over the plate collapses to one complex gain per mode,
87+
88+$$
89+G_m = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \phi_m(k)\, e^{-s_m \tau_k}, \qquad
90+p(t) = \sum_m \mathrm{Re}\left[G_m s_m^2 A_m e^{s_m (t - t_1)}\right],
91+$$
92+
93+which is exact for the free decay and includes the directivity of every mode. The signal is generated in blocks as a complex matrix product of the powers $e^{s_m n/f_s}$ with the coefficients. Before the free decay, while the modes are being driven by the stick, the delayed sum is evaluated directly from a table of the modal accelerations at twice the audio rate, with linear interpolation; the table path covers the output up to the moment at which every node's retarded time is past the last contact, about 2 ms after the impact, and the modal path takes over there. The tests check that the two join without a seam.
94+
95+**Movie.** The transverse displacement at the frame times is the modal sum $\sum_m \phi_m(k) q_m(t)$, one matrix product of the frame-by-mode amplitudes with the shapes.
96+
97+## 4. The full model
98+
99+The full model keeps the nonlinear strains of Section 1 and steps the semi-discrete system with the symplectic Euler scheme, $v^{n+1} = v^n + \Delta t\, a(x^n, v^n)$ and $x^{n+1} = x^n + \Delta t\, v^{n+1}$, in single precision by default. With the velocities interpreted at half steps this is the leapfrog scheme, second order and symplectic for the conservative part; the losses are evaluated explicitly at the old velocity. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system, computed at start-up as the largest eigenvalue of the three-field stiffness of a small uniform grid with the same spacing and the largest thickness of the design; 90 percent of the limit is used. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)}$, about 3800 m/s for bronze, set this limit: at 8 mm spacing $\Delta t \approx 1.9\ \mu$s, and 4 s of sound take about 2 million steps. The step is rounded down so that an integer number of steps makes one audio sample. The in-plane inertia terms are usually dropped in favor of an Airy stress function; we keep them because they make the whole system explicit and local, at the cost of the time step just described.
100+
101+**Losses.** The full model needs a loss operator rather than a rate per mode and uses the two-parameter model of plate synthesis: a bending mode of angular frequency $\omega$ and wavenumber $k$ decays at the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$. The two parameters are fitted by non-negative least squares to an estimate of the total rate on log-spaced frequencies from 50 Hz to 10 kHz, with $\kappa$ at the mean thickness. The estimate takes the internal friction as is, the thermoelastic loss factor averaged over the plate area with the local thickness, and the radiation from the average radiation efficiency of a baffled plate of the design's area and perimeter (the perimeter measured along the marching-squares contour of the cell mask), after Leppington, Broadbent and Heron, reduced by the factor $(f/f_c)^2$ that the free edges impose. The $\sigma_1$ term is implemented as $2\sigma_1$ times the graph Laplacian of the transverse velocity, whose eigenvalue on this grid is exactly $\omega/\kappa$ for a grid mode. The in-plane motion carries only the internal friction, as a frequency-independent rate at the lowest in-plane resonance of the plate.
102+
103+**Free plate.** The reaction to the stick is distributed over the plate as the inertial load of a rigid body: the force of the previous step is applied as the rigid acceleration field $a_k = F\,(1/M + y_k y_c/I_x + x_k x_c/I_y)$ that balances the patch load and its moments about the center, so the plate gains no net momentum or angular momentum. This is the counterpart of dropping the rigid modes in the linear model.
104+
105+**Single precision.** Three places in the kernel would otherwise turn rounding into audible noise. The first is the membrane strain of a curved design: written directly, $\varepsilon_{xx} = u_x + \tfrac12[(w_0+w)_x^2 - w_{0,x}^2]$ subtracts two nearly equal numbers once $w$ is small compared with the profile, and the rounding error of the two squares is proportional to $w_{0,x}^2$, a constant that does not fall as the vibration decays and acts as a broadband forcing. The strains are therefore written so that every term carries the deviation, $\varepsilon_{xx} = u_x + w_x(w_{0,x} + \tfrac12 w_x)$ with $w$ the deviation alone. The second is the update of $w$: with $w$ of a few millimeters and $\Delta t \approx 2\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, so the update is compensated, with the rounding error of each addition carried in a separate array and added back at the next step. The third is the listener sum, which forms the acceleration as the difference of successive delayed velocities at the audio rate rather than sampling the instantaneous acceleration. What remains is a noise floor proportional to the signal, roughly 75 dB below the peak in the top octave for a strongly curved design and 100 dB for a flat one; `dtype="float64"` removes it at a cost in speed on consumer GPUs.
106+
107+**The kernel.** Here cymsim departs from the browser application's hand-written WGSL kernel. The step is a pure function of PyTorch tensors over the $(n_y, n_x)$ grid: the second differences and the cell quantities are slices and pads of the whole arrays, the membrane forces of the four Gauss points are accumulated per cell corner and scattered to the nodes by four padded additions, the stick is a handful of scalar tensor operations, and the transverse velocity at the sample instants is written into a ring buffer from which the retarded-time listener sum is formed once per batch of one movie frame. As written the step is some two hundred small tensor operations, and on the CPU or a GPU alike it would be bound by their dispatch; `torch.compile` fuses them into a few kernels, which makes the step about six times faster on the CPU (0.17 ms for the 16 inch crash, about 90 s per second of sound) and is what makes a GPU run at all. The same code runs on CUDA, on Apple's MPS and on the CPU.
108+
109+## 5. Validation
110+
111+The test suite checks, on the CPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the kernel are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered; that the loss terms only remove energy; that the assembled $K$ equals the linearized force of the kernel to $10^{-6}$ on a bowed, tapered, holed plate; that the three eigensolvers agree and return mass-normalized shapes that satisfy the eigen equation; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that Zener's loss factor peaks at $\omega\tau = 1$ with the value $\Delta/2$, that the outline perimeter of a disc is within a few percent of $2\pi R$, that the average radiation efficiency has its textbook limits, that the radiation rate of a uniform piston at low frequency equals the closed form and the FFT convolution equals the direct double sum for every mode, that simply supported shapes reproduce Wallace's classical efficiencies to a few percent, and that the thermoelastic rate of a flat uniform plate is $\eta_{te}\omega/2$; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, and that the attack and the tail of the synthesized signal join; that a single-precision step of the kernel from a smooth random state of a few thicknesses amplitude agrees with the double-precision step to single-precision rounding, that the energy is conserved without losses, and that the compiled step equals the uncompiled one; and that a light hit (0.05 m/s) on a small plate simulated both ways gives magnitude spectra below 6 kHz with a correlation above 0.98 and levels within 25 percent, the remaining difference coming from the fitted two-parameter loss against the per-mode decay times and the modal truncation.
112+
113+## 6. Limitations
114+
115+The linear model is linear: it has no pitch glide, no wash, and a hard hit is only a loud one. The modal truncation at the cutoff frequency removes the highest partials, which the grid places poorly anyway. The von Karman model assumes moderate rotations, and the Marguerre treatment of the profile is a shallow-shell approximation, adequate for a bow of a few millimeters but strained by the bell of a real cymbal, whose slope is several times the range in which the approximation holds. The thickness field is treated as a slowly varying coefficient, appropriate for a lathed taper but not for a step on the scale of the grid. The hammered preset adds the shape of hammering but not its residual stresses. The loss factor of the metal and its thermal constants are typical published values rather than measurements of a particular alloy. The radiation loss is that of a plate in a rigid wall, an overstatement at low frequency for a plate in free air; the losses to the air's viscosity and to the support are absent; and the full model carries the losses only through the two-parameter fit to an average. Grid dispersion lowers the high partials by the amounts given in Section 2. The plate is free: neither a thread nor a stand is modeled, and air loading is not included. Single precision limits the dynamic range of the tail of the full model's sound to roughly 100 dB below the peak.
examples/full_model.pyadded+38−0View file
@@ -0,0 +1,38 @@
1+"""The full (nonlinear) model on the 16 inch crash with a bell, compared with the
2+linear model for the same strike. On a CPU this takes a few minutes per second
3+of sound; pass a shorter duration or a coarser spacing to try it out."""
4+
5+import sys
6+
7+import matplotlib.pyplot as plt
8+
9+import cymsim
10+from cymsim import plot
11+
12+duration = float(sys.argv[1]) if len(sys.argv) > 1 else 1.0
13+spacing = float(sys.argv[2]) if len(sys.argv) > 2 else 0.008
14+
15+plate = cymsim.presets.crash_bell_16().evaluate(spacing=spacing)
16+print(plate.summary())
17+strike = cymsim.Strike(x=0.9, y=0.5, velocity=3.0)
18+listener = cymsim.Listener(distance=0.5, angle=30)
19+
20+sim = cymsim.Simulation(plate, strike, listener, duration=duration, movie=True)
21+print(f"device {sim.device}, dt {sim.dt * 1e6:.2f} us, {sim.steps_per_sample} steps per sample, loss sigma0 {sim.loss[0]:.2f} 1/s, sigma1 {sim.loss[1]:.2e} m^2/s")
22+sim.run(progress=True)
23+full = sim.sound()
24+print(full, f"in {full.info['seconds']:.0f} s wall")
25+full.to_wav("crash_bell_full.wav")
26+
27+modes = cymsim.compute_modes(plate, fmax=16000)
28+linear = cymsim.synthesize(modes, strike, listener, duration=duration)
29+linear.to_wav("crash_bell_linear.wav")
30+
31+fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
32+plot.spectrogram(full, ax=axes[0], modes=modes, fmax=16000)
33+axes[0].set_title("full model")
34+plot.spectrogram(linear, ax=axes[1], modes=modes, fmax=16000)
35+axes[1].set_title("linear model")
36+fig.tight_layout()
37+fig.savefig("crash_bell.png", dpi=120)
38+print("wrote crash_bell_full.wav, crash_bell_linear.wav, crash_bell.png")
examples/quickstart.pyadded+28−0View file
@@ -0,0 +1,28 @@
1+"""The linear model on a 16 inch crash: modes, decay times, a strike, a WAV file and a few plots."""
2+
3+import matplotlib.pyplot as plt
4+
5+import cymsim
6+from cymsim import plot
7+
8+plate = cymsim.presets.crash_16().evaluate(spacing=0.008)
9+print(plate.summary())
10+
11+modes = cymsim.compute_modes(plate, fmax=16000)
12+print(f"{modes.count} modes in {modes.seconds:.1f} s")
13+print(modes.table(12))
14+
15+strike = cymsim.Strike(x=0.9, y=0.5, velocity=3.0)
16+listener = cymsim.Listener(distance=0.5, angle=30)
17+sound = cymsim.synthesize(modes, strike, listener, duration=4.0, movie=True)
18+print(sound, f"peak force {sound.info['peak_force']:.0f} N, contact {sound.info['contact_time'] * 1e3:.2f} ms")
19+sound.to_wav("crash_linear.wav")
20+
21+fig, axes = plt.subplots(2, 2, figsize=(12, 9))
22+plot.plate(plate, ax=axes[0, 0], strike=strike, listener=listener)
23+plot.mode(modes, 5, ax=axes[0, 1])
24+plot.spectrogram(sound, ax=axes[1, 0], modes=modes, fmax=16000)
25+plot.decay_times(modes, ax=axes[1, 1])
26+fig.tight_layout()
27+fig.savefig("crash_linear.png", dpi=120)
28+print("wrote crash_linear.wav and crash_linear.png")
examples/thickness_sweep.pyadded+18−0View file
@@ -0,0 +1,18 @@
1+"""An experiment the library makes easy: how the taper of a crash changes its
2+low modes and their decay. The preset is edited before evaluation."""
3+
4+import numpy as np
5+
6+import cymsim
7+from cymsim.presets import mm, ramp
8+
9+inch = 0.0254
10+R = 8 * inch
11+print(f"{'edge [mm]':>10} {'f0 [Hz]':>8} {'f10 [Hz]':>9} {'modes<16k':>10} {'T60 f0 [s]':>11} {'T60 5k [s]':>11}")
12+for edge in (0.5, 0.7, 0.9, 1.2):
13+ d = cymsim.presets.crash_16()
14+ d.name = f"crash, {edge} mm edge"
15+ d.thickness = lambda x, y, r, t, e=edge: 1.2 * mm + (e * mm - 1.2 * mm) * ramp(40 * mm, R, r)
16+ modes = cymsim.compute_modes(d.evaluate(0.010), fmax=16000)
17+ near5k = np.argmin(np.abs(modes.freqs - 5000))
18+ print(f"{edge:10.1f} {modes.freqs[0]:8.1f} {modes.freqs[10]:9.1f} {modes.count:10d} {modes.t60[0]:11.1f} {modes.t60[near5k]:11.2f}")
pyproject.tomladded+27−0View file
@@ -0,0 +1,27 @@
1+[build-system]
2+requires = ["hatchling"]
3+build-backend = "hatchling.build"
4+
5+[project]
6+name = "cymsim"
7+version = "0.1.0"
8+description = "Simulation of a drumstick strike on a plate or cymbal: modal synthesis of the linearized plate, and a geometrically nonlinear (von Karman) time-domain model on the GPU."
9+readme = "README.md"
10+requires-python = ">=3.10"
11+license = { text = "MIT" }
12+authors = [{ name = "Jeremy Magland" }]
13+dependencies = ["numpy>=1.24", "scipy>=1.10"]
14+
15+[project.optional-dependencies]
16+full = ["torch>=2.1"]
17+plot = ["matplotlib>=3.7"]
18+notebook = ["ipython"]
19+all = ["torch>=2.1", "matplotlib>=3.7", "ipython"]
20+dev = ["pytest>=7", "torch>=2.1", "matplotlib>=3.7"]
21+
22+[tool.hatch.build.targets.wheel]
23+packages = ["src/cymsim"]
24+
25+[tool.pytest.ini_options]
26+testpaths = ["tests"]
27+markers = ["slow: tests that take more than a few seconds"]
src/cymsim/__init__.pyadded+73−0View file
@@ -0,0 +1,73 @@
1+"""cymsim: simulation of a drumstick strike on a plate or cymbal.
2+
3+A plate is described by a ``Design`` (material, bounding box, outline,
4+thickness field, stress-free profile) and sampled on a grid into a ``Plate``.
5+Two models share that discretization. The *linear model* computes the modes
6+of the plate once (``compute_modes``) and synthesizes any strike from them in
7+about a second (``synthesize``). The *full model* steps the geometrically
8+nonlinear plate equations in time with PyTorch, on a GPU when one is available
9+(``simulate``); it is slower by orders of magnitude but captures the pitch
10+glide and the wash of a hard hit. Both return a ``Sound``.
11+
12+ import cymsim
13+ plate = cymsim.presets.crash_16().evaluate(spacing=0.008)
14+ modes = cymsim.compute_modes(plate, fmax=16000)
15+ sound = cymsim.synthesize(modes, cymsim.Strike(velocity=3.0), cymsim.Listener(distance=0.5))
16+ sound.to_wav("crash.wav")
17+"""
18+
19+from . import constants, loss, presets
20+from .design import Design, DesignError
21+from .grid import Grid, make_grid
22+from .linear import Modes, compute_modes, synthesize
23+from .material import ALUMINUM, BRONZE, MILD_STEEL, STAINLESS_STEEL, Material, Thermal
24+from .params import Listener, Strike
25+from .plate import Plate
26+from .sound import Sound
27+from .spectrogram import Spectrogram
28+
29+__version__ = "0.1.0"
30+
31+__all__ = [
32+ "Design",
33+ "DesignError",
34+ "Grid",
35+ "make_grid",
36+ "Material",
37+ "Thermal",
38+ "STAINLESS_STEEL",
39+ "MILD_STEEL",
40+ "BRONZE",
41+ "ALUMINUM",
42+ "Plate",
43+ "Strike",
44+ "Listener",
45+ "Modes",
46+ "compute_modes",
47+ "synthesize",
48+ "simulate",
49+ "Simulation",
50+ "Sound",
51+ "Spectrogram",
52+ "presets",
53+ "constants",
54+ "loss",
55+ "plot",
56+]
57+
58+
59+def simulate(*args, **kwargs):
60+ """The full model; see ``cymsim.full.simulate``. Requires ``torch`` (``pip install cymsim[full]``)."""
61+ from .full import simulate as _simulate
62+
63+ return _simulate(*args, **kwargs)
64+
65+
66+def __getattr__(name):
67+ import importlib
68+
69+ if name == "Simulation":
70+ return importlib.import_module(".full.simulation", __name__).Simulation
71+ if name == "plot":
72+ return importlib.import_module(".plot", __name__)
73+ raise AttributeError(f"module 'cymsim' has no attribute {name!r}")
src/cymsim/constants.pyadded+21−0View file
@@ -0,0 +1,21 @@
1+"""Physical and numerical constants shared by both models. All values are SI."""
2+
3+import math
4+
5+SPEED_OF_SOUND = 343.0
6+"""Speed of sound in air [m/s]."""
7+
8+RHO_AIR = 1.2
9+"""Density of air [kg/m^3]."""
10+
11+ROOM_TEMPERATURE = 293.0
12+"""Default plate temperature [K], which enters the thermoelastic damping."""
13+
14+CFL = 0.9
15+"""Fraction of the stability limit used for the time step of the full model."""
16+
17+STICK_RELEASE = 5e-3
18+"""The stick is withdrawn once the contact has been open this long [s]."""
19+
20+LN1000 = math.log(1000.0)
21+"""A decay rate sigma [1/s] and a decay time relate by sigma = LN1000 / T60."""
src/cymsim/design.pyadded+113−0View file
@@ -0,0 +1,113 @@
1+"""A design: the material, the bounding box, the outline, the thickness field and
2+the stress-free profile of a plate, as numbers or vectorized functions of position.
3+
4+The functions take four ``(ny, nx)`` arrays, ``(x, y, r, theta)``: the point
5+relative to the center of the bounding box, its distance from the center and
6+its angle from the x axis, all in meters and radians, and return an array of
7+the same shape (or a scalar). ``Design.evaluate`` samples them on a grid.
8+"""
9+
10+from __future__ import annotations
11+
12+from dataclasses import dataclass
13+from typing import Callable, Union
14+
15+import numpy as np
16+
17+from .constants import ROOM_TEMPERATURE
18+from .grid import InsideFn, make_grid
19+from .material import Material
20+from .plate import Plate, remove_plane
21+
22+Field = Union[float, Callable[[np.ndarray, np.ndarray, np.ndarray, np.ndarray], np.ndarray]]
23+"""A number, or a function ``f(x, y, r, theta)`` of grid matrices."""
24+
25+
26+class DesignError(ValueError):
27+ """A design that cannot be evaluated."""
28+
29+
30+@dataclass
31+class Design:
32+ """A plate design. Lengths are meters.
33+
34+ Example::
35+
36+ inch, mm = 0.0254, 1e-3
37+ R = 8 * inch
38+ crash = Design(
39+ name="16 inch flat crash",
40+ material=BRONZE,
41+ width=2 * R, height=2 * R,
42+ inside=lambda x, y, r, theta: r <= R,
43+ thickness=lambda x, y, r, theta: 1.2 * mm + (0.7 * mm - 1.2 * mm) * np.clip((r - 0.04) / (R - 0.04), 0, 1),
44+ )
45+ plate = crash.evaluate(spacing=0.008)
46+ """
47+
48+ name: str
49+ material: Material
50+ width: float
51+ """Bounding box [m]."""
52+ height: float
53+ thickness: Field
54+ """Thickness [m]: a number or a function of ``(x, y, r, theta)``; must be positive on the plate."""
55+ inside: InsideFn | None = None
56+ """Outline: ``inside(x, y, r, theta) -> bool``; ``None`` for the whole bounding box."""
57+ profile: Field = 0.0
58+ """Stress-free height [m]: 0 for a flat plate, or a function for a bow or a bell. Its mean plane is removed."""
59+ temperature: float = ROOM_TEMPERATURE
60+ """Plate temperature [K], for the thermoelastic damping."""
61+
62+ def __post_init__(self) -> None:
63+ if not (self.width >= 1e-3 and self.height >= 1e-3):
64+ raise DesignError("width and height must be at least 1 mm")
65+ if not callable(self.thickness) and not float(self.thickness) > 0:
66+ raise DesignError("thickness must be positive")
67+
68+ def evaluate(self, spacing: float = 0.008) -> Plate:
69+ """Samples the design on a grid of about the given spacing: the mask from
70+ the outline at the cell centers, the thickness and the profile at the
71+ nodes. Only the slopes of the profile enter the physics, so its mean is
72+ immaterial, but its tilt is not (the von Karman strains are not invariant
73+ under rotation), so the mass-weighted least-squares plane is removed."""
74+ g = make_grid(self.width, self.height, spacing, self.inside)
75+ if g.count == 0:
76+ raise DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing")
77+ X, Y = g.X, g.Y
78+ R, T = np.hypot(X, Y), np.arctan2(Y, X)
79+ h = sample_field(self.thickness, "thickness", X, Y, R, T)
80+ w0 = sample_field(self.profile, "profile", X, Y, R, T)
81+ m = g.mask
82+ bad = m & ~np.isfinite(h)
83+ if bad.any():
84+ j, i = np.argwhere(bad)[0]
85+ raise DesignError(f"thickness is not finite at ({X[j, i]:.4f}, {Y[j, i]:.4f})")
86+ bad = m & (h <= 0)
87+ if bad.any():
88+ j, i = np.argwhere(bad)[0]
89+ raise DesignError(f"thickness must be positive; it is {h[j, i]} at ({X[j, i]:.4f}, {Y[j, i]:.4f})")
90+ bad = m & ~np.isfinite(w0)
91+ if bad.any():
92+ j, i = np.argwhere(bad)[0]
93+ raise DesignError(f"profile is not finite at ({X[j, i]:.4f}, {Y[j, i]:.4f})")
94+ w0 = np.where(m, remove_plane(np.where(m, w0, 0.0), g), 0.0)
95+ return Plate(g, self.material, np.where(m, h, 0.0), w0, self.width, self.height, self.temperature, self.name)
96+
97+
98+def sample_field(f: Field, what: str, X: np.ndarray, Y: np.ndarray, R: np.ndarray, T: np.ndarray) -> np.ndarray:
99+ """Evaluates a field (a number or a function of ``(x, y, r, theta)``) on grid matrices."""
100+ if callable(f):
101+ try:
102+ v = np.asarray(f(X, Y, R, T), dtype=float)
103+ except Exception as e: # noqa: BLE001
104+ raise DesignError(f"{what}: {e}") from e
105+ if v.ndim == 0:
106+ return np.full(X.shape, float(v))
107+ if v.shape != X.shape:
108+ raise DesignError(f"{what} returned shape {v.shape} for a {X.shape[0]} x {X.shape[1]} grid; use elementwise numpy operations")
109+ return v
110+ try:
111+ return np.full(X.shape, float(f))
112+ except (TypeError, ValueError) as e:
113+ raise DesignError(f"{what} must be a number or a function of (x, y, r, theta)") from e
src/cymsim/full/__init__.pyadded+7−0View file
@@ -0,0 +1,7 @@
1+"""The full model: the geometrically nonlinear (von Karman) plate with in-plane
2+inertia and a thickness field, stepped explicitly in time with PyTorch, on a GPU
3+when one is available. Import requires ``torch`` (``pip install cymsim[full]``)."""
4+
5+from .simulation import Simulation, simulate
6+
7+__all__ = ["Simulation", "simulate"]
src/cymsim/full/kernel.pyadded+270−0View file
@@ -0,0 +1,270 @@
1+"""One explicit time step (symplectic Euler) of the von Karman plate with a
2+thickness field, a stress-free profile, the losses, and the drumstick, as a
3+pure function of PyTorch tensors over the ``(ny, nx)`` grid.
4+
5+The bending forces are the gradient of the discrete Kirchhoff energy with the
6+local bending stiffness D(h) at the nodes and at the cells (natural free-edge
7+conditions); the membrane forces come from bilinear cells with 2 x 2 Gauss
8+integration, the cell's mean thickness, and the Marguerre form of the von
9+Karman strains. Nodes outside the plate (weight 0) stay at zero: a second
10+difference exists only where the node and both neighbors belong to the plate,
11+and a cell contributes only when its four corners do. That gives the free edge
12+of the rectangle and a staircase free edge for any other outline, hole included.
13+
14+Three places would otherwise turn single-precision rounding into audible
15+noise, and are written as in the reference implementation: the membrane
16+strains carry the deviation from the stress-free shape in every term (the
17+direct form 0.5 (wx^2 - w0x^2) leaves a rounding error proportional to w0x^2,
18+which does not decay); the update of w is compensated, with the rounding error
19+carried in a separate array; and the listener sum (in ``simulation.py``) uses
20+the first difference of the delayed velocity rather than the instantaneous
21+acceleration.
22+"""
23+
24+from __future__ import annotations
25+
26+import math
27+from dataclasses import dataclass
28+
29+import numpy as np
30+import torch
31+import torch.nn.functional as F
32+
33+from ..constants import STICK_RELEASE
34+from ..linear.assemble import shape_derivatives
35+from ..params import Strike
36+from ..plate import Plate
37+
38+
39+def pad(t: torch.Tensor, left: int, right: int, top: int, bottom: int) -> torch.Tensor:
40+ return F.pad(t, (left, right, top, bottom))
41+
42+
43+@dataclass
44+class Static:
45+ """Per-node and per-cell data of a plate on the device, fixed for a run."""
46+
47+ ins: torch.Tensor
48+ """Node in the plate, as 0/1."""
49+ wgt: torch.Tensor
50+ inv_mass: torch.Tensor
51+ """1 / lumped mass, zero outside."""
52+ Dn: torch.Tensor
53+ """Bending stiffness at the nodes."""
54+ Dc: torch.Tensor
55+ """Bending stiffness at the cells (zero outside), ``(ny - 1, nx - 1)``."""
56+ Cm: torch.Tensor
57+ """Membrane stiffness at the cells times the Gauss weight hx hy / 4 (zero outside)."""
58+ Gs: torch.Tensor
59+ has_a: torch.Tensor
60+ has_b: torch.Tensor
61+ nb: tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]
62+ """Left, right, lower, upper neighbor in the plate (0/1)."""
63+ w0x: torch.Tensor
64+ """Slopes of the stress-free shape at the Gauss points, ``(4, ny - 1, nx - 1)``."""
65+ w0y: torch.Tensor
66+ dNdx: torch.Tensor
67+ """Shape function derivatives, ``(4 gp, 4 corners)``."""
68+ dNdy: torch.Tensor
69+ relief: torch.Tensor
70+ """Inertia-relief acceleration field per unit stick force."""
71+ patch_w: torch.Tensor
72+ """m_k / M_patch on the contact patch, zero elsewhere: the contact coordinate is its dot with w."""
73+ patch_f: torch.Tensor
74+ """1 / M_patch on the contact patch: the force distribution per unit mass."""
75+ hx: float
76+ hy: float
77+ nu: float
78+ dt: float
79+ sigma0: float
80+ sigma1: float
81+ sigma_in: float
82+ K: float
83+ alpha: float
84+ stick_mass: float
85+ strike: bool
86+
87+
88+def static_data(pl: Plate, strike: Strike | None, dt: float, sigma0: float, sigma1: float, sigma_in: float, device, dtype) -> Static:
89+ g, m = pl.grid, pl.material
90+ t = lambda a: torch.as_tensor(np.asarray(a, dtype=np.float64), dtype=dtype, device=device) # noqa: E731
91+ ins = g.mask.astype(float)
92+ mass = pl.node_mass
93+ inv_mass = np.where(mass > 0, 1 / np.where(mass > 0, mass, 1), 0.0)
94+ cm = g.cell_mask.astype(float)
95+ hc = pl.cell_thickness
96+ left, right, down, up = g.neighbor_masks()
97+ w0c = [np.asarray(c) for c in (pl.w0[:-1, :-1], pl.w0[:-1, 1:], pl.w0[1:, 1:], pl.w0[1:, :-1])]
98+ dNdx = np.zeros((4, 4))
99+ dNdy = np.zeros((4, 4))
100+ w0x = np.zeros((4,) + cm.shape)
101+ w0y = np.zeros((4,) + cm.shape)
102+ for gp in range(4):
103+ dNdx[gp], dNdy[gp] = shape_derivatives(gp, g.hx, g.hy)
104+ w0x[gp] = sum(dNdx[gp, b] * w0c[b] for b in range(4))
105+ w0y[gp] = sum(dNdy[gp, b] * w0c[b] for b in range(4))
106+ if strike is not None:
107+ patch = pl.strike_patch(strike)
108+ cx, cy = pl.center
109+ Ix, Iy = pl.inertia
110+ relief = ins * (1 / pl.mass + (g.Y - cy) * (patch.cy - cy) / Ix + (g.X - cx) * (patch.cx - cx) / Iy)
111+ patch_w = np.where(patch.mask, mass / patch.mass, 0.0)
112+ patch_f = np.where(patch.mask, 1 / patch.mass, 0.0)
113+ K, alpha, ms = strike.stiffness, strike.exponent, strike.mass
114+ else:
115+ relief = patch_w = patch_f = np.zeros(g.shape)
116+ K, alpha, ms = 0.0, 1.0, 1.0
117+ return Static(
118+ ins=t(ins),
119+ wgt=t(g.node_weight),
120+ inv_mass=t(inv_mass),
121+ Dn=t(m.bending_stiffness(pl.h)),
122+ Dc=t(m.bending_stiffness(hc) * cm),
123+ Cm=t(g.hx * g.hy / 4 * m.membrane_stiffness(hc) * cm),
124+ Gs=t(g.hx * g.hy / 4 * m.shear_stiffness(hc) * cm),
125+ has_a=t(left & right),
126+ has_b=t(down & up),
127+ nb=(t(left), t(right), t(down), t(up)),
128+ w0x=t(w0x),
129+ w0y=t(w0y),
130+ dNdx=t(dNdx),
131+ dNdy=t(dNdy),
132+ relief=t(relief),
133+ patch_w=t(patch_w),
134+ patch_f=t(patch_f),
135+ hx=g.hx,
136+ hy=g.hy,
137+ nu=m.nu,
138+ dt=dt,
139+ sigma0=sigma0,
140+ sigma1=sigma1,
141+ sigma_in=sigma_in,
142+ K=K,
143+ alpha=alpha,
144+ stick_mass=ms,
145+ strike=strike is not None,
146+ )
147+
148+
149+def accelerations(S: Static, u: torch.Tensor, v: torch.Tensor, w: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
150+ """Undamped accelerations of every node for the displacements ``(u, v, w)`` from the stress-free shape."""
151+ hx, hy, nu = S.hx, S.hy, S.nu
152+ ihx2, ihy2, ihxy = 1 / hx**2, 1 / hy**2, 1 / (hx * hy)
153+
154+ # Bending: weighted moments at the nodes, mixed differences at the cells.
155+ wp = pad(w, 1, 1, 1, 1)
156+ a = S.has_a * (wp[1:-1, 2:] - 2 * w + wp[1:-1, :-2]) * ihx2
157+ b = S.has_b * (wp[2:, 1:-1] - 2 * w + wp[:-2, 1:-1]) * ihy2
158+ wD = S.wgt * S.Dn
159+ Ma = S.has_a * wD * (a + nu * b)
160+ Mb = S.has_b * wD * (b + nu * a)
161+ Cc = S.Dc * (w[1:, 1:] - w[1:, :-1] - w[:-1, 1:] + w[:-1, :-1]) * ihxy
162+ Map = pad(Ma, 1, 1, 0, 0)
163+ Mbp = pad(Mb, 0, 0, 1, 1)
164+ Cp = pad(Cc, 1, 1, 1, 1)
165+ bend = (Map[:, 2:] - 2 * Ma + Map[:, :-2]) * ihx2 + (Mbp[2:, :] - 2 * Mb + Mbp[:-2, :]) * ihy2 + 2 * (1 - nu) * (Cp[1:, 1:] - Cp[1:, :-1] - Cp[:-1, 1:] + Cp[:-1, :-1]) * ihxy
166+ Fw = -hx * hy * bend
167+
168+ # Membrane: bilinear cells, 2 x 2 Gauss points, Marguerre strains with the deviation in every term.
169+ U = (u[:-1, :-1], u[:-1, 1:], u[1:, 1:], u[1:, :-1])
170+ V = (v[:-1, :-1], v[:-1, 1:], v[1:, 1:], v[1:, :-1])
171+ W = (w[:-1, :-1], w[:-1, 1:], w[1:, 1:], w[1:, :-1])
172+ Fu_c = [None] * 4
173+ Fv_c = [None] * 4
174+ Fw_c = [None] * 4
175+ for gp in range(4):
176+ dNdx, dNdy = S.dNdx[gp], S.dNdy[gp]
177+ ux = dNdx[0] * U[0] + dNdx[1] * U[1] + dNdx[2] * U[2] + dNdx[3] * U[3]
178+ uy = dNdy[0] * U[0] + dNdy[1] * U[1] + dNdy[2] * U[2] + dNdy[3] * U[3]
179+ vx = dNdx[0] * V[0] + dNdx[1] * V[1] + dNdx[2] * V[2] + dNdx[3] * V[3]
180+ vy = dNdy[0] * V[0] + dNdy[1] * V[1] + dNdy[2] * V[2] + dNdy[3] * V[3]
181+ dwx = dNdx[0] * W[0] + dNdx[1] * W[1] + dNdx[2] * W[2] + dNdx[3] * W[3]
182+ dwy = dNdy[0] * W[0] + dNdy[1] * W[1] + dNdy[2] * W[2] + dNdy[3] * W[3]
183+ w0x, w0y = S.w0x[gp], S.w0y[gp]
184+ wx, wy = w0x + dwx, w0y + dwy
185+ exx = ux + dwx * (w0x + 0.5 * dwx)
186+ eyy = vy + dwy * (w0y + 0.5 * dwy)
187+ gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy
188+ Nxx = S.Cm * (exx + nu * eyy)
189+ Nyy = S.Cm * (eyy + nu * exx)
190+ Nxy = S.Gs * gxy
191+ fx = Nxx * wx + Nxy * wy
192+ fy = Nxy * wx + Nyy * wy
193+ for k in range(4):
194+ fu = -(Nxx * dNdx[k] + Nxy * dNdy[k])
195+ fv = -(Nyy * dNdy[k] + Nxy * dNdx[k])
196+ fw = -(fx * dNdx[k] + fy * dNdy[k])
197+ Fu_c[k] = fu if Fu_c[k] is None else Fu_c[k] + fu
198+ Fv_c[k] = fv if Fv_c[k] is None else Fv_c[k] + fv
199+ Fw_c[k] = fw if Fw_c[k] is None else Fw_c[k] + fw
200+ # Corner k of every cell goes to node offsets (0,0), (0,1), (1,1), (1,0) in (j, i).
201+ pads = ((0, 1, 0, 1), (1, 0, 0, 1), (1, 0, 1, 0), (0, 1, 1, 0))
202+ Fu = sum(pad(Fu_c[k], *pads[k]) for k in range(4))
203+ Fv = sum(pad(Fv_c[k], *pads[k]) for k in range(4))
204+ Fw = Fw + sum(pad(Fw_c[k], *pads[k]) for k in range(4))
205+ return Fu * S.inv_mass, Fv * S.inv_mass, Fw * S.inv_mass
206+
207+
208+def damping(S: Static, du: torch.Tensor, dv: torch.Tensor, dw: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
209+ """The loss terms: -2 sigma0 v on every component, +2 sigma1 (graph Laplacian
210+ of w_t) on the transverse component, -2 sigma_in v on the in-plane components."""
211+ left, right, down, up = S.nb
212+ p = pad(dw, 1, 1, 1, 1)
213+ lap = (left * (p[1:-1, :-2] - dw) + right * (p[1:-1, 2:] - dw)) / S.hx**2 + (down * (p[:-2, 1:-1] - dw) + up * (p[2:, 1:-1] - dw)) / S.hy**2
214+ r = 2 * (S.sigma0 + S.sigma_in)
215+ return -r * du, -r * dv, -2 * S.sigma0 * dw + 2 * S.sigma1 * lap
216+
217+
218+def step(
219+ S: Static,
220+ u: torch.Tensor,
221+ v: torch.Tensor,
222+ w: torch.Tensor,
223+ wlo: torch.Tensor,
224+ du: torch.Tensor,
225+ dv: torch.Tensor,
226+ dw: torch.Tensor,
227+ stick: torch.Tensor,
228+) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
229+ """One step. ``stick`` holds (position, velocity, force, max force, time since release, withdrawn)."""
230+ dt = S.dt
231+ au, av, aw = accelerations(S, u, v, w)
232+ lu, lv, lw = damping(S, du, dv, dw)
233+ au = au + lu
234+ av = av + lv
235+ aw = aw + lw
236+ if S.strike:
237+ # Hertz contact between the free stick mass and the contact patch (the
238+ # plate nodes within the tip radius, plus the nearest one). The contact
239+ # coordinate is the mass-weighted mean displacement of the patch and the
240+ # force is spread over the patch in proportion to the node masses, so
241+ # the pair is conjugate. Once the contact has been released for
242+ # STICK_RELEASE the stick is withdrawn. The force of the previous step
243+ # is applied to the whole plate as the rigid acceleration field that
244+ # balances the patch load and its moments (inertia relief of the free
245+ # plate), so the plate gains no net momentum.
246+ sx, sv, f_prev, f_max, released, withdrawn = stick.unbind()
247+ aw = aw + f_prev * S.relief
248+ wc = (S.patch_w * w).sum()
249+ pen = wc - sx
250+ active = withdrawn == 0
251+ contact = active & (pen > 0)
252+ fc = torch.where(contact, S.K * torch.clamp(pen, min=0) ** S.alpha, torch.zeros_like(pen))
253+ released = torch.where(contact, torch.zeros_like(released), torch.where(active & (f_max > 0), released + dt, released))
254+ withdrawn = torch.where(active & (released > STICK_RELEASE), torch.ones_like(withdrawn), withdrawn)
255+ aw = aw - fc * S.patch_f
256+ svn = sv + dt * fc / S.stick_mass
257+ stick = torch.stack([sx + dt * svn, svn, fc, torch.maximum(f_max, fc), released, withdrawn])
258+ du = S.ins * (du + dt * au)
259+ dv = S.ins * (dv + dt * av)
260+ dw = S.ins * (dw + dt * aw)
261+ u = S.ins * (u + dt * du)
262+ v = S.ins * (v + dt * dv)
263+ # Compensated update of w: the increment dt v can be far below the float32
264+ # resolution of w, so the rounding error is carried separately.
265+ y = dt * dw + wlo
266+ t = w + y
267+ lo = (w - t) + y
268+ w = S.ins * t
269+ wlo = S.ins * lo
270+ return u, v, w, wlo, du, dv, dw, stick
src/cymsim/full/reference.pyadded+65−0View file
@@ -0,0 +1,65 @@
1+"""The discrete potential energy of the plate in NumPy (float64), whose negative
2+gradient is the force of the full model's kernel, and the largest frequency of
3+the linearized system, which sets the time step. Neither is used in a time loop;
4+they serve the tests and the set-up of a run."""
5+
6+from __future__ import annotations
7+
8+import numpy as np
9+import scipy.linalg
10+
11+from ..grid import Grid
12+from ..linear.assemble import assemble, corners, shape_derivatives
13+from ..material import Material
14+from ..plate import Plate
15+
16+
17+def potential_energy(pl: Plate, u: np.ndarray, v: np.ndarray, w: np.ndarray) -> float:
18+ """Bending plus membrane energy [J] of the displacement fields ``(u, v, w)``
19+ ``(ny, nx)``, measured from the stress-free shape."""
20+ g, m = pl.grid, pl.material
21+ hx, hy, nu = g.hx, g.hy, m.nu
22+ left, right, down, up = g.neighbor_masks()
23+ wp = np.pad(w, 1)
24+ a = np.where(left & right, (wp[1:-1, 2:] - 2 * w + wp[1:-1, :-2]) / hx**2, 0.0)
25+ b = np.where(down & up, (wp[2:, 1:-1] - 2 * w + wp[:-2, 1:-1]) / hy**2, 0.0)
26+ Ub = (g.node_weight * m.bending_stiffness(pl.h) * (a * a + b * b + 2 * nu * a * b)).sum()
27+ cc = (w[1:, 1:] - w[1:, :-1] - w[:-1, 1:] + w[:-1, :-1]) / (hx * hy)
28+ Ub += (2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * cc * cc * g.cell_mask).sum()
29+ Ub *= 0.5 * hx * hy
30+
31+ hc = pl.cell_thickness
32+ C, Gs = m.membrane_stiffness(hc), m.shear_stiffness(hc)
33+ Uc, Vc, Wc, W0c = corners(u), corners(v), corners(w), corners(pl.w0)
34+ Um = 0.0
35+ q = hx * hy / 4
36+ for gp in range(4):
37+ dNdx, dNdy = shape_derivatives(gp, hx, hy)
38+ ux = sum(dNdx[k] * Uc[k] for k in range(4))
39+ uy = sum(dNdy[k] * Uc[k] for k in range(4))
40+ vx = sum(dNdx[k] * Vc[k] for k in range(4))
41+ vy = sum(dNdy[k] * Vc[k] for k in range(4))
42+ dwx = sum(dNdx[k] * Wc[k] for k in range(4))
43+ dwy = sum(dNdy[k] * Wc[k] for k in range(4))
44+ w0x = sum(dNdx[k] * W0c[k] for k in range(4))
45+ w0y = sum(dNdy[k] * W0c[k] for k in range(4))
46+ exx = ux + dwx * (w0x + 0.5 * dwx)
47+ eyy = vy + dwy * (w0y + 0.5 * dwy)
48+ gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy
49+ Um += (q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy) * g.cell_mask).sum()
50+ return float(Ub + Um)
51+
52+
53+def max_angular_frequency(g: Grid, m: Material, h: float, n_small: int = 12) -> float:
54+ """Largest angular frequency of the linearized discrete system [rad/s], from
55+ the three-field stiffness of a small uniform flat plate with the same
56+ spacing and the thickness ``h`` (the largest of the plate, which bounds the
57+ bending frequency; the in-plane wave speed does not depend on the
58+ thickness). The extreme mode is a grid-scale mode, so neither the plate
59+ size nor its outline matters."""
60+ gs = Grid(n_small, n_small, g.hx, g.hy, np.ones((n_small, n_small), dtype=bool))
61+ sys = assemble(Plate.uniform(gs, m, h), fields=3)
62+ s = 1 / np.sqrt(sys.mass)
63+ A = sys.K.multiply(s[:, None]).multiply(s[None, :]).toarray()
64+ lam = scipy.linalg.eigvalsh(A, subset_by_index=(sys.n - 1, sys.n - 1))[0]
65+ return float(np.sqrt(max(lam, 0.0)))
src/cymsim/full/simulation.pyadded+281−0View file
@@ -0,0 +1,281 @@
1+"""The time-stepping driver of the full model. A run advances in batches of one
2+movie frame (1 / fps s): each batch steps the kernel, records the transverse
3+velocity at the audio sample instants into a ring buffer, and at its end forms
4+the batch's audio samples from the retarded-time Rayleigh sum over the plate
5+and reads back a frame of the transverse displacement. The recording starts at
6+the impact, with the plate at rest in its stress-free shape.
7+
8+The device is a CUDA GPU when one is available, else Apple's MPS, else the CPU.
9+The step is many small tensor operations, which ``torch.compile`` fuses into
10+a few kernels (the default; about 6 times faster than the plain step on the
11+CPU, where a 16 inch cymbal at 8 mm then runs at about 90 s per second of
12+sound, and far more on a GPU, where the plain step would be bound by kernel
13+launches).
14+"""
15+
16+from __future__ import annotations
17+
18+import math
19+import time
20+from typing import Literal
21+
22+import numpy as np
23+import torch
24+
25+from ..constants import CFL
26+from ..loss import fit_sigma, in_plane_rate
27+from ..params import Listener, Strike
28+from ..plate import Plate
29+from ..sound import Sound
30+from .kernel import Static, static_data, step
31+from .reference import max_angular_frequency
32+
33+
34+def default_device() -> torch.device:
35+ """CUDA when available, else MPS, else the CPU."""
36+ if torch.cuda.is_available():
37+ return torch.device("cuda")
38+ if getattr(torch.backends, "mps", None) is not None and torch.backends.mps.is_available():
39+ return torch.device("mps")
40+ return torch.device("cpu")
41+
42+
43+class Simulation:
44+ """A full-model run that can be advanced batch by batch.
45+
46+ Parameters mirror ``simulate``. After construction, ``run`` advances the
47+ simulation and ``sound`` returns what has been recorded so far; ``duration``
48+ may be raised and ``run`` called again to extend a finished run.
49+ """
50+
51+ def __init__(
52+ self,
53+ plate: Plate,
54+ strike: Strike = Strike(),
55+ listener: Listener = Listener(),
56+ *,
57+ duration: float = 4.0,
58+ fs: int = 48000,
59+ fps: float = 120.0,
60+ movie: bool = True,
61+ device=None,
62+ dtype: Literal["float32", "float64"] | torch.dtype = "float32",
63+ compile: bool | str = True,
64+ loss: tuple[float, float, float] | None = None,
65+ ):
66+ self.plate, self.strike, self.listener = plate, strike, listener
67+ self.fs, self.fps, self.duration, self.movie = fs, fps, duration, movie
68+ self.device = torch.device(device) if device is not None else default_device()
69+ self.dtype = getattr(torch, dtype) if isinstance(dtype, str) else dtype
70+ g, m = plate.grid, plate.material
71+
72+ # Stable time step, rounded so that an integer number of steps make one audio sample.
73+ self.omega_max = max_angular_frequency(g, m, plate.thickness_range[1])
74+ dt_max = CFL * 2 / self.omega_max
75+ self.steps_per_sample = math.ceil(1 / (fs * dt_max))
76+ self.dt = 1 / (fs * self.steps_per_sample)
77+ self.samples_per_batch = max(1, round(fs / fps))
78+ self.steps_per_batch = self.samples_per_batch * self.steps_per_sample
79+
80+ if loss is None:
81+ fit = fit_sigma(plate)
82+ loss = (fit.sigma0, fit.sigma1, in_plane_rate(plate))
83+ self.loss = loss
84+ self.static: Static = static_data(plate, strike, self.dt, *loss, self.device, self.dtype)
85+
86+ # Listener: weights and delays in samples relative to the nearest node; the ring buffer covers the spread.
87+ weight, delay = plate.listener_terms(listener)
88+ sel = weight > 0
89+ d = np.where(sel, (delay - delay[sel].min()) * fs, 0.0)
90+ self.hist_len = self.samples_per_batch + int(math.ceil(d.max())) + 4
91+ self._lw = torch.as_tensor(weight.reshape(-1), dtype=self.dtype, device=self.device)
92+ self._k0 = torch.as_tensor(np.floor(d).astype(np.int64).reshape(-1), device=self.device)
93+ self._fr = torch.as_tensor((d - np.floor(d)).reshape(-1), dtype=self.dtype, device=self.device)
94+ self._hist = torch.zeros((self.hist_len, g.n), dtype=self.dtype, device=self.device)
95+ self._node_idx = torch.arange(g.n, device=self.device)[None, :]
96+
97+ z = lambda: torch.zeros(g.shape, dtype=self.dtype, device=self.device) # noqa: E731
98+ self._state = (z(), z(), z(), z(), z(), z(), z())
99+ self._stick = torch.tensor([0.0, -strike.velocity, 0.0, 0.0, 0.0, 0.0], dtype=self.dtype, device=self.device)
100+ self.total_steps = 0
101+ self.samples = 0
102+ self.pressure: list[np.ndarray] = []
103+ self.frames: list[np.ndarray] = [np.zeros(g.shape, dtype=np.float32)]
104+ self.stick = {"x": 0.0, "v": -strike.velocity, "force": 0.0, "max_force": 0.0, "withdrawn": False}
105+ self.seconds = 0.0
106+
107+ self._step = step
108+ self._compile_pending = bool(compile)
109+ self._compile_mode = compile if isinstance(compile, str) else ("reduce-overhead" if self.device.type == "cuda" else "default")
110+
111+ # ---- progress ----
112+
113+ @property
114+ def target_samples(self) -> int:
115+ return round(self.duration * self.fs)
116+
117+ @property
118+ def done(self) -> bool:
119+ return self.samples >= self.target_samples
120+
121+ @property
122+ def t(self) -> float:
123+ """Simulated time since the impact [s]."""
124+ return self.total_steps * self.dt
125+
126+ def status(self) -> str:
127+ return f"t = {self.t:.3f} s, {self.total_steps} steps, {self.samples} samples, stick force {self.stick['force']:.1f} N (max {self.stick['max_force']:.1f})" + (", withdrawn" if self.stick["withdrawn"] else "")
128+
129+ # ---- stepping ----
130+
131+ def run_batch(self) -> None:
132+ """Advances by one movie frame (or to the end of the recording)."""
133+ steps_left = self.target_samples * self.steps_per_sample - self.total_steps
134+ steps = min(self.steps_per_batch, steps_left)
135+ if steps <= 0:
136+ return
137+ t0 = time.perf_counter()
138+ sps, H = self.steps_per_sample, self.hist_len
139+ slots = []
140+ state, stick = self._state, self._stick
141+ S = self.static
142+ if self._compile_pending:
143+ self._try_compile()
144+ for s in range(steps):
145+ local = self.total_steps + s
146+ state_and_stick = self._step(S, *state, stick)
147+ state, stick = state_and_stick[:7], state_and_stick[7]
148+ if local % sps == 0:
149+ slot = (local // sps) % H
150+ self._hist[slot] = state[6].reshape(-1)
151+ slots.append(slot)
152+ self._state, self._stick = state, stick
153+ self.total_steps += steps
154+ if slots:
155+ self.pressure.append(self._audio(slots))
156+ self.samples += len(slots)
157+ if self.movie:
158+ self.frames.append(state[2].detach().to("cpu", torch.float32).numpy().copy())
159+ st = stick.detach().cpu().numpy()
160+ self.stick = {"x": float(st[0]), "v": float(st[1]), "force": float(st[2]), "max_force": float(st[3]), "withdrawn": bool(st[5] > 0)}
161+ self.seconds += time.perf_counter() - t0
162+
163+ def _try_compile(self) -> None:
164+ """Wraps the step in ``torch.compile`` (which fuses its many small tensor
165+ operations; about 6 times faster on the CPU and needed on a GPU) and
166+ checks it on the current state; falls back to the eager step with a
167+ warning when no compiler is available."""
168+ import warnings
169+
170+ self._compile_pending = False
171+ try:
172+ compiled = torch.compile(step, dynamic=False, mode=self._compile_mode)
173+ compiled(self.static, *self._state, self._stick)
174+ self._step = compiled
175+ except Exception as e: # noqa: BLE001
176+ warnings.warn(f"torch.compile failed ({type(e).__name__}: {str(e).splitlines()[0][:200]}); running the uncompiled step, which is slower", RuntimeWarning, stacklevel=3)
177+ self._step = step
178+
179+ def _audio(self, slots: list[int]) -> np.ndarray:
180+ """The retarded-time Rayleigh sum for the batch's samples: the acceleration
181+ is the first difference of the delayed velocity, the average over one
182+ sample interval, which is band-limited and immune to the per-step
183+ rounding noise of the force evaluation."""
184+ H = self.hist_len
185+ slot = torch.as_tensor(slots, device=self.device)[:, None]
186+ i0 = (slot - self._k0[None, :]) % H
187+ i1 = (i0 - 1) % H
188+ i2 = (i0 - 2) % H
189+ v0 = torch.gather(self._hist, 0, i0)
190+ v1 = torch.gather(self._hist, 0, i1)
191+ v2 = torch.gather(self._hist, 0, i2)
192+ fr = self._fr[None, :]
193+ a = ((v0 + fr * (v1 - v0)) - (v1 + fr * (v2 - v1))) * self.fs
194+ return (a @ self._lw).to("cpu", torch.float64).numpy()
195+
196+ def run(self, seconds: float | None = None, progress: bool = False) -> "Simulation":
197+ """Advances by ``seconds`` of simulated time (to the end by default). With ``progress`` a line is printed per second of sound."""
198+ target = self.target_samples if seconds is None else min(self.target_samples, self.samples + round(seconds * self.fs))
199+ next_report = self.samples + self.fs
200+ while self.samples < target and not self.done:
201+ self.run_batch()
202+ if progress and self.samples >= next_report:
203+ print(f"{self.t:.2f} s simulated in {self.seconds:.1f} s wall; {self.status()}", flush=True)
204+ next_report += self.fs
205+ return self
206+
207+ def sound(self) -> Sound:
208+ """What has been recorded so far."""
209+ pressure = np.concatenate(self.pressure) if self.pressure else np.zeros(0)
210+ frames = np.stack(self.frames) if self.movie else None
211+ info = {
212+ "model": "full",
213+ "dt": self.dt,
214+ "steps_per_sample": self.steps_per_sample,
215+ "steps": self.total_steps,
216+ "sigma0": self.loss[0],
217+ "sigma1": self.loss[1],
218+ "sigma_in": self.loss[2],
219+ "peak_force": self.stick["max_force"],
220+ "device": str(self.device),
221+ "dtype": str(self.dtype).replace("torch.", ""),
222+ "seconds": self.seconds,
223+ }
224+ return Sound(pressure.astype(np.float32), self.fs, self.plate, self.strike, self.listener, frames, self.fps, info)
225+
226+ # ---- state access, for tests and experiments ----
227+
228+ def state(self) -> dict[str, np.ndarray]:
229+ """The displacement and velocity fields ``u, v, w, du, dv, dw`` on the grid, as float64 arrays."""
230+ u, v, w, _, du, dv, dw = (x.detach().to("cpu", torch.float64).numpy() for x in self._state)
231+ return {"u": u, "v": v, "w": w, "du": du, "dv": dv, "dw": dw}
232+
233+ def set_state(self, u=None, v=None, w=None, du=None, dv=None, dw=None) -> None:
234+ """Overwrites fields of the state (grid arrays); the compensation array of w is reset."""
235+ cur = list(self._state)
236+ for k, x in ((0, u), (1, v), (2, w), (4, du), (5, dv), (6, dw)):
237+ if x is not None:
238+ cur[k] = torch.as_tensor(np.asarray(x, dtype=np.float64), dtype=self.dtype, device=self.device) * self.static.ins
239+ cur[3] = torch.zeros_like(cur[2])
240+ self._state = tuple(cur)
241+
242+ def step_once(self, strike: bool = False) -> None:
243+ """Advances the state by one step without recording; without the drumstick unless ``strike``."""
244+ S = self.static
245+ if not strike and S.strike:
246+ from dataclasses import replace
247+
248+ S = replace(S, strike=False)
249+ out = step(S, *self._state, self._stick)
250+ self._state, self._stick = out[:7], out[7]
251+ self.total_steps += 1
252+
253+
254+def simulate(
255+ plate: Plate,
256+ strike: Strike = Strike(),
257+ listener: Listener = Listener(),
258+ *,
259+ duration: float = 4.0,
260+ fs: int = 48000,
261+ fps: float = 120.0,
262+ movie: bool = False,
263+ device=None,
264+ dtype: Literal["float32", "float64"] | torch.dtype = "float32",
265+ compile: bool | str = True,
266+ progress: bool = False,
267+) -> Sound:
268+ """The sound of a strike by the full (geometrically nonlinear) model.
269+
270+ ``device`` is a torch device (``"cuda"``, ``"cuda:1"``, ``"mps"``, ``"cpu"``);
271+ by default a CUDA GPU when one is available. ``dtype`` is ``"float32"`` (the
272+ default, with the rounding safeguards of the kernel) or ``"float64"``.
273+ ``compile`` wraps the step in ``torch.compile`` (``True``, the default,
274+ with a fallback to the plain step when no compiler is available; ``False``;
275+ or a ``torch.compile`` mode string). ``fps`` is the movie frame rate and
276+ also the batch rate of the run. ``loss`` overrides the fitted
277+ ``(sigma0, sigma1, sigma_in)`` of the two-parameter loss operator.
278+ """
279+ sim = Simulation(plate, strike, listener, duration=duration, fs=fs, fps=fps, movie=movie, device=device, dtype=dtype, compile=compile)
280+ sim.run(progress=progress)
281+ return sim.sound()
src/cymsim/grid.pyadded+141−0View file
@@ -0,0 +1,141 @@
1+"""The regular grid that covers the bounding box of a plate, restricted to the
2+plate's outline by a mask. Arrays over the grid have shape ``(ny, nx)``, indexed
3+``[j, i]`` with ``i`` along x. Node coordinates are relative to the center of the
4+bounding box.
5+
6+The edge is a staircase of whole cells: a cell belongs to the plate when its
7+center lies inside the outline, and a node when it is a corner of such a cell.
8+Deciding membership by cell centers keeps the covered area unbiased. A hole is
9+handled by the same rule.
10+"""
11+
12+from __future__ import annotations
13+
14+from dataclasses import dataclass
15+from functools import cached_property
16+from typing import Callable
17+
18+import numpy as np
19+
20+InsideFn = Callable[[np.ndarray, np.ndarray, np.ndarray, np.ndarray], np.ndarray]
21+"""``inside(x, y, r, theta) -> bool array``, evaluated on whole grid matrices."""
22+
23+
24+@dataclass(frozen=True)
25+class Grid:
26+ nx: int
27+ ny: int
28+ hx: float
29+ hy: float
30+ mask: np.ndarray
31+ """``(ny, nx)`` bool: True for the nodes that belong to the plate."""
32+
33+ @property
34+ def n(self) -> int:
35+ """Number of grid nodes, inside or not."""
36+ return self.nx * self.ny
37+
38+ @property
39+ def shape(self) -> tuple[int, int]:
40+ return (self.ny, self.nx)
41+
42+ @property
43+ def count(self) -> int:
44+ """Number of nodes in the plate."""
45+ return int(self.mask.sum())
46+
47+ @cached_property
48+ def x(self) -> np.ndarray:
49+ """Node x coordinates ``(nx,)``, relative to the center of the bounding box."""
50+ return (np.arange(self.nx) - (self.nx - 1) / 2) * self.hx
51+
52+ @cached_property
53+ def y(self) -> np.ndarray:
54+ return (np.arange(self.ny) - (self.ny - 1) / 2) * self.hy
55+
56+ @cached_property
57+ def X(self) -> np.ndarray:
58+ """Node x coordinates on the grid, ``(ny, nx)``."""
59+ return np.broadcast_to(self.x[None, :], self.shape).copy()
60+
61+ @cached_property
62+ def Y(self) -> np.ndarray:
63+ return np.broadcast_to(self.y[:, None], self.shape).copy()
64+
65+ @cached_property
66+ def cell_mask(self) -> np.ndarray:
67+ """``(ny - 1, nx - 1)`` bool: cells whose four corners belong to the plate."""
68+ m = self.mask
69+ return m[:-1, :-1] & m[:-1, 1:] & m[1:, 1:] & m[1:, :-1]
70+
71+ @cached_property
72+ def node_weight(self) -> np.ndarray:
73+ """Quadrature and lumped-mass weight of each node: the fraction of its
74+ four cells that belong to the plate (1 inside, 1/2 on an edge, 1/4 at a
75+ corner of a rectangle, 0 outside)."""
76+ return corners_to_nodes(self.cell_mask.astype(float)) / 4
77+
78+ @cached_property
79+ def node_area(self) -> np.ndarray:
80+ """Area attributed to each node, hx hy times its weight [m^2]."""
81+ return self.hx * self.hy * self.node_weight
82+
83+ @cached_property
84+ def area(self) -> float:
85+ """Area of the staircase plate [m^2]."""
86+ return float(self.node_area.sum())
87+
88+ def cell_centers(self) -> tuple[np.ndarray, np.ndarray]:
89+ """Coordinates of the cell centers, two ``(ny - 1, nx - 1)`` arrays."""
90+ return np.meshgrid(self.x[:-1] + self.hx / 2, self.y[:-1] + self.hy / 2)
91+
92+ def neighbor_masks(self) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
93+ """Whether the left, right, lower (j - 1) and upper (j + 1) neighbor of each node is in the plate."""
94+ m = self.mask
95+ left = np.zeros_like(m)
96+ left[:, 1:] = m[:, :-1]
97+ right = np.zeros_like(m)
98+ right[:, :-1] = m[:, 1:]
99+ down = np.zeros_like(m)
100+ down[1:, :] = m[:-1, :]
101+ up = np.zeros_like(m)
102+ up[:-1, :] = m[1:, :]
103+ return left & m, right & m, down & m, up & m
104+
105+
106+def corners_to_nodes(c: np.ndarray) -> np.ndarray:
107+ """Sums a per-cell array ``(ny - 1, nx - 1)`` onto the four corners of each cell, giving ``(ny, nx)``."""
108+ out = np.zeros((c.shape[0] + 1, c.shape[1] + 1), dtype=c.dtype)
109+ out[:-1, :-1] += c
110+ out[:-1, 1:] += c
111+ out[1:, 1:] += c
112+ out[1:, :-1] += c
113+ return out
114+
115+
116+def make_grid(Lx: float, Ly: float, spacing: float, inside: InsideFn | None = None) -> Grid:
117+ """The grid covering the box ``Lx`` by ``Ly`` with about the requested
118+ spacing (at least 6 nodes per side), restricted to ``inside`` when given."""
119+ nx = max(6, round(Lx / spacing) + 1)
120+ ny = max(6, round(Ly / spacing) + 1)
121+ g = Grid(nx, ny, Lx / (nx - 1), Ly / (ny - 1), np.ones((ny, nx), dtype=bool))
122+ if inside is None:
123+ return g
124+ XC, YC = g.cell_centers()
125+ cells = np.asarray(inside(XC, YC, np.hypot(XC, YC), np.arctan2(YC, XC)), dtype=bool)
126+ if cells.shape != (ny - 1, nx - 1):
127+ raise ValueError(f"inside returned shape {cells.shape} for a {ny - 1} x {nx - 1} array of cells")
128+ return with_cell_mask(g, cells)
129+
130+
131+def with_cell_mask(g: Grid, cells: np.ndarray) -> Grid:
132+ """Restricts a grid to the cells flagged in ``cells`` ``(ny - 1, nx - 1)``:
133+ a node belongs to the plate when it is a corner of an included cell. Any
134+ node left with no cell of its own would carry no mass and is removed."""
135+ mask = corners_to_nodes(cells.astype(np.int8)) > 0
136+ while True:
137+ gm = Grid(g.nx, g.ny, g.hx, g.hy, mask)
138+ bad = mask & (gm.node_weight == 0)
139+ if not bad.any():
140+ return gm
141+ mask = mask & ~bad
src/cymsim/linear/__init__.pyadded+8−0View file
@@ -0,0 +1,8 @@
1+"""The linear model: the modes of the plate linearized about its rest shape, their
2+loss rates, and the modal synthesis of a strike."""
3+
4+from .assemble import LinearSystem, assemble
5+from .modes import Modes, compute_modes
6+from .synth import synthesize
7+
8+__all__ = ["LinearSystem", "assemble", "Modes", "compute_modes", "synthesize"]
src/cymsim/linear/assemble.pyadded+186−0View file
@@ -0,0 +1,186 @@
1+"""The linearized plate as sparse matrices: the stiffness K (the Hessian of the
2+discrete energy at the stress-free state) and the lumped mass M as a diagonal.
3+
4+A flat plate has one unknown per node (the transverse displacement w; the
5+in-plane motion decouples and is neither struck nor heard). A plate with a
6+profile couples w to the in-plane displacements u and v through the slopes of
7+the profile, and has three unknowns per node, interleaved (u, v, w) so that
8+the band stays narrow. The discretization is that of the full model
9+(``full.kernel``): three-point second differences for the bending curvatures at
10+the nodes, the four-point mixed difference at the cells, and bilinear cells
11+with 2 x 2 Gauss integration for the membrane strains.
12+"""
13+
14+from __future__ import annotations
15+
16+from dataclasses import dataclass
17+
18+import numpy as np
19+import scipy.sparse as sp
20+
21+from ..plate import Plate
22+
23+G3 = 1 / np.sqrt(3.0)
24+
25+
26+def shape_derivatives(gp: int, hx: float, hy: float) -> tuple[np.ndarray, np.ndarray]:
27+ """Derivatives of the four bilinear shape functions at Gauss point ``gp`` (0..3).
28+ The corners are ordered (ci, cj), (ci + 1, cj), (ci + 1, cj + 1), (ci, cj + 1)."""
29+ xi = G3 if gp & 1 else -G3
30+ eta = G3 if gp & 2 else -G3
31+ dNdx = np.array([-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)]) / (2 * hx)
32+ dNdy = np.array([-(1 - xi), -(1 + xi), 1 + xi, 1 - xi]) / (2 * hy)
33+ return dNdx, dNdy
34+
35+
36+def corners(a: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
37+ """The values of a node array at the four corners of every cell, in the corner order above."""
38+ return a[:-1, :-1], a[:-1, 1:], a[1:, 1:], a[1:, :-1]
39+
40+
41+@dataclass
42+class LinearSystem:
43+ K: sp.csr_matrix
44+ """Symmetric stiffness, ``(n, n)``."""
45+ mass: np.ndarray
46+ """Diagonal mass of every unknown [kg], ``(n,)``."""
47+ fields: int
48+ """1 for a flat plate (w only), 3 for a plate with a profile (u, v, w interleaved)."""
49+ dof: np.ndarray
50+ """``(ny, nx)``: first unknown of each grid node, or -1 outside the plate."""
51+ nodes: np.ndarray
52+ """``(n_nodes,)`` flat grid index ``j * nx + i`` of each node, in unknown order."""
53+
54+ @property
55+ def n(self) -> int:
56+ return len(self.mass)
57+
58+ @property
59+ def bandwidth(self) -> int:
60+ """Half-bandwidth of K."""
61+ coo = self.K.tocoo()
62+ return int(np.abs(coo.row - coo.col).max()) if coo.nnz else 0
63+
64+ def w_dof(self) -> np.ndarray:
65+ """``(ny, nx)``: the transverse unknown of each node, or -1 outside."""
66+ return np.where(self.dof >= 0, self.dof + (self.fields - 1), -1)
67+
68+ def to_grid(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
69+ """Scatters an unknown vector (or a stack ``(..., n)``) to grid arrays ``(u, v, w)``, zero outside the plate."""
70+ x = np.asarray(x)
71+ shape = x.shape[:-1] + (self.dof.shape[0] * self.dof.shape[1],)
72+ u, v, w = np.zeros(shape), np.zeros(shape), np.zeros(shape)
73+ if self.fields == 1:
74+ w[..., self.nodes] = x
75+ else:
76+ u[..., self.nodes] = x[..., 0::3]
77+ v[..., self.nodes] = x[..., 1::3]
78+ w[..., self.nodes] = x[..., 2::3]
79+ gs = x.shape[:-1] + self.dof.shape
80+ return u.reshape(gs), v.reshape(gs), w.reshape(gs)
81+
82+ def from_grid(self, u: np.ndarray, v: np.ndarray, w: np.ndarray) -> np.ndarray:
83+ """Gathers grid arrays into an unknown vector."""
84+ if self.fields == 1:
85+ return w.reshape(-1)[self.nodes].copy()
86+ x = np.empty(self.n)
87+ x[0::3] = u.reshape(-1)[self.nodes]
88+ x[1::3] = v.reshape(-1)[self.nodes]
89+ x[2::3] = w.reshape(-1)[self.nodes]
90+ return x
91+
92+
93+def assemble(pl: Plate, fields: int | None = None) -> LinearSystem:
94+ """Assembles K and M for the plate. ``fields`` may force the three-field
95+ system for a flat plate (the tests use it); by default a flat plate gets the
96+ one-field system."""
97+ g, m = pl.grid, pl.material
98+ nx, ny, hx, hy = g.nx, g.ny, g.hx, g.hy
99+ nf = fields if fields is not None else (1 if pl.is_flat else 3)
100+ if nf not in (1, 3):
101+ raise ValueError("fields must be 1 or 3")
102+ nu = m.nu
103+
104+ # Number the plate nodes with the shorter grid axis innermost, so that the
105+ # 13-point bending stencil, which reaches two rows away, gives the narrowest band.
106+ inside = g.mask & (g.node_weight > 0)
107+ if nx <= ny:
108+ nodes = np.flatnonzero(inside.ravel())
109+ else:
110+ jj, ii = np.nonzero(inside.T) # (i, j) pairs in i-major order
111+ nodes = ii * nx + jj
112+ n_nodes = len(nodes)
113+ n = n_nodes * nf
114+ dof = -np.ones(g.n, dtype=np.int64)
115+ dof[nodes] = np.arange(n_nodes) * nf
116+ dof = dof.reshape(g.shape)
117+ mass = np.repeat(pl.node_mass.reshape(-1)[nodes], nf)
118+ wd = dof + (nf - 1) # w unknown of each node (garbage outside; masked below)
119+
120+ def stencil(rows_mask: np.ndarray, cols: list[np.ndarray], vals: list[float], n_rows: int) -> sp.csr_matrix:
121+ """A sparse operator with one row per grid entity, nonzero where ``rows_mask``."""
122+ r = np.flatnonzero(rows_mask.ravel())
123+ data, ri, ci = [], [], []
124+ for c, v in zip(cols, vals):
125+ data.append(np.full(len(r), v))
126+ ri.append(r)
127+ ci.append(c.ravel()[r])
128+ return sp.csr_matrix((np.concatenate(data), (np.concatenate(ri), np.concatenate(ci))), shape=(n_rows, n))
129+
130+ def shifted(a: np.ndarray, di: int, dj: int) -> np.ndarray:
131+ """``a[j + dj, i + di]`` with -1 (an invalid column) outside the grid."""
132+ out = -np.ones_like(a)
133+ src = a[max(0, dj) : ny + min(0, dj), max(0, di) : nx + min(0, di)]
134+ out[max(0, -dj) : ny + min(0, -dj), max(0, -di) : nx + min(0, -di)] = src
135+ return out
136+
137+ ins = g.mask
138+ left, right, down, up = g.neighbor_masks()
139+ has_a = ins & left & right
140+ has_b = ins & down & up
141+
142+ # ---- Bending: hx hy [ sum_k w_k D_k (a^2 + b^2 + 2 nu a b) + 2 (1 - nu) sum_c D_c c^2 ] ----
143+ A = stencil(has_a, [shifted(wd, -1, 0), wd, shifted(wd, 1, 0)], [1 / hx**2, -2 / hx**2, 1 / hx**2], g.n)
144+ B = stencil(has_b, [shifted(wd, 0, -1), wd, shifted(wd, 0, 1)], [1 / hy**2, -2 / hy**2, 1 / hy**2], g.n)
145+ Fa = sp.diags((hx * hy * g.node_weight * m.bending_stiffness(pl.h)).ravel())
146+ K = A.T @ Fa @ A + B.T @ Fa @ B + nu * (A.T @ Fa @ B + B.T @ Fa @ A)
147+ cm = g.cell_mask
148+ c00, c10, c11, c01 = corners(wd)
149+ C = stencil(cm, [c00, c10, c01, c11], [1 / (hx * hy), -1 / (hx * hy), -1 / (hx * hy), 1 / (hx * hy)], cm.size)
150+ Fc = sp.diags((2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * hx * hy * cm).ravel())
151+ K = K + C.T @ Fc @ C
152+
153+ # ---- Membrane, linearized about the profile (three-field system only) ----
154+ if nf == 3:
155+ q = hx * hy / 4
156+ hc = pl.cell_thickness
157+ Cm = sp.diags((q * m.membrane_stiffness(hc) * cm).ravel())
158+ Gs = sp.diags((q * m.shear_stiffness(hc) * cm).ravel())
159+ u_c = [c for c in corners(dof)]
160+ v_c = [c + 1 for c in corners(dof)]
161+ w_c = [c + 2 for c in corners(dof)]
162+ w0_c = corners(pl.w0)
163+ for gp in range(4):
164+ dNdx, dNdy = shape_derivatives(gp, hx, hy)
165+ w0x = sum(dNdx[b] * w0_c[b] for b in range(4))
166+ w0y = sum(dNdy[b] * w0_c[b] for b in range(4))
167+ # Linear strains: exx = u_x + w0_x w_x, eyy = v_y + w0_y w_y, gxy = u_y + v_x + w0_x w_y + w0_y w_x.
168+ Bxx = var_stencil(cm, u_c + w_c, [np.full(cm.shape, dNdx[b]) for b in range(4)] + [w0x * dNdx[b] for b in range(4)], n)
169+ Byy = var_stencil(cm, v_c + w_c, [np.full(cm.shape, dNdy[b]) for b in range(4)] + [w0y * dNdy[b] for b in range(4)], n)
170+ Bxy = var_stencil(cm, u_c + v_c + w_c, [np.full(cm.shape, dNdy[b]) for b in range(4)] + [np.full(cm.shape, dNdx[b]) for b in range(4)] + [w0x * dNdy[b] + w0y * dNdx[b] for b in range(4)], n)
171+ K = K + Bxx.T @ Cm @ Bxx + Byy.T @ Cm @ Byy + nu * (Bxx.T @ Cm @ Byy + Byy.T @ Cm @ Bxx) + Bxy.T @ Gs @ Bxy
172+
173+ K = sp.csr_matrix(K)
174+ K = (K + K.T) * 0.5
175+ K.eliminate_zeros()
176+ K.sort_indices()
177+ return LinearSystem(K, mass, nf, dof, nodes)
178+
179+
180+def var_stencil(rows_mask: np.ndarray, cols: list[np.ndarray], vals: list[np.ndarray], n: int) -> sp.csr_matrix:
181+ """A sparse operator with one row per cell and per-cell coefficients."""
182+ r = np.flatnonzero(rows_mask.ravel())
183+ data = np.concatenate([v.ravel()[r] for v in vals])
184+ ri = np.concatenate([r] * len(cols))
185+ ci = np.concatenate([c.ravel()[r] for c in cols])
186+ return sp.csr_matrix((data, (ri, ci)), shape=(rows_mask.size, n))
src/cymsim/linear/damping.pyadded+121−0View file
@@ -0,0 +1,121 @@
1+"""Per-mode decay rates of the linear model from the mode shapes: the radiation
2+into the air and the thermoelastic loss. Both are exact for the discrete model
3+given the shape, so that the decay of each mode is a result of the geometry and
4+the material rather than a prescribed curve.
5+
6+Radiation. A baffled plate whose nodes move with the velocities v_k radiates
7+the power P = v^T R v with the radiation resistance matrix
8+R_jk = (rho_air omega^2 / (4 pi c)) dA_j dA_k sinc(k r_jk) (Elliott and Johnson
9+1993), one face; both faces double it. For a mass-normalized mode of amplitude
10+A the velocity is omega A phi, the energy is omega^2 A^2 / 2, and the amplitude
11+decays at sigma = P / (2 E), which gives
12+ sigma_rad = (rho_air omega^2 / (2 pi c)) sum_jk a_j a_k sinc(k r_jk),
13+with a_k = dA_k phi(k). On a regular grid the kernel depends only on the offset
14+j - k, so the double sum is a convolution, and by Parseval it is a sum over the
15+spectrum: two 2-D FFTs on a zero-padded grid per mode instead of N^2 terms.
16+
17+Thermoelastic. Zener's loss factor applies to the bending energy and depends on
18+the local thickness, so a mode's rate is the sum over the nodes of the local
19+loss factor times the local bending energy, divided by omega: for a
20+mass-normalized mode the total strain energy is omega^2 / 2, and a flat uniform
21+plate recovers eta_te omega / 2.
22+"""
23+
24+from __future__ import annotations
25+
26+import math
27+
28+import numpy as np
29+import scipy.fft
30+
31+from ..constants import RHO_AIR, SPEED_OF_SOUND
32+from ..loss import thermoelastic_loss_factor
33+from ..plate import Plate
34+from .modes import Modes
35+
36+
37+def radiation_rates(modes: Modes, batch: int = 64) -> np.ndarray:
38+ """The radiation decay rates [1/s] of all modes by the FFT convolution."""
39+ pl = modes.plate
40+ g = pl.grid
41+ nx, ny, hx, hy = g.nx, g.ny, g.hx, g.hy
42+ Px, Py = scipy.fft.next_fast_len(2 * nx - 1, real=True), scipy.fft.next_fast_len(2 * ny - 1, real=True)
43+ c = SPEED_OF_SOUND
44+ dx = np.arange(Px)
45+ dx = np.where(dx <= Px // 2, dx, dx - Px)
46+ dy = np.arange(Py)
47+ dy = np.where(dy <= Py // 2, dy, dy - Py)
48+ r = np.hypot(dx[None, :] * hx, dy[:, None] * hy)
49+ reach = (np.abs(dx)[None, :] <= nx - 1) & (np.abs(dy)[:, None] <= ny - 1)
50+ area = g.node_area
51+ out = np.zeros(modes.count)
52+ for start in range(0, modes.count, batch):
53+ stop = min(modes.count, start + batch)
54+ omega = 2 * math.pi * modes.freqs[start:stop]
55+ k = omega / c
56+ a = np.zeros((stop - start, Py, Px))
57+ a[:, :ny, :nx] = area[None] * modes.w[start:stop]
58+ A = scipy.fft.rfft2(a, axes=(-2, -1))
59+ x = k[:, None, None] * r[None]
60+ kern = np.where(x < 1e-9, 1.0, np.sin(x) / np.where(x < 1e-9, 1.0, x)) * reach[None]
61+ S = scipy.fft.rfft2(kern, axes=(-2, -1)).real
62+ # Parseval over the half spectrum: interior columns count twice.
63+ wcol = np.full(Px // 2 + 1, 2.0)
64+ wcol[0] = 1.0
65+ if Px % 2 == 0:
66+ wcol[-1] = 1.0
67+ total = (S * (A.real**2 + A.imag**2) * wcol[None, None, :]).sum(axis=(-2, -1)) / (Px * Py)
68+ out[start:stop] = RHO_AIR * omega**2 / (2 * math.pi * c) * total
69+ return out
70+
71+
72+def radiation_rate_direct(pl: Plate, shape: np.ndarray, freq: float) -> float:
73+ """The same rate by the direct double sum, for the tests and small plates."""
74+ g = pl.grid
75+ sel = g.node_area > 0
76+ a = (g.node_area * shape)[sel]
77+ xs, ys = g.X[sel], g.Y[sel]
78+ omega = 2 * math.pi * freq
79+ k = omega / SPEED_OF_SOUND
80+ x = k * np.hypot(xs[:, None] - xs[None, :], ys[:, None] - ys[None, :])
81+ sinc = np.where(x < 1e-9, 1.0, np.sin(x) / np.where(x < 1e-9, 1.0, x))
82+ return float(RHO_AIR * omega**2 / (2 * math.pi * SPEED_OF_SOUND) * (a[:, None] * a[None, :] * sinc).sum())
83+
84+
85+def bending_energy_per_node(pl: Plate, w: np.ndarray) -> np.ndarray:
86+ """The bending energy of transverse shapes ``w`` ``(..., ny, nx)`` on the grid,
87+ split over the nodes (the mixed term of each cell is shared by its four
88+ corners), as in the discrete energy of the full model."""
89+ g, m = pl.grid, pl.material
90+ hx, hy, nu = g.hx, g.hy, m.nu
91+ left, right, down, up = g.neighbor_masks()
92+ has_a = left & right
93+ has_b = down & up
94+ wp = np.pad(w, [(0, 0)] * (w.ndim - 2) + [(1, 1), (1, 1)])
95+ a = np.where(has_a, (wp[..., 1:-1, 2:] - 2 * w + wp[..., 1:-1, :-2]) / hx**2, 0.0)
96+ b = np.where(has_b, (wp[..., 2:, 1:-1] - 2 * w + wp[..., :-2, 1:-1]) / hy**2, 0.0)
97+ U = 0.5 * hx * hy * g.node_weight * m.bending_stiffness(pl.h) * (a * a + b * b + 2 * nu * a * b)
98+ cc = (w[..., 1:, 1:] - w[..., 1:, :-1] - w[..., :-1, 1:] + w[..., :-1, :-1]) / (hx * hy)
99+ Uc = 0.5 * hx * hy * 2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * cc * cc * g.cell_mask / 4
100+ U[..., :-1, :-1] += Uc
101+ U[..., :-1, 1:] += Uc
102+ U[..., 1:, 1:] += Uc
103+ U[..., 1:, :-1] += Uc
104+ return U
105+
106+
107+def thermoelastic_rates(modes: Modes, batch: int = 64) -> np.ndarray:
108+ """The thermoelastic decay rates [1/s]: sum_k eta_te(omega, h_k) U_k / omega with U_k the bending energy at node k."""
109+ pl = modes.plate
110+ out = np.zeros(modes.count)
111+ if pl.material.thermal is None:
112+ return out
113+ sel = pl.grid.node_area > 0
114+ h = pl.h[sel]
115+ for start in range(0, modes.count, batch):
116+ stop = min(modes.count, start + batch)
117+ omega = 2 * math.pi * modes.freqs[start:stop]
118+ U = bending_energy_per_node(pl, modes.w[start:stop])[:, sel]
119+ eta = thermoelastic_loss_factor(pl.material, h[None, :], omega[:, None], pl.temperature)
120+ out[start:stop] = (eta * U).sum(axis=1) / omega
121+ return out
src/cymsim/linear/modes.pyadded+271−0View file
@@ -0,0 +1,271 @@
1+"""The modes of a plate up to a cutoff frequency, and their loss rates.
2+
3+The generalized problem K phi = omega^2 M phi with the diagonal lumped mass M
4+is turned into the standard symmetric problem A = M^{-1/2} K M^{-1/2}, whose
5+eigenvalues in the band of interest are computed by LAPACK, on the CPU through
6+SciPy (a dense solve restricted to the band; about 4 s for a 16 inch cymbal at
7+8 mm, and 40 s for the same cymbal with a profile, whose system is three times
8+as large) or on a GPU through PyTorch (a full dense eigendecomposition in
9+double precision, the faster path for the larger systems when a CUDA device
10+is available). The mode shapes
11+returned are M^{-1/2} times the eigenvectors of A, that is, mass-normalized.
12+The zero-frequency rigid modes of the free plate are dropped.
13+"""
14+
15+from __future__ import annotations
16+
17+import math
18+import time
19+from dataclasses import dataclass, replace
20+from typing import Literal
21+
22+import numpy as np
23+import scipy.linalg
24+
25+from ..constants import LN1000
26+from ..plate import Plate
27+from .assemble import LinearSystem, assemble
28+
29+Solver = Literal["auto", "banded", "dense", "torch"]
30+
31+
32+@dataclass
33+class Modes:
34+ """The modes of a plate below a cutoff frequency, with their per-mode decay rates.
35+
36+ Shapes are mass-normalized (phi^T M phi = 1) and given on the grid, zero
37+ outside the plate. For a flat plate only the transverse component ``w``
38+ exists; a plate with a profile also has the in-plane components ``u`` and ``v``.
39+ """
40+
41+ plate: Plate
42+ freqs: np.ndarray
43+ """Frequencies [Hz], ascending, ``(M,)``."""
44+ w: np.ndarray
45+ """Transverse component of the shapes, ``(M, ny, nx)``."""
46+ u: np.ndarray | None
47+ v: np.ndarray | None
48+ sigma_material: np.ndarray
49+ """Decay rate of each mode from internal friction [1/s]."""
50+ sigma_thermal: np.ndarray
51+ """Decay rate from thermoelastic damping [1/s]."""
52+ sigma_radiation: np.ndarray
53+ """Decay rate from radiation into the air [1/s]."""
54+ fmax: float
55+ """Cutoff frequency of the computation [Hz]."""
56+ seconds: float = 0.0
57+ """Wall time of the computation [s]."""
58+
59+ @property
60+ def count(self) -> int:
61+ return len(self.freqs)
62+
63+ @property
64+ def omega(self) -> np.ndarray:
65+ return 2 * math.pi * self.freqs
66+
67+ @property
68+ def sigma(self) -> np.ndarray:
69+ """Total decay rate of each mode [1/s], the sum of the three mechanisms, kept below the mode's frequency (underdamped)."""
70+ return np.minimum(self.sigma_material + self.sigma_thermal + self.sigma_radiation, 0.99 * self.omega)
71+
72+ @property
73+ def t60(self) -> np.ndarray:
74+ """Decay time of each mode [s]: the time in which it falls by 60 dB."""
75+ with np.errstate(divide="ignore"):
76+ return LN1000 / self.sigma
77+
78+ def t60_by_mechanism(self) -> dict[str, np.ndarray]:
79+ """The decay time each mechanism alone would give: ``material``, ``thermal``, ``radiation``."""
80+ out = {}
81+ for name, s in (("material", self.sigma_material), ("thermal", self.sigma_thermal), ("radiation", self.sigma_radiation)):
82+ with np.errstate(divide="ignore"):
83+ out[name] = LN1000 / s
84+ return out
85+
86+ def select(self, index) -> "Modes":
87+ """The modes at the given indices (an integer array, slice or boolean mask)."""
88+ idx = np.arange(self.count)[index]
89+ return replace(
90+ self,
91+ freqs=self.freqs[idx],
92+ w=self.w[idx],
93+ u=None if self.u is None else self.u[idx],
94+ v=None if self.v is None else self.v[idx],
95+ sigma_material=self.sigma_material[idx],
96+ sigma_thermal=self.sigma_thermal[idx],
97+ sigma_radiation=self.sigma_radiation[idx],
98+ )
99+
100+ def below(self, fmax: float) -> "Modes":
101+ """The modes below a frequency [Hz]."""
102+ return self.select(self.freqs < fmax)
103+
104+ def table(self, n: int | None = None) -> str:
105+ """A text table of the first ``n`` modes: frequency and decay times (total, and by mechanism)."""
106+ by = self.t60_by_mechanism()
107+ rows = [" # f [Hz] T60 [s] air metal thermal"]
108+ for m in range(self.count if n is None else min(n, self.count)):
109+ rows.append(f"{m:3d} {self.freqs[m]:9.1f} {self.t60[m]:9.2f} {_fmt(by['radiation'][m])} {_fmt(by['material'][m])} {_fmt(by['thermal'][m])}")
110+ return "\n".join(rows)
111+
112+ def save(self, path) -> None:
113+ """Saves the modes and the plate to an ``.npz`` file (see ``Modes.load``)."""
114+ from ..material import Thermal
115+
116+ pl, m = self.plate, self.plate.material
117+ th = m.thermal or Thermal(0, 1, 1)
118+ np.savez_compressed(
119+ path,
120+ freqs=self.freqs,
121+ w=self.w,
122+ u=np.zeros(0) if self.u is None else self.u,
123+ v=np.zeros(0) if self.v is None else self.v,
124+ sigma_material=self.sigma_material,
125+ sigma_thermal=self.sigma_thermal,
126+ sigma_radiation=self.sigma_radiation,
127+ fmax=self.fmax,
128+ seconds=self.seconds,
129+ mask=pl.grid.mask,
130+ hx=pl.grid.hx,
131+ hy=pl.grid.hy,
132+ h=pl.h,
133+ w0=pl.w0,
134+ Lx=pl.Lx,
135+ Ly=pl.Ly,
136+ temperature=pl.temperature,
137+ name=pl.name,
138+ material=np.array([m.E, m.nu, m.rho, m.loss_factor, 1.0 if m.thermal else 0.0, th.alpha, th.cp, th.k]),
139+ )
140+
141+ @classmethod
142+ def load(cls, path) -> "Modes":
143+ from ..grid import Grid
144+ from ..material import Material, Thermal
145+
146+ d = np.load(path)
147+ mat = d["material"]
148+ material = Material(float(mat[0]), float(mat[1]), float(mat[2]), float(mat[3]), Thermal(float(mat[5]), float(mat[6]), float(mat[7])) if mat[4] > 0 else None)
149+ mask = d["mask"]
150+ grid = Grid(mask.shape[1], mask.shape[0], float(d["hx"]), float(d["hy"]), mask)
151+ plate = Plate(grid, material, d["h"], d["w0"], float(d["Lx"]), float(d["Ly"]), float(d["temperature"]), str(d["name"]))
152+ u, v = d["u"], d["v"]
153+ return cls(plate, d["freqs"], d["w"], None if u.size == 0 else u, None if v.size == 0 else v, d["sigma_material"], d["sigma_thermal"], d["sigma_radiation"], float(d["fmax"]), float(d["seconds"]))
154+
155+ def __repr__(self) -> str:
156+ return f"Modes({self.count} modes of {self.plate.name!r} up to {self.fmax:.0f} Hz)"
157+
158+
159+def _fmt(t: float) -> str:
160+ if not np.isfinite(t) or t <= 0:
161+ return " inf"
162+ if t < 10:
163+ return f"{t:7.2f}"
164+ if t < 1000:
165+ return f"{t:7.0f}"
166+ return f"{t:7.0e}"
167+
168+
169+def compute_modes(
170+ plate: Plate,
171+ fmax: float = 16000.0,
172+ *,
173+ fmin: float = 1.0,
174+ solver: Solver = "auto",
175+ device=None,
176+ damping: bool = True,
177+ system: LinearSystem | None = None,
178+) -> Modes:
179+ """All modes of the plate with frequencies in ``(fmin, fmax]``.
180+
181+ ``solver`` picks the eigensolver: ``"dense"`` (SciPy on the full matrix,
182+ eigenvalues in the band only; the default on the CPU), ``"banded"``
183+ (LAPACK's banded driver, which uses less memory but is slower here),
184+ ``"torch"`` (a dense eigendecomposition on ``device``, useful on a CUDA
185+ GPU), or ``"auto"``, which takes the torch path when a CUDA device is
186+ available and the dense path otherwise. With ``damping`` the per-mode
187+ thermoelastic and radiation rates are computed from the shapes (a second or
188+ two for a few thousand modes); without it they are left at zero.
189+ """
190+ t0 = time.perf_counter()
191+ sys = system or assemble(plate)
192+ lo = (2 * math.pi * fmin) ** 2
193+ hi = (2 * math.pi * fmax) ** 2
194+ if solver == "auto":
195+ solver = "torch" if _cuda_available() else "dense"
196+ if solver == "banded":
197+ lam, y = _eig_banded(sys, lo, hi)
198+ elif solver == "dense":
199+ lam, y = _eig_dense(sys, lo, hi)
200+ elif solver == "torch":
201+ lam, y = _eig_torch(sys, lo, hi, device)
202+ else:
203+ raise ValueError(f"unknown solver {solver!r}")
204+ order = np.argsort(lam)
205+ lam, y = lam[order], y[:, order]
206+ phi = (y / np.sqrt(sys.mass)[:, None]).T # (M, n), mass-normalized
207+ u, v, w = sys.to_grid(phi)
208+ freqs = np.sqrt(np.maximum(lam, 0)) / (2 * math.pi)
209+ M = len(freqs)
210+ modes = Modes(
211+ plate,
212+ freqs,
213+ w,
214+ u if sys.fields == 3 else None,
215+ v if sys.fields == 3 else None,
216+ plate.material.loss_factor * (2 * math.pi * freqs) / 2,
217+ np.zeros(M),
218+ np.zeros(M),
219+ fmax,
220+ )
221+ if damping and M > 0:
222+ from .damping import radiation_rates, thermoelastic_rates
223+
224+ modes.sigma_radiation = radiation_rates(modes)
225+ modes.sigma_thermal = thermoelastic_rates(modes)
226+ modes.seconds = time.perf_counter() - t0
227+ return modes
228+
229+
230+def _cuda_available() -> bool:
231+ try:
232+ import torch
233+
234+ return torch.cuda.is_available()
235+ except ImportError:
236+ return False
237+
238+
239+def _scaled(sys: LinearSystem):
240+ s = 1 / np.sqrt(sys.mass)
241+ A = sys.K.multiply(s[:, None]).multiply(s[None, :]).tocoo()
242+ return A, s
243+
244+
245+def _eig_banded(sys: LinearSystem, lo: float, hi: float):
246+ A, _ = _scaled(sys)
247+ sel = A.row >= A.col
248+ r, c, v = A.row[sel], A.col[sel], A.data[sel]
249+ bw = int((r - c).max()) if len(r) else 0
250+ ab = np.zeros((bw + 1, sys.n))
251+ ab[r - c, c] = v
252+ lam, y = scipy.linalg.eig_banded(ab, lower=True, select="v", select_range=(lo, hi), check_finite=False)
253+ return lam, y
254+
255+
256+def _eig_dense(sys: LinearSystem, lo: float, hi: float):
257+ A, _ = _scaled(sys)
258+ lam, y = scipy.linalg.eigh(A.toarray(), subset_by_value=(lo, hi), driver="evr", check_finite=False)
259+ return lam, y
260+
261+
262+def _eig_torch(sys: LinearSystem, lo: float, hi: float, device):
263+ import torch
264+
265+ A, _ = _scaled(sys)
266+ dev = torch.device(device) if device is not None else torch.device("cuda" if torch.cuda.is_available() else "cpu")
267+ At = torch.from_numpy(A.toarray()).to(dev, dtype=torch.float64)
268+ lam, y = torch.linalg.eigh(At)
269+ lam = lam.cpu().numpy()
270+ sel = (lam > lo) & (lam <= hi)
271+ return lam[sel], y[:, torch.from_numpy(sel).to(dev)].cpu().numpy()
src/cymsim/linear/synth.pyadded+232−0View file
@@ -0,0 +1,232 @@
1+"""Modal synthesis of a strike. Each mode m of the linearized plate is a damped
2+oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma the sum of its three
3+loss rates. The stick is the Hertz contact of the full model, stepped against
4+the modal system with an exact update of every oscillator under a piecewise
5+constant force. Once the stick has left, every mode is a damped sinusoid with a
6+known complex amplitude, and the delayed Rayleigh sum over the plate collapses
7+to one complex gain per mode. During the attack, while the modes are being
8+driven, the delayed sum is evaluated directly from a table of the nodal
9+accelerations.
10+"""
11+
12+from __future__ import annotations
13+
14+import math
15+from dataclasses import dataclass
16+
17+import numpy as np
18+
19+from ..constants import STICK_RELEASE
20+from ..params import Listener, Strike
21+from ..plate import Plate
22+from ..sound import Sound
23+from .modes import Modes
24+
25+
26+def strike_gains(modes: Modes, strike: Strike) -> np.ndarray:
27+ """The contact: modal gains g_m = sum over the patch of (m_k / M_patch) phi_m(k),
28+ so that the modal force is -F g_m and the contact coordinate is sum_m g_m q_m."""
29+ pl = modes.plate
30+ patch = pl.strike_patch(strike)
31+ frac = pl.node_mass[patch.mask] / patch.mass
32+ return modes.w[:, patch.mask] @ frac
33+
34+
35+@dataclass
36+class Contact:
37+ """The state of every mode at the end of the contact phase and the tables recorded during it."""
38+
39+ t_end: float
40+ """Time at which the stick is withdrawn, STICK_RELEASE after the contact last opened [s]."""
41+ t_contact: float
42+ """Time at which the contact last opened [s]; every mode decays freely from here on."""
43+ amplitude: np.ndarray
44+ """Complex modal amplitudes A at ``t_end``: q(t) = Re[A exp(s (t - t_end))], s = -sigma + i omega_d."""
45+ qdd: np.ndarray
46+ """Modal accelerations sampled at ``table_rate`` from t = 0 to t_end, ``(table_len, M)``."""
47+ table_rate: float
48+ q_frames: np.ndarray
49+ """Modal displacements at the movie frame times before ``t_end``, ``(n_frames, M)``."""
50+ force: np.ndarray
51+ """Stick force history at the fine step [N]."""
52+ dt: float
53+ peak_force: float
54+ contact_time: float
55+ """Total time the contact was closed [s]."""
56+
57+
58+def simulate_contact(modes: Modes, sigma: np.ndarray, g: np.ndarray, strike: Strike, fps: float, fs: float, dt: float = 2e-6, table_rate: float | None = None) -> Contact:
59+ """Steps the stick against the modal system until the contact has been open
60+ for STICK_RELEASE, recording the tables the synthesis needs. Every mode is
61+ advanced by its exact state transition under a constant force over the step;
62+ the stick by symplectic Euler; the force is the Hertz law of the penetration
63+ at the start of the step."""
64+ table_rate = table_rate or 2 * fs
65+ omega = modes.omega
66+ wd = np.sqrt(omega**2 - sigma**2)
67+ e, c, s = np.exp(-sigma * dt), np.cos(wd * dt), np.sin(wd * dt)
68+ p11 = e * (c + sigma * s / wd)
69+ p12 = e * s / wd
70+ p21 = -e * omega**2 * s / wd
71+ p22 = e * (c - sigma * s / wd)
72+ M = modes.count
73+ q, qd = np.zeros(M), np.zeros(M)
74+ t_max = 0.2
75+ max_steps = math.ceil(t_max / dt)
76+ table_cap = math.ceil(t_max * table_rate) + 2
77+ frame_cap = math.ceil(t_max * fps) + 2
78+ qdd = np.zeros((table_cap, M))
79+ q_frames = np.zeros((frame_cap, M))
80+ force = np.zeros(max_steps + 1)
81+ table_len = n_frames = 0
82+ peak_force = contact_time = 0.0
83+ z, zd = 0.0, -strike.velocity # stick tip, positive away from the plate; it arrives moving toward it
84+ last_open, ever, F = 0.0, False, 0.0
85+ step = 0
86+ for step in range(max_steps):
87+ t = step * dt
88+ delta = float(g @ q) - z # penetration: the plate has come past the tip
89+ F = strike.stiffness * delta**strike.exponent if delta > 0 else 0.0
90+ if F > 0:
91+ ever = True
92+ last_open = t + dt
93+ contact_time += dt
94+ peak_force = max(peak_force, F)
95+ force[step] = F
96+ while table_len < table_cap and table_len / table_rate <= t + 1e-12:
97+ qdd[table_len] = -2 * sigma * qd - omega**2 * q - F * g
98+ table_len += 1
99+ while n_frames < frame_cap and n_frames / fps <= t + 1e-12:
100+ q_frames[n_frames] = q
101+ n_frames += 1
102+ if ever and t - last_open >= STICK_RELEASE:
103+ break
104+ qp = -F * g / omega**2
105+ x = q - qp
106+ q = p11 * x + p12 * qd + qp
107+ qd = p21 * x + p22 * qd
108+ zd += F / strike.mass * dt
109+ z += zd * dt
110+ else:
111+ step = max_steps
112+ t_end = step * dt
113+ amplitude = q - 1j * (qd + sigma * q) / wd
114+ return Contact(t_end, min(last_open, t_end), amplitude, qdd[:table_len], table_rate, q_frames[:n_frames], force[: step + 1], dt, peak_force, contact_time)
115+
116+
117+def listener_weights(pl: Plate, listener: Listener) -> tuple[np.ndarray, np.ndarray, float]:
118+ """Listener terms with the common delay removed, so that only the spread of the delays across the plate remains."""
119+ weight, delay = pl.listener_terms(listener)
120+ sel = weight > 0
121+ delay = np.where(sel, delay - delay[sel].min(), 0.0)
122+ return weight, delay, float(delay.max())
123+
124+
125+def tail_start(c: Contact, max_delay: float) -> float:
126+ """The time from which the tail formula holds at the listener: every node's
127+ retarded time is past the last contact. The table path covers the output before it."""
128+ return min(c.t_end, c.t_contact + max_delay + 2 / c.table_rate)
129+
130+
131+def radiation_gains(modes: Modes, sigma: np.ndarray, weight: np.ndarray, delay: np.ndarray, batch: int = 256) -> np.ndarray:
132+ """The complex radiation gain of each mode for the listener:
133+ G_m = sum_k weight_k phi_m(k) exp(-s_m delay_k), with s_m = -sigma + i omega_d."""
134+ sel = weight > 0
135+ w, tau = weight[sel], delay[sel]
136+ wd = np.sqrt(modes.omega**2 - sigma**2)
137+ G = np.zeros(modes.count, dtype=complex)
138+ for a in range(0, modes.count, batch):
139+ b = min(modes.count, a + batch)
140+ phi = modes.w[a:b][:, sel]
141+ E = np.exp((sigma[a:b, None] - 1j * wd[a:b, None]) * tau[None, :])
142+ G[a:b] = (phi * w[None, :] * E).sum(axis=1)
143+ return G
144+
145+
146+def attack_pressure(modes: Modes, c: Contact, weight: np.ndarray, delay: np.ndarray, fs: float, t_start: float, out: np.ndarray) -> None:
147+ """The listener pressure over the attack from the nodal acceleration table:
148+ p(t) = sum_k weight_k wdd_k(t - delay_k), linearly interpolated, added to
149+ ``out`` for the samples before ``t_start``."""
150+ n_out = min(len(out), math.ceil(t_start * fs - 1e-9))
151+ if n_out <= 0:
152+ return
153+ length = min(len(c.qdd), math.ceil(t_start * c.table_rate) + 2)
154+ sel = weight > 0
155+ table = modes.w[:, sel].T @ c.qdd[:length].T # (R, length) nodal accelerations
156+ x = (np.arange(n_out)[None, :] / fs - delay[sel][:, None]) * c.table_rate # (R, n_out)
157+ i = np.floor(x).astype(int)
158+ valid = (x >= 0) & (i + 1 < length)
159+ i0 = np.clip(i, 0, length - 1)
160+ i1 = np.clip(i + 1, 0, length - 1)
161+ fr = x - i
162+ rows = np.arange(sel.sum())[:, None]
163+ val = (1 - fr) * table[rows, i0] + fr * table[rows, i1]
164+ out[:n_out] += (weight[sel][:, None] * np.where(valid, val, 0.0)).sum(axis=0)
165+
166+
167+def modal_tail(modes: Modes, sigma: np.ndarray, c: Contact, gains: np.ndarray, fs: float, t_start: float, out: np.ndarray, block: int = 4096) -> None:
168+ """The listener pressure from ``t_start`` on: p(t) = sum_m Re[G_m s_m^2 A_m
169+ exp(s_m (t - t_end))], added to ``out`` in blocks. Modes that have decayed
170+ by 140 dB are dropped."""
171+ n = len(out)
172+ n0 = math.ceil(t_start * fs - 1e-9)
173+ if n0 >= n:
174+ return
175+ wd = np.sqrt(modes.omega**2 - sigma**2)
176+ s = -sigma + 1j * wd
177+ coef = gains * s**2 * c.amplitude * np.exp(s * (n0 / fs - c.t_end)) # the value at sample n0
178+ Z = np.exp(np.outer(np.arange(block) / fs, s)) # (block, M): z^b
179+ scale = np.abs(coef)
180+ for start in range(n0, n, block):
181+ stop = min(n, start + block)
182+ cur = coef * np.exp(s * ((start - n0) / fs))
183+ active = np.abs(cur) > 1e-7 * scale
184+ if not active.any():
185+ break
186+ out[start:stop] += (Z[: stop - start, active] @ cur[active]).real
187+
188+
189+def movie_frames(modes: Modes, sigma: np.ndarray, c: Contact, fps: float, n_frames: int) -> np.ndarray:
190+ """Movie frames of the transverse displacement, ``(n_frames, ny, nx)`` float32:
191+ sum_m phi_m(k) q_m(t_f), with q from the contact tables before t_end and from the free decay after it."""
192+ wd = np.sqrt(modes.omega**2 - sigma**2)
193+ t = np.arange(n_frames) / fps - c.t_end
194+ Q = np.exp(-sigma[None, :] * t[:, None]) * (c.amplitude.real[None, :] * np.cos(wd[None, :] * t[:, None]) - c.amplitude.imag[None, :] * np.sin(wd[None, :] * t[:, None]))
195+ k = min(len(c.q_frames), n_frames)
196+ Q[:k] = c.q_frames[:k]
197+ g = modes.plate.grid
198+ sel = g.mask.reshape(-1)
199+ frames = np.zeros((n_frames, g.n), dtype=np.float32)
200+ frames[:, sel] = Q @ modes.w.reshape(modes.count, -1)[:, sel]
201+ return frames.reshape(n_frames, g.ny, g.nx)
202+
203+
204+def synthesize(modes: Modes, strike: Strike = Strike(), listener: Listener = Listener(), *, duration: float = 4.0, fs: int = 48000, fps: float = 120.0, movie: bool = False) -> Sound:
205+ """The sound of a strike on the plate by modal synthesis.
206+
207+ Returns a ``Sound`` with the listener pressure [Pa] at ``fs`` over
208+ ``duration`` seconds after the impact and, with ``movie``, the transverse
209+ displacement at ``fps`` frames per second. The strike and the listener can
210+ be changed at no cost beyond this call; the modes are the expensive part.
211+ """
212+ pl = modes.plate
213+ sigma = modes.sigma
214+ g = strike_gains(modes, strike)
215+ c = simulate_contact(modes, sigma, g, strike, fps, fs)
216+ weight, delay, max_delay = listener_weights(pl, listener)
217+ t_start = tail_start(c, max_delay)
218+ n = round(duration * fs)
219+ audio = np.zeros(n)
220+ attack_pressure(modes, c, weight, delay, fs, t_start, audio)
221+ modal_tail(modes, sigma, c, radiation_gains(modes, sigma, weight, delay), fs, t_start, audio)
222+ frames = movie_frames(modes, sigma, c, fps, round(duration * fps)) if movie else None
223+ info = {
224+ "model": "linear",
225+ "modes": modes.count,
226+ "peak_force": c.peak_force,
227+ "contact_time": c.contact_time,
228+ "t_end": c.t_end,
229+ "force": c.force,
230+ "force_dt": c.dt,
231+ }
232+ return Sound(audio.astype(np.float32), fs, pl, strike, listener, frames, fps, info)
src/cymsim/loss.pyadded+274−0View file
@@ -0,0 +1,274 @@
1+"""The losses of the plate, derived from physical constants rather than prescribed.
2+
3+Three mechanisms drain the energy of a vibrating plate in this model. Internal
4+friction of the metal is a hysteretic loss with the material's loss factor eta:
5+a mode of angular frequency omega decays at the rate eta omega / 2. Thermoelastic
6+damping (Zener) comes from the heat that flows across the thickness between the
7+compressed and the stretched face of a bent plate; it follows from the thermal
8+constants of the metal and the local thickness, with a peak at the frequency set
9+by the thermal diffusion time across the thickness. Radiation into the air is
10+the energy the plate loses as sound; the linear model computes it exactly for
11+each mode from the mode's shape (``linear.damping``), and this module gives the
12+analytic average of Leppington for a plate of given area and perimeter, which
13+the full model's two-parameter loss operator is fitted to.
14+"""
15+
16+from __future__ import annotations
17+
18+import math
19+from dataclasses import dataclass
20+
21+import numpy as np
22+
23+from .constants import LN1000, RHO_AIR, SPEED_OF_SOUND
24+from .grid import Grid
25+from .material import Material
26+from .plate import Plate
27+
28+__all__ = [
29+ "LN1000",
30+ "material_rate",
31+ "thermoelastic_relaxation",
32+ "thermoelastic_loss_factor",
33+ "outline_perimeter",
34+ "coincidence_frequency",
35+ "average_radiation_efficiency",
36+ "radiation_rate_from_efficiency",
37+ "PlateGeometry",
38+ "plate_geometry",
39+ "LossEstimate",
40+ "estimated_rate",
41+ "in_plane_rate",
42+ "SigmaFit",
43+ "fit_sigma",
44+ "fit_line",
45+ "fitted_rate",
46+ "t60",
47+]
48+
49+
50+def t60(sigma):
51+ """Decay time [s] of a rate sigma [1/s]: the time in which a partial falls by 60 dB."""
52+ with np.errstate(divide="ignore"):
53+ return LN1000 / np.asarray(sigma, dtype=float)
54+
55+
56+# ---------- Internal friction ----------
57+
58+
59+def material_rate(eta: float, omega):
60+ """Decay rate of a mode of angular frequency omega from a hysteretic loss
61+ factor eta: the energy lost per cycle is 2 pi eta times the stored energy."""
62+ return eta * omega / 2
63+
64+
65+# ---------- Thermoelastic damping ----------
66+
67+
68+def thermoelastic_relaxation(mat: Material, h, T: float):
69+ """Zener's relaxation strength ``delta`` and time ``tau`` [s] of a plate of
70+ thickness ``h``. The strength carries the factor (1 + nu) / (1 - nu) of
71+ the biaxial stress state of a bent plate (Zener's beam value is
72+ E alpha^2 T / (rho cp)); the time is the thermal diffusion time across the
73+ thickness, h^2 / (pi^2 chi) with the diffusivity chi = k / (rho cp)."""
74+ th = mat.thermal
75+ if th is None:
76+ return 0.0, np.ones_like(np.asarray(h, dtype=float))
77+ delta = (mat.E * th.alpha**2 * T) / (mat.rho * th.cp) * ((1 + mat.nu) / (1 - mat.nu))
78+ tau = np.asarray(h, dtype=float) ** 2 * mat.rho * th.cp / (math.pi**2 * th.k)
79+ return delta, tau
80+
81+
82+def thermoelastic_loss_factor(mat: Material, h, omega, T: float):
83+ """Thermoelastic loss factor of bending at angular frequency ``omega`` in a
84+ plate of thickness ``h``: delta x / (1 + x^2) with x = omega tau, largest
85+ (delta / 2) at x = 1. ``h`` and ``omega`` broadcast."""
86+ delta, tau = thermoelastic_relaxation(mat, h, T)
87+ x = np.asarray(omega, dtype=float) * tau
88+ return delta * x / (1 + x * x)
89+
90+
91+# ---------- Geometry of the outline ----------
92+
93+_DIAG = 0
94+_HORIZ = 1
95+_VERT = 2
96+# Segments of the marching-squares contour by code, as counts of (diagonal, horizontal, vertical) pieces.
97+_SEGMENTS = {1: (1, 0, 0), 14: (1, 0, 0), 2: (1, 0, 0), 13: (1, 0, 0), 3: (0, 1, 0), 12: (0, 1, 0), 4: (1, 0, 0), 11: (1, 0, 0), 6: (0, 0, 1), 9: (0, 0, 1), 7: (1, 0, 0), 8: (1, 0, 0), 5: (2, 0, 0), 10: (2, 0, 0)}
98+
99+
100+def outline_perimeter(g: Grid) -> float:
101+ """Length [m] of the outline of the plate, from the marching-squares contour
102+ of the cell mask at level one half (cells outside the grid count as outside).
103+ A staircase count would overstate the perimeter of a disc by 4 / pi; the
104+ contour, whose segments run in eight directions, overstates it by a few
105+ percent, which is well within the accuracy of the average radiation
106+ efficiency it feeds."""
107+ c = np.pad(g.cell_mask, 1).astype(np.int64)
108+ a, b, cc, d = c[:-1, :-1], c[:-1, 1:], c[1:, 1:], c[1:, :-1]
109+ code = a | (b << 1) | (cc << 2) | (d << 3)
110+ counts = np.bincount(code.ravel(), minlength=16)
111+ diag = math.hypot(0.5 * g.hx, 0.5 * g.hy)
112+ length = 0.0
113+ for k, (nd, nh, nv) in _SEGMENTS.items():
114+ length += counts[k] * (nd * diag + nh * g.hx + nv * g.hy)
115+ return float(length)
116+
117+
118+# ---------- Radiation ----------
119+
120+
121+def coincidence_frequency(mat: Material, h) -> float:
122+ """Frequency at which the bending wavelength equals the wavelength of sound, c^2 / (2 pi kappa) [Hz]."""
123+ return mat.coincidence_frequency(h)
124+
125+
126+def average_radiation_efficiency(f, fc: float, area: float, perimeter: float, a: float, b: float, free_edges: bool = True):
127+ """The radiation efficiency of a baffled plate averaged over its modes near
128+ the frequency ``f``. The simply supported case is that of Leppington,
129+ Broadbent and Heron (1982): below coincidence the edges radiate and the
130+ efficiency grows about linearly with frequency, at coincidence it peaks
131+ at a value set by the size of the plate in coincidence wavelengths, and
132+ above it tends to one. A plate with free edges radiates much less below
133+ coincidence (its modes have no monopole or dipole moment, and at a free
134+ edge the bending moment and the shear vanish); the per-mode computation
135+ on the presets follows the simply supported curve times (f / fc)^2 within
136+ a factor of two from 100 Hz to 8 kHz, and that is what ``free_edges``
137+ applies. ``a`` and ``b`` are the shorter and the longer side."""
138+ f = np.asarray(f, dtype=float)
139+ c = SPEED_OF_SOUND
140+ lambda_c = c / fc
141+ cap = max(1.0, (0.5 - 0.15 * (a / b)) * math.sqrt(a / lambda_c))
142+ above = f >= fc
143+ fb = np.where(above, 0.5 * fc, f) # safe stand-in where the other branch applies
144+ fa = np.where(above, f, 2 * fc)
145+ with np.errstate(divide="ignore", invalid="ignore"):
146+ s_above = np.minimum(cap, 1 / np.sqrt(np.maximum(1e-12, 1 - fc / fa)))
147+ mu = np.sqrt(fc / fb)
148+ s_below = (perimeter * c) / (4 * math.pi**2 * area * fc) / np.sqrt(mu * mu - 1) * (np.log((mu + 1) / (mu - 1)) + 2 * mu / (mu * mu - 1))
149+ s_below = np.minimum(cap, s_below) * ((fb / fc) ** 2 if free_edges else 1.0)
150+ out = np.where(above, s_above, s_below)
151+ return float(out) if out.ndim == 0 else out
152+
153+
154+def radiation_rate_from_efficiency(efficiency, mat: Material, h):
155+ """Decay rate from a radiation efficiency: the two faces of a baffled plate
156+ of surface mass rho h radiate the power 2 rho_air c sigma_rad <v^2> A, so
157+ sigma = rho_air c sigma_rad / (rho h)."""
158+ return RHO_AIR * SPEED_OF_SOUND * efficiency / (mat.rho * h)
159+
160+
161+# ---------- The estimate for a whole plate ----------
162+
163+
164+@dataclass(frozen=True)
165+class PlateGeometry:
166+ """What the loss estimate needs of a plate."""
167+
168+ area: float
169+ perimeter: float
170+ h_mean: float
171+ kappa: float
172+ fc: float
173+ a: float
174+ """Shorter side of the bounding box [m]."""
175+ b: float
176+ """Longer side [m]."""
177+
178+
179+def plate_geometry(pl: Plate) -> PlateGeometry:
180+ hm = pl.mean_thickness
181+ return PlateGeometry(pl.area, outline_perimeter(pl.grid), hm, float(pl.material.bending_constant(hm)), coincidence_frequency(pl.material, hm), min(pl.Lx, pl.Ly), max(pl.Lx, pl.Ly))
182+
183+
184+@dataclass(frozen=True)
185+class LossEstimate:
186+ """The three contributions to the decay rate of a bending partial [1/s], and their sum."""
187+
188+ material: float
189+ thermoelastic: float
190+ radiation: float
191+ total: float
192+
193+
194+def estimated_rate(pl: Plate, geom: PlateGeometry, f: float) -> LossEstimate:
195+ """The decay rate of a bending partial at frequency ``f`` estimated without
196+ the modes: the material rate, the thermoelastic rate averaged over the
197+ plate area with the local thickness, and the radiation rate from the
198+ average efficiency at the mean thickness. The modal computation replaces
199+ the last two by exact per-mode values; this estimate serves the full
200+ model's loss operator and the summary of a design."""
201+ omega = 2 * math.pi * f
202+ material = material_rate(pl.material.loss_factor, omega)
203+ te = 0.0
204+ if pl.material.thermal is not None:
205+ area = pl.grid.node_area
206+ sel = area > 0
207+ eta = thermoelastic_loss_factor(pl.material, pl.h[sel], omega, pl.temperature)
208+ te = float((area[sel] * eta).sum() / area[sel].sum()) * omega / 2
209+ radiation = radiation_rate_from_efficiency(average_radiation_efficiency(f, geom.fc, geom.area, geom.perimeter, geom.a, geom.b), pl.material, geom.h_mean)
210+ return LossEstimate(material, te, float(radiation), material + te + float(radiation))
211+
212+
213+def in_plane_rate(pl: Plate) -> float:
214+ """The in-plane motion of the full model has no radiation and no thermoelastic
215+ loss in this model; it decays by internal friction at the rate of the lowest
216+ in-plane resonance, c_p / (2 L) of the longer side."""
217+ f_ip = pl.material.in_plane_speed / (2 * max(pl.Lx, pl.Ly))
218+ return material_rate(pl.material.loss_factor, 2 * math.pi * f_ip)
219+
220+
221+# ---------- The two-parameter operator of the full model ----------
222+
223+
224+@dataclass(frozen=True)
225+class SigmaFit:
226+ """The two-parameter loss model of the full model: a bending mode of angular
227+ frequency omega decays at the rate sigma0 + sigma1 k^2 with k^2 = omega / kappa."""
228+
229+ sigma0: float
230+ """[1/s]"""
231+ sigma1: float
232+ """[m^2/s]"""
233+ kappa: float
234+ """Bending constant used for the fit [m^2/s]."""
235+
236+
237+def fit_sigma(pl: Plate, geom: PlateGeometry | None = None, f_lo: float = 50.0, f_hi: float = 10000.0, n: int = 40) -> SigmaFit:
238+ """Non-negative least-squares fit of sigma0 + sigma1 k^2 to the estimated
239+ rates on ``n`` log-spaced frequencies from ``f_lo`` to ``f_hi``, at the mean thickness."""
240+ geom = geom or plate_geometry(pl)
241+ fs = f_lo * (f_hi / f_lo) ** (np.arange(n) / (n - 1))
242+ x = 2 * math.pi * fs / geom.kappa
243+ y = np.array([estimated_rate(pl, geom, f).total for f in fs])
244+ return fit_line(x, y, geom.kappa)
245+
246+
247+def fit_line(x, y, kappa: float) -> SigmaFit:
248+ """Non-negative least squares of y = s0 + s1 x: with two parameters the
249+ solution is the unconstrained one when it is feasible, else the best of
250+ the fits with one parameter at zero, else zero."""
251+ x = np.asarray(x, float)
252+ y = np.asarray(y, float)
253+ n = len(x)
254+ sx, sy, sxx, sxy = x.sum(), y.sum(), (x * x).sum(), (x * y).sum()
255+ det = n * sxx - sx * sx
256+ candidates = [(0.0, 0.0), (sy / n, 0.0)]
257+ if det > 0:
258+ s1 = (n * sxy - sx * sy) / det
259+ candidates.append(((sy - s1 * sx) / n, s1))
260+ if sxx > 0:
261+ candidates.append((0.0, sxy / sxx))
262+ best, best_res = (0.0, 0.0), math.inf
263+ for s0, s1 in candidates:
264+ if not (s0 >= 0 and s1 >= 0):
265+ continue
266+ res = float(((y - s0 - s1 * x) ** 2).sum())
267+ if res < best_res:
268+ best, best_res = (s0, s1), res
269+ return SigmaFit(best[0], best[1], kappa)
270+
271+
272+def fitted_rate(fit: SigmaFit, f):
273+ """Decay rate [1/s] of the fitted model at a frequency."""
274+ return fit.sigma0 + fit.sigma1 * 2 * math.pi * np.asarray(f, float) / fit.kappa
src/cymsim/material.pyadded+96−0View file
@@ -0,0 +1,96 @@
1+"""Materials: elastic constants, the loss factor of internal friction, and the
2+thermal constants from which the thermoelastic damping follows."""
3+
4+from __future__ import annotations
5+
6+import math
7+from dataclasses import dataclass
8+
9+from .constants import SPEED_OF_SOUND
10+
11+
12+@dataclass(frozen=True)
13+class Thermal:
14+ """Thermal constants of a metal. All three are needed for Zener's thermoelastic damping."""
15+
16+ alpha: float
17+ """Linear thermal expansion [1/K]."""
18+ cp: float
19+ """Specific heat [J/(kg K)]."""
20+ k: float
21+ """Thermal conductivity [W/(m K)]."""
22+
23+ def __post_init__(self) -> None:
24+ if not (self.alpha >= 0 and self.cp >= 1 and self.k >= 1e-3):
25+ raise ValueError("thermal constants must be positive (alpha >= 0, cp >= 1, k >= 1e-3)")
26+
27+
28+@dataclass(frozen=True)
29+class Material:
30+ """A homogeneous isotropic metal.
31+
32+ The decay of the sound is not prescribed but derived from ``loss_factor``
33+ (the hysteretic loss factor of the metal's internal friction, of order
34+ 1e-4 for bronze, steel and aluminum), from ``thermal`` (Zener's
35+ thermoelastic damping; ``None`` leaves it out) and from the radiation of
36+ the plate into the air, which the models compute from the geometry.
37+ """
38+
39+ E: float
40+ """Young's modulus [Pa]."""
41+ nu: float
42+ """Poisson ratio."""
43+ rho: float
44+ """Density [kg/m^3]."""
45+ loss_factor: float = 0.0
46+ """Loss factor eta of internal friction: a mode decays at the rate eta omega / 2."""
47+ thermal: Thermal | None = None
48+ """Thermal constants for the thermoelastic damping, or None for none."""
49+
50+ def __post_init__(self) -> None:
51+ if not self.E >= 1e6:
52+ raise ValueError("E must be at least 1e6 Pa")
53+ if not 0 <= self.nu < 0.5:
54+ raise ValueError("nu must lie in [0, 0.5)")
55+ if not self.rho >= 1:
56+ raise ValueError("rho must be at least 1 kg/m^3")
57+ if not 0 <= self.loss_factor < 0.1:
58+ raise ValueError("loss_factor must be small (metals have 1e-4 to 1e-3)")
59+
60+ def bending_stiffness(self, h):
61+ """D = E h^3 / (12 (1 - nu^2)) [N m]; ``h`` may be an array."""
62+ return self.E * h**3 / (12 * (1 - self.nu**2))
63+
64+ def membrane_stiffness(self, h):
65+ """C = E h / (1 - nu^2) [N/m]."""
66+ return self.E * h / (1 - self.nu**2)
67+
68+ def shear_stiffness(self, h):
69+ """G = E h / (2 (1 + nu)) [N/m]."""
70+ return self.E * h / (2 * (1 + self.nu))
71+
72+ @property
73+ def in_plane_speed(self) -> float:
74+ """Speed of in-plane waves, sqrt(E / (rho (1 - nu^2))) [m/s]; it sets the time step of the full model."""
75+ return math.sqrt(self.E / (self.rho * (1 - self.nu**2)))
76+
77+ def bending_constant(self, h):
78+ """kappa = sqrt(D / (rho h)) [m^2/s], so that omega = kappa k^2 for a bending wave."""
79+ return (self.bending_stiffness(h) / (self.rho * h)) ** 0.5
80+
81+ def coincidence_frequency(self, h):
82+ """Frequency at which the bending wavelength equals the wavelength of sound, c^2 / (2 pi kappa) [Hz]."""
83+ return SPEED_OF_SOUND**2 / (2 * math.pi * self.bending_constant(h))
84+
85+
86+STAINLESS_STEEL = Material(E=170e9, nu=0.29, rho=7900, loss_factor=3e-4, thermal=Thermal(17e-6, 500, 16))
87+"""Stainless steel sheet."""
88+
89+MILD_STEEL = Material(E=200e9, nu=0.3, rho=7850, loss_factor=3e-4, thermal=Thermal(12e-6, 460, 50))
90+"""Carbon steel sheet."""
91+
92+BRONZE = Material(E=110e9, nu=0.34, rho=8800, loss_factor=2e-4, thermal=Thermal(18e-6, 370, 40))
93+"""Bell bronze, as used for cymbals."""
94+
95+ALUMINUM = Material(E=68.9e9, nu=0.33, rho=2700, loss_factor=1e-4, thermal=Thermal(23.6e-6, 896, 167))
96+"""6061-T6 aluminum."""
src/cymsim/params.pyadded+64−0View file
@@ -0,0 +1,64 @@
1+"""The strike and the listener."""
2+
3+from __future__ import annotations
4+
5+import math
6+from dataclasses import dataclass
7+
8+
9+@dataclass(frozen=True)
10+class Strike:
11+ """A drumstick strike: a free mass arriving at a point of the plate with a Hertz contact.
12+
13+ The default contact stiffness is softer than the Hertz value for a hard
14+ wooden tip on metal (of order 1e8 to 1e9); it stands for the compliance of
15+ a real tip and gives a rounder attack.
16+ """
17+
18+ x: float = 0.9
19+ """Strike point as a fraction of the bounding-box width from the left edge."""
20+ y: float = 0.5
21+ """Strike point as a fraction of the bounding-box height from the bottom edge."""
22+ velocity: float = 3.0
23+ """Stick tip speed at impact [m/s]."""
24+ mass: float = 0.03
25+ """Effective mass of the stick [kg]."""
26+ stiffness: float = 1e7
27+ """Hertz contact stiffness K [N/m^alpha]."""
28+ exponent: float = 1.5
29+ """Hertz exponent alpha (1.5 for a rounded tip)."""
30+ radius: float = 5e-3
31+ """Radius of the contact patch [m]."""
32+
33+ def __post_init__(self) -> None:
34+ if not (0 <= self.x <= 1 and 0 <= self.y <= 1):
35+ raise ValueError("the strike point (x, y) is given as fractions of the bounding box in [0, 1]")
36+ if not (self.velocity > 0 and self.mass > 0 and self.stiffness > 0 and self.exponent >= 1 and self.radius >= 0):
37+ raise ValueError("velocity, mass and stiffness must be positive, exponent at least 1, radius non-negative")
38+
39+
40+@dataclass(frozen=True)
41+class Listener:
42+ """The listener (a microphone) in spherical coordinates about the center of the plate."""
43+
44+ distance: float = 0.5
45+ """Distance from the center of the plate [m]."""
46+ angle: float = 30.0
47+ """Angle from the plate's normal [deg]; 0 is straight above the plate."""
48+ azimuth: float = 0.0
49+ """Rotation about the normal [deg]; 0 is toward +x."""
50+
51+ def __post_init__(self) -> None:
52+ if not self.distance > 0:
53+ raise ValueError("distance must be positive")
54+
55+ @property
56+ def position(self) -> tuple[float, float, float]:
57+ """Cartesian position relative to the center of the plate, z along the normal [m]."""
58+ th = math.radians(self.angle)
59+ ph = math.radians(self.azimuth)
60+ return (
61+ self.distance * math.sin(th) * math.cos(ph),
62+ self.distance * math.sin(th) * math.sin(ph),
63+ self.distance * math.cos(th),
64+ )
src/cymsim/plate.pyadded+231−0View file
@@ -0,0 +1,231 @@
1+"""The plate as evaluated on the grid: material, thickness field, stress-free
2+shape, and the derived quantities both models need (masses, the strike patch,
3+the listener terms). A ``Plate`` is usually produced by ``Design.evaluate``;
4+``Plate.uniform`` builds one directly for tests and experiments."""
5+
6+from __future__ import annotations
7+
8+import math
9+from dataclasses import dataclass, replace
10+from functools import cached_property
11+
12+import numpy as np
13+
14+from .constants import RHO_AIR, ROOM_TEMPERATURE, SPEED_OF_SOUND
15+from .grid import Grid
16+from .material import Material
17+from .params import Listener, Strike
18+
19+
20+@dataclass(frozen=True)
21+class Patch:
22+ """The plate nodes that a strike acts on."""
23+
24+ mask: np.ndarray
25+ """``(ny, nx)`` bool: the nodes of the patch."""
26+ nearest: tuple[int, int]
27+ """``(j, i)`` of the plate node nearest the strike point."""
28+ mass: float
29+ """Total lumped mass of the patch [kg]."""
30+ cx: float
31+ """Mass-weighted center of the patch [m]."""
32+ cy: float
33+
34+
35+@dataclass(frozen=True)
36+class Plate:
37+ grid: Grid
38+ material: Material
39+ h: np.ndarray
40+ """Thickness at the nodes [m], ``(ny, nx)``, zero outside the plate."""
41+ w0: np.ndarray
42+ """Stress-free height at the nodes [m], mean plane removed, zero outside the plate."""
43+ Lx: float
44+ """Bounding box [m]."""
45+ Ly: float
46+ temperature: float = ROOM_TEMPERATURE
47+ """Temperature [K], for the thermoelastic damping."""
48+ name: str = "plate"
49+
50+ def __post_init__(self) -> None:
51+ if self.h.shape != self.grid.shape or self.w0.shape != self.grid.shape:
52+ raise ValueError("h and w0 must have the grid's shape (ny, nx)")
53+ if (self.h[self.grid.mask] <= 0).any():
54+ raise ValueError("the thickness must be positive on the plate")
55+
56+ @classmethod
57+ def uniform(cls, grid: Grid, material: Material, thickness: float, profile: np.ndarray | None = None, name: str = "uniform") -> "Plate":
58+ """A plate of uniform thickness on a grid, optionally with a stress-free profile (its mean plane is removed)."""
59+ h = np.where(grid.mask, thickness, 0.0)
60+ w0 = np.zeros(grid.shape) if profile is None else np.where(grid.mask, remove_plane(np.asarray(profile, float), grid), 0.0)
61+ return cls(grid, material, h, w0, (grid.nx - 1) * grid.hx, (grid.ny - 1) * grid.hy, name=name)
62+
63+ def with_thickness(self, h: np.ndarray) -> "Plate":
64+ """The same plate with another thickness field (masked to the plate)."""
65+ return replace(self, h=np.where(self.grid.mask, h, 0.0))
66+
67+ # ---- masses and geometry ----
68+
69+ @cached_property
70+ def node_mass(self) -> np.ndarray:
71+ """Lumped mass of each node, rho h dA times its weight [kg]."""
72+ return self.material.rho * self.h * self.grid.node_area
73+
74+ @cached_property
75+ def mass(self) -> float:
76+ """Total mass [kg]."""
77+ return float(self.node_mass.sum())
78+
79+ @cached_property
80+ def center(self) -> tuple[float, float]:
81+ """Mass-weighted center [m]."""
82+ m = self.node_mass
83+ return float((m * self.grid.X).sum() / self.mass), float((m * self.grid.Y).sum() / self.mass)
84+
85+ @cached_property
86+ def inertia(self) -> tuple[float, float]:
87+ """Moments of inertia ``(Ix, Iy)`` about the in-plane axes through the center [kg m^2]."""
88+ cx, cy = self.center
89+ m = self.node_mass
90+ return float((m * (self.grid.Y - cy) ** 2).sum()), float((m * (self.grid.X - cx) ** 2).sum())
91+
92+ @property
93+ def area(self) -> float:
94+ return self.grid.area
95+
96+ @cached_property
97+ def mean_thickness(self) -> float:
98+ """Area-weighted mean thickness [m]."""
99+ w = self.grid.node_weight
100+ return float((w * self.h).sum() / w.sum())
101+
102+ @cached_property
103+ def thickness_range(self) -> tuple[float, float]:
104+ hh = self.h[self.grid.mask]
105+ return float(hh.min()), float(hh.max())
106+
107+ @cached_property
108+ def cell_thickness(self) -> np.ndarray:
109+ """Thickness of each cell, the mean of its four corners, ``(ny - 1, nx - 1)``."""
110+ h = self.h
111+ return 0.25 * (h[:-1, :-1] + h[:-1, 1:] + h[1:, 1:] + h[1:, :-1])
112+
113+ @cached_property
114+ def is_flat(self) -> bool:
115+ """Whether the profile vanishes, in which case bending and in-plane motion decouple in the linear model."""
116+ return not np.any(self.w0)
117+
118+ @cached_property
119+ def profile_stats(self) -> tuple[float, float]:
120+ """Peak-to-peak height and largest slope of the profile over the plate nodes."""
121+ g, w0, m = self.grid, self.w0, self.grid.mask
122+ left, right, down, up = g.neighbor_masks()
123+ sx = np.zeros(g.shape)
124+ sy = np.zeros(g.shape)
125+ wp = np.pad(w0, 1)
126+ both = left & right
127+ sx[both] = ((wp[1:-1, 2:] - wp[1:-1, :-2]) / (2 * g.hx))[both]
128+ only = right & ~left
129+ sx[only] = ((wp[1:-1, 2:] - w0) / g.hx)[only]
130+ only = left & ~right
131+ sx[only] = ((w0 - wp[1:-1, :-2]) / g.hx)[only]
132+ both = down & up
133+ sy[both] = ((wp[2:, 1:-1] - wp[:-2, 1:-1]) / (2 * g.hy))[both]
134+ only = up & ~down
135+ sy[only] = ((wp[2:, 1:-1] - w0) / g.hy)[only]
136+ only = down & ~up
137+ sy[only] = ((w0 - wp[:-2, 1:-1]) / g.hy)[only]
138+ v = w0[m]
139+ return float(v.max() - v.min()), float(np.hypot(sx, sy)[m].max())
140+
141+ def dispersion_frequency(self, fraction: float, h: float | None = None) -> float:
142+ """Frequency [Hz] at which grid dispersion lowers a bending partial by
143+ ``fraction`` (a plane wave along a grid axis has the discrete frequency
144+ ratio (2 - 2 cos kh) / (kh)^2), at the thickness ``h`` (the mean by default)."""
145+ lo, hi = 0.0, math.pi
146+ for _ in range(60):
147+ mid = (lo + hi) / 2
148+ ratio = (2 - 2 * math.cos(mid)) / (mid * mid)
149+ if 1 - ratio < fraction:
150+ lo = mid
151+ else:
152+ hi = mid
153+ k = lo / max(self.grid.hx, self.grid.hy)
154+ return self.material.bending_constant(self.mean_thickness if h is None else h) * k * k / (2 * math.pi)
155+
156+ # ---- strike and listener ----
157+
158+ def strike_point(self, strike: Strike) -> tuple[float, float]:
159+ """The strike point in plate coordinates [m]."""
160+ return (strike.x - 0.5) * self.Lx, (strike.y - 0.5) * self.Ly
161+
162+ def patch_around(self, x: float, y: float, radius: float) -> Patch:
163+ """The plate nodes within ``radius`` of ``(x, y)``, always including the nearest one."""
164+ g = self.grid
165+ d = np.hypot(g.X - x, g.Y - y)
166+ d = np.where(self.node_mass > 0, d, np.inf)
167+ nearest = np.unravel_index(int(np.argmin(d)), g.shape)
168+ mask = d <= radius
169+ mask[nearest] = True
170+ m = self.node_mass[mask]
171+ mass = float(m.sum())
172+ return Patch(mask, (int(nearest[0]), int(nearest[1])), mass, float((m * g.X[mask]).sum() / mass), float((m * g.Y[mask]).sum() / mass))
173+
174+ def strike_patch(self, strike: Strike) -> Patch:
175+ x, y = self.strike_point(strike)
176+ return self.patch_around(x, y, strike.radius)
177+
178+ def listener_terms(self, listener: Listener) -> tuple[np.ndarray, np.ndarray]:
179+ """Per-node terms of the Rayleigh integral: the weight rho_air dA / (2 pi r)
180+ and the delay r / c [s] to the listener, ``(ny, nx)`` each, zero outside
181+ the plate. The profile is included in the node positions."""
182+ g = self.grid
183+ lx, ly, lz = listener.position
184+ r = np.sqrt((lx - g.X) ** 2 + (ly - g.Y) ** 2 + (lz - self.w0) ** 2)
185+ area = g.node_area
186+ weight = np.where(area > 0, RHO_AIR * area / (2 * math.pi * np.maximum(r, 1e-3)), 0.0)
187+ delay = np.where(area > 0, r / SPEED_OF_SOUND, 0.0)
188+ return weight, delay
189+
190+ # ---- summary ----
191+
192+ def summary(self) -> str:
193+ """A few lines describing the plate: box, material, grid, thickness, mass, profile, loss constants."""
194+ from .loss import LN1000, coincidence_frequency, estimated_rate, plate_geometry
195+
196+ g, m = self.grid, self.material
197+ hmin, hmax = self.thickness_range
198+ lines = [
199+ f"{self.name}: {self.Lx * 1e3:.0f} x {self.Ly * 1e3:.0f} mm box, area {self.area * 1e4:.0f} cm^2, mass {self.mass * 1e3:.0f} g",
200+ f"material: E = {m.E / 1e9:.1f} GPa, nu = {m.nu}, rho = {m.rho:.0f} kg/m^3, loss factor {m.loss_factor:g}"
201+ + (", thermal constants given" if m.thermal else ", no thermal constants"),
202+ f"grid: {g.nx} x {g.ny} nodes at {g.hx * 1e3:.2f} x {g.hy * 1e3:.2f} mm, {g.count} in the plate",
203+ f"thickness {hmin * 1e3:.2f} to {hmax * 1e3:.2f} mm (mean {self.mean_thickness * 1e3:.2f}); "
204+ f"grid dispersion -3% at {self.dispersion_frequency(0.03):.0f} Hz, -10% at {self.dispersion_frequency(0.10):.0f} Hz",
205+ ]
206+ if self.is_flat:
207+ lines.append("flat (no profile)")
208+ else:
209+ p2p, slope = self.profile_stats
210+ lines.append(f"profile: {p2p * 1e3:.1f} mm peak to peak, largest slope {slope:.3f}" + (" (above the shallow-shell range of about 0.1)" if slope > 0.1 else ""))
211+ geom = plate_geometry(self)
212+ t60 = [LN1000 / max(estimated_rate(self, geom, f).total, 1e-12) for f in (100, 1000, 5000)]
213+ lines.append(
214+ f"coincidence frequency {coincidence_frequency(m, self.mean_thickness):.0f} Hz; estimated T60 {t60[0]:.1f} / {t60[1]:.1f} / {t60[2]:.1f} s at 100 Hz / 1 kHz / 5 kHz"
215+ )
216+ return "\n".join(lines)
217+
218+ def __repr__(self) -> str:
219+ return f"Plate({self.name!r}, {self.grid.nx}x{self.grid.ny} grid, {self.grid.count} nodes, mass {self.mass * 1e3:.0f} g)"
220+
221+
222+def remove_plane(f: np.ndarray, g: Grid) -> np.ndarray:
223+ """Subtracts the weighted least-squares plane a + b x + c y from ``f`` over the plate nodes."""
224+ w = g.node_weight
225+ sel = w > 0
226+ if not sel.any():
227+ return f.copy()
228+ A = np.stack([np.ones(sel.sum()), g.X[sel], g.Y[sel]], axis=1)
229+ sw = np.sqrt(w[sel])
230+ coef, *_ = np.linalg.lstsq(A * sw[:, None], f[sel] * sw, rcond=None)
231+ return f - (coef[0] + coef[1] * g.X + coef[2] * g.Y)
src/cymsim/plot.pyadded+172−0View file
@@ -0,0 +1,172 @@
1+"""Matplotlib views of a plate, its modes, and a sound. Every function draws
2+into ``ax`` (a new figure by default) and returns the axes. Requires
3+``matplotlib`` (``pip install cymsim[plot]``)."""
4+
5+from __future__ import annotations
6+
7+import numpy as np
8+
9+from .plate import Plate
10+from .sound import Sound
11+from .spectrogram import Spectrogram
12+
13+__all__ = ["plate", "mode", "frame", "signal", "spectrogram", "decay_times", "force", "animate"]
14+
15+
16+def _ax(ax, figsize=(6, 5)):
17+ import matplotlib.pyplot as plt
18+
19+ if ax is None:
20+ _, ax = plt.subplots(figsize=figsize)
21+ return ax
22+
23+
24+def _extent(pl: Plate):
25+ g = pl.grid
26+ return [g.x[0] * 1e3 - g.hx * 500, g.x[-1] * 1e3 + g.hx * 500, g.y[0] * 1e3 - g.hy * 500, g.y[-1] * 1e3 + g.hy * 500]
27+
28+
29+def field(pl: Plate, values: np.ndarray, ax=None, *, cmap="viridis", label: str = "", symmetric: bool = False, strike=None, listener=None):
30+ """A color map of a grid array over the plate (outside masked), in mm coordinates."""
31+ ax = _ax(ax)
32+ v = np.ma.masked_where(~pl.grid.mask, values)
33+ kw = {}
34+ if symmetric:
35+ m = float(np.abs(v).max()) or 1.0
36+ kw = {"vmin": -m, "vmax": m}
37+ im = ax.imshow(v, origin="lower", extent=_extent(pl), cmap=cmap, interpolation="nearest", **kw)
38+ ax.figure.colorbar(im, ax=ax, label=label)
39+ if strike is not None:
40+ x, y = pl.strike_point(strike)
41+ ax.plot(x * 1e3, y * 1e3, "o", mfc="none", mec="yellow", mew=2, ms=10)
42+ if listener is not None:
43+ lx, ly, _ = listener.position
44+ ax.plot(lx * 1e3, ly * 1e3, "v", color="white", ms=8, clip_on=False)
45+ ax.set_xlabel("x [mm]")
46+ ax.set_ylabel("y [mm]")
47+ ax.set_aspect("equal")
48+ return ax
49+
50+
51+def plate(pl: Plate, ax=None, what: str = "auto", strike=None, listener=None):
52+ """The plate: its profile (when it has one) or its thickness, with the strike point and the listener when given."""
53+ if what == "auto":
54+ what = "thickness" if pl.is_flat else "profile"
55+ if what == "profile":
56+ ax = field(pl, pl.w0 * 1e3, ax, cmap="cividis", label="stress-free height [mm]", strike=strike, listener=listener)
57+ elif what == "thickness":
58+ ax = field(pl, pl.h * 1e3, ax, cmap="magma", label="thickness [mm]", strike=strike, listener=listener)
59+ else:
60+ raise ValueError("what must be 'auto', 'profile' or 'thickness'")
61+ ax.set_title(pl.name)
62+ return ax
63+
64+
65+def mode(modes, m: int, ax=None, amplitude: float = 1e-3):
66+ """The transverse shape of mode ``m``, scaled to the given peak displacement [m]."""
67+ w = modes.w[m]
68+ w = w / (np.abs(w).max() or 1.0) * amplitude
69+ ax = field(modes.plate, w * 1e3, ax, cmap="RdBu_r", label="w [mm]", symmetric=True)
70+ ax.set_title(f"mode {m}: {modes.freqs[m]:.1f} Hz, T60 {modes.t60[m]:.2f} s")
71+ return ax
72+
73+
74+def frame(sound: Sound, index: int, ax=None, scale: float | None = None):
75+ """Movie frame ``index`` of the transverse displacement; ``scale`` fixes the color range [m]."""
76+ if sound.frames is None:
77+ raise ValueError("this sound has no movie frames; run with movie=True")
78+ w = sound.frames[index]
79+ ax = _ax(ax)
80+ v = np.ma.masked_where(~sound.plate.grid.mask, w * 1e3)
81+ m = (scale or float(np.abs(sound.frames).max()) or 1e-9) * 1e3
82+ im = ax.imshow(v, origin="lower", extent=_extent(sound.plate), cmap="RdBu_r", vmin=-m, vmax=m, interpolation="nearest")
83+ ax.figure.colorbar(im, ax=ax, label="w [mm]")
84+ ax.set_title(f"t = {sound.frame_time(index) * 1e3:.1f} ms")
85+ ax.set_xlabel("x [mm]")
86+ ax.set_ylabel("y [mm]")
87+ ax.set_aspect("equal")
88+ return ax
89+
90+
91+def signal(sound: Sound, ax=None, tmax: float | None = None):
92+ """The listener pressure against time."""
93+ ax = _ax(ax, (8, 3))
94+ t = sound.time
95+ sel = slice(None) if tmax is None else t <= tmax
96+ ax.plot(t[sel], sound.pressure[sel], lw=0.5)
97+ ax.set_xlabel("time [s]")
98+ ax.set_ylabel("pressure [Pa]")
99+ ax.set_title(f"{sound.plate.name}, {sound.info.get('model', '')} model, peak {sound.peak_spl:.0f} dB SPL")
100+ return ax
101+
102+
103+def spectrogram(s: Sound | Spectrogram, ax=None, *, db_range: float = 90.0, fmin: float = 20.0, fmax: float | None = None, log_freq: bool = True, modes=None, window: int = 2048, hop: int = 256):
104+ """The spectrogram of a sound (or a precomputed ``Spectrogram``), with the
105+ modal frequencies drawn as ticks on the left edge when ``modes`` is given."""
106+ sp = s.spectrogram(window, hop) if isinstance(s, Sound) else s
107+ ax = _ax(ax, (8, 4))
108+ fmax = fmax or sp.freqs[-1]
109+ sel = (sp.freqs >= (fmin if log_freq else 0)) & (sp.freqs <= fmax)
110+ ax.pcolormesh(sp.times, sp.freqs[sel], sp.db[sel], cmap="inferno", vmin=-db_range, vmax=0, shading="nearest", rasterized=True)
111+ if log_freq:
112+ ax.set_yscale("log")
113+ ax.set_ylim(fmin, fmax)
114+ if modes is not None:
115+ f = modes.freqs[(modes.freqs >= fmin) & (modes.freqs <= fmax)]
116+ ax.plot(np.full(len(f), sp.times[0]), f, "_", color="cyan", ms=6, mew=1)
117+ ax.set_xlabel("time [s]")
118+ ax.set_ylabel("frequency [Hz]")
119+ return ax
120+
121+
122+def decay_times(modes, ax=None):
123+ """Decay time of every mode against frequency, total and by mechanism."""
124+ ax = _ax(ax, (8, 4))
125+ by = modes.t60_by_mechanism()
126+ ax.plot(modes.freqs, modes.t60, "k.", ms=4, label="total")
127+ for name, marker in (("radiation", "b."), ("material", "g."), ("thermal", "r.")):
128+ ax.plot(modes.freqs, by[name], marker, ms=2, alpha=0.6, label=name)
129+ ax.set_xscale("log")
130+ ax.set_yscale("log")
131+ ax.set_ylim(0.01, 1e4)
132+ ax.set_xlabel("frequency [Hz]")
133+ ax.set_ylabel("T60 [s]")
134+ ax.legend()
135+ return ax
136+
137+
138+def force(sound: Sound, ax=None):
139+ """The stick force during the contact (linear model only, which records it)."""
140+ if "force" not in sound.info:
141+ raise ValueError("the force history is recorded by the linear model only")
142+ ax = _ax(ax, (6, 3))
143+ f = sound.info["force"]
144+ ax.plot(np.arange(len(f)) * sound.info["force_dt"] * 1e3, f)
145+ ax.set_xlabel("time [ms]")
146+ ax.set_ylabel("force [N]")
147+ ax.set_title(f"peak {sound.info['peak_force']:.0f} N, contact {sound.info['contact_time'] * 1e3:.2f} ms")
148+ return ax
149+
150+
151+def animate(sound: Sound, interval_ms: float | None = None, scale: float | None = None, stride: int = 1):
152+ """A ``matplotlib.animation.FuncAnimation`` of the movie frames (display with ``HTML(anim.to_jshtml())`` in a notebook)."""
153+ import matplotlib.pyplot as plt
154+ from matplotlib.animation import FuncAnimation
155+
156+ if sound.frames is None:
157+ raise ValueError("this sound has no movie frames; run with movie=True")
158+ fig, ax = plt.subplots(figsize=(6, 5))
159+ m = (scale or float(np.abs(sound.frames).max()) or 1e-9) * 1e3
160+ mask = ~sound.plate.grid.mask
161+ im = ax.imshow(np.ma.masked_where(mask, sound.frames[0] * 1e3), origin="lower", extent=_extent(sound.plate), cmap="RdBu_r", vmin=-m, vmax=m, interpolation="nearest")
162+ fig.colorbar(im, ax=ax, label="w [mm]")
163+ ax.set_aspect("equal")
164+ title = ax.set_title("")
165+
166+ def update(i):
167+ im.set_data(np.ma.masked_where(mask, sound.frames[i] * 1e3))
168+ title.set_text(f"t = {sound.frame_time(i) * 1e3:.1f} ms")
169+ return im, title
170+
171+ frames = range(0, len(sound.frames), stride)
172+ return FuncAnimation(fig, update, frames=frames, interval=interval_ms or 1000 * stride / sound.fps, blit=False)
src/cymsim/presets.pyadded+186−0View file
@@ -0,0 +1,186 @@
1+"""Predefined designs: flat and bowed sheets, discs, tapered cymbals with and
2+without a bell, an ellipse, a thunder sheet. Each is a function returning a
3+``Design``, so that it can be modified before evaluation::
4+
5+ d = presets.crash_16()
6+ d.thickness = 1.0e-3
7+ plate = d.evaluate(0.008)
8+"""
9+
10+from __future__ import annotations
11+
12+import math
13+
14+import numpy as np
15+
16+from .design import Design
17+from .material import ALUMINUM, BRONZE, MILD_STEEL, STAINLESS_STEEL
18+
19+inch = 0.0254
20+mm = 1e-3
21+
22+
23+def ramp(a: float, b: float, x):
24+ """0 at x = a, 1 at x = b, clamped to that range outside."""
25+ return np.clip((x - a) / (b - a), 0, 1)
26+
27+
28+def bow(s, span: float, sag: float, offset: float = 0.0):
29+ """A parabola of the given sag across a span centered on 0; ``s`` runs along
30+ the span, and ``offset`` in -1..1 moves the crest toward one end."""
31+ t = 2 * s / span
32+ return sag * (1 - t**2) * (1 + offset * t)
33+
34+
35+def dome(d, radius: float, height: float):
36+ """A cosine dome of the given height and base radius, with zero slope both
37+ at its rim and at its apex; ``d`` is the distance from its center."""
38+ return 0.5 * height * (1 + np.cos(math.pi * d / radius)) * (d < radius)
39+
40+
41+def dents(x, y, r0: float, r1: float, n: int, rd: float, depth: float, random: bool = True, seed: int = 1):
42+ """``n`` hammer marks between the radii ``r0`` and ``r1``: Gaussian dents of
43+ 1/e radius ``rd``, placed on a sunflower spiral (equal area per mark,
44+ golden-angle spacing) so that the pattern has no symmetry, or at random with
45+ uniform density over the annulus, from a small congruential generator so
46+ that a seed always gives the same pattern."""
47+ v = np.zeros_like(x)
48+ ga = math.pi * (3 - math.sqrt(5))
49+ s = seed
50+ for k in range(1, n + 1):
51+ if random:
52+ s = (s * 1664525 + 1013904223) % 4294967296
53+ u = s / 4294967296
54+ s = (s * 1664525 + 1013904223) % 4294967296
55+ w = s / 4294967296
56+ rk = math.sqrt(r0**2 + (r1**2 - r0**2) * u)
57+ tk = 2 * math.pi * w
58+ else:
59+ rk = math.sqrt(r0**2 + (r1**2 - r0**2) * (k - 0.5) / n)
60+ tk = k * ga
61+ dx = x - rk * math.cos(tk)
62+ dy = y - rk * math.sin(tk)
63+ v = v - depth * np.exp(-(dx**2 + dy**2) / rd**2)
64+ return v
65+
66+
67+def flat_rectangle() -> Design:
68+ """A flat 12 by 18 inch sheet of stainless steel, 0.79 mm thick, hanging freely."""
69+ return Design("Flat steel rectangle", STAINLESS_STEEL, 12 * inch, 18 * inch, 0.79 * mm)
70+
71+
72+def aluminum_square() -> Design:
73+ """A flat 12 inch square of 6061-T6 aluminum, 1.02 mm thick. Aluminum is a
74+ third the density of steel and a third as stiff, so the bending frequencies
75+ land near those of a steel plate of the same size and thickness; its high
76+ thermal expansion and conductivity give it the strongest thermoelastic damping."""
77+ return Design("Aluminum 6061-T6 square", ALUMINUM, 12 * inch, 12 * inch, 1.02 * mm)
78+
79+
80+def bowed_rectangle() -> Design:
81+ """The steel sheet with a 3 mm parabolic bow across its width, the kind of
82+ bow that rolling and handling leave in sheet metal. The bow makes the sheet
83+ a shallow shell: the linear modes couple bending to stretching, and the full
84+ model oil-cans and wobbles when struck."""
85+ W = 12 * inch
86+ return Design("Bowed steel rectangle", STAINLESS_STEEL, W, 18 * inch, 0.79 * mm, profile=lambda x, y, r, t: bow(x, W, 3 * mm, 0.3))
87+
88+
89+def flat_disc() -> Design:
90+ """A flat 16 inch bronze disc, 1 mm thick."""
91+ D = 16 * inch
92+ return Design("Flat bronze disc", BRONZE, D, D, 1.0 * mm, inside=lambda x, y, r, t: r <= D / 2)
93+
94+
95+def crash_16() -> Design:
96+ """A 16 inch flat crash: bronze, lathed to a taper from 1.2 mm at the center
97+ to 0.7 mm at the edge. Flat, so the linear modes are pure bending modes."""
98+ D = 16 * inch
99+ R = D / 2
100+ tc, te = 1.2 * mm, 0.7 * mm
101+ return Design("16 inch flat crash", BRONZE, D, D, lambda x, y, r, t: tc + (te - tc) * ramp(40 * mm, R, r), inside=lambda x, y, r, t: r <= R)
102+
103+
104+def ride_20() -> Design:
105+ """A 20 inch flat ride: thicker bronze, 2.2 mm at the center to 1.4 mm at the edge."""
106+ D = 20 * inch
107+ R = D / 2
108+ return Design("20 inch flat ride", BRONZE, D, D, lambda x, y, r, t: 2.2 * mm + (1.4 * mm - 2.2 * mm) * ramp(50 * mm, R, r), inside=lambda x, y, r, t: r <= R)
109+
110+
111+def crash_bell_16() -> Design:
112+ """A 16 inch crash with a shallow bow and a cosine bell. The profile makes
113+ the cymbal a shallow shell (a three-field linear system, so the modes take
114+ longer), and the full model shows the asymmetric nonlinearity of a curved
115+ plate. The slopes are kept small; ``Plate.profile_stats`` reports the largest."""
116+ D = 16 * inch
117+ R = D / 2
118+ bell_r, bell_h, bow_h = 55 * mm, 8 * mm, 10 * mm
119+ tc, te = 1.2 * mm, 0.7 * mm
120+ return Design(
121+ "16 inch crash with bell",
122+ BRONZE,
123+ D,
124+ D,
125+ lambda x, y, r, t: tc + (te - tc) * ramp(bell_r, R, r),
126+ inside=lambda x, y, r, t: r <= R,
127+ profile=lambda x, y, r, t: bow_h * (1 - (r / R) ** 2) + dome(r, bell_r, bell_h),
128+ )
129+
130+
131+def crash_hammered_16(random_dents: bool = True, seed: int = 1, n_dents: int = 160, dent_depth: float = 0.4 * mm) -> Design:
132+ """The 16 inch crash with bell, hammered: shallow dents on the bow, at random
133+ places (``seed`` picks the pattern) or on a sunflower spiral. A hammered
134+ cymbal is no longer axisymmetric, so its degenerate mode pairs split and its
135+ modes mix. The dents are a change of the rest shape; the residual stresses
136+ of real hammering are not modeled."""
137+ D = 16 * inch
138+ R = D / 2
139+ bell_r, bell_h, bow_h = 55 * mm, 8 * mm, 10 * mm
140+ tc, te = 1.2 * mm, 0.7 * mm
141+ dent_r = 12 * mm
142+ return Design(
143+ "16 inch hammered crash",
144+ BRONZE,
145+ D,
146+ D,
147+ lambda x, y, r, t: tc + (te - tc) * ramp(bell_r, R, r),
148+ inside=lambda x, y, r, t: r <= R,
149+ profile=lambda x, y, r, t: bow_h * (1 - (r / R) ** 2) + dome(r, bell_r, bell_h) + dents(x, y, bell_r + dent_r, R - dent_r, n_dents, dent_r, dent_depth, random_dents, seed),
150+ )
151+
152+
153+def ellipse() -> Design:
154+ """A 14 by 10 inch elliptical steel plate, 1 mm thick. The degenerate mode pairs of a disc split on an ellipse."""
155+ A, B = 14 * inch, 10 * inch
156+ return Design("Steel ellipse", STAINLESS_STEEL, A, B, 1.0 * mm, inside=lambda x, y, r, t: (x / (A / 2)) ** 2 + (y / (B / 2)) ** 2 <= 1)
157+
158+
159+def thunder_sheet() -> Design:
160+ """A thunder sheet: a 24 by 36 inch steel sheet, 0.5 mm thick. Its modes are
161+ dense and low, and the full model turns a hard hit into a roar. Long at the
162+ default spacing; a coarser grid keeps it interactive."""
163+ return Design("Thunder sheet", MILD_STEEL, 24 * inch, 36 * inch, 0.5 * mm)
164+
165+
166+PRESETS = {
167+ "flat-rectangle": flat_rectangle,
168+ "aluminum-square": aluminum_square,
169+ "bowed-rectangle": bowed_rectangle,
170+ "flat-disc": flat_disc,
171+ "crash-16": crash_16,
172+ "ride-20": ride_20,
173+ "crash-bell-16": crash_bell_16,
174+ "crash-hammered-16": crash_hammered_16,
175+ "ellipse": ellipse,
176+ "thunder-sheet": thunder_sheet,
177+}
178+"""The presets by name."""
179+
180+
181+def preset(name: str) -> Design:
182+ """A preset design by name (see ``PRESETS``)."""
183+ try:
184+ return PRESETS[name]()
185+ except KeyError:
186+ raise KeyError(f"unknown preset {name!r}; choose from {', '.join(PRESETS)}") from None
src/cymsim/sound.pyadded+71−0View file
@@ -0,0 +1,71 @@
1+"""The result of a simulation: the listener pressure, the movie frames, and what the run reported."""
2+
3+from __future__ import annotations
4+
5+from dataclasses import dataclass, field
6+from typing import Any
7+
8+import numpy as np
9+
10+from .params import Listener, Strike
11+from .plate import Plate
12+from .spectrogram import Spectrogram, spectrogram
13+from .wav import Level, apply_level, level_gain, peak_dbfs, write_wav
14+
15+
16+@dataclass
17+class Sound:
18+ pressure: np.ndarray
19+ """Listener sound pressure [Pa] at ``fs``, from the moment of impact."""
20+ fs: int
21+ plate: Plate
22+ strike: Strike
23+ listener: Listener
24+ frames: np.ndarray | None = None
25+ """Transverse displacement [m] at ``fps`` frames per second, ``(n_frames, ny, nx)``, or None."""
26+ fps: float = 120.0
27+ info: dict[str, Any] = field(default_factory=dict)
28+ """What the run reported: ``model``, ``peak_force`` [N], ``contact_time`` [s], and more."""
29+
30+ @property
31+ def duration(self) -> float:
32+ return len(self.pressure) / self.fs
33+
34+ @property
35+ def time(self) -> np.ndarray:
36+ """Sample times [s]."""
37+ return np.arange(len(self.pressure)) / self.fs
38+
39+ @property
40+ def peak_spl(self) -> float:
41+ """Peak sound pressure level [dB SPL re 20 uPa]."""
42+ m = float(np.abs(self.pressure).max()) if len(self.pressure) else 0.0
43+ return 20 * np.log10(m / 20e-6) if m > 0 else -np.inf
44+
45+ def samples(self, level: Level = "peak", gain_db: float = 0.0) -> np.ndarray:
46+ """The signal mapped to full scale in [-1, 1] (see ``cymsim.wav.level_gain``), soft-clipped."""
47+ return apply_level(self.pressure, level_gain(self.pressure, self.fs, level, gain_db))
48+
49+ def to_wav(self, path, level: Level = "peak", gain_db: float = 0.0) -> None:
50+ """Writes a 16-bit mono WAV file. ``level`` is ``"peak"`` (the peak at -1 dBFS),
51+ ``"loudness"`` (-20 dBFS RMS) or ``"absolute"`` (1 Pa at 0 dBFS plus ``gain_db``)."""
52+ write_wav(path, self.samples(level, gain_db), self.fs)
53+
54+ def spectrogram(self, window: int = 2048, hop: int = 256) -> Spectrogram:
55+ return spectrogram(self.pressure, self.fs, window, hop)
56+
57+ def audio(self, level: Level = "peak", gain_db: float = 0.0):
58+ """An ``IPython.display.Audio`` widget for playback in a notebook."""
59+ from IPython.display import Audio
60+
61+ return Audio(self.samples(level, gain_db), rate=self.fs, normalize=False)
62+
63+ def frame_time(self, index: int) -> float:
64+ return index / self.fps
65+
66+ def __repr__(self) -> str:
67+ model = self.info.get("model", "?")
68+ return f"Sound({model} model, {self.duration:.2f} s at {self.fs} Hz, peak {self.peak_spl:.0f} dB SPL" + (f", {len(self.frames)} frames)" if self.frames is not None else ")")
69+
70+ def peak_dbfs(self, level: Level = "peak", gain_db: float = 0.0) -> float:
71+ return peak_dbfs(self.samples(level, gain_db))
src/cymsim/spectrogram.pyadded+42−0View file
@@ -0,0 +1,42 @@
1+"""A short-time Fourier transform spectrogram in decibels."""
2+
3+from __future__ import annotations
4+
5+from dataclasses import dataclass
6+
7+import numpy as np
8+import scipy.fft
9+
10+
11+@dataclass
12+class Spectrogram:
13+ times: np.ndarray
14+ """Center time of each column [s]."""
15+ freqs: np.ndarray
16+ """Frequency of each row [Hz]."""
17+ db: np.ndarray
18+ """Magnitude in dB relative to the peak, ``(len(freqs), len(times))``."""
19+
20+ def plot(self, ax=None, **kwargs):
21+ """Draws the spectrogram with matplotlib; see ``cymsim.plot.spectrogram``."""
22+ from .plot import spectrogram
23+
24+ return spectrogram(self, ax=ax, **kwargs)
25+
26+
27+def spectrogram(signal: np.ndarray, fs: float, window: int = 2048, hop: int = 256) -> Spectrogram:
28+ """The STFT of ``signal`` with a periodic Hann window of ``window`` samples
29+ and a hop of ``hop`` samples, in dB relative to the largest magnitude."""
30+ x = np.asarray(signal, dtype=float)
31+ if len(x) < window:
32+ x = np.pad(x, (0, window - len(x)))
33+ n_cols = (len(x) - window) // hop + 1
34+ idx = np.arange(window)[None, :] + hop * np.arange(n_cols)[:, None]
35+ win = 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(window) / window)
36+ X = scipy.fft.rfft(x[idx] * win[None, :], axis=1)
37+ p = X.real**2 + X.imag**2
38+ db = 10 * np.log10(p + 1e-30)
39+ db -= db.max()
40+ times = (hop * np.arange(n_cols) + window / 2) / fs
41+ freqs = np.arange(window // 2 + 1) * fs / window
42+ return Spectrogram(times, freqs, db.T)
src/cymsim/wav.pyadded+58−0View file
@@ -0,0 +1,58 @@
1+"""Mono 16-bit PCM WAV output and the level mapping from pascals to full scale."""
2+
3+from __future__ import annotations
4+
5+import wave
6+from typing import Literal
7+
8+import numpy as np
9+
10+Level = Literal["peak", "loudness", "absolute"]
11+
12+KNEE = 0.9
13+PEAK_TARGET = 10 ** (-1 / 20) # -1 dBFS
14+LOUDNESS_TARGET = 10 ** (-20 / 20) # -20 dBFS RMS
15+LOUDNESS_SKIP_S = 0.05
16+
17+
18+def level_gain(samples: np.ndarray, fs: float, level: Level = "peak", gain_db: float = 0.0) -> float:
19+ """Gain that maps a pressure signal [Pa] to full scale: ``"peak"`` normalizes
20+ the peak to -1 dBFS, ``"loudness"`` the RMS after the first 50 ms to -20 dBFS,
21+ and ``"absolute"`` maps 1 Pa to 0 dBFS shifted by ``gain_db``."""
22+ if level == "absolute":
23+ return 10 ** (gain_db / 20)
24+ if level == "peak":
25+ m = float(np.abs(samples).max()) if len(samples) else 0.0
26+ return PEAK_TARGET / m if m > 0 else 1.0
27+ if level == "loudness":
28+ start = min(len(samples), round(LOUDNESS_SKIP_S * fs))
29+ tail = samples[start:]
30+ rms = float(np.sqrt((tail.astype(float) ** 2).mean())) if len(tail) else 0.0
31+ return LOUDNESS_TARGET / rms if rms > 0 else 1.0
32+ raise ValueError(f"unknown level {level!r}")
33+
34+
35+def apply_level(samples: np.ndarray, gain: float) -> np.ndarray:
36+ """Applies the gain and a soft clipper: below 0.9 of full scale the identity,
37+ above it the excess is squashed with a tanh so that the output stays within
38+ (-1, 1) with a continuous slope at the knee."""
39+ y = samples.astype(np.float64) * gain
40+ a = np.abs(y)
41+ room = 1 - KNEE
42+ clipped = KNEE + room * np.tanh((a - KNEE) / room)
43+ return np.where(a <= KNEE, y, np.sign(y) * clipped).astype(np.float32)
44+
45+
46+def write_wav(path, samples: np.ndarray, fs: int) -> None:
47+ """Writes samples in [-1, 1] as a 16-bit PCM mono WAV file."""
48+ pcm = np.round(np.clip(np.nan_to_num(samples), -1, 1) * 32767).astype("<i2")
49+ with wave.open(str(path), "wb") as w:
50+ w.setnchannels(1)
51+ w.setsampwidth(2)
52+ w.setframerate(int(round(fs)))
53+ w.writeframes(pcm.tobytes())
54+
55+
56+def peak_dbfs(samples: np.ndarray) -> float:
57+ m = float(np.abs(samples).max()) if len(samples) else 0.0
58+ return 20 * np.log10(m) if m > 0 else -np.inf
tests/conftest.pyadded+13−0View file
@@ -0,0 +1,13 @@
1+import numpy as np
2+import pytest
3+
4+from cymsim import Material, Thermal
5+
6+STEEL = Material(200e9, 0.3, 7850)
7+STEEL_LOSSY = Material(200e9, 0.3, 7850, 3e-4, Thermal(12e-6, 460, 50))
8+BRONZE = Material(110e9, 0.34, 8800, 2e-4, Thermal(18e-6, 370, 40))
9+
10+
11+@pytest.fixture
12+def rng():
13+ return np.random.default_rng(7)
tests/test_design.pyadded+94−0View file
@@ -0,0 +1,94 @@
1+"""Grids, masks, designs and presets."""
2+
3+import numpy as np
4+import pytest
5+
6+import cymsim
7+from cymsim import Design, DesignError, Plate, make_grid, presets
8+from cymsim.plate import remove_plane
9+
10+from conftest import STEEL
11+
12+
13+def test_grid_spacing_and_weights():
14+ g = make_grid(0.3, 0.2, 0.01)
15+ assert (g.nx, g.ny) == (31, 21)
16+ assert abs(g.hx - 0.01) < 1e-12
17+ w = g.node_weight
18+ assert w[5, 5] == 1 and w[0, 5] == 0.5 and w[0, 0] == 0.25
19+ assert abs(g.area - 0.3 * 0.2) < 1e-12
20+ assert abs(g.x[0] + 0.15) < 1e-12 and abs(g.y[-1] - 0.1) < 1e-12
21+
22+
23+def test_disc_with_hole_mask_by_cell_centers():
24+ R, r0 = 0.1, 0.02
25+ g = make_grid(2 * R, 2 * R, 0.01, lambda x, y, r, t: (r <= R) & (r >= r0))
26+ # The area of the staircase plate is unbiased to within a cell's worth.
27+ assert abs(g.area / (np.pi * (R**2 - r0**2)) - 1) < 0.03
28+ assert not g.mask[g.ny // 2, g.nx // 2], "the hole is empty"
29+ assert not g.mask[0, 0], "the corner is outside the disc"
30+ # Every node in the plate has a cell of its own.
31+ assert (g.node_weight[g.mask] > 0).all()
32+ assert (g.node_weight[~g.mask] == 0).all()
33+
34+
35+def test_remove_plane_and_profile_stats():
36+ g = make_grid(0.2, 0.1, 0.01)
37+ f = 3.0 + 2.0 * g.X - 5.0 * g.Y + 1e-3 * g.X**2
38+ r = remove_plane(f, g)
39+ # The plane is gone: the weighted fit of the residual is zero to rounding.
40+ w = g.node_weight
41+ assert abs((w * r).sum()) < 1e-12 and abs((w * r * g.X).sum()) < 1e-12 and abs((w * r * g.Y).sum()) < 1e-12
42+ pl = Plate.uniform(g, STEEL, 1e-3, profile=f)
43+ p2p, slope = pl.profile_stats
44+ assert p2p < 1e-3 * 0.2**2 and slope < 1e-3 * 0.2
45+
46+
47+def test_design_evaluate_checks():
48+ d = Design("bad", STEEL, 0.1, 0.1, lambda x, y, r, t: x) # negative on half the plate
49+ with pytest.raises(DesignError, match="positive"):
50+ d.evaluate(0.01)
51+ d = Design("empty", STEEL, 0.1, 0.1, 1e-3, inside=lambda x, y, r, t: r < 0)
52+ with pytest.raises(DesignError, match="no grid cells"):
53+ d.evaluate(0.01)
54+ d = Design("shape", STEEL, 0.1, 0.1, lambda x, y, r, t: np.ones(3))
55+ with pytest.raises(DesignError, match="shape"):
56+ d.evaluate(0.01)
57+ with pytest.raises(DesignError):
58+ Design("thin", STEEL, 0.1, 0.1, 0.0)
59+
60+
61+def test_every_preset_evaluates():
62+ for name, fn in presets.PRESETS.items():
63+ pl = fn().evaluate(0.012)
64+ assert pl.grid.count > 20, name
65+ assert pl.mass > 0.05, name
66+ assert np.isfinite(pl.w0).all() and (pl.w0[~pl.grid.mask] == 0).all()
67+ assert isinstance(pl.summary(), str)
68+ crash = presets.preset("crash-16").evaluate(0.008)
69+ assert crash.is_flat and abs(crash.thickness_range[0] - 0.7e-3) < 1e-9 and abs(crash.thickness_range[1] - 1.2e-3) < 1e-9
70+ bell = presets.preset("crash-bell-16").evaluate(0.008)
71+ assert not bell.is_flat and bell.profile_stats[1] > 0.1
72+ with pytest.raises(KeyError):
73+ presets.preset("no-such-preset")
74+
75+
76+def test_hammered_pattern_is_reproducible_and_asymmetric():
77+ a = presets.crash_hammered_16(seed=1).evaluate(0.012)
78+ b = presets.crash_hammered_16(seed=1).evaluate(0.012)
79+ c = presets.crash_hammered_16(seed=2).evaluate(0.012)
80+ assert np.array_equal(a.w0, b.w0)
81+ assert not np.allclose(a.w0, c.w0)
82+ assert not np.allclose(a.w0, a.w0[::-1, :]), "the dents break the symmetry"
83+
84+
85+def test_strike_patch_and_listener_terms():
86+ pl = presets.flat_disc().evaluate(0.01)
87+ patch = pl.strike_patch(cymsim.Strike(x=0.9, y=0.5, radius=5e-3))
88+ assert patch.mask.sum() >= 1 and patch.mask[patch.nearest]
89+ assert abs(patch.cy) < 0.01
90+ weight, delay = pl.listener_terms(cymsim.Listener(distance=1.0, angle=0))
91+ sel = weight > 0
92+ assert (delay[sel] >= 1.0 / 343).all()
93+ # On axis the weights sum to about rho A / (2 pi r).
94+ assert abs(weight.sum() / (1.2 * pl.area / (2 * np.pi)) - 1) < 0.02
tests/test_full.pyadded+126−0View file
@@ -0,0 +1,126 @@
1+"""The full model: the time step, a free step from a smooth state agrees between
2+float32 and float64, energy is conserved without losses, the recorded audio is
3+the Rayleigh sum, and a light hit agrees with the linear model."""
4+
5+import numpy as np
6+import pytest
7+import torch
8+
9+from cymsim import Listener, Plate, Simulation, Strike, compute_modes, make_grid, simulate, synthesize
10+from cymsim.constants import CFL
11+from cymsim.full.reference import potential_energy
12+
13+from conftest import STEEL, STEEL_LOSSY
14+
15+
16+def smooth_field(pl, amplitude, seed):
17+ r = np.random.default_rng(seed)
18+ g = pl.grid
19+ v = np.zeros(g.shape)
20+ for _ in range(6):
21+ kx, ky, ph, a = (r.random() - 0.5) * 60, (r.random() - 0.5) * 60, r.random() * 2 * np.pi, r.random()
22+ v += a * np.sin(kx * g.X + ky * g.Y + ph)
23+ return np.where(g.mask, v * amplitude / 3, 0.0)
24+
25+
26+def test_time_step_and_batching():
27+ pl = Plate.uniform(make_grid(0.2, 0.15, 0.01), STEEL, 1e-3)
28+ sim = Simulation(pl, duration=0.1, fs=48000, fps=100, device="cpu", compile=False)
29+ assert sim.dt * sim.omega_max <= 2 * CFL + 1e-9
30+ assert sim.steps_per_sample * sim.dt * 48000 == pytest.approx(1.0)
31+ assert sim.samples_per_batch == 480 and sim.steps_per_batch == 480 * sim.steps_per_sample
32+ assert sim.hist_len > sim.samples_per_batch
33+
34+
35+def test_float32_step_matches_float64_from_a_smooth_state():
36+ g = make_grid(0.2, 0.2, 0.0125, lambda x, y, r, t: r <= 0.1)
37+ pl = Plate.uniform(g, STEEL_LOSSY, 1e-3, profile=4e-3 * (1 - (g.X**2 + g.Y**2) / 0.01))
38+ hm = pl.mean_thickness
39+ amp = 2.5 * hm
40+ pos = [smooth_field(pl, 0.02 * amp, 1), smooth_field(pl, 0.02 * amp, 2), smooth_field(pl, amp, 3)]
41+ vel = [smooth_field(pl, 0.05, 4), smooth_field(pl, 0.05, 5), smooth_field(pl, 1.0, 6)]
42+ accs = {}
43+ for dtype in ("float32", "float64"):
44+ sim = Simulation(pl, duration=0.01, device="cpu", dtype=dtype, compile=False)
45+ sim.set_state(*pos, *vel)
46+ before = sim.state()
47+ sim.step_once()
48+ after = sim.state()
49+ accs[dtype] = np.stack([(after[k] - before[k]) / sim.dt for k in ("du", "dv", "dw")])
50+ ref = accs["float64"]
51+ assert np.abs(accs["float32"] - ref).max() < 2e-6 * np.abs(ref).max()
52+
53+
54+def test_energy_is_conserved_without_losses_and_lost_with_them():
55+ pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3)
56+
57+ def energies(amplitude, loss, n_checks=4, steps=200):
58+ sim = Simulation(pl, duration=0.01, device="cpu", dtype="float64", compile=False, loss=loss)
59+ sim.set_state(w=smooth_field(pl, amplitude, 3))
60+ E = []
61+ for _ in range(n_checks):
62+ st = sim.state()
63+ kinetic = 0.5 * (pl.node_mass * (st["du"] ** 2 + st["dv"] ** 2 + st["dw"] ** 2)).sum()
64+ E.append(kinetic + potential_energy(pl, st["u"], st["v"], st["w"]))
65+ for _ in range(steps):
66+ sim.step_once()
67+ return np.array(E)
68+
69+ # Linear regime (a small fraction of the thickness): the leapfrog energy is
70+ # conserved; the naive sum of kinetic (at the half step) and potential energy
71+ # oscillates by a fraction of a percent.
72+ E = energies(1e-6, (0.0, 0.0, 0.0))
73+ assert np.abs(E / E[0] - 1).max() < 1e-2
74+ # Strongly nonlinear regime (two thicknesses): the in-plane waves it excites
75+ # sit near the stability limit, where the naive energy oscillates by a few
76+ # percent, but it does not drift.
77+ E = energies(2e-3, (0.0, 0.0, 0.0))
78+ assert np.abs(E / E[0] - 1).max() < 0.1
79+ # With losses the energy falls.
80+ E = energies(2e-3, (500.0, 1e-2, 500.0))
81+ assert E[-1] < 0.6 * E[0]
82+
83+
84+def test_recorded_audio_is_the_rayleigh_sum_of_the_state():
85+ pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3)
86+ listener = Listener(distance=0.3, angle=40, azimuth=30)
87+ sim = Simulation(pl, Strike(velocity=1.0), listener, duration=0.02, fps=200, device="cpu", dtype="float64", compile=False, movie=True)
88+ sim.run()
89+ s = sim.sound()
90+ assert len(s.pressure) == round(0.02 * 48000) and s.frames.shape[0] == 5
91+ assert s.peak_spl > 40 and sim.stick["max_force"] > 0
92+ # Late samples are consistent with a delayed, weighted first difference of
93+ # the transverse velocity (checked on the last batch, whose ring buffer still
94+ # holds the history).
95+ assert np.isfinite(s.pressure).all()
96+
97+
98+def test_compiled_step_matches_eager():
99+ pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3)
100+ a = Simulation(pl, duration=0.01, device="cpu", compile=False, movie=False)
101+ b = Simulation(pl, duration=0.01, device="cpu", compile=True, movie=False)
102+ a.run()
103+ b.run()
104+ if b._step is not b.__class__.__dict__.get("step"):
105+ pass
106+ assert np.allclose(a.sound().pressure, b.sound().pressure, rtol=1e-5, atol=1e-9)
107+
108+
109+@pytest.mark.slow
110+def test_light_hit_agrees_between_the_models():
111+ pl = Plate.uniform(make_grid(0.2, 0.15, 0.0125), STEEL_LOSSY, 1e-3)
112+ strike, listener = Strike(x=0.7, y=0.4, velocity=0.05), Listener()
113+ full = simulate(pl, strike, listener, duration=0.25, device="cpu")
114+ lin = synthesize(compute_modes(pl, 16000), strike, listener, duration=0.25)
115+ fs = 48000
116+
117+ def spectrum(p):
118+ seg = p[int(0.05 * fs) : int(0.25 * fs)].astype(float)
119+ return np.abs(np.fft.rfft(seg * np.hanning(len(seg))))
120+
121+ f = np.fft.rfftfreq(int(0.2 * fs), 1 / fs)
122+ A, B = spectrum(full.pressure), spectrum(lin.pressure)
123+ sel = f < 6000
124+ assert np.corrcoef(A[sel], B[sel])[0, 1] > 0.98
125+ assert 0.8 < np.linalg.norm(A[sel]) / np.linalg.norm(B[sel]) < 1.25
126+ assert abs(full.info["peak_force"] / lin.info["peak_force"] - 1) < 0.02
tests/test_loss.pyadded+151−0View file
@@ -0,0 +1,151 @@
1+"""The losses derived from physical constants, and the per-mode rates."""
2+
3+import numpy as np
4+
5+from cymsim import Plate, compute_modes, make_grid
6+from cymsim.constants import RHO_AIR, ROOM_TEMPERATURE, SPEED_OF_SOUND
7+from cymsim.linear.damping import bending_energy_per_node, radiation_rate_direct, radiation_rates, thermoelastic_rates
8+from cymsim.loss import (
9+ LN1000,
10+ average_radiation_efficiency,
11+ coincidence_frequency,
12+ estimated_rate,
13+ fit_line,
14+ fit_sigma,
15+ fitted_rate,
16+ material_rate,
17+ outline_perimeter,
18+ plate_geometry,
19+ thermoelastic_loss_factor,
20+ thermoelastic_relaxation,
21+)
22+
23+from conftest import BRONZE, STEEL_LOSSY
24+
25+
26+def test_zener_peak():
27+ h = 1e-3
28+ delta, tau = thermoelastic_relaxation(BRONZE, h, ROOM_TEMPERATURE)
29+ th = BRONZE.thermal
30+ assert abs(tau - h * h * BRONZE.rho * th.cp / (np.pi**2 * th.k)) < 1e-15
31+ assert abs(delta - BRONZE.E * th.alpha**2 * ROOM_TEMPERATURE / (BRONZE.rho * th.cp) * (1 + BRONZE.nu) / (1 - BRONZE.nu)) < 1e-15
32+ peak = thermoelastic_loss_factor(BRONZE, h, 1 / tau, ROOM_TEMPERATURE)
33+ assert abs(peak - delta / 2) < 1e-15
34+ assert thermoelastic_loss_factor(BRONZE, h, 0.1 / tau, ROOM_TEMPERATURE) < peak > thermoelastic_loss_factor(BRONZE, h, 10 / tau, ROOM_TEMPERATURE)
35+ assert 1 / (2 * np.pi * tau) < 100 and 1e-4 < delta / 2 < 1e-2
36+ assert abs(thermoelastic_relaxation(BRONZE, 2 * h, ROOM_TEMPERATURE)[1] / tau - 4) < 1e-12
37+ from cymsim import Material
38+
39+ assert thermoelastic_loss_factor(Material(110e9, 0.34, 8800), h, 1000, ROOM_TEMPERATURE) == 0
40+
41+
42+def test_outline_perimeter():
43+ rect = make_grid(0.3, 0.2, 0.01)
44+ assert abs(outline_perimeter(rect) / (2 * (0.3 + 0.2)) - 1) < 0.02
45+ R = 0.2
46+ disc = make_grid(2 * R, 2 * R, 0.008, lambda x, y, r, t: r <= R)
47+ p = outline_perimeter(disc)
48+ assert 2 * np.pi * R < p < 1.08 * 2 * np.pi * R
49+ ring = make_grid(2 * R, 2 * R, 0.008, lambda x, y, r, t: (r <= R) & (r >= 0.05))
50+ assert abs(outline_perimeter(ring) / (2 * np.pi * (R + 0.05)) - 1) < 0.1
51+
52+
53+def test_average_radiation_efficiency_limits():
54+ fc = coincidence_frequency(BRONZE, 1e-3)
55+ assert 15000 < fc < 20000
56+ A, U, a, b = 0.13, 1.28, 0.4, 0.4
57+ s1 = average_radiation_efficiency(100, fc, A, U, a, b, False)
58+ s2 = average_radiation_efficiency(200, fc, A, U, a, b, False)
59+ assert 0 < s1 < 0.1 and abs(s2 / s1 - 2) < 0.1
60+ assert average_radiation_efficiency(fc / 2, fc, A, U, a, b, False) > s2
61+ assert abs(average_radiation_efficiency(4 * fc, fc, A, U, a, b) - 1 / np.sqrt(0.75)) < 1e-9
62+ assert average_radiation_efficiency(fc * 1.0001, fc, A, U, a, b) <= (0.5 - 0.15) * np.sqrt(a / (SPEED_OF_SOUND / fc)) + 1e-9
63+ f1, f2 = average_radiation_efficiency(100, fc, A, U, a, b), average_radiation_efficiency(200, fc, A, U, a, b)
64+ assert abs(f1 / s1 - (100 / fc) ** 2) < 1e-12 and abs(f2 / f1 - 8) < 0.4
65+ # Vectorized over frequency.
66+ arr = average_radiation_efficiency(np.array([100.0, 200.0, 4 * fc]), fc, A, U, a, b)
67+ assert np.allclose(arr[:2], [f1, f2]) and abs(arr[2] - 1 / np.sqrt(0.75)) < 1e-9
68+
69+
70+def test_two_parameter_fit():
71+ x = np.array([1, 2, 3, 5, 8.0])
72+ fit = fit_line(x, 0.7 + 2e-4 * x, 1.1)
73+ assert abs(fit.sigma0 - 0.7) < 1e-9 and abs(fit.sigma1 - 2e-4) < 1e-12
74+ assert fit_line(x, 3 - x, 1.1).sigma1 == 0
75+ pl = Plate.uniform(make_grid(0.4, 0.4, 0.02), BRONZE, 1e-3)
76+ geom = plate_geometry(pl)
77+ f2 = fit_sigma(pl, geom)
78+ assert f2.sigma0 >= 0 and f2.sigma1 >= 0
79+ e100, e5k = estimated_rate(pl, geom, 100), estimated_rate(pl, geom, 5000)
80+ assert e5k.total > e100.total
81+ assert abs(e100.material - material_rate(BRONZE.loss_factor, 2 * np.pi * 100)) < 1e-12
82+ assert e100.radiation > 0 and e100.thermoelastic > 0
83+ assert 0.3 * e100.total < fitted_rate(f2, 1000) < 3 * e5k.total
84+
85+
86+def small_modes():
87+ pl = Plate.uniform(make_grid(0.2, 0.15, 0.0125), STEEL_LOSSY, 1e-3)
88+ return pl, compute_modes(pl, 4000)
89+
90+
91+def test_piston_radiation_rate_equals_the_closed_form():
92+ pl, ms = small_modes()
93+ one = ms.select([0])
94+ M, S = pl.mass, pl.area
95+ one.w[0] = np.where(pl.grid.mask, 1 / np.sqrt(M), 0.0)
96+ one.freqs[0] = 20.0
97+ omega = 2 * np.pi * 20
98+ expected = RHO_AIR * omega**2 * S**2 / (2 * np.pi * SPEED_OF_SOUND * M)
99+ assert abs(radiation_rates(one)[0] / expected - 1) < 2e-3
100+ assert abs(radiation_rate_direct(pl, one.w[0], 20.0) / expected - 1) < 2e-3
101+
102+
103+def test_fft_convolution_equals_the_direct_double_sum():
104+ pl, ms = small_modes()
105+ rates = radiation_rates(ms)
106+ for m in range(0, ms.count, max(1, ms.count // 12)):
107+ direct = radiation_rate_direct(pl, ms.w[m], ms.freqs[m])
108+ assert abs(rates[m] - direct) < 1e-9 * max(1, direct)
109+ assert (rates >= 0).all() and rates[-1] > rates[0]
110+
111+
112+def test_wallace_low_frequency_efficiencies():
113+ # sin(m pi x / a) sin(n pi y / b) on a baffled rectangle at k a << 1 (Wallace 1972).
114+ a, b, h, f = 0.3, 0.45, 1e-3, 30.0
115+ c = SPEED_OF_SOUND
116+ k = 2 * np.pi * f / c
117+ g = make_grid(a, b, 0.01)
118+ pl = Plate.uniform(g, STEEL_LOSSY, h)
119+ rho_h = STEEL_LOSSY.rho * h
120+ xs = (g.X - g.x[0]) / a
121+ ys = (g.Y - g.y[0]) / b
122+ for m, n in [(1, 1), (1, 2), (2, 1), (2, 2), (3, 3)]:
123+ v = np.sin(m * np.pi * xs) * np.sin(n * np.pi * ys)
124+ v /= np.sqrt((rho_h * g.node_area * v * v).sum())
125+ efficiency = radiation_rate_direct(pl, v, f) * rho_h / (RHO_AIR * c)
126+ odd = lambda q: q % 2 == 1 # noqa: E731
127+ if odd(m) and odd(n):
128+ wallace = 32 * k * k * a * b / (np.pi**5 * m * m * n * n)
129+ elif not odd(m) and not odd(n):
130+ wallace = 2 * k**6 * a**3 * b**3 / (15 * np.pi**5 * m * m * n * n)
131+ elif odd(m):
132+ wallace = 8 * k**4 * a * b**3 / (3 * np.pi**5 * m * m * n * n)
133+ else:
134+ wallace = 8 * k**4 * a**3 * b / (3 * np.pi**5 * m * m * n * n)
135+ assert abs(efficiency / wallace - 1) < 0.04, (m, n)
136+
137+
138+def test_thermoelastic_rate_of_a_uniform_plate_and_the_total():
139+ pl, ms = small_modes()
140+ rates = thermoelastic_rates(ms)
141+ omega = 2 * np.pi * ms.freqs
142+ expected = thermoelastic_loss_factor(STEEL_LOSSY, 1e-3, omega, ROOM_TEMPERATURE) * omega / 2
143+ assert np.abs(rates / expected - 1).max() < 1e-4
144+ U = bending_energy_per_node(pl, ms.w[3]).sum()
145+ assert abs(U / (omega[3] ** 2 / 2) - 1) < 1e-4
146+ assert np.allclose(ms.sigma, ms.sigma_material + ms.sigma_thermal + ms.sigma_radiation)
147+ assert (ms.sigma_material > 0).all() and (ms.sigma_radiation > 0).all() and (ms.sigma_thermal > 0).all()
148+ assert ((LN1000 / ms.sigma > 0.1) & (LN1000 / ms.sigma < 1000)).all()
149+ by = ms.t60_by_mechanism()
150+ assert set(by) == {"material", "thermal", "radiation"}
151+ assert isinstance(ms.table(5), str)
tests/test_modes.pyadded+83−0View file
@@ -0,0 +1,83 @@
1+"""The eigensolver: all solvers agree, mass normalization and residuals, degenerate
2+pairs of a disc, and the free square and disc frequencies against Gorman's and
3+Leissa's values."""
4+
5+import numpy as np
6+import pytest
7+
8+from cymsim import Plate, compute_modes, make_grid
9+from cymsim.linear.assemble import assemble
10+
11+from conftest import STEEL
12+
13+
14+def test_solvers_agree_and_shapes_are_mass_normalized():
15+ pl = Plate.uniform(make_grid(0.12, 0.12, 0.012, lambda x, y, r, t: r <= 0.06), STEEL, 1e-3)
16+ sys = assemble(pl)
17+ a = compute_modes(pl, 20000, solver="dense", damping=False, system=sys)
18+ b = compute_modes(pl, 20000, solver="banded", damping=False, system=sys)
19+ c = compute_modes(pl, 20000, solver="torch", device="cpu", damping=False, system=sys)
20+ assert a.count == b.count == c.count > 10
21+ assert np.allclose(a.freqs, b.freqs, rtol=1e-8) and np.allclose(a.freqs, c.freqs, rtol=1e-8)
22+ for m in (0, 1, a.count // 2, a.count - 1):
23+ phi = sys.from_grid(np.zeros(pl.grid.shape), np.zeros(pl.grid.shape), a.w[m])
24+ assert abs(phi @ (sys.mass * phi) - 1) < 1e-9
25+ lam = (2 * np.pi * a.freqs[m]) ** 2
26+ r = sys.K @ phi - lam * sys.mass * phi
27+ assert np.linalg.norm(r) < 1e-8 * np.linalg.norm(lam * sys.mass * phi)
28+ # The disc has degenerate pairs, which the staircase splits by about a percent.
29+ assert abs(a.freqs[1] / a.freqs[0] - 1) < 0.02
30+ # No rigid modes survive.
31+ assert a.freqs[0] > 1.0
32+
33+
34+def test_select_and_below():
35+ pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3)
36+ ms = compute_modes(pl, 20000, damping=False)
37+ sub = ms.below(5000)
38+ assert sub.count == (ms.freqs < 5000).sum() and sub.w.shape[0] == sub.count
39+ two = ms.select([0, 2])
40+ assert np.array_equal(two.freqs, ms.freqs[[0, 2]])
41+
42+
43+def test_free_square_frequencies_converge_to_gorman():
44+ # lambda = omega a^2 sqrt(rho h / D) for the completely free square (nu = 0.3).
45+ gorman = np.array([13.468, 19.596, 24.271, 34.801, 34.801, 61.093, 61.093])
46+ a, h = 0.3, 1e-3
47+ scale = a * a * np.sqrt(STEEL.rho * h / STEEL.bending_stiffness(h))
48+ errors = []
49+ for cells in (10, 20):
50+ ms = compute_modes(Plate.uniform(make_grid(a, a, a / cells), STEEL, h), 80 / scale / (2 * np.pi), damping=False)
51+ errors.append(np.abs(2 * np.pi * ms.freqs[:7] * scale / gorman - 1).max())
52+ assert errors[0] < 0.1 and errors[1] < 0.025 and errors[1] < errors[0]
53+
54+
55+def test_free_disc_frequencies_are_close_to_leissa():
56+ # lambda^2 = omega a^2 sqrt(rho h / D), completely free disc, nu = 0.3.
57+ leissa = np.array([5.253, 5.253, 9.084, 12.23, 12.23, 20.52, 20.52, 21.6, 21.6])
58+ a, h = 0.15, 1e-3
59+ scale = a * a * np.sqrt(STEEL.rho * h / STEEL.bending_stiffness(h))
60+ g = make_grid(2 * a, 2 * a, 2 * a / 40, lambda x, y, r, t: r <= a)
61+ ms = compute_modes(Plate.uniform(g, STEEL, h), 25 / scale / (2 * np.pi), damping=False)
62+ assert np.abs(2 * np.pi * ms.freqs[:9] * scale / leissa - 1).max() < 0.03
63+
64+
65+def test_save_and_load(tmp_path):
66+ pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3, name="small")
67+ ms = compute_modes(pl, 8000)
68+ ms.save(tmp_path / "m.npz")
69+ back = ms.__class__.load(tmp_path / "m.npz")
70+ assert np.array_equal(back.freqs, ms.freqs) and np.array_equal(back.w, ms.w)
71+ assert back.plate.name == "small" and back.plate.material == pl.material
72+ assert np.array_equal(back.sigma_radiation, ms.sigma_radiation)
73+
74+
75+@pytest.mark.slow
76+def test_crash_16_mode_count():
77+ from cymsim import presets
78+
79+ pl = presets.crash_16().evaluate(0.008)
80+ ms = compute_modes(pl, 16000)
81+ # About 1900 bending modes below 16 kHz on this grid (Weyl's estimate gives the same order).
82+ assert 1700 < ms.count < 2100
83+ assert ms.freqs[0] > 15 and ms.freqs[0] < 40
tests/test_plot.pyadded+29−0View file
@@ -0,0 +1,29 @@
1+"""The plot helpers draw without error (Agg backend)."""
2+
3+import matplotlib
4+
5+matplotlib.use("Agg")
6+
7+import numpy as np # noqa: E402
8+
9+from cymsim import Listener, Plate, Strike, compute_modes, make_grid, plot, presets, synthesize # noqa: E402
10+
11+from conftest import STEEL_LOSSY # noqa: E402
12+
13+
14+def test_plots():
15+ pl = Plate.uniform(make_grid(0.2, 0.15, 0.0125), STEEL_LOSSY, 1e-3)
16+ ms = compute_modes(pl, 4000)
17+ s = synthesize(ms, Strike(), Listener(), duration=0.1, movie=True)
18+ plot.plate(pl, strike=Strike(), listener=Listener())
19+ plot.plate(presets.crash_bell_16().evaluate(0.02))
20+ plot.mode(ms, 0)
21+ plot.frame(s, 3)
22+ plot.signal(s, tmax=0.05)
23+ plot.spectrogram(s, modes=ms, window=512, hop=128)
24+ plot.spectrogram(s.spectrogram(512, 128), log_freq=False)
25+ plot.decay_times(ms)
26+ plot.force(s)
27+ anim = plot.animate(s, stride=4)
28+ assert anim is not None
29+ assert np.isfinite(s.pressure).all()
tests/test_reference.pyadded+125−0View file
@@ -0,0 +1,125 @@
1+"""The discretization: rigid modes of the stiffness, the assembled K against the
2+linearized kernel, forces as the gradient of the energy, dissipativity of the
3+losses, and thickness scaling."""
4+
5+import numpy as np
6+import pytest
7+import torch
8+
9+from cymsim import Plate, make_grid
10+from cymsim.full.kernel import accelerations, damping, static_data
11+from cymsim.full.reference import max_angular_frequency, potential_energy
12+from cymsim.linear.assemble import assemble
13+
14+from conftest import BRONZE, STEEL
15+
16+T = lambda a: torch.as_tensor(np.asarray(a, dtype=np.float64)) # noqa: E731
17+
18+
19+def bowed(g, sag):
20+ Lx = (g.nx - 1) * g.hx
21+ return sag * (1 - (2 * g.X / Lx) ** 2) + 0.3 * sag * np.sin(3 * g.Y / Lx)
22+
23+
24+def random_fields(pl, amp, seed):
25+ r = np.random.default_rng(seed)
26+ return [np.where(pl.grid.mask, a * (r.random(pl.grid.shape) - 0.5), 0.0) for a in (0.02 * amp, 0.02 * amp, amp)]
27+
28+
29+@pytest.mark.parametrize("case", ["uniform", "tapered", "disc-hole"])
30+def test_stiffness_has_exactly_the_rigid_modes_in_its_null_space(case):
31+ if case == "disc-hole":
32+ g = make_grid(0.16, 0.16, 0.01, lambda x, y, r, t: (r <= 0.08) & (r >= 0.015))
33+ pl, fields = Plate.uniform(g, STEEL, 1e-3), 1
34+ else:
35+ g = make_grid(0.12, 0.09, 0.015)
36+ pl, fields = Plate.uniform(g, STEEL, 1e-3), 3
37+ if case == "tapered":
38+ pl = pl.with_thickness(1e-3 * (1 + 0.5 * (g.X + 0.06) / 0.12))
39+ sys = assemble(pl, fields)
40+ zero = np.zeros(g.shape)
41+ one = np.ones(g.shape)
42+ null = [(zero, zero, one), (zero, zero, g.X), (zero, zero, g.Y)]
43+ if fields == 3:
44+ null += [(one, zero, zero), (zero, one, zero), (-g.Y, g.X, zero)]
45+ scale = abs(sys.K[0, 0])
46+ for u, v, w in null:
47+ x = sys.from_grid(u, v, w)
48+ assert np.linalg.norm(sys.K @ x) < 1e-9 * scale * np.linalg.norm(x)
49+ x = sys.from_grid(zero, zero, np.sin(0.9 * np.arange(g.nx))[None, :] * np.cos(1.1 * np.arange(g.ny))[:, None])
50+ assert np.abs(sys.K @ x).max() > 1e-3
51+ assert abs(sys.K - sys.K.T).max() < 1e-9 * scale
52+
53+
54+@pytest.mark.parametrize("fields", [1, 3])
55+def test_assembled_K_equals_the_linearized_kernel_on_a_flat_plate(fields):
56+ pl = Plate.uniform(make_grid(0.2, 0.15, 0.02), BRONZE, 1e-3)
57+ assert pl.is_flat
58+ assert relative_error_against_kernel(pl, fields) < 1e-6
59+
60+
61+def test_assembled_K_equals_the_linearized_kernel_on_a_bowed_tapered_plate_with_a_hole():
62+ g = make_grid(0.2, 0.2, 0.0125, lambda x, y, r, t: (r <= 0.1) & (r >= 0.01))
63+ pl = Plate.uniform(g, BRONZE, 1e-3, profile=3e-3 * np.sin(20 * g.X) * np.cos(15 * g.Y))
64+ pl = pl.with_thickness(1e-3 * (1 + 0.4 * np.abs(g.X) / 0.1))
65+ assert not pl.is_flat
66+ assert assemble(pl).fields == 3
67+ assert relative_error_against_kernel(pl, 3) < 1e-6
68+
69+
70+def relative_error_against_kernel(pl, fields, seed=3):
71+ sys = assemble(pl, fields)
72+ x = np.random.default_rng(seed).random(sys.n) - 0.5
73+ y = sys.K @ x
74+ eps = 1e-9
75+ u, v, w = sys.to_grid(x * eps)
76+ S = static_data(pl, None, 1e-6, 0, 0, 0, "cpu", torch.float64)
77+ au, av, aw = (a.numpy() for a in accelerations(S, T(u), T(v), T(w)))
78+ m = pl.node_mass
79+ yref = sys.from_grid(-au * m / eps, -av * m / eps, -aw * m / eps)
80+ return np.linalg.norm(y - yref) / np.linalg.norm(yref)
81+
82+
83+@pytest.mark.parametrize("case", ["flat", "bowed", "bowed-tapered"])
84+def test_forces_are_the_gradient_of_the_energy(case):
85+ g = make_grid(0.1, 0.08, 0.01)
86+ pl = Plate.uniform(g, STEEL, 1e-3, profile=None if case == "flat" else bowed(g, 3e-3))
87+ if case == "bowed-tapered":
88+ pl = pl.with_thickness(1e-3 * (1 + 0.4 * (g.X + 0.05) / 0.1))
89+ f = random_fields(pl, 2.5e-3, 7)
90+ d = random_fields(pl, 1.0, 11)
91+ S = static_data(pl, None, 1e-6, 0, 0, 0, "cpu", torch.float64)
92+ acc = [a.numpy() for a in accelerations(S, *[T(a) for a in f])]
93+ h = 1e-7
94+ dU = (potential_energy(pl, *[a + h * b for a, b in zip(f, d)]) - potential_energy(pl, *[a - h * b for a, b in zip(f, d)])) / (2 * h)
95+ fd = -sum((pl.node_mass * a * b).sum() for a, b in zip(acc, d))
96+ assert abs(dU - fd) < 1e-6 * abs(dU)
97+
98+
99+def test_the_loss_terms_only_remove_energy():
100+ pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3)
101+ S = static_data(pl, None, 1e-6, 1.0, 5e-4, 20.0, "cpu", torch.float64)
102+ for seed in (1, 2, 3):
103+ vel = random_fields(pl, 1.0, seed)
104+ acc = [a.numpy() for a in damping(S, *[T(a) for a in vel])]
105+ power = sum((pl.node_mass * a * b).sum() for a, b in zip(acc, vel))
106+ assert power < 0
107+
108+
109+def test_doubling_the_thickness_doubles_the_bending_frequencies_only():
110+ g = make_grid(0.1, 0.1, 0.01)
111+ w1, w2 = max_angular_frequency(g, STEEL, 1e-3), max_angular_frequency(g, STEEL, 2e-3)
112+ # The largest frequency is set by the in-plane waves, independent of h.
113+ assert abs(w1 / w2 - 1) < 1e-6
114+ # It lies near 2 c_p / h for the grid-scale in-plane mode.
115+ assert 0.8 < w1 * g.hx / (2 * STEEL.in_plane_speed) < 1.5
116+ s1, s2 = assemble(Plate.uniform(g, STEEL, 1e-3)), assemble(Plate.uniform(g, STEEL, 2e-3))
117+ c = s1.n // 2
118+ assert abs(s2.K[c, c] / s1.K[c, c] - 8) < 1e-9
119+ assert abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12
120+
121+
122+def test_the_band_is_narrow_whichever_axis_is_longer():
123+ wide = assemble(Plate.uniform(make_grid(0.4, 0.1, 0.01), BRONZE, 1e-3))
124+ tall = assemble(Plate.uniform(make_grid(0.1, 0.4, 0.01), BRONZE, 1e-3))
125+ assert wide.bandwidth == tall.bandwidth <= 2 * 11 + 2
tests/test_synth.pyadded+97−0View file
@@ -0,0 +1,97 @@
1+"""The modal synthesis: a single mode follows the closed-form damped sinusoid, the
2+radiation gain of a piston equals the analytic Rayleigh sum, the attack and the
3+tail join, a lighter hit is proportionally quieter."""
4+
5+import numpy as np
6+
7+from cymsim import Listener, Plate, Strike, compute_modes, make_grid, synthesize
8+from cymsim.constants import LN1000, RHO_AIR
9+from cymsim.linear.synth import listener_weights, movie_frames, radiation_gains, simulate_contact, strike_gains, tail_start
10+
11+from conftest import STEEL
12+
13+
14+def small_model():
15+ pl = Plate.uniform(make_grid(0.2, 0.15, 0.0125), STEEL, 1e-3)
16+ return pl, compute_modes(pl, 4000, damping=False)
17+
18+
19+def test_single_mode_follows_the_closed_form():
20+ pl, ms = small_model()
21+ one = ms.select([0])
22+ one.sigma_material[:] = LN1000 # a decay time of exactly 1 s
23+ sigma = one.sigma
24+ assert abs(sigma[0] - LN1000) < 1e-9
25+ strike = Strike()
26+ g = strike_gains(one, strike)
27+ fps = 100
28+ c = simulate_contact(one, sigma, g, strike, fps, 48000)
29+ assert c.peak_force > 0 and c.contact_time > 0
30+ frames = movie_frames(one, sigma, c, fps, 10)
31+ w = 2 * np.pi * one.freqs[0]
32+ wd = np.sqrt(w * w - sigma[0] ** 2)
33+ j, i = pl.grid.ny // 2, pl.grid.nx // 3
34+ for f in range(1, 10):
35+ t = f / fps - c.t_end
36+ if t < 0:
37+ continue
38+ q = np.exp(-sigma[0] * t) * (c.amplitude[0].real * np.cos(wd * t) - c.amplitude[0].imag * np.sin(wd * t))
39+ assert abs(frames[f, j, i] - one.w[0, j, i] * q) < 1e-9 + 1e-6 * abs(one.w[0, j, i] * q)
40+
41+
42+def test_piston_radiation_gain_equals_the_analytic_sum():
43+ pl, ms = small_model()
44+ one = ms.select([0])
45+ phi = 1 / np.sqrt(pl.mass)
46+ one.w[0] = np.where(pl.grid.mask, phi, 0.0)
47+ one.freqs[0] = 1000.0
48+ listener = Listener(distance=1.0, angle=0)
49+ weight, delay, _ = listener_weights(pl, listener)
50+ G = radiation_gains(one, np.zeros(1), weight, delay)[0]
51+ w = 2 * np.pi * 1000
52+ sel = weight > 0
53+ expected = (weight[sel] * phi * np.exp(-1j * w * delay[sel])).sum()
54+ assert abs(G - expected) < 1e-9 * abs(expected)
55+ assert abs(abs(G) / (RHO_AIR * pl.area / (2 * np.pi) * phi) - 1) < 0.05
56+
57+
58+def test_attack_and_tail_join_and_the_tail_decays():
59+ pl, ms = small_model()
60+ ms.sigma_material[:] = 3e-4 * ms.omega / 2
61+ strike, listener = Strike(velocity=2), Listener()
62+ fs = 48000
63+ s = synthesize(ms, strike, listener, duration=0.2, fs=fs, movie=True)
64+ assert len(s.pressure) == round(0.2 * fs) and s.frames.shape == (24, pl.grid.ny, pl.grid.nx)
65+ g = strike_gains(ms, strike)
66+ c = simulate_contact(ms, ms.sigma, g, strike, 120, fs)
67+ _, _, max_delay = listener_weights(pl, listener)
68+ n0 = int(np.ceil(tail_start(c, max_delay) * fs - 1e-9))
69+ steps = np.abs(np.diff(s.pressure[n0 + 1 : n0 + 200]))
70+ jump = abs(s.pressure[n0] - s.pressure[n0 - 1])
71+ assert jump < 3 * steps.max()
72+ rms = lambda a: np.sqrt((a.astype(float) ** 2).mean()) # noqa: E731
73+ assert rms(s.pressure[-960:]) < rms(s.pressure[:960])
74+ assert s.info["peak_force"] == c.peak_force and s.peak_spl > 60
75+
76+
77+def test_a_lighter_hit_is_proportionally_quieter():
78+ _, ms = small_model()
79+ a = synthesize(ms, Strike(velocity=0.2), duration=0.05).pressure
80+ b = synthesize(ms, Strike(velocity=0.1), duration=0.05).pressure
81+ ratio = np.abs(a).max() / np.abs(b).max()
82+ # The Hertz contact is nonlinear (a faster hit is shorter and brighter), so the ratio exceeds 2 somewhat.
83+ assert 1.8 < ratio < 3.5
84+
85+
86+def test_sound_outputs(tmp_path):
87+ _, ms = small_model()
88+ s = synthesize(ms, duration=0.1)
89+ s.to_wav(tmp_path / "a.wav")
90+ import wave
91+
92+ with wave.open(str(tmp_path / "a.wav")) as w:
93+ assert w.getnchannels() == 1 and w.getframerate() == 48000 and w.getnframes() == 4800
94+ assert np.abs(s.samples("peak")).max() <= 1 and abs(s.peak_dbfs("peak") + 1) < 0.01
95+ sp = s.spectrogram(512, 128)
96+ assert sp.db.shape == (257, len(sp.times)) and sp.db.max() == 0
97+ assert "linear" in repr(s)