9512.net
甜梦文库
当前位置:首页 >> >>

An augmented Lagrangian-based approach to the Oseen problem


SIAM J. SCI. COMPUT. Vol. 28, No. 6, pp. 2095–2113

c 2006 Society for Industrial and Applied Mathematics

AN AUGMENTED LAGRANGIAN-BASED APPROACH TO THE OSEEN PROBLEM?
MICHELE BENZI? AND MAXIM A. OLSHANSKII? Abstract. We describe an e?ective solver for the discrete Oseen problem based on an augmented Lagrangian formulation of the corresponding saddle point system. The proposed method is a block triangular preconditioner used with a Krylov subspace iteration like BiCGStab. The crucial ingredient is a novel multigrid approach for the (1,1) block, which extends a technique introduced by Sch¨berl for elasticity problems to nonsymmetric problems. Our analysis indicates that this approach o results in fast convergence, independent of the mesh size and largely insensitive to the viscosity. We present experimental evidence for both isoP2-P0 and isoP2-P1 ?nite elements in support of our conclusions. We also show results of a comparison with two state-of-the-art preconditioners, showing the competitiveness of our approach. Key words. Navier–Stokes equations, ?nite element, iterative methods, multigrid, preconditioning AMS subject classi?cations. 65F10, 65N22, 65F50 DOI. 10.1137/050646421

1. Introduction. We consider the numerical solution of the steady Navier– Stokes equations governing the ?ow of a Newtonian, incompressible viscous ?uid. Let Ω ? Rd (d = 2, 3) be a bounded, connected domain with a piecewise smooth boundary ?Ω. Given a force ?eld f : Ω → Rd and boundary data g : ?Ω → Rd , the problem is to ?nd a velocity ?eld u : Ω → Rd and a pressure ?eld p : Ω → R such that (1.1) (1.2) (1.3) ?νΔu + (u · ?) u + ?p = f in Ω,

div u = 0 in Ω, u = g on ?Ω,

where ν > 0 is the kinematic viscosity coe?cient (inversely proportional to the Reynolds number Re), Δ is the Laplace operator in Rd , ? denotes the gradient, and div is the divergence. To determine p uniquely we may impose some additional condition, such as (p, 1) =
Ω

p dx = 0.

Equation (1.1) represents conservation of momentum, while (1.2) represents the incompressibility condition, or mass conservation. Due to the presence of the convective term (u·?) u in the momentum equations, the Navier–Stokes system is nonlinear. It can be linearized in various ways. A widely used linearization scheme is the one
? Received by the editors November 30, 2005; accepted for publication (in revised form) May 30, 2006; published electronically December 5, 2006. http://www.siam.org/journals/sisc/28-6/64642.html ? Department of Mathematics and Computer Science, Emory University, Atlanta, GA 30322 (benzi@mathcs.emory.edu). The work of this author was supported in part by the National Science Foundation grant DMS-0511336. ? Department of Mechanics and Mathematics, Moscow State M. V. Lomonosov University, Moscow 119899, Russia (Maxim.Olshanskii@mtu-net.ru). The work of this author was supported in part by the Russian Foundation for Basic Research and the Netherlands Organization for Scienti?c Research grants NWO-RFBR 047.016.008, RFBR 05-01-00864, RFBR-DFG 06-01-04000.

2095

2096

MICHELE BENZI AND MAXIM A. OLSHANSKII

based on Picard’s iteration; see, e.g., [10, section 7.2.2]. Starting with an initial guess u(0) (with div u(0) = 0) for the velocity ?eld, Picard’s iteration constructs a sequence of approximate solutions (u(k) , p(k) ) by solving the linear Oseen problem: (1.4) (1.5) (1.6) ?νΔu(k) + (u(k?1) · ?) u(k) + ?p(k) = ? f div u(k) = 0 u
(k)

in Ω, in Ω, on ?Ω

=g

(k = 1, 2, . . . ). Note that no initial pressure needs to be speci?ed. Under certain conditions on ν (which should not be too small) and ? (which should not be too f large in an appropriate norm), the steady Navier–Stokes equations (1.1)–(1.3) have a unique solution (u? , p? ), and the iterates (u(k) , p(k) ) converge to it as k → ∞ for any choice of the initial velocity u(0) ; see [10, Chapter 7] and the references therein. Hence, at each Picard iteration one needs to solve an Oseen problem of the form (1.7) (1.8) (1.9) ?νΔu + (w · ?) u + ?p = ? f div u = 0 u=g in Ω, in Ω, on ?Ω

with a known, divergence-free coe?cient w. Discretization of (1.7)–(1.9) using LBBstable ?nite elements [10, 26] results in a linear system of the form (1.10) A B BT O u p = f 0 , or A x = b,

in which u represents the discrete velocities and p the discrete pressure. Here A = diag (A1 , . . . , Ad ) is a block diagonal matrix, where each block corresponds to a discrete convection-di?usion operator with the appropriate boundary conditions. Note that A is nonsymmetric. The rectangular matrix B T represents the discrete gradient operator, while B represents its adjoint, the (negative) divergence operator. The solvability of this system is discussed, e.g., in [2, section 3.2]; here we assume that the coe?cient matrix in (1.10) is nonsingular. Linear systems of the form (1.10) are often referred to as generalized saddle point systems. In recent years, a great deal of e?ort has been invested in solving systems of this form. Most of the work has been aimed at developing e?ective preconditioning techniques; see [10], and [2] for an extensive survey. In spite of these e?orts, there is still considerable interest in preconditioning techniques that are truly robust, i.e., techniques which result in convergence rates that are largely independent of problem parameters such as mesh size and viscosity. In this paper we describe a promising approach based on an augmented Lagrangian (AL) formulation. The success of this method crucially depends on the availability of a robust multigrid solver for the (1,1) block (submatrix) in the AL formulation of the saddle point system; we develop such a method by building on previous work by Sch¨berl [29], together with appropriate o smoothers for convection-dominated ?ows. Rather than as a solver, this multigrid iteration will be used to de?ne a block preconditioner for the outer iteration on the coupled saddle point system. We will show that this approach is especially appropriate for discretizations based on discontinuous pressure approximations, but can be used to construct preconditioners for other discretizations using continuous pressures. As an example of the ?nite element (FE) method with discontinuous pressures we will use the isoP2-P0 pair. In this case, our numerical experiments demonstrate a robust

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2097

behavior of the solver with respect to h and ν for some typical wind vector functions w in (1.7). Further, the isoP2-P1 FE pair is used for the continuous pressure-based approximation. For this case numerical experiments show h-independent convergence rates, with mild dependence on ν when the viscosity becomes very small. The remainder of the paper is organized as follows. In section 2 we present the (standard) abstract formulation of the Oseen problem and discuss the FE approximations. In section 3 we describe the AL approach and the resulting block preconditioner and brie?y discuss related approaches that can be found in the literature. Section 4 is devoted to a study of the spectrum of the preconditioned system, assuming exact solves for the (1,1) block. The multigrid solver for the (1,1) block is described in section 5. Numerical experiments on some two-dimensional (2D) model problems are presented in section 6. Here we also present a comparison with two of the best available preconditioning techniques, showing that our method is quite competitive in terms of convergence rates, robustness, and e?ciency. Some concluding remarks are given in section 7. 2. Finite element method. For simplicity, we assume the boundary conditions in (1.9) to be Dirichlet homogeneous (i.e., g = 0). The weak formulation of the Oseen problem reads as follows: Given f ∈ H?1 (Ω), ?nd u ∈ H1 (Ω) and p ∈ L2 (Ω) such 0 0 that ? v ∈ H1 (Ω), q ∈ L2 (Ω), 0 0 L(u, p; v, q) := ν(?u, ?v) + ((w · ?) u, v) ? (p, div v) + (q, div u). Given conforming FE spaces Vh ? H1 (Ω) and Qh ? L2 (Ω), the Galerkin FE dis0 0 cretization of (1.7)–(1.9) is based on the following weak formulation: Find {uh , ph } ∈ Vh × Qh such that (2.1) L(uh , ph ; vh , qh ) = f , vh ? vh ∈ Vh , qh ∈ Qh . L(u, p; v, q) = f , v

