Algebraic multigrid method with skew-Hermitian smoothers

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 producttriangular 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 twodimensional 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 eficiency of the algebraic multigrid method with the skew-Hermitian smoothers for this problem.

Keywords

• eol>algebraic multigrid method
• skew-Hermitian smoothers
• unsteady Navier-Stokes problem

1. Introduction

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 diference approximation of the 2-D incompressible unsteady Navier-Stokes equations governing the flow of viscous Newtonian lfuids. We use two iterative methods based on the Hermitian and skew-Hermitian splitting of the matrix 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) (2) (4)

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, where 0 = 1 ( + ∗), 1 = 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 coeficient matrix of the SLAE ( 1 ) is real.

Then, we can split the skew-Hermitian part 1 of the matrix ∈ ℂ × into

1 = + , where and are the strictly lower and the strictly upper triangular parts of 1, respectively. Obviously, that = − ∗ . where where and are two acceleration parameters, and ( ) is defined by ( +1) = (,

) ( ) + ( )−1, (, ) = ( )

−1( ( ) − ), ( ) = + ((1 + ) + (1 − ) ), = ±1.

( +1) = (,

) ( ) + ( )−1, (, ) = ( )

−1( ( ) − ), ( ) = ( + 2 ̂ ) −1 ( + 2 ̂ ), . For = 0, 1, 2, ... until { ( )} convergence, compute

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 coeficient given in [9], the PSTS iteration method has been investigated in [8]. matrix . These methods have been further developed in [9]. Moreover, based on the methodology

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 efective 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 where

( ) = 0( ) + 1( ), { 0( ) = 12 ( ( ) + ( ) ) = + ( 2 )2 ̂ −1 ̂ , 1( ) = 21 ( ( ) − ( ) ) = 2 ( ̂ + ̂ ) = 2 1.

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 efect 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

defined by ( 6 ). That is to say, Let 0( ) and 1( ) be, respectively, the Hermitian and the skew-Hermitian parts of the matrix ( ) and are two acceleration parameters, and ( ) is defined by with ∈ ℂ × a prescribed Hermitian matrix. . For = 0, 1, 2, ... until { ( )} convergence, compute

The PSTS iteration method [8]: Given an initial guess (0), and two positive parameters and ( 5 ) ( 6 ) ( 7 )

Suppose that

0 < ℎ ≤ 0 ≤ ℎ , ≤ ̂ −1 ̂ ≤ , 0 < ≤ ≤ , ≤ 0( ) ≤ , Notice, that and

satisfy the following inequalities: values of the corresponding matrix, respectively [8]. ̂ −1 ̂ ℎ, , , where is an identity matrix. Since 0

and are Hermitian positive definite matrices, 0( ) and 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 = ,

We need to require positive definiteness of the matrix ( ) from ( 7 ). Since the matrix ̂ −1 ̂ is negative definite, then 0( ) > 0, if the parameter ∈ (0, ), where

√ convergent, i.e., the spectral radius ( (, the acceleration parameters and satisfy and Theorem 1. [8] Let the matrices and ( ) be positive definite. Then the PSTS iteration method is )) of its iteration matrix (, ) is less than 1, provided that ≥ + ( 2 ) , ≤ + ( 2 ) .

2

2 0 < < ,

0 < < , 0 < <

Θ = 2 Θ ,

ℎ + ( 2 )

We use some denotations and theoretical results from [11, 12].

We are interested in solution of the problem (1) in

and the inner product is denoted by (.,.) with ∥ . ∥ corresponding norm in , = 1, 2...

Let 1 ⊂ 2 ⊂ ... be a family of nested finite-dimensional linear spaces. The dimension of is 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 = −1 ̃ is satisfied ( ) = 1.

Using the operators ̃ and

and the inner product

Then we may define the following norms on : we define the energy norm ‖ ‖ ̃ = ( ̃ , )1/2, ∀ ∈ , (, ) = ( , ), ∀, ∈ . ∥ ∥, =∥ ∥ = ( ,

∥ ∥1, = ( ̃ , )1/2 =∥ ∥ ∥0, = ( , )1/2 = (, )1/2 . 1/2 ) , = 0, 1, ...

∥ ̃ ,

Now we make three basic assumptions [11, 12]:

Define the subspace

⊂ , = { ∈ ∶ ( , ) = 0, ∀ ∈ −1} Assumption 1. There exists ∈ −1, , 0 < ≤ 1 and < ∞ such that for ∀ ∈ Assumption 2. There exists , → 0 ( → ∞) such that ∥ − ∥1, ≤ ∥ ∥1+,

2 2 ∣ ( ̂ , ) ∣< ∥ ∥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 suficiently small.

In [12] it was shown that assumptions 1-3 were satisfied for suficiently wide class of boundary problems in two dimension-limited domains with diferent 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 .

MGM-iteration. We need to estimate the energy-norm of the contraction number defined by

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 = sup ∥ 2 − ∗ ∥1 ∥ 0 − ∗ ∥1 , 0 ≠ ∗.

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 ∥2 1 ( 12 ) with = from Assumption 3. Then ≤ ̃ for the two-grid contraction number where ⎧ ⎪ ⎪ ⎪ ⎩ ⎪ 1− ̃ ≡ ̃( ) ≡ sup ⎨ ( 1+ −1 ( 1+ )2/ 2/ (1 + Δ))

2+ 2 2+2 ⎪ ⎪ 2 + 2 − 2

≤ 1, , ≥ 0 = ∥ 1 + ∗ ∥1 ∥ 1 ∥1 , = ∥ ∗ ∥1 ∥ 1 ∥1 , 1/2 ∶ ⎪⎪ ⎫ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ assumptions 1, 2, 3 respectively, while is calculation accuracy.

where constants and Δ depend on properties of smoothing system and constants , , , taken from

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]. then the smoothing assumption ( 12 ) is satisfied with

