Quasi-optimal and pressure robust discretizations of the Stokes equations by new augmented Lagrangian formulations

TL;DR

This paper introduces new discretization methods for the stationary Stokes equations that are quasi-optimal and pressure robust, using augmented Lagrangian formulations to improve accuracy without requiring divergence-free finite element pairs.

Contribution

The paper develops a novel augmented Lagrangian approach that achieves pressure robustness and quasi-optimality for various finite element pairs in Stokes discretizations, without divergence-free constraints.

Findings

Velocity error proportional to best approximation error

Pressure error bounded by approximation errors

Applicable to conforming and nonconforming finite element pairs

Abstract

We approximate the solution of the stationary Stokes equations with various conforming and nonconforming inf-sup stable pairs of finite element spaces on simplicial meshes. Based on each pair, we design a discretization that is quasi-optimal and pressure robust, in the sense that the velocity $H^1$-error is proportional to the best $H^1$-error to the analytical velocity. This shows that such a property can be achieved without using conforming and divergence-free pairs. We bound also the pressure $L^2$-error, only in terms of the best approximation errors to the analytical velocity and the analytical pressure. Our construction can be summarized as follows. First, a linear operator acts on discrete velocity test functions, before the application of the load functional, and maps the discrete kernel into the analytical one. Second, in order to enforce consistency, we employ a new augmented Lagrangian formulation, inspired by Discontinuous Galerkin methods.

Figures

Tables

Table 1. Table 5.1. Test case § 5.2 . Velocity H 1 superscript 𝐻 1 H^{1} -errors of the new ℙ 2 / ℙ 0 subscript ℙ 2 subscript ℙ 0 \mathbb{P}_{2}/\mathbb{P}_{0} discretization and corresponding EOCs with composite (left) or standard (right) quadrature rules for α ∈ { 1 , 10 3 } 𝛼 1 superscript 10 3 \alpha\in\{1,10^{3}\} .

	$\alpha =1$	$\alpha ={10}^3$
N	$H^1$-error   EOC	$H^1$-error   EOC
4	3.32e-04	3.32e-04
5	8.31e-05    1.00	8.31e-05    1.00
6	2.08e-05    1.00	2.08e-05    1.00
7	5.19e-06    1.00	5.19e-06    1.00
8	1.30e-06    1.00	1.30e-06    1.00

Equations

$$-\mu \Delta u+\nabla p=f\text{and}\mathrm{div}u=0\text{in }\Omega ,u=0\text{on }\partial \Omega$$

$$\Vert \nabla (u-u_h){\Vert }_{L^2(\Omega )}\le c\underset{w_h\in V_h}{\inf }\Vert \nabla (u-w_h){\Vert }_{L^2(\Omega )}$$

$$\Vert \nabla (u-u_h){\Vert }_{L^2(\Omega )}\le c\left(\underset{w_h\in V_h}{\inf }\Vert \nabla (u-w_h){\Vert }_{L^2(\Omega )}+\frac{1}{\mu }\underset{q_h\in Q_h}{\inf }\Vert p-q_h{\Vert }_{L^2(\Omega )}\right).$$

$$\Vert u-u_h{\Vert }_h\le c\left(\underset{w_h\in V_h}{\inf }\Vert u-w_h{\Vert }_h+\frac{1}{\mu }\underset{q_h\in Q_h}{\inf }\Vert p-q_h{\Vert }_{L^2(\Omega )}+|||(u,p)||{|}_h\right)$$

$$\Vert u-u_h{\Vert }_h\le c\left(\underset{w_h\in V_h}{\inf }\Vert u-w_h{\Vert }_h+|||(u,0)||{|}_h\right)$$

$$\forall v\in H_0^1(\Omega {)}^d$$

$$\forall q\in L_0^2(\Omega )$$

$$\forall v_h\in V_h$$

$$\forall q_h\in Q_h$$

$$\forall q_h\in Q_h\beta \Vert q_h{\Vert }_{L^2(\Omega )}\le \underset{v_h\in V_h}{\sup }\frac{b_h(v_h,q_h)}{\Vert v_h{\Vert }_h}$$

$$\mu \Vert u-u_h{\Vert }_h+\Vert p-p_h{\Vert }_{L^2(\Omega )}\le C\left(\mu \underset{w_h\in V_h}{\inf }\Vert u-w_h{\Vert }_h+\underset{q_h\in Q_h}{\inf }\Vert p-q_h{\Vert }_{L^2(\Omega )}\right)$$

$$\begin{array}{c}{V}_h\subseteq H_0^1(\Omega {)}^d{E}_h={\mathrm{Id}}_{V_h}\\ a_h(w_h,v_h)={\int }_{\Omega }\nabla w_h:\nabla v_h{b}_h(v_h,q_h)=-{\int }_{\Omega }{q}_h\mathrm{div}{v}_h\end{array}$$

$$a_h({\phi }_i,{\phi }_j)b({\phi }_i,{\psi }_k)⟨f,E_h{\phi }_i⟩$$

$$Z:=\{z\in H_0^1(\Omega {)}^d\mid \mathrm{div}z=0\}.$$

$$\forall z\in Z\mu {\int }_{\Omega }\nabla u:\nabla z=⟨f,z⟩$$

$$\forall q_h\in Q_h{\int }_{\Omega }{q}_h{\underline{\mathrm{div}}}_{\text{h}}{v}_h=-b_h(v_h,q_h)$$

$$Z_h:=\{z_h\in V_h\mid {\underline{\mathrm{div}}}_{\text{h}}{z}_h=0\}.$$

$$\forall z_h\in Z_h\mu a_h(u_h,z_h)=⟨f,E_h{z}_h⟩$$

$$\underset{z_h\in Z_h}{\inf }\Vert u-z_h{\Vert }_h\le \left(1+{\beta }^{-1}\right)\underset{w_h\in V_h}{\inf }\Vert u-w_h{\Vert }_h$$

$$\Vert u-u_h{\Vert }_h\le C\underset{w_h\in V_h}{\inf }\Vert u-w_h{\Vert }_h$$

$$\begin{array}{c}{V}_h\subseteq H_0^1(\Omega {)}^d\mathrm{div}{V}_h=Q_h{E}_h={\mathrm{Id}}_{V_h}\\ a_h(w_h,v_h)={\int }_{\Omega }\nabla w_h:\nabla v_h{b}_h(v_h,q_h)=-{\int }_{\Omega }{q}_h\mathrm{div}{v}_h\end{array}$$

$$\forall v_h\in V_h,p\in Q_h$$

$$\forall u\in Z\cap Z_h,v_h\in V_h$$

$$E_h(Z_h)\subseteq Z$$

$$\forall u\in Z\cap Z_h,z_h\in Z_h{a}_h(u,z_h)={\int }_{\Omega }\nabla u:\nabla E_h{z}_h.$$

$$\forall v_h\in V_h\mathrm{div}{E}_h{v}_h={\underline{\mathrm{div}}}_{\text{h}}{v}_h.$$

$$S_{\ell }^k:=\{v\in H^k(\Omega )\mid \forall K\in \mathcal{M}{v}_{|K}\in {\mathbb{P}}_{\ell }(K)\}$$

$${\mathring{S}}_{\ell }^1:=S_{\ell }^1\cap H_0^1(\Omega )\text{and}{\hat{S}}_{\ell }^k:=S_{\ell }^k\cap L_0^2(\Omega ).$$

$$V_h=({\mathring{S}}_{\ell }^1{)}^2\text{and}{Q}_h={\hat{S}}_{\ell -2}^0,b_h(v_h,q_h)=-{\int }_{\Omega }{q}_h\mathrm{div}{v}_h$$

$${\underline{\mathrm{div}}}_{\text{h}}{v}_h={\Pi }_{\ell -2}^K\mathrm{div}{v}_h\text{in}K$$

$$\forall v_h\in ({\mathring{S}}_{\ell }^1{)}^2$$

$$\forall q_h\in {\hat{S}}_{\ell -2}^0$$

$${\mathring{S}}_{\ell }^1({\mathcal{M}}_K)\text{and}{\hat{S}}_{\ell -1}^0({\mathcal{M}}_K)$$

Full text

1. Introduction

We consider the discretization of the stationary Stokes equations

[TABLE]

with viscosity $\mu>0$, in a bounded domain $\Omega\subseteq\mathbb{R}^{d}$, $d\in\{2,3\}$. According to the classical approach of Brezzi [12], we approximate the analytical velocity $u$ and the analytical pressure $p$ by means of discrete spaces $V_{h}$ and $Q_{h}$, which are required to fulfill the so-called inf-sup condition. We additionally assume that $V_{h}$ and $Q_{h}$ are finite element spaces on a simplicial mesh of $\Omega$.

To motivate our work, let us focus on the velocity $H^{1}$-error, i.e. the error between $u$ and the discrete velocity $u_{h}$, measured in the $H^{1}$-norm. We refer to [8, Chapter 5] for the proof of the results listed hereafter. The Céa’s-type quasi-optimal estimate

[TABLE]

is well-known for standard discretizations (see (2.2) and (2.11) below) with conforming and divergence-free pairs, i.e. under the assumptions $V_{h}\subseteq H^{1}_{0}(\Omega)^{d}$ and $\operatorname{div}V_{h}=Q_{h}$. Such pairs have attracted a growing interest in recent years; see [18, 19, 32, 37] and the references therein. Owing to (1.2), this class of discretizations seems particularly attractive, because it fully exploits, up to a constant, the approximation properties of the space $V_{h}$ in the $H^{1}$-norm. This prevents, in particular, from the following issues.

For standard discretizations with general conforming pairs (see (2.2) and (2.5) below) one typically has

[TABLE]

Thus, if $\operatorname{div}V_{h}\neq Q_{h}$, the right-hand side suggests that the velocity $H^{1}$-error may be not robust with respect to the pressure. This is indeed the case and such effect is known in the literature as poor mass conservation. It becomes extreme for purely irrotational loads or for small values of the viscosity; see, for instance, [24]. Poor mass conservation discourages, in particular, from the use of unbalanced pairs, i.e. pairs $V_{h}/Q_{h}$ so that the approximation power of $V_{h}$ in the $H^{1}$-norm is higher than the one of $Q_{h}$ in the $L^{2}$-norm; cf. Remark 3.1.

Recall also that, in the nonconforming case $V_{h}\nsubseteq H^{1}_{0}(\Omega)^{d}$, estimates in the form

[TABLE]

