A Splitting Preconditioner for the Incompressible Navier–Stokes Equations

In this paper, a splitting preconditioner based on the relaxed dimensional factorization (RDF) preconditioner and the modified augmented Lagrangian (MAL) preconditioner for the incompressible Navier–Stokes equations is presented. The preconditioned matrix is analyzed, and similar results arising from the RDF and the MAL preconditioners are obtained. The corresponding details of the spectrum analysis are given. Finally, we compare the three preconditioners and numerical experiments are implemented by using the IFISS package.


Introduction
We consider the following incompressible Navier-Stokes equations describing the flow of viscous Newtonian fluids: where Ω ⊂ R d (d = 2, 3) is an open bounded domain with boundary ∂Ω, [0, T ] is the time interval, the unknown velocity fields u(x, t) and pressure fields p(x, t), ν is the kinematic viscosity, ∆ is the vector Laplacian operator, ∇ is the gradient operator, div is the divergence, f , g and u 0 are given functions. After implicit time discretization and linearization of the Navier-Stokes system of equations by Picard fixed-point iteration, we get a sequence of the Oseen problems. Discretization of the Oseen problems using finite element methods results in a sequence of large sparse linear systems of equations. These equations are expressed as with u and p representing the discrete velocity and pressure, respectively. A denotes the discretization of the diffusion, convection and time-dependent terms.
A is a diagonal block matrix, e.g. A = A1 0 0 A2 in 2D. B T is the discrete gradient, B denotes the (negative) discrete divergence, C is a stabilization matrix and depends on the discretization stability condition, f and g contain the forcing and boundary terms. If we use the LBB-stable finite element to discretize this problem and use a simple transformation J = I A 0 0 −I C , where I A and I C are the identity matrices, then (1.5) can be rewritten in the mathematically equivalent system as where the spectrum of the coefficient matrix of (1.6) is entirely contained in the right half plane (see [3]). These systems can be solved by direct methods, but they require extensive resources in terms of computational time and memory. For 3D and large 2D problems, iterative methods, in combination with suitable preconditioners, are the methods of choice. Benzi et. al have presented the RDF preconditioner to effectively solve the system (1.6) in [5]. The RDF preconditioner is one kind of dimensional splitting (DS) preconditioners [4,8,10] and is more effective than the preconditioner derived in [4]. The experimentally optimal parameter within the RDF preconditioner always strongly depends on mesh size. Although the parameter values can be estimated by Fourier Analysis, see [5], they become smaller and smaller for finer grids and cannot perfectly approximate the experimentally optimal values. However, according to the structure of the RDF preconditioner, practical implementation of the RDF preconditioner is good. The augmented Lagrangian (AL) preconditioner has been presented by Benzi et. al and it is based on the block triangular preconditioner for the augmented system of (1.6), see [3,6,7,9,14,15].
For the AL preconditioner, approximating the Schur complement BA −1 B T is relatively easy, and the main issue is how to effectively compute solution of linear systems associated with the augmented block; see Section 4. The modified AL preconditioner is presented to decrease computational cost of the AL preconditioner, but it may be an inferior position comparing with the RDF preconditioner because more matrix-vector multiplications within the augmented system make costs increase. However, the MAL preconditioner retains a portion of features of the AL preconditioner. The parameter of the MAL preconditioner is almost steady and robust, it is largely insensitive to problem parameters including grid size, viscosity ν, non-uniform meshes, etc. The confidence interval of this parameter is wide, the optimal parameter value obtained by Fourier analysis [5,9] can be close to the experimentally optimal value of the MAL preconditioner. Our starting point is to investigate the sensitivity of the parameter of these preconditioners, and to search a preconditioner satisfying these advantages. That is, we consider whether there is a preconditioner which can be combined with both RDF and MAL preconditioners. In this paper we propose a splitting preconditioner, whose ideas come from the RDF preconditioner and the MAL preconditioner. Moreover, we analyze the spectrum distribution of the corresponding preconditioned matrix and the choice of the optimal parameter. The rest of the paper is organized as follows. In Sections 2 and 3, we briefly introduce the RDF preconditioner and the MAL preconditioner respectively. In Section 4, we present our splitting preconditioner and derive some properties. Then we show how Fourier analysis can be used to select the parameter. In Section 5, we present results of numerical experiments including comparisons with the RDF and the MAL preconditioners. Finally, some conclusions are drawn in Section 6.

