<figure><img src="image_1.png" /></figure>

HAL
open science

# Explicit T-coercivity for the Stokes problem: a coercive finite element discretization

Patrick Ciarlet, Erell Jamelot

<figure><img src="image_2.png" /></figure>

## To cite this version:

Patrick Ciarlet, Erell Jamelot. Explicit T-coercivity for the Stokes problem: a coercive finite element discretization. 2025. hal-04414789v6

HAL Id: hal-04414789

https://inria.hal.science/hal-04414789v6

Preprint submitted on 20 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Explicit *T*-coercivity for the Stokes problem: A coercive finite element discretization ☆

<figure><img src="image_5.png" /></figure>

Patrick Ciarlet Jr<sup>a, ID</sup>, Erell Jamelot<sup>b, ID,*</sup>

## ARTICLE INFO

MSC:

65N30
35J57
76D07

Keywords:

Stokes problem

$T$-coercivity

Stabilized $P^1 \times P^0$ pair

## ABSTRACT

Using the *T*-coercivity theory as advocated in Chesnel and Ciarlet (2013) [25], we propose a new variational formulation of the Stokes problem which does not involve nonlocal operators. With this new formulation, unstable finite element pairs are stabilized. In addition, the numerical scheme is easy to implement, and a better approximation of the velocity and the pressure is observed numerically when the viscosity is small.

## 1. Introduction

The Stokes problem describes the steady state of incompressible Newtonian flows. It is derived from the Navier-Stokes equations [33]. With regard to numerical analysis, the study of the Stokes problem helps to build an appropriate approximation of the Navier-Stokes equations. We propose here to write a new variational formulation of the Stokes problem using the *T*-coercivity theory, following Section 2 in [25].

In Section 2, we recall the *T*-coercivity theory as written in [26,25].

In Section 3 we apply this theory to the generalized Stokes problem, that is possibly with $div \mathbf{u} \neq 0$. In particular, solving this problem with homogeneous boundary condition helps solving the incompressible Stokes problem with inhomogeneous Dirichlet boundary conditions. We prove basic *T*-coercivity by finding an operator *T* such that the global variational formulation of the generalized Stokes problem is *T*-coercive.

From that point on, we assume that $div \mathbf{u} = 0$ and solve the incompressible Stokes problem. In Section 4, we build and analyse a new variational formulation with this operator *T* (explicit *T*-coercivity). Then in Section 5 we propose a numerical algorithm based on the new variational formulation and, in the Section after, we introduce discretizations and we study more specifically stability properties. As a particular case, it is shown that the unstable finite element pair $P^1 \times P^0$ yields stability. Finally, we provide some numerical experiments in Section 7 to illustrate our points.

Some concluding remarks are made in Section 8.

<footnote><sup>a</sup> POEMS, CNRS, Inria, ENSTA, Institut Polytechnique de Paris, 828 Boulevard des Maréchaux, 91120 Palaiseau, France</footnote>

<footnote><sup>b</sup> Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191, France</footnote>

<footnote>☆ Explicit *T*-coercivity for Stokes.</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: patrick.ciarlet@ensta-paris.fr (P. Ciarlet), erell.jamelot@cea.fr (E. Jamelot).</footnote>
## 2. T-coercivity

We recall here the T-coercivity theory as written in [26,25]. Consider the variational problem, where V and W are two Hilbert spaces and $\ell \in W'$:

$$ \text{Find } u \in V \text{ such that } \forall w \in W, a(u, w) = \ell(w). \quad (2.1) $$

Classically, we know that Problem (2.1) is well-posed if $a(\cdot, \cdot)$ satisfies the stability and the solvability conditions of the so-called Banach-Nečas-Babuška (BNB) Theorem (see e.g. [31, Theorem 25.9]). For some models, one can also prove the well-posedness using the T-coercivity theory. Whereas the BNB theorem relies on an abstract inf-sup condition, T-coercivity uses explicit inf-sup operators. We refer for instance to [18,47,14,25,15,13] for problems with sign-changing coefficients, to [23,22,52] for integral equations, to [16] for interior transmission problems, to [40,21,26,41,38] for Helmholtz-like problems, to [44,28,34] for the neutron diffusion equation, to [27] for the magnetostatic problem.

**Definition 1 (Basic T-coercivity).** Let V and W be two Hilbert spaces and $a(\cdot, \cdot)$ be a continuous bilinear form over $V \times W$. It is T-coercive if

$$ \exists T \in \mathcal{L}(V, W), \text{ bijective, } \exists \alpha > 0, \forall v \in V, |a(v, Tv)| \ge \alpha \|v\|_V^2. \quad (2.2) $$

Obviously, if (2.2) holds for some $T \in \mathcal{L}(V, W)$, then T is injective.

**Theorem 1 (Well-posedness [26,25]).** Let $a(\cdot, \cdot)$ be a continuous bilinear form on $V \times W$. The Problem (2.1) is well-posed if, and only if, the form $a(\cdot, \cdot)$ is T-coercive in the sense of Definition 1.

Two extensions of T-coercivity are proposed in [25, §2.3.2]. The first one is, when one is interested in solving the discretized problem, to mimic the design of the operator T at the discrete level to obtain a uniform discrete inf-sup condition: we call it the discrete T-coercivity. Among the above cited references, this is developed in [23,40,18,26,25,41,28,13,38,52,43]. The second one is the use of the operator T directly in the variational formulation: an equivalent variational formulation is obtained, which reads

$$ \text{Find } u \in V \text{ such that } \forall v \in V, a(u, Tv) = \ell(Tv). \quad (2.3) $$

By design, $a(\cdot, \cdot)$ is a continuous bilinear form, coercive on $V \times V$. We call it the explicit T-coercivity, see [27].

## 3. The Stokes problem

Let $\Omega$ be a connected bounded domain of $\mathbb{R}^d$, $d = 2, 3$, with a polygonal ($d = 2$) or Lipschitz polyhedral ($d = 3$) boundary $\partial\Omega$. We consider the generalized Stokes problem:

$$ \text{Find } (u, p) \text{ such that } \quad -\nu \Delta u + \text{grad } p = \mathbf{f}, \\ \text{div } u = g, \quad (3.1) $$

with Dirichlet boundary conditions for u and a normalization condition for p: $\int_{\Omega} p = 0$. If the Dirichlet boundary conditions are homogeneous, we write (3.1)$_H$, and (3.1)$_NH$ else.

The vector field u represents the velocity of the fluid and the scalar field p represents its pressure divided by the fluid density which is supposed to be constant. The first equation of (3.1) corresponds to the momentum balance equation and the second one corresponds to the conservation of the mass. The constant parameter $\nu > 0$ is the kinematic viscosity of the fluid.

Let us provide some definitions and reminders. Let us set $L^2(\Omega) := (L^2(\Omega))^d$, $H_0^1(\Omega) := (H_0^1(\Omega))^d$, $H^{-1}(\Omega) := (H^{-1}(\Omega))^d$ the dual space of $H_0^1(\Omega)$ and $L^2_{zmv}(\Omega) := \{q \in L^2(\Omega) | \int_{\partial\Omega} q = 0\}$, $L^2_{zmv}(\partial\Omega) := \{q \in L^2(\partial\Omega) | \int_{\partial\Omega} q = 0\}$.

The natural function space for the velocity is $H^1(\Omega)$ and, if homogeneous boundary conditions are prescribed, it is $H_0^1(\Omega)$, while for the pressure the natural function space is $L^2_{zmv}(\Omega)$. The data is $(f, g)$. The vector field $f \in H^{-1}(\Omega)$ represents a body force divided by the fluid density. The scalar field $g \in L^2_{zvm}(\Omega)$ is some abstract data that will be useful for the analysis of the case of nonhomogeneous Dirichlet boundary conditions.

Let us first recall the Poincaré-Steklov inequality:

$$ \exists C_{PS} > 0 |\forall v \in H_0^1(\Omega), \quad \|v\|_{L^2(\Omega)} \le C_{PS} \|\text{grad } v\|_{L^2(\Omega)}. \quad (3.2) $$

Thanks to this result, in $H_0^1(\Omega)$ the semi-norm is equivalent to the full norm, so we choose an inner product equal to $(v, w)_{H_0^1(\Omega)} = (\text{grad } v, \text{grad } w)_{L^2(\Omega)}$ with associated norm $\|v\|_{H_0^1(\Omega)} = \|\text{grad } v\|_{L^2(\Omega)}$. Let $v, w \in H_0^1(\Omega)$. We denote by $(v_i)_{i=1}^d$ (resp. $(w_i)_{i=1}^d$) the components of v (resp. w), and we set $\text{Grad } v = (\partial_j v_i)_{i,j=1}^d \in L^2(\Omega)$, where $L^2(\Omega) := [L^2(\Omega)]^{d \times d}$. We have:
$$ (\mathbf{Grad} \mathbf{v}, \mathbf{Grad} \mathbf{w})_{L^2(\Omega)} = (\mathbf{v}, \mathbf{w})_{H_0^1(\Omega)} = \sum_{i=1}^{d} (v_i, w_i)_{H_0^1(\Omega)} $$

and:

$$ \|\mathbf{v}\|_{H_0^1(\Omega)} = \left( \sum_{j=1}^{d} \|v_j\|_{H_0^1(\Omega)}^2 \right)^{1/2} = \|\mathbf{Grad} \mathbf{v}\|_{L^2(\Omega)} $$

Let us set $\mathbf{V} := \{\mathbf{v} \in H_0^1(\Omega) \mid \text{div } \mathbf{v} = 0\}$. The vector space $\mathbf{V}$ is a closed subset of $H_0^1(\Omega)$. We denote by $\mathbf{V}^{\perp}$ the orthogonal complement of $\mathbf{V}$ in $H_0^1(\Omega)$. We recall that there exists a (nonlocal) right inverse of the divergence operator:

**Proposition 1** ([3, Theorem 3.1 (e)]) . The operator div : $H_0^1(\Omega) \to L^2(\Omega)$ is an isomorphism of $\mathbf{V}^{\perp}$ onto $L_{zmv}^2(\Omega)$. We call $C_{\text{div}} \ge 1$ the smallest constant such that:

$$ \forall p \in L_{zmv}^2(\Omega), \exists! \tilde{v}_p \in \mathbf{V}^{\perp} \mid \text{div } \tilde{v}_p = p \text{ and } \|\tilde{v}_p\|_{H_0^1(\Omega)} \le C_{\text{div}} \|p\|_{L^2(\Omega)}. \quad (3.3) $$

In the statement of Proposition 1, we note that since

$$ \forall \mathbf{v} \in H_0^1(\Omega), \quad \|\mathbf{v}\|_{H_0^1(\Omega)}^2 = \|\mathbf{curl} \mathbf{v}\|_{L^2(\Omega)}^2 + \|\text{div } \mathbf{v}\|_{L^2(\Omega)}^2, $$

one has necessarily $C_{\text{div}} \ge 1$.

The variational formulation of Problem (3.1)$_H$ reads:

Find $(\mathbf{u}, p) \in H_0^1(\Omega) \times L_{zmv}^2(\Omega)$ such that

$$ \begin{aligned} \mathbf{v}(\mathbf{u}, \mathbf{v})_{H_0^1(\Omega)} - (p, \text{div } \mathbf{v})_{L^2(\Omega)} &= \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} & \forall \mathbf{v} \in H_0^1(\Omega); \\ (q, \text{div } \mathbf{u})_{L^2(\Omega)} &= (g, q)_{L^2(\Omega)} & \forall q \in L_{zmv}^2(\Omega). \end{aligned} \quad (3.4) $$

This is a saddle point-problem. Using classical theory, one proves that Problem (3.4) is well-posed with the help of Poincaré-Steklov inequality (3.2) and Proposition 1. Check for instance the proof of [33, Theorem 1.5.1].

Let us set $\mathcal{X} := H_0^1(\Omega) \times L_{zmv}^2(\Omega)$ which is a Hilbert space which we endow with the following norm:

$$ \| (\mathbf{v}, q) \|_{\mathcal{X}, \nu} = \left( \| \mathbf{v} \|_{H_0^1(\Omega)}^2 + \nu^{-2} \| q \|_{L^2(\Omega)}^2 \right)^{1/2}. \quad (3.5) $$

This norm is chosen to account for physical phenomena due to small viscosity. Typically, if the result of a physics experiment is $(\mathbf{u}, p)$, the ratio of the two components of the norm $\| (\mathbf{u}, p) \|_{\mathcal{X}, \nu}$, respectively equal to $\| \mathbf{u} \|_{H_0^1(\Omega)}$ and $\nu^{-1} \| p \|_{L^2(\Omega)}$, varies linearly with $\mathbf{v}$.

We define the following continuous symmetric bilinear form:

$$ \left\{ \begin{array}{l} a: \mathcal{X} \times \mathcal{X} \to \mathbb{R} \\ (\mathbf{u}', p') \times (\mathbf{v}, q) \mapsto \mathbf{v}(\mathbf{u}', \mathbf{v})_{H_0^1(\Omega)} - (p', \text{div } \mathbf{v})_{L^2(\Omega)} - (q, \text{div } \mathbf{u}')_{L^2(\Omega)}. \end{array} \right. \quad (3.6) $$

We also define the linear continuous form:

$$ \ell: \mathcal{X} \to \mathbb{R} \\ (\mathbf{v}, q) \mapsto \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} - (g, q)_{L^2(\Omega)}. \quad (3.7) $$

We can write Problem (3.1)$_H$ in an equivalent way as follows:

$$ \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ s.t. } a((\mathbf{u}, p), (\mathbf{v}, q)) = \ell((\mathbf{v}, q)) \quad \forall (\mathbf{v}, q) \in \mathcal{X}. \quad (3.8) $$

Let us prove that Problem (3.8) is well-posed using basic T-coercivity.

**Proposition 2.** The bilinear form $a(\cdot, \cdot)$ is T-coercive.

**Proof.** We follow here the proof given in [6,17]. Let us consider $(\mathbf{u}', p') \in \mathcal{X}$ and let us build $(\mathbf{v}^*, q^*) = T(\mathbf{u}', p') \in \mathcal{X}$ satisfying (2.2) (with $V = W = \mathcal{X}$). We need three main steps.

