H Mathis

C Canc Ès

E Godlewski

N Seguin

DYNAMIC MODEL ADAPTATION FOR MULTISCALE SIMULATION OF HYPERBOLIC SYSTEMS WITH RELAXATION

Keywords: Hyperbolic system, finite volume method, relaxation, multiscale method, dynamic model adaptation, Chapman-Enskog expansion, two-phase flows. Mathematics Subject Classification. 35L45, 65M08, 65M55, 35C20, 76T10

In numerous industrial CFD applications, it is usual to use two (or more) different codes to solve a physical phenomenon: where the flow is a priori assumed to have a simple behavior, a code based on a coarse model is applied, while a code based on a fine model is used elsewhere. This leads to a complex coupling problem with fixed interfaces. The aim of the present work is to provide a numerical indicator to optimize to position of these coupling interfaces. In other words, thanks to this numerical indicator, one could verify if the use of the coarser model and of the resulting coupling does not introduce spurious effects. In order to validate this indicator, we use it in a dynamical multiscale method with moving coupling interfaces. The principle of this method is to use as much as possible a coarse model instead of the fine model in the computational domain, in order to obtain an accuracy which is comparable with the one provided by the fine model. We focus here on general hyperbolic systems with stiff relaxation source terms together with the corresponding hyperbolic equilibrium systems. Using a numerical Chapman-Enskog expansion and the distance to the equilibrium manifold, we construct the numerical indicator. Based on several works on the coupling of different hyperbolic models, an original numerical method of dynamic model adaptation is proposed. We prove that this multiscale method preserves invariant domains and that the entropy of the numerical solution decreases with respect to time. The reliability of the adaptation procedure is assessed on various 1D and 2D test cases coming from two-phase flow modeling.

Introduction

Interface coupling of existing numerical codes. The problem we address enters the framework of modeling of complex flows arising in industry. The phenomena we consider imply that we have to handle different scales both in time and space. It leads to the use of a hierarchy of models which differ according to these scales. As an example, one may think about water circuits in pressurized water reactors (PWR). It is clear that a complex model of phase transition has to be used in the steam generator or in the condenser while simpler models are sufficient to describe the flow in most parts of the pipes. Besides, in the case of loss-of-coolant accident, only much more accurate models should be able to describe these highly heterogeneous flows.

