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

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

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

Files#

• run.py, compare.py, plot.py: run one backend, compare the results, and plot a field.
• problem.py: the showcase problem, shared by both backends: cases, source positions, right-hand sides and checks.
• fmm3dbie_gpu/surface.py: patch discretization and the rounded cube. It is pure NumPy.
• 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.