Schwarz methods for a preconditioned WOPSIP ... - Semantic Scholar

Schwarz methods for a preconditioned WOPSIP method for elliptic problems Paola F. Antonietti∗ MOX–Laboratory for Modeling and Scientific Computing, Dipartimento di Matematica Politecnico di Milano, Piazza Leondardo da Vinci 32, 20133 Milano, Italy

Blanca Ayuso de Dios† Centre de Recerca Matem´ atica Campus de Bellaterra, 08193 Bellaterra (Barcelona), Spain

Susanne C. Brenner‡ Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA

Li-yeng Sung§ Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA

June 28, 2012 Abstract We propose and analyze several two-level non-overlapping Schwarz methods for a preconditioned weakly over-penalized symmetric interior penalty (WOPSIP) discretization of a second order boundary value problem. We show that the preconditioners are scalable and that the condition number of the resulting preconditioned linear systems of equations is independent of the penalty parameter and is of order Hh−1 , where H and h represent the mesh sizes of the coarse and fine partitions, respectively. Numerical experiments that illustrate the performance of the proposed two-level Schwarz methods are also presented.

1

Introduction

The weakly over-penalized symmetric interior penalty (WOPSIP) method is a Discontinuous Finite Element method that was first introduced in [10] for the approximation of second order elliptic problems, such as the model Poisson problem of finding u ∈ H01 (Ω) such that ∫ ∫ ∇v · ∇v dx = f v dx ∀ v ∈ H01 (Ω), (1) Ω

Ω

where Ω ⊂ R2 is a polygonal domain and f a given source term in L2 (Ω). The WOPSIP method is an Interior Penalty method that preserves the symmetric positive-definiteness of the continuous problem ∗ [email protected] † [email protected]

‡ [email protected] § [email protected]

1

without requiring any tuning of a penalty parameter. Furthermore, although the method is inconsistent (in the sense of [4]) the weak over-penalization of the jumps allows for achieving optimal convergence in the energy and L2 -norms. Moreover, since the penalty term only involves the averages of the jumps across the edges, the performance of the WOPSIP method is not sensitive to regular hanging nodes. Another nice feature of the method is its simplicity and the structure of the matrix in the algebraic formulation of the method, which makes it an attractive option for practical computations. However, due to the over-penalization, the condition number of the stiffness matrix is of order O(h−4 ), h being the mesh-size. To remedy the very ill-conditioning of the WOPSIP method, a simple blockdiagonal preconditioner was introduced in [10] that reduces the condition number of the preconditioned system to O(h−2 ), which is typical for discretizations of second order elliptic problems. A particular construction of this preconditioner was later shown in [8] to retain the intrinsic built-in parallelism of the WOPSIP method, rendering the resulting preconditioned WOPSIP an ideal method for parallel computations. The goal of this paper is to construct and analyze non-overlapping Schwarz preconditioners for the (already) preconditioned WOPSIP discretization of a second order elliptic problem. Research in Domain Decomposition preconditioners for DG methods, in particular for approximating second order elliptic problems, has gained considerable attention in recent years [2, 3, 16, 17, 18, 19]. Non-overlapping and overlapping Schwarz methods were first introduced in [19] for the classical Interior Penalty (IP) scheme. These overlapping preconditioners were extended in [20] for the IP approximations of convection-diffusion problems and in [13] to the C 0 -IP approximation of a fourth order problem. A slightly different construction of non-overlapping Schwarz methods was introduced in a unified framework in [2], for all the stable and consistent DG methods in [4]. These preconditioners differ from those proposed in [19] (and all the other work in the subject) in the design of the local solvers: inexact in [2] versus exact in [19]. Both solvers have been shown to perform similarly, although the analysis in the former case turns out to be more involved. In this paper, we propose some non-overlapping preconditioners within the Schwarz framework for the preconditioned (first order) WOPSIP method. The construction with both exact and inexact local solvers is considered and we allow for many possible choices of the coarse solver. We perform the convergence analysis of the proposed Schwarz preconditioners, following the classical framework [21]. They are shown to be scalable and the condition number of the resulting preconditioned systems is shown to be of order O(Hh−1 ) (H denoting the mesh size of the coarse partition), similar to other Schwarz methods for DG discretizations of second order problems. It is worth mentioning that for the original WOPSIP method (without preconditioner), overlapping Schwarz methods have been previously considered in [6]. An important difference with these Schwarz methods, is that ours require the solution of local problems whose conditioning is O(h−2 ) in contrast to the much more ill-conditioned O(h−4 ) considered in [6]. In our convergence analysis, however, an additional difficulty arises since the matrix for the preconditioned WOPSIP method in previous work [8, 10] is non-symmetric. A first step (and contribution) in our analysis, is to identify a new DG method, P-WOPSIP, by characterizing the variational formulation that gives rise to an operator that is spectrally equivalent to the non-symmetric matrix associated with the preconditioned WOPSIP. We analyze the spectral properties of the P-WOPSIP method, showing that it is indeed spectrally equivalent to all the classical consistent and stable DG methods for second order problems. Its convergence and approximation properties can then be studied following the classical framework [4]. Furthermore, special care is taken to ensure that the resulting symmetric DG method (P-WOPSIP), inherits the simplicity and special structure of the original preconditioned WOPSIP [8, 10]. A particular attractive feature of the P-WOPSIP method is that its spectral properties are completely independent of the penalty parameter α. As a result, the performance of the Schwarz preconditioners constructed here is also independent of the penalty parameter. Finally, we mention that the WOPSIP method and the results in this paper can be extended to three

2

dimensional problems in a straightforward manner. High order WOPSIP methods have been recently introduced in [11]. The extension of the results presented here to the higher order methods will be a subject of future research. The rest of the paper is organized as follows. In the next section, we review the WOPSIP method, its block-diagonal preconditioner and characterize the bilinear form defining the P-WOPSIP method. In Section 3 we introduce the Schwarz methods for the P-WOPSIP discretization and discuss some computational issues. The convergence analysis is carried out in Section 4. Numerical experiments validating the theory are presented Section 5. Finally, the proofs of some technical results needed in our theoretical analysis are given in the Appendices. Throughout the paper, we shall use standard notation for Sobolev spaces (cf. [1]), and x . y will mean that there exists a generic constant C > 0 (that may not be the same at different occurrences but is always independent of the mesh and penalty parameters) so that x ≤ C y. Analogously, x ≈ y will mean that C −1 y ≤ x ≤ C y, for a constant C > 0.

2

Problem setting and WOPSIP discretization

In this section, we introduce some notation, recall the WOPSIP approximation and present some of the properties of the formulation. Let {Th }h>0 be a family of quasi-uniform triangulations of Ω. The mesh size is defined by h := maxT ∈Th diam T . We denote by Vh the first order discontinuous finite element space associated with Th , defined by Vh := {v ∈ L2 (Ω) : v|T ∈ P1 (T ) ∀ T ∈ Th } , where P1 (T ) is the space of linear polynomials in T . The set of all the edges in Th is denoted by Eh ; the set of internal edges by Eh◦ and the set of boundary edges by Eh∂ , so that Eh := Eh◦ ∪ Eh∂ . For any e ∈ Eh , he will denote the length of the edge e. We use standard notation for trace operators [4] to define the jumps [[v]], [[τ ]] and averages {{v}}, {{τ }} of (sufficiently regular) scalar and vector–valued functions v and τ . For each interior edge e ∈ Eh◦ such that e = ∂T + ∩ ∂T − we define − − + + − − [[v]] := ve+ n+ e + ve ne , [[τ ]] := τ e · ne + τ e · ne ,

{{v}} := (ve+ + ve− )/2,

− {{τ }} := (τ + e + τ e )/2,

where ve+ (respectively ve− ) denotes the trace of v on e taken within the interior of T + (respectively − + T − ), and n+ (respectively e (respectively ne ) is the unit normal of e pointing towards the outside of T − ∂ T ). For e ∈ Eh , we define {{τ }} := τ e , [[v]] := ve n. We do not need either [[τ ]] or {{v}} on boundary edges, and we leave them undefined. The WOPSIP approximation to the solution of the model problem (1) reads: Find uh ∈ Vh such that ∫ Ah (uh , v) = f v dx ∀ v ∈ Vh , (2) Ω

where Ah (·, ·) : Vh × Vh −→ R is the bilinear form defined in [10, 8]: ∑ α ∫ ∑ ∫ ∇w · ∇v dx + Ah (w, v) := Π0 ([[w]]) · Π0e ([[v]]) ds h3e e e T T ∈Th

e∈Eh

3

∀ w, v ∈ Vh ,

(3)

Here, α denotes the penalty parameter which we assume to be ≥ 1 and Π0e : L2 (e) −→ P0 (e) is the L2 -orthogonal projection onto the space P0 (e) of constant functions on e: ∫ 1 Π0e (v) := v ds = v(me ) ∀ e ∈ Eh ∀v ∈ Vh , (4) he e where in the last step we have used the midpoint rule for integration and me is the midpoint of the edge e ∈ Eh . For vector valued functions Π0e (·) is defined componentwise. By considering the energy norm: ∥v∥2h :=

∑ T ∈Th

∥∇v∥20,T +

∑ 1 ∥Π0 ([[v]])∥20,e h3e e

∀ v ∈ Vh ,

e∈Eh

(observe that ∥v∥2h = Ah (v, v) for α = 1), it can be shown that for α > 0, the bilinear form defining the WOPSIP method is coercive and continuous in Vh : Ah (v, v) ≥ ∥v∥2h Ah (v, w) ≤ max (1, α)∥v∥h ∥w∥h

∀v ∈ Vh , ∀v, w ∈ Vh .

Also, optimal rates of convergence in the ∥ · ∥h and L2 -norms can be proved for the WOPSIP approximation to problem (1) (i.e., the solution of (2)). For details see [9, 10].

2.1

An efficient preconditioner for the WOPSIP method

We recall that, given a basis of Vh , any function v ∈ Vh is uniquely determined by a set of degrees of freedom (dofs). If Ah is the stiffness matrix associated with the bilinear form Ah (·, ·) and the given basis, problem (2) can be rewritten as the linear system of equations Ah u = f, with Ah symmetric and positive definite. Due to the over-penalization of the method, it can be easily seen that the condition number of Ah is of order κ(Ah ) = O(h−4 ). To effectively compute the approximation with the WOPSIP method, the bilinear form ∑ α ∫ ∑ ∑ wT (me )vT (me ) + Bh (w, v) := Π0 ([[w]]) · Π0e ([[v]]) ds ∀w, v ∈ Vh , h3e e e T ∈Th e⊂∂T

e∈Eh

was introduced in [10], where wT := w|T for all T ∈ Th . Denoting by Bh the matrix associated with the bilinear form Bh and the given basis, it was proved in [10] that vT Bh v . vT Ah v . h−2 vT Bh v ∀ v ∈ Rn ,

(5)

