Preconditioners for higher order finite element discretizations of H(div)-elliptic problem

Preconditioners for higher order finite element discretizations of H(div)-elliptic problem

Junxian Wang; Liuqiang Zhong; Shi Shu

School of Mathematical and Computational Sciences, Xiangtan University, Hunan 411105, China. E-mails: [e-mail] / [e-mail] / [e-mail]

ABSTRACT

In this paper, we are concerned with the fast solvers for higher order finite element discretizations of H(div)-elliptic problem. We present the preconditioners for the first family and second family of higher order divergence conforming element equations, respectively. By combining the stable decompositions of two kinds of finite element spaces with the abstract theory of auxiliary space preconditioning, we prove that the corresponding condition numbers of our preconditioners are uniformly bounded on quasi-uniform grids.

Mathematical subject classification: Primary: 65F10; Secondary: 65N22.

Key words: preconditioner, higher order finite element, stable decomposition, H(div)-elliptic problem.

1 Introduction

Let Ω be a simply connected polyhedron in

³ with boundary Γ and unit outward normal ν. We define the Hilbert spaces H₀(div; Ω) as follows

H₀ (div; Ω) = {u ∈ (L² (Ω))³ | ∇ . u∈L² (Ω), ν . u = 0 on Γ}

with the inner product

(u, ν)_div = (u, ν) + (∇. u, ∇.ν),

where (·, ·) denotes the inner product in (L²(Ω))³ or L²(Ω).

In this paper, we consider the following variational problem: Find u∈ H₀(div; Ω) such that

where ƒ ∈ H₀(div; Ω)' is a given data and

with the constant τ > 0.

The bilinear form a(·, ·)induces the energy norm

Variational problem of the form (1) arises in numerous problems of practical import. Typical examples include the mixed method for second order elliptic problems, the least squares method of the form discussed in [3], and the sequential regularization method for the time dependent Navier-Stokes equation discussed in [6]. For a more detailed discussion of applications, we refer to [1].

To avoid the repeated use of generic but unspecified constants, following [9], we will use the following short notation: x y means x < Cy, x y means x > cy, and x ≈ y means cx < y < Cy, where c and C are generic positiveconstants independent of the variables that appear in the inequalities and especially the mesh parameters.

Outline. The remainder of this article is organized as follows. In the next section, we introduce two kinds of higher order finite element equations, and present the corresponding frame of constructing preconditioner. We construct the preconditioners for two kinds of higher order divergence conforming element equations, and prove that their corresponding condition number is uniformly bounded in Section 3 and Section 4, respectively.

2 Finite element equations and framework of preconditioner

Let

_h be a shape regular tetrahedron meshes of Ω, where h is the maximum diameter of the tetrahedra in _h. Now, we present two families of divergence conforming finite elements spaces (see [7])

where

_k denote the standard space of polynomials of total degree less than or equal to k, and _k denote the space of homogeneous polynomials of order k.

We consider the solution of systems of linear algebraic equations which arise from the finite element discretization of variational problems (1): Find (k > 1, l = 1, 2) such that

Their algebraic systems can be described as

Since is symmetric positive definite, we use precondition conjugate gradient (PCG) methods to solve algebraic systems (5). In this paper, we will construct the preconditioners for the cases of higher order finite equations, and present some estimates of the corresponding condition numbers.

For this purpose, we need to introduce some auxiliary spaces and corresponding operators.

Let V = with inner product a(·, ·) given by (2).

Let

₁, ..., _J, J ∈ , be Hilbert spaces endowed with inner products _j(., .), j = 1, ..., J. The operators _j: _j

are isomorphisms induced by _j(., .), namely

_j(_j, _j) = < _j

_j, _j > ∀_j, _j∈_j,

here we tag dual spaces by ' and use angle brackets for duality pairings. For each

_j, there exist continuous transfer operators Π_j : _j

V. Then we can construct the preconditioner for operator as follows:

where

_j :

_j are given preconditioners for _j, and are adjoint operators of Π_j.

Now, we present the following theorem of an estimate for the spectral condition number of the preconditioner given by (6).

Theorem 2.1. Assume that there exist constants c_j, such that

and for ∀ u ∈ V, there exist _j∈ _j such that u = and

then for the preconditioner B given by (6), we have the following estimate for the spectral condition number

Proof. We define the space

= ₁ × ₂ × ... × _j

with the inner product

and the following two operators

Π = (Π₁,Π₂, ..., Π_j) :

V,

= diag (₁, ₂, ..., _j) :

,

