fmm3dbie-webgpu-poc#

A limitation of the CUDA proof of concept for fmm3dbie's interior Helmholtz impedance solver is that it needs an NVIDIA GPU, CUDA and a Python environment, which makes it hard to try and hard to share. Here we present a port of the same solver to WebGPU, as a web page. It runs on whatever GPU the browser can reach, including integrated graphics, and solves the same showcase problem interactively.

Live page: https://fmm3dbie-webgpu-poc--proof-of-concept.vault1-sites.magland.org/

The page solves either of the CUDA version's two showcase problems. The first is a point source at a fixed frequency in a closed room with partially absorbing walls. The room is the superellipsoid |x|^p + |y|^p + |z|^p = 1, and with time dependence e^{-iωt} the pressure amplitude satisfies Δu + k²u = −δ(x − x₀) inside and ∂u/∂n − ikβu = 0 on the walls, where β > 0 is the wall admittance. The page lets the user set k, β, the discretization, the room shape and the source position, and shows the field on three planes through the source and on the walls, in 3D, animated over one period.

The second is a listener near a stationary, rigid cymbal in free space. The source is placed at the listener, so that by reciprocity the field u(x) is the pressure at the listener due to a unit source at x. The cymbal is 40 cm across, with a raised bell; since fmm3dbie's equations need a closed surface, it is a thin closed body, 20 mm thick at the center and 10 mm (1/40 of the diameter) at its rounded rim. The equation is the exterior mode of the same impedance equation, with λ = 0. The page lets the user set the frequency, the patch size and order, the rim thickness and the listener's position, and shows the field on three planes through a 1 m box around the cymbal, with the level relative to the free field. ?problem=cymbal opens the page on this problem.

Method#

The equation, representation and discretization are those of the CUDA version: the correction v = u − G is represented as S_k[σ] + iD_k S_{i|k|}[σ], on 6m² curved quadrilateral patches with (order + 1)² Gauss–Legendre nodes each. The near part of each layer potential is computed once with adaptive quadrature (polar coordinates with a sinh substitution on the target's own patch, a quadtree on neighbouring patches) and stored as weights on the GPU. The far part is a direct O(N²) sum over oversampled nodes, recomputed at every GMRES iteration. GMRES runs on the CPU in double precision; only the operator is applied on the GPU.

Unlike the CUDA version, which uses double precision for the near quadrature, everything on the GPU is single precision, since WebGPU has no double precision type. Three changes make this sufficient:

1. For a target on the wall, the vector from the target to a nearby quadrature node is formed in the coordinates of the patch's parametrization rather than as a difference of two rounded positions. It then keeps full relative precision however short it is, which the double-layer kernels need, since their numerators (y − x)·n are O(r²). For the room this is the superellipsoid's analytic parametrization. For the cymbal, and for any other surface, it is the Legendre interpolant of the node data (the geometry fmm3dbie uses), and the vector comes from divided differences of the Legendre polynomials, P_n(u) − P_n(u_t) = (u − u_t) δP_n, with δP_n from its own three-term recurrence. For a target on a neighboring patch, the CPU finds its preimage in the patch's (extended) interpolant by Gauss–Newton, and passes the small remaining gap along with it.
2. The weakly singular difference D′k − D′{i|k|} is evaluated from its power series when kr < 1, where the direct formula subtracts two nearly equal O(1/r³) terms.
3. sin and cos are evaluated with our own range reduction and polynomials. WGSL only requires them to be accurate to 2⁻¹¹ in absolute terms, and some implementations (SwiftShader among them) are no better than that, which limited the near weights to about 1e-4.

Accuracy#

npm run test:gpu compares the near weights with a brute-force double precision reference, checks Green's identity through the full operator, and solves a problem with a known solution (the field of a source outside the room). On 54 patches of order 6 at k = 4 with the SwiftShader software adapter, the near weights agree with the reference to 7e-6 or better, the Green's identity residual is 3.4e-6 (3.6e-6 in double precision on the CPU for the same discretization), and the self-test error is 1.6e-6 in the interior and 9.0e-6 on the walls. For comparison, the CUDA version reports 2.2e-6 on its small case.

The suite then runs the same checks on a coarse cymbal (192 patches of order 5, with a 20 mm rim). On an RTX 6000 Ada the near weights agree with the double precision reference to 1.7e-5 or better, and the exterior self-test error, with sources spread inside the cymbal, is 2.0e-4 around the cymbal and 6.9e-4 on it, which is the discretization error of that coarse surface.

test.html?bench=small|medium|large runs the named cymbal cases of the CUDA version, on the same surface and unknowns. On the same RTX 6000 Ada (Chrome under Xvfb, Vulkan), with GMRES tolerance 1e-6:

The two codes agree to about 1e-6, so on these cases single precision costs no accuracy that the discretization would otherwise give. The WebGPU code is two to three times slower than CUDA on the same GPU. fmm3dbie on the workstation's 32 CPU threads takes 27 to 81 s for these solves.

For sharper rooms the discretization, not the precision, limits the accuracy. At p = 30 the same surface gives a Green's identity residual of 1.1e-4 on the GPU and 1.1e-4 in double precision; refining to 96 patches of order 8 brings the self-test error down to 1.7e-5.

Running it#

npm install
npm run dev          # the app at http://localhost:5173
npm run test:node    # CPU checks of the surface and the reference quadrature
npm run test:gpu     # the GPU validation suite in headless Chrome
npm run smoke        # build, then solve the quick case in the app, headless
node scripts/smoke.mjs quick cymbal   # the same for the cymbal, after a build

The page accepts ?problem=cymbal, ?case=<name> and ?autosolve, which solves at once. After a solve, Export GIF writes the three planes side by side as one GIF, animated over one period for the pressure (24 frames) or a single image for the level.

The headless tests use a hardware WebGPU adapter when Chrome finds one and SwiftShader otherwise. On a Linux machine without a display, headless Chrome may only find SwiftShader even when there is a GPU; a windowed Chrome under Xvfb (xvfb-run -a google-chrome --ignore-gpu-blocklist --enable-unsafe-webgpu --enable-features=Vulkan ...) found the NVIDIA adapter in our tests. Note that in headless Chrome, presenting a frame to a WebGPU canvas makes later asynchronous WebGPU calls fail ("a valid external Instance reference no longer exists"), so the smoke test opens the page with ?noroom, which skips the 3D view. We have not seen this in a browser with a window.

Limitations#

The near weights take 6 × (order + 1)² × 8 bytes per pair of a node and a near patch. At order 8 this is several hundred MB for the room's medium case (17,496 points) and more than most browsers allow in one buffer for its large case (48,600 points). The thin cymbal has more near pairs per node, and its small, medium and large cases need 640, 730 and 780 MB. Splitting the weights across buffers, or recomputing them, would remove this limit.

The far field is a direct sum, so each operator application costs O(N²); a GPU FMM would replace it for larger problems without changing the near weights.

The interpolated geometry with divided differences works for any surface of quadrilateral patches, but only the rounded cube and the cymbal are wired into the page. The far-field energy check for the cymbal runs on the CPU and is skipped for the largest cases.

Timings on the software adapter are not representative, and we have not yet measured the named cases on hardware GPUs.