where n := dim Vh . From (5), it immediately follows that −2 κ(B−1 ). h Ah ) = O(h

The issue of the efficiency of the preconditioner Bh was further explored in [8], where the authors showed that if a suitable ordering of the dofs is employed the resulting matrix Bh (and so its action) turns out to be block-diagonal with 1 × 1 and 2 × 2 blocks and therefore can be computed in parallel. The aim of this paper is to design a Schwarz method for the efficient solution of the linear system −1 B−1 h Ah u = Bh f .

4

6

7

15

9

16

4

5

23

22

17

30

18 31

33

28

40

29

47

46 41

2

42

3

6

18 17

19

20

30

8

29

31 32

38

11

37

39

40

48 47

8

24

32

48

21

26

41

44

21

11

45

35

22

25

42

43

14 19

20

2

13

1

12

10 39

38

3

43

37

44

26

25 36

34

5

27

27 28

16

15

4

13

14

1

23 24

45

46

10

36

35

33 34

12

9

7

Figure 1: Sample of elementwise (left) and edgewise (right) ordering of the degrees of freedom. Note that, although Ah is a symmetric and positive definite (s.p.d.) matrix, B−1 h Ah is no longer symmetric in general. Hence, to avoid the non-symmetry and the resulting difficulties, we consider the equivalent linear system of equations −1/2 Dh y = Bh f, (6) 1/2

where y := Bh u and

−1/2

Dh := Bh

−1/2

Ah Bh

, 1/2

which is well-defined since Bh is s.p.d. and therefore it admits a unique s.p.d. square root Bh . From now on, we focus on the construction of Schwarz preconditioners for the s.p.d. system (6). Clearly, it still holds that −1/2 −1/2 κ(Bh Ah Bh ) = O(h−2 ), −1/2

since we can take v = Bh

w in (5) for any w ∈ Rn .

So far, we have not said anything about the selection of the basis or the location of the dofs of Vh . In [8], it was shown that the use of the Crouziex-Raviart basis for P1 (T ) on each T ∈ Th and the choice of the dofs at the midpoints of the edges in each T have some advantages. More precisely, the authors showed that by using an edgewise ordering of dofs (that is, the dofs associated to the midpoints of an interior edge are always consecutive, cf. Figure 1 for an example), the matrix Bh , and consequently B−1 h , turns out to be block-diagonal with 1 × 1 and 2 × 2 blocks, and therefore the preconditioned WOPSIP method has an intrinsic highly parallel structure. In the next section we show that by using −1/2 the same special ordering, also the action of Bh retains the same highly parallel structure and can be efficiently computed. Moreover, we shall also show that this ordering facilitates our analysis of the Schwarz methods for the preconditioned WOPSIP discretization. Hence, throughout the rest of the paper it is assumed that the edgewise ordering is employed (see Section 3.2 for details on the implementation).

2.2

−1/2

Construction of Bh

As shown in [8], by ordering the dofs in an edgewise manner (cf. Figure 1 (right)) the matrix representing Bh is block-diagonal, with either 2 × 2 blocks (corresponding to an interior edge) or 1 × 1 blocks (corresponding to a boundary edge). Denoting by Beh the block of the matrix Bh corresponding

5

to the dofs associated to the edge e ∈ Eh , we have ]  [ 1 1 + θ −1  e  if e ∈ Eh◦ ,  −1 1 + θe θe e Bh = ]  1 [   1 + θe if e ∈ Eh∂ , θe where θe = h2e /α for all edges e ∈ Eh . Observe that for any e ∈ Eh◦ , since Beh is s.p.d. it can be diagonalized as follows: ] [ ] [ ] [ 1 0 1 1 1 1 1 1 Beh = QΛQT = · . · e 0 2+θ 1 −1 2 2 1 −1 θe And so, we obtain an explicit expression for (Beh )−1/2 , ] [ [ ] [ 1 0 1 1 1 1 1 √ e −1/2 −1/2 T (Bh ) · = QΛ Q = · θe 1 0 2 2 1 −1 2+θe √

where we have set βe :=

θe 2 + θe

1 −1

]

1 = 2

[

1 + βe 1 − βe

1 − βe 1 + βe

] , (7)

∀e ∈ Eh◦ .

For e ∈ Eh∂ , we simply have (Beh )−1/2

=

[

βe∂

]

√ ,

βe∂

:=

θe 1 + θe

∀e ∈ Eh∂ .

Let β be the piecewise constant function on Eh defined by √  θe   if e ∈ Eh◦ ,  βe := h2 2 + θe β|e := θe := e ∀e ∈ Eh . √  α θe   β ∂ := if e ∈ Eh∂ , e 1 + θe √ Observe that β|e ≈ h/ α for all e ∈ Eh and ( )  α 1 1   if e ∈ Eh◦ , ≤ h 2 2α + h 2h α e e e ) ( (β|e )2 3 =  he α 1 1   ≤ if e ∈ Eh∂ . he α + h2e he

(8)

(9)

Having found an explicit expression for each block (Beh )−1/2 , we look at its action on the vector of degrees of freedom associated to an edge e ∈ Eh . Let e ∈ Eh◦ be an arbitrary edge shared by the elements T + and T − , e = T + ∩ T − , and let ue := [u+ , u− ]T denote the nodal values of the trace u|e of u at the midpoint of the edge e. It follows from (7) that ] [ n+ e · [[u]] { {u} } + β e 2 . (Beh )−1/2 ue = (10) n+ {{u}} − βe 2e · [[u]] If e ∈ Eh∂ is a boundary edge, we have (Beh )−1/2 ue =

[ √

θe 1+θe

6

]

ue = βe∂ ue .

(11)

−1/2

We can then define Bh −1/2 Bh .

: Vh −→ Vh to be the operator represented by the block-diagonal matrix

Finally, we introduce the bilinear form Dh (·, ·) : Vh × Vh −→ R defined by −1/2

Dh (u, v) := Ah (Bh

−1/2

u, Bh

v) =

∑ ∫

−1/2

T

T ∈Th

∇(Bh

−1/2

u) · ∇(Bh

v) dx

∑ α ∫ −1/2 −1/2 + Π0 ([[Bh u]]) · Π0e ([[Bh v]]) ds , h3e e e

(12)

e∈Eh

and the norm ∥u∥2DG :=

∑

∥∇u∥20,T +

T ∈Th

∑ 1 ∥Π0 ([[u]])∥20,e he e

∀ u ∈ Vh .

(13)

e∈Eh

The next result shows that Dh (·, ·) is continuous and coercive in Vh with respect to the above DG norm, provided h is small enough (see Remark 2.2). Lemma 2.1. The bilinear form Dh (·, ·) defined by (12) is continuous in the DG norm (13), and it is also coercive for all h ≤ h0 with √ α h0 := , (14) 16Ct2 + 1 where Ct is a trace inequality constant. More precisely, there exist constants Cc , Cs > 0 depending only on the shape regularity of the partition but independent of the penalty parameter α and the mesh size h, such that Continuity: Coercivity:

Dh (u, w) ≤ Cc ∥u∥DG ∥w∥DG Dh (u, u) ≥

Cs ∥u∥2DG

∀ u, w ∈ Vh ;

(15)

∀ u ∈ Vh .

(16)

The proof of Lemma 2.1 can be found in Appendix A. We also define the following bilinear forms ∑ ∫

∑ 1 ∫ [[w]] · [[v]] ds, he e T ∈Th T e∈Eh ∑ ∫ ∑ 1 ∫ Sh∗ (·, ·) : Vh × Vh −→ R, Sh∗ (w, v) := ∇w · ∇v dx + Π0 ([[w]]) · Π0e ([[v]]) ds. he e e T Sh (·, ·) : Vh × Vh −→ R,

Sh (w, v) :=

∇w · ∇v dx +

T ∈Th

(17)

e∈Eh

Remark 2.2. The restriction on h in Lemma 2.1 is necessary for guaranteeing the coercivity in the DG norm ∥·∥DG . Note from the definition of h0 and taking into account our assumption α ≥ 1 together with the fact that for piecewise linear polynomials on triangles Ct2 ≈ 1/3 (see for instance [22]), the above restriction on h is a very mild one. Remark 2.3. Notice that since Dh (·, ·) is symmetric, Lemma 2.1 implies in particular that Dh (·, ·), Sh (·, ·) and Sh∗ (·, ·) are spectrally equivalent.

3

Schwarz methods for the P-WOPSIP discretization

In this section we introduce the Schwarz methods and provide some technical tools needed in the analysis. ∪N We denote by TN a partition of Ω into N non-overlapping subdomains, i.e., Ω = i=1 Ωi , and by {TH }H>0 and {Th }h>0 two families of coarse and fine partitions, respectively, with mesh sizes H > 0 7

and h > 0. All the partitions are assumed to be regular and quasi-uniform and we shall always proceed under the assumption that Th , TH and TN are nested: TN ⊆ TH ⊆ Th , i.e., each Ωi , i = 1, . . . , N , can be written as the union of some elements D ∈ TH , each of which is the union of elements of the finer partition Th ; that is ∪ D= Ti ∀D ∈ TH . Ti ∈Th Ti ⊂D

For each subdomain Ωi ∈ TN , i = 1, . . . , N , we define the local DG spaces Vhi as Vhi := {u ∈ L2 (Ωi ) : v|T ∈ P1 (T )

∀ T ∈ Th , T ⊂ Ωi },

and denote by RTi : Vhi −→ Vh the standard inclusion operator from Vhi to Vh , and by Ri its transpose with respect to the canonical bilinear form. We observe that Vh = RT1 Vh1 ⊕ . . . ⊕ RTN VhN . Finally, we define Γ :=

N ∪

Γij

Γij := {e ∈ Eh , such that e ⊂ ∂Ωi ∩ ∂Ωj , i ̸= j} .

(18)

i,j=1

We now introduce the local solvers, for which we consider two classes: exact local solvers (as those proposed in [19]) and inexact local solvers (as those introduced in [2, 3]). (i) Exact local solvers: For each subdomain Ωi ∈ TN , i = 1, . . . , N , the local bilinear form DiE (·, ·) : Vhi × Vhi −→ R is defined as the restriction of the P-WOPSIP bilinear form (12) to the space RTi Vhi : −1/2 T −1/2 Ri ui , Bh RTi vi )

