Analysis of some projection method based preconditioners for models of incompressible flow
Abstract
In this paper, several projection method based preconditioners for various incompressible flow models are studied. In particular, we are interested in the theoretical analysis of a pressurecorrection projection method based preconditioner [14]. For both the steady and unsteady Stokes problems, we will show that the preconditioned systems are well conditioned. Moreover, when the flow model degenerates to the mixed form of an elliptic operator, the preconditioned system is an identity no matter what type of boundary conditions are imposed; when the flow model degenerates to the steady Stokes problem, the multiplicities of the nonunitary eigenvalues of the preconditioned system are derived. These results demonstrate the effects of boundary treatments and are related to the stability of the staggered grid discretization. To further investigate the effectiveness of these projection method based preconditioners, numerical experiments are given to compare their performances. Generalizations of these preconditioners to other saddle point problems will also be discussed.
rojection method, preconditioning, Stokes equations, saddle point problems
1 Introduction
We consider to solve the incompressible Stokes equations
(1) 
subject to suitable boundary conditions on . Here, is the velocity, is the pressure, is the density function, is the viscosity, may include an external force and the nonlinear term , is 0 or a source term. For simplicity, we assume that , and are constants.
Different methods have been proposed and studied for the solution of the above problem. The original projection method, introduced by Chorin and Temam [7, 39], decouples the computation of velocity and pressure into solving several Poisson problems. It yields time stepping schemes for the incompressible Stokes equations which require relatively simple linear solvers. However, this approach complicates the specification of physical boundary conditions and may cause artificial numerical boundary layers [9]. Appropriate ”artificial” boundary conditions must be prescribed for the intermediate variables and , and obtaining highorder accuracy for and may not be possible in some cases [9, 15]. Compared with the original projection method, implicit schemes which couple the computations of and together can provide accurate numerical solution and possess better stability [15, 29]. However, implicit coupled schemes require efficient solvers for the resulted saddle point systems, which are usually difficult to solve. Using preconditioned Krylov methods provides a good remedy to overcome the difficulties [2, 20, 28]. Consequently, designing effective preconditioner is of vital importance. In a recent work [14], a pressurecorrection projection method [33, 14, 34] was adopted as the preconditioner for the resulting saddle point system. Numerical experiments [14, 6] show that the preconditioner is effective and efficient. However, we note that the preconditioner in [14] is restricted to the time dependent Stokes case and the corresponding theoretical justification is absent. In this work, we are interested in the following aspects: generalizing the preconditioner in [14] to steady Stokes problem; providing a theoretical analysis for both the steady and unsteady problems; studying some related preconditioners which are based on exact or spectralequivalent approximations of the Schur complement; building up the connections between our work with some existing works; extending the projection method based preconditioners to other saddle point problems.
In this paper, we use an implicit time scheme and a staggered grid spatial discretization for (1). The resulting saddle point system y using a preconditioned GMRES method [37, 36]. One of our interests is analyzing the pressurecorrection projection method based preconditioner [14]. The obtained results are as follows. For both the steady and unsteady Stokes problems, it is shown that the preconditioned systems are well conditioned. More precisely, the nontrivial eigenvalues of the preconditioned systems have uniform lower and upper bounds which are independent of the mesh refinement and the physical parameters. Moreover, for steady Stokes equations, the multiplicities of the nonunitary eigenvalues of the preconditioned system are derived based on the analysis of the commutator difference operator. Specifically, for two dimensional Dirichlet boundary value problem, if there are cells along each direction, it is shown that there are at most eigenvalues not equal to 1. For periodic boundary value problem, the preconditioned operator is an identity operator. Meanwhile, motivated by other related works [2, 3, 11, 20, 25, 26, 27, 28, 31, 32], we propose and investigate several projection method based preconditioners. These preconditioners use exact or approximate inverse of the Schur complement. Their effectiveness, advantages and disadvantages, and the connections with those existing works will be highlighted. In addition, we will discuss how to generalize our work to other saddle point problems.
This paper is organized as follows. In section 2, we present the time and the spatial discretization of the Stokes problem, the projection method based preconditioners, the corresponding matrix representations, and other related preconditioners. In section 3, eigenvalue analysis of the preconditioned systems are presented. In section 4, numerical experiments are given to compare the performances of different preconditioners. Concluding remarks are drawn in the last section.
2 Discretizations and precondtioners
2.1 Time and spatial discretizations
To have better stability and accuracy properties, an implicit time stepping scheme is adopted. In this paper, we apply the following backward Euler scheme.
(2) 
Here, and are the approximate solution of and at .
The staggered grid (i.e. markerandcell or MAC [17]) method is applied for the spatial discretization. In this work, we assume that the computational domain is a two dimensional rectangle and there are and cells along the direction and the direction respectively. For simplicity, let us further assume that and (and therefore ). To present the staggered grid discretization, we introduce the following sets.
As shown in Figure 2.1, , and are the point sets for , and respectively. The divergence of is approximated at cell centers by with
is well defined at all points in . The gradient of is approximated at the and edges of the grid cells by with
We see that is well defined in the interior of and is well defined in the interior of . There are three different approximations to the Laplace operator required in the staggered grid scheme, one defined at the cell centers, one defined at the edges of the grid, and one defined at the edges of the grid. All employ the standard 5point finite difference stencil. The Laplacian of is approximated at cell centers by
The approximations to the Laplacians of and are
The finite difference approximation to the vector Laplacian of is denoted as . We see that is well defined in and is well defined in .
To evaluate the above finite difference approximations near the boundaries of , we need the specification of boundary values along and ”ghost values” located outside of . For periodic boundary conditions, the ghost values at boundaries are copies of the corresponding interior values. For example, the ghost values at the west boundary are copied from the corresponding function values which are evaluated at the nodes nearest to the east boundary. To demonstrate the treatment of Dirichlet boundary condition, we consider the discretization of near the south boundary. By using the standard stencils, we evaluate
(3) 
Here, is the  component of the velocity at the node . Similar treatments are applied to near the west and east boundaries. This type of boundary treatment for keeping the symmetry of the discrete system can be interpreted as follows. If is defined on (i=1,2), then we denote the function which is defined on and coincides with on as . If , we then set . If , we set with being the nearest node of to . By this way, the 0 Dirichlet boundary condition is weakly imposed and for the points in , and there holds . For Dirichlet boundary condition, the discretization based on (3) for approximating the second order derivatives has local truncation order while the overall discretization is spatially globally secondorder accurate [18, 29].
After the time discretization and the staggered grid spatial discretization, the fully discrete system is of the following saddle point form.
(4) 
Here, is the discrete operator for , for , and for , includes both the force terms and the velocity terms from the previous time step, and represents the source term and the boundary data. Moreover, no matter what type of boundary conditions, we note that the staggered grid discretization always satisfies the compatibility: . That is, the two operators are adjoint to each other. Moreover, with being the discrete Laplacian operator applied to the pressure variable. As previously introduced, is welldefined at cell centers. In addition, the boundary condition for is naturally determined by the boundary condition of the global system: when the boundary conditions for the Stokes system are Dirichlet type, is imposed with Neumann boundary conditions, when boundary conditions for the Stokes system are periodic, is also imposed with periodic boundary conditions.
2.2 Projection method based preconditioners
The pressurecorrection projection algorithm [22, 14, 33, 34] for the constant density and constant viscosity Stokes equations is as follows.

