On the initial higher-order pressure convergence in equal-order finite element discretizations of the Stokes system

The finite element discretization of pressure and velocity in incompressible flow systems can be done through either stable spaces or stabilized pairs. When using stabilized formulations with lowest-order piecewise linear spaces, the classical theory guarantees only linear convergence in L for the pressure approximation, although a higher order is often observed – but seldom discussed – in numerical practice. Such experimental observations may, in the absence of a sound a priori error analysis, mislead the selection of discretization spaces in practical applications. Therefore, we present herein a numerical analysis demonstrating that an initial higher-order pressure convergence may in fact occur under certain conditions. Moreover, our numerical experiments clearly indicate that whether and for how long this behaviour holds is a problem-dependent matter. Our findings confirm that an optimal pressure convergence can in general not be expected when using unbalanced velocity-pressure discretizations.


Introduction
The Stokes system of incompressible creeping flows is a mixed problem having a scalar pressure p and a velocity vector u as primal unknowns. The stability and convergence of finite element discretisations for this saddle-point problem are subject to the well-known Babuška-Brezzi theory [1,2], which is violated for instance when both u and p are approximated using the same polynomial degree. Unique solvability can be attained by using either stable pairs such as in Taylor-Hood [3] and MINI [4] elements, or stabilized equalorder discretisations. The main weakness of equal-order (as well as of MINI) elements is the suboptimal pressure convergence resulting from the unbalanced approximation properties of the velocity and pressure spaces [5]. When using piecewise linear elements, for example, the standard theory only guarantees first-order convergence for the pressure in L 2 [6], which is one order lower than the corresponding best approximation error estimate. Nevertheless, the order observed in numerical experiments [7,8,9,10,11,12,13,14,15] is often higher than one even for very fine meshes, which usually goes undiscussed or at least unexplained.
Around a decade ago, Eichel et al. [12] proved a half-an-order-higher pressure convergence for low-order elements in problems with very smooth solutions and uniform, orthogonal meshes. To the best of our knowledge, however, there is currently no available analysis on a possible higher-order pressure convergence for equal-order elements of arbitrary degree and for more general meshes. In fact, Cioncolini and Boffi [14] have very recently carried out an empirical numerical investigation on unstructured meshes to assess whether the assumptions in [12] might have been too strict. Despite providing insightful numerical evidence on the performance of the MINI element, the results in [14] do not allow a definitive conclusion regarding this apparent "superconvergence": in most cases a 3/2 slope was indeed verified, while in some examples there was an apparent degradation of this higher slope at finer levels.
In this context, we present new theoretical and numerical results on the pressure convergence of unbalanced discretisations of the Stokes system. Our theory includes equal-order pressure stabilization methods, as well as the MINI element. Developed by Arnold et al. [4], the MINI element provides one of the simplest inf-sup stable spaces for incompressible flow problems and consists of enriching a piecewise linear velocity space with bubble functions. For the Stokes system, the bubble degrees of freedom can be completely decoupled from the first-order part of the velocity, which results in a lowestorder discretisation with a simple matrix added to the pressure-pressure block (cf. [16]). This means that the Stokes system discretised with MINI elements can be considered as either a stable or a low-order stabilized formulation. For this reason, we include herein both MINI and equal-order elements under the umbrella of pressure-stabilized methods.
Using a Schur complement formulation we show, under standard mesh regularity assumptions, that an initial higher-order pressure convergence may take place depending on certain problem-and discretisation-dependent constants, similarly as observed in boundary element methods [17] when using equal-order elements approximating the Cauchy data on the boundary. In comparison to the existing theory [12], ours requires less regularity on the (exact) velocity, allows unstructured meshes and accommodates higher-order elements. Finally, we provide numerical examples demonstrating that our estimate is sharp: indeed such a superconvergence cannot be expected to hold undefinitely, eventually breaking down in most cases -even for structured, uniformly refined meshes.

The Stokes system
As a model problem we consider the Dirichlet boundary value problem for the Stokes system where Ω ⊂ R n , n = 2, 3, is a connected, bounded domain with Lipschitz boundary Γ = ∂Ω. The standard varational formulation of (1) is to find (u, p) ∈ H 1 0 (Ω) × L 2 (Ω) such that is satisfied for all (v, q) ∈ H 1 0 (Ω) × L 2 (Ω), which ensures the scaling condition Ω p dx = 0 for any solution of (2). When using the Riesz representation we can define linear bounded operators A : We also define Lp, q Ω := Ω p dx Ω q dx for all p, q ∈ L 2 (Ω).
Hence we can write the variational formulation (2) in operator form as Note that we have for all u, v ∈ H 1 0 (Ω), and Alternatively, we can consider a perturbed system by adding to L a stabilization operator to result in an invertible pressure-pressure operator D. This is the basis for stabilized formulations, where consistency terms can also be added to the right-hand side in (3) and to the pressure-velocity block B . In this setting, the system takes the more general form with D non-negative and C bounded: Since A is invertible, we can solve the first equation in (4) to get the Schur complement system In other words, we have the variational formulation to find p ∈ L 2 (Ω) such that From the properties of A, B, C and D, at the continuous level we immediately get that the operator S : L 2 (Ω) → L 2 (Ω) is bounded and elliptic, that is, Therefore, we conclude unique solvability of the variational problem (6), which will be the basis for deriving our theory.

