Citation
Adaptive-multilevel BDDC

Material Information

Title:
Adaptive-multilevel BDDC
Creator:
Sousedik, Bedrich
Place of Publication:
Denver, Colo.
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
xvii, 109 leaves : ; 28 cm

Thesis/Dissertation Information

Degree:
Doctorate ( Doctor of Philosophy)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Mathematical and Statistical Sciences, CU Denver
Degree Disciplines:
Applied Mathematics
Committee Chair:
Mandel, Jan
Committee Members:
Bennethum, Lynn S.
Burda, Pavel
Welch, Samuel W. J.

Subjects

Subjects / Keywords:
Decomposition method ( lcsh )
Boundary value problems ( lcsh )
Boundary value problems ( fast )
Decomposition method ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 100-109).
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Bedrich Sousedik.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
655258948 ( OCLC )
ocn655258948
Classification:
LD1193.L622 2010d S68 ( lcc )

Full Text
ADAPTIVE MULTILEVEL BDDC
by
Bedrich Sousedfk
M. Eng., Czech Technical University in Prague, 2001
M.S., University of Colorado Denver, 2008
Ph.D., Czech Technical University in Prague, 2008
A thesis submitted to the
University of Colorado Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
2010


This thesis for the Doctor of Philosophy
degree by
Bedrich Sousedik
has been approved
by
Date


Sousedfk, Bedrich (Ph.D., Applied Mathematics)
Adaptive Multilevel BDDC
Thesis directed by Professor Jan Mandel
ABSTRACT
The Balancing Domain Decomposition by Constraints (BDDC) method
by Dohrmann (2003) is the most advanced method from the Balancing fam-
ily of iterative substructuring methods by Mandel (1993) for the solution of
large systems of linear algebraic equations arising from discretizations of ellip-
tic boundary value problems. The method is closely related to FETI-DP by
Farhat et. al. (2001), and it is the same as other two methods proposed inde-
pendently by Fragakis and Papadrakakis (2003) and by Cros (2003).
Because these are two-level methods, solving the coarse problem exactly
becomes a bottleneck when the number of substructures becomes large. The
coarse problem in BDDC has the same structure as the original problem, so it
is straightforward to apply the BDDC method recursively to solve the coarse
problem only approximately. In the first part we formulate a new family of
abstract Multispace BDDC methods and give a condition number bound from
the abstract additive Schwarz preconditioning theory. The Multilevel BDDC is
then treated as a special case of the Multispace BDDC, and it is also shown that
the original, two-level, BDDC can be written as a multispace method.


In the second part we propose a method for adaptive selection of the coarse
space for the original two-level BDDC method. The method works by adding
coarse degrees of freedom constructed from eigenvectors associated with intersec-
tions of selected pairs of adjacent substructures. It is assumed that the starting
coarse degrees of freedom are already sufficient to prevent relative rigid body
motions in any selected pair of adjacent substructures. A heuristic indicator of
the condition number is developed and a minimal number of coarse degrees of
freedom is added to decrease the indicator under a given threshold.
In the third part we combine the advantages of both approaches to propose
a new method called Adaptive Multilevel BDDC that preserves both parallel
scalability with increasing number of subdomains and very good convergence
properties. Performance of the method is illustrated by several numerical ex-
amples in two and three spatial dimensions.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Jan Mandel


DEDICATION
Dedicated to Alyssa Heberton-Morimoto.


ACKNOWLEDGMENT
I would like to give my deepest thank to my advisor Professor Jan Mandel.
He was always a good source of interesting new problems, and let me work on
topics that I liked. Without his guidance and mentorship this thesis would have
been impossible. A large debt of gratitude goes to the author of the BDDC
method Dr. Clark R. Dohrmann for bringing our attention to its multilevel
extension, and to Professors Marian Brezina, Jaroslav Kruis, Jaroslav Novotny
and Jakub Sfstek who provided the data for numerical experiments and helped
me with visualization of finite element meshes. In particular, I have been truly
impressed by Marians own visualization package My Vis. I am grateful to Prof.
Pavel Burda and Dr. Chris Harder for proof-reading of a draft of this thesis.
I would also like to thank the faculty and my fellow students at the Uni-
versity of Colorado Denver for creating excellent working conditions. Last but
not least I thank my parents for bringing me up the way they did and for their
continuing support. Finally, I thank my family and friends for support and
distraction from science.
The doctoral dissertation has been kindly supported by the National Sci-
ence Foundation under grants CNS-0325314, CNS-0719641, and DMS-0713876.
The research completed during my visits at the Institute of Thermomechanics,
Academy of Sciences of the Czech Republic has been supported by the Grant
Agency of the Czech Republic GA CR 106/08/0403.


CONTENTS
Figures .............................................................. viii
Tables................................................................ xiii
Chapter
1. Introduction......................................................... 1
2. Problem Setting and Preliminaries................................... 11
2.1 Additive Schwarz preconditioners................................... 12
2.2 Abstract variational BDDC preconditioner........................... 14
3. Abstract Multispace BDDC............................................ 17
3.1 Finite element setting ............................................ 20
3.2 Two-level BDDC as Multispace BDDC.................................. 25
4. Multilevel BDDC..................................................... 36
5. Adaptive Coarse Degrees of Freedom.................................. 48
5.1 Preconditioned LOBPCG ............................................. 59
6. Adaptive Multilevel BDDC.......................................... 65
6.1 Implementation remarks............................................. 67
7. Numerical Examples.................................................. 70
7.1 Two-dimensional results............................................ 70
7.2 Three-dimensional results.......................................... 77
8. Conclusion.......................................................... 98
References............................................................. 100
vii


FIGURES
Figure
3.1 An example of a domain decomposition for the two-level (top) and
the three-level (bottom) BDDC methods............................. 35
4.1 Spaces, embeddings and projections in the Multilevel BDDC. ... 39
5.1 Comparison of the non-preconditioned (top) vs. preconditioned
LOBPCG (bottom) for one of the faces with large jumps in coef-
ficients of the composite cube problem, cf. Chapter 7 and Fig. 7.2.
Estimated eigenvalue errors are in the panels on the left-hand side,
and Euclidean norms of residuals for different eigenpairs are in the
panels on the right-hand side. We see that LOBPCG without a
preconditioner showed essentially no convergence, and with the pre-
conditioner we have reached convergence in less than 30 iterations. 61
5.2 A pair of substructures of the mining reel problem from Figure 7.8,
obtained from the automatic decomposition by METIS 4.0. We see
that one of the substructures has 4 spurious rigid body modes. . .
vni
62


5.3 Comparison of the non-preconditioned (top) vs. preconditioned
LOBPCG (bottom) for the detection of spurious rigid-body modes
of the pair of subdomains from Fig. 5.2. Estimated eigenvalue er-
rors are in the panels on the left-hand side, and Euclidean norms of
residuals for different eigenpairs are in the panels on the right-hand
side. We see that LOBPCG without a preconditioner essentially did
not detect the nullspace, and the application of preconditioner led
to relatively fast detection of the nullspace................... 63
5.4 Convergence of the preconditioned LOBPCG for the pair of subdo-
mains from Fig. 5.2 with the nullspace detection (Fig. 5.3). Esti-
mated eigenvalue errors are in the panels on the left-hand side, and
Euclidean norms of residuals for different eigenpairs are in the panels
on the right-hand side.......................................... 64
7.1 Two remote corners of the two-level decomposition into 48 x 48 (=
2304) subdomains (top), and the decomposition into 9 subdomains
for the three-level method (bottom). The jagged edge from the lower
decomposition level (top) is indicated on the second-level decompo-
sition (bottom) by the thick line............................... 72
IX


7.2 Finite element discretization and substructuring of the cube with
jumps in coefficients, consisting of 107811 degrees of freedom, dis-
tributed into 8 substructures with 30 corners, 16 edges and 15
faces (the bars cut the substructures only through faces). Notably,
similar problems are solved in practice to determine numerically
(anisotropic) properties of composite materials [72]. Courtesy of
Jakub Sfstek........................................................ 78
7.3 Finite element discretization and substructuring of the large cube
with jumps in coefficients, consisting of 823 875 degrees of freedom,
distributed into 512 substructures with 721 corners, 1176 edges and
1 344 faces on the first decomposition level (top), and 4 substructures,
6 corners, one edge and 4 faces on the second decomposition level
(bottom). Courtesy of Jakub Sfstek............................... 81
7.4 Finite element discretization and substructuring of the dam, consist-
ing of 2 006 748 degrees of freedom, distributed into 400 substructures
with 3 990 corners, 3 070 edges and 2 274 faces. Courtesy of Jaroslav
Kruis................................................................. 84
7.5 Correspondence of finite elements and the subdomains on the sec-
ond decomposition level. The dam problem with 3800080 finite
elements, 400 substructures on the first level and 8 substructures on
the second level...................................................... 85
x


7.6 Finite element discretization and substructuring of the dam, consist-
ing of 2 006 748 degrees of freedom, distributed into 1024 substruc-
tures with 10 693 corners, 7 713 edges and 6 182 faces. Courtesy of
Jaroslav Kruis...................................................... 87
7.7 Correspondence of finite elements and the subdomains on the sec-
ond decomposition level. The dam problem with 3800080 finite
elements, 1024 substructures on the first level and 32 substructures
on the second level................................................. 88
7.8 Finite element discretization and substructuring of the mining reel
problem, consisting of 1 739 211 degrees of freedom, distributed into
400 substructures with 4010 corners, 831 edges and 1906 faces.
Courtesy of Jan Lestina, Jaroslav Novotny and Jakub Sfstek. ... 91
7.9 Correspondence of finite elements on the zero decomposition level
and the subdomains on the second decomposition level. Mining reel
with 140 816 finite elements, 400 substructures on the first level and
8 substructures on the second level................................. 92
7.10 Finite element discretization and substructuring of the mining reel
problem, consisting of 1 739 211 degrees of freedom, distributed into
1024 substructures with 7864 corners, 1197 edges, and 3895 faces.
Courtesy of Jan Lestina, Jaroslav Novotny and Jakub Sfstek. ... 92
7.11 Correspondence of finite elements on the zero decomposition level
and the subdomains on the second decomposition level. Mining reel
with 140 816 finite elements, 1024 substructures on the first level and
32 substructures on the second level.............................. 93
xi


7.12 Correspondence of finite elements on the zero decomposition level
and the subdomains on the second decomposition level. The bridge
problem with 880 000 finite elements, 1 024 substructures on the first
level and 8 substructures on the second level (a scaled view).... 95
7.13 Finite element discretization of the bridge construction with 3 173 664
degrees of freedom, distributed into 1 024 subdomains with 6051 cor-
ners, 2 099 edges and 3 034 faces. Courtesy of Jaroslav Kruis. ... 96
Xll


TABLES
Table
7.1 Results for the planar elasticity from Fig. 7.1 obtained using non-
adaptive 2-level method. Constraints are corners, or corners and
arithmetic averages over edges, denoted as c, c+f, resp., and Nc
is number of constraints (coarse degrees of freedom), C is size of
the coarse problem related to size of a subdomain problem, k is the
approximate condition number computed from the Lanczos sequence
in conjugate gradients, and it is the number of iterations for relative
residual tolerance 10-8............................................. 73
7.2 Results for the planar elasticity from Fig. 7.1 obtained using non-
adaptive 3-level method. Nc is the number of coarse degrees of
freedom on the first (+ the second) decomposition level, C is the
relative coarsening with respect to the size of substructures on the
first level (the size of the coarse problem for the three-level method
is negligible). Other headings are the same as in Table 7.1......... 73
7.3 The largest eigenvalues \st,k of the local eigenvalue problems for sev-
eral pairs of subdomains s and t of the 2-level elasticity problem from
Fig. 7.1 (top), with (s,t) = (2,50) being the jagged edge.
xm
74


7.4 The largest eigenvalues of the local problems for several pairs
of subdomains s, t on the level i = 2, cf. Fig. 7.1 (lower panel), with
t 2 (top) and with r = 10 (bottom). The jagged edge is between
subdomains 2 and 5.................................................. 75
7.5 Results for the planar elasticity from Fig. 7.1 obtained using the
adaptive 2-level method. Here, r is condition number target, oj is
condition number indicator, and the other headings are the same as
in Table 7.1.......................................................... 76
7.6 Results for the planar elasticity from Fig. 7.1 obtained using the
adaptive 3-level method. Headings are the same as in Tables 7.2
and 7.5............................................................... 76
7.7 Results for the cube from Fig. 7.2 obtained using a preconditioning
by incomplete Cholesky factorization. The global stiffness matrix
has size 107 811 with 7 737 641 nonzeros. Here, nnz(R) is the number
of nonzeros in the upper triangular Cholesky factor R, fill-in is the
relative fill-in computed as 2 times nnz(R) divided by the number of
nonzeros in the global stiffness matrix, k is the approximate condition
number computed from the Lanczos sequence in conjugate gradients,
it is number of iterations for relative residual tolerance 10-8. With
the zero-level of fill-in the method did not converge............... 78
xiv


7.8 Results for the cube from Fig. 7.2 obtained using the non-adaptive
2-level method. Constraints are corners, or corners and arithmetic
averages over edges and faces denoted as c, c+e, c+e+f resp., and
c+e+f (3eigv), corresponding to corner constraints, arithmetic av-
erages, and 3 weighted averages over each face obtained using the
adaptive method. Next, Nc is the number of constraints, C is the
size of the coarse problem related to size of a subdomain problem,
k is the approximate condition number computed from the Lanczos
sequence in conjugate gradients, it is number of iterations for relative
residual tolerance 10-8......................................... 79
7.9 Results for the cube from Fig. 7.2 obtained using the adaptive 2-
level method. Here, r is the condition number target, u is the con-
dition number indicator. An approximate number of nonzeros of the
Cholesky factor of a substructure problem is 8 500 000 for all values
of r, and the number of nonzeros in the Cholesky factor of the coarse
problem is denoted by nnz(c). The other headings are the same as
in Table 7.8.................................................... 79
7.10 Results for the cube from Fig. 7.3 obtained using the non-adaptive
2- level method. The headings are the same as in Table 7.8...... 80
7.11 Results for the cube from Fig. 7.3 obtained using the adaptive 2-level
method. The headings are the same as in Table 7.9............... 80
7.12 Results for the cube from Fig. 7.3 obtained using the non-adaptive
3- level method. The headings are the same as in Tables 7.2 and 7.8. 82
xv


