Numerical Solution of the Navier–Stokes Equations Using Multigrid Methods with HSS-Based and STS-Based Smoothers

Multigrid methods (MGMs) are used for discretized systems of partial differential equations (PDEs) which arise from finite difference approximation of the incompressible Navier–Stokes equations. After discretization and linearization of the equations, systems of linear algebraic equations (SLAEs) with a strongly non-Hermitian matrix appear. Hermitian/skew-Hermitian splitting (HSS) and skew-Hermitian triangular splitting (STS) methods are considered as smoothers in the MGM for solving the SLAE. Numerical results for an algebraic multigrid (AMG) method with HSS-based smoothers are presented.


Introduction
Mathematical modeling of hydrodynamics is the base for research of various natural phenomena, technological processes, and environmental problems. The main equations describing this problem are the Navier-Stokes equations. Development and research of effective numerical algorithms for solving these equations and their practical realization is an actual task. The use of the MGM for the numerical solution of the Navier-Stokes equations describing the motion of an incompressible viscous fluid is discussed. Currently, various discretization methods for the corresponding differential model are known. However, with any choice of the discretizing method, the problem of constructing effective methods for solving large systems of algebraic equations-to which the discrete model is reduced-arises. This problem is especially relevant in the nonstationary case, when multiple solutions of the systems of algebraic equations are required at each discrete time step.
To discretize the system of two-dimensional Navier-Stokes equations on regular grids, we use the finite difference method. The equations are considered in the natural variables "velocity-pressure": where P/ρ is replaced by P (i.e., ρ is normalized at 1), P is the static pressure, V is the velocity vector, and ν is the kinematic viscosity coefficient. At the initial moment of time and at the boundary of the domain, the initial and boundary conditions are set, respectively. Flow simulating is accompanied by a number of mathematical difficulties. One of the problems in solving this system is the nonlinearity associated with convective terms in the equations, which can lead to the appearance of oscillations of the solution in regions with large gradients. The main efforts of the researchers were directed at overcoming the difficulties associated with the nonlinearity of the Navier-Stokes system of equations.
One of the most time-consuming stages of the computational procedure is finding the solution of the system of linear algebraic equations (SLAE). Modern application packages usually use the linearization of the original equations, and Krylov subspace methods are used to solve the resulting SLAEs. Despite the fact that these methods have proven themselves well, they have some problems in cases of significant nonsymmetry of the SLAEs-associated, for example, with variable coefficients in differential equations or using complex numerical boundary conditions. For time discretization of the unsteady problem, we use an implicit difference scheme. Here, we do not specifically consider the stages of discretization and linearization of the Navier-Stokes equations, but focus on solving SLAEs. Given that the SLAEs resulting from the use of the implicit time schemes have a large dimension and a sparse nonsymmetric matrix, we propose using the MGM to solve them.
Thus, we consider the iterative solution of the large sparse SLAE where A ∈ C n×n is a non-Hermitian and positive definite matrix. Naturally, the matrix A can be split as where and A * denotes the conjugate transpose of the matrix A. Positive definiteness of the matrix A means that for all x ∈ C n \ {0}, x * A 0 x > 0. Here, x * denotes the conjugate transpose of the complex vector x. Let in some matrix norm ||| · |||, |||A 0 ||| << |||A 1 |||, then the matrix A is called a strongly non-Hermitian one. This situation occurs in many real applications, such as the discretization of the Navier-Stokes equations. The Hermitian and skew-Hermitian splitting (HSS) iteration methods, based on HS splitting (3) and (4), for solving large sparse non-Hermitian positive definite SLAE were firstly proposed in [1]. The HSS iteration method has been widely developed in [2][3][4][5] and others.
Then, we can split the skew-Hermitian part A 1 of the matrix A ∈ C n×n into where K L and K U are the strictly lower and the strictly upper triangular parts of A 1 , respectively. Obviously, that K L = −K * U . Based on the splitting (3)-(5) in [6][7][8] classes of skew-Hermitian triangular splitting (STS), iteration methods for solving SLAE (2) have been proposed. The triangular operator of the STS uses only the skew-Hermitian part of the coefficient matrix A. These methods have been further developed in [9][10][11][12].
The use of the multigrid method (MGM) with the STS-based smoothers for solving convection-diffusion problems has been studied in [13]. The convergence of the MGM with the STS-based smoothers has also been proved in this research. The local Fourier analysis of the MGM with the triangular skew-symmetric smoothers has been performed in [14]. The results of numerical experiments for convection-diffusion problems with large Peclet numbers by the geometric MGM have been presented in both researches.
In [15], it was shown that the MGMs with the HSS-based smoothers converge uniformly for second-order nonselfadjoint elliptic boundary value problems. This happens if the mesh size of the coarsest grid is sufficiently small, but independent of the number of the multigrid levels.

