Coarse Solvers for Exascale Solution of Poisson Problems

Abstract: We present a two-level Schwarz method as an alternative to Algebraic Multigrid method(AMG) used as the last level (coarse) solver of the p-multigrid pMG preconditioner for pressure Poisson equation resulting from Spectral/Finite element descretization of incompressible Navier-Stokes equation. Proposed Schwarz method consits of a local problem in the original pMG coarse space and a global coarse problem. Main contribution of the paper is a novel, structured and a non-nested coarse space for the global coarse problem. Structured nature of the proposed global coarse space enable communication-free interpolation between the original p-multgrid coarse space and the global coarse problem. We demonstrate the effectiveness of the proposed method compared to the state of the art AMG solver BoomerAMG by a series of experiments performed using Nek5000/RS, a suite of highly scalable incompressible Navier-Stokes solvers, on Summit/Frontier supercomputers at Oak Ridge Leadership Computing Facility.

• The paper introduces a two-level Schwarz method with a structured, nonnested coarse space to overcome coarse grid bottlenecks in exascale Poisson solvers.
• It combines local overlapping subdomain solves and a communication-free global coarse solve, yielding 1.6x–2.7x faster coarse-grid performance compared to AMG.
• Scalability tests on Summit and Frontier demonstrate improved parallel efficiency and reduced time-to-solution for large-scale Navier-Stokes simulations.

Two-Level Schwarz Methods for Parallel Exascale Poisson Solvers

Background and Motivation

The paper addresses the computational bottleneck in p-multigrid (pMG) solvers for elliptic PDEs, notably the pressure Poisson equation from spectral/finite element discretizations of the incompressible Navier-Stokes equations. On exascale architectures (e.g., GPU-based supercomputers), the bottleneck emerges at the coarse grid level (lowest polynomial order, $N=1$), where the problem size reduces to the number of elements but the resulting system matrix $A$ is dense and communication-intensive. Standard approaches, such as Algebraic Multigrid (AMG) (e.g., BoomerAMG), exhibit sub-optimal scaling—especially as process count $P$ increases, limiting achievable strong scaling.

The motivation is driven by production-scale applications (e.g., full-core reactor simulations with $E\sim10^8$–$10^9$ elements, $n\sim10^{10}$–$10^{11}$ dofs), where traditional AMG-based coarse solves dominate overall runtime and restrict scalability. The coarse grid solve can account for 45% of flow simulation time at peak Summit runs.

Methodology

Two-Level Overlapping Schwarz Scheme

The core contribution is the introduction of a two-level Schwarz method as an alternative coarse grid solver within pMG. The Schwarz approach decomposes the problem into:

• Local subdomain solves: Each MPI process handles a local problem (~$E/P$ elements), extended by overlap of width $\delta = \mathcal{O}(h)$.
• Global coarse problem: A reduced coarse space with only a few dofs per process (typically 4–10), designed via structured, nonnested interpolation.

The additive Schwarz preconditioner is formalized as:

$$M^{-1} = \sum_{p=1}^{P} R_p A_p^{-1} R_p^T + J A_r^{-1} J^T$$

where $A$0 handles restriction/prolongation to overlapping local domains, and $A$1 interpolates between the original coarse mesh and the structured coarse problem.

Nonnested Structured Coarse Space

The global reduced coarse space is constructed as a tensor product of piecewise linear interpolants spanning a rectangular box covering the domain. This structured space enables communication-free interpolation/prolongation—both $A$2 and $A$3 can be computed locally, avoiding global communication in interpolation.

The nonnested design decouples the reduced space from the fine mesh. Any geometric mismatches (holes/inclusions, etc.) are handled via matrix assembly, ensuring robustness to irregular domains and boundary conditions.

The global coarse problem $A$4 is solved exactly using the $A$5 algorithm (fast parallel direct solver), with communication complexity near-optimal: $A$6 for small systems, where $A$7 is the network latency.

