Fast solution of incompressible flow problems with two-level pressure approximation

This paper develops efficient preconditioned iterative solvers for incompressible flow problems discretised by an enriched Taylor-Hood mixed approximation, in which the usual pressure space is augmented by a piecewise constant pressure to ensure local mass conservation. This enrichment process causes over-specification of the pressure when the pressure space is defined by the union of standard Taylor-Hood basis functions and piecewise constant pressure basis functions, which complicates the design and implementation of efficient solvers for the resulting linear systems. We first describe the impact of this choice of pressure space specification on the matrices involved. Next, we show how to recover effective solvers for Stokes problems, with preconditioners based on the singular pressure mass matrix, and for Oseen systems arising from linearised Navier-Stokes equations, by using a two-stage pressure convection-diffusion strategy. The codes used to generate the numerical results are available online.


Introduction
Reliable and efficient iterative solvers for models of steady incompressible flow emerged in the early 1990s.Strategies based on block preconditioning of the underlying matrix operators using (algebraic or geometric) multigrid components have proved to be the key to realising mesh independent convergence (and optimal complexity) without the need for tuning parameters, particularly in the context of classical mixed finite element approximation, see Elman et al. [8, chap. 9].The focus of this contribution is on efficient solver strategies in cases where (an inf-sup) stable Taylor-Hood mixed approximation is augmented by a piecewise constant pressure in order to guarantee local conservation of mass.The augmentation leads to over-specification of the pressure solution when the pressure space is defined in the natural way via the union of Taylor-Hood pressure basis functions and basis functions for the piecewise constant pressure space, requiring a redesign of the established solver technology.
The idea of adding a piecewise constant pressure to the standard rectangular biquadratic velocity, bilinear pressure (Q 2 -Q 1 ) approximation was originally suggested during discussion around a blackboard at a conference on finite elements in fluids held in Banff in 1980; see Gresho et al. [9].The need for local mass conservation was motivated by competition from finite volume methods (such as the MAC scheme) in the design of effective strategies for modelling buoyancy-driven flow in the atmosphere.
The extension of the augmentation idea to Taylor-Hood (P 2 -P 1 ) triangular approximation was proposed in a paper presented at a conference in Reading in 1982; see Griffiths [10].A proof of stability of the augmented P 2 -P * 1 approximation on triangular meshes was constructed by Thatcher & Silvester [24] in 1987.An extended version of this manuscript included discussion of Q 2 -Q * 1 hexahedral elements [23].A rigorous assessment of the augmentation strategy was undertaken by Boffi et al. two decades later [4].The strategy of augmenting a continuous pressure approximation to give a locally mass-conserving strategy can also be generalised to higher-order Q k -Q k−1 and P k -P k−1 Taylor-Hood approximations.Inf-sup stability is assured for k ≥ 2 in two dimensions and for P k -P k−1 with k ≥ 3 in three dimensions; see [3, p. 506].
A closely related, yet fundamentally different Stokes approximation strategy is developed in [27].Their idea is to stabilise the lowest-order P 1 -P 0 approximation space defined on triangles or tetrahedra by augmenting the continuous velocity approximation by a piecewise constant approximation.The extended mixed approximation is stable but nonconforming.A similar approach is taken in [6], but an alternative weak form is used.Note that the design of linear solvers for the resulting discrete equations is relatively straightforward.A special feature of the extended approximation is that it can be augmented by postprocessing to give a pressure-robust approximation; see [13].
More generally, such an extended Galerkin (EG) approximation can be viewed as an intermediate between continuous Galerkin (CG) and discontinuous Galerkin (DG).This interpretation of piecewise constant augmentation was originally put forward by Sun & Liu [22] and is motivated by the fact that it gives local flux conservation when modelling transport in porous media flow problems, but with fewer degrees of freedom compared to vanilla DG.
The novel contribution in this work lies in the linear algebra aspects of two-field pressure approximation.The immediate issue that needs to be dealt with is the fact that, when the pressure space is specified by the frame formed by adding usual finite element basis functions for the underlying Taylor-Hood space to basis functions for the discontinuous pressure space, the mass matrix that determines the stability of the resulting mixed approximation is singular.This aspect is noted in the context of EG approximation in [15,Remark 4.1] and is discussed in section 2. The main issue that has to be dealt with in practical flow simulation is the over-specification and associated ill-conditioning of the discrete operators that arise in the preconditioning of the linearised Navier-Stokes operator.This is the focus of the discussion in section 3. The conclusion of this study is that optimal complexity convergence rates can be recovered when using two-field pressure approximation, but only after a careful redesign of the preconditioning components.