Multigrid Methods
The MGMs are proving themselves as very successful tools for solving the SLAE associated with discretization of partial differential equations (PDEs).
The main idea of the MGM has been proposed by R.P. Fedorenko in [16]. Then, A. Brandt [17], W. Hackbusch [18], and other researchers showed the efficiency of the multigrid approach and extended Fedorenko's idea.
The multigrid technique is based on two principles: error smoothing and coarse grid correction. The smoothing property is fundamental for the MGM. It is connected with fast damping high-frequency Fourier components of an initial error in decomposition on the basis from eigenvectors.
There exist two approaches in the MGM: geometric multigrid and algebraic multigrid methods. Geometric multigrid methods were critical to the early development of the MGM and still play an important role today. Nevertheless, there are classes of problems for which geometric techniques are too difficult to apply or cannot be used at all. These classes can be solved by the algebraic multigrid (AMG) methods, as introduced in [19,20].
The MGM is not a fixed algorithm. Rather, there is a multigrid technique that defines its scope. The efficiency of the MGM depends on the adjustment of its components to the considered problem [21].
The key to this is the correct choice of its components and effective interaction between smoothing and coarse-grid correction [22]. We need to use special iteration methods as smoothers for the MGM and nonstandard course-grid correction to a good approximation of the smooth error components.
The smoothing method is the central component of the multigrid algorithm; it is the most dependent part of the MGM on the type of the problem being solved. The role of smoothing methods is that they should not so much reduce the total error as smooth it (namely, suppress the high-frequency harmonics of the error) so that the error can be well approximated on a coarse grid.
Standard smoothing methods are linear iteration methods, for example, the Gauss-Seidel method. An alternative is the following methods: The MGMs can be used as solvers as well as preconditioners. The MGMs have been widely used for complicated nonsymmetric and nonlinear systems, like the Lame equations of elasticity or the Navier-Stokes problems.

Smoothers Based on the HSS and the STS Iteration Methods
A particular problem when using the MGM is the choice of smoothers. There are a number of iteration methods that can be used as smoothers, but not all of them are effective for solving strongly non-Hermitian SLAEs. The behavior of the HSS and the STS iteration methods is similar to the behavior of the Gauss-Seidel method, which quickly damps the high-frequency harmonics of the error, slowing down in the future. We give the formulas of these iteration methods.
The HSS iteration method [1]: Given an initial guess v (0) , for k = 0, 1, 2, ... until {v (k) } convergence, compute where α is a given positive constant and I is an identity matrix. Bai, Golub, and Ng [1] proved that the HSS iteration method converges unconditionally to the exact solution of the SLAE (2). Moreover, the upper bound of the contraction factor depends on the spectrum of A 0 but is independent of the spectrum of A 1 .
We can rewrite the HSS iteration method in the following form: and The STS iteration method [6,8]: Given an initial guess v (0) and two positive parameters ω and τ.
ω and τ are two acceleration parameters, and B(ω) is defined by For the STS method a convergence analysis, optimal choice of parameters and an accelerating procedure have presented in [8]. As it was mentioned above, smoothers in the MGMs should have a smoothing effect on the error of approximation. It was shown in [14] that the skew-Hermitian triangular iteration methods have such properties. Therefore, these methods can be used as smoothers in the MGMs.

