1"""Configuration dataclasses and their YAML representation.
3Every physical and numerical parameter of a run lives in one of the dataclasses below, grouped
4by what it describes: the material, the cymbal profile, the seeded imperfections, the mesh, the
5list of hits (each with its own exciter), the listener, the numerics and the output. Defaults
6reproduce the 18-inch medium-thin crash of the original prototype. A SimulationConfig can be
7built directly in Python or loaded from a YAML file whose keys mirror the field names; a hit's
8exciter is chosen by its ``kind`` key ("stick", "mallet" or "force").
9"""
10from __future__ import annotations
12import dataclasses
13import math
14import types
15import typing
16from dataclasses import dataclass, field
17from pathlib import Path
19import yaml
22# ------------------------------------------------------------------------------- material
24class Material:
25 """Isotropic material of the shell (default: B20 bronze) and its two loss parameters."""
26 E: float = 1.1e11 # Young's modulus [Pa]
27 nu: float = 0.34 # Poisson ratio
28 rho: float = 8700.0 # density [kg/m^3]
29 sigma0: float = 0.9 # frequency-independent loss [1/s]; T60 of the lows ~ 6.9/sigma0
30 sigma1: float = 5e-4 # loss proportional to omega [m^2/s]; sets the ceiling of the cascade
33# -------------------------------------------------------------------------------- profile
35class Profile:
36 """Axisymmetric part of the cymbal: bow, bell and thickness taper. Mean rim height is z0 = 0."""
37 radius: float = 0.2286 # outer radius [m] (18 inch)
38 bell_height: float = 0.013 # bell height above the bow at the shoulder [m]
39 bell_top_radius: float = 0.025 # radius of the flat bell plateau [m]
40 bell_shoulder_radius: float = 0.080 # mean radius of the bell/bow junction [m]
41 bow_height: float = 0.030 # bow rise from the rim to the shoulder [m]
42 bow_power: float = 1.4 # bow shape exponent: 1 = cone, > 1 = curved, flat at the shoulder
43 h_rim: float = 1.0e-3 # thickness at the rim [m]
44 h_shoulder: float = 1.7e-3 # thickness at the shoulder [m]
45 h_bell: float = 2.6e-3 # thickness in the bell [m]
46 max_slope: float = 0.4 # warn when max |grad z0| exceeds this (shallow-shell bound)
50class Imperfections:
51 """Seeded symmetry breaking: hammer dents, rim warp, shoulder ovality and a bell offset."""
52 seed: int = 1
53 n_dents: int = 400 # Gaussian hammer marks on the bow, random sign
54 dent_depth: float = 0.25e-3 # typical dent depth [m] (x 0.5 to 1.5 jitter)
55 dent_sigma: float = 7e-3 # dent Gaussian sigma [m] (x 0.7 to 1.3 jitter); keep >= element size
56 dent_margin: float = 0.01 # dents stay this far from the shoulder and the rim [m]
57 rim_warp: float = 0.6e-3 # m = 2 rim-warp amplitude [m]; mode m gets 2 * rim_warp / m
58 rim_warp_modes: tuple[int, ...] = (2, 3, 4, 5)
59 shoulder_ovality: float = 0.03 # relative m = 2 modulation of the shoulder radius
60 bell_offset: float = 1.5e-3 # bell centre offset from the geometric centre [m]
64class MeshConfig:
65 element_size: float = 6e-3 # target edge length of the triangles [m]
66 smoothing: int = 1 # Laplacian smoothing passes; more improves shape, shrinks the smallest elements
69# ------------------------------------------------------------------------------- exciters
71class StickExciter:
72 """Hard stick: Hertz contact (alpha = 1.5) along a short line, as in a shoulder crash."""
73 kind: str = "stick"
74 mass: float = 0.030 # effective mass [kg]
75 k: float = 5e8 # total contact stiffness [N/m^alpha]
76 alpha: float = 1.5 # contact exponent
77 patch_length: float = 0.03 # length of the contact line [m]
78 patch_points: int = 7 # springs along the line
79 patch_angle: float = 0.6 # direction of the line relative to the radial direction [rad]
80 gap: float = 1e-4 # initial clearance above the surface at release [m]
84class MalletExciter:
85 """Soft felt mallet: power-law contact on a small disk; k is set from a target compression."""
86 kind: str = "mallet"
87 mass: float = 0.30 # head plus effective arm mass [kg]
88 alpha: float = 2.5 # felt exponent (1.5 = Hertz, 2 to 3 = felt)
89 compression: float = 12e-3 # felt compression at full stop on a rigid wall [m]
90 radius: float = 0.012 # head contact radius [m]; 7 springs (centre + hexagon)
91 gap: float = 1e-4
95class ForceExciter:
96 """Prescribed force at one point, independent of the motion: a raised cosine
97 F(t) = peak (1 - cos(2 pi t / duration)) / 2 for 0 <= t <= duration after the hit time,
98 pushing the surface down. There is no contact and no rebound. This is the excitation of
99 modal plate models such as Ducceschi and Touze (2015), and is meant for comparison with them;
100 ``Hit.velocity`` is ignored."""
101 kind: str = "force"
102 peak: float = 100.0 # peak force [N]
103 duration: float = 2e-3 # total length of the pulse [s]
106Exciter = typing.Union[StickExciter, MalletExciter, ForceExciter]
107EXCITER_KINDS: dict[str, type] = {"stick": StickExciter, "mallet": MalletExciter, "force": ForceExciter}
111class Hit:
112 """One strike: when, where, how fast, and with what."""
113 time: float = 0.0 # release time of the exciter [s]
114 r: float = 0.137 # radial position of the contact centre [m]
115 phi: float = 0.3 # azimuth of the contact centre [rad]
116 velocity: float = 8.0 # approach speed [m/s]
117 exciter: Exciter = field(default_factory=MalletExciter)
120# ------------------------------------------------------------------------------- listener
122class Listener:
123 """Listener point and air; the mean rim plane is z = 0 with the cymbal centre at the origin."""
124 position: tuple[float, float, float] = (0.40, 0.30, 0.90) # [m]
125 c_air: float = 343.0
126 rho_air: float = 1.2
127 back_diffraction_gain: float = 1.0 # weight of the diffracted back-side term (0 = baffled only)
130# ------------------------------------------------------------------------------- numerics
132class Numerics:
133 fs: int = 44100 # audio sample rate [Hz]; the sub-step rate is an integer multiple
134 duration: float = 2.0 # simulated time [s]
135 cfl_safety: float = 0.6 # fraction of the linear stability limit actually used
136 stability_bound: str = "global" # "global": lambda_max(M^-1 K) by power iteration; "element": the Irons element bound (larger, used before 2026-10)
137 rot_inertia_factor: float = 1 / 64 # lumped rotational inertia I = m (h^2/12 + f A_e); see docs
138 inplane_mass_scale: float = 1.0 # > 1 slows in-plane waves and relaxes the time step; see docs
139 steps_per_block: int = 20000 # sub-steps per compiled block
140 dtype: str = "float32" # simulation precision: "float32" or "float64"
141 nonlinear: bool = True # False drops the quadratic strain terms (linear reference run)
142 tangent: bool = False # also integrate the tangent: d/d delta of the run when every prescribed peak force is x (1 + delta)
143 n_progress: int = 20 # number of progress/energy reports during the run
147class Output:
148 prefix: str = "crash" # output files are <prefix>_<name>.wav / .png / .npz
149 wav: bool = True
150 npz: bool = True
151 cascade_plot: bool = True # needs matplotlib
152 record_acceleration: bool = False # keep the normal acceleration of every node at the audio rate
153 displacement_rate: float = 60.0 # frames per second of transverse displacement kept for display; 0 = none
154 probes: list[tuple[float, float]] = field(default_factory=list) # points (x, y) [m] where the normal acceleration is kept at the audio rate
157# The fields in which a continued run may differ from the run it continues: the length of the run,
158# how often energy is reported, and the naming of the output files. Two configurations that agree
159# on everything else describe runs that are prefixes of one another (Simulator.run(resume=...)),
160# and share a cache key (cache.py).
161CONTINUATION_FIELDS = (("numerics", "duration"), ("numerics", "n_progress"), ("output", "prefix"),
162 ("output", "wav"), ("output", "npz"), ("output", "cascade_plot"))
165# ----------------------------------------------------------------------------- simulation
167class SimulationConfig:
168 material: Material = field(default_factory=Material)
169 profile: Profile = field(default_factory=Profile)
170 imperfections: Imperfections = field(default_factory=Imperfections)
171 mesh: MeshConfig = field(default_factory=MeshConfig)
172 hits: list[Hit] = field(default_factory=lambda: [Hit()])
173 listener: Listener = field(default_factory=Listener)
174 numerics: Numerics = field(default_factory=Numerics)
175 output: Output = field(default_factory=Output)
177 # ---- validation
178 def validate(self) -> None:
179 p, n = self.profile, self.numerics
180 if not self.hits:
181 raise ValueError("at least one hit is required")
182 if n.dtype not in ("float32", "float64"):
183 raise ValueError(f"numerics.dtype must be 'float32' or 'float64', got {n.dtype!r}")
184 if n.stability_bound not in ("global", "element"):
185 raise ValueError(f"numerics.stability_bound must be 'global' or 'element', got {n.stability_bound!r}")
186 if not (0 < n.cfl_safety <= 1):
187 raise ValueError("numerics.cfl_safety must lie in (0, 1]")
188 if self.mesh.element_size <= 0 or self.mesh.element_size > p.radius / 2:
189 raise ValueError("mesh.element_size must be positive and well below the radius")
190 if n.tangent and self.output.record_acceleration:
191 raise ValueError("numerics.tangent does not support output.record_acceleration")
192 if n.tangent and not any(isinstance(h.exciter, ForceExciter) for h in self.hits):
193 raise ValueError("numerics.tangent differentiates the prescribed peak forces: it needs a force exciter")
194 if self.output.displacement_rate < 0:
195 raise ValueError("output.displacement_rate must be non-negative")
196 for i, h in enumerate(self.hits):
197 if h.time < 0:
198 raise ValueError(f"hit {i}: time must be non-negative")
199 if h.velocity <= 0:
200 raise ValueError(f"hit {i}: velocity must be positive (the exciter approaches from above)")
201 ex = h.exciter
202 extent = (ex.patch_length / 2 if isinstance(ex, StickExciter)
203 else ex.radius if isinstance(ex, MalletExciter) else 0.0)
204 if h.r + extent >= p.radius - self.mesh.element_size:
205 raise ValueError(f"hit {i}: contact patch reaches the rim (r = {h.r}, extent = {extent})")
206 if ex.kind not in EXCITER_KINDS:
207 raise ValueError(f"hit {i}: unknown exciter kind {ex.kind!r}")
208 if isinstance(ex, ForceExciter) and (ex.peak < 0 or ex.duration <= 0):
209 raise ValueError(f"hit {i}: a force pulse needs peak >= 0 and duration > 0")
210 for i, (x, y) in enumerate(self.output.probes):
211 if math.hypot(x, y) >= p.radius:
212 raise ValueError(f"output.probes[{i}] = ({x}, {y}) lies outside the cymbal")
214 # ---- dict / YAML
215 def to_dict(self) -> dict:
216 return _plain(dataclasses.asdict(self))
218 def identity_dict(self) -> dict:
219 """``to_dict()`` without CONTINUATION_FIELDS: what a run and its continuation share."""
220 d = self.to_dict()
221 for section, name in CONTINUATION_FIELDS:
222 d[section].pop(name, None)
223 return d
226 def from_dict(cls, d: dict) -> "SimulationConfig":
227 return _from_dict(cls, d)
229 def to_yaml(self, path: str | Path) -> None:
230 Path(path).write_text(yaml.safe_dump(self.to_dict(), sort_keys=False))
233 def from_yaml(cls, path: str | Path) -> "SimulationConfig":
234 with open(path) as f:
235 d = yaml.safe_load(f) or {}
236 return cls.from_dict(d)
239# --------------------------------------------------------------- generic dataclass <-> dict
240def dict_differences(a, b, path: str = "") -> list[str]:
241 """Dotted paths at which two plain nested structures differ."""
242 if isinstance(a, dict) and isinstance(b, dict):
243 return [d for k in sorted(set(a) | set(b), key=str)
244 for d in dict_differences(a.get(k), b.get(k), f"{path}.{k}" if path else str(k))]
245 return [] if a == b else [path or "(root)"]
248def _plain(obj):
249 """Recursively turn tuples into lists so that yaml.safe_dump accepts the structure."""
250 if isinstance(obj, dict):
251 return {k: _plain(v) for k, v in obj.items()}
252 if isinstance(obj, (list, tuple)):
253 return [_plain(v) for v in obj]
254 return obj
257def _build(tp, value):
258 """Construct a value of annotated type ``tp`` from plain YAML data."""
259 origin = typing.get_origin(tp)
260 if isinstance(tp, type) and dataclasses.is_dataclass(tp):
261 return _from_dict(tp, value)
262 if origin is list:
263 (item,) = typing.get_args(tp)
264 return [_build(item, v) for v in value]
265 if origin is tuple:
266 return tuple(value)
267 if origin in (typing.Union, getattr(types, "UnionType", typing.Union)): # UnionType: Python >= 3.10
268 # The only union in the schema is the exciter, discriminated by its "kind" key.
269 if isinstance(value, dict) and "kind" in value:
270 try:
271 return _from_dict(EXCITER_KINDS[value["kind"]], value)
272 except KeyError:
273 raise ValueError(f"unknown exciter kind {value['kind']!r}; "
274 f"expected one of {sorted(EXCITER_KINDS)}") from None
275 if any(isinstance(value, t) for t in typing.get_args(tp)):
276 return value
277 raise ValueError("an exciter mapping needs a 'kind' key ('stick', 'mallet' or 'force')")
278 if tp is float and isinstance(value, (int, float)):
279 return float(value)
280 return value
283def _from_dict(cls, d):
284 if isinstance(d, cls):
285 return d
286 if not isinstance(d, dict):
287 raise ValueError(f"expected a mapping for {cls.__name__}, got {type(d).__name__}")
288 hints = typing.get_type_hints(cls)
289 names = {f.name for f in dataclasses.fields(cls)}
290 unknown = set(d) - names
291 if unknown:
292 raise ValueError(f"unknown field(s) for {cls.__name__}: {sorted(unknown)}")
293 return cls(**{k: _build(hints[k], v) for k, v in d.items()})
296def default_config(exciter: str = "mallet") -> SimulationConfig:
297 """The prototype's default run for the given exciter kind."""
298 if exciter == "mallet":
299 hit = Hit(time=0.0, r=0.60 * Profile.radius, phi=0.3, velocity=8.0, exciter=MalletExciter())
300 elif exciter == "stick":
301 hit = Hit(time=0.0, r=0.85 * Profile.radius, phi=0.3, velocity=7.0, exciter=StickExciter())
302 else:
303 raise ValueError(f"unknown exciter kind {exciter!r}")
304 return SimulationConfig(hits=[hit])
307def t60_estimate(material: Material, frequency: float) -> float:
308 """Decay time to -60 dB implied by sigma0 and sigma1 at a given frequency, for a flat plate
309 with bending wavenumber k^2 = omega / sqrt(D / rho h). Rough, but useful for tuning losses."""
310 # The sigma1 term acts as 2 sigma1 k^2; for a thin plate k^2 = omega sqrt(rho h / D) which
311 # depends on thickness; we quote the value for h = 1 mm.
312 h = 1e-3
313 D = material.E * h**3 / (12 * (1 - material.nu**2))
314 omega = 2 * math.pi * frequency
315 k2 = omega * math.sqrt(material.rho * h / D)
316 sigma = material.sigma0 + material.sigma1 * k2
317 return 6.91 / sigma