# Development of discontinuous Galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint

Vincent Perrier

► To cite this version:

Vincent Perrier. Development of discontinuous Galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint. 2024. hal-04564886

HAL Id: hal-04564886

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

Preprint submitted on 6 May 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Development of discontinuous Galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint: the case of linear systems

Vincent Perrier

► To cite this version:

Vincent Perrier. Development of discontinuous Galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint: the case of linear systems. Journal of Computational Physics, 2025, 544, pp.114445. 10.1016/j.jcp.2025.114445. hal-04564886v4

HAL Id: hal-04564886

https://inria.hal.science/hal-04564886v4

Submitted on 13 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Development of discontinuous Galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint

Vincent Perrier

Team Cagire, INRIA Bordeaux Sud-Ouest.

Laboratoire de Mathématiques et de leurs applications

Bâtiment IPRA, Université de Pau et des Pays de l'Adour,

Avenue de l'Université, 64013 Pau Cedex

May 6, 2024

## Abstract

Some hyperbolic systems are known to include implicit preservation of differential constraints: these are for example the time conservation of the curl or the divergence of a vector that appear as an implicit constraint. In this article, we show that this kind of constraint can be easily conserved at the discrete level with the classical discontinuous Galerkin method, provided the right approximation space is used for the vectorial space, and under some mild assumption on the numerical flux. For this, we develop a discrete differential geometry framework for some well chosen piece-wise polynomial vector approximation space. More precisely, we define the discrete Hodge star operator, the exterior derivative, and their adjoints. The discrete adjoint divergence and curl are proven to be exactly preserved by the discontinuous Galerkin method under a small assumption on the numerical flux. Numerical tests are performed on the wave system, the two dimensional Maxwell system and the induction equation, and confirm that the differential constraints are preserved at machine precision while keeping the high order of accuracy.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Continuous and discrete de-Rham complex</td><td>4</td></tr><tr><td>2.1</td><td>Differential forms and the de-Rham complex</td><td>4</td></tr><tr><td>2.1.1</td><td>Space of alternate forms</td><td>4</td></tr><tr><td>2.1.2</td><td>Space of differential forms and the de-Rham complex</td><td>6</td></tr><tr><td>2.2</td><td>The discrete case</td><td>7</td></tr><tr><td>2.2.1</td><td>The classical conformal discrete de-Rham complex</td><td>7</td></tr><tr><td>2.2.2</td><td>Some complexes involving a discrete discontinuous space for velocities</td><td>7</td></tr><tr><td>3</td><td>Discrete Hodge-star operator, codifferentials and Hodge Laplacian</td><td>10</td></tr><tr><td>3.1</td><td>Discrete Hodge star operator</td><td>10</td></tr><tr><td>3.2</td><td>Definition of the codifferential operators</td><td>11</td></tr><tr><td>3.2.1</td><td>Definition of δ<sup>1</sup></td><td>11</td></tr><tr><td>3.2.2</td><td>Definition of δ<sup>2</sup></td><td>11</td></tr><tr><td>3.3</td><td>Hodge Laplacian</td><td>12</td></tr></table>
# Development of discontinuous Galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint: the case of linear systems

Vincent Perrier

Team Cagire, INRIA Bordeaux Sud-Ouest.

Laboratoire de Mathématiques et de leurs applications

Bâtiment IPRA, Université de Pau et des Pays de l'Adour,

Avenue de l'Université, 64013 Pau Cedex

October 13, 2025

## Abstract

Some hyperbolic systems are known to include implicit preservation of differential constraints: these are for example the time conservation of the curl or the divergence of a vector that appear as an implicit constraint. In this article, we show that this kind of constraint can be easily conserved at the discrete level with the classical discontinuous Galerkin method, provided the right approximation space is used for the vectorial space, and under some mild assumption on the numerical flux. Even if the approximation space that is used does not always match with the one used classically for the discontinuous Galerkin method, the approximation space is *fully discontinuous*, i.e. the method is not a staggered method. We start by recalling a discrete de-Rham framework in which discontinuous approximation spaces for vectors fits. The discrete adjoint divergence and curl, which are *global* divergence and curl, are proven to be exactly preserved by the discontinuous Galerkin method under a small assumption on the numerical flux. Numerical tests are performed in dimension 2 on the wave system, the Maxwell system and the induction equation, on Cartesian meshes, and unstructured triangular and quadrangular meshes and confirm that the differential constraints are preserved at machine precision while keeping the high order of accuracy.

## Contents

1 Introduction 2

2 Approximation spaces 4

2.1 Continuous de-Rham complex 4

2.2 Conforming discrete de-Rham complex 5

2.3 Nonconforming discrete de-Rham complex 5

2.4 Scalar products and adjoint curl and divergence 7

2.5 Projection operator in the $C_k$ space 8

3 Discrete preservation of the curl or the divergence 9

3.1 Numerical scheme 9

3.2 Discrete preservation of the vorticity 9

3.3 Discrete preservation of the irrotational component 11

4 Preservation of initially divergence free field for the induction equation 12
<table><tr><td>4</td><td>Discrete preservation of the curl or the divergence</td><td>13</td></tr><tr><td>4.1</td><td>Numerical scheme</td><td>13</td></tr><tr><td>4.2</td><td>Discrete preservation of the vorticity</td><td>13</td></tr><tr><td>4.3</td><td>Discrete preservation of the irrotational component</td><td>15</td></tr></table>

<table><tr><td>5</td><td>Preservation of initially curl or divergence free fields for Lie derivative based advection equation</td><td>15</td></tr><tr><td>5.1</td><td>The direct and adjoint continuous equations</td><td>15</td></tr><tr><td>5.2</td><td>Discretization</td><td>17</td></tr></table>

<table><tr><td>6</td><td>Numerical results</td><td>18</td></tr><tr><td>6.1</td><td>Discrete conservation of the divergence: Maxwell system</td><td>20</td></tr><tr><td>6.1.1</td><td>Conservation of the divergence of a stationary solution</td><td>20</td></tr><tr><td>6.1.2</td><td>Convergence test</td><td>20</td></tr><tr><td>6.2</td><td>Discrete conservation of a curl: Wave system</td><td>26</td></tr><tr><td>6.2.1</td><td>Conservation of the curl of a stationary solution</td><td>26</td></tr><tr><td>6.2.2</td><td>Convergence test</td><td>27</td></tr><tr><td>6.3</td><td>Induction equation</td><td>32</td></tr><tr><td>6.3.1</td><td>Conservation of the divergence free field</td><td>33</td></tr><tr><td>6.3.2</td><td>Convergence test</td><td>33</td></tr></table>

<table><tr><td>7</td><td>Conclusion</td><td>35</td></tr></table>

<table><tr><td>A</td><td>Summary of the operators on differential forms and on their proxies</td><td>41</td></tr></table>

# 1 Introduction

In this article, we are interested in the discrete conservation of differential constraints that appear implicitly in a hyperbolic system of conservation law. Suppose for example that a hyperbolic system includes a vectorial unknown **u** for which the conservation law is

$$ \partial_t \mathbf{u} + \nabla \cdot \bar{\mathbf{G}} = 0, $$

where $\bar{\mathbf{G}}$ is a matrix. If $\bar{\mathbf{G}}$ is proportional to the identity matrix, namely if a scalar function $g$ exists such that $\bar{\mathbf{G}} = gI_d$, then the conservation law becomes

$$ \partial_t \mathbf{u} + \nabla g = 0, $$

and taking the curl of the equation on **u** gives formally the conservation of the curl of **u**:

$$ \partial_t (\nabla \times \mathbf{u}) = \nabla \times (\partial_t \mathbf{u}) = -\nabla \times (\nabla g) = 0. $$

In the same manner, if $\bar{\mathbf{G}}$ is antisymmetric, then the divergence of **u** is constant.

This kind of implicit differential constraint appears in a large number of systems including incompressible Navier-Stokes system, Maxwell system, magnetohydrodynamics (MHD), and wave and elastodynamic problems written in first order formulation. Preservation of such constraints at the discrete level has been addressed by several strategies, which can be gathered as

* **Staggering of unknowns.** Staggering unknowns consists in non collocated distribution of unknowns, for example in defining scalar unknowns in each cells, and vector unknowns in the sides of the mesh. The staggering of unknowns was proposed for the incompressible Navier-Stokes system in the MAC scheme [44], and for the Yee scheme [68] for Maxwell equations, see also [5, 6, 7] for the MHD system. Staggering of unknowns has been the object of a huge body of papers, including see e.g. [52, 51, 25, 28, 26] for the analysis of the MAC scheme. Among the large family of staggered schemes enters also the work on compatible
<table><tr><td><b>5</b></td><td><b>Numerical results</b></td><td><b>14</b></td></tr><tr><td>5.1</td><td>Discrete conservation of the divergence: Maxwell system</td><td>16</td></tr><tr><td>5.1.1</td><td>Conservation of the divergence of a stationary solution</td><td>16</td></tr><tr><td>5.1.2</td><td>Convergence test</td><td>18</td></tr><tr><td>5.2</td><td>Discrete conservation of a curl: Wave system</td><td>18</td></tr><tr><td>5.2.1</td><td>Conservation of the curl of a stationary solution</td><td>18</td></tr><tr><td>5.2.2</td><td>Convergence test</td><td>23</td></tr><tr><td>5.3</td><td>Induction equation</td><td>29</td></tr><tr><td>5.3.1</td><td>Conservation of the divergence free field</td><td>29</td></tr><tr><td>5.3.2</td><td>Convergence test</td><td>31</td></tr><tr><td>5.3.3</td><td>Stability test</td><td>31</td></tr><tr><td><b>6</b></td><td><b>Conclusion</b></td><td><b>34</b></td></tr></table>

# 1 Introduction

In this article, we are interested in the discrete conservation of differential constraints that appear implicitly in a hyperbolic system of conservation law. Suppose for example that a hyperbolic system includes a vectorial unknown **u** for which the conservation law is

$$ \partial_t \mathbf{u} + \nabla \cdot \bar{\mathbf{G}} = 0, $$

where $\bar{\mathbf{G}}$ is a matrix. If $\bar{\mathbf{G}}$ is proportional to the identity matrix, namely if a scalar function $g$ exists such that $\bar{\mathbf{G}} = gI_d$, then the conservation law becomes

$$ \partial_t \mathbf{u} + \nabla g = 0, $$

and taking the curl of the equation on **u** gives formally the conservation of the curl of **u**:

$$ \partial_t (\nabla \times \mathbf{u}) = \nabla \times (\partial_t \mathbf{u}) = -\nabla \times (\nabla g) = 0. $$

In the same manner, if $\bar{\mathbf{G}}$ is antisymmetric, then the divergence of **u** is constant.

This kind of implicit differential constraint appears in a large number of systems including incompressible Navier-Stokes system, Maxwell system, magnetohydrodynamics (MHD), and wave and elastodynamic problems written in first order formulation. Preservation of such constraints at the discrete level has been addressed by several strategies, which can be gathered as

* **Staggering of unknowns.** Staggering unknowns consists in non collocated distribution of unknowns, for example in defining scalar unknowns in each cells, and vector unknowns in the sides of the mesh. The staggering of unknowns was proposed for the incompressible Navier-Stokes system in the MAC scheme [54], and for the Yee scheme [86] for Maxwell equations, see also [5, 6, 8] for the MHD system. Staggering of unknowns has been the object of a huge body of papers, including see e.g. [68, 67, 31, 35, 32] for the analysis of the MAC scheme. Among the large family of staggered schemes enters also the work on compatible discretization [49, 60], or the staggered discontinuous Galerkin method [77]. This includes also the work on discrete exterior calculus with finite elements, including theoretical work around the Hodge Laplacian and its well posed mixed formulation [2, 3, 4], see also the book [1], work on electromagnetism [16, 17, 18, 47, 48] and discretization of the Lie-advection equation [44, 45, 46] and application [69]. Staggered discretizations were extended also to polytopal meshes, mainly for incompressible systems see e.g. the virtual element method [11], the Hybrid High order (HHO) methods [26] or the Compatible Discrete Operators (CDO) method [13, 14, 63]. The main difficulty of staggered discretizations for nonlinear hyperbolic systems is to keep a correct definition of the conservation; especially, keeping conservation while shock-limiting staggered data is challenging.
discretization [41, 46], or the staggered discontinuous Galerkin method [61]. This includes also the work on discrete exterior calculus with finite elements, including theoretical work around the Hodge Laplacian and its well posed mixed formulation [2, 3, 4], see also the book [1], work on electromagnetism [15, 16, 17, 39, 40] and discretization of the Lie-advection equation [36, 37, 38] and application [53]. Staggered discretizations were extended also to polytopal meshes, mainly for incompressible systems see e.g. the virtual element method [10], the Hybrid High order (HHO) methods [22] or the Compatible Discrete Operators (CDO) method [12, 13, 47]. The main difficulty of staggered discretizations for nonlinear hyperbolic systems is to keep a correct definition of the conservation; especially, keeping conservation while shock-limiting staggered data is challenging.

* **Projection method.** The projection method was mainly used for divergence cleaning in [18] for the MHD system and in [11] for incompressible Navier-Stokes system. The projection method is a predictor-corrector method. Suppose for example that the divergence of a vector $v$ should be preserved. Suppose also that a predictor step provided a candidate update for $v^{n+1}$ of $v^n$, but that the divergence is not preserved. Then a potential $\varphi$ such that

$$ \Delta\varphi = \nabla \cdot (v^{n+1} - v^n), $$

may be computed. Then $\tilde{v}^{n+1} := v^{n+1} - \nabla\varphi$ is a projection of $v^{n+1}$ for which the divergence is preserved. The main drawback of this method is the cost of the inversion of an elliptic system at each time step, and the definition of the boundary conditions for this elliptic system.

* **Generalized Lagrange Multiplier method.** This method, first developed in [49, 20] consists in considering the divergence to be preserved as an additional variable. Then an additional equation for this variable, and a relaxation process ensure that the divergence is asymptotically preserved. This method was extended to curl preservation in [23]. The difficulty with this system is its higher computational cost because of the additional variables, and the tuning of the numerical relaxation parameters.

Another category of method was especially designed for the conservation of the zero divergence of the magnetic field in the Maxwell system, the MHD system or the induction model system; these systems are slightly different from the ones previously discussed, because the divergence of the magnetic field is not directly preserved by the system. Instead, the divergence is solution of a conservative transport equation, so that an initially divergence free magnetic field is divergence free for all time. These two categories of methods are

* **the Godunov-Powell method.** This method consists in relying on the formulation of Godunov [30, 29] and Powell [55, 56] for the MHD equation, and in trying to control the divergence of the magnetic field [27, 48] (note however that the magnetic field is not divergence free, the aim of the method is only to keep this divergence "low").

* **Constrained Transport Method.** This method was originally proposed in [24] as an alternative of the Yee scheme [68], and under a staggered fashion. This method is based on staggered ideas, and in the Maxwell context, the challenge consists in the computation of a reliable staggered electric field based on a magnetic field that, still may not be staggered: indeed, in [65], several methods were compared, including different versions of the constrained transport method, and it was shown that the method can be collocated, see also [64, 63, 35]. Concerning the constrained transport method, we also refer to the review [62] and references therein, [66] for a high order constrained transport method based on spectral differences, and [42] for an extension to the preservation of a curl.

Apart from these schemes that are especially designed for preserving exactly a discrete curl or a discrete divergence, some collocated numerical schemes seem to be naturally able to preserve these constraints. This is for example the case of numerical schemes developed within the low Mach
* **Projection method.** The projection method was mainly used for divergence cleaning in [19] for the MHD system and in [12] for incompressible Navier-Stokes system, see also [23, the globally divergence free method]. The projection method is a predictor-corrector method. Suppose for example that the divergence of a vector **v** should be preserved. Suppose also that a predictor step provided a candidate update for **v**<sup>n+1</sup> of **v**<sup>n</sup>, but that the divergence is not preserved. Then a potential ϕ such that

$$ \Delta\varphi = \nabla \cdot (\mathbf{v}^{n+1} - \mathbf{v}^n), $$

may be computed. Then $\tilde{\mathbf{v}}^{n+1} := \mathbf{v}^{n+1} - \nabla\varphi$ is a projection of $\mathbf{v}^{n+1}$ for which the divergence is preserved. The main drawback of this method is the cost of the inversion of an elliptic system at each time step, and the definition of the boundary conditions for this elliptic system.

* **Generalized Lagrange Multiplier method.** This method, first developed in [65, 24] consists in considering the divergence to be preserved as an additional variable. Then an additional equation for this variable, and a relaxation process ensure that the divergence is asymptotically preserved. This method was extended to curl preservation in [29]. The difficulty with this system is its higher computational cost because of the additional variables, and the tuning of the numerical relaxation parameters.

Another category of method was especially designed for the conservation of the zero divergence of the magnetic field in the Maxwell system, the MHD system or the induction model system; these systems are slightly different from the ones previously discussed, because the divergence of the magnetic field is not directly preserved by the system. Instead, the divergence is solution of a conservative transport equation, so that an initially divergence free magnetic field is divergence free for all time. These two categories of methods are

* **the Godunov-Powell method.** This method consists in relying on the formulation of Godunov [37, 36] and Powell [73, 74] for the MHD equation, and in trying to control the divergence of the magnetic field [34, 64] (note however that the magnetic field is not divergence free, the aim of the method is only to keep this divergence "low").

* **Constrained Transport Method.** This method was originally proposed in [30] as an alternative of the Yee scheme [86], and under a staggered fashion see also [58] for an application to the high order discretization of the MHD system. This method is based on staggered ideas, and in the Maxwell context, the challenge consists in the computation of a reliable staggered electric field based on a magnetic field that, still may not be staggered: indeed, in [81], several methods were compared, including different versions of the constrained transport method, and it was shown that the method can be collocated, see also [80, 79, 43]. Concerning the constrained transport method, we also refer to the review [78] and references therein, [82] for a high order constrained transport method based on spectral differences, and [50] for an extension to the preservation of a curl.

Apart from these schemes that are especially designed for preserving exactly a discrete curl or a discrete divergence, some collocated numerical schemes seem to be naturally able to preserve these constraints. This is for example the case of numerical schemes developed within the low Mach number community [25, 39, 40, 52], which inspired the present article, and which was extended for low order quadrangular meshes in [53]. Even if these references do no state any relation with curl preservation, they are actually preserving a curl that will be defined in this article. Another family of schemes that seem to preserve exactly divergence or curl constraints are node based numerical schemes, see e.g. [15] for the divergence preservation or [9, 10] for curl preservation. Considering their stencils, these schemes seem to have a close relationship with [79, 50], which are based on the unstaggered constrained transport method.

In this article, we wish to develop two-dimensional discontinuous Galerkin methods that naturally preserve a global divergence or a curl exactly, which match with [25, 39, 40, 52] for straight
number community [21, 32, 33, 43], which inspired the present article. Even if these references do not state any relation with curl preservation, they are actually preserving a curl that will be defined in this article. Another family of schemes that seem to preserve exactly divergence or curl constraints are node based numerical schemes, see e.g. [14] for the divergence preservation or [8, 9] for curl preservation. Considering their stencils, these schemes seem to have a close relationship with [63, 42], which are based on the unstaggered constrained transport method. In this article, we wish to develop discontinuous Galerkin methods that naturally preserve a divergence or a curl exactly, which match with [21, 32, 33, 43] for straight triangular meshes but which differ on quadrangular meshes because we use alternative approximation space for vectors, proposed in [54].

This article is organized as follows. In section 2, the basics of the two-dimensional de-Rham complex is recalled. The classical conformal finite element approximation based on the continuous / Nédélec / Raviart-Thomas complex is recalled, and the nonconforming approach proposed in [54] is reviewed. The classical discrete operators (Hodge-star operator, codifferentials) are then defined in section 3. Based on these approximation spaces and matching discrete operators, we prove in section 4 that the classical discontinuous Galerkin method is able to preserve a curl or a divergence constraint provided the right approximation space is used, and under a small hypothesis regarding the diffusion direction of the numerical flux. Then in section 5, we explain how the numerical scheme can be also extended to the induction equation; we especially prove that if the vector field is correctly initialized, then it is divergence free (still in the adjoint sense) for all time. In section 6, the numerical scheme is tested for the preservation of a divergence with the two dimensional Maxwell system, the preservation of a curl with the two dimensional first order formulation of the wave system, and the preservation of the divergence free field with the induction equation. In this numerical section, convergence tests are also performed. This article finishes with the conclusion in section 7.

## 2 Continuous and discrete de-Rham complex

### 2.1 Differential forms and the de-Rham complex

#### 2.1.1 Space of alternate forms

Differential forms are a powerful tool which may be used also in the partial differential equations context for exhibiting and analyzing geometrical structure of these equations.

This theory starts with alternate k-multilinear forms on $\mathbb{R}^2$ that we denote by $\Lambda^k$. In this paper, continuous and discrete linear forms will be denoted by gothic letters, e.g. $f, u$ or $g$. Scalars will be denoted by italic small letters (e.g. $f, g$), and vectors by roman bold letters (e.g. $u, v$). On $\mathbb{R}^2$, $\Lambda^k$ is the null space for $k \ge 3$. For $0 \le k \le 2$, we have:

* $\Lambda^0$ is the set of 0-linear forms on $\mathbb{R}^2$; a 0 multilinear form $f \in \Lambda^0$ is completely determined by a constant $f$:

$$f() = f.$$

* $\Lambda^1$ is the set of 1-linear forms. They are usually represented in two manner

- Either by the scalar product by a vector $u$:

$$\forall u \in \Lambda^1 \quad \exists! u \in \mathbb{R}^2 \quad \forall v_1 \in \mathbb{R}^2 \quad u(v_1) = u \cdot v_1, \quad (1)$$

- or by the determinant with a vector $\tilde{u}$:

$$\forall u \in \Lambda^1 \quad \exists! \tilde{u} \in \mathbb{R}^2 \quad \forall v_1 \in \mathbb{R}^2 \quad u(v_1) = \det(v_1, \tilde{u}). \quad (2)$$

Note that $u \cdot v_1 = \det(v_1, \tilde{u})$ if and only if $\tilde{u}$ is the image of $u$ by the $\pi/2$ rotation which we denote by a $\perp$ exponent: $\tilde{u} = u^\perp$.
triangular meshes but which differ on quadrangular meshes because we use alternative approximation space for vectors, proposed in [71].

This article is organized as follows. In section 2, the approximation spaces developed in [71] are recalled, and their relation with the de-Rham complex is recalled. Based on these approximation spaces and matching discrete operators, we prove in section 3 that the classical discontinuous Galerkin method is able to preserve a curl or a divergence constraint provided the right approximation space is used, and under a small hypothesis regarding the diffusion direction of the numerical flux. Then in section 4, we explain how the numerical scheme can be also extended to the induction equation; we especially prove that if the vector field is correctly initialized, then it is divergence free (still in the adjoint sense) for all time. In section 5, the numerical scheme is tested for the preservation of a divergence with the two dimensional Maxwell system, the preservation of a curl with the two dimensional first order formulation of the wave system, and the preservation of the divergence free field with the induction equation. In this numerical section, convergence tests are also performed. This article finishes with the conclusion in section 6.

# 2 Approximation spaces

## 2.1 Continuous de-Rham complex

The most general context for introducing the gradient, curl and divergence operator relies on exterior calculus. In dimension 2, this consists in considering the de-Rham complex:

$$ \Lambda^0(\Omega) \xrightarrow{d^0} \Lambda^1(\Omega) \xrightarrow{d^1} \Lambda^2(\Omega), \qquad (1) $$

where $\Lambda^k(\Omega)$ is the set of smooth $k$-differential forms on $\Omega$, and $d$ is the exterior derivative (see [21, Chap. 1.4]). When dealing with partial differential equations, it is usually more convenient to describe (1) in terms of *proxies*. In dimension 2, a differential form of $\Lambda^1$ is usually represented in two manners:

* Either by the scalar product by a vector field $\mathbf{u}$:

$$ \forall \mathbf{u} \in \Lambda^1 \quad \exists! \mathbf{u} \in \mathbb{R}^2 \quad \forall \mathbf{v}_1 \in \mathbb{R}^2 \quad \mathbf{u}(\mathbf{v}_1) = \mathbf{u} \cdot \mathbf{v}_1, \qquad (2) $$

* or by the determinant with a vector field $\tilde{\mathbf{u}}$:

$$ \forall \mathbf{u} \in \Lambda^1 \quad \exists! \tilde{\mathbf{u}} \in \mathbb{R}^2 \quad \forall \mathbf{v}_1 \in \mathbb{R}^2 \quad \mathbf{u}(\mathbf{v}_1) = \det(\mathbf{v}_1, \tilde{\mathbf{u}}). \qquad (3) $$

Note that $\mathbf{u} \cdot \mathbf{v}_1 = \det(\mathbf{v}_1, \tilde{\mathbf{u}})$ if and only if $\tilde{\mathbf{u}}$ is the image of $\mathbf{u}$ by the $\pi/2$ rotation which we denote by a $\perp$ exponent: $\tilde{\mathbf{u}} = \mathbf{u}^{\perp}$.

From now on, we denote by $\nabla$ the vector operator $(\partial_x, \partial_y)$ and $\nabla^{\perp}$ the vector operator $(-\partial_y, \partial_x)$. Then, depending on the choice of proxy, (1) becomes, if (2) is chosen:

$$ C^{\infty}(\Omega) \xrightarrow{\nabla} (C^{\infty}(\Omega))^2 \xrightarrow{\nabla^{\perp}} C^{\infty}(\Omega), \qquad (4) $$

whereas if (3) is chosen, (1) becomes:

