1"""The eigensolver: all solvers agree, mass normalization and residuals, degenerate
2pairs of a disc, and the free square and disc frequencies against Gorman's and
3Leissa's values."""
5import numpy as np
6import pytest
8from cymsim import Plate, compute_modes, make_grid
9from cymsim.linear.assemble import assemble
11from conftest import STEEL
14def 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
34def 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]])
43def 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]
55def 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
65def 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)
76def test_crash_16_mode_count():
77 from cymsim import presets
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