7.13 Results for the cube from Fig. 7.3 obtained using the adaptive 3-level
method. The headings are the same as in Tables 7.6 and 7.9....... 82
7.14 Results for the dam (Fig. 7.4, 400 substructures) obtained using
the non-adaptive 2-level method. The headings are the same as in
Table 7.8..................................................... 85
7.15 Results for the dam from Figs. 7.4 and 7.5 (400 + 8 substructures)
obtained using the non-adaptive 3-level method. The headings are
the same as in Tables 7.2 and 7.8................................ 86
7.16 Results for the dam (Fig. 7.6, 1024 substructures) obtained using
the non-adaptive 2-level method. The headings are the same as in
Table 7.8..................................................... 86
7.17 Results for the dam from Figs. 7.6 and 7.7 (1024 + 32 substructures)
obtained using the non-adaptive 3-level method. The headings are
the same as in Tables 7.2 and 7.8................................... 86
7.18 Results for the mining reel (Fig. 7.8, 400 substructures) obtained
using the adaptive 2-level method. The headings are the same as in
Table 7.9..................................................... 90
7.19 Results for the mining reel from Figs. 7.8 and 7.9 (400 and 8 sub-
domains) obtained using teh adaptive 3-level method. The headings
are the same as in Tables 7.6 and 7.9............................ 90
7.20 Results for the mining reel (Fig. 7.10, 1 024 substructures) obtained
by the non-adaptive 2-level method. The headings are the same as
in Table 7.8........................................................ 93
xvi


7.21 Results for the mining reel (Fig. 7.10, 1024 substructures) obtained
using the adaptive 2-level method. The headings are the same as in
Table 7.9....................................................... 94
7.22 Results for the mining reel from Figs. 7.10 and 7.11 (1024 + 32 sub-
structures) obtained using the adaptive 3-level method. The head-
ings are the same as in Tables 7.6 and 7.9...................... 94
7.23 Results for the bridge construction from Fig. 7.13 obtained using the
adaptive 2-level method. The headings are the same as in Table 7.9. 97
7.24 Results for the bridge construction from Figs. 7.13 and 7.12 (1024 +
8 substructures) obtained using the adaptive 3-level method. The
headings are the same as in Tables 7.6 and 7.9.................. 97
xvu


1. Introduction
A dramatic increase in availability of multicore, parallel computers over the
last few decades has motivated a progressive research in many areas, including
domain decomposition methods. Out of many possible interpretations of the
domain decomposition cf., e.g., monographs [74, 77] it will be in our context
understood as a separation of a physical domain into subdomains, called alter-
natively (with the same meaning) substructures. These methods then aim to
allow for an efficient and scalable solution of systems of linear equations arising
from discretizations of partial differential equations (PDEs), in particular by
finite element methods. In this thesis we will restrict our attention to elliptic
problems such as Laplace equation or problems of linear elasticity.
The important class of domain decomposition methods that we will focus on
is known as iterative substructuring. From the point of view of linear algebra, the
solution of a large problem is replaced by repeated solutions of a number of inde-
pendent subproblems, to be solved in parallel. The algorithms are formulated as
preconditioned Krylov iterative methods such as conjugate gradients (CG), or
GMRES. All these methods access the specific problem through only two sub-
routines, namely the system matrix-vector multiplication, and preconditioner
matrix-vector multiplication. There is often a pre-processing step consisting of
a transformation of variables before an iterative method is applied; in that case,
subroutines to pre-process the right-hand side and to recover the solution also
need to be provided.
1


The following method components form the system and the preconditioner
presented to the iterative method:
The solution space or, as we will see later, some larger space is split into
a number of subspaces, and independent problems to be solved in parallel
are set up in each one of them, and solved in each iteration.
Correspondence between the subspaces and the solution space is estab-
lished to build a global solution from the local solutions.
A global coarse problem with one or as few as possible variables per sub-
space is set up and then solved in each iteration to coordinate the solution
between all subspaces at once.
Careful splitting into subspaces and preconditioning in the subspaces are
needed for scalability with subdomain size. A coarse problem is necessary for
scalability with increasing number of subdomains, and hence, processors. We
now briefly review some domain decomposition methods important in the con-
text of the BDDC method the main topic of this thesis. The comprehen-
sive overview of the field has been the subject of numerous monographs, e.g.,
[47, 74, 77, 81, 87, 91], cf., also the standard finite element monograph [5]. An
important role in the present development is also played by another thesis of
the author [79], where some of the most popular methods were formulated and
compared under so-called minimalist assumptions. This revealed that some of
the methods are closely related, equivalent or even the same. The theory has
been used in the final chapter to develop an adaptive method that significantly
improves convergence of two most advanced methods, BDDC and FETI-DP.
2


We note that the basic idea of the adaptive method is the same in the two
theses. However, the implementation in [79] uses global spaces, projections,
and a change of variables to preserve sparsity of global operators, whereas the
formulation presented here builds explicitly all subspaces (including the coarse).
Perhaps the earliest primal domain decomposition methods were based on
the concept of substructuring. The domain is divided into non-overlapping sub-
structures, and there is one subspace formed by the degrees of freedom within
a substructure. The subspaces overlap in the degrees of freedom on substruc-
ture interfaces. The internal degrees of freedom in each substructure are then
eliminated; this alone results in a significant decrease of the condition number,
typically from 0(h~2) to 0(h.~l). So even simple diagonal preconditioners and
no coarse problem result in good practical parallel methods [35]. The system
matrix-vector multiplication is implemented by solving an independent Dirich-
let problem in each substructure, which is the same as multiplication by the
Schur complement in each substructure. More sophisticated diagonal precondi-
tioners need to know the diagonal entries of the Schur complements, which are
not available, but can be estimated from the matrix-vector multiplication [11],
Asymptotically optimal (for small h) diagonal preconditioning of the problem
reduced to substructure interfaces, along with various coarse problems, were in-
troduced in pioneering theoretically oriented papers [2, 3, 19, 89], The coarse
problems are constructed from piecewise linear or piecewise constant functions
on substructure scale. These methods have condition numbers that, for geomet-
rically regular problems, provably do not grow faster than log2^, where H is
the typical substructure size and h is the mesh step or element size.
3


Methods preconditioning the system reduced to interfaces by a weighted sum
of inverses of the reduced substructure matrices [14, 15, 30] have become known
as Neumann-Neumann methods [20], because multiplication of a vector by the
inverse is the same as numerically solving the differential equations on both sides
of a substructure interface with the Neumann boundary conditions. A scalable
version of the Neumann-Neumann method requires a suitable coarse problem.
The Balancing Domain Decomposition (BDD) by Mandel [53] constructs the
coarse problem from the natural nullspace of the problem (constants for scalar
problems, rigid body modes for elasticity). The BDD has become quite popular,
because it can be formulated purely algebraically, requiring only substructure
matrix-vector multiplications and the substructure nullspace vectors. For ex-
ample, in application of BDD to mixed discretization of a flow problem [12],
the matrix-vector multiplications are provided by pressure to flux and flux to
pressure mappings on substructure interfaces. The BDD method and its theory
were further developed for problems with coefficient jumps [55, 57], and also for
plates and shells [48, 49].
Currently the most advanced version of the primal substructuring method is
the Balancing Domain Decomposition by Constraints (BDDC) by Dohrmann [17].
It uses (similarly to [49]) a coarse space defined by energy minimization on each
substructure separately, subject to given values of coarse degrees of freedom
(such as values at substructure corners or edge or face averages), which act
as constraints, but the BDDC method uses a careful additive approach: the
substructure spaces and the coarse space form an energy orthogonal decomposi-
tion of a space of vectors that are not continuous across substructure interfaces;
4


the only approximation comes from enforcing the continuity across substructure
interfaces by averaging. This results in lower condition numbers and faster con-
vergence than other methods. The first proofs that BDDC is asymptotically
optimal (the condition number does not grow faster than log2^) were given by
Mandel, Dohrmann and Tezaur in [58, 59], Unlike other primal substructuring
methods, BDDC requires only the substructure matrices.
The dual substructuring methods enforce the continuity of the solution on
substructure interfaces by Lagrange multipliers; for the case of 2 substructures,
see [30]. In general, some substructure matrices will have a nontrivial nullspace,
and guaranteeing that all local problems have right-hand sides orthogonal to the
nullspace gives the original Finite Element Tearing and Interconnecting (FETI)
method by Farhat and Roux [26], later called FETI-1. Farhat, Mandel and
Roux have shown that resolving the nullspace acts as a natural coarse prob-
lem [25], but the original FETI method still had only a diagonal preconditioner,
with resulting lack of scalability. An asymptotically optimal preconditioner
involving a Dirichlet problem in each substructure and the first proof of the
polylogarithmic growth of the condition number were given in [66]. Further
developments of the FETI method include more complicated coarse spaces for
plates and shells [24, 22], which also proved to have polylogarithmic bound on
the condition number growth [68]. Currently the most advanced method is Fi-
nite Element Tearing and Interconnecting Dual, Primal (FETI-DP) by Farhat
et al. [23, 21]. The method enforces continuity of the values of the solution
at substructure corners and, in 3D, of the values of face averages, cf. [23, 44],
by choosing those values as (primal) coarse degrees of freedom, and enforcing
5


the equality of the rest of the degrees of freedom on substructure interfaces by
Lagrange multipliers. All original and coarse degrees of freedom are eliminated
(this involves setting up a coarse problem and solving it in every iteration) and
the remaining dual system in terms of the Lagrange multipliers is solved itera-
tively. Polylogarithmic condition number bounds for FETI-DP were first proved
by Mandel and Tezaur [67] and generalized to the case of coefficient jumps be-
tween substructures by Klawonn, Widlund and Dryja [43],
The FETI method is conceptually dual to the BDD. In the FETI method,
evaluation of the system matrix-vector multiplication involves solving a Neu-
mann problem on each substructure, and the preconditioner involves solving a
Dirichlet problem on each substructure. In the BDD, the role of the Dirichlet
and the Neumann problems is reversed. The rest of the algebra is somewhat dif-
ferent so the base methods are not dual in any exact sense. However, algebraic
relations between the mesh transfer operators in primal and dual methods were
proved by Rixen at al. [75] and Klawonn and Widlund [41] have reduced the
asymptotic analysis of BDD and FETI to essentially a single inequality. Fra-
gakis and Papadrakakis [29] then classified common components of BDD and
FETI type methods and observed numerically that the preconditioned operator
in one version [1] of the FETI method, with certain improvements to make it
more robust, has the same eigenvalues as the BDD method, and proposed several
primal methods derived from FETI class methods. Proof that the eigenvalues
of BDDC and FETI-DP are identical except for eigenvalue equal to one was
given by Mandel, Dohrmann and Tezaur [59]. These proofs were followed by
a surge of important simplified proofs, alternative frameworks, and extensions
6


to various problems by the best analysts in the substructuring field and their
students and coworkers, e.g., [4, 6, 7] and [28, 40, 43, 52, 51, 69]. The most
popular methods were also analyzed by Mandel and Sousedik under so-called
minimalist assumptions [62, 79, 80]. This setting allowed to show that in the
case of corner constraints only, methods identical to BDDC were independently
derived as primal versions of FETI-DP by Cros [13] and by Fragakis and Pa-
padrakakis [29]. The BDDC and, equivalently, FETI-DP are quite robust. It
can be proved that the condition number remains bounded even for large classes
of subdomains with rough interfaces in 2D [39, 90] as well as in many cases of
strong discontinuities of coefficients, including some configurations when the
discontinuities cross substructure boundaries [70, 71].
Solving the coarse problem exactly in the original BDDC method becomes
a bottleneck as the number of unknowns and, in particular, the number of
substructures gets too large. Since the coarse problem in BDDC, unlike in the
FETI-DP, has the same structure as the original problem, it is straightforward
to apply the method recursively to solve the coarse problem only approximately
[17]. The original, two-level, BDDC has been extended into three-levels by Tu
[83, 82, 85, 84] and into a general multilevel method by Mandel, Sousedik and
Dohrmann [63, 64], However, the abstract condition number bounds reveal
deteriorating convergence with increasing number of levels. We also note that
recently Tu extended BDDC into three-level methods for saddle point problems
[86], and with Kim for mortar discretizations [36], Another research direction
includes applications of inexact substructure solvers in the BDDC methods [18,
52], and inexact coarse problem solvers in the FETI methods [33, 37, 38].
7


Moreover, despite their robustness, the condition number of the BDDC and,
equivalently, FETI-DP does deteriorate in many situations of practical impor-
tance and an adaptive method is warranted. Enriching the coarse space so
that the iterations run in a subspace devoid of difficult modes has been a
long-standing trick in iterative substructuring methods and used, e.g., in the
development of BDD and FETI for plates from the base BDD and FETI meth-
ods [24, 48, 49, 68]. Methods that build a coarse space adaptively from lo-
cal eigenvalue calculations were also devised in a number of other contexts
[8, 27, 54, 56, 73]. Adaptive enrichment for BDDC and FETI-DP was pro-
posed by Mandel and Sousedfk in [60, 61], with the added coarse functions built
from eigenproblems based on adjacent pairs of substructures in 2D. However,
the adaptive method was formulated in terms of FETI-DP operators and it was
quite complicated. Later, the algorithm has been developed directly in terms of
BDDC operators and extended to 3D by Mandel, Sousedfk and Sfstek [65, 79],
resulting in a much simplified formulation and implementation. This implemen-
tation framework operates on global matrices, builds no explicit coarse problem,
and gets much of its parallelism through the direct solver used to solve an aux-
iliary decoupled system. To preserve sparsity, the authors used a variant of the
change of variables from [51], extended to an arbitrary number of constraints.
It has been shown on numerical examples that the heuristic eigenvalue-based
estimates work reasonably well and that the adaptive approach can result in the
concentration of computational work in a small troublesome part of the prob-
lem, which leads to a good convergence behavior at a small added cost. The
only requirement for the adaptive algorithm is that there is a sufficient num-
8


