S Diot

R Loubère

S Clain

The MOOD method in the three-dimensional case: Very-High-Order Finite Volume Method for Hyperbolic Systems

Keywords: Finite Volume, high-order, conservation law, polynomial reconstruction, 3D, unstructured, Euler, MOOD, positivity-preserving

The Multi-dimensional Optimal Order Detection (MOOD) method for two-dimensional geometries has been introduced in "A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)", J. Comput. Phys. 230 (2011), and enhanced in "Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials", Comput. & Fluids 64 (2012). We present in this paper the extension to 3D mixed meshes composed of tetrahedra, hexahedra, pyramids and prisms. In addition, we simplify the u2 detection process previously developed and show on a relevant set of numerical tests for both the convection equation and the Euler system that the optimal high-order of accuracy is reached on smooth solutions while spurious oscillations near singularities are prevented. At last, the intrinsic positivity-preserving property of the MOOD method is confirmed in 3D and we provide simple optimizations to reduce the computational cost such that the MOOD method is very competitive compared to existing high-order Finite Volume methods.

First introduced in [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF], the Multi-dimensional Optimal Order Detection (MOOD) method proposes a new strategy to provide third-order approximations to hyperbolic scalar or vectorial problems for two-dimensional geometry with structured meshes. The author then gave an extension in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] to general unstructured 2D meshes where they achieved a sixth-order convergence in space introducing new detection-limitation procedure. The issue we address in the present paper is to extend the MOOD method to three-dimensional geometries with general polyhedral unstructured meshes for the scalar advection equation and the hydrodynamics Euler system. The method casts in the generic framework of the finite volume method but fundamentally differs from the traditional techniques by the specific detection-limitation procedure implemented by the authors. Indeed, classical high-order polynomial reconstruction-based schemes such as MUSCL [START_REF] Kolgan | Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme[END_REF][START_REF] Kolgan | Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a flow over a body under an angle of attack[END_REF][START_REF] Hubbard | Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids[END_REF][START_REF] Buffard | Monoslope and Multislope MUSCL Methods for unstructured meshes[END_REF][START_REF] Kolgan | Application of the principle of minimizing the derivative to the construction of finite-difference schemes fr computing discontinuous solutions of gas dynamics[END_REF] or ENO/WENO methods [START_REF] Harten | Uniformly highorder accurate nonoscillatory schemes I[END_REF][START_REF] Harten | Uniformly highorder accurate nonoscillatory schemes III[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing scheme[END_REF][START_REF] Zhang | Third-order WENO scheme on three dimensional tetrahedral meshes[END_REF][START_REF] Shu | High-order Weighted Essentially Non-Oscillatory schemes for convection dominated problems[END_REF] rely on an a priori limiting procedure to achieve some stability properties. For instance, in MUSCL-like methods unlimited slopes are reduced through the use of limiters to respect some Discrete Maximum Principle or Total Variation Diminishing properties. In the same way, ENO/WENO-like methods employ an essentially non-oscillatory polynomial which provides an accurate solution while preventing undesirable oscillations from appearing.

We state that such limitation strategies are a priori in the sense that only the data at time t n are used to first perform the limitation procedure and then compute an approximation at time t n+1 . Generally, the "worst case scenario" (speculative approach) has to be considered as plausible and, consequently a "precautionary principle" is applied. It results that most of the time the limitation mechanism unnecessarily operates and may reduce the scheme accuracy due to restrictive assessments. The MOOD principle lies in a different approach since we first compute a candidate solution for time t n+1 and use this a posteriori information to check if the proposed approximation is valid. Roughly speaking, we compute a candidate solution without any limitation using local polynomial reconstructions to provide accurate approximation of the flux (the degree is set to a prescribed maximal value). We then detect if this solution locally fails to fulfill some stability criteria (detection of problematic cells) and further decrement polynomial degree only on problematic regions (limitation step) before recomputing a new candidate solution. An iterative procedure (the MOOD algorithm) is carried out by successively decrementing the degree to provide the optimal local polynomial reconstruction for each cell to satisfy the given stability criteria. At the end of the MOOD algorithm, the candidate solution is eligible and turns out to be the approximation at time t n+1 . The a posteriori strategy brings new benefits. We dramatically reduce the number of polynomial reconstructions regarded to the ENO/WENO method since our technique only requires one polynomial function for each cell. Most of the time, the polynomial with maximal degree is employed since the limitation mechanism is only activated for problematic cells (objective approach). From a physical point of view, the positivity preserving property (for the Euler equations as instance) is simply guaranteed by the a posteriori strategy applying a simple detection procedure which checks the physical admissibility of the solution.

The paper is organized as follows. In section 2, we detail the concept of the MOOD method, while the detection criteria are developed in section 3 both for the advection equation and the hydrodynamics Euler system of equations. Numerical tests are proposed in section 4 to prove the efficiency of the method: we first consider the scalar advection equation and show effective high-order of accuracy for regular solutions with the fourth-and sixth-order schemes considering meshes made of hexahedra and pyramids. We then propose an H-shaped discontinuous profile in rotation to verify the non-oscillatory property of the MOOD method. Finally, numerical simulations are carried out for the Euler system to test the method with a nonlinear vectorial problem. As preliminary experiments, the classical 1D test cases, namely the Sod and Lax shock tubes and the Shu-Osher and Woodward-Collela problems, are run on 3D tetrahedral and pyramidal meshes. Then the numerical order of accuracy is checked on the 2D isentropic vortex in motion extended by invariance for a 3D mesh and a realistic 2D test case (introduced in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]), namely the impact of a shock wave on a cylindrical cavity, is carried out on a mesh made of triangular and quadrangular prisms. At last, we present results for the 3D explosion problem on a pyramidal mesh and the interaction of a shock wave with a quarter of cone on a mesh of 1.1 millions of tetrahedra. We moreover provide computational cost (CPU and memory storage) for the 3D explosion problem for the MOOD method for different polynomial degrees. We conclude with section 5 and delineate some future perspectives.

THE MOOD CONCEPT

We consider the generic hyperbolic equation defined on a domain Ω ⊂ R 3 , t > 0 cast in the conservative form

∂ t U + ∇ • F (U ) = 0, (1a) U (•, 0) = U 0 , (1b)
where U = U (x, t) is the vector of unknown functions depending on x = (x, y, z) ∈ Ω and on the time t. We denote by F the so-called physical flux where we shall consider the autonomous case F = F (U ) (Euler system as instance) and the non-autonomous situation F = F (x, U ) such as ∇ x .F (x, .) = 0 (scalar convection case). Function u 0 stands for the initial condition while the boundary conditions will be prescribed in section devoted to the numerical simulations.

Framework

In order to design the numerical scheme, we introduce the following notation illustrated in Figure 1. The computational domain Ω is assumed to be a polyhedron bounded set of R 3 divided into polyhedral cells K i , i ∈ E el where E el is the cell index set. For each cell K i , we denote by c i the cell centroid, and define the set ν(i) of all the indexes j ∈ E el such that elements K j share a common face f ij with K i and the set ν(i) of all the indexes j ∈ E el such that K i ∩ K j = ∅ (see illustrations in Figure 2). Moreover for each face f ij = K i ∩ K j , n ij stands for the unit normal vector going from K i to K j and we denote by ξ ij,r , q ij,r , r = 1, ..., R ij the quadrature rule for the numerical integration on f ij where ξ ij,r is the weight associated to the r th quadrature point q ij,r with R r=1 ξ ij,r = 1, ∀i ∈ E el and ∀j ∈ ν(i) (see Figure 1).

Figure 1. Notation for a three-dimensional mesh: exploded view of the face f ij between two cells K i and K j . Centroids are respectively denoted by c i and c j . Three quadrature points q ij,r , r = 1, 2, 3, on f ij are drawn for illustration. The unit normal vector pointing from K i to K j is denoted n ij .

S. DIOT

To avoid a specific treatment of the boundary faces we introduce the notion of virtual cell. To this end, assuming that cell K i has a face f ie = K i ∩ ∂Ω on ∂Ω, we introduce the virtual cell K j where j / ∈ E el obtained by symmetrical transformation of the original cell K i which represents the exterior side of Ω. We shall denote by E bd the index set of all virtual cells such that E el = E el ∪ E bd is the index set of all cells (including the virtual ones). The generic first-order Finite Volume scheme associated to equation (1) writes

