A unified framework for a posteriori error estimation for the Stokes problem

Abstract

In this paper, a unified framework for a posteriori error estimation for the Stokes problem is developed. It is based on $[H^1_0(\Omega )]^d$-conforming velocity reconstruction and $\underline{\boldsymbol{H}}(\mathrm{div},\Omega )$-conforming, locally conservative flux (stress) reconstruction. It gives guaranteed, fully computable global upper bounds as well as local lower bounds on the energy error. In order to apply this framework to a given numerical method, two simple conditions need to be checked. We show how to do this for various conforming and conforming stabilized finite element methods, the discontinuous Galerkin method, the Crouzeix–Raviart nonconforming finite element method, the mixed finite element method, and a general class of finite volume methods. The tools developed and used include a new simple equilibration on dual meshes and the solution of local Poisson-type Neumann problems by the mixed finite element method. Numerical experiments illustrate the theoretical developments.

Appendix

1.1 Characterization of the inf–sup constants

In this section we will show a proof of Lemma 3.1, following the ideas of [45]. We start by the following well-known result, cf. [37, 54]:

Lemma 9.1

(Characterization of the $\inf$–$\sup$ constant) The $\inf$–$\sup$ constant $\beta$ of (2.4) is the square root of the smallest eigenvalue to the following generalized eigenvalue problem

$$\begin{aligned} a(\mathbf{{u}},\mathbf{{v}}) + b(\mathbf{{v}},p)&= 0\quad \forall \mathbf{{v}}\in \mathbf{{V}}, \end{aligned}$$

(9.1a)

$$\begin{aligned} b(\mathbf{{u}},q)&= -\lambda (p,q) \quad \forall q \in Q. \end{aligned}$$

(9.1b)

Proof

Define the following operators

$$\begin{aligned} A&:=-\Delta : \mathbf{{V}}\rightarrow \mathbf{{V}}^{*},\\ B&:=\nabla : Q \rightarrow \mathbf{{V}}^{*}, \end{aligned}$$

so that it holds

$$\begin{aligned} a(\mathbf{{v}},\mathbf{{v}})= (A^{1/2}\mathbf{{v}},A^{1/2}\mathbf{{v}}) \quad \text{ and} \quad b(\mathbf{{v}},q)=(B^{*}\mathbf{{v}},q) \end{aligned}$$

for all $\mathbf{{v}}\in \mathbf{{V}}$ and all $q \in Q$. With this notation, we have

$$\begin{aligned} \frac{b(\mathbf{{v}},q)}{\Vert \nabla \mathbf{{v}}\Vert \, \Vert q\Vert } =\frac{(B^{*}\mathbf{{v}},q)}{\Vert A^{1/2}\mathbf{{v}}\Vert \, \Vert q\Vert }. \end{aligned}$$

Substituting $\mathbf{{z}}:=A^{1/2}\mathbf{{v}}$ gives

$$\begin{aligned} \frac{b(\mathbf{{v}},q)}{\Vert \nabla \mathbf{{v}}\Vert \, \Vert q\Vert } =\frac{(B^{*}A^{-1/2}\mathbf{{z}},q)}{\Vert \mathbf{{z}}\Vert \, \Vert q\Vert } =\frac{(\mathbf{{z}}, A^{-1/2} Bq)}{\Vert \mathbf{{z}}\Vert \, \Vert q\Vert } \end{aligned}$$

and hence the supremum is obtained by choosing

$$\begin{aligned} \mathbf{{z}}= \frac{A^{-1/2}Bq}{\Vert A^{-1/2}Bq \Vert } \end{aligned}$$

and we come to

$$\begin{aligned} \beta = \inf _{q \in Q} \frac{\Vert A^{-1/2} Bq\Vert }{ \Vert q\Vert } . \end{aligned}$$

Squaring gives

$$\begin{aligned} \beta ^{2} = \inf _{q \in Q} \frac{\Vert A^{-1/2} Bq\Vert ^{2}}{ \Vert q\Vert ^{2}} . \end{aligned}$$