Two-field pressure mass matrix
In this section we consider the Stokes problem where u and p are the fluid velocity and pressure, respectively, and , is a polygonal or polyhedral domain.
Throughout this section, we assume that are an inf-sup stable pair of finite element spaces.Then the corresponding finite element approximation problem is to find where (•, •) denotes the usual L 2 inner product, right-hand side functions f h and g h are associated with the specified boundary velocity (g h is zero for enclosed flow), and Solving the finite element problem (1) is then equivalent to solving the linear system where A ∈ R nu×nu is symmetric positive definite and B ∈ R np×nu [8, chap. 3].The matrix A is of well-known saddle point type, and solvers for this sort of system have been extensively studied, see, e.g., Benzi et al. [1].Since A is large and sparse, the system is typically solved by an iterative method, with preconditioned MINRES [18] a popular choice.
An ideal preconditioner for A is [14,17] since the eigenvalues of the preconditioned matrix are 0, 1, and (1 ± √ 5)/2, with the zero eigenvalue appearing only if the preconditioned matrix is singular.This preconditioner is usually too costly to apply, and so efficient block diagonal preconditioners for A are typically based on spectrally equivalent preconditioners for A and the negative Schur complement S = BA −1 B T .For Stokes problems, the solve with A can be replaced by, e.g., an algebraic or geometric multigrid solver, while for stable discretisations S is spectrally equivalent to M Q , the pressure mass matrix [8, chap. 3], i.e., there exist constants γ and Γ, independent of h, such that for all vectors q except those corresponding to the function q h ≡ 1 on Ω.For certain element pairs, M Q itself is easily inverted, e.g., for discontinuous P 0 or Q 0 pressures, the mass matrix is diagonal.Otherwise, M Q can be replaced by its diagonal or by a fixed number of steps of Chebyshev semi-iteration acceleration of a Jacobi iteration when, as is common, diag(M Q ) −1 M Q has eigenvalues within a small interval and this interval lies away from the origin [25,26].
To summarise, for inf-sup stable finite element pairs, an effective preconditioner for (2) is where M ∈ R nu×nu is A or an approximation, and M S ∈ R np×np is M Q or an approximation.(See, e.g., Elman et al. [8,chap. 4] for results with Q 2 -Q 1 approximation on quadrilaterals.)In this section, our aim is to determine effective preconditioners for the enriched Taylor-Hood element.As we will see, although enriching the pressure space results in better mass conservation properties, it can also present challenges for solving (2).
2.1.Augmented Taylor-Hood Elements.We see from (1) that the mass conservation condition ∇ • u = 0 is imposed only in a weak sense, and if Taylor-Hood elements are employed we can only guarantee that (1b) will hold.However, by augmenting the pressure space by piecewise constant pressures it is possible to obtain local mass conservation, so that the average of the divergence is zero on each individual element.The natural choice is to define this space by a frame consisting of the union of the continuous Taylor-Hood pressure basis functions and the discontinuous pressure basis functions.In this case, as we will discuss later in this section, constant pressures have multiple representations, and this results in certain challenges for solving (2).Although it is certainly possible to instead compute a basis for this augmented pressure space, we show that the linear algebra challenges associated with using a frame can be overcome, which simplifies the implementation of the method.
Let us now describe the enriched Taylor-Hood finite element spaces.We first introduce a shape-regular family of simplicial, quadrilateral (in 2D) or hexahedral (in 3D) decompositions of the domain Ω.We assume that any two elements have at most a common face, edge, or vertex and denote by h the maximum diameter of any element.The total number of elements in the resulting mesh is n el .
We denote the usual Taylor-Hood finite element space by In the latter case, we additionally assume that the polynomial degree, k, satisfies k ≥ d − 1.The corresponding enriched Taylor-Hood space is and Q 0 h is the space of discontinuous pressures that are constant on each element.Thus, we see that the velocity approximation space is identical to that of the corresponding Taylor-Hood element, while augmented with piecewise constant pressures.It follows that functions in Q ⋆ h may be discontinuous across inter-element boundaries.We stress that the enriched Taylor-Hood space is well defined, and inf-sup stable.From a linear algebra perspective, however, a critical point is the representation of Q ⋆ h .Ideally, we would like to represent functions in Q ⋆ h as linear combinations of basis functions.In this case the resulting linear system would be nonsingular, and it is likely that the approaches described at the stat of this section could be applied directly.However, the most natural choice is define Q ⋆ h by a frame, i.e., to let where {φ k } n k k=1 and {φ k } np k=n k +1 are Lagrange bases of Q T H h and Q 0 h (see (5)).This approach to specifying the pressure space has been used in, e.g., [4,20,23].
Although this makes specification of the pressure straightforward, with this choice any constant function on Ω can be represented by a function in Q T H h or a function in Q 0 h .This has profound consequences for the linear algebra: the pressure mass matrix M Q that determines M S in (4) becomes singular, and the rank of the matrix B T is reduced.In the rest of this section we first establish these properties, before showing that it is still possible to solve (2) by preconditioned MINRES in this case.
Using (6), we can relate vectors p ∈ R np and functions p = p k + p 0 ∈ Q ⋆ h , where As mentioned above, any constant function has multiple representations when specified via (6).In particular, p = p k + p 0 ≡ 0 with p k = α and p 0 = −α, for any α ∈ R. We see from (7) that this representation of the zero function corresponds to the vector αk, where and n 0 = n p − n k .A direct consequence of the correspondence between k and the zero function on Ω is that M Q k = 0 and B T k = 0, as we now show.
be the pressure mass matrix for the enriched Taylor-Hood pressure space Q ⋆ h in (5), specified by the frame (6).Then, null(M Q ) = span{k}, where k is given in (8).
Proof.Let p ∈ R np , p = 0.Then, using (7), we find that which implies that p ≡ 0 in Ω since p is continuous on each element.Since p ≡ 0 corresponds to vectors of the form αk, we find that span{k} ⊆ null(M Q ).
We now show that there are no other vectors in the nullspace.Since p = p k + p 0 ≡ 0, But p 0 is piecewise constant, and p k is continuous on Ω.It follows that p 0 ≡ α and p k ≡ −α on Ω for some constant α ∈ R. Such functions correspond to vectors of the form αk, which shows that null(M Q ) = span{k}, as required.

