Multilevel Spectral Domain Decomposition

06/11/2021 ∙ by Peter Bastian, et al. ∙ University of Heidelberg 0

Highly heterogeneous, anisotropic coefficients, e.g. in the simulation of carbon-fibre composite components, can lead to extremely challenging finite element systems. Direct solvers for the resulting large and sparse linear systems suffer from severe memory requirements and limited parallel scalability, while iterative solvers in general lack robustness. Two-level spectral domain decomposition methods can provide such robustness for symmetric positive definite linear systems, by using coarse spaces based on independent generalized eigenproblems in the subdomains. Rigorous condition number bounds are independent of mesh size, number of subdomains, as well as coefficient contrast. However, their parallel scalability is still limited by the fact that (in order to guarantee robustness) the coarse problem is solved via a direct method. In this paper, we introduce a multilevel variant in the context of subspace correction methods and provide a general convergence theory for its robust convergence for abstract, elliptic variational problems. Assumptions of the theory are verified for conforming, as well as for discontinuous Galerkin methods applied to a scalar diffusion problem. Numerical results illustrate the performance of the method for two- and three-dimensional problems and for various discretization schemes, in the context of scalar diffusion and linear elasticity.



There are no comments yet.


page 19

page 22

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this paper we are concerned with the solution of large and sparse linear systems



is a symmetric and positive definite (SPD) matrix arising from the discretization of an elliptic (system of) partial differential equation(s) (PDE) on a bounded domain

, with the spatial dimension.

Direct methods for solving (1) are very effective for relatively small problems but suffer from severe memory requirements and limited parallel scalability [21, 26]. We focus on iterative methods instead. In PDE applications, the mesh parameter needs to be chosen sufficiently small to control the error in the numerical solution, leading to very large systems . A variety of methods have been developed in the past that converge (almost) independently of the mesh size as well as of the number of processors – assuming a parallel partitioning of into subdomains with diameter bounded by . These include in particular multigrid (MG) and domain decomposition (DD) methods [28, 30, 11]. These methods require a “coarse solver component” providing global information transfer. In this paper we focus on extensions of the two-level overlapping Schwarz DD method.

While robustness with respect to (w.r.t.) and is achieved by many methods, robustness w.r.t. mesh anisotropy or to large variations in problem parameters, such as the permeability coefficient in porous media flow or the Lamé parameters in linear elasticity, are more difficult to achieve. Construction of robust coarse spaces was a theme in multigrid early on [4] and lead to the development of algebraic multigrid (AMG) methods [33]. Rigorous convergence bounds for aggregation-based AMG applied to nonsingular symmetric M-matrices with nonnegative row sums are provided in [23]. In the context of overlapping DD methods, it was shown in [27] based on weighted Poincaré inequalities [25] that the standard method can be robust w.r.t. strong coefficient variation inside subdomains, but this puts hard constraints on the domain decomposition. Construction of coarse spaces based on multiscale basis functions [1, 17]

can be very effective, but also leads to no rigorous robustness w.r.t. arbitrary coefficient variations. A breakthrough was achieved by constructing coarse spaces based on solving certain local generalized eigenvalue problems (GEVP). In

[24], a local eigenvalue problem involving the Dirichlet-to-Neumann map was introduced, and later analysed in [12] based on [25]. Different GEVP were proposed in [16, 13, 29]

together with an analysis that is solely based on approximation properties of the chosen local eigenfunctions.

In this paper, we consider extensions of the GenEO (Generalized Eigenproblems in the Overlap) method introduced in [29]. The method works on general elliptic systems of PDEs and is rigorously shown to be robust w.r.t. coefficient variations when all eigenfunctions w.r.t. eigenvalues below a threshold are included into the coarse space. This number can be related to the number of isolated high-conductivity regions, [16], but also depends on the specific GEVP used. The local GEVP can either be constructed from local stiffness matrices or in a fully algebraic way from the global stiffness matrix through a symmetric positive semi-definite (SPSD) splitting [2]. An extension of the GenEO approach to nonoverlapping domain decomposition methods as well as non-selfadjoint problems can be found in [18].