Numerical Experiments
A wide class of CFD (Computational Fluid Dynamics) problems is associated with solving the equations of motion of a viscous incompressible fluid with a predominance of convective transfer. As a model, we consider the problem of internal single-phase chemically homogeneous flows, which are described by the unsteady Navier-Stokes equations in the domain Ω with a solid boundary Γ. At the initial stages of the development of CFD, preference was given to explicit methods that were used to solve stationary and nonstationary Navier-Stokes equations. Recently, increased attention has been paid to implicit methods. This is primarily due to the insufficient computational efficiency of explicit methods in solving the equations of motion of a viscous fluid using small difference grids. From the point of view of computational linear algebra, the matrices obtained at each time step when integrating unsteady equations using implicit schemes (after linearization) are nonselfadjoint and require special iterative methods for their effective solution. Therefore, in this research, we suggest using the AMG with special smoothers to solve such SLAEs.
So, we consider the model unsteady Navier-Stokes problem or ∂u ∂x where ν is the kinematic viscosity coefficient; Re = UL/ν is the Reynolds number, where U is a characteristic velocity of the flow and L is a characteristic length scale; V = (u(x, y, t), v(x, y, t)) is the velocity vector; P is the static pressure; the initial pressure distribution is given by a linear function. The initial conditions are taken to be zero. At the boundary, no-slip conditions are accepted. It means that at a solid boundary, the fluid will have zero velocity relative to the boundary. There are no mass forces in the formulation; motion is determined only by the boundary and initial conditions for the velocity field as well as the initial pressure distribution. For convenience, only square domain Ω = (0, 1) × (0, 1) will be considered. We assume that the fluid motion occurs in the time interval [0, T]. Therefore, the equations are considered in the domain Ω × (0, T) with the boundary Γ × [0, T]. The Navier-Stokes equations with the introduced boundary conditions have a solution determined up to an arbitrary constant for pressure, therefore, an agreement was adopted on the next normalization Ω P(x, y, t)dxdy = 0, ∀t. The most common approach to solving the Navier-Stokes equations in natural variables essentially uses the replacement of the difference continuity equation by the difference Poisson equation for pressure. Following this approach, first the difference equations are constructed that approximate the mass and momentum conservation equations and then, by algebraic transformations, the Poisson equation for determining the pressure is derived. This equation is used in the calculations instead of the continuity equation.
First, the equations of motion and continuity (6) and (7) are rewritten in schematic form [23]: where R contains all convective and diffusive forces, We fix the time step δt and introduce a discrete time grid t n = nδt, n ≥ 0 and denote the approximation to f (t n , x, y) as f (n) . Then, the fully implicit scheme will have the form The Poisson equation for pressure is obtained by taking the divergence from both sides of the Equation (14), taking into account the continuity Equation (16): But following [23], at this moment, the Poisson Equation (18) does not need to be created. Instead, we need to do a discretization. In addition, the continuity Equation (16) is first discretized before substituting the discrete version of (14). To approximate the problem in space, the finite difference method is used. Let the equations in discrete form be given by where D h and G h are the discrete div and ∇ operator, respectively. Then, V h , P h and R h are the discrete grid functions corresponding with V, P, and R. After discretization of (16), the number of velocity unknowns equals the number of discrete momentum equations. The number of pressure unknowns is equal to the number of discrete continuity equations, since both are equal to the number of grid cells [23]. Our approach uses the idea of [23], but it differs in implementation. The uniform grid Ω is introduced in the domain Ω with steps h 1 and h 2 ; h 1 = 1/N 1 , h 2 = 1/N 2 , where N 1 , N 2 are the number of cells in each direction. The grid cells are positioned such that the cell faces coincide with the boundary Γ of Ω. The discretization in space of the Navier-Stokes equations is performed on MAC (Marker and Cell) [24] (staggered) grids when pressure P and velocities in two-dimensional problems are determined on three grids shifted relative to each other. So, P is located in the center of each cell, the x-component velocity u is on the middle points of vertical faces, the y-component velocity v is on the middle points of horizontal faces. For the MAC-method, the solution advanced in time by solving the momentum equation with the best current estimate of pressure distribution. Such a solution initially would not satisfy the continuity equation unless the correct pressure distribution was used. The pressure is improved by numerically solving the Poisson equation with estimated velocity field. We rewrite the equation for pressure in the following form: We now introduce the grid sets and the corresponding spaces: h be the linear space of vector functions defined on D 1 × D 2 and vanishing at the corresponding grid boundaries, and P h is the space of functions defined on D 3 and orthogonal to unity. Thus, Variables are denoted by a single set of indices, despite the fact that different variables are calculated at different grid nodes. As a result, the indices i, j refer to a set of three mismatched points.
The term R (n+1) in (15) contains the nonlinear terms. So, for treating this nonlinearity, Newton linearization around the old time level is used. For example, we want to linearize a nonlinear term The expression in the right-hand side of (23)  be set equal to zero. That is, the correction of pressure is required to compensate for nonzero dilation at the n iterative level. The Poisson equation is then solved for the revised pressure field. The improved pressure is then used in the momentum equation for better solution at time step. If the dilation (divergence of velocity field) is not zero, the cyclic process of solving the momentum equation and Poisson equation is repeated until the velocity field is divergence free.
Thus, our computational scheme can be represented as follows: 1. Velocity field components u = u (n+1) and v = v (n+1) are determined by solving the implicit momentum equation with P , and for treating nonlinearity, the Newton linearization around the old time level is used.
2. The Poisson equation with estimated velocity field components u = u (n+1) and v = v (n+1) is solved for the revised pressure field P = P (n+1) .
Revised velocity field components u and v are determined by solving the implicit momentum equation with revised pressure P. Process of solving the momentum equation and Poisson equation is repeated until the velocity field is divergence free. Thus, at each time step in solving the Navier-Stokes equation, we need to solve SLAE with nonsymmetric matrices that are solved by the AMG method with HSS smoothers.
There are two coarsening approaches in the AMG: RS and PMIS algorithms. Coarsening splits initial grid on C-points and F-points-coarse and fine grid points, respectively. The RS (Ruge-Stuben) algorithm [25] is a traditional coarsening approach. The RS algorithm is based on two heuristic criteria that achieve optimal convergence and minimal computational cost. The first criterion provides the achievement of good convergence, as the effective coarsening scheme should allow to accurately interpolate a smooth error. Then, it is desirable that each F-point (Fine-grid point) has as many strongly influencing C-points (Coarse-grid points) as possible [26]. The criterion, provided minimal computational cost for different levels of V-cycle, requires that the set of C-points is the maximum subset of all F-points, to obtain more accurate interpolation, provided that no C-point is strongly dependent on another C-point (the set is maximum and independent), since such points would have increased the computational costs without providing visible benefits of interpolation [26]. In general, as the convergence is increased, the computational costs of the V-cycle decrease. Therefore, the first criterion is strictly observed and the second one is guidance. The RS algorithm has two passes. The first pass splits the full grid in C and F points; the second one ensures strict implementation of the first criterion [26]. PMIS (parallel changes independent set), the algorithm of coarsening, is based on the same principles as the RS algorithm except that a heuristic criterion is not strictly observed, i.e., F-F connections without a common C-point are permitted. Unlike the RS coarsening, the PMIS is not sequential. However, the precision may be deteriorated because an insufficient number of points reduces the accuracy of interpolation [26].
Numerical experiments have been done using the PMIS-algorithm. In Tables 1 and 2 , we give the number of AMG-iterations with the HSS-based smoother on the different grids, where α is the parameter of the HSS iteration method. For comparison, we give the AMG calculations when the Gauss-Seidel method is used as the smoothing procedure. In our implementations, all iterations are started from the zero vector, and terminated when where r (p) = b − Av (p) is the residual vector of the SLAE (2) at the current iterate v (p) and r (0) is the initial residual. Our comparisons are done for the number of iteration steps and the elapsed CPU time (in seconds, in parentheses). The abbreviation "n.c." in Table 2  From the data shown in the Table 1, an increase in the number of iterations with an increase in the mesh size follows. However, this relates to some features of the algebraic approach in MGM (more precisely, the PMIS algorithm in the AMG). The traditional (a scalable) approach in the AMG (RS algorithm) works well for problems arising from the discretization of PDEs in two spatial dimensions. For many two-dimensional problems, a solver can be obtained with the number of iterations, regardless of the size of the problem n, as well as the solution time per iteration, linearly proportional to n. For the RS algorithm, the convergence factor is separated from unity and does not depend on the size of the problem n. But when using regular AMG interpolation in combination with PMIS, AMG convergence worsens depending on the size of the problem. This results in a loss of scalability [27]. However, when traditional AMG algorithms are applied to three-dimensional (3D) problems, numerical tests show [27] that in many cases scalability is lost. However, the number of iterations may remain constant. The computational complexity and size of the stencil can increase significantly, which will lead to an increase in execution time and memory usage. In addition, the PMIS algorithm allows for natural parallelization, unlike the RS algorithm. These properties of the PMIS algorithm seem promising to us for the further study of the three-dimensional Navier-Stokes equations using parallel computing.  Table 3 shows the number of iteration steps and CPU time of the AMG+HSS method depending on the value α, when ν = 10 −5 . For the AMG+HSS method, the optimal (experimental) parameter value that reduces the number of iterations depends on the size of the grid. As the grid size increases, the value of α, which provides the best convergence, decreases. Numerical experiments showed that for parameter values less than 0.2, the AMG+HSS method diverges. Thus, the numerical experiments have showed that the HSS-based smoothers can be effectively used for the AMG, in which the stage of coarse-grid correction can be considered as a kind of accelerating procedure of the HSS methods.

Conclusions
In our previous theoretical and numerical studies of the MGM with the STS-based smoothers, the stationary (and nonstationary) linear diffusion-convection equation with dominant convection was considered as a test problem [13,14]. All theoretical results and calculations were performed using geometric MGM. Here, we first use the HSS-method as the smoother in the algebraic MGM for solving the unsteady Navier-Stokes equations. It is supposed to further prove the theoretically smoothing properties of the HSS iteration methods and to prove the convergence of the MGM with the corresponding smoothers. In addition, theoretical and numerical results should be obtained for the MGM with the STS-based smoothers for the Navier-Stokes problem. The PMIS algorithm was not chosen by us by chance. Preliminary testing of it on this model problem showed its robustness. In addition, the PMIS algorithm allows for natural parallelization, unlike the RS algorithm. These properties of the PMIS algorithm seem promising to us for the further study of the three-dimensional Navier-Stokes equations using parallel computing.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.