Remark 2. A very similar argument shows that k ∈ null(B T ).
What do these results mean for the solution of (2) by preconditioned MINRES, when the coefficient matrix is constructed using the frame in (6)?First, if B T k = 0 then it follows that A is always singular, with Aw = 0, where If the linear system (2) is consistent then this does not pose a problem for preconditioned MINRES.However, as we shall now see, the proposed block diagonal preconditioner (4), with M S = M Q , is also singular, and so effective preconditioning requires some care.

9:
Solve Pz (j+1) = v (j+1) 10: 11: : 13: x (j) = x (j−1) + c j+1 ηw (j+1) 18: <Test for convergence> 20: end for 2.2.Preconditioning considerations.Knowing that the enriched Taylor-Hood element is inf-sup stable implies that the preconditioner (4) will be effective when solving (2).However, with the common choice of frame (6), the matrix M Q , and hence the preconditioner P, are singular, since Pw = 0, where w is given in (9).(Note that this implies that A and P have a common nullspace.)Here, we will show that using a singular preconditioner causes no difficulty for the preconditioned MINRES method in Algorithm 3 in exact arithmetic.At the end of section 2, we present numerical results with two different preconditioners to demonstrate our approach.
We start by noting that at each iteration step we need to solve a linear system of the form Pz (j) = v (j) .If this system is consistent there are infinitely many solutions, which take the form z (j) = z (j) +ζ j w, with z (j) ⊥ w and ζ j ∈ R. Our next step is to examine the effect of |ζ j | on the scalars and vectors in Algorithm 3. Observe that if the linear system Ax = b is consistent then b ⊥ w and so v (1) ⊥ w.It then follows by induction that v (j) ⊥ w, j = 1, 2, . . ., and that the systems Pz (j) = v (j) are all consistent.Furthermore, because of the orthogonality of v (j) and w, so that ζ j does not affect γ j .Similarly, δ j = Az (j) , z (j) = Az (j) , z (j) , which shows that δ j is similarly unaffected by ζ j .Indeed, the only quantities that are affected by the nullspace components ζ j w, j = 1, 2, . . ., are the vectors w (j) and x (j) .In exact arithmetic, this is not a problem, since solutions of Ax = b may certainly contain a component in the direction of w.However, in unlucky cases, the size of the nullspace component of x (j) may be so large as to dominate the approximate solution.Alternatively, in finite precision v (j) and w may not be exactly orthogonal.Hence, it may be wise to explicitly ensure that z ⊥ w.One option is to orthogonalise z against w after each preconditioner solve, but if |ζ j | is large then the result may be inaccurate.A more robust approach is to note that solutions of Pz = r are minimisers of the quadratic form 1 2 z T Pz − z T r, since P is positive semidefinite.Constraining z to be orthogonal to w is then equivalent to the following optimisation problem: Applying a Lagrange multiplier approach results in the augmented system where λ is the Lagrange multiplier, and solving this augmented system gives a vector z that is orthogonal to w.Moreover, since we see that only the solve with M S needs to be modified.