= diag (₁, ₂, ...,_j) :

.

Thus we can rewrite the definition of operator B given by (6):

B = Π Π*.

Using the definitions of inner product in , operators Π and , and conditions (7)-(8), then there exists a constant := , such that

and for ∀ u ∈ V, there exists ∈ , such that u = Π and

||||< c₀||u||_A

From Corollary 2.3 of [5], we immediately get an estimate for the spectral condition number of the preconditioned operator B

The desired estimates then follow by combining the above inequality and the following fact

□

The principal challenge confronted in the development of preconditioners by applying Theorem 2.1 is to construct some appropriate spaces and operators which satisfy (7) and (8). In the following two sections, we present the corresponding spaces and operators for two kinds of divergence conforming element spaces, respectively.

3 Preconditioner for finite element equations of first kind

We first introduce Sobolev functional space

H₀(curl; Ω) = {u ∈ (L²(Ω))³|∇ × u ∈ (L²(Ω))³, ν × u = 0 on Γ}

with the norm

There exist two families of edge finite element spaces for the space H₀(curl; Ω) (see [2, 4, 7]).

1. k order Nédélec element of first kind:

where

_k = (_k-1)³ ⊕ {p ∈ (_k)³ | p (x) · x = 0}

2. k order Nédélec element of second kind:

We also need to introduce the following space of piecewise k–degree discontinuous scalar elements on _h:

The Sobolev spaces H₀(div; Ω), H₀(curl; Ω) and the corresponding finite element spaces possess the exceptional exact sequence properties (see [4, 7])

Assuming that u has the necessary smoothness, we can define two kinds of interpolants: and , such that u ∈and u ∈(more details refer to [4, 7]). Especially, the interpolation is not defined for a general function in H₀(div; Ω). Here let us quote a slightly simplified version (see Theorem 5.25 of [7]).

Lemma 3.1. Suppose that there are constants δ > 0 such that u ∈ (H^1/2+δ(K))³ for each K in _h. Then u is well-defined, and we have

with a constant only depending on the shape regularity of

_h.

The finite element spaces is equipped with bases (k, 1) comprising locally supported functions. These bases are L² stable in the sense that

with constant only depending on the shape-regularity of

_h.

Lemma 3.2. The interpolation operator is bounded on ((Ω))³ and satisfies

with a constant only depending on the shape regularity of

_h.

Furthermore, all above operators possess the following commuting diagram property (see [7])

We may apply the quasi-interpolation operators for Lagrangian finite element space introduced in [8] to the components of vector fields separately. This gives rise to the projectors

_h:((Ω)) ³

, which inherits the continuity

and satisfies the local projection error esitmate

Now, we present the stable decomposition of , k > 2.

Lemma 3.3. For any ∈ , there exist Span{b}, ∈ , such that

and

where the constant

₀ only depends on Ω and the shape regularity of _h.

Proof. For any given ∈ , using the continuous Helmholtz decomposition, there exist Ψ∈ ((Ω)) ³, p ∈ H₀(curl; Ω) such that

and

with constants only depending on Ω.

Taking the div of both sides of (23) and using (14), we get

Owing to Lemma

Ψ is well defined. Furthermore, the commuting diagram property (18) implies

This confirms that the third term in the splitting

actually belongs to the kernel of div. By (12), then there esists q ∈ H₀(curl; Ω) such that

Noting that

_h Ψ∈⊂ , which leads to

Substituting (25), (26) and (27) into (23), we have

Since , (Id – _h ) Ψ ,_hΨ ∈ , we obtain ∇ × (q + p) ∈ (div0) by using (28), then observing (13), there exists q_h ∈ , such that

Let

It's easy to obtain ∈ by noting that _hΨ∈ ⊂ and ∇ × q_h ∈∇ × ⊂ . Substituting (29), (30) and (31) into (28), we conclude

which completes the proof of (21).

Using (30), triangular inequality, Lemma 3.2, (20) and (24), we have

which leads to

It follows readily from inverse estimate and (16) that

Using inverse estimate again yields

By means of (33) and inverse estimate, we get

In view of (32), triangular inequality (34), (35) and (36), we have

which completes the proof of (22). □

We rely on the stable decomposition for V = in Lemma 3.3 and apply the abstract theory in Section 2 to define the preconditioner for finite element equations of first kind.

Let V = and choose two auxiliary spaces and the corresponding transfer operators as follows.

1. = , with inner product ₁(·, ·) which is defined by

where

The transfer operator is Π₁ = Id.

2.

