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

Rodolphe Turpault
email: rodolphe.turpault@u-bordeaux.fr

Alexandre Bourriaud
email: alexandre.bourriaud@u-bordeaux.fr

A MOOD-like compact high order Finite Volume scheme with Adaptive Mesh Refinement

Keywords: Finite Volume scheme, High accuracy, Hyperbolic systems, MOOD, AMR

In this paper, a novel Finite Volume (FV) scheme for obtaining high order approximations of solutions of multidimensional hyperbolic systems of conservation laws within an Adaptive Mesh Refinement framework is proposed. It is based on a point-wise polynomial reconstruction that avoids the recalculation of reconstruction stencils and matrices whenever a mesh is refined or coarsened. It also couples both the limiting of the FV scheme and the refinement procedure, taking advantage of the Multi-dimensional Optimal Order Detection (MOOD) detection criteria. The resulting computational procedure is employed to simulate test cases of increasing difficulty using two models of Partial Differential Equations: the Euler system and the radiative M 1 model, thus demonstrating its efficiency.

Introduction

In this work we design a novel multi-dimensional high order explicit Finite Volume (FV) scheme within an Adaptive Mesh Refinement (AMR) context. This method is dedicated to solving general time dependent hyperbolic systems of non-linear conservation laws in multi-dimensions. The Finite Volume method is one of the preferred methods for solving possibly discontinuous solutions of conservation laws [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. It evolves a set of piece-wise constant averaged cell centered values on structured and unstructured meshes. However, because of Lax-Wendroff's theorem [START_REF] Lax | Systems of conservation laws[END_REF], any high order numerical method must be conservative, and, because of Godunov's theorem [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] it must also be non-linear. Several major improvements have been made to the accuracy of the FV method since its introduction. In particular, second order accurate FV methods are nowadays available in many commercial codes. Arbitrary High Order (HO) FV methods have still not gained this popularity, despite their availability in many research laboratories. Some of the main ingredients of an arbitrary HO FV scheme on general unstructured fixed grids are 1. First, an HO time discretization: for instance considering the classical Runge-Kutta time discretization [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods[END_REF][START_REF] Gottlieb | Total Variation Diminishing Runge-Kutta Schemes[END_REF] or even more recent ones like the ADER technique [START_REF] Toro | Towards very high order Godunov schemes[END_REF][START_REF] Titarev | ADER: Arbitrary high order Godunov approach[END_REF]. 2. Second, the HO space discretization: often based on cell-centered polynomial reconstruction by considering a local stencil, for instance using piece-wise linear reconstruction (Monotonic Upstream-centered Scheme for Conservation Laws, MUSCL, [START_REF] Van Leer | Towards the ultimate conservative difference scheme. v. a second-order sequel to Godunov's method[END_REF], Total Variation Diminishing/Bounded TVD, TVB), or higher than second order (Piecewise Parabolic Method PPM [START_REF] Colella | The piecewise parabolic method (ppm) for gas-dynamical simulations[END_REF], or arbitrary higher orders [START_REF] Barth | Higher order solution of the euler equations on unstructured grids using quadratic reconstruction[END_REF] to cite only two). 3. Third, a stabilization/limiting strategy: for instance by blending polynomials in an Essentially Non-Oscillatory (ENO) manner [START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: Analysis and implementation[END_REF][START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF], with WENO (Weighted ENO) probably being the most popular one [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF]; or adapting the notion of slope limiter to polynomials (hierarchical slope limiter [START_REF] Michalak | Accuracy preserving limiter for the high-order accurate solution of the euler equations[END_REF][START_REF] Tsoutsanis | Extended bounds limiter for high-order finite-volume schemes on unstructured meshes[END_REF]); or using some more recent a posteriori Multi-dimensional Optimal Order Detection (MOOD) strategy [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional Optimal Order Detection (MOOD)[END_REF][START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes with high-order MOOD reconstruction[END_REF], to cite only few 1 .

Of the three previous ingredients, the stabilization/limiting is still not entirely agreed upon amongst researchers. Also for HO FV schemes, the main cost is due to the reconstruction/limiting procedure even in a parallel framework, see [START_REF] Dumbser | Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems[END_REF][START_REF] Farmakis | Weno schemes on unstructured meshes using a relaxed a posteriori MOOD limiting approach[END_REF]. Indeed, the polynomial reconstruction requires that one or more stencils of cell neighbors are chosen, to compute several polynomials, and, further, that those polynomials are blended together with non-linear weights. The larger the stencils, the more expensive the reconstruction step is, even if the mesh connectivity remains fixed. Moreover, if some cells happen to change connectivity, then the impacted stencils must be reconstructed as well as the reconstruction matrices. This increases even further the cost of this already expensive procedure. In addition, when steep gradients (shocks, contacts, jumps or high gradients with respect to the space step) develop, then the properties of the numerical solution such as robustness, accuracy, admissibility, stability, are almost entirely dependent on the quality of the limiting strategy. In this work we rely on an a posteriori MOOD strategy, which has produced some interesting results, by changing somewhat the paradigm of limiting HO schemes [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional 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 Multidimensional Optimal Order Detection method in the three-dimensional case: very high-order finite volume method for hyperbolic systems[END_REF][START_REF] Clain | The MOOD method for the non-conservative shallow-water system[END_REF][START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes with high-order MOOD reconstruction[END_REF][START_REF] Coudière | A domain decomposition strategy for a very high-order finite volumes scheme applied to cardiac electrophysiology[END_REF][START_REF] Blanchard | High order accurate conservative remapping scheme on polygonal meshes using a posteriori MOOD limiting[END_REF][START_REF] Nogueira | High-accurate sph method with multidimensional optimal order detection limiting[END_REF][START_REF] Boscheri | A second-order cell-centered lagrangian ADER-MOOD finite volume scheme on multidimensional unstructured meshes for hydrodynamics[END_REF][START_REF] Siengdy | Solution property preserving reconstruction for finite volume scheme: a boundary variation diminishing+multidimensional optimal order detection framework[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] Dumbser | A posteriori subcell limiting of the discontinuous galerkin finite element method for hyperbolic conservation laws[END_REF]. In a MOOD loop, instead of a priori inferring the best-fitted reconstructed polynomials, unlimited HO reconstructions are first employed. Then, locally for the troubled cells, one decreases the polynomial degree of the reconstruction to satisfy the so-called 'Detection criteria'. In the worst case scenario, a first order FV scheme is locally used for those troubled cells, while in others a higher accurate unlimited reconstruction could be employed. The detection criteria are based on physics (positivity, admissibility) and computational (no NaN) requirements, as well as numerical consideration such as ENO behavior. The main originality of a MOOD approach is its a posteriori feature: a candidate solution is computed at the next time step, but if this solution in a cell is observed to be inadequate, then only this cell and its neighbors are sent back for re-computation. Usually we observe that a large majority of the cells does not require any limiting [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]. Among those requiring some intervention, many drop to the first order scheme. This way of limiting produces almost fail-safe schemes, and therefore extremely robust HO FV codes [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] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes with high-order MOOD reconstruction[END_REF][START_REF] Siengdy | Solution property preserving reconstruction for finite volume scheme: a boundary variation diminishing+multidimensional optimal order detection framework[END_REF][START_REF] Farmakis | Weno schemes on unstructured meshes using a relaxed a posteriori MOOD limiting approach[END_REF][START_REF] Semplice | Adaptive-mesh-refinement for hyperbolic systems of conservation laws based on a posteriori stabilized high order polynomial reconstructions[END_REF].

A complementary approach to improving the limiting strategy is to adapt the mesh. For instance, one can gather the cells where the flow presents interesting physical features and, conversely, discard them in flat areas. This behavior can be obtained either by moving the cells towards high gradients (Lagrangian schemes [START_REF] Maire | A high-order cell-centered lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF], Arbitrary-Lagrangian-Eulerian schemes [START_REF] Boscheri | Direct arbitrary-lagrangian-eulerian ADER-MOOD finite volume schemes for multidimensional hyperbolic conservation laws[END_REF], mesh moving techniques [START_REF] Weiming | A moving mesh method based on the geometric conservation law[END_REF], etc.), or, by maintaining a base mesh, but refining and coarsening the cells where appropriate (local AMR [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Semplice | Adaptive-mesh-refinement for hyperbolic systems of conservation laws based on a posteriori stabilized high order polynomial reconstructions[END_REF]). These two techniques may occasionally be coupled together [START_REF] Morrell | A cell by cell anisotropic adaptive mesh ale scheme for the numerical solution of the euler equations[END_REF][START_REF] Waltz | A coupled ALE-AMR method for shock hydrodynamics[END_REF]. The local AMR technique has become popular because it allows saving of a large amount of computational resources when the solution presents only localised features and large zones away from perturbations. As such, AMR FV schemes may solve problems that are intractable on a uniform grid. Moreover with a good and appropriate refinement/coarsening strategy, we can hope to automatically adapt the mesh and the resolution to the underlying unsteady flow. Many AMR libraries are nowadays available including patch-based (Clawpack [START_REF]Clawpack software[END_REF]), cell-based (p4est [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF]), and block-based [START_REF] Keppens | Parallel, grid-adaptive approaches for relativistic hydro and magnetohydrodynamics[END_REF] ones. The AMR procedure generally acts serially after the HO FV scheme has produced a solution. Usually "indicators" and "thresholds" drive the AMR code to refine/coarsen and re-update the solution in some area. The goal is to use a fine mesh on interesting physical features and a coarse one on flat area, thus saving computational resources. These physics-based and smoothness-based indicators are of paramount importance to ensure the quality of the refinement and, as a consequence, the accuracy of the solution.

Coupling an AMR strategy and a HO FV scheme may seem simple at first glance. First the HO FV scheme produces an accurate ENO numerical solution, then the AMR indicators refine/coarsen where appropriate and third, the solution is interpolated/remapped on the new mesh. And the cycle continues. However, making the HO FV scheme work in pair with the AMR strategy is not an easy task. It usually requires tailoring the refinement/coarsening thresholds. Some examples of such a coupling is employed in several codes starting probably from [START_REF] Marsha | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF] and then others like [START_REF] Quirk | A parallel adaptive grid algorithm for computational shock hydrodynamics[END_REF][START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. a new high resolution code called ramses[END_REF][START_REF] Keppens | Adaptive mesh refinement for conservative systems: multi-dimensional efficiency evaluation[END_REF][START_REF] Plewa | Adaptive mesh refinement theory and applications : proceedings of the chicago workshop on adaptive mesh refinement methods[END_REF][START_REF] Baeza | Adaptive mesh refinement techniques for high-order shock capturing schemes for multi-dimensional hydrodynamic simulations[END_REF][START_REF] Van Der Holst | Crash: A block-adaptive-mesh code for radiative shock hydrodynamics-implementation and verification[END_REF][START_REF] Dumbser | Ader-weno finite volume schemes with space-time adaptive mesh refinement[END_REF][START_REF] Buchmüller | Finite volume weno methods for hyperbolic conservation laws on cartesian grids with adaptive mesh refinement[END_REF][START_REF] Semplice | Adaptive-mesh-refinement for hyperbolic systems of conservation laws based on a posteriori stabilized high order polynomial reconstructions[END_REF][START_REF] Freret | High-order finite-volume method with block-based amr for magnetohydrodynamics flows[END_REF] to cite but a few. Stated differently, it is important to design the numerical scheme and the AMR procedure to ensure that they work appropriately together.

In this work we rely on p4est AMR environment [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF] to manage the mesh data-structure and its evolution, and, improve two aspects of an HO FV scheme under AMR framework:

• first, we plan to reduce the need to rebuild reconstruction stencils and matrices in the HO FV scheme when the mesh changes connectivity after the AMR step;

• second, we couple the HO FV MOOD limiting procedure with the AMR refinement/coarsening strategy.

The article is organized as follows. In the second section we briefly present the context of our HO FV scheme within an AMR framework, with a discussion on some weaknesses of the overall approach. The third section presents our solution, consisting of a HO FV scheme with point-wise reconstruction with a stability analysis on advection equation and its embedding into our AMR environment. Numerical experiments are gathered in the fourth section. Two main hyperbolic models are tested, an hydrodynamics model (Euler equations) and a radiative transfer one (M 1 model).

Several tests are chosen to highlight specific properties of the overall HO FV AMR scheme: accuracy, ability to capture isolated and complex interacting waves, ability to appropriately refine the mesh, physical admissibility of the solutions, etc. Finally, conclusions are drawn and various perspectives are discussed.

Context of High Order Finite Volume schemes and Adaptive Mesh Refinement framework

In this section we introduce the hyperbolic systems of conservation laws under consideration. Next the family of HO FV methods of Godunov type is described to setup notation. Then the AMR framework is briefly introduced and the overall scheme is summarized. A final discussion of some weaknesses concludes this section.

Governing equations

2D hyperbolic systems of conservation laws modeling gas dynamics and radiative transfer are considered in this work. They are generically written as

∂U ∂t + ∇ • F(U) = 0, (1)
or

∂U ∂t + ∂F(U) ∂x + ∂G(U) ∂y = 0, ( 2
)
where U is the vector of conservative variables, F(U) = (F(U), G(U)) the associated flux vector, t the time and x = (x, y) the space variable. U and F, G depend on the space and time variable, but we omit this dependency to lighten the notation. The dimension in space is 2 while the dimension of the system is denoted by s. Two models of Partial Differential Equations (PDEs) in 2D are considered in this work. First, the Euler hydrodynamics system of conservation laws, which is a model of dimension s = 4 where

U = (ρ, ρu, ρv, ρE), F(U) = (ρu, ρu 2 + p, ρuv, ρuE + pu), G(U) = (ρv, ρuv, ρv 2 + p, ρvE + pv), (3)
where ρ is the gas density, u = (u, v) t its velocity vector and E the total energy. The system is closed by invoking a complete equation of state (EOS) linking two thermodynamics variables, for instance the density ρ and the specific internal energy ε = E -1 2 u 2 . The EOS computes the pressure p and deduces the sound-speed a as

p ≡ p(ρ, ε) = (γ -1)ρε, a 2 = γ p ρ , ( 4
)
where γ is the ratio of specific heats, a constant which determines the kind of gas considered. This system of PDEs is hyperbolic if the states are admissible and belong to, see [START_REF] Toro | Riemann Sovlers and Numerical Methods for Fluid Dynamics[END_REF],

A = U = (ρ, ρu, ρv, ρE) ∈ R 4 / ρ > 0, ε > 0 . (5)
The second model is the so-called M 1 model of radiative transfer without source term of dimension s = 3 where

U = (E R , F R,x , F R,y ), F(U) = (F R,x , c 2 P R,xx E R , c 2 P R,xy E R ), G(U) = (F R,y , c 2 P R,yx E R , c 2 P R,yy E R ), (6)
where c 2.998 × 10 8 m.s -1 is the speed of the light in vacuum. E R is the radiative energy, while

F R = (F R,x , F R,y ) t is the radiative flux. The anisotropy factor, f = F R cE R
, is an important derived quantity. The system is further closed by an entropic closure relationship, which leads to the following expression of the radiative pressure tensor

P R = 1 -χ 2 I + 3χ -1 2 F R ⊗ F R F R 2 = P R,xx P R,xy P R,yx P R,yy , with χ ≡ χ ( f ) = 3 + 4 f 2 5 + 2 4 -3 f 2 ,
where I is the identity matrix. This system of conservation laws is also hyperbolic provided that U ∈ A, see [START_REF] Dubroca | Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif[END_REF] A

= U = (E R , F R,x , F R,y ) ∈ R 3 / E R > 0, F R < cE R . (7)
We can remark that the last physical criteria in A can be recast in terms of the anisotropy factor as f < 1. Moreover, the M 1 model is a stiff system compared to Euler hydrodynamics system of conservation laws. These systems are further equipped with some initial data U(t = 0) ≡ U 0 and boundary conditions.

Finite volume discretization

The fixed computational domain Ω ∈ R 2 is a polygon meshed with cells with neither gaps nor holes.

Mesh

Recall that we want to employ AMR technique later on. As a consequence, each cell Ω c is square, of size ∆x × ∆y and is referred to with a single index 1 ≤ c ≤ M. The mesh has M cells and M may evolve in time. The sizes ∆x and ∆y may be different from cell to cell, and, rigorously should be denoted ∆x c and ∆y c , but we omit the c index whenever possible. The cell center is denoted (x c , y c ) and the cell interfaces are horizontal and vertical segments. Let us introduce the notation

x ± c = x c ± ∆x 2 and y ± c = y c ± ∆x 2 , as such Ω c = [x - c ; x + c ] × [y - c ; y + c ].
The cell area (also referred to as its 'volume' by an abuse of semantics) is |Ω c | = ∆x 2 . Although each cell has a square shape, due to the AMR nature of the mesh, it may be a degenerate polygon. In figure 1 denoted by ∂Ω c consists of 4 to 8 vertical and horizontal edges. Each edge of cell c is common to a 'generic' neighbor cell d and this edge is denoted e cd of length cd . The unit outward pointing normal to ∂Ω c is n c . The set of neighbor cells sharing an edge with cell c is N(c), while the set of its neighbors is V(c). The time domain T = [0, T ] is split into time steps [t n , t n+1 ], of size ∆t = t n+1t n subject to a CFL (Courant-Friedrichs-Lewy) like condition.

Numerical flux and discretization

In a FV context, the data within a cell is a constant approximating the mean value of the solution U(x, t n ) over the cell at time t n , that is

U n c 1 |Ω c | Ω c U(x, t n ) dv = 1 ∆x 2 x + c x - c y + c y - c U(x, y, t n ) dx dy. (8)
The vector of all mean values at time t n is denoted as

U n = U n c c=1,...,M . (9)
Integrating (1) in space over Ω c and time over [t n , t n+1 ], knowing U n c for all cell c yields

t n+1 t n Ω c ∂U ∂t + ∇ • F(U) dv dt = 0, (10)
which, given the definition (8), and considering an explicit time discretization, becomes

|Ω c | U n+1 c -U n c + ∆t ∂Ω c F(U n (s)) • n c ds = 0. ( 11
)
Next the boundary ∂Ω c is split into its horizontal and vertical edges e cd shared by cells c and d, and

|Ω c | U n+1 c -U n c + ∆t d∈N(c) e cd F(U n (s)) • n cd ds = 0. ( 12
)
Designing a FV numerical scheme consists of computing the above integral accurately. At cell interface the flux may not be uniquely defined, i.e in general U n c U n d , therefore F(U n c ) F(U n d ). Therefore we introduce the numerical fluxes F and G in the horizontal and vertical directions respectively. Then, along cell edges e cd , we choose a 1D quadrature formula made of G > 0 quadrature points and

e cd F(U n (s)) • n cd ds G g=1 cd w g F U n c , U n d ; (x g ) , (13)
where e cd is a vertical or horizontal edge of length cd between cell c and d, x g = (x g , y g ) is a quadrature point along it, and w g is its associated weight. The numerical flux F depends on the reconstructed polynomials U n c in cell c evaluated at point x g and its counterpart U n d in cell d. The reconstruction operators are described in the next section. Substituting ( 13) into [START_REF] Boscheri | Direct arbitrary-lagrangian-eulerian ADER-MOOD finite volume schemes for multidimensional hyperbolic conservation laws[END_REF], we write the generic FV scheme for any cell c as

U n+1 c = U n c -∆t d∈N(c) cd |Ω c | G g=1 w g F U n c , U n d ; x g . (14)
In this work the Rusanov, HLL and HLLC numerical fluxes are considered [START_REF] Toro | Riemann Sovlers and Numerical Methods for Fluid Dynamics[END_REF]. For instance, the Rusanov flux at any point z on edge e cd is generically of the form

F U n c , U n d ; z = F U n c (z) + F U n d (z) 2 • n cd - |Λ cd | 2 U n d (z) -U n c (z) • n cd . ( 15
)
The first term is a centered physical flux while the second one is the dissipative term scaled by factor |Λ cd |, which must be larger than the largest eigenvalue of the flux's Jacobian. For the Euler system we pick

|Λ cd | = max(|u c ± c c |, |u d ± c d |)
and for the M 1 model we set |Λ cd | = c . Both terms employ polynomial extrapolated values U n c (z) and U n c (z) at quadrature point z. The simplest first-order accurate scheme does not reconstruct any polynomial, that is U n c (x) = U n c for all point x ∈ Ω c , and considers one quadrature point located at the edge center. Notice that the scheme ( 14) could be written in a somewhat simpler form by using the fact that the normal is either (1, 0) or (0, 1). However, the unstructured point of view in [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF] allows for more generality when non-conformal AMR meshes with more than four edges are considered (see figure 1). Finally, let us recall that such an explicit Euler scheme is stable under the following CFL constraint

∆t expE ≤ CFL ∆x 4 max cd |Λ cd | , (16)
where CFL number is a user-given constant usually 0 < CFL < 1.

Time discretization

The high order time discretization employs SSP-RK method (Strong Stability Preserving Runge-Kutta) [START_REF] Gottlieb | Total Variation Diminishing Runge-Kutta Schemes[END_REF]. The same accuracy is chosen for time and space approximations. Because SSP schemes preserve the convex sets, the limiting can essentially be undertaken in the space discretization, as outlined in the following sections. Here our third order SSP-RK scheme is the SSP-RK(3,3) scheme

U (1) h = U n h + ∆t L(U n h ), (17)
U (2) h = U (1) h + ∆t L(U (1) h ) → U * h = 3 4 U n h + 1 4 U (2) h , (18)
U (3) h = U (2) h + ∆t L(U * h ) → U n+1 h = 1 3 U n h + 2 3 U (3) h . ( 19
)
where L is our spatial FV discrete operator. For 2nd and 4th orders respectively, we rely on the SSP-RK(2,2) and SSP-RK [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF] schemes from [START_REF] Gottlieb | Total Variation Diminishing Runge-Kutta Schemes[END_REF][START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF]. These schemes are stable provided that

∆t ≤ κs∆t expE , ( 20
)
where s is the number of stages and κ ∈ (0, 1) is the efficiency coefficient of the SSP-RK scheme.

Polynomial reconstruction from cell mean values

We have seen that the accuracy of flux ( 15) is linked to the quadrature formula, which we always assume to be sufficiently accurate, and, more importantly, to the polynomial extrapolated states U n c (z), U n d (z). In order to compute these extrapolated states seen in ( 14), a polynomial is reconstructed for each primitive/conservative variable and cell at time t n . Notice that a reconstruction based on characteristic variables may improve the results for simple models. However, for more complex ones, these variables are expensive to compute [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] Miyoshi | A short note on reconstruction variables in shock capturing schemes for magnetohydrodynamics[END_REF], hence our choice to rely on conservative/primitive ones. We search for a polynomial of degree d > 0 in each cell denoted

U c ∈ P d (Ω c ) such that 1 |Ω c | Ω c U c (x) dv = U c , (21)
and U c best fits the mean values of neighbor cells in the least-square sense. The stencil of neighbor cells is denoted by S d c . In other words, the polynomial coefficients of U c ∈ P d (Ω c ) are computed by minimizing the cost function

J( U c ) = 1 2 d∈S d c 1 |Ω d | Ω d U c (x) dv -U d 2 . ( 22
)
The stencil S d c should contain at least N = (d+1)(d+2) 2 cells to solve ( 21)-( 22) for degree d > 0. In practice, more cells are considered, usually about 3 2 N is employed [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional Optimal Order Detection (MOOD)[END_REF]. The same degree is used for all physical variables in a cell, and the quadrature formula is always assumed to be of at least the same order of accuracy than d. As already mentioned, the reconstruction for d = 0 consists of U c (x) = U c for all x in Ω c . Once the reconstructed polynomial U c in cell c of degree d > 0 is computed, it is evaluated for each quadrature point z to compute [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods[END_REF] and then deduce [START_REF] Buchmüller | Finite volume weno methods for hyperbolic conservation laws on cartesian grids with adaptive mesh refinement[END_REF]. In practice solving the least-square reconstruction problem consists in solving a linear system of size d × d. The matrix depends on the location of each neighboring cell Ω d , d ∈ S d c relative to the current one. Therefore the choice of the neighborhood S d c is crucial. For instance, in figure 2-c we present the 13 cells neighborhood constructed by a two-layer face-based neighborhood. However, this choice may lead to a geometrical bias in direction of large cells. More importantly, if the mesh changes connectivity then these stencils have to be recomputed. In fact, the stencil selection is still a somewhat subjective procedure for several polynomial based numerical schemes, like WENO, but also to determine the minimal/maximal bounds used within classical slope limiters. In summary, the classical polynomial reconstruction operator requires: a cell c, a polynomial degree d > 0, an appropriately large enough set of neighbor mean values

A posteriori MOOD limiting strategy

Any high order scheme is subject to the Gibbs phenomenon, illustrated by the creation of spurious oscillations in the vicinity of steep gradients. In this work we rely on a MOOD loop [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional 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 Multidimensional Optimal Order Detection method in the three-dimensional case: very high-order finite volume method for hyperbolic systems[END_REF] to stabilize the solution. The MOOD paradigm assesses a posteriori the admissibility of a candidate numerical solution at time t n+1 . If a cell does not fulfill some physical and numerical validity criteria, it is recomputed with a more robust scheme, as illustrated in figure 3.

Cascade and parachute scheme. In order to design a MOOD loop the user must define an ordered list of schemes. This ordered list is called the 'cascade'. At a minimum the cascade must be composed of two schemes: the most accurate, often oscillatory scheme with reconstructions of degree d max > 0, and, a dissipative, robust and fail-safe one called the 'parachute scheme', in our case with d = 0. The unlimited high accurate scheme of (d max + 1)th order is the first one to be tried, while the first-order accurate one is the last one, that is the parachute scheme (d = 0). Some intermediate schemes could be tried in-between. In this work we will employ the following cascade:

d max = 2 → 1 → 0 although different cascades such as d max = 2 → 0, or d max = 2 → 1 + slope limiter → 0 would be possible.
Detection criteria: Computer, Physical and Numerical admissibility. To assess if a numerical solution is accepted, we design so called Detection criteria, see figure 3. If a candidate solution U n+1, * i in cell Ω i does not pass the criteria, it is sent back at t n for re-computation following [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional 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 Multidimensional Optimal Order Detection method in the three-dimensional case: very high-order finite volume method for hyperbolic systems[END_REF]. We design three sets of detection criteria Computational -CAD. This criterion is employed to check for floating point exceptions, in other words the numerical solution is assessed for the Computer Admissible Detection criteria if it does not present any NaN (Not-a-Number), Inf (Infinite), etc. values.

Physical -PAD. The numerical solution in a given cell is assessed against the PAD criterion if it remains in the set of admissible states A. The PAD criteria is strongly related to the PDE model being solved.

Numerical -NAD. The numerical solution in a given cell Ω i is numerically admissible if it presents an ENO behavior.

Here we rely on a Relaxed Discrete Maximum Principle (RDMP) on the conservative variables, which writes for any variable A in cell number i as

-δ + m d i ≤ A n+1, * i ≤ M d i + δ, (23)
where

m d i = min k∈S d i (A n k ), M d i = max k∈S d i (A n k ) and A n k is the mean value of variable A in cell Ω k ∈ S d i at time t n . The so-called relaxation parameter δ is δ = min 10 -4 , 10 -3 |M d i -m d i | .
Note that if δ = 0, then a strict DMP is forced for which a 2nd order accuracy bound would be observed in L ∞ norm.

These criteria are of paramount importance for the detection of invalid cells, and, as such, should not be too permissive nor too preventive.

MOOD loop. As seen in figure 3 the MOOD loop includes the FV scheme solver. At t n the maximal degree d max is set. The d max + 1 order FV scheme is used to produce a candidate solution. Then the detection procedure splits the admissible cells from the problematic/troubled ones. These troubled cells and their direct neighbors are then re-calculated by the next and more robust scheme in the cascade. Once the next scheme in the cascade produces a new candidate solution, it is then processed by the detection criteria. The new troubled cells are further corrected in another iteration of the MOOD loop using an even more robust scheme. When the parachute scheme is reached in the cascade, then the candidate solution is declared as admissible. The MOOD loop is rather inexpensive if few cells require correction, and, more importantly, it always converges. The sketch from figure 3 represents the whole high accurate a posteriori MOOD limited numerical scheme.

General context for Adaptive Mesh Refinement -Quadtree, indicators and projection

The AMR technique consists of refining the mesh to maintain important physical gradients in the finest cell level, and, also coarsening it in flat/unimportant areas.

In the AMR context of p4est [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF], the meshes are stored as quadtrees. An illustration is presented in figure 2(b). The level of refinement of a cell corresponds to its level in the tree. Each cell is a square which is split into 4 smaller squares when refined. When four neighbor cells, emanating from the same mother cell, are identified for coarsening, they are merged back onto the mother cell. Usually a minimal and maximal level denoted l ≤ L are given by the user to avoid extreme refinement or coarsening. As such, in 2D, the coarsest mesh has a minimum of (2 l ) 2 cells, while the finest possible one is (2 L ) 2 .

In order to refine or coarsen a cell, for classical cell-based AMR procedure, two indicators are in general defined.

Given the solution on the previous AMR mesh, a projection operator computes/interpolates/remaps the solution on these newly refined/coarsened cells.

Refining/Coarsening indicators. We must determine locally computable and consistent refinement and coarsening indicators R c (U) and C c (U) for cell c, and their associated thresholds R and C, so that

If R c ≥ R =⇒ Cell c requires refinement. If C c ≤ C =⇒ Cell c requires coarsening.
Usually the gradient of a variable (density, entropy, etc.) is used with ad-hoc thresholds. There is no agreement on an universal indicator, and, it is not the purpose of this work to optimize this choice. Here, instead of designing R c , we propose to take advantage of the detection criterion used by the MOOD loop to define an innovative refinement strategy. The idea is to use the MOOD detection criterion to detect a troubled cell which will first be refined. If the daughter cells are still troubled then limiting will occur. Only the coarsening strategy remains classical, that is the threshold C and diagnostic C c need to be chosen. To do so, we pick a classical and simple indicator related to the gradient. The default setting is to consider the gradient of the first variable. Typically, in the case of Euler equations, we therefore choose

C c = ∇ρ 2 . (24)
Note that a cell is coarsened only if the three other sister cells also agree for coarsening.

Projection/Remapping operator. When a cell c is refined, it is split into four children square cells onto which the solution must be projected. First one computes the reconstruction polynomial of degree d = d max , U n c ∈ P d in cell c. Recall that d max + 1 is the maximal order of accuracy fixed by the user at the beginning of the computation. Then, for any child cell c k , k = 1, 2, 3, 4, such that 4

k=1 Ω c k = Ω c , one defines

∀1 ≤ k ≤ 4, U n c k = 1 |Ω c k | Ω c k U n c (x) dv. ( 25
)
The conservation is respected by construction because

4 k=1 |Ω c k |U n c k = 4 k=1 Ω c k U n c (x) dv = Ω c U n c (x) dv = U n c , (26)
the last equality corresponding to equation [START_REF] Diot | Improved detection criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]. Notice that the projection in ( 25) is accurate but may lead to nonphysical values, for instance negative density. In this unlikely situation, the "reconstruction" is restricted to U n c = U n c ∈ P 0 . We do not expect such situation to be encountered often as the operations of refining and coarsening are supposed to occur in smooth regions. When four child cells c k , k = 1, 2, 3, 4, are coarsened into their parent cell indexed by m, Ω m = 4

k=1 Ω c k , then the new conservative mean value is computed as

U n m = 1 |Ω m | 4 k=1 |Ω c k |U n c k . (27)
In this case no reconstruction is needed, and we notice that coarsening four cells, inexorably leads to some loss of information. As such, coarsening must be done with care.

This concludes our ultra-simplified AMR context, recalling that more subtle and clever procedures could be designed to better adapt the mesh to the flow.

Summary

As a summary we can state that the FV numerical method presented in the previous section enjoys 1. an arbitrary order of accuracy in time thanks to SSP-RK discretization with a well defined CFL condition; 2. an arbitrary order of accuracy in space thanks to the polynomial reconstructions and the high order quadrature formula for the integration along the edges; 3. an a posteriori MOOD limiting strategy to damp spurious oscillations and ensure physical admissibility; 4. an AMR procedure to refine/coarsen the mesh where appropriate according to ad-hoc indicators.

As already mentioned, due to the use of an AMR procedure the reconstruction stencils and matrices in step 2. need to be re-updated when a cell changes its neighbors or is created by the AMR procedure.

Discussion: revamping stencils and reconstruction matrices

Inexorably the AMR procedure creates new cells and discards some others. Therefore the stencils S d i for all cells i in the vicinity of a new or discarded cell must be reconstructed. Usually stencils in transition from coarse to fine levels are not uniform, which does not help the general stability and the conditioning of the associated linear systems.

In figure 2 we present some examples of such non-uniformity. In fact, for a second-order scheme, refining one cell impacts all its neighbours' stencils: from 4 to 12 stencils may need modification, 3 new ones must be designed from scratch. For a third or fourth order reconstructions, the stencils are wider, so that the number of stencils to be revamped highly grows. On the other hand, coarsening four cells also leads to equivalent stencil reshaping. Accordingly, the modification of the stencils implies that the associated reconstruction matrices must be re-assembled and inverted. Designing the stencils and their associated reconstruction matrices is one of the main costs for a large number of high order FV schemes on fixed grid, especially when truly high orders of accuracy are targeted, i.e. beyond 3rd order. On an AMR grid, this cost is even higher because the mesh adapts dynamically to the flow without any possibility to anticipate the evolution of the stencils 2Fortunately, this extra work could be drastically reduced if the polynomial reconstruction matrix would act on the same "stencil of points" independently of the underlying mesh configuration. This is the goal of the next section: to design such a reconstruction procedure acting on a stencil of fixed point values rather than on a stencil of cells. The primary goal of this work is to design a new point-wise FV scheme which is agnostic to cell neighborhood changes, avoiding de facto the redefinition of stencils for the polynomial reconstruction within an AMR configuration.

High Order Finite Volume schemes with point-wise reconstruction (HOp)

General principles

The main idea is motivated by the observation that the polynomial reconstruction procedure can be designed from highly accurate point values instead of cell-centered mean values, unlike in [START_REF] Diot | The Multidimensional Optimal Order Detection method in the three-dimensional case: very high-order finite volume method for hyperbolic systems[END_REF]. The local polynomial reconstruction problem is then stated as finding U c ∈ P d (Ω c ) defined in cell c such that ( 21) is fulfilled, and, U c fits at best P >

J( U c ) = 1 2 p∈P d c w p U c ( xp ) -U p 2 , (28)
under the constraint

1 |Ω c | Ω c U n c (x) dv = U c . (29)
Notice that in [START_REF] Freret | High-order finite-volume method with block-based amr for magnetohydrodynamics flows[END_REF] one evaluates the unknown polynomial at point xp , which we refer to as an "evaluation point", instead of computing the mean value over neighbor cells. With such a procedure, one must choose the number of evaluation points P ≥ 1 2 (d + 1)(d + 2), their positions xp and weights w p and values U p to be matched with sufficient accuracy. Then, one searches for the reconstruction polynomial U c of degree d (i.e (d + 1)th order accurate) which has the same mean value as U n c in cell c, and matches at best the point values ( xp , U p ) p∈P d c . The weights are intended to balance the relative importance of each point. The values U p should be accurate enough for the final scheme to have the expected order, and must be constructed elsewhere. The procedure is referred to as the P d,h point-wise reconstruction procedure: d is the degree of the searched polynomial U c , while h + 1 is the order of accuracy of the point-wise data. In figure 4 main advantage of this approach is the fact that the evaluations points are chosen at the beginning of the simulation for each degree and fixed during the entire AMR simulation. Any cell is always a square (homothetic to the unit square), therefore only one reconstruction matrix per set of evaluation points must be stored. In addition, these matrices are assembled once, at initialisation. Stated differently, this approach allows the AMR calculation to be in the same reconstruction environment as when the mesh is fixed.

Evaluation point positions and number. The number of points is chosen close but superior to the optimal choice, which is 1 2 (d + 1)(d + 2) for a degree d in 2D. Moreover, it may be wise to choose a multi-symmetric distribution, since all cells are squares of size ∆x × ∆x. For simplicity, we distribute the points on the boundary of the homothetic square of size (∆x + ξ) × (∆x + ξ) with 0 ≤ ξ ≤ 1 2 ∆x. The first set of N = 8 points consists of the four points located at corners (points 1,5,9,13) and the mid-edge ones (points 3,7,11,15), for degrees less than d = 2, as shown in figure 4-left. For degrees d ≤ 3 we may choose the distribution made of N = 16 points, as in figure 4-middle. In our work the choice of the distance ξ is set to ξ = ∆x 2 to ensure that the evaluation points remain inside a neighboring cell. Finally, the cell into which an evaluation point is located is identified by two indices c and p, such as Ωp,c . Notice that Ωp,c is necessarily a neighbor cell of Ω c (node-wise). This feature simplifies the points localization's procedure since only the knowledge of the direct neighbors and their size are required. An evaluation point located on an edge is associated with the two neighbouring cells, for instance see point 2 in figure 4-middle.

In the present study, we restrict ourselves to polynomials of degree 3, and thus to schemes of order 4. Going beyond 4th order would probably require increasing the number of evaluation points, but, even then, it is not clear that there would be a suitable set of weights, nor that the benefit in an AMR code would be cost effective. The extension of this method to order higher than 4 is therefore an open question.

Evaluation point values. The point values U p to be matched should naturally be h ≥ d accurate values. Due to symmetries in the points' distributions, we expect that h = d -1 is enough. The numerical simulations in section 4 support this statement. For the sake of completeness, both situations are analyzed in the following. The simplest first-order accurate version consists of setting the point value to the mean cell centered value, so that the 'associated reconstruction' is first-order accurate with h = 0:

U p = U d ,
where

Ω d = Ωp,c , (30)
leading to a procedure P d,0 . For higher accuracy one can rely on local polynomial reconstruction U c (x) in cell c and set

U p = U c ( xp ), where U c ∈ P h (Ω c ). ( 31
)
The reconstruction U c usually involves only the direct neighbors of Ω c , which are easily accessible. If an evaluation point lies on the edge between two cells, says, c 1 and c 2 , the point-wise value is set to be the average:

U p = 1 2 U c 1 ( xp ) + U c 2 ( xp )
, for instance the mid-edge points 3,7,11,15 in figure 4.

Evaluation of point weights. The stability region of this reconstruction procedure can be increased by correctly determining the weights. In practice the weights may modify the matrix conditioning, but do not change the matrix structure. In the next subsections we propose a stability analysis to estimate some appropriate weight values for two second-order procedures P 1,0 and P 1,1 , two third-order ones P 2,1 and P 2,2 and a fourth-order one P 3,2 .

Stability analysis of High Order Finite Volumes point-wise reconstruction schemes for advection equation

The stability analysis is conducted on the 2D advection equation on a uniform 3×3 grid centered at the origin, meaning that the current cell center is x c = (0, 0). The expansion of the polynomial writes for any degree d > 0:

U c (x, y) = α 00 + d i, j=0 i+ j≤d α i j x i y j -I i j ,
where

I i j = 1 |Ω c | Ω c x i y j dv, (32)
where α = (α i j ) t i, j is the vector of the unknown polynomial coefficients of size

K = (d+1)(d+2) 2 in 2D. Matching U c (x)
on the P evaluation points ( xp , w p , U p ) p∈P d c leads to the following linear system:

                                           w 1         α 00 + i, j α i j xi 1 ŷ j 1 -I i j         = w 1 U 1 , w 2         α 00 + i, j α i j xi 2 ŷ j 2 -I i j         = w 2 U 2 , . . . . . . . . . w P         α 00 + i, j α i j xi P ŷ j -I i j         = w P U P , (33)
which in matrix form yields

W M α = W Û, (34)
where

M =                    1
                   , W = diag(w), and U = U 1 , U 2 , • • • U P t , (35)
where diag(w) is a diagonal matrix, the terms of which are the components of w = (w 1 , . . . , w P ) t . The linear system is a priori rectangular with more points than unknown coefficients: P > K. One solves

M t W M α = M t W Û =⇒ α = M t W M -1 M t W Û. ( 36
)
We carry out a stability analysis in order to choose the weights and point positions. The linear stability study uses these reconstructions in an upwind FV scheme solving the 2D advection equation at constant velocity (1, 1). Then a von Neumann stability analysis is performed. To ease the study, the weights are split into three families: centered weights arbitrarily set to w c = 1 (points 3,7,11,15 in figure 4), corners weights w cn (points 1,5,9,13), and mid-edge weights w e > 0 (points 2,4,6,8,10,12,14,16). As such, only two parameters remain: w cn > 0 and w e > 0. We resort to numerical sampling to check under which ∆t constraint all Fourier modes are damped. Therefore there exists a stability region bounded by a maximal CFL number. Some results for varying weights are summarized in Table 1. It can be seen that the influence of the weights on the stability region is important for the second order schemes. If the sensitivity is high for the P 1,1 scheme, some values make it possible to reach a CFL value close to one. On the other hand, it does not seem possible to achieve this with the P 1,0 scheme, which is much less stable in general. Therefore, the latter is not used in practice, and reference will be made to P 1,1 when referring to a second order scheme. On the other hand, any reasonable choice of weights gives close to optimal stability for orders 3 and 4.

Since these schemes are to be integrated in a MOOD structure, and for the sake of simplicity in the programming, the same set of weights is chosen in the following to obtain the highest possible stability for all orders: (10, 5, 1). Naturally, it would be possible to choose two different sets of weights, but we have chosen the same ones for both reconstructions. In addition to simplifying the stability study, it allows us to use the coefficients previously obtained in the computation of the cell polynomial to compute the point values at the next time step.

High order FV numerical scheme with point-wise reconstruction and integration in an AMR framework

The previous sections have introduced a FV numerical scheme for which the reconstructed states at the quadrature points used to estimate the flux functions are given at 'fixed' points independent of the shape and the number of surrounding neighbor cells. Recall that the cell neighborhood often changes in an AMR simulation. The goal of our approaches is to avoid exploring each modification of neighborhood to reconstruct the cell stencils and matrices, as it is classically done for cell-based polynomial reconstructions. The HOp numerical scheme can be written as [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF], except that the polynomial reconstructions U n c are replaced by point values U c ( x p ) at fixed positions x p . This procedure is called 'point-wise polynomial reconstruction'. We adopt the 16 points depicted in figure 4 along with the associated weights 10, 5 and 1 obtained from the stability analysis in section 3.2. This set of points allows computation of second, third and fourth order accurate polynomial reconstructions, hence d ≤ 3. The rest of the numerical method remain the same. At the beginning of any time-step the cell mean values in a given mesh are required and also the point-wise value at each evaluation points. We employ the initial condition to set up the cell mean values and the point-wise ones for n = 0. Within an AMR environment and considering the MOOD limiting strategy, we perform the updates as follows.

• We take advantage of the MOOD criterion to detect the cells that must be recomputed with a lower order, and those that must be refined. First, an unlimited candidate solution is computed and checked against the detection criteria in order to identify the troubled cells. These troubled cells are first refined, and, an unlimited candidate solution is re-computed. We then enter the classical MOOD loop illustrated in figure 3. This use of the MOOD criterion avoids the tailoring of any refinement parameters.

• In order to simplify the computational procedure, the update of the point values is performed at the end of the loop using the coefficients already computed for the polynomial reconstruction in the cells. To put it more clearly, in each cell the point values necessary for the next iteration are computed by evaluating the polynomial as in [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] where α 00 = U s+1 c and the other coefficients α i j of the cell reconstruction at Runge-Kutta stage s are reused to avoid additional computations. These coefficients are thus shifted by one stage (i.e. one stage of the Runge-Kutta method). This shift does not affect the order of the method as will be shown in the numerical results. A sketch of our AMR strategy for the current HO FV scheme based on MOOD limiting is depicted in figure 5. As expected, after the troubled cells have been refined, we retrieve the classical FV MOOD based scheme, as in figure 3.

The blue boxes, 'Refine' and 'Coarsen', include the computation of the new values in the refined/coarsen cells. As usual in a MOOD-like strategy, if no cell is detected as troubled,then the solution at t n+1 is unlimited with dth order of accuracy. On the other hand, if some cells are troubled, they are first refined and, then, recomputed. The second detection pass states if they are accepted or sent back to be recomputed with a lower order scheme in the cascade.

To complete this description of the computational procedure, we mention three practical points:

• The coarsening is not necessarily performed at each update, but only after a given number of time steps (3 by default in the numerical simulations presented later). It seems indeed unreasonable to coarsen after each stage of the Runge-Kutta method. Furthermore, the objective in this article is not to obtain an optimal coarsened mesh but to always ensure an accurate approximation.

• The reconstructions for cells on the boundaries are computed using the same stencils as inside the domain without taking into account exterior "ghost" points. The weights are of course different, and, it is therefore necessary to pre-compute an adapted reconstruction for the boundary cells. The approach remains the same as inside the domain except that there are not enough points left to obtain a 4th order reconstruction. This is why we restrict ourselves to 3rd order accurate schemes in the numerical section.

• The choice to refine before lowering the order of the scheme is made for the sake of consistency and in the spirit of trying to obtain the most accurate approximation possible. It would indeed be possible to refine after the MOOD loop though in this case another detection criterion should be designed.

This concludes the description of the numerical simulation. It is validated on several numerical experiments in the following section.

Numerical experiments

In this paper we restrict the maximal polynomial degree to d = 2. Several key points are tested:

• the point-wise reconstruction accuracy on smooth test case;

• the FV scheme with point-wise reconstruction on fixed mesh with limiting for Euler and M 1 models;

• the coupling of the HO FV scheme with point-wise reconstruction in the AMR environment. Recall that our strategy consists of refining the troubled cells and coarsening only at the end of the time-step.

In order to test these key points we propose several appropriate problems. First, the order of convergence and the accuracy of the second-order (P 1,1 ) and third-order (P 2,2 ) schemes are monitored for the Euler model employing the vortex test case on fixed Cartesian and AMR meshes. Being regular, this problem does not need any limiting. We also gives the results obtained by a classical cell-centered reconstruction for comparison purposes. Next, the results of the FV scheme under MOOD limiting with point-wise reconstruction on a fixed Cartesian mesh are presented for the M 1 model on a filamenting test case. Then, the behavior of the complete AMR simulation code is tested on the cylindrical Sod shock tube to show the impact of the refinement/coarsening strategy coupled with the limiting. Next, on two 2D Riemann problems with 4 shocks and 4 contact discontinuities, we pursue the testing campaign of the AMR code on Euler model. In the latter case, the influence of the selection of the first order solver used as a basis for the method is investigated by comparing the results achieved with HLL and HLLC Riemann solvers. Finally, a 2D Riemann problem with highly twisted and distorted behavior for the M 1 model is considered. This test reveals the flexibility of the overall AMR method to complex and intricate situations.

In the following, the AMR mesh is defined with the two user-specified parameters (l, L): the coarsest possible level of refinement l and the finest ones L ≥ l. The number of allowed refinement steps is therefore N r = Ll + 1, and, if l = L then the mesh can not be refined nor coarsened. On the unit square, the (l, L) AMR mesh corresponds to cells of minimal and maximal lengths in the range 1/2 l ≤ ∆x = ∆y ≤ 1/2 L . For instance the parameters (l, L) = (4, 7) allows N r = 4 levels of refinements with the coarsest cell length ∆x = 2 -4 = 1/16, and finest one ∆x = 2 -7 = 1/128.

Euler model: vortex test case

This first test case involves the Euler system of conservation laws, with the goal of checking the numerical accuracy and orders of convergence of the method. It consists of the transport of a vortex in the computational domain Ω = [-5, 5] × [-5, 5] with periodic boundary conditions. The following initial condition is used

H 0 (x, y) = (ρ 0 (x, y), u 0 (x, y), v 0 (x, y), p 0 (x, y)) t =                        (1 -25(γ-1)e 1-r 2 8γπ 2 ) 1 γ-1 1 -5ye 0.5(1-r 2 ) 2π 1 + 5xe 0.5(1-r 2 ) 2π (1 -25(γ-1)e 1-r 2 8γπ 2 ) γ γ-1                        , r = x 2 + y 2 , γ = 1.4. ( 37
)
Here, the simulations are performed on a hierarchy of fixed meshes until the final time t final = 5 × 10 -2 is reached. The numerical approximations obtained with both the point-wise and classical cell-wise reconstructions are then compared with the exact solution. The errors are gathered in table 2, where we can observe that the expected 2nd and 3rd orders of convergence are attained for the cell-wise and point-wise schemes. The errors given by the two types of reconstructions are comparable, but marginally better for the classical reconstruction on coarse meshes. Notice that for the point-wise reconstruction the correct orders of convergence are also obtained with different weights, even if the errors may be slightly different (not shown here).

Classical cell-wise reconstruction

Point-wise reconstruction Expected order → 3 3

N c L 1 error L 2 error L 1 order L 2 order L 1 error L 2 error L 1 order L
Table 2: L 1 and L 2 errors and convergence rates for the unlimited HO Finite Volume scheme equipped either with the classical cell-wise reconstruction or the current point-wise one for second and third order versions.

Next, we present the errors and convergence rates for the 3rd order scheme employing 3 AMR grids with 3 levels of refinement. They are compared with the fixed grid schemes. The grid is initially refined around the vortex where the gradients are large and the mesh is coarse elsewhere. Notice that the AMR parameters are not optimised for this test, but the refined grid follows the vortex. We observe in table 3 that the optimal rates of convergence are reached also when AMR grid is used. Moreover the AMR errors are about ∼ 1.7 -1.8 times larger than the fixed grid ones for the finest resolution, that is comparing 3.31 × 10 -4 for the AMR grid and 1.91 × 10 -4 for the fixed grid (128 × 128 mesh). Compared to the coarsest fixed grid results (32 × 32), the AMR errors are about 32 times smaller.

Point-wise reconstruction AMR Fixed Fixed grid results AMR grid results Ratio 3rd-order scheme ∼1.75 ∼1.77

(l, L) N max c L 1 err L 2 err L 1 ord L 2 ord L 1 err L 2 err L 1 ord L 2 ord L 1 L 2 (8,
Table 3: L 1 and L 2 errors and convergence rates for the unlimited third-order HO Finite Volume scheme equipped with the current point-wise reconstruction with fixed or AMR grids.

M 1 model: 2D Riemann problem with shocks

The test case considered here is a 2D 4-state Riemann problem for the M 1 model ( 6), which was first proposed in The solution is given by

U 0 (x, y) = (E 0 R (x, y), F 0 R,x (x, y), F 0 R,y (x, y)) =                  (E 0 , F 0 , 0) if x ≤ 0.5 and y ≤ 0.5, (E 0 , 0, F 0 )
if x ≤ 0.5 and y ≥ 0.5, (E 0 , 0, -F 0 ) if x ≥ 0.5 and y ≤ 0.5, (E 0 , -F 0 , 0) if x ≥ 0.5 and y ≥ 0.5, [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] where E 0 = 7.56 × 10 -4 J • m -3 , and

F 0 = (1 -10 -8 )cE 0 W • m -2 .
It is a particularly demanding test case for several reasons. On the one hand, the solution remains close to the limit of admissible states in a large part of the domain throughout the simulation. On the other hand, it presents complex structures, including filamentary shocks which are delicate and difficult to capture with low order schemes or brutal limiting strategy. To achieve a good approximation, a numerical method must be highly accurate and robust. This problem is tested on fixed Cartesian meshes and the results are shown in Figure 6. One represents at t final = 2 × 10 -9 s the values of the anisotropy factor f = F R cE R and the radiative energy E R and the distribution of the polynomial degrees used for the reconstruction. The results obtained using the point-wise reconstruction on meshes with 128 2 and 256 2 cells are compared to those given by the classical scheme using cell-centered polynomial reconstruction on a mesh with 256 2 cells. It can be seen that the schemes capture all the characteristics of the solution with 256 2 cells, which is still relatively coarse. No important difference is visible between the scheme with point-wise reconstruction and the classical one, which confirms the results observed in the previous section. Moreover, only a few cells switch back to the first order in figure 6-bottom. The highest order scheme is therefore used for most cells. This example shows that the point-wise reconstruction provides approximations which are comparable to those given by the classical cell-based FV scheme. HOp scheme can be tested with confidence in an AMR simulation code.

Euler model: cylindrical Sod shock

This problem is used as a sanity test to validate the ability of the scheme to solve a 1D cylindrical problem with separate waves unaligned with the mesh. AMR is activated in all remaining experiments allowing testing the quality of the full simulation code. The initial condition on the domain Ω = [-1, 1] × [-1, 1] is split into the following inside and outside cylindrical states H 0 (r) = (ρ 0 (r), u 0 (r), v 0 (r), p 0 (r)) =

(1, 0, 0, 1) if r ≤ 0.5, (0.125, 0, 0, 0.1) if r > 0.5

C = x 2 + y 2 , γ = 1.4. (39)
These two states generate moving cylindrical waves: a rarefaction traveling inward, and a contact discontinuity plus a shock wave traveling outward. The computation is performed until final time t final = 0.25. For this test we present the density profile in colors and the associated AMR mesh.

Comparison between different levels of refinement. First of all we present the behavior of the current HOp scheme when a (l, L) = (4, 7) and (5, 8) sets of parameters are employed. In figure 7 we present the density (colors) and the mesh at final time. The coarsening acts every other iteration, with parameter C = 8 × 10 -1 . We observe that the AMR mesh is correctly refined on steep gradients and coarsened on regular regions. Both simulations capture the main waves and, as expected the solution is visually sharper for the (5, 8) mesh. As expected, a smaller coarsening parameter C = 8 × 10 -2 generated less coarsening (bottoms panels) showing that the threshold parameters in a AMR context have an important impact although fixing them a priori is not an easy task.

Mesh adaptation in time. We present the case (l, L) = (5, 8) and C = 8 × 10 -1 in figure 8-top. We observe that the mesh is appropriately refined and coarsened, following the main waves in time. The final time result of this simulation is plotted on the middle panels of figure 7. Moreover, on the right panels we present the cell polynomial degrees employed showing that the maximal degree 2 is used on safe zones, while 0 is selected around steep gradients. There are also some cells that switch back to first order in the areas where the solution is flat. In this case however, the quality of the approximation is not negatively impacted. Finally, in figure 9 we present the number of cells as a function of time with comparison to the maximal and minimal allowed numbers related to the (l, L) = (5, 8) parameters. We observe that the AMR strategy gently increases the number of cells, without any spurious effects, to reach about 37% of the maximal value and an average of 26% for the entire simulation. These results on the cylindrical Sod problem show that the scheme manages to adapt the mesh to the underlying solution, also maintaining an almost rotation-invariant mesh. Moreover the mesh is refined everywhere the solution presents some features. The results become more precise as the number of levels of refinement is increased. Notice that the 1D profiles of density confirm these statements which we omit to save space. We have observed that the coarsening parameter plays a major role in obtaining a more or less adapted mesh, i.e a mesh that is appropriately coarsened in space/time according to the underlying solution. In our case these parameters have been chosen only to illustrate the behavior of the AMR numerical scheme. More tailored parameters would certainly produce even better AMR meshes. However, the focus of this work is not to discuss how to achieve such optimal parameters, should they exist. At last we have observed that the distribution of polynomial degrees follows the main waves, and, the number of AMR cells remains at the order of 33% of the maximal allowed number of cells for this specific test case.

Euler model: a 4-shock 2D Riemann problem

This test and the one in section 4.5 are four state 2D Riemann problems for the Euler system of conservation laws (see [START_REF] Schulz-Rinne | Classification of the riemann problem for two-dimensional gas dynamics[END_REF] or [START_REF] Kurganov | Solution of two-dimensional riemann problems for gas dynamics without riemann problem solvers[END_REF] for a numerical study). Both act on the computational domain

Ω = [-1, 1] × [-1, 1]
. This first 2D Riemann problem considers the following initial condition H 0 (x, y) = (ρ 0 (x, y), u 0 (x, y), v 0 (x, y), p 0 (x, y))

=                 
(0.138, 1.206, 1.206, 0.029) if x ≤ 0.5 and y ≤ 0.5, (0.5323, 1.206, 0.0, 0.3) if x ≤ 0.5 and y ≥ 0.5, (0.5323, 0.0, 1.206, 0.3) if x ≥ 0.5 and y ≤ 0.5, (1.5, 0.0, 0.0, 1.5) if x ≥ 0.5 and y ≥ 0.5. (40) The final time is set to t final = 0.3. This initialization corresponds to the 3rd configuration in [START_REF] Schulz-Rinne | Classification of the riemann problem for two-dimensional gas dynamics[END_REF][START_REF] Kurganov | Solution of two-dimensional riemann problems for gas dynamics without riemann problem solvers[END_REF] which develops four main shocks and a central zone of complex interactions. Three AMR configurations are run: (l, L) = (4, 7), [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes[END_REF] and [START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes[END_REF][START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes[END_REF], where (l, L) = (8, 8) corresponds to a fixed mesh made of 256 2 square cells. In figure 10 we present the density variable and the AMR meshes colored by density. The bottom left panel shows the fixed fine mesh result. The middle panels show the AMR results with [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2d unstructured meshes[END_REF], while the top ones present the (4, 7) AMR ones. The AMR simulations have a well-adapted mesh to the underlying solution. As expected one more level of refinement helps to capture more structures in the central part. The bottom-right panel presents the evolution of the number of cells as a function of time. We observe that less than 10% of the maximal number of cells (at level 8) is attained, and, on average only 7% are used. Even with the extra-cost of the AMR procedure the gain is largely in favour of the AMR simulation compared to the fixed fine mesh simulation. Next, in figure 11 we focus on the polynomial degree distribution at the final time for the three simulations. Interestingly most of the cells are updated with the unlimited highest order scheme. The limited ones are aligned with the main shocks, and, appropriately located just ahead and behind. We also observe that some cells are troubled even if they seem located on "flat" regions. But this does not deteriorate the accuracy. Finally, in figure 12, we report the distribution of cells in each level of refinement (left) and the percentage of covered surface (right) for the (l, L) = (4, 7) simulation for which the maximal number of cells is 128 2 = 16384 while the minimal one is 16 2 = 256. We gather some of these percentages in table 4. We observe that the surfaces covered by each level are well balanced at about 25% at the end of the simulation. This solution is self-similar and complex, therefore it is expected that the number of fine cells increases over time.

Euler: 4-contact-wave 2D Riemann problem

This next 2D Riemann problem for Euler's equations corresponds to the 6th configuration in [START_REF] Schulz-Rinne | Classification of the riemann problem for two-dimensional gas dynamics[END_REF][START_REF] Kurganov | Solution of two-dimensional riemann problems for gas dynamics without riemann problem solvers[END_REF], its solution contains 4 planar contact discontinuities which interact in the central region. The initial condition on the computational domain is H 0 (x, y) = (ρ 0 (x, y), u 0 (x, y), v 0 (x, y), p 0 (x, y)) Two different sets of simulations are performed. The first set employs the HLL [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] numerical flux while the second one uses HLLC one [START_REF] Toro | Riemann Sovlers and Numerical Methods for Fluid Dynamics[END_REF]. The later is known to better capture contacts. Indeed, in figure 13 we present the specific internal energy at final time t final = 0.3 for the AMR meshes (l, L) = (4, 7) (left panels) or (l, L) = (5, 8) (right panels). The top/bottom panels present the HLL/HLLC results. For this test case the coarsening does not truly act, hence the results are mainly those obtained by the scheme on a fixed finest mesh. As expected one more level of refinement produces more accurate results. Also the HLLC Riemann solver produces sharper and finer central structures which was expected. Next, in order to understand why the coarsening does not seem to act much, we propose in figure 14 the AMR (5, 8) results with HLLC Riemann solver for several intermediate times. A bird's eye view is plotted, that is a 3D view seen from top with elevation and color representing the value of the variable. Here, the specific internal energy (top panels) and the AMR mesh (bottom panels) are presented. For a different perspective, figure 15 shows the AMR (4, 7) meshes colored by the specific internal energy with HLLC Riemann solver for several intermediate times. There are several reasons why the mesh remains at high levels of refinement over a significant portion of the domain. In addition to the initial discontinuities, the central structure is found to be relatively large and contains high gradients. It is therefore expected that the mesh be refined in the corresponding area. In addition, extra 1D waves of low amplitude can be seen propagating. These waves are artifacts due to the fact that the initial states only approximately lie on the curves of a contact. A good HO scheme must capture and follow those waves which are indeed the correct solution for this problem. Any extra coarsening on these waves would represent a failure of the AMR strategy to maintain the highest accuracy possible on physical waves. A similar phenomenon also occurs in the 4-shock case studied earlier.

=                  (
One could probably obtain a little more coarsening by relaxing the MOOD criterion a little less, however one does not expect the gain to be substantial. This test case is therefore not genuinely favorable to an AMR strategy. It is interesting to look at the impact on the computation time on this unfavorable test case and this is the subject of paragraph 4.7.

M 1 model: 2D Riemann problem with contact discontinuities

The second test-case for the M 1 model is original and called twist in the following. The initial solution is given by

U 0 (x, y) = (E 0 R (x, y), F 0 R,x (x, y), F 0 R,y (x, y)) =                  (E 0 , F 0 , 0) if x ≤ 0.
5 and y ≤ 0.5, (E 0 , 0, -F 0 ) if x ≤ 0.5 and y ≥ 0.5, (E 0 , 0, F 0 ) if x ≥ 0.5 and y ≤ 0.5, (E 0 , -F 0 , 0) if x ≥ 0.5 and y ≥ 0.5, [START_REF] Morrell | A cell by cell anisotropic adaptive mesh ale scheme for the numerical solution of the euler equations[END_REF] where E 0 = 7.56 × 10 -4 SI and F 0 = (1 -10 -8 )cE 0 SI and the final time is set to t final = 1.5 × 10 -9 s. This test contains relatively narrow structures and steep gradients that evolve over time entrained by a rotation like motion. Moreover, some of these structures are associated with contact discontinuities. Likewise for Euler equations, the contact discontinuities of the M 1 model are intermediate waves which are somewhat difficult to capture for a Rusanov numerical flux. We consider three meshes: a fixed Cartesian one with (l, L) = (8, 8) corresponding to 256 2 cells, a (l, L) = (4, 7) AMR mesh and a (l, L) = (5, 8) one. The HO FV scheme with point-wise reconstruction is considered for the AMR meshes whereas the classical HO FV with cell-centered reconstruction is used for the fixed grid. The results can be seen in the figure 16 for the radiative energy variable E R on the left panels and the anisotropy factor f on the right ones. We observe that the HOp scheme with the same nominal mesh accuracy (middle panels) produces comparable results with the fixed mesh and cell-centered reconstruction (bottom panels) of figure 16. As expected, if one less level of refinement is allowed (top panels), then the sharpness of the waves is less pronounced and fewer small scale features are captured. Nonetheless the overall solution is well resolved. Next, in figure 17 we plot the last polynomial degrees used and observe that most of the cells are updated with high order. Only few cells in the plateaus and on the surrounding area of the main waves are still limited. It does not seem that the HOp scheme increases the number of troubled cells compared to the classical HO FV scheme with cell-centered reconstruction. Finally, in figure 18 we present the AMR mesh and the radiative flux vector (F R,x , F R,y ) for the (4, 7) case for three intermediate times. The mesh and the flux are colored with respect to the magnitude of the anisotropy factor. We can clearly see the rotational structure in the center and the steep gradients which are appropriately anticipated and covered by the finest mesh level. The overall AMR procedure based on MOOD criteria is therefore able to maintain the highest accuracy where appropriate. Of course, tailored coarsening parameters would have increased the number of coarsened zones but at the price of lengthy search for these parameters. While feasible, this would not strengthen this proof of concept.

CPU time

In this section we present a rough comparison of CPU times for the Riemann problem 3 for Euler model and the twist problem for M 1 model. Our AMR simulation code has only been designed to be a proof of concept for the HO FV scheme with point-wise reconstruction and the use of MOOD criteria as an AMR refinement strategy. Therefore our sequential AMR code is by no means optimal or even designed towards efficiency, which may be the topic of a future work. Nonetheless we can still compare the relative cost of the AMR simulations against the fixed grid ones, recalling that those 'efficiency' records could drastically be improved. The results are gathered in Table 5, where we express the CPU time ratio with respect to the fixed grid simulation time. The two problems illustrate either a relatively favorable situation or a genuinely unfavorable one for the AMR efficiency. For the Riemann problem we see that for a comparable quality of approximation, see figure 10, the use of AMR allows to speed up the computation by a factor of 2.3. It should be noted that this value is far from being optimal since on the one hand the code can still be improved, and, on the other hand, the refinement and coarsening parameters have been fixed to default values for all test cases. Even in this context using an AMR code significantly reduces the CPU time. For the M 1 Twist test case, see figure 16, for which an excessive refinement is observed, one notices a gain in computation time of about 18%. This gain is obviously less important than for the previous case, but still significant. Even if one expects to save time on the first iterations where only few cells are refined, the interaction with the AMR library usually lowers the expected gain and, consequently, lowers the overall efficiency. These results could be improved by choosing well suited AMR parameters for each test case, in particular the coarsening ones. However we observe that even for the unfavorable test cases and our non optimal code, the AMR strategy improves the overall efficiency. Table 5: CPU time comparisons between AMR and fixed grid simulations.

Conclusion and perspectives

In this paper we have introduced an innovative point-wise polynomial reconstruction procedure to pair up with a highorder Finite Volume (FV) scheme stabilized by an a posteriori Multi-dimensional Optimal Order Detection (MOOD) paradigm within a local cell-based Adaptive Mesh Refinement (AMR) framework, called HOp. This scheme is devoted to solve hyperbolic systems of partial differential equations in multi-dimensions, such as Euler or radiative M 1 systems in this work.

When the mesh is modified during the AMR stage of an HO FV scheme, the cell-centered reconstruction matrices and associated stencils must be rebuilt. This is an expensive cost of the overall method especially when frequent refinement/coarsening stages occur. We have designed a new point-wise polynomial reconstruction which is essentially independent of the surrounding stencil because it is based on point-wise values at fixed points in the vicinity of the current cell. As such it does not require such an expensive stencil and matrix rebuilding. We have analysed the stability of some point-wise based reconstruction schemes by varying the weight distribution. Next we have included this reconstruction within a HO FV scheme limited by a MOOD procedure. Some tests for Euler and M 1 models have shown that this scheme produces equivalent accurate solutions than the classical cell-centered reconstruction based scheme on fixed grids. Additionally, we have revamped the AMR refinement procedure. Classically it requires the user to choose a diagnostic, for instance the norm of the density gradient, and one threshold value. Both the diagnostic and the threshold are not obvious to determine. Worse, they may depend on the test case. We have re-used the MOOD detection criteria to determine which cells must be refined. In short, a detected troubled cell is first refined. Then, only in the second MOOD iteration it is limited. Coarsening is applied at the very end of the time-step on an admissible solution using a classical diagnostic/threshold procedure. The coarsening step does discard some information by averaging, and, as such must be carefully manipulating. In our scheme only genuinely flat regions are allowed for coarsening. On the contrary, refinement must occur to ensure that the high accuracy is maintained wherever the HO FV scheme is acting. This HOp FV scheme has been implemented within the p4est AMR framework [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF] on 2D quad-tree mesh. Several schemes up to order 3 -to keep a consistent expression of the boundary conditions-have been implemented and tested on a series of examples of increasing difficulty. The test cases assess that the AMR mesh appropriately adapts to the underlying solution. Moreover, the solution seems to be captured with high accuracy and no positivity issue as expected. As such it validates the use of the MOOD detection criteria as a refinement criterion. In addition the troubled cells detected and limited by MOOD seem appropriately located. We have also observed that the total number of cells remains on average largely smaller than the one of the finest possible mesh. Although the implementation is far from being optimal, some gains in CPU time are observed. The gain is more or less important for the test cases depending on the surface covered by the fine scale structures, and, the finest cell level. Obviously with a fine tuning of the coarsening parameters we could improve the CPU time by reducing the number of fine cells, but this is not the goal of this article to play with parameters. Here we focus on maintaining the high accuracy of the FV scheme in an AMR framework simplifying the refinement criteria and reducing the cost of reconstruction. Notice that the scheme is also flexible in terms of the system of PDEs and baseline numerical flux (Rusanov, HLL, HLLC. . . ). The results show the applicability of the point-wise reconstruction and the MOOD-based refinement technique.

From this point on, several extensions are envisioned: A 3D version, a parallel and optimized code being the most natural possibilities without any additional limitation. One can also consider other types of models, for example taking into account certain source terms and for which numerical schemes based on approximate Riemann solvers do exist.

  two AMR meshes are plotted in panels (b) and (c). On panel (b) the central cell c is a degenerate heptagon while in (c) it is a degenerate hexagon. The boundary of any cell d Valid AMR mesh, c is a degenerate heptagon. c (c) Invalid AMR mesh. Two levels of refinement for blue/black cells.

Figure 1 :

 1 Figure 1: Notation and some valid/invalid AMR meshes. The invalidity comes from the fact that blue cells can not share an edge with black ones.

  U d c = U d , s.t d ∈ S d c from a stencil S d c . From this, the operator reconstructs a polynomial of degree d which is conservative and corresponds to the least square approximation of the cell mean values U d c . c (a) Set of neighbors S d i (red crosses) for polynomial reconstruction in cell c.

  Quadtree of an AMR mesh where the color corresponds to the level of refinement. c (c) Example neighbor cell stencil (red+blue crosses) based on face-neighbors.

Figure 2 :

 2 Figure 2: (a) Example of the set of neighbors S d i (red crosses) adapted for polynomial reconstruction of degree d = 1 or 2 in central cell c for an AMR mesh -(b) Example of quadtree mesh -(c) Example of the set of neighbors S d i (red and blue crosses) based on face-neighbors. Geometrical bias can be created towards large cells.

Figure 3 :

 3 Figure 3: Sketch of a FV scheme stabilized by a MOOD loop.

1 2 (

 2 d + 1)(d + 2) given point values ( xp , U p , w p ) 1≤p≤P located at position xp close enough to cell c and with weight w p > 0 for 1 ≤ p ≤ P. Let us denote the set of points by P d c = xp , for 1 ≤ p ≤ P . The least-square minimization procedure (22) is replaced by

Figure 4 :

 4 Figure 4: Left: 8 point distribution P d c around cell c usable for 1 ≤ d ≤ 2 on a uniform mesh -Middle: 16 point distribution for 1 ≤ d ≤ 4 on a uniform mesh -Right: the same 16 point distribution on an AMR mesh.

  x1 -I 10 ŷ1 -I 01 x1 ŷ1 -I 11 . . . 1 x2 -I 10 ŷ2 -I 01 x2 ŷ2 -I 11 . . . . . . . . . . . . 1 xP -I 10 ŷP -I 01 xP ŷ2 -I 11 . . .

Figure 5 :

 5 Figure 5: Sketch of the FV scheme stabilized by a MOOD loop within our AMR framework.

Figure 6 :

 6 Figure 6: M 1 model -Filament test case -Simulation at t final = 2 × 10 -9 s -Comparison of the HO FV point-wise scheme results with fixed meshes (l, L) = (7, 7) (left column) and (8, 8) (middle) versus the classical cell-centered based reconstruction FV scheme with ∆x = 1/2 8 (right) -From top to bottom: anisotropy factor, energy and polynomial degrees. The same legend applies to each line.

Figure 7 :

 7 Figure 7: Euler equation -Cylindrical Sod problem -Simulation at t final = 0.25 -Density and mesh -Comparison of the results with levels of refinement (l, L) = (4, 7) (top), (5, 8) with C = 8 × 10 -1 (middle) and (5, 8) with C = 8 × 10 -2 (bottom).

Figure 8 : 9 :

 89 Figure 8: Euler equation -Cylindrical Sod problem -Density and mesh -Comparison of the results for (l, L) = (5, 8) levels of refinement and coarsening parameter 8 × 10 -1 -Top to bottom: different intermediate times starting at t = 0 -Left: density and mesh -Right: cell polynomial degrees.

Figure 10 :

 10 Figure 10: Euler equation -4-shocks Riemann problem -Simulation at t final = 0.3 -Comparison of the results with different levels of refinement (l, L) = (4, 7) (top), (4, 8) (middle) and (8, 8) (bottom) Left column: Density. Right: Mesh colored by density -Bottom-right: evolution of the number of cells for the (4, 7) simulation.

Figure 11 :

 11 Figure 11: Euler equation -4-shocks Riemann problem -Simulation at t final = 0.3 -Polynomial degrees -Comparison of the results with (l, L) = (4, 7) (left), (4, 8) (middle) and (8, 8) (right).

Figure 12 :

 12 Figure 12: Euler equation -4-shocks Riemann problem -Simulation at t final = 0.3 -Left: Number of cells in each level as function of iteration number -Right: Percentage of surface covered by each level of refinement.

Figure 13 :

 13 Figure 13: Euler model -4-contact-waves Riemann problem -Simulation at t final = 0.3 -Specific internal energy -Comparison of the results with AMR (l, L) = (4, 7) (left) or (l, L) = (5, 8) (right) and HLL (top) or HLLC Riemann solver (bottom).

Figure 14 :

 14 Figure 14: Euler model -4-contact-waves Riemann problem -Intermediate times -Specific internal energy (color and elevation) (top panels) and mesh (bottom panels) -AMR (l, L) = (5, 8) mesh and HLLC Riemann solver.

Figure 15 :

 15 Figure 15: Euler model -4-contact-waves Riemann problem -Intermediate times -Specific internal energy and AMR mesh -Results with AMR (l, L) = (4, 7) and HLLC Riemann solver.

Figure 16 :

 16 Figure 16: M 1 model -Twist test case -Simulation at t final = 1.5 × 10 -9 s -Comparison of the results with AMR meshes (l, L) = (4, 7) (top) and (5, 8) (middle) versus the classical cell-centered based reconstruction FV scheme with ∆x = 1/2 8 (bottom) -Left: Energy, right: anisotropy factor.

Figure 17 :

 17 Figure 17: M 1 model -Twist test case -Simulation at t final = 1.5 × 10 -9 s -Comparison of the results with AMR meshes (l, L) = (4, 7) (left) and (5, 8) (middle) versus the classical cell-centered based reconstruction FV scheme with ∆x = 1/2 8 (right) -Polynomial degrees.

Figure 18 :

 18 Figure 18: M 1 model -Twist test case -Simulation at t final = 1.5 × 10 -9 s -Results with AMR mesh (l, L) = (4, 7) at three intermediate times -Top: Radiative flux vectors -Bottom: Meshes. Both are colored with respect to the magnitude of the anisotropy factor.

  L) = (4, 7) (l, L) = (5, 8) (l, L) = (8, 8)

Table 1 :

 1 Maximal CFL number (16) computed by the von Neumann analysis for the upwind HOp FV scheme on 2D avection equation.

			P 1,0						P 1,1						P 2,1
	Points w cn w e w c CFL max		Points w cn w e w c CFL max		Points w cn w e w c CFL max
	12	0	10	1	0.0		12	0	10	1	0.92		16	0.1 0.1	1	0.96
	12	0	0.1	1	0.0		12	0	0.1	1	0.52		16	10	0.1	1	0.96
	16	0.1 0.1	1	0.04		16	0.1 0.1	1	0.68		16	10	5	1	1.00
	16	10	0.1	1	0.20		16	10	0.1	1	0.0		16	1	1	1	0.96
	16	10	5	1	0.20		16	10	5	1	0.96		16	0.5 0.8	1	0.92
						P 2,2						P 3,2
				Points w cn w e w c CFL max		Points w cn w e w c CFL max
				16	10	0.1	1	0.96		16	10	0.1	1	1.04
				16	10	5	1	0.96		16	10	5	1	1.00
				16	1	1	1	0.96		16	1	1	1	1.00
				16	0.5 0.8	1	1.00		16	0.5 0.8	1	1.00

Table 4 :

 4 Euler equation -4-shocks Riemann problem -AMR simulation with (l, L) = (4, 7) levels -Number of cells per level and percentage of surface covered by each level for several iterations. The maximal number of cells is 16384.

	Iteration	= 4	Number of cells per level = 5 = 6 = 7 Total	Percentage of surface covered by level = 4 = 5 = 6 = 7 Total
	1	194	122	248	1024	1588	76%	12%	6%	6%	100%
	300	131	144	585	3356	4216	51%	14%	15%	20%	100%
	700	106	164	748	3984	5002	42%	16%	18%	24%	100%
	1000	84	204	814	4488	5590	33%	20%	20%	27%	100%
	1397	65	229	891	4996	6181	25%	22%	22%	31%	100%

Notice that one could consider listing exhaustively all possible types of neighborhood and their associated matrix for any polynomial degree.Perhaps tractable in 2D, this approach is not reasonable in

3D.

Acknowledgments

This work has been partly produced within the "Laboratoire de Recherche Conventionné" Anabase (Applied numerical analysis to (complex) flow simulation in Bordeaux), anabase.math.u-bordeaux.fr. The authors would like to thank the editors and reviewers of the special issue of AMC.