DiE (ui , vi ) := Dh (RTi ui , RTi vi ) = Ah (Bh

∀ui , vi ∈ Vhi .

(19)

(ii) Inexact local solvers: Following [2], we first consider the model problem (1) set in the subdomain Ωi : ∫ ∫ ∇ui · ∇vi dx = f vi ∀vi ∈ H01 (Ωi ). (20) Ωi

Ωi

th

The i -local solver (that is associated to the subdomain Ωi ) is defined as the P-WOPSIP approximation to (20). Hence, the local bilinear form DiI (·, ·) : Vhi × Vhi −→ R is given by: −1/2

DiI (ui , vi ) := Ai (Bi where Ai (·, ·) is given by: Ai (wi , vi ) :=

∑ ∫ T ∈Th T ⊂Ωi

−1/2

ui , Bi

∇wi · ∇vi dx +

T

vi )

∀ui , vi ∈ Vhi ,

∑ α ∫ Π0 ([[wi ]]) · Π0e ([[vi ]]) ds. h3e e e

e∈Eh e⊂Ωi

(21)

(22)

and Bi : Vhi −→ (Vhi )′ refers to the operator associated to the bilinear form Bi (·, ·) defined by ∑ α ∫ ∑ ∑ wiT (me )viT (me ) + Π0 ([[wi ]]) · Π0e ([[vi ]]) ds ∀wi , vi ∈ Vhi . Bi (wi , vi ) := h3e e e T ∈Th e⊂∂T T ⊂Ωi

e∈Eh e⊂Ωi

8

Note that, the edges e ∈ Eh◦ such that e ⊂ ∂Ωi although interior edges in the global partition Th , are however boundary edges with respect to the local partitioning induced in the subdomain Ωi . On these edges the definition of the jump operator on boundary edges applies, i.e., [[wi ]] = wi |e n. −1/2 Consequently the action of Bi on the functions restricted to these edges is given by (11). A key issue in the analysis of the non-overlapping Schwarz methods is the relation between the global bilinear form Dh (·, ·) and the sum of the local solvers. To study such a relation, we need first to introduce some additional notation. Recalling the definition (18) of the interface Γ, we define the strip ΩΓ as N ∪ (23) ΩΓ = ΩΓij , ΩΓij = {T ∈ Th | T has one edge in Γij } . i,j=1

Following [19, 2] we have the following result whose proof can be found in Appendix B. Lemma 3.1. For any u ∈ Vh , let ui ∈ Vhi , i = 1, . . . , N , be the unique functions such that u = ∑N T i=1 Ri ui . Then the following identities hold: Dh (u, u) =

N ∑

DiE (ui , ui ) + IhE (u, u) ,

(24)

DiI (ui , ui ) + IhI (u, u) ,

(25)

i=1

Dh (u, u) =

N ∑ i=1

where IhE (u, u) = 2

[ ∑ ∫ T ∈ΩΓ

−1/2 T Ri ui )

T

∇(Bh

−1/2 T Rj uj ) dx

· ∇(Bh

−

∑

βe2

e∈Γ

α h3e

∫

] Π0e (ui )Π0e (uj ) ds ,

(26)

e

IhI (u, u) = IhE (u, u) + GIh (u, u),

(27)

with βe defined as in (8), and GIh (u, u)

=

N [ ∑ ∫ ∑ i=1

T ∈ΩΓ

T

−1/2 ∇(Bh RTi ui )

·

∑ ∫

−1/2 ∇(Bh RTi ui ) dx −

T ∈ΩΓ T ⊂Ωi

−

−1/2

∇(Bi

−1/2

ui ) · ∇(Bi

] ui ) dx,

T

∑ α ∫ ( ) ηe2 [Π0e (ui )]2 + [Π0e (uj )]2 ds. (28) 3 he e e∈Γ

Here, ηe is defined as: ηe2 := −[βe ]2 + [βe∂ ]2 =

θe > 0, (1 + θe )(2 + θe )

θe = h2e /α .

(29)

The last ingredient in the construction of the Schwarz methods is the coarse solver. We consider a coarse partition TH and we take for ℓ = 0, 1 VH ≡ Vh0 := {v ∈ L2 (Ω) : v|T ∈ Pℓ (T )

∀D ∈ TH },

i.e., the functions in the coarse space can be either piecewise constant or piecewise linear. We denote by RT0 : Vh0 −→ Vh the standard inclusion operator from Vh0 to Vh , by R0 its transpose with respect to the canonical bilinear forms, and define the following three coarse solvers: D0 (u0 , v0 ) := Dh (RT0 u0 , RT0 v0 )

∀u0 , v0 ∈ Vh0 , 9

(30)

S0 (u0 , v0 ) := Sh (RT0 u0 , RT0 v0 )

∀u0 , v0 ∈ Vh0 ,

(31)

S0∗ (u0 , v0 )

∀u0 , v0 ∈

(32)

:=

Sh∗ (RT0 u0 , RT0 v0 )

Vh0 ,

where Sh (·, ·), Sh∗ (·, ·) : Vh × Vh −→ R are defined in (17). Remark 3.2. As in [19, 2], the coarse solver D0 (·, ·) is defined as the restriction of the original method to the coarse finite element space Vh0 . However, it should be noted that D0 (u0 , v0 ) := Dh (RT0 u0 , RT0 v0 ) ̸= DH (u0 , v0 )

∀ u0 , v0 ∈ VH .

In particular to ensure the performance of the resulting Schwarz method it turns out to be essential to choose the penalty parameter αH in the definition of AH (·, ·) as αH = α(H/h)3 . Remark 3.3. Since the coarse solvers (30), (31) and (32) are defined as the restriction of Dh (·, ·), Sh (·, ·) and Sh∗ (·, ·), respectively, to the coarse space Vh0 , we can immediately conclude that all the coarse solvers are spectrally equivalent thanks to Remark 2.3.

3.1

Schwarz operators

We now define the Schwarz operators and show that they can be viewed as preconditioners for the original (preconditioned) system of equations (6). For the exact local solvers, let PEi : Vh −→ RiT Vhi be defined as −1/2

Dh (PEi u, RiT vi ) := Dh (u, RTi vi ) = Ah (Bh

−1/2 T Ri vi )

u, Bh

∀vi ∈ Vhi .

(33)

For the inexact local solvers we set PIi := RTi PeiI : Vh −→ RTi Vhi ⊂ Vh , where PeiI : Vh −→ Vhi is defined as −1/2 −1/2 (34) DiI (PeiI u, vi ) := Dh (u, RTi vi ) = Ah (Bh u, Bh RTi vi ) ∀vi ∈ Vhi . We observe that the operators PEi and PIi are well-defined since the local bilinear forms DiE (·, ·) and DiI (·, ·) are coercive. We also define the operators P0 , Q0 , T0 : Vh −→ RT0 Vh0 as follows: −1/2

∀v0 ∈ Vh0 ,

Sh (Q0 u, RT0 v0 ) := Dh (u, RT0 v0 ) =

−1/2 T R 0 v0 ) −1/2 −1/2 T Ah (Bh u, Bh R0 v0 )

Sh∗ (T0 u, RT0 v0 ) := Dh (u, RT0 v0 ) =

−1/2 −1/2 Ah (Bh u, Bh RT0 v0 )

∀v0 ∈ Vh0 .

