**4. 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

$$\frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla)\mathbf{V} = -\nabla P + \nu \Delta \mathbf{V}\_{\prime} \tag{6}$$

$$
\operatorname{div} \mathbf{V} = 0,\tag{7}
$$

or

$$
\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \frac{\partial \boldsymbol{u}}{\partial \mathbf{x}} + \boldsymbol{v} \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{y}} + \frac{\partial P}{\partial \mathbf{x}} - \frac{1}{R\varepsilon} \left( \frac{\partial^2 \boldsymbol{u}}{\partial \mathbf{x}^2} + \frac{\partial^2 \boldsymbol{u}}{\partial \boldsymbol{y}^2} \right) = \mathbf{0},\tag{8}
$$

$$
\frac{
\partial \upsilon
}{
\partial t
} + u \frac{
\partial \upsilon
}{
\partial \mathbf{x}
} + v \frac{
\partial \upsilon
}{
\partial y
} + \frac{
\partial P
}{
\partial y
} - \frac{1}{Re} \left( \frac{
\partial^2 \upsilon
}{
\partial \mathbf{x}^2
} + \frac{
\partial^2 \upsilon
}{
\partial y^2
} \right) = 0,\tag{9}
$$

$$
\frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial y} = 0,
$$

$$\begin{aligned} u\left(\mathbf{x},y,t\right)&=0, & v\left(\mathbf{x},y,t\right)&=0 & on & \Gamma, \\ u\left(\mathbf{x},y,0\right)&=0, & v\left(\mathbf{x},y,0\right)&=0, \\ P\left(\mathbf{x},y,0\right)&=\mathbf{j}\mathbf{x} - \frac{1}{2}\mathbf{j}', & \mathbf{j}' &= const,\end{aligned}$$

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 <sup>Ω</sup> *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]:

$$\frac{\partial \mathbf{V}}{\partial t} + \nabla P = \mathbf{R}\_{\prime} \tag{11}$$

where **R** contains all convective and diffusive forces,

$$\mathbf{R} = -(\mathbf{V} \cdot \nabla)\mathbf{V} + \frac{1}{Re}\Delta \mathbf{V},\tag{12}$$

$$
\operatorname{div} \mathbf{V} = 0.\tag{13}
$$

We fix the time step *δt* and introduce a discrete time grid *t <sup>n</sup>* = *<sup>n</sup>δt*, *<sup>n</sup>* ≥ 0 and denote the approximation to *f*(*t <sup>n</sup>*, *x*, *y*) as *f* (*n*). Then, the fully implicit scheme will have the form

$$\frac{1}{\delta t}(\mathbf{V}^{(n+1)} - \mathbf{V}^{(n)}) + \nabla P^{(n+1)} = \mathbf{R}^{(n+1)},\tag{14}$$

*Symmetry* **2020**, *12*, 233

$$\mathbf{R}^{(n+1)} = -(\mathbf{V}^{(n+1)} \cdot \nabla)\mathbf{V}^{(n+1)} + \frac{1}{Re}\Delta\mathbf{V}^{(n+1)},\tag{15}$$

$$
\operatorname{div} \mathbf{V}^{(n+1)} = 0,\tag{16}
$$

$$|\mathbf{V}^{(n+1)}|\_{\Gamma} = 0.\tag{17}$$

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):

$$
\Delta P^{(n+1)} - \operatorname{div} \mathbf{R}^{(n+1)} = \operatorname{div} \frac{\mathbf{V}^{(n)}}{\delta t}.\tag{18}
$$

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

$$D\_h \mathbf{V}\_h^{(n+1)} = 0,\tag{19}$$

$$\frac{1}{\delta t}(\mathbf{V}\_h^{(n+1)} - \mathbf{V}\_h^{(n)}) + G\_h P\_h^{(n+1)} = \mathbf{R}\_h^{(n+1)},\tag{20}$$

$$\mathbf{V}\_h^{(n+1)}|\_{\Gamma} = 0,\tag{21}$$