There are several critical issues associated with the Galerkin FE method for the Oseen problem. One is the compatibility of Vh and Qh , i.e., the validation of the LBB (infsup) stability condition, which guarantees that the FE velocity space is “rich enough” relative to the FE pressure space. This ensures well-posedness and full approximation order for the FE linear problem. For the numerical experiments in the paper we used the isoP2-P0 FE pair (piecewise constant pressure and piecewise linear continuous velocity on two-times ?ner triangulation) and the isoP2-P1 FE pair (piecewise linear pressure and piecewise linear continuous velocity on two-times ?ner triangulation). Both pairs are LBB stable: There exists a mesh-independent constant μ(Ω) such that (2.2) inf sup (qh , div vh ) = μh ≥ μ(Ω) > 0. ?vh qh

qh ∈Qh vh ∈Vh

We note that the chosen FE pairs have stable higher order variants P2-P0 and P2-P1 (the latter is also known as Taylor–Hood element). All considerations and conclusions in the paper remain valid for these higher order ?nite elements. Piecewise linear velocity-based elements were used solely for the ease of implementation. Another potential source of instabilities in (2.1) is the presence of dominating convection terms. This necessitates stabilization of the discrete system if the mesh is not su?ciently ?ne to resolve the sti? behavior of the system. We outline below the SUPG-type method we used. Many more details on the family of SUPG stabilization

2098

MICHELE BENZI AND MAXIM A. OLSHANSKII

methods can be found in, e.g., [6, 28, 31]. Using (2.1) as the starting point, a weighted residual for the FE solution multiplied by a solution-dependent test function is added: (2.3) L(uh , ph ; vh , qh ) +
τ ∈Th

στ (?νΔuh + w·?uh + ?ph ? f , w·?vh )τ = (f , vh ) ? vh ∈ Vh , qh ∈ Qh .

The “new” term in (2.3) is evaluated elementwise for each element τ of a triangulation Th . The parameters στ are element- and problem-dependent. The general idea behind the choice of στ is to add almost no additional stabilization terms in regions of small mesh Reynolds numbers, so as to recover optimality of the Galerkin method, but to add them in regions of large mesh Reynolds numbers. Several recipes for the particular choice of the stabilization parameters can be found in the literature. We use the one from [22]: (2.4) στ = σ ? hτ w 2Reτ , 1 + Reτ L2 (τ ) Reτ := w
L2 (τ ) hτ

ν

.

This formula for the parameter with σ ∈ [0.2, 1] was successfully used in many nu? merical experiments in [31]. In our experiments we take σ = 0.3. ? Remark 2.1. The particular choice of stabilization parameters στ can be quite important for the discrete solution accuracy. On the other hand, using high order ?nite elements [20] or su?ciently ?ne meshes may dispense with the need to stabilize the problem. We found, however, that in order to have an e?cient solver it is still important to add stabilization to the coarse grid problems in the multigrid method (see details in section 5). Additionally, the e?ciency of the iterative solver described below was found to be essentially independent of the particular choice of σ = O(1) or of ? di?erent formulae for στ known from the literature for convection-di?usion problems (see, e.g., [11]). We can also assume that any pressure-dependent part in the stabilization term enters the nonlinear residual in the Picard iterations, thus preserving the saddle point structure (1.10) of the problem. 3. Augmented Lagrangian formulation. Some of the most e?ective solvers for the Oseen problem found in the literature are based on block triangular preconditioners of the form (3.1) P= ? A O BT ? S ,

? ? where A and S are approximations to A and to the Schur complement S = ?BA?1 B T , ? respectively; see, e.g., [2, 10]. The approximation A often consists of an inexact solver for linear systems with coe?cient matrix A, such as one or more iterations of an appropriate multigrid solver for convection-di?usion equations. The search for ? good approximations S to the Schur complement (or its inverse) has generated much research in recent years [10]. Note that S = ?BA?1 B T is a dense matrix. When ν is relatively large (moderate Reynolds numbers), then S is well conditioned and can be well approximated by the pressure mass matrix Mp , provided that the discretization is LBB-stable. When ν is small, however, ?nding good and inexpensive approximations to S is not easy. While the best available methods exhibit h-independent convergence of the preconditioned iterations [10], some degradation in the rate of convergence is

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2099

observed as ν → 0. For large Reynolds numbers, the number of iterations can be high. Here we follow an approach that bypasses the need for a good approximation to the dense matrix S. Suppose we replace the original problem (1.10) with the “regularized” one, (3.2) A B BT 1 ?γ W u p = f 0 ,

where γ > 0 is a parameter and W is a positive de?nite matrix. The solution [u(γ), p(γ)] of (3.2) tends to the solution of (1.10) for γ → ∞. Observing that the Schur complement of the (1,1) block in the coe?cient matrix of (3.2) (sometimes referred to as the “primal” or “velocity” Schur complement) is A + γB T W ?1 B, we expect (3.3) P(γ) = A + γB T W ?1 B O BT 1 ?γ W

(or, in practice, some inexact version of it) to be an e?ective preconditioner for the regularized problem (3.2). Note that if W is diagonal, or block diagonal with small blocks, the matrix A + γB T W ?1 B is also going to be sparse. Our choice of W is discussed in the next section. There is, however, a well-known di?culty with this approach. The matrix in the regularized problem is close to the original one only for γ large, and for large γ the (1,1) block A + γB T W ?1 B becomes increasingly ill conditioned (note that B T B has a null space of dimension n ? m for n = dim Vh , m = dim Qh ). Hence, ?nding an e?ective inexact solver for the (1,1) block becomes di?cult. It is therefore preferable to keep the value of γ moderate, say, γ = O(1) or, for reasons of scaling, γ ≈ w . This is the idea behind the classical AL approach. In its simplest form, this method consists of replacing the original system (1.10) with the equivalent (3.4) A + γB T W ?1 B B BT O u p = f 0 , or A x = b.

Clearly, this system has precisely the same solution as the original one (1.10). We refer to [13] for an extensive discussion of AL techniques and their applications. We propose to precondition the AL system (3.4) with the block triangular preconditioner (3.5) P= ? Aγ O BT ? S ,

? where the action of A?1 is computed via an appropriate (inexact) solver for linear γ systems involving the matrix A + γB T W ?1 B, (3.6) ? Aγ ≈ A + γB T W ?1 B,

? and S is implicitly de?ned through its inverse, (3.7) ? ?1 ? S ?1 := ?ν Mp ? γW ?1 .

