Multispace and Multilevel BDDC

The BDDC (Balancing Domain Decomposition by Constraints) method by Dohrmann [4] is the most advanced method from the BDD family introduced by Mandel [12]. It is a Neumann-Neumann iterative substructuring method of Schwarz type [5] that iterates on the system of primal variables reduced to the interfaces between the substructures. The BDDC method is closely related to the FETI-DP method (Finite Element Tearing and Interconnecting - Dual, Primal) by Farhat et al. [6, 7]. FETI-DP is a dual method that iterates on a system for Lagrange multipliers that enforce continuity on the interfaces, with some “coarse” variables treated as primal, and it is a further development of the FETI method by Farhat and Roux [8]. Polylogarithmic condition number estimates for BDDC were obtained in [13, 14] and a proof that the eigenvalues of BDDC and FETI-DP are actually the same except for eigenvalue equal to one was given in Mandel et al. [14]. Simpler proofs of the equality of eigenvalues were obtained by Li and Widlund [10], and also by Brenner and Sung [3], who also gave an example when BDDC has an eigenvalue equal to one but FETI-DP does not. In the case of many substructures, solving the coarse problem exactly becomes a bottleneck. However, since the coarse problem in BDDC has the same form as the original problem, the BDDC method can be applied recursively to solve the coarse problem only approximately. This leads to a multilevel form of BDDC in a straightforward manner, see Dohrmann [4]. Polylogarithmic condition number bounds for three-level BDDC (BDDC with two coarse levels) were proved in two and three spatial dimensions by Tu [20, 19].

In this paper, we present a new abstract Multispace BDDC method. The method extends a simple variational setting of BDDC from Mandel and Sousedík [15], which could be understood as an abstract version of BDDC by partial subassembly in Li and Widlund [11]. However, we do not adopt their change of variables, which does not seem to be suitable in our abstract setting. We provide a condition number estimate for the abstract Multispace BDDC method, which generalizes the estimate for a single space from [15]. The proof is based on the abstract additive Schwarz theory by Dryja and Widlund [5]. Many BDDC formulations (with an explicit treatment of substructure interiors, after reduction to substructure interfaces, with two levels, and multilevel) are then viewed as abstract Multispace BDDC with a suitable choice of spaces and operators, and abstract condition number estimates for those BDDC methods follow. This result in turn gives a polylogarithmic condition number bound for Multilevel BDDC applied to a second-order scalar elliptic model problems, with an arbitrary number of levels. A brief presentation of the main results of the paper, without proofs and with a simplified formulation of the Multispace BDDC estimate, is contained in the conference paper [16].

The paper is organized as follows. In Sec. 2 we introduce the abstract problem setting. In Sec. 3 we formulate an abstract Multispace BDDC as an additive Schwarz preconditioner. In Sec. 4 we introduce the settings of a model problem using finite element discretization. In Sec. 5 we recall the algorithm of the original (two-level) BDDC method and formulate it as a Multispace BDDC. In Sec. 6 we generalize the algorithm to obtain Multilevel BDDC and we also give an abstract condition number bound. In Sec. 7 we derive the condition number bound for the model problem. Finally, in Sec. 8, we report on numerical results.

We wish to solve an abstract linear problem

where $X$ is a finite dimensional linear space, $a\left(\cdot,\cdot\right)$ is a symmetric positive definite bilinear form defined on $X$, $f_{X}\in X^{\prime}$ is the right-hand side with $X^{\prime}$ denoting the dual space of $X$, and $\left\langle\cdot,\cdot\right\rangle$ is the duality pairing. The form $a\left(\cdot,\cdot\right)$ is also called the energy inner product, the value of the quadratic form $a\left(u,u\right)$ is called the energy of $u$, and the norm $\left\|u\right\|_{a}=a\left(u,u\right)^{1/2}$ is called the energy norm. The operator $A_{X}:X\mapsto X^{\prime}$ associated with $a$ is defined by

A preconditioner is a mapping $B:X^{\prime}\rightarrow X$ and we will look for preconditioners such that $\left\langle r,Br\right\rangle$ is also symmetric and positive definite on $X^{\prime}$. It is well known that then $BA_{X}:X\rightarrow X$ has only real positive eigenvalues, and convergence of the preconditioned conjugate gradients method is bounded using the condition number

which we wish to bound above.

All abstract spaces in this paper are finite dimensional linear spaces and we make no distinction between a linear operator and its matrix.

To introduce abstract Multispace BDDC preconditioner, suppose that the bilinear form $a$ is defined and symmetric positive semidefinite on some larger space $W\supset X$. The preconditioner is derived from the abstract additive Schwarz theory, however we decompose some space between $X$ and $W$ rather than $X$ as it would be done in the additive Schwarz method: In the design of the preconditioner, we choose spaces $V_{k}$, $k=1,\ldots,M$, such that

