Keywords: numerical analysis, second order schemes, hyperbolic problems, residual distribution, explicit schemes numerical analysis, second order schemes, hyperbolic problems, residual distribution, explicit schemes

Explicit Runge-Kutta Residual Distribution schemes for Time Dependent Problems : second order case

Résumé : In this paper we construct spatially consistent second order explicit discretizations for time dependent hyperbolic problems, starting from a given Residual Distribution (RD) discrete approximation of the steady operator. We explore the properties of the RD mass matrices necessary to achieve consistency in space, and finally show how to make use of second order mass lumping to obtain second order explicit schemes. The discussion is particularly relevant for schemes of the residual distribution type which we will use for all our numerical experiments. However, similar ideas can be used in the context of residual based finite volume discretizations.

Introduction

The aim of this study is to understand, given a residual based discretization of steady limit of a hyperbolic conservation law, how to contruct fully explicit consistent discretizations for time dependent problems. As a case study we consider schemes of the residual distribution (RD) type [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF][START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF]. While well understood in the steady case, their formulation in the time dependent case has never been completely clarified, due to the lack of a rigorous formulation allowing a natural extension. In particular, the lack of sufficient constraints on the discretization has led in time to a number of different formulations all featuring different mass matrices [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for navier stokes simulations[END_REF][START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF][START_REF] Deconinck | Status of multidimensional upwind residual distribution schemes and applications in aeronautics[END_REF][START_REF] Palma | A secondorder accurate monotone implicit fluctuation splitting scheme for unsteady problems[END_REF]. In this paper we show that all these formulations are equivalent, up to a dissipation operator. We also show the existence of entire families of additional consistent mass matrices.

Finally, we show how to combine all of the above discretization with high order (second order) mass-lumping to obtain fully explicit schemes. Note that in the case of nonlinear RD discretizations based on second order time integration in time, positivity preservation is obtained only under an explicit CFL-type condition [START_REF] Abgrall | Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. This fact, related to the properties of the underlying ODE integrator [START_REF] Bolley | Conservation de la positivité lors de la discétization des problèmes d'èvolution paraboliques[END_REF], and the highly implicit nature of the schemes leads to poor efficiency. The explicit formulation propose here is one possible solution to this issue.

The work discussed here somehow generalizes the initial work of [START_REF] Rossiello | Second-orderaccurate explicit fluctuation splitting schemes for unsteady problems[END_REF] where only central Lax-Wendroff type discretizations are considered. Moreover, in the paper we will show that the ideas presented here also apply to other classes of schemes, such as the ones proposed in [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF].

The structure of the paper is as follows. We start by introducing the notations used throughout the discussion. Then, in section 3 we review different formulations of RD for time dependent problems. We discuss their relations, and show how other families of consistent formulations exist. We elaborate further on these ideas in 4 where we finally show how to obtain fully explicit formulations which still retain the same formal accuracy of the fully implicit formulations. In section 5 we give a summary of the different numerical schemes that we use in the numerical tests which are discussed in sections 6 and 7. We end the paper with some conclusive remarks and some thoughts for further developments.

Mathematical problem and notation

We seek approximations of solutions of the time dependent hyperbolic problem

r(u) = 0 , r(u) = ∂ t u + ∇ • F (u) (1)
on some spatial domain Ω, and on some temporal domain [0, t f ]. We will mainly focus on the two-dimensional case Ω ∈ R 2 , but the generalization to three spatial dimension is trivial. We discretize Ω by an unstructured triangulation denoted by T h , with T denoting the generic element of the mesh, and h the mesh parameter (characteristic mesh size). When no confusion is generated we denote the nodes of T by {1, 2, 3}. In every element, we denote by n j the inward pointing vector normal to the edge facing node j, scaled by the length of the edge. Denoting by ϕ i the P 1 Lagrange basis function corresponding ot node i ∈ T h , we have

∇ϕ i T = n i 2|T | (2)
The P 1 approximation of u will be denoted by u h , and it is given by

u h = i∈T h u i ϕ i = T ∈T h j∈T u j ϕ j T (3)
The temporal domain is discretized by a set of non-overlapping time slabs [t n , t n+1 ]. We denote by ∆t = t n+1 -t n the time step.

To simplify the presentation of the next sections we also introduce here the element fluctuation defined as

φ(u h ) = T ∇ • F h (u h ) dx dy = ∂T F h (u h ) • n dl , (4)
INRIA the element residual

Φ(u h ) = T r(u h ) dx dy = T (∂ t u h + ∇ • F h (u h )) dx dy = j∈T |T | 3
du j dt + φ(u h )
(5) and the local Galerkin residual

φ G i (u h ) = T ϕ i ∇ • F h (u h ) dx dy (6)
In the expressions above u h represents the P 1 numerical approximation of the unknown, and F h (u h ) a discrete approximation of the flux. Note that all of the above quantities depend on time. Moreover, to simplify the notation, we do not introduce a super-or sub-script indicating the element T over which they are evaluated, this being always clear from the context.

Second order RD : the proliferation of mass matrices

Let us for the moment consider the particular case of [START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF] given by the linear constant advection problem

∂ t u + a • ∇u = 0 (7)
We focus our attention on discrete counterparts of (7) that, on a slab T h × [t n , t n+1 ] can be written as

T |i∈T    j∈T m T ij du j dt + β i φ(u h )    = 0 ∀ i ∈ T h (8)
Last definitions give a scheme that, requires the solution of a (generally) nonlinear system if the mass matrix m ij is non-diagonal. Moreover, introducing the nodal residuals

Φ i (u h ) = j∈T m T ij du j dt + β i φ(u h ) (9)
we also require that, j∈T Φ j (u h ) = Φ(u h ) [START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF] with Φ(u h ) given by [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for steady flow problems[END_REF]. The prototype [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for navier stokes simulations[END_REF] is meant to be a consistent generalization to the time dependent case of the fluctuation splitting/residual distribution discretization which approximates the steady limit of [START_REF] Bolley | Conservation de la positivité lors de la discétization des problèmes d'èvolution paraboliques[END_REF] as

T |i∈T β i φ(u h ) = 0 ∀ i ∈ T h (11)
with j∈T β j = 1 [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF] In order to distinguish the steady advective operator from the time dependent equation, we have chosen to keep a distinction between the fluctuation (4) and the residual [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for steady flow problems[END_REF], the latter representing the integral of the whole equation. Historically, the first consistent approaches to obtain such a generalization were based on two different points of view. In the first approach [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], one simply replaces the fluctuations in the discrete equations [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF] with the full residual [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for steady flow problems[END_REF]. This residual is then distributed exactly as in [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF] (with the LDA scheme in the original reference [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][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]), leading to

0 = T |i∈T β i Φ(u h ) = T |i∈T   j∈T m F1 ij du j dt + β i φ(u h )   , m F1 ij = |T | 3 β i ( 13
)
with δ ij Kroenecker's delta, and F1 standing for Formulation 1.

A second approach [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] uses an analogy with stabilized Galerkin finite element schemes in which the discrete equations [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF] are obtained as

T |i∈T β i φ(u h ) = T |i∈T φ G i (u h )+ T |i∈T δφ i = Ω ϕ i a•∇u h dx dy+ T |i∈T T δ ϕi a•∇u h dx dy
with the perturbation to the test function δ ϕi depending on the distribution coefficients β T i . In particular, for constant advection and a P 1 variable approximation, one can assume δ ϕi to be a constant, to find easily (in two space dimensions) δ ϕi | T = β i -1/3. In the time dependent case this naturally leads to 0 = Ω ϕ i r(u h ) dx dy + T |i∈T T δ ϕi r(u h ) dx dy

= T |i∈T   j∈T m F2 ij du j dt + β i φ(u h )   , m F2 ij = |T | 36 ( 3
δ ij + 12β i -1) (14)
with δ ij Kroenecker's delta, and F2 standing for Formulation 2.

A different approach has instead been proposed in [START_REF] Palma | A secondorder accurate monotone implicit fluctuation splitting scheme for unsteady problems[END_REF], where the authors use the idea that in every T ∈ T h the combination of terms arising from the multiplication of the mass matrix with the nodal time derivatives should give back an integral of the time derivative of u h over a dual sub-element C j ∈ T . Consistency is guaranteed by the requirement |C j | = β j |T | 1 . In particular, in the paper the authors require the node j to belong to C j (cf. figure 1). Conditions for second order of accuracy are shown to be

i∈T m ij = |T | 3 , j∈T m ij = |T |β i (15)
In the reference, the authors ultimately arrive to the following formulation

0 = T |i∈T Ti r(u h ) dx dy = T |i∈T   j∈T m F3 ij du j dt + β i φ(u h )   , m F3 ij = |T | 3 β i (δ ij + 1 -β j ) (16)
1 which implicitly assumes β i ≥ 0 ∀i

INRIA

with δ ij Kroenecker's delta, and F3 standing for Formulation 3. This formulation is actually proposed mainly for schemes with a multidimensional upwind characher [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF], and in particular for those schemes that, when the advection speed points toward an edge of element T , give β j ≥ 0 only if j belongs to this edge. In this case, only the rows of m F3 ij relative to these nodes contain non-zero elements.

The idea of [START_REF] Palma | A secondorder accurate monotone implicit fluctuation splitting scheme for unsteady problems[END_REF] can actually be used to derive still another member to the family of consistent mass matrices. It suffices to follow the exact same developments done in the reference, except that we allow the sub-triangle j not to contain node j itself. In particular, whenever β i ≥ 0 ∀ i, we note that we can find a unique point, say M ∈ T , such that ϕ i (M ) = β i ∀ i. This means that the β i coefficients represent the area coordinates of M (see figure 1). With the notation of figure 1, using the fact that u h (M ) = β 1 u 1 + β 2 u 2 + β 3 u 3 we find in the time dependent case :

0 = T |i∈T Ti r(u h ) dx dy = T |i∈T   j∈T m F4 ij du j dt + β T i φ(u h )   , m F4 ij = |T | 3 β i (1 -δ ij + β j ) (17)
with δ ij Kroenecker's delta, and F4 standing for Formulation 4. Note that the difference with respect to the matrix proposed in [START_REF] Palma | A secondorder accurate monotone implicit fluctuation splitting scheme for unsteady problems[END_REF] is that here j / ∈ T j . One easily checks that conditions (15) are verified. All of the above construction can be still generalized by noting that, given the distribution coefficients β T i , the only constranits available are given by [START_REF] Csík | Space-time residual distribution schemes for hyperbolic conservation laws[END_REF]. These constraints actually correspond to the conservation requirement [START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF], and to the requirement that, whenever the differential operator

1 1 2 2 3 3 M M T 1 T 2 T 3 C 1 C 2 C 3
r(u h ) = ∂ t u h + a • ∇u h is (as in the steady case) locally constant, say r(u h ) T = r T , then Φ i = β T i |T | r T
These properties are always verified if we can find a Petrov-Galerkin test function ω i such that on every T

Φ i (u h ) = T ω i r(u h ) dx dy (18)
with

j∈T ω j = 1 , 1 |T | T ω i dx dy = β T i ( 19
)
As we will show later, the formal second order of accuracy is guaranteed as long as ω i is locally bounded. Note also that

m T ij = T ω i ϕ j dx dy (20)
The number of functions that verify these constraints is infinite. For example, to obtain the formulation F1 (cf. equation ( 13)) one can choose on each T

ω F1 i T = β T i χ T
having denoted by χ T the characteristic function

χ T (x, y) = 1 if (x, y) ∈ T 0 if (x, y) / ∈ T
Conversely, formulation F2 (cf. equation ( 14)) is obtained for

ω F2 i = ϕ i + T |i∈T δ ϕi χ T
Formulations F3 and F4 are instead obtained by taking for example (cf. equations ( 16) and ( 17))

ω F3/F4 i T = χ Ti
Moreover, for any given test function ω i verifying all the consistency, conservation, and accuracy constraints, we can easily come up with a modified function, say ω i with all the desirable properties. For example, if we can find three bounded functions, say f 1 , f 2 , and f 3 such that

3 j=1 f j = C f
with C f a constant, we can modify ω i as follows

ω i = ω i + K(f i -f i ) , f i = 1 |T | T f i dx dy
with K an arbitrary parameter ! Note that this term does not affect consistency or conservation, due to the fact that (using for the nodes of T the local renumbering {i, j, k} → {1, 2, 3})

3 j=1 f j -f j = 0 , T f j -f j dx dy = 0 INRIA
nor it does pollute the accuracy of the discretization, as long as the each f i is bounded. Moreover, in the P 1 case ∇u h is constant per element, so that

T f j -f j a • ∇u h dx dy = 0
so that the extra term only affects the form of the mass matrix. This leads clearly to quite a large number of consistent mass matrices, and extra constraints are needed to make sure one does the right thing. For example stability is not at all taken into account in the analysis performed so far.

The last observation leads to very interesting consequences if we take f i = ϕ i . In this case we have

ϕ i = 1 |T | T ϕ i = 1 3
So any mass matrix can be modified as

m T ij = m T ij + K T (ϕ i -ϕ i )ϕ j dx dy = T ω i ϕ j dx dy + K T (ϕ i -ϕ i )ϕ j dx dy
that leads to the semi-discrete scheme

T |i∈T   j∈T m T ij + Kδm ij du j dt + β i φ(u h )   = 0 , δm ij = |T | 36 (3δ ij -1) (21)
with δ ij Kroenecker's delta. As already noted in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF], the matrix δm ij is symmetric, and defines a dissipation operator, that is

v T [δm ij ] v ≥ 0 , ∀ v ∈ R 3
In the last reference, this term has been used to provide further stabilization to a nonlinear second-order variant of a Lax-Friedrich's scheme. The interesting observation is that if we take K = 3 and apply the modification to the Galerkin scheme we obtain :

m ij = Galerkin |T | 12 (δ ij + 1) + 3δmij |T | 12 (3δ ij -1) = |T | 3 δ ij
Which is just another way to show that mass lumping for the Galerkin scheme does not reduce the accuracy in the P 1 case but it does introduce a degree of dissipation.

It is also worth noting that comparing ( 13) and ( 14) with [START_REF] Cohen | High order triangular finite elements with mass lumping for the wave equation[END_REF], one immediately sees that

m F2 ij = m F1 ij + δm ij (22)
The third and first formulations are linked by a very similar relation :

m F3 ij = m F1 ij + δm ij , δm ij = |T | 3 β i δ ij -β i β T j ( 23
)
RR n 6998

Provided that β i ≥ 0 ∀ i, then the symmetric matrix δm ij also defines a dissipation operator. In particular, ∀ v ∈ R 3 we have

v T δm ij v = |T | 3 β T 1 β T 2 (v 1 -v 2 ) 2 + |T | 3 β T 1 β T 3 (v 1 -v 3 ) 2 + |T | 3 β T 3 β T 2 (v 3 -v 2 ) 2 ≥ 0
A similar relation holds for the last formulation, only this time we have

m F1 ij = m F4 ij + δm ij (24)
This last relation allows finally to show that all the formulations are equivalent up to a dissipation term :

m F1 ij = m F4 ij + δm ij m F2 ij = m F1 ij + δm ij = m F4 ij + δm ij + δm ij m F3 ij = m F1 ij + δm ij = m F4 ij + 2 δm ij (25)
Remark 3.1. The last relations show that, for a fixed approximation of the advection operator, the formulation F4 is the least dissipative of all.

Remark 3.2. In one space dimension, if the spatial discretization is given by the classical 1d upwind scheme, the formulations F1, F3, and F4 become identical. The formulation F2, instead, reduces to the 1D SUPG scheme obtained by defining the τ SUPG parameter as [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/leastsquares algorithms[END_REF][START_REF] Hauke | Fourier analysis of semi-discrete and spacetime stabilized methods for the advectivediffusivereactive equation: I[END_REF] τ = ∆ x 2|a|

The objective of the following sections is to show how to make use of the properties discussed above to contruct fully explicit residual based schemes.

Mass lumping and bubble stabilization

In this section we make use of the elements already discussed to propose a fully explicit variant of the residual distribution schemes described in the previous sections. The start with we make the assumption that, however complex the definition of the β i coefficients and of the mass matrix, there exists a uniformly bounded and locally differentiable function γ i , such that we can rewrite the discretization [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for navier stokes simulations[END_REF] as

Ω ϕ i (∂ t u h + a • ∇u h ) dx dy + T |i∈T T γ i (∂ t u h + a • ∇u h ) dx dy = 0 ( 26
)
where γ i plays the role of a "stabilizing" bubble function, and satisfies

                     j∈T γ j = 0 T (ϕ i + γ i )∂ t u h dx dy = j∈T m ij du j dt 1 |T | T (ϕ i + γ i )dx dy = β i (27)

INRIA

The last three relations guarantee the satisfaction of the conservation property, and of the consistency with the (given) spatial discretization, so that ultimately (cf. equation ( 9))

Φ i (u h ) = j∈T m ij du j dt + β i φ(u h ) = T ϕ i r(u h ) dx dy + T γ i r(u h ) dx dy (28)
The role of the bubble is of course that of providing stabilization to the otherwise unstable Galerkin scheme. Note that it is not necesary to actually show particular forms of such function which, as we shall see in the following, is just an artifact allowing to analyze the accuracy of the schemes proposed in the paper. Nevertheless, whenever we can exhibit the existence of a Petrov-Galerkin test function ω i such that (18) holds, we can simply set

γ i | T = ω i | T -ϕ i | T (29)
For example, for the formulations seen in the previous sections we have :

γ F1 i T = β i -ϕ i | T ; γ F2 i T = δ ϕi = β i - 1 3 ; γ F3/F4 i T = χ Ti -ϕ i | T . (30)

Construction of explicit schemes

The starting point of the contruction is to choose an explicit time-stepping scheme. We will focus here on Runge-Kutta (RK) schemes, however other possibilities exist, and will be studied in the future. Let us denote by δu k = u k -u n the increment of the k-th step of a given explicit RK scheme. Similarly, let f k be the k-th step evolution operator so that for the problem

du dt + f (u) = 0
we can rewrite each RK step as

δu k ∆t + f k = 0
In particular, in the following we will denote by r k the quantity

r k = δu k ∆t + f k (31)
For example for the classical TVD RK2 scheme we have

         r 1 = δu 1 ∆t + f 1 = 0 , f 1 = a • ∇u n r 2 = δu 2 ∆t + f 2 = 0 , f 2 = 1 2 a • ∇u n + 1 2 a • ∇u 1 (32)
RR n 6998

Similarly, the TVD RK3 scheme gives

                     r 1 = δu 1 ∆t + f 1 = 0 , f 1 = a • ∇u n r 2 = δu 2 ∆t + f 2 = 0 , f 2 = 1 4 a • ∇u n + 1 4 a • ∇u 1 r 3 = δu 3 ∆t + f 3 = 0 , f 3 = 1 6 a • ∇u n + 1 6 a • ∇u 1 + 2 3 a • ∇u 2 (33)
With this notation we can write the k-th step of the RK time integrator as

r k = 0
Its Galerkin discretization writes

Ω ϕ i r k (u h ) dx dy = Ω ϕ i δu k h ∆t dx dy + Ω ϕ i f k (u h ) dx dy = 0
The next step is to add the contribution of the bubble. The standard approach would be to write this contribution as

T |i∈T T γ i r k (u h ) dx dy = T |i∈T T γ i δu k h ∆t + f k (u h ) dx dy ,
however, even when lumping the Galerkin component of the mass matrix, this would lead to a scheme with a non-diagonal mass matrix, still requiring the solution of an a-priori nonlinear system at each RK step. What we propose is to replace only in the bubble contribution the k-th step residual r k (u h ) by a modified residual r k (u h ), which makes use of a different approximation of the time derivative. In practice, we will look for r k (u h )s differing from r k (u h ) only in the definition of the time increment, that is

r k (u h ) = δu k ∆t + f k (34)
The constraints on r k (u h ) guaranteeing that the overall accuracy of the discretization is not deteriorated will be given in the next section. For the moment

INRIA

we observe that when adding this contribution we obtain :

0 = Ω ϕ i r k (u h ) dx dy + T |i∈T T γ i r k (u h ) dx dy = Ω ϕ i δu k h ∆t dx dy + Ω ϕ i f k (u h ) dx dy + T |i∈T T γ i δu k h ∆t + f k (u h ) dx dy = Ω ϕ i δu k h ∆t dx dy + Ω ϕ i δu k h ∆t + f k (u h ) dx dy + T |i∈T T γ i δu k h ∆t + f k (u h ) dx dy - Ω ϕ i δu k h ∆t dx dy = Ω ϕ i δu k h ∆t + T |i∈T T (ϕ i + γ i ) δu k h ∆t + f k (u h ) dx dy - Ω ϕ i δu k h ∆t dx dy = Ω ϕ i δu k h ∆t + T |i∈T Φ RK(k) i - Ω ϕ i δu k h ∆t dx dy (35)
The last relations are obtained by first adding and substracting the Galerkin integral of the approximate time increment δu k /∆t, then using the properties of the bubble function γ i (cf. equations ( 9) and ( 30)), and finally introducing the fully discrete split residuals :

Φ RK(k) i = T (ϕ i + γ i )r k (u h ) dx dy = j∈T m ij δu k j ∆t + β i φ RK(k) (u h ) (36)
with

φ RK(k) = T f k (u h ) dx dy and j∈T Φ RK(k) j = T r h (u h ) dx dy = Φ RK(k)
(37) At this point two possibilities exist, leading to two distinct classes of methods.

Selectively Lumped (SL) schemes. If in the last line of [START_REF] Jund | Méthodes d'éléments finis d'ordre élevé pour la simulations numérique de la propagation d'ondes[END_REF] only massmatrix corresponding to the first Galerkin integral is lumped we obtain the following explicit formulation :

|S i | δu k i ∆t = - T |∈T   Φ RK(k) i - T ϕ i δu k h ∆t dx dy   (38)
In this case, the effect of the lumping simply leads to the following modification of the mass matrix. If by m G ij we denote the Galerkin mass matrix, we have for the selectively lumped schemes :

m SL ij = m T ij -m G ij , m G ij = |T | 12 (δ ij + 1) (39)
Globally Lumped (GL) schemes. If we lump all the Galerkin integrals we obtain the following explicit formulation :

|S i | δu k i -δu k i ∆t = - T |∈T Φ RK(k) i (40)
In this case, there is no modification at all on the residual distribution formulation, however, the lumping modifies the explicit iterations that now depend on the definition of δu k . All that remains to do is to properly define r k and δu k , such that we can still keep the desired accuracy.

Accuracy and time-stepping

We want to derive a sufficient condition on the r k guaranteeing that the accuracy of the Runge-Kutta Galerkin approximation is not lost when adding the bubble contribution. To do this we use a truncation error analysis, following the approach of [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF]. All the details of the analysis are given in two appendices at the end of the paper. The general idea of the proof is, given a classical solution w, and a smooth function ψ ∈ C 1 0 (Ω), to verify under which conditions the truncation error

E n = i∈T h ψ i Ω ϕ i δw n+1 h ∆t + f n+1 (w h ) dx dy+ i∈T h ψ i T |i∈T T γ i δw n+1 h ∆t + f n+1 (w h ) dx dy (41)
is of an order O(h p ). The analysis reported in appendix 1 and 2 is done for the general case p ≥ 2, even though the paper focuses only on the case p = 2. Note that in the definition of the error we have used the notation introduced in the previous section. This means that f n+1 (w h ) represents the discrete evolution operator of the last RK step, which atually makes use of flux values at known time-steps (cf. equations ( 32) and ( 33)). The analysis makes use of the following two hypotheses.

Hypothesis 4.1 (RK truncation error). Given a smooth classical solution w such that ∂ t w + ∇ • F (w) = 0, a p-th order RK scheme verifies the truncation error estimate

r n+1 (w) = δw n+1 ∆t + f n+1 (w) = C RK ∆t p Hypothesis 4.2 (Approximate semi-discrete residual estimate). Given a smooth classical solution w such that ∂ t w + ∇ • F (w) = 0, the approximate semidiscrete residual r verifies the estimate r n+1 (w) = δw n+1 ∆t + f n+1 (w) = C RK ∆t l for some l ≤ p.
The main result is summarized by the following proposition.

Proposition 4.3 (Accuracy and time-stepping). Given a p-th order spatial approximation and a p-th order RK scheme verifying hypothesis 4.1, the truncation error (41) verifies an estimate of the type

E n ≤ C h p provided that INRIA 1.
the bubble γ i is uniformly bounded 2. the approximate semi-discrete residual verifies hypothesis 4.2 with

l ≥ p -1
In particular, in the second order case of interest here, it is enough to provide definitions of the approximate time increments yielding a first order semidiscrete operator.

Remark 4.4 (Accuracy, time-stepping, and distribution coefficients). As seen in section 3, for all the known consistent formulations of RD we can provide define the bubble fuction as γ i = ω i -ψ i is always bounded. For the formulations recalled in section 3, ω i , and hence γ i , is bounded whenever the distribution coefficients β T i are. To end the construction we give particular definitions of r k that satisfy hypothesis 4.2 (see appendix 2) :

RK2 scheme δu 1 = 0 ⇒ r 1 = ∇ • F (u n ) δu 2 = δu n+1 = u 1 -u n ⇒ r 2 = u 1 -u n ∆t + ∇ • F (u n ) + ∇ • F (u 1 ) 2 (42)
When combining this definition with the updates [START_REF] Woodward | The numerical simulation of twodimensional flows with strong shocks[END_REF], and (40), we obtain for the SL schemes

             |S i | u 1 i -u n i ∆t = - T |i∈T β i φ(u n h ) |S i | u n+1 i -u n i ∆t = - T |i∈T Φ RK2(2) i - j∈T m G ij u 1 j -u n j ∆t (43)
with m G ij as in (39), and with

Φ RK2(2) i = j∈T m ij u 1 j -u n j ∆t + 1 2 β i φ(u n h ) + φ(u 1 h )
The update for the GL schemes is somewhat simpler and given by

             |S i | u 1 i -u n i ∆t = - T |i∈T β i φ(t n ) |S i | u n+1 i -u 1 i ∆t = - T |i∈T Φ RK2(2) i (44)
RK3 scheme

δu 1 = 0 ⇒ r 1 = ∇ • F (u n ) δu 2 = u 1 -u n 2 ⇒ r 2 = u 1 -u n 2∆t + ∇ • F (u n ) + ∇ • F (u 1 ) 2 δu 3 = δu n+1 = 2(u 2 -u n ) ⇒ r 3 = 2(u 2 -u n ) ∆t + ∇ • F (u n ) + ∇ • F (u 1 ) + 4∇ • F (u 2 ) 6 (45)
Note that in this case the coefficients involved in the definition take into account the fact that u 1 and u 2 are initial guesses for the solution at times t n + ∆t and t n + ∆t/2, respectively. When combining this definition with the updates [START_REF] Woodward | The numerical simulation of twodimensional flows with strong shocks[END_REF], and (40), we obtain for the selectively lumped schemes

                           |S i | u 1 i -u n i ∆t = - T |i∈T β i φ(u n h ) |S i | u 2 i -u n i ∆t = - T |i∈T Φ RK3(2) i - j∈T m G ij u 1 j -u n j 2∆t |S i | u n+1 i -u n i ∆t = - T |i∈T Φ RK3(3) i - j∈T m G ij 2 u 2 j -u n j ∆t (46)
with m G ij as in (39), and with

Φ RK3(2) i = j∈T m ij u 1 j -u n j 2∆t + 1 4 β i φ(u n h ) + φ(u 1 h ) and Φ RK3(3) i = j∈T m ij 2 u 2 j -u n j ∆t + β i 1 6 φ(u n h ) + 1 6 φ(u 1 h ) + 2 3 φ(u 2 h )
As before, the update for the globally lumped schemes is somewhat simpler and given by

                         |S i | u 1 i -u n i ∆t = - T |i∈T β i φ(u n h ) |S i | ∆t u 2 i - u 1 i + u n i 2 = - T |i∈T Φ RK3(2) i 2|S i | ∆t u n+1 i + u n i 2 -u 2 i = - T |i∈T Φ RK3(3) i (47)
Remark 4.5 (Fluctuations/signals). Both formulations, the one based on selective lumping and the one based on global lumping, allow to see the RD component of the discretization as an error between two different approximations of the unknown at certain time levels. When using the formulation F1 of the RD discretization (cf. section 3, equation (13)), the second step of the RK2 scheme with selective lumping can be recast as

|S i | u n+1 i -u n i ∆t - Ω ϕ i u 1 h -u n h ∆t = - T |i∈T β i Φ RK2(2) ( 48
)
where

Φ RK2(2) = T u 1 h -u n h ∆t + 1 2 ∇ • F h (u n h ) + 1 2 ∇ • F h (u 1 h ) dx dy

INRIA

Clearly equation (48) expresses the error between two local approximations of the time variation of the unknown as a function of signals proportional to elemental errors represented by the residual Φ RK2 (2) . This is even more apparent in the case of the globally lumped scheme which reads, in RK2 case :

|S i | u n+1 i -u 1 i ∆t = - T |i∈T β i Φ RK2(2) (49)
The RK3 version of the last equation is obtained immediately from equation (47). In this case the RD wheighted average on the left expresses the between the two different approximations of the unknown at time t n+1 . The same remarks applies of course to the case of the RK3 schemes. In some way the explicit formulations proposed here lead us back to the original ideas of P.L.Roe [START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF] in which the nodal error is proportional to the signals sent by surrounding elements.

Remark 4.6 (Relations with explicit predictor-corrector). The explicit formulation proposed here is also related to the explicit predictor/multi-corrector formulation of the SUPG scheme used for example in [START_REF] Hauke | Fourier analysis of semi-discrete and spacetime stabilized methods for the advectivediffusivereactive equation: I[END_REF][START_REF] Hauke | Fourier analysis of semi-discrete and spacetime stabilized methods for the advectivediffusivereactive equation: II[END_REF][START_REF] Hauke | Fourier analysis of semi-discrete and spacetime stabilized methods for the advectivediffusivereactive equation: III[END_REF] (see also [START_REF] Tezduyar | Development of time-accurate finite element techniques for first order hyperbolic systems with emphasys on the compressible euler equations[END_REF][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] Scovazzi | Multi-scale lagrangian shock hydrodynamics on q1/p0 finite elements: Theoretical framework and twodimensional computations[END_REF]). In the simplest setting, in this formulation on replaces an implicit time integrator by a finite number of explicit steps. In the case of the Crank-Nicholson time integrator for example the idea is to rewrite the SUPG scheme as

|S i | u 1 i -u n i ∆t = - Ω ϕ i a • ∇u n h dx dy + T |i∈T T a • ∇ϕ i τ a • ∇u n h dx dy |S i | u k i -u n i ∆t = - Ω ϕ i a • ∇ u k-1 h + u n h 2 dx dy + T |i∈T T a • ∇ϕ i τ u k-1 h -u n h ∆t + a • ∇ u k-1 h + u n h 2 dx dy
where k ≥ 2 and u n+1 i = u kmax i . The second relation can immediately be recast as

|S i | u k i -u n i ∆t = - T |i∈T   T |i∈T Φ SUPG(k) i - T ϕ i u k-1 h -u n h ∆t dx dy   with Φ SUPG(k) i = T (ϕ i + a • ∇ϕ i τ ) u k-1 h -u n h ∆t + a • ∇ u k-1 h + u n h 2 dx dy and j∈T Φ SUPG(k) j = Φ k = T u k-1 h -u n h ∆t + a • ∇ u k-1 h + u n h 2 dx dy
This is basically the selectively lumped formulation of the SUPG scheme in RD form. In particular, when using only one correction step we end up exactly with the RK2 scheme (43).

Remark 4.7 (Explicit Residual Based Finite Volume formulation). The approach presented here finds application also in the case of Finite Volume discretizations where the stabilization operator is proportional to some local approximation of the residual, rather than to local variations of the solution. Such schemes have been proposed for example in [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF][START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF] and, in a different spirit, in [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF]. The schemes of [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF][START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF] in their basic formulation can be rewritten as

|C i | du i dt + ∂Ci H C • n dl - 1 2 j h j Ψ ij Φ ij = 0 (50)
where H C is a centered finite volume numerical flux, while the last term represent stabilization terms. These terms are function of a local residual Φ ij computed on the staggered cell C ij (cf. figure 2) and defined as [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF][START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF]]

Φ ij = Cij du h dt + ∇ • F h (u h ) dx dy
where now u h is a polynomial approximation of the unknown recontructed starting from cell averages. We refer to [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF][START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF] for further details, and in particular for the definition of the local mesh size h j , and of the Ψ ij parameter in (50).

The important point is that the residual Φ ij has to include the time derivative of the numerical unknown to attain consistency in the spatial discretization. In [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF][START_REF] Corre | A residual-based scheme for computing compressible flows on unstructured grids[END_REF] the authors use the same discrete operator to approximate both the time derivative of u i in (50), and in Φ ij . this naturally leads to the appearance of a mass matrix rendering the scheme implicit in space. The approach proposed here allows to overcome this limitation allowing the contruction of an explicit RK schemes in which the time derivative in Φ ij is approximated by time increments using known values of the discrete solution. In the RK2 case for example the scheme would read : This section is devoted to the description of the schemes actually used in the numerical tests discussed in the following paragraphs, and of some details relative to their implementation. We will discuss the results obtained with four well known schemes : the LDA scheme, the blended LDAN scheme or B scheme, the Streamline Upwind scheme, or SU scheme for short, and a centered Blended scheme, Bc for short, constructed starting from the limited stabilized Laxfriedrich's scheme of [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF]. The results obtained with other RD schemes are very similar in nature. An importat remark is that so far we still have not worked on the adaptation of nonlinear RD discretizations to the construction proposed in tha paper. This means that we limited ourselves to code the schemes as they are presented in the literature. Improvements will be made in the future concerning strict preservation of positivity. Even so, as we will see, the numerical results are excellent, and confirm our theoretical analysis.

|C i | u 1 i -u n i ∆t + ∂Ci H n C • n dl - 1 2 j h j Ψ ij Cij ∇ • F n h dx dy = 0 |C i | u n+1 i -u n i ∆t + ∂Ci H n C + H 1 C 2 • n dl - 1 2 j h j Ψ ij Cij u 1 h -u n h ∆t + ∇ • F n h + ∇ • F 1 h 2 dx dy = 0 with H 1 C = H C (u 1 h ) and F 1 h = F h (u 1 h ). C i C j C ij ∂C i

LDA scheme

We test our construction on the well known second order linear multidimensional upwind LDA scheme, defined by the distribution coefficients [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF] :

β LDA i = k + i   j∈T k + j   -1 (51)
where, using the notation of equation ( 2), we define ∀T ∈ T h

k i = 1 2 ∂F(u) ∂u • n i (52)
with u the arithmetic average of the values of u h in the nodes of T . Note that in the case of a system of conservation laws, the k i s are matrices, and their sign in (51) is computed in the standard matrix sense, via eigenvalue decomposition.

For more details on the definition and properties of the LDA scheme, the reader is referred to [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF][START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF].

In the scalar case we will compare the results obtained when using the different formulations recalled in section 3. In particular, the scheme has been coded exactly as decribed in equations ( 43) and ( 46), for the selectively lumped scheme, and in equations ( 44) and (47) for the global lumped scheme. In both cases, we replace the quantities Φ RK2(k) i and Φ RK3(k) i by (see equations ( 37), (42), and (45) for the notation)

Φ LDA(k) i = j∈T m LDA ij δu k j + ∆t β LDA i φ RK(k) (53)
In particular, the form of the mass matrix m LDA ij will depend of the formulation chosen (cf. section 3). To shorten the text we will lump together the acronyms when referring to a scheme. For example, we shall speak of the LDA-F1-SL-RK2 when referring to the scheme obtained using the LDA distribution coefficients, the mass matrix of the formulation 1, selective lumping, and the RK2 scheme in time. Similarly for all the other combinations.

Blended LDA-N scheme

As suggested by its name, the Blended LDA-N scheme, or B scheme for short, is a blending between the LDA scheme of section 5.1 with the first order positive multidimensional upwind N scheme defined by the spatial splitting [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF][START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF][START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF]

φ N i = k + i (u i -u in ) , u in =   j∈T k + j   -1   -φ(u h ) + j∈T k + j u j  
In particular, following [START_REF] Abgrall | Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF], we set for the B scheme

Φ B(k) i = (1 -l(u h ))Φ LDA(k) i + l(u h )Φ N(k) (54) with Φ LDA(k) i
given by ( 53) and with

Φ N(k) i = |T | 3 δu k i + ∆t φ N(k) i having denoted by φ N(k) i
the spatial contribution of the N scheme corresponding to the k-the RK step. Expression (54) is used in (43) , ( 46), ( 44) and (47) to replace Φ

RK2(k) i and Φ RK3(k) i .
Concerning the blending paramenter l(u h ) we have used the standard definition of Deconinck et al. [START_REF] Deconinck | Status of multidimensional upwind residual distribution schemes and applications in aeronautics[END_REF][START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF] (cf. also equation ( 37)) :

l(u h ) = |Φ RK(k) | j∈T |Φ N(k) j |
for systems of equations, the blending procedure has been performed on residuals projected in characteristic directions, as explained in [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for steady flow problems[END_REF][START_REF] Abgrall | Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF].

As a last remark, we note that only when using global lumping for l(u h ) = 1 does the B scheme defined by (54) reduce to the N scheme with RK time integration. In the selective lumping case, for l(u h ) = 1 we get (cf. equation [START_REF] Woodward | The numerical simulation of twodimensional flows with strong shocks[END_REF])

|S i | δu k i ∆t + T |i∈T φ N(k) i = |S i | δu k i ∆t - Ω ϕ i δu k h ∆t dx dy
where the left hand side corresponds to the k-th RK step of the N scheme, while the right hand side contains some kind of anti-diffusive correction (cf. section 3.).

The SU scheme

To test the behavior of our formulation with different type of discretizations, we also consider centered schemes. The first is referred to in the RD literature either as SUPG scheme or as LW scheme. It is defined by the distribution coefficients

β SU i = 1 3 + k i τ (55)

INRIA

Independently on the definition of the scaling parameter τ , the second term on the last definition introduces some Streamline Upwinding in the distribution [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF], which is why we refer to this scheme as to the SU scheme. In our computations we have taken

τ =   j∈T |k j |   -1 (56)
For the Euler equations, last expression is meant in the usual matrix sense. Finally, we replace the quantities Φ

RK2(k) i and Φ RK3(k) i
by (see equations ( 37), (42), and (45) for the notation)

Φ SU(k) i = j∈T m SU ij δu k j + ∆t β SU i φ RK(k) (57)
As for the LDA scheme, also for the SU scheme the form of the mass matrix m SU ij depends of the formulation chosen (cf. section 3).

Central blended scheme

In [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] the authors introduce a centered discretization based on a nonlinear variant of a Lax-Friedrich's scheme. This limited stabilized Lax-Friedrich's scheme, or LLFs scheme, is obtained starting from the positive first order Lax-Friedrich's (LF) splitting

φ LF i = 1 3   φ(u h ) + α LF j∈T (u i -u j )   ( 58
)
where α LF is the Lax-Friedrich's dissipation coefficient which we set to

α LF = 1 2 a T h T , a T = max j]∈T
∂F (u j ) ∂u in the scalar case, while for the Euler equations we have set

α LF = 1 2 max j]∈T ( u j + a j ) h T
with u the flow speed, a the speed of sound, and h T a reference length for element T .

The LF scheme is only first order. To obtain a formally second order nonlinear splitting we proceed as follows. First we define the LF-RK splitting

Φ LF(k) i = |T | 3 δu k i + ∆t φ LF(k) i having denoted by φ LF(k) i
the k-the RK step of the spatial operator (58). Next, we compute bounded distribution coefficients by applying a sign preserving nonlinear mapping. Several ways of doing this exist, and we refer to [START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF][START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF] for a discussion. Here, we set (cf. equation (37)) :

β LLF i = max 0 , Φ LF(k) i Φ RK(k) j∈T max 0 , Φ LF(k) j Φ RK(k) (59) RR n 6998
The limited LF scheme is then defined by

Φ LLF(k) i = β LLF i Φ RK(k)
As shown in previous work [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | Stable and convergent residual distribution for time-dependent conservation laws[END_REF], the limiter (60) not taking into account the directional propagation of the information typical of hyperbolic problems, the LLF scheme shows mild spurious modes that eventually reduce its accuracy to first order. This is cured as in the above references by adding an upwind bias inspired by the SU scheme :

β LLFs i = β LLF i + δ(u h )k i τ (60)
with τ as in (56). We refer the reader to [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | Stable and convergent residual distribution for time-dependent conservation laws[END_REF] for more details on the theoretical background leading to this choice. We limit ouselves to recall that δ(u h ) is a smoothness sensor such that δ(u h ) = 1 in smooth areas, while δ = O(h T ) in presence of discontinuities. In our computations we have set in the scalar case [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | Stable and convergent residual distribution for time-dependent conservation laws[END_REF]]

δ(u h ) = min 1, ∆t h 2 T a T |u| T |Φ RK(k) | (61)
where |u| T is the maximum of the absolute value of the solution over the element.

For the Euler equations, the extension is done following [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF] : the limiter (60) is evaluated on residual projected on local characteristic directions, while the |Φ RK(k) | in (61) is replaced by the scalar entropy component of Φ RK(k) . This is computed as

ϕ s = l 0 • Φ RK(k)
where l 0 is the left eigenvector of the flux Jacobian corresponding to the entropy wave. For the Euler equations δ(u h ) is then the scalar quantity (see [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF] for more)

δ(u h ) = min 1, ∆t h 2 T |ϕ s | (62)
Normally, we would set

Φ LLFs(k) i = β LLFs i Φ RK(k) (63)
and replace Φ 46), ( 44) and ( 47) by ( 63). However, we found that much better results are obtained, at negligible extra cost, by using the central blended scheme, or Bc scheme for short, defined by

RK2(k) i and Φ RK3(k) i in (43) , (
Φ Bc(k) i = β Bc i Φ RK(k) , β Bc i = δ(u h )β SU i + (1 -δ(u h ))β LLF i ( 64
)
From definitions (55) and (60) we immediately see that

Φ Bc(k) i = Φ LLFs(k) i + δ(u h ) 1 3 Φ RK(k)
Compared to the LLFs scheme, the only extra cost to evaluate (64) is the addition of the term δ(u h )Φ RK(k) /3. Being δ(u h ) always a scalar, this cost is clearly negligible.

INRIA

Computation of the time step

All the numerical results presented in the folllowing section have been obtained by computing the timestep as (cf. section 5.3) :

∆t = min i∈T h |S i | T |i∈T α LF (65)
For all the nonlinear problems considered, α LF is evaluated using solution values at the last known time step.

A fourier analysis on structured triangulations is under way to have a better estimate of the time step stability limit for the linear schemes.

Scalar results

The scalar tests we present have two objectives : verify the accuracy of our explicit formulation for different forms of the mass matrix, and for schemes of different nature (multidimensional upwind, and centered) ; test the non-oscillatory nature of the results obtained with the nonlinear schemes, when no modifications are introduced to take into account the additional terms introduced by RK formulation.

Unless stated, all the numerical tests, including the Euler tests, have been performed on unstructured triangulations with the topology shown on figure 3. We use this test to study the accuracy of the different schemes discussed in the paper. The accuracy is monitored by the convergence of the L 1 norm of the error with respect to the exact solution. The behaviour of the L ∞ and L 2 norms is qualitatively and quantitatively very similar.

The first exercice is to verify that indeed our RK formulation leads to second order discretizations, independently on the starting form of the (consistent) mass matrix. We perform the test for all the mass matrix formulations for the LDA scheme, which is the most popular multidimensional upwind RD scheme.

The results are summarized in figures 4 and 5, where we report the grid convergence history and the rate of convergence history, respectively. The first remark we can make is that our explicit formulation does lead to a second order discretization. This is clear especially from the rates of convergence observed. What is more interesting is that the RK2 schemes all yield the same accuracy, while the RK3 scheme with global lumping seem to actually be less and less accurate as the mesh is refined. We believe this might be the consequence of a (mild) linear stability problem. We are currently performing a Fourier analysis on structured grids to better understand this behaviour. There are minor differences between the different mass matrix forms which, in our opinion, do not justify the use of the more complex formulations F3 and F4 (cf. section 3), especially in view of the extension to systems.

We repeat the same exercise with the SU scheme, only this time we only test the mass matrix formulations F1 and F2 (cf. section 3). The results are shown on figure 6. The same remarks made for the LDA scheme apply also to the SU distribution : second order of accuracy is obtained already with the RK2 scheme, independently of the mass matrix and lumping choices ; the RK3 scheme with global lumping suffers from a drop in the convergence rate, which might be caused by the presence of a linear instability.

We now come to the nonlinear schemes. We firt test the B scheme, using either formulation F1, or formulation F2 for the LDA mass matrix. The results are displayed on figure 7. The asymptotic rate of convergence obtained is about 1.75-1.8, independently on the formulation. Clearly, when using global lumping, the drop in convergence speed of the LDA affacts the B scheme as well. Lastly, on figure 8 we report the results obtained with the Bc scheme. Once more, we observe asymptotic convergence rates ranging from 1.7 to 1.9, with the exception of the RK3 scheme in conjunction with global lumping.

We believe these tests confirm our theoretical contruction. In particular the fact that with the RK2 scheme one already obtains a second order discretization. Moreover, the fact that different forms of the mass matrix lead to very similar accuracy properties leads us to the conclusion that the choice of the form of the mass matrix should be done on the basis of stability (or positivity eventually) considerations. This is the objective of our current investigations.

INRIA -2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -4 -3.5 -3 -2.5 -2 -1.5 log 10 (h) 1 1 2 1.6 log 10 ( ǫ L 1 ) F1-SL-RK2 F1-GL-RK2 F1-SL-RK3 F1-GL-RK3 -2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -4 -3.5 -3 -2.5 -2 -1.5 log 10 (h) 1 1 2 1.6 log 10 ( ǫ L 1 ) F2-SL-RK2 F2-GL-RK2 F2-SL-RK3 F2-GL-RK3 -2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -4 -3.5

Discontinuous solutions : 2d Burger's equation

We consider now the nonlinear 2d Burger's equation

∂ t u + ∂ x u 2 2 + ∂ y u 2 2 = 0
We solve the problem on the square [-1, 1] 2 with the discontinuous initial solution

u 0 = 1 if x ∈ [-0.6, -0.1] × [-0.35, 0.15] 0 otherwise
The problem is solved up to the final time t = 1 on an unstructured triangulation with the topology shown on figure 3, and reference size h ≈ 1/80. We compare the results of all the schemes considered. Only the simplest mass matrix form F1 (cf. section 3) has been used.

We first consider the multidimensional upwind LDA and B schemes. The results for different RK schemes and lumping strategy are shown on figures from 9 to 16. Concerning the LDA scheme, as one would expect, the solution exhibits oscillations near the discontinuities. These oscillations are much more pronounced when using selective lumping. More interesting are however the results of the B scheme, shown on figures 13 to 16. From all the contour plots we can see that the solution is smoother (the kinks close to the shock are less pronounced) when compared to the LDA scheme. When using selective lumping oscillations still appear close to the discontinuity. This, as observed at the end of section 5.2, is a consequence of the non-positive coefficients introduced by the Galerkin integral present when lumping selectively. A mixed formulation, in which these terms are also multiplied by the blending coefficient, might be used to cure the problem, but this is beyond the scopes of this paper and left for future work. The results obtained with global lumping show the expected monotone resolution of the discontinuities. Probably, the small negative undershoots can also be avoided by properly redefining the blending. Again, this is beyon the scopes of this paper and left for the future.

The results obtained with the centered distributions are displayed on figures from 17 to 24. The qualitative bahavior of these schemes is similar to the one of the multidimensional upwind discretizations. The linear SU scheme gives oscillations near the discontinuities. Milder oscillations are obtained when using global lumping. Concerning the Bc scheme, the results show smoother contours (the kinks close to the shock are less pronounced), and less oscillations. However, only the results obtained using global lumping show a monotone resolution of the moving shock.

We judge the results obtained on this nonlinear problem very encouraging : even without modifying the basic RD distribution, the nonlinear second order explicit RK-RD schemes can yield monotone solutions. This is further confirmed by the Euler results discussed hereafter.

Euler equations

As already said in section 5, the extension of the schemes to the system of Euler equations is performed formally. As in the scalar case, the nonlinear schemes are not modified to take into account the additional terms coming from the explicit RK formulation and to improve their behavior close to discontinuities. For simplicity, only the schemes based on the simple mass matrix formulation F1 (cf. section 3) are tested. The objective of the tests is to asses the accuracy of the discretizations, and the behavior of the nonlinear schemes in presence of a strong moving planar shock, and for more complex flow structures involving several contact lines and interactions between shocks and expansions.

Advection of a vortex : grid convergence

The accuracy of the schemes is measured on the advection of a constant density vortex. The test has been initially proposed in [START_REF] Dobes | Second order blended multidimensional upwind residual distribution scheme for steady and unsteady computations[END_REF], to which we refer for all the details concerning its implementation. The solution involves the advection of a vortex with a constant density profile, and a smooth pressure variation of which the analytical form is known [START_REF] Dobes | Second order blended multidimensional upwind residual distribution scheme for steady and unsteady computations[END_REF]. We solve the problem on a set of 5 unstructured grids with the topology shown on figure 3. The coarsest grid as a reference size of h ≈ 1/20. The other meshes are obtained by means of 4 successive steps of conformal refinement. We measure the accuracy by means of the L 2 norm of the relative pressure error

ǫ p = p -p exact p ∞
see [START_REF] Dobes | Second order blended multidimensional upwind residual distribution scheme for steady and unsteady computations[END_REF] for the definition of p ∞ and of -p exact . The behaviour of the L ∞ and L 1 norms is qualitatively and quantitatively very similar. The results are displayed on figures from 25 to 28 in terms of error convergence history, and rate of convergence history. The results are qualitatively very similar to the ones discussed in section 6.1. With the exception of the first refining step, we do obtain roughly second order of convergence with all the schemes except the RK3 ones when using global lumping. These schemes, exactly as in the scalar case, show a more or less evident decrease in accuracy, as the mesh is refined.

The poor convergence rate at the first refinement step might be explained by the coarseness of the first meshes : the starting mesh only has 10 points through the vortex core, the second one 20 points. The drop in convergence rates of the RK3-GL schemes might be a consequence of a linear stability problem. This is under investigation. The main difference between the distribution strategies is that the B scheme gives a slighly smaller asymptotic accuracy of about 1.7, while all the others attain convergence rates closer to 2. While this is expected for the linear schemes, we believe the face that the Bc scheme shows a better convergence is due to the definition of the entropy smoothness sensor δ(u h ) proposed in [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF] and used for the blending, which really is turned on only very close to discontinuities. Once more, the improvement of the nonlinear schemes is a topic for future work. Nevertheless, the results obtained confirm once more our theoretical analysis. Left column : convergence history. Right column : convergence rates.

-2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 - 1.2 -5 -4.5 -4 -3.5 -3 -2.5 -2 log 10 (h) 1 1 2 1.6 SL
-RK2 GL-RK2 SL-RK3 GL-RK3 log 10 ( ǫp L 2 ) -2.4 -2.2 -2 -1.8 -1.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2
-2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 - 1.2 -5 -4.5 -4 -3.5 -3 -2.5
-2

log 10 (h) 1 1 2 1.6 SL-RK2 GL-RK2 SL-RK3 GL-RK3 log 10 ( ǫp L 2 ) -2.4 -2.2 -2 -1.8 -1.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2
INRIA -2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 - 1.2 -5 -4.5 -4 -3.5 -3 -2.5
-2

log 10 (h) 1 1 2 1.6 SL-RK2 GL-RK2 SL-RK3 GL-RK3 log 10 ( ǫp L 2 ) -2.4 -2.2 -2 -1.8 -1.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2
-RK2 GL-RK2 SL-RK3 GL-RK3 log 10 ( ǫp L 2 ) -2.4 -2.2 -2 -1.8 -1.6 0.8 1.2 -5 -4.5 -4 -3.5 -3 -2.5 -2 log 10 (h) 1 1 2 1.6 SL

Double Mach reflection

In this section we check the behavior of the nonlinear schemes in presence of a strong moving planar shock. The test case is that of a reflection of a Mach 10 oblique shock over a ramp proposed by Woodward and Colella in [START_REF] Woodward | The numerical simulation of twodimensional flows with strong shocks[END_REF], to which we refer for details concerning the implementation. The computations have been run on an untructured triangulation with the topology shown on figure 3 and reference mesh size h ≈ 1/100.

We display on figures 29 and 30 the density contours obtained with the B and Bc scheme, respectively. The first remark we can make is that even in presence of a strong moving shock both nonlinear discretizations yield quite smooth and non oscillatory results. To confirm this we report on figures 31 and 32 the density and entropy distributions across the shock and on the wall, respectively. From figure 31 we see that the shock is resolved very sharply. Only a small overshoot in its vicinity is observed in almost all the solutions. An exception to this is the B scheme with RK3 time stepping and global lumping. In this case, as it can be also seen from the contour plot on figure 29 (bottomright), we obtain some strange behavior in correspondence of the compression region where the bent incoming shock changes curvature. This seems to be a feature propagating from the upper boundary, where the exact shock movement is strongly imposed. This affects both the shock profile, as seen on figure 31, and the structures on the lower wall, as seen on figure 32. So far we have not been able to explain this behavior.

Apart from the above remarks, the solutions obtained are very satisfactory. The minumim and maximum values of the density, reported on table 1 also show that the minumim of the density is always very close to 1.4 (its analytical value), again with the exception of the B-GL-RK3 scheme.

Mach 3 wind tunnel with a step

This final test is also taken from [START_REF] Woodward | The numerical simulation of twodimensional flows with strong shocks[END_REF] and involves the formation and evolution of a moving shock in a Mach 3 wind tunnel with a step. We refer to [START_REF] Woodward | The numerical simulation of twodimensional flows with strong shocks[END_REF] for details concerning the implementation of the test case. The mesh used for the computations is the same used in [START_REF] Csík | Space time residual distribution schemes for hyperbolic conservation laws on unstructured linear finite elements[END_REF][START_REF] Csík | Space-time residual distribution schemes for hyperbolic conservation laws[END_REF]. A close up view in vicinity of the corner of the step is displayed on figure 33. The reference mesh size far from Left : B scheme. Right : Bc scheme the corner is h ≈ 1/80. The mesh is refined at the corner to attain a mimium size of h ≈ 1/1000. No particular numerical treatment has been used near the corner to handle the supersonic expansion taking place during the transient.

The solutions obtained at times t = 0.5, t = 1.5, and t = 4.0 with the B and Bc schemes are shown on figures from 34 to 41. All the figures show a monotone and sharp resolution of the shocks, and of the contact lines. The nonoscillatory character of the results is confirmed by the line plots of the solution along the upper wall of the step (line y = 0.2 containing the corner singularity). We never obtained negative densities. Note that this is a test case where the explicit formulation does give an advantage with respect to the implicit schemes based on Crank-Nicholson time integration [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]. Even if implicit in time, the positivity of the schemes proposed in the last references is still guaranteed by an explicit type time step restriction which, in presence of mesh refinement, renders the implicit formulation extremely time consuming. As a last remark, we note that the Bc scheme with selective lumping yields a much better resolution of the flow, as seen for example from the kinks of the initial shock (top-left on figures 38 and 39), and from the resolution of the contact emanating from the interaction of the corner expansion with the reflected shock (middle-left on figures 38 and 39).

Conclusions

In this paper we have provided a construction of explicit second order Residual Distribution schemes based on Runge-Kutta time integration. We used second order mass lumping and a finite element interpretation to achieve a discretization where the consistent RD mass matrix does not multiply the solution at the new time level thus allowing for a truly explicit solution procedure. All the theoretical arguments justifying our contruction have been thouroughly exposed, and strong numerical evidence has been given to confirm them.

The results obtained are very encouraging both concerning accuracy, and monotonicity. We think this work paves the way for a different class of RD schemes based on explicit, or mixed, time integration where the RD mass ma- trix does not necessarily need to be inverted.

Concerning the developments of the work reported in this paper, we mention the following points :

We are currently performing a Fourier analysis on structured triangulations to better undertand the linear stability properties of the linear schemes, when using different forms of hte mass matrix, and also to understand the influence of the type of lumping on the stability of the resulting scheme ;

Even though we judge the results presented here quite satisfactory (especially for the Euler equations) there is definitely space to improve the INRIA nonlinear schemes, taking into account the fully discrete RK time stepping and the terms arising from mass lumping in the positivity analysis ; An immediate application of the explicit RK-RD schemes is given by te Shallow Water equations where the preservation of the positivity of the depth leads, for the standard implicit RD, to a strict contraint on the time step [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. We think the RK-RD approach proposed would represent an improvement, still preserving most of the nice properties of the RD discretization ;

The extension to more than second order should be relatively stright forward when making use of higher order elements allowing higher order mass lumping. From this point of view we will profit of the work that has been done on the wave equation (see e.g. [START_REF] Cohen | High order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Jund | Arbitrary high order finite element schemes and high order mass lumping[END_REF] and references therein). that the modified semi-discrete operator used in the bubble function verifies the estimate

r n+1 (w) = δw n+1 ∆t + f n+1 (w) = C RK (w n )∆t p
where certainly l ≤ p. Both hypotheses are verified in appendix 2 for the RK2 and RK3 schemes considered in the paper.

Next we define, for the stabilized Galerkin scheme, the following truncation error :

E n = i∈T h ψ i Ω ϕ i r n+1 (w h ) dx dy + i∈T h ψ i T |i∈T T γ i r n+1 (w h ) dx dy (66)
where the ψ i s are nodal values of a C 1 0 (Ω) function, which is assumed to verify [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF]

ψ L ∞ (Ω) ≤ C ψ , ψ h L ∞ (Ω) ≤ C ψ h ; ∇ψ L ∞ (Ω) ≤ C ∇ψ , ∇ψ h L ∞ (Ω) ≤ C ∇ψ h (
67) having denoted by ψ h the p-th order polynomial approximation of ψ corresponding to the approximation space chosen. Similarly, w h represents the spatial interpolant of the given smooth exact solution w. As a second step, we can immediately rewite the error as

E n = I Ω ψ h r n+1 (w h ) dx dy + II T ∈T h j∈T T γ j ψ j r n+1 (w h ) dx dy
Next we estimate the two terms in the error. For I we immediately make use of hypothesis 4.1 :

I = Ω ψ h (∂ t w n h + ∇ • F h (w n h )) dx dy + Ω ψ h C RK (w n h )∆t p dx dy
Now, being w a smooth exact solution, we have ∂ t w n + ∇ • F (u n ) = 0, hence I can be rewritten as

I = Ω ψ h ∂ t (w n h -w n )dx dy + Ω ψ h ∇ • (F h (w n h ) -F (w n )) dx dy + Ω ψ h C RK (w n h )∆t p dx dy = Ω ψ h ∂ t (w n h -w n )dx dy - Ω (F h (w n h ) -F (w n )) • ∇ψ h dx dy + Ω ψ h C RK (w n h )∆t p dx dy
where, following [START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF], we have broken the second integral over elements, integrated by parts over each element, re-assembled, and used the fact that ψ ∈ C 1 0 (Ω). Due to the assumptions on ψ (cf. equation ( 67)), we can now use the properties of the approximation to estimate all terms, ending up with |I| ≤ C w (T h , w n )h p +C F (T h , w n )h p +C RK (T h , w n )∆t p = C 0 (T h , w n )h p +C RK (T h , w n )∆t p (68) This term is nothing else than the truncation error of the Galerkin scheme. As expected, it is of an order dictated purely by the spatial and temporal approximations.

INRIA

We now estimate the term II. First of all we note that since j γ j = 0, then we can write II = 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i ) r n+1 (w h ) dx dy where we racall that K denotes the number of DoF in an element. Next we use hypothesis 4.2 to get II = 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i )∂ t w n h dx dy

+ 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i )∇ • F h (w n h ) dx dy + 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i )C RK (w n )∆t l dx dy
Using again the fact that w is a classical solution we have

II = 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i )∂ t (w n h -w n ) dx dy + 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i )∇ • (F h (w n h ) -F (w n )) dx dy + 1 K T ∈T h j∈T i∈T T γ j (ψ j -ψ i )C RK (w n )∆t l dx dy
In order to give an upper bound to the last expression, we make use of the fact that in 2D the total number of elements in the mesh can be bounded by h -2 , the properties of ψ to deduce that ψ j -ψ i can be bounded by ∇ψ L ∞ (Ω) h, the fact that |T | ≤ C 0 h 2 , and the properties of the approximation. This leads to