where *Dh* and *Gh* are the discrete *div* and ∇ operator, respectively. Then, **V***h*, *Ph* and **R***<sup>h</sup>* 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*<sup>1</sup> and *h*2; *h*<sup>1</sup> = 1/*N*1, *h*<sup>2</sup> = 1/*N*2, where *N*1, *N*<sup>2</sup> 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:

$$
\Delta P = \frac{d}{dx}\left(-\left(u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}\right) + \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)\right) + \tag{22}
$$

$$
\frac{d}{dy}\left(-\left(u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y}\right) + \frac{1}{Re}\left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right)\right) - \frac{\partial}{\partial t}\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right).
$$

We now introduce the grid sets and the corresponding spaces:

+

$$\overline{D}\_1 = \{ \mathbf{x}\_{ij} = ((i + 1/2)h\_1, jh\_2) : i = 0, \ldots, N\_1 - 1, j = 0, \ldots, N\_2 \},$$

$$\overline{D}\_2 = \{ \mathbf{x}\_{ij} = (ih\_1 \iota \,(j + 1/2)h\_2) : i = 0, \ldots, N\_1, j = 0, \ldots, N\_2 - 1 \},$$

$$D\_3 = \{ \mathbf{x}\_{ij} = (ih\_1 \iota \, jh\_2) : i = 1, \ldots, N\_1 - 1, j = 1, \ldots, N\_2 - 1 \}.$$

Let **V***<sup>h</sup>* = *V*1,*<sup>h</sup>* × *V*2,*<sup>h</sup>* be the linear space of vector functions defined on *D*<sup>1</sup> × *D*<sup>2</sup> and vanishing at the corresponding grid boundaries, and *Ph* is the space of functions defined on *D*<sup>3</sup> and orthogonal to unity. Thus,

$$V\_{1,l} = \{ u\_{ij} = u(\mathbf{x}\_{ij}) : \mathbf{x}\_{ij} \in \overline{D}\_1, \mathbf{u}\_{0,j} = \mathbf{u}\_{N\_1 - 1, j} = u\_{i,0} = \boldsymbol{u}\_{i, N\_2} = 0 \},$$

$$V\_{2,l} = \{ v\_{ij} = v(\mathbf{x}\_{ij}) : \mathbf{x}\_{ij} \in \overline{D}\_2, v\_{0,j} = v\_{N\_1, j} = v\_{i,0} = v\_{i, N\_2 - 1} = 0 \},$$

$$P\_{\mathcal{h}} = \{ P\_{\mathcal{ij}} = P(\mathbf{x}\_{ij}) : \mathbf{x}\_{ij} \in D\_3, \sum\_{ij} h\_1 h\_2 P\_{\mathcal{ij}} = 0 \}.$$

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 *<sup>u</sup>*(*n*+1)*φ*(*n*+1) *<sup>x</sup>* , then

$$
\mu^{(n+1)}\phi\_x^{(n+1)} = \mu^{(n)}\phi\_x^{(n+1)} + \mu^{(n+1)}\phi\_x^{(n)} - \mu^{(n)}\phi\_x^{(n)} + \mathcal{O}(\delta t^2).\tag{23}
$$

The expression in the right-hand side of (23) is linear in the variables at the new time level and possesses a discretization error *O*(*δt* 2).

Let *<sup>D</sup>* <sup>=</sup> *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂v <sup>∂</sup><sup>y</sup>* in (22) be the local dilation term, and other terms with velocity field determined

from the solution of momentum equation with a provisional estimate of pressure *P* = ( *f* 1 *f* 2 ), counter, and *D*(*n*+1) *ij* 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.

$$\frac{u\_{ij}^{(n+1)} - u\_{ij}^{(n)}}{\delta t} + \left( u\_{ij}^{(n)} \left( \frac{u\_{ij}^{(n+1)} - u\_{i-1,j}^{(n+1)}}{h\_1} \right) + u\_{ij}^{(n+1)} \left( \frac{u\_{ij}^{(n)} - u\_{i-1,j}^{(n)}}{h\_1} \right) - u\_{ij}^{(n)} \left( \frac{u\_{ij}^{(n)} - u\_{i-1,j}^{(n)}}{h\_1} \right) \right) + \dots$$

$$+ \left( v\_{ij}^{(n)} \left( \frac{u\_{ij}^{(n+1)} - u\_{i,j-1}^{(n+1)}}{h\_2} \right) + v\_{ij}^{(n+1)} \left( \frac{u\_{ij}^{(n)} - u\_{i,j-1}^{(n)}}{h\_2} \right) - v\_{ij}^{(n)} \left( \frac{u\_{ij}^{(n)} - u\_{i,j-1}^{(n)}}{h\_2} \right) \right) - \dots \tag{24}$$

$$- \frac{1}{\text{Re}} \left( \frac{u\_{i+1,j}^{(n+1)} - 2u\_{ij}^{(n+1)} + u\_{i-1,j}^{(n+1)}}{h\_1^2} + \frac{u\_{i,j+1}^{(n+1)} - 2u\_{ij}^{(n+1)} + u\_{i,j-1}^{(n+1)}}{h\_2^2} \right) = f\_{1\prime}^{\prime}$$

$$\frac{\upsilon\_{ij}^{(n+1)} - \upsilon\_{ij}^{(n)}}{\delta t} + \left( u\_{ij}^{(n)} \left( \frac{\upsilon\_{ij}^{(n+1)} - \upsilon\_{i-1,j}^{(n+1)}}{h\_1} \right) + u\_{ij}^{(n+1)} \left( \frac{\upsilon\_{ij}^{(n)} - \upsilon\_{i-1,j}^{(n)}}{h\_1} \right) - u\_{ij}^{(n)} \left( \frac{\upsilon\_{ij}^{(n)} - \upsilon\_{i-1,j}^{(n)}}{h\_1} \right) \right) + \dots$$

$$+ \left( \upsilon\_{ij}^{(n)} \left( \frac{\upsilon\_{ij}^{(n+1)} - \upsilon\_{i,j-1}^{(n+1)}}{h\_2} \right) + \upsilon\_{ij}^{(n+1)} \left( \frac{\upsilon\_{ij}^{(n)} - \upsilon\_{i,j-1}^{(n)}}{h\_2} \right) - \upsilon\_{ij}^{(n)} \left( \frac{\upsilon\_{ij}^{(n)} - \upsilon\_{i,j-1}^{(n)}}{h\_2} \right) \right) - \tag{25}$$

$$-\frac{1}{\mathcal{R}c} \left( \frac{\upsilon\_{i+1,j}^{(n+1)} - 2\upsilon\_{ij}^{(n+1)} + \upsilon\_{i-1,j}^{(n+1)}}{h\_1^2} + \frac{\upsilon\_{i,j+1}^{(n+1)} - 2\upsilon\_{ij}^{(n+1)} + \upsilon\_{i,j-1}^{(n+1)}}{h\_2^2} \right) = f\_{2,i}'$$

$$\frac{u\_{ij}^{(n+1)} - u\_{i-1,j}^{(n+1)}}{h\_1} + \frac{\upsilon\_{ij}^{(n+1)} - \upsilon\_{i,j-1}^{(n+1)}}{h\_2} = 0. \tag{26}$$

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).

1 *h*2 1 *P*(*n*+1) *<sup>i</sup>*+1,*<sup>j</sup>* <sup>−</sup> <sup>2</sup>*P*(*n*+1) *<sup>i</sup>*,*<sup>j</sup>* <sup>+</sup> *<sup>P</sup>*(*n*+1) *i*−1,*j* + 1 *h*2 2 *P*(*n*+1) *<sup>i</sup>*,*j*+<sup>1</sup> <sup>−</sup> <sup>2</sup>*P*(*n*+1) *<sup>i</sup>*,*<sup>j</sup>* <sup>+</sup> *<sup>P</sup>*(*n*+1) *i*,*j*−1 = <sup>=</sup> <sup>1</sup> *δt* ⎛ ⎝*u*(*n*+1) *ij* <sup>−</sup> *<sup>u</sup>*(*n*+1) *i*−1,*j h*1 + *v* (*n*+1) *ij* − *v* (*n*+1) *i*,*j*−1 *h*2 ⎞ <sup>⎠</sup> <sup>+</sup> + 1 *h*1 ⎛ ⎝−*u*(*n*+1) *ij* ⎛ ⎝*v* (*n*+1) *ij* − *v* (*n*+1) *i*−1,*j h*1 ⎞ <sup>⎠</sup> − *<sup>v</sup>* (*n*+1) *ij* ⎛ ⎝*v* (*n*+1) *ij* − *v* (*n*+1) *i*,*j*−1 *h*2 ⎞ ⎠ ⎞ <sup>⎠</sup> <sup>+</sup> + 1 *Reh*<sup>1</sup> ⎛ ⎝*u*(*n*+1) *<sup>i</sup>*+1,*<sup>j</sup>* <sup>−</sup> <sup>2</sup>*u*(*n*+1) *ij* <sup>+</sup> *<sup>u</sup>*(*n*+1) *i*−1,*j h*2 1 <sup>+</sup> *<sup>u</sup>*(*n*+1) *<sup>i</sup>*,*j*+<sup>1</sup> <sup>−</sup> <sup>2</sup>*u*(*n*+1) *ij* <sup>+</sup> *<sup>u</sup>*(*n*+1) *i*,*j*−1 *h*2 2 ⎞ <sup>⎠</sup> <sup>+</sup> + 1 *h*2 ⎛ ⎝−*u*(*n*+1) *ij* ⎛ ⎝*v* (*n*+1) *ij* − *v* (*n*+1) *i*−1,*j h*1 ⎞ <sup>⎠</sup> − *<sup>v</sup>* (*n*+1) *ij* ⎛ ⎝*v* (*n*+1) *ij* − *v* (*n*+1) *i*,*j*−1 *h*2 ⎞ ⎠ ⎞ <sup>⎠</sup> <sup>+</sup> + 1 *Reh*<sup>2</sup> ⎛ ⎝*v* (*n*+1) *<sup>i</sup>*+1,*<sup>j</sup>* − 2*v* (*n*+1) *ij* + *v* (*n*+1) *i*−1,*j h*2 1 + *v* (*n*+1) *<sup>i</sup>*,*j*+<sup>1</sup> − 2*v* (*n*+1) *ij* + *v* (*n*+1) *i*,*j*−1 *h*2 2 ⎞ ⎠ .

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

$$\frac{||r^{(p)}||\_2}{||r^{(0)}||\_2} \le 10^{-6} \text{ .}$$

where *<sup>r</sup>*(*p*) <sup>=</sup> *<sup>b</sup>* <sup>−</sup> *Av*(*p*) is the residual vector of the SLAE (2) at the current iterate *<sup>v</sup>*(*p*) and *<sup>r</sup>*(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 means "no convergence".

The experiments are run in MATLAB (version R2018b) with a machine precision of 10<sup>−</sup>16.

From Tables 1 and 2, it follows that the AMG methods with the HSS-smoother have fast convergence speed for all tested values of the viscosity coefficient (*<sup>ν</sup>* = <sup>10</sup>−<sup>1</sup> ÷ <sup>10</sup>−5) on all used grids, while the AMG with the Gauss–Seidel smoother does not converge for *ν* = 10−4, 10−<sup>5</sup> on all considered grids, and does not converge on the grids 260 × 260 and 520 × 520 nodes for all values of the viscosity coefficient. For all tests, the AMG+HSS (Algebraic multigrid method with Hermitian/ Skew-Hermitian Splitting smoother) outperforms the AMG+GS (Algebraic multigrid method with Gauss–Seidel smoother) with respect to both number of iteration steps and CPU time. Moreover, the number of iteration steps and CPU time increase with increasing grid size for both methods.

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 1.** Algebraic multigrid (AMG)+(HSS) Hermitian/skew-Hermitian splitting iterations with different *ν*.



**Table 2.** Algebraic multigrid (AMG)+(GS) Gauss–Seidel iterations with different *ν*.

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.

**Table 3.** (AMG)+(HSS) iterations with different *α*, *ν* = 10<sup>−</sup>5.


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.