Finite element error analysis
Let us assume a shape-regular triangulation of the domain Ω into simplicial elements Ω T , and two finite element spaces Π h × X h ⊂ L 2 (Ω) × H 1 0 (Ω) for the discretization of pressure and velocity. We denote by h T := n |Ω T | the size of Ω T , and by h := max{h T } the global mesh size. For a conforming finite element space Π h = S ν h (Ω) of piecewise polynomial basis functions of degree ν, we consider the Galerkin formulation to find Using standard arguments we arrive at Cea's lemma and from the approximation property of Π h we finally conclude the error estimate when assuming p ∈ H ν+1 (Ω).
Since the composed operator S = CA −1 B + D in general does not allow a direct evaluation, we construct a suitable approximation by defining, for any p ∈ L 2 (Ω), a vector w = A −1 Bp ∈ H 1 0 (Ω), which is the unique solution of the variational formulation Let X h := [S µ h (Ω)∩H 1 0 (Ω)] n be a second finite element space of polynomial basis functions with degree µ, for which we consider the Galerkin formulation to find w h ∈ X h such that From the ellipticity of A and the boundedness of B we have . Moreover, using standard arguments, we obtain the a priori error estimate when assuming w ∈ H µ+1 (Ω). Instead of we now define the approximate operator where we have Moreover,S : L 2 (Ω) → L 2 (Ω) is bounded: Let us assume thatS is elliptic in Π h , i.e., which is satisfied when using either inf-sup stable finite elements or appropriate stabilization operators. Then, we consider the perturbed variational formulation to findp h ∈ Π h such that We now recall the well-known Strang lemma.
Theorem 3.1. LetS be the approximate Schur complement operator as defined in (8), andp h the corresponding pressure approximation from (9). Then, under the assumptions of (7) and, additionally, the ellipticity ofS, there holds the error estimate with α and β independent of h.
Proof. From the triangle inequality and the error estimate (7) we have From the ellipticity ofS in Π h we conclude that is, which concludes the proof with α = (1 + c 4 /cS)c 1 and β = c 3 /cS .
The error estimate (10) implies the choice µ = ν + 1 to ensure an optimal order of convergence. On the other hand, the choice of equal-order elements, i.e., µ = ν, will asymptotically not result in an improved order of convergence. However, we are now in a position to show an initial higher-order convergence for the pressure approximation under certain conditions. From the triangle and Young's inequalities and the estimate (7), we As in the previous proof we have for µ = ν. We also have Thus, due to (12) we get

Let us write
So, as long as is satisfied for some γ < 1, we get , so that we can finally conclude This estimate provides an explanation for the higher (than one) order often observed for the pressure convergence in numerical practice, see e.g. Refs. [7,8,9,10,11,12,13,14,15]. Although most assumptions made towards proving (14) are rather standard, this is not the case for condition (13). In fact, how can we interpret such a condition? For small h, the expression on the left-hand side of (13) behaves (on a logarithmic scale) as a line with slope ν. If the right-hand side p −p h is assumed to be of order ν + 1/2 for some range of h, we see that even if condition (13) is satisfied initially, the two lines (left-and right-hand side) will eventually intersect as h decreases, and thus (13) will no longer hold. This is why estimate (14) shows only a possible initial higher-order convergence. As an alternative scenario, the curve with the higher slope may already start below the other curve, so that (13) will not hold even for the coarsest possible h, see Figure 1 for a graphical illustration. Besides, due to the several constants in (13), whether and for how long this higher-order convergence holds may depend on various factors such as the problem (domain and solution) and the discretisation. Also for this reason, the higher slope might not break down at all within a practical range of mesh sizes.

Pressure-based stabilization methods
As realization of the perturbed (stabilized) system (4), we now give three popular examples of pressure-stabilization methods, which we will also use in the numerical experiments. For simplicity of presentation, we will consider purely pressure-based stabilization, i.e., when C = B and g = 0 in (4). The pressure-pressure operator D can be expressed through Dp, q Ω := Ω p dx Ω q dx + s(p, q) for all p, q ∈ L 2 (Ω) , where s(·, ·) is a stabilizing bilinear form. The first and probably simplest stabilization method ever conceived is the pressure Poisson equation (PPE) by Brezzi and Pitkäranta [18], in which with the optimal parameter α = 1/12 for linear elements [5]. Another classical example is the polynomial pressure projection (PPP) method by Dohrmann and Bochev [9]: in which π ν−1 is an operator projecting p h ∈ S ν h (Ω) locally onto a space with reduced polynomial degree ν − 1 (refer to [9] for implementation details).
A slightly more complicated case is that of the MINI element, in which Π h = S 1 h (Ω) and X h is constructed by enriching a first-order velocity space n with a space X b h of standard bubble functions, for stabilization. The discretization of the variational problem (4) leads to a linear algebraic system with subscript h used for indicating the discrete counterparts of the respective operators/quantities in the infinite-dimensional case (4). It is simple to show that the spaces X 1 h (Ω) and X b h (Ω) are A-orthogonal, i.e., Thus, by splitting u h = u 1 h + u b h , we get the system which can be rewritten as That is, the bubble degrees of freedom can be easily eliminated (as A b h is diagonal), so that we are effectively left with an equal-order formulation with the stabilizing term [16] for further details). Although we shall stick to these three methods in our numerical examples, most popular pressure stabilization methods fit into the theory presented in Section 3.

