Add the cymbal problem: a listener near a rigid cymbal in free space
A reciprocity solve with the exterior mode of the impedance solver (lambda = 0)
around a thin closed cymbal with a raised bell and a tapered rim. The problems
move to problems/ (room, cymbal) and run.py selects one with --problem. Both
backends run it; the self-test uses sources spread inside the cymbal, and the
energy balance compares the source power with the power radiated to infinity.
14 changed files+940−205
.gitignoremodified+1−0View file
@@ -1,3 +1,4 @@
11 results/
22 __pycache__/
33 *.pyc
4+fort.*
README.mdmodified+43−2View file
@@ -1,11 +1,13 @@
11 # fmm3dbie on the GPU: a proof of concept
22
3-fmm3dbie solves boundary integral equations on surfaces in three dimensions with high-order quadrature and the fast multipole method. On a CPU, most of its run time for the Helmholtz equation goes into computing the near-field quadrature corrections, and the rest into FMM-accelerated matrix-vector products inside GMRES. Here we present a GPU implementation of one fmm3dbie solver, the interior Helmholtz impedance problem (`helm_rpcomb_imp_solver`). It uses the same integral equation and the same surface discretization, together with a showcase problem that runs on both codes so that their speed and accuracy can be compared directly.
3+fmm3dbie solves boundary integral equations on surfaces in three dimensions with high-order quadrature and the fast multipole method. On a CPU, most of its run time for the Helmholtz equation goes into computing the near-field quadrature corrections, and the rest into FMM-accelerated matrix-vector products inside GMRES. Here we present a GPU implementation of one fmm3dbie solver, the Helmholtz impedance problem (`helm_rpcomb_imp_solver`), in its interior and exterior modes. It uses the same integral equation and the same surface discretization, together with a showcase problem that runs on both codes so that their speed and accuracy can be compared directly.
44
55 The showcase is a point source at a fixed frequency in a closed room with partially absorbing walls. The room is a slightly rounded cube, the superellipsoid |x|⁸ + |y|⁸ + |z|⁸ = 1. With time dependence e^{-iωt}, the pressure amplitude satisfies Δu + k²u = −δ(x − x₀) inside the room and ∂u/∂n + ikλu = 0 on the walls, with outward normal n and λ = −β. Here β > 0 is the normalized wall admittance, and smaller β gives more reflective walls.
66
77 
88
9+A second showcase, a listener near a rigid cymbal in free space, is described in [The cymbal problem](#the-cymbal-problem). It uses the exterior mode of the same solver, and both backends run it.
10+
911 A single precision port of the same solver to WebGPU, which runs in a browser without CUDA, is in [fmm3dbie-webgpu-poc](https://vault1.magland.org/proof-of-concept/fmm3dbie-webgpu-poc).
1012
1113 ## What is being compared
@@ -35,6 +37,10 @@ The fmm3dbie backend needs gfortran, make and OpenBLAS; the GPU backend needs an
3537 python3 compare.py
3638 python3 plot.py results/small-gpu
3739
40+ python3 run.py --problem cymbal --backend gpu --case large --selftest --reciprocity
41+ python3 run.py --problem cymbal --backend fmm3dbie --case large --selftest --reciprocity
42+ python3 plot.py results/cymbal-large-gpu
43+
3844 Each run writes `results/<case>-<backend>/summary.json`, with timings, iteration counts, errors and machine information, and `field.npz`. `compare.py` tabulates all runs and, for every case run with both backends, reports the speedup and the relative difference of the two wall fields. It also writes `results/comparison.md`. The runs do not have to be on the same machine: copy the result directories together before comparing. For a fair speedup, though, run fmm3dbie on the same machine as the GPU and with all of its cores (`OMP_NUM_THREADS`).
3945
4046 The cases are
@@ -67,11 +73,46 @@ The two codes need the same number of GMRES iterations, which is consistent with
6773
6874 At the default GMRES tolerance both codes' self-test errors are close to the tolerance. With a GMRES tolerance of 1e-10 the medium-case GPU error is 2.7e-7 with the single precision far field and 2.9e-10 with `--far64`. In an earlier run on a laptop, fmm3dbie reached 1.9e-8 there with eps 1e-7, at about four times the cost of eps 1e-4 (400 s against 92 s).
6975
76+## The cymbal problem
77+
78+The room problem puts the unknowns on a closed surface around the region of interest. A common acoustic question has the opposite structure: a small obstacle in free space, and a listener. Here we consider a listener near a cymbal, which we treat as stationary and sound-hard. We place a point source at the listener's position L and solve Δu + k²u = −δ(x − L) outside the cymbal, with ∂u/∂n = 0 on it and the radiation condition at infinity. By reciprocity, u(x) is also the pressure at the listener due to a unit source at x, so one solve gives the listener's response to a source anywhere around the cymbal. We evaluate u on a 3D grid in a bounding box (a 1 m cube around the cymbal) and on three planes through it. The box only selects where the field is evaluated: the integral representation satisfies the radiation condition, so, unlike a volume discretization, the problem needs no truncation or absorbing layer.
79+
80+fmm3dbie's integral equations need a closed surface. A zero-thickness cymbal would need an open-surface formulation, with hypersingular operators and a density that is singular at the edge, which fmm3dbie does not provide. We therefore model the cymbal as a thin closed body (`fmm3dbie_gpu/cymbal.py`). An axisymmetric midsurface with a raised bell and a shallow bow is offset by ±T/2 along its normal, where the thickness T decreases from 20 mm at the center to 10 mm at the rim, and a half-torus of diameter T closes the body at the rim. The cymbal is 40 cm across, so the rim thickness is 1/40 of the diameter; a real cymbal is closer to 1/400, and `--t-rim` sets it. The faces are meshed with an O-grid (an equiangular central square of at least 2 × 2 patches and rings of patches around it), graded toward the rim and refined across the bell. Every patch lies on one face or on the rim, so each is analytic, although the curvature jumps where the faces meet the rim. The rings are not conforming, which neither code requires. Node derivatives are computed by the complex step method, and at order 8 the surface reproduces the body's volume to 3e-9.
81+
82+The equation is the exterior mode of the same solver, `helm_rpcomb_imp_solver` with ifinout = 1 and λ = 0, and the GPU code needed no change to its operator for it. The GPU code uses the interpolated geometry for the cymbal, as fmm3dbie does, since no exact chart has been written for it in CUDA. The cases are in meters, with c = 343 m/s.
83+
84+| case | frequency | diameter | target patch size | order | patches | points |
85+|---|---|---|---|---|---|---|
86+| small | 1 kHz | 1.2 wavelengths | 6 cm | 6 | 388 | 19,012 |
87+| medium | 3 kHz | 3.5 wavelengths | 4 cm | 8 | 444 | 35,964 |
88+| large | 8 kHz | 9.3 wavelengths | 3 cm | 8 | 524 | 42,444 |
89+
90+`--freq`, `--h`, `--norder`, `--t-rim` and `--listener` override them. The listener is at (0.25, 0, 0.45) m by default, above the cymbal and outside its rim.
91+
92+There are three checks. The self-test uses about 5,000 point sources spread 3 mm apart over the midsurface out to r = 12 cm, with smoothly tapering strengths; their field is an exact solution outside. A single interior source would be at most 1 cm from the surface, and its data would vary on that scale, which the main problem's data does not. The energy balance compares the power emitted by the source, k/(4π) + Im v(L), with the power radiated to infinity, k∫|u∞|²dΩ, with the far-field pattern u∞ computed from the densities at the nodes. A rigid cymbal absorbs nothing, so the two must agree. `--reciprocity` adds a solve with the source at a point B below the cymbal and compares the scattered field at B due to the listener with the field at the listener due to B.
93+
94+Both backends were run on the same workstation as the room problem (Xeon w5-3435X with all 32 threads for fmm3dbie, eps = 1e-4; RTX 6000 Ada for the GPU), with GMRES tolerance 1e-6.
95+
96+| case | fmm3dbie solve | GPU setup + solve | speedup | GMRES iterations | self-test error, fmm3dbie / GPU | reciprocity, fmm3dbie / GPU | surface field difference |
97+|---|---|---|---|---|---|---|---|
98+| small | 27.0 s | 0.78 + 0.99 s | 15x | 25 | 4.3e-6 / 1.5e-5 | 2.0e-5 / 7.7e-5 | 6.6e-5 |
99+| medium | 71.9 s | 1.31 + 1.89 s | 23x | 32 | 9.2e-7 / 1.1e-6 | 2.0e-6 / 5.2e-6 | 2.4e-6 |
100+| large | 80.6 s | 1.50 + 1.82 s | 24x | 24 | 2.4e-5 / 4.6e-6 | 3.5e-6 / 7.0e-7 | 9.7e-6 |
101+
102+The energy mismatch was below 1.1e-6 for every run. The GMRES iteration counts are the same for both codes and do not grow with the frequency, which suggests that at this thickness the thin body does not degrade the conditioning. On the small case the GPU's discretization error is about three times fmm3dbie's. It does not depend on any of the GPU's quadrature parameters (near-field radius, quadtree refinement, rule orders, far-field oversampling) or on the far-field precision, and on a finer mesh (h = 4 cm, order 8, at 1 kHz) the two codes agree to 3e-6, so we take it to be a difference between two consistent discretizations rather than a quadrature error; we have not found its source. On the large case the GPU is the more accurate of the two. Note that the first version of the mesh had a single central patch at h = 6 cm, which left thin slivers of the first ring at the square's corners; with at least 2 × 2 central patches the small-case errors of both codes fell by a factor of 5 to 20.
103+
104+
105+
106+The figure shows the large case. The bottom row is the level relative to the free field: 0 dB where the cymbal has no effect, about +6 dB in the interference fringes above it, and a shadow of more than 15 dB below it.
107+
108+Several limitations are specific to this problem. The rim thickness is 1/40 of the diameter, ten times a real cymbal's; thinner bodies need more near-quadrature refinement, and we have not checked how the iteration counts behave. The rim is rounded, whereas a real cymbal's edge is closer to square. The cymbal is rigid, so its own vibration, which is how it makes sound, is not modeled.
109+
70110 ## Files
71111
72112 - `run.py`, `compare.py`, `plot.py`: run one backend, compare the results, and plot a field.
73-- `problem.py`: the showcase problem, shared by both backends: cases, source positions, right-hand sides and checks.
113+- `problems/room.py`, `problems/cymbal.py`: the two showcase problems, shared by both backends: cases, surfaces, sources, self-tests and checks. `problems/common.py` has the Green's function and other shared helpers.
74114 - `fmm3dbie_gpu/surface.py`: patch discretization and the rounded cube. It is pure NumPy.
115+- `fmm3dbie_gpu/cymbal.py`: the cymbal's shape and its patch mesh, also pure NumPy. Both backends build the cymbal from it.
75116 - `fmm3dbie_gpu/cuda_ops.py`: the CUDA kernels: near pairs, quadtrees, near weights, far-field sum and near-field product.
76117 - `fmm3dbie_gpu/solver.py`: the operator, evaluation, and GMRES, following fmm3dbie's `zgmres_guru`.
77118 - `reference/`: the ctypes interface to fmm3dbie and the same rounded-cube surface built with fmm3dbie's routines. We use it because the Python bindings shipped with fmm3dbie need `numpy.distutils`, which NumPy 2 no longer has.
compare.pymodified+13−12View file
@@ -3,8 +3,8 @@
33 python compare.py [results_dir]
44
55 Prints a table of all runs found in results_dir (default results/) and, for
6-every case run with both backends, the speedup and the relative difference of
7-the wall fields. The table is also written to results_dir/comparison.md.
6+every problem and case run with both backends, the speedup and the relative
7+difference of the surface fields. The table is also written to results_dir/comparison.md.
88 """
99
1010 import glob
@@ -25,30 +25,31 @@ def main():
2525 for f in sorted(glob.glob(os.path.join(rdir, "*", "summary.json"))):
2626 s = json.load(open(f))
2727 s["dir"] = os.path.dirname(f)
28- runs[(s["case"], s["backend"])] = s
28+ runs[(s.get("problem", "room"), s["case"], s["backend"])] = s
2929 if not runs:
3030 print(f"no results in {rdir}")
3131 return
3232
33- lines = ["| case | k | points | backend | machine | setup (s) | solve (s) | setup + solve (s) | "
33+ lines = ["| problem | case | k | points | backend | machine | setup (s) | solve (s) | setup + solve (s) | "
3434 "iterations | self-test error | energy mismatch |",
35- "|---|---|---|---|---|---|---|---|---|---|---|"]
35+ "|---|---|---|---|---|---|---|---|---|---|---|---|"]
3636 order = {"small": 0, "medium": 1, "large": 2}
37- for (case, backend), s in sorted(runs.items(), key=lambda kv: (order.get(kv[0][0], 9), kv[0][1])):
37+ sort_key = lambda pc: (pc[0] != "room", pc[0], order.get(pc[1], 9))
38+ for (problem, case, backend), s in sorted(runs.items(), key=lambda kv: (sort_key(kv[0][:2]), kv[0][2])):
3839 m = s["machine"]
3940 where = m.get("gpu", "").split(",")[0] if backend == "gpu" else f"{m.get('cpu', '?')}, {m.get('cpu_count')} threads"
40- lines.append(f"| {case} | {s['k']:g} | {s['npts']:,} | {backend} | {where} | {s['setup_s']:.2f} | "
41+ lines.append(f"| {problem} | {case} | {s['k']:.4g} | {s['npts']:,} | {backend} | {where} | {s['setup_s']:.2f} | "
4142 f"{s['solve_s']:.2f} | {s['setup_plus_solve_s']:.2f} | {s['niter']} | "
4243 f"{fmt(s.get('selftest_err'), '.1e')} | {s['energy_mismatch']:.1e} |")
4344
44- lines += ["", "| case | speedup (setup + solve) | wall field difference |", "|---|---|---|"]
45- for case in sorted({c for c, _ in runs}, key=lambda c: order.get(c, 9)):
46- if (case, "gpu") in runs and (case, "fmm3dbie") in runs:
47- g, c = runs[(case, "gpu")], runs[(case, "fmm3dbie")]
45+ lines += ["", "| problem | case | speedup (setup + solve) | surface field difference |", "|---|---|---|---|"]
46+ for problem, case in sorted({pc[:2] for pc in runs}, key=sort_key):
47+ if (problem, case, "gpu") in runs and (problem, case, "fmm3dbie") in runs:
48+ g, c = runs[(problem, case, "gpu")], runs[(problem, case, "fmm3dbie")]
4849 ug = np.load(os.path.join(g["dir"], "field.npz"))["u_wall"]
4950 uc = np.load(os.path.join(c["dir"], "field.npz"))["u_wall"]
5051 diff = np.linalg.norm(ug - uc) / np.linalg.norm(uc) if ug.shape == uc.shape else float("nan")
51- lines.append(f"| {case} | {c['setup_plus_solve_s'] / g['setup_plus_solve_s']:.0f}x | {diff:.1e} |")
52+ lines.append(f"| {problem} | {case} | {c['setup_plus_solve_s'] / g['setup_plus_solve_s']:.0f}x | {diff:.1e} |")
5253 text = "\n".join(lines)
5354 print(text)
5455 with open(os.path.join(rdir, "comparison.md"), "w") as f:
docs/cymbal_slices_large.pngadded+0−0View file
Binary file not shown.
fmm3dbie_gpu/cymbal.pyadded+246−0View file
@@ -0,0 +1,246 @@
1+"""A cymbal as a thin closed body built from quadrilateral patches (pure NumPy).
2+
3+fmm3dbie's integral equations need a closed surface, so the cymbal is not a
4+zero-thickness shell but a thin body around an axisymmetric midsurface
5+z = F(rho), rho = x^2 + y^2, with thickness T(rho). The midsurface has a raised
6+bell at the center and a shallow bow outside it; the thickness decreases from
7+the center to the rim, with zero slope at the rim. The top and bottom faces
8+are the offset surfaces M +- (T/2) N, with M the midsurface and N its upward
9+unit normal, and a half-torus of diameter T closes the body at the rim. Every
10+patch lies entirely on one face or on the rim, so every patch is analytic even
11+though the curvature jumps where the faces meet the rim.
12+
13+The faces are meshed with an O-grid in the (x, y) plane: an equiangular central
14+square of m_c x m_c patches, then rings of patches between the square and a
15+circle and between successive circles. Rings are graded toward the rim, where
16+the density varies on the scale of the thickness. Neighboring rings need not
17+have the same number of patches around, so the mesh is not conforming; neither
18+fmm3dbie nor the GPU code requires that it be.
19+
20+Node derivatives dx/du, dx/dv are computed by the complex step method, which is
21+exact to rounding for this analytic parametrization.
22+"""
23+
24+import numpy as np
25+
26+from . import surface
27+
28+CSTEP = 1e-30
29+
30+
31+class CymbalShape:
32+ """Axisymmetric cymbal shape; lengths in meters.
33+
34+ radius: outer radius, including the rim.
35+ t_rim: thickness at the rim; taper: thickness at the center / t_rim.
36+ bell_radius, bell_height, bell_width: the bell is a logistic step in rho,
37+ of height bell_height, centered at r = bell_radius, with width about
38+ bell_width in r.
39+ bow_height: fall of the midsurface from the center to the rim, apart from the bell.
40+ """
41+
42+ def __init__(self, radius=0.2, t_rim=0.01, taper=2.0, bell_radius=0.06, bell_height=0.025,
43+ bell_width=0.01, bow_height=0.015):
44+ self.radius, self.t_rim, self.taper = radius, t_rim, taper
45+ self.bell_radius, self.bell_height, self.bell_width = bell_radius, bell_height, bell_width
46+ self.bow_height = bow_height
47+ self.r_mid = radius - t_rim / 2 # radius of the midsurface, where the rim is attached
48+
49+ # midsurface height, its rho-derivative, and the thickness, all functions of rho
50+ def _bell_arg(self, rho):
51+ return (rho - self.bell_radius ** 2) / (2 * self.bell_radius * self.bell_width)
52+
53+ def F(self, rho):
54+ return (self.bow_height * (1 - rho / self.r_mid ** 2)
55+ + self.bell_height / (1 + np.exp(self._bell_arg(rho))))
56+
57+ def dF(self, rho):
58+ e = np.exp(self._bell_arg(rho))
59+ return (-self.bow_height / self.r_mid ** 2
60+ - self.bell_height * e / (1 + e) ** 2 / (2 * self.bell_radius * self.bell_width))
61+
62+ def T(self, rho):
63+ return self.t_rim * (1 + (self.taper - 1) * (1 - rho / self.r_mid ** 2) ** 2)
64+
65+ def face_point(self, x, y, side):
66+ """Point of the top (side = +1) or bottom (side = -1) face above/below (x, y) of
67+ the midsurface. Complex-safe."""
68+ rho = x * x + y * y
69+ g = 2 * self.dF(rho) # grad F = g (x, y)
70+ q = np.sqrt(1 + g * g * rho)
71+ s = side * self.T(rho) / 2 / q
72+ return np.array([x - s * g * x, y - s * g * y, self.F(rho) + s])
73+
74+ def rim_point(self, theta, phi):
75+ """Point of the rim at azimuth theta; phi from -pi/2 (bottom face) to pi/2 (top face)."""
76+ R = self.r_mid
77+ rho = R * R
78+ fp = 2 * R * self.dF(rho) # dz/dr of the midsurface at the rim
79+ q = np.sqrt(1 + fp * fp)
80+ c, s = np.cos(theta), np.sin(theta)
81+ h = self.T(rho) / 2
82+ a, b = h * np.cos(phi) / q, h * np.sin(phi) / q # along the tangent (1, fp) and normal (-fp, 1)
83+ rr = R + a - b * fp
84+ return np.array([rr * c, rr * s, self.F(rho) + a * fp + b])
85+
86+ def cross_section(self, n=400):
87+ """Closed outline of the body in the (r, z) half plane, (2, m), counterclockwise."""
88+ r = np.linspace(0, self.r_mid, n)
89+ top = self.face_point(r, 0 * r, 1)
90+ bot = self.face_point(r, 0 * r, -1)
91+ phi = np.linspace(-np.pi / 2, np.pi / 2, n // 4)
92+ rim = self.rim_point(0 * phi, phi)
93+ rz = np.hstack([bot[[0, 2]], rim[[0, 2]], top[[0, 2]][:, ::-1]])
94+ return rz
95+
96+ def inside(self, pts, margin=0.0):
97+ """True for points (3, n) inside the body or within margin of its surface.
98+
99+ Every vertical line meets the cross section in one interval, so we offset the
100+ outline outward by margin, split it at its outermost point into a lower and an
101+ upper chain, and compare z with both chains at the point's radius."""
102+ rz = self.cross_section(4000)
103+ if margin > 0:
104+ t = np.roll(rz, -1, axis=1) - np.roll(rz, 1, axis=1) # counterclockwise tangent
105+ t /= np.linalg.norm(t, axis=0)
106+ rz = rz + margin * np.array([t[1], -t[0]]) # outward normal
107+ rz[0] = np.maximum(rz[0], 0)
108+ i = int(np.argmax(rz[0]))
109+ lo, up = rz[:, :i + 1], rz[:, i:][:, ::-1]
110+ r = np.hypot(pts[0], pts[1])
111+ z_lo = np.interp(r, np.maximum.accumulate(lo[0]), lo[1])
112+ z_up = np.interp(r, np.maximum.accumulate(up[0]), up[1])
113+ return (r <= rz[0, i]) & (pts[2] > z_lo) & (pts[2] < z_up)
114+
115+ def volume(self, n=20000):
116+ """Volume of the body by Pappus' theorem on the cross section (a check)."""
117+ rz = self.cross_section(n)
118+ r, z = rz
119+ r1, z1 = np.roll(r, -1), np.roll(z, -1)
120+ # int r dA over the cross section, by Green's theorem: (1/2) oint r^2 dz
121+ return 2 * np.pi * np.sum((r * r + r * r1 + r1 * r1) / 6 * (z1 - z))
122+
123+
124+def ring_layout(shape, h, grade=1.5, c_frac=0.12):
125+ """Radii and patch counts of the O-grid for a target patch size h.
126+
127+ Returns (m_c, c, radii, counts): the central square has half-width c and
128+ m_c x m_c patches (m_c >= 2); ring i lies between radii[i] and radii[i + 1]
129+ (radii[0] = c means the square) and has counts[i] patches around (a multiple
130+ of 4). The ring around the square reaches r = 2c, so that it is 0.59c wide at
131+ the square's corners. Radial widths are at most h; they grow geometrically
132+ from t_rim at the rim by the factor grade, and they are at most twice the
133+ bell width across the bell.
134+ """
135+ R = shape.r_mid
136+ c = c_frac * R
137+ m_c = max(2, int(np.ceil(2 * c / h)))
138+ first = 2 * c # outer radius of the ring around the square
139+ b0 = shape.bell_radius - 2.5 * shape.bell_width
140+ b1 = shape.bell_radius + 2.5 * shape.bell_width
141+
142+ def wmax(edge):
143+ w = min(h, shape.t_rim + (grade - 1) * (R - edge))
144+ if edge > b0 and edge - w < b1:
145+ w = min(w, 2 * shape.bell_width)
146+ return w
147+
148+ radii = [R]
149+ while True:
150+ w = wmax(radii[-1])
151+ if radii[-1] - w < first + 0.5 * w:
152+ break
153+ radii.append(radii[-1] - w)
154+ radii = np.array([c, first] + radii[::-1])
155+ counts = []
156+ for i in range(len(radii) - 1):
157+ width = radii[i + 1] - radii[i]
158+ arc = min(h, 2.5 * width)
159+ n = 4 * max(m_c, int(np.ceil(2 * np.pi * radii[i + 1] / (4 * arc))))
160+ counts.append(n)
161+ counts[0] = 4 * m_c # the first ring matches the square
162+ # rings never have fewer patches around than the ring inside them
163+ counts = list(np.maximum.accumulate(counts))
164+ return m_c, c, radii, counts
165+
166+
167+def _face_charts(m_c, c, radii, counts):
168+ """Planar charts (u, v) in [-1, 1]^2 -> (x, y) for all face patches, as closures.
169+ Each chart is oriented so that d/du x d/dv points up."""
170+ charts = []
171+ q = np.pi / 4
172+ # central square, equiangular so that its edges match the first ring
173+ br = np.linspace(-1, 1, m_c + 1)
174+ for i in range(m_c):
175+ for j in range(m_c):
176+ def ch(u, v, a0=br[i], a1=br[i + 1], b0=br[j], b1=br[j + 1]):
177+ s = a0 + 0.5 * (u + 1) * (a1 - a0)
178+ t = b0 + 0.5 * (v + 1) * (b1 - b0)
179+ return c * np.tan(q * s), c * np.tan(q * t)
180+ charts.append(ch)
181+ # rings; u is radial, v is the angle (counterclockwise)
182+ for ir in range(len(radii) - 1):
183+ r0, r1, n = radii[ir], radii[ir + 1], counts[ir]
184+ dth = 2 * np.pi / n
185+ for k in range(n):
186+ th0 = -q + k * dth # patch angles start at the square corner (-pi/4)
187+ if ir == 0:
188+ def ch(u, v, th0=th0, r1=r1, dth=dth):
189+ th = th0 + 0.5 * (v + 1) * dth
190+ sec = np.floor((th.real + q) / (2 * q)) # sector index 0..3
191+ loc = th - sec * 2 * q # angle relative to the sector's center
192+ rin = c / np.cos(loc)
193+ a = 0.5 * (u + 1)
194+ rad = (1 - a) * rin + a * r1
195+ return rad * np.cos(th), rad * np.sin(th)
196+ else:
197+ def ch(u, v, th0=th0, r0=r0, r1=r1, dth=dth):
198+ th = th0 + 0.5 * (v + 1) * dth
199+ rad = r0 + 0.5 * (u + 1) * (r1 - r0)
200+ return rad * np.cos(th), rad * np.sin(th)
201+ charts.append(ch)
202+ return charts
203+
204+
205+def _node_data(fun, uvs):
206+ """xyz, r_u, r_v, unit normal (12, n) of the map fun(u, v) -> (3, n), with complex step derivatives."""
207+ u, v = uvs[0].astype(complex), uvs[1].astype(complex)
208+ x = np.real(fun(u, v))
209+ du = np.imag(fun(u + 1j * CSTEP, v)) / CSTEP
210+ dv = np.imag(fun(u, v + 1j * CSTEP)) / CSTEP
211+ nrm = np.cross(du.T, dv.T).T
212+ nrm /= np.linalg.norm(nrm, axis=0)
213+ return np.vstack([x, du, dv, nrm])
214+
215+
216+def patch_maps(shape, h, grade=1.5):
217+ """All patch maps (u, v) -> xyz (3, n) with outward orientation, and the ring layout."""
218+ m_c, c, radii, counts = ring_layout(shape, h, grade)
219+ charts = _face_charts(m_c, c, radii, counts)
220+ maps = []
221+ for ch in charts: # top faces
222+ maps.append(lambda u, v, ch=ch: shape.face_point(*ch(u, v), 1))
223+ for ch in charts: # bottom faces: swapping u and v flips the normal
224+ maps.append(lambda u, v, ch=ch: shape.face_point(*ch(v, u), -1))
225+ nrim = counts[-1]
226+ for k in range(nrim): # rim: u is the azimuth, v the angle over the rim
227+ th0 = -np.pi / 4 + 2 * np.pi * k / nrim
228+ maps.append(lambda u, v, th0=th0: shape.rim_point(th0 + 0.5 * (u + 1) * 2 * np.pi / nrim,
229+ 0.5 * np.pi * v))
230+ return maps, dict(m_c=m_c, c=c, radii=radii, counts=counts, nrim=nrim)
231+
232+
233+def cymbal_srcvals(shape, h, norder, grade=1.5):
234+ """Node data (12, npts) in the fmm3dbie iptype 11 convention, and the layout."""
235+ uvs, _ = surface.disc_nodes(norder)
236+ maps, layout = patch_maps(shape, h, grade)
237+ return np.hstack([_node_data(f, uvs) for f in maps]), layout
238+
239+
240+def cymbal(shape, h, norder, grade=1.5):
241+ """A surface.Surface for the cymbal, with the interpolated geometry inside the
242+ integrals, as in fmm3dbie."""
243+ srcvals, layout = cymbal_srcvals(shape, h, norder, grade)
244+ s = surface.Surface(norder, srcvals)
245+ s.layout = layout
246+ return s
plot.pymodified+95−2View file
@@ -1,7 +1,9 @@
1-"""Plot the room field saved by run.py (python plot.py results/<case>-<backend>): slices of Re(u) and level in dB,
2-plus an animation of the real pressure field Re(u exp(-i omega t)) over one period."""
1+"""Plot a field saved by run.py (python plot.py results/<case>-<backend>): slices of Re(u) and level in dB,
2+plus an animation of the real pressure field Re(u exp(-i omega t)) over one period. Works for the room
3+and the cymbal problems."""
34
45 import argparse
6+import json
57 import os
68
79 import matplotlib
@@ -25,6 +27,9 @@ def main():
2527 args = ap.parse_args()
2628
2729 d = np.load(os.path.join(args.out, "field.npz"))
30+ with open(os.path.join(args.out, "summary.json")) as f:
31+ if json.load(f).get("problem", "room") == "cymbal":
32+ return plot_cymbal(d, args)
2833 g, x0, k = d["grid"], d["x0"], float(d["k"])
2934 ext = [g[0], g[-1], g[0], g[-1]]
3035 us = {n: d[f"u_{n}"] for n, *_ in SLICES}
@@ -94,5 +99,93 @@ def main():
9499 print(f"wrote {args.out}/slices.png and {args.out}/animation.gif")
95100
96101
102+CYMBAL_SLICES = ["xz", "yz", "xy"]
103+
104+
105+def plot_cymbal(d, args):
106+ """Slices through the evaluation box around the cymbal. u(x) is the pressure at the
107+ listener due to a unit source at x (reciprocity). The level is shown relative to
108+ the free field, 20 log10 |u / G|, which is 0 dB where the cymbal has no effect."""
109+ x0, k, freq = d["x0"], float(d["k"]), float(d["freq"])
110+ rz, R = d["outline"], float(d["radius"])
111+ names = [n for n in CYMBAL_SLICES if f"u_{n}" in d]
112+ panels = []
113+ for n in names:
114+ ga, gb = d[f"axes_{n}"]
115+ fixed, val = d[f"plane_{n}"]
116+ fixed = int(fixed)
117+ ia, ib = [i for i in range(3) if i != fixed]
118+ A, B = np.meshgrid(ga, gb, indexing="xy")
119+ pts = np.zeros((3, A.size))
120+ pts[ia], pts[ib], pts[fixed] = A.ravel(), B.ravel(), val
121+ r = np.linalg.norm(pts - x0[:, None], axis=0)
122+ G = (np.exp(1j * k * r) / (4 * np.pi * r)).reshape(A.shape)
123+ panels.append((n, ga, gb, ia, ib, fixed, val, d[f"u_{n}"], G, r.reshape(A.shape)))
124+ far = np.concatenate([np.abs(u[np.isfinite(u) & (r > np.pi / k)]) for *_, u, G, r in panels])
125+ umax = np.percentile(far, 99)
126+
127+ def outline(ax, ia, ib, fixed, val):
128+ if fixed == 2: # horizontal plane: the cymbal's rim, projected
129+ t = np.linspace(0, 2 * np.pi, 200)
130+ ax.plot(R * np.cos(t), R * np.sin(t), color="0.2", lw=0.8, ls="--")
131+ elif abs(val) < R: # vertical plane: the cross section, at distance |val| from the axis
132+ r, z = rz
133+ m = r >= abs(val)
134+ s = np.sqrt(np.maximum(r[m] ** 2 - val ** 2, 0))
135+ for sign in (1, -1):
136+ ax.fill(sign * s, z[m], color="0.25", lw=0)
137+
138+ fig, axes = plt.subplots(2, len(panels), figsize=(4.4 * len(panels), 8.4), constrained_layout=True,
139+ squeeze=False)
140+ for col, (n, ga, gb, ia, ib, fixed, val, u, G, r) in enumerate(panels):
141+ ext = [ga[0], ga[-1], gb[0], gb[-1]]
142+ ax = axes[0, col]
143+ im0 = ax.imshow(masked(u.real), origin="lower", extent=ext, cmap="RdBu_r",
144+ vmin=-umax, vmax=umax, interpolation="bilinear")
145+ ax.set_title(f"Re u, {n}-plane at {'xyz'[fixed]} = {val:.2f} m")
146+ ax = axes[1, col]
147+ with np.errstate(divide="ignore", invalid="ignore"):
148+ level = 20 * np.log10(np.abs(u) / np.abs(G))
149+ im1 = ax.imshow(masked(level), origin="lower", extent=ext, cmap="PuOr_r",
150+ vmin=-15, vmax=15, interpolation="bilinear")
151+ ax.set_title(f"level re free field (dB), {n}-plane")
152+ for ax in axes[:, col]:
153+ outline(ax, ia, ib, fixed, val)
154+ if abs(x0[fixed] - val) < 1e-9:
155+ ax.plot(x0[ia], x0[ib], marker="o", ms=5, mfc="none", mec="0.15", mew=1)
156+ ax.set_xlabel(f"{'xyz'[ia]} (m)")
157+ ax.set_ylabel(f"{'xyz'[ib]} (m)")
158+ ax.set_aspect("equal")
159+ fig.colorbar(im0, ax=axes[0, :], shrink=0.85, label="Re u")
160+ fig.colorbar(im1, ax=axes[1, :], shrink=0.85, label="dB")
161+ wl = 2 * np.pi / k
162+ fig.suptitle(f"Listener at ({x0[0]:.2f}, {x0[1]:.2f}, {x0[2]:.2f}) m near a rigid cymbal, "
163+ f"{freq:g} Hz (diameter = {2 * R / wl:.2f} wavelengths)\nBy reciprocity u(x) is the "
164+ f"pressure at the listener due to a unit source at x")
165+ fig.savefig(os.path.join(args.out, "slices.png"), dpi=130)
166+ plt.close(fig)
167+
168+ # Animation of the vertical plane through the listener
169+ n, ga, gb, ia, ib, fixed, val, u, G, r = panels[0]
170+ fig, ax = plt.subplots(figsize=(6, 5.4), constrained_layout=True)
171+ im = ax.imshow(masked(u.real), origin="lower", extent=[ga[0], ga[-1], gb[0], gb[-1]], cmap="RdBu_r",
172+ vmin=-umax, vmax=umax, interpolation="bilinear")
173+ outline(ax, ia, ib, fixed, val)
174+ ax.set_xlabel(f"{'xyz'[ia]} (m)")
175+ ax.set_ylabel(f"{'xyz'[ib]} (m)")
176+ fig.colorbar(im, ax=ax, shrink=0.85, label="Re(u e^{-iωt})")
177+ title = ax.set_title("")
178+
179+ def frame(i):
180+ im.set_data(masked((u * np.exp(-2j * np.pi * i / args.frames)).real))
181+ title.set_text(f"{n}-plane, t/T = {i / args.frames:.2f}")
182+ return [im]
183+
184+ anim = FuncAnimation(fig, frame, frames=args.frames, blit=False)
185+ anim.save(os.path.join(args.out, "animation.gif"), writer=PillowWriter(fps=12), dpi=80)
186+ plt.close(fig)
187+ print(f"wrote {args.out}/slices.png and {args.out}/animation.gif")
188+
189+
97190 if __name__ == "__main__":
98191 main()
problem.pydeleted+0−114View file
@@ -1,114 +0,0 @@
1-"""The showcase problem, shared by both backends.
2-
3-A 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
5-exp(-i omega t) the pressure amplitude u satisfies
6-
7- (Delta + k^2) u = -delta(x - x0) in the room,
8- du/dn + i k lambda u = 0 on the walls,
9-
10-with n the outward normal and lambda = -beta, where beta > 0 is the wall
11-admittance. We write u = G + v with G the free-space field of the source, and
12-find v from the right-preconditioned combined field integral equation of
13-fmm3dbie (helm_rpcomb_imp_solver, interior mode).
14-
15-This module has no dependence on either backend.
16-"""
17-
18-import json
19-import os
20-import platform
21-import socket
22-import subprocess
23-
24-import numpy as np
25-
26-# Named cases. Both backends discretize the surface identically, so the same
27-# case gives the same unknowns on both.
28-CASES = {
29- "small": dict(k=10.0, beta=0.2, m=4, norder=6), # 3.2 wavelengths, 4,704 points
30- "medium": dict(k=10.0, beta=0.2, m=6, norder=8), # 3.2 wavelengths, 17,496 points
31- "large": dict(k=40.0, beta=0.05, m=10, norder=8), # 12.7 wavelengths, 48,600 points
32-}
33-A, P = 1.0, 8 # room half-width and superellipsoid exponent
34-SOURCE = np.array([-0.35, 0.2, 0.1]) # the point source
35-EXTERIOR_SOURCE = np.array([1.7, 0.6, -0.5]) # for the self-test (outside the room)
36-ALPHA = 1.0
37-
38-
39-def green(x, x0, k):
40- d = x - x0[:, None]
41- r = np.linalg.norm(d, axis=0)
42- return np.exp(1j * k * r) / (4 * np.pi * r)
43-
44-
45-def green_dn(x, nrm, x0, k):
46- d = x - x0[:, None]
47- r = np.linalg.norm(d, axis=0)
48- g = np.exp(1j * k * r) / (4 * np.pi * r)
49- return g * (1j * k - 1 / r) * np.sum(d * nrm, axis=0) / r
50-
51-
52-def inside(x, a=A, p=P, margin=0.0):
53- return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p
54-
55-
56-def room_rhs(xs, nrm, k, lam):
57- """Impedance data for v = u - G: -(dG/dn + i k lambda G)."""
58- return -(green_dn(xs, nrm, SOURCE, k) + 1j * k * lam * green(xs, SOURCE, k))
59-
60-
61-def selftest_rhs(xs, nrm, k, lam):
62- """Impedance data of the field of a source outside the room."""
63- return green_dn(xs, nrm, EXTERIOR_SOURCE, k) + 1j * k * lam * green(xs, EXTERIOR_SOURCE, k)
64-
65-
66-def selftest_points(n=50):
67- return np.random.default_rng(0).uniform(-0.8 * A, 0.8 * A, size=(3, n))
68-
69-
70-def slice_points(ngrid):
71- """Points on three planes through the source; returns grid, {name: (pts, mask, shape)}."""
72- g = np.linspace(-A, A, ngrid)
73- out = {}
74- for name, (i, j, fixed) in {"xy": (0, 1, 2), "xz": (0, 2, 1), "yz": (1, 2, 0)}.items():
75- Ag, Bg = np.meshgrid(g, g, indexing="xy")
76- pts = np.zeros((3, Ag.size))
77- pts[i], pts[j], pts[fixed] = Ag.ravel(), Bg.ravel(), SOURCE[fixed]
78- out[name] = (pts, inside(pts, margin=1e-3), Ag.shape)
79- return g, out
80-
81-
82-def energy_balance(k, lam, v_at_source, wts, u_wall):
83- """Power radiated by the source, k/(4 pi) + Im v(x0), and power absorbed by the
84- walls, -k Re(lambda) int |u|^2 dS. They must agree."""
85- p_src = k / (4 * np.pi) + v_at_source.imag
86- p_wall = -k * np.real(lam) * np.sum(wts * np.abs(u_wall) ** 2)
87- return p_src, p_wall
88-
89-
90-def machine_info():
91- info = dict(host=socket.gethostname(), platform=platform.platform(), cpu_count=os.cpu_count())
92- try:
93- with open("/proc/cpuinfo") as f:
94- for line in f:
95- if line.startswith("model name"):
96- info["cpu"] = line.split(":", 1)[1].strip()
97- break
98- except OSError:
99- pass
100- try:
101- gpu = subprocess.run(["nvidia-smi", "--query-gpu=name,memory.total", "--format=csv,noheader"],
102- capture_output=True, text=True, timeout=10)
103- if gpu.returncode == 0:
104- info["gpu"] = gpu.stdout.strip().splitlines()[0]
105- except (OSError, subprocess.SubprocessError):
106- pass
107- return info
108-
109-
110-def save(outdir, summary, **arrays):
111- os.makedirs(outdir, exist_ok=True)
112- np.savez(os.path.join(outdir, "field.npz"), **arrays)
113- with open(os.path.join(outdir, "summary.json"), "w") as f:
114- json.dump(summary, f, indent=2)
problems/__init__.pyadded+11−0View file
@@ -0,0 +1,11 @@
1+"""The showcase problems, shared by both backends: room (interior impedance
2+problem in a rounded cube) and cymbal (exterior sound-hard problem around a thin
3+cymbal)."""
4+
5+import importlib
6+
7+PROBLEMS = ["room", "cymbal"]
8+
9+
10+def get(name):
11+ return importlib.import_module(f"problems.{name}")
problems/common.pyadded+71−0View file
@@ -0,0 +1,71 @@
1+"""Helpers shared by the problems: the Green's function, boundary data, machine
2+information and saving results. No dependence on either backend."""
3+
4+import json
5+import os
6+import platform
7+import socket
8+import subprocess
9+
10+import numpy as np
11+
12+
13+def green(x, x0, k):
14+ d = x - x0[:, None]
15+ r = np.linalg.norm(d, axis=0)
16+ return np.exp(1j * k * r) / (4 * np.pi * r)
17+
18+
19+def green_dn(x, nrm, x0, k):
20+ d = x - x0[:, None]
21+ r = np.linalg.norm(d, axis=0)
22+ g = np.exp(1j * k * r) / (4 * np.pi * r)
23+ return g * (1j * k - 1 / r) * np.sum(d * nrm, axis=0) / r
24+
25+
26+def impedance_data(xs, nrm, x0, k, lam):
27+ """du/dn + i k lambda u on the surface for u the field of a unit source at x0."""
28+ return green_dn(xs, nrm, x0, k) + 1j * k * lam * green(xs, x0, k)
29+
30+
31+def source_sum(x, x0s, ws, k, nrm=None, lam=0.0, chunk=64, xp=np):
32+ """Field sum_j ws[j] G(x, x0s[:, j]) at points x (3, n), or with nrm given its
33+ impedance data du/dn + i k lambda u. xp is numpy or cupy; the result is NumPy."""
34+ x, x0s, ws = xp.asarray(x), xp.asarray(x0s), xp.asarray(ws)
35+ nrm = None if nrm is None else xp.asarray(nrm)
36+ out = xp.zeros(x.shape[1], dtype=complex)
37+ for i in range(0, x0s.shape[1], chunk):
38+ d = x[:, None, :] - x0s[:, i:i + chunk, None] # (3, chunk, n)
39+ r = xp.sqrt(xp.sum(d * d, axis=0))
40+ g = xp.exp(1j * k * r) / (4 * np.pi * r)
41+ if nrm is not None:
42+ g = g * ((1j * k - 1 / r) * xp.einsum("icn,in->cn", d, nrm) / r + 1j * k * lam)
43+ out += ws[i:i + chunk] @ g
44+ return out.get() if xp is not np else out
45+
46+
47+def machine_info():
48+ info = dict(host=socket.gethostname(), platform=platform.platform(), cpu_count=os.cpu_count())
49+ try:
50+ with open("/proc/cpuinfo") as f:
51+ for line in f:
52+ if line.startswith("model name"):
53+ info["cpu"] = line.split(":", 1)[1].strip()
54+ break
55+ except OSError:
56+ pass
57+ try:
58+ gpu = subprocess.run(["nvidia-smi", "--query-gpu=name,memory.total", "--format=csv,noheader"],
59+ capture_output=True, text=True, timeout=10)
60+ if gpu.returncode == 0:
61+ info["gpu"] = gpu.stdout.strip().splitlines()[0]
62+ except (OSError, subprocess.SubprocessError):
63+ pass
64+ return info
65+
66+
67+def save(outdir, summary, **arrays):
68+ os.makedirs(outdir, exist_ok=True)
69+ np.savez(os.path.join(outdir, "field.npz"), **arrays)
70+ with open(os.path.join(outdir, "summary.json"), "w") as f:
71+ json.dump(summary, f, indent=2)
problems/cymbal.pyadded+245−0View file
@@ -0,0 +1,245 @@
1+"""The cymbal problem: a reciprocity solve for a listener near a rigid cymbal.
2+
3+A point source at the listener's position L radiates in free space; the only
4+obstacle is a stationary, sound-hard cymbal. With time dependence
5+exp(-i omega t) the pressure amplitude u satisfies
6+
7+ (Delta + k^2) u = -delta(x - L) outside the cymbal,
8+ du/dn = 0 on the cymbal,
9+
10+and the Sommerfeld radiation condition. By reciprocity, u(x) is also the
11+pressure at the listener due to a unit point source at x, so one solve gives
12+the listener's response to a source anywhere; we evaluate it on a 3D grid in a
13+bounding box, and on three planes through the box. The box only selects where
14+the field is evaluated: the integral representation satisfies the radiation
15+condition, so the problem needs no truncation.
16+
17+The cymbal is a thin closed body (fmm3dbie_gpu/cymbal.py): fmm3dbie's integral
18+equations need a closed surface. We write u = G + v, with G the free-space
19+field of the source, and find v from fmm3dbie's right-preconditioned combined
20+field 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.
22+Lengths are in meters and the speed of sound is 343 m/s.
23+
24+This module has no dependence on either backend.
25+"""
26+
27+import time
28+
29+import numpy as np
30+
31+from fmm3dbie_gpu.cymbal import CymbalShape, cymbal, cymbal_srcvals
32+
33+from .common import green, impedance_data
34+
35+NAME = "cymbal"
36+IFINOUT = 1 # exterior problem
37+ALPHA = 1.0
38+C_SOUND = 343.0
39+
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.
42+CASES = {
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+}
47+RADIUS = 0.2 # a 40 cm cymbal
48+T_RIM = 2 * RADIUS / 40 # rim thickness diameter/40; twice that at the center
49+LISTENER = np.array([0.25, 0.0, 0.45]) # above the cymbal and off its axis
50+BOX = np.array([[-0.5, 0.5], [-0.5, 0.5], [-0.35, 0.65]]) # evaluation region (1 m cube)
51+RECIPROCITY_POINT = np.array([-0.15, 0.05, -0.12]) # below the cymbal
52+SLICE_BELOW = -0.1 # height of the horizontal slice, below the cymbal
53+NGRID = 201 # default grid size of the slices, 5 mm spacing
54+
55+
56+def 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)")
70+
71+
72+def 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
82+
83+
84+def results_name(case_name, backend):
85+ return f"cymbal-{case_name}-{backend}"
86+
87+
88+def shape(case):
89+ return CymbalShape(radius=RADIUS, t_rim=case["t_rim"])
90+
91+
92+def 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})")
97+
98+
99+def gpu_surface(case, args):
100+ return cymbal(shape(case), case["h"], case["norder"])
101+
102+
103+def 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)
107+
108+
109+def source(case):
110+ return np.array(case["listener"])
111+
112+
113+def 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
129+
130+
131+def 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]
136+
137+
138+def 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)
164+
165+
166+def 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()])
170+
171+
172+def 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
184+
185+
186+def 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
192+
193+
194+def 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)
200+
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])
209+
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)
220+
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")
problems/room.pyadded+136−0View file
@@ -0,0 +1,136 @@
1+"""The room problem.
2+
3+A 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
5+exp(-i omega t) the pressure amplitude u satisfies
6+
7+ (Delta + k^2) u = -delta(x - x0) in the room,
8+ du/dn + i k lambda u = 0 on the walls,
9+
10+with n the outward normal and lambda = -beta, where beta > 0 is the wall
11+admittance. We write u = G + v with G the free-space field of the source, and
12+find v from the right-preconditioned combined field integral equation of
13+fmm3dbie (helm_rpcomb_imp_solver, interior mode).
14+
15+This module has no dependence on either backend.
16+"""
17+
18+import time
19+
20+import numpy as np
21+
22+from .common import green
23+
24+NAME = "room"
25+IFINOUT = 0 # interior problem
26+
27+# Named cases. Both backends discretize the surface identically, so the same
28+# case gives the same unknowns on both.
29+CASES = {
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+}
34+A, P = 1.0, 8 # room half-width and superellipsoid exponent
35+SOURCE = np.array([-0.35, 0.2, 0.1]) # the point source
36+EXTERIOR_SOURCE = np.array([1.7, 0.6, -0.5]) # for the self-test (outside the room)
37+ALPHA = 1.0
38+NGRID = 101 # default grid size of the field slices
39+
40+
41+def 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")
46+
47+
48+def 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
55+
56+
57+def results_name(case_name, backend):
58+ return f"{case_name}-{backend}"
59+
60+
61+def describe(case):
62+ return (f"k = {case['k']}, beta = {case['beta']}, "
63+ f"room width {2 * A * case['k'] / (2 * np.pi):.2f} wavelengths")
64+
65+
66+def 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)
70+
71+
72+def fmm3dbie_surface(case):
73+ import rounded_cube as rc
74+ return rc.rounded_cube(a=A, p=P, m=case["m"], norder=case["norder"])
75+
76+
77+def source(case):
78+ return SOURCE
79+
80+
81+def 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)
85+
86+
87+def selftest_points(case, n=50):
88+ return np.random.default_rng(0).uniform(-0.8 * A, 0.8 * A, size=(3, n))
89+
90+
91+def inside(x, a=A, p=P, margin=0.0):
92+ return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p
93+
94+
95+def 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
105+
106+
107+def 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
113+
114+
115+def 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")
reference/fmm3dbie_ctypes.pymodified+6−4View file
@@ -149,11 +149,13 @@ def helm_rpcomb_eval(surf, zk, alpha, sigma, siksigma, targs, eps=1e-6):
149149 return pot
150150
151151
152-def helm_rpcomb_eval_boundary(surf, zk, alpha, sigma, siksigma, eps=1e-6):
153- """Interior boundary limit of u = S_k[s] + i alpha D_k S_{i|k|}[s] at the nodes.
152+def helm_rpcomb_eval_boundary(surf, zk, alpha, sigma, siksigma, eps=1e-6, ifinout=0):
153+ """Boundary limit of u = S_k[s] + i alpha D_k S_{i|k|}[s] at the nodes, from the
154+ inside (ifinout = 0) or the outside (ifinout = 1).
154155
155156 The on-surface evaluation gives the principal value of D_k; with outward
156- normals the interior limit adds -1/2 times the double-layer density.
157+ normals the interior limit adds -1/2 times the double-layer density and the
158+ exterior limit +1/2 times it.
157159 """
158160 uvs = get_disc_nodes(surf.norder, surf.iptype)
159161 ipatch_id = np.repeat(np.arange(1, surf.npatches + 1), surf.npols).astype(np.int64)
@@ -165,4 +167,4 @@ def helm_rpcomb_eval_boundary(surf, zk, alpha, sigma, siksigma, eps=1e-6):
165167 _lib.helm_rpcomb_eval_(*surf.args(), _ref_i8(12), _ref_i8(surf.npts), _ptr(surf.srcvals),
166168 _ptr(ipatch_id), _ptr(uvs_targ), _ref_f8(eps), _ptr(zpars),
167169 _ptr(sigma), _ptr(siksigma), _ptr(pot))
168- return pot - 0.5j * alpha * siksigma
170+ return pot + (-1) ** (ifinout + 1) * 0.5j * alpha * siksigma
run.pymodified+72−70View file
@@ -1,11 +1,15 @@
1-"""Solve the showcase problem with one backend and save timings and fields.
1+"""Solve a showcase problem with one backend and save timings and fields.
22
33 python run.py --backend gpu --case small --selftest
44 python run.py --backend fmm3dbie --case small --selftest
5-
6-Results go to results/<case>-<backend>/ (summary.json, field.npz). The two
7-backends solve the same discrete problem on the same surface nodes; compare.py
8-puts their results side by side.
5+ python run.py --problem cymbal --backend fmm3dbie --case small --selftest
6+
7+The problems are `room` (the default; a point source in a room with impedance
8+walls) and `cymbal` (a listener near a rigid cymbal in free space). Results go
9+to results/<case>-<backend>/ for the room and results/cymbal-<case>-<backend>/
10+for the cymbal (summary.json, field.npz). The two backends solve the same
11+discrete problem on the same surface nodes; compare.py puts their results side
12+by side. `python run.py --problem cymbal --help` lists the cymbal's options.
913 """
1014
1115 import argparse
@@ -18,7 +22,8 @@ import numpy as np
1822 HERE = os.path.dirname(os.path.abspath(__file__))
1923 sys.path.insert(0, HERE)
2024
21-import problem as pb # noqa: E402
25+import problems # noqa: E402
26+from problems.common import impedance_data, machine_info, save, source_sum # noqa: E402
2227
2328
2429 class GPUBackend:
@@ -26,24 +31,22 @@ class GPUBackend:
2631
2732 name = "gpu"
2833
29- def __init__(self, case, args):
34+ def __init__(self, pb, case, args):
3035 import cupy as cp
31- from fmm3dbie_gpu import cuda_ops, solver, surface
32- self.cp = cp
36+ from fmm3dbie_gpu import cuda_ops, solver
37+ self.cp = self.xp = cp
3338 t0 = time.time()
3439 cuda_ops.module() # compile or load the CUDA kernels once, outside the timings
3540 self.kernel_load_s = time.time() - t0
36- self.surf = surface.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"],
37- exact=not args.interp_geometry)
38- self.lam = -case["beta"]
41+ self.surf = pb.gpu_surface(case, args)
3942 t0 = time.time()
40- self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, self.lam, ifinout=0,
43+ self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, case["lam"], ifinout=pb.IFINOUT,
4144 far_dtype=np.float64 if args.far64 else np.float32, verbose=False)
4245 cp.cuda.Device().synchronize()
4346 self.setup_s = time.time() - t0
4447 self.eps_gmres = args.eps_gmres
4548 self.options = dict(far_field="float64" if args.far64 else "float32",
46- geometry="interpolated" if args.interp_geometry else "exact")
49+ geometry="exact" if self.surf.exact is not None else "interpolated")
4750
4851 def solve(self, rhs):
4952 sigma, phi1, niter, rres = self.op.solve(rhs, eps_gmres=self.eps_gmres, verbose=False)
@@ -55,27 +58,33 @@ class GPUBackend:
5558 def eval_wall(self, state):
5659 return self.op.eval(None, *state, on_surface=True).get()
5760
61+ def to_numpy(self, q):
62+ return self.cp.asnumpy(q)
63+
5864 def sync(self):
5965 self.cp.cuda.Device().synchronize()
6066
67+ def timed_solve(self, rhs):
68+ return timed(self, self.solve, rhs)
69+
6170
6271 class Fmm3dbieBackend:
6372 """fmm3dbie on the CPU (helm_rpcomb_imp_solver), through ctypes."""
6473
6574 name = "fmm3dbie"
6675
67- def __init__(self, case, args):
76+ def __init__(self, pb, case, args):
6877 sys.path.insert(0, os.path.join(HERE, "reference"))
6978 import fmm3dbie_ctypes as fb
70- import rounded_cube as rc
7179 self.fb = fb
80+ self.xp = np
7281 self.log = os.path.join(args.outdir, "fortran.log")
7382 os.makedirs(args.outdir, exist_ok=True)
7483 open(self.log, "w").close()
7584 self.kernel_load_s = 0.0
7685 with fb.fortran_output_to(self.log):
77- self.surf = rc.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"])
78- self.k, self.lam = case["k"], -case["beta"]
86+ self.surf = pb.fmm3dbie_surface(case)
87+ self.k, self.lam, self.alpha, self.ifinout = case["k"], case["lam"], pb.ALPHA, pb.IFINOUT
7988 self.setup_s = 0.0 # the near quadrature is computed inside each solver call
8089 self.eps, self.eps_gmres = args.eps, args.eps_gmres
8190 self.options = dict(eps=args.eps, threads=os.environ.get("OMP_NUM_THREADS", "all"))
@@ -83,21 +92,28 @@ class Fmm3dbieBackend:
8392 def solve(self, rhs):
8493 with self.fb.fortran_output_to(self.log):
8594 sigma, sik, niter, errs, rres = self.fb.helm_rpcomb_imp_solver(
86- self.surf, self.k, pb.ALPHA, self.lam * np.ones(self.surf.npts), rhs, ifinout=0,
87- eps=self.eps, eps_gmres=self.eps_gmres)
95+ self.surf, self.k, self.alpha, self.lam * np.ones(self.surf.npts), rhs,
96+ ifinout=self.ifinout, eps=self.eps, eps_gmres=self.eps_gmres)
8897 return (sigma, sik), int(niter), float(rres)
8998
9099 def eval_points(self, state, pts):
91100 with self.fb.fortran_output_to(self.log):
92- return self.fb.helm_rpcomb_eval(self.surf, self.k, pb.ALPHA, *state, pts, eps=self.eps)
101+ return self.fb.helm_rpcomb_eval(self.surf, self.k, self.alpha, *state, pts, eps=self.eps)
93102
94103 def eval_wall(self, state):
95104 with self.fb.fortran_output_to(self.log):
96- return self.fb.helm_rpcomb_eval_boundary(self.surf, self.k, pb.ALPHA, *state, eps=self.eps)
105+ return self.fb.helm_rpcomb_eval_boundary(self.surf, self.k, self.alpha, *state,
106+ eps=self.eps, ifinout=self.ifinout)
107+
108+ def to_numpy(self, q):
109+ return q
97110
98111 def sync(self):
99112 pass
100113
114+ def timed_solve(self, rhs):
115+ return timed(self, self.solve, rhs)
116+
101117
102118 def timed(be, f, *a):
103119 t0 = time.time()
@@ -107,80 +123,66 @@ def timed(be, f, *a):
107123
108124
109125 def main():
126+ pre = argparse.ArgumentParser(add_help=False)
127+ pre.add_argument("--problem", choices=problems.PROBLEMS, default="room")
128+ pb = problems.get(pre.parse_known_args()[0].problem)
129+
110130 ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
131+ ap.add_argument("--problem", choices=problems.PROBLEMS, default="room")
111132 ap.add_argument("--backend", choices=["gpu", "fmm3dbie"], required=True)
112133 ap.add_argument("--case", choices=list(pb.CASES), default="small")
113- ap.add_argument("--k", type=float, help="override the case's wavenumber")
114- ap.add_argument("--beta", type=float, help="override the case's wall admittance")
115- ap.add_argument("--m", type=int, help="override the patches per face edge")
116- ap.add_argument("--norder", type=int, help="override the order per patch")
134+ pb.add_arguments(ap)
117135 ap.add_argument("--selftest", action="store_true",
118136 help="also solve a problem with a known solution (a second solve)")
119137 ap.add_argument("--eps-gmres", type=float, default=1e-6)
120138 ap.add_argument("--eps", type=float, default=1e-4, help="fmm3dbie quadrature and FMM tolerance")
121139 ap.add_argument("--far64", action="store_true", help="gpu: double precision far field")
122140 ap.add_argument("--interp-geometry", action="store_true",
123- help="gpu: interpolated patch geometry inside the integrals, as fmm3dbie")
124- ap.add_argument("--ngrid", type=int, default=101, help="grid size of the field slices (0 to skip)")
125- ap.add_argument("--outdir", help="default results/<case>-<backend>")
141+ help="gpu: interpolated patch geometry inside the integrals, as fmm3dbie "
142+ "(the cymbal always uses it)")
143+ ap.add_argument("--ngrid", type=int, default=pb.NGRID,
144+ help=f"grid size of the field slices (default {pb.NGRID}; 0 to skip)")
145+ ap.add_argument("--outdir", help="default results/<case>-<backend> (room), "
146+ "results/cymbal-<case>-<backend> (cymbal)")
126147 args = ap.parse_args()
127148
128- case = dict(pb.CASES[args.case])
129- for key in ["k", "beta", "m", "norder"]:
130- if getattr(args, key) is not None:
131- case[key] = getattr(args, key)
132- args.outdir = args.outdir or os.path.join(HERE, "results", f"{args.case}-{args.backend}")
133- k, lam = case["k"], -case["beta"]
149+ case = pb.make_case(args.case, args)
150+ args.outdir = args.outdir or os.path.join(HERE, "results", pb.results_name(args.case, args.backend))
151+ k, lam = case["k"], case["lam"]
134152
135- be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(case, args)
153+ be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(pb, case, args)
136154 surf = be.surf
137155 xs, nrm = np.asarray(surf.srcvals)[0:3], np.asarray(surf.srcvals)[9:12]
138- wts = np.asarray(surf.wts)
139- print(f"[{args.backend}] case {args.case}: k = {k}, beta = {case['beta']}, "
140- f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points, "
141- f"room width {2 * pb.A * k / (2 * np.pi):.2f} wavelengths")
156+ print(f"[{args.backend}] {pb.NAME} {args.case}: {pb.describe(case)}; "
157+ f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points")
142158 if args.backend == "gpu":
143159 print(f" kernel load {be.kernel_load_s:.2f} s, operator setup {be.setup_s:.2f} s")
144160
145- summary = dict(backend=args.backend, case=args.case, **case, npts=int(surf.npts),
146- eps_gmres=args.eps_gmres, options=be.options, machine=pb.machine_info(),
161+ summary = dict(problem=pb.NAME, backend=args.backend, case=args.case,
162+ **{key: val for key, val in case.items()}, npts=int(surf.npts),
163+ eps_gmres=args.eps_gmres, options=be.options, machine=machine_info(),
147164 kernel_load_s=be.kernel_load_s, setup_s=be.setup_s)
148165
149166 if args.selftest:
150- (state, niter, rres), t = timed(be, be.solve, pb.selftest_rhs(xs, nrm, k, lam))
151- xt = pb.selftest_points()
167+ x0s, ws = pb.selftest_sources(case)
168+ (state, niter, rres), t = be.timed_solve(source_sum(xs, x0s, ws, k, nrm=nrm, lam=lam, xp=be.xp))
169+ xt = pb.selftest_points(case)
152170 ut = be.eval_points(state, xt)
153171 ub = be.eval_wall(state)
154- uex, ubex = pb.green(xt, pb.EXTERIOR_SOURCE, k), pb.green(xs, pb.EXTERIOR_SOURCE, k)
172+ uex, ubex = source_sum(xt, x0s, ws, k, xp=be.xp), source_sum(xs, x0s, ws, k, xp=be.xp)
155173 err = float(np.linalg.norm(ut - uex) / np.linalg.norm(uex))
156174 err_wall = float(np.linalg.norm(ub - ubex) / np.linalg.norm(ubex))
157- print(f" self-test: {niter} iterations, {t:.2f} s; relative error {err:.2e} at 50 interior "
158- f"points, {err_wall:.2e} on the walls")
175+ print(f" self-test: {niter} iterations, {t:.2f} s; relative error {err:.2e} at "
176+ f"{xt.shape[1]} points off the surface, {err_wall:.2e} on the surface")
159177 summary.update(selftest_err=err, selftest_err_wall=err_wall, selftest_solve_s=t)
160178
161- (state, niter, rres), t_solve = timed(be, be.solve, pb.room_rhs(xs, nrm, k, lam))
162- print(f" room problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s")
163- v0 = complex(be.eval_points(state, pb.SOURCE[:, None])[0])
164- u_wall = pb.green(xs, pb.SOURCE, k) + be.eval_wall(state)
165- p_src, p_wall = pb.energy_balance(k, lam, v0, wts, u_wall)
166- mismatch = abs(p_src - p_wall) / p_src
167- print(f" energy balance: source {p_src:.8f}, walls {p_wall:.8f}, relative mismatch {mismatch:.2e}")
168- summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve,
169- p_src=p_src, p_wall=p_wall, energy_mismatch=mismatch)
170-
171- arrays = dict(x0=pb.SOURCE, k=k, lam=lam, a=pb.A, p=pb.P, xs=xs, wts=wts, u_wall=u_wall)
172- if args.ngrid > 0:
173- g, sl = pb.slice_points(args.ngrid)
174- t0 = time.time()
175- for name, (pts, mask, shape) in sl.items():
176- u = np.full(mask.size, np.nan + 1j * np.nan)
177- u[mask] = pb.green(pts[:, mask], pb.SOURCE, k) + be.eval_points(state, pts[:, mask])
178- arrays[f"u_{name}"] = u.reshape(shape)
179- be.sync()
180- summary["slices_s"] = time.time() - t0
181- arrays["grid"] = g
182- print(f" field on three {args.ngrid} x {args.ngrid} slices: {summary['slices_s']:.2f} s")
183- pb.save(args.outdir, summary, **arrays)
179+ (state, niter, rres), t_solve = be.timed_solve(-impedance_data(xs, nrm, pb.source(case), k, lam))
180+ print(f" main problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s")
181+ summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve)
182+
183+ arrays = {}
184+ pb.report(be, state, case, args, summary, arrays)
185+ save(args.outdir, summary, **arrays)
184186 print(f" saved {args.outdir}")
185187
186188
tests/geometry_test.pymodified+1−1View file
@@ -13,7 +13,7 @@ sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
1313 sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
1414
1515 from fmm3dbie_gpu import surface as S # noqa: E402
16-from problem import green, green_dn # noqa: E402
16+from problems.common import green, green_dn # noqa: E402
1717 from test_quadrature import kernels6 # noqa: E402
1818
1919 a_, p_, m_, norder = 1.0, 8, 4, 6