This is the Rayleigh quotient for the eigenvalue problem

$$\begin{aligned} B^{*}A^{-1}B p=\lambda p. \end{aligned}$$

(9.2)

Denoting $\mathbf{{u}}=-A^{-1}B p$ this is written as

$$\begin{aligned} A\mathbf{{u}}+ B p&= 0, \\ B^{*}\mathbf{{u}}&= -\lambda p, \end{aligned}$$

i.e., the operator form of (9.1a)–(9.1b). $\square$

We are now ready to prove Lemma 3.1.

Proof

(Proof of Lemma 3.1) In complete analogy to the preceding proof, with $B, A$ replaced by

$$\begin{aligned} \mathcal{B }=\left( \begin{array}{c@{\quad }c} A&B\\ B^{*}&0 \end{array} \right) \quad \text{ and} \quad \mathcal{A }= \left( \begin{array}{c@{\quad }c}A&0\\ 0&\beta ^{2} I\end{array} \right) , \end{aligned}$$

respectively, the $\inf$–$\sup$ constant is the square root of the smallest eigenvalue $\mu$ of

$$\begin{aligned} \mathcal{B }^{*} \mathcal{A }^{-1} \mathcal{B } \mathcal{V } = \mu \mathcal{A } \mathcal{V }. \end{aligned}$$

Written out explicitly, with $\mathcal{V }^T=(\mathbf{{u}},p)^{T}$, this is

$$\begin{aligned} \left( \begin{array}{c@{\quad }c}A&B\\ B^{*}&0 \end{array}\right) \left( \begin{array}{c@{\quad }c}A^{-1}&0\\ 0&\beta ^{-2} I \end{array} \right) \left( \begin{array}{c@{\quad }c}A&B\\ B^{*}&0 \end{array}\right) \left( \begin{array}{c} \mathbf{{u}}\\ p \end{array} \right) = \mu \left( \begin{array}{c@{\quad }c}A&0\\ 0&\beta ^{2} I \end{array} \right) \left( \begin{array}{c} \mathbf{{u}}\\ p \end{array} \right). \end{aligned}$$

From here we see that $\mu = \nu ^{2}$, where $\nu$ is the eigenvalue to

$$\begin{aligned} \left( \begin{array}{c@{\quad }c} A&B \\ B^{*}&0 \end{array}\right) \left( \begin{array}{c} \mathbf{{u}}\\ p \end{array} \right) = \nu \left( \begin{array}{c@{\quad }c}A&0\\ 0&\beta ^{2} I \end{array} \right) \left( \begin{array}{c} \mathbf{{u}}\\ p \end{array} \right). \end{aligned}$$

Explicitly,

$$\begin{aligned} A\mathbf{{u}}+ B p&= \nu A\mathbf{{u}},\\ B^{*}\mathbf{{u}}&= \beta ^{2}\nu p. \end{aligned}$$

From here, we see that $\nu =1$ is an eigenvalue. Suppose next that $\nu \not =1$. Solving for $(\nu -1)\mathbf{{u}}$ in the first equation and substituting in the second one gives

$$\begin{aligned} B^{*}A^{-1}B p=\beta ^{2}\nu (\nu -1)p. \end{aligned}$$

Comparing with (9.2) shows that

$$\begin{aligned} \nu (\nu -1) = 1, \end{aligned}$$

giving

$$\begin{aligned} \nu =\frac{1\pm \sqrt{5}}{2}. \end{aligned}$$

The constant in the stability condition is thus $\min {|\nu |}$, i.e.

$$\begin{aligned} \min \left\{ 1, \frac{\sqrt{5}+1}{2}, \frac{\sqrt{5}-1}{2} \right\} =\frac{\sqrt{5}-1}{2}. \end{aligned}$$

This completes the proof. $\square$

1.2 Equilibration for higher-order conforming and conforming stabilized finite element methods on dual meshes