are often derived. Here $\|\cdot\|_{h}$ is an extension of the $H^{1}$-norm to $H^{1}_{0}(\Omega)^{d}+V_{h}$ and the semi-norm $\lvert\!\lvert\!\lvert{\cdot}\rvert\!\rvert\!\rvert_{h}$ is defined on (a subspace of) $H^{1}_{0}(\Omega)^{d}\times L^{2}(\Omega)$. Since the lack of smoothness in $V_{h}$ is commonly compensated by additional regularity of the load beyond $H^{-1}(\Omega)^{d}$, the semi-norm $\lvert\!\lvert\!\lvert{\cdot}\rvert\!\rvert\!\rvert_{h}$ cannot be extended to $H^{1}_{0}(\Omega)^{d}\times L^{2}_{0}(\Omega)$ and potentially dominates the right-hand side of (1.4) for rough solutions. Therefore, an estimate like (1.3) cannot be expected to hold, cf. Remark 2.3.

Several techniques are available in the literature to deal with the above mentioned difficulties. The discretization of [4, section 6] and the general framework in [33] indicate how to avoid the issue with $\lvert\!\lvert\!\lvert{\cdot}\rvert\!\rvert\!\rvert_{h}$ for nonconforming pairs. The over-penalized augmented Lagrangian formulation of [10] and the grad-div stabilization [28] may serve to mitigate the impact of poor mass conservation. More recently, Linke et al. [23, 24, 25] proposed a class of discretizations, which differ from standard ones only in the treatment of the load and enjoy the following pressure robust upper bound

[TABLE]

for several conforming and nonconforming pairs.

In this paper, we show that the quasi-optimal and pressure robust estimate (1.2) is not a prerogative of conforming and divergence-free pairs, but can be achieved also by (carefully designed) discretizations, based on general inf-sup stable pairs. In this way, we combine the advantages of the various techniques listed above. We also bound the pressure $L^{2}$-error only in terms of the best approximation errors to the analytical velocity and to the analytical pressure. To our best knowledge, similar error bounds were previously obtained only in [35] in the rather specific case of the lowest-order nonconforming Crouzeix-Raviart pair [14]. In particular, our results make unbalanced pairs a valuable option, if one is more interested in the analytical velocity rather than in the analytical pressure.

Our approach is guided by few simple necessary conditions and builds on two main ingredients. First, we discretise the load with the help of an operator which maps $V_{h}$ into $H^{1}_{0}(\Omega)^{d}$ and discretely divergence-free into exactly divergence-free functions. The importance of the latter property was first devised in [24]. For this purpose, we solve local Stokes problems with Scott-Vogelius elements on a barycentric refinement of the mesh, see [19, 30, 36]. Second, we discretise the weak form of the Laplace operator in a way inspired by Discontinuous Galerkin (DG) methods, in order to enforce the necessary consistency. The resulting discretization can be interpreted as a new augmented Lagrangian formulation, cf. Remark 3.7.

The rest of the paper is organized as follows. In section 2 we set up the abstract framework. In section 3 we illustrate our construction by means of a model example. Various generalizations are then discussed in section 4. Finally, in section 5 we complement our theoretical findings through some numerical experiments.

2. Abstract framework

This section introduces an abstract discretization of (1.1) and the properties in which we are interested. Two basic results are also proved. We use standard notations for Lebesgue and Sobolev spaces.

2.1. Quasi-optimal discretizations

Let $\Omega\subseteq\mathbb{R}^{d}$, $d\in\{2,3\}$, be an open and bounded polytopic domain with Lipschitz-continuous boundary. The weak formulation of the stationary Stokes equations in $\Omega$, with viscosity $\mu>0$ and load $f\in H^{-1}(\Omega)^{d}$, looks for $u\in H^{1}_{0}(\Omega)^{d}$ and $p\in L^{2}_{0}(\Omega)$ such that

[TABLE]

Here $\colon$ denotes the euclidean scalar product of $d\times d$ tensors and $\left\langle\cdot,\cdot\right\rangle$ is the dual pairing of $H^{-1}(\Omega)^{d}$ and $H^{1}_{0}(\Omega)^{d}$. Due to the boundary condition on the analytical velocity $u$, the analytical pressure $p$ belongs to $L^{2}_{0}(\Omega):=\{q\in L^{2}(\Omega)\mid\int_{\Omega}q=0\}$. Problem (2.1) is uniquely solvable, according to [8, Theorem 8.2.1].

Remark 2.1* (Alternative formulation).*

Most of our subsequent results remain unchanged in case the gradient is replaced by the symmetric gradient in the first equation of (2.1) and the homogeneous Neumann condition is imposed on (a portion of) $\partial\Omega$. The only remarkable difference is that a piecewise Korn’s inequality may fail to hold for some of the nonconforming pairs mentioned in section 4.1, see [2, 11]. This problem, however, can be overcome e.g. by an additional jump penalization in the spirit of [34, Section 3.3].

We consider discretizations that mimic the variational structure of problem (2.1). More precisely, we approximate $u$ and $p$ in finite-dimensional linear spaces $V_{h}$ and $Q_{h}$. We require $Q_{h}\subseteq L^{2}_{0}(\Omega)$ and measure the pressure error in the $L^{2}$-norm $\|\cdot\|_{L^{2}(\Omega)}$. Instead, we allow for nonconforming discrete velocity spaces $V_{h}\nsubseteq H^{1}_{0}(\Omega)^{d}$. In order to measure the velocity error, we assume that an extension $\|\cdot\|_{h}$ of the $H^{1}$-norm $\|\operatorname{\nabla}\cdot\|_{L^{2}(\Omega)}$ to $H^{1}_{0}(\Omega)^{d}+V_{h}$ is at our disposal. We replace the bilinear forms in (2.1) with discrete surrogates $a_{h}:V_{h}\times V_{h}\to\mathbb{R}$ and $b_{h}:V_{h}\times Q_{h}\to\mathbb{R}$. Moreover, we let $E_{h}:V_{h}\to H^{1}_{0}(\Omega)^{d}$ be a linear operator. Hence, we look for a discrete velocity $u_{h}\in V_{h}$ and a discrete pressure $p_{h}\in Q_{h}$ such that

[TABLE]

To ensure that this problem is uniquely solvable, we assume hereafter that $a_{h}$ is coercive on $V_{h}$ and that the pair $V_{h}/Q_{h}$ is inf-sup stable, i.e.

[TABLE]

for some constant $\beta>0$, see [8, Corollary 4.2.1]. Note, in particular, that the duality $\left\langle f,E_{h}v_{h}\right\rangle$ is well-defined for all $f\in H^{-1}(\Omega)^{d}$ and $v_{h}\in V_{h}$, also in the nonconforming case.

We shall pay special attention to the following property, which guarantees that $(u_{h},p_{h})$ is a near-best approximation of $(u,p)$ in $V_{h}\times Q_{h}$.

Definition 2.2** (Quasi-optimality).**

Denote by $(u,p)$ and $(u_{h},p_{h})$ the solutions of (2.1) and (2.2), respectively, with load $f$ and viscosity $\mu$. We say that (2.2) is a quasi-optimal discretization of (2.1) when there is a constant $C\geq 1$ such that

[TABLE]

for all $f\in H^{-1}(\Omega)^{d}$ and $\mu>0$. We denote by $C_{\mathrm{qo}}$ the smallest such constant.

According to [8, Theorem 5.2.5], the discretization (2.2) is quasi-optimal if

[TABLE]

i.e. if $V_{h}/Q_{h}$ is a conforming pair and $a_{h}$, $b_{h}$ and $E_{h}$ are simple restrictions of their conforming counterparts in (2.1). In sections 3 and 4 we show that quasi-optimality can be achieved also with nonconforming pairs and/or for different choices of $a_{h}$ and $E_{h}$.

Remark 2.3* (Smoothing by $E_{h}$).*

Since $V_{h}$ is finite-dimensional, the operator $E_{h}$ is bounded and the solution of (2.2) depends continuously on the $H^{-1}$-norm of $f$. This property, in turn, prevents the issue pointed out in the introduction concerning the semi-norm $\lvert\!\lvert\!\lvert{\cdot}\rvert\!\rvert\!\rvert_{h}$ in (1.4). Of course, such observation is of practical interest only if the norm of $E_{h}$ is of moderate size, so that it does not affect too much the stability constant of (2.2). We call $E_{h}$ ”smoothing” operator, because it increases the smoothness of the elements of $V_{h}$ whenever $V_{h}\nsubseteq H^{1}_{0}(\Omega)^{d}$. For conforming pairs, one can let $E_{h}$ be the identity as in (2.5). This choice is compatible with quasi-optimality but, possibly, it is not pressure robust; compare with section 2.2 below.

Remark 2.4* (Computational feasibility).*

It is highly desirable that there are bases $\{\varphi_{1},\dots,\varphi_{N}\}$ and $\{\psi_{1},\dots,\psi_{M}\}$ of $V_{h}$ and $Q_{h}$, respectively, such that the scalars

[TABLE]

can be computed or approximated, up to a prescribed tolerance, with $O(1)$ operations, for all $i,j=1,\dots,N$ and $k=1,\dots,M$. This ”computational feasibility” is not necessary for quasi-optimality but guarantees that the solution of (2.2) can be computed with optimal complexity.

2.2. Quasi-optimal and pressure robust discretizations

The analytical velocity $u$ solving (2.1) can be equivalently characterized as the solution of an elliptic problem. In fact, the second equation imposes that $u$ is divergence-free or, in other words, that it is an element of the kernel

[TABLE]

Then, testing the first equation with an arbitrary element of $Z$, we obtain the reduced problem

[TABLE]

which is uniquely solvable, according to the Lax-Milgram lemma and the Friedrichs inequality.

The same structure can be observed at the discrete level. To see this, we first introduce the discrete divergence $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}:V_{h}\to Q_{h}$ by

[TABLE]

for all $v_{h}\in V_{h}$. The second equation of (2.2) imposes that $u_{h}$ is discretely divergence-free, i.e. it is an element of the discrete kernel

[TABLE]

Then, testing the first equation with an arbitrary element of $Z_{h}$, we derive the discrete reduced problem

[TABLE]

which is uniquely solvable, since $a_{h}$ is coercive on $V_{h}$. In the vein of [12, Remark 2.1], it is worth recalling that this is a (possibly) nonconforming discretization of (2.6), because $Z_{h}$ may fail to be a subspace of $Z$, even if $V_{h}\subseteq H^{1}_{0}(\Omega)^{d}$.

Similarly as in Definition 2.2, we will be interested in the question whether $u_{h}$ is a near-best approximation of $u$ in $Z_{h}$. This actually amounts to ask whether $u_{h}$ is near-best in $V_{h}$, because the inf-sup condition (2.3) implies

[TABLE]

according to [8, Proposition 5.1.3] and [29, Lemma 2.1].

Definition 2.5** (Quasi-optimality and pressure robustness).**

Denote by $u$ and $u_{h}$ the solutions of (2.6) and (2.8), respectively, with load $f$ and viscosity $\mu$. We say that (2.2) is a quasi-optimal and pressure robust discretization of (2.1) when there is a constant $C\geq 1$ such that

[TABLE]

