R Abgrall

M Ricchiuto

High order methods for CFD

We provide a review of high order methods for CFD. Besides recalling some classical methods, we show a framework allowing on one hand to see and work with these methods under a different light, and on the other to provide a different path to construct numerical methods for flow equations. In particular, we focus on Residual Based techniques, and Residual Distribution methods, as a framework to construct schemes of arbitrary order. The somewhat classical second-order multidimensional upwind fluctuationsplitting/residual-distribution schemes are reviewed in the chapter by Deconinck and Ricchiuto.

Introduction

In this chapter we discuss with some extend high order methods for advection dominated problems. Typical examples are the advection diffusion equation (with large Peclet number), the Euler equations, the Navier Stokes equations, the Shallow water equations and many problems in geophysical flow, to mention just a few. The list is immense.

In this kind of problems, one has to deal with contransting constraints. First the solution must be accurate. Wherever the solution is smooth, the truncation error must scale as h r+1 where h is the typical size of the mesh elements, and r is an integer. However, it is also well known that in many case, the solution is not globally smooth, and that it may admit local very large gradients. For example, these may be shock waves (or slip lines) for the Euler equations, and boundary layers for the Navier-Stokes equations at very high Reynolds number. Other effects may come into play, as e.g. dispersion, as in some geophysical and environnemental applications (wave propagation, capillary flows).

In this chapter, we address these issues and give several examples of successful modern high order methods. The literature on this topic has exploded since the mid 90's, and it is not possible to give an exhaustive view of all what has been achieved, so we had to make choices, which, of course, are biased by our own work.

We will start by reviewing to some extent possibly the two most popular methods today: WENO finite volume schemes, and the Discontinuous Galerkin (DG) methods. We give some details, in particular indicate the main principles. However, we do not discuss all the issues related to these methods. In particular, set aside the choice of the approximation of the viscous terms: this can easily be found in the many papers that have appeared on the topics or in monographs such as [START_REF] Abgrall | Chapters for the Handbook of Numerical Methods for Hyperbolic Problems[END_REF][START_REF] Di Pietro | Mathematical aspects of Discontinuous Galerkin methods[END_REF]. Our main focus is on a less successful approach known as the residual distribution method. As also discussed in the companion chapter by Deconinck and Ricchiuto, these schemes share a lot with continuous Finite Element Methods, such as the Streamline diffusion method, but also embed properties typical of the finite volume method: a lot of emphasis is put on L ∞ stability constraints, allowingto avoid spurious oscillations at discontinuities. In this contribution, however, we extend this analysis, showing how the Residual-Based philosophy underlying these schemes provides a framework which allows to embed most (or all) other arbitrary order methods, and work with them under a different light, thus providing mode insight in these methods, and perhaps new, alternative contructions. Of course it also provides a setting to construct different arbitrary order schemes, and we will review the main challenges encountered when doing that, as well as some of the solutions proposed so far to overcome these challenges.

This chapter is organized into 4 sections. The first one gives a review of the WENO and DG methods. The second section develop in details the Residual distribution method, both from an historical perspective and its most recent achievements. The third section provides several applications, both for compressible flows and aerodynamics, and for some geophysical flows. A conclusion follows. We hope that the bibliography is rich enough to cover and complete all the topics we have mentioned in the text.

A review of existing methods

We start by considering the following problem:

∂w ∂t + div f (w) -div fv(w, ∇w) = 0 (1a)
defined in Ω ⊂ R d with d = 1, 2, 3, and w : Ω × R + → D ⊂ R m . We need to set up an initial condition w(x, 0) = u0(x),

x ∈ Ω (1b)

and boundary conditions on ∂Ω. In (1a), the flux f = (f1, . . . , f d ) is assumed to be defined on D and to be smooth enough. The viscous flux fv is assumed to satisfy similar hypothesis (however see below).

Here, we are mostly interested in fluid mechanics problems where the Navier Stokes equation is the canonical example. In that case, the state variable w needs to satisfy additional constraints: the density and the internal energy need to be positive; this is what is expressed in saying that w ∈ D.

In this section, we review some of the high order methods that have been developed in the recent years. The research activity in this field has been very intense over the last years, so that it is impossible to make an exhaustive review of what has been done. We had to make choices, this choices are biased.

The second choice we have made, at least for this section, is to deal with a simplified problem. Instead of (1a), we will deal with ∂w ∂t + div f (w) = 0 [START_REF] Di Pietro | Mathematical aspects of Discontinuous Galerkin methods[END_REF] with the initial condition (1b) and relevant boundary conditions. In the case of fluid mechanics, the simplest system of this kind is the Euler system. There, the state variable is

w = (ρ, ρu, E) T
where ρ is the density, u is the fluid velocity, and E is the total energy, i.e. the sum of the internal energy e and the quinetic energy. The flux is given by: f = (ρu, ρu ⊗ u + p Id d , (E + p)u) T .

To close the system, we need to define the pressure p = p(ρ, e). General assumptions on this function can be found in [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], a typical example is that of perfect gaz:

p = (γ -1)e.
It is well known that only weak solutions of (2) can be considered because there is no hope to have regular solutions in general. Hence, we need to consider weak solutions, i.e. measurable functions in L ∞ (Ω × R + ) m ∩ L 1 (Ω × R + ) m such that for any compactly supported regular test function ϕ ∈ C 1 0 (Ω × R + ) m , we have

Ω×R + ∂ϕ ∂t (x, t) • w(x, t)dxdt + Ω×R + ∇ϕ(x, t) • f (w)dxdt - Ω ϕ(x, 0) • w0(x)dx = 0. (3)
Of course the initial condition needs to be also in L ∞ (Ω) m ∩ L 1 (Ω) m . Note we have not taken into account the boundary conditions. This is a complex problem (for a rigorous treatment and also from a practical point of view), we refer to [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF] for systems. It is also well known that the definition of weak solution is not enough. Even in the scalar case, this does not guaranty the uniqueness of solution and some selection mechanism is needed. For the system case, the solution is way more complex. To go towards this, one classical consider an entropy, i.e. a strictly convex function S defined on D such that there exists an entropy flux G such that:

∇uG = ∇wS • ∇wf .
An entropy solution should satisfy, in the sens of distribution, the following inequality:

∂S ∂t + div G ≤ 0. (4)
In the case of the Euler equations, the (mathematical) entropy is given by S = -ρs where s is the (physical) entropy. The entropy flux is G = Su. The reader may consult [START_REF] Harten | On the symmetric form of conservation laws with entropy[END_REF][START_REF] Th | Finite element formulation for Computational Fluid Dynamics : I symmetric forms of the compressible Euler and Navier Stokes equations and the secoond law of thermodynamics[END_REF] for further considerations. The form (3) is the origin of all possible forms of numerical approximation of the system [START_REF] Abgrall | Chapters for the Handbook of Numerical Methods for Hyperbolic Problems[END_REF]. The first thing to do is to approximate the domain Ω. For the sake of simplicity, we assume here that Ω is polygonal. Then we discretize Ω using meshes. With the vocabulary of unstructured meshes, we consider a tessaletion T h . The domain Ω is Ω = ∪K∈T h K.

As usual, we assume that the elements K are non overlapping. The elements K will be triangles or quadrangles in 2D, tetrahedrons, hexahedrons, pyramid, etc in 3D, or may have more complicated forms. All depends on the choices made for approximating the solution and the choices of test functions ϕ.

• Finite volume schemes: one considers that w is constant in each cell,

wK (t) ≈ 1 |K| K w(x, t)dx,
and the test functions are also constant. On can nevertheless get high order accuracy of the averaged value. This is the topic of section 2.1.1.

• Continuous finite elements. Here we assume a globally continuous approximation w h of w. Typicaly, for any element, w h|K is a polynomial of degree k. Because of the continuity requirement, this imposes constraints on the mesh: the intersection of two elements is either empty, or reduced to a (complete) face or they are identical. The mesh is said to be conformal. The elements need also, in general, be simplices because of the polynomial approximation. These methods are sketched in 2.1.3.

• Discontinuous Galerkin methods. Here the continuity requirement is dropped. This enables a lot of freedom: the mesh needs not be conformal, the element can be general, so that mesh refinement becomes simple in principle. These methods are sketched in 2.1.2.

In the rest of this section, we first consider the spatial approximation (hence using a semi-discrete formulation), and then we discuss a bit the temporal approximation.

Space discretizations

ENO and WENO

Here, we consider the finite volume formulation of (3). The states are described by {wK (t)}K∈T h . Starting from (3), we get d dt K w(x, t)dx

+ ∂K f (w) • nd∂K = 0
Here n is the outward unit normal to the boundary ∂K or K. This can be obtained from (3) by first regularizing via mollification ϕε the characteristic function of K, and taking the limit when ε → 0, we see that ∇ϕε → n.

Since K is polygonal, denoting by e a generic face/edge of K, we see that an approximation of (3) is

d dt wK (t) + 1 |K| e∈∂K e f (w h ) • nde = 0. (5)
This relation has not yet a meaning since u h is discontinuous accross edges. In the normal direction to e, we need to solve the following Riemann problem:

∂w ∂t + ∂f (w) • n ∂n = 0
with the initial condition:

w(x, 0) = wK if x • n < 0 w K + else.
Here, u K + is the state on the other side of e. This problem is solved either exactly or in an approximated way. A meaning of the edge integral is given thanks to the use of numerical flux f (uK , u K + ), see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] for an extensive discussion about numerical flux and Riemann solvers. Hence, the finite volume method in its simplest form is:

d dt wK (t) + 1 |K| e∈∂K e
f (w h|K , w h|K + , n)de = 0.

Note that the edges integrals are evaluated via quadrature formula. As this, only first order accuracy can be achieved. Formal high order accuracy can be obtained by using the MUSCL method due to van Leer [START_REF] Van Leer | Toward the ultimate conservative difference scheme V. a second order sequel to Godunov's method[END_REF]. The idea is to consider a polynomial reconstruction of degree p, R(u h ), within each cell K and to replace [START_REF] Th | Finite element formulation for Computational Fluid Dynamics : I symmetric forms of the compressible Euler and Navier Stokes equations and the secoond law of thermodynamics[END_REF] by

d dt wK (t) + 1 |K| e∈∂K e f (R(w h )K , R(w h ) K + , n)de = 0. ( 7
)
and the quadrature formula need to be of sufficient order, typicaly exact for polynomials of degree p.