The RDF Preconditioner
In this section, we make a brief description of the RDF preconditioner; for further details, see [4,5]. For simplicity, we only consider the 2D case. The reformulation of (1.6) is expressed as The RDF preconditioner in [5] is defined as follows where τ is the positive parameter. The RDF preconditioner is also scaled, a scaling is applied to the coefficient matrix before forming the preconditioner. The behaviour of RDF preconditioning can be improved by diagonal scaling. Unless otherwise specified, we always perform a preliminary symmetric scaling of the system Hx = b in the form D − 1 2 HD − 1 2 y = D − 1 2 b, with y = D 1 2 x and D = diag(D 1 , D 2 , I m ), where diag(D 1 , D 2 ) is the main diagonal of the velocity mass matrix for 2D problems. Obviously, this diagonal scaling is regarded as a simple preconditioner.

The Modified Augmented Lagrangian Preconditioner
For the steady-state problems (1.6), the equivalent AL formulation [14,15] is given by Thus the ideal AL system and the ideal AL preconditioner [6,7,9] are whereŜ al is the approximate Schur complement. In order to retain the transformation property of the system (1.6), in this paper, the proposed P al -type AL preconditioner is the block triangular matrix as follows whereŜ is the approximate Schur complement and usually implicitly defined byŜ where M p denotes the approximate pressure mass matrix, ν is the viscosity. A good choice of W is the pressure mass matrixM p , see [6,14,16]. In practice, we always use the main diagonal of the pressure mass matrix. In many cases, W is also replaced by M p . For decreasing calculation cost, in practice, M p is a diagonal matrix or is replaced by spectrally equivalent diagonal matrix. Considering 2D case, we have A = diag(A 1 , A 2 ) and B = ( Then the (1, 2) block is dropped, the modified augmented Lagrangian preconditioner is obtained. Readers can refer to [6,7,9] for further details including the choice of matrix W . Throughout this paper, we always consider the MAL preconditioner.

A Splitting Preconditioner
In this section, we define a preconditioner via the splitting scheme. Our splitting preconditioner is formulated as follows where W x is SPD and α > 0, By comparing (4.1) with (2.2) and (3.2), it is shown that the block structure of these preconditioners is similar. We can find that W x and W −1 x are added into the RDF preconditioner on the one hand, then can get (4.1); on the other hand, dealing with the MAL preconditioner also yields (4.1). Obviously, (4.1) is not a type of the dimensional factorization preconditioner, it is just a generic splitting type. Based on this relationship, we consider W x = W , and don't pursue other choices of W x . Now we analyze some properties of new preconditioner. Firstly, we consider the spectral properties of the new preconditioned matrix HP −1 , then we resort into Fourier analysis for guiding in the choice of the parameter α. More details which are similar to properties of the RDF preconditioner and the MAL preconditioner are given; see [5,7,9].
According to the special structure of (4.1), the new preconditioner can be factorized into some factors. It is convenient to implement numerical experiments and analyze its properties. In [5,17], there are some well known results. We extend those propositions for the preconditioned matrix of (4.1) as follows.

Lemma 2. Assume
Then V 22 has zero eigenvalues with multiplicity at least n 2 , and the remaining eigenvalues Proof. Referring to [5], this proposition can be easily shown .
Theorem 1. The preconditioned matrix HP −1 has eigenvalues equal to 1 with multiplicity at least n 1 + n 2 , the remaining eigenvalues are the eigenvalues λ i of the matrix Z α .
Proof. Referring to [5], this theorem can be easily proved.
We note that where the eigenvalues µ i satisfy the generalized eigenvalue problem Proof. Applying the Sherman-Morrison-Woodbury formulation toŜ k giveŝ Then we get similar results to [5].
These propositions show that there exist inseparable relationships between the RDF preconditioner and our splitting preconditioner. Assume all the idealized conditions are satisfied; see the parameter analysis part of this section. The matrices A 1 , A 2 , B 1 , B 2 and W are all diagonalizable by the discrete Fourier transform, i.e., Here a i , b i , w are vectors containing the eigenvalues of the corresponding matrices, and the unitary matrix U is composed of Fourier components e i2πhθk √ l with i = √ −1, k = 1, 2, . . . , l, θ = 1, 2, . . . , l and grid size h = 1/l, see [5]. Refer to the RDF preconditioner, let the matricesŜ 1 ,Ŝ 2 and Z α can be diagonalized by U , then the eigenvalues of Z α can be simply expressed as where s 1 , s 2 are the eigenvalues ofŜ 1 ,Ŝ 2 and are respectively given by Readers can refer to [5] for further details. According to [5], we can compare the corresponding eigenvalues formulation with respect to the RDF preconditioner, i.e., z τ = 1 . Let α ≤ 1, τ ≤ 1 and w ≤ 1, it is not hard to see that the magnitude of 1 τ and 1 τ 2 increase more quickly than α and α 2 . Then in conjunction with the above equation, we also can find the corresponding optimal value of (α, w) and τ so that the eigenvalues z α have a more clustering effect than the eigenvalues z τ , the worst case is that the clustering effect of the both eigenvalues is the same.
We also propose another method to analyze the spectrum of the preconditioned matrix in order to observe the relationship between our preconditioner and the MAL preconditioner. Z α is expressed as , then Since proofs of Lemmas 3 and 4 are easy, readers can refer to [6]. The conditions of Lemma 4 are satisfied if we assume that B has full row rank and A α is positive definite. Hence the remaining m eigenvalues λ i are solutions of the generalized eigenproblem αBÂ −1 Hence, the non-unit eigenvalues of the preconditioned matrix are given as follows: . Now we analyze boundary of λ i , some corresponding results and relations are obtained. Writing λ i = a λ + ib λ andμ i = a µ + ib µ , it is easy to obtain the following expressions of the real and the imaginary parts of λ i : Theorem 4. The remaining m eigenvalues λ i are given by Theorem 2, wherê µ i = a µ + ib µ satisfies (4.2). The following estimates in [2,6,13]: are also valid for the eigenvalues of the matrix Z α .
It is shown that our preconditioner contains some features of the MAL preconditioner. We consider the matrix , but the purely imaginary eigenvalues of K and the parameter factors may lead to the eigenvalues of A α more decentralized, i.e., the eigenvalues of the corresponding preconditioned matrix have less clustering effect. According to analysis (Theorem 3 and [7]) and experiments, the eigenvalues with respect to our preconditioner tend to less clustering effect than the MAL preconditioner, this feature may cause the convergence rates of iterative method become slow. In Figures 1 and 2, the conclusion of the foregoing analysis are confirmed. The SPP denotes our splitting preconditioner, and we show the spectrum plots of the preconditioned matrices obtained from the RDF, the SPP , the MAL and the AL preconditioners under the experimentally optimal values of τ , α, and γ respectively. We can see that the eigenvalues of the preconditioned matrix of the new preconditioner tend to more clustering effect compared with the RDF preconditioner, clustering effect of the eigenvalues of the preconditioned matrix of the AL preconditioner is the best. From plots, we can see that our preconditioner and the MAL preconditioner are slightly similar, this explains the fact that the convergence features of our preconditioner tend to the MAL preconditioner. We note that the following numerical experiments verify this point. Thus our preconditioner is a compromise preconditioner. Allowing α → 0+ or α → ∞ doesn't make the norm H − P 2 very small, then it suggests that there is an optimal α for our preconditioner. Hence, we resort to Fourier analysis (FA) for approximating the optimal value of the parameter α; see [5,9]. As usual, the use of this technique needs some rather drastic simplifications and assumptions about the problem. We assume that the Oseen problem has constant coefficients, it is defined on the unit square with periodic boundary conditions, and is discretized on a uniform, with grid size h = 1/l. Moreover, we assume the matrices A 1 , A 2 , B 1 , B 2 and W are all square and commute. Though these assumptions are virtually never met in real problems, we emphasize that these assumptions are made to guide the choice of the parameter, and the parameter obtained by Fourier analysis often gives good results for more general problems [5,9,11].
To further simplify, we resort to the symbols of the corresponding operators. The block matrix A i is a discrete scalar convection-diffusion-reaction operator we mainly consider νL+N i operator, L is discrete Laplacians and N i is discrete convection operators. Then the discrete 2D steady state convection-diffusion operator ℵ = I l ⊗ (νL x + N x ) + (νL y + N y ) ⊗ I l , where ⊗ denotes the tensor product, I l is the identity matrix with l order, L x and L y are discrete 1D Laplacians, N x and N y are discrete 1D convection operators in the x and y directions. Similarly, the matrix B i represent discrete partial derivatives with respect to x and y. Noting the discretization of the ordinary derivatives d dx and d dy by R x and R y , then B 1 = I l ⊗ R x and B 2 = R y ⊗ I l . Let us assume that diffusion and convection terms are discretized by centered finite differences and the divergence is discretized by one-sided differences, and observe that W scales as h 2 , then we get the following correspondence between the operators and their symbols, θ = (θ x , θ y ): We note that the symbol respect the scaling of the matrices discretized by finite element methods, then L, N , R can be expressed as diagonal matrices, whose diagonal entries are the corresponding eigenvalues [9]. Hence matrices A i and B i can be represented by symbols as well. Moreover, we express B 1 A −1 , respectively. From [9], matrix Z α = α(S 1 W −1 + S 2 W −1 ) − 2α 2 S 1 W −1 S 2 W −1 can be expressed as then the eigenvalues of Z α are given by Therefore, we want the eigenvalues of Z α to be around 1, it is equivalent to solving the following optimization problem min α>0 λ(Z α ) − 1 , subject to θ = 1, 2, . . . , l.
Comparing with the choice of the parameter in the MAL preconditioner [9], it is shown that these both minimum problems are equivalent. Thus we can turn to the choice of the parameter of the MAL preconditioner to obtain the parameter of our splitting preconditioner.

Numerical Experiments
In this section, we shall present numerical experiments for the linear systems coming from the two dimensional Oseen models of the incompressible flow discretized by finite element methods to test the performance of our preconditioner. Two main test problems are generated by IFISS software package [12]: Lid-driven cavity (LCD) problem and backward facing-step (BFS) problem. These experiments were performed in MATLAB 7 on an Intel (R) Pentium (R) 4 CPU with 3.00 GHz and 1GB of memory. We should mention using restarted GMRES [18] as the Krylov subspace solver, stopped when the relative residual norm is reduced below 10 −6 , the maximum subspace dimension is set to 20, in the tests it always uses zero initial guesses. Unless otherwise specified, the corresponding parameter of discretization is default in the IFISS software package. We also set W = M p = diag(M p ), whereM p denotes the pressure mass matrix in our numerical experiments, exact solves are performed on the subsystems involvingÂ 1 andÂ 2 by means of direct LU factorization (for Oseen problems) after proceeding by AMD reordering technique [1,5,7,9].
We consider the 2D leaky-lid driven cavity problem discretized by Q2-Q1 and Q2-P1 finite elements on the uniform and stretched grids. We are interested in observing the effect of the backward facing step problem on a non-square domain. In the given tables and figures, RDF denotes the RDF preconditioner, MAL is the modified AL preconditioner, SPP denotes our splitting preconditioner.
Example 1. The steady cavity Oseen problem. In this example, we perform the diagonal scaling technique for all preconditioning. In Tables 1-5, we present preconditioned GMRES iteration counts for the leaky lid driven cavity problem discretized by Q2-Q1 and Q2-P1 finite elements. In Table 1, we use experimental optimal parameter values to obtain the corresponding iteration counts and timings. We can see that the iteration counts of the three preconditioners are more or less similar (just a difference of 1-2 iterations in many cases), thus this data is not sufficient to conclude which preconditioner is better in the case of uniform grids. However, at the same time, it is shown that the cost of the MAL preconditioner is much larger than the other two preconditioners from Table  1, and is about three times the cost of the RDF or the SPP preconditioner, here Setup time mainly shows the time of direct LU factorization of two linear subsystems. Indeed, see [5], the RDF preconditioner can be factorized into four simple matrix factors, then it requires solving two linear systems at each iterative step. Similarly, our preconditioner can be factored as follows; see Lemma 1: At each step, the SPP preconditioner requires solving two linear systems with the coefficient matricesÂ 1 =A 1 +αB T 1 W −1 B 1 andÂ 2 =A 2 +αB T 2 W −1 B 2 . If necessary, we also perform diagonal scaling to improve the convergence rates. Compared with the RDF preconditioner, the SPP preconditioner just additionally needs some matrix W -vector multiplications, if matrix W is simple, e.g., the diagonal matrix, then we can neglect their computation. This feature is reflected in Table 1. For the MAL preconditioner, the coefficient matrix requires more matrix-vector (M-V) multiplications obviously during iterative process, i.e., O(( where v is a vector and O indicates computational complexity. In GMRES algorithm, there are two M-V multiplication, one is initial residual obtained and the other is in the orthogonalization. Table 2 shows the average CPU times of matrix-vector multiplications, carrying out the preconditioner and orthogonalization for a single iterative process. We can find that number of nonzeros of the coefficient matrix A + αB T W −1 B are much more, B T W −1 B makes sparsity of the coefficient matrix smaller. This may be more algebraic operation so that the computation time of total iteration increases. Since the linear subsystems are similar for all three preconditioners, Table 1 has shown the more or less costs for using direct LU factorization method. In Table 2, we can easily find that the coefficient matrix of the augmented system lead to the more costs.  In Figures 3 and 4, the sensitivity of the parameter is shown for different mesh sizes and the viscosity. We can see that the sensitivity of the parameter of the MAL and the SPP preconditioners are similar, and have the same change tendency. At the same times, it's shown that the parameter of the MAL and the SPP preconditioners are robust. According to the choice of the parameter of our preconditioner and [9], in Tables 3 and 4, we know that the parameter values of the SPP preconditioner and the MAL preconditioner are equivalent, which are obtained by Fourier Analysis, i.e. γ = α. For the mild viscosity, the SPP preconditioner and the MAL preconditioner are better than the RDF preconditioner, the convergence rate of GMRES of the MAL preconditioner is faster than the SPP preconditioner for the small viscosity. In Figure 5, the   parameters τ , α, γ obtained by Fourier Analysis and found experimentally are illustrated respectively. The parameters α and γ obtained by Fourier Analysis can work well. We can see that the SPP preconditioner has almost the same effect of the MAL preconditioner without the relatively slow convergence rate. In Table 5, we display the number of GMRES iterations for steady problem discretized by Q2-P1 element on a stretched grid, the experimental optimal values of the parameter τ, α, γ are presented respectively in parentheses. It is shown that the RDF preconditioner converges slowly while the MAL precondi- The parameter FA γ,α optimal γ optimal α FA τ optimal τ (b) Figure 5. The parameters τ , α, γ change curves with respect to FA optimal and experimentally optimal vs. the mesh size in the RDF, SPP, and MAL preconditioners for steady Oseen problem (Q2-Q1, uniform grids).   tioner and the SPP preconditioner converge quickly, and the convergence rate of the MAL preconditioner is slightly faster than that of the SPP preconditioner.
From Tables 1-5, we can see that the performance of the SPP preconditioner is good for steady problems, it has the similar convergence behaviour with the MAL preconditioner, and the similar structure with the RDF preconditioner.
Example 2. The unsteady cavity Oseen problem. Similarly, we present the unsteady problems discretized by Q2-Q1 element on a stretched grid. Linear systems of this type tend to be easier to solve than the ones arising from the steady case, since the presence of the additional positive definite term σM matrix makes the A block more diagonally dominant, where M is the velocity mass matrix; the parameter σ ≥ 0 is typically proportional to the reciprocal of the time step, and is zero in the steady case. In our following experiments we also let σ = 1/h, where h is the mesh size (see [5,7]).
Referring to [7], since the MAL preconditioner can also work well when the simple choice of the parameter γ = 1, thus in this example we also set the same parameter value in order to observe the convergence behavior of the SPP preconditioner. In Table 9, we show the rates of convergence for these preconditioners (RDF, MAL, and SPP). It is shown that the convergence rates of the RDF preconditioner are mildly dependent of mesh size, but they are independent of the viscosity when the parameter is obtained experimentally. Table 9. Number of GMRES(20) iterations with the RDF, MAL and SPP preconditioners for unsteady Oseen problems (Q2-Q1 FEM stretched grids). The parameter τ is experimentally optimal, the parameter α, γ equal 1, the Such complementŜ = (ν + γ) −1 W .  (1) 26(1) 25(0.0008) 55(1) 55(1) The MAL preconditioner and the SPP preconditioner work well in the case of ν = 0.1, α = γ = 1 via using the diagonal scaling technique. However, the convergence rates of the SPP and the MAL preconditioners mildly depend on the mesh size if the parameters equal one. The smaller ν becomes, the more noticeable the deterioration gets. Nevertheless, it can be still accepted. Here we use the implicit Schur complementŜ −1 = (ν +γ)W −1 instead ofŜ = γ −1 W to improve the convergence behavior of GMRES for the MAL preconditioner.
From Tables 6 to 9, considering the convergence rate of preconditioner for the unsteady problems, the SPP preconditioner works better than the MAL preconditioner on uniform grids, and it is similar to the MAL preconditioner on stretched girds. If the approximate parameter values can be precise enough to tend to the experimentally optimal values, then the RDF preconditioner presents the best performance. Otherwise, the SPP preconditioner may be a nice choice for the unsteady problems. ied its behaviour for Navier-Stokes equations. In addition, we have compared our preconditioner with the RDF preconditioner and the MAL preconditioner, and have used Fourier analysis tool for guiding the choice of the parameter. For steady problems, our preconditioner show out the similar behaviour to the MAL preconditioner, and our preconditioner and the MAL preconditioner converge more quickly than the RDF preconditioner, especially on stretched grids. Though the convergence rate of the MAL preconditioner is slightly faster than our preconditioner, the computational cost of our preconditioner is much less than the MAL preconditioner, and is more or less compared with the RDF preconditioner. Meanwhile, Our preconditioner has the features of the RDF and the MAL preconditioners, but it more tends to be the MAL preconditioner. For unsteady case, our preconditioner works very well on uniform grids. Despite the convergence behaviour of our preconditioner mildly depends on the mesh size on stretched grids, these results are still acceptable, the performance of our preconditioner is better than the MAL preconditioner in many circumstances. All experimental results show that our preconditioner still has a certain competitiveness.
Finally, using direct methods for the solution of inner linear systems is certainly feasible, but it is not a good idea to solve larger 2D or 3D problems due to memory and time constraints, then exact solvers can be replaced with inexact solvers. Therefore, inexact solve method and 3D case need to be further considered in the future.