Stéphane Clain
email: clain@math.uminho.pt

Raphaël Loubère
email: raphael.loubere@u-bordeaux.fr

Gaspar J Machado

a posteriori stabilized sixth-order finite volume scheme for one-dimensional steady-state hyperbolic equations

We propose a new family of finite volume high-accurate numerical schemes devoted to solve one-dimensional steady-state hyperbolic systems. High-accuracy (up to the sixth-order presently) is achieved thanks to polynomial reconstructions while stability is provided with an a posteriori MOOD method which control the cell polynomial degree for eliminating non-physical oscillations in the vicinity of discontinuities. Such a procedure demands the determination of a chain detector to discriminate between troubled and valid cells, a cascade of polynomial degrees to be successively tested when oscillations are detected, and a parachute scheme corresponding to the last, viscous, and robust scheme of the cascade. Experimented on linear, Burgers', and Euler equations, we demonstrate that the schemes manage to retrieve smooth solutions with optimal order of accuracy but also irregular solutions without spurious oscillations.

Introduction

Numerical approximations of the steady-state Euler system date back to the early seventies with the NASA Ames group [START_REF] Pulliam | Early development of implicit methods for computational fluid dynamics at nasa ames[END_REF]. The development of implicit time algorithm [START_REF] Beam | An implicit finite-difference algorithm for hyperbolic systems in conservation law form[END_REF] with high-order finite difference methods [START_REF] Beam | Upwind second-order difference schemes and applications in aerodynamic flows[END_REF] led to the first simulations for two-and three-dimensional complex geometries calculated on Illiac IV Computer [START_REF] Steger | Implicit finite difference simulation of flow about arbitrary geometries with application to airfoils[END_REF][START_REF] Pulliam | Implicit finite-difference simulations of three dimensional compressible flow[END_REF][START_REF] Pulliam | Simulation of three-dimensional compressible viscous flow on the illiac iv computer[END_REF]. The steady-state was achieved as the limit stage of the non-stationary problem using ficticious time step (time marching method) [START_REF] Venkatakrishnan | Implicit solvers for unstructured meshes[END_REF]. The main difficulties are, on one hand, to achieve an accurate approximation where the solution is smooth enough, and, on the other hand, to produce a robust solution without non-physical oscillations.

