Journal of Computational and Applied Mathematics 205 (2007) 453 – 457 www.elsevier.com/locate/cam
Letter to the Editor
Eigenvalue estimates for a preconditioned Galerkin matrix arising from mixed finite element discretizations of viscous incompressible flows Paul Deuring Laboratoire de Mathématiques Pures et Appliquées, Université du Littoral “Côte d’Opale”, la Mi-Voix, B.P. 699, 62228 Calais cedex, France Received 14 March 2006
Abstract The article deals with Galerkin matrices arising with finite element discretizations of the Navier–Stokes system. Usually these matrices are indefinite and nonsymmetric. They have to be preconditioned if a related linear system is to be solved efficiently by an iterative method. We consider preconditioning by a pressure mass matrix. It is shown how upper and lower bounds of the eigenvalues of a preconditioned Galerkin matrix may be found by variational arguments. © 2006 Elsevier B.V. All rights reserved. MSC: 65F10; 65F35; 65N22; 65N25 Keywords: Nonsymmetric systems; Indefinite systems; Preconditioning; Finite element methods; Navier–Stokes equations
1. Introduction, description of the problem and main results This article deals with systems of linear equations which arise when the incompressible Navier–Stokes system is discretized in space by mixed finite element methods. In the case of the evolutionary Navier–Stokes system, space discretization is supposed to be preceded by a fully- or semi-implicit time discretization. In the case of the fully-implicit approach, or if the stationary Navier–Stokes system is considered, we assume that the discrete problem is linearized by a Picard iteration. Such approximations lead to discrete systems of linear equations which may be written as follows: A BT U F · = , (1) B −C P G with the blocks A, B and B T in the system matrix of (1) (“Galerkin matrix”) corresponding to a discrete convection–diffusion, gradient, and divergence operator, respectively. In the evolutionary case, the block A additionally contains a component which is due to time discretization. Usually, A is positive definite but not symmetric. Typically the matrix B has full rank in the case of nonenclosed flow, and is rank deficient by 1 else. When LBB-stable mixed finite elements are used, the matrix C is zero. Otherwise, in the case of stabilized mixed finite element methods, E-mail address:
[email protected]. 0377-0427/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2006.05.011
454
P. Deuring / Journal of Computational and Applied Mathematics 205 (2007) 453 – 457
C frequently is positive semi-definite and symmetric. The Galerkin matrix as a whole is indefinite [4, p. 285] and nonsymmetric. If the size of the matrices involved is large, Eq. (1) is often solved by iterative methods. In that case, it is imperative to precondition the Galerkin matrix in (1). As explained in [4, Section 8.1], a good choice of a right preconditioner is a BT A block triangular matrix of the form 0 − S , where A is an approximation of A, and S an approximation of the Schur A 0 complement matrix C + B · A−1 · B T . Here we will consider analogous left preconditioners of the form B − S , but our arguments remain valid in the case of right-oriented preconditioning. = A, that is, we assume that a fast solver for a discrete convection–diffusion problem is available; We suppose that A see [4, Sections 4 and 8.5] in this respect. Thus, we only consider the effect of the preconditioner S of the Schur complement. We take S = Q, where Q denotes the pressure mass matrix. This choice of S is well adapted to the case of Reynolds numbers of order 10, as indicated in [4, p. 346] and [6]. The preceding assumptions mean that we study the following preconditioned version of (1):
A B
0 −Q
−1 A · B
BT −C
U A · = P B
0 −Q
−1 F · . G
(2)
When an iterative method like GMRES is applied to a linear system, the convergence behaviour of the method depends on whether the eigenvalues of the system matrix are clustered or dispersed; compare [4, Sections 4.2, 8.3]. Therefore, in order to evaluate the efficiency of an iterative method applied to (2), it is of interest to look for upper and lower bounds of the eigenvalues of the system matrix of (2). Usually such bounds are found by matrix algebraic methods. In this respect, we mention Refs. [1–6], which in addition report on results of test computations. Many more references related to the preconditioning of Galerkin matrices as in (1) may be found in [4]. In the work at hand, we present an alternative approach for finding bounds of eigenvalues of the preconditioned Galerkin matrix in (2). This approach is based on variational arguments and has two main advantages: it is simple, and it leads to bounds which are completely explicit with respect to parameters characterizing the given flow field and the finite element method involved. Let us describe this approach and state our results. Our starting point is a variational problem arising with the discretization and linearization process mentioned at the beginning of this article. Written in an abstract form, this problem reads as follows: Find v ∈ V , ∈ M such that (v, w) + n( w , v, w) + b(w, ) = F(w) for w ∈ V , b(v, ) + c(, ) = G() for ∈ M,
(3)
where V and M are Hilbert spaces of finite dimension, , b and c are bilinear forms, n is a trilinear form, F a bounded functional on V , G a bounded functional on M, and w a fixed element of V . We suppose there are constants K1 , . . . , K6 such that |(w, z)| K1 · wV · zV , (w, w) K2 · w2V , |b(w, )|K3 · wV · M , |c(, )|K4 · M · M , |n( w , w, z)| K5 · w 0 · wV · zV , M K6 · (sup{b(y, ) : y ∈ V , yV = 1} + c(, )1/2 )
(4) (5) (6) (7)
. It is further assumed that the form c is symmetric, for w, z ∈ V , , ∈ M, where w 0 may be any norm of w c(, )0 for ∈ M, and n( w , w, z) = −n( w , z, w) for z, w ∈ V , in particular n( w, z, z) = 0(z ∈ V ). We remark that w corresponds to the velocity approximation from the preceding nonlinear iteration step. Inequality (7) constitutes a generalized LBB-condition, which reduces to the usual LBB-condition if c = 0. As an example for a finite element method which verifies (7) with c = 0, we mention the stabilized P1–P1 finite element method proposed by Rebollo [7], which is based on a variant of internal condensation of bubble functions; see [1, Sections 2 and 4.1] for a proof of (7) in that case. Concerning the LBB-case (c = 0), we refer to [4, (7.34), (7.35)]. The trilinear form c in that reference must be supplemented by the term ( 21 ) · div uh · ( uh+1 · v) dx (notations as in [4]), so as to obtain a skew-symmetric form n as in (3). Adding such an additional term in order to obtain skew-symmetry is a common trick in finite element approximations of Navier–Stokes equations.
P. Deuring / Journal of Computational and Applied Mathematics 205 (2007) 453 – 457
455
If the Navier–Stokes system is considered in normalized form, and is thus characterized by the Reynolds number, this number enters in the constant K5 in (6). On the other hand, if problem (1) arises from the Navier–Stokes system written in terms of physical velocity and pressure, the viscosity coefficient enters into K1 . Typically, it even coincides with K1 when the stationary Navier–Stokes system is discretized by mixed finite element methods. In order to establish a link between (1) and (3), let us choose a base (1 , . . . , n ) of V , and a base (1 , . . . , m ) of M, where n := dim V , m := dim M. Then define A ∈ Rn×n , B ∈ Rm×n , C, Q ∈ Rm×m by w , j , i ))1 i,j n , A := ((j , i ) + n( C := (c(i , j ))1 i,j m ,
B := (b(j , i ))1 i m,1 j n ,
Q := ((i , j )M )1 i,j m ,
where ( , )M denotes the scalar product of the Hilbert space M. Note that our definition of the matrix Q means this matrix indeed corresponds to the pressure mass matrix in a finite element context. With the preceding choice of A, B and C, variational problem (3) is equivalent to (1), and thus to (2). Abbreviate S := C + B · A−1 · B T . Since for every invertible matrix S ∈ Rm×m , in particular for S = Q, we have −1 A 0 A BT · B − S B −C A−1 I A−1 · B T A BT 0 = −1 = , · S −1 B C 0 S −1 · S S · B · A−1 − we see that unity is an eigenvalue of geometric multiplicity n of the system matrix of (2). Moreover it follows there are another m eigenvalues which may equivalently be characterized as generalized eigenvalues of the Schur complement operator in the sense that each such generalized eigenvalue verifies the equation (C + B · A−1 · B T ) · P = · Q · P
(8)
for some P ∈ Cm \{0}. Our main result concerns this Eq. (8). In fact, we will show Theorem 1.1. Suppose that B has full rank. Let ∈ C be such that Eq. (8) is valid for some P ∈ Cm \{0}. Then R 2 · (K3 + K4 ) · (1 + K3 /K2 ), R K2 /(2 · K62 · [4 · (K1 + K5 · w 0 )2 + K2 ]),
|I | 2 · K32 /K2 ,
where the constants K1 , . . . , K6 where introduced in (4)–(7). We also consider the case that the operator B is rank deficient by 1 (see Section 3), but this deficiency will not be expressed in terms of the elements of B, but via certain properties of the bilinear forms b and c. These properties correspond to the situation arising with usual mixed finite element discretizations of enclosed flows; see [4, Section 5.3], for example. In this case the bounds stated in Theorem 1.1 remain valid for nonzero eigenvalues of (8): Theorem 1.2. Suppose there is an element m0 ∈ M such that c(m0 , ) = 0 for ∈ M and b(w, m0 ) = 0 for w ∈ V . := {m ∈ M : (m, m0 )M = 0}. Suppose that inequality (7) is valid only for ∈ M. Let ∈ C be such that = 0 Put M m and Eq. (8) is valid for some P ∈ C \{0}. Then the inequalities in Theorem 1.1 remain valid. We remark that in standard mixed finite element approximations of enclosed flows, the element m0 corresponds to a constant function. Of course, variational problem (3) has to be modified in order to correspond to such approximations: instead of M, and the second equation is to be valid only the component of the solution (v, ) of (3) is to belong to M for ∈ M. 2. Proof of Theorem 1.1 Let ∈ C, P ∈ C\{0} be such that Eq. (8) is valid. Put 1 := c( 1 , ) + b(u1 , ) = R · ( 1 , )M − I · ( 2 , )M c( 2 , ) + b(u2 , ) = R · ( 2 , )M + I · ( 1 , )M
m
j =1 RPj
for ∈ M, for ∈ M,
· j and 2 :=
m
j =1 IPj
· j . Then (9) (10)
456
P. Deuring / Journal of Computational and Applied Mathematics 205 (2007) 453 – 457
where ui , for i ∈ {1, 2}, verifies the equation w , ui , w) = b(w, i ) (ui , w) + n(
for w ∈ V .
(11)
Without loss of generality, we may suppose that 1 2M + 2 2M = 1. Adding (9) with = 1 and (10) with = 2 , we obtain by (5) R = R · ( 1 2M + 2 2M ) =
2
i=1
2
(b(ui , i ) + c( i , i ))
i=1
(K3 · ui V · i M + K4 · i 2M ) (K3
+ K4 ) · 2 +
2
ui V
.
(12)
1=1
We further find by using (4), (11), (5) and our assumption on n: w , ui , ui )) = K2−1 · b(ui , i ) ui 2V K2−1 · (ui , ui ) = K2−1 · ((ui , ui ) + n( (K3 /K2 ) · ui V · i K (K3 /K2 ) · ui V , so that ui V K3 /K2
for i ∈ {1, 2}.
(13)
Thus, we may conclude from (12) that the first inequality stated in Theorem 1.1 is valid. When we take = 2 in (9) and = 1 in (10), we get by a subtraction and by (5) |I | = |I · (( 1 , 1 )M + ( 2 , 2 )M )| = | − b(u1 , 2 ) + b(u2 , 1 )| K3 · (u1 V · 2 M + u2 V · 1 M ) K3 · (u1 V + u2 V ) 2 · K32 /K2 , where the last inequality follows from (13). Thus, we have shown the last estimate in Theorem 1.1. Turning to the second one, we choose vi ∈ V such that vi V = 1
and b( vi , i ) (1/2) · sup{b(v, i ) : v ∈ V , vV = 1} for i ∈ {1, 2}.
Then we obtain with (7): 1=
2
i=1
i 2M K62 ·
2 · K62 ·
2
2
(2 · b( vi , i ) + c( i , i )1/2 )2
i=1
(4 · b( vi , i )2 + c( i , i )).
(14)
i=1
On the other hand, we have for i ∈ {1, 2}: |b( vi , i )| = |(ui , vi ) + n( w , ui , vi )|(K1 + K5 · w 0 ) · ui V ,
(15)
where we used (11), (4), (6) and the assumption vi V = 1. But by (4) and (9)–(11): 2
i=1
ui 2V K2−1 · ((ui , ui ) + n( w , ui , ui )) = K2−1 · K2−1 ·
2
i=1
2
b(ui , i )
i=1
(b(ui , i ) + c( i , i )) = K2−1 · R ·
2
i=1
( i , i )M = K2−1 · R .
(16)
P. Deuring / Journal of Computational and Applied Mathematics 205 (2007) 453 – 457
457
This means in particular that 2
b(ui , i )0.
(17)
i=1
Combining (15) and (16) yields 2
b( vi , i )2 (K1 + K5 · w 0 )2 · K2−1 · R .
(18)
i=1
We further find with (17), (9) and (10): 2
c( i , i )
2
(c( i , i ) + b(ui , i )) = R ·
i=1
i=1
2
( i , i )M = R .
(19)
i=1
Thus, by (18) and (19): 2
(4 · b( vi , i )2 + c( i , i ))[4 · (K1 + K5 · w 0 )2 · K2−1 + 1] · R .
i=1
This estimate and inequality (14) imply the second estimate stated in Theorem 1.1.
3. The case of enclosed flows: proof of Theorem 1.2 as in Theorem 1.2, and we require the same assumptions as in that theorem. Then the We choose m0 and M coefficient matrix in (1) is rank deficient by 1, as follows by the arguments below. In order to prove Theorem 1.2, we take 1 , . . . , m ∈ R with m0 = m · i=1 i i . Then C · = (c(i , m0 ))1 i m = 0
and B T · = (b(j , m0 ))1 j n = 0,
so that (C + B · (A−1 )T · B T ) · = 0. Now let ∈ C\{0} and P ∈ Cm \{0} be such that Eq. (8) is satisfied. Define 1 and 2 as at the beginning of Section 2. Then (m0 , 1 )M + i · (m0 , 2 )M = T · Q · RP + i · T · Q · IP = T · Q · P = −1 · T · · Q · P = −1 · T · (C + B · A−1 · B T ) · P = −1 · P T · (C + B · (A−1 )T · B T ) · = 0,
(20)
where we used the relation (C + B · (A−1 )T · B T ) · = 0 in the last equation. It follows from (20) that (m0 , i ) = 0 This means that the reasoning in Section 2 remains valid under the assumptions of for i ∈ {1, 2}, hence 1 , 2 ∈ M. Theorem 1.2. Therefore R and I are bounded as stated in Theorem 1.1. References [1] C. Calgaro, P. Deuring, D. Jennequin, A preconditioner for generalized saddle point problems: application to 3D stationary Navier–Stokes equations, Numer. Methods Partial Differential Equations, to appear. [2] Z.-H. Cao, Fast iterative solution of stabilized Navier–Stokes systems, Appl. Math. Comput. 157 (2004) 219–241. [3] H. Elman, D. Silvester, Fast nonsymmetric iterations and preconditioning for Navier–Stokes equations, SIAM J. Sci. Comput. 17 (1996) 33–46. [4] H. Elman, D. Silvester, A. Wathen, Finite Elements and Fast Iterative Solvers, with Applications in Incompressible Fluid Dynamics, Oxford University Press, Oxford, 2005. [5] A. Klawonn, G. Starke, Block triangular preconditioners for saddle point problems with a penalty term, SIAM J. Sci. Comput. 19 (1998) 172–184. [6] P. Krzyzanowski, On block preconditioners for nonsymmetric saddle point problems, SIAM J. Sci. Comput. 23 (2001) 157–169. [7] T.C. Rebollo, A term by term stabilization algorithm for finite element solution of incompressible flow problems, Numer. Math. 79 (1998) 283–319.