for all $f\in H^{-1}(\Omega)^{d}$ and $\mu>0$. We denote by $\mathrm{C_{\mathrm{qopr}}}$ the smallest such constant.

Problem (2.6) reveals that the analytical velocity $u$ is independent of the pressure $p$ and depends on the load $f$ only through its restriction to $Z$. This implies, for instance, that $u$ is invariant with respect to irrotational perturbations of $f$, see Linke [24]. The near-best estimate (2.10) guarantees that $u_{h}$ reproduces such invariance property at the discrete level and justifies the designation ”pressure robust”.

The discretization (2.2) is known to be quasi-optimal and pressure robust if

[TABLE]

i.e. if $V_{h}/Q_{h}$ is a conforming and divergence-free pair and $a_{h}$, $b_{h}$ and $E_{h}$ are simple restrictions of their continuous counterparts in (2.1). In fact, in this case, we have $Z_{h}\subseteq Z$ and (2.8) is a conforming Galerkin discretization of (2.6). Therefore, Céa’s lemma and (2.9) imply $\mathrm{C_{\mathrm{qopr}}}\leq(1+\beta^{-1})$. It is our purpose to show that quasi-optimality and pressure robustness can be achieved also by other discretizations than (2.11).

2.3. Necessary consistency conditions

The left- and the right-hand sides of (2.4) are seminorms on $Z\times L^{2}_{0}(\Omega)$ and the kernel of the latter is $(Z\cap Z_{h})\times Q_{h}$, as a consequence of (2.9). Quasi-optimality actually prescribes that such seminorms are equivalent, because the converse of (2.4) immediately follows from the inclusion $(u_{h},p_{h})\in Z_{h}\times Q_{h}$. Hence, a simple necessary condition is that the kernels of the two seminorms coincide. In other words, whenever the solution $(u,p)$ of (2.1) is in $Z_{h}\times Q_{h}$, it must solve also (2.2). This is an algebraic consistency condition, which can be rephrased in terms of the forms $a_{h}$ and $b_{h}$ and of the operator $E_{h}$, in the spirit of [33, Definition 2.7].

Lemma 2.6** (Consistency for quasi-optimality).**

Assume that (2.2) is a quasi-optimal discretization of (2.1). Then, necessarily we have

[TABLE]

Proof.

Denote by $(u,p)$ the solution of (2.1) and assume first $u=0$ and $p\in Q_{h}$. Quasi-optimality implies that the solution $(u_{h},p_{h})$ of (2.2) satisfies $u_{h}=0$ and $p_{h}=p$. Comparing the first equations of (2.1) and (2.2), we derive the identity $b_{h}(v_{h},p)=-\int_{\Omega}p\operatorname{div}E_{h}v_{h}$ for all $v_{h}\in V_{h}$. Condition (2.12a) then follows from the definition of $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ in (2.7). Next, assume $u\in Z\cap Z_{h}$ and $p=0$. Since quasi-optimality implies $u_{h}=u$ and $p_{h}=0$, condition (2.12b) can be derived comparing the first equations of (2.1) and (2.2) as before. ∎

The conforming discretization (2.5) is a simple option to fulfill (2.12), but not the only possible. Examples with nonconforming discrete velocity space can be found in [4, Section 6] and [35]. Standard nonconforming discretizations, like the one of Crouzeix and Raviart [14], do not fulfill (2.12), because they do not employ a smoothing operator. It is also worth noticing that (2.12) involves the interplay of $a_{h}$ and $b_{h}$ with $E_{h}$. This indicates that the discretization of the differential operator in (1.1) and the one of the corresponding load should not be regarded as independent tasks.

Proceeding similarly as in Lemma 2.6, we derive necessary conditions for quasi-optimality and pressure robustness.

Lemma 2.7** (Consistency for quasi-optimality and pressure robustness).**

Assume that (2.2) is a quasi-optimal and pressure robust discretization of (2.1). Then, necessarily we have

[TABLE]

Proof.

Let $z_{h}\in Z_{h}$ be such that $\operatorname{div}E_{h}z_{h}\neq 0$. Assuming that $(u,p)=(0,\operatorname{div}E_{h}z_{h})$ solves (2.1), we infer $\left\langle f,E_{h}z_{h}\right\rangle=-\|\operatorname{div}E_{h}z_{h}\|_{L^{2}(\Omega)}^{2}\neq 0$. Inserting this information in (2.8), we obtain $u_{h}\neq 0$. Therefore, we have $\|u-u_{h}\|_{h}>\inf_{v_{h}\in V_{h}}\|u-v_{h}\|_{h}=0$, which contradicts quasi-optimality and pressure robustness. This proves (2.13a). Assertion (2.13b) may be checked similarly to (2.12b) in Lemma 2.6. ∎

Condition (2.13b) is clearly necessary for (2.12b), while (2.13a) is neither necessary nor sufficient for (2.12a). We mention also that (2.13a) differs from the condition exploited in [25] to achieve pressure robustness, in that here $E_{h}$ is required to map into $H^{1}_{0}(\Omega)^{d}$ and not only into $H_{\mathrm{div}}(\Omega)$, cf. Remark 2.3.

Remark 2.8* (Failure of $E_{h}=\mathrm{Id}_{V_{h}}$).*

If $V_{h}/Q_{h}$ is a conforming and divergence-free pair, the abstract discretization (2.2) with (2.11) verifies the first necessary condition in Lemma 2.7. If, instead, the pair is conforming but not divergence-free, we have $Z_{h}\nsubseteq Z$. In this case, the operator $E_{h}$ cannot coincide with the identity on $Z_{h}$.

In the next sections, we design some new discretizations proceeding as follows. Given an inf-sup stable pair $V_{h}/Q_{h}$, together with the corresponding bilinear form $b_{h}$, we construct $a_{h}$ and $E_{h}$ so that the necessary conditions in Lemmas 2.6 and 2.7 hold true. Then, we use standard techniques from the analysis of saddle point problems to verify (2.4) and (2.10) and to bound the constants $C_{\mathrm{qo}}$ and $\mathrm{C_{\mathrm{qopr}}}$. Alternatively, one could exploit [33, Theorem 4.14], which guarantees that (2.13) is a sufficient condition for quasi-optimality and pressure robustness. Such result provides also a formula for $\mathrm{C_{\mathrm{qopr}}}$. Analogously, generalizing the framework of [33], one could show also that (2.12) is a sufficient condition for quasi-optimality and derive a formula for $C_{\mathrm{qo}}$. We prefer to proceed as indicated, to make sure this paper can be read independently of [33].

3. A paradigmatic discretization

Assume that we are given an inf-sup stable pair $V_{h}/Q_{h}$, together with the corresponding bilinear form $b_{h}$. A possible strategy to fulfill the necessary conditions (2.12a) and (2.13a) is to employ a ”divergence-preserving” smoothing operator, i.e.

[TABLE]

Once such operator is given, conditions (2.12b) and (2.13b) prescribe the restriction of $a_{h}$ on $(Z\cap Z_{h})\times V_{h}$. Then, inspired by [1] and [34], we extend the resulting form to $V_{h}\times V_{h}$, in a way that additionally ensures symmetry and coercivity. In order to keep the exposition as clear as possible, we first exemplify this idea in a model setting. We postpone various generalizations to the next section.

3.1. The unbalanced $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair

We consider hereafter pairs of finite element spaces on a face-to-face simplicial mesh $\mathcal{M}$ of $\Omega$ in the sense of [15, Definition 1.36]. We write $c$ for a nondecreasing and nonnegative function of the shape parameter of $\mathcal{M}$, which possibly depends also on different quantities (like, e.g., the space dimension), but neither on other properties of $\mathcal{M}$ nor on the viscosity $\mu$. Such constant may change at different occurrences. We occasionally abbreviate $a\leq cb$ as $a\lesssim b$ and $c^{-1}b\leq a\leq cb$ as $a\eqsim b$.

For all integers $\ell\geq 0$, we denote by $\mathbb{P}_{\ell}(S)$ the space of polynomials with total degree $\leq\ell$ on a simplex $S\subseteq\mathbb{R}^{d}$. The space of $H^{k}$-conforming element-wise polynomials on $\mathcal{M}$ then reads

[TABLE]

with $k\in\{0,1\}$ and the convention $H^{0}(\Omega):=L^{2}(\Omega)$. Motivated by the homogeneous boundary condition in (1.1), we consider the subspaces

[TABLE]

To exemplify our construction, we assume $d=2$ for the remaining part of this section. We consider the conforming $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair, which is given by

[TABLE]

with $\ell\geq 2$. The inf-sup condition (2.3) holds with $\beta^{-1}\leq c$, see [8, Remark 8.6.2].

Remark 3.1* (Unbalanced pairs).*

The $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair is unbalanced, in the sense that the approximation power $\ell-1$ of the discrete pressure space in the $L^{2}$-norm is strictly less than the approximation power $\ell$ of the discrete velocity space in the $H^{1}$-norm. Other examples can be obtained enriching the velocity space of any inf-sup stable pair. The use of conforming unbalanced pairs, in combination with the standard discretization (2.5), is discouraged by the error estimate (1.3) and Remark 2.8; see also [8, Remark 8.6.2]. Still, quasi-optimal and pressure robust discretizations based on such pairs would be a valuable option, if one is more interested in the analytical velocity rather than in the analytical pressure.

The discrete divergence $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ in the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair coincides with the $L^{2}$-orthogonal projection of the analytical divergence onto $\widehat{S}_{\ell-2}^{0}$. Since (2.7) actually holds for all discrete pressures in $S_{\ell-2}^{0}$, we can compute $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ element-wise as follows

[TABLE]

for all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $K\in\mathcal{M}$, where $\Pi^{K}_{\ell-2}$ is the $L^{2}$-orthogonal projection onto $\mathbb{P}_{\ell-2}(K)$. Therefore, denoting by $Z_{h}^{ub}$ the discrete kernel, we conclude $Z_{h}^{{ub}}\nsubseteq Z$.111The superscript $"{ub}"$ stands for ”unbalanced”. Along this section, we use it to label spaces, forms and operators related to the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair. This confirms that the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair is conforming but not divergence-free.

The abstract discretization (2.2) with (2.5), based on the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair, states $u_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $p_{h}\in\widehat{S}_{\ell-2}^{0}$ such that

[TABLE]

3.2. Local inversion of the divergence

Proceeding as in [35], we enforce (3.1) with the help of local right inverses of the divergence. Such operators can be defined through discrete Stokes-like problems on the barycentric refinement of each element. To see this, fix $K\in\mathcal{M}$ and let $\mathcal{M}_{K}$ denote the triangulation of $K$ obtained connecting each vertex with the barycenter; cf. Figure 3.1. For $\ell\in\mathbb{N}$, we define the local spaces

[TABLE]

