Coarse Solvers for Exascale Solution of Poisson Problems
arXiv:2606.20496v1 [math.NA] 18 Jun 2026
Thilina Ratnayakaa , Paul Fischerb and Luke Olsonc ARTICLE INFO
ABSTRACT
Keywords: Schwarz Multigrid Exascale
We present a two level Schwarz method as an alternative to Algebraic Multigrid method (AMG) used as the last level (coarse) solver of the 𝑝-multigrid 𝑝MG preconditioner for pressure Poission equation resulting from Spectral/Finite element descretization of incompressible Navier-Stokes eqaution. Proposed Schwarz method consits of a local problem in the original 𝑝MG coarse space and a global coarse problem. Main contribution of the paper is a novel, structured and a nonnested coarse space for the global coarse problem. Structured nature of the proposed global coarse space enable communication-free interpolation between the original 𝑝-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.
1. Introduction 𝑝-multigrid (𝑝MG) is an effective preconditioning strategy for the scalable solution of elliptic problems that are discretized by high-order finite-element (FEM) or spectral-element methods (SEM) [22, 3, 21, 28]. For example, using a local polynomial order 𝑁 = 7 with the SEM, a multigrid smooth/restrict/prolongate schedule of polynomial orders 𝑁 = 7 → 5 → 3 → 1 is effective and reduces the number of degrees-of-freedom (dofs) per element in 3D from 343 → 125 → 27 → 1. Through local (fast) tensor-product operations that require only 𝑂(𝑁 3 ) memory references and 𝑂(𝑁 4 ) work per element, the method is highly efficient [9]. While this communication-minimal relaxation process leads to a two-order-of-magnitude reduction in problem size, the resulting coarse grid problem for 𝑁 = 1, denoted 𝐴𝑐 𝑢𝑐 = 𝑏𝑐 , is communication intensive because 𝐴−1 𝑐 is completely dense. This is expected as the Green’s functions for the Poisson problem have slow decay; as a result, data in any entry of the distributed coarse input, 𝑏𝑐 , has an immediate impact on all values of the distributed solution, 𝑢𝑐 , which implies that the communication is all-to-all. In this paper, we consider development of coarse grid solvers tailored to exascale architectures (specifically, GPUbased platforms with high node counts). The target application is the Poisson problem that arises when simulating unsteady incompressible flow with spatial discretization based on the SEM [9, 27]. The SEM is a heterogeneous discretization that features locally structured degrees-of-freedom 𝑢𝑒𝑖𝑗𝑘 , (𝑖, 𝑗, 𝑘) ∈ [0, … 𝑁 − 1]3 , which are ordered lexicographically within each of 𝐸 elements, Ω𝑒 , 1 ≤ 𝑒 ≤ 𝐸. To accommodate complex geometries, the computational mesh comprising the elements is globally unstructured. The mesh restrictions are simply that each Ω𝑒 should be an ̂ ∶= [−1, 1]3 , that the elements are nonoverlapping,1 invertible (and well-conditioned) map of the reference element Ω and that their union covers the computational domain; that is, Ω = ∪𝑒 Ω𝑒 . For 𝑝MG [28] or multilevel Schwarz [18, 22] preconditioners, the heterogeneous nature of the SEM leads to a natural decomposition of structured (tensor-product) work on the local mesh, which accounts for the majority of the floating-point operations (flops), and communicationintensive operations on the unstructured element-vertex (i.e., 𝑁 = 1) mesh, referred to as the coarse grid. It is well known that the coarse grid solve2 is challenging in a parallel setting when the number of processes, 𝑃 , is large. It has been observed that the coarse grid solve is the only operation in the solution of elliptic partial differential equations (PDEs) whose cost increases with 𝑃 [17]. The other terms scale like (𝑛∕𝑃 )𝛾 for some 𝛾 ≥ 0, where 𝑛 is the total number of dofs. While terms with 𝛾 < 1 will impede ideal speed-up in the case of strong-scaling, where 𝑃 is increased for fixed 𝑛, it is only the coarse grid solve that impacts weak scaling, where 𝑃 is increased while the ∗ Thilina Ratnayaka
[email protected] (T. Ratnayaka); [email protected] (P. Fischer); [email protected] (L. Olson) ORCID (s): 1 Overlapping SEM meshes in the spirit of overset methods are, however, feasible [24, 26]. 2 We refer to the coarse grid solve throughout the text, but the coarse grid problem is rarely solved exactly because the level-𝑁 outer iteration
is embedded in a Krylov-subspace projection (KSP) such as conjugate gradients or GMRES. For the coarse grid, a simple relaxation or other error-reduction scheme often suffices.
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 1 of 14
Fast Coarse Solvers
number of dofs local to each process, 𝑛∕𝑃 , is fixed. As 𝑃 increases, the coarse grid solve invariably sets the limits on strong-scaling, which directly impacts time-to-solution in many applications. The importance of the coarse grid solve for exascale applications is illustrated in Min et al. [25], where it accounts for 45% of the flow simulation time when running on 𝑃 = 27 648 Nvidia V100 GPUs on the Summit supercomputer at the Oak Ridge Leadership Computing Facility (OLCF). In this case, the Navier-Stokes problem has 𝐸 = 98 million elements of order 𝑁 = 8 (leading to 𝑛 = 51 billion dofs) and the coarse solve in the 𝑝MG-preconditioned pressure Poisson problem uses a single 𝑉 -cycle of BoomerAMG [20]. For 𝑝MG on tensor-product elements, the coarse grid size is 𝑛𝑐 ≈ 𝐸. Current exa- and pre-exascale platforms are routinely enabling 𝐸 = 108 –109 , resulting in large “coarse” systems. In [25], the coarse grid overhead is indrectly reduced by adding more smoothing sweeps at the finer 𝑝MG levels, thereby reducing the number of outer GMRES iterations and consequent visits to the 𝑁 = 1 level of 𝑝MG. In this paper, we consider a direct approach to reducing communication for the 𝑁 = 1 coarse solve 𝐴𝑐 𝑢𝑐 = 𝑏𝑐 by introducing a low-communication two-level overlapping Schwarz smoother with a novel nonnested coarse space, as an alternative to AMG. The overlapping systems (size ≈ 𝐸∕𝑃 ) are highly localized and well suited to computation on the GPU, yielding full 𝑃 -fold parallelism with minimal off-device data transfer. The nonnested (reduced) coarse space requires only a few (e.g., 4–10) dofs per process and the corresponding reduced system, 𝐴𝑟 𝑢𝑟 = 𝑏𝑟 (size 𝑛𝑟 = 𝑂(𝑃 ) ≪ 𝐸), is solved using the communication-minimal 𝑋𝑋 𝑇 algorithm developed in [15, 31]. Rest of the paper is organized as follows. In Section 2, we review past coarse grid solve strategies and underscore the need for a novel approach in the context of GPU-based exascale platforms. In Section 3, we present the proposed two-level Schwarz method, including a coarser system, 𝐴𝑏 𝑢𝑏 = 𝑏𝑏 , which is based on a simple nonnested approximation space. We discuss the background and details of the nonnested coarse space in Section 4. Numerical results with both Nek5000 and NekRS using four different meshes are shown in Section 5. Conclusions and future directions are given in Section 6.3.
2. Survey of Coarse Grid Solvers Solving the coarse problem using the explicit inverse of 𝐴𝑐 becomes prohibitive in 3D for large 𝑛𝑐 due to the −1 storage cost of 𝐴−1 𝑐 . Even if 𝐴𝑐 is distributed across the number of processes 𝑃 , this method has a storage cost of 2 𝑛𝑐 ∕𝑃 . However, the fact that we can simply solve the system using 𝑢 = 𝐴−1 𝑐 𝑏 is quite appealing since all the dot products can be performed in parallel. Approach described in [1] try to exploit this feature by inverting the Cholesky factor 𝐿 of 𝐴𝑐 using a partitioned inverse approach. Any unit lower triangular matrix of size 𝑛𝑐 can be expressed as a ∏𝑛𝑐 𝐿𝑖 , 𝐿𝑖 is of the form 𝐼 + 𝑚𝑖 𝑒𝑇𝑖 where 𝑒𝑖 is the 𝑖𝑡ℎ unit vector and product of 𝑛𝑐 elementary matrices 𝐿𝑖 i.e., 𝐿 = 𝑖=1 𝑚𝑖 is a vector with first 𝑖 values being zeros. Then these 𝐿𝑖 elementary matrices are grouped together to form 𝑘 factors ∏ such that 𝐿 = 𝑘𝑖=1 𝑃𝑖 where each factor 𝑃𝑖 has the property that 𝑃𝑖−1 can be represented using the same space as 𝑃𝑖 . ∏ The solution 𝑢 = 𝐿−1 𝑏 = 1𝑖=𝑘 𝑃𝑖−1 𝑏 can be calculated by doing 𝑘 matrix vector multiplications. Communication for this approach scales as 𝑂((log2 𝑃 )2 ) and is not optimal. During time stepping of a computational fluid dynamics simulation, system 𝐴𝑐 𝑢𝑐 = 𝑏𝑐 is solved multiple times with different right hand sides for the same constant matrix 𝐴𝑐 . In time transient problems, successive right hand sides (RHS) often share enough information so a good initial guess based on previously generated Krylov subspaces can yield significant reductions in the global solve costs [14, 16]. Farhat and Chen [13] suggest using this approach to solve 𝐴𝑐 𝑢𝑐 = 𝑏𝑐 for successive RHS by using a slightly modified version of Conjugate Gradient (CG) iteration to make sure that the search directions generated for new RHS are 𝐴−orthogonal to previously generated Krylov subspace. When a solution 𝑢𝑐 for a new RHS 𝑏𝑐 has to be found, the initial guess is computed as the solution for a restricted versions of 𝐴𝑐 𝑢𝑐 = 𝑏𝑐 in the previously constructed Krylov subspace. The dominant cost for this method is this projection step to calculate the initial guess since the cost for subsequent CG iterations are negligible. Each of the CG iterations following this initial projection incurs a latency cost of 4𝛼 log2 𝑃 due to the reductions associated with CG. 𝑇 An alternative approach is to use the sparse direct factorization 𝐴−1 𝑐 = 𝑋𝑋 proposed in [31]. Main idea is to find a set of 𝑘, 𝐴𝑐 -conjugate sparse basis vectors 𝑋𝑘 = [𝑥1 𝑥2 … 𝑥𝑘 ] satisfying 𝑥𝑇𝑖 𝐴𝑐 𝑥𝑗 = 𝛿𝑖𝑗 , where 𝛿𝑖𝑗 is the Kronecker delta function. Then an approximation to 𝑥𝑐 , 𝑥̄ 𝑘 can be calculated by projecting 𝑏𝑐 in to this space i.e., 𝑥̄ 𝑘 = 𝑋𝑘 𝑋𝑘𝑇 𝑏𝑐 . If 𝑘 = 𝑛𝑐 , this approach yields the exact solution, 𝑥𝑐 = 𝑥̄ 𝑛 . By using a nested-disection ordering, the 𝑐 authors showed that the communication complexity for a system arising from a discrete Poisson problem in 𝑑 space (𝑑−1)∕𝑑 dimensions is 2(𝛼 ∗ + 𝛽 ∗ 𝐶𝑛𝑐 ) log2 𝑃 , where 𝐶 is a constant and (𝛼 ∗ , 𝛽 ∗ ) are interprocess latency and inverse bandwidth respectively. For sufficiently small systems this cost is ≈ 2𝛼 ∗ log2 𝑃 , which is near optimal. T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 2 of 14
Fast Coarse Solvers
3. Two-Level Coarse Solver In this section, we develop a low-communication preconditioner, 𝑀𝑐 , to approximate the coarse-level (𝑁 = 1) system matrix, 𝐴𝑐 . Given that the full pMG system is embedded in a Krylov subspace projection (KSP), an exact solution for the 𝑁=1 system is not required; a single sweep of a multilevel preconditioner is often sufficient. A potential consequence of an inexact solve, however, is an increase in the number of KSP iterations, so any economization must be weighed against the overall cost of the outer KSP solver, as discussed in the examples of Section 5. We base 𝑀𝑐 on the two-level additive Schwarz method (ASM) of Dryja and Widlund [10, 11, 12]. To simplify notation, we drop the subscript 𝑐 for the coarse problem and simply refer to the system as 𝐴𝑢 = 𝑏 and to the preconditioner as 𝑀, unless otherwise needed. In this context, 𝐴 is assumed to be a sparse 𝑛 × 𝑛 symmetric positive definite (SPD) matrix derived from an unstructured nodal finite element (FEM) discretization of the Poisson equation with homogenous Dirichlet boundary conditions. (In our case, the 𝐴𝑐 system is based on a Galerkin restriction of the originating order-𝑁 SEM operator.) Many features of this two-level method would work equally well for other PDEs, however, and for other nodal-based discretizations such as finite differences or finite volumes. We reiterate that our target problem size will involve only 𝑛∕𝑃 ≈ 103 –104 dofs per GPU and that this problem is thus intrinsically communication bound, despite the potentiality that 𝑛 might be as large as 109 . To set the stage for later developments, we introduce two representations for FEM basis functions on the domain, Ω, having boundary 𝜕Ω = 𝜕Ω𝑁 ∪ 𝜕Ω𝐷 . We assume that the boundary conditions for the PDE are homogeneous Neumann on 𝜕Ω𝑁 and homogeneous Dirichlet on 𝜕Ω𝐷 . Let 𝑋0ℎ =span{𝜙𝑗 (𝐱)} be the set of 𝐶 0 continuous basis functions satisfying 𝜙𝑗 (𝐱𝑖 ) = 𝛿𝑖𝑗 (the Kronecker delta function) that forms the FEM basis on nodal points 𝐱𝑖 . We assume that all basis functions vanish on 𝜕Ω𝐷 and that any trial solution 𝑢(𝐱) ∈ 𝑋0ℎ can be written as 𝑢(𝐱)
=
𝑛 ∑
𝜙𝑗 (𝐱) 𝑢𝑗 =
𝑗=1
𝑛𝑣 𝐸 ∑ ∑
𝑙𝑘𝑒 (𝐱) 𝑢𝑒𝑘 .
(1)
𝑒=1 𝑘=1
The first expression on the right is the nodal (or global) form involving unknown basis coefficients, 𝑢𝑗 , often referred to as global degrees-of-freedom (dofs). The second expression is the elemental (or local) form. We assume that Ω = ∪𝑒 Ω𝑒 represents a decomposition of Ω into 𝐸 nonoverlapping trianglular or quadrilateral elements (tetrahedral or hexahedral in 3D) of size 𝑂(ℎ), each having 𝑛𝑣 vertices, 𝐱𝑘𝑒 . For each element, there are 𝑛𝑣 local Lagrange interpolants, 𝑙𝑘𝑒 (𝐱), 𝑘 = 1, … , 𝑛𝑣 , which vanish outside of Ω𝑒 . There are a total of 𝑛𝑙 = 𝐸 ⋅ 𝑛𝑣 local basis coefficients, 𝑢𝑒𝑘 . To impose 𝐶 0 ′ ′ continuity on 𝑢(𝐱), we additionally have a constraint on these coefficients, 𝐱𝑘𝑒 = 𝐱𝑘𝑒 ′ ⟹ 𝑢𝑒𝑘 = 𝑢𝑒𝑘′ . The standard way to impose this constraint is to start with the global representation and to map the coefficients 𝑢𝑗 to their local counterparts, 𝐱𝑘𝑒 = 𝐱𝑗 ⟹ 𝑢𝑒𝑘 = 𝑢𝑗 . The transformation from global to local coefficients can be expressed as a matrix-vector product, 𝑢𝐿 = 𝑄𝑢, where 𝑄 is an 𝑛𝑙 × 𝑛 Boolean matrix, 𝑢𝐿 = {𝑢𝑒𝑘 } is the set of local unknowns, and 𝑢 = {𝑢𝑗 } is the set of global unknowns. Application of 𝑄 is a global-to-local copy operation, whereas application of 𝑄𝑇 involves summation (contraction) of local basis coefficients to their global counterparts. The FEM discretization of the Poisson problem, −∇2 𝑢 = 𝑓 , is based on the weak form. Let (𝑣, 𝑢) ∶= ∫Ω 𝑣 𝑢 𝑑𝑉 be the 𝐿2 inner-product on Ω, perhaps approximated by a suitable quadrature rule. Then the standard Galerkin projection approach leads to the FEM system matrix, 𝐴, having entries 𝑎𝑖𝑗 = (∇𝜙𝑖 , ∇𝜙𝑗 ), and right-hand side, 𝑏, with entries 𝑏𝑖 = (𝜙𝑖 , 𝑓 ). It is easy to show that 𝐴 = 𝑄𝑇 𝐴𝐿 𝑄, where 𝐴𝐿 =block-diag(𝐴𝑒 ) comprises local system matrices having entries 𝑎𝑒𝑖𝑗 =
∫Ω𝑒
∇𝑙𝑖𝑒 ⋅ ∇𝑙𝑗𝑒 𝑑𝑉
(2)
Application of 𝑄𝑇 and 𝑄 to 𝐴𝐿 is referred to as the matrix assembly process, which we will use later in the development of our reduced coarse-space operator. For applications on 𝑃 distributed-memory compute units (i.e., processes or MPI ranks), it is natural to cluster the elements into 𝑃 contiguous subdomains using, say, recursive spectral bisection [29] or other partitioning strategy that minimizes the number of shared dofs on the subdomain interfaces. We denote these subdomains as Ω𝑝 , 𝑝 = 1, … , 𝑃 , and map data and solution values associated with Ω𝑝 to process 𝑝. Nodal data on 𝜕Ω𝑝 , the boundary of Ω𝑝 , is redundantly represented on any process that shares those boundary nodes. With this element-based partition, nearestneighbor communication in our implementation typically arises when computing matrix-vector products of the form, T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 3 of 14
Fast Coarse Solvers
𝑤 = 𝐴𝑢. In local form, this product is expressed as 𝑟𝐿 = 𝐴𝐿 𝑢𝐿 , which is communication free because 𝐴𝐿 is blockdiagonal, followed by direct-stiffness summation, 𝑤𝐿 = 𝑄𝑄𝑇 𝑟𝐿 , in which shared nodal values are summed by 𝑄𝑇 and redistributed by 𝑄 [9]. For the Schwarz method, we extend each subdomain Ω𝑝 by one or more layers of elements from neighboring processes to form a set of overlapping subdomains, Ω𝑝 , with overlap 𝛿 = 𝑂(ℎ). For each domain, we define a Boolean restriction operator, 𝑅𝑝 , such that 𝑢𝑝 ∶= 𝑅𝑝 𝑢 returns the vector of nodal values that are interior to Ω𝑝 (i.e., excluding values on 𝜕Ω𝑝 ). Note that 𝑤𝑝 ∶= 𝑅𝑇𝑝 𝑢𝑝 extends (by zero) a vector of local subdomain values on Ω𝑝 to a global vector of length 𝑛. With this decomposition, the additive Schwarz preconditioner is 𝑧
=
𝑃 ∑
̃𝑇 𝐴−1 𝑅𝑝 𝑟 + 𝐽 𝐴−1 𝐽 𝑇 𝑟, 𝑅 𝑝 𝑝 𝑟
(3)
𝑝=1
which is the sum of local subdomain operators (subscript 𝑝) and a reduced coarse grid operator (subscript 𝑟), which we introduce later. The method is naturally parallel because the subdomain problems can be solved independently. The reduced coarse-space system, 𝐴𝑟 ∶= 𝐽 𝑇 𝐴𝐽 , with 𝐽 a coarse-to-fine interpolation matrix, is not trivially solved in parallel but it is of modest size compared to 𝐴 and can be solved using an 𝑋𝑋 𝑇 factorization, as we discuss below. The local systems are 𝐴𝑝 ∶= 𝑅𝑝 𝐴𝑅𝑇𝑝 , which constitute the principal submatrices of 𝐴 that correspond to the 𝑇
̃ , which is equal to 𝑅𝑇 for the classic interior nodes of Ω𝑝 . In (3), we also introduce the prolongation operator, 𝑅 𝑝 𝑝 ASM. Once the vectors are formally extended by zero, they can be added. Thus, summation of the local solutions, 𝑧̄ 𝑝 ∶= 𝑅𝑇𝑝 𝐴−1 𝑝 𝑅𝑝 𝑟, amounts to summing shared nodal values in regions of overlap. These local solution values must be communicated so that each neighboring process can complete the sum for all nodes in Ω𝑝 . Alternatively, one can ̃𝑇 such that each process ignores all nonlocal components from the local Schwarz solve. This approach yields define 𝑅 𝑝 the restricted additive Schwarz (RAS) method of Cai and Sarkis [6], which is in some cases more rapidly convergent ̃𝑇 is the same as 𝑅𝑇 , save that columns associated with Ω𝑝 ∖Ω𝑝 are null, so less than the traditional ASM. For RAS, 𝑅 𝑝 𝑝 communication is required. One still needs communication to enforce continuity for the redundantly-stored surface nodes, however. It is clear that the local solutions of the Schwarz method require a minimal amount of communication. Only surface values need to be communicated (particularly in the minimal-overlap case). On GPUs, one can solve (exactly or approximately) the local problem on the device. There is no need to communicate 𝑂(𝑛∕𝑃 ) data between the device and the host. Unlike multigrid, two-level Schwarz incurs nearest neighbor communication only on the fine level. The remaining communication is deferred to the coarse problem. As is well known, Schwarz methods must be augmented with a coarse space approximation to realize convergence rates that are independent of 𝑃 (i.e., that are scalable) [30]. The most common approaches are the additive variant (3) or a hybrid method that has a locally additive step, ASM1 :
𝑧loc
=
𝑃 ∑
̃𝑇 𝐴−1 𝑅𝑝 𝑟, 𝑅 𝑝 𝑝
(4)
𝑝=1
followed by a coarse-space correction, ASM2 :
=
𝑧
𝑇 𝑧loc + 𝐽 𝐴−1 𝑟 𝐽 (𝑟 − 𝐴𝑧loc ).
(5)
The multiplicative approach (4)–(5) serializes the local and coarse solve steps but is often more effective than (3) in reducing the overall number of iterations in the full system. In (3) and (5), 𝐽 is an interpolation matrix from the reduced coarse-space basis functions to the nodal points, 𝐱𝑗 . The main idea is that the interpolated range of 𝐽 should provide an efficient representation of the low-wave number components of the Poisson operator on Ω. For nested coarse spaces with exact local and coarse solves, the condition number for the preconditioned system is 𝜅(𝑀 −1𝐴)
=
𝑂(1 + 𝐻∕𝛿),
(6)
where 𝛿 is the amount of domain overlap and 𝐻 is the diameter of the subdomains, which is assumed to be the same as the diameter of the coarse-space elements. Cai [5] gives an alternative estimate for the nonnested case, T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 4 of 14
Fast Coarse Solvers
where 𝐻 a characteristic size of the coarse-space elements (i.e., boxes, in the notation of the next section), which is independent of the subdomain size. Under slightly different constraints on the coarse space, Chan, Smith, and Zou [7] give the following bound for additive Schwarz with nonnested spaces, ) ( 𝐻 2 . (7) 𝜅(𝑀 −1 𝐴) ≤ 𝐶 1 + 𝛿 The authors remark (Remark 5 [7]) that this bound can be improved to ( ) 𝐻max −1 𝜅(𝑀 𝐴) ≤ 𝐶 1 + , 𝛿
(8)
if the subdomains Ω𝑝 form a quasi-uniform triangulation of Ω and if 𝐻 ≤ 𝛽𝐻max for some fixed constant 𝛽 where 𝐻max = max diamΩ𝑝 .
4. Nonnested Coarse Spaces As noted earlier, a coarse-space is essential for efficient iterative solution of large-scale problems. The role of the reduced coarse space is to provide a mechanism for generating long wavelength approximations to the solution with relatively few degrees of freedom. Without such an approximation, many iterations may be required for lowwavenumber components of the solution to emerge in the iteration process, particularly if one is using restartedGMRES rather than a projection-based iteration. For unstructured mesh problems, it is not immediately obvious how to construct a coarse space that utilizes the underlying approximation space, 𝑋0ℎ . Aggregation (either smoothed or unsmoothed) and fine-coarse (FC) splittings are two of the most popular approaches in AMG. Nonnested coarse spaces have been discussed in both multigrid and Schwarz contexts [5, 4, 32]. Here, we follow the nonnested approach as it affords significant flexibility and the desired level of granularity for our targeted exascale platforms. Moreover, we can leverage the communication-minimal 𝑋𝑋 𝑇 coarse grid solver that is designed for problems at this scale [15, 31]. Without the restriction of nested spaces, one can develop a reduced coarse space 𝑋 𝑟 based on relatively simple interpolants. Let 𝑋̂ 𝑟 =span{Φ𝑗 (𝐱)}, 𝑗 = 1, … , 𝑛𝑟 , define a set reduced-space interpolation functions on Ω𝑟 ⊃ Ω and define the interpolation matrix 𝐽 as having entries 𝐽𝑖𝑗 ∶= Φ𝑗 (𝐱𝑖 ). Then the reduced approximation space is given by 𝑋 𝑟 =span{Ψ𝑖 (𝐱)} ⊂ 𝑋0ℎ , where Ψ𝑗 (𝐱)
=
𝑛 ∑
𝜙𝑖 (𝐱)𝐽𝑖𝑗 .
(9)
𝑖=1
Note that for any nodal point 𝐱𝑖 in the originating FEM mesh we have Ψ𝑗 (𝐱𝑖 ) = 𝐽𝑖𝑗 = Φ𝑗 (𝐱𝑖 ). However, the functions Ψ𝑗 are properly in 𝑋0ℎ —they satisfy the original continuity requriments and boundary conditions—whereas the Φ𝑗 s are not so constrained. We can choose the basis for 𝑋̂ 𝑟 based solely on efficiency and ease-of-implementation considerations. In fact, for advection-diffusion problems there are stability advantages to admitting discontinuous bases in 𝑋̂ 𝑟 , as we illustrate in the Appendix. Discontinuous bases for 𝑋̂ 𝑟 are also attractive for local mesh refinement in the reduced space as one requires no special treatment for hanging nodes. With the definition (9), the Galerkin statement for the Poisson problem in the reduced space is, Find 𝑢𝑟 (𝐱) ∈ 𝑋 𝑟 ⊂ ℎ 𝑋0 such that (∇𝑣, ∇𝑢𝑟 )
=
(∇𝑣, ∇𝑢) ∀ 𝑣 ∈ 𝑋 𝑟 .
(10)
In matrix form, the reduced-space basis coefficients 𝑢𝑟 = [𝑢𝑟,1 … 𝑢𝑟,𝑛𝑟 ]𝑇 are governed by the reduced system, 𝐴𝑟 𝑢𝑟 = 𝐽 𝑇 𝑏, with the SPD system matrix, 𝐴𝑟 ∶= 𝐽 𝑇 𝐴𝐽 . Here, we take 𝑋̂ 𝑟 to be the space of tensor products of piecewise linear interpolants that cover a rectangular box in 𝑑 lR . The 1D linear bases give rise to box-like elements in 2D and 3D, as illustrated in Fig. 1. With this basis, it is easy to construct a structured coarse-space to interpolate to the fine nodes. Information regarding topology or boundary conditions associated with Ω is not required. Aside from some solvability constraints (discussed below), the domain of 𝑋 𝑟 can be a simple superset of Ω, which implies that minimal geometric information is required to construct a viable set of interpolants, Φ𝑗 (𝐱). We simply find the limits of Ω in each of the principal directions and partition the space into T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 5 of 14
Fast Coarse Solvers
a set of 𝑛𝑏 = 𝑏𝑥 × 𝑏𝑦 (×𝑏𝑧 ) boxes that span these dimensions. We choose the number of boxes to be a small multiple of 𝑃 (e.g., 4–10) so that there are just a few dofs per process and to ensure that the box size, 𝐻, is greater than the characteristic finite element size, ℎ. The number of reduced-space unknowns is 𝑛𝑟 ≤ (𝑏𝑥 + 1)(𝑏𝑦 + 1)(𝑏𝑧 + 1), with the upper-bound corresponding to the total number of box vertices. Inequality holds for cases when there are empty boxes or, more precisely, when a larger value for 𝑛𝑟 would yield linearly dependent columns in 𝐽 . In the example of Fig. 1, the reduced-space problem has 16 unknowns, one for each vertex associated with the coarse bilinear interpolants, here denoted with a dual subscript, Φ𝑖𝑗 , 𝑖, 𝑗 ∈ {0, … , 3}2 . Figure 1(b) illustates the vertex-oriented interpretation of the basis functions, with the support for Φ22 shown in blue. The support for Ψ22 , which includes the blue and the red triangles, can be seen to extend beyond the squares that surround Φ22 . Figure 1(c) gives an “element” based interpretation of the support for the center-most box. Here, every triangular element marked in blue is influenced only by the four basis functions, Φ11 , Φ21 , Φ12 , and Φ22 . The red triangles indicate border elements whose support spans both the blue region and the adjacent boxes. A consequence of these overlapping elements is to increase the overall support of the individual basis functions, Ψ𝑖𝑗 (𝐱), which leads to reduced sparsity in 𝐴𝑟 . We address this issue in the next subsection. Note that there is no need to integrate or differentiate the interpolants Φ𝑖𝑗 . All of the physics is contained in 𝐴 and consequently embedded in 𝐴𝑟 ∶= 𝐽 𝑇 𝐴𝐽 . If 𝐽 is of full rank and 𝐴 is SPD, then 𝐴𝑟 will 𝑇 be SPD and the reduced coarse problem, 𝑧𝑟 = 𝐽 𝐴−1 𝑟 𝐽 𝑟, will be solvable. This system will provide the best 𝐴-norm −1 approximation to the solution 𝐴 𝑟 in the range of 𝐽 . By construction, 𝑧𝑟 will provide the desired long-wavelength approximation to the unknown error. Note that the holes in the domain of Fig. 1 pose no particular difficulty. If the boundaries of the holes are part of 𝜕Ω𝐷 the solution will vanish on those boundaries and the reduced approximation, 𝑧𝑟 , will respond accordingly. Similarly, if the boundary for the hole is part of 𝜕Ω𝑁 the interpolated solution will float according to the nontrivial values in 𝑧𝑟 associated with nodes 𝐱𝑖 ∈ 𝜕Ω𝑁 . A significant attraction of the nonnested spaces is that one can generate a simple coarse space that can be applied without detailed attention to geometry or boundary conditions. Situations where concern for geometry is warranted, however, include external flows such as flow over aircraft. In these cases, one might find that the far-field elements are larger than 𝐻, which can result in boxes that are interior to Ω but which contain no fine-scale vertices, potentially leading to null columns in 𝐽 since those bases would have no target interpolation points. Unlike vertices exterior to Ω, interior Φ𝑗 s cannot be trivially discarded because they are needed to formally provide function continuity across Ω. In such situations it is necessary to coarsen the reduced space and perhaps augment it with local mesh refinement in the near-field regions. The flexibility of the current approach provides a significant benefit in this case because one does not need continuity of Φ𝑗 at hanging nodes. In a high-level language like Matlab, implementation of the reduced space solver is straightforward. One begins with a sparse matrix 𝑛 × 𝑛𝑟 matrix 𝐽 that is void and then marches through each nodal point, 𝐱𝑘 , 𝑘 = 1, … , 𝑛, to identify its bounding box. Given the bounding box, one can identify each of the 2𝑑 box vertices and construct the linear interpolants that map the box vertex values to the node in question. Consider for example a 2D case with a 𝑏𝑥 × 𝑏𝑦 array of boxes. Let the columns of 𝐽 , which represent the discrete support of Φ𝑖𝑗 , be ordered lexicographically. Then 𝐽
=
[Φ00 (𝐱) Φ10 (𝐱)
⋯
Φ𝑏𝑥 𝑏𝑦 (𝐱)],
(11)
where 𝐱 = [𝐱1 𝐱2 ⋯ 𝐱𝑛 ]𝑇 is the vector of nodal points. Because the support of each Φ𝑖𝑗 is compact, 𝐽 will be sparse with at most 2𝑑 𝑛 nonzeros for a 𝑑-dimensional problem. (Each nodal point 𝐱𝑖 is influenced by the corner vertices of its bounding box, so each row of 𝐽 will have 2𝑑 nonzeros.) Construction of 𝐽 proceeds by finding and evaluating the four bounding interpolants for each 𝐱𝑖 . One can flag the vertices in each box that has a nontrivial entry. Any bounding vertex that is not flagged is compressed out of 𝐽 . In our production 𝑋𝑋 𝑇 solver that is used to solve the 𝐴𝑟 system, the setup phase does not require contiguous index sets, so columns that are not passed into the setup call are automatically suppressed. It is helpful to recognize that the action of 𝐽 is to copy the reduced vertex values to the 2𝑑 cells in the box domain that surrounds each vertex. That vertex value is interpolated onto each FEM node within a given cell. The action of 𝐽 𝑇 , therefore, is to multiply the FEM nodal values by the same weights and to then sum these onto the 2𝑑 indvidual vertex values. One can thus cast the formation of 𝐴𝑟 as a matrix assembly process. Consider the 2D case as an example and let 𝑆𝑖𝑗 denote the cell subdomain covered by the support of {Φ𝑖−1,𝑗−1 Φ𝑖,𝑗−1 Φ𝑖−1,𝑗 Φ𝑖,𝑗 }. Let 𝐴𝑖𝑗 denote the 4 × 4 matrix associated with cell 𝑆𝑖𝑗 that accounts for finite elements interior to 𝑆𝑖𝑗 . (This matrix is the cell equivalent of 𝐴𝑒 introduced in the FEM context.) For elements Ω𝑒 that are interior to 𝑆𝑖𝑗 , one can form a local 𝑛𝑣 × 4 interpolant, T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 6 of 14
Fast Coarse Solvers
𝐽𝑖𝑗𝑒 , that maps the four vertex values for cell 𝑆𝑖𝑗 to the 𝑛𝑣 nodes in Ω𝑒 . Let 𝐸𝑖𝑗 be the set of elements that are interior to 𝑆𝑖𝑗 and define ∑ 𝐴𝑖𝑗 = (𝐽𝑖𝑗𝑒 )𝑇 𝐴𝑒 𝐽𝑖𝑗𝑒 . (12) 𝑒∈𝐸𝑖𝑗
In an ideal situation, one could simply assemble each of these 4 × 4 matrices to form 𝐴𝑟 . Unfortunately, this is not quite the case because the procedure above applies only to elements interior to 𝑆𝑖𝑗 (i.e., the blue triangles in Fig. 1(b)). If we account for the red triangles then the simple assembly procedure breaks down because of the extended support of the Ψ𝑖𝑗 s. Possible mitigation options are to ignore the red triangles or to assign them decisvely to one cell or another (e.g., by the location of their centroid). These ideas are illustrated in Fig. 2. Ignoring the red triangles is roughly equivalent to using an incomplete quadrature rule that skips elements cut by cell boundaries. The problem with this approach is that one might have a chain of isolated triangles with vertices on 𝜕Ω𝑁 , each cut by a cell boundary. Representation of those vertex values would be lost by the discard rule and the reduced solution would not be propagated to those vertices. Assigning elements (as opposed to vertices) to cells, on the other hand has several advantages. First, all nodal values are represented. Second, any cell that has at least one element will be assigned at least 𝑛𝑣 rows in the corresponding columns of 𝐽 , thereby reducing the likelihood that columns of 𝐽 will be linearly dependent. One downside of the centroid assignment approach is that the output is no longer continuous. This situation can be remedied, however, by applying 𝑊 𝑄𝑐 𝑄𝑇𝑐 to 𝑧𝑟 = 𝐽 𝐴−1 𝑟 𝐽 𝑟, where 𝑄𝑐 is the coarse (𝑁 = 1) space assembly operator, 𝑊 = diag(𝑤−1 ) is a diagonal weight matrix comprising the inverse entries of the 𝑖 𝑇 𝑇 counting vector, 𝑤𝐿 = 𝑄𝑐 𝑄𝑐 𝑒𝐿 , and 𝑒𝐿 = [1 1 … 1] is the local unit vector of length 𝐸𝑛𝑣 . (This same weighting procedure is used to enforce continuity at the end of the Schwarz step, (3) or (4), and can be combined in the additive case.) By using either of these approaches, one can build 𝐴𝑟 through matrix assembly. The cell-based localization leads to a significant gain in sparsity, particularly in 3D. With the full support for Ψ, 𝐴𝑟 will have 5𝑑 nonzeros per row, as illustrated in Fig. 2(a) for 𝑑 = 1. With either of the compact schemes, the number of nonzeros per row is only 3𝑑 , which yields almost a 5× savings for 𝑑 = 3. Finally, switching from a width-5 stencil to a width-3 stencil reduces by a factor of two the size of the separator sets when 𝐴𝑟 is permuted by a nested-dissection ordering. For Cholesky factorization, this 5-to-3 switch reduces the number of nonzeros by a factor of four. For the 𝑋𝑋 𝑇 factorization discussed below, the number of nonzeros and the communication overhead, both of which scale with the aggregate size of the separators, are reduced by a factor of two.
5. Results In this section, we explore the performance of the new coarse grid solver in the context of the SEM-based incompressible flow solvers, Nek5000 and NekRS. A critical question in the current development was whether the outer, 𝑝MG-preconditioned, GMRES iteration counts for the pressure solves would change when the default AMG coarse solver was replaced by the proposed two-level Schwarz solve. If the reduction of the coarse solve time using the new coarse grid solver can’t compensate for the increased cost due to potentially large outer iterations, then the new approach will not be competitive. To answer this question, we compare AMG and the new solver as a coarse grid solver in four production-level Navier-Stokes simulations. The example cases used are illustrated in Fig. 3. The first is flow in a T-junction, which has a mesh comprising 𝐸 = 62176 spectral elements of order 𝑁 = 7. The reduced-space domain is based on a box of having 24 × 2 × 7 cells. Only a small fraction of these are used, however, as only those cells that contain spectral element centroids are activated. These cells are outlined in the figure. Note that the majority of the active reducedspace vertices are exterior to the domain in this case. The second test case is flow through a 146-pebble bed having 𝐸 = 62138 elements of order 𝑁 = 7, which is illustrated in the exploded view in Fig. 3. The fluid flows around the spherical pebbles and the sphere interiors are external to Ω. Here, we use a 7 × 7 × 12 box mesh for the reduced space, which nominally has 832 DOFs save that corner cells in the 𝑥-𝑦 plane fall outside the disk, so those corner DOFs are not included. This example shows that complex domains with void inclusions present no particular difficulty for the new coarse-space approach. The third test is also a pebble bed but with ≈ 45000 pebbles with 𝐸 = 13032440 elements of order 𝑁 = 7. Reduced space is based on a box mesh of 46 × 46 × 23 for Nek5000 and 28 × 28 × 14 for NekRS. The forth test is vertical flow through an annular bed of 352625 pebbles with 𝐸 = 98782067 elements of order 𝑁 = 8 (𝑛 = 51B). The reduced space is based on a box mesh of 47 × 47 × 35 for Nek5000 and 32 × 32 × 25 for NekRS. Again, any box that is fully exterior to Ω is automatically dropped by the implementation, which only adds T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 7 of 14
Fast Coarse Solvers
active vertices to the set of reduced-space unknowns. For each test case, the pressure-Poisson boundary conditions are Neumann everywhere except at the domain outflow (the end of the long branch for the T-junction case and the upper boundary for the pebble bed cases), where the solution has homogeneous Dirichlet conditions. Although these boundary conditions are not applied to the reduced-space basis functions, Φ𝑗 , they are automatically imposed on Ψ𝑗 . The original coarse-space dimension is 𝑛𝑐 ≈ 𝐸 for all the four cases since it is the space associated with the last level of 𝑝MG (𝑁 = 1) which consists of element vertices. We run all the tests in the context of full Navier-Stokes simulations, which require the solution of a pressure Poisson problem at each time step, 𝑡𝑚 . For iterative solvers, it pays to solve only for the change in the solution, 𝛿𝑝 ∶= 𝑝𝑚 − 𝑝̄ , ∑ where 𝑝̄ = 𝑝𝑚−1 or, better, 𝑝̄ = 𝑙𝑗=1 𝛽𝑗 𝑝𝑚−𝑗 is the 𝐴-orthognal (best fit) projection onto 𝑙 prior solutions [16]. This projection process tends to remove low-wavenumber content from the initial residual and typically yields a two- to three-fold savings in Navier-Stokes run time [16, 25]. All Nek5000 runs were performed on the Summit supercomputer at OLCF, which has 4608 nodes, each with two IBM POWER9 CPUs and six NVIDIA Volta V100 GPUs. Nek5000 only utilized the CPU cores since it does not support running on the GPUs. NekRS runs were performed on the Frontier supercomputer at the same facility which has ≈ 9400 nodes, each with a single AMD EPYC CPU and four MI250X GPUs. NekRS was run on the GPUs using HIP backend of OCCA [23] portability library except for the coarse solver which was run on the CPUs in case of two level Schwarz solver and the default AMG solver using BoomerAMG. The first set of results, shown in Fig. 4, are for Nek5000, which runs on CPUs and uses multilevel 𝑝MG with Schwarz-based smoothing for the finest two levels, 𝑁 = 7 and 𝑁 = 3. The default for the coarse 𝑁 = 1 problem is to solve it directly with the 𝑋𝑋 𝑇 factorization or, for problems with 𝐸 > 350000, to solve it approximately with a single sweep of an AMG 𝑉 -cycle using Hypre [20].The Nek5000 preconditioner, described in [18, 22], is additive between levels and takes the form ( ) ) ( −1 + 𝐽3 𝑀3−1 + 𝐽1 𝑀1−1 𝐽1𝑇 𝐽3𝑇 𝑟, (13) 𝑧 = 𝑀𝑁 where 𝑀1−1 = 𝐴−1 𝑐 represents the coarse solve and −1 𝑀𝑁
=
𝑊𝑁
𝐸 ∑
𝑅𝑇𝑁,𝑒 𝐴̃ −1 𝑁,𝑒 𝑅𝑁,𝑒
(14)
𝑒=1
represents the element-local additive Schwarz smoother for polynomial order 𝑁 > 1. The local Schwarz problems are solved (approximately) using the tensor-product based fast diagonalization methods described in [19, 22]. The diagonal weight matrix, 𝑊𝑁 , is used to average the solution in the element overlap regions, which was found to significantly improve the smoothing properties of the ASM [22]. A central idea in the original Nek5000 ASM implementation was to ensure that the result of each matvec was projected. Consequently, only a single matvec is used per GMRES iteration, but the overall iteration count is higher than that of NekRS, which uses multiple smoothing sweeps at each 𝑝MG level. Figure 4 shows the pressure iteration counts (or number of coarse grid solves) per time step for each of the examples when using several different coarse solvers for the multilevel 𝑝MG preconditioner in Nek5000. The top row shows the results as number of iterations per step, whereas the bottom shows the cumulative iteration counts, which are easier to read given the step-to-step variability in iteration counts. (This variability is largely due to the pressure projection, which extracts all temporal regularity out of the solution such that the initial residual is highly variable [16] .) Shown 𝑇 are five different coarse solvers: using just the global reduced-space (box) solver, 𝑀𝑟−1 ∶= 𝐽 𝐴−1 𝑟 𝐽 ; using just the local Schwarz solve, ASM1 (4); combining these into a two-level Schwarz solve, ASM2 (5); using a single 𝑉 -cycle of Hypre; and using the direct 𝑋𝑋 𝑇 -based solve. First off, we note that both parts of the ASM2 solver are crucial for the success of the method. Neither the local solve nor the reduced-space box-solve are sufficient to yield a fast 𝑝MG solver. Moreover, we see that ASM2 , with the sparse box-solver for the coarse grid problem, is almost as effective as 𝑋𝑋 𝑇 and slightly more effective than a single AMG 𝑉 -cycle. It is remarkable just how effective the low-communication coarse space is for these applications.
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 8 of 14
Fast Coarse Solvers Table 1 Comparison of solver times and pressure iterations per timestep with NekRS when using AMG and two level Schwarz solver as a coarse grid solve. Case T-Junction T-Junction 146 Pebbles 146 Pebbles 45K Pebbles 45K Pebbles 350K Pebbles 350K Pebbles
Solver Schwarz AMG Schwarz AMG Schwarz AMG Schwarz AMG
P 8 8 8 8 1680 1680 27648 27648
E/P 7772 7772 7767 7767 7758 7758 3573 3573
n/P 2.66e+06 2.66e+06 2.66e+06 2.66e+06 2.66e+06 2.66e+06 1.83e+06 1.83e+06
NS Time(s) 1.37e-01 1.21e-01 1.14e-01 1.14e-01 2.60e-01 2.69e-01 2.17e-01 2.18e-01
Pres. Time(s) 8.59e-02 6.97e-02 6.72e-02 6.63e-02 1.95e-01 2.05e-01 1.59e-01 1.61e-01
Coarse Time(s) 2.57e-02 1.78e-02 1.53e-02 1.68e-02 4.66e-02 7.35e-02 1.70e-02 4.56e-02
Pres. Iter. 4.32 3.84 3.45 3.17 8.35 7.31 6.18 4.80
Second set of results, shown in Figure 5 are for the same set of cases simulated using NekRS running on GPUs except for the coarse solver which is running on the CPUs. We only plot the cumulative iterations counts for the cases with NekRS in Figure 5. NekRS uses a multilevel 𝑝MG with a Schwarz based smoother and a schedule 𝑁 = 7 → 3 → 1 similar to Nek5000 settings used in experiments in Figure 4. Default coarse solver in NekRS is a single sweep of an AMG 𝑉 -cycle on the original system 𝐴𝑐 (𝑁 = 1) using Hypre [20] in single precision on the CPU. Similar to what we saw with Nek5000, both parts of ASM2 are required for the success of the new coarse solver with NekRS as well. One main difference compared to Nek5000 experiments is that a single 𝑉 -cycle with Hypre is more effective in terms of the cumulative pressure iteration counts than the Schwarz based coarse solver for the NekRS experiments. Our primary interest is the time per timestep at large scale simulations rather than the cumulative iteration counts and the whole point of the design of the new Schwarz based coarse solver was to have a smaller cumulative coarse grid solve time compared to AMG albeit a slightly larger cumulative number of pressure iterations. Table 1 compare the solve times and pressure iterations per timestep when using AMG and two level Schwarz method as the coarse grid solver in NekRS for four cases including three cases used in experiments in Figure 4-5 and the annular bed of 352625 pebbles. The latter case uses a schedule 𝑁 = 8 → 6 → 4 → 1 for the multilevel 𝑝MG preconditioner. The results show that the two level Schwarz method is competitive (and does slightly better) at higher process counts where the communication cost of AMG becomes significant due to multiple levels in the 𝑉 -cycle. This can be seen for the two larger meshes run with 𝑃 = 1680 and 𝑃 = 27648 in Table 1 where the time spent in coarse grid solve per timestep is about 1.6 and 2.7 times faster than default AMG solver respectively. These speedups eventually translate to an overall speedup of full Navier-Stokes solve time per timestep. While, the number of pressure iterations per timestep is lower with AMG for all the cases, reduced coarse times eventually helped the two level method catch up with AMG at higher process counts for the two large cases with high number of MPI processes. Figure 6, which shows timing data for a strong scaling study of the 45000 pebbles mesh with NekRS further demonstrate this point. The plots were generated by collecting timing data by increasing the number of processes, 𝑃 , from 840 to 10080 (or decreasing 𝑛∕𝑃 ). We can see that Navier-Stokes solve time, pressure solve time and coarse grid solve time per timestep goes down faster with the two level Schwarz solver compared to AMG as the number of processes increase in the leftmost plot in Figure 6. Crossover point where the two level Schwarz solver becomes faster than the default AMG approach is 𝑛∕𝑃 ≈ 3.5 × 106 . Middle plot of Figure 6 shows the parallel efficiency for the two solvers as 𝑛∕𝑃 decreases. Reference point for the parallel efficiency is the NekRS run with Hypre AMG solver with 𝑃 = 840 which is assumed to have 100% parallel efficiency. It is evident from the middle plot that the new two level Schwarz solver has better parallel efficiency than AMG for 𝑛∕𝑃 values which falls in our original design range for 𝑛∕𝑃 . Right most plot in Figure 6 shows the breakdown of the coarse grid solve time for the two level Schwarz solver. Cost of the operations which don’t involve communication goes down noticeably as 𝑛∕𝑃 decreases. These include the local Schwarz solve (𝐴−1 𝑝 ), right hand side update (𝑟 − 𝐴𝑧𝑙𝑜𝑐 ) required for doing reduced solve in a multiplicative manner and interpolation (𝐽 ) from the original coarse space 𝐴𝑐 to the reduced structued space 𝐴𝑟 . Note that in the case of 𝐽 and 𝐽 𝑇 , it is the structured nature of 𝐴𝑟 which allows us to compute the interpolation and prolongation in parallel with no communication. The fact that the local Schwarz solve is the most expensive operation in the two level Schwarz for most of the range presents an opportunity for further optimizations which are discussed breifly in Section 6.3. Cost T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 9 of 14
Fast Coarse Solvers Table 2 Navier-Stokes solve times per timestep with different Chebyshev orders for 350K pebbles mesh with NekRS using 𝑃 = 27648 processes with ≈ 2𝑀 grid points per process when using AMG and two level Schwarz solver as coarse grid solvers. Solver AMG AMG AMG Schwarz Schwarz Schwarz
Cheb. Order 3 2 1 3 2 1
NS Time(s) 2.18e-01 2.60e-01 3.65e-01 2.17e-01 2.48e-01 3.31e-01
Description Latency of the Network Short-long message demarcation Arithmetic time (inverse FLOPS)
Pres.Time(s) 1.61e-01 2.02e-01 3.04e-01 1.59e-01 1.89e-01 2.75e-01 Symbol 𝛼 𝑚2 = 𝛼∕𝛽 𝑡𝑎
Coarse Time(s) 4.56e-02 7.02e-02 1.38e-01 1.70e-02 2.77e-02 5.68e-02
Pres. Iter. 4.80 7.28 14.50 6.18 9.78 20.90
Measured or Estimated value 4e-6 ≈ 5000 10−9 s
Table 3 System parameters for Frontier supercomputer at OLCF
of solving the reduced system 𝐴−1 𝑟 stays more or less the same as 𝑛∕𝑃 decreases. This is not suprising since the size of this system is not affected by 𝑛∕𝑃 but fixed for the entire study. Althoug a slight decrease in the cost of solving 𝐴−1 𝑟 is observed as 𝑛∕𝑃 decreases which hints at the fact that it is the local part of 𝑋𝑋 𝑇 which is dominating the cost of solving 𝐴−1 𝑟 . Table 2 illustrate the importance of visiting the coarse grid system as few times as possible in order to reduce the overall Navier-Stokes solve time. This is the strategy adapted in [25] in order to reduce the coarse grid overhead indirectly by doing more Chebyshev smoothing steps in order to reduce the number of pressure solves. Increasing number of Chebyshev smoothing steps from 1 to 3 decreases the number of pressure iterations by more than a factor of three.
6. Discussion In this section, first we are going to do a complexity analysis of both AMG and two level Schwarz method to see if we are on the right ballpark in terms of the timing data collected and to estimate the costs for future exascale runs. Then we are going to analyse the affect of the new two level Schwarz coarse solver on strong scaling of Nek5000/RS. Finally, we are going to discuss the future work that we are planning to do in order to improve the performance of the solver.
6.1. Complexity Analysis Here we develop time complexity estimates for parallel evaluation of the two-level Schwarz and AMG in the context of a coarse grid solver for pMG preconditioner. We are considering a problem size consitent with the strong scaling limit for GPUs, which is about 𝑛∕𝑃 ≈ 2.66 × 106 or 𝐸∕𝑃 ≈ 8000 with 𝑁 = 7. For Schwarz, we assume that the local and reduced coarse-space solves are evaluated serially in a multiplicative way, so that one might consider using a hybrid Schwarz approach in which the additive local solves are followed by a coarse grid correction. We define the required system parameters, along with estimated/measured values for Frontier supercomputer at OLCF in Table 3. Note that both AMG and Schwarz are run on the CPU in the context of the coarse solver, so we only need the parameters for the CPU and the network. For ASM1 , we have pre- and post-solve (𝑅𝑝 and 𝑅̃ 𝑇 ) near-neighbor communication requiring ≈ 50 messages of size ≤ (𝐸∕𝑃 )2∕3 . For GPUs at the strong-scale limit, we anticipate 𝐸∕𝑃 ≈ 8000 and therefore message sizes 𝑚 ≤ (𝐸∕𝑃 )2∕3 = 400 < 𝑚2 , which implies that the messages are latency bound, each with a cost ≤ 2𝛼. Thus, anticipated communication cost for ASM1 is 𝑇𝑐,𝐴𝑆𝑀1
=
50 ⋅ (2𝛼) = 100 ⋅ 𝛼 ≈ .0002s
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
(15) Page 10 of 14
Fast Coarse Solvers
Above communication time is the sum of the time for 𝑅𝑝 and 𝑅̃ 𝑇 shown in the rightmost plot of Figure 6. Cost of 𝑅𝑝 and 𝑅̃ 𝑇 are more or less the same, and thus should be half of the above estimate. We can see that the timing data for 𝑅𝑝 and 𝑅̃ 𝑇 for 𝐸∕𝑃 ≈ 8000 in rightmost Figure 6 is in the same order of magnitude as the estimate above. Also, there is a local computation to be done for 𝑅𝑝 and 𝑅̃ 𝑇𝑝 in addition to the communication which explains why the measured time goes down slightly as local problem size 𝑛∕𝑃 goes down. Local solve for ASM1 is implemented as a Cholesky solve on the CPU, so we can estimate the time using the CPU arithmetic time 𝑡𝑎 and the number of operations performed during the for the local solve. Assuming a Cholesky factorization of 𝐴𝑝 = 𝐿𝐿𝑇 , the number of operations in forward and backward solve is roughly 2𝑛𝑧 where 𝑛𝑧 is the number of nonzeros in 𝐿 (we have to access each non-zero value in 𝐿 or 𝐿𝑇 once during the solve and update the right hand side by substracting its contribution after the multiplication by relevant solution component). We have noticed that 𝑛𝑧 could be as high as 5 × 106 for the local problems when 𝐸∕𝑃 ≈ 8000 (or 𝑛∕𝑃 ≈ 2.66 × 106 for 𝑁 = 7). =
𝑇𝑎,𝐴𝑆𝑀1
2𝑛𝑧 𝑡𝑎 ≈ 0.01s
(16)
With both communication in Equation 15 and arithmetic time in Equation 16, we can estimate the total time for ASM1 as: 𝑇𝐴𝑆𝑀1
=
𝑇𝑐,𝐴𝑆𝑀1 + 𝑇𝑎,𝐴𝑆𝑀1 ≈ 0.01s
(17)
For the reduced-space solve 𝐴𝑟 𝑢𝑟 = 𝑏𝑟 , we follow the complexity estimate for the 𝑋𝑋 𝑇 solver developed in [31]. For a regular 3D mesh (which provides a conservative bound), the estimated communication and arithmetic times for 𝑋𝑋 𝑇 are: ) ] [( 2.5 2∕3 log2 𝑃 𝛼 (18) 𝑛𝑟 𝑇𝑐,𝐴−1 = 2 1 + 𝑟 𝑚2 ( ) 5∕3 𝑇𝑎,𝐴−1 = 2 2.5 𝑛𝑟 ∕𝑃 𝑡𝑎 (19) 𝑟
For the 45000 Pebbles problem, reduced system size 𝑛𝑟 is 28 × 28 × 14 ≈ 11000 for NekRS. If we consider 𝑃 = 1680 (for which 𝐸∕𝑃 ≈ 8000 and 𝑛∕𝑃 = 2.66 × 106 ), we can calculate the arithmetic, communication and total time for the reduced-space solve as follows: [( ) ] 2.5 𝑇𝑐,𝐴−1 = 2 1 + × 110002∕3 × log2 1680 × 4 × 10−6 (20) 𝑟 5000 [( ) ] 2.5 = 2 1+ × 500 × 11 × 4 × 10−6 = 2 × 13.5 × 4 × 10−6 = 0.0001s (21) 5000 ( ) 𝑇𝑎,𝐴−1 = 2 2.5 × 110005∕3 ∕1680 × 1 × 10−9 = 2 × (2.5 × 5.4𝑒6∕1680) × 1 × 10−9 (22) 𝑟
𝑇𝐴−1 𝑟
=
2 × (2.5 × 3200) × 1 × 10−9 = 0.00001s
(23)
=
𝑇𝑐,𝐴−1 + 𝑇𝑎,𝐴−1 = 0.0001s
(24)
𝑟
𝑟
From leftmost plot in Figure 6, we can see that the measured time for the reduced-space is about 0.0003s which is the same order of magnitude as the complexity estimate. Also, we expect the actual run time to be higher than the estimate since we don’t have a nested dissection ordering for the reduced system and the Equation 19 is derived assuming a nested dissection ordering. We contrast these times with estimated and measured times for AMG applied to the full system 𝐴𝑢 = 𝑏 (i.e., 𝐴𝑐 𝑢𝑐 = 𝑏𝑐 in the pMG solver). BoomerAMG based coarse solver used in the experiments performed a single 𝑉 -cycle with a single smoothing step in each level using a Chebyshev smoother with degree 2. Thus, each smoothing step consists of two nearest neighbor communications with about 26 messages per communication. These messages are also latency bound similar to 𝑅𝑝 and 𝑅̃ 𝑇𝑝 in the ASM1 solver with a cost of 2𝛼. If there are 𝐿 levels in the 𝑉 −cycle, the communication cost of the full 𝑉 −cycle is given by: 𝑇𝑐,𝐴𝑀𝐺
=
2 ⋅ 𝐿 ⋅ (2 ⋅ 2𝛼 × 26) ≈ 200 ⋅ 𝐿 ⋅ 𝛼
(25)
The total arithmetic time is estimated as twice the cost of the smoothing step on the fine grid which is a sparse matvec with 𝐸∕𝑃 rows and ≈ 27 non-zeros per row. 𝑇𝑎,𝐴𝑀𝐺
=
2 ⋅ (27𝐸∕𝑃 ) ⋅ 𝑡𝑎 ≈ 50 ⋅ 𝐸∕𝑃 ⋅ 𝑡𝑎
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
(26) Page 11 of 14
Fast Coarse Solvers
Thus, the total time for AMG solver with a sing 𝑉 −cycle is given by: 𝑇𝐴𝑀𝐺
=
𝑇𝑐,𝐴𝑀𝐺 + 𝑇𝑎,𝐴𝑀𝐺 = 200 ⋅ 𝐿 ⋅ 𝛼 + 50 ⋅ 𝐸∕𝑃 ⋅ 𝑡𝑎
(27)
For the 45000 Pebbles problem, at 𝑃 = 1680, BoomerAMG had 𝐿 = 10 levels in the 𝑉 −cycle and the local problem 𝐸∕𝑃 ≈ 8000 similar to the Schwarz based approach. Communication and arithmetic times for the 𝑉 −cycle can then be calculated as: 𝑇𝑐,𝐴𝑀𝐺
=
200 ⋅ 10 ⋅ 𝛼 = 0.008s
𝑇𝑎,𝐴𝑀𝐺
=
50 ⋅ 8000 ⋅ 1 × 10−9 = 0.0004s
(29)
𝑇𝐴𝑀𝐺
=
𝑇𝑐,𝐴𝑀𝐺 + 𝑇𝑎,𝐴𝑀𝐺 = 0.0084s
(30)
(28)
Using the data in Table 1 for 45K Pebbles problem with AMG, we can calculate the time for a single coarse solve by dividing coarse time by the number of pressure iterations which comes out to 7.35𝑒 − 2∕7.31 ≈ 0.01s. We can see that the measured time for AMG matches pretty well with the estimated time.
6.2. Strong Scaling In this section, we investigate strong scalability of Schwarz and AMG based coarse-grid solvers for pMG both experimentally and using the timing estimates developed in Section 5 and show that the Schwarz has better strong scaling properties than AMG. As shown in Section 5, strong scaling experiments done on 45000 pebbles case shown in Figure 6 show that, as we increae the number of processes 𝑃 (thus decreasing 𝑛∕𝑃 ), time spent in both Schwarz and AMG (leftmost plot in Figure 6) decreases with Schwarz time decreasing faster than AMG. Middle plot in Figure 6 shows the parallel efficiency for the same experiment assuming BoomerAMG run with 𝑃 = 840 has 100% of effiiciency. An HPC user running production level simulations has to decide on the trade-off between parallel efficiency and the time to solution. It is clear from the plot that as we keep adding processors, time to solution keep decreasing till the parallel efficiency reaches 30-40%. But a low level of parallel efficiency may not be acceptable due to underutilization of computing resources. Let’s assume that users are willing to accept a parallel efficiency of 80%. Then, from the middle plot in Figure 6, we can see that we can use larger number of processors with Schwarz to achieve a smaller time to solution than AMG. We can see from the middle plot in Figure 6 that AMG needs about 3 × 106 gridpoints for 80% parallel efficiency whereas Schwarz needs ≈ 2 × 106 . We can verify the above observation by looking at the timing estimates developed in Section 5. Note that we are operating in the regime where local problem size, 𝐸∕𝑃 ≈ 8000 and 𝑃 ≈ 105 . For estimating time with AMG, we can use Equation 27 with 𝐿 ≈ 10 since we expect the number of levels in AMG 𝑉 -cycle to be around 10 for problems of size 𝐸 = 109 . For the Schwarz solver, our local solve times are going to stay constant since both local problem size and number of neighbors more or less stays constant as the number of processors grow. Only the reduced solve time will increase as we increase the number of processors. We can calculate the time spent in the reduced solve as the sum of Equation 18 and Equation 19. The reduced solve time is dependent on the number of reduced-space points 𝑛𝑟 . By construction, 𝑛𝑟 = 𝛾𝑃 , where 𝛾 = 3–10 is the target number of reduced-space points per rank. We assume 𝛼, and 𝑡𝑎 values from Table 3 for Frontier at OLCF will also hold true for future exascale systems. Figure 7 show the estimated solve times for Schwarz and AMG for different 𝛾 and 𝐿 values.
6.3. Conclusion We presented a two level Schwarz method as a coarse solver for the 𝑝-multigrid preconditioner in SEM/FEM simulations with a novel nonnested coarse space for the global coarse grid of the Schwarz method. Structured nature of this coarse space enables communication-free interpolation between the local problem and global coarse grid. We used CHOLMOD [8], a sparse Cholesky solver for solving the local problem exactly. CHOLMOD uses Approximate Minimum Degree (AMD) [2] ordering to reduce fill-in in the Cholesky factor. Global coarse problem is solved exactly using 𝑋𝑋 𝑇 solver. Although, these two solves can be performed either in an adiditive or a multiplicative manner, we found that the reduction of pressure iterations resulting from the latter is essential for the new solver to be competitive with state of the art AMG solvers. There are multiple aspects in the current solver that can be improved in future research work to make it more robust and faster. Currently, we are solving both local and global coarse problems of the Schwarz method in double precision. Switching to single precision will cut the communication cost by a factor of two thus could result in a considerable T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 12 of 14
Fast Coarse Solvers
speed up in global coarse solve time. Similarly, local Cholesky solve time can see a speedup too due to the low memory footprint of the Cholesky factor. Another effective way to reduce memory footprint and communication cost is to reduce nonzero entries in both the local Cholesky factor as well as the 𝑋𝑋 𝑇 factorization with a reordering of the respective sparse systems. Furthermore, reduction of the nonzeros in the factors reduce the number of floating point operations that has to be performed during a solve. Another avenue worth exploring is the use of inexact solvers for both local Schwarz problem as well as the global coarse grid problem. This could be a winning strategy if the increase of pressure iterations and subsequent increase of the time spent in finer levels of 𝑝-multigrid can be compensated by the reduction in time spent in the coarse solve. For the local problem, incomplete Cholesky on the CPU or an iterative method like CG with a suitable preconditioner on GPU can be used as an inexact solver.
Acknowledgments This material is based upon work supported by the U.S. Department of Energy, Office of Science, under contract DE-AC02-06CH11357 and by the Exascale Computing Project (17-SC-20-SC). The research used resources at the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC05-00OR22725.
References [1] Alvarado, F.L., Pothen, A., Schreiber, R., 1992. Highly parallel sparse triangular solution. Pennsylvania State University, Department of Computer Science. [2] Amestoy, P.R., Davis, T.A., Duff, I.S., 2004. Algorithm 837: Amd, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software (TOMS) 30, 381–388. [3] Anderson, R., Andrej, J., Barker, A., Bramwell, J., Camier, J.S., Cerveny, J., Dobrev, V., Dudouit, Y., Fisher, A., Kolev, T., et al., 2021. Mfem: A modular finite element methods library. Computers & Mathematics with Applications 81, 42–74. [4] Bramble, J.H., Pasciak, J.E., Xu, J., 1991. The analysis of multigrid algorithms with nonnested spaces or noninherited quadratic forms. Math. of Comp. 56, 1–34. [5] Cai, X.C., 1995. The use of pointwise interpolation in domain decomposition methods with nonnested meshes. SIAM J. Sci. Comput. 16, 250–256. [6] Cai, X.C., Sarkis, M., 1999. A restricted additive Schwarz preconditioner for general sparse linear systems. SIAM J. Sci. Comput 21, 792–797. [7] Chan, T.F., Smith, B.F., Zou, J., 1996. Overlapping schwarz methods on unstructured meshes using non-matching coarse grids. Num. Math. 73, 149–167. [8] Chen, Y., Davis, T.A., Hager, W.W., Rajamanickam, S., 2008. Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate. ACM Transactions on Mathematical Software (TOMS) 35, 1–14. [9] Deville, M., Fischer, P., Mund, E., 2002. High-order methods for incompressible fluid flow. Cambridge University Press, Cambridge. [10] Dryja, M., Widlund, O., 1987. An Additive Variant of the Schwarz Alternating Method for the Case of Many Subregions. Technical Report TR 339. Courant Inst., NYU. Dept. Comp. Sci. [11] Dryja, M., Widlund, O., 1989. An additive Schwarz algorithm for two- and three-dimensional finite element elliptic problems, in: Chan, T., Glowinski, R., Périaux, J., Widlund, O. (Eds.), Domain Decomposition Methods, SIAM. [12] Dryja, M., Widlund, O., 1994. Domain decomposition algorithms with small overlap. SIAM J. Sci. Comput. 15, 604–620. [13] Farhat, C., Chen, P., 1994. Tailoring domain decomposition methods for efficient parallel coarse grid solution and for systems with many right hand sides. Contemporary Mathematics 180, 401–406. [14] Fischer, P., 1993. Projection techniques for iterative solution of 𝐴𝑥 = 𝑏 with successive right-hand sides. Technical Report 93–90. ICASE. Hampton,Va. [15] Fischer, P., 1996. Parallel multi-level solvers for spectral element methods, in: Ilin, A., Scott, L. (Eds.), Third Int. Conference on Spectral and High Order Methods, Houston J. of Mathematics. pp. 595–604. [16] Fischer, P., 1998. Projection techniques for iterative solution of 𝐴𝑥 = 𝑏 with successive right-hand sides. Comput. Methods Appl. Mech. Engrg. 163, 193–204. [17] Fischer, P., Heisey, K., Min, M., 2015. Scaling limits for PDE-based simulation (invited), in: 22nd AIAA Computational Fluid Dynamics Conference, AIAA Aviation, AIAA 2015-3049. [18] Fischer, P., Lottes, J., 2004. Hybrid Schwarz-multigrid methods for the spectral element method: Extensions to Navier-Stokes, in: Kornhuber, R., Hoppe, R., Périaux, J., Pironneau, O., Widlund, O., Xu, J. (Eds.), Domain Decomposition Methods in Science and Engineering Series, Springer, Berlin. [19] Fischer, P., Miller, N., Tufo, H., 2000. An overlapping Schwarz method for spectral element simulation of three-dimensional incompressible flows, in: Bjørstad, P., Luskin, M. (Eds.), Parallel Solution of Partial Differential Equations, Springer, Berlin. pp. 158–180. [20] Henson, V., Yang, U., 2002. BoomerAMG: a parallel algebraic multigrid solver and preconditioner. Applied Numerical Mathematics 41, 155–177. [21] Hiptmair, R., 2002. Finite elements in computational electromagnetism. Acta Numerica 11, 237–339. [22] Lottes, J.W., Fischer, P.F., 2005. Hybrid multigrid/Schwarz algorithms for the spectral element method. J. Sci. Comput. 24, 45–78. [23] Medina, D.S., St-Cyr, A., Warburton, T., 2014. OCCA: A unified approach to multi-threading languages. arXiv preprint arXiv:1403.0968 .
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 13 of 14
Fast Coarse Solvers [24] Merrill, B., Peet, Y., Fischer, P., Lottes, J., 2016. A spectrally accurate method for overlapping grid solution of incompressible Navier-Stokes equations. J. Comput. Phys. 307, 60–93. [25] Min, M., Lan, Y., Fischer, P., Merzari, E., Kerkemeier, S., Phillips, M., Rathnayake, T., Novak, A., Gaston, D., Chalmers, N., Warburton, T., 2022. Optimization of full-core reactor simulations on Summit, in: Proc. Conf. on Supercomput., IEEE. [26] Mittal, K., Dutta, S., Fischer, P., 2020. Multirate time-stepping for the incompressible Navier-Stokes equations in overlapping grids. J. Comput. Phys. 437, 110335. [27] Patera, A., 1984. A spectral element method for fluid dynamics : laminar flow in a channel expansion. J. Comput. Phys. 54, 468–488. [28] Phillips, M., Kerkemeier, S., Fischer, P., 2022. Tuning spectral element preconditioners for parallel scalability on GPUs, in: Proc. of the 2022 SIAM Conf. on Par. Proc. for Sci. Comp., SIAM. pp. 37–48. [29] Pothen, A., Simon, H., Liou, K., 1990. Partitioning sparse matrices with eigenvectors of graphs. SIAM J. Matrix Anal. Appl. 11, 430–452. [30] Smith, B., Bjørstad, P., Gropp, W., 1996. Domain Decomposition: Parallel Multilevel Methods for Elliptic PDEs. Cambridge University Press, Cambridge. [31] Tufo, H., Fischer, P., 2001. Fast parallel direct solvers for coarse grid problems. J. Parallel Distrib. Comput. 61, 151–177. [32] Xu, J., 1996. The auxiliary space method and optimal multigrid preconditioning techniques for unstructured grids. Computing 56, 215–235.
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 14 of 14
Fast Coarse Solvers
(a)
Φ03
Φ13
Φ23
Φ33
Φ32
Φ02
Φ12
Φ22
Φ21
Φ31
Φ01
Φ11
Φ20
Φ30
Φ00
Φ10
Φ03
Φ13
Φ23
Φ33
Φ02
Φ12
Φ22
Φ01
Φ11
Φ00
Φ10
(b)
Φ03
Φ13
Φ23
Φ33
Φ32
Φ02
Φ12
Φ22
Φ32
Φ21
Φ31
Φ01
Φ11
Φ21
Φ31
Φ20
Φ30
Φ00
Φ10
Φ20
Φ30
(c)
Figure 1: Two-level Schwarz illustration. (a) Parallel partition (𝑃 = 2) of the 𝑁 = 1 mesh: each rank solves a local Poisson problem on their respective shaded region (green or blue), including a one- or two-element overlap extension, shown in gray. (b) Vertex-based support of the coarse grid interpolants: shown in blue is the support for Φ22 . The red triangles indicate border elements, for which support of Ψ22 would interact with Ψ0𝑘 and Ψ𝑘0 , 𝑘 = 1 and 2, which would in general give rise to a 5 × 5 stencil for these lexicographically-ordered coarse-space functions. (c) Element-based illustration of the coarse decomposition: if the (red) overlapping elements are eliminated then the coarse stencil will only be 3 × 3 and the coarse-space operator, 𝐴𝑟 , can be formed using standard FEM assembly techniques.
Φ𝑗−1
Φ𝑗+1
Φ
Ψ
Ψorig
Ψ𝑗−1
Ψ𝑗+1
𝑗-1
𝑗
𝑗+1
(a)
(b)
Ψ𝑗−1
Ψ𝑗
Ψ𝑗+1
(c)
⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟
𝑆(Ψ𝑗−1 )
𝑆(Ψ𝑗+1 )
Figure 2: Reduced-space basis functions in 1D: (a) standard, following Eq. (9), showing overlapping support of Ψ𝑗−1 and Ψ𝑗+1 ; (b) gap-based support; (c) element-centroid-based support. The overlap in case (a) leads to a 5-point stencil while (b) and (c) yield a 3-point stencil such that 𝐴𝑟 is tridiagonal.
Figure 3: Test cases (left to right): T-junction (𝐸=62176) configuration, including coarse-space cells; 146-pebble (𝐸=62138) mesh illustrating an eight-way mesh partition; 45000-pebble (𝐸=13032440) mesh, and 352625-pebble (𝐸=98782067) mesh.
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 15 of 14
Fast Coarse Solvers
Figure 4: Nek5000 test results for different coarse solvers: (top) iteration count-per-step; (bottom) cumulative iteration counts.
Figure 5: NekRS results for different coarse solvers: cumulative iteration counts.
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 16 of 14
Fast Coarse Solvers
Figure 6: Strong scaling study of 45000 pebbles mesh with NekRS. Left to right: Navier-Stokes, Pressure and Coarse grid solve time per timestep for AMG and two level Schwarz method; Parallel efficiency for AMG and two level Schwarz method; Breakdown of the coarse grid solve time for the two level Schwarz method.
Figure 7: Estimated solve times for Schwarz and AMG for different 𝛾 and 𝐿 values.
T. Ratnayaka, P. Fischer, L. Olson: Preprint submitted to Elsevier
Page 17 of 14