|II| ≤ C(Ω, T h ) h -2 γ L ∞ (Ω) C 0 h 2 ∇ψ L ∞ (Ω) h C w (T h , w n ) h p + C(Ω, T h ) h -2 γ L ∞ (Ω) C 0 h 2 ∇ψ L ∞ (Ω) h C ∇F (T h , w n ) h p-1 + C(Ω, T h ) h -2 γ L ∞ (Ω) C 0 h 2 ∇ψ L ∞ (Ω) h C RK (w n )∆t l With γ L ∞ (Ω) = max T ∈T h max j∈T γ j L ∞ (T ) . Setting C 1 (Ω, T h , w n ) = C(Ω, T h )C 0 C ∇ψ max(C w (T h , w n ), C ∇F (T h , w n )
and C RK (Ω, T h , w n ) = C(Ω, T h )C 0 C ∇ψ C RK (w n ) we get the estimate

|II| ≤ γ L ∞ (Ω) C 1 (Ω, T h , w n ) h p + C RK (Ω, T h , w n ) h ∆t l
Assembling the Galerkin error and the error associated to the bubble, we get finally

E n ≤ C 0 (T h , w n ) h p + C RK (T h , w n ) ∆t p + γ L ∞ (Ω) C 1 (Ω, T h , w n ) h p + γ L ∞ (Ω) C RK (Ω, T h , w n ) h ∆t l (69)
This immediately shows that, provided that the bubble functions are uniformly bounded, we are allowed to have l ≤ p, in particular, it is enough to take l = p-1 to retain the accuracy of the Galerkin approximation.

In particular, if, as it is always the case for explicit schemes, we can find two positive bounded constants C h/∆t and C ∆t/h such that

C h/∆t ≤ ∆t h ≤ C ∆t/h
then we have for l = p -1

E n ≤ C h p (70)
with

C = C 0 (T h , w n )+C RK (T h , w n ) C p ∆t/h + γ L ∞ (Ω) C 1 (Ω, T h , w n ) + C RK (Ω, T h , w n )C p-1 ∆t/h
Note that, mass lumping is kept out of the analysis. However, as shown in section 3, it can be included in the definition of the (bounded) bubble function, at least in the P 1 case. For the higher order case, we refer to [START_REF] Cohen | High order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Jund | Arbitrary high order finite element schemes and high order mass lumping[END_REF][START_REF] Jund | Méthodes d'éléments finis d'ordre élevé pour la simulations numérique de la propagation d'ondes[END_REF] for more.

Appendix 2 : hypotheses 4.1 and 4.2

In this appendix we justify the choice of the approximate time incerements δu k (cf. section 4.1 and appendix 1) for the RK2 and RK3 schemes of section 4.1.

We recall that the constraint to respect is that for the last RK step

r k = δu k ∆t + f k = O(∆t l )
with l ≥ p -1, where with p we denote the (desired) overall accuracy of the scheme. The analysis will be performed for the autonomous ODE :

∂ t u + f (u) = 0 ( 71
)
RK2 scheme. Let us start by verifying hypothesis 4.1 for the RK2 scheme defined by

u 1 =u n -∆tf (u n ) u n+1 =u n - ∆t 2 f (u n ) - ∆t 2 f (u 1 )
When replacing u n and u n+1 by the values at t n and t n+1 of an exact solution w(t), and using the fact that w 1 = w n -∆tf (w n ), we can write

f (w 1 ) =f (w n + (w 1 -w n )) = f (w n ) + (w 1 -w n )∂ u f (u n ) + (w 1 -w n ) 2 2 ∂ uu f (u n ) + O (w 1 -w n ) 3 =f (w n ) -∆tf (u n )∂ u f (u n ) + ∆t 2 2 f (u n ) 2 ∂ uu f (u n ) + O(∆t 3 )
which immediately leads to

w n+1 -w n ∆t + 1 2 f (u n ) + f (u 1 ) =∂ t w n + ∆t 2 ∂ tt w n + ∆t 2 6 ∂ ttt w n +f (w n ) - ∆t 2 f (u n )∂ u f (u n ) + ∆t 2 4 f (u n ) 2 ∂ uu f (u n ) + O(∆t 3 ) =∂ t w n + f (w n ) - ∆t 2 3 1 2 f (w n )∂ u f (w n ) 2 + f (w n ) 2 ∂ uu f (w n ) + O(∆t 3 ) =
INRIA having used the relations

∂ t w n = -f (w n ) ∂ tt w n =f (w n )∂ u f (w n ) ∂ ttt w n = -f (w n )∂ u f (w n ) 2 -f (w n ) 2 ∂ uu f (w n ) ∂ tttt w n =f (w n )∂ u f (w n ) 3 + 4f (w n ) 2 ∂ u f (w n )∂ uu f (w n ) + f (w n ) 3 ∂ uuu f (w n ) ( 72
)
To verify hypothesis 4.2, we perform a similar exercise :

r n+1 = w 1 -w n ∆t + 1 2 f (w n ) + f (w 1 ) = -f (w n ) + 1 2 f (w n ) + f (w 1 ) = - 1 2 f (w n ) + 1 2 f (w n ) - ∆t 2 f (w n )∂ u f (w n ) + O(∆t 2 ) = - ∆t 2 f (w n )∂ u f (w n ) + O(∆t 2 ) = O(∆t)
proving that this definition of r n+1 is enough for use in the stabilization term in second order schemes.

RK3 scheme. We repeat the same exercise for the RK3 scheme defined by

u 1 =u n -∆tf (u n ) u 2 =u n - ∆t 4 f (u n ) - ∆t 4 f (u 1 ) u n+1 =u n - ∆t 6 f (u n ) - 2∆t 3 f (u 2 ) - ∆t 6 f (u 1 )
When replacing u n and u n+1 by the values at t n and t n+1 of an exact solution w(t), we can easily prove the following developments

f (w 1 ) =f (w n ) -∆tf (w n )∂ u f (w n ) + ∆t 2 2 f (w n ) 2 ∂ uu f (w n ) - ∆t 3 6 f (w n ) 3 ∂ uuu f (w n ) + O(∆t 4 ) f (w 2 ) =f (w n ) - ∆t 2 f (w n )∂ u f (w n ) + ∆t 2 4 f (w n )∂ u f (w n ) 2 + 1 2 f (w n ) 2 ∂ uu f (w n ) - ∆t 3 8 2f (w n ) 2 ∂ u f (w n )∂ uu f (w n ) + 1 6 f (w n ) 3 ∂ uuu f (w n ) + O(∆t 4 )
These developments can be readily used to show that

w n+1 -w n ∆t + 1 6 f (w n ) + 1 6 f (w 1 ) + 2 3 f (w 2 ) = ∂ t w n + ∆t 2 ∂ tt w n + ∆t 2 6 ∂ ttt w n + ∆t 3 24 ∂ tttt w n + 1 6 f (w n ) + 1 6 f (w n ) - ∆t 6 f (w n )∂ u f (w n ) + ∆t 2 12 f (w n ) 2 ∂ uu f (w n ) - ∆t 3 36 f (w n ) 3 ∂ uuu f (w n ) + 2 3 f (w n ) - ∆t 3 f (w n )∂ u f (w n ) + ∆t 2 6 1 2 f (w n ) 2 ∂ uu f (w n ) + f (w n )∂ u f (w n ) 2 - ∆t 3 12 2f (w n ) 2 ∂ u f (w n )∂ uu f (w n ) + 1 6 f (w n ) 3 ∂ uuu f (w n ) + O(∆t 4 )
which, using (72), leads immediately to

w n+1 -w n ∆t + 1 6 f (w n ) + 1 6 f (w 1 ) + 2 3 f (w 2 ) = ∂ t w n + f (w n ) + ∆t 3 12 f (w n )∂ u f (w n ) 3 + O(∆t 4 ) = O(∆t 3 )
A similar exercise can be used now to show that

r 2 (w) = w 1 -w n 2∆t + 1 4 f (w n ) + 1 4 f (w 1 ) = - 1 2 f (w n ) + 1 4 f (w n ) + 1 4 f (w n ) - ∆t 4 f (w n )∂ u f (w n ) + O(∆t 2 ) = - ∆t 4 f (w n )∂ u f (w n ) + O(∆t 2 ) = O(∆t)
and more importantly that

r n+1 (w) =2 w 2 -w n ∆t + 1 6 f (w n ) + 1 6 f (w 1 ) + 2 3 f (w 2 ) = 2 ∆t - ∆t 2 f (w n ) + ∆t 2 4 f (w n )∂ u f (w n ) - ∆t 3 8 f (w n ) 2 ∂ uu f (w n ) + 1 6 f (w n ) + 1 6 f (w n ) - ∆t 6 f (w n )∂ u f (w n ) + ∆t 2 12 f (w n ) 2 ∂ uu f (w n )+ 2 3 f (w n ) - ∆t 3 f (w n )∂ u f (w n ) + ∆t 2 6 f (w n )∂ u f (w n ) 2 + 1 2 f (w n ) 2 ∂ uu f (w n ) + O(∆t 3 ) = ∆t 2 6 f (w n )∂ u f (w n ) 2 - 1 2 f (w n ) 2 ∂ uu f (w n ) + O(∆t 3 ) = O(∆t 2 )
which shows that also for the RK3 scheme, our definitions of the r k do verify the accuracy constraint.

Figure 1 :

 1 Figure 1: Left. Formulation 3 : dual areas C j , |C j | = β j |T |. Right. Formulation 4 : area coordinates of the distribution point M ; |T | j = β j |T |

Figure 2 :

 2 Figure 2: Residual Based Finite Volume. Cells C i and C j , and staggered cell C ij .

Figure 3 :

 3 Figure 3: Typical topology of the meshes used in the numerical tests

Figure 4 :

 4 Figure 4: Scalar advection : grid convergence for the LDA scheme. Top-left : formulation F1. Top-right : formulation F2. Bottom-left : formulation F3. Bottom-right : formulation F4.

-

Figure 9 :

 9 Figure 9: 2d Burger's equation : LDA-F1-SL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 10 :Figure 11 :

 1011 Figure 10: 2d Burger's equation : LDA-F1-GL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 12 :

 12 Figure 12: 2d Burger's equation : LDA-F1-GL-RK3 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 13 :

 13 Figure 13: 2d Burger's equation : B-F1-SL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 14 :

 14 Figure 14: 2d Burger's equation : B-F1-GL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 15 :Figure 16 :

 1516 Figure 15: 2d Burger's equation : B-F1-SL-RK3 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 17 :

 17 Figure 17: 2d Burger's equation : SU-F1-SL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 18 :

 18 Figure 18: 2d Burger's equation : SU-F1-GL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 19 :Figure 20 :

 1920 Figure 19: 2d Burger's equation : SU-F1-SL-RK3 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 21 :

 21 Figure 21: 2d Burger's equation : Bc-SL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 22 :

 22 Figure 22: 2d Burger's equation : Bc-GL-RK2 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

Figure 23 :Figure 24 :

 2324 Figure 23: 2d Burger's equation : Bc-SL-RK3 scheme. Left : contours at time t = 1. Middle : solution along the line y = 0.3 and along the symmetry line. Right : minimum and maximum values of the solution.

. 6 log 10 Figure 25 :

 61025 Figure 25: Vortex advection : grid convergence for the LDA scheme with F1. Left column : convergence history. Right column : convergence rates.

. 6 log 10 Figure 26 :

 61026 Figure 26: Vortex advection : grid convergence for the B scheme with F1. Left column : convergence history. Right column : convergence rates.

. 6 log 10 Figure 28 :

 61028 Figure 28: Vortex advection : grid convergence for the Bc scheme. Left column : convergence history. Right column : convergence rates.

Figure 29 :Figure 30 :

 2930 Figure 29: Double Mach reflection. Density contours for the B scheme. 30 equally spaced contours from 1 to 24. Top row : RK2. Bottom row : RK3. Left column : selective lumping. Right column : global lumping.

Figure 31 :

 31 Figure 31: Double Mach reflection. Solution across the oblique shock (y = 0.6). Top row : density. Bottom row : entropy. Left column : B scheme. Right column : Bc scheme. ρ min ρ max

Figure 32 :

 32 Figure 32: Double Mach reflection. Solution along the wall (y = 0). Top row : density. Bottom row : entropy. Left column : B scheme. Right column : Bc scheme.

Figure 33 :ρ

 33 Figure 33: Mach 3 wind tunnel : close-up view of the mesh around the corner (h = 1/80 far from the corner, h = 10 -3 at the corner)

Figure 34 :

 34 Figure 34: Mach 3 wind tunnel : B-SL-RK2 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ

  min ρ max t = 0.5 0.065 5.75 t = 1.5 0.067 8.19 t = 4.0 0.071 6.31

Figure 35 :

 35 Figure 35: Mach 3 wind tunnel : B-SL-RK3 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ

  min ρ max t = 0.5 0.062 5.73 t = 1.5 0.067 7.89 t = 4.0 0.072 6.29

Figure 36 :

 36 Figure 36: Mach 3 wind tunnel : B-GL-RK2 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ

  min ρ max t = 0.5 0.062 5.73 t = 1.5 0.067 7.86 t = 4.0 0.072 6.30

Figure 37 :

 37 Figure 37: Mach 3 wind tunnel : B-GL-RK3 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ

  min ρ max t = 0.5 0.049 5.69 t = 1.5 0.065 8.04 t = 4.0 0.059 6.25

Figure 38 :

 38 Figure 38: Mach 3 wind tunnel : Bc-SL-RK2 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ 28 Figure 39 :

 2839 Figure 39: Mach 3 wind tunnel : Bc-SL-RK3 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ 31 Figure 40 :

 3140 Figure 40: Mach 3 wind tunnel : Bc-GL-RK2 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

ρ 33 Figure 41 :

 3341 Figure 41: Mach 3 wind tunnel : Bc-GL-RK3 scheme. Left : density contours at time t = 0.5 (top), t = 1.5 (middle), and t = 4.0 (bottom) ; 30 equally spaced contours between 0.5 and 8. Right : density distribution along the line y = 0.2, and minimum and maximum values of the density.

Table 1 :

 1 Double Mach reflection : minimum and maximum values of the density.

				ρ min ρ max
		1.40 22.30	Bc-SL-RK2	1.398 24.12
	B-SL-RK3	1.29 22.30	Bc-SL-RK3	1.4	24.07
	B-GL-RK2	1.37 22.30	Bc-GL-RK2	1.397 24.00
	B-GL-RK3	0.77 22.20	Bc-GL-RK3	1.396 23.98

RR n 6998

Ricchiuto & Abgrall

Acknowledgements

We are grateful to G. Scovazzi (SANDIA NL, Albuquerque) for his numerous remarks on the solution of the nonlinear Newton loop in RD discretizations, and for pointing out the relations with predictor-multicorrector schemes used in SUPG discretizations.

Appendix 1 : proof of proposition 4.3 We derive a sufficient condition on the r k guaranteeing that the accuracy of the Runge-Kutta Galerkin approximation is not lost when adding the bubble contribution. To do this we use a truncation error analysis, following the approach of [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF]. What we want to do is to show that for a p-th order spatial approximation, and when employing a p-th order RK scheme, the solution at time t n verifies a truncation error equation of the type E n ≤ C n h p . A global space-time truncation error estimate can then be obtained by integrating over each time interval and adding up over all the time slabs.

To do this we consider a more general polynomial approximation in space. The degrees of freedom are still approximations of the values of the unknown in some nodal locations of the mesh, except that, differently from the P 1 case, these are the element vertices plus other locations, as for example in standard P k elements, or in more "exotic" polynomial spaces, such as the ones proposed in [START_REF] Cohen | High order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Jund | Arbitrary high order finite element schemes and high order mass lumping[END_REF][START_REF] Jund | Méthodes d'éléments finis d'ordre élevé pour la simulations numérique de la propagation d'ondes[END_REF]. As before, we denote by ϕ i the basis functions spanning the polynomial space, by u h the spatial approximation of a function u, and with K we denote the total number of degrees of freedom (DoF) contained in an element of the mesh.

We start by recalling that, given a smooth classical solution of the problem w, hypothesis 4.1 guarantees that a p-th order RK scheme verifies the truncation error estimate

where, with the notation of section 4.1, we have explicitely used the fact that for the last stage of the RK scheme r k = r n+1 . Similarly, hypothesis 4.2 ensures