3688ceacymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman modelJeremy Magland 1"""The discretization: rigid modes of the stiffness, the assembled K against the
2linearized kernel, forces as the gradient of the energy, dissipativity of the
3losses, and thickness scaling."""
5import numpy as np
6import pytest
7import torch
9from cymsim import Plate, make_grid
10from cymsim.full.kernel import accelerations, damping, static_data
11from cymsim.full.reference import max_angular_frequency, potential_energy
12from cymsim.linear.assemble import assemble
14from conftest import BRONZE, STEEL
16T = lambda a: torch.as_tensor(np.asarray(a, dtype=np.float64)) # noqa: E731
19def bowed(g, sag):
20 Lx = (g.nx - 1) * g.hx
21 return sag * (1 - (2 * g.X / Lx) ** 2) + 0.3 * sag * np.sin(3 * g.Y / Lx)
24def random_fields(pl, amp, seed):
25 r = np.random.default_rng(seed)
26 return [np.where(pl.grid.mask, a * (r.random(pl.grid.shape) - 0.5), 0.0) for a in (0.02 * amp, 0.02 * amp, amp)]
30def test_stiffness_has_exactly_the_rigid_modes_in_its_null_space(case):
31 if case == "disc-hole":
32 g = make_grid(0.16, 0.16, 0.01, lambda x, y, r, t: (r <= 0.08) & (r >= 0.015))
33 pl, fields = Plate.uniform(g, STEEL, 1e-3), 1
34 else:
35 g = make_grid(0.12, 0.09, 0.015)
36 pl, fields = Plate.uniform(g, STEEL, 1e-3), 3
37 if case == "tapered":
38 pl = pl.with_thickness(1e-3 * (1 + 0.5 * (g.X + 0.06) / 0.12))
39 sys = assemble(pl, fields)
40 zero = np.zeros(g.shape)
41 one = np.ones(g.shape)
42 null = [(zero, zero, one), (zero, zero, g.X), (zero, zero, g.Y)]
43 if fields == 3:
44 null += [(one, zero, zero), (zero, one, zero), (-g.Y, g.X, zero)]
45 scale = abs(sys.K[0, 0])
46 for u, v, w in null:
47 x = sys.from_grid(u, v, w)
48 assert np.linalg.norm(sys.K @ x) < 1e-9 * scale * np.linalg.norm(x)
49 x = sys.from_grid(zero, zero, np.sin(0.9 * np.arange(g.nx))[None, :] * np.cos(1.1 * np.arange(g.ny))[:, None])
50 assert np.abs(sys.K @ x).max() > 1e-3
51 assert abs(sys.K - sys.K.T).max() < 1e-9 * scale
55def test_assembled_K_equals_the_linearized_kernel_on_a_flat_plate(fields):
56 pl = Plate.uniform(make_grid(0.2, 0.15, 0.02), BRONZE, 1e-3)
57 assert pl.is_flat
58 assert relative_error_against_kernel(pl, fields) < 1e-6
61def test_assembled_K_equals_the_linearized_kernel_on_a_bowed_tapered_plate_with_a_hole():
62 g = make_grid(0.2, 0.2, 0.0125, lambda x, y, r, t: (r <= 0.1) & (r >= 0.01))
63 pl = Plate.uniform(g, BRONZE, 1e-3, profile=3e-3 * np.sin(20 * g.X) * np.cos(15 * g.Y))
64 pl = pl.with_thickness(1e-3 * (1 + 0.4 * np.abs(g.X) / 0.1))
65 assert not pl.is_flat
66 assert assemble(pl).fields == 3
67 assert relative_error_against_kernel(pl, 3) < 1e-6
70def relative_error_against_kernel(pl, fields, seed=3):
71 sys = assemble(pl, fields)
72 x = np.random.default_rng(seed).random(sys.n) - 0.5
73 y = sys.K @ x
74 eps = 1e-9
75 u, v, w = sys.to_grid(x * eps)
76 S = static_data(pl, None, 1e-6, 0, 0, 0, "cpu", torch.float64)
77 au, av, aw = (a.numpy() for a in accelerations(S, T(u), T(v), T(w)))
78 m = pl.node_mass
79 yref = sys.from_grid(-au * m / eps, -av * m / eps, -aw * m / eps)
80 return np.linalg.norm(y - yref) / np.linalg.norm(yref)
84def test_forces_are_the_gradient_of_the_energy(case):
85 g = make_grid(0.1, 0.08, 0.01)
86 pl = Plate.uniform(g, STEEL, 1e-3, profile=None if case == "flat" else bowed(g, 3e-3))
87 if case == "bowed-tapered":
88 pl = pl.with_thickness(1e-3 * (1 + 0.4 * (g.X + 0.05) / 0.1))
89 f = random_fields(pl, 2.5e-3, 7)
90 d = random_fields(pl, 1.0, 11)
91 S = static_data(pl, None, 1e-6, 0, 0, 0, "cpu", torch.float64)
92 acc = [a.numpy() for a in accelerations(S, *[T(a) for a in f])]
93 h = 1e-7
94 dU = (potential_energy(pl, *[a + h * b for a, b in zip(f, d)]) - potential_energy(pl, *[a - h * b for a, b in zip(f, d)])) / (2 * h)
95 fd = -sum((pl.node_mass * a * b).sum() for a, b in zip(acc, d))
96 assert abs(dU - fd) < 1e-6 * abs(dU)
99def test_the_loss_terms_only_remove_energy():
100 pl = Plate.uniform(make_grid(0.1, 0.08, 0.01), STEEL, 1e-3)
101 S = static_data(pl, None, 1e-6, 1.0, 5e-4, 20.0, "cpu", torch.float64)
102 for seed in (1, 2, 3):
103 vel = random_fields(pl, 1.0, seed)
104 acc = [a.numpy() for a in damping(S, *[T(a) for a in vel])]
105 power = sum((pl.node_mass * a * b).sum() for a, b in zip(acc, vel))
106 assert power < 0
109def test_doubling_the_thickness_doubles_the_bending_frequencies_only():
110 g = make_grid(0.1, 0.1, 0.01)
111 w1, w2 = max_angular_frequency(g, STEEL, 1e-3), max_angular_frequency(g, STEEL, 2e-3)
112 # The largest frequency is set by the in-plane waves, independent of h.
113 assert abs(w1 / w2 - 1) < 1e-6
114 # It lies near 2 c_p / h for the grid-scale in-plane mode.
115 assert 0.8 < w1 * g.hx / (2 * STEEL.in_plane_speed) < 1.5
116 s1, s2 = assemble(Plate.uniform(g, STEEL, 1e-3)), assemble(Plate.uniform(g, STEEL, 2e-3))
117 c = s1.n // 2
118 assert abs(s2.K[c, c] / s1.K[c, c] - 8) < 1e-9
119 assert abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12
122def test_the_band_is_narrow_whichever_axis_is_longer():
123 wide = assemble(Plate.uniform(make_grid(0.4, 0.1, 0.01), BRONZE, 1e-3))
124 tall = assemble(Plate.uniform(make_grid(0.1, 0.4, 0.01), BRONZE, 1e-3))
125 assert wide.bandwidth == tall.bandwidth <= 2 * 11 + 2