$$ C^{\infty}(\Omega) \xrightarrow{\nabla^{\perp}} (C^{\infty}(\Omega))^2 \xrightarrow{-\nabla} C^{\infty}(\Omega). \qquad (5) $$

A fundamental result on the de-Rham complex, is the link between the dimension of the cohomology groups of (1) and the Betti numbers [42], which are characteristic of the topology of the domain:

$$ \left\{ \begin{array}{l} b_0 = \dim (\ker \nabla), \\ b_1 = \dim (\ker (\nabla^{\perp}.) / \operatorname{Range}(\nabla)), \\ b_2 = \dim (C^{\infty}(\Omega) / \operatorname{Range}(\nabla^{\perp}.)). \end{array} \right. \qquad (6) $$

In this article, we focus on the two-dimensional torus $\mathbb{T}^2$, for which the Betti numbers are $b_0 = b_2 = 1$ and $b_1 = 2$.
<table><caption>Table 1: Expression of the effect of the exterior product on the scalar and vector proxies. This effect is independent of the choice of the proxy for $\Lambda^2$ ((1) or (2)).</caption><thead><tr><th>u</th><th>v</th><th>Proxy (u &wedge; v)</th></tr></thead><tbody><tr><td>&Lambda;<sup>0</sup></td><td>&Lambda;<sup>0</sup></td><td>uv</td></tr><tr><td>&Lambda;<sup>0</sup></td><td>&Lambda;<sup>1</sup></td><td>uv</td></tr><tr><td>&Lambda;<sup>0</sup></td><td>&Lambda;<sup>2</sup></td><td>uv</td></tr><tr><td>&Lambda;<sup>1</sup></td><td>&Lambda;<sup>1</sup></td><td>det (u, v)</td></tr></tbody></table>

* $\Lambda^2$ is the set of alternate 2-linear forms. It is known that all the alternate 2-linear forms are proportional to the determinant:

$$ \forall f \in \Lambda^2 \quad \exists! f \in \mathbb{R} \quad \forall \mathbf{v}_1, \mathbf{v}_2 \in \mathbb{R}^2 \quad f(\mathbf{v}_1, \mathbf{v}_2) = f \det(\mathbf{v}_1, \mathbf{v}_2). $$

The choice of the proxy (1) or (2) will have an impact on the whole list of operators that will be defined. We will denote by Proxy the proxy of a linear form when (1) is chosen, and by Proxy$^\perp$ the proxy of a linear form when (2) is chosen. A large number of operators on linear forms will be introduced all along this article, and the matches between operators on proxies and operators on linear forms is summarized in Appendix A. The scalar $f$ is called the scalar proxy of the alternate 0 or 2 linear form $f$, whereas the vectors $u$ or $\tilde{u}$ are called the vector proxies of the alternate 2-linear form $u$. Several additional tools are needed on these sets of alternate linear forms:

* The exterior product $\wedge$ is defined in general as [19, Chap. 1.4]

$$ \begin{array}{ccc} \wedge : & \Lambda^i \times \Lambda^j & \mapsto \Lambda^{i+j} \\ (u, v) & \mapsto u \wedge v, & \end{array} $$

with

$$ u \wedge v (\mathbf{v}_1, \dots, \mathbf{v}_i, \mathbf{v}_{i+1}, \dots, \mathbf{v}_{i+j}) = \sum_{\sigma \in \mathcal{S}_{i+j}} \varepsilon(\sigma) u (\mathbf{v}_{\sigma(1)}, \dots, \mathbf{v}_{\sigma(i)}) v (\mathbf{v}_{\sigma(i+1)}, \dots, \mathbf{v}_{\sigma(i+j)}), $$

where $\mathcal{S}_{i+j}$ is the set of $(i+j)$-permutations such that

$$ \sigma(1) < \dots < \sigma(i) \quad \text{and} \quad \sigma(i+1) < \dots < \sigma(i+j), $$

and $\varepsilon(\sigma)$ is the signature of the permutation $\sigma$. The exterior product ensures $u \wedge v = (-1)^{ij} v \wedge u$. The effect of the exterior product on the proxies is recalled in Table 1.

* A scalar product may be defined on $\Lambda^k$, which in our case is equivalent to the classical scalar or vector scalar product on the proxies, and which we denote by $\langle \cdot | \cdot \rangle_{\Lambda^k}$

* Last, we denote by vol the volume form, the element of $\Lambda^2$ such that

$$ \forall \mathbf{v}_1, \mathbf{v}_2 \in \mathbb{R}^2 \quad \mathrm{vol}(\mathbf{v}_1, \mathbf{v}_2) := \det(\mathbf{v}_1, \mathbf{v}_2), $$

and we define the Hodge-star operator, that will be denoted by $\star_{\{k\}}$, as the operator from $\Lambda^k$ to $\Lambda^{2-k}$ such that

$$ \forall u \in \Lambda^k \quad \forall v \in \Lambda^{2-k} \quad u \wedge v = \langle \star_{\{k\}} (u) | v \rangle_{\Lambda^{2-k}} \mathrm{vol}. $$

Concerning the proxies, $\star_{\{0\}}$ and $\star_{\{2\}}$ are the identity, whereas $\star_{\{1\}}$ is the $\pi/2$ rotation.
## 2.2 Conforming discrete de-Rham complex

From an approximation point of view, the study of discrete counterpart of the complex of (1) has been an intensive research topic over the last forty years, including the work of Whitney [83], Bossavit and Hiptmair on electromagnetism [16, 17, 18, 48], and the work on formalization of Arnold and collaborators of finite element exterior calculus [2, 3] which led to the reference book [1]. A classical conformal discrete counterpart of the de-Rham diagram of (1) on triangles relies on the family of $P_r^- \Lambda^k$ of [1, Chapter 7], for which a key property is preservation of the equalities (6) at the discrete level. The family $P_r^- \Lambda^k$ can then be translated back into continuous/Nédélec/Raviart-Thomas complex [75, 76, 66] for the proxies. For example, the discrete counterpart of (4) is:

$$ P_{k+1} \xrightarrow{\nabla} N_k \xrightarrow{\nabla^{\perp}} dP_k, \qquad (7) $$

where $P_{k+1}$ is the continuous finite elements space of degree $k+1$, $N_k$ is two-dimensional Nédélec approximation space (vectorial finite element space that is tangential conforming), and $dP_k$ is the discontinuous Galerkin approximation space. The discrete counterpart of (5) is:

$$ P_{k+1} \xrightarrow{\nabla^{\perp}} RT_k \xrightarrow{-\nabla} dP_k, \qquad (8) $$

where $RT_k$ is the two-dimensional Raviart-Thomas approximation space (vectorial finite element space that is normal conforming).

## 2.3 Nonconforming discrete de-Rham complex

We begin by providing some definitions on the mesh we are working on. We work on the two-dimensional torus $T^2$, which is supposed to be oriented. This domain is divided into a set of cells which are either triangles or quadrangles. The set of sides of the mesh is denoted by $S$. Each side $S$ of the mesh is oriented once for all, and we denote by $n_S$ the clockwise normal to the side. Each side $S$ has two neighbouring cells, the one for which $n_S$ is outgoing is the left cell, and the one for which $n_S$ is ingoing is the right cell. The jump of a scalar $f$ that is continuous on the cells and discontinuous across the faces is defined as

$$ \forall S \in S \quad [f]_S := f_L - f_R, $$

where $f_L$ is the value on the left cell and $f_R$ is the value on the right cell.

Recently, it was proposed to relax the continuity constraints induced by the conformal hypothesis [59]. Contrarily to the classical conformal approximation, the approximation of the vector spaces proposed in [59] are completely discontinuous, without any hypothesis on the continuity of the normal component (which holds for $RT$) or the tangential component (which holds for $N$). This idea was extended in [71] for finding basis with fewer degrees of freedom. The different approximation space for vectors proposed in [71] will be denoted by $dB_k^{\text{div}}$ and $dB_k^{\text{curl}}$. In this case, the classical derivation operators cannot be applied to these approximation space, and the derivative in the sense of distributions is used instead

$$ \forall \mathbf{u} \in dB_k^{\text{div}} \qquad \nabla_{\mathcal{D}'} \cdot \mathbf{u} := \begin{cases} \forall C \in \mathcal{C} & (\nabla_{\mathcal{D}'} \cdot \mathbf{u})^C := \nabla \cdot \mathbf{u} \\ \forall S \in \mathcal{S} & (\nabla_{\mathcal{D}'} \cdot \mathbf{u})^S := -[\mathbf{u}] \cdot \mathbf{n}_S. \end{cases} $$

and

$$ \forall \mathbf{u} \in dB_k^{\text{curl}} \qquad \nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{u} := \begin{cases} \forall C \in \mathcal{C} & (\nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{u})^C := \nabla^{\perp} \cdot \mathbf{u} \\ \forall S \in \mathcal{S} & (\nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{u})^S := [\mathbf{u}^{\perp}] \cdot \mathbf{n}_S. \end{cases} $$

We see that the image of the operators $(\nabla_{\mathcal{D}'} \cdot)$ and $(\nabla_{\mathcal{D}'}^{\perp} \cdot)$ is in a Cartesian product of cell space and face space, that we will denote in general $C_k$. In [71], several such vectorial approximation spaces were proposed, for which we proved that the discrete and continuous cohomology spaces are matching on $T^2$. The proposed approximation spaces may then be put in the following complex:
### 2.1.2 Space of differential forms and the de-Rham complex

We denote by $C^\infty\Lambda^k$ the space of differential $k$-forms, namely the functions $T^2 \mapsto \Lambda^k$ that are infinitely differentiable, where $T^2$ is the two-dimensional torus. It is equipped with the norm

$$ \langle \mathbf{u} | \mathbf{v} \rangle_{L^2 \Lambda^k} := \int_{T^2} \langle \mathbf{u} | \mathbf{v} \rangle_{\Lambda^k} . $$

If $\mathbf{u} \in C^\infty\Lambda^k$, then its derivative $\mathbf{u}'$ is a linear application from $\mathbb{R}^2$ to $\Lambda^k$, which may be seen as a $(k+1)$-linear application. Denoting by

$$ \partial_{\mathbf{v}_j} \mathbf{u} := \mathbf{u}' \cdot \mathbf{v}_j, $$

the exterior derivative [19, Chap. 2.3] $d^k$ is defined as

$$ \begin{array}{ccc} d^k : C^\infty\Lambda^k & \mapsto & C^\infty\Lambda^{k+1} \\ \mathbf{u} & \mapsto & d^k \mathbf{u}, \end{array} $$

with

$$ (d^k \mathbf{u}) (\mathbf{v}_1, \dots, \mathbf{v}_{k+1}) = \sum_{j=1}^{k+1} (-1)^{j-1} \partial_{\mathbf{v}_j} (\mathbf{v}_1, \dots, \mathbf{v}_{j-1}, \mathbf{v}_{j+1}, \dots, \mathbf{v}_k). $$

Denoting by $\nabla$ the vector operator $(\partial_x, \partial_y)^T$, then

* If the choice of proxy is (1), then $d^0$ is $\nabla$, whereas $d^1$ is $\nabla^\perp$.

* If the choice of proxy is (2), then $d^1$ is $\nabla^\perp$, whereas $d^1$ is $-\nabla^\perp$.

Regarding the relations between the $\pi/2$ rotation and the different differential operators, we wish to recall the following relations that will be useful in section 3:

$$ \left\{ \begin{array}{l} \nabla^\perp \cdot \mathbf{u} = -\nabla \cdot \mathbf{u}^\perp \\ \nabla \cdot (\mathbf{u}^\perp) = -\nabla^\perp \cdot \mathbf{u}, \\ \mathbf{u} \cdot \nabla^\perp f = -\mathbf{u}^\perp \cdot \nabla f. \end{array} \right. $$

Last, we denote by $\delta^k$ the coderivative which are defined in two dimensions as

$$ \delta^1 := -\star_{\{2\}} d^1 \star_{\{1\}}, \qquad (3) $$

and

$$ \delta^2 := -\star_{\{1\}} d^0 \star_{\{2\}}. \qquad (4) $$

Concerning the proxies, the $\delta$ operators are

* If the choice of proxy is (1), then $\delta^2$ is $(-\nabla^\perp)$, whereas $\delta^1$ is $(-\nabla)$.

* If the choice of proxy is (2), then $\delta^2$ is $\nabla$, whereas $\delta^1$ is $(-\nabla^\perp)$.

$C^\infty\Lambda^k$ may be equipped with the graph norm $\| \cdot \|_k + \|d^k \cdot \|_{k+1}$, and the closure of $C^\infty\Lambda^k$ is denoted by $HA^k$. This gives rise to a complex similar to the de-Rham complex, but with Sobolev spaces $HA^k$, called Hilbert complex [1, Chapter 4].

In Figure 1, the two-dimensional Hilbert complex was written, both in terms of differential forms spaces $HA^k$, and in terms of their vector proxies. A fundamental result on the de-Rham complex, which holds also for the Hilbert complex, is the link between the dimension of the cohomology groups and the Betti numbers[34], which are characteristic of the topology of the domain:

$$ \left\{ \begin{array}{l} b_0 = \dim (\ker \nabla), \\ b_1 = \dim (\ker (\nabla^\perp \cdot) / \operatorname{Range}(\nabla)), \\ b_2 = \dim (L^2 / \operatorname{Range}(\nabla^\perp \cdot)). \end{array} \right. \qquad (5) $$

In our case, on $T^2$, the Betti numbers are $b_0 = b_2 = 1$ and $b_1 = 2$.
• If the choice of proxy is (2),

$$ A_{k+1} \xrightarrow{\nabla} \mathbf{d}\mathcal{B}_k^{\text{curl}} \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} C_k, \qquad (9) $$

where $A_{k+1}$ is the space of continuous finite element spaces, i.e. either $\mathbb{P}_{k+1}$ on triangular meshes, or $\mathbb{Q}_{k+1}$ on quadrangular meshes, and $\mathbf{d}\mathcal{B}_k^{\text{curl}}$ is discussed in the remaining part of the section.

• If the choice of proxy is (3), then

$$ A_{k+1} \xrightarrow{\nabla^{\perp}} \mathbf{d}\mathcal{B}_k^{\text{div}} \xrightarrow{-\nabla_{\mathcal{D}'}^{\perp}} C_k, \qquad (10) $$

where $\mathbf{d}\mathcal{B}_k^{\text{div}}$ is discussed in the remaining part of the section.

Note that in (10) and (9), the spaces $C_k$ and $A_{k+1}$ are the same, only $\mathbf{d}\mathcal{B}_k^{\text{div}}$ and $\mathbf{d}\mathcal{B}_k^{\text{curl}}$ are different. In [71], the following approximation spaces were proposed; the first family is based on the conformal case, whereas the two other proposed approximation spaces are optimal in the number of degrees of freedom.

1. **Approximation spaces directly based on discontinuous versions of Raviart-Thomas/Nédélec finite elements.** In this version of the discontinuous spaces, proposed in [59], the polynomial basis are exactly the same as for the conformal approximation on each cell, but all the continuity constraints are relaxed. These spaces were denoted by $\mathbf{d}\mathbb{R}^k$ and $\mathbf{d}\mathbb{N}_k$. This leads to

$$ C_k = \mathrm{d}\mathbb{P}_k(\mathcal{C}) \times \mathrm{d}\mathbb{P}_k(\mathcal{S}), $$

for triangular meshes and to

$$ C_k = \mathrm{d}\mathbb{Q}_k(\mathcal{C}) \times \mathrm{d}\mathbb{P}_k(\mathcal{S}), $$

for quadrangular meshes. These approximation space work fine, but some degrees of freedom are useless for ensuring the matching between the discrete and continuous cohomology groups. Also the Raviart-Thomas/Nédélec finite element basis should be generated which is not always straightforward on triangular meshes.

2. **Optimal approximation spaces on triangles.** On triangles, the classical approximation space for vectors, namely the one obtained by tensorizing the classical approximation space for scalars can be put in the diagrams (9) and (10), as addressed in [71]; in this case, $\mathbf{d}\mathcal{B}_k^{\text{curl}} = \mathbf{d}\mathcal{B}_k^{\text{div}} = \mathbf{d}\mathbb{P}_k$, and

$$ C_k = \mathrm{d}\mathbb{P}_{k-1}(\mathcal{C}) \times \mathrm{d}\mathbb{P}_k(\mathcal{S}). $$

3. **Optimal approximation spaces on quadrangles.** The case of quadrangular meshes is quite surprising because the lowest order approximation space (namely with $A_{k+1} = \mathbb{Q}_1$) does not give the classical finite volume vector space, but rather an enriched version including three basis vector instead of two:

$$ \mathbf{d}\mathcal{B}_0^{\text{div}} = \mathrm{Span} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} -x \\ y \end{pmatrix} \right), $$

and

$$ \mathbf{d}\mathcal{B}_0^{\text{curl}} = \mathrm{Span} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} y \\ x \end{pmatrix} \right). $$

In the general case, following [71] we denote by $\mathrm{d}\mathbb{Q}_{i,j}$ the space of polynomials of degree lower than $i$ in $x$ and lower than $j$ in $y$, and we set

$$ \mathbf{d}\mathcal{B}_k^{\text{div}} = \left( \begin{array}{c} \mathrm{d}\mathbb{Q}_{k,k} + \mathrm{d}\mathbb{Q}_{k+1,k-1} \\ \mathrm{d}\mathbb{Q}_{k,k} + \mathrm{d}\mathbb{Q}_{k-1,k+1} \end{array} \right) \oplus \mathrm{Span} \left( \begin{array}{c} -x^{k+1}y^k \\ x^k y^{k+1} \end{array} \right) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Two dimensional Hilbert complex. On the top, the differential forms version, with the exterior differential d and the codifferential δ. On the bottom part, the vector proxy version, with the different choices of proxy: (1), leading to the complex ∇ / ( ∇ ⊥ · ) and (2), leading to the complex ∇ ⊥ / ( − ∇ · ) versions, and their adjoints. The Figure can be also read column-wise, for example the space of differential forms H Λ 1 matches with H curl or H div depending on the choice of proxy representation of the linear forms, and the exterior differential d 0 matches with ∇ or ∇ ⊥ depending on the choice of proxy.</figcaption></figure>

2.2 The discrete case

2.2.1 The classical conformal discrete de-Rham complex

From an approximation point of view, the study of discrete counterpart of the Hilbert complex of Figure 1 has been an intensive research topic over the last forty years, including the work of Whitney [67], Bossavit and Hiptmair on electromagnetism [15, 16, 17, 40], and the work on formalization of Arnold and collaborators of finite element exterior calculus [2, 3] which led to the reference book [1]. The classical conformal discrete counterpart of the de-Rham diagram of Figure 1 relies on the continuous/Nédélec/Raviart-Thomas complex [4, 57, 58, 50], where H 1 is replaced by the continuous finite elements P$_{k}$$_{+1}$ , H div by the two-dimensional Raviart-Thomas approximation space RT$_{k}$ , H curl by the two-dimensional Nédélec approximation space N$_{k}$ , and L 2 by the discontinuous Galerkin approximation space d P$_{k}$ , in which a key property from the preservation point of view relies on the equalities (5) which hold also at the discrete level.

2.2.2 Some complexes involving a discrete discontinuous space for velocities

We begin by providing some definitions on the mesh we are working on. The domain T 2 is divided into a set of cells C . The set of sides of the mesh is denoted by S . Each side of the mesh is oriented, and with respect to a given side S ∈ S , the left cell is the one for which the canonical normal is outgoing, whereas the right cell is the one in which the canonical normal is ingoing. The jump of a scalar f that is continuous on the cells and discontinuous across the faces is defined as

$$
\forall S \in \mathcal{S} \quad [f]_S := f_L - f_R,
$$

where f$_{L}$ is the value on the left cell and f$_{R}$ is the value on the right cell.
and

$$ \mathbf{d}\mathbb{B}_k^{\text{curl}} = \left( \begin{array}{c} \mathbf{d}\mathbb{Q}_{k,k} + \mathbf{d}\mathbb{Q}_{k-1,k+1} \\ \mathbf{d}\mathbb{Q}_{k,k} + \mathbf{d}\mathbb{Q}_{k+1,k-1} \end{array} \right) \oplus \text{Span} \left( \begin{array}{c} x^k y^{k+1} \\ x^{k+1} y^k \end{array} \right). $$

Last, the space $\mathbb{C}_k$ is

$$ \mathbb{C}_k = \mathrm{d}\check{\mathbb{Q}}_k(\mathcal{C}) \times \mathrm{d}\mathbb{P}_{k-1}(\mathcal{S}), $$

with

$$ \mathrm{d}\check{\mathbb{Q}}_k := \mathbf{d}\mathbb{Q}_{k,k-1} + \mathbf{d}\mathbb{Q}_{k-1,k}, $$

or more clearly, $\mathrm{d}\check{\mathbb{Q}}_k$ contains all canonical monomials of $\mathbf{d}\mathbb{Q}_k$ except for $x^k y^k$. As remarked in [71], the difference with the discontinuous Nédélec and Raviart-Thomas elements is, in this case, only of a single element, and the benefits of this optimal basis with respect to the classical one is less evident.

## 2.4 Scalar products and adjoint curl and divergence

The scalar product on the finite element spaces will be denoted with $\langle \cdot | \cdot \rangle$ with the finite element space as index, e.g. $\langle \cdot | \cdot \rangle_{\mathbb{A}_{k+1}}$. For the spaces $\mathbb{A}_{k+1}$, $\mathbf{d}\mathbb{B}_k^{\text{curl}}$ and $\mathbf{d}\mathbb{B}_k^{\text{div}}$, the scalar product is the standard $L^2$ product.

As the space $\mathbb{C}_k$ includes in general both components in the cell and components in the face, the scalar product in $\mathbb{C}_k$ is the graph scalar product:

$$ \langle f|g \rangle_{\mathbb{C}_k} := \sum_{c \in \mathcal{C}} \int_C f^C g^C + \sum_{S \in \mathcal{S}} \int_S f^S g^S. $$

The definition of these scalar products allow to define the adjoint differential operators as follows:

**Definition 1** (Adjoint differential operators). *Based on the complex (9), the adjoint gradient $\nabla^*$ is defined as*

$$ \begin{array}{l} \nabla^* : \mathbf{d}\mathbb{B}_k^{\text{curl}} \mapsto \mathbb{A}_{k+1} \\ \text{such that } \forall \mathbf{u} \in \mathbf{d}\mathbb{B}_k^{\text{curl}} \quad \forall f \in \mathbb{A}_{k+1} \qquad \langle \nabla^* \mathbf{u} | f \rangle_{\mathbb{A}_{k+1}} = \langle \mathbf{u} | \nabla^* f \rangle_{\mathbf{d}\mathbb{B}_k^{\text{curl}}}, \end{array} $$

whereas the adjoint curl $(\nabla^\perp)^*$ is defined as

$$ \begin{array}{l} (\nabla^\perp)^* : \mathbb{C}_k \mapsto \mathbf{d}\mathbb{B}_k^{\text{curl}} \\ \text{such that } \forall \mathbf{u} \in \mathbf{d}\mathbb{B}_k^{\text{curl}} \quad \forall f \in \mathbb{C}_k \qquad \langle (\nabla^\perp)^* f | \mathbf{u} \rangle_{\mathbf{d}\mathbb{B}_k^{\text{curl}}} = \langle f | \nabla^\perp \cdot \mathbf{u} \rangle_{\mathbb{C}_k}. \end{array} $$

Based on the complex (10), the adjoint curl $(\nabla^\perp)^*$ is defined as

$$ \begin{array}{l} (\nabla^\perp)^* : \mathbf{d}\mathbb{B}_k^{\text{div}} \mapsto \mathbb{A}_{k+1} \\ \text{such that } \forall \mathbf{u} \in \mathbf{d}\mathbb{B}_k^{\text{div}} \quad \forall f \in \mathbb{A}_{k+1} \qquad \langle (\nabla^\perp)^* \mathbf{u} | f \rangle_{\mathbb{A}_{k+1}} = \langle \mathbf{u} | \nabla^\perp f \rangle_{\mathbf{d}\mathbb{B}_k^{\text{div}}}, \end{array} $$

whereas the adjoint divergence $(\nabla^\cdot)^*$ is defined as

$$ \begin{array}{l} (\nabla^\cdot)^* : \mathbb{C}_k \mapsto \mathbf{d}\mathbb{B}_k^{\text{div}} \\ \text{such that } \forall \mathbf{u} \in \mathbf{d}\mathbb{B}_k^{\text{div}} \quad \forall f \in \mathbb{C}_k \qquad \langle (\nabla^\cdot)^* f | \mathbf{u} \rangle_{\mathbf{d}\mathbb{B}_k^{\text{div}}} = \langle f | \nabla^\cdot \mathbf{u} \rangle_{\mathbb{C}_k}. \end{array} $$

Note that at the continuous level, the gradient and the divergence are adjoint up to a sign, this means that at the discrete level, $-\nabla^*$ is an approximation of the divergence on $\mathbf{d}\mathbb{B}_k^{\text{curl}}$. In the same manner, $-(\nabla^\perp)^*$ is an approximation of the curl on $\mathbf{d}\mathbb{B}_k^{\text{div}}$.

These adjoint differential operators play a key role in this article, because they are the ones that will be proven to be naturally preserved by the discontinuous Galerkin method.

Last, we remark that it is possible to define these adjoint operators as the composition of exterior derivative with Hodge-star operators. This was done in the first version of this article [70], but was removed in the new version for clarity purpose.
<figure><img src="image_3.png" /><figcaption>Figure 2: Two dimensional discrete de-Rham complex on the proxies proposed in [54]. The approximation space $A_{k+1}$ is always a continuous finite element approximation space, namely $P_{k+1}$ for triangles and $Q_{k+1}$ for quadrangles. Depending on the properties needed, the choice of approximation space for $dB_k^{curl}$ and $dB_k^{div}$ may be very diverse. As the spaces $dB_k^{curl}$ and $dB_k^{div}$ include discontinuous functions, the application of the exterior derivative in the distribution sense (operators $\nabla_{\mathcal{D}'} \cdot$ and $\nabla_{\mathcal{D}'} \cdot$) includes the classical derivative located on the cells and the contribution of the normal or tangential jumps multiplied by a Dirac distribution located on the sides of the mesh. Therefore $C_k$ appears as a Cartesian product of a cell-based finite element space and a side-based finite element space.</figcaption></figure>

Recently, it was proposed to relax the continuity constraints induced by the conformal hypothesis [45]. Contrarily to the classical conformal approximation, the approximation of the vector spaces $H^{curl}$ and $H^{div}$ proposed in [45] are completely discontinuous, without any hypothesis on the continuity of the normal component (which holds for $RT$) or the tangential component (which holds for $N$). This idea was extended in [54] for finding basis with fewer degrees of freedom. The different approximation space for vectors proposed in [54] will be denoted by $dB_k^{div}$ and $dB_k^{curl}$. In this case, the classical derivation operators cannot be applied to these approximation space, and the derivative in the sense of distributions is used instead

$$ \forall \mathbf{u} \in dB_k^{div} \quad \nabla_{\mathcal{D}'} \cdot \mathbf{u} := \begin{cases} \forall C \in \mathcal{C} & (\nabla_{\mathcal{D}'} \cdot \mathbf{u})^C := \nabla \cdot \mathbf{u} \\ \forall S \in \mathcal{S} & (\nabla_{\mathcal{D}'} \cdot \mathbf{u})^S := -[\mathbf{u} \cdot \mathbf{n}] \end{cases} $$

and

$$ \forall \mathbf{u} \in dB_k^{curl} \quad \nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{u} := \begin{cases} \forall C \in \mathcal{C} & (\nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{u})^C := \nabla^{\perp} \cdot \mathbf{u} \\ \forall S \in \mathcal{S} & (\nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{u})^S := -[\mathbf{u}^{\perp} \cdot \mathbf{n}] \end{cases} $$

We see that the image of the operators $(\nabla_{\mathcal{D}'} \cdot)$ and $(\nabla_{\mathcal{D}'}^{\perp} \cdot)$ is in a Cartesian product of cell space and face space, that we will denote in general $C_k$. In [54], several such vectorial approximation spaces were proposed, for which the harmonic gap property [1, Chap. 5.2.3], which may be summarized as

**Definition 1** (Harmonic gap property). A discrete diagram ensures the harmonic gap property if the discrete and continuous cohomology spaces are isomorphic, which induces that (5) holds.

was proven on $\mathbb{T}^2$. The proposed approximation spaces may all be put in the diagram of Figure 2. In [54], the following approximation spaces were proposed; the first family is based on the conformal case, whereas the two other proposed approximation spaces are optimal in the number of degrees of freedom.

1. Approximation spaces directly based on discontinuous versions of Raviart-Thomas/Nédélec finite elements. In this version of the discontinuous spaces, proposed in [45], the polynomial basis are exactly the same as for the conformal approximation on each
## 2.5 Projection operator in the $\mathbb{C}_k$ space

The spaces $A$, $d\mathbb{B}_k^{\text{curl}}$ and $d\mathbb{B}_k^{\text{div}}$ are rather usual approximation spaces, and the classical $L^2$ projection is used on these spaces. For these spaces, we denote by $\pi$ with the approximation space as an exponent the $L^2$ projection: $\pi^{h_k}$, $\pi^{d\mathbb{B}_k^{\text{curl}}}$ and $\pi^{d\mathbb{B}_k^{\text{div}}}$.

The approximation space $\mathbb{C}_k$ is less classical, as it is a Cartesian product of degrees of freedom on the cells and degrees of freedom on the faces, still it is similar as the approximation space of the Hybrid High Order methods [27, 28, 26]. Therefore, following [27, 28, 26], we define the projection on $\mathbb{C}_k$ as follows

**Definition 2** (Projection operator on $\mathbb{C}_k$). Given a smooth function $g \in C^\infty(\Omega)$, the projection operator $\pi^{\mathbb{C}_k}$ is defined as

$$ \sum_{C \in \mathcal{C}} \int_C g \varphi^C + \sum_{S \in \mathcal{S}} \int_S g \varphi^S = \sum_{C \in \mathcal{C}} \int_C (\pi^{\mathbb{C}_k} g)^C \varphi^C + \sum_{S \in \mathcal{S}} \int_S (\pi^{\mathbb{C}_k} g)^S \varphi^S. $$

The projection operator defined in **Definition 2** ensures then the following proposition:

**Proposition 1** (Commutation property of the projection operator on $\mathbb{C}_k$). The following diagram

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

commutes, and so does the following diagram

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

*Proof.* The proof is performed for the first diagram, the proof for the second following the same lines. We denote by $g$ an element of $C^\infty(\Omega)$. Then for any $\mathbf{v} \in (C^\infty(\Omega))^2$, we have

$$ \int_{\Omega} \nabla \cdot (g \mathbf{v}) = \int_{\Omega} g \nabla \cdot \mathbf{v} + \int_{\Omega} \nabla g \cdot \mathbf{v}, $$

but as $\Omega$ is the two dimensional torus, the left hand side is zero. This directly gives

$$ (\nabla \cdot)^* g = -\nabla g. $$

The projection of $(-\nabla g)$ on $d\mathbb{B}_k^{\text{div}}$ ensures

$$ \forall \Psi \in d\mathbb{B}_k^{\text{div}} \quad -\sum_{C \in \mathcal{C}} \int_C \nabla g \cdot \Psi = -\sum_{C \in \mathcal{C}} \int_C \pi^{d\mathbb{B}_k^{\text{div}}} (\nabla g) \cdot \Psi = \sum_{C \in \mathcal{C}} \int_C \pi^{d\mathbb{B}_k^{\text{div}}} \circ \nabla^*(g) \cdot \Psi. \quad (11) $$

On the other hand, the definition of the projection of $g$ on $\mathbb{C}_k$ ensures

$$ \forall \varphi \in \mathbb{C}_k \quad \sum_{C \in \mathcal{C}} \int_C g \varphi^C + \sum_{S \in \mathcal{S}} \int_S g \varphi^S = \sum_{C \in \mathcal{C}} \int_C (\pi^{\mathbb{C}_k} g)^C \varphi^C + \sum_{S \in \mathcal{S}} \int_S (\pi^{\mathbb{C}_k} g)^S \varphi^S. $$
cell, but all the continuity constraints are relaxed. These spaces were denoted by $dRT_k$ and $dN_k$. This leads to

$$C_k = d\mathbb{P}_k(C) \times d\mathbb{P}_k(S),$$

for triangular meshes and to

$$C_k = d\mathbb{Q}_k(C) \times d\mathbb{P}_k(S),$$

for quadrangular meshes. These approximation space work fine, but some degrees of freedom are useless for ensuring the harmonic gap property. Also the Raviart-Thomas/Nédélec finite element basis should be generated which is not always straightforward on triangular meshes.

2. **Optimal approximation spaces on triangles.** On triangles, the classical approximation space for vectors, namely the one obtained by tensorizing the classical approximation space for scalars can be put in the diagram Figure 2, as addressed in [54]; in this case, $d\mathbb{B}_k^{\text{curl}} = d\mathbb{B}_k^{\text{div}} = d\mathbb{P}_k$, and

$$C_k = d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(S).$$

3. **Optimal approximation spaces on quadrangles.** The case of quadrangular meshes is quite surprising because the lowest order approximation space (namely with $A_{k+1} = Q_1$) does not give the classical finite volume vector space, but rather an enriched version including three basis vector instead of two:

$$d\mathbb{B}_0^{\text{div}} = \text{Span} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} -x \\ y \end{pmatrix} \right),$$

