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

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.

1. Introduction

We study the numerical solution methods of the incompressible viscous fluid problems with the following form: 𝜕𝐮()(],[][]𝐮()𝜕𝑡−𝜈Δ𝐮+𝐮⋅∇𝐮+∇𝑝=𝐟inΩ×0,Γ∇⋅𝐮=0inΩ×0,Γ,ℬ𝐮=𝐠on𝜕Ω×0,Γ,𝐱,0=𝐮₀inΩ. (1.1)(1.2)(1.3)(1.4)

Equations (1.1) to (1.4) are also known as the Navier-Stokes equations. Here Ω is an open set of ℝ^𝑑, where 𝑑=2, or 𝑑=3, with boundary 𝜕Ω; the variable 𝐮=𝐮(𝐱,𝑡)∈ℝ^𝑑 is a vector-valued function representing the velocity of the fluid, and the scalar function 𝑝=𝑝(𝐱,𝑡)∈ℝ represents the pressure. The pressure field, 𝑝, is determined up to an additive constant. To uniquely determine 𝑝, we may impose some additional condition, such as_Ω𝑝𝑑𝐱=0. (1.5)

The source function 𝐟 is given on Ω. Here 𝜈>0 is a given constant called the kinematic viscosity, which is 𝜈=𝑂(Re⁻¹). Re is the Reynolds number: Re=𝑉𝐿/𝜈, where 𝑉 denotes the mean velocity, and 𝐿 is the diameter of Ω, see [1]. Also, Δ is the (vector) Laplacian operator in 𝑑 dimensions, ∇ is the gradient operator, and ∇⋅ is the divergence operator. In (1.3) ℬ is some boundary operator; for example, the Dirichlet boundary condition 𝐮=𝐠; Neumann boundary condition 𝜕𝐮/𝜕𝐧=𝐠, where 𝐧 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 𝛼𝐮−𝜈Δ𝐮+(𝐯⋅∇)𝐮+∇𝑝=𝐟inΩ,∇⋅𝐮=0inΩ,ℬ𝐮=𝐠on𝜕Ω. (1.6)(1.7)(1.8)

Here 𝐯 is a known divergence-free vector obtained from the previous linearized step (e.g., 𝐯=𝐮_𝑘). We call the vector 𝐯 the wind function. In addition, 𝛼=𝑂(1/𝛿𝑡) where 𝛿𝑡 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 𝒜𝐮=𝐛 has the form 𝐀𝐁^𝑇𝐮𝑝=𝐟𝐁−𝐂𝑔. (1.9)

In this paper, we are interested in an alternative linearization of the steady-state Navier-Stokes equation. Based on the identity ()1𝐮⋅∇𝐮=2∇(𝐮⋅𝐮)+(∇×𝐮)×𝐮. (1.10) In order to linearize it, we replace 𝐮 in one place with a known divergence-free vector 𝐯 which can be the solution obtained from the previous Picard iteration. In this case we have ()1𝐯⋅∇𝐮≈2∇(𝐮⋅𝐮)+(∇×𝐯)×𝐮. (1.11)

After substituting the right-hand side into (1.6), we find that the corresponding linearized equations have the following form:𝜕𝐮(],[][]𝐮()𝜕𝑡−𝜈Δ𝐮+𝐰×𝐮+∇𝑃=𝐟inΩ×0,𝑇∇⋅𝐮=0inΩ×0,𝑇,ℬ𝐮=𝐠on𝜕Ω×0,𝑇,𝐱,0=𝐮₀inΩ, (1.12)(1.13)(1.14)(1.15) where 𝑃=𝑝+(1/2)‖𝐮‖²₂ is the so-called Bernoulli pressure. For the two-dimensional case 𝐰×=0𝑤−𝑤0, (1.16) where 𝑤=∇×𝐯=−𝜕𝑣₁/𝜕𝑥₂+𝜕𝑣₂/𝜕𝑥₁ is a scalar function.

