Research Article  Open Access
The Mixed Finite Element Multigrid Method for Stokes Equations
Abstract
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 BraessSarazin type smoothers. A comparative study of the smoothers shows that the BraessSarazin smoothers enhance good performance of the multigrid method. We study the problem in a twodimensional domain using stable HoodTaylor 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.
1. 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: where is a velocity field, represents pressure, and is an external force field. The problem is considered with (1)â€“(3) defined on the domain 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 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 infsup stable. In our experiment we use HoodTaylor 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 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â€“6] and elliptic partial differential equations [7â€“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â€“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], Vankatype [9, 12, 13, 15], BraessSarazintype ([13, 16â€“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: BraessSarazin, 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.
2. 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 and the second equation (2) with , subsequently integrating over the domain , applying the Gauss theorem, and incorporating the boundary condition (3), we obtain the variational form.
Find and such thatwhere and are continuous bilinear forms defined aswhere represents a componentwise scalar product that is and and . The wellposedness follows from the coercivity of in the LaxMilgram theorem [19, 20] and partly from the infsup condition [7, 8, 16, 21â€“23]. Below is a sketch of the analysis of the existence uniqueness and stability of the solution of mixed problem (5):(i)the bilinear form is bounded or continuous if(ii)the bilinear form is coercive on ; that is, there exists a positive constant : (iii)the bilinear form is symmetric and nonnegative if (iv)the bilinear form is bounded if (v)the bilinear form satisfies the infsup condition; that is, there exists a constant :For instance, in [22], it is shown that in our concrete case fulfills the infsup condition; thus we can combine (i)â€“(v) to give the following theorem.
Theorem 1. The variational problem (5) is uniquely solvable provided properties (i)â€“(v) are all satisfied.
The proof relies on the closed range theorem and on the LaxMilgram theorem. The details can be found in [22, 24].
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â€“21, 25].
We assume that . We define the finite dimensional spaces. Let and be subspaces of and , respectively.
Now we can formulate a discrete version of problem (5).
Find a couple such that The finite element discretization should satisfy the infsup 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 independent ellipticity constant) and that there exists a constant (independent of ) such that the discrete infsup condition holds. Then the associated (discretized, steady state) Stokes problem has a unique solution , and there exists a constant such that
If the basis of is given by and of is given by , thenwhere is the number of inner nodes and is the number of boundary nodes so the coefficients 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 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 HoodsTaylor square finite elements which are known to be stable.
We specifyAn element is uniquely determined by specifying components of 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 byThe entries of the right hand side vector areThe linear algebraic system can be represented aswhere , , and .
The solution vectors from (15) are the mixed finite element solution. The system (17)â€“(19) is called the discrete Stokes problem.
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].
3. 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 be a sequence of subspaces of the finite dimensional subspaces and defined on sequence of grids with mesh sizes with . We define a hierarchy/family of nested finite element subspaces for the velocity and pressure:where subspace which corresponds to is the refinement of with subspace such that . At the discrete level with the defined discrete spaces and bases, the linear algebraic system is defined by where , , and .
The main goal is to find the pair 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 , where(i) is the output of velocity and pressure after one step of the multigrid algorithm at level ;(ii) is the input velocity at level ;(iii) is the input pressure at level ;(iv) and are restriction operators for velocity and pressure, respectively, from level to level ;(v) and are prolongation operators for velocity and pressure, respectively, from level to level .
Algorithm 3 (multigrid algorithm). if (coarsest grid) else define â€‰(1)â€‰PreSmoothing: Smoothing operator starting with, with smoothing steps, producing ,(a)â€‰defect/residual (2)â€‰restrict the defect(3)â€‰approximate solution(4)â€‰applying one/two iterations of at the recursive call:(a)â€‰apply steps of (b)â€‰Set (c)â€‰Set (d)â€‰computeâ€‰for â€‰end(5)â€‰Correction Step define the new iterate:
Postsmoothing. Starting with perform smoothing steps using smoothing operator to produce 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 represents such methods; in this study it represents BraessSarazin, 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 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 is defined bywhere the restriction operators and for velocity and pressure, respectively. The prolongation between levels and is defined again as where the prolongation operators and are representations of the following relations for the quadratic interpolation of the velocity () and for the linear interpolation of the pressure ().(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 . 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.
3.1. The Smoothers
The most crucial part is the proper choice of a smoothing technique. Usually, the wellknown smoothing iterations for the scalar problems (damped Jacobi or GaussSeidel 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 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 Vankatype 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 BraessSarazintype 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.
3.1.1. BraessSarazinType Smoother
The BraessSarazin smoothers proposed in [17] and used in [13, 18] solve a large saddle point problem in each smoothing step. This BraessSarazin or SIMPLEtype 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:with and given. The smoothing BraessSarazin iteration (35) solves the auxiliary problem with and . Inherent in the system system (36) is the problem of the auxiliary pressure variable This system is solved approximatively by an iterative solver. From the system we get approximately which can be used to approximately determine from
3.1.2. Inexact Uzawa Type Smoothers
The variant of the inexact Uzawa iteration used as a smoother is outlined.
Algorithm 4. (1) For : smoothing steps.
(2) Compute the residual .
(3) Compute the residual .
(4) Solve .
(5) Solve .
(6) Solve .
(7) Update the velocity and pressure End.
Step (6) in the outline may be rearranged as with obtained as a byproduct of step (5). This saves the application of 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â€“28].
3.1.3. 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 (22) and with . We solve the transformed residual equationwith and being invertible approximations of and , respectively. A single iteration with the update through a distributive matrix is performed by the following algorithm.
Algorithm 5 (). (1)â€‰ â€‰Smooth momentum equations (2) Smooth the transformed continuity equation (3) Transform the correction back to the original variablesThe 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.
3.1.4. 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 with . The block preconditioning requires the solution of two systems of equations with matrices and at each MINRES iteration. If 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 , the mass matrix can be further replaced by its diagonal matrix [8].
3.2. 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 and are the number of pre and 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 BraessSarazin smoother with being the iteration matrix of the smoother (34) and the 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 twogrid method, with and the pre and postsmoothing steps, respectively, applied to (22) results in the iteration matrixThe 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 define the following norms, Euclidean norm by , and on the following norm is applied:Furthermore we introduceUsing the norms defined above and taking and above we obtainThe 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 regular. Let be the coefficient stiffness matrix and and the prolongation and the restriction operators. Then there exists a constant independent of and using scaling induced by thenwhere .
The smoothing property is dependent on the smoother used. It varies from one smoother to another. In this work we used the BraessSarazin in which we solve the system (37) exactly and sufficiently accurate inexact inner solver.
Theorem 7 (smoothing property). Let be the coefficient stiffness matrix and the smoothing operator . Then where for and is a decreasing function with .
Combining the approximation property Theorem 6 with the smoothing property Theorem 7 produces a twogrid convergence result.
Theorem 8. Assume that and that is such that the problem (5) is regular. Then for the twogrid method the following holds:with a constant independent of l and m.
Using this twogrid contraction number bound the multigrid cycle method convergence results can be derived using ideas in [1, 2].
4. 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 [11]. The Dirichlet noslip boundary condition is applied on the side and bottom boundaries. The mixed finite element method was used to discretize the cavity domain .
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 preconditioners 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 and one cycle is used to approximate Schur compliment matrix . The second case is to use one cycle for both evaluations of and . The method becomes IUZAWAMG.(iii)BraessSarazin smoother (BS): the two cases are conspired for the evaluation of the preconditioners. Firstly, the approximation and one cycle is used to solve the approximate Schur compliment matrix . The second case is to use one cycle for both evaluations of and . The method becomes BSMG.(iv)PMINRES smoother: the first case is to use diagonal preconditioner for and and the second case is one cycle for committing the inversion of the Laplacian operator for velocity as one 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 HoodTaylor 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 programming language and the MATLAB builtin Minres functions are used for the smoother. For the discretization we start with a uniform square grid with 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 and being presmoothing and postsmoothing steps, respectively. The smoothers are determined by specifying approximations and as highlighted in (i)â€“(vi) and in all cases where one cycle inner multigrid iteration is used with and 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 , where is satisfied. The schemes converge if the stopping criteria are satisfied. The results show the first case of the evaluation, the preconditioners of the smoothers and with , and for all cases evaluation of by one cycle inner multigrid iteration with and being Gauss Seidel iteration of the presmoothing and postsmoothing steps, respectively.
Tables 2 and 3 show the number of iterations and computing time to demonstrate the effects of different cycle (1, 2, 3, and 4) and cycle (1, 2, 3, and 4) being pre and postsmoothing steps with BraessSarazin (BS) smoother . We compare the performance of the cycle and cycle multigrid iterations with various smoothing steps at different grid levels using one of the smoothers, BraessSarazin.


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 as and for BraessSarazin, 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 BraessSarazin 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 cycle method (). The multigrid is a wellknown 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 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 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 norm of the divergence error; that is, where is the velocity error estimate and and are the global error estimator, using different smoothers from one level to the other.