on $\mathcal{M}_{K}$ similarly to the global spaces $\mathring{S}_{\ell}^{1}$ and $\widehat{S}_{\ell-1}^{0}$ in (3.3). In particular, all $v_{k}\in\mathring{S}_{\ell}^{1}(\mathcal{M}_{K})$ vanish on $\partial K$ and all $q_{K}\in\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{K})$ are such that $\int_{K}q_{K}=0$. The pair $\mathring{S}_{\ell}^{1}(\mathcal{M}_{K})^{2}/\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{K})$ is conforming and divergence-free in $K$.

According to [19, Theorem 3.1], we have the local inf-sup stability

[TABLE]

This entails that we can define a linear operator $R_{\ell}^{K}:L^{2}(\Omega)\to H^{1}_{0}(\Omega)^{2}$ as follows. Given $q\in L^{2}(\Omega)$, let $u_{K}=u_{K}(q)\in\mathring{S}_{\ell}^{1}(\mathcal{M}_{K})^{2}$ and $p_{K}=p_{K}(q)\in\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{K})$ solve

[TABLE]

Hence, we set

[TABLE]

Proposition 3.2** (Local right inverses).**

Let $K\in\mathcal{M}$ be a mesh element and $\ell\in\mathbb{N}$. The operator $R_{\ell}^{K}$ is well-defined and, for all $q\in L^{2}(\Omega)$, we have

[TABLE]

Proof.

The operator $R_{\ell}^{K}$ is well-defined and satisfies (3.9a) in view of the local inf-sup (3.7) and [8, Corollary 4.2.1]. The property in (3.9b) directly follows from the second equation of problem (3.8), because $\operatorname{div}u_{K}\in\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{K})$. ∎

Remark 3.3* (Computation of the local right inverses).*

In what follows, we shall need to compute $R_{\ell}^{K}q$ for all $K\in\mathcal{M}$ and various $q\in S_{\ell-1}^{0}$. To this end, a possible strategy is to precompute the solution of (3.8) on a reference triangle $K_{\mathrm{ref}}$, for all possible loads $q_{\mathrm{ref}}$ in a basis of $\mathbb{P}_{\ell-1}(K_{\mathrm{ref}})$. The computational complexity of this task only depends on $\ell$. Then, the solution of (3.8) in $K$ can be obtained in terms of the corresponding solution in $K_{\mathrm{ref}}$, by means of the contravariant Piola transformation; see [8, Section 2.1.3].

We have considered here the two-dimensional case only to be consistent with the simplification introduced in section 3.1. The same construction is actually possible in any space dimension $d\geq 2$.

3.3. A new augmented Lagrangian formulation

We now propose a new discretization of the Stokes equations, based on the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair. The first ingredient of our construction is a linear operator $E_{h}^{ub}:(\mathring{S}_{\ell}^{1})^{2}\to H^{1}_{0}(\Omega)^{2}$ fulfilling (3.1). In view of $Z_{h}^{{ub}}\nsubseteq Z$ and Remark 2.8, the identity on $(\mathring{S}_{\ell}^{1})^{2}$ cannot accommodate this property. Therefore, we introduce a ”divergence correction” $R_{h}^{ub}:(\mathring{S}_{\ell}^{1})^{2}\to H^{1}_{0}(\Omega)^{2}$

[TABLE]

Proposition 3.4** (Divergence-preserving smoothing operator).**

The linear operator $E_{h}^{ub}:(\mathring{S}_{\ell}^{1})^{2}\to H^{1}_{0}(\Omega)^{2}$ given by

[TABLE]

fulfills (3.1) and is such that, for all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$,

[TABLE]

Proof.

For all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $K\in\mathcal{M}$, it holds

[TABLE]

In view of (3.5), we have $\int_{K}(\operatorname{\underline{\operatorname{div}}_{\textit{h}}}v_{h}-\operatorname{div}v_{h})=0$. Since the inclusion $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ implies also $(\operatorname{\underline{\operatorname{div}}_{\textit{h}}}v_{h}-\operatorname{div}v_{h})_{|K}\in\mathbb{P}_{\ell-1}(K)$, Proposition 3.2 and the identity above ensure that $E_{h}^{ub}$ fulfills (3.1). This, in turn, easily implies the lower bound ”$\gtrsim$” in (3.11). The corresponding upper bound ”$\lesssim$” is a consequence of the identity $\|\operatorname{\nabla}(v_{h}-E_{h}^{ub}v_{h})\|_{L^{2}(K)}=\|\operatorname{\nabla}R^{K}_{\ell}v_{h}\|_{L^{2}(K)}$, $K\in\mathcal{M}$, combined with (3.9a). ∎

The second ingredient of our construction is a suitable bilinear form $a_{h}$. Accounting for the definition of $E_{h}^{ub}$ in (3.10), the necessary conditions (2.12b) and (2.13b) prescribe

[TABLE]

for all $u\in Z\cap Z_{h}^{{ub}}$ and $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$. A simple option would be to let the right-hand side define $a_{h}$ on $(\mathring{S}_{\ell}^{1})^{2}\times(\mathring{S}_{\ell}^{1})^{2}$. Still, it has to be noticed that the second summand $\int_{\Omega}\operatorname{\nabla}u\colon\operatorname{\nabla}R_{h}^{ub}v_{h}=-\sum_{K\in\mathcal{M}}\int_{K}\operatorname{\Delta}u\cdot R_{h}^{ub}v_{h}$ cannot be expected to vanish. Therefore, it obstructs the symmetry and, possibly, also the nondegeneracy of $a_{h}$. To overcome this problem, we observe that $R_{h}^{ub}$ vanishes on $Z\cap Z_{h}^{{ub}}$, according to (3.11). This suggests to re-establish symmetry and nondegeneracy mimicking the construction of the Symmetric Interior Penalty (DG-SIP) discretization of second-order problems, see [1] or [15, section 4.2.1]. Thus, we set $a_{h}=a_{h}^{ub}$, where

[TABLE]

where $\eta>0$ is a penalty parameter. Note that $a_{h}^{ub}$ fulfills (3.12).

The abstract discretization (2.2) with the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair, $a_{h}=a_{h}^{ub}$ and $E_{h}=E_{h}^{ub}$ reads as follows: Find $u_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $p_{h}\in\widehat{S}_{\ell-2}^{0}$ such that

[TABLE]

We begin our discussion on the new discretization by checking that a solution $(u_{h},p_{h})$ exists and is unique. In view of the above-mentioned inf-sup stability of the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair, it suffices to prove that $a_{h}^{ub}$ is coercive on $(\mathring{S}_{\ell}^{1})^{2}$. We proceed similarly as in [15, Lemma 4.1.2].

Lemma 3.5** (Coercivity of $a_{h}^{ub}$).**

The bilinear form $a_{h}^{ub}$ is coercive on $(\mathring{S}_{\ell}^{1})^{2}$ for all $\eta>1$ and we have

[TABLE]

for all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$.

Proof.

Let $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$. Setting $w_{h}=v_{h}$ in (3.13), we obtain

[TABLE]

The Cauchy-Schwartz and the weighted Young’s inequality further provide the upper bound $2\left\lvert\int_{\Omega}\operatorname{\nabla}v_{h}\colon\operatorname{\nabla}R_{h}^{ub}v_{h}\right\rvert\leq\eta^{-1}\|\operatorname{\nabla}v_{h}\|_{L^{2}(\Omega)}^{2}+\eta\|\operatorname{\nabla}R_{h}^{ub}v_{h}\|_{L^{2}(\Omega)}^{2}$. Inserting this inequality into the previous identity concludes the proof. ∎

Let us comment on the cost for assembling and solving the new discretization.

Remark 3.6* (Feasibility of the new discretization).*

Assume that $\{\varphi_{1},\dots,\varphi_{N}\}$ and $\{\psi_{1},\dots,\psi_{M}\}$ are nodal bases of $(\mathring{S}_{\ell}^{1})^{2}$ and $\widehat{S}_{\ell-2}^{0}$, respectively. All functions $\varphi_{i}$ and $\psi_{k}$, with $i=1,\dots,N$ and $k=1,\dots,M$, are locally supported. Hence, the construction of $E_{h}^{ub}\varphi_{i}$ involves the solution of a limited number of local problems (3.8) and we have $\mathrm{supp}(E_{h}^{ub}\varphi_{i})\subseteq\mathrm{supp}(\varphi_{i})$. Moreover, thanks to the local characterization of the discrete divergence (3.5), the entire computation of $E_{h}^{ub}\varphi_{i}$ requires $O(1)$ operations. This entails that the bilinear forms $a_{h}^{ub}(\varphi_{i},\varphi_{j})$ and $\int_{\Omega}\psi_{k}\operatorname{div}\varphi_{i}$ and the linear form $\langle f,E_{h}^{ub}\varphi_{i}\rangle$ can be evaluated with $O(1)$ operations for all $i,j=1,\dots,N$ and $k=1,\dots,M$. Thus, the discretization (3.14) is computationally feasible, in the sense of Remark 2.4. Let us mention also that the stiffness matrices associated with $a_{h}^{ub}$ and its counterpart in (3.6) are of course different but, for all $\eta>1$, their condition numbers differ, at most, by the ratio of the continuity and the coercivity constants of $a_{h}^{ub}$. This ratio is bounded by $c\eta^{2}(\eta-1)^{-1}$, as a consequence of Proposition 3.4 and Lemma 3.5.

The following remarks connect (3.14) with other existing discretizations.

Remark 3.7* (Connection with augmented Lagrangian formulations).*

In view of (3.11), the last summand $\eta\int_{\Omega}\operatorname{\nabla}R_{h}^{ub}w_{h}\colon\operatorname{\nabla}R_{h}^{ub}v_{h}$ in the definition of $a_{h}^{ub}$ penalizes the functions that are in the discrete kernel $Z_{h}^{{ub}}$ and not in $Z$. More precisely, the penalization is equivalent to $\eta\int_{\Omega}\operatorname{div}w_{h}\operatorname{div}v_{h}$ on $Z_{h}^{{ub}}$. This indicates that (3.14) can be interpreted as a new augmented Lagrangian formulation for the Stokes problem; see [8, Section 6.1]. The additional terms enforcing consistency and symmetry distinguish our formulation from previous ones.

Remark 3.8* (Connection with DG discretizations).*

The DG-SIP bilinear form in [1] consists of four terms. The first two terms serve to accommodate consistency, see [15, Section 4.2] or [34]. In particular, the second one arises due to the use of possibly nonconforming, i.e. discontinuous, functions. The two remaining terms are designed to further enforce symmetry and coercivity, respectively, still preserving consistency. The same structure can be observed in the form $a_{h}^{ub}$. Here nonconformity has to be intended in the sense that $Z_{h}^{ub}\nsubseteq Z$, i.e. discretely divergence-free functions are possibly not divergence-free. A remarkable difference from the DG-SIP bilinear form is that the coercivity of $a_{h}^{ub}$ can be guaranteed for all $\eta>1$ and not only for sufficiently large $\eta$.