A direct consequence of the use of different models and numerical codes in disjoint parts of the computational domain is the development of theoretical and numerical techniques of multiscale coupling. The coupling has to be non intrusive because of the complexity of the codes under study, leading to methods which only make use of boundary conditions. This has been the subject of a series of works where several methods of coupling have been proposed for hyperbolic systems of partial differential equations [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF][START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. II. The case of systems[END_REF][START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Numerical coupling of two-phase flows[END_REF][START_REF] Ambroso | Coupling of general Lagrangian systems[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF][START_REF] Galié | Couplage interfacial de modles en dynamique des fluides[END_REF][START_REF] Boutin | Étude mathématique et numérique d'équations hyperboliques non-linéaires : couplage de modèles et chocs non classiques[END_REF]. In all these works, the interface of coupling which separates two different models is fixed. Since the models are different, there is not a consensual way of coupling, it depends on conservation and continuity principles which may be in contradiction (see for instance the examples in [START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF][START_REF] Ambroso | Coupling of multiphase flow models[END_REF]). As a consequence, the results vary according to the method of coupling and the more the flux is disturbed at the coupling interface, the greater is the sensibility to the coupling method. In order to avoid the dependence of the results on the interface coupling, it is natural to locate the interface where its impact is the most reduced.

This approach aims at reducing the complexity of the overall computation because solving the most accurate model is often out of reach or inappropriate. For example, in the reactor vessel, 3D models of multiphase flows with phase transition can be used, but, due to their complexity, such models might provide unexploitable results in most regions of the primary loop, where the classical one-phase 1D Euler equations with variable cross-section is more relevant. We are therefore concerned by using as much as possible the coarser model. Thus, in some cases, this method of model adaptation could lead to a decrease of the CPU cost, but it is not yet our purpose here.

Optimization of the position of the coupling interfaces. Our goal is to provide new tools to optimize the location of the interface, following the preliminary work [START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF]. The context is the following: assume that the overall flow can be fully described by an accurate model, the so-called fine model, which involves small scales. We aim at detecting the regions of the computational domain where this fine model can be replaced by a reduced one, the so-called coarse model, without deteriorating the accuracy of the results. This model can be obtained by setting the small scales to zero in the fine model. The whole procedure relies on the computation of a (numerical) indicator, depending on time and space, that allows to perform the parting of the computational domain into fine and coarse sub-domains. Such an indicator has to measure the difference between the solutions of the fine model and of the coarse one. Therefore, being given a threshold, it suffices to check whether the indicator is less than this threshold or not. In the first case, this means that replacing the fine model by the coarse model is possible, without introducing too much error. In the latter case, we continue to use the fine model. This leads to a multiscale method which dynamically selects the best model to use locally.

Of course, the development of such an indicator depends on the models under investigation. As mentioned above, we want to apply our method to compressible models of compressible two-phase flows. Due to the complexity of the underlying PDE's, see for instance [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF][START_REF] Drew | Theory of Multicomponent Fluids[END_REF], we are not able to derive robust and guaranteed criteria to position the interface, but the validation comes from the numerical experiments.

Let us begin by introducing the general theoretical framework and the so-called fine and coarse models.

Nonlinear systems of balance laws and hyperbolic relaxation. We focus in this paper on hyperbolic systems with relaxation. The fine model basically takes the form (more details are provided in Section 2.1)

(1)

∂ t W + d α=1 ∂ α F α (W ) = 1 ε R(W ),
with the initial condition W (0, x) = W 0 (x), which governs the evolution of the state vector W (x, t) : R d × R + → Ω ⊂ R n . Let us stress that ε is a positive constant (independent of space, time and W ) which is assumed to be small (but it does not tend to 0). Therefore, system (1) is said to be stiff. Following the pioneer works on hyperbolic systems with stiff relaxation [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Yong | Singular perturbations of first-order hyperbolic systems with stiff source terms[END_REF], we assume that there exists a linear operator P 1 : R n → R k of rank k < n such that P 1 R(W ) = 0 for all W ∈ Ω. In other words, if we denote u = P 1 W , system (1) includes k independent conservation laws

∂ t u + d α=1 ∂ α P 1 F α (W ) = 0.
We also assume that the equilibrium manifold {W ∈ Ω | R(W ) = 0} can be parameterized in terms of u, that is to say there exists M : R k → R n such that for all u we have R(M (u)) = 0. This leads to the definition of the equilibrium system

(2)

∂ t u + d α=1 ∂ α P 1 F α (M (u)) = 0,
which is the so-called coarse model. Additional assumptions ensure that the solutions of the fine model tend to solutions of the coarse model when ε → 0, see for instance [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Yong | Singular perturbations of first-order hyperbolic systems with stiff source terms[END_REF][START_REF] Tzavaras | Relative entropy in hyperbolic relaxation[END_REF], or in other words, that the leading order with respect to ε of the solutions of the fine model ( 1) is governed by the coarse model [START_REF] Allaire | A strictly hyperbolic equilibrium phase transition model[END_REF].

As mentioned above, we aim at defining an indicator which measures the difference between the solution W of the fine model [START_REF] Abdulle | The heterogeneous multiscale method[END_REF] and the solution u of the coarse model [START_REF] Allaire | A strictly hyperbolic equilibrium phase transition model[END_REF]. The main tool we use is based on the classical first-order Chapman-Enskog expansion

(3) W = M (P 1 W ) + εW 1 + O(ε 2 )
which provides a representation formula of the first order term W 1 . As a result, for some norm to be precised, εW 1 appears to be the right quantity for estimating the difference between the fine model and the coarse model. At the numerical level, the same method can be applied. On the basis of classical Finite Volume schemes for the fine model (1), we perform a numerical Chapman-Enskog expansion which provides the numerical counterpart of εW 1 . Note that, since the source term in (1) will be approximated in an implicit way, the computation of the numerical indicator εW 1 should not require to solve [START_REF] Abdulle | The heterogeneous multiscale method[END_REF]. In practice, this indicator measures the local smoothness of the solution and, in order to complete it, we also take into account the distance to the equilibrium manifold. Other indicators could be derived from error estimates of [START_REF] Tzavaras | Relative entropy in hyperbolic relaxation[END_REF] between general balance laws and their equilibrium counterparts; in [START_REF] Jin | A domain decomposition method for semilinear hyperbolic systems with two-scale relaxations[END_REF], a more academic case is studied, but accounting for the error due to the interface coupling.

Application to dynamic model adaptation. Let us now come back to our original problem of coupling two models through a thin interface. If we have such an indicator at our disposal, we may think that we are able to efficiently guess the position of the interface. In order to assess the relevance of our indicator, we insert it in a completely adaptive algorithm which automatically selects on-the-fly the "right" model to use. More precisely, being given a threshold θ > 0, we compare in each cell and at each time step our local indicator with θ: if it is greater than θ, we conclude that the difference between the two models is large, then the fine model is locally used; otherwise, we conclude that the coarse model can be used. Thus, at each time step, the computational domain is divided in several parts where the fine or the coarse model is used and at each interface, the two models are coupled. This method can be too complex to hope any CPU time saving in the present form, in particular when the dimensions of the fine and the coarse models are not very different. However, this numerical indicator and the adaptive algorithm allow to optimize the position of the coupling interfaces.

Comparisons with related methods. Many numerical methods exist to address model adaptation and multiscale methods. However, our motivations prevent us from applying standard techniques to our problem since we are in a situation where the numerical methods are already implemented in existing codes and cannot be modified.

In particular, asymptotic preserving schemes (let us only refer to the pioneer works [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Klar | An asymptotic preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF] and the reviews [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF]) are designed to provide a numerical scheme valid for both the fine model ( 1) and the coarse model [START_REF] Allaire | A strictly hyperbolic equilibrium phase transition model[END_REF] in such a way that consistency, stability (thus convergence) and accuracy of this numerical scheme are achieved independently of ε. Other similar techniques belong to the class of heterogeneous multiscale methods and micro-macro decompositions. See for instance [START_REF] Abdulle | The heterogeneous multiscale method[END_REF][START_REF] Kissling | The computation of nonclassical shock waves with a heterogeneous multiscale method[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Lemou | Micro-macro schemes for kinetic equations including boundary layers[END_REF] and also [START_REF] Fatone | Multimodels for incompressible flows: iterative solutions for the Navier-Stokes/Oseen coupling[END_REF][START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF][START_REF] Braack | A posteriori control of modeling errors and discretization errors[END_REF][START_REF] Perotto | Adaptive modeling for free-surface flows[END_REF][START_REF] Golse | A domain decomposition analysis for a two-scale linear transport problem[END_REF][START_REF] Degond | Macroscopic fluid models with localized kinetic upscaling effects[END_REF][START_REF] Bernardi | Automatic insertion of a turbulence model in the finite element discretization of the Navier-Stokes equations[END_REF][START_REF] Jin | A domain decomposition method for semilinear hyperbolic systems with two-scale relaxations[END_REF] for related works on coupling models with different scales. In particular, let us mention the work by Degond, Dimarco and Mieussens [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF] (and references therein, in particular [START_REF] Degond | A moving interface method for dynamic kineticfluid coupling[END_REF]) where similar ideas are used to design adaptive methods in a kinetic-fluid context.

Moreover, the authors with Frédéric Coquel propose in [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF] a complete analysis for a simpler model to understand how to balance the errors due to the coupling of models with the errors due to the dynamic and local replacement of the fine model by the coarse model. This analysis enables to design a multiscale method with model adaptation.

Outline of the paper. In Section 2 we present the hyperbolic model with relaxation we consider. Starting from the theory of Chen, Levermore and Liu [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF], we depict the general framework of hyperbolic systems with stiff relaxation. Performing a Chapman-Enskog expansion, we show that the smooth solutions of the relaxation system solve the associated equilibrium model up to second order terms. Thus we obtain an explicit formula of the first order term of the expansion, which is an element of the indicator we use thereafter to perform the adaptive method. The remaining part of this section is dedicated to the numerical schemes for the fine model and for the coarse model.

In Section 3 is presented the algorithm of model adaptation. The first step concerns the construction of the numerical indicator. It is obtained, on the one hand, by performing the same Chapman-Enskog expansion as in the previous section, but on the numerical scheme. On the other hand, we take into account the distance to the equilibrium manifold. As the coarse model is the hyperbolic limit of the fine relaxation model, we assume that the numerical scheme we use preserves this asymptotic property, at least for the algorithm we provide in this section (see Remark 1.1 for more details). This section ends with a thorough description of our multiscale method based on a dynamic selection of the model. We also prove that this method preserves admissible states and is entropy decreasing.

Section 4 is dedicated to numerical illustrations of model adaptation. We present two different applications, where we compare the numerical solutions obtained using the fine model, the coarse model and the result of model adaptation. The first example is a direct application of the previous sections. The systems are conservative and the entropy is strictly convex and satisfies the assumptions of [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]. The compatibility of the numerical scheme (asymptotic preserving property) allows us to perform exactly the adaptive method presented in Section 3. The second example models flow with phase transition. The entropy of the model is not strictly convex which leads to a degeneracy of the indicator for pure phases. In Section 5, we present the last example which concerns compressible two-phase flows described by the Baer-Nunziato model [START_REF] Baer | A two phase mixture theory for the deflagration to detonation (ddt) transition in reactive granular materials[END_REF]. This test case does not enter the frame of the previous sections since the fine model is not strictly hyperbolic and non conservative. It also involves a nonlinear source term. The coarse model is the so-called homogeneous two-fluid model, which is composed by four conservation laws (the Euler equations endowed with a transport equation on the fraction of volume of one phase). Here, we choose to drop the compatibility between the two numerical schemes. Thus we provide a numerical scheme for the coarse model that is not derived from the fine one. This implies that the numerical indicator cannot be deduced from a Chapman-Enskog expansion on the fine scheme, we thus modify our strategy by providing a classical discretization of the first order term obtained in the continuous Chapman-Enskog expansion.

Remark 1.1. For the sake of conciseness, the models and methods are presented in the standard relaxation framework defined in reference works, such as in [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]. But our approach easily extends to the applications we have in mind, which come from nuclear thermohydraulics, leading to quite complex systems (non conservative, non convex entropy. . . ). These models do not fulfill the assumptions we make in Section 2. However, we will see that an adaptation procedure can be implemented and provides interesting results.

Hyperbolic system with relaxation

2.1. Algebraic structure of relaxation. We consider the system of hyperbolic equations with relaxation terms (4)

∂ t W + d α=1 ∂ α F α (W ) = 1 ε R(W ),
with the initial condition W (0, x) = W 0 (x). The state vector W (t, x) : R + ×R d → Ω belongs to the convex set Ω ⊂ R n of admissible states. The flux and the source term F α , R : R n → R n , α = 1, . . . , d are supposed to be smooth. System (4) will be referred to in the sequel as the fine model. The chosen framework is that of Chen, Levermore and Liu [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] and we briefly list the different assumptions, referring to [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] for their justification. There exists a linear operator P 1 : R n → R k of rank k, 1 ≤ k < n, such that

P 1 R(W ) = 0, ∀W ∈ Ω.
The operator P 1 defines k conserved quantities u = P 1 W that satisfy

∂ t u + d α=1 ∂ α P 1 F α (W ) = 0.
There also exists another linear operator P 2 : R n → R n-k of rank nk such that the operator

P = P 1 P 2 is invertible in R n . Defining for all W ∈ Ω, α = 1, . . . , d, f α 1 (P 1 W, P 2 W ) = P 1 F α (W ), f α 2 (P 1 W, P 2 W ) = P 2 F α (W ), r(P 1 W, P 2 W ) = P 2 R(W ),
and setting v = P 2 W we can rewrite the system (4) as

∂ t u + d α=1 ∂ α f α 1 (u, v) = 0, (5)
∂ t v + d α=1 ∂ α f α 2 (u, v) = 1 ε r(u, v). (6)
The states W satisfying R(W ) = 0 are called equilibrium states or Maxwellians. The equilibrium manifold M = {W ∈ Ω | R(W ) = 0} is parameterized in terms of the k conserved quantities u = P 1 W and we can define a one-to-one mapping M : P 1 Ω → M such that u uniquely determines a local equilibrium value W eq = M (u) ∈ M and P 1 M (u) = u for all u ∈ R k . Let us also introduce the smooth map v eq : P 1 Ω → P 2 Ω defined by [START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF] v eq (u) = P 2 M (u)

and the operator P from R n to R n , satisfying

P : Z = (u, v) T → P M (u) = (u, v eq (u)) T . (8)
Then r(u, v eq (u)) = 0, for all u ∈ P 1 Ω. Furthermore we assume that the mapping

(9) v -→ r(u, v) is a C 1 -diffeomorphism on a neighborhood of v eq (u).
When ε is small with respect to the characteristic quantities of the system (4), the solutions of (4) are near equilibrium and we may consider the reduced system

∂ t P 1 W + d α=1 ∂ α P 1 F α (M (P 1 W )) = 0,
that can also be written [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF] ∂ t u

+ d α=1 ∂ α f α 1 (u, v eq (u)) = 0.
System [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF] is referred to in the sequel as the coarse model. Assume now that the relaxation system (4) is endowed with a convex entropy Φ : Ω → R satisfying:

(i) ∇ 2 Φ(W )∇F (W ) • ξ is symmetric for all W ∈ Ω and ξ ∈ R d , (ii) ∇Φ(W ) T R(W )≤0, ∀W ∈ Ω.
The condition (i) is the classical Lax condition for hyperbolic conservation laws (see e.g. [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]). It ensures the existence of an entropy flux Ψ :

Ω → R d such that (11) ∇Φ(W ) T ∇F (W ) = ∇Ψ(W ) T , ∀W ∈ Ω.
Every classical solution of (4) satisfies

(12) ∂ t Φ(W ) + d α=1 ∂ α Ψ α (W ) = 1 ε ∇Φ(W ) T R(W ),
and the second condition (ii) implies that the relaxation system is entropy dissipative since the right-hand side of ( 12) is non-positive. Besides the restriction of the entropy pair (Φ, Ψ) on the equilibrium manifold M

φ(u) := Φ(M (u)), ψ(u) := Ψ(M (u)),
gives an entropy pair (φ, ψ) for the system [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF], so that the equilibrium system is hyperbolic (see [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] for a detailed proof).

2.2. Chapman-Enskog expansion and dissipation. We now focus on the relaxation system [START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF]. Let us recall the classical result based on the Chapman-Enskog expansion.

Proposition 2.1. Up to terms of order ε 2 , the smooth solutions of the relaxation system (5-6) satisfy (at least formally)

∂ t u + d α=1 ∂ α f α 1 (u, v eq (u)) = -ε d α=1 ∂ α (∇ v f α 1 (u, v eq (u))v 1 ) , ( 13
) v = v eq (u) + εv 1 , (14)
where (15)

v 1 = ∇ v r(u, v eq (u)) -1 d α=1 ∂ α f α 2 (u, v eq (u)) -∇v eq (u) T d α=1 ∂ α f α 1 (u, v eq (u)) .
Proof. Let us consider the Chapman-Enskog expansion

v ε = v eq (u) + εv 1 + O(ε 2 ).
Plugging it into (5-6) leads to

∂ t u + d α=1 ∂ α f α 1 (u, v eq (u)) = -ε d α=1 ∂ α (∇ v f α 1 (u, v eq (u))v 1 ) + O(ε 2 ), ∂ t v eq (u) + d α=1 ∂ α f α 2 (u, v eq (u)) = -∇ v r(u, v eq (u))v 1 + O(ε).
The first system of equations is exactly (13) up to ε 2 terms. If we multiply it by ∇v eq (u) T and combine it with the second system of equations, we obtain the expression of the first order correction term v 1 [START_REF] Caetano | Sur certains problèmes de linéarisation et de couplage pour les systèmes hyperboliques non linéaires[END_REF].

Note that the first order term v 1 is merely an explicit function of u.

The system (13-15) is a closed system of order 2 which can be seen as an intermediate model between the fine model [START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF] and the equilibrium model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF]: formally, smooth solutions of (13-15) solve [START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF] up to terms of order ε 2 and, on the other hand, when ε tends to 0 one recovers the equilibrium model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF].

For stability reasons, one may expect that the second order term in (13-15) is dissipative. Note that this hypothesis is not verified without further assumptions, but following [START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservation laws[END_REF] and [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF], the existence of a uniformly convex entropy satisfying (i) and (ii) implies that the second order correction term of [START_REF] Boutin | Étude mathématique et numérique d'équations hyperboliques non-linéaires : couplage de modèles et chocs non classiques[END_REF][START_REF] Braack | A posteriori control of modeling errors and discretization errors[END_REF][START_REF] Caetano | Sur certains problèmes de linéarisation et de couplage pour les systèmes hyperboliques non linéaires[END_REF] is dissipative (this is the so-called sub-characteristic condition).

Remark 2.1. Let us recall that the Chapman-Enskog expansion is not valid near discontinuities, the above calculus requires the solution to be smooth.

Numerical schemes.

We now detail the scheme associated with the fine model [START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF].

Let us consider a mesh T h of the computational domain D ⊂ R d made of cells K, that are polygonal disjoint open subsets. If K and L are two neighboring cells (in the sequel we will use the notation L ∈ N (K)) their common face is denoted e KL = K ∩ L and n KL is the outgoing unit normal to K. The notations |K| and |e KL | refer to the volume of the cell K and the surface of the edge e KL respectively. We denote ∆t the time step and the sequence of time steps (t n ) n such that t n+1 = t n + ∆t, for all n ∈ N. We define the approximation of the initial data by

W 0 K = 1 |K| K W 0 (x)dx.
The finite volume formulation is obtained by integrating the system (5-6) on the space-time domain [t n , t n+1 ] × K and we introduce the two numerical fluxes F 1 and F 2 , respectively consistent with the fluxes f 1 and f 2 in (5-6). The relaxation system is approximated using a splitting strategy between the convective part and the source term. Suppose the approximation W n K , for all K ∈ T h , is known at time t n and let us introduce the notation

Z n K = (u n K , v n K ) T = P W n K .
For simplicity, the numerical schemes are expressed in terms of Z n K instead of W n K . In a first step, from t n to an intermediate time t n+1,-, the convective part is approximated by

u n+1,- K = u n K - ∆t |K| L∈N (K) |e KL |F 1 (Z n K , Z n L , n KL ), (16) v n+1,- K = v n K - ∆t |K| L∈N (K) |e KL |F 2 (Z n K , Z n L , n KL ). ( 17
)
Then the value Z n+1,- K is taken as the initial data for solving the source term:

u n+1 K = u n+1,- K , ( 18
)
v n+1 K = v n+1,- K + ∆t ε r(u n+1 K , v n+1 K ). ( 19
)
The classical implicit Euler scheme is chosen in order to ensure the unconditional stability of the second step.

Remark 2.2. Somehow formally, we assume that the drawback of this method is due to the resolution of the source term, so that the underlying resolution may be costly and complex. For instance, alternative methods can be found in [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], and in [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] for high order time discretizations.

One may easily deduce a numerical scheme for the coarse model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF]. The first step remains unchanged. The second step reduces to take v at equilibrium. The resulting numerical scheme is simply

u n+1 K = u n K - ∆t |K| L∈N (K) |e KL |F 1 (Z n K , Z n L , n KL ), (20)
v n+1 K = v eq (u n+1 K ). ( 21
)
The use of this numerical scheme avoids any computation of v: [START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF] is not used and the discretization of the stiff source term ( 19) is replaced by [START_REF] Coquel | Well-posedness and singular limit of a semilinear hyperbolic relaxation system with a two-scale discontinuous relaxation rate[END_REF].

In practice, the numerical scheme for the coarse model ( 10) can be different (see Section 4). Nonetheless, requiring the above compatibility between the two numerical schemes enables us to define the adaptive method in a simpler way.

Model adaptation

As mentioned previously the problematic we are interested in is to guess where the fine model can be replaced by the coarse model.

The decomposition of the computational domain into a fine and a coarse subdomain relies on the use of an indicator, the computation of which we now describe. This indicator relies both on a discretization of the first order correction term v 1 resulting from the Chapman-Enskog expansion (see Section 2.2, proposition 2.1), and on the distance to the equilibrium manifold, which is more precisely vv eq (u). This indicator allows us to realize the dynamic decomposition of the space and time computational domain.

At the end of this section, we precise the adaptive method based on the coupling of the relaxation model and of the equilibrium one at each sub-domain interfaces.

3.1. General ideas. At each time step, the strategy is to detect a sub-domain (non necessarily connected) where the fine model can be replaced by the coarse one. Let D f (t) (resp. D c (t)) be the sub-domain where the fine model (resp. the coarse model) has to be solved. We impose that the sub-domains do not intersect:

D f (t) ∪ D c (t) = D, D f (t) ∩ D c (t) = ∅.
At the interfaces between the sub-domains D f (t) ∩ D c (t) a coupling strategy has to be applied.

The first part of the indicator is based on the Chapman-Enskog expansion described in Proposition 2.1. We make the assumption that the difference between the solution of the fine model [START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF] and the solution of the coarse model ( 10) can be measured by some norm of εv 1 , defined by [START_REF] Braack | A posteriori control of modeling errors and discretization errors[END_REF][START_REF] Caetano | Sur certains problèmes de linéarisation et de couplage pour les systèmes hyperboliques non linéaires[END_REF]. Therefore, if εv 1 is less than a given threshold, then the coarse model is solved and if εv 1 is larger, the fine model is solved.

It is well-known that the Chapman-Enskog expansion is only valid for smooth solutions and in the neighborhood of the equilibrium manifold M. The quantity εv 1 measures the local smoothness of the solution. Since it vanishes for constant in space solutions, we have to consider also the distance to the equilibrium manifold M, by simply computing vv eq (u) , which in some sense measures the error with respect to time (see Appendix A for more details). The resulting indicator formally writes E = max(ε v 1 , vv eq (u) ), and depends on t and x.

A tolerance θ being given, the parting procedure reads as follows:

• The region where E ≤ θ is chosen to be D c (t). In that domain the error between the equilibrium model and the relaxation one is assumed to be negligible, so that the coarse model ( 10) is applied. • The domain D f (t) corresponds to the region where E > θ and the relaxation model (5-6) is solved inside. • At the interfaces separating the sub-domains D c (t) and D f (t), a numerical coupling method, like the ones developed in [START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF][START_REF] Caetano | Sur certains problèmes de linéarisation et de couplage pour les systèmes hyperboliques non linéaires[END_REF], is used.

Numerical indicator.

Let us now build the numerical indicator. To do so, we first use the Chapman-Enskog expansion at the discrete level on the numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF].

Proposition 3.1. Let us consider the numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF]. Then, up to terms of order ε 2 , one has

(22) v n+1 K = v eq (u n+1 K ) + εv n+1 1,K , where (23) v n+1 1,K = ∇ v r(u n+1 K , v eq (u n+1 K )) -1 |K| L∈N (K) |e KL | F 2 (P(Z n K ), P(Z n L ), n KL ) +∇v eq (u) T F 1 (P(Z n K ), P(Z n L ), n KL ) ,
where P is defined in (8) and u satisfies

∇v eq (u) T (u n+1 -u n ) = v eq (u n+1 ) -v eq (u n ).
The term v n+1 1,K is a first order approximation of v 1 at time t n+1 in cell K.

Proof. We plug the ansatz [START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF] in the numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF] and, noting

