Algebraic multigrid method with skew-Hermitian smoothers Tatiana S. Martynovaa , Galina V. Muratovaa and Zeng-Qi Wangb a Vorovich Institute of Mathematics, Mechanics and Computer Science, Southern Federal University, 200/1 Stachki Ave., Bld. 2, Rostov-on-Don, 344090, Russian Federation b School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai 200240, P.R. China Abstract An algebraic multigrid method with a special type of smoothers is used for solving large sparse systems of linear algebraic equations with a strongly non-Hermitian matrix. Namely, skew-Hermitian triangular and product- triangular iterative methods have been used as the smoothers. When implementing the algebraic multigrid method, the Parallel Modified Independent Set (PMIS) algorithm for coarsening grid has been used. The two- dimensional unsteady Navier-Stokes problem for viscous incompressible fluid, written in primitive variables "velocity-pressure", is considered as a model one. Numerical experiments have shown high efficiency of the algebraic multigrid method with the skew-Hermitian smoothers for this problem. Keywords algebraic multigrid method, skew-Hermitian smoothers, unsteady Navier-Stokes problem 1. Introduction Multigrid methods (MGMs) are called scalable or optimal because they can solve a linear system with 𝑁 unknowns with only 𝑂(𝑁 ) work. Because this work can be effectively distributed across a parallel machine, MGMs can solve ever larger problems on proportionally larger parallel computers in essentially constant time, making them an ideal solver for large-scale scientific simulation [1]. Nowadays MGMs are the most efficient tools for solving a large class of problems including bound- ary value problems for partial differential equations (PDEs) or systems of PDEs. The MGM efficiency depends on the adjustment of its components to the problem to be solved. Important parts of the MGM are a smoothing procedure and the coarse-grid correction. We suggest to use special iterative meth- ods as smoothers for the MGM and non-standard course-grid correction to a good approximation of the smooth error components. A successful use of MGMs for solving hard practical problems depends also on their robustness against bad or small parameters (e.g. geometrical parameters, singularities, anisotropies etc.) arising in many applications. It is known that some iteration methods have a smoothing effect on the error of approximation [2, 3]. This property is fundamental for the multigrid idea and is connected with fast damping high- frequency Fourier components of an initial error in decomposition on the basis from eigenvectors. 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 Far Eastern Workshop on Computational Technologies and Intelligent Systems, March 2–3, 2021, Khabarovsk, Russia " martynova@sfedu.ru (T.S. Martynova); muratova@sfedu.ru (G.V. Muratova); wangzengqi@sjtu.edu.cn (Z. Wang) ~ https://sfedu.ru/en/person/martynova (T.S. Martynova); https://sfedu.ru/en/person/muratova (G.V. Muratova); http://math.sjtu.edu.cn/faculty/zqwang1/zqwang.html (Z. Wang)  0000-0002-7073-1162 (T.S. Martynova); 0000-0002-9816-6759 (G.V. Muratova); 0000-0001-6087-4023 (Z. Wang) Β© 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) 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 [2, 3]. There exist two approaches in the MGM: geometric multigrid and algebraic multigrid methods (AMG) [1, 3]. AMG is one of the most important issues from a practical point of view. In contrast to the geometrical multigrid method, its algebraic version recovers algebraically the ingredients for a complete multigrid algorithm from the fine grid data only. AMG determines coarse grids, inter-grid transfer operators, and coarse-grid equations based solely on the matrix entries. There are two coarsening approaches in the AMG: RS (Ruge-Stuben) and PMIS (parallel changes independent set) algorithms. The RS algorithm [4] is a traditional coarsening approach, it is based on two heuristic criteria that achieve optimal convergence and minimal computational cost. The PMIS is based on the same principles as the RS algorithm except that a heuristic criterion is not strictly observed [5]. An advantage of the PMIS is the possibility of natural parallelization, which makes it applicable for 3-D problems, as well as for problems with a huge number of unknowns. We are interested in numerical solution of the SLAEs that arise from finite difference approximation of the 2-D incompressible unsteady Navier-Stokes equations governing the flow of viscous Newtonian fluids. We use two iterative methods based on the Hermitian and skew-Hermitian splitting of the ma- trix 𝐴 of the SLAE as smoothers for the MGM. These are the skew-Hermitian triangular splitting (STS) [6, 7] and the product-type skew-Hermitian triangular splitting (PSTS) [8] iteration methods. Based on the smoothing and approximation properties convergence of the MGM with the PSTS smoothers has been proved. Numerical experiments were carried out using the AMG (PMIS algorithm) with STS-, PSTS- and Gauss-Seidel based smoothers. 2. Smoothers based on the STS and the PSTS iteration methods Thus, we consider iterative solution of the large sparse SLAE 𝐴𝑒 = 𝑏, 𝑒, 𝑏 ∈ ℂ𝑛 , (1) where 𝐴 ∈ ℂ𝑛×𝑛 is a non-Hermitian and positive definite matrix. Naturally, the matrix 𝐴 can be split into its Hermitian and skew-Hermitian parts as 𝐴 = 𝐴0 + 𝐴1 , (2) where 1 1 𝐴0 = (𝐴 + π΄βˆ— ), 𝐴1 = (𝐴 βˆ’ π΄βˆ— ), (3) 2 2 and π΄βˆ— denotes the conjugate transpose of the matrix 𝐴. Positive definiteness of the matrix 𝐴 means that for all π‘₯ ∈ ℂ𝑛 ⧡ {0}, π‘₯ βˆ— 𝐴0 π‘₯ > 0. Here π‘₯ βˆ— denotes the conjugate transpose of the complex vector π‘₯. When 𝐴 is positive definite then its Hermitian part 𝐴0 is Hermitian positive definite, and the diagonal of its skew-Hermitian part 𝐴1 is purely imaginary. Let in some matrix norm ||| β‹… |||, |||𝐴0 ||| << |||𝐴1 |||, then the matrix 𝐴 is called strongly non-Hermitian one. This situation occurs in many real applications, such as the discretization of the Navier-Stokes equations. Notice, that π‘‘π‘–π‘Žπ‘”(𝐴1 ) = 0 when the coefficient matrix 𝐴 of the SLAE (1) is real. Then, we can split the skew-Hermitian part 𝐴1 of the matrix 𝐴 ∈ ℂ𝑛×𝑛 into 𝐴1 = 𝐾𝐿 + πΎπ‘ˆ , (4) where 𝐾𝐿 and πΎπ‘ˆ are the strictly lower and the strictly upper triangular parts of 𝐴1 , respectively. Obviously, that 𝐾𝐿 = βˆ’πΎπ‘ˆβˆ— . Based on the splitting (2)-(4) in [6, 7] classes of the STS iteration methods for solving SLAE (1) have been proposed. The triangular operator of the STS uses only skew-Hermitian part of the coefficient matrix 𝐴. These methods have been further developed in [9]. Moreover, based on the methodology given in [9], the PSTS iteration method has been investigated in [8]. 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 STS and the PSTS 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 STS iteration method [6, 7]: Given an initial guess 𝑒 (0) , and two positive parameters πœ” and 𝜏 . For π‘˜ = 0, 1, 2, ... until {𝑒 (π‘˜) } convergence, compute 𝑒 (π‘˜+1) = 𝐺(πœ”, 𝜏 )𝑒 (π‘˜) + 𝜏 𝐡(πœ”)βˆ’1 𝑏, where 𝐺(πœ”, 𝜏 ) = 𝐡(πœ”)βˆ’1 (𝐡(πœ”) βˆ’ 𝜏 𝐴), πœ” and 𝜏 are two acceleration parameters, and 𝐡(πœ”) is defined by 𝐡(πœ”) = 𝐡𝑐 + πœ”((1 + 𝑗)𝐾𝐿 + (1 βˆ’ 𝑗)πΎπ‘ˆ ), 𝑗 = Β±1. with 𝐡𝑐 ∈ ℂ𝑛×𝑛 a prescribed Hermitian matrix. The PSTS iteration method [8]: Given an initial guess 𝑒 (0) , and two positive parameters πœ” and 𝜏 . For π‘˜ = 0, 1, 2, ... until {𝑒 (π‘˜) } convergence, compute 𝑒 (π‘˜+1) = 𝐺(πœ”, 𝜏 )𝑒 (π‘˜) + 𝜏 𝐡(πœ”)βˆ’1 𝑏, (5) where 𝐺(πœ”, 𝜏 ) = 𝐡(πœ”)βˆ’1 (𝐡(πœ”) βˆ’ 𝜏 𝐴), πœ” and 𝜏 are two acceleration parameters, and 𝐡(πœ”) is defined by πœ” Μ‚ βˆ’1 πœ” 𝐡(πœ”) = (𝐡𝑐 + 𝐾𝐿 )𝐡𝑐 (𝐡𝑐 + πΎΜ‚π‘ˆ ), (6) 2 2 where 𝐾̂𝐿 = 𝐾𝐿 + 𝐹0 , πΎΜ‚π‘ˆ = πΎπ‘ˆ βˆ’ 𝐹0 , 𝐹0 ∈ ℂ𝑛×𝑛 is a Hermitian matrix, 𝐡𝑐 ∈ ℂ𝑛×𝑛 is a prescribed Hermitian positive definite matrix. Obviously that 𝐾̂𝐿 = βˆ’πΎΜ‚π‘ˆβˆ— , 𝐴1 = (𝐾𝐿 + 𝐹0 ) + (πΎπ‘ˆ βˆ’ 𝐹0 ) = 𝐾̂𝐿 + πΎΜ‚π‘ˆ . It was shown in [10] that the skew-Hermitian iteration methods have smoothing effect on the error of approximation. Therefore, these methods can be used as smoothers in the MGMs. Let us give a proof of the MGM convergence with the PSTS smoothers. 3. MGM convergence with the PSTS smoothers Let 𝐡0 (πœ”) and 𝐡1 (πœ”) be, respectively, the Hermitian and the skew-Hermitian parts of the matrix 𝐡(πœ”) defined by (6). That is to say, 𝐡(πœ”) = 𝐡0 (πœ”) + 𝐡1 (πœ”), (7) where { 𝐡0 (πœ”) = 21 (𝐡(πœ”) + 𝐡(πœ”)𝑇 ) = 𝐡𝑐 + ( πœ”2 )2 𝐾̂𝐿 π΅π‘βˆ’1 πΎΜ‚π‘ˆ , 𝐡1 (πœ”) = 12 (𝐡(πœ”) βˆ’ 𝐡(πœ”)𝑇 ) = πœ”2 (𝐾̂𝐿 + πΎΜ‚π‘ˆ ) = πœ”2 𝐴1 . Suppose that 0 < π›Όβ„Ž 𝐼 ≀ 𝐴0 ≀ π›½β„Ž 𝐼 , 0 < 𝛼𝑐 𝐼 ≀ 𝐡𝑐 ≀ 𝛽𝑐 𝐼 , (8) 𝛼𝑙 𝐼 ≀ 𝐾̂𝐿 π΅π‘βˆ’1 πΎΜ‚π‘ˆ ≀ 𝛽𝑙 𝐼 , 𝛼𝑠 𝐼 ≀ 𝐡0 (πœ”) ≀ 𝛽𝑠 𝐼 , (9) where 𝐼 is an identity matrix. Since 𝐴0 and 𝐡𝑐 are Hermitian positive definite matrices, 𝐡0 (πœ”) and 𝐾̂𝐿 π΅π‘βˆ’1 πΎΜ‚π‘ˆ are Hermitian matrices, and 𝐴1 is a skew-Hermitian matrix, the bounds π›Όπœ“ and π›½πœ“ , πœ“ = β„Ž, 𝑐, 𝑙, 𝑠, can be easily expressed with respect to the smallest and the largest eigenvalues or singular values of the corresponding matrix, respectively [8]. We need to require positive definiteness of the matrix 𝐡(πœ”) from (7). Since the √matrix 𝐾𝐿 𝐡𝑐 πΎπ‘ˆ is Μ‚ βˆ’1 Μ‚ negative definite, then 𝐡0 (πœ”) > 0, if the parameter πœ” ∈ (0, πœ”π‘šπ‘Žπ‘₯ ), where πœ”π‘šπ‘Žπ‘₯ = 2 (βˆ’ 𝛼𝛼𝑐𝑙 ). Notice, that 𝛼𝑙 and 𝛽𝑙 satisfy the following inequalities: πœ” 2 πœ” 2 𝛼𝑠 β‰₯ 𝛼𝑐 + ( ) 𝛼𝑙 , 𝛽𝑠 ≀ 𝛽𝑐 + ( ) 𝛽𝑙 . (10) 2 2 Theorem 1. [8] Let the matrices 𝐴 and 𝐡(πœ”) be positive definite. Then the PSTS iteration method is convergent, i.e., the spectral radius 𝜌(𝐺(πœ”, 𝜏 )) of its iteration matrix 𝐺(πœ”, 𝜏 ) is less than 1, provided that the acceleration parameters πœ” and 𝜏 satisfy 0 < 𝜏 < πœ”, 0 < πœ” < πœ”π‘šπ‘Žπ‘₯ , and 2 π›½β„Ž 0<𝜏 < , Θ= . (11) Θ 𝛼𝑐 + ( πœ”2 )2 𝛼𝑙 Let us study the convergence of the MGMs with the PSTS smoothers. We use some denotations and theoretical results from [11, 12]. Let 𝐻1 βŠ‚ 𝐻2 βŠ‚ ... be a family of nested finite-dimensional linear spaces. The dimension of π»π‘š is π‘›π‘š and the inner product is denoted by (.,.)π‘š with βˆ₯ . βˆ₯π‘š corresponding norm in π»π‘š , π‘š = 1, 2... We are interested in solution of the problem (1) in π»π‘š π΄π‘š 𝑒 = 𝑏. Let this problem have a unique solution for some 𝑏 ∈ π»π‘š and Μƒπ‘š + π΄Μ‚π‘š , π΄π‘š = 𝐴 where 𝐴 Μƒπ‘š is a symmetric positive definite operator in π»π‘š . Let π‘„π‘š is the other symmetric positive definite operator in π»π‘š and the condition for the spectral radius of the operator πΊπ‘š = π‘„π‘š Μƒπ‘š is βˆ’1 𝐴 satisfied 𝜌 (πΊπ‘š ) = 1. Using the operators π΄Μƒπ‘š and π‘„π‘š we define the energy norm Μƒπ‘š 𝑒, 𝑒)1/2 , β€–π‘’β€–π΄Μƒπ‘š = (𝐴 βˆ€π‘’ ∈ π»π‘š , and the inner product (𝑒, 𝑣)π‘„π‘š = (π‘„π‘š 𝑒, 𝑣), βˆ€π‘’, 𝑣 ∈ π»π‘š . Then we may define the following norms on π»π‘š : βˆ₯ 𝑒 βˆ₯𝑠,π‘š =βˆ₯ 𝑒 βˆ₯πΊπ‘šπ‘  = (πΊπ‘šπ‘  𝑒, 𝑒)1/2 , 𝑠 = 0, 1, ... π‘„π‘š βˆ₯ 𝑒 βˆ₯1,π‘š = (π΄Μƒπ‘š 𝑒, 𝑒)1/2 =βˆ₯ 𝑒 βˆ₯ Μƒ , π΄π‘š βˆ₯ 𝑒 βˆ₯0,π‘š = (π‘„π‘š 𝑒, 𝑒)1/2 = (𝑒, 𝑒)1/2 π‘„π‘š . Define the subspace πΉπ‘š βŠ‚ π»π‘š , πΉπ‘š = {𝑒 ∈ π»π‘š ∢ (π΄π‘š 𝑒, 𝑣) = 0, βˆ€π‘£ ∈ π»π‘šβˆ’1 } Now we make three basic assumptions [11, 12]: Assumption 1. There exists 𝑣 ∈ π»π‘šβˆ’1 , 𝛾 , 0 < 𝛾 ≀ 1 and 𝛿 < ∞ such that for βˆ€ 𝑒 ∈ π»π‘š βˆ₯ 𝑒 βˆ’ 𝑣 βˆ₯21,π‘š ≀ 𝛿 𝛾 βˆ₯ 𝑒 βˆ₯21+𝛾 ,π‘š Assumption 2. There exists πœ‚π‘š , πœ‚π‘š β†’ 0 (π‘š β†’ ∞) such that ∣ (π΄Μ‚π‘š 𝑒, 𝑣)π‘š ∣< πœ‚π‘š βˆ₯ 𝑒 βˆ₯1,π‘š βˆ₯ 𝑣 βˆ₯1,π‘š for βˆ€ 𝑒 ∈ πΉπ‘š , βˆ€π‘£ ∈ π»π‘š , π‘š is a positive integer. Assumption 3. There exists πœ‡π‘š , πœ‡π‘š β†’ 0 (π‘š β†’ ∞) such that ∣ (π΄Μ‚π‘š 𝑒, 𝑣 ) βˆ£β‰€ πœ‡π‘š βˆ₯ 𝑒 βˆ₯1,π‘š βˆ₯ 𝑣 βˆ₯0,π‘š π‘š for βˆ€π‘’, 𝑣 ∈ π»π‘š , π‘š is a positive integer. It is assumed that πœ‚π‘š , πœ‡π‘š are sufficiently small. In [12] it was shown that assumptions 1-3 were satisfied for sufficiently wide class of boundary problems in two dimension-limited domains with different boundary conditions. Assumption 1 means a generalization of the usual approximation property known for the symmetric case. Assumptions 2 and 3 restrict the non-symmetric part π΄Μ‚π‘š of the operator π΄π‘š . Let us consider two-grid algorithm. We denote the exact solution of the problem (1) by 𝑒 βˆ— . Let 𝑒 0 is an initial guess, 𝑒 1 is a problem solution after smoothing and 𝑒 2 is a problem solution after MGM-iteration. We need to estimate the energy-norm of the contraction number defined by βˆ₯ 𝑒 2 βˆ’ 𝑒 βˆ— βˆ₯1 𝜎 = sup , 𝑒0 β‰  π‘’βˆ—. βˆ₯ 𝑒 0 βˆ’ 𝑒 βˆ— βˆ₯1 Denote the error by 𝑒 𝑖 = 𝑒 𝑖 βˆ’ 𝑒 βˆ— , 𝑖 = 0, 1, 2. Theorem 2. [12] Let the three basic assumptions for the two-grid method be satisfied. Furthermore, let the following smoothing proposition also be satisfied: there exists Ξ” (0 < Ξ” < ∞) and πœ— > 0 such that βˆ₯ 𝑒 1 βˆ₯21 +πœ— βˆ₯ 𝑒 1 βˆ₯22 ≀ (1 + πœ‡Ξ”) βˆ₯ 𝑒 0 βˆ₯21 (12) with πœ‡=πœ‡π‘˜ from Assumption 3. Then 𝜎 ≀ πœŽΜƒ for the two-grid contraction number where ⎧ βŽͺ 1/2 ⎫ βŽͺ βŽͺ βŽͺ πœ‰ 2 +πœ– 2 𝜍 2 +2πœ–πœ‚πœ‰ 𝜍 βŽͺ βŽͺ 1βˆ’πœ‚ (1 + πœ‡Ξ”) ∢ πœŽΜƒ ≑ πœŽΜƒ (πœ–) ≑ sup ⎨ ( 1+𝛿 βˆ’1 πœ— ( 1+πœ‚ )2/𝛾 πœ‰ 2/𝛾 ) ⎬ βŽͺ βŽͺ βŽͺ βŽͺ πœ‰ 2 + 𝜁 2 βˆ’ 2πœ‚πœ‰ 𝜁 ≀ 1, 𝜁 , πœ‰ β‰₯ 0 βŽͺ βŽͺ ⎩ ⎭ βˆ₯ 𝑒 1 + 𝑒 βˆ— βˆ₯1 βˆ₯ 𝑒 βˆ— βˆ₯1 πœ‰ = , 𝜍= , βˆ₯ 𝑒 1 βˆ₯1 βˆ₯ 𝑒 1 βˆ₯1 where constants πœ— and Ξ” depend on properties of smoothing system and constants 𝛾 , 𝛿, πœ‚, πœ‡ taken from assumptions 1, 2, 3 respectively, while πœ– is calculation accuracy. In [12] the proof of the same theorem for the multigrid method is given. For further considerations it is necessary to have one more theorem from [11]. Theorem 3. [11] Let the problem (1) be solved by iterative method (5) - (6) written in the form of 𝑒 (𝑛+1) = 𝑒 (𝑛) βˆ’ π΅Μƒβˆ’1 (𝐴𝑒 (𝑛) βˆ’ 𝑏 ) , Μƒ = 𝐴0 , 𝐴 and let the three basic assumptions be satisfied. If there exists the constant πœ— > 0 so that inequality is satisfied 𝐡̃ + π΅Μƒβˆ— βˆ’ 𝐴 Μƒ βˆ’ 𝐡̃)βˆ— 𝑄 βˆ’1 (𝐴 Μƒ β‰₯ πœ— (𝐴 Μƒ βˆ’ 𝐡̃) , (13) then the smoothing assumption (12) is satisfied with Ξ” = (1 + πœ— ) 𝛽 (πœ‡π›½ + 2) , (14) where Μƒ (π΅Μƒβˆ’1 )βˆ— 𝐴 𝛽 = (𝜌 (π΅Μƒβˆ’1 𝐴 Μƒ))1/2 . Now we can prove the convergence of the method. Theorem 4. For the method (5) - (6) there exists the constant πœ— > 0 so that inequality (13) is satisfied. Proof. In case of using method (5)-(6) 𝐡̃ = 1/𝜏 𝐡, 𝑄 = 𝐴0 . Consider inequality (13) and transform left and right parts of it. Using (7) we obtain 𝐡̃ + π΅Μƒβˆ— βˆ’ 𝐴 Μƒ = 2/𝜏 𝐡0 βˆ’ 𝐴0 , (15) Μƒ βˆ’ 𝐡̃)βˆ— 𝐴 (𝐴 Μƒβˆ’1 (𝐴 Μƒ βˆ’ 𝐡̃) = (𝐴0 βˆ’ 1 π΅βˆ— ) π΄βˆ’1 1 0 (𝐴0 βˆ’ 𝜏 𝐡 ) = 𝜏 (16) = 𝐴0 βˆ’ 𝜏1 π΅βˆ— βˆ’ 𝜏1 𝐡 + 𝜏12 π΅βˆ— π΄βˆ’1 2 1 βˆ— βˆ’1 0 𝐡 = 𝐴0 βˆ’ 𝜏 𝐡0 + 𝜏 2 𝐡 𝐴0 𝐡 Denote 2 𝑆= 𝐡0 βˆ’ 𝐴0 , (17) 𝜏 then from (8) - (11) it’s easy to obtain that π›½β„Ž 2 πœ” 2 < , 𝛼𝑐 + ( 2 ) 𝛼𝑙 𝜏 thus, 2/𝜏 𝐡0 > 𝐴0 for acceleration parameters πœ”, 𝜏 satisfying the conditions of the Theorem 1, and 𝑆 = 𝑆 βˆ— > 0. From (15) - (17) we obtain that for inequality (13): 1 βˆ— βˆ’1 𝑆 β‰₯ πœ— (βˆ’π‘† + 𝐡 𝐴0 𝐡 ) . 𝜏2 Multiplying the left and the right parts of this inequality on 𝑆 βˆ’1/2 we obtain 1 βˆ’1/2 βˆ— βˆ’1 βˆ’1/2 𝐼 β‰₯ πœ— (βˆ’πΌ + 𝑆 𝐡 𝐴0 𝐡𝑆 ) . (18) 𝜏2 Denote 𝐿 = 𝜏1 π΄βˆ’1/2 0 𝐡𝑆 βˆ’1/2 . Then (18) is transformed to 𝐼 β‰₯ πœ— (πΏβˆ— 𝐿 βˆ’ 𝐼 ) . If we take 1 πœ—= , (19) βˆ₯ πΏβˆ— 𝐿 βˆ’ 𝐼 βˆ₯ then inequality (13) is satisfied. The Theorem 4 is proved. Thus, if the conditions of the Theorem 4 are satisfied then two-grid method converges with the PSTS-smoothers (6), πœ— is calculated by (19) and Ξ” is calculated by (14) correspondingly. The results obtained for the two-grid method can be extended to the multigrid method. 4. Numerical experiments Let us consider the 2-D unsteady Navier-Stokes equations of viscous incompressible fluids. The prim- itive variables mathematical formulation of this problem is: find V and 𝑃 (up to constant) so that πœ•V βˆ’ πœˆΞ”V + (V β‹… βˆ‡)V + βˆ‡π‘ƒ = f π‘œπ‘› Ξ© Γ— (0, 𝑇 ], πœ•π‘‘ 𝑑𝑖𝑣V = 0 π‘œπ‘› Ξ© Γ— [0, 𝑇 ], V=g π‘œπ‘› πœ•Ξ© Γ— [0, 𝑇 ], V(x, 0) = V0 (x) π‘œπ‘› Ξ© where 𝜈 is the kinematic viscosity coefficient; Ξ” is the Laplacian, βˆ‡ is the gradient, 𝑑𝑖𝑣 is the diver- gence; V = (𝑀(π‘₯, 𝑦, 𝑑), 𝑣(π‘₯, 𝑦, 𝑑)) is the velocity vector, 𝑃 is the pressure, g is the velocity prescribed on the boundary πœ•Ξ© and V0 is the divergence-free initial velocity field. We assume that the fluid motion occurs in the time interval [0, 𝑇 ]. One of the main difficulties in using the AMG is the constructing smoothers that are robust over a wide range of the viscosity coefficient, especially for small values of 𝜈. Classical MGMs were de- veloped for elliptic PDEs. When applied to nonelliptic and singular perturbation problems such as high-Reynolds flows, performance seems to deteriorate significantly. As a model problem, we consider the standard lid-driven cavity problem [13] in the square domain Ξ© = (0, 1) Γ— (0, 1). Let us write this problem in scalar dimensionless form πœ•π‘€ πœ•π‘€ πœ•π‘€ πœ•π‘ƒ 1 πœ•2𝑀 πœ•2𝑀 +𝑀 +𝑣 + βˆ’ + = 0, (20) πœ•π‘‘ πœ•π‘₯ πœ•π‘¦ πœ•π‘₯ 𝑅𝑒 ( πœ•π‘₯ 2 πœ•π‘¦ 2 ) πœ•π‘£ πœ•π‘£ πœ•π‘£ πœ•π‘ƒ 1 πœ•2𝑣 πœ•2𝑣 +𝑀 +𝑣 + βˆ’ + = 0, (21) πœ•π‘‘ πœ•π‘₯ πœ•π‘¦ πœ•π‘¦ 𝑅𝑒 ( πœ•π‘₯ 2 πœ•π‘¦ 2 ) πœ•π‘€ πœ•π‘£ + = 0, (22) πœ•π‘₯ πœ•π‘¦ 𝑀 = 0, 𝑣 = 0, (π‘₯ = 0, π‘₯ = 1, 𝑦 = 0), 𝑀 = 1, 𝑣 = 0, (𝑦 = 1), (23) 𝑀 (π‘₯, 𝑦, 0) = 0, 𝑣 (π‘₯, 𝑦, 0) = 0. where 𝑅𝑒 is the Reynolds number (𝑅𝑒 = π‘ˆ 𝐿/𝜈, π‘ˆ is a characteristic velocity of the flow, 𝐿 is a char- acteristic length scale). Homogeneous Dirichlet boundary conditions are prescribed for all velocity components with the exception of a positive unit horizontal velocity along the top face. The initial condition assumes the fluid to be at rest. We consider this problem as a test for study the effect of the suggested smoothers for the MGM. The fully implicit scheme has been used for the time discretization. The use of implicit schemes re- moves restrictions on the time integration step, which is selected based on the required computational accuracy. We fix the time step 𝛿𝑑 and introduce a discrete time grid 𝑑 𝑛 = 𝑛𝛿𝑑, 𝑛 β‰₯ 0: 1 (𝑛+1) 1 (V βˆ’ V(𝑛) ) + βˆ‡π‘ƒ (𝑛+1) = βˆ’(V(𝑛+1) β‹… βˆ‡)V(𝑛+1) + Ξ”V(𝑛+1) , (24) 𝛿𝑑 𝑅𝑒 𝑑𝑖𝑣V(𝑛+1) = 0, (25) The boundary conditions are taken at 𝑑 = 𝑑 𝑛+1 : 𝑀 (𝑛+1) = 0, 𝑣 (𝑛+1) = 0, (π‘₯ = 0, π‘₯ = 1, 𝑦 = 0), 𝑀 (𝑛+1) = 1, 𝑣 (𝑛+1) = 0, (𝑦 = 1). The uniform grid Ξ© is introduced in the domain Ξ© with steps β„Ž1 and β„Ž2 ; β„Ž1 = 1/𝑁1 , β„Ž2 = 1/𝑁2 , where 𝑁1 , 𝑁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 spatial discretization of the Navier-Stokes equations is performed on MAC (Marker-and-Cell) [14] (staggered) grids when pressure and velocities in two- dimensional problems are determined on three grids shifted relative to each other. So, the pressure 𝑃 is located in the center of each cell, the π‘₯-component velocity 𝑀 is on the middle points of vertical faces, the 𝑦-component velocity 𝑣 is on the middle points of horizontal faces. The continuity equation (25) is discretized at cell centers using central difference schemes. Thus we introduce the grid sets and the corresponding spaces: 𝐷 1 = {π‘₯𝑖𝑗 = ((𝑖 + 1/2)β„Ž1 , π‘—β„Ž2 ) ∢ 𝑖 = 0, ..., 𝑁1 βˆ’ 1, 𝑗 = 0, ..., 𝑁2 }, 𝐷 2 = {π‘₯𝑖𝑗 = (π‘–β„Ž1 , (𝑗 + 1/2)β„Ž2 ) ∢ 𝑖 = 0, ..., 𝑁1 , 𝑗 = 0, ..., 𝑁2 βˆ’ 1}, 𝐷3 = {π‘₯𝑖𝑗 = (π‘–β„Ž1 , π‘—β„Ž2 ) ∢ 𝑖 = 1, ..., 𝑁1 βˆ’ 1, 𝑗 = 1, ..., 𝑁2 βˆ’ 1}. Let π•β„Ž = 𝑉1,β„Ž Γ— 𝑉2,β„Ž be the linear space of vector functions defined on 𝐷 1 Γ— 𝐷 2 and vanishing at the corresponding grid boundaries, and π‘ƒβ„Ž is the space of functions defined on 𝐷3 and orthogonal to unity. Thus, 𝑉1,β„Ž = {𝑀𝑖𝑗 = 𝑀(π‘₯𝑖𝑗 ) ∢ π‘₯𝑖𝑗 ∈ 𝐷 1 , 𝑀0,𝑗 = 𝑀𝑁1 βˆ’1,𝑗 = 𝑀𝑖,0 = 𝑀𝑖,𝑁2 = 0}, 𝑉2,β„Ž = {𝑣𝑖𝑗 = 𝑣(π‘₯𝑖𝑗 ) ∢ π‘₯𝑖𝑗 ∈ 𝐷 2 , 𝑣0,𝑗 = 𝑣𝑁1 ,𝑗 = 𝑣𝑖,0 = 𝑣𝑖,𝑁2 βˆ’1 = 0}, π‘ƒβ„Ž = {𝑃𝑖𝑗 = 𝑃(π‘₯𝑖𝑗 ) ∢ π‘₯𝑖𝑗 ∈ 𝐷3 , βˆ‘ β„Ž1 β„Ž2 𝑃𝑖𝑗 = 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 𝑖, 𝑗 refer to a set of three mismatched points. The equations contain the discrete nonlinear terms. For treating this non-linearity Newton lin- earization around the old time level is used. If we want to linearize a nonlinear term 𝑀 (𝑛+1) πœ™π‘₯(𝑛+1) , then 𝑀 (𝑛+1) πœ™π‘₯(𝑛+1) = 𝑀 (𝑛) πœ™π‘₯(𝑛+1) + 𝑀 (𝑛+1) πœ™π‘₯(𝑛) βˆ’ 𝑀 (𝑛) πœ™π‘₯(𝑛) + 𝑂(𝛿𝑑 2 ). (26) Table 1 AMG+STS Iterations and CPU with Different 𝜈. Grid 𝜈 = 10βˆ’1 𝜈 = 10βˆ’2 𝜈 = 10βˆ’3 𝜈 = 10βˆ’4 𝜈 = 10βˆ’5 60 Γ— 60 36 (38.18) 34 (36.42) 31 (32.20) 28 (30.51) 20 (22.11) 120 Γ— 120 47 (62.32) 40 (55.54) 37 (53.60) 33 (50.20) 22 (29.31) 180 Γ— 180 49 (87.61) 44 (84.52) 38 (79.62) 34 (76.85) 24 (48.41) 260 Γ— 260 54 (122.71) 47 (129.33) 43 (132.42) 35 (110.26) 27 (83.33) 520 Γ— 520 59 (251.51) 56 (248.58) 54 (235.20) 44 (215.64) 36 (151.19) Table 2 AMG+PSTS Iterations and CPU with Different 𝜈. Grid 𝜈 = 10βˆ’1 𝜈 = 10βˆ’2 𝜈 = 10βˆ’3 𝜈 = 10βˆ’4 𝜈 = 10βˆ’5 60 Γ— 60 42 (46.82) 38 (45.70) 37 (44.24) 32 (42.36) 22 (26.82) 120 Γ— 120 49 (64.38) 42 (57.22) 39 (60.52) 35 (54.81) 24 (36.55) 180 Γ— 180 51 (97.51) 46 (94.32) 40 (87.26) 36 (82.64) 26 (54.81) 260 Γ— 260 56 (157.84) 49 (150.57) 45 (142.81) 37 (124.83) 29 (91.29) 520 Γ— 520 62 (287.86) 59 (283.35) 57 (280.55) 49 (232.36) 40 (187.65) The expression in the right-hand side of (26) is linear in the variables at the new time level and possesses a discretization error 𝑂(𝛿𝑑 2 ). After discretization and linearization of the problem (20)-(23) we need to solve large sparse strongly nonsymmetric SLAEs at each time step. The choice of an appropriate iterative method and its imple- mentation largely determine the overall efficiency of the computational algorithm. In practice, clas- sical iterative methods are used, such as BiCG, GMRES, GMRES(m) [15] and others. We use the AMG with the special STS- and PSTS- smoothers. Both of the smoothers are compared with the Gauss- Seidel (GS) one which is the classic smoothing method [3]. The number of the pre-smoothing and post-smoothing steps [2] equals 15 for all methods. Calculations are carried out using the V-cycle [2, 3] as the least demanding from the computational point of view. For the PSTS method, 𝐡𝑐 = 𝐼 and 𝐹0 was chosen such that the matrix (𝐹0 + 𝐾𝐿 ) was unitary [8]. In the numerical experiments, the time step is equal to 0.01 and the Reynolds number is 𝑅𝑒 = 𝜈 βˆ’1 . In Tables 1, 2, 3 we compare the AMG-iterations with the STS-, PSTS- and GS-based smoothers on the different grids for small viscosity values. In our implementations, all iterations are started from the zero vector, and terminated either when β€–π‘Ÿ (𝑝) β€–2 ≀ 10βˆ’6 , β€–π‘Ÿ (0) β€–2 where π‘Ÿ (𝑝) = 𝑏 βˆ’ 𝐴𝑒 (𝑝) is the residual vector of the SLAE (1) at the current iterate 𝑒 (𝑝) , and π‘Ÿ (0) is the initial residual. Our comparisons are done for the number of iteration steps (denoted by β€œIT”) and the elapsed CPU time (in seconds and denoted by β€œCPU”). The abbreviation β€œn.c.” in the tables means β€œno convergence”. The experiments are run in MATLAB (version R2018b) with a machine precision 10βˆ’16 . From Tables 1, 2, 3 an increase in the number of iterations follows with an increase in the mesh size for all tested methods. This is due to some features of the algebraic approach in the multigrid method (namely, the PMIS algorithm in the AMG). The traditional scalable approach (RS algorithm) works well for many 2-D problems. In this case, the solver can be obtained with the number of iterations that does not depend on the size of the problem. But when using AMG interpolation in combination Table 3 AMG+GS Iterations and CPU with Different 𝜈. Grid 𝜈 = 10βˆ’1 𝜈 = 10βˆ’2 𝜈 = 10βˆ’3 𝜈 = 10βˆ’4 𝜈 = 10βˆ’5 60 Γ— 60 39 (39.57) 46 (42.51) 54 (114.85) n.c. n.c. 120 Γ— 120 57 (83.82) 64 (122.61) 83 (160.50) n.c. n.c. 180 Γ— 180 59 (162.36) 82 (185.38) 85 (192.20) n.c. n.c. 260 Γ— 260 n.c. n.c. n.c. n.c. n.c. 520 Γ— 520 n.c. n.c. n.c. n.c. n.c. with the PMIS, the AMG convergence deteriorates depending on the problem size [5]. However, the RS algorithm has obvious drawbacks, such as poor efficiency for 3-D problems [5]. Moreover, from Tables 1, 2, 3 it follows that the AMG method with the STS- and the PSTS- smoothers has fast convergence speed for all tested values of the viscosity coefficient (𝜈 = 10βˆ’1 Γ· 10βˆ’5 ) on all used grids, while the AMG with the Gauss–Seidel smoother does not converge for 𝜈 = 10βˆ’4 , 10βˆ’5 on all considered grids, and does not converge on the grids 260 Γ— 260 and 520 Γ— 520 nodes for all val- ues of the viscosity coefficient. For all tests, the AMG+STS and AMG+PSTS methods outperform the AMG+GS method with respect to both number of iteration steps and CPU time. The AMG+STS method slightly ahead of the AMG+PSTS method, apparently due to the fact that the first method belongs to the class of triangular, and the second one relates to the class of two-step iterative methods. In the numerical tests the iterative parameters for all methods was chosen exper- imentally. The decrease in the number of the AMG+STS and AMG+PSTS iterations with a decrease in the viscosity coefficient values is explained that these methods converge the faster, the more pro- nounced the nonsymmetry of the SLAE is. The number of iterations and CPU time of the AMG+GS method increases with decreasing 𝜈. 5. Conclusions Numerical experiments carried out for the 2-D model incompressible unsteady Navier-Stokes problem have shown efficiency of the STS- and the PSTS- smoothers for the AMG and their advantage over the classical Gauss-Seidel smoother. The use of these smoothers for the AMG allows solving this problem for small values of the viscosity coefficient. The lack of convergence of the AMG+STS and the AMG+PSTS methods for the smallest values of 𝜈 is not observed, which could be expected for the unsteady problem. The obtained SLAE is simpler than for the stationary problem, since the spectrum of this system lies in the right half-plane. This is due to a shift in the spectrum when the time derivative is approximated. It is supposed to solve the corresponding 3-D problem in the future using these methods. Acknowledgments This research was funded by the Government of the Russian Federation (Project No.075-15-2019-1928), Russia. References [1] R. D. Falgout, An introduction to algebraic multigrid, Computing in Science & Engineering 8 (2006) 24–33. doi:10.1109/MCSE.2006.105. [2] W. Hackbusch, Multi-grid methods and applications, Springer-Verlag, Berlin, 1985. [3] U. Trottenberg, C. Oosterlee, A. Schuller, Multigrid, 1st. ed., Academic Press, 2001. [4] K. Stuben, Algebraic multigrid (amg): experiences and comparisons, Applied Mathematics and Computation 13 (1983) 419–451. doi:10.1016/0096-3003(83)90023-1. [5] H. de Sterck, U. M. Yang, J. J. Heys, Reducing complexity in parallel algebraic multigrid preconditioners, SIAM Journal on Matrix Analysis and Applications 27 (2006) 1019–1039. doi:10.1137/040615729. [6] L. A. Krukier, Convergence acceleration of triangular iterative methods based on the skew- symmetric part of the matrix, Applied Numerical Mathematics 30 (1999) 281–290. [7] L. A. Krukier, L. G. Chikina, T. V. Belokon, Triangular skew-symmetric iterative solvers for strongly nonsymmetric positive real linear system of equations, Applied Numerical Mathematics 41 (2002) 89–105. [8] L. A. Krukier, T. S. Martynova, Z.-Z. Bai, Product-type skew-hermitian triangular splitting iter- ation methods for strongly non-hermitian positive definite linear systems, Journal of Computa- tional and Applied Mathematics 232 (2009) 3–16. doi:10.1016/j.cam.2008.10.033. [9] L. Wang, Z.-Z. Bai, Skew-hermitian triangular splitting iteration methods for non-hermitian positive definite linear systems of strong skew-hermitian parts, BIT Numerical Mathematics 44 (2004) 363–386. doi:10.1023/B:BITN.0000039428.54019.15. [10] G. V. Muratova, L. A. Krukier, E. M. Andreeva, Fourier analysis of multigrid method with trian- gular skew-symmetric smoothers, Communication on Applied Mathematics and Computation 27 (2013) 355–362. doi:10.3969/j.issn.1006-6330.2013.03.007. [11] Z.-H. Cao, Convergence of multigrid methods for nonsymmetric, indefinite problems, Applied Mathematics and Computation 28 (1988) 269–288. doi:10.1016/0096-3003(88)90076-8. [12] J. Mandel, Multigrid convergence for nonsymmetric, indefinite variational problems and one smoothing step, Applied Mathematics and Computation 19 (1986) 201–216. doi:10.1016/ 0096-3003(86)90104-9. [13] J.-L. Guermond, C. Migeon, G. Pineau, L. Quartapelle, Start-up flows in a three-dimensional rectangular driven cavity of aspect ratio 1:1:2 at re=1000, Journal of Fluid Mechanics 450 (2002) 169–199. doi:10.1017/S0022112001006383. [14] S. McKee, M. F. Tome, V. G. Ferreira, J. A. Cuminato, A. Castelo, F. Sousa, N. Mangiavacchi, The mac method, Computers & Fluids 37 (2008) 907–930. doi:10.1016/j.compfluid.2007.10. 006. [15] Y. Saad, Iterative methods for sparse linear systems, 2nd. ed., SIAM, 2003.