3688ceacymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman modelJeremy Magland 1"""The modes of a plate up to a cutoff frequency, and their loss rates.
3The generalized problem K phi = omega^2 M phi with the diagonal lumped mass M
4is turned into the standard symmetric problem A = M^{-1/2} K M^{-1/2}, whose
5eigenvalues in the band of interest are computed by LAPACK, on the CPU through
6SciPy (a dense solve restricted to the band; about 4 s for a 16 inch cymbal at
78 mm, and 40 s for the same cymbal with a profile, whose system is three times
8as large) or on a GPU through PyTorch (a full dense eigendecomposition in
9double precision, the faster path for the larger systems when a CUDA device
10is available). The mode shapes
11returned are M^{-1/2} times the eigenvectors of A, that is, mass-normalized.
12The zero-frequency rigid modes of the free plate are dropped.
13"""
15from __future__ import annotations
17import math
18import time
19from dataclasses import dataclass, replace
20from typing import Literal
22import numpy as np
23import scipy.linalg
25from ..constants import LN1000
26from ..plate import Plate
27from .assemble import LinearSystem, assemble
29Solver = Literal["auto", "banded", "dense", "torch"]
33class Modes:
34 """The modes of a plate below a cutoff frequency, with their per-mode decay rates.
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 """
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]."""
60 def count(self) -> int:
61 return len(self.freqs)
64 def omega(self) -> np.ndarray:
65 return 2 * math.pi * self.freqs
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)
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
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
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 )
100 def below(self, fmax: float) -> "Modes":
101 """The modes below a frequency [Hz]."""
102 return self.select(self.freqs < fmax)
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)
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
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 )
142 def load(cls, path) -> "Modes":
143 from ..grid import Grid
144 from ..material import Material, Thermal
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"]))
155 def __repr__(self) -> str:
156 return f"Modes({self.count} modes of {self.plate.name!r} up to {self.fmax:.0f} Hz)"
159def _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}"
169def 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]``.
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
224 modes.sigma_radiation = radiation_rates(modes)
225 modes.sigma_thermal = thermoelastic_rates(modes)
226 modes.seconds = time.perf_counter() - t0
227 return modes
230def _cuda_available() -> bool:
231 try:
232 import torch
234 return torch.cuda.is_available()
235 except ImportError:
236 return False
239def _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
245def _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
256def _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
262def _eig_torch(sys: LinearSystem, lo: float, hi: float, device):
263 import torch
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()