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

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 superellipsoid's analytic 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²).
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.

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

The headless tests use a hardware WebGPU adapter when Chrome finds one and SwiftShader otherwise. 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 medium case (17,496 points) and more than most browsers allow in one buffer for the large case (48,600 points). 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.

Only the rounded-cube family of rooms is available, because the single precision near quadrature relies on its exact parametrization. A general surface would need the interpolated geometry and a different way of forming short displacements accurately.

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