# Boundary conditions in approximate commutator preconditioners for the Navier-Stokes equations.

1. Introduction. Consider the Navier-Stokes equations(1.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

on [OMEGA] [subset] [R.sup.d], d = 2 or 3. Here, u is the d-dimensional velocity field, p is the pressure, and v is the kinematic viscosity, which is inversely proportional to the Reynolds number. The value [lambda] = 0 corresponds to the steady-state problem and [lambda] = 1 to the case of unsteady flow. It is assumed that u satisfies suitable boundary conditions on [partial derivative][OMEGA], which is subdivided as follows:

[partial derivative][[OMEGA].sub.i] = {x [member of] [partial derivative][OMEGA] | u x n < 0}, the inflow boundary, [partial derivative][[OMEGA].sub.c] = {x [member of] [partial derivative][OMEGA] | u x n = 0}, the characteristic boundary, [partial derivative][[OMEGA].sub.o] = {x [member of] [partial derivative][OMEGA] | u * n > 0}, the outflow boundary.

Linearization and discretization of (1.1) by finite elements, finite differences or finite volumes leads to a sequence of linear systems of equations of the form

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

These systems, which are the focus of this paper, must be solved at each step of a nonlinear (Picard or Newton) iteration, or at each time step. Here, B and [B.sup.T] are matrices corresponding to discrete divergence and gradient operators, respectively and F operates on the discrete velocity space. In this paper, we focus only on div-stable discretizations where the corresponding (2, 2) block entry is identically zero.

In recent years, there has been considerable activity in the development of efficient iterative methods for the numerical solution of the stationary and fully-implicit versions of this problem. These are based on new preconditioning methods derived from the structure of the linearized discrete problem given in (1.2). An overview of the ideas under consideration can be found in the monograph of Elman, Silvester and Wathen [5]. A survey of solver algorithms and issues associated with saddle point systems appears in [1]. The key to attaining fast convergence lies with the effective approximation of the Schur complement operator

(1.3) S = [BF.sup.-1][B.sup.T],

which is obtained by algebraically eliminating the velocities from the system.

Two approaches of interest are the pressure convection-diffusion (PCD) preconditioner proposed by Kay, Loghin and Wathen [8] and Silvester et al. [12], and the least squares commutator (LSC) preconditioner developed by Elman et al. [3]. We will describe (variants of) these in Section 2. They are derived using a certain commutator associated with convection-diffusion and divergence operators, which we introduce here. Consider the linear convection-diffusion operator

(1.4) F = -v[[nabla].sup.2] + w x [nabla],

where the convection coefficient w is a vector field in [R.sup.d]. This operator is defined on the velocity space (we will not be precise about function spaces) and derived from linearization of the convection term in (1.1) via Picard iteration. We will also assume that there is an analogous operator [F.sup.(p)] defined on the pressure space. Let B denote the divergence operator on the velocity space. For two-dimensional problems,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now, define the commutator

(1.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where [B.sub.x] = [partial derivative]/[partial derivative]x and [B.sub.y] = [partial derivative]/[partial derivative]y. When [w.sub.1] and [w.sub.2] are constant and boundary effects are ignored, E = 0 in (1.5). This observation leads to the PCD and LSC preconditioners, as shown in Section 2.

This discussion does not take into account any effects that boundary conditions may have on the outcome. For example, the PCD preconditioner requires a discrete approximation [F.sub.p] to the operator [F.sup.(p)], and for this it is necessary that boundary conditions associated with [F.sup.(p)] be specified. In previous work, decisions about defining boundary conditions within the preconditioner have been made in an ad hoc manner. These decisions have been primarily guided by experimentation and they lack a solid basis. A poor choice of boundary conditions can critically affect the performance of the preconditioners.

In this study, we systematically study the effects of boundary conditions on commutators of the form (1.5), and we use these observations to define new versions of the PCD and LSC preconditioners. We show that for one-dimensional problems, it is possible to specify boundary conditions in such a way that the commutator is zero even at domain boundaries. In particular, a Robin boundary condition is used at inflow boundaries for the convection-diffusion operator [F.sup.(p)]. It is interesting to note that Robin conditions appear in several somewhat related preconditioning contexts. In [15, Section 11.5.1, p. 329] (see also the references therein), Robin conditions are discussed in terms of domain decomposition and how they enhance coercivity. In our context, the Robin condition leads directly to discrete (matrix) operators for which a corresponding discrete commutator is zero, and it results in a perfect PCD preconditioner in a one-dimensional setting. This basic idea is then generalized to higher dimensions by splitting the differential operators into components based on coordinate directions. This split leads to several equations associated with the commutator that can be analyzed in a fashion similar to the one-dimensional case. Although not all the commutator equations can be satisfied simultaneously, it is possible to satisfy an important subset of them. This again leads to an appropriate set of boundary conditions within the PCD method for coordinate-aligned domains. The resulting new PCD preconditioner, when combined with Krylov subspace solvers, displays significantly improved convergence properties for solving a problem with inflow and outflow boundary conditions. For enclosed flow problems (containing only characteristic boundaries), our analysis shows that the Neumann condition previously used [8] to define the convection-diffusion operator [F.sup.(p)] is, in fact, the right choice.

Our analysis also leads naturally to a modification of the LSC preconditioner that emphasizes important relationships for components of the commutator associated with boundaries.

The convergence of the resulting LSC-preconditioned GMRES iteration appears to be independent of the discretization mesh parameter. This contrasts with the (mild) mesh dependence seen for the standard LSC method and it also helps explain that this dependence is caused by boundary effects.

An outline of the paper is as follows. In Section 2, we briefly review the derivation of the PCD and LSC preconditioners. In Section 3, we present some preliminary results for a commutator associated with a one-dimensional model. In Section 4, we show how to split the differential commutator into components and examine the effects of directionality on properties of the components and of discrete versions of them. In Section 5, we analyze the effects of satisfying different relationships associated with the commutator. Specifically, we show that a component of the commutator corresponding to the direction orthogonal to a Dirichlet boundary plays a special role in the qualities of the commutator. New versions of the PCD and LSC method that take these considerations into account are derived in Section 6. In Section 7, we demonstrate the improved performance of the new preconditioners on two benchmark problems, the backward-facing step, whose boundary contains inflow, outflow, and characteristic components, and the driven cavity problem, which has only characteristic boundaries. Finally, in Section 8, we make some concluding remarks.

2. Review: preconditioners derived from approximate commutators. If finite element methods are used to discretize the component operators of E in (1.5), then a discrete version of the commutator takes the form

(2.1) E = ([Q.sup.-1.sub.p]B)([Q.sup.-1.sub.v]F) - ([Q.sup.-1.sub.p][F.sub.p])([Q.sup.-1.sub.p]B),

where [Q.sub.v] and [Q.sub.p] are the velocity mass matrix and pressure mass matrix, respectively. Assuming that this matrix version of the commutator is also small (i.e., E [approximately equal to] 0), a straightforward algebraic manipulation leads to the approximation

[(B[Q.sup.-1.sub.v][B.sup.T]).sup.-1] [F.sub.p][Q.sup.-1.sub.p]B[F.sup.-1][B.sup.T] [approximately equal to] I.

That is, the inverse of the Schur complement can be approximated by

(2.2) [(B[F.sup.-1][B.sup.T]).sup.-1] [approximately equal to] [A.sup.-1.sub.p][F.sub.p][Q.sup.-1.sub.p],

where

(2.3) [A.sub.p] = B[Q.sup.-1.sub.v][B.sup.T]

is a discrete Laplacian. If the commutator is small, we expect [A.sup.-1.sub.p][F.sub.p][Q.sup.-1.sub.p] to be an effective preconditioner for the Schur complement. In an implementation, a diagonal matrix (spectrally equivalent to [Q.sub.p]) can be used to approximate the action of [Q.sup.-1.sub.p], and similarly, [Q.sup.-1.sub.v] can be replaced by a spectrally equivalent approximation in (2.3) [17]. In the following, we will assume that [Q.sub.p] and [Q.sub.v] represent these diagonal approximate mass matrices. Moreover, one or two multigrid cycles applied to a Poisson equation can be used to approximate the action of [A.sup.-1.sub.p]. A complete specification of the PCD preconditioner (2.2) requires that the matrix [F.sub.p] be defined as though it comes from a differential operator (a convection-diffusion operator) having some boundary conditions associated with it; this is a primary focus of this paper

For the LSC method, [F.sub.p] is defined in an alternative fashion. The basic idea is to formally solve a least squares problem designed to make a discrete version of the commutator small. In particular, we solve a weighted least squares problem for the ith row of [F.sub.p][Q.sup.-1.sub.p] (equivalently, the ith column of [Q.sup.-1.sub.p][F.sup.T.sub.p]) via

(2.4) min [parallel] [[[F.sup.T][Q.sup.-1.sub.v][B.sup.T]].sub.:,i] - [B.sup.T][[[X.sup.T]].sub.:,i] [[parallel].sub.H],

where MATLAB-style notation is used to refer to matrix columns. Here, H is a positive-definite matrix defining an inner product and induced norm,

(2.5) [<q, q>.sub.H] = (H q,q), [[parallel]q[parallel].sub.H] = [<q, q>.sup.1/2.sub.H].

X is to be determined and produces an approximation to [F.sub.p][Q.sup.-1.sub.p]. Note that (2.4) is derived for the commutator with B; it is expressed in terms of BT only so that it looks like a least squares problem for a column vector.

The choice of the weighting matrix H and its relation to the commutatorwill be described in Section 6. Briefly here, the problem (2.4) is derived by multiplying E of (2.1) by [Q.sub.p], transposing, and then attempting to make the result small in the least squares sense with respect to the H-inner product. The normal equations associated with this problem are

B H [B.sup.T] [[[X.sup.T]].sub.:,i] = [[B H [F.sup.T] [Q.sup.-1.sub.v][B.sup.T]].sub.j].

This leads to the definition

(2.6) [F.sub.p] = (B[Q.sup.-1.sub.v] F H[B.sup.T])[(B H [B.sup.T]).sup.-1][Q.sub.p].

Substitution of this expression into (2.2) and using (2.3) then gives an approximation to the inverse of the Schur complement matrix,

(2.7) [(B[Q.sup.-1.sub.v][B.sup.T]).sup.-1](B[Q.sup.-1.sub.v]F H [B.sup.T])[(B H [B.sup.T]).sup.-1].

It is important to recognize that the matrix [F.sub.p] is never explicitly computed. Moreover, as above, (2.7) is used with [Q.sup.-1.sub.v] approximated by a diagonal matrix, and for practical computations the actions of [(B[Q.sup.-1.sub.v][B.sup.T]).sup.-1] and [(B H [B.sup.T]).sup.-1] are approximated by one or two multigrid cycles.

These (PCD and LSC) preconditioners differ in several ways from the "standard" ones described, for example, in [5]. Most notably, previously the commutator was defined using the gradient operator, as F[nabla] - [nabla][F.sup.(p)]. The difference from (1.5) appears innocuous since (for constant w) the components of the composite differential operators are the same. However, the commutator for the gradient leads to different requirements on boundary conditions than that for the divergence operator, and there are advantages to using the divergence operator. We will elaborate on this at the end of this section.

In addition to this essential difference, concerning boundary conditions, there are also a few structural differences between the new variants derived here and those previously developed:

* PCD preconditioning. Previously, manipulation of the commutator that uses the gradient operator led to the approximation

(2.8) [(B [F.sup.-1] [B.sup.T]).sup.-1] [approximately equal to] [Q.sup.-1.sub.p][F.sub.p][A.sup.-1.sub.p].

Thus the order in which the operators appear differs from (2.2).

* LSC preconditioning. The scaling matrix H has previously been taken to be the same as [Q.sup.-1.sub.v], a diagonal approximation to the inverse velocity matrix in [3] and it is taken as the reciprocal of the diagonal of F in [9] for highly variable viscosity Stokes problems. Here, we consider other diagonal forms of H in (2.7) to allow for different treatment of boundary effects. A derivation of (2.7) using the gradient would also change the order in which [Q.sub.v] and H appear.

* Both new methods use (2.3) explicitly for the discrete pressure Laplacian [A.sub.p], so that no decision is needed for boundary conditions for this operator.

These differences are negligible from the point of view of computational requirements.

3. One-dimensional analysis. Consider the one-dimensional operators

(3.1) F = [F.sub.p] = -v [d.sup.2]/d[x.sup.2] + w d/dx, B = d/dx

applied on an interval, say [OMEGA] = (0,1), where v and w are positive constants. We are interested here in the impact of boundary conditions on differential and discrete commutators, and we need not think of these operators as being associated with a problem such as (1.1) It is easy to see that the commutator BF - [F.sup.(p)] B is identically zero in [omega], with

(3.2) BF = [F.sup.(p)]B = -v [d.sup.3]/d[x.sup.3] + w [d.sup.2]/d[x.sup.2].

Let us assume that a Dirichlet value u(0) = 0 is specified at the inflow boundary x = 0 for functions u defined on [OMEGA], and that at the outflow boundary x = 1 we have a condition in which vu' = 0; the latter requirement is intended to be consistent with a standard outflow condition for (1.1); see (4.1). In addition, one boundary condition on B is required for the first-order equation Bu = 0 to have a unique solution. To be consistent with the left-to-right nature of the flow, we also take this to be a Dirichlet condition at the left boundary.

Now, given u, let v = Fu. For v to be appropriate as an argument of B on the left side of (3.2), we require v(0) = 0. Then

v = Fu = (-v d/dx + w) du/dx = (-v d/dx + w) p, p = du/dx.

That is, specifying the inflow value of v is the same as specifying the inflow value of(-v d/dx + w)p. But p is the argument of the operator [F.sup.(p)] in the second term of the commutator, [F.sup.(p)] d/dx. So, to make the commutator 0 at the inflow, we define [F.sub.p] applied to (any) p to satisfy the Robin boundary condition

(3.3) -vp' + wp = 0

at x = 0. This boundary condition more closely resembles a Dirichlet or a Neumann condition depending on the size of convection relative to diffusion. This is fundamentally different from what has been done previously [5], where either a Neumann or Dirichlet condition is imposed on a boundary depending on whether it corresponds to inflow, outflow, or no slip. The Robin condition is consistent with the analysis in [11] where it is shown that Dirichlet conditions on inflow are preferred as v approaches 0. Also, notice that with our choice of boundary condition for [F.sub.p], the third-order operator in (3.2) has two inflow boundary conditions specified at x = 0, u = 0 and -v [d.sup.2]u/d[x.sup.2] + w du/dx = 0.

A similar argument can be made at outflow, where F is defined with a Neumann condition and no additional boundary conditions are needed for B. Thus, the composite differential operator BFu is only required to satisfy a Neumann condition at outflow, i.e., [u.sub.x] = 0. For the composite differential operator [F.sup.(p)]Bu, we have p = Bu = [u.sub.x]. This implies that [F.sup.(p)]p must be equipped with a Dirichlet condition p = 0 at the outflow in order for the composite differential operator [F.sup.(p)]Bu to satisfy the same Neumann condition as BFu.

[FIGURE 3.1 OMITTED]

Now consider the discrete case, using the staggered mesh given in Figure 3. 1 corresponding to a one-dimensional version of a mesh that might be used for "marker-and-cell" (MAC) finite differences [7]. Given a subdivision of (0,1), by analogy with higher dimensions, take the discrete "velocities" (arguments of the matrix approximations F and B to F and B, respectively) to be defined on interval endpoints, and the discrete "pressures" (arguments of [F.sub.p] and [B.sup.T]) to be defined on interval centers. The discrete operators take the form of n x n matrices

(3.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

(3.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

We will explain in detail below how the first and last rows of F are determined. We show here how, given F, the requirement that the commutator BF - [F.sub.p]B be zero defines the first and last diagonal entries of [F.sub.p]. Suppose these entries are unknowns to be determined,

[[[F.sub.p]].sub.11] = [[xi].sub.1], [[[F.sub.p]].sub.nn] = [[xi].sub.n].

It is easy to show that requiring the first row of the products BF and [F.sub.p]B to be equal leads to the condition on the (1,1)-entries

[[BF].sub.11] = d(a + b) = [[xi].sub.1]d + bd = [[[F.sub.p]B].sub.11],

giving [[xi].sub.1] = a as in (3.5). The situation is slightly different for the right, outflow, boundary. Here, the last rows of BF and [F.sub.p]B are the same only if two conditions hold:

(3.6) (a + [[xi].sub.n])d = (a + b + c)d, [[xi].sub.n]d = (b + c)d.

Fortuitously, these conditions are compatible and give [[xi].sub.n] = b + c as in (3.5). Thus, it is possible to construct [F.sub.p] so that the discrete commutator BF - [F.sub.p]B is identically 0. It follows from this that

B[F.sup.-1] [B.sup.T] = [F.sup.-1.sub.p] B[B.sup.T],

i.e., the matrix [F.sub.p] can be used to produce a "perfect" preconditioner for B[F.sup.-1][B.sup.T].

The entries of F, B and the interior rows of [F.sub.p] are determined in a standard way. For concreteness, we will describe the case where centered differences are used for all convection terms. After scaling by [h.sup.2], this gives a = v + wh/2 and b = v - wh/2. The first row of F corresponds to an equation centered at [x.sub.1] with a Dirichlet value for [x.sub.0]. In the last row of F, the choice c = 2v comes from using a central discretization to the PDE and eliminating the ghost point via the outflow Neumann boundary condition. With d =1, B corresponds to a centered difference approximation, scaled by [h.sup.2]. The entries of [F.sub.p] centered at [x.sub.1/2] are determined using

(3.7) [[[F.sub.p]P].sub.1/2] = -(v + 2h/2)[p.sub.-1/2] + 2v[p.sub.1/2] - {y - wh/2)[p.sub.3/2].

Special treatment is needed for the value [p.sub.-1/2] at the "ghost point" [x.sub.-1/2] = -h/2. If we use the Robin condition (3.3), approximated at [x.sub.0] by

(3.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

then setting the difference operator on the right to zero leads to

(3.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Substitution into (3.7) gives

[[[F.sub.p]p].sub.1/2] = (v + wh/2)[p.sub.1/2] - (v - wh/2) [p.sub.3/2] = a[p.sub.1/2] - b[p.sub.3/2].

That is, the discrete Robin condition (3.3) is exactly what is needed to make the discrete commutator zero at the inflow boundary. It is also possible to interpret the discrete outflow condition as an approximation that incorporates a Dirichlet assumption. Specifically, the entries of [F.sub.p] centered at [x.sub.n-1/2] are

(3.10) [[[F.sub.p]p].sub.n-1/2] = -(v + wh/2) [p.sub.n-3/2] + (3v - wh/2) [p.sub.n-1/2].

This can be obtained by assuming a standard interior stencil and eliminating the ghost point using 1/2 {[p.sub.n-1/2] + [p.sub.n+1/2]) = 0 which is a discrete approximation to the Dirichlet condition p = 0. It should be noted that the MAC discretization is special in that the pressure grid contains no points on the boundary. Thus, each row of [F.sub.p] is an approximation of the differential operator (making use of boundary conditions for the first and last row). However, most discretizations have points on the boundary. Consider, for example, the discretization BF of the composite differential operator BF. The last row of F can be written as an approximation to F in (3.1) where a ghost point is removed using the outflow condition vu' = 0. Application of B then entails finite difference approximation to the derivative operator, applied along this row. Thus, BF would effectively be an approximation to (3.2) where a ghost point is removed using the outflow condition on F. It is not generally possible to exactly satisfy a discrete commuting relationship by taking a Dirichlet approximation for [F.sub.p] (e.g., [F.sub.p](k, n) = 0 for k [not equal to] n, [F.sub.p](n, n) = 1 where n is the dimension of [F.sub.p]). This Dirichlet condition implies that the last row of [F.sub.p]B is simply the last row of B which is an approximation to [u.sub.x]. Further, the scaling of a simple Dirichlet condition (e.g., taking instead [F.sub.p](n,n) = [alpha]) has an effect on how much the commuting relationship is violated, i.e., the size of [[parallel]BF - [F.sub.p]B[parallel].sub.2]. This follows from the fact that BF is independent of this scaling while [F.sub.p]B obviously depends on the scaling. This scaling issue does not arise in a standard discretization context so long as the right hand side is scaled appropriately. In our case, however, it is the scaling of [F.sub.p] which must be consistent with F. This will be discussed further in Section 7.

To summarize this discussion: the discrete commutator equation is solved exactly for a one-dimensional constant wind model problem with a MAC discretization by using a Robin condition at inflow and a Dirichlet condition at outflow. Although exact commuting is not always possible in other situations, the differential commuting relationships provide justification for using these conditions as a guide for higher dimensional settings.

Before moving on to higher dimensional problems, we briefly discuss the commutator with the gradient, which in the discrete one-dimensional setting has the form F [B.sup.T] - [B.sup.T] [F.sub.p]. With F as in (3.4), a derivation identical to that above produces a zero commutator at the left boundary with the choice [[[F.sub.p]].sub.11] = b. We can interpret this in terms of centered finite differences. Starting from (3.7), let the ghost point be defined using a discrete approximation to the boundary condition p'(0) = 0,

[p.sub.1/2] - [p.sub.-1/2]/h = 0 or [p.sub.-1/2] = [p.sub.1/2].

Substitution into (3.7) gives

[[[F.sub.p]p].sub.1/2] = {v - wh/2)[p.sub.1/2] - (v - wh/2) [p.sub.3/2] = b [p.sub.1/2] - b [p.sub.3/2].

Thus, the discrete commutator with the gradient is zero at the left boundary if [F.sub.p] is a discrete version of [F.sup.(p)] obtained using a Neumann condition at the inflow. In contrast to the case for the divergence, however, it is not possible to also make the discrete commutator zero at the outflow boundary. Moreover, with the Neumann inflow condition, the continuous operator [F.sup.(p)] becomes degenerate in the case of pure convection; that is, in the limit v [right arrow] 0, [F.sup.(p)] applied to any constant function is zero. No such difficulty occurs when the divergence is used for the commutator, where the Robin condition (3.3) induces a Dirichlet condition at the inflow.

4. Analysis for higher dimensions. Commuting ideas are now extended to problems in higher dimensions. We consider a two-dimensional rectangular domain [OMEGA] aligned with the coordinate axes (see Figure 4.1), where inflow and outflow conditions hold on the left and right vertical boundaries of [OMEGA], respectively. This corresponds to Dirichlet values for u on the left that satisfy [u.sub.1] > 0, and the solution satisfies

(4.1) v [partial derivative][u.sub.1]/[partial derivative]x - p = 0, [partial derivative][u.sub.2]/dx = 0,

on the right boundary with [u.sub.1] > 0. To simplify the discussion, we use periodic boundary conditions on the top and bottom.

[FIGURE 4.1 OMITTED]

Consider a splitting of the convection-diffusion operators F and [F.sup.(p)] into components associated with coordinate directions,

F = [F.sub.x] + [F.sub.y], [F.sup.(p)] = [F.sup.(p).sub.x] + [F.sup.(p).sub.y],

where

(4.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and [F.sup.(p)] is split analogously. We can use these to split the commutator (1.5) into four components,

(4.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where the first pair comes from a splitting of the first block of (1.5), and the second pair comes from the second block. In particular, E = [(1) + (2), (3) + (4)]. The expressions in (4.3) can be categorized by component orientations in relation to the vertical boundaries:

(1) orthogonal-orthogonal,

(2) orthogonal-tangential,

(3) tangential-orthogonal,

(4) tangential-tangential.

The first direction is associated with a component of B, and the second is associated with a component of F. For example, (2) contains the horizontal-oriented part of B ([partial derivative]/[partial derivative]x) which is orthogonal to the left and right boundaries, together with the vertical-oriented part of F, tangent to the left and right boundaries.

Our aim is to understand how the choice of boundary conditions for [F.sup.(p)] affects these commutators, and ultimately, to understand the impact this choice has on discrete versions of the commutators and the preconditioners. Indeed, although (4.2)-(4.3) correspond to differential operators, we are primarily interested in the discrete ones, and F, B and [F.sup.(p)] can be thought of as continuous approximations to their discrete analogues.

We begin with the discrete setting. The discrete analogue of (4.3) is

(4.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

For marker-and-cell finite differences, assume [OMEGA] is subdivided into an m x n grid; an example is shown in Figure 4.1 for m = 5 and n = 4. The component matrices are then structured as follows:

(4 5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and

(4.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [F.sup.(p).sub.x] are finite difference discretizations of [F.sub.x]. [F.sub.1], [F.sub.2] and [F.sub.3] are tridiagonal matrices corresponding to discretization along a single horizontal grid line. [F.sub.1] is identical to F of (3.4). [F.sub.2] is somewhat different due to the location with respect to the boundaries of the grid points labeled "[cross product]." [F.sub.3] near the boundaries is to be determined. [B.sub.1] also corresponds to discretization along a single horizontal grid line and is identical to B of (3.4). [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and [F.sup.(p).sub.y] are finite difference approximations to [F.sub.y]. They are identical to each other because of the periodic boundary conditions. All the approximations to [F.sub.x] and [F.sub.y] are block matrices where the block order is n and each of the individual blocks is of order m. It is easy to see that the entries of the commutators in (4.4) corresponding to interior points of [OMEGA] are zero when w is constant in (4.2).

Let us examine in detail what happens along boundaries. The discrete commutator (1) is the block-diagonal matrix

(4.7) diag([B.sub.1][F.sub.1] - [F.sub.3][B.sub.1], ..., [B.sub.1][F.sub.1] - [F.sub.3][B.sub.1]).

As [F.sub.1] and [B.sub.1] are identical to the matrices in the one-dimensional scenario described in Section 3, it follows that this expression is zero when [F.sub.3] at inflow includes a discrete version of the Robin condition

(4.8) -v[partial derivative]p/[partial derivative]x + [w.sub.1p] = 0.

Specifically, [[[F.sub.3]].sub.11] = a and [[[F.sub.3]].sub.12] = -b. Similarly, [F.sub.3] should include a Dirichlet condition at outflow. The complication not seen in the one-dimensional case is that [F.sub.3] also appears in the commutator (3). In particular, this discrete commutator is

(4.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

which is only zero if [F.sub.3] = [F.sub.2]. This implies that [F.sub.3] should include a Dirichlet condition at inflow and a Neumann condition at outflow, as these are the boundary conditions for [F.sub.x] (and so consequently for [F.sub.2]). For example, at inflow [[[F.sub.2]].sub.11] = 2a + b and [[[F.sub.2]].sub.12] = -b, and so taking [[[F.sub.3]].sub.11] = 2a + b and [[[F.sub.3]].sub.12] = -b makes (3) equal zero along the inflow boundary. (1) Thus, the conditions required to make (3) equal to zero are incompatible with the conditions required to make (1) equal zero. [F.sub.3] can be chosen to make either (1) or (3) in (4.4) equal to zero at the inflow boundary. However, they cannot both be zero simultaneously.

Now consider the commutators (2) and (4) of (4.4). Observe that for (2),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Thus, no special requirements are needed along the vertical boundaries in order for (2) and (4) to be zero; they are compatible with each other and they do not affect (1) and (3) as they do not depend on [F.sub.3]. This implies that it is possible for three of the four commutators in (4.4) to be zero simultaneously at each vertical boundary, but not all four. We note, however, that for centered finite differences, in the limit v [right arrow] 0, a = -b = [w.sub.1]h/2, so that the discrete commutator is identically zero at the inflow in the hyperbolic limit.

We summarize these observations with the assertion that the discrete operators (2), (3), and (4) of (4.4) are "self-commuting," in the sense that these commutators are zero when the discrete pressure convection-diffusion matrix [F.sup.(p)] is defined from the same boundary conditions used to specify the velocities. Furthermore, the property of self-commuting depends only on the special matrix structures and not on values in the particular stencils, as [F.sub.2], d, q, and l are arbitrary. In fact, self-commuting does not even depend on the specific boundary conditions per se; instead it is a consequence of the fact that (2), (3), and (4) contain at most one discrete difference operator associated with a direction orthogonal to the vertical boundaries. In particular, (4) comes from only tangential differencing, which, in light of the periodic boundary conditions at top and bottom, are represented by circulant matrices. It is well known that circulant matrices commute, which makes (4) equal to zero. Commutators (2) and (3) come from an orthogonal/tangential pair. The tangential differencing leads to circulant matrices with scaled identity submatrices (2) whereas the orthogonal differencing leads to a block diagonal matrix with the same submatrix in each block diagonal entry. Here, self-commuting relies on the fact that scalar multiplication with a matrix commutes. To summarize, the order in which the individual operators appear in (2)-(4) is not important in the discrete setting. If one views the continuous commutators (4.3) as approximations to the discrete commutators, one can loosely make an argument that the order of operators is not important in a continuous setting either (though the continuous commutators may not be identically zero). Further, (1) (corresponding to the one-dimensional problem) is unique in that it has an orthogonal-orthogonal characterization and it is the only discrete commutator that does not involve self-commuting.

This discussion has focused on boundaries assumed to be of either inflow or outflow type. We conclude this section with an observation for the case where the velocities satisfy a characteristic Dirichlet boundary condition. This holds along top and bottom boundaries in a simple channel flow application which is similar to the problem under discussion here, where w = [(1,0).sup.T], and it also applies to the driven cavity problem. In this case, it is condition (4) that is of (orthogonal,orthogonal) type with respect to both the top and bottom boundaries, and a zero commutator is obtained using the Robin boundary condition analogous to (4.8) to define [F.sup.(p).sub.y] along these boundaries. For [w.sub.2] = 0, this Robin condition is

(4.10) -v[partial derivative]p/[partial derivative] + [w.sub.2]p = -v[partial derivative]p/[partial derivative]y = 0,

at the bottom, a pure Neumann condition. (The top is the same, with the opposite sign for [partial derivative]p/[partial derivative]y.) This is what has been used previously for characteristic boundaries [5, 8], and this observation provides justification for this choice. As above, this would not be compatible with the Dirichlet condition needed to make the commutator (2) equal to zero. Our computational experience, which we report in Section 7, is that the Robin condition is to be preferred.

To summarize, three out of four commutator equations can be satisfied along the boundaries of [OMEGA], but there is a conflict between the commutators of "orthogonal-orthogonal" and "tangential-orthogonal" types. When Dirichlet conditions are specified for the Navier-Stokes equations, such as at an inflow boundary, Robin boundary conditions for [F.sup.(p)] enables the first of these commutator types to be zero. We will explore the effect of this choice in the following sections.

We comment on the analogous question for three-dimensional problems, where there are nine commutator equations due to three-way splittings of both B and F corresponding to the three coordinate directions. Six commutator equations (similar to (2) and (4) in 2D) are satisfied, as they do not contain differentiation orthogonal to the inflow boundary within the split F operator. Two of the remaining three conditions can be satisfied using a Neumann condition or, alternatively, one of the remaining three can be satisfied with a Robin condition. Further study would be required to fully assess this choice; our intuition is that the Robin condition remains the most important.

5. Perturbation analysis. As not all commutators in (4.4) can be zero simultaneously, we now seek to understand the impact of nonzero commutators, using a combination of analysis and empirical results for the MAC discretization. Once again, we consider a rectangular domain where the velocities satisfy periodic boundary conditions on the top and bottom, a Dirichlet inflow condition on the left boundary, and an outflow condition (4.1) on the right boundary. The PCD preconditioner (2.2) approximates the inverse of the Schur complement S = B[F.sup.-1][B.sup.T]. For MAC discretization on uniform grids, the mass matrices [Q.sub.v] and [Q.sub.p] are both diagonal of the form [h.sup.2]I, and they cancel each other in (2.2). Hence, the inverse of the preconditioning operator is

[M.sup.-1] = [(B[B.sup.T]).sup.-1][F.sub.p] = ([B.sub.x][B.sup.T.sub.x] + [B.sub.y][B.sup.T.sub.y])([F.sup.(p).sub.x] + [F.sup.(p).sub.y]).

The convergence of a preconditioned iterative method is (for "right-oriented" preconditioning) governed by properties of [SM.sup.-1]. We explore the variant obtained from a similarity transformation,

(5.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

are the commutator errors The similarity transformation will allow us to decouple the analysis of different components of the boundary

Further simplification of (5.1) yields

(5.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

If [E.sub.xx] = [E.sub.yy] = [E.sub.xy] = [E.sub.yx] = 0, then Y = J and M is an ideal preconditioner. We will explore the size of the perturbations from the identity as a measure of the effectiveness of M as an approximation to S, i.e., as a preconditioning operator.

As shown in Section 4, we can choose [F.sup.(p).sub.y] so that [E.sub.xy] and [E.sub.yy] are zero near the inflow boundary. Let [F.sup.(p, [perpendicular to]).sub.x] be the version of [F.sup.(p).sub.x] determined from Robin boundary conditions that makes [E.sub.xx] (i.e., (1)) equal zero along the inflow, let [F.sup.(p,|).sub.x] similarly be the operator obtained from Dirichlet inflow conditions, for which [E.sub.yx] ((3)) is zero, let [delta][F.sup.(p).sub.x] = [F.sup.(p, [perpendicular to]).sub.x] - [F.sup.(p,|).sub.x]. Let [Y.sup.[perpendicular to]] be the version of Y obtained when [F.sup.(p).sub.x] = [F.sup.(p,[perpendicular to]).sub.x] and let [Y.sup.|] be the version obtained when [F.sup.(p).sub.x] = [F.sup.(p,|).sub.x]. It follows that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Notice that only the factor [delta][F.sup.(p).sub.x] is tied directly to the size of errors in the commutators. Furthermore, each row of the perturbations is associated with a commutator error at one pressure grid point. That is, the commutator error at a particular pressure point only affects the row in the perturbed matrix associated with that point. This means that we can explore the effect of each boundary in isolation. For example, errors in the commutator associated with points adjacent to the inflow (outflow) boundary do not have any effect on entries in rows of the perturbation associated with points adjacent to outflow (inflow) boundaries.

To examine these perturbations, we consider a 40 x 40 MAC mesh with [w.sub.1] = 1, [w.sub.2] = 0, v = .25, and centered finite differences for both the convection and diffusion terms. The computed condition numbers of the preconditioned Schur complement using the orthogonal commutator (from [F.sup.(p,[perpendicular to]).sub.x]) and the tangential commutator (from [F.sup.(p,|).sub.x]) are 45.2 and 4404.1 respectively. Thus, the system is much better conditioned using the orthogonal commutator. This correlates well with the size of the inflow perturbations associated with using each of the two perturbation formulas:

(5.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

[FIGURE 5.1 OMITTED]

Specifically, Figure 5.1 illustrates a single row of [P.sup.[perpendicular to]] and [P.sup.|] associated with the 21st vertical pressure point adjacent to the inflow boundary. In this and several other figures below, this matrix row is displayed as a grid function on the underlying 40 x 40 grid. (3) It is obvious that the perturbation associated with the orthogonal commutator is small and localized. In addition, the Euclidean norm of the row vector from [P.sup.[perpendicular to]] depicted in the left side of Figure 5.1 is approximately 1/2. In contrast, the perturbations associated with the tangential commutator, shown on the right side of the figure, are much larger and do not decay quickly to zero. The Euclidean norm of the row vector from [P.sup.|] is 7.34.

Comparison of the two perturbation formulas reveals that the only differences are the appearance of either [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the appearance of either [B.sub.x] or [B.sub.y], and the sign of the perturbation term. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are almost identical as they correspond to the same convection-diffusion operator and both [u.sub.1] and [u.sub.2] satisfy Neumann conditions on outflow and Dirichlet conditions on all other boundaries. (4) Thus, the difference in magnitude of the two perturbations [P.sup.[perpendicular to]] and [P.sup.|] must be due to the presence of one or the other of [B.sub.x] or [B.sub.y].

We can get a clear understanding of [P.sup.[perpendicular to]]. Using the fact that [B.sub.y] and [B.sup.T.sub.y] commute with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and B[B.sup.T] due to the periodic boundary conditions on the top and bottom boundaries, we can rewrite this perturbation as

(5.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Each inflow row of [delta][F.sup.(p).sub.x] is nonzero in precisely one entry, corresponding to a grid point next to the inflow boundary; this difference comes from the difference in boundary conditions. Let [r.sup.T] denote one such row, so that the corresponding row of [P.sup.[perpendicular to]] is given by (the transpose of)

(5.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where r can now be interpreted as a discrete point source. Figure 5.2 plots the intermediate quantity s [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] as a grid function, where r comes from the 21st vertical point adjacent to the inflow boundary. It can be seen that this function has a large variation in the x-direction, and variations in the y-direction that are small everywhere, largest near the inflow boundary, and becoming smoother as one moves away from the inflow boundary. Consequently, application of the vertical discrete Laplacian [B.sub.y][B.sup.T.sub.y] to s (see (5.5)) largely eliminates the variation in the x-direction, and the resulting row of [P.sup.[perpendicular to]] is small. Our speculation is that the shape of the function depicted in Figure 5.2 is tied directly to the fact that an inflow row of [delta][F.sup.(p).sub.x] is a point source (centered adjacent to the inflow boundary) and that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] enforces a Neumann condition on the left boundary and a Dirichlet condition on the right boundary.

[FIGURE 5.2 OMITTED]

An expression analogous to (5.4) for [P.sup.|] is

(5.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Because the horizontal variation of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is large, application of the horizontal discrete Laplacian [B.sub.x][B.sup.T.sub.x] is not likely to produce a small result. It should be noted however that the "commuting trick" used to obtain (5.4) cannot be used in this case, so that [P.sup.|] is not in fact equal to the expression of (5.6), and this is merely a heuristic observation.

The conclusion reached from this discussion is that Robin boundary conditions for [F.sup.(p)] at the inflow make the preconditioned operator more resemble the identity than Dirichlet conditions, suggesting that Robin conditions are to preferred.

One can also compare perturbations for outflow by choosing [F.sup.(p).sub.x] to make either (1) (Dirichlet conditions) or (3) (Neumann) zero at the outflow boundary. Figure 5.3 illustrates perturbations for a single outflow row. There is no clear sense as to which perturbation is superior to the other. The Euclidean norm of the orthogonal perturbation is a little larger (2.16 vs. 1.05 for the tangential perturbation), the signs of the perturbations are different, and [P.sup.[perpendicular to]] extends somewhat further into the domain away from the outflow boundary than [P.sup.|]. This suggests that there might be a slight advantage for Neumann conditions at the outflow boundary, although the differences here are less conclusive.

[FIGURE 5.3 OMITTED]

In the figures discussed above, v was fixed at 1/4. Figure 5.4 shows [P.sup.[perpendicular to]] and -[P.sup.|] along the single vertical lines next to either the inflow or outflow boundaries in the two-dimensional grid, for different values of v. Specifically, rows adjacent to the inflow (or outflow) boundaries are considered, and to compress the presentation, only the values along the leftmost (for inflow) or rightmost (for outflow) grid line are plotted. One can see that as v decreases, the difference between the orthogonal and tangential inflow perturbations also decreases. This is most likely due to the fact that the Robin condition (from the orthogonal commutator) approaches a Dirichlet boundary condition (which is obtained with the tangential commutator) as v is decreased.

In summary, we conclude that at inflow boundaries, the orthogonal commutator (coming from Robin conditions) gives rise to smaller perturbations than the tangential commutator. The difference between these two commutators becomes less significant as the Reynolds number increases (i.e., as v decreases). The perturbations at outflow boundaries are roughly of equal size. Numerical experiments that augment these observations are presented in Section 7.

6. New variants of the PCD and LSC preconditioners. We use the results of Sections 4 and 5 to develop new variants of the pressure convection-diffusion and least squares commutator preconditioners. Both the PCD and LSC preconditioners are defined by the preconditioning matrix

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

When S is the true Schur complement, the preconditioned GMRES method converges in two iterations [10]. For the PCD method we approximate [S.sup.-1] via (2.2)-(2.3), where [F.sub.p] is a discrete convection-diffusion operator on the pressure space, together with the specification of boundary conditions for [F.sub.p]. The results above identify new criteria for these conditions. If the Navier-Stokes equations (1.1) are posed with Dirichlet conditions of either inflow or characteristic type, then we will define [F.sub.p] with Robin conditions

(6.2) -v[partial derivative]p/[partial derivative]n + (w x n)p,

as in (4.8) and (4.10). This reflects a decision to favor the commutator of (orthogonal-orthogonal) type, (1) in (4.4), as suggested by the results of Section 4. It constitutes a new strategy at inflows, where previously a Dirichlet condition was used for [F.sub.p], and which can now be seen to favor the commutator of tangential-orthogonal type, (3) in (4.4)). We will give additional evidence of the superiority of the new approach in Section 7. For characteristic boundaries, the new approach (Neumann conditions for [F.sub.p]) is the same as what was done previously. We also note that although the discrete operators discussed above were obtained for MAC finite differences, in the following we will use the new guidelines in a finite element setting. Details for specifying [F.sub.p] in this setting are given in [5, Ch. 8].

[FIGURE 5.4 OMITTED]

The new variant of the LSC preconditioner is defined by (6.1) where [S.sup.-1] is given by (2.7) and the scaling matrix H is chosen in a such a way that the discrete commutator (1) in (4.4) is given appropriate emphasis. To motivate this strategy, we begin with a version of the commutator of (2.1), scaled by [Q.sub.p]:

(6.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where the two terms on the right here could be further split into sums of the form (1)+(2) and (3)+(4). This scaling is used to make the derivation similar to what was done for finite differences in the previous two sections. We seek an approximation X to [F.sub.p][Q.sup.-1.sub.p]. For each row i of X, we will attempt to minimize the difference

[[B[Q.sup.-1.sub.v]F].sub.i,:] - [X.sub.i,:]B

in a least squares sense using a weighted norm that forces (1) to be more heavily emphasized than (3) in (6.3).

This will be achieved using a weighting matrix. We have

(6.4) [parallel][([B[Q.sup.-1.sub.v]F].sub.i,:] - [X.sub.i,:]B) [H.sup.1/2][parallel] = [[parallel][B.sup.T][X.sub.:,i] - F[Q.sup.-1.sub.v][B.sup.T.sub.i,:][parallel].sub.H],

where the norm in the expression on the left is the vector Euclidean norm, and H is a symmetric positive-definite matrix that induces an H-norm as in (2.5). We choose H = [W.sup.1/2][Q.sup.-1.sub.v][W.sup.1/2] where, as above, [Q.sup.-1.sub.v] is a diagonal approximation to the velocity mass matrix and W is a diagonal weighting matrix. Note that that the jth column of the commutator on the left of (6.4) is scaled by the jth diagonal entry of W. When the domain boundaries are aligned with the coordinate axis, we take [W.sub.jj] = 1 for all j except those corresponding to velocity components that are near the boundary (i.e., defined on a node contained in an element adjacent to [partial derivative][OMEGA]) and oriented in a direction tangent to the boundary. For those indices, we take [W.sub.jj] = [epsilon], a small parameter. Specifically,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The effect of this is that for each row i corresponding to a pressure unknown on a boundary, the weight is small in all columns j corresponding to the velocity component tangent to that boundary, and the associated component (3) is deemphasized in the least squares problem (6.4). (5)

This strategy is intended to mimic the treatment of boundary conditions used for the PCD preconditioner. We also point out its impact in the Stokes limit. For Stokes problems, a good approximation to the inverse of the Schur complement is [Q.sup.-1.sub.p] which gives rise to iterative solvers with convergence rates independent of discretization mesh size [16, 13]. In light of (2.2), this suggests that [F.sub.p] should be the same as [A.sub.p] (including boundary conditions) in the Stokes case. The boundary conditions of [A.sub.p] of (2.3) are completely determined by those associated with B. It is easy to show that Dirichlet conditions in the original system give rise to Neumann conditions for [A.sub.p]. Emphasis on (3) makes LSC's version in [F.sub.p] of (2.6) more like a Dirichlet condition (which obviously does not match the Neumann condition for [A.sub.p]). Emphasis on (1) in the new LSC preconditioner has the effect of forcing [F.sub.p] to more closely correspond to an operator defined with Neumann conditions (which properly matches the Neumann condition for [A.sub.p]).

Finally, we note that [A.sub.p] specified by (2.3) is non-degenerate only when the underlying discretization of (1.1) is div-stable [5, 6], so that B is of full rank (or rank-deficient by one in the case of enclosed flows). Techniques for defining [A.sub.p] and other operators arising in the LSC preconditioning for finite element discretizations that require stabilization are discussed in [2]. We expect these ideas to carry over to the techniques discussed in this study, although we do not consider this issue here.

7. Numerical results. We now show the results of numerical experiments with the new variants of the PCD and LSC preconditioners, addressing the following issues:

* For the PCD preconditioner, different choices of boundary conditions in the definition of [F.sup.(p)] enable different components of the commutators (4.3) and (for MAC discretization) (4.4) to be zero. We compare the different variants of the new PCD preconditioner defined using various choices of boundary conditions at inflow and outflow boundaries.

* We compare the new versions of both the PCD and LSC preconditioners with the original versions of them discussed, for example, in [5].

* We compare the new versions of the PCD and LSC preconditioners.

[FIGURE 7.1 OMITTED]

Two benchmark problems were used, a backward facing step, which is an example of an inflow/outflow problem, and a regularized driven cavity, which contains an enclosed flow with only characteristic boundaries. Figure 7.1 shows examples of the streamlines for the two problems. For the step, the inflow is at x = -1 and the outflow is at x = 5 for Reynolds numbers 10 and 100, and at x = 10 and x = 20 for Reynolds numbers 200 and 400, respectively. The examples were generated using the IFISS software package [14] and discretized with with a div-stable [Q.sub.2]-[Q.sub.1] finite element discretization (biquadratic velocities, bilinear pressures) on a uniform mesh of velocity element width 2/[2.sup.l-1] (giving velocity nodal mesh width 2/[2.sup.l]). The new algorithms were also tested with the same benchmark problems and a [Q.sub.2]-[P.sub.-1] element (which contains discontinuous pressures); the results were qualitatively the same as those shown below. Additional details concerning these problems are given in [5, Ch. 8].

The nonlinear discrete systems were solved using either a Picard or Newton iteration, which was stopped when the relative accuracy in the nonlinear residual satisfied a tolerance of [10.sup.-5]. Each step of the nonlinear iteration requires a linear system solve, where the coefficient matrix is a discrete Oseen operator for Picard iteration or a Jacobian matrix for Newton iteration. The results reported below show the performance of preconditioned GMRES for solving the last system arising during the course of the nonlinear iteration.

Consider the performance of the new PCD preconditioner. At inflow and outflow boundaries, the matrix [F.sub.p] used by the PCD preconditioner can be defined from four possible combinations of boundary conditions, corresponding to either Robin or Dirichlet conditions at inflow boundaries and Neumann or Dirichlet conditions at outflow boundaries. Figure 7.2 shows the performance of these four variants for four systems arising from the backward facing step. Neumann conditions are used along the top and bottom characteristic boundaries.

For comparison, the figure also contains results for the "original PCD" method of (2.8), which uses Dirichlet conditions for [F.sub.p] and [A.sub.p] at the inflow and Neumann conditions on all other boundaries [5, p. 348]. In these tests, the discretization parameter was l = 6, giving nodal velocity grids of size 64 x 96 for Re [less than or equal to] 100 (top of figure), 64 x 192 (Re = 200, bottom left) and 64 x 384 (Re = 400, bottom right). It can be seen that a choice of [F.sub.p] with Robin conditions at the inflow boundary is more effective than when Dirichlet conditions are used. This means that we are favoring commutator (1) at the inflow. The situation at the outflow is less clear, and for small Re, a Dirichlet condition at the outflow (which favors commutator (3)) appears to offer some advantage; we will return to this point in a moment. The difference between the two choices becomes negligible as Re increases. This can be attributed to the fact that as v [right arrow] 0, the Robin condition (6.2) is close to a Dirichlet condition. For large Re, the new variants are more efficient than the original version of the PCD preconditioner, in large part also because of a significantly shorter transient period of slow convergence.

[FIGURE 7.2 OMITTED]

Returning to the question of Dirichlet conditions for [F.sub.p] at the outflow, we note that implementation of Dirichlet conditions in a boundary value problem entails adjusting the coefficient matrix so that the known boundary values are obtained. A typical strategy is to force appropriate rows of the coefficient matrix to contain only diagonal entries, with values equal to 1. Of course, the value 1 is arbitrary and other values can be chosen for solving discrete PDEs, as long as the right-hand side is also properly adjusted. Here, we are interested in commuting and there is no right-hand side, so the choice for this value is less clear. We have found that with a strategy of this type for specifying [F.sub.p], the performance of the PCD preconditioner is sensitive to the scale of these diagonal entries. The results of Figure 7.2 come from taking these entries to be the average of the diagonal values of [F.sub.p] from all rows not corresponding to inflow and outflow boundaries. Intuitively, choosing scalings that are roughly comparable with diagonal entries of F (or of [F.sub.p]) is consistent with trying to make discrete commutators small. Figure 7.3 shows what happens when other scalings are used at the outflow, for one example from Figure 7.2 (discrete Jacobian, Re = 100). It is evident that performance is sensitive to this choice, and with other scalings, it is not better than when a Neumann condition are used at the outflow. Since the latter strategy does not entail a parameter, we use that in subsequent tests. Thus, for problems with inflow and outflow boundaries, we define [F.sub.p] with Robin conditions at the inflow (favoring commutator (1) of (orthogonal,orthogonal) type), and Neumann conditions at the outflow (favoring (3) of (tangential,orthogonal) type). These choices are consistent with the conclusions about perturbations reached in Section 5.

[FIGURE 7.3 OMITTED]

We have seen that with the new ordering of operators in the preconditioner (contrast (2.2) and (2.8)), it is important to also use the new Robin boundary conditions for [F.sub.p] at inflow boundaries. (Dirichlet conditions are the "old" choice.) One could ask the opposite question, how the new boundary conditions would work with the original ordering. To explore this, we show in Table 7.1 the iteration counts for three versions of each of the preconditioners, for the backward facing step and Re = 200; the stopping criterion was [[parallel][r.sub.k][parallel].sub.2]/[[parallel]b[parallel].sub.2] < [10.sup.-6] where [r.sub.k] and b are the residual and right-hand side vectors, respectively. The results for the new ordering of operators and the new treatment of boundaries are in the right column for each preconditioner. The middle column (for each) uses the old ordering with the new boundary conditions. For LSC, this means the inverse of the preconditioner is

[(B H [B.sup.T]).sup.-1](B H F [Q.sup.-1.sub.v] [B.sup.T])[(B[Q.sup.-1.sub.v][B.sup.T]).sup.-1],

which is derived using weighted least squares and commuting with the gradient operator (compare with (2.7). The results here are mixed: for the PCD preconditioner with the old ordering, the (new) Robin boundary condition is less effective than the (old) Dirichlet condition, whereas the opposite conclusion holds for the LSC preconditioner. However, the best combination is the new methodology. We note, however, that for the driven cavity problem, which has only characteristic boundaries, the old and new boundary conditions are the same, and we have also found that performance is not sensitive to operator ordering.

In general we have found the Jacobian systems arising from a Newton iteration to be slightly more costly than the Oseen systems, and this trend holds here as well. For the remainder of this section, we restrict our attention to the Oseen problem.

[FIGURE 7.4 OMITTED]

Figure 7.4 shows the performance of the PCD preconditioner for fixed values of the Reynolds number (denoted Re) and various grid refinements. The two graphs each show the results for both the new and original versions of the preconditioner. The graph on the left is for the driven cavity problem, where the original PCD preconditioner is known to exhibit mesh independent behavior [5]. Here, the eight curves (four for the new preconditioner and four for the original) are largely indistinguishable. The graph on the right is for the step problem. Here, the convergence behavior of GMRES with the new method is essentially independent of the the mesh, whereas the initial transient exhibited by the original method increases as the mesh is refined.

This issue was previously not well understood. The results shown in Figures 7.2 and 7.4 clearly demonstrate that boundary conditions are responsible for the difference between the original and new methods. For the cavity problem, with characteristic boundaries, Neumann boundary conditions are the right choice for defining [F.sub.p]. For the inflow/outflow (step) problem, with the wrong choice of boundary conditions for [F.sub.p], GMRES exhibits a long period of slow convergence in its initial steps. The new method produces mesh-independent convergence for both types of problems; in contrast, previously, only the asymptotic convergence rate was known to be independent of the mesh [4].

[FIGURE 7.5 OMITTED]

Figure 7.5 shows analogous performance results for the LSC preconditioner. It can be seen that with the new version of this method, GMRES also exhibits mesh-independent performance. (Note that for the step, performance improves as the grid is refined, until the convergence rate appears to settle to a constant for grid parameters l = 6 and 7.) In contrast, for previous versions, the asymptotic convergence factor degrades slightly as the grid is refined.

Finally, Figure 7.6 compares the performance of both the new preconditioners for two sets of problems with fixed grid parameter l = 6 and a variety of Reynolds numbers. The iteration counts for the LSC preconditioner are somewhat lower than those for the PCD preconditioner (although the former requires one more Poisson solve). The trends with respect to Reynolds number are largely the same for the two methods.

8. Conclusions. We have analyzed the role of boundary conditions within the pressure convection-diffusion preconditioner for saddle point systems arising from the incompressible Navier-Stokes equations. This examination has led to alternative formulations for the pressure convection-diffusion preconditioner and the LSC preconditioner that in effect emphasize certain commuting relationships near the boundary. Computationally, the number of required GMRES iterations is noticeably better than with the original versions of these preconditioners on two model benchmark problems. Further, the measured GMRES convergence rate with the new preconditioners now appears to be independent of the mesh resolution.

[FIGURE 7.6 OMITTED]

Acknowledgment. We thank David Silvester for a careful reading of this paper and several helpful remarks.

REFERENCES

[1] M. BENZI, G. GOLUB, AND J. LIESEN, Numerical solution of saddle point problems, Acta Numer., 14 (2005), pp. 1-137.

[2] H. ELMAN, V. E. HOWLE, J. SHADID, D. SILVESTER, AND R. TUMINARO, Least squares preconditioners for stabilized discretizations of the Navier-Stokes equations, SIAM J. Sci. Comput., 30 (2007), pp. 290-311.

[3] H. C. ELMAN, V. E. HOWLE, J. SHADID, R. SHUTTLEWORTH, AND R. TUMINARO, Block preconditioners based on approximate commutators, SIAM J. Sci. Comput., 27 (2006), pp. 1651-1668.

[4] H. C. ELM AN, D. LOGHIN, AND A.J. WATHEN, Preconditioning techniques for Newton's method for the incompressible Navier-Stokes equations, BIT, 43 (2003), pp. 961-974.

[5] H. C. ELM AN, D.J. SILVESTER, AND A. J. WATHEN, Finite Elements and Fast Iterative Solvers, Oxford University Press, Oxford, 2005.

[6] M. D. GUNZBURGER, Finite Element Methods for Viscous Incompressible Flows, Academic Press, San Diego, 1989.

[7] F. H. HARLOW AND J. E. WELCH, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids, 8 (1965), pp. 2182-2189.

[8] 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.

[9] D. A. MAY AND L. MORESI, Preconditioned iterative methods for Stokes flow problems arising in computational geodynamics, Phys. Earth Planet. Inter., 171 (2008), pp. 33-47.

[10] M. F. MURPHY, G. H. GOLUB, AND A. J. WATHEN, A note on preconditioning for indefinite linear systems, SIAM J. Sci. Comput., 21 (2000), pp. 1969-1972.

[11] M. A. OLSHANSKII AND Y. V. VASSILEVSKI, Pressure Schur complement preconditioners for the discrete Oseenproblem, SIAM J. Sci. Comput., 29 (2007), pp. 2686-2704.

[12] D. SILVESTER, H. ELMAN, D. KAY, AND A. WATHEN, Efficient preconditioning of the linearized Navier-Stokes equations for incompressible flow, J. Comput. Appl. Math., 128 (2001), pp. 261-279.

[13] D. SILVESTER AND A. WATHEN, Fast iterative solution of stabilised Stokes systems. Part II: Using general block preconditioners, SIAM J. Numer. Anal., 31 (1994), pp. 1352-1367.

[14] D. J. SILVESTER, H. C. ELMAN, AND A. RAMAGE, IFISS: Incompressible Flow Iterative Solution Software. Available at http : //www.manchester.ac.uk/ifiss.

[15] A. TOSELLI AND O. WIDLUND, Domain Decomposition Methods--Algorithms and Theory, Springer Series in Computational Mathematics, Springer, New York, 2004.

[16] R. VERFURTH, A combined conjugate gradient-multigrid algorithm for the numerical solution of the Stokes problem, IMA J. Numer. Anal., 4 (1984), pp. 441-455.

[17] A.J. WATHEN, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457.

(1) The (1, 1)-entry of [F.sub.2] is larger than the (1, 1)-entry of [F.sub.1] due to the fact that the leftmost grid point "[cross product]" is closer to the boundary than the leftmost grid point "x."

(2) By taking d = 1, l = -1, and q = 0, [B.sub.y] can be viewed as a special case of the structure for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and [F.sup.(p).sub.y].

(3) All perturbation rows associated with points adjacent to the inflow boundary have the same values due to the periodic boundary conditions. The only difference is that grid functions associated with different rows are shifted so that their peak corresponds to the row location.

(4) The only differences come from specific aspects of the MAC discretization, as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is defined on vertical edges whereas [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is defined on horizontal edges.

(5) From this discussion, one might conclude that [epsilon] = 0 is an appropriate choice, although it is evident that in this case H is singular. We have not found performance to be overly sensitive to [epsilon] and we have fixed it to be [epsilon] = .1 in experiments described below.

HOWARD C. ELMAN ([dagger]) AND RAY S. TUMINARO ([double dagger])

* Received April 3, 2009. Accepted for publication November 2, 2009. Published online December 31, 2009. Recommended by M. Benzi.

([dagger]) Department of Computer Science and Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742 (elman@cs.umd.edu). The work of this author was supported by the U. S. Department of Energy under grant DEFG0204ER25619 and by the U. S. National Science Foundation under grant CCF0726017.

([double dagger]) Sandia National Laboratories, PO Box 969, MS 9159, Livermore, CA 94551 (rstumin@sandia.gov). This author was supported in part by the DOE Office of Science ASCR Applied Math Research program. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

Table 7.1 Iteration counts for different combinations of operator orderings and boundary conditions. Backward facing step, Re = 200. PCD Old order Old order New order l Old b.c. New b.c. New b.c. 4 46 53 30 5 42 56 25 6 47 72 28 7 63 96 32 LSC Old order Old order New order l Old b.c. New b.c New b.c. 4 29 34 34 5 22 25 25 6 22 22 17 7 32 20 15

Printer friendly Cite/link Email Feedback | |

Author: | Elman, Howard C.; Tuminaro, Ray S. |
---|---|

Publication: | Electronic Transactions on Numerical Analysis |

Article Type: | Report |

Geographic Code: | 1USA |

Date: | Jan 1, 2009 |

Words: | 11453 |

Previous Article: | An efficient generalization of the rush-Larsen method for solving electro-physiology membrane equations. |

Next Article: | Convergence of a lattice numerical method for a boundary-value problem with free boundary and nonlinear Neumann boundary conditions. |

Topics: |