₂ = with inner product ₂(·, ·) = a(·, ·) in the sense that

which concludes that

₂ = . The transfer operator is Π₂ = Id.

Making use of (6), the auxiliary space preconditioner for reads

where is the preconditioner of , ₁ is the preconditioners of ₁.

Noting that

₁ denotes the diagonal matrix of , in the practical application, we will take ₁ as the Jacobi (or Gauss-Seidel) smoothing operator for . Obviously, this special choose satisfies

where the constant

₁ is independent of the mesh parameters.

First, we prove that the above transfer operators satisfy the condition (7).

Due to the definitions of inner product and transfer operator in space

₁, for any given , where α_b∈ , we have

where the constant M bounds the number of basis functions whose support overlaps with a single element K.

For any given

₂∈ ₂, it's easy to obain

Combining (39) with (40), we conclude that (7) holds with the constants c₁ = M and c₂ = 1.

Secondly, the above spaces and operators satisfy the condition (8) by using the Lemma 3.3.

Summing up, we obtain the following theorem by using Theorem 2.1.

Theorem 3.4. Forgiven by (37), and ₁satisfies the condition of (38), then we have

with a constant only depending on the constants₀,₁and the shape regularity of _h.

4 Preconditioner for finite element equations of second kind

Now, we present the another stable decomposition of with k > 2.

Lemma 4.1. For any ∈ , there are ∈ and φ_h ∈ such that

and

where the constant c₀only depends on Ω and the shape regularity of _h.

Proof. For any ∈ , we can interpolate by Lemma 3.1. Thus, using (18), we have

In view of (14), we have

Making use of (45) and noting that = Id in (44), we get

namely

Noting that - ∈ , then by (46) and (13), there exists φ_h∈ , such that

where = , which completes the proof of (42).

Using (47), (15) with δ = 1/2, and the inverse estimate, we obtain

Squaring and summing over all the elements, we get

In view of (3) and (48), we find

Making use of (47), triangular inequality and (48), we have

A direct manipulation of (47) gives that

A combination of (49), (50) and (51) concludes (43). □

In this case, let V = . We choose the following two auxiliary spaces and the corresponding transfer operator.

1.

₁ = with inner product ₁(., .) = a(., .) in the sense that

which concludes that

₁ = . The corresponding transfer operator is Π₁ = Id.

2.

₂ = with inner product

The corresponding transfer operator is Π₂ = curl.

Then by using (6), we obtain the auxiliary space preconditioner for as follows

where is the preconditioner of , and ₂ is the preconditioners of ₂ given by (52).

Especially, we adopt the preconditioner

₂ in [10], this choice satisfy

where the constant C₁ is independent of the mesh parameters.

It is easy to prove that the above transfer operators satisfy the conditions (7). In fact, using the definitions of inner products and transfer operators in spaces

_l (l = 1,2), we have

namely, the conditions (7) of Theorem 2.1 hold with the constants c₁ = c₂ = 1.

Applying Theorem 2.1 and using Lemma 4.1, we have the following Theorem.

Theorem 4.2. Forgiven by (53), and ₂satisfies the condition of (54), then we have

with a constant only depending on the constants c₀and C₁and the shape regularity of _h.

Combining Theorem 3.4 and Theorem 4.2, by using a Jacobi (or Gauss-Seidel) smoothing, we can translate the construction of preconditioner for into the one of. Furthermore, by using the preconditioner of H(curl; Ω)-elliptic problem, we can translate the preconditioner for into the one for . Since Hiptmair and Xu [5] have constructed an efficient preconditioner for , we construct the efficient precondtioners for (k = 1, l = 2 or k > 2, l = 1,2) and prove the corresponding spectral condition numbers are uniformly bounded and independent of mesh size h and the parameter τ by this recursive form.

5 Implementation of algorithm and numerical experiments

For simplicity, we only give the description of the preconditioning algorithm defined by (53) when k = 2.

Note that when k = 2, (53) turn to

In the following, we first discuss the description of algorithm about the preconditioner . For this purpose, we introduce the following operators

and

then, the algorithm about the operator can be described by (see [5] formore details)

Algorithm 5.1. For a given g ∈ , then u_g = g ∈ can be obtained as follows:

Step 1: Applying m ₁ times symmetric Gauss_Seidel iterations in variational problem

with a zero initial guess to get₁, whereƒ = g.

Step 2: Computing₂∈ by

Step 3:Computing₃∈ by

which can be obtained by

1. Applying m₂times symmetric Gauss_Seidel iterations in (59) with a zero initial guess to get ₄.

