A dimension-oblivious domain decomposition method based on space-filling curves

1. Problem Statement

The paper addresses the challenge of solving high-dimensional elliptic partial differential equations (PDEs), specifically the Poisson equation, on modern high-performance computing (HPC) systems. The core difficulties identified are:

• The Curse of Dimensionality: Standard discretizations (full grids) become computationally infeasible as the dimension $d$ increases. The authors utilize the Sparse Grid Combination Technique, which decomposes the high-dimensional problem into many independent, lower-dimensional subproblems on anisotropic grids.
• Scalability and Load Balancing: To utilize exascale systems (millions of cores), the solver must handle a massive number of independent subproblems simultaneously. Each subproblem has a different number of degrees of freedom (DOFs) due to anisotropic mesh sizes. The solver must partition these subproblems efficiently across arbitrary numbers of processors ($P$) regardless of the dimension $d$.
• Fault Tolerance: Future exascale systems require fault tolerance. The authors aim to develop solvers that can recover from hardware failures via data redundancy. This necessitates large overlaps between subdomains, which contradicts traditional domain decomposition (DD) goals of minimizing overlap to reduce communication costs.
• Dimension-Obliviousness: Conventional geometric DD methods rely on fixed-dimensional spatial information, making them difficult to apply to the varying anisotropic grids generated by the combination technique. A method that operates purely algebraically, independent of $d$, is required.

2. Methodology

The authors propose a two-level algebraic domain decomposition solver based on Space-Filling Curves (SFCs).

A. Space-Filling Curve Partitioning

Instead of geometric partitioning, the method maps the $d$-dimensional grid points to a one-dimensional interval $[0, 1]$ using a discrete space-filling curve (specifically the Hilbert curve).

• Mapping: Grid points are ordered by their SFC index.
• Disjoint Partitioning: The ordered list of $N$ points is split into $P$ contiguous segments of roughly equal size ($N/P$), creating disjoint subdomains $\tilde{\Omega}_i$. This ensures perfect load balancing regardless of dimension $d$ or anisotropy.
• Overlap Construction: To achieve data redundancy for fault tolerance, the disjoint sets are enlarged to create overlapping subdomains $\Omega_i$. The overlap is controlled by a parameter $\gamma$ in the 1D parameter space of the SFC.
  • $\Omega_i$ includes $\tilde{\Omega}_i$ plus neighboring sets $\tilde{\Omega}_{i \pm k}$ based on $\gamma$.
  • This approach creates overlaps that are "dimension-oblivious" because they are defined by the 1D curve ordering, not by geometric dimensions.

B. Algebraic Two-Level Preconditioner

The solver uses an additive Schwarz method with a coarse grid correction, constructed purely algebraically:

• Fine Level: Local problems are solved on the overlapping subdomains $\Omega_i$.
• Coarse Level: A coarse space is defined algebraically without geometric coarse grids. For each disjoint subdomain $\tilde{\Omega}_i$, $q$ degrees of freedom are selected. The restriction operator $R_0$ maps fine grid values to these coarse values.
• Weighting: To handle the overcounting of unknowns in overlapping regions, a partition of unity is applied. The authors introduce a specific weighting strategy where, for $\gamma = n/2$, the weights become uniform constants ($\omega_i = 1/(2\gamma+1)$), ensuring the resulting operator is symmetric.
• Implementation: The coarse problem is solved redundantly on every processor (Bank-Holst technique) to avoid communication bottlenecks and enhance fault tolerance.

C. Iterative Solvers

The preconditioner is applied within:

1. Damped Richardson Iteration: A linear iterative scheme.
2. Preconditioned Conjugate Gradient (PCG): For symmetric positive definite systems.

3. Key Contributions

1. Dimension-Oblivious Partitioning: The method successfully partitions anisotropic, high-dimensional grids into balanced subdomains using only 1D SFC ordering, eliminating the need for geometric dimension-dependent logic.
2. Large Overlap for Fault Tolerance: Unlike traditional DD methods that minimize overlap, this method intentionally uses large overlaps (controlled by $\gamma$) to store redundant data, enabling data recovery in the event of processor failures.
3. Algebraic Coarse Grid: The coarse problem is constructed algebraically using the SFC ordering, making it applicable to arbitrary dimensions and anisotropic meshes where geometric coarse grids are ill-defined.
4. Symmetry via Specific Overlap: The authors prove that choosing $\gamma$ as integer multiples of $0.5$(e.g.,$\gamma = 0.5, 1.0$) results in a symmetric preconditioner even with the large overlaps required for fault tolerance.

4. Numerical Results

The authors conducted extensive numerical experiments on dimensions $d=1$ to $6$ using up to one million processors.

• Convergence Behavior:

  • The solver exhibits optimal weak scaling: The number of iterations required to converge remains constant as the number of processors $P$ increases (keeping the subdomain size fixed).
  • Dimension Independence: The convergence rate is largely independent of the dimension $d$. In fact, the 1D case was the most difficult (requiring the most iterations), while higher dimensions ($d \ge 6$) converged faster or similarly, stabilizing at a low iteration count.
  • PCG vs. Richardson: The Preconditioned Conjugate Gradient method significantly outperformed the Richardson iteration, reducing iteration counts by a factor of roughly $\sqrt{\kappa}$.
  • Overlap Parameter: An overlap parameter of $\gamma = 0.5$ provided the best balance between convergence speed and redundancy. Larger $\gamma$ improved fault tolerance but slightly increased iteration counts (unless weighted correctly).

• Scaling:

  • Weak Scaling: The method maintained constant iteration counts for $P$ up to $256$ subdomains per subproblem, and the overall combination technique utilized over 1 million processors effectively.
  • Strong Scaling: The method showed expected behavior, with iteration counts decreasing as $P$ increased (until the subdomain size became too small relative to the coarse grid).

• Combination Technique Application:

  • When applied to the full sparse grid combination method for $d=2$ to $6$, the solver successfully managed the massive parallelism (millions of tasks).
  • Runtimes per subproblem remained consistent across different dimensions and discretization levels, demonstrating that the algorithm's performance depends only on the local problem size, not the global dimension.

5. Significance

This paper presents a foundational step toward exascale computing for high-dimensional PDEs.

• Scalability: It demonstrates that sparse grid combination techniques can be scaled to utilize millions of cores efficiently, overcoming the limitations of traditional solvers.
• Fault Tolerance: By integrating large overlaps directly into the algebraic solver structure, it paves the way for resilient computing on unreliable hardware, a critical requirement for future exascale systems.
• Generality: The "dimension-oblivious" nature of the solver makes it applicable to a wide range of problems involving anisotropic meshes and high dimensions, where standard geometric domain decomposition fails.

In summary, the authors have developed a robust, algebraic domain decomposition solver that leverages space-filling curves to achieve optimal convergence and extreme scalability in arbitrary dimensions, specifically tailored for the demands of fault-tolerant exascale computing.