We formulate the condition number bound first in the full strength allowed by the proof. The bound used in the rest of this paper will be a corollary.

Dryja and Widlund [5] proved that if there exist constants $C_{0},$ $\omega,$ and a symmetric matrix $\mathcal{E}=(e_{ij})_{i,j=1}^{M}$, such that

then

where $\rho$ is the spectral radius.

Now the idea of the proof is essentially to map the assumptions of the abstract additive Schwarz estimate from the decomposition (7) of the space $X$ to the decomposition (2). Define the spaces

We will show that the preconditioner (3) satisfies (8), where $b_{k}$ is defined by

with the operators $G_{k}:X\rightarrow V_{k}^{\mathcal{M}}$ defined by

First, from the definition of operators $G_{k}$, spaces $X_{k},$ and because $a$ is positive definite on $V_{k}$ by Assumption (1), it follows that $G_{k}x$ and $G_{k}z$ in (12) exist and are unique, so $b_{k}$ is defined correctly. To prove (8), let $v_{k}$ be as in (3) and note that $v_{k}$ is the solution of

Consequently, the preconditioner (3) is an abstract additive Schwarz method and we only need to verify the inequalities (9)–(11). To prove (9), let $u\in X$. Then, with $v_{k}$ from the assumption (4) and with $u_{k}=Q_{k}G_{k}v_{k}$ as in (12), it follows that

Next, let $u_{k}\in X_{k}$. From the definitions of $X_{k}$ and $V_{k}^{\mathcal{M}}$, it follows that there exist unique $v_{k}\in V_{k}^{\mathcal{M}}$ such that $u_{k}=Q_{k}v_{k}$. Using the assumption (5) and the definition of $b_{k}$ in (12), we get

which gives (10). Finally, (6) is the same as (11).   

The next Corollary was given without proof in [16, Lemma 1]. This is the special case of Theorem 3 that will be actually used. In the case when $M=1$, this result was proved in [15].

Let $\Omega$ be a bounded domain in $\mathbb{R}^{d}$, $d=2$ or $3$, decomposed into $N$ nonoverlapping subdomains $\Omega^{s}$, $s=1,...,N$, which form a conforming triangulation of the domain $\Omega$. Subdomains will be also called substructures. Each substructure is a union of Lagrangian $P1$ or $Q1$ finite elements with characteristic mesh size $h$, and the nodes of the finite elements between substructures coincide. The nodes contained in the intersection of at least two substructures are called boundary nodes. The union of all boundary nodes is called the interface $\Gamma$. The interface$~\Gamma$ is a union of three different types of open sets: faces, edges, and vertices. The substructure vertices will be also called corners. For the case of regular substructures, such as cubes or tetrahedrons, we can use standard geometric definition of faces, edges, and vertices; cf., e.g., [9] for a more general definition.

In this paper, we find it more convenient to use the notation of abstract linear spaces and linear operators between them instead of the space $\mathbb{R}^{n}$ and matrices. The results can be easily converted to the matrix language by choosing a finite element basis. The space of the finite element functions on $\Omega$ will be denoted as $U$. Let $W^{s}$ be the space of finite element functions on substructure $\Omega^{s}$, such that all of their degrees of freedom on $\partial\Omega^{s}\cap\partial\Omega$ are zero. Let

and consider a bilinear form arising from the second-order scalar elliptic problem as

Now $U\subset W$ is the subspace of all functions from $W$ that are continuous across the substructure interfaces. We are interested in the solution of the problem (1) with $X=U$,

where the bilinear form $a$ is associated on the space $U$ with the system operator$~A$, defined by

and $f\in U^{\prime}$ is the right-hand side. Hence, (17) is equivalent to

Define $U_{I}\subset U$ as the subspace of functions that are zero on the interface $\Gamma$, i.e., the “interior” functions. Denote by $P$ the energy orthogonal projection from $W$ onto $U_{I}$,

Functions from $\left(I-P\right)W$, i.e., from the nullspace of $P,$ are called discrete harmonic; these functions are $a$-orthogonal to $U_{I}$ and energy minimal with respect to increments in $U_{I}$. Next, let $\widehat{W}$ be the space of all discrete harmonic functions that are continuous across substructure boundaries, that is

In particular,

A common approach in substructuring is to reduce the problem to the interface. The problem (17) is equivalent to two independent problems on the energy orthogonal subspaces $U_{I}$ and $\widehat{W}$, and the solution $u$ satisfies $u=u_{I}+\widehat{u}$, where