ber of corner constraints to prevent rigid body motions between any pair of
adjacent substructures. This requirement has been recognized recently in other
contexts [9, 50], and in the context of BDDC by Burda et al. [10].
The main goal of this thesis is to combine the adaptive and multilevel ap-
proaches to the BDDC method in order to develop its variant that would preserve
parallel scalability with an increasing number of subdomains and also show ex-
cellent convergence properties. Some of the material has already been published.
This thesis is organized as follows. In Chapter 2 we establish the notation and
introduce problem settings and preliminaries. Following [63, 64], in Chapter 3
we formulate the Abstract Multispace BDDC. In Chapter 4 we then present
the Multilevel BDDC as a particular instance of the Multispace BDDC. We also
derive an abstract condition number bound. Chapter 5 is based on a series of pa-
pers by Mandel, Sousedfk and Sfstek [61, 65, 79, 78], except for Section 5.1 which
is new. We describe the adaptive selection of constraints in terms of BDDC op-
erators only ([65, 79]). However, the presented formulation includes an explicit
coarse space correction, outlined for 2D in [78], which allows for a multilevel
extension of the adaptive algorithm. Chapter 6 is new and contains the main
result of this thesis. We combine the adaptive and multilevel approaches and
formulate the new method called Adaptive Multilevel BDDC which allows for
the adaptive selection of constraints on each decomposition level. Numerical
experiments are presented in Chapter 7. It appears that the adaptive algorithm
is able to effectively detect troublesome parts on each decomposition level and
decrease the number of iterations at a small added cost. Finally, Chapter 8
contains the summary and concluding remarks.
9


The Multispace and Multilevel BDDC (Chapters 3-4) resulted from a joint
research with Jan Mandel and Clark R. Dohrmann. The first idea of the adaptive
algorithm (Chapter 5) can be traced back to a joint paper with Jan Mandel [61].
Its extension into 3D was originally published in the authors first thesis [79].
However, this work contains the first publication of this formulation with an
explicit coarse space, along with a preconditioner for LOBPCG (Section 5.1).
The main result of the thesis, the Adaptive-Multilevel BDDC algorithm is also
an original contribution of the author; and only a 2D outline of the algorithm
has been recently submitted in the conference proceedings [78]. The 3D version
of the algorithm appears here for the first time.
10


2. Problem Setting and Preliminaries
We wish to solve an abstract linear problem
u E X : a(u,v) = (fx,v), Vv E X, (2.1)
where X is a finite dimensional linear space, a(-,-) is a symmetric positive
definite bilinear form defined on X, fx X' is the right-hand side with X'
denoting the dual space of X, and (, ) is the duality pairing. The form a (, )
is also called the energy inner product, the value of the quadratic form a (u, u)
is called the energy of u, and the norm ||u||a = [a (u,u)\1^2 is called the energy
norm. The operator Ax : X 1-4 X' associated with a is defined by
a(u, v) = (Axu, v), Vu, v E X,
and (2.1) can be equivalently written as a system of linear algebraic equations
Axu = fx,
which we would like to solve by a preconditioned conjugate gradient method.
Here, a preconditioner is a mapping B : X' > X and we will look for
preconditioners such that (r, Br) is also symmetric and positive definite on X'.
In iteration k the method computes the residual
r = Axv, -fxe X',
and the preconditioner computes the increment to the approximate solution u^
as a linear combination of the preconditioned residual Br^ E X with precon-
ditioned residuals in earlier iterations.
11


It is well known that BA\ : X > X has only real positive eigenvalues, and
convergence of the preconditioned conjugate gradients method can be estab-
lished from the eigenvalues A of the preconditioned operator BAx; the condition
number
__ ^max(BAx)
K ~ Amin(BAx)'
gives a well-known bound on the error reduction, cf. e.g. [31],
where e^ u is the error of the solution in iteration k.
2.1 Additive Schwarz preconditioners
Because of their importance in the forthcoming development, we briefly
review the theory of abstract additive Schwarz preconditioners. Such a precon-
ditioner is specified by a decomposition of the solution space X into subspaces,
X = Xi + ... + XM, (2.2)
and by symmetric positive definite bilinear forms bi on X{. The preconditioner
is a linear operator
B-.X'^X, B-.r^u,
defined by solving the following variational problems on the subspaces and
adding the results:
M
B : r u = ^ uk Xk : bk(uk, yk) = (r, yk), Xk. (2.3)
fc=i
12


Theorem 2.1 (Dryja and Widlund [20]) If there exist constants C0, uj, and
a symmetric matrix E = (ey)jj=1 such that
M M
Vug X 3uk e Xk, k = 1,..., M : u = Uk, (2.4)
fc=i *:=1
Vk = l,...,MVukeXk:\\uk\\2a Vwfc e Xk, k 1,..., A/ a(ui, uf) ^ e^j ||U,|ja | Ma (2.6)
then
k = < CoUjp(E).
This theorem with proof can be found in Dryja and Widlund [20, Theo-
rem 1], or in the monograph by Toselli and Widlund [81, Theorem 2.7].
Remark 2.2 Because p(E) < ||£||00 and we can choose e^ = 1 if X^Xj ^ {0},
eij = 0 otherwise, we have the easy estimate
p(E)< max \{je{l,...,M}:XinXj^{0}}\ In the case when M = 1, the preconditioner simplifies to
r i> u e X : b(u, v) = (r, v), Vu X,
and Theorem 2.1 reduces to the statement that if
1 "uWl^Ml^CoWuWl Mu ex,
U)
then k < ujC0.
(2.8)
13


2.2 Abstract variational BDDC preconditioner
Next, we present the original BDDC preconditioner in a particularly simple
form originally formulated in [61] which has been inspired by a view of the
Neumann-Neumann methods, going back once more to Dryja and Widlund [20].
Suppose that the bilinear form a is defined and symmetric positive semidef-
inite on a larger space W D X. The abstract version of the original BDDC
preconditioner [17, 58] is characterized by a selection of an intermediate space V,
X C V C W. (2.9)
Assumption 2.3 The form a (, ) is positive definite on V.
Algorithm 2.4 (Abstract BDDC) Given the space V and a projection Q,
such that a (, ) is positive definite on V, and
X C V, Q : V X,
define the preconditioner B : r X' i u X by
B : r i> u = Qv, v G V : a(v,z) = (r, Qz), Vz G V. (2.10)
This formulation is remarkably simple, involving only the bilinear form a,
the space V, and the projection Q. We note that in typical substructuring
applications the spaces W and X are given, and the space V and the projection Q
are to be chosen.
14


At this point we would like to premise that the main goal of this dissertation
is to address several different ways of selecting V and Q. But before proceeding
any further, let us recall the abstract condition number bound of the original
two-level BDDC method, cf. also [61, Theorem 2].
Theorem 2.5 Let Q be a projection onto X, let us define a space VM by
VM = {v e V : Vz V : Qv = Qz = ||i;||^ < ||z||*} ,
and let u be a minimal constant such that
VveVM: \\Qv\\l <
UJ V
Then the preconditioner from Algorithm 2-4 satisfies
\\Qv\fc
K ||w||a
Proof. Define an operator G : X t VM by
(2.11)
G : u i-> v, -a(v,v) > min, s.t. v VM, u = Qv. (2.12)
Since a is positive definite on VM, the operator G is well defined. Define the
bilinear form 6 on A by b(u,v) = a(Gu,Gv). Now let u and v be as in (2.10).
Since v is the solution of
(v, v) (r,Qv) min, s.t. v VM,
it follows that u is the solution of
u X : b(u,x) (r,x), Vx X. (2-13)
15


It remains to compare ||u|| and ||u||^. Let u X and define v = Gu. Then,
from the minimization property (2.12) and the fact that Qu = u 6 X C VM, it
follows that
IMI* = IMI < IMI.
On the other hand,
IMIo = ll^ll^wlHI^WblHI*,
which concludes the proof.
16


3. Abstract Multispace BDDC
We will derive the abstract Multispace BDDC preconditioner from the ab-
stract additive Schwarz theory. However, unlike in this theory, we will decom-
pose some space between X and W rather than X. Suppose again that the
bilinear form a is defined and symmetric positive semidefinite on a larger space
W D X. In the design of the multispace preconditioner, we choose spaces 14,
k = 1,..., A/, such that
M
Xc^VkcW, (3.1)
fc=l
so this can also be viewed as replacing the space V in (2.9) by a sum ^4-
We do not wish to relate the two choices of the intermediate space at this point.
Assumption 3.1 The form a(-,-) is positive definite on each 14 separately.
Algorithm 3.2 (Abstract Multispace BDDC) Given spaces 14 and linear
operators Qk, k 1,..., M such that a (, ) is positive definite on each space 14,
and
M
X C ^ 14, Qk 14 i X,
k= 1
define the preconditioner B : r X' i> u G X by
M
= ^ Qkvk, vk£Vk\ a (vk, zk) = (r, Qkzk), Vzk G 14- (3.2)
k= 1
We first formulate the condition number bound in the full strength allowed
by the proof. We then find the bound used for the Multilevel BDDC as a
corollary.
17


Theorem 3.3 Define for all k 1,..., M the spaces Vf* by
Vk4 = {* G Vk : \/zk e Vk : Qkvk = Qkzk => ||ufc||^ < |j^fc||^} .
If there exist constants Co, u, and a symmetric matrix £ = (eij)(j=1, such that
M M
VuX 3vk eVk, k l,..., M : u = QkVk, £ INI < C0 |M|* (3.3)
fc=i
fc=i
VA = 1,...,A/ Vvfc e : ||Q*t;fc||= < w |
G Qk^ki k 1,..., A/ a (zj, ^ ejj ||2j||a \\zj||a ,
(3.4)
(3.5)
then the preconditioner from Algorithm 3.2 satisfies
Ax )
^min (DAx )
< C0up{£).
Proof. We interpret the Multispace BDDC preconditioner as an abstract ad-
ditive Schwarz method. So the essential idea of the proof is to map the assump-
tions of the abstract additive Schwarz estimate from the decomposition (2.2) of
the space X to the decomposition (3.1). Define the spaces
Xk = QkVk.
We will show that the preconditioner (3.2) satisfies (2.3), where bk is defined by
h(uk, Vk) = a (Gkx, Gkz), x, z e X, uk = QkGkx, yk = QkGkz (3.6)
with the operators Gk : X > Vjf4 defined by
1 M
Gk-.u^f vk, -a (vk, vk) -> min, s.t. vk Vf1, u = ^ Qkvk. (3.7)
fc=i
First, from the definition of operators Gk, spaces Xk, and because a is positive
definite on Vk by Assumption (3.1), it follows that Gkx and Gkz in (3.6) exist
18


and are unique, so bk is defined correctly. To prove (2.3), let vk be as in (3.2)
and note that vk is the solution of
(vk, vk) (r, Qkvk) -> min, vk G Vk.
Consequently, the preconditioner (3.2) is an abstract additive Schwarz method
and we only need to verify the inequalities (2.4)-(2.6). To prove (2.4), let u G X.
Then, with vk from the assumption (3.3) and with uk = QkGkvk as in (3.6), it
follows that
M M M
u = J2Uk' Y KHk* = Y ll^lio ^ C IMI
k=l fe=1 fc=l
Next, let uk G Xk. From the definitions of Xk and V^4, it follows that there
exist unique vk G such that uk = Qkvk. Using the assumption (3.4) and the
definition of bk in (3.6), we get
IMIa = WQkVkWl < W IMIa = W IMItfc -
which gives (2.5). Finally, (3.5) is the same as (2.6).
The next Corollary was given without proof in [63, Lemma 1]. This is the
special case of Theorem 3.3 that will actually be used. In the case when M = 1,
this result was proved in [61].
Corollary 3.4 Assume that the subspaces Vk are energy orthogonal, the opera-
tors Qk are projections, a (, ) is positive definite on each space Vk, and
Vu G X :
M
U
= ^ vk, vk G Vk
fc= l
M
U
^ ^ Qk^k
(3.8)
fc=i
Then the abstract Multispace BDDC preconditioner from Algorithm 3.2 satisfies
l|2
Amax{BAX) ^ \\QkVk\\a
K = ---)-z:r < U) = max SUp --s-2-
Amin^^x) k VfcgVfc ||^*:||a
(3.9)
19


Proof. We only need to verify the assumptions of Theorem 3.3. Let u X
and choose Vk as the energy orthogonal projections of u on 14- First, since the
spaces Vk are energy orthogonal, u = Y(,vk, Qk are projections, and from (3.8)
u = YlQkvk, we get that ||u||^ = )T) ||n/t||^ which proves (3.3) with Co = 1. Next,
the assumption (3.4) becomes the definition of the optimal uj in (3.9). Finally,
(3.5) with £ = I follows from the orthogonality of subspaces V*,.
Remark 3.5 The assumption (3.8) can be written as
M
'^QkPk =1,
k= 1 X
where Pk is the a-orthogonal projection from V) onto Vk. Hence, the prop-
erty (3.8) is a type of decomposition of unity. In the case when M 1, (3.8)
means that the projection Qj is onto X.
3.1 Finite element setting
Let fi be a bounded domain in Kd, d = 2 or 3, decomposed into N nonover-
lapping subdomains fi3, s = 1,..., N, which form a conforming triangulation of
the domain Q. Subdomains will be also called substructures. Each substructure
is a union of Lagrangian PI 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 T. The interface T
is a union of three different types of open sets: faces, edges, and vertices. The
substructure vertices will also be called corners. For the case of regular sub-
structures, such as cubes or tetrahedrons, we can use the standard geometric
definition of faces, edges, and vertices; cf., e.g., [42] for a more general definition.
20


