The Mixed Finite Element Multigrid Method for Stokes Equations

The stable finite element discretization of the Stokes problem produces a symmetric indefinite system of linear algebraic equations. A variety of iterative solvers have been proposed for such systems in an attempt to construct efficient, fast, and robust solution techniques. This paper investigates one of such iterative solvers, the geometric multigrid solver, to find the approximate solution of the indefinite systems. The main ingredient of the multigrid method is the choice of an appropriate smoothing strategy. This study considers the application of different smoothers and compares their effects in the overall performance of the multigrid solver. We study the multigrid method with the following smoothers: distributed Gauss Seidel, inexact Uzawa, preconditioned MINRES, and Braess-Sarazin type smoothers. A comparative study of the smoothers shows that the Braess-Sarazin smoothers enhance good performance of the multigrid method. We study the problem in a two-dimensional domain using stable Hood-Taylor Q 2-Q 1 pair of finite rectangular elements. We also give the main theoretical convergence results. We present the numerical results to demonstrate the efficiency and robustness of the multigrid method and confirm the theoretical results.


Introduction
This study considers the numerical solution of the large scale linear algebraic system arising from the discretization of the partial differential equations. The discretization is achieved by the finite element method. For positive definite linear systems, linked to the Poisson equations, the multigrid (MGM) methods are proven to be the most effective and fast methods [1,2]. However it is more challenging for linear indefinite algebraic systems. In this paper we consider multigrid methods for solving linear indefinite algebraic system of equations arising from the mixed finite element discretization of the steady state Stokes problem: div u = 0, in Ω, where u is a velocity field, represents pressure, and f is an external force field. The problem is considered with (1)-(3) defined on the domain Ω ⊆ R 2 with boundary Ω.
The main goal of this work is to construct and analyze numerical methods that produce an appropriate solution to the Stokes problem. The main thrust is to apply an iterative method, the multigrid method, to solve the linear system of equations that arise from the discretization of the Stokes equations. The MFEM applied to (1)-(3) with carefully chosen finite element spaces results in the algebraic system which must be solved. The velocity variable u together with the pressure variable is the solutions of the system. We discretize the domain of the Stokes problem by the rectangular grids with a pair of conforming mixed finite element spaces that are inf-sup stable. In our experiment we use Hood-Taylor 2 -1 pair as used by [3]. The process produces a symmetric indefinite system of linear algebraic equations. In this paper we study an efficient solver for this system. This work on multigrid method has been motivated 2 The Scientific World Journal by the need to effectively and efficiently solve large application problems. The multigrid method has been shown to be very efficient and successful in solving control problems [1,2,[4][5][6] and elliptic partial differential equations [7][8][9] in an accurate and computationally efficient way. The multigrid method has been applied to problems discretized by the finite difference method and widely by finite element method [3,8,[10][11][12][13][14]. The effectiveness of the multigrid method depends on the correct choice of the smoothers. Various smoothers have been suggested in literature weighted Jacobi, Gauss Seidel [11], Ilu [8], Vanka-type [9,12,13,15], Braess-Sarazin-type ( [13,[16][17][18]), Semi implicit method for pressure linked equations [15], SOR/Richardson [18], and inexact Uzawa [18]. It is the purpose of this study to apply the multigrid solver to the Stokes problem with the following iterative solvers as smoothers: Braess-Sarazin, inexact Uzawa, preconditioned MINRES, and the distributed Gauss Seidel. The inner solver of these smoothers can also be taken as the multigrid method for the definite subsystems. There is no work known where a comparative study is made on the effects of these four smoothers on the performance of the multigrid method for indefinite systems. The first step is to transform the continuous problem to the discrete system and apply the MFEM that produces the linear algebraic system on which the multigrid method is developed, analyzed, and finally numerically and computationally implemented.
The key features and ingredients of the multigrid method are smoothing and coarse grid correction that involves the intergrid transfers and a solution correction step. The main results of the work are the convergence of the multigrid method in calculating the velocity and pressure variables in an appropriate norm which is based on the smoothing and approximation properties [9,18]. The rest of the paper is organized as follows. In Section 2 we give the discrete system of the Stokes problem by mixed finite element method. In Section 3 the iterative solution technique, the geometric method, and smoothers are outlined. The known theoretical convergence analysis results are also outlined. In Section 4 a numerical experimental and comparative analysis on the effects of smoothers on the performance multigrid method is presented and discussed and the conclusion is given.

The Stokes Discrete System
For the discretization of the Stokes equations in the domain Ω we need to transform the system (1)-(3) to the weak variational form. For the weak variational formulation of the Stokes equations we define the following solution and test spaces: By multiplication of the first equation (1) with v ∈ 1 0 and the second equation (2) with ∈ 2 (Ω), subsequently integrating over the domain Ω, applying the Gauss theorem, and incorporating the boundary condition (3), we obtain the variational form.
The Scientific World Journal 3 2.1. The Mixed Finite Element Discretization. The mixed finite element discretization of the weak formulation of the Stokes equations produces a linear algebraic system of equations. The finite element method described here is based on [7,16,19,21,23,25]. We will introduce the concept of mixed finite element methods. Details can be found in [7,[19][20][21]25].
We assume that Ω ⊆ R 2 . We define the finite dimensional spaces. Let ℎ and V ℎ be subspaces of and V, respectively. Now we can formulate a discrete version of problem (5).
The finite element discretization should satisfy the inf-sup condition. The following theorem shows that again the infsup condition is of major importance (for the proof we refer to [22]).

Theorem 2.
Assume that is ℎ -elliptic (with h independent ellipticity constant) and that there exists a constant > 0 (independent of ℎ) such that the discrete inf-sup condition holds. Then the associated (discretized, steady state) Stokes problem has a unique solution (u ℎ , ℎ ), and there exists a constant 1 such that If the basis of ℎ is given by { 1 , . . . , } and of ℎ is given by { 1 , . . . , }, then where is the number of inner nodes and is the number of boundary nodes so the coefficients u : = + 1, . . . , + interpolate the boundary data and = + . The mixed finite element entails partitioning of the solution domain Ω into quadrilaterals; in our case that is Ω = ∪ we denote a set of rectangular (square) elements by ℎ = { 1 , 2 , 3 , . . .} and on each element and we denote the space ( ) of degree less than or equal to . There are a variety of finite element pairs whose effectiveness is through stabilization [26]. In this work we are going to use Hoods-Taylor 2 -1 square finite elements which are known to be stable.
We specify An element (u ℎ , ℎ ) ∈ ℎ × ℎ is uniquely determined by specifying components of u ℎ on the nodes and on the midpoints of the edges of the elements and the values of ℎ on the nodes of the elements. The mixed finite element method results in the coupled linear algebraic system which has to be solved by the appropriate solvers. The resulting system is with ℎ being a block Laplacian matrix and ℎ being the divergence matrix whose entries are given by The entries of the right hand side vector are The linear algebraic system can be represented as where The discretization and assembling of matrices are done by the MATLAB implementation of the mixed finite element method [8]. ℎ is stiffness matrix resulting from the discretization of the Laplacian. The resultant coefficient matrix is large, sparse, indefinite and the system must be solved iteratively, in this case by multigrid solvers. The multigrid solver is a well known fast solver for the elliptic partial differential equations [2,5].

Multigrid Method
The main focus of this section is the construction of the multigrid solver to find the approximate solution of (20) at the finest mesh/discretization. Let (V × ) be a sequence of subspaces of the finite dimensional subspaces ℎ and ℎ defined on sequence of grids ∈ {0, 1, 2, 3, . . . , max } with mesh sizes ℎ 0 , ℎ 1 , ℎ 2 , . . . , ℎ with ℎ +1 := (1/2)ℎ . We define a hierarchy/family of nested finite element subspaces for the velocity and pressure: where The Scientific World Journal At the discrete level with the defined discrete spaces and bases, the linear algebraic system is defined by where M := [ ], := [u ], and := [ f ]. The main goal is to find the pair = (u , ) of the discrete velocity and the discrete pressure variables at the finest level . Now we introduce the multigrid iteration for solving the discretized equation (22) on grid . We define the multigrid algorithm at level as is the output of velocity and pressure after one step of the multigrid algorithm at level ; (ii) u old is the input velocity at level ; (iii) old is the input pressure at level ; (iv) u, , −1 and , , −1 are restriction operators for velocity and pressure, respectively, from level to level − 1; (v) u, −1, and , −1, are prolongation operators for velocity and pressure, respectively, from level − 1 to level .

end (5) Correction
Step define the new iterate: Postsmoothing. Starting with (u * , * ) perform 2 smoothing steps using smoothing operator S to produce (u new , new ) : The multigrid method described above belongs to a class of optimal order methods for solving linear systems emanating from the discretization techniques like the finite element method. Its known convergence speed does not deteriorate when the discretization is refined whereas classical iterative solvers slow down for the decreasing mesh size [1,2,5,6]. The starting point of the multigrid concept is the observation that classical iteration methods have some smoothing properties. The operator S represents such methods; in this study it represents Braess-Sarazin, inexact Uzawa, distributed Gauss Seidel, and the preconditioned minimum residual method. These methods are characterized by poor/slow global convergence rates and for errors whose length scales are comparable to mesh sizes, they provide rapid damping leaving behind smooth, longer wave length errors. These smooth parts of the error are responsible for the poor convergence. A geometric multigrid method involves a hierarchy of meshes and related discretization. A quantity that is smooth on a certain grid can also be approximated on a coarser grid. Low frequency error components can be effectively reduced by a coarse grid correction procedure. Since the action of a smoothing iteration leaves only smooth error components, it is possible to represent them as the solution of an appropriate coarser system. Once this coarser problem is solved, the solution is interpolated back to the fine grid to correct the fine approximation for its low frequency errors.
The Scientific World Journal 5 The most essential ingredients of the multigrid method are the smoothing operator, for which using a wrong smoother will destroy the efficiency of the entire multigrid method, and the coarse grid correction which involves the prolongation and the restriction operators. In multigrid methods we have to transform information from one grid to another and for that purpose we use prolongations and restrictions operators.
Restriction transfers values from fine grid to the next coarse grid. Prolongation transfers values from the coarse grid to the next fine grid. Next we discuss the key components of the multigrid method.
(a) Intergrid transfer operators: the intergrid transfer operators are the restriction and prolongation between different grid levels. The restriction operator maps the residual from the finer grid to a coarser grid while the prolongation operator transfers vectors from coarse grid to fine grid. The restriction between levels and − 1 is defined by where the restriction operators (u, , −1) : R → R −1 and ( , , −1) : R → R −1 for velocity and pressure, respectively. The prolongation between levels − 1 and is defined again as where the prolongation operators (u, −1, ) : R −1 → R and ( , −1, ) : R −1 → R are representations of the following relations V −1 ⊂ V for the quadratic interpolation of the velocity ( 2 ) and −1 ⊂ for the linear interpolation of the pressure ( 1 ).
(b) Coarse grid correction: the other key ingredient of the multigrid method is the coarse grid correction.
In the multigrid solution process we need to solve the problem at the finest define level = max . The problem is defined on the coarser grid levels and on the coarsest grid level the problem is solved exactly.
There are very few situations in which a grid can be coarsened to the extent that it is not practical to solve the problem using a direct method but iteratively. In this work the iterative solver used as a smoother is applied to solve the problem at the coarsest level.

The Smoothers.
The most crucial part is the proper choice of a smoothing technique. Usually, the well-known smoothing iterations for the scalar problems (damped Jacobi or Gauss-Seidel relaxation) are not appropriate for saddle point problems or are even not defined, for example, in saddle point systems like (22). There are natural ways to generalize scalar smoothing schemes to systems of PDEs. The smoothing process is the main ingredient of the multigrid method. The convergence of the multigrid method is influenced by the smoothing process [11,14,18]. We perform a number of iterations of an iterative solver to smooth the residual. The main goal is to compare the effectiveness of different iterative schemes as smoothers of the multigrid methods. On each level of a multigrid method, a system involving operator S has to be solved approximately. The smoother dumps out highly oscillating error modes of the systems. In this paper we consider the following smoothing process: Several smoothers have been proposed and applied in literature. Brandt [4] advocates for the use of the distributed Gauss Seidel smoothing. The Vanka-type smoother is widely used with a coupled Gauss Seidel scheme [13,14] that introduces the idea of transforming smoothers and combines with incomplete factorization to develop an efficient smoothing. John and Tobska [14] and Pernice [15] used the Braess-Sarazin-type smoother with the Schur complement schemes as smoothers which exhibit wonderful smoothing properties.
The following algorithms describe the iterative schemes that are used as smoothers in this study.

Braess-Sarazin-Type
Smoother. The Braess-Sarazin smoothers proposed in [17] and used in [13,18] solve a large saddle point problem in each smoothing step. This Braess-Sarazin or SIMPLE-type iteration uses (̂) as a smoother for the saddle point problem (22). The smoother as presented in [17] and generalised in [18] consisted of constant application of the smoothing iteration: witĥ= diag( ) and = 2 given. The smoothing Braess-Sarazin iteration (35) solves the auxiliary problem with r = u + − f and = u − . Inherent in the system system (36) is the problem of the auxiliary pressure variablê̂̂=̂− This system is solved approximatively by an iterative solver.
From the system we get̂approximately which can be used to approximately determineû fromû End.
This saves the application of̃− 1 at the end of every outer iteration and hence improves the efficiency of the algorithm. The other variants of the inexact Uzawa method are analysed in [26][27][28].

The Distributed Gauss Seidel Type Smoothers (DGS).
The standard smoothing iteration schemes like Jacobi and Gauss Seidel smoothers are not applicable to the system (22) because of the nature of the coefficient matrix; particularly the zero block in the diagonal hampers the smoothing process. The distributive smoother transforms the vital operators to the main diagonal and applied as a decoupled smoother. The DGS was introduced in [4] is related to a successive application of standard Gauss Seidel applied to the matrix operator M (22) and G = ( − ) with MG = ( ).
We solve the transformed residual equation witĥand̂being invertible approximations of and := , respectively. A single iteration with the update through a distributive matrix G is performed by the following algorithm.
(2) Smooth the transformed continuity equation (3) Transform the correction back to the original variables u +1 = w + , The DGS has been widely used as a smoother for the finite difference discretization. In this paper the DGS type smoothers are used for finite element discretization of the Stokes problem.

The Preconditioned Minimum Residual Smoother.
The preconditioned minimum residual method is a Krylov subspace method for solving symmetric indefinite systems and uses popular block preconditioners. This method is used as a smoother for the multigrid method of the Stokes problem in this paper. For the Stokes equations, the classical blockdiagonal preconditioner for MINRES method [8] is witĥ=̂− 1 . The block preconditioning requires the solution of two systems of equations with matriceŝand at each MINRES iteration. If −1 is computed exactly, the preconditioned Krylov methods converge in two or three steps [10]. For practical implementations, the Schur complement̂is replaced by the mass matrix of the pressure space. For discontinuous pressure space, is block diagonal and easy to invert. For continuous pressure space, say 1 , the mass matrix can be further replaced by its diagonal matrix [8].

Multigrid Convergence.
The convergence analysis of the multigrid method relies on the two properties, namely, the approximation and the smoothing. The general convergence rates are independent of ℎ (the mesh size), is the level of discretization, and 1 and 2 are the number of preand postsmoothing iterations [1,12,18]. The results for the convergence of the multigrid method for the scalar elliptic problems cannot apply to the Stokes equations. We provide a snapshot of the available convergence results of the multigrid method for Stokes equations. The ideas presented in this paper are based on the work in [12,16,18]. An iteration of single multigrid step consists of a combination of smoothing step and a coarse grid correction step. We will consider the multigrid convergence with the Braess-Sarazin smoother with S being the iteration matrix of the smoother (34) and the M being the Stokes stiffness matrix in (22). The operator and its adjoint are intergrid transfer operators, prolongation, and the restriction, respectively. The convergence analysis of the multigrid method begins with the analysis of the two-grid method, with 1 and 2 the pre-and post-smoothing steps, respectively, applied to (22) results in the iteration matrix The key point on the analysis of the multigrid method is that the error can be split into two components. That is the one produced by the smoothing process and the one produced by the coarse grid correction. The coarse grid error consists of the low frequency components and the smoothing consists of the high frequency components of error. The ability to cope with the low frequency components is called the approximation property and with the second is called the smoothing property. For the analysis of the multigrid convergence [2,29] used the framework based on the smoothing and approximation property. For analysis we The Scientific World Journal 7 define the following norms, Euclidean norm by ‖ ⋅ ‖, and on R + the following norm is applied: with Θ := ( ) . (46) Furthermore we introducê Using the norms defined above and taking 1 = and 2 = 0 above we obtain The theorems below state the two properties and the multigrid convergence. For detailed proof we refer to [12,18].

Theorem 6 (approximation property). Assume that Ω is such that the problem (5) is 2 -regular. Let M be the coefficient stiffness matrix and and the prolongation and the restriction operators. Then there exists a constant M independent of and using ℎ-scaling induced by M then
The smoothing property is dependent on the smoother used. It varies from one smoother to another. In this work we used the Braess-Sarazin in which we solve the system (37) exactly and sufficiently accurate inexact inner solver.
Combining the approximation property Theorem 6 with the smoothing property Theorem 7 produces a two-grid convergence result. Theorem 8. Assume that 2 = 0 and that Ω is such that the problem (5) is 2 -regular. Then for the two-grid method the following holds: with a constant M independent of l and m.
Using this two-grid contraction number bound the multigrid -cycle method convergence results can be derived using ideas in [1,2].

Numerical Results
In this section we present the numerical solution of classical Stokes problem (1)-(3) using the solver presented above. The solver is denoted by MGM (Algorithm 3). We present the results of this method as outlined above to run the traditional test problem, the driven cavity flow problem [11,12,27,28]. It is a model of the flow in a square cavity (the domain is Ω ◻ ) with the top lid moving from left to right in our case the regularized cavity model { = 1 : −1 ≤ ≤ 1 | = 1 − 4 } [11]. The Dirichlet no-slip boundary condition is applied on the side and bottom boundaries. The mixed finite element method was used to discretize the cavity domain Ω = (−1, 1) 2 .
We pay particular attention to the computational performance of the multigrid method on the system (22) at different grid levels. We compare the effectiveness of different smoothing/relaxation methods in the performance of the multigrid method and different approximations for the preconditionerŝand̂of the smoothers. The following setup of the smoothers listed is considered.
(i) Distributed Gauss Seidel (DGS) smoother: we use one Gauss Seidel iteration for the evaluations of̂and one Gauss Seidel iteration for the computation of̃. The method becomes DGSMG.
(ii) Inexact Uzawa smoother (IUzawa): the two cases are considered for the evaluation of the preconditioners. Firstly, the approximation̂= diag( ) and one V(1, 1)-cycle is used to approximate Schur compliment matrix̂. The second case is to use one V(1, 1)-cycle for both evaluations of̂and̂. The method becomes IUZAWAMG.
(iii) Braess-Sarazin smoother (B-S): the two cases are conspired for the evaluation of the preconditioners. Firstly, the approximation̂= diag( ) and one V(1, 1)-cycle is used to solve the approximate Schur compliment matrix̂. The second case is to use one V(1, 1)-cycle for both evaluations of̂and̂. The method becomes B-SMG.
(iv) PMINRES smoother: the first case is to use diagonal preconditioner for̂and̂and the second case is one V(1, 1)-cycle for committing the inversion of the Laplacian operator for velocity as one V(1, 1) cycle is used to approximate the Schur compliment using the pressure mass matrix to accelerate the MINRES. The method becomes PMINRESMG.
The comparison is made on the performance of the multigrid schemes with different smoothers (i)-(iv) and cases involving different approximations of preconditioners in terms of iterative counts and CPU time. The numerical treatment is given to the discrete Stokes problem which resulted from the mixed finite Hood-Taylor stable elements consisting of biquadratic elements for the velocities and bilinear elements for the pressure, on a uniform grid. Implementation of our algorithms was performed on a Windows 7 platform with 2.13 GHz speed intel dual core processor by using MATLAB 7.14 8 The Scientific World Journal   programming language and the MATLAB built-in Minres functions are used for the smoother. For the discretization we start with a uniform square grid with ℎ 0 = 1/2 and we apply regular refinements to this starting discretization to obtain the finest grid level. The discretized equations are solved using the multigrid iteration with the -cycle and -cycle and 1 and 2 being presmoothing and postsmoothing steps, respectively. The smoothers are determined by specifying approximationŝand̂as highlighted in (i)-(vi) and in all cases where one V-cycle inner multigrid iteration is used with 1 and 2 being Gauss Seidel iteration steps of the presmoothing and the postsmoothing, respectively.
In this work we use the structured mesh and regular refinements. The finite element matrices on the rectangular grids are assembled and the meshes are generated by the MATLAB IFISS toolbox [3] in a hierarchy of grids which are produced by successive regular refinements. We need to choose the coarse mesh (the starting mesh), the finest mesh which corresponds to the maximum level of refinement on which the final approximate solution is considered.
The assembled matrices are stored for each refinement level for the system (22). Table 1 shows an example of the refinement levels, we use the coarsest (starting) level to have 9 nodes for velocity and 4 nodes for pressure variables (level 1) but we start the computation at level 2. Table 1 shows the refinement levels and the number of grid points (nodes) for each level.
The zero initial guess is chosen for all the tests. In all the tests the iterations are repeated until the tolerance   pre-and postsmoothing steps with Braess-Sarazin (B-S) smoother (diag( ), V(1, 1)). We compare the performance of the -cycle and -cycle multigrid iterations with various smoothing steps at different grid levels using one of the smoothers, Braess-Sarazin. From Tables 2 and 3 we observe that the number of iterations decreases when the smoothing steps decrease and the CPU time increases as expected with the increase in smoothing steps. Tables 4 and 5 show the numerical results obtained of the multigrid solver at different grid levels. The number of -cycle and -cycle multigrid iterations and CPU time are shown, respectively. All the results presented underline the efficiency of the multigrid solver to indefinite systems of equations. In both tables we present results of the four studied smoothers of the multigrid solver. In Tables 4 and 5 we choose the approximation of the smoother preconditioners aŝ= 2 diag( ) and̂= V-cycle(1, 1) for Braess-Sarazin, IUzawa, and PMINRES. For the DGS we use the one Gauss Seidel for both and . We fix the number of smoothing steps to (3,3) for all the results in the tables.
Comparing the performance of the -cycle and -cycle multigrid solver, we observe that the smoothers have different effects on the performance of the multigrid solver. The multigrid solver is optimal and the iterations are bound for all the grid levels. In Tables 4 and 5 we compare all smoothers and we observe that the Braess-Sarazin smoother leads to faster convergence of the multigrid with fewer iterations and less CPU ahead of other smoothers. The inexact Uzawa did not disappoint in relaxing the error but the DGS and the PMINRES lead to more iterations and computing times. The other observation in Tables 4 and 5 is that the -cycle converges in fewer iterations than the -cycle though it has more computing times for all smoothers.
In Tables 6 and 7 we use different approximations for the preconditioner of the smoothers of the multigrid -cycle and -cycle, respectively. In applying the preconditioners, we approximate the preconditioner̂of the Laplacian stiffness and sparse matrix and̂by a geometric multigrid V(1, 1)cycle method ( ). The multigrid is a well-known fast solver for such systems. The multigrid solver is an inner iteration of the smoothers. The results in Tables 6 and 7 also show that the one iteration of the multigrid V-cycle is a suitable approximation of the smoothers since the multigrid solver has improved iterations from the ones in Tables 4 and 5. In both tables the multigrid method is optimal in solving the indefinite systems and the number of iterations is bounded for all smoothers independent of the mesh size or grid level. Table 8 shows the changes in the estimated a posteriori errors for regularized driven cavity flow using 2 -1 approximation for the flow: using the strategy built in IFISS [3,8] that, for every element error, the local error estimation is given by the combination of the energy norm of the velocity error and the 2 norm of the divergence error; that is,   where e is the velocity error estimate and = ‖∇ ⋅ u‖ and := (∑ ∈ ℎ 2 ) 1/2 are the global error estimator, using different smoothers from one level to the other. From Table 8 we note that the velocity divergence is clearly converging at a faster rate to (ℎ 3 ), which means that the estimated global error is increasingly dominated by the velocity error component as ℎ → 0. Figure 1 shows the sample grid output at the levels 1/64 and 1/128, the sample velocity solution (exponential streamlines), and the pressure plot at the same level with the same smoother.

Conclusion
The purpose of this study was to explore the multigrid solver for the Stokes equations. We have introduced four smoother iterative methods for both multigrids -cycle and -cycle to solve the indefinite systems emanating for the mixed finite element discretization of the Stokes problem. We analyse the construction of the multigrid solver, construction of the smoothers, computation costs, and CPU time as an indicator of the performance of each smoother at all grid levels. Numerical experimental results are given for both -cycle and -cycle for the smoothers at different grid levels.
The Scientific World Journal We have found out that for both cases and for all smoothers used in this study the multigrid solver is optimal and the number of iterations is bounded for all the grid levels. For the steady Stokes equations and the choices of the smoothers used the Braess-Sarazin like smoother became the best iteration to relax the error of the multigrid solver. All the numerical results show that the one V-cycle multigrid iteration is also a suitable preconditioner for the smoothers used.