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  and lead to the development of algebraic multigrid (AMG) methods . Rigorous convergence bounds for aggregation-based AMG applied to nonsingular symmetric M-matrices with nonnegative row sums are provided in . In the context of overlapping DD methods, it was shown in  based on weighted Poincaré inequalities  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, a local eigenvalue problem involving the Dirichlet-to-Neumann map was introduced, and later analysed in  based on . 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 . 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, , 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 . An extension of the GenEO approach to nonoverlapping domain decomposition methods as well as non-selfadjoint problems can be found in .
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  has been proposed in . 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  based on the SPSD splitting introduced in .
In this paper, we formulate a natural extension of GenEO  to an arbitrary number of levels. In contrast to , our method is formulated in a variational framework based on subspace correction . The convergence theory of  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 . Condition number bounds are derived for an iterated two-level method as in  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  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 , 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:
On the finest level , decompose into overlapping sets
by first partitioning into nonoverlapping sets of elements using a graph partitioner such es ParMetis  and then adding a user defined overlap in terms of layers of elements.
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
The subdomains are now defined from the mesh decomposition as
On the coarsest level 0 we employ always only one subdomain.
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.
A partition of unity on level is a family of operators such that
is zero on for all and
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:
Set . Set .
For each subdomain solve a GEVP
with appropriately defined local bilinear forms and
detailed below. Let eigenvalues and corresponding eigenvectors be ordered s.t..
With a user defined parameter define the coarse space as
Set . If goto step 2, otherwise stop.
The subspaces for the subspace correction method based on are then given by
together with .
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
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)
We say the multilevel domain decomposition admits a level-wise coloring with colors, if for each level there exists a map such that
The multilevel domain decomposition is called admissible, if on each level every interior face is interior to at least one subdomain.
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
If the subspaces and , , admit a stable splitting with constant , parallel subspace correction satisfies the lower bound
See e.g. .
3.2 Abstract Schwarz Theory for Spectral Domain Decomposition
We now generalize the theory in  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 :
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
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
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
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 .
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
Let be positive semi-definite bilinear forms on with and , , eigenpairs of (12) with , , -orthonormal and , , -orthogonal. Then the projection satisfies the stability estimates
In  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
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 .
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.
Follows from any being in at most subdomains on any level .