3688ceacymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman modelJeremy Magland 1"""The full model: the time step, a free step from a smooth state agrees between
2float32 and float64, energy is conserved without losses, the recorded audio is
3the Rayleigh sum, and a light hit agrees with the linear model."""
5import numpy as np
6import pytest
7import torch
9from cymsim import Listener, Plate, Simulation, Strike, compute_modes, make_grid, simulate, synthesize
10from cymsim.constants import CFL
11from cymsim.full.reference import potential_energy
13from conftest import STEEL, STEEL_LOSSY
16def 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)
26def 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
35def 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()
54def 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)
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)
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]
84def 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()
98def 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)
110def 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
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))))
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