fmm3dbie on the GPU: a proof of concept#

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.

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.

A second showcase, a listener near a rigid cymbal in free space, is described in The cymbal problem. It uses the exterior mode of the same solver, and both backends run it.

A single precision port of the same solver to WebGPU, which runs in a browser without CUDA, is in fmm3dbie-webgpu-poc.

What is being compared#

Both backends solve for the correction v = u − G, where G is the free-space field of the source. The representation is v = S_k[σ] + iα D_k S_{i|k|}[σ], and the equation is fmm3dbie's discretized impedance equation. Both use the same 6m² quadrilateral patches (m × m on each face of the cube) of order p, with (p+1)² Gauss–Legendre nodes each, so the unknowns are the same.

• fmm3dbie runs fmm3dbie itself through ctypes (reference/), with the quadrature and FMM tolerance --eps (default 1e-4).
• gpu runs fmm3dbie_gpu/, a CuPy program with CUDA kernels compiled at run time. It does not use fmm3dbie.

The GPU code differs from fmm3dbie in four ways, which also explain its speed and accuracy.

1. Self-patch quadrature. fmm3dbie uses precomputed generalized Gaussian rules. We use polar coordinates on the four triangles fanned from the target to the patch edges, with a sinh substitution along each edge. This handles targets close to an edge, needs no tables, and computes every (target, patch) pair in one thread block.
2. Near-patch quadrature. We use a quadtree on the reference square, refined by distance and by wavelength, with an 8 × 8 Gauss–Legendre rule per leaf. This is fmm3dbie's distance-based strategy rather than its adaptive one.
3. Far field. Patches that are not near a target are summed directly, O(N²), rather than with an FMM. Near patches are excluded from this sum, so it has no large cancelling terms and runs in single precision, with partial sums accumulated in double precision.
4. Geometry. For the rounded cube we evaluate the exact surface inside the integrals. fmm3dbie uses polynomial interpolants of the node positions and tangents. The interpolants are slightly inconsistent with each other, and the double-layer kernel amplifies this at nodes next to patch edges, by up to 0.6% of the local integral in our tests. With the exact geometry the GPU solution is more accurate than fmm3dbie's on the same nodes. --interp-geometry switches to the interpolated model for a like-for-like comparison.

Quick start#

The fmm3dbie backend needs gfortran, make and OpenBLAS; the GPU backend needs an NVIDIA GPU with a CUDA 12 driver. Both need Python 3.9 or later with NumPy, and the plots need Matplotlib.

scripts/build_fmm3dbie.sh                 # clones and builds ~/src/fmm3dbie (or set FMM3DBIE_LIB)
scripts/setup_gpu.sh                      # CuPy and CUDA libraries as pip wheels, if no CUDA toolkit
source scripts/gpuenv.sh

python3 run.py --backend gpu --case small --selftest
python3 run.py --backend fmm3dbie --case small --selftest
python3 compare.py
python3 plot.py results/small-gpu

python3 run.py --problem cymbal --backend gpu --case large --selftest --reciprocity
python3 run.py --problem cymbal --backend fmm3dbie --case large --selftest --reciprocity
python3 plot.py results/cymbal-large-gpu

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).

The cases are

and --k, --beta, --m and --norder override them. Run python3 run.py --help for the other options.

What each run reports#

• Timings. For the GPU, setup is the near quadrature for the operator and solve is GMRES. The one-time loading of the compiled CUDA kernels is reported separately and excluded. fmm3dbie computes its near quadrature inside the solver call, so its setup is zero and its solve includes it. Compare setup + solve.
• Self-test (--selftest). This solves a second problem whose solution is known, the field of a source outside the room, and reports the relative error at 50 interior points and on the walls. This is the accuracy measure. Note that the GMRES tolerance (--eps-gmres, default 1e-6) bounds it from below.
• Energy balance. The power radiated by the source, k/(4π) + Im v(x₀), must equal the power absorbed by the walls, −k Re(λ)∫|u|² dS. Their relative mismatch is a check that needs no exact solution.

Results#

Both backends were run with run.py, at the default settings, on the same workstation. That machine has an Intel Xeon w5-3435X (16 cores, 32 threads, all used by fmm3dbie) and an NVIDIA RTX 6000 Ada (48 GB). fmm3dbie was built at commit 94cb732 and run with eps = 1e-4. The GMRES tolerance was 1e-6 for both. The table is results/comparison.md as written by compare.py.

The two codes need the same number of GMRES iterations, which is consistent with their discretizing the same operator. Their wall fields agree to about the accuracy of each. The speedup grows with the problem size, because the fixed costs of the GPU code become less important. On the large case, 85% of the GPU time is GMRES. Most of that is presumably the O(N²) far-field sum, which is the part an FMM would reduce. The CUDA kernels are compiled or loaded once per process, taking about 2 s, and this is excluded from the GPU times. fmm3dbie's first solve in a process can be slower than its second: on the medium case the self-test solve took 81 s and the room solve 18 s. We report the room solve.

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).

The cymbal problem#

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.

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.

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.

--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.

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.

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.

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.

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.

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.

Files#

• run.py, compare.py, plot.py: run one backend, compare the results, and plot a field.
• 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.
• fmm3dbie_gpu/surface.py: patch discretization and the rounded cube. It is pure NumPy.
• fmm3dbie_gpu/cymbal.py: the cymbal's shape and its patch mesh, also pure NumPy. Both backends build the cymbal from it.
• fmm3dbie_gpu/cuda_ops.py: the CUDA kernels: near pairs, quadtrees, near weights, far-field sum and near-field product.
• fmm3dbie_gpu/solver.py: the operator, evaluation, and GMRES, following fmm3dbie's zgmres_guru.
• 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.
• tests/test_quadrature.py: compares the GPU near weights with brute-force integration. For example, python3 tests/test_quadrature.py 8 40 m=10 tp=45 runs at order 8 and k = 40.
• tests/geometry_test.py: checks Green's identity at nodes next to patch edges with the interpolated, derivative-consistent and exact geometry models. This is the evidence for point 4 above.

Limitations#

The far field is a direct sum, so each product costs O(N²). At 48,600 points one product takes about 40 ms. Somewhere between 10⁵ and 2·10⁵ points a GPU FMM becomes necessary; FMM3D's Helmholtz scheme is the obvious model. The near weights would not change, since the FMM would only replace the dense sum.

The exact geometry exists only for the rounded cube. General surfaces use the interpolated geometry, and with it fmm3dbie's sensitivity at nodes next to patch edges.

The self-patch rule has been checked against brute-force integration for orders 6 and 8 and for k up to 40, but not for higher orders or more strongly curved patches.

Only this one solver is ported. The two codes are compared on one family of problems, and the relative speed will differ for other geometries, orders and tolerances.