The solution of the interior problem (22) decomposes into independent problems, one per each substructure. The reduced problem (23) is then solved by preconditioned conjugate gradients. The reduced problem (23) is usually written equivalently as

where $s$ is the form $a$ restricted on the subspace $\widehat{W}$, and $g$ is the reduced right hand side, i.e., the functional $f$ restricted to the space $\widehat{W}$. The reduced right-hand side $g$ is usually written as

because $a(u_{I},\widehat{v})=0$ by (21). In the implementation, the process of passing to the reduced problem becomes the elimination of the internal degrees of freedom of the substructures, also known as static condensation. The matrix of the reduced bilinear form $s$ in the basis defined by interface degrees of freedom becomes the Schur complement, and (24) becomes the reduced right-hand side. For details on the matrix formulation, see, e.g., [17, Sec. 4.6] or [18, Sec. 4.3].

The BDDC method is a two-level preconditioner characterized by the selection of certain coarse degrees of freedom, such as values at the corners and averages over edges or faces of substructures. Define $\widetilde{W}\subset W$ as the subspace of all functions such that the values of any coarse degrees of freedom have a common value for all relevant substructures and vanish on $\partial\Omega,$ and $\widetilde{W}_{\Delta}\subset W$ as the subspace of all function such that their coarse degrees of freedom vanish. Next, define $\widetilde{W}_{\Pi}$ as the subspace of all functions such that their coarse degrees of freedom between adjacent substructures coincide, and such that their energy is minimal. Clearly, functions in $\widetilde{W}_{\Pi}$ are uniquely determined by the values of their coarse degrees of freedom, and

We assume that

That is the case when $a$ is positive definite on the space $U$, where the problem (1) is posed, and there are sufficiently many coarse degrees of freedom. We further assume that the coarse degrees of freedom are zero on all functions from $U_{I}$, that is,

In other words, the coarse degrees of freedom depend on the values on substructure boundaries only. From (25) and (27), it follows that the functions in $\widetilde{W}_{\Pi}$ are discrete harmonic, that is,

Next, let $E$ be a projection from $\widetilde{W}$ onto $U$, defined by taking some weighted average on substructure interfaces. That is, we assume that

Since a projection is the identity on its range, it follows that $E$ does not change the interior degrees of freedom,

since $U_{I}\subset U$. Finally, we show that the operator $\left(I-P\right)E$ is a projection. From (30) it follows that $E$ does not change interior degrees of freedom, so $EP=P$. Then, using the fact that $I-P$ and $E$ are projections, we get

In the next section, the space $X$ will be either $U$ or $\widehat{W}$.

We show several different ways the original, two-level, BDDC algorithm can be interpreted as multispace BDDC. We consider first BDDC applied to the reduced problem (23), that is, (1) with $X=\widehat{W}$. This was the formulation considered in [14]. Define the space of discrete harmonic functions with coarse degrees of freedom continuous across the interface

Because we work in the space of discrete harmonic functions and the output of the averaging operator $E$ is not discrete harmonic, denote

In an implementation, discrete harmonic functions are represented by the values of their degrees of freedom on substructure interfaces, cf., e.g. [18]; hence, the definition (33) serves formal purposes only, so that everything can be written in terms of discrete harmonic functions without passing to the matrix formulation.

Using the decomposition (32), we can split the solution in the space $\widetilde{W}_{\Gamma}$ into the independent solution of two subproblems: mutually independent problems on substructures as the solution in the space $\widetilde{W}_{\Gamma\Delta}=(I-P)\widetilde{W}_{\Delta}$, and a solution of global coarse problem in the space $\widetilde{W}_{\Pi}$. The space $\widetilde{W}_{\Gamma}$ has a decomposition

the same as the decomposition (32), and Algorithm 7 can be rewritten as follows.

To verify the assumptions of Corollary 4, note that the decomposition (36) is $a-$orthogonal, $a(\cdot,\cdot)$ is positive definite on both $\widetilde{W}_{\Gamma\Delta}$ and $\widetilde{W}_{\Pi}$ as subspaces of $\widetilde{W}_{\Gamma}$ by (26), and $E_{\Gamma}$ is a projection by (31).   

Next, we present a BDDC formulation on the space $U$ with explicit treatment of interior functions in the space $U_{I}$ as in [4, 13], i.e., in the way the BDDC algorithm was originally formulated.

The interior corrections (40) and (44) decompose into independent Dirichlet problems, one for each substructure. The substructure correction (42) decomposes into independent constrained Neumann problems, one for each substructure. Thus, the evaluation of the preconditioner requires three problems to be solved in each substructure, plus solution of the coarse problem (43). In addition, the substructure corrections can be solved in parallel with the coarse problem.