Dh (P0 u, RT0 v0 ) := Dh (u, RT0 v0 ) = Ah (Bh

u, Bh

∀v0 ∈ Vh0 ,

(35)

Since the coarse bilinear forms Dh (·, ·), Sh (·, ·) and Sh∗ (·, ·) are coercive, the operators P0 , Q0 and T0 are well defined. We are now ready to define the following additive Schwarz operators: PE :=

N ∑

PEi + P0 ,

QE :=

i=1

PI :=

N ∑ i=1

N ∑

PEi + Q0 .

TE :=

i=1

PIi + P0 ,

QI :=

N ∑

N ∑

PEi + T0 .

(36)

PIi + T0 .

(37)

i=1

PIi + Q0 .

i=1

TI :=

N ∑ i=1

In the case of exact local solvers, the matrix representation of the additive Schwarz operators PE , QE

10

and TE is given by

( PE =

N ∑

) −1 RTi (DE Ri + i )

RT0 D−1 0 R0

i=1

( Q = E

( T = E

N ∑

Dh := ME1 Dh , )

−1 RTi (DE Ri i )

+

RT0 S−1 0 R0

+

RT0 (S∗0 )−1 R0

Dh := ME2 Dh ,

i=1 N ∑

) −1 RTi (DE Ri i )

Dh := ME3 Dh ,

i=1

S∗0

where D0 , S0 and are the matrix representations of the bilinear forms D0 (·, ·), S0 (·, ·) and S0∗ (·, ·), respectively. We observe that the preconditioners differ by the choice of the coarse solver (cf. Table 1). ME1 employs as coarse solver the restriction of the P-WOPSIP bilinear form to the finite element coarse space whereas for ME2 and ME3 the coarse solver is defined as the restriction to the coarse space of the bilinear form Sh (·, ·) and Sh∗ (·, ·), respectively. The same kind of representation holds for inexact local solvers. Details are given in Table 1.

Preconditioner ME1 := ME2 := ME3 := MI1 := MI2 := MI3 :=

∑N

Coarse Component

Local Components (i = 1, . . . , N )

i=1

−1 Ri + RT0 D−1 RTi (DE i ) 0 R0

D0 := R0 Dh RT0

T DE i := Ri Dh Ri

i=1

−1 Ri + RT0 S−1 RTi (DE i ) 0 R0

S0 := R0 Sh RT0

T DE i := Ri Dh Ri

i=1

−1 Ri + RT0 (S∗0 )−1 R0 RTi (DE i )

S∗0 := R0 S∗h RT0

T DE i := Ri Dh Ri

i=1

RTi (AIi )−1 Ri + RT0 D−1 0 R0

D0 := R0 Dh RT0

see (21)

i=1

RTi (AIi )−1 Ri + RT0 S−1 0 R0

S0 := R0 Sh RT0

see (21)

i=1

RTi (AIi )−1 Ri + RT0 (S∗0 )−1 R0

S∗0 := R0 S∗h RT0

see (21)

∑N ∑N ∑N ∑N ∑N

Table 1: Coarse and local components for the preconditioners ME1 − ME2 − ME3 and MI1 − MI2 − MI3 .

3.2

Computational issues

Let v be a vector representing a finite element function v in the edgewise ordering, and let P be the permutation matrix so that Pv becomes the vector representing v in the elementwise ordering. We define J to be the matrix representing the jumps term ∑ α ∫ T w Jv = Π0 ([[w]]) · Π0e ([[v]]) ds, h3e e e e∈Eh

and G to be the matrix representing the volume term ∑ ∫ ∇w · ∇v dx. wT Gv = T ∈Th

T

We remark that in the elementwise ordering the matrix G is block-diagonal with 3 × 3 blocks, whereas in the edgewise ordering the matrix J is block diagonal and therefore the preconditioner Bh = I + J is −1/2 is block diagonal as well and the 2 × 2 block diagonal, with I the identity matrix. Therefore, Bh 11

−1/2

Algorithm 1 Compute z = Bh

−1/2

Ah Bh

v

1/2 Bh z

Solve =v Compute x := Jz Compute y := PT GPz 1/2 Solve Bh z = x + y

blocks can be computed directly with (7). Algorithm 1 computes the action of the stiffness matrix of −1/2 the WOPSIP method and the action of the preconditioner Bh on a vector (cf. [8]). Next, we also describe the action of the additive Schwarz preconditioner ME1 on a vector v edgewise ordered (cf. Algorithm 2). The routines for the other preconditioners can be written exactly in the same way with only notational changes involved. Note that, for the application of the preconditioner, it is more convenient to employ the elementwise ordering of the dofs, and to number first the dofs corresponding to elements in the first subdomain, then the dofs corresponding to elements in the second subdomain and so on. With such an ordering, the local solvers turn out to be a block Jacobi preconditioner where each block corresponds to the dofs in a subdomain. Algorithm 2 Compute z = ME1 v T Solve z = RT0 D−1 0 R0 Pv for i = 1, . . . , N do T z ⇐ z + RTi D−1 i Ri z end for z ⇐ PT z.

4

Convergence analysis

In this section we present the convergence analysis of the proposed Schwarz methods for the P-WOPSIP scheme. We start by stating the main result of this section: Theorem 4.1. Let P be any of the Schwarz operators defined in (36) and (37). Then the condition number of P satisfies ( H) κ(P) ≤ C 1 + , h where the positive constant C is independent of H, h, α and the number of subdomains. Remark 4.2. Theorem 4.1 guarantees that the condition number of the preconditioned system is independent of the parameter α. This result is in contrast to those for the Schwarz methods for standard and classical DG discretizations of second order problems [4], for which is shown that the condition number of the preconditioned system depends linearly on the penalty parameter α (see [2, 3] for details). The rest of the section is devoted to the proof of the above theorem. We follow the classical abstract convergence theory of Schwarz methods [15, 14] (cf. also [21, Chapter 2] and [12, Chapter 7]), and therefore, we only have to verify the following three assumptions. Assumption A 1 (Stable decomposition). There exists C0 > 0 such that every u ∈ Vh admits a ∑N decomposition u = i=0 RTi ui , with u0 ∈ Vh0 , and ui ∈ Vhi , i = 1, . . . , N , that satisfies N ∑

Dir (ui , ui ) + γ0 (u0 , u0 ) ≤ C02 Dh (u, u),

i=1

12

r = E or I ,

where γ0 (·, ·) is one of the coarse bilinear forms defined in (30)–(32). Assumption A2 (Strengthened Cauchy–Schwarz inequalities). There exist 0 ≤ εij ≤ 1, 1 ≤ i, j ≤ N , such that Dh (RTi ui , RTj uj ) ≤ εij Dh (RTi ui , RTi ui )1/2 Dh (RTj uj , RTj uj )1/2 for all vi ∈ Vhi , uj ∈ Vhj . Define ρ(E) to be the spectral radius of E := {εij }i,j=1,...,N . Assumption A3 (Local stability). There exists ωr > 0 such that Dh (RTi ui , RTi ui ) ≤ ωr Dir (ui , ui )

∀ ui ∈ Vhi , r = {E, I}.

(38)

Under Assumptions A1–A3, we have the following condition number estimate: κ(P) ≤ C02 ωr (ρ(E) + 1). Therefore we can prove Theorem 4.1 by showing that (i) ρ(E) ≤ Na + 1, where Na is the maximum number of adjacent subdomains that a given subdomain might have, (ii) ωr . 1, and (iii) C02 . 1 + (H/h). We start by verifying the bound for ρ(E). Following [19, 2], it is straightforward to see that εii = 1 −1/2 T −1/2 T Rj uj ) ̸= 0 only if for i = 1, . . . , N . For i ̸= j, we note that Dh (RTi ui , RTj uj ) = Ah (Bi Ri ui , B j ∂Ωi ∩ ∂Ωj ̸= ∑ ∅, therefore εij = 1 in those cases, and εij = 0 otherwise. Then, ρ(E) can be bounded by ρ(E) ≤ maxi j |εij | ≤ 1 + Na . In the next sections we verify the bounds for ωr and C0 .

4.1

Local stability

We now prove that the local solvers satisfy a local stability property with ωr . 1. Observe that for the exact local solvers defined in (19), it follows from their definition that (38) holds with ωE ≡ 1. Before showing that Assumption A3 holds also for the inexact local solvers, we define the norm ∥ · ∥DG,Ωi according to (13) but at the subdomain level, i.e., ∥ui ∥2DG,Ωi :=

∑ T ∈Th T ⊂Ωi

∥∇ui ∥20,T +

∑ 1 ∑ 1 ∥Π0e ([[ui ]])∥20,e + ∥Π0 (ui n)∥20,e he he e

e∈Eh e⊂Ωi

∀ui ∈ Vhi ,

e∈Eh e⊂∂Ωi

and observe that the coercivity (16) holds also at the subdomain level for h ≤ h0 with h0 given in (14), by the definition of DiI (·, ·). Lemma 4.3. For i = 1, . . . N , let DiI : Vhi × Vhi −→ R be the bilinear form defined by (21). Then, there exists a positive number ωI . 1 such that the following local stability property holds: Dh (RTi ui , RTi ui ) ≤ ωI DiI (ui , ui )

∀ ui ∈ Vhi

∀ i = 1, . . . , N .

Proof. Observe that ∥RiT ui ∥2DG = ∥ui ∥2DG,Ωi . It then follows from (15) and (16) (for DiI (·, ·)) that Dh (RiT ui , RTi ui ) ≤ Cc ∥RiT u∥2DG ≤ Cc ∥ui ∥2DG,Ωi ≤ CDiI (ui , ui ) .

13

4.2

Stable decomposition

In this section we finally show that the decomposition underlying the definition of the additive Schwarz operator is indeed stable with respect to the energy norm defined by Dh (·, ·). We first state an auxiliary result needed in the proof of Proposition 4.6. This result provides an estimate for the interface bilinear forms IhE (·, ·) and IhI (·, ·). Lemma 4.4. For any u ∈ Vh , it holds that ∑ ∑ |Ihr (u, u)| . ∥u∥2DG + h−1 ∥u∥20,E , r = E or I . (39) D∈TH E⊂∂D

Proof. We start by proving the bound for IhE (·, ·). From the definition (26) of IhE given in Lemma 3.1 and the standard triangle inequality, we have |IhE (u, u)| ≤ 2 |F1 | + 2 |F2 |, where ∑ ∫ −1/2 −1/2 ∇(Bh RTi ui ) · ∇(Bh RTj uj ) dx, F1 := T ∈ΩΓ

F2 :=

∑ e∈Γ

βe2

T

α h3e

∫ Π0e (ui )Π0e (uj ) ds. e

We next estimate the two terms separately, starting with F2 . By recalling the definition (8) of βe on e ∈ Eh◦ and using (9), the Cauchy-Schwarz inequality, the arithmetic-geometric inequality and the stability of the projection Π0e (·), we find (

)1/2 ( )1/2 ∑ 1 0 2 ∥Π (uj )∥0,e he e e∈Γ e∈Γ ) ∑ 1 ∑ 1 ∑( 1 1 ≤ ∥Π0e (ui )∥20,e + ∥Π0e (uj )∥20,e ≤ ∥ui ∥20,e + ∥uj ∥20,e . he he he he

∑ 1 2 |F2 | ≤ 2 ∥Π0 (ui )∥20,e he e

e∈Γ

e∈Γ

e∈Γ

Observe that each subdomain Ωi is the union of some elements D ∈ TH and the mesh is quasi-uniform. Denoting by E the edges of D, we have ∑ ∑ ∑ 2 −1 2 −1 ∥u∥20,E . 2 |F2 | ≤ (h−1 (40) e ∥ui ∥0,e + he ∥uj ∥0,e ) . h e∈Γ

D∈TH E⊂∂D

Next, we estimate the term |F1 |. Using the Cauchy-Schwarz inequality we have N ∑ ∑ ∫ −1/2 T −1/2 T |F1 | = ∇(Bh Ri ui ) · ∇(Bh Rj uj ) dx T i,j=1 T ∈ΩΓij i̸=j

.

N ( ∑ ∑ i,j=1 i̸=j

−1/2 T Ri ui ∥0,T

∥∇Bh

−1/2 T Rj uj ∥0,T

∥∇Bh

) .

(41)

T ∈ΩΓij

Observe that for any fixed j ̸= i and T ∈ ΩΓij with T ⊆ Ωi , the divergence theorem, the CauchySchwarz and the trace inequalities together with the stability of the projection Π0e (·) give ∫ ∫ −1/2 T −1/2 T −1/2 T −1/2 −1/2 2 ∥∇Bh Rj uj ∥0,T = − ∆(Bh Rj uj ) Bh Rj uj dx + ∇(Bh RTj uj ) · n Bh RTj uj ds ∂T ∫ T −1/2 T −1/2 T 0 = ∇(Bh Rj uj ) · n Πe (Bh Rj uj ) ds ∂T

14

−1/2 T Rj uj )∥0,∂T

≤ ∥∇(Bh

−1/2 T Rj uj )∥0,T

−1/2 T Rj uj )∥0,∂T

