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 FETIDP 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 twolevel 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, twolevel, 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 twolevel 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 HebertonMorimoto.
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 proofreading 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 CNS0325314, CNS0719641, and DMS0713876.
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 Twolevel 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 Twodimensional results............................................ 70
7.2 Threedimensional results.......................................... 77
8. Conclusion.......................................................... 98
References............................................................. 100
vii
FIGURES
Figure
3.1 An example of a domain decomposition for the twolevel (top) and
the threelevel (bottom) BDDC methods............................. 35
4.1 Spaces, embeddings and projections in the Multilevel BDDC. ... 39
5.1 Comparison of the nonpreconditioned (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 lefthand side,
and Euclidean norms of residuals for different eigenpairs are in the
panels on the righthand 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 nonpreconditioned (top) vs. preconditioned
LOBPCG (bottom) for the detection of spurious rigidbody modes
of the pair of subdomains from Fig. 5.2. Estimated eigenvalue er
rors are in the panels on the lefthand side, and Euclidean norms of
residuals for different eigenpairs are in the panels on the righthand
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 lefthand side, and
Euclidean norms of residuals for different eigenpairs are in the panels
on the righthand side.......................................... 64
7.1 Two remote corners of the twolevel decomposition into 48 x 48 (=
2304) subdomains (top), and the decomposition into 9 subdomains
for the threelevel method (bottom). The jagged edge from the lower
decomposition level (top) is indicated on the secondlevel 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 2level 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 108............................................. 73
7.2 Results for the planar elasticity from Fig. 7.1 obtained using non
adaptive 3level 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 threelevel 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 2level 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 2level 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 3level 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, fillin is the
relative fillin 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 108. With
the zerolevel of fillin the method did not converge............... 78
xiv
7.8 Results for the cube from Fig. 7.2 obtained using the nonadaptive
2level 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 108......................................... 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 nonadaptive
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 2level
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 nonadaptive
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 3level
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 nonadaptive 2level 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 nonadaptive 3level 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 nonadaptive 2level 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 nonadaptive 3level 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 2level 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 3level 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 nonadaptive 2level 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 2level 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 3level 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 2level 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 3level 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 matrixvector multiplication, and preconditioner
matrixvector multiplication. There is often a preprocessing step consisting of
a transformation of variables before an iterative method is applied; in that case,
subroutines to preprocess the righthand 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 socalled 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 FETIDP.
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 nonoverlapping 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
matrixvector 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 matrixvector 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 NeumannNeumann 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 NeumannNeumann 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
matrixvector multiplications and the substructure nullspace vectors. For ex
ample, in application of BDD to mixed discretization of a flow problem [12],
the matrixvector 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 righthand sides orthogonal to the
nullspace gives the original Finite Element Tearing and Interconnecting (FETI)
method by Farhat and Roux [26], later called FETI1. 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 (FETIDP) 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 FETIDP 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 matrixvector 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 FETIDP 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 socalled
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 FETIDP by Cros [13] and by Fragakis and Pa
padrakakis [29]. The BDDC and, equivalently, FETIDP 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
FETIDP, 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, twolevel, BDDC has been extended into threelevels 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 threelevel 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, FETIDP 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
longstanding 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 FETIDP 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 FETIDP 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 eigenvaluebased
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 34) 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 AdaptiveMultilevel 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 righthand 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 ua = [a (u,u)\1^2 is called the energy
norm. The operator Ax : X 14 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 wellknown 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
NeumannNeumann 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
twolevel 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 24 satisfies
\\Qv\fc
K
wa
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. (213)
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 2ja \\zja ,
(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
l2
Amax{BAX) ^ \\QkVk\\a
K = )z:r < U) = max SUp s2
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 aorthogonal 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 secondorder 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, (311)
and / U' is the righthand 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 aorthogonal 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 = {IP)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 righthand
side. For details on the matrix formulation, see, e.g., [77, Sec. 4.6] or [81, Sec.
4.3],
The BDDC method is a twolevel 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 = (IP)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 Twolevel BDDC as Multispace BDDC
We show several different ways the original twolevel 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 34 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 precorrection
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 postcorrection
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 = (IP)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, (342)
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 aorthogonal
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 aorthogonal 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 aorthogonal.
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 un Â£ 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
Aii 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 offdiagonal 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
\\(IP)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
\\{IP)Ew\\l ^ \\{IP)E(IP)w\\l
IMIa " \\(IP)M\l
from E (/ P) = E, which follows from (3.23), and from u;a > (7 F)u;a,
which follows from the aorthogonality of the projection P.
Before proceeding into the Multilevel BDDC section, let us write concisely
the spaces and operators involved in the twolevel 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 twolevel (top)
and the threelevel (bottom) BDDC methods.
35
4. Multilevel BDDC
In this chapter, we generalize the twolevel 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 threelevel 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 aorthogonal 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 = (IPi)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 = (IPi)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 precorrection 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 = mnL1,
otherwise set up the righthand 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 postcorrection on level i,
vn Un : a (vu, zti) = a (uBi, zH), Vzn G Un. (4.12)
Apply the combined corrections,
Ui = un + uBi vn. (413)
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 41 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)
F2L4 = {1 Pl2)WAC2, V2L3 U]L1,
y2L2 {I Pli)WaL1, y2Ll = WnLl,
41
Qi I, Q2 Q3 (I Pi) Ei,
Q4 = Qs = (/ Pi) Ex (/ P2) E2,
(4.15)
Q2L4 Q2L3 (I Pi) Ei (/ P12) El2,
Qil2 = Q2L1 = (I Pi) Ex (/  Pli) ELi,
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 v2l3 = Uili
So, let zn Â£ Uji, i = 2,..., L 1, be arbitrary. From (4.6) using (4.10) and
(47),
a(un, zn) = (ri, zH) = (rBi1, Ei^zn) = (4.16)
= (rj_i, EiiZn) aiuni^iizn)
Since from (4.6), using the fact that PiiEiiZn Un1, it follows that
(rxi,PiiEiiZn) a(%in\,Pi\Ei\Zji) 0,
we get from (4.16),
(^/u Zn) (rj_i, (/ Pi1) EiiZji) a (Ujii, (/ Pii)EiiZn),
and because a (u{i1, (I Pii)EiiZn) = 0 by orthogonality, we obtain
a{un, zn) = (r<_i, (/ Pt1) EiXzH)
42
Repeating this process recursively using (4.16), we finally get
a(uIi,zIi) = (rii, (/ Pi~i) EiiZu) = ...
= (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
V2L2 = u>ali 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 aorthogonal projection, it holds that
a {PiWAi, PwAi) = a (wAi, PiWAi) = (rBi, EiPiWAi) = 0,
43
where we have used EJJa 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\ WULi is defined by (4.9),
then i>2li = whli 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) Eli {wAli + tdiLi)] >
and, in particular for U\,
U\ = un +
+ (I Pi) Ei [u>ai 4... 4 (/ Pl 1) El 1 (wAli + ^ni,1)]
= QlVl 4 Q2V2 4 ... 4 Q2L2V2L2 + Q2L\V2L\
It remains to verify the assumptions of Corollary 3.4.
The spaces Wm and WAi, for alii = 1,..., L 1, are aorthogonal by (4.2)
and from (4.3),
(i Pi)wAi cwAi,
thus (I Pi) WAi is aorthogonal 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 aorthogonal.
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
= (IPi)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, (419)
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 + ... + Wali + Wnci,
un+1 (//i+i, wal\ (I Pc,i) Â£ WnLi)
and in (4.21) as
Ui+\ un+1 + ... + (/ Pj+i) Ej+i (/ Pl_i) Pl_i (?uali + f^nLi),
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 ~ Plx)ELXwULX\\1
then the Multilevel BDDC preconditioner (Algorithm 41) 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^nti = ^2Li V2L1 = Wnti*
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 twolevel 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 41) 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 twolevel BDDC method. It was presented in [61] starting
from corner constraints only, formulated in terms of FETIDP, 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 onetoone,
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
twolevel BDDC satisfies k < u, where
CR = RcQp.
(5.3)
uc = QpV, w = Rv.
U) = sup
wW
\\(IP)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 righthand side
of (54) 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 (CourantFisherWeyl 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 T7T
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 lefthand side of (5.6) is symmetric positive
semidefinite and the bilinear form on the righthand 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,
'\\(IP)(IE)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 zeroone 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(ie)(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 righthand 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 righthand 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 lefthand side is symmetric positive semidefinite and if the
pair of substructures intersects boundary with Dirichlet boundary conditions, the
matrix on the righthand 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 righthand
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 twolevel 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 semidefinite 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 inversefree 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 matrixfree format, i.e., as functions for matrixvector
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 righthand side only positive semidefinite. Initial experiments reveal that
the nonpreconditioned 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 RayleighRitz 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 righthand 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 nonpreconditioned (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 lefthand side, and Euclidean norms of residuals
for different eigenpairs are in the panels on the righthand 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 noilpreconditioned (top) vs. preconditioned
LOBPCG (bottom) for the detection of spurious rigidbody modes of the pair of
subdomains from Fig. 5.2. Estimated eigenvalue errors are in the panels on the
lefthand side, and Euclidean norms of residuals for different eigenpairs are in the
panels on the righthand 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 lefthand side, and Euclidean norms of residuals
for different eigenpairs are in the panels on the righthand 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 41 satisfies
uj = nffui, (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 40 [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
apriori 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 twolevel routines without any modifications.
68
Finally, we remark that instead of interior pre and postcorrection on the
lowest decomposition level, cf. equations (4.6)(4.7) and (4.12)(4.13), we reduce
the problem to interfaces in the preprocessing 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
twolevel 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 twodimensional results are reproduced
from [78], the threedimensional results appear here for the first time.
7.1 Twodimensional 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 thirdlevel.
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 twolevel 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 threelevel decomposition (Fig. 7.1). Because the
number of coarse degrees of freedom depends on apriori 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 twolevel method (Table 7.5) with the threelevel 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 twolevel 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 secondlevel decomposition (bottom) by the thick line.
72
Table 7.1: Results for the planar elasticity from Fig. 7.1 obtained using non
adaptive 2level 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 108.
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 3level 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 threelevel 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 2level 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 2level 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 3level 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 Threedimensional 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 noiladaptive 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 fillin of the Cholesky factor was quite high when compared with
the fillin 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, fillin is the relative fillin 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. zerolevel of: nnz(R) illin th fillin e method c cond id not iter
no fillin 3922 726 1.01  oo
1 103 9 784 734 2.53 4.14 106 331
1 io4 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 nonadaptive
2level 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 2level method. When the
nonadaptive 2level is replaced by the 3level method, Tables 7.10 and 7.12, the
condition number estimate as well as the number of iterations grows. However,
with the adaptive 3level approach (Table 7.13) we were able to achieve nearly
the same convergence properties as in the adaptive 2level method (Table 7.11).
Table 7.10: Results for the cube from Fig. 7.3 obtained using the nonadaptive
2level 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
2level 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 nonadaptive
3level 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
3level 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 nonadaptive 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 nonadaptive constraints for the de
composition into 1024 substructures are summarized in Tables 7.16 and 7.17 for
the two and threelevel 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 3level
method double, or even triple, compared to the 2level 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.47.7.
83