? ?1 Here Mp denotes an approximate solve with the pressure mass matrix. This choice ? is motivated by Lemma 4.1 below and will be analyzed in the next section. The of S

2100

MICHELE BENZI AND MAXIM A. OLSHANSKII

pressure mass matrix Mp is diagonal for FE pairs based on piecewise constant pressure approximations. For higher order pressure approximations Mp is no longer diagonal ? but it is spectrally equivalent to a diagonal matrix Mp ; see [34]. Thus, linear systems with Mp can be solved easily, either exactly or approximately. It follows from the identity (3.8) P ?1 = ?γ A?1 O O Im In O BT ?Im In O O ? ?S ?1

that the action of the preconditioner on a given vector requires one application of ? ? A?1 , one of S ?1 , and one sparse matrix-vector product with B T . The question γ of the adequacy of the augmented Schur complement preconditioner (3.7) will be ? addressed in sections 4 and 6. Besides this, the major issue is the de?nition of A?1 , γ the approximate (1,1) block solver. Developing such a solver is a nontrivial task. Note that the introduction of the additional term γB T W ?1 B in the (1,1) block of the saddle point matrix introduces a coupling between the components of the velocity vector, which is originally avoided in the Oseen linearization compared to the Newton one. For vanishing convection the (1,1) block essentially reduces to the one for the classical elasticity problem. For moderate and high Reynolds numbers, however, it becomes increasingly nonsymmetric. It is nevertheless possible to develop e?cient iterative solvers for the (1,1) block, at least for some choices of the FE discretization. In section 5 we describe an e?ective multigrid method for the (1,1) block. We conclude this section with a brief discussion of related work. An AL approach to the Oseen problem can already be found in [13, Chapter 2], in the context of the classical Uzawa’s method. However, the crucial aspect of how to e?ciently solve linear systems involving A + γB T W ?1 B was not discussed there. Among more recent papers, we mention [7], where a preconditioner based on (3.2) is applied to (symmetric) linear incompressible elasticity problems. The inexact solution of systems involving (3.2) is accomplished by a “primal” Schur complement approach, leading to an overall preconditioning strategy which resembles the approach taken here. There are, however, important di?erences. Our preconditioner is applied to the AL formulation (3.4), whereas the authors of [7] work with the original saddle point formulation. A more important di?erence is the fact that the problem tackled in [7] is symmetric; thus, the authors are able to leverage existing techniques for linear elasticity problems when developing an e?cient solver for the (1,1) block. We mention that a related approach has been used in [17] to construct block preconditioners for saddle point systems arising from the curl-curl formulation of the time-harmonic Maxwell’s equations, also a symmetric problem. The paper [14] applies a preconditioning technique similar to the one in [7] to Stokes and linearized Navier–Stokes problems. The numerical evidence reported in [14] indicates that as the Reynolds number increases, the rate of convergence of the inner iterative solver used for the (1,1) block (a two-level overlapping additive Schwarz preconditioner with a multiplicative coarse grid correction) tends to deteriorate as the mesh is re?ned. The deterioration of this solver is already noticeable for a 2D driven cavity problem with Reynolds number Re = 400, a fairly easy problem (close to Stokes ?ow). Furthermore, no analysis of the preconditioner is given. On the positive side, the solver in [14] appears to be well suited for massively parallel architectures, an aspect which we do not address in this paper. Finally, adding the product γ(div uh , div vh ) on the left-hand side of the FE method formulation (2.1) is well known in the literature on stabilized FE methods

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2101

for ?uid problems and sometimes called “the ? div-stabilization.” Using this stabilization may change the FE solution. However, the e?ect on the iterative solution of (1.10) is similar to the one resulting from adding γB T W ?1 B (with W = Mp ) to the (1,1) block of the linear algebraic system. This linear algebra aspect of the ? div-stabilization was studied for the Stokes problem in [24] and for the linearized rotation form of the Navier–Stokes problem in [25]. In particular, numerical experiments in [25] demonstrate ν-independent convergence of the preconditioned BiCGstab method if γ = O(1) and if exact solves are used for systems involving the (1,1) block. However, in [24] and [25] the question of ?nding a robust solver or preconditioner for the (1,1) block with respect to the ratio ν/γ was left open. 4. Analysis of the preconditioner. It is well known that characterizing the rate of convergence of nonsymmetric preconditioned iterations can be a di?cult task. In particular, eigenvalue information alone may not be su?cient to give meaningful estimates of the convergence rate of a method like preconditioned GMRES [16]. The situation is even more complicated for a method like BiCGStab, for which virtually no convergence theory exists. Nevertheless, experience shows that for many linear systems arising in practice, a well-clustered spectrum (away from zero) usually results in rapid convergence of the preconditioned iteration. In this section we show that for the preconditioner (3.5), the eigenvalues of the preconditioned matrix M = P ?1 A are enclosed in a rectangular region contained in the right half-plane (z) > 0, the sides of which do not depend on the mesh size h and for su?ciently large γ do not depend on ν. Here we assume exact solves for the ? (1,1) block, i.e., Aγ = A + γB T W ?1 B. Hence, at least for this “ideal” version of the preconditioner, the eigenvalues of the preconditioned matrix are bounded away from zero, independently of h and ν. Our analysis makes use of the following simple lemma, which is a straightforward consequence of [12, Exercise 12.12] or [15, Proposition 2.1]. Lemma 4.1. Let A ∈ Rn×n and B ∈ Rm×n (m ≤ n). Let γ ∈ R, and assume that matrices A, A + γB T W ?1 B, BA?1 B T , and B(A + γB T W ?1 B)?1 B T are all invertible. Then (4.1) B(A + γB T W ?1 B)?1 B T
?1

= BA?1 B T

?1

+ γW ?1 .

We note that the conditions of Lemma 4.1 are satis?ed if we assume that B has full row rank and that A is positive real (i.e., the symmetric part of A is positive de?nite); see [2, section 3]. These conditions are satis?ed, in particular, for the matrices considered in this paper. Recall that P ?1 A and A P ?1 have the same eigenvalues. A simple calculation shows that (4.2) A P ?1 = In B(A + γB T W ?1 B)?1 O ? ?B(A + γB W ?1 B)?1 B T S ?1
T

.

Hence, the preconditioned matrix has the eigenvalue 1 of multiplicity n, with the remaining eigenvalues λi (1 ≤ i ≤ m) being solutions of the generalized eigenproblem ? B(A + γB T W ?1 B)?1 B T p = ?λS p. ? Setting W = Mp and Mp = Mp , Lemma 4.1 implies (ν + γ)λ?1 = γ + μ?1 , i i

2102

MICHELE BENZI AND MAXIM A. OLSHANSKII

where μi satis?es the generalized eigenproblem (4.3) BA?1 B T q = μMp q.

Hence, the nonunit eigenvalues of P ?1 A are given by (4.4) λi = γ+ν , γ + μ?1 i 1 ≤ i ≤ m.

It has been shown in [8] that the eigenvalues of (4.3) are enclosed in a rectangle contained in the half-plane (z) > 0; using this result and the relation (4.4), we can conclude that the same is true of the eigenvalues of P ?1 A. If we denote by ai and bi the real and imaginary parts of μi , respectively, easy manipulations result in the following expressions for the real and the imaginary parts of λi : (4.5) (λi ) = (ν + γ)(ai + γ(a2 + b2 )) i i (γai + 1)2 + (γbi )2 and (λi ) = (ν + γ)bi . (γai + 1)2 + (γbi )2