U n+1 i = U n i -∆t j∈ν(i) |f ij | |K i | F(U n i , U n j , n ij ), (2)
where U n i is an approximation of the mean value of U at time t n on K i , F(U n i , U n j , n ij ) is a numerical flux which satisfies the properties of consistency and monotonicity for the scalar case, ∆t stands for the time step while |f ij | and |K i | are the area of face f ij and the volume of cell K i respectively.

To provide high-order finite volume schemes, we use convex combinations of the initial buildingblock [START_REF] Abgrall | Essentially non-oscillatory Residual Distribution schemes for hyperbolic problems[END_REF] with better approximations at the quadrature points to compute the numerical flux (see [START_REF] Clain | A high-order finite volume method for hyperbolic systems: 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] for instance). The high-order schemes are thus obtained from an original first-order Finite Volume and that is of crucial importance from a computational and implementation point of view (re-use of the original first-order code to achieve high-order approximations). We substitute the first-order approximations of the flux integral by higher-order versions, the scheme then writes

U n+1 i = U n i -∆t j∈ν(i) |f ij | |K i | Rij r=1 ξ ij,r F(U n ij,r , U n ji,r , n ij ), (3)
where U n ij,r , U n ji,r are high-order approximations of U at quadrature points q ij,r on both side of f ij . For meshes constituted of tetrahedral cells, all faces are triangles. Consequently R ij and ξ ij,r are independent of i and j and the previous scheme rewrites as a convex combination of the first-order scheme (2)

U n+1 i = R r=1 ξ r   U n i -∆t j∈ν(i) |f ij | |K i | F(U n ij,r , U n ji,r , n ij )   . ( 4
)
Remark 2.1. When dealing with general polyhedral cells, the quadrature rules may be different from a face to another and such a convex combination is not valid anymore. However, since each polygonal face can be split into triangles, one can recover equation ( 4) by considering each polyhedron as a polyhedron only constituted by triangular faces (and consequently with more faces than the original one).

Let denote by

U n h = i∈E el U n i 1I Ki the constant piecewise representation of approximation (U n i ) i∈E el , we introduce operator H R (U n h ) such that relation (3) rewrites as U n+1 h = U n h + ∆tH R (U n h ). (5)
From the original forward Euler discretization in time [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] we derive a high-order approximation in time using a Runge-Kutta 3 TVD method:

U n+1 h = U n h + 2U (3) h 3 with        U (1) h = U n h +∆t H R (U n h ) U (2) h = U (1) h +∆t H R (U (1) h ) U (3) h = U (2) h +∆t H R ( U (2) h ) (6)
where U (2) h is the convex combination

(3U n h + U (2) h )/4.
The time discretization introduces a 3 rd -order error which makes the whole scheme to be formally 3 rd -order accurate. However setting ∆t = ∆x r/3 where r is the spatial order of accuracy and ∆x is a characteristic length provide same order for spatial and time errors.

Reconstruction

We have formally defined an arbitrary high-order accurate Finite Volume scheme, providing that U ij,r is a high-order accurate point-wise approximation of U (q ij,r ) computed from cell K i . In this subsection, we briefly describe the technique to produce such approximations and we refer to [START_REF] Clain | A high-order finite volume method for hyperbolic systems: 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] and references herein for details. Let us consider a scalar variable u and denote by u i (•, d) a local polynomial approximation of degree d reconstructed on cell K i from the mean values of function u on a set of neighboring cells S d i called stencil. For the sake of conservation, i.e.

1
|Ki| Ki u i (x; d) dx = u i , we assume that the polynomial has the following structure

u i (x; d) = u i + 1≤|α|≤d R α i (x -c) α - 1 |K| K (x -c) α dx , (7)
where the polynomial coefficients R α i are fixed by solving a least-squares problem equivalent to minimizing the functional

E = j∈S d i 1 |K j | Kj u i (x; d) dx -u j 2 .
In practice, the polynomial coefficients are obtained by multiplying the pseudoinverse of the least-square problem matrix (that we store in memory) with the vector of mean values on the stencil, see [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] for details. We moreover recall that for the vectorial case, the reconstructions are performed for all the conservative components independently.

Finally, considering that polynomial reconstructions u i (x; d) are provided for all cells K i , i ∈ E el , we compute the approximation at each quadrature point of each face f ij by u ij,r = u i (q ij,r ; d). The so-called (d-)unlimited scheme (3) is thus defined by employing the reconstructed values in the numerical flux without any restriction (i.e. no limitation). Remark 2.2. We recall that the reconstruction process is very time and memory consuming and would like to emphasize that contrarily to WENO methods we consider only one reconstruction stencil per cell and per degree, so that a lot of computational resources are saved.

The MOOD concept

It is well-known that the first-order scheme (2) is robust but tremendously diffusive, while unlimited schemes of higher-order produce spurious oscillations in the vicinity of steep gradients. Limitation mechanisms have been developed to prevent the oscillations from appearing, such as the slope or flux limitation in MUSCL methods [START_REF] Kolgan | Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme[END_REF][START_REF] Kolgan | Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a flow over a body under an angle of attack[END_REF][START_REF] Hubbard | Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids[END_REF][START_REF] Buffard | Monoslope and Multislope MUSCL Methods for unstructured meshes[END_REF][START_REF] Kolgan | Application of the principle of minimizing the derivative to the construction of finite-difference schemes fr computing discontinuous solutions of gas dynamics[END_REF] or the computation of an Essentially Non-Oscillatory polynomial reconstruction in WENO methods [START_REF] Harten | Uniformly highorder accurate nonoscillatory schemes I[END_REF][START_REF] Harten | Uniformly highorder accurate nonoscillatory schemes III[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing scheme[END_REF][START_REF] Zhang | Third-order WENO scheme on three dimensional tetrahedral meshes[END_REF][START_REF] Shu | High-order Weighted Essentially Non-Oscillatory schemes for convection dominated problems[END_REF]. As we mention in the introduction, all the classical techniques act a priori in the sense that we determine the limitation process in function of the current data (i.e. the solution at time t n ). As a consequence, the a priori strategy imposes very drastic accuracy reduction due to strong and unnecessary limitations (the worst case scenario has to be considered). Moreover, computational resources are allocated to perform the limitation process where most of the time it is useless. In the ENO/WENO case for instance, several polynomial reconstructions are required even if the solution is locally regular and can be approximated with only one polynomial function.

In two recent papers [START_REF] Clain | A high-order finite volume method for hyperbolic systems: 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], we have introduced a new approach based on an a posteriori evaluation of the solution to determine if the limitation procedure has to be applied and where. The technique is a posteriori in the sense that we compute a candidate solution (a potential approximation at time t n+1 ) and we use the data of the candidate solution to determine if the solution is valid. More precisely, the detection-limitation mechanism operates in several steps. A candidate solution is first computed with the highest-order unlimited scheme (the polynomials with maximal degree). Then a detection procedure is performed to determine the problematic cells, i.e. all cells where the approximation does not respect some given criteria (see next section). For problematic cells, the solution is recomputed with a lower-order unlimited scheme (using polynomials with lower degree) and we repeat the procedure detection-degree decrementing (the MOOD algorithm) till the cell satisfy the detection criteria or the polynomial degree is zero. In the last case, a robust first-order scheme (2) is triggered and a meaningful solution is thus provided. Note that we need to guarantee that the MOOD algorithm stops after a finite number of iterations.

We now set some fundamental notions to define the MOOD method. We name Cell Polynomial Degree, shortened as CellPD and denoted by d i , the degree of the polynomial reconstruction on cell K i . We name Face Polynomial Degrees, shortened as FacePD and denoted by d ij and d ji , the degrees of the polynomial reconstructions actually used to compute approximations, U ij,r and U ji,r , of the solution on face f ij at quadrature points q ij,r respectively from K i and K j . The computation of d ij and d ji , named FacePD strategy, consists in evaluating the FacePD d ij , d ji that we employ on both sides of the interface f ij with respect to the CellPD of the neighboring cells. In previous studies (see [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] for details), we have proposed and experimented several strategies and introduced the upper-limiting property for a FacePD strategy which states that for any degree d, the following property holds

d i = d =⇒ d ij ≤ d and d ji ≤ d, ∀j ∈ ν(i).
This guarantees (see [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF]) that the MOOD algorithm stops after a finite number of iterations. In practice, we use the simple rule

