1"""The room problem.
3A point source at x0 inside a rounded cube (the superellipsoid
4|x|^p + |y|^p + |z|^p = a^p) with impedance walls. With time dependence
5exp(-i omega t) the pressure amplitude u satisfies
7 (Delta + k^2) u = -delta(x - x0) in the room,
8 du/dn + i k lambda u = 0 on the walls,
10with n the outward normal and lambda = -beta, where beta > 0 is the wall
11admittance. We write u = G + v with G the free-space field of the source, and
12find v from the right-preconditioned combined field integral equation of
13fmm3dbie (helm_rpcomb_imp_solver, interior mode).
15This module has no dependence on either backend.
16"""
18import time
20import numpy as np
22from .common import green
24NAME = "room"
25IFINOUT = 0 # interior problem
27# Named cases. Both backends discretize the surface identically, so the same
28# case gives the same unknowns on both.
29CASES = {
30 "small": dict(k=10.0, beta=0.2, m=4, norder=6), # 3.2 wavelengths, 4,704 points
31 "medium": dict(k=10.0, beta=0.2, m=6, norder=8), # 3.2 wavelengths, 17,496 points
32 "large": dict(k=40.0, beta=0.05, m=10, norder=8), # 12.7 wavelengths, 48,600 points
33}
34A, P = 1.0, 8 # room half-width and superellipsoid exponent
35SOURCE = np.array([-0.35, 0.2, 0.1]) # the point source
36EXTERIOR_SOURCE = np.array([1.7, 0.6, -0.5]) # for the self-test (outside the room)
37ALPHA = 1.0
38NGRID = 101 # default grid size of the field slices
41def add_arguments(ap):
42 ap.add_argument("--k", type=float, help="override the case's wavenumber")
43 ap.add_argument("--beta", type=float, help="override the case's wall admittance")
44 ap.add_argument("--m", type=int, help="override the patches per face edge")
45 ap.add_argument("--norder", type=int, help="override the order per patch")
48def make_case(name, args):
49 case = dict(CASES[name])
50 for key in ["k", "beta", "m", "norder"]:
51 if getattr(args, key, None) is not None:
52 case[key] = getattr(args, key)
53 case["lam"] = -case["beta"]
54 return case
57def results_name(case_name, backend):
58 return f"{case_name}-{backend}"
61def describe(case):
62 return (f"k = {case['k']}, beta = {case['beta']}, "
63 f"room width {2 * A * case['k'] / (2 * np.pi):.2f} wavelengths")
66def gpu_surface(case, args):
67 from fmm3dbie_gpu import surface
68 return surface.rounded_cube(a=A, p=P, m=case["m"], norder=case["norder"],
69 exact=not args.interp_geometry)
72def fmm3dbie_surface(case):
73 import rounded_cube as rc
74 return rc.rounded_cube(a=A, p=P, m=case["m"], norder=case["norder"])
77def source(case):
78 return SOURCE
81def selftest_sources(case):
82 """A source outside the room, (positions (3, n), strengths (n,)); its field is the
83 exact solution inside."""
84 return EXTERIOR_SOURCE[:, None], np.ones(1)
87def selftest_points(case, n=50):
88 return np.random.default_rng(0).uniform(-0.8 * A, 0.8 * A, size=(3, n))
91def inside(x, a=A, p=P, margin=0.0):
92 return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p
95def slice_points(ngrid):
96 """Points on three planes through the source; returns grid, {name: (pts, mask, shape)}."""
97 g = np.linspace(-A, A, ngrid)
98 out = {}
99 for name, (i, j, fixed) in {"xy": (0, 1, 2), "xz": (0, 2, 1), "yz": (1, 2, 0)}.items():
100 Ag, Bg = np.meshgrid(g, g, indexing="xy")
101 pts = np.zeros((3, Ag.size))
102 pts[i], pts[j], pts[fixed] = Ag.ravel(), Bg.ravel(), SOURCE[fixed]
103 out[name] = (pts, inside(pts, margin=1e-3), Ag.shape)
104 return g, out
107def energy_balance(k, lam, v_at_source, wts, u_wall):
108 """Power radiated by the source, k/(4 pi) + Im v(x0), and power absorbed by the
109 walls, -k Re(lambda) int |u|^2 dS. They must agree."""
110 p_src = k / (4 * np.pi) + v_at_source.imag
111 p_wall = -k * np.real(lam) * np.sum(wts * np.abs(u_wall) ** 2)
112 return p_src, p_wall
115def report(be, state, case, args, summary, arrays):
116 """Energy balance and field slices of the room solution."""
117 k, lam = case["k"], case["lam"]
118 xs, wts = np.asarray(be.surf.srcvals)[0:3], np.asarray(be.surf.wts)
119 v0 = complex(be.eval_points(state, SOURCE[:, None])[0])
120 u_wall = green(xs, SOURCE, k) + be.eval_wall(state)
121 p_src, p_wall = energy_balance(k, lam, v0, wts, u_wall)
122 mismatch = abs(p_src - p_wall) / p_src
123 print(f" energy balance: source {p_src:.8f}, walls {p_wall:.8f}, relative mismatch {mismatch:.2e}")
124 summary.update(p_src=p_src, p_wall=p_wall, energy_mismatch=mismatch)
125 arrays.update(x0=SOURCE, k=k, lam=lam, a=A, p=P, xs=xs, wts=wts, u_wall=u_wall)
126 if args.ngrid > 0:
127 g, sl = slice_points(args.ngrid)
128 t0 = time.time()
129 for name, (pts, mask, shape) in sl.items():
130 u = np.full(mask.size, np.nan + 1j * np.nan)
131 u[mask] = green(pts[:, mask], SOURCE, k) + be.eval_points(state, pts[:, mask])
132 arrays[f"u_{name}"] = u.reshape(shape)
133 be.sync()
134 summary["slices_s"] = time.time() - t0
135 arrays["grid"] = g
136 print(f" field on three {args.ngrid} x {args.ngrid} slices: {summary['slices_s']:.2f} s")