1. According to Proposition 1, there exists $\tilde{v}_{p'} \in \mathbf{V}^{\perp}$ such that:

$$ \text{div } \tilde{v}_{p'} = p' \text{ in } \Omega \text{ and } \| \tilde{v}_{p'} \|_{H_0^1(\Omega)} \le C_{\text{div}} \| p' \|_{L^2(\Omega)}. \quad (3.9) $$

Let us set $(\mathbf{v}^*, q^*) := (\lambda \mathbf{u}' - \nu^{-1} \tilde{v}_{p'}, -\lambda p')$, with $\lambda > 0$ to be fixed. We obtain:

$$ a((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) = \nu \lambda \|\mathbf{u}'\|_{H_0^1(\Omega)}^2 + \nu^{-1} \|p'\|_{L^2(\Omega)}^2 - (\mathbf{u}', \tilde{v}_{p'})_{H_0^1(\Omega)}. \quad (3.10) $$
2. In order to bound the last term of (3.10), we use the Young inequality and then inequality (3.9), and it follows that for all $\eta > 0$:

$$ (\mathbf{u}', \bar{\mathbf{v}}_{p'})_{\mathbf{H}_0^1(\Omega)} \le \frac{\eta}{2} \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + \frac{\eta^{-1}}{2} (C_{\text{div}})^2 \|p'\|_{L^2(\Omega)}^2. \quad (3.11) $$

3. Using the bound (3.11) in (3.10):

$$ a((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge (\nu\lambda - \frac{\eta}{2}) \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + (\nu^{-1} - \frac{\eta^{-1}}{2} (C_{\text{div}})^2) \|p'\|_{L^2(\Omega)}^2. $$

We look for $\eta > 0$ such that $2\nu\lambda > \eta$ and $\eta > \frac{\nu}{2} (C_{\text{div}})^2$, which amounts to requiring

$$ \lambda > \frac{1}{4} (C_{\text{div}})^2. $$

According to the above, provided that $\lambda > \frac{1}{4} (C_{\text{div}})^2$, we have proved that the operator $T_\lambda \in \mathcal{L}(\mathcal{X})$ defined by $T_\lambda((\mathbf{u}', p')) = (\lambda\mathbf{u}' - \nu^{-1}\bar{\mathbf{v}}_{p'}, -\lambda p')$ is such that:

$$ \exists \alpha_\lambda > 0, \forall (\mathbf{u}', p') \in \mathcal{X}, a((\mathbf{u}', p'), T_\lambda((\mathbf{u}', p'))) \ge \alpha_\lambda \| (\mathbf{u}', p') \|_{\mathcal{X},v}^2. $$

As noted in Section 2, the injectivity of the operator $T_\lambda$ follows. Given $(\mathbf{v}^*, q^*) \in \mathcal{X}$, choosing $(\mathbf{u}', p') = ((\lambda^{-1}\mathbf{v}^* - \nu^{-1}\lambda^{-2}\bar{\mathbf{v}}_{q^*}, -\lambda^{-1}q^*))$ yields $T_\lambda((\mathbf{u}', p')) = (\mathbf{v}^*, q^*)$. Hence, the operator $T_\lambda \in \mathcal{L}(\mathcal{X})$ is bijective. $\square$

**Remark 1.** In the above proof, we established that the bijective operator $T_\lambda((\mathbf{u}', p')) = (\lambda\mathbf{u}' - \nu^{-1}\bar{\mathbf{v}}_{p'}, -\lambda p')$ leads to T-coercivity as soon as $\lambda > \frac{1}{4} (C_{\text{div}})^2$. Observe that one can have even more flexibility in the choice of T by choosing a different factor in front of $\bar{\mathbf{v}}_{p'}$, and then choosing $\lambda$ accordingly.

**Remark 2.** Interestingly, another operator T has been proposed in [4], when $\Omega$ is a 2D convex domain. The main difference is that another right inverse of the divergence operator $p' \mapsto \bar{\mathbf{v}}_{p'}$ is chosen. Unfortunately, some orthogonality properties are lost with this choice, so difficulties arise. We refer to Appendix A for a discussion.

We can now prove the following result for the generalized Stokes problem:

**Theorem 2.** Problem (3.8) is well-posed, so it admits one and only one solution for any $(\mathbf{f}, g) \in \mathbf{H}^{-1}(\Omega) \times L_{zmv}^2(\Omega)$. Writing $\mathbf{u} = \mathbf{u}_0 + \mathbf{u}_\perp$ with $\mathbf{u}_0 \in \mathbf{V}$ and $\mathbf{u}_\perp \in \mathbf{V}^\perp$, the solution is such that:

$$ \forall \mathbf{f} \in \mathbf{H}^{-1}(\Omega), \forall g \in L_{zmv}^2(\Omega) \quad \left\{ \begin{array}{l} \|\mathbf{u}_\perp\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|g\|_{L^2(\Omega)}, \\ \|\mathbf{u}_0\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}, \\ \|p\|_{L^2(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} + \nu C_{\text{div}}^2 \|g\|_{L^2(\Omega)}. \end{array} \right. \quad (3.12) $$

**Proof.** According to Proposition 2, the continuous bilinear form $a(\cdot, \cdot)$ is T-coercive. Hence, according to Theorem 1, Problem (3.8) is well-posed. Let us now derive (3.12). Consider $(\mathbf{u}, p)$ the unique solution of Problem (3.8), where $\mathbf{u} = \mathbf{u}_0 + \mathbf{u}_\perp$ with $\mathbf{u}_0 \in \mathbf{V}$ and $\mathbf{u}_\perp \in \mathbf{V}^\perp$. Choosing $\mathbf{v} = 0$, we obtain that $\text{div} \mathbf{u}_\perp = g$, so that $\|\mathbf{u}_\perp\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|g\|_{L^2(\Omega)}$. Now, choosing $\mathbf{v} = \mathbf{u}_0$ and using orthogonality, we have: $\nu \|\mathbf{u}_0\|_{\mathbf{H}_0^1(\Omega)}^2 = \langle \mathbf{f}, \mathbf{u}_0 \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} \le \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|\mathbf{u}_0\|_{\mathbf{H}_0^1(\Omega)}$, so that: $\|\mathbf{u}_0\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}$. Next, we choose, in (3.8), $q = 0$ and $\mathbf{v} = -\bar{\mathbf{v}}_p \in \mathbf{V}^\perp$, where $\text{div} \bar{\mathbf{v}}_p = p$ and $\|\bar{\mathbf{v}}_p\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|p\|_{L^2(\Omega)}$ (see Proposition 1). Since $\mathbf{u}_0 \in \mathbf{V}$, it holds that $(\mathbf{u}_0, \bar{\mathbf{v}}_p)_{\mathbf{H}_0^1(\Omega)} = 0$. This gives:

$$ \begin{aligned} \|p\|_{L^2(\Omega)}^2 &= (p, \text{div} \bar{\mathbf{v}}_p)_{L^2(\Omega)} = -\langle \mathbf{f}, \bar{\mathbf{v}}_p \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} + \nu \langle \Delta \mathbf{u}_\perp, \bar{\mathbf{v}}_p \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)}, \\ &\le \left( \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} + \nu \|\mathbf{u}_\perp\|_{\mathbf{H}_0^1(\Omega)} \right) \|\bar{\mathbf{v}}_p\|_{\mathbf{H}_0^1(\Omega)} \\ &\le C_{\text{div}} \left( \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} + \nu C_{\text{div}} \|g\|_{L^2(\Omega)} \right) \|p\|_{L^2(\Omega)}, \end{aligned} $$

so that: $\|p\|_{L^2(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} + \nu C_{\text{div}}^2 \|g\|_{L^2(\Omega)}$. $\square$

In the following sections, we will consider that $g = 0$, i.e. we focus on the classical incompressible Stokes model.

## 4. New variational formulations

We solve the classical incompressible Stokes model, with $g = 0$.
4.1. Explicit T-coercivity

Let $\lambda > \frac{1}{4}(C_{\text{div}})^2$. In Remark 1, we introduced the operator $T_{\lambda} \in \mathcal{L}(\mathcal{X})$ defined by:

$$ T_{\lambda}: \begin{cases} \mathcal{X} \to \mathbb{R} \\ (\mathbf{v}, q) \mapsto (\lambda \mathbf{v} - \nu^{-1} \tilde{\mathbf{v}}_q, -\lambda q) \end{cases}, \qquad (4.1) $$

where $\tilde{\mathbf{v}}_q \in \mathbf{V}^{\perp}$ is given by (3.3): $\text{div} \tilde{\mathbf{v}}_q = q$ and $\|\tilde{\mathbf{v}}_q\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|q\|_{L^2(\Omega)}$. We now write the variational formulation (3.8) with test function $T_{\lambda}((\mathbf{v}, q))$ instead of $(\mathbf{v}, q)$. Let us define $\tilde{a}_{\lambda}((\mathbf{u}', p'), (\mathbf{v}, q)) = a((\mathbf{u}', p'), T(\mathbf{v}, q))$. We have:

$$ \begin{aligned} \tilde{a}_{\lambda}((\mathbf{u}', p'), (\mathbf{v}, q)) ={}& \nu \lambda (\mathbf{u}', \mathbf{v})_{\mathbf{H}_0^1(\Omega)} - (\mathbf{u}', \tilde{\mathbf{v}}_q)_{\mathbf{H}_0^1(\Omega)} \\ & - \lambda(p', \text{div } \mathbf{v})_{L^2(\Omega)} + \nu^{-1}(p', q)_{L^2(\Omega)} \\ & + \lambda(q, \text{div } \mathbf{u}')_{L^2(\Omega)}. \end{aligned} \qquad (4.2) $$

The term with $\tilde{\mathbf{v}}_q$ requires some knowledge of the nonlocal right inverse of the divergence operator. We will see in the next section that it can be removed.

According to Proposition 2, we have the...

**Proposition 3.** The bilinear form (4.2) is coercive.

Introducing $\ell_{\lambda}((\mathbf{v}, q)) := \ell(T_{\lambda}(\mathbf{v}, q))$, we can propose a first new variational formulation to Problem (3.1)$_H$, which reads:

$$ \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ s.t. } \tilde{a}_{\lambda}((\mathbf{u}, p), (\mathbf{v}, q)) = \ell_{\lambda}((\mathbf{v}, q)) \quad \forall (\mathbf{v}, q) \in \mathcal{X}. \qquad (4.3) $$

Using the bijectivity of $T_{\lambda}$, it is obvious that (4.3) is equivalent to (3.8), so well-posedness of (4.3) is a direct consequence of the well-posedness of (3.8), and vice versa.

Let us derive an “explicit” expression of $\ell_{\lambda}((\mathbf{v}, q))$, which will prove useful later on. We recall that we have in the sense of distributions:

$$ -\Delta(\cdot) = \text{curl curl}(\cdot) - \text{grad div}(\cdot), \qquad (4.4) $$

so that

$$ \forall \mathbf{f}' \in \mathbf{H}^{-1}(\Omega), \quad \exists! (z_{\mathbf{f}'}, \mathbf{w}_{\mathbf{f}'}) \in L^2_{zmv}(\Omega) \times \mathbf{V} \mid \mathbf{f}' = \text{grad } z_{\mathbf{f}'} - \Delta \mathbf{w}_{\mathbf{f}'} = \text{grad } z_{\mathbf{f}'} + \text{curl curl } \mathbf{w}_{\mathbf{f}'}, \qquad (4.5) $$

where $(\mathbf{w}_{\mathbf{f}'}, z_{\mathbf{f}'})$ solves (3.1)$_H$ with $\nu = 1$ and data $\mathbf{f} = \mathbf{f}'$ and $g = 0$.

**Proposition 4.** Let $\mathbf{f}' \in \mathbf{H}^{-1}(\Omega)$. Given, $q \in L^2_{zmv}$, let $\tilde{\mathbf{v}}_q \in \mathbf{V}^{\perp}$ be defined by (3.3). We have:

$$ \langle \mathbf{f}', \tilde{\mathbf{v}}_q \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} = -(z_{\mathbf{f}'}, q)_{L^2(\Omega)}. \qquad (4.6) $$

**Proof.** Let $\mathbf{f}'$ be decomposed as in (4.5).

On the one hand, integrating by parts twice and using (4.4), we get:

$$ \begin{aligned} -\langle \text{curl curl } \mathbf{w}_{\mathbf{f}'}, \tilde{\mathbf{v}}_q \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} &= -\langle \text{curl curl } \tilde{\mathbf{v}}_q, \mathbf{w}_{\mathbf{f}'} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)}, \\ &= \langle \Delta \tilde{\mathbf{v}}_q, \mathbf{w}_{\mathbf{f}'} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} - \langle \text{grad } q, \mathbf{w}_{\mathbf{f}'} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)}, \\ &= -(\tilde{\mathbf{v}}_q, \mathbf{w}_{\mathbf{f}'} \mathbf{H}_0^1(\Omega)) + (q, \text{div } \mathbf{w}_{\mathbf{f}'})_{L^2(\Omega)}, \\ &= 0, \end{aligned} $$

resp. since $\text{div } \tilde{\mathbf{v}}_q = q$, $\tilde{\mathbf{v}}_q \in \mathbf{V}^{\perp}$, $\mathbf{w}_{\mathbf{f}'} \in \mathbf{V}$, and $\text{div } \mathbf{w}_{\mathbf{f}'} = 0$.

On the other hand, integrating by parts once, we have:

$$ \langle \text{grad } z_{\mathbf{f}'}, \tilde{\mathbf{v}}_q \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} = -(z_{\mathbf{f}'}, q)_{L^2(\Omega)}, $$

so the claim follows. $\square$

So, we have the
**Corollary 1.** The right-hand-side $\ell_\lambda$ is equal to

$$ (v, q) \mapsto \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} + v^{-1}(z_f, q)_{L^2(\Omega)}. $$

**Proof.** This is an obvious consequence of (4.6). $\square$

## 4.2. New variational formulation using orthogonality

Going back to the original Problem (3.1)$_H$, we note that all solutions $(\mathbf{u}', p')$ to Problem (4.3) are such that $\mathbf{u}' \in \mathbf{V}$. So, in the statement of Problem (4.3), the term $(\mathbf{u}', \tilde{\mathbf{v}}_q)_{\mathbf{H}_0^1(\Omega)}$ can be removed by orthogonality from the expression (4.2) of the bilinear form $\tilde{a}_\lambda$: as a matter of fact, according to (3.3), for all $q \in L^2_{mv}(\Omega)$, $(\mathbf{u}', \tilde{\mathbf{v}}_q)_{\mathbf{H}_0^1(\Omega)} = 0$ because $\tilde{\mathbf{v}}_q \in \mathbf{V}^\perp$. Interestingly, in this manner one improves the stability in the proof of Proposition 2, because the cross term vanishes in the expression (3.10): no treatment is required and the resulting bilinear form is coercive for all $\lambda > 0$. So, we propose to remove the second term in the expression (4.2) of the bilinear form with operator $T_\lambda$, that is, we introduce the bilinear form $a_\lambda(\cdot, \cdot)$ such that:

$$ \left\{ \begin{array}{l} a_\lambda : \quad \mathcal{X} \times \mathcal{X} \to \mathbb{R} \\ ((\mathbf{u}', p'), (\mathbf{v}, q)) \mapsto v \lambda (\mathbf{u}', \mathbf{v})_{\mathbf{H}_0^1(\Omega)} + v^{-1} (p', q)_{L^2(\Omega)} \\ \phantom{((\mathbf{u}', p'), (\mathbf{v}, q)) \mapsto} + \lambda [(q, \operatorname{div} \mathbf{u}')_{L^2(\Omega)} - (p', \operatorname{div} \mathbf{v})_{L^2(\Omega)}] \end{array} \right. \qquad (4.7) $$

**Proposition 5.** The bilinear form (4.7) is coercive.

**Proof.** We have:

$$ \begin{aligned} a_\lambda ((\mathbf{u}', p'), (\mathbf{u}', p')) &= v \lambda \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + v^{-1} \|p'\|_{L^2(\Omega)}^2, \\ &\geq v \min(1, \lambda) \|(\mathbf{u}', p')\|_{\mathcal{X}, \mathbf{v}}^2. \quad \square \end{aligned} $$

With the help of explicit $T$-coercivity and using orthogonality, we can now propose a second new variational formulation to Problem (3.1)$_H$, which reads:

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ such that} \\ a_\lambda((\mathbf{u}, p), (\mathbf{v}, q)) = \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)} + v^{-1} (z_f, q)_{L^2(\Omega)} \quad \forall (\mathbf{v}, q) \in \mathcal{X}. \end{array} \right. \qquad (4.8) $$

**Theorem 3.** For all $\lambda > 0$, Problem (4.8) is well-posed and is equivalent to Problem (3.1)$_H$ with $g = 0$.

By contrast with Remark 1, the result here holds for all $\lambda > 0$.

**Proof.** The bilinear form $a_\lambda(\cdot, \cdot)$ is continuous and coercive. Let $\mathbf{f} \in \mathbf{H}^{-1}(\Omega)$, and let $z_f$ be given by (4.5), the linear form $\ell_\lambda(\cdot)$ is continuous. According to Lax-Milgram Theorem, Problem (4.8) is well-posed. It exists a unique solution $(\mathbf{u}, p)$ which depends continuously on the data.

Regarding the equivalence with Problem (3.1)$_H$ with $g = 0$, we already observed that solving (4.3) is equivalent to solving (3.8), and that both of them are equivalent to solving Problem (3.1)$_H$ with $g = 0$. Then, if $(\mathbf{u}, p)$ solves Problem (3.1)$_H$ with data $\mathbf{f} \in \mathbf{H}^{-1}(\Omega)$ and $g = 0$, one has in particular that $\mathbf{u} \in \mathbf{V}$. Hence it follows from the above that $(\mathbf{u}, p)$ solves (4.8) with data $\mathbf{f}$ and $g = 0$, resp. $z_f$ given by (4.5).

