The notion of conservation for residual distribution schemes (or fluctuation splitting schemes), with some applications

TL;DR

This paper explores the concept of discrete conservation in hyperbolic conservation laws through fluctuation splitting schemes, demonstrating how existing schemes can be reformulated and extended to satisfy additional conservation laws.

Contribution

It introduces fluctuation splitting schemes for residual distribution and shows how to reformulate existing schemes in flux form to achieve conservation.

Findings

Most known schemes can be reformulated in flux form

New schemes can satisfy additional conservation laws

Provides a unified framework for residual distribution schemes

Abstract

In this paper, we discuss the notion of discrete conservation for hyperbolic conservation laws. We introduce what we call a fluctuation splitting schemes (or residual distribution, also RDS) and show on several examples how this cal lead to new development. In particular, we show that most, if not all known schemes can be rephrased in flux form, and also show how to satisfy additional conservation laws. This review paper is built on \cite{AbgrallConservation,Abgrall2017,paola,svetlana,ABGRALL2018640}.

Figures

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\text{ div }\mathbf{f}(\mathbf{u})=0\text{ for }\mathbf{x}\in {\mathbb{R}}^d$$

$$u(\mathbf{x},0)=u_0(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d.$$

$$\int_{\mathbb{R}^{d}\times\mathbb{R}^{+}}\bigg{(}\dfrac{\partial\varphi}{\partial t}(\mathbf{x},t)\mathbf{u}(\mathbf{x},t)+\nabla\varphi(\mathbf{x},t)\cdot{\mathbf{f}}(\mathbf{u}(\mathbf{x},t))\;d\mathbf{x}\;dt\bigg{)}+\int_{\mathbb{R}^{d}}\varphi(\mathbf{x},0)\mathbf{u}_{0}(\mathbf{x})\;d\mathbf{x}=0$$

$$u(x,t)=\left\{\begin{array}[]{ll}u^{+}(\mathbf{x},t)&\text{ if }(\mathbf{x},t)\in\Omega^{+}\\ u^{-}(\mathbf{x},t)&\text{ if }(\mathbf{x},t)\in\Omega^{-}\end{array}\right.,$$

$$\mathbf{f}(u^{+}(\mathbf{x},t))\cdot \mathbf{n}-\mathbf{f}(u^{-}(\mathbf{x},t))\cdot \mathbf{n}=\sigma (\mathbf{x},t)(u^{+}(\mathbf{x},t)-u^{-}(\mathbf{x},t))$$

$${\nabla }_{\mathbf{u}}{E}^T{\nabla }_{\mathbf{u}}{f}_j={\nabla }_{\mathbf{u}}{g}_j.$$

$$\frac{\partial E}{\partial t}+\text{ div }\mathbf{g}(\mathbf{u})=0.$$

$$\int_{\mathbb{R}^{d}\times\mathbb{R}^{+}}\bigg{(}\dfrac{\partial\varphi}{\partial t}(\mathbf{x},t)E(\mathbf{u})+\nabla\varphi(\mathbf{x},t)\cdot\mathbf{g}(\mathbf{u}(\mathbf{x},t))\;d\mathbf{x}\;dt\bigg{)}+\int_{\mathbb{R}^{d}}\varphi(\mathbf{x},0)E(\mathbf{u}_{0}(\mathbf{x}))\;d\mathbf{x}\leq 0.$$

$$\frac{1}{\Delta x}{\int }_{C_j}{u}^0(x)dx.$$

$$u_{j}^{n+1}=u_{j}^{n}-\dfrac{\Delta t}{\Delta x}\big{(}\hat{f}_{j+1/2}-\hat{f}_{j-1/2}\big{)},$$

$$u_{\Delta }(x,t)=u_j^n\text{ if }(x,t)\in C_{j+1/2}\times [t_n,t_{n+1}[.$$

$$E(u_{j}^{n+1})-E(u_{j}^{n})-\dfrac{\Delta t}{\Delta x}\big{(}\hat{g}_{j+1/2}-\hat{g}_{j-1/2}\big{)}\leq 0$$

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0$$

$$\frac{\partial u}{\partial t}+\frac{1}{2}\frac{\partial u^2}{\partial x}=0,f(u)=\frac{u^2}{2}.$$

$$u_{i}^{n+1}=u_{i}^{n}-\dfrac{\Delta t}{\Delta x}u_{i}^{n}\big{(}u_{i}^{n}-u_{i-1}^{n}\big{)},$$

$$u_{i}^{n+1}=u_{i}^{n}-\dfrac{\Delta t}{\Delta x}\big{(}f(u_{i}^{n})-f(u_{i-1}^{n})\big{)}.$$

$$u_{i}^{n+1}=u_{i}^{n}-\dfrac{\Delta t}{\Delta x}a_{i-1/2}\big{(}u_{i}^{n}-u_{i-1}^{n}\big{)}.$$

$$u_i^{n+1}=(1-\lambda a_{i-1/2})u_i^n+\lambda a_{i-1/2}{u}_{i-1}^n,$$

$$v_{i}^{n}\big{(}u_{i}^{n+1}-u_{i}^{n}\big{)}+\lambda\big{(}\hat{\mathbf{f}}_{i+1/2}-\hat{\mathbf{f}}_{i-1/2}\big{)}=0,\qquad\hat{\mathbf{f}}_{j+1/2}=\frac{(u_{i}^{n})^{2}}{2}.$$

$$g={\int }_{u}v{f}_{v}dv=vf-{\int }_{v}fdv=vf-\theta (v),$$

$${\hat{g}}_{i+1/2}={\bar{v}}_{i+1/2}{\hat{\mathbf{f}}}_{i+1/2}-{\bar{\theta }}_{i+1/2}$$

$$\begin{split}v_{i}^{n}\big{(}u_{i}^{n+1}-u_{i}^{n}\big{)}+\lambda\big{(}\hat{g}_{i+1/2}-\hat{g}_{i-1/2}\big{)}=\lambda\int_{v_{i}}^{v_{i+1}}\big{(}f(v)-\hat{f}_{i+1/2}\big{)}dv+\lambda\int_{v_{i-1}}^{v_{i}}\big{(}f(v)-\hat{f}_{i-11/2}\big{)}dv.\end{split}$$

$$\begin{split}E_{i}^{n+1}&=E_{i}^{n}+v_{i}^{n}(u_{i}^{n+1}-u_{i}^{n})+\frac{1}{2}\int_{0}^{1}E_{uu}((1-s)u_{i}^{n}+su_{i}^{n+1})(u_{i}^{n+1}-u_{i}^{n})^{2}ds\\ &=E_{i}^{n}+v_{i}^{n}(u_{i}^{n+1}-u_{i}^{n})+\frac{\lambda^{2}}{8}\int_{0}^{1}E_{uu}((1-s)u_{i}^{n}+su_{i}^{n+1})\bigg{(}\big{(}u_{i}^{n}\big{)}^{2}-\big{(}u_{i-1}^{n}\big{)}^{2}\bigg{)}^{2}ds\\ &=E_{i}^{n}+v_{i}^{n}(u_{i}^{n+1}-u_{i}^{n})+\frac{\lambda^{2}}{2}a_{i-1/2}^{2}(u_{i}^{n}-u_{i-1}^{n})^{2}\;\int_{0}^{1}E_{uu}((1-s)u_{i}^{n}+su_{i}^{n+1})ds\end{split}$$

$$\begin{split}E_{i}^{n+1}-E_{i}^{n}&+\lambda\big{(}\hat{g}_{i+1/2}-\hat{g}_{i-1/2}\big{)}=\lambda\int_{v_{i}}^{v_{i+1}}\big{(}f(v)-\hat{f}_{i+1/2}\big{)}dv+\lambda\int_{v_{i-1}}^{v_{i}}\big{(}f(v)-\hat{f}_{i-1/2}\big{)}dv\\ &\qquad+\frac{\lambda^{2}}{2}a_{i-1/2}^{2}(u_{i}^{n}-u_{i-1}^{n})^{2}\;\int_{0}^{1}E_{uu}((1-s)u_{i}^{n}+su_{i}^{n+1})ds\end{split}$$

$$\int_{v_{i}}^{v_{i+1}}\big{(}f(v)-\hat{f}_{i+1/2}\big{)}dv<0,\qquad\int_{v_{i-1}}^{v_{i}}\big{(}f(v)-\hat{f}_{i-1/2}\big{)}dv<0.$$

$$\frac{\partial u}{\partial t}+\frac{\partial f(u)}{\partial x}=\mu$$

$$\frac{\partial u}{\partial t}+\frac{1}{2}\frac{\partial u^2}{\partial x}=0.$$

$$\frac{\partial v}{\partial t}+\frac{1}{4}\frac{\partial v^4}{\partial x}=0.$$

$$\sigma =\frac{u^{+}+u^{-}}{2},$$

$${\sigma }^{\prime }=\frac{1}{4}\frac{(v^{+}{)}^4-(v^{-}{)}^4}{v^{+}-v^{-}}=\frac{(v^{+}{)}^3+3(v^{+}{)}^2{v}^{-}+3v^{+}(v^{-}{)}^2+(v^{-}{)}^3}{4},$$

$$\mathbf{f}(\mathbf{u})=\left(\begin{array}{c}\mathbf{m}\\ \frac{\mathbf{m}\otimes \mathbf{m}}{\rho }+p\text{Id}\\ \mathbf{m}\frac{E+p}{\rho }\end{array}\right).$$

Full text

This paper is also a direct consequence of the work of P.L. Roe, in particular [6, 7] where the notion of conservation I will discussed is first introduced. In [8], P.L. Roe mentions the Hermes project, and the role of Dassault Aviation in it. I was suggested by Bruno Stoufflet, now Vice-President R&D and advanced business in this company, to have a detailed look at [7]. To be honnest, at the time, I did not understood anything, and this was the case for several years. I was lucky to work with Katherine Mer, at the time a postdoc, now research engineer at CEA, and she helped me a lot in starting to understand this notion of conservation. The present contribution can be seen as what I managed to understand after many years playing around the very productive notion of residual distribution schemes (or fluctuation splitting schemes), introduced by P. L. Roe.

1 Introduction

The aim of this paper is to discuss some aspects related to the notion of weak solutions of

[TABLE]

Since at least the work of P. Lax, we know that the correct setting to define the notion of solution to (1) is the following: $\mathbf{u}:\mathbb{R}^{d}\times\mathbb{R}^{+}\rightarrow\mathcal{D}\subset\mathbb{R}^{m}$ is a weak solution of (1) is $\mathbf{u}\in L^{1}(\mathbb{R}^{d}\times\mathbb{R}^{+})^{m}\cap L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}^{+})^{m}$, the initial condition in (1b) belongs to $L^{1}(\mathbb{R}^{d})^{m}\cap L^{\infty}(\mathbb{R}^{d})^{m}$ and for any $\varphi\in C^{1}_{0}(\mathbb{R}^{d}\times\mathbb{R}^{+})$, we have

[TABLE]

In (1)-(2), the flux is ${\mathbf{f}}=(f_{1},\ldots,f_{d})$; the functions $f_{i}$ map the open subset $\mathcal{D}$ of $\mathbb{R}^{m}$ to $\mathbb{R}^{m}$, and they are assumed to be $C^{1}$ for simplicity.

If one is looking for piecewise $C^{1}$ solutions, one sees that the solution must satisfy the Rankine-Hugoniot relations. More precisely, if $u$ is defined by $\mathbb{R}^{d}\times\mathbb{R}^{+}=\Omega^{+}\cup\Omega^{-}$ where the boundary $\Sigma=\Omega^{+}\cap\Omega^{-}$ is a regular hypersurface. To make things simple, we assume that $\Sigma=\{(\mathbf{x},t)\in\mathbb{R}^{d}\times\mathbb{R}^{+},\mathbf{x}=\psi(t)\}$ with $\psi$ $C^{1}$. Then if

[TABLE]

with $u^{\pm}$ smooth in $\Omega^{\pm}$, $u$ is a weak solution of (1) if it is a classical solution in the interior of $\Omega^{+}$ and $\Omega^{-}$, and for any point $(\mathbf{x},t)\in\Sigma$, we have

[TABLE]

where ${\mathbf{n}}$ is a normal of $\Sigma$ at $(\mathbf{x},t)$, $\mathbf{u}^{\pm}(\mathbf{x},t)=\lim\limits_{\mathbf{y}\rightarrow\mathbf{x},\mathbf{y}\in\Omega^{\pm}}u^{\pm}(\mathbf{y},t)$ and $\sigma=\dfrac{d\psi}{dt}$.

An important notion is that of entropy. An entropy $E$ is a convex function defined on $\mathcal{D}$ (hence this set is assumed to be convex) such that there exists $\mathbf{g}=(g_{1},\ldots,g_{d})$ $C^{1}$ with, for any $j=1,\ldots d$,

[TABLE]

Hence, if $\mathbf{u}$ is $C^{1}$, we also have that

[TABLE]

It is well known that the weak solutions of (1) are not smooth nor continuous in general, so that the above equality cannot be met in general for weak solutions. It is said that a weak solution $\mathbf{u}$ is an entropy solution if for any positive $\varphi\in C^{1}_{0}(\mathbb{R}^{d}\times\mathbb{R}^{+})$, we have

[TABLE]

Details can be obtained in classical references such as [9, 10].

The whole purpose is to define a suitable numerical framework for approximating (1) such that, when a sequence of meshes is considered, with a spatial characteristic size that is converging to zero, the sequence of numerical solution will converge to a weak solution and, if one or more entropies are also considered, to a weak entropy solutions for each of these entropies.

The format of this paper is as follows. I start by recalling the classical notion of discrete conservation introduce by P. Lax and B. Wendroff in the early 60’s, and recall what may happen when the approximation does not exactly fit this framework. I also recall that not all scheme fits in that framework, dispites their success. Then I introduce the notion or residual distribution scheme, and give a Lax-Wendroff like theorem. Using this notion of conservation, I show that any residual distribution scheme is also a finite volume scheme, with non standard flux functions. I also show, using the same concepts, how to satisfy more than one conservation relation, and how to discretise conservative systems not written in conservation form, such as the Euler equation in primitive variables. A conclusion follows.

2 Classical setting: the Lax Wendroff theorem

The answer, or one answer to this question, has been given by Lax and Wendroff [11]. We formulate it in one spatial dimension, for simplicity, and provide references for the extension in several dimensions.

Theorem 2.1** (Lax-Wendroff).**

Consider the problem (1) for $d=1$. Consider a mesh $\{x_{j}=j\Delta x\}_{j\in\mathbb{Z}}$, the control volumes $C_{j}=]x_{j-1/2},x_{j+1/2}]$, and $\lambda>0$ Let $u_{j}^{0}$ be an approximation of