Remark 3.9* (Connection with R-FEM discretizations).*

Rearranging terms in (3.13), we see that the form $a_{h}^{ub}$ can be rewritten as follows

[TABLE]

This sheds additional light on the condition $\eta>1$ in Lemma 3.5 and provides an interesting connection with the Recovered Finite Element Method (R-FEM) of Georgoulis and Pryer [17].

3.4. Error estimates

We now aim at showing that, unlike (3.6), (3.14) is a quasi-optimal and pressure robust discretization of (2.1). As a preliminary step, we bound the consistency error generated by the last two terms in the definition of $a_{h}^{ub}$. Such terms can be expected to generate a consistency error, as they were artificially added to the right-hand side of (3.12).

Lemma 3.10** (Consistency error).**

Let $\eta>1$ be given. We have

[TABLE]

for all $z_{h}\in Z_{h}^{{ub}}$ and $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$.

Proof.

The definitions of $a_{h}^{ub}$ and $E_{h}^{ub}$ imply

[TABLE]

The equivalence (3.11) reveals, in particular, $\|\operatorname{\nabla}R_{h}^{ub}z_{h}\|_{L^{2}(\Omega)}\lesssim\|\operatorname{\nabla}(z-z_{h})\|_{L^{2}(\Omega)}$ for all $z\in Z$. The characterization (3.5) of the discrete divergence $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ and (3.11) entail also $\|\operatorname{\nabla}(v_{h}+\eta R_{h}^{ub}v_{h})\|_{L^{2}(\Omega)}\lesssim\eta\|\operatorname{\nabla}v_{h}\|_{L^{2}(\Omega)}$. Inserting these bounds into the identity above concludes the proof. ∎

Recall from section 2.2 that the discrete velocity $u_{h}$ solving (3.14) is in the discrete kernel $Z_{h}^{{ub}}$ and can be equivalently characterized through the reduced problem

[TABLE]

Theorem 3.11** (Quasi-optimality and pressure robustness).**

For all $\eta>1$, problem (3.14) is a quasi-optimal and pressure robust discretization of (2.1) with constant $\mathrm{C_{\mathrm{qopr}}}\leq c\eta^{2}(\eta-1)^{-1}$.

Proof.

Denote by $u\in Z$ and $u_{h}\in Z_{h}^{{ub}}$ the solutions of problems (2.6) and (3.17), respectively, with load $f\in H^{-1}(\Omega)^{2}$ and viscosity $\mu>0$. Let $z_{h}\in Z_{h}^{{ub}}$ be arbitrary and define $v_{h}:=u_{h}-z_{h}$. Lemma 3.5 and problem (3.17) reveal

[TABLE]

Since $v_{h}\in Z_{h}^{{ub}}$, we have $E_{h}^{ub}v_{h}\in Z$ as a consequence of Proposition 3.4. Hence, problem (2.6) yields $\mu^{-1}\langle f,E_{h}^{ub}v_{h}\rangle=\int_{\Omega}\operatorname{\nabla}u\colon\operatorname{\nabla}E_{h}^{ub}v_{h}$. We insert this identity into the previous inequality and invoke Proposition 3.4 and Lemma 3.10. Owing to the inclusion $u\in Z$, it results

[TABLE]

We conclude taking the infimum over all $z_{h}\in Z_{h}$ and recalling (2.9). ∎

Let us mention that a better bound of the constant $\mathrm{C_{\mathrm{qopr}}}$ in terms of $\eta$, namely $\mathrm{C_{\mathrm{qopr}}}\leq c\eta(\eta-1)^{-1/2}$, could be obtained with the help of [33, Theorem 4.14]. Both, this estimate and the one in Theorem 3.11, suggest to set $\eta=2$. The next remark additionally confirm that we may have $\mathrm{C_{\mathrm{qopr}}}\to+\infty$ as $\eta\to+\infty$, thus pointing out the importance of explicitly knowing a safe value of the penalty parameter.

Remark 3.12* (Locking effect).*

The penalization in $a_{h}^{ub}$ imposes that the solution $u_{h}^{ub}$ of (3.17) approaches the subspace $Z\cap Z_{h}^{ub}$ for $\eta\to+\infty$, as a consequence of Proposition 3.4. This entails that the constant $\mathrm{C_{\mathrm{qopr}}}$ in Theorem 3.11 remains bounded in the limit $\eta\to+\infty$ only if the equivalence

[TABLE]

holds for all $z\in Z$. Conversely, if (3.18) holds, we can assume that the function $z_{h}$ in the proof of Theorem 3.11 varies only in $Z\cap Z_{h}^{ub}$. This, in turn, provides a robust upper bound of $\mathrm{C_{\mathrm{qopr}}}$ in the limit $\eta\to+\infty$. Whenever condition (3.18) fails, a locking effect may occur, in the sense of [3]. We illustrate this in section 5.3 by means of a numerical experiment.

Theorem 3.11 states that the discretization (3.14) enjoys a better velocity $H^{1}$-error estimate than the standard one (3.6), cf. Remark 2.8. The next result additionally ensures that the two discretizations are actually comparable if one considers the sum of the velocity $H^{1}$-error times viscosity plus the pressure $L^{2}$-error. Thus, in other words, the modifications introduced in (3.14) do not impair the quasi-optimality of (3.6).

Theorem 3.13** (Quasi-optimality).**

For all $\eta>1$, problem (3.14) is a quasi-optimal discretization of (2.1) with constant $C_{\mathrm{qo}}\lesssim\eta^{3}/(\eta-1)$.

Proof.

Denote by $(u,p)$ and $(u_{h},p_{h})$ the solutions of problems (2.1) and (3.14), respectively, with load $f\in H^{-1}(\Omega)^{2}$ and viscosity $\mu>0$. In view of Theorem 3.11, it suffices to bound the pressure error $\|p-p_{h}\|_{L^{2}(\Omega)}$. To this end, let $q_{h}\in\widehat{S}_{\ell-2}^{0}$ be arbitrary and recall that the discrete divergence $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ is given by (2.7). The inf-sup stability of the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair and Proposition 3.4 yield

[TABLE]

For all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$, a comparison of (2.1) and (3.14) entails

[TABLE]

where we have made use again of Proposition 3.4. The last summand in the right-hand side vanishes if we let $q_{h}$ be the $L^{2}$-orthogonal projection of $p$. Hence, invoking Lemma 3.10 and proceeding as in the proof of Theorem 3.11, we infer

[TABLE]

The triangle inequality and Theorem 3.11 conclude the proof. ∎

3.5. Inhomogeneous continuity equation

It is worth having a look at the case when the incompressibility constraint $\operatorname{div}u=0$ of (1.1) is replaced by the inhomogeneous continuity condition $\operatorname{div}u=g$ with $g\in L^{2}_{0}(\Omega)$. The corresponding weak formulation reads as follows: Find $u\in H^{1}_{0}(\Omega)^{2}$ and $p\in L^{2}_{0}(\Omega)$ such that

[TABLE]

A possible extension of the discretization (3.14) with the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair consists in finding $u_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $p_{h}\in\widehat{S}_{\ell-2}^{0}$ such that

[TABLE]

The second equations of (3.20) and (3.21) impose $u\in Z(g)$ and $u_{h}\in Z_{h}^{{ub}}(g)$, respectively, where

[TABLE]

and $\Pi_{\ell-2}$ is the $L^{2}$-orthogonal projection onto $\widehat{S}_{\ell-2}^{0}$.

Lemma 3.10 states that the consistency error in the left hand side of (3.16) vanishes whenever $z_{h}\in Z\cap Z_{h}^{{ub}}$. If, instead, we assume $z_{h}\in Z(g)\cap Z_{h}^{{ub}}(g)$ for some $g\in L^{2}_{0}(\Omega)$ with $g\neq\Pi_{\ell-2}g$, the consistency error may not vanish. In fact, we possibly have $R_{h}^{ub}z_{h}\neq 0$, as a consequence of Proposition 3.4. This suggests that a bound of the consistency error solely in terms of the best approximation $H^{1}$-error to $z_{h}$ by elements of $Z(g)$ is likely not possible. Therefore, we do not expect that the discrete velocity $u_{h}$ solving (3.21) is a near-best approximation of the analytical velocity in $(\mathring{S}_{\ell}^{1})^{2}$, with respect to the $H^{1}$-norm.

Still, combining the equivalence (3.11) and the $L^{2}$-orthogonality of $\Pi_{\ell-2}$, we obtain the following generalization of Lemma 3.10

[TABLE]

for all $z_{h}\in Z_{h}^{{ub}}(g)$ and $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$, with $g\in L^{2}_{0}(\Omega)$. Apart from the additional term in the right-hand side of this estimate, the technique in the proof of Theorem 3.11 can be still applied, with the help of [8, Proposition 5.1.3], and we finally derive

[TABLE]

for any fixed $\eta>1$. Similarly as in (1.3), here the approximation power of the discrete pressure space in the $L^{2}$-norm may impair the velocity $H^{1}$-error, because the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair is unbalanced. We confirm this suspicion by means of a numerical experiment in section 5.4. Still, we remark that this estimate, unlike (1.3), is pressure robust, i.e. independent of the analytical pressure. A corresponding bound of the pressure error can be derived arguing as in the proof of Theorem 3.13.

The nonconforming discretization proposed in section 4.1 has the remarkable property that the consistency error can always be bounded solely in terms of the best approximation $H^{1}$-error to the analytical velocity; cf. Remark 4.3. Therefore, in that case, we achieve quasi-optimality and pressure robustness even if an inhomogeneous continuity condition is imposed.

4. Generalizations of the paradigmatic discretization

The idea illustrated in the previous section can be generalized in various directions. An immediate observation is that the same construction applies to any other conforming and inf-sup stable pair $V_{h}/Q_{h}$ such that

• $(i)$

$\widehat{S}_{0}^{0}$ is a subset of $Q_{h}$ and

• $(ii)$

the discrete divergence $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ can be computed element-wise.

The first condition is needed in Proposition 3.4 to ensure that the smoothing operator $E_{h}^{ub}$ fulfills (3.1). The second one guarantees that the divergence correction $R^{ub}$ can be computed element-wise. As a consequence, the proposed discretization is computationally feasible, cf. Remark 3.6. Conditions (i) and (ii) are verified, for instance, by the following generalization of the $\mathbb{P}_{\ell}/\mathbb{P}_{\ell-2}$ pair

[TABLE]

where $d\leq k\leq\ell$ and $d\in\{2,3\}$. Another possibility is to consider the conforming Crouzeix-Raviart pairs described in [8, Sections 8.6.2 and 8.7.2]. Stable pairs with continuous pressure, i.e. $Q_{h}\subseteq C^{0}(\Omega)$, do not fulfill (i), while (ii) is violated, for instance, by the modified Hood-Taylor pairs of Boffi et al. [9].