Numerical examples
We now present a series of numerical examples in order to verify our error analysis. Relative pressure errors will be measured as When having (u, p) ∈ H 3 (Ω) × H 2 (Ω), Eichel et al. [12] proved a O(h 3/2 ) pressure convergence in L 2 (Ω) for first-order elements in uniform triangular meshes diagonally refined from a tensor-product mesh. Hence, we aim to verify here what may happen in a more general setup. For this, we consider first-order discretizations of problems having (u, p) ∈ H 2 (Ω) × H 2 (Ω) but u / ∈ H 3 (Ω), which fits our theory but not the existing one [12]. For the refinement studies we have Ω = (0, 1) 2 , starting from the mesh depicted in Figure 2 and then applying several levels of standard red uniform refinement. The linear algebraic system is solved directly to avoid the influence of iterative solver tolerances. The first example has an exact solution: , p = 4π cos 4πx sin 4πy , where the second-order gradient of the velocity field has a discontinuity at x = 1/2 due to the discontinuous right-hand side vector f = −16π 2 sin 4πx sin 4πy 8π 2 cos 4πx [2 cos 4πy − 1 + sign(1/2 − x)] . Table 1: Stokes flow with non-homogeneous boundary conditions: relative pressure error and estimated order of convergence (eoc) for different stabilization methods. The uniform refinement study starts from a coarse, non-orthogonal grid (cf. Figure 2 We solve the problem using three different stabilization methods, as described in Section 4. The results of the refinement study, displayed in Table 1 , confirm our a priori estimates: for each of the three methods, a slope between 1.5 and 2 holds for several levels, but after nine or ten levels of refinement the convergence starts slowing down towards a linear behaviour. For comparison, we present in Table 2 the velocity convergence in the H 1 (Ω) semi-norm, showing that the predicted first order is already reached around the sixth level of refinement.
Then, to illustrate the problem-dependent nature of the initial higher-order pressure convergence, we consider a problem with different solution. In the same unit square as before, with homogeneous Dirichlet boundary conditions and discontinuous forcing term we get the analytical solution again with u ∈ H 2 (Ω)\H 3 (Ω). We first perform a refinement study starting from the same non-orthogonal grid shown in Figure 2. The results, depicted in Table 3, highlight the dependence of the initial higher-order behaviour upon the chosen stabilization method: the PPE and PPP methods display higher orders for a few levels (breaking down sooner than in the previous example), while for the MINI elements the 1.5 slope is never reached.  Table 3: Stokes flow with homogeneous boundary conditions: relative pressure error and estimated order of convergence (eoc) for different stabilization methods. The uniform refinement study starts from a coarse, non-orthogonal grid (cf. Figure 2  Finally, we perform a refinement study with orthogonal grids: the initial mesh is similar to the one in Figure 2, but with the inner node now centralized at (x, y) = (1/2, 1/2). The results are very interesting: this time, the higher-order convergence does not break down within the present range of h, and in fact the MINI elements only reach the 1.5 slope at the finest levels. Whether the convergence would eventually become linear at finer levels cannot be predicted from our theory.

Conclusions
In this work, we have presented a numerical analysis on the pressure convergence of some classical finite element discretizations of the Stokes system. Although it is widely known that optimal convergence can be attained by going one degree higher in the velocity discretization, very little has been published to date on the pressure convergence of unbalanced pairs such as in MINI or equal-order elements. Thus, our main goal has been to answer a rather old question surrounding finite-element-based incompressible flow approximations: do we really lose one full order in the pressure convergence by not using balanced, Taylor-Hood-like pairs? By considering the pressure Schur complement formulation arising after eliminating the velocity, we have been able to show that, depending on certain constants, the pressure may in fact converge one or half an order faster than predicted by standard mixed finite element theory -but not necessarily for long. We can thus speak of a conditional and initial higher-order pressure convergence in unbalanced approximations. Taking stabilized first-order elements as a model discretization, our numerical examples confirm these results: a higher slope than one may indeed occur, but whether and for how long depends on the exact solution, the triangulation, the stabiliza-tion method, among other factors. In some cases the higher slope may not break down even after several levels of refinement, but our numerical counter-examples clearly confirm that this cannot be expected in general. Although we have considered the Stokes system as a model setting, we expect similar results to hold for Navier-Stokes flows. We hope and expect that our investigation can bring some clarity into the selection of finite element spaces for incompressible flow simulations, an important part of conceptual software design.