[TABLE]

Consider the numerical scheme (with $\tfrac{\Delta t}{\Delta x}=\lambda$)

[TABLE]

with the numerical flux: $\hat{f}_{j+1/2}=\hat{f}_{j+1/2}(v_{j-k},\ldots,v_{j},\ldots v_{j+k}).$ We define $u_{\Delta}$ by:

[TABLE]

Assume that:

$\hat{f}_{j+1/2}$* is consistant with ${\mathbf{f}}$: for any $v$, $\hat{f}_{j+1/2}(v,\ldots,v\ldots v)={\mathbf{f}}(v)$,* 2. 2.

$\hat{f}$* is a continuous function of its arguments,* 3. 3.

The scheme is stable: there exists a constant $C(u_{0},f)$ such that $\max\limits_{j}|u_{j}^{n}|\leq C$ for all $n$ 4. 4.

There exists a subsequence of $u_{\Delta}$ that converges to $v$ in $L^{2}(\mathbb{R}\times\mathbb{R}^{+})$

Then $v$ is a weak solution of the problem.

Proof.

Use the form of the scheme to do integration by part (Abel summation procedure)+Lebesgue’s dominated convergence thm. ∎

Corollary 2.2**.**

If $E$ is an entropy, and the numerical scheme satisfies the following inequalities