We now aim at addressing more substantial generalizations. We mainly focus on the necessary modifications and, in particular, we omit all proofs that are similar to the ones in the previous section.

4.1. Nonconforming pairs

Assume that $V_{h}/Q_{h}$ is a nonconforming pair, i.e. $V_{h}\nsubseteq H^{1}_{0}(\Omega)^{d}$. In this case, it does not seem appropriate to define the smoothing operator $E_{h}$ as in (3.10), because of the condition $E_{h}(V_{h})\subseteq H^{1}_{0}(\Omega)^{d}$. A possible fix for this problem is to replace $v_{h}$ with $M_{h}v_{h}$, where $M_{h}:V_{h}\to H^{1}_{0}(\Omega)^{d}$ is a linear operator. To make sure that a counterpart of Proposition 3.4 holds, we require that $\operatorname{div}M_{h}v_{h}$ has element-wise the same mean as $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}v_{h}$ for all $v_{h}\in V_{h}$. Therefore, we resort to a element-wise ”mean mass preserving” operator; cf. Proposition 4.1.

As before, we illustrate this idea by means of a model example, namely the two-dimensional nonconforming Crouzeix-Raviart pair of degree $\ell\geq 2$. We do not consider the lowest-order case $\ell=1$, as it is rather specific and it is already covered by [35], cf. Remark 4.2. A similar technique can be applied, for instance, with the modified Crouzeix-Raviart pairs of [27] or with the three-dimensional generalizations of the Kouhia-Stenberg pair from [21]. The original two-dimensional pair of Kouhia and Stenberg [22] can be treated as indicated in Remark 4.2.

Let the mesh $\mathcal{M}$ be as in section 3 and denote by $\mathcal{F}$ the faces of $\mathcal{M}$. A subscript to $\mathcal{F}$ indicates that we consider only those faces that are contained in the set specified by the subscript. We orient each interior face $F\in\mathcal{F}_{\Omega}$ with a normal unit vector $n_{F}$. We denote by $\llbracket\cdot\rrbracket_{|F}$ the jump on $F$ in the direction of $n_{F}$. For boundary faces $F\in\mathcal{F}_{\partial\Omega}$, we orient $n_{F}$ so that it points outside $\Omega$ and let $\llbracket\cdot\rrbracket_{|F}$ coincide with the trace on $F$, cf. [15, Section 1.2.3]. We use the subscript $\mathcal{M}$ to indicate the broken version of a differential operator on $\mathcal{M}$. For instance, the broken gradient of an element-wise $H^{1}$-function $v$ is given by $(\operatorname{\operatorname{\nabla}_{\mathcal{M}}}v)_{|K}:=\operatorname{\nabla}(v_{|K})$ for all $K\in\mathcal{M}$.

The nonconforming Crouzeix-Raviart space of degree $\ell\in\mathbb{N}$ on $\mathcal{M}$, with homogeneous boundary conditions, can be defined as follows

[TABLE]

Notice that the integral $\int_{F}v$ is well-defined for all $v\in\mathring{CR}_{\ell}$ and $F\in\mathcal{F}$ and vanishes if $F\in\mathcal{F}_{\partial\Omega}$. Yet, the jumps on mesh faces are not vanishing in general.

We assume hereafter $\ell\geq 2$. The two-dimensional nonconforming Crouziex-Raviart pair of degree $\ell$ is

[TABLE]

Results concerning the inf-sup stability can be found in [6, 13, 16]. Since the broken divergence $\operatorname{\operatorname{div}_{\mathcal{M}}}$ maps $V_{h}$ into $Q_{h}$, it coincides with the discrete divergence from (2.7), i.e. $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}=\operatorname{\operatorname{div}_{\mathcal{M}}}$. We measure the velocity error in the broken $H^{1}$-norm, augmented with scaled jumps. Thus, in the notation of section 2, we set

[TABLE]

where $h_{F}$ is the diameter of $F$. An equivalent alternative would be to consider only the broken $H^{1}$-norm. Both options extend the $H^{1}$-norm to $H^{1}_{0}(\Omega)^{2}+(\mathring{CR}_{\ell})^{2}$.

Let $\mathcal{V}_{\ell,\Omega}$ be the set of interior Lagrange nodes of degree $\ell$ in $\mathcal{M}$. For all $\nu\in\mathcal{V}_{\ell,\Omega}$, we denote by $\Phi_{\ell}^{\nu}$ the Lagrange basis function of $\mathring{S}_{\ell}^{1}$ associated with the evaluation at $\nu$, i.e. $\Phi_{\ell}^{\nu}(\nu^{\prime})=\delta_{\nu\nu^{\prime}}$ for all $\nu^{\prime}\in\mathcal{V}_{\ell,\Omega}$. Fix also an element $K_{\nu}\in\mathcal{M}$ with $\nu\in K_{\nu}$. We define a ”simplified nodal averaging” operator $A_{h}^{cr}:(\mathring{CR}_{\ell})^{2}\to(\mathring{S}_{\ell}^{1})^{2}$ by

[TABLE]

Next, let $m_{F}$ be the midpoint of any interior face $F\in\mathcal{F}_{\Omega}$. Consider the bubble function $\Phi_{2}^{F}:=3(2\left\lvert F\right\rvert)^{-1}\Phi_{2}^{m_{F}}$, where $\Phi_{2}^{m_{F}}$ is the Lagrange basis function of $\mathring{S}_{2}^{1}$ associated with the evaluation at $m_{F}$. The normalization implies $\int_{F^{\prime}}\Phi_{2}^{F}=\delta_{FF^{\prime}}$ for all $F^{\prime}\in\mathcal{F}$, according to the Simpson quadrature formula. We introduce a ”bubble” operator $B_{h}^{cr}:(\mathring{CR}_{\ell})^{2}\to(\mathring{S}_{\ell}^{1})^{2}$ by

[TABLE]

We combine $A_{h}^{cr}$ and $B_{h}^{cr}$ to obtain the announced element-wise mean mass preserving operator $M_{h}^{cr}$. Roughly speaking, we use $B_{h}^{cr}$ to enforce the first part of (4.2) below, while $A_{h}^{cr}$ is responsible for the second part.

Proposition 4.1** (Element-wise mean mass preserving operator).**

The linear operator $M_{h}^{cr}:(\mathring{CR}_{\ell})^{2}\to(\mathring{S}_{\ell}^{1})^{2}$ given by

[TABLE]

is such that

[TABLE]

for all $v_{h}\in(\mathring{CR}_{\ell})^{2}$ and $K\in\mathcal{M}$.

Proof.

Let $v_{h}\in(\mathring{CR}_{\ell})^{2}$ and $F^{\prime}\in\mathcal{F}_{\Omega}$ be given. The normalization of the functions $\{\Phi_{2}^{F}\}_{F\in\mathcal{F}_{\Omega}}$ reveals

[TABLE]

The same identities hold also for boundary faces $F^{\prime}\in\mathcal{F}_{\partial\Omega}$, in view of the boundary conditions in $\mathring{CR}_{\ell}$ and $\mathring{S}_{\ell}^{1}$. Rearranging terms, we obtain $\int_{F^{\prime}}M_{h}^{cr}v_{h}=\int_{F^{\prime}}v_{h}$ for all $F^{\prime}\in\mathcal{F}$. Then, for all $K\in\mathcal{M}$, the Gauss theorem yields the first part of (4.2)

[TABLE]

A detailed proof of the second part of (4.2) can be found in [34, Section 3], where a similar, actually more involved, operator is considered. For this reason, we only sketch the proof. Let $K\in\mathcal{M}$ be given. Owing to the triangle inequality, we initially bound $\|\operatorname{\nabla}(v_{h}-A_{h}^{cr}v_{h})\|_{L^{2}(K)}$ and $\|\operatorname{\nabla}B_{h}^{cr}(v_{h}-A_{h}^{cr}v_{h})\|_{L^{2}(K)}$. The scaling of the functions $\{\Phi_{2}^{F}\}_{F\in\mathcal{F}_{\partial K}}$ and the trace inequality imply

[TABLE]

where $h_{K}$ is the diameter of $K$. Next, for all $\nu\in\mathcal{V}_{\ell,K}$, we have $v_{h|K}(\nu)=A_{h}^{cr}v_{h}(\nu)$ if $\nu\in\mathrm{int}(K)$, otherwise $\left\lvert v_{h|K}(\nu)-A_{h}^{cr}v_{h}(\nu)\right\rvert\lesssim\sum_{F\ni\nu}h_{F}^{-1/2}\|\llbracket v_{h}\rrbracket\|_{L^{2}(F)}$, where $F$ varies in $\mathcal{F}$. This estimate and the scaling of the Lagrange basis functions entail that the right-hand side of (4.3) is bounded by $\sum_{F\cap K\neq\emptyset}h_{F}^{-1/2}\|\llbracket v_{h}\rrbracket\|_{L^{2}(F)}$. Squaring and summing over all $K\in\mathcal{M}$, we finally obtain

[TABLE]

We conclude recalling the definition of the norm $\|\cdot\|_{{cr}}$. ∎

According to the first part of (4.2), we can now construct a smoothing operator similarly to $E_{h}^{ub}$ in Proposition 3.4. Recalling the local operators $R_{\ell}^{K}$ introduced in section 3.2, we define $E_{h}^{cr}:(\mathring{CR}_{\ell})^{2}\to H^{1}_{0}(\Omega)^{2}$ by

[TABLE]

Owing to the identity $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}=\operatorname{\operatorname{div}_{\mathcal{M}}}$, we see that $E_{h}^{cr}$ fulfills condition (3.1), as a consequence of Propositions 3.2 and 4.1. Moreover, the stability of the operators $R_{\ell}^{K}$ and the second part of (4.2) provide a strengthened counterpart of (3.11) in that, for all $v_{h}\in(\mathring{CR}_{\ell})^{2}$, we have

[TABLE]

Next, inspired by the definition of $a_{h}^{ub}$ in (3.13) as well as by identity (3.15), we introduce the following bilinear form $a_{h}^{cr}$ on $(\mathring{CR}_{\ell})^{2}$

[TABLE]

where $R_{h}^{cr}:=(E_{h}^{cr}-\mathrm{Id})$ and $\eta>1$ is a penalty parameter. The above-mentioned properties of $E_{h}^{cr}$ imply that the necessary conditions in Lemmas 2.6 and 2.7 are fulfilled if we set $a_{h}=a_{h}^{cr}$ and $E_{h}=E_{h}^{cr}$. In this setting, the abstract discretization (2.2) reads as follows: Find $u_{h}\in(\mathring{CR}_{\ell})^{2}$ and $p_{h}\in\widehat{S}_{\ell-1}^{0}$ such that