Conversely, assume that $(\mathbf{u}, p)$ solves (4.8), with data $\mathbf{f} \in \mathbf{H}^{-1}(\Omega)$, and $z_f$ given by (4.5). Denote by $(\mathbf{u}^\dagger, p^\dagger)$ the solution to Problem (3.1)$_H$ with data $\mathbf{f}$. According to the above and by linearity, $(\mathbf{u} - \mathbf{u}^\dagger, p - p^\dagger)$ solves (4.8) with vanishing data, hence to is equal to $(0, 0)$ by uniqueness: in other words, $(\mathbf{u}, p)$ solves Problem (3.1)$_H$ with data $\mathbf{f}$. $\square$

Notice that we can write Problem (4.8) with two equations as:

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ s.t. for all } (\mathbf{v}, q) \in \mathcal{X} \\ \text{(i)} \quad v \lambda (\mathbf{u}, \mathbf{v})_{\mathbf{H}_0^1(\Omega)} - \lambda(p, \operatorname{div} \mathbf{v})_{L^2(\Omega)} = \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)}, \\ \text{(ii)} \quad \lambda(q, \operatorname{div} \mathbf{u})_{L^2(\Omega)} + v^{-1}(p, q)_{L^2(\Omega)} = v^{-1}(z_f, q)_{L^2(\Omega)}. \end{array} \right. \qquad (4.9) $$

This new variational formulation appears as a stabilized variational formulation, in the sense of §II.1.2 in [33], pages 120-123. It can be solved once $z_f$, or a suitable approximation of $z_f$, is available. This suggests to use (4.9) as a post-processing step as follows:

* Compute $z_f$ by solving numerically the Stokes problem (3.4) with $v = 1$, and data $\mathbf{f}$;
• Solve (4.9) with the data **f** and the computed z<sub>f</sub>.

Notice that we can recover that the solution **u** to (4.9) belongs to V in a simple manner: let us split **u** = **u**<sub>0</sub> + **u**<sub>⊥</sub>, where (**u**<sub>0</sub>, **u**<sub>⊥</sub>) ∈ V×V<sub>⊥</sub>, so that div **u** = div **u**<sub>⊥</sub>. Choosing **v** = **u**<sub>⊥</sub> in (4.9)-(i), and **q** = div **u** in (4.9)-(ii):

$$
\left\{
\begin{array}{ll}
\text{(i)} & v \|\mathbf{u}_{\perp}\|_{H_0^1(\Omega)}^2 - (p, \operatorname{div} \mathbf{u})_{L^2(\Omega)} = \langle \mathbf{f}, \mathbf{u}_{\perp} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = -(z_f, \operatorname{div} \mathbf{u})_{L^2(\Omega)}, \\
\text{(ii)} & \lambda \|\operatorname{div} \mathbf{u}\|_{L^2(\Omega)}^2 + v^{-1}(p, \operatorname{div} \mathbf{u})_{L^2(\Omega)} = v^{-1}(z_f, \operatorname{div} \mathbf{u})_{L^2(\Omega)}.
\end{array}
\right.
\tag{4.10}
$$

Summing (4.10)-(i) and (4.10)-(ii) times v, we obtain:

$$
v \left( \| \mathbf{u}_{\perp} \|_{H_0^1(\Omega)}^2 + \lambda \| \operatorname{div} \mathbf{u} \|_{L^2(\Omega)}^2 \right) = 0.
$$

Hence, **u**<sub>⊥</sub> = 0 and div **u** = 0.

In Appendix A, we present another variational formulation (not relying on orthogonality) that can be obtained with the help of explicit T-coercivity. In Appendix B, we address the case of the classical incompressible Stokes model, with nonhomogeneous Dirichlet boundary conditions.

5. Numerical algorithms

We solve the classical incompressible Stokes model, with g = 0. We propose below two strategies depending on whether z<sub>f</sub>
is known or not. Below, we consider homogeneous Dirichlet boundary conditions (see Appendix B for nonhomogeneous Dirichlet
boundary conditions).

If z<sub>f</sub> is known, one solves (4.9) directly. If this is not the case, we propose to start from some approximate value of z<sub>f</sub>, and then
to iterate by solving a series of problems like (4.9).

We propose here to study the numerical algorithms derived from these strategies. The discretization will be detailed in Section 7.
Let $N_u$ (resp. $N_p$) be the number of discrete velocity (resp. pressure) unknowns. Let $\mathcal{A} \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_u}$ be the velocity stiffness matrix,
$\mathcal{B} \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_u}$ be the velocity-pressure coupling matrix and $\mathcal{M} \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_p}$ be the pressure mass matrix. So, after discretization, if
$z_f$ is known, one solves a linear system like

$$
\left\{
\begin{array}{l}
\text{Find } (\underline{U}, \underline{P}) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\
\nu \lambda \mathcal{A} \underline{U} - \lambda \mathcal{B}^T \underline{P} = \lambda \underline{F} \\
\lambda \mathcal{B} \underline{U} + \nu^{-1} \mathcal{M} \underline{P} = \nu^{-1} \mathcal{M} \underline{Z}
\end{array}
\right.
\tag{5.1}
$$

where $\underline{U} \in \mathbb{R}^{N_u}$ represents the discrete velocity, resp. $\underline{P} \in \mathbb{R}^{N_p}$ the discrete pressure, while $\underline{F} \in \mathbb{R}^{N_u}$ stands for $\mathbf{f}$, and $\underline{Z} \in \mathbb{R}^{N_p}$ stands for $z_f$. Classically, $\underline{U} \mapsto (\underline{A}\underline{U}|\underline{U})$ measures discrete velocities, $\underline{P} \mapsto (\underline{M}\underline{P}|\underline{P})$ measures discrete pressures and, as a consequence, $\underline{U} \mapsto (\underline{B}\underline{U}|\underline{M}^{-1}\underline{B}\underline{U})$ measures the divergence of the discrete velocities.

On the other hand, if $z_f$ is not known, starting from an initial guess $\underline{P}^{-1} \in \mathbb{R}^{N_p}$, for $n=0,1,\dots$, one solves linear systems like

$$
\left\{
\begin{array}{l}
\text{Find } (\underline{U}^n, \underline{P}^n) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\
\nu \lambda \mathcal{A} \underline{U}^n - \lambda \mathcal{B}^T \underline{P}^n = \lambda \underline{F}^n \\
\lambda \mathcal{B} \underline{U}^n + \nu^{-1} \mathcal{M} \underline{P}^n = \nu^{-1} \mathcal{M} \underline{P}^{n-1}
\end{array}
\right.
\tag{5.2}
$$

In this situation, one needs to prove that the iterative solver is actually converging, and to provide a stopping criterion.

**Remark 3.** The matrix of algorithm (5.2) is similar to the one of the first-order artificial compressibility algorithm [50,35,36]. In [35] and [36], the space discretization is done with the MAC approximation on a Cartesian grid [39]. In [36], it is also done with the Taylor-Hood finite element [51].

**Theorem 4 (Convergence).** Assume that the matrices $\mathcal{A} \in \mathbb{R}^{N_u} \times N_u$ and $\mathcal{M} \in \mathbb{R}^{N_p} \times N_p$ are symmetric positive-definite. Then, the sequence $(\underline{U}^n, \underline{P}^n)_n$ of solutions to (5.2) is converging, and it converges to $(\underline{U}^\infty, \underline{P}^\infty) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p}$ which is governed by

$$
\left\{
\begin{array}{l}
\text{Find } (\underline{U}^{\infty}, \underline{P}^{\infty}) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\
\nu \mathcal{A} \underline{U}^{\infty} - \mathcal{B}^T \underline{P}^{\infty} = \underline{F}^{\infty} \\
\mathcal{B} \underline{U}^{\infty} = 0 \\
\mathcal{M} (\underline{P}^{\infty} - \underline{P}^{n-1}) \in \operatorname{im}(\mathcal{B})
\end{array}
\right.
\tag{5.3}
$$

**Proof.** In the proof, we use the change of variable $\underline{Y} := \mathcal{M}^{\frac{1}{2}} \underline{P} \in \mathbb{R}^{N_p}$, and specifically $\underline{Y}^n = \mathcal{M}^{\frac{1}{2}} \underline{P}^n$ for $n \ge -1$. Then let $\underline{\delta} \underline{P}^n = \underline{P}^n - \underline{P}^{n-1}$ and $\underline{\delta} \underline{Y}^n := (\mathcal{M}^{\frac{1}{2}} \underline{\delta} \underline{P}^n)$ for $n \ge 0$.
1/ We notice first that the following recursive relation holds:

$$ \forall n \ge 1, \quad (M + \lambda BA^{-1}B^T)\delta \underline{P}^n = M\delta \underline{P}^{n-1}, \qquad (5.4) $$

where the matrix $\lambda BA^{-1}B^T$ is a symmetric positive matrix. By construction, $\delta \underline{Y}^n = \nu \lambda M^{-\frac{1}{2}} \underline{B} \underline{U}^n \in \text{im}(M^{-\frac{1}{2}}\underline{B})$. Since $\mathbb{R}^{N_p} = \text{im}(M^{-\frac{1}{2}}\underline{B}) \oplus \ker(B^T M^{-\frac{1}{2}})$, one has $\delta \underline{Y}^n = \ker(B^T M^{-\frac{1}{2}})^\perp$ and (5.4) reads:

$$ \forall n \ge 1, \quad (I + \lambda M^{-\frac{1}{2}}BA^{-1}B^T M^{-\frac{1}{2}})\delta \underline{Y}^n = \delta \underline{Y}^{n-1}. \qquad (5.5) $$

The iteration matrix $(I + \lambda M^{-\frac{1}{2}}BA^{-1}B^T M^{-\frac{1}{2}})^{-1}$ restricted to $\ker(B^T M^{-\frac{1}{2}})^\perp$ is a symmetric positive definite matrix, whose eigenvalues are all strictly smaller than 1. Hence, its spectral radius $\rho_\lambda$ is strictly smaller than 1.

2/ We deduce that it exists $C_\lambda \in (0, 1)$ independent of the data and of $n$ such that for all $n \ge 1$, $|\delta \underline{Y}^n|_2 \le C_\lambda |\delta \underline{Y}^{n-1}|_2$. Thus, the series $\sum_{n \ge 1} \delta \underline{Y}^n$ is a convergent series. Noticing that $\sum_{n=0}^N \delta \underline{Y}^n = \underline{Y}^N - \underline{Y}^{-1}$, we conclude that the sequence $(\underline{Y}^n)_n$ is a convergent sequence and we set $\underline{Y}^\infty = \lim_{n \to \infty} \underline{Y}^n$. By construction $\underline{Y}^\infty - \underline{Y}^{-1} \in \text{im}(M^{-\frac{1}{2}}\underline{B}) = \ker(B^T M^{-\frac{1}{2}})^\perp$, and for all $n$, the algorithm doesn't change the orthogonal projection of $\underline{Y}^n$ onto $\ker(B^T M^{-\frac{1}{2}})$, which is equal to that of $\underline{Y}^{-1}$ onto $\ker(B^T M^{-\frac{1}{2}})$.

3/ Because $M^{-\frac{1}{2}}$ is invertible, we obtain that the sequence $(\underline{P}^n)_n$ is also a convergent sequence and we call $\underline{P}^\infty = M^{-\frac{1}{2}}\underline{Y}^\infty$ its limit. In the same way, noting that $\underline{A}$ is invertible, we obtain that the sequence $(\underline{U}^n)_n$ is a convergent sequence and we call $\underline{U}^\infty = \nu^{-1}\underline{A}^{-1}(\underline{F}_u + \underline{B}^T\underline{P}^\infty)$ its limit.

Passing to the limit, it holds that $\underline{B}\underline{U}^\infty = 0$ and $M(\underline{P}^\infty - \underline{P}^{-1}) \in \text{im}(\underline{B})$. $\square$

**Remark 4.** In the above proof, we notice that the mapping $\lambda \mapsto \rho_\lambda$ is monotonically decreasing. In section 7 dedicated to the numerical experiments, we set the value of $\lambda$ equal to 10.

**Remark 5.** There are two exclusive cases for Problem (5.3):

* either, $\text{im}(\underline{B}) = \mathbb{R}^{N_p}$: a discrete inf-sup condition holds. Uniqueness of $\underline{P}^\infty$ is guaranteed and the last line of (5.3) is trivial;

* or, $\text{im}(\underline{B}) \subset \mathbb{R}^{N_p}$: there is no discrete inf-sup condition. However in this case, the last line of (5.3) guarantees the uniqueness of $\underline{P}^\infty$. As a matter of fact, the algorithm does not change the orthogonal projection of $\underline{P}^n$ onto $\ker(B^T)$, which is equal to that of $\underline{P}^{-1}$ onto $\ker(B^T)$, where orthogonality is understood with respect to the inner product $(M \cdot |\cdot|)$.

Let us consider that the assumptions of Theorem 4 hold. A critical ingredient is to define the stopping criterion, in particular what are the relevant quantities of interest. By definition, $(M\delta \underline{P}^n|\delta \underline{P}^n) = |\delta \underline{Y}^n|_2^2$ is strictly monotonically decreasing. Hence, denoting by $|\cdot|_A: \underline{U} \mapsto (\underline{A}\underline{U}|\underline{U})^{1/2}$ the norm of the discrete velocity, resp. $|\cdot|_M: \underline{P} \mapsto (\underline{M}\underline{P}|\underline{P})^{1/2}$ the norm of the discrete pressure, we define the stopping criterion by comparing the adimensionalized quantities $|\delta \underline{P}^n|_M$ and $|\underline{U}^n|_A$ (we note that both quantities are easily computable). This amounts to setting the stopping criterion as

$$ |\delta \underline{P}^n|_M \le \epsilon |\underline{U}^n|_A \qquad (5.6) $$

for ad hoc $\epsilon > 0$. In our numerical experiments, we fixed $n = 1$ or $n = 8$ with $\epsilon = 10^{-12}$.

Finally, we recall that $\underline{B}\underline{U}^n$ stands for the divergence of the discrete velocity div $\underline{u}_h^n$, and that one aims at diminishing the value of its norm. Considering again the auxiliary variable $\delta \underline{Y}^n$ introduced in the proof of Theorem 4, we find:

$$ \lambda^2 |M^{-1/2} \underline{B} \underline{U}^n|_2^2 = \nu^{-2} |\delta \underline{Y}^n|_2^2. $$

As a consequence of the proof of the previous theorem, we infer that the sequence of norms $(|M^{-1/2} \underline{B} \underline{U}^n|_2)_n$ is strictly monotonically decreasing, and so is $(\| \text{div } \underline{u}_h^n \|_{L^2(\Omega)})_n$.

# 6. Discretization and stability estimates

We solve numerically the new variational formulation, written like in Problem (4.9). Below, we consider homogeneous Dirichlet boundary conditions (see Appendix B for nonhomogeneous Dirichlet boundary conditions).

## 6.1. Discretizations

Consider $(\mathcal{T}_h)_h$ a simplicial triangulation sequence of $\Omega$. For all $D \subset \mathbb{R}^d$, and $k \in \mathbb{N}$, we call $P^k(D)$ the set of order $k$ polynomials on $D$, $\mathbf{P}^k(D) = (P^k(D))^d$, and we consider the broken polynomial space:

$$ P_{\text{disc}}^k(\mathcal{T}_h) := \{q \in L^2(\Omega); \forall T \in \mathcal{T}_h, q|_T \in P^k(T)\}. $$

We use the notation $P^0(\mathcal{T}_h)$ for $P_{\text{disc}}^0(\mathcal{T}_h)$. We define the conforming spaces of $P^k$-Lagrange functions:
$$V_h^{k,cg} := \{v_h \in H^1(\Omega); \forall T \in \mathcal{T}_h, v_{h|T} \in P^k(T)\}$$

$$V_{0,h}^{k,cg} := \{v_h \in V_h^{k,cg}; v_{h|\partial\Omega} = 0\}.$$

We set $V_h^{k,cg} := (V_h^{k,cg})^d$ and $V_{0,h}^{k,cg} := (V_{0,h}^{k,cg})^d$. We call $Q_h^k := P_{disc}^k(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$. Notice that $V_{0,h}^{k,cg} \times Q_h^{k'} \subset H_0^1(\Omega) \times L_{zmv}^2(\Omega)$ and $\text{div } V_{0,h}^{k,cg} \subset Q_h^{k-1}$.

Let $k \ge 1$. The (conforming) discretization of Problem (4.9) with the Scott-Vogelius $P^k \times P_{disc}^{k-1}$ finite element [49,12] reads:

$$ \begin{cases} \text{Find } (\mathbf{u}_h, p_h) \in \mathbf{V}_{0,h}^{k,cg} \times Q_h^{k-1} \text{ s.t. for all } (\mathbf{v}_h, q_h) \in \mathbf{V}_{0,h}^{k,cg} \times Q_h^{k-1}: \\ \nu \lambda(\mathbf{u}_h, \mathbf{v}_h) \mathbf{H}_0^1(\Omega) - \lambda(p_h, \text{div } \mathbf{v}_h) L^2(\Omega) = \lambda(\mathbf{f}, \mathbf{v}_h) \mathbf{H}_0^{-1}(\Omega), \\ \lambda(q_h, \text{div } \mathbf{u}_h) L^2(\Omega) + \nu^{-1}(p_h, q_h) L^2(\Omega) = \nu^{-1}(z_{\mathbf{f}}, q_h) L^2(\Omega). \end{cases} \quad (6.1) $$

Writing $\mathbf{f} = -\nu\Delta\mathbf{u} + \text{grad } z_{\mathbf{f}}$ (recall (3.1) with $z_{\mathbf{f}} = p$), Problem (6.1) reads:

$$ \begin{cases} \text{Find } (\mathbf{u}_h, p_h) \in \mathbf{V}_{0,h}^{k,cg} \times Q_h^{k-1} \text{ s.t. for all } (\mathbf{v}_h, q_h) \in \mathbf{V}_{0,h}^{k,cg} \times Q_h^{k-1}: \\ \nu \lambda(\mathbf{u}_h, \mathbf{v}_h) \mathbf{H}_0^1(\Omega) - \lambda(p_h, \text{div } \mathbf{v}_h) L^2(\Omega) = \nu \lambda(\mathbf{u}, \mathbf{v}_h) \mathbf{H}_0^1(\Omega) \\ \qquad \qquad \qquad \qquad - \lambda(z_{\mathbf{f}}, \text{div } \mathbf{v}_h) L^2(\Omega), \\ \lambda(q_h, \text{div } \mathbf{u}_h) L^2(\Omega) + \nu^{-1}(p_h, q_h) L^2(\Omega) = \nu^{-1}(z_{\mathbf{f}}, q_h) L^2(\Omega). \end{cases} \quad (6.2) $$

**Remark 6.** It is well known that discretizing (3.4) with the Scott-Vogelius $P^k \times P_{disc}^{k-1}$ finite element, $1 \le k \le 3$, is not stable on all shape regular meshes [37]. A wide range of strategies to get round this problem have been explored for years. Below is a quick overview of these strategies:

* The discrete velocity space can be enriched see [10,5,8,24,45] for $k=1$.

* The discrete variational formulation can be stabilized, see [11,9,54,2] for $k=1$, [1] for $k=2$.

* The mesh can be designed in such a way that the uniform discrete inf-sup condition applies, see [32] for $k=1$, [49,53,7] for $k>1$ and [37] for $k \ge 4$.

## 6.2. Error estimate when $z_{\mathbf{f}}$ is known

Let $\Pi_{h,cg}$ be the $L^2(\Omega)$-orthogonal projection operator from $L^2(\Omega)$ to $V_{0,h}^k$ and $\Pi_{h,dg}$ be the $L^2(\Omega)$-orthogonal projection operator from $L^2(\Omega)$ to $P_{disc}^{k-1}(\mathcal{T}_h)$. We have the

**Proposition 6.** Let $(\mathbf{u}, p) \in \mathcal{X}$ be the exact solution and $(\mathbf{u}_h, p_h) \in \mathcal{X}_h$ be the solution to Problem (6.1). We have the following error estimate:

$$ \| (\mathbf{u}_h - \mathbf{u}, p_h - \Pi_{h,dg} z_{\mathbf{f}}) \|_{\mathcal{X}, \sqrt{\lambda\nu}} \le \sqrt{1 + d\lambda} \| \mathbf{u} - \Pi_{h,cg} \mathbf{u} \|_{H_0^1(\Omega)}. \quad (6.3) $$

Suppose that $\mathbf{u} \in P^k(\Omega)$. Then $\mathbf{u}_h = \mathbf{u}$ and $p_h = \Pi_{h,dg} z_{\mathbf{f}}$.

Notice that the error estimate for the velocity is independent of the pressure. Hence, the (conforming) discretization of Problem (4.9) with the Scott-Vogelius $P^k \times P_{disc}^{k-1}$ finite element and for which $z_{\mathbf{f}}$ is known is a so called pressure robust method for which the discrete velocity $\mathbf{u}_h$ tends to $\mathbf{u}$ independently of $\nu$. The discrete pressure $p_h$ tends to $\Pi_{h,dg} z_{\mathbf{f}}$ all the faster the smaller $\nu$ is. Let us now prove Proposition 6.

**Proof.** Setting $z_{\mathbf{f},h} = \Pi_{h,dg} z_{\mathbf{f}}$, we have: $(z_{\mathbf{f}}, \text{div } \mathbf{v}_h)_{L^2(\Omega)} = (z_{\mathbf{f},h}, \text{div } \mathbf{v}_h)_{L^2(\Omega)}$ and $(z_{\mathbf{f}}, q_h)_{L^2(\Omega)} = (z_{\mathbf{f},h}, q_h)_{L^2(\Omega)}$ for all $\mathbf{v}_h \in V_{0,h}^k$, and all $q_h \in Q_h^{k-1}$. Summing both equations of Problem (6.2) and reshuffling terms, it comes:

$$ \begin{cases} \text{Find } (\mathbf{u}_h, p_h) \in \mathbf{V}_{0,h}^k \times Q_h^{k-1} \text{ s.t. for all } (\mathbf{v}_h, q_h) \in \mathbf{V}_{0,h}^k \times Q_h^{k-1}: \\ \nu \lambda(\mathbf{u}_h - \mathbf{u}, \mathbf{v}_h) \mathbf{H}_0^1(\Omega) + \nu^{-1}(p_h - z_{\mathbf{f},h}, q_h)_{L^2(\Omega)} \\ \qquad = \lambda(p_h - z_{\mathbf{f},h}, \text{div } \mathbf{v}_h)_{L^2(\Omega)} - \lambda(q_h, \text{div } \mathbf{u}_h)_{L^2(\Omega)}. \end{cases} \quad (6.4) $$

Choosing $\mathbf{v}_h = \mathbf{u}_h - \Pi_{h,cg} \mathbf{u} = (\mathbf{u}_h - \mathbf{u}) + (\mathbf{u} - \Pi_{h,cg} \mathbf{u})$ and $q_h = p_h - z_{\mathbf{f},h}$ in (6.4), and dividing by $\lambda\nu$, we obtain:

$$ \begin{aligned} & \|\mathbf{u}_h - \mathbf{u}\|_{H_0^1(\Omega)}^2 + (\mathbf{u}_h - \mathbf{u}, \mathbf{u} - \Pi_{h,cg} \mathbf{u})_{H_0^1(\Omega)} + \lambda^{-1} \nu^{-2} \|p_h - z_{\mathbf{f},h}\|_{L^2(\Omega)}^2 \\ & \qquad = -\nu^{-1}(p_h - z_{\mathbf{f},h}, \text{div } \Pi_{h,cg} \mathbf{u})_{L^2(\Omega)}. \end{aligned} $$

Hence, using Cauchy-Schwarz inequality, we get:
$$
\begin{align*}
& \|\mathbf{u}_h - \mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}^2 + \lambda^{-1} v^{-2} \|p_h - z_{\mathbf{f},h}\|_{L^2(\Omega)}^2 \\
& \le \|\mathbf{u}_h - \mathbf{u}\|_{\mathbf{H}_0^1(\Omega)} \|\mathbf{u} - \Pi_{h,\text{cg}} \mathbf{u}\|_{\mathbf{H}_0^1(\Omega)} + v^{-1} \|p_h - z_{\mathbf{f},h}\|_{L^2(\Omega)} \|\text{div} \Pi_{h,\text{cg}} \mathbf{u}\|_{L^2(\Omega)} \\
& \le (\|\mathbf{u}_h - \mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}^2 + \lambda^{-1} v^{-2} \|p_h - z_{\mathbf{f},h}\|_{L^2(\Omega)}^2)^{1/2} \times \\
& \qquad (\|\mathbf{u} - \Pi_{h,\text{cg}} \mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}^2 + \lambda \|\text{div} \Pi_{h,\text{cg}} \mathbf{u}\|_{L^2(\Omega)}^2)^{1/2}
\end{align*}
$$

Finally $\| \operatorname{div} \Pi_{h,cg} \mathbf{u} \|_{L^2(\Omega)} \le \sqrt{d} \| \mathbf{u} - \Pi_{h,cg} \mathbf{u} \|_{\mathbf{H}_0^1(\Omega)}$, so we have:

$$
\|(\mathbf{u}_h - \mathbf{u}, p_h - z_{\mathbf{f},h})\|_{X, \sqrt{\lambda v}} \le \sqrt{1 + d\lambda} \|\mathbf{u} - \Pi_{h,cg} \mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}.
$$

