On a MIC(0) preconditioning of non-conforming mixed FEM elliptic problems

Mathematics and Computers in Simulation 76 (2007) 149–154

On a MIC(0) preconditioning of non-conforming mixed FEM elliptic problems S. Margenov a,∗ , P. Minev b b

a Institute for Parallel Processing, Bulgarian Academy of Sciences, Sofia, Bulgaria Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, Canada T6G 2G1

Available online 23 January 2007

Abstract In this paper we analyze a preconditioner for mixed finite element systems arising in the approximation of a second order elliptic problem with Neumann boundary conditions by triangular non-conforming elements. This problem stems from the so called projection methods for the unsteady Navier–Stokes equations and is one of the most computationally intensive parts of the method. The present study is focused on the efficient implementation of the modified incomplete LU factorization MIC(0) as a preconditioner in the PCG iterative method for the reduced Schur complement system following from the mixed formulation of the pressure Poisson problem. The presented model analysis and the related set of numerical tests illustrate well the potential of the proposed approach. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. AMS Classification: 65F10; 65N30 Keywords: Non-conforming mixed FEM; Preconditioning

1. Introduction The projection methods are probably the most popular methods for solving the unsteady Navier–Stokes equations. Their popularity is due to the fact that they reduce the solution of the generalized Stokes problem to the solution of a vectorial advection-diffusion equation for the velocity and a second order elliptic problem for the pressure. In a recent paper (see [4]) it was proposed to use P1 elements for the discretization of the advection-diffusion equation that follows from the projection approach, and first order nonconforming Crouzeix–Raviart elements for the discretization of the projection step. As a result, the projection step takes the following mixed form: ¯ u − ∇p = f in Ω,

∇ ·u=0

in Ω,

u·n=0

on ∂Ω,

(1)

where Ω is a convex polygonal domain in R2 and f is a given smooth vector valued function. In the paper we use the CFD terminology and refer to the vector u as a velocity and p as a pressure. The weak formulation of the problem (1) reads: 2 For a given f ∈ (L2 (Ω)) find (u, p) ∈ (V, Q), satisfying (u, v) + (p, ∇ · v) = (f, v) ∀v ∈ V,

(∇ · u, q) = 0

∀q ∈ Q

2

where V = {v ∈ (L2 (Ω)) ; ∇ · v ∈ L2 (Ω); v · n = 0 on ∂Ω} and Q = {q ∈ L2 (Ω); ∗

Corresponding author. E-mail addresses: [email protected] (S. Margenov), [email protected] (P. Minev).

0378-4754/$32.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.01.021

(2) 

Ω q dΩ

= 0}.

150

S. Margenov, P. Minev / Mathematics and Computers in Simulation 76 (2007) 149–154

Elliptic problems can be discretized by a variety of methods and each of them has its advantages and disadvantages. It is usually the particular application that favours one choice versus another. For CFD applications, in particular, the use of nonconforming (discontinuous pressure) elements for the projection yields a locally divergence-free velocity which is a desirable feature in many occasions. The discretization of (2) proceeds as follows. Suppose that Ω is partitioned using triangular elements. The partition is denoted by Th and is assumed to be quasi-uniform with a characteristic mesh-size h (for a precise definition of a quasi-uniform mesh see [10], p. 106). Then the discretization of V is given by Vh = 2 {vh ∈ (L2 (Ω)) ; vh ∈ (P1 (e))2 , ∀e ∈ T; vh is continuous at mi,e ∈ Ω; vh (mj,e ) · n = 0, ∀mj,e ∈ ∂Ω}, where mj,e is the midpoint of the i-th side of the element e (i = 1, 2, 3). The pressure is approximated by piecewise constants on T  and the corresponding discrete space is Qh = {qh ∈ L2 (Ω); qh |e = P0 (e), ∀e ∈ T; Ω qh = 0}. The spaces (Vh , Qh ) satisfy the inf-sup stability condition associated with saddle point problems. In general, this is not enough to guarantee the convergence of the projection algorithm as demonstrated in [15]. However, it was shown in [4] that if the Crouzeix–Raviart elements are used to discretize the projection step (1) then the corresponding approximation for the velocity converges to the exact solution of the original generalized Stokes problem. The finite element discretization of the problem leads to the system Pwh = b,