∥Π0e (Bh

h−1/2 ∥RTj uj ∥0,e e

. ∥∇(Bh

(e = ∂T ∩ Γij ) ,

−1/2 T Rj uj )

where in the last step we have used the fact that Π0e (Bh due to (53) (see also Figure 4), and hence −1/2 T Rj uj ∥0,T

∥∇Bh

. h−1/2 ∥uj ∥0,e e

e = ∂T ∩ Γij

(42)

̸= 0 only on the edge e = ∂T ∩ Γij T ⊂ Ωi ∩ ΩΓij

i ̸= j.

After inserting the estimate (42) into (41) and using the continuity (15) of Dh (·, ·) given in Lemma 2.1, we finally obtain |F1 | .

N ( ∑ ∑ i,j=1 i̸=j

. .

−1

+h

−1/2 T Ri uj ∥0,T h−1/2 ∥ui ∥0,e e

∥∇Bh

)

T ∈ΩΓij T ⊂Ωj

∥∇u∥20,T

T ∈ΩΓ

∑

+

T ∈ΩΓij T ⊂Ωi

∑ ( ∥u∥2DG

−1/2 T Ri ui ∥0,T h−1/2 ∥uj ∥0,e e

∥∇Bh

∑

+

∑

2 h−1 e ∥Π0 ([[u]])∥0,e

)

−1

+h

∑

(e = ∂T ∩ Γij ) ∥u∥20,E

D∈TH E⊂∂D

e⊂∂T

∑

∑

∥u∥20,E .

D∈TH E⊂∂D

The above estimate together with (40) concludes the proof for IhE (·, ·). To bound IhI (·, ·) we observe that, thanks to Lemma 3.1 IhI (u, u) = IhE (u, u) + GIh (u, u)

∀ u ∈ Vh ,

and therefore it is enough to bound GIh (·, ·) which we recall is defined as

GIh (u, u) =

N ∑ ∫ ∑ i=1 T ∈ΩΓ

−

T

N ∑ ∫ ∑ i=1 T ∈ΩΓ T ⊂Ωi

−

−1/2 T Ri ui )

∇(Bh

−1/2

∇(Bi

−1/2 T Ri ui ) dx

· ∇(Bh

−1/2

ui ) · ∇(Bi

(F 3)

ui ) dx

(F 4)

T

∫ N ∑ ∑ ) α 2 ( 0 η [Πe (ui )]2 + [Π0e (uj )]2 ds . e 3 he e i=1

(F 5)

e∈Γ

We start with the last term F5 . Recalling the definition (29) of ηe and using the fact that α2 /(α + h2e )(2α + h2e ) ≤ 1, we have α 2 1 1 α θe α2 = ≤ . η = e 3 3 2 2 he he (1 + θe )(2 + θe ) he (α + he )(2α + he ) he Then, taking into account the stability of the projection Π0e (·) and arguing as we did for F2 , we obtain |F5 | ≤

N ∑ N ∑ ∑ ∑ ∑ ) α 2( 1 2 2 ∥ui ∥20,e . h−1 η ∥u ∥ + ∥u ∥ . 2 i j e 0,e 0,e 3 he he i=1 i=1 e∈Γ

e∈Γ

15

∑

D∈TH E⊂∂D

∥u∥20,E .

We now estimate the other terms. The estimate for F3 is similar to the estimate for |F1 |: |F3 | ≤

N ∑ ∑

−1/2 T Ri ui )∥20,T

∥∇(Bh

i=1 T ∈ΩΓ T ⊂Ωi

.

∑ (

∥∇u∥20,T +

T ∈ΩΓ

.

∥u∥2DG

−1

+h

∑

∑

+

N ∑ ∑

−1/2 T Ri ui )∥20,T

∥∇(Bh

i=1 T ∈ΩΓ T ̸⊂Ωi

) ∑ 2 −1 h−1 e ∥Π0 ([[u]])∥0,e + h

∥u∥20,E

D∈TH E⊂∂D

e⊂∂T

∑

∑

∥u∥20,E .

D∈TH E⊂∂D

Finally, the term F4 is readily estimated by the continuity of the bilinear form DiI (·, ·): |F4 | ≤

N ∑ ∑

−1/2 ui )∥20,T ∥∇(Bi

i=1 T ∈ΩΓ T ⊂Ωi

.

N ∑

∥ui ∥2DG,Ωi . ∥u∥2DG + h−1

i=1

∑

∥u∥20,E .

D∈TH E⊂∂D

The last preliminary result concerns the coarse solver. Lemma 4.5. For any u ∈ Vh , let u0 ∈ Vh0 = VH be defined as ∫ 1 u0 |D := u dx ∀D ∈ TH . |D| D Then it holds that

∑

( ) γ0 (u0 , u0 ) . 1 + Hh−1 Dh (u, u),

(43)

(44)

where γ0 (·, ·) is one of the coarse bilinear forms defined in (30)–(32). Proof. It is sufficient to show the bound in the case γ0 (·, ·) = S0∗ (·, ·); the other two cases follow from the observation made in Remark 3.3. Let u ∈ Vh and let u0 be defined as in (43). Note that u0 is piecewise constant (by definition) on TH . Then it follows from the definition of S0∗ (·, ·), adding and subtracting u and the stability of the projection Π0 (·) that ∑ α ∫ ∗ Π0e ([[R0T u0 ]]) 2 ds S0 (u0 , u0 ) = he e e∈Eh ∑ 1 ∫ ∑ 1 ∫ 2 T 0 Π0 ([[u]]) 2 ds . Π ([[R0 u0 − u]]) ds + he e e he e e e∈Eh e∈Eh ∑ 1 ∫ 2 [[R0T u0 − u]] ds + Dh (u, u) , . he e e∈Eh

where in the last step we have also used the coercivity of Dh (·, ·) (cf. (16)). We now observe that last term can be estimated exactly following [19] and [2, Lemma 4.3]: ∑ 1 ∫ [[R0T u0 − u]] 2 ds . Hh−1 ∥u∥2DG . Hh−1 Dh (u, u) . he e e∈Eh

We close the section with the proof of Assumption A1.

16

∑N Proposition 4.6 (Stable decomposition). For any u ∈ Vh , let u = i=0 RTi ui , ui ∈ Vhi , i = 0, . . . , N , where u0 ∈ Vh0 is defined by ∫ 1 u0 |D := u dx ∀D ∈ TH , |D| D and u1 , . . . , uN are (uniquely) determined by u − RT0 u0 = RT1 u1 + · · · + RTN uN . Then, there exists C02 . 1 + (H/h) such that N ∑

Dir (ui , ui ) + γ0 (u0 , u0 ) ≤ C02 Dh (u, u), r = {E, I} ,

i=1

where γ0 (·, ·) is one of the coarse bilinear forms defined in (30)–(32). Proof. The proof follows those given in [19, 2]. We set γ0 (·, ·) = D0 (·, ·). Given u ∈ Vh , we decompose ∑N u − RT0 u0 uniquely as i=1 RTi ui . Taking into account Lemma 3.1 we can write D0 (u0 , u0 ) +

N ∑

Dir (ui , ui ) = D0 (u0 , u0 ) + Dh (u − R0T u0 , u − R0T u0 ) − Ihr (u − R0T u0 , u − R0T u0 ) , (45)

i=1

then we just need to estimate each term on the right hand side. The first term is readily estimated by using Lemma 4.5: ( ) D0 (u0 , u0 ) . 1 + Hh−1 Dh (u, u).

(46)

For the second term on the right hand side of (45), triangle inequality, the continuity of Dh (·, ·) (cf. (15)) together with (46) and the definition of the coarse solver gives, Dh (u − R0T u0 , u − R0T u0 ) . Dh (u, u) + Dh (R0T u0 , R0T u0 ) ( ) = Dh (u, u) + D0 (u0 , u0 ) . 1 + Hh−1 Dh (u, u).

(47)

For the last term, it follows from (16), Lemma 4.4 and (47) that ∑ ∑ r Ih (u − R0T u0 , u − R0T u0 ) . ∥u − R0T u0 ∥2DG + h−1 ∥u − R0T u0 ∥20,E D∈TH E⊂∂D

∑ ( ) . 1 + Hh−1 Dh (u, u) + h−1

∑

∥u − R0T u0 ∥20,E .

D∈TH E⊂∂D

Noting now that the trace inequality together with a Poincar´e-Friedrichs inequality [7, 12] gives [ ] ∑ ∑ ∑ −1 −1 h−1 ∥u − R0T u0 ∥20,E . HD h ∥u − R0T u0 ∥20,D + HD h−1 ∥∇h (u − R0T u0 )∥20,D D∈TH E⊂∂D

D∈TH

. Hh−1 ∥u∥2DG + Hh−1

∑

∥∇u∥20,T . Hh−1 Dh (u, u) ,

T ∈Th

we finally obtain the estimate r ( ) Ih (u − R0T u0 , u − R0T u0 ) . 1 + Hh−1 Dh (u, u) for the last term in (45). Substituting this estimate together with (47) and (46) into (45), the proof is completed. For the other coarse solvers, the proof follows exactly the same steps, replacing the bound in (46) by the corresponding one.

17

5

Numerical results

In this section we present a series of numerical experiments to highlight the practical performance of our non-overlapping Schwarz preconditioners. We restrict ourselves to two-dimensional model problems: we let Ω = (0, 1) × (0, 1) and choose f such that the analytical solution of the model problem (1) is given by u(x, y) = exp(xy)(x − x2 )(y − y 2 ). Throughout Sections 5.1 and 5.2 we take the stability constant α appearing in the formulation of the WOPSIP method (3) to be 1; numerical results with different choices of the penalty parameter are discussed in Section 5.3. We employ a uniform subdomain partition consisting of N = 4, 16 squares, and consider initial coarse and fine refinements as depicted in Figure 2 (for N = 4 (top) and N = 16 (bottom)). We denote by H0 and h0 the corresponding initial coarse and fine mesh sizes, respectively, and consider j = 1, 2, 3, successive uniform refinements of the initial grids.

Figure 2: Initial coarse and fine refinements, respectively, on N = 4 subdomain partitions (top) and N = 16 subdomain partitions. We solved the preconditioned linear systems of equations by the conjugate gradient (CG) iterative solver with a (relative) tolerance set equal to 10−9 allowing a maximum of 100 iterations. For the solution of the linear system (6) we employed the CG iterative solver with the same relative tolerance but allowing a maximum of 1100 iterates.

5.1

Exact local solvers

In this section we test the performance of the Schwarz preconditioners ME1 =

N ∑

−1 RTi (DE Ri + RT0 D−1 i ) 0 R0 ,

i=1

ME2 =

N ∑

−1 RTi (DE Ri + RT0 S−1 i ) 0 R0 ,

i=1

ME3 =

N ∑

−1 RTi (DE Ri + RT0 (S∗0 )−1 R0 , i )

i=1

18

applied to the original symmetric (preconditioned) systems of equations (6). In the first set of experiments we have considered a coarse space constructed from piecewise linear discontinuous elements. The condition number estimates together with the corresponding iteration counts needed to reach convergence (in parenthesis) for all the considered preconditioners are reported in Table 2 on a partition with 16 subdomains. For the sake of comparison, we also report (last but one row of Table 2) the con−1/2 −1/2 together with the iteration counts needed for the dition number estimate of the matrix Bh Ah Bh solution of the linear system of equations (6). The last row of Table 2 shows the condition number and Preconditioner ME1 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

5.1815 (23) -

5.9422 (25) 5.2452 (23) -

9.0724 (33) 5.9113 (26) 5.3608 (23) -

h0 /8 17.1532 8.9409 5.9822 5.4379

(45) (33) (26) (23)

Preconditioner ME2 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

6.6250 (27) -

7.6977 (30) 7.5231 (29) -

11.6106 (37) 8.0398 (31) 7.9609 (30) -

h0 /8 21.7983 11.9791 8.3793 8.2273

(49) (39) (32) (31)

Preconditioner ME3 H↓ h→

h0

h0 /2

h0 /4

6.6459 (27) -

7.7096 (30) 7.5603 (29) -

11.6134 (37) 8.0511 (31) 7.9947 (30) -

Ah Bh

1.2229e+2 ( 52)

4.7212e+2 (100)

1.8726e+3 (202)

7.4752e+3 ( 410)

Ah

2.5173e+4 (115)

3.9949e+5 (229)

6.3800e+6 (486)

1.0309e+8 (1040)

H0 H0 /2 H0 /4 H0 /8 −1/2

Bh

−1/2

h0 /8 21.7993 11.9826 8.3903 8.2664

(49) (39) (32) (31)

Table 2: Preconditioners ME1 , ME2 and ME3 (N = 16, α = 1): condition number estimates and iteration counts. Piecewise linear discontinuous coarse space. the corresponding iteration counts of the original unpreconditioned system of equations Ah x = f. The numerical results confirm the theoretical estimates provided in Theorem 4.1: the condition number of the preconditioned system √ behaves asymptotically as H/h, and, consequently, the iteration counts behaves asymptotically as H/h. By a comparison with the computed condition number of the matrix Ah it is clear that the application of all the preconditioners drastically reduce the condition number of the system, and consequently, the iteration counts needed for convergence. Next, we investigate the scalability of the preconditioners, i.e., the independence of the performance on the number of subdomains. To this end we repeated the same set of experiments decreasing the number of subdomains from N = 16 to N = 4: the results are reported in Table 3. As predicted from our theoretical estimates, the condition number of the preconditioned system is independent of the number of subdomains.

19

Preconditioner ME1 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

4.8722 (21) -

5.7042 (24) 5.0710 (22) -

8.9088 (30) 5.8873 (25) 5.3561 (23) -

h0 /8 16.5838 9.0200 5.9792 5.4471

(40) (32) (26) (23)

Preconditioner ME2 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

6.8843 (27) -

7.7219 (29) 7.4641 (29) -

11.8249 (37) 8.2226 (31) 8.0459 (30) -

h0 /8 21.5324 12.1736 8.4641 8.2667

(50) (38) (32) (31)

Preconditioner ME3 H↓ h→

h0

h0 /2

h0 /4

6.9179 (27) -

7.7294 (29) 7.5241 (29) -

11.8275 (37) 8.2336 (31) 8.0865 (30) -

Ah Bh

1.2229e+2 ( 52)

4.7212e+2 (100)

1.8726e+3 (202)

7.4752e+3 ( 410)

Ah

2.5173e+4 (115)

3.9949e+5 (229)

6.3800e+6 (486)

1.0309e+8 (1040)

H0 H0 /2 H0 /4 H0 /8 −1/2

Bh

−1/2

h0 /8 21.5322 12.1763 8.4768 8.3095

(50) (38) (32) (31)

Table 3: Preconditioners ME1 , ME2 and ME3 (N = 4, α = 1): condition number estimates and iteration counts. Piecewise linear discontinuous coarse space.

20

Next, we investigate the effect of the coarse space on the performance of our preconditioners. To this end, we ran the same set of experiments as before (on a partition with 16 subdomains) employing a piecewise constant coarse space. Table 4 reports the condition number estimates and the corresponding iteration counts. Note that whenever we employ a piecewise constant coarse space the bilinear forms Sh (·, ·) and Sh∗ (·, ·) turn out to be identical, and therefore the preconditioners ME2 and ME3 coincide. For this reason, in Table 4 we only report the results obtained with the preconditioners ME1 and ME2 . We observe that the performance of the preconditioners are consistently poorer. On the other hand with this choice of coarse space only one degree of freedom per coarse element is required. Preconditioner ME1 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

9.1821 (27) -

19.1929 (38) 10.6855 (31) -

39.5807 (50) 22.1790 (44) 11.8751 (36) -

h0 /8 80.6716 45.2087 24.4631 12.7832

(70) (64) (50) (38)

Preconditioner ME2 H↓ h→

h0

h0 /2

h0 /4

9.4006 (29) -

20.2353 (40) 11.6003 (35) -

42.7406 (54) 23.9082 (48) 13.2304 (38) -

Ah Bh

1.2229e+2 ( 52)

4.7212e+2 (100)

1.8726e+3 (202)

7.4752e+3 ( 410)

Ah

2.5173e+4 (115)

3.9949e+5 (229)

6.3800e+6 (486)

1.0309e+8 (1040)

H0 H0 /2 H0 /4 H0 /8 −1/2

Bh

−1/2

h0 /8 88.3442 49.1192 26.8473 14.4078

(76) (68) (55) (41)

Table 4: Preconditioners ME1 and ME2 (N = 16, α = 1): condition number estimates and iteration counts. Piecewise constant discontinuous coarse space.

5.2

Inexact local solvers

In this section we test the performance of the Schwarz preconditioners (with inexact local solvers) MI1 =

N ∑

RTi (AIi )−1 Ri + RT0 D−1 0 R0 ,

i=1

MI2 =

N ∑

RTi (AIi )−1 Ri + RT0 S−1 0 R0 ,

i=1

MI3 =

N ∑

RTi (AIi )−1 Ri + RT0 (S∗0 )−1 R0 ,

i=1

applied to the original symmetric (preconditioned) systems of equations (6). We ran the same set of experiments as before. More precisely, in Table 5 and Table 6 we compare the condition number estimates and the iteration counts obtained on a subdomain partition made of N = 16 and N = 4 subdomains, respectively, employing a piecewise linear discontinuous coarse space. As expected, the preconditioners with inexact local solvers are also scalable, and the condition number 21

Preconditioner MI1 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

4.2990 (21) -

7.3572 (28) 4.7076 (22) -

14.3928 (40) 7.7532 (30) 4.9095 (23) -

h0 /8 29.5718 14.4282 8.1509 5.0121

(57) (41) (30) (23)

Preconditioner MI2 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

5.5549 (25) -

9.4431 (32) 6.2441 (26) -

18.7975 (45) 10.1556 (34) 6.5799 (28) -

h0 /8 38.8577 18.7271 10.6361 6.6286

(64) (47) (35) (29)

Preconditioner MI3 H↓ h→

h0

h0 /2

h0 /4

5.5855 (25) -

9.4545 (32) 6.2788 (27) -

18.8013 (45) 10.1688 (35) 6.6139 (28) -

Ah Bh

1.2229e+2 ( 52)

4.7212e+2 (100)

1.8726e+3 (202)

7.4752e+3 ( 410)

Ah

2.5173e+4 (115)

3.9949e+5 (229)

6.3800e+6 (486)

1.0309e+8 (1040)

H0 H0 /2 H0 /4 H0 /8 −1/2

Bh

−1/2

h0 /8 38.8591 18.7303 10.6507 6.6689

(64) (47) (35) (29)

Table 5: Preconditioners MI1 , MI2 and MI3 (N = 16, α = 1): condition number estimates and iteration counts. Piecewise linear discontinuous coarse space.

22

estimates of the preconditioned system are in agreement with Theorem 4.1: the computed condition number seems to behave as O(H/h). Preconditioner MI1 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

4.4817 (20) -

7.3013 (26) 4.7381 (21) -

13.8371 (36) 7.8524 (28) 4.7748 (22) -

h0 /8 27.7972 14.4631 8.1640 4.9653

(51) (38) (29) (23)

Preconditioner MI2 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

5.9734 (24) -

9.7880 (31) 6.4031 (26) -

18.2750 (42) 10.4022 (34) 6.6069 (28) -

h0 /8 36.4445 18.9555 10.7379 6.5152

(58) (45) (35) (28)

Preconditioner MI3 H↓ h→

h0

h0 /2

h0 /4

6.0084 (24) -

9.8017 (31) 6.4362 (27) -

18.2784 (42) 10.4172 (34) 6.6491 (28) -

Ah Bh

1.2229e+2 ( 52)

4.7212e+2 (100)

1.8726e+3 (202)

7.4752e+3 ( 410)

Ah

2.5173e+4 (115)

3.9949e+5 (229)

6.3800e+6 (486)

1.0309e+8 (1040)

H0 H0 /2 H0 /4 H0 /8 −1/2

Bh

−1/2

h0 /8 36.4448 18.9593 10.7544 6.5739

(58) (45) (34) (28)

Table 6: Preconditioners MI1 , MI2 and MI3 (N = 4, α = 1): condition number estimates and iteration counts. Piecewise linear discontinuous coarse space. Finally, we test again the performance of the preconditioners where the coarse spaces are constructed from piecewise constant polynomials. The computed condition number estimates and the corresponding iteration counts obtained by employing the preconditioners MI1 and MI2 are shown in Table 7. By comparing the results with the analogous ones presented in the previous Section 5.1 it can been inferred that employing exact local solvers improves the performance of the preconditioner slightly.

23

Preconditioner MI1 H↓ h→ H0 H0 /2 H0 /4 H0 /8

h0

h0 /2

h0 /4

10.6099 (29) -

24.6977(43) 12.2398 (34) -

54.5351 (63) 26.7205 (49) 13.2809 (37) -

h0 /8 114.8625 56.2146 28.0078 14.0256

(91) (71) (51) (38)

Preconditioner MI2 H↓ h→

h0

h0 /2

h0 /4

12.0117 (31) -

28.4619 (46) 14.0306 (36) -

62.7619 (67) 30.7275 (52) 15.1082 (40) -

Ah Bh

1.2229e+2 ( 52)

4.7212e+2 (100)

1.8726e+3 (202)

7.4752e+3 ( 410)

Ah

2.5173e+4 (115)

3.9949e+5 (229)

6.3800e+6 (486)

1.0309e+8 (1040)

H0 H0 /2 H0 /4 H0 /8 −1/2

Bh

−1/2

h0 /8 132.3684 64.7140 31.8235 15.8052

(95) (74) (55) (42)

Table 7: Preconditioners MI1 and MI2 (N = 16, α = 1): condition number estimates and iteration counts. Piecewise constant discontinuous coarse space.

5.3

Variable penalty parameter

The aim of this section is to validate the independence on the penalty parameter of the estimates for the condition number of the preconditioned system proved in Theorem 4.1. For the sake of brevity throughout this section we employ a piecewise constant coarse solver. In Figure 3 we report the condition number estimates of the preconditioned system for different values of α. Although our theory requires α ≥ 1, for the sake of completeness we report here the results obtained with α = 10−2 , . . . , 104 , and different mesh configurations. Results obtained with exact local solvers, i.e., the preconditioners ME1 , ME2 and ME3 are shown in Figure 3(a), whereas Figure 3(b) shows the analogous results obtained with inexact local solvers, i.e., the preconditioners MI1 , MI2 and MI3 . As expected, our preconditioner is fairly insensitive on the choice of the penalty parameter..

Acknowledgments The work of P.F. Antonietti and B. Ayuso de Dios was partially supported by Azioni Integrate Italia– Spagna through the projects IT097ABB10 and HI2008-0173. B. Ayuso de Dios was also partially supported by MEC through the grant MTM2011-27739-C04-04. The work of S.C. Brenner and L.-Y. Sung was supported in part by the National Science Foundation under Grant DMS-10-16332.

A

Appendix

We prove Lemma 2.1 in this appendix. For that purpose, we first recall a result that provides a natural splitting of the DG linear functions. Proposition A.1. [5, Proposition 3.1] For any u ∈ Vh there exist a unique v ∈ VhCR and a unique z ∈ Zh such that u = v + z. That is: Vh = VhCR ⊕ Zh , where VhCR is the classical Crouziex-Raviart

24

(a) Exact local solvers 6

6

10

10

k(ME 1 Dh) k(ME 2 Dh) k(ME 3 Dh)

k(ME 1 Dh) k(ME 2 Dh) k(ME 3 Dh) k(Dh )

k(Dh ) 4

4

10

10

2

2

10

10 −2

10

0

10

2

4

10

10

−2

10

0

10

2

10

4

10

(b) Inexact local solvers 6

6

10

10

k(MI1 D h ) k(MI2 D h ) k(MI3 D h )

k(MI1 D h ) k(MI2 D h ) k(MI3 D h ) k(Dh )

k(Dh ) 4

4

10

10

2

2

10

10 −2

10

0

10

2

10

4

10

−2

10

0

10

2

10

Figure 3: Condition number estimates as a function of the penalty parameter α (N = 16): h = h0 /8, H = H0 (left) and h = h0 /8, H = H0 /2 (right). Piecewise constant coarse solver.

25

4

10

space defined by { VhCR := v ∈ L2 (Ω) : v|T ∈ P1 (T ) and the space Zh is defined by: { Zh := v ∈ L2 (Ω) : v|T ∈ P1 (T )

∀T ∈ Th

∀T ∈ Th

and

and

Π0e ([[v]]) = 0

} ∀ e ∈ Eh◦ .

Π0e ({{v}}) = 0

} ∀ e ∈ Eh◦ .

A natural set of basis functions associated to midpoints of edges can be given for both spaces VhCR and Zh , i.e., Zh Zh ◦ }e∈Eh◦ ⊕ span{ψe,T }e∈Eh∂ . Zh = span{ψe,T VhCR = span{φCR (48) e }e∈Eh Therefore, an edgewise ordering of the dofs of any u ∈ Vh facilitates the use of the above splitting. For any u ∈ Vh , let v ∈ VhCR and z ∈ Zh such that u = v + z. Now, we fix an interior edge e = ∂T + ∩ ∂T − , e ∈ Eh◦ and, we denote by ve (resp. ze ) the vector containing the degrees of freedom of v|e∩T + and v|e∩T − (resp. z|e∩T + and z|e∩T − ). Using the definition of the spaces VhCR and Zh , it follows that, [ ] [ + ] [ ] ] [ ve |ze | u ve + |ze | ve = , ze = , ue = = = ve + ze , ve −|ze | u− ve − |ze | Therefore, we have that (Beh )−1/2 ue =

[

ve + βe |ze | ve − βe |ze |

] = ve + βe ze

∀e ∈ Eh◦ ,

where βe is defined as in (8). Therefore, with such a decomposition the action of the operator (Beh )−1/2 can be read as the one that, on each interior edge, leaves untouched the Crouziex-Raviart part of the DG function and acts only on its highly oscillatory component z. Analogously, on each boundary edge e ∈ Eh∂ , we have (Beh )−1/2 QT ue = βe∂ ze ∀e ∈ Eh∂ , where, following [5], we have assigned the dofs corresponding to the boundary edges (of the Dirichlet problem) in Zh (or, analogously, the Dirichlet boundary conditions with Crouzeix-Raviart elements are imposed strongly). Summarizing, we have { ve + βe ze if e ∈ Eh◦ , e −1/2 (Bh ) ue = (49) βe ze if e ∈ Eh∂ , where βe is defined in (8). We also recall the following result from [5]. We report the proof for the sake of completeness. Lemma A.2. For any z ∈ Zh it holds that ∑

∥∇z∥20,T ≤ Ct2

T ∈Th

∑ 1 ∥Πe ([[z]])∥20,e , he 0

(50)

e∈Eh

where Ct is the constant in the trace inequality (and therefore depends only on the shape regularity of the mesh). Proof. Integrating by parts, recalling that since z ∈ Zh is piecewise linear then ∆z = 0 on each T ∈ Th , and the definition (4) of the operator Πe0 (·) yield ∑∫ ∑∫ ∑ ∫ ∑ 2 {{∇z}} · [[z]] ds + [[∇z]] · {{z}} ds (∆z)z dx + ∥∇z∥0,T = (∇z, ∇z)0,T = − T ∈Th

T ∈Th

T

e∈Eh

26

e

◦ e∈Eh

e

=

∑∫ e∈Eh

=

∑∫

e∈Eh

≤

∑

∑∫

{{∇z}} · Πe0 ([[z]]) ds + e

[[∇z]] · Πe0 ({{z}}) ds

e

◦ e∈Eh

{{∇z}} · Πe0 ([[z]]) ds e

−1/2 h1/2 ∥Πe0 ([[z]])∥0,e . e ∥ {{∇z}} ∥0,e he

e∈Eh

The estimate (50) follows by employing the standard trace inequality. Finally, we are ready to prove Lemma 2.1. Proof of Lemma 2.1. From the decomposition of the space Vh given in Proposition A.1, any u ∈ Vh can be decomposed uniquely as u = ucr + uz , with ucr ∈ VhCR and uz ∈ Zh . Recalling now that −1/2 Bh u = ucr + βuz , with β|e = βe defined as in (8), we find ∑ α ∑ α ∑ α −1/2 0 0 cr z 2 2 ∥Π ([[B = ∥Π ([[u + β u ]])∥ = βe2 3 ∥Π0 ([[uz ]])∥20,e u]])∥ e 0,e 0,e h h3e h3e he e∈Eh e∈Eh e∈Eh ∑ 1 ≤ ∥Π0 ([[u]])∥20,e , he e∈Eh

where the last bound follows from estimate (9). We now show the continuity (15). For any u, w ∈ Vh , we write u = ucr + uz and w = wcr + wz with ucr , wcr ∈ VhCR and uz , wz ∈ Zh . Then, the Cauchy-Schwarz inequality and the above inequality, gives ( Dh (u, w) ≤

∑

)1/2 ( ∥∇(u

cr

+ βu

z

)∥20,T

T ∈Th

( +

)1/2

∑

∥∇(w

cr

+ βw

T ∈Th

∑ 1 ∥Π0 ([[u]])∥20,e he

z

)∥20,T

)1/2 (

e∈Eh

∑ 1 ∥Π0 ([[w]])∥20,e he

)1/2 . (51)

e∈Eh

We now estimate each term on the right hand side separately (it is enough to do this for u = ucr + βuz ). The triangle inequality, the arithmetic-geometric inequality together with Lemma A.2, and the definition (8) of β yield the estimate ∑ ∑ ( ) ∥∇(ucr + βuz )∥20,T . ∥∇ucr ∥20,T + β 2 ∥∇uz ∥20,T T ∈Th

T ∈Th

=

∑ (

∥∇(ucr + uz − uz )∥20,T + β 2 ∥∇uz ∥20,T

T ∈Th

.

∑

∥∇(ucr + uz )∥20,T +

T ∈Th

.

∑

T ∈Th

)

∑ (1 + β 2 ) e ∥Πe0 ([[uz ]])∥20,e he

e∈Eh

∥∇u∥20,T

∑ 1 + ∥Πe ([[u]])∥20,e . he 0

(52)

e∈Eh

which yields (15). We now prove the coercivity. Cauchy-Schwarz inequality, the arithmetic-geometric inequality and estimate (50) from Lemma A.2 imply ∫ 2 β∇h ucr · ∇h uz dx ≤2∥β∇h ucr ∥0,Ω ∥∇h uz ∥0,Ω Ω

27

1 ∥∇h uz ∥20,Ω , 8Ct2 1 ∑ 1 ≤8Ct2 β 2 ∥∇h ucr ∥20,Ω + ∥Π0 ([[uz ]])∥20,e , 8 he ≤8Ct2 β 2 ∥∇h ucr ∥20,Ω +

e∈Eh

where Ct denotes the constant for the trace inequality. Moreover the scaling of βe given in (9) and the assumption α ≥ 1 yield ( ) ( ) 1 1 α 1 1 1 α 2 α ≥ βe 3 = ≥ ≥ ≥ , he he he kα + h2e he kα + 1 2khe 4he with k = 2 if e ∈ Eh◦ and k = 1 if e ∈ Eh∂ . These observations together finally give ∫ ∑ cr 2 2 z 2 2 α 0 z 2 cr z Dh (u, u) ≥ ∥∇h u ∥0,Ω + β ∥∇h u ∥0,Ω + βe 3 ∥Π ([[u ]])∥0,e − 2β ∇h u · ∇h u dx he Ω e∈Eh ∑ ( ) 1 ≥ 1 − 8Ct2 β 2 ∥∇h ucr ∥20,Ω + β 2 ∥∇h uz ∥20,Ω + ∥Π0 ([[u]])∥20,e , 8he e∈Eh

Hence, by taking h so that 1 − 8Ct2 β 2 ≥ 1/2 > 0, we have Dh (u, u) ≥

∑ ∑ 1 1 ∥∇h ucr ∥20,Ω + β 2 ∥∇uz ∥20,T + ∥Π0 ([[u]])∥20,e , 2 8he T ∈Th

e∈Eh

and therefore, discarding the low order terms we finally obtain (because of Lemma A.2) Dh (u, u) & ∥∇h ucr ∥20,Ω +

∑ 1 ∥Π0 ([[u]])∥20,e & ∥u∥2DG 8he

e∈Eh

for all h ≤

B

√

α 16Ct2 +1

and the proof is complete.

Appendix

We provide the proof of Lemma 3.1 in this appendix. For simplicity we present the proof in the case of N = 2 subdomains, that is Ω = Ω1 ∪ Ω2 . The extension to the case of N subdomains is straightforward and we omit the details. We first show (24). Taking into account the definition (12) of Dh (·, ·), the linearity and symmetry of Dh (·, ·) and of the exact local solvers DiE (·, ·) (cf. (19)), it is easy to see that IhE (u, u) = Dh (u, u) − D1E (u1 , u1 ) − D2E (u2 , u2 ) = Dh (R1T u1 , R2T u2 ) + Dh (R2T u2 , R1T u1 ) = 2 Dh (RT1 u1 , RT2 u2 ) −1/2

−1/2

= 2 Ah (Bh RT1 u1 , Bh RT2 u2 ) ∑ ∫ −1/2 −1/2 ∇(Bh RT1 u1 ) · ∇(Bh RT2 u2 ) dx =2 T ∈Th

T

+2

∑ α ∫ −1/2 −1/2 Π0 ([[Bh RT1 u1 ]]) · Π0e ([[Bh RT2 u2 ]]) ds . h3e e e

e∈Eh

To give a more explicit expression of the last two terms on the right hand side, we take a closer look at the support of the terms involved. We first observe that supp(RT1 u1 ) ∩ supp(RT2 v2 ) ⊆ Γ, therefore it

28

−1/2

is enough to consider the action of Bh on e ∈ Γ. Fix an edge e ∈ Γ shared by the elements T1 ⊆ Ω1 and T2 ⊆ Ω2 , and recall that (see Figure 4(a)), { { u1 = u|Ω1 in Ω1 , 0 in Ω1 , T T R1 u1 = R2 u2 = 0 in Ω2 , u2 = u|Ω2 in Ω2 , −1/2

on internal edges (since e ∈ Γ so e ∈ Eh◦ ) we have   1 + βe [ ] u 1 u1   −1/2 −1/2 (Bh RT1 u1 )|e = Bh =  1 −2 β , 0 e u1 2   1 − βe ] [ u 2  0  −1/2 −1/2 (Bh RT2 u2 )|e = Bh =  1 +2 β , u2 e u2 2

Taking into account the action of Bh

(53)

−1/2

where βe is defined as in (8). Therefore, under the action of Bh , the support of RT1 u1 (resp. RT2 u2 ) expands into the Ω2 (resp. Ω1 ) across Γ, with an additional dof at the midpoint of the edge e (see −1/2 Figure 4(b)). Next, we have to further consider the actions of the operators ∇ and [[·]] on Bh RTi ui . −1/2 T R1 u1

(a) dofs of RT 1 u1

(b) dofs of Bh

−1/2 T R1 u1 )

−1/2 T R1 u1 ]]

(c) dofs of ∇(Bh

(d) dofs of [[Bh

−1/2

−1/2

−1/2

Figure 4: Degrees of freedom of RT1 u1 , Bh RT1 u1 , ∇(Bh RT1 u1 ) and [[Bh RT1 u1 ]], respectively, on a N = 2 subdomain partition. The dofs marked with • are different from zero; those marked with ◦ are equal to zero.

29

For the gradient term, it is clear that the resulting support expands along the strip of elements that touch Γ (see Figure 4(c)), that is: −1/2 T R1 u1 )

supp(∇Bh

−1/2 T R2 u2 )

∩ supp(∇Bh

⊆ ΩΓ ,

(54)

where the set ΩΓ is defined in (23). For the penalty term, it can be seen that (cf. Figure 4(d)), −1/2 T R1 u1 ]]))

supp(Π0e ([[Bh

−1/2 T R2 u2 ]]))

∩ supp(Π0e ([[Bh

⊆Γ.

(55)

Using the definition of the jump operator on interior edges and (53) we have −1/2 T R1 u1 ]]

[[Bh

−1/2 T R2 u2 ]]

= βe [[RT1 u1 ]] = βe u1 n1e ,

[[Bh

= βe [[RT2 u2 ]] = βe u2 n2e ,

(56)

and taking into account the definition (4) we obtain ∑ α −1/2 ∑ α ∫ −1/2 −1/2 −1/2 Π0e ([[Bh RT1 u1 ]]) · Π0e ([[Bh RT2 u2 ]]) ds = [[B RT1 u1 ]](me )[[Bh RT2 u2 ]](me ) 3 he e h2e h e∈Γ e∈Γ ∫ ∑ α ∑ α 2 + − 2 = (β u1 (me )u2 (me )ne · ne ) = − β Π0 (u1 ) · Π0e (u2 ) ds . h2e e h3e e e e e∈Γ

e∈Γ

Therefore, we finally have ] [ ∫ ∑ ∫ ∑ −1/2 T −1/2 T 0 0 E 2 α Ih (u, u) = 2 ∇(Bh R1 u1 ) · ∇(Bh R2 u2 ) dx − βe 3 Πe (u1 )Πe (u2 ) ds , he e T T ∈ΩΓ

e∈Γ

which establishes (26) and hence (24). We now turn to the case of inexact local solvers and the proof of (25). We first note that, when acting on (the restriction of the functions to) interior edges e ∈ Eh◦ that do not belong to the interface Γ, we −1/2 −1/2 have that Bh ≡ Bi . Hence, we can write: IhI (u, u) = Dh (u, u) − D1I (u1 , u1 ) − D2I (u2 , u2 ) = W1 + W2 , where W1 :=

∑ ∫ T ∈ΩΓ

−1/2

T

∇(Bh

−1/2

u) · ∇(Bh

u) dx −

2 ∑ ∫ ∑ i=1 T ∈ΩΓ T ⊂Ωi

−1/2

∇(Bi

−1/2

ui ) · ∇(Bi

ui ) dx,

T

and W2 :=

∑ α ∫ [ −1/2 −1/2 −1/2 −1/2 Π0 ([[Bh u]]) · Π0e ([[Bh u]]) − Π0e ([[B1 u1 ]]) · Π0e ([[B1 u1 ]]) h3e e e e∈Γ

−1/2

−Π0e ([[B2

−1/2

u2 ]]) · Π0e ([[B2

] u2 ]]) ds.

We first observe that the main difference with respect to the case of exact solvers is that the action −1/2 −1/2 of Bh on a function restricted to an edge e ∈ Γ differs from the action of the local operator Bi entering in the definition of DiI (·, ·). In the former case e ∈ Γ is an interior edge, while for the latter e is a boundary edge. In fact, in view of (11) we have ] [ ∂ ] [ u1 βe u1 −1/2 −1/2 (B1 u1 )|e = B1 = , 0 0 ] ] [ [ (57) 0 0 −1/2 −1/2 , (B2 u2 )|e = B2 = u2 βe∂ u2 30

−1/2

and therefore in this case the support of Bi W1 = 2

∑ ∫ T

T ∈ΩΓ



−1/2 T R1 u1 )

∇(Bh

ui remains in Ωi . For W1 , (54) together with (57) gives

−1/2 T R2 u2 ) dx

· ∇(Bh

 ∫ ∫ 2 ∑ ∑ ∑   −1/2 −1/2 −1/2 −1/2 + ∇(Bh RTi ui ) · ∇(Bh RTi ui ) dx − ∇(Bi ui ) · ∇(Bi ui ) dx . (58)  i=1

T ∈ΩΓ

T

T ∈ΩΓ T ⊂Ωi

T

For the term W2 , using (57) and (55) we have −1/2

[[Bi

ui ]] = βe∂ ui ,

where βe∂ is defined in (8). Hence, taking into account the above identity together with (56) and (55) we find, ) ]2 [ ]2 Π0e (u1 ) + Π0e (u2 ) − 2Π0e (u1 )Π0e (u2 ) , [ ]2 −1/2 −1/2 Π0e ([[Bi ui ]]) · Π0e ([[Bi ui ]]) = [βe∂ ]2 Π0e (ui ) −1/2

Π0e ([[Bh

−1/2

u]]) · Π0e ([[Bh

u]]) = [βe ]2

([

i = 1, 2 .

Thus we have

∫ ∫ ∑ α ( ∑ α ) ) ( 0 0 2 0 2 ∂ 2 [Πe (u1 )]2 + [Π0e (u2 )]2 ds W2 = − 2 [βe ] Πe (u1 )Πe (u2 ) ds + [βe ] − [βe ] 3 h3e h e e e∈Γ e e∈Γ ∫ ∑ α ∑ α ∫ ( ) =−2 [βe ]2 Π0e (u1 )Π0e (u2 ) ds − ηe2 [Π0e (u1 )]2 + [Π0e (u2 )]2 ds, (59) 3 3 he he e e e∈Γ

e∈Γ

where ηe is defined as in (29). Putting together (58) and (59) we finally obtain IhI (u, u)

=

IhE (u, u)

+

i=1

e∈Γ

 ∫ ∫ ∑  ∑ −1/2 −1/2 −1/2 −1/2 ∇(Bi ui ) · ∇(Bi ui ) dx , ∇(Bh RTi ui ) · ∇(Bh RTi ui ) dx −  

2 ∑

∑ α ∫ ( ) ηe2 [Π0e (u1 )]2 + [Π0e (u2 )]2 ds − 3 he e

T ∈ΩΓ

T

T ∈ΩΓ T ⊂Ωi

T

which is (27), and concludes the proof.

References [1] R. A. Adams. Sobolev spaces. Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1975. Pure and Applied Mathematics, Vol. 65. [2] P. F. Antonietti and B. Ayuso de Dios. Schwarz domain decomposition preconditioners for discontinuous Galerkin approximations of elliptic problems: non-overlapping case. M2AN Math. Model. Numer. Anal., 41(1):21–54, 2007. [3] P. F. Antonietti and B. Ayuso de Dios. Multiplicative Schwarz methods for discontinuous Galerkin approximations of elliptic problems. M2AN Math. Model. Numer. Anal., 42(3):443–469, 2008.

31

[4] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 39(5):1749–1779 (electronic), 2001/02. [5] B. Ayuso de Dios and L. T. Zikatanov. Uniformly convergent iterative methods for discontinuous Galerkin discretizations. J. Sci. Comput., 40(1–3):4–36, 2009. [6] A.T. Barker, S.C. Brenner, E.-H. Park, and L.-Y. Sung. Two-level additive Schwarz preconditioners for a weakly over-penalized symmetric interior penalty method. J. Sci. Comput., 47:27–49, 2011. [7] S.C. Brenner. Poincar´e-Friedrichs inequalities for piecewise H 1 functions. SIAM J. Numer. Anal., 41(1):306–324 (electronic), 2003. [8] S.C. Brenner, T. Gudi, L. Owens, and L.-Y. Sung. An intrinsically parallel finite element method. J. Sci. Comput., 42(1):118–121, 2010. [9] S.C. Brenner and L. Owens. A weakly over-penalized non-symmetric interior penalty method. J. Numer. Anal. Ind. Appl. Math., 2(1-2):35–48, 2007. [10] S.C. Brenner, L. Owens, and L.-Y. Sung. A weakly over-penalized symmetric interior penalty method. Electron. Trans. Numer. Anal., 30:107–127, 2008. [11] S.C. Brenner, L. Owens, and L.-Y. Sung. Higher order weakly over-penalized symmetric interior penalty methods. J. Comp. Appl. Math., 236:2883–2894, 2012. [12] S.C. Brenner and L.R. Scott. The Mathematical Theory of Finite Element Methods (Third Edition). Springer-Verlag, New York, 2008. [13] S.C. Brenner and K. Wang, Two -level additive Schwarz preconditioners for C0 interior penalty methods. Numer. Math. 102 (2005) 231-255. [14] M. Dryja and O. B. Widlund. Some domain decomposition algorithms for elliptic problems. In Iterative methods for large linear systems (Austin, TX, 1988), pages 273–291. Academic Press, Boston, MA, 1990. [15] M. Dryja and O. B. Widlund. Towards a unified theory of domain decomposition algorithms for elliptic problems. In Third International Symposium on Domain Decomposition Methods for Partial Differential Equations (Houston, TX, 1989), pages 3–21. SIAM, Philadelphia, PA, 1990. [16] M. Dryja, J. Galvis, and M. Sarkis. BDDC methods for discontinuous Galerkin discretization of elliptic problems. J. Complexity, 23(4-6):715-739, 2007. [17] M. Dryja, J. Galvis, and M. Sarkis. Neumann-neumann methods for a DG discretization of elliptic problems with discontinuous coefficients on geometrically nonconforming substructures. In Domain Decomposition methods in science and engineering XIX, pages 27–38, Lect. Notes Comput. Sci. Eng., Springer, Heidelberg, 2011. [18] M. Dryja and M. Sarkis Additive average Schwarz methods for discretization of elliptic problems with highly discontinuous coefficients. Comput. Methods Appl. Math., 10:164–176, 2010. [19] X. Feng and O. A. Karakashian. Two-level additive Schwarz methods for a discontinuous Galerkin approximation of second order elliptic problems. SIAM J. Numer. Anal., 39(4):1343–1365 (electronic), 2001. [20] C. Lasser and A. Toselli, An overlapping domain decomposition preconditioner for a class of discontinuous Galerkin approximations of advection-diffusion problems. Math. Comp. 72 (2003) 1215-1238 [21] A. Toselli and O. Widlund. Domain decomposition methods—algorithms and theory, volume 34 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2005.

32

[22] T. Warburton and J. S. Hesthaven. On the constants in hp-finite element trace inverse inequalities. Comput. Methods Appl. Mech. Engrg., 192(25):2765–2773, 2003.

33