Two-Level Neumann–Neumann Preconditioner

• The paper introduces a two-level Neumann–Neumann preconditioner that leverages local Neumann solves and adaptive spectral enrichment to ensure robust convergence for high-contrast elliptic and Maxwell problems.
• It employs overlapping Schwarz methods with localized Neumann-type boundary conditions and a global coarse correction to efficiently reduce errors across subdomains.
• Adaptive coarse spaces, constructed via localized eigenvalue problems on interfaces, yield contrast-independent performance and balance computational cost with convergence speed.

The two-level Neumann–Neumann preconditioner is a scalable iterative preconditioning framework for large, sparse linear systems arising from conforming finite element discretizations of elliptic and Maxwell-type PDEs. It extends classical Schwarz domain decomposition methods by incorporating both local solves on overlapping subdomains with homogeneous Neumann-type boundary conditions and a global coarse space for low-frequency error correction. Adaptive coarse spaces, constructed via local generalized eigenvalue problems on interfaces, are critical for contrast-independent and robust performance, especially in the presence of variable coefficients and non-convex domains. The approach enables fully algebraic construction, making it suitable for problems where geometric or coefficient information is only available through assembled system matrices (Bootland et al., 2020, Heinlein et al., 2022).

1. Formulation of the Model Problem

Consider a domain $\Omega \subset \mathbb{R}^d$ (typically $d=2,3$) discretized by conforming finite elements, such as lowest-order Nédélec edge elements for Maxwell systems or Lagrange P1/Q1 for scalar elliptic equations. The variational problem for the positive Maxwell system is to find $u: \Omega \to \mathbb{R}^3$ such that

$$\begin{aligned} \mathrm{curl} (\mu^{-1} \mathrm{curl}\;u) + \sigma u = f \quad & \text{in } \Omega, \ n \times u = 0 \quad & \text{on } \partial \Omega. \end{aligned}$$

This leads to a symmetric positive-definite (SPD) system $A u = f$ with $A_{ij} = a(\phi_j, \phi_i)$ where $a(u,v)$ incorporates the relevant bilinear forms and $\phi_i$ are basis functions of the finite element space $V_h$.

For scalar elliptic problems, the standard form involves $\alpha(x) \nabla u \cdot \nabla v$, potentially with highly heterogeneous and high-contrast coefficients $\alpha(x)$ (Heinlein et al., 2022).

2. Subdomain Decomposition and Local Neumann Problems

The global domain $\Omega$ is partitioned into overlapping subdomains $\{\Omega_i\}_{i=1}^N$, each with diameter $\sim H$ and overlap width $\delta$. Degrees of freedom (DoFs) are restricted to each subdomain via Boolean restriction operators $R_i$, with $R_i^T$ the natural zero-extension. Local finite element spaces $V_i = R_i V_h$ and inner products induced by restricted stiffness matrices $A_i = R_i A R_i^T$ are defined.

On each $\Omega_i$, local problems enforce Neumann-type boundary conditions:

• For Maxwell: $n \times (\mu^{-1} \mathrm{curl}\;u_i) = 0$ on $\partial \Omega_i \setminus \partial \Omega$, $n \times u_i = 0$ on $\partial \Omega_i \cap \partial \Omega$.
• For scalar: standard homogeneous Neumann or Dirichlet conditions as appropriate.

Local solvers $A_i^{-1}$ appear in the preconditioner as $R_i^T A_i^{-1} R_i$.

3. Construction of the Two-level Coarse Space

The coarse space $Z$ is constructed to ensure robust and scalable convergence. Its structure includes:

a) Near-kernel Space (Gradient Components)

Fields of the form $u = \nabla \psi$ with $\psi \in H_0^1(\Omega)$ span the null-space of $\mathrm{curl}\,\mathrm{curl}$ for Maxwell; in discretization, gradients of P1 scalar FE functions form the basis $G_h$, assembled into $Z_0$.

b) Adaptive Spectral Enrichment

Key enrichment is achieved by solving generalized eigenvalue problems (EVPs) on subdomain interfaces:

• For Maxwell (Bootland et al., 2020): Find $\phi$ supported on $\Gamma_i$ such that

$$\int_{\Gamma_i} \mu^{-1} \mathrm{curl}\;\phi \cdot \mathrm{curl}\;\psi\, dS = \lambda \int_{\Gamma_i} \sigma\; \phi\;\psi\, dS$$

