ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 1"""Reference implementation (Python, cymsim) of the modal precompute that cymplay does in the browser.
3Used only to validate the browser numbers; the app itself needs no Python.
5The in-plane DOFs (u, v) of the linearised shell are eliminated by static condensation, which
6leaves a dense stiffness on the (w, theta_x, theta_y) DOFs. Its mass-normalised eigenmodes below
7``--fmax`` give, for each mode m: the frequency f_m, the decay rate sigma_m (sigma0 plus the
8projected sigma1 term), the complex listener gains (front Rayleigh sum and diffracted back sum,
9with the element delays folded into a phase at f_m), and the transverse mode shape on the nodes,
10quantised to int16 with one scale per mode. The mesh, the nodal surface height, the full config
11and the exciter defaults are stored alongside so that the app can locate strike points, run the
12contact law and write a config back out.
14Output format: magic ``CYMB``, a uint32 little-endian header length, a UTF-8 JSON header, then
15the raw little-endian arrays at the byte offsets the header lists.
17Usage (on a machine with enough memory for a dense (3N)^2 float64 matrix, 1.5 GB at 6 mm):
18 python make_modes.py ../cym2/examples/crash_mallet.yaml public/modes/crash_mallet.bin
19"""
20from __future__ import annotations
22import argparse
23import dataclasses
24import json
25import struct
26import sys
27import time
29import numpy as np
30import scipy.linalg as sla
31import scipy.sparse as sp
32import scipy.sparse.linalg as spla
34try:
35 import shim # noqa: F401 (office machine: Python 3.9 compatibility, adds cym2/src to the path)
36except ImportError:
37 pass
39import jax
40import jax.numpy as jnp
42jax.config.update("jax_enable_x64", True)
44from cymsim.config import MalletExciter, SimulationConfig, StickExciter # noqa: E402
45from cymsim.elements import element_energy # noqa: E402
46from cymsim.geometry import CymbalProfile # noqa: E402
47from cymsim.mesh import disk_mesh # noqa: E402
48from cymsim.radiation import RadiationModel # noqa: E402
49from cymsim.shell import N_DOF, ShellModel # noqa: E402
51FORMAT_VERSION = 2
54def log(msg):
55 print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
58def build(cfg: SimulationConfig, fmax: float, condense_rotations: bool = False):
59 t0 = time.time()
60 prof = CymbalProfile(cfg.profile, cfg.imperfections)
61 mesh = disk_mesh(cfg.profile.radius, cfg.mesh.element_size, cfg.mesh.smoothing)
62 sh = ShellModel(mesh, prof, cfg.material, cfg.numerics)
63 N, E, tri = mesh.n_nodes, mesh.n_elements, mesh.tri
64 log(f"mesh: {N} nodes, {E} elements; shell built in {time.time() - t0:.0f} s")
66 # ---- global tangent stiffness at rest from the element Hessians (float64)
67 Kh = np.asarray(jax.vmap(jax.hessian(element_energy), in_axes=(0, 0, 0, 0, 0, 0, 0, 0, None))(
68 jnp.zeros((E, 3, N_DOF)), jnp.asarray(mesh.Gx), jnp.asarray(mesh.Gy), jnp.asarray(mesh.area),
69 jnp.asarray(sh.z0x), jnp.asarray(sh.z0y), jnp.asarray(sh.Cm_e), jnp.asarray(sh.Ke),
70 sh.nl)).reshape(E, 15, 15)
71 gi = (tri[:, :, None] * 5 + np.arange(5)[None, None, :]).reshape(E, 15)
72 rows = np.repeat(gi, 15, axis=1).ravel()
73 cols = np.tile(gi, (1, 15)).ravel()
74 K = sp.coo_matrix((Kh.ravel(), (rows, cols)), shape=(5 * N, 5 * N)).tocsr()
75 M = sh.M.ravel()
76 log(f"assembled K: {K.nnz} nonzeros")
78 # ---- static condensation of the in-plane DOFs. Three in-plane DOFs are pinned to remove
79 # the in-plane rigid motions; the condensed operator then has exactly three zero modes
80 # (transverse translation and the two tilts), which we drop.
81 if condense_rotations:
82 ip = np.sort(np.concatenate([np.arange(N) * 5 + d for d in (0, 1, 3, 4)]))
83 else:
84 ip = np.sort(np.concatenate([np.arange(N) * 5, np.arange(N) * 5 + 1]))
85 wb = np.setdiff1d(np.arange(5 * N), ip)
86 # pinned in-plane DOFs u0, v0, v1 in the local numbering of ``ip``
87 pins = [int(np.searchsorted(ip, g)) for g in (0, 1, 5 + 1)]
88 Kuu = K[ip][:, ip].tolil()
89 Kuw = K[ip][:, wb].toarray()
90 Kww = K[wb][:, wb].toarray()
91 for p in pins: # u0, v0, v1
92 Kuu[p, :] = 0
93 Kuu[:, p] = 0
94 Kuu[p, p] = 1.0
95 Kuw[p, :] = 0
96 t0 = time.time()
97 lu = spla.splu(Kuu.tocsc())
98 X = lu.solve(Kuw) # (2N, 3N)
99 S = Kww - Kuw.T @ X
100 S = 0.5 * (S + S.T)
101 del Kww, Kuw, X
102 log(f"condensed stiffness {S.shape} in {time.time() - t0:.0f} s")
104 # ---- mass-normalised eigenmodes below fmax (with some margin, trimmed afterwards)
105 Mw = M[wb]
106 s = 1.0 / np.sqrt(Mw)
107 S *= s[:, None]
108 S *= s[None, :]
109 t0 = time.time()
110 lam, V = sla.eigh(S, subset_by_value=[-np.inf, (2 * np.pi * fmax) ** 2], driver="evr",
111 overwrite_a=True, check_finite=False)
112 del S
113 f = np.sqrt(np.maximum(lam, 0)) / (2 * np.pi)
114 keep = f > 1.0 # drop the rigid modes
115 f, V = f[keep], V[:, keep]
116 Phi = V * s[:, None] # Phi^T M Phi = I on the (w, thx, thy) DOFs
117 Nm = len(f)
118 log(f"eigh: {Nm} modes between {f[0]:.1f} Hz and {f[-1]:.0f} Hz in {time.time() - t0:.0f} s; "
119 f"dropped {int((~keep).sum())} rigid modes")
120 stride = 1 if condense_rotations else 3
121 Phi_w = np.ascontiguousarray(Phi[0::stride, :]) # (N, Nm), the w rows
123 # ---- modal damping: sigma_m = sigma0 + sigma1 phi^T L phi with L = sum_e (rho h A)_e G^T G
124 rhoA = cfg.material.rho * sh.h_e * mesh.area
125 Lrows, Lcols, Lvals = [], [], []
126 for i in range(3):
127 for j in range(3):
128 Lrows.append(tri[:, i]); Lcols.append(tri[:, j])
129 Lvals.append(rhoA * (mesh.Gx[:, i] * mesh.Gx[:, j] + mesh.Gy[:, i] * mesh.Gy[:, j]))
130 L = sp.coo_matrix((np.concatenate(Lvals), (np.concatenate(Lrows), np.concatenate(Lcols))),
131 shape=(N, N)).tocsr()
132 damp_q = np.einsum("nm,nm->m", Phi_w, L @ Phi_w)
133 sigma = cfg.material.sigma0 + cfg.material.sigma1 * damp_q
134 log(f"damping: sigma from {sigma.min():.2f} to {sigma.max():.1f} 1/s "
135 f"(T60 {6.91 / sigma.max() * 1e3:.0f} ms at the top mode)")
137 # ---- listener gains: Rayleigh front sum and diffracted back sum per mode, with the element
138 # delays as phases at the mode frequency. pressure = Re(G_m * complex acceleration).
139 rad = RadiationModel(mesh, prof, sh.z0_node, cfg.listener, 1e-6, 1000, jnp.float64)
140 phi_e = Phi_w[tri].mean(1) # (E, Nm) element-mean mode shape
141 c = cfg.listener.c_air
142 omega = 2 * np.pi * f
143 G_front = np.einsum("em,e,em->m", phi_e, np.asarray(rad.coef_f),
144 np.exp(-1j * np.outer(rad.r_front, omega) / c))
145 G_back = np.einsum("em,e,em->m", phi_e, np.asarray(rad.coef_b),
146 np.exp(-1j * np.outer(rad.r_back, omega) / c))
147 log("radiation gains done")
149 # ---- quantise the mode shapes: int16 with one scale per mode
150 scale = np.abs(Phi_w).max(0) / 32767.0
151 Phi_q = np.round(Phi_w / scale[None, :]).astype(np.int16)
153 arrays = {
154 "xy": mesh.xy.astype(np.float32), # (N, 2)
155 "tri": tri.astype(np.int32), # (E, 3)
156 "z0": sh.z0_node.astype(np.float32), # (N,)
157 "h": sh.h_node.astype(np.float32), # (N,)
158 "freq": f.astype(np.float64), # (Nm,)
159 "damp_q": damp_q.astype(np.float64), # (Nm,) sigma = sigma0 + sigma1 * damp_q
160 "gain_front": np.stack([G_front.real, G_front.imag], 1).astype(np.float32), # (Nm, 2)
161 "gain_back": np.stack([G_back.real, G_back.imag], 1).astype(np.float32), # (Nm, 2)
162 "phi_scale": scale.astype(np.float32), # (Nm,)
163 "phi_q": Phi_q, # (N, Nm) int16, node-major
164 }
165 meta = {
166 "format_version": FORMAT_VERSION,
167 "n_nodes": int(N), "n_elements": int(E), "n_modes": int(Nm), "fmax": float(fmax),
168 "radius": float(cfg.profile.radius),
169 "E0": float(cfg.material.E), "rho0": float(cfg.material.rho),
170 "mass": float(sh.mass),
171 "config": cfg.to_dict(),
172 "playback_element_size": float(cfg.mesh.element_size),
173 "condense_rotations": bool(condense_rotations),
174 "exciters": {"stick": dataclasses.asdict(StickExciter()),
175 "mallet": dataclasses.asdict(MalletExciter())},
176 "notes": "phi = phi_q * phi_scale[m] is the mass-normalised transverse mode shape on the "
177 "nodes; modal EOM: a'' + 2 sigma a' + (2 pi f)^2 a = phi(p) . F; pressure = "
178 "Re(sum_m G_m * complex acceleration_m), unbaffled = front - back.",
179 }
180 return meta, arrays
183def write_bin(path: str, meta: dict, arrays: dict[str, np.ndarray]) -> None:
184 offset = 0
185 entries = {}
186 blobs = []
187 for name, a in arrays.items():
188 a = np.ascontiguousarray(a)
189 if a.dtype.byteorder == ">":
190 a = a.astype(a.dtype.newbyteorder("<"))
191 b = a.tobytes()
192 entries[name] = {"dtype": a.dtype.name, "shape": list(a.shape), "offset": offset, "nbytes": len(b)}
193 blobs.append(b)
194 offset += len(b)
195 header = json.dumps({**meta, "arrays": entries}).encode("utf-8")
196 with open(path, "wb") as fh:
197 fh.write(b"CYMB")
198 fh.write(struct.pack("<I", len(header)))
199 fh.write(header)
200 for b in blobs:
201 fh.write(b)
202 log(f"wrote {path}: {(8 + len(header) + offset) / 1e6:.1f} MB")
205def read_bin(path: str):
206 """Reader for the same format (used by the validation script)."""
207 with open(path, "rb") as fh:
208 data = fh.read()
209 assert data[:4] == b"CYMB", "not a cymplay modes file"
210 (hlen,) = struct.unpack("<I", data[4:8])
211 meta = json.loads(data[8:8 + hlen].decode("utf-8"))
212 base = 8 + hlen
213 arrays = {}
214 for name, e in meta.pop("arrays").items():
215 arrays[name] = np.frombuffer(data, dtype=np.dtype(e["dtype"]).newbyteorder("<"),
216 count=int(np.prod(e["shape"])), offset=base + e["offset"]).reshape(e["shape"])
217 return meta, arrays
220def main(argv=None):
221 ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
222 ap.add_argument("config", help="cymsim YAML config")
223 ap.add_argument("out", help="output .bin file")
224 ap.add_argument("--fmax", type=float, default=20000.0, help="keep modes below this frequency [Hz]")
225 ap.add_argument("--element-size", type=float, help="override mesh.element_size [m]")
226 ap.add_argument("--condense-rotations", action="store_true",
227 help="also condense theta_x, theta_y (w-only modes, as the browser does)")
228 ap.add_argument("--no-imperfections", action="store_true",
229 help="zero every seeded imperfection (deterministic geometry, comparable across RNGs)")
230 args = ap.parse_args(argv)
231 cfg = SimulationConfig.from_yaml(args.config)
232 if args.element_size:
233 cfg.mesh.element_size = args.element_size
234 if args.no_imperfections:
235 im = cfg.imperfections
236 im.n_dents, im.rim_warp, im.shoulder_ovality, im.bell_offset = 0, 0.0, 0.0, 0.0
237 cfg.validate()
238 meta, arrays = build(cfg, args.fmax, args.condense_rotations)
239 write_bin(args.out, meta, arrays)
242if __name__ == "__main__":
243 sys.exit(main())