2.3.
Approximating the two-level pressure mass matrix.Now let us consider approximations of the matrix M Q which, because of the two-level pressure approximation, has 2 × 2 block structure: where Note that Q k is the standard Taylor-Hood pressure mass matrix, and Q 0 is the standard discontinuous pressure mass matrix, while R represents cross terms between the spaces Q T H h and Q 0 h (see (5)).Naively, one may wish to approximate Q k by its diagonal, but this results in slow convergence rates; it can be shown, using the method described by Wathen [26], that diag(Q k ) −1 Q k has eigenvalues in an interval [0, µ] where, for example, for P 1 elements on triangles µ = 3. Accordingly, applying a fixed number of iterations of a Chebyshev semi-iteration based on a Jacobi iteration, as is standard for nonsingular Galerkin mass matrices [25], leads to large iteration counts here.Replacing Q k or Q 0 in M Q by their diagonals results in similarly poor performance.
However, it is possible to replace M Q by an approximation designed for symmetric positive semidefinite matrices, see, e.g., [5,7,16] and the references therein.For example, symmetric Gauss-Seidel is convergent for such matrices [7].Here, we have applied a fixed number of iterations of Chebyshev semi-iteration based on symmetric Gauss-Seidel [5].
The results are still mesh-dependent, because the largest non-unit eigenvalue of the underlying symmetric Gauss-Seidel iteration matrix approaches 1 as the mesh is refined but, for large problems, the cost per iteration is much lower than applying M Q exactly.
2.4.Reliable computation of the discrete inf-sup constant.The inf-sup constant depends on the shape of the domain Ω.Thus the constant for a step domain with a long channel is much smaller than the constant for a square domain.Estimation of the infsup constant dynamically provides useful information regarding the connection between the velocity error and the pressure error.The inf-sup constant also features in the norm equivalence between the residual norm of the saddle-point system and the natural "energy" norm ∇u L 2 (Ω) + p L 2 0 (Ω) , as discussed in the motivating paper [21].A typical strategy to estimate the inf-sup constant for (1) is to find the largest γ that satisfies (3), i.e., to find the smallest nonzero eigenvalue of the generalised eigenvalue problem where A is prescribed by the velocity basis functions, M Q by the pressure functions, and B by both the velocity and pressure functions.
Although the inf-sup constant is independent of the choice of these functions, in practice, its computation may be affected.For example, with the common choice to define Q ⋆ h by the frame (6), Proposition 1 and Remark 2 show that k lies in the nullspaces of both B T and M Q , which means that this generalised eigenvalue problem is singular, i.e., any λ ∈ R satisfies (10) when v = k.It is known that generalised eigenvalue problems with singular pencils are challenging to solve numerically [11], and additional checks must be performed to ensure that an estimate of γ is not associated with the eigenvector k.In practice we find that standard methods for computing eigenvalues of sparse matrices may struggle to accurately compute these eigenvalues, precisely because BA −1 B T and M Q are singular.
An alternative is to define Q ⋆ h by a basis instead of the commonly-used frame (6), but it turns out that this is not necessary.The EST-MINRES approach proposed in [21] is much more robust for this singular eigenvalue problem and gives consistently reliable results for certain preconditioners.The intuition is that by ensuring that any solves with M Q are orthogonal to k within the preconditioned MINRES method, (see Section 2.2), we instead solve (10) for v ⊥ k.
We illustrate the EST-MINRES inf-sup constant estimates using two representative test problems that we describe below.All numerical results were obtained using T-IFISS [2] (for triangular elements) and IFISS3D [19] (for cubic elements).The stopping criterion for preconditioned MINRES is a reduction of the norm of the preconditioned residual by eight orders of magnitude, i.e., r k P −1 / r 0 P −1 < 10 −8 .We apply two different block diagonal preconditioners (4).The first is Although this preconditioner is expensive to apply, because it involves exact solves with A, we have used it here to more clearly illustrate the key findings of this section, namely, that (4) with M S = M Q is an effective preconditioner for augmented Taylor-Hood problems, despite the singularity of M Q , and that EST-MINRES provides reliable approximations to γ, the discrete inf-sup constant in (3).However, we note that it is possible to replace A by, e.g., an algebraic multigrid method such as the HSL code MI20 [12].With this AMG approximation, we are able to solve a 3D problem with nearly 10 6 degrees of freedom in under 10 seconds on a quad-core, 62 GB RAM, Intel i7-6700 CPU with 3.20GHz.Moreover, the effect on the inf-sup constant is fairly small.
Our second, cheaper, preconditioner is where A AMG is one AMG V-cycle, with the default parameters in T-IFISS and IFISS3D and M cheb is 20 iterations of a Chebyshev semi-iteration method based on symmetric Gauss-Seidel.
Test problem 1 (two-dimensional enclosed flow).Our first example is a classical drivencavity flow in the square domain D = [−1, 1] 2 .A Dirichlet no-flow condition is imposed on the bottom and side boundaries, while on the lid the nonzero tangential velocity is u y = 1 − x 4 .The domain is subdivided uniformly into n 2 bisected squares.We use the standard P 2 -P 1 Taylor-Hood mixed approximation and the augmented P 2 -P * 1 approximation.The two components of the pressure solution for the P 2 -P * 1 approximation, computed for n = 32, are illustrated in Fig. 1.The centroid pressure field is concentrated in the two corners where the pressure is singular, and the centroid pressures are an order of magnitude smaller than the vertex pressure in all elements.Consequently, the overall pressure field is visually identical to the P 1 pressure field shown in the left plot in Fig. 1.
Test problem 2 (three-dimensional enclosed flow).Our second problem is a threedimensional version of driven-cavity flow.The domain is now D = [−1, 1] 3 .As in the previous example, the flow is enclosed, but now the nonzero tangential velocity u y = (1 − x 4 )(1 − z 4 ) is specified on the top of the cavity.The domain is subdivided uniformly into n 3 cubic elements, and we use Table 1 shows preconditioned MINRES iteration counts and EST-MINRES discrete inf-sup constant approximations for Example 1 with preconditioner P 1 .We first note that for both the P 2 -P 1 and P 2 -P * 1 approximations the iteration counts are quite similar and are mesh independent.In both cases the inf-sup constant approximations also appear to be converging from above.The approximation for P 2 -P * 1 elements appears to rapidly converge to four digits, indicating that even a relatively coarse grid is sufficient to obtain an approximation to the discrete inf-sup constant.However, the approximation for P 2 -P 1 elements appears to converge more slowly.We also note that the two approaches give different inf-sup constant estimates, at least for the grids shown here.This is not so surprising as the matrices in (3) depend on the choice of finite element spaces.
When we replace the preconditioner by the cheaper P 2 , the situation is largely unchanged for P 2 -P 1 elements (see Table 2).However, for P 2 -P * 1 elements, the preconditioned MINRES convergence rate and inf-sup constant estimate degrade as the mesh is refined.The effect on the convergence rate is caused by the largest non-unit eigenvalue of the symmetric Gauss-Seidel iteration matrix approaching 1 as the mesh is refined; this appears to also impact inf-sup constant approximation.Fig. 2 plots the inf-sup approximations at each iteration of preconditioned MINRES for Example 1 for preconditioner P 1 .We see that a good approximation of the inf-sup constant is obtained after 20-25 iterations.It is again clear that for the enriched Taylor-Hood approximations we obtain very similar approximations for all grids.
The results for Example 2 are given in Tables 3 and 4 and Fig. 3. Preconditioned MINRES iteration counts for P 1 are again mesh independent and, for all but the coarsest mesh, are almost identical for the two element pairs.The discrete inf-sup constant approximations again appear to converge from above.Now, at least for the grids presented EST-MINRES estimates of the discrete inf-sup constant γ 2 at each iteration for Example 1 and preconditioner P 1 with P 2 -P 1 (left) and P 2 -P * 1 (right) elements for the grids specified in Table 1.here, the estimates of γ 2 for Taylor-Hood elements seem to "track" the augmented Taylor-Hood estimates, i.e., the Taylor-Hood approximation on grid j is almost the same as the augmented Taylor-Hood estimate on grid j − 1. Table 4 shows that using the cheaper P 2 also results in mesh-independent convergence for Q 2 -Q 1 elements.Additionally, the inf-sup constant approximations seem to "track" the Taylor-Hood estimates in Table 3.Now, however, there is more modest growth in the iteration counts when this approximate preconditioner is used for Q 2 -Q * 1 elements and the inf-sup approximations are much closer to those in Table 3.
We also see from Fig. 2, which plots the inf-sup approximations at each iteration of preconditioned MINRES, for preconditioner P 1 , that the discrete inf-sup constant is already reasonably well approximated after 25 iterations.
elements for the grids specified in Table 3.