Z n 1,K = (0 R k , v n 1,K ) T for all K ∈ T h , we obtain (24)
u n+1 K = u n K + ∆t |K| L∈N (K) |e KL | F 1 (P(Z n K ), P(Z n L ), n KL ) + ε ∇ 1 F 1 (P(Z n K ), P(Z n L ), n KL ) T Z n 1,K + ε ∇ 2 F 1 (P(Z n K ), P(Z n L ), n KL ) T Z n 1,L + O(ε 2 ),
where

∇ β F 1 (Z 1 , Z 2 , n), β = 1, 2, is the partial derivative of F 1 with respect to Z β .
We also obtain (25)

v eq (u n+1 K ) = v eq (u n K ) + ∆t |K| L∈N (K) |e KL |F 2 (P(Z n K ), P(Z n L ), n KL ) + ∆t ∇ v r(u n+1 K , v eq (u n+1 K )) T v n+1 1,K + O(ε)
. Since v eq is smooth, there exists u(., .) such that:

∀u 1 , u 2 ∈ P 1 Ω, v eq (u 2 ) -v eq (u 1 ) = [∇v eq (u(u 1 , u 2 ))] T (u 2 -u 1 ).
Multiplying [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF] by ∇v eq (u(u n K , u n+1 K ))

T and dropping terms of order ε leads to

v eq (u n+1 K ) = v eq (u n K ) - ∆t |K| ∇v eq (u) T L∈N (K) |e KL |F 1 (P(Z n K ), P(Z n L ), n KL ).
Combining with [START_REF] Degond | Macroscopic fluid models with localized kinetic upscaling effects[END_REF] and dropping terms of order ε provides [START_REF] Degond | A moving interface method for dynamic kineticfluid coupling[END_REF]. Using the consistency of the numerical fluxes, it is easy to check that ( 23) is a consistent approximation of the continuous formula [START_REF] Caetano | Sur certains problèmes de linéarisation et de couplage pour les systèmes hyperboliques non linéaires[END_REF].

Note that replacing the terms v n 1,K and v n 1,L into (24) using the expression (23) could allow us to determine the discrete counterpart of the parabolic system [START_REF] Boutin | Étude mathématique et numérique d'équations hyperboliques non-linéaires : couplage de modèles et chocs non classiques[END_REF]. Moreover, one can check that the zeroth order part of ( 24) coincides with the coarse scheme [START_REF] Coquel | Relaxation of fluid systems[END_REF][START_REF] Coquel | Well-posedness and singular limit of a semilinear hyperbolic relaxation system with a two-scale discontinuous relaxation rate[END_REF], which is consistent with what holds at the continuous level.

In practice, the terms ∇v eq (u(u n K , u n+1 K )) and (∇ v r(u n+1 K , v eq (u n+1 K ))) -1 are respectively approximated by ∇v eq (u n K ) and (∇ v r(u n K , v eq (u n K ))) -1 to make the computation of the indicator fully explicit. This leads to the following modified expression of v n+1 1,K :

(26)

