ISOGEOMETRIC BDDC PRECONDITIONERS WITH DELUXE ...

ISOGEOMETRIC BDDC PRECONDITIONERS WITH DELUXE SCALING ˜ DA VEIGA ∗ , LUCA F. PAVARINO ∗ , SIMONE SCACCHI LOURENCO BEIRAO OLOF B. WIDLUND † , AND STEFANO ZAMPINI ‡

∗

,

TR2013-955 Abstract. A BDDC (Balancing Domain Decomposition by Constraints) preconditioner with a novel deluxe scaling, introduced by Dohrmann for problems with more than one variable coefficient, is extended to Isogeometric Analysis of scalar elliptic problems. This new scaling turns out to be much more powerful than BDDC with standard ρ- and stiffness scaling studied in a previous study of isogeometric BDDC. Our h-analysis shows that the condition number of the resulting deluxe BDDC preconditioner is scalable with a quasi-optimal polylogarithmic bound which is also independent of coefficient discontinuities across subdomain interfaces. Extensive numerical experiments support the theory and show that the deluxe scaling can yield a remarkable improvement over the older scalings, in particular, for large polynomial degree and high regularity of the isogeometric elements. Key words. domain decomposition, BDDC preconditioners, isogeometric analysis, elliptic problems AMS subject classifications. 65F08, 65N30, 65N35, 65N55

