• Keine Ergebnisse gefunden

Algebraic multilevel preconditioning in

N/A
N/A
Protected

Academic year: 2022

Aktie "Algebraic multilevel preconditioning in"

Copied!
26
0
0

Wird geladen.... (Jetzt Volltext ansehen)

Volltext

(1)

www.oeaw.ac.at

www.ricam.oeaw.ac.at

Algebraic multilevel preconditioning in

three-dimensional H(div) space

S. Tomar

RICAM-Report 2013-04

(2)

R IC A M re p o rt 2 0 1 2 -2 7 (v 3 ) an d 2 0 1 3 -0 4 (v 2 ), 8 Ju ly 2 0 1 3

S. K. TOMAR

Abstract. An algebraic multilevel iteration method for solving system of linear algebraic equations arising inHpcurlqandHpdivqspaces are presented. The algorithm is developed for the discrete problem obtained by using the space of lowest order Nedelec and Raviart-Thomas-Nedelec elements. The theoretical analysis of the method is based only on some algebraic sequences and generalized eigenvalues of local (element- wise) problems. In the hierarchical basis framework, explicit recursion formulae are derived to compute the element matrices and the constantγ(which measures the quality of the space splitting) at any given level.

It is proved that the proposed method is robust with respect to the problem parameters, and is of optimal order complexity. Supporting numerical results, including the case when the parameters have jumps, are also presented.

1. Introduction

Consider the finite element discretization of variational problems related to the bilinear form (1.1) Apu,vq:“αpu,vq `βpXu,Xvq, α, βPR`,

defined on the Hilbert space

(1.2) HpΩ,Xq:“ tv P pL2pΩqqd :XvPL2pΩqu.

HereΩ Ă Rd,d “ 2,3, is a Lipschitz domain, andXis the curl operator for d “ 2 and div operator for d “ 3. Note that div v “ Bxv1 ` Byv2 ` Bzv3 is the divergence of a three-dimensional vector v “ rv1,v2,v3sT, curl v “ Bxv2´ Byv1 is the scalar curl of a two-dimensional vector v “ rv1,v2sT, and p¨,¨qdenotes the inner-product in L2pΩq. Forα “ β “ 1, the bilinear form (1.1) is precisely the inner-product inHpΩ,Xq.

The adjoint of operatorXis defined by Xa

#curl forX“curl,d“2

´grad forX“div,d“3 ,