Step 1: Solve for an intermediate velocity using an implicit viscosity step:
(5) 
Step 2: Note that does not generally satisfy the discrete incompressibility constraint, we project to the divergence free space by solving
(6) This leads to a Poisson equation for the intermediate variable :
(7) 
Step 3: Update the velocity and correct the pressure term by
(8) (9)
This algorithm is analogous to the second order projection method whic is based on the CrankNicolson scheme [22, 14]. From (5)(9), the above projection algorithm can be written into a factorization form.
(10) 
Here is rank one deficient. is in the sense that . As pointed out in [6], the last 2 steps can be combined together and therefore leads to the following more efficient implementation.
(11) 
Here,
(12) 
is an approximation of the inverse of the Schur complement
(13) 
The approximation property of (12) to (13) was firstly found in [5] and has been justified using the Fourier analysis approach [5, 6, 25] and the operator mapping theory [3, 26, 27, 31, 32, 40].
Note that the pressure correction (9) is based on the following arguments. Multiplying by on both sides of (6), and noting that and satisfy the discrete Stokes equation (2) and satisfies (5), we derive that
(14) 
Multiplying by on both sides of (14), the least square solution of (14) satisfies
Therefore, we see that the pressure correction can be
(15) 
By assuming that there holds the commutation property , we obtain (9). Generally, the commutation property does not hold and only in special cases, such as the boundary conditions are periodic.
From the above discussions, if the pressure correction step is modified to be (15), then we obtain
(16) 
Furthermore, for steady Stokes equation, . By setting in (10), one obtains the projection method based preconditioner for the steady Stokes problem:
(17) 
From our study [6], the projection algorithm expressed as (17) can not be directly applied as a solver for the stationary Stokes. For instance, if (17) is used for solving the driven cavity problem, one can not resolve the corner vortices. However, we note that (17) is still a robust preconditioner for (4).
Motivated by the theory developed in [28, 20] and the above discussion, we see that one can use the exact inverse of the Schur complement. Then, we obtain the following projection type preconditioner.
(18) 
Remark 1. In this work, the derivations of the preconditioners are based on the backward Euler scheme (2). The derivation based on the CrankNicolson scheme can be found in [14]. In addition, we note that is similar to the BFBt preconditioner in [11] if the Oseen equations degenerates to the Stokes equations.
2.3 Other related preconditioners
As is a good approximation of the inverse of the Schur complement, we consider to extend the block diagonal preconditioners discussed in [5, 26, 27, 31, 32] to the following lower and upper triangular preconditioners. We study
(19) 
and
(20) 
These preconditioners are actually the approximations of the following block triangular preconditioners
(21) 
respectively.
In the following, we will study all the preconditioners mentioned above. In particular, we will discuss the advantage, the disadvantages, and the algebraic properties of all the precondtioners.
2.4 Matrix representations
In this subsection, we present the matrix representations of the discrete differential operators. If the lexicographic order [12] (indexing the unknowns from the left to the right and from the bottom to the top) of the degrees of freedom is used, from the staggered grid discretization and the boundary treatments, then the resulting linear system for the 2D Stokes problem with Dirichlet boundary on a square domain is a saddle point problem of the form (4). The corresponding matrices are [11]
Here, and are the matrix representations of the discrete and the discrete according to the above variables ordering, and
If periodic boundary conditions are imposed in both the  and  directions, then
Here and are circulant matrices [8, 38].
Moreover, we have
Here, is a Fourier matrix and is a diagonal matrix [38] with the eigenvalues as its diagonal entries. Therefore, when boundary conditions are purely periodic, there holds the commutation property: . Moreover, a Fourier matrix times a vector can be realized by applying the Fast Fourier transform to the vector. Therefore, FFT is usually used for the resulted system [38].
Based on the above discussions and explanations, if the boundary condition is periodic in the  direction and is Dirichlet in the  direction, it is not difficult to write out the corresponding matrix representations either.
3 Analysis
To link with the original saddle point form (4), we formally calculate
(22)  
(23) 
Here, we assume that is invertible. In implementation, times a pressure vector is evaluated in a way which ensures that the nullspace is handled consistently and the mean pressure value is kept constant. Note that , when the boundary condition is purely periodic, the basic operators , , and are commutative, is exactly the discrete identity operator. Most importantly, from (22), we see that the preconditioned system is block upper triangular. Therefore,
(24) 
It means that the eigenvalues of the preconditioned system are either 1 or the eigenvalues of the block of (22).
Similarly, for the other preconditioners, it is easy to derive that
(25) 
and (or calculate )
(26) 
From (22)(26), we see that , and actually share the same spectral. Moreover, the eigenvalues of the preconditioned systems are either 1 or the eigenvalues of the . Therefore, we shall analyze the spectral of under the staggered grid discretization. It should be pointed out that, in (22), the block contains a projection operator . The projection operator, maps a velocity vector to . Intuitively, (22) has better approximation properties than (25) or (26). We will check that whether does provide better performance or not.
For the exact inverse Schur complement based precondtioner, , we have the following proposition, which is an extension of Remark 2 in [28].
Let be the preconditioner defined in (18), we have
(27) 
The preconditioned system have the minimal polynomial .
3.1 Analysis of the preconditioners for the stationary Stokes problem
By setting in (10), (16) and (18), we obtain the projection method based precondtioner (17) for the stationary Stokes problem. From (24) and noting that in the stationary Stokes case, we see that the eigenvalues of the preconditioned system are either 1 or the eigenvalues of the Schur complement (13). In this subsection, we mainly discuss the preconditioned system which uses (17) as the preconditioner. In particular, we will focus on the analysis for the Dirichlet boundary problem.
Let be the linear space of vector functions defined on and be the linear space defined on satisfying mean value 0. For homogeneous Dirichlet boundary problem, the boundary treatments are stated in section 2.1. Moreover, we endow and with the following inner products and norms.
In the definition of norm, we have implicitly used the 0 Dirichlet boundary conditions. We see that the functional spaces and are approximations of and respectively [23, 24].
The staggered grid discretization is uniformly divstable [13, 16, 23, 30], it means that the following lemma holds [23].
There exists a constant independent of the mesh refinement such that
(28) 
The proof of this lemma in finite element terminology can be found in [13, 16]. The constant depends on the shape of the domain (for example, the ratio of the length along the  direction and the length along the  direction [30]) but independent of the mesh refinement. This lemma states that is a surjective from to and the staggered grid discretization is uniformly divstable [4].
To study the spectral of the Schur complement, we introduce the following lemma. It points out that the spectrum of the Schur complement is equivalent to the solution of a generalized eigenvalue problem.
If and are an nonzero eigenvalue and the corresponding eigenvector of the Schur complement, i.e.,
(29) 
then is the eigenvalue of the generalized eigenvalue problem:
(30) 
Furthermore, if is nonzero, the inverse of the above argument also holds. {proof} Multiplying both sides of (29) by and denoting , we see that
Moreover, as and are similar to each other, they share the same spectral. Therefore, it is sufficient to investigate the eigenvalue problem
This eigenvalue problem is equivalent to the generalized eigenvalue problem (30) if is nonzero. We therefore get the conclusion.
From Lemma 3.1, we see that the nonzero eigenvalues of the Schur complement is equivalent to the generalized eigenvalue of (30). Moreover, because both and are symmetric, by using (30), we have
Here, and are the maximum and the minimum eigenvalues of the Schur complement respectively. Thus, to estimate the bounds of the spectral of the Schur complement, it is sufficient to study the generalized Rayleigh quotient,
(31) 
Remark 2. When (1) degenerates to the stationary Stokes problem, one can absorb the inverse of the viscosity into the pressure term or one can scale in approximating the inverse of the Schur complement. Therefore, without loss of the generality, we will assume that and in this section.
The following theorem states the eigenvalue analysis of the preconditioned system for the steady Stokes problem.
For the steady Stokes problem with purely Dirichlet boundary condition, we have (i) . The multiplicity of 0 eigenvalue is 1. (ii). There are at most eigenvalues not equal to 1. (If , there are at most eigenvalues not equal to 1.) For boundary condition with the direction periodic and the direction Dirichlet, there are at most eigenvalues not equal to 1.
(i) The multiplicity of the 0 eigenvalue is easy to understand because the dimension of the null space of the Schur complement is 1. The lower bound of the eigenvalues of the Schur complement is , which is independent of the mesh refinement [12, 23, 30]. This result is derived from the definition of the discrete infsup condition [12, 23, 30]. To see this, we check
(32)  
(33)  
(34)  
(35)  
(36) 
Here means the positive square root of (that is ).
To estimate the upper bound of the eigenvalues, we study the properties of the Rayleigh quotient (31). Firstly, we note that
To check the eigenvalues of the Schur complement are not greater than 1, it is sufficient to prove
(37) 
By using the matrix representations in the previous section and the properties of the Kronecker product, we have
(38) 
For any block symmetric operator, there holds
Thus, to prove is symmetric semipositive definite, it is sufficient to prove (SPD) and are symmetric positive semidefinite. By using the properties of the Kronecker product,
Since is positive semidefinite, we see that is positive semidefinite. Therefore (37) holds.
(ii). To show that there are at most nonzero eigenvalues, we study the commutation difference operator [15].
From (38), we see that
(39) 
The right hand side is a matrix with its rank equals to . It means that are all zeros in the interior of the domain. Moreover, for any nontrivial eigenvector , there exists such that . Therefore, from Lemma 3.1 and (39), we see that there are at most eigenvalues are not equal to 1.
If the boundary conditions are purely periodic, then all the nonzero eigenvalues of are equal to 1. From the commutation property and Lemma 3.2, the nontrivial eigenvalues of are equivalent to the generalized eigenvalues of . As , all the nontrivial eigenvalues of are equal to 1. For the case with the boundary condition which is periodic in  direction and Dirichlet in  direction, the proof is similar to the above discussion.
Remark 3. From Theorem 3.1 and the infsup condition, we conclude that
(40) 
The ideas of the proof of Theorem 3.1 actually originate from some differential operator identities. Note that in continuous level, if the function is smooth enough, the following identities of differential operators [1] hold.
(41) 
and
(42) 
as
(43) 
The communication difference operator in (38) is actually the matrix representation of except that there are some modifications near the boundaries (because of the boundary treatment). Moreover, the operators and are the actually discrete analogies of (41)(43).
BC type  DOF  No. nonunitary eigs  