[TABLE]

Similarly as $a_{h}^{ub}$ in Lemma 3.5, the form $a_{h}^{cr}$ is coercive on $(\mathring{CR}_{\ell})^{2}$, for $\eta>1$, with constant $\geq(1-\eta^{-1})$. Moreover, in view of (4.5), we can estimate the consistency error of (4.6) by the following counterpart of Lemma 3.10

[TABLE]

for all $w_{h},v_{h}\in(\mathring{CR}_{\ell})^{2}$. Hence, we conclude that (4.6) is a quasi-optimal and pressure robust discretization of (2.1) in the norm $\|\cdot\|_{{cr}}$ and the constant $\mathrm{C_{\mathrm{qopr}}}$ from Definition 2.5 solely depends on $\eta$ and the shape parameter of $\mathcal{M}$.

Whenever the pair $(\mathring{CR}_{\ell})^{2}/\widehat{S}_{\ell-1}^{0}$ is inf-sup stable, an estimate of the pressure $L^{2}$-error, only in terms of the best approximation errors to the analytical velocity and the analytical pressure, can also be established similarly as in Theorem 3.13. Thus, problem (4.6) is also a quasi-optimal discretization of (2.1).

Locally supported basis functions of $\mathring{CR}_{\ell}$ are described in [5, section 3]. With this basis and the standard nodal basis of $S_{\ell-1}^{0}$, we see that (4.6) is computationally feasible in the sense of Remark 2.4, cf. Remark 3.6.

Remark 4.2* (The pair $(\mathring{CR}_{1})^{2}/\widehat{S}_{0}^{0}$).*

In principle, the approach described for $\ell\geq 2$ applies also with $\ell=1$, up to observing that $R_{2}^{K}$ (and not $R_{1}^{K}$) should be used in (4.4). The point is that, in this case, an element-wise integration by parts and the identity $\int_{F^{\prime}}E_{h}^{cr}v_{h}=\int_{F^{\prime}}M_{h}^{cr}v_{h}=\int_{F^{\prime}}v_{h}$, with $F^{\prime}\in\mathcal{F}$, reveal $\int_{\Omega}\operatorname{\operatorname{\nabla}_{\mathcal{M}}}w_{h}\colon\operatorname{\nabla}R_{h}^{cr}v_{h}=0$ for all $w_{h},v_{h}\in(\mathring{CR}_{1})^{2}$. Hence, the form $a_{h}^{cr}$ is given by $a_{h}^{cr}(w_{h},v_{h})=\int_{\Omega}\operatorname{\operatorname{\nabla}_{\mathcal{M}}}w_{h}\colon\operatorname{\operatorname{\nabla}_{\mathcal{M}}}v_{h}+\eta\int_{\Omega}\operatorname{\operatorname{\nabla}_{\mathcal{M}}}R_{h}^{cr}w_{h}\colon\operatorname{\operatorname{\nabla}_{\mathcal{M}}}R_{h}^{cr}v_{h}$, showing that the penalization is actually not needed. Setting $\eta=0$ annihilates the consistency error and corresponds to the discretization proposed in [35].

Remark 4.3* (Inhomogeneous continuity equation).*

The infimum in the right-hand side of (4.7) is taken over $H^{1}_{0}(\Omega)^{2}$ and not only over $Z$, unlike Lemma 3.10. This prevents the issue pointed out in section 3.5. Therefore, the nonconforming Crouzeix-Raviart pair can be used to design a quasi-optimal and pressure robust discretization of problem (3.20) with the inhomogeneous continuity condition $g\neq 0$.

4.2. Conforming pairs with continuous pressure

Another class of pairs still not covered by our discussion are conforming pairs with continuous pressure. In fact, the following observations obstruct the construction of a smoothing operator as indicated in Proposition 3.4.

• $(i)$

Since $\widehat{S}_{0}^{0}$ is not a subspace of $Q_{h}$, the identity $\int_{K}\operatorname{\underline{\operatorname{div}}_{\textit{h}}}v_{h}=\int_{K}\operatorname{div}v_{h}$ may fail to hold for some $v_{h}\in V_{h}$ and $K\in\mathcal{M}$.

• $(ii)$

The computation of $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$ is likely unfeasible in the sense of Remark 2.4.

Item (i) entails that we cannot correct the divergence element-wise by means of the operators $R_{\ell}^{K}$ from section 3.2. The shape functions of the lowest-order continuous space $S_{1}^{1}$ suggest to work on patches of elements sharing a vertex, instead. Item (ii) further indicates that we should never require a direct computation of $\operatorname{\underline{\operatorname{div}}_{\textit{h}}}$. The construction of a quasi-optimal and pressure robust discretization of the Stokes equations is still possible under these constraints, but it is more involved than the ones in the previous sections. We mainly adapt ideas by Lederer et al. [23].

As an example, we let the mesh $\mathcal{M}$ be as in section 3 and consider the two-dimensional Hood-Taylor pair

[TABLE]

with $\ell\geq 2$. The inf-sup condition (2.3) holds with $\beta^{-1}\leq c$ under mild assumptions on $\mathcal{M}$, see [7]. The discrete divergence coincides with the $L^{2}$-orthogonal projection of the analytical divergence onto $\widehat{S}_{\ell-1}^{1}$. We denote by $Z_{h}^{ht}$ the discrete kernel.

Let $\mathcal{V}:=\mathcal{V}_{1}$ denote the set of all vertices of $\mathcal{M}$. For each $\nu\in\mathcal{V}$, let $\Phi_{1}^{\nu}$ be the Lagrange basis function of $S_{1}^{1}$ associated with the evaluation at $\nu$, i.e. $\Phi_{1}^{\nu}(\nu^{\prime})=\delta_{\nu\nu^{\prime}}$ for all $\nu^{\prime}\in\mathcal{V}$. Recall that $\Phi_{1}^{\nu}$ is supported on the patch $\omega_{\nu}:=\{K\in\mathcal{M}\mid\nu\in K\}$. Consider the barycentric refinement $\mathcal{M}_{\nu}$ of $\omega_{\nu}$, i.e. the mesh obtained connecting the vertices and the barycenter of any triangle in $\omega_{\nu}$, cf. Figure 4.1. The space $S_{\ell}^{0}(\mathcal{M}_{\nu})$ and the subspaces

[TABLE]

are defined on $\mathcal{M}_{\nu}$ analogously to $S_{\ell}^{0}$ in (3.2) and $\mathring{S}_{\ell}^{1}$ and $\widehat{S}_{\ell-1}^{0}$ in (3.3), respectively. The element-wise local Lagrange interpolant $I_{\ell}^{\nu}:S_{\ell}^{0}(\mathcal{M}_{\nu})\to S_{\ell-1}^{0}(\mathcal{M}_{\nu})$ is given by

[TABLE]

where $\mathcal{V}_{\ell-1,K}$ is the set of Lagrange nodes of degree $\ell-1$ in $K$ and $\Phi_{\ell-1}^{\nu^{\prime},K}$ is the Lagrange basis function of $\mathbb{P}_{\ell}(K)$ associated with the evaluation at $\nu^{\prime}$ and extended to zero outside $K$. Consider also the simplified local averaging $A_{\ell}^{\nu}:S_{\ell}^{0}(\mathcal{M}_{\nu})\to S_{\ell-1}^{1}$

[TABLE]

where $K_{\nu}\in\mathcal{M}$ is a fixed element such that $\nu\in K_{\nu}$ and $v$ is extended to zero outside $\omega_{\nu}$. As before, $\Phi_{\ell-1}^{\nu}$ denotes the Lagrange basis function of $\mathring{S}_{\ell-1}^{1}$ associated with the evaluation at $\nu$.

We are now ready to define the operators $R_{\ell}^{\nu}:L^{2}(\Omega)\to H^{1}_{0}(\Omega)^{2}$ that will be used to correct the divergence in each patch $\omega_{\nu}$, $\nu\in\mathcal{V}$. Here $R_{\ell}^{\nu}$ plays the same role as $R_{\ell}^{K}$ in section 3.2. Given $q\in L^{2}(\Omega)$, let $u_{\nu}=u_{\nu}(q)\in\mathring{S}_{\ell}^{1}(\mathcal{M}_{\nu})^{2}$ and $p_{\nu}=p_{\nu}(q)\in\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{\nu})$ be such that

[TABLE]

This problem is uniquely solvable, according to [19, Corollary 6.2]. Then, we set

[TABLE]

Remark 4.4* (Local problems).*

The use of the barycentric refinement $\mathcal{M}_{\nu}$ is a main difference compared to [23]. This ensures that the pair $\mathring{S}_{\ell}^{1}(\mathcal{M}_{\nu})^{2}/\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{\nu})$ is inf-sup stable. In fact, it is known that the stability of the Scott-Vogelius pair on $\omega_{\nu}$ (without the barycentric refinement) may be impaired if $\nu$ is a singular or nearly singular vertex, see [32]. The partition of unity $\{\Phi_{1}^{\nu}\}_{\nu\in\mathcal{V}}$ and the interpolants $\{I_{\ell}^{\nu}\}_{\nu\in\mathcal{V}}$ account for the overlapping of the patches, while the averaging operators $\{A_{\ell}^{\nu}\}_{\nu\in\mathcal{V}}$ are used to avoid a direct computation of the discrete divergence in (4.9).

We define a global divergence correction $R_{h}^{ht}:(\mathring{S}_{\ell}^{1})^{2}\to H^{1}_{0}(\Omega)^{2}$

[TABLE]

In contrast to $E_{h}^{ub}$ and $E_{h}^{cr}$ from (3.10) and (4.4), respectively, we now make use of a smoothing operator $E_{h}^{ht}$ which is not guaranteed to be divergence-preserving, i.e. (3.1) may fail to hold. We shall see, however, that it still satisfies the necessary conditions in Lemmas 2.6 and 2.7. In the following proposition we only prove a basic stability estimate, for the sake of simplicity.

Proposition 4.5** (Smoothing operator for the Hood-Taylor pair).**

The linear operator $E_{h}^{ht}:(\mathring{S}_{\ell}^{1})^{2}\to H^{1}_{0}(\Omega)^{2}$ given by

[TABLE]

satisfies (2.12a) and (2.13a) and is such that, for all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$,

[TABLE]

Proof.

For all $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $q_{h}\in\widehat{S}_{\ell-1}^{1}$, we have

[TABLE]

The first identity follows from the second equation of (4.8), which actually holds for all $q_{\nu}$ in $S_{\ell-1}^{0}(\mathcal{M}_{\nu})$ (and not only in $\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{\nu})$), as both sides vanish if $q_{\nu}$ is constant. To check the second identity, observe that $A_{\ell}^{\nu}(q_{h}\Phi_{1}^{\nu})=I_{\ell}^{\nu}(q_{h}\Phi_{1}^{\nu})$ for all $\nu\in\mathcal{V}$, due to the continuity of $q_{h}\Phi_{1}^{\nu}$. Thus, we derive the identity