and

$$d\mathbb{B}_0^{\text{curl}} = \text{Span} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} y \\ x \end{pmatrix} \right).$$

In the general case, we denote by $d\mathbb{Q}_{i,j}$ the space of polynomials of degree lower than $i$ in $x$ and lower than $j$ in $y$, and we set

$$d\mathbb{B}_k^{\text{div}} = [(d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1}) \times (d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1})] \oplus \text{Span} \left( \begin{array}{c} -x^{k+1}y^k \\ x^k y^{k+1} \end{array} \right)$$

and

$$d\mathbb{B}_k^{\text{curl}} = [(d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1}) \times (d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1})] \oplus \text{Span} \left( \begin{array}{c} x^k y^{k+1} \\ x^{k+1} y^k \end{array} \right).$$

Last, the space $C_k$ is

$$C_k = d\tilde{\mathbb{Q}}_k(C) \times d\mathbb{P}_{k-1}(S),$$

with

$$d\tilde{\mathbb{Q}}_k := d\mathbb{Q}_{k,k-1} + d\mathbb{Q}_{k-1,k},$$

or more clearly, $d\tilde{\mathbb{Q}}_k$ contains all canonical monomials of $d\mathbb{Q}_k$ except for $x^k y^k$. As remarked in [54], the difference with the discontinuous Nédélec and Raviart-Thomas elements is, in this case, only of a single element, and the benefits of this optimal basis with respect to the classical one is less evident.

The scalar product on the finite element spaces will be denoted with $\langle \cdot | \cdot \rangle$ with the finite element space as index, e.g. $\langle \cdot | \cdot \rangle_{A_{k+1}}$. When the space of differential forms will be needed, we will take the same notation, but with gothic letters, e.g. $\langle \cdot | \cdot \rangle_{\mathfrak{A}_{k+1}}$ is the scalar product on the discrete space of differential forms which proxies are in $A_{k+1}$.

Last, as the space $C_k$ includes in general both components in the cell and components in the face, the scalar product in $C_k$ is the graph scalar product:

