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

HAL
open science

# A curl preserving finite volume scheme by space velocity enrichment. Application to the low Mach number accuracy problem

Jonathan Jung, Vincent Perrier

► To cite this version:

Jonathan Jung, Vincent Perrier. A curl preserving finite volume scheme by space velocity enrichment. Application to the low Mach number accuracy problem. 2024. hal-04563838

HAL Id: hal-04563838

https://inria.hal.science/hal-04563838

Preprint submitted on 30 Apr 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A curl preserving finite volume scheme by space velocity enrichment. Application to the low Mach number accuracy problem.

Jonathan Jung\* and Vincent Perrier†

July 18, 2024

## Abstract

In this article, we address the problem of accuracy of finite volume schemes in the low Mach number limit. It has been known for years that collocated finite volume schemes are naturally correctly behaving in this limit on triangular meshes [21, 22, 16], but fail in general on other types of mesh. We are first interested in the general problem of the conservation of vorticity for the wave system. By enriching the approximation space for vectors, we prove that the Hodge-Helmholtz context developed for triangular meshes in [16] can be recovered in the quadrangular mesh case. This leads to a numerical scheme for the wave system that naturally preserves the vorticity under mild assumption on the numerical flux. The new approximation space is then used with the barotropic Euler system. Numerical tests show that the new numerical scheme is accurate for both steady and acoustic problems at low Mach number.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Cartesian and periodic case</td><td>3</td></tr><tr><td>2.1</td><td>A new finite element space for vectorial field</td><td>4</td></tr><tr><td>2.2</td><td>Approximation space on Cartesian mesh</td><td>4</td></tr><tr><td>2.3</td><td>Discrete Hodge-Helmholtz decomposition</td><td>5</td></tr><tr><td>2.4</td><td>First order wave system</td><td>8</td></tr><tr><td>2.4.1</td><td>Discretization</td><td>8</td></tr><tr><td>2.4.2</td><td>Godunov stabilization and stationary solution</td><td>9</td></tr><tr><td>2.4.3</td><td>Conservation of the adjoint curl with the Godunov' flux</td><td>10</td></tr><tr><td>3</td><td>Quadrangular case with boundary conditions</td><td>11</td></tr><tr><td>3.1</td><td>Discretization and boundary conditions</td><td>11</td></tr><tr><td>3.2</td><td>Discrete Hodge-Helmholtz decomposition</td><td>12</td></tr><tr><td>3.2.1</td><td>The two Piola transformations and the choice of the approximation space</td><td>12</td></tr><tr><td>3.2.2</td><td>Determination of u<sub>h</sub><sup>φ</sup></td><td>13</td></tr><tr><td>3.2.3</td><td>Existence and uniqueness of the decomposition</td><td>14</td></tr><tr><td>3.3</td><td>Structure preserved and long time behaviour with Godunov' flux</td><td>17</td></tr><tr><td>4</td><td>Discretization for the barotropic Euler system</td><td>19</td></tr><tr><td>5</td><td>Numerical results</td><td>19</td></tr><tr><td>5.1</td><td>Order of accuracy on the adjoint curl</td><td>20</td></tr><tr><td>5.2</td><td>Wave equation</td><td>20</td></tr><tr><td>5.2.1</td><td>Periodic vortex</td><td>20</td></tr><tr><td>5.2.2</td><td>Cylinder scattering</td><td>21</td></tr><tr><td>5.3</td><td>Euler equation</td><td>23</td></tr><tr><td>5.3.1</td><td>Cylinder scattering</td><td>23</td></tr><tr><td>5.3.2</td><td>Propagation of a low Mach number acoustic wave over a steady vortex</td><td>25</td></tr></table>

<footnote>\*LMA, UPPA, Pau, France and Cagire team, Inria Bordeaux Sud-Ouest, France jonathan.jung@univ-pau.fr.</footnote>

<footnote>†Cagire team, Inria Bordeaux Sud-Ouest, France and LMA, UPPA, Pau, France, vincent.perrier@inria.fr.</footnote>
# A curl preserving finite volume scheme by space velocity enrichment. Application to the low Mach number accuracy problem.

Jonathan Jung* and Vincent Perrier†

March 21, 2024

## Abstract

In this article, we address the problem of accuracy of finite volume schemes in the low Mach number limit. It has been known for years that collocated finite volume schemes are naturally correctly behaving in this limit on triangular meshes [20, 21, 15], but fail in general on other types of mesh. We are first interested in the general problem of the conservation of vorticity for the wave system. By enriching the approximation space for vectors, we prove that the Hodge-Helmholtz context developed for triangular meshes in [15] can be recovered in the quadrangular mesh case. This leads to a numerical scheme for the wave system that naturally preserves the vorticity under mild assumption on the numerical flux. The new approximation space is then used with the barotropic Euler system. Numerical tests show that the new numerical scheme is accurate for both steady and acoustic problems at low Mach number.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Cartesian and periodic case</td><td>3</td></tr><tr><td>2.1</td><td>A new finite element space for vectorial field</td><td>3</td></tr><tr><td>2.2</td><td>Approximation space on Cartesian mesh</td><td>4</td></tr><tr><td>2.3</td><td>Discrete Hodge-Helmholtz decomposition</td><td>5</td></tr><tr><td>2.4</td><td>First order wave system</td><td>8</td></tr><tr><td>2.4.1</td><td>Discretization</td><td>8</td></tr><tr><td>2.4.2</td><td>Godunov stabilization and stationary solution</td><td>9</td></tr><tr><td>2.4.3</td><td>Conservation of the adjoint curl with the Godunov' flux</td><td>10</td></tr><tr><td>3</td><td>Quadrangular case with boundary conditions</td><td>10</td></tr><tr><td>3.1</td><td>Discretization and boundary conditions</td><td>11</td></tr><tr><td>3.2</td><td>Discrete Hodge-Helmholtz decomposition</td><td>12</td></tr><tr><td>3.2.1</td><td>The two Piola transformations and the choice of the approximation space</td><td>12</td></tr><tr><td>3.2.2</td><td>Determination of u<sub>h</sub><sup>c</sup></td><td>13</td></tr><tr><td>3.2.3</td><td>Existence and uniqueness of the decomposition</td><td>14</td></tr><tr><td>3.3</td><td>Structure preserved and long time behaviour with Godunov' flux</td><td>17</td></tr><tr><td>4</td><td>Discretization for the barotropic Euler system</td><td>19</td></tr><tr><td>5</td><td>Numerical results</td><td>19</td></tr><tr><td>5.1</td><td>Order of accuracy on the adjoint curl</td><td>20</td></tr><tr><td>5.2</td><td>Wave equation</td><td>20</td></tr><tr><td>5.2.1</td><td>Periodic vortex</td><td>20</td></tr><tr><td>5.2.2</td><td>Cylinder scattering</td><td>21</td></tr><tr><td>5.3</td><td>Euler equation</td><td>22</td></tr><tr><td>5.3.1</td><td>Cylinder scattering</td><td>22</td></tr><tr><td>5.3.2</td><td>Propagation of a low Mach number acoustic wave over a steady vortex</td><td>23</td></tr></table>

<footnote>*LMA, UPPA, Pau, France and Cagire team, Inria Bordeaux Sud-Ouest, France jonathan.jung@univ-pau.fr.</footnote>

<footnote>†Cagire team, Inria Bordeaux Sud-Ouest, France and LMA, UPPA, Pau, France, vincent.perrier@inria.fr.</footnote>
# 1 Introduction

In this article, we are interested in two types of problems that are strongly interconnected: the accuracy problem at low Mach number of upwind finite volume numerical schemes for the compressible barotropic Euler system, and the conservation of the vorticity in the first order formulation of the wave system. The strong connection between these two problems was for example discussed in [27], and we mainly focus this state of the art section on the long time accuracy problem of the wave system:

$$ \left\{ \begin{aligned} \partial_{\tau} p + \frac{1}{\rho_0} \mathrm{div}_{\mathbf{x}} \mathbf{u} &= 0, \\ \partial_{\tau} \mathbf{u} + \kappa_0 \nabla_{\mathbf{x}} p &= 0, \end{aligned} \right. \qquad (1) $$

where $\mathbf{u}$ is the velocity and $p$ the pressure. System (1) depends on two strictly non-negative parameters, $\kappa_0$ and $\rho_0$. The wave velocity is $c_0$, linked with the parameters of the system by $c_0^2 = \kappa_0/\rho_0$. By taking the curl of the velocity equation of (1), we formally find $\partial_{\tau} (\nabla_{\mathbf{x}} \times \mathbf{u}) = 0$, which means that the vorticity $\nabla_{\mathbf{x}} \times \mathbf{u}$ is conserved. Conservation of this quantity is a necessary condition for having long time accuracy.

The problem of conservation of the vorticity for (1) has been tackled by several strategies, some of them being linked with the conservation of the divergence of a vector $\mathbf{v}$ when $\mathbf{v}$ ensures the following type of equation

$$ \partial_{\tau} \mathbf{v} + \nabla_{\mathbf{x}} \times \mathbf{g} = 0. \qquad (2) $$

These different strategies may be gathered into the three following families

1. **Projection method.** Projection methods have been widely used in the context of preservation of a divergence for incompressible flows [6] but also for magnetohydrodynamic system [10]. This method was rather designed for (2), and consists in computing a candidate update $\bar{\mathbf{v}}^{n+1}$ which is then corrected. For ensuring preservation of the divergence, a potential $\varphi$ is computed such that

$$ \Delta\varphi = \mathrm{div}_{\mathbf{x}} (\bar{\mathbf{v}}^{n+1} - \mathbf{v}^n), \qquad (3) $$

and the update is projected as $\mathbf{v}^{n+1} = \bar{\mathbf{v}}^{n+1} - \nabla_{\mathbf{x}}\varphi$, which leads to the conservation of the divergence of $\mathbf{v}$. As far as we know, it has never been proposed for the conservation of the curl for (1), but could be easily adapted. This method raises several problems: first, this requires to solve the elliptic problem (3) at each time step, which is costly, and second, (3) should be equipped with appropriate boundary conditions, which is not always straightforward.

2. **Discretizations based on staggered data.** This type of discretization consists in having data located on the faces or edges of the cells, whereas other data are located in the center of the cells. For example, the MAC scheme [30] consists in discretizing the pressure at the center of cells whereas the velocity components are normal to each face. The MAC scheme has been thoroughly analyzed in [35, 36], and extended to compressible flows in [18, 20, 19]. Similar approach were designed in [24, 15, 23, 32], in which general definitions of discrete gradient, divergence and curl were provided with a staggered design. Application of staggered schemes to hyperbolic problems can be found for example in [8]. A high order discontinuous Galerkin version of staggered schemes was proposed in [41], and high order finite volume schemes were proposed in [2]. Still linked with staggered ideas, the scheme [7] was proposed, for which the data are collocated, but the divergence is cleaned on a staggered grid. The main drawback of staggered data is its difficulty of implementation on unstructured meshes at high order.

3. **Generalized Lagrange multipliers.** The general Lagrange multiplier method was designed originally for the conservation of the divergence for the Maxwell system [34] and for the MHD system [14]. This method consists in adding an extra variable replacing the divergence to be preserved, and to define an extra equation on this additional variable, which includes a relaxation between the divergence of the vector and the additional variable. The extension to the curl constraint preservation was proposed in [17]. The main drawbacks of this method are the increase of the number of unknowns to solve, and the addition of numerical parameters, that should be tuned (propagation speed of the additional variable, stiffness of the relaxation).
# 1 Introduction

In this article, we are interested in two types of problems that are strongly interconnected: the accuracy problem at low Mach number of upwind finite volume numerical schemes for the compressible barotropic Euler system, and the conservation of the vorticity in the first order formulation of the wave system. The strong connection between these two problems was for example discussed in [25], and we mainly focus this state of the art section on the long time accuracy problem of the wave system:

$$ \left\{ \begin{aligned} \partial_{\tau} p + \frac{1}{\rho_0} \mathrm{div}_{\mathbf{x}} \mathbf{u} &= 0, \\ \partial_{\tau} \mathbf{u} + \kappa_0 \nabla_{\mathbf{x}} p &= 0, \end{aligned} \right. \qquad (1) $$

where $\mathbf{u}$ is the velocity and $p$ the pressure. System (1) depends on two strictly non-negative parameters, $\kappa_0$ and $\rho_0$. The wave velocity is $c_0$, linked with the parameters of the system by $c_0^2 = \kappa_0/\rho_0$. By taking the curl of the velocity equation of (1), we formally find $\partial_{\tau} (\nabla_{\mathbf{x}} \times \mathbf{u}) = 0$, which means that the vorticity $\nabla_{\mathbf{x}} \times \mathbf{u}$ is conserved. Conservation of this quantity is a necessary condition for having long time accuracy.

The problem of conservation of the vorticity for (1) has been tackled by several strategies, some of them being linked with the conservation of the divergence of a vector $\mathbf{v}$ when $\mathbf{v}$ ensures the following type of equation

$$ \partial_{\tau} \mathbf{v} + \nabla_{\mathbf{x}} \times \mathbf{g} = 0. \qquad (2) $$

These different strategies may be gathered into the three following families

1. **Projection method.** Projection methods have been widely used in the context of preservation of a divergence for incompressible flows [6] but also for magnetohydrodynamic system [10]. This method was rather designed for (2), and consists in computing a candidate update $\bar{\mathbf{v}}^{n+1}$ which is then corrected. For ensuring preservation of the divergence, a potential $\varphi$ is computed such that

$$ \Delta \varphi = \mathrm{div}_{\mathbf{x}} (\bar{\mathbf{v}}^{n+1} - \mathbf{v}^n), \qquad (3) $$

and the update is projected as $\mathbf{v}^{n+1} = \bar{\mathbf{v}}^{n+1} - \nabla_{\mathbf{x}}\varphi$, which leads to the conservation of the divergence of $\mathbf{v}$. As far as we know, it has never been proposed for the conservation of the curl for (1), but could be easily adapted. This method raises several problems: first, this requires to solve the elliptic problem (3) at each time step, which is costly, and second, (3) should be equipped with appropriate boundary conditions, which is not always straightforward.

2. **Discretizations based on staggered data.** This type of discretization consists in having data located on the faces or edges of the cells, whereas other data are located in the center of the cells. For example, the MAC scheme [28] consists in discretizing the pressure at the center of cells whereas the velocity components are normal to each face. The MAC scheme has been thoroughly analyzed in [33, 34], and extended to compressible flows in [18, 19, 17]. Similar approach were designed in [22, 30], in which general definitions of discrete gradient, divergence and curl were provided with a staggered design. Application of staggered schemes to hyperbolic problems can be found for example in [8]. A high order discontinuous Galerkin version of staggered schemes was proposed in [37], and high order finite volume schemes were proposed in [2]. Still linked with staggered ideas, the scheme [7] was proposed, for which the data are collocated, but the divergence is cleaned on a staggered grid. The main drawback of staggered data is its difficulty of implementation on unstructured meshes at high order.

3. **Generalized Lagrange multipliers.** The general Lagrange multiplier method was designed originally for the conservation of the divergence for the Maxwell system [32] and for the MHD system [14]. This method consists in adding an extra variable replacing the divergence to be preserved, and to define an extra equation on this additional variable, which includes a relaxation between the divergence of the vector and the additional variable. The extension to the curl constraint preservation was proposed in [16]. The main drawbacks of this method are the increase of the number of unknowns to solve, and the addition of numerical parameters, that should be tuned (propagation speed of the additional variable, stiffness of the relaxation).
Apart from these types of schemes that are especially designed for conserving divergence or curl constraints, some cell-centered finite volume schemes seem to be able to naturally conserve these constraint. Even if it was not expressed in these terms, [16, 21], which was then extended to high order in [29] are examples of finite volumes and discontinuous Galerkin methods that are long time accurate for (1), and so must preserve the curl in a sense that was not yet defined. Note that these articles deal with triangular or tetrahedral meshes, and are known to fail if other type of meshes are used (see e.g. [29], in which the high order discontinuous Galerkin method on quadrangular and triangular meshes are compared). A second family of collocated schemes seem to be able to preserve correctly divergence or curl constraints: the node-based solvers, such as [9] for divergence preservation or [3, 4, 5] for vorticity preservation and application to the low Mach number flows. These schemes may be seen as an extension of [25, 44, 43] to non Cartesian meshes.

In this article, we wish to find an extension of the schemes studied in [16, 21, 29] to quadrangular meshes. Note that as we are dealing with the two-dimensional case, the definition of the curl may be tricky. We therefore define two operators:

* the rotational operator, denoted as $\mathrm{rot}_x$ which maps vector fields to scalar fields as

$$ \mathrm{rot}_x \mathbf{u} := \frac{\partial \mathbf{u}_y}{\partial x} - \frac{\partial \mathbf{u}_x}{\partial y}, $$

* the curl operator, denoted as $\mathrm{curl}_x$, which maps scalar fields to vector fields as

$$ \mathrm{curl}_x f = \left( -\frac{\partial f}{\partial y}, \frac{\partial f}{\partial x} \right)^T . $$

These two operators are such that $\mathrm{rot}_x$ is the opposite of the adjoint of $\mathrm{curl}_x$ and inversely. Also, the quantity that is preserved for (1) is $\mathrm{rot}_x \mathbf{u}$. The theoretical study of first order schemes [16, 21] can rely on a discrete Hodge-Helmholtz decomposition that reads on periodic domains as

$$ d\mathbb{P}_0 = \mathbb{R}^2 \oplus \mathrm{curl}_x \mathbb{P}_1 \oplus \nabla_x \mathrm{CR}, \qquad (4) $$

where $d\mathbb{P}_0$ is the set of piecewise constant vectors, $\mathbb{P}_1$ is the space of continuous finite element of degree 1, and $\mathrm{CR}$ is the Crouzeix-Raviart finite element space [13]. The focus of this article is on finding a similar discrete Hodge-Helmholtz decomposition as (4), but on quadrangular meshes. For this, we propose to enrich the approximation space of velocity.

This article is organized as follows. In section 2, the new approximation space is introduced for Cartesian meshes, and the Cartesian meshes counterpart of (4) is proven. A finite volume numerical scheme for (1) is proposed with this new approximation space, and we prove that this numerical scheme preserves a curl. An originality of this curl is that it is defined in an adjoint sense. Then the scheme is extended in section 3 to the quadrangular case with non periodic boundary conditions in the spirit of [27]. In section 4, the numerical scheme is extended to the barotropic Euler system. Because the Euler system is nonlinear, and because the new basis for vectors is not piecewise constant, the numerical scheme is no more a purely finite volume scheme, but includes a cell integral as in the discontinuous Galerkin method. The section 5 is dedicated to numerical results where tests are performed on the wave system with periodic Cartesian meshes and with general quadrangular meshes. Tests are also performed for the accuracy problem at low Mach number for stationary problems and for propagation of acoustic waves in a low Mach number flow. This article finishes with the section 6, a conclusion.

## 2 Cartesian and periodic case

On triangular mesh with periodic boundary conditions, the discrete Hodge-Helmholtz decomposition (4) holds. Moreover, it is adapted to the wave system (1) in the sense that its divergence free component is preserved over time using a Godunov numerical scheme [16]. On Cartesian mesh, this decomposition does not exist anymore for piecewise constant functions. In this section, a new finite element space that is richer than piecewise constant functions is defined and allows to recover a discrete Hodge-Helmholtz decomposition. Using this vectorial approximation space, a numerical discretization of the wave system (1) is proposed. Then, it is proved that this discretization coupled with a Godunov numerical flux preserves the divergence free component of the discrete Hodge-Helmholtz decomposition but also the adjoint curl over time.
Apart from these types of schemes that are especially designed for conserving divergence or curl constraints, some cell-centered finite volume schemes seem to be able to naturally conserve these constraint. Even if it was not expressed in these terms, [15, 20], which was then extended to high order in [27] are examples of finite volumes and discontinuous Galerkin methods that are long time accurate for (1), and so must preserve the curl in a sense that was not yet defined. Note that these articles deal with triangular or tetrahedral meshes, and are known to fail if other type of meshes are used (see e.g. [27], in which the high order discontinuous Galerkin method on quadrangular and triangular meshes are compared). A second family of collocated schemes seem to be able to preserve correctly divergence or curl constraints: the node-based solvers, such as [9] for divergence preservation or [3, 4, 5] for vorticity preservation and application to the low Mach number flows. These schemes may be seen as an extension of [23, 40, 39] to non Cartesian meshes.

In this article, we wish to find an extension of the schemes studied in [15, 20, 27] to quadrangular meshes. The theoretical study of first order schemes [15, 20] can rely on a discrete Hodge-Helmholtz decomposition that reads on periodic domains as

$$d\mathbb{P}_0 = \mathbb{R}^2 \oplus \nabla_{\mathbf{x}}^{\perp} \mathbb{P}_1 \oplus \nabla_{\mathbf{x}} \mathbb{C}\mathbb{R}, \qquad (4)$$

where $d\mathbb{P}_0$ is the set of piecewise constant vectors, $\mathbb{P}_1$ is the space of continuous finite element of degree 1, and $\mathbb{C}\mathbb{R}$ is the Crouzeix-Raviart finite element space [13]. The focus of this article is on finding a similar discrete Hodge-Helmholtz decomposition as (4), but on quadrangular meshes. For this, we propose to enrich the approximation space of velocity.

This article is organized as follows. In section 2, the new approximation space is introduced for Cartesian meshes, and the Cartesian meshes counterpart of (4) is proven. A finite volume numerical scheme for (1) is proposed with this new approximation space, and we prove that this numerical scheme preserves a curl. An originality of this curl is that it is defined in an adjoint sense. Then the scheme is extended in section 3 to the quadrangular case with non periodic boundary conditions in the spirit of [25]. In section 4, the numerical scheme is extended to the barotropic Euler system. Because the Euler system is nonlinear, and because the new basis for vectors is not piecewise constant, the numerical scheme is no more a purely finite volume scheme, but includes a cell integral as in the discontinuous Galerkin method. The section 5 is dedicated to numerical results where tests are performed on the wave system with periodic Cartesian meshes and with general quadrangular meshes. Tests are also performed for the accuracy problem at low Mach number for stationary problems and for propagation of acoustic waves in a low Mach number flow. This article finishes with the section 6, a conclusion.

## 2 Cartesian and periodic case

On triangular mesh with periodic boundary conditions, the discrete Hodge-Helmholtz decomposition (4) holds. Moreover, it is adapted to the wave system (1) in the sense that its divergence free component is preserved over time using a Godunov numerical scheme [15]. On Cartesian mesh, this decomposition does not exist anymore for piecewise constant functions. In this section, a new finite element space that is richer than piecewise constant functions is defined and allows to recover a discrete Hodge-Helmholtz decomposition. Using this vectorial approximation space, a numerical discretization of the wave system (1) is proposed. Then, it is proved that this discretization coupled with a Godunov numerical flux preserves the divergence free component of the discrete Hodge-Helmholtz decomposition but also the adjoint curl over time.

### 2.1 A new finite element space for vectorial field

Let $\hat{K} = [0; 1]^2$ the unit square element. The finite dimensional space vector fields $\hat{V}_h$ on $\hat{K}$ considered matches with the space $\hat{S}_0(\hat{K})$ defined in [1] by

$$\hat{S}_0(\hat{K}) = \mathbb{Q}_0(\hat{K})^2 + \operatorname{span}\left(\begin{pmatrix} \hat{x} - 1/2 \\ -(\hat{y} - 1/2) \end{pmatrix}\right) = \operatorname{span}\left(\begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} \hat{x} - 1/2 \\ -(\hat{y} + 1/2) \end{pmatrix}\right). \quad (5)$$
## 2.1 A new finite element space for vectorial field

Let $\hat{K} = [0; 1]^2$ the unit square element. The finite dimensional space vector fields $\hat{V}_h$ on $\hat{K}$ considered matches with the space $\hat{S}_0(\hat{K})$ defined in [1] by

$$ \hat{S}_0(\hat{K}) = \mathbb{Q}_0(\hat{K})^2 + \operatorname{span} \left( \begin{pmatrix} \hat{x} - 1/2 \\ -(\hat{y} - 1/2) \end{pmatrix} \right) = \operatorname{span} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} \hat{x} - 1/2 \\ -\hat{y} + 1/2 \end{pmatrix} \right). \quad (5) $$

We denote by $\Sigma = \{\sigma_1, \sigma_2, \sigma_3\}$ the set of linear forms defined on $L^2(\hat{K})^2$ by

$$ \sigma_1(\hat{u}) = \int_{\hat{K}} \hat{u}(\hat{x}) \cdot \begin{pmatrix} 1 \\ 0 \end{pmatrix} d\hat{x}, \quad (6) $$

$$ \sigma_2(\hat{u}) = \int_{\hat{K}} \hat{u}(\hat{x}) \cdot \begin{pmatrix} 0 \\ 1 \end{pmatrix} d\hat{x}, \quad (7) $$

$$ \sigma_3(\hat{u}) = 12 \int_{\hat{K}} \hat{u}(\hat{x}) \cdot \begin{pmatrix} \hat{x} - 1/2 \\ -\hat{y} + 1/2 \end{pmatrix} d\hat{x}. \quad (8) $$

We directly obtain that $(\hat{K}, \Sigma, \hat{S}_0(\hat{K}))$ is a finite element. The finite element $(\hat{K}, \Sigma, \hat{S}_0(\hat{K}))$ is then strictly included in the Raviart-Thomas $RT_0$ finite element of degree zero [40, 1].

## 2.2 Approximation space on Cartesian mesh

In this section, we are interested in the Cartesian mesh case. For the sake of simplicity, the domain considered is a square of size $l \times l$, divided into $N$ cells in each direction. We denote by $(x_0, y_0)$ the bottom left point of the domain. Then the domain is $\Omega = [x_0; x_0 + l] \times [y_0; y_0 + l]$, the nodes of the mesh are given by the couples $(x_i, y_j)$ where

$$ x_i = x_0 + ih, \quad i = 0, 1, \dots, N $$

$$ y_j = y_0 + jh, \quad i = 0, 1, \dots, N $$

with $h = \frac{l}{N}$. All cell $c_{i,j} = [x_i; x_{i+1}] \times [y_j; y_{j+1}]$ is the image of $\hat{K} = [0; 1]^2$ under the following linear application

$$ \mathbf{F}_{c_{i,j}} : \begin{pmatrix} \hat{x} \\ \hat{y} \end{pmatrix} \mapsto \begin{pmatrix} x_i + \hat{x}h \\ y_j + \hat{y}h \end{pmatrix}. $$

We denote by $\mathcal{C}$ the set of cells of the mesh and by $\mathcal{S}$ the set of sides. For a given vectorial approximation space $\hat{V}$ on the reference square, the approximation space of vector fields on $\Omega$ is defined as

$$ \mathbf{V}_h = \{ \mathbf{u}_h \in L^2(\Omega)^2 \mid \forall c \in \mathcal{C}, \mathbf{u}_h|_c \circ \mathbf{F}_c \in \hat{\mathbf{V}} \}. $$

The average and jump of $\mathbf{u}_h \in \mathbf{V}_h$ on an interior side $S \in \mathcal{S}$ separating two cells $c_L$ and $c_R$ are defined by

$$ \llbracket \mathbf{u}_h \rrbracket = \mathbf{u}_h|_{c_L} - \mathbf{u}_h|_{c_R}, \quad \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = (\mathbf{u}_h|_{c_L} - \mathbf{u}_h|_{c_R}) \cdot \mathbf{n}^S, \\ \{ \mathbf{u}_h \} = \frac{\mathbf{u}_h|_{c_L} + \mathbf{u}_h|_{c_R}}{2}, \quad \{ \mathbf{u}_h \cdot \mathbf{n}^S \} = \frac{\mathbf{u}_h|_{c_L} + \mathbf{u}_h|_{c_R}}{2} \cdot \mathbf{n}^S. $$

where $\mathbf{n}^S$ is the normal vector to $S$ going from cell $c_L$ to cell $c_R$. We denote by $\mathcal{C}(S)$ the set of the two adjacent cells to a side $S \in \mathcal{S}$ and by $\mathcal{S}(c)$ the set of the four sides of a cell $c \in \mathcal{C}$.

It was proven in [1] that the best approximation of $\mathbf{u}$ in $\mathbf{V}_h$ is of order 1 for the $L^2(\Omega)$ norm if and only if $\hat{V}$ contains $\hat{S}_0$. This leads to consider the following approximation space

$$ d\mathbf{S}_h(\Omega) := \{ \mathbf{u}_h \in L^2(\Omega)^2 \mid \forall c \in \mathcal{C}, \mathbf{u}_h|_c \circ \mathbf{F}_c \in \hat{\mathbf{S}}_0 \}. \quad (9) $$