[TABLE]

where the entropy flux $\hat{g}_{l+1/2}=\hat{g}(v_{j-k},\ldots,v_{j},\ldots v_{j+k})$ is consistent with the entropy flux $\mathbf{g}$, then under the assumptions of the Lax Wendroff theorem, the function $v$ is a weak entropy solution for $E$.

Proof.

The proof is similar, using positive test functions. ∎

Extension of this results for several spatial dimensions exists, see for example [12].

Since this result, researchers have tried to improve the quality of the numerical approximation by designing more accurate, more robust flux functions, and to encapsulate in better type of time stepping approximation. But in all cases by strickly respecting the Lax Wendroff theorem. Indeed there are many excellent reasons for that, and let us show a couple of counter examples showing what occurs when this framework is violated.

2.1 Violating the flux form

Let us consider the Burgers equation,

[TABLE]

that can also be rewritten in conservation form (if one assumes that the solution is smooth),

[TABLE]

Assuming that $u$ stays positive, two reasonable approximations are

• •

Non conservation form:

[TABLE]

• •

Conservation form

[TABLE]

Both schemes can be rewritten in a very similar form. In both cases, we have

[TABLE]

where in the case (5), $a_{j+1/2}=u_{j}^{n}$ and in the case (6), $a_{j+1/2}=\tfrac{u_{j+1}+u_{j}}{2}$. Only the speed is modified.

Since we assume that $u_{j}^{n}\in[0,A]$, and since we can write in both cases

[TABLE]

we see that $u_{i}^{n+1}\in[0,A]$ if $\lambda|a_{i+1/2}|\leq\lambda A\leq 1.$ Similarly, we have easily that if in addition that $u_{j}^{n}\in[0,A]$ $\Delta x\sum_{j}|u_{j}^{n}|^{2}\leq B$ then $\Delta x\sum_{j}|u_{j}^{n+1}|^{2}\leq B$ under the same constraint. We also have that $\sum_{j\in\mathbb{Z}}|u_{j+1}^{n+1}-u_{j}^{n+1}|\leq\sum_{j\in\mathbb{Z}}|u_{j+1}^{n}-u_{j}^{n}|$ under the same constraints, i.e. the schemes are both total variation diminishing. using Helly’s theorem, we see that in both cases a subsequence in converging in $L^{2}$ to some function. Now, doing numerical simulations we see that we do not converge to the same solution. The one obtained from (5) is a weak solution, thanks to Lax Wendroff theorem. If $E$ is any entropy, denoting by $v$ the gradient of the entropy with respect to $u$, we see that for the conservative case, we get:

[TABLE]

The entropy flux (when we express the flux in term of the entropy variable) is with some abuse of language

[TABLE]

and following Tadmor, we introduce the entropy flux

[TABLE]

where $\bar{f}$ is the arithmetic average between $f_{i}$ and $f_{i+1}$. After some calculations, we get

[TABLE]

Since

[TABLE]

we get

[TABLE]

In our case, we always have $a_{j+1/2}>0$, since $u_{i}^{n}>0$ for all $i$ and $n$, we see that for $\lambda$ small enough the right hand side is negative because

[TABLE]

So the weak solution of the conservative scheme is an entropy solution for any entropy, and by uniqueness, this is the solution.

In [13], it is shown that a scheme that is stable in $L^{\infty}$ and in $BV$ and for which there exists a subsequence that converges in $L^{1}_{loc}$ to a solution, in the sense of distribution, of

[TABLE]

where $\mu$ is a locally bounded real-valued Borel measure defined on $\mathbb{R}\times\mathbb{R}^{+}$. In the case that the initial condition has a finite number of monotonicity changes, they also show that $\mu$ is concentrated on the union of the curves of discontinuity, at least for small times. It seems quite difficult to provide a more constructive description of $\mu$ and this does not prevent that $\mu\equiv 0$…! Though it seems to be seldom the case.

However, there are many cases where one would like to deal with the non conservative version of a model, see bellow for a practical example.

2.2 Change of variables

It is well know that nonlinear change of variables are not permitted. The classical example is again the Burgers equation

[TABLE]

If one sets $u=v^{3}$, this is a one-to-one change of variable, and the smooth solutions in $v$ will satisfy

[TABLE]

However the weak solutions are not similar, because the Rankine-Hugoniot relations are different: For the original Burgers equation, the velocity of a shock between the states $u^{\pm}$ is

[TABLE]

while for the second problem it is

[TABLE]

and even if $u^{\pm}=(v^{\pm})^{3}$, there is no chance that the two velocity match, in general.

2.3 Some questions

However, in many case, one would like to change variables and/or work with non conservative models. A good example is the the Euler system for fluid mechanics.

The conserved version of this system is, setting $\mathbf{u}=(\rho,\mathbf{m}=\rho\mathbf{v},E)^{T}$ where here $\mathbf{v}$ represents the velocity, $\rho$ is the density and $E=e+\tfrac{1}{2}\rho\mathbf{v}^{2}$ is the total energy, with $e$ the internal energy, the flux is

