# An Alternative HSS Preconditioner for the Unsteady Incompressible Navier-Stokes Equations in Rotation Form

Authors: Jia Liu
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We study the preconditioned iterative method for the unsteady Navier-Stokes equations. The rotation form of the Oseen system is considered. We apply an efficient preconditioner
which is derived from the Hermitian/Skew-Hermitian preconditioner to the Krylov subspace-iterative method. Numerical experiments show the robustness of the preconditioned iterative methods with respect to the mesh size, Reynolds numbers, time step, and algorithm parameters. The
preconditioner is efficient and easy to apply for the unsteady Oseen problems in rotation form.


---

## Body

## 1. Introduction

We study the numerical solution methods of the incompressible viscous fluid problems with the following form:(1.1)∂u∂t-νΔu+(u⋅∇)u+∇p=finΩ×(0,Γ],(1.2)∇⋅u=0inΩ×[0,Γ],(1.3)Bu=gon⁡∂Ω×[0,Γ],(1.4)u(x,0)=u0inΩ.Equations (1.1) to (1.4) are also known as the Navier-Stokes equations. Here Ω is an open set of ℝd, where d=2, or d=3, with boundary ∂Ω; the variable u=u(x,t)∈ℝd is a vector-valued function representing the velocity of the fluid, and the scalar function p=p(x,t)∈ℝ represents the pressure. The pressure field, p, is determined up to an additive constant. To uniquely determine p, we may impose some additional condition, such as(1.5)∫Ωpdx=0.The source functionf is given on Ω. Here ν>0 is a given constant called the  kinematic viscosity, which is ν=O(Re-1). Re is the  Reynolds number: Re=VL/ν, where V denotes the mean velocity, and L is the diameter of Ω, see [1]. Also, Δ is the (vector) Laplacian operator in d dimensions, ∇ is the gradient operator, and ∇· is the divergence operator. In (1.3) ℬ is some boundary operator; for example, the Dirichlet boundary condition u=g; Neumann boundary condition ∂u/∂n=g, where n denotes the outward-pointing normal to the boundary, or a mixture of the two.We use fully implicit time discretization and picard linearization to obtain a sequence of Oseen problems, that is, linear problems of the form(1.6)αu-νΔu+(v⋅∇)u+∇p=finΩ,(1.7)∇⋅u=0inΩ,(1.8)Bu=gon⁡∂Ω.Herev is a known divergence-free vector obtained from the previous linearized step (e.g., v=uk). We call the vector v the wind function. In addition, α=O(1/δt) where δt is the time step. Equations (1.6)–(1.8) are referred to as the Oseen problem.We can use either finite element or finite different methods to discretize (1.6)–(1.8). The resulting discrete system 𝒜u=b has the form
(1.9)[ABTB-C][up]=[fg].In this paper, we are interested in an alternative linearization of the steady-state Navier-Stokes equation. Based on the identity(1.10)(u⋅∇)u=12∇(u⋅u)+(∇×u)×u.
In order to linearize it, we replace u in one place with a known divergence-free vector v which can be the solution obtained from the previous Picard iteration. In this case we have
(1.11)(v⋅∇)u≈12∇(u⋅u)+(∇×v)×u.After substituting the right-hand side into (1.6), we find that the corresponding linearized equations have the following form:(1.12)∂u∂t-νΔu+w×u+∇P=finΩ×(0,T],(1.13)∇⋅u=0inΩ×[0,T],(1.14)Bu=gon⁡∂Ω×[0,T],(1.15)u(x,0)=u0inΩ,
where P=p+(1/2)∥u∥22 is the so-called Bernoulli pressure. For the two-dimensional case
(1.16)w×=(0w-w0),
where w=∇×v=-∂v1/∂x2+∂v2/∂x1 is a scalar function.In the three-dimensional case, we have(1.17)w×=(0-w3w2w30-w1-w2w10),
here (w1,w2,w3)=w=∇×v, where wi denotes the ith component of ∇×v. Assume v=(v1,v2,v3), then we have the formal expression of w(1.18)∇×v=|ijk∂∂x∂∂y∂∂xv1v2v3|.Here the divergence-free vector fieldv again denotes the approximate velocity from the previous Picard iteration. Note that when the “wind” function v is irrotational (∇×v=0), (1.12)–(1.14) reduce to the Stokes problem. It is not difficult to see that the linearizations (1.6)–(1.8) and (1.12)–(1.14), although both conservative, are not mathematically equivalent. The momentum equation (1.12) is called the rotation form. We can see that no first-order terms in the velocities appear in (1.12); on the other hand, the velocities in the d scalar equations comprising (1.12) are now coupled due to the presence of the term w×u. The disappearance of the convective terms suggests that the rotation form (1.12) of the momentum equations may be advantageous over the standard form (1.6) from the linear solution point of view. This observation was first made by Olshanskii and his coworkers in [2–5]. In their papers, they showed the advantages of the rotation form over the standard convection form in several aspects.After we discretize the Oseen problem in rotation form (1.12)–(1.14), we obtain the sparse linear system 𝒜x=b, where(1.19)A=[ABTB0].HereA=L+K, where L is the discretization of the operator α-νΔ, and matrix K is the discretization of the term w×, where w=∇×v. In the 2D case,(1.20)K=[0D-D0].
The rectangular matrix B is the discretization of the negative divergence, and BT is the discretization of the gradient.If we use a finite difference method, like Mac and Cell (MAC), see [6], then D is a diagonal matrix where its diagonal elements are the values of w evaluated at the grid edges. Matrix D is a weighted mass matrix if a finite element method is used. In the 3D case, we have
(1.21)K=[0-D3D2D30-D1-D2D10].Again matricesD1,D2, and D3 are all diagonal matrices or weighted mass matrices. Typical sparsity patterns for 𝒜 in the 2D and 3D case are displayed in Figures 1(a) and 1(b).Sparsity patterns for different types of the Oseen problem in rotation form.
(a)
2D case using MAC(b)
3D case using MACFor some discretization methods, a stabilization matrix needs to be added to the (2,2) block of𝒜, namely, a matrix -C, where C is a symmetric positive semidefinite diagonal or scaled mass matrix, or scaled Laplacian with small norm. Figures 2(a) and 2(b) show the sparsity pattern for the coefficient matrix 𝒜 with or without stabilization term in the 2D case. Such a stabilization is not necessary for the MAC discretization.Sparsity patterns for different types of the 2D Oseen problem in convection form.
(a)
without stabilization(b)
with stabilizationTo solve the system𝒜x=b, we can consider the Krylov subspace methods with the preconditioning. Many powerful preconditioning techniques have been explored for the generalized Oseen problems, for example, Uzawa-type preconditioner, block and approximate schur complement preconditioner, pressure preconditioner, and so forth, see [7–11] for more details. However, there is no “best” preconditioner for the saddle point system. To find the “best” preconditioner, we would like to find a preconditioner P, such that the rate of convergence of the preconditioned Krylov subspace matrix is low and bounded independent of the mesh size, viscosity ν and time step α. In addition, the cost of the preconditioning steps must be low. In this paper, we describe such a new preconditioner that satisfies the above requirements in most of cases and demonstrate its utility.A summary of the paper is as follows. Section2 demonstrates the Alternative Hermitian and Skew-Hermitian (AHSS) preconditioner; studies some of its convergence properties and the application of the HSS preconditioner for Krylov subspace methods; Section 3 shows the results of a series of numerical experiments. Finally, section 4 summarizes the approach and future work.