In the three-dimensional case, we have 𝐰×=0−𝑤₃𝑤₂𝑤₃0−𝑤₁−𝑤₂𝑤₁0, (1.17) here (𝑤₁,𝑤₂,𝑤₃)=𝐰=∇×𝐯, where 𝑤_𝑖 denotes the 𝑖th component of ∇×𝐯. Assume 𝐯=(𝑣₁,𝑣₂,𝑣₃), then we have the formal expression of 𝑤𝜕∇×𝐯=𝐢𝐣𝐤𝜕𝜕𝑥𝜕𝜕𝑦𝑣𝜕𝑥₁𝑣₂𝑣₃. (1.18)

Here the divergence-free vector field 𝐯 again denotes the approximate velocity from the previous Picard iteration. Note that when the “wind” function 𝐯 is irrotational (∇×𝐯=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 𝑑 scalar equations comprising (1.12) are now coupled due to the presence of the term 𝐰×𝐮. 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 𝒜𝐱=𝐛, where𝒜=𝐀𝐁^𝑇𝐁𝟎. (1.19)

Here 𝐀=𝐋+𝐊, where 𝐋 is the discretization of the operator 𝛼−𝜈Δ, and matrix 𝐊 is the discretization of the term 𝐰×, where 𝐰=∇×𝐯. In the 2D case,𝐊=𝟎𝐃−𝐃𝟎. (1.20) The rectangular matrix 𝐁 is the discretization of the negative divergence, and 𝐁^𝑇 is the discretization of the gradient.

If we use a finite difference method, like Mac and Cell (MAC), see [6], then 𝐃 is a diagonal matrix where its diagonal elements are the values of 𝐰 evaluated at the grid edges. Matrix 𝐃 is a weighted mass matrix if a finite element method is used. In the 3D case, we have 𝐊=𝟎−𝐃₃𝐃₂𝐃₃𝟎−𝐃₁−𝐃₂𝐃₁𝟎. (1.21)

Again matrices 𝐃₁,𝐃₂, and 𝐃₃ 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).

(a) 2D case using MAC

(b) 3D case using MAC

Figure 1 

Sparsity patterns for different types of the Oseen problem in rotation form.

For some discretization methods, a stabilization matrix needs to be added to the (2,2) block of 𝒜, namely, a matrix −𝐂, where 𝐂 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.

(a) without stabilization

(b) with stabilization

Figure 2 

Sparsity patterns for different types of the 2D Oseen problem in convection form.

To solve the system 𝒜𝐱=𝐛, 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 𝑃, 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. Section 2 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) 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 ℜ(𝑧)>0. (Here we use ℜ(𝑧) to denote the real part of 𝑧∈ℂ).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 𝐻≡(1/2)(𝐴+𝐴^𝑇) and 𝐾≡(1/2)(𝐴−𝐴^𝑇), we have the following splitting of 𝒜 into two parts: 𝒜=𝐀𝐁^𝑇=−𝐁𝟎𝐇𝐁^𝑇+−𝐁𝟎𝐊𝟎𝟎𝟎. (2.2)

We denoteℋ=𝐇𝐁^𝑇−𝐁𝟎,𝒦=𝐊𝟎𝟎𝟎. (2.3)

Therefore, we defined the preconditioner as the following:𝑃_𝜌=12𝜌ℋ+𝐼_𝑛+𝑚𝒦+𝐼_𝑛+𝑚. (2.4) Here 𝐼_𝑛+𝑚 denotes the identity matrix of order 𝑛+𝑚, 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) Here ℐ denotes the identity matrix of order 𝑛+𝑚. Note that ℋ+𝜌ℐ=𝐇+𝜌𝐈_𝑛𝐁^𝑇𝐁𝜌𝐈_𝑚 (2.6) is the shifted discretized Stokes problem, where 𝐼_𝑛 denotes the identity matrix of order 𝑛, and 𝐼_𝑚 denotes the identity matrix of order 𝑚. We obtian that 𝒦+𝜌ℐ=𝐊+𝜌𝐈_𝑛𝟎𝟎𝜌𝐈_𝑚 (2.7) is nonsingular and has positive definite symmetric part.