Most of the high-order technologies developed for the non-stationary case were adapted to the steady-state case with the time marching algorithm using an explicit scheme [START_REF] Jameson | Numerical solution of Euler Equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] or an implicit formulation [START_REF] Yee | Implicit total variation diminishing (tvd) schemes for steady-state calculations[END_REF][START_REF] Maciel | Explicit and implicit tvd high resolution schemes in 2d[END_REF] equipped with a Newton-Krylov iterative procedure [START_REF] Venkatakrishnan | Implicit solvers for unstructured meshes[END_REF][START_REF] Venkatakrishnan | On the accuracy of limiters and convergence to steady state solutions[END_REF][START_REF] Nielsen | Application of newton-krylov methodology to a threedimensional unstructured Euler code[END_REF]. Taking the time step to infinity leads to a Newton-Krylov method for the steady-state problem which represents an alternative to the time marching method [START_REF] Knoll | Jacobian-free newton-krylov methods: a survey of approaches and applications[END_REF]. Second-order finite volume approximations with MUSCL limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme II: Monotonicity and conservation combined in a second order scheme[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme V: A second order sequel to Godunov's method[END_REF][START_REF] Venkatakrishnan | On the accuracy of limiters and convergence to steady state solutions[END_REF][START_REF] Venkatakrishnan | Convergence to steady state solutions of the Euler equations on unstructured grids with limiters[END_REF] on structured [START_REF] Yee | Application of second-order-accurate total variation diminishing (tvd) schemes to the Euler equations in general geometries[END_REF][START_REF] Jameson | Numerical solution of Euler Equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] or unstructured meshes [START_REF] Dervieux | Steady Euler simulation using unstructured meshes[END_REF][START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF][START_REF] Jameson | Finite volume solution of the two-dimensional Euler equations on a regular triangular mesh[END_REF][START_REF] Mavriplis | Multigrid solution of the two-dimensional Euler equations on unstructured triangular meshes[END_REF] revealed to be the standard approach in the nineties and have been widely developed. Unfortunatly, the accuracy remains below the optimal value two due to an excessive limiting even for smooth solution where limiting should not act toprovide the maximum order. Moreover, higher-order methods are appealing because they may dramatically reduce the computational time since they provide an excellent approximation with coarse meshes for smooth solutions. Of course, shocks and vortices still require an adequate local mesh to catch all the structures but high-order schemes prove to be efficient in reducing the numerical diffusion. For instance the k-exact reconstruction [START_REF] Barth | Recent developments in high order k-exact reconstruction on unstructured meshes[END_REF][START_REF] Barth | High order solution of the Euler equations on unstructured grids quadratic reconstruction[END_REF][START_REF] Mitchell | K-exact reconstruction for the navier-stokes equation on arbitrary grids[END_REF] increases the method accuracy using quadratic or cubic polynomial approximations [START_REF] Michalak | Matrix-explicit gmres for a higher-order accurate inviscid compressible flow solver[END_REF][START_REF] Nejat | Effect of discretization order on preconditioning and convergence of a high-order unstructured newton-gmres solver for the Euler equations[END_REF][START_REF] Nejat | A high-order accurate unstructured finite volume newton-krylov algorithm for inviscid compressible flows[END_REF]. Nevertheless, traditional TVD (Total Variation Diminishing) limiting procedures drastically reduce the order of accuracy despite the construction of alternative limiters [START_REF] Venkatakrishnan | On the accuracy of limiters and convergence to steady state solutions[END_REF][START_REF] Michalak | Limiters for unstructured higher-order accurate solutions of the Euler equations[END_REF] to enhance the quality of the solution.

More recently, Weighted Essentially Non-Oscillatory (WENO) technology [START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation[END_REF][START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF][START_REF] Jiang | Efficient implementation of weighted eno[END_REF][START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF] delivers more accurate approximations while preserving the robustness for the numerical solution and it enables to achieve higher-orders than traditionnal TVD methods. The technique was to a later time adapted for the steady-state situation [START_REF] Zhu | Runge-Kutta discontinuous galerkin method using WENO limiter II: Unstructured meshes[END_REF][START_REF] Li | A multigrid block lower-upper symmetric gauss-seidel algorithm for steady Euler equation on unstructured grids[END_REF]. The main difficulties lie on the non-differentiability of the numerical scheme leading to the creation of slight oscillations nearby shock discontinuities. More adapted WENO schemes then have been developed [START_REF] Zhang | A new smoothness indicator for the weno schemes and its effect on the convergence to steady state solution[END_REF][START_REF] Hu | A robust weno type finite volume solver for steady Euler equations on unstructured grids[END_REF][START_REF] Hu | An adaptive finite volume method for 2d steady Euler equations with weno reconstruction[END_REF] to overcome the problem of convergence toward the steady-state solution. We also mention some interesting developments for the steady-state case based on the Discontinuous Galerkin method [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF][START_REF] Bassi | High-order accurate discontinuous finite element solution of the 2D Euler equations[END_REF][START_REF] Pesch | A discontinuous galerkin finite element discretization of the Euler equations for compressible and incompressible fluids[END_REF], the Spectral Volume technique [START_REF] Breviglieri | Implicit high-order spectral finite volume method for inviscid compressible flows[END_REF][START_REF] Haga | An implicit lu-sgs scheme for the spectral volume method on unstructured tetrahedral grids[END_REF], and the Residual Distribution schemes [START_REF] Chou | High order residual distribution conservative finite difference weno schemes for steady state problems on non-smooth meshes[END_REF][START_REF] Abgrall | Development of residual distribution schemes for discontinuous galerkin method: the scalar case with linear elements[END_REF][START_REF] Hu | A robust high-order residual distribution type scheme for steady Euler equations on unstructured grids[END_REF].

All the technologies listed above employ an a priori limiting procedures to reduce the accuracy or the degree of the polynomial reconstructions in the presence of non-physical oscillations. Indeed, the reconstructed values are modified before feeding the numerical flux for the subsequent update stage. Consequently, unexpected over-limiting usually happens since the correction is carried out with the spirit to eliminate any potential risks. Moreover, the limiting procedure is only based on the numerical regularity of the approximation and rarely on the physics of the problem. For example, the positivity preserving is a consequence of the TVD restriction on the reconstruction but is not directly imposed as a restriction per se. The Multidimensional Optimal Order Method (MOOD) has been designed in [START_REF] Clain | A high-order finite volume method for systems of conservation lawsmultidimensional optimal order detection (MOOD)[END_REF][START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The MOOD method in the three-dimensional case: Very-high-order finite volume method for hyperbolic systems[END_REF] on different paradigms. It is an a posteriori method since the modification of the polynomial reconstruction is performed after calculating a candidate solution. The main avantage is that we only modify the solution where it is necessary leading to an improvement of the accuracy with respect to other methods. Another point is the capacity of the method to directly take into account the physics of the problem (positivity, entropy production, etc.). The MOOD paradigm has been developed for non-stationary hyperbolic system [START_REF] Zanotti | A posteriori subcell limiting for discontinuous galerkin finite element method for hyperbolic system of conservation laws[END_REF][START_REF] Loubère | A new family of high order unstructured mood and ader finite volume schemes for multidimensional systems of hyperbolic conservation laws[END_REF][START_REF] Clain | The mood method for the non-conservative shallow-water system[END_REF] and we aim at demonstrating that the methodology is adapted to steady-state situations. The key point is the introduction of an additional unknown vector which represents the maximum admissible polynomial degree on each cell. We seek for the numerical solution (density, the velocity, and the pressure for Euler system for instance) as well as the cell polynomial degree that produce the best approximation exempt from non-physical oscillations. In this work we only tackle the one-dimensional situation in order to focus on the MOOD algorithm and prove the ability of the method to provide high-order approximations for classical hyperbolic equation or system of equations. We also skip the question of the iterative solver efficiency since most of the implicit or explicit techniques used in the WENO, DG, or SV methods do apply in our context. This document is divided in seven sections. After the introduction, we present in section 2 the framework (notation, mesh, schemes) while section 3 is devoted to the MOOD methodology. We start by considering the simple scalar linear equations in section 4, then the non-linear Burgers' equation case in section 5, and dedicate section 6 to the Euler system. Numerical experiments involving regular and irregular solutions are presented for each equation or system of equations. We end the document with the conclusions of this work and present some perspectives.

Framework

Let Ω =]x L , x R [ be the computational domain discretized into a regular mesh T h constituted of cells K i = [x i- 1 2 , x i+ 1 2 ] for i = 1, . . . , I (see Fig. 1) of constant length h = (x R -x L )/I with center

x i = 1 2 (x i+ 1 2 + x i-1 2
). The boundaries of Ω are denoted by x 1 2 = x L and x

I+ 1 2 = x R . h x1 2 K 1 x 1 x3 2 x i-1 2 K i x i x i+ 1 2 x I-1 2 K I x I x I+ 1 2
Fig. 1 Notation: point/interface, mesh, and cells.

For any bounded function φ defined on Ω, φ i stands for an approximation of the mean value φ ex i of function φ over cell K i , that is

φ i ≈ φ ex i = 1 h x i+ 1 2 x i-1 2 φ(x) dx, i = 1, . . . , I,
while vector Φ = (φ i ) i=1,...,I ∈ R I gathers the values of the approximation.

Polynomial reconstruction

Let K i be a cell and d ∈ N. The stencil S i is composed of the n i closest neighbor cells (excluding cell K i ) such that n i > d. The polynomial of degree d associated to cell K i and approximation φ i is defined as

φ i (x; d) = φ i + d k=1 R i,k (x -x i ) k -X i,k , with R i = (R i,k ) k=1,.
..,d the vector that gathers the unknown polynomial coefficients. We set

X i,k = 1 h Ki (x -x i ) k dx to achieve the conservative property 1 h Ki φ i (x; d) dx = φ i . (1)
For a given stencil S i , we consider the quadratic functional

E i (R i ) = j∈Si 1 h Kj φ i (x; d) dx -φ j 2 .
We denote by R i = ( R i,k ) k=1,...,d the unique vector which minimizes the quadratic functional and set φ i (x; d) the associated polynomial that corresponds to the best approximation in the least squares sense of the data of the stencil

φ i (x; d) = φ i + d k=1 R i,k (x -x i ) k -X i,k .
Such a reconstruction satisfies the conservation property (1).

Remark 1

We have adopted the following notations: φ i represents an approximation of φ over cell K i while φ i (x; d) is a generic polynomial function of degree d associated to K i which satisfies the conservation property (1).

Generic finite volume scheme

We consider the steady state one-dimensional equation with source term

dF(φ) dx = S, in Ω, (2)
where F(φ) stands for the physical flux, S represents a regular source term, and Dirichlet boundary conditions φ D are prescribed on the boundaries. Integration of (2) over a cell K i gives

x i+ 1 2 x i-1 2 dF(φ(x)) dx dx = x i+ 1 2 x i-1 2 S(x) dx,
and integration by parts yields

F(φ(x i+ 1 2 ) -F(φ(x i-1 2 ) -h   1 h x i+ 1 2 x i-1 2 S(x) dx   = 0.
We substitute the exact relation with the numerical scheme

F(φ i+ 1 2 ,-, φ i+ 1 2 ,+ ; x i+ 1 2 ) -F(φ i-1 2 ,-, φ i-1 2 ,+ ; x i-1 2 ) -h S i = 0, (3)
where F(φ i+ 1 2 ,-, φ i+ 1 2 ,+ ; x i+ 1 2 ) is the two-point numerical flux evaluated at point x i+ 1 2 , scalar values φ i+ 1 2 ,-and φ i+ 1 2 ,+ stand for approximations on the left and right sides of interface x i+ 1 2 respectively, and S i is an approximation of the source term computed over the cell with Gaussian quadrature formula. Boundary conditions are taken into account for the first and last cell equations taking φ 1 2 ,-= φ D (x L ) and φ I+ 1 2 ,+ = φ D (x R ).

Cell polynomial degree and residual formulation

Polynomial reconstruction is an efficient tool to provide accurate approximations. Nonetheless, one has to pay attention to the local regularity and choose the appropriate degree for the representation. For a cell K i , we denote by d i its associated degree and M = (d i ) i=1,...,I stands for the Cell Polynomial Degree (CPD) map, that is, the vector which indicates the degree of the reconstruction for each cell. Moreover, we denote by

d i+ 1 2 = min(d i , d i+1
) the effective polynomial degree we use to compute approximations on interface x i+ 1 2 and set

φ i+ 1 2 ,-= φ i (x i+ 1 2 ; d i+ 1 2 ), φ i+ 1 2 ,+ = φ i+1 (x i+ 1 2 ; d i+ 1 2
). We introduce the residual G i (Φ, M) at cell K i that is given by

G i (Φ, M) = F(φ i+ 1 2 ,-, φ i+ 1 2 ,+ ; x i+ 1 2 ) -F(φ i-1 2 ,-, φ i-1 2 ,+ ; x i-1 2 ) -h S i , (4)
and the global residual G(Φ, M) ∈ R I gathers the residual expressions G i (Φ, M) for all cell i. For a given CPD map M, we introduce the residual operator

Φ ∈ R I → G(Φ, M) ∈ R I
and seek a solution Φ M for the system of equations G(Φ, M) = 0 using an iterative solver and given an initial guess. We also notice that the case where d i = 0 corresponds to the first-order approximation. On the other hand, if we are dealing with smooth solutions, let us choose d i = d with d given by the user and solving G(Φ, M) = 0 provides a d + 1th-order approximation. The point is to deal with solutions involving discontinuities. In this case the CPD map can not be constant and requires a low degree where the solution is rough and high degree where the solution is smooth. Since we do not know the regularity of the solution, M also turns to be an unknown of the problem and one search a solution (Φ , M ) such that G(Φ , M ) = 0.

Steady-state iterative solvers

Even dealing with linear equations, the problem is nonlinear due to the choice of M. Therefore iterative procedures have to be designed to produce a succession of approximations (Φ k , M k ) that converge to a satisfying solution (Φ , M ). Furthermore, if the system is nonlinear (Burgers equation, Euler system), an extra inner loop is required to solve the nonlinear problem for a fixed CPD map.

One way to proceed requires a solver which is problem dependent Φ → G(Φ, M) such that, being given a CPD map M k and an initial guess Φ k , the (non)linear ((N)L) solver computes Φ k+1 solution of G(Φ, M k ) = 0. Then a new CPD map M k+1 is provided on the base of the data of Φ k+1 (evaluation of M k+1 will be explained in the next section) and the iterative process goes one until we reach to a limit CPD map, i.e. M k = M k+1 since in that case, we have Φ k = Φ k+1 and the steady-state approximation is obtained.

A second strategy consists in reformulating equation (2) by introducing a fake time dependency

∂ φ ∂t + ∂F( φ) ∂x = S (5)
with prescribed initial and boundary conditions. The time parameter t is devoted to tend to infinity providing the steady state solution as φ = lim t→∞ φ(•, t) . Equation ( 5) is solved using an explicit high accurate finite volume code as presented in [START_REF] Diot | The MOOD method in the three-dimensional case: Very-high-order finite volume method for hyperbolic systems[END_REF] and we tag the solver as "Time Marching solver" (TM).

The MOOD method

The dynamical determination of the CPD map M is of paramount importance to ensure stability, robustness, and accuracy where appropriate. In the case of a discontinuous solution, high-order schemes produce numerical instabilities resulting into non-physical oscillations in the vicinity of the shocks. Consequently, in the presence of a discontinuous solution or a steep gradient, the cell polynomial degree d i must drop to zero to avoid spurious oscillations whereas the maximal userdefined degree d max can be employed in the zones where the solution is regular enough, ensuring de facto a local high accuracy of the approximate solution.

MOOD loop

We dynamical determine the CPD map using a a posteriori MOOD-like approach developed in [START_REF] Clain | A high-order finite volume method for systems of conservation lawsmultidimensional optimal order detection (MOOD)[END_REF][START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The MOOD method in the three-dimensional case: Very-high-order finite volume method for hyperbolic systems[END_REF]. Namely, for a given stage k and its associated map M k , a candidate solution Φ k is computed.

We then perform a detection procedure to identify which cells present numerical artifacts that would demand more dissipation. For those cells only we reduce the polynomial degree following a cascade, i.e. a decreasing sequence of polynomial degrees, keeping the good cell polynomial degrees alike. This determines the updated CPD map M k+1 which is further employed to compute the new candidate solution Φ k+1 . We depict in Fig. 2 the MOOD loop embedding the solver that provides solution Φ k for a fixed CPD map Starting from an initial guess Φ 0 and an initial CPD map M 0 (k = 0), the solver furnishes a candidate solution Φ k+1 which is tested (detection) and subsequent decrementing of the polynomial degree map may occur giving a new CPD map M k+1 . If so, or if convergence is not reached then another iterate is performed.

M k . M Φ k k Φ k+1 M k+1 M k+1 M k = M Φ 0 0 k=0 Initial guess POLYNOMIAL RECONSTRUCTION OUT ELSE if k = k+1

Detection chain

The detection chain is the procedure by which the algorithm analyses a candidate solution. For each cell K i , we apply the chain to assess if the value is considered as acceptable (with respect to the neighbor cells). If not, the solution needs to be recomputed with more dissipation locally by reducing the degree of the polynomial reconstruction associated to K i . The detection chain is composed of a succession of elementary detector procedures that focus on a specific potential problem of the solution. We list in the following the elementary detectors used in the applications.

Basic detectors

The basic detectors are designed to check if the candidate solution is physically admissible and to prevent the appearance of oscillations characterized by the creation of local extrema.

Physical Admissible Detector (PAD)

The candidate solution satisfies the PAD condition on cell K i if φ k i is a physical state associated with the solved system of PDEs. Such a condition is crucial for the Euler system and it consists of the positivity of the density and pressure values.

Extrema Detector (ED) One states that the candidate solution does not present an extremum at cell

K i if min(φ k i-1 , φ k i+1 ) ≤ φ k i ≤ max(φ k i-1 , φ k i+1 ).
When such a condition is violated, we have detected a potential oscillation of the discrete candidate solution.

Relaxation detectors

New extrema may derive from local oscillations associated to the Gibbs phenomenon, but may also correspond to smooth physical extrema. Another false detection situation arises when the solution is locally constant and extremely small perturbations due to the real number IEEE format may be interpreted as oscillations. Relaxation detectors aim at eliminating non-critical situations that have been falsely detected by ED. To this end, we compute an approximation of the local curvature setting C i (Φ) = φ i (x; 2), i = 1, . . . , I.

We then define the following local curvature indicators

χ i,m = min(C i-1 , C i , C i+1 ), χ i,M = max(C i-1 , C i , C i+1 ), i = 2, . . . , I -1,
where we omit dependency on the vector Φ for the sake of simplicity. Notice that one can have

|χ i,m | > |χ i,M |.
Plateau Detector (PD) Let ε PD be a given tolerance parameter. Then,

d i is not modified if max(|χ i,m |, |χ i,M |) ≤ ε PD .
Such condition means that the local curvature is so small that the numerical solution is locally linear and therefore the degree of the polynomial reconstruction should not be altered. Local Oscillation Detector (LOD) We must reduce the cell polynomial degree d i if one has

χ i,m χ i,M ≤ 0,
for any new extrema. This condition detects a local oscillation due to the variation of the curvature signs.

Smoothness Detector (SD) This detector determines if the minimum and the maximum curvatures are close enough with respect to the threshold parameter, meaning that the numerical solution is considered locally smooth. Let ε SD be a given tolerance parameter. The numerical solution is considered locally smooth if

1 ≥ min(|χ i,m |, |χ i,M |) max(|χ i,m |, |χ i,M |) ≥ 1 -ε SD .
If that is the case we do not modify the polynomial degree, otherwise the degree should be reduced.

In Fig. 3 we present a sketch of the chain of detectors proposed in this work. Note that some detectors may be removed depending on the system of equations (for instance the PAD may not make sense in the case of linear advection of a scalar). Note also that the PD permits to validate and accept a cell which is on a "plateau" and not wasting time trying to compute and compare zero valued curvatures in the subsequent detectors of the chain. The top/bottom arrows permit to exit the chain sooner if possible. In this work we fix ε PD = h and ε SD = 1/2.

The cascade

The last ingredient of a MOOD approach is to choose the cascade of polynomial degrees that the algorithm will successively try. A costly but systematic cascade consists in decrementing the polynomial degree by one from d max to 0. A less expensive choice consists in ignoring some intermediate stages. For example, d max = 5 → 2 → 0 is a less computational demanding cascade since we ignore the intermediate stages 4, 3, and 1. In this work we mainly use the cascade

d max = 5 → 2 → (1 →)0 or d max = 1 → 0.

Advection problem

We first consider the scalar linear convective steady-state problem

d dx u(x)φ(x) = S(x), x ∈ Ω, (6)
where u(x) is a given velocity function. We also prescribe boundary function φ D , which is just relevant at the boundary interfaces with inflow velocities. Since u or S may present some irregularities, the solution is not necessarily smooth enough to perform all the polynomial reconstructions with d max hence the CPD map will be adapted with respect to the numerical approximation.

Numerical flux

The numerical flux F(φ i+ 1 2 ,-, φ i+ 1 2 ,+ ; x i+ 1 2 ) is the classical upwind one. For a given CPD map M and noticing that the problem G(Φ, M) = 0 is affine, we use a GMRES procedure to nullify the residual and tagged this first strategy as "Linear solver" (L).

Next, the chain detector associated to this linear equation follows the chain depicted in Fig. 4 where all detectors have been defined in the previous section. We consider only two cascades:

d max = 1 → 0, or d max = 5 → 2 → 1 → 0,
that is second-order and sixth-order accurate schemes. The first-order accurate scheme with d max = 0 is also employed for comparison purposes.

Numerical experiments

The L 1 and L ∞ errors calculated for the whole mesh are given by

E 1 = I i=1 |φ i -φ ex i |h, E ∞ = I max i=1 |φ i -φ ex i |. (7)
Note that the previous norms explicitely depend on the numerical solution Φ and mesh size I, as E 1 (Φ, I) and E ∞ (Φ, I). Next we define the rate of convergence between two solutions/meshes (Φ k , I k ), for k = 1, 2 where I 1 < I 2 as

O α ≡ O α ((Φ 1 , I 1 ); (Φ 2 , I 2 )) = | log E α (Φ 1 , I 1 )/E α (Φ 2 , I 2 )| | log I 1 /I 2 | , α = 1, ∞. (8)
In order to derive an exact solution, we systematicaly choose φ(x) = 1/u(x), where u(x) is non-zero on the considered computational domain Ω = [0; 1]. This leads to u(x)φ(x) = 1, hence u(x)φ(x) = 0, for all x in Ω and S(x) = 0. Consequently, an exact solution associated to a given velocity field u is easily accessible. Moreover boundary conditions are given by φ 1 2 ,-= 1/u(0) and φ I+ 1 2 ,+ = 1/u(1), where u is known. We experiment two algorithms: the linear (L) and the time-marching (TM) solvers which both employ a MOOD like approach to stabilize the computation. The methodology of testing is based on the following test cases with the previous manufactured solutions.

Ω = [0; 1] -φ(x) = 1/u(x) -Left: regular solution u(x) = 2 sin(2πx) + 3 -Right: irregular solution u(x) = x + 2, if x < 1/2 or u(x) = 50(x -1/2) + 5/2, otherwise.
Regular solution (see Fig. 5-left). With this test, one tries to retrieve a regular steady-state solution of ( 6) with high accuracy. The unlimited version of the schemes are used to determine the maximal possible accuracy, then the MOOD-stabilized version of the schemes are employed and it is shown that optimal order of accuracy is attained. In other words we expect that no bad cells should be detected, leading de facto to a strict equivalence between unlimited and a posteriori MOOD schemes. Irregular solution (see Fig. 5-right). An irregular solution is simulated to prove the oscillatory character of the unlimited schemes, which, however is cured by the MOOD approach leading to a stable, yet more accurate, solution than the one provided by the first-order scheme.

Regular solution

Let us consider the case where u(x) = 2 sin(2πx) + 3 (see Fig. 5-left for a plot of the solution and the velocity) and the Linear and Time Marching schemes employ polynomial reconstructions P 0 , P 1 , and P 5 . That is to say first-, second-, and sixth-order accurate schemes. The boundary conditions are set to 1/3 at x = 0 and an outflow condition at x = 1. In Table 1 we present E 1 the errors in L 1 and E ∞ the L ∞ norms, the rates of convergence O 1 and O ∞ for L and TM schemes. Successively refined meshes made of I = 70 to 200 cells are considered. We also report the number of bad cells detected by the chain detector during the whole simulation (column labeled with the tag "B"). From the numbers in this table we observe that the optimal convergence rates are systematicaly attained. Moreover no bad cells are detected meaning that the a posteriori MOOD scheme would produce exactly the same results. As a consequence, because no bad cells were detected, a a posteriori MOOD L or TM scheme also reaches the optimal order of accuracy associated to its d max for this smooth test case.

Irregular solution

In this second test we consider velocity

u(x) = x + 2, if x ∈ [0; 1/2], 50(x -1/2) + 5/2, if x ∈ [1/2; 1],
leading to the solution depicted in Fig. 5 on right panel. This solution is continuous but not regular at location x = 1/2. The boundary conditions are set to 1/2 at x = 0 and an outflow condition at x = 1. The first set of simulations considers unlimited schemes and the results are presented in top panel of Fig. 6. Note that the L and TM solvers are producing the very same solutions, therefore we skip TM results (only the number of iterations, and, consequently the CPU time is larger for the TM scheme).

From the figures we deduce that the schemes are capturing the correct solution but with spurious oscillations for the high accurate schemes. On the same figure, on bottom panels, we plot the a posteriori MOOD limited solutions when the following cascades are employed: P 0 (left panel, red color), P 1 → P 0 (middle panel, blue color) and P 5 → P 2 → P 1 → P 0 (right panel, green color). We can observe that the numerical solutions are now exempt from spurious oscillations in the vicinity of the irregularity. There, cells are updated with lower order schemes employing either P 1 or P 0 reconstructions. Note that the color of the cell corresponds to its cell polynomial degree. Consequently we observe that for the cascade P 1 → P 0 two cells are updated with first order of accuracy, while, for the nominaly sixth-order accurate scheme two extra cells are updated with P 1 reconstructions. As such the convergence rate can not exceed a first order.

Nonetheless, it may still be of interest to employ high accurate schemes. Indeed, in Fig. 7, we present the logscale of the cell error for the first-, second-, and sixth-order accurate L schemes and three meshes I = 80, 240, and 320 cells. (The results for the TM scheme are equivalent.) As expected the schemes behave alike close to the discontinuity. There, the high accurate schemes are not more accurate than first order as expected. Contrarily, in the smooth parts of the flow, the high accurate schemes perform with an appreciable high accuracy. Note that, even if the discontinuity is located on one point, its effect on the numerical error spreads along a cone, the width of which is larger with the nominal accuracy of the scheme. Away from this cone the error produced by a high accurate scheme dramaticaly drops towards machine precision.

The Burgers' equation

Burger's equation is an attractive benchmark equation since we are dealing with a very simple scalar non-linear model which may provide discontinuous solutions. The equation writes 1 2 where S(x, φ) is a nonlinear source depending on the position and the solution.

dφ 2 dx = S(x, φ), x ∈ Ω, (9)

Numerical flux

The numerical flux F(φ i+ 1 2 ,-, φ i+ 1 2 ,+ ; x i+ 1 2 ) is the classical Rusanov one and the chain detector associated to the Burgers' equation follows the chain depicted in Fig. 4 for the linear case. For Burgers equation we consider only two cascades:

d max = 1 → 0, d max = 5 → 2 → 1 → 0,
that is second-order and sixth-order accurate schemes. Moreover the first-order accurate scheme is carried out for comparison purposes.

Numerical experiments

In this section we again experiment with the two algorithms, the nonlinear (NL) and the explicit time-marching (TM) solvers which both employ a MOOD like approach to stabilize the computation. The methodology of testing is based on simulating a regular solution (see Fig. 8-left) and an irregular one (see Fig. 8-right) following the experiments carried out for the linear equation. We adopt the same error norms and rate of convergence as for the linear equation, refer to ( 7)-(8).

Regular solution

In order to test the implementation of the numerical method we manufacture an exact solution.

Let us consider the regular function φ(x) = sin (3πx) exp x + 2 [START_REF] Breviglieri | Implicit high-order spectral finite volume method for inviscid compressible flows[END_REF] and formaly compute the source term

S(x) = 1 2
dφ 2 dx 2 = (sin (3πx) exp x + 2) (3π cos (3πx) exp x + sin (3πx)) .

We set the domain to be Ω = [0, 1] and compute the boundary conditions according to the exact solution, namely φ 1 2 ,-= φ(0) and φ I+ 1 2 ,+ = φ(1). In Table 2 we present E 1 the errors in L 1 and E ∞ the L ∞ norms, and the rates of convergence O 1 and O ∞ for the NL and TM schemes employing P 0 , P 1 , and P 5 polynomial reconstructions. Successively refined meshes made of 70, 80, 90, and 100 cells are considered. We also report the number of bad cells detected by the chain detector during the whole simulation with value B the number of bad cells detected by the MOOD loop.

From the table we can conclude that both solvers provide the optimal expected order of accuracy. The NL solver converges in about 10 outer iterations for all simulation with a residual at machine tolerance. The TM scheme demands from 4000 to 5000 iterations to converge at machine tolerance for P 0 and P 1 polynomial degrees due to its explicit nature while the P 5 reconstructions demands about 8000 to 10000 iterations to reduce the residual around 5 × 10 -5 , hence explaining the non optimal order of convergence. We illustrate this behavior in Fig. 9 where we plot on the bottom panels the evolution of the residuals as a function of the nonlinear solver iterations (NL scheme) or as a function of the time cycles (TM scheme). From this plot we can observe that the evolution of the residual for the TM scheme presents a plateau at 0.1 during few thousands of cycles before dropping to another plateau at ∼ 10 -6 . For all other simulations, that is with different d max = 0, 1, the same plateau at 0.1 is systematicaly observed, but the dropping reaches machine precision. On top panel of Fig. 9 we plot the solution provided by both schemes for 100 cells and P 5 reconstructions. As expected, even if machine precision has not been reached for the TM scheme, the solution is visually indistinguishible from the NL converged one. Note that the solution is visually improved when the mesh is refined and/or the maximal polynomial degree d max increases (not shown). Also we have observed that the MOOD detector does not flag any cell during the whole simulation, B = 0 systematically in Table 2, proving that a limiting procedure is not needed. As such MOOD limited TM or NL schemes provide exactly the same results as the unlimited ones. Consequently we do not show these tables and curves.

We have also tested other regular solutions such as polynominals, exponentials, etc., and the same conclusions hold. Consequently we also omit those results.

Irregular solution

Next, following [START_REF] Chen | Fast sweeping methods for steady state problems for hyperbolic conservation laws[END_REF] for instance, we will consider a discontinuous solution. Considering the source term S(x) = -π cos(πx)φ(x) and the boundary conditions φ 1 2 ,-= 1 and φ I+ 1 2 ,+ = -0.1, we can determine the exact solution with shock

φ(x) = 1 -sin(πx), if x ∈ [0; x s ], -0.1 -sin(πx), if x ∈ [x s ; 1], (12)
where x s is the location of the shock. Two shocks are possible -x s = 0.1486 and x s = 0.8514but only the first shock is stable for small perturbation. Exact solution is imposed at the boundaries of the domain. We plot in Fig. 8 the solutions (black line) along with the source term (red line).

Unlimited solvers. Let us first observe the behavior of unlimited TM and NL solvers. We represent in Fig. 10 the solutions produced by the NL and TM solvers when d max = 0, 1, and 5 and I = 70 cells are considered. We observe that all simulations are capable of capturing the shock wave but only the schemes without reconstructions (P 0 ) produce monotone solutions. Any high accurate unlimited scheme, P 1 , P 5 for TM or NL solvers does produce parasitical oscillations in the vicinity of the shock as expected. The same behavior is observed when the mesh is refined up to I = 100 cells. This behavior is confirmed by the number of bad cells detected which is of the order 5 for P 1 and P 5 simulations, and, of course 0 for P 0 simulation. All simulations have converged with machine precision residual apart from the TM P 5 simulations for which the residual saturates at about 10 -4 .

Both solvers, NL and TM, are performing alike in terms of accuracy and, from now on we will only use the NL solver to pursue our test campaign.

2nd order limited MOOD NL solver. In this paragraph we only test the 2nd order limited MOOD NL solver, that is, with the cascade d max = 1 → 0. The MOOD loop does produce several candidate solutions, each at convergence of the residual and an associated cell polynomial degree map that we plot in Fig. 11. As observed, for I = 120 cells, five MOOD iterations are needed to produce a valid solution in terms of the detection citeria. In red we emphasize the bad cells, updated with a first order accurate scheme, while in blue the cells updated with maximal 2nd order accuracy.

6th order limited MOOD NL solver. Next we test the 6th order limited MOOD NL solver with the cascade d max = 5 → 2 → 0. The MOOD loop produces several candidate solutions. Each candidate solution is obtained at convergence of the residual, and, has an associated cell polynomial degree map that we plot in Fig. 12 (colours). For I = 120 cells, ten MOOD iterations are needed to produce a valid solution in terms of the detection citeria. In red we emphasize the bad cells, updated with a first order accurate scheme, while in purple the cells updated with maximal 3rd order accuracy and in green the cells updated with the 6th order scheme. We present nine iterations in the figure, the final solution being plotted on the larger bottom panel. We observe that the bad cells are located in the vicinity of the shock wave as expected. About five bad cells for the final iterations need limiting. Note that the first MOOD iteration produces a solution which is not correct, but, when used as initial guess for the second iteration, it permits to capture the correct solution in the smooth parts. Then the other iterations fix the shock region.

Mesh refinement for the limited MOOD NL solver. Next, in Fig. 13, top panels, we present the behavior of the 6th order scheme when the mesh is refined (80, 120, 240, and 320 cells). The number of MOOD iterations is respectively 8, 10, 10, and 13. From this figure we can clearly observe that no spurious oscillations are present, while they were produced in previous iterations, and the number of corrected cells in the vicinity of the shock is about 4. For comparison purposes we have added on bottom panels the results obtained with a different cascade adding P 1 reconstructions, i.e. Last, we display in Fig. 14 the absolute value of the difference between the exact and the approximate solution for each cell for three of the previous meshes in logscale and for cascades P 1 → P 0 and P 5 → P 2 → P 0 . Due to the discontinuous nature of the solution we can not expect more than first-order accuracy. Nonetheless we observe that the errors visually decrease with the iteration number and the mesh size. The maximal error is appropriately located in the vicinity of the shock wave, but the smooth parts of the flow are genuinely better approximated for the last iterations.

P
The a posteriori MOOD limiting is preserving the non-oscillatory character of the numerical solution for the highest order schemes. We observe that the procedure demands several iterations to reach a steady cell polynomial degree map. Nevertheless, once found, this map ensures the correctness of the solution by having the cell polynomial degrees reduced in the vicinity of the shock wave but nowhere else. We remind that the MOOD limiting de facto produces a (very) good guess. A very last step of the algorithm consists in using this guess and reseting the cell polynomial degree map to its maximal value.

In the case of the time-marching solver the MOOD procedure has more difficulties to converge as a "blinking" phenomena occurs due to the non-differential character of the method with respect to the CPD map (we jump from one degree to another one). Some pairs of cells alternatively drop then increase their degree so that the cell polynomial degree map can never attain a steady situation which is mandatory for convergence. This drawback could be solved (unsatisfactory) by setting a maximal number of iterations after which the map is kept fixed. Doing so the scheme converges at the price of fixing such a sensitive parameter.

The Euler system

We seek V ≡ V (x) = (ρ, u, p) T ≡ (ρ(x), u(x), p(x)) T , the vector of state variables density, velocity, and pressure, solution of the steady-state one-dimensional Euler equations with source terms given by dF (V ) dx = S, in Ω, [START_REF] Clain | A high-order finite volume method for systems of conservation lawsmultidimensional optimal order detection (MOOD)[END_REF] where the flux vector is given by F (V ) = (ρu, ρu 2 + p, u(E + p)) T and the regular source vector term is given by S ≡ S(x) = ( dD dx , dF dx , dH dx ) T ≡ ( dD dx (x), dF dx (x), dH dx (x)) T . Total energy E, specific internal energy e, and sound-speed a are expressed as

E = 1 2 ρu 2 + ρe, e = p ρ(γ -1) , a = γp ρ , (14)
respectively, where γ > 1 is the ratio of specific heat. System ( 13) is completed with Dirichlet boundary conditions

V (x L ) = V L ≡ (ρ L , u L , p L ) T and V (x R ) = V R ≡ (ρ R , u R , p R ) T ,
depending on the flow regime at x L and x R . On the other hand, U = (ρ, ρu, E) T stands for the the conservative variables vector.

We deduce the algebraic system to be solved

ρu = D, (15)
ρu 2 + p = F, (16)
u(E + p) = H, ( 17
) 1 2 ρu 2 + p γ -1 = E. (18)
Notice that D, F , and H depend on additional constants we shall fix with the boundary conditions.

Solution for the Euler system

We aim to derive exact solutions for system (15)-( 18) that we shall employ as manufactured solutions in the numerical section. First of all, functions F and H have to respect some restrictions to provide a physical solution. From ( 15) and ( 16), we deduce that p = F -D 2 ρ . The pressure and density non-negativity yields that

F ≥ 0, ρ ≥ D 2 F . ( 19
)
On the other hand from relations ( 17) and ( 14), we deduce that

p = γ -1 γ H u - 1 2 ρu 2 .
The non-negativity of the pressure and ( 15) implies

p = γ -1 γ H u - 1 2 ρu 2 ≥ 0 =⇒ Hu ≥ 1 2 ρu 4 =⇒ DH ≥ 1 2 D 4 ρ 2 .
(20)

Supercritical and subcritical solutions

Eliminating p and u from system ( 15)-( 18) provides an implicit relation between D, F , H, and ρ given by

F = γ + 1 2γ D 2 ρ + γ -1 γ ρH D , (21)
such that ρ is a positive solution. Equation ( 21) is rewritten as a quadratic polynomial in ρ of the form

(2(γ -1)H) ρ 2 -(2γDF ) ρ + (γ + 1)D 3 = 0. ( 22
)
Existence of solution demands that

4D 2 (γ 2 F 2 -2(γ -1)(γ + 1)DH) ≥ 0 =⇒ DH F 2 ≤ γ 2 2(γ 2 -1) (23)
leading to an additional restriction. From restrictions ( 19), [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF], and ( 23) we define the admissible space

Λ = (D, F, H), F ≥ 0, DH ≥ 0, γ 2 F 2 -2(γ 2 -1)DH ≥ 0 . (24)
Assuming (D, F, H) ∈ Λ, the supersonic solution ρsup (D, F, H) and the subsonic solution ρsub (D, F, H) correspond to the solution of ( 22) and are given by

ρsup (D, F, H) = γF -γ 2 F 2 -2(γ 2 -1)DH 2(γ -1) H D , and
ρsub (D, F, H) = γF + γ 2 F 2 -2(γ 2 -1)DH 2(γ -1) H D , respectively.
Remark 2 Solutions of the steady-state problem are constituted of branches of the supersonic and subsonic solutions with possible jump transitions. Consequently, since boundary conditions on the left and right sides have to satisfy the same relations ( 15)- [START_REF] Haga | An implicit lu-sgs scheme for the spectral volume method on unstructured tetrahedral grids[END_REF], they also correspond to subcritical or subcritical states.

Remark 3

We also deduce an important relation between the total energy and the density. Indeed, eliminating ρu 2 between ( 16) and ( 18) and mutipling ( 17) by ρ, we obtain

F -p 2 + p γ -1 = E, E + p = H D ρ.
Eliminating the pressure between the two last equations gives an affine relation between the total energy and the density

E = γ -1 γ + 1 F + 3 -γ γ + 1 H D ρ. (25)

Shock condition

We now analyse the condition to provide an admissible entropic steady-state shock.

Theorem 4 Assume that a steady-state solution admits an entropic stationary genuinely nonlinear shock at point x c . Then there only exist two admissible situations:

(a) the solution is supercritical on the left and subcritical on the right with D > 0 and H > 0;

(b) the solution is subcritical on the left and supercritical on the right with D < 0 and H < 0.

Proof. We present the proof for the case D > 0 and H > 0, i.e. the velocity is positive. A steadystate genuinely nonlinear shock requires that the flow regime on the left and righ sides of x c are differents. We state that a subcritical state on the left and supercritical state on the right is not possible. Indeed, a supercritical flow on the right side u > 0 means that three simple waves go out of the shock while only two simple waves enter on the left side (u and u + c). Consequently the shock is not entropic in the Lax sense. On the contrary, if the supercritical state lies on the left side, three simple waves enter in the shock while only two waves go out in accordance to the entropy principle. The entropic steady-state shock is then associated to the simple wave u -c.

One difficulty when dealing with steady-state solution is that the shock can take place at any point x c of the domain. The transition location is not contained in system ( 15)-( 18) and requires an additional restriction. To this end, we consider the steady-state situation as the limit case of the non-stationary problem solution U = ( ρ, ρu, E) T in conservative variables with initial conditions ρ 0 , ρ 0 u 0 , and E 0 . We prescribe the same Dirichlet boundary conditions as in the steady-state situation.

Theorem 5 Assume that density ρ(x, t), the mass flow ρu(x, t), and the total energy E are continuous at the boundary for t > 0, that is, there are no jumps between the boundary conditions and the solution at x L and x R (as a limit from the inner domain) for ρ, ρu, and E. Then for any time t > 0, we have Similar operations provide the two other relations.

Assuming that the unstationary solution converges to a limit steady-state solution, we get

xR xL ρ(x) dx = xR xL ρ 0 (x) dx, xR xL ρ(x)u(x) dx = xR xL ρ 0 (x)u 0 (x) dx = xR xL D(x) dx, xR xL E(x) dx = xR xL E 0 (x) dx.

Remark 6

The second relation is automatically satisfied if one assumes that for ρ 0 > 0 given, we define u 0 such that ρ 0 u 0 = D. In the following, we assume that the velocity u 0 satisfies relation [START_REF] Dervieux | Steady Euler simulation using unstructured meshes[END_REF]. In the same way, we assume that E 0 is compute with the relation [START_REF] Jameson | Numerical solution of Euler Equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] and provides

E 0 = γ -1 γ + 1 F + 3 -γ γ + 1 H D ρ 0 . (26)
Corollary 7 Assume that the steady-state solution is the limit of the unstationary problem with a unique discontinuity at point x c . Moreover, assume that the initial density ρ 0 satisfies the condition

Assume at last that there exists C ∈ R such that on the whole domain

H(x) D(x) = C. (28)
Then we have only two situations.

Proof. We only consider the first case. Steady-state solutions have to satisfy relation [START_REF] Hu | An adaptive finite volume method for 2d steady Euler equations with weno reconstruction[END_REF], hence are composed of subsonic and supersonic branches we join with discontinuities. Under the assumptions given in Theorem 5, mass conservation reads

xR xL ρ(x) dx = xR xL ρ 0 (x) dx.
From Theorem 4, the entropic shock requires that the supercritical solution lies on the left side of the shock situated at point x c . Hence ρ = ρ sup for x < x c while ρ = ρ sub for x > x c leading to relation [START_REF] Loubère | A new family of high order unstructured mood and ader finite volume schemes for multidimensional systems of hyperbolic conservation laws[END_REF]. Let

A sup = xR xL ρ sup (x) dx, A sub = xR xL ρ sub (x) dx.
Since ρ sup < ρ sub , function

x c → xc xL ρ sup (x) dx + xR xc ρ sub (x) dx ∈ [A sup , A sub ]
is strictly increasing with respect to x c . Consequently condition ( 27) is mandatory to provide existence and uniqueness of x c such that (29) holds. Relation [START_REF] Maciel | Explicit and implicit tvd high resolution schemes in 2d[END_REF] is obtained in a similar way.

The last point we have to check is that we must respect the additional constraint for the total energy. Since the initial total energy is defined by ( 26), we have, using relation ( 25)

xR xL E(x) dx = xR xL γ -1 γ + 1 F + 3 -γ γ + 1 H D ρ dx = γ -1 γ + 1 xR xL F + 3 -γ γ + 1 C xR xL ρ dx = γ -1 γ + 1 xR xL F + 3 -γ γ + 1 C xR xL ρ 0 dx = xR xL E 0 (x) dx.
Consequently, the additional constraint is automatically fulfilled.

Remark 8 Notice that the corollary requires that H D is constant otherwise, the compatibility does not hold. In that case, continuity of the density and total energy at the boundary no longer holds. Corollary 7 suggests the additional condition to fix the shock position by imposing the total mass. Therefore, the system (13) will be augmented with the relation

xR xL ρ(x) dx = A, (31)
for a given A ∈ [A sup , A sub ].

Remark 9 Since the boundary conditions have to correspond to one of the two states (subsonic or supersonic), we have 4 possible situations. Attending to Theorem 4 and Corollary 7, and assuming D > 0 and H > 0, we have (a) supersonic states on both boundary points -the solution for this case will be the regular function ρ sup ; (b) subsonic states on both boundary points -the solution for this case will be the regular function ρ sub ;

(c) supersonic state on the left boundary and subsonic state on the right boundary point -one shock sharing the supersonic solution situated on the left and the subsonic one on the righ; (d) subsonic state on the left boundary point and supersonic state on the right boundary point -the case is impossible since no entropic shock transition is possible from a subcritical to a supercritical branch.

High accurate numerical scheme

The high-accurate scheme is derived from the one proposed for Burgers' equation. Namely, one has to define the numerical flux to plug into the residual formulation. Let us denote by U = (U i ) i=1,...,I the 3×I matrix containing the approximations of the mean-values for ρ, ρu, and E while U k stands for the data associated to stage k. The numerical flux F(U -, U + ; x) is the HLL one proposed in [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF] where U -and U + represent the left and right states at interface x.

For a given CPD map M k , we introduce the nonlinear operator

U ∈ R 3×I → G(U , M k ) ∈ R 3×I
and we seek for U k+1 such that G(U k+1 , M k ) = 0. Again we consider either a nonlinear solver (NL) and the time-marching (TM) technique to solve the set of equations. Notice that for the NL case, an additional condition is considered to fix the shock position when dealing with two different regimes on the left and right boundaries. For the TM method, the initial condition for density will fix the position of the shock. The chain detector associated to the Euler equation follows the chain depicted in Fig. 15. Here the PAD box checks the positivity of density and pressure. The MOOD procedure is adapted to this system of equations. Because we consider only one polynomial degree per cell for each component of U i.e. density, mass flow and total energy, we determine only one polynomial degree map by only testing the density variable for the ED, LOD, and SD detection criteria. Consistently with Burgers' equation experiments we consider only two cascades:

d max = 1 → 0 or d max = 5 → 2 → 1 → 0. The MOOD loop is escaped if M k+1 = M k .

Numerical experiments

To assess the validity of the method, we consider the following data. The domain is given by x L = 0, x R = 1 and we take D(x) = 1, F (x) = 0.027x + 0.6137, and H(x) = 0.375. The first two conditions for the admissible space ( 24) are trivially satisfied. The third condition, DH

F 2 ≤ γ 2 2(γ 2 -1)
is also satisfied since max

x∈[0,1] DH F 2 = max x∈[0,1]
0.375 (0.027x + 0.6137) 2 = 0.9957 ≤ γ 2 2(γ 2 -1) = 1.0208.

The expressions for the supersonic solution and the subsonic solution are given by ρ sup (x) = 0.126x -3.3333 1.96(0.027x + 0.6137) 2 -0.72 + 2.8639, [START_REF] Michalak | Matrix-explicit gmres for a higher-order accurate inviscid compressible flow solver[END_REF] ρ sub (x) = 0.126x + 3.3333 1.96(0.027x + 0.6137) 2 -0.72 + 2.8639, [START_REF] Michalak | Limiters for unstructured higher-order accurate solutions of the Euler equations[END_REF] respectively and depicted in Fig. 16. Fig. 16 Plots of the supersonic and subsonic branches ρsup (red) and ρ sub (blue) from [START_REF] Michalak | Matrix-explicit gmres for a higher-order accurate inviscid compressible flow solver[END_REF], and the exact solution of the irregular density (shock located at x = 0.6) employed in this paper (black).

As mentioned in Remark 9, only three situations have to be considered. We have two regular situations if considering supersonic or subsonic states on both boundary points. The case where we prescribe subsonic state at the left boundary point and supersonic state at the right does not provide an entropic solution and must be discarded. At last, prescribing supersonic state on the left boundary point and subsonic state on the right boundary produces a solution with a shock wave located at x c . We adopt the same error norms and rate of convergence as for the linear and Burgers equations, refer to ( 7)- [START_REF] Beam | An implicit finite-difference algorithm for hyperbolic systems in conservation law form[END_REF], and only monitor the error for the density variable.

Regular case

We first try to retrieve a regular steady-state solution of ( 13) with high accuracy. The unlimited versions of the schemes are used by canceled the MOOD stage to determine the maximal possible accuracy. Then the MOOD-stabilized version of the scheme is employed and it is shown that optimal order of accuracy is attained. We consider the regular supersonic/supersonic case which corresponds to the bottom-red curve of Fig. 16 and regular meshes made of I = 25, 50, 75, and 100 cells and report in Table 3 the errors when the first-, second-, and sixth-order accurate MOOD NL and explicit TM schemes are used. Note that no bad cells are detected for this test case which implies that the unlimited or a posteriori MOOD limited schemes behave strictly equivalently. In the same way, the simulations for the regular solution on the subsonic/subsonic branch lead to the same conclusions and we omit those tables.

Discontinuous case

A solution with a steady shock wave is simulated to prove that the oscillatory character of the unlimited scheme is cured by the MOOD approach leading to a stable and more accurate solution than the first order scheme. To fix the position at x c = 0.6, we set

A = xR xL ρ 0 (x) dx = xR xL (ρ L -ρ R )(1 -x) 0.4550 + ρ R dx = 2.8975.
The exact solution for the density variable is depicted in Fig. 16. We enforce the position by using the additional restriction [START_REF] Mavriplis | Multigrid solution of the two-dimensional Euler equations on unstructured triangular meshes[END_REF]. Notice that the restriction is not necessary for the time marching method if we use ρ 0 and u 0 = D/ρ 0 as initial conditions. p 0 derives from relation ( 16) and E 0 is deduced from relation [START_REF] Haga | An implicit lu-sgs scheme for the spectral volume method on unstructured tetrahedral grids[END_REF].

Parachute and first-order accurate schemes In Fig. 17 we present the numerical solutions for density variable and I = 25, 50, and 100 cells for the NL scheme with polynomial degree 0. From this plot we observe a convergence towards the exact solution when the mesh size increases without any spurious oscillation. This validates the robust and least accurate scheme that is used in the following as the parachute scheme of the cascade. On the right panel of Fig. 17 are displayed the logscale of the cell error for meshes of size I = 50, 75, 100, and 200. Because the scheme is nominaly of first order of accuracy we observe on the smooth regions that the error decreases by a factor 2 from the 50 to 100 sizes. Obviously the maximal error is located in the vicinity of the shock and do not decrease when the mesh is refined, only the "numerical region of influence" reduces in size, nonetheless the number of cells in this region remains constant, about 4. The same study is carried out for the first-order accurate scheme with cascade P 1 → P 0 and the results are reported in Fig. 18. The same conclusions applies and the error in the smooth parts of the flow are about two orders of magnitude smaller.

a posteriori MOOD NL schemes. Next we test the nominaly sixth-, second-, and first-order accurate NL scheme, that is, schemes using a maximal polynomial degree of 5, 1, and 0 and a MOOD loop to drop the degree to 0 if needed. The cascade of degrees are P 5 → P 2 → P 1 → P 0 for the sixth-order scheme and P 1 → P 0 for the second-order one. The numerical densities obtained for 75 cells are shown in Fig. 19 and the cells are colored according to their polynomial degree. The first-order accurate scheme needs one iteration while the second-and sixth-order accurate schemes demands two and nine respectively. We conclude that all the schemes produce a non-oscillatory solution, capture the correct shock location, employ the most accurate reconstructions far from the shock region, and decrement the polynomial degree of only few cells in the vicinity of the shock wave. We observe that the first iteration of any high-accurate scheme presents spurious oscillations. Then, the MOOD loop further decreases the polynomial degree up to capturing a monotone solution. The exact same behavior is observed when the mesh size is increased or decreased, so we skip these figures. The sixth-order accurate scheme can reproduce the sharp shock wave obtained by the second-order accurate scheme, by restricting its spread over one cell, while the first-order accurate scheme demands about three cells. An interesting behavior can be observed when the cascade skips the P 1 reconstruction, that is when the choices are restricted to P 5 → P 2 → P 0 . Formaly the scheme is still nominally sixth order accurate. The results are displayed in Fig. 20 and we observe that the spread of the shock is ressembling the first-order scheme from top-panel of Fig. 19. We conclude that it is rather important to let the numerical method to choose P 1 reconstructions when appropriate.

Next, we present in Fig. 21 the logscale of the cell error produced by the three schemes (first-, second-, and sixth-order accurate) for 75, 100, and 200 cells at convergence, that is for the final MOOD iteration. The left panel results are associated to the right-most curves of Fig. 20. We observe that the highest the maximal polynomial degree, the lower the error in smooth regions, Mesh convergence for the sixth-order NL scheme. In Fig. 22 we present the numerical solutions for density variable and I = 25, 50, and 100 cells for the NL scheme a posteriori limited with maximal degree 5 and an a posteriori MOOD correction loop. From this plot we observe a convergence towards the exact solution when the mesh size increases and no spurious oscillations are observed. This shows that the nominaly sixth-order accurate numerical method is able to converge towards the exact solution. Obvisouly the convergence rate is still limited to the first-order one but, as previously seen, away from the shock wave, better accuracy can be expected. On the right panel of Fig. 22 are displayed the logscale of the cell error for meshes of size I = 50, 75, 100, and 200. Because the scheme is nominally of sixth order of accuracy we observe on the smooth region to the left away from the cone of influence, a drop of error of a factor ∼ 100 when the mesh size is doubled (from 50 to 100 and 100 to 200). This roughly corresponds to 2 6 , the factor expected for a sixth-order accurate scheme.

Conclusion and perspectives

A new strategy to construct high accurate finite volume numerical schemes has been presented to solve the one-dimensional steady-state linear, Burgers, and Euler equations. High accuracy is gained via polynomial reconstructions whereas, in the presence of steep gradients or shock waves, the numerical schemes are stabilized with an a posteriori MOOD strategy. The technique consists in detecting bad cells on a candidate solution, and further recomputing the solution using lower degree polynomial reconstructions in the vicinity of the bad cells. Doing so we define: (1) a chain of detectors which determines if a cell is valid or not, (2) a cascade of polynomial degree, that is a list of ordered polynomial degrees that one wishes to try, and (3) a parachute scheme that always produces a valid solution. The chain detector is based on classical physical admissibility, local smoothness, plateau, and extrema detection. The parachute scheme is always the first-order finite volume scheme. The cascades are of two kinds, a nominaly second-order accurate one P 1 → P 0 , and a sixth-order accurate one P 5 → P 2 (→ P 1 ) → P 0 .

To perform the simulations, two solvers are considered: a non-linear implicit solver and a so-called explicit pseudo-time marching solver. The one-dimensional scalar linear equation, the Burgers scalar non-linear equation, and the one-dimensional Euler system with source terms were simulated. For each system of equation, several exact regular and irregular solutions are derived. The unknows of these problems are not only the cell-centered finite volume values of the solution, but also the cell polynomial degree (CPD) map.

In the vicinity of a discontinuous solution, a low-order scheme must be employed (i.e. for instance using a polynomial of degree 0) to avoid Gibbs phenomenon, whereas the maximal degree is expected to produce a valid and acurate solution on smooth parts of the flow. An iterative MOOD loop is used to compute a candidate solution associated to its CPD map, then the chain detection checks for bad cells and modifies the CPD map accordingly. A new candidate solution is then recomputed and again checked for possible detected troubled cells.

Numerical experiments show that the optimal order of accuracy is systematicaly reached for smooth solutions. For non-smooth solutions, the a posteriori MOOD stabilization leads to nonoscillatory solutions while maintaining reasonable accuracy on the smooth parts of the flow. As expected the troubled cells are mostly located in the vicinity of discontinuities. We have numerically observed that the stabilization technique performs for a nominaly second-order scheme but also for the sixth-order one without visible spurious effects.

The goal of the study is to present a proof of concept that the a posteriori MOOD stabilization strategy can be applied to the steady-state case. In the near future we plan to extend this approach to multi-dimensional situations taking into account classical boundary conditions (inflow, outflow, wall, symetry) and to efficient solvers. The results for the one-dimensional geometry indicate that the method performs well and motivate to tackle the more complex situations involving two-and also three-dimensional geometries.

Fig. 2

 2 Fig.2Sketch of the MOOD loop. Starting from an initial guess Φ 0 and an initial CPD map M 0 (k = 0), the solver furnishes a candidate solution Φ k+1 which is tested (detection) and subsequent decrementing of the polynomial degree map may occur giving a new CPD map M k+1 . If so, or if convergence is not reached then another iterate is performed.

Fig. 4

 4 Fig.4Chain of detectors used for Linear and Burgers' equation to check if a cell value is valid and accepted, or, if it is discarded and need to be recomputed after reduction of the local cell polynomial degree. ED: Extrema Detector, PD: Plateau Detector, LOD: Local Oscillation Detector, SD: Smoothness Detector.

Fig. 5

 5 Fig. 5 Exact solutions (black) and associate velocities (red) for the linear equation on Ω = [0; 1] -φ(x) = 1/u(x) -Left: regular solution u(x) = 2 sin(2πx) + 3 -Right: irregular solution u(x) = x + 2, if x < 1/2 or u(x) = 50(x -1/2) + 5/2, otherwise.

Fig. 6 Fig. 7

 67 Fig.6Numerical solutions for linear equation and I = 80 cells -Top: unlimited schemes with reconstruction degrees 0 (left), 1 (middle), and 5 (right) -Bottom: a posteriori MOOD limited schemes with cascade P 0 (left), P 1 → P 0 (middle), and P 5 → P 2 → P 1 → P 0 (right). (L and TM schemes provide similar pictures so only L solver's results are shown.)

Fig. 8

 8 Fig. 8 Exact solutions (black) and source term (red) for each of the proposed problem for Burgers' equation -Left: regular solution -Right: irregular solution.

5 Fig. 9 Fig. 10

 5910 Fig.9Simulation of regular solution φ(x) = sin(3πx) exp x + 2 of Burgers equation for the Nonlinear (NL, circle sybols) and Time-Marching (TM, cross symbols) schemes -100 cells and P 5 polynomial reconstructions -Top: numerical solutions. Bottom: evolution of the residuals as a function of the nonlinear solver iterations (NL scheme) or as a function of the time cycles (TM scheme).

Fig. 11

 11 Fig. 11 Simulation of irregular solution of Burgers equation for the Nonlinear (NL) scheme with a posteriori MOOD limiting -120 cells and P 1 → P 0 polynomial reconstructions -Top: numerical solutions for the first 4 iterations. Bottom: last iteration and valid solution.

Fig. 12

 12 Fig.12Simulation of irregular solution of Burgers equation for the Nonlinear (NL) scheme with a posteriori MOOD limiting -120 cells and P 5 → P 2 → P 0 polynomial reconstructions -Top: numerical solutions for 8 iterations. Bottom: last iteration and valid solution.

Fig. 13

 13 Fig.[START_REF] Clain | A high-order finite volume method for systems of conservation lawsmultidimensional optimal order detection (MOOD)[END_REF] Simulation of irregular solution of Burgers equation for the Nonlinear (NL) scheme with a posteriori MOOD limiting -Top: P → P 2 → P 0 polynomial reconstructions -Bottom: P 5 → P 2 → P 1 → P 0 polynomial reconstructions -From left to right: last and valid solutions for 80, 120, 240, and 320 cells.

Fig. 14

 14 Fig.14Simulation of irregular solution of Burgers equation for the Nonlinear (NL) scheme with a posteriori MOOD limiting -Top: cascade P 1 → P 0 , Bottom: cascade P 5 → P 2 → P 0 -From left to right: logscale of the errors for 120, 240, and 320 cells for the first iteration (red), intermediate iterations (green/purple), and the last iteration (blue).

  x) dx, φ = ρ, ρu, E. Moreover, if the solution converges to a steady-state solution in time, denoted by ρ, u, and E, then we have xR xL ρ 0 (x) dx = xR xL ρ(x) dx, xR xL ρ 0 (x)u 0 (x) dx = xR xL D(x) dx, xR xL E 0 (x) dx = xR xL E(x) dx. Proof. Integrating the non-stationary mass conservation equation on the whole domain provides ( ρu) |x=xR -( ρu) |x=xL = D(x R ) -D(x L ). Since the boundary conditions satisfy the algebric relation (15) with no discontinuity, we deduce that ( ρu) |x=xL = D(x L ), ( ρu) |x=xR = D(x R ) and so d dt xR xL ρ dx = 0. Integration in time over the interval [0, t] gives xR xL ρ(x, t) dx = xR xL ρ 0 (x) dx.

  (a) Case D > 0, H > 0. The steady-state solution is constituted of the supercritical solution in [x L , x c [ and the subcritical solution on ]x c , x R ] with a shock located at the unique point x c such that Case D < 0, H < 0. The steady-state solution is constituted of the subcritical solution in [x L , x c [ and the supercritical solution on ]x c , x R ] with a shock located at the unique point x c such that xc xL ρ sub (x) dx + xR xc ρ sup (x) dx = xR xL ρ 0 (x) dx.

Fig. 15

 15 Fig.[START_REF] Dervieux | Steady Euler simulation using unstructured meshes[END_REF] Chain of detectors used for Euler system of equations to check if a cell value is valid and accepted, or, if it is discarded and need to be recomputed after reduction of the local cell polynomial degree. PAD: Physical Admissible Detector, ED: Extrema Detector, PD: Plateau Detector, LOD: Local Oscillation Detector, SD: Smoothness Detector.

Fig. 18

 18 Fig. 18 Numerical solutions for density variablefor the NL scheme a posteriori limited with maximal degree 1 -Left panel: density for I = 25, 50 and 100 cells -Right panel: error for each cell for I = 50, 75, 100 and 200 cells.

Fig. 21

 21 Fig. 21 Numerical errors for density variable and I = 75, 100, and 200 cells for the NL scheme a posteriori limited with maximal degree 0 (red), 1 (blue), and 5 (green). Only the errors at final iteration are shown. The 75 cell results correspond to the right-most curves of Fig. 19.

Fig. 22

 22 Fig. 22 Numerical solutions for density variablefor the NL scheme a posteriori limited with maximal degree 5 with cascade P 5 → P 2 → P 1 → P 0 -Left panel: density for I = 25, 50, and 100 cells -Right panel: error for each cell for I = 50, 75, 100, and 200 cells.

Table 1

 1 Errors and convergence rates for a regular solution of linear equation for the Linear (L) and Time-Marching (TM) schemes. B represents the number of bad cells detected by MOOD loop.

				Linear scheme			Time Marching scheme
		I	E 1	O 1	E∞	O∞	B	E 1	O 1	E∞	O∞	B
		70	1.1E-02	-	2.8E-02	-	0	1.1E-02	-	2.8E-02	-	0
		80	1.0E-02	1.0	2.4E-02	1.0	0	1.0E-02	1.0	2.4E-02	1.0	0
	P 0	90	8.9E-03	1.0	2.2E-02	1.0	0	8.9E-03	1.0	2.2E-02	1.0	0
		100 8.0E-03	1.0	2.0E-02	1.0	0	8.0E-03	1.0	2.0E-02	1.0	0
		200 4.0E-03	1.0	9.8E-03	1.0	0	4.0E-03	1.0	9.8E-03	1.0	0
		70	2.9E-04	-	1.4E-03	-	0	2.9E-04	-	1.4E-03	-	0
		80	2.2E-04	2.1	1.1E-03	2.0	0	2.2E-04	2.1	1.1E-03	2.0	0
	P 1	90	1.7E-04	2.1	8.2E-04	2.1	0	1.7E-04	2.1	8.2E-04	2.1	0
		100 1.4E-04	2.1	6.7E-04	2.0	0	1.4E-04	2.1	6.7E-04	2.0	0
		200 3.4E-05	2.0	1.7E-04	2.0	0	3.4E-05	2.0	1.7E-04	2.0	0
		70	3.7E-07	-	3.2E-06	-	0	3.6E-07	-	3.1E-06	-	0
		80	1.7E-07	6.0	1.4E-06	6.1	0	1.6E-07	5.8	1.4E-06	6.2	0
	P 5	90	8.2E-08	6.1	7.4E-07	5.6	0	8.2E-08	5.9	7.0E-07	5.7	0
		100 4.3E-08	6.0	3.9E-07	6.1	0	4.5E-08	5.7	3.6E-07	6.2	0
		200 6.9E-10	6.0	6.4E-09	5.9	0	8.1E-10	5.8	5.6E-09	6.0	0

Table 2

 2 Errors and convergence rates for a regular solution of Burgers' equation for the Nonlinear (NL) and Time-Marching (TM) schemes. B represents the number of bad cells detected by MOOD loop.

				Nonlinear scheme			Time Marching scheme
		I	E 1	O 1	E∞	O∞	B	E 1	O 1	E∞	O∞	B
		70	7.8E-02	-	1.9E-01	-	0	7.8E-02	-	1.9E-01	-	0
	P 0	80 90	6.8E-02 6.0E-02	1.1 1.0	1.6E-01 1.5E-01	1.0 1.0	0 0	6.8E-02 6.0E-02	1.1 1.0	1.6E-01 1.5E-01	1.0 1.0	0 0
		100 5.4E-02	1.0	1.3E-01	1.0	0	5.4E-02	1.0	1.3E-01	1.0	0
		70	1.8E-03	-	7.4E-03	-	0	1.8E-03	-	7.4E-03	-	0
	P 1	80 90	1.4E-03 1.1E-03	2.0 2.0	5.4E-03 4.1E-03	2.4 2.3	0 0	1.4E-03 1.1E-03	2.0 2.0	5.4E-03 4.1E-03	2.4 2.3	0 0
		100 8.6E-04	2.0	3.2E-03	2.3	0	8.6E-04	2.0	3.2E-03	2.3	0
		70	1.1E-07	-	8.8E-07	-	0	1.8E-07	-	6.2E-07	-	0
	P 5	80 90	4.9E-08 2.4E-08	6.0 6.1	3.8E-07 1.8E-07	6.2 6.2	0 0	1.0E-07 5.7E-08	4.3 4.8	3.6E-07 2.1E-07	4.0 4.7	0 0
		100 1.2E-08	6.2	1.1E-07	5.2	0	3.3E-08	5.1	1.2E-07	5.0	0

  5 → P 2 → P 1 → P 0 . The results are almost equivalent even if some cells are updated with P 1 reconstructions. The number of iterations are respectively 5, 12, 15, and 17.

	1											1												1												1
										Exact										Exact											Exact											Exact
									MOOD iteration 1									MOOD iteration 3									MOOD iteration 3										MOOD iteration 4
										Cells: P5									Cells: P5											Cells: P5											Cells: P5
	0.5											0.5								Cells: P2		0.5									Cells: P2		0.5									Cells: P2
																																	Cells: P0											Cells: P0
	0											0												0												0
	-0.5											-0.5												-0.5												-0.5
	-1											-1												-1												-1
	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	1											1												1												1
										Exact										Exact											Exact											Exact
									MOOD iteration 5									MOOD iteration 7									MOOD iteration 8										MOOD iteration 9
										Cells: P5									Cells: P5											Cells: P5											Cells: P5
	0.5									Cells: P2	0.5								Cells: P2		0.5									Cells: P2		0.5									Cells: P2
										Cells: P0																						Cells: P0											Cells: P0
	0											0												0												0
	-0.5											-0.5												-0.5												-0.5
	-1											-1												-1												-1
	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	-1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
														1
																														Exact
																											MOOD iteration 10
																														Cells: P 5
														0.5																Cells: P 2
																														Cells: P 0
														0
														-0.5
														-1
														-1.5	0		0.1	0.2	0.3		0.4		0.5		0.6		0.7		0.8	0.9		1

Table 3

 3 Errors and convergence rates for the numerical density -Regular solution of Euler system for the Nonlinear (NL) and Time-Marching (TM) schemes. B represents the number of bad cells detected by MOOD loop.

				Nonlinear scheme			Time Marching scheme
		I	E 1	O 1	E∞	O∞	B	E 1	O 1	E∞	O∞	B
		25	7.8E-03	-	1.3E-02	-	0	7.8E-02	-	1.9E-01	-	0
	P 0	50 75	3.9E-03 2.6E-03	1.0 1.0	6.6E-03 4.4E-03	1.0 1.0	0 0	6.8E-02 6.0E-02	1.1 1.0	1.6E-01 1.5E-01	1.0 1.0	0 0
		100 2.0E-03	1.0	3.3E-03	1.0	0	5.4E-02	1.0	1.3E-01	1.0	0
		25	9.0E-05	-	9.3E-04	-	0	1.8E-03	-	7.4E-03	-	0
	P 1	50 75	1.9E-05 7.7E-06	2.3 2.2	2.6E-04 1.2E-04	1.8 1.9	0 0	1.4E-03 1.1E-03	2.0 2.0	5.4E-03 4.1E-03	2.4 2.3	0 0
		100 4.1E-06	2.2	6.9E-05	1.9	0	8.6E-04	2.0	3.2E-03	2.3	0
		25	1.5E-07	-	3.4E-06	-	0	1.8E-07	-	6.2E-07	-	0
	P 5	50 75	2.9E-09 2.4E-10	5.7 6.1	1.3E-07 1.6E-08	4.7 5.2	0 0	1.0E-07 5.7E-08	4.3 4.8	3.6E-07 2.1E-07	4.0 4.7	0 0
		100 4.0E-11	6.3	3.4E-09	5.4	0	3.3E-08	5.1	1.2E-07	5.0	0

  Fig. 17 Numerical solutions for density variable for the NL scheme a posteriori limited with maximal degree 0 -Left panel: density for I = 25, 50, and 100 cells -Right panel: error for each cell for I = 50, 75, 100, and 200 cells.

	4							10 0
			Exact solution							200 cells -P 0
	3.8		P 0 100 cells							100 cells -P 0
			P 0 50 cells							75 cells -P 0
	3.6		P 0 25 cells							50 cells -P 0
								10 -1
	3.4
	3.2
	3							10 -2
	2.8
	2.6
								10 -3
	2.4
	2.2
	2	0	0.2	0.4	0.6	0.8	1	10 -4	0	0.2	0.4	0.6	0.8	1
	4							10 -1
			Exact solution							200 cells -P 1 -P 0
	3.8		P 1 -P 0 100 cells P 1 -P 0 50 cells					10 -2		100 cells -P 1 -P 0 75 cells -P 1 -P 0
	3.6		P 1 -P 0 25 cells							50 cells -P 1 -P 0
	3.4							10 -3
	3.2							10 -4
	3
	2.8							10 -5
	2.6							10 -6
	2.4
								10 -7
	2.2
	2	0	0.2	0.4	0.6	0.8	1	-8 10	0	0.2	0.4	0.6	0.8	1

Acknowledgments

The material of this research has been partly built during the SHARK workshops taking place in Ofir, Portugal, http://www.math.univ-toulouse.fr/SHARK-FV/. The authors would like to acknowledge the financial support of Campus France through the PHC Pessoa labeled 26922YH and entitled "Investigation on very high-order finite volume numerical schemes for fluid hydrodynamics simulation". The research was also supported by the Research Center CMAT of the University of Minho, Portugal, with the Portuguese Funds from the Fundação para a Ciência e a Tecnologia, through the Project PEstOE/MAT/UI0013/2014.

i.e. away from the shock wave. We notice a "characteristics error value" of ∆x k , where k = 1, 2, and 6 for the three schemes, for the error produced in the smooth regions when enough cells are employed. Obvisouly, the maximal error is systematicaly reached in the vicinity of the shock wave, and, there, all schemes behave alike. However, moving away of 5 cells from this region, the high-accurate schemes are performing much better, and, the influence of the presence of a shock follows some sort of a "light cone". The width of this cone does not seem to be related to the number of iterations needed to reach the steady-state solution.