**Proposition 1** (Properties $d\mathbf{S}_h(\Omega)$). $d\mathbf{S}_h(\Omega)$ is a finite element space approximating $L^2(\Omega)^2$ at order one, and its dimension is

$$ \dim(d\mathbf{S}_h(\Omega)) = 3\#\mathcal{C} = 3N^2. $$

Also, the following properties hold for all $\mathbf{u}_h \in d\mathbf{S}_h$

$$ \bullet \quad \forall c \in \mathcal{C}, \quad \operatorname{div}_x (\mathbf{u}_h|_c) = 0 $$

$$ \bullet \quad \forall c \in \mathcal{C}, \quad \forall S \in \mathcal{S}(c), \quad \mathbf{u}_h|_c \cdot \mathbf{n}^S \text{ is constant along } S. \quad \text{This also means that the jump } \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket \text{ is constant along each side.} $$

*Proof.* All these properties are direct consequences of straightforward computations. $\square$
We denote by $\Sigma = \{\sigma_1, \sigma_2, \sigma_3\}$ the set of linear forms defined on $L^2(\hat{K})^2$ by

$$ \sigma_1(\hat{\mathbf{u}}) = \int_{\hat{K}} \hat{\mathbf{u}}(\hat{\mathbf{x}}) \cdot \begin{pmatrix} 1 \\ 0 \end{pmatrix} d\hat{\mathbf{x}}, \qquad (6) $$

$$ \sigma_2(\hat{\mathbf{u}}) = \int_{\hat{K}} \hat{\mathbf{u}}(\hat{\mathbf{x}}) \cdot \begin{pmatrix} 0 \\ 1 \end{pmatrix} d\hat{\mathbf{x}}, \qquad (7) $$

$$ \sigma_3(\hat{\mathbf{u}}) = 12 \int_{\hat{K}} \hat{\mathbf{u}}(\hat{\mathbf{x}}) \cdot \begin{pmatrix} \hat{x} - 1/2 \\ -\hat{y} + 1/2 \end{pmatrix} d\hat{\mathbf{x}}. \qquad (8) $$

We directly obtain that $(\hat{K}, \Sigma, \hat{\mathbf{S}}_0(\hat{K}))$ is a finite element. The finite element $(\hat{K}, \Sigma, \hat{\mathbf{S}}_0(\hat{K}))$ is then strictly included in the Raviart-Thomas $RT_0$ finite element of degree zero [36, 1].

## 2.2 Approximation space on Cartesian mesh

In this section, we are interested in the Cartesian mesh case. For the sake of simplicity, the domain considered is a square of size $l \times l$, divided into $N$ cells in each direction. We denote by $(x_0, y_0)$ the bottom left point of the domain. Then the domain is $\Omega = [x_0; x_0 + l] \times [y_0; y_0 + l]$, the nodes of the mesh are given by the couples $(x_i, y_j)$ where

$$ x_i = x_0 + ih, \quad i = 0, 1, \dots, N $$

$$ y_j = y_0 + jh, \quad i = 0, 1, \dots, N $$

with $h = \frac{l}{N}$. All cell $c_{i,j} = [x_i; x_{i+1}] \times [y_j; y_{j+1}]$ is the image of $\hat{K} = [0; 1]^2$ under the following linear application

$$ \mathbf{F}_{c_{i,j}} : \begin{pmatrix} \hat{x} \\ \hat{y} \end{pmatrix} \mapsto \begin{pmatrix} x_i + \hat{x}h \\ y_j + \hat{y}h \end{pmatrix}. $$

We denote by $\mathcal{C}$ the set of cells of the mesh and by $\mathcal{S}$ the set of sides. For a given vectorial approximation space $\hat{\mathbf{V}}$ on the reference square, the approximation space of vector fields on $\Omega$ is defined as

$$ \mathbf{V}_h = \{ \mathbf{u}_h \in L^2(\Omega)^2 \mid \forall c \in \mathcal{C}, \mathbf{u}_h|_c \circ \mathbf{F}_c \in \hat{\mathbf{V}} \}. $$

The average and jump of $\mathbf{u}_h \in \mathbf{V}_h$ on an interior side $S \in \mathcal{S}$ separating two cells $c_L$ and $c_R$ are defined by

$$ \llbracket \mathbf{u}_h \rrbracket = \mathbf{u}_h|_{c_L} - \mathbf{u}_h|_{c_R}, \quad \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = (\mathbf{u}_h|_{c_L} - \mathbf{u}_h|_{c_R}) \cdot \mathbf{n}^S, \\ \llbracket \mathbf{u}_h \rrbracket = \frac{\mathbf{u}_h|_{c_L} + \mathbf{u}_h|_{c_R}}{2}, \quad \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = \frac{\mathbf{u}_h|_{c_L} + \mathbf{u}_h|_{c_R}}{2} \cdot \mathbf{n}^S. $$

where $\mathbf{n}^S$ is the normal vector to side $S$ going from cell $c_L$ to cell $c_R$. We denote by $\mathcal{C}(S)$ the set of the two adjacent cells to a side $S \in \mathcal{S}$ and by $\mathcal{S}(c)$ the set of the four sides of a cell $c \in \mathcal{C}$.

It was proven in [1] that the best approximation of $\mathbf{u}$ in $\mathbf{V}_h$ is of order 1 for the $L^2(\Omega)$ norm if and only if $\hat{\mathbf{V}}$ contains $\hat{\mathbf{S}}_0$. This leads to consider the following approximation space

$$ d\mathbf{S}_h(\Omega) := \{ \mathbf{u}_h \in L^2(\Omega)^2 \mid \forall c \in \mathcal{C}, \mathbf{u}_h|_c \circ \mathbf{F}_c \in \hat{\mathbf{S}}_0 \}. \qquad (9) $$

**Proposition 1** (Properties $d\mathbf{S}_h(\Omega)$). $d\mathbf{S}_h(\Omega)$ is a finite element space approximating $L^2(\Omega)^2$ at order one, and its dimension is

$$ \dim(d\mathbf{S}_h(\Omega)) = 3\#\mathcal{C} = 3N^2. $$

Also, the following properties hold for all $\mathbf{u}_h \in d\mathbf{S}_h$

$$ \bullet \quad \forall c \in \mathcal{C}, \quad \mathrm{div}_x (\mathbf{u}_h|_c) = 0 $$

$$ \bullet \quad \forall c \in \mathcal{C}, \forall S \in \mathcal{S}(c), \mathbf{u}_h|_c \cdot \mathbf{n}^S \text{ is constant along } S. \text{ This also means that the jump } \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket \text{ is constant along each side.} $$

*Proof.* All these properties are direct consequences of straightforward computations. □
<figure><img src="image_1.png" /><figcaption>Figure 1: Labels of the vertices (1, 2, 3, 4) and the sides of the references square (B,R,T,L for Bottom, Right, Top, Left).</figcaption></figure>

## 2.3 Discrete Hodge-Helmholtz decomposition

In this section, a discrete Hodge-Helmholtz decomposition of the new approximation space $dS_h(\Omega)$ is proposed. This result is the Cartesian version of the decomposition (4) that holds on triangular mesh. We denote by $Q_1(\Omega)$ the set of continuous, piecewise bilinear Lagrange finite element

$$Q_1(\Omega) = \{ \psi \in C^0(\Omega) \mid \forall c \in \mathcal{C}, \hat{\psi}_c := \psi|_c \circ \mathbf{F}_c \in Q_1(\hat{K}) \},$$

where $Q_1(\hat{K})$ is spanned by the following basis

$$\left\{ \begin{array}{l} \hat{\psi}_1(\hat{x}, \hat{y}) = (\hat{x} - 1)(\hat{y} - 1), \\ \hat{\psi}_2(\hat{x}, \hat{y}) = -\hat{x}(\hat{y} - 1), \\ \hat{\psi}_3(\hat{x}, \hat{y}) = \hat{x}\hat{y}, \\ \hat{\psi}_4(\hat{x}, \hat{y}) = -(\hat{x} - 1)\hat{y}, \end{array} \right. \qquad (10)$$

see Figure 1 for the labelling of the vertices. On a Cartesian mesh with periodic boundary conditions, the dimension of $Q_1(\Omega)$ is the number of vertices, namely

$$\dim Q_1(\Omega) = N^2. \qquad (11)$$

We denote by $\mathbb{R}aTu(\Omega)$ the Rannacher-Turek [39] finite element space

$$\mathbb{R}aTu(\Omega) = \{ \varphi \in L^2(\Omega) \mid \forall c \in \mathcal{C}, \hat{\varphi}_c := \varphi|_c \circ \mathbf{F}_c \in \mathbb{R}aTu(\hat{K}) \text{ and } \forall S \in \mathcal{S}, \int_S [\varphi] = 0 \}.$$

The space $\mathbb{R}aTu(\hat{K})$ is spanned by the following basis functions

$$\left\{ \begin{array}{l} \hat{\varphi}_B(\hat{x}, \hat{y}) = \frac{3}{4} - \hat{y} - \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \\ \hat{\varphi}_R(\hat{x}, \hat{y}) = -\frac{1}{4} + \hat{x} + \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \\ \hat{\varphi}_T(\hat{x}, \hat{y}) = -\frac{1}{4} + \hat{y} - \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \\ \hat{\varphi}_L(\hat{x}, \hat{y}) = \frac{3}{4} - \hat{x} + \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \end{array} \right. \qquad (12)$$

and its degrees of freedom are the integrals on each side. On a Cartesian mesh with periodic boundary conditions, the dimension of $\mathbb{R}aTu(\Omega)$ is the number of sides, namely

$$\dim \mathbb{R}aTu(\Omega) = 2N^2. \qquad (13)$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Labels of the vertices (1, 2, 3, 4) and the sides of the references square (B,R,T,L for Bottom, Right, Top, Left).</figcaption></figure>

## 2.3 Discrete Hodge-Helmholtz decomposition

In this section, a discrete Hodge-Helmholtz decomposition of the new approximation space $dS_h(\Omega)$ is proposed. This result is the Cartesian version of the decomposition (4) that holds on triangular mesh. We denote by $Q_1(\Omega)$ the set of continuous, piecewise bilinear Lagrange finite element

$$Q_1(\Omega) = \{ \psi \in C^0(\Omega) \mid \forall c \in \mathcal{C}, \hat{\psi}_c := \psi|_c \circ \mathbf{F}_c \in Q_1(\hat{K}) \},$$

where $Q_1(\hat{K})$ is spanned by the following basis

$$\left\{ \begin{array}{l} \hat{\psi}_1(\hat{x}, \hat{y}) = (\hat{x} - 1)(\hat{y} - 1), \\ \hat{\psi}_2(\hat{x}, \hat{y}) = -\hat{x}(\hat{y} - 1), \\ \hat{\psi}_3(\hat{x}, \hat{y}) = \hat{x}\hat{y}, \\ \hat{\psi}_4(\hat{x}, \hat{y}) = -(\hat{x} - 1)\hat{y}, \end{array} \right. \qquad (10)$$

see Figure 1 for the labelling of the vertices. On a Cartesian mesh with periodic boundary conditions, the dimension of $Q_1(\Omega)$ is the number of vertices, namely

$$\dim Q_1(\Omega) = N^2. \qquad (11)$$

We denote by $\mathbb{R}aTu(\Omega)$ the Rannacher-Turek [35] finite element space

$$\mathbb{R}aTu(\Omega) = \{ \varphi \in L^2(\Omega) \mid \forall c \in \mathcal{C}, \hat{\varphi}_c := \varphi|_c \circ \mathbf{F}_c \in \mathbb{R}aTu(\hat{K}) \text{ and } \forall S \in \mathcal{S}, \int_S [\varphi] = 0 \}.$$

The space $\mathbb{R}aTu(\hat{K})$ is spanned by the following basis functions

$$\left\{ \begin{array}{l} \hat{\varphi}_B(\hat{x}, \hat{y}) = \frac{3}{4} - \hat{y} - \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \\ \hat{\varphi}_R(\hat{x}, \hat{y}) = -\frac{1}{4} + \hat{x} + \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \\ \hat{\varphi}_T(\hat{x}, \hat{y}) = -\frac{1}{4} + \hat{y} - \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \\ \hat{\varphi}_L(\hat{x}, \hat{y}) = \frac{3}{4} - \hat{x} + \frac{3}{2} \left( \left( \hat{x} - \frac{1}{2} \right)^2 - \left( \hat{y} - \frac{1}{2} \right)^2 \right), \end{array} \right. \qquad (12)$$

and its degrees of freedom are the integrals on each side. On a Cartesian mesh with periodic boundary conditions, the dimension of $\mathbb{R}aTu(\Omega)$ is the number of sides, namely

$$\dim \mathbb{R}aTu(\Omega) = 2N^2. \qquad (13)$$
**Proposition 2** (Exact discrete Hodge-Helmholtz decomposition of $d\mathbf{S}_h(\Omega)$ in the Cartesian periodic case). On a uniform Cartesian mesh $\mathcal{M}_h$ of $\Omega$, the following decomposition holds

$$d\mathbf{S}_h(\Omega) = \mathbb{R}^2 \oplus \nabla_{\mathbf{x}} (\mathrm{RaTu}(\Omega)) \oplus \mathrm{curl}_{\mathbf{x}} (\mathbb{Q}_1(\Omega)).$$

More precisely, a velocity field $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$ can be uniquely written as

$$\mathbf{u}_h = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} + \mathcal{P}_h^{\varphi}[\mathbf{u}_h] + \mathcal{P}_h^{\psi}[\mathbf{u}_h]$$

where $(\alpha, \beta) \in \mathbb{R}^2$ and

* $\mathcal{P}_h^{\varphi}[\mathbf{u}_h]$ is the gradient of a Rannacher-Turek scalar potential, i.e. $\mathcal{P}_h^{\varphi}[\mathbf{u}_h] = \nabla_{\mathbf{x}}\varphi_h$ with $\varphi_h \in \mathrm{RaTu}(\Omega)$ where $\nabla_{\mathbf{x}}$ is applied cellwise:

$$\forall c \in \mathcal{C} \quad (\nabla_{\mathbf{x}}\varphi_h)|_c := \nabla_{\mathbf{x}}(\varphi_h)|_c$$

* $\mathcal{P}_h^{\psi}[\mathbf{u}_h]$ is the curl of a continuous $\mathbb{Q}_1$ scalar function, i.e. $\mathcal{P}_h^{\psi}[\mathbf{u}_h] = \mathrm{curl}_{\mathbf{x}}(\psi_h)$ with $\psi_h \in \mathbb{Q}_1(\Omega)$ where $\mathrm{curl}_{\mathbf{x}}(\psi) = (-\partial_y\psi, \partial_x\psi)^T$ is applied cellwise:

$$\forall c \in \mathcal{C} \quad (\mathrm{curl}_{\mathbf{x}} \psi_h)|_c = \mathrm{curl}_{\mathbf{x}} (\psi_h)|_c.$$

Moreover, this decomposition is orthogonal for the $L^2(\Omega)$ scalar product

$$\langle \mathbf{u}, \mathbf{v} \rangle = \int_{\Omega} \mathbf{u} \cdot \mathbf{v}.$$

*Proof.* We first prove that if $\psi \in \mathbb{Q}_1(\Omega)$, $\mathrm{curl}_{\mathbf{x}} \psi \in d\mathbf{S}_h$. Since for all cell $c$, we have $\psi|_c = \hat{\psi}_c \circ \mathbf{F}_c^{-1}$ where $\hat{\psi}_c \in \mathbb{Q}_1(\hat{K})$, we have

$$\mathrm{curl}_{\mathbf{x}} (\psi|_c) = \frac{1}{h} \mathrm{curl}_{\hat{\mathbf{x}}} \hat{\psi}_c. \qquad (14)$$

Then it is sufficient to prove that the curl of each basis function of $\mathbb{Q}_1(\hat{K})$ is included in $\hat{\mathbf{S}}_0$. Taking the curl of each of the functions defined in (10), we get

$$\left\{ \begin{array}{l} \mathrm{curl}_{\hat{\mathbf{x}}} (\hat{\psi}_1) = \begin{pmatrix} -\hat{x} + 1 \\ \hat{y} - 1 \end{pmatrix} \in \hat{\mathbf{S}}_0, \\ \mathrm{curl}_{\hat{\mathbf{x}}} (\hat{\psi}_2) = \begin{pmatrix} \hat{x} \\ -\hat{y} + 1 \end{pmatrix} \in \hat{\mathbf{S}}_0, \\ \mathrm{curl}_{\hat{\mathbf{x}}} (\hat{\psi}_3) = \begin{pmatrix} -\hat{x} \\ \hat{y} \end{pmatrix} \in \hat{\mathbf{S}}_0, \\ \mathrm{curl}_{\hat{\mathbf{x}}} (\hat{\psi}_4) = \begin{pmatrix} \hat{x} - 1 \\ -\hat{y} \end{pmatrix} \in \hat{\mathbf{S}}_0. \end{array} \right.$$

We now prove that if $\varphi \in \mathrm{RaTu}(\Omega)$, we have $\nabla_{\mathbf{x}}\varphi \in d\mathbf{S}_h$. Since for all cell $c$, $\varphi|_c = \hat{\varphi}_c \circ \mathbf{F}_c^{-1}$ where $\hat{\varphi}_c \in \mathrm{RaTu}(\hat{K})$, we have

$$\nabla_{\mathbf{x}} (\varphi|_c) = \frac{1}{h} \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_c. \qquad (15)$$

Then, it is also sufficient to prove that the gradient of each basis function of the Rannacher-Turek finite element on $\hat{K}$ is included in $\hat{\mathbf{S}}_0(\hat{K})$. Taking the derivative of the basis functions $(\hat{\varphi}_B, \hat{\varphi}_R, \hat{\varphi}_T, \hat{\varphi}_L)$ defined in (12) leads to

$$\left\{ \begin{array}{l} \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_B = \begin{pmatrix} -3\hat{x} + 3/2 \\ 3\hat{y} - 5/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}), \\ \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_R = \begin{pmatrix} 3\hat{x} - 1/2 \\ -3\hat{y} + 3/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}), \\ \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_T = \begin{pmatrix} -3\hat{x} + 3/2 \\ 3\hat{y} - 1/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}), \\ \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_L = \begin{pmatrix} 3\hat{x} - 5/2 \\ -3\hat{y} + 3/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}). \end{array} \right.$$
**Proposition 2** (Exact discrete Hodge-Helmholtz decomposition of $d\mathbf{S}_h(\Omega)$ in the Cartesian periodic case). On a uniform Cartesian mesh $\mathcal{M}_h$ of $\Omega$, the following decomposition holds

$$d\mathbf{S}_h(\Omega) = \mathbb{R}^2 \oplus \nabla_{\mathbf{x}} (\mathbb{R}\mathrm{aTu}(\Omega)) \oplus \nabla_{\mathbf{x}}^{\perp} (\mathbb{Q}_1(\Omega)).$$

More precisely, a velocity field $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$ can be uniquely written as

$$\mathbf{u}_h = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} + \mathcal{P}_h^{\varphi}[\mathbf{u}_h] + \mathcal{P}_h^{\psi}[\mathbf{u}_h]$$

where $(\alpha, \beta) \in \mathbb{R}^2$ and

* $\mathcal{P}_h^{\varphi}[\mathbf{u}_h]$ is the gradient of a Rannacher-Turek scalar potential, i.e. $\mathcal{P}_h^{\varphi}[\mathbf{u}_h] = \nabla_{\mathbf{x}}\varphi_h$ with $\varphi_h \in \mathbb{R}\mathrm{aTu}(\Omega)$ where $\nabla_{\mathbf{x}}$ is applied cellwise:

$$\forall c \in \mathcal{C} \quad (\nabla_{\mathbf{x}}\varphi_h)|_c := \nabla_{\mathbf{x}}(\varphi_h)|_c$$

* $\mathcal{P}_h^{\psi}[\mathbf{u}_h]$ is the curl of a continuous $\mathbb{Q}_1$ scalar function, i.e. $\mathcal{P}_h^{\psi}[\mathbf{u}_h] = \nabla_{\mathbf{x}}^{\perp}(\psi_h)$ with $\psi_h \in \mathbb{Q}_1(\Omega)$ where $\nabla_{\mathbf{x}}^{\perp}(\psi) = (-\partial_y\psi, \partial_x\psi)^T$ is applied cellwise:

$$\forall c \in \mathcal{C} \quad (\nabla_{\mathbf{x}}^{\perp}\psi_h)|_c = \nabla_{\mathbf{x}}^{\perp}(\psi_h)|_c.$$

Moreover, this decomposition is orthogonal in $L^2(\Omega)^2$.

*Proof.* We first prove that if $\psi \in \mathbb{Q}_1(\Omega)$, $\nabla_{\mathbf{x}}^{\perp}\psi \in d\mathbf{S}_h$. Since for all cell $c$, we have $\psi|_c = \hat{\psi}_c \circ \mathbf{F}_c^{-1}$ where $\hat{\psi}_c \in \mathbb{Q}_1(\hat{K})$, we have

$$\nabla_{\mathbf{x}}^{\perp}(\psi|_c) = \frac{1}{h} \nabla_{\hat{\mathbf{x}}}^{\perp} \hat{\psi}_c. \qquad (14)$$

Then it is sufficient to prove that the curl of each basis function of $\mathbb{Q}_1(\hat{K})$ is included in $\hat{\mathbf{S}}_0$. Taking the curl of each of the functions defined in (10), we get

$$\left\{ \begin{aligned} \nabla_{\hat{\mathbf{x}}}^{\perp}(\hat{\psi}_1) &= \begin{pmatrix} -\hat{x} + 1 \\ \hat{y} - 1 \end{pmatrix} \in \hat{\mathbf{S}}_0, \\ \nabla_{\hat{\mathbf{x}}}^{\perp}(\hat{\psi}_2) &= \begin{pmatrix} \hat{x} \\ -\hat{y} + 1 \end{pmatrix} \in \hat{\mathbf{S}}_0, \\ \nabla_{\hat{\mathbf{x}}}^{\perp}(\hat{\psi}_3) &= \begin{pmatrix} -\hat{x} \\ \hat{y} \end{pmatrix} \in \hat{\mathbf{S}}_0, \\ \nabla_{\hat{\mathbf{x}}}^{\perp}(\hat{\psi}_4) &= \begin{pmatrix} \hat{x} - 1 \\ -\hat{y} \end{pmatrix} \in \hat{\mathbf{S}}_0. \end{aligned} \right.$$

We now prove that if $\varphi \in \mathbb{R}\mathrm{aTu}(\Omega)$, we have $\nabla_{\mathbf{x}}\varphi \in d\mathbf{S}_h$. Since for all cell $c$, $\varphi|_c = \hat{\varphi}_c \circ \mathbf{F}_c^{-1}$ where $\hat{\varphi}_c \in \mathbb{R}\mathrm{aTu}(\hat{K})$, we have

$$\nabla_{\mathbf{x}}(\varphi|_c) = \frac{1}{h} \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_c. \qquad (15)$$

Then, it is also sufficient to prove that the gradient of each basis function of the Rannacher-Turek finite element on $\hat{K}$ is included in $\hat{\mathbf{S}}_0(\hat{K})$. Taking the derivative of the basis functions $(\hat{\varphi}_B, \hat{\varphi}_R, \hat{\varphi}_T, \hat{\varphi}_L)$ defined in (12) leads to

$$\left\{ \begin{aligned} \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_B &= \begin{pmatrix} -3\hat{x} + 3/2 \\ 3\hat{y} - 5/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}), \\ \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_R &= \begin{pmatrix} 3\hat{x} - 1/2 \\ -3\hat{y} + 3/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}), \\ \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_T &= \begin{pmatrix} -3\hat{x} + 3/2 \\ 3\hat{y} - 1/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}), \\ \nabla_{\hat{\mathbf{x}}} \hat{\varphi}_L &= \begin{pmatrix} 3\hat{x} - 5/2 \\ -3\hat{y} + 3/2 \end{pmatrix} \in \hat{\mathbf{S}}_0(\hat{K}). \end{aligned} \right.$$

We now prove that the sum is orthogonal. For $\psi_h \in \mathbb{Q}_1(\Omega)$ and $\varphi_h \in \mathbb{R}\mathrm{aTu}(\Omega)$, we have

$$\begin{aligned} \int_{\Omega} (\nabla_{\mathbf{x}}\varphi_h) \cdot (\nabla_{\mathbf{x}}^{\perp}\psi_h) &= \sum_{c \in \mathcal{C}} \int_{c} (\nabla_{\mathbf{x}}\varphi_h)|_c \cdot (\nabla_{\mathbf{x}}^{\perp}\psi_h)|_c \\ &= \sum_{c \in \mathcal{C}} \sum_{S \in \mathcal{S}(c)} \int_{S} (\varphi_h)|_c \nabla_{\mathbf{x}}^{\perp}((\psi_h)|_c) \cdot \mathbf{n}^S \\ &\quad - \sum_{c \in \mathcal{C}} \int_{c} (\varphi_h)|_c \mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}^{\perp}((\psi_h)|_c)) \\ &= \sum_{S \in \mathcal{S}} \int_{S} [\varphi_h \nabla_{\mathbf{x}}^{\perp}(\psi_h) \cdot \mathbf{n}^S]. \end{aligned}$$
We now prove that the sum is orthogonal. For $\psi_h \in Q_1(\Omega)$ and $\varphi_h \in \mathbb{R}\mathit{aTu}(\Omega)$, we have

$$
\begin{align*}
\int_{\Omega} (\nabla_{\mathbf{x}} \varphi_h) \cdot (\mathbf{curl}_{\mathbf{x}} \psi_h) &= \sum_{c \in \mathcal{C}} \int_{c} (\nabla_{\mathbf{x}} \varphi_h)|_{c} \cdot (\mathbf{curl}_{\mathbf{x}} \psi_h)|_{c} \\
&= \sum_{c \in \mathcal{C}} \sum_{S \in \mathcal{S}(c)} \int_{S} (\varphi_h)|_{c} \mathbf{curl}_{\mathbf{x}} ((\psi_h)|_{c}) \cdot \mathbf{n}^{S} \\
&\quad - \sum_{c \in \mathcal{C}} \int_{c} (\varphi_h)|_{c} \mathrm{div}_{\mathbf{x}} (\mathbf{curl}_{\mathbf{x}} ((\psi_h)|_{c})) \\
&= \sum_{S \in \mathcal{S}} \int_{S} [\varphi_h \mathbf{curl}_{\mathbf{x}} (\psi_h) \cdot \mathbf{n}^{S}],
\end{align*}
$$

where the Gauss's theorem was used for the second computing line. As $\psi_h \in Q_1$, it is linear along each side, so that $\mathbf{curl}_x \psi_h \cdot \mathbf{n}^S$ is constant along each side. This leads to

$$
\int_{\Omega} (\nabla_{\mathbf{x}} \varphi_h) \cdot (\mathbf{curl}_{\mathbf{x}} \psi_h) = \sum_{S \in \mathcal{S}} \mathbf{curl}_{\mathbf{x}} (\psi_h) \cdot \mathbf{n}^S \int_{S} [\varphi_h].
$$

Last, by definition of the degrees of freedom of $\mathbb{R}\mathit{aTu}(\Omega)$, we have

$$
\forall S \in \mathcal{S} \quad \forall \varphi_h \in \mathbb{R}\mathit{aTu}(\Omega) \qquad \int_S [\varphi_h] = 0,
$$

which proves that $\nabla_{\mathbf{x}}\varphi_h$ and $\mathbf{curl}_{\mathbf{x}} \psi_h$ are orthogonal.

Following exactly the same lines, we can prove that any $(\alpha, \beta) \in \mathbb{R}^2$ is orthogonal to any element of $\nabla_{\mathbf{x}}\mathbb{R}\mathit{aTu}(\Omega)$.

It remains to deal with the orthogonality of $(\alpha, \beta) \in \mathbb{R}^2$ and any element of $\mathbf{curl}_x Q_1(\Omega)$. We consider $\psi_h$ an element of $Q_1(\Omega)$, and compute as follows

$$
\begin{align*}
\int_{\Omega} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \cdot \mathbf{curl}_{\mathbf{x}} \psi_h &= \sum_{c \in \mathcal{C}} \int_{c} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \cdot \mathbf{curl}_{\mathbf{x}} \psi_h \\
&= - \sum_{c \in \mathcal{C}} \int_{c} \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \nabla_{\mathbf{x}} \psi_h \\
&= - \sum_{c \in \mathcal{C}} \int_{c} \mathrm{div}_{\mathbf{x}} \left( \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \psi_h \right) \\
&= - \sum_{c \in \mathcal{C}} \sum_{S \in \mathcal{S}(c)} \int_{S} \psi_h \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \mathbf{n}^S \\
\int_{\Omega} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \cdot \mathbf{curl}_{\mathbf{x}} \psi_h &= - \sum_{S \in \mathcal{S}} \int_{S} [\psi_h \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \mathbf{n}^S].
\end{align*}
$$

