Efficient Domain Decomposition for the Helmholtz Equation on GPUs

Abstract: The Helmholtz equation governs wave propagation in acoustics, electromagnetics, and seismology, but its indefinite nature makes it difficult to solve with iterative methods. Domain decomposition methods are a natural fit for massively parallel architectures, yet mapping efficient Helmholtz solvers onto modern GPUs remains a challenge. We address both with two key contributions: (1) a block-level domain decomposition scheme, in which each subdomain is assigned to a single thread block and all solves run concurrently in a single kernel launch, and (2) WaveHoltz as the subdomain solver. WaveHoltz is a fixed-point iteration that is uniquely well-suited to the GPU execution model due to its minimal memory footprint and no reduction operations. Together, these eliminate device-level synchronizations and replace global memory traffic with shared memory and register-level operations, keeping subdomain data largely resident in L1 and L2 cache. We explore two threading strategies: one degree of freedom per thread for small subdomains, and multiple degrees of freedom per thread for larger ones. Benchmarks of our CUDA based implementation on a NVIDIA A100 show that WaveHoltz achieves 2x-25x speedup over MINRES, with the advantage growing with subdomain size. Crucially, evaluating the subdomain solver in single rather than double precision yields an additional 2x-10x speedup--a benefit largely unattainable by MINRES due to loss of Krylov vector orthogonality under reduced precision.

• The paper introduces an innovative block-level domain decomposition that maps each subdomain to a GPU thread block using the WaveHoltz iterative method.
• It leverages spectral element discretization and optimized GPU memory hierarchies to achieve significant speedups, with benchmarks showing up to 25x faster performance over MINRES.
• Mixed precision, minimal synchronization, and efficient cache usage enable scalable 2D and 3D Helmholtz simulations on modern GPU architectures.

Efficient Domain Decomposition for the Helmholtz Equation on GPUs

Introduction

The Helmholtz equation, prevalent in modeling wave propagation across acoustics, electromagnetics, and seismology, is notably challenging for iterative solvers due to its indefinite nature. While domain decomposition methods (DDMs) provide natural scalability for parallel architectures, efficient GPU-based Helmholtz solvers remain nontrivial. This paper introduces a block-level domain decomposition approach wherein each subdomain is solved using the WaveHoltz iterative algorithm, mapped efficiently onto GPU thread blocks. The design eliminates device-level synchronization and leverages shared memory, maximizing L1/L2 cache residency and parallelism. Benchmarks demonstrate strong performance and scalability, particularly when using mixed precision.

Algorithmic Framework

Domain Decomposition and Transmission Conditions

The paper adopts a non-overlapping domain decomposition paradigm, partitioning the computational domain into subdomains assigned to individual thread blocks. The transmission condition between subdomains is set as $B_{e,r} = \partial_n - ik$, a zero-order approximation, which avoids complications with non-local Dirichlet-to-Neumann operators while maintaining algorithmic simplicity on GPUs.

Spectral Element Discretization

Each subdomain is discretized using spectral element methods (SEM), exploiting tensor product bases and GLL quadrature. The resulting system is solved in matrix-free fashion, further optimizing GPU performance by reducing memory overhead and utilizing basis structure for efficient evaluation.

WaveHoltz Iteration

WaveHoltz is employed as the subdomain solver—a fixed-point iteration designed for GPU execution. Its salient features are:

• Minimal memory footprint: Requires only one auxiliary vector compared to five in MINRES.
• No reduction operations: Avoids thread synchronizations typical in Krylov methods.
• Explicit time-stepping: Leveraging mass-lumped diagonal matrices, WaveHoltz is executed as a modified leapfrog scheme, ensuring arithmetic intensity and enabling efficient cache usage.

Benchmarks indicate WaveHoltz converges within 5–20 iterations for typical subdomain sizes, facilitating rapid fixed-point updates in domain decomposition iterations.

GPU Kernel Strategies

Two threading strategies are explored:

One DOF per Thread: Each thread in a block handles a single degree of freedom, optimal for subdomains with ≤1024 DOFs. All computation, except matrix-vector products, is confined to registers, resulting in minimal global memory traffic.

T DOFs per Thread: Each thread manages T DOFs (T = 2 or 4), extending support to larger subdomains at the cost of increased cache pressure. Here, partial results are serialized but efficiently buffered, balancing throughput and memory access.

The domain decomposition cycle is fused into a single kernel launch per subdomain, thus minimizing global memory access latency (300 cycles) relative to shared memory (20–30 cycles). This architectural mapping is validated by profiling cache traffic, which shows substantial confinement to L1/L2 cache.

Numerical Benchmarks

Two-Dimensional Performance

• The DDM cycle (WaveHoltz iteration) is 2x–4x the cost of operator evaluation for low-order SEM (P=1), and up to 200 MDOFs/sec throughput for higher orders (P=7).
• Optimal performance is achieved with smaller subdomains and lower-order discretizations. Cache hit rates are >98% in ideal configurations.
• The use of single precision delivers 1.5x–3.5x speedup on NVIDIA A100, and even higher on consumer GPUs where double precision is emulated.

Three-Dimensional Performance

• Throughput is halved due to cubic scaling of DOFs per element, but the approach retains practical utility for low-order problems.
• The DDM cycle is 10x slower than operator evaluation for P=1, and 30x–40x for higher orders.

Comparison with MINRES

• WaveHoltz achieves 2x–25x speedup over MINRES in single precision, with the gap increasing as subdomain size grows.
• MINRES demonstrates minimal benefit from mixed precision, as Krylov vector orthogonality deteriorates under reduced floating point accuracy.
• In double precision, WaveHoltz remains consistently faster (up to 9x), except for small subdomains in three dimensions.
• These results underscore the suitability of WaveHoltz for GPU architectures, particularly when memory hierarchy and precision are leveraged efficiently.

Implications and Future Directions

The proposed methodology exhibits strong practical and theoretical advantages:

• Practical performance: The single-kernel-per-block approach maximizes throughput, minimizes synchronization, and exploits GPU memory hierarchy optimally.
• Robustness to mixed precision: WaveHoltz's stability in single precision enables further acceleration, a property not shared by Krylov subspace methods on indefinite systems.
• Extensibility: Higher order transmission conditions (e.g., optimized Schwarz or ABCs) can be incorporated with minimal algorithmic change. The use of coarse-space corrections, such as plane wave ansatz or eigenvector-based augmentation, is a feasible extension to address adverse scaling in three dimensions.
• Multi-GPU potential: The independence of subdomain solves makes distribution across GPUs trivial, requiring only interface communication.

Theoretical implications include the improved algorithmic mapping of fixed-point DDM cycles onto massively parallel hardware, opening avenues for scalable, high-frequency wave simulations in large domains. The variable tolerance strategy, wherein initial solves are relaxed and tightened as iterations advance, promises further runtime reduction and warrants investigation.

Conclusion

This work establishes an efficient domain decomposition procedure for the Helmholtz equation on GPUs. By allocating each subdomain to a single thread block and employing the WaveHoltz iteration, synchronization and memory traffic are minimized, resulting in substantial performance gains over traditional Krylov solvers, especially in reduced precision. The method is robust and extensible, with clear pathways for algorithmic enhancement and deployment across multiple GPUs. The findings inform both practical system design for high-performance wave simulation and theoretical strategies for indefinite PDEs in parallel environments.

(2606.21695)