Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

A fast spectral overlapping domain decomposition method with discretization-independent conditioning bounds

Simon Dirckx, Anna Yesypenko, Per-Gunnar Martinsson

TL;DR

This work develops a fast, discretization-independent domain decomposition method for general variable-coefficient elliptic PDEs by tessellating the domain into overlapping thin slabs and forming a reduced interface system whose blocks are smooth Dirichlet-to-Dirichlet maps. The method leverages high-order spectral local solves (HPS) and randomized HBS compression to represent dense interface blocks data-sparsely, yielding a well-conditioned equilibrium operator that behaves like a second-kind Fredholm operator. Key contributions include discretization-independent conditioning bounds, robust GMRES performance scaling as O(1/H), and demonstrated capability to solve oscillatory 2D and 3D problems with up to 28 million degrees of freedom. The approach offers scalable parallelism, supports high-order discretizations, and opens pathways to a full direct solver for large-scale elliptic problems, with potential extensions to varied boundary conditions and HPC implementations.

Abstract

A domain decomposition method for the solution of general variable-coefficient elliptic partial differential equations on regular domains is introduced. The method is based on tessellating the domain into overlapping thin slabs or shells, and then explicitly forming a reduced linear system that connects the different domains. Rank-structure ('H-matrix structure') is exploited to handle the large dense blocks that arise in the reduced linear system. Importantly, the formulation used is well-conditioned, as it converges to a second kind Fredholm equation as the precision in the local solves is refined. Moreover, the dense blocks that arise are far more data-sparse than in existing formulations, leading to faster and more efficient H-matrix arithmetic. To form the reduced linear system, black-box randomized compression is used, taking full advantage of the fact that sparse direct solvers are highly efficient on the thin sub-domains. Numerical experiments demonstrate that our solver can handle oscillatory 2D and 3D problems with as many as 28 million degrees of freedom.

A fast spectral overlapping domain decomposition method with discretization-independent conditioning bounds

TL;DR

This work develops a fast, discretization-independent domain decomposition method for general variable-coefficient elliptic PDEs by tessellating the domain into overlapping thin slabs and forming a reduced interface system whose blocks are smooth Dirichlet-to-Dirichlet maps. The method leverages high-order spectral local solves (HPS) and randomized HBS compression to represent dense interface blocks data-sparsely, yielding a well-conditioned equilibrium operator that behaves like a second-kind Fredholm operator. Key contributions include discretization-independent conditioning bounds, robust GMRES performance scaling as O(1/H), and demonstrated capability to solve oscillatory 2D and 3D problems with up to 28 million degrees of freedom. The approach offers scalable parallelism, supports high-order discretizations, and opens pathways to a full direct solver for large-scale elliptic problems, with potential extensions to varied boundary conditions and HPC implementations.

Abstract

A domain decomposition method for the solution of general variable-coefficient elliptic partial differential equations on regular domains is introduced. The method is based on tessellating the domain into overlapping thin slabs or shells, and then explicitly forming a reduced linear system that connects the different domains. Rank-structure ('H-matrix structure') is exploited to handle the large dense blocks that arise in the reduced linear system. Importantly, the formulation used is well-conditioned, as it converges to a second kind Fredholm equation as the precision in the local solves is refined. Moreover, the dense blocks that arise are far more data-sparse than in existing formulations, leading to faster and more efficient H-matrix arithmetic. To form the reduced linear system, black-box randomized compression is used, taking full advantage of the fact that sparse direct solvers are highly efficient on the thin sub-domains. Numerical experiments demonstrate that our solver can handle oscillatory 2D and 3D problems with as many as 28 million degrees of freedom.

Paper Structure

This paper contains 28 sections, 1 theorem, 37 equations, 24 figures, 2 algorithms.

Table of Contents

  1. Introduction
  2. A model problem
  3. Constructing the reduced system
  4. Advantages of the reduced system formulation
  5. Connection to prior work and classical domain decompositions
  6. Derivation of the reduced linear system
  7. A numerical derivation
  8. An analytic derivation
  9. Discretization of the local problem
  10. Rank structure and randomized compression
  11. Summary of the Proposed Method
  12. Condition number estimates
  13. The discrete case
  14. Continuum analysis
  15. GMRES iterations
  16. ...and 13 more sections

Key Result

Theorem 1

Let $\rho(\bm{\mathsf{S}})$ and $\rho(\bm{\mathsf{S}}^{-1})$ denote the spectral radius of the discretized equilibrium operator and the spectral radius of its inverse respectively. Then there is a $C\in\mathbb{R}$, independent of the local slab discretizations, but possibly depending on the positive

Figures (24)

  • Figure 1: Examples of domains that can naturally be tessellated into thin slabs or shells. From left to right, the number of double-wide slabs, $N_{\rm ds}$, is $3$, $2$ and $3$.
  • Figure 2: The model problem considered in Section \ref{['sec:LA']}: A linear system $\bm{\mathsf{A}}$ results from discretization on a computational grid (gray nodes hold Dirichlet data and are not 'active'). The domain is split into thin slabs, separated by the nodes in $J:=\{J_{j}\}_{j=1}^{4}$ (red). The reduced linear system (middle, see (\ref{['eq:modelT']})) has the Schur complement coefficient matrix $\bm{\mathsf{T}} = \bm{\mathsf{A}}(J,J) - \bm{\mathsf{A}}(J,J^c)\,\bm{\mathsf{A}}(J^c,J^c)^{-1}\,\bm{\mathsf{A}}(J^c,J)$. The matrix $\bm{\mathsf{S}}$ from (\ref{['eq:modelS']}) shown on the right is obtained by block diagonal preconditioning of $\bm{\mathsf{T}}$.
  • Figure 3: Continuum domain decomposition described in Section \ref{['sec:analytic']}. (a) BVP on the domain $\Omega$, with known Dirichlet data on $\Gamma$ (green). (b) Local problem on the double-wide strip $\Psi_{j}$, with known Dirichlet data on $\Gamma\cap\partial\Psi_{j}$ (green). The unknown data on $\Gamma_{j-1}$ and $\Gamma_{j+1}$ (blue) is $u_{j-1}$ and $u_{j+1}$, respectively. For fixed $\{u_{j-1},u_{j+1}\}$, there is a unique solution $u_{j}$ on $\Gamma_{j}$ (red).
  • Figure 4: Illustration of the solution operator principle for general configurations, with nonzero boundary conditions (blue) and load (gray).
  • Figure 5: Illustration of the discretization technique described in Section \ref{['sec:HPS']} for discretizing the local boundary value problems introduced in Section \ref{['sec:analytic']}.
  • ...and 19 more figures

Theorems & Definitions (9)

  • Remark 1
  • Remark 2
  • Remark 3
  • Theorem 1
  • proof
  • Remark 4
  • Remark 5
  • Remark 6: Weak scaling in the parallel setting
  • Remark 7: Complexity of the Hierarchical Poincaré–Steklov discretization