where, for a scalar functionw,gradw“ rBxw,Byw,BzwsT(for three-dimensional problem), andcurlw“ rByw,´BxwsT (for two-dimensional problem). Associated with the inner-productA, there exists a linear operatorA:“αI`βXaX, which mapsHpΩ,Xqonto its dual space, and is determined by the relation (1.3) pAu,vq “Apu,vq, @v PHpΩ,Xq,

Given a finite element space Vh ofHpΩ,Xq, the symmetric and positive-definite (SPD) operator Ah : Vh Ñ Vh, which is the discretization of the operatorAtogether with natural boundary conditions, is determined by

(1.4) pAhuh,vhq “Apuh,vhq, @vh PVh.

The operator equationAu“f, forf P pL2pΩqqd, then leads to the following discrete problem

(1.5) Ahuh“ fh,

which is uniquely solvable. ForHpΩ,curlq, such problems frequently occur in various contexts in elec- tromagnetism, e.g., low-frequency time-harmonic Maxwell equations [33], or some formulations of the (Navier -) Stokes equations [18], and forHpΩ,divqsuch problems frequently occur in, e.g., mixed formu- lations of elliptic problems, least-squares formulations of elliptic problems, part of fluid flow problems, and in functional-type a posteriori error estimates, see, e.g., [2,3,30,39,40] and the reference therein.

Therefore, developing fast solvers for large system of equations (1.5) is of significant importance.

Date: July 8, 2013.

1991Mathematics Subject Classification. 65N30, 65N22, 65N55.

Key words and phrases. Algebraic multilevel iteration method, lowest-order Nedelec and Raviart-Thomas-Nedelec spaces, optimal order complexity,HpdivqandHpcurlqspaces.

1

(3)

Preconditioning methods for such linear systems in Hpcurlq within the framework of domain de- composition methods, multigrid methods, and auxiliary space methods have been proposed by several authors, see e.g., [2, 3, 19, 23, 43, 44] and the references therein. The first results for multigrid in Hpdivq(based on smoothing and approximation property) was presented in [13] for triangular elements.

The first results for multigrid in Hpcurlq(within the framework of overlapping Schwarz methods) were obtained by Hiptmair in [20]. A unified treatment of multigrid methods for Hpcurlq and Hpdivq was presented by Hiptmair and Toselli in [21]. However, the condition number estimates of their precondi- tioned system were not robust with respect to the parameters αandβ. Arnold et al. [3] employed the multigrid framework by developing necessary estimates for mixed finite element methods (FEM) based on discretizations ofHpdivqandHpcurlq, and thereby obtained parameter independent condition number estimates of the preconditioned system. Pasciak and Zhao studied the overlapping Schwarz methods for Hpcurlqin polyhedral domains in [37], and Reitzinger and Schoeberl studied algebraic multigrid meth- ods for edge elements in [38]. Auxiliary space preconditioning, proposed by Xu in [46], was studied for H0pΩ,curlq (the space Hpcurlqwith zero tangential trace) by Hiptmair et al. [22]. Nodal auxiliary space preconditioning inHpcurlqandHpdivqwas studied by Hiptmair and Xu in [23], and the proposed preconditioner was robust with respect to the parametersαandβ.

The main principles in constructing efficient multigrid and multilevel solvers for (1.5) are projections into spaces of divergence-free vector fields, see [44], or, alternatively, a discrete version of the Helmholtz decomposition, see e.g., [2], and/or the construction of a proper auxiliary space, see e.g., [23]. Moreover, an effective error reduction generally demands to complement the coarse-grid correction by an appropri- ate smoother, e.g., additive or multiplicative Schwarz smoother, cf. [3]. The simple scalar (point-wise) smoothers, in general, do not work satisfactorily for this class of problems. All of these methods may be viewed as subspace correction methods [45,47], where different choices of specific components result in different methods (which also applies to the method presented in this paper).

Algebraic multilevel iteration (AMLI) methods were introduced by Axelsson and Vassilevski in a se- ries of papers [7,8,9,10]. The AMLI methods, which are recursive extensions of two-level methods for FEM [5], have been extensively analyzed in the context of conforming and nonconforming FEM (including discontinuous Galerkin methods), see [11,12,16,24,26,27,28,34,35,36]. For a detailed systematic exposition of AMLI methods, see the monographs [25,42]. These methods utilize a sequence of coarse-grid problems that are obtained from repeated application of a natural (and simple) hierarchi- cal basis transformation, which is computationally advantageous. The underlying technique of these methods often requires only a few minor adjustments (mainly two-level hierarchical basis transforma- tion) even if the underlying problem changes significantly. This is evident from the two different kind of problems considered in this paper, where the same algorithms (see Section4) are used. Furthermore, the AMLI methods are robust with respect to the jumps in the operator coefficients (where classical multigrid methods suffer), and are computationally advantageous than classical algebraic multigrid methods.

In this paper, we first derive the results for two-dimensional Hpcurlq problem. Note that, in two- dimensions, the lowest-order Nedelec space can be obtained by a 90 degrees rotation of lowest-order Raviart-Thomas space. Therefore, the space splitting presented in [28] also applies in this case (and vice-versa). However, we present a unified treatment of the element matrices arising from pu,vq and pXu,Xvq, which helps in deriving the explicit recursion formulae in simpler forms and without any undetermined constants. Moreover, with the unified treatment we are able to extend the results to three- dimensional lowest-order Raviart-Thomas-Nedelec elements in a straight-forward manner. Our analysis is based only on some algebraic sequences and the generalized eigenvalues of local (element-wise) prob- lems. In hierarchical setting, we derive explicit recursion formulae to compute the element matrices and the constantγ(which measures the quality of the space splitting) at any given level. The method is shown to be robust with respect to the parameters, i.e., the results hold uniformly for 0ăα, βă 8.

The remainder of this paper is organized as follows. In Section2we briefly discuss the finite element discretization of the model problem (1.1) using the lowest-order Nedelec and Raviart-Thomas-Nedelec spaces. Section3starts with a brief description of the AMLI procedure (in Section3.1). After presenting hierarchical basis transformations in Section 3.2, the construction of the hierarchical splitting of the lowest-order Nedelec and Raviart-Thomas-Nedelec spaces is presented in Section3.3. In Section3.4a local two- and multi-level analysis is then presented and the main result is proved. The algorithms used in this paper are provided in Section4. Finally, in Section 5we present numerical experiments. These

(4)

include the cases with known analytical solution (α“β“1), fixing one of the parameters and varying other from 10´6to 106, and the case of jumping coefficients. The conclusions are drawn in Section6.

2. Finite element discretization

In this section we briefly discuss the finite element discretization using lowest order Nedelec space in two-dimensions and lowest-order Raviart-Thomas-Nedelec space in three-dimensions, respectively.

2.1. Finite element discretization using Nedelec elements. We consider the tessellation of Ω Ă R2 using square elements, and choose the reference element ˆK asr´1,1s ˆ r´1,1s. LetPrx,rypKˆqdenote the space of polynomials of degree ď rx in x and ď ry iny. Also, let PrpBKˆq denote the space of polynomials of degreeď r on BK. For the construction ofˆ Vh, we use the space of lowest-order edge elements (Nedelec space of first kind), which is denoted by N0. The space N0pKˆqis defined as

N0pKˆq “P0,1pKˆq ˆP1,0pKˆq “

"

vpx,ˆ yˆq “

„ v1`v2yˆ v3`v4

*

(2.1) .

Thus, the local basis for N0has dimension 4. Moreover, forv0 PN0pKˆqwe have (2.2) curlv0PP0,0, v0¨t|BKˆ PP0pBKˆq,

where tdenotes the unit tangential vector to the element boundaries. For further details the reader is referred to, e.g., [33].

Now letF: ˆKÑR2be a diffeomorphism of the reference element ˆKonto a physical elementK, i.e., K“FpKˆq. ByJ we denote the Jacobian matrix of the mapping, and byJDits determinant, which are defined as

J “

ˆ Bxˆx Byˆx Bxˆy Byˆy

˙

, JD“ |detJ| “ BxˆxByˆy´ ByˆxBxˆyą0.

Then we have the following transformation relations:

w“J´Tw;ˆ curlw“JD´1curlw,ˆ @wPHpK,curlq,wˆPHpK,ˆ curlq. (2.3)

The vector transformationwÑJ´Twˆis called the covariant transformation, and curlw“JD´1curlwˆ is obtained via the well known Piola transformation wÑJD´1Jw.ˆ

We denote the element matrix forş

Ku¨v byLK, and forş

Kcurl ucurlv byCK. For the N0space based on uniform mesh composed of square elements, the element matricesLKandCKhave the following structure

LK “ 1 6

»

— –

2 1 0 0

1 2 0 0

0 0 2 1

0 0 1 2

fi ffi ffi fl

, CK “ 1

h2

»

— –

1 ´1 ´1 1

´1 1 1 ´1

´1 1 1 ´1

1 ´1 ´1 1

fi ffi ffi fl (2.4) .

The overall element matrixAK,C :“αLK`βCK, is thus given by

AK,C“ 1 6h2

»

— –

2αh2`6β αh2´6β ´6β 6β

αh2´6β 2αh2`6β 6β ´6β

´6β 6β 2αh2`6β αh2´6β

6β ´6β αh2´6β 2αh2`6β

fi ffi ffi fl (2.5) .