Here, we find it more convenient to use the notation of abstract linear spaces
and linear operators between them instead of the space Mn and matrices. The
results can be easily converted to matrix language by choosing a finite element
basis. The space of the finite element functions on Q will be denoted as U. Let
Wa be the space of finite element functions on substructure Qs, such that all of
their degrees of freedom on dQs fl dil are zero. Let
W = W1 x x WN,
and consider a bilinear form a (, ) arising from the second-order elliptic problem
such as Poisson equation or a problem of linear elasticity.
Now U C W is the subspace of all functions from W that axe continuous
across the substructure interfaces. We are interested in the solution of the
problem (2.1) with X = U,
u U : a(u,v) = (f,v), Vv U, (3.10)
where the bilinear form a is associated on the space U with the system
operator A, defined by
A : U -> [/', a(u,v) = (Au,v) for all u,v E. U, (3-11)
and / U' is the right-hand side. Hence, (3.10) is equivalent to
Au = f. (3.12)
Define Uj C U as the subspace of functions that are zero on the interface
T, i.e., the interior functions. Denote by P the energy orthogonal projection
from W onto Uj,
P :w G W i Vi Uj : a (vj, Zi) = a (w, Zj), Vz/ £ f7/.
21


Functions from (I P) W, i.e., from the nullspace of P, are called discrete har-
monic; these functions are a-orthogonal to Uj and energy minimal with respect
to increments in Uj. Next, let W be the space of all discrete harmonic functions
that are continuous across substructure boundaries, that is
W = {I-P)U. (3.13)
In particular,
U = Uj W, UiaW. (3.14)
A common approach in substructuring is to reduce the problem to the inter-
face. Problem (3.10) is equivalent to two independent problems on the energy
orthogonal subspaces Uj and W, and the solution u satisfies u uj + u, where
Uj G Uj : a(uj,vj) = (f,vj), Vv/ G f//, (3.15)
ueW : a(u,v) = (/, v), Vv G W. (3.16)
The solution of the interior problem (3.15) decomposes into independent prob-
lems, one for each substructure. The reduced problem (3.16) is then solved
by preconditioned conjugate gradients. The reduced problem (3.16) is usually
written equivalently as
u G W : s(u, v) = {g, v), Vv G W,
where s is the form a restricted on the subspace W, and g is the reduced right-
hand side, i.e., the functional / restricted to the space W. The reduced right-
hand side g is usually written as
(9,v) = (f,v) ~ a{ui,v), WgVT, (3.17)
22


because a(ui,v) 0 by (3.14). 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 S, and (3.17) becomes the reduced right-hand
side. For details on the matrix formulation, see, e.g., [77, Sec. 4.6] or [81, Sec.
4.3],
The BDDC method is a two-level preconditioner characterized by the se-
lection of certain coarse degrees of freedom, such as values at the corners and
averages over edges or faces of substructures. Define W C 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 dVl, and WA C W as
the subspace of all functions such that their coarse degrees of freedom vanish.
Next, define Wn 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 Wn are uniquely determined by the values of
their coarse degrees of freedom, and
WA a W, and W = WA VTn (3.18)
We assume that
a is positive definite on W. (3.19)
This will be the case when a is positive definite on the space U, for which
problem (2.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
23


from {//, that is,
£// C WA. (3.20)
In other words, the coarse degrees of freedom depend on the values on substruc-
ture boundaries only. From (3.18) and (3.20), it follows that the functions in
Wfi are discrete harmonic, that is,
Wn = (I-P)Wn. (3.21)
Next, let E be a projection from W onto U, defined by taking some weighted
average on substructure interfaces. That is, we assume that
E :W U, EU U, E2 = E. (3.22)
Since a projection is the identity on its range, it follows that E does not change
the interior degrees of freedom,
EUi = Ui, (3.23)
since Uj C U. Finally, we show that the operator (I P) E is a projection.
From (3.23) 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
[(/ P) E}2 = (I P) E (I P) E
= (I P) (E P) E (3.24)
= (I P) (I P) E = (I P) E.
Remark 3.6 In [59, 61], the whole analysis was done in spaces of discrete
harmonic functions after eliminating Ui, and the space W was the solution space.
24


In particular, W consisted of discrete harmonic functions only, while the same
space here would be (/ P)W. In our context, the decomposition of this space
used in [59, 61] would be written as
(/ P)W = (/ P)WA Wn, (/ P)WA a Wn. (3.25)
In the next section, the space X will be either U or W.
3.2 Two-level BDDC as Multispace BDDC
We show several different ways the original two-level BDDC algorithm can
be interpreted as multispace BDDC. We first consider BDDC applied to the
reduced problem (3.16), that is, (2.1) with X = W. This was the formulation
considered in [59]. Define the space of discrete harmonic functions with coarse
degrees of freedom coinciding across the interface
Wr = (I P) W.
Because we work in the space of discrete harmonic functions and the output of
the averaging operator E is not discrete harmonic, denote
Er = (/ P) E. (3.26)
In an implementation, discrete harmonic functions are represented by the
values of their degrees of freedom on substructure interfaces, cf., e.g. [81]; hence,
the definition (3.26) serves formal purposes only, so that everything can be
written in terms of discrete harmonic functions without passing to the matrix
formulation.
25


Algorithm 3.7 ([61], BDDC on the reduced problem) Define the precon-
ditioner r G W' i> ueW by
u = EpWp, wr G Wp : a (wp, Zp) = (r, £rzr), Vzp G Wr. (3.27)
Proposition 3.8 ([61]) The BDDC preconditioner on the reduced problem in
Algorithm 3.7 is the abstract Multispace BDDC from Algorithm 3.2 with M = 1
and the space and operator given by
X = W, Vi = Wr, Qx = Ep. (3.28)
Also, the assumptions of Corollary 3.4 are satisfied.
Proof. We only need to note that the bilinear form a(, ) is positive definite
on Wp C W by (3.19), and the operator Ep defined by (3.26) is a projection
by (3.24). The projection Ep is onto W because E is onto U by (3.22), and
I P maps U onto W by the definition of W in (3.13).
Using the decomposition (3.25), we can split the solution in the space Wp
into the independent solution of two subproblems: mutually independent prob-
lems on substructures as the solution in the space WpA (/ P)W&, and a
solution of the global coarse problem in the space Wu- The space Wr has the
same decomposition as in (3.25), i.e.,
Wr = WrAWn, and WpAaWn, (3.29)
and Algorithm 3.7 can be rewritten as follows.
26


Algorithm 3.9 ([59], BDDC on the reduced problem) Define the precon-
ditioner r W' i> u e W by u = Er (tcrA + wn), where
wrA £ WrA : a (^rA, ^a) = (u ^r-^rA) > ^zrA WpA, (3.30)
wn lUn ' (%, ^rn) = (u ErZrn), Vzrn £ Wm (3.31)
Proposition 3.10 The BDDC preconditioner on the reduced problem in Algo-
rithm 3.9 is the abstract Multispace BDDC from Algorithm 3.2 with M = 2 and
the spaces and operators given by
X = W, V^WrA, V2 = Wn, Ql=Q2 = Er. (3.32)
Also, the assumptions of Corollary 3-4 are satisfied.
Proof. Let r W'. Define the vectors Vi, i 1,2 in Multispace BDDC
by (3.2) with Vt and Qi given by (3.32). Let u, WrA, be the quantities in
Algorithm 3.9, defined by (3.30)-(3.31). Using the decomposition (3.29), any
wr Wr can be written uniquely as wr = wpA+wn for some WrA and tt>n
corresponding to (3.2) as V\ = Wr a and v2 = w n, and u = Er (icrA + wn).
To verify the assumptions of Corollary 3.4, note that the decomposition
(3.29) is aorthogonal, a(-, ) is positive definite on both UVa and Vbn as sub-
spaces of Wr by (3.19), and Er is a projection by (3.24).
Next, we present a BDDC formulation on the space U with explicit treat-
ment of interior functions in the space f// as in [17, 58], i.e., in the way the
BDDC algorithm was originally formulated.
27


Algorithm 3.11 ([17, 58], original BDDC) Define the preconditioner r G
U' i> u G U as follows. Compute the interior pre-correction
uj G U] : a (uj, Zj) = (r, z/), Vz/ G f//. (3.33)
Set up the updated residual
rBeU\ (rB,v) = (r,v) a(ui,v), Vv G U. (3.34)
Compute the substructure correction
uA = EwA, wa G Wa : a (wA, zA) = (rB, EzA), VzA WA. (3.35)
Compute the coarse correction
un = Ewn, wn G Wn : a (wn, zu) = (rB, Ezn), Vzn G Wn- (3.36)
Add the corrections
uB = uA + un.
Compute the interior post-correction
Vi G Ui : a(vi,Zj) = a(uB,zj), Vz/G f//. (3.37)
Apply the combined corrections
u = uB vj 4- /. (3.38)
The interior corrections (3.33) and (3.37) decompose into independent
Dirichlet problems, one for each substructure. The substructure correction (3.35)
decomposes into independent constrained Neumann problems, one for each sub-
structure. Thus, the evaluation of the preconditioner requires three problems to
28


be solved in each substructure, plus solution of the coarse problem (3.36). In
addition, the substructure corrections can be solved in parallel with the coarse
problem.
Remark 3.12 As it is well known [17], the first interior correction (3.33) can be
omitted in the implementation by starting the iterations from an initial solution
such that the residual in the interior of the substructures is zero,
a (u, Zj) (fx,Zi) = 0, 'izj C//,
i.e., such that the error is discrete harmonic. Then the output of the precondi-
tioner is discrete harmonic and thus the errors in all the CG iterations (which are
linear combinations of the original error and outputs from the preconditioner)
are also discrete harmonic by induction.
The following proposition will be the starting point for the multilevel case.
Proposition 3.13 The original BDDC preconditioner in Algorithm 3.11 is the
abstract Multispace BDDC from Algorithm 3.2 with M 3 and the spaces and
operators given by
X = u, Vi = [/;, V2 = (I- P)WA, K3 = Wn, (3.39)
Qi = I, Q2 = Q3 = (I-P)E, (3.40)
and the assumptions of Corollary 3.f are satisfied.
Proof. Let r E U'. Define the vectors vt, i = 1,2,3, in Multispace BDDC
by (3.2) with the spaces V) given by (3.39) and with the operators Qi given
29


by (3.40). Let uj, rB, w&, wB, uB, V[, and u be the quantities in Algorithm 3.11,
defined by (3.33)-(3.38).
First, with V) = t//, the definition of V\ in (3.2) with k = 1 is identical to
the definition of uj in (3.33), so Uj = iq.
Next, consider wA G WA defined in (3.35). We show that wA satisfies (3.2)
with k = 2, i.e., v2 wA. So, let 2a WA be arbitrary. From (3.35) and (3.34),
a (wA, za) = (rB, Eza) = (r, EzA) a (u/, EzA). (3.41)
Now from the definition of U; by (3.33) and the fact that PEzA Uj, we get
(r, PEza) a (u/, PEzA) = 0, (3-42)
and subtracting (3.42) from (3.41) gives
a (wA, 2a) = (r, (I P) Eza) a (uj, (I P) EzA)
= (r,(I~P) Eza),
because a (u/, (/ P) EzA) = 0 by orthogonality. To verify (3.2), it is enough
to show that PwA = 0; then wA (/ P)WA = V2. Since P is an a-orthogonal
projection, it holds that
a (PwA, PwA) = a (wA, PwA) = (rB, EPwA) 0, (3.43)
where we have used EUi C Ui following the assumption (3.23) and the equality
(rB, zi) = (r, zj) a (uh 2/) = 0,
for any zj G Uj, which follows from (3.34) and (3.33). Since a is positive definite
on W D U[ by assumption (3.19), it follows from (3.43) that PwA = 0.
30


In exactly the same way, from (3.36) (3.38), we get that if wn 6 Wn is
defined by (3.36), then u3 = w;n satisfies (3.2) with k = 3. (The proof that
Pwn = 0 can be simplified but there is nothing wrong with proceeding exactly
as for wa-)
Finally, from (3.37), V[ = P (Ewa + Ewn), so
u = U/ + (uB Vi)
= uj + (/ P) Ewa + {I P) Ewn
= Q\V\ + Q2V2 + Q3V3.
It remains to verify the assumptions of Corollary 3.4.
First, the spaces Wn and Wa are a-orthogonal by (3.18) and, from (3.20),
(/ P) Wa C Wa,
thus (I P) Wa -La Wn- Clearly, (I P) Wa -La Uj. Since Wn consists of dis-
crete harmonic functions from (3.21), so Wn -La Uj, it follows that the spaces Vi,
i = 1,2,3, given by (3.39), axe a-orthogonal.
Next, (/ P) E is a projection by (3.24), and so are the operators Qi
from (3.40).
It remains to prove the decomposition of unity (3.8). Let
u' = Uj + Wa + wn G U, ut G Uj, wA (I P) WA, wn Wn, (3.44)
and let
v = ui + (I P) Ewa + (I P) Ewn-
31


From (3.44), w& + wn £ U since v! £ U and ii/ 6 (// C U. Then E (w& + u>n) =
wa + wn by (3.22), so
v uj + (I P) E (wa + wn)
= U! + (/ P) (wA + inn)
= uj + wa + wu ~ u',
because both wa and wx\ are discrete harmonic.
The next theorem shows an equivalence of the three algorithms introduced
above.
Theorem 3.14 The eigenvalues of the preconditioned operators from Algo-
rithm 3.7, and Algorithm 3.9 are exactly the same. They are also the same
as the eigenvalues from Algorithm 3.11, except possibly for multiplicity of the
eigenvalue equal to one.
Proof. From the decomposition (3.29), we can write any w £ ITp uniquely as
w = wa+w\\ for some u>a £ WVa and u-n £ VFn, so the preconditioned operators
from Algorithms 3.7 and 3.9 are spectrally equivalent and we need only to show
their spectral equivalence to the preconditioned operator from Algorithm 3.11.
First, we note that the operator A : U t-> U' defined by (3.11), and given in the
block form as
An Ait
A =
Ati ^4rr
32


with blocks
A-ii Uj > Uj, Air '-Uj ¥ W',
An : W -> U'n Apr : W -> W1,
is block diagonal and Ari = Air = 0 for any u G U, written as u = uj + w,
because U] _La W. Next, we note that the block .App : IT' > W is the Schur
complement operator corresponding to the form s. Finally, since the block
An is used only in the preprocessing step, the preconditioned operator from
Algorithms 3.7 and 3.9 is simply MrrArr t W' > u W.
Let us now turn to Algorithm 3.11. Let the residual r U be written as
r = ri + rr, where 77 U', and rp W'. Taking rp = 0, we get r rj, and it
follows that rs ub = vj = 0, so u = uj. On the other hand, taking r = rp
gives uj = 0, rB = rr, Vj Pub and finally u = (/ P) E(wA + ^n), so u G W.
This shows that the off-diagonal blocks of the preconditioner M are zero, and
therefore it is block diagonal
0 A/pp
Next, let us take u = ui, and consider rp = 0. The algorithm returns rg = ub =
vj = 0, and finally u = uj. This means that MnAnUj = U], so Mu = Ajj.
The operator A : U {/', and the block preconditioned operator MA : r £
U' > u G U from Algorithm 3.11 can be written, respectively, as
An 0 , MA = I 0
0 .App 0 MrrArr
33


where the right lower block MrrArr r E IT' -> u G W is exactly the same as
the preconditioned operator from Algorithms 3.7 and 3.9.
The BDDC condition number estimate is well known from [58], Following
Theorem 3.14 and Corollary 3.4, we only need to estimate ||(/ P) Ew\\a on W.
Theorem 3.15 ([58])
satisfies k < lj, where
The condition number of the original BDDC algorithm
UJ = sup
wW
\\(I-P)Ew\\2a
(3.45)
Remark 3.16 In [58], the theorem was formulated by taking the supremum over
the space of discrete harmonic functions (I P)W. However, the supremum
remains the same by taking the larger space W D (/ P)W, since
\\{I-P)Ew\\l ^ \\{I-P)E(I-P)w\\l
IMIa " \\(I-P)M\l
from E (/ P) = E, which follows from (3.23), and from ||u;||a > ||(7 F)u;||a,
which follows from the a-orthogonality of the projection P.
Before proceeding into the Multilevel BDDC section, let us write concisely
the spaces and operators involved in the two-level preconditioner as
U, £ U £ WA Wn = W C W.
We axe now ready to extend this decomposition into the multilevel case.
34


H

Figure 3.1: An example of a domain decomposition for the two-level (top)
and the three-level (bottom) BDDC methods.
35


4. Multilevel BDDC
In this chapter, we generalize the two-level BDDC preconditioner to multiple
levels, using the abstract Multispace BDDC framework from Algorithm 3.2. The
substructuring components from Section 3.2 will be denoted by an additional
subscript 1? as D, s = 1 etc., and called level 1. The level 1 coarse
problem (3.36) will be called the level 2 problem. It has the same finite element
structure as the original problem (2.1) on level 1, so we put U? = Wm. Level 1
substructures are level 2 elements and level 1 coarse degrees of freedom axe level 2
degrees of freedom. Repeating this process recursively, level i 1 substructures
become level i elements, and the level i substructures are agglomerates of level i
elements. Level i substructures are denoted by Q;s, s = l,Ni} and they are
assumed to form a conforming triangulation with a characteristic substructure
size Hi. For convenience, we denote by Qq the original finite elements and put
H0 = h. The interface Tt on level i is defined as the union of all level i boundary
nodes, i.e., nodes shared by at least two level i substructures, and we note
that Tj C r,_i. Level i 1 coarse degrees of freedom become level i degrees
of freedom. The shape functions on level i are determined by minimization
of energy with respect to level i 1 shape functions, subject to the value of
exactly one level i degree of freedom being one and the other level i degrees of
freedom being zero. The minimization is done on each level i element (level i 1
substructure) separately, so the values of level i 1 degrees of freedom are in
general discontinuous between level i 1 substructures, and only the values of
36


level i degrees of freedom between neighboring level i elements coincide. For an
example of a decomposition for two and a three-level method, see Figure 3.2.
The development of the spaces on level i now parallels the finite element
setting in Section 3.1. Denote Ui = Wni_i. Let W* be the space of functions on
the substructure fl, such that all of their degrees of freedom on dfl- n dCl are
zero, and let
Wi = Wl x x WtNi.
Then Ui C Wt is the subspace of all functions from W that are continuous
across the interfaces Tj. Define Un C Ui as the subspace of functions that are
zero on Tj, i.e., the functions interior to the level i substructures. Denote by
Pi the energy orthogonal projection from Wi onto UIt,
Pi : Wi G Wi i vH eUn:a (vu, zh) = a (wu zH), V27i G UH.
Functions from (I Pi) Wi, i.e., from the nullspace of Pt, are called discrete
harmonic on level i; these functions are a-orthogonal to Un and energy minimal
with respect to increments in Un. Denote by Wi C Ui the subspace of discrete
harmonic functions on level i, that is
Wi = (I-Pi)Ui. (4.1)
In particular, Un _La Wx. Define W{ C W, as the subspace of all functions such
that the values of any coarse degrees of freedom on level i have a common value
for all relevant level i substructures and vanish on n <9Q, and WAi C Wt as
the subspace of all functions such that their level i coarse degrees of freedom
vanish. Define Wni as the subspace of all functions such that their level i coarse
37


degrees of freedom between adjacent substructures coincide, and such that their
energy is minimal. Clearly, functions in Wm are uniquely determined by the
values of their level i coarse degrees of freedom, and
w* -La Wuu Wi = W We assume that the level i coarse degrees of freedom are zero on all functions
from Uji, that is,
Un C WAi. (4.3)
In other words, level i coarse degrees of freedom depend on the values on level i
substructure boundaries only. From (4.2) and (4.3), it follows that the functions
in W\\t are discrete harmonic on level i, that is
Wni = (I-Pi)Wili. (4.4)
Let E be a projection from IT) onto Ui, defined by taking some weighted average
on Ti
Ef.Wi^Ui, EiUIi = UIi, Ef = Ei.
Since projection is the identity on its range, Ei does not change the level i
interior degrees of freedom, in particular
EiUn = UH. (4.5)
The hierarchy of spaces and operators is shown concisely in Figure 4.1.
The Multilevel BDDC method is defined recursively [17, 63, 64] by solving
the coarse problem on level i only approximately, by one application of the
preconditioner on level i + 1. Eventually, at the top level L 1, the coarse
38


Figure 4.1: Spaces, embeddings and projections in the Multilevel BDDC.
39


problem, which is the level L problem, is solved exactly. We need a more formal
description of the method here, which is provided by the following algorithm.
Algorithm 4.1 (Multilevel BDDC) Define the preconditioner r\ E U[ i>
U\ E U\ as follows:
for i = 1,..., L 1,
Compute interior pre-correction on level i,
ua E Uh : a (uH, zH) = (ru zH), \/zH E Un. (4.6)
Get an updated residual on level i,
rm Ui, (rm, v{) = (n, - a (-uH, v,), Vuj E Ui. (4.7)
Find the substructure correction on level i:
WAi WAi : a {wAi, zAi) = (rBi, EiZAi), VzAi E WAi. (4.8)
Formulate the coarse problem on level i,
Wni : a (wm, zUi) = (rBi, EiZni), Vzni Wni, (4.9)
If i = L 1, solve the coarse problem directly and set ul = mnL-1,
otherwise set up the right-hand side for level i + 1,
ri+1 G Wm, (n+i,Zi+i) = (rBi,EiZi+1), Vzi+i G Wm = Ui+U (4.10)
end.
40


for i = L
Average the approximate corrections on substructure interfaces on level i,
uBi = Ei(wAi + ui+i). (4.11)
Compute the interior post-correction on level i,
vn Un : a (vu, zti) = a (uBi, zH), Vzn G Un. (4.12)
Apply the combined corrections,
Ui = un + uBi vn. (4-13)
end.
We can now show that the Multilevel BDDC can be cast as the Multi-
space BDDC on energy orthogonal spaces, using the hierarchy of spaces from
Figure 4.1.
Lemma 4.2 The Multilevel BDDC preconditioner in Algorithm 4-1 is the ab-
stract Multispace BDDC preconditioner from Algorithm 3.2 with M = 2L 1,
and the spaces and operators
X = UU Vi = Un, V2 = (I ~ Pi)WA1, V3 = UI2,
V, = (/ P2)WA2, Vb = UI3, ... (4.14)
F2L-4 = {1 Pl-2)WAC-2, V2L-3 U]L-1,
y2L-2 {I Pl-i)WaL-1, y2L-l = WnL-l,
41


Qi I, Q2 Q3 (I Pi) Ei,
Q4 = Qs = (/ Pi) Ex (/ P2) E2,
(4.15)
Q2L-4 Q2L-3 (I Pi) Ei (/ P1-2) El-2,
Qil-2 = Q2L-1 = (I Pi) Ex (/ - Pl-i) EL-i,
and the assumptions of Corollary 3.4 are satisfied.
Proof. Let r\ U[. Define the vectors Vk, k = 1,..., 2L 1 by (3.2) with
the spaces and operators given by (4.14)-(4.15), and let un, rgt, wm, ri+j,
Ubi, vH, and ux be the quantities in Algorithm 4.1, defined by (4.6)-(4.13).
First, with V\ Un, the definition of v\ in (3.2) is (4.6) with i = 1 and
un V\ . We show that in general, for level i 1,..., L1, and space k 2i l,
we get (3.2) with 14 = Un, so that vl = un and in particular v2l-3 = Uil-i-
So, let zn £ Uji, i = 2,..., L 1, be arbitrary. From (4.6) using (4.10) and
(4-7),
a(un, zn) = (ri, zH) = (rBi-1, Ei-^zn) = (4.16)
= (rj_i, Ei-iZn) aiun-i^i-izn)
Since from (4.6), using the fact that Pi-iEi-iZn Un-1, it follows that
(rxi,PiiEiiZn) a(%in\,Pi\Ei\Zji) 0,
we get from (4.16),
(^/u Zn) (rj_i, (/ Pi1) EiiZji) a (Ujii, (/ Pi-i)Ei-iZn),
and because a (u{i-1, (I Pi-i)Ei-iZn) = 0 by orthogonality, we obtain
a{un, zn) = (r<_i, (/ Pt-1) Ei-XzH)
42


Repeating this process recursively using (4.16), we finally get
a(uIi,zIi) = (ri-i, (/ Pi~i) Ei-iZu) = ...
= (n, (I Pi)Ex (I /*_!) E^h) .
Next, consider w&i £ Wa* defined by (4.8). We show that for = 1,..., L1,
and k = 2i, we get (3.2) with 14 = WAu so that vk wAi, and in particular
V2L-2 = u>al-i- So, let ZAi £ WAi be arbitrary. From (4.8) using (4.7),
a (wau zAi) = (rBi, EiZAi) = (n, EtZAi) a (uIU EiZAi) (4.17)
From the definition of Un by (4.6) and since P{ExzAi G Un it follows that
(Tj, PiExZAi) U (u/i, EiExZAi) 0,
so (4.17) gives
0, (WAi: ZAi) (^"ij (I Pi) EiZAi) 0, (u/j, (/ Px) ExZAi)
Next, because a (un, (I Pi) Exz Ai) = 0 by orthogonality, and using (4.10),
a (u>A, ZAi) = iXii (I Pi) ExZAi) if Bih Eil (/ Pi) ExZAi)
Repeating this process recursively, we finally get
a(wAi,ZAi) = (ri, (I ~ Pi) EiZAi) = ...
= (r1? (/ Pi) Ei (I Pi) EiZ^).
To verify (3.2), it remains to show that PiWAi = 0; then uiAl & (I Pi)WAi = V*.
Since Pi is an a-orthogonal projection, it holds that
a {PiWAi, PwAi) = a (wAi, PiWAi) = (rBi, EiPiWAi) = 0,
43


where we have used E-JJa C Uu following the assumption (4.5) and the equality
(j"Bit Zli) (^"i) Zli) ^ Zji) 0
for any zIt e Uu, which follows from (4.6) and (4.7).
In exactly the same way, we get that if w\\l-\ WUL-i is defined by (4.9),
then i>2l-i = whl-i satisfies (3.2) with k = 2L 1.
Finally, from (4.11 )-(4.13) for any i = L 2,..., 1, we get
Ui U/j 4 U[3i U/j
U/j 4- (/ P^) Ei {w&i + tq+1)
= tlji 4- (/ Pi) Ei 4- UIi+i + (I Pi+i) Ei+x (WAi+l + Ui+2)]
= Uii +
+ (I Pi) Ei [uAi 4-... 4- {I Pl- 1) El-i {wAl-i + tdiL-i)] >
and, in particular for U\,
U\ = un +
+ (I Pi) Ei [u>ai 4-... 4- (/ Pl- 1) El- 1 (wAl-i + ^ni,1)]
= QlVl 4- Q2V2 4- ... 4- Q2L-2V2L-2 + Q2L-\V2L-\-
It remains to verify the assumptions of Corollary 3.4.
The spaces Wm and WAi, for alii = 1,..., L 1, are a-orthogonal by (4.2)
and from (4.3),
(i Pi)wAi cwAi,
thus (I Pi) WAi is a-orthogonal to W\n. Since Wm = Ui+i consists of dis-
crete harmonic functions on level i from (4.4), and Un+1 C Ui+i, it follows by
induction that the spaces 14, given by (4.14), are a-orthogonal.
44


We now show that the operators Qk defined by (4.15) are projections. From
our definitions, coarse degrees of freedom on substructuring level i (from which
we construct the level i + 1 problem) depend only on the values of degrees of
freedom on the interface Tj and F^ C F, for j > z. Then,
(I PJEjV PJEiil Pj)Ej = (I- PJEiV Pj)Ej. (4.18)
Using (4.18), and since (/ Pi) is a projection by (3.24), we get
[(/ -Pi)Ei---{I- P^ Eif = (/ Pi) Ei (/ Pi) Ei (/ Pi) Ei
= (I-Pi)Ei---(I- Pi) Eu
so the operators Qk from (4.15) are projections.
It remains to prove the decomposition of unity (3.8). Let it* G Ui, such that
u'f = uh + wAi + iti+i, (4-19)
uh e UH, wAi G (/ Pi) WAl, ui+i Ui+i (4.20)
and
Vi = un + (I Pi) EiWAi + (I Pi) EiUi+i. (4.21)
From (4.19), wAi + ul+\ £/, since Ui Ui and uIt G Un C Ui. Then
Ei [wAi + i+i] = wAi + ui+i by (4.5), so
Vi = uji + (I Pi)Ei [wAi + ui+1] = uji + (/ - Pi) [wAi + Uj+i] =
= uh + wAi + ui+i = uh + wAi + ui+i = a',
because wAi and ui+j are discrete harmonic on level i. The fact that ui+l in
(4.19) and (4.21) are the same on arbitrary level i can be proved in exactly the
45


same way using induction and putting ui+1 in (4.19) as
Ui+1 u/i+i + ... + Wal-i + Wnc-i,
un+1 (//i+i, wal-\ (I Pc,-i) £ WnL-i)
and in (4.21) as
Ui+\ un+1 + ... + (/ Pj+i) Ej+i (/ Pl_i) Pl_i (?ual-i + f^nL-i),
which concludes the proof.
The following bound follows from writing of the Multilevel BDDC as Mul-
tispace BDDC in Lemma 4.2 and the estimate for Multispace BDDC in Corol-
lary 3.4.
Lemma 4.3 If for some ui > 1,
||(/ PJE^Wl < uj IK^ V^ai (/ A) WAU
||(/ Pi)ElUl2\\l < u \\uI2\\2a Vui2 UI2,
(4.22)
II(/-Pi)Pi {I ~ Pl-x)EL-XwUL-X\\1 then the Multilevel BDDC preconditioner (Algorithm 4-1) satisfies k < uj.
Proof. Choose the spaces and operators as in (4.14)-(4.15) so that Un = V\ G
Vi = Un, wAi v2 V2 = (I Pi) WA\, u^nt-i = ^2L-i V2L-1 = Wnt-i*
Tlie bound now follows from Corollary 3.4.
The following lemma is an immediate consequence of Lemma 4.3, and it can
be viewed as a multilevel analogy of Theorem 3.15. In fact, in the same way as
46


Theorem 3.15, formulated as Theorem 5.2, will serve as a starting point for the
adaptive selection of constraints for the two-level BDDC method, the following
Lemma, formulated as Theorem 6.1, will serve as a starting point for adaptive
selection of constraints for the Multilevel BDDC method.
Lemma 4.4 If for some cu* > 1,
W(I -Pi^mWl^LJiWwiWl, vWiewu i = \,...,L-\, (4.23)
then the Multilevel BDDC preconditioner (Algorithm 4-1) satisfies k < cut.
Proof. Note from
Wi, and generally (I
Lemma 4.3 that (I Pi) WAi C Hai C Wi, Uj2 C Wui C
- Pi) c WAi C Wi, UIi+i C Wm C Wi. m
47


5. Adaptive Coarse Degrees of Freedom
We formulate an algorithm for adaptive selection of the coarse degrees of
freedom for the two-level BDDC method. It was presented in [61] starting
from corner constraints only, formulated in terms of FETI-DP, with the result
translated to BDDC. Later the method has been extended in [65, 79] to the case
of a general space W and implemented in BDDC directly using a projection on
the subspace W. The current formulation allows for an explicit coarse space and
hence for a multilevel extension.
The space W is constructed using coarse degrees of freedom. These can
be, e.g., values at corners, and averages over edges or faces. The space W is
then given by the requirement that the coarse degrees of freedom on adjacent
substructures coincide; for this reason, the terms coarse degrees of freedom and
constraints axe used interchangeably. The edge (or face) averages are necessary
in 3D problems to obtain scalability with subdomain size. Ideally, one can prove
the polylogarithmic condition number bound
k < const ^1 + log , (5.1)
where H is the subdomain size and h is the finite element size.
Remark 5.1 The initial selection of constraints in the proposed adaptive ap-
proach will be done such that (5.1) is satisfied. See, e.g., [43] for theoretical
reasoning for these constraints.
48


To choose the space W, cf. [61, Section 2.3], suppose we are given a linear
operator C : W > X and define,
W = {w G W : C (I E) w = 0} .
(5.2)
The values Cw will be called local coarse degrees of freedom, and the space W
consists of all functions w whose local coarse degrees of freedom on adjacent
substructures have zero jumps. To represent their common values, i.e., the
global coarse degrees of freedom of vectors u G W, suppose there is a space Uc
and linear operators Qj, : U > Uc, Rc : Uc > X such that Rc is one-to-one,
and injection R:U^W such that
The space W then satisfies
W = {w G W : 3uc G Uc: Cw = Rcuc},
and from (5.3), for w G W, the unique uc that satisfies Cw = Rcuc is given by
In order to formulate the adaptive algorithm, we first restate the condition
number bound from (2.11) in a way suitable for our purposes, cf. Theorem 3.15.
Theorem 5.2 ([61, Theorem 3]) The condition number bound (3.45) of the
two-level BDDC satisfies k < u, where
CR = RcQp.
(5.3)
uc = QpV, w = Rv.
U) = sup
wW
\\(I-P)Ew\\l = gup ||(7 (/ P)EH1*
(5.4)
49


With respect to Remark 3.16, we can conveniently look for the condition
number bound w only in the subspace Wr = (/ P)W. Next, observe that
(/ E) Pv 0 for all v W, so we can define the space W in (5.2) using
discrete harmonic functions w (/ P) W, for which
(I (I P) E)w = (I P) (I E) w, (5.5)
because Pw = 0 if w G (I P) W. Clearly, the bound (5.4) can be found as a
maximum eigenvalue of an associated eigenvalue problem, using (5.5) written as
((/ P) (I E) w, (/ P) (I E) z)a = A (w, z)a Vz Wr. (5.6)
Remark 5.3 The eigenvalue problem (5.6) corresponds to the right-hand side
of (5-4) combined with (5.5). This is motivated by the definition (5.2) of the
space W: in the adaptive algorithm we will prescribe certain weighted jumps
of functions to be zero across substructure interfaces.
The following is a well known result from linear algebra, cf., e.g., [16, The-
orem 5.2].
Lemma 5.4 (Courant-Fisher-Weyl minmax principle) Letc{-,-) be sym-
metric positive semidefinite bilinear form on vector space V of dimension n
and b(-,-) symmetric positive definite bilinear form on V. Then the generalized
eigenvalue problem
w G V : c (w, u) = Ab (w, u) Vu V,
has n linearly independent eigenvectors W( and the corresponding eigenvalues are
real and nonnegative and the eigenvectors are stationary points of the Rayleigh
50


quotient c(w,w) /b(w,w), with the stationary values equal to A*.
^2 > > An > 0. Then, for any subspace Vk C V of dimension n
c (w, w)
max T7------T
wevklw^o b (w, w)
^ Xk+1,
Order Aj >
-k,
with equality if
Vk = {w V : c(we, w) = 0, W = 1,..., k} .
Since the bilinear form on the left-hand side of (5.6) is symmetric positive
semidefinite and the bilinear form on the right-hand side is symmetric positive
definite, Lemma 5.4, using (5.5), implies
Corollary 5.5 ([61]) The generalized eigenvalue problem (5.6) has eigenvalues
Ai > A2 > ... > A > 0. Denote the corresponding eigenvectors by W(. Then,
for any k 1,..., n 1, and any linear functionals Le on Wr, t = 1,..., k,
'\\(I-P)(I-E)w\\l _
max
II^L
: w £ Wr, Le (w) = 0 = 1,..., k > > Xk+i,
with equality if
Le (w) a((I P) (/ E) we, (/ P) (/ E) w). (5.7)
Therefore, because (/ E) is a projection, the optimal decrease of the con-
dition number bound (5.4) can be achieved by adding to the constraint matrix
C in the definition of W the rows ct defined by cjw = Le (w). However, solving
the global eigenvalue problem (5.6) is expensive, and the vectors q are not of the
form required for substructuring, i.e., each q with nonzero entries corresponding
to only one corner, an edge or a face at a time.
51


For these reasons, we replace (5.6) by a collection of local problems, each
defined by considering only two adjacent subdomains Qs and Q1. All quantities
associated with such pairs will be denoted by the superscript st. Here, subdo-
mains are considered adjacent if they share an edge in 2D, or a face in 3D. We
note that edges in 2D will be regarded as faces. Using also (5.5), the generalized
eigenvalue problem (5.7) becomes a problem to find w E Wp such that
ast pstj (J Est) _ pst) (7 pet) = Xaet ^ ^ ^ £ yfst
(5.8)
Assumption 5.6 The comer constraints are already sufficient to prevent rela-
tive rigid body motions of any pair of adjacent substructures, so
\/w E Wst : Astw = 0 = (/ Est) w = 0,
i.e., the comer degrees of freedom are sufficient to constrain the rigid body modes
of the two substmctures into a single set of rigid body modes, which are contin-
uous across the interface Tst.
The maximal eigenvalue ust of (5.8) is finite due to Assumption 5.6, and we
define the heuristic condition number indicator
u = max {u/ : and 0* are adjacent} . (5.9)
Considering two adjacent subdomains Qs and Ql only, we get the added
constraints Lp (w) = 0 from (5.7) as
ast ((I Pst) (I Est) we, (I Pst) (I Est) w) = 0 W=1 (5.10)
where wp are the eigenvectors corresponding to the k largest eigenvalues from (5.8)
52


To formulate a numerical algorithm, we need to write the generalized eigen-
value problem (5.8) and the added constraints (5.10) in terms of matrices and
vectors. To avoid complicated notation, we now drop the superscripts st, or,
equivalently, let us consider a domain which consists of only two substruc-
tures. For convenience, we will also identify finite element functions with vectors
formed by their degrees of freedom, and we will also identify linear operators
with their matrices, in bases that will be clear from the context.
The vectors of the local substructure degrees of freedom ws Ws and the
vector of the global degrees of freedom u G U are related by ws = Rsu, where
Rs is the restriction operator (a zero-one matrix), so that
RS:U-+WS, RsRsT = I, (5.11)
and R : U W. The Schur complement matrices Sa are assumed to be
symmetric and positive semidefinite. Let us consider the vectors and matrices
to be given in a block form
ws Ss Rs
w wl , 5 = Sl , R = Rl
(5.12)
We will need a more specific construction of the matrix C in the substruc-
turing framework. We build a block diagonal matrix C satisfying (5.3) by
C =
cl
Cs = RscQTpRsT.
(5.13)
Then (5.3) follows from (5.13) and (5.11).
Here is an interpretation. The matrix Cs maps a vector of local degrees of
freedom on substructure i to a vector of local coarse degrees of freedom on the
53


substructure, and R* restricts a vector of all global coarse degree of freedom to
a vector of local coarse degree of freedom on substructure s. A global coarse
degree of freedom is given by a row of Qp. The operator Qj, acts on vectors of
global degrees of freedom in U and it selects global coarse degrees of freedom
in Uc as linear combinations of global degrees of freedom. In our problems,
there are corner coarse degrees of freedom, which are values at corners, and
edge (face) coarse degrees of freedom, which are linear combinations of values
on edges (faces).
Consider the space W given some initial constraint matrix C containing at
least corner constraints. Let us denote D C (7 E) and define the orthogonal
projection onto null D by
U = I -Dt (DDt)~1 D.
The generalized eigenvalue problem (5.6) now becomes
n (/ p)T (i e)t s(i-e)(i p) nw = Ansiiw. (5.14)
Since
nulinsn C null n (7 pf (I Ef S (7 E) (7 P) II, (5.15)
the eigenvalue problem (5.14) reduces in the factorspace modulo null IISII to the
problem with the operator on the right-hand side positive definite. In some of
the computations, we have used the subspace iteration method LOBPCG [45] to
find the dominant eigenvalues and their eigenvectors. The LOBPCG iterations
then simply run in the factorspace. To use standard eigenvalue solvers, (5.14)
is converted to a matrix eigenvalue problem by penalizing the components in
54


null D and rigid body modes, already described in [61, 79], and recalled here for
completeness. We can formulate 5.14 as a generalized eigenvalue problem with
the matrix on the right-hand side positive definite using the following procedure.
Theorem 5.7 ([61]) Let a > 0. Then the nonzero eigenvalues and the eigen-
vectors of (5.14) are same as those of the generalized eigenvalue problem
n (/ P)T (/ E)t S(I E){I P)Uw = A (USU + a (/ II)) w, (5.16)
The matrix on the left-hand side is symmetric positive semidefinite and if the
pair of substructures intersects boundary with Dirichlet boundary conditions, the
matrix on the right-hand side is symmetric positive definite.
In practice, we choose a to be roughly the same magnitude as S. Note
that if the eigenvalues are computed approximately, the result will in general
depend on a. Also, for subspace iteration methods the matrices, in particular
the Schur complement S, need not to be formulated explicitly, and only matrix-
vector products are evaluated. However, the matrices on both sides may be still
singular because of rigid body modes that move substructures s and t as a whole.
To reduce (5.16) to an eigenvalue problem with the matrix on the right-hand
side positive definite, we use matrix Z that generates a superspace of rigid body
modes of the two substructures
null S C range Z.
The matrix Z can be available from finite element software or can be easily
computed from the geometry of the finite element mesh with rigid body modes
as its columns, e.g., [81, Chapter 8]. To avoid using any information other
55


than the system matrices, we can instead use as Z a block diagonal matrix of
the coarse basis functions of the two subtructures because their span contains
the rigid body modes. However, the computations in this case will be more
expensive because there are typically more coarse basis functions for the two
substructures than the number of the rigid body modes.
Remark 5.8 We have used in the computations the true rigid body modes com-
puted from the mesh geometry for the two-level method, and the coarse basis
functions on higher levels of Multilevel BDDC.
We find a basis of null (nsn + a (I n))by computing the nullspace of a
much smaller symmetric positive semi-definite matrix,
null (ZT (nsn + a (I n)) Z) = range K,
and applying the QR decomposition
ZK = QR, QtQ = I,
which gives
range Q null (nsn + a (/ n)).
Consequently,
n = / qqt,
is the orthogonal projection onto range (nSTl + a (I n)).
The following theorem follows similarly as Theorem 5.7.
56


Theorem 5.9 ([61]) The nonzero eigenvalues and the corresponding eigenvec-
tors of
n (/ Ef s(i -e) n* = a* (nsn + <*(/- n)) wk,
are the same as nonzero eigenvalues and corresponding eigenvectors of
Xwk = XkYwk, (5.17)
where
X = U{I Ef S(I E)Yl,
v = (TT (nsn + a (/ n)) IT + a (/ n)).
In addition, Y is symmetric positive definite.
Remark 5.10 The automatic decomposition can result, due to the vagaries of a
mesh partitioner, in very irregular substructures including spurious mechanisms,
see Figs. 5.2, or 7.8. It such cases the nullspace of the Schur complement is in
general unknown. See Section 5.1 for further discussion.
From the matrix form (5.15) of the eigenvalue problem, the constraints to
be added are
Le H = wjn (/ pf (/ Ef s(i -e)(i- p) n; = o.
That is, we wish to add to the constraint matrix C the rows
ce = wJU (I P)T (/ Ef S(I -E)(I P) n. (5.18)
Proposition 5.11 The vectors q, constructed for a domain consisting of only
two substructures fls and have matching entries on the interface between the
two substructures, with opposite signs.
57


Proof: Consider the vector w W that has two entries equal to 1,
corresponding to a degree of freedom on the interface, and all other entries
equal to 0. Using the definition of q and because (/ E) u = 0 for all u U,
we get Cgw L( (w) = 0.
It remains to construct the augmentation of the primal constraint matrix QP
from the augmentation q. Due to Proposition 5.11, each row of q can be split
into two blocks and written as
Q
The augmentation of Qp is then constructed by simply taking c£, and computing
RsTcf
(5.19)
Because Ra is a 0 1 matrix, it means that columns q are formed by a scattering
of the entries in c)T. Each column of q defines a coarse degree of freedom, which
is used to augment the matrix QP as
[Qpq]. (5.20)
Unfortunately, the added columns will generally have nonzero entries over
all of the interface of 0s and Ql, including the edges in 3D where ils and Ul
intersect other substructures. Consequently, the added q are not of the form
required for substructuring, i.e., each q with nonzeros in one edge or face only.
In the computations reported in Section 7, we drop the adaptively generated edge
constraints in 3D. Then it is no longer guaranteed that the condition number
indicator uj < t. However, the method is still observed to perform well.
58


Now Corollary 5.5 and the formulation of the constraints (5.18)-(5.20) sug-
gest a way to decrease the indicator u and the proposed adaptive algorithm
follows.
Algorithm 5.12 (Adaptive BDDC [61]) Find the smallest k for every two
adjacent substructures to guarantee that Xk+i < t, where t is a given tolerance,
and add the constraints (5.10) to the definition ofW.
5.1 Preconditioned LOBPCG
The most important step towards future parallel implementation of the
adaptive method seems to be an efficient solution of local generalized eigenvalue
problems. An attractive approach is to use an inverse-free method, such as the
method by Golub and Ye [32], or LOBPCG by Knyazev [45]. These methods al-
low the matrices to be in a matrix-free format, i.e., as functions for matrix-vector
multiplication, which are readily available in our implementation. In particular,
LOBPCG might be more suitable because it allows for the resolution of more
eigenpairs at once and it can simply run in the factorspace with the operator on
the right-hand side only positive semi-definite. Initial experiments reveal that
the non-preconditioned LOBPCG as well as the method of Golub and Ye work
well for reasonably hard problems [79]. Unfortunately, it turns out that many
iterations are required for problems with extremely irregular structures and/or
high jumps in coefficients, and preconditioning of the local eigenproblems seems
to be necessary [65].
One desirable property of a preconditioner is that it must be invariant on
the null 5. Also, unless the component of the solution in the direction of the
59


nullspace is small, the errors will accumulate, which may eventually result in
instability of the code at the Rayleigh-Ritz step, but only after a large number
of steps [46]. Because the Schur complement operator S plays a central role in
the operators on both sides of the local generalized eigenvalue problems (5.14),
(5.16) , and (5.17), a straightforward idea is to use a local version of the BDDC
preconditioner, denoted here as Mloc. The only difference is that this precondi-
tioner acts on the larger space W, unlike the preconditioners in Section 3.2, and
therefore in place of averaging E we have used an injection R instead. The coarse
space correction is obtained using corner (and in 3D edge) constraints, shared
by the two substructures. We restrict the action of the preconditioner in the
suitable subspaces using the two projections II resp. II introduced previously,
milocu resp.. nmflocim. (5.21)
Unfortunately, as we have pointed out already in Remark 5.10, the auto-
matic decomposition of the finite element mesh can result, due to the vagaries
of a mesh partitioner (in the experiments we have used METIS 4.0 [34]), in very
irregular substructures including spurious mechanisms, see Figs. 5.2, or 7.8. In
such cases the nullspace of the Schur complement is in general unknown, and the
LOBPCG iterations with the preconditioner (5.21) can fail. To detect the eigen-
vectors in the nullspace of the operator on the right-hand side in (5.17), resp.
(5.16) , we have used again LOBPCG with Mloc as a preconditioner. However
because the Schur complement is singular, so is Mloc. To circumvent this, we
have applied a shift and the action of the preconditioner was in this case given
as Mloc + I. Once the nullspace is detected, we can enrich the nullspace basis,
reconstruct the projection n and rerun the preconditioned LOBPCG iterations.
60


Figure 5.1: Comparison of the non-preconditioned (top) vs. preconditioned
LOBPCG (bottom) for one of the faces with large jumps in coefficients of the
composite cube problem, cf. Chapter 7 and Fig. 7.2. Estimated eigenvalue
errors are in the panels on the left-hand side, and Euclidean norms of residuals
for different eigenpairs are in the panels on the right-hand side. We see that
LOBPCG without a preconditioner showed essentially no convergence, and with
the preconditioner we have reached convergence in less than 30 iterations.
61


Figure 5.2: A pair of substructures of the mining reel problem from Figure 7.8,
obtained from the automatic decomposition by METIS 4.0. We see that one of
the substructures has 4 spurious rigid body modes.
62


Iteration number
50 100 150 200 250 300 350
Iteration number
Figure 5.3: Comparison of the noil-preconditioned (top) vs. preconditioned
LOBPCG (bottom) for the detection of spurious rigid-body modes of the pair of
subdomains from Fig. 5.2. Estimated eigenvalue errors are in the panels on the
left-hand side, and Euclidean norms of residuals for different eigenpairs are in the
panels on the right-hand side. We see that LOBPCG without a preconditioner
essentially did not detect the nullspace, and the application of preconditioner
led to relatively fast detection of the nullspace.
63


Iteration number Iteration number
Figure 5.4: Convergence of the preconditioned LOBPCG for the pair of subdo-
mains from Fig. 5.2 with the nullspace detection (Fig. 5.3). Estimated eigenvalue
errors are in the panels on the left-hand side, and Euclidean norms of residuals
for different eigenpairs are in the panels on the right-hand side.
64


6. Adaptive Multilevel BDDC
We build on the previous two chapters to propose a new variant of the
Multilevel BDDC with adaptive selection of constraints on each level. The
starting point is Lemma 4.4, formulated as a multilevel analogy to Theorem 5.2.
Theorem 6.1 The condition number bound k < uu of Multilevel BDDC from
Algorithm 4-1 satisfies
uj = nf-fui, (6.i)
where
= SUp _ sup !!(/-(/-**<
wWi ll^illa weWi ll^illa
The development of adaptive selection of constraints in Multilevel BDDC
now proceeds similarly as in Chapter 5. We formulate (6.1) as a set of eigen-
value problems for each decomposition level. On each level we solve for every
two adjacent substructures a generalized eigenvalue problem and we add the
constraints to the definitions of W*.
The heuristic condition number indicator is defined as
ui = nfjj1 max {u>tst : D? and Sl\ are adjacent} . (6.2)
We now describe the Adaptive Multilevel BDDC in more detail. The al-
gorithm consists of two main steps: (i) setup (adaptive selection of constraints),
and (ii) loop of the preconditioned conjugate gradients with the Multilevel
BDDC from Algorithm 4.1 as a preconditioner. The setup can be summarized
as follows (cf. [78, Algorithm 4] for the 2D case):
65


Algorithm 6.2 (Setup of Adaptive Multilevel BDDC) Adding of coarse
degrees of freedom to guarantee that the condition number indicator w < rL_1,
for a given a target value r:
for levels i = 1 : L 1,
Create substructures with roughly the same numbers of degrees of freedom (one
can use a graph partitioner, e.g., METIS 4-0 [34])-
Find a set of initial constraints (in particular sufficient number of comers),
and set up the BDDC structures for the adaptive algorithm (the next loop
over faces).
for all faces Ei on level i,
Compute the largest local eigenvalues and corresponding eigenvectors, until
the first mst is found such that A^st < r.
Compute the constraint weights and update the global coarse degrees of
freedom selection matrix.
end.
Setup the BDDC structures for level i and check size of the coarse problem:
if small enough, call this level L problem, factor it directly, break the loop.
end.
66


6.1 Implementation remarks
The matrices of the averaging operator E were constructed with entries pro-
portional to the diagonal entries of the substructure matrices before elimination
of interiors, which is also known as a stiffness scaling [39].
6.1.1 Initial constraints
Following Remark 5.1, in order to satisfy the polylogarithmic condition num-
ber bounds, we have used corners, and in 3D corners with arithmetic averages
over edges as initial constraints. It is essential (Assumption 5.6) to generate a
sufficient number of corners as initial constraints to prevent rigid body motions
between any pair of adjacent substructures. This topic has been addressed in
the literature several times cf., e.g., [9, 50] in a different context, or a recent con-
tribution in the context of BDDC by Burda et al. [10]. The selection of corners
in our implementation follows the original implementation by Dohrmann [17].
Let Afst denote the set of all nodes shared by substructures 0s and Ql. The first
corner cf £ Afst is chosen as a node shared by the largest number of substruc-
tures. The second corner cf £ Afst is chosen as a node with greatest distance
from cf. For problems in three dimensions, a third corner cf £ Afst is chosen as
a node for which the area of the triangle connecting cf, cf, and cf is maximized
from vector cross product. However if all nodes in Afst are shared only by the
two substructures, the algorithm starts by a random selection of an initial node
in Afst and cf is identified as a node maximizing the distance from the initial
node, as suggested in [10].
67


6.1.2 Algebraic coarse elements
The substructures in engineering applications were obtained using a graph
partitioner METIS 4.0 [34]. The connectivity graph has been weighted in both
vertices and edges in order to minimize the number of cuts. The vertex weights
were given by the total number of degrees of freedom in the substructure and
the weights of graph edges were determined as the numbers of the degrees of
freedom identified on faces by the adaptive algorithm.
The substructures on higher levels were then treated in an algebraic way,
unlike the geometric substructures illustrated in Figure 3.2, as (coarse) elements
with energy minimal basis functions. The routines for multilevel algorithm must
allow for (coarse) elements with variable number of nodes, and they must also
allow for variable number of degrees of freedom on each node (corresponding to
a face) the number of nodes and the number of their degrees of freedom is
a-priori unknown due to their adaptive selection. For this purpose we transform,
renumber and reorder all globs so that the nodes corresponding to corners have
the lowest numbers, followed by nodes corresponding to edges, and finally by
nodes corresponding to faces identified by the adaptive algorithm.
It is also convenient to use the same assembling routines for higher levels. To
this end, the globs that do not consist of single nodes (edges or faces) axe replaced
by virtual nodes with coordinates given by arithmetic averages of coordinates
of all nodes belonging to that particular glob. We note that the coordinates of
these nodes are in fact not used by the algorithm, so the only purpose is to allow
the use of the basic two-level routines without any modifications.
68


Finally, we remark that instead of interior pre- and post-correction on the
lowest decomposition level, cf. equations (4.6)-(4.7) and (4.12)-(4.13), we reduce
the problem to interfaces in the pre-processing step, cf. also Remark 3.12.
6.1.3 Adaptive constraints
The adaptive algorithm uses matrices and operators that are readily avail-
able in an implementation of the BDDC method with an explicit coarse space,
with one exception: in order to satisfy the local partition of unity, cf. [62, eq. (9)],
EfRf = I,
we need to generate locally the weight matrices Ef.
In the computations reported in Chapter 7, we drop the adaptively generated
edge constraints in 3D. Then, it is no longer guaranteed that the condition
number indicator ui < rL_1. However, the method is still observed to perform
well. Since the constraint weights are thus supported only on faces, and the
entries corresponding to edges are set to be zero, we orthogonalize and normalize
the vectors of constraint weights to preserve numerical stability.
69


7. Numerical Examples
The main purpose of this chapter is to compare performance of the standard
two-level BDDC with the adaptive and multilevel extensions. For consistency
with our previous research [61, 65, 79], we first present results for several aca-
demic examples. Results of Multilevel BDDC for scalar problems can be found in
[64], The computations were done in Matlab (version 7.8.0.347 (R2009a)). The
generalized eigenproblems on pairs of substructures were solved by setting up
the matrices and using standard methods for the symmetric eigenvalue problem
in Matlab, and we have also tested LOBPCG by Knyazev [45] with a precon-
ditioner described in Section 5.1. The two-dimensional results are reproduced
from [78], the three-dimensional results appear here for the first time.
7.1 Two-dimensional results
The method has been tested on planar elasticity (with A = 1, p = 2) on
a square domain discretized by Lagrange bilinear finite elements with 1 182 722
degrees of freedom. The domain was decomposed into 48 x 48(= 2304) subdo-
mains on the second level and into 3 x 3(= 9) subdomains on the third-level.
Such a decomposition leads to the coarsening ratio //*///*_! = 16, for i 1,2.
In order to test the adaptive selection of constraints, one single edge is jagged
on both levels, see Fig. 7.1. Recall that edges in 2D are regarded as faces.
In the first set of experiments, we have compared performance of the non-
adaptive BDDC method with 2 and 3 decomposition levels. The results are
presented in Tables 7.1 and 7.2. As suggested by Lemma 4.4, the convergence
of the algorithm deteriorates when additional levels are introduced.
70


In the next set of experiments, we have tested the adaptive algorithm for
the two-level BDDC. The results are summarized in Table 7.5. The algorithm
performs consistently with our previous formulation in [61]. The eigenvalues as-
sociated with edges between substructures clearly distinguish between the single
problematic edge and the others (Table 7.3). Adding the coarse dofs created
from the associated eigenvectors decreases the value of the condition number
indicator u and improves convergence at the cost of increasing the number of
coarse dofs.
Finally, we have tested the performance of the Adaptive Multilevel BDDC
for the model problem with three-level decomposition (Fig. 7.1). Because the
number of coarse degrees of freedom depends on a-priori chosen value of r and
the coarse basis functions on level i become shape basis functions on level i + 1,
the solutions of local eigenvalue problems will depend on r as well. This fact
is illustrated by Table 7.4 for r = 2, and r = 10 (the local eigenvalues for
t = 3 were essentially same as for r = 2). Comparing the values in the two
panels of this table, we see that lower values of r result in worse conditioning of
the local eigenvalue problems on higher decomposition level. This immediately
raises the conjecture that it might not be desirable to decrease the values of
r arbitrarily low in order to achieve a better convergence of the method. On
the other hand, for the model problem, comparing the convergence results for
the two-level method (Table 7.5) with the three-level method (Table 7.6), we
see that with the adaptive constraints we were able to achieve nearly the same
convergence properties of both methods.
71


90
92
96
100
94
2300
2251
2204
-
193 194 195 196
145 146 147 148
97 98 99 100 101
49 50 51 52 53
1 2 3 4 5
0 2 4 6 6 10
2301
2253
2205
2157
2109
2302
2254
2206
2158
2110
2303
2255
2207
2159
2111
2304
2256
2208
2160
2112
96
96
94
99
90
100|
90 -
80-
70-
60-
50-
40-
30
20 -
10
0
'UirUUUUlr
20
40 60 80 1 00
Figure 7.1: Two remote corners of the two-level decomposition into 48 x 48 (=
2304) subdomains (top), and the decomposition into 9 subdomains for the three-
level method (bottom). The jagged edge from the lower decomposition level
(top) is indicated on the second-level decomposition (bottom) by the thick line.
72


Table 7.1: Results for the planar elasticity from Fig. 7.1 obtained using non-
adaptive 2-level method. Constraints are corners, or corners and arithmetic
averages over edges, denoted as c, c+f, resp., and Nc is number of constraints
(coarse degrees of freedom), C is size of the coarse problem related to size of a
subdomain problem, k is the approximate condition number computed from the
Lanczos sequence in conjugate gradients, and it is the number of iterations for
relative residual tolerance 10-8.
constraint Nc C K it
c 4794 9.3 18.41 43
c+f 13818 26.9 18.43 32
Table 7.2: Results for the planar elasticity from Fig. 7.1 obtained using non-
adaptive 3-level method. Nc is the number of coarse degrees of freedom on
the first (+ the second) decomposition level, C is the relative coarsening with
respect to the size of substructures on the first level (the size of the coarse
problem for the three-level method is negligible). Other headings are the same
as in Table 7.1.
constraint Nc C K it
c 4794 + 24 1.0 67.5 74
c+f 13818 + 48 3.0 97.7 70
73


Table 7.3: The largest eigenvalues Xst several pairs of subdomains s and t of the 2-level elasticity problem from Fig. 7.1
(top), with (s,t) = (2,50) being the jagged edge.
s t At,i Kt, 2 Kt, 3 ^st, 4 Xst, 5 Kt, 7 Kt, 8
1 2 3.8 2.4 1.4 1.3 1.2 l.i 1.1 1.1
1 49 6.0 3.5 2.7 1.4 1.3 l.i 1.1 1.1
2 3 5.4 2.6 1.6 1.3 1.2 l.i 1.1 1.1
2 50 24.3 18.4 18.3 16.7 16.7 14.7 13.5 13.1
3 4 3.4 2.4 1.4 1.3 1.1 1.1 1.1 1.1
3 51 7.4 4.6 3.7 1.7 1.4 1.3 1.2 1.1
49 50 12.6 5.1 4.3 1.9 1.6 1.3 1.2 1.2
50 51 8.7 4.8 3.9 1.8 1.5 1.3 1.2 1.2
50 98 7.5 4.6 3.7 1.7 1.4 1.3 1.2 1.1
74


Table 7.4: The largest eigenvalues X3t,k of the local problems for several pairs
of subdomains s, t on the level i = 2, of. Fig. 7.1 (lower panel), with r = 2 (top)
and with r = 10 (bottom). The jagged edge is between subdomains 2 and 5.
s t ^s,l Kt, 2 Kt, 3 ^st, 4 Xst, 5 Kt, 7 Kt, 8
1 2 16.5 9.0 5.4 2.6 2.1 1.4 1.3 1.3
1 4 6.5 4.7 1.9 1.7 1.3 1.2 1.2 1.1
2 3 23.1 9.4 4.6 3.2 2.1 1.6 1.4 1.3
2 5 84.3 61.4 61.4 55.9 55.8 49.3 48.0 46.9
3 6 13.7 8.8 4.4 2.2 1.9 1.4 1.3 1.2
4 7 6.5 4.7 1.9 1.7 1.3 1.2 1.2 1.1
5 6 18.9 13.1 11.3 3.8 2.6 2.1 1.9 1.5
5 8 17.3 12.9 10.8 3.6 2.3 2.0 1.8 1.4
8 9 13.7 8.8 4.4 2.2 1.9 1.4 1.3 1.2
1 2 7.7 4.5 2.7 1.6 1.4 1.2 1.2 1.1
1 4 3.6 3.0 1.5 1.5 1.2 1.2 1.1 1.1
2 3 10.9 4.8 2.7 1.7 1.5 1.2 1.2 1.1
2 5 23.2 17.2 13.7 13.7 12.7 12.4 11.0 10.9
3 6 6.1 4.2 2.5 1.5 1.3 1.2 1.1 1.1
4 7 3.6 3.0 1.5 1.5 1.2 1.2 1.1 1.1
5 6 9.8 6.2 4.1 2.1 1.6 1.5 1.3 1.2
5 8 8.6 5.9 3.9 2.0 1.5 1.4 1.2 1.2
8 9 6.1 4.2 2.5 1.5 1.3 1.2 1.1 1.1
75


Table 7.5: Results for the planar elasticity from Fig. 7.1 obtained using the
adaptive 2-level method. Here, r is condition number target, lj is condition
number indicator, and the other headings are the same as in Table 7.1.
T Nc C U) K it
oo(=c) 4794 9.3 - 18.41 43
10 4805 9.4 8.67 8.34 34
3 18110 35.3 2.67 2.44 15
2 18305 35.7 1.97 1.97 13
Table 7.6: Results for the planar elasticity from Fig. 7.1 obtained using the
adaptive 3-level method. Headings are the same as in Tables 7.2 and 7.5.
r Nc C U) K it
oo(=c) 4794 + 24 1.0 - 67.5 74
10 4805 + 34 1.0 > (9.80)2 37.42 60
3 18110 + 93 3.9 > (2.95)2 3.11 19
2 18305 + 117 4.0 > (1.97)2 2.28 15
76


7.2 Three-dimensional results
The performance of the Adaptive BDDC in the presence of jumps in material
coefficients has been tested on a cube with material parameters E = 106 Pa and
v 0.45 (rubber), penetrated by four bars with parameters E 2.1* 1011 Pa and
v = 0.3 (steel), consisting of 107811 degrees of freedom, and distributed into
8 substructures with 30 corners, 16 edges and 15 faces, see Fig. 7.2, see also [65].
Comparing the results in Tables 7.8 and 7.9 we see that with r = 10 000 only 10
additional averages over faces decrease the number of iterations nearly 3 times,
and with r = 2 the number of iterations decreased more than 13 times compared
to the noil-adaptive algorithm with arithmetic averages over all globs (c+e+f)
whereas the number of constraints increased approximately 2.5 times. Results
(and numbers of nonzeros in the action) of the BDDC preconditioner in Table 7.9
can be compared with results obtained by incomplete Cholesky factorization
applied to the global matrix in Table 7.7. We see that for a lower number of
iterations the fill-in of the Cholesky factor was quite high when compared with
the fill-in of the subdomain and the coarse problems in the BDDC method.
The performance of the Adaptive Multilevel BDDC in the presence of
jumps in material coefficients has been tested on a cube designed similarly as the
problem above (also with the same material parameters), this time consisting of
823875 degrees of freedom and distributed into 512 substructures, 721 corners,
1176 edges and 1 344 faces on the first decomposition level, and 4 substructures,
6 corners, one edge and 4 faces on the second decomposition level, see Fig. 7.3.
77


Figure 7.2: Finite element discretization and substructuring of the cube with
jumps in coefficients, consisting of 107811 degrees of freedom, distributed into
8 substructures with 30 corners, 16 edges and 15 faces (the bars cut the sub-
structures only through faces). Notably, similar problems are solved in practice
to determine numerically (anisotropic) properties of composite materials [72],
Courtesy of Jakub Sfstek.
Table 7.7: Results for the cube from Fig. 7.2 obtained using a precondition-
ing by incomplete Cholesky factorization. The global stiffness matrix has size
107811 with 7 737641 nonzeros. Here, nnz(R) is the number of nonzeros in
the upper triangular Cholesky factor R, fill-in is the relative fill-in computed
as 2 times nnz(R) divided by the number of nonzeros in the global stiffness
matrix, k is the approximate condition number computed from the Lanczos se-
quence in conjugate gradients, it is number of iterations for relative residual
tolerance 10~8
With the drop tol. zero-level of: nnz(R) ill-in th fill-in e method c cond id not iter
no fill-in 3922 726 1.01 - oo
1 103 9 784 734 2.53 4.14 106 331
1 io-4 30968 534 8.00 2.25 106 170
1 lO"5 88 125 845 22.78 119.12 37
1 10~6 194448 707 50.26 3.63 15
1 10"7 273649916 70.73 1.88 10
78


Table 7.8: Results for the cube from Fig. 7.2 obtained using the non-adaptive
2-level method. Constraints are corners, or corners and arithmetic averages over
edges and faces denoted as c, c+e, c+e+f resp., and c+e+f (3eigv), corresponding
to corner constraints, arithmetic averages, and 3 weighted averages over each face
obtained using the adaptive method. Next, Nc is the number of constraints, C
is the size of the coarse problem related to size of a subdomain problem, k
is the approximate condition number computed from the Lanczos sequence in
conjugate gradients, it is number of iterations for relative residual tolerance 10~8.
constraint Nc K it
c 90 408114 455
c+e 138 125 378 307
c+e+f 183 18915.1 211
c+e+f (3eigv) 183 1267.61 81
Table 7.9: Results for the cube from Fig. 7.2 obtained using the adaptive 2-
level method. Here, r is the condition number target, tJ is the condition number
indicator. An approximate number of nonzeros of the Cholesky factor of a
substructure problem is 8 500 000 for all values of r, and the number of nonzeros
in the Cholesky factor of the coarse problem is denoted by nnz(c). The other
headings are the same as in Table 7.8.
T Nc nnz(c) CJ K it
oo(=c+e) 138 6618 268 390.00 125378.00 307
10000 148 7402 5096.07 1 843.70 104
1000 159 8 271 368.78 173.57 38
100 162 8448 5.94 6.42 24
5 198 13029 4.99 4.55 21
2 465 87012 < 2 2.80 16
79


Comparing the results in Tables 7.10 and 7.11 we see that similar to the pre-
vious problem, a relatively small number of (additional) constraints leads to
a dramatic decrease in number of iterations in the 2-level method. When the
non-adaptive 2-level is replaced by the 3-level method, Tables 7.10 and 7.12, the
condition number estimate as well as the number of iterations grows. However,
with the adaptive 3-level approach (Table 7.13) we were able to achieve nearly
the same convergence properties as in the adaptive 2-level method (Table 7.11).
Table 7.10: Results for the cube from Fig. 7.3 obtained using the non-adaptive
2-level method. The headings are the same as in Table 7.8.
constraint Nc K it
c 2163 312 371 > 3 000
c+e 5 691 45849 1521
e+e+f 9 723 16384 916
c+e+f (3eigv) 9 723 3848 367
Table 7.11: Results for the cube from Fig. 7.3 obtained using the adaptive
2-level method. The headings are the same as in Table 7.9.
T Nc C u K it
oo(=c+e) 5691 3.54 o(104) 45 848.60 1521
10000 5883 3.66 8 776.50 5 098.60 441
1000 6027 3.75 5.33 9.92 32
10 6149 3.82 6.25 6.66 28
5 9119 5.67 < 5 4.79 24
2 25009 15.54 < 2 2.92 18
80


Figure 7.3: Finite element discretization and substructuring of the large cube
with jumps in coefficients, consisting of 823 875 degrees of freedom, distributed
into 512 substructures with 721 corners, 1176 edges and 1 344 faces on the first
decomposition level (top), and 4 substructures, 6 corners, one edge and 4 faces
on the second decomposition level (bottom). Courtesy of Jakub Sfstek.
81


Table 7.12: Results for the cube from Fig. 7.3 obtained using the non-adaptive
3-level method. The headings are the same as in Tables 7.2 and 7.8.
constraint Nc C K it
c 2 163 + 18 0.34 + 0.01 o(107) 3000
c+e 5691 + 21 0.88 + 0.01 o(106) > 3000
c+e+f 9 723 + 33 1.51 + 0.02 461 750 1573
c+e+f (3eigv) 9 723 + 33 1.51 + 0.02 125 305 981
Table 7.13: Results for the cube from Fig. 7.3 obtained using the adaptive
3-level method. The headings are the same as in Tables 7.6 and 7.9.
r Nc C u K it
oo(=c+e) 5691 + 21 0.88 + 0.01 - o(106) > 3000
10000 5883 + 28 0.91 + 0.02 8 776.50 26874.40 812
1000 6027 + 34 0.94 + 0.02 766.82 1 449.50 145
100 6027 + 53 0.94 + 0.03 99.05 100.89 59
10 6 149 + 65 0.96 + 0.04 7.93 7.91 30
5 9119 + 67 1.42 + 0.04 < 5 6.18 25
2 25009+ 122 3.89 + 0.08 < 2 3.08 18
82


7.2.1 Applications to engineering problems
7.2.1.1 Application of Multilevel BDDC to a dam
The performance of the Multilevel BDDC has been tested on the realistic
engineering problem of a dam discretized using 3 800 080 tetrahedral finite ele-
ments with 668916 nodes and 2 006 748 degrees of freedom, with two variants
of substructuring: first decomposed into 400 substructures with 3990 corners,
3 070 edges and 2 274 faces, see Fig. 7.4, and then decomposed into 1 024 sub-
structures with 10693 corners, 7713 edges and 6 182 faces, see Fig. 7.6.
The results with non-adaptive constraints for the decomposition into 400
substructures are summarized in Tables 7.14 and 7.15 for the two- and three-
level methods, respectively. Results with non-adaptive constraints for the de-
composition into 1024 substructures are summarized in Tables 7.16 and 7.17 for
the two- and three-level methods, respectively. At the first glance, comparing
the values in Tables 7.14 and 7.15, it might appear that for the decomposition
into 400 substructures, the increase in the number of iterations is not too signif-
icant if one uses corners and arithmetic averages over edges (or faces). However,
comparing the values in Tables 7.16 and 7.17, it turns out that for the decom-
position into 1 024 substructures the number of iterations needed for the 3-level
method double, or even triple, compared to the 2-level method.
Nevertheless, we see that for this problem the simple arithmetic averages
already work well enough as there are no interfaces that require extra work -
the quality of the decomposition is uniform, as seen in Figures 7.4-7.7.
83