[TABLE]

The system is closed if the pressure is given in term of $\mathbf{u}$, $p=p(e,\rho)$. Defining the enthalpy $H=\tfrac{E+p}{\rho}$, under the condition that

[TABLE]

the system is hyperbolic. The speed of sound is given by

[TABLE]

If $p=(\gamma-1)e$, i.e. for a perfect gas, the speed of sound is given by the classical

[TABLE]

Another way to write the system is:

[TABLE]

and this form is interesting because the pressure is a variable that can be tracked directly from the system, and not as a consequene of the total energy equation. Unfortunately the system (7) is not in conservation form, and then starting from that relation is a priori not a good idea.

2.4 Many schemes are not naturally written in a finite volume form

To solve the problem (1), there are other methods than finite volume schemes, this does not mean they are not good. However, it seems that one of the essential feature of finite volume schemes, i.e. local conservation, is lost. In order to simplify, we consider the steady problem

[TABLE]

The domain $\Omega$ is assumed to be bounded, and regular. We assume for simplicity that its boundary is never characteristic. We also assume that it has a polygonal shape and thus any triangulation that we consider covers $\Omega$ exactly, only for simplicity. In (8b), ${\mathbf{n}}(\mathbf{x})$ is the outward unit vector at $\mathbf{x}\in\partial\Omega$ and $\mathbf{u}_{b}$ is a regular enough function. The weak formulation of (8) is: $\mathbf{u}\in L^{\infty}(\Omega)$ is a weak solution of (8) if for any $\varphi\in C^{1}_{0}(\Omega)$,

[TABLE]

where $\mathbf{\mathcal{F}}_{\mathbf{n}}$ is a flux that is almost everywhere the upwind flux:

[TABLE]

Let us show some examples, and before let us introduce some notations.

We denote by $\mathcal{E}_{h}$ the set of internal edges/faces of $\mathcal{T}_{h}$, and by $\mathcal{F}_{h}$ those contained in $\partial\Omega$. $\mathcal{K}$ stands either for an element $K$ or a face/edge $e\in\mathcal{E}_{h}\cup\mathcal{F}_{h}$. The boundary faces/edges are denoted by $\Gamma$. The mesh is assumed to be shape regular, $h_{K}$ represents the diameter of the element $K$. Similarly, if $e\in\mathcal{E}_{h}\cup\mathcal{F}_{h}$, $h_{e}$ represents its diameter.

We follow Ciarlet’s definition [14, 15] of a finite element approximation: we have a set of degrees of freedom $\Sigma_{K}$ of linear forms acting on the set $\mathbb{P}^{k}$ of polynomials of degree $k$ such that the linear mapping

[TABLE]

is one-to-one. The space $\mathbb{P}^{k}$ is spanned by the basis function $\{\varphi_{\sigma}\}_{\sigma\in\Sigma_{K}}$ defined by

[TABLE]

We have in mind either Lagrange interpolations where the degrees of freedom are associated to points in $K$, or other type of polynomials approximation such as Bézier polynomials where we will also do the same geometrical identification. Considering all the elements covering $\Omega$, the set of degrees of freedom is denoted by $\mathcal{S}$ and a generic degree of freedom by $\sigma$. We note that for any $K$,

[TABLE]

For any element $K$, $\#K$ is the number of degrees of freedom in $K$. If $\Gamma$ is a face or a boundary element, $\#\Gamma$ is also the number of degrees of freedom in $\Gamma$.

The integer $k$ is assumed to be the same for any element. We define

[TABLE]

The solution will be sought for in a space $v_{h}$ that is:

• •

Either $v_{h}=\mathcal{V}^{h}$. In that case, the elements of $v_{h}$ can be discontinuous across internal faces/edges of $\mathcal{T}_{h}$. There is no conformity requirement on the mesh.

• •

Or $v_{h}=\mathcal{V}_{h}\cap C^{0}(\Omega)$ in which case the mesh needs to be conformal.

We also need to integrate functions. This is done via quadrature formula, and the symbol $\oint$ used in volume integrals

[TABLE]

or boundary integrals

[TABLE]

means that these integrals are done via user defined numerical quadratures.

If $e\in\mathcal{E}_{h}$, represents any internal edge, i.e. $e\subset K\cap K^{+}$ for two elements $K$ and $K^{+}$, we define for any function $\psi$ the jump $[\nabla\psi]=\nabla\psi_{|K}-\nabla\psi_{|K^{+}}$. Here the choice of $K$ and $K^{+}$ is important and will become clear in each example. Similarly, \{\mathbf{v}\}=\tfrac{1}{2}\big{(}\mathbf{v}_{|K}+\mathbf{v}_{|K^{+}}\big{)}.

If $\mathbf{x}$ and $\mathbf{y}$ are two vectors of $\mathbb{R}^{q}$, for $q$ integer, $\langle\mathbf{x},\mathbf{y}\rangle$ is their scalar product. In some occasions, it can also be denoted as $\mathbf{x}\cdot\mathbf{y}$ or $\mathbf{x}^{T}\mathbf{y}$. We also use $\mathbf{x}\cdot\mathbf{y}$ when $\mathbf{x}$ is a matrix and $\mathbf{y}$ a vector: it is simply the matrix-vector multiplication.

The first example is the SUPG scheme, originally designe by T. J. Hughes and collaborators. A variant of it is used as a production code in Dassault, where shocks need to be considered …The formulation is: find $\mathbf{u}_{h}\in v_{h}=\mathcal{V}_{h}\cap C^{0}(\Omega)$ such that for any $\mathbf{w}_{h}\in v_{h}$,

[TABLE]

here $\tau_{K}$ is a strictly positive parameter, and it specific design is at the core of the method (for stability reasons). In (10), the first term of the right hand side is the Galerkin term. It is obtained by multiplying (8) by the test function $\mathbf{w}_{h}$, apply the divergence theorem and take into account the continuity of the elements of $v_{h}$ accross the edges of the mesh. The last term corresponds to the boundary conditions. The second term is a stability term: if $\mathbf{w}_{h}=\mathbf{u}_{h}$, this term is positive, but cancel if $\nabla\mathbf{f}(\mathbf{u}_{h})\cdot\nabla\mathbf{u}_{h}=0$, so that the residual property is met: the Galerkin term also vanishes (up-to the boundary terms), if $\nabla\mathbf{f}(\mathbf{u}_{h})\cdot\nabla\mathbf{u}_{h}=0$.

A variant of the SUPG scheme is obtained when one changes the stabilisation term:

[TABLE]

Here, the residual property is kept, up to boundary terms, for smooth solutions, since the jump term cancels, see [16] for details.

In these two cases, there is no clear flux formulation. One can also say the same for the discontinuous Galerkin schemes, where we look for $\mathbf{u}_{h},\mathbf{v}_{h}\in v_{h}=\mathcal{V}^{h}$ such that

[TABLE]