## 2. The Alternative HSS Preconditioner

The alternative HSS preconditioner is based on the nonsymmetric formulation(2.1)[ABT-B0][up]=[f-g].
We have analyzed the advantages of the nonsymmetric formulation in [12]. We gain positive (semi)-definiteness in this case. By changing the sign in front of the (2,1) and (2,2) blocks, we obtain an equivalent linear system with a matrix whose spectrum is entirely contained in the half-plane ℜ(z)>0. (Here we use ℜ(z) to denote the real part of z∈ℂ).The spectra of the nonsymmetric formulation is more friendly to the convergence of Krylov subspace iterations. For an example, GMRES methods, see [13, 14].We have investigated the preconditioner based on the Hermitian and Skew-Hermitian splitting methods for the Navier-Stokes problem, see [12, 15, 16]. However, the HSS preconditioner still has some problems. When the time step is not small enough or the viscosity is relative larger, the iteration number increases a lot. Therefore, we are trying to find another preconditioner which works better than the HSS preconditioner. We find out that if we use a different splitting of the coefficient matrix, we can get a very good results. The following splitting is the new preconditioner we will introduce in the paper. Letting H≡(1/2)(A+AT) and K≡(1/2)(A-AT), we have the following splitting of 𝒜 into two parts:
(2.2)A=[ABT-B0]=[HBT-B0]+[K000].We denote(2.3)H=[HBT-B0],K=[K000].Therefore, we defined the preconditioner as the following:(2.4)Pρ=12ρ(H+In+m)(K+In+m).
Here In+m denotes the identity matrix of order n+m, and ρ>0 is a parameter.Similar in spirit to the classical ADI (alternating-direction implicit) method, we consider the following two splittings of𝒜̂:
(2.5)A=(H+ρI)-(ρI-K),A=(K+ρI)-(ρI-H).
Here ℐ denotes the identity matrix of order n+m. Note that
(2.6)H+ρI=[H+ρInBTBρIm]
is the shifted discretized Stokes problem, where In denotes the identity matrix of order n, and Im denotes the identity matrix of order m. We obtian that
(2.7)K+ρI=[K+ρIn00ρIm]
is nonsingular and has positive definite symmetric part.Alternating between these two splittings leads to the the following iteration:(2.8)(H+ρI)uk+1/2=(ρI-K)uk+b̂,(K+ρI)uk+1=(ρI-H)uk+1/2+b̂,
(k=0,1,…). Here b̂ denotes the right-hand side of (1.9); the initial guess u0 is chosen arbitrarily. Elimination of uk+1/2 from (2.8) leads to a stationary (fixed-point) iteration of the form
(2.9)uk+1=Tρuk+c,k=0,1,…,
where 𝒯ρ=(𝒦+ρℐ)-1(ρℐ-ℋ)(ℋ+ρℐ)-1(ρℐ-𝒦) is the iteration matrix and c:=(𝒮+ρℐ)-1(ρℐ-ℋ)(ℋ+ρℐ)-1(ρℐ-𝒮). The iteration converges for arbitrary initial guesses u0 and right-hand sides b to the solution u*=𝒜̂-1 if and only if ϱ(𝒯ρ)<1, where ϱ(𝒯ρ) denotes the spectral radius of 𝒯ρ.Theorem 2.1.
Consider the problem (2.1), that is,
(2.10)[ABT-B0][up]=[f-g].
We assume that A is positive real, and B has full rank. Then the iteration (2.9) from the splitting (2.3) is unconditionally convergent; that is, ϱ(𝒯ρ)<1 for all ρ>0 and α≥0.Proof.
Consider the splitting (2.3). The iteration matrix 𝒯ρ is similar to
(2.11)T̂ρ≔(ρI-H)(H+ρI)-1(ρI-K)(K+ρI)-1,=RU,
where ℛ:=(ρℐ-ℋ)(ℋ+ρℐ)-1 is symmetric and 𝒰:=(ρℐ-𝒦)(𝒦+ρℐ)-1.
By Kellogg's lemma,∥ℛ∥=∥(ρℐ-ℋ)(ℋ+ρℐ)-1∥≤1 since ℋ is positive semidefinite. 𝒰 is the unitary matrix so ∥𝒰∥=∥(ρℐ-𝒦)(𝒦+ρℐ)-1∥=1. Therefore,
(2.12)ϱ(Tρ)≤1.
We claim thatϱ(𝒯ρ)≠1.
Assume thatλ is one eigenvalue of the preconditioned linear system P-1𝒜x=P-1b. We have
(2.13)Ax=λ12ρ(H+ρI)(K+ρI)x=λ2ρ(HK+ρA+ρ2I)x=λ2ρHKx+λ2Ax+ρλ2x.
Therefore,
(2.14)(1-λ2)Ax=ρλ2(I+1ρ2HK)x.
We claim that 1-λ/2≠0, otherwise, (ℐ+(1/ρ2)ℋ𝒦)x=0. It turns out that
(2.15)[1ρ2HS+I0-BSI]x=0.
However,ρ((1/ρ2)HS+I)=1+ρ((1/ρ2)HS), and HS is orthogonal similar with the matrix H-1/2SH1/2 which is a skew symmetric matrix with only pure imaginary eigenvalues. Thus, if B is full rank, we claim that λ≠2.
We define thatθ=λρ/(2-λ). Since λ≠2, θ is welldefined. Thus, with λ=2θ/(θ+2), we consider the following equation:
(2.16)Ax=θ(I+1ρ2HK)x.
If |λ|<1, then, |2θ/(θ+2)|<1. Since |1-(2θ/(θ+2))|=|(ρ-θ)/(θ+ρ)|≤1, we need to show |ρ-θ|/|ρ+θ|≠1, which means that θ is not a pure imaginary number.
Next we will prove thatθ is not a pure imaginary number.
Consider the system(2.17)[ABT-B0][u-p]=[1ρ2HS+I0-1ρ2BSI][u-p].
We can obtain the following system of the equations:
(2.18)Au+BTp=θu+θρ2HSu,-Bu=-θρ2BSu+θp.
We solvep from the second equation p=(1/θ)B((θ/ρ2)S-I)u. Plug in p into the first equation, we have
(2.19)Au+1ρ2BTBSu-θu-θρ2HSu=1θBTBu.
Applying θuH to the both sides, we can obtain the following equation:
(2.20)θuHAu+θρ2uHBTBSu-θ2uHu-θ2ρ2uHHSu=uHBTBu=‖Bu‖2≥0.
Therefore, θA+(θ/ρ2)BTBS-θ2I-(θ2/ρ2)HS is a Hermitian matrix. Suppose that  Re(θ)=0, that is, θ=ti, where t≠0. We denote that G=θA+(θ/ρ2)BTBS-θ2I-(θ2/ρ2)HS=tiA+(ti/ρ2)BTBS+t2I+(t2/ρ2)HS. While GH=-tiAH+(ti/ρ2)BTBS+t2I+(t2/ρ2)HS=G, which leads to A=AH. A contradiction. Because A is the discretization of the Oseen problem which is not Hermitian.
Thus,θ is not a pure imaginary number.