vn+1 1,K = (∇ v r(u n K , v eq (u n K ))) -1 |K| L∈N (K) |e KL | [F 2 (P(Z n K ), P(Z n L ), n KL )
+∇v eq (u n K ) T F 1 (P(Z n K ), P(Z n L ), n KL ) . Moreover, we include in the numerical indicator the distance to the equilibrium at time t n , which writes v n Kv eq (u n K ) . We deduce the following numerical indicator :

(27) E n+1 K := max(ε vn+1 1,K , v n K -v eq (u n K )
). This provides an estimate of the difference of the solutions given by the two schemes at time t n+1 , as an explicit function of the discrete solution (u n K , v n K ), for all K ∈ T h , (at time t n ) that is known within the whole computational domain. The first argument in the right-hand side [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] measures the local smoothness due to the relaxation source term. Even if εv n+1

1,K provides an estimate of the error at time t n+1 , it is not sufficient: in the case of constant in space solutions away from equilibrium, this quantity vanishes. However, in such a case, the second argument is positive and therefore the numerical indicator E n+1 K remains relevant (see Appendix A for more details).

3.3.

Model adaptation and coupling. The model adaptive method is performed at each time step t n . Assume that the subdivision of the computational domain into D f (t n ) × (t n , t n+1 ) and D c (t n ) × (t n , t n+1 ) is known. The adaptive method we propose consists in solving the following spatially coupled problem for all n ∈ N:

     (4) is solved in D f (t n ) × (t n , t n+1 ), (10) is solved in D c (t n ) × (t n , t n+1 ), + coupling conditions in D f (t n ) ∩ D c (t n ) × (t n , t n+1 ).
Several coupling conditions can be considered. In the following we favor state coupling, as presented for instance in [START_REF] Caetano | Sur certains problèmes de linéarisation et de couplage pour les systèmes hyperboliques non linéaires[END_REF][START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF]. Let K be a cell belonging to the fine domain D f (t n ) and L its neighbor belonging to the coarse domain D c (t n ). Then at the edge e KL , we impose (u n L , v eq (u n L )) as Dirichlet boundary condition for the fine model and u n K as Dirichlet boundary condition for the coarse model: if

K ∈ D f (t n ) and L ∈ D c (t n ), the numerical fluxes are (C1) F 1 (Z n K , P(Z n L , )n KL ), F 2 (Z n K , P(Z n L , )n KL ) for the computation of Z n+1,- K by the fine scheme (16-17), (C2) F 1 (P(Z n L , )Z n K , n LK ) for the computation of Z n+1
L by the coarse scheme (20), written on cell L. This amounts to impose in a weak sense the continuity of u and v = v eq (u) through the coupling interface (see [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF][START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. II. The case of systems[END_REF]).

Remark 3.1. One may impose other coupling conditions. For instance the equation on u being solved everywhere, one may impose the continuity of the flux f 1 through the interface e KL . This strategy of flux coupling is not tested here but may have special interest in some cases (see for instance Section 4). The algorithm will put the interfaces so that, in their neighborhood, the fine and the coarse models are mostly equivalent, the fine model being close to equilibrium, thus all the coupling strategies may be supposed to provide similar results.

3.4. Adaptive method. We now detail the general algorithm of the dynamic coupling between the fine and the coarse models. Let Z n K be the solution in the cell K known at time t n to be updated to time t n+1 . The algorithm reads as follows (recall that θ is a given threshold for defining the partition):

A) For all cell K ∈ T h , compute the numerical indicator E n+1 K given by (27) B) For all cell K ∈ T h , if [|E n+1 K | > θ] then -K ∈ D f (t n ) Else -K ∈ D c (t n ) C) At this stage, D f (t n ) ∪ D c (t n ) = D. For all cell K ∈ T h : -If [K ∈ D f (t n ) and for all L ∈ N (K), L ∈ D f (t n )] -Compute Z n+1 K using the numerical scheme (16-19). -If [K ∈ D c (t n ) and for all L ∈ N (K), L ∈ D c (t n )]
-Compute Z n+1 K using the numerical scheme [START_REF] Coquel | Relaxation of fluid systems[END_REF][START_REF] Coquel | Well-posedness and singular limit of a semilinear hyperbolic relaxation system with a two-scale discontinuous relaxation rate[END_REF]. -Else -Compute Z n+1 K using the state coupling method (C1-C2).

This adaptive method satisfies some basic properties:

Proposition 3.2. (i)
The adaptive method (A-B-C) is conservative with respect to u. (ii) If the initial condition W 0 is a constant equilibrium (i.e. W 0 ∈ M), then for all K ∈ T h and n ≥ 0, Z n K = P W 0 . In other words, the adaptive method exactly provides the constant solution of the original relaxation system (1).

These properties (also discussed in [START_REF] Degond | A moving interface method for dynamic kineticfluid coupling[END_REF] for instance) can be easily proved. The first one directly follows from the conservation assumption

F 1 (Z n K , P(Z n L ), n KL ) = -F 1 (P(Z n L ), Z n K , n LK )
which is verified by the coupling conditions (C1-C2). The second one is deduced by the consistency of the numerical fluxes and by the use of the state coupling (C1-C2).

In order to investigate the nonlinear stability of the adaptive method, let us focus on an equivalent formulation.

3.5. Alternative formulation and properties. Using notations f = (f 1 , f 2 ) T and S = (0, r) T , system (5-6) can be written as ( 28)

∂ t Z + d α=1 ∂ α f α (Z) = 1 ε S(Z).
Following [START_REF] Gosse | Time-splitting schemes and measure source terms for a quasilinear relaxing system[END_REF] and denoting by Z(t n+ ) the right-value of t → Z(t) at t = t n , the classical time-splitting discretization of ( 28) can be written (29)

∂ t Z + d α=1 ∂ α f α (Z) = ∆t ε ∞ n=1 S(Z(t n+ ))δ(t -t n ),
where δ stands for the classical Dirac measure. System ( 29) is the continuous-inspace version of the fine numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF]. Indeed, writing the Rankine-Hugoniot jump relations for (29) through {t = t n } yields

Z(t n+ ) -Z(t n-) = ∆t ε S(Z(t n+ ))
where Z(t n-) is the left-value of t → Z(t) at t = t n , which is exactly the implicit Euler method for ∂ t Z = 1 ε S(Z). On the other hand, since the coarse numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF] and (20-21) fits the relaxation process discretized in [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF], the coarse scheme can be interpreted as [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] ∂

t Z + d α=1 ∂ α f α (Z) = ∞ n=0 P(Z(t n-)) -Z(t n-) δ(t -t n ),
where P is defined in [START_REF] Ambroso | Coupling of general Lagrangian systems[END_REF]. Indeed, [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] yields

Z(t n+ ) -Z(t n-) = P(Z(t n-)) -Z(t n-) =⇒ Z(t n+ ) = P(Z(t n-)).
Now, let us recall that the decomposition of the computational domain is performed according to the sign of [E(Z(t n-))θ]. Hence the semi-discrete method of model selection can be formally written as (31)

∂ t Z + d α=1 ∂ α f α (Z) = ∆t ε ∞ n=1 S(Z(t n+ ))δ(t -t n )1 {E(Z(t n-))>θ} (x) + ∞ n=0 P(Z(t n-)) -Z(t n-) δ(t -t n )1 {E(Z(t n-))≤θ} (x),
the system (31) being complemented with the initial data Z(t 0-) = Z 0 . At this stage, in order to describe thoroughly the adaptive method, it remains to:

(i) define the coupling at each interface between the two sub-domains, (ii) introduce the space discretization (for instance, in the context of the Godunov scheme, it suffices to use the L 2 -projection on the set of constant-by-cell functions).

As far as the definition of the coupling is concerned, the use of the classical state coupling deeply simplifies the adaptive method since, as mentioned before, the same numerical flux F 1 is used in the coupling conditions (C1-C2) for both cells at the coupling interface, with the result that the whole computation of u is conservative.

Assuming that Z(t n , •), D f (t n ) and D c (t n ) are defined (which is the case after steps A and B), step C can be split in the following way:

C'.a) Project onto the equilibrium manifold in the coarse sub-domain:

(32) Z a (t n , x) = P(Z(t n , x)) if x ∈ D c (t n ), Z(t n , x) if x ∈ D f (t n ).
C'.b) Solve the Cauchy problem of the homogeneous part of ( 31):

(33)

   ∂ t Z b + d α=1 ∂ α f (Z b ) = 0, t ∈ (t n , t n+1 ), x ∈ R d , Z b (t n , x) = Z a (t n , x), x ∈ R d ,
and obtain Z b (t n+1 , •). C'.c) Solve the source terms of (31):

• In D f (t n ), solve the Cauchy problem (34)

∂ t Z c = ∆t ε S(Z c ), t ∈ (t n , t n+1 ), x ∈ D f (t n ), Z c (t n , x) = Z b (t n+1 , x), x ∈ D f (t n ). • In D c (t n ), compute (35) Z c (t n+1 , x) = P Z b (t n+1 , x) , x ∈ D c (t n ).
C'.d) Update the value of Z:

(36) Z(t n+1 , x) = Z c (t n+1 , x), x ∈ R d .
If the homogeneous system ( 33) is solved by the fine numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF]) and if the system of ODE's ( 34) is solved by the implicit Euler scheme [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF], then the steps (C'.a-C'.d) are equivalent to step C. This decomposition of the adaptive method enables us to state the following fundamental stability property:

Proposition 3.3. The continuous-in-space version of adaptive method (A-B-C') is entropy decreasing, which means that its solution W = P Z satisfies

(37) d R Φ(W (t n+1 , x))dx ≤ d R Φ(W (t n , x))dx.
Moreover, if the fine numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF] is entropy satisfying in the sense of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], then the fully discrete adaptive method (A-B-C) is entropy decreasing, that is to say

(38) K∈T h Φ(W n+1 K ) ≤ K∈T h Φ(W n K )
where

W n K = P Z n K .
The proof is straightforward since each step C'.a and C'.c makes the entropy decrease pointwise, that is to say for almost every x ∈, Φ(P Z a (t n , x)) ≤ Φ(P Z(t n , x)) and Φ(P Z c (t n+1 , x)) ≤ Φ(P Z c (t n , x)).