If a flux formulation is obvious for the averages of the conserved variables, this corresponds to only one degree of freedom. One can easily construct isomorphisms between $\mathbb{P}^{k}$ and $\mathbb{R}^{\text{dim }\mathbb{P}^{k}}$: for this, one can split the elements $K$ into $\text{dim }\mathbb{P}^{k}$ non overlapping control volumes and in general the mapping between $\mathbb{P}^{k}$ and the $\text{dim }\mathbb{P}^{k}$ dimensional vector consisting of the average of the elements of $\mathbb{P}^{k}$ on these control volume will be one-to one. This means that one can reformulate the discontinuous Galerkin scheme as a scheme actioning not on $\mathcal{V}^{h}$ but on a direct sum of copies of $\mathbb{R}^{\text{dim }\mathbb{P}^{k}}$. One one hand one would expect a natural finite volume formulation of the method, from geometry, on the other side, this formulation is not clear, though this kind of idea have already be used for hexahedral meshes in [17] or the DGSEM schemes, see for example [18].

3 A different point of view

Looking again at the schemes (11) and (10), we see that we can rewrite them as:

[TABLE]

where the element and boundary residuals $\Phi_{\sigma}^{K}(\mathbf{u}_{h})$ and $\Psi_{\sigma}^{f}(\mathbf{u}_{h})$ are defined as follows:

• •

Case of the SUPG scheme (10)

[TABLE]

with $\tau_{K}>0$.

• •

Case of the scheme (11)

[TABLE]

with $\theta_{e}>0$.

• •

In both cases,

[TABLE]

In addition, since $\sum\limits_{\sigma\in K}\varphi_{\sigma}=1$, we see that in both cases, for the internal elements $K$, we have

[TABLE]

These are not the only schemes that can be rewritten in the form (13). For example, the dG scheme (12) can be rewritten as such with

[TABLE]

where (14a) has to be slightly modified into

[TABLE]

Any finite volume also has a similar structure. Let us start with the scheme

[TABLE]

defined for the mesh $\{x_{j}\}_{j\in Z}$. The control volumes are $[x_{j-1/2},x_{j+1/2}]$ with $x_{j+1/2}=\tfrac{x_{j}+x_{j+1}}{2}$. Here $|C_{j}$ stands for the measure of $C_{j}$, $|C_{j}|=x_{j+1/2}-x_{j-1/2}$. As in [8], I introduce the quantities defined for the element $K_{j+1/2}=[x_{j},x_{j+1}]$ whatever $j\in\mathbb{Z}$ by:

[TABLE]

where $u_{h}$ is the piecewise linear interpolant of the $\{u_{l}\}_{l\in\mathbb{Z}}$ at the mesh nodes (hence we change interpretation) we see that the finite volume scheme can be rewritten as

[TABLE]

In a way, the quantity $\overrightarrow{\Phi}_{K_{j+1/2}}(u_{h})$ is the ”amount” of information sent by $K_{j+1/2}$ to the vertex $j$, while $\overleftarrow{\Phi}_{K_{j+1/2}}(u_{h})$ is the ”amount” of information sent by $K_{j+1/2}$ to the vertex $j+1$. Since the vertex $j$ belongs to $K_{j+1/2}$ and $K_{j-1/2}$, we just add the two pieces of information. Now going back to the element $K_{j+1/2}$, adding together the two pieces of informations, we get the total information, i.e.

[TABLE]

with some abuse of language.

This construction can be extended to any volume, and the key fact is that the volumes are closed, so that the integral of the outward unit normal to the control volume vanishes. Let us be more explicit, and we choose a 2D example.

Consider a conformal mesh, the vertices are $\mathbf{x}_{i}$, and the elements are generically denoted by $K$. For simplicity, we assume that $K$ is a simplex, so it is convex and we can consider its centroid. For any face, we consider again the centroid, and we connect all this in the same way as on

Here, we rephrase [19]. The notations are defined in Figure 2.

Again, we specialize ourselves to the case of triangular elements, but exactly the same arguments can be given for more general elements, provided a conformal approximation space can be constructed. This is the case for triangle elements, and we can take $k=1$.

The control volumes in this case are defined as the median cell, see figure 2. We concentrate on the approximation of $\text{div }{\mathbf{f}}$, see equation (8). Since the boundary of $C_{\sigma}$ is a closed polygon, the scaled outward normals ${\mathbf{n}}_{\gamma}$ to $\partial C_{\sigma}$ sum up to 0:

[TABLE]

where $\gamma$ is any of the segment included in $\partial C_{\sigma}$, such as $PG$ on Figure 2. Hence

[TABLE]

To make things explicit, in $K$, the internal boundaries are $PG$, $QG$ and $RG$, and those around $\sigma\equiv 1$ are $PG$ and $RG$. We set

[TABLE]

The last relation uses the consistency of the flux and the fact that $C_{\sigma}\cap K$ is a closed polygon. The quantity $\Phi_{\sigma}^{K}(\mathbf{u}_{h})$ is the normal flux on $C_{\sigma}\cap K$. If now we sum up these three quantities and get:

[TABLE]

where ${\mathbf{n}}_{j}$ is the scaled inward normal of the edge opposite to vertex $\sigma_{j}$, i.e. twice the gradient of the $\mathbb{P}^{1}$ basis function $\varphi_{\sigma_{j}}$ associated to this degree of freedom. Thus, we can reinterpret the sum as the boundary integral of the Lagrange interpolant of the flux. The finite volume scheme is then a residual distribution scheme with residual defined by (15) and a total residual defined by

[TABLE]

3.1 The residual distribution point of view

As we said, this point of view was first introduced by P.L. Roe in his seminal 1981 paper, [6], and then further developped for several dimensions in [7]. However, from an historical view point, one of the very first multidimensional residual distribution paper was written by Ni, an engineer at Bombardier, see [20].

Definition 3.1** (Residual distribution schemes).**

Considering (1), and a mesh of $\Omega$ made of simplices $K$, we will say that a scheme is a residual distribution scheme if one approximates the solution $\mathbf{u}$ of (1) by $\mathbf{u}_{h}\in v_{h}$ $v_{h}$ is the set of functions that are polynomials of degree $k$ on each element $K$ and globally continuous or not by the scheme (13) where the residuals $\Phi_{\sigma}^{K}(\mathbf{u}_{h})$ and the boundary residuals $\Psi_{\sigma}^{f}(\mathbf{u}_{h})$ satisfy the conservation relations (14a) and (14b).

One can show and is a generalisation of the classical Lax-Wendroff theorem, see [21].

Theorem 3.2**.**

Assume the family of meshes $\mathcal{T}=(\mathcal{T}_{h})$ is shape regular. We assume that the residuals $\{\Phi_{\sigma}^{\mathcal{K}}\}_{\sigma\in\mathcal{K}}$, for $\mathcal{K}$ an element or a boundary element of $\mathcal{T}_{h}$, satisfy:

• •

For any $M\in\mathbb{R}^{+}$, there exists a constant $C$ which depends only on the family of meshes $\mathcal{T}_{h}$ and $M$ such that for any $\mathbf{u}_{h}\in v_{h}$ with $||\mathbf{u}_{h}||_{\infty}\leq M$, then

[TABLE]

• •

The conservation relations (14a) and (14b).