## 3. Application of the Preconditioner

To solve the preconditionerPαzk=rk, we first solve the system(3.1)Hwk=rk,
for wk, followed by(3.2)Kzk=wk.The first system requires solving systems with coefficient matrixℋ, which is a system from discretized Stokes problem. We have many efficient solvers to solve this type of the system, see [17–19].The second system requires solving the sparse tridiagonal matrixK+ρIn. This can be done by sparse LU factorization, preconditioned GMRES method. Notice that since K+ρI is a tridiagnoal matrix, it is very easy to solve this system. In practice, we can solve (3.1) and (3.2) with inexact solvers. Our experience is that the rate of convergence of the outer Krylov subspace iteration is scarcely affected by the use of inexact inner solves. We can only use 1 step of pcg method for (3.1) and 1 step of gmres method for (3.2).


## 4. Numerical Experiments

In this section we report on several numerical experiments meant to illustrate the behavior of the HSS preconditioner on a wide range of model problems. We consider both Stokes and Oseen-type problems, steady and unsteady, in 2D. The use of inexact solves is also discussed. Our results include iteration counts, as well as some timings and eigenvalue plots.All results were computed in  Matlab 7.6.0 on one processor of an Intel core i7 with 8 GB of memory.In all experiments, a symmetric diagonal scaling was applied before forming the preconditioner, so as to have all the nonzeros diagonal entries of the saddle point matrix equal to 1. We found that this scaling is beneficial to convergence, and it makes finding (nearly) optimal values of the shiftρ easier. Of course, the right-hand side and the solution vector were scaled accordingly. We found, however, that without further preconditioning Krylov subspace solvers converge extremely slowly on all the problems considered here. Right preconditioning was used in all cases.Here we consider linear systems arising from the discretization of the linearized Navier-Stokes equations in rotation form. The computational domain is the unit squareΩ=[0,1]×[0,1]. Homogeneous Dirichlet boundary conditions are imposed on the velocities; experiments with different boundary conditions were also performed, with results similar to those reported below. We experimented with different forms of the divergence-free field v appearing (via w=∇×v) in the rotation form of the unsteady Oseen problem. Here we present results for the choice w=16x(x-1)+16y(y-1) (2D case) and w=(-4z(1-2x),4z(2y-1),2y(1-y)) (3D case). The 2D case corresponds to the choice of v. The equations were discretized with a Marker-and-Cell (MAC) scheme with a uniform mesh size h. The outer iteration (full GMRES) was stopped when (4.1)‖rk‖2‖r0‖2<10-6,
where rk denotes the residual vector at step k. For the results presented in this section, we use the exact solver to solve the two systems.In Tables1, 2, and 3, we present results for unsteady problems with α=10 to α=100 for different values of ν. We can see from the table that AHSS preconditioning results in fast convergence in all cases, and that the rate of convergence is virtually h-independent. Here as in all other unsteady (or quasi-steady) problems that we have tested, the rate of convergence is not overly sensitive to the choice of ρ, especially for small ν. A good choice is ρ≈0.01 for the two most grids.Table 1
Results for 2D unsteady Oseen problem different values ofα (exact solves) and ν=0.1.
Gridα=10α=20α=50α=1008×8754316×16864432×32865464×648654128×1288654Table 2
Results for 2D unsteady Oseen problem different values ofα (exact solves) and ν=0.05.
Gridα=10α=20α=50α=1008×8443316×16653332×32864364×649755128×12810755Table 3
Results for 2D unsteady Oseen problem different values ofα (exact solves) and ν=0.01.
Gridα=10α=20α=50α=1008×8332216×16433232×32543364×648533128×1289643


## 5. Conclusions

In this paper, we have considered preconditioned iterative methods applied to discretizations of the Navier-Stokes equations in rotation form. We focus on the unsteady case of the linearized Navier-Stokes problem. We have compared the performance of the alternative HSS (AHSS) preconditioners with regard to the mesh size, the Reynolds number, the time step, and other problem parameters. We find that the AHSS preconditioner has a robust behavior especially for unsteady Oseen problems. Although our computational experience has been limited to uniform MAC discretizations and simple geometries, the preconditioner should be applicable to more complicated problems and discretizations, including unstructured grids.Compared with HSS (Hermitian and Skew-Hermitian preconditioner), the AHSS preconditioner works better for relative large viscosity. For an example,ν>0.05. For the smaller viscosity, ν<0.01, HSS preconditioner will be recommended.In the future study, we will investigate the performance the AHSS preconditioner based using the inexact solvers for the inner iteration. Also the picard's iteration will be tested.



---

*Source: 307939-2012-04-05.xml*