Two-level domain decomposition methods traditionally employ direct solvers in the subdomains and on the coarse level. While iterative solvers could be employed, they then need to be robust with respect to coefficient variations as well. If such a solver would be available, it could be used instead of the domain decomposition method. Thus we assume that such solver is not available. The use of direct solvers in the subdomains and on the coarse level puts an upper limit on the size of the subdomain/coarse problem and thus on the total problem size due to run-time and memory requirements. The more severe penalty is typically imposed by the memory requirements, in particular for todays supercomputers which have many cores with relatively little memory per core. In order to give a concrete example, consider the HAWK system of HLRS, Stuttgart, Germany. It consists of 5632 nodes, each containing 2 CPUs with 64 cores each and 256 GB of memory, i.e. 2GB of memory per core. This limits the fine level subdomain size to about degrees of freedom (dof) in scalar, three-dimensional problems. In order to maintain a good speedup (in the setup phase), the coarse system size is then also limited to 250000 dof or about 12500 subdomains (or cores) when we assume 20 basis vectors per subdomain from the GEVP. Thus the GenEO method would not scale to the full Hawk machine. This estimate could be improved by employing parallel direct solvers, but the scalability of these solvers is limited [21, 26, 3].

Scalability beyond cores can be achieved by employing more than two levels. A robust multilevel method employing spectral coarse spaces based on the work [16] has been proposed in [31]. It employs a hierarchy of finite element meshes, GEVPs that are very similar to the ones that we propose in this paper and a nonlinear AMLI-cycle. Robustness and level independence is achieved with a -cycle, which is not desirable in a parallel method. While robustness with respect to coefficient variations is proven and demonstrated numerically, the problem sizes are rather small. A multilevel version of GenEO was proposed in [3] based on the SPSD splitting introduced in [2].

In this paper, we formulate a natural extension of GenEO [29] to an arbitrary number of levels. In contrast to [3], our method is formulated in a variational framework based on subspace correction [32]. The convergence theory of [29] is generalized to nonconforming discretizations as well as to multiple levels and several different GEVPs. This enables us to prove robust convergence also for discontinuous Galerkin methods suitable for high coefficient contrast [15]. Condition number bounds are derived for an iterated two-level method as in [3] but also for a fully additive multilevel preconditioner. Numerical results up to subdomains and more than dofs demonstrate the effectiveness of the approach.

The paper is organized as follows. In Section 2 we formulate the multilevel spectral domain decomposition method in a variational framework. In Section 3 we provide the analysis of the method and briefly describe our implementation in Section 4. Numerical results follow in Section 5 and we end with conclusions in Section 6.

2 Multilevel Spectral Domain Decomposition

2.1 Subspace Correction

Throughout the paper we assume the linear system (1) arises by inserting a basis representation into the variational problem


where is a finite element space (a finite dimensional vector space), is a symmetric positive definite bilinear form and is a linear form. The subscript denotes that and are defined on a shape regular mesh consisting of elements with diameter at most , partitioning the domain . Elements are assumed to be open and the image of a reference element under the diffeomorphism . Reference elements are either the reference simplex or the reference cube in dimension .

Subspace correction methods [32] are based on a splitting

of into possibly overlapping subspaces . Any such splitting gives rise to the iterative parallel subspace correction method


Here, is a suitably chosen damping factor. Sequential subspace correction, see [32], typically converges faster but offers less parallelism. Parallel subspace correction is also called additive subspace correction and sequential subspace correction is called multiplicative subspace correction due to the form of the error propagation operator. Hybrid variants may offer a good compromise in practice.

Practical implementation of subspace correction employs a basis representation

The rectangular matrices represent the basis of in terms of the basis of . Expanding leads to the algebraic form

with , and the preconditioner . is typically used as a preconditioner in the conjugate gradient method.

Multilevel spectral domain decomposition methods are introduced below in the framework of subspace correction. They employ a decomposition


where is the number of levels with being the finest level and the coarsest level. Identifying we have the nested level-wise spaces for . On each level , is split into subspaces .

2.2 Hierarchical Domain Decomposition

The construction of the subspaces is related to a hierarchical decomposition of the domain into subdomains for and . “Hierarchical” means that each subdomain is the union of subdomains on the next finer level . In particular, our multilevel method employs a single fine mesh given by the user. All subdomains are unions of elements of the mesh . Figure 1 shows a decomposition of a triangular mesh employing three levels.

The domain decomposition is obtained from a decomposition of as follows:

  1. On the finest level , decompose into overlapping sets

    by first partitioning into nonoverlapping sets of elements using a graph partitioner such es ParMetis [19] and then adding a user defined overlap in terms of layers of elements.

  2. On levels , determine decompositions of the subdomain index sets

    The sets may overlap but are not required to. Such decomposition is obtained by a graph partitioner using the subdomain graph instead of the mesh. With this we obtain the mesh partitioning on level as

  3. The subdomains are now defined from the mesh decomposition as

  4. On the coarsest level 0 we employ always only one subdomain.