The following result is an immediate consequence of (4.5). Theorem 4.2. Assume W = Mp . The preconditioned matrix P ?1 A has the eigenvalue 1 of multiplicity n. The remaining m eigenvalues λi are given by (4.4), where μi = ai + i bi satis?es (4.3). The following estimates hold: 1? 1 ≤ (λi ) < 1 + νγ ?1 γai + 1 and | (λi )| < (1 + νγ ?1 ) min γbi , 1 γbi .

In particular, all the eigenvalues tend to 1 for γ → ∞. The bounds in Theorem 4.2, while independent of h, will generally depend on ν in an implicit way through ai and bi . For the case of the Galerkin FE method (2.1) (with no stabilization), the results in [8] provide us with the following estimates for the eigenvalues of (4.3): (4.6) c ν ≤ ai ≤ C ν ?1 , |bi | ≤ C1 ν ?1

for some h- and ν-independent positive constants c, C, C1 (here we assume 0 < ν ≤ 1). Together with the estimates from Theorem 4.2, this result implies that it is su?cient to set γ = O(ν ?1 ) to ensure that all nonunit eigenvalues of P ?1 A are contained in a box [a, A] × [?b, b], a > 0, in the complex plane, with a, A, b independent of ν and h. We note that this theoretical result is of somewhat limited use, since in practice we found that it is unnecessary to use large values of γ (!), or to solve linear systems involving the (1,1) block “exactly” (the latter being perhaps less surprising). Remark 4.3. In numerical experiments we found that using γ ≈ w already provides ν- and h-independent convergence of the Krylov subspace method for a wide range of values of ν and h (it should be noted that we also use streamline di?usion? type velocity stabilization). Furthermore, assume that A?1 is not a solver but a γ preconditioner for the (1,1) block and that its e?ciency deteriorates with ν/γ → 0. Then the optimal value of γ in the Krylov subspace method will be quite modest for ν 1 (see Table 6.2 in section 6). ? Remark 4.4. The result of the theorem remains valid if one considers W = Mp ? and with Mp replacing Mp in (4.3). Following the arguments in [8], it can be easily veri?ed that (4.6) still holds. In general, W is not necessarily an approximation to a mass matrix or identity; see the discussion in [15] in conjunction with other saddle point problems. In this paper, however, we do not pursue other choices of W .

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2103

5. A multigrid method for the (1,1) block. In this section we develop a multigrid method to (approximately) solve the system of equations for the discrete velocity at each outer iteration. For the purpose of analysis in this section we consider ? W = Mp ; however, in the implementation it is more practical to set W = Mp . By a simple scaling argument one can assume w ∞ = 1; thus we can let γ = 1. The system one needs to solve is (5.1)
?1 A u + B T Mp B u = f.

To de?ne an FE counterpart of (5.1) we introduce the orthogonal projector Ph : L2 → Qh . The corresponding FE problem takes the form 0 (5.2) ν(?uh , ?vh ) +
τ ∈Th

στ (w · ?uh , w · ?vh )τ + ((w · ?) uh , vh ) + (Ph div uh , div vh ) = (f, vh ) ? vh ∈ Vh .

An ideal multigrid method for solving (5.1) should be robust with respect to ν and di?erent types of wind w and have optimal computational complexity. There are several di?culties in developing such a method. First, note that matrix B has a kernel of dimension n ? m. In the case of ν 1 and if w is small in part of Ω, some vectors ?1 from ker(B) may nearly vanish after multiplication by A+B T Mp B; see Figure 5.1 for an example of locally supported functions from ker(B) for a particular FE pair. Thus, common relaxation iterations like Gauss–Seidel could fail to smooth out oscillatory components from ker(B) that may be present in the error. At the same time, if a local Reynolds number is large, the problem has convection-dominated character at least on the kernel of B. Building robust multigrid methods for convection-dominated problems is a well-known challenge; see the overview in [23]. To de?ne the necessary components of the multigrid method, let us consider two pairs of ?nite element spaces Vh × Qh and VH × QH , corresponding to ?ne and coarse grids, respectively. It was shown in [27, 10, 23] that a stable direct discretization on the coarse grid is the preferred choice in geometric multigrid methods for convection-di?usion problems. Additionally, in the case of elasticity-type problems (w = 0 in (2.1)), stability of the coarse grid problem means that velocity and auxiliary pressure approximations on the coarser grid remain compatible, so as to avoid the so-called locking phenomenon [1]. Thus, to build a matrix for a coarse grid problem in our multigrid method we consider the FE formulation ν(?uH , ?vH ) +
τ ∈TH

στ (w · ?uH , w · ?vH )τ + ((w · ?) uH , vH ) + (PH div uH , div vH ) = (f, vH ) ? vH ∈ VH ,

where the new set of parameters στ corresponds to the coarse grid and PH projects on QH . We further de?ne the problem-dependent “energy” norms on H1 , 0 v
2 h

= ν(?v, ?v) + (Ph div v, div v)

and

v

2 H

= ν(?v, ?v) + (PH div v, div v),

and set ah (u, v) = ν(?u, ?v) + (Ph div u, div v). Consider the standard prolongation operator p : VH → Vh induced by the embedding VH ? Vh . This is a common choice in a multigrid method for solving elliptic problems. However, it appears not

2104

MICHELE BENZI AND MAXIM A. OLSHANSKII

to work well for (5.1) in the case of small ν. The reason is that vH ∈ ker(PH div) does not yield in general vH ∈ ker(Ph div). Therefore, the standard prolongation is not uniformly stable in the energy norm; i.e., the constant C in the estimate (5.3) p vH
h

≤ C vH

H

? vH ∈ VH

could grow unboundedly for ν → 0. In [29] it is suggested that to enforce the stability of the prolongation one can make a correction of vH on a ?ne grid in a way involving solution of local problems. To present this idea in a general setting we consider the L2 -orthogonal decomposition (5.4) ? Qh = QH ⊕ Qh .

We prove the following lemma. ? Lemma 5.1. Assume the subspace Vh ? Vh is such that there is a constant c independent of h with (5.5) inf sup ? (?h , div vh ) q ≥c>0 ?? h qh v ?

qh ∈Qh vh ∈Vh ? ? ? ?

? ? and Vh ? ker(PH div). For some ?xed uH let uh be the solution of the problem (5.6) ? u ? ah (? h , vh ) = ah (uH , vh ) ? ? ? vh ∈ V h .

De?ne the prolongation p : VH → Vh by (5.7) ? p uH = uH ? uh .

Then p is stable; i.e., (5.3) holds with a constant C independent of h and ν. Proof. De?ne the following bilinear form on H1 × L2 : 0 0 B(u, p; v, q) := ν(u, v) + (p, div v) + (q, div u) ? (p, q). For given uH denote (5.8) pH = PH div uH .

? ? ? ? De?ne uh ∈ Vh and ph ∈ Qh by solving (5.9) ? ? ? ? ? ? B(? h , ph ; vh , qh ) = B(uH , pH ; vh , qh ) ? vh ∈ Vh , qh ∈ Qh . u ? ? ?

? This problem is well-posed thanks to (5.5). Taking vh = 0 in (5.9) we get (5.10) ? ? ? ? ? ? (div(uH ? uh ), qh ) = (pH ? ph , qh ) ? qh ∈ Qh .