Implementation and Numerical Experiments

Experiments were performed on Summit and Frontier supercomputers, using the highly scalable Nek5000 (CPU) and NekRS (GPU/CPU hybrid) solvers. The coarse grid solvers were implemented using:

• CHOLMOD for local exact Cholesky solves with AMD ordering,
• $A$8 factorization for global coarse solves,
• Hypre/BoomerAMG for comparison.

Four production Navier-Stokes cases were considered:

• T-junction (E=62176)
• 146-pebble bed (E=62138)
• 45,000-pebble bed (E=13,032,440)
• 352,625-pebble annular bed (E=98,782,067, $A$9B)

Detailed timings, iteration counts, and strong scaling studies were conducted.

Numerical Results

Iteration Counts and Solver Timings

The two-level Schwarz method exhibits comparable iteration counts per timestep to AMG and $P$0 (Table 1, Fig. 4-5). Notably:

• Combination of local Schwarz solves and structured coarse space is crucial; neither alone produces competitive results.
• Slightly higher GMRES iterations with Schwarz vs. AMG are offset by lower coarse solve timings at high process counts.

Strong Scaling and Parallel Efficiency

Strong scaling studies (Fig. 6, Table 1) demonstrate superior parallel efficiency for the two-level Schwarz solver, especially as $P$1 increases; time-to-solution decreases faster than AMG at fixed parallel efficiency thresholds.

Speedup factors for coarse-grid solve times were observed (1.6x–2.7x faster than AMG for largest meshes and process counts), translating to overall Navier-Stokes simulation time reductions.

Communication-free interpolation enabled by the structured coarse space is pivotal at large scale, as it eliminates bottlenecks associated with AMG multilevel communication.

Complexity Analysis

Theoretical and measured complexity estimates align closely. For pMG preconditioner systems at strong scaling limit ($P$2), Schwarz's local and global solves maintain nearly constant cost per process, while AMG's multilevel communications (with $P$3 levels) drive up latency-bound time per solve.

Implications and Future Directions

Practical Implications

The structured, nonnested coarse space in two-level Schwarz methods provides a scalable alternative to AMG for exascale elliptic problems—offering superior parallel efficiency and reduced time-to-solution as systems and process counts grow. The communication-minimal design is suited for next-generation heterogeneous architectures (GPU-dominated HPC systems).

Integration in production CFD codes (Nek5000/NekRS) demonstrates practical competitiveness and robustness, even for complex geometries and domain inclusions.

Theoretical Implications

The approach bypasses the limitations of nested coarse spaces and multilevel communication in AMG, opening avenues for scalable domain decomposition methods. The theoretical bounds on the condition number with nonnested Schwarz spaces suggest comparable asymptotic convergence properties.

Future Research

Several avenues for performance optimization are suggested:

• Switching local/global solves to single precision to halve communication costs.
• Reducing fill-in and nonzeros in Cholesky and $P$4 factors via improved sparse reorderings.
• Using inexact solvers (incomplete Cholesky, iterative CG with preconditioning) for local/global solves, subject to maintaining acceptable outer iteration counts.

Further exploration of mesh coarsening and local refinement strategies in the reduced space, especially for highly irregular domains, is warranted.

Conclusion

This paper introduces a two-level Schwarz method with a novel structured, nonnested coarse space as an exascale-compatible coarse solver for pMG elliptic PDEs. Empirical and theoretical results demonstrate its competitiveness with AMG in terms of iteration counts, parallel efficiency, and overall time-to-solution, particularly at high process counts and large mesh sizes. The method's communication-minimal properties suit it well for current and future heterogeneous supercomputing platforms. Further optimizations (precision reduction, sparse reordering, inexact solves) are feasible and comprise avenues for subsequent research, with broad implications for scalable PDE solvers in computational science and engineering.

Reference: "Coarse Solvers for Exascale Solution of Poisson Problems" (2606.20496)