Figure 1: Multilevel domain decomposition in two dimensions. A triangular mesh with 1.6 million vertices is partitioned into 4096 (level 3), 256 (level 2) and 16 (level 1) nested subdomains using ParMetis [19]. No overlap is used on the coarse levels.

For finite volume and discontinuous Galerkin methods we need to introduce notation for mesh faces. is an interior face if it is the intersection of two elements and has dimension . All interior faces are collected in the set . Likewise, is a boundary face if it is the intersection of some element with of dimension . All boundary faces make up the set . With each we associate a unit normal vector oriented from to . For a boundary face its unit normal coincides with the unit normal to . Related to the submeshes we define the corresponding sets of faces equationparentequation


2.3 Spectral Coarse Space Construction

Based on the hierarchical domain decomposition we can now formulate the construction of the coarse spaces and introduced in (4). First define the auxiliary local spaces


The restriction operator , , restricts the domain of a finite element function. For being zero on , the extension operator extends functions by zero outside of . Another major ingredient is the partition of unity.

Definition 1

A partition of unity on level is a family of operators such that

  1. is zero on for all and

  2. for all ,

Restriction and extension operators as well as the partition of unity operators with the required properties can be defined for conforming finite element spaces as well as discontinuous Galerkin finite element spaces.

The construction of the coarse spaces is now recursive over the levels:

  1. Set . Set .

  2. For each subdomain solve a GEVP


    with appropriately defined local bilinear forms and

    detailed below. Let eigenvalues and corresponding eigenvectors be ordered s.t.


  3. With a user defined parameter define the coarse space as

  4. Set . If goto step 2, otherwise stop.

The subspaces for the subspace correction method based on are then given by


together with .

Remark 1

The most important properties of this construction are:

  • Within each level all GEVPs can be solved in parallel.

  • Basis functions have support only in .

  • Functions in are zero at their subdomain boundary (except on the global Neumann boundary). Functions in are not necessarily zero at their subdomain boundary (except on the global Dirichlet boundary).

3 Convergence Theory

3.1 Standard Subspace Correction Theory

For the standard analysis of subspace correction methods, the following -orthogonal projections and are introduced [32, 30]:

This allows to write the error propagation operator of the iteration (3) as

Taking , the convergence factor of the iteration (3) is with the spectral condition number . The goal of the analysis below is to provide upper and lower bounds of the form

which in turn give a bound on the condition number .

The analysis is based on two major properties.

Definition 2 (Coloring and domain decomposition)
  1. We say the multilevel domain decomposition admits a level-wise coloring with colors, if for each level there exists a map such that

  2. The multilevel domain decomposition is called admissible, if on each level every interior face is interior to at least one subdomain.

Lemma 1

Let the multilevel domain decomposition have a finite coloring with colors. Then parallel subspace correction satisfies the upper bound


See [28, p. 182].

Definition 3 (Stable splitting)

The subspaces and , , , admit a stable splitting if there exists and for each a decomposition , , , such that

Lemma 2

If the subspaces and , , admit a stable splitting with constant , parallel subspace correction satisfies the lower bound


See e.g. [30].

Theorem 1

If the subspaces and , , have a finite coloring with colors and admit a stable splitting with constant , parallel subspace correction satisfies the bound


Follows immediately from Lemma 1 and Lemma 2.

3.2 Abstract Schwarz Theory for Spectral Domain Decomposition

We now generalize the theory in [29] to the multilevel case and to more general discretization schemes. For each application, the following three definitions have to be verified.

Definition 4 (Strengthened triangle inequality under the square)

The domain decomposition allows a strengthened triangle inequality under the square if there exists independent of and such that for any collection of :

Here is the norm induced by .

Definition 5 (Positive semi-definite splitting)

The bilinear forms introduced in (7) provide a positive semi-definite splitting of if there exists independent of and such that for each level :

Here is the semi-norm induced by . Symmetric positive semi-definite splittings on the algebraic level were introduced in [2, 3].

Definition 6 (Local stability)

The subspaces and are called locally stable if there exists a constant idependent of and for each , , a decomposition with and such that the following inequalities hold:

Lemma 3 (Levelwise stability)

Let the spectral coarse spaces admit strengthened triangle inequalities under the square, let the bilinear forms provide a symmetric positive semi-definite splitting and let the subspaces be locally stable. Then, for each level there exists a decomposition with and such that



Adaptation of [29, Lemma 2.9]. From Def. 6 and Def. 5 we get

Thus, it follows from and from Def. 4 that

Combining both results yields