d ij = d ji = min(d i , d j ).
As mentioned above, the detection mechanism is performed on the candidate solution U ⋆ h and criteria have to be set to specify what is a good solution. To this end, we denote by A the set of detection criteria (e.g. positivity of a variable or a maximum principle) that the numerical approximation has to respect on each cell and we say that a candidate solution is A-eligible if it fulfills all the criteria of A. If the candidate solution is not A-eligible on cell K i , then we decrement the polynomial degree. However the solution may not be A-eligible regardless of the set A even if the polynomial degree is zero for the cell. Consequently, we shall consider the solution acceptable on the cell if either it is A-eligible, or is a first-order solution (i.e. CellPD has been decremented to zero). To sum up, the MOOD algorithm for the explicit discretization in time consists in the following stages depicted in Figure 3: Following [START_REF] Clain | A high-order finite volume method for hyperbolic systems: 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] we extend the MOOD algorithm initially designed for a one-time step scheme to the RK3-TVD scheme by applying it to each sub-step of the RK3-TVD ( 6) procedure. The MOOD method is now completely defined except from the detection criteria that have to be suited to the problem we intend to solve. Such a difficult task requires the complete next section.

0. Initialize d i = d max , ∀i ∈ E el . 1.

DETECTION CRITERIA

The crucial point of the MOOD method is the elaboration of the detection criteria set A which characterizes the properties we want the numerical solution to fulfill. A fundamental purpose of the detection criteria is to obtain higher-order of accuracy for regular solutions while preventing numerical oscillations in the vicinity of discontinuous profiles. This would consequently provide an efficient and robust method. We face several difficulties to design such a set since accuracy and robustness are antagonist objectives. Moreover, in the Euler problem, a physically admissible solution is mandatory since the positivity of the density and the pressure is required to compute the numerical flux. It results that the detection criteria would cover a wide spectrum of properties and restrictions. A key point we shall detail in the following is the notion of "numerical regularity" in the sense that we have to determine if, for a local stencil and a set of data (for instance the mean values), we can associate a regular or a irregular function. This point is really important since the choice of the reconstruction (namely the polynomial degree) depends on it.

The present section intends to extend and improve detection criteria initially introduced in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] to evaluate the local "numerical regularity" of the approximation. We first begin the study for the advection equation in section 3.1 and address the hydrodynamics Euler system in section 3.2.

Advection equation

The scalar advection problem is characterized by the physical flux F (U ) = V U where V ∈ R 3 stands for the velocity that we assume to be a regular function on Ω and satisfies ∇ x V (x) = 0 while U = U (t, x) ∈ R is the passive scalar quantity transported by the fluid.

When dealing with a constant velocity, the exact solution is simply given by U (x, t) = U 0 (x -V t) and clearly fulfills a maximum principle, e.g. the minimum of the solution can not be lower than the initial condition minimum (and a similar property for the maximum). Consequently, it seems natural to impose such a condition at the numerical level and, as proposed in [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF], we integrate in the set A the Discrete Maximum Principle (DMP) on mean values for the candidate solution U ⋆ h formulated like this: min

j∈ν(i) (U n i , U n j ) ≤ U ⋆ i ≤ max j∈ν(i) (U n i , U n j ). (8)
A solution is A-eligible if condition ( 8) is satisfied for all the cells and we have proved in [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] that the scheme equipped with such set A provides a numerical solution which, under first-order scheme CFL condition, satisfies the DMP. However a strict application of relation ( 8) at smooth extrema unavoidably reduces the scheme accuracy to two. It suggests that relation ( 8) is too restrictive and should be relaxed.

In [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] we have relaxed the condition on cells which violate the DMP. More specifically, the relation ( 8) has been supplemented with a new criteria, the so-called u2 detection criteria which provides an effective arbitrary high-order of accuracy. As mention in the beginning of the section, the key point is to determine if the numerical solution is regular enough to be approximated by a highorder polynomial reconstruction and avoid the Gibbs phenomena. To this end, let assume that the candidate solution does not satisfy the DMP criteria on cell K i . A first step consists in reconstructing quadratic polynomials on K i denoted by U i and on its neighbors K j for j ∈ ν(i) denoted by U j . In a second step, we define approximations to the local minimal and maximal curvatures, namely

X min i = min j∈ν(i) ∂ xx U i , ∂ xx U j and X max i = max j∈ν(i) ∂ xx U i , ∂ xx U j , (9)
Y min i = min j∈ν(i) ∂ yy U i , ∂ yy U j and Y max i = max j∈ν(i) ∂ yy U i , ∂ yy U j , (10)
Z min i = min j∈ν(i) ∂ zz U j , ∂ zz U j and Z max i = max j∈ν(i) ∂ zz U i , ∂ zz U j , (11)
where we emphasize that the second derivatives are constant and naturally referred to as curvatures.

The u2 detection criterion holds in the following definition.

Definition 3.1 (u2 detection criterion). A candidate solution U ⋆ i in cell K i which violates the DMP is nonetheless eligible if the following holds

X max i X min i > 0 and X min i X max i ≥ 1 -ε,
and

Y max i Y min i > 0 and Y min i Y max i ≥ 1 -ε,
and

Z max i Z min i > 0 and Z min i Z max i ≥ 1 -ε,
where ε is a smoothness parameter.

The definition derives from the idea that the comparison of local second derivatives of the quadratic reconstructions on a neighborhood provides a relevant information on the numerical smoothness of the underlying solution. More precisely, we consider that the underlying solution (characterized by the piecewise constant mean value) is (ε-)smooth if for each direction the curvatures have the same sign (no oscillation or inflection point) and are (ε-)close enough to eachother. Such a definition lies in a fitting of the parameter ε the value of which defines the threshold between what is considered as smooth extrema or as discontinuity. Therefore the determination of ε is of crucial importance since it rules the decrementing process activation.

From a practical point of view, the u2 detection criteria operates in two stages. First the test on the sign of curvatures (left inequalities) is performed. If oscillations are detected, i.e. the product is negative, the cell is considered as problematic and the decrementing procedure must be applied. The second stage is performed only if the product is positive. It consists in computing the ratios between minimal and maximal curvatures and comparing it to 1ε (right inequalities). If the curvatures ratio does not respect the inequality, the cell is considered as problematic and the decrementing process must be applied.

In [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF], we propose a parameter ε depending on a local characteristic length and on the spatial dimension of the domain. This was a first attempt to the determination of ε and deeper investigations have shown that a simpler definition provides same quality results. To set the ε value, we extend the parameter as a new function ε

x = ε x X min i X max i
(for the x-direction) with respect to the curvatures which have to satisfy the restriction

X min i X max i ≥ 1 -ε x . ( 12
)
The goal is to determine a relevant function ε x which enables high-order approximation and robustness. We first note that the curvatures ratio ranges between zero and one so that ε x must range in [0, 1] to make sense. Moreover, the ratio is expected to be close to zero on discontinuities and close to one on smooth functions which are the two extreme cases. For a non-smooth function, we expect that the limiting procedure operates and since the closer to zero ε is, the less smooth the underlying function is considered, ε x is expected to be close to zero on discontinuities, i.e. lim r→0 + ε x (r) = 0 + . On the other hand, a ratio close to one indicates smooth functions, so we expect lim

r→1 - ε x (r) = 1 -
to relax the restriction. We thus propose to define ε x as a continuous increasing function of the curvatures ratio such that ε x (0) = 0 and ε x (1) = 1. After several attempts, it appears that the simple function ε x (r) = r is an excellent choice. When substituting expression of ε x = X min i /X max i in relation [START_REF] Dumbser | ADER schemes on unstructured meshes for nonconservative hyperbolic systems: Applications to geophysical flows[END_REF], the x-direction curvatures criterion becomes

X min i X max i ≥ 1 - X min i X max i ,
and yields

X min i X max i ≥ 1/2.
Finally we apply the same reasoning for yand z-directions and obtain

Y min i Y max i ≥ 1/2 and Z min i Z max i ≥ 1/2.
The linearity of function ε x simplifies the final inequalities and leads to the constant value ε = 1/2 in definition 3.1.

Remark 3.2. The definition of ε is really simpler than the one proposed in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]. However numerous numerical test cases have been carried out and no change in the quality of results have been reported.