Then if there exists a constant $C_{max}$ such that the solutions of the scheme (13) satisfy $||\mathbf{u}_{h}||_{\infty}\leq C_{max}$ and a function $\mathbf{v}\in L^{2}(\Omega)$ such that $(\mathbf{u}_{h})_{h}$ or at least a sub-sequence converges to $\mathbf{v}$ in $L^{2}(\Omega)$, then $\mathbf{v}$ is a weak solution of (8)

An immediate side result is the following result on entropy inequalities:

Proposition 3.3**.**

Let $(U,\mathbf{g})$ be a entropy-flux couple for (8) and $\hat{\mathbf{g}}_{\mathbf{n}}$ be a numerical entropy flux consistent with $\mathbf{g}\cdot{\mathbf{n}}$. Assume that the residuals satisfy: for any element $K$,

[TABLE]

Then, under the assumptions of theorem 3.2, the limit weak solution also satisfies the following entropy inequality: for any $\varphi\in C^{1}(\overline{\Omega})$, $\varphi\geq 0$,

[TABLE]

Instead of considering conservation at the level of the internal faces of the mesh, we consider it at the level of the elements. This opens new perspectives, and we will show some of them in the sequel

3.2 New examples

Using this point of view, and following the pionneering work of P.L Roe, one can define the limited Residual Distributive Schemes, see [22, 23], namely

[TABLE]

or

[TABLE]

or

[TABLE]

where the parameters $\beta_{\sigma}$ are defined to guarantee conservation,

[TABLE]

and such that (19) without the streamline term and (20) without the jump terms satisfy a discrete maximum principle. The streamline term and jump term are introduced because one can easily see that spurious modes may exist, but their role is very different compared to (10) and (11) where they are introduced to stabilize the Galerkin scheme: if formally the maximum principle is violated, experimentally the violation is extremely small, if nonexistant. See [24, 22] for more details.

A similar construction can be done starting from a discontinuous Galerkin scheme without non-linear stabilisation such as limiting, has been applied. This has been done in [25] and developped further in [26].

The non-linear stability is provided by the coefficient $\beta_{\sigma}$ which is a non-linear function of $\mathbf{u}_{h}$. Possible values of $\beta_{\sigma}$ are described in the appendix 13.

Here we consider a globally continuous approximation: $\mathbf{u}_{h}\in\mathcal{V}_{h}\cap C^{0}(\Omega)$.

Consider one element $K$. Since there is no ambiguity, the drop, for the residuals, any reference to $K$ in the following. The total residual is defined by

[TABLE]

and we assume to have monotone residuals $\{\Phi_{\sigma}^{L}\}_{\sigma\in K}$. By this we mean

[TABLE]

with $c_{\sigma\sigma^{\prime}}\geq 0$ that also satisfies

[TABLE]

It can easily be shown that the condition $c_{\sigma\sigma^{\prime}}\geq 0$ garanties that the scheme is monotone under a CFL like condition. One example is given by the Rusanov residuals:

[TABLE]

where $\bar{\mathbf{u}}$ is the arithmetic average of of the $\mathbf{u}_{\sigma}^{\prime}s$ on $K$ and $\alpha$ satisfies:

[TABLE]

Here $\#K$ is the number of degrees of freedom in $K$. Indeed, this residual can be rewritten as

[TABLE]

with

[TABLE]

Under the condition above, $c_{\sigma\sigma^{\prime}}\geq 0$ and hence we have a maximum principle.

The coefficients $\beta_{\sigma}$ introduced in the relations (19) and (20) are defined by:

[TABLE]

and can be shown to be always defined, to guaranty a local maximum principle for (19) and (20), see [22].

Remark 3.4** (About the coefficients $c_{\sigma\sigma^{\prime}}^{L}$).**

All the examples of monotone residual we are aware of are such that for linear problems, te $c_{\sigma\sigma^{\prime}}^{L}$ are independant of $\mathbf{u}_{h}$. Then one can show that for any $\sigma\in K$,

[TABLE]

This relation implies the conservation relation (14a).

4 RDS as finite volume schemes

In this section, we show how to interpret RD schemes as finite volume schemes. This amounts to defining control volumes and flux functions. We first have to define what is a flux in this context and to adapt the notion of consistency.

Let us consider any common edge or face $\Gamma$ of $K^{+}$ and $K^{-}$, two elements. Let ${\mathbf{n}}$ be the normal to $\Gamma$, see Figure 3. Depending on the context, ${\mathbf{n}}$ is a scaled normal or $||{\mathbf{n}}||=1$.

To each edge $\Gamma$ is associated a set of states $S=\{\mathbf{u}_{1},\ldots,\mathbf{u}_{l})$. A flux $\hat{\mathbf{f}}_{\mathbf{n}}(S)$ between $K^{+}$ and $K^{-}$ has to satisfy

[TABLE]

The consistency property that stands for the consistency is that if all the states are identical in an element, then each of the residuals vanishes. Hence, we define a multidimensional flux as follows: A multidimensional flux

[TABLE]

is consistent if, when $\mathbf{u}_{1}=\mathbf{u}_{2}=\ldots=\mathbf{u}_{N}=\mathbf{u}$ then

[TABLE]

The results of this section apply to any finite element method but also to discontinuous Galerkin methods. There is no need for exact evaluation of integral formula (surface or boundary), so that these results apply to schemes as they are implemented.

Let $K$ is a polytope contained in $\mathbb{R}^{d}$ with degrees of freedoms on the boundary of $K$. The set $\mathcal{S}$ is the set of degrees of freedom in $K$. We consider a triangulation $\mathcal{T}_{K}$ of $K$ whose vertices are exactly the elements of $\mathcal{S}$. Choosing an orientation of $K$, it is propagated on $\mathcal{T}_{K}$: the edges are oriented. This is illustrated in figure 4 for a $\mathbb{P}^{2}$ triangle and a $\mathbb{Q}^{2}$ quad

The problem is to find normals ${\mathbf{n}}_{\sigma\sigma^{\prime}}$ with ${\mathbf{n}}_{\sigma\sigma^{\prime}}=-{\mathbf{n}}_{\sigma^{\prime}\sigma}$ and quantities $\hat{\mathbf{f}}_{{\mathbf{n}}_{\sigma\sigma^{\prime}}}:=\hat{\mathbf{f}}_{\sigma,\sigma^{\prime}}$ for any edge $[\sigma,\sigma^{\prime}]$ of $\mathcal{T}_{K}$ such that:

[TABLE]

Note that (22b) implies the conservation relation

[TABLE]

In short, we will consider

[TABLE]

but other examples can be considered provided the consistency (22c) relation holds true: the choice is a bit arbitrary, provided (22c) holds true.

Any edge $[\sigma,\sigma^{\prime}]$ is either direct or, if not, $[\sigma^{\prime},\sigma]$ is direct. Because of (22b), we only need to know $\hat{\mathbf{f}}_{\sigma,\sigma^{\prime}}$ for direct edges. Thus we introduce the notation $\hat{\mathbf{f}}_{\{\sigma,\sigma^{\prime}\}}$ for the flux assigned to the direct edge whose extremities are $\sigma$ and $\sigma^{\prime}$. We can rewrite (22a) as, for any $\sigma\in\mathcal{S}$,