? Also, since (5.8) implies (pH ?div uH , qH ) = 0 and the assumptions Vh ? ker(PH div) ? h ⊥ QH imply (?h , qH ) = (div uh , qH ) = 0 for any qH ∈ QH , we get ? p and Q (5.11) ? ? (div(uH ? uh ), qH ) = (pH ? ph , qH ) ? qH ∈ QH .

? Since Qh = QH ⊕ Qh , relations (5.10) and (5.11) yield (5.12) ? ? Ph div(uH ? uh ) = pH ? ph .

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2105

? Using (5.12) and (5.9) with qh = 0 we can see that uH and uh found from (5.9) satisfy ? ? ? ? ? ? ν(?? h , ?? h ) + (Ph div uh , div vh ) = ν(?uH , ?? h ) + (Ph div uH , div vh ) ? vh ∈ Vh . u v v This is the same as (5.6). Thanks to the stability assumption (5.5) the bilinear form B is stable with respect 1 to the norm (ν u 2 + ν ?1 p 2 )? 2 (this can be easily veri?ed by rescaling the problem to a Stokes problem with the penalty term ν(p, q)). Thus, for the solution of the problem (5.9) we obtain the estimate (5.13) ? (ν uh
2

+ ν ?1 ph 2 )? 2 ≤ C ?
1

? q vh ,?h ∈Vh ×Qh

sup

=C

? q vh ,?h ∈Vh ×Qh

sup

