Christian Tenaud
email: christian.tenaud@limsi.fr

Olivier Roussel
email: olivier.roussel@coria.fr

Linda Bentaleb
email: linda.bentaleb@yahoo.fr

Unsteady compressible flow computations using an adaptive multiresolution technique coupled with a high-order one-step shock-capturing scheme

Keywords: adaptive mesh refinement, multiresolution, finite volume, shock-capturing scheme, compressible flow, shock-wave PACS: 47.11.-j, 47.11.Df, 47.11.St, 47.40.-x, 47.40.Nm

published or not. The documents may come

Unsteady compressible flow computations using an adaptive multiresolution technique coupled with a

high-order one-step shock-capturing scheme

Introduction

In the high speed flow regime, many aerodynamic configurations involve interactions between shock waves and turbulence such as, for instance, within air intakes, compressor or turbine configurations where shock wave/turbulent shear layer (e.g. boundary layer) interactions occur. An accurate prediction of such interactions is of importance in effective design of transonic or supersonic vehicles since they greatly affect the aerodynamic loads. At the present time, it is commonly admitted that advanced numerical simulations (mainly LES) are powerful tools for accurate predictions of shock wave turbulent shear layer interactions, including large-scale flow phenomena [START_REF] Ducros | Large eddy simulation of shock/turbulence interaction[END_REF][START_REF] Knight | Advances in CFD prediction of shock wave turbulent boundary layer interactions[END_REF]. In these approaches, the quality of the solutions depends not only on the capability of the numerical scheme associated with the sub-grid modelings in LES but also on the ability of the computational grid to capture the governing dynamical process. In fact, when dealing with shock waves, LES computations must however use numerical schemes which can both represent small scale structures with the minimum of numerical dissipation, mainly to minimize the interaction with the sub-grid scale model, and capture discontinuities with robustness [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF][START_REF] Titarev | Ader: Arbitrary high order Godunov approach[END_REF]. Nevertheless, some phenomena could not be accounted by sub-grid modeling and accurate schemes coupled with locally very fine grid are needed to recover a high quality of the solution. For instance, according to the theoretical developments in the Linearized Interaction Approximation, the shock-wave/turbulence interaction phenomenon requires the correct prediction of the shock wave deformation occurring at small scales. These small scale shock deformations which could not be accounted by LES modeling, need locally very fine grid. The production of vorticity through baroclinic effect is also a phenomenon largely encountered in real flow physics that could not be accounted by sub-grid modeling and needs accurate numerical scheme on grid tightened in the production regions. These examples show that it is necessary to dynamically refine the grid locally, in the regions where the unsteady phenomena occur. It then motivates the introduction of self-adaptive discretization, as the solution may be over-resolved in large subsets of the computational domain when using equidistant fine grids. Therefore, to be efficient in terms of CPU time and memory usage, a mesh refinement method must be employed to save grid points in smooth regions and to concentrate them in the regions where phenomena (discontinuity, vor-Coupled time-space schemes are usually developed following the Lax-Wendroff approach. Among them, the One Step (OS) schemes have been developed, first for 1D (linear and non-linear) scalar equations, and then extended to multidimensional systems of non-linear equations (Daru and Tenaud [START_REF] Daru | Evaluation of TVD high resolution schemes for unsteady viscous shocked flows[END_REF][START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF]). Such schemes have a minimal stencil, and optimal non-oscillatory conditions, based on Monotonicity-Preserving constraints, can easily be implemented. These accurate numerical schemes offer a compromise between high accuracy in smooth regions and an efficient shock capturing technique. They provide very accurate results, which compare well to high-order separate time-space classical schemes, at a lower cost [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF].

Besides the numerical scheme, the quality of solutions also depends on the capability of the computational grid to capture the governing dynamical mechanisms. In that sense, Adaptive techniques for problems exhibiting locally steep gradients or shock-like structures have been developed since the end of the 1970s. Historically, adaptive methods like Multi-Level Adaptive Techniques (MLAT) (Brandt [9]) or Adaptive Mesh Refinement (AMR) (Berger et al [START_REF] Bell | Three-dimensional adaptive mesh refinement for hyperbolic conservation laws[END_REF][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]) were the first to achieve this goal, using a set of locally refined grids where steep gradients of high truncation errors are found. However, the data compression rate is high where the solution is almost constant, but remains low where the solution is regular. To overcome this difficulty, adaptive multiresolution methods, based on Harten's pioneering work [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF], have been developed for 1D and 2D hyperbolic conservation laws (Cohen et al [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF], Müller et al [START_REF] Gottschlich-Müller | Adaptive finite volume schemes for conservation laws based on local multiresolution techniques[END_REF]). They have then been extended to 3D parabolic problems (Roussel et al [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF]). First simulations of 3D supersonic flows in the laminar regime using adaptive multiresolution methods were performed by Bramkamp et al [START_REF] Bramkamp | H-adaptive multiscale schemes for compressible Navier-Stokes equations -polyhedral discretization, data compression and mesh generation[END_REF][START_REF] Bramkamp | An adaptive multiscale finite volume solver for unsteady and steady flow computations[END_REF], with separate RK/ENO time-space discretizations. It has been shown in these papers that a high compression rate can be reached for solutions with inhomogeneous regularity. For an overview on adaptive multiresolution techniques, we refer to the books of Cohen [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF] and Müller [START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF].

This paper aims at evaluating in practical situations the capability of the multiresolution adaptive technique coupled with a one-step shock capturing scheme to recover elementary physical mechanisms by achieving gains in both CPU time and memory usage compared to single grid computations. Numerical simulations are conducted on both inviscid and viscous compressible flows with high frequency oscillations in the transonic and supersonic regimes. As far as there could exist a competition between the discretization error of the scheme and the perturbation error introduced by the MR technique, we use several approximation orders of the OSMP scheme on several grids. The question that arises is: is it better to employ a low order (at least 2nd order) scheme on a very refined grid than use a high order scheme on a coarse grid? Through comparisons with 2nd and 3rd order schemes, we then explore the efficiency of a high order scheme coupled with the MR technique. We then propose an appropriate MR threshold parameter value that ensures accurate results, while achieving drastic gains on the CPU time and memory usage.

The paper is organized as follows. After a brief review of the governing equations ( § 2), we present in section 3 the numerical approach used in this work: the so-called OSMP scheme (based on coupled time and space integration with MP constraint). Section 4 is dedicated to detail the multiresolution procedure. The evaluation of the method is then presented in section 5 on several numerical results, 2D and 3D inviscid configurations and 2D configurations where viscous effects are present. Finally we conclude and present perspectives for future works.

Governing equations

We consider the Navier-Stokes equations expressed in dimensionless form using Cartesian coordinates:

w t + ∇ • f E (w) -f V (w, ∇w) = 0, (1)
where w = (ρ, ρu, ρE) t is the vector of the conservative variables, using the classical notations, and f E (w) and f V are the Euler and the viscous fluxes respectively:

f E =         ρ u ρ u ⊗ u + P γM 2 0 I ρ u E + u P γM 2 0         , (2)
f V =        0 σ u • σ + Ψ        , (3)
with the strain rate tensor

σ = µ Re ∇u + ∇ t u - 2 3 ∇ • u I ,
and the heat flux

Ψ = µ (γ -1) Re Pr M 2 0 ∇T.
In addition, a perfect gas law is needed:

P γ M 2 0 = (γ -1) ρ E - 1 2 ρ u • u , (4)
T = P ρ , (5)
with ρ the fluid density, u the velocity vector, P the static pressure, T the static temperature, E the total energy per unit of mass and µ the dimensionless dynamic viscosity.

These equations are written in dimensionless form using the reference values of the density (ρ 0 ), the velocity (v 0 ), and the length scale (L 0 ). The Reynolds number is based on the reference values: Re = ρ 0 v 0 L 0 /µ(T 0 ). The Mach number is M 0 = v 0 / γRT 0 with R = 287 J.kg -1 .K -1 (for air) the gas constant.

For simplicity, this study considers an ideal gas (air) with constant specific heat ratio γ = 1.4 and constant Prandtl number Pr = 0.73. The dimensionless dynamic viscosity µ is expressed with the Sutherland's law

µ = T 3 2 1 + T s T + T s , (6)
where the Sutherland's constant T s = 110.4 T 0 .

Coupled time and space integration

The resolution of the Navier-Stokes equations (1, 2, 3) is based on a finite volume approach on Cartesian grids. Let us consider a computational domain (Ω) partitioned in subsets (V j as control volumes) dense in Ω.

We denote by wn j the cell-averaged value of the discrete quantity w(x, t), estimated at a grid point x j = j δx centered in the control volume (V j ) and at time t n = n δt, δt and δx being, respectively, the time step and control volume size. The cell-averaged definition is the following:

wj (t) = 1 |V j | V j w(x, t) dx, (7)
where

|V j | = V j
dx is the measure of the control volume.

The Navier-Stokes equations (1, 2, 3) are solved using an operator splitting approach: as a one-step approach is used to treat the convective terms, convection and diffusion operators are treated separately to recover at least the second order accuracy. Let us note that the operator splitting procedure could allow us to implement an implicit integration of the viscous fluxes to relax the Von-Neumann constraint on the time step. We implement the operator splitting as follows:

w n+1 j = L E δt/2 • L V δt/2 • L V δt/2 • L E δt/2 w n j , (8)
where L E δt is the discrete approximation of the operator of convection; L E δt :

w t = -∇ • f E (w) and L V
δt the discrete approximation of the diffusion operator; L V δt :

w t = +∇ • f V (w, ∇w) .
The integration of the viscous part of equation ( 1) is performed using a classical 2nd-order centered scheme coupled with a 2nd-order Runge-Kutta time integration. Regarding the convective fluxes, we used the high-order coupled time and space scheme proposed in [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF][START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF][START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF] that gives very accurate results on various laminar test-cases at very low-cost in terms of CPU time [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF][START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF][START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF].

High-order one-Step scheme for convection

The integration of the convective terms is performed using a one-step (OS) scheme which is of the Lax-Wendroff type. The construction of such schemes can also be obtained by correcting the successive error terms of the modified equations to increase the order of accuracy of the schemes. In this way, we obtain high-order accurate schemes relative to both time and space. The implementation of these schemes is very simple because the increase of accuracy is obtained via a change of an accuracy function (noted Φ o , where o is the order of accuracy) applied to a classic incremental scheme (see [START_REF] Daru | Evaluation of TVD high resolution schemes for unsteady viscous shocked flows[END_REF][START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF][START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF] for more details on the procedure).

As far as the 1D system of Euler equations (w t + ∇ • f E (w) = 0) is considered, the One Step (OS) scheme reads in the Finite-Volume framework:

w n+1 j = w n j - δt δx (F o j+1/2 -F o j-1/2 ), (9)
where F o j+1/2 is the oth-order accurate numerical flux of the scheme at the cell interface (j + 1/2), which can be written:

F o j+1/2 = F Roe j+1/2 + 1 2 k (ψ o r) k,j+1/2 . (10)
For clarity, the superscript n has been omitted in the expression of the fluxes.

As far as non-linear system is concerned, while the accuracy order cannot be maintained easily, a simple local linearization on the eigenvector basis at the interfaces of the cells has been preferred to the expansive Cauchy-Kowalewski procedure which might be required along the derivation of the successive modified equations. Here, we denote by λ k and r k , respectively, the eigenvalues and right eigenvectors of the Roe-averaged jacobian matrix A = df dw j+1/2 evaluated on the cell interface (j + 1/2). The k-th Riemann invariant is then evaluated by

δα k,j+1/2 = r k • w n j+1 -w n j . (11)
Here, F Roe j+1/2 denotes the first order Roe flux, defined as

F Roe j+1/2 = 1 2 (f j + f j+1 ) - 1 2 k (δ|f |r) k,j+1/2 , (12)
with

δ|f | k,j+1/2 = |λ| k,j+1/2 δα k,j+1/2 .
The function (ψ o k,j+1/2 ) drives the order of accuracy and extends to the oth order the accuracy of the basic first-order Roe scheme. A general form of accuracy functions ψ can be found in [START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF]. Here, computations are performed using the 2 nd , 3 rd and 7 th -order accuracy functions.

The present unlimited One-Step 7-th order accurate scheme (called OS7) requires a stencil of only nine points to achieve the accuracy order in both time and space. It has a classical CFL stability condition 0 ≤ |ν k | ≤ 1. Note that this type of scheme yields the exact solution for CFL number equal to 1, in the scalar case.

To prevent spurious oscillations in the vicinity of strong gradient regions, an ad hoc limiting process must be employed. Total Variation Diminishing (TVD) schemes are generally considered to be well suited for the capture of sharp discontinuities without oscillations. Nevertheless, TVD constraints are known to clip the extrema, which appears as a serious drawback in a limiting procedure. To avoid this loss of accuracy near extrema, it is necessary to satisfy the Monotonicity Preserving (MP) criteria, introduced by Suresh and Huynh [START_REF] Suresh | Accurate monotonicity-preserving schemes with runge-kutta time stepping[END_REF], that enlarge the TVD intervals to provide room for the numerical flux to maintain an accurate value. Daru and Tenaud [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF] derived a Monotonicity Preserving version of the present scheme (named OSMP) that preserves accuracy near extrema. In particular, the original MP constraints of Suresh and Huynh [START_REF] Suresh | Accurate monotonicity-preserving schemes with runge-kutta time stepping[END_REF] has been recasted in the TVD framework without CFL restriction to generalize the MP conditions in terms of flux limiting. The MP conditions that preserve accuracy can be expressed directly as constraints on ψ o functions [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF]. The MP constraints are written for each characteristic field and can be found in [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF][START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF].

By taking the limiting function ψ o-M P k in [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF], the resulting scheme will be high-order accurate almost everywhere, except around discontinuities where it becomes first order accurate as is the case of all TVD schemes. Let us note that the essential difference is that MP constraints act for non-monotone data and so far are TVD for monotone data such that the scheme is not oscillatory around discontinuities.

Extension to the multidimensional case

The extension in the multidimensional case is delicate as far as a one-step approach is used. In fact, we need to consider cross derivative terms that appear in the second and higher order terms, which are left uncontrolled if one applies a direction by direction MP correction to a Lax-Wendroff unsplit scheme. We also need to guarantee that the resulting scheme will be nonoscillatory. The simplest way to avoid the problem of cross derivatives and to recover the good properties of the one-dimensional scheme is to use a Strang directional splitting strategy [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF][START_REF] Strang | On the construction and comparison of difference schemes[END_REF], which is unfortunately only second order accurate. While the order of accuracy is lowered compared to the tensorial multistage approach, the OSMP scheme with the Strang algorithm provides with results with very small error level at low cost [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF].

In two dimensions, the splitting of the system of equations (1) writes

w n+2 j = L δx L δy L δy L δx w n j . (13)
Here L δx and L δy are discrete approximations of operators in each space direction. When directional operators do not commute, the second order accuracy is recovered every two time steps with however the symmetric property of the solution [START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF]. In the same way, a three dimensional splitting can be constructed that gives second order accuracy every six time steps (when spatial operators do not commute):

w n+6 j = (L δx L δy L δz ) (L δx L δz L δy ) (L δy L δz L δx ) (14) (L δy L δx L δz ) (L δz L δy L δx ) (L δz L δx L δy ) w n j ,
with L δz the operator discrete approximation in the third space direction.

Adaptive multiresolution method

Since the finite volume computation is based on cell-average values of the conservative variables, the cell-average multiresolution analysis [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF] is here employed. In the context of adaptive mesh refinement, the principle of the multiresolution analysis is to represent the data on a set of nested dyadic grids. We denote by l = 0, 1, ...., L the grid level from the coarsest (l = 0) to the finest (l = L). Here for simplicity, we only consider embedded Cartesian grids while the method has already been developed on both unstructured grids [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF] and general coordinate system [START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF]. At each grid level (l), we consider dense partitions V l j (control volumes) of the computational domain Ω. A tree data structure is used to encode the multiresolution analysis technique. An example of an embedded grid on a tree data structure is presented in 1-D on Figure 1. The root denotes the basis of the tree, nodes are elements of the tree and leaves are upper elements enhanced in dark in Figure 1. In d dimensions, a parent-cell at a level l has always 2 d children cells at the level l + 1. The cell-average value on each control volume at level l is defined as follows (for instance in 3-D):

wl i,j,k = 1 |V l i,j,k | V i,j,k w(x, t) dx, (15)
where

V l i,j,k := [2 -l i, 2 -l (i + 1)] × [2 -l j, 2 -l (j + 1)] × [2 -l k, 2 -l (k + 1)] with i, j, k ∈ 0, ..., 2 l-1 and |V l i,j,k | = V l i,j,k
dx is the measure of the control volume.

Cell-average multiresolution analysis

Relations between cell-average values at different levels are defined by the multiresolution analysis. They are obtained by projection and prediction operators to progress from a fine grid (at level l + 1) to a coarser one (at level l) and vice versa.

The projection operator (noted P l+1→l ) allows to compute the cell-average value of the solution on a node at grid level l from cell-average values of the solution on children nodes at grid level l + 1. As far as grids are nested, the projection operator is exact and unique [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF], given that cell-average values at two successive grid levels are related by:

wl j = p∈C l j |V l+1 j | |V l j | wl+1 p , (16)
where C l j denotes the index ensemble of the 2 d children at grid level (l + 1) of the cell V l j . Such a way, knowing the solution on leaves (at a grid level l), cell-average values can be calculated from grid level (l) down to the root-cell (l = 0). Table 1 Coefficients of centered linear polynomial interpolations [START_REF] Postel | Approximations multiéchelles. Ecole de printemps de mécanique des fluides numérique[END_REF] The prediction operator (noted P l→l+1 ) maps wl to an approximation ŵl+1 of wl+1 . In contrast with the projection operator, there is an infinite number of choices for the definition of P l→l+1 . Nevertheless, in order to be applicable in a graded tree structure, it needs to be local, i.e. based on an interpolation using the s nearest neighbors in each direction, and consistent with the projection operator, i.e. P l+1→l • P l→l+1 = Id. The former property implies that the s-nearest neighbors must include the parent-cell [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF][START_REF] Postel | Approximations multiéchelles. Ecole de printemps de mécanique des fluides numérique[END_REF]. Here to predict the approximated value, we used centered linear polynomial interpolation:

order (r) s ξ 1 ξ 2 ξ 3 1 0 0 0 0 3 1 -1 8 0 0 5 2 -22
           ŵl+1 2j = wl j + s p=1 ξ p wl j+p -wl j-p , ŵl+1 2j+1 = wl j - s p=1 ξ p wl j+p -wl j-p , (17)
where ξ p are coefficients of the linear polynomial interpolation of order r = 2 s + 1, given in table 1 for s ≤ 3. Since a Cartesian mesh is used, extension to multidimensional polynomial interpolations is easily obtained by a tensorial product of the 1-D operator [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF][START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF].

Multiresolution transform : The error of prediction at a grid level l is estimated by evaluating "details" ( dl j ) defined as the difference between numerical solutions and interpolated values at this grid level l:

dl j = wl j -ŵl j . (18)
Depending on the accuracy order (r) of the centered polynomial interpolation, details recover null values for smooth solutions with locally bounded r-th order derivatives [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF]. Moreover, when the solution is regular, it was shown [START_REF] Bramkamp | An adaptive multiscale finite volume solver for unsteady and steady flow computations[END_REF][START_REF] Cohen | Multiresolution finite volume schemes on triangles[END_REF][START_REF] Postel | Approximations multiéchelles. Ecole de printemps de mécanique des fluides numérique[END_REF] that details decay with a rate at least of 2 -l . Therefore, for smooth solutions, the higher the grid level, the smaller the details. On the opposite, details recover significantly high values in regions where singularities of the solution occur. Thanks to the consistency assumption [START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF][START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF], the sum of details on the children of a parent cell is equal to zero [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF]. Therefore, in d dimensions, the knowledge of the 2 d children cell-averages of a given parent cell is equivalent to the knowledge of the parent cell-average and 2 d -1 details. If we note Dl the vector of all details at a grid level l, the vector of the solution ( Wl ) at a grid level l can be represented by the vector of the solution ( Wl-1 ) at the grid level l -1 and the vector of details (D l ), i.e. Wl ←→ ( Wl-1 , Dl ). Recursively on all L grid levels, one gets the so-called multiresolution transform [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF] that maps the vector of the solution on the finest grid (L-level) to the solution on the root-cell plus a series of vectors of details from grid level 1 to the highest level (L):

M : WL -→ ( W0 , D1 , . . . , DL ). ( 19
)
This multiresolution transform (M) must only be seen as a change of basis, keepping the number of degrees of freedom unchanged. This new data format is however more convenient for data compression because detail magnitudes are zero when solutions have locally bounded r-th order derivatives and decay with the grid level for smooth solutions. Nevertheless, to be competitive for computations, the multiresolution transform should not increase the computational complexity related to the number of floating point operations required by the algorithm. To decrease the computational complexity, one of the concepts is to adapt locally the mesh to the behavior of the solution through the use of a local grid refinement, which should decrease CPU time and memory requirements.

Local grid refinement

Indicator functions are needed to locate regions where the solution requires a mesh refinement for capturing small scale dynamical structures. Here, assuming that this basis [START_REF] Gottschlich-Müller | Adaptive finite volume schemes for conservation laws based on local multiresolution techniques[END_REF] is stable [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF], the indicator is based on the details because they decay with the grid level for smooth solutions and recover sufficiently large values in discontinuous regions. The idea is to only retain details that a measure (for instance, absolute value for scalar case) is greater than a threshold parameter (ε l ) on each grid level (l) and set others to zero. The threshold operation is then obtained through the following algorithm by taking into account a measure of details scaled by global maximum:

       if ∀j, | dl j | max j | wl j | < ε l =⇒ dl j = 0. otherwise dl j ∈ D (20)
Therefore, the computational mesh is formed by leaves that the normalized L 1 -norm of their details is greater than ε l . Leaves with details set to zero are discarded from the tree-data structure.

The threshold parameter (ε l ) has a crucial influence on the error introduced by the multiresolution procedure [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF], called the perturbation error, calculated as the difference between the finite-volume solution evaluated on the finest grid and the solution obtained on the refined grid using the multiresolution algorithm. If one wants that perturbation errors (here in L 1 -norm) recover the same magnitude (ε) on each grid level l, the threshold parameter (ε l ) must be prescribed at

ε l = 2 d(l-L) ε. ( 21
)
where d denotes the dimension. This has been first introduced by Harten [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF] as a heuristic hypothesis which has been corroborated by theoretical studies [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF][START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF]. That way, the value of ε drives the efficiency of the grid refinement. We will show later (see § 5) the influence of the ε value on the solution accuracy and on gains in both the memory usage and the CPU time.

As far as we deal with hyperbolic conservation laws, discontinuities could appear from a regular solution after a finite time. Therefore, to foresee the formation of discontinuities, the mesh can locally be refined by adding children of a leaf. Assuming that the forthcoming loss of regularity can be detected by detail values estimated on coarser grids and considering a leaf at a grid level l, if the following criterion:

dl j max j wl j ≥ 2 (2.s-1) .ε l , (22)
is satisfied, children at a grid level l + 1 are added to the tree.

In order to be graded, the tree must verify that, for each leaf at a level l, the prediction operator (P l→l+1 ) can be evaluated and the numerical fluxes can be calculated ensuring conservativity. Adjacent cells of level at least equal to l-1 must be added in each direction (Figure 1), the diagonal being included in multi-dimension. The number of adjacent cells depends on both the number (s) of the nearest neighbor cells of the polynomial interpolation and the number of points in the stencil of the numerical scheme. For instance, using the OSMP7 scheme, s + 2 nearest neighbor cells must be added in each direction including the diagonals in multi-D.

Conservative flux computation

To ensure conservativity, numerical fluxes on a cell face must have the same value whatever the cell considered adjacent to the face. Regarding a leaf of a locally refined grid, adjacent cells could however not exist at the same grid level as illustrated on Figure 2. As shown on this figure, the outgoing total flux from the right face of both V l+1 2i+1,2j and V l+1 2i+1,2j+1 should be balanced with the outgoing flux on the left face of cell V l i+1,j . In such a situation, the numerical fluxes are always computed at the highest grid level (l + 1) using virtual cells (named ghost or phantom cells elsewhere) that have been added at the grid level l + 1 (see figures 1 in 1D (dashed blue cells) and 2 in 2D (dashed cells)).

Since we use an explicit algorithm, let us note that the solution is not integrated on these virtual cells but just evaluated using the prediction operator from the solution at time n δt on the grid level l. To ensure a strict conservativity in the flux computation between cells at different grid levels, without increasing significantly the number of costly flux evaluations, the ingoing flux on the leaf at the grid level l is equal to the sum of the outgoing fluxes on the leaves at the grid level l + 1 (Figure 2) [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF], i.e.

F l i,j→i+1,j = F l+1 2i+1,2j→2i+2,2j + F l+1 2i+1,2j+1→2i+2,2j+1 .

Algorithm

In what follows, we present the adaptive MR algorithm that allows to code and let evolve a solution on a hybrid graded tree. For more details we refer to [START_REF] Postel | Approximations multiéchelles. Ecole de printemps de mécanique des fluides numérique[END_REF][START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF][START_REF] Roussel | An adaptive multiresolution scheme for combustion problems: application to flame ball-vortex interaction[END_REF][START_REF] Tenaud | Tutorials on adaptive multiresolution for mesh refinement applied to fluid dynamics and reactive media problems[END_REF]. From now on, we consider v L as the solution on the locally most refined cells, not necessarily on the finest grid.

(1) Initialization: (a) Initialize nested grids: we first build nested grids for grid-levels l = 0, 1, . . . , L from the coarsest to the finest grid. There are regular disjoint partitions (cells) V l j of an open subset Ω ⊂ R d (15), such that each V l j , j ∈ S l , is the union of a finite number of cells V l k , k ∈ S l+1 , where S l and S l+1 are ensembles of indices of consecutive embedded grids. (b) Initial solution: the initial solution is given either on cells of the finest grid V L j , j ∈ S L or on leaves of an initial graded tree, if it exists. (c) Threshold parameter: we also define a unique value of the threshold parameter ε. The threshold value for each grid level (l) is then given by equation ( 21); (2) Beginning of the time integration: (3) For time = 0 up to n time × δt (where δt is the time-step prescribed assuming a CFL condition), Do (a) MR algorithm: (i) Encoding: Knowing the cell-average values (v L j ) of the solution on the locally most refined leaves, propagate the solution from the finest grid to the coarser using the projection operator [START_REF] Daru | Numerical simulation of the shock wave / boundary layer interaction in a shock tube by using a high resolution monotonicitypreserving scheme[END_REF]. Then, encode the solution M L = M v L using the multiresolution transform [START_REF] Gottschlich-Müller | Adaptive finite volume schemes for conservation laws based on local multiresolution techniques[END_REF] through detail evaluations (18 and 17) from coarse grids to the finest. (ii) Predictive Harten's thresholding process: An effective data compression is performed through the thresholding process since some cells are literally discarded from the tree data structure. However, we have seen before that a graded tree data structure must be respected. Hence, before alleged useless cells are discarded, they are first marked using a logical function following (20, 21 and 22). (iii) Building a graded tree: marking cells or creating new cells that must belong to the tree to be graded; (iv) Pruning the graded tree: once all cells needed to form a graded tree have been marked, useless nodes are literally deleted from the tree data structure, with a prescribed tolerance ε. (v) Virtual leaves are added for conservativity constraints (b) Perform the time evolution of the solution:

(i) While a "true" (not virtual) leaf is selected, the time evolution of the solution is then performed by solving the discrete approximations of the Navier-Stokes equations (13 or 14). (4) End For (5) End of time integration (6) Save the solution at the final time (7) END Program

Numerical results

The chosen test-cases deal with interactions between wave propagations or shear layers and shock phenomena. It is a first step to recover very elementary phenomena that contribute to the turbulence production that could not be accounted by sub-grid modeling. The aim here is to perform computations in several configurations to assess the efficiency and the accuracy of the adaptive multiresolution method coupled with a high-order discretization scheme. In the following, the classical finite volume scheme on a regular finest grid is denoted by (FV) whereas the adaptive multiresolution scheme is denoted by (MR).

For the MR computations, we perform a parametric study on the threshold parameter (ε) to assess its optimal value, i.e. to find the largest value for which the discretization error of the numerical scheme is balanced with the accumulated threshold error, that is the so-called perturbation error [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF][START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF]. The performance analyses enable us to show if the choice of a MR algorithm is justified for such computations.

Convergence study on a 2D Euler case: advection of a strong vortex

The propagation at 45 o to the grid lines of a strong vortex at a supersonic Mach number is considered. This test-case has previously been treated in [START_REF] Balsara | Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy[END_REF]. The vortex is initially centered in a domain (x × y) ∈ [-5, 5] × [-5, 5]. The strong vortex is superimposed as a perturbation over an initial flow (ρ, p, u, v) = (1, 1, 1, 1), given by: (δu, δv) = 2π e 0.5(1-r 2 ) (-y, x); δT = -(γ -1) 2 8 γ π 2 e (1-r 2 ) ; δS = 0 where r 2 = x 2 + y 2 , S = p/ρ γ is the entropy and = 5 is the vortex strength. Periodic boundary conditions are imposed in both directions. Computations have been performed using several grids from the coarsest 32×32 using 5 gridlevels to the finest grid 512 × 512 using 9 grid-levels with a dyadic evolution.

The exact solution of this problem is just the passive convection of the vortex. The density contours of the initial solution are plotted on Figure 3. At t = 10, when the vortex recovers its initial position, the solution obtained with the MR approach with ε = 10 -3 and 9 grid levels (the finest grid is equivalent to 512 × 512 grid points), compares very well with the initial solution. The adapted grid is encoded on 5 grid levels (from 5 to 9 grid levels, Figure 4) and is well localized in the vicinity of the vortex. Errors in the L 1 -norm relative to the exact solution, reported in Figure 5, are calculated on the density at t = 10 when the vortex recovers its initial position. The OSMP7 scheme on a single grid (FV approach) recovers the fifth-order of accuracy in the L 1 norm (Figure 5-left). Although the same accuracy order as the MPWENO(r=5) (from Balsara and Shu [START_REF] Balsara | Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy[END_REF]) is recovered, the magnitude of L 1 -errors from the OSMP7 scheme is much lower than the one recovered by the MPWENO(r=5). This may mainly be attributed to the directional Strang splitting procedure [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF][START_REF] Strang | On the construction and comparison of difference schemes[END_REF], given that the boundary conditions are periodic. Besides, computations have been performed using the MR algorithm and 2nd, 3rd and 7th order OSMP schemes, with several threshold ε values ranging from 0 to 10 -5 . Let us remark that when ε = 0 the solution is computed everywhere Fig. 5. Convergence curves for FV approach (on the left) using the OSMP7 scheme compared to the MPWENO5 scheme [START_REF] Balsara | Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy[END_REF], and (on the right) for the MR approach using several ε values.

Fig. 6. Convergence curves for the MR approach using 2nd and 3rd order OSMP schemes for several ε values.

in the domain on the highest grid-level (the finest grid) with however the use of the MR algorithm. L 1 -errors relative to the exact solution for several ε values are reported on Figure 5-right for OMSP7, and on Figure 6 for both OMSP2 and OSMP3 schemes. For ε ≤ 10 -4 , the accuracy of the FV solution is recovered whatever the scheme. Let us note that, as the results for ε ≤ 10 -4 completely fit the FV results for OSMP2 and OSMP3 schemes, they are not reported on Figure 6 to allow reading. When ε = 10 -3 or 10 -2 , the accuracy of the solution is spoiled by the perturbation error induced by the multiresolution algorithm. Regarding the OSMP7 scheme, solutions for ε = 10 -3 or 10 -2 could however be considered as acceptable since error magnitudes are very low and at most decrease with a second-order of accuracy when ε = 10 -2 , for instance. Fig. 7. CPU-time versus memory usage for the MR approach using 9 grid levels (the finest grid is equivalent to 512 × 512 grid points) compared to the FV approach: on the left, the 7th order OSMP scheme using ε = 0, 10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 , (the same legend as in Figure 5 is used for ε values) and on the right, the 2nd and 3rd OSMP schemes using ε = 0, 10 -2 , 10 -3 , 10 -4 , 10 -5 .

This is not the case when the 2nd and 3rd order OSMP schemes are used, since error magnitudes are not low enough to consider solutions as reliable.

Even with very fine grids, the 2nd and 3rd order error magnitudes are much greater than the ones obtained on the coarsest grid with the OSMP7 scheme.

In Figure 7, we present the CPU time versus memory usage for the MR approach with 9 grid levels (the finest grid is equivalent to 512 × 512 grid points) to compute the solution with several ε values. The CPU time is reported w.r.t. the CPU time of the FV on the same highest grid-level. The memory usage evaluates the number of active leaves w.r.t. the number of grid points in the FV grid. The FV values are reported on these figures to exhibit performances (gains or losses) of the MR approach on CPU-time and memory usage relative to the value of the threshold parameter (ε, same legend as in Figure 5).

We can first remark that the MR algorithm, when applied everywhere on the finest grid (ε = 0), uses the same memory usage as the FV approach due to the one to one multiresolution transform [START_REF] Gottschlich-Müller | Adaptive finite volume schemes for conservation laws based on local multiresolution techniques[END_REF]. Gains on the memory usage are then reached when leaves can be discarded from the tree data structure. Secondly, when applied everywhere on the finest grid (ε = 0), the MR algorithm costs approximately 60 % to 70 % more than the one of the FV procedure. Consequently, gains on the CPU time can only be achieved when at least 40 % of leaves can be discarded from the tree-data structure. For the highest resolution (i.e. 512 × 512 grid points). This is obtained (Figure 7) when ε ≤ 10 -3 , for the OSMP7 scheme. The lower order OSMP schemes use less CPU time and memory usage than the corresponding OSMP7 scheme for the same ε value. Although efficiency seems to be slightly better than for the 7th order scheme, the middling accuracy recovered disqualifies the OSMP2 and OSMP3 schemes.

To go further in the analysis, gains and losses compared to the FV approach on both CPU-time and memory usage are presented in Table 2, for MR calculations using OSMP7 and several ε values on several grid levels. It is clear that, for grids with sufficiently resolved phenomena, gains on CPU-time and memory usage can only be achieved if ε ≤ 10 -3 . Though these solutions exhibit very low error magnitudes (Figure 5), the accuracy is significantly degraded by the perturbation error for ε = 10 -2 . Hence, ε = 10 -3 seems a good compromise that noticeably saves CPU-time and memory usage without spoiling the excellent accuracy.

Since the threshold parameter (ε) has a crucial influence on the error introduced by the multiresolution procedure [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF], it is important to check how MR solutions converge with ε. On figure 8, we plot the perturbation error in L 1 -norm, calculated as the difference between the FV solution and the MR solution for several ε values, obtained on 9 grid levels (i.e the finest grid is 512 × 512) at t = 10. The perturbation error recovers a linear fit versus ε which is in complete agreement with the heuristic MR approach (eq. 21) from Harten's work [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF]. This test-case is a good candidate to exhibit the complete power capacity of the procedure by enhancing accuracy performances and the perturbation error influence.

Interaction between a weak shock and a spot of temperature.

The goal is here to check the ability of the MR approach to predict a basic process of vorticity production induced by baroclinic effects [START_REF] Picone | Theory of vorticity generation by shock wave and flame interactions[END_REF][START_REF] Picone | Vorticity generation by shock propagation through bubbles in a gas[END_REF][START_REF] Pernaud-Thomas | Méthodes d'Ordre Elevé Appliquées au Calcul d'Ecoulements de Fluides Compressibles[END_REF]. This phenomenon can be observed during the interaction between a shock wave and a flame or a bubble. As it is well known in 2D compressible flows, without initial vorticity only the baroclinic effect can be responsible for the production of vorticity through the baroclinic torque, appearing in the vorticity transport equation:

-

1 ρ 2 ∇P × ∇ρ. ( 23
)
The flow configuration is the one studied by B. Pernaud-Thomas [START_REF] Pernaud-Thomas | Méthodes d'Ordre Elevé Appliquées au Calcul d'Ecoulements de Fluides Compressibles[END_REF] The spot of temperature is defined as an isobaric perturbation of temperature:

∆T (r) T 0 = 1 α 4 r 2 -α 2 2 e -r 2 /σ 2 , ( 24
)
with α = 7, σ = 0.07 and r = (x -x 0 ) 2 + (y -y 0 ) 2 . The center of the spot is located upstream of the shock wave at x 0 = 1/2, y 0 = 1/2. Initially, this perturbation of temperature is superimposed to the basic flow. On the upstream boundary (x = 0), since the inflow is supersonic, all the conservative quantities are prescribed. At the outlet (x = 2), non-reflecting conditions are applied, based on a characteristic approach. Periodicity is assumed in the vertical direction (y). We refer the reader to [START_REF] Tenaud | Evaluation of some high-order shock capturing schemes for direct numerical simulation of unsteady two-dimensional free flows[END_REF] for more details.

To ensure that the mesh has the same grid step in both directions in the rectangular domain ((x, y) ∈ [0, 2] × [0, 1]), two dyadic trees from two roots distributed in the x-direction (i.e. the first root at (x 0 1 , y 0 1 ) = (0.5, 0.5) and the second one at (x 0 2 , y 0 2 ) = (1.5, 0.5)) have been employed. If we use L grid levels in each part of the tree, the number of grid points on the finest grid will therefore be N x = 2 × 2 L , N y = 2 L . FV computations have been performed using 2 uniform grids: the coarsest grid 512 × 256 and the finest grid 1024 × 512. MR computations have also been performed on the same uniform grids, respectively the coarsest grid using 8 grid-levels and the finest grid using 9 Fig. 9. Interaction of a weak shock with a spot of temperature: density contours at t = 0.5 on the left, and at t = 1 on the right, computed with OSMP7 scheme and the MR approach using 9 grid levels and ε = 10 -3 (21 contours from 0.6 to 1.29). Fig. 10. Interaction of a weak shock with a spot of temperature: Pressure contours at t = 0.5 on the left, and at t = 1 on the right, computed with OSMP7 scheme and the MR approach using 9 grid levels and ε = 10 -3 (21 contours from 0.53 to 0.75).

grid-levels. Several threshold values ε= 10 -5 , 10 -4 , 10 -3 and 10 -2 , have been used for MR computations. The solution is also computed with ε = 0 using the highest grid-level everywhere in the domain with however the use of the MR algorithm. Simulations were conducted with the 2nd, 3rd and 7th order OSMP schemes, on a dimensionless time of unity using CF L = 0.5. Results are compared to the reference solution [START_REF] Tenaud | Evaluation of some high-order shock capturing schemes for direct numerical simulation of unsteady two-dimensional free flows[END_REF] obtained through a fully resolved simulation performed on a very fine grid (801x101) using the Hermit-6 scheme with about 10 points within the shock wave. To exhibit the nice results, we show (Figures 9, 10 and 11) the dimensionless density, pressure and vorticity contours at the intermediate dimensionless time (t = 0.5) and at the final time (t = 1), obtained with the OSMP7 scheme coupled with the MR technique using 9 grid levels and a threshold value ε = 10 -3 . These results are in complete agreement with the reference solution [START_REF] Tenaud | Evaluation of some high-order shock capturing schemes for direct numerical simulation of unsteady two-dimensional free flows[END_REF].

During the spot/shock interaction, the integral over the computational domain (Ω) of the vorticity (ω), the baroclinic torque [START_REF] Knight | Advances in CFD prediction of shock wave turbulent boundary layer interactions[END_REF], and their modulus are computed:

Ω ωdΩ, Ω |ω| dΩ, Ω ∇P × ∇ρ ρ 2 dΩ, Ω ∇P × ∇ρ ρ 2 dΩ. (25)
Fig. 11. Interaction of a weak shock with a spot of temperature: vorticity contours at t = 0.5 on the left, and at t = 1 on the right, computed with OSMP7 scheme and the MR approach using 9 grid levels and ε = 10 -3 (20 contours from -1 to +1).

Fig. 12. Interaction of a weak shock with a spot of temperature: adapted grid at t = 0.5 on the left, and at t = 1 on the right, obtained with OSMP7 scheme and the MR approach using 9 grid levels and ε = 10 -3 .

As far as the initial state does not exhibit initial vorticity and initial baroclinic torque, integrals of these quantities over the computation domain [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] must remain at zero. The maxima of these integral values, coming from global errors composed with approximation errors and perturbation errors, are reported in Table 3. The computer double-precision is easily recovered by the OSMP7 scheme when FV (ε = 0) is used on both grid levels. It is clear that the accuracy order of the scheme is of importance on the vorticity and the baroclinic productions since errors of lower order schemes (OSMP2, OSMP3) are far away from the computer accuracy even if the highest grid level is used. When a MR technique is employed with rather high threshold values (i.e., ε = 10 -2 and 10 -3 ), errors of the OSMP7 scheme are at least one order of magnitude lower than with OSMP2 and OSMP3 schemes. For lower ε values (ε = 10 -4 and 10 -5 ), OSMP7 errors drastically decrease. As expected, OSMP7 errors decrease more rapidly than for low order schemes as ε decreases.

Time histories of integrals of the vorticity and baroclinic torque magnitudes are plotted on Figure 13. The MR results obtained on 9 grid levels (i.e. using 1024 × 512 grid points on the highest level) with a threshold value ε = 10 -3 favorably compare to the reference solution (Hermit-6 scheme on a fine grid) [START_REF] Tenaud | Evaluation of some high-order shock capturing schemes for direct numerical simulation of unsteady two-dimensional free flows[END_REF]. The spot/shock interaction starts when the core of the vortex reaches the shock location, at about a dimensionless time t = 0.3. The maximum of the vorticity production is reached when the center of the spot coincides with 4.0219e-07 4.0007e-07 3.0446e-09 1.00652e-06 1.25717e-06 2.98157e-12 Fig. 13. Time history of the integral over the computational domain of |ω|, on the left, and the baroclinic torque, on the right, obtained by MR approach (9 grid levels with ε = 10 -3) with several orders of the OSMP scheme compared to the reference solution using a Hermit 6 scheme [START_REF] Tenaud | Evaluation of some high-order shock capturing schemes for direct numerical simulation of unsteady two-dimensional free flows[END_REF].

the shock location (Figure 13). As we can see on Figure 11, the vorticity is produced within and just downstream the shock as a pair of counter-rotating vortices. Low pressure structures associated with these vortices are visible in Figure 10. The interaction also produces a modification of the shock curvature which could be seen on density (Figure 9) and pressure (Figure 10) contours. When the spot is convected downstream of the weak shock, the intensity of the baroclinic torque decreases rapidly and the magnitude of vorticity consequently remains rather constant (Figure 13). At the final time, the shock wave that has been bent for the interaction still moves to recover its equilibrium position, and the spot located at x = 1.4 still interacts with the shock (Figure 10). This remaining interaction is responsible for the residual value of the vorticity production at the final time of the computation (Figure 13). The vorticity is concentrated in two counter-rotating vortices. A sign of a very good quality of the MR solution is that no residual vorticity is visible in regions where the curvature of the shock is still pronounced after the interaction (Figure 11).

The influence of the order of accuracy of the OSMP schemes is shown on Figure 13. The higher the order of accuracy, the smaller the production of vorticity. It is noticeable that OSMP7 gives very accurate results, since the maximum of baroclinic torque coincides with the hermit-6 scheme and since the intensity of vorticity after the interaction is lower than with the hermit-6 scheme. In fact, as far as the grid is much refined in the vertical direction than for the hermit-6 simulation, the vorticity relaxation and the remaining vorticity is much better captured by the OSMP7 scheme. We must also note that great discrepancies occur between the 2nd and 3rd order OSMP schemes and the OSMP7 scheme on both remaining vorticity after the interaction and baroclinic torque maximum production. As a conclusion on the influence of the order of accuracy, we could claim that the middling accuracy recovered disqualifies the OSMP2 and OSMP3 schemes. In the followings, only results obtained with the OSMP7 scheme will be discussed.

The corresponding adapted grids are also shown in Figure 12. Note that the grid adaptivity, based on both detail and predicting criteria (Eqs. 20 and 22), is capable of following travelling waves even though their magnitudes are weak. At t = 0.5, the shock wave, the spot of temperature and the weak circular pressure wave coming from their interaction are clearly enhanced with a locally high grid level. Since the initial solution does not satisfy the Navier-Stokes equations, plane acoustic waves are initially produced at the shock location and convected downstream. These acoustic waves are responsible for the high density of grid points visible close to the outlet of the domain. At t = 1, the grid refinement also educes acoustic wave interactions mainly allowed by the periodicity condition on the upper and lower boundaries. Otherwise, coarse grids are selected when the solution is smooth enough, especially in front of the shock wave where the solution recovers the uniform flow conditions after the interaction.

We report in Table 4 the CPU-time and memory usage of MR solutions using several ε values. The percentage of CPU time is relative to the FV values. As mentioned previously, the percentage of memory usage is the ratio between the number of leaves in the graded tree and the grid points of the equivalent FV mesh. Gains on memory usage are then reached as soon as leaves can be discarded from the tree-data structure. When ε = 0, the finest grid is used everywhere in the computational domain and we can then measure the specific cost of the MR algorithm. The MR algorithm cost is approximately 65 % greater than the one of the FV procedure. Gains on the CPU time can then be achieved, as far as 40 % of leaves are at least discarded from the tree-data structure. This is obtained when ε ≤ 10 -3 . To exhibit the convergence of the MR method, we plotted the perturbation error in L 1 -norm versus the threshold parameter (ε), for the simulation obtained on 9 grid levels (i.e the finest grid is 1024 × 512) at t = 1 (Figure 14); recalling that the perturbation error is calculated as the difference between the FV solution and the MR solution for several ε values. Following the heuristic MR approach (eq. 21) from Harten's work [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF], the perturbation error recovers a linear fit versus ε though a weak shock is present in the domain. Let say however that, using 9 grid levels, the weak shock is well resolved since more than 10 grid points are located in the thickness of the shock wave. In conclusion, as in the previous test-case, ε = 10 -3 seems the most convenient value since it saves both CPU-time and memory usage while giving an accuracy similar to the one of the FV solution on the regular mesh.

The viscous shock tube problem.

We consider the test-case studied in [START_REF] Daru | Evaluation of TVD high resolution schemes for unsteady viscous shocked flows[END_REF][START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF], i.e. a square shock tube with a unit side length and insulated walls. The diaphragm is initially located in the middle of the tube, i.e. x = 0.5. The initial states, in terms of dimensionless quantities, are:

on the left-side of the diaphragm: ρ L = 120,

P L = ρ L γ , u L = v L = 0,
and on the right-side:

ρ R = 1.2, P R = ρ R γ , u R = v R = 0.
The reference velocity is based on the initial speed of sound, corresponding to a reference Mach number M 0 = 1. The Reynolds number based on right-side initial quantities is prescribed at Re R = 200. At the initial time, the diaphragm is broken. A shock wave, followed by a contact discontinuity, moves toward the low-pressure region (i.e. to the right-side, in the present case), while a sonic rarefaction wave moves toward the high-pressure region on the left-side. The shock Mach number is equal to 2.37. Interactions between shock, rarefaction and expansion waves are well described in [START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF].

In the 2D viscous case, the propagating incident shock wave and contact discontinuity interact with the horizontal wall, creating a thin boundary layer.

After its reflection on the right end-wall, the shock wave interacts with this boundary layer. As the stagnation pressure in the boundary layer is lower than the one within the outflow region, a separation region, named bubble, appears over a large extent within the boundary layer, resulting in a major modification of the flow pattern and the formation of a lambda-shape like shock pattern, as one can see in Figure 15. The triple point emerging from the lambda-shape like shock pattern generates a slip line that rolls up in the right end corner.

To compare with the original solution obtained with the FV approach using the OSMP7 scheme (Figure 8-a of [START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF]), we also present a zoom in the contours of the density gradient in the vicinity of the interaction (Figure 16). Using ε = 10 -3 on L = 9 grid levels, the solution compares very well with the FV solution [START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF]. To emphasize these good results, a comparison on density distributions along the lower-wall of the shock-tube is presented in Figure 17.

When ε = 10 -2 , the slight perturbation errors introduce discrepancies when compared to the FV solution. More accurate results are obtained with ε = 10 -3 , which are in very good agreement with the reference solution.

Thanks to the MR approach, the adapted grid is tightened in regions where intense gradients occur. Note that the adapted mesh also improves the vertical waves emitted downstream the reflected shock-wave that however have a weak intensity (Figure 15). The gain in the memory usage achieved by the MR technique is calculated on the ratio between the number of active leaves used by the MR method and the number of grid points of the finest grid. The compression rate of grid points is very important for the first stages of the computation, since the solution is mostly composed with planar waves and a thin boundary layer. This compression rate rapidly evolves after the reflection of the shock wave occurs on the right end wall.

For ε = 10 -3 and 9 grid levels, the gain in the memory usage reaches its maximum value at the final time (t = 1). There, only 35 % of the finest grid points are used to encode the solution. Consequently, gains on the CPU time can be achieved. Only 80 % of the CPU time used by the FV computation on an equivalent single grid is here required to obtain the MR solution. This CPU-time seems to be lower than in the previous test-cases w.r.t. the memory saved. This may be explained by the few number of points used in first part of the computation until the shock wave reflects on the right end-wall. As before, the MR approach with ε = 10 -3 gives very accurate results at very low cost that is really competitive compared to the FV solution. Here we consider an inviscid spherical Riemann problem studied in [START_REF] Langseth | A wave propagation method for threedimensional hyperbolic conservation laws[END_REF], i.e. a spherical shock tube between two parallel solid plates situated at z = 0 and z = 1. Initially the gas is at rest everywhere. Inside a sphere with a radius r sp = 0.2, centered at x = 0, y = 0, and z = 0.4, density and pressure are initially prescribed at: ρ in = 1, P in = 5 while elsewhere density and pressure are: ρ out = 1, P out = 1. As in the original publication [START_REF] Langseth | A wave propagation method for threedimensional hyperbolic conservation laws[END_REF], by considering symmetry in (0, x, y) planes parallel to the end plates, the computational domain is limited to 1 4 of the whole domain: (x, y, z) ∈ [0, 1.5] × [0, 1.5] × [0, 1]; the z-axis being the symmetry axis (See Figure 18). The initial condition on the spherical Riemann problem is distributed on a Cartesian mesh using a distance function in order to avoid spurious oscillations. A FV solution is obtained on a fine grid with (N x × N y × N z ) = (150 × 150 × 100) equally spaced points. Two MR computations are also performed using respectively 5 and 6 grid levels on a tree forest corresponding respectively to a number of grid points at the highest level: (N x × N y × N z ) = (160 × 160 × 96) using 75 sub-trees (5 roots in x and y directions and 3 roots in the z direction), and (N x × N y × N z ) = (192 × 192 × 128) using 18 sub-trees (3 roots in x and y directions and 2 roots in the z direction). The initial solution and the corresponding adapted grid are shown on Figure 18. Using ε = 10 -3 , the adapted grid uses less than 1 % of the maximum grid points at the initial state.

Solutions are analyzed at the final dimensionless time t = 0.7, as in [START_REF] Langseth | A wave propagation method for threedimensional hyperbolic conservation laws[END_REF]. A shock wave followed by a contact discontinuity initially moves outwards to the low-pressure region and finally interacts with the solid walls. Expansion waves moving inwards to the high-pressure region converge to the center of the initial sphere, causing a local implosion. A secondary outwards moving shock wave is then created during the reflexion at the focal point. The main features of this test-case are the interaction between these two shock waves at different speeds and the reflection of these waves on solid end-walls, resulting in complex shock-shock interaction patterns.

Figure 19 shows the density contours at the dimensionless time t = 0.7. The adapted grid is also materialized by the center of leaves superimposed to the density contours. Setting ε = 10 -3 , the MR solution at t = 0.7 uses 50 % of grid points of the finest grid, distributed on 5 grid levels (from level 3 to level 7). To validate this MR solution, it is compared to the FV solution obtained with the same numerical scheme on a 150 x 150 x 100 grid. The pressure contours in a radial plane normal to the end-plates are presented in Figure 20.

Although very small discrepancies are visible, the MR solution compares very well with the FV solution we obtained using the OSMP7 scheme, but also with the solution found in [START_REF] Langseth | A wave propagation method for threedimensional hyperbolic conservation laws[END_REF]. Even though the adapted grid uses 50 % of the maximum grid points, the MR approach is here very competitive, since only 70 % of the FV CPU time is required to compute the final MR solution. This low CPU-time value w.r.t. the memory usage at the end of the computation is a consequence of the few number of points used in first part of the computation, i.e. until shock waves interact. Finally, to illustrate the good results obtained by the MR approach, the kinetic energy histories are compared between both MR on 2 adapted grids and FV solutions (Figure 21). The grid convergence seems to be obtained with the MR approach, leading to reliable solutions on both graded trees. Let us note that, compared to the FV approach, a very good agreement is achieved, meaning that the MR technique is very competitive compared to the classical FV method on a single grid.

Conclusion

In the present paper, we assess the capability of a new adaptive multiresolution technique coupled with an accurate one-step shock-capturing scheme, in order to speed up and improve the quality of transonic and supersonic flow computations with high frequency oscillations. We demonstrated its efficiency and showed its accuracy on several 2D and 3D prototype test-cases that are representative of elementary physical mechanisms in fluid mechanics. Both inviscid and viscous compressible flow phenomena are reviewed through the convection of an intense vortex and interactions between shock waves or shear layers. The starting point of the method is an accurate one-step (OS) finite volume scheme, developed following the Lax Wendroff approach. A monotony preserving (MP) criterion is added to preserve the accuracy near extrema. The dynamical adaptive strategy is based on the cell-average multiresolution analysis and a graded tree data structure that evolves with time. It is shown that the adaptive multiresolution enables to significantly speed up the computations with a very small loss in accuracy, thanks to the OSMP7 scheme, as long as the grid at the highest level is fine enough to capture the phenomenon. This loss of accuracy mainly dominated by the perturbation error is mainly related to the value of the threshold parameter (ε), which has a crucial influence on results. Hence, we checked how MR solutions converge with ε and what could be an adequate value of this threshold parameter. Through the prototype testcases, we showed that ε = 10 -3 is a good compromise that noticeably saves CPU-time and memory usage without spoiling the excellent accuracy provided by the OSMP7 scheme. With this threshold value, very good agreements with reference solutions obtained by means of the FV method on a single grid, are achieved by the MR technique, meaning that it is very competitive compared to single grid approaches.

Future work will deal with the simulation of turbulent flows using the MR approach. Current works are still in progress to compute accurately and with a significant CPU time reduction shock-turbulence interactions, mainly in the context of shock-turbulent boundary layer configurations, in transonic and supersonic flow regimes.

Fig. 1 .

 1 Fig. 1. Graded tree data structure in 1D: dyadic tree in red dash-dotted line, leaves in bold black line and virtual cells in blue dashed line.

Fig. 2 .

 2 Fig. 2. Ingoing and outgoing flux computations (arrows) in 2D on faces between adjacent cells at two different grid levels.

Fig. 3 .Fig. 4 .

 34 Fig. 3. 2D Euler test-case: advection of a strong vortex at a supersonic Mach number. On the left, density contours at the initial state and on the right, density contours obtained at t = 10 (after one turn) by MR (9 grid levels, ε = 10 -3 )

Fig. 8 .

 8 Fig. 8. L 1 -norm of the perturbation error versus the threshold parameter (ε) obtained for 9 grid levels (i.e. the finest grid is 512 × 512) at t = 10. The solid black line figures the linear regression.

Fig. 14 .

 14 Fig. 14. L 1 -norm of the perturbation error versus the threshold parameter (ε) obtained for 9 grid levels at t = 1.

Fig. 15 .

 15 Fig.[START_REF] Daru | High order one-step monotonicity preserving schemes for unsteady flow calculations[END_REF]. Viscous shock tube problem at t = 1 for Re = 200: solution obtained by the MR approach using 9 grid-levels, i.e. the finest grid has 1024 x 512 grid points, with ε = 10 -3 . On the left, iso-contours of the density (21 contours from 20 to 120). On the right, corresponding adapted grid using 5 encoding levels (from level 5 to level 9).

Fig. 16 .

 16 Fig. 16. Zoom in the interaction region for the viscous shock tube problem at t = 1 for Re = 200: On the top, adapted grid. On the bottom, contours of the density gradient.

Fig. 17 .

 17 Fig. 17. Distribution of the density on the lower-wall of the shock tube: MR solutions for two threshold values, i.e; ε = 10 -2 , red dashed line and ε = 10 -3 , blue dashed line, on 9 grid levels are compared with the FV solution (1001 × 501, black line).

Fig. 18 .

 18 Fig. 18. Initial state of the spherical Riemann problem: density contours and cell centers of the adapted grid, materialized with black dots (5 grid levels, ε = 10 -3 ). The lower and upper walls are enhanced in blue and red colors.

5. 4 A

 4 spherical Riemann problem.

Fig. 19 .

 19 Fig. 19. The spherical Riemann problem at the dimensionless time t = 0.7 obtained with the MR technique using 5 grid levels and ε = 10 -3 : density contours (on the right) and center of leaves of the adapted grid materialized with red dots (on the left). The lower and upper walls are enhanced in blue and red colors.

Fig. 20 .

 20 Fig. 20. Pressure contours, at the dimensionless time t = 0.7, in a radial plane of the spherical Riemann problem: FV solution obtained on a unique mesh using 150 × 150 × 100 grid (on the left) and the MR solution using 5 grid levels (the finest grid uses 160 × 160 × 96) and ε = 10 -3 (on the right).

Fig. 21 .

 21 Fig. 21. Time history of the kinetic energy over the domain: comparison between the FV solution obtained on a unique mesh using 150 × 150 × 100 grid and MR solutions using ε = 10 -3 on two adapted grids encoded on 5 grid levels (the finest grid uses 160 × 160 × 96) and 6 grid levels (the finest grid uses 192 × 192 × 128).

Table 2

 2 Advection of a strong vortex by a supersonic flow: gains and losses on both CPUtime and memory usage compared to the FV approach, for MR calculations using OSMP7 scheme and several ε values and several maximum grid levels.

	Method	N x × N y % CPU time % memory usage
	FV method	-	100 %	100 %
	MR -ε = 0	5 grid levels: 32 × 32	116 %	100 %
		6 grid levels: 64 × 64	137 %	100 %
		7 grid levels: 128 × 128	147 %	100 %
		8 grid levels: 256 × 256	148 %	100 %
		9 gris levels: 512 × 512	158 %	100 %
	MR -ε = 10 -2	5 grid levels: 32 × 32	109 %	83.59 %
		6 grid levels: 64 × 64	93 %	58.91 %
		7 grid levels: 128 × 128	75 %	47.52 %
		8 grid levels: 256 × 256	77 %	47.98 %
		9 grid levels: 512 × 512	62 %	35.97 %
	MR -ε = 10 -3	5 grid levels: 32 × 32	117 %	100 %
		6 grid levels: 64 × 64	125 %	83.37 %
		7 grid levels: 128 × 128	103 %	65.17 %
		8 grid levels: 256 × 256	92 %	60.60 %
		9 grid levels: 512 × 512	93 %	59.19 %
	MR -ε = 10 -4	5 grid levels: 32 × 32	117 %	100 %
		6 grid levels: 64 × 64	144 %	100 %
		7 grid levels: 128 × 128	130 %	91.10 %
		8 grid levels: 256 × 256	110 %	74.15 %
		9 grid levels: 512 × 512	106 %	70.93 %
	MR -ε = 10 -5	5 grid levels: 32 × 32	117 %	100 %
		6 grid levels: 64 × 64	143 %	100 %
		7 grid levels: 128 × 128	144 %	100 %
		8 grid levels: 256 × 256	146 %	95.03 %
		9 grid levels: 512 × 512	143 %	84.15 %

Table 3

 3 Maximum values of the vorticity and baroclinic torque productions using 2nd, 3rd and 7th order OSMP schemes for the MR approach using several ε values and both 8 and 9 grid levels.

	ε	level		max |	ω ∂Ω|		max |	∇P × ∇ρ ρ 2	∂Ω|
				Ω				Ω
			2 nd order	3 rd order	7 th order	2 nd order	3 rd order	7 th order
	0	8	1.1411e-09 9.2137e-10 3.5799e-16 5.02732e-09	2.7265e-09	1.1366e-15
		9	1.0644e-09 1.3290e-10 9.1135e-16 1.61482e-10 1.12989e-10 4.00384e-15
	10 -2	8	4.6253e-05 3.3874e-05 8.0392e-06 2.01097e-04	2.1252e-04	1.4123e-05
		9	4.3897e-05 3.9883e-05 3.2958e-06 3.09335e-04 2.66831e-04 3.30163e-05
	10 -3	8	7.0449e-06 6.6538e-06 4.3513e-07 5.35605e-05	5.8540e-05	1.1484e-06
		9	6.7396e-06 5.1923e-06 9.4462e-07 8.21467e-05 9.79718e-05 4.57065e-06
	10 -4	8	2.0455e-06 1.5565e-06 2.3525e-08 7.37053e-06	7.1631e-06	1.0664e-09
		9	3.

6096e-06 3.7981e-06 8.9881e-08 1.61768e-05 1.11815e-05 3.33593e-08 10 -5 8 2.2491e-07 2.0350e-07 3.9251e-09 5.31054e-09 2.8386e-09 4.1756e-12 9

Table 4

 4 Interaction of a spot of temperature and a weak shock for two meshes (512 × 256 and 1024 × 512) at the dimensionless time t = 1: CPU-time and memory usage for the MR approach with the OSMP7 scheme using several ε values. Gains and losses in CPU-time and memory usage refer to the FV approach.

	Method	N x × N y % CPU time % memory usage
	MR -ε = 0	8 grid levels: 512 × 256	165 %	100 %
		9 grid levels: 1024 × 512	164 %	100 %
	MR -ε = 10 -2	8 grid levels: 512 × 256	71 %	42.46 %
		9 grid levels: 1024 × 512	35 %	20.68 %
	MR -ε = 10 -3	8 grid levels: 512 × 256	97 %	66.40 %
		9 grid levels: 1024 × 512	74 %	51.10 %
	MR -ε = 10 -4	8 grid levels: 512 × 256	118 %	76.02 %
		9 grid levels: 1024 × 512	98 %	69.67 %
	MR -ε = 10 -5	8 grid levels: 512 × 256	118%	79.77 %
		9 grid levels: 1024 × 512	118 %	77.82 %

Acknowledgments

The authors would like to thank Marie Postel for her fruitful discussions and valuable comments. This work was partially supported by a scholarship from the French Ministère de l'Enseignement Supérieur et de la Recherche.