# A GenEO Domain Decomposition method for Saddle Point problems

Frédéric Nataf, Pierre-Henri H Tournier

► To cite this version:

Frédéric Nataf, Pierre-Henri H Tournier. A GenEO Domain Decomposition method for Saddle Point problems. Comptes Rendus. Mécanique, 2023, 351 (S1), pp.1-18. 10.5802/crmeca.175. hal-02343808v6

HAL Id: hal-02343808

https://hal.science/hal-02343808v6

Submitted on 16 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A GenEO Domain Decomposition method
for Saddle Point problems

F. Nataf *¹ and P.H. Tournier †²

¹,²Sorbonne Université, CNRS, Université de Paris, Inria Equipe Alpines, Laboratoire Jacques-Louis Lions, F-75005 Paris, France,

November 8, 2022

Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Preconditioning of the primal problem</td><td>3</td></tr><tr><td>3</td><td>Schur complement preconditioning</td><td>4</td></tr><tr><td>3.1</td><td>First spectrally equivalent preconditioner</td><td>4</td></tr><tr><td>3.2</td><td>Preconditioning of S<sub>1</sub></td><td>6</td></tr><tr><td>3.2.1</td><td>One-level DD for S<sub>1</sub></td><td>6</td></tr><tr><td></td><td>Surjectivity of R</td><td>7</td></tr><tr><td></td><td>Continuity of R</td><td>7</td></tr><tr><td></td><td>Stable decomposition</td><td>7</td></tr><tr><td>3.2.2</td><td>Two-level DD for S<sub>1</sub></td><td>8</td></tr><tr><td>3.3</td><td>Final Preconditioner for the Schur complement</td><td>9</td></tr><tr><td>4</td><td>Recap</td><td>9</td></tr><tr><td>4.1</td><td>Setup for the Schur complement preconditioner</td><td>9</td></tr><tr><td>4.2</td><td>DD solver for the saddle point system</td><td>9</td></tr><tr><td>5</td><td>Numerical experiments</td><td>10</td></tr><tr><td>5.1</td><td>Software, hardware, implementation details</td><td>10</td></tr><tr><td>5.2</td><td>Parameters of the method</td><td>12</td></tr><tr><td>5.3</td><td>Weak scalability test for heterogeneous steel and rubber beam</td><td>12</td></tr><tr><td></td><td>Iteration counts.</td><td>13</td></tr><tr><td></td><td>Timings.</td><td>13</td></tr><tr><td>5.4</td><td>Strong scalability test for heterogeneous steel and rubber beam</td><td>13</td></tr><tr><td></td><td>Iteration counts.</td><td>13</td></tr><tr><td></td><td>Timings.</td><td>14</td></tr><tr><td>5.5</td><td>Comparison with other algorithms</td><td>14</td></tr><tr><td>5.5.1</td><td>Comparison with a direct solver.</td><td>14</td></tr><tr><td>5.5.2</td><td>Comparison with other Schwarz methods</td><td>15</td></tr><tr><td>5.5.3</td><td>Comparison with algebraic multigrid method</td><td>15</td></tr></table>

<footnote>*frederic.nataf@sorbonne-universite.fr</footnote>

<footnote>†pierre-henri.tournier@sorbonne-universite.fr</footnote>
Abstract

We introduce an adaptive element-based domain decomposition (DD) method for solving saddle point problems defined as a block two by two matrix. The algorithm does not require any knowledge of the constrained space. We assume that all sub matrices are sparse and that the diagonal blocks are spectrally equivalent to a sum of positive semi definite matrices. The latter assumption enables the design of adaptive coarse space for DD methods that extends the GenEO theory [40] to saddle point problems. Numerical results on three dimensional elasticity problems for steel-rubber structures discretized by a finite element with continuous pressure are shown for up to one billion degrees of freedom.

# 1 Introduction

Solving saddle point problems with parallel algorithms is very important for many branches of scientific computing: fluid and solid mechanics especially for incompressible or nearly incompressible materials, computational electromagnetism, multi-physics problems, inverse problems and optimization. They also arise when analyzing realistic problems with the incorporation of a set of linear multipoint constraints (MPC) into the linear solver. Examples of the problems described by the MPCs are the non-penetration condition in the contact problem, connection of different types of finite elements, modelling of rigid bodies, and so on. Using Lagrange multipliers to enforce them, yields linear saddle point systems.

For small enough problems, direct solvers are the method of choice since they are robust and have a predictable efficiency. But for large problems the memory requirements are very large and at some point the solver runs out of memory, see e.g. the tests performed in § 5.5.1. Then iterative solvers are needed. Domain decomposition methods are, alongside multigrid methods, one of the dominant paradigms for defining efficient and robust preconditioners in modern large-scale applications dealing with partial differential equations and algebraic systems of equations.

We are interested in domain decomposition (DD) methods since they are naturally well-fitted to modern parallel architectures. For specific systems of partial differential equations with a saddle point formulation, efficient DD methods have been designed, see e.g., [36, 31, 37] and [42] references therein. Also in [23], a GenEO coarse space is introduced for the P.L. Lions' algorithm and its efficiency is mathematically proved for symmetric definite positive problems. In the above article, numerical experiments are conducted on three dimensional elasticity problems for steel-rubber structures discretized by a finite element with continuous pressure. Although the method works well in practice, it lacks theoretical convergence guarantees and also demands the design of specific absorbing conditions as interface conditions.

As for a convergence rate analysis for a discretization with a continuous pressure, the recent article [45] generalizes the theory developed in [44] to the case of nonzero pressure block but under the assumption that the discontinuities are resolved by the subdomains. Compared to the above mentioned works, the method we propose has a provable control on the condition number for zero or non zero pressure block with a continuously discretized pressure also in the case arbitrary heterogeneities and, contrarily to [23], bypasses the need for absorbing boundary conditions.

Here as in [34, 5, 11, 38], we consider the problem in the form of a two by two block matrix. Let $m$ and $n$ be two integers with $m < n$. Let $A n \times n$ be a SPD matrix, $B$ a sparse $m \times n$ full rank matrix of constraints and $C$ a $m \times m$ non negative matrix (in particular, $C = 0$ is allowed), we consider the following saddle point matrix:

$$ \mathcal{A} := \begin{pmatrix} A & B^T \\ B & -C \end{pmatrix}. \qquad (1) $$
When the kernel of matrix B is known, very efficient multigrid methods have been designed in the context of finite element methods, see e.g., [9, 26, 25, 2, 39, 16]. Without this knowledge, it is nevertheless also possible to design efficient geometric multigrid methods as in [13] where the fine mesh is obtained by several uniform mesh refinements. In the case where matrix C is zero, augmented Lagrangian methods [18, 20, 17] may be used. Applications to incompressible fluid dynamics and incompressible elasticity are numerous, see e.g. in [16, 19] and references therein.

Here we do not assume any knowledge on the kernel of matrix B nor that C is zero and we work with arbitrary meshes. The following three factor factorization, see e.g., [6]:

$$ \begin{pmatrix} A & B^T \\ B & -C \end{pmatrix} = \begin{pmatrix} I & 0 \\ BA^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & -(C + BA^{-1}B^T) \end{pmatrix} \begin{pmatrix} I & A^{-1}B^T \\ 0 & I \end{pmatrix}, $$

shows that solving the linear system with A can be performed by solving sequentially linear systems with A and one with the Schur complement $C + BA^{-1}B^T$. In order to build a scalable method, we assume that all three matrices A, B and C are sparse and that A and C are the sum of positive semi definite matrices. This is easily achieved in finite element or finite volume contexts for partial differential equations. The latter assumption enables the design of adaptive coarse space for DD methods, see [12].

The highlights of the article are:

* new provably scalable DD method that is a blend of two-level Schwarz and Neumann-Neumann domain decomposition methods.