Remark 3.3. Numerical experiments show that the choice of the neighborhood where the curvatures are computed should define a convex hull which contains the reference cell K i . To constitute such a stencil, we used the index set of cells ν(i) in 2D (see [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]) but this choice is not relevant for threedimensional meshes and we use the index set ν(i) in equations ( 9)- [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] to provide the expected results even for large form factor meshes. To conclude the section we propose in Figure 4 an algorithmic view of the complete detection process [DMP → u2] for the advection equation constituted of the DMP of equation ( 8) relaxed by the u2 detection criteria of definition 3.1. We emphasize that the algorithm is given in the case of a cell K i with U ⋆ i its associated candidate solution mean value.

Hydrodynamics Euler system

The Euler system for three-dimensional geometries writes

∂ t      ρ ρu ρv ρw E      + ∂ x      ρu ρu 2 + p ρuv ρuw u(E + p)      + ∂ y      ρv ρuv ρv 2 + p ρvw v(E + p)      + ∂ z      ρv ρuw ρvw ρw 2 + p w(E + p)      = 0, (13)
where ρ stands for the density, u, v and w for the velocity components in the x, y and z directions respectively, p for the pressure and E for the total energy. This system is closed by the Equation Of State (EOS) of a perfect gas p = (γ -1)ρǫ, where ǫ is the specific internal energy, γ the ratio of specific heats and the total energy is constituted of the internal and kinetic energy

E = ρ (u 2 + v 2 + w 2 )/2 + ǫ .
At last, vector U = (ρ, ρu, ρv, ρw, E) represents the conservative variables of the system while W = (ρ, u, v, w, p) are the primitive ones. Note that contrarily to WENO methods we do not use the characteristic variables.

To provide accurate and oscillation-free solutions we use on the one hand polynomial reconstruction and apply, on the other hand the MOOD algorithm. We mention that all the polynomial reconstructions are performed on conservative variables and only one CellPD is used for all variables (see [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] for the motivations and justifications). We now turn to the detection-limitation procedure and we have to design specific detection criteria for the Euler problem.

Following [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF], a first and also mandatory detection criteria corresponds to ensuring the physical meaningfulness of the primitive variables. We then introduce the Physical Admissibility Detection (PAD in short) which considers that the candidate solution on a cell K i is not valid if we have ρ ⋆ i or p ⋆ i are negative (after having computed pressure p ⋆ i ). We underline the important property that a high-order scheme (whichever the degree of the polynomial reconstruction) equipped with the PAD and a first-order scheme which preserves the positivity (of density and pressure) under a CFL condition is automatically positivity preserving. This property straightforwardly derives from the a posteriori nature of the MOOD method and has been proved in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF].

However the PAD detection process does not prevent spurious oscillations from appearing and we turn to the adaptation of the [DMP → u2] detection process proposed in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]. Initially defined for scalar quantity, we apply the [DMP → u2] on the density ρ (detection) and recall that the decrementing is performed for all variables (limitation). Note that the smoothness parameter ε is still set as 1/2 in the u2 definition as in previous section.

The set of constraints A for Euler system is thus constituted by the PAD followed by the [DMP → u2] detection process applied to the density variable since we first check the PAD and if the cell is valid we continue with the [DMP → u2]detection. In Figure 5 we give an algorithmic view of the complete detection process [PAD → DMP → u2] for the hydrodynamics Euler system constituted of the PAD detection criteria, the DMP of equation ( 8) on the density relaxed by the u2 detection criteria of definition 3.1. We emphasize that the algorithm is given in the case of a cell K i with U ⋆ i = ρ ⋆ i , (ρu) ⋆ i , (ρv) ⋆ i , (ρw) ⋆ i , E ⋆ i its associated candidate solution mean value and that the candidate pressure p ⋆ i has to be computed.

We now highlight some implementation aspects about the detection process which enable to improve the solution accuracy. Actually in the above algorithm, the [PAD → DMP → u2] performs well but does not, in some cases, fully reach the optimal order of accuracy for smooth solutions. Deeper investigations on the isentropic vortex in motion problem have shown that the detection process inappropriately decrements some cells of the flat region while it operates well in the area where curvatures are not negligible. The undesirable limitation derives from the extra-small curvatures treatment by the u2 detection where some spurious micro-oscillations take place on the flat area and wrongly activate the curvature sign detection. It results that the sign criterion is not relevant when all the curvatures sizes are too small with respect to a mesh parameter δ. To overcome the over-detection phenomena, we introduce a relaxation parameter in the u2 criterion to fix the problem.

Definition 3.4 (u2 detection criterion). A candidate solution U ⋆

i in cell K i for which the density ρ ⋆ i violates the DMP is nonetheless eligible if

X max i X min i > -δ and max |X max i |, |X min i | < δ or X min i X max i ≥ 1/2 ,
and

Y max i Y min i > -δ and max |Y max i |, |Y min i | < δ or Y min i Y max i ≥ 1/2 , S. DIOT and Z max i Z min i > -δ and max |Z max i |, |Z min i | < δ or Z min i Z max i ≥ 1/2 ,
where δ is the greatest length of geometrical entity of dimension one defined by the length of the cells in R, the maximal length of the cell interfaces in R 2 and the maximal length of edges of the cell interface for three-dimensional meshes.

The correction only damps extra-small oscillations such that minimal and maximal curvatures product satisfies the left condition. When maximal curvatures are larger than δ, the condition on the ratios of curvatures implies that the underlying function will be considered as non-smooth.

Remark 3.5. The value of δ has been determined after numerous simulation experiments. It enables to fully reach the optimal order for the Euler system but does not affect the method in wisely capturing discontinuous profiles. The correction has even been tested for the convection equation and accuracy losses have not been reported.

In the same way, we slightly relax the DMP criteria to reduce the computational effort to avoid the waste of resources when performing the u2 detection criterion on plateaus. We consider that a DMP violation is not relevant if

max j∈ν(i) (ρ RK i , U RK j ) -min j∈ν(i) (ρ RK i , U RK j ) < δ 3
where index RK corresponds to one of the Runge-Kutta sub-steps.

The MOOD method for the Euler hydrodynamics system is now completely defined and numerical simulations are carried out for three-dimensional geometries presented in section 4.

Implementation and optimizations

To conclude the section we detail two important and simple optimizations that we apply to drastically improve the efficiency of the MOOD method.

Local re-updating. The MOOD method may seem computationally expensive since the MOOD algorithm we run for each time step, recompute the candidate solution several times whereas polynomial degrees have been only modified for a small number cells. At the first stage, an initial candidate solution is computed on all cells. Then the MOOD algorithm successively detects and limits the problematic cells. The evaluation of a new candidate solution during the MOOD algorithm by means of scheme (4) only involves the fluxes at the interfaces of corrected cells. Consequently only problematic cells and their neighbors by face must be recomputed. It drastically reduces the computational effort since in most cases the solution is acceptable on more than 80-90% of cells, even when shocks are present.

Reduced polynomial degree decrementing. The original decrementing procedure consists in dropping one-by-one polynomial degrees until zero is reached. Such an approach may be both costly in CPU and memory resources since reconstruction matrices must be stored for all degrees.

It would nonetheless still be less memory consuming than for the WENO method due to the large number of polynomial functions involved in the WENO technique. Moreover numerical experiments suggest the following alternative: whether the solution is very smooth, whether the solution presents some discontinuities. To take advantage of it, we change the decrementing strategy by starting from the highest degree, reducing to degree 2 if any and setting degree equal to 0 if the candidate solution is still not A-eligible. We then manage to reduce the number of decrementing stages and save computational resources. We point out that the size of the reconstruction stencil is also an important parameter since a large stencil (required for the maximal degree) will be influenced by a discontinuity located in the second or third layer of cells around the reference one while a more compact one (for a P 2 reconstruction) still preserves the local regularity of the underlying function. Another reason to use the P 2 reconstruction is that it is also used for the u2 detection process and always has to be stored. Therefore in practice, we only store two reconstruction matrices per cell, one for the maximal degree and one for the degree two. It is thus important to remark that the storage cost of the matrix for degree two is always much lower than the one for the maximal degree. Indeed for two-dimensional situations, the memory cost of the pseudoinverse matrix associated to polynomial of degree 2 represents about 10 times 5 elements, about 16 times 9 for P 3 reconstruction and around 28 times 20 for P 5 . Analogically for three-dimensional situations, the P 2 reconstruction matrix represents about 16 times 9 elements while it is about 38 times 19 for P 3 and 110 times 55 for P 5 .

