1"""The cymbal problem: a reciprocity solve for a listener near a rigid cymbal.
3A point source at the listener's position L radiates in free space; the only
4obstacle is a stationary, sound-hard cymbal. With time dependence
5exp(-i omega t) the pressure amplitude u satisfies
7 (Delta + k^2) u = -delta(x - L) outside the cymbal,
8 du/dn = 0 on the cymbal,
10and the Sommerfeld radiation condition. By reciprocity, u(x) is also the
11pressure at the listener due to a unit point source at x, so one solve gives
12the listener's response to a source anywhere; we evaluate it on a 3D grid in a
13bounding box, and on three planes through the box. The box only selects where
14the field is evaluated: the integral representation satisfies the radiation
15condition, so the problem needs no truncation.
17The cymbal is a thin closed body (fmm3dbie_gpu/cymbal.py): fmm3dbie's integral
18equations need a closed surface. We write u = G + v, with G the free-space
19field of the source, and find v from fmm3dbie's right-preconditioned combined
20field equation for the impedance problem in exterior mode with lambda = 0
21(helm_rpcomb_imp_solver, ifinout = 1), the same solver as the room problem.
22Lengths are in meters and the speed of sound is 343 m/s.
24This module has no dependence on either backend.
25"""
27import time
29import numpy as np
31from fmm3dbie_gpu.cymbal import CymbalShape, cymbal, cymbal_srcvals
33from .common import green, impedance_data
35NAME = "cymbal"
36IFINOUT = 1 # exterior problem
37ALPHA = 1.0
38C_SOUND = 343.0
40# h is the target patch size in meters; the mesh is also refined toward the rim
41# and across the bell, independently of the frequency.
42CASES = {
43 "small": dict(freq=1000.0, h=0.06, norder=6),
44 "medium": dict(freq=3000.0, h=0.04, norder=8),
45 "large": dict(freq=8000.0, h=0.03, norder=8),
46}
47RADIUS = 0.2 # a 40 cm cymbal
48T_RIM = 2 * RADIUS / 40 # rim thickness diameter/40; twice that at the center
49LISTENER = np.array([0.25, 0.0, 0.45]) # above the cymbal and off its axis
50BOX = np.array([[-0.5, 0.5], [-0.5, 0.5], [-0.35, 0.65]]) # evaluation region (1 m cube)
51RECIPROCITY_POINT = np.array([-0.15, 0.05, -0.12]) # below the cymbal
52SLICE_BELOW = -0.1 # height of the horizontal slice, below the cymbal
53NGRID = 201 # default grid size of the slices, 5 mm spacing
56def add_arguments(ap):
57 ap.add_argument("--freq", type=float, help="override the case's frequency (Hz)")
58 ap.add_argument("--h", type=float, help="override the case's target patch size (m)")
59 ap.add_argument("--norder", type=int, help="override the order per patch")
60 ap.add_argument("--t-rim", type=float, help=f"rim thickness in m (default {T_RIM})")
61 ap.add_argument("--listener", type=lambda s: [float(v) for v in s.split(",")],
62 help="listener (= source) position x,y,z in m")
63 ap.add_argument("--reciprocity", action="store_true",
64 help="also solve with the source at a point below the cymbal and check "
65 "that the field there due to the listener equals the field at the "
66 "listener due to it (a second solve)")
67 ap.add_argument("--n3d", type=int, default=-1,
68 help="grid points per axis of the 3D grid in the box (default: a quarter "
69 "wavelength spacing, 41 to 121 points; 0 to skip)")
72def make_case(name, args):
73 case = dict(CASES[name])
74 for key in ["freq", "h", "norder", "t_rim"]:
75 if getattr(args, key, None) is not None:
76 case[key] = getattr(args, key)
77 case.setdefault("t_rim", T_RIM)
78 case["listener"] = list(args.listener) if getattr(args, "listener", None) else list(LISTENER)
79 case["k"] = 2 * np.pi * case["freq"] / C_SOUND
80 case["lam"] = 0.0
81 return case
84def results_name(case_name, backend):
85 return f"cymbal-{case_name}-{backend}"
88def shape(case):
89 return CymbalShape(radius=RADIUS, t_rim=case["t_rim"])
92def describe(case):
93 wl = C_SOUND / case["freq"]
94 return (f"f = {case['freq']:g} Hz, k = {case['k']:.2f}/m, cymbal diameter "
95 f"{2 * RADIUS / wl:.2f} wavelengths, rim thickness {1000 * case['t_rim']:.1f} mm "
96 f"(diameter/{2 * RADIUS / case['t_rim']:.0f})")
99def gpu_surface(case, args):
100 return cymbal(shape(case), case["h"], case["norder"])
103def fmm3dbie_surface(case):
104 import fmm3dbie_ctypes as fb
105 srcvals, _ = cymbal_srcvals(shape(case), case["h"], case["norder"])
106 return fb.Surface(case["norder"], 11, srcvals)
109def source(case):
110 return np.array(case["listener"])
113def selftest_sources(case, r_max=0.12, spacing=0.003):
114 """Point sources inside the cymbal, (positions (3, n), strengths (n,)); their field
115 is the exact solution outside. A single source would be at most T/2 = 1 cm from
116 the surface, and its data would vary on that scale, which the main problem's
117 data does not. Instead we spread sources 3 mm apart over the midsurface out to
118 r = 12 cm, where they are at least 7 mm from the surface, with strengths that
119 taper smoothly to zero at r_max. The data is then smooth on the scale of the
120 patches, apart from a ripple of relative size about exp(-2 pi 7 / 3) ~ 4e-7."""
121 sh = shape(case)
122 g = np.arange(-r_max, r_max + spacing / 2, spacing)
123 X, Y = np.meshgrid(g, g, indexing="ij")
124 rho = (X * X + Y * Y).ravel()
125 m = rho < r_max ** 2
126 x0s = np.stack([X.ravel()[m], Y.ravel()[m], sh.F(rho[m])])
127 ws = (1 - rho[m] / r_max ** 2) ** 3 * spacing ** 2 / (np.pi * r_max ** 2 / 4) # sum ~ 1
128 return x0s, ws
131def selftest_points(case, n=50):
132 """Points around the cymbal, at least 5 mm from its surface."""
133 sh = shape(case)
134 pts = np.random.default_rng(0).uniform([-0.3, -0.3, -0.15], [0.3, 0.3, 0.2], size=(4 * n, 3)).T
135 return pts[:, ~sh.inside(pts, margin=0.005)][:, :n]
138def far_field_power(k, x0, xs, nrm, wts, sigma, phi1, alpha=ALPHA, xp=np):
139 """Power radiated to infinity, k int |u_inf|^2 dOmega, for u = G(., x0) + S_k sigma
140 + i alpha D_k phi1, with u ~ u_inf exp(ikr)/r. The far field pattern is computed
141 with the smooth rule at the surface nodes, and the sphere integral with Gauss-
142 Legendre in cos(theta) and the trapezoid rule in phi, exact for the band limit.
143 xp is numpy or cupy."""
144 rmax = max(np.linalg.norm(xs, axis=0).max(), np.linalg.norm(x0))
145 L = int(np.ceil(k * rmax)) + 20
146 ct, wct = np.polynomial.legendre.leggauss(L + 1)
147 nphi = 2 * L + 2
148 ph = 2 * np.pi * np.arange(nphi) / nphi
149 st = np.sqrt(1 - ct ** 2)
150 dirs = np.stack([np.outer(st, np.cos(ph)), np.outer(st, np.sin(ph)),
151 np.outer(ct, np.ones(nphi))]).reshape(3, -1)
152 wdir = np.outer(wct, np.full(nphi, 2 * np.pi / nphi)).ravel()
153 xs, nrm, dirs = xp.asarray(xs), xp.asarray(nrm), xp.asarray(dirs)
154 qs, qd = xp.asarray(wts * sigma), xp.asarray(1j * alpha * wts * phi1)
155 uinf = xp.exp(-1j * k * (dirs.T @ xp.asarray(x0)))
156 step = 64 if xp is np else 1024
157 for i in range(0, dirs.shape[1], step):
158 d = dirs[:, i:i + step]
159 e = xp.exp(-1j * k * (d.T @ xs))
160 uinf[i:i + step] += e @ qs + ((-1j * k * (d.T @ nrm)) * e) @ qd
161 uinf = uinf.get() if xp is not np else uinf
162 uinf /= 4 * np.pi
163 return k * np.sum(wdir * np.abs(uinf) ** 2)
166def box_grid(n):
167 axes = [np.linspace(lo, hi, n) for lo, hi in BOX]
168 X, Y, Z = np.meshgrid(*axes, indexing="ij")
169 return axes, np.stack([X.ravel(), Y.ravel(), Z.ravel()])
172def slice_points(ngrid, x0):
173 """Three planes through the box: xz through the listener, yz through the cymbal's
174 axis, and xy below the cymbal. Returns {name: (pts, (a-axis, b-axis), fixed)}."""
175 out = {}
176 for name, (i, j, fixed, val) in {"xz": (0, 2, 1, x0[1]), "yz": (1, 2, 0, 0.0),
177 "xy": (0, 1, 2, SLICE_BELOW)}.items():
178 ga, gb = np.linspace(*BOX[i], ngrid), np.linspace(*BOX[j], ngrid)
179 Ag, Bg = np.meshgrid(ga, gb, indexing="xy")
180 pts = np.zeros((3, Ag.size))
181 pts[i], pts[j], pts[fixed] = Ag.ravel(), Bg.ravel(), val
182 out[name] = (pts, (ga, gb), (fixed, val))
183 return out
186def eval_total(be, state, sh, x0, k, pts, spacing):
187 """u = G + v at the points, NaN inside the cymbal and at the source."""
188 mask = ~sh.inside(pts, margin=1e-3) & (np.linalg.norm(pts - x0[:, None], axis=0) > 1e-3 * spacing)
189 u = np.full(mask.size, np.nan + 1j * np.nan)
190 u[mask] = green(pts[:, mask], x0, k) + be.eval_points(state, pts[:, mask])
191 return u
194def report(be, state, case, args, summary, arrays):
195 k, x0, sh = case["k"], source(case), shape(case)
196 srcvals, wts = np.asarray(be.surf.srcvals), np.asarray(be.surf.wts)
197 xs, nrm = srcvals[0:3], srcvals[9:12]
198 v0 = complex(be.eval_points(state, x0[:, None])[0])
199 u_wall = green(xs, x0, k) + be.eval_wall(state)
201 t0 = time.time()
202 sigma, phi1 = (np.asarray(be.to_numpy(q)) for q in state)
203 p_src = k / (4 * np.pi) + v0.imag
204 p_far = far_field_power(k, x0, xs, nrm, wts, sigma, phi1, xp=be.xp)
205 mismatch = abs(p_src - p_far) / p_src
206 print(f" energy balance: source {p_src:.8f}, far field {p_far:.8f}, relative mismatch "
207 f"{mismatch:.2e} ({time.time() - t0:.1f} s)")
208 summary.update(p_src=p_src, p_far=p_far, energy_mismatch=mismatch, v_listener=[v0.real, v0.imag])
210 if args.reciprocity:
211 xb = RECIPROCITY_POINT
212 (state_b, niter_b, _), t = be.timed_solve(-impedance_data(xs, nrm, xb, k, 0.0))
213 v_lb = complex(be.eval_points(state, xb[:, None])[0])
214 v_bl = complex(be.eval_points(state_b, x0[:, None])[0])
215 rdiff = abs(v_lb - v_bl) / abs(v_lb)
216 print(f" reciprocity: scattered field at {xb} due to the listener {v_lb:.6e}, at the "
217 f"listener due to {xb} {v_bl:.6e}, relative difference {rdiff:.2e} "
218 f"({niter_b} iterations, {t:.2f} s)")
219 summary.update(reciprocity_diff=rdiff, reciprocity_solve_s=t)
221 arrays.update(x0=x0, k=k, freq=case["freq"], xs=xs, wts=wts, u_wall=u_wall, box=BOX,
222 outline=sh.cross_section(), radius=RADIUS)
223 wl = 2 * np.pi / k
224 if args.ngrid > 0:
225 t0 = time.time()
226 for name, (pts, (ga, gb), (fixed, val)) in slice_points(args.ngrid, x0).items():
227 u = eval_total(be, state, sh, x0, k, pts, (ga[1] - ga[0]))
228 arrays[f"u_{name}"] = u.reshape(gb.size, ga.size)
229 arrays[f"axes_{name}"] = np.stack([ga, gb])
230 arrays[f"plane_{name}"] = np.array([fixed, val])
231 be.sync()
232 summary["slices_s"] = time.time() - t0
233 print(f" field on three {args.ngrid} x {args.ngrid} slices: {summary['slices_s']:.2f} s")
234 n3d = args.n3d
235 if n3d < 0:
236 n3d = int(np.clip(np.ceil((BOX[0, 1] - BOX[0, 0]) / (wl / 4)) + 1, 41, 121))
237 if n3d > 0:
238 t0 = time.time()
239 axes, pts = box_grid(n3d)
240 u = eval_total(be, state, sh, x0, k, pts, axes[0][1] - axes[0][0])
241 arrays["u_grid"] = u.reshape(n3d, n3d, n3d)
242 arrays["grid_axes"] = np.stack(axes)
243 be.sync()
244 summary.update(grid_s=time.time() - t0, n3d=n3d)
245 print(f" field on a {n3d}^3 grid in the box: {summary['grid_s']:.2f} s")