Two-field pressure preconditioning strategies for Oseen flow.
In this section we consider the Oseen problem that arises from applying fixed point iteration to the Navier-Stokes equations modelling steady flow in a channel domain Ω ⊂ R 2 with viscosity parameter ν.The Oseen problem is a linear elliptic system of PDEs and assuming a divergence-free convection field w and a smooth (e.g., polygonal) domain Ω, has a unique weak solution for all positive values of the viscosity parameter. 1The inf-sup stability of our two-field mixed approximation is a sufficient condition for convergence of the mixed approximation to the weak solution of the Oseen problem as h → 0. Our focus will be on the associated discrete matrix system, where the unknown coefficient vector involves the discrete velocity vector u ∈ R nu and the two-field pressure vector p ∈ R np .The right-hand side vectors f and g are associated with the specified boundary velocity (as in the Stokes flow case).References to the "residual norm" throughout this section will refer to the standard Euclidean norm for finite dimensional vectors.Thus, for the system (11) the convergence of our iterative solver strategies is assessed by monitoring The nonsymmetry of the matrix F means that the iterative solver of choice is GMRES (see [8, section 9.1]) together with a block preconditioning operator of the form where M is an optimal complexity (multigrid) operator effecting the action of the inverse of the matrix F and M S is an optimal complexity approximation of the Schur complement matrix BF −1 B T .We will discuss results for two representative flow problems herein.resulting system (11) is singular with the one-dimensional pressure nullspace described in section 2. The two components of the pressure solution computed on a representative P 2 -P * 1 mesh are illustrated in Fig. 4. The centroid pressure values are an order of magnitude smaller than the vertex pressure in all the elements-they provide the "corrections" to the vertex pressures that are needed to ensure local (elementwise) conservation of mass.
Test problem 4 (two-dimensional enclosed flow).We consider the classical driven-cavity enclosed flow problem defined on the domain D = [−1, 1] 2 with the viscosity parameter ν set to 1/100 and a nonzero tangential velocity u y = 1 − x 4 specified on the top of the cavity.We take P 2 -P * 1 augmented Taylor-Hood mixed approximation with the domain subdivided uniformly into n 2 bisected squares.We consider the discrete system (11) that arises after 5 fixed-point iterations starting from the corresponding Stokes flow solution.The discrete system is singular with a two-dimensional pressure nullspace, corresponding to a constant vertex pressure and a constant centroid pressure.The contour plot on the left in Fig. 5 shows the element divergence errors ∇ • u h L 2 (△) computed on a coarse mesh (n = 32).The associated centroid "correction" pressure field contours shown on the plot in the right can be seen to provide a good indication of the regions where the divergence error is concentrated.As in Fig. 4 the centroid pressure values are an order of magnitude smaller than the vertex pressure in all elements, so the overall pressure field is visually identical to the P 1 pressure field (not shown here).