To conclude this section, we would like to draw some remarks about the potentiality of the MOOD method to be parallelized. Within the MOOD algorithm, only classical unlimited schemes are used without modification so that the parallelization of this part of the method can be done as efficiently as the state-of-the-art methods (WENO method for instance). The only novelty brought by the MOOD method is the iterative process constituting the MOOD algorithm. A potential difficulty comes from the fact that the number of cells on which the numerical scheme acts changes from an iteration of the MOOD algorithm to another, since the procedure is only applied to problematic cells. However it may not dramatically affect the parallelization efficiency: firstly, because an efficient treatment of the list of problematic cells can be achieved and secondly, because the time spent to recompute new candidate solutions is negligible compared to the time to compute the initial one since the number of problematic cells is (in general) very low compared to the total number of cells. The parallelization capacity of the MOOD method is thus as good as the state-of-the-art higher-order finite volume methods.

S. DIOT

NUMERICAL RESULTS

The MOOD method has been implemented into a 3D unstructured code dealing with polyhedra having coplanar faces: tetrahedron, hexahedron, pyramid and prism. The polynomial reconstruction procedure is implemented independently of the degree d max and we provide in the present paper numerical results up to d max = 5. Following [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] and remarks in section 3.3, we use the decrementing sequence P dmax -P 2 -P 0 . The reconstruction matrices are computed and stored in a preprocessing step since they only depend on geometry. Moreover fluxes across faces are approximated by the mean of Gaussian quadrature formulae on a triangular decomposition of the faces (see Figure 1). At last concerning the time discretization, the first-order time step ∆t is controlled by a CFL coefficient equal to 0.5. For the convergence studies on smooth solutions we use the time step ∆t = ∆x r/3 to achieve a global r th -order of accuracy and compute the relative L 1 and L ∞ errors for a a bounded, L 1 function ϕ by

L 1 error: i∈E el |ϕ N i -ϕ ex i ||K i | i∈E el |ϕ ex i ||K i |
and L ∞ error:

max i∈E el |ϕ N i -ϕ ex i | max i∈E el |ϕ ex i | ,
where (ϕ ex i ) i∈E el and (ϕ N i ) i∈E el are respectively the exact and the approximated cell mean values at final time t = t final .

Advection equation

For the scalar advection equation, the MOOD method is employed with the [DMP → u2] detection process and two test cases are carried out: the Triple Sine Translation (TST) to assess the effective very high-order of accuracy and the rotation of a discontinuous H-like shape to test its ability to damp the spurious oscillations.

Triple Sine Translation

Let Ω be the unit cube. We consider a constant translation velocity V = (1, 1, 1) and the C ∞ initial condition U 0 (x, y, z) = sin(2πx) sin(2πy) sin(2πz).

The final time is t final = 2.0 and periodic boundary conditions imply that the exact final solution coincides with the initial one. The computations are first carried out on a series of successively refined regular hexahedral meshes from 8 3 to 64 3 cells. To underline the capability of the MOOD method to handle mixed element meshes, we also consider a series of meshes built from a series of regular hexahedral meshes from 4 3 to 48 3 cells into which we regularly split half of cells into 6 pyramids (see top line of Figure 6).

In Figure 6, we display the convergence curves for the L 1 and L ∞ errors of the MOOD-P 2 , MOOD-P 3 and MOOD-P 5 methods and give in Table I, the corresponding errors and rates of convergence. As expected, the optimal rate of convergence is achieved. Notice that on the coarsest meshes the initial mean values are not representative of the underlying smooth function and are coherently handled by the method as discontinuous profiles. As such the sine function is underresolved; for instance in 1D, averaging the function sin(2πx) or an Heaviside-like function on [0; 1] using four cells provides to same mean values. Table I. L 1 , L 2 and L ∞ errors and convergence rates for the TST problem with the MOOD-P 2 , MOOD-P 3 and MOOD-P 5 methods. Top lines: hexahedral meshes. Bottom lines: mixed hexahedral/pyramidal meshes.

H-like shape rotation

We now turn to the rotation of an H-like shape in the unit cube Ω. The initial shape is given by

U 0 (x, y, z) = 1 if (|x -0.5| > 0.1)
or (|y -0.5| < 0.1), 0 elsewhere, in the cube [0.2; 0.8] 3 and 0 elsewhere. The rotation axis is the diagonal line joining the origin (0, 0, 0) and the point (1, 1, 1). We stop the simulation after one full rotation when the shape is back to its original position. Note that the velocity depends on the spatial position but is divergence-free so that the maximum principle also applies in that case. Numerical simulations are carried out on a 86215 tetrahedra mesh generated by the free mesher Gmsh. Results are displayed with an extruded view on the cut plane z = 1/2. Initialization details are illustrated in Figure 7. We plot in Figure 8 the solution on the cut plane z = 1/2 for the unlimited P 3 and P 5 schemes and the MOOD-P 3 , MOOD-P 5 methods. We notice that the unlimited schemes produce oscillations, depicted in green in the figure, whereas the MOOD method provides an oscillation-free solution even for polynomials of degree 5. It highlights the capacity of the [DMP → u2] detection process to correctly treat discontinuous shapes on genuinely unstructured 3D meshes.

Euler system