This appendix concerns conforming and conforming stabilized finite element methods of Sect. 7.2. More precisely, for higher-order continuous pressure elements of Sect. 7.2.2, we show how to, from (7.29), obtain new normal flux functions ${\varvec{\varUpsilon }}_{F}(\mathbf{{u}}_h,p_h)$ for which (7.30) holds. This can be seen as an equivalent of the equilibration procedure of [4] on dual meshes.

Let $D \in \mathcal D ^{\mathrm{int}}_h$, $V$ be the associated vertex, $T \in \mathfrak T _V$, and $i = 1,\ldots ,d$. Denote the contribution to the correction terms of the right-hand side of (7.29) by

$$\begin{aligned}&m_{V,T,i} \nonumber \\&:={} - (\mathbf{{f}}+ \Delta \mathbf{{u}}_h - \nabla p_h, \mathbf e _i)_{T \cap D} + (\mathbf{{f}}+ \Delta \mathbf{{u}}_h - \nabla p_h, {\varvec{\psi }}_{V,i})_{T} \nonumber \\&\quad - \frac{1}{2}\sum _{F \in \mathcal F _T^\mathrm{int}} \langle \left[\!\left[\nabla \mathbf{{u}}_h \mathbf{{n}}_F\right]\!\right], {\varvec{\psi }}_{V,i}\rangle _F + \frac{1}{2} \sum _{F \in \mathcal F _T^\mathrm{int}} \langle \left[\!\left[\nabla \mathbf{{u}}_h \mathbf{{n}}_F\right]\!\right], \mathbf e _i\rangle _{F \cap D}. \end{aligned}$$

(9.3)

We will speak about these quantities as of “normal fluxes” $m_{V,T,i}$. Remark that $\left[\!\left[\nabla \mathbf{{u}}_h \mathbf{{n}}_F\right]\!\right] = 0$ on such sides $F \in \partial \mathcal S _D^\mathrm{int}$ which are not contained in $\partial \mathcal T _h$, cf. Fig. 1. Thus, from (7.29) and the above formula, we have

$$\begin{aligned}&\sum _{F \in \mathcal{F }_D} \langle {\varvec{\varUpsilon }}_{F}(\mathbf{{u}}_h) \mathbf{{n}}_D {\cdot } \mathbf{{n}}_F, \mathbf e _i \rangle _F - (\nabla p_h, \mathbf e _i)_D + (\mathbf{{f}}, \mathbf e _i)_D + \sum _{T \in \mathfrak T _V} m_{V,T,i} = 0, \nonumber \\&\quad i = 1,\ldots ,d, \quad \forall D \in \mathcal D ^{\mathrm{int}}_h. \end{aligned}$$

(9.4)

For the sake of simplicity, let us define $m_{V,T,i}$ in the same way also for $D \in \mathcal D ^{\mathrm{ext}}_h$ and the associated vertex $V$.

Consider a fixed $T \in \mathcal T _h$ and $i = 1,\ldots ,d$. We have associated the normal flux $m_{V_j,T,i}$ to each of the vertices $V_{j}$ of $T$, $j = 1,\ldots ,d + 1$, cf. Fig. 15. We now want to equilibrate the normal fluxes $m_{V_j,T,i}$: the purpose is to associate to each of the sides $F_m \subset T$, $m = 1,\ldots ,d + 1$, $F_m \in \partial \mathcal S _h^\mathrm{int}$ such that $F \subset \partial D$ for some $D \in \mathcal D _h$, a correction normal flux $\upsilon _{F_m,i}$ (in the direction of the fixed normal $\mathbf{{n}}_F$) such that the following holds (we give an example for $d=2$, corresponding to Fig. 15):