Lettinge“κh2, withκ“α{β, the element matrix can be written as

AK,C “ β 6h2

»

— –

2e`6 e´6 ´6 6

e´6 2e`6 6 ´6

´6 6 2e`6 e´6

6 ´6 e´6 2e`6

fi ffi ffi fl (2.6) .

Clearly, for allα, βPR`, and thusκPR`, we haveeą0. Note that for fixedκ, andhÑ0, the element matrixAK,Cis dominated by the matrixCK (which has a non-zero kernel), whereas for moderate values ofhit is a regular matrix. The near-nullspace of the matrixAK,C is given by the nullspace of the matrix CK, which is associated with the local bilinear formCKpu,vq :“ pcurlu,curl vqK.As we shall see in the analysis, the proposed method is of optimal order for all 0ăα, βă 8.

The following result can now be easily shown using [28, Lemma 2.1].

(5)

Lemma 2.1. (Near-nullspace of matrix AK,C). The element matrix AK,C given in (2.5) is symmetric positive definite (SPD). Moreover, the nullspace of the matrix CK for a general element K with nodal coordinatespxi,yiq, iP t1,2,3,4uis given by

(2.7) kerpCKq “spantp1,1,0,0qT,p0,0,1,1qT,px1,x2,y3,y4qTu.

Furthermore, in case of a uniform mesh composed of squareN0elements, the matrix CKis same for each element K and its nullspace is given by

kerpCKq “spantp1,1,0,0qT,p0,0,1,1qT,p´1,0,0,1qTu.

Remark 2.2. When using the lowest order Nedelec elements, the matrixCKis always of rank one. In the global assembly this yields a matrixCwhose rank equals the number of elements in the mesh. That is, the kernel of the global matrixChas dimension dimpkerpCqq “nE ´nK, wherenE denotes the number of faces andnK denotes the number of elements in the finite element mesh. Thereby, the dimension of the kernel is slightly more than half of the total number of degrees of freedom.

2.2. Finite element discretization using Raviart-Thomas-Nedelec elements. We consider the tessel- lation ofΩĂR3using cubic elements, and choose the reference element ˆK asr´1,1s3. LetPrx,ry,rzpKˆq denote the space of polynomials of degree ď rx in x, ď ry iny and ď rz inz, respectively. Also, let Pr1,r2pBKˆq denote the space of polynomials of degrees ď r1 and ď r2 in the respective dimensions on BK. For the construction ofˆ Vh, we use the space of lowest-order Raviart-Thomas-Nedelec elements, which is denoted by RTN0. The space RTN0pKˆqis defined as

RTN0pKˆq “P1,0,0pKˆq ˆP0,1,0pKˆq ˆP0,0,1pKˆq “

$

&

%

vpx,ˆ y,ˆ zˆq “

» –

v1`v2xˆ v3`v4yˆ v5`v6

fi fl

, . - (2.8) .

Thus, the local basis for RTN0has dimension 6. Moreover, forv0 PRTN0pKˆqwe have (2.9) divv0PP0,0,0, v0¨n|BKˆ PP0,0pBKˆq,

wherendenotes the unit normal vector to the element faces. For further details the reader is referred to, e.g., [14].

Now letF: ˆKÑR3be a diffeomorphism of the reference element ˆKonto a physical elementK, i.e., K“FpKˆq. ByJ we denote the Jacobian matrix of the mapping, and byJDits determinant, which are defined as

J “

¨

˝

Bxˆx Byˆx Bˆzx Bxˆy Byˆy Bzˆy Bxˆz Byˆz Bˆzz

˛

‚, JD“ |detJ| ą0.

Then we have the following relations:

w“JD´1Jw;ˆ divw“JD´1divw,ˆ @wPHpK,divq,wˆPHpK,ˆ divq, (2.10)

by the well known Piola transformation, see e.g., [14].

We denote the element matrix forş

Ku¨v byLK, and forş

Kdivudivv byDK. For the RTN0space based on uniform mesh composed of cubic elements, the element matricesLKandDKhave the following structure

LK “ 1 6h

»

— –

2 1 0 0 0 0

1 2 0 0 0 0

0 0 2 1 0 0

0 0 1 2 0 0

0 0 0 0 2 1

0 0 0 0 1 2

fi ffi ffi ffi ffi ffi ffi fl

, DK “ 1

h3

»

— –

1 ´1 1 ´1 1 ´1

´1 1 ´1 1 ´1 1

1 ´1 1 ´1 1 ´1

´1 1 ´1 1 ´1 1

1 ´1 1 ´1 1 ´1

´1 1 ´1 1 ´1 1 fi ffi ffi ffi ffi ffi ffi fl (2.11) .

The overall element matrixAK,D:“αLK`βDK, is thus given by

AK,D“ 1 6h3

»

— –

2αh2`6β αh2´6β 6β ´6β 6β ´6β

αh2´6β 2αh2`6β ´6β 6β ´6β 6β

6β ´6β 2αh2`6β αh2´6β 6β ´6β

´6β 6β αh2´6β 2αh2`6β ´6β 6β

6β ´6β 6β ´6β 2αh2`6β αh2´6β

´6β 6β ´6β 6β αh2´6β 2αh2`6β fi ffi ffi ffi ffi ffi ffi fl (2.12) .

(6)

With the definition ofeintroduced before (2.6), the element matrix can be written as

AK,D“ β 6h3

»

— –

2e`6 e´6 6 ´6 6 ´6

e´6 2e`6 ´6 6 ´6 6

6 ´6 2e`6 e´6 6 ´6

´6 6 e´6 2e`6 ´6 6

6 ´6 6 ´6 2e`6 e´6

´6 6 ´6 6 e´6 2e`6 fi ffi ffi ffi ffi ffi ffi fl (2.13) .

Note again that for fixedκ, andhÑ 0, the element matrixAK,Dis dominated by the matrixDK (which has a non-zero kernel), whereas for moderate values ofhit is a regular matrix. The near-nullspace of the matrixAK,Dis given by the nullspace of the matrixDK, which is associated with the local bilinear form DKpu,vq :“ pdivu,divvqK.As we shall see in the analysis, the proposed method is of optimal order for all 0ăα, βă 8.

Proposition 2.3. (Near-nullspace of matrixAK,D).The element matrix AK,Dgiven in(2.12)is symmetric positive definite (SPD). Furthermore, in case of a uniform mesh composed of cubicRTN0elements, the matrix DKis same for each element K and its nullspace is given by

kerpDKq

“spantp1,1,0,0,0,0qT,p´1,0,1,0,0,0qT,p1,0,0,1,0,0qT,p´1,0,0,0,1,0qT,p1,0,0,0,0,1qTu. Proof. Since the coefficientsαand β in (2.12) are positive, it follows from equation (1.5) thatAK,D is SPD for a general elementK. Moreover, for a uniform mesh composed of cubic RTN0elements, since the vectorp1,´1,1,´1,1,´1qT is orthogonal to the kernel ofDK, it is clear that the rank-one matrixDK is of the formc¨ p1,´1,1,´1,1,´1qT ¨ p1,´1,1,´1,1,´1q, for some constantc.

Remark 2.4. When using the lowest order Raviart-Thomas-Nedelec elements, the matrixDKis always of rank one. In the global assembly this yields a matrixDwhose rank equals the number of elements in the mesh. That is, the kernel of the global matrix Dhas dimension dimpkerpDqq “nF´nK, wherenF

denotes the number of faces andnKdenotes the number of elements in the finite element mesh. Thereby, the dimension of the kernel is slightly more than two-third of the total number of degrees of freedom.

3. Algebraic multilevel iteration

For the solution of the linear system arising from (1.5), we describe and analyze the AMLI method in the remainder of this section. Our presentation follows [28].

3.1. The AMLI procedure. In what follows we will denote byMpℓqa preconditioner for a finite element (stiffness) matrixApℓq corresponding to aℓtimes refined meshp0 ďℓ ď Lq. We will also make use of the corresponding ℓth level hierarchical matrix ˆApℓq, which is related toApℓq via a two-level hierarchical basis (HB) transformation Jpℓq, i.e.,

(3.1) Aˆpℓq “JpℓqApℓqpJpℓqqT.

The transformation matrix Jpℓqspecifies the space splitting, and will be described in detail in Section3.2.

ByApℓqi j and ˆApℓqi j , 1 ďi, jď2, we denote the blocks ofApℓq and ˆApℓqthat correspond to the fine-coarse partitioning of degrees of freedom (DOF) where the DOF associated with the coarse mesh are numbered last.

The aim is to build a multilevel preconditioner MpLqfor the coefficient matrixApLq:“Ahat the level of the finest mesh that has a uniformly bounded (relative) condition number

κpMpLq´1ApLqq “Op1q,

and an optimal computational complexity, that is, linear in the number of degrees of freedom NLat the finest mesh (grid). In order to achieve this goal hierarchical basis methods can be combined with various types of stabilization techniques. One particular purely algebraic stabilization technique is the so-called Algebraic Multi-Level Iteration (AMLI) method, which is presented below.

(7)

We have the following two-level hierarchical basis representation at levelℓ

(3.2) Aˆpℓq

«Aˆpℓq11pℓq12pℓq21pℓq22

«Aˆpℓq11pℓq12pℓq21 Apℓ´1q

ff .

Starting at level 0 (associated with the coarsest mesh), on which a complete LU factorization of the matrixAp0qis performed, we define

(3.3) Mp0q:“Ap0q.

Given the preconditioner Mpℓ´1qat levelℓ´1, the preconditioner Mpℓqat levelℓis then defined by

(3.4) Mpℓq:“LpℓqUpℓq,

where

(3.5) Lpℓq :“

« Cpℓq11 0 Aˆpℓq21 Cpℓq22

, Upℓq :“

» –

I C11pℓq´1pℓq12

0 I

fi fl.

HereC11pℓqis a preconditioner for the pivot blockApℓq11, and

(3.6) C22pℓq:“Apℓ´1q

´

I´ppℓqpMpℓ´1q´1Apℓ´1q´1

is an approximation to the Schur complement S “ Apℓ´1q ´Aˆpℓq21C11pℓq´1pℓq12, where Apℓ´1q “ Aˆpℓq22 is the stiffness matrix at the coarse levelℓ´1, and ppℓqis a certain stabilization polynomial of degreeν satisfying the condition

(3.7) 0ď ppℓqpxq ă1, @0ăxď1, and ppℓqp0q “1.

It is easily seen that (3.6) is equivalent to

(3.8) Cpℓq22´1“Mpℓ´1q´1qpℓqpApℓ´1qMpℓ´1q´1q, where the polynomialqpℓqpxqis given by

(3.9) qpℓqpxq “ 1´ppℓqpxq

x .

We note that the multilevel preconditioner defined via (3.4) is getting close to a two-level method when qpℓqpxqclosely approximates 1{x, in which caseCpℓq22´1 « Apℓ´1q´1. In order to construct an efficient multilevel method the action ofC22pℓq´1 on an arbitrary vector should be much cheaper to compute (in terms of the number of arithmetic operations) than the action ofApℓ´1q´1. Optimal order solution algo- rithms typically require that the arithmetic work for one application of Cpℓq22´1 is of the orderOpNℓ´1q whereNℓ´1denotes the number of unknowns at levelℓ´1.

To reduce the overall complexity of AMLI methods (to achieve optimal computational complexity), various stabilization techniques can be used. It is well known from the theory introduced in [7,8] that a properly shifted and scaled Chebyshev polynomial ppℓq :“ pν of degreeνcan be used to stabilize the condition number ofMpℓq´1Apℓq(and thus obtain optimal order computational complexity). Other poly- nomials such as the best polynomial approximation of 1{xin uniform norm also qualify for stabilization, see, e.g., [29]. This approach requires the computation of polynomial coefficients which depends on the bounds of the eigenvalues of the preconditioned system. Alternatively, a few inner flexible conjugate gradient (FCG) type iterations are performed at coarse levels to stabilize (or freeze the residual reduction factor of) the outer FCG iteration, which lead to parameter-free AMLI methods [9,10,24,34,35,36].

In general, the resulting nonlinear (variable step) multilevel preconditioning method is almost equally efficient as linear AMLI method, and, because its realization does not rely on any spectral bounds, it is easier to implement than the linear AMLI method (based on a stabilization polynomial). For a conver- gence analysis of nonlinear AMLI see, e.g., [24,25,42].

(8)

Typically, the iterative solution process is of optimal order of computational complexity if the degree ν “νof the matrix polynomial (or alternatively, the number of inner iterations for nonlinear AMLI) at levelℓsatisfies the optimality condition

1{ b

p1´γ2q ă νă τ, (3.10)

whereτ«τ“N{Nℓ´1denotes the reduction factor of the number of degrees of freedom (DOF), and γ denotes the constant in the strengthened Cauchy-Bunyakowski-Schwarz (CBS) inequality. In case of standard (full) coarsening, the value of τis approximately 4 for the sequence of N0 spaces, and 8 for the sequence of RTN0spaces. These sequences will be constructed in the next subsections. For a more detailed discussion of AMLI methods, including implementation issues see, e.g., [25,42].

Remark 3.1. The commonly used AMLI algorithm was originally introduced and studied in a multi- plicative form (3.4), see [7,8]. However, the preconditioner can also be constructed in the additive form, which is defined as follows [4,6,25]

(3.11) MApℓq:“

«

C11pℓq 0 0 Cpℓq22

ff .

In this case the optimal order of computational complexity demands that the matrix polynomial degree (or the number of inner iterations of nonlinear AMLI) satisfy the following relation

b

p1`γq{p1´γq ă νă τ.

(3.12)

3.2. Hierarchical basis forVh. The AMLI methods we are considering here, for the solution of (1.5), are based on a proper splitting of the spaceVh.

4

8

1 2

3 6 4

I

5 7

II

III

9 10

11 IV 12

1

3

2

Figure1. Macro-element obtained after one regular mesh-refinement step

For N0 subspace of Hpcurlq, the particular two-level HB transformation that induces this splitting was introduced in the context of linear nonconforming Crouzeix-Raviart (CR) elements in [11,12]. It was later studied for quadrilateral rotated bilinear (Rannacher-Turek) type elements in [16]. Note that the similarities of the HB transformation when using CR elements and Nedelec elements is due to the algebraic nature of the problem. For the discretization based on linear elements (for meshes consisting of triangles) or bilinear elements (for meshes consisting of squares), similar HB transformation matrix can be used. However, suitable changes will be required when working with meshes consisting of general quadrilaterals.

Consider two consecutive discretizations TH (coarse level) andTh(fine level). Figure1illustrates a macro-elementG(at fine level) obtained from a coarse element by one regular mesh-refinement step. Let ϕG “ tφipx,yqu12i“1be the macro-element vector of the nodal basis functions. Using the local numbering of DOF, as shown in Figure1(right picture), a macro-element level (local) transformation matrix JG is constructed based on differences and aggregates of each pair of basis functionsφiandφjthat correspond

(9)

to a macro element edge, i.e.,

(3.13) JG “ 1

2

»

— –

2 2

2 2

1 ´1

1 ´1

1 ´1

1 ´1

1 1

1 1

1 1

1 1

fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl .

Figure2. Macro-element obtained after one regular mesh-refinement step

For RTN0subspace ofHpdivq, the particular two-level HB transformation, that induces this splitting, was introduced in the context of Rannacher-Turek elements for three-dimensional elliptic problems [17].

Consider two consecutive discretizations TH (coarse level) and Th (fine level). Figure 2 illustrates a

(10)

macro-element G(at fine level) obtained from a coarse element by one regular mesh-refinement step.

The colors green, magenta and blue represent the face directions and face DOFs forx,yandzdirections, respectively. LetϕG “ tφipx,yqu36i“1 be the macro-element vector of the nodal basis functions. Using the local numbering of DOF, as shown in Figure2(second, third and fourth row of pictures), a macro- element level (local) transformation matrixJGis constructed based on differences and aggregates of basis functionsφiandφjthat correspond to a macro element face, i.e.,

(3.14a) JG “ 1

4

„ 4I JG;22

 , whereIis the 12ˆ12 identity matrix and

(3.14b) JG;22

»

— –

PD PD

PD PD

PD PD P1A P2A P3A P4A P5A P6A

fi ffi ffi ffi ffi ffi ffi ffi ffi fl .

Here each blockPiA,i“1,2, . . . ,6, which reflects the basis functions obtained by aggregates, is a 6ˆ4 matrix with all zeros except ith-row which has all ones. The block PD, which reflects the orthogonal transformation to aggregates, and obtained by suitable combination of differences, is given by

(3.14c) PD

» –

1 ´1 1 ´1

1 1 ´1 ´1

1 ´1 ´1 1

fi fl.

The transformations (3.13)-(3.14) define a two-level hierarchical basis ˆϕGlocally, namely, ˆϕG“JGϕG. 3.3. Hierarchical splitting. LetAGbe the macro-element stiffness matrix corresponding toG PT “ Th. The global stiffness matrixAhcan be written as

Ah“ ÿ

GPT

RGTAGRG,

where RG denotes the natural inclusion (canonical injection) of the matrix AG for all G in T. Note that the matrix AG is of size 12ˆ12 for two-dimensional Hpcurlq problem, and of size 36ˆ36 for three-dimensional Hpdivqproblem. Then the hierarchical two-level macro-element matrix is given by

G“ JGAGJGT,

and the related global two-level matrix can be obtained via assembling, i.e., ˆAh “ ř

GPT RTGGRG. Alternatively, one can compute the matrix ˆAhvia the triple matrix product

(3.15) Aˆh“JAhJT,

where the global transformation matrix Jis induced by the local transformations, i.e., J|G“JG, @GPT.

In other words, global and local transformations are compatible in the sense that restricting Jto the DOF of any macro-elementGwe obtainJG. Now, if we number those DOF that correspond to interior nodes of the macro elements first, the global two-level stiffness matrix ˆAhhas the 2ˆ2 block structure

(3.16) Aˆh

„ Aˆ1112

2122

 ,

where ˆA11 corresponds to theinterior unknowns. We follow the first reduce (FR) approach, see e.g., [11,12,16,17], where these interior unknowns are first eliminatedexactly. This static condensation step can be written in the form

(3.17) Aˆh

„ Aˆ11 0 Aˆ21 B

 „ I1´11112

0 I2

 ,

(11)

with the Schur complement B “Aˆ22´Aˆ21´11112. Next, the matrix Bis partitioned into 2ˆ2 blocks, i.e.,

(3.18) B“

„ B11 B12 B21 B22

 ,

where B11 and B22 correspond to the differences and aggregates of basis functions (associated with one macro-element edge or face), respectively. The matrix B22 at level ℓ then defines the coarse-grid matrixApℓ´1qin the AMLI hierarchy, cf. (3.2). This algorithm can be applied recursively on each level ℓ “L,L´1, . . . ,1. The resulting algorithm is then of optimal computational complexity, see e.g., [28, Remark 3.1].

3.4. Local analysis. In the two-level framework we denote byV1 and V2the subspaces of the finite element spaceVh. The spaceV2is spanned by the coarse-space basis functions (aggregates) andV1is the complement ofV2inVh, i.e.,Vhis a direct sum ofV1andV2:

(3.19) Vh“V1‘V2.

A measure for the quality of this splitting is the constantγin the strengthened CBS inequality, which is defined by the relation

γ “cospV1, V2q:“ sup uPV1, vPV2

Apu,vq aApu,uqApv,vq.

It is well known (see, e.g., [5]) that γ can be estimated locally over each macro elementG, and that γ “maxGγG,where

γG :“ sup

uPV1pGq, vPV2pGq

AGpu,vq aAGpu,uqAGpv,vq.

The spacesV1pGq,V2pGq, and the bilinear formAGpu,vqcorrespond to the restriction ofV1,V2, and Apu,vq, respectively, to the macro elementG.

We perform this local analysis on the matrix level, where the splitting (3.19) is obtained via the two- level hierarchical basis transformation described in Section 3.2, and the space Vh corresponds to the choice of lowest order Nedelec or Raviart-Thomas-Nedelec elements. In this setting the upper left block of ˆAhis block-diagonal. Note that, for two-dimensionalHpcurlqproblem, the diagonal blocks of ˆA11are of size 4ˆ4, which can be associated with the interior nodest1,2, . . . ,4uin the right picture of Figure1, and for three-dimensional Hpdivqproblem, the diagonal blocks of ˆA11 are of size 12ˆ12, which can be associated with the interior nodest1,2, . . . ,12uin the center column of second, third and fourth row of pictures in Figure 2. Therefore, we first compute the local Schur complements arising from static condensation of the interior DOF and obtain the matrices BG. Next we split each matrixBGas

BG

„BG,11 BG,12

BG,21 BG,22

 udifferences uaggregates ,

written again in two-by-two block form. For two-dimensional Hpcurlq problem, the block BG,11 and BG,22are both of size 4ˆ4, and for three-dimensionalHpdivqproblem the blockBG,11is of size 18ˆ18 and the block BG,22is of size 6ˆ6. We have thus reduced the problem of estimating the CBS constant of the splitting (3.19) to a small-sized local problem that involves the matrix BG. Following the general theory, see [5,15], to estimate the CBS constantγ, it suffices to compute the minimal eigenvalue of the generalized eigenproblem

(3.20) SGvG“λG,minBG,22vG, @vG,

where SG “BG,22´BG,21B´1G,11BG,12. The CBS constantγcan then be estimated as follows

(3.21) γ2ďmax

GPT γG2 “max

GPTp1´λG,minq.

Note that the matrix BG,11 is a well conditioned matrix, see Figure3, and therefore, it can be inverted cheaply, either by an iterative process or by, for example, an incomplete LUfactorization [41], which is denoted byBi11in Figure3.

We now first prove some auxiliary (stand-alone) results on algebraic sequences, which we will use to bound the CBS constantγ.

(12)

3 4 5 6 7 8 9 2

2.5 3 3.5 4 4.5 5

Number of levels

Condition number

κ(B11) κ((B11i )−1B11)

(a) Two-dimensionalHpcurlq

2 2.5 3 3.5 4 4.5 5 5.5 6

4 4.5 5 5.5 6 6.5 7 7.5 8

Number of levels

Condition number

κ(B11) κ((B11i )−1B11)

(b) Three-dimensionalHpdivq

Figure3. Condition number of the matrixBG,11 Lemma 3.2. For all eą0, consider the coupled sequences

b0“e´6, a0 “2e`6“2pb0`9q, (3.22a)

bℓ`1“ ´b2{a, aℓ`1“2a`bℓ`1, ℓ“0,1,2, . . . . (3.22b)

Let r“b{a. Then we have

bℓ`1{a “ ´r2, aℓ`1{a“2´r2, rℓ`1“ ´r2{p2´r2q, (3.23a)

aℓ`1´bℓ`1 “2a, aℓ`1`bℓ`1“ 2

apa2 ´b2q “2ap1´r2q, aℓ`1`bℓ`1

aℓ`1´bℓ`1 “1´r2. (3.23b)

Moreover, the following bound holds

´1ăr0ă1{2, and ´1ăr ď0 @ℓ“1,2, . . . , (3.24a)

aą. . .a1 ąa0ą6, 0ďr2ď. . .ďr21 ďr02ă1,@ℓ“0,1,2, . . . . (3.24b)

Proof. Using the definition ofrin (3.22b), we getbℓ`1{a “ ´r2, and thusaℓ`1{a “2´r2. The last relation of (3.23a) then immediately follows. The relations (3.23b) are also easily obtained from (3.22b) and (3.23a).

Clearly, fore ą 0, we havea0 ą 6, and sincer0 “b0{a0 “ pe´6q{p2e`6q, it is easy to see that

´1 ă r0 ă 1{2. The latter also implies that 0 ď r20 ă 1. We now prove the remaining bounds using induction.

ℓ“0. Sincea1{a0 “2´r02 ą1, we havea1 ąa0 ą6. Moreover,r1 “ ´r20{p2´r20q. This implies that´1ăr1ď0, and thus 0ďr21 ă1. Furthermore, whenr0,0, we have

r21

˜ ´r02 2´r20

¸2

ñ r21

r20 “ r20

p2´r02q2 ă1.

And, sincer1 “0 ifr0“0, we haver12ďr20ă1.

ℓ“n. Assume that the relations (3.24) hold for ℓ “ n. Since an`1{an “ 2 ´rn2 ą 1, we have an`1 ąan ą6. Moreover,rn`1 “ ´r2n{p2´r2nq. This implies that ´1 ărn`1 ď 0, and thus 0ďrn`12 ă1. Also, whenrn,0, we have

r2n`1

ˆ ´rn2 2´r2n

˙2

ñ rn`12

r2n “ rn2

p2´r2nq2 ă1.

And, sincern`1 “0 ifrn “0, we havern`12 ďr2n ă1.

This concludes the proof.

(13)

Lemma 3.3. Let eą0and the sequences aand b be as defined in Lemma3.2. Then for

(3.25) c2ℓ,C “ 36pa`bq

pa2 ´36qpa´bq, the following bounds hold for allℓ“0,1,2, . . .

(3.26) c2ℓ,C ăc2ℓ´1,C ă. . .ăc21,C ăc20,Că3{8.

Proof. Froma0 “2e`6 andb0“e´6, we havea0´b0 “e`12,a0`b0 “3e,a0´6“2e, and a0`6“2pe`6q. Substituting these relations in the definition ofc20,C, we get

(3.27) c20,C “ 27

pe`6qpe`12q ă3{8.

Now

c21,C´c20,C “ 36`

pa1`b1qpa20´36qpa0´b0q ´ pa0`b0qpa21´36qpa1´b1q˘ pa21´36qpa1´b1qpa20´36qpa0´b0q .

Substituting the values ofa0,a1,b0andb1, and after some lengthy, but simple calculations, we find that c21,C´c20,C “ 108e`

´9e2p312`80e`5e2

pe`3qpa21´36qpa1´b1qpa20´36qpa0´b0q. Since the denominator is a positive quantity, we getc21,C´c20,C ă0, and thus

(3.28) c21,C ă3{8.

For remaining bounds, we again use induction. Note that, using (3.23b) we get (3.29) c2ℓ`1,C“ 36paℓ`1`bℓ`1q

pa2ℓ`1´36qpaℓ`1´bℓ`1q “ 36p1´r2q pa2ℓ`1´36q. Therefore, to show thatc2ℓ`1,Că3{8, it suffices to show that

(3.30) a2ℓ`1´36ą96p1´r2q.

Sincec21,C ă3{8, we clearly havea21´36ą96p1´r20q. Now assume that the relation (3.30) holds for ℓ“n´1, i.e.,

(3.31) a2n´36ą96p1´r2n´1q.

Multiplying (3.31) byp2´r2nq2and subtracting 36 from both sides we get

p2´r2nq2a2n´36ą36p2´rn2q2`96p1´r2n´1qp2´r2nq2´36 ñ a2n`1´36ą96`

p2´r2nq2p11{8´r2n´1q ´3{8˘ . We need to show thatp2´rn2q2p11{8´r2n´1q ´3{8ą1´rn2, i.e.,

gn :“ p2´rn2q2p11{8´r2n´1q `r2n´11{8ą0.

(3.32)

From the recurrence relation onrnfrom (3.23a), we have rn2“ rn´14

p2´r2n´1q2, 2´rn2“ pr4n´1´8r2n´1`8q p2´r2n´1q2 . Substituting these relations ingn, and after some lengthy calculations we obtain (3.33) gn “ p1´r2n´1q2

p2´r2n´1q4p´r6n´1`15r4n´1´64r2n´1`66q. Now forrn´12 P r0,1q, we have

1´r2n´1ą0, 2´rn´12 ą0, 66´64r2n´1ą0, 15rn´14 ´r6n´1 ě0,

which proves thatgn ą 0, and thata2n`1´36 ą96p1´r2nq. Therefore, the inequality (3.30) holds for allℓ“0,1, . . ..

(14)

To prove the monotonicity ofc2ℓ,C, we show that

(3.34) f:“c2ℓ`1,C{c2ℓ,Că1.

Using (3.29) we get

f “ p1´r2qpa2 ´36q p1´r2ℓ´1qpa2ℓ`1´36q. Multiplying numerator and denominator byp2´r2q2, we obtain

f“ p1´r2q`

p2´r2q2a2 ´36p2´r2q2˘ p1´rℓ´12 qpa2ℓ`1´36qp2´r2q2

“ p1´r2q p1´r2ℓ´1q

´

a2ℓ`1´36`36p1´ p2´r2q2q¯ pa2ℓ`1´36qp2´r2q2

“ p1´r2q

p1´r2ℓ´1qp2´r2q2 ` 36p1´r2qp1´ p2´r2q2q p1´r2ℓ´1qpa2ℓ`1´36qp2´r2q2. Now sincec2ℓ`1,C ă3{8, we havep1´r2q{pa2ℓ`1´36q ă1{96 from (3.30). Therefore,

fă p1´r2q

p1´r2ℓ´1qp2´r2q2 ` 36p1´ p2´r2q2q 96p1´r2ℓ´1qp2´r2q2

“ p1´r2q `3

8p1´ p2´r2q2q p1´r2ℓ´1qp2´r2q2

11{8´r2´3

8p2´r2q2 p1´rℓ´12 qp2´r2q2 . This gives

f´1ă

11{8´r2 ´3

8p2´r2q2´ p1´r2ℓ´1qp2´r2q2 p1´r2ℓ´1qp2´r2q2

“ 11{8´r2 ` p2´r2q2p´11{8`r2ℓ´1q p1´rℓ´12 qp2´r2q2 . Using (3.32) we therefore get

f´1ă ´g

p1´r2ℓ´1qp2´r2q2 ă0,

sinceg ą0, 1´r2ℓ´1ą0, andp2´r2q2ą0. This proves (3.34) and concludes the proof.

Lemma 3.4. Let eą0and the sequences aand b be as defined in Lemma3.2. Then for

(3.35) c2ℓ,D“ 72pa`bq

pa`12qpa´6qpa´bq, the following bounds hold for allℓ“0,1,2, . . .

(3.36) c2ℓ,Dăc2ℓ´1,Dă. . .ăc21,Dăc20,Dă1{2.

Proof. Substituting the relations fora0,b0,a0´b0,a0`b0,a0´6, anda0`6 from Lemma3.3in the definition ofc20,D, we get

(3.37) c20,D“ 54

pe`9qpe`12q ă1{2.

Now substituting the values ofa0,a1,b0andb1, and after some lengthy, but simple calculations, we find that

(3.38) c21,D´c20,D“ ´486ep5e2`88e`372q

pe`9qpe`12qp7e`48qp7e2`84e`108q ă0.

(15)

For remaining bounds, we use induction and proceed as follows. Let tm :“ 1{2´c2m,D and tm`1 :“ 1{2´c2m`1,D. Then, expandingam`1andbm`1in terms ofamandbm, and dropping the subscripts ofam

andbmfor brevity reasons, we get

tm:“1{2´c2m,D“ ´216a`6a2`a3´72b´6ab´a2b 2pa´6qpa`12qpa´bq “: nm

dm

(3.39a) ,

tm`1:“1{2´c2m`1,D“ ´216a2`12a3`4a4`144b2´6ab2´4a2b2`b4

2p´6a`2a2´b2qp12a`2a2´b2q “:nm`1 dm`1

(3.39b) ,

wherenmandnm`1are the numerators oftmandtm`1, respectively, anddmanddm`1are the denominators oftmandtm`1, respectively. Assume that the relation (3.36) holds forℓ“mě1, i.e.,tm ą0. We need to show thattm`1ą0. Sinceaą6,aą |b|, andbă0 formě1, we see thatdmanddm`1are positive.

Therefore, it suffices to show thatnm`1is positive whenevernmis positive. Givena{2ą1, we consider nm`1´a

2nm. We have 2pnm`1´a

2nmq “ pa`bqp7a3`18a2´6a2b`288b`2b3´216a´12ab´2ab2q

“ pa`bqp´6bpa2`2a´48q `3apa2`6a´72q `2pa3`b3q `2apa2´b2qq ą0,

sinceaą6,aą |b|, andbă0 formě1. This proves thatnm`1ą0, and hence,tm`1ą0.

The monotonicity ofc2ℓ,Dcan be shown by using (3.38) and showing the induction thatc2m`1,D´c2m,Dă 0 wheneverc2m,D´c2m´1,D ă 0. The details are omitted here (the results can also be verified by using algebraic cylindrical decomposition in a computer algebra system like Mathematica [31]).

The sequences a, b, and r are plotted in Figure 4, and the sequences c2ℓ,C and c2ℓ,D are plotted in Figure5.

0 5 10 15 20 25 l

100 000 200 000 300 000 400 000 a

0 5 10 15 20 25 l

-15 -10 0 5 b

0 5 10 15 20 25 l

-0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 r

Figure4. a,b, andrfore“10m0, wherem0“ t2,1,0, . . .´11,´12u(left to right)

0 5 10 15 20 25 l

0.05 0.10 0.15 0.20 0.25 0.30 0.35 c2

(a)c2ℓ,C

0 5 10 15 20 25 l

0.1 0.2 0.3 0.4 0.5 c2

(b)c2ℓ,D

Figure5. c2 fore“10m0, wherem0 “ t2,1,0, . . .´11,´12u(left to right)

We are now in a position to prove the following theorem which provides a theoretical estimate that holds on all levels of recursive splitting of the N0subspace ofHpcurlq, and RTN0subspace ofHpdivq.

Referenzen

ÄHNLICHE DOKUMENTE

At the same time, the standard Tikhonov method with a regularization parameter chosen in accordance with the balancing principle (2.5) implemented for ρ(u, v) = ku − vk L 2 does

For two dimensional simple manifolds with boundary, injectivity with in the solenoidal tensor fields has been establish fairly recent: in the non-attenuated case for 0- and 1-tensors

We consider an algebraic multilevel preconditioning method for SPD matrices resulting from finite element discretization of elliptic PDEs.. In particu- lar, we focus on

Moreover, numerical tests with a nonlinear algebraic multilevel iteration (AMLI) method demonstrate that the presented two-level method can be applied successfully in the recursive

The following remarks are taken from the literature and concern the Cauchy problem for linear systems of partial differential equations with constant coefficients for function

An optimal order algebraic multilevel iterative (AMLI) method for solving systems of linear alge- braic equations arising from the finite element discretization of certain

The criteria for the inclusion of participants and the making of decisions in EGAs are not generally compatible with the conventional norms for democratic legitimation used within

AWBET Cross-border shareholders and participations – transactions [email protected] AWBES Cross-border shareholders and participations – stocks