We now consider the three-dimensional hydrodynamics Euler system on unstructured meshes. The first test cases proposed in section 4.2.1 deal with the Sod and Lax shock tubes following the Ox axis (invariant with respect to the other directions). The simulations are carried out on a tetrahedral mesh to study the MOOD method capacity to handle simple waves. Section 4.2.2 is dedicated to the Shu-Osher and Blastwave problems approximated on pyramidal cells which respectively involve a complex oscillatory solution and strong interactions between simple waves along the Ox direction. We address in section 4.2.3 the effective numerical accuracy of the method with the isentropic vortex problem for which an exact smooth solution exists. Then in section 4.2.4, we assess the ability of the MOOD method to simulate complex realistic physics on a mesh of triangular and quadrangular prisms by carrying out the impact of a shock wave on a cylindrical cavity proposed in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]. At last, we provide the MOOD method results for two genuinely threedimensional test cases. First we compare the behavior and computational cost (CPU and memory storage) of the MOOD method with different degrees and detection processes by simulating the so-called explosion problem [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, 3 rd revision[END_REF] using unstructured pyramidal meshes in section 4.2.5; Then in section 4.2.6, we consider the interaction of a shock wave with a quarter of cone on a mesh of 1.1 millions of tetrahedra with the 4 th -order MOOD method.

Sod and Lax shock tubes

The original Sod [START_REF] Sod | A survey of several finite difference methods for systems of non-linear hyperbolic conservation laws[END_REF] and Lax [START_REF] Lax | Weak solutions of nonlinear hyperbolic equations and their numerical computation[END_REF] problems concern onedimensional Riemann shock tubes whose solutions consist of a left-moving rarefaction fan, a right-moving contact discontinuity and a right-moving shock wave. In the three-dimensional context, we reproduce the expansion following the Ox axis setting initial condition invariant in y, z and we prescribe reflecting boundary conditions on the cylinder sides. The domain is filled with an ideal gas with γ = 1.4 and the discontinuity is located in x = 0.5 at t = 0. The initial density/velocity/pressure values and final time t f inal are given by The computational domain we consider is a cylinder of unit length and radius R = 0.025 with 0x line as symmetry axis which is paved with 7517 unstructured tetrahedra as shown in figure 9.

Figure 9. Mesh constituted of 7517 tetrahedra used for the Sod and Lax problems. Some cells are drawn non-opaque to see some interior tetrahedra.

We display in Figure 10 the numerical approximations of the density computed with the MOOD-P 3 method using the [PAD → DMP → u2] detection process and the exact solution (red line). In order to provide a clear and relevant representation of the solution along the Ox axis, we slice the whole cylinder in 100 uniform cylinders (since the average characteristic length is 10 -2 ) and plot the average of the solution on each of them. As expected the MOOD-P 3 method provides a very good approximation of the solution and maintains sharp discontinuities. In particular, we underline the very few numbers of points in the contact discontinuity.

Shu-Osher and Blastwave problems

The Shu-Osher problem has been introduced in [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing scheme[END_REF] to test the ability of a scheme to capture both small-scale smooth flow along with shock wave. The one-dimensional computational domain is Ω = [-5; 5] and the final time is t final = 1.8. An initial x-directional shock wave located at x = -0.4 separates the domain into a left post-shock state We consider a perfect gas with γ = 1.4. Reflecting boundary conditions are used to preserve the invariance following axis Oy, Oz except from the left boundary condition which is an inflow one.

The Blastwave problem has been introduced by Collela and Woodward in [START_REF] Colella | The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations[END_REF] to test the performance of numerical schemes on problems involving strong and thin shock structures. The initial conditions consist of two parallel planar flow discontinuities on domain Ω = [0, 1] separated by the planes x 1 = 0.1 and x 2 = 0.9. The density is unity on the whole domain and the gas is assumed initially at rest. The pressure is given by p L = 1000 on the left, p C = 0.01 in the center and p R = 100 on the right. Reflecting boundary conditions are prescribed and the final time is t final = 0.038.

We consider a 21600 regular pyramids mesh (see Figure 12-top right for a pattern example) obtained from a 400 × 3 × 3 regular hexahedral mesh for which each cell is split into six pyramids. The original hexahedral mesh is built by setting ∆x = ∆y = ∆z with ∆x = 0.075 for Shu-Osher problem and ∆x = 0.0075 for the Blastwave problem. Since there is no exact solution for both tests we have computed reference solutions using a first-order finite volume scheme with very fine meshes. As in the previous simulations, the solutions are plotted following the Ox direction considering an underlying 400 points uniform one-dimensional mesh and circles in Figure 11 represent the mean density on three-dimensional slices of thickness ∆x.

Density approximations obtained with the MOOD-P 3 are presented in Figure 11 and compared to the reference solution (red line). For the Shu-Osher problem (left) we report that the [PAD → DMP → u2] detection criteria does not over-smooth the oscillatory solution and accurately capture the high-frequencies waves. On the other hand, for the Blastwave problem (right) we observe sharp contact discontinuities and shock waves are well-preserved. No spurious oscillations are generated and the central structure of the solution is very well approximated.

Isentropic vortex

The isentropic vortex problem was initially introduced for the twodimensional space [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperin Advanced Numerical Approximation of Nonlinear Hyperbolic Equations[END_REF][START_REF] Yee | Entropy Splitting and Numerical Dissipation[END_REF] to test the accuracy of numerical methods since the exact solution is smooth and has an analytical expression. We simply extend the original problem for the threedimensional situation taking the two-dimensional solution invariant following Oz. Let us consider the computational domain Ω = [-5, 5] × [-5, 5] × [0, z max ] and an ambient flow characterized with ρ ∞ = 1.0, u ∞ = 1.0, v ∞ = 1.0, w ∞ = 1.0, p ∞ = 1.0, with a normalized ambient temperature A z-invariant vortex is centered on the axis line x vortex = (x vortex , y vortex , z) = (0, 0, z) z ∈ R and supplemented to the ambient gas at the initial time t = 0 with the following conditions

u = u ∞ + δu, v = v ∞ + δv, T * = T * ∞ + δT * where δu = -y ′ β 2π exp 1 -r 2 2 , δv = x ′ β 2π exp 1 -r 2 2 , δT * = - (γ -1)β 8γπ 2 exp 1 -r 2 ,
with r = x ′ 2 + y ′ 2 and x ′ = xx vortex , y ′ = yy vortex . The vortex strength is given by β = 5.0 and the initial density follows relation

ρ = ρ ∞ T * T * ∞ 1 γ-1 = 1 - (γ -1)β 8γπ 2 exp 1 -r 2 1 γ-1 . ( 14
)
The domain is paved either with N × N × 4 hexahedra , N = 20, 40, 60, 80, 120 or with N × N × 24 pyramids (each hexahedron from the previous mesh is split into 6 pyramids, see Figure 12). To reduce the computational effort, only four cells are considered in the zdirection and z max is taken such that ∆x = ∆y = ∆z, that is to say z max = 4∆x = 40/N . The minimal/maximal number of cells is 1600/57600 hexahedra and 9600/153600 pyramids. We prescribe periodic boundary conditions everywhere.

In Figure 12 we display the convergence curves for the L 1 and L ∞ errors on the density approximations for MOOD-P 2 , MOOD-P 3 , MOOD-P 5 methods, while we provide in Table II the corresponding errors and convergence rates. We report effective orders corresponding to the expected optimal rates of convergence for both types of meshes and underline the MOOD method capacity to provide effective high-order of accuracy on a smooth but non-trivial solution for the three-dimensional Euler system.

Impact of a shock wave on a cylindrical cavity

Based on the experiment proposed in [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF], we have introduced this test case in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] for the two-dimensional case. We here extend it to 3D by invariance along the z-direction. It consists in a planar shock wave which impacts a cylindrical cavity creating complex structures and instabilities. The original purpose of this stringent numerical test is to prove the ability of the MOOD method to capture physics in realistic conditions. In this paper it moreover assesses the capacity of the MOOD method to deal with mixed triangular and quadrangular prisms since the mesh is obtained by extrusion (only two layers) along the Oz axis of a 2D mesh containing 101127 cells (triangles and quadrangles). We moreover point out that important differences between cell sizes are present in the domain, since the largest characteristic length is 0.008 and the smallest one is 0.00015. At last, we run the simulation on the lower half part of the domain but plot a full domain using a symmetry argument. Details of the mesh are provided in Figure 13. The detailed configuration and boundary conditions are provided in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF], and we recall that we consider the case of a nominal incident shock Mach number of 1.33 in ambient air (with γ = 1.4) at 0.95 bar pressure and that the variables initialization consists in the pre-shock values (ρ, u, v, w, p) = (1.1175, 0.0, 0.0, 0.0, 95000.0) and the post-shock ones (ρ, u, v, w, p) = (1.7522, 166.3435, 0.0, 0.0, 180219.75) leading to conditions of [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF] at temperature T = 296.15K. In Figure 14, we plot the magnitude of the density gradient computed with the MOOD P 2 method equipped with the [PAD → DMP → u2] detection process at different times of the simulation in order to give an overview of the physical phenomena. We emphasize that the instabilities along the cylindrical wall are very well captured. Finally in Figure 15, we provide a zoom of the final solution on the created instabilities which perfectly match the experimental results of [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF].

The explosion problem

We consider the so-called explosion problem [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, 3 rd revision[END_REF] given by a gas initially at rest in the unit cube where a quarter of the ball of radius r c = 0.4 centered at the origin has a density ρ b = 1.0, a pressure p b = 1.0 whereas the exterior is characterized by ρ e = 0.125, p e = 0.1. The domain is partitioned into 20 3 hexahedral cells for which each hexahedron is split into 6 pyramids leading to a mesh of 48000 pyramids. Simulations are carried out till the final time t final = 0.25. A reference solution has been computed with a two-dimensional cylindrical staggered numerical Lagrangian scheme [START_REF] Loubère | Staggered Lagrangian discretization based on cell-centered Riemann solver and associated hydrodynamics scheme[END_REF].

We report in Figure 16 the density approximations in function of the radius for a classical MUSCL scheme [START_REF] Park | Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids[END_REF], the MOOD-P 2 , MOOD-P 3 , MOOD-P 5 methods equipped with the [PAD → DMP → u2] detection and the reference solution. Note that we use the same type of representation than for the previous test cases by slicing the radius in 100 uniform cells.

The solution shape is well reproduced by all methods and the higher the polynomial degree is the sharper the contact discontinuity and the shock wave are. In a three-dimensional context with discontinuous solutions, the MOOD-P 3 method seems to be the right balance between accuracy and cost. The slight improvement gained by the MOOD-P 5 compared to MOOD-P 3 may not justify the computational over-cost (see further). We also notice that the head of the rarefaction wave is badly resolved by the MUSCL method whereas the MOOD-P 2 and especially the MOOD-P 3 method give accurate approximations. To compare the different detection strategies, we present in Figure 17 the final solutions obtained by the MOOD method with the PAD alone and the [PAD → DMP → u2] detection processes using P 5 polynomial reconstructions. Note that contrary to previous figures, we plot the density values for all cells by associating them with the radius corresponding to the cell centroid. As expected, the PAD detection process does not damp spurious oscillations close to the shock wave (see the zoom panel) and extra oscillations are also visible on the head of the rarefaction. We recall that the numerical approximation using the PAD detection process is the most accurate one on smooth solutions since only the physical admissibility of the solution is required so that few numerical diffusion is produced. On the opposite, the [PAD → DMP → u2] detection process damps the oscillations close S. DIOT to the shock and to the head of rarefaction but also maintain a very good accuracy with a slight non monotonic behavior. Finally we provide in Tables III the computational cost of the MOOD method for this test case in this particular configuration when running on a single core of the three following machines (using -O3 flag for gfortran compiler):

M1: server with two Intel Xeon E5335 (4 cores) @ 2.00Ghz, 8MB of L2 Cache, 16GB of RAM M2: laptop with Intel Core2Duo P7550 (2 cores) @ 2.26GHz, 3MB of L2 Cache, 8GB of RAM M3: desktop with Intel Core i5 2500 (4 cores) @ 3.30GHz, 6MB of L2 Cache, 8GB of RAM Note that the same three machines have been used in [START_REF] Diot | Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF] to assess the computational cost of the MOOD method for two-dimensional geometries. We first observe that the memory storage doubles when the polynomial degree is increased by one: 0.4 for P 2 , 0.8 for P 3 , 1.6 for P 4 (not presented in the table) and 3.0 for P 5 . Notice that the memory consumption is very low since only two reconstruction pseudoinverse matrices (for P 2 and P dmax ) per cell are effectively stored. The CPU cost increases by a factor about 2.4 from P 2 to P 3 and about 2.7 from P 3 to P 5 .

By extrapolation of these results we estimate the cost of the MOOD method for larger meshes. As instance for one million cells mesh and 1000 time steps the method cost should be: Consequently simulations with nowadays sequential computers with a one million cells mesh (assuming one thousand time steps) can be obtained for about one day of computation with MOOD-P 3 method. The MOOD method is thus a very competitive very high-order finite volume method, and these results shall be improved by an efficient parallelization.

Interaction of a shock wave with a quarter of cone

To conclude the numerical tests section, we run the test case named interaction of a shock wave with a quarter of cone with the 4 th -order MOOD-P 3 method equipped with the [PAD → DMP → u2] detection process. This 3D extension of the so-called interaction of a shock wave with a wedge has been proposed in [START_REF] Dumbser | Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems[END_REF] as instance.

The domain consists in a quarter on cylinder of radius R = 2.25 centered on the Ox axis which covers the interval [-1.1; 3.0] in the x-direction. Note that three modifications have been made in comparison to [START_REF] Dumbser | Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems[END_REF] in order to reduce the computational cost: the test is run on a quarter of cylinder instead of a half one, the initial interface is placed at x = -0.2 instead of x = -1.0 and the domain covers in the x-direction the interval [-1.1; 3.0] instead of [-1.5; 3.0]. Finally the mesh obtained by the free mesher Gmsh contains 1161854 tetrahedra in three refinement zones and exactly matches the initial interface, see top of Figure 19.

We recall that the circular cone under consideration is such that its length is 1, its tip and foot radii are 0.02 and 0.5 respectively while its tip is placed at the origin. Moreover wall boundary conditions are prescribed everywhere except from the top and bottom of the quarter of cylinder where the exact solution according to the Rankine-Hugoniot conditions is imposed. At last the initial pre-and post shock conditions are given by (ρ, u, v, w, p) = (2.122, 0.0, 0.0, 0.0, 1.805) and (ρ, u, v, w, p) = (1.4, 0.0, 0.0, 0.0, 1.0) respectively with γ = 1.4 and the final time is chosen such that it corresponds to the final time of [START_REF] Dumbser | Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems[END_REF].

In Figure 18, we propose numerical Schlieren-type images on the solution in the Ox -Oy and Ox -Oz planes. We remark that the symmetry is very well conserved since both images are almost identical and that all waves that are present in results of [START_REF] Dumbser | Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems[END_REF] are also resolved here although much less cells (more than 3.5 times less) are considered. This proves that the MOOD method performs very well on 3D unstructured meshes. Finally on bottom of Figure 19, we provide a 3D view for which isosurfaces have been chosen to represent the principal waves in the whole domain. It is thus clear that the method properly reproduces the cylindrical symmetry even on this fully unstructured 3D tetrahedral mesh.

CONCLUSION

In this paper we have proposed the three-dimensional extension of the so-called MOOD method [START_REF] Clain | A high-order finite volume method for hyperbolic systems: 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]. The Multi-dimensional Optimal Order Detection expression refers to an original way of determining the optimal local polynomial degree to be used in the reconstruction step of a classical high-order unlimited scheme. To each cell corresponds a polynomial reconstruction for which we a posteriori determine the degree according to given criteria (positivity as instance) against which we test each candidate solution. The detection criteria is based on a relaxed version of the discrete maximum principle (DMP) associated with a so-called u2 detection procedure which analyses the numerical curvatures in the neighborhood of a DMP violating cell and determine if the underlying function is regular or not. In the latter case the polynomial degree in the associated cell is decremented and the solution is locally recomputed. We have detailed the numerical method for three-dimensional unstructured meshes and improved the detection criteria both for the advection equation and the Euler system. Moreover some optimizations for the three dimensional case have been provided to significantly improve the efficiency of the method.