(3)

where wTh = {uTh , ph } with uTh , ph being the vectors of the velocity and pressure unknowns. The matrix P results from a symmetric saddle point problem and therefore is symmetric and indefinite, i.e., ⎤ ⎡ M B1 ⎥ ⎢ M B2 ⎦ . (4) P =⎣ B1T

B2T

Here, M is the mass matrix corresponding to the scalar FEM space used to approximate the velocity components. Note, that M is a diagonal matrix. This follows directly from the fact, that the quadrature formula on a triangle with nodes in the midpoints of the edges is exact for second degree polynomials. The solution of linear systems resulting from saddle point problems has been discussed by many authors (see [3] and the references therein). For this system we can either eliminate the pressure unknowns and end up with a system for uh that corresponds to a divergence-free basis in Vh (see [9]) or, since the mass matrix is diagonal, we can eliminate the velocity unknowns and derive a system for the pressure with a symmetric and positive semidefinite matrix S = B1T M −1 B1 + B2T M −1 B2 .

(5)

In this paper we discuss the second approach and derive a preconditioner for the Schur complement (5) based on the modified Cholesky factorization MIC(0). The rest of the paper is organized as follows. In Sections 2 we present the MIC(0) algorithm. A model analysis of the proposed PCG algorithm for the considered mixed FEM problem is given in Section 3. Finally, in Section 4 we present numerical results on a test problem, that illustrate the computational efficiency of the considered iterative solution method. 2. MIC(0) factorization Let us first recall some known facts about the modified incomplete factorization MIC(0) preconditioner applied to a matrix A. Our presentation at this point follows the setting of [6]; see also [12,13] for some more details. Let us rewrite the real N × N matrix A = (aij ) in the form A = D − L − LT ,

(6)

where D is the diagonal and (−L) is the strictly lower triangular part of A. Then we consider the approximate factorization of A which has the following form: CMIC(0) (A) = (X − L)X−1 (X − L)T ,

S. Margenov, P. Minev / Mathematics and Computers in Simulation 76 (2007) 149–154

151

where X = diag(x1 , . . . , xN ) is a diagonal matrix determined by the condition of equal rowsums: CMIC(0) (A)e = Ae,

e = (1, . . . , 1)T ∈ RN .

For the purpose of preconditioning, we are interested in the case when X > 0 and thus CMIC(0) (A) is positive definite. If this holds, we call the MIC(0) factorization stable. We have the following theorem about the stability of the MIC(0) factorization. Theorem 1. Let A = (aij ) be a symmetric real N × N m atrix and let A = D − L − Lt be the splitting (6) of A. Let us assume that L ≥ 0,

Ae ≥ 0,

Ae + Lt e > 0,

e = (1, . . . , 1)t ∈ RN ,

i.e. that A is a weakly diagonally dominant matrix with non-positive off-diagonal entries and that A + Lt = D − L is strictly diagonally dominant. Then the relation xi = aii −

i−1 N aik k=1

xk

akj

j=k+1

gives positive values xi and the diagonal matrix X = diag(x1 , . . . , xN ) defines a stable MIC(0) factorization of A. For a complete proof the reader is referred to [12]. Remark 2. The numerical tests presented in the last section are performed using the perturbed version of the ˜ = A + D. ˜ The diagonal perturbation MIC(0)algorithm, where the incomplete factorization is applied to the matrix A ˜ ˜ ˜ ˜ D = D(ξ) = diag(d1 , . . . dN ) is defined as follows:

ξaii if aii ≥ 2wi d˜ i = , ξ 1/2 aii if aii < 2wi  where 0 < ξ < 1 is a constant and wi = j>i − aij . Let us focus on the preconditioned conjugate gradient (PCG) solution of the reduced system for the Schur complement S, using the preconditioner C = CMIC(0) (S). Each PCG iteration consists of one solution of problem with the preconditioning matrix C, one matrix vector multiplication with the original matrix S, two inner products, and three linked vector triads of the form v:=αv + u. It is easy to show, that the matrix S remains sparse in the considered setting. Moreover, there are at most four nonzero entries in each of its rows. Therefore, the computational complexity of one PCG iteration is given by NitPCG ≈ N(C−1 v) + N(Sv) + 10N. Then, N(C−1 v) ≈ 9N,

N(Sv) ≈ 7N,

and finally CCNitPCG ≈ 26N.

(7)

Obviously, the algorithm under consideration is relatively cheap, since the solution of the preconditioning system takes only about one third of the total cost. 3. Model analysis This model analysis concerns mostly the structure of the Schur complement matrix S and its properties and tries to determine if it is possible to apply MIC(0) factorization for a PCG iterative solution of the reduced problem for the pressure Ph . Before we start, we should note that it is quite difficult to impose the condition Ω qh dΩ = 0 on the space for the pressure Qh . Therefore, the usual approach is to relax this requirement on Qh , impose a homogeneous

152

S. Margenov, P. Minev / Mathematics and Computers in Simulation 76 (2007) 149–154

Fig. 1. Schur complement four point stencil for the pressure.

Dirichlet condition to one, let say the first, of the components of the vector ph , and then impose the integral condition above on the so computed solution. This procedure yields a symmetric and positive definite Schur complement S. Let us consider the problem in the unit square Ω = (0, 1)2 . Let us also assume that a uniform mesh is used with a mesh size parameter h = 1/n. Then, the mass matrix is easy to compute, namely M=

h2 I. 3

The Schur complement is a discrete elliptic operator acting on Qh . As it was already mentioned, due to the diagonal structure of M, the matrix S remains sparse. In fact, S could be interpreted as a stiffness matrix corresponding to a T-shaped four point stencil. The particular coefficients for the model problem are shown in Fig. 1. The conclusion from Fig. 1 is that the Schur complement is an M-matrix. Obviously, the first two conditions for a stable MIC(0) factorization (see Theorem 1) are satisfied. They are independent of the ordering of the unknowns, corresponding to S. This is not the case with the last condition of the theorem. Here we consider two possible element numberings presented in Fig. 2. Both are natural from the FEM point of view. The first one follows the vertical columns of elements, and is denoted by TW (triangle-wise). The second variant follows the vertical columns of square cells, and is therefore denoted by SW (square-wise). In both cases, the Schur complement S has a block-tridiagonal form, i.e., S = blocktridiag (Si,i−1 , Si,i , Si,i+1 ). The diagonal blocks Si,i are diagonal matrices for variant TW which could be very attractive from algorithmic point of view. This could in particular be very useful for the parallel implementation of the algorithm. Unfortunately, D − L is not strictly diagonally dominant in this case, i.e., the last condition from Theorem 1 is not fulfilled. The structure of S does not look so nicely if elements’ numbering SW is applied, but the conditions for a stable MIC(0) factorization are fully satisfied in this case. We have used the variant SW to implement the resulting algorithm.

Fig. 2. Element numbering of the pressure unknowns (a) variant TW; (b) variant SW.

S. Margenov, P. Minev / Mathematics and Computers in Simulation 76 (2007) 149–154

153

Table 1 CG and PCG iterations: MIC(0) factorization of S n

N

CG

PCG

16 32 64 128 256 512

512 2048 8192 32768 131072 524288

86 177 358 737 1487 3010

13 19 28 41 59 87