Pressure convection-diffusion preconditioning for Oseen flow.
There are two alternative ways of approximating the key matrix BF −1 B T in the case that B T is generated by a two-field pressure approximation so that . The focus will be on pressure convection-diffusion (PCD) preconditioning in this section.Both of the Schur complement approximations can be motivated by starting with the Oseen matrix operator (11) and observing that the diagonal blocks of F are discrete representations of the convection-diffusion operator defined on the velocity space.In practical calculations ν > 0 is proportional to the inverse of the flow Reynolds number and w h is the discrete approximation to the flow velocity computed at the most recent nonlinear iteration.The PCD approximation supposes that there is an analogous operator to (14), namely defined on the two components of the augmented pressure space.
To this end, defining {φ j } n 1 j=1 to be the basis for the C 0 pressure discretisation, we construct matrices Q 1 and F 1 so that We then note that if the commutator with the divergence operator is small then we have the approximation where M is the diagonal of the mass matrix associated with the basis representation of the velocity space.2Rearranging ( 16) gives the first Schur complement approximation A discrete version of L p for the piecewise constant pressure space can be generated by considering the jumps in pressure across inter-element boundaries; see [8, pp. 268-370].To this end, defining {ϕ j } n 0 j=1 to be the (indicator function) basis for the discontinuous pressure, we construct matrices Q 0 and F 0 via and note that if the commutator with the divergence operator is small then we have suggesting the second Schur complement approximation Combining ( 17) with ( 19) then gives a two-field PCD approximation Two features of the PCD approximation (20) are worth noting.The first point is that the coupling between the pressure components is represented by the 2 × 2 block matrix BM −1 B T rather than by the pressure mass matrix M Q .The second key point is that the matrices M S and BF −1 B T have the same nullspace, independent of the nature of underlying flow problem that is being solved.The PCD approximation in (20) is imperfect in practice.To illustrate this, representative convergence histories that arise in solving the inflow-outflow problem using P 2 -P * 1 approximation (shown in Fig. 4) are presented in Fig. 6.Taking the solution from the previous Picard iteration as the initial guess, we note that the initial residual norm of the target matrix system ( 11) is close to 10 −4 independent of the spatial discretisation.Convergence plots are presented for two preconditioning strategies, namely with M S defined in (20).The first strategy is the block triangular extension of the Stokes preconditioning strategy discussed in section 2. We note that the resulting convergence is very slow-around 50 iterations are required to reduce the residual norm by an order of magnitude-but is independent of the discretisation level.In contrast we see that the PCD preconditioning strategy has two distinctive phases of convergence behaviour.An initial phase of relatively fast convergence is followed by a secondary phase where GMRES stagnates.We hypothesise that this stagnation is a consequence of the ill-conditioning of the eigenvectors of the matrix operator M S .A notable feature is that the onset of the stagnation is delayed when solving the same problem on a finer grid.The two-phase convergence behaviour is ubiquitous-the same pattern is seen using rectangular elements and the stagnation does not go away when the viscosity parameter is increased from 1/50 to 1/5.An alternative strategy is clearly needed!One way of designing a more robust PCD preconditioning strategy for a two-field pressure approximation is to exploit the fast convergence of the PCD approximation for the unaugmented Taylor-Hood approximation.The starting point for such a strategy is to rewrite the system (11) in the form  The proposed solution algorithm is then a two-stage process.
• Input: residual reduction tolerance η and reduced system residual factor c Step I Generate a PCD solution to the reduced system using the Schur complement approximation (16), stopping the GMRES iteration when the residual is reduced by a factor of cη.
Step II Generate a solution to the target system (22) with residual tolerance η using preconditioning strategy M 1 in ( 21) with the refined initial guess  Sample results generated using this strategy with tolerance η set to 10 −4 and c set to 10 are presented in Fig. 7. Sample results for the cavity test problem are presented in Fig. 8.The results in Fig. 7 and Fig. 8 are representative of the two-stage convergence profiles that are generated when solving these test problems at other Reynolds numbers.Our experience is that the level of residual reduction is perfectly robust with regards to the spatial discretisation-typically giving smaller iteration counts when the mesh resolution is increased (a known feature of PCD preconditioning).The convergence rates of both stages of the algorithm deteriorate slowly when the Reynolds number is increased.Our strategy for terminating the first stage of the iteration is motivated by the following result.1 , q * 1 , 0] to the discrete system (22) satisfies the bound where f is the initial residual error associated with a zero initial vector.
Proof.The vector z * associated with the intermediate solution is the three-component vector The stopping test for solving the reduced system ensures that Thus we have The bound (24) has two terms on the right-hand side.While the first term can be controlled by reducing η and/or c, the second term measures the local incompressibility of the intermediate Taylor-Hood solution B 0 u * , where u * h is the expansion of the coefficient vector u * 1 in the basis of the velocity approximation space.Setting η = 10 −4 and c = 10 we see that the second term saturates the residual error and the residual error jumps up when the switch is made from the first to the second step of the algorithm.This "transition" jump in the residual norm is clearly evident in the convergence plots.The convergence in the second step is rapid initially but eventually mirrors the rate observed for M 1 approximation with a standard starting guess.