A similar proposition holds if one assumes that Ω is a convex invariant domain, that is to say: if W 0 (x) ∈ Ω for all x ∈ R d , then the solution W of (4) satisfies W (t, x) ∈ Ω for all (t, x) ∈ R + × R d . Proposition 3.4. Assume that Ω is a convex invariant domain. The continuousin-space version of adaptive method (A-B-C') preserves the invariant domain Ω. Moreover, if the fine numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF] preserves the invariant domain Ω, then the fully discrete adaptive method (A-B-C) preserves the invariant domain Ω.

The proof is also straightforward since one may easily check that each step of the adaptive methods preserves the invariant domain.

Numerical applications

We provide in this section some numerical examples to illustrate the reliability of the indicator [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] and the relevance of its use for dynamic model adaptation.

In particular, we do not restrict ourselves to models which strictly fulfill the assumptions we made in Section 2, extending our method to more complex systems, of great interest for the applications we have in mind (which come from thermohydraulics in the context of nuclear reactors). In what follows, only the Chaplygin gas model presented in Section 4.1 falls in this theoretical framework.

First we consider the Chaplygin gas model, that consists in the barotropic Euler equations with relaxation, corresponding to rate-type fluids of Suliciu [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions[END_REF]. The simple form of the source term allows us to exhibit a dissipative convex entropy. We present one test case, inspired by [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF] (we also provide another test case for this system in Appendix A which concerns the study of the different parts of the numerical indicator).

Secondly, we consider a phase transition model. In that case the entropy is no longer strictly convex and the Chapman-Enskog expansion may not be valid. In particular, the Chapman-Enskog part of the indicator [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] vanishes in pure phases.

We will address in the last section the approximation of the solutions of the compressible seven-equation two-phase flow model. Although this model does not enter the frame of the section 2 (non convex entropy, non-strictly hyperbolic and non conservative model, incompatible numerical schemes for the fine and coarse models. . . ), the adaptation process is carried out, leading to a slightly different adaptive method.

In all the following tests, we compare the results provided by the numerical scheme [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF] for the fine model [START_REF] Ambroso | Coupling of multiphase flow models[END_REF][START_REF] Ambroso | Homogeneous models with phase transition: coupling by finite volume methods[END_REF], which corresponds to the reference solution, with the results provided by the adaptive method developed in Section 3. We also show the results obtained by the numerical scheme [START_REF] Coquel | Relaxation of fluid systems[END_REF][START_REF] Coquel | Well-posedness and singular limit of a semilinear hyperbolic relaxation system with a two-scale discontinuous relaxation rate[END_REF] for the coarse model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF] in order to verify that they are very different from those of the reference solution and that model adaptation is really needed. 4.1. Chaplygin gas model. We first consider a relaxation model which describes Chaplygin gas [START_REF] Serre | Multidimensional shock interaction for a Chaplygin gas[END_REF]. The one-dimensional problem we address corresponds to a fluid flow governed by the relaxation system (39)

∂ t τ -∂ x u = 0, ∂ t u + ∂ x Π(τ, T ) = 0, ∂ t T = 1 ε (T eq (τ ) -T ),
where T eq (τ ) = τ . This model is also derived from the works of Suliciu [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions[END_REF] on rate-type fluids. The state variable τ and u stand for the specific volume and the velocity while T is a perturbed specific volume. The extended pressure law Π is defined by

Π(τ, T ) = p(T ) + a 2 (T -τ ),
where p follows a perfect gas law p(T ) = T -γ , with γ = 1.4 and a = 1.5 in the tests hereafter. The associated equilibrium system is obtained setting T = T eq (τ ) and corresponds to the barotropic Euler equations in Lagrangian coordinates: [START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF] ∂ t τ -∂ x u = 0,

∂ t u + ∂ x p(τ ) = 0.
An admissible entropy for the Suliciu system is (see [START_REF] Coquel | Relaxation of fluid systems[END_REF])

Φ(τ, u, T ) = 1 2 |u| 2 + 1 1 -γ T 1-γ + a 2 2 (T 2 -τ 2 ) + (T -γ + a 2 T )(τ -T ).
This entropy is strictly convex and dissipative with respect to the source term under the so-called Whitham's condition

a 2 > max s (-p ′ (s)).
This condition also implies that the Chapman-Enskog expansion is dissipative, see for instance [START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservation laws[END_REF]. We consider Rusanov numerical fluxes [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves with barriers[END_REF] and, using the previous computations [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF], we obtain

vn+1 1,K = - 1 2∆x (p(τ n K+1 ) + p(τ n K ) -a K+1/2 (u n K+1 -u n K ) -p(τ n K ) -p(τ n K-1 ) + a K-1/2 (u n K -u n K-1 )
) and the numerical indicator is defined by [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF]. (1, 0, 1) if x < 0, (0.8, 0, 0.8) else.

For the numerical simulations, the mesh is composed by 200 cells and the threshold θ for the selection of the appropriate model is equal to 5.10 -5 as in the previous case. The numerical results are gathered in Figures 1 and2. The solution is composed by two fronts which start from the center of the domain, one going to the left and the other to the right. The solution of the coarse model ( 40) is composed by a shock wave which moves to the left and by a rarefaction wave which moves to the right.

Numerical results at T 1 = 0.2. The reference solution provided by the fine model corresponds to the red curve (which is actually hidden by the blue points). One can see the formation of the two fronts, with an intermediate state. The difference between the profiles of τ and T shows that the fronts are far from the equilibrium T = T eq (τ ). In particular, the coarse model provides sharper fronts.

The results of the adaptive method are superposed with the results of the fine model. In particular, the diffused fronts are very well approximated for all the variables. The numerical indicator is plotted in Figure 2-right and we can see that the variations are localized at the fronts. It is very interesting to note that the numerical indicator detects the equilibrium region in the neighborhood of x = 0 and that the coarse model is used there while the fine model is used at the fronts.

Numerical results at T 2 = 0.4. The two fronts continue to propagate and the difference between the results for the fine model and for the coarse model are similar to those at the previous time and are located at the extremities of the front. The intermediate part is expending and is composed to a constant equilibrium state.

The regions where the numerical indicator varies still correspond to the fronts and the profiles provided by the adaptive method are very similar to the reference profiles. The fine model is only used by the adaptive method near the fronts and one can check that the the intermediate state is solved using the coarse model, see Figure 2 Numerical results at T 3 = 0.6. At this time, the two fronts begin to exit from the domain. Let us recall that we impose free boundary conditions. They do not affect the solutions, no reflection is detected.

Since the numerical indicator is local, any boundary condition which is imposed by a "ghost cell" method can be handled without any difficulty. Indeed, the results provided by the adaptive method are very close to those of the fine model and one can check that no spurious effect appears at the boundaries.

Numerical results at T 4 = 0.8. Here, the two fronts are almost out of the domain. The solution is composed by the end parts of the fronts and a large constant equilibrium state.

The partition obtained by the adaptive method is once again in agreement with the solution. The fine model is only used at the boundaries, where the fronts are still present. The domain where the coarse model is used now represents 70% of the whole domain and corresponds to the intermediate constant state. 4.1.2. Shock tube test case: convergence with respect to the threshold. We now illustrate the convergence of the solution of the adaptation algorithm towards the solution provided by the fine model. To do so, we compute the L 1 difference between the perturbed specific volume T obtained by the adaptation algorithm and the one obtained by the fine model, for several values of the threshold θ.

In Fig. 3, we plot in a Log-Log scale this difference with respect to θ. We can observe a monotone decrease of the error, of order 1. Moreover, the error is bounded by the difference between the fine and the coarse model (blue curve in Fig. 3), as expected.

4.2.

Phase transition model. We now address the numerical approximation of a liquid-vapor compressible flow (see [START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF]). Each phase, denoted β = 1 or 2, is described by its own Equation of State (EoS) and thermodynamical quantities: ρ β denotes the density, τ β = 1/ρ β the specific volume, e β the internal energy. We assume that both phases are described by a perfect gas law, such that the pressure p β , the temperature T β and the entropy of each phase are classically given by:

p β (ρ β , e β ) = (γ β -1)ρ β e β , T β = e β , s β = ln e β ρ γ β -1 β
, where γ β > 1 stands for the polytropic exponent. We also introduce the chemical potential

g β = p β ρ β -T β s β + e β .
The density of the fluid is related to those of both phases by the following relations

ρ = αρ 1 + (1 -α)ρ 2 , (ρ -1 =) τ = ϕτ 1 + (1 -ϕ)τ 2
, where α and ϕ denote respectively the volume fraction and the mass fraction of the gaseous phase. The thermodynamical equilibrium between both phases corresponds to the equality of pressures, temperatures and chemical potentials. It is characterized by the following fractions at equilibrium:

ϕ eq (ρ) = τ -τ * 2 τ * 1 -τ * 2 , α eq (ρ) = ρ -ρ * 2 ρ * 1 -ρ * 2 , τ * 1 = 1 ρ * 1 , τ * 2 = 1 ρ * 2 ,
together with

ρ * 1 = exp(-1) γ 2 -1 γ 1 -1 γ 2 γ 2 -γ 1 , ρ * 2 = exp(-1) γ 2 -1 γ 1 -1 γ 1 γ 2 -γ 1 .
∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(ρu ⊗ u) + ∇p = 0, ∂ t (ρE) + div((ρE + p)u) = 0, ∂ t (ρϕ) + div(ρuϕ) = 1 ε (ϕ eq (ρ) -ϕ), (41)
where E = e + |u| 2 /2 stands for the total energy of the fluid and the complete equilibrium mass fraction is defined by

ϕ eq (ρ) =      1 if ρ ≤ ρ * 1 , ϕ e (ρ) if ρ * 1 ≤ ρ ≤ ρ * 2 , 0 if ρ * 2 ≤ ρ.
To close the system, we consider the EoS

p = p(ρ, e, ϕ) = (γ(ϕ) -1)ρe, where γ(ϕ) = γ 1 ϕ + γ 2 (1 -ϕ).
Remark 4.1. The entropy of the system is not strictly convex, as shown in [START_REF] Jaouen | Étude mathématique et numérique de stabilité pour des modèles hydrodynamiques avec transition de phase[END_REF][START_REF] Helluy | Relaxation models of phase transition flows[END_REF][START_REF] Allaire | A strictly hyperbolic equilibrium phase transition model[END_REF][START_REF] Helluy | Pressure laws and fast Legendre transform[END_REF]. Moreover, the source term in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]

vanishes when ρ ∈ [ρ * 1 , ρ * 2 ]
. As a consequence, the map (9) is not invertible and the Chapman-Enskog expansion does not provide a dissipative system in the whole domain. Nonetheless, we still use the same numerical indicator in order to see the influence of this degeneracy.

Let us now present the equilibrium model. If α = 0 (resp. α = 1) then ϕ = 0 (resp. ϕ = 1) and we directly recover the Euler equations for the pure phase. If 0 < α < 1, the limit ε → 0 leads to g 1 = g 2 : the thermodynamical equilibrium is reached. This asymptotic defines the coarse model: (42)

∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(ρu ⊗ u) + ∇p = 0, ∂ t (ρE) + div((ρE + p)u) = 0, Fine model (4)
Coarse model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF] Adaptation with θ = 10 -2 Adaptation with θ = 1 where p = p(ρ, e, ϕ eq (ρ)) that is

(43) p = p(ρ, e) =      (γ 1 -1)ρe if ρ ≤ ρ * 1 , (γ 1 -1)ρ * 1 e if ρ * 1 ≤ ρ ≤ ρ * 2 , (γ 2 -1)ρe if ρ * 2 ≤ ρ.
Remark 4.2. The fine model (41) can be expected to be well-posed, when the two phases co-exist, following [START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF][START_REF] Yong | Entropy and global existence for hyperbolic balance laws[END_REF]. On the other hand, the result of ill-posedness by Chiodaroli, De Lellis and Kreml [START_REF] Chiodaroli | Global ill-posedness of the isentropic system of gas dynamics[END_REF] could be applied to the coarse model [START_REF] Helluy | Pressure laws and fast Legendre transform[END_REF], in space dimension greater that 1. These results are not in contradiction since the solutions of relaxation models are smooth while the ill-posedness result for the system of barotropic gas dynamics is related to discontinuous solutions.

The numerical schemes we use are based on the HLLC scheme [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] and the numerical indicator is given by [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF]. Note that the indicator vn+1 1,K is not a continuous function since ϕ ′ eq admits discontinuities at ρ * 1 and ρ * 2 .

Shock bubble interactions test case. The adaptive method is tested on a 2D weak interaction of a bubble with a planar shock wave. The domain is [-0.5, 0.5] 2 and the pressure laws are perfect gas equations of state with γ 1 = 1.

(ρ, u, v, p, ϕ)(t, x = -0.5, y) = (ρ * 2 , 0.1, 0, 1, 0). The initial data is (ρ, u, v, p, ϕ)(0, x, y) = (ρ * 2 , 0, 0, 1, 0) if x 2 + y 2 < 10 -2 , (ρ * 1 , 0, 0, 1 , 1) else,
which corresponds to a bubble of vapor surrounded by the liquid phase. The relaxation time ε is equal to 1. The solution is composed by a shock wave which impacts the bubble, rebounds on the right boundary and impacts once again the bubble. The solution is plotted after the first interaction at t 1 = 0.6 and during the second interaction at t 2 = 1.26. The mesh is composed of 5906 triangular cells and 3054 vertices. The mass fraction ϕ is depicted in Figures 4 and6 while the density ρ is depicted in Figures 5 and7. In all these figures, we compare the approximate solutions associated with

• the fine model ( 41),

• the coarse model ( 42),

• the adaptive method with a threshold θ = 1 and • the adaptive method with a threshold θ = 10 -2 .

Fine model (4)

Coarse model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF] Adaptation with θ = 10 -2 Adaptation with θ = 1 The coarse model provides solutions very different from those of the fine model; in particular, the bubble fully liquefies during the second shock-bubble interaction, see the mass fraction in Fig. 6 (top right).

The results with θ = 10 -2 are very accurate compared to the solution provided by the fine model. Only a qualitative difference can be seen in Fig. 6 (bottom left), where the interface between the two phases provided by the adaptive method is slightly sharper. On the contrary, the adaptive method with θ = 1 gives results much closer to those of the coarse model. In particular, in Fig. 6 (bottom right), the bubble has almost disappeared. We present in Fig. 8 the characteristic function of the coarse domain. For θ = 10 -2 , the fine domain is restricted to the bubble and one can see that for θ = 1, it only remains four cells in the fine domain at time t 2 , located at the center of the (former) bubble. This clearly explains the discrepancies between the results of the adaptive method with θ = 1 and those of the fine model.

Shock bubble interactions test case: convergence with respect to the threshold. We now present the convergence of the solution of the adaptation algorithm towards the solution provided by the fine model. As in the previous test case, we compute the difference between the mass fraction ϕ given by the adaptation algorithm and the one obtained by the fine model.

Once again, one can see in Fig. 9 that this difference tends to 0 when θ tends to 0 (actually, for θ < 10 -10 , the difference is exactly 0, which explains the vertical Fine model [START_REF] Ambroso | Numerical coupling of two-phase flows[END_REF] Coarse model [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF] Adaptation with θ = 10 -2 Adaptation with θ = 1

Extension to a compressible two-phase flow model

We now consider the computation of a one-dimensional compressible two-fluid flow: each phase is considered as a single phase separated from the other. The two phases have their own thermodynamics and distinct velocities. Thus the balance equations can be given for both phases adding exchange terms between the two phases through the interfaces. This model was first proposed by Baer and Nunziato [START_REF] Baer | A two phase mixture theory for the deflagration to detonation (ddt) transition in reactive granular materials[END_REF] and then widely studied, see for instance, in a non exhaustive way, Abgrall and Saurel [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase[END_REF], Gallouët et al [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF], Ambroso et al [START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF], Saleh [START_REF] Saleh | Analyse et simulation numérique par relaxation d'écoulements diphasiques compressibles. Contribution au traitement des phases évanescentes[END_REF].

Considering our adaptation procedure for such a model is a very challenging benchmark test because this two-fluid flow model does not enter the frame of the previous sections. Indeed we will see that:

• The fine model is non-strictly hyperbolic and non-conservative.

• The source term is non-zero over 3 equations.

• The entropy is not strictly convex.

• The Chapman-Enskog expansion is rather complicated and requires adding some assumptions to the system (in particular on the thermodynamics of the two phases).

Adaptation with θ = 10 -2 Adaptation with θ = 1

Adaptation with θ = 10 -2 Adaptation with θ = 1 • The numerical scheme used for the coarse model is not derived from the one of the fine model. Thus vn+1 1,K is not deduced from a Chapman-Enskog expansion applied to the numerical scheme but from a direct discretization of the continuous expression of the first order term v 1 .

• The adaptive method has to be modified in order to take into account the non-compatibility of the numerical schemes.

5.1. Fine and coarse models. Following [START_REF] Dellacherie | Relaxation schemes for the multicomponent Euler system[END_REF] for instance, the governing set of equations reads in 1D ( 44)

                                     ∂ t α 1 + u I ∂ x α 1 = λ p (p 1 -p 2 ), ∂ t (α 1 ρ 1 ) + ∂ x (α 1 ρ 1 u 1 ) = 0, ∂ t (α 1 ρ 1 u 1 ) + ∂ x (α 1 ρ 1 u 2 1 + α 1 p 1 ) -p I ∂ x α 1 = -λ u (u 1 -u 2 ), ∂ t (α 1 ρ 1 E 1 ) + ∂ x ((α 1 ρ 1 E 1 + α 1 p 1 )u 1 ) -p I u I ∂ x α 1 = -λ T (T 1 -T 2 ) -u I λ u (u 1 -u 2 ) -p I λ p (p 1 -p 2 ), ∂ t (α 2 ρ 2 ) + ∂ x (α 2 ρ 2 u 2 ) = 0, ∂ t (α 2 ρ 2 u 2 ) + ∂ x (α 2 ρ 2 u 2 2 + α 2 p 2 ) -p I ∂ x α 2 = λ u (u 1 -u 2 ), ∂ t (α 2 ρ 2 E 2 ) + ∂ x ((α 2 ρ 2 E 2 + α 2 p 2 )u 2 ) -p I u I ∂ x α 2 = λ T (T 1 -T 2 ) + u I λ u (u 1 -u 2 ) + p I λ p (p 1 -p 2 ). Here ρ k , u k , p k , E k = e k + 1 2 u 2
k and T k denote respectively the density, the velocity, the pressure, the total energy (e k being the internal energy) and the temperature of each phase with k = 1, 2. The volume fractions α k satisfy α 1 + α 2 = 1, so that the two phases are immiscible. Each phase is characterized by its own EoS. We assume for simplicity that both phases are described by perfect gas laws

p k = π k (e k , ρ k ) = (γ k -1)ρ k e k ,
where γ k > 1 denotes the polytropic exponent. In that case the temperatures satisfy (after normalization C v = 1)

T k = T k (e k , ρ k ) = e k .
The interaction between phases is described with both differential terms and relaxation terms arising in the mass, momentum and total energy equations, where positive relaxation coefficients λ p , λ u and λ T appear. The differential terms are non-conservative products that involve an interfacial pressure p I and an interfacial velocity u I . These terms depict the coupled evolution of the phases: if ∂ x α k = 0 and the source terms are neglected then the two phases are totally decoupled and each one follows the classical gas dynamics equations. A usual choice for the interfacial pressure and velocity is

u I = u 1 , p I = p 2 .
Other interfacial closure relations can be chosen, we refer to [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF] for detailed explanations.

Let us rewrite the seven-equation model under a form similar to (5-6):

(45)

∂ t U + ∂ x f 1 (U, V ) = 0, ∂ t V + ∂ x f 2 (U, V ) + g(U, V )∂ x α 1 = r(U, V ),
where

U = (α 1 ρ 1 , ρ, ρu, ρE) T , V = (α 1 , α 1 ρ 1 u 1 , α 1 ρ 1 E 1 ) T , r(U, V ) =   λ p (p 1 -p 2 ) -λ u (u 1 -u 2 ) -λ T (T 1 -T 2 ) -u I λ u (u 1 -u 2 ) -p I λ p (p 1 -p 2 )   , f 1 (U, V ) =       α 1 ρ 1 u 1 ρu ρu 2 + p + α 1 ρ 1 α 2 ρ 2 ρ (u 1 -u 2 ) 2 (ρE + p)u + α 1 ρ 1 α 2 ρ 2 ρ (u 1 -u 2 ) E 1 + p1 ρ1 -E 2 -p2 ρ2       , f 2 (U, V ) =   0 α 1 ρ 1 u 2 1 + α 1 p 1 (α 1 ρ 1 E 1 + α 1 p 1 )u 1   , g(U, V ) =   u 1 -p 2 -p 2 u 1   , and (46)
ρ = α 1 ρ 1 + α 2 ρ 2 , ρu = α 1 ρ 1 u 1 + α 2 ρ 2 u 2 , ρE = α 1 ρ 1 E 1 + α 2 ρ 2 E 2 , p = α 1 p 1 + α 2 p 2 .
The relaxation terms reflect the fact that the phases are not at mechanical, kinematic and thermodynamical equilibrium. When the relaxation times go to 0 (i.e. when the time-scales λ p , λ T and λ u go to +∞), the source term r(U, V ) vanishes and the system reaches the equilibrium characterized by [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF]

p 1 = p 2 = p, u 1 = u 2 = u, T 1 = T 2 = T.
Because we consider a perfect gas mixture, the equilibrium of temperatures is equivalent to [START_REF] Jin | A domain decomposition method for semilinear hyperbolic systems with two-scale relaxations[END_REF] e 1 = e 2 = e.

Thus the map V eq (such that r(U, V ) = 0 ⇔ V = V eq (U )) can be written as an explicit function of U :

(49) V eq (U ) =        (γ 1 -1)α 1 ρ 1 (γ 2 -1)(ρ -α 1 ρ 1 ) + (γ 1 -1)α 1 ρ 1 (α 1 ρ 1 ) ρu ρ (α 1 ρ 1 ) ρE ρ        .
This asymptotic limit defines the coarse model ( 50)

         ∂ t (α 1 ρ 1 ) + ∂ x (α 1 ρ 1 u) = 0, ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = 0, ∂ t (ρE) + ∂ x ((ρE + p)u) = 0,
which corresponds to the classical Euler equations coupled with the evolution equation of the volume fraction α 1 . Combining ( 47) and ( 48), we deduce that the equilibrium pressure can also be written as a γ-gas law p = (γ(U ) -1)ρe with [START_REF] Lemou | Micro-macro schemes for kinetic equations including boundary layers[END_REF] γ

(U ) -1 = (γ 2 -1) ρ -α 1 ρ 1 ρ + (γ 1 -1) α 1 ρ 1 ρ .
5.2. Discretization of each model. We now address the approximation of the fine model [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF]. The approximation of the seven-equation model has been the subject of many contributions. The numerical scheme we use is adapted from [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF]. It consists in treating the convective terms by the Rusanov scheme [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves with barriers[END_REF] that handles the non-conservative terms and the relaxation terms by a fractional step approach. Let Z n K = (U n K , V n K ) be the one dimensional solution in the cell K at time t n which we want to update, until time t n,1 . The non-conservative Rusanov scheme writes ( 52)

∆x (Z) n,1 K -(Z) n K + ∆t (F f ) n K+1/2 -(F f ) n K-1/2 -∆t(φ) n K (α 1 ) n K+1/2 -(α 1 ) n K-1/2 = 0,
where

F f = f 1 f 2
, f 1 and f 2 being defined in [START_REF] Jaouen | Étude mathématique et numérique de stabilité pour des modèles hydrodynamiques avec transition de phase[END_REF], and

φ = (0, 0, 0, 0, u 1 , p 2 , p 2 u 1 ) T , 2(F f ) n K+1/2 = (F f ) n K + (F f ) n K+1 -r K+1/2 ((Z) n K+1 -(Z) n K ), r K+1/2 = max(r K , r K+1 ), r K = max(|(u 1 ) n K |, |(u 1 ) n K | + (c 1 ) n K , |(u 2 ) n K | + (c 2 ) n K ), 2(α 1 ) n K+1/2 = (α 1 ) n K + (α 1 ) n K+1 .
Here c k , k = 1, 2, denotes the speed of sound of the phase k, which reads

c k = γ k p k ρ k ,
the phase k following a perfect gas law. The second step consists in computing the source terms, separating the velocity, the pressure and the temperature relaxations and thus introducing three intermediate times t n,1 , t n,u and t n,p between t n and t n+1 . We recall the main ideas, the whole procedure being widely depicted in [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF]. Let Z n,u be the approximation of the solution after the velocity relaxation with Z n,1 as initial condition. The velocity relaxation only acts on the velocities u k and the total energies E k , k = 1, 2:

∂ t (α 1 ρ 1 ) = ∂ t ρ = ∂ t (ρu) = ∂ t (ρE) = 0, ∂ t (α 1 ρ 1 u 1 ) = -λ u (u 1 -u 2 ), ∂ t (α 1 ρ 1 E 1 ) = -λ u (u 1 -u 2 ), which is equivalent, for k = 1, 2, to α k ρ k ∂ t u k = (-1) k λ u (u 1 -u 2 ), α k ρ k ∂e k ∂p k ∂p k ∂t = (-1) k λ u u 1 (u 1 -u 2 ).
Applying the implicit Euler method to approximate the previous ordinary differential equations leads to the following approximation of the state vector in the cell K between time t n,1 and t n,u

(α k ) n,u K = (α k ) n,1 K , (α k ρ k ) n,u K = (α k ρ k ) n,1 K , (u k ) n,u K = ((α l ρ l ) n,1 K + λ u ∆t)(α k ρ k u k ) n,1 K + λ u ∆t(α l ρ l u l ) n,1 K (α 1 ρ 1 ) n,1 K (α 2 ρ 2 ) n,1 K + ρ n,1 K λ u ∆t (p 1 ) n,u K = (p 1 ) n,1 K , (p 2 ) n,u K = (p 2 ) n,1 K + λ u ∆t (γ 2 -1)(α 2 ) n,1 K ((u 1 ) n,u K -(u 2 ) n,u K ) 2 .
The pressure relaxation term is now taken into account. From time t n,u to time t n,p the solution is computed by solving the following ODE system

∂ t (α 1 ρ 1 ) = ∂ t (α 1 ρ 1 u 1 ) = 0, ∂ t ρ = ∂ t (ρu) = ∂ t (ρE) = 0, ∂ t α 1 = λ p (p 1 -p 2 ), α 1 ρ 1 ∂ t E 1 = -p 2 λ p (p 1 -p 2 ).
Following [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF] we use an explicit form of λ p λ p = 1 λp

α 1 α 2 p 1 + p 2 ,
where λp is constant.

Integrating the equations on the total energies gives

(p k -p l )(t) = (p k -p l ) n,u exp -λ p t 0 (A k -A l ) (τ )dτ , (p 1 p 2 )(t) = (p 1 p 2 ) n,u exp -λ p t 0 A 1 p 1 (p 1 -p 2 ) + A 2 p 2 (p 2 -p 1 ) (τ )dτ with A k = ∂p k ∂ρ k ρ k α k - ∂p k ∂e k p 2 α k ρ k .
Then multiplying the mass fraction equation by 1 α1α2 and integrating leads to

α 1 1 -α 1 (t) = α 1 1 -α 1 n,u exp 1 λp t 0 p 1 -p 2 p 1 + p 2 (τ )dτ .
Finally we address the relaxation in temperature and solve the ODE system below between t n,p and t

n+1 ∂ t α 1 = ∂ t (α 1 ρ 1 ) = ∂ t (α 1 ρ 1 u 1 ) = 0, ∂ t ρ = ∂ t (ρu) = ∂ t (ρE) = 0, α 1 ρ 1 ∂ t E 1 = -λ T (T 1 -T 2 ).
The two phases follow a perfect gas law so that T k = e k , k = 1, 2. Thus the previous equations are equivalent to

α k ρ k ∂ t e k = (-1) k λ T (e k -e l ), k = l.
These ordinary differential equations are approximated by

(e k ) n+1 K = ((α l ρ l ) n,p K + λ T ∆t)(α k ρ k e k ) n,p K + λ T ∆t(α l ρ l e l ) n,p K (α 1 ρ 1 ) n,p K (α 2 ρ 2 ) n,p K + λ T ∆tρ n,p K , for k = l.
Remark 5.1. Note that this strategy preserves the positivity of the pressure p k , the temperature T k and the maximum principle on the mass fraction α k ∈ [0, 1] (see [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]).

In Section 3 the numerical scheme of the coarse model was deduced from the one of the fine model using the asymptotic preserving property of the fine scheme. In the present case the coarse model has only four equations and is conservative contrary to the fine model. The idea is to use a different numerical scheme to approximate the coarse model [START_REF] Klar | An asymptotic preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF]. In the following we consider the classical Rusanov [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves with barriers[END_REF] scheme for the approximation of the coarse model. Let U n K be the one dimensional solution in the cell K at time t n to be advanced at time t n+1 . The finite volume scheme for the coarse model ( 49) is

(53) ∆x (U ) n+1 K -(U ) n K + ∆t (F c ) n K+1/2 -(F c ) n K-1/2 = 0, together with F c =     α 1 ρ 1 u 1 ρu ρu 2 + p (ρE + p)u     , and (54) 2
(F c ) n K+1/2 = (F c ) n K + (F c ) n K+1 -s K+1/2 ((U ) n K+1 -(U ) n K ), s K+1/2 = max(s K , s K+1 ), s K = max(|(u) n K |, |(u) n K | + (c) n K )
, where c = γp/ρ denotes the speed of sound of the fluid. 5.3. Indicator and adaptive method. Up to now, the adaptive method relied on a numerical indicator which is partly deduced from a discrete Chapman-Enskog expansion. Here the complexity of the fine scheme prevents such computations. To overcome the problem we propose to use a direct discretization of the calculations of [START_REF] Dellacherie | Relaxation schemes for the multicomponent Euler system[END_REF]. The author performs the Chapman-Enskog expansion on a two-phase model close to [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF]. Let us recall the result, detailed calculations are given in [START_REF] Dellacherie | Relaxation schemes for the multicomponent Euler system[END_REF].

Proposition 5.1. Assume that the relaxation process has only one time scale: λ u = λ p = λ T = 1 ε , ε > 0 fixed. Up to terms of order ε 2 , the smooth solutions of (44) satisfy

A = ρ(Y 1 ) 2 Y 2 ρ ρ 1 -1 ∂ x p, B = e(Y 1 -Y 2 ) 2 ∂ x u, C = ρeY 1 Y 2 ∂ x p γ 1 Y 1 ρ ρ 1 -1 + γ 2 Y 2 ρ ρ 2 -1 + uB. (56)
Note that the non-conservative terms p 2 ∂ x α k , p 2 u 1 ∂ x α k and the temperature relaxation term are mandatory to compute the Chapman-Enskog expansion. Moreover the first order terms in ε are explicit functions of U . decreasing (as soon as the underlying numerical schemes possess these features).

On the basis of these properties, one may hope to prove convergence results, for instance in the simplified framework of [START_REF] Coquel | Well-posedness and singular limit of a semilinear hyperbolic relaxation system with a two-scale discontinuous relaxation rate[END_REF].

Mainly three test cases are presented in order to assess the relevance of the indicator, in the context of dynamic model adaptation. When the model completely fits the classical assumptions on hyperbolic systems with relaxation, our strategy provides very satisfying results. The use of the fine model only coincides with the strong variations of the solution. The two other test cases are much more challenging since they do not fully enter the initial framework. But here again, the indicator allows an accurate dynamic partition of the domain and the adaptive method provides in most cases very close results to the reference solution of the fine model. However, let us recall that our foremost interest in this work concerns the development of an indicator to determine when and where a fine model can be replaced by an associated coarse model, and not yet to construct a multiscale method to save CPU time. In order to hope any decrease of the computational cost, the gap between the fine and the coarse must be important, with complex or tabulated laws to be solved by iterative methods, for the fine model.

The main direction of improvement of our adaptive algorithm concerns the development of error estimates. Indeed, we do not take into account the error due to the coupling and do not really compare the error of model adaptation with the numerical error. Secondly, the final goal of adaptation is to be able to control the threshold in order to obtain a final accumulated error which is a priori given by the user. This would necessitate a posteriori estimates, which are very difficult to have for hyperbolic systems of balance laws. Several works on these topics are under investigation. In particular, a similar adaptive method in the scalar setting is studied, for which rigorous convergence results can be obtained, see [START_REF] Cancès | Error analysis of a dynamic model adaptation procedure for nonlinear hyperbolic equations[END_REF].

Solutions for ε = 10 to detect the different waves, but it vanishes in the left part of the domain. The main consequence is that the results are very similar to those with smaller ε and are far from the expected solution. On the contrary, the indicator (59) works well at the left of the domain where the solution is constant but far from the equilibrium, but poorly resolves the front at the center of the domain, see Fig. 13 (top left).

Clearly, the adaptive method with the full indicator [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] takes advantage of both parts, to produce results which are very similar to those of the fine model.

4. 1 . 1 .

 11 Shock tube test case. We are focusing on a shock tube for the model (39) with ε = 2.10 -2 , with classical homogeneous Neumann conditions at the boundaries of the domain[-1, 1]. The initial data is (τ, u, T )(0, x) =

Figure 1 .

 1 Figure 1. Shock tube: specific volume τ (left), perturbed specific volume T (right).

Solutions at time T 1 = 6 -

 16 Characteristic function of the fine domainSolutions at time T 3 = 0.

Figure 2 .

 2 Figure 2. Shock tube: velocity u (left), characteristic function of the fine model and numerical indicator (right).

Figure 3 .

 3 Figure 3. Shock tube: L 1 error (in Log-Log scale) between the fine model and the adaptation algorithm w.r.t. θ: perturbed specific volume T .

Figure 4 .

 4 Figure 4. Shock bubble interactions. Solution after the first interaction: mass fraction ϕ

1 Figure 5 .

 15 Figure 5. Shock bubble interactions. Solution after the first interaction: density ρ

Figure 6 .

 6 Figure 6. Shock bubble interactions. Solution during the second interaction: mass fraction ϕ

Figure 7 .

 7 Figure 7. Shock bubble interactions. Solution during the second interaction: density ρ

Figure 8 .

 8 Figure 8. Partition of the domain between the coarse domain (red) and the fine domain (blue) after the first interaction (up) and during the second interaction (down)

Figure 9 .

 9 Figure 9. Shock bubble interactions: L 1 error (in Log-Log scale) between the fine model and the adaptation algorithm w.r.t. θ: mass fraction ϕ.

( 55 )

 55 ∂ t (α 1 ρ 1 ) + ∂ x (α 1 ρ 1 u) = ε∂ x A, ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = ε∂ x B, ∂ t (ρE) + ∂ x ((ρE + p)u) = ε∂ x C, where, defining Y k = α k ρ k ρ ,

Solutions at time T 1 =Figure 10 .

 110 Figure 10. Relaxation test case: relative pressure (left), relative temperature, zoom on [0, 0.3] (right).

Figure 12 .

 12 Figure 12. Relaxation test case. L 1 error (in Log-Log scale) between the fine model and the adaptation algorithm w.r.t. θ: V -V eq (U ) .

Figure 13 .

 13 Figure 13. Solutions with the indicators (27), (58) and (59), for different values of ε: perturbed specific volume T (left) and the indicators (58) and (59) (right).

  -right.

				Solutions at time T 1 = 0.2
	1			Fine model	1				Fine model
	0.98			Coarse model Adaptation algorithm	0.98				Coarse model Adaptation algorithm
	0.96				0.96
	0.94				0.94
	0.92				0.92
	0.9				0.9
	0.88				0.88
	0.86				0.86
	0.84				0.84
	0.82				0.82
	0.8				0.8
	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8	1	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8
				Solutions at time T 2 = 0.4
	1			Fine model	1				Fine model
	0.98			Coarse model Adaptation algorithm	0.98				Coarse model Adaptation algorithm
	0.96				0.96
	0.94				0.94
	0.92				0.92
	0.9				0.9
	0.88				0.88
	0.86				0.86
	0.84				0.84
	0.82				0.82
	0.8				0.8
	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8	1	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8
				Solutions at time T 3 = 0.6
	1			Fine model	1				Fine model
	0.98			Coarse model Adaptation algorithm	0.98				Coarse model Adaptation algorithm
	0.96				0.96
	0.94				0.94
	0.92				0.92
	0.9				0.9
	0.88				0.88
	0.86				0.86
	0.84				0.84
	0.82				0.82
	0.8				0.8
	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8	1	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8
				Solutions at time T 4 = 0.8
	0.94			Fine model	0.94				Fine model
				Coarse model					Coarse model
	0.93			Adaptation algorithm	0.93				Adaptation algorithm
	0.92				0.92
	0.91				0.91
	0.9				0.9
	0.89				0.89
	0.88				0.88
	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8	1	-1	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8

Acknowledgments. The authors would like to thank Frédéric Coquel who initially pointed out the issue of model adaptation to them and for fruitful discussions. They also thank the reviewers, in particular for the helpful suggestion which

Then we propose to directly use the L 1 -norm of a discrete approximation of the corrector of first order in ε (56) as the numerical indicator (without using the fine numerical scheme as it is done in Section 3). The numerical indicator thus reads [START_REF] Saleh | Analyse et simulation numérique par relaxation d'écoulements diphasiques compressibles. Contribution au traitement des phases évanescentes[END_REF] vn+1

together with

where

We now address the adaptive method. Because the coarse scheme is not derived from the fine scheme, the adaptive method has to be modified.

We recall that P 1 is the linear operator such that P 1 Z = U and M is the map such that M (U ) belongs to the equilibrium manifold M (see Section 2.1). Let

) be the solution in the cell K known at time t n which we want to update, until time t n+1 . Let B n K be the balance of fluxes in the cell K at time t n , i.e.

The 1D adaptive method reads as follows: A) For all cell K, compute the numerical error E n+1 K using ( 27) with ( 57

) n K+1/2 using the numerical scheme [START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF]. Update the balance

-Else -Compute (F f ) n K+1/2 using the numerical scheme [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. Update the balance

5.4. Relaxation test case. We apply this algorithm in the case of the evolution of a two-phase flow in a duct of uniform section and 7m length. The simulation is performed on 1000 cells. The two phases are depicted by a perfect gas law with γ 1 = γ 2 = 1.4. The flow is initially at equilibrium, the initial conditions are α 1 = 0.9, u 1 = u 2 = 5, p 1 = p 2 = 10 5 , ρ 1 = ρ 2 = 1. An inlet Dirichlet boundary condition is applied on the left side using the initial condition except that p 1 = 1.005 × 10 5 and the relaxation parameters are λ u = λ p = λ T = ε -1 = 10 4 . An outlet Neumann boundary condition is applied on the right side. These boundary conditions are taken into account using the classical ghost cell method. Finally, the adaptation parameter is θ = 10 -1 . This test case allows us to study the asymptotic behavior of the two-phase model. Due to the relaxation terms in the fine model, the relative quantities tend to zero as x increases, leading to a relaxation boundary layer on the left. It is worth noting that only the boundary relative pressure is non-zero but the boundary layer also appears for the relative velocity and for the relative temperature.

On the other hand, when considering this test case with the coarse model, all the relative quantities are null, by definition of the equilibrium manifold [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF].

We present the results obtained for four successive times: t = 0.05, 0.1, 0.15,and 0.2. The relative pressure |p 1p 2 |/max(p 1 , p 2 ) and the relative temperature

and the characteristic function of the fine domain with the numerical indicator are plotted in Fig. 11. Due to the remark here-above, the results given by the coarse model are always constant in space and equal to 0.

Once again, the results provided by the adaptation algorithm are very close to those of the fine model. However, one can observe a small overshoot on the relative pressure at time 0.05, see Fig. 10 (top left). Note that the results for the relative temperature are zoomed in space since the associated boundary layer is very small, see Fig. 10 (right). When considering the decomposition of the domain in Fig. 11 (right), one can see that the fine model is only solved at the left part of the computational domain, which corresponds to the boundary layer.

Relaxation test case: convergence with respect to the threshold.

As for the previous test cases, we study the convergence of the solutions of the adaptation algorithm towards those of the fine model. Contrary to the previous cases, the full convergence is not reached. Indeed, one can see in Fig. 12 that for θ < 10 -7 , the error does not decrease anymore. This means that the numerical indicator is not fully relevant in this case, which is not totally surprising since it is based on a direct discretization of the first order term v 1 of the Chapman-Enskog expansion. Therefore, it is not completely compatible with the numerical schemes we use for the fine and the coarse models.

Conclusion

We propose an original indicator which can be used to place in an optimal way coupling interfaces between asymptotically compatible models. Due to the applications under investigation, we are not allowed to modify the numerical schemes for computing both models, which prevents us from using either classical asymptotic preserving schemes on the whole domain or sophisticated micro-macro decompositions. Our method relies on a numerical Chapman-Enskog expansion and with an estimate of the distance to the equilibrium manifold. We then develop an adaptive method which dynamically and locally selects the right model to solve. Using an appropriate coupling method at the interface between the different models, we are able to prove that our method exactly maintains constant equilibrium states, but also provides numerical solutions in the set of admissible states and is entropy

Appendix A. Comparison of different indicators

In this appendix, we provide a numerical comparison of the different components of the numerical indicator ( 27), ( 27)

To do so, we compare on a simple test case the numerical solutions provided by the adaptation algorithm with different indicators. We compare the indicator ( 27) with

E n+1 dv,K := v n Kv eq (u n K ) . On the one hand, for data which do not belong to the equilibrium manifold M and are constant (w.r.t. space), the indicator E n+1 CE,K is not sufficient, since:

, where W 0 is a constant state which is not at equilibrium. In such case, if the adaptation algorithm is only based on the indicator E n+1 CE,K , then only the coarse model will be used, on the whole computational domain. Therefore, at the next time step, the solution will directly be at equilibrium:

, which is of course not consistent with the true solution.

On the other hand, for data which belong to the equilibrium manifold M, the indicator E n+1 dv,K is not relevant, since: [START_REF] Tzavaras | Relative entropy in hyperbolic relaxation[END_REF] ∀K, W n K ∈ M =⇒ ∀K, E n+1 dv,K = 0. As a consequence, the adaptation algorithm based on the indicator E n+1 dv,K would only use the coarse model, whatever the smoothness of the solution is, and then would allow shock formation. Since in many cases relaxation terms have a dissipative effect [START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF][START_REF] Yong | Entropy and global existence for hyperbolic balance laws[END_REF], such discontinuous solutions are impossible. Therefore, the indicator ( 59) is not relevant either.

In order to illustrate this discussion, we provide here some numerical tests with these different indicators. Moreover, we compare different values of ε, since the adaptive method is a priori more useful when ε is large (i.e. when the solutions provided by the fine and the coarse models are very different).

We consider the Suliciu model with Chaplygin gas, which is presented in Section 4.1. The domain is [-1, 1] with homogeneous Neumann conditions and the results are plotted at t = 0.2. The initial data is composed by two constant states:

Note that only the state at the right-hand side is at equilibrium. We present three results, associated with three values of ε: 1, 10 -2 and 10 -4 . We plot the results obtained with the adaptive method with the different indicators. The threshold θ is set to a value for which the solution of the adaptive method is close to the solution of the fine model (respectively 0.2, 10 -3 and 5.10 -6 ). We plot in Fig. 13 the variable T obtained by the fine model and the adaptive methods with each indicator. We also present at the right-hand side the indicators (58) and (59) (the full indicator [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] can be deduced by taking the max). First of all, one can check that for the smallest value of ε, all the numerical indicators lead to similar results. This is not surprising since, in this case, the solution provided by the fine model and the one provided by the coarse model are very close. For the intermediate value of ε (= 10 -2 ), Even if the shape of the solutions are comparable, the indicator (58) leads to underestimated results, while the indicator [START_REF] Serre | Multidimensional shock interaction for a Chaplygin gas[END_REF] gives greater results that those of the full indicator and of the fine model. One can also see that the indicators are very similar in Fig. 13 (middle right). In the last case with ε = 1, one recover the behaviors described by the two examples [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions[END_REF] and [START_REF] Tzavaras | Relative entropy in hyperbolic relaxation[END_REF]. Indeed, the indicator based on the Chapman-Enskog expansion is able led us to improve the numerical indicator [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF]. This work has been supported by the LRC Manon (Modélisation et approximation numérique orientées pour l'énergie nucléaire -CEA/DM2S-LJLL).