In this last equality, $\begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \mathbf{n}^S$ is constant along the side, and so can be put out of both the jump and the integral. But then it remains the jump of $\psi_h$, which vanishes because $\psi_h$ is continuous. We have thus proven that $\mathbf{curl}_x Q_1(\Omega)$ is orthogonal to any $(\alpha, \beta)^T \in \mathbb{R}^2$, and this ends the proof of orthogonality.

For the moment, we have proven that

$$
\mathbb{R}^2 \oplus \nabla_{\mathbf{x}} (\mathbb{R}\mathit{aTu}(\Omega)) \oplus \mathbf{curl}_{\mathbf{x}} (Q_1(\Omega)) \subset \mathrm{d}S_h(\Omega).
$$

It remains to prove the equality, and for this, we rely on the dimension of the different approxima-
tion spaces.

First, we consider $\nabla_{\mathbf{x}} (\mathbb{R}\mathit{aTu}(\Omega))$. If $\varphi_h \in \mathbb{R}\mathit{aTu}(\Omega)$ is such that $\nabla_{\mathbf{x}} \varphi_h = 0$, then $\varphi_h$ is constant on each cell. Still, as $\mathbb{R}\mathit{aTu}(\Omega)$ is such that the average of the jump is zero along each face, $\varphi_h$ is actually constant on the whole domain. This means that $\dim (\ker \nabla_{\mathbf{x}} (\mathbb{R}\mathit{aTu}(\Omega))) = 1$, so that using (13):

$$
\dim (\nabla_{\mathbf{x}} (\mathbb{R}\mathit{aTu}(\Omega))) = \dim \mathbb{R}\mathit{aTu}(\Omega) - 1 = 2N^2 - 1.
$$

Doing the same for $\mathbf{curl}_x Q_1(\Omega)$, and using (11), we find

$$
\dim (\mathbf{curl}_x (Q_1(\Omega))) = \dim Q_1(\Omega) - 1 = N^2 - 1.
$$

This leads to

$$
\dim S_h(\Omega) - \dim (\nabla_{\mathbf{x}} (\mathbb{R}\mathit{aTu}(\Omega))) - \dim (\mathbf{curl}_{\mathbf{x}} (Q_1(\Omega))) - \dim \mathbb{R}^2 = 0,
$$

and proves the equality of dimensions, which ends the proof.

□
As $\psi_h \in \mathbb{Q}_1$, it is linear along each side, so that $\nabla_{\mathbf{x}}^{\perp} \psi_h \cdot \mathbf{n}^S$ is constant along each side. This leads to

$$ \int_{\Omega} (\nabla_{\mathbf{x}} \varphi_h) \cdot (\nabla_{\mathbf{x}}^{\perp} \psi_h) = \sum_{S \in \mathcal{S}} \nabla_{\mathbf{x}}^{\perp} (\psi_h) \cdot \mathbf{n}^S \int_S [[\varphi_h]]. $$

Last, by definition of the degrees of freedom of $\mathbb{R}aTu(\Omega)$, we have

$$ \forall S \in \mathcal{S} \quad \forall \varphi_h \in \mathbb{R}aTu(\Omega) \qquad \int_S [[\varphi_h]] = 0, $$

which proves that $\nabla_{\mathbf{x}} \varphi_h$ and $\nabla_{\mathbf{x}}^{\perp} \psi_h$ are orthogonal.

Following exactly the same lines, we can prove that any $(\alpha, \beta) \in \mathbb{R}^2$ is orthogonal to any element of $\nabla_{\mathbf{x}} \mathbb{R}aTu(\Omega)$.

It remains to deal with the orthogonality of $(\alpha, \beta) \in \mathbb{R}^2$ and any element of $\nabla_{\mathbf{x}}^{\perp} \mathbb{Q}_1(\Omega)$. We consider $\psi_h$ an element of $\mathbb{Q}_1(\Omega)$, and compute as follows

$$ \begin{aligned} \int_{\Omega} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \cdot \nabla_{\mathbf{x}}^{\perp} \psi_h &= \sum_{c \in \mathcal{C}} \int_c \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \cdot \nabla_{\mathbf{x}}^{\perp} \psi_h \\ &= - \sum_{c \in \mathcal{C}} \int_c \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \nabla_{\mathbf{x}} \psi_h \\ &= - \sum_{c \in \mathcal{C}} \int_c \mathrm{div}_{\mathbf{x}} \left( \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \psi_h \right) \\ &= - \sum_{c \in \mathcal{C}} \sum_{S \in \mathcal{S}(c)} \int_S \psi_h \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \mathbf{n}^S \end{aligned} $$ $$ \int_{\Omega} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \cdot \nabla_{\mathbf{x}}^{\perp} \psi_h = - \sum_{S \in \mathcal{S}} \int_S [[\psi_h \begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp}] \cdot \mathbf{n}^S]. $$

In this last equality, $\begin{pmatrix} \alpha \\ \beta \end{pmatrix}^{\perp} \cdot \mathbf{n}^S$ is constant along the side, and so can be put out of both the jump and the integral. But then it remains the jump of $\psi_h$, which vanishes because $\psi_h$ is continuous. We have thus proven that $\nabla_{\mathbf{x}}^{\perp} \mathbb{Q}_1(\Omega)$ is orthogonal to any $(\alpha, \beta)^T \in \mathbb{R}^2$, and this ends the proof of orthogonality.

For the moment, we have proven that

$$ \mathbb{R}^2 \oplus \nabla_{\mathbf{x}} (\mathbb{R}aTu(\Omega)) \oplus \nabla_{\mathbf{x}}^{\perp} (\mathbb{Q}_1(\Omega)) \subset d\mathbf{S}_h(\Omega). $$

It remains to prove the equality, and for this, we rely on the dimension of the different approximation spaces.

First, we consider $\nabla_{\mathbf{x}} (\mathbb{R}aTu(\Omega))$. If $\varphi_h \in \mathbb{R}aTu(\Omega)$ is such that $\nabla_{\mathbf{x}} \varphi_h = 0$, then $\varphi_h$ is constant on each cell. Still, as $\mathbb{R}aTu(\Omega)$ is such that the average of the jump is zero along each face, $\varphi_h$ is actually constant on the whole domain. This means that $\dim (\ker \nabla_{\mathbf{x}} (\mathbb{R}aTu(\Omega))) = 1$, so that using (13):

$$ \dim (\nabla_{\mathbf{x}} (\mathbb{R}aTu(\Omega))) = \dim \mathbb{R}aTu(\Omega) - 1 = 2N^2 - 1. $$

Doing the same for $\nabla_{\mathbf{x}}^{\perp} \mathbb{Q}_1(\Omega)$, and using (11), we find

$$ \dim (\nabla_{\mathbf{x}}^{\perp} (\mathbb{Q}_1(\Omega))) = \dim \mathbb{Q}_1(\Omega) - 1 = N^2 - 1. $$

This leads to

$$ \dim \mathbf{S}_h(\Omega) - \dim (\nabla_{\mathbf{x}} (\mathbb{R}aTu(\Omega))) - \dim (\nabla_{\mathbf{x}}^{\perp} (\mathbb{Q}_1(\Omega))) - \dim \mathbb{R}^2 = 0, $$

and proves the equality of dimensions, which ends the proof. $\square$

**Corollary 2.1.** For all $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$, we have

$$ \mathcal{P}_h^{\varphi} [\mathbf{u}_h] = 0 \quad \Leftrightarrow \quad \forall S \in \mathcal{S}, \quad [[\mathbf{u}_h \cdot \mathbf{n}^S]] = 0. $$

*Proof.* If $\mathcal{P}_h^{\varphi} [\mathbf{u}_h] = 0$, from Proposition 2 we can write

$$ \mathbf{u}_h = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} + \nabla_{\mathbf{x}}^{\perp} \psi $$
**Corollary 2.1.** For all $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$, we have

$$ \mathcal{P}_h^\varphi[\mathbf{u}_h] = 0 \quad \Leftrightarrow \quad \forall S \in \mathcal{S}, \quad \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0. $$

*Proof.* If $\mathcal{P}_h^\varphi[\mathbf{u}_h] = 0$, from Proposition 2 we can write

$$ \mathbf{u}_h = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} + \mathrm{curl}_{\mathbf{x}} \psi $$

for some $(\alpha, \beta) \in \mathbb{R}^2$ and $\psi \in \mathbb{Q}_1(\Omega)$. Moreover, as recalled in the proof of Proposition 2, we have for all $S \in \mathcal{S}$, $\llbracket \mathrm{curl}_{\mathbf{x}} \psi \cdot \mathbf{n}^S \rrbracket = 0$ and then $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0$.

We now assume that $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$ is such that for all $S \in \mathcal{S}$, $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0$. From Proposition 2, we know that $(\alpha, \beta) \in \mathbb{R}^2$, $\varphi \in \mathbb{R}aTu(\Omega)$ and $\psi \in \mathbb{Q}_1(\Omega)$ exist such that

$$ \mathbf{u}_h = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} + \nabla_{\mathbf{x}}\varphi + \mathrm{curl}_{\mathbf{x}} \psi. $$

As $\llbracket \mathrm{curl}_{\mathbf{x}} \psi \cdot \mathbf{n} \rrbracket = \llbracket (\alpha, \beta)^T \cdot \mathbf{n} \rrbracket = 0$, $\llbracket \nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} \rrbracket = \llbracket \mathbf{u}_h \cdot \mathbf{n} \rrbracket = 0$. Then

$$ \begin{aligned} \int_{\Omega} \|\nabla_{\mathbf{x}}\varphi\|^2 &= \sum_{c \in \mathcal{C}} \int_{c} \nabla_{\mathbf{x}}\varphi \cdot \nabla_{\mathbf{x}}\varphi \\ &= \sum_{c \in \mathcal{C}} \left( \int_{c} \mathrm{div}_{\mathbf{x}}(\varphi\nabla_{\mathbf{x}}\varphi) - \int_{c} \varphi\mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}\varphi) \right) \\ &= \sum_{c \in \mathcal{C}} \int_{c} \varphi\nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} - \sum_{c \in \mathcal{C}} \int_{c} \varphi\mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}\varphi) \\ \int_{\Omega} \|\nabla_{\mathbf{x}}\varphi\|^2 &= \sum_{S \in \mathcal{S}} \int_{S} \llbracket \varphi\nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} \rrbracket - \sum_{c \in \mathcal{C}} \int_{c} \varphi\mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}\varphi). \end{aligned} $$

The second sum is zero, because as stated in Proposition 1, the divergence of all the basis functions is zero. Concerning the first sum, still following Proposition 1, $\nabla_{\mathbf{x}}\varphi \cdot \mathbf{n}$ is constant on each side, which gives

$$ \int_{\Omega} \nabla_{\mathbf{x}}\varphi \cdot \nabla_{\mathbf{x}}\varphi = \sum_{S \in \mathcal{S}} \nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} \int_{S} \llbracket \varphi \rrbracket, $$

and by hypothesis, this last integral vanishes on all sides, so that

$$ \int_{\Omega} \|\nabla_{\mathbf{x}}\varphi\|^2 = 0, $$

which ends the proof.

□

## 2.4 First order wave system

In this section, a numerical discretization of the first order waves system (1), based on the new approximation space $d\mathbf{S}_h$, is proposed. Then, it is proved that this numerical scheme, coupled with a Godunov numerical flux, preserves the divergence free component of the Hodge-Helmholtz decomposition and preserves the adjoint curl.

### 2.4.1 Discretization

The discontinuous Galerkin discretization of system (1) on a Cartesian periodic mesh consists in finding $(p_h, \mathbf{u}_h)$ in some approximation space $V_h = V_h \times V_h$ such that for all $v \in V_h$ and all $\mathbf{w} \in V_h$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_{c} v \partial_{\tau} p_h - \sum_{c \in \mathcal{C}} \int_{c} \frac{1}{\rho_0} \mathbf{u}_h \cdot \nabla v + \sum_{S \in \mathcal{S}} \int_{S} \llbracket v \rrbracket \left( \frac{1}{\rho_0} \left\{ \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket + d_{11} \llbracket p_h \rrbracket \right\} \right) = 0, \\ \sum_{c \in \mathcal{C}} \int_{c} \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h - \sum_{c \in \mathcal{C}} \int_{c} \kappa_0 p_h \mathrm{div}_{\mathbf{x}} \mathbf{w} + \sum_{S \in \mathcal{S}} \int_{S} \llbracket \mathbf{w} \rrbracket \cdot \left( \kappa_0 \left\{ \llbracket p_h \mathbf{n}^S \rrbracket + \llbracket \mathbf{D}_{22}(\mathbf{n}^S) \mathbf{u}_h \rrbracket \right\} \right) = 0 \end{array} \right. \quad (16) $$

where

$$ \mathbf{D}(\mathbf{n}) = \begin{pmatrix} d_{11} & 0 \\ 0 & \mathbf{D}_{22}(\mathbf{n}) \end{pmatrix} \quad (17) $$

is the stabilization matrix. We recall that Godunov stabilization is given by

$$ \mathbf{D}^{\text{Godunov}}(\mathbf{n}) = \frac{c_0}{2} \begin{pmatrix} 1 & 0 \\ 0 & \mathbf{n}\mathbf{n}^T \end{pmatrix} \quad (18) $$
for some $(\alpha, \beta) \in \mathbb{R}^2$ and $\psi \in \mathbb{Q}_1(\Omega)$. Moreover, from the proof Proposition 2, we have for all $S \in \mathcal{S}$, $\llbracket \nabla_{\mathbf{x}}^{\perp} \psi \cdot \mathbf{n}^S \rrbracket = 0$ and then $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0$.

We now assume that $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$ is such that for all $S \in \mathcal{S}$, $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0$. From Proposition 2, we know that $(\alpha, \beta) \in \mathbb{R}^2$, $\varphi \in \mathbb{R}aTu(\Omega)$ and $\psi \in \mathbb{Q}_1(\Omega)$ exist such that

$$ \mathbf{u}_h = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} + \nabla_{\mathbf{x}}\varphi + \nabla_{\mathbf{x}}^{\perp}\psi. $$

As $\llbracket \nabla_{\mathbf{x}}^{\perp} \psi \cdot \mathbf{n} \rrbracket = \llbracket (\alpha, \beta)^T \cdot \mathbf{n} \rrbracket = 0$, $\llbracket \nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} \rrbracket = \llbracket \mathbf{u}_h \cdot \mathbf{n} \rrbracket = 0$. Then

$$ \begin{aligned} \int_{\Omega} \|\nabla_{\mathbf{x}}\varphi\|^2 &= \sum_{c \in \mathcal{C}} \int_{c} \nabla_{\mathbf{x}}\varphi \cdot \nabla_{\mathbf{x}}\varphi \\ &= \sum_{c \in \mathcal{C}} \left( \int_{c} \mathrm{div}_{\mathbf{x}}(\varphi\nabla_{\mathbf{x}}\varphi) - \int_{c} \varphi\mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}\varphi) \right) \\ &= \sum_{c \in \mathcal{C}} \int_{\mathcal{S}(c)} \varphi\nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} - \sum_{c \in \mathcal{C}} \int_{c} \varphi\mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}\varphi) \\ \int_{\Omega} \|\nabla_{\mathbf{x}}\varphi\|^2 &= \sum_{S \in \mathcal{S}} \int_{S} \llbracket \varphi\nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} \rrbracket - \sum_{c \in \mathcal{C}} \int_{c} \varphi\mathrm{div}_{\mathbf{x}}(\nabla_{\mathbf{x}}\varphi). \end{aligned} $$

The second sum is zero, because as stated in Proposition 1, the divergence of all the basis functions is zero. Concerning the first sum, still following Proposition 1, $\nabla_{\mathbf{x}}\varphi \cdot \mathbf{n}$ is constant on each side, which gives

$$ \int_{\Omega} \nabla_{\mathbf{x}}\varphi \cdot \nabla_{\mathbf{x}}\varphi = \sum_{S \in \mathcal{S}} \nabla_{\mathbf{x}}\varphi \cdot \mathbf{n} \int_{S} \llbracket \varphi \rrbracket, $$

and by hypothesis, this last integral vanishes on all sides, so that

$$ \int_{\Omega} \|\nabla_{\mathbf{x}}\varphi\|^2 = 0, $$

which ends the proof.

□

## 2.4 First order wave system

In this section, a numerical discretization of the first order waves system (1), based on the new approximation space $d\mathbf{S}_h$, is proposed. Then, it is proved that this numerical scheme, coupled with a Godunov numerical flux, preserves the divergence free component of the Hodge-Helmholtz decomposition and preserves the adjoint curl.

### 2.4.1 Discretization

The discontinuous Galerkin discretization of system (1) on a Cartesian periodic mesh consists in finding $(p_h, \mathbf{u}_h)$ in some approximation space $V_h = V_h \times V_h$ such that for all $v \in V_h$ and all $\mathbf{w} \in V_h$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_{c} v \partial_{\tau} p_h - \sum_{c \in \mathcal{C}} \int_{c} \frac{1}{\rho_0} \mathbf{u}_h \cdot \nabla v + \sum_{S \in \mathcal{S}} \int_{S} \llbracket v \rrbracket \left( \frac{1}{\rho_0} \{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \} + d_{11} \llbracket p_h \rrbracket \right) = 0, \\ \sum_{c \in \mathcal{C}} \int_{c} \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h - \sum_{c \in \mathcal{C}} \int_{c} \kappa_0 p_h \mathrm{div}_{\mathbf{x}} \mathbf{w} + \sum_{S \in \mathcal{S}} \int_{S} \llbracket \mathbf{w} \rrbracket \cdot (\kappa_0 \{ \{ p_h \mathbf{n}^S \} \} + \llbracket \mathbf{D}_{22}(\mathbf{n}^S) \mathbf{u}_h \rrbracket) = 0 \end{array} \right. \quad (16) $$

where

$$ \mathbf{D}(\mathbf{n}) = \begin{pmatrix} d_{11} & 0 \\ 0 & \mathbf{D}_{22}(\mathbf{n}) \end{pmatrix} \quad (17) $$

is the stabilization matrix. We recall that Godunov stabilization is given by

$$ \mathbf{D}^{\text{Godunov}}(\mathbf{n}) = \frac{c_0}{2} \begin{pmatrix} 1 & 0 \\ 0 & \mathbf{n}\mathbf{n}^T \end{pmatrix} \quad (18) $$

while Rusanov stabilization is given by

$$ \mathbf{D}^{\text{Rusanov}}(\mathbf{n}) = \frac{c_0}{2} \begin{pmatrix} 1 & 0 \\ 0 & I_2 \end{pmatrix}. \quad (19) $$

As a consequence of the Proposition 1, using the approximation space $V_h(\Omega) = d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$, the cell integrals in (16) vanish and discretization (16) can be seen as a finite volume
while Rusanov stabilization is given by

$$ \mathbf{D}^{\text{Rusanov}}(\mathbf{n}) = \frac{c_0}{2} \begin{pmatrix} 1 & 0 \\ 0 & I_2 \end{pmatrix}. \qquad (19) $$

As a consequence of the Proposition 1, using the approximation space $V_h(\Omega) = d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$, the cell integrals in (16) vanish and discretization (16) can be seen as a finite volume scheme that consists in finding $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ such that for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h(\Omega)$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_{c} v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}} \int_{S} \llbracket v \rrbracket \left( \frac{1}{\rho_0} \{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \} + d_{11} \llbracket p_h \rrbracket \right) = 0, \\ \\ \sum_{c \in \mathcal{C}} \int_{c} \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_{S} \llbracket \mathbf{w} \rrbracket \cdot ( \kappa_0 \{ \{ p_h \mathbf{n}^S \} \} + \llbracket \mathbf{D}_{22}(\mathbf{n}^S) \mathbf{u}_h \rrbracket ) = 0. \end{array} \right. \qquad (20) $$

Note that, in the case of a Cartesian mesh, the mass matrix is diagonal.

### 2.4.2 Godunov stabilization and stationary solution

**Proposition 3.** Using the discretization (20) of the wave system with Godunov stabilization (18), a state $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0 \times d\mathbf{S}_h$ is a stationary solution of (20) if and only if $p_h$ is uniform and for all side $S \in \mathcal{S}$, $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0$ (divergence free velocity field).

The following proof is an adaptation of the proof performed in [21] for triangular mesh. It can also be found in [22].

*Proof.* If $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0 \times d\mathbf{S}_h$ is a stationary solution of (20) with Godunov stabilization (18), it satisfies for all $v \in d\mathbb{P}_0$ and all $\mathbf{w} \in d\mathbf{S}_h$,