Least square commutator approximation for Oseen flow.
A second way of approximating the key matrix BF −1 B T in the case that B T is generated by a two-field pressure approximation is given by the least-squares commutator (LSC) preconditioner The attractive feature of LSC is that the construction of M S is completely algebraic.The only technical issue is the need to make adjustments on rows and columns associated with tangential velocity degree of freedom adjacent to inflow and fixed wall boundaries.These adjustments are associated with a diagonal scaling matrix D so that H = D −1/2 M D −1/2 .Full details can be found in [8, pp. 376-379]).
The LSC approximation (28) is also far from perfect when using triangular elements. 3o illustrate this, representative convergence histories that arise in solving the inflowoutflow problem using P 2 -P * 1 approximation are presented in Fig. 9. Convergence plots are presented for two preconditioning strategies, namely the refined PCD from the previous section and with M S defined in (28).The PCD histories in Fig. 9 display mesh independent convergence and are comparable with those observed using square elements in Fig. 7.The associated cpu times for the solution on the intermediate grid with 2×11264 elements were 13 seconds for the first step (28 iterations) and 26 seconds for the second step (53 iterations).The LSC preconditioning strategy is not robust-the convergence rate deteriorates with increasing grid refinement.The cpu time for generating the LSC solution on the intermediate grid was over 100 seconds (61 iterations).