The following proposition will be the starting point for the multilevel case.

First, with $V_{1}=U_{I}$, the definition of $v_{1}$ in (3) with $k=1$ is identical to the definition of $u_{I}$ in (40), so $u_{I}=v_{1}$.

Next, consider $w_{\Delta}\in\widetilde{W}_{\Delta}$ defined in (42). We show that $w_{\Delta}$ satisfies (3) with $k=2$, i.e., $v_{2}=w_{\Delta}$. So, let $z_{\Delta}\in\widetilde{W}_{\Delta}$ be arbitrary. From (42) and (41),

Now from the definition of $u_{I}$ by (40) and the fact that $PEz_{\Delta}\in U_{I}$, we get

and subtracting (49) from (48) gives

because $a\left(u_{I},\left(I-P\right)Ez_{\Delta}\right)=0$ by orthogonality. To verify (3), it is enough to show that $Pw_{\Delta}=0;$ then $w_{\Delta}\in(I-P)\widetilde{W}_{\Delta}=V_{2}$. Since $P$ is an $a$-orthogonal projection, it holds that

where we have used $EU_{I}\subset U_{I}$ following the assumption (30) and the equality

for any $z_{I}\in U_{I}$, which follows from (41) and (40). Since $a$ is positive definite on $\widetilde{W}\supset U_{I}$ by assumption (26), it follows from (50) that $Pw_{\Delta}=0$.

In exactly the same way, from (43) – (45), we get that if $w_{\Pi}\in$ $\widetilde{W}_{\Pi}$ is defined by (43), then $v_{3}=w_{\Pi}$ satisfies (3) with $k=3$. (The proof that $Pw_{\Pi}=0$ can be simplified but there is nothing wrong with proceeding exactly as for $w_{\Delta}$.)

Finally, from (44), $v_{I}=P\left(Ew_{\Delta}+Ew_{\Pi}\right)$, so

It remains to verify the assumptions of Corollary 4.

First, the spaces $\widetilde{W}_{\Pi}$ and $\widetilde{W}_{\Delta}$ are $a$-orthogonal by (25) and, from (27),

thus $\left(I-P\right)\widetilde{W}_{\Delta}\perp_{a}\widetilde{W}_{\Pi}$. Clearly, $\left(I-P\right)\widetilde{W}_{\Delta}\perp_{a}U_{I}$. Since $\widetilde{W}_{\Pi}$ consists of discrete harmonic functions from (28), so $\widetilde{W}_{\Pi}\perp_{a}U_{I}$, it follows that the spaces $V_{i}$, $i=1,2,3$, given by (46), are $a$-orthogonal.

Next, $\left(I-P\right)E$ is by (31) a projection, and so are the operators $Q_{i}$ from (47).

It remains to prove the decomposition of unity (14). Let

and let

From (51), $w_{\Delta}+w_{\Pi}\in U$ since $u^{\prime}\in U$ and $u_{I}\in U_{I}\subset U$. Then $E\left(w_{\Delta}+w_{\Pi}\right)=w_{\Delta}+w_{\Pi}$ by (29), so

because both $w_{\Delta}$ and $w_{\Pi}$ are discrete harmonic.   

The next Theorem shows an equivalence of the three Algorithms introduced above.

Let us now turn to Algorithm 11. Let the residual $r\in U$ be written as $r=r_{I}+r_{\Gamma}$, where $r_{I}\in U_{I}^{\prime}$ and $r_{\Gamma}\in\widehat{W}^{\prime}$. Taking $r_{\Gamma}=0$, we get $r=r_{I}$, and it follows that $r_{B}=u_{B}=v_{I}=0$, so $u=u_{I}$. On the other hand, taking $r=r_{\Gamma}$ gives $u_{I}=0$, $r_{B}=r_{\Gamma}$, $v_{I}=Pu_{B}$ and finally $u=\left(I-P\right)E(w_{\Delta}+w_{\Pi})$, so $u\in\widehat{W}$. This shows that the off-diagonal blocks of the preconditioner $M$ are zero, and therefore it is block diagonal

Next, let us take $u=u_{I},$ and consider $r_{\Gamma}=0$. The algorithm returns $r_{B}=u_{B}=v_{I}=0$, and finally $u=u_{I}$. This means that $M_{II}\mathcal{A}_{II}u_{I}=u_{I},$ so $M_{II}=\mathcal{A}_{II}^{-1}$. The operator $A:U\rightarrow U^{\prime}$, and the block preconditioned operator $MA:r\in U^{\prime}\rightarrow u\in U$ from Algorithm 11 can be written, respectively, as