$$ \left\{ \begin{array}{l} \sum_{S \in \mathcal{S}} \int_{S} \llbracket v \rrbracket \left( \frac{1}{\rho_0} \{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \} + \frac{c_0}{2} \llbracket p_h \rrbracket \right) = 0, \\ \\ \sum_{S \in \mathcal{S}} \int_{S} \llbracket \mathbf{w} \rrbracket \cdot \left( \kappa_0 \{ \{ p_h \mathbf{n}^S \} \} + \frac{c_0}{2} \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket \mathbf{n}^S \right) = 0. \end{array} \right. \qquad (21) $$

Moreover, for $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0 \times d\mathbf{S}_h$, $v \in d\mathbb{P}_0$, and $\mathbf{w} \in d\mathbf{S}_h$, we have

$$ \forall c \in \mathcal{C}, \quad 0 = \int_c \frac{1}{\rho_0} \mathrm{div}_x(\mathbf{u}_h v) = \sum_{S \in \mathcal{S}(c)} \int_S \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S v $$

and

$$ \forall c \in \mathcal{C}, \quad 0 = \int_c \kappa_0 \mathrm{div}_x(p_h \mathbf{w}) = \sum_{S \in \mathcal{S}(c)} \int_S \kappa_0 p_h \mathbf{w} \cdot \mathbf{n}^S. $$

Using test functions vanishing everywhere except on one cell $c \in \mathcal{C}$, and subtracting the above equalities from (21) gives

$$ \forall c \in \mathcal{C}, \quad \forall \mathbf{w} \in \hat{\mathbf{S}}_0(c), \quad \left\{ \begin{array}{l} \sum_{S \in \mathcal{S}(c)} \int_S \frac{1}{2} \left( -\frac{1}{\rho_0} \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket + c_0 \llbracket p_h \rrbracket \right) = 0, \\ \\ \sum_{S \in \mathcal{S}(c)} \int_S \frac{\mathbf{w} \cdot \mathbf{n}^S}{2} \left( -\kappa_0 \llbracket p_h \rrbracket + c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket \right) = 0. \end{array} \right. \qquad (22) $$

Last, multiplying the first line of (22) by $\rho_0 c_0$ and using $\kappa_0 = \rho_0 c_0^2$, and multiplying the second line of (22) by $-1$ lead to

$$ \forall c \in \mathcal{C}, \quad \forall \mathbf{w} \in \hat{\mathbf{S}}_0(c), \quad \left\{ \begin{array}{l} \sum_{S \in \mathcal{S}(c)} \int_S \frac{1}{2} (\kappa_0 \llbracket p_h \rrbracket - c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket) = 0, \\ \\ \sum_{S \in \mathcal{S}(c)} \int_S \frac{\mathbf{w} \cdot \mathbf{n}^S}{2} (\kappa_0 \llbracket p_h \rrbracket - c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket) = 0. \end{array} \right. \qquad (23) $$

From Proposition 1, $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket$ is constant along each side $S \in \mathcal{S}$, and so

$$ \gamma_S := \kappa_0 \llbracket p_h \rrbracket - c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket, $$
scheme that consists in finding $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ such that for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h(\Omega)$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_{c} v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}} \int_{S} \llbracket v \rrbracket \left( \frac{1}{\rho_0} \{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \} + d_{11} \llbracket p_h \rrbracket \right) = 0, \\ \\ \sum_{c \in \mathcal{C}} \int_{c} \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_{S} \llbracket \mathbf{w} \rrbracket \cdot (\kappa_0 \{ \{ p_h \mathbf{n}^S \} \} + \llbracket \mathbf{D}_{22}(\mathbf{n}^S) \mathbf{u}_h \rrbracket) = 0. \end{array} \right. \qquad (20) $$

Note that, in the case of a Cartesian mesh, the mass matrix is diagonal.

### 2.4.2 Godunov stabilization and stationary solution

**Proposition 3.** Using the discretization (20) of the wave system with Godunov stabilization (18), a state $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0 \times d\mathbf{S}_h$ is a stationary solution of (20) if and only if $p_h$ is uniform and for all side $S \in \mathcal{S}$, $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket = 0$ (divergence free velocity field).

The following proof is an adaptation of the proof performed in [20] for triangular mesh. It can also be found in [21].

*Proof.* If $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0 \times d\mathbf{S}_h$ is a stationary solution of (20) with Godunov stabilization (18), it satisfies for all $v \in d\mathbb{P}_0$ and all $\mathbf{w} \in d\mathbf{S}_h$,

$$ \left\{ \begin{array}{l} \sum_{S \in \mathcal{S}} \int_S \llbracket v \rrbracket \left( \frac{1}{\rho_0} \{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \} + \frac{c_0}{2} \llbracket p_h \rrbracket \right) = 0, \\ \\ \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{w} \rrbracket \cdot \left( \kappa_0 \{ \{ p_h \mathbf{n}^S \} \} + \frac{c_0}{2} \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket \mathbf{n}^S \right) = 0. \end{array} \right. \qquad (21) $$

Moreover, for $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0 \times d\mathbf{S}_h$, $v \in d\mathbb{P}_0$, and $\mathbf{w} \in d\mathbf{S}_h$, we have

$$ \forall c \in \mathcal{C}, \quad 0 = \int_c \frac{1}{\rho_0} \mathrm{div}_x (\mathbf{u}_h v) = \sum_{S \in \mathcal{S}(c)} \int_S \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S v $$

and

$$ \forall c \in \mathcal{C}, \quad 0 = \int_c \kappa_0 \mathrm{div}_x (p_h \mathbf{w}) = \sum_{S \in \mathcal{S}(c)} \int_S \kappa_0 p_h \mathbf{w} \cdot \mathbf{n}^S. $$

Using test functions vanishing everywhere except on one cell $c \in \mathcal{C}$, and subtracting the above equalities from (21) gives

$$ \forall c \in \mathcal{C}, \quad \forall \mathbf{w} \in \hat{\mathbf{S}}_0(c), \quad \left\{ \begin{array}{l} \sum_{S \in \mathcal{S}(c)} \int_S \frac{1}{2} \left( -\frac{1}{\rho_0} \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket + c_0 \llbracket p_h \rrbracket \right) = 0, \\ \\ \sum_{S \in \mathcal{S}(c)} \int_S \frac{\mathbf{w} \cdot \mathbf{n}^S}{2} \left( -\kappa_0 \llbracket p_h \rrbracket + c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket \right) = 0. \end{array} \right. $$

and so

$$ \forall c \in \mathcal{C}, \quad \forall \mathbf{w} \in \hat{\mathbf{S}}_0(c), \quad \left\{ \begin{array}{l} \sum_{S \in \mathcal{S}(c)} \int_S \frac{1}{2} (\kappa_0 \llbracket p_h \rrbracket - c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket) = 0, \\ \\ \sum_{S \in \mathcal{S}(c)} \int_S \frac{\mathbf{w} \cdot \mathbf{n}^S}{2} (\kappa_0 \llbracket p_h \rrbracket - c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket) = 0. \end{array} \right. \qquad (22) $$

From Proposition 1, $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket$ is constant along each side $S \in \mathcal{S}$, and so

$$ \gamma_S := \kappa_0 \llbracket p_h \rrbracket - c_0 \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket, $$

is also constant on each side $S \in \mathcal{S}$. For a square cell $c \in \mathcal{C}$, denoting by $(\gamma_R, \gamma_T, \gamma_L, \gamma_B) \in \mathbb{R}^4$ the unknown on the right, top, left and bottom faces of the cell $c$, (22) leads to the following linear system to solve

$$ \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \\ 1 & -1 & 1 & -1 \end{bmatrix} \begin{bmatrix} \gamma_R \\ \gamma_T \\ \gamma_L \\ \gamma_B \end{bmatrix} = 0. $$
is also constant on each side $S \in \mathcal{S}$. For a square cell $c \in \mathcal{C}$, denoting by $(\gamma_R, \gamma_T, \gamma_L, \gamma_B) \in \mathbb{R}^4$ the unknown on the right, top, left and bottom faces of the cell $c$, (23) leads to the following linear system to solve

$$ \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \\ 1 & -1 & 1 & -1 \end{bmatrix} \begin{bmatrix} \gamma_R \\ \gamma_T \\ \gamma_L \\ \gamma_B \end{bmatrix} = 0. $$

Since the matrix is invertible, we get $\gamma_R = \gamma_T = \gamma_L = \gamma_B = 0$ and then

$$ \forall c \in \mathcal{C}, \quad \forall S \in \mathcal{S}(c), \quad \kappa_0[\mathbf{p}_h] - c_0[\mathbf{u}_h \cdot \mathbf{n}^S] = 0. $$

If $S \in \mathcal{S}$ is a side, denoting by $c'$ the neighboring cell to $c$ with respect to the normal $\mathbf{n}^S$ and performing the same study for cell $c'$ we obtain

$$ \kappa_0[\mathbf{p}_h] + c_0[\mathbf{u}_h \cdot \mathbf{n}^S] = 0. $$

Then, for all side $S \in \mathcal{S}$, we obtain $[\mathbf{p}_h] = 0 = [\mathbf{u}_h \cdot \mathbf{n}^S]$ and the result is proven. The reciprocal is obvious. $\square$

### 2.4.3 Conservation of the adjoint curl with the Godunov' flux

In Proposition 2, the following map is considered:

$$ \begin{array}{rcl} \mathrm{curl}_{\mathbf{x}} : & \mathbb{Q}_1 & \mapsto & \mathrm{d}\mathbf{S}_h(\Omega) \\ & \psi & \mapsto & \mathrm{curl}_{\mathbf{x}} \psi. \end{array} $$

In this subsection, we consider the adjoint of this operator defined as

$$ \forall \mathbf{u} \in \mathrm{d}\mathbf{S}_h(\Omega) \quad \forall \psi \in \mathbb{Q}_1 \quad \int_{\Omega} \psi (\mathrm{curl}_{\mathbf{x}})^* \mathbf{u} = \int_{\Omega} \mathbf{u} \cdot \mathrm{curl}_{\mathbf{x}} \psi. \quad (24) $$

From a practical point of view, $(\mathrm{curl}_{\mathbf{x}})^*$ can be computed by inverting the $\mathbb{Q}_1$ mass matrix. As the operator $\mathrm{rot}_{\mathbf{x}}$ is the opposite of the adjoint of $\mathrm{curl}_{\mathbf{x}}$, $(\mathrm{curl}_{\mathbf{x}})^*$ can also be seen as the approximation of the opposite of $\mathrm{rot}_{\mathbf{x}}$.

**Proposition 4** (Conservation of the adjoint curl). Consider the numerical scheme (20) of the wave system (1) with Godunov' stabilization (18). Then

$$ \forall \psi \in \mathbb{Q}_1 \quad \int_{\Omega} \mathrm{curl}_{\mathbf{x}} \psi \cdot \partial_{\tau} \mathbf{u}_h = 0, $$

which means that $\partial_{\tau} ((\mathrm{curl}_{\mathbf{x}})^* \mathbf{u}_h) = 0$.

*Proof.* Taking the equation on the velocity evolution of (20) gives

$$ \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_S [\mathbf{w}] \cdot (\kappa_0 \{\{ p_h \mathbf{n}^S \}\} + [\mathbf{D}_{22}(\mathbf{n}^S) \mathbf{u}_h ] ) = 0. $$

With the Godunov' numerical flux, we have $\mathbf{D}_{22}(\mathbf{n}^S) = \frac{c_0}{2} \mathbf{n}^S (\mathbf{n}^S)^T$, which gives

$$ \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_S [\mathbf{w}] \cdot (\kappa_0 \{\{ p_h \mathbf{n}^S \}\} + \frac{c_0}{2} [\mathbf{n}^S]^T \mathbf{u}_h ] \mathbf{n}^S = 0, $$

which can be rewritten

$$ \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_S [\mathbf{w} \cdot \mathbf{n}^S] (\kappa_0 \{\{ p_h \} \} + \frac{c_0}{2} [\mathbf{u}_h \cdot \mathbf{n}^S]) = 0. $$

Suppose now that $\mathbf{w} = \mathrm{curl}_{\mathbf{x}} \psi$ for $\psi \in \mathbb{Q}_1$. As $\psi$ is linear on each side, its gradient is continuous in the side direction, which means that $[\mathrm{curl}_{\mathbf{x}} \psi \cdot \mathbf{n}^S] = 0$ on each side, which ends the proof. $\square$

Note that from an implementation point of view, the quadrature formula used for the mass matrix should match with the one used for computing the discrete operator $(\mathrm{curl}_{\mathbf{x}})^*$ from (24).
Since the matrix is invertible, we get $\gamma_R = \gamma_T = \gamma_L = \gamma_B = 0$ and then

$$ \forall c \in \mathcal{C}, \quad \forall S \in \mathcal{S}(c), \quad \kappa_0[\llbracket p_h \rrbracket] - c_0[\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket] = 0. $$

If $S \in \mathcal{S}$ is a side, denoting by $c'$ the neighboring cell to $c$ with respect to the normal $\mathbf{n}^S$ and performing the same study for cell $c'$ we obtain

$$ \kappa_0[\llbracket p_h \rrbracket] + c_0[\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket] = 0. $$

Then, for all side $S \in \mathcal{S}$, we obtain $[\llbracket p_h \rrbracket] = 0 = [\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket]$ and the result is proven. The reciprocal is obvious. $\square$

### 2.4.3 Conservation of the adjoint curl with the Godunov' flux

In Proposition 2, the following map is considered:

$$ \begin{aligned} \nabla^{\perp} : \quad & \mathbb{Q}_1 \quad \mapsto \quad \mathrm{d}\mathbf{S}_h(\Omega) \\ & \psi \quad \mapsto \quad \nabla^{\perp}\psi. \end{aligned} $$

In this subsection, we consider the adjoint of this operator defined as

$$ \forall \mathbf{u} \in \mathrm{d}\mathbf{S}_h(\Omega) \quad \forall \psi \in \mathbb{Q}_1 \quad \int_{\Omega} \psi (\nabla^{\perp})^{*} \mathbf{u} = \int_{\Omega} \mathbf{u} \cdot \nabla^{\perp} \psi. \quad (23) $$

From a practical point of view, $(\nabla^{\perp})^*$ can be computed by inverting the $\mathbb{Q}_1$ mass matrix. As the curl is a self-adjoint operator, $(\nabla^{\perp})^*$ can be seen also as a discrete curl.

**Proposition 4 (Conservation of the adjoint curl).** Consider the numerical scheme (20) of the wave system (1) with Godunov' stabilization (18). Then

$$ \forall \psi \in \mathbb{Q}_1 \quad \int_{\Omega} \nabla^{\perp} \psi \cdot \partial_{\tau} \mathbf{u}_h = 0, $$

which means that $\partial_{\tau} ((\nabla^{\perp})^* \mathbf{u}_h) = 0$.

*Proof.* Taking the equation on the velocity evolution of (20) gives

$$ \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_S [\llbracket \mathbf{w} \rrbracket] \cdot (\kappa_0 \{\{ p_h \mathbf{n}^S \}\} + [\llbracket \mathbf{D}_{22}(\mathbf{n}^S) \mathbf{u}_h \rrbracket]) = 0. $$

With the Godunov' numerical flux, we have $\mathbf{D}_{22}(\mathbf{n}^S) = \frac{c_0}{2} \mathbf{n}^S (\mathbf{n}^S)^T$, which gives

$$ \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_S [\llbracket \mathbf{w} \rrbracket] \cdot (\kappa_0 \{\{ p_h \mathbf{n}^S \}\} + \frac{c_0}{2} [\llbracket (\mathbf{n}^S)^T \mathbf{u}_h \rrbracket] \mathbf{n}^S) = 0, $$

which can be rewritten

$$ \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}} \int_S [\llbracket \mathbf{w} \cdot \mathbf{n}^S \rrbracket] (\kappa_0 \{\{ p_h \} \} + \frac{c_0}{2} [\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket]) = 0. $$

Suppose now that $\mathbf{w} = \nabla^{\perp}\psi$ for $\psi \in \mathbb{Q}_1$. As $\psi$ is linear on each side, its gradient is continuous in the side direction, which means that $[\nabla^{\perp}\psi \cdot \mathbf{n}^S] = 0$ on each side, which ends the proof. $\square$

Note that from an implementation point of view, the quadrature formula used for the mass matrix should match with the one used for computing the discrete operator $(\nabla^{\perp})^*$ from (23).

## 3 Quadrangular case with boundary conditions

In this section, we still consider the wave system (1), but the domain $\Omega$ is no more periodic, and we consider two types of boundary conditions:

* Inlet/Outlet boundary conditions, in which a pressure $p_b$ and a velocity $\mathbf{u}_b$ are weakly imposed

* Wall boundary conditions, in which $\mathbf{u} \cdot \mathbf{n} = 0$ is weakly imposed, and the imposed pressure flux is 0.
### 3 Quadrangular case with boundary conditions

In this section, we still consider the wave system (1), but the domain $\Omega$ is no more periodic, and we consider two types of boundary conditions:

* Inlet/Outlet boundary conditions, in which a pressure $p_b$ and a velocity $\mathbf{u}_b$ are weakly imposed

* Wall boundary conditions, in which $\mathbf{u} \cdot \mathbf{n} = 0$ is weakly imposed, and the imposed pressure flux is 0.

We suppose that $p_b$ is uniform, and that $\int_{\partial\Omega_{\text{inlet}/\text{outlet}}} \mathbf{u}_b \cdot \mathbf{n} = 0$, which is necessary for ensuring a long time limit of the wave system.

When dealing with the bounded case, it is important to understand that (1) does not always preserve the curl of the velocity, because of boundary conditions. The structure that is preserved is a bit more complicated, and relies on a special Hodge-Helmholtz decomposition that was proposed in [27], and is recalled here

**Proposition 5** (Hodge-Helmholtz decomposition adapted to the wave system, see [27, Prop. 1 & 2]). Consider the problem (1), with inlet/oulet boundary conditions ($p_b, \mathbf{u}_b$) or wall boundary conditions, and with an initial pressure $p_0$ and velocity $\mathbf{u}_0$. For any $\mathbf{u}$, consider $\varphi \in H^1(\Omega)$ a solution of the variational problem

$$ \forall g \in H^1(\Omega) \qquad \int_{\Omega} \nabla g \cdot \nabla \varphi = \int_{\Omega} \mathbf{u} \cdot \nabla g - \int_{\Omega} g \mathbf{u}_b \cdot \mathbf{n}, $$

and which is unique up to a constant, and define $\mathbf{u}^{\varphi} = \nabla \varphi$, and $\mathbf{u}^{\psi} = \mathbf{u} - \mathbf{u}^{\varphi}$. Then the Hodge-Helmholtz decomposition

$$ \mathbf{u} = \mathbf{u}^{\varphi} + \mathbf{u}^{\psi}, $$

is adapted to the wave system in the sense that the component $\mathbf{u}^{\psi}$ (divergence free component) is constant in time. Moreover, the long time limit of the wave system exists and is given by $(p_b, \mathbf{u}^{\psi}(0))$ where $\mathbf{u}^{\psi}(0)$ corresponds to the divergence free component of the initial velocity field $\mathbf{u}_0$.

Note that in Proposition 5, due to boundary conditions, the set of all the $\mathbf{u}^{\varphi}$ and of all the $\mathbf{u}^{\psi}$ are not vectorial spaces, but affine spaces, which means that the decomposition of 0 provided by Proposition 5 is not 0 + 0.

The aim of this section is to develop a discrete counterpart of the structure depicted in Proposition 5 on quadrangular meshes.

#### 3.1 Discretization and boundary conditions

We denote by $M_h$ a conformal quadrangular mesh, by $\mathcal{P}$ the set of points of $M_h$, $\mathcal{C}$ the set of cells of $M_h$, by $\mathcal{S}_i$ the set of interior faces of $M_h$, and by $\mathcal{S}_b$ the set of boundary faces. Boundary faces are supposed to be oriented, such that the normal is outgoing. Each interior side $S$ is arbitrarily oriented, and denoting by $\mathbf{n}^S$ its normal, we will call the left cell the one from which $\mathbf{n}^S$ is outgoing, and the right cell the one in which $\mathbf{n}^S$ is ingoing. The discontinuous Galerkin discretization of the wave system (1) consists in finding $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ such that for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h(\Omega)$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_c v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S [v] \left( \frac{1}{\rho_0} \left\{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \right\} + d_{11} [p_h] \right) \\ \qquad + \sum_{S \in \mathcal{S}_b} \int_S v^L \left[ \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \right]_b = 0, \\[1.5em] \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{w}] \cdot (\kappa_0 \left\{ p_h \mathbf{n}^S \right\} + [\mathbf{D}_{22}(\mathbf{n}^S)\mathbf{u}_h]) \\ \qquad + \sum_{S \in \mathcal{S}_b} \int_S \mathbf{w}^L \cdot [\kappa_0 p_h \mathbf{n}^S]_b = 0 \end{array} \right. \quad (25) $$

where the boundary fluxes in (25) are given by

$$ \left[ \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \right]_{\text{wall}} = \left( \begin{array}{c} 0 \\ \kappa_0 p_h \cdot \mathbf{n}^S + c_0 (\mathbf{u}_h \cdot \mathbf{n}^S) \mathbf{n}^S \end{array} \right) \quad (26) $$
We suppose that $p_b$ is uniform, and that $\int_{\partial\Omega_{\text{inlet/outlet}}} \mathbf{u}_b \cdot \mathbf{n} = 0$, which is necessary for ensuring a long time limit of the wave system.

When dealing with the bounded case, it is important to understand that (1) does not always preserve the curl of the velocity, because of boundary conditions. The structure that is preserved is a bit more complicated, and relies on a special Hodge-Helmholtz decomposition that was proposed in [25], and is recalled here

**Proposition 5** (Hodge-Helmholtz decomposition adapted to the wave system, see [25, Prop. 1 & 2]). Consider the problem (1), with inlet/oulet boundary conditions ($p_b, \mathbf{u}_b$) or wall boundary conditions, and with an initial pressure $p_0$ and velocity $\mathbf{u}_0$. For any $\mathbf{u}$, consider $\varphi \in H^1(\Omega)$ a solution of the variational problem

$$ \forall g \in H^1(\Omega) \qquad \int_{\Omega} \nabla g \cdot \nabla \varphi = \int_{\Omega} \mathbf{u} \cdot \nabla g - \int_{\Omega} g \mathbf{u}_b \cdot \mathbf{n}, $$

and which is unique up to a constant, and define $\mathbf{u}^{\varphi} = \nabla \varphi$, and $\mathbf{u}^{\psi} = \mathbf{u} - \mathbf{u}^{\varphi}$. Then the Hodge-Helmholtz decomposition

$$ \mathbf{u} = \mathbf{u}^{\varphi} + \mathbf{u}^{\psi}, $$

is adapted to the wave system in the sense that the component $\mathbf{u}^{\psi}$ (divergence free component) is constant in time. Moreover, the long time limit of the wave system exists and is given by $(p_b, \mathbf{u}^{\psi}(0))$ where $\mathbf{u}^{\psi}(0)$ corresponds to the divergence free component of the initial velocity field $\mathbf{u}_0$.

Note that in Proposition 5, due to boundary conditions, the set of all the $\mathbf{u}^{\varphi}$ and of all the $\mathbf{u}^{\psi}$ are not vectorial spaces, but affine spaces, which means that the decomposition of 0 provided by Proposition 5 is not 0 + 0.

The aim of this section is to develop a discrete counterpart of the structure depicted in Proposition 5 on quadrangular meshes.

## 3.1 Discretization and boundary conditions

We denote by $\mathcal{M}_h$ a conformal quadrangular mesh, by $\mathcal{P}$ the set of points of $\mathcal{M}_h$, $\mathcal{C}$ the set of cells of $\mathcal{M}_h$, by $\mathcal{S}_i$ the set of interior faces of $\mathcal{M}_h$, and by $\mathcal{S}_b$ the set of boundary faces. Boundary faces are supposed to be oriented, such that the normal is outgoing. Each interior side $S$ is arbitrarily oriented, and denoting by $\mathbf{n}^S$ its normal, we will call the left cell the one from which $\mathbf{n}^S$ is outgoing, and the right cell the one in which $\mathbf{n}^S$ is ingoing. The discontinuous Galerkin discretization of the wave system (1) consists in finding $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ such that for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h(\Omega)$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_c v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{v}] \left( \frac{1}{\rho_0} \{ \{ \mathbf{u}_h \cdot \mathbf{n}^S \} \} + d_{11} [p_h] \right) \\ \qquad + \sum_{S \in \mathcal{S}_b} \int_S v^L \left[ \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \right]_b = 0, \\[1.5em] \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{w}] \cdot (\kappa_0 \{ \{ p_h \mathbf{n}^S \} \} + [\mathbf{D}_{22}(\mathbf{n}^S)\mathbf{u}_h]) \\ \qquad + \sum_{S \in \mathcal{S}_b} \int_S \mathbf{w}^L \cdot [\kappa_0 p_h \mathbf{n}^S]_b = 0 \end{array} \right. \quad (24) $$

where the boundary fluxes in (24) are given by

$$ \left[ \begin{array}{c} \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \\ \kappa_0 p_h \mathbf{n}^S \end{array} \right]_{\text{wall}} = \left( \begin{array}{c} 0 \\ \kappa_0 p_h \cdot \mathbf{n}^S + c_0 (\mathbf{u}_h \cdot \mathbf{n}^S) \mathbf{n}^S \end{array} \right) \quad (25) $$

for wall boundary condition and

$$ \left[ \begin{array}{c} \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \\ \kappa_0 p_h \mathbf{n}^S \end{array} \right]_{\text{SW}} = \left( \begin{array}{c} \frac{1}{\rho_0} \frac{\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} \\ \kappa_0 \frac{p_h + p_b}{2} \mathbf{n}^S + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) \mathbf{n}^S \end{array} \right) \quad (26) $$

for inlet/outlet boundary condition.
for wall boundary condition and

$$ \left[ \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \right]_{\text{SW}} = \left( \frac{1}{\rho_0} \frac{\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} \right) \left( \kappa_0 \frac{p_h + p_b}{2} \mathbf{n}^S + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) \mathbf{n}^S \right) \quad (27) $$

for inlet/outlet boundary condition. Note that in the case of inlet/outlet boundary condition, a full state $(p_h, \mathbf{u}_b)$ is enforced. Still, this enforcement is performed weakly through Riemann problems between the interior state and the imposed boundary state. Details on the derivation of expressions (27) and (26) can be found in [29, Equation (2.9)] and [27, Appendix A.].

## 3.2 Discrete Hodge-Helmholtz decomposition

For triangular meshes, a discrete counterpart of the decomposition of Proposition 5 was proven in [27]. It is based on the definition of $\varphi$ as follows

$$ \text{Find } \varphi \in \mathbb{R} \quad \forall g \in \mathbb{R} \quad \sum_{c \in \mathcal{C}} \int_c \nabla_{\mathbf{x}} \varphi \cdot \nabla_{\mathbf{x}} g = \sum_{S \in \mathcal{S}_i} \int_S g[\mathbf{u} \cdot \mathbf{n}^S] + \sum_{S \in \mathcal{S}_b} \int_S g(\mathbf{u} \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S), \quad (28) $$

where $\mathbb{R}$ is the Crouzeix-Raviart finite element space. Then taking $\mathbf{u}_h^\varphi = \nabla_{\mathbf{x}} \varphi$ (where the gradient is computed cellwise) and $\mathbf{u}_h^\psi = \mathbf{u} - \mathbf{u}_h^\varphi$ gives a discrete counterpart of Proposition 5 for triangular meshes.

For Cartesian meshes with periodic boundary conditions, from Proposition 2 and Theorem 2.1, we can write all $\mathbf{u}_h \in d\mathbf{S}_h$, as

$$ \mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi $$

where $\mathbf{u}_h^\varphi \in \nabla_{\mathbf{x}} (\mathbb{R}aTu(\Omega))$ and $\mathbf{u}_h^\psi$ ensures for all $S \in \mathcal{S}_i$, $[\mathbf{u}^\psi \cdot \mathbf{n}^S] = 0$.

We wish to extend this discrete decomposition for general quads and for bounded domain with boundary conditions of type (27) and/or (26), which first raises the problem of the approximation space to choose for the velocity space.

### 3.2.1 The two Piola transformations and the choice of the approximation space

The contravariant Piola transformation takes $\hat{\mathbf{u}}: \hat{K} \mapsto \mathbb{R}^2$ to $\mathbf{u}: c \mapsto \mathbb{R}^2$ defined by

$$ \mathbf{P}_{\mathbf{F}_c}^{\text{div}}(\hat{\mathbf{u}}) := \frac{1}{\det(\mathbf{D}_{\hat{\mathbf{x}}}\mathbf{F}_c)} (\mathbf{D}_{\hat{\mathbf{x}}}\mathbf{F}_c) \hat{\mathbf{u}} \circ \mathbf{F}_c^{-1}. \quad (29) $$

Using this definition gives

$$ \mathrm{curl}_{\mathbf{x}} \psi = \mathbf{P}_{\mathbf{F}_c}^{\mathrm{div}} (\mathrm{curl}_{\hat{\mathbf{x}}} \hat{\psi}) $$

On the other hand, the covariant Piola transformation takes $\hat{\mathbf{u}}: \hat{K} \mapsto \mathbb{R}^2$ to $\mathbf{u}: c \mapsto \mathbb{R}^2$ defined by

$$ \mathbf{P}_{\mathbf{F}_c}^{\mathrm{curl}}(\hat{\mathbf{u}}) := (\mathbf{D}_{\hat{\mathbf{x}}} \mathbf{F}_c)^{-T} \hat{\mathbf{u}} \circ \mathbf{F}_c^{-1}. \quad (30) $$

Using this definition gives

$$ \nabla_{\mathbf{x}} \varphi = \mathbf{P}_{\mathbf{F}_c}^{\mathrm{curl}} (\nabla_{\hat{\mathbf{x}}} \hat{\varphi}). $$

As we are considering the general quadrangular case, the transformations $\mathbf{P}_{\mathbf{F}_c}^{\mathrm{div}}$ and $\mathbf{P}_{\mathbf{F}_c}^{\mathrm{curl}}$ are no more linear inside the elements, and the vector approximation space generated by the parametric finite elements $\mathrm{curl}_{\mathbf{x}} \mathbb{Q}_1$ does not match any more with the parametric $\nabla_{\mathbf{x}} \mathbb{R}aTu$. Based on the proofs of the properties derived in section 2, we consider that the properties ensured by $\mathrm{curl}_{\mathbf{x}} \mathbb{Q}_1$ are important to preserve, which leads us to define the approximation space $d\mathbf{S}_h$ of vector fields on $\Omega$ as

$$ d\mathbf{S}_h(\Omega) = \{ \mathbf{u}_h \in L^2(\Omega)^2 \mid \forall c \in \mathcal{C}, \mathbf{u}_h|_c \in \mathbf{P}_{\mathbf{F}_c}^{\mathrm{div}}(\hat{\mathbf{S}}_0) \}. \quad (31) $$

The contravariant Piola transformation gives for all $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$, for all $p \in d\mathbb{P}_0(\Omega)$, and all cell $c \in \mathcal{C}$,

$$ \int_c p \operatorname{div}_{\mathbf{x}}(\mathbf{u}) = \int_{\hat{K}} \hat{p} \operatorname{div}_{\hat{\mathbf{x}}} (\hat{\mathbf{u}}) \quad \text{and} \quad \forall S \in \mathcal{S}(c), \quad \int_S p \mathbf{u} \cdot \mathbf{n} = \int_{\hat{S}} \hat{p} \hat{\mathbf{u}} \cdot \hat{\mathbf{n}} \quad (32) $$

where $\mathbf{n}$ and $\hat{\mathbf{n}}$ denote the unit outward normals on $S$ and $\hat{S}$.
## 3.2 Discrete Hodge-Helmholtz decomposition

For triangular meshes, a discrete counterpart of the decomposition of Proposition 5 was proven in [25]. It is based on the definition of $\varphi$ as follows

$$ \text{Find } \varphi \in \mathbb{R} \quad \forall g \in \mathbb{R} \quad \sum_{c \in \mathcal{C}} \int_c \nabla_{\mathbf{x}} \varphi \cdot \nabla_{\mathbf{x}} g = \sum_{S \in \mathcal{S}_i} \int_S g[\mathbf{u} \cdot \mathbf{n}^S] + \sum_{S \in \mathcal{S}_b} \int_S g(\mathbf{u} \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S), \qquad (27) $$

where $\mathbb{R}$ is the Crouzeix-Raviart finite element space. Then taking $\mathbf{u}_h^\varphi = \nabla_{\mathbf{x}}\varphi$ (where the gradient is computed cellwise) and $\mathbf{u}_h^\psi = \mathbf{u} - \mathbf{u}_h^\varphi$ gives a discrete counterpart of Proposition 5 for triangular meshes.

For Cartesian meshes with periodic boundary conditions, from Proposition 2 and Theorem 2.1, we can write all $\mathbf{u}_h \in d\mathbf{S}_h$, as

$$ \mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi $$

where $\mathbf{u}_h^\varphi \in \nabla_{\mathbf{x}}(\mathbb{R}aTu(\Omega))$ and $\mathbf{u}_h^\psi$ ensures for all $S \in \mathcal{S}_i$, $[\mathbf{u}^\psi \cdot \mathbf{n}^S] = 0$.

We wish to extend this discrete decomposition for general quads and for bounded domain with boundary conditions of type (26) and/or (25), which first raises the problem of the approximation space to choose for the velocity space.

### 3.2.1 The two Piola transformations and the choice of the approximation space

The contravariant Piola transformation takes $\hat{\mathbf{u}}: \hat{K} \mapsto \mathbb{R}^2$ to $\mathbf{u}: c \mapsto \mathbb{R}^2$ defined by

$$ \mathbf{P}_{\mathbf{F}_c}^{\text{div}}(\hat{\mathbf{u}}) := \frac{1}{\det(\mathbf{D}_{\hat{\mathbf{x}}} \mathbf{F}_c)} (\mathbf{D}_{\hat{\mathbf{x}}} \mathbf{F}_c) \hat{\mathbf{u}} \circ \mathbf{F}_c^{-1}. \qquad (28) $$

Using this definition gives

$$ \nabla_{\mathbf{x}}^{\perp} \psi = \mathbf{P}_{\mathbf{F}_c}^{\text{div}} (\nabla_{\hat{\mathbf{x}}}^{\perp} \hat{\psi}) $$

On the other hand, the covariant Piola transformation takes $\hat{\mathbf{u}}: \hat{K} \mapsto \mathbb{R}^2$ to $\mathbf{u}: c \mapsto \mathbb{R}^2$ defined by

$$ \mathbf{P}_{\mathbf{F}_c}^{\text{curl}}(\hat{\mathbf{u}}) := (\mathbf{D}_{\hat{\mathbf{x}}} \mathbf{F}_c)^{-T} \hat{\mathbf{u}} \circ \mathbf{F}_c^{-1}. \qquad (29) $$

Using this definition gives

$$ \nabla_{\mathbf{x}}\varphi = \mathbf{P}_{\mathbf{F}_c}^{\text{curl}}(\nabla_{\hat{\mathbf{x}}}\hat{\varphi}). $$

As we are considering the general quadrangular case, the transformations $\mathbf{P}_{\mathbf{F}_c}^{\text{div}}$ and $\mathbf{P}_{\mathbf{F}_c}^{\text{curl}}$ are no more linear inside the elements, and the vector approximation space generated by the parametric finite elements $\nabla_{\mathbf{x}}^{\perp} \mathbb{Q}_1$ does not match any more with the parametric $\nabla_{\mathbf{x}} \mathbb{R}aTu$. Based on the proofs of the properties derived in section 2, we consider that the properties ensured by $\nabla_{\mathbf{x}}^{\perp} \mathbb{Q}_1$ are important to preserve, which leads us to define the approximation space $d\mathbf{S}_h$ of vector fields on $\Omega$ as

$$ d\mathbf{S}_h(\Omega) = \{ \mathbf{u}_h \in L^2(\Omega)^2 \mid \forall c \in \mathcal{C}, \mathbf{u}_h|_c \in \mathbf{P}_{\mathbf{F}_c}^{\text{div}}(\hat{\mathbf{S}}_0) \}. \qquad (30) $$

The contravariant Piola transformation gives for all $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$, for all $p \in d\mathbb{P}_0(\Omega)$, and all cell $c \in \mathcal{C}$,