If $\mathbf{u} \in \mathbf{P}^k(\Omega)$, then $\Pi_{h,cg}\mathbf{u} = \mathbf{u}$ and the right-hand side is equal to 0. Hence, $\mathbf{u}_h = \mathbf{u}$ and $p_h = z_{\mathbf{f},h}$. $\square$

6.3. Error estimate when z<sub>f</sub> is not known

If $z_f$ is not known explicitly, let us assume that we have at hand an approximation $p_h^{\text{old}} \in Q_h^{k-1}$ and use it to solve the following approximation of Problem (4.9):

$$
\begin{equation}
\left\{
\begin{array}{@{}l@{}}
\text{Find } (\tilde{\mathbf{u}}, \tilde{\rho}) \in \mathcal{X} \text{ s.t. for all } (\mathbf{v}, q) \in \mathcal{X}: \\
\nu \lambda (\tilde{\mathbf{u}}, \mathbf{v})_{\mathbf{H}_0^1(\Omega)} - \lambda(\tilde{\rho}, \operatorname{div} \mathbf{v})_{L^2(\Omega)} = \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)}, \\
\lambda(q, \operatorname{div} \tilde{\mathbf{u}})_{L^2(\Omega)} + \nu^{-1} (\tilde{\rho}, q)_{L^2(\Omega)} = \nu^{-1} (p_h^{\text{old}}, q)_{L^2(\Omega)}.
\end{array}
\right.
\tag{6.5}
\end{equation}
$$

Let us write again **f** = -νΔ**u** + **grad** z<sub>f</sub>, so that Problem (6.5) reads:

$$
\left\{
\begin{array}{l}
\text{Find } (\tilde{\mathbf{u}}, \tilde{\rho}) \in \mathcal{X} \text{ s.t. for all } (\mathbf{v}, q) \in \mathcal{X}: \\
\nu \lambda (\tilde{\mathbf{u}} - \mathbf{u}, \mathbf{v})_{\mathbf{H}_0^1(\Omega)} - \lambda (\tilde{\rho} - z_{\mathbf{f}}, \operatorname{div} \mathbf{v})_{L^2(\Omega)} = 0, \\
\lambda (q, \operatorname{div} \tilde{\mathbf{u}})_{L^2(\Omega)} + \nu^{-1} (\tilde{\rho} - p_h^{\text{old}}, q)_{L^2(\Omega)} = 0.
\end{array}
\right.
\qquad (6.6)
$$

**Proposition 7.** Let (u, p) ∈ 𝒳 be the exact solution and (tilde{u},tilde{p}) ∈ 𝒳 be the solution to Problem (6.5). We have the following estimates:

$$
\left\{
\begin{array}{l}
\| \tilde{\rho} - z_{\mathbf{f}} \|_{L^2(\Omega)} \le \| z_{\mathbf{f}} - p_h^{\text{old}} \|_{L^2(\Omega)}, \\
\| \tilde{\mathbf{u}} - \mathbf{u} \|_{\mathbf{H}_0^1(\Omega)} \le (\sqrt{\lambda \nu})^{-1} \| z_{\mathbf{f}} - p_h^{\text{old}} \|_{L^2(\Omega)}.
\end{array}
\right.
\qquad (6.7)
$$

**Proof.** Choosing **v** = **u** − **u** and **q** = **p** − **z**<sub>f</sub> in (6.6), summing both equations and dividing by **λν**, it comes:

$$
\| \tilde{\mathbf{u}} - \mathbf{u} \|_{\mathbf{H}_0^1(\Omega)}^2 + \lambda^{-1} v^{-2} (\tilde{\rho} - p_h^{\text{old}}, \tilde{\rho} - z_{\mathbf{f}})_{L^2(\Omega)} = 0.
$$

Writing $\tilde{\rho} - p_h^{\text{old}} = (\tilde{\rho} - z_f) + (z_f - p_h^{\text{old}})$, we obtain:

$$
\begin{align*}
\|(\tilde{\mathbf{u}} - \mathbf{u}, \tilde{\rho} - z_{\mathbf{f}})\|_{\mathcal{X}, \sqrt{\lambda \nu}}^2 &= -\lambda^{-1} v^{-2} (z_{\mathbf{f}} - p_h^{\text{old}}, \tilde{\rho} - z_{\mathbf{f}})_{L^2(\Omega)} \\
&\le \lambda^{-1} v^{-2} \|z_{\mathbf{f}} - p_h^{\text{old}}\|_{L^2(\Omega)} \|\tilde{\rho} - z_{\mathbf{f}}\|_{L^2(\Omega)}.
\end{align*}
$$

We obtain successively estimates (6.7).

<figure><img src="image_6.png" /></figure>

The discretization of Problem (6.5) with the Scott-Vogelius $\mathbf{P}^k \times \mathbf{P}_{\text{disc}}^{k-1}$ finite element reads:

$$
\left\{
\begin{array}{l}
\text{Find } (\tilde{\mathbf{u}}_h, p_h^{\text{new}}) \in \mathbf{V}_{0,h}^k \times \mathbf{Q}_h^{k-1} \text{ s.t. for all } (\mathbf{v}_h, q_h) \in \mathbf{V}_{0,h}^k \times \mathbf{Q}_h^{k-1}: \\
\nu \lambda (\tilde{\mathbf{u}}_h, \mathbf{v}_h)_{\mathbf{H}_0^1(\Omega)} - \lambda (p_h^{\text{new}}, \operatorname{div} \mathbf{v}_h)_{L^2(\Omega)} = \lambda \langle \mathbf{f}, \mathbf{v}_h \rangle_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^1(\Omega)}, \\
\lambda (q_h, \operatorname{div} \tilde{\mathbf{u}}_h)_{L^2(\Omega)} + \nu^{-1} (p_h^{\text{new}}, q_h)_{L^2(\Omega)} = \nu^{-1} (p_h^{\text{old}}, q_h)_{L^2(\Omega)}.
\end{array}
\right.
\qquad (6.8)
$$

**Theorem 5.** Let (u, p) ∈ 𝒳 be the exact solution and (tilde{u}_h, p_h^{new}) ∈ 𝒳 be the solution to Problem (6.8). We have the following error estimate:

$$
\| (\tilde{\mathbf{u}}_h - \mathbf{u}, p_h^{\text{new}} - \Pi_{h,dg} z_{\mathbf{f}}) \|_{X, \sqrt{\lambda \nu}} \le \sqrt{1 + d\lambda} \| \mathbf{u} - \Pi_{h,cg} \mathbf{u} \|_{\mathbf{H}_0^1(\Omega)} \\
+ (\sqrt{\lambda \nu})^{-1} \| \Pi_{h,dg} z_{\mathbf{f}} - p_h^{\text{old}} \|_{L^2(\Omega)} \tag{6.9}
$$

Suppose that $\mathbf{u} \in \mathbf{P}^k(\Omega)$. Then $\Pi_{h,cg}\mathbf{u} = \mathbf{u}$ and we obtain the estimate below:

$$
\| (\tilde{\mathbf{u}}_h - \mathbf{u}, p_h^{\text{new}} - \Pi_{h,dg} z_{\mathbf{f}}) \|_{X, \sqrt{\lambda \nu}} \le (\sqrt{\lambda \nu})^{-1} \| \Pi_{h,dg} z_{\mathbf{f}} - p_h^{\text{old}} \|_{L^2(\Omega)} \quad (6.10)
$$
In particular, if $p_h^{old}$ is a good approximation of $\Pi_{h,dg} z_f$, the solution $(\tilde{u}_h, p_h^{new})$ is also a good approximation of $(u, \Pi_{h,dg} z_f)$. Interestingly, the above with $k=1$ corresponds to the $P^1 \times P^0$ finite element pair, which is known to be unstable for the discretization of the usual variational formulation (3.4) of the Stokes problem. Let us prove Theorem 5.

**Proof.** Setting $z_{\mathbf{f},h} = \Pi_{h,dg} z_{\mathbf{f}}$, we have: $\langle \mathbf{f}, v_h \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = \nu(u, v_h)_{H_0^1(\Omega)} - (z_{\mathbf{f},h}, \text{div } v_h)_{L^2(\Omega)}$ for all $v_h \in V_{0,h}^k$. Hence, Problem (6.8) can be written as:

$$ \left\{ \begin{array}{l} \text{Find } (\tilde{u}_h, p_h^{new}) \in V_{0,h}^k \times Q_h^{k-1} \text{ s.t. for all } (v_h, q_h) \in V_{0,h}^k \times Q_h^{k-1} : \\ \nu \lambda (\tilde{u}_h - u, v_h)_{H_0^1(\Omega)} - \lambda (p_h^{new} - z_{\mathbf{f},h}, \text{div } v_h)_{L^2(\Omega)} = 0, \\ \lambda (q_h, \text{div } \tilde{u}_h)_{L^2(\Omega)} + \nu^{-1} (p_h^{new} - p_h^{old}, q_h)_{L^2(\Omega)} = 0. \end{array} \right. \qquad (6.11) $$

Choosing $(v_h, q_h) = (\tilde{u}_h - \Pi_{h,cg} u, p_h^{new} - z_{\mathbf{f},h})$ in Problem (6.11), we have:

$$ \left\{ \begin{array}{l} \nu \lambda (\tilde{u}_h - u, \tilde{u}_h - \Pi_{h,cg} u)_{H_0^1(\Omega)} - \lambda (p_h^{new} - z_{\mathbf{f},h}, \text{div}(\tilde{u}_h - \Pi_{h,cg} u))_{L^2(\Omega)} = 0, \\ \lambda (p_h^{new} - z_{\mathbf{f},h}, \text{div} \tilde{u}_h)_{L^2(\Omega)} + \nu^{-1} (p_h^{new} - p_h^{old}, p_h^{new} - z_{\mathbf{f},h})_{L^2(\Omega)} = 0. \end{array} \right. $$

Summing both equations and dividing by $\lambda\nu$, it now comes:

$$ \begin{aligned} & (\tilde{u}_h - u, \tilde{u}_h - \Pi_{h,cg} u)_{H_0^1(\Omega)} + \lambda^{-1} \nu^{-2} (p_h^{new} - p_h^{old}, p_h^{new} - z_{\mathbf{f},h})_{L^2(\Omega)} \\ & \qquad + \nu^{-1} (p_h^{new} - z_{\mathbf{f},h}, \text{div} \Pi_{h,cg} u)_{L^2(\Omega)} = 0. \end{aligned} $$

Noticing that $\tilde{u}_h - \Pi_{h,cg} u = (\tilde{u}_h - u) + (u - \Pi_{h,cg} u)$ and $p_h^{new} - p_h^{old} = (p_h^{new} - z_{\mathbf{f},h}) + (z_{\mathbf{f},h} - p_h^{old})$, we get:

$$ \begin{aligned} & \| (\tilde{u}_h - u, p_h^{new} - z_{\mathbf{f},h}) \|_{\lambda, \sqrt{\lambda\nu}}^2 + (\tilde{u}_h - u, u - \Pi_{h,cg} u)_{H_0^1(\Omega)} \\ & \qquad + \lambda^{-1} \nu^{-2} (z_{\mathbf{f},h} - p_h^{old}, p_h^{new} - z_{\mathbf{f},h})_{L^2(\Omega)} + \nu^{-1} (p_h^{new} - z_{\mathbf{f},h}, \text{div} \Pi_{h,cg} u)_{L^2(\Omega)} = 0. \end{aligned} $$

Using Cauchy-Schwarz inequality, we deduce that:

$$ \begin{aligned} & \| (\tilde{u}_h - u, p_h^{new} - z_{\mathbf{f},h}) \|_{\lambda, \sqrt{\lambda\nu}}^2 \le \| \tilde{u}_h - u \|_{H_0^1(\Omega)} \| u - \Pi_{h,cg} u \|_{H_0^1(\Omega)} \\ & \qquad + (\sqrt{\lambda\nu})^{-1} \| p_h^{new} - z_{\mathbf{f},h} \|_{L^2(\Omega)} \left( (\sqrt{\lambda\nu})^{-1} \| z_{\mathbf{f},h} - p_h^{old} \|_{L^2(\Omega)} + \sqrt{\lambda} \| \text{div} \Pi_{h,cg} u \|_{L^2(\Omega)} \right). \end{aligned} $$

Using again Cauchy-Schwarz inequality as in the proof of Proposition 6, together with $\| \text{div} \Pi_{h,cg} u \|_{L^2(\Omega)} \le \sqrt{\nu} \| u - \Pi_{h,cg} u \|_{H_0^1(\Omega)}$, we obtain the estimates (6.9) and (6.10). $\square$

## 6.4. Discussion on the error indicators

As we noted after stating Theorem 5, if the error indicator on the pressure $\|z_{\mathbf{f},h} - p_h^{old}\|_{L^2(\Omega)}$ is small, so is $\|z_{\mathbf{f},h} - p_h^{new}\|_{L^2(\Omega)}$. Indeed, assuming that the errors on the velocities are negligible compared to these indicators, one derives from (6.9) the bound $\|z_{\mathbf{f},h} - p_h^{new}\|_{L^2(\Omega)} \le \|z_{\mathbf{f},h} - p_h^{old}\|_{L^2(\Omega)}$. Interestingly, this situation is likely to occur when the viscosity is small, since there is a multiplicative factor equal to $(\sqrt{\lambda\nu})^{-1}$ in front of the error indicators on the pressure. This further justifies the use of the iterative algorithm presented in Section 5 in this case.

## 7. Numerical results

The numerical algorithms presented here are written using the same notation as in Section 5.

We propose some numerical experiments with $g=0$, depending on whether or not $z_{\mathbf{f}}$ is known explicitly. In the latter case, we first compute some approximation $p_h^{old}$. In principle, either a classical, conforming or nonconforming, discretization can be used to compute $p_h^{old}$. Hence, we apply an iterative approach with a nonconforming initial computation followed by (several iterations of) our new formulation. We call the second part the post-processing. For each test-case, $(u, p)$ is given, whereas the viscosity $\nu$ can vary (see the discussion after (3.5)). The numerical results are obtained on a github platform, implemented in Octave language, see [42].

### 7.1. Resolution algorithm

Consider a nonconforming discretization of the classical variational formulation (3.4), cf. [29]. Let $V_h^{k,nc} \subset H^1(\Omega)$ (resp. $V_{0,h}^{k,nc} \subset H_0^1(\Omega)$) be the nonconforming discretization space of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$). We set $V_h^{k,nc} = (V_h^{k,nc})^d$ and $V_{0,h}^{k,nc} = (V_{0,h}^{k,nc})^d$. The discrete velocity space is $V_{0,h}^{k,nc}$ and the discrete pressure space is $Q_h^{k-1}$. Let $(\psi_i)_{i=1}^{N_u}$ be a basis of $V_{0,h}^{k,nc}$, and $(\phi_i)_{i=1}^{N_p}$ be a basis of $Q_h^{k-1}$. We set: $\underline{U} := (\underline{U}_i)_{i=1}^{N_u}$ where $u_h := \sum_{i=1}^{N_u} \underline{U}_i \psi_i$, $\underline{P} := (\underline{P}_i)_{i=1}^{N_p}$ where $p_h := \sum_{i=1}^{N_p} \underline{P}_i \phi_i$, and $\underline{F}_u := (\underline{F}_{u,t})_{t=1}^{N_u}$. We let $\ell_{\mathbf{f}} \in \mathcal{L}(V_{0,h}^{k,nc}, \mathbb{R})$ be such that $\forall v_h \in V_{0,h}^{k,nc}$, $\ell_{\mathbf{f}}(v_h) = (\mathbf{f}, v_h)_{L^2(\Omega)}$ if $\mathbf{f} \in L^2(\Omega)$, $\ell_{\mathbf{f}}(v_h) = (\mathbf{f}, J_h(v_h))_{H^{-1}(\Omega), H_0^1(\Omega)}$ if $\mathbf{f} \notin L^2(\Omega)$, where $J_h: V_{0,h}^{k,nc} \to V_{0,h}^k$ is for instance
an averaging operator [30, §22.4.1]. Then $\underline{F}_{u,i} = \ell_f(\psi_i)$. Let $\mathcal{A} \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_u}$ be the velocity stiffness matrix, $\mathbb{B} \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_p}$ be the velocity-pressure coupling matrix and $\mathcal{M} \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_p}$ be the pressure mass matrix. The linear system to be solved is