2. Computing ₅ ∈ by

3. Set ₃ = ₄+ ()^T

₅.

Step 4:Set u_g = ₁+()^T

₂+()^T

₃.

By [5], the preconditioner defined by Algorithm 5.1 satisfy

K () < C₁

where the constant C₁ is independent of the mesh size h and parameter τ.

Next, we give the description of algorithm for the operator curl

₂curl^*. Firstly, let

n = dim (), dim (),

and

then we introduce the transfer matrix(or operator)

By using , we can define the following matrix(or operator)

In view of (4.1) in [10], we can construct the preconditioner

₂ for , and its spectral condition number satisfy

where the constant C₂ is independent of the mesh size h and parameter τ.

Noting that the operator

₂ can be divided into three parts: the first part is to use the Jacobi (or Gauss-Seidel) smoothing for (52) in space , the second part is to solve the restriction of (52) in , the third part is to solve the restriction of (52) in ∇ . We can drop the second and third parts by using the fact that the second part is the same as (60) and curl grad 0. Hence the operator curl

₂curl^* can be simplified.

Summing up, we can obtain the following algorithm of the preconditioner .

Algorithm 5.2.For g ∈ , the solution u_g = g ∈ can be gotten as follows:

Step 1:Computing u₁∈ by Algorithm 5.1.

Step 2:Applying m₃times symmetric Gauss_Seidel iterations to get u₂∈ V^2,1by

( u₂, v₂) = (g, ∇ x v₂), ∀ v₂∈ V^2,1.

Step 3: Set

u_g = u₁+ u₂

For variational problem (4), we apply Algorithm 5.2 to the following two examples:

Example 5.1.The computational domain is Ω = [0,1] × [0, 1] × [0, 1] and the corresponding structured grids can be seen in Figure 1. For the convenience of computing the exact errors, we construct an exact solution u = (u₁, u₂, u₃) as

Example 5.2. The computational domain is the spheres of radius 1 and the corresponding unstructured grids can be seen in Figure 2, the exact solution u = (u₁, u₂, u₃) is

Now, we present some numerical experiments with m₁ = m₂ = m₃ = 3.

Table 1 gives the L₂ and H(div) error estimates for Example 5.1 when τ = 1, which shows that is the optimal convergence.

The condition number estimates and iteration counts for Example 5.1 and Example 5.2 are listed in Tables 2 - 5 for different values of the mesh size h and the scaling parameter τ. By these Tables, we find that the condition number and iteration counts are independent of the mesh size h and weakly dependent on the parameter τ.

Acknowledgements. The authors are partially supported by the National Natural Science Foundation of China (Grant No. 10771178), NSAF (Grant No. 10676031), the National Basic Research Program of China (973 Program) (Grant No. 2005CB321702). Especially, the first author is supported by Hunan Provincial Innovation Foundation For Postgraduate (Grant No. CX2009B121).

Received: 27/II/09.

Accepted: 30/VII/09.

#CAM-68/09.

• [1] D.N. Arnold, R.S. Falk and R Winther, Preconditioning in H(div) and applications, Math. Comp., 66 (1997), 957-984.
• [2] A. Bossavit, Computational Electromagnetism. Variational Formulation, Complementarity, Edge Elments San Diego: Academic Press (1998).
• [3] Z. Cai, R. Lazarov, T. Manteuffel and S. McCormick, First-order system least squares for second-order partial differential equations: Part I, SIAM J. Numer. Anal., 31 (1994), 1785-1799.
• [4] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numer., 11 (2002), 237-339.
• [5] R. Hiptmair and J. Xu, Nodal auxiliary space preconditioning in H(curl) and H(div) spaces, SIAM J. Numer. Anal., 45 (2007), 2483-2509.
• [6] P. Lin, A sequential regularization method for time-dependent incompressible Navier-Stokes equations, SIAM J. Numer. Anal., 34 (1997), 1051-1071.
• [7] P. Monk, Finite Element Methods for Maxwell Equations Oxford University Press, Oxford (2003).
• [8] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp., 54 (1990), 483-493.
• [9] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Review, 34 (1992), 581-613.
• [10] L. Zhong, S. Shu, D. Sun and L. Tan, Preconditioners for higher order edge finite element discretizations of Maxwell's equations, Sci. China Ser. A, 50 (2008), 1537-1548.

Publication Dates

• Publication in this collection
  19 Mar 2010
• Date of issue
  2010

History

• Received
  27 Feb 2009
• Accepted
  30 July 2009