$$\langle f | g \rangle_{C_k} := \sum_{c \in C_k} \int_C f^C g^C + \sum_{S \in S_k} \int_S f^S g^S.$$
Using the definition of $(\nabla_{\mathcal{G}'\cdot})^*$ gives

$$ \forall \Psi \in \mathbf{d}\mathcal{B}_k^{\text{div}} \quad \sum_{C \in \mathcal{C}} \int_C (\nabla_{\mathcal{G}'\cdot})^* \circ \pi^{\mathcal{C}_k}(g) \cdot \Psi = \sum_{C \in \mathcal{C}} \int_C (\pi^{\mathcal{C}_k} g)^C \nabla \cdot \Psi - \sum_{S \in \mathcal{S}} \int_S (\pi^{\mathcal{C}_k} g)^S \llbracket \Psi \cdot n_S \rrbracket. $$

Combining the two last equalities provides

$$ \forall \Psi \in \mathbf{d}\mathcal{B}_k^{\text{div}} \quad \sum_{C \in \mathcal{C}} \int_C (\nabla_{\mathcal{G}'\cdot})^* \circ \pi^{\mathcal{C}_k}(g) \cdot \Psi = \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \Psi - \sum_{S \in \mathcal{S}} \int_S g \llbracket \Psi \cdot n \rrbracket. \quad (12) $$

It remains to apply the Stokes formula on each cell C for getting

$$ \begin{aligned} -\sum_{C \in \mathcal{C}} \int_C \nabla g \cdot \Psi &= -\sum_{C \in \mathcal{C}} \int_C \nabla \cdot (g\Psi) + \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \Psi \\ &= -\sum_{C \in \mathcal{C}} \sum_{S \in \mathcal{S}} \int_S g \Psi \cdot n^{\text{out}} + \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \Psi \\ &= -\sum_{S \in \mathcal{S}} \int_S \llbracket g \Psi \cdot n_S \rrbracket + \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \Psi \\ -\sum_{C \in \mathcal{C}} \int_C \nabla g \cdot \Psi &= -\sum_{S \in \mathcal{S}} \int_S g \llbracket \Psi \cdot n_S \rrbracket + \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \Psi. \end{aligned} $$

This last equality provides the equality of the left hand side of (11) with the right hand side of (12), which induces the equality of the right hand side of (11) and the right hand side of (12), which ends the proof. □

Note that the second commutative diagram of Proposition 1 will be relevant for initializing divergence free fields in the discrete adjoint sense in section 4; especially, starting from a $g \in C^\infty(\Omega)$ the path $\pi^{\mathbf{d}\mathcal{B}_k^{\text{curl}}} \circ (\nabla^\perp \cdot)^*(g)$ ensures a high order approximation of $\nabla^\perp g$, whereas the path $(\nabla_{\mathcal{G}'\cdot})^* \circ \pi^{\mathcal{C}_k}(g)$ ensures that the discrete field is divergence free in the discrete adjoint sense.

# 3 Discrete preservation of the curl or the divergence

## 3.1 Numerical scheme

In this section, we consider a system of conservation law in which one of the unknowns is a vector $\mathbf{u}$. We denote by $\bar{\mathbf{G}}$ the flux of the conservation law involving $\mathbf{u}$.

$$ \partial_t \mathbf{u} + \nabla \cdot \bar{\mathbf{G}} = 0. \quad (13) $$

We consider the numerical resolution of (13) on a periodic domain. The discontinuous Galerkin method reads

$$ \text{Find } \mathbf{u} \in \mathcal{B} \quad \forall \mathbf{v} \in \mathcal{B} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C \bar{\mathbf{G}} : \nabla \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{v} \rrbracket \cdot \tilde{\mathbf{G}} = 0, \quad (14) $$

where $\tilde{\mathbf{G}}$ is the numerical flux, and $\mathcal{B}$ is the approximation space for vectors that is not defined for the moment.

## 3.2 Discrete preservation of the vorticity

In this section, we consider the case in which $\bar{\mathbf{G}}$ is proportional to the identity $\bar{\mathbf{G}} = gI_d$, which means that we have $\nabla \cdot \bar{\mathbf{G}} = \nabla g$. This means that the conservation law on $\mathbf{u}$ (13) may be simplified as

$$ \partial_t \mathbf{u} + \nabla g = 0, \quad (15) $$
# 3 Discrete Hodge-star operator, codifferentials and Hodge Laplacian

## 3.1 Discrete Hodge star operator

Our discrete de-Rham complex includes the 0, 1 and 2-forms. For completing the discrete counterpart of the continuous de-Rham complex, we wish to define the discrete Hodge-star operators, which we denote by $\star_{\{k\}}$ for the Hodge-star operator mapping the $k$ forms to the $d-k$ forms. In dimension 2, the Hodge star operator $\star_{\{1\}}$ maps the 1-forms to themselves, and is translated in terms of proxy as a $\pi/2$ rotation. We define the discrete operator $\star_{\{1\}}$ also as a $\pi/2$ rotation.

We now define the discrete Hodge star operators $\star_{\{0\}}$ and $\star_{\{2\}}$ between the discrete spaces $A_{k+1}$ and $C_k$. If the approximation degrees are considered, we see that the space $A_{k+1}$ is of a higher degree per cell and per side, which means that a projection should be a good option for the definition of the Hodge star operator $\star_{\{0\}}$:

**Definition 2** (Definition of $\star_{\{0\}}$). The Hodge star operator $\star_{\{0\}}$ that maps $A_{k+1}$ to $C_k$ is defined as

$$ \begin{array}{rcl} \star_{\{0\}} : & A_{k+1} & \mapsto C_k \\ & \varphi & \mapsto \star_{\{0\}}(\varphi) := \begin{cases} \mathcal{P}_C(\varphi) & \forall C \in \mathcal{C} \\ \mathcal{P}_S(\varphi) & \forall S \in \mathcal{S} \end{cases} \end{array} \quad (6) $$

where $\mathcal{P}$ is the $L^2$ projection.

We now define the discrete Hodge star operator $\star_{\{2\}}$. For a given $\Psi \in C_k$, we wish to define $\star_{\{2\}}$ such that

$$ \forall \varphi \in A_{k+1} \quad \langle \star_{\{2\}} (\Psi) | \varphi \rangle_{A_{k+1}} = \langle \Psi | \star_{\{0\}} \varphi \rangle_{C_k} $$

This can be developed on one hand as

$$ \langle \star_{\{2\}} (\Psi) | \varphi \rangle_{A_{k+1}} = \sum_{C \in \mathcal{C}} \int_C \star_{\{2\}} (\Psi) \varphi, $$

and on the other hand as

$$ \langle \Psi | \star_{\{0\}} \varphi \rangle_{C_k} = \sum_{C \in \mathcal{C}} \int_C \Psi \star_{\{0\}} (\varphi) + \sum_{S \in \mathcal{S}} \int_S \Psi \star_{\{0\}} (\varphi). $$

By definition of the projection, this may be simplified as

$$ \langle \Psi | \star_{\{0\}} (\varphi) \rangle_{C_k} = \sum_{C \in \mathcal{C}} \int_C \Psi^C \varphi + \sum_{S \in \mathcal{S}} \int_S \Psi^S \varphi, $$

and this leads to the following definition:

**Definition 3** (Definition of $\star_{\{2\}}$). The Hodge star operator $\star_{\{2\}}$ that maps $C_k$ to $A_{k+1}$ is defined as

$$ \forall \varphi \in A_{k+1} \quad \forall \psi \in C_k \quad \sum_{C \in \mathcal{C}} \int_C \star_{\{2\}} (\psi) \varphi = \sum_{C \in \mathcal{C}} \int_C \psi^C \varphi + \sum_{S \in \mathcal{S}} \int_S \psi^S \varphi. $$

Note that the computation of $\star_{\{2\}}$ as in Definition 3 requires the inversion of the mass matrix of $A_{k+1}$.

**Remark 1.** Following [54], we always have

$$ \dim C_k > \dim A_{k+1}, $$

which means that neither $\star_{\{0\}}$ nor $\star_{\{2\}}$ can be invertible.

Also, following the rank nullity theorem, we have

$$ \dim (\ker \star_{\{2\}}) = \dim C_k - \dim (\star_{\{2\}}) \ge \dim C_k - \dim A_{k+1} > 0, $$

which means that $\star_{\{2\}}$ has always a non trivial kernel.

It is however easy to check that $\star_{\{2\}} \circ \star_{\{0\}} = I_d$.
where $g$ is a scalar that may depend nonlinearly on the variables of the system. Then the discontinuous Galerkin discretization reads

$$ \text{Find } \mathbf{u} \in \mathbf{d}\mathbf{B}_k^{\text{div}} \quad \forall \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{div}} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{v} \rrbracket \cdot \tilde{\mathbf{G}} = 0, \quad (16) $$

where $\tilde{\mathbf{G}}$ is the numerical flux. Equation (15) induces formally the conservation of $\nabla^{\perp} \cdot \mathbf{u}$, and the discontinuous Galerkin method (16) may have a similar property summarized in the following proposition.

**Proposition 2** (Conservation of $(\nabla^{\perp} \mathbf{u})^*$). Consider the numerical scheme (16). If $\tilde{\mathbf{G}}$ is parallel to $\mathbf{n}_S$, then the numerical scheme (16) induces $\partial_t ((\nabla^{\perp})^* \mathbf{u}) = 0$.

*Proof.* If $\tilde{\mathbf{G}}$ is parallel to $\mathbf{n}_S$, then the numerical flux may be rewritten

$$ \tilde{\mathbf{G}} = \tilde{g} \mathbf{n}_S, $$

where $\tilde{g}$ is a scalar, so that the numerical scheme (16) may be rewritten

$$ \text{Find } \mathbf{u} \quad \forall \mathbf{v} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{v} \cdot \mathbf{n}_S \rrbracket \tilde{g} = 0. $$

In this last equation, suppose that $\mathbf{v} = \nabla^{\perp} f$ for $f \in A_{k+1}$. Then the normal jump vanishes at each side:

$$ \forall S \in \mathcal{S} \quad \llbracket \mathbf{v} \cdot \mathbf{n}_S \rrbracket = 0, $$

and on each cell, its strong divergence is 0:

$$ \nabla \cdot \mathbf{v} = 0, $$

so that

$$ \forall f \in A_{k+1} \quad \sum_{C \in \mathcal{C}} \int_C (\nabla^{\perp} f) \cdot \partial_t \mathbf{u} = 0. $$

This last equation means that

$$ \forall f \in A_{k+1} \quad \langle \nabla^{\perp} f | \partial_t \mathbf{u} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = 0. $$

By definition of $(\nabla^{\perp})^*$, we get

$$ \forall f \in A_{k+1} \quad \langle f | (\nabla^{\perp})^* (\partial_t \mathbf{u}) \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = 0. $$

It remains to prove that $(\nabla^{\perp})^* (\partial_t \mathbf{u}) = \partial_t ((\nabla^{\perp})^* \mathbf{u})$. Denoting by $f$ an element of $A_{k+1}$, we have

$$ \begin{aligned} \langle \partial_t ((\nabla^{\perp})^* \mathbf{u}) | f \rangle_{A_{k+1}} &= \int_{\Omega} \partial_t ((\nabla^{\perp})^* \mathbf{u}) f \\ &= \partial_t \left( \int_{\Omega} ((\nabla^{\perp})^* \mathbf{u}) f \right) \\ &= \partial_t \left( \int_{\Omega} \mathbf{u} \cdot \nabla^{\perp} f \right) \\ &= \int_{\Omega} (\partial_t \mathbf{u}) \cdot \nabla^{\perp} f \\ &= \langle \partial_t \mathbf{u} | \nabla^{\perp} f \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} \\ \langle \partial_t ((\nabla^{\perp})^* \mathbf{u}) | f \rangle_{A_{k+1}} &= \langle ((\nabla^{\perp})^* (\partial_t \mathbf{u})) | f \rangle_{A_{k+1}}, \end{aligned} $$

which ends the proof.

□
## 3.2 Definition of the codifferential operators

### 3.2.1 Definition of $\delta^1$

The continuous codifferential $\delta^1$ was defined in (3). We wish to find a formula for the discrete codifferential on the proxies. We focus on the proxy choice (2).

**Proposition 1** (Definition of $(\nabla^{\perp})^*$. $(\nabla^{\perp})^*$ is defined as in (3) but with the discrete operators. It ensures

$$ \forall \mathbf{u} \in \mathbf{dB}_{k}^{\text{div}} \quad \forall f \in A_{k+1} \quad \langle (\nabla^{\perp})^{*} \mathbf{u} | f \rangle_{A_{k+1}} = \langle \mathbf{u} | \nabla^{\perp} f \rangle_{\mathbf{dB}_{k}^{\text{div}}}. $$

*Proof.* The composition of $d^1$ with $\star_{\{1\}}$ appearing in (3) is a $(\nabla^{\perp})$ on the proxy, which is translated into a $(\nabla^{\perp}_{\mathcal{G}', \cdot})$ at a discrete level. Then for $\mathbf{u} \in \mathbf{dB}_{k}^{\text{div}}$,

$$ \text{Proxy} (d^1 \star_{\{1\}} \mathbf{u}) = \begin{cases} \forall C \in \mathcal{C} & \nabla^{\perp} \cdot \mathbf{u} \\ \forall S \in \mathcal{S} & \llbracket \mathbf{u}^{\perp} \cdot \mathbf{n} \rrbracket \end{cases} $$

Then

$$ \begin{aligned} \langle (\nabla^{\perp})^* \mathbf{u} | f \rangle_{A_{k+1}} &= \langle \star_{\{2\}} d^1 \star_{\{1\}} \mathbf{u} | f \rangle_{A_{k+1}} \\ &= \langle \star_{\{2\}} \text{Proxy} (d^1 \star_{\{1\}} \mathbf{u}) | f \rangle_{A_{k+1}} \\ &= \sum_{C \in \mathcal{C}} \int_C f \nabla^{\perp} \cdot \mathbf{u} + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{u}^{\perp} \cdot \mathbf{n} \rrbracket f \\ &= - \sum_{C \in \mathcal{C}} \int_C f \nabla \cdot \mathbf{u}^{\perp} + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{u}^{\perp} \cdot \mathbf{n} \rrbracket f \\ &= - \sum_{C \in \mathcal{C}} \int_C \nabla \cdot (f \mathbf{u}^{\perp}) - \sum_{C \in \mathcal{C}} \int_C \mathbf{u}^{\perp} \cdot \nabla f + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{u}^{\perp} \cdot \mathbf{n} \rrbracket f \\ &= - \sum_{S \in \mathcal{S}} \int_S f \llbracket \mathbf{u}^{\perp} \cdot \mathbf{n} \rrbracket - \sum_{C \in \mathcal{C}} \int_C \mathbf{u}^{\perp} \cdot \nabla f + \sum_{S \in \mathcal{S}} \int_S \llbracket \mathbf{u}^{\perp} \cdot \mathbf{n} \rrbracket f \\ &= - \sum_{C \in \mathcal{C}} \int_C \mathbf{u}^{\perp} \cdot \nabla f \\ &= \sum_{C \in \mathcal{C}} \int_C \mathbf{u} \cdot \nabla^{\perp} f \\ \langle (\nabla^{\perp})^* \mathbf{u} | f \rangle_{A_{k+1}} &= \langle \mathbf{u} | \nabla^{\perp} f \rangle_{\mathbf{dB}_{k}^{\text{div}}}, \end{aligned} $$

which ends the proof.

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

### 3.2.2 Definition of $\delta^2$

The continuous codifferential $\delta^2$ was defined in (4). Still, $\delta^2$ depends on $\star_{\{2\}}$, and as stated in Remark 1, $\star_{\{2\}}$ has a non trivial kernel, which may lead to a mismatch between the properties of the discrete and continuous codifferentials. This is why we define the discrete codifferential $\delta^2$ as follows

**Definition 4** (Definition of $(-\nabla_{\mathcal{G}', \cdot})^*$. $(-\nabla_{\mathcal{G}', \cdot})^*$ is defined as

$$ \forall \mathbf{u} \in \mathbf{dB}_{k}^{\text{div}} \quad \forall f \in C_k \quad \langle (-\nabla_{\mathcal{G}', \cdot})^* f | \mathbf{u} \rangle_{\mathbf{dB}_{k}^{\text{div}}} = \langle f | -\nabla_{\mathcal{G}', \cdot} \mathbf{u} \rangle_{C_k}. $$

We however still have an identity similar to (4) at the discrete level

**Proposition 2.** If the discrete proxy of $\delta^2$ is defined as in Definition 4, then

$$ \delta^2 \circ \star_{\{0\}} = - \star_{\{1\}} \circ d^0. $$
**Remark 1** (full discretization in space: effect of quadrature formula). The Proposition 2 was proven in the non fully spatial discretized case, namely without the use of quadrature formula. It is important to note that both the cell and face integral of the spatial discretisation are vanishing pointwise, therefore they vanish whatever the quadrature formula used. Concerning the adjoint curl that is preserved, the numerical computation of $(\nabla^{\perp})^*$ should take exactly the same quadra- ture formula as the one used for computing the mass matrix. For example, if the mass matrix is lumped, then the adjoint curl should be defined with the same lumped mass matrix.

**Remark 2** (Time discretization). Proposition 2 addresses the semi-discrete case (discrete in space, but not in time). A proof for the fully discrete case can be done by following exactly the same lines as the proof of Proposition 2, by replacing the time derivative by its discretization, still this leads to a cumbersome proof. Actually, the extension to the time discrete case is a simple consequence of [62, Section 3], which states that if the semi-discrete scheme preserves the discrete curl, then all B-series (Butcher series) [41, 61] time integrators (including all Runge-Kutta ones), also preserve the discrete curl.

**Remark 3** (Lax-Friedrich flux with normal diffusion). A widely used numerical flux is the Lax-Friedrich flux which reads

$$ \tilde{\mathbf{G}}(\overline{\mathbf{G}}_L, \mathbf{u}_L, \overline{\mathbf{G}}_R, \mathbf{u}_R, \mathbf{n}_S) = \frac{\overline{\mathbf{G}}_L \mathbf{n}_S + \overline{\mathbf{G}}_R \mathbf{n}_S}{2} + \frac{\lambda}{2}(\mathbf{u}_L - \mathbf{u}_R), \quad (17) $$

where $\lambda$ is the maximum absolute value of the eigenvalues of the system. Dealing with the system (15), the numerical flux (17) can be simplified as

$$ \tilde{\mathbf{G}}(g_L, \mathbf{u}_L, g_R, \mathbf{u}_R, \mathbf{n}_S) = \frac{g_L \mathbf{n}_S + g_R \mathbf{n}_S}{2} + \frac{\lambda}{2}(\mathbf{u}_L - \mathbf{u}_R). $$

The centered part of the flux is clearly parallel to the normal, and should not be modified because it ensures the consistency of the numerical scheme. The diffusive part may be decomposed into its normal and tangential part as

$$ (\mathbf{u}_L - \mathbf{u}_R) = \mathbf{n}_S \mathbf{n}_S^T (\mathbf{u}_L - \mathbf{u}_R) + (\mathbf{I}_d - \mathbf{n}_S \mathbf{n}_S^T) (\mathbf{u}_L - \mathbf{u}_R). $$

Therefore an easy way to ensure the hypothesis of Proposition 2 is to use the following Lax-Friedrich flux with purely normal diffusion

$$ \tilde{\mathbf{G}}(g_L, \mathbf{u}_L, g_R, \mathbf{u}_R, \mathbf{n}_S) = \frac{g_L \mathbf{n}_S + g_R \mathbf{n}_S}{2} + \frac{\lambda \mathbf{n}_S \mathbf{n}_S^T}{2} (\mathbf{u}_L - \mathbf{u}_R). \quad (18) $$

**Remark 4** (Link with the Hodge Laplacian). Putting diffusion operator to the equation (13) consists, at least formally, in adding a vector Laplacian $\Delta$ to the right hand side:

$$ \Delta \mathbf{u} = \begin{pmatrix} \partial_{xx} \mathbf{u}_x + \partial_{yy} \mathbf{u}_x \\ \partial_{xx} \mathbf{u}_y + \partial_{yy} \mathbf{u}_y \end{pmatrix} = \nabla (\nabla \cdot \mathbf{u}) + \nabla^{\perp} (\nabla^{\perp} \cdot \mathbf{u}). \quad (19) $$

For ensuring preservation of the vorticity, it is sufficient to put a diffusion that preserves also the vorticity, namely only the $\nabla (\nabla \cdot \mathbf{u})$ component. In [70, Remark 4], we proved that the diffusion induced by (18) matches with a discrete version of the operator $\nabla (\nabla \cdot \mathbf{u})$.

### 3.3 Discrete preservation of the irrotational component

In this section, we consider the case in which $\overline{\mathbf{G}}$ is antisymmetric; in dimension 2, $\overline{\mathbf{G}}$ may be written as

$$ \begin{pmatrix} 0 & -g \\ g & 0 \end{pmatrix}, $$

where $g$ is a scalar that may depend nonlinearly on the variables of the system. This means that $\mathbf{u}$ ensures the following conservation law

$$ \partial_t \mathbf{u} + \nabla^{\perp} g = 0. \quad (20) $$
Proof. We denote by $f \in A_{k+1}$ and $u \in \mathfrak{B}_k$. Then

$$
\begin{align*}
\langle \delta^2 \circ \star_{\{0\}} f | u \rangle_{\mathfrak{B}_k} &= \langle (-\nabla_{\mathscr{G}'\cdot})^* \star_{\{0\}} f | u \rangle_{\mathbf{dB}_k^{\text{div}}} \\
&= - \langle \star_{\{0\}} f | \nabla_{\mathscr{G}'\cdot} u \rangle_{\mathscr{C}_k} \\
&= - \sum_{c \in C} \int_C f \nabla \cdot u + \sum_{S \in S} \int_S f [[u \cdot n]] \\
&= - \sum_{c \in C} \int_C \nabla \cdot (f u) + \sum_{c \in C} \int_C u \cdot \nabla f + \sum_{S \in S} \int_S f [[u \cdot n]] \\
&= \sum_{c \in C} \int_C u \cdot \nabla f \\
&= - \sum_{c \in C} \int_C u \cdot (\nabla^{\perp} f)^{\perp}
\end{align*}
\[
\langle \delta^2 \circ \star_{\{0\}} f | u \rangle_{\mathfrak{B}_k} = - \langle \star_{\{1\}} \circ d^0 f | u \rangle_{\mathfrak{B}_k},
\]
$$

which ends the proof.

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

## 3.3 Hodge Laplacian

Even if we are dealing with hyperbolic systems, the Laplace operator plays a role as a regularization operator, this is why a small section is dedicated to this operator. The Hodge-Laplacian is defined as $\delta^{i+1} \circ d^i + d^{i-1} \circ \delta^i$. We are interested in this operator only on $\mathfrak{B}_k$. As we are dealing with Galerkin methods, we are rather interested in the $L^2$ product of the Hodge Laplacian of a $u \in \mathfrak{B}_k$ with a given $v \in \mathfrak{B}_k$.

**Proposition 3** (Expression of $\delta^2 \circ d^1$). For (1) as a choice of proxy, for all $u, v \in \mathfrak{B}_k$

$$
\begin{align}
\langle \delta^2 \circ d^1 u | v \rangle_{\mathfrak{B}_k} &= \langle (\nabla^{\perp})^* \nabla^{\perp} \cdot u | v \rangle_{\mathbf{dB}_k^{\text{curl}}} \tag{7} \\
&= \sum_{c \in C} \int_C (\nabla^{\perp} \cdot u) (\nabla^{\perp} \cdot v) + \sum_{S \in S} \int_S [[u]^T (I_2 - n^T n) [[v]]]. \nonumber
\end{align}
$$

For (2) as a choice of proxy, for all $u, v \in \mathfrak{B}_k$

$$
\begin{align}
\langle \delta^2 \circ d^1 u | v \rangle_{\mathfrak{B}_k} &= \langle (\nabla \cdot)^* \nabla \cdot u | v \rangle_{\mathbf{dB}_k^{\text{curl}}} \nonumber \\
&= \sum_{c \in C} \int_C (\nabla \cdot u) (\nabla \cdot v) + \sum_{S \in S} \int_S [[u]^T (n^T n) [[v]]]. \tag{8}
\end{align}
$$

Proof. We first deal with the choice (1) of proxies:

$$
\begin{align*}
\langle \delta^2 \circ d^1 u | v \rangle_{\mathfrak{B}_k} &= \langle (\nabla^{\perp})^* \nabla^{\perp} \cdot u | v \rangle_{\mathbf{dB}_k^{\text{curl}}} \\
&= \langle \nabla^{\perp} \cdot u | \nabla^{\perp} \cdot v \rangle_{\mathbf{dB}_k^{\text{curl}}} \\
&= \sum_{c \in C} \int_C (\nabla^{\perp} \cdot u) (\nabla^{\perp} \cdot v) + \sum_{S \in S} \int_S [[u \cdot n^{\perp}]] [[v \cdot n^{\perp}]] \\
\langle \delta^2 \circ d^1 u | v \rangle_{\mathfrak{B}_k} &= \sum_{c \in C} \int_C (\nabla^{\perp} \cdot u) (\nabla^{\perp} \cdot v) + \sum_{S \in S} \int_S [[u]^T n^{\perp} (n^{\perp})^T [[v]]],
\end{align*}
$$

and we have

$$
\mathbf{n}^{\perp} (\mathbf{n}^{\perp})^{\mathrm{T}} = \begin{pmatrix} -\mathbf{n}_y \\ \mathbf{n}_x \end{pmatrix} (-\mathbf{n}_y \quad \mathbf{n}_x) = \begin{pmatrix} \mathbf{n}_y^2 & -\mathbf{n}_x \mathbf{n}_y \\ -\mathbf{n}_x \mathbf{n}_y & \mathbf{n}_x^2 \end{pmatrix} = I_2 - \mathbf{n}^T \mathbf{n},
$$

which finally gives (7). (8) may be obtained in the same manner.

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

The case of $d^0 \circ \delta^1$ is harder to address, because of the computation of $\delta^1$. As it maps to $A_{k+1}$, its computation requires the inversion of the mass matrix of $A_{k+1}$ which is a continuous finite element space and so is not local. It is probably possible to simplify it by using mass lumping, still for the remaining of the article, only the expression of $\delta^2 \circ d^1$ matters. This is why the case of $d^0 \circ \delta^1$ is not further investigated.
Then the discontinuous Galerkin discretization reads

$$ \text{Find } \mathbf{u} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \quad \forall \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C g \nabla^{\perp} \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v}] \cdot \tilde{\mathbf{G}} = 0, \quad (21) $$

where $\tilde{\mathbf{G}}$ is the numerical flux.

**Proposition 3** (Adjoint conservation of the divergence). Consider the numerical scheme (21). If $\tilde{\mathbf{G}}$ is orthogonal to $\mathbf{n}_S$, then $(\nabla)^* \mathbf{u}$ is conserved by the numerical scheme.

The proof of Proposition 3 follows exactly the same lines as Proposition 2 and is not repeated.

**Remark 5** (Lax-Friedrich flux with purely tangential diffusion). The Lax-Friedrich numerical flux (17) for (20) may be simplified as

$$ \tilde{\mathbf{G}}(\mathbf{g}_L, \mathbf{u}_L, \mathbf{g}_R, \mathbf{u}_R, \mathbf{n}_S) = \frac{g_L \mathbf{n}_S^{\perp} + g_R \mathbf{n}_S^{\perp}}{2} + \frac{\lambda}{2} (\mathbf{u}_L - \mathbf{u}_R). $$

The centered part of the flux, which is clearly orthogonal to $\mathbf{n}_S$, ensures the consistency of the numerical scheme. Still relying on the decomposition of the diffusive part into its normal and tangential part, an easy way to ensure the hypothesis of Proposition 3 consists in using the following *Lax-Friedrich flux with purely tangential diffusion*

$$ \tilde{\mathbf{G}}(\mathbf{g}_L, \mathbf{u}_L, \mathbf{g}_R, \mathbf{u}_R, \mathbf{n}_S) = \frac{g_L \mathbf{n}_S^{\perp} + g_R \mathbf{n}_S^{\perp}}{2} + \frac{\lambda (\mathbf{I}_d - \mathbf{n}_S \mathbf{n}_S^T)}{2} (\mathbf{u}_L - \mathbf{u}_R). \quad (22) $$

Similarly to Remark 4, we proved in [70] that the diffusion of (22) matches with a discrete version of the operator $\nabla^{\perp} (\nabla^{\perp} \cdot)$, which is the component of the vector Laplacian (19) that preserves the divergence.

# 4 Preservation of initially divergence free field for the induction equation

In this section, we are interested in the induction equation, which reads, in dimension 3:

$$ \partial_t \mathbf{u} + \nabla_{3d} \times (\mathbf{b} \times \mathbf{u}) + \mathbf{b} \nabla_{3d} \cdot \mathbf{u} = 0, $$

where $\mathbf{u}$ is the unknown, $\mathbf{b}$ is a given vector field, and $\nabla_{3d}$ denotes the three dimensional nabla operator. We suppose that the two vectors $\mathbf{b}$ and $\mathbf{u}$ have only a $x$ and a $y$ component, and we suppose also that these components depend only on the variables $x$ and $y$. Then

$$ \begin{aligned} \nabla_{3d} \times (\mathbf{b} \times \mathbf{u}) &= \nabla_{3d} \times \left( \begin{pmatrix} \mathbf{b}_x \\ \mathbf{b}_y \\ 0 \end{pmatrix} \times \begin{pmatrix} \mathbf{u}_x \\ \mathbf{u}_y \\ 0 \end{pmatrix} \right) \\ &= \nabla_{3d} \times \begin{pmatrix} 0 \\ 0 \\ \mathbf{b}_x \mathbf{u}_y - \mathbf{b}_y \mathbf{u}_x \end{pmatrix} \\ &= \begin{pmatrix} \partial_y (\mathbf{b}_x \mathbf{u}_y - \mathbf{b}_y \mathbf{u}_x) \\ -\partial_x (\mathbf{b}_x \mathbf{u}_y - \mathbf{b}_y \mathbf{u}_x) \\ 0 \end{pmatrix} \\ \nabla_{3d} \times (\mathbf{b} \times \mathbf{u}) &= \begin{pmatrix} \nabla^{\perp} (\mathbf{b} \cdot \mathbf{u}^{\perp}) \\ 0 \end{pmatrix}. \end{aligned} $$

As a consequence, the two dimensional induction equation reads

$$ \partial_t \mathbf{u} + \nabla^{\perp} (\mathbf{b} \cdot \mathbf{u}^{\perp}) + \mathbf{b} \nabla \cdot \mathbf{u} = 0. \quad (23) $$
# 4 Discrete preservation of the curl or the divergence

## 4.1 Numerical scheme

In this section, we consider a system of conservation law in which one of the unknowns is a vector **u**. We denote by $\overline{\mathbf{G}}$ the flux of the conservation law involving **u**

$$ \partial_t \mathbf{u} + \nabla \cdot \overline{\mathbf{G}} = 0. \qquad (9) $$

We consider the numerical resolution of (9) on a periodic domain. The discontinuous Galerkin method reads

$$ \text{Find } \mathbf{u} \in \mathbf{B} \quad \forall \mathbf{v} \in \mathbf{B} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C \overline{\mathbf{G}} : \nabla \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v}] \cdot \tilde{\mathbf{G}} = 0, \quad (10) $$

where $\tilde{\mathbf{G}}$ is the numerical flux, and $\mathbf{B}$ is the approximation space for vectors that is not defined for the moment.

## 4.2 Discrete preservation of the vorticity

In this section, we consider the case in which $\overline{\mathbf{G}}$ is proportional to the identity $\overline{\mathbf{G}} = gI_d$, which means that we have $\nabla \cdot \overline{\mathbf{G}} = \nabla g$. This means that the conservation law on $\mathbf{u}$ (9) may be simplified as

$$ \partial_t \mathbf{u} + \nabla g = 0, \qquad (11) $$

where $g$ is a scalar that may depend nonlinearly on the variables of the system. Then the discontinuous Galerkin discretization reads

$$ \text{Find } \mathbf{u} \in \mathbf{d}\mathbf{B}_k^{\text{div}} \quad \forall \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{div}} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v}] \cdot \tilde{\mathbf{G}} = 0, \quad (12) $$

where $\tilde{\mathbf{G}}$ is the numerical flux. Equation (11) induces formally the conservation of $\nabla^{\perp} \cdot \mathbf{u}$, and the discontinuous Galerkin method (12) may have a similar property summarized in the following proposition.

**Proposition 4 (Conservation of $(\nabla^{\perp} \mathbf{u})^*)$.** Consider the numerical scheme (12). If $\tilde{\mathbf{G}}$ is parallel to $\mathbf{n}$, then the numerical scheme (12) induces $\partial_t ((\nabla^{\perp})^* \mathbf{u}) = 0$.

*Proof.* If $\tilde{\mathbf{G}}$ is parallel to $\mathbf{n}$, then the numerical flux may be rewritten

$$ \tilde{\mathbf{G}} = \tilde{g}\mathbf{n}, $$

where $\tilde{g}$ is a scalar, so that the numerical scheme (12) may be rewritten

$$ \text{Find } \mathbf{u} \quad \forall \mathbf{v} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C g \nabla \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v} \cdot \mathbf{n}] \tilde{g} = 0. $$

In this last equation, suppose that $\mathbf{v} = \nabla^{\perp} f$ for $f \in A_{k+1}$. Then the normal jump vanishes at each side:

$$ \forall S \in \mathcal{S} \quad [\mathbf{v} \cdot \mathbf{n}] = 0, $$

and on each cell, its strong divergence is 0:

$$ \nabla \cdot \mathbf{v} = 0, $$

so that

$$ \forall f \in A_{k+1} \quad \sum_{C \in \mathcal{C}} \int_C (\nabla^{\perp} f) \cdot \partial_t \mathbf{u} = 0. $$
This equation formally ensures the following advection equation on $\nabla \cdot \mathbf{u}$:

$$ \partial_t (\nabla \cdot \mathbf{u}) + \nabla \cdot ((\nabla \cdot \mathbf{u}) \mathbf{b}) = 0, \qquad (24) $$

which means that if $\nabla \cdot \mathbf{u}$ is initially equal to 0, it is always equal to 0. We are now interested in the discretization of (23). In Proposition 3, we proved that $(\nabla)^*$ could be preserved by the discontinuous Galerkin method provided the approximation space is $\mathbf{d}\mathbf{B}_k^{\text{curl}}$ and the numerical flux is tangential to the normal. This suggests to propose the following numerical scheme:

$$ \text{find } \mathbf{u} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \text{ such that for all } \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \\ \sum_{c \in \mathcal{C}} \int_c \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{c \in \mathcal{C}} \int_c (\mathbf{b} \cdot \mathbf{u}^{\perp}) \nabla^{\perp} \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v}] \cdot \tilde{\mathbf{G}} + \sum_{c \in \mathcal{C}} \int_c \mathbf{v} \cdot \mathbf{b} (-\nabla^{\star} \mathbf{u}) = 0, \qquad (25) $$

for which we can prove

**Proposition 4** (Equation on $-\nabla^*\mathbf{u}$). If $\mathbf{u}$ is solution of (25), and if the numerical flux is orthogonal to $\mathbf{n}_S$, then $-\nabla^*\mathbf{u}$ is solution of the continuous Galerkin discretization for the advection equation on the divergence (24)

$$ \text{Find } D \in \mathbb{A}_{k+1} \quad \forall f \in \mathbb{A}_{k+1} \quad \sum_{c \in \mathcal{C}} \int_c f \partial_t D - \sum_{c \in \mathcal{C}} \int_c D \mathbf{b} \cdot \nabla f = 0. \qquad (26) $$

*Proof.* We consider (25), which is tested with $\mathbf{v} = \nabla f$, where $f \in \mathbb{A}_{k+1}$. Then

$$ \nabla^{\perp} \cdot \mathbf{v} = \nabla^{\perp} \cdot (\nabla f) = 0. $$

As $\tilde{\mathbf{G}}$ is orthogonal to $\mathbf{n}_S$, it may be rewritten $\tilde{g}\mathbf{n}_S^{\perp}$. This leads to

$$ [\mathbf{v}] \cdot \tilde{\mathbf{G}} = [\mathbf{v}] \cdot \tilde{g}\mathbf{n}_S^{\perp} = [\mathbf{v} \cdot \mathbf{n}_S^{\perp}] \tilde{g}, $$

and if $f \in \mathbb{A}_{k+1}$, then $[\nabla f \cdot \mathbf{n}_S^{\perp}] = 0$. Therefore (25) tested with $\nabla f$ for $f \in \mathbb{A}_{k+1}$ gives

$$ \forall f \in \mathbb{A}_{k+1} \quad \sum_{c \in \mathcal{C}} \int_c \nabla f \cdot \partial_t \mathbf{u} + \sum_{c \in \mathcal{C}} \int_c \nabla f \cdot \mathbf{b} (-\nabla^* \mathbf{u}) = 0. \qquad (27) $$

The time derivative term may be transformed as (the commutation between $\nabla^*$ and $\partial_t$ can be proven as in the proof of Proposition 2):

$$ \begin{aligned} \sum_{c \in \mathcal{C}} \int_c \nabla f \cdot \partial_t \mathbf{u} &= \langle \nabla f | \partial_t \mathbf{u} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{curl}}} \\ &= \langle f | \nabla^* (\partial_t \mathbf{u}) \rangle_{\mathbb{A}_{k+1}} \\ &= \langle f | \partial_t \nabla^* (\mathbf{u}) \rangle_{\mathbb{A}_{k+1}} \\ \sum_{c \in \mathcal{C}} \int_c \nabla f \cdot \partial_t \mathbf{u} &= - \langle f | \partial_t (-\nabla^* (\mathbf{u})) \rangle_{\mathbb{A}_{k+1}} . \end{aligned} $$

Equation (27) becomes finally

$$ \forall f \in \mathbb{A}_{k+1} \quad \sum_{c \in \mathcal{C}} \int_c f \partial_t (-\nabla^* (\mathbf{u})) - \sum_{c \in \mathcal{C}} \int_c \nabla f \cdot \mathbf{b} (-\nabla^* \mathbf{u}) = 0, $$

which means that $(-\nabla^* \mathbf{u})$ is solution of (26).

□

Note that $-\nabla^* \mathbf{u}$ is solution of the numerical scheme (26), which is known to be unstable with explicit time stepping, as it leads to central-difference type approximations of differential operators [20]. However, it will be used in practical applications with initial condition $\nabla^* \mathbf{u}(t=0) = 0$, so that this instability should not raise any problem provided the initialization is carefully performed.
This last equation means that

$$ \forall f \in A_{k+1} \quad \langle \nabla^{\perp} f | \partial_t \mathbf{u} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = 0. $$

By definition of $(\nabla^{\perp})^*$, we get

$$ \forall f \in A_{k+1} \quad \langle f | (\nabla^{\perp})^* (\partial_t \mathbf{u}) \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = 0. $$

It remains to prove that $(\nabla^{\perp})^* (\partial_t \mathbf{u}) = \partial_t ((\nabla^{\perp})^* \mathbf{u})$. Denoting by $f$ an element of $A_{k+1}$, we have

$$ \begin{aligned} \langle \partial_t ((\nabla^{\perp})^* \mathbf{u}) | f \rangle_{A_{k+1}} &= \int_{\Omega} \partial_t ((\nabla^{\perp})^* \mathbf{u}) f \\ &= \partial_t \left( \int_{\Omega} ((\nabla^{\perp})^* \mathbf{u}) f \right) \\ &= \partial_t \left( \int_{\Omega} \mathbf{u} \cdot \nabla^{\perp} f \right) \\ &= \int_{\Omega} (\partial_t \mathbf{u}) \cdot \nabla^{\perp} f \\ &= \langle \partial_t \mathbf{u} | \nabla^{\perp} f \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} \\ \langle \partial_t ((\nabla^{\perp})^* \mathbf{u}) | f \rangle_{A_{k+1}} &= \langle (\nabla^{\perp})^* (\partial_t \mathbf{u}) | f \rangle_{A_{k+1}}, \end{aligned} $$

which ends the proof.

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

**Remark 2.** The Proposition 4 was proven in the semi-discrete case. It is however easy to prove the same in the discrete case. Still, the definition of the adjoint curl $(\nabla^{\perp})^*$ should take into account how the mass matrix is computed. For example, if the mass matrix is lumped, then the adjoint curl should be defined with the same lumped mass matrix.

**Remark 3** (Lax-Friedrich flux with normal diffusion). A widely used numerical flux is the Lax-Friedrich flux which reads

$$ \tilde{\mathbf{G}}(\overline{\mathbf{G}}_L, \mathbf{u}_L, \overline{\mathbf{G}}_R, \mathbf{u}_R, \mathbf{n}) = \frac{\overline{\mathbf{G}}_L \mathbf{n} + \overline{\mathbf{G}}_R \mathbf{n}}{2} + \frac{\lambda}{2}(\mathbf{u}_L - \mathbf{u}_R), \quad (13) $$

where $\lambda$ is the maximum absolute value of the eigenvalues of the system. Dealing with the system (11), the numerical flux (13) can be simplified as

$$ \tilde{\mathbf{G}}(g_L, \mathbf{u}_L, g_R, \mathbf{u}_R, \mathbf{n}) = \frac{g_L \mathbf{n} + g_R \mathbf{n}}{2} + \frac{\lambda}{2} (\mathbf{u}_L - \mathbf{u}_R). $$