Alternating between these two splittings leads to the the following iteration: ℋ+𝜌ℐ𝐮_𝑘+1/2=𝜌ℐ−𝒦𝐮_𝑘+𝐮𝐛,𝒦+𝜌ℐ_𝑘+1=𝜌ℐ−ℋ𝐮_𝑘+1/2+𝐛, (2.8) (𝑘=0,1,…). Here 𝐛 denotes the right-hand side of (1.9); the initial guess 𝐮₀ is chosen arbitrarily. Elimination of 𝐮_𝑘+1/2 from (2.8) leads to a stationary (fixed-point) iteration of the form 𝐮_𝑘+1=𝒯_𝜌𝐮_𝑘+𝐜,𝑘=0,1,…, (2.9) where 𝒯_𝜌=(𝒦+𝜌ℐ)⁻¹(𝜌ℐ−ℋ)(ℋ+𝜌ℐ)⁻¹(𝜌ℐ−𝒦) is the iteration matrix and 𝐜∶=(𝒮+𝜌ℐ)⁻¹(𝜌ℐ−ℋ)(ℋ+𝜌ℐ)⁻¹(𝜌ℐ−𝒮). The iteration converges for arbitrary initial guesses 𝐮₀ and right-hand sides 𝐛 to the solution 𝐮_∗=𝒜⁻¹ if and only if 𝜚(𝒯_𝜌)<1, where 𝜚(𝒯_𝜌) denotes the spectral radius of 𝒯_𝜌.