1. Introduction. The design of efficient and scalable iterative solvers for Isogeometric Analysis (IGA, see e.g. [?, ?, ?]) is a far from routine due to the integration of Finite Element Analysis and Computer Aided Design techniques which are required to build smooth and high-order discretizations based on nonuniform rational B-splines (NURBS) representations for both the domain geometry and finite element basis functions. In particular, there are complications arising from the fact that the basis functions are not nodal, which leads to wide (fat) interfaces. The main goal of this paper is to design and analyze a BDDC (Balancing Domain Decomposition by Constraints, see [?, ?]), preconditioner for Isogeometric Analysis based on a novel type of interface averaging, which we will denote by deluxe scaling. This variant was recently introduced by Dohrmann and Widlund in a study of H(curl) problems, see [?] and also [?] for its application to problems in H(div) and for Reissner–Mindlin plates [?]. In our previous work on isogeometric BDDC [?], standard BDDC scalings were employed with weights for the averaging built directly from the values of the elliptic coefficients in each subdomain (ρ-scaling) or from the values of the diagonal elements of local and global stiffness matrices (stiffness scaling). The novel deluxe scaling, originally developed to deal with elliptic problems with more than one variable coefficient, is instead based on solving local problems built from local Schur complements associated with sets of what are known as the dual variables. This new scaling turns out to be much more powerful than the standard ρ- and stiffness scaling even for scalar elliptic problems with one variable coefficient. The main result of our h-analysis shows that the condition number of the resulting deluxe ∗ Dipartimento di Matematica, Universit` a di Milano, Via Saldini 50, 20133 Milano, Italy. E-mail: [email protected], [email protected], [email protected]. This work was supported by grants of M.I.U.R. (PRIN 200774A7LH 003). † Courant Institute of Mathematical Sciences, 251 Mercer Street, New York, NY 10012. [email protected],http://cs.nyu.edu/cs/faculty/widlund/index.html. This work was supported in part by the U.S. Department of Energy under contract DE-FG02-06ER25718 and in part by the National Science Foundation Grant DMS-1216564. ‡ CINECA, Italy. E-mail: [email protected].

1

2

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

BDDC preconditioner satisfies the same quasi-optimal polylogarithmic bound in the ratio H/h of subdomain to element diameters, as in [?], and that it is independent of the number of subdomains and jumps of the coefficients of the elliptic problem across subdomain interfaces. Moreover, our numerical experiments with deluxe scaling show a remarkable improvement, in particular, for increasing polynomial degree p of the isogeometric elements, regardless of the element regularity k. In particular, for 2D problems, the convergence rate of deluxe BDDC appears to be independent of p, while for 3D problems, we observe equally strong results for p ≤ 5; even for larger p, we find several orders of magnitude better results than what has been obtained with stiffness and ρ-scaling. Recent work on IGA preconditioners have focused on overlapping Schwarz preconditioners, [?, ?, ?], multigrid methods, [?], and nonoverlapping preconditioners, [?, ?]. Among the recent extensions of BDDC methods, we mention the work on mortar discretizations, [?, ?], discontinuous Galerkin methods, [?], advection-diffusion and indefinite problems, [?, ?], inexact solvers, [?, ?], Reissner–Mindlin plates, [?, ?], spectral elements, [?], and multilevel algorithms, [?, ?]. We remark that we could also consider FETI-DP algorithms, see, e.g., [?, ?] defined with the same set of primal constraints as our BDDC algorithm, since it is known that then the BDDC and FETI-DP operators have the same eigenvalues with the exception of at most two; see [?, ?, ?]. The rest of this paper is organized as follows. We recall the basics on IGA discretizations of elliptic problems in Sec. ??. In Sec. ??, we introduce the domain and space decompositions in the isogeometric context, the required restriction and interpolation operators, the deluxe scaling, and construct the BDDC preconditioner. In Sec. ??, we prove a condition number bound for the deluxe BDDC preconditioned operator. In Sec. ??, the results of serial and parallel numerical tests in two and three dimensions are presented, confirming our theoretical estimates. 2. Isogeometric discretization of scalar elliptic problems. We consider the model elliptic problem: −∇ · (ρ∇u) = f in Ω,

Lorenzo: could you add a multi-patch comment?

u = 0 on ∂Ω,

(2.1)

where ρ is a scalar field satisfying 0 < ρmin ≤ ρ(x) ≤ ρmax ∀x ∈ Ω and where Ω ⊂ Rd , d = 2, 3, is a bounded and connected CAD domain. We discretize (??) with IGA based on B-splines and NURBS basis functions; see, e.g., [?]. When we describe our problem and the iterative method, we will confine our discussion, for simplicity, mostly to the two-dimensional single-patch case but we will briefly comment on extensions to three dimensions and multi-patch domains. The bivariate B-spline discrete space is defined by p,q Sbh := span{Bi,j (ξ, η), i = 1, . . . , n, j = 1, . . . , m},

(2.2)

p,q where the bivariate B-spline basis functions Bi,j (ξ, η) = Nip (ξ) Mjq (η) are defined by tensor products of one-dimensional B-splines functions Nip (ξ) and Mjq (η) of degree p and q, respectively. Analogously, the NURBS space is the span of NURBS basis functions defined in one dimension by

N p (ξ)ωi N p (ξ)ωi = i , Rip (ξ) := Pn i p w(ξ) ˆı=1 Nˆı (ξ)ωˆı

(2.3)

3

Deluxe BDDC Methods for Isogeometric Analysis

with the weight function w(ξ) := tensor product p,q Ri,j (ξ, η)

Pn

ˆı=1

Nˆıp (ξ)ωˆı ∈ Sbh , and in two dimensions by a

p,q p,q Bi,j (ξ, η)ωi,j Bi,j (ξ, η)ωi,j P P = := n , m p,q w(ξ, η) B (ξ, η)ω ˆ ˆı=1 ˆı,ˆ j j=1 ˆı,ˆ j

(2.4)

where w(ξ, η) is the weight function and ωi,j = (Cω i,j )3 the weights associated with a n × m net of control points Ci,j . The discrete space of NURBS functions on the domain Ω is defined as the span of the push-forward of the NURBS basis functions (??) (see, e.g., [?, ?]) p,q Nh := span{Ri,j ◦ F−1 , with i = 1, . . . , n; j = 1, . . . , m},

(2.5)

b → Ω, the geometrical map between parameter and physical spaces with F : Ω F(ξ, η) =

n X m X

p,q Ri,j (ξ, η)Ci,j .

(2.6)

i=1 j=1

For simplicity, we will consider the case with a Dirichlet boundary condition imposed on all of ∂Ω and can then define the spline space in the parameter space by b 2 = [span{B p,q (ξ, η), i = 2, . . . , n − 1, j = 2, . . . , m − 1}]2 . Vbh := [Sbh ∩ H01 (Ω)] i,j and the NURBS space in physical space as p,q Uh := [Nh ∩ H01 (Ω)]2 = [span{Ri,j ◦ F−1 , with i = 2, . . . , n − 1; j = 2, . . . , m − 1}]2 .

The IGA formulation of problem (??) then reads: ( Find uh ∈ Uh such that: a(uh , vh ) =< f, vh > with the bilinear form a(uh , vh ) =

R Ω

∀v ∈ Uh ,

(2.7)

ρ∇uh ∇vh dx.

3. BDDC preconditioners for the Schur complement system. When using iterative substructuring methods, such as BDDC, we first reduce the problem to one on the interface by implicitly eliminating the interior degrees of freedom, a process known as static condensation; see, e.g., Toselli and Widlund [?, Ch. 4]. 3.1. Knots and subdomain decomposition. A decomposition is first built for the underlying space of spline functions in the parametric space, and is then easily extended to the NURBS space in the physical domain. From the full set of knots, {ξ1 = 0, ..., ξn+p+1 = 1}, we select a subset {ξik , k = 1, . . . , N + 1} of non-repeated knots with ξi1 = 0, ξiN +1 = 1. The interface knots are given by ξik for k = 2, .., N and they define a decomposition of the closure of the reference interval into subdomains 
Ib = [0, 1] =

[

 Ibk , with Ibk = (ξik , ξik+1 ),

k=1,..,N

that we assume to have a similar length, i.e., Hk := diam(Ibk ) ≈ H.

4

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

In more than one dimension, we just use tensor products. Thus, in two dimension, we define the subdomains by Ibk = (ξik , ξik+1 ),

b kl = Ibk × Ibl , Ω

Ibl = (ηjl , ηjl+1 ),

1 ≤ k ≤ N1 , 1 ≤ l ≤ N2 . (3.1)

For simplicity, we reindex the subdomains using only one index to obtain the decomposition of our domain [

b= Ω

bk, Ω

k=1,..,K

into K = N1 N2 subdomains; analogously, into K = N1 N2 N3 subdomains in three dimensions. Throughout this paper, we assume that both the subdomains and elements defined by the coarse and full sets of knot vectors are shape regular and with quasi-uniform characterstic diameters H and h, respectively. 3.2. The Schur complement system. As in classical iterative substructuring, we reduce the problem to one on the interface Γ :=

K [

 b k \∂ Ω b ∂Ω

k=1

by static condensation, i.e., by eliminating the interior degrees of freedom associated with the basis functions with support in each subdomain. The resulting Schur b k and its local interface Γk := ∂ Ω b k \ ∂Ω b will be denoted by S (k) . complement for Ω In the sequel, we will use the following sets of indices: ΘΩ = {(i, j) ∈ N2 : 2 ≤ i ≤ n − 1, 2 ≤ j ≤ m − 1}, p,q ΘΓ = {(i, j) ∈ ΘΩ : supp(Bi,j ) ∩ Γ 6= ∅}, p,q ΘC = {(i, j) ∈ ΘΓ : supp(Bi,j ) ∩ C 6= ∅},

where C denotes the set of vertices on Γ. We note that ΘΓ typically consists of several layers of knots reflecting the fact that a knot on Γ belongs to the support of a number of B-spline basis functions. We therefore talk about a fat interface and fat edges, etc. The discrete interface space is defined as p,q VbΓ := span{Bi,j , (i, j) ∈ ΘΓ },

b k are defined as and the local spaces of functions with support in Ω (k)

VI

b k ). := Vbh ∩ H01 (Ω

The space Vbh can be decomposed as (k)

⊕K k=1 VI

⊕ H(VbΓ ),

(3.2)

where H : VbΓ → Vbh , is the piece-wise discrete spline harmonic extension operator, which provides the minimal energy extension of values given on VbΓ . The interface component of the discrete solution satisfies the reduced system s(uΓ , vΓ ) =< fb, vΓ >,

∀vΓ ∈ VbΓ ,

(3.3)

5

Deluxe BDDC Methods for Isogeometric Analysis

with suitable right-hand side fb and Schur complement bilinear form is defined by s(wΓ , vΓ ) := a(H(wΓ ), H(vΓ )).

(3.4)

For simplicity, in the sequel we will drop the subscript Γ for functions in VbΓ . In matrix form, (??) is the Schur complement system SbΓ w = fb,

(3.5)

−1 b where SbΓ = AΓΓ − AΓI A−1 II AΓI , f = fΓ − AΓI AII fI , are obtained from the original discrete problem by reordering the spline basis functions into sets of interior (subscript I) and interface (subscript Γ) basis functions      AII AΓI wI fI = . (3.6) AΓI AΓΓ wΓ fΓ

The Schur complement system (??) is solved by a Preconditioned Conjugate Gradient (PCG) iteration, where only the action of SbΓ on a vector is needed and SbΓ is never explicitly formed. In fact, this action can be assembled from actions of the subdomain Schur complements S (k) on subvectors. The preconditioned Schur complement system solved by PCG is then −1 b −1 MBDDC SΓ w = MBDDC fb,

(3.7)

−1 where MBDDC is the BDDC preconditioner, defined in (??) below. Before we can con−1 struct MBDDC , we need to introduce some restriction and scaling operators associated with the subspace decompositions.

3.3. Subspace decompositions. Let V (k) be the local discrete space defined b k and with functions that vanish on ∂ Ω b k ∩ ∂ Ω. b Analogously to on the subdomain Ω the space splitting (??), we split this local space into a direct sum of its interior (I) (k) L (k) and interface (Γ) subspaces V (k) = VI VΓ , where (k)

VI (k)

VΓ

(k)

(k)

p,q := span{Bi,j , (i, j) ∈ ΘI }, ΘI

(k)

p,q b k }, := {(i, j) ∈ ΘΩ : supp(Bi,j )⊂Ω

(k)

p,q p,q b k ∩ Γk 6= ∅}, := span{Bi,j , (i, j) ∈ ΘΓ }, ΘΓ := {(i, j) ∈ ΘΓ : supp(Bi,j ) ∩ (∂ Ω

and we define the associated product spaces by VI :=

K Y k=1

(k)

VI ,

VΓ :=

K Y

(k)

VΓ .

k=1

The functions in VΓ are generally discontinuous (multi-valued) across Γ, while our isogeometric approximations belong to VbΓ , the subspace of VΓ consisting of functions continuous across Γ. We will select some interface basis functions as primal (subscript Π), that will be made continuous across the interface and will be subassembled between their supporting elements, and we will call dual (subscript ∆) the remaining interface degrees of freedom that can be discontinuous across the interface and which vanish at the primal degrees of freedom. This splitting allows us to decompose each

6

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini (k)

local interface space into primal and dual subspaces VΓ the associated product spaces by V∆ :=

K Y k=1

(k)

V∆ ,

VΠ :=

K Y

(k)

= VΠ

L

(k)

V∆ , and define

(k)

VΠ .

k=1

We also need an intermediate subspace VeΓ ⊂ VΓ of partially continuous basis functions M VeΓ := V∆ VbΠ , where the product space V∆ has been defined above and VbΠ is a global subspace consisting of the selected primal variables. Particular choices of primal set are given in Sec. ?? below. In our analysis, we will focus on the simplest primal set consisting of basis functions associated with the subdomain vertices, but in our numerical experiments we will also study richer choices using sets of edge/face basis functions with constant values at the knots of the associated edge/face, see Sec. ??. In order to define our preconditioners, we will need the following restriction and interpolation operators represented by matrices with elements in the set {0, 1}: eΓ∆ : VeΓ −→ V∆ , R (k) (k) R∆ : V∆ −→ V∆ ,

eΓΠ : VeΓ −→ VbΠ , R (k) (k) RΠ : VbΠ −→ VbΠ

bΠ : VbΓ −→ VbΠ R b(k) : VbΓ −→ V (k) . R ∆

(3.8)

∆

For any edge/face F, we also define RF as the restriction matrix to F. 3.4. Deluxe scaling (see Dohrmann and Widlund [?]). Let Ωk be any subdomain in the partition, k = 1, 2, ..., K. We will indicate by Ξk the index set of all the Ωj , j 6= k, that share an edge F with Ωk . For regular quadrilateral subdomain partitions in two dimensions, the cardinality of Ξk is 4. In BDDC, the average w ¯ := ED w of an element in w ∈ VeΓ , is computed separately for the sets of interface degrees of freedom of certain different equivalence classes. Each of these classes is defined by the set of indices of the subdomains boundaries to which the degrees of freedom belong. In order to define the deluxe scaling, we first consider a class F with only two elements, k, j, as for an edge in two dimensions or a face in three dimensions. We (k) (j) define two principal minors, SF and SF , obtained from S (k) and S (j) by removing all rows and columns which do not belong to the degrees of freedom which are common to the (fat) boundary of Ωk and Ωj . (k) Let wF := RF w(k) ; the deluxe average across F is then defined as  −1   (k) (j) (k) (k) (j) (j) w ¯F = SF + SF SF wF + SF wF .

(3.9)

If the Schur complements of an equivalence class have small dimensions, they can be computed explicitly, otherwise their action can be computed by solving Dirichlet problems on the union of the relevant subdomains. In this subsection, we will not consider cases where there are additional equivalence classes of dual unknowns that are associated with more than two subdomains; see Subsec. ?? for the general case. Each of the relevant equivalence classes, which involve the subdomain Ωk , will contribute to the values of w. ¯ Each of these contributions will belong to VbΓ , after being extending by zero to Γ \ F; the resulting element is

Deluxe BDDC Methods for Isogeometric Analysis

7

T given by RF w ¯F . We then add the contributions from the different equivalence classes to obtain

w ¯ = ED w = wΠ +

X

T RF w ¯F ,

F

and its complementary projection PD w := (I − ED )w := PD w = w∆ −

X

T RF w ¯F .

F

In order to rewrite ED in matrix form, we define, for each subdomain Ωk , the associated scaling matrix    D(k) =   

(k)



DF j1

(k)

DF j2 ..

. (k)

  ,  

DF jk where j1 , j2 , . . . , jk ∈ Ξk and the diagonal blocks are given by the deluxe scaling  −1 (k) (k) (j) (k) DF := SF + SF SF . In case of more than two subdomains, the diagonal blocks are replaced by the analogous deluxe scaling expressions given in Sec. ??. (k) (k) e (k) We can then define the scaled local operators by RD,Γ := D(k) RΓ , R D,∆ := (k)

(k)

RΓ,∆ RD,Γ and the global scaled operator eD,Γ := the direct sum R b Π ⊕k R e(k) , R D,∆

(3.10)

so that the averaging operator is T eΓ R eD,Γ ED = R ,

eΓ := R b Π ⊕k R e(k) . where R ∆ 3.5. Choice of primal constraints. The choice of primal degrees of freedom is fundamental for the construction of efficient BDDC preconditioners. We will consider the following three primal choices. p,q i) VbΠ = VbΠC : set of coefficient values of Bi,j , (i, j) ∈ ΘC associated with the subdomain corners. This choice is not always sufficient to obtain scalable and fast preconditioners, in particular for problems in three dimensions, and this has motivated the search for richer primal sets that may yield faster preconditioners. ii) VbΠ = VbΠC+E : previous set augmented with the subdomain edge averages. iii) VbΠ = Vb C+E+F : previous set augmented with the subdomain face averages. Π

When working with the primal sets VbΠC+E and VbΠC+E+F employing edge and/or face averages we will assume that, after a change of basis, each primal variable corresponds to an explicit degree of freedom, see [?], and Sec. ?? for implementation details. Our theoretical analysis in Sec. ?? will focus on the simplest primal set VbΠC .

8

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

3.6. The BDDC preconditioner. We denote by A(k) the local stiffness matrix b k . By partitioning the local degrees of freedom into those restricted to subdomain Ω interior (I) and those interface (Γ), as before, and by further partitioning the latter into dual (∆) and primal (Π) degrees of freedom, then A(k) can be written as  (k)  (k)T (k)T " # AII A∆I AΠI (k) (k)T T  AII AΓI  A(k) = =  A(k) A(k) A(k)  . (k) (k) ∆I ∆∆ Π∆ AΓI AΓΓ (k) (k) (k) AΠI AΠ∆ AΠΠ Using the scaled restriction matrices, defined in (??) and (??), the BDDC preconditioner can be written as −1 T eD,Γ eD,Γ , MBDDC =R SeΓ−1 R

(3.11)

where  T  eΓ∆ SeΓ−1 = R

K h X

T

0

(k)

R∆

i

"

k=1

(k) AII (k) A∆I

(k)T A∆I (k) A∆∆

#−1 

0 (k) R∆

  eΓ∆ + ΦS −1 ΦT . R ΠΠ

(3.12) b k , with a The first term in (??) is the sum of local solvers on each subdomain Ω Neumann condition on the local ∆ edges and with the coarse degrees of freedom constrained to vanish. The second term is a coarse solver for the primal variables, that we implemented as in [?, ?] by using the coarse matrix  # " #−1 " K (k)T h i (k) (k)T X A A A (k)T  (k) (k) ΠI II ∆I  R(k) SΠΠ = RΠ AΠΠ − A(k) AΠ∆ Π (k) (k) ΠI (k)T A∆I A∆∆ AΠ∆ k=1 and a matrix Φ mapping primal degrees of freedom to interface variables, given by # " #−1 " K h (k)T i (k) (k)T X T A A A (k) T T (k) ΠI II ∆I e −R RΠ . Φ=R 0 R∆ ΓΠ Γ∆ (k) (k) (k)T A A AΠ∆ ∆I ∆∆ k=1 The columns of Φ represent the coarse basis functions defined as the minimum energy extension, using the original bilinear form, into the subdomains and subject to the chosen set of primal constraints. 4. Condition number bounds. In this section, we will derive a condition number bound for our preconditioned operator. We will make use of the following preliminary result, that is an immediate combination of Theorems 6.1 and 6.2 in [?], re-written in the present notation. Theorem 4.1. Let the dimension of the problem be equal to 2. Then for all (k) subdomains Ωk it exists a boundary seminorm | · |W k such that for all w(k) ∈ VΓ X |RF w(k) |2W F , |w(k) |2W = (4.1) k

F ∈∂Ω

|w(k) |2W k

≤ C (w(k) )T S (k) w(k) ,

(4.2)

with the positive constant C independent of h, H and where |w|W F is a seminorm on the space of discrete functions associated with the edge F. Moreover, for all w(k) ∈ (k) VΓ that vanish at the subdomain primal (corner) degrees of freedom, it holds 2 (w(k) )T S (k) w(k) ≤ C(1 + log2 (H/h))|w(k) |W , k

(4.3)

9

Deluxe BDDC Methods for Isogeometric Analysis

with the positive constant C independent of h, H. For a proof and an explicit form of this seminorm, see [?]. A bound on the condition number of the preconditioned operator can be obtained eΓ R eT . We recall by bounding the SeΓ -norm of the average operator defined by ED = R D,Γ −1 eT Se−1 R eD,Γ and Sb = R eT SeΓ R eΓ . that MBDDC =R D,Γ Γ Γ Lemma 4.2. (Lower bound) uT MBDDC u ≤ uT SbΓ u, ∀u ∈ VbΓ .

(4.4)

eT R e Proof. Let w = MBDDC u. Since, as is easy to show, R Γ D,Γ = I, we have eΓT R eD,Γ w = uT R eΓT SeΓ Se−1 R eD,Γ w uT MBDDC u ≤ uT w = uT R Γ 1/2

1/2

Γ

Γ

eΓ u, R eΓ u) (Se−1 R eD,Γ w, Se−1 R eD,Γ w) ≤ (R Γ Γ e e S S

T eD,Γ w)1/2 eΓT SeΓ R eΓ u)1/2 (wT R eD,Γ = (uT R SeΓ−1 SeΓ SeΓ−1 R

= (uT SbΓ u)1/2 (uT MBDDC u)1/2 . Lemma 4.3. (Upper bound) If |ED u|2Se ≤ CE |u|2Se Γ

Γ

∀u ∈ VeΓ , then

uT SbΓ u ≤ CE uT MBDDC u.

(4.5)

eΓ M −1 eT SeΓ R eΓ u = uT R eT SeΓ R Proof. uT SbΓ u = uT R Γ Γ BDDC MBDDC u T eΓT SeΓ R eΓ R eD,Γ eD,Γ w =uT R SeΓ−1 R

eD,Γ w)1/2 eD,Γ w, ED Se−1 R eΓ u, R eΓ u)1/2 (ED Se−1 R ≤(R Γ Γ e e S S Γ

Γ

eT SeΓ R eD,Γ w, Se−1 R eD,Γ w)1/2 eΓ u)1/2 C 1/2 (Se−1 R ≤(u R Γ E Γ Γ e S T

Γ

1/2 =CE (uT SbΓ u)1/2 (uT MBDDC u)1/2 ,

where the last step follows as in the proof of Lemma ??. With these two lemmas, we can prove our main theoretical result. Theorem 4.4. Let the dimension of the problem be equal to 2 and the primal set be given by the subdomain corner set VbΠC . Then it holds   −1 Cond MBDDC SbΓ ≤ C(1 + log(H/h))2 , with the positive constant C independent of h, H and the  jumps of the coefficient ρ. −1 Proof. Lemma ?? provides the lower bound λmin MBDDC SbΓ ≥ 1. Lemma ??   −1 provides the upper bound λmax MBDDC SbΓ ≤ CE if we prove the bound on ED required in the hypothesis of Lemma (??). We prove the quivalent bound on the PK ¯ 2 complementary operator PD = I − ED , and since |ED w|2Se = k=1 |R Γ ED w|S (k) , we Γ ¯ Γ PD w|2 (k) . will focus on proving a bound for |R S We recall that the deluxe average across an edge/face F shared by two subdomains (k) (j) (k) (k) (j) (j) (k) has been defined in (??) by w ¯F = (SF + SF )−1 (SF wF + SF wF ), where wF :=

10

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

RF w(k) and RF is the restriction matrix to the edge/face F. This element w ¯F is then T extended by zero to Γk \ F to obtain the element RF w ¯F in VbΓ . Recalling that Ξk denotes the index set of the subdomains sharing an edge or a face with Ωk , by adding the relevant contributions from the different edges/faces, we obtain X X T T w ¯ = ED w = wΠ + RF w ¯F , hence PD w = w∆ − RF w ¯F , so that F ∈Ξk

F ∈Ξk

X

¯ Γ PD w|2 (k) ≤ |Ξk | |R S

(k)

T |RF (wF − w ¯F )|2S (k) ,

(4.6)

F ∈Ξk

where |Ξk | = 4 in our special two dimensional case; in the general case we will use that the number of neighbors is finite. Still focusing, for simplicity, on the case where F is shared by two subdomains (an edge in two dimensions or a face in three), we find, by simple algebra, that (k)

(k)

(j)

(j)

(k)

(j)

wF − w ¯F = (SF + SF )−1 SF (wF − wF ), (k)

(k)

(j)

(j)

so that (k)

(j)

T T ¯F )|2S (k) = |RF |RF (wF − w (SF + SF )−1 SF (wF − wF )|2S (k) .

(4.7)

By adding and subtracting a suitable coarse function wΠ ∈ VbΓ (a specific choice will be given below), we have (k)

(j)

(k)

(j)

(k)

(k)

(j)

(j)

wF −wF = wF −RF wΠ −(wF −RF wΠ ) = wF −RF wΠ −(wF −RF wΠ ). (4.8) (k)

T Using (??) in (??) and noting that RF S (k) RF = SF , we then find that (k)

(k)

(k)

T T T |RF (wF − w ¯F )|2S (k) = (RF (wF − w ¯F ))T S (k) (RF (wF − w ¯F )) (k)

(j)

(j)

(k)

(j)

(k)

(k)

(j)

(j)

(k)

(j)

= (wF − wF )T SF (SF + SF )−1 SF (SF + SF )−1 SF (wF − wF )) (k)

(k)

(j)

(k)

(j)

(k)

(k)

(j)

(j)

(k)

(k)

≤ 2(wF − RF wΠ )T SF (SF + SF )−1 SF (SF + SF )−1 SF (wF − RF wΠ ) (j)

(j)

(j)

(k)

(j)

(k)

(k)

(j)

(j)

(j)

(j)

+ 2(wF − RF wΠ )T SF (SF + SF )−1 SF (SF + SF )−1 SF (wF − RF wΠ ). We can simplify this expression by using the two bounds (j)

(k)

(j)

(k)

(k)

(j)

(j)

(k)

(j)

(k)

(j)

(k)

(k)

(j)

(j)

(j)

SF (SF + SF )−1 SF (SF + SF )−1 SF ≤SF , SF (SF + SF )−1 SF (SF + SF )−1 SF ≤SF , which follow easily by considering the action of these operators on any eigenvector of (k) (j) the generalized eigenvalue problem SF φ = λSF φ and just using that all eigenvalues are strictly positive. Hence, we obtain (k)

(k)

(k)

(k)

(j)

(j)

(j)

(k)

(k)

T |RF (wF − w ¯F )|2S (k) ≤ 2(wF − RF wΠ )T SF (wF − RF wΠ ) + (j)

(j)

2(wF − RF wΠ )T SF (wF − RF wΠ ).

to Olof: ok reference to standard edge lemma?

(4.9)

What remains, and this is where analysis rather than linear algebra enters the picture, is to establish an edge lemma (see [?] for a finite element analog)

11

Deluxe BDDC Methods for Isogeometric Analysis (k)

(k)

(k)

(k)

(k)

(wF − RF wΠ )T SF (wF − RF wΠ ) ≤ C(1 + log(H/h))2 (w(k) )T S (k) w(k) , (4.10) (k)

where the values of w(k) shares the values on F with wF but are otherwise arbitrary. Note that this is where the logarithmic factors enters; the energy of a minimal extension will be smaller than that of the extension by zero. To obtain the estimate (??), we select the function wΠ introduced in (??) as the function of VbΓ with the same primal values as w and with values on the edges of (k) each Ωk such that the discrete seminorm |wΠ |W k introduced in Theorem ?? (see also [?, eq. (6.13)]) is minimized. Note that, since the norm | · |W k is the sum of edge-wise norms (??) and the values at the corner (primal) dofs are assigned, it is immediate to check that the minimization problem above corresponds to an edgeby-edge minimization of the seminorms | · |W F , for all F edges of the subdomain (k) partition. Therefore, wΠ is continuous across edges, i.e. wΠ ∈ VbΓ . The definition above guarantees that (k)

|wΠ |W k ≤ |w(k) |W k . (k)

(k)

(4.11)

(k)

Denoting zF = wF − RF wΠ , the left-hand side of (??) can be bounded by (k)T

zF

(k) (k)

(k)

(k)

T T SF zF =(RF zF )T S (k) (RF zF ) (k)

T 2 ≤C1 (1 + log(H/h))2 |RF zF |W ,

(4.12)

k

(k)

(k)

(k)

(k)

T T T RF (w(k) − wΠ ) (wF − RF wΠ ) = RF zF = RF where we have used (??). Since RF is nonzero only on the edge F and the discrete seminorm |·|W k is defined edge-by-edge as a sum of four terms, only the one associated with the common edge F is nonzero. Hence, thanks to (??), we have (k)

(k)

(k)

T 2 |RF zF |W ≤ |w(k) − wΠ |2W ≤ 2(|w(k) |2W + |wΠ |2W ) ≤ 4|w(k) |2W . k

k

k

k

k

By (??), we can then return to the local Schur bilinear form T

|w(k) |2W ≤ Cw(k) S (k) w(k) , k

and obtain the edge lemma (??) We can also compute, in exactly the same way, the contribution of the other subdomain Ωj as well as contributions from other relevant subdomain edges or faces. Summing over the edges of Ωk , we obtain from (??), (??), and (??), ¯ Γ PD w|2 (k) ≤ C(1 + log(H/h))2 |R S

X

|w(j) |2S (j) ,

j∈Ξk ∪{k}

¯ Γ ED w|2 (k) , and summing over the subdomains hence the same bound on |R S |ED w|2Se ≤ C(1 + log(H/h))2 |w|2Se , Γ

Γ

(4.13)

i.e. the ED bound required by Lemma ?? with CE = C(1 + log(H/h))2 , which proves the theorem.

to Olof: what do you want to do with the internal note? It is now in comment

12

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

4.1. The case of more than two subspaces. We have shown above how to (k) reduce an estimate of the norm of wF − w ¯F for the case of an equivalence class associated with just two subdomains. In many cases, we also need to work with equivalence classes with three or more subdomains; this is, e.g., the case if some of the degrees of freedom associated with the set ΩC are not primal, as for the reduced primal set VbΠC4 described below. This is also necessary for three dimensional problems unless we make the many degrees of freedom associated with all fat subdomain edges primal. We will first show that the same algebraic ideas can be used to take care of the case of three subdomains, Ωk1 , Ωk2 , and Ωk3 and an edge E. (k ) (k ) (k ) (k ) Letting SE 123 := SE 1 + SE 2 + SE 3 , our averaging operator now has the form (k1 )

(k123 ) −1

w ¯E

:= (SE

)

(k1 )

(SE

(k1 )

wE

(k2 )

+ SE

(k2 )

wE

(k3 )

+ SE

(k3 )

wE

),

As before, the required Schur complements can be computed explicitly if they have small dimensions, otherwise their action can be computed by solving Dirichlet problems on the union of the relevant subdomains. We find that (k1 )

wE

(k123 ) −1

(k1 )

−w ¯E

(k1 )

Since RE S (k1 ) RET = SE (k1 )

|RET (wE

)

= (SE

(k2 )

((SE

(k3 )

+ SE

(k1 )

(k2 )

− SE

)wE

(k2 )

wE

(k3 )

− SE

(k3 )

wE

).

, then the analog of (??) for three subdomains becomes (k123 ) −1

−w ¯E )|2S (k1 ) = |(SE

)

(k2 )

((SE

(k3 )

+SE

(k1 )

)wE

(k2 )

−SE

(k2 )

wE

(k3 )

−SE

(k3 )

wE

)|2 (k1 ) . SE

This norm can be bounded from above by the sum of the three terms: (k1 )T

3wE

(k2 )

(SE

(k2 )T

SE

(k3 )T

SE

+ 3wE + 3wE

(k3 )

+ SE

(k2 )

(SE

(k3 )

(SE

(k123 ) −1

)(SE

)

(k123 ) −1

SE

(k123 ) −1

SE

)

)

(k1 )

SE

(k1 )

(SE

(k1 )

(SE

(k123 ) −1

(SE

)

(k123 ) −1

SE

(k123 ) −1

SE

)

)

(k )T

(k )

(k2 )

(k3 )

(SE

(k2 )

wE

(k3 )

wE

+ SE

(k1 )

)wE

(k2 )

(k3 )

(4.14) .

(k )

The first term can be bounded by 3wE 1 SE 1 wE 1 by using the same argument (k ) as in the case of two subdomains and by working with the two matrices SE 1 and (k2 ) (k3 ) SE + SE . The second expression can also be reduced to an argument with the two matrices, (k2 ) (k ) (k ) SE and SE 1 + SE 3 , after first bounding the expression (??) from above by (k2 )T

3wE

(k )T

(k2 )

SE

(k )

(k123 ) −1

(SE

)

(k1 )

(SE

(k3 )

+ SE

(k123 ) −1

)(SE

)

(k2 )

SE

(k2 )

wE

.

(k )

A bound of 3wE 2 SE 2 wE 2 results. (k )T (k ) (k ) The third expression can be bounded by 3wE 3 SE 3 wE 3 in the same way. We then obtain the analog of (??) (k1 )

|RET (wE

paragaph can be eliminated?

(k1 )T

−w ¯E )|2S (k1 ) ≤ 3wE

(k1 )

(k1 )

(k2 )T

(k2 )

(k2 )

(k3 )T

(k3 )

(k )

wE 3 . (4.15) As in the case of two subdomains, we can also subtract a common element of the (k ) (k ) (k ) primal space from wE 1 , wE 2 , and wE 3 , to obtain (k1 )

|RET (wE

SE

wE

+ 3wE

SE

wE

(k )

(k1 )

(wE

(k1 )

− RE wΠ 1 ) +

(k )

(k2 )

(wE

(k2 )

− RE wΠ 2 ) +

(k )

(k3 )

(wE

(k3 )

− RE wΠ 3 ).

(k1 )

− RE wΠ 1 )T SE

(k2 )

− RE wΠ 2 )T SE

(k3 )

− RE wΠ 3 )T SE

−w ¯E )|2S (k1 ) ≤ 3(wE 3(wE 3(wE

+ 3wE

SE

(k )

(k )

(k )

(4.16)

13

Deluxe BDDC Methods for Isogeometric Analysis

Each of these terms can then be bounded by a counter part of the edge lemma (??). Turning to the four subdomain case, we have (k1 )

w ¯E

(k123 )

where SE

(k1234 ) −1

:= (SE

)

(k1 )

:= SE

(k1234 ) −1

(SE

)

(k2 )

((SE

(k1 )

(SE

(k2 )

+ SE

(k3 )

+ SE

(k1 )

wE

(k3 )

+ SE

(k4 )

+ SE

(k2 )

+ SE

(k4 )

+ SE (k1 )

)wE

(k2 )

wE

(k3 )

+ SE

(k3 )

wE

(k4 )

wE

(k1 )

−w ¯E

, and we find that wE (k2 )

− SE

(k2 )

wE

(k3 )

− SE

(k4 )

+ SE

(k3 )

wE

),

(k1 )

(k4 )

− SE

equals (k4 )

wE

).

(k )

The norm |RET (wE 1 − w ¯E )|2S (k1 ) can then be bounded, by using the same arguments as for three subdomains, by the sum of four terms, each with a coefficient 4. 4.2. A reduced primal set. In case of maximal regularity k = p − 1 of the isogeometric basis functions, the fat interface leads to a rich primal set VbΠC with p2 primal degrees of freedom for each subdomain vertex in 2D. We can then consider a reduced primal set VbΠC4 where out of these p2 degrees of freedom per vertex only the 4 corner degrees of freedom are retained as primal and the other p2 − 4 are considered dual. For this reduced primal set, we expect the same scalable convergence bound of Theorem ??, since the presence of two primal corners per edge allow us to still prove an edge lemma and the additional p2 − 4 dual variables per vertex are now shared by 4 subdomains and can the treated with the techniques of the previous Sec. ??. The numerical results presented in Sec. ?? confirm this scalable bound but show that the BDDC convergence rate with VbΠC4 primal set rapidly deteriorates with increasing p. 4.3. The three-dimensional case. We conjecture that a scalable convergence bound as in Theorem ?? holds also in 3D, but a complete proof is beyond the scope of this paper. We only note that the basic tools required are isogeometric edge and face lemmas in 3D, which we believe can be obtained by extending the 2D isogeometric techniques of our previous work [?], and the deluxe estimates in the case of more than two subdomains considered in Sec. ??. Remark. If we use a conventional averaging procedure, we start with a piecewise discrete harmonic function which is discontinuous across the interface and which is the result of solving a system with S˜ as matrix. By averaging point-wise across the interface, we introduce nonzero residuals next to the interface. We then have to remove them by solving a Dirichlet problem on each subdomain. If we use the new deluxe averaging, we obtain one contribution from each subdomain edge/face and in other applications additional contributions, e.g., associated with the subdomain vertices. Each of them will be discrete harmonic in the subdomains and therefore the local Dirchlet solves no longer will be necessary. 5. Numerical results. In this section, we report on numerical experiments with the isogeometric BDDC deluxe preconditioner for two and three dimensional elliptic model problems (??) on both the parametric (reference square or cube) and physical domains, discretized with isogeometric NURBS spaces with a mesh size h, polynomial degree p = q and regularity k. The domain is decomposed into K nonoverlapping subdomains of characteristic size H, as described in Sec. ??. In the following, we will denote by kΓ ≤ k the regularity at the subdomain interface. The discrete Schurcomplement problems are solved by the PCG method with the isogeometric BDDC preconditioner, with a zero initial guess and and a stopping criterion of a 10−6 reduction of the PCG residual. In the tests, we study how the convergence rate of the BDDC preconditioner depends on h, K, p, k, kΓ , and jumps in the coefficient of the

to Olof: should we drop this remark? If not, where should we place it?

14

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

a) homogeneous

b) central jump

c) checkerboard

Fig. 5.1. 2D quarter-ring domains used in the numerical tests. Examples with 4 × 4 subdomains

K 2×2 4×4 8×8 16 × 16

h = 1/16 κ2 it. 1.24 5

h = 1/32 κ2 it. 1.42 6 2.02 8

h = 1/64 κ2 it. 1.65 6 2.68 10 2.39 10

h = 1/128 κ2 it. 1.92 6 3.46 11 3.29 12 2.64 11

Table 5.1 BDDC deluxe preconditioner for a 2D quarter-ring domain: condition number κ2 and iteration counts it. as a function of the number of subdomains K and mesh size h. Fixed p = 3, k = 2.

details should be added

elliptic problem. In all tests, the BDDC condition number is essentially the maximum eigenvalue of the preconditioned operator, since its minimum eigenvalue is always very close to 1. The 2D tests have been performed with a MATLAB code based on the GeoPDEs library [?], while the 3D tests have been performed with a FORTRAN parallel code based on the PETSc library [?], in particular on the PCBDDC subroutine (contributed to the PETSC library by S. Zampini) and run on the BlueGene/Q cluster of the CINECA Laboratory (http://www.hpc.cineca.it/hardware/ibm-bgq-fermi). We start with the results of 2D experiments with deluxe BDDC with primal set VbΠC consisting of subdomain vertices. In the 3D experiments, we will consider primal sets consisting of subdomain vertices (VbΠC ), or subdomain vertices and edge averages (VbΠC+E ), or subdomain vertices and edge and face averages (VbΠC+E+F ). 5.1. Scalability in K and quasi-optimality in H/h in 2D. The condition number κ2 and iteration counts it. of the BDDC deluxe preconditioner are reported in Table ?? for a quarter-ring domain (see ?? a)), as a function of the number of subdomains K and mesh size h, for fixed p = 3, k = 2. The results show that the proposed preconditioner is scalable, since moving along the diagonals of the table the condition number appears to be bounded from above by a constant independent of K. Analogous results for a higher degree p = 5 and regularity k = 4 are reported in Table ??, where the condition numbers and iteration counts are even better than the ones for the lower degree case of Table ??. In contrast to the scalings proposed in our previous work [?], the BDDC deluxe preconditioner appears to retain a very good performance in spite of the increase of the polynomial degree p. To better understand this issue, we investigate how the performance depend on p in the next test.

15

Deluxe BDDC Methods for Isogeometric Analysis

K 2×2 4×4 8×8 16 × 16

h = 1/16 κ2 it. 1.19 5

h = 1/32 κ2 it. 1.35 6 1.62 8

h = 1/64 κ2 it. 1.55 6 2.19 9 1.77 8

h = 1/128 κ2 it. 1.78 6 2.86 10 2.55 10 1.87 8

Table 5.2 BDDC deluxe preconditioner for a 2D quarter-ring domain: condition number κ2 and iteration counts it. as a function of the number of subdomains K and mesh size h. Fixed p = 5, k = 4.

p κ2 it.

2 3.22 10

3 2.68 10

4 2.41 9

5 2.19 9

6 2.04 9

7 1.91 8

8 1.80 8

9 1.72 8

10 1.62 9

Table 5.3 BDDC deluxe dependence on p in the 2D quarter-ring domain: condition number κ2 and iteration counts it. as a function of the spline degree p. Fixed h = 1/64, K = 4 × 4, k = p − 1.

5.2. Dependence on p in 2D. In this test, we compare the BDDC deluxe performance with respect to the polynomial degree p = q and the regularity k. We recall that our theoretical results in Sec. ?? are only an h-analysis and do not cover the convergence rate dependence on p and k. The domain considered is the quarter-ring of Fig. ?? a) discretized with a mesh size h = 1/64 and K = 4 × 4 subdomains, while the degree p varies from 2 to 10 and the regularity k = p − 1 is maximal everywhere. The results in Table ?? show that the condition numbers and iteration counts are bounded independently of the degree p and actually improve slightly for increasing p. This remarkable results are a considerable improvement over the results in p with ρ and stiffness scaling of our previous work as in Tables 4, 5 of [?], especially taking into account that the latter were obtained only for low regularity at the subdomain interface kΓ = 0, 1, while here the regularity k = p − 1 is maximal everywhere. 5.3. Robustness with respect to discontinuous coefficients in 2D. We next investigate the robustness of the BDDC deluxe preconditioner with respect to jump discontinuities of the coefficient ρ of the elliptic problem (??). We consider three different classical tests, that we call ”central jump”, ”random mix” and ”checkerboard” in a 2D quarter-ring decomposed into 4 × 4 subdomains, see Fig. ??. In the central jump test, the coefficient ρ varies of 8 orders of magnitude (from 10−4 to 104 ) in the 2 × 2 central subdomains, while it equals 1 in the surrounding subdomains. In the random mix test, ρ has random values varying by 8 orders of magnitude in the different subdomains, see Fig. ?? (upper part) for a depiction of ρ in parametric space. In the checkerboard test, ρ = 104 in the black subdomains, while ρ = 1 in the white subdomains. We fix h = 1/64, H/h = 16, p = 3, k = 2, and the regularity at the subdomain interfaces kΓ = 0, since in the presence of jumps in the elliptic coefficient this is the correct choice for approximation reasons, see e.g. [?]. The condition numbers and iteration counts reported in the tables of Fig. ?? clearly show the robustness of BDDC deluxe, since κ2 and it. are essentially independent of the jumps in ρ. 5.4. Performance of the reduced primal set VbΠC4 in 2D. We study the optimality, scalabilty and dependance on growing p of the BDDC deluxe preconditioner

16

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

central 1 1 1 ρ 1 ρ 1 1

jump 1 1 ρ 1 ρ 1 1 1

10−3 101 10−2 100

central jump κ2 it. ρ 10−4 5.41 13 10−2 5.42 13 1 5.94 13 7.18 13 102 104 7.26 13

random mix 102 10−4 −1 10 100 3 10 102 4 10 10−3

102 104 10−4 101

random mix κ2 it. 5.25 11 checkerboard κ2 it. 5.11 14

Fig. 5.2. BDDC deluxe robustness with respect to jump discontinuities in the coefficient ρ of the elliptic problem. Central jump, random mix and checkerboard tests. Condition number κ2 and CG iteration counts it. Fixed h = 1/64, K = 4 × 4, H/h = 16.

OPTIMALITY p = 3, N = 2 × 2 H/h κ2 it. 8 3.67 8 3.67 8 16 32 3.66 9 64 4.19 9

SCALABILITY p = 3, H/h = 8 N κ2 it. 2×2 3.67 8 4×4 43.49 16 8×8 46.02 27 16 × 16 46.58 29

DEPENDANCE ON p 1/h = 64, N = 4 × 4 p κ2 it. 2 16.03 15 3 56.70 14 4 967.40 20 5 11415.81 22

Table 5.4 Optimality, scalability and dependance on p of the BDDC deluxe preconditioner with reduced primal set VbΠC4 in the 2D quarter of ring domain: condition number κ2 and CG iteration counts it. The spline regularity k is always maximal, i.e. k = p − 1.

with the reduced primal set VbΠC4 in the 2D quarter of ring domain. The spline regularity k is kept maximal, i.e. k = p − 1. The results reported in Table ?? show that, according to the theoretical estimates, the BDDC deluxe preconditioner is optimal and scalable also with the reduced primal set VbΠC4 , but the behaviour with respect to increasing p is much worse than the with the rich primal set VbΠC . 5.5. Scalability in K in 3D. We now perform 3D parallel tests on a BlueGene/Q machine, assigning one subdomain per processor. We start with weak scalability tests on the unit cube: we increase the number of subdomains K = 23 , · · · , 103 , keeping the local size H/h = 8 fixed, polynomial degree p = 3 and maximum regularity k = 2. We consider three primal spaces of increasing sizes consisting of subdomain vertices (VbΠC ), or subdomain vertices and edge averages (VbΠC+EDeli ), or subdomain vertices and edge and face averages (VbΠC+E+F ). The results reported in Table ?? show that BDDC with both deluxe scaling (top) and stiffness scaling (bottom) are scalable, since the condition numbers and iteration counts are bounded from above by constants independent of the number of subdomains. As expected, larger primal spaces yield faster convergence rates. Deluxe scaling performs better than stiffness scaling in all tests, requiring about half the iteration counts. We remark that BDDC deluxe yields quite small condition numbers, bounded by 9 for the primal space C, by 2 for VbΠC+E and by 1.4 for VbΠC+E+F .

17

Deluxe BDDC Methods for Isogeometric Analysis

1 0.8 0.6 0.4 0.2 0 −0.2 2 1.5 1 1.5

0.5 1 0

0.5 0

Fig. 5.3. 3D domains used in the numerical tests.

63 73 Deluxe scaling κ2 8.96 8.38 8.44 8.38 8.35 8.35 it. 20 21 23 24 23 23 κ2 2.06 2.01 1.98 1.98 1.98 1.98 it. 10 11 11 10 10 10 1.40 1.40 κ2 1.42 1.40 1.41 1.40 it. 8 8 8 8 8 8 Stiffness scaling κ2 20.09 19.24 19.16 19.16 19.16 19.16 it. 26 33 38 39 39 39 κ2 6.04 6.08 6.08 6.10 6.09 6.10 it. 21 22 22 22 22 23 κ2 6.04 6.08 6.08 6.10 6.09 6.10 it. 21 22 22 22 22 23 K

VbΠC VbΠC+E VbΠC+E+F

VbΠC VbΠC+E VbΠC+E+F

23

33

43

53

83

93

103

8.35 24 1.98 10 1.40 8

8.36 24 1.98 10 1.40 8

8.35 24 1.98 10 1.40 8

19.16 39 6.09 22 6.09 22

19.16 39 6.10 23 6.10 23

19.17 39 6.10 22 6.10 22

Table 5.5 BDDC weak scalability on the unit cube with different coarse spaces. Deluxe scaling (top), stiffness scaling (bottom). Condition number κ2 and iteration counts it. as a function of the number of subdomains K = n × n × n; fixed k = 2, p = 3, H/h = 8.

5.6. Quasi-optimality in H/h in 3D. Table ?? reports the results of parallel tests with increasing values of H/h = 4, 8, · · · , 24 for fixed p = 3, k = 2, 4 × 4 × 4 subdomains. Since the domain and its subdivision are fixed (H = 1/4), here we are varying the mesh size h. We consider again the primal spaces VbΠC , VbΠC+E , and VbΠC+E+F and deluxe scaling (top), stiffness scaling (bottom). The results show a linear dependence of the condition number on H/h for the primal space C, both for stiffness and deluxe scaling. The addition of edge averages (VbΠC+E ) seems sufficient to obtain a quasi-optimal method. Again deluxe scaling requires half (or less) the iterations of stiffness scaling. The void columns could not be run due to memory limitations. 5.7. Dependence on p in 3D. In this test, we study the BDDC convergence rate for increasing polynomial degree p = 2, 3 · · · , 7 and maximal regularity k = p − 1. The domain considered is the unit cube, with mesh size h = 1/32, subdivided into K = 4 × 4 × 4 subdomains. The results reported in Table ?? show that the condition

18

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini

H/h VbΠC VbΠC+E VbΠC+E+F

VbΠC VbΠC+E VbΠC+E+F

4

8

κ2 it. κ2 it. κ2 it.

2.62 12 1.54 9 1.54 9

κ2 it. κ2 it. κ2 it.

7.03 24 6.35 23 6.35 23

12 16 Deluxe scaling 6.13 10.10 14.19 18 21 24 1.80 2.03 2.21 10 11 12 1.37 1.46 1.62 8 8 9 Stiffness scaling 10.59 21.30 34.64 26 34 40 6.09 6.13 6.26 22 22 23 6.09 6.13 6.16 22 22 21

20

24

17.91 27 2.35 12 1.75 9

— — — — — —

46.97 46 8.15 26 6.19 22

66.86 54 10.28 29 6.21 22

Table 5.6 BDDC dependence on H/h, unit cube. Condition number κ2 and PCG iteration counts it. Fixed p = 3, k = 2, K = 4 × 4 × 4.

VbΠC VbΠC+E VbΠC+E+F

VbΠC VbΠC+E VbΠC+E+F

p

2

3

κ2 it. κ2 it. κ2 it.

14.81 31 1.43 12 1.43 8

8.42 23 1.41 11 1.41 8

κ2 it. κ2 it. κ2 it.

29.91 37 5.43 19 1.81 10

19.18 38 6.09 22 6.08 22

4 5 6 Deluxe scaling 6.73 5.10 6.28 21 18 21 1.76 2.87 6.03 10 14 20 1.76 2.87 6.03 10 13 20 Stiffness scaling 38.42 469.73 4705.34 61 184 588 33.02 440.48 4258.72 55 180 547 33.12 440.19 4257.88 56 180 554

7 53.15 58 50.29 57 50.26 55 97419.31 >1000 89459.12 >1000 89229.54 >1000

Table 5.7 BDDC dependence on p for the 3D unit cube with different coarse spaces. Deluxe scaling (top), stiffness scaling (bottom). Condition number κ2 and iteration counts it. as a function of the spline degree p. Fixed h = 1/32, K = 4 × 4 × 4, k = p − 1.

numbers and iteration counts using deluxe scaling (top table) are almost independent of p for p ≤ 5; when p > 5, the convergence rate of deluxe BDDC begin to degrade but is still orders of magnitude better than the convergence rate of BDDC with stiffness scaling (bottom table). Indeed, the condition numbers with stiffness scaling are already high with p > 3 and rapidly degenerate for higher p. Interestingly, the condition number of the vertex only (C) coarse space become closer to the other coarse spaces when using deluxe scaling with higher polynomial degrees. The addition of faces averages (F) to the primal space does not improve the condition number in

19

Deluxe BDDC Methods for Isogeometric Analysis

VbΠC

VbΠC+E

VbΠC+E+F

ρ 10−4 10−2 1 102 104 10−4 10−2 1 102 104 10−4 10−2 1 102 104

central jump κ2 it. 117.37 44 118.40 44 134.04 48 137.15 50 137.40 52 5.33 18 5.33 18 5.27 18 4.92 16 4.88 16 1.98 10 1.99 10 2.05 10 2.05 10 2.05 10

checkerboard κ2 it. — — — — 134.04 48 102.11 43 104.31 44 — — — — 5.27 18 4.19 16 4.20 16 1.98 10 1.99 10 2.05 10 2.05 10 2.05 10

random κ2 — — 134.04 126.53 123.63 — — 5.27 4.83 4.87 — — 2.05 2.00 2.00

mix it. — — 48 47 46 — — 18 16 16 — — 10 10 9

Table 5.8 BDDC deluxe robustness with respect to jump discontinuities in the coefficient ρ of the elliptic problem. Central jump, random mix and checkerboard tests. Condition number κ2 and PCG iteration counts it. Fixed h = 1/32, p = 3, C 0 continuity at the interface, 4 × 4 × 4 subdomains.

any of the cases considered (except for p = 2 and stiffness scaling). The number of subproblems which must be solved in the deluxe case is 252; the adjacency graph of such problems can be colored with 18 colors. The total solving time using deluxe scaling becomes lower than the solving time with stiffness scaling starting from p = 6. 5.8. Robustness with respect to discontinuous coefficients in 3D. Finally, we report the results of parallel tests with discontinuous coefficient ρ for a unit cube discretized with H/h = 8, p = 3, kΓ = 0 continuity at the interface, 4×4×4 subdomains. The distribution of values of ρ is the 3D analog of the 2D case, i.e. ”central jump”, ”checkerboard” and ”random mix”. The random setting is somewhat different from the 2D case; given a value of ρ, every subdomain generates a random number r in the interval (− log10 ρ, log10 ρ) and the scaling factor for the subdomain is then chosen as 10r . In the checkerboard case, the local matrices are scaled by ρ on black subdomains and by 1/ρ on white subdomains. Note that the rows associated with a value of ρ < 1 for the checkerboard and random test cases are void since these cases are already covered by the rows corresponding to the reciprocal of the ρ value. The results clearly show the BDDC robustness wit respect to jumps in ρ. The convergence rate with the primal space C is somewhat unsatisfactory, but VbΠC+E and especially VbΠC+E+F yield remarkably small condition numbers and low iteration counts. REFERENCES [1] S. Balay, J. Brown, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. Curfman McInnes, B. F. Smith, H. Zhang. PETSc Web page, http://www.mcs.anl.gov/petsc (2012). [2] Y. Bazilevs, L. Beir˜ ao da Veiga, J.A. Cottrell, T.J.R. Hughes, G. Sangalli. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. Math. Mod. Meth. Appl. Sci., 16, 1–60, 2006. [3] L. Beir˜ ao da Veiga, C. Chinosi, C. Lovadina and L. F. Pavarino. Robust BDDC preconditioners

20

[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

[16]

[17] [18]

[19] [20] [21] [22] [23] [24] [25] [26] [27]

[28] [29] [30] [31]

L. Beir˜ ao da Veiga, L. F. Pavarino, S. Scacchi, O. B. Widlund, S. Zampini for Reissner-Mindlin plate bending problems and MITC elements. SIAM J. Numer. Anal., 47 (6): 4214–4238, 2010. L. Beir˜ ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. Overlapping Schwarz methods for Isogeometric Analysis. SIAM J. Numer. Anal., 50 (3), 1394-1416, 2012. L. Beir˜ ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. BDDC preconditioners for Isogeometric Analysis. Math. Mod. Meth. Appl. Sci. 23 (6): 1099–1142, 2013. L. Beir˜ ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. Isogeometric Schwarz preconditioners for linear elasticity systems. Comp. Meth. Appl. Mech. Engrg., 253: 439–454, 2013. L. Beir˜ ao da Veiga, D. Cho and G. Sangalli. Anisotropic NURBS approximation in isogeometric analysis. Comput. Meth. Appl. Mech. Engrg., 209-212: 1–11, 2012. M. Bercovier, I. Soloveichik. Additive Schwarz Decomposition methods applied to isogeometric analysis. Submitted. S. C. Brenner and L.-Y. Sung. BDDC and FETI-DP without matrices or vectors, Comput. Methods Appl. Mech. Engrg., 196 (2007), pp. 1429–1435. J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs. Isogeometric Analysis. Towards integration of CAD and FEA. Wiley, 2009. C. De Falco, A. Reali, R. Vazquez. GeoPDEs: a research tool for Isogeometric Analysis of PDEs. Advances in Engineering Software, 42 (12), 1020-1034, 2011. M. Dryja, J. Galvis, and M. Sarkis, BDDC methods for discontinuous Galerkin discretization of elliptic problems, J. Complexity, 23 (2007), pp. 715–739. C.R. Dohrmann. A Preconditioner for Substructuring Based on Constrained Energy Minimization. SIAM J. Sci. Comput., 25: 246–258, 2003. C. R. Dohrmann, An approximate BDDC preconditioner, Numer. Linear Algebra Appl., 14 (2007), pp. 149–168. C.R. Dohrmann and O.B. Widlund. Some Recent Tools and a BDDC Algorithm for 3D Problems in H(curl). Proceedings of the 20th International Conference on Domain Decomposition Methods. Held in San Diego, CA, February 7-11, 2011. Springer LNCSE. To appear. C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, and D. Rixen, FETI-DP: a dual-primal unified FETI method – part I. A faster alternative to the two-level FETI method, Internat. J. Numer. Meth. Engrg., 50 (2001), pp. 1523–1544. K. Gahalaut, J. Kraus, S. Tomar. Multigrid Methods for Isogeometric Discretization. Comput. Methods Appl. Mech. Engrg., in press. doi: 10.1016/j.cma.2012.08.015. T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement. Comp. Meth. Appl. Mech. Engrg., 194, 4135–4195, 2005. H. H. Kim, M. Dryja, and O. B. Widlund, A BDDC method for mortar discretizations using a transformation of basis, SIAM J. Numer. Anal., 47 (2008/09), pp. 136–157. H. H. Kim and X. Tu, A three-level BDDC algorithm for mortar discretizations, SIAM J. Numer. Anal., 47 (2009), pp. 1576–1600. S. K. Kleiss, C. Pechstein, B. J¨ uttler, S. Tomar. IETI - Isogeometric Tearing and Interconnecting, Comput. Methods Appl. Mech. Engrg., 247248: 201215, 2012. A. Klawonn and O. Widlund. Dual-primal FETI methods for linear elasticity. Comm. Pure Appl. Math., 59: 1523–1572, 2006. J. H. Lee. Domain Decomposition Methods for Reissner-Mindlin Plates Discretized with the Falk-Tu Elements. Ph.D. Thesis, Courant Institute TR2011-937, 2011. J. H. Lee, Balancing Domain Decomposition by Constraints for Numerically Thin ReissnerMindlin Plates Approximated by Falk–Tu Elements. In preparation. J. Li and O. B. Widlund. FETI-DP, BDDC, and block Cholesky methods. Int. J. Numer. Meth. Engrg., 66 (2): 250–271, 2006. J. Li and O. B. Widlund, On the use of inexact subdomain solvers for BDDC algorithms, Comput. Methods Appl. Mech. Engrg., 196 (2007), pp. 1415–1428. J. Li and X. Tu, Convergence analysis of a balancing domain decomposition method for solving a class of indefinite linear systems, Numer. Linear Algebra Appl., 16 (2009), pp. 745–773. J. Mandel and C. R. Dohrmann. Convergence of a balancing domain decomposition by constraints and energy minimization. Numer. Lin. Alg. Appl., 10 (7): 639–659, 2003. J. Mandel, C. R. Dohrmann and R. Tezaur. An algebraic theory for primal and dual substructuring methods by constraints. Appl. Numer. Math., 54 (2): 167–193, 2005. J. Mandel, B. Soused´ık, and C. R. Dohrmann, Multispace and multilevel BDDC, Computing, 83 (2008), pp. 55–85. D.-S. Oh, O. B. Widlund and C. Dohrmann. A BDDC algorithm for Raviart-Thomas vector fields. TR2013-951, Courant Institute, NYU, 2013.

Deluxe BDDC Methods for Isogeometric Analysis

21

[32] L. F. Pavarino, O. B. Widlund and S. Zampini. BDDC Preconditioners for Spectral Element Discretizations of Almost Incompressible Elasticity in Three Dimensions. SIAM J. Sci. Comp., 32 (6): 3604–3626, 2010. [33] A. Toselli, O. B. Widlund. Domain Decomposition Methods: Algorithms and Theory. Computational Mathematics, Vol. 34. Springer-Verlag, Berlin, 2004. [34] X. Tu. Three-level BDDC in three dimensions. SIAM J. Sci. Comp., 29 (4): 1759-1780, 2007. [35] X. Tu and J. Li, A balancing domain decomposition method by constraints for advectiondiffusion problems, Commun. Appl. Math. Comput. Sci. 3 (2008), pp. 2560. [36] O. B. Widlund, Iterative substructuring methods: algorithms and theory for elliptic problems in the plane. In First Inter. Symp. on Domain Decomposition Methods for Partial Differential Equations, R. Glowinski et al. Eds., pp. 113–128, SIAM, 1988.