The centered part of the flux is clearly parallel to the normal, and should not be modified because it ensures the consistency of the numerical scheme. The diffusive part may be decomposed into its normal and tangential part as

$$ (\mathbf{u}_L - \mathbf{u}_R) = \mathbf{n}\mathbf{n}^T (\mathbf{u}_L - \mathbf{u}_R) + (\mathbf{I}_d - \mathbf{n}\mathbf{n}^T) (\mathbf{u}_L - \mathbf{u}_R). $$

Therefore an easy way to ensure the hypothesis of Proposition 4 is to use the following Lax-Friedrich flux with purely normal diffusion

$$ \tilde{\mathbf{G}}(g_L, \mathbf{u}_L, g_R, \mathbf{u}_R, \mathbf{n}) = \frac{g_L \mathbf{n} + g_R \mathbf{n}}{2} + \frac{\lambda \mathbf{n}\mathbf{n}^T}{2} (\mathbf{u}_L - \mathbf{u}_R). \quad (14) $$

**Remark 4** (Link with the Hodge Laplacian). The diffusion induced by (14) is exactly the face component of the component $\delta^2 \circ d^1$ of the Hodge Laplacian that was found in (8). Formally, the numerical resolution of (11) with the numerical flux (14) may then be seen as

$$ \partial_t \mathbf{u} - \delta^2 \tilde{\mathbf{f}} = \lambda \delta^2 \circ d^1 \mathbf{u}, $$