for all $\psi$ supported on $\Gamma_i$. Eigenvectors $\phi_{ik}$ with $\lambda_{ik} \leq \tau$ (threshold) are selected, extended discretely into $\Omega_i$ (Neumann extension), and assembled into $Z_1$.

• For scalar elliptic (Heinlein et al., 2022): Two algebraic EVPs are used per decomposition edge $e$:
  • Dirichlet EVP (AGDSW-type): Modes identified via $S_e \psi = \mu A_{ee} \psi$ where $S_e = A_{ee} - A_{eR} A_{RR}^{-1} A_{Re}$.
  • Transfer EVP (Multiscale/Optimal Local Approximation): Harmonic extension-based operators yield $T^T A_{ee} T v = \lambda M_e v$ with $M_e$ based on minimum coefficients and mesh geometry.

The full coarse space is $Z = [Z_0\ |\ Z_1]$, with dimension $n_0 = N_G + N_{\tau}$.

4. Definition and Application of the Two-level Preconditioner

The two-level Neumann–Neumann (overlapping Schwarz) preconditioner is defined as:

$$M^{-1} = R_0^T A_0^{-1} R_0 + \sum_{i=1}^N R_i^T A_i^{-1} R_i$$

where $A_0 = Z^T A Z$, and the prolongation $R_0^T = Z$. For a residual $r$, application involves:

• Coarse correction: Solve $w_0 = A_0^{-1} (Z^T r)$
• Local corrections: For each $i$, $w_i = A_i^{-1} (R_i r)$
• Update: $M^{-1} r = Z w_0 + \sum_{i=1}^N R_i^T w_i$

During Krylov iteration (CG, GMRES), these steps are performed at each stage, with the coarse space enabling global low-frequency error reduction (Bootland et al., 2020, Heinlein et al., 2022).

5. Spectral Estimates and Condition Number Bounds

Convergence analysis relies on stable decomposition and strengthened Cauchy–Schwarz arguments. For appropriate choices of overlap $\delta$ and coarse mode threshold $\tau$, there exist $c, C$ independent of mesh size $h$ and number of subdomains $N$ so that:

$c\ a(u,u) \leq a(M^{-1} A u, u) \leq C\ a(u,u)$

leading to condition number bound:

$$\kappa(M^{-1}A) \leq C'/c' \leq C'' (1 + \log(H/\delta))^2$$

provided $\tau = O(1)$. The bound is uniform with respect to coefficient jumps and subdomain irregularity if local interface EVPs use local coefficients. For highly heterogeneous elliptic problems, coarse spaces constructed from adaptive eigenvalue problems give contrast-independent bounds (Heinlein et al., 2022).

6. Algorithmic Procedure for Implementation

Algorithmic realization is fully algebraic, requiring only the global assembled matrix $A$:

1. Partition $A$ into subdomain index sets (e.g., METIS), extract local matrices $A_i$.
2. Identify interface nodes, split into edges $e$ and vertices $V$.
3. Build classical GDSW vertex and constant-edge spaces.
4. For each edge:
   • Define oversampling regions, extract local matrices,
   • Solve Dirichlet and transfer EVPs, select modes,
   • Extend interface patterns interior via harmonic extension.
5. Orthogonalize per edge via small POD to remove near dependencies.
6. Assemble coarse prolongation $E_0$, restriction $R_0 = E_0^T$, and coarse matrix $A_0$.
7. In Krylov solves, apply $M^{-1}$ by coarse solve, local solves, and update (Heinlein et al., 2022).

All manipulations (block extraction, extension, eigenmode computation) are executed via index-set and sparse-matrix operations, making geometric or mesh knowledge unnecessary.

7. Remarks on Effectiveness and Flexibility

The near-kernel gradient space $Z_0$ is essential for representing the null-space (e.g., for Maxwell, when $\sigma \to 0$), and for bounding the condition number. Adaptive spectral enrichment via $Z_1$ provides robustness to coefficient jumps and irregular subdomain shapes. The overall computational cost balances coarse-solve size $n_0$ against overlap $\delta$ and desired convergence. Coarse-space enrichment thresholds ($\tau$, $\text{tol}_{dir}$, $\text{tol}_{tr}$) control the trade-off between cost and robustness. This preconditioner is applicable to a wide range of problems, including those with high-contrast or oscillatory coefficients, and in scenarios where mesh or material information is limited to matrix data (Bootland et al., 2020, Heinlein et al., 2022).