The two-level decomposition from Definition 6 implies a multilevel decomposition as follows: For any given function :


We can now formulate the first major result of our paper.

Lemma 4 (Multilevel stability)

Let the spectral coarse spaces admit strengthened triangle inequalities under the square, let the bilinear forms provide a symmetric positive semi-definite splitting and let the subspaces be locally stable. Then the multilevel decomposition (10) satisfies for any

with .


For any level , we obtain as in the proof of Lemma 3

Furthermore, Lemma 3 implies

for . Together we obtain with setting

Theorem 2

Let the assumptions of Lemma 4 hold and let the domain decomposition admit a coloring with colors on each level. Then the parallel multilevel subspace correction method satisfies the condition number bound

with .


Follows from Lemma 4 and Theorem 1.

Remark 2

The upper bound in Lemma 3 predicts an exponential increase of the condition number with the number of levels . Since we are only interested in a moderate number of levels or , this may be acceptable. Moreover, our numerical results below suggest that the bound is pessimistic. In our experiments, we observe , which may actually be due to the lower bound.

3.3 Generalized Eigenproblems

The stability estimates in Definition 6 are closely linked to properties of the GEVP solved in each subdomain. This subsection establishes the necessary results. In this section, and denote generic symmetric and positive semi-definite bilinear forms on a -dimensional vector space . For a symmetric and positive semi-definite bilinear form denote by

the kernel of . In the following, we make the important assumption (to be verified below for the bilinear forms in (7)) that the bilinear forms satisfy

Definition 7 (Generalized Eigenvalue Problem)

We call with , an eigenpair of , if either and


or and . For such a pair, is called an eigenvalue and an eigenvector of . The collection of all eigenvalues (counted according to their geometric multiplicity) is called the spectrum of .

Lemma 5

The generalized eigenvalue problem (12) is non-defective, i.e. it has a full set of eigenvectors with either or .


Assumption (11) implies that is a symmetric and positive definite bilinear form. Employing a spectral transformation yields

with a symmetric and positive definite bilinear form on the right. Now standard spectral theory for this problem shows that there exists a full set of eigenvectors with real nonnegative eigenvalues. It also follows that . Now corresponds to , and corresponds to , .

From now on we assume that eigenpairs are ordered by size, i.e. . With and the spectrum reads

The eigenvectors , corresponding to the first eigenvalues, can be chosen to be simultaneously -orthonormal and -orthogonal. The remaining eigenvectors are chosen to be -orthogonal and are also -orthogonal, such that forms a basis of .

The projection to the first eigenvectors is defined by


Furthermore, we also introduce the induced semi-norms

Lemma 6

Let be positive semi-definite bilinear forms on with and , , eigenpairs of (12) with , , -orthonormal and , , -orthogonal. Then the projection satisfies the stability estimates


and for


The case is proved in [29, Lemma 2.11]. The proof of (14) in [29] extends without any modification to the case , since the are still -orthogonal and . For the proof of (15) a small modification is necessary. Let be arbitrary. Then,

Remark 3

In [29] the bilinear form in the GenEO method is positive semi-definite but the authors did not include this case in their Lemma 2.11 but rather treat this fact in Lemmata 3.11, 3.16 and 3.18 for a special case. We think the assumption is more natural and easy to prove in applications. In [3, Lemma 2.3] the authors consider the GEVP for two symmetric positive semi-definite matrices , with nontrivial intersection . This is not required in our analysis in the variational setting. However, in the implementation a basis for the space is required on each level. It turns out, it is prohibitively expensive to construct such a basis on levels and only a generating system is available. This then results in a GEVP with . We refer to Section 4.2 how to overcome this problem.

3.4 Application to Continuous Galerkin

In this section, we consider the application to the solution of the scalar elliptic boundary value problem equationparentequation

on (16c)

with Dirichlet and Neumann boundary conditions in a domain . The diffusion coefficient is symmetric and positive definite with eigenvalues bounded uniformly from above and below for all .

We discretize (16) with conforming finite elements on simplicial or hexahedral meshes [14] resulting in the weak formulation


The local bilinear forms on the left side of the GEVP (7) are then defined as


The following Lemma shows that the define a symmetric positive semi-definite splitting and the strengthened triangle inequality under the square holds.

Lemma 7

Let the domain decomposition satisfy definition 2 with . Then the local bilinear forms (17) satisfy definitions 4 and 5 with .


Follows from any being in at most subdomains on any level .

Extending the results in [29], we show local stability (Def. 6) for three different right hand sides in the GEVPs: equationparentequation



The choice (18b) defines the original GenEO method fro