where $\delta^2 \circ d^1$ is the face component of the operator $\delta^2 \circ d^1$, and $\tilde{\mathbf{f}}$ is the projection of $\mathbf{f}$ on the space $\mathcal{C}_k$. The conservation of the adjoint curl is then translated into $\partial_t (\delta^1 \mathbf{u}) = 0$.
**Proposition 5** (Exact preservation of $\nabla^*u = 0$). We consider equation (23). Suppose that $u(t=0) = u^0$ derives from a potential $f^0$: $u^0 = \nabla^\perp f^0$, and we denote by $f_h^0$ the projection of $f^0$ in $\mathbb{C}_k$ as defined in Definition 2. Then $u_h^0 := -(\nabla_{\mathcal{D}'}^\perp)^* f_h^0$ is such that $\nabla^*u_h^0 = 0$, and under hypothesis of Proposition 4, the constraint $\nabla^*u_h = 0$ is preserved for all time by the numerical scheme (25).

*Proof.* The beginning of the proof relies only on relations between the kernel and range of an operator and its adjoint. Note that we are working here in finite dimension, so that these relations do not include closures of range or kernel.

If $f_h^0 \in \mathbb{C}_k$, then $(\nabla_{\mathcal{D}'}^\perp)^* f_h^0 \in (\ker(\nabla_{\mathcal{D}'}^\perp))^{\perp}$. But as Range $\nabla \subset \ker(\nabla_{\mathcal{D}'}^\perp)$, we also have $(\nabla_{\mathcal{D}'}^\perp)^* f_h^0 \in (\text{Range} \nabla)^\perp$. As $(\text{Range} \nabla)^\perp = \ker \nabla^*$, we have $(\nabla_{\mathcal{D}'}^\perp)^* f_h^0 \in \ker \nabla^*$, and so $\nabla^*u_h^0 = 0$.

Then by Proposition 4, $\nabla^*u_h$ ensures the numerical scheme (26) which preserves $\nabla^*u_h = 0$ for all time, which ends the proof. $\square$

**Remark 6** (Divergence free initialization). We provide some details about how the initial condition is computed. Usually, the computation of the initial condition $u_h^0$ is computed by projection of the initial condition $u^0$ on the finite element space. This leads to the following system to solve

$$M u_h^0 = \mathbf{RHS},$$

where $M$ is the mass matrix of the finite element space of the velocity space: if $v$ denotes a basis of this finite element basis, the mass matrix is

$$M_{i,j} = \int_{\Omega} v_i \cdot v_j, \qquad (28)$$

and RHS is the right hand side, equal to

$$\mathbf{RHS}_i = \int_{\Omega} u^0 \cdot v_i. \qquad (29)$$

Suppose now that the initial condition derives from a potential vector $f$: $u^0 = \nabla^\perp f^0$. Then at the discrete level, we wish to have $u_h^0 = -(\nabla_{\mathcal{D}'}^\perp)^* f_h^0$, namely

$$\forall v \in \mathbf{d}\mathbb{B}_k^{\text{curl}} \quad \langle u_h^0 | v \rangle_{\mathbf{d}\mathbb{B}_k^{\text{curl}}} = - \langle (\nabla_{\mathcal{D}'}^\perp)^* f_h^0 | v \rangle_{\mathbf{d}\mathbb{B}_k^{\text{curl}}} = - \langle f_h^0 | \nabla_{\mathcal{D}'}^\perp \cdot v \rangle_{\mathbb{C}_k}.$$

This gives

$$\forall v \in \mathbf{d}\mathbb{B}_k^{\text{curl}} \quad \langle u_h^0 | v \rangle_{\mathbf{d}\mathbb{B}_k^{\text{curl}}} = - \sum_{c \in \mathcal{C}} \int_c f_h^0 \nabla^\perp \cdot v + \sum_{S \in \mathcal{S}} \int_S f_h^0 \llbracket v \cdot n_S^\perp \rrbracket.$$

Therefore, the matrix of the system to solve is still (28), but the right hand side is no more (29), but

$$\mathbf{RHS}_i = - \sum_{c \in \mathcal{C}} \int_c f_h^0 \nabla^\perp \cdot v_i + \sum_{S \in \mathcal{S}} \int_S f_h^0 \llbracket v_i \cdot n_S^\perp \rrbracket.$$

# 5 Numerical results

For all the numerical tests, the computational domain is the periodic square domain $[0, 1]^2$. We will consider three types of mesh: Cartesian, unstructured quadrangular and unstructured triangular. The meshes on which the conservation tests will be performed are represented in Figure 1.

The convergence tests will be performed on Cartesian meshes with uniform size of the mesh, with $h = 0.1$, $h = 0.05$, $h = 0.025$ and $h = 0.0125$. The convergence tests on triangular meshes will be performed on a set of unstructured meshes with a number of cells equal to 268, 1036, 4186
## 4.3 Discrete preservation of the irrotational component

In this section, we consider the case in which $\overline{G}$ is antisymmetric; in dimension 2, $\overline{G}$ may be written as

$$ \begin{pmatrix} 0 & -g \\ g & 0 \end{pmatrix}, $$

where $g$ is a scalar that may depend nonlinearly on the variables of the system. This means that $\mathbf{u}$ ensures the following conservation law

$$ \partial_t \mathbf{u} + \nabla^\perp g = 0. \qquad (15) $$

Then the discontinuous Galerkin discretization reads

$$ \text{Find } \mathbf{u} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \quad \forall \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \quad \sum_{C \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{C \in \mathcal{C}} \int_C g \nabla^\perp \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v}] \cdot \tilde{\mathbf{G}} = 0, \quad (16) $$

where $\tilde{\mathbf{G}}$ is the numerical flux.

**Proposition 5** (Adjoint conservation of the divergence). *Consider the numerical scheme (16). If $\tilde{\mathbf{G}}$ is orthogonal to $\mathbf{n}$, then $(\nabla)^* \mathbf{u}$ is conserved by the numerical scheme.*

The proof of Proposition 5 follows exactly the same lines as Proposition 4 and is not repeated.

**Remark 5** (Lax-Friedrich flux with purely tangential diffusion). The Lax-Friedrich numerical flux (13) for (15) may be simplified as

$$ \tilde{\mathbf{G}}(\mathbf{g}_L, \mathbf{u}_L, \mathbf{g}_R, \mathbf{u}_R, \mathbf{n}) = \frac{g_L \mathbf{n}^\perp + g_R \mathbf{n}^\perp}{2} + \frac{\lambda}{2} (\mathbf{u}_L - \mathbf{u}_R). $$

The centered part of the flux, which is clearly orthogonal to $\mathbf{n}$, ensures the consistency of the numerical scheme. Still relying on the decomposition of the diffusive part into its normal and tangential part, an easy way to ensure the hypothesis of Proposition 5 consists in using the following *Lax-Friedrich flux with purely tangential diffusion*

$$ \tilde{\mathbf{G}}(\mathbf{g}_L, \mathbf{u}_L, \mathbf{g}_R, \mathbf{u}_R, \mathbf{n}) = \frac{g_L \mathbf{n}^\perp + g_R \mathbf{n}^\perp}{2} + \frac{\lambda (I_d - \mathbf{n}\mathbf{n}^T)}{2} (\mathbf{u}_L - \mathbf{u}_R). \quad (17) $$

The diffusion of (17) matches with the face component of the $\delta^2 \circ d^1$ component of the Hodge Laplacian (7).

Last, we remark by translating back the equation into differential forms, the same equation as in Remark 4 holds.

# 5 Preservation of initially curl or divergence free fields for Lie derivative based advection equation

## 5.1 The direct and adjoint continuous equations

In this section, we address a slightly different problem, in which a vector field is seen as a proxy of the advection equation based on the *Lie derivative*. We first introduce an operator on the differential forms $\Lambda^k$ for $k \ge 1$ the contraction by a vector $\mathbf{b}$. It is defined as

$$ \begin{aligned} i_{\mathbf{b}} : & \Lambda^k \mapsto \Lambda^{k-1} \\ \mathbf{u} & \mapsto i_{\mathbf{b}} \mathbf{u}, \end{aligned} $$

with

$$ \forall \mathbf{v}_1, \dots, \mathbf{v}_{k-1} \in \mathbb{R}^2 \quad i_{\mathbf{b}} \mathbf{u} (\mathbf{v}_1, \dots, \mathbf{v}_{k-1}) := \mathbf{u} (\mathbf{b}, \mathbf{v}_1, \dots, \mathbf{v}_{k-1}). $$
<figure><img src="image_4.png" /><figcaption>Figure 1: The meshes on which the conservation of divergence or curl are performed. Top left: Cartesian mesh, top right: unstructured quadrangular mesh, bottom: unstructured triangular mesh.</figcaption></figure>
We also define the co-contraction $j_b$ as

$$j_b : \begin{array}{ccc} \Lambda^k & \mapsto & \Lambda^{k+1} \\ \mathbf{u} & \mapsto & j_b \mathbf{u}, \end{array}$$

such that $\star_{\{k+1\}} \circ j_b = (-1)^k i_b \circ \star_{\{k\}}$.

* If the choice (1) is done, then

$$\left\{ \begin{array}{lll} \forall f \in \Lambda^2 & \text{Proxy}(f) = f & \text{Proxy}(i_b f) = f b^{\perp} \\ \forall u \in \Lambda^1 & \text{Proxy}(u) = u & \text{Proxy}(i_b u) = u \cdot b, \\ \forall f \in \Lambda^0 & \text{Proxy}(f) = f & \text{Proxy}(j_b f) = f b, \\ \forall u \in \Lambda^1 & \text{Proxy}(u) = u & \text{Proxy}(j_b u) = -u^{\perp} \cdot b, \end{array} \right.$$

* If the choice (2) is done, then

$$\left\{ \begin{array}{lll} \forall f \in \Lambda^2 & \text{Proxy}^{\perp}(f) = f & \text{Proxy}^{\perp}(i_b f) = -f b \\ \forall u \in \Lambda^1 & \text{Proxy}^{\perp}(u) = u & \text{Proxy}^{\perp}(i_b u) = \det(b, u) \\ \forall f \in \Lambda^0 & \text{Proxy}^{\perp}(f) = f & \text{Proxy}^{\perp}(j_b f) = f b^{\perp}, \\ \forall u \in \Lambda^1 & \text{Proxy}^{\perp}(u) = u & \text{Proxy}^{\perp}(j_b u) = -u \cdot b, \end{array} \right.$$

The Lie derivative $L_b$ is defined as

$$L_b : \begin{array}{ccc} C^{\infty} \Lambda^k & \mapsto & C^{\infty} \Lambda^k \\ \mathbf{u} & \mapsto & L_b \mathbf{u} = i_b \circ d^k \mathbf{u} + d^{k-1} \circ i_b \mathbf{u}. \end{array}$$

We consider the following advection equation discussed in [36, 38, 37, 53]

$$\partial_t u + L_b u = 0. \qquad (18)$$

It is easy to check that the Lie derivative commutes with the exterior derivative, so that if $u$ follows (18), then

$$\partial_t (d^k u) + L_b (d^k u) = 0.$$

Especially, if $d^k u(t=0) = 0$, then for all time, $d^k u = 0$.

Considering now the Lie co-derivative $\mathcal{L}_b$, defined as

$$\mathcal{L}_b : \begin{array}{ccc} C^{\infty} \Lambda^k & \mapsto & C^{\infty} \Lambda^k \\ \mathbf{u} & \mapsto & \mathcal{L}_b \mathbf{u} = j_b \circ \delta^k \mathbf{u} + \delta^{k+1} \circ j_b \mathbf{u}. \end{array}$$

the adjoint Lie advection problem

$$\partial_t u - \mathcal{L}_b u = 0, \qquad (19)$$

can be considered. Then similar properties as for (18) can be proven for $\delta^k u$ instead of $d^k u$.

Written with proxies on $C^{\infty}\Lambda^1$, equations (18) and (19) lead to two equations

* The first one is

$$\partial_t u + \nabla (b \cdot u) + b^{\perp} \nabla^{\perp} \cdot u = 0,$$

and is (18) with the choice of proxy (1), and is also (19) with the choice of proxy (2). $\nabla^{\perp} \cdot u$ follows the conservative equation

$$\partial_t (\nabla^{\perp} \cdot u) + \nabla \cdot ((\nabla^{\perp} \cdot u) b) = 0.$$

If $\nabla^{\perp} \cdot u(t=0) = 0$, then $\nabla^{\perp} \cdot u = 0$ for any $t > 0$. It can be encountered in the equation on $w$ in the two-phase flows model of [60], which was for example numerically addressed in [59].

* The second one is the two dimensional induction equation

$$\partial_t u + \nabla^{\perp} (\det(b, u)) + b \nabla \cdot u = 0, \qquad (20)$$

obtained with (19) with the choice of proxy (1), and is also (18) with the choice of proxy (2). Then $\nabla \cdot u$ follows the conservative equation

$$\partial_t (\nabla \cdot u) + \nabla \cdot ((\nabla \cdot u) b) = 0. \qquad (21)$$

If $\nabla \cdot u(t=0) = 0$, then $\nabla \cdot u = 0$ for any $t > 0$.
and 16682, and minimal size (defined as the square root of the surface of the smallest element) approximately equal to 0.0844, 0.0419, 0.0223 and 0.0109.

For all the test cases, the time stepping is done with classical Strong-Stability Preserving explicit time schemes [38]. If k is the polynomial order of the space discretization, the order of the time integration scheme is k + 1. The CFL number is equal to 0.5 for k = 0, 0.33 for k = 1 and 0.2 for k = 2.

Concerning the optimal order of convergence we may get, we make the following remark

**Remark 7 (Optimal order of convergence).** It has been known for a while that the convergence order of the discontinuous Galerkin method for advection problems is of order k + 1/2 [51, 72] in general. On Cartesian meshes, the order of convergence can be improved to k + 1 by using tensor products of one dimensional basis [56]. Here, we are working with basis for vectors that are not tensor products of one dimensional basis. Therefore, the optimal order we can expect is only k + 1/2, even though we will frequently observe k + 1.

## 5.1 Discrete conservation of the divergence: Maxwell system

In this section, we are interested in the two-dimensional Maxwell system in the vacuum, which reads

$$ \left\{ \begin{array}{l} \partial_t b + \nabla^{\perp} \cdot \mathbf{e} = 0 \\ \partial_t \mathbf{e} + c^2 \nabla^{\perp} b = 0, \end{array} \right. \qquad (30) $$

where e is the two-dimensional electric field and b is the z-component of the magnetic field.

### 5.1.1 Conservation of the divergence of a stationary solution

For this first test case, the following initial condition is imposed

$$ \left\{ \begin{array}{l} b(\mathbf{x}) = 0 \\ \mathbf{e}_x(\mathbf{x}) = \bar{x} e^{-\bar{r}^2/2} \\ \mathbf{e}_y(\mathbf{x}) = \bar{y} e^{-\bar{r}^2/2}, \end{array} \right. $$

with $\bar{x} = \frac{x - x_c}{r_0}$, $\bar{y} = \frac{y - y_c}{r_0}$, $r^2 = (x - x_c)^2 + (y - y_c)^2$, and $\bar{r} = \frac{r}{r_0}$. This initial condition is such that $\nabla^{\perp} \cdot \mathbf{e} = 0$, and so is clearly a stationary solution of (30). This solution was built by considering the potential $e^{-\bar{r}^2/2}$ and by taking its gradient. The numerical parameters are $r_0 = 0.15$, $x_c = y_c = 0.5$, and $c^2 = 1$. The computation is led until $t = 3$ on the different meshes represented in Figure 1, and The L² difference between $\nabla^* \mathbf{e}$ and its initial value is computed along the time. Two configurations of finite element spaces for approximating e are used:

* the finite element spaces $\mathbf{d}\mathbf{B}_k^{\text{curl}}$,

* the classical finite element spaces for discontinuous Galerkin methods obtained by tensorization of the scalar basis (note however that in the triangular case, these two approximation spaces match),

and with two different numerical flux:

* the Lax-Friedrich flux,

* the Lax-Friedrich flux with purely tangential diffusion (22), which matches in the case of the two-dimensional Maxwell system (30) with the Godunov' flux.

Numerical results are represented in Figure 2 for degree k = 0, 1, 2, and show that the only combination that preserves correctly the divergence is the one with the finite element space $\mathbf{d}\mathbf{B}_k^{\text{curl}}$ and with the Godunov' numerical flux.
## 5.2 Discretization

We are interested in the discretization of (20). This equation can be seen either as the equation (18) with the proxy (2), or as the equation (19) with the proxy (1). In Proposition 5, we proved that $(\nabla)^*$ could be preserved by the discontinuous Galerkin method provided the approximation space is $\mathbf{d}\mathbf{B}_k^{\text{curl}}$ and the numerical flux is tangential to the normal. This suggests to see the equation (20) as the equation (19) with the proxy (1). This leads to the following numerical scheme:

$$ \text{find } \mathbf{u} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \text{ such that for all } \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \\ \sum_{c \in \mathcal{C}} \int_C \mathbf{v} \cdot \partial_t \mathbf{u} - \sum_{c \in \mathcal{C}} \int_C \det(\mathbf{b}, \mathbf{u}) \nabla^\perp \cdot \mathbf{v} + \sum_{S \in \mathcal{S}} \int_S [\mathbf{v}] \cdot \tilde{\mathbf{G}} + \sum_{c \in \mathcal{C}} \int_C \mathbf{v} \cdot \mathbf{b} (-\nabla^* \mathbf{u}) = 0, \quad (22) $$

for which we can prove

**Proposition 6** (Equation on $-\nabla^*\mathbf{u}$). If $\mathbf{u}$ is solution of (22), and if the numerical flux is orthogonal to $\mathbf{n}$, then $-\nabla^*\mathbf{u}$ is solution of the continuous Galerkin discretization for the advection equation on the divergence (21)

$$ \text{Find } D \in \mathbb{A}_{k+1} \quad \forall f \in \mathbb{A}_{k+1} \quad \sum_{c \in \mathcal{C}} \int_C f \partial_t D - \sum_{c \in \mathcal{C}} \int_C D \mathbf{b} \cdot \nabla f = 0. \quad (23) $$

*Proof.* We consider (22), which is tested with $\mathbf{v} = \nabla f$, where $f \in \mathbb{A}_{k+1}$. Then

$$ \nabla^\perp \cdot \mathbf{v} = \nabla^\perp \cdot (\nabla f) = 0. $$

As $\tilde{\mathbf{G}}$ is orthogonal to $\mathbf{n}$, it may be rewritten $\tilde{g}\mathbf{n}^\perp$. This leads to

$$ [\mathbf{v}] \cdot \tilde{\mathbf{G}} = [\mathbf{v}] \cdot \tilde{g}\mathbf{n}^\perp = [\mathbf{v} \cdot \mathbf{n}^\perp] \tilde{g}, $$

and if $f \in \mathbb{A}_{k+1}$, then $[\nabla f \cdot \mathbf{n}^\perp] = 0$. Therefore (22) tested with $\nabla f$ for $f \in \mathbb{A}_{k+1}$ gives

$$ \forall f \in \mathbb{A}_{k+1} \quad \sum_{c \in \mathcal{C}} \int_C \nabla f \cdot \partial_t \mathbf{u} + \sum_{c \in \mathcal{C}} \int_C \nabla f \cdot \mathbf{b} (-\nabla^* \mathbf{u}) = 0. \quad (24) $$

The time derivative term may be transformed as (the commutation between $\nabla^*$ and $\partial_t$ can be proven as in the proof of Proposition 4):

$$ \sum_{c \in \mathcal{C}} \int_C \nabla f \cdot \partial_t \mathbf{u} = \langle \nabla f | \partial_t \mathbf{u} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{curl}}} \\ = \langle f | \nabla^* (\partial_t \mathbf{u}) \rangle_{\mathbb{A}_{k+1}} \\ = \langle f | \partial_t \nabla^* (\mathbf{u}) \rangle_{\mathbb{A}_{k+1}} \\ \sum_{c \in \mathcal{C}} \int_C \nabla f \cdot \partial_t \mathbf{u} = - \langle f | \partial_t (-\nabla^* (\mathbf{u})) \rangle_{\mathbb{A}_{k+1}}. $$

Equation (24) becomes finally

$$ \forall f \in \mathbb{A}_{k+1} \quad \sum_{c \in \mathcal{C}} \int_C f \partial_t (-\nabla^* (\mathbf{u})) - \sum_{c \in \mathcal{C}} \int_C \nabla f \cdot \mathbf{b} (-\nabla^* \mathbf{u}) = 0, $$

which means that $(-\nabla^* \mathbf{u})$ is solution of (23). $\square$

Note that $-\nabla^* \mathbf{u}$ is solution of the numerical scheme (23), which is known to be unstable. However, it will be used in practical applications with initial condition $\nabla^* \mathbf{u}(t=0) = 0$, so that this instability should not raise any problem provided the initialization is carefully performed.

**Proposition 7** (Exact preservation of $\nabla^* \mathbf{u} = 0$). We consider equation (20). Suppose that $\mathbf{u}(t=0) = \mathbf{u}^0$ derives from a potential $f^0$: $\mathbf{u}^0 = \nabla^\perp f^0$, and we denote by $f_h^0$ the projection of $f^0$ in $\mathbb{C}_k$. Then $\mathbf{u}_h^0 := -(\nabla^\perp \tilde{\mathbf{u}})^* f_h^0$ is such that $\nabla^* \mathbf{u}_h^0 = 0$, and under hypothesis of Proposition 6, the constraint $\nabla^* \mathbf{u}_h = 0$ is preserved for all time by the numerical scheme (22).
<figure><img src="image_5.png" /><figcaption>Figure 2: Plot of $||\nabla^*e - \nabla^*e_0||_2$ with respect to time for different degree, approximation space, type of meshes and numerical flux. In the left column, Cartesian and unstructured quadrangular meshes are considered, with the $d\mathcal{Q}_k$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. In the right column, Cartesian, unstructured quadrangular and triangular meshes are considered, with the $dB_k^{curl}$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. On triangular meshes, the space $dB_k^{curl}$ and $dP_k$ are the same, and the computations on these spaces are represented on the right column. Note that the y scaling of the top right figure ($dB_k^{curl}$ with Godunov's scheme) is $10^{-12}$, whereas it is of the order of 1 for the other plots.</figcaption></figure>
*Proof.* The beginning of the proof relies only on relations between the kernel and range of an operator and its adjoint. Note that we are working here in finite dimension, so that these relations do not include closures of range or kernel.

If $f \in \mathbb{C}_k$, then $(\nabla_{\mathcal{D}'}^{\perp})^* f \in (\ker (\nabla_{\mathcal{D}'}^{\perp}))^{\perp}$. But as $\text{Range} \nabla \subset \ker (\nabla_{\mathcal{D}'}^{\perp})$, we also have $(\nabla_{\mathcal{D}'}^{\perp})^* f \in (\text{Range} \nabla)^{\perp}$. As $(\text{Range} \nabla)^{\perp} = \ker \nabla^*$, we have $(\nabla_{\mathcal{D}'}^{\perp})^* f \in \ker \nabla^*$, and so $\nabla^* \mathbf{u}_h^0 = 0$.

Then by Proposition 6, $\nabla^* \mathbf{u}_h$ ensures the numerical scheme (23) which preserves $\nabla^* \mathbf{u}_h = 0$ for all time, which ends the proof. □

**Remark 6 (Divergence free initialization).** We provide some details about how the initial condition is computed. Usually, the computation of the initial condition $\mathbf{u}_0^h$ is computed by projection of the initial condition $\mathbf{u}_0$ on the finite element space. This leads to the following system to solve

$$M \mathbf{u}_0^h = \mathbf{RHS},$$

where $M$ is the mass matrix of the finite element space of the velocity space: if $\mathbf{v}$ denotes a basis of this finite element basis, the mass matrix is

$$M_{i,j} = \int_{\Omega} \mathbf{v}_i \cdot \mathbf{v}_j, \qquad (25)$$

and RHS is the right hand side, equal to

$$\mathbf{RHS}_i = \int_{\Omega} \mathbf{u}_0 \cdot \mathbf{v}_i. \qquad (26)$$

Suppose now that the initial condition derives from a potential vector $f$: $\mathbf{u}_0 = \nabla_{\mathcal{D}'}^{\perp} f$. Then at the discrete level, we wish to have $\mathbf{u}_0^h = -(\nabla_{\mathcal{D}'}^{\perp})^* f$, namely

$$\forall \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \quad \langle \mathbf{u}_0^h | \mathbf{v} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = - \langle (\nabla_{\mathcal{D}'}^{\perp})^* f | \mathbf{v} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = - \langle f | \nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{v} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}}. $$

This gives

$$\forall \mathbf{v} \in \mathbf{d}\mathbf{B}_k^{\text{curl}} \quad \langle \mathbf{u}_0^h | \mathbf{v} \rangle_{\mathbf{d}\mathbf{B}_k^{\text{div}}} = - \sum_{c \in \mathcal{C}_c} \int_c f \nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{v} + \sum_{S \in \mathcal{S}_S} \int_S f \llbracket \mathbf{v} \cdot \mathbf{n}^{\perp} \rrbracket. $$

Therefore, the matrix of the system to solve is still (25), but the right hand side is no more (26), but

$$\mathbf{RHS}_i = - \sum_{c \in \mathcal{C}_c} \int_c f \nabla_{\mathcal{D}'}^{\perp} \cdot \mathbf{v}_i + \sum_{S \in \mathcal{S}_S} \int_S f \llbracket \mathbf{v}_i \cdot \mathbf{n}^{\perp} \rrbracket. $$

# 6 Numerical results

For all the numerical tests, the computational domain is the periodic square domain $[0, 1]^2$. We will consider three types of mesh: Cartesian, unstructured quadrangular and unstructured triangular. The meshes on which the conservation tests will be performed are represented in Figure 3.

The convergence tests will be performed on Cartesian meshes with uniform size of the mesh, with $h = 0.1$, $h = 0.05$, $h = 0.025$ and $h = 0.0125$. The convergence tests on triangular meshes will be performed on a set of unstructured meshes with a number of cells equal to 268, 1036, 4186 and 16682, and minimal size (defined as the square root of the surface of the smallest element) approximately equal to 0.0844, 0.0419, 0.0223 and 0.0109.

For all the test cases, the time stepping is done with classical Strong-Stability Preserving explicit time schemes [31]. If $k$ is the polynomial order of the space discretization, the order of the time integration scheme is $k+1$. The CFL number is equal to 0.5 for $k = 0$, 0.33 for $k = 1$ and 0.2 for $k = 2$.
### 5.1.2 Convergence test

The numerical test is taken from [65, Section 5.1]. The initial condition is

$$ egin{cases} b(\mathbf{x}) = \frac{\omega}{c^2} \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x) \\ \mathbf{e}_x(\mathbf{x}) = -k_{\perp}\pi \sin(k_{\perp}\pi y) \cos(k_{\parallel}\pi x) \\ \mathbf{e}_y(\mathbf{x}) = k_{\parallel}\pi \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x) \end{cases} \tag{31} $$

and the exact solution is

$$ egin{cases} b(\mathbf{x}, t) = \frac{\omega}{c^2} \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x - \omega t) \\ \mathbf{e}_x(\mathbf{x}, t) = -k_{\perp}\pi \sin(k_{\perp}\pi y) \cos(k_{\parallel}\pi x - \omega t) \\ \mathbf{e}_y(\mathbf{x}, t) = k_{\parallel}\pi \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x - \omega t) \end{cases} $$

The longitudinal and transverse wave numbers $k_{\parallel}$ and $k_{\perp}$ are data of the test case. The frequency $\omega$ is such that

$$ k_{\parallel}^2 + k_{\perp}^2 = \frac{\omega^2}{\pi^2 c^2}. $$

The numerical parameters are $c = 1$, and $k_{\parallel} = k_{\perp} = 2$.

A second test case is performed, in which a stationary non divergence free solution is added to the initial solution (31) (it is similar to the stationary solution of subsubsection 5.1.1, but regular), equal to

$$ egin{cases} b(\mathbf{x}) = 0 \\ \mathbf{e}_x(\mathbf{x}) = 2K_0 \alpha \bar{x} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2} \\ \mathbf{e}_y(\mathbf{x}) = 2K_0 \alpha \bar{y} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2}, \end{cases} \tag{32} $$

if $r < r_0$, and 0 otherwise. The numerical parameters are $K_0 = 100$, $r_0 = 0.35$, $x_c = y_c = 0.5$, $\alpha = 4$.

The convergence curves obtained on the two test cases on each variable is plotted in Figure 3. The errors obtained and rate of convergence for Cartesian meshes on the variable $e_x$ are gathered in Table 1 for the case with the initial condition (31) and in Table 2 where (32) was added to the initial condition (31). With the single initial condition (31), the order of convergence obtained is not very sensitive to the change of approximation basis or of the change of numerical flux. This is probably due to the fact that (31) is already divergence free, and so the numerical divergence is already low when computing the initial condition, and so the benefit of preserving exactly the divergence by using the space $dB_{k}^{\text{curl}}$ with Godunov flux is low.

The same plots are performed on the series of triangular meshes, shown in Figure 4, and the errors obtained for the variable $e_x$ are shown in Table 3 and Table 4. The same observation as for Cartesian meshes hold.

## 5.2 Discrete conservation of a curl: Wave system

In this section, the wave system

$$ egin{cases} \partial_t p + \nabla \cdot \mathbf{u} = 0 \\ \partial_t \mathbf{u} + c^2 \nabla p = 0, \end{cases} \tag{33} $$

which couples the pressure $p$ and the velocity $\mathbf{u}$ is considered. The wave velocity $c$ is a parameter of the system, and will be always equal to 1 in the numerical applications.

### 5.2.1 Conservation of the curl of a stationary solution

For this first test case, the following initial condition is imposed

$$ egin{cases} p(\mathbf{x}) = 0 \\ \mathbf{u}_x(\mathbf{x}) = -\bar{y} e^{-\bar{r}^2/2} \\ \mathbf{u}_y(\mathbf{x}) = \bar{x} e^{-\bar{r}^2/2}, \end{cases} $$
<figure><img src="image_8.png" /><figcaption>Figure 3: The meshes on which the conservation of divergence or curl are performed. Top left: Cartesian mesh, top right: unstructured quadrangular mesh, bottom: unstructured triangular mesh.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3: Error obtained on the test case described in subsection 5.1.2 with initial condition (31) on the left, and with the initial condition (32) added to (31) on the right, on a series of Cartesian meshes. For each of the test cases, the error obtained for the variables $b$ (top row), $e_x$ (middle row) and $e_y$ (bottom row) is shown for different approximation spaces and for the Lax-Friedrich and Godunov flux.</figcaption></figure>
## 6.1 Discrete conservation of the divergence: Maxwell system

In this section, we are interested in the two-dimensional Maxwell system in the vacuum, which reads

$$ \left\{ \begin{array}{l} \partial_t b + \nabla^{\perp} \cdot \mathbf{e} = 0 \\ \partial_t \mathbf{e} + c^2 \nabla^{\perp} b = 0, \end{array} \right. \qquad (27) $$

where $\mathbf{e}$ is the two-dimensional electric field and $b$ is the $z$-component of the magnetic field.

### 6.1.1 Conservation of the divergence of a stationary solution

For this first test case, the following initial condition is imposed

$$ \left\{ \begin{array}{l} b(\mathbf{x}) = 0 \\ \mathbf{e}_x(\mathbf{x}) = \bar{x} e^{-\bar{r}^2/2} \\ \mathbf{e}_y(\mathbf{x}) = \bar{y} e^{-\bar{r}^2/2}, \end{array} \right. $$

with $\bar{x} = \frac{x - x_c}{r_0}$, $\bar{y} = \frac{y - y_c}{r_0}$, $r^2 = (x - x_c)^2 + (y - y_c)^2$, and $\bar{r} = \frac{r}{r_0}$. This initial condition is such that $\nabla^{\perp} \cdot \mathbf{e} = 0$, and so is clearly a stationary solution of (27). This solution was built by considering the potential $e^{-\bar{r}^2/2}$ and by taking its gradient. The numerical parameters are $r_0 = 0.15$, $x_c = y_c = 0.5$, and $c^2 = 1$. The computation is led until $t = 3$ on the different meshes represented in Figure 3, and The $L^2$ difference between $\nabla^* \mathbf{e}$ and its initial value is computed along the time. Two configurations of finite element spaces for approximating $\mathbf{e}$ are used:

* the finite element spaces $\mathbf{dB}_k^{\text{curl}}$,

* the classical finite element spaces for discontinuous Galerkin methods obtained by tensorization of the scalar basis (note however that in the triangular case, these two approximation spaces match),

and with two different numerical flux:

* the Lax-Friedrich flux,

* the Lax-Friedrich flux with purely tangential diffusion (17), which matches in our case with the Godunov' flux.

Numerical results are represented in Figure 4 for degree $k = 0, 1, 2$, and show that the only combination that preserves correctly the divergence is the one with the finite element space $\mathbf{dB}_k^{\text{curl}}$ and with the Godunov' numerical flux.

### 6.1.2 Convergence test

The numerical test is taken from [49, Section 5.1]. The initial condition is

$$ \left\{ \begin{array}{l} b(\mathbf{x}) = \frac{\omega}{c^2} \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x) \\ \mathbf{e}_x(\mathbf{x}) = -k_{\perp}\pi \sin(k_{\perp}\pi y) \cos(k_{\parallel}\pi x) \\ \mathbf{e}_y(\mathbf{x}) = k_{\parallel}\pi \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x) \end{array} \right. \qquad (28) $$

and the exact solution is

$$ \left\{ \begin{array}{l} b(\mathbf{x}, t) = \frac{\omega}{c^2} \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x - \omega t) \\ \mathbf{e}_x(\mathbf{x}, t) = -k_{\perp}\pi \sin(k_{\perp}\pi y) \cos(k_{\parallel}\pi x - \omega t) \\ \mathbf{e}_y(\mathbf{x}, t) = k_{\parallel}\pi \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x - \omega t) \end{array} \right. $$

The longitudinal and transverse wave numbers $k_{\parallel}$ and $k_{\perp}$ are data of the test case. The frequency $\omega$ is such that

$$ k_{\parallel}^2 + k_{\perp}^2 = \frac{\omega^2}{\pi^2 c^2}. $$
<table><caption>Table 1: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 5.1.2 with initial condition (31), on a series of quadrangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>curl</sup><sub>0</sub></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>curl</sup><sub>0</sub></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.76e+00</td><td></td><td>2.76e+00</td><td></td><td>2.99e+00</td><td></td><td>2.99e+00</td><td></td></tr><tr><td>0.05</td><td>2.02e+00</td><td>0.45</td><td>2.02e+00</td><td>0.45</td><td>2.44e+00</td><td>0.30</td><td>2.44e+00</td><td>0.30</td></tr><tr><td>0.025</td><td>1.26e+00</td><td>0.68</td><td>1.26e+00</td><td>0.68</td><td>1.65e+00</td><td>0.56</td><td>1.65e+00</td><td>0.56</td></tr><tr><td>0.0125</td><td>7.08e-01</td><td>0.83</td><td>7.08e-01</td><td>0.83</td><td>9.76e-01</td><td>0.76</td><td>9.76e-01</td><td>0.76</td></tr></tbody><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>curl</sup><sub>1</sub></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>curl</sup><sub>1</sub></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>1.84e-01</td><td></td><td>1.54e-01</td><td></td><td>1.75e-01</td><td></td><td>1.54e-01</td><td></td></tr><tr><td>0.05</td><td>4.23e-02</td><td>2.12</td><td>3.57e-02</td><td>2.11</td><td>3.88e-02</td><td>2.17</td><td>3.57e-02</td><td>2.11</td></tr><tr><td>0.025</td><td>1.04e-02</td><td>2.03</td><td>8.72e-03</td><td>2.03</td><td>9.28e-03</td><td>2.06</td><td>8.72e-03</td><td>2.03</td></tr><tr><td>0.0125</td><td>2.60e-03</td><td>2.00</td><td>2.17e-03</td><td>2.01</td><td>2.29e-03</td><td>2.02</td><td>2.17e-03</td><td>2.01</td></tr></tbody><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>curl</sup><sub>2</sub></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>curl</sup><sub>2</sub></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>6.41e-03</td><td></td><td>4.90e-03</td><td></td><td>7.17e-03</td><td></td><td>5.98e-03</td><td></td></tr><tr><td>0.05</td><td>8.16e-04</td><td>2.97</td><td>6.08e-04</td><td>3.01</td><td>9.32e-04</td><td>2.94</td><td>7.61e-04</td><td>2.97</td></tr><tr><td>0.025</td><td>1.03e-04</td><td>2.98</td><td>7.59e-05</td><td>3.00</td><td>1.19e-04</td><td>2.97</td><td>9.57e-05</td><td>2.99</td></tr><tr><td>0.0125</td><td>1.30e-05</td><td>2.99</td><td>9.48e-06</td><td>3.00</td><td>1.49e-05</td><td>2.99</td><td>1.20e-05</td><td>3.00</td></tr></tbody></table>

<table><caption>Table 2: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 5.1.2 with initial condition (31), on a series of Cartesian meshes. Results show a high benefit in using the space $dB_{k}^{curl}$ with the Godunov flux, which is the only one to always reach the optimal order.</caption><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>curl</sup><sub>0</sub></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>curl</sup><sub>0</sub></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.85e+00</td><td></td><td>2.78e+00</td><td></td><td>3.12e+00</td><td></td><td>3.11e+00</td><td></td></tr><tr><td>0.05</td><td>2.10e+00</td><td>0.44</td><td>2.03e+00</td><td>0.46</td><td>2.58e+00</td><td>0.27</td><td>2.56e+00</td><td>0.28</td></tr><tr><td>0.025</td><td>1.33e+00</td><td>0.66</td><td>1.26e+00</td><td>0.69</td><td>1.81e+00</td><td>0.51</td><td>1.78e+00</td><td>0.52</td></tr><tr><td>0.0125</td><td>7.59e-01</td><td>0.81</td><td>7.09e-01</td><td>0.83</td><td>1.14e+00</td><td>0.67</td><td>1.10e+00</td><td>0.70</td></tr></tbody><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>curl</sup><sub>1</sub></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>curl</sup><sub>1</sub></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.45e-01</td><td></td><td>1.71e-01</td><td></td><td>2.52e-01</td><td></td><td>1.98e-01</td><td></td></tr><tr><td>0.05</td><td>7.29e-02</td><td>1.75</td><td>3.94e-02</td><td>2.12</td><td>5.84e-02</td><td>2.11</td><td>4.11e-02</td><td>2.26</td></tr><tr><td>0.025</td><td>2.28e-02</td><td>1.68</td><td>9.70e-03</td><td>2.02</td><td>1.31e-02</td><td>2.15</td><td>9.83e-03</td><td>2.07</td></tr><tr><td>0.0125</td><td>6.96e-03</td><td>1.71</td><td>2.42e-03</td><td>2.00</td><td>2.87e-03</td><td>2.20</td><td>2.43e-03</td><td>2.02</td></tr></tbody><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>curl</sup><sub>2</sub></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>curl</sup><sub>2</sub></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>3.15e-02</td><td></td><td>1.12e-02</td><td></td><td>2.94e-02</td><td></td><td>1.38e-02</td><td></td></tr><tr><td>0.05</td><td>6.47e-03</td><td>2.29</td><td>2.07e-03</td><td>2.43</td><td>6.48e-03</td><td>2.18</td><td>2.72e-03</td><td>2.35</td></tr><tr><td>0.025</td><td>1.15e-03</td><td>2.50</td><td>2.53e-04</td><td>3.04</td><td>1.14e-03</td><td>2.50</td><td>3.96e-04</td><td>2.78</td></tr><tr><td>0.0125</td><td>1.93e-04</td><td>2.57</td><td>3.17e-05</td><td>3.00</td><td>1.93e-04</td><td>2.57</td><td>5.41e-05</td><td>2.87</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Figure 4: Plot of $||∇^*e − ∇^*e_0||_2$ with respect to time for different degree, approximation space, type of meshes and numerical flux. In the left column, Cartesian and unstructured quadrangular meshes are considered, with the $dQ_k$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. In the right column, Cartesian, unstructured quadrangular and triangular meshes are considered, with the $dB_k^{curl}$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. On triangular meshes, the space $dB_k^{curl}$ and $dP_k$ are the same, and the computations on these spaces are represented on the right column. Note that the y scaling of the top right figure ($dB_k^{curl}$ with Godunov's scheme) is $10^{-12}$, whereas it is of the order of 1 for the other plots.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Error obtained on the test case described in subsection 5.1.2 with initial condition (31) on the left, and with the initial condition (32) added to (31) on the right, on a series of unstructured triangular meshes. For each of the test cases, the error obtained for the variables $b$ (top row), $e_x$ (middle row) and $e_y$ (bottom row) is shown for different degrees, and for the Lax-Friedrich and Godunov numerical flux.</figcaption></figure>
<table><caption>Table 2: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 6.1.2 with initial condition (28), on a series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sub>0</sub><sup>curl</sup></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sub>0</sub><sup>curl</sup></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.76e+00</td><td></td><td>2.76e+00</td><td></td><td>2.99e+00</td><td></td><td>2.99e+00</td><td></td></tr><tr><td>0.05</td><td>2.02e+00</td><td>0.45</td><td>2.02e+00</td><td>0.45</td><td>2.44e+00</td><td>0.30</td><td>2.44e+00</td><td>0.30</td></tr><tr><td>0.025</td><td>1.26e+00</td><td>0.68</td><td>1.26e+00</td><td>0.68</td><td>1.65e+00</td><td>0.56</td><td>1.65e+00</td><td>0.56</td></tr><tr><td>0.0125</td><td>7.08e-01</td><td>0.83</td><td>7.08e-01</td><td>0.83</td><td>9.76e-01</td><td>0.76</td><td>9.76e-01</td><td>0.76</td></tr></tbody><tbody><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sub>1</sub><sup>curl</sup></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sub>1</sub><sup>curl</sup></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>1.84e-01</td><td></td><td>1.54e-01</td><td></td><td>1.75e-01</td><td></td><td>1.54e-01</td><td></td></tr><tr><td>0.05</td><td>4.23e-02</td><td>2.12</td><td>3.57e-02</td><td>2.11</td><td>3.88e-02</td><td>2.17</td><td>3.57e-02</td><td>2.11</td></tr><tr><td>0.025</td><td>1.04e-02</td><td>2.03</td><td>8.72e-03</td><td>2.03</td><td>9.28e-03</td><td>2.06</td><td>8.72e-03</td><td>2.03</td></tr><tr><td>0.0125</td><td>2.60e-03</td><td>2.00</td><td>2.17e-03</td><td>2.01</td><td>2.29e-03</td><td>2.02</td><td>2.17e-03</td><td>2.01</td></tr></tbody><tbody><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sub>2</sub><sup>curl</sup></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sub>2</sub><sup>curl</sup></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>6.41e-03</td><td></td><td>4.90e-03</td><td></td><td>7.17e-03</td><td></td><td>5.98e-03</td><td></td></tr><tr><td>0.05</td><td>8.16e-04</td><td>2.97</td><td>6.08e-04</td><td>3.01</td><td>9.32e-04</td><td>2.94</td><td>7.61e-04</td><td>2.97</td></tr><tr><td>0.025</td><td>1.03e-04</td><td>2.98</td><td>7.59e-05</td><td>3.00</td><td>1.19e-04</td><td>2.97</td><td>9.57e-05</td><td>2.99</td></tr><tr><td>0.0125</td><td>1.30e-05</td><td>2.99</td><td>9.48e-06</td><td>3.00</td><td>1.49e-05</td><td>2.99</td><td>1.20e-05</td><td>3.00</td></tr></tbody></table>

The numerical parameters are $c = 1$, and $k_\parallel = k_\perp = 2$.

A second test case is performed, in which a stationary non divergence free solution is added to the initial solution (28) (it is similar to the stationary solution of subsection 6.1.1, but regular), equal to

$$ \left\{ \begin{array}{l} b(\mathbf{x}) = 0 \\ \mathbf{e}_x(\mathbf{x}) = 2K_0 \alpha \bar{x} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2} \\ \mathbf{e}_y(\mathbf{x}) = 2K_0 \alpha \bar{y} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2}, \end{array} \right. \qquad (29) $$

if $r < r_0$, and 0 otherwise. The numerical parameters are $K_0 = 100$, $r_0 = 0.35$, $x_c = y_c = 0.5$, $\alpha = 4$.

The convergence curves obtained on the two test cases on each variable is plotted in Figure 5. The errors obtained and rate of convergence for Cartesian meshes on the variable $e_x$ are gathered in Table 2 for the case with the initial condition (28) and in Table 2 where (29) was added to the initial condition (28). With the single initial condition (28), the order of convergence obtained is very few sensitive to the change of approximation basis or of the change of numerical flux. This is probably due to the fact that (28) is already divergence free, and so the numerical divergence is already low when computing the initial condition, and so the benefit of preserving exactly the divergence by using the space $dB_k^{\text{curl}}$ with Godunov flux is low.

The same plots are performed on the series of triangular meshes, shown in Figure 6, and the errors obtained for the variable $e_x$ are shown in Table 4 and Table 5. The same observation as for Cartesian meshes hold.
<table><caption>Table 3: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 5.1.2 with initial condition (31), on the series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.39e+00</td><td></td><td>2.70e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.58e+00</td><td>0.59</td><td>1.93e+00</td><td>0.47</td></tr><tr><td>0.02226197</td><td>9.25e-01</td><td>0.85</td><td>1.19e+00</td><td>0.77</td></tr><tr><td>0.01091023</td><td>4.94e-01</td><td>0.88</td><td>6.58e-01</td><td>0.83</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.43e-01</td><td></td><td>1.54e-01</td><td></td></tr><tr><td>0.04187874</td><td>2.92e-02</td><td>2.27</td><td>3.15e-02</td><td>2.27</td></tr><tr><td>0.02226197</td><td>6.47e-03</td><td>2.38</td><td>7.43e-03</td><td>2.29</td></tr><tr><td>0.01091023</td><td>1.53e-03</td><td>2.02</td><td>1.98e-03</td><td>1.85</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>7.53e-03</td><td></td><td>8.83e-03</td><td></td></tr><tr><td>0.04187874</td><td>9.22e-04</td><td>3.00</td><td>1.28e-03</td><td>2.75</td></tr><tr><td>0.02226197</td><td>1.10e-04</td><td>3.37</td><td>1.76e-04</td><td>3.14</td></tr><tr><td>0.01091023</td><td>1.30e-05</td><td>2.99</td><td>2.49e-05</td><td>2.75</td></tr></tbody></table>

<table><caption>Table 4: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 5.1.2 with initial condition (31), on the series of triangular meshes. Results show a higher benefit than in Table 3 in using the Godunov flux for preserving exactly the divergence.</caption><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.41e+00</td><td></td><td>2.83e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.59e+00</td><td>0.60</td><td>2.08e+00</td><td>0.44</td></tr><tr><td>0.02226197</td><td>9.28e-01</td><td>0.85</td><td>1.34e+00</td><td>0.70</td></tr><tr><td>0.01091023</td><td>4.95e-01</td><td>0.88</td><td>7.88e-01</td><td>0.74</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.53e-01</td><td></td><td>1.94e-01</td><td></td></tr><tr><td>0.04187874</td><td>3.36e-02</td><td>2.17</td><td>4.74e-02</td><td>2.01</td></tr><tr><td>0.02226197</td><td>7.67e-03</td><td>2.34</td><td>1.33e-02</td><td>2.01</td></tr><tr><td>0.01091023</td><td>1.84e-03</td><td>2.00</td><td>4.04e-03</td><td>1.67</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.39e-02</td><td></td><td>2.10e-02</td><td></td></tr><tr><td>0.04187874</td><td>2.15e-03</td><td>2.66</td><td>4.94e-03</td><td>2.07</td></tr><tr><td>0.02226197</td><td>2.57e-04</td><td>3.37</td><td>7.07e-04</td><td>3.08</td></tr><tr><td>0.01091023</td><td>3.03e-05</td><td>3.00</td><td>9.99e-05</td><td>2.74</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure 5: Error obtained on the test case described in subsection 6.1.2 with initial condition (28) on the left, and with the initial condition (29) added to (28) on the right, on a series of Cartesian meshes. For each of the test cases, the error obtained for the variables $b$ (top row), $e_x$ (middle row) and $e_y$ (bottom row) is shown for different approximation spaces and for the Lax-Friedrich and Godunov flux.</figcaption></figure>
with $\bar{x} = \frac{x-x_c}{r_0}$, $\bar{y} = \frac{y-y_c}{r_0}$, $r^2 = (x-x_c)^2 + (y-y_c)^2$, and $\bar{r} = \frac{r}{r_0}$. This initial condition is such that $\nabla \cdot \mathbf{u} = 0$, and so is clearly a stationary solution of (33). This solution was built by considering the potential $e^{-\bar{r}^2/2}$ and by taking its rotated gradient $\nabla^{\perp}$. The numerical parameters are $r_0 = 0.15$, and $x_c = y_c = 0.5$. The computation is led until $t = 3$ on the different meshes represented in Figure 1, and The $L^2$ difference between $(\nabla^{\perp})^* \mathbf{u}$ and its initial value is computed along the time. Two configurations of finite element spaces for approximating $\mathbf{u}$ are used:

* the finite element spaces $\mathbf{dB}_{k}^{\text{div}}$,

* the classical finite element spaces for discontinuous Galerkin methods obtained by tensorization of the scalar basis (note however that in the triangular case, these two approximation spaces match),

and with two different numerical flux:

* the Lax-Friedrich flux,

* the Lax-Friedrich flux with purely normal diffusion (18), which matches in our case with the Godunov' flux.

Numerical results are represented in Figure 5 for degree $k = 0, 1, 2$, and show that the only combination that preserves correctly the curl is the one with the finite element space $\mathbf{dB}_{k}^{\text{div}}$ and with the Godunov' numerical flux.

### 5.2.2 Convergence test

The numerical test is adapted from the one for the divergence convergence test case of subsubsection 5.1.2. The initial condition is

$$ 
\left\{ 
\begin{aligned}
p(\mathbf{x}) &= \frac{\omega}{c^2} \sin(k_{\parallel}\pi y) \cos(k_{\parallel}\pi x) \\
\mathbf{u}_x(\mathbf{x}) &= k_{\parallel}\pi \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x) \\
\mathbf{u}_y(\mathbf{x}) &= k_{\perp}\pi \sin(k_{\perp}\pi y) \cos(k_{\parallel}\pi x)
\end{aligned} 
\right. 
\qquad (34) 
$$

and the exact solution is

$$ 
\left\{ 
\begin{aligned}
p(\mathbf{x}) &= \frac{\omega}{c^2} \sin(k_{\perp}\pi y - \omega t) \cos(k_{\parallel}\pi x) \\
\mathbf{u}_x(\mathbf{x}) &= k_{\parallel}\pi \cos(k_{\perp}\pi y - \omega t) \sin(k_{\parallel}\pi x) \\
\mathbf{u}_y(\mathbf{x}) &= k_{\perp}\pi \sin(k_{\perp}\pi y - \omega t) \cos(k_{\parallel}\pi x).
\end{aligned} 
\right. 
$$

The longitudinal and transverse wave numbers $k_{\parallel}$ and $k_{\perp}$ are data of the test case. The frequency $\omega$ is such that

$$ k_{\parallel}^2 + k_{\perp}^2 = \frac{\omega^2}{\pi^2 c^2}. $$

The numerical parameters are $c = 1$, and $k_{\parallel} = k_{\perp} = 2$.

A second test case is performed, in which a stationary non curl free solution is added to the initial solution (34) (it is similar to the stationary solution of subsubsection 5.2.1, but regular), equal to

$$ 
\left\{ 
\begin{aligned}
b(\mathbf{x}) &= 0 \\
\mathbf{u}_x(\mathbf{x}) &= -2K_0 \alpha \bar{y} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2} \\
\mathbf{u}_y(\mathbf{x}) &= 2K_0 \alpha \bar{x} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2},
\end{aligned} 
\right. 
\qquad (35) 
$$

if $r < r_0$, and 0 otherwise. The numerical parameters are $K_0 = 100$, $r_0 = 0.35$, $x_c = y_c = 0.5$, $\alpha = 4$. The convergence curve for this test case with the two initial conditions are represented in Figure 6 for Cartesian meshes, and the matching table with the convergence errors and rates are summarized in Table 5 and Table 6. For triangular meshes, the convergence curves are plotted in Figure 7, and the errors and convergence rates are summarized in the Table 7 and Table 8. Same comments as for the convergence test case of the previous section hold.
<table><caption>Table 3: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 6.1.2 with initial condition (28), on a series of Cartesian meshes. Results show a high benefit in using the space $dB_k^{curl}$ with the Godunov flux, which is the only one to always reach the optimal order.</caption><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>curl</sup><sub>0</sub></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>curl</sup><sub>0</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.85e+00</td><td></td><td>2.78e+00</td><td></td><td>3.12e+00</td><td></td><td>3.11e+00</td><td></td></tr><tr><td>0.05</td><td>2.10e+00</td><td>0.44</td><td>2.03e+00</td><td>0.46</td><td>2.58e+00</td><td>0.27</td><td>2.56e+00</td><td>0.28</td></tr><tr><td>0.025</td><td>1.33e+00</td><td>0.66</td><td>1.26e+00</td><td>0.69</td><td>1.81e+00</td><td>0.51</td><td>1.78e+00</td><td>0.52</td></tr><tr><td>0.0125</td><td>7.59e-01</td><td>0.81</td><td>7.09e-01</td><td>0.83</td><td>1.14e+00</td><td>0.67</td><td>1.10e+00</td><td>0.70</td></tr></tbody><tbody><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>curl</sup><sub>1</sub></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>curl</sup><sub>1</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>2.45e-01</td><td></td><td>1.71e-01</td><td></td><td>2.52e-01</td><td></td><td>1.98e-01</td><td></td></tr><tr><td>0.05</td><td>7.29e-02</td><td>1.75</td><td>3.94e-02</td><td>2.12</td><td>5.84e-02</td><td>2.11</td><td>4.11e-02</td><td>2.26</td></tr><tr><td>0.025</td><td>2.28e-02</td><td>1.68</td><td>9.70e-03</td><td>2.02</td><td>1.31e-02</td><td>2.15</td><td>9.83e-03</td><td>2.07</td></tr><tr><td>0.0125</td><td>6.96e-03</td><td>1.71</td><td>2.42e-03</td><td>2.00</td><td>2.87e-03</td><td>2.20</td><td>2.43e-03</td><td>2.02</td></tr></tbody><tbody><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>curl</sup><sub>2</sub></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>curl</sup><sub>2</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>3.15e-02</td><td></td><td>1.12e-02</td><td></td><td>2.94e-02</td><td></td><td>1.38e-02</td><td></td></tr><tr><td>0.05</td><td>6.47e-03</td><td>2.29</td><td>2.07e-03</td><td>2.43</td><td>6.48e-03</td><td>2.18</td><td>2.72e-03</td><td>2.35</td></tr><tr><td>0.025</td><td>1.15e-03</td><td>2.50</td><td>2.53e-04</td><td>3.04</td><td>1.14e-03</td><td>2.50</td><td>3.96e-04</td><td>2.78</td></tr><tr><td>0.0125</td><td>1.93e-04</td><td>2.57</td><td>3.17e-05</td><td>3.00</td><td>1.93e-04</td><td>2.57</td><td>5.41e-05</td><td>2.87</td></tr></tbody></table>

<table><caption>Table 4: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 6.1.2 with initial condition (28), on the series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th colspan="2"></th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.39e+00</td><td></td><td>2.70e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.58e+00</td><td>0.59</td><td>1.93e+00</td><td>0.47</td></tr><tr><td>0.02226197</td><td>9.25e-01</td><td>0.85</td><td>1.19e+00</td><td>0.77</td></tr><tr><td>0.01091023</td><td>4.94e-01</td><td>0.88</td><td>6.58e-01</td><td>0.83</td></tr></tbody><tbody><tr><th colspan="2"></th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.08439823</td><td>1.43e-01</td><td></td><td>1.54e-01</td><td></td></tr><tr><td>0.04187874</td><td>2.92e-02</td><td>2.27</td><td>3.15e-02</td><td>2.27</td></tr><tr><td>0.02226197</td><td>6.47e-03</td><td>2.38</td><td>7.43e-03</td><td>2.29</td></tr><tr><td>0.01091023</td><td>1.53e-03</td><td>2.02</td><td>1.98e-03</td><td>1.85</td></tr></tbody><tbody><tr><th colspan="2"></th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.08439823</td><td>7.53e-03</td><td></td><td>8.83e-03</td><td></td></tr><tr><td>0.04187874</td><td>9.22e-04</td><td>3.00</td><td>1.28e-03</td><td>2.75</td></tr><tr><td>0.02226197</td><td>1.10e-04</td><td>3.37</td><td>1.76e-04</td><td>3.14</td></tr><tr><td>0.01091023</td><td>1.30e-05</td><td>2.99</td><td>2.49e-05</td><td>2.75</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Figure 5: Plot of $\| (\nabla^{\perp})^* \mathbf{u} - (\nabla^{\perp})^* \mathbf{u}_0 \|_2$ with respect to time for different degree, approximation space, type of meshes and numerical flux. In the left column, Cartesian and unstructured quadrangular meshes are considered, with the $d\mathbf{Q}_k$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. In the right column, Cartesian, unstructured quadrangular and triangular meshes are considered, with the $dB_k^{\text{div}}$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. On triangular meshes, the space $dB_k^{\text{div}}$ and $d\mathbf{P}_k$ are the same, and the computations on these spaces are represented only on the right column. Note that the y scaling of the top right figure ($dB_k^{\text{div}}$ with Godunov' scheme) is $10^{-12}$, whereas it is of the order of 1 for the other plots.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 6: Error obtained on the test case described in subsection 6.1.2 with initial condition (28) on the left, and with the initial condition (29) added to (28) on the right, on a series of unstructured triangular meshes. For each of the test cases, the error obtained for the variables $b$ (top row), $e_x$ (middle row) and $e_y$ (bottom row) is shown for different degrees, and for the Lax-Friedrich and Godunov numerical flux.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: Error obtained on the test case described in subsection 5.2.2 with initial condition (34) on the left, and with the initial condition (35) added to (34) on the right, on a series of Cartesian meshes. For each of the test cases, the error obtained for the variables $p$ (top row), $u_x$ (middle row) and $u_y$ (bottom row) is shown for different approximation spaces and for the Lax-Friedrich and Godunov flux.</figcaption></figure>
<table><caption>Table 5: Errors and convergence rates obtained on the variable $e_x$ with the test case described in subsection 6.1.2 with initial condition (28), on the series of triangular meshes. Results show a higher benefit than in Table 4 in using the Godunov flux for preserving exactly the divergence.</caption><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.41e+00</td><td></td><td>2.83e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.59e+00</td><td>0.60</td><td>2.08e+00</td><td>0.44</td></tr><tr><td>0.02226197</td><td>9.28e-01</td><td>0.85</td><td>1.34e+00</td><td>0.70</td></tr><tr><td>0.01091023</td><td>4.95e-01</td><td>0.88</td><td>7.88e-01</td><td>0.74</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.53e-01</td><td></td><td>1.94e-01</td><td></td></tr><tr><td>0.04187874</td><td>3.36e-02</td><td>2.17</td><td>4.74e-02</td><td>2.01</td></tr><tr><td>0.02226197</td><td>7.67e-03</td><td>2.34</td><td>1.33e-02</td><td>2.01</td></tr><tr><td>0.01091023</td><td>1.84e-03</td><td>2.00</td><td>4.04e-03</td><td>1.67</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.39e-02</td><td></td><td>2.10e-02</td><td></td></tr><tr><td>0.04187874</td><td>2.15e-03</td><td>2.66</td><td>4.94e-03</td><td>2.07</td></tr><tr><td>0.02226197</td><td>2.57e-04</td><td>3.37</td><td>7.07e-04</td><td>3.08</td></tr><tr><td>0.01091023</td><td>3.03e-05</td><td>3.00</td><td>9.99e-05</td><td>2.74</td></tr></tbody></table>

## 6.2 Discrete conservation of a curl: Wave system

In this section, the wave system

$$ \left\{ \begin{array}{l} \partial_t p + \nabla \cdot \mathbf{u} = 0 \\ \partial_t \mathbf{u} + c^2 \nabla p = 0, \end{array} \right. \qquad (30) $$

which couples the pressure $p$ and the velocity $\mathbf{u}$ is considered. The wave velocity $c$ is a parameter of the system, and will be always equal to 1 in the numerical applications.

### 6.2.1 Conservation of the curl of a stationary solution

For this first test case, the following initial condition is imposed

$$ \left\{ \begin{array}{l} p(\mathbf{x}) = 0 \\ \mathbf{u}_x(\mathbf{x}) = -\bar{y} e^{-\bar{r}^2/2} \\ \mathbf{u}_y(\mathbf{x}) = \bar{x} e^{-\bar{r}^2/2}, \end{array} \right. $$

with $\bar{x} = \frac{x - x_c}{r_0}$, $\bar{y} = \frac{y - y_c}{r_0}$, $r^2 = (x - x_c)^2 + (y - y_c)^2$, and $\bar{r} = \frac{r}{r_0}$. This initial condition is such that $\nabla \cdot \mathbf{u} = 0$, and so is clearly a stationary solution of (30). This solution was built by considering the potential $e^{-\bar{r}^2/2}$ and by taking its rotated gradient $\nabla^{\perp}$. The numerical parameters are $r_0 = 0.15$, and $x_c = y_c = 0.5$. The computation is led until $t = 3$ on the different meshes represented in Figure 3, and The $L^2$ difference between $(\nabla^{\perp})^* \mathbf{u}$ and its initial value is computed along the time. Two configurations of finite element spaces for approximating $\mathbf{u}$ are used:

* the finite element spaces $\mathbf{d}\mathbf{B}_k^{\text{div}}$,

* the classical finite element spaces for discontinuous Galerkin methods obtained by tensorization of the scalar basis (note however that in the triangular case, these two approximation spaces match),

and with two different numerical flux:
<table><caption>Table 5: Errors and convergence rates obtained on the variable $u_x$ with the test case described in subsection 5.2.2 with initial condition (34), on a series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>div</sup><sub>0</sub></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>div</sup><sub>0</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.76e+00</td><td></td><td>2.76e+00</td><td></td><td>2.99e+00</td><td></td><td>2.99e+00</td><td></td></tr><tr><td>0.05</td><td>2.02e+00</td><td>0.45</td><td>2.02e+00</td><td>0.45</td><td>2.44e+00</td><td>0.30</td><td>2.44e+00</td><td>0.30</td></tr><tr><td>0.025</td><td>1.26e+00</td><td>0.68</td><td>1.26e+00</td><td>0.68</td><td>1.65e+00</td><td>0.56</td><td>1.65e+00</td><td>0.56</td></tr><tr><td>0.0125</td><td>7.08e-01</td><td>0.83</td><td>7.08e-01</td><td>0.83</td><td>9.76e-01</td><td>0.76</td><td>9.76e-01</td><td>0.76</td></tr><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>div</sup><sub>1</sub></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>div</sup><sub>1</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>1.84e-01</td><td></td><td>1.54e-01</td><td></td><td>1.75e-01</td><td></td><td>1.54e-01</td><td></td></tr><tr><td>0.05</td><td>4.23e-02</td><td>2.12</td><td>3.57e-02</td><td>2.11</td><td>3.88e-02</td><td>2.17</td><td>3.57e-02</td><td>2.11</td></tr><tr><td>0.025</td><td>1.04e-02</td><td>2.03</td><td>8.72e-03</td><td>2.03</td><td>9.28e-03</td><td>2.06</td><td>8.72e-03</td><td>2.03</td></tr><tr><td>0.0125</td><td>2.60e-03</td><td>2.00</td><td>2.17e-03</td><td>2.01</td><td>2.29e-03</td><td>2.02</td><td>2.17e-03</td><td>2.01</td></tr><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>div</sup><sub>2</sub></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>div</sup><sub>2</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>6.41e-03</td><td></td><td>4.90e-03</td><td></td><td>7.17e-03</td><td></td><td>5.98e-03</td><td></td></tr><tr><td>0.05</td><td>8.16e-04</td><td>2.97</td><td>6.08e-04</td><td>3.01</td><td>9.32e-04</td><td>2.94</td><td>7.61e-04</td><td>2.97</td></tr><tr><td>0.025</td><td>1.03e-04</td><td>2.98</td><td>7.59e-05</td><td>3.00</td><td>1.19e-04</td><td>2.97</td><td>9.57e-05</td><td>2.99</td></tr><tr><td>0.0125</td><td>1.30e-05</td><td>2.99</td><td>9.48e-06</td><td>3.00</td><td>1.49e-05</td><td>2.99</td><td>1.20e-05</td><td>3.00</td></tr></tbody></table>

<table><caption>Table 6: Errors and convergence rates obtained on the variable $u_x$ with the test case described in subsection 5.2.2 with initial condition (34), on a series of Cartesian meshes. Results show a high benefit in using the space $dB_{k}^{curl}$ with the Godunov flux, which is the only one to always reach the optimal order.</caption><thead><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>div</sup><sub>0</sub></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>div</sup><sub>0</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.85e+00</td><td></td><td>2.78e+00</td><td></td><td>3.12e+00</td><td></td><td>3.11e+00</td><td></td></tr><tr><td>0.05</td><td>2.10e+00</td><td>0.44</td><td>2.03e+00</td><td>0.46</td><td>2.58e+00</td><td>0.27</td><td>2.56e+00</td><td>0.28</td></tr><tr><td>0.025</td><td>1.33e+00</td><td>0.66</td><td>1.26e+00</td><td>0.69</td><td>1.81e+00</td><td>0.51</td><td>1.78e+00</td><td>0.52</td></tr><tr><td>0.0125</td><td>7.59e-01</td><td>0.81</td><td>7.09e-01</td><td>0.83</td><td>1.14e+00</td><td>0.67</td><td>1.10e+00</td><td>0.70</td></tr><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>div</sup><sub>1</sub></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>div</sup><sub>1</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>2.45e-01</td><td></td><td>1.71e-01</td><td></td><td>2.52e-01</td><td></td><td>1.98e-01</td><td></td></tr><tr><td>0.05</td><td>7.29e-02</td><td>1.75</td><td>3.94e-02</td><td>2.12</td><td>5.84e-02</td><td>2.11</td><td>4.11e-02</td><td>2.26</td></tr><tr><td>0.025</td><td>2.28e-02</td><td>1.68</td><td>9.70e-03</td><td>2.02</td><td>1.31e-02</td><td>2.15</td><td>9.83e-03</td><td>2.07</td></tr><tr><td>0.0125</td><td>6.96e-03</td><td>1.71</td><td>2.42e-03</td><td>2.00</td><td>2.87e-03</td><td>2.20</td><td>2.43e-03</td><td>2.02</td></tr><tr><th colspan="4">Godunov</th><th colspan="4">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>div</sup><sub>2</sub></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>div</sup><sub>2</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td>0.1</td><td>3.15e-02</td><td></td><td>1.12e-02</td><td></td><td>2.94e-02</td><td></td><td>1.38e-02</td><td></td></tr><tr><td>0.05</td><td>6.47e-03</td><td>2.29</td><td>2.07e-03</td><td>2.43</td><td>6.48e-03</td><td>2.18</td><td>2.72e-03</td><td>2.35</td></tr><tr><td>0.025</td><td>1.15e-03</td><td>2.50</td><td>2.53e-04</td><td>3.04</td><td>1.14e-03</td><td>2.50</td><td>3.96e-04</td><td>2.78</td></tr><tr><td>0.0125</td><td>1.93e-04</td><td>2.57</td><td>3.17e-05</td><td>3.00</td><td>1.93e-04</td><td>2.57</td><td>5.41e-05</td><td>2.87</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 7: Plot of $|| (\nabla^{\perp})^* \mathbf{u} - (\nabla^{\perp})^* \mathbf{u}_0 ||_2$ with respect to time for different degree, approximation space, type of meshes and numerical flux. In the left column, Cartesian and unstructured quadrangular meshes are considered, with the $d\mathbf{Q}_k$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. In the right column, Cartesian, unstructured quadrangular and triangular meshes are considered, with the $d\mathbf{B}_k^{\text{div}}$ approximation space with the Godunov (top figure) and Lax-Friedrich (bottom figure) numerical flux. On triangular meshes, the space $d\mathbf{B}_k^{\text{div}}$ and $d\mathbf{P}_k$ are the same, and the computations on these spaces are represented only on the right column. Note that the y scaling of the top right figure ($d\mathbf{B}_k^{\text{div}}$ with Godunov' scheme) is $10^{-12}$, whereas it is of the order of 1 for the other plots.</figcaption></figure>

* the Lax-Friedrich flux,

* the Lax-Friedrich flux with purely normal diffusion (14), which matches in our case with the Godunov' flux.

Numerical results are represented in Figure 7 for degree $k = 0, 1, 2$, and show that the only combination that preserves correctly the curl is the one with the finite element space $d\mathbf{B}_k^{\text{div}}$ and with the Godunov' numerical flux.

### 6.2.2 Convergence test

The numerical test is adapted from the one for the divergence convergence test case of sub-section 6.1.2. The initial condition is

$$ \left\{ \begin{array}{l} p(\mathbf{x}) = \frac{\omega}{c^2} \sin(k_{\perp}\pi y) \cos(k_{\parallel}\pi x) \\ \mathbf{u}_x(\mathbf{x}) = k_{\parallel}\pi \cos(k_{\perp}\pi y) \sin(k_{\parallel}\pi x) \\ \mathbf{u}_y(\mathbf{x}) = k_{\perp}\pi \sin(k_{\parallel}\pi y) \cos(k_{\perp}\pi x) \end{array} \right. \qquad (31) $$
<figure><img src="image_10.png" /><figcaption>Figure 7: Error obtained on the test case described in subsection 5.2.2 with initial condition (34) on the left, and with the initial condition (35) added to (34) on the right, on a series of unstructured triangular meshes. For each of the test cases, the error obtained for the variables $p$ (top row), $u_x$ (middle row) and $u_y$ (bottom row) is shown for different degrees, and for the Lax-Friedrich and Godunov numerical flux.</figcaption></figure>
<table><caption>Table 6: Errors and convergence rates obtained on the variable $\mathbf{u}_x$ with the test case described in subsection 6.2.2 with initial condition (31), on a series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sub>0</sub><sup>div</sup></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sub>0</sub><sup>div</sup></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.76e+00</td><td></td><td>2.76e+00</td><td></td><td>2.99e+00</td><td></td><td>2.99e+00</td><td></td></tr><tr><td>0.05</td><td>2.02e+00</td><td>0.45</td><td>2.02e+00</td><td>0.45</td><td>2.44e+00</td><td>0.30</td><td>2.44e+00</td><td>0.30</td></tr><tr><td>0.025</td><td>1.26e+00</td><td>0.68</td><td>1.26e+00</td><td>0.68</td><td>1.65e+00</td><td>0.56</td><td>1.65e+00</td><td>0.56</td></tr><tr><td>0.0125</td><td>7.08e-01</td><td>0.83</td><td>7.08e-01</td><td>0.83</td><td>9.76e-01</td><td>0.76</td><td>9.76e-01</td><td>0.76</td></tr></tbody></table><table><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sub>1</sub><sup>div</sup></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sub>1</sub><sup>div</sup></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>1.84e-01</td><td></td><td>1.54e-01</td><td></td><td>1.75e-01</td><td></td><td>1.54e-01</td><td></td></tr><tr><td>0.05</td><td>4.23e-02</td><td>2.12</td><td>3.57e-02</td><td>2.11</td><td>3.88e-02</td><td>2.17</td><td>3.57e-02</td><td>2.11</td></tr><tr><td>0.025</td><td>1.04e-02</td><td>2.03</td><td>8.72e-03</td><td>2.03</td><td>9.28e-03</td><td>2.06</td><td>8.72e-03</td><td>2.03</td></tr><tr><td>0.0125</td><td>2.60e-03</td><td>2.00</td><td>2.17e-03</td><td>2.01</td><td>2.29e-03</td><td>2.02</td><td>2.17e-03</td><td>2.01</td></tr></tbody></table><table><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sub>2</sub><sup>div</sup></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sub>2</sub><sup>div</sup></th></tr><tr><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>6.41e-03</td><td></td><td>4.90e-03</td><td></td><td>7.17e-03</td><td></td><td>5.98e-03</td><td></td></tr><tr><td>0.05</td><td>8.16e-04</td><td>2.97</td><td>6.08e-04</td><td>3.01</td><td>9.32e-04</td><td>2.94</td><td>7.61e-04</td><td>2.97</td></tr><tr><td>0.025</td><td>1.03e-04</td><td>2.98</td><td>7.59e-05</td><td>3.00</td><td>1.19e-04</td><td>2.97</td><td>9.57e-05</td><td>2.99</td></tr><tr><td>0.0125</td><td>1.30e-05</td><td>2.99</td><td>9.48e-06</td><td>3.00</td><td>1.49e-05</td><td>2.99</td><td>1.20e-05</td><td>3.00</td></tr></tbody></table>

and the exact solution is

$$ egin{cases} p(\mathbf{x}) = \frac{\omega}{c^2} \sin(k_{\parallel}\pi y - \omega t) \cos(k_{\parallel}\pi x) \\ \mathbf{u}_x(\mathbf{x}) = k_{\parallel}\pi \cos(k_{\parallel}\pi y - \omega t) \sin(k_{\parallel}\pi x) \\ \mathbf{u}_y(\mathbf{x}) = k_{\perp}\pi \sin(k_{\perp}\pi y - \omega t) \cos(k_{\perp}\pi x). \end{cases} $$

The longitudinal and transverse wave numbers $k_{\parallel}$ and $k_{\perp}$ are data of the test case. The frequency $\omega$ is such that

$$ k_{\parallel}^2 + k_{\perp}^2 = \frac{\omega^2}{\pi^2 c^2}. $$

The numerical parameters are $c = 1$, and $k_{\parallel} = k_{\perp} = 2$.

A second test case is performed, in which a stationary non curl free solution is added to the initial solution (31) (it is similar to the stationary solution of subsection 6.2.1, but regular), equal to

$$ egin{cases} b(\mathbf{x}) = 0 \\ \mathbf{u}_x(\mathbf{x}) = -2K_0 \alpha \bar{y} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2} \\ \mathbf{u}_y(\mathbf{x}) = 2K_0 \alpha \bar{x} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r}^2)^2}, \end{cases} \quad (32) $$

if $r < r_0$, and 0 otherwise. The numerical parameters are $K_0 = 100$, $r_0 = 0.35$, $x_c = y_c = 0.5$, $\alpha = 4$. The convergence curve for this test case with the two initial conditions are represented in Figure 8 for Cartesian meshes, and the matching table with the convergence errors and rates are summarized in Table 6 and Table 7. For triangular meshes, the convergence curves are plotted in Figure 9, and the errors and convergence rates are summarized in the Table 8 and Table 9. Same comments as for the convergence test case of the previous section hold.
<table><caption>Table 7: Errors and convergence rates obtained on the variable $u_x$ with the test case described in subsection 5.2.2 with initial condition (34), on the series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.39e+00</td><td></td><td>2.70e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.58e+00</td><td>0.59</td><td>1.94e+00</td><td>0.47</td></tr><tr><td>0.02226197</td><td>9.29e-01</td><td>0.84</td><td>1.20e+00</td><td>0.77</td></tr><tr><td>0.01091023</td><td>4.93e-01</td><td>0.89</td><td>6.58e-01</td><td>0.84</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.43e-01</td><td></td><td>1.55e-01</td><td></td></tr><tr><td>0.04187874</td><td>2.94e-02</td><td>2.26</td><td>3.15e-02</td><td>2.27</td></tr><tr><td>0.02226197</td><td>6.48e-03</td><td>2.39</td><td>7.41e-03</td><td>2.29</td></tr><tr><td>0.01091023</td><td>1.53e-03</td><td>2.03</td><td>1.97e-03</td><td>1.86</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>7.62e-03</td><td></td><td>8.87e-03</td><td></td></tr><tr><td>0.04187874</td><td>9.12e-04</td><td>3.03</td><td>1.28e-03</td><td>2.76</td></tr><tr><td>0.02226197</td><td>1.09e-04</td><td>3.35</td><td>1.78e-04</td><td>3.12</td></tr><tr><td>0.01091023</td><td>1.30e-05</td><td>2.98</td><td>2.49e-05</td><td>2.76</td></tr></tbody></table>

<table><caption>Table 8: Errors and convergence rates obtained on the variable $u_x$ with the test case described in subsection 5.2.2 with initial condition (34), on the series of triangular meshes. Results show a higher benefit than in Table 7 in using the Godunov flux for preserving exactly the divergence.</caption><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.41e+00</td><td></td><td>2.83e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.59e+00</td><td>0.59</td><td>2.08e+00</td><td>0.43</td></tr><tr><td>0.02226197</td><td>9.33e-01</td><td>0.85</td><td>1.34e+00</td><td>0.69</td></tr><tr><td>0.01091023</td><td>4.95e-01</td><td>0.89</td><td>7.87e-01</td><td>0.75</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.55e-01</td><td></td><td>1.94e-01</td><td></td></tr><tr><td>0.04187874</td><td>3.37e-02</td><td>2.18</td><td>4.63e-02</td><td>2.05</td></tr><tr><td>0.02226197</td><td>7.66e-03</td><td>2.34</td><td>1.34e-02</td><td>1.96</td></tr><tr><td>0.01091023</td><td>1.84e-03</td><td>2.00</td><td>4.13e-03</td><td>1.65</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.40e-02</td><td></td><td>1.99e-02</td><td></td></tr><tr><td>0.04187874</td><td>2.14e-03</td><td>2.68</td><td>4.94e-03</td><td>1.99</td></tr><tr><td>0.02226197</td><td>2.54e-04</td><td>3.37</td><td>7.14e-04</td><td>3.06</td></tr><tr><td>0.01091023</td><td>3.06e-05</td><td>2.97</td><td>9.98e-05</td><td>2.76</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 8: Error obtained on the test case described in subsection 6.2.2 with initial condition (31) on the left, and with the initial condition (32) added to (31) on the right, on a series of Cartesian meshes. For each of the test cases, the error obtained for the variables $p$ (top row), $u_x$ (middle row) and $u_y$ (bottom row) is shown for different approximation spaces and for the Lax-Friedrich and Godunov flux.</figcaption></figure>
<table><caption>Table 9: $L^2$ norm of the initial divergence when the initial condition is computed as in Remark 6 for the different meshes of Figure 1 and for degree 0, 1 and 2.</caption><thead><tr><th>degree</th><th>0</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>Triangular</td><td>0</td><td>2.9602e-14</td><td>9.43413e-14</td></tr><tr><td>Cartesian</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Unstructured quad</td><td>0</td><td>4.26275e-14</td><td>1.69086e-13</td></tr></tbody></table>

## 5.3 Induction equation

In this section, we are interested in the system (23). The test case is taken from [82, Section 4.3] Rotating discontinuous magnetic field loop, still, it was modified in order to ensure that the magnetic loop is regular because we wish to do a convergence study. The computational domain is still $[0, 1]^2$, and the vector field $v$ is an orthonormal velocity with respect to the center of the computational domain (0.5, 0.5): $v = -e_\theta$.

The initial condition is defined by

$$ \mathbf{u}^0(\mathbf{x}) = \begin{cases} \text{if } \bar{r} < r_0 \\ \mathbf{u}_x = -2K_0\alpha\bar{y} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r})^2} \\ \mathbf{u}_y = 2K_0\alpha\bar{x} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r})^2} \end{cases} \quad (36) $$