$$ \left\{ \begin{array}{l} \text{Find } (\underline{U}, \underline{P}) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\ \nu \mathcal{A} \underline{U} - \mathbb{B}^T \underline{P} = \underline{F}_{u} \\ \mathbb{B} \underline{U} = 0 \end{array} \right. . \qquad (7.1) $$

Let us set $\mathcal{K} = \mathbb{B} \mathcal{A}^{-1} \mathbb{B}^T \in \mathbb{R}^{N_p} \times \mathbb{R}^{N_p}$. The matrix $\mathcal{K}$ is a symmetric matrix, furthermore it is positive definite as soon as the kernel of $\mathbb{B}^T$ is reduced to $\{0\}$. When it is the case, the coupled velocity-pressure problem (7.1) can be solved using the three + one steps below (the fourth step being straightforward):

$$ \begin{array}{ll} \text{Prediction:} & \text{Solve in } \underline{U}_{\star} \text{ such that } \nu \mathcal{A} \underline{U}_{\star} = \underline{F}_{u} \\ \text{Pressure solver:} & \text{Solve in } \underline{P} \text{ such that } \mathcal{K} \underline{P} = -\nu \mathbb{B} \underline{U}_{\star} \\ \text{Correction:} & \text{Solve in } \delta \underline{U} \text{ such that } \nu \mathcal{A} \delta \underline{U} = \mathbb{B}^T \underline{P}. \\ \text{Update:} & \underline{U} = \delta \underline{U} + \underline{U}_{\star}. \end{array} \qquad (7.2) $$

One can check easily that the above computed solution $(\underline{U}, \underline{P})$ solves (7.1). The pressure solver with matrix $\mathcal{K}$ is based on the Uzawa algorithm, which is the conjugate gradient algorithm in the context of the Stokes problem. It can be preconditioned by the inverse of the mass matrix associated to the discrete pressure $\mathcal{M}$ (see e.g. [48, Lemma 5.9]). Thanks to the uniform discrete inf-sup condition, the number of iterations of the conjugate gradient algorithm is independent of the meshsize. The matrices $\mathcal{M}$ and $\mathbb{B}$ are kept with all the $P_{\text{disc}}^{k-1}$ degrees of freedom for the discrete pressure. To take into account the zero mean value constraint, at each iteration of the preconditioned conjugate gradient algorithm, the discrete pressure is first computed in $L^2(\Omega)$, then orthogonally projected in $L^2_{\text{zmv}}(\Omega)$. We use the Cholesky factorization (computed once and for all) to solve linear systems with matrix $\mathcal{A}$.

Consider next the variational formulation (6.1) or (6.8). Let $(\psi_i)_{i=1}^{N_u}$ be the Lagrange basis of $\mathbf{V}_{0,h}^k$, and $(\phi_i)_{i=1}^{N_p}$ be a basis of $Q_h^{k-1}$. We set: $\underline{U} := (\underline{U}_i)_{i=1}^{N_u}$ where $\mathbf{u}_h := \sum_{i=1}^{N_u} \underline{U}_i \psi_i$ and $\underline{P} := (\underline{P}_i)_{i=1}^{N_p}$ where $p_h := \sum_{i=1}^{N_p} \underline{P}_i \phi_i$. We set $\underline{F}_{u} := (\underline{F}_{u,i})_{i=1}^{N_u}$ where $\underline{F}_{u,i} = (\underline{f}_i, \psi_i)_{\mathbf{H}^{-1}(\Omega), \mathbf{H}_0^{1}(\Omega)}$ and $\underline{F}_{p} := (\underline{F}_{p,i})_{i=1}^{N_p}$ where $\underline{F}_{p,i} = (z_f, \phi_i)_{L^2(\Omega)}$, cf. (6.1), or $\underline{F}_{p,i} = (p_h^{old}, \phi_i)_{L^2(\Omega)}$, where the discrete pressure $p_h^{old}$ is an approximation of $z_f$, cf. (6.8). The linear system to be solved is

$$ \left\{ \begin{array}{l} \text{Find } (\underline{U}, \underline{P}) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\ \nu \lambda \mathcal{A} \underline{U} - \lambda \mathbb{B}^T \underline{P} = \lambda \underline{F}_{u} \\ \lambda \mathbb{B} \underline{U} + \nu^{-1} \mathcal{M} \underline{P} = \nu^{-1} \underline{F}_{p} \end{array} \right. . \qquad (7.3) $$

In that case, we set $\mathcal{A}_{\lambda} = \mathcal{A} + \lambda \mathbb{B}^T \mathcal{M}^{-1} \mathbb{B} \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_u}$, which is automatically a symmetric positive definite matrix. To solve the coupled velocity-pressure problem (7.3), one relies on the two steps below:

$$ \begin{array}{ll} \text{Velocity solver:} & \text{Solve in } \underline{U} \text{ such that } \nu \mathcal{A}_{\lambda} \underline{U} = \underline{F}_{u} + \mathbb{B}^T \mathcal{M}^{-1} \underline{F}_{p} \\ \text{Pressure solver:} & \text{Solve in } \underline{P} \text{ such that } \mathcal{M} \underline{P} = \underline{F}_{p} - \nu \lambda \mathbb{B} \underline{U}. \end{array} \qquad (7.4) $$

One can check easily that the above computed solution $(\underline{U}, \underline{P})$ solves (7.3). The matrix $\mathcal{A}_{\lambda}$ can be assembled easily since $\mathcal{M}$ is a block-diagonal matrix. Using $\mathbf{P}^1 \times \mathbf{P}^0$ discretization, the computation of $\underline{P}$ is straightforward. Again, matrices $\mathcal{M}$ and $\mathbb{B}$ are kept with all the $P_{\text{disc}}^{k-1}$ degrees of freedom for the discrete pressure. We proceed as before to take into account the zero mean value constraint.

**Remark 7.** We recall that solving the linear system (7.1) via the solver (7.2) is not so straightforward, even with $\nu = 1$. Solving the linear system (7.3) is much easier.

## 7.2. Settings

We consider Problem (3.1) with homogeneous or nonhomogeneous boundary conditions in $\Omega = (0, 1)^2$. Let $(\mathbf{u}, p)$ be the exact solution, and $(\mathbf{u}_h, p_h)$ be the numerical solution. We compare numerical methods, showing how the coercive $\mathbf{P}^1 \times \mathbf{P}^0$ formulation can be used in a post-processing step (i.e. solving (6.8) with $p_h^{old}$ known), improving then the approximation of the discrete velocity. On Figs. 2 to 11, and in Tables 3 to 7 we give the following names to our numerical methods:

* Method with Crouzeix-Raviart (CR): computations are made with the nonconforming Crouzeix-Raviart $\mathbf{P}_{nc}^1 \times \mathbf{P}^0$ formulation [29, Example 4], which is not a pressure robust method. We call $p_{nc}$ the resulting discrete pressure.

* Method with exact pressure (EP): computations are made with the coercive $\mathbf{P}^1 \times \mathbf{P}^0$ formulation (6.1) and $\lambda = 1$, knowing the exact pressure $z_f$ (i.e. we solve Problem (6.1)).

* Method with post-processing (Post): computations are made in two steps. In a first step, we approximate the pressure by the CR-method. Then, in a post-processing step, we use this numerical pressure as the source term in the EP-method (i.e. we solve Problem (6.8) with $p_h^{old} = p_{nc}$) and $\lambda = 10$. Unless the stopping criterion (5.6) is achieved, we iterate eight times the second step, updating $p_h^{old}$ at each new iteration. The algorithm corresponds to Algorithm (5.2).
<table><caption>**Table 1**
Number of unknowns: first two test cases (left), last test case (right).</caption><thead><tr><th>h</th><th># dof CR</th><th># dof EP</th><th># dof CR</th><th># dof EP</th></tr></thead><tbody><tr><td>1.00 × 10<sup>-1</sup></td><td>1 048</td><td>566</td><td>2 184</td><td>1 114</td></tr><tr><td>5.00 × 10<sup>-2</sup></td><td>4 376</td><td>2 270</td><td>8 064</td><td>4 074</td></tr><tr><td>2.50 × 10<sup>-2</sup></td><td>17 368</td><td>8 846</td><td>30 464</td><td>15 314</td></tr><tr><td>1.25 × 10<sup>-2</sup></td><td>67 816</td><td>34 230</td><td>117 664</td><td>58 994</td></tr><tr><td>6.25 × 10<sup>-3</sup></td><td>272 624</td><td>136 954</td><td>464 448</td><td>232 866</td></tr></tbody></table>

Let $N_T$ be the number of triangles. For the velocity, the number of unknowns is of order $N_u \approx 3 N_T$ for the CR-method and $N_u \approx N_T$ for the EP-method. For the pressure, the number of unknowns is $N_p = N_T - 1$ for both methods. As a consequence, there are roughly twice as many unknowns for the CR-method than there are for the EP-method. We report in Table 1 the number of unknowns "# dof" for the numerical tests.

We propose three numerical examples based on manufactured solutions. Let $I_{h,cg}$ be the interpolation operator from $C^0(\bar{\Omega})$ to $V_h^{1,cg}$ and $I_{h,nc}$ be the interpolation operator from $C^0(\bar{\Omega})$ to $V_h^{1,nc}$. We use the discrete $L^2$-error for the velocity and for the pressure:

$$ 
\begin{aligned}
\varepsilon_0^v(\mathbf{u}_h) &:= \| I_h \mathbf{u} - \mathbf{u}_h \|_{L^2(\Omega)} / \| (\mathbf{u}, p) \|_{X,v} \\
\varepsilon_0^v(p_h) &:= v^{-1} \| \Pi_{h,dg} p - p_h \|_{L^2(\Omega)} / \| (\mathbf{u}, p) \|_{X,v}
\end{aligned} 
\quad (7.5) $$

where $I_h$ is either equal $I_{h,nc}$ or $I_{h,cg}$. For the EP-method and the Post-method, we moreover compare the indicators $\varepsilon_D^v$ and $\varepsilon_1^v$ defined by:

$$ 
\begin{aligned}
\varepsilon_1^v &:= \| \mathrm{Grad}(I_{h,cg} \mathbf{u} - \mathbf{u}_h) \|_{L^2(\Omega)} / \| (\mathbf{u}, p) \|_{X,v} \\
\varepsilon_D^v &:= \| \mathrm{div} \, \mathbf{u}_h \|_{L^2(\Omega)} / \| (\mathbf{u}, p) \|_{X,v}
\end{aligned} 
\quad (7.6) $$

We first plot the results as error curves for both the velocity and the pressure as a function of the meshsize, and we give the numerically observed average convergence rates. Second, we report the elapsed CPU times for the CR-method and the Post-method. In the Tables, column "overhead" indicates the ratio between the cost of the second (post-processing) step for the Post-method and the cost of the first step (CR-method). Finally, we plot the errors of the two methods as a function of the elapsed CPU times of the resolution algorithm.

For the Post-method "Post" plots represent the errors after a single iteration of the second step and "Post-8" plots represent the errors after iterating the second step eight times, updating $p_h^{old}$ at each new iteration, cf. Algorithm (5.2).

**Remark 8.** The use of discrete mixed formulations with $H(\mathrm{div})$-conforming projection of the test function on the right-hand side leads to a pressure robust discrete velocity [46,20]. It has been proven to be an efficient numerical solution to solve Stokes problem when the parameter $v$ is small, however the $H(\mathrm{div})$-conforming discrete spaces must be chosen with care. The Post-method can also be used in that case to obtain an $H^1$-conforming reconstruction of the discrete velocity, see §7.5.

### 7.3. Regular manufactured solutions

We postulate that for the EP-method and the Post-method, $\varepsilon_0^v(\mathbf{u}_h) \le C_0 h^2$, where $C_0$ is independent of the meshsize. The error estimates for the CR-method are given by [29, Theorems 3, 4, 6]. Notice that the norm $\|(\cdot, \cdot)\|_{X,v}$ depends on the viscosity $v$ (3.5) and that the error estimates on the pressure and the velocity are linked. We check the convergence rates and the viscosity dependency of the error estimates. The first test enters the framework of §6.2, see the estimate (6.3).

#### Test case with a linear velocity: Figs. 1-4, Tables 2-3

We consider Problem (3.1)$_{NH}$ with $f = \mathrm{grad} p$, where: $\mathbf{u} = (-y, x)^T$ and $p = x^3 + y^3 - 1/2$.

The number of unknowns are reported on Table 1 (left).

As expected in (6.3), the EP-method returns $\varepsilon_0^v(\mathbf{u}_h) = \mathcal{O}(10^{-15})$, $\varepsilon_0^v(p_h) = \mathcal{O}(10^{-14})$ for $v = 1$; and $\varepsilon_0^v(\mathbf{u}_h) = \mathcal{O}(10^{-18})$, $\varepsilon_0^v(p_h) = \mathcal{O}(10^{-17})$ for $v = 10^{-6}$, so the errors are not reported.

Fig. 1 (resp. 2) shows the discrete error values $\varepsilon_0^v(\mathbf{u}_h)$ (left) and $\varepsilon_0^v(p_h)$ (right) plotted against the meshsize for $v = 1$ (resp. $v = 10^{-6}$). Notice that the Post-method improves the approximation initially computed with the CR-method.

Fig. 3 shows the discrete error values $\varepsilon_0^v$ and $\varepsilon_D^v$ plotted against the meshsize for $v = 10^{-6}$. Notice that the ratio $\varepsilon_D^v/\varepsilon_1^v$ is of order 0.5 for the Post-method with a single iteration and of order 0.1 for the Post-method with eight iterations.

Table 2 shows the average convergence rates for the velocity, $\tau_\mathbf{u}$ and the pressure, $\tau_p$. For the CR-method, the average convergence rate for the pressure, $\tau_p$ is better than expected, possibly because the source term is a polynomial of degree 3 which gradient is numerically exactly integrated.

Table 3 shows the CPU times for the CR-method and the Post-method with either one, or eight, iterations. With our stopping criterion, around 30 iterations of the preconditioned conjugate gradient algorithm are needed to solve the pressure solver of algorithm (7.2) with the CR-method: this is consistent with Remark 7. By design, the Post-method which includes the CR-method as its first step
<figure><img src="image_7.png" /><figcaption>Fig. 1. Linear velocity, $v = 1$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right).</figcaption></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 2. Linear velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right).</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 3. Linear velocity, $v = 10^{-6}$. Plots of $\epsilon_1^v$ (left) and $\epsilon_D^v$ (right).</figcaption></figure>

<figure><img src="image_12.png" /></figure>

requires more CPU time. However, we notice that the CPU time of the second (post-processing) step, the overhead, is only a small fraction of the first step, and also that it decreases dramatically as the meshsize decreases. This can be explained by the fact that on the one hand, there are fewer unknowns and, on the other hand, algorithm (7.4) is faster than algorithm (7.2) (cf. Remark 7). For the Post-method, it seems worth doing eight iterations, especially when the meshsize is small.

Fig. 4 shows the discrete error values $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right) plotted against the CPU time for $v = 10^{-6}$. To reach $\epsilon_0^v(\mathbf{u}_h) \lesssim 2 \times 10^{-5}$, the required CPU time is 0.5 s for the Post-method with one iteration and less than 0.01 s for the Post-method with eight iterations, to be compared with 20 s for the CR-method.
<table><caption>**Table 2**
Linear velocity. Average convergence rates.</caption><thead><tr><th>ν</th><th>τ</th><th>CR</th><th>Post</th><th>Post-8</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>τ<sub>u</sub></td><td>2.03</td><td>1.99</td><td>1.93</td></tr><tr><td>τ<sub>p</sub></td><td>1.55</td><td>1.59</td><td>1.71</td></tr><tr><td>10<sup>-6</sup></td><td>τ<sub>u</sub></td><td>2.03</td><td>1.97</td><td>2.53</td></tr><tr><td></td><td>τ<sub>p</sub></td><td>1.55</td><td>1.68</td><td>1.78</td></tr></tbody></table>

<table><caption>**Table 3**
Linear velocity, $v = 10^{-6}$. CPU time (s).</caption><thead><tr><th>h</th><th>CPU CR</th><th>CPU Post</th><th>overhead</th><th>CPU Post-8</th><th>overhead</th></tr></thead><tbody><tr><td>1.00 × 10<sup>-1</sup></td><td>4.34 × 10<sup>-3</sup></td><td>5.08 × 10<sup>-3</sup></td><td>17%</td><td>7.40 × 10<sup>-3</sup></td><td>70%</td></tr><tr><td>5.00 × 10<sup>-2</sup></td><td>2.32 × 10<sup>-2</sup></td><td>2.61 × 10<sup>-2</sup></td><td>13%</td><td>4.46 × 10<sup>-2</sup></td><td>92%</td></tr><tr><td>2.50 × 10<sup>-2</sup></td><td>2.51 × 10<sup>-1</sup></td><td>2.69 × 10<sup>-1</sup></td><td>7.2%</td><td>3.65 × 10<sup>-1</sup></td><td>45%</td></tr><tr><td>1.25 × 10<sup>-2</sup></td><td>3.00 × 10<sup>+0</sup></td><td>3.07 × 10<sup>+0</sup></td><td>2.3%</td><td>3.54 × 10<sup>+0</sup></td><td>10%</td></tr><tr><td>6.25 × 10<sup>-3</sup></td><td>4.89 × 10<sup>+1</sup></td><td>4.94 × 10<sup>+1</sup></td><td>1.0%</td><td>5.20 × 10<sup>+1</sup></td><td>6.3%</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 4.** Linear velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right) against CPU time (s).</figcaption></figure>

• Test case with a sinusoidal solution: Figs. 5-8, Tables 5-4

We consider Problem (3.1)$_H$ with $f = -v \Delta \mathbf{u} + \text{grad } p$, where:

$$ \mathbf{u} = \begin{pmatrix} (1 - \cos(2 \pi x)) \sin(2 \pi y) \\ (\cos(2 \pi y) - 1) \sin(2 \pi x) \end{pmatrix} \quad \text{and} \quad p = \sin(2 \pi x) \sin(2 \pi y). $$

The number of unknowns are reported on Table 1 (left).

Fig. 5 shows $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right) plotted against the meshsize, for $v = 1$. In that case, for the Post-method, there is no need doing eight iterations, so we give numerical results for computations with a single iteration only. The Post-method and the EP-method give similar velocity errors. The CR-method and the Post-method give similar pressure errors. For a given meshsize $h$, the velocity error is smaller for the EP-method and the Post-method than for the CR-method; and the pressure error is smaller for the CR-method than for the EP-method, but it is obtained at a higher cost. Notice that the Post-method reduces the velocity error without worsening the pressure error.

Fig. 6 shows $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right) plotted against the meshsize, now for $v = 10^{-6}$. Both post-processings improve the initial computation. Even with eight iterations, the overhead cost remains affordable, especially when the meshsize is small.

The EP-method gives much smaller velocity and pressure errors than the CR-method since it is a pressure robust method. For a given meshsize $h$, we note that the Post-method allows to reduce the velocity error by a factor larger than 10. The pressure error of the Post-method using one iteration is close to that of the CR-method, while it is greatly improved using eight iterations.

Fig. 7 shows the discrete error values $\epsilon_1^v$ and $\epsilon_D^v$ plotted against the meshsize for $v = 10^{-6}$. Iterating several times yields again a significant decrease of the indicators.

Table 4 shows the average convergence rates for the velocity, $\tau_u$ and the pressure, $\tau_p$. In the case $v = 1$, the average convergence rates are as expected. In the case $v = 10^{-6}$, the average convergence rate for the velocity, $\tau_u$, is better than expected for the EP-method and the Post-method, and the average convergence rate for the pressure, $\tau_p$, is better than expected for the Post-method, probably because the asymptotic convergence regime has not been reached.
<figure><img src="image_14.png" /><figcaption>Fig. 5. Sinusoidal velocity, $v = 1$. Plots of $\epsilon_0^\nu(\mathbf{u}_h)$ (left) and $\epsilon_0^\nu(p_h)$ (right).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 6. Sinusoidal velocity, $v = 10^{-6}$. Plots of $\epsilon_0^\nu(\mathbf{u}_h)$ (left) and $\epsilon_0^\nu(p_h)$ (right).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 7. Sinusoidal velocity, $v = 10^{-6}$. Plots of $\epsilon_1^\nu$ (left) and $\epsilon_D^\nu$ (right).</figcaption></figure>

Table 5 shows the CPU times for the CR-method and the Post-method with either one, or eight, iterations. As we have used the same meshes, the computation times are similar to those in Table 3, showing once again that the overhead cost decreases dramatically as the mesh size decreases.

Fig. 8 shows $\epsilon_0^\nu(\mathbf{u}_h)$ (right) and $\epsilon_0^\nu(p_h)$ (left) plotted against the CPU time for the CR-method and the Post-method with a single iteration (“Post” plot) or eight iterations (“Post-8” plot). In the same way as for the first test, we note that to reach $\epsilon_0^\nu(\mathbf{u}_h) \lesssim 5 \times 10^{-6}$, the required CPU time is $0.2$ s for the Post-method with one iteration and less than $0.02$ s for the Post-method with eight iterations, to be compared with the CR-method which does not reach this threshold.
<table><caption>**Table 4**
Sinusoidal velocity. Average convergence rates.</caption><thead><tr><th>ν</th><th>τ</th><th>CR</th><th>EP</th><th>Post</th><th>Post-8</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>τ<sub>u</sub></td><td>1.99</td><td>2.07</td><td>2.06</td><td>—</td></tr><tr><td>τ<sub>p</sub></td><td>1.03</td><td>1.06</td><td>1.11</td><td>—</td></tr><tr><td rowspan="2">10<sup>-6</sup></td><td>τ<sub>u</sub></td><td>2.05</td><td>2.17</td><td>2.47</td><td>2.36</td></tr><tr><td>τ<sub>p</sub></td><td>1.31</td><td>1.07</td><td>1.47</td><td>1.74</td></tr></tbody></table>