pure Diri  16  16  736  60 
pure Diri  32  32  3008  124 
xperiodic yDiri  16  32  1520  31 
xperiodic yDiri  32  64  6112  63 
To verify the conclusions in Theorem 3.1, we present some numerical experiments. In Table 1, we list the number of the nonunitary eigenvalues under different types of the boundary conditions and different numbers of the grid points. The computation is based on uniform mesh partitions (with the meshsize equals to 1) on rectangular domains. We form the Schur complement and call Matlab function to calculate the eigenvalues. For purely Dirichlet boundary condition, the total DOF is equal to . If the boundary condition along  direction is periodic, the total DOF is equal to . If the boundary conditions are purely periodic, the total DOF is equal to . The eigenvalue calculations are based on (22) and (24). From Table 3.1, we see that the results verify the theoretical predictions. We comment here that there are similar results in [25] and the reference therein [24]. The discussions therein are based on the comparisons of the dimensions of functionals spaces under different boundary condition types.
3.2 Analysis of the preconditioners for the unsteady models
For the time dependent Stokes problem, we can again get started from (22) and (24). As and share the same spectral with that of . The following spectral analysis is valid for all these three preconditioners.
To simplify our presentation, we analyze an equivalent system like that in [26]. Multipling (4) by , denoting as and absorbing the scaling into the pressure term, we see that the resulting system is again a saddle point problem of the form (4) with and . The corresponding projection method based preconditioner is
(44) 
and the preconditioned system is (22) with
(45)  
(46) 
As mentioned before, when boundary condition is purely periodic, the ”=” holds in (45).
By using the connections between the eigenvalues of and a generalized eigenvalue problem, we will reduce the spectral analysis to the analysis of the coefficients of a quadratic polynomial. We have the following theorem for the unsteady Stokes problem.
For the unsteady Stokes problem, the eigenvalues of the preconditioned system have uniform lower and upper bounds, which are independent of the mesh refinement and the parameter . More precisely, . {proof} To investigate the eigenvalues of , we consider the following generalized eigenvalue problem [12].
(47) 
If , then . Substituting it into the first equation of (47), we have
(48) 
For (48), taking the inner product with and dividing by , we obtain the following quadratic polynomial.
We note that satisfies a quadratic polynomial equation and therefore its values can be found by using Vieta’s formulas. To prove has uniform lower and upper bounds, it is sufficient to prove the generalized Rayleigh quotient
has uniform lower and upper bounds.
Firstly, there holds
(49) 
To prove (49), we note that
Using Theorem 3.1 and (40), is symmetric positive semidefinite. Moveover, because , there exists a such that . Therefore,
It follows that
(50) 
Therefore, (49) holds.
To estimate the lower bound, we check
(51)  
(52)  