where $r^2 = (x - x_c)^2 + (y - y_c)^2$, $\bar{r} = r/r_0$, $\bar{x} = (x - x_c)/r_0$, $\bar{y} = (y - y_c)/r_0$, and $K_0, \alpha, x_c, y_c$ and $r_0$ are numerical parameters of the test case. The numerical parameters are $K_0 = 2, \alpha = 4, x_c = 0.5, y_c = 0.75$, and $r_0 = 0.125$. This solution is such that $\mathbf{u}^0 = \nabla^\perp f^0$, where $f^0$ is the regular function

$$ f^0(\mathbf{x}) := \begin{cases} -K_0 \bar{r}_0 e^{-\alpha/(1-\bar{r}^2)} & \text{if } \bar{r} < \bar{r}_0 \\ 0 & \text{otherwise.} \end{cases} $$

The initial solution is a vortex rotating around the point $(x_c, y_c)$. The field $v$ induces a rotation of the solution around (0.5, 0.5), which means that the exact solution is a vortex rotating around its center, and the center of the vortex rotates around (0.5, 0.5) with angular velocity equal to 1, namely

$$ \mathbf{u}(\mathbf{x}, t) = R(-t)\mathbf{u}^0(R(t)\mathbf{x}), $$

where $R(t)$ is the matrix of rotation around (0.5, 0.5) of angle $t$:

$$ R(t) := \begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix}. $$

The computations are led with the approximation space $\mathbf{d}\mathcal{B}_k^{\text{curl}}$, and with the Godunov' flux. The initial condition is computed with the method explained in Remark 6.

### 5.3.1 Conservation of the divergence free field

We first check that the divergence of $\mathbf{u}$ is preserved, equal to 0. For this, the test case is run until time $t = \pi$ on the different meshes shown in Figure 1. The initial $L^2$ norm of $\nabla^* \mathbf{u}^0$ is summarized in Table 9. Results show that the method described in Remark 6 ensures the zero divergence for the initial condition.

The difference with respect to the initial divergence is plotted in Figure 8, and show preservation up to round-off errors of the divergence free field.
<table><caption>Table 7: Errors and convergence rates obtained on the variable $u_x$ with the test case described in subsection 6.2.2 with initial condition (31), on a series of Cartesian meshes. Results show a high benefit in using the space $dB_k^{\text{curl}}$ with the Godunov flux, which is the only one to always reach the optimal order.</caption><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>div</sup><sub>0</sub></th><th colspan="2">dQ<sub>0</sub></th><th colspan="2">dB<sup>div</sup><sub>0</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.85e+00</td><td></td><td>2.78e+00</td><td></td><td>3.12e+00</td><td></td><td>3.11e+00</td><td></td></tr><tr><td>0.05</td><td>2.10e+00</td><td>0.44</td><td>2.03e+00</td><td>0.46</td><td>2.58e+00</td><td>0.27</td><td>2.56e+00</td><td>0.28</td></tr><tr><td>0.025</td><td>1.33e+00</td><td>0.66</td><td>1.26e+00</td><td>0.69</td><td>1.81e+00</td><td>0.51</td><td>1.78e+00</td><td>0.52</td></tr><tr><td>0.0125</td><td>7.59e-01</td><td>0.81</td><td>7.09e-01</td><td>0.83</td><td>1.14e+00</td><td>0.67</td><td>1.10e+00</td><td>0.70</td></tr></tbody><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>div</sup><sub>1</sub></th><th colspan="2">dQ<sub>1</sub></th><th colspan="2">dB<sup>div</sup><sub>1</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>2.45e-01</td><td></td><td>1.71e-01</td><td></td><td>2.52e-01</td><td></td><td>1.98e-01</td><td></td></tr><tr><td>0.05</td><td>7.29e-02</td><td>1.75</td><td>3.94e-02</td><td>2.12</td><td>5.84e-02</td><td>2.11</td><td>4.11e-02</td><td>2.26</td></tr><tr><td>0.025</td><td>2.28e-02</td><td>1.68</td><td>9.70e-03</td><td>2.02</td><td>1.31e-02</td><td>2.15</td><td>9.83e-03</td><td>2.07</td></tr><tr><td>0.0125</td><td>6.96e-03</td><td>1.71</td><td>2.42e-03</td><td>2.00</td><td>2.87e-03</td><td>2.20</td><td>2.43e-03</td><td>2.02</td></tr></tbody><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>div</sup><sub>2</sub></th><th colspan="2">dQ<sub>2</sub></th><th colspan="2">dB<sup>div</sup><sub>2</sub></th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.1</td><td>3.15e-02</td><td></td><td>1.12e-02</td><td></td><td>2.94e-02</td><td></td><td>1.38e-02</td><td></td></tr><tr><td>0.05</td><td>6.47e-03</td><td>2.29</td><td>2.07e-03</td><td>2.43</td><td>6.48e-03</td><td>2.18</td><td>2.72e-03</td><td>2.35</td></tr><tr><td>0.025</td><td>1.15e-03</td><td>2.50</td><td>2.53e-04</td><td>3.04</td><td>1.14e-03</td><td>2.50</td><td>3.96e-04</td><td>2.78</td></tr><tr><td>0.0125</td><td>1.93e-04</td><td>2.57</td><td>3.17e-05</td><td>3.00</td><td>1.93e-04</td><td>2.57</td><td>5.41e-05</td><td>2.87</td></tr></tbody></table>

<table><caption>Table 8: Errors and convergence rates obtained on the variable $u_x$ with the test case described in subsection 6.2.2 with initial condition (31), on the series of triangular meshes. Results show a low benefit in using the Godunov flux, namely in exactly preserving the divergence.</caption><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.39e+00</td><td></td><td>2.70e+00</td><td></td><td>2.70e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.58e+00</td><td>0.59</td><td>1.94e+00</td><td>0.47</td><td>1.94e+00</td><td></td></tr><tr><td>0.02226197</td><td>9.29e-01</td><td>0.84</td><td>1.20e+00</td><td>0.77</td><td>1.20e+00</td><td></td></tr><tr><td>0.01091023</td><td>4.93e-01</td><td>0.89</td><td>6.58e-01</td><td>0.84</td><td>6.58e-01</td><td></td></tr></tbody><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.43e-01</td><td></td><td>1.55e-01</td><td></td><td>1.55e-01</td><td></td></tr><tr><td>0.04187874</td><td>2.94e-02</td><td>2.26</td><td>3.15e-02</td><td>2.27</td><td>3.15e-02</td><td></td></tr><tr><td>0.02226197</td><td>6.48e-03</td><td>2.39</td><td>7.41e-03</td><td>2.29</td><td>7.41e-03</td><td></td></tr><tr><td>0.01091023</td><td>1.53e-03</td><td>2.03</td><td>1.97e-03</td><td>1.86</td><td>1.97e-03</td><td></td></tr></tbody><thead><tr><th colspan="3">Godunov</th><th colspan="3">Lax-Friedrich</th></tr><tr><th><i>h</i></th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>7.62e-03</td><td></td><td>8.87e-03</td><td></td><td>8.87e-03</td><td></td></tr><tr><td>0.04187874</td><td>9.12e-04</td><td>3.03</td><td>1.28e-03</td><td>2.76</td><td>1.28e-03</td><td></td></tr><tr><td>0.02226197</td><td>1.09e-04</td><td>3.35</td><td>1.78e-04</td><td>3.12</td><td>1.78e-04</td><td></td></tr><tr><td>0.01091023</td><td>1.30e-05</td><td>2.98</td><td>2.49e-05</td><td>2.76</td><td>2.49e-05</td><td></td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 8: Time evolution of $||∇^*u - ∇^*u^0||_2$ with respect to the time for a coarse Cartesian, unstructured quadrangular and a triangular mesh (all shown in Figure 1). Note that the scale is $10^{-11}$, so that the figure shows exact (up to round-off errors) preservation of the divergence.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 9: Error obtained on the test case described in subsection 6.2.2 with initial condition (31) on the left, and with the initial condition (32) added to (31) on the right, on a series of unstructured triangular meshes. For each of the test cases, the error obtained for the variables $p$ (top row), $u_x$ (middle row) and $u_y$ (bottom row) is shown for different degrees, and for the Lax-Friedrich and Godunov numerical flux.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9: Error curves for the test case described in subsection 5.3.2. The top figures are the error on $\mathbf{u}_x$, whereas the bottom figures are the error on $\mathbf{u}_y$. The left column matches with triangular meshes whereas the right column is for Cartesian meshes.</figcaption></figure>

### 5.3.2 Convergence test

We wish now to perform a convergence test on the series of Cartesian and unstructured triangular meshes described in the beginning of the section. In this test, the computation is led until $t = 0.5$. The convergence curves for the two variables $\mathbf{u}_x$ and $\mathbf{u}_y$ are plotted in Figure 9, whereas the convergence rate computations are shown in Table 10. Both the table and the figure represent a convergence close of the optimal order of convergence. Note however that in agreement with Remark 7, the order of accuracy for the variable $\mathbf{u}_x$, degree 1 on Cartesian mesh is 3/2.

### 5.3.3 Stability test

In Proposition 4, the divergence free preservation is obtained for the numerical flux (22), which has purely tangential diffusion. Numerical flux (22) is less diffusive than the Lax-Friedrich numerical flux, which has diffusion in all the directions. The aim of this test is to assess that the tangential diffusion is sufficient for ensuring $L^2$ stability.

For this, we do the same test as in [82, Section 4.2], *Discontinuous magnetic field loop*. The initial condition is

$$ \mathbf{u}^0(\mathbf{x}) = \begin{cases} \text{if } \bar{r} < r_0 & \begin{cases} \mathbf{u}_x = -K_0 \bar{y} \\ \mathbf{u}_y = K_0 \bar{x} \end{cases} \\ \text{if } \bar{r} \ge r_0 & 0, \end{cases} $$

with the same notations as in (36). The numerical parameters are $x_c = y_c = 0.5$, $K_0 = 0.01$. Taking a discontinuous initial condition is well suited with the stability study, because a discontinuous solution contains a lot of Fourier modes, and these Fourier modes are decreasing slowly. Following [82], the velocity vector field is $\mathbf{v} = (1, 1)$ for allowing a long time integration until $t = 2$. The boundary conditions are periodic. The CFL number is equal to 0.5 for $k = 0$, 0.33 for $k = 1$ and 0.2 for $k = 2$, as for the classical discontinuous Galerkin method. The stability is assessed
<table><caption>Table 9: Errors and convergence rates obtained on the variable $\mathbf{u}_x$ with the test case described in subsection 6.2.2 with initial condition (31), on the series of triangular meshes. Results show a higher benefit than in Table 8 in using the Godunov flux for preserving exactly the divergence.</caption><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>2.41e+00</td><td></td><td>2.83e+00</td><td></td></tr><tr><td>0.04187874</td><td>1.59e+00</td><td>0.59</td><td>2.08e+00</td><td>0.43</td></tr><tr><td>0.02226197</td><td>9.33e-01</td><td>0.85</td><td>1.34e+00</td><td>0.69</td></tr><tr><td>0.01091023</td><td>4.95e-01</td><td>0.89</td><td>7.87e-01</td><td>0.75</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.55e-01</td><td></td><td>1.94e-01</td><td></td></tr><tr><td>0.04187874</td><td>3.37e-02</td><td>2.18</td><td>4.63e-02</td><td>2.05</td></tr><tr><td>0.02226197</td><td>7.66e-03</td><td>2.34</td><td>1.34e-02</td><td>1.96</td></tr><tr><td>0.01091023</td><td>1.84e-03</td><td>2.00</td><td>4.13e-03</td><td>1.65</td></tr></tbody><thead><tr><th rowspan="2">h</th><th colspan="2">Godunov</th><th colspan="2">Lax-Friedrich</th></tr><tr><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td>0.08439823</td><td>1.40e-02</td><td></td><td>1.99e-02</td><td></td></tr><tr><td>0.04187874</td><td>2.14e-03</td><td>2.68</td><td>4.94e-03</td><td>1.99</td></tr><tr><td>0.02226197</td><td>2.54e-04</td><td>3.37</td><td>7.14e-04</td><td>3.06</td></tr><tr><td>0.01091023</td><td>3.06e-05</td><td>2.97</td><td>9.98e-05</td><td>2.76</td></tr></tbody></table>

## 6.3 Induction equation

In this section, we are interested in the system (20). The test case is taken from [66, Section 4.3] Rotating discontinuous magnetic field loop, still, it was modified in order to ensure that the magnetic loop is regular because we wish to do a convergence study. The computational domain is still $[0, 1]^2$, and the vector field $\mathbf{v}$ is an orthonormal velocity with respect to the center of the computational domain $(0.5, 0.5)$: $\mathbf{v} = -\mathbf{e}_\theta$.

The initial condition is defined by

$$ \mathbf{u}^0(\mathbf{x}) = \begin{cases} \text{if } \bar{r} < r_0 \\ \mathbf{u}_x = -2K_0\alpha\bar{y} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r})^2} \\ \mathbf{u}_y = 2K_0\alpha\bar{x} \frac{e^{-\alpha/(1-\bar{r}^2)}}{(1-\bar{r})^2} \\ \text{if } \bar{r} \ge r_0 \end{cases}, \quad 0, $$

where $r^2 = (x - x_c)^2 + (y - y_c)^2$, $\bar{r} = r/r_0$, $\bar{x} = x/r_0$, $\bar{y} = y/r_0$, and $K_0$, $\alpha$, $x_c$, $y_c$ and $r_0$ are numerical parameters of the test case. This solution is such that $\mathbf{u}^0 = \nabla^\perp f^0$, where $f^0$ is the regular function

$$ f^0(\mathbf{x}) := \begin{cases} -K_0 \bar{r}_0 e^{-\alpha/(1-\bar{r}^2)} & \text{if } \bar{r} < \bar{r}_0 \\ 0 & \text{otherwise.} \end{cases} $$

The initial solution is a vortex rotating around the point $(x_c, y_c)$. The field $\mathbf{v}$ induces a rotation of the solution around $(0.5, 0.5)$, which means that the exact solution is a vortex rotating around its center, and the center of the vortex rotates around $(0.5, 0.5)$ with angular velocity equal to 1, namely

$$ \mathbf{u}(\mathbf{x}, t) = R(-t)\mathbf{u}^0(R(t)\mathbf{x}), $$

where $R(t)$ is the matrix of rotation around $(0.5, 0.5)$ of angle $t$:

$$ R(t) := \begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix}. $$
<table><caption>Table 10: Errors obtained for the test case of the Rotating regular magnetic loop on triangular and quadrangular meshes. The error obtained is close of the optimal rate of convergence.</caption><thead><tr><th colspan="2">degree 0</th><th colspan="4">Cartesian</th><th colspan="4">Triangle</th></tr><tr><th colspan="2"></th><th colspan="2">u<sub>x</sub></th><th colspan="2">u<sub>y</sub></th><th colspan="2">u<sub>x</sub></th><th colspan="2">u<sub>y</sub></th></tr><tr><th></th><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr></thead><tbody><tr><td></td><td>1.00e-01</td><td>6.15e-03</td><td></td><td>5.84e-03</td><td></td><td>8.44e-02</td><td>6.71e-03</td><td></td><td>5.43e-03</td><td></td></tr><tr><td></td><td>5.00e-02</td><td>5.07e-03</td><td>0.28</td><td>4.57e-03</td><td>0.35</td><td>4.19e-02</td><td>4.72e-03</td><td>0.50</td><td>4.67e-03</td><td>0.22</td></tr><tr><td></td><td>2.50e-02</td><td>3.15e-03</td><td>0.69</td><td>2.83e-03</td><td>0.69</td><td>2.23e-02</td><td>2.84e-03</td><td>0.81</td><td>2.72e-03</td><td>0.85</td></tr><tr><td></td><td>1.25e-02</td><td>1.69e-03</td><td>0.90</td><td>1.50e-03</td><td>0.92</td><td>1.09e-02</td><td>1.61e-03</td><td>0.80</td><td>1.49e-03</td><td>0.84</td></tr><tr><td></td><td>6.25e-03</td><td>8.72e-04</td><td>0.95</td><td>7.79e-04</td><td>0.94</td><td>5.60e-03</td><td>9.44e-04</td><td>0.80</td><td>7.50e-04</td><td>1.03</td></tr><tr><th colspan="2">degree 1</th><th colspan="4">Cartesian</th><th colspan="4">Triangle</th></tr><tr><th colspan="2"></th><th colspan="2">u<sub>x</sub></th><th colspan="2">u<sub>y</sub></th><th colspan="2">u<sub>x</sub></th><th colspan="2">u<sub>y</sub></th></tr><tr><th></th><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td></td><td>1.00e-01</td><td>2.42e-03</td><td></td><td>2.47e-03</td><td></td><td>8.44e-02</td><td>3.80e-03</td><td></td><td>3.85e-03</td><td></td></tr><tr><td></td><td>5.00e-02</td><td>1.29e-03</td><td>0.90</td><td>8.24e-04</td><td>1.58</td><td>4.19e-02</td><td>1.93e-03</td><td>0.97</td><td>1.48e-03</td><td>1.36</td></tr><tr><td></td><td>2.50e-02</td><td>6.01e-04</td><td>1.11</td><td>1.93e-04</td><td>2.10</td><td>2.23e-02</td><td>1.02e-03</td><td>1.01</td><td>4.95e-04</td><td>1.73</td></tr><tr><td></td><td>1.25e-02</td><td>2.26e-04</td><td>1.41</td><td>4.17e-05</td><td>2.21</td><td>1.09e-02</td><td>3.11e-04</td><td>1.66</td><td>1.71e-04</td><td>1.49</td></tr><tr><td></td><td>6.25e-03</td><td>7.89e-05</td><td>1.52</td><td>9.64e-06</td><td>2.11</td><td>5.60e-03</td><td>1.17e-04</td><td>1.46</td><td>6.24e-05</td><td>1.51</td></tr><tr><th colspan="2">degree 2</th><th colspan="4">Cartesian</th><th colspan="4">Triangle</th></tr><tr><th colspan="2"></th><th colspan="2">u<sub>x</sub></th><th colspan="2">u<sub>y</sub></th><th colspan="2">u<sub>x</sub></th><th colspan="2">u<sub>y</sub></th></tr><tr><th></th><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th><th>h</th><th>Error</th><th>rate</th><th>Error</th><th>rate</th></tr><tr><td></td><td>1.00e-01</td><td>1.57e-03</td><td></td><td>1.20e-03</td><td></td><td>8.44e-02</td><td>3.15e-03</td><td></td><td>2.28e-03</td><td></td></tr><tr><td></td><td>5.00e-02</td><td>3.39e-04</td><td>2.21</td><td>2.09e-04</td><td>2.52</td><td>4.19e-02</td><td>8.42e-04</td><td>1.88</td><td>5.59e-04</td><td>2.00</td></tr><tr><td></td><td>2.50e-02</td><td>5.53e-05</td><td>2.62</td><td>2.59e-05</td><td>3.01</td><td>2.23e-02</td><td>1.32e-04</td><td>2.93</td><td>8.35e-05</td><td>3.01</td></tr><tr><td></td><td>1.25e-02</td><td>5.29e-06</td><td>3.39</td><td>2.41e-06</td><td>3.42</td><td>1.09e-02</td><td>2.42e-05</td><td>2.38</td><td>1.55e-05</td><td>2.36</td></tr><tr><td></td><td>6.25e-03</td><td>4.91e-07</td><td>3.43</td><td>2.69e-07</td><td>3.17</td><td>5.60e-03</td><td>3.49e-06</td><td>2.90</td><td>2.06e-06</td><td>3.02</td></tr></tbody></table>

by computing the norm of the magnetic energy $u \cdot u/2$ at each time step, and by checking that it is decreasing. This magnetic energy is adimensioned by its value at time 0. Results obtained on the triangular mesh of Figure 1 are shown in Figure 10 for different numerical flux and different initialization:

* If the divergence free initialization is done as in Remark 6, and the flux (22) is used, the scheme is stable. This can be explained by the fact that the divergence free component of $u$ is preserved, while the flux (22) is sufficient for stabilizing the evolution of the curl free component of $u$. At the fourth revision step of this article, one reviewer claimed that the numerical scheme was unstable in long time, and so, on the bottom right of this Figure, results for a final time of 20 are also depicted, and show stability too.

* If the initial condition is computed by projection, without the process described in Remark 6, and the flux (22) is used, then the scheme is not stable in general, in agreement with Proposition 4: both the divergence free and curl free component of $u$ evolve, and the flux (22) is not sufficient for stabilizing the evolution of the two components.

* If the divergence free initialization is done as in Remark 6, and the Lax-Friedrich flux is used, the scheme is stable, still, the divergence of $u$ is not kept to 0.

In all our numerical test, no instability was observed, still, the thorough analysis of the numerical scheme remains to be theoretically performed. We foresee also an instability that may come from the accumulation of round-off errors:

**Remark 8 (Impact of accumulation of round-off errors).** As seen in Figure 8, the norm of the divergence is low, but seems to be increasing. We infer that it may come from the numerical scheme (26) ensured by the divergence. The numerical scheme (26) preserves an initially zero