<table><caption>**Table 5**
Sinusoidal velocity, $v = 10^{-6}$. CPU time (s).</caption><thead><tr><th>h</th><th>CPU CR</th><th>CPU Post</th><th>overhead</th><th>CPU Post-8</th><th>overhead</th></tr></thead><tbody><tr><td>1.00 × 10<sup>-1</sup></td><td>3.85 × 10<sup>-3</sup></td><td>4.64 × 10<sup>-3</sup></td><td>21%</td><td>7.06 × 10<sup>-3</sup></td><td>83%</td></tr><tr><td>5.00 × 10<sup>-2</sup></td><td>2.22 × 10<sup>-2</sup></td><td>2.51 × 10<sup>-2</sup></td><td>13%</td><td>4.27 × 10<sup>-2</sup></td><td>92%</td></tr><tr><td>2.50 × 10<sup>-2</sup></td><td>2.43 × 10<sup>-1</sup></td><td>2.58 × 10<sup>-1</sup></td><td>6.2%</td><td>3.45 × 10<sup>-1</sup></td><td>42%</td></tr><tr><td>1.25 × 10<sup>-2</sup></td><td>4.27 × 10<sup>+0</sup></td><td>4.34 × 10<sup>+0</sup></td><td>1.6%</td><td>4.78 × 10<sup>+0</sup></td><td>12%</td></tr><tr><td>6.25 × 10<sup>-3</sup></td><td>4.53 × 10<sup>+1</sup></td><td>4.57 × 10<sup>+1</sup></td><td>0.7%</td><td>4.82 × 10<sup>+1</sup></td><td>6.4%</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 8.** Sinusoidal velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(u_h)$ (left) and $\epsilon_0^v(p_h)$ (right) against CPU time (s).</figcaption></figure>

## Comment on the algorithm

At this stage, we noticed that eight iterations was a good compromise. In order to design an optimized algorithm, another stopping criterion could be set by comparing two successive computations of $p_h^{old}$, i.e. comparing $|\delta P^n|_{\Omega_0}$ to $|\delta P^{n-1}|_{\Omega_0}$. The choice of $\lambda$ could also be further optimized, depending on $v$ and on the meshsize.

## Some observations

The coercive $P^1 \times P^0$ formulation (EP-method) gives pressure robust results, the obvious limitation being that it requires to know explicitly the potential of the gradient part of the source term. If it is not known, the two step method greatly reduces the velocity error, compared with the calculation carried out using the Crouzeix-Raviart $P^1$ formulation (CR-method). Moreover, using the second (post-processing) step iteratively greatly improves the initial result. Finally, the reduction factor is greater the smaller $v$ is.

## 7.4. Low regularity manufactured solution: Figs. 9-11, Tables 6-7

Last, we consider Problem (3.1)$_{NH}$ with a low regularity solution. Let $(\rho, \theta)$ be the polar coordinates centered in $(0.5, 0.5)$. Let $\alpha = 0.45$. We set $\mathbf{f} = -v\Delta\mathbf{u} + \text{grad} p$ where $(\mathbf{u}, p) = (\rho^\alpha e_\theta, \rho - \int_\Omega \rho)$. Results are given for $v = 10^{-6}$. The number of unknowns are reported on Table 1 (right). We used a refined mesh around $(0.5, 0.5)$, where the solution is of low regularity.

Fig. 9 shows the discrete error values $\epsilon_0^v(u_h)$ (left) and $\epsilon_0^v(p_h)$ (right) plotted against the meshsize. We remark that EP-method shows far better results than the CR-method, and that the Post-method allows again to improve the approximation of the CR-method.

Fig. 10 shows the discrete error values $\epsilon_1^v$ (left) and $\epsilon_D^v$ (right) plotted against the meshsize. We have $\epsilon_D^v/\epsilon_1^v \approx 0.5$ for the Post-method with a single iteration while $\epsilon_D^v/\epsilon_1^v \approx 0.1$ or the Post-method with eight iterations.

Table 6 shows the averaged convergence rates between the successive meshes. For the EP-method and the CR-method, we postulate that, asymptotically, $\tau_u = 1 + \alpha$ and $\tau_p = \alpha$. In both cases, the obtained convergence rates are better than expected, which suggest that the asymptotic convergence regime is not reached.
<figure><img src="image_18.png" /><figcaption>Fig. 9. Low regularity velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right).</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 10. Low regularity velocity, $v = 10^{-6}$. Plots of $\epsilon_1^v$ (left) and $\epsilon_p^v$ (right).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Table 6</figcaption></figure>

Low regularity velocity, $v = 10^{-6}$. Averaged convergence rates.

<table><caption>Table 7</caption><thead><tr><th>τ</th><th>CR</th><th>EP</th><th>Post-8</th></tr></thead><tbody><tr><td>τ<sub>u</sub></td><td>2.0</td><td>2.2</td><td>2.8</td></tr><tr><td>τ<sub>p</sub></td><td>1.6</td><td>0.7</td><td>1.9</td></tr></tbody></table>

Low regularity velocity, $v = 10^{-6}$. CPU time (s).

<table><thead><tr><th>h</th><th>CPU CR</th><th>CPU Post-8</th><th>overhead</th></tr></thead><tbody><tr><td>1.00 × 10<sup>-1</sup></td><td>1.12 × 10<sup>-2</sup></td><td>1.92 × 10<sup>-2</sup></td><td>71%</td></tr><tr><td>5.00 × 10<sup>-2</sup></td><td>2.77 × 10<sup>-1</sup></td><td>3.31 × 10<sup>-1</sup></td><td>19%</td></tr><tr><td>2.50 × 10<sup>-2</sup></td><td>3.15 × 10<sup>+0</sup></td><td>3.40 × 10<sup>+0</sup></td><td>7.9%</td></tr><tr><td>1.25 × 10<sup>-2</sup></td><td>1.06 × 10<sup>+1</sup></td><td>1.20 × 10<sup>+1</sup></td><td>13%</td></tr><tr><td>6.25 × 10<sup>-3</sup></td><td>1.32 × 10<sup>+2</sup></td><td>1.39 × 10<sup>+2</sup></td><td>5%</td></tr></tbody></table>

Table 7 shows the CPU times for the CR-method and the Post-method with eight iterations. Again, the overhead cost decreases sharply with the meshsize.

Fig. 11 shows the discrete error values $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right) against the CPU time for the CR-method and the Post-method with eight iterations. To reach $\epsilon_0^v(\mathbf{u}_h) \lesssim 10^{-5}$, the required CPU time is 0.02 s for the Post-method with eight iterations, to be compared with more than 100 s for the CR-method.
<figure><img src="image_22.png" /><figcaption>Fig. 11. Low regularity velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right) against CPU time (s).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 12. Linear velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 13. Sinusoidal velocity, $v = 10^{-6}$. Plots of $\epsilon_0^v(\mathbf{u}_h)$ (left) and $\epsilon_0^v(p_h)$ (right).</figcaption></figure>

## 7.5. Numerical results using Raviart-Thomas projection

The use of Crouzeix-Raviart $P^1_{nc} \times P^0$ formulation with $H(\text{div})$-conforming projection of the test function on the right-hand side leads to a pressure robust discrete velocity [46,20]. In this case too, we can use the Post-method to obtain a precise $H^1$-conforming approximation of the discrete velocity. On Figs. 12-14, we represent $\epsilon_0^v(\mathbf{u}_h)$ and $\epsilon_0^v(p_h)$ against the meshsize for our three tests. For the CR-method, we use the lowest order Raviart-Thomas projection [46, §3.2] to compute $p_h^{old}$, which is then used in the Post-method with $\lambda = 1$. For the first test (linear velocity), we obtain $\epsilon_0^v(\mathbf{u}_h) \lesssim 10^{-13}$ and $\epsilon_0^v(p_h) \lesssim 10^{-12}$, and regardless, we observe that the errors (close to machine precision for the CR-method) do not deteriorate when the Post-method is applied. The pressure error $\epsilon_0^v(p_h)$ remains unchanged for the other two tests. The approximation of the velocity is improved for all three tests.
<figure><img src="image_25.png" /><figcaption>Fig. 14. Low regularity velocity, $v = 10^{-6}$. Plots of $\epsilon_0^+(\mathbf{u}_h)$ (left) and $\epsilon_0^-(p_h)$ (right).</figcaption></figure>

## 8. Conclusion

We proposed and analysed a new variational formulation of the Stokes problem based on T-coercivity theory. This variational formulation is coercive, and can be discretized with the Scott-Vogelius $P^k \times P^{k-1}_{disc}$ finite element for all $k \ge 1$. To solve the linear system resulting from the discretization, we need to know the pressure, or at least some approximation of it, which in our numerical tests is the discrete pressure obtained using the classical nonconforming method with the Crouzeix-Raviart $P^1_{nc} \times P^0$ finite element. This two step method improves the numerical results by notably reducing the errors obtained after the use of the classical method, especially when the viscosity is small. More significantly, the two step method consistently outperforms the classical method in terms of precision with respect to CPU time.

## Acknowledgements

Erell Jamelot would like to thank CEA SIMU/SITHY project.

The authors would like to thank the reviewers for their careful assessment.

## Appendix A. New variational formulation without orthogonality

Let us briefly consider what would happen if we were to consider a nonlocal right inverse of the divergence operator, that is different from the one proposed in Proposition 1. For example, with values in the whole of $H_0^1(\Omega)$ (not restricted to $V^\perp$). Such an operator $q \mapsto \bar{v}_q$ is considered in Section 3.1 of [19] (2D domain with a smooth boundary) and in [4] (2D convex domain). For this operator, let $\bar{C}_{div}$ be a constant such that:

$$ \forall q \in L^2_{zmv}(\Omega), \exists \bar{v}_q \in H_0^1(\Omega) \mid \operatorname{div} \bar{v}_q = q \text{ and } \|\bar{v}_q\|_{H_0^1(\Omega)} \le \bar{C}_{div} \|q\|_{L^2(\Omega)}. \quad (A.1) $$

For both operators, we emphasize that $\bar{v}_q$ does not belong to $V^\perp$ in general.

Let $\lambda > \frac{1}{4}(\bar{C}_{div})^2$ (cf. Remark 1). T-coercivity can be obtained as before, with the operator

$$ \left\{ \begin{array}{l} \bar{T} : \lambda \to \mathbb{R} \\ (\mathbf{v}, q) \mapsto (\lambda \mathbf{v} - v^{-1}\bar{v}_q, -\lambda q) \end{array} \right. \quad (A.2) $$

Define $\bar{a}_\lambda((\mathbf{u}', p'), (\mathbf{v}, q)) = a((\mathbf{u}', p'), \bar{T}(\mathbf{v}, q))$. We have:

$$ \begin{aligned} \bar{a}_{\lambda}((\mathbf{u}', p'), (\mathbf{v}, q)) ={}& v \lambda(\mathbf{u}', \mathbf{v})_{H_0^1(\Omega)} - (\mathbf{u}', \bar{v}_q)_{H_0^1(\Omega)} \\ & - \lambda(p', \operatorname{div} \mathbf{v})_{L^2(\Omega)} + v^{-1}(p', q)_{L^2(\Omega)} \\ & + \lambda(q, \operatorname{div} \mathbf{u}')_{L^2(\Omega)}. \end{aligned} \quad (A.3) $$

However, for $\mathbf{u}' \in V$ (that is for solutions $(\mathbf{u}', p)$ to Problem (4.3)), the term $(\mathbf{u}', \bar{v}_q)_{H_0^1(\Omega)}$ can no longer be removed from the expression (4.2) of the bilinear form $\bar{a}_\lambda$ because orthogonality is lost. As a matter of fact one has

$$ (\mathbf{u}', \bar{v}_q)_{H_0^1(\Omega)} = (\operatorname{curl} \mathbf{u}', \operatorname{curl} \bar{v}_q)_{L^2(\Omega)}. $$

Hence, in the non-orthogonal case, one has to choose an ad hoc constant $\lambda > \frac{1}{4}(\bar{C}_{div})^2$ to ensure T-coercivity (cf. Proposition 2). Similarly in the expression of the right-hand side $\bar{\mathcal{E}}_\lambda$.
**Proposition 8.** Let $f' \in H^{-1}(\Omega)$ be decomposed as in (4.5). Given, $q \in L^2_{zmv}(\Omega)$, let $\bar{v}_q \in H_0^1(\Omega)$ be defined by (A.1). We have:

$$ \langle \mathbf{f}', \bar{v}_q \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = -(z_{f'}, q)_{L^2(\Omega)} + (\operatorname{curl} \mathbf{w}_{f'}, \operatorname{curl} \bar{v}_q)_{L^2(\Omega)}. \quad (A.4) $$

To conclude on the use of explicit T-coercivity in the non-orthogonal case, we observe that the problem, if split as in (4.9), leads to a more intricate variational formulation, which reads

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ s.t. for all } (\mathbf{v}, q) \in \mathcal{X} \\ \text{(i)'} \quad \nu \lambda(\mathbf{u}, \mathbf{v})_{H_0^1(\Omega)} - \lambda(p, \operatorname{div} \mathbf{v})_{L^2(\Omega)} = \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)}, \\ \text{(ii)'} \quad -(\mathbf{u}, \bar{v}_q)_{H_0^1(\Omega)} + \lambda(q, \operatorname{div} \mathbf{u})_{L^2(\Omega)} + \nu^{-1}(p, q)_{L^2(\Omega)} \\ \qquad = \nu^{-1}(z_{f'}, q)_{L^2(\Omega)} - \nu^{-1}(\operatorname{curl} \mathbf{w}_{f'}, \operatorname{curl} \bar{v}_q)_{L^2(\Omega)}. \end{array} \right. $$

Regarding the right-hand side in (ii'), it requires the knowledge of both $z_f$ and of $w_f$ (or of suitable approximations): a two-step procedure could be used as before (see the end of §4.2). One needs to evaluate also $(\mathbf{u}, \bar{v}_q)_{H_0^1(\Omega)} = (\operatorname{curl} \mathbf{u}, \operatorname{curl} \bar{v}_q)_{L^2(\Omega)}$ in the left-hand side of (ii'). The latter part requires some knowledge of the nonlocal right inverse of the divergence operator [4].

## Appendix B. Nonhomogeneous Dirichlet boundary conditions

We now solve the classical incompressible Stokes model, with non zero Dirichlet boundary conditions. Hence, we focus on Problem (3.1)$_{NH}$ with $g = 0$. The variational formulation reads:

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{u}, p) \in H^1(\Omega) \times L^2_{zmv}(\Omega) \text{ such that} \\ \nu(\mathbf{u}, \mathbf{v})_{H_0^1(\Omega)} - (p, \operatorname{div} \mathbf{v})_{L^2(\Omega)} = \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} \quad \forall \mathbf{v} \in H_0^1(\Omega), \\ (q, \operatorname{div} \mathbf{u})_{L^2(\Omega)} = 0 \quad \forall q \in L^2_{zmv}(\Omega), \\ \mathbf{u} = \mathbf{g} \text{ on } \partial\Omega, \end{array} \right. \quad (B.1) $$

where $g \in H_0^1(\partial\Omega)$ is some boundary data such that $g \cdot n \in L^2_{zmv}(\partial\Omega)$. Let $\mathbf{u}_g \in H^1(\Omega)$ be such that $\mathbf{u}_g = \mathbf{g}$ on $\partial\Omega$ and $-\Delta\mathbf{u}_g = 0$ in $\Omega$. Then consider $(\mathbf{u}_0, p) \in H_0^1(\Omega) \times L^2_{zmv}(\Omega)$ the unique solution to (3.1)$_H$ with data $(\mathbf{f} + \nu\Delta\mathbf{u}_g, -\operatorname{div}\mathbf{u}_g)$: Find $(\mathbf{u}_0, p) \in \mathcal{X}$ such that for all $(\mathbf{v}, q) \in \mathcal{X}$