$$ \int_c p \operatorname{div}_{\mathbf{x}}(\mathbf{u}) = \int_{\hat{K}} \hat{p} \operatorname{div}_{\hat{\mathbf{x}}} (\hat{\mathbf{u}}) \quad \text{and} \quad \forall S \in \mathcal{S}(c), \quad \int_S p \mathbf{u} \cdot \mathbf{n} = \int_{\hat{S}} \hat{p} \hat{\mathbf{u}} \cdot \hat{\mathbf{n}} \qquad (31) $$

where $\mathbf{n}$ and $\hat{\mathbf{n}}$ denote the unit outward normals on $S$ and $\hat{S}$.

**Proposition 6** (Finite element for vectors on quadrangular mesh). $d\mathbf{S}_h$ defined by (30) is a finite element space approximating $L^2(\Omega)^2$ at order one with $\dim(d\mathbf{S}_h) = 3\#\mathcal{C}$. Moreover, all $\mathbf{u}_h \in d\mathbf{S}_h$ defined by (30) satisfies for all cell $c \in \mathcal{C}$

$$ \operatorname{div}_{\mathbf{x}} (\mathbf{u}_{h}|_c) = 0 \quad \text{and} \quad \forall S \in \mathcal{S}(c), \quad \mathbf{u}_{h}|_c \cdot \mathbf{n}^S \text{ is constant along } S. $$

Then, the jump $[\mathbf{u}_h \cdot \mathbf{n}^S]$ is also constant along an interior face $S \in \mathcal{S}_i$.

*Proof.* This is a direct consequence of Proposition 1 and (31). □
**Proposition 6** (Finite element for vectors on quadrangular mesh). $d\mathbf{S}_h$ defined by (31) is a finite element space approximating $L^2(\Omega)^2$ at order one with $\dim(d\mathbf{S}_h) = 3\#\mathcal{C}$. Moreover, all $\mathbf{u}_h \in d\mathbf{S}_h$ defined by (31) satisfies for all cell $c \in \mathcal{C}$

$$ \mathrm{div}_{\mathbf{x}} (\mathbf{u}_h|_c) = 0 \quad \text{and} \quad \forall S \in \mathcal{S}(c), \quad \mathbf{u}_h|_c \cdot \mathbf{n}^S \text{ is constant along } S. $$

Then, the jump $\llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket$ is also constant along an interior face $S \in \mathcal{S}_i$.

*Proof.* This is a direct consequence of Proposition 1 and (32). □

### 3.2.2 Determination of $\mathbf{u}_h^\varphi$

The determination of the $\mathbf{u}_h^\varphi$ component is slightly more complicated than in the Cartesian or in the triangular case. This comes from our choice of approximation space (31), which does not contain the gradient of the parametric Rannacher-Turek elements. Note that the parametric Rannacher-Turek finite element has a lot of drawbacks (for example, it is not optimal order on general quads [39]), and the definition of nonconforming finite elements on general quads is still an active research topic [12, 33, 26, 31].

In this article we propose to define directly a subset of the vectorial space (31), based on the properties on the sides that we wish, instead of relying on a variant of the Rannacher-Turek finite element. For this we first remark that the solution of (28) can actually be explicitly computed. Indeed, denoting by $\varphi_S$ the Crouzeix-Raviart basis function associated to the side $S$, such that

$$ \int_S \varphi_S = 1, $$

and denote by $\mathbf{u}^S := \nabla \varphi_S$. Then changing $\nabla \varphi$ by $\mathbf{u}_h^\varphi$, $g$ by $\varphi_S$ and $\nabla g$ by $\mathbf{u}^S$ in (28) gives:

$$ \sum_{c \in \mathcal{C}(S)} \int_c \mathbf{u}_h^\varphi \cdot \mathbf{u}^S = \begin{cases} \int_S \llbracket \mathbf{u} \cdot \mathbf{n}^S \rrbracket & \text{if } S \in \mathcal{S}_i \\ \int_S (\mathbf{u} \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) & \text{if } S \in \mathcal{S}_b. \end{cases} \quad (33) $$

In the quadrangular case, based on (33), for each side $S$, we define $\mathbf{u}^S \in d\mathbf{S}_h$, which support is in $\mathcal{C}(S)$ (the cells adjacent to the side $S$), such that

$$ \begin{align} \forall S \in \mathcal{S}_i \quad & \forall \mathbf{u} \in d\mathbf{S}_h(\Omega) && \int_{\mathcal{C}(S)} \mathbf{u} \cdot \mathbf{u}^S = \int_S \llbracket \mathbf{u} \cdot \mathbf{n}^S \rrbracket, \tag{34} \\ \forall S \in \mathcal{S}_b \quad & \forall \mathbf{u} \in d\mathbf{S}_h(\Omega) && \int_{\mathcal{C}(S)} \mathbf{u} \cdot \mathbf{u}^S = \int_S \mathbf{u} \cdot \mathbf{n}^S. \nonumber \end{align} $$

Such a function can be built by just solving the mass-matrix system based on the equalities of (34) on each cell in $\mathcal{C}(S)$ for each basis member of $d\mathbf{S}_h$ of $\mathcal{C}(S)$. Then equation (34) holds for each $\mathbf{u} \in d\mathbf{S}_h$ by linearity.

**Definition 3.1** (Definition of the space in which $\mathbf{u}_h^\varphi$ will be searched).

$$ d\mathbf{S}_h^\varphi = \mathrm{Span} \{ \mathbf{u}^S, \ S \in \mathcal{S} \}. $$