$$\begin{aligned} \left(\begin{array}{r@{\quad }r@{\quad }r} 1&1&0 \\ 0&-1&1 \\ -1&0&-1 \\ \end{array}\right) \left(\begin{array}{c} \upsilon _{F_1,i} \\ \upsilon _{F_2,i} \\ \upsilon _{F_3,i}\\ \end{array}\right) = \left(\begin{array}{c} m_{V_1,T,i} \\ m_{V_2,T,i} \\ m_{V_3,T,i} \\ \end{array}\right). \end{aligned}$$

(9.5)

The value $m_{V_1,T,i}$ represents the total normal flux from the element $T \cap D_1$ to the elements $T \cap D_2$ and $T \cap D_3$ (where $D_i$ are the dual volumes associated with the vertices $V_i$). We clearly want to keep this total normal flux but to split it into the side normal fluxes $\upsilon _{F_1,i}$ and $\upsilon _{F_2,i}$; we proceed similarly for $m_{V_2,T,i}$ and $m_{V_3,T,i}$. The essential feature is that the corrections normal fluxes $\upsilon _{F_m,i}$ are univocally defined for each side $F_m$, $m = 1,\ldots ,d + 1$, cf. once again Fig. 15.

It turns out that the system matrix in (9.5) is singular, as the sum of all the row vectors equals zero. It is, however, easy to check that its rank is equal to $d$. Fortunately, the right-hand side in (9.5) is compatible: by the fact that the basis functions ${\varvec{\psi }}_{V_j,i}$ form a partition of unity on the chosen element $T \in \mathcal T _h$,

$$\begin{aligned} \sum _{j = 1}^{d+1} {\varvec{\psi }}_{V_j,i}|_T = \mathbf e _i|_T, \end{aligned}$$

we easily get

$$\begin{aligned} \sum _{j = 1}^{d+1} m_{V_j,T,i} = 0, \end{aligned}$$

$i = 1,\ldots ,d$. Thus, there exists a solution to (9.5). Note that (9.5) is always a system of a fixed small size $(d+1) \times (d+1)$ on each $T \in \mathcal T _h$, for approximations (7.11a)–(7.11b) or (7.12a)–(7.12b) of any order $k$.

Fig. 15

Equilibration of the correction terms inside each triangle

Using $\upsilon _{F_m,i}$ for each $T \in \mathcal T _h$, we can now define new normal flux functions ${\varvec{\varUpsilon }}_{F}(\mathbf{{u}}_h,p_h)$ for sides $F \in \partial \mathcal S _h^\mathrm{int}$ such that $F \subset \partial D$ for some $D \in \mathcal D _h$, in a way that (7.30) holds. More precisely, let

$$\begin{aligned} ({\varvec{\upsilon }}_F(\mathbf{{u}}_h,p_h))^i :=|F|^{-1} \upsilon _{F,i}, \quad i = 1,\ldots ,d. \end{aligned}$$

(9.6)

Note that, consequently, (9.5) gives

$$\begin{aligned} \sum _{T \in \mathfrak T _V} m_{V,T,i} = \sum _{F \in \mathcal{F }_D} \upsilon _{F,i} \mathbf{{n}}_D {\cdot } \mathbf{{n}}_F = \sum _{F \in \mathcal{F }_D} \langle {\varvec{\upsilon }}_F(\mathbf{{u}}_h,p_h) \mathbf{{n}}_D {\cdot } \mathbf{{n}}_F, \mathbf e _i \rangle _F \end{aligned}$$

(9.7)

for every $D \in \mathcal D ^{\mathrm{int}}_h$ and the associated vertex $V$, $i = 1,\ldots ,d$. Let $F \in \partial \mathcal S _h^\mathrm{int}$ such that $F \subset \partial D$ for some $D \in \mathcal D _h$ and set

$$\begin{aligned} {\varvec{\varUpsilon }}_{F}(\mathbf{{u}}_h,p_h) :=(\nabla \mathbf{{u}}_h \mathbf{{n}}_F)|_F + {\varvec{\upsilon }}_F(\mathbf{{u}}_h,p_h). \end{aligned}$$

(9.8)

We then see that (9.4) together with (9.7) and (9.8) implies (7.30).