B(? h , ph ; vh , qh ) u ? ? ? 1 ? (ν vh 2 + ν ?1 qh 2 )? 2 ? ? ? B(uH , pH ; vh , qh ) ≤ 2C ν uH . 2 + ν ?1 q 2 )? 1 2 ? (ν vh ?h

Using the triangle inequality and (5.13) we get ? ν uH ? uh
2

+ pH ? ph ?

2

≤ C ((ν + ν 2 ) uH

2

+ pH

2

).

Thanks to (5.8) and (5.12) the last inequality is equivalent to ? ν uH ? uh
2

? + Ph div(uH ? uh )

2

≤ C (ν uH

2

+ PH div uH

2

),

which is (5.3) for the prolongation de?ned in (5.7). ? Note that for isoP2-P0 or P2-P0 ?nite elements, the spaces Qh and the appropriate ? h are easy to describe. Indeed, let TH be the coarse grid (pressure) triangulation. V Then ? Qh = qh ∈ Qh : ?
τ

qh dx = 0 ? τ ∈ TH ?

,

? ? Vh = {? h ∈ Vh : vh |?τ = 0 ? τ ∈ TH } . v The assumptions of Lemma 5.1 can be readily checked. For this choice of the sub? space Vh , solving (5.6) requires the solution of Np independent problems of small dimension (for regular subdivisions the dimension is 6 in two dimensions, and 9 in three dimensions), where Np is the number of elements in TH . On the other hand, ? for isoP2-P1 or P2-P1 ?nite elements we are unable to ?nd a simple choice for Vh leading to easily solved local problems (5.6). Therefore, in this case we compute the ? ? correction vh in (5.7) by solving locally (5.6) on each element τ ∈ TH , taking for Vh all vh ∈ Vh vanishing on ?τ . From experiments we saw that such a prolongation still improves the robustness of the multigrid method in the case of continuous pressure discretization compared to the canonical prolongation (? h = 0). Moreover, some furu ther improvement in this case was observed if the w-dependent terms were included in the de?nition of the bilinear form ah in (5.6). However, the observed improvement in this case was not as dramatic as for isoP2-P0 elements; see also Remark 5.2. Typically the restriction operator r : Vh → VH is taken to be the adjoint of the prolongation. This choice is also convenient for theoretical analysis [29]. However, an operator adjoint to (5.7) is not easily computable. Therefore, as a restriction we use the L2 -orthogonal projection from Vh to VH . It is well known that the L2 -orthogonal projection on VH is stable in the H 1 -norm (see, e.g., [4]) under some assumptions on TH . Although we did not ?nd in the literature any results which would allow us

2106

MICHELE BENZI AND MAXIM A. OLSHANSKII

 I 

I ? R



I 6 -

Fig. 5.1. Basis of weakly divergence-free functions for isoP2-P 0 and P 2-P 0. Pressures from Qh are constants on each triangle; dots denote velocity nodes.

to conclude the stability of the L2 -orthogonal projection with respect to the norms · h , · H if ν → 0, this choice works very well as the restriction operator in our multigrid method. Finally, we describe the smoother. As already mentioned, an e?ective smoother should e?ectively eliminate the oscillatory components in the kernel of Ph div. For piecewise constant pressure-based FE pairs one ?nds that uh ∈ ker(Ph div) if and only if for any τ ∈ Th the total ?ux for uh through ?τ equals zero. Thus, the kernel of Ph div can be described explicitly. For example, basis functions for ker(Ph div) in the case of P2-P0 ?nite elements are shown in Figure 5.1 [5, 29] (these are actually three of nine basis functions “captured” inside the oval on the right). Note that all basis functions have local support. Therefore an e?cient smoother can be built as a block relaxation procedure, where degrees of freedom (DOFs) supporting each basis function from ker(Ph div) are organized in one block. The blocks can be overlapping. This is described more formally below. We used a block Gauss–Seidel method, which can be written as a sequential subspace correction method (in the terminology of Xu [36]) based on the decomposition l Vh = i=0 Vi , Vi ? Vh , where the sum is not necessarily a direct sum. Let Ai be the sti?ness matrix for (5.2) on Vi . Then one relaxation iteration consists of the following steps: Let u0 = uold , then compute (5.14)
?1 ui+1 = ui ? pi A?1 ri ((A + B T Mp B) ui ? f ) i

for i = 0, . . . , l,

T and put unew = ul+1 . Restriction ri in (5.14) is the simple nodal one and pi = ri . For isoP2-P0 or P2-P0 ?nite elements the natural choice is to gather nodal DOFs for velocity inside ovals similar to the one in Figure 5.1. To be more precise, let si , i = 1, . . . , l, be the set of all interior vertices of Th . For each i de?ne τ i as a union of all triangles τ ∈ Th sharing si . Then

Vi := {vh ∈ Vh : vh = 0

in Ω \ τ i }.

We note that in the case of piecewise linear pressure FE pairs (e.g., isoP2-P1) it is harder to describe the kernel of Ph div. Thus, in this case we simply use the same choice of Vi to de?ne the relaxation step (5.14). As we will see, the resulting smoother remains quite e?ective. As usual for Gauss–Seidel type methods, the performance of the iterations (5.14) depends on the ordering of blocks of unknowns. Since the de?nition of Vi is solely based on the choice of a node si , the relaxation method (5.14) is uniquely de?ned by the ordering of nodes in Th . It is well known that smoothing iterations for convectiondominated problems are sensitive to the numbering of unknowns and that the fastest

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2107

convergence is achieved for the streamwise ordering of nodes; see, e.g., [3, 18, 23]. Finding appropriate orderings may be a di?cult task, especially for complex threedimensional (3D) ?ow ?elds w. However, our numerical experiments suggest that the multigrid method with the above-described block Gauss–Seidel smoothing applied to (5.1) is not very sensitive to the numbering of nodes. Because of this, in our implementation we used a lexicographical ordering of the blocks of unknowns. One smoothing iteration consists of two relaxation steps (5.14): one with left-right/top-bottom ordering, the other with top-bottom/left-right ordering. We conclude the section with some remarks about the implementation and performance of the multigrid scheme. Remark 5.2. While the multigrid method described above turns out to be very robust for (iso)P2-P0 ?nite elements, in the case of (iso)P2-P1 discretization some ν degradation of convergence is still observed as γ → 0. This results in a large number of outer iterations. One way to improve convergence is to increase the number of ν smoothing steps on each level as γ → 0. Another, more e?cient way (followed in the numerical experiments; see section 6) is to take somewhat smaller γ for ν 1, keeping the number of smoothing steps constant and low. Doing this, one ensures that the cost of one outer iteration does not increase. At the same time one sacri?ces slightly the quality of the Schur complement preconditioner in (3.7): Indeed, for γ = 0 the preconditioner from (3.7) is poor for small ν (see [8]). Nevertheless, the total number of outer iterations stays quite low for all reasonable values of ν (see section 6). Remark 5.3. In general we found the W-cycle to be more robust than the V-cycle. For the problems we tested, the V-cycle needs approximately twice or three times as many smoothing steps as the W-cycle to demonstrate similar convergence rates. Thus, the results in the next section are those obtained with the W-cycle implementation. However, for problems with locally re?ned grids using the V -cycle instead may be preferable in order to preserve the optimal complexity of the multigrid method. Remark 5.4. The dimension of the subspace Vi equals 14 in the 2D case (two velocity components for each node inside the oval in Figure 5.1). In the 3D case and for a regular subdivision we get dim Vi = 39. It is a common practice in using block smoothers for ?ow problems (such as Vanka smoothing [33]) that instead of applying A?1 in (5.14) one considers a simple, e.g., diagonal or triangular, approximation of i Ai [31, section 2.4.1]. This helps reduce signi?cantly the computational cost of one smoothing iteration. In this paper we solve subproblems with Ai exactly. A multigrid scheme based on inexact block inverses will be studied elsewhere. 6. Numerical experiments. In this section we present results for both isoP2P0 and isoP2-P1 elements. We focus on two commonly used model problems on the unit square. In the ?rst problem the wind function is constant and is given by (6.1) w= 1 0 .

The second problem is a classical recirculating ?ow problem. Here the wind is a “rotating vortex,” i.e., the vector ?eld (6.2) w= 4(2y ? 1)(1 ? x)x ?4(2x ? 1)(1 ? y)y .

In the experiments we use our preconditioner with the Krylov subspace method BiCGStab [32]. Inexact solves involving the (1,1) block consist of a single W(1,1)cycle with the multigrid method described in the previous section, that is, one W-cycle

2108

MICHELE BENZI AND MAXIM A. OLSHANSKII Table 6.1 Results for AL approach, isoP2-P 0 ?nite elements. Mesh size h 1 0.1 w is constant wind (6.1) 1/16 7 5 1/32 7 5 1/64 5 5 1/128 5 5 w is rotating vortex (6.2) 1/16 1/32 1/64 1/128 5 4 4 4 5 4 4 5 Viscosity ν 0.01 5 6 6 5 6 5 5 5 10?3 6 7 5 5 10 10 9 7 10?4 6 8 7 6 15 21 18 14

Number of preconditioned BiCGStab iterations ? (A?1 is one W(1,1)-cycle, γ = 1). γ

with one presmoothing and one postsmoothing step on each level. Each smoothing step consists of two iterations of block Gauss–Seidel relaxation with the alternating ordering described at the end of the previous section. The number of levels is such that the coarsest grid problem corresponds to h = 1 . The coarsest grid problem is 2 solved exactly. ? We set W = Mp in all experiments. Solves involving the (2,2) block are inexpen? sive. The action of S ?1 on a vector is computed via (3.7). For isoP2-P0 elements this amounts to a diagonal scaling operation. For the case of isoP2-P1 elements, the pressure mass matrix is not diagonal. Linear systems involving Mp are approximately solved by performing a ?xed number (?fteen) of stationary Richardson-type relaxations with parameter α = 1.75. Note that the condition number of Mp is bounded by a constant independent of h (and of ν). Since we use a ?xed number of stationary iterations for both the (1,1) and the (2,2) blocks, the overall preconditioner remains constant from one outer iteration to the next. Hence, we can use a standard Krylov subspace method such as BiCGStab for the outer iteration. In all runs we stop the outer iteration once the 2-norm of the initial residual has been reduced by at least six orders of magnitude, and we always use the zero vector as the initial guess. The results for isoP2-P0 elements are reported in Table 6.1. For both choices of the wind function, the convergence rate is h-independent. For the constant wind problem, the convergence rate is also independent of ν, even for values of ν as small as 10?4 . For the recirculating ?ow problem, the results are virtually the same except for an increase in the number of iterations for the smallest value of the viscosity, ν = 10?4 . Even in this extreme case, however, the convergence is quite fast. Results for the isoP2-P1 discretization are reported in Table 6.2. The results are still quite good. For the smallest values of the viscosity we found that γ must be reduced in order to retain fast convergence (see Remarks 4.3 and 5.2 for a discussion). The results are, in all cases, h-independent; some increase in the iteration count occurs for the smallest values of ν, but overall the convergence properties remain excellent. 6.1. Other methods. In this subsection we report the results obtained when solving the same model problems with the pressure convection-di?usion (PCD) preconditioner developed in [21] and with a full Vanka multigrid method as a preconditioner (see, e.g., [19, 33]). Both techniques are based on the original, nonaugmented

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM Table 6.2 Results for AL approach, isoP2-P 1 ?nite elements. Mesh size h 1 0.1 Viscosity ν 0.01 10?3 10?4

2109

Parameter γ 1 1 w is constant wind (6.1) 1/16 6 6 1/32 7 6 1/64 7 6 1/128 7 6 w is rotating vortex (6.2) 1/16 1/32 1/64 1/128 6 5 4 4 6 6 5 4 1 7 10 8 9 7 9 10 10 0.1 8 8 7 9 13 11 11 12 0.02 24 22 19 18 25 32 37 34

Number of preconditioned BiCGStab iterations ? (A?1 is one W(1,1)-cycle). γ

formulation (γ = 0). The PCG preconditioner is built upon separately solving velocity and pressure subproblems, while Vanka’s approach is a geometric multigrid method for a coupled problem (1.10); see the review paper [35] discussing this class of methods. The PCD method is widely regarded as one of the best available preconditioners in the class of block solvers for the steady Oseen problem, at least in the case of continuous pressure approximations. The preconditioner is of the form (3.1), where ? S, the approximate pressure Schur complement, is de?ned through the action of its inverse: (6.3)
?1 ? S ?1 = Mp Fp A?1 . p

In (6.3) Mp denotes the pressure mass matrix, Fp denotes a convection-di?usion operator acting on the pressure space, and Ap represents the pressure Laplacian; for the latter two operators, care should be taken to enforce the appropriate boundary conditions. For the given model problems we used Neumann conditions on the entire boundary. In the case of piecewise linear pressure approximations Fp was constructed following [10, page 348]. For piecewise constant pressure elements we set ?1 Ap := BMu B T (so-called mixed approximation of the pressure Poisson problem) and Fp := ru→p Ax rp→u , where Ax is the x-subblock of A, ru→p and rp→u are suitable mappings from Qh to Vh , and vice versa. Note that in two dimensions each application of the PCD preconditioner requires solving two convection-di?usion problems for the velocities and one Poisson-type problem for the pressure. These solves are usually performed inexactly in order to reduce costs. In our implementation we used exact solves for both the (1,1) block and for the inverses appearing in (6.3). As recommended in [10], we use the PCD preconditioner with BiCGStab(2) [30]. This Krylov subspace method requires four matrix-vector multiplies and four applications of the preconditioner per iteration. In order to facilitate the comparison with the results obtained with BiCGStab and our own preconditioner, we report in Table 6.3 the number of “BiCGStab-equivalent iterations” (i.e., the number of BiCGStab(2)

2110

MICHELE BENZI AND MAXIM A. OLSHANSKII Table 6.3 Results for PCD preconditioner, isoP2-P 0/isoP 2-P 1 ?nite elements. Mesh size h 1 0.1 w is constant wind (6.1) 1/16 6 / 12 8 / 16 1/32 6 / 10 10 / 16 1/64 6 / 10 8 / 14 1/128 6 / 10 8 / 12 w is rotating vortex (6.2) 1/16 1/32 1/64 1/128 6/8 6/8 4/6 4/6 10 / 12 10 / 12 8 / 12 8 / 10 Viscosity ν 0.01 12 / 24 14 / 24 16 / 24 16 / 26 30 / 40 30 / 40 26 / 40 22 / 44 10?3 30 / 34 24 / 28 22 / 32 24 / 36 > 400 / 188 > 400 / 378 > 400 / > 400 228 / > 400 10?4 100 / 80 86 / 92 64 / 66 64 / 58

Number of preconditioned BiCGStab-equivalent iterations (exact solves).

iterations multiplied by 2). We also experimented with BiCGStab, but the results were not better. It appears from these results that the PCD preconditioner is generally less robust and e?ective than the AL-based preconditioner, especially for small ν. We note, however, that PCD preconditioners can be built on a basis of various solvers for the (1,1) block and for the pressure Poisson problem; in particular, algebraic multigrid methods can be used, while the e?cient implementation of the AL-based preconditioner described in this paper requires a hierarchy of grids to be available, as usual for a geometric multigrid method. Also, it should be mentioned that the PCD preconditioner results in preconditioned operators for which some of the eigenvalues have an imaginary part growing like O(ν ?1 ); see [9]. It is well known that BiCG-type methods tend to perform poorly for such problems. Thus, it is possible that better results could be obtained for small ν if PCD preconditioning is used with GMRES or with BiCGStab( ) with larger values of . The Vanka preconditioner is built on a geometric multigrid method with block Gauss–Seidel smoothing for system (1.10). One relaxation step consists of l sequential updates: (6.4) ui+1 pi+1 = ui pi
T ?1 ? ri Ci ri A

ui pi

?

f 0

for i = 0, . . . , l.

Restriction ri is the simple nodal one. Entries of matrices Ci are those from A corresponding to unknowns in the ith block. For approximations based on discontinuous pressure, unknowns in one block consist of all pressure and velocity DOFs tied to one cell (triangle). For isoP2-P0 elements this results in a matrix of the form (6.5) Ci = Ai Bi
T Bi 0

,

where the block Ai is 12 × 12 and Bi is 1 × 12. Since for a regular subdivision in two dimensions the number of cells is approximately twice the number of vertices and the dimension of Ai in (5.14) is 14 × 14, we conclude that for isoP2-P0 elements one smoothing step in Vanka’s method is more computationally expensive than in the multigrid method from section 5.

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM Table 6.4 Vanka multigrid preconditioner, isoP2-P 0/isoP 2-P 1 ?nite elements. Mesh size h 1 0.1 w is constant wind (6.1) 1/16 4/6 4/6 1/32 4/6 4/6 1/64 4/6 4/6 1/128 4/5 4/5 w is rotating vortex (6.2) 1/16 1/32 1/64 1/128 4/4 4/4 4/4 4/4 5/5 4/5 4/4 4/4 Viscosity ν 0.01 3 / 10 3/8 4/8 4/6 5/7 6/6 5/5 4/5 10?3 4 / 18 4 / 22 3 / 21 3 / 14 8 / 19 10 / 20 8 / 18 9 / 13 10?4 5 / 20 4 / 33 5 / 56 4 / 64 13 / 30 17 / 38 40 / 59 > 400 / 67

2111

Number of preconditioned BiCGStab iterations.

For approximations based on the continuous pressure it is recommended in [19] that unknowns in a block should consist of all the pressure and velocity DOFs tied to one vertex. For isoP2-P1 elements this results in a matrix of the same structure as (6.5). However, now each matrix Ci has dimension 39 × 39 for a regular 2D triangulation (196 × 196 in three dimensions). This makes one relaxation step considerably more expensive than in (5.14). To facilitate comparison we include in a block only velocity DOFs in the corresponding vertex and next to it (i.e., exactly the same as in the multigrid method described in section 5). Taking into account one pressure unknown in the vertex, this results in a 15 × 15 matrix Ci . Thus the computational cost of one smoothing step is nearly the same as in the multigrid method used in the ALbased preconditioner. One problem with this reduction is that now the Ci matrix for isoP2-P1 is LBB-unstable and the relaxation fails to smooth the error in the pressure. As typical for stabilized methods, the solution is to replace the zero (2,2) “block” in (6.5) (its dimension is actually 1×1) by the matrix ?σ(x)I, where σ(x) is an elementdependent stabilization parameter. Numerical experiments reveal that multigrid convergence is very sensitive to the value of σ(x). An ad hoc choice which we found to be close to optimal for various values of h and ν is σ(x) = (ν/2 + 2 w(x) h2 )?1 , where w(x) is the Euclidean norm of vector w(x) in the corresponding vertex. Results for the Vanka multigrid preconditioner and both FE pairs are shown in Table 6.4. Similar to the multigrid used in the AL-based method, one smoothing iteration in Vanka consists of two relaxation steps with alternating ordering directions. From the results reported in the table, it appears that the Vanka-type multigrid preconditioner gives good results for a wide range of values of h and ν. For very small values of the viscosity, however, the performance of the method may in some cases deteriorate, and the preconditioned iteration may even fail to converge. 7. Conclusions. The paper presents a new preconditioning technique for ?nite element discretizations of the Oseen problem arising from Picard linearizations of the steady Navier–Stokes equations. The approach is based on an AL formulation of the resulting saddle point system. The preconditioner is block triangular and requires an e?cient approximate solver for the (1,1) block; the second block is essentially a scaled mass matrix and is easy to invert. The (1,1) block can be seen as the discretization of a nonconventional, nonsymmetric operator consisting of a linear elasticity operator plus a convection term. We have developed an e?ective inexact solver for this opera-

2112

MICHELE BENZI AND MAXIM A. OLSHANSKII

tor by combining a multigrid method of Sch¨berl (originally designed for symmetric o problems) with a suitable smoother. We have shown that assuming exact solves for the (1,1) block, the preconditioner clusters the eigenvalues in a small rectangular region of the complex plane around the point (1, 0), and that the sides of this rectangle are independent of both h and ν. Our numerical experiments on some 2D model problems indicate that in practice, a single W(1,1)-cycle is su?cient to achieve fast, h-independent convergence of the outer BiCGStab iteration. Furthermore, the rate of convergence appears to be essentially ν-independent, except for a slight degradation for extremely low values of ν. The preconditioner is especially e?ective for isoP2-P0 and P2-P0 elements, but excellent performance is also observed for discretizations using piecewise linear pressure approximations. In the latter case, a good choice of the parameter γ is important for ν 1. Finally, we mention that while in this paper we have focused on the steady 2D case, our approach can be extended to the solution of unsteady and 3D problems. Acknowledgments. We would like to thank Clark Dohrmann for his helpful comments on an earlier version of the paper. Thanks also to David Silvester, Andy Wathen, and an anonymous referee for their careful reading of the manuscript. This work was done while the second author was visiting the Department of Mathematics and Computer Science at Emory University during fall 2005. The hospitality and support of Emory University are gratefully acknowledged.
REFERENCES [1] I. Babuˇka and M. Suri, Locking e?ects in the ?nite element approximation of elasticity s problems, Numer. Math., 62 (1992), pp. 439–463. [2] M. Benzi, G. H. Golub, and J. Liesen, Numerical solution of saddle point problems, Acta Numer., 14 (2005), pp. 1–137. [3] J. Bey and G. Wittum, Downwind numbering: Robust multigrid for convection-di?usion problems, Appl. Numer. Math., 23 (1997), pp. 177–192. [4] J. H. Bramble, J. E. Pasciak, and O. Steinbach, On the stability of the L2 projection in H 1 (Ω), Math. Comp., 71 (2001), pp. 147–156. [5] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, New York, 1991. [6] A. N. Brooks and T. J. R. Hughes, Streamline upwind/Petrov–Galerkin formulation for convection dominated ?ows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg., 32 (1982), pp. 199–259. [7] C. R. Dohrmann and R. B. Lehoucq, A primal-based penalty preconditioner for elliptic saddle point systems, SIAM J. Numer. Anal., 44 (2006), pp. 270–282. [8] H. Elman and D. Silvester, Fast nonsymmetric iterations and preconditioning for the Navier–Stokes equations, SIAM J. Sci. Comput., 17 (1996), pp. 33–46. [9] H. C. Elman, D. J. Silvester, and A. J. Wathen, Performance and analysis of saddle point preconditioners for the discrete steady-state Navier–Stokes equations, Numer. Math., 90 (2002), pp. 665–688. [10] H. C. Elman, D. J. Silvester, and A. J. Wathen, Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics, Numer. Math. Sci. Comput., Oxford University Press, Oxford, UK, 2005. [11] B. Fischer, A. Ramage, D. J. Silvester, and A. J. Wathen, On parameter choice and iterative convergence for stabilised discretisations of advection-di?usion problems, Comput. Methods Appl. Mech. Engrg., 179 (1999), pp. 179–195. [12] R. Fletcher, Practical Methods of Optimization, 2nd ed., Wiley & Sons, Chichester, UK, 1987. [13] M. Fortin and R. Glowinski, Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, Stud. Math. Appl. 15, North–Holland, Amsterdam, New York, Oxford, 1983. [14] D. K. Gartling and C. R. Dohrmann, Quadratic ?nite elements and incompressible viscous ?ow, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 1692–1708.

AUGMENTED LAGRANGIAN APPROACH TO THE OSEEN PROBLEM

2113

[15] G. H. Golub and C. Greif, On solving block-structured inde?nite linear systems, SIAM J. Sci. Comput., 24 (2003), pp. 2076–2092. [16] A. Greenbaum, Iterative Methods for Solving Linear Systems, Frontiers Appl. Math. 17, SIAM, Philadelphia, 1997. ¨ [17] C. Greif and D. Schotzau, Preconditioners for the discretized time-harmonic Maxwell equations in mixed form, Numer. Linear Algebra Appl., to appear. [18] W. Hackbusch and T. Probst, Downwind Gauss–Seidel smoothing for convection dominated problems, Numer. Linear Algebra Appl., 4 (1997), pp. 85–102. [19] V. John, Higher order ?nite element methods and multigrid solvers in a benchmark problem for the 3D Navier–Stokes equations, Internat. J. Numer. Methods Fluids, 40 (2002), pp. 775–798. [20] V. John and G. Matthies, Higher order ?nite element discretizations in a benchmark problem for incompressible ?ows, Internat. J. Numer. Methods Fluids, 37 (2001), pp. 885–903. [21] D. Kay, D. Loghin, and A. Wathen, A preconditioner for the steady-state Navier–Stokes equations, SIAM J. Sci. Comput., 24 (2002), pp. 237–256. [22] G. Lube, Stabilized Galerkin ?nite element methods for convection dominated and incompressible ?ow problems, in Numerical Analysis and Mathematical Modelling, Banach Center Publ. 29, Polish Academy of Sciences, Warsaw, 1994, pp. 85–104. [23] M. A. Olshanskii and A. Reusken, Convergence analysis of a multigrid method for a convection-dominated model problem, SIAM J. Numer. Anal., 42 (2004), pp. 1261–1291. [24] M. A. Olshanskii and A. Reusken, Grad-div stabilization for Stokes equations, Math. Comp., 73 (2004), pp. 1699–1718. [25] M. A. Olshanskii, A low order Galerkin ?nite element method for the Navier–Stokes equations of steady incompressible ?ow: A stabilization issue and iterative methods, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 5515–5536. [26] A. Quarteroni and A. Valli, Numerical Approximation of Partial Di?erential Equations, Springer Ser. Comput. Math. 23, Springer, Berlin, 1994. [27] A. Ramage, A multigrid preconditioner for stabilised discretisations of advection-di?usion problems, J. Comput. Appl. Math., 110 (1999), pp. 187–203. [28] H. G. Roos, M. Stynes, and L. Tobishka, Numerical Methods for Singularly Perturbed Equations: Convection-Di?usion and Flow Problems, Springer, Berlin, 1996. ¨ [29] J. Schoberl, Multigrid methods for a parameter dependent problem in primal variables, Numer. Math., 84 (1999), pp. 97–119. [30] G. L. G. Sleijpen and D. R. Fokkema, BiCGStab( ) for linear equations involving unsymmetric matrices with complex spectrum, Electron. Trans. Numer. Anal., 1 (1993), pp. 11–32. [31] S. Turek, E?cient Solvers for Incompressible Flow Problems. An Algorithmic and Computational Approach, Lect. Notes Comput. Sci. Eng. 6, Springer, Berlin, Heidelberg, New York, 1999. [32] H. A. van der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 13 (1992), pp. 631–644. [33] S. P. Vanka, Block-implicit multigrid solution of Navier–Stokes equations in primitive variables, J. Comput. Phys., 65 (1986), pp. 138–158. [34] A. Wathen, Realistic eigenvalues bounds for the Galerkin mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449–457. [35] P. Wesseling and C. W. Oosterlee, Geometric multigrid with applications to computational ?uid dynamics, J. Comput. Appl. Math., 128 (2001) pp. 311–334. [36] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Rev., 34 (1992), pp. 581–613.


赞助商链接

更多相关文章:
更多相关标签:

All rights reserved Powered by 甜梦文库 9512.net

copyright ©right 2010-2021。
甜梦文库内容来自网络,如有侵犯请联系客服。zhit325@126.com|网站地图