**Proposition 7** (Dimension of $d\mathbf{S}_h^\varphi$). If a vector $\alpha \in \mathbb{R}^{\#\mathcal{S}}$ is such that

$$ \sum_{S \in \mathcal{S}} \alpha_S \mathbf{u}^S = 0, $$

then

$$ \alpha_1 = \alpha_2 = \dots = \alpha_{\#\mathcal{S}}. $$

This leads to

$$ \dim d\mathbf{S}_h^\varphi = \#\mathcal{S} - 1. $$

*Proof.* Consider a linear combination of $\mathbf{u}^S$ that is zero

$$ \sum_{S \in \mathcal{S}} \alpha_S \mathbf{u}^S = 0. $$

Consider now a cell $c$ and an element $\mathbf{u}_c$ of $d\mathbf{S}_h(\Omega)$, which support is reduced to $c$, such that its trace is 0 along two of the sides, and 1 along one side denoted $S_i$, and -1 along the other side,
### 3.2.2 Determination of $u_h^\varphi$

The determination of the $u_h^\varphi$ component is slightly more complicated than in the Cartesian or in the triangular case. This comes from our choice of approximation space (30), which does not contain the gradient of the parametric Rannacher-Turek elements. Note that the parametric Rannacher-Turek finite element has a lot of drawbacks (for example, it is not optimal order on general quads [35]), and the definition of nonconforming finite elements on general quads is still an active research topic [12, 31, 24, 29].

In this article we propose to define directly a subset of the vectorial space (30), based on the properties on the sides that we wish, instead of relying on a variant of the Rannacher-Turek finite element. For this we first remark that the solution of (27) can actually be explicitly computed. Indeed, denoting by $\varphi_S$ the Crouzeix-Raviart basis function associated to the side $S$, such that

$$ \int_S \varphi_S = 1, $$

and denote by $u^S := \nabla \varphi_S$. Then changing $\nabla \varphi$ by $u_h^\varphi$, $g$ by $\varphi_S$ and $\nabla g$ by $u^S$ in (27) gives:

$$ \sum_{c \in \mathcal{C}(S)} \int_c u_h^\varphi \cdot u^S = \begin{cases} \int_S [\mathbf{u} \cdot \mathbf{n}^S] & \text{if } S \in \mathcal{S}_i \\ \int_S (\mathbf{u} \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) & \text{if } S \in \mathcal{S}_b. \end{cases} \quad (32) $$

In the quadrangular case, based on (32), for each side $S$, we define $u^S \in d\mathbf{S}_h$, which support is in $\mathcal{C}(S)$ (the cells adjacent to the side $S$), such that

$$ \begin{aligned} \forall S \in \mathcal{S}_i \quad & \forall \mathbf{u} \in d\mathbf{S}_h(\Omega) && \int_{\mathcal{C}(S)} \mathbf{u} \cdot \mathbf{u}^S = \int_S [\mathbf{u} \cdot \mathbf{n}^S], \\ \forall S \in \mathcal{S}_b \quad & \forall \mathbf{u} \in d\mathbf{S}_h(\Omega) && \int_{\mathcal{C}(S)} \mathbf{u} \cdot \mathbf{u}^S = \int_S \mathbf{u} \cdot \mathbf{n}^S. \end{aligned} \quad (33) $$

Such a function can be built by just solving the mass-matrix system based on the equalities of (33) on each cell in $\mathcal{C}(S)$ for each basis member of $d\mathbf{S}_h$ of $\mathcal{C}(S)$. Then equation (33) holds for each $\mathbf{u} \in d\mathbf{S}_h$ by linearity.

**Definition 3.1** (Definition of the space in which $u_h^\varphi$ will be searched).

$$ d\mathbf{S}_h^\varphi = \operatorname{Span}\{\mathbf{u}^S, \ S \in \mathcal{S}\}. $$

**Proposition 7** (Dimension of $d\mathbf{S}_h^\varphi$). If a vector $\alpha \in \mathbb{R}^{\#\mathcal{S}}$ is such that

$$ \sum_{S \in \mathcal{S}} \alpha_S \mathbf{u}^S = 0, $$

then

$$ \alpha_1 = \alpha_2 = \dots = \alpha_{\#\mathcal{S}}. $$

This leads to

$$ \dim d\mathbf{S}_h^\varphi = \#\mathcal{S} - 1. $$

*Proof.* Consider a linear combination of $u^S$ that is zero

$$ \sum_{S \in \mathcal{S}} \alpha_S \mathbf{u}^S = 0. $$

Consider now a cell $c$ and an element $\mathbf{u}_c$ of $d\mathbf{S}_h(\Omega)$, which support is reduced to $c$, such that its trace is 0 along two of the sides, and 1 along one side denoted $S_i$, and $-1$ along the other side, denoted $S_j$. Such an element $\mathbf{u}_c$ exists because of Proposition 6: the three degrees of freedom in the cell $c$ allow to define a $\mathbf{u}_c$ such that it is 0 on two sides and 1 on a third side, and as $\mathbf{u}_c$ is divergence free, its trace is $-1$ on the last side. Then

$$ \int_{\Omega} \mathbf{u}_c \cdot \left( \sum_{S \in \mathcal{S}} \alpha_S \mathbf{u}^S \right) = \sum_{S \in \mathcal{S}} \alpha_S \int_{\Omega} \mathbf{u}_c \cdot \mathbf{u}^S = \alpha_{S_i} \int_{S_i} [\mathbf{u}_c \cdot \mathbf{n}^{S_i}] + \alpha_{S_j} \int_{S_j} [\mathbf{u}_c \cdot \mathbf{n}^{S_j}] = \alpha_{S_i} - \alpha_{S_j}, $$

which leads to $\alpha_{S_i} = \alpha_{S_j}$, and gives $\dim d\mathbf{S}_h^\varphi \le \#\mathcal{S} - 1$. We now consider the vector

$$ \mathbf{s} := \sum_{S \in \mathcal{S}} \mathbf{u}^S. $$
denoted $S_j$. Such an element $\mathbf{u}_c$ exists because of Proposition 6: the three degrees of freedom in the cell $c$ allow to define a $\mathbf{u}_c$ such that it is 0 on two sides and 1 on a third side, and as $\mathbf{u}_c$ is divergence free, its trace is $-1$ on the last side. Then

$$ \int_{\Omega} \mathbf{u}_c \cdot \left( \sum_{S \in \mathcal{S}} \alpha_S \mathbf{u}^S \right) = \sum_{S \in \mathcal{S}} \alpha_S \int_{\Omega} \mathbf{u}_c \cdot \mathbf{u}^S = \alpha_{S_i} \int_{S_i} \left[ \left[ \mathbf{u}_c \cdot \mathbf{n}^{S_i} \right] \right] + \alpha_{S_j} \int_{S_j} \left[ \left[ \mathbf{u}_c \cdot \mathbf{n}^{S_j} \right] \right] = \alpha_{S_i} - \alpha_{S_j}, $$

which leads to $\alpha_{S_i} = \alpha_{S_j}$, and gives dim $d\mathbf{S}_h^\varphi \le \#\mathcal{S} - 1$. We now consider the vector

$$ \mathbf{s} := \sum_{S \in \mathcal{S}} \mathbf{u}^S. $$

Then for all $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$

$$ \begin{aligned} \int_{\Omega} \mathbf{u}_h \cdot \mathbf{s} &= \sum_{S \in \mathcal{S}_i} \int_S \left[ \left[ \mathbf{u}_h \cdot \mathbf{n}^S \right] \right] + \sum_{S \in \mathcal{S}_b} \int_S \mathbf{u}_h \cdot \mathbf{n}^S \\ &= \sum_{c \in \mathcal{C}} \int_c \operatorname{div}_{\mathbf{x}} \mathbf{u}_h \\ &= 0, \end{aligned} $$

because all the elements of $d\mathbf{S}_h(\Omega)$ have a zero divergence on each cell. This ends the proof. $\square$

Now, based on the variational formulation (28), we propose to define $\mathbf{u}_h^\varphi$ as follows

**Proposition 8 (Definition of $\mathbf{u}_h^\varphi$).** We denote by $\mathbf{u}_h$ and element of $d\mathbf{S}_h(\Omega)$. Suppose that $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$. Then the problem to find $\beta \in \mathbb{R}^{\#\mathcal{S}}/(1, 1...1)$ such that $\mathbf{v}_h = \sum_S \beta_S \mathbf{u}^S$, ensuring

$$ \forall \alpha \in \mathbb{R}^{\#\mathcal{S}}/(1, 1...1) \qquad \sum_S \alpha_S \int_{\Omega} \mathbf{v} \cdot \mathbf{u}^S = \sum_{S \in \mathcal{S}_i} \alpha_S \int_S \left[ \left[ \mathbf{u}_h \cdot \mathbf{n}^S \right] \right] + \sum_{S \in \mathcal{S}_b} \alpha_S \int_S (\mathbf{u}_h - \mathbf{u}_b) \cdot \mathbf{n}^S, \qquad (35) $$

has a unique solution $\beta^\varphi$, and we denote by

$$ \mathbf{u}_h^\varphi := \sum_S \beta_S^\varphi \mathbf{u}^S. $$

*Proof.* The proof relies on the Lax-Milgram theorem. The right hand side is clearly a linear form. The left hand side is clearly a positive bilinear form. If the problem (35) is considered for $\alpha \in \mathbb{R}^{\#\mathcal{S}}$, then the problem is singular because of Proposition 7. Note however that as $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$, the problem (35) considered for $\alpha \in \mathbb{R}^{\#\mathcal{S}}$ has an infinite number of solutions. Last, the problem (35) is coercive because we are working on $\mathbb{R}^{\#\mathcal{S}}/(1, 1...1)$, on which $\mathbf{u}_h^\varphi = 0$ if and only if $\beta = 0$. $\square$

### 3.2.3 Existence and uniqueness of the decomposition

We are now able to define a Hodge-Helmholtz decomposition on general quadrangular mesh with boundary condition.

**Proposition 9 (Discrete Hodge-Helmholtz decomposition).** Let $\Omega$ a connected bounded domain with possibly $r \in \mathbb{N}$ hole(s) and boundary conditions of type (26) or/and (27). Suppose that a boundary value $\mathbf{u}_b \cdot \mathbf{n}$ is known on the boundary such that $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$. Then a velocity $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$ on a general quadrangular mesh $\mathcal{M}_h$ can be uniquely written as

$$ \mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi $$

where

$$ \bullet \quad \mathbf{u}_h^\varphi \in d\mathbf{S}_h^\varphi $$

$$ \bullet \quad \mathbf{u}_h^\psi \in d\mathbf{S}_h^\psi \text{ such that} $$

$$ d\mathbf{S}_h^\psi := \{ \mathbf{v}_h \in d\mathbf{S}_h \mid \forall S \in \mathcal{S}_i, \left[ \left[ \mathbf{v}_h \cdot \mathbf{n}^S \right] \right] = 0 \quad \text{and} \quad \forall S \in \mathcal{S}_b, \mathbf{v}_h \cdot \mathbf{n}^S = \mathbf{u}_b \cdot \mathbf{n}^S \}. $$
Then for all $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$

$$ 
\begin{aligned}
\int_{\Omega} \mathbf{u}_h \cdot \mathbf{s} &= \sum_{S \in \mathcal{S}_i} \int_S [[\mathbf{u}_h \cdot \mathbf{n}^S]] + \sum_{S \in \mathcal{S}_b} \int_S \mathbf{u}_h \cdot \mathbf{n}^S \\
&= \sum_{c \in \mathcal{C}} \int_c \mathrm{div}_{\mathbf{x}} \mathbf{u}_h \\
&= 0,
\end{aligned}
$$

because all the elements of $d\mathbf{S}_h(\Omega)$ have a zero divergence on each cell. This ends the proof. $\square$

Now, based on the variational formulation (27), we propose to define $\mathbf{u}_h^\varphi$ as follows

**Proposition 8 (Definition of $\mathbf{u}_h^\varphi$).** We denote by $\mathbf{u}_h$ and element of $d\mathbf{S}_h(\Omega)$. Suppose that $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$. Then the problem to find $\beta \in \mathbb{R}^{\#\mathcal{S}}/(1, 1...1)$ such that $\mathbf{v}_h = \sum_S \beta_S \mathbf{u}^S$, ensuring

$$ \forall \alpha \in \mathbb{R}^{\#\mathcal{S}}/(1, 1...1) \qquad \sum_S \alpha_S \int_{\Omega} \mathbf{v} \cdot \mathbf{u}^S = \sum_{S \in \mathcal{S}_i} \alpha_S \int_S [[\mathbf{u}_h \cdot \mathbf{n}^S]] + \sum_{S \in \mathcal{S}_b} \alpha_S \int_S (\mathbf{u}_h - \mathbf{u}_b) \cdot \mathbf{n}^S, \qquad (34) $$

has a unique solution $\beta^\varphi$, and we denote by

$$ \mathbf{u}_h^\varphi := \sum_S \beta_S^\varphi \mathbf{u}^S. $$

*Proof.* The proof relies on the Lax-Milgram theorem. The right hand side is clearly a linear form. The left hand side is clearly a positive bilinear form. If the problem (34) is considered for $\alpha \in \mathbb{R}^{\#\mathcal{S}}$, then the problem is singular because of Proposition 7. Note however that as $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$, the problem (34) considered for $\alpha \in \mathbb{R}^{\#\mathcal{S}}$ has an infinite number of solutions. Last, the problem (34) is coercive because we are working on $\mathbb{R}^{\#\mathcal{S}}/(1, 1...1)$, on which $\mathbf{u}_h^\varphi = 0$ if and only if $\beta = 0$. $\square$

### 3.2.3 Existence and uniqueness of the decomposition

We are now able to define a Hodge-Helmholtz decomposition on general quadrangular mesh with boundary condition.

**Proposition 9 (Discrete Hodge-Helmholtz decomposition).** Let $\Omega$ a connected bounded domain with possibly $r \in \mathbb{N}$ hole(s) and boundary conditions of type (25) or/and (26). Suppose that a boundary value $\mathbf{u}_b \cdot \mathbf{n}$ is known on the boundary such that $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$. Then a velocity $\mathbf{u}_h \in d\mathbf{S}_h(\Omega)$ on a general quadrangular mesh $\mathcal{M}_h$ can be uniquely written as

$$ \mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi $$

where

$$ \bullet \quad \mathbf{u}_h^\varphi \in d\mathbf{S}_h^\varphi $$

$$ \bullet \quad \mathbf{u}_h^\psi \in d\mathbf{S}_h^\psi \text{ such that} $$

$$ d\mathbf{S}_h^\psi := \{ \mathbf{v}_h \in d\mathbf{S}_h \mid \forall S \in \mathcal{S}_i, [[\mathbf{v}_h \cdot \mathbf{n}^S]] = 0 \quad \text{and} \quad \forall S \in \mathcal{S}_b, \mathbf{v}_h \cdot \mathbf{n}^S = \mathbf{u}_b \cdot \mathbf{n}^S \}. $$

Moreover

$$ \dim(d\mathbf{S}_h^\psi) = 3\#\mathcal{C} - \#\mathcal{S} + 1. $$

We note that unlike the periodic case of Proposition 2, due to the boundary conditions, the decomposition of Proposition 9 is no more orthogonal.

*Proof.* We first address the existence of the decomposition. Given $\mathbf{u}_h$, we define $\mathbf{u}_h^\varphi$ as in Proposition 8, and define $\mathbf{u}_h^\psi$ as

$$ \mathbf{u}_h^\psi := \mathbf{u}_h - \mathbf{u}_h^\varphi. $$

It remains to check that $\mathbf{u}_h^\psi \in d\mathbf{S}_h^\psi$. For this, the formula (34) is tested with $\alpha_S = 0$ for all $S$, except for one, which we denote as $S_k$, for which $\alpha_{S_k} = 1$. Then
Moreover

$$ \dim(d\mathbf{S}_h^\psi) = 3\#\mathcal{C} - \#\mathcal{S} + 1. $$

We note that unlike the periodic case of Proposition 2, due to the boundary conditions, the decomposition of Proposition 9 is no more orthogonal.

**Proof.** We first address the existence of the decomposition. Given $\mathbf{u}_h$, we define $\mathbf{u}_h^\varphi$ as in Proposition 8, and define $\mathbf{u}_h^\psi$ as

$$ \mathbf{u}_h^\psi := \mathbf{u}_h - \mathbf{u}_h^\varphi. $$

It remains to check that $\mathbf{u}_h^\psi \in d\mathbf{S}_h^\psi$. For this, the formula (35) is tested with $\alpha_S = 0$ for all $S$, except for one, which we denote as $S_k$, for which $\alpha_{S_k} = 1$. Then

* If $S_k \in \mathcal{S}_i$ then the left hand side is

$$ \int_{\Omega} \mathbf{u}_h^\varphi \cdot \mathbf{u}_{S_k} = \int_{S_k} [\mathbf{u}_h^\varphi \cdot \mathbf{n}^{S_k}], $$

where (34) was used. The right hand side is

$$ \int_{S_k} [\mathbf{u}_h \cdot \mathbf{n}^{S_k}]. $$

This leads to

$$ \int_{S_k} [\mathbf{u}_h^\psi \cdot \mathbf{n}^{S_k}] = \int_{S_k} [(\mathbf{u}_h - \mathbf{u}_h^\varphi) \cdot \mathbf{n}^{S_k}] = 0, $$

and as all the elements of $d\mathbf{S}_h$ have a constant trace on each side, this gives $\mathbf{u}_h^\psi \cdot \mathbf{n}^{S_k} = 0$.

* If $S_k \in \mathcal{S}_b$ then the left hand side is still

$$ \int_{\Omega} \mathbf{u}_h^\varphi \cdot \mathbf{u}^{S_k} = \int_{S_k} [\mathbf{u}_h^\varphi \cdot \mathbf{n}^{S_k}], $$

whereas the right hand side is

$$ \int_{S_k} [(\mathbf{u}_h - \mathbf{u}_b) \cdot \mathbf{n}^{S_k}]. $$

This leads to

$$ \int_{S_k} [(\mathbf{u}_h^\psi - \mathbf{u}_b) \cdot \mathbf{n}^{S_k}] = \int_{S_k} [(\mathbf{u}_h - \mathbf{u}_h^\varphi) \cdot \mathbf{n}^{S_k}] = 0, $$

and so $\mathbf{u}_h^\psi \cdot \mathbf{n}^{S_k} = \mathbf{u}_b \cdot \mathbf{n}^{S_k}$.

Therefore, $\mathbf{u}_h^\psi \in d\mathbf{S}_h^\psi$, and this proves the existence of the decomposition.

We now would like to prove uniqueness of the decomposition. Suppose that a given $\mathbf{u}_h$ can be decomposed as

$$ \mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi = \mathbf{v}_h^\varphi + \mathbf{v}_h^\psi. $$

Then we define $\mathbf{d}_h$ such that

$$ \mathbf{d}_h := \mathbf{u}_h^\varphi - \mathbf{v}_h^\varphi = \mathbf{v}_h^\psi - \mathbf{u}_h^\psi, $$

As $[\mathbf{d}_h] = 0$ on all interior and boundary sides (because $\mathbf{v}_h^\psi, \mathbf{u}_h^\psi \in d\mathbf{S}_h^\psi$), this gives

$$ \sum_{S} \alpha_S \int_{\Omega} \mathbf{d}_h \cdot \mathbf{u}^S = 0. $$

As $\mathbf{d}_h = \mathbf{u}_h^\varphi - \mathbf{v}_h^\varphi \in d\mathbf{S}_h^\varphi$, and by uniqueness of the solution of (35), we have $\mathbf{d}_h = 0$, which proves uniqueness.

As $\dim d\mathbf{S}_h = 3\#\mathcal{C}$ and $\dim d\mathbf{S}_h^\varphi = \#\mathcal{S} - 1$, the existence and uniqueness of the decomposition gives

$$ \dim d\mathbf{S}_h^\psi = \dim d\mathbf{S}_h - \dim d\mathbf{S}_h^\varphi = 3\#\mathcal{C} - \#\mathcal{S} + 1, $$

which ends the proof.

□

**Proposition 10** (Characterization of $d\mathbf{S}_h^\psi$). We denote by $\bar{\mathbf{u}}^\varphi$ and $\bar{\mathbf{u}}^\psi$ the decomposition of the null vector given by Proposition 9. Then
• If $S_k \in \mathcal{S}_i$ then the left hand side is

$$ \int_{\Omega} \mathbf{u}_h^{\varphi} \cdot \mathbf{u}_{S_k} = \int_{S_k} [\mathbf{u}_h^{\varphi} \cdot \mathbf{n}^{S_k}], $$

where (33) was used. The right hand side is

$$ \int_{S_k} [\mathbf{u}_h \cdot \mathbf{n}^{S_k}]. $$

This leads to

$$ \int_{S_k} [\mathbf{u}_h^{\psi} \cdot \mathbf{n}^{S_k}] = \int_{S_k} [(\mathbf{u}_h - \mathbf{u}_h^{\varphi}) \cdot \mathbf{n}^{S_k}] = 0, $$

and as all the elements of $d\mathbf{S}_h$ have a constant trace on each side, this gives $\mathbf{u}_h^{\psi} \cdot \mathbf{n}^{S_k} = 0$.

• If $S_k \in \mathcal{S}_b$ then the left hand side is still

$$ \int_{\Omega} \mathbf{u}_h^{\varphi} \cdot \mathbf{u}^{S_k} = \int_{S_k} [\mathbf{u}_h^{\varphi} \cdot \mathbf{n}^{S_k}], $$

whereas the right hand side is

$$ \int_{S_k} [(\mathbf{u}_h - \mathbf{u}_b) \cdot \mathbf{n}^{S_k}]. $$

This leads to

$$ \int_{S_k} [ (\mathbf{u}_h^{\psi} - \mathbf{u}_b) \cdot \mathbf{n}^{S_k} ] = \int_{S_k} [ (\mathbf{u}_h - \mathbf{u}_h^{\varphi}) \cdot \mathbf{n}^{S_k} ] = 0, $$

and so $\mathbf{u}_h^{\psi} \cdot \mathbf{n}^{S_k} = \mathbf{u}_b \cdot \mathbf{n}^{S_k}$.

Therefore, $\mathbf{u}_h^{\psi} \in d\mathbf{S}_h^{\psi}$, and this proves the existence of the decomposition.

We now would like to prove uniqueness of the decomposition. Suppose that a given $\mathbf{u}_h$ can be decomposed as

$$ \mathbf{u}_h = \mathbf{u}_h^{\varphi} + \mathbf{u}_h^{\psi} = \mathbf{v}_h^{\varphi} + \mathbf{v}_h^{\psi}. $$

Then we define $\mathbf{d}_h$ such that

$$ \mathbf{d}_h := \mathbf{u}_h^{\varphi} - \mathbf{v}_h^{\varphi} = \mathbf{v}_h^{\psi} - \mathbf{u}_h^{\psi}, $$

As $[\mathbf{d}_h] = 0$ on all interior and boundary sides (because $\mathbf{v}_h^{\psi}, \mathbf{u}_h^{\psi} \in d\mathbf{S}_h^{\psi}$), this gives

$$ \sum_{S} \alpha_S \int_{\Omega} \mathbf{d}_h \cdot \mathbf{u}^S = 0. $$

As $\mathbf{d}_h = \mathbf{u}_h^{\varphi} - \mathbf{v}_h^{\varphi} \in d\mathbf{S}_h^{\varphi}$, and by uniqueness of the solution of (34), we have $\mathbf{d}_h = 0$, which proves uniqueness.

As dim $d\mathbf{S}_h = 3\#\mathcal{C}$ and dim $d\mathbf{S}_h^{\varphi} = \#\mathcal{S} - 1$, the existence and uniqueness of the decomposition gives

$$ \dim d\mathbf{S}_h^{\psi} = \dim d\mathbf{S}_h - \dim d\mathbf{S}_h^{\varphi} = 3\#\mathcal{C} - \#\mathcal{S} + 1, $$

which ends the proof.

□

**Proposition 10** (Characterization of $d\mathbf{S}_h^{\psi}$). We denote by $\bar{\mathbf{u}}^{\varphi}$ and $\bar{\mathbf{u}}^{\psi}$ the decomposition of the null vector given by Proposition 9. Then

• if $\psi_p$ denotes an element of the canonical $\mathbb{Q}_1$ basis function for an interior point $p \in \mathcal{P}_i$, then $\nabla^{\perp} \psi_p + \bar{\mathbf{u}}^{\psi} \in d\mathbf{S}_h^{\psi}$

• if $\psi_{\partial_k}$ denotes the element of $\mathbb{Q}_1$ that is 1 on the whole kth connected component of $\partial\Omega$ and 0 elsewhere, then $\nabla^{\perp} \psi_{\partial_k} + \bar{\mathbf{u}}^{\psi} \in d\mathbf{S}_h^{\psi}$.

Reciprocally, any element of $d\mathbf{S}_h^{\psi}$ can be written as a linear combination of $\bar{\mathbf{u}}^{\psi}, \nabla^{\perp} \psi_p$ for $p \in \mathcal{P}_i$ and $\nabla^{\perp} \psi_{\partial_k}$.
* if $\psi_p$ denotes an element of the canonical $\mathbb{Q}_1$ basis function for an interior point $p \in \mathcal{P}_i$, then $\mathrm{curl}_x \psi_p + \bar{\mathbf{u}}^\psi \in d\mathbf{S}_h^\psi$

* if $\psi_{\partial_k}$ denotes the element of $\mathbb{Q}_1$ that is 1 on the whole kth connected component of $\partial\Omega$ and 0 elsewhere, then $\mathrm{curl}_x \psi_{\partial_k} + \bar{\mathbf{u}}^\psi \in d\mathbf{S}_h^\psi$.

Reciprocally, any element of $d\mathbf{S}_h^\psi$ can be written as a linear combination of $\bar{\mathbf{u}}^\psi$, $\mathrm{curl}_x \psi_p$ for $p \in \mathcal{P}_i$ and $\mathrm{curl}_x \psi_{\partial_k}$.

*Proof.* All the $\mathrm{curl}_x \psi_p$ and $\mathrm{curl}_x \psi_{\partial_k}$ are such that the jump is 0 along all boundary and interior sides. Then it is clear that by adding $\bar{\mathbf{u}}^\psi$, they all belong to $d\mathbf{S}_h^\psi$. Also, we remark that all these functions are not linearly independent, as the function of $\mathbb{Q}_1$ equal to 1 on all points belong to this space. This gives

$$ \dim \operatorname{span}(\mathrm{curl}_x \psi_p, \mathrm{curl}_x \psi_{\partial_k}) = \#\mathcal{P}_i + r. $$

It remains to prove that

$$ \dim d\mathbf{S}_h^\psi = \#\mathcal{P}_i + r. $$

For this, we first use the Euler relation in the mesh

$$ \#\mathcal{P} - \#\mathcal{S} + \#\mathcal{C} = 1 - r. \qquad (36) $$

We then remark that the sum

$$ \sum_{c \in \mathcal{C}} \sum_{S \in \mathcal{S}(c)} 1, $$

can be computed in two ways, and gives $4\#\mathcal{C}$ on one hand, and $2\#\mathcal{P}_i + \#\mathcal{P}_b$ on the other hand:

$$ 4\#\mathcal{C} = 2\#\mathcal{P}_i + \#\mathcal{P}_b. $$

Using this last formula for eliminating $\#\mathcal{P} = \#\mathcal{P}_i + \#\mathcal{P}_b$ from (36) leads to

$$ \#\mathcal{P} - (4\#\mathcal{C} - \#\mathcal{P}_i) + \#\mathcal{C} = 1 - r, $$

which provides the equality

$$ 3\#\mathcal{C} = \#\mathcal{P} + \#\mathcal{P}_i + r - 1. $$

Last, remarking that $\#\mathcal{P}_b = \#\mathcal{P}_b$, we get

$$ 3\#\mathcal{C} - \#\mathcal{P} + 1 = \#\mathcal{P}_i + r, $$

which ends the proof.

□

**Proposition 11** (Characterization of $(d\mathbf{S}_h^\varphi)^\perp$). We have

$$ \mathrm{span}(\mathrm{curl}_x \psi_p, \mathrm{curl}_x \psi_{\partial_k}) = (d\mathbf{S}_h^\varphi)^\perp $$

*Proof.* If $p \in \mathcal{P}_i$, we have for all $\alpha \in \mathbb{R}^{\#\mathcal{S}}$,

$$ \begin{aligned} \int_{\Omega} \mathrm{curl}_x \psi_p \cdot \left( \sum_S \alpha_S \mathbf{u}^S \right) &= \sum_{c \in \mathcal{C}} \int_c \mathrm{curl}_x \psi_p \cdot \left( \sum_S \alpha_S \mathbf{u}^S \right) \\ &= \sum_S \alpha_S \sum_{c \in \mathcal{C}(S)} \int_c \mathrm{curl}_x \psi_p \cdot \mathbf{u}^S \\ &= \sum_{S \in \mathcal{S}_i} \alpha_S \int_S \left[ \mathrm{curl}_x \psi_p \cdot \mathbf{n}^S \right] + \sum_{S \in \mathcal{S}_b} \alpha_S \int_S \mathrm{curl}_x \psi_p \cdot \mathbf{n}^S \\ &= 0. \end{aligned} $$

Similarly, if $\partial_k$ is a connected component of $\partial\Omega$, we have for all $\alpha \in \mathbb{R}^{\#\mathcal{S}}$,

$$ \int_{\Omega} \mathrm{curl}_x \psi_{\partial_k} \cdot \left( \sum_S \alpha_S \mathbf{u}^S \right) = 0, $$

so that $\mathrm{span}(\mathrm{curl}_x \psi_p, \mathrm{curl}_x \psi_{\partial_k}) \subset (d\mathbf{S}_h^\varphi)^\perp$. Since

$$ \dim \operatorname{span}(\mathrm{curl}_x \psi_p, \mathrm{curl}_x \psi_{\partial_k}) = \#\mathcal{P}_i + r = 3\#\mathcal{C} - \#\mathcal{S} + 1 = \dim (d\mathbf{S}_h^\varphi)^\perp, $$

it ends the proof.

□
Proof. All the $\nabla^{\perp}\psi_p$ and $\nabla^{\perp}\psi_{\partial_k}$ are such that the jump is 0 along all boundary and interior sides. Then it is clear that by adding $\bar{u}^{\psi}$, they all belong to $d\mathbf{S}_h^{\psi}$. Also, we remark that all these functions are not linearly independent, as the function of $Q_1$ equal to 1 on all points belong to this space. This gives

$$ \dim \operatorname{span} (\nabla^{\perp} \psi_p, \nabla^{\perp} \psi_{\partial_k}) = \# \mathcal{P}_i + r. $$

It remains to prove that

$$ \dim d\mathbf{S}_h^{\psi} = \# \mathcal{P}_i + r. $$

For this, we first use the Euler relation in the mesh

$$ \# \mathcal{P} - \# \mathcal{S} + \# \mathcal{C} = 1 - r. \qquad (35) $$

We then remark that the sum

$$ \sum_{c \in \mathcal{C}} \sum_{S \in \mathcal{S}(c)} 1, $$

can be computed in two ways, and gives $4\#\mathcal{C}$ on one hand, and $2\#\mathcal{P}_i + \#\mathcal{P}_b$ on the other hand:

$$ 4\#\mathcal{C} = 2\#\mathcal{P}_i + \#\mathcal{P}_b. $$

Using this last formula for eliminating $\#\mathcal{P} = \#\mathcal{P}_i + \#\mathcal{P}_b$ from (35) leads to

$$ \# \mathcal{P} - (4\#\mathcal{C} - \#\mathcal{P}_i) + \#\mathcal{C} = 1 - r, $$

which provides the equality

$$ 3\#\mathcal{C} = \#\mathcal{P} + \#\mathcal{P}_i + r - 1. $$

Last, remarking that $\#\mathcal{P}_b = \#\mathcal{P}_i$, we get

$$ 3\#\mathcal{C} - \#\mathcal{P} + 1 = \#\mathcal{P}_i + r, $$

which ends the proof.

□

**Proposition 11** (Characterization of $(d\mathbf{S}_h^{\varphi})^{\perp}$). We have

$$ \operatorname{span} (\nabla^{\perp} \psi_p, \nabla^{\perp} \psi_{\partial_k}) = (d\mathbf{S}_h^{\varphi})^{\perp} $$

Proof. If $p \in \mathcal{P}_i$, we have for all $\alpha \in \mathbb{R}^{\#\mathcal{S}}$,

$$ \begin{aligned} \int_{\Omega} \nabla^{\perp} \psi_p \cdot \left( \sum_S \alpha_S \mathbf{u}^S \right) &= \sum_{c \in \mathcal{C}} \int_c \nabla^{\perp} \psi_p \cdot \left( \sum_S \alpha_S \mathbf{u}^S \right) \\ &= \sum_S \alpha_S \sum_{c \in \mathcal{C}(S)} \int_c \nabla^{\perp} \psi_p \cdot \mathbf{u}^S \\ &= \sum_{S \in \mathcal{S}_i} \alpha_S \int_S \left[ \left\| \nabla^{\perp} \psi_p \cdot \mathbf{n}^S \right\| \right] + \sum_{S \in \mathcal{S}_b} \alpha_S \int_S \nabla^{\perp} \psi_p \cdot \mathbf{n}^S \\ &= 0. \end{aligned} $$

Similarly, if $\partial_k$ is a connected component of $\partial\Omega$, we have for all $\alpha \in \mathbb{R}^{\#\mathcal{S}}$,

$$ \int_{\Omega} \nabla^{\perp} \psi_{\partial_k} \cdot \left( \sum_S \alpha_S \mathbf{u}^S \right) = 0, $$

so that $\operatorname{span} (\nabla^{\perp} \psi_p, \nabla^{\perp} \psi_{\partial_k}) \subset (d\mathbf{S}_h^{\varphi})^{\perp}$. Since

$$ \dim \operatorname{span} (\nabla^{\perp} \psi_p, \nabla^{\perp} \psi_{\partial_k}) = \# \mathcal{P}_i + r = 3\#\mathcal{C} - \#\mathcal{P} + 1 = \dim (d\mathbf{S}_h^{\varphi})^{\perp}, $$

it ends the proof.

□
### 3.3 Structure preserved and long time behaviour with Godunov' flux

Discretization (25) with Godunov stabilization can be written as finding $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ such that for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h$, we have

$$ 
\left\{
\begin{aligned}
& \sum_{c \in \mathcal{C}} \int_C v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{v}] \left( \frac{1}{\rho_0} \{\{ \mathbf{u}_h \cdot \mathbf{n}^S \}\} + \frac{c_0}{2} [\mathbf{p}_h] \right) \\
& \qquad + \sum_{S \in \mathcal{S}_{SW}} \int_S v \left( \frac{1}{\rho_0} \frac{\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} + \frac{c_0}{2} (p_h - p_b) \right) = 0, \quad (37a) \\
\\
& \sum_{c \in \mathcal{C}} \int_C \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{w} \cdot \mathbf{n}^S] \left( \kappa_0 \{\{ p_h \}\} + \frac{c_0}{2} [\mathbf{u}_h \cdot \mathbf{n}^S] \right) \\
& \qquad + \sum_{S \in \mathcal{S}_{Wall}} \int_S \mathbf{w} \cdot \mathbf{n}^S (\kappa_0 p_h + c_0 \mathbf{u}_h \cdot \mathbf{n}^S) \\
& \qquad + \sum_{S \in \mathcal{S}_{SW}} \int_S \mathbf{w} \cdot \mathbf{n}^S \left( \kappa_0 \frac{p_h + p_b}{2} + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) \right) = 0. \quad (37b)
\end{aligned}
\right.
$$

The following propositions are the discrete version of Proposition 5.

**Proposition 12** (Discrete Hodge-Helmholtz decomposition adapted to the wave system). Denote $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ the solution of the wave system discretization (25) with Godunov stabilization (18). The Hodge-Helmholtz decomposition $\mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi$ of Proposition 9 is adapted to the wave system in the sense that $\mathbf{u}_h^\psi$ is constant in time.

*Proof.* For all $p \in \mathcal{P}_i$, since $[\mathbf{curl}_x \psi_p \cdot \mathbf{n}^S] = 0$ for all $S \in \mathcal{S}_i$ and $\mathbf{curl}_x \psi_p \cdot \mathbf{n}^S = 0$ for all $S \in \mathcal{S}_b$, (37b) with forward or backward Euler method for the time integration gives

$$ \forall p \in \mathcal{P}_i \quad \sum_{c \in \mathcal{C}} \int_C \mathbf{curl}_x \psi_p \cdot \partial_t \mathbf{u}_h = 0 $$

so that

$$ \forall p \in \mathcal{P}_i \quad \sum_{c \in \mathcal{C}} \int_C \mathbf{curl}_x \psi_p \cdot (\partial_t \mathbf{u}_h^\varphi + \partial_t \mathbf{u}_h^\psi) = 0. $$

Using Proposition 11, it gives

$$ \forall p \in \mathcal{P}_i \quad \sum_{c \in \mathcal{C}} \int_C \mathbf{curl}_x \psi_p \cdot (\partial_t \mathbf{u}_h^\psi) = 0. \quad (38) $$

In the same way, we obtain for all $\partial_k$ connected component of $\partial\Omega$ that

$$ \sum_{c \in \mathcal{C}} \int_C \mathbf{curl}_x \psi_{\partial_k} \cdot (\partial_t \mathbf{u}_h^\psi) = 0. \quad (39) $$

Then, using the characterization of $d\mathbf{S}_h^\psi$ of Proposition 10, we obtain with (38) and (39) that

$$ \partial_t \mathbf{u}_h^\psi = 0 $$

which ends the proof.

□

**Proposition 13** (Long time limit). We consider the discretization (25) of the wave system (1) with Godunov stabilization (18). The initial conditions are denoted by $(p_h^0, \mathbf{u}_h^0)$. Suppose that the boundary conditions ensure

* $p = p_b$ where $p_b$ is uniform.

* $\mathbf{u} = \mathbf{u}_b$ where $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$.

Then the long time limit exists and is denoted by $(p_h^\infty, \mathbf{u}_h^\infty)$, it satisfies $p_h^\infty$ is uniform (equals to $p_b$) and $\mathbf{u}_h^\infty = \mathbf{u}_h^\psi(0)$ where $\mathbf{u}_h^\psi(0)$ corresponds to the divergence free component of the initial velocity field $\mathbf{u}_h^0$ obtained with Proposition 9. Moreover, defining the relative energy as

$$ E_{\mathbf{u}_h^\psi(0)}(\mathbf{U}_h) := \frac{\rho_0^2 c_0^2}{2} (p_h - p_b)^2 + \frac{1}{2} \| \mathbf{u}_h - \mathbf{u}_h^\psi(0) \|^2, $$

the volumic average of the relative energy is a Liapunov functional.
### 3.3 Structure preserved and long time behaviour with Godunov' flux

Discretization (24) with Godunov stabilization can be written as finding $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ such that for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h$, we have

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_C v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{v}] \left( \frac{1}{\rho_0} \{\{ \mathbf{u}_h \cdot \mathbf{n}^S \}\} + \frac{c_0}{2} [\mathbf{p}_h] \right) \\ \qquad + \sum_{S \in \mathcal{S}_{SW}} \int_S v \left( \frac{1}{\rho_0} \frac{\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} + \frac{c_0}{2} (p_h - p_b) \right) = 0, \end{array} \right. \quad (36a) $$

$$ \left\{ \begin{array}{l} \sum_{c \in \mathcal{C}} \int_C \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}_i} \int_S [\mathbf{w} \cdot \mathbf{n}^S] (\kappa_0 [\{p_h\}] + \frac{c_0}{2} [\mathbf{u}_h \cdot \mathbf{n}^S]) \\ \qquad + \sum_{S \in \mathcal{S}_{\text{wall}}} \int_S \mathbf{w} \cdot \mathbf{n}^S (\kappa_0 p_h + c_0 \mathbf{u}_h \cdot \mathbf{n}^S) \\ \qquad + \sum_{S \in \mathcal{S}_{SW}} \int_S \mathbf{w} \cdot \mathbf{n}^S (\kappa_0 \frac{p_h + p_b}{2} + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S)) = 0. \end{array} \right. \quad (36b) $$

The following propositions are the discrete version of Proposition 5.

**Proposition 12** (Discrete Hodge-Helmholtz decomposition adapted to the wave system). Denote $(p_h, \mathbf{u}_h) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ the solution of the wave system discretization (24) with Godunov stabilization (18). The Hodge-Helmholtz decomposition $\mathbf{u}_h = \mathbf{u}_h^\varphi + \mathbf{u}_h^\psi$ of Proposition 9 is adapted to the wave system in the sense that, using a forward or a backward Euler method for the time integration, $\mathbf{u}_h^\psi$ is constant in time.

*Proof.* For all $p \in \mathcal{P}_i$, since $[\nabla_{\mathbf{x}}^{\perp} \psi_p \cdot \mathbf{n}^S] = 0$ for all $S \in \mathcal{S}_i$ and $\nabla_{\mathbf{x}}^{\perp} \psi_p \cdot \mathbf{n}^S = 0$ for all $S \in \mathcal{S}_b$, (36b) with forward or backward Euler method for the time integration gives

$$ \forall p \in \mathcal{P}_i \quad \sum_{c \in \mathcal{C}} \int_C \nabla_{\mathbf{x}}^{\perp} \psi_p \cdot \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\delta t} = 0 $$

so that

$$ \forall p \in \mathcal{P}_i \quad \sum_{c \in \mathcal{C}} \int_C \nabla_{\mathbf{x}}^{\perp} \psi_p \cdot \left( \left( \mathbf{u}_h^{\varphi} + \mathbf{u}_h^{\psi} \right)^{n+1} - \left( \mathbf{u}_h^{\varphi} + \mathbf{u}_h^{\psi} \right)^n \right) = 0. $$

Using Proposition 11, it gives

$$ \forall p \in \mathcal{P}_i \quad \sum_{c \in \mathcal{C}} \int_C \nabla_{\mathbf{x}}^{\perp} \psi_p \cdot \left( \left( \mathbf{u}_h^{\psi} \right)^{n+1} - \left( \mathbf{u}_h^{\psi} \right)^n \right) = 0. \quad (37) $$

In the same way, we obtain for all $\partial_k$ connected component of $\partial\Omega$ that

$$ \sum_{c \in \mathcal{C}} \int_C \nabla_{\mathbf{x}}^{\perp} \psi_{\partial_k} \cdot \left( \left( \mathbf{u}_h^{\psi} \right)^{n+1} - \left( \mathbf{u}_h^{\psi} \right)^n \right) = 0. \quad (38) $$

Then, using the characterization of $d\mathbf{S}_h^\psi$ of Proposition 10, we obtain with (37) and (38) that

$$ \left( \mathbf{u}_h^{\psi} \right)^{n+1} - \left( \mathbf{u}_h^{\psi} \right)^n = 0 $$

which ends the proof.

□

**Proposition 13** (Long time limit). We consider the discretization (24) of the wave system (1) with Godunov stabilization (18). The initial conditions are denoted by $(p_h^0, \mathbf{u}_h^0)$. Suppose that the boundary conditions ensure

* $p = p_b$ where $p_b$ is uniform.

* $\mathbf{u} = \mathbf{u}_b$ where $\int_{\partial\Omega} \mathbf{u}_b \cdot \mathbf{n} = 0$.

Then the long time limit exists and is denoted by $(p_h^\infty, \mathbf{u}_h^\infty)$, it satisfies $p_h^\infty$ is uniform (equals to $p_b$) and $\mathbf{u}_h^\infty = \mathbf{u}_h^\psi(0)$ where $\mathbf{u}_h^\psi(0)$ corresponds to the divergence free component of the initial velocity field $\mathbf{u}_h^0$ obtained with Proposition 9. Moreover, defining the relative energy as

$$ E_{\mathbf{u}_h^\psi(0)}(\mathbf{U}_h) := \frac{\rho_0^2 c_0^2}{2} (p_h - p_b)^2 + \frac{1}{2} \left\| \mathbf{u}_h - \mathbf{u}_h^\psi(0) \right\|^2, $$

the volumic average of the relative energy is a Liapunov functional.
Proof. Since for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h$, we have

$$
\begin{align*}
0 &= \sum_{c \in \mathcal{C}} \int_c \mathrm{div}_{\mathbf{x}}(v\mathbf{w}) \\
&= \sum_{S \in \mathcal{S}_i} \int_S [[v\mathbf{w} \cdot \mathbf{n}^S]] + \sum_{S \in \mathcal{S}_b} \int_S v\mathbf{w} \cdot \mathbf{n}^S \\
&= \sum_{S \in \mathcal{S}_i} \int_S [[v]] \{\{ \mathbf{w} \cdot \mathbf{n}^S \}\} + \{\{ v \}\} [[\mathbf{w} \cdot \mathbf{n}^S]] + \sum_{S \in \mathcal{S}_b} \int_S v\mathbf{w} \cdot \mathbf{n}^S,
\end{align*}
$$

so that (37a) can be rewritten as

$$
\begin{align}
& \sum_{c \in \mathcal{C}} \int_c v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S \left( -\frac{1}{\rho_0} \{\{ v \}\} \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket + \frac{c_0}{2} \llbracket v \rrbracket \llbracket p_h \rrbracket \right) \nonumber \\
& \quad - \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S v \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \nonumber \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S v \left( \frac{1}{\rho_0} \frac{-\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} + \frac{c_0}{2} (p_h - p_b) \right) = 0. \tag{40}
\end{align}
$$

Moreover, $(p_b, \mathbf{u}_h^\psi(0)) \in \mathbb{P}_0(\Omega) \times d\mathbf{S}_h^\psi$ satisfies for all $v \in d\mathbb{P}_0(\Omega)$

$$
\sum_{S \in \mathcal{S}_i} \int_S \left( -\frac{1}{\rho_0} \{\{ v \}\} \llbracket \mathbf{u}_h^\psi(0) \cdot \mathbf{n}^S \rrbracket + \frac{c_0}{2} \llbracket v \rrbracket \llbracket p_b \rrbracket \right) = 0, \quad (41)
$$

and for all $\mathbf{w} \in d\mathbf{S}_h$,

$$
\begin{align}
& \sum_{S \in \mathcal{S}_i} \int_S [[\mathbf{w} \cdot \mathbf{n}^S]] \left( \kappa_0 p_b + \frac{c_0}{2} [[\mathbf{u}_h^\psi(0) \cdot \mathbf{n}^S]] \right) + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \kappa_0 p_b + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \kappa_0 p_b \tag{42} \\
&= \sum_{S \in \mathcal{S}_i} \int_S \left( -\kappa_0 \{\{ \mathbf{w} \cdot \mathbf{n}^S \}\} [[p_b]] + \frac{c_0}{2} [[\mathbf{w} \cdot \mathbf{n}^S]] [[\mathbf{u}_h^\psi(0) \cdot \mathbf{n}^S]] \right) \nonumber \\
&=0. \nonumber
\end{align}
$$

Denoting by $p_h^{\text{Rel}} = p_h - p_b$ and $\mathbf{u}_h^{\text{Rel}} = \mathbf{u}_h - \mathbf{u}_h^\psi(0)$ and subtracting (41) to (40) and (42) to (37b), we obtain for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h$,

$$
\left\{
\begin{aligned}
& \sum_{c \in \mathcal{C}} \int_c v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S \left( -\frac{1}{\rho_0} \{\{ v \}\} \llbracket \mathbf{u}_h^{\text{Rel}} \cdot \mathbf{n}^S \rrbracket + \frac{c_0}{2} \llbracket v \rrbracket \llbracket p_h^{\text{Rel}} \rrbracket \right) \\
& \quad - \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S v \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S v \left( \frac{1}{\rho_0} \frac{-\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} + \frac{c_0}{2} (p_h - p_b) \right) = 0 \tag{43a} \\
& \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}_i} \int_S \llbracket \mathbf{w} \cdot \mathbf{n}^S \rrbracket \left( \kappa_0 \{\{ p_h^{\text{Rel}} \}\} + \frac{c_0}{2} \llbracket \mathbf{u}_h^{\text{Rel}} \cdot \mathbf{n}^S \rrbracket \right) \\
& \quad + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S \mathbf{w} \cdot \mathbf{n}^S (\kappa_0 (p_h - p_b) + c_0 \mathbf{u}_h \cdot \mathbf{n}^S) \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \left( \frac{\kappa_0}{2} (p_h - p_b) + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) \right) = 0 \tag{43b}
\end{aligned}
\right.
$$

Multiplying (43a) with $v = p_h^{\text{Rel}}$ by $\rho_0^2 c_0^2$ and adding to (43b) with $\mathbf{w} = \mathbf{u}_h^{\text{Rel}}$, we obtain

$$
\begin{align}
& \sum_{c \in \mathcal{C}} \int_c \partial_{\tau} E_{\mathbf{U}_h^\psi(0)} (\mathbf{U}_h) + \frac{1}{2} \sum_{S \in \mathcal{S}_i} \int_S \rho_0^2 c_0^3 \llbracket p_h^{\text{Rel}} \rrbracket^2 + c_0 \llbracket \mathbf{u}_h^{\text{Rel}} \cdot \mathbf{n}^S \rrbracket^2 \nonumber \\
& \quad + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S c_0 (\mathbf{u}_h \cdot \mathbf{n}^S)^2 \tag{44} \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \frac{\rho_0^2 c_0^3}{2} (p_h - p_b)^2 + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S)^2 = 0. \nonumber
\end{align}
$$

so that

$$
\sum_{c \in \mathcal{C}} \int_c \partial_{\tau} E_{\mathbf{U}_h^\psi(0)} (\mathbf{U}_h) \le 0.
$$
Proof. Since for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h$, we have

$$
\begin{align*}
0 &= \sum_{c \in \mathcal{C}} \int_c \mathrm{div}_{\mathbf{x}}(v\mathbf{w}) \\
&= \sum_{S \in \mathcal{S}_i} \int_S [[v\mathbf{w} \cdot \mathbf{n}^S]] + \sum_{S \in \mathcal{S}_b} \int_S v\mathbf{w} \cdot \mathbf{n}^S \\
&= \sum_{S \in \mathcal{S}_i} \int_S [[v]] \{\{ \mathbf{w} \cdot \mathbf{n}^S \}\} + \{\{ v \}\} [[\mathbf{w} \cdot \mathbf{n}^S]] + \sum_{S \in \mathcal{S}_b} \int_S v\mathbf{w} \cdot \mathbf{n}^S,
\end{align*}
$$

so that (36a) can be rewritten as

$$
\begin{align}
& \sum_{c \in \mathcal{C}} \int_c v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S \left( -\frac{1}{\rho_0} \{\{ v \}\} \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket + \frac{c_0}{2} \llbracket v \rrbracket \llbracket p_h \rrbracket \right) \nonumber \\
& \quad - \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S v \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \nonumber \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S v \left( \frac{1}{\rho_0} \frac{-\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} + \frac{c_0}{2} (p_h - p_b) \right) = 0. \tag{39}
\end{align}
$$

Moreover, $(p_b, \mathbf{u}_h^\psi(0)) \in \mathbb{P}_0(\Omega) \times d\mathbf{S}_h^\psi$ satisfies for all $v \in d\mathbb{P}_0(\Omega)$

$$
\sum_{S \in \mathcal{S}_i} \int_S \left( -\frac{1}{\rho_0} \{\{ v \}\} \llbracket \mathbf{u}_h^\psi(0) \cdot \mathbf{n}^S \rrbracket + \frac{c_0}{2} \llbracket v \rrbracket \llbracket p_b \rrbracket \right) = 0, \quad (40)
$$

and for all $\mathbf{w} \in d\mathbf{S}_h$,

$$
\begin{align}
& \sum_{S \in \mathcal{S}_i} \int_S [[\mathbf{w} \cdot \mathbf{n}^S]] \left( \kappa_0 p_b + \frac{c_0}{2} [[\mathbf{u}_h^\psi(0) \cdot \mathbf{n}^S]] \right) + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \kappa_0 p_b + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \kappa_0 p_b \tag{41} \\
&= \sum_{S \in \mathcal{S}_i} \int_S \left( -\kappa_0 \{\{ \mathbf{w} \cdot \mathbf{n}^S \}\} [[p_b]] + \frac{c_0}{2} [[\mathbf{w} \cdot \mathbf{n}^S]] [[\mathbf{u}_h^\psi(0) \cdot \mathbf{n}^S]] \right) \nonumber \\
&= 0. \nonumber
\end{align}
$$

Denoting by $p_h^{\text{Rel}} = p_h - p_b$ and $\mathbf{u}_h^{\text{Rel}} = \mathbf{u}_h - \mathbf{u}_h^{\psi}(0)$ and subtracting (40) to (39) and (41) to (36b), we obtain for all $v \in d\mathbb{P}_0(\Omega)$ and all $\mathbf{w} \in d\mathbf{S}_h$,

$$
\left\{
\begin{aligned}
& \sum_{c \in \mathcal{C}} \int_c v \partial_{\tau} p_h + \sum_{S \in \mathcal{S}_i} \int_S \left( -\frac{1}{\rho_0} \{\{ v \}\} \llbracket \mathbf{u}_h^{\text{Rel}} \cdot \mathbf{n}^S \rrbracket + \frac{c_0}{2} \llbracket v \rrbracket \llbracket p_h^{\text{Rel}} \rrbracket \right) \\
& \quad - \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S v \frac{1}{\rho_0} \mathbf{u}_h \cdot \mathbf{n}^S \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S v \left( \frac{1}{\rho_0} \frac{-\mathbf{u}_h \cdot \mathbf{n}^S + \mathbf{u}_b \cdot \mathbf{n}^S}{2} + \frac{c_0}{2} (p_h - p_b) \right) = 0 \tag{42a} \\
& \sum_{c \in \mathcal{C}} \int_c \mathbf{w} \cdot \partial_{\tau} \mathbf{u}_h + \sum_{S \in \mathcal{S}_i} \int_S \llbracket \mathbf{w} \cdot \mathbf{n}^S \rrbracket \left( \kappa_0 \{\{ p_h^{\text{Rel}} \}\} + \frac{c_0}{2} \llbracket \mathbf{u}_h^{\text{Rel}} \cdot \mathbf{n}^S \rrbracket \right) \\
& \quad + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \left( \kappa_0 (p_h - p_b) + c_0 \mathbf{u}_h \cdot \mathbf{n}^S \right) \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \mathbf{w} \cdot \mathbf{n}^S \left( \frac{\kappa_0}{2} (p_h - p_b) + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S) \right) = 0 \tag{42b}
\end{aligned}
\right.
$$

Multiplying (42a) with $v = p_h^{\text{Rel}}$ by $\rho_0^2 c_0^2$ and adding to (42b) with $\mathbf{w} = \mathbf{u}_h^{\text{Rel}}$, we obtain

$$
\begin{align}
& \sum_{c \in \mathcal{C}} \int_c \partial_{\tau} E_{\mathbf{U}_h^\psi(0)} (\mathbf{U}_h) + \frac{1}{2} \sum_{S \in \mathcal{S}_i} \int_S \rho_0^2 c_0^3 \llbracket p_h^{\text{Rel}} \rrbracket^2 + c_0 \llbracket \mathbf{u}_h^{\text{Rel}} \cdot \mathbf{n}^S \rrbracket^2 \nonumber \\
& \quad + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S c_0 (\mathbf{u}_h \cdot \mathbf{n}^S)^2 \tag{43} \\
& \quad + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \frac{\rho_0^2 c_0^3}{2} (p_h - p_b)^2 + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S)^2 = 0. \nonumber
\end{align}
$$

so that

$$
\sum_{c \in \mathcal{C}} \int_c \partial_{\tau} E_{\mathbf{U}_h^\psi(0)} (\mathbf{U}_h) \le 0.
$$
Denoting by $(p_h^\infty, \mathbf{u}_h^\infty) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ the long time limit, (44) gives

$$ 
\left\{ 
\begin{array}{lll}
\llbracket p_h^\infty \rrbracket = 0 & \text{and} & \llbracket \mathbf{u}_h^\infty \cdot \mathbf{n}^S \rrbracket = 0, \quad \forall S \in \mathcal{S}_i, \\
\mathbf{u}_h^\infty \cdot \mathbf{n}^S = 0, & & \forall S \in \mathcal{S}_{\text{wall}}, \\
p_h^\infty = p_b & \text{and} & \mathbf{u}_h^\infty \cdot \mathbf{n}^S = \mathbf{u}_b \cdot \mathbf{n}^S, \quad \forall S \in \mathcal{S}_{\text{SW}},
\end{array}
\right.
$$

so that $p_h^\infty = p_b$ and $\mathbf{u}_h^\infty = \mathbf{u}_h^\psi(0)$ by uniqueness of the Hodge decomposition of Proposition 9. $\square$

# 4 Discretization for the barotropic Euler system

The isentropic Euler system is shortly noted as

$$ \partial_t \mathbf{W} + \nabla \cdot \mathbf{f}(\mathbf{W}) = 0 \qquad (45) $$

where $\mathbf{W} = (\rho, \rho\mathbf{u})^T$ is the vector of conservative variables, $\rho$ is the density and $\mathbf{u}$ the velocity. $\mathbf{f}$ is the flux defined as

$$ \mathbf{f}(\mathbf{W}) = \begin{pmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + pI \end{pmatrix}. $$

$p$ is the pressure, and is linked with the density through the equation of state: $p = p(\rho)$, which is supposed to be convex and strictly increasing. The sound velocity $c$ is defined as $\sqrt{p'(\rho)}$, and the Mach number is defined as $|\mathbf{u}|/c$.

The numerical solutions of (45) will be searched in the finite element space

$$ \mathcal{V}_h = d\mathbb{P}_0 \times d\mathbf{S}_h. $$

We consider the discontinuous Galerkin method for (45):

$$ \text{Find } \mathbf{W}_h \in \mathcal{V}_h \quad \forall \varphi \in \mathcal{V}_h \quad \sum_{c \in \mathcal{C}_i} \int_c (\varphi \cdot \partial_t \mathbf{W}_h - \mathbf{f}(\mathbf{W}_h) \cdot \nabla \varphi) \\ + \sum_{S \in \mathcal{S}_i} \int_S \llbracket \varphi \rrbracket \cdot \tilde{\mathbf{f}}(\mathbf{W}_h, \mathbf{n}^S) + \sum_{S \in \mathcal{S}_b} \int_S \varphi^{\text{Left}} \cdot \mathbf{f}^b(\mathbf{W}_h) \cdot \mathbf{n}^S = 0. \qquad (46) $$

The numerical flux $\tilde{\mathbf{f}}(\mathbf{W}_h, \mathbf{n}^S)$ may be any of known numerical flux [42]: Roe, Lax-Friedrich, HLL, HLLC, exact Godunov... Note that in the case of the wave system (16) and (25), the cell integration term vanishes, because of the linearity of the flux and the fact that the flux of velocity of (1) includes only a gradient, which leads to a weak formulation involving only the divergence of the test functions which is 0 for all test functions. This is no longer the case for the barotropic Euler system, and the cell integral must be taken into account.

The link between the behaviour of (45) at low Mach number and the wave system was extensively discussed in [27, 29]. The main result, which may be obtained by a two-scale asymptotic expansion of the numerical scheme (46) [29, Section 3.2], is the classification of the numerical flux for Euler with respect to the possible numerical flux for the wave system:

* A first family of flux is asymptotically consistent with the Godunov' flux for the wave system. This family includes for example the HLLC scheme, the Roe scheme, and the Osher scheme. For this family of numerical flux, Proposition 12 ensures the preservation of the divergence free component of the discrete Hodge-Helmholtz decomposition and Proposition 13 ensures to reach the right long time limit (according to continuous case of Proposition 5) for the wave system. Then, the numerical scheme (46) is expected to be low Mach number accurate.

* A second family of flux is asymptotically consistent with the Rusanov scheme for the wave system. This family includes for example the Rusanov, HLL or Lax-Friedrich flux for the barotropic Euler system. In this case, the numerical scheme (46) is expected to be inaccurate at low Mach number.

# 5 Numerical results

In this section, numerical tests are performed to illustrate the results of section 2 and section 3 concerning the discretization (20) of the wave system (1) obtained the new approximation space $d\mathbf{S}_h$. Finally, the low Mach number behavior of the Euler discretization (46) using the new approximation space is studied.
Denoting by $(p_h^\infty, u_h^\infty) \in d\mathbb{P}_0(\Omega) \times d\mathbf{S}_h(\Omega)$ the long time limit, (43) gives

$$
\left\{
\begin{array}{ll}
\llbracket p_h^\infty \rrbracket = 0 & \text{and} \quad \llbracket u_h^\infty \cdot \mathbf{n}^S \rrbracket = 0, \quad \forall S \in \mathcal{S}_i, \\
u_h^\infty \cdot \mathbf{n}^S = 0, & \text{and} \quad \forall S \in \mathcal{S}_{\text{wall}} \\
p_h^\infty = p_b & \text{and} \quad u_h^\infty \cdot \mathbf{n}^S = u_b \cdot \mathbf{n}^S, \quad \forall S \in \mathcal{S}_{\text{SW}},
\end{array}
\right.
$$

so that $p_h^\infty = p_b$ and $u_h^\infty = u_h^\psi(0)$ by uniqueness of the Hodge decomposition of Proposition 9. $\square$

4 Discretization for the barotropic Euler system

The isentropic Euler system is shortly noted as

$$
\partial_t \mathbf{W} + \nabla \cdot \mathbf{f}(\mathbf{W}) = 0 \qquad (44)
$$

where $\mathbf{W} = (\rho, \rho\mathbf{u})^T$ is the vector of conservative variables, $\rho$ is the density and $\mathbf{u}$ the velocity. $\mathbf{f}$ is the flux defined as

$$
\mathbf{f}(\mathbf{W}) = \begin{pmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + pI \end{pmatrix}.
$$

$p$ is the pressure, and is linked with the density through the equation of state: $p = p(\rho)$, which is supposed to be convex and strictly increasing. The sound velocity $c$ is defined as $\sqrt{p'(\rho)}$, and the Mach number is defined as $|\mathbf{u}|/c$.

The numerical solutions of (44) will be searched in the finite element space

$$
V_h = d\mathbb{P}_0 \times d\mathbf{S}_h.
$$

We consider the discontinuous Galerkin method for (44):

$$
\text{Find } \mathbf{W}_h \in \mathcal{V}_h \quad \forall \varphi \in \mathcal{V}_h \quad \sum_{c \in \mathcal{C}_i} \int_c (\varphi \cdot \partial_t \mathbf{W}_h - \mathbf{f}(\mathbf{W}_h) \cdot \nabla \varphi) \\
+ \sum_{S \in \mathcal{S}_i} \int_S \llbracket \varphi \rrbracket \cdot \tilde{\mathbf{f}}(\mathbf{W}_h, \mathbf{n}^S) + \sum_{S \in \mathcal{S}_b} \int_S \varphi^{\text{Left}} \cdot \mathbf{f}^b(\mathbf{W}_h) \cdot \mathbf{n}^S = 0. \tag{45}
$$

The numerical flux $\tilde{f}(\mathbf{W}_h, \mathbf{n}^S)$ may be any of known numerical flux [38]: Roe, Lax-Friedrich, HLL, HLLC, exact Godunov... Note that in the case of the wave system (16) and (24), the cell integration term vanishes, because of the linearity of the flux and the fact that the flux of velocity of (1) includes only a gradient, which leads to a weak formulation involving only the divergence of the test functions which is 0 for all test functions. This is no longer the case for the barotropic Euler system, and the cell integral must be taken into account.

The link between the behaviour of (44) at low Mach number and the wave system was exten-
sively discussed in [25, 27]. The main result, which may be obtained by a two-scale asymptotic
expansion of the numerical scheme (45) [27, Section 3.2], is the classification of the numerical flux
for Euler with respect to the possible numerical flux for the wave system:

* A first family of flux is asymptotically consistent with the Godunov' flux for the wave system. This family includes for example the HLLC scheme, the Roe scheme, and the Osher scheme. For this family of numerical flux, Proposition 12 ensures the preservation of the divergence free component of the discrete Hodge-Helmholtz decomposition and Proposition 13 ensures to reach the right long time limit (according to continuous case of Proposition 5) for the wave system. Then, the numerical scheme (45) is expected to be low Mach number accurate.

* A second family of flux is asymptotically consistent with the Rusanov scheme for the wave system. This family includes for example the Rusanov, HLL or Lax-Friedrich flux for the barotropic Euler system. In this case, the numerical scheme (45) is expected to be inaccurate at low Mach number.

5 Numerical results

In this section, numerical tests are performed to illustrate the results of section 2 and section 3 concerning the discretization (20) of the wave system (1) obtained the new approximation space dS<sub>h</sub>. Finally, the low Mach number behavior of the Euler discretization (45) using the new approximation space is studied.
From a numerical point of view, cells integrals of (20) and (46) are performed using a 2 × 2 = 4 points Gauss' quadrature formula while side integrals of (20) and (46) are performed using midpoint approximation. Concerning time integration, all the results presented below were obtained with a forward Euler scheme. Last, all the initializations are performed with a L² projection on the finite element space. Note that none of the propositions (and especially the main ones, namely Proposition 4, Proposition 12, Proposition 13) depend on the way the initial condition is done.

## 5.1 Order of accuracy on the adjoint curl

The aim of this test case is to assess the order of approximation given by the adjoint curl, $(\mathbf{curl}_x)^*$, defined in (24). For this, we first consider a function $f$ defined in polar coordinate $r \mapsto f(r)$. We then consider the vector field defined as

$$ \mathbf{u} = \mathbf{curl}_x f = \frac{\partial f}{\partial r} \mathbf{e}_{\theta}. $$

We have then

$$ \mathrm{rot}_{\mathbf{x}} \mathbf{u} = -\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial f}{\partial r} \right). $$

For assessing the order, we choose a function that is regular and has a compact support in $[0, 1]^2$:

$$ f : r \mapsto f(r) = \begin{cases} r_0 e^{-\frac{r}{1-(r/r_0)^2}} & \text{if } r \le r_0 \\ 0 & \text{otherwise,} \end{cases} $$

where $r$ is the $L^2$ distance defined with respect to a point $[x_c, y_c] \in [0, 1]^2$, and $r_0$ is such that the circle of center $[x_c, y_c]$ and of radius $r_0$ is strictly inside $[0, 1]^2$. Then, denoting by $\bar{r} = \frac{r}{r_0}$,

$$ \mathbf{u} = \begin{pmatrix} \frac{2\alpha y}{r_0} \frac{e^{-\frac{\alpha}{1-\bar{r}^2}}}{(1-\bar{r}^2)^2} \\ -\frac{2\alpha x}{r_0} \frac{e^{-\frac{\alpha}{1-\bar{r}^2}}}{(1-\bar{r}^2)^2} \end{pmatrix}, \qquad (47) $$

and

$$ \mathrm{rot}_{\mathbf{x}} \mathbf{u} = \frac{-4\alpha (\alpha \bar{r}^2 + \bar{r}^4 - 1) e^{-\frac{\alpha}{1-\bar{r}^2}}}{r_0^2 (1-\bar{r}^2)^4}. \qquad (48) $$

The numerical test is led as follows. We first project $\mathbf{u}$ of (47) on the finite element space $d\mathbf{S}_h$, for finding $\mathbf{u}_h$. Then $(\mathbf{curl}_x)^* \mathbf{u}_h \in \mathbb{Q}_1$ is computed using (24). Last, the error with respect to the exact solution (48) is computed. The computation is led on four regular Cartesian meshes with 10 × 10, 20 × 20, 40 × 40, 80 × 80, and 160 × 160 cells. The data used are $\alpha = 4$ (for avoiding a too sharp behaviour near $r_0$), $r_0 = 0.15$, and $x_c = y_c = 0.5$. The error obtained is presented in Table 1, and shows a second order of accuracy of the adjoint curl. This second order of accuracy of the adjoint curl first surprised us, because the projection of the vector field $\mathbf{u}$ is done on a piecewise constant approximation space, for which the order of accuracy is 1, and we could expect that the error on the adjoint curl be bounded by the initial projection. Still, the approximation space of $(\mathbf{curl}_x)^* \mathbf{u}$ is piecewise linear, and so a second order of approximation can be obtained. The observed order of accuracy probably comes from a commutation property between the projection operators and the continuous and discrete $(\mathbf{curl}_x)^* \mathbf{u}$. The study of this commutation property is out of the scope of this paper.

## 5.2 Wave equation

### 5.2.1 Periodic vortex

The aim of this section is to assess the result of Proposition 4. For this, we consider the following vortex initial solution on the quad $[0, 1]^2$

$$ p_0 = 0, \qquad \mathbf{u}_0 = \begin{pmatrix} -\frac{y}{r_0} e^{-\bar{r}^2} \\ \frac{x}{r_0} e^{-\bar{r}^2} \end{pmatrix}, $$
<table><caption>Table 1: Order of accuracy obtained after a projection on $dS_h$ and a computation of the adjoint curl $(\nabla^{\perp})^*$.</caption><thead><tr><th>h</th><th>Error</th><th>Local order</th></tr></thead><tbody><tr><td>0.1</td><td>0.111664</td><td>—</td></tr><tr><td>0.05</td><td>0.0588788</td><td>0.923</td></tr><tr><td>0.025</td><td>0.0159502</td><td>1.884</td></tr><tr><td>0.0125</td><td>0.00379819</td><td>2.070</td></tr><tr><td>0.00625</td><td>0.000855902</td><td>2.149</td></tr><tr><td>0.003125</td><td>0.000208047</td><td>2.040</td></tr></tbody></table>

From a numerical point of view, cells integrals of (20) and (45) are performed using a $2 \times 2 = 4$ points Gauss' quadrature formula while side integrals of (20) and (45) are performed using midpoint approximation. Concerning time integration, all the results presented below were obtained with a forward Euler scheme.

## 5.1 Order of accuracy on the adjoint curl

The aim of this test case is to assess the order of approximation given by the adjoint curl $(\nabla^{\perp})^*$ defined in (23). For this, we first consider a function $f$ defined in polar coordinate $r \mapsto f(r)$. We then consider the vector field defined as

$$ \mathbf{u} = \nabla^{\perp} f = \frac{\partial f}{\partial r} \mathbf{e}_{\theta}. $$

We have then

$$ \nabla^{\perp} \cdot \mathbf{u} = -\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial f}{\partial r} \right). $$

For assessing the order, we choose a function that is regular and has a compact support in $[0, 1]^2$:

$$ f : r \mapsto f(r) = \begin{cases} r_0 e^{-\frac{r}{1-\left(\frac{r}{r_0}\right)^2}} & \text{if } r \le r_0 \\ 0 & \text{otherwise,} \end{cases} $$

where $r$ is the $L^2$ distance defined with respect to a point $[x_c, y_c] \in [0, 1]^2$, and $r_0$ is such that the circle of center $[x_c, y_c]$ and of radius $r_0$ is strictly inside $[0, 1]^2$. Then, denoting by $\bar{r} = \frac{r}{r_0}$,

$$ \mathbf{u} = \begin{pmatrix} \frac{2\alpha y}{r_0} \frac{e^{-\frac{\alpha}{1-\bar{r}^2}}}{(1-\bar{r}^2)^2} \\ -\frac{2\alpha x}{r_0} \frac{e^{-\frac{\alpha}{1-\bar{r}^2}}}{(1-\bar{r}^2)^2} \end{pmatrix}, \quad (46) $$

and

$$ \nabla^{\perp} \cdot \mathbf{u} = \frac{-4\alpha (\alpha \bar{r}^2 + \bar{r}^4 - 1) e^{-\frac{\alpha}{1-\bar{r}^2}}}{r_0^2 (1-\bar{r}^2)^4}. \quad (47) $$

The numerical test is led as follows. We first project $\mathbf{u}$ of (46) on the finite element space $dS_h$, for finding $\mathbf{u}_h$. Then $(\nabla^{\perp})^* \mathbf{u}_h \in \mathbb{Q}_1$ is computed using (23). Last, the error with respect to the exact solution (47) is computed. The computation is led on four regular Cartesian meshes with $10 \times 10$, $20 \times 20$, $40 \times 40$, $80 \times 80$, and $160 \times 160$ cells. The data used are $\alpha = 4$ (for avoiding a too sharp behaviour near $r_0$), $r_0 = 0.15$, and $x_c = y_c = 0.5$. The error obtained is presented in Table 1, and shows a second order of accuracy of the adjoint curl.

## 5.2 Wave equation

### 5.2.1 Periodic vortex

The aim of this section is to assess the result of Proposition 4. For this, we consider the following vortex initial solution on the quad $[0, 1]^2$

$$ p_0 = 0, \quad \mathbf{u}_0 = \begin{pmatrix} -\frac{y}{r_0} e^{-\bar{r}^2} \\ \frac{x}{r_0} e^{-\bar{r}^2} \end{pmatrix}, $$
<table><caption>Table 1: Order of accuracy obtained after a projection on $dS_h$ and a computation of the adjoint curl, $(\mathbf{curl}_x)^*\mathbf{u}$.</caption><thead><tr><th>h</th><th>Error</th><th>Local order</th></tr></thead><tbody><tr><td>0.1</td><td>0.111664</td><td>—</td></tr><tr><td>0.05</td><td>0.0588788</td><td>0.923</td></tr><tr><td>0.025</td><td>0.0159502</td><td>1.884</td></tr><tr><td>0.0125</td><td>0.00379819</td><td>2.070</td></tr><tr><td>0.00625</td><td>0.000855902</td><td>2.149</td></tr><tr><td>0.003125</td><td>0.000208047</td><td>2.040</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 2: Quadrangular meshes used for the conservation of the adjoint of the curl test described in subsection 5.2.1: Cartesian on the left and quadrangular unstructured on the right.</figcaption></figure>

with the same notations as in the previous subsection, and with periodic boundary conditions. The vortex is centered around $(x_c, y_c) = (0.5, 0.5)$, and we take $r_0 = 0.15$. The test is run until time $t = 1$ on a Cartesian mesh of $10 \times 10$ cells, and on an unstructured quadrangular mesh with similar resolution represented in Figure 2. At each time step, the difference between the adjoint curl of the velocity, $(\mathbf{curl}_x)^*\mathbf{u}$ and its initial value $(\mathbf{curl}_x)^*\mathbf{u}_0$ is computed. The $L^2$ norm of this difference is plotted in Figure 3. Note that $(d\mathbb{Q}_0)^2 \subset dS_h$, so that the $(\mathbf{curl}_x)^*$ defined in (24) is also defined in $(d\mathbb{Q}_0)^2$. As proven in Proposition 4, $(\mathbf{curl}_x)^*\mathbf{u}$ is conserved when the approximation space is $dS_h$ and the numerical flux is the Godunov' flux. In the other cases, namely either when the velocity approximation space is $(d\mathbb{Q}_0)^2$ or when the Lax-Friedrich numerical flux is used, the adjoint curl is not preserved.

### 5.2.2 Cylinder scattering

The purpose of this test case is to illustrate Proposition 13. The domain is an annulus $[r_0, r_1] \times [0; 2\pi]$ where $r_0 = 0.5$ and $r_1 = 5.5$. Wall boundary condition (26) is applied in $r = r_0$ and

<figure><img src="image_3.png" /><figcaption>Figure 3: Plot of $L^2$ norm of the difference between the adjoint curl of the velocity $(\mathbf{curl}_x)^*\mathbf{u}$ and the initial adjoint curl of the velocity $(\mathbf{curl}_x)^*\mathbf{u}_0$ for the periodic test case described in subsection 5.2.1. The left plot matches with a Cartesian mesh and the right plot with a quadrangular unstructured mesh (the two meshes are represented in Figure 2).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Quadrangular meshes used for the conservation of the adjoint of the curl test described in subsection 5.2.1: Cartesian on the left and quadrangular unstructured on the right.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Plot of $L^2$ norm of the difference between the adjoint curl of the velocity $(\nabla^{\perp})^* \mathbf{u}$ and the initial adjoint curl of the velocity $(\nabla^{\perp})^* \mathbf{u}_0$ for the periodic test case described in subsection 5.2.1. The left plot matches with a Cartesian mesh and the right plot with a quadrangular unstructured mesh (the two meshes are represented in Figure 2).</figcaption></figure>

with the same notations as in the previous subsection, and with periodic boundary conditions. The vortex is centered around $(x_c, y_c) = (0.5, 0.5)$, and we take $r_0 = 0.15$. The test is run until time $t = 1$ on a Cartesian mesh of $10 \times 10$ cells, and on an unstructured quadrangular mesh with similar resolution represented in Figure 2. At each time step, the difference between the adjoint curl of the velocity, $(\nabla^{\perp})^* \mathbf{u}$ and its initial value $(\nabla^{\perp})^* \mathbf{u}_0$ is computed. The $L^2$ norm of this difference is plotted in Figure 3. As proven in Proposition 4, $(\nabla^{\perp})^* \mathbf{u}$ is conserved when the approximation space is $d\mathbf{S}_h$ and the numerical flux is the Godunov' flux. In the other cases, namely either when the velocity approximation space is $(d\mathbb{Q}_0)^2$ or when the Lax-Friedrich numerical flux is used, the adjoint curl is not preserved.

### 5.2.2 Cylinder scattering

The purpose of this test case is to illustrate Proposition 13. The domain is an annulus $[r_0, r_1] \times [0; 2\pi[$ where $r_0 = 0.5$ and $r_1 = 5.5$. Wall boundary condition (25) is applied in $r = r_0$ and inlet/outlet boundary condition (26) is applied in $r = r_1$ with $p_b = 0$ and $\mathbf{u}_b = (1, 0)^T$. Initial data are given by $p^0 = 0$ and $\mathbf{u}^0 = 0$.

In Figure 4, we plot the relative energy residual defined in (43) by

$$ \frac{1}{2} \sum_{S \in \mathcal{S}_i} \int_S (\rho_0^2 c_0^3 [p_h]^2 + c_0 [\mathbf{u}_h \cdot \mathbf{n}^S]^2) + \sum_{S \in \mathcal{S}_{\text{wall}}} \int_S c_0 (\mathbf{u}_h \cdot \mathbf{n}^S)^2 \\ + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \frac{\rho_0^2 c_0^3}{2} (p_h - p_b)^2 + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S)^2, $$

as a function of the computational time for the numerical solution $(p_h, \mathbf{u}_h)$. A mesh containing $10 \times 20$ cells, where the first number corresponds to radial discretization and the second one to orthoradial discretization, is used. Using $(d\mathbb{Q}_0)^2$ or the new approximation space $d\mathbf{S}_h$, the numerical schemes converge to a long-time limit $(p_h^\infty, \mathbf{u}_h^\infty)$. This result was proved in [26] for $(d\mathbb{Q}_0)^2$ with Rusanov and Godunov stabilization. Concerning $d\mathbf{S}_h$, no numerical difficulty was observed.
<figure><img src="image_4.png" /><figcaption>Figure 4: Wave cylinder scattering - Energy residual as a function of the computational time.</figcaption></figure>

inlet/outlet boundary condition (27) is applied in $r = r_1$ with $p_b = 0$ and $\mathbf{u}_b = (1, 0)^T$. Initial data are given by $p^0 = 0$ and $\mathbf{u}^0 = 0$.

In Figure 4, we plot the relative energy residual defined in (44) by

$$ \frac{1}{2} \sum_{S \in \mathcal{S}_i} \int_S (\rho_0^2 c_0^3 [\mathbf{p}_h]^2 + c_0 [\mathbf{u}_h \cdot \mathbf{n}^S]^2) + \sum_{S \in \mathcal{S}_{\text{Wall}}} \int_S c_0 (\mathbf{u}_h \cdot \mathbf{n}^S)^2 \\ + \sum_{S \in \mathcal{S}_{\text{SW}}} \int_S \frac{\rho_0^2 c_0^3}{2} (p_h - p_b)^2 + \frac{c_0}{2} (\mathbf{u}_h \cdot \mathbf{n}^S - \mathbf{u}_b \cdot \mathbf{n}^S)^2, $$

as a function of the computational time for the numerical solution $(p_h, \mathbf{u}_h)$. A mesh containing 10 × 20 cells, where the first number corresponds to radial discretization and the second one to orthoradial discretization, is used. Using $(d\mathbb{Q}_0)^2$ or the new approximation space $d\mathbf{S}_h$, the numerical schemes converge to a long-time limit $(p_h^\infty, \mathbf{u}_h^\infty)$. This result was proved in [28] for $(d\mathbb{Q}_0)^2$ with Rusanov and Godunov stabilization. Concerning $d\mathbf{S}_h$, no numerical difficulty was observed. In agreement with the result of Proposition 13, the relative energy residual of the long time limit obtained with $d\mathbf{S}_h$ and Godunov stabilization is zero such that the long time limit satisfies $p_h^\infty = p_b$ and $\mathbf{u}_h^\infty = \mathbf{u}_b^\psi(0)$ by uniqueness of the Hodge-Helmholtz decomposition of Proposition 9. Using $(d\mathbb{Q}_0)^2$ or $d\mathbf{S}_h$ with Rusanov stabilization, the energy residual of the long time limit is not zero.

For this test case, the exact long time wave solution can be computed analytically and is given by

$$ \left\{ \begin{array}{l} p_{\text{ex}}^{\infty}(r, \theta) = p_b = 0 \\ \mathbf{u}_{\text{ex}}^{\infty}(r, \theta) = \mathbf{u}^{\psi}(0) = \frac{r_1^2}{r_1^2 - r_0^2} \begin{pmatrix} 1 - \frac{r_0^2}{r_1^2} \cos(2\theta) \\ -\frac{r_0^2}{r_1^2} \sin(2\theta) \end{pmatrix}. \end{array} \right. \qquad (49) $$

In Figure 5, a mesh convergence study is performed. The computation is performed on five meshes containing respectively 6 × 12, 10 × 20, 20 × 40, 40 × 80 and 80 × 160 cells. Using the new approximation space $d\mathbf{S}_h$ with Godunov stabilization, the long time pressure is exactly zero and the long time velocity converges to the exact solution with a rate of one. This means that we have a commutation between the mesh convergence and the long time limit convergence. In the two other approximation type, namely with $(d\mathbb{Q}_0)^2$ with the Godunov' scheme, and with $d\mathbf{S}_h$ and the Rusanov scheme, we know that for each mesh, we have a long time limit as proven in [28], but the long time limit does not converge to the correct solution when the mesh is refined. This means that in these cases, the limits in long time and the mesh convergence do not commute.

The results found in Figure 5 deserves additional comments. In [29, Section 5.1.1], the numerical results obtained on a different test case show that the order of accuracy obtained on quads is generally one order lower than the optimal one. This means that a first order of accuracy should be obtained with the classical DG¹ numerical method, still at the price of having eight degrees of freedom per cell for the velocity (four by velocity component), whereas the numerical method proposed here has only three degrees of freedom per cell for the velocity. The enrichment procedure for the velocity space was further discussed in the two submitted article [38, 37] for the higher order case.
<figure><img src="image_6.png" /><figcaption>Figure 4: Wave cylinder scattering - Energy residual as a function of the computational time.</figcaption></figure>

In agreement with the result of Proposition 13, the relative energy residual of the long time limit obtained with dS_h and Godunov stabilization is zero such that the long time limit satisfies $p_h^\infty = p_b$ and $u_h^\infty = u_h^\psi(0)$ by uniqueness of the Hodge-Helmholtz decomposition of Proposition 9. Using $(d\mathbb{Q}_0)^2$ or dS_h with Rusanov stabilization, the energy residual of the long time limit is not zero.

For this test case, the exact long time wave solution can be computed analytically and is given by

$$ \left\{ \begin{array}{l} p_{\text{ex}}^{\infty}(r, \theta) = p_b = 0 \\ \mathbf{u}_{\text{ex}}^{\infty}(r, \theta) = \mathbf{u}^{\psi}(0) = \frac{r_1^2}{r_1^2 - r_0^2} \left( \begin{array}{l} 1 - \frac{r_0^2}{r_1^2} \cos(2\theta) \\ -\frac{r_0^2}{r_1^2} \sin(2\theta) \end{array} \right). \end{array} \right. \quad (48) $$

In Figure 5, a mesh convergence study is performed. The computation is performed on five meshes containing respectively 6 × 12, 10 × 20, 20 × 40, 40 × 80 and 80 × 160 cells. Using the new approximation space dS_h with Godunov stabilization, the long time pressure is exactly zero and the long time velocity converges to the exact solution with a rate of one.

## 5.3 Euler equation

### 5.3.1 Cylinder scattering

We illustrate the good behavior at low Mach numbers of the numerical scheme obtained with the new approximation space $(\rho, \rho\mathbf{u})_h \in d\mathbb{P}_0 \times dS_h$ and a Roe numerical flux. Note that similar results are obtained with any other numerical flux that degenerates to a Godunov stabilization for the asymptotic wave system. The domain and meshes used are the same as in subsubsection 5.2.2. Wall boundary condition is applied in $r = r_0$ while Steger-Warming boundary condition is applied in $r = r_1$ with a state characterized by its density at infinity $\rho_b$ and its velocity at infinity $\mathbf{u}_b = (u_b, 0)^T$. For all the computations, $\rho_b$ is set to $\rho_b = 2$ and $u_b$ is deduced from the Mach number $M_b$ by $u_b = M_b \sqrt{p'(\rho_b)}$. Initial data are uniform and set equal to $\rho^0 = \rho_b$ and $\mathbf{u}^0 = 0$. The exact incompressible solution is given by

$$ \left\{ \begin{array}{l} \mathbf{u}_{\text{ex}}^{\text{Incomp}}(r, \theta) = u_b \mathbf{u}_{\text{ex}}^{\infty}(r, \theta), \\ \rho_{\text{ex}}^{\text{Incomp}}(r, \theta) = \rho_b + \rho_b \left( \frac{r_1^2}{r_1^2 - r_0^2} \right)^2 \left( \frac{r_0^2}{r_1^2} \cos(2\theta) - \frac{r_0^4}{2r_1^4} \right) M_b^2, \end{array} \right. $$

where $\mathbf{u}_{\text{ex}}^\infty$ is given by (48).

In Figure 6, the iso-contours of the steady velocity field are plotted for $M_b = 10^{-4}$ with the mesh containing 20 × 40 cells. The solution obtained with the new approximation space dS_h and Roe' numerical flux seems to be in agreement with the exact incompressible solution. However, the study of isolines is not sufficient to conclude on low Mach number behavior. A more detailed study needs to be carried out not only on the mesh convergence of the compressible low Mach number solution to the incompressible solution, but also on the order of magnitude of the density gradient and the momentum divergence with respect to the Mach number.

In Figure 7, a mesh convergence at Mach number $M_b = 10^{-4}$ is performed. The steady solution obtained with the new approximation space dS_h and Roe' numerical flux converges towards the incompressible exact solution with a rate close to one.
<figure><img src="image_5.png" /><figcaption>Figure 5: Wave cylinder scattering - L² norm of the error between the exact long time solution and the numerical long time solution. A log-log plot is used.</figcaption></figure>

## 5.3 Euler equation

### 5.3.1 Cylinder scattering

We illustrate the good behavior at low Mach numbers of the numerical scheme obtained with the new approximation space $(\rho, \rho\mathbf{u})_h \in d\mathbb{P}_0 \times d\mathbf{S}_h$ and a Roe numerical flux. Note that similar results are obtained with any other numerical flux that degenerates to a Godunov stabilization for the asymptotic wave system. The domain and meshes used are the same as in subsection 5.2.2. Wall boundary condition is applied in $r = r_0$ while Steger-Warming boundary condition is applied in $r = r_1$ with a state characterized by its density at infinity $\rho_b$ and its velocity at infinity $\mathbf{u}_b = (u_b, 0)^T$. For all the computations, $\rho_b$ is set to $\rho_b = 2$ and $u_b$ is deduced from the Mach number $M_b$ by $u_b = M_b \sqrt{p'(\rho_b)}$. Initial data are uniform and set equal to $\rho^0 = \rho_b$ and $\mathbf{u}^0 = 0$. The exact incompressible solution is given by

$$ \left\{ \begin{array}{l} \mathbf{u}_{\text{ex}}^{\text{Incomp}}(r, \theta) = u_b \mathbf{u}_{\text{ex}}^{\infty}(r, \theta), \\ \rho_{\text{ex}}^{\text{Incomp}}(r, \theta) = \rho_b + \rho_b \left( \frac{r_1^2}{r_1^2 - r_0^2} \right)^2 \left( \frac{r_0^2}{r^2} \cos(2\theta) - \frac{r_0^4}{2r^4} \right) M_b^2, \end{array} \right. $$

where $\mathbf{u}_{\text{ex}}^{\infty}$ is given by (49).

In Figure 6, the iso-contours of the steady velocity field are plotted for $M_b = 10^{-4}$ with the mesh containing 20 × 40 cells. The solution obtained with the new approximation space $d\mathbf{S}_h$ and Roe' numerical flux seems to be in agreement with the exact incompressible solution. However, the study of isolines is not sufficient to conclude on low Mach number behavior. A more detailed study needs to be carried out not only on the mesh convergence of the compressible low Mach number solution to the incompressible solution, but also on the order of magnitude of the density gradient and the momentum divergence with respect to the Mach number.

In Figure 7, a mesh convergence at Mach number $M_b = 10^{-4}$ is performed. The steady solution obtained with the new approximation space $d\mathbf{S}_h$ and Roe' numerical flux converges towards the incompressible exact solution with a rate close to one.

In Figure 8, we study the order of the dimensionless density gradient and the dimensionless momentum divergence with respect to the Mach number. Then, we plot the semi-norms $|\tilde{\rho}_h|_{H^1}$
<figure><img src="image_7.png" /><figcaption>Figure 5: Wave cylinder scattering - $L^2$ norm of the error between the exact long time solution and the numerical long time solution. A log-log plot is used.</figcaption></figure>

In Figure 8, we study the order of the dimensionless density gradient and the dimensionless momentum divergence with respect to the Mach number. Then, we plot the semi-norms $|\tilde{\rho}_h|_{H^1}$ and $|\tilde{\rho}\tilde{u}_h|_{H^{\text{div}}}$, defined by

$$ 
\begin{aligned}
|p_h|^2_{H^1} &= \sum_{c \in \mathcal{C}} \int_c \|\nabla p_h\|^2 + \sum_{S \in \mathcal{S}_i} \int_S \llbracket p_h \rrbracket^2 \\
|\mathbf{u}_h|^2_{H^{\text{div}}} &= \sum_{c \in \mathcal{C}} \int_c (\text{div}_{\mathbf{x}} \mathbf{u}_h)^2 + \sum_{S \in \mathcal{S}_i} \int_S \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket^2
\end{aligned}
$$

as a function of the Mach number for a fixed mesh containing 10 × 20 cells. Using the new approximation space $d\mathbf{S}_h$ and Roe' numerical flux, we observe that $|\tilde{\rho}_h|_{H^1} = \mathcal{O}(M^2)$ and $|\tilde{\rho}\tilde{u}_h|_{H^{\text{div}}} = \mathcal{O}(M)$ and so

$$ \nabla \tilde{\rho} = \mathcal{O}(M^2) \quad \text{and} \quad \text{div}_{\mathbf{x}}(\tilde{\rho}\tilde{u}) = \mathcal{O}(M) $$

as expected.

### 5.3.2 Propagation of a low Mach number acoustic wave over a steady vortex

In this test case, we evaluate the ability of the numerical scheme to handle both incompressible phenomena and acoustic wave propagation at low Mach numbers. This test case was proposed in [11], because it had been remarked that some low Mach number fixes spoil the propagation of acoustic waves. The domain is the rectangle $[-0.1; 1.1] \times [0; 1]$. Periodic conditions are applied to top and bottom boundaries, and input-output conditions are used for left and right boundaries. The vortex is centered around $(x_c, y_c) = (0.5, 0.5)$, characterized by its reference Mach number $M_{\text{ref}}$, and is given by

$$ 
\left\{ 
\begin{array}{l}
\rho(\mathbf{x}) = \rho_0 \left( 1 - \frac{M_{\text{ref}}^2}{\lambda_{\text{max}}^2} e^{-\frac{2\alpha^2}{1-\bar{r}^2}} \right), \\
\mathbf{u}(\mathbf{x}) = u_0 \left( \begin{array}{l} y \\ -x \end{array} \right) \frac{2\alpha}{\lambda_{\text{max}} r_0 (1 - \bar{r}^2)} e^{-\frac{\alpha^2}{1-\bar{r}^2}},
\end{array}
\right. 
$$
<figure><img src="image_6.png" /><figcaption>Figure 6: Euler cylinder scattering - Iso-contours of the norm of the velocity obtained at Mach number $M_b = 10^{-4}$. Twenty equally reparted contours between $8 \times 10^{-6}$ and $3 \times 10^{-4}$ are plotted.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Euler cylinder scattering - Iso-contours of the norm of the velocity obtained at Mach number $M_b = 10^{-4}$. Twenty equally reparted contours between $8 \times 10^{-6}$ and $3 \times 10^{-4}$ are plotted.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7: Euler cylinder scattering - L² norm of the error between the exact incompressible solution and the numerical solution. Results are shown for a Mach number of $M_b = 10^{-4}$. A log-log plot is used.</figcaption></figure>

and $|\tilde{\rho}\tilde{u}_h|_{H^{\text{div}}}$, defined by

$$ \begin{aligned} |p_h|_{H^1}^2 &= \sum_{c \in \mathcal{C}} \int_c \|\nabla p_h\|^2 + \sum_{S \in \mathcal{S}_i} \int_S \llbracket p_h \rrbracket^2 \\ |\mathbf{u}_h|_{H^{\text{div}}}^2 &= \sum_{c \in \mathcal{C}} \int_c (\text{div}_{\mathbf{x}} \mathbf{u}_h)^2 + \sum_{S \in \mathcal{S}_i} \int_S \llbracket \mathbf{u}_h \cdot \mathbf{n}^S \rrbracket^2 \end{aligned} $$

as a function of the Mach number for a fixed mesh containing 10 × 20 cells. Using the new approximation space $d\mathbf{S}_h$ and Roe' numerical flux, we observe that $|\tilde{\rho}_h|_{H^1} = \mathcal{O}(M^2)$ and $|\tilde{\rho}\tilde{u}_h|_{H^{\text{div}}} = \mathcal{O}(M)$ and so

$$ \nabla \tilde{\rho} = \mathcal{O}(M^2) \quad \text{and} \quad \text{div}_{\mathbf{x}}(\tilde{\rho}\tilde{u}) = \mathcal{O}(M) $$

as expected.

### 5.3.2 Propagation of a low Mach number acoustic wave over a steady vortex

In this test case, we evaluate the ability of the numerical scheme to handle both incompressible phenomena and acoustic wave propagation at low Mach numbers. This test case was proposed in [11], because it had been remarked that some low Mach number fixes spoil the propagation of acoustic waves. The domain is the rectangle $[-0.1; 1.1] \times [0; 1]$. Periodic conditions are applied to top and bottom boundaries, and input-output conditions are used for left and right boundaries. The vortex is centered around $(x_c, y_c) = (0.5, 0.5)$, characterized by its reference Mach number $M_{\text{ref}}$, and is given by

$$ \left\{ \begin{array}{l} \rho(\mathbf{x}) = \rho_0 \left( 1 - \frac{M_{\text{ref}}^2}{\lambda_{\text{max}}^2} e^{-\frac{2\alpha^2}{1-\bar{r}^2}} \right), \\ \mathbf{u}(\mathbf{x}) = u_0 \left( \begin{array}{l} y \\ -x \end{array} \right) \frac{2\alpha}{\lambda_{\text{max}} r_0 (1 - \bar{r}^2)} e^{-\frac{\alpha^2}{1-\bar{r}^2}}, \end{array} \right. $$

where $\rho_0 = 2$, $u_0 = M_{\text{ref}}\sqrt{p'(\rho_0)}$, $\alpha = 2$, $r_0 = 0.45$, $\bar{r} = r/r_0$, $r = \sqrt{(x-x_c)^2 + (y-y_c)^2}$ and $\lambda_{\text{max}}$ is chosen such that the maximal velocity norm is $u_0$, i.e. $\lambda_{\text{max}} = \frac{2\alpha\bar{r}_{\text{max}}}{1-\bar{r}_{\text{max}}^2} e^{-\frac{\alpha^2}{1-\bar{r}_{\text{max}}^2}}$ with
<figure><img src="image_9.png" /><figcaption>Figure 7: Euler cylinder scattering - $L^2$ norm of the error between the exact incompressible solution and the numerical solution. Results are shown for a Mach number of $M_b = 10^{-4}$. A log-log plot is used.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Euler cylinder scattering - Semi-norms $|\tilde{\rho}|_{H^1}$ and $|\tilde{\rho}\tilde{u}|_{H^{\text{div}}}$ with respect to the Mach number for $M_b = 10^{-1}$ to $M_b = 10^{-7}$. A log-log plot is used.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: Euler cylinder scattering - Semi-norms $|\tilde{\rho}|_{H^1}$ and $|\tilde{\rho}\tilde{u}|_{H^{\text{div}}}$ with respect to the Mach number for $M_b = 10^{-1}$ to $M_b = 10^{-7}$. A log-log plot is used.</figcaption></figure>

$\bar{r}_{\text{max}} = \sqrt{-\alpha^2 + \sqrt{1 + \alpha^4}}$. The low Mach number acoustic wave is centered in $x = 0$ and given for $x \in [-0.05; 0.05]$ by

$$ \left\{ \begin{array}{l} \rho(\mathbf{x}) = \rho_0 \left(1 + M_{\text{ref}} e^{1 - \frac{1}{1-\alpha^2}}\right), \\ \mathbf{u}(\mathbf{x}) = u_0 + \frac{2}{\gamma-1} \left(\sqrt{p'(\rho(\mathbf{x}))} - \sqrt{p'(\rho_0)}\right), \end{array} \right. $$

where $\bar{x} = x/0.05$.

In Figure 9, the Mach number is plotted at different time with a mesh containing 480 × 400 cells for a reference Mach number $M_{\text{ref}} = 10^{-2}$. Using the new approximation space $d\mathbf{S}_h$ with Roe' numerical flux, the acoustic wave propagates correctly and the vortex is preserved over time. Using $d\mathbf{S}_h$ with Rusanov' numerical flux or $(d\mathbb{Q}_0)^2$, the acoustic wave propagates correctly but the vortex is quickly diffused over time.

# 6 Conclusion

In this article, the problem of low Mach number accuracy was addressed through the problem of discrete curl preservation. We proposed to use a special approximation space for the velocity, which is not the tensor product of the approximation space for scalars. This idea is partly inspired by the Raviart-Thomas or Nédéléc finite element basis, which are vector finite element spaces, but is different because all the degrees of freedom are located within the cells. With this new approximation space for velocities, we were able to design a curl, defined in the adjoint sense, that is preserved by the numerical scheme under mild assumption on the numerical flux. Going further, in the general quadrangular case, with boundary conditions, we were also able to prove the existence and uniqueness of a discrete Hodge-Helmholtz decomposition, and also to prove that this Hodge-Helmholtz decomposition is preserved by finite volume numerical scheme under the same assumption on the numerical flux. All these properties were confirmed with numerical experiments on the wave system.

The new approximation space for vectors was then used with the barotropic Euler system for the momentum variable, and was thoroughly tested in different low Mach number configurations. All numerical tests show that this single change of approximation space provides a strong benefit for solving the low Mach number accuracy on quadrangular meshes.

This paper deals with the case of dimension two. The new approximation space defined by (1) for Cartesian mesh and by (31) for general quadrangular mesh can be extended to dimension three. For dimension three, $\hat{\mathbf{S}}_0$ defined by (5) must be replaced by $\mathbb{Q}_0^3$ enriched by the two basis functions $(\hat{x} - 1/2, -\hat{y} + 1/2, 0)^T$ and $(\hat{x} - 1/2, 0, -\hat{z} + 1/2)^T$. All the results of the paper can be extended to dimension three.

Current investigations are focused on the higher order extension of this method, on extensions to the problem of conservation of divergence constraints, and also on the three dimensional extension of the scheme.
where $\rho_0 = 2$, $u_0 = M_{\text{ref}}\sqrt{p'(\rho_0)}$, $\alpha = 2$, $r_0 = 0.45$, $\bar{r} = r/r_0$, $r = \sqrt{(x-x_c)^2 + (y-y_c)^2}$ and $\lambda_{\text{max}}$ is chosen such that the maximal velocity norm is $u_0$, i.e. $\lambda_{\text{max}} = \frac{2\alpha\bar{r}_{\text{max}}}{1-\bar{r}_{\text{max}}^2}e^{-\frac{\alpha^2}{1-\bar{r}_{\text{max}}^2}}$ with $\bar{r}_{\text{max}} = \sqrt{-\alpha^2 + \sqrt{1+\alpha^4}}$. The low Mach number acoustic wave is centered in $x=0$ and given for $x \in [-0.05; 0.05]$ by

$$ \left\{ \begin{array}{l} \rho(\mathbf{x}) = \rho_0 \left(1 + M_{\text{ref}} e^{1-\frac{1}{1-\bar{x}^2}}\right), \\ \mathbf{u}(\mathbf{x}) = u_0 + \frac{2}{\gamma-1} \left(\sqrt{p'(\rho(\mathbf{x}))} - \sqrt{p'(\rho_0)}\right), \end{array} \right. $$

where $\bar{x} = x/0.05$.

In Figure 9, the Mach number is plotted at different time with a mesh containing 480 × 400 cells for a reference Mach number $M_{\text{ref}} = 10^{-2}$. Using the new approximation space $dS_h$ with Roe' numerical flux, the acoustic wave propagates correctly and the vortex is preserved over time. Using $dS_h$ with Rusanov' numerical flux or $(d\mathbb{Q}_0)^2$, the acoustic wave propagates correctly but the vortex is quickly diffused over time.

# 6 Conclusion

In this article, the problem of low Mach number accuracy was addressed through the problem of discrete curl preservation. We proposed to use a special approximation space for the velocity, which is not the tensor product of the approximation space for scalars. This idea is partly inspired by the Raviart-Thomas or Nédéléc finite element basis, which are vector finite element spaces, but is different because all the degrees of freedom are located within the cells. With this new approximation space for velocities, we were able to design a curl, defined in the adjoint sense, that is preserved by the numerical scheme under mild assumption on the numerical flux. Going further, in the general quadrangular case, with boundary conditions, we were also able to prove the existence and uniqueness of a discrete Hodge-Helmholtz decomposition, and also to prove that this Hodge-Helmholtz decomposition is preserved by finite volume numerical scheme under the same assumption on the numerical flux. All these properties were confirmed with numerical experiments on the wave system.

The new approximation space for vectors was then used with the barotropic Euler system for the momentum variable, and was thoroughly tested in different low Mach number configurations. All numerical tests show that this single change of approximation space provides a strong benefit for solving the low Mach number accuracy on quadrangular meshes.

This paper deals with the case of dimension two. The new approximation space defined by (1) for Cartesian mesh and by (30) for general quadrangular mesh can be extended to dimension three. For dimension three, $\hat{S}_0$ defined by (5) must be replaced by $\mathbb{Q}_0^3$ enriched by the two basis functions $(\hat{x} - 1/2, -\hat{y} + 1/2, 0)^T$ and $(\hat{x} - 1/2, 0, -\hat{z} + 1/2)^T$. All the results of the paper can be extended to dimension three.

Current investigations are focused on the higher order extension of this method, on extensions to the problem of conservation of divergence constraints, and also on the three dimensional extension of the scheme.

# References

[1] Douglas N. Arnold, Daniele Boffi, and Richard S. Falk. Quadrilateral H (div) finite elements. *SIAM Journal on Numerical Analysis*, 42(6):2429–2451, 2005.

[2] Dinshaw S Balsara, Roger Käppeli, Walter Boscheri, and Michael Dumbser. Curl constraint-preserving reconstruction and the guidance it gives for mimetic scheme design. *Communications on Applied Mathematics and Computation*, pages 1–60, 2021.

[3] Wasilij Barsukow. Stationarity preserving schemes for multi-dimensional linear systems. *Mathematics of Computation*, 88(318):1621–1645, 2019.

[4] Wasilij Barsukow. Truly multi-dimensional all-speed schemes for the euler equations on cartesian grids. *Journal of Computational Physics*, 435:110216, 2021.

[5] Wasilij Barsukow, Pierre-Henri Maire, and Raphaël Loubère. A node-conservative vorticity-preserving finite volume method for linear acoustics on unstructured grids. *Math. of Comp.*, 2023. Submitted.
<figure><img src="image_9.png" /><figcaption>Figure 9: Euler equations - Propagation of a low Mach acoustic wave over a steady vortex - Mach number obtained at different times for a reference Mach number $M_{\text{ref}} = 10^{-2}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Euler equations - Propagation of a low Mach acoustic wave over a steady vortex - Mach number obtained at different times for a reference Mach number $M_{\text{ref}} = 10^{-2}$.</figcaption></figure>
References

[1] Douglas N. Arnold, Daniele Boffi, and Richard S. Falk. Quadrilateral H (div) finite elements. *SIAM Journal on Numerical Analysis*, 42(6):2429–2451, 2005.

[2] Dinshaw S Balsara, Roger Käppeli, Walter Boscheri, and Michael Dumbser. Curl constraint-preserving reconstruction and the guidance it gives for mimetic scheme design. *Communications on Applied Mathematics and Computation*, pages 1–60, 2021.

[3] Wasilij Barsukow. Stationarity preserving schemes for multi-dimensional linear systems. *Mathematics of Computation*, 88(318):1621–1645, 2019.

[4] Wasilij Barsukow. Truly multi-dimensional all-speed schemes for the Euler equations on Cartesian grids. *Journal of Computational Physics*, 435:110216, 2021.

[5] Wasilij Barsukow, Pierre-Henri Maire, and Raphaël Loubère. A node-conservative vorticity-preserving finite volume method for linear acoustics on unstructured grids. *Math. of Comp.*, 2023. Submitted.

[6] John B. Bell, Phillip Colella, and Harland M. Glaz. A second-order projection method for the incompressible Navier-Stokes equations. *Journal of computational physics*, 85(2):257–283, 1989.

[7] Walter Boscheri, Giacomo Dimarco, and Lorenzo Pareschi. Locally structure-preserving div-curl operators for high order discontinuous Galerkin schemes. *Journal of Computational Physics*, 486:112130, 2023.

[8] Walter Boscheri, Michael Dumbser, Matteo Ioriatti, Ilya Peshkov, and Evgeniy Romenski. A structure-preserving staggered semi-implicit finite volume scheme for continuum mechanics. *Journal of Computational Physics*, 424:109866, 2021.

[9] Walter Boscheri, Raphaël Loubère, and Pierre-Henri Maire. An unconventional divergence preserving finite-volume discretization of Lagrangian ideal MHD. *Communications on Applied Mathematics and Computation*, pages 1–55, 2023.

[10] Jeremiah U. Brackbill and Daniel C. Barnes. The effect of nonzero $\nabla \cdot B$ on the numerical solution of the magnetohydrodynamic equations. *Journal of Computational Physics*, 35(3):426–430, 1980.

[11] Pascal Bruel, Simon Delmas, Jonathan Jung, and Vincent Perrier. A low Mach correction able to deal with low Mach acoustics. *Journal of Computational Physics*, 378:723–759, 2019.

[12] Zhiqiang Cai, Jim Douglas, and Xiu Ye. A stable nonconforming quadrilateral finite element method for the stationary Stokes and Navier-Stokes equations. *Calcolo*, 36:215–232, 1999.

[13] Michel Crouzeix and P-A Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes equations. *Revue française d'automatique informatique recherche opérationnelle. Mathématique*, 7(R3):33–75, 1973.

[14] Andreas Dedner, Friedemann Kemm, Dietmar Kröner, C-D Munz, Thomas Schnitzer, and Matthias Wesenberg. Hyperbolic divergence cleaning for the MHD equations. *Journal of Computational Physics*, 175(2):645–673, 2002.

[15] Sarah Delcourte, Komla Domelevo, and Pascal Omnes. A discrete duality finite volume approach to hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes. *SIAM Journal on Numerical Analysis*, 45(3):1142–1174, 2007.

[16] Stéphane Dellacherie, Pascal Omnes, and Felix Rieper. The influence of cell geometry on the Godunov scheme applied to the linear wave equation. *Journal of Computational Physics*, 229(14):5315–5338, 2010.

[17] Michael Dumbser, Francesco Fambri, Elena Gaburro, and Anne Reinarz. On GLM curl cleaning for a first order reduction of the CCZ4 formulation of the Einstein field equations. *Journal of Computational Physics*, 404:109088, 2020.

[18] Robert Eymard, Thierry Gallouët, Raphaele Herbin, and Jean-Claude Latché. Convergence of the MAC scheme for the compressible Stokes equations. *SIAM Journal on Numerical Analysis*, 48(6):2218–2246, 2010.