Theorem 2.1. Consider the problem (2.1), that is, 𝐀𝐁^𝑇𝐮𝑝=𝐟−𝐁𝟎−𝑔. (2.10) We assume that 𝐀 is positive real, and 𝐁 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) where ℛ∶=(𝜌ℐ−ℋ)(ℋ+𝜌ℐ)⁻¹ is symmetric and 𝒰∶=(𝜌ℐ−𝒦)(𝒦+𝜌ℐ)⁻¹.
By Kellogg's lemma, ‖ℛ‖=‖(𝜌ℐ−ℋ)(ℋ+𝜌ℐ)⁻¹‖≤1 since ℋ is positive semidefinite. 𝒰 is the unitary matrix so ‖𝒰‖=‖(𝜌ℐ−𝒦)(𝒦+𝜌ℐ)⁻¹‖=1. Therefore, 𝜚𝒯_𝜌≤1. (2.12)
We claim that 𝜚(𝒯_𝜌)≠1.
Assume that 𝜆 is one eigenvalue of the preconditioned linear system 𝑃⁻¹𝒜𝑥=𝑃⁻¹𝐛. We have 1𝒜𝐱=𝜆𝐱=𝜆2𝜌ℋ+𝜌ℐ𝒦+𝜌ℐ2𝜌ℋ𝒦+𝜌𝒜+𝜌²ℐ𝐱=𝜆𝜆2𝜌ℋ𝒦𝐱+2𝒜𝐱+𝜌𝜆2𝐱. (2.13) Therefore, 𝜆1−2𝒜𝐱=𝜌𝜆21ℐ+𝜌²ℋ𝒦𝐱. (2.14) We claim that 1−𝜆/2≠0, otherwise, (ℐ+(1/𝜌²)ℋ𝒦)𝐱=0. It turns out that 1𝜌²𝐇𝐒+𝐈𝟎−𝐁𝐒𝐈𝐱=0. (2.15)
However, 𝜌((1/𝜌²)𝐇𝐒+𝐈)=1+𝜌((1/𝜌²)𝐇𝐒), and 𝐇𝐒 is orthogonal similar with the matrix 𝐇^−1/2𝐒𝐇^1/2 which is a skew symmetric matrix with only pure imaginary eigenvalues. Thus, if 𝐁 is full rank, we claim that 𝜆≠2.
We define that 𝜃=𝜆𝜌/(2−𝜆). Since 𝜆≠2, 𝜃 is welldefined. Thus, with 𝜆=2𝜃/(𝜃+2), we consider the following equation: 1𝒜𝐱=𝜃ℐ+𝜌²ℋ𝒦𝐱. (2.16) 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 𝐀𝐁^𝑇𝐮=1−𝐁0−𝑝𝜌²−1𝐇𝐒+𝐈𝟎𝜌²𝐮𝐁𝐒𝐈−𝑝. (2.17) We can obtain the following system of the equations: 𝐀𝐮+𝐁^𝑇𝜃𝑝=𝜃𝐮+𝜌²𝜃𝐇𝐒𝐮,−𝐁𝐮=−𝜌²𝐁𝐒𝐮+𝜃𝑝. (2.18)
We solve 𝑝 from the second equation 𝑝=(1/𝜃)𝐁((𝜃/𝜌²)𝐒−𝐈)𝐮. Plug in 𝑝 into the first equation, we have 1𝐀𝐮+𝜌²𝐁^𝑇𝜃𝐁𝐒𝐮−𝜃𝐮−𝜌²1𝐇𝐒𝐮=𝜃𝐁^𝑇𝐁𝐮. (2.19) Applying 𝜃𝐮^𝐻 to the both sides, we can obtain the following equation: 𝜃𝐮^𝐻𝜃𝐀𝐮+𝜌²𝐮^𝐻𝐁^𝑇𝐁𝐒𝐮−𝜃²𝐮^𝐻𝜃𝐮−²𝜌²𝐮^𝐻𝐇𝐒𝐮=𝐮^𝐻𝐁^𝑇𝐁𝐮=‖𝐁𝐮‖₂≥0. (2.20) Therefore, 𝜃𝐀+(𝜃/𝜌²)𝐁^𝑇𝐁𝐒−𝜃²𝐈−(𝜃²/𝜌²)𝐇𝐒 is a Hermitian matrix. Suppose that Re(𝜃)=0, that is, 𝜃=𝑡𝐢, where 𝑡≠0. We denote that 𝐆=𝜃𝐀+(𝜃/𝜌²)𝐁^𝑇𝐁𝐒−𝜃²𝐈−(𝜃²/𝜌²)𝐇𝐒=𝑡𝐢𝐀+(𝑡𝑖/𝜌²)𝐁^𝑇𝐁𝐒+𝑡²𝐈+(𝑡²/𝜌²)𝐇𝐒. While 𝐆^𝐻=−𝑡𝐢𝐀^𝐻+(𝑡𝑖/𝜌²)𝐁^𝑇𝐁𝐒+𝑡²𝐈+(𝑡²/𝜌²)𝐇𝐒=𝐆, which leads to 𝐀=𝐀^𝐻. A contradiction. Because 𝐀 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 preconditioner 𝑃_𝛼𝑧_𝑘=𝑟_𝑘, we first solve the systemℋ𝑤_𝑘=𝑟_𝑘, (3.1) for 𝑤_𝑘, followed by𝒦𝑧_𝑘=𝑤_𝑘. (3.2)

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 matrix 𝐾+𝜌𝐼_𝑛. This can be done by sparse LU factorization, preconditioned GMRES method. Notice that since 𝐾+𝜌𝐼 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 𝐯 appearing (via 𝐰=∇×𝐯) in the rotation form of the unsteady Oseen problem. Here we present results for the choice 𝑤=16𝑥(𝑥−1)+16𝑦(𝑦−1) (2D case) and 𝐰=(−4𝑧(1−2𝑥),4𝑧(2𝑦−1),2𝑦(1−𝑦)) (3D case). The 2D case corresponds to the choice of 𝐯. The equations were discretized with a Marker-and-Cell (MAC) scheme with a uniform mesh size ℎ. The outer iteration (full GMRES) was stopped when 𝑟_𝑘2𝑟₀₂<10⁻⁶, (4.1) where 𝑟_𝑘 denotes the residual vector at step 𝑘. For the results presented in this section, we use the exact solver to solve the two systems.

In Tables 1, 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 ℎ-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.

Table 2 

Results for 2D unsteady Oseen problem different values of 𝛼 (exact solves) and 𝜈=0.05.

Table 3 

Results for 2D unsteady Oseen problem different values of 𝛼 (exact solves) and 𝜈=0.01.

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.

Acknowledgment

The authors would like to thank Michele Benzi for helpful suggestions.