4. Numerical tests The numerical tests presented below illustrate the PCG convergence rate when the size of the discrete problem varies. A relative stopping criterion (C−1 rnit , rnit ) <ε (C−1 r0 , r0 ) is used in the PCG algorithm, where r i stands for the residual at the ith iteration step, (·, ·) is the standard Euclidean inner product, and ε = 10−6 . The computational domain is the unit square Ω = (0, 1)2 . A uniform mesh is used, where h = 1/n, and the size of the discrete problem is N = 2n2 . The qualitative analysis of the results given in Table 1 clearly confirms that: (a) the number of iterations for CG 1/2 1/4 method is O(n) = O(N 1/2 √ ); (b) the number of iterations for MIC(0) PCG method is O(n ) = O(N ), namely, it grows proportionally to n. 5. Concluding remarks In this paper we propose a preconditioning technique for the iterative solution of mixed finite element systems arising in approximations of second order elliptic boundary value problems by non-conforming finite elements. Our study was motivated by more general applications in CFD, where the Crouzeix–Raviart non-conforming elements are shown to be very a attractive approximation tool because it yields a velocity field which is locally divergence-free. The considered approach is based on the elimination of the velocity unknowns in the system that results from the mixed FE discretization of the problem. It was shown, that the Schur complement matrix for the reduced pressure problem preserves a strong sparsity satisfying the conditions for a stable MIC(0) factorization. The presented numerical tests well illustrate the potential of the reported method, algorithm and code, for the efficient solution of large scale problems. Our further plans include a generalization to 3D problems on tetrahedral meshes. We plan also to consider a twolevel modification of our algorithm based on a first step of aggregation/differentiation of the Schur complement (see, e.g., [7,14]). Depending on the aggregation strategy, we could get for the second (elliptic) diagonal block of the obtained hierarchical presentation of S: (a) a more regular structure, like the standard five point (seven point in 3D case) stiffness matrix, or (b) a better for parallel implementation structure, which is similar to what we have in case of Rannacher–Turek rotated quadrilateral elements [5,16]. Acknowledgments This work was supported in part by BiS21++ EC Grant, Bulgarian NSF Grant I-1402/2004 and a Discovery Grant of NSERC. References [3] O. Axelsson, M. Neytcheva, Preconditioned methods for linear systems arising in constrained optimization problems, Numer. Linear Algebra Appl. 10 (2003) 3–31.

154

S. Margenov, P. Minev / Mathematics and Computers in Simulation 76 (2007) 149–154

[4] B. Bejanov, J.-L. Guermond, P. Minev, A locally div-free projection scheme for incompressible flows based on non-conforming finite elements, Int. J. Numer. Meth. Fluids 49 (2005) 239–258. [5] G. Bencheva, S. Margenov, Parallel incomplete factorization preconditioning of rotated linear FEM systems, J. Comput. Appl. Mech. 4 (2) (2003) 105–117. [6] R. Blaheta, Displacement decomposition-incomplete factorization preconiditioning techniques for linear elasticity problems, Numer. Linear Algebra Appl. 1 (1994) 107–126. [7] R. Blaheta, S. Margenov, M. Neytcheva, Uniform estimate of the constant in the strengthened CBS inequality for anisotropic non-conforming FEM systems, Numer. Linear Algebra Appl. 11 (2004) 309–326. [9] S. Brenner, A nonconforming multigrid method for the stationary Stokes equations, Math. Comp. 55 (1990) 411–437. [10] S. Brenner, L. Scott, The mathematical theory of finite element methods, Springer, 1994. [12] I. Gustafsson, Modified incomplete Cholesky (MIC) factorization, in: D.J. Evans (Ed.), Preconditioning Methods; Theory and Applications, Gordon and Breach, 1984, pp. 265–293. [13] I. Gustafsson, An incomplete factorization preconditioning method based on modification of element matrices, BIT 36 (11) (1996) 86–100. [14] R. Lazarov, S. Margenov, On two-level MIC(0) preconditioning of Crouzeix–Raviart nonconforming FEM systems, Springer Lecture Notes Comp. Sci. 2542 (2003) 192–201. [15] K.A. Mardal, X.-C. Tai, R. Winther, A robust finite element method for Darcy–Stokes flow, SIAM J. Numer. Anal. 40 (5) (2002) 1605–1631. [16] R. Rannacher, S. Turek, Simple nonconforming quadrilateral Stokes Element, Numer. Meth. PDEs 8 (2) (1992) 97–112.