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

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.


Introduction
We study the numerical solution methods of the incompressible viscous fluid problems with the following form: Bu g on ∂Ω × 0, Γ , p p x, t ∈ R 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 Ω p dx 0. 1.5 The source function f is given on Ω.Here ν > 0 is a given constant called the kinematic viscosity, which is ν O Re −1 .Re is the Reynolds number: Re V L/ν, 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 B 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 Bu g on ∂Ω. 1.8 Here v is a known divergence-free vector obtained from the previous linearized step e.g., v u k .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 Au b has the form In this paper, we are interested in an alternative linearization of the steady-state Navier-Stokes equation.Based on the identity 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 After substituting the right-hand side into 1.6 , we find that the corresponding linearized equations have the following form: Bu g on ∂Ω × 0, T , 1.14 where P p 1/2 u 2 2 is the so-called Bernoulli pressure.For the two-dimensional case where w ∇ × v −∂v 1 /∂x 2 ∂v 2 /∂x 1 is a scalar function.
In the three-dimensional case, we have here w 1 , w 2 , w 3 w ∇ × v, where w i denotes the ith component of ∇ × v. Assume v v 1 , v 2 , v 3 , then we have the formal expression of w Here the divergence-free vector field v 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 Ax b, where Here A 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, The rectangular matrix B is the discretization of the negative divergence, and B T 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 For some discretization methods, a stabilization matrix needs to be added to the 2,2 block of A, 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 A with or without stabilization term in the 2D case.Such a stabilization is not necessary for the MAC discretization.To solve the system Ax 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.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.

The Alternative HSS Preconditioner
The alternative HSS preconditioner is based on the nonsymmetric formulation 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 ∈ C .The spectra of the nonsymmetric formulation is more friendly to the convergence of Krylov subspace iterations.
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 A T and K ≡ 1/2 A − A T , we have the following splitting of A into two parts: We denote Therefore, we defined the preconditioner as the following:

2.4
Here I n 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 A: Here I denotes the identity matrix of order n m.Note that is the shifted discretized Stokes problem, where I n denotes the identity matrix of order n, and I m denotes the identity matrix of order m.We obtian that is nonsingular and has positive definite symmetric part.
Alternating between these two splittings leads to the the following iteration: 2.8 k 0, 1, . . . .Here b denotes the right-hand side of 1.9 ; the initial guess u 0 is chosen arbitrarily.Elimination of u k 1/2 from 2.8 leads to a stationary fixed-point iteration of the form where

2.10
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, T ρ < 1 for all ρ > 0 and α ≥ 0.
Proof.Consider the splitting 2.3 .The iteration matrix T ρ is similar to
Assume that λ is one eigenvalue of the preconditioned linear system P −1 Ax P −1 b.We have
We define that θ λρ/ 2 − λ .Since λ / 2, θ is welldefined.Thus, with λ 2θ/ θ 2 , we consider the following equation: 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

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 r k 2 r 0 2 < 10 −6 , 4.1 where r k 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 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 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.

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.
1 to 1.4 are also known as the Navier-Stokes equations.Here Ω is an open set of R d , where d 2, or d 3, with boundary ∂Ω; the variable u u x, t ∈ R d is a vector-valued function representing the velocity of the fluid, and the scalar function Journal of Applied Mathematics

Figure 1 :
Figure 1: Sparsity patterns for different types of the Oseen problem in rotation form.
and D 3 are all diagonal matrices or weighted mass matrices.Typical sparsity patterns for A in the 2D and 3D case are displayed in Figures 1 a and 1 b .

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

Theorem 2 . 1 .
− K is the iteration matrix and c : S ρI −1 ρI − H H ρI −1 ρI − S .The iteration converges for arbitrary initial guesses u 0 and right-hand sides b to the solution u * A −1 if and only if T ρ < 1, where T ρ denotes the spectral radius of T ρ .Consider the problem 2.1 , that is,

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.