* successfully applied to the notorious difficult problem of large scale (1 billion d.o.f's) nearly incompressible elasticity with heterogeneities not resolved by the domain decomposition and discretized with a continuous pressure.

The paper is organized as follows. In § 2, we recall the two-level additive Schwarz method denoted $M_A$ used to precondition the matrix A (the primal-primal block of the saddle point problem). Then in § 3, we introduce the operator $P_S := C + BM_A^{-1}B^T$ which is spectrally equivalent to the Schur complement S. Its preconditioning is studied in § 3.2. In § 4, we combine these different components to define in a compact way the parallel saddle point preconditioner. In § 5, we present weak and strong scaling experiments on large scale elasticity problems for steel-rubber structures discretized by a finite element with continuous pressure. Comparisons with direct and multigrid solvers are also given. These problems are highly heterogeneous since the Lamé-Poisson coefficients of the rubber are $(E_1, \nu_1) = (1 \times 10^7, 0.4999)$ and those of the steel are $(E_2, \nu_2) = (2 \times 10^9, 0.35)$.

## 2 Preconditioning of the primal problem

The sparse $n \times n$ SPD matrix A is preconditioned by a two-level Schwarz type DD method:

$$ M_A^{-1} := R_0^T (R_0 A R_0^T)^{-1} R_0 + \sum_{i=1}^{N} R_i^T (R_i A R_i^T)^{-1} R_i, \quad (2) $$

where $R_0$ is full rank $dim(V_0) \times n$ with $V_0$ denoting the space spanned by the columns of $R_0^T$. The following assumptions are crucial to ensure the final method is scalable:

### Assumption 2.1 (dimension and structure of the coarse space)

* The coarse space dimension, $dim(V_0)$, is $O(N)$ typically 10-20 times N.

* The coarse space is made of extensions by zero of local vectors.
Using the GenEO method [40], it is possible to specify in advance two constants $0 < \lambda_m < 1 < \lambda_M$ and then build a coarse space $V_0$ such that $M_A^{-1}$ is spectrally equivalent to $A^{-1}$:

$$ \frac{1}{\lambda_M} M_A^{-1} \le A^{-1} \le \frac{1}{\lambda_m} M_A^{-1}. \qquad (3) $$

The dimension of the coarse space $V_0$ is typically proportional to the number of subdomains. This corresponds to Assumption 2.1. More precisely, for each subdomain $1 \le i \le N$, let $D_i$ be a non negative diagonal matrix that defines a discrete partition of unity, i.e.:

$$ \sum_{i=1}^{N} R_{i}^{T} D_{i} R_{i} = I, $$

and $A_i^{Neu}$ be a symmetric semi-definite positive matrix such that for the maximum multiplicity of the intersection of subdomains denoted $k_1$, we have:

$$ \sum_{i=1}^{N} R_{i}^{T} A_{i}^{Neu} R_{i} \le k_{1} A. \qquad (4) $$

Then, the GenEO eigenvalue problem is local to each subdomain and reads: Find $(\lambda_{ik}, V_{ik}) \in \mathbb{R} \times \mathbb{R}^{rank(R_i)}$ such that:

$$ (D_i R_i A R_i^T D_i) V_{ik} = \lambda_{ik} A_i^{Neu} V_{ik}. \qquad (5) $$

Let $\tau > 0$ be a positive threshold, the coarse space is the vector space spanned by the vectors $R_i^T D_i V_{ik}$ for all $\lambda_{ik} > \tau$. Then inequality (3) holds with $\lambda_m := (1 + k_1 \tau)^{-1}$ and $\lambda_M := k_0$ where $k_0$ is the maximal number of neighbours of a subdomain including itself.

Our aim in the next section is to precondition the Schur complement $-S$ of matrix $A$ (eq. (1)) where

$$ S := C + B A^{-1} B^{T}. \qquad (6) $$

This is achieved by a series of spectrally equivalent matrices or preconditioners, the first one being $P_S$ defined as follows:

$$ P_S := C + B M_A^{-1} B^T \qquad (7) $$

and the final one being $N_S^{-1}$ introduced in § 3.3, see eq. (24). Finally in § 4.2 we will introduce the preconditioner of the saddle point matrix $\mathcal{A}$.

# 3 Schur complement preconditioning

## 3.1 First spectrally equivalent preconditioner

Note that $P_S$ is by definition a sum of $N + 2$ positive semi definite matrices

$$ P_S := B R_0^T (R_0 A R_0^T)^{-1} R_0 B^T + C + \sum_{i=1}^{N} B R_i^T (R_i A R_i^T)^{-1} R_i B^T. \qquad (8) $$

Since $B$ is a sparse matrix, it is interesting to introduce, for all $0 \le i \le N$, $\tilde{R}_i$ the restriction operator on the support of $\mathfrak{S}(B R_i^T)$ (the image of $B R_i^T$) so that:

$$ \tilde{R}_i^T \tilde{R}_i B R_i^T = B R_i^T. \qquad (9) $$

Then, by defining for $0 \le i \le N$,

$$ \tilde{B}_i := \tilde{R}_i B R_i^T, $$
the operator $P_S$ is rewritten as

$$P_S := \tilde{R}_0^T \tilde{B}_0 (R_0 A R_0^T)^{-1} \tilde{B}_0^T \tilde{R}_0 + C + \sum_{i=1}^{N} \tilde{R}_i^T \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T \tilde{R}_i. \quad (10)$$

We consider a partition of unity on $H := \mathbb{R}^m$ defined with local diagonal matrices $(\tilde{D}_i)_{1\le i\le N} \in \mathbb{R}^{\dim(\Im(B R_i^T)) \times \dim(\Im(B R_i^T))}$:

$$\sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i \tilde{R}_i = I_H.$$

We make the following assumption

**Assumption 3.1** There exist symmetric positive semidefinite matrices $(\tilde{C}_i)_{1\le i\le N}$ such that for some constant $\tilde{k}_1$

$$C \le \sum_{i=1}^{N} \tilde{R}_i^T \tilde{C}_i \tilde{R}_i \le \tilde{k}_1 C. \quad (11)$$

This assumption is not so restrictive. Indeed, for a minimization problem with constraints enforced exactly without penalization nor relaxation, we have $C = 0$ and the assumption is automatically satisfied. Moreover, we have:

**Lemma 3.1** If $C$ is a diagonal matrix, Assumption 3.1 is satisfied with $\tilde{k}_1 = 1$.

**Proof** If $C$ is a diagonal matrix, it suffices to take

$$\tilde{C}_i := \tilde{R}_i C \tilde{R}_i^T \tilde{D}_i,$$

which is a diagonal non negative matrix. Indeed, we then have:

$$C = \sum_{i=1}^{N} C \tilde{R}_i^T \tilde{D}_i \tilde{R}_i = \sum_{i=1}^{N} \tilde{R}_i^T (\tilde{R}_i C \tilde{R}_i^T \tilde{D}_i) \tilde{R}_i.$$

<figure><img src="image_2.png" /></figure>

**Remark 3.1** Note also that in the finite element case it suffices to restrict the variational form that defines $C$ to the subdomains. In this case $\tilde{k}_1$ is the multiplicity of the intersections of the subdomains used to define the $\tilde{C}_i$'s.

Let us define the operator $M_S$ as the sum of a non local but low rank matrix $S_0$:

$$S_0 := \tilde{R}_0^T \tilde{B}_0 (R_0 A R_0^T)^{-1} \tilde{B}_0^T \tilde{R}_0,$$

and of $S_1$ which is a sum of $N$ local positive semi definite matrices:

$$S_1 := \sum_{i=1}^{N} \tilde{R}_i^T (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \tilde{R}_i,$$

that is

$$M_S := S_0 + S_1.$$

By Assumption 3.1, the operator $M_S$ is spectrally equivalent to $P_S$ which is also spectrally equivalent to $S$. Note that we may assume that $S_1$ is invertible whereas it does not make sense for $S_0$. Note that if it is not the case, since we build a preconditioner, $S_1$ can be regularized by a small diagonal term with little effect on the efficiency of the preconditioner.

We consider next the construction of a preconditioner $M_{S_1}^{-1}$ to $S_1$ leveraging the fact that $S_1$ is a sum of symmetric semidefinite positive matrices. Let us stress that this property stems from the domain decomposition structure of the preconditioner for matrix $A$ which apart from the coarse level is block diagonal.
## 3.2 Preconditioning of $S_1$

It is well known that one level domain decomposition methods are in most cases not scalable. Nevertheless, the study of a one-level method in § 3.2.1 enables the identification of a suitable coarse space that will be efficiently embedded in a scalable two-level domain decomposition method in § 3.2.2.

Our studies of the spectrum of the DD preconditioners are based on the Fictitious Space lemma which is recalled here, see [35] for the original paper and [21] for a modern presentation.

**Lemma 3.2 (Fictitious Space Lemma, Nepomnyaschikh 1991)** Let $H$ and $H_D$ be two Hilbert spaces, with the scalar products denoted by $(\cdot, \cdot)$ and $(\cdot, \cdot)_D$. Let the symmetric positive bilinear forms $a: H \times H \to \mathbb{R}$ and $b: H_D \times H_D \to \mathbb{R}$, generated by the s.p.d. operators $A: H \to H$ and $B: H_D \to H_D$, respectively (i.e. $(Au, v) = a(u, v)$ for all $u, v \in H$ and $(Bu_D, v_D)_D = b(u_D, v_D)$ for all $u_D, v_D \in H_D$). Suppose that there exists a linear operator $\mathcal{R}: H_D \to H$ that satisfies the following three assumptions:

(i) $\mathcal{R}$ is surjective.

(ii) Continuity of $\mathcal{R}$: there exists a positive constant $c_R$ such that

$$a(\mathcal{R}u_D, \mathcal{R}u_D) \le c_R \cdot b(u_D, u_D) \quad \forall u_D \in H_D. \qquad (12)$$

(iii) Stable decomposition: there exists a positive constant $c_T$ such that for all $u \in H$ there exists $u_D \in H_D$ with $\mathcal{R}u_D = u$ and

$$c_T \cdot b(u_D, u_D) \le a(\mathcal{R}u_D, \mathcal{R}u_D) = a(u, u). \qquad (13)$$

We introduce the adjoint operator $\mathcal{R}^*$: $H \to H_D$ by $(\mathcal{R}u_D, u) = (u_D, \mathcal{R}^*u)_D$ for all $u_D \in H_D$ and $u \in H$.

Then, we have the following spectral estimate

$$c_T \cdot a(u, u) \le a(\mathcal{R}B^{-1}\mathcal{R}^*Au, u) \le c_R \cdot a(u, u), \quad \forall u \in H, \qquad (14)$$

which proves that the eigenvalues of operator $\mathcal{R}B^{-1}\mathcal{R}^*A$ are bounded from below by $c_T$ and from above by $c_R$.

The Fictitious Space Lemma (FSL) can also conveniently be related to the book [42]: the first assumption corresponds to equation (2.3), page 36 where the global Hilbert space is assumed to satisfy a decomposition into subspaces, the second assumption is related to Assumptions 2.3 and 2.4, page 40 and the third assumption corresponds to the Stable decomposition Assumption 2.2 page 40.

### 3.2.1 One-level DD for $S_1$

As in [12] chapter 7, we begin with a one-level Neumann-Neumann [7] type DD method defined in terms of the Fictitious Space Lemma (FSL). This study will be the basis for constructing a scalable two-level preconditioner in § 3.2.2. Recall the formula for the one-level preconditioner $M_{S_1,1}^{-1}$ for $S_1$:

$$M_{S_1,1}^{-1} := \sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T)^{\dagger} \tilde{D}_i \tilde{R}_i, \qquad (15)$$

where the superscript $\dagger$ denotes a pseudo inverse in case the operator in brackets is not invertible. For sake of simplicity, we assume that they are invertible so that the following framework enables the study of $M_{S_1,1}$ with the fictitious space lemma. Let

$$H := \mathbb{R}^m$$
and let $a$ be the following bilinear form:

$$a: H \times H \to \mathbb{R} \quad a(\mathbf{P}, \mathbf{Q}) := (S_1\mathbf{P}, \mathbf{Q}).$$

Let

$$H_D := \prod_{i=1}^{N} \mathbb{R}^{\mathrm{rank}(\tilde{B}_i)},$$

and $b$ be the following bilinear form:

$$b: H_D \times H_D \to \mathbb{R} \quad b((\mathbf{P}_i)_{1\le i \le N}, (\mathbf{Q}_i)_{1\le i \le N}) := \sum_{i=1}^{N} ((\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \mathbf{P}_i, \mathbf{Q}_i). \qquad (16)$$

We define $\mathcal{R}$:

$$\mathcal{R}: \quad H_D \to H \\ (\mathbf{P}_i)_{1\le i \le N} \mapsto \sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i \mathbf{P}_i.$$

We now check the three assumptions of the FSL.

**Surjectivity of $\mathcal{R}$** For any $\mathbf{P} \in H$, we have using the partition of unity that:

$$\mathbf{P} = \sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i \tilde{R}_i \mathbf{P},$$

so that

$$\mathbf{P} = \mathcal{R}((\tilde{R}_i \mathbf{P})_{1\le i \le N}). \qquad (17)$$

**Continuity of $\mathcal{R}$** On one hand, we have using $k_0$ the number of neighbours of a subdominant (including itself), $k_0 := \max_{1\le i \le N} \#\mathcal{O}(i)$ where $\mathcal{O}(i) := \{1 \le j \le N \mid \tilde{R}_i \tilde{D}_i S_1 \tilde{D}_j \tilde{R}_j^T \ne 0\}$:

$$a(\mathcal{R}(\mathbf{P}), \mathcal{R}(\mathbf{P})) = \| (\sum_{i=1}^N \tilde{R}_i^T \tilde{D}_i \mathbf{P}_i) \|_a^2 \le k_0 \sum_{i=1}^N \| \tilde{R}_i^T \tilde{D}_i \mathbf{P}_i \|_a^2 \\ \le k_0 \left( (\sum_{j \in \mathcal{O}(i)} \tilde{R}_j^T (\tilde{C}_j + \tilde{B}_j (R_j A R_j^T)^{-1} \tilde{B}_j^T) \tilde{R}_j) \tilde{R}_i^T \tilde{D}_i \mathbf{P}_i, \tilde{R}_i^T \tilde{D}_i \mathbf{P}_i \right).$$

On the other hand, we have by definition:

$$b(\mathcal{P}, \mathcal{P}) = \sum_{i=1}^{N} ((\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \mathbf{P}_i, \mathbf{P}_i).$$

We can take:

$$c_R := \max_{1 \le i \le N} \max_{\mathbf{P}_i \in \mathbb{R}^{\mathrm{rank}(\tilde{B}_i)}} \frac{(\sum_{j \in \mathcal{O}(i)} \tilde{R}_i \tilde{R}_j^T (\tilde{C}_j + \tilde{B}_j (R_j A R_j^T)^{-1} \tilde{B}_j^T) \tilde{R}_j \tilde{R}_i^T \tilde{D}_i \mathbf{P}_i, \tilde{D}_i \mathbf{P}_i)}{((\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \mathbf{P}_i, \mathbf{P}_i)}, \qquad (18)$$

but we have no control on $c_R$ which may be large. This motivates the introduction of a spectral coarse space in § 3.2.2 with the generalized eigenvalue problem (19).

**Stable decomposition** Let $\mathbf{P} \in H$, we start from its decomposition (17) and estimate its $b$-norm

$$b(\mathcal{P}, \mathcal{P}) = \sum_{i=1}^{N} ((\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \tilde{R}_i \mathbf{P}, \tilde{R}_i \mathbf{P}) = a(\mathbf{P}, \mathbf{P}),$$

so that we can take $c_T = 1$ which is an optimal value.
### 3.2.2 Two-level DD for $S_1$

In order to control the value of $c_R$ defined above, we introduce a two-level preconditioner similarly to what is done for Schur complement methods in [12, § 7.8.3, page 197] or in [41]. The generalized eigenvalue problem in each subdomain $1 \le i \le N$ to be solved to build the coarse space is inferred from the definition of the constant $c_R$ in eq. (18):

$$ \begin{aligned} \tilde{D}_i \tilde{R}_i \left( \sum_{j \in O(i)} \tilde{R}_j^T (\tilde{C}_j + \tilde{B}_j (R_j A R_j^T)^{-1} \tilde{B}_j^T) \tilde{R}_j \right) \tilde{R}_i^T \tilde{D}_i \mathbf{P}_{ik} \right) &= \lambda_{ik} (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \mathbf{P}_{ik}. \end{aligned} \quad (19) $$

This generalized eigenvalue problem contains inverses of some local matrices on both sides and in order to solve it via e.g. Arpack, we have to factorize one of them. This difficulty can be overcome since the right matrix is the inverse of the Schur complement of the extended sparse matrix (23). It is thus amenable to a factorization using only sparse matrix factorizations. Note that matrix-vector products on the left hand side involve neighbor to neighbor communications and thus synchronization points. Consequently, these generalized eigenvalue problem solves are a little more involved to code than in the standard GenEO.

The coarse space is defined as follows. Let $\tau_{S_1}$ be a user-defined threshold; for each subdomain $1 \le i \le N$, we introduce a subspace $\tilde{W}_i \subset \mathbb{R}^{rank(\tilde{B}_i)}$:

$$ \tilde{W}_i := \operatorname{Span}\{\mathbf{P}_{ik} \mid \lambda_{ik} > \tau_{S_1}\}. \quad (20) $$

Then the coarse space $\tilde{W}_0$ is defined (with some abuse of notation) as

$$ \tilde{W}_0 := \bigoplus_{1 \le i \le N} \tilde{R}_i^T \tilde{D}_i \tilde{W}_i. $$

Let $Z_{S_1}$ be a rectangular matrix whose columns span the coarse space $\tilde{W}_0$. Let $\tilde{P}_0$ be the $S_1$ orthogonal projection from $\mathbb{R}^m$ on $\tilde{W}_0$ whose formula is

$$ \tilde{P}_0 = Z_{S_1} (Z_{S_1}^T S_1 Z_{S_1})^{-1} Z_{S_1}^T S_1. \quad (21) $$

In order to avoid a too cumbersome analysis, we make the following assumption:

**Assumption 3.2** We assume that for all subdomains $1 \le i \le N$, $\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T$ is invertible.

Finally, the preconditioner for $S_1$ reads

$$ \begin{aligned} M_{S_1}^{-1} &:= Z_{S_1} (Z_{S_1}^T S_1 Z_{S_1})^{-1} Z_{S_1}^T + (I - \tilde{P}_0) \\ &\times \left( \sum_{i=1}^N \tilde{R}_i^T \tilde{D}_i (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T)^\dagger \tilde{D}_i \tilde{R}_i \right) (I - \tilde{P}_0^T). \end{aligned} \quad (22) $$

If Assumption 3.2 is not satisfied for some subdomain $i$, we should incorporate the kernel of $\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T$ into the coarse space and make use of a pseudo inverse in the definition of the preconditioner as it is done for the FETI method, see [15] or [12] § 7.8.2 and references therein. Recall that from [12] chapter 7, we have for $\alpha := \max(1, \frac{k_0}{\tau_{S_1}})$:

$$ \frac{1}{\alpha} M_{S_1}^{-1} \le S^{-1} \le M_{S_1}^{-1}. $$

A careful implementation of (22) requires two coarse solves. Note that the application of $M_{S_1}^{-1}$ can be done using only sparse solvers since solving a linear system with a local Schur complement

$$ (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \mathbf{P}_i = \mathbf{G}_i, $$

amounts to solving an augmented sparse system which has the form of a local saddle point system:

$$ - \begin{pmatrix} R_i A R_i^T & \tilde{B}_i^T \\ \tilde{B}_i & -\tilde{C}_i \end{pmatrix} \begin{pmatrix} \mathbf{U}_i \\ \mathbf{P}_i \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ \mathbf{G}_i \end{pmatrix}. \quad (23) $$
### 3.3 Final Preconditioner for the Schur complement

From the spectrally equivalent preconditioner $M_{S_1}$ to $S_1$, we define $N_S$ a spectrally equiv- alent preconditioner to $M_S$ and thus to $S$ as well:

$$N_S := S_0 + M_{S_1}. \qquad (24)$$

We now consider the application of the preconditioner $N_S$, that is for some right hand side $G \in \mathbb{R}^m$, the solving in $P$ of the following system:

$$N_S P = G, \qquad (25)$$

by a Krylov solver with $M_{S_1}^{-1}$ as a preconditioner.

## 4 Recap

### 4.1 Setup for the Schur complement preconditioner

We have a setup phase which is composed of:

1. Build the two-level preconditioner $M_A^{-1}$ for $A$, see eq. (2),

2. Build the two-level preconditioner $M_{S_1}^{-1}$ for $S_1$, see eq. (22).

Once the setup is complete, applying preconditioner $N_S^{-1}$ can be performed following Algorithm 1

**Algorithm 1** $N_S^{-1}$ matvec product

INPUT: $G \in \mathbb{R}^m$ OUTPUT: $P = N_S^{-1} G$

1. Solve eq. (25) in $P$ by a Krylov method with $M_S^{-1}$ as preconditioner.

### 4.2 DD solver for the saddle point system

**Algorithm 2** DD saddle point solver

INPUT: $(\mathbf{F}_U, \mathbf{F}_P) \in \mathbb{R}^{n+m}$ OUTPUT: $(\mathbf{U}, \mathbf{P})$ the solution to (26).

1. Solve $A\mathbf{G}_U = \mathbf{F}_U$ by a PCG with $M_A^{-1}$ as a preconditioner

2. Compute $\mathbf{G}_P := \mathbf{F}_P - B\mathbf{G}_U$

3. Solve $(C + BA^{-1}B^T)\mathbf{P} = -\mathbf{G}_P$ by a PCG with $N_S^{-1}$ as a preconditioner, see Algorithm 1

4. Compute $\mathbf{G}_U := \mathbf{F}_U - B^T\mathbf{P}$

5. Solve $A\mathbf{U} = \mathbf{G}_U$ by a PCG with $M_A^{-1}$ as a preconditioner

The solution to the saddle point problem:

$$\begin{pmatrix} A & B^T \\ B & -C \end{pmatrix} \begin{pmatrix} \mathbf{U} \\ \mathbf{P} \end{pmatrix} = \begin{pmatrix} \mathbf{F}_U \\ \mathbf{F}_P \end{pmatrix} \qquad (26)$$

is obtained by Algorithm 2. Note that Step 3 demands the matrix-vector product with matrix $C + BA^{-1}B^T$ which is done by an iterative solve for matrix $A^{-1}$.
# 5 Numerical experiments

In this section, we perform 3D experiments to illustrate the theory and the performance of the method. We are interested in a heterogeneous elasticity problem with nearly incompressible material, typically rubber-steel structures. First, we recall the definition of the coefficients and the corresponding variational formulation.

The mechanical properties of a solid are characterized by its Young modulus $E$ and Poisson ratio $\nu$, or alternatively by its Lamé coefficients $\lambda$ and $\mu$. They verify the following relations:

$$ \lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)} \quad \text{and} \quad \mu = \frac{E}{2(1 + \nu)}. \qquad (27) $$

For the discretization, we choose a continuous pressure space and take the lowest order Taylor-Hood finite element $C0P2 - C0P1$ whose stability is proved, see e.g., [8]. The domain $\Omega$ is a beam and the variational problem consists in finding $(u_h, p_h) \in \mathcal{V}_h := \mathbb{P}_2^3 \cap (H_0^1(\Omega))^3 \times \mathbb{P}_1$ with Dirichlet boundary conditions on the four lateral faces and Neumann boundary conditions on the other two faces such that for all $(v_h, q_h) \in \mathcal{V}_h$:

$$ \left\{ \begin{array}{ll} \int_{\Omega} 2\mu\underline{\varepsilon}(\mathbf{u}_h) : \underline{\varepsilon}(\mathbf{v}_h)dx & - \int_{\Omega} p_h \mathrm{div}(\mathbf{v}_h)dx = \int_{\Omega} \mathbf{f}\mathbf{v}_h dx \\ - \int_{\Omega} \mathrm{div}(\mathbf{u}_h)q_h dx & - \int_{\Omega} \frac{1}{\lambda}p_h q_h = 0. \end{array} \right. \qquad (28) $$

Letting $\mathbf{u}$ denote the degrees of freedom of $\mathbf{u}_h$ and $p$ that of $p_h$, the problem can be written in matrix form as:

$$ \left( \begin{array}{cc} A & B^T \\ B & -C \end{array} \right) \left( \begin{array}{c} \mathbf{u} \\ p \end{array} \right) = \left( \begin{array}{c} f \\ 0 \end{array} \right). \qquad (29) $$

The matrix $C$ is a mass matrix arising from the discretization of a variational form. This enables us to satisfy Assumption 3.1 with $\tilde{C}_i$ the corresponding mass matrix but only defined on subdomain $\tilde{\Omega}_i$ which is the extension by a layer of direct neighbors of subdomain $\Omega_i$.

In the following numerical experiments, we consider a heterogeneous beam composed of 10 alternating layers of rubber material $(E_1, \nu_1) = (1 \times 10^7, 0.4999)$ and steel material $(E_2, \nu_2) = (2 \times 10^9, 0.35)$, see Fig. 1. Note that the discontinuities are not resolved by the domain decomposition.

## 5.1 Software, hardware, implementation details

In the following numerical experiments, the iteration counts assess the control of the condition number via the adaptive coarse spaces. We also report timings along with iteration counts. We also showcase the fact that the size of the GenEO coarse space adapts automatically to the difficulty of the problem at hand, for example when going from a homogeneous to a heterogeneous problem. We illustrate the efficiency of the method by performing weak and strong scalability tests, using the automatic graph partitioner Metis [30] for the subdomain partitioning.

The problem is discretized and solved with the open-source parallel finite element software FreeFEM [24]. FreeFEM is a domain specific language (DSL) where the problem to be solved is defined in terms of its variational formulation. Then the local matrices $(A_i^{Neu})_{1 \le i \le N}$ (see eq. (4)) and $(\tilde{C}_i)_{1 \le i \le N}$ (see eq. (11)) are easily obtained by restricting the corresponding variational formulations to adequate local subdomains. Note that these matrices are different from the restriction of the global matrices $A$ and $C$ to the local degrees
of freedom. The domain decomposition algorithm presented in this paper is implemented on top of the *ffddm* framework, a set of parallel *FreeFEM* scripts implementing Schwarz domain decomposition methods. *ffddm* already implements the GenEO method [40] for SPD problems, and its building blocks are designed to simplify the implementation and prototyping of new domain decomposition methods such as the saddle point solver presented in this paper. The *ffddm* documentation is available on the FreeFem.org web page, see [43].

Numerical results are obtained on the french GENCI supercomputer *Occigen*, on the Haswell partition composed of 50544 cores of Intel Xeon E5-2690V3 processors clocked at 2.6 GHz. The interconnect is an InfiniBand FDR 14 pruned fat tree. We use Intel compilers, the Intel Math Kernel Library and Intel MPI.

As is usually done in domain decomposition methods, we assign one subdomain per MPI process. Our implementation is pure MPI and no multithreading is done ; we assign one MPI process per computing core. The mesh of the computational domain is partitioned using the automatic graph partitioner *Metis* [30] (see Figure 1). Local subdomain matrices are factorized by the sparse direct solver *MUMPS* [1]. Local eigenvalue problems are solved with *Arpack* [32] ; both libraries are interfaced with *FreeFEM*. GenEO coarse space matrices $R_0AR_0^T$ in (2) and $Z_{S_1}^T S_1 Z_{S_1}$ in (21) are assembled and factorized in a distributed manner on a few cores (24 in most of the experiments) using the parallel solver *MUMPS*.

<figure><img src="image_3.png" /><figcaption>Figure 1: Heterogeneous beam composed of 10 alternating layers of rubber and steel. Coefficient distribution (left) and mesh partitioning into 16 subdomains by the automatic graph partitioner *Metis* (right).</figcaption></figure>

For illustration purposes, we represent in Figure 2 (top) the eigenvalues of the local GenEO eigenvalue problems for both coarse spaces, $V_0$ for A and $\tilde{W}_0$ for $S_1$ (the former corresponding to eq. (5) and the latter to eq. (19)), for the heterogeneous beam problem with 16 million degrees of freedom, corresponding to the first row of Table 1. The figures show the inverse of the first 40 largest eigenvalues for 10 of the subdomains for the experiment corresponding to the first row of Table 1, so that the eigenvectors corresponding to the smallest values on the graphs (below the dashed line) will be selected to enter the coarse space. For comparison, we also solve the constant coefficient problem corresponding to an homogeneous steel (and compressible) beam and show the eigenvalues in Figure 2 (bottom).

We can see the effect of heterogeneities on the spread of eigenvalues for different subdomains compared to the homogeneous case. In addition, we retrieve the 6 eigenvalues corresponding to the rigid body modes for A in the homogeneous case, and we see that we need a larger set of eigenvectors in order to build a robust coarse space in the heterogeneous case. Figure 2 also shows that there is no need for a coarse space for $S_1$ in the compressible homogeneous case. A strong feature of the GenEO method is that relevant
<figure><img src="image_4.png" /><figcaption>Figure 2: Top: Heterogeneous steel and rubber beam. Bottom: Homogeneous steel only beam. Inverse of the eigenvalues of the local GenEO eigenvalue problems for both coarse spaces, $V_0$ for A (left) and $\tilde{W}_0$ for $S_1$ (right), for 10 of the subdomains.</figcaption></figure>

eigenvectors to enter the coarse space are selected automatically, adapting to the problem at hand and its spatial heterogeneity. Moreover, the robustness of the preconditioner does not rely on a specific partitioning, which allows the use of automatic graph partitioners such as *Metis* [30] or *Scotch* [10].

## 5.2 Parameters of the method

The method has a few parameters in play:

* The number of layers of mesh elements in the overlap region between subdomains is 2 for the velocity blocks (corresponding to $R_i$ in (2)) and 4 for the pressure blocks (corresponding to $\tilde{R}_i$ in (10)). This corresponds to the minimum overlap that satisfies relation (9) for a symmetric construction of the overlap between subdomains.

* For the heterogeneous beam problem, we set the threshold $\tau_A$ for selecting the local eigenvectors entering the coarse space $V_0$ to $\tau_A = 10$ (corresponding to $1/\lambda = 0.1$ on Figure 2, left). The threshold $\tau_{S_1}$ for selecting the local eigenvectors entering the coarse space $\tilde{W}_0$ is set to $\tau_{S_1} = 3.33$ (corresponding to $1/\lambda = 0.3$ on Figure 2, right).

The selection of these thresholds is based on the fact that the iteration counts are very stable in both weak and strong scaling experiments. It means that for a given physics the parameters $\tau_A$ and $\tau_{S_1}$ can be tuned for a small test case and then used in large scale experiments.

## 5.3 Weak scalability test for heterogeneous steel and rubber beam

Here we present weak scalability results for the heterogeneous beam composed of 10 alternating layers of rubber ($E_1, \nu_1$) = (1 × 10⁷, 0.4999) and steel ($E_2, \nu_2$) = (2 × 10⁹, 0.35).
Local problem size is kept roughly constant as $N$ grows, and the total number of dofs $n$ goes from 16 million on 262 cores to 1 billion on 16800 cores.

We report in Table 1 the iteration counts and computing times for the DD saddle point solver Algorithm 2. Note that in Algorithm 2 we replace PCG by right-preconditioned GMRES for step 3. The stopping criterion is a tolerance smaller than $10^{-5}$. Moreover, we use flexible GMRES, as we solve (25) inexacty using GMRES with a tolerance of $10^{-2}$ in order to apply $N_S^{-1}$. In order, columns correspond to: number of cores, number of dofs $n$, size of the coarse space for A $dim(V_0)$, size of the coarse space for $S_1 dim(W_0)$, setup time corresponding to the assembly and factorization of the various local and coarse operators, number of outer GMRES iterations, GMRES computing time, total computing time (setup + GMRES) and average number of inner GMRES iterations for each solution of (25). All timings are reported in seconds.

<table><caption>Table 1: Weak scaling experiment for 3D heterogeneous elasticity: beam with 10 alternating layers of steel and rubber. Reported iteration counts and timings for DD saddle point algorithm 2.</caption><thead><tr><th>#cores</th><th>n</th><th>dim(V<sub>0</sub>)</th><th>dim(W&#x0303;<sub>0</sub>)</th><th>setup(s)</th><th>#It</th><th>gmres(s)</th><th>total(s)</th><th>#It N<sub>S</sub><sup>-1</sup></th></tr></thead><tbody><tr><td>262</td><td>15 987 380</td><td>5 383</td><td>3 319</td><td>710.7</td><td>24</td><td>631.6</td><td>1342.3</td><td>11</td></tr><tr><td>525</td><td>27 545 495</td><td>9 959</td><td>2 669</td><td>526.6</td><td>21</td><td>519.5</td><td>1046.1</td><td>12</td></tr><tr><td>1 050</td><td>64 982 431</td><td>17 837</td><td>4 587</td><td>675.2</td><td>22</td><td>665.9</td><td>1341.1</td><td>11</td></tr><tr><td>2 100</td><td>126 569 042</td><td>32 361</td><td>7 995</td><td>689.2</td><td>25</td><td>733.8</td><td>1423.0</td><td>10</td></tr><tr><td>4 200</td><td>218 337 384</td><td>59 704</td><td>13 912</td><td>593.0</td><td>27</td><td>705.4</td><td>1298.4</td><td>10</td></tr><tr><td>8 400</td><td>515 921 881</td><td>141 421</td><td>25 949</td><td>735.8</td><td>32</td><td>1152.5</td><td>1888.3</td><td>10</td></tr><tr><td>16 800</td><td>1 006 250 208</td><td>260 348</td><td>41 341</td><td>819.2</td><td>29</td><td>1717.9</td><td>2537.1</td><td>12</td></tr></tbody></table>

**Iteration counts.** We first discuss iteration counts. We see that outer iteration count remains stable, between 21 and 32. The inner iteration count is also stable and remains around 11. We also observed (figures are not reported here) than the inner GMRES tolerance of $10^{-2}$ does not affect the outer iteration count compared to an accurate solution with a stricter tolerance of $10^{-5}$, and allows a significant reduction in inner iteration count. For example, 11 iterations on average instead of 28 on 1050 cores for the same outer iteration count of 22, leading to a decrease from 1178.2 to 665.9 seconds in GMRES timing.

**Timings.** In terms of setup timings, the computing time remains relatively stable, with roughly 15% increase for a factor of 64 in problem size. Around 60% of the setup time is spent in the solution of the eigenvalue problems (19) for $S_1$.

The solution time stays relatively stable up to 4200 cores, where it starts to degrade. This can be related to the increased cost of the coarse space solves with matrices $R_0AR_0^T$ and $Z_{S_1}^T S_1 Z_{S_1}$ as their size increases: total time spent in coarse space solves is 14.7, 62.1 and 679.3 seconds on 262, 4200 and 16800 cores respectively. A possible improvement would be to use a multi-level method to solve the coarse problems iteratively.

## 5.4 Strong scalability test for heterogeneous steel and rubber beam

Strong scalability results for the same test case as above consisting of an heterogeneous beam composed of 10 alternating layers of rubber and steel are presented in Table 2. The problem size is 27.5 million and the strong scaling test ranges from 525 to 4200 cores. Iteration counts and computing times for the DD saddle point solver Algorithm 2 are reported.

**Iteration counts.** Outer iteration count remains stable, with a slight increase from 21 to 23. Inner iteration counts remains also stable, even slightly decreasing from 12 to 9.
<table><caption>Table 2: Strong scaling experiment for 3D heterogeneous elasticity: beam with 10 alternating layers of steel and rubber. Reported iteration counts and timings for DD saddle point Algorithm 2.</caption><thead><tr><th>#cores</th><th>n</th><th>dim(V<sub>0</sub>)</th><th>dim(W&#x0303;<sub>0</sub>)</th><th>setup(s)</th><th>#It</th><th>gmres(s)</th><th>total(s)</th><th>#It N<sub>S</sub><sup>-1</sup></th></tr></thead><tbody><tr><td>525</td><td>27 545 495</td><td>9 959</td><td>2 669</td><td>526.6</td><td>21</td><td>519.5</td><td>1046.1</td><td>12</td></tr><tr><td>1 050</td><td>27 545 495</td><td>15 078</td><td>4 082</td><td>265.7</td><td>21</td><td>224.7</td><td>490.4</td><td>11</td></tr><tr><td>2 100</td><td>27 545 495</td><td>23 172</td><td>6 453</td><td>168.8</td><td>23</td><td>131.1</td><td>299.9</td><td>10</td></tr><tr><td>4 200</td><td>27 545 495</td><td>37 768</td><td>11 152</td><td>103.8</td><td>23</td><td>91.3</td><td>195.1</td><td>9</td></tr></tbody></table>

**Timings.** We see that the setup timing decreases accordingly as the subdomains shrink in size, from 526.6 seconds on 525 cores to 103.8 seconds on 4200 cores; the speedup efficiency with respect to 525 cores ranges from 99% on 1050 cores to 63% on 4200 cores.

We see a similar trend for the solution time, ranging from 519.5 seconds on 525 cores to 91.3 seconds on 4200 cores; the speedup efficiency with respect to 525 cores ranges from 116% on 1050 cores to 71% on 4200 cores. This decrease in efficiency can be explained by the increased relative cost of coarse space solves as subdomains get smaller: from 3% on 525 cores to 30% on 4200 cores. The added overlap also plays a greater role in the loss of efficiency as subdomains get smaller.

## 5.5 Comparison with other algorithms

### 5.5.1 Comparison with a direct solver

<table><caption>Table 3: Comparison with the parallel sparse direct solver *MUMPS* for 3D heterogeneous elasticity: beam with 10 alternating layers of steel and rubber. Reported timings for four discretization levels while also varying the number of cores (OOM means the computation ran out of available memory).</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">#cores</th><th colspan="3">MUMPS</th><th colspan="4">DD saddle point solver</th></tr><tr><th>setup(s)</th><th>solve(s)</th><th>total(s)</th><th>setup(s)</th><th>#It</th><th>gmres(s)</th><th>total(s)</th></tr></thead><tbody><tr><td>139 809</td><td>16</td><td>7.1</td><td>0.1</td><td>7.2</td><td>27.1</td><td>18</td><td>19.7</td><td>46.8</td></tr><tr><td>1 058 312</td><td>32</td><td>85.7</td><td>0.8</td><td>86.5</td><td>166.2</td><td>20</td><td>137.2</td><td>303.4</td></tr><tr><td>1 058 312</td><td>65</td><td>71.0</td><td>0.6</td><td>71.6</td><td>91.0</td><td>21</td><td>77.1</td><td>168.1</td></tr><tr><td>1 058 312</td><td>131</td><td>63.2</td><td>0.5</td><td>63.7</td><td>59.7</td><td>24</td><td>49.7</td><td>109.4</td></tr><tr><td>3 505 582</td><td>55</td><td>477.8</td><td>3.7</td><td>481.5</td><td>404.1</td><td>24</td><td>430.1</td><td>834.2</td></tr><tr><td>3 505 582</td><td>110</td><td>392.3</td><td>2.3</td><td>394.6</td><td>242.5</td><td>23</td><td>212.8</td><td>455.3</td></tr><tr><td>3 505 582</td><td>221</td><td>387.0</td><td>2.1</td><td>389.1</td><td>134.8</td><td>23</td><td>109.4</td><td>244.2</td></tr><tr><td>3 505 582</td><td>442</td><td>453.9</td><td>2.2</td><td>456.1</td><td>88.2</td><td>24</td><td>68.6</td><td>156.8</td></tr><tr><td>8 235 197</td><td>262</td><td>OOM</td><td>/</td><td>/</td><td>278.5</td><td>25</td><td>264.3</td><td>542.8</td></tr><tr><td>8 235 197</td><td>525</td><td>1622.1</td><td>6.1</td><td>1628.2</td><td>172.1</td><td>24</td><td>136.0</td><td>308.1</td></tr><tr><td>8 235 197</td><td>1050</td><td>1994.3</td><td>7.4</td><td>2001.7</td><td>136.5</td><td>25</td><td>99.7</td><td>236.2</td></tr></tbody></table>

In Table 3, we compare the performance of the solver to the parallel sparse direct solver *MUMPS* for the heterogeneous steel and rubber beam test case with four discretization levels, while also varying the number of cores. As we can see, *MUMPS* is comparatively more efficient for smaller problems, with for example a total time of 86.5 seconds compared to 303.4 seconds for our saddle point solver for 1 million unknowns on 32 cores. However, as expected, we see a large increase in memory and computational cost as the size of the system gets larger: for 8.2 million unknowns, *MUMPS* runs out of memory on 262 cores and solves the problem in 1628.2 seconds on 525 cores, compared to 308.1 seconds for the
DD solver. Moreover, we can see from Table 3 that the DD saddle point solver offers much better strong scalability.

### 5.5.2 Comparison with other Schwarz methods

In Figure 3, we plot the convergence history of GMRES for both compressible homogeneous and heterogeneous steel-rubber cases with 27.5 million unknowns on 525 cores (left), with a comparison to the standard one-level Additive Schwarz Method (ASM) from PETSc [4] on the global problem (right), illustrating the difficulty of the test case at hand. Our saddle point solver needs 15 and 21 iterations to converge for the homogeneous and heterogeneous cases respectively, compared to 856 and 2880 for the one-level method, with significant undesired plateaux.

Note also that for the same physical test case, iteration counts are better than in [23] but timings are not as good. As mentioned in the introduction, our method has the advantage of a provable convergence estimate and to not depend on the design of specific absorbing conditions for the elasticity system.

<figure><img src="image_5.png" /><figcaption>Figure 3: GMRES convergence history of the saddle point solver (left) compared to the one-level Additive Schwarz Method (right) for the homogeneous steel beam and heterogeneous rubber/steel beam problems discretized with 27.5 million unknowns (corresponding to the first row of Table 2), on 525 cores.</figcaption></figure>

### 5.5.3 Comparison with algebraic multigrid method

We also performed comparisons with the Geometric Algebraic Multigrid (GAMG) preconditioner from PETSc. We were not able to find a suitable tuning of parameters for GAMG for the saddle point formulation. However, we performed comparisons between GAMG and standard GenEO for the velocity formulation on the homogeneous beam, varying the Poisson ratio $\nu$ from 0.48 to 0.499. The near kernel of the problem which corresponds to rigid body modes is provided to the multigrid method. The GenEO threshold $\tau$ is set to 3.33, and we select at most 80 eigenvectors in each subdomain. Even though GAMG is faster for $\nu \le 0.49$, we can see that GenEO is more robust as $\nu$ increases. In particular, GAMG fails to converge in 2000 iterations for $\nu \ge 0.495$.
<table><caption>Table 4: GAMG versus standard GenEO for the velocity formulation on the homogeneous beam discretized with 7.9 million unknowns, using 131 and 525 cores. Reported iteration counts and timings for different values of the Poisson ratio $\nu$ ranging from 0.48 to 0.499.</caption><thead><tr><th colspan="2">131 cores</th><th colspan="2">GAMG</th><th colspan="4">GenEO</th></tr><tr><th>ν</th><th>#It</th><th>total(s)</th><th>dim(V<sub>0</sub>)</th><th>setup(s)</th><th>#It</th><th>gmres(s)</th><th>total(s)</th></tr></thead><tbody><tr><td>0.48</td><td>60</td><td>67.1</td><td>10 480</td><td>200.7</td><td>24</td><td>11.2</td><td>212.0</td></tr><tr><td>0.485</td><td>109</td><td>89.0</td><td>10 480</td><td>199.5</td><td>27</td><td>12.7</td><td>212.2</td></tr><tr><td>0.49</td><td>210</td><td>137.0</td><td>10 480</td><td>202.0</td><td>32</td><td>15.0</td><td>217.0</td></tr><tr><td>0.495</td><td>>2000</td><td>/</td><td>10 480</td><td>199.9</td><td>43</td><td>20.2</td><td>220.1</td></tr><tr><td>0.499</td><td>>2000</td><td>/</td><td>10 480</td><td>199.2</td><td>99</td><td>48.6</td><td>247.7</td></tr></tbody></table><table><thead><tr><th colspan="2">525 cores</th><th colspan="2">GAMG</th><th colspan="4">GenEO</th></tr><tr><th>ν</th><th>#It</th><th>total(s)</th><th>dim(V<sub>0</sub>)</th><th>setup(s)</th><th>#It</th><th>gmres(s)</th><th>total(s)</th></tr></thead><tbody><tr><td>0.48</td><td>56</td><td>25.5</td><td>41 766</td><td>60.4</td><td>18</td><td>5.0</td><td>65.4</td></tr><tr><td>0.485</td><td>60</td><td>26.1</td><td>41 984</td><td>60.9</td><td>20</td><td>5.3</td><td>66.2</td></tr><tr><td>0.49</td><td>116</td><td>33.3</td><td>42 000</td><td>60.4</td><td>23</td><td>5.9</td><td>66.3</td></tr><tr><td>0.495</td><td>>2000</td><td>/</td><td>42 000</td><td>60.4</td><td>32</td><td>7.6</td><td>68.1</td></tr><tr><td>0.499</td><td>>2000</td><td>/</td><td>42 000</td><td>60.6</td><td>95</td><td>20.3</td><td>81.0</td></tr></tbody></table>

# 6 Conclusion and outlook

Under the assumption that the diagonal block matrices of a saddle point problem are spectrally equivalent to a sum of positive semi definite matrices, we have introduced an adaptive domain decomposition (DD) method. The key ingredients of the algorithm are a two level Schwarz method for the primal unknowns block and a two-level Neumann-Neuman for the Schur complement problem on the dual unknowns.

Two coarse spaces are built by solving generalized eigenvalue problems, one for the primal unknowns and the second one for the dual unknowns. The robustness of the method was assessed on a notoriously difficult three dimensional elasticity problem for a steel-rubber structure discretized with continuous pressure.

Several issues deserve further investigations. First a multilevel method with more than two levels would enable even larger and possibly faster simulations. Also, the tests were performed with *FreeFEM* scripts using the standalone *ffddm* [43] framework written in the *FreeFEM* language. The integration of the method in the C++/MPI library *hpddm* [28], which already implements the GenEO methods [40, 22], could lead to faster codes and a more general diffusion of the saddle point preconditioner since *hpddm* is callable [29] from the PETSc library [3]. In a different setting, the design of adaptive coarse space is strongly connected to multiscale finite element (MFE) methods (see e.g., [14, 27, 33] and references therein) and this work could be used in designing MFE methods for saddle point problems.

# Acknowledgment

This work was granted access to the HPC resources of OCCIGEN@CINES under the allocation 2021-067730 granted by GENCI.

# Code reproducibility

Our numerical results can be reproduced running the script https://github.com/FreeFem/FreeFem-sources/blob/master/examples/ffddm/elasticity_saddlepoint.edp available in the FreeFem distribution starting from version 4.10.
References

[1] Patrick R. Amestoy, Iain S. Duff, Jean-Yves L'Excellent, and Jacko Koster. A fully asynchronous multifrontal solver using distributed dynamic scheduling. *SIAM J. Matrix Analysis and Applications*, 23(1):15-41, 2001.

[2] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. Multigrid in h (div) and h (curl). *Numerische Mathematik*, 85(2):197-217, Apr 2000.

[3] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, and Hong Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.5, Argonne National Laboratory, 2014.

[4] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, *Modern Software Tools in Scientific Computing*, pages 163-202. Birkhäuser Press, 1997.

[5] Michele Benzi, Gene H. Golub, and Jörg Liesen. Numerical solution of saddle point problems. *Acta Numer.*, 14:1-137, 2005.

[6] Michele Benzi and Andrew J. Wathen. Some preconditioning techniques for saddle point problems. In *Model order reduction: theory, research aspects and applications*, volume 13 of *Math. Ind.*, pages 195-211. Springer, Berlin, 2008.

[7] Jean-François Bourgat, Roland Glowinski, Patrick Le Tallec, and Marina Vidrascu. Variational formulation and algorithm for trace operator in domain decomposition calculations. In Tony Chan, Roland Glowinski, Jacques Périaux, and Olof Widlund, editors, *Domain Decomposition Methods*, pages 3-16, Philadelphia, PA, 1989. SIAM.

[8] Franco Brezzi and Michel Fortin. *Mixed and hybrid finite element methods*, volume 15. Springer Science & Business Media, 2012.

[9] J Cahouet and J-P Chabard. Some fast 3d finite element solvers for the generalized stokes problem. *International Journal for Numerical Methods in Fluids*, 8(8):869-895, 1988.

[10] C. Chevalier and F. Pellegrini. PT-SCOTCH: a tool for efficient parallel graph ordering. *Parallel Computing*, 6-8(34):318-331, 2008.

[11] Eric de Sturler and Jörg Liesen. Block-diagonal and constraint preconditioners for nonsymmetric indefinite linear systems. I. Theory. *SIAM J. Sci. Comput.*, 26(5):1598-1619, 2005.

[12] Victorita Dolean, Pierre Jolivet, and Frédéric Nataf. *An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation*. SIAM, 2015.

[13] Daniel Drzisga, Lorenz John, Ulrich Rude, Barbara Wohlmuth, and Walter Zulehner. On the analysis of block smoothers for saddle point problems. *SIAM Journal on Matrix Analysis and Applications*, 39(2):932-960, 2018.

[14] Yalchin Efendiev, Juan Galvis, and Thomas Y Hou. Generalized multiscale finite element methods (gmsfem). *Journal of Computational Physics*, 251:116-135, 2013.

[15] Charbel Farhat and Francois-Xavier Roux. A method of Finite Element Tearing and Interconnecting and its parallel solution algorithm. *Int. J. Numer. Meth. Engrg.*, 32:1205-1227, 1991.
[16] Patrick E. Farrell, Lawrence Mitchell, and Florian Wechsung. An Augmented Lagrangian Preconditioner for the 3D Stationary Incompressible Navier-Stokes Equations at High Reynolds Number. *SIAM J. Sci. Comput.*, 41(5):A3073-A3096, 2019.

[17] Michel Fortin and Roland Glowinski. *Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems*. Elsevier, 2000.

[18] R Glowinski and P Le Tallec. Numerical solution of problems in incompressible finite elasticity by augmented lagrangian methods. i. two-dimensional and axisymmetric problems. *SIAM Journal on Applied Mathematics*, 42(2):400-429, 1982.

[19] R. Glowinski and P. Le Tallec. Numerical solution of problems in incompressible finite elasticity by augmented lagrangian methods ii. three-dimensional problems. *SIAM Journal on Applied Mathematics*, 44(4):710-733, 1984.

[20] Roland Glowinski and Patrick Le Tallec. *Augmented Lagrangian and operator-splitting methods in nonlinear mechanics*. SIAM, 1989.

[21] M. Griebel and P. Oswald. On the abstract theory of additive and multiplicative Schwarz algorithms. *Numer. Math.*, 70(2):163-180, 1995.

[22] R. Haferssas, P. Jolivet, and F. Nataf. An Additive Schwarz Method Type Theory for Lions's Algorithm and a Symmetrized Optimized Restricted Additive Schwarz Method. *SIAM J. Sci. Comput.*, 39(4):A1345-A1365, 2017.

[23] Ryadh Haferssas, Pierre Jolivet, and Frédéric Nataf. An additive schwarz method type theory for lions's algorithm and a symmetrized optimized restricted additive schwarz method. *SIAM Journal on Scientific Computing*, 39(4):A1345-A1365, 2017.

[24] F. Hecht. New development in Freefem++. *J. Numer. Math.*, 20(3-4):251-265, 2012.

[25] R. Hiptmair. Multigrid method for Maxwell's equations. *SIAM J. Numer. Anal.*, 36(1):204-225, 1998.

[26] Ralf Hiptmair. Multigrid method for h (div) in three dimensions. *Electron. Trans. Numer. Anal*, 6(1):133-152, 1997.

[27] Patrick Jenny, Seong H Lee, and Hamdi A Tchelepi. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. *Multiscale Modeling & Simulation*, 3(1):50-64, 2005.

[28] Pierre Jolivet and Frédéric Nataf. HPDDM: High-Performance Unified framework for Domain Decomposition methods, MPI-C++ library. https://github.com/hpddm/hpddm, 2014.

[29] Pierre Jolivet, Jose E. Roman, and Stefano Zampini. KSPHPDDM and PCHPDDM: Extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. *Computers & Mathematics with Applications*, 84:277-295, 2021.

[30] G. Karypis and V. Kumar. METIS: A software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices. Technical report, Department of Computer Science, University of Minnesota, 1998. http://glaros.dtc.umn.edu/gkhome/views/metis.

[31] Axel Klawonn. An optimal preconditioner for a class of saddle point problems with a penalty term. *SIAM Journal on Scientific Computing*, 19(2):540-552, 1998.
[32] Richard B Lehoucq, Danny C Sorensen, and Chao Yang. ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods, volume 6. SIAM, 1998.

[33] Chupeng Ma, Robert Scheichl, and Tim Dodwell. Novel design and analysis of generalized FE methods based on locally optimal spectral approximations. *arXiv preprint arXiv:2103.09545*, 2021.

[34] Malcolm F. Murphy, Gene H. Golub, and Andrew J. Wathen. A note on preconditioning for indefinite linear systems. *SIAM J. Sci. Comput.*, 21(6):1969–1972, 2000.

[35] Sergey V. Nepomnyaschikh. Mesh theorems of traces, normalizations of function traces and their inversions. *Sov. J. Numer. Anal. Math. Modeling*, 6:1–25, 1991.

[36] Joseph E Pasciak and Jun Zhao. Overlapping schwarz methods in h (curl) on polyhedral domains. *Journal of Numerical Mathematics*, 10(3):221–234, 2002.

[37] Luca F. Pavarino and Olof B. Widlund. Balancing Neumann-Neumann methods for incompressible Stokes equations. *Comm. Pure Appl. Math.*, 55(3):302–335, 2002.

[38] T. Rees and M. Wathen. An element-based preconditioner for mixed finite element problems. *SIAM Journal on Scientific Computing*, 2020.

[39] S. Reitzinger and J. Schöberl. An algebraic multigrid method for finite element discretizations with edge elements. *Numerical Linear Algebra with Applications*, 9(3):223–238, 2002.

[40] Nicole Spillane, Victorita Dolean, Patrice Hauret, Frédéric Nataf, Clemens Pechstein, and Robert Scheichl. Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps. *Numer. Math.*, 126(4):741–770, 2014.

[41] Nicole Spillane and Daniel Rixen. Automatic spectral coarse spaces for robust finite element tearing and interconnecting and balanced domain decomposition algorithms. *Internat. J. Numer. Methods Engrg.*, 95(11):953–990, 2013.

[42] Andrea Toselli and Olof Widlund. *Domain Decomposition Methods - Algorithms and Theory*, volume 34 of Springer Series in Computational Mathematics. Springer, 2005.

[43] Pierre-Henri Tournier and Frédéric Nataf. FFDDM: Freefem domain decomposition method. https://doc.freefem.org/documentation/ffddm/index.html, 2019.

[44] Xuemin Tu and Jing Li. A FETI-DP type domain decomposition algorithm for three-dimensional incompressible Stokes equations. *SIAM J. Numer. Anal.*, 53(2):720–742, 2015.

[45] O Widlund, Stefano Zampini, S Scacchi, and L Pavarino. Block feti-dp/bddc preconditioners for mixed isogeometric discretizations of three-dimensional almost incompressible elasticity. *Mathematics of Computation*, 90(330):1773–1797, 2021.