[TABLE]

with

[TABLE]

$\mathcal{E}^{+}$ represents the set of direct edges.

Hence the problem is to find a vector $\hat{\mathbf{f}}=(\hat{\mathbf{f}}_{\{\sigma,\sigma^{\prime}\}})_{\{\sigma,\sigma^{\prime}\}\text{ direct edges}}$ such that

[TABLE]

where $\Psi=(\Psi_{\sigma})_{\sigma\in\mathcal{S}}$ and $A_{\sigma\sigma^{\prime}}=\varepsilon_{\sigma,\sigma^{\prime}}$.

We have the following lemma which shows the existence of a solution.

Lemma 4.1**.**

For any couple $\{\Phi_{\sigma}\}_{\sigma\in\mathcal{S}}$ and $\{\hat{\mathbf{f}}_{\sigma}^{b}\}_{\sigma\in\mathcal{S}}$ satisfying the condition (22c), there exists numerical flux functions $\hat{\mathbf{f}}_{\sigma,\sigma^{\prime}}$ that satisfy (22). Recalling that the matrix of the Laplacian of the graph is $L=AA^{T}$, we have

The rank of $L$ is $|\mathcal{S}|-1$ and its image is \big{(}\text{span}\{\mathbf{1}\})^{\bot}. We still denote the inverse of $L$ on \big{(}\text{span}\{\mathbf{1}\})^{\bot} by $L^{-1}$, 2. 2.

With the previous notations, a solution is

[TABLE]

The proof can be found in [1]. The computation of $A^{T}L^{-1}$ is easy in practice: since $L$ is symetric, the range of $L$ is orthogonal to its kernel, spanned by $x_{0}=(1,\ldots,1)^{T}$. Hence, for any $\lambda\neq 0$, the matrix $L+\lambda\frac{x_{0}\otimes x_{0}}{||x_{0}||^{2}}$ is invertible, and the matrix written as $L^{-1}$ with some abuse of language is $(L+\lambda\frac{x_{0}\otimes x_{0}}{||x_{0}||^{2}})^{-1}-\frac{x_{0}\otimes x_{0}}{\lambda||x_{0}||^{2}}$: the computation can be done by any standard matrix inversion package.

If in addition, the boundary flux satisfy for any $\mathbf{u}$

[TABLE]

as for example (22d), then this set of flux are consistent and the normals ${\mathbf{n}}_{\sigma,\sigma^{\prime}}$ are given by

[TABLE]

We can state:

Proposition 4.2**.**

If the residuals $(\Phi_{\sigma})_{\sigma\in K}$ and the boundary fluxes $(\hat{\mathbf{f}}_{\sigma}^{b})_{\sigma\in K}$ satisfy (22c) and (25), then we can find a set of consistent flux $(\hat{\mathbf{f}}_{\sigma,\sigma^{\prime}})_{[\sigma,\sigma^{\prime}]}$ satisfying (22). They are given by (24). In addition, for a constant state,

[TABLE]

for the normals defined by (26).

This also defines the control volumes since we know their normals.

We can state a couple of general remarks:

Remark 4.3**.**

The flux $\hat{\mathbf{f}}_{\sigma,\sigma^{\prime}}$ depend on the $\Psi_{\sigma}$ and not directly on the $\hat{\mathbf{f}}_{\sigma}^{b}$. We can design the fluxes independently of the boundary flux, and their consistency directly comes from the consistency of the boundary fluxes. 2. 2.

The residuals depends on more than 2 arguments. For stabilized finite element methods, or the non linear stable residual distribution schemes, see e.g. **[27, 28, 22]**, the residuals depend on all the states on $K$. Thus the formula (24) shows that the flux depends on more than two states in contrast to the 1D case. In the finite volume case however, the support of the flux function is generally larger than the three states of $K$, think for example of an ENO/WENO method, or a simpler MUSCL one. 3. 3.

The formula (24) are influenced by the form of the total residual (16). We show in the next paragraph how this can be generalized. 4. 4.

The formula (24) make no assumption on the approximation space $v_{h}$: they are valid for continuous and discontinuous approximations. The structure of the approximation space appears only in the total residual. 5. 5.

Quadrature formula: all the relations we use are obtained by quadrature formula. This means that the integration does not need to be exact.

To end this paragraph, let us give one example. We consider a triangle and a quadratic approximation, see figure 5

The one can easily show that

[TABLE]

[TABLE]

Then we choose the boundary flux:

[TABLE]

and get:

[TABLE]

The normals are given by:

[TABLE]

[TABLE]

5 Satisfaction of constraints

In this section, we want to show that this notion of conservation at the level of elements can also be useful to construct, from a known scheme, a new one that satisfy additional constraints. More explicitly, let us consider the two different problems:

• •

Assume there is a function $\mathbf{u}\in\mathbb{R}^{m}\mapsto\mathbf{U}(\mathbf{u})\in\mathbb{R}^{m\times m}$ such that the (possibly) non conservative non linear PDE

[TABLE]

An example, taken from fluid mechanics, is

[TABLE]

If we take

[TABLE]

then we recover the conservative form of the Euler equations.

• •

Assume formally that the conservative system

[TABLE]

satisfies an additional conservation relation,

[TABLE]

where $E$ is a function of the state $\mathbf{u}$. Assume in addition, that we satisfy the second system by some algebraic manipulations, for example there exists a mapping $\mathbf{u}\in\mathbb{R}^{m}\mapsto\mathbf{v}(\mathbf{u})\in\mathbb{R}^{m}$ such that

[TABLE]

In other words,

[TABLE]

An example is the entropy $E$.

We first show how, starting from a known scheme, we can modify it so that the new scheme will satisfy the additional constraints. Since both problems are unsteady, we first explain how we discretise them in space and time, then we show how to modify the schemes.

5.1 Discretisation in space and time

Let us consider (1) and follow [29, 2]. If one wants to discretize this problem starting from the schemes of the type (13) where the residuals are given by (10), (11) or (18), and if we want to keep the automatic consistency with the original PDE that is provided by the residual formulation, we are led to method where there is a mass matrix. This is also true with (12), but here the mass matrix is block diagonal (so easy to invert)contrarily to the other cases where it is only sparse. In addition, in the case (18), it is very unclear that this mass matrix (which is formal in that case) is …invertible: the trick introduced in [29] is precisely done to avoid to invert the mass matrix, to the price of adding some dissipative term. In [2], this trick was reinterpreted as a Defect Correction approach and then generalised to any order. In detail the formal scheme will write as: find $\mathbf{u}_{h}\in v_{h}$ such that for any test function $v_{h}$m, we have

[TABLE]

where $a$ is the form defined in (10), (11) or (18) that we can write as

[TABLE]

