Revert "Optional slope-corrected Rayleigh quadrature (listener.quadrature: slope) against the spurious high-frequency radiation of the element-mean sum"
This reverts commit 401ef929b72c6d905a50ef3c892e389cd8cd5d77.
5 changed files+3−78
docs/model.mdmodified+0−8View file
@@ -210,14 +210,6 @@ $$
210210
211211 where $\bar a_e$ is the mean over the element's three nodes of the *normal acceleration* $a_i = n_i \cdot (\ddot u_i, \ddot v_i, \ddot w_i)$, $n_i$ being the unit normal of the rest surface at node $i$ (the area-weighted mean of the adjacent elements' constant slopes, pointing to the $+z$ side), and $r_e$ the distance from the element centroid, placed on the curved surface at height $z_0$, to the listener. The acceleration is taken from the integrator directly: for the leapfrog scheme $(v^{n+1} - v^n) / \Delta t$ is the exact central second difference of $q$ at step $n$, so no differentiation of stored data is involved. This is exact for a rigid piston in an infinite baffle and an approximation here: the surface is shallow rather than flat, and the air does not load the shell (radiation damping is part of what $\sigma_0$ and $\sigma_1$ lump). At the bell, where the slope reaches 0.39, the normal component differs from $\ddot w$ by several percent, which is why we use it rather than $\ddot w$.
212212
213-The element mean $\bar a_e$ integrates exactly a field that is linear over the element, and that is not accurate enough here. Below the coincidence frequency, where bending waves travel slower than sound (about 24 kHz for a 0.5 mm steel plate, lower for thicker plates), the true radiation is a small residue of a near-complete cancellation over the surface, and only a few modes per band radiate strongly. The quadrature error leaves a residue of its own for every mode, with an effectively random weight, which in a cascade adds a dense broadband noise to the listener signal. On a flat 0.5 mm steel plate with 6 mm elements, this error is about 10 dB above the true radiated power in every band, and it falls at second order with the element size. With `listener.quadrature: slope` the element source is instead
214-
215-$$
216-\bar a_e = \frac13 \sum_{i \in e} a_i + \frac18 \sum_{i \in e} \nabla \ddot w_i \cdot (x_c - x_i), \qquad \nabla \ddot w_i = (-\ddot\theta_{y,i},\, \ddot\theta_{x,i}),
217-$$
218-
219-which uses the rotational accelerations the integrator already has and is exact for quadratic fields. On the same plate it lowers the error by 15 to 18 dB above 2 kHz. For a curved shell the slope of the normal acceleration differs slightly from that of $\ddot w$; we use the latter. The default remains `mean`, because the stored acceleration record holds no rotations, so a reader can reproduce the listener signal from the record only for `mean`, and because the web companion implements only `mean`.
220-
221213 A cymbal has no baffle. Its back face radiates with the opposite sign, and that sound reaches the listener by diffracting around the rim; at low frequency the two contributions nearly cancel, which is why an unbaffled plate is a weak low-frequency source. We estimate the back contribution with the same sum but with $r_e$ replaced by the shortest path over the rim, $r_e^{back} = \min_j \big( |x_e - x_j^{rim}| + |x_j^{rim} - x_L| \big)$ over 720 rim points, weighted by `back_diffraction_gain`, and subtract it: $p_{unbaffled} = p_{front} - g\, p_{back}$. This is a Kirchhoff-style edge-diffraction estimate, not a boundary-element solution, and it ignores shadowing of the front by the shell itself. Both signals are returned; the baffled one is a cleaner measure of the shell motion, the unbaffled one sounds more like a cymbal in a room.
222214
223215 Retarded times are handled at sub-step resolution: the delay of element $e$ in sub-steps, $d_e = r_e / (c \Delta t)$, is split into an integer part and a fraction, and each element's contribution is added to two adjacent slots of a ring buffer with linear interpolation. After each block the host reads the slots of the finished sub-steps, clears them and decimates them to the audio rate with the kernel of section 9, the same kernel that filters the acceleration record. Since delay and filtering are both linear and time-invariant, the stored listener signal follows from the stored accelerations up to the linear-interpolation error at the sub-step scale, which is negligible; the test suite rebuilds the baffled signal from the record to $1.5 \times 10^{-4}$ relative rms.
src/cymsim/config.pymodified+0−5View file
@@ -125,9 +125,6 @@ class Listener:
125125 c_air: float = 343.0
126126 rho_air: float = 1.2
127127 back_diffraction_gain: float = 1.0 # weight of the diffracted back-side term (0 = baffled only)
128- quadrature: str = "mean" # per-element Rayleigh source: "mean" of the nodal accelerations, or "slope"
129- # (adds the nodal slopes; exact for quadratic fields, far less spurious
130- # radiation below coincidence; see docs/model.md, section 7)
131128
132129
133130 # ------------------------------------------------------------------------------- numerics
@@ -178,8 +175,6 @@ class SimulationConfig:
178175 raise ValueError("numerics.cfl_safety must lie in (0, 1]")
179176 if self.mesh.element_size <= 0 or self.mesh.element_size > p.radius / 2:
180177 raise ValueError("mesh.element_size must be positive and well below the radius")
181- if self.listener.quadrature not in ("mean", "slope"):
182- raise ValueError(f"listener.quadrature must be 'mean' or 'slope', got {self.listener.quadrature!r}")
183178 if self.output.displacement_rate < 0:
184179 raise ValueError("output.displacement_rate must be non-negative")
185180 for i, h in enumerate(self.hits):
src/cymsim/radiation.pymodified+0−22View file
@@ -8,12 +8,6 @@ contribution reaches the listener by diffracting around the rim. We approximate
88 Rayleigh sum with the path length replaced by the shortest route over the rim (a Kirchhoff-type
99 edge-diffraction estimate, not a boundary-element solution) and subtract it from the front term.
1010
11-Each element's source strength is the mean of its nodal normal accelerations, or, with
12-``listener.quadrature = "slope"``, that mean corrected with the nodal slopes of the transverse
13-acceleration (the rotational accelerations), which makes the element integral exact for quadratic
14-fields. Below the coincidence frequency the true radiation is a small residue of a near-complete
15-cancellation over the surface, and the plain mean leaves a spurious residue of its own.
16-
1711 Retarded times are handled at sub-step resolution with linear interpolation into two delay
1812 buffers. The buffers are rings that hold a little more than two blocks plus the longest delay;
1913 after each block the host reads the finished sub-steps, clears their slots and decimates them to
@@ -31,20 +25,6 @@ from .geometry import CymbalProfile
3125 from .mesh import Mesh
3226
3327
34-def element_sources(a_node, slope_x, slope_y, tri, xc):
35- """Per-element source strength (mean acceleration over the element) for the Rayleigh sum.
36-
37- With ``slope_x`` = None this is the mean of the three nodal values, exact for a field that is
38- linear on the element. Otherwise the nodal slopes are added,
39- (1/A) int_e a dA ~ (1/3) sum_i a_i + (1/8) sum_i grad a_i . (x_c - x_i),
40- which is exact for quadratic fields. ``xc``: (E, 3, 2) offsets x_c - x_i on the projected domain.
41- """
42- acc_e = a_node[tri].mean(1)
43- if slope_x is None:
44- return acc_e
45- return acc_e + 0.125 * (slope_x[tri] * xc[..., 0] + slope_y[tri] * xc[..., 1]).sum(1)
46-
47-
4828 class RadiationModel:
4929 def __init__(self, mesh: Mesh, profile: CymbalProfile, z0_node: np.ndarray, listener: Listener,
5030 dt: float, steps_per_block: int, dtype, n_rim: int = 720):
@@ -76,8 +56,6 @@ class RadiationModel:
7656 self.coef_f, self.coef_b = f(coef_f), f(coef_b)
7757 self.frac_f, self.frac_b = f(df - np.floor(df)), f(db - np.floor(db))
7858 self.i0f, self.i0b = jnp.asarray(i0f), jnp.asarray(i0b)
79- self.slope = listener.quadrature == "slope"
80- self.xc = f(mesh.centroid[:, None, :] - mesh.xy[mesh.tri]) # (E, 3, 2)
8159 self.dtype = dtype
8260
8361 def initial_buffers(self):
src/cymsim/simulate.pymodified+3−7View file
@@ -36,7 +36,7 @@ from .decimation import DECIMATION_SPAN, StreamingDecimator, decimation_taps
3636 from .exciter import HitSchedule, build_patch
3737 from .geometry import CymbalProfile
3838 from .mesh import Mesh, disk_mesh
39-from .radiation import RadiationModel, element_sources
39+from .radiation import RadiationModel
4040 from .shell import N_DOF, ShellModel
4141
4242 log = logging.getLogger("cymsim")
@@ -341,12 +341,8 @@ class Simulator:
341341 q1 = q + dt * v1
342342 zs1, vs1 = sched.advance(zs, vs, F, dt)
343343 # normal acceleration at step n: central second difference of q, exact for leapfrog
344- acc = (v1 - v) / dt
345- a_node = (acc[:, :3] * normal).sum(1)
346- if rad.slope: # slopes of the transverse acceleration: grad w = (-theta_y, theta_x)
347- acc_e = element_sources(a_node, -acc[:, 4], acc[:, 3], tri, rad.xc)
348- else:
349- acc_e = element_sources(a_node, None, None, tri, rad.xc)
344+ a_node = (((v1 - v) / dt)[:, :3] * normal).sum(1)
345+ acc_e = a_node[tri].mean(1)
350346 bf, bb = rad.scatter(n, acc_e, bf, bb)
351347 if n_probes:
352348 pbuf = pbuf.at[n % spb].set((a_node[p_nodes] * p_w).sum(1)) # blocks start at multiples of spb
tests/test_radiation.pydeleted+0−36View file
@@ -1,36 +0,0 @@
1-import copy
2-
3-import numpy as np
4-
5-from cymsim.mesh import disk_mesh
6-from cymsim.radiation import element_sources
7-from cymsim.simulate import Simulator
8-
9-
10-def test_slope_corrected_sources_are_exact_for_quadratics():
11- m = disk_mesh(0.2, 0.02)
12- x, y = m.xy[:, 0], m.xy[:, 1]
13- a = 1.0 + 2 * x - 3 * y + 5 * x * x - 4 * x * y + 7 * y * y
14- gx, gy = 2 + 10 * x - 4 * y, -3 - 4 * x + 14 * y
15- xc = m.centroid[:, None, :] - m.xy[m.tri]
16- # exact element means of the quadratic: value at the centroid plus the second-moment term
17- p = m.xy[m.tri] - m.centroid[:, None, :]
18- xx, xy_, yy = [(p[:, :, i] * p[:, :, j]).sum(1) / 12 for i, j in ((0, 0), (0, 1), (1, 1))]
19- cx, cy = m.centroid[:, 0], m.centroid[:, 1]
20- exact = 1 + 2 * cx - 3 * cy + 5 * (cx * cx + xx) - 4 * (cx * cy + xy_) + 7 * (cy * cy + yy)
21- corrected = np.asarray(element_sources(a, gx, gy, m.tri, xc))
22- plain = np.asarray(element_sources(a, None, None, m.tri, xc))
23- assert np.abs(corrected - exact).max() < 1e-12
24- assert np.abs(plain - exact).max() > 1e-5 # the plain mean is not exact
25-
26-
27-def test_slope_quadrature_run(coarse_config):
28- cfg = copy.deepcopy(coarse_config)
29- cfg.listener.quadrature = "slope"
30- res = Simulator(cfg).run()
31- ref = Simulator(coarse_config).run()
32- p, q = res.pressure_rayleigh, ref.pressure_rayleigh
33- assert np.isfinite(p).all()
34- # same motion, slightly different sum: the two signals agree at low frequency and differ a little
35- rel = np.sqrt(np.mean((p - q) ** 2) / np.mean(q ** 2))
36- assert 1e-6 < rel < 0.5