Now we can prove the convergence of the method.

̃+ ̃∗ − ̃ ≥ ( ̃ − ̃)∗ −1 ( ̃ − ̃),

Δ = (1 + ) ( + 2) , = ( ( ̃−1 ̃( ̃−1)∗ ̃)) .

1/2 satisfied where 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

Denote then from ( 8 ) - ( 11 ) it’s easy to obtain that

̃+ ̃∗ − ̃ = 2/ 0 − 0, (

̃ − ̃)∗ ̃−1 ( ̃ − ̃)= ( 0 − 1 ∗) 0−1 ( 0 − 1 = 0 − 1 ∗ − 1 + 12 ∗ 0 = 0 − 2 0 + 12 ∗ 0 −1 −1

)= 2 ℎ =

0 − 0, + ( 2 )2 < , 2 = ∗ > 0.

Theorem 4. For the method ( 5 ) - ( 6 ) there exists the constant > 0 so that inequality ( 13 ) is satisfied. left and right parts of it. Using ( 7 ) we obtain

Proof. In case of using method ( 5 )-( 6 ) ̃= 1/ , = 0. Consider inequality ( 13 ) and transform ( 13 ) ( 14 ) ( 15 ) (16) (17) (18) thus, 2/ 0 > 0 for acceleration parameters ,

satisfying the conditions of the Theorem 1, and

From (15) - (17) we obtain that for inequality (13):

Denote = 1

0−1/2 −1/2. Then (18) is transformed to Multiplying the left and the right parts of this inequality on −1/2 we obtain ≥ (− + 2 1 ∗ 0−1 ). ≥ (− + 2 1 −1/2 ∗ 0 −1/2 .

−1 ) If we take

≥ ( ∗ − ) . =

1 ∥ ∗ − ∥

, + + + + + + − − + = 0, 1 1 2 2 ( 2 + ( 2 + 2 2 ) 2 2 ) = 0, = 0, = 0, = 0, ( = 0, = 1, = 0),

= 1, = 0, ( = 1), (, , 0) = 0, (, , 0) = 0. 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 obtained for the two-grid method can be extended to the multigrid method.

PSTS-smoothers ( 6 ), is calculated by (19) and Δ is calculated by ( 14 ) correspondingly. The results

4. Numerical experiments

Let us consider the 2-D unsteady Navier-Stokes equations of viscous incompressible fluids. The primitive variables mathematical formulation of this problem is: find V and (up to constant) so that V − ΔV + (V ⋅ ∇)V + ∇ = f

Ω × (0, ],

V = 0 V = g Ω × [0, ], Ω × [0, ], V(x, 0) = V0(x) Ω gence; V = ( the boundary (, , Ω and V

0 ), (, , occurs in the time interval [0, ]. where is the kinematic viscosity coeficient; Δ is the Laplacian, ∇ is the gradient, is the diveris the divergence-free initial velocity field. We assume that the fluid motion

)) is the velocity vector, is the pressure, g is the velocity prescribed on

One of the main dificulties in using the AMG is the constructing smoothers that are robust over a wide range of the viscosity coeficient, especially for small values of . Classical MGMs were developed 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 (19) (20) (21) (22) (23) where is the Reynolds number ( = / , is a characteristic velocity of the flow, is a characteristic 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 efect of the suggested smoothers for the MGM.

The fully implicit scheme has been used for the time discretization. The use of implicit schemes removes 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 (V( +1) − V( )) + ∇ ( +1) = −(V( +1) ⋅ ∇)V( +1) + 1 ΔV( +1),

V( +1) = 0, (24) (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 twodimensional 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 diference 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 diferent variables are calculated at diferent 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 linearization around the old time level is used. If we want to linearize a nonlinear term ( +1) ( +1), then ( +1) ( +1) = ( ) ( +1) + ( +1) ( ) − ( ) ( ) + ( 2). (26) 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 implementation largely determine the overall eficiency of the computational algorithm. In practice, classical 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 GaussSeidel (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 diferent 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 with the PMIS, the AMG convergence deteriorates depending on the problem size [5]. However, the RS algorithm has obvious drawbacks, such as poor eficiency for 3-D problems [5].

Moreover, from Tables 1, 2, 3 it follows that the AMG method with the STS- and the PSTSsmoothers has fast convergence speed for all tested values of the viscosity coeficient ( = 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 values of the viscosity coeficient. 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 experimentally. The decrease in the number of the AMG+STS and AMG+PSTS iterations with a decrease in the viscosity coeficient values is explained that these methods converge the faster, the more pronounced 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 eficiency 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 coeficient. 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 skewsymmetric 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 iteration methods for strongly non-hermitian positive definite linear systems , Journal of Computational 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 triangular 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 .