, $J$ is a possible jump term (such as in (11) and (20) and $v_{h}$ a test function that is

• •

constant by element $K$ for (19) and (20),

• •

any element of $v_{h}$ for (11),

• •

of the form $w_{h}+h_{K}\tau_{K}\nabla_{u}{\mathbf{f}}(\mathbf{u}_{h})\cdot\nabla w_{h}$ for (10)

Then

[TABLE]

A priori, the idea is to start from

[TABLE]

and to use some quadrature formula. Note that the test function do not depend on time. This amount to subdivide the interval $[t_{n},t_{n+1}]$ with sub-timesteps $t_{n,0}=t_{n}<t_{n,1}<\ldots<t_{n,p-1}<t_{n,p}=t_{n+1}$ and to approximate the relation (30) at the sub-timesteps. If $\mathbf{u}_{h}$ is the vector $(\mathbf{u}_{h}^{n},\mathbf{u}_{h}^{n,1},\ldots,\mathbf{u}_{h}^{n,p-1},\mathbf{u}_{h}^{n+1})$ with $\mathbf{u}_{h}^{n,l}\approx\mathbf{u}(t_{n,l})$, we write the approximation as

[TABLE]

For example, the Crank-Nicholson method leads to

[TABLE]

and

[TABLE]

In [2], knowing $\mathbf{u}_{h}^{n}$, we compute $\mathbf{u}_{h}^{n+1}$ by the following algorithm (we just show the simplest version).

• •

Set $\mathbf{u}^{(0)}=\mathbf{u}_{h}^{n}$,

• •

for $p=1,\ldots,N$, compute $\mathbf{u}^{(p+1)}$ by solving

[TABLE]

where

[TABLE]

The forms $a$ are written in term of spatial residuals of the type (10), (11) or (19) and (20).

In practice $N$ is equal to the number of sub-time steps, and the accuracy is not spoiled if some technical condition described in [2] are met. They are met in particular is the degrees of freedoms $\sigma$ are the control points of the Bézier polynomials of degree $N$.

5.2 Corrected scheme: the example (27)

We have a set of residuals $\{\Phi_{\sigma}^{K}\}_{\sigma\in K}$ that satisfy the ’conservation’ relations:

[TABLE]

The surface integral is approximated by a quadrature formula such that the error is $|K|O(h^{k+1})$, where $k$ is the polynomial degree.

If there exists some average of $\mathbf{U}(\mathbf{u})$, say $\overline{\mathbf{U}}(\mathbf{u}^{(p+1)},\mathbf{u}^{(p)})$ such that

[TABLE]

then the conservation constraints are recovered if for any $K$, we have:

[TABLE]

where $\Phi_{\sigma}^{K}(\mathbf{u}^{n},\mathbf{u}^{n,1,(p)},\ldots,\mathbf{u}^{n+1,(p)})$ are the residual that are needed to define $A$ in (31) and $\mathcal{I}({\mathbf{f}}(\mathbf{u}^{n})\cdot{\mathbf{n}},{\mathbf{f}}(\mathbf{u}^{n,1,(p)})\cdot{\mathbf{n}},\ldots,{\mathbf{f}}(\mathbf{u}^{n+1,(p)})\cdot{\mathbf{n}})$ is a weighted average of the normal flux at the sub-timesteps after the $p$-th iteration, this translate the way (31) discretise (30). In the case of the Crank-Nicholson method, this is simply an arithmetic average between the states at $t_{n}$ and $t_{n+1}$. Of course these relations are in general not satisfied by the original residuals, and we modify them by adding a correction

[TABLE]

We only describe what happens for the elements $K$.

The correction $r_{\sigma}^{K}$ must satisfy

[TABLE]

We have one (vectorial) relations, and at least 2 unknowns (since an element has at least two degrees of freedom). So by a simple linear algebra argument, there might be a solution. The problem is to find the average in such a way that the explicit nature of the scheme is kept, and second to compute the correction: this is a case by case situation. We show with the example of fluid mechanics how this can be achieved and show it on a specific example: second order of accuracy.

Using a piecewise linear interpolation in time, the conservation relations are, for $p=0,1$,

[TABLE]

Here we haver written $\mathbf{u}^{(p)}$ for $\mathbf{u}^{n+1,(p)}$ since there is no ambiguity. Let us write explicitly the differences $\mathbf{w}^{(p)}-\mathbf{w}^{(n)}$ in term of the primitive variable. In the sequel, for any quantity $f$, $\Delta f=f^{(p+1)}-f^{(p)}$.

Since the density is a conservative and a primitive variable, there is nothing to write.

[TABLE]

i.e. in matrix form:

[TABLE]

The matrix is lower triangular, so that the scheme can be kept explicit: we first compute the density: we know $\rho^{(p+1)},\rho^{(p)},\mathbf{u}^{(p)},e^{(p)}$. Then we compute the velocity $\mathbf{u}^{(p+1)}$: we know $\rho^{(p)},\rho^{(p+1)},\mathbf{u}^{(p)},\mathbf{u}^{(p+1)},e^{(p)}$, and then the internal energy. The last question is how to evaluate the corrections. There is no correction on the density component of the residuals. For the velocity component we write

[TABLE]

In this relation, the right hand side can be explicitly computed from what is known: we have one vectorial relation and as many unknown as degrees of freedom in $K$. A priori, there is no reason why to weight differently the degrees of freedom, so we assume

[TABLE]

and then

[TABLE]

The same method is used for the energy correction, and we get

[TABLE]

where the $\Phi_{\rho\mathbf{u}}^{K}$ are corrected residuals and $N_{K}$ the number of degrees of freedom in the element.

5.3 Corrected scheme: the example 28

In that case, the same kind of trick is used, with a small difference: the corrections $r_{\sigma}$ must not destroy the initial conservation law, and we must have a relation of the type (33). In order to illustrate, we give the example of fluid mechanics with the mathematical entropy $S=p\rho^{-\gamma+1}$. In the sequel, we denote by $\kappa=\gamma-1$ and recall that $p=\kappa e$.

We can write that

[TABLE]

where

[TABLE]

Assuming we have a initial scheme for the variables $(\rho,\mathbf{u},p)$, combining the techniques of the previous example and the formula bellow, we see that the residual on the velocity should corrected as previously, while the residuals on the pressure should be corrected such that

[TABLE]

This leads to a linear system of the type

[TABLE]

where $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ are computable quantities. Since there is more than two degrees of freedom in $K$, if the density is not uniform, we can compute the corrections. If the densities are the same ($\rho^{(p)}_{\sigma}=\rho^{(p)}$), then one has to choose initial residuals such that \big{(}\rho^{(p)}\big{)}^{-\kappa}\mathcal{E}_{1}=\mathcal{E}_{2}.

Remark 5.1**.**

Using this entropy, we keep the explicit nature of the scheme. With other entropies, this is less clear.

6 Conclusions and perspectives

In this paper, we have discussed how the conservation property writes in the residual distribution framework. Instead of looking at what happens at the cell interfaces, we look at the element contributions. Using this concept, it is possible to reformulate most if not all the known schemes as finite volume schemes, with explicit formula for the flux. Using this notion, it is possible to construct schemes that start from a non conservative formulation of a conservative systems, or to enforce more than one conservation relation. We show the principles, and provide some examples. Other examples are possible, see for example [4].