[TABLE]

showing that condition (2.12a) holds. Next, let $z_{h}\in Z^{ht}$ be given and consider $q_{h}=\operatorname{div}E_{h}^{ht}z_{h}$. Recall that $\{\Phi_{1}^{\nu}\}_{\nu\in\mathcal{V}}$ is a partition of unity and extend $I_{\ell}^{\nu}(q_{h}\Phi_{1}^{\nu})$ to zero outside $\omega_{\nu}$. We infer $\sum_{\nu\in\mathcal{V}}I_{\ell}^{\nu}(q_{h}\Phi_{1}^{\nu})=q_{h}$. Then, since $z_{h}$ is discretely divergence-free, we have

[TABLE]

This reveals $\operatorname{div}E_{h}^{ht}z_{h}=0$ and confirms that condition (2.13a) holds. Finally, owing to the stability of $I_{\ell}^{\nu}$ and $A_{\ell}^{\nu}$ in the $L^{2}(\omega_{\nu})$-norm, we infer

[TABLE]

for all $\nu\in\mathcal{V}$ and $v_{h}\in(\mathring{S}_{\ell}^{1})^{2}$. This entails $\|\operatorname{\nabla}\operatorname{div}R_{\ell}^{\nu}v_{h}\|_{L^{2}(\omega_{\nu})}\lesssim\|\operatorname{div}v_{h}\|_{L^{2}(\omega_{\nu})}$, owing to [8, Corollary 4.2.1] and the inf-sup stability of the pair $\mathring{S}_{\ell}^{1}(\mathcal{M}_{\nu})^{2}/\widehat{S}_{\ell-1}^{0}(\mathcal{M}_{\nu})$ stated in [19, Corollary 6.2]. The definition of $R_{h}^{ht}$ in (4.9) then implies

[TABLE]

for all $K\in\mathcal{M}$, where $K^{\prime}$ varies in $\mathcal{M}$. We conclude summing over all elements of $\mathcal{M}$ and recalling the definition of $E_{h}^{ht}$. ∎

Next, for $\eta>1$, we introduce the following bilinear form on $(\mathring{S}_{\ell}^{1})^{2}$

[TABLE]

The abstract discretization (2.2) with $a_{h}=a_{h}^{ht}$ and $E_{h}=E_{h}^{ht}$ looks for $u_{h}\in(\mathring{S}_{\ell}^{1})^{2}$ and $p_{h}\in\widehat{S}_{\ell-1}^{1}$ such that

[TABLE]

This discretization is computationally feasible in the sense of Remark 2.4, cf. Remark 3.6. Yet, the implementation is more costly than the one of (3.14) and (4.6) because, in general, we cannot resort to one reference configuration for the solution of the local problems (4.8). The error analysis of (4.11) proceeds almost verbatim as in section 3.4, with the help of Proposition 4.5. The only remarkable difference is that estimate (3.19) in the proof of Theorem 3.13 should be replaced by the weaker one $\|p_{h}-q_{h}\|_{L^{2}(\Omega)}\lesssim\mu\eta\|\operatorname{\nabla}(u-u_{h})\|_{L^{2}(\Omega)}+\|p-q_{h}\|_{L^{2}(\Omega)}$, because identity (3.1) may fail to hold.

5. Numerical experiments with the unbalanced $\mathbb{P}_{2}/\mathbb{P}_{0}$ pair

In this section we restrict our attention to the two-dimensional Stokes equations, with unit viscosity, posed in the unit square. In the notation of section 2, this corresponds to

[TABLE]

We investigate numerically the new discretization (3.14), based on the unbalanced $\mathbb{P}_{2}/\mathbb{P}_{0}$ pair, i.e.

[TABLE]

If not specified differently, the penalty parameter is set to

[TABLE]

We shall consider the following families $(\mathcal{M}_{N}^{D})_{N\in\mathbb{N}_{0}}$ and $(\mathcal{M}_{N}^{C})_{N\in\mathbb{N}_{0}}$ of triangular meshes of $\Omega$. For $N\in\mathbb{N}_{0}$, we divide $\Omega$ into $2^{N}\times 2^{N}$ identical squares, with edges parallel to the $x_{1}$- and $x_{2}$-axis and with area $2^{-2N}$. We obtain the ”diagonal mesh” $\mathcal{M}_{N}^{D}$ dividing each square by the diagonal with positive slope. Similarly, we obtain the ”crisscross mesh” $\mathcal{M}_{N}^{C}$ drawing both diagonals of each square, cf. Figure 5.1. All experiments have been implemented in ALBERTA 3.0 [20, 31].

5.1. Smooth solution

To illustrate the quasi-optimality and pressure robustness of the new $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization, we first consider a test case with smooth analytical solution, given by

[TABLE]

where $\operatorname{curl}(w):=(\partial_{2}w,-\partial_{1}w)$. We compare the performances of the standard $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization (3.6) and the new one (3.14) on the crisscross meshes $\mathcal{M}_{N}^{C}$ with $N=0,\dots,8$. Figure 5.2 displays the respective balances of velocity $H^{1}$-error and pressure $L^{2}$-error versus $\#\mathcal{M}_{N}^{C}$, that is the number of triangles in the mesh.

We first observe that the pressure $L^{2}$-errors of both discretizations behave quite similarly and converge to zero with the maximum decay rate $(\#\mathcal{M}_{N}^{C})^{-0.5}$. The velocity $H^{1}$-error of the standard discretization converges to zero with the same decay rate, as suggested by estimate (1.3), according to the approximation power of the discrete pressure space in the $L^{2}$-norm. Note, however, that such rate is suboptimal with respect to the approximation power of the discrete velocity space in the $H^{1}$-norm. In contrast, the velocity $H^{1}$-error of the new discretization exhibits the maximum decay rate $(\#\mathcal{M}_{N}^{C})^{-1}$, as predicted by Theorem 3.11. The next experiments are intended to highlight some of the ingredients that contribute to make this optimal-order convergence possible.

5.2. Composite numerical quadrature

The evaluation of the duality $\left\langle f,E_{h}v_{h}\right\rangle$, $v_{h}\in(\mathring{S}_{2}^{1})^{2}$, in the new $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization requires, in particular, the evaluation of $\left\langle f,\widetilde{v}_{h}\right\rangle$ for test functions $\widetilde{v}_{h}$ that are element-wise quadratic on the barycentric refinement of the mesh at hand. This suggests that, for each triangle $K$ in the mesh, a composite quadrature rule, based on the barycentric refinement of $K$, should be used. If one, instead, uses a standard quadrature rule in $K$, the resulting quadrature error could be not negligible, due to the low regularity of $\widetilde{v}_{h}$. Moreover, since the quadrature error is potentially not pressure robust, as pointed out in [26, section  6.2], this may even affect the decay rate of the velocity $H^{1}$-error.

To illustrate such effect, we consider a test case with analytical solution

[TABLE]

For $\alpha\in\{1,10^{3}\}$, we apply the new $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization on the crisscross meshes $\mathcal{M}_{N}^{C}$ with $N=0,\dots,8$. We assemble the right-hand side both with a composite and a standard quadrature rule of degree $6$. For $N=4,\dots,8$, the corresponding velocity $H^{1}$-errors are reported in Table 5.1. In each case, we compute also the so-called experimental order of convergence (EOC), defined as

[TABLE]

where $e_{N}$ denotes the $H^{1}$-error on $\mathcal{M}_{N}^{C}$.

When the composite quadrature rule is applied, the results seem insensitive to the parameter $\alpha$ and we observe the maximum decay rate $(\#\mathcal{M}_{N}^{C})^{-1}$. In contrast, the use of the standard quadrature rule impairs the pressure robustness stated in Theorem 3.11. In fact, for sufficiently large $N$, the velocity $H^{1}$-error is essentially proportional to $\alpha$ and exhibits the suboptimal decay rate $(\#\mathcal{M}_{N}^{C})^{-0.5}$.

5.3. Locking

As mentioned in Remark 3.8, the bilinear form $a_{h}^{ub}$ in the new $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization has the same structure as the DG-SIP form of [1]. Still, one main difference is that Lemma 3.5 ensures the coercivity of the former for any penalty $\eta>1$ (and not only for sufficiently large $\eta$). Moreover, the coercivity constant is $\geq 0.5$ for $\eta=2$. Having an explicit and safe choice of the penalty parameter is particularly useful in this context, because we may have locking for large $\eta$, in view of Remark 3.12.

To illustrate this, we consider a test case with analytical solution

[TABLE]

We apply the new $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization for $\eta\in\{2,32,512\}$ both on diagonal meshes $\mathcal{M}_{N}^{D}$ and on crisscross meshes $\mathcal{M}_{N}^{D}$, with $N=0,\dots,7$. The velocity $H^{1}$-errors displayed in the right part of Figure 5.3 indicate that the new discretization is robust with respect to $\eta$ on crisscross meshes. This follows from the fact that condition (3.18) in Remark 3.12 holds for such meshes, as a consequence of [30, Theorem 4.3.1]. In contrast, adopting the terminology of [3], we observe on the left part of Figure 5.3 locking of order $(\mathcal{M}_{N}^{D})^{1/2}$ when diagonal meshes are used.

5.4. Inhomogeneous continuity equation

We finally point out that the quasi-optimality and pressure robustness of the new $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization, as stated in Theorem 3.11, hinges on the homogeneity of the continuity equation in the Stokes problem (2.1), cf. section 3.5.

To see this, we consider the more general problem (3.20) and approximate the analytical solution

[TABLE]

on the crisscross meshes $\mathcal{M}_{N}^{C}$ with $N=0,\dots,8$. Note, in particular, that $\operatorname{div}u$ is not element-wise constant on $\mathcal{M}_{N}^{C}$.

Comparing the velocity $H^{1}$-errors of the standard $\mathbb{P}_{2}/\mathbb{P}_{0}$ discretization (3.6) and the new one (3.14), we see that the former is slightly smaller than the latter and that both errors converge to zero with decay rate $(\mathcal{M}_{N}^{C})^{-0.5}$; cf. Figure 5.4. This confirms that inequality (3.22) captures the correct behavior of the new discretization. Thus, for this problem, we expect that the new discretization performs significantly better than the standard one only in case of large pressure $L^{2}$-errors.

Acknowledgements

We wish to thank Rüdiger Verfürth for reading some preliminary versions of this manuscript and for suggesting several improvements in the presentation.

Funding

The authors gratefully acknowledge partial support by the DFG research grant KR 3984/5-1 “Convergence Analysis for Adaptive Discontinuous Galerkin Methods”.