The MOOD method has been implemented on several kinds of unstructured meshes with P k polynomial reconstructions (k varying from 1 to 5). We have provided some sanity checks with simple configurations and performed more advanced full three-dimensional tests to assess the ability of the MOOD method to accurately capture waves on real unstructured meshes. For the scalar convection equation with a regular initial shape the method gives an effective high-order of accuracy corresponding to the optimal one and we have shown that spurious oscillations are damped when discontinuous profiles are convected. The results for unidirectional problems for the Euler system with three-dimensional unstructured meshes show that small scaled structures are captured while shock waves are resolved within few cells. For the isentropic vortex test case extended to the three-dimensional context with non-trivial exact solution, effective high-orders of accuracy are measured and optimal orders are reported for P 2 , P 3 and P 5 polynomials. We prove that the MOOD method is able to capture the realistic physics of the impact of a shock wave on a cylindrical cavity on a non trivial mesh made of a mix of triangular and quadrangular prisms. At last, the three-dimensional explosion problem has been carried out to show the improvement gained with the use of high-order MOOD methods and the slight numerical diffusion generated by the u2 detection process which enables to prevent spurious numerical oscillations from appearing. We have also provided the solution computed with the PAD detection process alone to support the intrinsic positivity-preserving property of the MOOD method and measures of the CPU cost to underline that the MOOD method is effective on nowadays personal computers. Finally the interaction of a shock wave on a quarter of cone with the 4 th -order MOOD method proves that the MOOD method provides a very good reproduction of the physics on a unstructured non-regular 3D mesh of 1.1 millions of tetrahedra.

In a near future, we plan to adapt the MOOD within an ADER technique to avoid the multiple time steps of the Runge-Kutta approach and overcome the third-order accuracy restriction. Furthermore although the MOOD method significantly reduces the necessary computational resources (CPU and memory storage), a parallelized version is of crucial importance to treat huge size simulations. Finally the application of the MOOD method to more complex physics (multimaterial, multi-phase, etc.) is also an important challenge that has to be tackled.

6 3 Detection Criteria 7 3 . 1

 631 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 The MOOD concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Advection equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 Hydrodynamics Euler system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Implementation and optimizations . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4 Numerical results 14 4.1 Advection equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2 Euler system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 .

 2 Figure 2. Illustrations for index sets ν(i) (left) and ν(i) (right) in 3D.

Figure 3 .

 3 Figure 3. Flowchart of the MOOD algorithm: Recon. stands for polynomial reconstruction, Eval. for highorder evaluations at quadrature points, H-O Update for high-order update of the solution and Det. Proc. for detection process.

Figure 4 .

 4 Figure 4. Algorithmic view of the [DMP → u2] detection process for the advection equation.

Figure 5 .

 5 Figure 5. Algorithmic view of the [PAD → DMP → u2] detection process for the Euler system.

Figure 6 .

 6 Figure 6. Triple sine translation: convergence curves for L 1 (middle) and L ∞ (bottom) errors for series of hexahedral (left) and hexahedral/pyramidal (right) meshes. Examples of such meshes are given on top line.