$$ a((\mathbf{u}_0, p), (\mathbf{v}, q)) = \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} + (q, \operatorname{div} \mathbf{u}_g)_{L^2(\Omega)}, \quad (B.2) $$

where we use that $(\mathbf{u}_g, \mathbf{v})_{H_0^1(\Omega)} = 0$. Defining $\mathbf{u} = \mathbf{u}_g + \mathbf{u}_0 \in H^1(\Omega)$, we find that $(\mathbf{u}, p)$ is solution to Problem (B.1). Uniqueness and continuous dependence with respect to the data are easily obtained. The next step is to replace $(\mathbf{v}, q)$ by $T((\mathbf{v}, q)) = (\lambda\mathbf{v} - \nu^{-1}\mathbf{v}_q, -\lambda q)$. One finds by orthogonality that for all $(\mathbf{v}, q) \in \mathcal{X}$:

$$ \begin{aligned} \nu \lambda(\mathbf{u}_0, \mathbf{v})_{H_0^1(\Omega)} - \lambda(p, \operatorname{div} \mathbf{v})_{L^2(\Omega)} + \nu^{-1}(p, q)_{L^2(\Omega)} &+ \lambda(q, \operatorname{div} \mathbf{u}_0)_{L^2(\Omega)} \\ & = \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} - \nu^{-1} \langle \mathbf{f}, \mathbf{v}_q \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} - \lambda(q, \operatorname{div} \mathbf{u}_g)_{L^2(\Omega)}. \end{aligned} $$

Hence, a new variational formulation for nonhomogeneous boundary Dirichlet conditions is: Find $(\mathbf{u}, p) \in H^1(\Omega) \times L^2_{zmv}(\Omega)$ such that for all $(\mathbf{v}, q) \in \mathcal{X}$

$$ \begin{aligned} \nu \lambda(\mathbf{u}, \mathbf{v})_{H_0^1(\Omega)} - \lambda(p, \operatorname{div} \mathbf{v})_{L^2(\Omega)} + \nu^{-1}(p, q)_{L^2(\Omega)} &+ \lambda(q, \operatorname{div} \mathbf{u})_{L^2(\Omega)} \\ & = \lambda \langle \mathbf{f}, \mathbf{v} \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} + \nu^{-1}(z_{f'}, q)_{L^2(\Omega)}. \end{aligned} \quad (B.3) $$

This is the same variational formulation as the one with homogeneous Dirichlet boundary conditions, except that $\mathbf{u} \in H^1(\Omega)$ and $\mathbf{u}_{|\partial\Omega} = \mathbf{g}$.

Regarding the numerical algorithms, if $z_f$ is known, one solves a linear system like

$$ \left\{ \begin{array}{l} \text{Find } (\underline{U}, \underline{P}) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\ \nu \lambda \underline{A} \underline{U} - \lambda \underline{B}^T \underline{P} = \lambda \underline{F} - \nu \lambda \underline{A} \underline{U}_{NH} \\ \lambda \underline{B} \underline{U} + \nu^{-1} \underline{M} \underline{P} = \nu^{-1} \underline{M} \underline{Z} - \lambda \underline{B} \underline{U}_{NH} \end{array} \right. \quad (B.4) $$

where, in the right-hand side, $\underline{U}_{NH} \in \mathbb{R}^{N_u}$ accounts for $\mathbf{u}_g$. This is completely similar to (5.1). While, if $z_f$ is not known, starting from an initial guess $\underline{P}^{-1} \in \mathbb{R}^{N_p}$, for $n = 0, 1, \dots$, one solves linear systems like

$$ \left\{ \begin{array}{l} \text{Find } (\underline{U}^n, \underline{P}^n) \in \mathbb{R}^{N_u} \times \mathbb{R}^{N_p} \text{ such that:} \\ \nu \lambda \underline{A} \underline{U}^n - \lambda \underline{B}^T \underline{P}^n = \lambda \underline{F}^n - \nu \lambda \underline{A} \underline{U}^N \\ \lambda \underline{B} \underline{U}^n + \nu^{-1} \underline{M} \underline{P}^n = \nu^{-1} \underline{M} \underline{P}^{n-1} - \lambda \underline{B} \underline{U}^N \end{array} \right. . \quad (B.5) $$
Interestingly, one recovers the same results as those of Theorem 4, because one finds identical iterating matrices, cf. (5.4) and (5.5). Let $\tilde{\Pi}_{h,cg}$ be the $L^2(\Omega)$-orthogonal projection operator from $L^2(\Omega)$ to $V_h^{k,cg}$. The (conforming) discretization of Problem (B.3) with $\mathbf{P}^k \times P_{disc}^{k-1}$ finite element reads:

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{u}_h, p_h) \in \mathbf{V}_h^{k,cg} \times \mathbf{Q}_h^{k-1} \text{ s.t. for all } (\mathbf{v}_h, q_h) \in \mathbf{V}_{0,h}^{k,cg} \times \mathbf{Q}_h^{k-1} \\ \qquad \nu \lambda(\mathbf{u}_h, \mathbf{v}_h) \mathbf{H}_0^{H^{-1}(\Omega)} - \lambda(p_h, \text{div } \mathbf{v}_h) L^2(\Omega) = \lambda(\mathbf{f}, \mathbf{v}_h) \mathbf{H}^{-1}(\Omega) \mathbf{H}_0^*(\Omega), \\ \qquad \lambda(q_h, \text{div } \mathbf{u}_h) L^2(\Omega) + \nu^{-1}(p_h, q_h) L^2(\Omega) = \nu^{-1}(z_f, q_h) L^2(\Omega), \\ \mathbf{u}_h|_{\partial\Omega} = (\tilde{\Pi}_{h,cg} \mathbf{u}_g)|_{\partial\Omega}. \end{array} \right. \qquad (B.6) $$

For nonhomogeneous boundary Dirichlet conditions, one can study the error estimates by introducing $\mathbf{u}_{0,h} = \mathbf{u}_h - \tilde{\Pi}_{h,cg} \mathbf{u}_g \in \mathbf{V}_{0,h}^{k,cg}$ and using the results that have been obtained for the homogeneous boundary conditions. First, in the estimate (6.7) of Proposition 7, we change the term $\|z_f - p_h^{old}\|_{L^2(\Omega)}$ in the right-hand side of the two equations into $\|z_f - p_h^{old}\|_{L^2(\Omega)} + \lambda\nu\|\text{div}(\mathbf{u}_g - \tilde{\Pi}_{h,cg} \mathbf{u}_g)\|_{L^2(\Omega)}$. Second, in the estimate (6.9) of Theorem 5, we add $\sqrt{\lambda}\|\text{div}(\mathbf{u}_g - \tilde{\Pi}_{h,cg} \mathbf{u}_g)\|_{L^2(\Omega)}$ to the right-hand side.

## Data availability

Data will be made available on request.

## References

[1] M. Ainsworth, A. Allendes, G. Barrenechea, R. Rankin, On the adaptive selection of the parameter in stabilized finite element approximations, SIAM J. Numer. Anal. 51 (2013) 1585–1609.

[2] A. Allendes, G. Barrenechea, C. Naranjo, A divergence-free low-order stabilized finite element method for a generalized steady state Boussinesq problem, Comput. Methods Appl. Mech. Eng. 340 (2018) 90–120.

[3] C. Amrouche, P.G. Ciarlet, C. Mardare, On a lemma of Jacques-Louis Lions and its relation to other fundamental results, J. Math. Pures Appl. 104 (2015) 207–226.

[4] R. Araya, G. Gatica, A new nonconforming Galerkin scheme for the Stokes problem: partially circumventing the discrete Babuska-Brezzi condition, Comput. Methods Appl. Mech. Eng. 107 (1993) 193–208.

[5] D. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337–344.

[6] M. Barré, P. Ciarlet Jr, The T-coercivity approach for mixed problems, C. R. Math. 362 (2024) 1051–1088.

[7] G. Barrenechea, E. Burman, E. Cáceres, J. Guzmán, Continuous interior penalty stabilization for divergence-free finite element methods, IMA J. Numer. Anal. 44 (2024) 980–1002.

[8] G. Barrenechea, F. Valentin, An unusual stabilized finite element method for a generalized Stokes problem, Numer. Math. 92 (2002) 653–677.

[9] G. Barrenechea, F. Valentin, Consistent local projection stabilized finite element methods, SIAM J. Numer. Anal. 48 (2010) 1801–1825.

[10] C. Bernardi, G. Raugel, Méthodes d'éléments finis mixtes pour les équations de Stokes et de Navier-Stokes dans un polygone non convexe, Calcolo 18 (1981) 255–291.

[11] P. Bochev, C. Dohrmann, M. Gunzburger, Stabilization of low-order mixed finite elements for the Stokes equations, SIAM J. Numer. Anal. 44 (2006) 82, 20.

[12] D. Boffi, F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods and Applications, Springer-Verlag, 2013.

[13] A.S. Bonnet-Ben Dhia, C. Carvalho, P. Ciarlet Jr, Mesh requirements for the finite element approximation of problems with sign-changing coefficients, Numer. Math. 138 (2018) 801–838.

[14] A.S. Bonnet-Ben Dhia, L. Chesnel, P. Ciarlet Jr, T-coercivity for scalar interface problems between dielectrics and metamaterials, Math. Modell. Numer. Anal. 46 (2012) 1363–1387.

[15] A.S. Bonnet-Ben Dhia, L. Chesnel, P. Ciarlet Jr, T-coercivity for the Maxwell problem with sign-changing coefficients, Commun. Partial Differ. Equ. 39 (2014) 1007–1031.

[16] A.S. Bonnet-Ben Dhia, L. Chesnel, H. Haddar, On the use of T-coercivity to study the interior transmission eigenvalue problem, C. R. Math. Acad. Sci. Paris 349 (2011) 647–651.

[17] A.S. Bonnet-Ben Dhia, P. Ciarlet Jr, Variational methods for the analysis of noncoercive problems, M.Sc. AMS Lecture Notes, Institut Polytechnique de Paris, France, 2024 (in French), https://hal.science/hal-03764649/.

[18] A.S. Bonnet-Ben Dhia, P. Ciarlet Jr, C. Zwölf, Time harmonic wave diffraction problems in materials with sign-shifting coefficients, J. Comput. Appl. Math. 234 (2010) 1912–1919, corrigendum p. 2616.

[19] F. Boyer, P. Fabrie, Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models, Springer-Verlag, 2013.

[20] C. Brennecke, A. Linke, C. Merdon, J. Schöberl, Optimal and pressure-independent $L^2$ velocity error estimates for a modified Crouzeix-Raviart Stokes element with BDM reconstructions, J. Comput. Math. 33 (2015) 191–208.

[21] A. Buffa, Remarks on the discretization of some noncoercive operator with applications to heterogeneous Maxwell equations, SIAM J. Numer. Anal. 43 (2005) 1–18.

[22] A. Buffa, S. Christiansen, The electric field integral equation on Lipschitz screens: definitions and numerical approximation, Numer. Math. 94 (2003) 229–267.

[23] A. Buffa, M. Costabel, C. Schwab, Boundary element methods for Maxwell's equations on non-smooth domains, Numer. Math. 92 (2002) 679–710.

[24] N. Chaabane, V. Girault, B. Riviere, T. Thompson, A stable enriched Galerkin element for the Stokes problem, Appl. Numer. Math. 132 (2018) 1–21.

[25] L. Chesnel, P. Ciarlet Jr, T-coercivity and continuous Galerkin methods: application to transmission problems with sign changing coefficients, Numer. Math. 124 (2013) 1–29.

[26] P. Ciarlet Jr, T-coercivity: application to the discretization of Helmholtz-like problems, Comput. Math. Appl. 64 (2012) 22–24.

[27] P. Ciarlet Jr, E. Jamelot, Variational methods for solving numerically magnetostatic systems, Adv. Comput. Math. 50 (2024), https://doi.org/10.1007/s10444-023-10089-1.

[28] P. Ciarlet Jr, E. Jamelot, F. Kpadonou, Domain decomposition methods for the diffusion equation with low-regularity solution, Comput. Math. Appl. 74 (2017) 2369–2384.

[29] M. Crouzeix, P.A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations, RAIRO. Anal. Numér. 7 (1973) 33–75.

[30] A. Ern, J.L. Guermond, Finite Elements I, Texts in Applied Mathematics, vol. 72, Springer, 2021.

[31] A. Ern, J.L. Guermond, Finite Elements II, Texts in Applied Mathematics, vol. 73, Springer, 2021.
[32] M. Fabien, J. Guzmán, M. Neilan, A. Zytoon, Low-order divergence-free approximations for the Stokes problem on Worsey-Farin and Powell-Sabin splits, Comput. Methods Appl. Mech. Eng. 390 (2022) 11444.

[33] V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer-Verlag, 1986.

[34] L. Giret, Non-Conforming Domain Decomposition for the Multigroup Neutron SPN Equation, Ph.D. thesis, Université Paris-Saclay, 2018.

[35] J.L. Guermond, P.D. Minev, High-order time stepping for the incompressible Navier-Stokes equations, SIAM J. Sci. Comput. 37 (2015) A2656–A2681.

[36] J.L. Guermond, P.D. Minev, High-order time stepping for the Navier-Stokes equations with minimal computational complexity, J. Comput. Appl. Math. 310 (2017) 92–103.

[37] J. Guzmán, L.R. Scott, The Scott-Vogelius finite elements revisited, Math. Comput. 88 (2019) 515–529.

[38] M. Halla, Galerkin approximation of holomorphic eigenvalue problems: weak T-coercivity and T-compatibility, Numer. Math. 148 (2021) 387–407.

[39] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface, Phys. Fluids 8 (1965) 2182–2189.

[40] R. Hiptmair, Finite elements in computational electromagnetics, Acta Numer. (2002) 237–339.

[41] T. Hohage, L. Nannen, Convergence of infinite element methods for scalar waveguide problems, BIT Numer. Math. 55 (2015) 215–254.

[42] E. Jamelot, Stokes_ExplicitTC code, https://github.com/cea-trust-platform/Stokes_ExplicitTC, 2024.

[43] E. Jamelot, Stability estimates for solving Stokes problem with nonconforming finite elements, C. R. Math. 363 (2025) 115–137.

[44] E. Jamelot, P. Ciarlet Jr, Fast non-overlapping Schwarz domain decomposition methods for solving the neutron diffusion equation, J. Comput. Phys. 241 (2013) 445–463.

[45] Y. Li, L.T. Zikatanov, New stabilized $P_1 \times P_0$ finite element methods for nearly inviscid and incompressible flows, Comput. Methods Appl. Mech. Eng. 393 (2022) 114815.

[46] A. Linke, On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime, Comput. Methods Appl. Mech. Eng. 278 (2014) 782–800.

[47] S. Nicaise, J. Venel, A posteriori error estimates for a finite element approximation of transmission problems with sign changing coefficients, J. Comput. Appl. Math. 235 (2011) 4272–4282.

[48] M.A. Olshanskii, E.E. Tyrtyshnikov, Iterative Methods for Linear Systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2014.

[49] L.R. Scott, M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in: Lectures in Applied Mathematics, Part 2, Amer. Math. Soc., 1985, pp. 221–244.

[50] J. Shen, On error estimates of the penalty method for unsteady Navier-Stokes equations, SIAM J. Numer. Anal. 32 (1995) 386–403.

[51] C. Taylor, T. Hood, Numerical solution of the Navier-Stokes equations using the finite element technique, Comput. Fluids 1 (1973) 73–100.

[52] G. Unger, Convergence analysis of a Galerkin boundary element method for electromagnetic resonance problems, Ser. Partial Differ. Equ. Appl. 39 (2021).

[53] S. Zhang, A new family of stable mixed finite elements for the 3D Stokes equations, Math. Comput. 74 (2005) 543–554.

[54] T. Zhang, L. Tang, A stabilized finite volume method for Stokes equations using the lowest order $P_1 - P_0$ element pair, Adv. Comput. Math. 41 (2015) 781–798.