Conclusions
Two-level pressure approximation for incompressible flow problems offer the prospect of accurate approximation with minimal computational overhead.Derived quantities of practical importance such as the mean wall shear stress are likely to be computed much more precisely if incompressibility is enforced locally. 4However, the augmented pressure space causes some challenges for the linear algebra when the pressure space is defined by combining the usual Taylor-Hood pressure basis functions with basis functions for the piecewise continuous pressure space.This is because constant functions can be expressed using either the usual (continuous) Taylor-Hood pressure space, or the augmented piecewise constant pressure space.Specifically, with this common choice for specifying the pressure space, the pressure mass matrix becomes singular, and care should be taken to carefully construct and apply preconditioners that involve this matrix.Care should also be exercised when approximating the discrete inf-sup constant, and we find that naive approaches are not always reliable.On the other hand, the approximation implemented in EST-MINRES is robust.Our computational experimentation indicates that our twostage PCD strategy could be the best way of iteratively solving two-level pressure discrete linear algebra systems in the sense of algorithmic reliability and computational efficiency.

Figure 1 .
Figure 1.Representative P 2 -P * 1 pressure field solution for the cavity flow in Example 1 computed on a uniform mesh with 2048 right-angled triangles and 1089 vertices.Contours are equally spaced between the maximum and minimum values in both plots.

Figure 2 .
Figure2.EST-MINRES estimates of the discrete inf-sup constant γ 2 at each iteration for Example 1 and preconditioner P 1 with P 2 -P 1 (left) and P 2 -P * 1 (right) elements for the grids specified in Table1.

Figure 3 .
Figure 3. EST-MINRES estimates of the discrete inf-sup constant γ 2 at each iteration for Example 2 and preconditionerP 1 with Q 2 -Q 1 (left) and Q 2 -Q * 1 (right)elements for the grids specified in Table3.

Figure 4 .Figure 5 .
Figure 4. Representative P 2 -P * 1 pressure field solution for flow over a step (R = 100) computed on a uniform mesh with 5632 right-angled triangles and 2945 vertices.

Figure 6 .
Figure 6.Absolute residual reduction for test problem 3 when computing P 2 -P * 1 solutions using preconditioners M 1 (x) or M 2 (o) on two nested meshes.

Figure 7 .
Figure 7. Absolute residual reduction for test problem 3 when computing Q 2 -Q * 1 solutions using using preconditioners M 1 (x) or two-stage PCD (Step I: o; Step II: o) on two nested meshes.

Figure 8 .
Figure 8. Absolute residual reduction for test problem 4 when computing P 2 -P * 1 solutions using using preconditioners M 1 (x) or two-stage PCD (Step I: o; Step II: o) on two nested meshes.

Proposition 4 .
The residual error z * associated with the intermediate solution [u *

Figure 9 .
Figure 9. Absolute residual reduction for test problem 3 when computing P 2 -P * 1 solutions using using refined PCD (left) or preconditioner M 3 (right) on three nested meshes.

Table 1 .
Preconditioned MINRES iterations and discrete inf-sup constant approximations for Example 1 and preconditioner P 1

Table 2 .
Preconditioned MINRES iterations and discrete inf-sup constant approximations for Example 1 and preconditioner P 2 .

Table 3 .
Preconditioned MINRES iterations and discrete inf-sup constant approximations for Example 2 and preconditioner P 1 .

Table 4 .
Preconditioned MINRES iterations and discrete inf-sup constant approximations for Example 2 and preconditioner P 2 .