Figure 7 .

 7 Figure 7. Initialization of the H-like shape rotation problem. Top left: interior view of the tetrahedral mesh. Top right: initialization of the H-like shape (isosurface 1/2, rotation axis is the red line). Bottom left: cut plane z = 1/2. Bottom right: extruded initial values from the cut plane.

Figure 8 .

 8 Figure 8. Results of the H-like shape rotation problem on the cut plane z = 1/2 for the unlimited P 3 and P 5 schemes (top line) and for the MOOD P 3 and P 5 methods (bottom line). The highlighted green cells correspond to values below 0 or above 1.

•

  Sod: (ρ, u, p) L = (1.0, 0.0, 1.0) and (ρ, u, p) R = (0.125, 0.0, 0.1), t final = 0.2, • Lax: (ρ, u, p) L = (0.445, 0.698, 3.528) and (ρ, u, p) R = (0.5, 0.0, 0.571), t final = 0.13.

Figure 10 .

 10 Figure 10. The MOOD-P 3 density results are displayed for the Sod (left) and Lax (right) problems on tetrahedral mesh vs the exact solution (red line).

Figure 11 .

 11 Figure 11. Results for the Shu-Osher (left) and Blastwave (right) problems on pyramids. MOOD-P 3 density results are displayed on the left and right columns respectively vs the reference solution (red line).

Figure 12 .

 12 Figure 12. Isentropic vortex in motion: convergence curves for L 1 (middle) and L ∞ (bottom) errors for series of hexahedral (left) and pyramidal (right) meshes. Examples of such meshes are given on top line.

Figure 13 .

 13 Figure 13. Impact of a shock on a cylindrical cavity: details of the mesh containing 202254 triangular and quadrangular prisms.

Figure 14 .

 14 Figure 14. Impact of a shock on a cylindrical cavity: magnitude of the density gradient at different times from left to right and top to bottom.

Figure 15 .

 15 Figure 15. Impact of a shock on a cylindrical cavity: Zoom on the created instabilities at final time.

Figure 16 .

 16 Figure 16. Density results for the explosion problem in 3D. Comparison between a classical MUSCL method and the MOOD-P 2 , MOOD-P 3 and MOOD-P 5 methods with [PAD → DMP → u2] detection process on tetrahedral mesh. The straight line corresponds to the reference solution.

Figure 17 .

 17 Figure 17. Density results for the explosion problem in 3D obtained with the MOOD-P 5 method. Comparison between the PAD alone and [PAD → DMP → u2] detection process. The straight line corresponds to the reference solution and the symbols represent mean values of all cells.

Figure 19 .

 19 Figure 19. Interaction of a shock wave with a half cone: on top, view of the interior of the tetrahedral mesh with the different zones of refinement; on bottom, isosurfaces corresponding to the principal waves.

  Compute polynomial reconstruction of degree d i , ∀i ∈ E el . 2. Compute FacePD d ij and d ji and evaluate high-order approximations at quadrature points on face f ij , ∀j ∈ ν(i), ∀i ∈ E el . 3. Compute candidate solution mean values through unlimited scheme (3), ∀i ∈ E el . 4. Detection process: decrement CellPD of cells where solution is not acceptable. 5. Stop if the solution is acceptable else go to stage 1.

Table II .

 II L 1 , L 2 and L ∞ errors and convergence rates for the isentropic vortex problem with the MOOD-P 2 , MOOD-P 3 and MOOD-P 5 methods. Top lines: hexahedral meshes. Bottom lines: pyramidal meshes.

			MOOD on hexahedra
	Deg.	Cell nb	L 1 error		L 2 error		L ∞ error
	P 2	20 × 20 × 4	1.2149125472e-2	-	3.1900165943e-2	-	2.4441600308e-1	-
		40 × 40 × 4	3.2245099895e-3 1.91 8.0410260168e-3 1.98 5.8071095773e-2 2.07
		60 × 60 × 4	1.2274201731e-3 2.38 3.0952199533e-3 2.35 2.2241472168e-2 2.37
		80 × 80 × 4	5.8449248920e-4 2.57 1.4891720330e-3 2.54 1.0333244490e-2 2.66
		120 × 120 × 4 1.8870676632e-4 2.78 4.8726522430e-4 2.75 3.2966795901e-3 2.82
	Expected order		3		3		3
	P 3	20 × 20 × 4	3.8426301161e-3	-	9.2866634542e-3	-	6.3860302401e-2	-
		40 × 40 × 4	7.2909293814e-4 2.39 1.6463735019e-3 2.49 1.3928173243e-2 2.19
		60 × 60 × 4	1.4537313954e-4 3.97 3.4725838178e-4 3.84 2.5316808689e-3 4.20
		80 × 80 × 4	4.3014601762e-5 4.23 1.1403422006e-4 3.87 8.9234157884e-4 3.62
		120 × 120 × 4 7.7653485653e-6 4.22 2.0827671718e-5 4.19 1.7793186459e-4 3.98
	Expected order		4		4		4
	P 5	20 × 20 × 4	2.8991068920e-3	-	4.8543664172e-3	-	3.2038381504e-2	-
		40 × 50 × 4	2.2151699683e-4 3.71 5.5851141683e-4 3.12 6.2194475329e-3 2.36
		60 × 60 × 4	2.8610132561e-5 5.04 7.5286576723e-5 4.94 4.9068468256e-4 6.26
		80 × 80 × 4	5.4168534310e-6 5.78 1.5519206048e-5 5.49 1.5955744462e-4 3.90
		120 × 120 × 4 4.0840597698e-7 6.38 1.1795674119e-6 6.36 1.0709587465e-5 6.66
	Expected order		6		6		6
			MOOD on pyramids
	Deg.	Cell nb	L 1 error		L 2 error		L ∞ error
	P 2	20 × 20 × 24	3.3660908651e-3	-	8.2020368268e-3	-	5.8966752971e-2	-
		40 × 40 × 24	6.0800306087e-4 2.47 1.4780372369e-3 2.47 1.2917288297e-2 2.19
		60 × 60 × 24	1.9831385885e-4 2.76 5.0256415975e-4 2.66 3.4489695638e-3 3.25
		80 × 80 × 24	7.9096059248e-5 3.19 2.0028509695e-4 3.19 1.3642153624e-3 3.22
	Expected order		3		3		3
	P 3	20 × 20 × 24	8.8005733635e-4	-	2.0405839361e-3	-	2.2060839273e-2	-
		40 × 40 × 24	6.4460987694e-5 3.77 1.4763173293e-4 3.78 1.3204082077e-3 4.06
		60 × 60 × 24	1.2809782719e-5 3.98 2.9223354775e-5 3.99 2.8960192576e-4 3.74
		80 × 80 × 24	4.0713141263e-6 3.98 9.3121356054e-6 3.97 8.6899534957e-5 4.18
	Expected order		4		4		4
	P 5	20 × 20 × 24	3.7944742185e-4	-	9.8016940506e-4	-	2.0273963181e-2	-
		40 × 40 × 24	9.7451977113e-6 5.28 2.4664732108e-5 5.31 2.4540502338e-4 6.36
		60 × 60 × 24	8.0304771455e-7 6.15 2.0569058735e-6 6.12 2.0022941712e-5 6.18
		80 × 80 × 24	1.2520658320e-7 6.46 3.1436119294e-7 6.53 3.2667224251e-6 6.30
	Expected order		6		6		6

Copyright c 2012 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids (2012) Prepared using fldauth.cls DOI: 10.1002/fld S. DIOT

Table III

 III

. CPU time in microseconds per iteration per cell and memory storage in Gigabytes for the MOOD-P k methods (k = 2, 3, 5) with the [PAD → DMP → u2] detection process on three different computers.

•

  MOOD-P 2 is 66000 seconds on M1, that is to say ∼ 18 hours (∼ 16 hours on M2 and ∼ 8.3 hours on M3) with about 8 Gb of memory storage, • MOOD-P 3 is 163000 seconds, ∼ 2 days on M1 (∼ 1.5 day and ∼ 19 hours on M2 and M3) with about 16 Gb of memory storage, • MOOD-P 5 is 439000 seconds, ∼ 5 days on M1 (∼ 4.5 days and ∼ 2.1 days on M2 and M3) with about 62 Gb of memory storage.

Copyright c 2012 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids (2012) Prepared using fldauth.cls DOI: 10.1002/fld

S. DIOT