The design of a reconstruction operator is a research field by itself since one wants to avoid the Gibbs phenomena where the solution becomes steep or discontinuous. After the seminal work of van Leer, a very large literature has been devoted to this problem. A large part of it is about Total Variation Diminishing schemes (see Sweby's paper [START_REF] Sweby | High resolution schemes using limiter for hyperbolic conservation laws[END_REF]), but it is rather difficult to reach higher than second order accuracy (see [START_REF] Chakravarthy | A new class of high accuracy tvd schemes for hyperbolic conservation laws[END_REF] for an attempt in 1D), and it can be shown that a TVD scheme in more than one dimension, even in the scalar case, is at most first order accurate, see [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF]. This negative result has motivated to look for criteria that are less strict than the TVD one, and the most successfull method is probably the Essentialy Non Oscillatory method, orginally due to Harten and co-workers [START_REF] Harten | Uniformly high-order accurate nonoscillatory schemes I[END_REF][START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes[END_REF], then refined by Shu and co-workers [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes, ii[END_REF]. Extension to unstructured meshes can be found in [START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: Analysis and implementation[END_REF]. Better stability properties are obtained by the so-called Weighted Essentially Non Oscillatory technique (WENO), see [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF] for the original reference, and further refined by Shu and co-workers, see [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF] for a review.

The principle can be explained assuming a regular mesh, in one dimension. Extension to 2D, and more general meshes can be found in [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF][START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF][START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using WENO limiters II: unstructured meshes[END_REF] for example. Taking a mesh {xj}j∈z, with xj = j∆x, we first define four approximation of a smooth function at x i+1/2 = x i +x i+1 2 by:

• using the stencil S (1) = {xi-2, xi-1, xi}, we have u

(1)

i+1/2 = 3 8 ui-2 -5 4 ui-1 + 15 8 ui = u(x i+1/2 ) + O(∆x 3 ) • Using the stencil With S (2) = {xi-1, xi, xi+1}, we get: u (2) = -1 8 ui-1 + 3 4 ui + 3 8 ui+1 = u(x i+1/2 ) + O(∆x 3 ) • With the stencil, With S (3) = {xi, xi+1, xi+2}, we have u (3) (x i+1/2 ) = 3 8 ui+ 3 4 ui+1-1 8 ui+2+O(∆x 3 ) • Last, with S = {xi-2, xi-1, xi, xi+1, xi+2}, we obtain u (4) i+1/2 = 3 128 ui-2 - 5 32 ui-1 + 45 64 ui + 15 32 ui+1 - 5 128 ui+2 = u(x i+1/2 ) + O(∆x 5 ).
Then we notice that u (4)

i+1/2 = γ1u (1)
i+1/2 + γ2u

(2)

i+1/2 + γ3u (3)
i+1/2 with γ1 = 1 16 , γ2 = 5 8 and γ3 = 5 16 . Note that γ1 +γ2 +γ3 = 1. This enables to approximate u(x i+1/2 ) by u i+1/2 = w1u

(1)

i+1/2 +w2u (2) i+1/2 +w3u (3) i+1/2
with wj ≈ γj and wj ≈ 0 if a discontinuity exists in S (j) . In order to achieve this, we define wj = wj w1 + w2 + w3 with wj = γj (ε + βj) 2 and the smoothner indicater βj is:

βj = 2 l=1 ∆x 2l-1 x i+1/2 x i-1/2 d l dx l pj(x) 2 dx.
Compared to the method we are going to discuss now, for a given formal accuracy they clearly need the lowest possible storage. The price to pay for this is that the computational stencil is quite large. In the case of an irregular mesh, many precautions need to be taken in order to effectively reach the formal accuracy. In the case of unstructured mesh, their extension is possible but very technical.

Discontinuous Galerkin methods.

Formulation and basic properties. This class of method has original been designed by W.H. Reed and T.R. Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF], the first analysis was done by Lesaint and Raviart [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF] and further refined by [START_REF] Johnson | An analysis of the discontinous galerkin method for a scalar hyperbolic equation[END_REF]. The references [START_REF] Chavent | The local projection p 0 p 1 -discontinuous galerkin finite element method for scalar conservation laws[END_REF] and mostly [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element methid for conservation laws[END_REF][START_REF] Cockburn | The local Runge-Kutta local projection discontinuous finite element method for scalar conservation laws[END_REF] and their sequel paved the way to the success of DG methods for hyperbolic problems and the Navier Stokes equations (among many other applications). The reference [START_REF] Cockburn | Discontinuous Galerkin Methods-Theory, computation and applications[END_REF] represents the state-of-the-art in the early 2000, the reference [START_REF] Di Pietro | Mathematical aspects of Discontinuous Galerkin methods[END_REF] is a more mathematical presentation of the theory, it also contains a lot of informations on how to approximate parabolic problems in that framework. It is impossible to give a complete survey of this field because the number of papers grows exponentially. Again, we will sketch the method for purely hyperbolic problems, and refer to the reference section of [START_REF] Di Pietro | Mathematical aspects of Discontinuous Galerkin methods[END_REF] for more informations on how to approximate the second order terms.

Again, we start from (3). We consider a tessellation of the computational domain Ω, like in the finite volume method, but here we look for solutions that are polynomial of degree r ≥ 0 in each element. More precisely, we want to approximate the solution in V h defined by:

V h = {w h ∈ (L ∞ (Ω)) m ∩ (L 1 (Ω)) m , for any element K, (w h ) |K ∈ (P r (K)) m }
No continuity requirement is needed, and more over the degree r may depend on the element. Then we apply the weak formulation, taking as test function any ϕ ∈ V h : for any element K,

K ∂ϕ ∂t • w h (x, t)dx - K ∇ϕ • f (w h (x, t))dx + ∂K ϕ • (f (w h (x, t)) • n)d∂K = 0.
However, as for the Finite Volume method, this formulation is meaningless because u h appearing in the boundary integral is in general multivalued, so the flux term cannot be given a meaning. The idea contained in [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element methid for conservation laws[END_REF][START_REF] Cockburn | The local Runge-Kutta local projection discontinuous finite element method for scalar conservation laws[END_REF] is, as for the finite volume method, to introduce a numerical flux f . The DG (semi discrete) formulation is thus: find w h ∈ V h such that for any K and any

ϕ h ∈ V h , K ∂ϕ ∂t • w h (x, t)dx - K ∇ϕ • f (w h (x, t))dx + ∂K ϕ • f (w h|K , w h|K + , n)d∂K = 0 (8)
In any element K, w h ∈ P r (K). This vector space is spanned by a finite set of polynomial functions:

w h = R j=0
w (j) ϕj, so that we arrive at the following form of the semi-discrete scheme: for any K,

MK d dt WK + F (w h|K ) = 0
where the mass matrix is (MK )ij = K ϕiϕjdx is clearly invertible. This is a block diagonal matrix, hence its inversion (needed for time discretization) is rather straightforward. The vector F is defined by its components:

Fj = - K ∇ϕj • f (w h (x, t))dx + ∂K ϕj • f (w h|K , w h|K + , n)d∂K.
The choice of the degree of freedom, i.e. the choice of the basis function, is an issue by itself. The choices are made depending whether to favor a geometrical interpretation (Lagrange basis), to facilitate the change of polynomial degree within elements (in the case of degree adaptivity), or if the element shape is or not completely general (for example in the case of Lagrangian hydrodynamics, [START_REF] Vilar | A discontinuous galerkin discretization for solving the two-dimensional gas dynamics equations written under total lagrangian formulation on general unstructured grids[END_REF]), etc.

Non linear stability. One can show, in the scalar case, that a global entropy inequality can be easily derived, see [START_REF] Jiang | On a cell entropy for Discontinuous Galerkin Methods[END_REF]. In the scalar case, a natural entropy is U (u) = u 2 2 : this is a convex function, and an entropy g = (gx, gy) flux satisfies

uf x = g x , u f y = g y .
We see that gx = ufx -u fxdu, gy = ufy -u fydu. In the following, we set hx = u fxdu and hy = u fydu. We wish to establish an inequality of the type: for any K,

K ∂U (u h ) ∂t dx + ∂K ĝ • ndx ≤ 0. (9)
Here ĝ • n is an entropy flux, i.e. a numerical flux consistant with g. This inequality simply states that we have a local L 2 energy bound. For any v h ,

K ∂u h ∂t v h dx - K f (u h )∇v h xdx + ∂K v h f (u h|K , u h|K -, n)d∂K = 0.
Then we choose v h = u h , so that

1 2 K ∂(u h ) 2 ∂t dx - K f (u h )∇u h dx + ∂K u h f (u h|K , u h|K -, n)d∂K = 0 (10) We get 1 2 K ∂(u h ) 2 ∂t dx + ∂K Ĝ(u h|K , u h|K -, n)d∂K + AK = 0 with Ĝj+1/2 = u h|K + u h|K - 2 f (u h|K , u h|K -, n) - 1 2 g(u h|K ) + g(u h|K -) • n
This flux is consistant with g • n and we also have set

AK = [u h|K ,u h|K -] f (u h|K , u h|K -, n) -f (v) • n dv.
Using the mean value theorem, we see that

AK = (u h|K -u h|K -) f (u h|K , u h|K -, n) -f (ξ) • n
for a suitable ξ between u h|K and u h|K -. If f is an E-Scheme (see [START_REF] Osher | Rieman solvers, the entropy condition and difference approximations[END_REF]), i.e. if for any ξ between u and v, f (u, v, n) -f (ξ) • n (u -v) ≤ 0, we see that ( 9) holds true for any E-scheme. Typical example are the exact Godunov solver, the Rusanov scheme:

f (u, v, n) = 1 2 f (u) • n + f (v) • n + α(u -v) for α ≥ max ξ∈[min(u,v),max(u,v)] f (ξ) • n .
The case of systems is of course more complex. One possible solution could be to rewrite the system (1) in term of entropy variables: v = ∇wS(w) where S is a (mathematical) entropy, in the following form:

w ∂w ∂t + div h(v) = 0
The change of variables v → w is one-to-one and does not affect the Rankine-Hugoniot relations. Instead of approximating the state variable w by piecewise polynomials, we can approximate the entropy with polynomials of degree r + 1 and define the approximation state as

V h = {v ∈ (L 1 (Ω)) m ∩ (L ∞ (Ω)) m , v ∈ P r+1 (K)}
and look for w(v), with v ∈ V h such that for any ϕ ∈ V h , for any

K K ϕ ∂w(v) ∂t - K ∇ϕ • h(v)dx + ∂K ϕ h(vK , v K -, n)d∂K = 0
Clearly, if h is an E-scheme, we have the entropy inequality:

K ∂S(v) ∂t + ∂K G(vK , v K -, n)d∂K ≤ 0
for a suitable consistant entropy flux G.

Controlling spurious oscillations. Another and somewhat related issue is to control the Gibbs phenomena: when the numerical solution develops steep gradients, either because the mesh resolution is not sufficient or because discontinuities are appearing, spurious oscillations will appear. One of the fundamental question is how to control them as automatically as possible. One of the solution is to get inspired by what has been done for finite volume schemes, taking into account the negative result of Goodman and LeVeque [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF]. In order to describe what has been achieved, let us turn back to the 1D case for the simplification of exposure. In that case, the scheme reduces to:

K j ϕ du h dt - K j ϕ f (u h )dx + ϕ(x j+1/2 ) f (u - h,j+1/2 , u + h,j+1/2 ) -ϕ(x j-1/2 ) f (u + h,j-1/2 , u - h,j-1/2 ) = 0 ( 11
)
where

Kj = (x j-1/2 , x j+1/2 ), xj = x j-1/2 ,x j+1/2 2 , ∆j = x j+1/2 -x j-1/2 and u ± h,j+1/2 = u h (xj ± δ j u 2 ) = u h (xj) ± δ ± j u.
One of the main remark that enables to understand the behavior of methods is what is called Harten's lemma. Instead of considering the semi-discrete case, let us use full discretized form, we will come back to this in section 2.2 Assume we have a sequence {u n j } j∈z,n∈N that satisfies (λ > 0):

u n+1 j = u n j -λ C j+1/2 (u n j+1 -u n j ) -D j-1/2 (u n j -u n j-1 ) . ( 12
)
If for any j ∈ Z, we have

C j+1/2 ≥ 0, D j+1/2 ≥ 0, λ C j+1/2 + D j+1/2 ≤ 1 (13)
then the sequence satisfies a L ∞ , a L 1 and a TVD bound, where the total variation of u = {uj}j∈z is:

T V (u) = j∈z |uj+1 -uj|.
There is no reason why the arguments u ± j are such that the sequence defined by [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF] are such that the

I j-1 I j I j+1 u 0,j-1 + u + j-1 u 0,j-1 u 0,j + u - j u 0,j+1 -u - j+1 u 0,j+1 u 0,j u 0,j -u - j 1
Figure 1: Geometrical representation of a P 1 approximation. In the element K j , u h = u 0,j + δu j (xx j ).

x j is the centroid of K j , δ j is its lengh. We have set δu ± j = u h (x j ± ∆j 2 )u 0,j .

conditions [START_REF] Harten | Uniformly high-order accurate nonoscillatory schemes I[END_REF] are true. To do so, one technique is to introduce a limiter. The simplest is the generalized minmod limiter: m(a1, a2, . . . , am) = s min(|a1|, |a2|, . . . , |am|) if s = sign(a1) = . . . sign(am) 0 else [START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes[END_REF] The arguments in the flux f [START_REF] Chakravarthy | A new class of high accuracy tvd schemes for hyperbolic conservation laws[END_REF] are modified as follows. We replace δ ± j by

δ ± j ) mod = m(δ ± j , ∆+u0,j, ∆-u0,j).
There is a huge litterature on this topic. One may quote among others: [START_REF] Cockburn | TVD Runge Kutta local projection discontinuous Galerkin Finite element method for conservation laws II : general framework[END_REF][START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF][START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF][START_REF] Krivodonova | Shock detection and limiting with discontinuous galerkin methods for hyperbolic conservation laws[END_REF][START_REF] Qiu | A comparison of trouble cell indicators for runge-kutta discontinuous galerkin method using weno limiters[END_REF][START_REF] Li | Hybrid weighted essentially non-oscillatory schemes with different indicators[END_REF]. Other kind of polynomial representation can also be used, such as Hermite approximation, see [START_REF] Qiu | Hermite weno schemes and their application as limiters for runge-kutta discontinuous galerkin method ii: Two dimensional case[END_REF][START_REF] Balsara | A sub-cell based indicator for troubled zones in RKDG schemes and a novel class of hybrid RKDG+HWENO schemes[END_REF] and their relation to limiting. Another approach is to combine WENO limiters and Discontinuous Galerkin method, see [START_REF] Dumbser | Arbitrary high order PNPM schemes on unstructured meshes for the compressible navier-stokes equations[END_REF][START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using WENO limiters II: unstructured meshes[END_REF].

Stabilized Continous FEM

Another approach for spatial approximation is to consider the following trial space:

V h = {w h ∈ (L ∞ (Ω)) m ∩ (L 1 (Ω)) m ∩ (C 0 (Ω)) m , for any element K, (w h ) |K ∈ (P r (K)) m }.
The fundamental difference is that we now require continuity. If one use elements of V h as test function, i.e. look for solutions w h such that: for any ϕ

∈ V h , Ω ϕ • ∂w h ∂t dx + a(w h , ϕ) = 0 (15a) where a(w, v) = - Ω ∇v • f (w h )dx + BC(w, ϕ). (15b)
The problem amounts to solving (with clear notations)

M dw h dt + A(w h ) = 0, Mij = Ω ϕiϕjdx.
The mass matrix M is also invertible. It is a sparse matrix but it is not block diagonal, contrarily to the Discontinuous Galerkin method. In (15b), the operator BC describes the approximation of the weakly enforced boundary conditions. We do not describe it because it depends on the nature of the boundary conditions: it is problem dependent. The method ( 15) is known have stability difficulties, so it is better to add to the Galerkin variational form a stabilization operator. There are several forms of this stabilization operator: the stream line operator [START_REF] Th | Finite element formulation for Computational Fluid Dynamics : I symmetric forms of the compressible Euler and Navier Stokes equations and the secoond law of thermodynamics[END_REF][START_REF] Johnson | An analysis of the discontinous galerkin method for a scalar hyperbolic equation[END_REF] and a jump operator [START_REF] Burman | Consistent supg-method for transient transport problems: Stability and convergence[END_REF]. Instead of solving [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF], we solve:

Ω ϕ • ∂w h ∂t + a(w h , ϕ) + aS(w h , ϕ) = 0 (16a)
where aS is a stabilisation operator.

In the case of the streamline operator, the choice is [START_REF] Th | Finite element formulation for Computational Fluid Dynamics : I symmetric forms of the compressible Euler and Navier Stokes equations and the secoond law of thermodynamics[END_REF]:

aS(w h , ϕ) = K hK K ∇wf (w h ) • ∇ϕ TK ∂w h ∂t + ∇uf (w h ) • ∇w h dx = 0, (16b)
where hK represents the diameter of K and TK ≥ 0 is a stabilisation parameter (or matrix).

In the case of the jump operator, we take [START_REF] Burman | Consistent supg-method for transient transport problems: Stability and convergence[END_REF] aS(w h , ϕ)

= internal edges γeh 2 e e ∇wf (w h ) [∇w h ][∇ϕ h ] (16c)
where γe ≥ 0, he is the measure of the edge e. The choice of the stabilisation operator is done so that if the exact solution also satisfies [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes, ii[END_REF]. Note that the structure of the mass matrix is affected in the case (16b), hence its invertibility is less obvious. In the case of (16c), the mass matrix is not changed, but the compactness of the computational stencil is slightly affected. Since these methods share a lot of similarities with the residual distribution methods (indeed, they can be seen as a particular case), we postpone the discussion later in the text.

Temporal discretizations

There are several standard ways of approximating in time, depending on how we look at time with respect to space. Either they are two unrelated parameters, so that one first approximate in space and then in time thanks to the method of lines. Or we consider the equation

∂w ∂t + div f (w) = 0
as a space-time divergence applied to the flux (w, f (w)). Here, we focus on the explicit method of lines. A typical example is the well knows method of line. After having discretized in space, we have to discretize a problem of the form:

∂w ∂t = L(w). (17)
Depending on the stiffness of the problem, more general speaking, of the properties we are looking for, one may consider explicit or implicit scheme. A very popular method is the so called Strong Stability

Preserving technique [START_REF] Gottlieb | Strong Stability-Preserving High-Order Time Discretisation Methods[END_REF]: If the Euler operator w → v = w -∆tL(w) preserves the L ∞ norm or the L 1 norm or the TVD semi norm for ∆t ≤ ∆t0, then the SSP Runge-Kutta method it is built on will also have the same property under a condition of the type ∆t ≤ C∆t0. In the cases we are interested in, ∆t0 is defined by mean of a CFL-type condition, and is approximation dependent. The explicit SSP RK schemes are written as Runge Kutta schemes

u (0) = u n u (i) = i-1 k=0 α i,k u (k) + ∆tβ i,k L(u (k) ) , i = 1, . . . , m u n+1 = u (m)
where the α i,k and β i,k are all positive. By consistency, i-1 k=1 α i,k = 1, s that the intermediate stages can be written as convex combination of the Euler operator. The integer m is the number of stages. Examples of such SSP RK methods are the following:

• Second order in time and C = 1 (no degradation of the time step)

u (1) = u n + ∆tL(u n ) u n+1 = 1 2 u n + 1 2 u (1) + ∆tL(u (1) )
• Third order in time and C = 1.

u (1) = u n + ∆tL(u n ) u (2) = 3 4 u n + 1 4 u (1) + ∆tL(u (1) ) u n+1 = 1 3 u n + 2 3 u (2) + ∆tL(u (2) )
In these examples, the number of stages is equal to the order of the scheme. It can be shown [START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF] that there exists no 4th order SSP RK scheme with four stages. Spiteri and Ruuth [START_REF] Steven | Two barriers on strong-stability-preserving time discretization methods[END_REF]developed fourth order methods with m = 5, 6, 7 and 8 stages: for example u (1) = un + 0.391752226571890 ∆tL(u n ), u (2) = 0.444370493651235 u n + 0.555629506348765 u (1) + 0.368410593050371 ∆tL( (1) ), u (3) = 0.620101851488403 u n + 0.379898148511597 u (2) + 0.251891774271694 ∆tL( (2) ), u (4) = 0.178079954393132 u n + 0.821920045606868 u (3) + 0.544974750228521 ∆tL( (3) ), u n+1 = 0.517231671970585 u (2) + 0.096059710526147 u (3) + 0.063692468666290 ∆tL( (3) ) +0.386708617503269 u (4) + 0.226007483236906 ∆tL( (4) )

for which C = 1.508. A rather complete discussion can be found in [START_REF] Gottlieb | Strong Stability-Preserving High-Order Time Discretisation Methods[END_REF][START_REF] Gottlieb | On high order strong stability preserving runge-kutta and multi step time discretizations[END_REF][START_REF] Steven | Two barriers on strong-stability-preserving time discretization methods[END_REF]. Error estimates for explicit Runge Kutta time stepping can be found in [START_REF] Burman | Explicit Runge-Kutta schemes and finite elements with symetric stabilization for first-order linear PDE systems[END_REF][START_REF] Zhang | Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws[END_REF][START_REF] Zhang | Stability analysis and a priori error estimates for the third order explicit runge-kutta discontinuous galerkin method for scalar conservation laws[END_REF][START_REF] Meng | Optimal error estimates for dicontinuous galerkin methods based on upwind-biased fluxes for linear hyperbolic equations[END_REF].

A different setting : residual distribution

We now consider the framework known today as residual distribution. Its roots can be found the seminal work of P.L. Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF] on fluctuation splitting, and in all the contributions of the 90s on wave decomposition, hyperbolic elliptic splitting, and multidimensional upwind methods [START_REF] Roe | Linear advection schemes on triangular meshes[END_REF][START_REF] Roe | optimum" upwind advection on a triangular mesh[END_REF][START_REF] Roe | Multidimensional upwinding : motivation and concepts[END_REF][START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF][START_REF] Sidilkover | Unification of some advection schemes in two dimensions[END_REF][START_REF] Struijs | Progress on multidimensional upwind Euler solvers for unstructured grids[END_REF][START_REF] Deconinck | A multidimensional generalization of Roe's difference splitter for the Euler equations[END_REF][START_REF] Nishikawa | A third-order fluctuation splitting scheme that preserves potential flow[END_REF], and see also [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF].

We present it here as a general framework to study and unify the "more classical" approaches recalled in section 2, while giving some additional flexibility to construct new "non classical" discretizations.

Steady hyperbolic problems

Consider the scalar steady state advection equation

a • ∇u = 0 on Ω ⊂ R 2 (18)
where ∇ • a = 0, and with boundary conditions To find a numerical approximation of the solution of ( 18)-( 19), on a tessellation of the spatial domain Ω h , we use a generalization of the fluctuation splitting strategy put forward by P.L. Roe. In particular, we start by considering u h , a continuous nodal finite element approximation of the solution

∂Ω ( a • n) -(g -u) = ∂Ω - a • n(g -u) = 0 (19)
u h = i∈Ω h ϕiui = K∈Ω h uiϕi K (20)
For a given degree of freedom i of the continuous collocated finite element expansion, let Ki be the set of elements sharing i as a node, and similarly, let Fi be the set of mesh faces sharing i. Given an initial guess for the degrees of freedom, we proceed as follows 1. For all elements K, compute the fluctuation/residual

φ K = K a • ∇u h | K dx ≈ - K ∂tu h dx (21)
2. For all elements K, distribute the fluctuation to the three nodes of K. Let φ K j denote the amount of fluctuation sent to node j ∈ K, then the conservation/consistency requirement is

j=j K j=1 φ K j = φ K (22)
3. For all nodes i ∈ Ω h , assemble signals from surrounding elements and evolve toward steady state by some iterative procedure such as e.g.

u n+1 i = u n i -ωi K∈K i φ K i ( 23
)
The method described by ( 21)-( 22)-( 23) aims at providing a solution to the discrete algebraic system

K∈K i φ K i = 0 , ∀ i ∈ Ω h (24)
As formulated, it does not include boundary conditions (BCs). The most general way to introduce them is to consider, for any face f ∈ ∂Ω h the face fluctuations.

φ f = f a • n(g * h -u h ) df (25)
where, to embed the compatibility condition implicit in [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF], we have introduced the numerical flux :

(g * a) • n = a • n 1 + sign( a • n) 2 u + 1 -sign( a • n) 2 g
Face fluctuations can be split to the degrees of freedom j ∈ f by means of distributed residuals φ f j such that

j f j=1 φ f j = φ f (26)
Finally, the complete discrete fluctuation splitting/residual distribution steady equations read

K∈K i φ K i + f ∈F i ∂Ω h φ f i = 0 (27)

Accuracy conditions

The first formulation of these schemes, on linear triangular elements, relied for the construction of second order discretizations on the so-called linearity preservation property [START_REF] Roe | Linear advection schemes on triangular meshes[END_REF][START_REF] Paillere | Compact cell vertex convection schemes on unstructured meshes[END_REF] defined as follows.

Definition 3.1 (Linearity preservation). Let {β K j }j∈K be a set of distribution coefficients uniformly bounded with respect to h, u h , φ K , and with respect to the data of the problem ( a, boundary data, etc), and verifying the consistency property

j K j=1 β K j = 1 (28)
A Linearity Preserving scheme is one for which

φ K i = β K i φ K (29)
Proposition 3.2. Linearity preserving schemes are second order accurate.

The simple property stated in definition 3.1 and in proposition 3.2 has been known since the late 80s, but it has taken more than a decade to be formally understood. A more general characterisation of the accuracy of these schemes, due to [START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF] and generalized in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Abgrall | An example of high order residual distribution scheme using non-lagrange elements[END_REF][START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF], is the following. Definition 3.3 (Truncation error and accuracy). Let ψ be smooth function, ψ ∈ C r+1 (Ω). Let Ω h be an unstructured grid composed of non-overlapping elements. On the generic element K ∈ Ω h consider the r-th degree continuous polynomial approximation [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF]. Let in particular ψ h = j∈K ψjϕj be the r-th degree polynomial approximation of type (20) of ψ, the values ψj being obtained by Galerkin projection. Consider now an exact smooth function u ∈ H r+1 verifying (18)-( 19) in a classical sense : a • ∇u = 0 in Ω, and u = g on ∂Ω -. Let u h be its polynomial approximation of degree r of type [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF], obtained by Galerkin projection. Let now φ K j (u h ) and φ f j (u h ) the value of the split residuals obtained when replacing the nodal values of the solution obtained with the scheme by the values uj of the Galerkin projection of u. We define the integral truncation error (u h , ψ)

(u h , ψ) = j∈Ω h ψj    K∈K j φ K j (u h ) + f ∈F j φ f j (u h )    = K∈Ω h j K j=1 ψjφ K j (u h ) + f ∈∂Ω h j f j=1 ψjφ f j (u h ) (30)
We shall say that a scheme is r + 1 order accurate if it verifies the truncation error estimate

| (u h , ψ)| ≤ C(Ω h ) h r+1
The following general characterization is possible.

Proposition 3.4. In d spatial dimensions, a sufficient condition for scheme [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element methid for conservation laws[END_REF] to be r + 1 order accurate in the sense of definition 3.3 is to simultaneously have

φ K i (u h ) ≤ CΩ h h r+d ∀K ∈ Ω h , ∀ i ∈ K φ f i (u h ) ≤ C ∂Ω h h r+d-1 ∀f ∈ ∂Ω h , ∀ i ∈ f (31)
The proof of this property is omitted for brevity. The interested reader can refer to [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Abgrall | An example of high order residual distribution scheme using non-lagrange elements[END_REF][START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] for details. The importance of this characterization is that it allows to provide some design conditions. To see this, one must first recall that for the solution u of definition 3.3, and for its Galerkin projection u h on the r-degree finite element polynomial space, we can use classical approximation results [START_REF] Ciarlet | General lagrange and hermite interpolation in R n with applications to finite element methods[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF] to show that a

• ∇u h -a • ∇u ≤ Cuh r in Ω h
and, provided that the boundary ∂Ω is also smooth enough, provided that ∂Ω h is a high order polynomial rendering of the exact boundary [START_REF] Abgrall | An example of high order residual distribution scheme using non-lagrange elements[END_REF][START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF], we also have1

(g * a • n) h -(g * a • n) ≤ Cu,nh r+1
where the norms used are standard L norms, such as the L 2 or the max norm, with no derivatives involved.

With the regularity hypotheses made on the mesh, we also have that

|K| = O(h d ) and |f | = O(h d-1
), for and K and any f . This, and the fact that a • ∇u = 0 and g * -u = 0 for the exact solution, leads to the conclusion that: a sufficient condition for a scheme to be r + 1 order accurate in the sense of definition 3.3 is that we can find for any K ∈ Ω h and for any f ∈ ∂Ω h sets of uniformly bounded test functions ω K i and ω f i , such that

j K j=1 ω K i = 1 and j f j=1 ω f i = 1
and that the distribution can be obtained as

φ K i (u h ) = K ω K i a • ∇u h dx and φ f i (u h ) = f ω f i a • n (g * h -u h ) df (32)
Clearly, the linearity preserving schemes of definition 3.1 are obtained as the particular case in which the test functions are constant within each element ! As we will see immediately, this consistency analysis, applies trivially to classical continuous Galerkin discretizations, as well as to their stabilized counterparts, as well as to discontinuous Galerkin methods.

Stability and convergence

The above consistency analysis gives conditions under which that if convergence with respect to the mesh parameter h is obtained, r + 1 convergence rates are expected w.r.t. h for a r-th degree polynomial approximation, and in correspondence of sufficiently smooth solutions. The missing piece of information is : how do we make sure that convergence is indeed achieved ? A classical finite element convergence analysis would need two main ingredients [START_REF] Ern | Theory and practice of finite elements[END_REF] : a consistency estimate, which we have provided, and a stability condition, which we have not. If we could provide a stability statement which ensures e.g. that ∀u h in our approximation space

| K j K j=1 ujφ K j (u h ) + f j f j=1 ujφ f j (u h )| ≥ C u h 2 , with 0 < C < ∞ , (33)
then using more or less classical arguments [START_REF] Ern | Theory and practice of finite elements[END_REF], we could infer the existence of the discrete solution, and derive more rigorous estimates for the error associated to this solution.

Unfortunately, to this day residual distribution schemes lack a framework for stability analysis. Some weaker results showing the decay of the solution energy (L 2 norm) during iterations [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] have been shown in several works [START_REF] Barth | An energy look at the N scheme[END_REF][START_REF] Abgrall | Residual distribution schemes for conservation laws via adaptive quadrature[END_REF][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. These conditions are however not sufficient to say more on the discrete solution.

On the other hand, we are able to rule out some schemes as the following property shows in two space dimensions2 . Proposition 3.5 (Fall of the βφ paradigm, 2d advection). Consider the solution of a • ∇u = 0 in two space dimensions, with a constant, and with ∂Ω a collection of straight sides. Any scheme of the form 0 =

K∈K i β K i φ K + f ∈F i ∂Ω h β f i φ f
cannot be freed of high frequency spurious modes whatever the form of

β K i , if K is a P k Lagrange triangle with k > 2 and if K is a Q k Lagrange quadrilateral ∀k ≥ 1.
Proof. For all elements considered, we explicitly show one spurious mode exact solution of the discrete problem with homogeneous boundary conditions. This mode can be added to any grid function without the scheme detecting its presence, thus preserving this unphysical perturbation.

First recall that for homogeneous boundary conditions and using the hypothesis on ∂Ω

φ K = ∂K a • nu h = f ∈∂K a • n f u h df and φ f = - f 1 -sign( a • n) 2 a • n u h df = - 1 -sign( a • n) 2 a • n f u h df
so we focus on the approximation of the integrals of u h over element faces. Denote by let the number of freedom on each face f ∈ ∂K be C f + 2. We consider the mode defined by uj

= 1 if j is a vertex, otherwise on each f ∈ ∂K we set ∀j = v uj = - 2 C f f ϕv df f ϕj df
having denoted with v one of the two vertices forming face f . The mode is compatible with the continuity of the representation, and with the adoption of hybrid meshes. For P k triangles with k ≥ 3 and Q k elements with k ≥ 2, the value of the solution at nodes within the elements remains arbitrary. For this mode, one easily checks that φ K = 0 , ∀K, and that

φ f = 0 , ∀ f ∈ ∂Ω h .
The only remaining element is the Q 1 quadrilateral which is easily checked to suffer from the checkerboard spurious mode in which u oscillates between -1 and 1 on every face.

The important consequence of proposition 3.5 is that we have to start looking for schemes exploiting the sub-elemental variation of the discrete solution. A well known example of such a scheme, perfectly fitting the framework presented, is the SUPG scheme of T.J.R. Hughes and collaborators [START_REF] Hughes | Streamline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF][START_REF] Franca | Stabilized finite element methods : I. application to the advective-diffusive model[END_REF][START_REF] Hughes | Multiscale and Stabilized Methods[END_REF] obtained by setting

φ K i = K ϕi a • ∇u h dx + Streamline Dissipation K a • ∇ϕi τK a • ∇u h dx and φ f i = f ϕi a • n (g * h -u h ) df (34)
Stability results for the SUPG scheme can be obtained in the classical sense discussed in the beginning of this section (see for example [START_REF] Szepessy | Convergence of a shock-capturing streamline diffusion finite element method for a scalar conservation law in two space dimensions[END_REF][START_REF] Johnson | On the convergence of shock capturing streamline diffusion finite element methods for hyperbolic conservation laws[END_REF][START_REF] Johnson | On the convergence of a finite element method for a nonlinear hyperbolic conservation law[END_REF][START_REF] Bochev | Stability of the supg finite element method for transient advection-diffusion problems[END_REF][START_REF] Burman | Consistent supg-method for transient transport problems: Stability and convergence[END_REF][START_REF] Hughes | Multiscale and Stabilized Methods[END_REF] and references therein), and are based on the positive-semidefinite nature of the bi-linear form associated to the Streamline-dissipation term.

Other examples of schemes allowing to overcome the flaw of proposition 3.5 will be given in the following. In general, guidelines to construct such methods can be obtained by considering the convergence of iteration [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF]. In the simplest setting of scalar advection, if we recast this iteration as the following update for the array of degrees of freedom U

U n+1 = U n -ω(A h U n -F )
convergence requires that for some r < 1 and for all V

(Id -ωA h )V 2 ≤ r V 2
which can be developed into

V t A h V ≥ 1 -r 2ω V 2 + ω 2 A h V 2 ≥ C h V 2 ≥ 0
leading back to a condition of type [START_REF] Cockburn | TVD Runge Kutta local projection discontinuous Galerkin Finite element method for conservation laws II : general framework[END_REF], and to the necessary (albeit not sufficient) condition

V t A h V ≥ 0 ( 35
)
which we will use in the following.

Embedding a discrete maximum principle

Monotonicity of the numerical solution is retained by the so-called local positivity contraints for the distribution. This property is related to positive coefficient theory which has replaced the TVD theory to construct high order schemes [START_REF] Goodman | On the accuracy of stable schemes for 2d scalar conservation laws[END_REF][START_REF] Spekreijse | Multigrid solution of monotone second-order discretizations of hyperbolic conservation laws[END_REF][START_REF] Barth | Numerical methods for conservation laws on structured and unstructured meshes[END_REF].

Definition 3.6 (Positive scheme). A (locally) positive scheme is one for which

φ K i = j∈K j =i cij(ui -uj) , cij ≥ 0 ∀j ∈ K (36)
Positivity is the key to the construction of non-oscillatory schemes [START_REF] Roe | Linear advection schemes on triangular meshes[END_REF][START_REF] Paillere | Compact cell vertex convection schemes on unstructured meshes[END_REF] :

min j∈K i u n j ≤ u n+1 i ≤ max j∈K i u n j ∀ i ∈ Ω h
under the following condition

min i∈Ω h ωi K∈K i j∈K j =i cij ≤ 1
Proof. The proof follows from the positivity of the cijs and time step restriction, and from

u n+1 i = 1 -ωi K∈K i j∈K j =i cij u n i + ωi j∈K i j =i K∈K i K j ciju n j
This characterization can be generalized to fully consistent time dependent discretizations as we will show later.

A general framework : relation with classical discretization approaches

The formalism introduced in the previous sections for the scalar advection equation can be easily generalized to (systems of) steady nonlinear conservation laws of the form div f (w) = 0 on Ω ⊂ R 2 [START_REF] Qiu | A comparison of trouble cell indicators for runge-kutta discontinuous galerkin method using weno limiters[END_REF] with the appropriate boundary conditions on ∂Ω. We now look for a solution satisfying

K∈K i φ φ φ K i + f ∈F i ∂Ω h φ φ φ f i = 0 ( 38
)
where in every element

K i∈K φ φ φ K i = ∂K f (w h ) • n d∂K = K div f (w h )dx (39)
while on a boundary face f we have

i∈f φ φ φ f i = f f (w h , g, n) -f (w h ) • n df (40)
where the numerical flux f (w h , g, n) accounts for the boundary conditions. This framework defines a sort of super class of methods, which allows to embed, and has relations, with all the discretization approaches introduced in section 2.

Continuous FEM as residual distribution

The simplest example is perhaps that of the stabilized finite elements discussed in section 2.1.3. In particular, given a continuous collocated finite element expansion for which we can write V h = span{ϕi}i∈Ω h , then the (un-stabilized) continuous Galerkin method is obtained simply by setting

φ φ φ K i = K ϕi div f (w h ) dx , φ φ φ f i = f ϕi f (w h , g, n) -f (w h ) • n df
For nodal finite elements, the relation i∈K ϕi = 1 ensures that consistency is satisfied. There is, however, a slight catch which is worth underlying. The relation between the last definitions, and the consistency condition [START_REF] Qiu | Hermite weno schemes and their application as limiters for runge-kutta discontinuous galerkin method ii: Two dimensional case[END_REF] requires that exact integration is performed, w.r.t. the assumed polynomial variation of w h and the definition of the nonlinear flux f . This is required to go from the integral of the flux divergence, to the boundary integral [START_REF] Qiu | Hermite weno schemes and their application as limiters for runge-kutta discontinuous galerkin method ii: Two dimensional case[END_REF], so that the variational formulation (15b) is recovered. In practice, exact quadrature is never used. The practical way to handle this issue, is to introduce a high order polynomial representation of the flux f h . Based on the accuracy of a given quadrature formula, we can uniquely identify the polynomial degree of such an expansion, built starting from the reconstructed values of w h at a sufficient number of flux evaluation points, exactly as done in the so called quadrature free approaches used in Discontinuous Galerkin [START_REF] Atkins | Quadrature-Free Implementation of Discontinuous Galerkin Method for Hyperbolic Equations[END_REF][START_REF] Lockard | Efficient implementations of the quadrature free Discontinuous Galerkin method[END_REF], and in the most recent flux reconstruction methods (cf [START_REF] Huynh | High-order methods for computational fluid dynamics: A brief review of compact differential formulations on unstructured grids[END_REF] and references therein). Based on the exact evaluation of the integrals of this polynomial flux, we can re-formulate our consistency conditions as being

i∈K φ φ φ K i = ∂K f h (w h ) • n d∂K = K div f h (w h )dx (41)
and

i∈f φ φ φ f i = f f h (w h , g, n) -f h (w h ) • n df (42)
The choice of the polynomial degree has to respect at least some accuracy constraints which are easily deduced from the analysis of section 3.1.1. In particular, for this analysis to apply in the nonlinear case, one must ensure that for a given smooth flux f , and for a r-th degree finite element approximation

div f h (w h ) -div f K ≤ Cuh r , f h • n -f • n f ≤ Cu,nh r+1
This requires the flux polynomial to be of degree r f ≥ r. With this modification, the (un-stabilised) continuous Galerkin approximation will read

φ φ φ K i = K ϕi div f h (w h ) dx , φ φ φ f i = f ϕi f h (w h , g, n) -f h (w h ) • n df
while a stabilised variant is readily obtained by including in φ φ φ K i the streamline dissipation term (cf. section 2.1.3), leading to a SUPG distribution :

φ φ φ K i = K ϕi div f h (w h ) dx + K ∇wf (w h ) • ∇ϕi TK ∇wf (w h ) • ∇w h dx ( 43
)
where the relation i∈K ϕi = 1 allows to show that ( 42)-( 41) are met.

Although different in spirit, the edge-stabilized schemes discussed e.g. in [START_REF] Burman | Explicit runge-kutta schemes and finite elements with symmetric stabilization for first-order linear pde systems[END_REF][START_REF] Burman | Interior penalty continuous and discontinuous finite element approximations of hyperbolic equations[END_REF] can be also recast in the formalism above by setting

φ φ φ K i = K ϕi div f h (w h ) dx + ∂K γ ∂K (w h )[∇w h ] • [∇ϕi] ( 44
)
where again consistency is a consequence of the partition of unity property, while one ca easily demonstrate that the accuracy conditions are met provided that γ ∂K (w h ) = O(h 2 ), as in (16c).

Finite volume vs residual distribution : local conservation and continuous FEM

We recall here the analogy between residual distribution node centered finite volume schemes on median dual cells. With the notation of section 2.1.1, and with reference to figure 3, we can write the semi-discrete evolution equation for wi, the average of w on cell Ci, as

|Ci| dwi dt = - j f ij f R(w h )i, R(w h )j, nij = - K∈K i j∈K j =i f K ij f R(w h )i, R(w h )j, nij = - K∈K i j∈K j =i f K ij •nij (45) Local conservation is equivalent now to the condition f K ij • nij + f K ji • nji = 0 ( 46
)
Since Ci is a closed polygon, we also have K∈K i j∈K j =i nij = 0, which allows to recast [START_REF] Steven | Two barriers on strong-stability-preserving time discretization methods[END_REF] as

|Ci| dwi dt = - K∈K i j∈K j =i ( f K ij -fi) • nij (47)
If we now set

φ φ φ K i = j∈K j =i ( f K ij -fi) • nij (48)
we find that the local conservation property [START_REF] Gottlieb | On high order strong stability preserving runge-kutta and multi step time discretizations[END_REF] implies

j∈K φ φ φ K i = 1 2 j∈K fj • nj = ∂K f h • n df = φ φ φ K with f h = j ϕjfj.
This shows that for any given higher order finite volume discretization we may define a residual distribution method consistent with a second order polynomial flux approximation. While this was known for a certain time, the reverse is not. In particular, given a definition of the split residuals {φ φ φ K j }j∈K we may ask if there exist a definition of consistent fluxes expressing local conservation over the median dual cell for the residual distribution method. If we require these fluxes to satisfy [START_REF] Zhang | Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws[END_REF], the we may write the system

f K ij • nij + f K il • n il = φ φ φ K i - 1 2 fi • ni f K ji • nji + f K jl • n jl = φ φ φ K j - 1 2 fj • nj f K li • n li + f K lj • n lj = φ φ φ K l - 1 2 f l • n l (49)
using local conservation, and setting Ψ Ψ

Ψ K i = φ φ φ K i -fi • ni/2, we obtain a linear system for ( f K ij • nij, f K jl • n jl , f K li • n li )   1 0 -1 -1 1 0 0 -1 1     f K ij • nij f K jl • n jl f K li • n li   =    Ψ Ψ Ψ K i Ψ Ψ Ψ K j Ψ Ψ Ψ K l    (50)
The associated matrix has rank 2. We can easily find particular solutions setting to zero one of the unknowns and solving the resulting sub system. Averaging out the three particular solutions obtained (for symmetry) er en up with the following multidimensional numerical fluxes w.r.t. which the residual distribution scheme is locally conservative on the median dual cell :

f K ij • nij = f K ij • nij(wi, wj, w l ) = 1 3 (Ψ Ψ Ψ K i -Ψ Ψ Ψ K j ) = 1 3 (φ φ φ K i -φ φ φ K j ) - 1 6 (fi • ni -fj • nj) f K jl • n jl = f K jl • n jl (wi, wj, w l ) = 1 3 (Ψ Ψ Ψ K j -Ψ Ψ Ψ K l ) = 1 3 (φ φ φ K j -φ φ φ K l ) - 1 6 (fj • nj -f l • n l ) f K li • n li = f K li • n li (wi, wj, w l ) = 1 3 (Ψ Ψ Ψ K l -Ψ Ψ Ψ K i ) = 1 3 (φ φ φ K l -φ φ φ K i ) - 1 6 (f l • n l -fi • ni) (51)
These are three states multidimensional numerical fluxes. Consistency can be formulated as

f K ij • nij(w, w, w) = f (w) • nj -ni 6
as for a constant state over the element we always have φ φ φ K i = 0 ∀ i. Other standard properties of numerical fluxes, e.g. Lipschitz continuity, are inherited from the properties of the physical flux f , and of the split residuals φ φ φ K i . A similar construction can be repeated for schemes written on high order finite elements. To understand how this works, we start from the P 2 case. We consider the set-up defined of figure 4 : the element is split first into 4 sub-triangles K1, K2, K3, and K4. From this sub-triangulation, we can construct a dual mesh as in the P 1 case. The dual mesh is a collection of cells Cj whose intersection with an element K defines 6 sub-zones, represented by dashed lines in the figure. The notation used in this case is similar to the one used before : in the sub-triangle Ki, we denote by n K i ij the normal to the portion of the face separating the median dual cells Ci and Cj, and by f K i ij • n K i ij the corresponding local numerical flux. We can now write down the finite volume equations for each control cell Cj, and then proceed as in the P 1 case to determine contributions associated to each sub-element Ki. To relate these sub-elemental residuals to the P distributed residuals, we sum for each node the contribution from the sub-elements to [START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF] having set

F K i = C i ∂K f (w h ) • n dΓ
with w h the P finite element solution, and with the obvious relation

K∈K i F K i = 0
due the the continuity of the flux. If now we set, to simplify the notation

f14 := f K 1 14 • n K 1 14 , f61 := f K 1 61 • n K 1 61 , f64 := f K 4 64 • n K 4 64 -f K 1 46 • n K 1 46 f42 := f K 2 42 • n K 2 42 , f25 := f K 2 25 • n K 2 25 , f45 := f K 4 45 • n K 4 45 -f K 2 54 • n K 2 54 f53 := f K 3 53 • n K 3 53 , f36 := f K 3 36 • n K 3 36 , f56 := f K 4 56 • n K 4 56 -f K 3 65 • n K 3 65 and Ψ Ψ Ψ K i = φ φ φ K i -F K i , we obtain f14 -f61 = Ψ Ψ Ψ K 1 -f42 + f25 = Ψ Ψ Ψ K 2 -f53 + f36 = Ψ Ψ Ψ K 3 -f14 -f64 + f45 + f42 = Ψ Ψ Ψ K 4 -f25 -f45 + f56 + f53 = Ψ Ψ Ψ K 5 -f36 -f56 + f64 + f61 = Ψ Ψ Ψ K 6 ( 53
)
System (53) has a very neat interpretation : the sub-triangulation of figure 4 defines a triangulation of element K associated to its degrees of freedom. For any edge between two degrees of freedom, say [i, j], we look for fluxes fij satisfying [START_REF] Roe | Linear advection schemes on triangular meshes[END_REF], with the constraint fij + fji = 0.

In the P case, one easily shows that the matrix associated to (53) has rank 5, which can be used to obtain a definition of the equivalent finite volume fluxes as done for linear elements. We only sketch the generalization to P k elements which relies on the following main elements

• construct a triangulation of K which vertices are the degrees of freedom of the interpolation ;

• associate to his sub-triangulation a dual tessellation to be used to define local conservation equations ;

• set Ψ Ψ Ψ K i = φ φ φ K i - C i ∂K f (w h ) • n dΓ
• write equations for conservative edge fluxes fij. Assemble a liner system for a subset F of the ordered couples (i, j) associated to the edges of the sub-triangulation, with F containing either (i, j), or (j, i) for any two fixed nodes. We have that 1. the matrix coefficients of the linear system are

θij =    0 (i, j) is no an edge 1 (i, j
) is an edge and (i, j) ∈ F -1

(i, j) is an edge and (i, j) / ∈ F 2. the i-th right hand side of the system is equal to Ψ Ψ Ψ K i 3. the rank of the system matrix is equal to n dof -1

Our analysis can be also generalized to three space dimensions, and to other type of elements, since it only relies on the possibility of constructing a set of connected dual cells, which is possible for any mesh. This analysis shows that whatever the type of element, the approximation of the residual distribution spatial discretisation can be reformulated by means of a finite volume approximation defined by multidimensional fluxes function of n dof states. Hence all continuous finite element schemes admitting a residual distribution re-formulation are locally conservative.

WENO-RD and bridge with DG

One can slightly extend the RD formalism. In what is written above, the main assumption is that the approximation is globally continuous. This assumption can be relaxed. Assume that, as for DG, the trial function space is made of functions that are polynomials on each elements, but we relax the continuity assumption. This problem has been studied in [START_REF] Abgrall | Development of residual distribution schemes for the discontinuous galerkin method: the scalar case with linear elements[END_REF][START_REF] Abgrall | A residual method using discontinuous elements for the computation of possibly non smooth flows[END_REF][START_REF] Hubbard | Discontinuous fluctuation distribution[END_REF]. In each element K, we assume that we have NK degrees of freedom, say {ii j , j = 1, NK }, and assume that we have constructed residuals Φ K i i j (u h ). The conservation relation must be relaxed into

N K j=1 Φ K i i j (w h ) = ∂K f (w h,+ , w h,-, n)d∂K ( 54
)
where f is a consistent numerical flux, w h,+ , w h,-are the states on the two sides of the faces that makes ∂K. In [START_REF] Abgrall | A residual method using discontinuous elements for the computation of possibly non smooth flows[END_REF], it is shown how to reformulate a DG method with P 1 elements in this framework. Though limited to P 1 element, this approach can easily be extended to higher order of approximation, see also [START_REF] Abgrall | A residual method using discontinuous elements for the computation of possibly non smooth flows[END_REF] for a more systematic (than [START_REF] Abgrall | Development of residual distribution schemes for the discontinuous galerkin method: the scalar case with linear elements[END_REF]) technique. In the discontinuous case, a simple variant can be found, see [START_REF] Abgrall | A residual method using discontinuous elements for the computation of possibly non smooth flows[END_REF]. This remark make it possible to use the technique that we describe in section 3.1.7. A completely different approach has been pursued in [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF]. Starting from a finite difference-like grid and using WENO reconstruction, these authors have been successful in developing RD-like approximation for hyperbolic problems.

A general Lax-Wendroff result

One of the key constraint a RD scheme must fullfil is that for any element or face, the sum of the sub-residual must equal the total residuals, these are the conditions [START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using WENO limiters II: unstructured meshes[END_REF] and [START_REF] Chavent | The local projection p 0 p 1 -discontinuous galerkin finite element method for scalar conservation laws[END_REF]. 3 These conditions guaranty a Lax-Wendroff like result, see [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF]. More precisely, if we assume that: Assumption 3.1. The mesh T h is conformal and regular. By regular we mean that all elements are roughly the same size, more precisely that there exist constants C1 and C2 such that for any element

K, C1 ≤ sup K∈T h h 2 |K| ≤ C2
We introduce the spaces:

V k h = {v h ∈ C 0 (R d ) p ; v h|K polynomial of degree k, ∀ K ∈ T h } X h = {v h ; v h|C constant ∈ R p , ∀ C ∈ C h }.
Here, f |K denotes the restriction of f to K. The second assumption is on the residuals.

Assumption 3.2. Let T h be a triangulation satisfying the assumption 3.1. For any C ∈ R + , there exists

C (C, T h ) ∈ R + which depends only on C and T h such that for any w ∈ X h , with ||w|| L ∞ (R 2 ) ≤ C we have ∀K, ∀i, ||Φ K i || ≤ C (C, T h ) h j∈K ||w(j) -w(i)||. ( 55
)
We assume that the residual Φ

K K i
and the numerical solution satisfy the following conditions.

Assumption 3.3. There exists an approximation f h of the flux f such that

(i) ∀w h ∈ X h , Φ K := K div f h (w h )dx = i∈K Φ K i (w h ), (ii) ∀w h ∈ X h , ∀K1, K2 neigbors, f h (w h ) |K 1 . n = f h (w h ) |K 2 . n a.e. on K1 ∩ K2
where n is a normal of K1 ∩ K2.

(iii) For any C > 0, there exists C (C) such that for any

w h ∈ X h with ||w h || L ∞ (R 2 ) ≤ C, one has for K ∈ T h and f h K = f h |K , || div f h K (u h )|| ≤ C h i,j
||w h i -w h j || a.e. on K.

(iv) For any sequence

(w h ) h bounded in L ∞ (R 2 × R + ) p independently of h and convergent in L 2 loc (R 2 × R + ) p to w, we have lim h ||f h (w h ) -f (w)|| L 1 loc (R d ×R + ) p = 0.
We have the following result, Theorem 3.4. Let be w0 ∈ L ∞ (R d ) p and w h the approximation given by

w n+1 i = w n i - ∆ |Ci| K i φ φ φ K i (w n+1 ) + F ∩∂Ω i φ φ φ F i w 0 i = w0(i).
We assume that the scheme satisfies the assumptions 3.2 and 3.3. We also assume there exists a constant C that depends only on C1, C2 and u0 and a function w

∈ (L 2 (R d × R + )) p such that sup h sup x,y,t |w h (x, y, t)| ≤ C lim h ||w -w h || L 2 loc (R d ×R + ) p = 0 Then w is a weak solution of ∂w ∂t + div f (w) = 0 w(x, 0) = w0(x)
The proof is given in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF].

Construction of non-classical high order schemes

Well posed linear schemes

The simplest method is obtained by splitting the element residuals in a symmetric manner, as e.g.

φ φ φ K i = 1 n dof φ φ φ K
This definition verifies trivially all the accuracy criteria, and the conditions for the Lax-Wendroff theorem. Nevertheless, it is flawed by the existence of spurious mode discussed in section 3.1.2 which is a clear symptom of a lack of stability. An example of a stabilized method is the Lax-Friedrich's type distribution

φ φ φ LF i = 1 n dof φ φ φ K + αK (wi -wK ) ( 56
)
with wK the arithmetic average of the solution values in K. In the scalar case, we can easily prove the stability of this method in both the L 2 and L ∞ norms (cf. section 3.1.3). This method does not verify the accuracy conditions of section 3.2.1.

To obtain a stable high order method we can follow the ideas of [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF], and add to an unstable high order method a streamline dissipation term

φ φ φ K i = β K i φ φ φ k + θK K (∇wf • ∇ϕi)TK (∇wf • ∇w h ) dx (57)
with θK a scalar coefficient, and where, for generality, we have replaced 1/n dof by a generic bounded distribution coefficient. Following [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF], we seek for rules to define the term θK TK for scalar advection. We start by recasting the method obtained with (57) as

- Ω h u h a • ϕi dx + K∈K i K (β K i -ϕi)( a • ∇u h ) dx + K∈K i θK K ( a • ∇ϕi)TK ( a • ∇u h ) dx = boundary cond.s
We can associate to this method the bi-linear form

a(v h , u h ) + K bK (v h , u h ) = l bc.s (v h ) where a(v h , u h ) = a Galerkin (v h , u h ) + K aK (v h , u h ) with aK (v h , u h ) = K (v β K -v h )( a • ∇u h ) dx and where bK (v h , u h ) is the streamline dissipation term bK (v h , u h ) = K ( a • ∇v h )θK TK ( a • ∇u h ) dx
For simplicity, and following the ideas of [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF], we now express the increment v β K -v h as a function of ∇v h , and of a (scheme dependent) element length h β K and direction

ξ β K aK (v h , u h ) + bK (v h , u h ) = K (ξ β K • ∇v h ) h β K ( a • ∇u h ) dx + K ( a • ∇v h ) θK TK ( a • ∇u h ) dx
Finally, we want to define the coefficient θK TK such that

aK (v h , u h ) + bK (u h , u h ) = K (ξ β K • ∇u h ) h β K ( a • ∇u h ) dx + K θK TK ( a • ∇u h ) 2 dx ≥ 0 (58)
In particular, to have [START_REF] Struijs | Progress on multidimensional upwind Euler solvers for unstructured grids[END_REF] satisfied in practice, we consider the fully discrete evaluation of the streamline dissipation term

K ( a • ∇ϕi) θK TK ( a • ∇u h ) dx ≈ θK TK |K| x quad ω quad ( a • ∇ϕi(x quad )) ( a • ∇u h (x quad )) (59)
having assumed for simplicity that a constant value of TK is used over each element. We seek now guidelines to choose a quadrature formula.

A necessary condition to have [START_REF] Struijs | Progress on multidimensional upwind Euler solvers for unstructured grids[END_REF] is that the quadratic form

qK (u h ) := |K| x quad ω quad ( a • ∇u h (x quad )) 2
must be positive whenever a • ∇u h = 0. A sufficient condition for this to be true is that

if a • ∇u h (x quad ) = 0 ∀ x quad then a • ∇u h = 0 (60)
In this case, we can find positive bounded constants such that

C1,qqK (u h ) ≤ K ( a • ∇u h ) 2 dx ≤ C2,qqK (u h )
Since V h = span{ϕi} is a finite-dimensional spaces, the discrete quantity

Q(u h ) = K qK (u h ) defines on V h a norm equivalent to u h → Ω ( a • ∇u h ) 2 dx
This allows to prove the following result.

Proposition 3.8 (Quadrature of the streamline dissipation). Independently on the values of the weights ω quad , provided that the number of points used ti evaluate (59) is large enough to guarantee (60) then we can find (θK TK )0 such that the scheme obtained with (57) is well posed whenever θK TK ≥ (θK TK )0.

Proof. See [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF].

In light of this analysis, the set of points used to evaluate [START_REF] Deconinck | A multidimensional generalization of Roe's difference splitter for the Euler equations[END_REF] does not need necessarily be a set of quadrature points, as the relevant condition is not to evaluate the streamline dissipation term exactly, but to ensure [START_REF] Nishikawa | A third-order fluctuation splitting scheme that preserves potential flow[END_REF]. In this light, term (59) can be seen as a sort of filtering term allowing to remove spurious modes and guarantee the well-posedness of the method. In particular, even for constant scalar advection, the number of points sufficient to have [START_REF] Nishikawa | A third-order fluctuation splitting scheme that preserves potential flow[END_REF] is smaller than that of most quadrature/cubature formulas providing an exact evaluation of the streamline dissipation term, and in any case simpler point values can be used, as e.g. the x dof (cf. [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF]).

Another path to avoid the flaw associated to proposition 3.5, is to "bring the distribution coefficient β under the integral". The classical definition associated to SUPG [START_REF] Gottlieb | Strong Stability-Preserving High-Order Time Discretisation Methods[END_REF] is one example of such a method. However, we can also provide similar generalizations of the so-called multidimensional upwind methods, which have constituted one of the elements of originality of the residual distribution approach (cf. [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF] and references therein). In particular, we consider the method defined in the scalar case by

φ K i = K ( a • ∇ϕi) + γK ( a • ∇u h ) dx (61)
with γK > 0. If, without loss of generality, we consider the linear advection problem admitting a solution u > 0 4 , we may assume that we seek a discrete solution

u h ∈ V + h = {u h ∈ span{ϕi}|u h > 0}, we can in this case write i∈Ω h ui K∈K i φ K i = i∈Ω h K ( a • ∇u h ) + γK a • ∇u h dx = i∈Ω h K ( a • ∇u h ) + γK ( a • ∇u h ) + dx ≥ 0
4 which we may always do in the linear case by rescaling the boundary data g by

M = | min x∈Ω g| + M , M > 0 as a • ∇u h = ( a • ∇u h ) + + ( a • ∇u h ) -
, and ( a • ∇u h ) + ( a • ∇u h ) -= 0. This shows that condition [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF] is met, giving an indication of the well-posedness of the method, confirmed by all numerical evidence. Note that for the method to satisfy (39), we need to set

γK = j∈K ( a • ∇ϕj) + -1
which can be show to reduce in the P 1 case leads to

φ K i = β K i φ K , β K i = ( a • ni) + / j∈K ( a • nj) +
which is nothing else that the well known multidimensional upwind LDA scheme introduced by Roe, Deconinck and co-workers in the 90's [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. The generalization ( 61) is obtained by formally replacing the so-called upwind parameters ki = a • ni/2 by the i-th streamline component of the solution gradient ki = a•∇ϕi, and by performing the distribution of the local residual, instead of distributing the integrated cell residual φ K . The extension to nonlinear problems is obtained by replacing in [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF] the residual a•∇u h , with ∇ • f h , with f h a higher order polynomial, of at degree at least k + 1 (one higher than the solution), built starting from the values of u h . Refer to [START_REF] D'angelo | Adjoint-based error estimation for adaptive petrovgalerkin finite element methods[END_REF][START_REF] Vymazal | Highorder residual distribution and error estimation for steady and unsteady compressible flow[END_REF] for more details.

Non-oscillatory methods

The analysis of section 3.1.3 constitutes the basic artillery used in the past years to construct methods allowing a non-oscillatory approximation of discontinuous solutions. The key element of these constructions is some low (first) order linear scheme, satisfying [START_REF] Krivodonova | Shock detection and limiting with discontinuous galerkin methods for hyperbolic conservation laws[END_REF]. A typical example of such a scheme is given by [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF]. For this scheme, and in the scalar case, we can readily prove that (36) holds in d space dimensions, as soon as

αK > h d-1 K ∇ bbw f L ∞ (K)
, with hK a characteristic length scale of the element. A neat way of producing a formally high order method starting from [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF], is to fabricate uniformly bounded distribution coefficients by applying a nonlinear mapping to the quantities φ LF i /φ K . An example of such a mapping is the well known "PSI limiter" [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]

β LLF i = (φ LF i /φ K ) + j∈K (φ LF j /φ K ) + = (φ LF i φ K ) + j∈K (φ LF j φ K ) + (62)
For the corresponding scheme, one easily shows that

φ LLF i = β LLF i φ K = γiφ LF i , γi = (φ LF i /φ K ) + j∈K (φ LF j /φ K ) + φ K /φ LF i ∈ [0, 1] (63)
Thus, this limited Lax-Friedrich's distribution is by construction stable in the L ∞ norm. Unfortunately, we also know from proposition 3.5 that this scheme will not in general converge, and anyways may not converge to the right solution. A cure to this problem has been suggested in [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF], and consists in adding the filtering term (59) in smooth regions. The resulting method reads

φ LLFs i = β LLF i φ K + θ(u h )|K| x quad ω quad ( a • ∇ϕi(x quad )) TK ( a • ∇u h (x quad )) (64)
where θ(u h ) is defined such that the conditions of proposition 3.8 are met in smooth regions, while θ < O(hK ) in vicinity of discontinuities. Practical definitions of this term can be found in [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF]. The extension of this construction to systems is performed by computing the limiter (62) either equation by equation, or by a prior projection of the residuals on characteristic directions, and by replacing the advection vector in (64) by the flux Jacobian matrices. A common definition of the scaling matrix TK is

TK = |K| v (∇wf (wv) • ∇ϕv(xv)) + -1
with v the vertices of element K.

An alternative construction consists in adding to a linear high order and stable scheme a local amount of shock capturing dissipation. This approach dates back a long way [START_REF] Hughes | A new finite element formulation for CFD IV: a discontinuitycapturing operator for multidimensional advective-diffusive systems[END_REF]. In the framework of residual distribution schemes it has been reformulated by means of a technique reminiscent of flux limiting in the finite volume context : the nonlinear blending of a linear high order method with a linear low (first) order positive coefficient one. For example, blending [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF] with a high order stabilized method would lead to

φ φ φ K i = 1 n dof φ φ φ K + δ(w h )αK (wi -wK ) + (Id -δ(w h ))θK |K| x quad ω quad ( a • ∇ϕi(x quad )) TK ( a • ∇u h (x quad )) (65
) where different forms of the stabilization as selected depending on whether w h is smooth, in which case δ(w h ) ≤ O(hK ), or discontinuous, in which case Id -δ(w h ) ≤ O(hK ). More involved constructions considering replacing [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF] in the blending by [START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF] have also been proposed e.g. in [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF].

Handling source terms

The extension of the above framework to the approximation of solutions of div f (w) + s(w, x) = 0 [START_REF] Abgrall | An example of high order residual distribution scheme using non-lagrange elements[END_REF] is based on the inclusion of the source in the re-definition of the local element residual, leading to the requirement

i∈K φ φ φ K i = ∂K f h (w h ) • n d∂K + K s h (w h , x) dx (67)
All of the methods described earlier can be extended to this more general setting.

Interesting results can be shown when s depends on some given data, say a given field f (x) :

s(w, x) = s(w, f (x))
This is the case in some environnemental applications (e.g. shallow water equations), or when considering the solution of the differential problem on a manifold (see e.g. [START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems on curved manifolds[END_REF] and references therein). Such problems often embed some particular solutions which are characterized by the existence of a set of invariants v = v(w, f ) constant throughout the spatial domain. Assuming a sufficient smoothness of f , of the solution, and of the mapping (v, f ) → w(v, f ), in this case we can write that div f (w) = (∇wf ∇vw) ∇v + (∇wf ∇ f w)

Λ(v,f ) ∇f Solutions characterized by the invariance relation v = v0 = c t ∀x, satisfy (cf. ( 66
)) Λ(v0, f )∇f + s(v0, f ) = 0 (68)
An interesting result concerning this class of solutions for schemes defined by

φ φ φ K i = K ω K i (∇ • f h (v h , f ) + s h (v h , f )) (69)
thus based on a direct approximation of the invariant states, instead of the conserved variables w. The following is shown in [START_REF] Ricchiuto | On the c-property and generalized c-property of residual distribution for the shallow water equations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Ricchiuto | Contributions to the development of residual discretizations for hyperbolic conservation laws with application to shallow water flows[END_REF].

Proposition 3.9 (Steady invariants and superconsistency). Under standard regularity assumptions on the mesh, provided that (69) is true for some test function ω K i which is uniformly bounded w.r.t. h, v h , element residuals, and w.r.t. to the data of the problem, then for exact integration the scheme defined by (69) preserves exactly the equilibrium [START_REF] Ciarlet | General lagrange and hermite interpolation in R n with applications to finite element methods[END_REF]. For approximate integration, assuming that a flux quadrature exact for approximate polynomial fluxes of degree p f is used, and a source quadrature exact for approximate polynomial sources of degree pv, and assuming that f ∈ H p+1 with ∇f ∈ H p , and p > min(p f , pv), then the scheme defined by (69) is super-cosnsistent w.r.t. solutions characterized by [START_REF] Ciarlet | General lagrange and hermite interpolation in R n with applications to finite element methods[END_REF], and in particular, its consistency is of order r = min(p f + 2, pv + 3).

This result shows the potential of residual framework considered here in guaranteeing the balance of the flux divergence with complex source terms. Similar and stronger results, especially on simple particular solutions encountered in shallow water flows, have already been recalled in the previous chapter [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. Some applications of this property will be shown in the results section.

Handling viscous terms

When considering the advection diffusion equation (with ∇ • a = 0):

a • ∇u -∇ • (K(u) • ∇u) = 0 ( 70
)
with K(u) a positive semi-definite diffusion matrix coefficient, the first idea is to look at it as a standard conservation relation with an enhanced flux, now

f (u, ∇u) = au -K(u) • ∇u,
to which, one could apply the same construction as before. However, there is a fundamental difference: if the approximation of the solution is sought to be piecewise polynomial and globally continuous, its gradient will still be piecewise polynomial, but will not be globally continuous anymore. One of the fundamental requirement of the previous developments is that the flux on the boundary of the element is single valued. This cannot be anymore the case here unless something is done.

There are two ways of solving this issue. Both are similar what is done in LDG methods. The first step is in each case to rewrite the partial differential equation into a, possibly hyperbolic, first order system of PDEs. For the two-dimensional advection-diffusion equation, setting K = ν Id, we consider the hyperbolic first order system

∂u ∂t + a • ∇u = ν (∂xp + ∂yq) ∂p ∂t = 1 Tr (∂xu -p) ∂q ∂t = 1 Tr (∂yu -q) (71)
where p and q are the gradient variables and Tr is a relaxation time. At the steady state the system (71) is equivalent to the original equation ( 70), independently of the parameter Tr, and p,q become equivalent to the derivatives of the unknown. The idea of reformulating a parabolic problem with second order derivatives as a hyperbolic system such as e.g. [START_REF] Abgrall | Residual distribution schemes for conservation laws via adaptive quadrature[END_REF], is not new, as it dates back to the work of Vertnotte and Cattaneo [START_REF] Vernotte | Les paradoxes de la theorie continue de l?equation de la chaleur[END_REF][START_REF] Cattaneo | A form of heat-conduction equations which eliminates the paradox of instantaneous propagation[END_REF] to study the heat equation. This idea has been efficiently exploited by H. Nishikawa and A. Mazaheri to construct schemes for the steady and time dependent diffusion, advectiondiffusion, advetcion-diffuction and Navier-Stokes equations, see e.g. [START_REF] Nishikawa | A first-order system approach for diffusion equation. i: Second-order residualdistribution schemes[END_REF][START_REF] Nishikawa | A first-order system approach for diffusion equation. ii: Unification of advection and diffusion[END_REF][START_REF] Nishikawa | Beyond interface gradient: A general principle for constructing diffusion scheme[END_REF][START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme I: Basic principle[END_REF] and [START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] for recent formulations of residual distribution and Discontinuous Galerkin based on this approach. There are two ways to approach system [START_REF] Abgrall | Residual distribution schemes for conservation laws via adaptive quadrature[END_REF]. The first is to make use of its hyperbolicity, and reuse all the artillery already available. In this case, the overhead of having to introduce the gradient variables, can be compensated, by a careful design of the scheme which may guarantee the same accuracy for both the solution and its derivatives. This may have an impact on the computation of e.g. forces and heat fluxes and allow the use of coarser meshes to provide accurate values of these quantities. This is a path being followed in [START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF], but the demonstration of its feasibility for practical application is still a work in progress.

Another way to exploit system (71) has been suggested in [START_REF] Nishikawa | Beyond interface gradient: A general principle for constructing diffusion scheme[END_REF][START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme I: Basic principle[END_REF], and developed from scalar advection diffusion up to laminar Navier-Stokes and RANS equations in [START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF][START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF][START_REF] Santis | High-order linear and non-linear residual distribution schemes for turbulent compressible flows[END_REF]. In this alternative approach only the first discrete equation or u h is kept. This is of course a discretization of ( 70), which however depends on values of p and q. These values are now replaced by an appropriate high order reconstruction of the solution derivatives starting from u h . Simple solutions are possible, as e.g. the use of simple arithmetic averages for the viscous fluxes on element boundaries, see for example [START_REF] Abgrall | Numerical approximation of parabolic problems by means of residual distribution schemes[END_REF]. However, these simple choices lead to suboptimal accuracy, mostly because one order of accuracy is lost in the evaluation of the gradient. In [START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF], a systematic study of possible recovery methods (arithmetic average, least square, ect) has been conducted, and the best solution is to take into advantage, via a local least square minimization algorithm, of the existence of super convergence points in element. At these points, as put forward by Zienkiewiscz and Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF], the gradient is approximated at full order. In the following, this reconstruction will be nicknamed as SPR-ZZ.

Example of a non-classical scheme. We want to show how to use these ideas to generalise the schemes of section §3.1.7 for the solution of [START_REF] Barth | An energy look at the N scheme[END_REF]. The schemes obtained are those used in the numerical results we will discuss later.

So we start from scheme (57), assuming for simplicity β K i = 1/n dof . If we assume to be in the purely diffusive case, apply this scheme to system (71), and we only look at the local distributed residuals only for the first equation we have

φ fos i = 1 n dof φ K,ν + K τν ∇ϕi • (∇u h -(p h , q h ))
having assumed for the stabilization matrix TK = δν Id, and where

φ K,ν = - ∂K ν(p h , q h ) • n
Note also that the effect of the relaxation time Tr has been embedded in the δν coefficient.

The trick is now to replace the nodal values of the gradients p and q by accurately reconstructed ones, which we obtain with the SPR-ZZ procedure recalled above. The important part is the definition of the total residual. From the Lax-Wendroff theorem recalled in section §3.1.6 however we know that the numerical approximations of bot these fluxes must be edge-continuous. The simplest way to achieve that is to use for the viscous flux the finite element approximation based on the reconstructed nodal gradients. We denote this quantity by ∇u h . So, for pure diffusion, the scheme is finally defined by

φ K,ν i = 1 n dof φ K,ν + K δν ∇ϕi • ∇u h -∇u h where φ K,ν = - ∂K ν ∇u h • n
and with δν having the dimensions of a diffusion coefficient.

For advection diffusion, we can apply the same procedure. Starting with the total residual

φ K = ∂K a u h -ν ∇u h • n d∂Ω,
we can deduce from the first order system formulation two types of regularization terms leading to local nodal residual5

φ K i = 1 n dof φ K + K τa a • ∇ϕi ( a • ∇u h -ν∇ • ∇u h ) + K δν ∇ϕi • ∇u h -∇u h
The optimal choice of the scaling parameters τa and δν has been shown to require some dependence on the elemental Re number Re = λ h ν , see for example [START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF][START_REF] Nishikawa | A first-order system approach for diffusion equation. ii: Unification of advection and diffusion[END_REF][START_REF] Ricchiuto | On uniformly high-order accurate residual distribution schemes for advection-diffusion[END_REF] and references therein. This is taken into account by setting

φ K i = 1 n dof φ K + ξ(Re) K a • ∇ϕi τ a • ∇u h -∇ • (ν∇u h ) dΩ + 1 -ξ(Re) K νδ 2 ∇u h -∇u h • ∇ϕi dΩ, (72)
where the function ξ(Re) is such that ξ(Re) → 0 in the diffusion limit (Re → 0) and ξ(Re) → 1 in the advection limit (Re → ∞).

To account for non-smooth solutions, one can use the same technique discussed in section §3.1.7, replace the centered contribution by a nonlinear limited residual, and pre-multiply the stabilisation terms by some smoothness sensor, so that the scheme can be generally written in the finale general form

φ K i =β K i φ K +θ(u h )ξ(Re) K a • ∇ϕi τ a • ∇u h -∇ • (ν∇u h ) dΩ +θ(u h ) 1 -ξ(Re) K νδ 2 ∇u h -∇u h • ∇ϕi dΩ, (73)
where β K i is computed following the limiting procedure discussed in section 3.1.7 in the non-smooth case.

The numerical scheme obtained for the advection-diffusion scalar equation is then extended to the compressible Navier-Stokes equations. The governing equations read

∂w ∂t + ∇ • f a (w) -∇ • f v (w, ∇w) = 0,
where w and f a (w) are the vector of the conservative variables and the advective flux function, respectively, as defined for the Euler equations, while are the components of the stress tensor, with µ the dynamic viscosity of the fluid, and qx, qy are the components of the heat flux q which is defined as

f v (w, ∇w) = ( f v x , f v y ) T is the viscous flux function f v x (w, ∇w) =       0 τxx τxy uτxx + vτxy -qx       , f v y (w, ∇w) =       0
q = k∇T,
where T is the temperature and k is the thermal conductivity coefficient. It is well know that the viscous flux function f v is homogeneous with respect to the gradient of the conservative variable ∇w

f v (w, ∇w) = K(w)∇w,
with the homogeneity tensor K(w) = ∂f d ∂w . The discretization of the Navier-Stokes equations is straightforward. The total residual on a generic element K is given by

Φ K = ∂K f a (w) -K(w) ∇w • n,
with ∇w the reconstructed gradient of the conservative variables and the boundary integral computed by the means of a quadrature rule. The total residual is first distributed to all the DOF of the element using the low order Rusanov scheme, subsequently the limitation procedure is applied to obtain an high order residual as described in the section 4.2.1. In the last step the filtering term is added together with the dumping term acting for the viscous part. The complete scheme reads

Φ K i = ΦK i + ξ(Re) K A • ∇ϕi τ A • ∇w h -∇ • K∇w h dΩ + 1 -ξ(Re) E 1 2 K ∇w h -∇w h • ∇ϕi dΩ. ( 74
)
with ΦK i denoting the (unfiltered) centered or nonlinear distribution.

Time dependent problems

In this chapter, we consider the approximation of time dependent solutions to a system of conservation laws reading

∂tw + ∇ • f (w) = 0 on Ω × [0 , T fin ] ⊂ R d × R + (75)
As shown in [START_REF] Struijs | A Multi-Dimensional Upwind Discretization Method for the Euler Equations on Unstructured Grids[END_REF], and then in [START_REF] Caraeni | Development of a Multidimensional Upwind Residual Distribution Solver for Large Eddy Simulation of Industrial Turbulent Flows[END_REF][START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF], and in [START_REF] Maerz | Improving time accuracy of residual distribution schemes[END_REF][START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF][START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF][START_REF] Ricchiuto | Residual distribution for general time dependent conservation laws[END_REF] (cf. also the chapter [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]), to obtain high order schemes for this case, one must carefully design a coupling between the stencil used to approximate the integral of the time derivative and the flux divergence. Some approaches to obtain this coupling are recalled, and a more general prototype is analyzed. The links with other methods are briefly recalled. The first part of the section is devoted to fully implicit methods. We then discuss a path allowing the contruction of explicit apporaches which do-not require the inversion of a mass matrix, or for which this matrix reduces to the symmetric positive-definite Galerkin one. Note that, beside classical stabilized finite element schemes (SUPG, GLS, etc), here the status of residual distribution type methods is less advanced. Here we discuss some of the most interesting ideas toward generalizing the methods presented for steady state. Some research directions to push the limits of the existing construtions will be discussed later.

Implicit prototype for time dependent solutions

We introduce the time discretized version of ( 75) by means of an r + 1th order time integration scheme

Γ n+1 (w) = p i=0 αi δw n+1-i ∆t + q j=0 θj∇ • f n+1-j (76)
where ∆t = minn(t n+1 -t n ), with ∆t n+1 = t n+1 -t n , with δw n+1 = w n+1 -w n , f n+1-j = f n+1-j (w n+1-j ), and with the αi and θj coefficients given by a time integration scheme of choice. This may be a generic stage of a multi-stage method, or a multi-step scheme. Space time schemes can be embedded in the analysis that follows by appropriate definitions of the αis, and of the δw n+1-i to embed eventually jumps in the time direction, when using discontinuous in time space-time elements. An important assumption, is that the time stepping verifies the conservation identity

N n=0 p i=0 αiδw n+1-i = w N -w 0 = w(T fin ) -w0 (77)
We set on every K ∈ Ω h

Φ K = K Γ n+1 (w h ) = K p i=0 αi δw n+1-i h ∆t + q j=0 θj∇ • f n+1-j h (78)
with w h and f h continuous finite element polynomial approximations of degree k and (at least) k respectively. Similarly, on each boundary face f we set

φ φ φ f = f q j=0 θj(ĝ -f h ) n+1-j • n (79)
with ĝ a numerical flux consistent with the BCs. Similarly to the previous sections, we consider the scheme that computes w h as the solution of

K∈K i Φ K i + f ∈F i φ φ φ f i = 0 ( 80
)
where ∀ K and ∀ f j∈K

Φ K j = Φ K and j∈f φ φ φ f j = φ φ φ f (81)

Consistency analysis

To begin with, we generalize the consistency conditions. To simplify the notation we consider the scalar case, and we neglect the boundary conditions, which can be easily embedded in the spatial operator as shown.We will assume some classical regularity properties for the mesh and the time stepping strategy, namely

C0 ≤ sup K∈Ω h h 2 |K| ≤ C1 , C 0 ≤ ∆t h ≤ C 1 (82)
Let now w ∈ C l+1 be an exact classical solution of (75), with l ≥ max(r, k), and such that

p i=0 αi δw n+1-i ∆t + q j=0 θj∇ • f n+1-j = ∂tw + ∇ • f + O(∆t r+1 ) (83)
We denote by w m h , the kth degree continuos finite element projection/interpolation of w m . Consider now ψ ∈ C 1 0 (Ω × [0, T fin ]), a smooth test function with ψ| ∂Ω = 0. Let ψ h be its k degree polynomial finite element projection/interpolation, with ψ n i the corresponding values at the choses degrees of freedom. It is also assumed that [START_REF] Ciarlet | General lagrange and hermite interpolation in R n with applications to finite element methods[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF] there exist constants C 0 , C 1 , C2 such that

∂tψ h L ∞ (Ω h ) ≤ C 0 , ψ h (t + ∆t) -ψ h (t) L ∞ (Ω h ) ≤ C 0 ∆t ψ h L ∞ (Ω h ) ≤ C 1 , |ψi -ψj| ≤ ∇ψ h L ∞ (Ω h ) h ≤ C2h (84)
We define the following truncation error for scheme ( 80)

(w h , ψ) := N n=0 i∈Ω h ∆t n+1 ψ n+1 i K∈K i Φ K i (w h ) = N n=0 K∈Ω h i∈K t n+1 t n ψ n+1 i Φ K i (w h ) (85)
We introduce the Galerkin splitting in space

Φ G i = K ϕiΓ n+1
and note that j∈K

(Φ K j -Φ G j ) = 0
This allows to recast the error as

(w h , ψ) = N n=0 t n+1 t n    Ω h ψ n+1 h Γ n+1 (w h ) + 1 CK K∈Ω h i,j∈K (ψi -ψj)(Φ K i -Φ G i )    (86)
Multiplying ( 83) by ψ h and integrating over space and time we can get

N n=0 t n+1 t n Ω h ψ n+1 h Γ n+1 (w) = N n=0 ∆tO(∆t r+1 ) = O(∆t r+1 )
So the error can be estimated as follows

(w h , ψ) = I + II + III + O(∆t r+1 ) I = N n=0 t n+1 t n Ω h ψ n+1 h p i=0 αi δ(w h -w) n+1-i ∆t II = N n=0 t n+1 t n q j=0 Ω h ψ n+1 h ∇ • (f h -f ) n+1-j III = N n=0 t n+1 t n 1 CK K∈Ω h i,j∈K (ψi -ψj)(Φ K i -Φ G i )
Estimating each of the terms we get to the conditions of the cell and boundary splittings allowing to preserve the O(∆t r+1 ) appearing on the right hand side. This is readily done by using the hypotheses on the regularity of u and standard interpolation results [START_REF] Ciarlet | General lagrange and hermite interpolation in R n with applications to finite element methods[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF]. In particular, for term I we can use hypothesis [START_REF] Johnson | On the convergence of a finite element method for a nonlinear hyperbolic conservation law[END_REF] to write

I = N n=0 t n+1 t n Ω h p i=0 αi δ(ψ h w h -ψ h w) n+1-i ∆t + N n=0 t n+1 t n Ω h p i=0 αi(ψ n+1 h -ψ n-i+1/2 h ) δ(w h -u) n+1-i ∆t - N n=0 t n+1 t n Ω h (w h -u) n-i+1/2 p i=0 αi δ(ψ h ) n+1-i ∆t = Ω h (ψ h (w h -u))(T fin ) - Ω h (ψ h (w h -u))0+ N n=0 t n+1 t n Ω h p i=0 αi(ψ n+1 h -ψ n-i+1/2 h ) δ(w h -u) n+1-i ∆t - N n=0 t n+1 t n Ω h (w h -u) n-i+1/2 p i=0 αi δ(ψ h ) n+1-i ∆t
Using [START_REF] Huynh | High-order methods for computational fluid dynamics: A brief review of compact differential formulations on unstructured grids[END_REF], and the regularity of u, we can now bound this term as

|I| = O(h k+1 ) + C T fin ∆t ∆t O(h k+1 )C 0 sup i=1,p |αi| = O(h k+1 )
The analysis of the remaining terms is practically identical to the one of section §3.1.1, and omitted for brevity (the interested reader can refer to [START_REF] Ricchiuto | Contributions to the development of residual discretizations for hyperbolic conservation laws with application to shallow water flows[END_REF] for details). The final result is the following.

Proposition 3.10 (Accuracy of RD, unsteady case). Under assuption (82) on the time stepping, given a k + 1th order continuous polynomial approximation of the unknown ad of the fluxes, and a r + 1th order accurate time integration scheme, scheme (80) verifies the truncation error estimate

| (w h , ψ)| ≤ O(h p+1 ) , p = min(k, r) provided that sup K∈Ω h sup i∈K |Φ K i (w h )| =O(h p+d ) ( 87
)
whenever w h is the finite element projection/interpolation of a smooth exact solution. In this case we say that the scheme is p + 1th order accurate.

Moreover we have the following estimate.

Lemma 3.11 (Consistency estimate, time dependent case). Under the hypotheses of proposition 3.10 the following consistency estimates hold.

Γ n+1 (w h ) = O(h k ) + O(∆t r+1 ) , Φ K (w h ) = O(h p+d ) (88)
Proof. The proof is easily obtained by considering that due to (83)

Γ n+1 (w h ) = O(∆t r+1 ) + Γ n+1 (w h ) -Γ n+1 (w)
By its definition, and under the hypotheses made, one easily checks that Γ

n+1 (w h ) -Γ n+1 (w) = O(h k ).
The estimate on Φ K (w h ) is trivially obtained upon integration of Γ n+1 .

As a consequence we have the following corollary.

Corollary 3.12 (High order residual schemes). Under the hypotheses of proposition 3.10, a sufficient condition for a scheme of the form (80) to be p + 1th order accurate if there exist a test function ωi uniformly bounded w.r.t. h, w h , Γ n+1 (w h ), and w.r.t the data of the problem, such that

Φ K i (w h ) = K ωiΓ n+1 (w h ) ( 89
)

Examples of implicit high order schemes

Typical examples of high order methods are obtained with the natural extension to the time dependent case case of SUPG-type methods (see e.g. [START_REF] Hughes | Development of time-accurate finite element techniques for first order hyperbolic systems with emphasys on the compressible euler equations[END_REF][START_REF] Shakib | A new finite element formulation for computational fluid dynamics: Ix. fourier analysis of space-time galerkin/least-squares algorithms[END_REF][START_REF] Hughes | Multiscale and Stabilized Methods[END_REF][START_REF] Chalot | Higher-order stabilized finite elements in an industrial navier-stokes code[END_REF] and references therein). Some notable examples of less classical high order schemes exploit ( 89) with ωi = β K i constant per element. The fist of such accuracy preserving schemes are found in the work of D. Caraeni [START_REF] Caraeni | Development of a Multidimensional Upwind Residual Distribution Solver for Large Eddy Simulation of Industrial Turbulent Flows[END_REF][START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF], up to third order of accuracy for the Navier-Stokes equations, and more recently in [START_REF] Rossiello | Third-order-accurate fluctuation splitting schemes for unsteady hyperbolic problems[END_REF] where the third order scheme of Caraeni has been blended with a monotone one via a FCT procedure to provide oscillation free high order solutions of the compressible Euler equations.

Other non-classical constructions have tried to exploit the similarities between stabilized finite elements, and RD methods with constant distribution coefficients. The objective of these works is to find clever definitions of mass matrices/test functions guaranteeing the satisfaction of t [START_REF] Hubbard | Discontinuous fluctuation distribution[END_REF]. This was the initial idea behind the work of Maerz and Ferrante at the von Karman Institute [START_REF] Maerz | Improving time accuracy of residual distribution schemes[END_REF][START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF] later pursued first in [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF][START_REF] Mezine | Conception de Schémas Distributifs pour l'aérodynamique stationnaire et instationnaire[END_REF], and [START_REF] Mezine | Monotone and stable residual distribution schemes on prismatic space-time elements for unsteady conservation laws[END_REF][START_REF] Ricchiuto | Conservative residual distribution schemes for general unsteady systems of conservation laws[END_REF][START_REF] Ricchiuto | Residual distribution for general time dependent conservation laws[END_REF][START_REF] Ricchiuto | Stable and convergent residual distribution for time-dependent conservation laws[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] (see also [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]). This has provided interesting results, but so far only for second order methods.

Finally, examples of space time RD schemes up to third order are discussed in [START_REF] Ricchiuto | Construction of very high order residual distribution schemes for unsteady advection: preliminary results[END_REF][START_REF] Koloszár | Residual distribution method for aeroacoustics[END_REF] with monotonicity preserving extensions discussed in [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF][START_REF] Abgrall | Towards very high-order accurate schemes for unsteady convection problems on unstructured meshes[END_REF][START_REF] Ricchiuto | Residual distribution for general time dependent conservation laws[END_REF][START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF].

All these works, use almost exactly the same techniques developed for steady problems, either treating the time derivative as a source term, or as an additional space direction. The potential of these methods is that they may allow preservation of monotonicity unconditionally w.r.t the time step size, which is very interesting when considering local mesh refinement (see e.g. [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF][START_REF] Sarmany | Unconditionally stable space-time discontinuous residual distribution for shallow water flows[END_REF]), or stiff problems (viscous terms, chemical reactions etc.). The drawback of this formulation, is that the nonlinear stabilization involved depends on the unknown solution at the new time level, thus ruling out a-priori simpler genuinely explicit time marching methods, often preferred in the hyperbolic case. Some exceptions to this rule exist, as for example Taylor-Galerkin, and Lax-Wendroff type methods which can also be recast in a residual distribution formalism (see e.g. [START_REF] Hubbard | Compact high resolution algorithms for time dependent advection problems on unstructured grids[END_REF][START_REF] Ricchiuto | Time accurate solution of hyperbolic partial differential equations using fct and residual distribution[END_REF][START_REF] Rossiello | Second-order-accurate explicit fluctuation splitting schemes for unsteady problems[END_REF], and [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]).

A technique to side step this issue, and construct some non-classical genuinely explicit monotone and high order residual methods is discussed in the next section.

Genuinely explicit time advancement for residual methods

The main idea here is to start from a prototype high order scheme, which we will write in general as (boundary conditions are neglected for simplicity)

Ω h ωi(w h ) (∂tw h + ∇ • f h (w h )) + K∈Ω h ∂K γ ∂K (w h )[∇w h ] • [∇ϕi] = 0 with [•]
a jump of a quantity, as in (16c). The weight ωi in the first term is better expressed as a composition of local restrictions ωi = K ω K i , and depends on the specific method chosen. For SUPGtype schemes we can write ω K i = ϕ K i + γ K i (w h ) with the first term only depending on the mesh. For other methods such as RD schemes, similar decompositions may be invoked, however these are not unique [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF]. Other definitions can be obtained by considering Variational Multi-Scale stabilization techniques, or bubble functions (see [START_REF] Hughes | Multiscale and Stabilized Methods[END_REF] for a review). The last term in the method is one of the possible forms of an edge stabilization [START_REF] Burman | Explicit runge-kutta schemes and finite elements with symmetric stabilization for first-order linear pde systems[END_REF][START_REF] Burman | Interior penalty continuous and discontinuous finite element approximations of hyperbolic equations[END_REF]. Due to the presence of the ∂tw h term in the residual r(w h ), and to the continuity of the approximation, the first term will led to a global mass matrix in the resulting system of ODEs. This matrix in general depends on the discrete solution w h , and, in the case of RD schemes, is not uniquely defined, nor guaranteed to be invertible [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF].

The first idea to simplify things, comes originally from [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF], and requires the introduction of some discrete approximation of the ODE system. As done before, we consider a semi-discretization in time, and the semi-discrete residual which we write here as

Γ n+1 = Γ n+1 (w n+1 h ; {w (s) h }) = α-1w n+1 h + S s=0 αsw (s) h + ∆t S s=0 θs∇ • f h (w (s) h ) ( 90
)
with the w (s) values being either those computed from previous time steps (multi-step scheme) or from some previous predictor stages (multi-stage). Note that the two summations on the right hand side are independent on the unknown w n+1 . As before, for a rth order accurate method in time, the local truncation error relation will be of the type Γ n+1 = O(∆t r+1 ) = O(h r+1 ), if (82) hold, as is the case always for explicit time schemes. If we proceeded as in the last section, we would plug Γ n+1 in the spatial discretization, and the term α-1w n+1 h would lead to the inversion of a (nonlinear) mass matrix. However, in [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF] it is proved that given a kthe order accurate approximation in space, and a rth order accurate approximation in time, provided that the ratio ∆t/h is uniformly bounded, the space-time discretization

Ω h ϕi(Γ n+1 (w n+1 h ; {w (s) h }) -Γn+1 ({w (s) h })) = - Ω h ωi(w h ) Γn+1 ({w (s) h }) -∆t S s=0 βs K∈Ω h ∂K γ ∂K (w (s) h )[∇w h ] • [∇ϕi]
verifies a truncation error/consistency estimate of the type = O(h p ), with p = min(k +1, r +1), provided that for a smooth exact solution, the modified semi-discrete residual Γn+1 verifies the consistency estimate

Γn+1 = O(h p-1 )
The first practical use of this reduced consistency requirement for Γn+1 is to modify a given time discretization to obtain residual expressions one order lower. For example, for the classical third order RK3 method one has [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF] first step Γ

(1)

RK3 =w (1) -

w n + ∇ • f (w n ) Γ(1) RK3 = ∇ • f (w n ) second step        Γ (2) RK3 = w (2) -w n + ∆t 4 (∇ • f (w n ) + ∇ • f (w (1) )) Γ(2) RK3 = w (1) -w n 2 + ∆t 4 (∇ • f (w n ) + ∇ • f (w (1) )) final step      Γ n+1 RK3 = w n+1 -w n + ∆t 6 (∇ • f (w n ) + 4∇ • f (w (2) ) + ∇ • f (w (1) )) Γn+1 RK3 = 2(w (2) -w n )+ ∆t 6 (∇ • f (u n ) + 4∇ • f (w (2) ) + ∇ • f (w (1) ))
For the extrapolated backward differencing method (eBDf3) one finds [START_REF] Klosa | Extrapolated BDF residual distribution schemes for the shallow water equations[END_REF][START_REF] Klosa | Multi-step and multi-stage explicit time stepping for residual distribution[END_REF]]

Γ n+1 eBDf 3 = 11 6 w n+1 -3w n + 3 2 w n-1 - 1 3 w n-2 +∆t(3∇ • f (w n ) -3∇ • f (w n-1 ) + ∇ • f (w n-2 )) Γn+1 eBDf 3 = 5 2 w n -4w n-1 + 3 2 w n-2 +∆t(3∇ • f (u n ) -3∇ • f (w n-1 ) + ∇ • f (w n-2 ))
Equation [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF], can be also seen as a defect correction method, in which a lower order residual is used as a means of approximating solutions of a high order one.

Note however, that relation [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF] still requires the inversion of the Galerkin mass matrix which, even though symmetric positive-definite, is not an inverse monotone matrix. This may destroy all the efforts made in the construction of a shock capturing mechanism in the method. The solution is to constrain the choice of finite element spaces to those allowing to lump this matrix. Several choices exist, either based on standard Lagrange elements on a cubature grids with strictly positive cubature weights [START_REF] Cohen | High order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Giraldo | A diagoanl mass-matrix triangular spectral element method based on cubature points[END_REF][START_REF] Xu | Gauss-Lobatto integration on the triangle[END_REF][START_REF] Mulder | New triangular mass-lumped finite elements of degree six for wave propagation[END_REF], or using non-Lagrange elements having a similar property, as the Bezier basis proposed in [START_REF] Abgrall | An example of high order residual distribution scheme using non-lagrange elements[END_REF] (see also [START_REF] Rogers | An Introduction to NURBS: with Historical Perspectives[END_REF] chapter 5). Whatever the choice, this approach leads to a fully explicit space-time discretization reading

|Vi| α-1w n+1 i + S s=0 αsw (s) i = - Ω h ωi(w h ) Γn+1 ({w (s) h }) -∆t S s=0 θs K∈Ω h ∂K γ ∂K (w (s) h )[∇w h ] • [∇ϕi]
with |Vi| a nodal volume depending on the areas of the surrounding elements and on the quadrature weights induced by the finite element basis, and with the αs obtained from the "defect-correction" in time Γ -Γ.

This construction provides genuinely explicit variants of all well known stabilized continuous finite elements (SUPG, GLS, VMS, etc), as well as of nonlinear residual distribution schemes discussed in this chapter. Thorough numerical validations have been reported in [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Klosa | Extrapolated BDF residual distribution schemes for the shallow water equations[END_REF][START_REF] Klosa | Multi-step and multi-stage explicit time stepping for residual distribution[END_REF]. Some examples will be provided in the comping sections.

Applications

Scalar examples

We start with a few scalar convergence tests, to check some of the theoretical aspects discussed in the chapter. consider the approximation of solutions of the steady scalar advection equation ( 18) on the domain Ω = [0, 1] 2 , with a = (0, 1), and with inlet condition u(x, 0) = sin 2 (κπx).

We start by a result taken from [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF]. The test aims at verifying the analysis of section §3.1.7. The grid convergence has been run for κ = 1 with the nonlinear LLFs scheme [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF], and with different evaluation strategies for the streamline dissipation, or filtering term. In particular, the discrete term in ( 59) is taken as the arithmetic average of its value in a certain set of points. Note that, with the exceptions of linear polynomials, this evaluation does not give in general any kextact quadrature formula. Table 1 shows the impact of under evaluating this term. For a P 2 finite element approximation, first order of accuracy us obtained, unless a three points stencil is used. Similarly, for P 3 finite element approximation, a stencil of at least 6 points is required. Provided that the number of points is large enough, we see that indeed we recover the expected second, third and fourth order rates, even though the expressions used to evaluate the streamline dissipation are not obtained form a high order quadrature formula.

k = 1 k = 2 k = 2 k = 3
k = 3 filter : P 0 dof filter : P 0 dof filter : P 1 dof filter : P 1 dof filter: 1: Scalar advection : grid convergence for the LLFs scheme [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF]. Verification of the analysis of section §3.1.7: impact of the number of evaluation points for the "filtering term" (from [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF]).

P 2 dof h L 2 L 2 L 2 L 2 L 2 1/
The next example, taken from [START_REF] Vymazal | Very high order residual distribution schemes via variable distribution coefficients[END_REF] (see also [START_REF] Vymazal | Highorder residual distribution and error estimation for steady and unsteady compressible flow[END_REF][START_REF] D'angelo | Adjoint-based error estimation for adaptive petrovgalerkin finite element methods[END_REF]) aims at verifying the convergence rates obtained with the "variable β" LDA [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. Polynomial approximations up to degree k = 7 are tested using meshes with roughly the same number of degrees of freedom in all cases (from ≈ 2000 for the coarsest mesh to ≈ 32, 000 for the finest). The simulations are run with κ = 5. The results, summarized on table 2, show that indeed the method converges with a rate in between k + 1/2 and k + 1. For k = 6 and k = 7, converging results have only been obtained by using the optimized collocation of the degrees of freedom based on the warp and blend procedure discussed in [START_REF] Warburton | An explicit construction of interpolation nodes on the simplex[END_REF]. Computations on standard Lagrange elements with equally spaced degrees of freedom did not converge for k > 5.

External aerodynamics

In this section, we report a couple of results from [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] for compressible fluids without viscous effect (Euler equations) and [START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF][START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF] for the Navier Stokes case. The interested reader may consult [START_REF] Santis | High-order linear and non-linear residual distribution schemes for turbulent compressible flows[END_REF] for informations and results for the turbulent case (Spalart and Allmaras model). [START_REF] Vymazal | Very high order residual distribution schemes via variable distribution coefficients[END_REF], see also [START_REF] Vymazal | Highorder residual distribution and error estimation for steady and unsteady compressible flow[END_REF][START_REF] D'angelo | Adjoint-based error estimation for adaptive petrovgalerkin finite element methods[END_REF]).

Euler equations

Method: from scalar to systems So far, we have only dealt with scalar problems: the computation of the residual distribution parameters is done via arithmetic and logical operations on scalar. This cannot be as simple for systems, vecause dividing vectors has no meaning. The method that is followed has been introduced in [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF]. The idea is as follows: given an element K, we first consider an average state w. The choice of thsi average states does not seem to be essential, and we take the arithmetic mean. From this, one can evaluate the Jacobians of the flux at this state, say A(bbu). The next step is to choose a direction d. Again the choice does not seem to be essential and for fluid dynamic problems, we consider the normalised velocity except when the velocity vanishes. In that case we take an arbitrary direction. Once this is done, we compute the eigenvectors {rj}j=1,••• ,m of the matrix A(w) • d. Any vector X can be decomposed on this basis as:

X = m i=1 j (X)rj.
The eigenvectors rj are often called the right eigenvectors, while the linear forms j are often called the left eigenvectors of A(w) • d.

Starting from the LLF residual, {Φj}j=1,K where K is the number of degrees of freedom in K. For any eigenvector ri, we consider the quantities { i(Φj )}j=1,...,K that clearly satisfy:

N j=1 i(Φj ) = (Φ).
Because of this we interpret these quantities as residual and we can apply the technique of section 3.1.7 to evaluate, for any j = 1, . . . K, ( i(Φj )) = β i j (Φ) where, for example β i j is evaluated via PSI "limiter" [START_REF] Paillere | Compact cell vertex convection schemes on unstructured meshes[END_REF]. Once this is done, we define

Φ j = m i=1 ( i(Φj )) ,
this satisfies the accuracy requirements. If needed (and this generally the case), one can add a least square filtering term,

Φ j = Φ j + θ(w h )|K| x quad ω quad (A(w) • ∇ϕi(x quad )) TK (A(w) • ∇w h (x quad )
where

T -1 K = N j=1 A(w) • ∇ϕi(x quad ) .
In [151] that the matrix N j=1 A(w) • ∇ϕi(x quad ) is always invertible except when the velocity defined by w is zero. However, it that case, the matrices A(w • ∇ϕi(x quad ) TK can always be defined, see [151] for details.

Applications. These results are taken from [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF]. The meshes use triangles only unless specified.

In our first example, the domain is a square Ω = [0, 1] 2 . The boundary conditions are :

• If y > 0.5 and x = 0, the Mach number is set to M∞ = 4, the density is ρ∞ = 0.5 and the velocity is (u∞ = M∞c∞, 0) with c∞ = γp∞/ρ∞.

• If y ≤ 0.5 and x = 0, the Mach number is set to 2.4, the velocity is (u∞, 0) and the density set to 1

• The other boundary are assumed to be supersonic.

In such a configuration, the flow is steady and supersonic. We have a shock wave on the bottom, followed by a slip line and then a fan, see figure 5. Since the flow is supersonic, the x-coordinate plays the role of time : if one makes a cross-section x = const, we have a self-similar solution of the same type as what one gets for a one dimensional shock tube. It is clear that there is no oscillation at all on the density. The same conclusion holds for the other variables (not displayed).

The next example is the classical flow at M∞ = 0.35 over a sphere. In that case, the flow is symmetric with respect to the x-axis of the domain, but also with respect to the y axis. We have run this case with a second order scheme, a third order scheme, and again the second order scheme on the mesh that has the same degrees of freedom as those of the P 2 scheme. In other words, we subdivide each triangle into 4 smaller triangles which vertices are those of the large triangle and the mid-edges points. The initial mesh has 2719 nodes, 5308 elements and 100 nodes on cylinder. It is displayed on figure 6.

We see on Figure 7 which displays the pressure coefficient isolines the improvement of the solution quality when the scheme is upgraded from second order to third order. More important, the same Figure indicates clearly that the second order scheme on the refined mesh gives less accurate results than the third order one. Note that we have the same degrees of freedom in both cases.

We have re-run this test case on an hybrid mesh using the second order and the third order schemes. In both cases, the same degrees of freedom are used (i.e. we use the DOFs of the sub-triangulation for the second order scheme). The results are shown on figure 8. The mesh use 81 points on the sphere. We get the same conclusions as before.

Our next examples is a flow over a NACA012 airfoil. It is transonic, and has the following conditions at infinity: M = 0.8, angle of attack of 1.25 • The mesh has 10959 points and 21591. This corresponds to 43509 degrees of freedom. On figure 9, we have displayed the Mach number, the pressure coefficients en relative entropy deviation for the third order version of the scheme.The solutions are fine. Note however a non physical overshoot in the entropy across the upper shock.

We have run many other tests as the following (results not shown). If we compare the second order solution run with a mesh constructed from the mesh we have used where the element is sub-triangulated so that we have the same number of degrees of freedom, we can see an excellent agreement between the solutions with a main difference however. In both cases, the shock with is one element, but one element for the third order solution is roughly twice as large as an element for the second order one. Hence, the shock look more diffused in the third order case. However, the entropy levels are much lower, as we have already seen in the two sphere subsonic case.

Another case is the Ringleb flow. It has been devised by F. Ringleb [START_REF] Ringleb | Exakte Lösungen der Differentialgleichungen einer abadiatischen Gassströmung[END_REF] in 1940, see [153] for a derivation of more general solutions. This is an isentropic, irrotational two dimensional flow. It is defined from the streamline function (θ is the velocity angle with respect to a given direction and v is the norm of the velocity) ψ = sin θ v . From this, it is possible to get the explicit form of the streamlines

x = 1 2 1 ρ 1 v 2 - 2 k 2 + J 2 y = ± 1 kρv 1 - q k 2 with k = 1 φ a constant on any stream line, J = 1 c + 1 3c 2 + 1 5c 2 -1 2 log 1 + c 1 -c c = 1 -γ-1 2 q 2 , ρ = c 2/(γ-1)
The pressure is determined by the equal entropy assumption. We see that the isotach lines are the circles

x - J 2 2 + y 2 = 1 4ρ 2 q 4
From this it is possible to determine the exact solution: given a point (x, y), we determine the speed of sound c such that (x, y) belongs to the circle of center (J(c)/2, 0) and radius 1/2(ρq 2 ). Once this is done, we can get all the other values. We have run this case in the (symmetric) domain defined by

• the circle q = 0.3 on the top and the bottom,

• the extreme stream lines k = 0.4 and k = 0.8.

The simulation has been conducted with two series of meshes. The first one is made of quads cut into two triangles, always in the same direction. The mesh is then made symetric. In the second one, we only consider the quads. In both cases, we have 2 × P points on the streamlines k = 0.3 and 0.8 and P points on the circles q = 0.3. Here we have taken P = 15, 30, 60 and 100. The error (in the L 2 norm for the density are shown on figure 10. We see a slope of -3 for the third order scheme and -1.5 for the second order scheme. We also note that though the formal accuracy in both case is as expected, the effective accuracy on the quad meshes is much superior to what is obtained for triangle meshes.

We have run the same scheme on a scramjet-like configuration using an hybrid mesh as shown on Figure 11. This example has already been run in [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF]. The inflow mach number is set to 3.5. The geometry is such that many waves coexist and interact in very complex flow patterns. This situation is particularly clear on the upper part of the internal body where shocks, fans and their reflection due to wall interact. Again, in both cases, the same number of degrees of freedom have been used. Once again, the scheme has been run starting from a uniform flow configuration. Figure 12 shows the Mach number isolines. As expected, there is no real difference between the solutions since the flow is basically made of shock, fans, slip lines and constant states : this is not an accuracy case, but a case that shows that, despite the flow complexity, the third order scheme is robust. However, one can see a small difference between the solutions : the slip line created by the interaction of two shocks after the blade is a little bit more twisted for the third order scheme than the second order one. We also see that the resolution of the discontinuities is in both case approximately one cell width.

Navier Stokes equations

We rereport here results taken from [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF]. The scheme and problems have already been discussed. For more details, the reader may also consult [START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF]. The filtering term has to be more elaborate in order to take into account the viscous terms.

The first example is the classical test case consisting of a subsonic viscous flow over a NACA-0012 airfoil at zero angle of attack. The free stream mach number is 0.5 and Reynolds number is 5 000. This is a widely used test case for two dimensional laminar flows; a distinctive feature of this test case is a steady separation bubble near the trailing edge of the airfoil. An example of computational grid is displayed in figure 13. The grid extends about 50 chords away from the airfoil. The airfoil boundary is considered adiabatic, no-slip and is represented by piece-wise quadratic elements, the far-field boundary condition is applied on the outer boundary of the domain, see [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF] for a precise description of the boundary conditions approximation, as well as details on the steady state solver. The steady state is considered to be reached when the L 2 norm of the density residual is drop by ten orders of magnitudes respect to the initial value.

In figure 14 are depicted the solutions computed with the linear scheme and the SPR-ZZ gradient reconstruction, for P 1 and P 2 elements. The solution with the Φ 1 elements has been computed on a grid obtained from that with P 2 elements (4 216 elements) and splitting each P 2 triangle with four P 1 triangles, in such a way the number of DOFs for the second and third order simulation is exactly the same. Note, in figure 14, that although there is not much difference in the Mach number contours between the second and the third order simulations, the streamlines near the trailing edge are very different, and only the third order scheme is able to reproduce the symmetric recirculation bubble. For the same simulations, in figure 15 and 16 are reported the pressure and skin friction coefficients profiles, respectively. Note the more regularity of the solution of the third order simulation respect to the second order one, for the same number of DOFs.

The second example is a steady, laminar flow at high angle of attack, around a delta wing with sharp edges. As the flow passes the leading edge, it rolls up and creates a big vortex structure which is convected far behind the wing, at the same time, near the leading edge a smaller secondary vortex appears. A free stream Mach number M = 0.5 is considered, the Reynolds number, based on the root chord of the wing is Re = 4 000, the angle of attack is α = 12.5 The geometry of the delta wing is depicted in figure 17, together with an example of a coarse grid used for the simulations. The grid consists of tetrahedra; finer levels of grids are obtained by uniformly splitting each tetrahedron of the coarser level with eight tetrahedra. Note the presence of very stretched elements on the wing. The wing surface is treated as no-slip adiabatic wall, the vertical plane intersecting the root of wing is treated as a symmetry plane, while far field boundary conditions are applied on the outer boundary of the domain.

The solution is initialized with an uniform flow, the lower order solution is used as initial solution for the third order computation. For this test case the linear scheme is used with the SPR-ZZ gradient reconstruction method; in figure 18 are reported the streamlines and Mach number contours, at different stations, of the third order solution on the finest grid.

In figure 19 are reported the drag and lift coefficients computed with linear and quadratic elements, on three uniformly refined grids. For comparison, are reported also the reference values computed in [START_REF] Leicht | Error estimation and anisotropic mesh refinement for 3d laminar aerodynamic flow simulations[END_REF] by extrapolating the results obtained with a higher order DG method. Observing the convergence of the drag coefficient in term of DOFs, it can be noted that there is no significant gain in using a higher order approximation, with respect to the second order. This behavior can be caused by the singularity at the leading edge of the wing, which might mask the benefits of a higher order approximation with an uniform mesh refinement. Regarding the convergence of the lift coefficient, it could be observed a clear benefit of using a higher order approximation, because the big vortex structure over the wing is better captured with higher order elements.

As last test example, the interaction of an oblique shock wave with a laminar boundary layer is considered. The aim of this test is to show the non-oscillatory properties of the non-linear scheme in presence of discontinuities of the solution and at the same time, the capability to maintain the accuracy required for the discretization of the boundary layer.

The test consists in a laminar boundary layer developing over a flat plate and an incident shock impinging the boundary layer. Since the flow is supersonic, a shock appears at the leading edge of the flat plate, that interacts with the oblique shock. Furthermore, at the impinging point, the incident shock produces a separation of the boundary layer, the shock is then reflected and an expansion fan appears, turning the flow toward the wall and causing a reattachment of the boundary layer, as it is depicted in figure 20.

In the numerical simulations, the oblique shock is generated by imposing the incoming supersonic flow state on the lower part of left boundary, while another supersonic state is imposed on the upper part of the left boundary and on the top boundary; this state is computed using the relations of the oblique shocks, such that the incident shock has a certain angle of incidence θs. The height of the computational domain is 0.94, while the range of the computational domain in the x direction is [-0.2, 2], the flat plate has length L = 2 with the leading edge of flat plate at x = 0. Along the plate, the no-slip adiabatic wall boundary condition is applied, while the symmetry boundary condition is applied on the remaining part of the bottom boundary. On the right boundary, the outflow boundary condition is applied, see figure 20. The inflow states are chosen such that the free-stream Mach number is M = 2.15 and the angle of the incident shock is θs = 30.8 • , in this case the impingement point would be at center of the plate for an inviscid fluid. The Reynolds number based on the free-stream values and the distance between the plate leading edge and the inviscid shock impingement point is 1 × 10 5 .

The non-linear scheme with the SPR-ZZ gradient recovery strategy is used to perform the numerical simulations at second and third order of accuracy. The computational domain is generated from the triangulation of a 90 × 85 structured grid; the first number refers to the number of elements on the horizontal boundaries, with 80 elements along the plates, the second number refers to the number of elements on the vertical boundaries. The element distribution is uniform on the x direction, while along the y direction a non-uniform distribution of the elements is used, with a mesh spacing ∆y = 0.5 × 10 -3 near the bottom boundary. For comparison, a second order simulation is also performed on a finer grid with the same number of DOFs of the third order simulation on the coarse grid. The simulation is initialized with an uniform solution, and the second order solution is used as initial solution for the third order approximation. Except the case of the second order simulation on the coarse grid, for which the initial residual is reduces by ten orders of magnitude, the residual for the third order simulation and the second order one on the finer grid could not be reduced by more than eight orders of magnitude.

In figure 21-(a) are shown the contours of the pressure for the third order simulation; all the features of this problem are well represented. In figure 21-(b) is reported a zoom of the solution where the incident shock impinges the boundary layer. Two features are evident: the reflection of the incident shock and the recirculation bubble as a consequence of the separation and subsequent reattachment of the boundary layer produced by the incident shock and the expansion fan, respectively.

The profiles of density, pressure and Mach number along the lines at y = 0.29 and y = 0.15 are reported in figure 22. Note that the third order scheme gives a very sharp and monotone representation of the discontinuities and also smooth portions of the solution are better represented compared to the second order solution. It is important to remember that smooth and discontinuous solutions are treated within Figure 21: Left: contours of the pressure obtained with the third order scheme for the shock/boundary layer interaction. Right: zoom of the solution near the impinging point of the shock with the boundary layer, streamlines are also reported to show the separation bubble.

the same non-linear scheme without any special treatment or tuning parameter. For a fair comparison, it is also reported the solution obtained with the second order scheme on a finer mesh; it is worth noticing that, although a mesh refinement produces an improvement of the numerical solution, the level of accuracy obtained with the second order scheme is still lower than that obtained with the third order scheme, for the same number of DOFs. Finally, in figure 23 are reported the values of the pressure and of the friction coefficient along the plate. The oscillations near the point x = 0 are due to the singularity of the solution at the leading edge of the flat plate, but they are limited only in small region around the leading edge. The third order scheme seems less sensitive to this singularity compared to the second order simulations. The separation bubble can easily detected by the negative values of friction coefficients, note also the pressure plateau in the detached zone.

Free surface flows

The framework presented in this chapter has proved quite interesting to construct discrete approximations of systems of PDEs modelling free surface flows, namely the shallow water equations and dispersive enhancements (Boussinesq and/or Green-Naghdi equations). Early work on steady hydrostatic flows had been done in the PhD of M. Hubbard (see e.g. [START_REF] Garcia-Navarro | Genuinely multidimensional upwinding for the 2d shallow-water equations[END_REF][START_REF] Hubbard | Conservative multidimensional upwinding for the steady twodimensional shallow water equations[END_REF], and also the paper by Brufau and Garcia-Navarro [START_REF] Brufau | Unsteady free surface flow simulation over complex topography with a multidimensional upwind technique[END_REF]). More recent work, combining high order of accuracy in space and time, the preservation of moving steady states, robust handing of dry areas, and dispersive extensions is discussed in [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | On the c-property and generalized c-property of residual distribution for the shallow water equations[END_REF][START_REF] Sarmany | Unconditionally stable space-time discontinuous residual distribution for shallow water flows[END_REF][START_REF] Ricchiuto | Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF].

Inundation of a complex three-dimensional beach

The first example we consider involves the solution of the shallow water equations, and is taken from [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]. In this paper, the nonlinear stabilized Lax-Friedrich's method has been combined with the fully explicit time marching strategy discussed in section §3.2, and modified to allow a (provable) preservation of the non-negativity of the water depth. To illustrate the capabilities of the method obtained we consider a standard benchmark in the oceanography community involving the Tsunami runup onto a complex threedimensional beach. The so-called Monai valley benchmark aims a simulating a scaled down laboratory experiment reproducing the impact of the tsunami wave that hit the Okushiri island in Japan in 1993. The bathymetry, and inlet data are available on the web page of the third international workshop on long wave runup models [160] (see also [161,[START_REF] Liu | Advanced Numerical Models for Simulating Tsunami Waves and Runup[END_REF]), with the data relative to the time series of the wavet level in three gauges close to the shore. The shape of the bathymetry and of the inlet wave, as well as the position of the three wave gauges are shown in the left, middle and right pictures in figure 24. In the observations [160,161,[START_REF] Liu | Advanced Numerical Models for Simulating Tsunami Waves and Runup[END_REF] the highest runup is of 32 [m], and it occurs in the region of the Monay valley where the bathymetry is steepest. For clarity, this region is encircled in the results presented in the following.

Approximation of moving steady states

The super consistency property discussed in section §3.1.8 also has applications in shallow water flows. In this case, the state vector w is defined by the quantities H, the water depth, and q, the volume flux q = H u, with u the depth-averaged flow velocity. A known steady state involving moving water, is the pseudo-one dimensional flow characterized by q = q0 = c t , and E = E0, with E the total energy g(H + b) + u • u/2, and b = b(x, y) the bathymetry. This solution allows to check numerically proposition 3.9. To do this, we consider the tests discussed in [START_REF] Ricchiuto | On the c-property and generalized c-property of residual distribution for the shallow water equations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF].

The first, involves a small perturbation of the initial steady state over a bathymetry with C 1 regularity obtained as a series of ribs defined by truncated sin 2 functions. The evolution of the perturbation on an irregular triangulation is studied. The typical result is shown on 27 showing a 3D view of the free surface level. The left picture is obtained with a standard scheme based on a P 1 approximation of the state vector w and of the bathymetry. The right result is instead obtained with the scheme based on a direct approximation of the total energy E and of the flux q, and with a higher order approximation and quadrature of the bathymetric gradient. The improvement is quite remarkable.

The second test consists in verifying the property of proposition 3.9 by computing, on irregular triangulations, the solution error at a finite time when starting from the exact nodal steady state. This is done with bathymetries of increasing smoothness, and with volume and edge quadrature strategies of increasing accuracy. The results are sumamrized on figure 28 the grid convergence obtained on bathymetries with different regularity when using quadrature strategies with errors of orders h 2 , h 4 , h 6 , and h 8 respectively. The last column show the error convergence on a fixed mesh when increasing the accuracy of the quadrature. In articular, picture (e) is obtained on the coarsest mesh used in the convergence study, while picture (f) on the finest. The underlying approximation is P 1 . Not only this result confirms the super consistency analysis, but it also shows that for exact quadrature, the residual approach would yield exact preservation of the steady state.

For additional examples involving other steady state solutions, the interested reader is referred to [START_REF] Ricchiuto | On the c-property and generalized c-property of residual distribution for the shallow water equations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Ricchiuto | Contributions to the development of residual discretizations for hyperbolic conservation laws with application to shallow water flows[END_REF].

Residual based stabilized methods for dispersive waves

Another challenging application in free surface flows in the inclusion of non-hydrostatic effects in depthaveraged models. The interested reader may consult the review papers [START_REF] Kirby | Chapter 1 boussinesq models and applications to nearshore wave propagation, surf zone processes and wave-induced currents[END_REF][START_REF] Brocchini | A reasoned overview on boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF] and the book [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] for an overview of the modelling issues. Concerning numerics, the typical form of a depth averaged Boussinesq-type model is ∂tK + ∇ • f (w) + s(w, x) = 0 [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF] where the quantity K(w) is related to the state vector by w -T (w) = K

where T (•) is a non-linear elliptic operator. Here physical dispersion is present in the PDE. The challenge is thus to design a numerical method with low dissipation and very low dispersion errors to allow long time integration of propagating waves, while however guaranteeing a sufficient degree of dissipation to avoid spurious modes. The use of some stability mechanism is also required as the term T (w) is often neglected locally to recover the hyperbolic shallow water equations, and model breaking regions as moving bores [START_REF] Tonelli | Simulation of wave breaking over complex bathymerties by a Boussinesq model[END_REF][START_REF] Bonneton | A spliting approach for the fully nonlinear and weakly dispersive green-naghdi model[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF]. The requirement is then to have a low dissipation/dispersion method, capable of handling both the parabolic Boussinesq equations, and the hyperbolic shallow water ones, with eventually capabilities for capturing of shocks and dry areas. This has led to the work presented in [START_REF] Ricchiuto | Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF][START_REF] Bacigaluppi | A 1d stabilized finite element model for nonhydrostatic wave breaking and run-up[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF] which has tried to extend upwind and multidimensional upwind residual based stabilization techniques to these systems. The main idea is to decouple the approximation of the two sub problems above. The elliptic step ( 92) is solved with a standard C 0 Galerkin method, while an upwind scheme is used in the evolution step [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF]. The work discussed in the references shows evidence that this approach is a sound one, and that provided that the hyperbolic step is solved with al least third order of accuracy, the elliptic phase can be solved with a second order method without affecting the dispersion accuracy. This generalizes on unstructured grids, and to residual based stabilized method, an idea proposed in the finite difference context by Wei and Kirby in [START_REF] Wei | A time-dependent numerical code for extended Boussinesq equations[END_REF]. The schemes obtained, all reduce in 1D to a streamline upwind method stabilizing the Galerkin approximation of the first order PDE [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF] with cell integrals depeding on the residual of (91), and on the sign of the shallow water Jacobians. In two space dimensions, both a standard stramline upwind formulation, and a multidimensional upwind variant based on the LDA method (61) have been proposed in [START_REF] Ricchiuto | Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF].

We present here three results. The first is the characterization of the accuracy of the methods obtained. Figure 29 shows results relative to a second order Galerkin approximation of [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF], and a third order stramline upwind (SU) or or fourth order Galerkin (cG) approximation of [START_REF] Abgrall | A simple construction of very high order nonoscillatory compact schemes on unstructured meshes[END_REF]. In particular, the left picture provides a numerical convergence study on a propagating solitary wave. Despite the second order treatment of the elliptic term, the overall accuracy measured for a propagating solution is three. More importantly, the middle and right pictures study the dispersion errors of the schemes and compare them to second and fourth order finite differencing. The result shows two important features: both the cG and SU are as good or better than the fourth order finite difference method; for propagating solutions, the upwind SU stabilization actually improves the dispersion properties of the scheme providing lower dispersion errors, especially for shorter waves.

The second result tests the ability of the proposed method to correctly reproduce the energy exchange between different harmonics when monochromatic waves shoal on a two dimensional circular shelf. This is a standard benchmark for multidimensional Boussinesq-type codes (see [START_REF] Sørensen | A new form of the boussinesq equations with improved linear dispersion characteristics. part 2: a slowing varying bathymetry[END_REF][START_REF] Beji | A formal derivation and numerical modeling of the improved boussinesq equations for varying depth[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended boussinesq equations[END_REF][START_REF] Sørensen | Boussinesq-type modelling using an unstructured finite element technique[END_REF][START_REF] Eskilsson | An unstructured spectral/hp element model for enhanced boussinesq-type euations[END_REF][START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2dh improved boussinesq equations[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2d boussinesq-type equations[END_REF] and references therein). In [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF] experiments have been conducted in several configurations involving values of period and amplitude for the incoming monochromatic wave. Here we discuss the results for: case (a) with T=2 s, A = 0.0075 m, h0/λ = 0.117; case (b) with T=1 s, A = 0.0195 m, h0/λ = 0.306. The first case has a relatively weak dispersive character, but presents an important energy transfer to higher harmonics. The second case is quite demanding as it involves a higher dispersion degree, outside the validity of the most simple Boussinesq models. Figure 30 summarizes the results obtained solving the enhanced Boussinesq equations of [START_REF] Schaffer | Further enhancements of boussinesq-type equations[END_REF] on unstrucutred triangulations. Both a 'classical' streamline upwind stabilization and a multidimensional one based on the LDA distribution (61) have been tested. The pictures clearly show that these multidimensional stabilized methods have a high potential in resolving the energy transfer between harmonics, also in the more demanding cases.

Finally, we show the results obtained on an experiment carried on [START_REF] Berkhoff | Verification of numerical wave propagation models for simple harmonic linear water waves[END_REF] and involving the refraction and diffraction of monochromatic waves over a complex bathymetry. A sketch of the experiment is reported on the top-left picture on figure 31. The bathymetry involves a shoal presenting a constant angle with the main incoming wave direction, with an elliptic bump which leads to a complex multidimensional wave pattern which involves dispersive effects both in the main wave direction, and along the orthogonal. As shown in the sketch on figure 31, the experiments provide the the normalized time-average of the water height in 8 different sections. Profiting from the general formulation used here, the problem is solved on an unstructured triangulation refined in correspondence of the sampling region, as shown in the picture on the top-right on figure 31.

On the same figure, the typical instantaneous wave pattern obtained is shown. One can clearly see the effect of the submerged feature in diffracting the incoming waves. To provide a more quantitative appreciation of the result, the comparison with the experiments is shown for thee of the eight sections on figure 32. The results, again obtained with two different upwind (and multidimensional upwind) stabilization approaches, confirm the potential of residual based methods in capturing complex dispersive wave phenomena.

Conclusion, Open challenges

Over the years, the Residual distribution technology has proved that continuous finite elements allow the same flexibility as discontinuous finite elements. The stencils are comparable, in particular for viscous calculations, and less degrees of freedom are always needed, eventhough the difference between completely discontinuous approxiamtion and continuous ones tend to become smaller and smaller as the polynomial degree increases. We have also shown that all these methods are locally conservative contrarily to a common belief. The techniques developed here shows that the schemes is very robust. We could not show all possible results, but simulation for hypersonic flows are possible without major difficulties. We have also shown (see reference), that iterative convergence to machine zero is possible even for turbulent flows, see [START_REF] Santis | High-order linear and non-linear residual distribution schemes for turbulent compressible flows[END_REF].

Figure 2 :

 2 Figure 2: Residual distribution

Proposition 3 . 7 (

 37 Local Positivity and discrete maximum principle). Locally positive schemes, combined with the evolution step (23) verify the discrete maximum principle

Figure 3 :

 3 Figure 3: Node centered finite volume

Figure 7 :

 7 Figure 7: Subsonic sphere problem : Isolines of the pressure coefficient. We have the same isolines on each figure.

Figure 8 :

 8 Figure 8: Subsonic sphere problem, hybrid mesh : Pressure coefficient and entropy variation on an hybrid mesh, M ∞ = 0.35.

Figure 22 :Figure 22 :Figure 22 :Figure 22 :Figure 9 :

 222222229 Figure 22: Transonic NACA012 problem. Isolines of the Mach number, pressure, density and entropy for the NACA012 case.

Figure 24 :

 24 Figure 24: Ringleb flow problem. L 2 error on the density for the Ringleb flow. Tri stands for triangle, Quad for quadrangle. O2 stands for second order, O3 for third order.

Figure 10 :

 10 Figure 10: Ringleb flow problem. L 2 error on the density for the Ringleb flow. Tri stands for triangle, Quad for quadrangle. O2 stands for second order, O3 for third order.

x 10 11 Figure 25 :

 1125 Figure 25: Zoom of the mesh for the scramjet problem.

Figure 26 :

 26 Figure 26: Scramjet problem. Mach number distribution. Top : the third order solution, bottom the second order solution. The same isolines are plotted.

Figure 11 :

 11 Figure 11: Zoom of the mesh for the scramjet problem.

x 10 11 Figure 25 :

 1125 Figure 25: Zoom of the mesh for the scramjet problem.

Figure 26 : 10 11 Figure 25 :

 26101125 Figure 26: Scramjet problem. Mach number distribution. Top : the third order solution, bottom the second order solution. The same isolines are plotted.

Figure 26 :Figure 12 :

 2612 Figure 26: Scramjet problem. Mach number distribution. Top : the third order solution, bottom the second order solution. The same isolines are plotted.

Figure 13 :

 13 Figure 13: An example of computational grid used for the NACA-0012 test case.

Figure 14 :

 14 Figure 14: Mach Number contours (top) and streamlines near the trailing edge (bottom) for the second (left) and third (right) order linear scheme.

Figure 15 :Figure 16 :Figure 17 :

 151617 Figure 15: Pressure coefficient along the whole NACA-0012 airfoil for the second and third order simulations, with the same number of DOFs.

Figure 18 :Figure 19 :Figure 20 :

 181920 Figure 18: Streamlines and slices of Mach number contours along and behind the delta wing, for a third order simulation on a fine grid.

Figure 22 :

 22 Figure 22: Density, pressure and Mach number profiles along the line y = 0.29 (a, c, e) and the line y = 0.15 (b, d, f) for the shock/boundary layer interaction problem.

Figure 23 :

 23 Figure 23: Pressure (a) and skin friction (b) profiles along the flat plate for the shock/boundary layer interaction problem.

Figure 24 :

 24 Figure 24: Monai valley benchmark. Left: bathymetry. Middle: incoming wave. Right: position of the wave gauges.

Figure 25 :

 25 Figure 25: Monai valley benchmark. 3D view of the inundation process.

Figure 26 :

 26 Figure 26: Monai valley benchmark. Time series in gauge 5 (left), and runup plot (right).

  in which figures (a), (b), (c), and (d) show

Figure 27 :

 27 Figure 27: Moving steady states: evolution of a small perturbation in a homo-energetic steady state. 3D plot of the free surface. Left: approximation in physical variables. Right: approximation in steady invariants.

Figure 28 :

 28 Figure 28: Moving steady states: super consistency of the scheme. Left and middle column: grid convergence for different quadrature strategies. Right column: quadrature convergence on the coarsest (top) and finest (bottom) grid.

Figure 29 :

 29 Figure 29: Accuracy of residual approximation of dispersive equations. Left: grid convergence study on an exact solitary wave solution. Middle and left: dispersion error in function of the number of nodes per wavelength for the upwind stabilized and unstabilized scheme, and comparison with second and fourth order finite differencing. Solid line: kh 0 = 0.5 (long wave). Circles: kh 0 = 2.6 ("short" wave).

Figure 30 :

 30 Figure 30: Wave diffraction on a circular shoal, case (a) (left) and case (b) (right).

Figure 31 :

 31 Figure 31: Wave scattering on an elliptic shoal. Problem sketch (top-left), close up of the grid (top-right), and instantaneous wave patterns (bottom).

Figure 32 :

 32 Figure 32: Wave scattering on an elliptic shoal. Time-average of the water height: simulations vs experiments.

Table

	25	0.50493E-02	0.25122E-01	0.32612E-04	2.17274E-02	0.12071E-05
	1/50	0.14684E-02	0.12935E-01	0.48741E-05	1.13486E-02	0.90642E-07
	1/100	0.41019E-03	0.83978E-02	0.66019E-06	5.83347 E-03 0.53860E-08
	average	1.790	0.7904	2.812	0.9292	3.914
	rate

Table 2 :

 2 Scalar advection : convergence for the variable β LDA (61) (courtesy of M. Vymazal

	k	N dof	h	L 2	rate
	1	2094 0.02185 3.49E-02	-
		8124 0.01109 7.44E-03 2.24
		32546 0.00554 1.36E-03 2.46
	2	2189 0.02137 1.37E-02	-
		8217 0.01103 2.19E-03 2.65
	(equi-spaced) 32181 0.00557 3.04E-04 2.88
	3	2113 0.02175 5.20E-03	-
		8347 0.01095 2.94E-04 4.16
	(equi-spaced) 33520 0.00546 2.11E-05 3.89
	4	2017 0.02227 2.57E-03	-
		8593 0.01079 9.94E-05 4.71
	(equi-spaced) 32553 0.00554 4.28E-06 4.55
	5	2381 0.02049 1.15E-03	-
		8611 0.01078 2.83E-05 5.36
	(equi-spaced) 33546 0.00546 5.32E-07 5.75
	6	2317 0.02077 6.68E-04	-
		8293 0.01098 7.30E-06 7.01
	(warp-blend) 33073 0.00550 7.29E-08 6.67
	7	2633 0.01949 3.82E-04	-
		9430 0.01030 2.44E-06 7.92
	(warp-blend) 34427 0.00539 9.86E-09 8.51

this aspect can be explained in more rigorous terms and made systematic if taken into account from the start

A similar explicit construction can be done in the three space dimensions including high prder tets, prisms, and hexas. Details are left out of this manuscript

The results above can easily be generalised to the conditions 54 provided f is Lipschitz continuous.

having set here T K =diag(τa, δν , δν )

Acknowledgements

This work would not have been possible without the contributions, suggestions and help of our collaborators, students, and friends (in alphabetical order): P. Congedo (INRIA), H. Deconinck (VKI), S. D'Angelo (VKI, now Credo Consulting), D. De Santis (INRIA, now Stanford University), A. Filippini (Inria) M. Hubbard (Nottingham University), A. Larat (INRIA, now Ecole Centrale Paris), M. Vymazal (VKI, nom Imperial College) and many others. RA has been conducting this work with the partial support of SNFS grant # 200021 153604. MR has been partly funded by the TANDEM project (reference ANR-11-RSNR-0023-01 French Programme Investissements d'Avenir.

which the node belongs. As before, this leads to a system of equations, which reads :

second order third order Figure 5: Jet problem : isolines of the density, second and third order LLxFf scheme. All the degrees of freedom are plotted. and the same isolines are plotted Figure 6: Subsonic sphere problem : Zoom of the mesh for the sphere problem. The mesh has no symetry.

However, all the problems have not been solved so far:

• High order and unsteady problems: this chapter has displayed a couple of solutions for geophysical flows. Other examples, related to compressible flow problems, can be found in [START_REF] Ricchiuto | Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case[END_REF] where a fully explicit method is described, or [START_REF] Sarmany | Unconditionally stable space-time discontinuous residual distribution for shallow water flows[END_REF][START_REF] Abgrall | Towards very high-order accurate schemes for unsteady convection problems on unstructured meshes[END_REF] for implicit technique. Considering now higher than second order in time, research is still needed, however see [START_REF] Abgrall | How to avoid mass matrix for linear hyperbolic problems[END_REF] for a fully explicit (i.e. mass matrix free) technique for linear problems. The same technique can be applied for non linear problems.

• Error estimation and adaptation. Some work on adjoint problems in the RD framework has been done by [START_REF] D'angelo | Adjoint-based error estimation for adaptive petrovgalerkin finite element methods[END_REF][START_REF] D'angelo | Generalized framework for adjoint error estimation of PG method in linear problems[END_REF].

• p-adaptation : continuous and discontinuous approximation. Some work in that direction has been done in [START_REF] Abgrall | p-adaptation is possible with continuous finite elements: the Euler equations case[END_REF][START_REF] Abgrall | p-adaptation is possible with continuous finite elements: the Navier Stokes equations case[END_REF] [