An Efficient Iterative Method for the Generalized Stokes Problem

Ahmed Samehgif  and Vivek Saringif


This paper presents an efficient iterative scheme for the generalized Stokes problem, which arises frequently in the simulation of time-dependent Navier-Stokes equations for incompressible fluid flow. The general form of the linear system is


where tex2html_wrap_inline271 is an tex2html_wrap_inline273 symmetric positive definite matrix, in which M is the mass matrix, T is the discrete Laplace operator, tex2html_wrap_inline279 and tex2html_wrap_inline281 are positive constants proportional to the inverses of the time-step tex2html_wrap_inline283 and the Reynolds number Re respectively, and B is the discrete gradient operator of size tex2html_wrap_inline287 (k;SPMlt;n). Even though the matrix A is symmetric and positive definite, the system is indefinite due to the incompressibility constraint ( tex2html_wrap_inline293 ). This causes difficulties both for iterative methods and commonly used preconditioners. Moreover, depending on the ratio tex2html_wrap_inline295 , A behaves like the mass matrix M at one extreme and the Laplace operator T at the other, causing problems for the common iterative methods employed to solve this system.

The generalized Stokes problem is one of the most time-consuming steps in the large scale simulation of incompressible fluid flows. Iterative methods, which are indispensable for solving this problem in realistic situations, rely heavily on effective preconditioners that can be efficiently implemented on multiprocessors. Therefore, the issues of efficient iterative algorithms and robust, effective and parallelizable preconditioners for the generalized Stokes problem must be resolved satisfactorily.

Previous efforts to solve the linear system in Eq. 1 can be broadly classified as follows:

  1. Uzawa-type methods. These methods solve the linear system:


    Each iteration requires the solution of the linear system Ax=b. If an iterative method is used to solve Ax=b, then we obtain a two-level solver with inner and outer iterations. The inner system may be solved in many ways including multigrid, spectral methods and preconditioned CG algorithm. Convergence may be considerably improved by replacing A with tex2html_wrap_inline309 or by accelerating the linear system by the CG method. Single-level methods which approximate the generalized Stokes problem by a nearly incompressible problem have also been developed. However, effective preconditioners for this approach often rely on expensive linear system solves using an inner iterative method. The reader is referred to [4, 1, 3, 7, 10, 9] for a variety of Uzawa-type methods.

  2. Projection methods. Projection methods (see [8, 2]) require the solution of the linear system:


    where P is the orthogonal projection onto the null space of tex2html_wrap_inline313 . In each iteration, one computes the action of the projector tex2html_wrap_inline315 on a vector u, which requires solving the system tex2html_wrap_inline319 , where tex2html_wrap_inline321 is analogous to a Laplace operator on the domain. For a large three dimensional problem, computing the projection of u onto tex2html_wrap_inline325 is almost as difficult as solving the system Au=b, and may require the use of an inner iterative method.

  3. Augmentation methods. Augmentation methods proposed in [5, 6] reformulate the system by augmenting the matrix B with a suitable matrix F, and use the CG method to solve the resulting symmetric and positive definite system. Here, one needs to solve systems of the type (B,F)x=b and tex2html_wrap_inline335 . At present however, the choice of F for which (B,F) is easy to invert and the system matrix is well-conditioned, has not yet been determined.

Both Uzawa-type and projection methods are expensive two-level nested iterative methods with rapid convergence assured only for certain extreme values of tex2html_wrap_inline295 . Further, these schemes suffer from the lack of good preconditioners, which have been elusive due to the complicated coefficient matrices arising in these methods.

In this paper, we propose a multilevel scheme for the solution of elliptic problems that has the desirable properties of effective preconditioning and efficient implementation on multiprocessors. We also extend this multilevel scheme to the generalized Stokes problem and present numerical experiments comparing the multilevel scheme with the Uzawa and projection methods.

Multilevel Scheme for the Generalized Stokes Problem A multilevel approach is used to compute a basis for tex2html_wrap_inline325 , where B is the discrete gradient operator. This null-space-basis, tex2html_wrap_inline347 is expressed as the product of a sequence of sparse matrices, s.t. it requires only O(n) operations to perform a matrix-vector product with tex2html_wrap_inline347 . In fact, we determine matrices P, D and Z such that


where P is a nonsingular tex2html_wrap_inline273 matrix, D is a tex2html_wrap_inline365 diagonal matrix and Zis a tex2html_wrap_inline365 orthogonal matrix. Further, tex2html_wrap_inline371 where tex2html_wrap_inline347 comprises of the last n-k columns of P.

Eq. 1 may be rewritten as


where tex2html_wrap_inline379 . The algorithm for computing the solution of the generalized Stokes problem is as follows:


To solve the generalized Stokes problem efficiently, it is imperative that the matrix tex2html_wrap_inline389 , in which tex2html_wrap_inline271 , is well-conditioned for a wide range of the parameters tex2html_wrap_inline279 and tex2html_wrap_inline281 . When tex2html_wrap_inline295 is large, A may be approximated by tex2html_wrap_inline401 , and when tex2html_wrap_inline295 is small, A approaches the Laplace operator tex2html_wrap_inline407 . It can be shown that for large values of tex2html_wrap_inline295 , standard preprocessing techniques like scaling the matrix A with a diagonal matrix obtained by lumping the mass matrix M yields a well-conditioned system. In cases where tex2html_wrap_inline295 is small, we expect tex2html_wrap_inline389 to be well-conditioned due to the ``inverse'' nature of the matrices A and tex2html_wrap_inline421 .

Numerical experiments for the solution of the generalized Stokes problem for a lid-driven cavity problem were performed on an IBM RS6000 with 66.5 MHz clock and 256 MB memory. We consider a two-dimensional unit square domain with tex2html_wrap_inline423 on the boundary except the top edge where tex2html_wrap_inline425 , and p=0 at the bottom-left corner. The domain is uniformly discretized with the Marker-and-Cell (MAC) scheme. Figs. 1 and  2 presents the number of iterations and time required to solve the linear system in Eq. 5. From these results, it may be concluded that in the worst case, the condition number of the system matrix tex2html_wrap_inline389 is tex2html_wrap_inline431 , which indicates effective preconditioning of the system matrix.

The number of iterations and time required by the outer iterative method in the Uzawa method are presented in Figs. 3 and  4. It may be noted that in general, each iteration of these methods is relatively expensive due to system solves.

Figure 1: Multilevel scheme for the solution of the generalized Stokes problem on a uniform tex2html_wrap_inline249 mesh for tex2html_wrap_inline251 .

Figure 2: Multilevel scheme for the solution of the generalized Stokes problem on a uniform tex2html_wrap_inline249 mesh for tex2html_wrap_inline251 .

Figure 3: Uzawa methods for the solution of the generalized Stokes problem on a uniform tex2html_wrap_inline249 mesh for tex2html_wrap_inline251 .

Figure 4: Uzawa methods for the solution of the generalized Stokes problem on a uniform tex2html_wrap_inline249 mesh for tex2html_wrap_inline251 .

These results clearly indicate that the multilevel scheme is a robust and effective preconditioning strategy for the solution of the linear system of the generalized Stokes problem for a wide choice of tex2html_wrap_inline295 . This makes the multilevel scheme an attractive approach especially since convergence is relatively independent of tex2html_wrap_inline451 , and therefore doesn't constrain the choice of the time-step tex2html_wrap_inline283 . Furthermore, it is a single-level scheme with O(n) computation per iteration, unlike the two-level Uzawa-type and projection methods which require expensive linear system solves in each iteration. It must also be pointed out that the effort required in generating the basis tex2html_wrap_inline347 can be amortized over a number of system solves since B is invariant over numerous time steps for the nonlinear Navier-Stokes equations. Current work includes numerical experiments with systems obtained using the finite elements method on unstructured grids for two and three-dimensional domains. We expect to present comparative results of numerical experiments for these problems.


R.E. Bank, B.D. Welfert, and H. Yserentant. A class of iterative methods for solving saddle point problems. Numer. Math., 56:645-666, 1990.

R. Bramley. An orthogonal projection algorithm for generalized Stokes problems. Technical Report 1190, CSRD, Univ. of Illinois Urbana-Champaign, 1992.

N. Dyn and Jr. W.E. Ferguson. The numerical solution of equality-constrained quadratic programming problems. Math. Comp., 41(163):165-170, 1983.

M. Fortin and R. Glowinski. Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems. North-Holland, New York, 1983.

F. G. Lou and A. Sameh. An expansion method for solving saddle-point problems. Technical Report 92-083, AHPCRC, Univ. of Minnesota, 1992.

F. G. Lou, A. Sameh, and V. Sarin. An augmentation method for solving saddle-point problems. Technical Report 95/29, MSI, Univ. of Minnesota, 1995.

T. Rusten and R. Winther. A preconditioned iterative method for saddle-point problems. SIAM J. Matrix Anal. Appl., 13(3):887-904, 1992.

A. Sameh and J. A. Wisniewski. A trace minimization algorithm for the generalized eigenvalue problem. SIAM J. Numer. Anal., 19(6):1243-1259, 1982.

D. Silvester and A. Wathen. Fast iterative solution of stabilised Stokes systems Part 2: Using general block preconditioners. SIAM J. Numer. Anal., 31(5):1352-1367, 1994.

A. Wathen and D. Silvester. Fast iterative solution of stabilised Stokes systems Part 1: Using simple diagonal preconditioners. SIAM J. Numer. Anal., 30(3):630-649, 1993.

This research has been supported in part by the NSF under the grant NSF/CDA 9396332-001.

Dept. of Computer Science, Univ. of Minnesota, Twin Cities

Dept. of Computer Science, Univ. of Illinois, Urbana-Champaign

Vivek Sarin
Tue Sep 22 08:12:33 EST 1998

| AEM Home | Institute of Technology |

| Academics | Research | People | Information | Contact AEM |
| Major Projects | by Specialty | Other Programs |

Solid-Liquid Flows Home | Publications

Last updated October 16, 2000