FR ÉD ÉRIC COQUEL Benjamin Boutin

Frédéric Coquel
email: coquel@cmap.polytechnique.fr

AND Philippe G Lefloch

COUPLING TECHNIQUES FOR NONLINEAR HYPERBOLIC EQUATIONS. IV. WELL-BALANCED SCHEMES FOR SCALAR MULTIDIMENSIONAL AND MULTI-COMPONENT LAWS

Coupling techniques for nonlinear hyperbolic equations. IV. Well-balanced schemes for scalar multi-dimensional and multi-component laws

Introduction

Objective of this paper. This is a continuation of [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF][START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF][START_REF] Boutin | Coupling techniques nonlinear hyperbolic equations. III. A regularization method based on thick interfaces[END_REF] devoted to coupling techniques for nonlinear hyperbolic equations. In the present paper, we deal with the coupling of multi-dimensional hyperbolic equations, based on an arbitrary partition of the physical domain. The main motivation stems from the study of complex systems resulting from the combination of elementary components modeled by different equations. Indeed, each component may be subject to physical phenomena involving fairly different time and space scales. Tackling this multiscale problem with sufficient accuracy and efficiency requires to consider distinct physical models for the description of each component, so as to end up with a suitable description of the whole physical system. For instance, large-scale power plants provide a typical example of interest [START_REF] Guelfi | NEPTUNE: A new software platform for advanced nuclear thermal hydraulic[END_REF]. Describing the evolution in time requires the exchange of transient informations at each physical boundary separating two distinct hyperbolic models. These transient informations or boundary conditions are referred hereafter to as coupling conditions.

This problem seems to be rather new in the applied mathematical community. Its analysis was initiated by Godlewski and Raviart [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF] for scalar equations in one space variable. Therein, the coupling problem is formulated in terms of two initial boundary value problems (IBVP) supplemented with coupling boundary conditions at a given (infinitely thin) interface. These boundary conditions are stated in such a way that in "most cases" they ensure the continuity of the main unknown, at least, roughly speaking, as long as no wave from the left-and right-hand problems interact at the interface. If this condition does not hold, one says that the interface is resonant. In Ambroso et al. [START_REF] Ambroso | Coupling of general Lagrangian systems[END_REF][START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF][START_REF] Ambroso | Relaxation methods and coupling procedures[END_REF], quite general continuity conditions based on a nonlinear transformation of the unknown were investigated. Following earlier investigations by LeFloch and collaborators [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF][START_REF] Joseph | Boundary layers in weak solutions to hyperbolic conservation laws[END_REF][START_REF] Lefloch | Propagating phase boundaries. Formulation of the problem and existence via the Glimm scheme[END_REF][START_REF] Lefloch | Hyperbolic Systems of Conservation Laws. The theory of classical and nonclassical shock waves[END_REF][START_REF] Lefloch | Kinetic relations for undercompressive shock waves. Physical, mathematical, and numerical issues[END_REF][START_REF] Lefloch | theory for nonlinear hyperbolic systems in nonconservative form[END_REF] on undercompressive shocks and interfaces, nonconservative hyperbolic systems, and boundary value problems, we stress that additional information coming from physical modeling is necessary in order to single out the relevant continuity conditions (or transmission condition) at the interfaces. Various conditions were introduced and studied in a variety of physical frameworks, ranging from gas dynamics [START_REF] Ambroso | Coupling of general Lagrangian systems[END_REF] to multiphase flows [START_REF] Ambroso | The drift-flux asymptotic limit of barotropic two-phase two-pressure models[END_REF][START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF].

Thin interface versus thick interface. We briefly mention some transmission conditions of interest when the coupling invoves two Euler systems with distinct pressure laws. Typically, one imposes the continuity of the density ρ, velocity component u, and pressure p, or else the continuity of the convervative variables (ρ, ρu, ρE) (where E denotes the total energy). These conditions determine the class of constant solutions in the time and the space variables, and have either constant density, velocity, and pressure, or else constant density, momentum, and total energy. In both cases, the proposed coupling conditions are nonconservative, since the total mass of density, momentum, and total energy do vary with time. A fully conservative coupling may turn relevant in some applications, as was addressed in [START_REF] Ambroso | Relaxation methods and coupling procedures[END_REF] (following [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF]) via suitable a relaxation method.

The resonance phenomena, likely to take place around thin interfaces, brings a main difficulty in the mathematical analysis of coupled initial boundary value problems. Solutions can be shown to exist under general conditions but resonance generally comes at the expense of uniqueness. We refer the reader to [START_REF] Boutin | Existence result for the coupling problem of two scalar conservation laws with Riemann initial data[END_REF] for a discussion of scalar equation and to [START_REF] Ambroso | Coupling of general Lagrangian systems[END_REF] for a distinct behavior exhibited for characteristic but non-resonant interfaces. A selection criterion for discontinuous solutions, therefore, is required. Recall that, for the fully conservative coupling, several distinct entropy criteria have been proposed, each selecting a distinct weak solution in agreement with the physical context. (See [START_REF] Bürger | Conservation laws with discontinuous flux: a short introduction[END_REF] for a review and [START_REF] Seguin | Analysis and approximation of a scalar conservation law with a flux function with discontinuous coefficients[END_REF][START_REF] Helluy | Relaxation models of phase transition flows, M2AN Math[END_REF][START_REF] Bachmann | Existence uniqueness of entropy solution of scalar conservation laws with a flux function involving discontinuous coefficients[END_REF]).

To deal with general transmission conditions, a macroscopic selection principle analogous to the entropy inequalities is not available and one needs a detailed description of microscopic mechanisms coming with suitable regularizing procedures. In [START_REF] Boutin | Dafermos's regularization for interface coupling of conservation laws[END_REF][START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF][START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF], we introduced an alternative modeling for the coupling problem associated with two hyperbolic equations in one space variable. This alternative method relies on the introduction of an augmented PDE (partial differential equations) formulation that avoids the need of a detailed description of the interfaces. The proposed formalism is based on an additional unknown, the color function which takes values in the range [0, 1]. Extreme values 0 and 1 are devoted to restore the left-and right-problems to be coupled, while intermediate values may serve to model a smooth transition from one problem to the other.

Outline of this paper. The interest in this augmented formulation comes from its very capability to support various regularization mechanisms. Viscous perturbations were introduced by the authors [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF][START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF] for scalar problems and, specifically, a self-similar approach was developed, which allows for the study of the existence and uniqueness of solutions to the coupled Riemann problem in the limit of vanishing viscosity. The analysis has been carried out for a general class of systems [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF] and led to an existence theory under fairly general assumptions. In [START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF], the analysis of the internal structure of resonant interfaces was performed and led us to a characterization of the set of admissible Riemann solutions. Despite of the viscous mechanisms a failure of uniqueness may be observed for resonant infinitely thin interfaces.

Riemann solutions may be indeed understood as describing the long time asymptotic of the solutions of the Cauchy problem. Failure of uniqueness for thin interfaces just reflects the property that distinct regularizations of thin interfaces may give rise to different solutions and thus with a distinct long time behavior. This observation has motivated a second regularization procedure based on thick interfaces.

Thick interfaces within the augmented PDE framework are based on a regularization of the discontinuous color function, considered in the thin regime. This approach has been introduced by the authors first within the framework of two coupled conservation laws in one space variable [START_REF] Boutin | Coupling techniques nonlinear hyperbolic equations. III. A regularization method based on thick interfaces[END_REF]. Existence and uniqueness for the coupled Cauchy problem was proven for general initial data with bounded sup-norm. One of the main ingredients of proof was the design of a well-balanced finite volume method. The well-balanced property means that the exact constant solutions selected by a given transmission condition are exactly preserved at the discrete level, whatever choice is made for the regularized color function. This consistency property is of central importance.

In the present paper, we introduce a framework which covers coupling problems in several space variables and with distinct scalar hyperbolic equations, allowing for possible covering in space. An outline of this paper is as follows.

• In Section 2, we show how to extend the two existing coupling frameworks in one space variable to the coupling of two distinct hyperbolic equations in several space variables. We then show how to extend the augmented PDE formalism to encompass the case of several hyperbolic equations with possible covering. In our approach, a vector-valued color map is introduced so that each component is associated with one of the equations and takes values in the interval [0, 1]. The specific definition of the regularized color function provides us with a transition from an equation to another (possibly more than one). • We check the existence and uniqueness of entropy solutions to the coupled Cauchy problem (with initial data in L ∞ ) under fairly general assumptions on the transmission conditions and the equations under consideration. • Next, in Section 3, we design a robust and flexible finite volume framework based on general triangulations. Importantly, by construction, the proposed method is well-balanced and our strategy for achieving the well-balanced property is an extension of the subcell reconstruction approach (analyzed by Bouchut in a different context [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]). In particular, we introduce two distinct meshes: the first one, the primal mesh, describes the main coupled unknown. The second mesh, referred to as the dual one, is built from the primal mesh and carries the approximation of the color function. A comprehensive derivation of this dual mesh is also proposed.

• In Sections 4 and 5, we then derive a sup-norm estimate, and observe that a uniform estimate on the total variation seems to be out of reach, due to the subcell reconstruction procedure. Consequently, we propose to use DiPerna's framework based on entropy measure-valued solutions and, by deriving suitable entropy inequalities and entropy dissipation bounds, we establish the strong convergence of the proposed method. • Finally, in Section 6, numerical experiments are presented which concern problems with covering in space and, therefore, highlight the interest of the new coupling strategy.

2.

A framework for multi-dimensional coupling 2.1. Coupling of two systems.

Pasting together two initial boundary value problems. In this section, we introduce the coupling problem associated with two hyperbolic equations coupled at a given interface. At this stage, it suffices to think of an hyperplane, say {x 1 = 0}. We will extend two distinct coupling strategies that have been developed in a single space variable. The first procedure consists in modeling the coupling problem as two initial boundary value problems (IBVP) with time dependent boundary conditions prescribing the evolution of traces of the coupled solutions on both sides of the hyperplane {x 1 = 0}. In contrast, the second strategy introduced in [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF][START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF][START_REF] Boutin | Coupling techniques nonlinear hyperbolic equations. III. A regularization method based on thick interfaces[END_REF] is based on augmented PDE systems, and handles the coupling problem as an initial data problem written over the entire space R d . This new framework brings mathematical and numerical advantages, pointed out at the end of this section. Consider an hyperplane of R d with unit normal vector ν ∈ R d , we denote

H = {x ∈ R d /x.ν = 0}, partitioning R d into two half-domains D -= {x ∈ R d /x.ν < 0} and D + = {x ∈ R d /x.ν > 0}.
In each open subdomain, a distinct conservation law is prescribed:

(2.1) ∂ t w + d i=1 ∂ xi a ± i (w) = 0, w(t, x) ∈ R, t > 0, x ∈ D ± ,
where the flux-functions A ± : R → R d , with components (a ± i ) i=1,...,d , are assumed to be twice differentiable for definiteness. An initial data w(0, x) = w 0 (x) supplements this formulation, but obviously some extra-condition, the coupling condition, must be prescribed at the interface H. For simplicity, we restrict ourselves in this introductory section to piecewise smooth solutions w with bounded left and right traces at the interface H: w(t, y±) := lim z→0+ w(t, y ± zν), y ∈ H. Then, it sounds natural that the coupling condition we seek should relate these traces

(2.2) C(w(t, y-), w(t, y+)) = 0, t > 0, y ∈ H,
for some nonlinear mapping C to be specified. The implicit function theorem is assumed to apply so as to recast (2.2) in the more tractable form w(t, y-) = c(w(t, y+)), t > 0, y ∈ H, for some function c mapping R onto R. Assuming from now on c to be strictly monotone, we re-express the above coupling condition in terms of two nonlinear monotone functions θ -and θ

+ with c = θ -1 -• θ + : (2.3) θ -(w(t, y-)) = θ + (w(t, y+)), t > 0, y ∈ H.
Here and without loss of generality, θ -and θ + are assumed to be strictly increasing and to map R onto R, and their inverse functions are denoted by γ -and γ + . On the basis of this pair of functions, we introduce the following useful change of unknown:

(2.4) u(t, x) = θ -(w(t, x)), t > 0, x ∈ D -, θ + (w(t, x)), t > 0, x ∈ D + ,
so that the coupling condition (2.3) resumes to:

(2.5) u(t, y-) = u(t, y+), y ∈ H.

Observe that in the new unknown, (2.5) just reads as a continuity condition for u.

It is worth underlining that (2.5) defines the constant solutions of the coupling problem (2.1)-(2.3), i.e. time independent functions w = w(x) which solve (2.1) and (2.5). Such functions clearly obey

(2.6) u(w(x)) = u , t > 0, x ∈ R d \ H,
for some real u ∈ R (here, the notation u(w) reminds of the change of unknown (2.4)). This observation actually just opens a path toward the mathematical study of perturbed solutions of the trivial solution (2.6). We refer the reader to the work [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF] devoted to the existence of self-similar coupled solutions for systems.

Observe that the coupling condition (2.3) plays the role of a pair of transient boundary conditions for the interface H. In other words, the coupling framework we address merely takes the form of two nonlinear hyperbolic IBVPs linked via the transient boundary condition (2.3). It becomes clear that the coupling condition (2.3) is actually expressed in a strong sense, since it is formulated without reference to the signature of the wave speeds at the interface H. It is nevertheless well-known that the sign of the wave velocities at a boundary directly affects the boundary condition to be prescribed. Hence, the coupling condition (2.3) or its equivalent form (2.5) must be stated in a weak sense.

We follow the approach for coupled problem in one space variables, originally developed in Godlevski, Raviart, and collaborators (cf. [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. II. The case of systems[END_REF][START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF] and [START_REF] Ambroso | Coupling of general Lagrangian systems[END_REF][START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF][START_REF] Chalons | The interface coupling of the gas dynamics equations[END_REF]). In these papers, a weak form of the coupling condition (2.3) is formulated in terms of an admissible boundary set, proposed by Dubois and LeFloch [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF] and based on the notion of Riemann solutions. Such a notion here readily extends since the coupling condition expressed in (2.3) just links the traces of the coupled solution w in the normal direction ν and thus essentially concerns the quasi-one dimensional form of (2.3) written for plane wave solution in the ν-direction. Thus it turns natural to consider the coupled problem in one space variable (up to some shift in the space variable z)

(2.7) ∂ t w + ∂ z A ± ν (w) = 0, t > 0, ±z > 0, where we have set A ± ν (w) = A ± (w) • ν.
In order to state the weak form of the boundary condition θ -(w(t, y-)) = θ + (w(t, y+)), y ∈ H, we first recall the Dubois-LeFloch framework for say the right IBVP:

∂ t w + ∂ z A + ν (w) = 0, t > 0, z > 0, (2.8) w(t, 0+) = b, t > 0, (2.9)
for some prescribed real b. Following Dubois and LeFloch, a weak formulation of (2.9) is stated in terms of Riemann solutions associated with (2.8), that is, W(•; w L , w R ) (for left-and right-hand states w L , w R ):

(2.10)

w(t, 0+) ∈ O + ν (b) = W(0+; b, w), w ∈ R .
Observe that the analogous of (2.10) for the left IBVP built from A - ν would read w(t, 0-) ∈ O - ν (b) = W(0-; w, b), w ∈ R . These considerations naturally yield us to the following coupled boundary conditions (2.3) at any point y ∈ H and for t > 0:

(2.11)

w(t, y+) ∈ O - ν (θ -1 + • θ -(w(t, y-))), w(t, y-) ∈ O + ν (θ -1 -• θ + (w(t, y+))
). This simple problem, based of two coupled equations at a given hyperplane, can be easily extended to more general interfaces resulting from a partition of R d into two non-overlapping open sets D + and D -such that D -∪D + = R d , separated by a smooth boundary ∂D = D -∩ D + . Smoothness allows to define without ambiguity an unit normal vector ν(y) for all y ∈ ∂D so that left and right traces at ∂D for piecewise smooth solutions of the coupled problem (2.3) may be defined as follows:

w(t, y±) = lim z→0+ w(t, y ± zν(y)), y ∈ ∂D.
The expected coupling condition just takes the weak form (2.11).

Coupling technique based on an augmented PDE's system. As already emphasized, an alternative coupling framework has been introduced by the authors in [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF]. Instead of dealing with two IBVPs coupled at a given interface via boundary conditions, our new approach treats the coupling problem as a single initial value problem, over the entire space R d via an augmented PDE formulation. This strategy was introduced by the authors [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF][START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF][START_REF] Boutin | Coupling techniques nonlinear hyperbolic equations. III. A regularization method based on thick interfaces[END_REF] for problems in one space variable. In order to encompass problems in several space variables, we perform hereafter a comprehensive derivation.

The derivation starts from the characteristic functions of the two open sets D - and D + , we denote by v -= χ D-, v + = χ D+ . It heavily makes use of the change of unknown u introduced in (2.4), we rephrase as:

u(t, x) = θ -(w(t, x)), if v -(x) = 1, θ + (w(t, x)), if v -(x) = 0, i.e. if v + (x) = 1, t > 0, x / ∈ ∂D.
Equipped with these notation, we recast the two distinct hyperbolic equations in D ± in terms of u:

γ ± (u)∂ t u + d i=1 γ ± (u)a ± i (γ ± (u))∂ xi u = 0, t > 0,
x ∈ D ± , restricting ourselves to smooth solutions in a first stage. Recall that γ + (respectively γ -) denotes the inverse function of θ + (resp. θ -). We further proceed by rewritting the above two equations in term of a single equation in

x ∈ R d \ ∂D: v -γ -(u)+v + γ + (u) ∂ t u+ d i=1 v -γ -(u)a - i (γ -(u))+v + γ + (u)a + i (γ + (u)) ∂ xi u = 0,
At this stage, it must be noticed that the two characteristic functions v -and v + in the above equation may be replaced by a single function say v, by setting for instance v -(x) = 1v(x) and v + (x) = v(x) for x ∈ R d \ ∂D with v = χ D+ . In the following, such a function v will be refered to as a color function. For the moment v is nothing but a step function taking values in {0, 1} but it is important to conceive v as a function taking values in the interval [0, 1] so that the value 0 restores the equation set in D -while the value 1 restores the equation set in D + . Intermediate values of v then may be thought as modeling a smooth shift from one problem to the other. Keeping this in mind we now recast the equations above in the form of an augmented PDE system with unknown u and v, for t > 0 and x ∈ R d \ ∂D:

(2.12)

(1 -v)γ -(u) + vγ + (u) ∂ t u + (1 -v)γ -(u)∇A -(γ -(u)) + vγ + (u)∇A + (γ + (u)) • ∇ x u = 0, ∂ t v = 0.
We stress that the 1-dimensional form of these equations written for plane wave solutions in the direction ν reads (t > 0, x ∈ R d \ ∂D, or ±z > 0):

(2.13) (1 -v)γ -(u) + vγ + (u) ∂ t u + (1 -v)γ -(u)∇A -(γ -(u)) • ν + vγ + (u)∇A + (γ + (u)) • ν ∂ z u = 0, ∂ t v = 0.
This system is easily seen to be hyperbolic if (and only if) the following quantity is not zero

(2.14) (1 -v)γ -(u)∇A -(γ -(u)) • ν + vγ + (u)∇A + (γ + (u)) • ν = 0.
For such states, the standing wave associated with the additional unknown v can be seen to admit u as a Riemann invariant. In other words, as long as the nondegeneracy condition (2.14) is valid, u stays continuous at the jumps of the color function v, namely across the coupling boundary ∂D at which the value of v shifts from 0 to 1. In other words and whenever (2.14) is valid, the coupling condition (2.5) is satisfied in the strong sense across the standing wave (2.15) u(t, y-) = u(t, y+), y ∈ ∂D.

Violation of the condition (2.14) at a point of jump for v, namely at the interface ∂D, expresses that waves from the left and right propagate with opposite sign at the interface; the first order system (2.13) is then only weakly hyperbolic. This is the resonance phenomena for which we refer the reader to, for instance, Goatin and LeFloch [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] and the references cited therein. As far as the coupling issue is concerned, the continuity condition (2.15) is no longer satisfied and the weak form (2.11) of the coupling condition must be addressed. Turning considering the augmented formulation (2.13), resonance phenomena has been studied in depth in [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF] in the scalar case thanks to a self-similar viscous perturbation. The Riemann solutions for (2.13) defined in the limit of vanishing viscosity satisfy (2.11) when resonance takes place. To sum up, weak solutions of the augmented equations (2.13) and thus their multi-dimensional form (2.12) naturally encode the weak form of the coupling condition. We now generalize the rather special form of the augmented equation and adopt the general framework introduced by the authors in [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF] (which also applies to systems in one space variable). We thus introduce coupling functions C 0 : R × [0, 1] → R and C i : R × [0, 1] → R with i ∈ {1, . . . , d} satisfying the following consistency properties:

(2.16) lim v→0 C 0 (u, v) = γ -(u), lim v→1 C 0 (u, v) = γ + (u), lim v→0 C i (u, v) = a - i (γ -(u)), lim v→1 C i (u, v) = a + i (γ + (u)),
so as to consider in place of (2.12) the general augmented equations:

(2.17)

∂ u C 0 (u, v)∂ t u + d i=1 ∂ u C i (u, v)∂ xi u = 0, ∂ t v = 0, t > 0, x ∈ R d ,
which equivalently recasts as:

∂ t C 0 (u, v) + d i=1 ∂ xi C i (u, v) - d i=1 ∂ v C i (u, v)∂ xi v = 0, ∂ t v = 0. t > 0, x ∈ R d .
In the following, the coupling functions C 0 and (C i ) 1≤i≤d are smooth and

C 0 , (C i ) 1≤i≤d ∈ C 2 (R × [0, 1]),
and C 0 , in addition, obeys

∂ u C 0 (u, v) > 0, u ∈ R, v ∈ [0, 1]
, which is a nondegeneracy condition for the time arrow in (2.17). The resonance phenomenon is the main difficulty in the coupling problematic and has made the matter of previous works especially in the one-dimensional case [START_REF] Boutin | Coupling techniques for nonlinear hyperbolic equations. I. Self-similar diffusion for thin interfaces[END_REF][START_REF] Boutin | Coupling techniques for nonlinear equations[END_REF][START_REF] Boutin | Coupling techniques nonlinear hyperbolic equations. III. A regularization method based on thick interfaces[END_REF]. In this one-dimensional setting, the analysis proves that if the resonance occurs for (2.17) the self-similar weak solutions obtained via self-similar regularization satisfy the coupling relation (2.15). Nevertheless in the general case where resonance may appear, uniqueness then generaly fails for the initial value problem.

The central interest of the augmented formulation (2.17) over more classical coupling approaches built from a collection of IBVPs stems from the fact it can be supplemented with a variety of regularizing mechanisms at the coupling interfaces. These regularization mechanisms are intended to handle the resonance phenomena which is likely to take place at the interfaces. A first regularization procedure relies on introduction of suitable viscous mechanisms. Such mechanisms yield a non trivial internal structure to resonant interfaces which proves to be useful in the selection of discontinuous solutions. It turns that discontinuous solutions may not be unique for thin interfaces. The augmented formulation (2.17) actually allows for another regularization mechanism based on thick interfaces. The color function which is naturally discontinuous (for the description of thin interfaces) is regularized in the thick regime. Such a regularization technique has been analyzed in one space variable, and existence and uniqueness of a solution for the Cauchy problem has been established. In the next section, we show how to extend this regularization procedure to several space variables.

Remark 2.1. A canonical example of coupling functions satisfying the above conditions is

C 0 (u, v) = (1 -v)γ -(u) + vγ + (u), C i (u, v) = (1 -v)a - i (γ -(u)) + va + i (γ + (u)), 1 ≤ i ≤ d.

A framework for multi-component coupling problems.

Multi-component coupling of initial boundary value problems. We are in a position to present the general coupling framework we intend to analyze in this paper. The proposed extension treats the coupling of (L + 1), L ≥ 1, distinct conservation laws in several space dimensions, with possible covering. The coupling modeling via augmented PDEs relies on a partition of the space R d in a finite number of non-overlapping, non-empty and open sets (D l ) 0≤l≤L :

(2.18)

L l=0 D l = R d .
The set of boundaries B are given by

(2.19) B = k =l D k ∩ D l .
An interface H kl is by definition the part of the boundary of D k which is only shared with D l (see also Fig. 1 for an example with N = 2 and L = 3):

(2.20)

H kl = (D k ∩ D l ) \ i =k,l D i .
These interfaces H kl are supposed to be smooth enough so that they admit an unit normal vector ν kl (y), which is well-defined except at some "exceptional" points (like corners, etc.). We suppose the set of boundaries B to be of d-dimensional Lebesgue measure zero, and, more precisely, the remaining set B \ (∪ k =l H kl ) has only components of Hausdorff dimension less than or equal to (d -2) (see for example the four points underlined in Figure 1). In each domain D l , the unknown w is governed by a specific conservation law with flux-function

D 0 D 0 D 1 D 2 H 02
A l = (a l i ) 1≤i≤d : w ∈ R → A l (w) ∈ R d : (2.21) ∂ t w + d i=1 ∂ xi a l i (w) = 0, w(t, x) ∈ R, t > 0, x ∈ D l .
Following the description introduced in the previous section, we start focusing the discussion on the definition of constant states (2.4)-(2.5)-(2.6) for the global problem set on the whole space R d . These solutions are recovered through a certain change of variable in each subdomain D l , for l = 0, . . . , L,

(2.22) u(t, x) = θ l (w(t, x)), t > 0, x ∈ D l ,
so that the stationnary solutions w(x) for the coupled problem (2.21) are the real constants u in the u variable:

(2.23) u(w(x)) = u , x ∈ R d \ B.
The coupling functions θ l are supposed to map increasingly R onto itself and we denote once again γ l the inverse functions:

(2.24)

γ l = θ -1 l , l = 0, . . . , L.
Observe that a different outlook where the coupling functions would be associated to the interfaces H kl rather than to the domains themselves could only be local in space and therefore would not allow a matching of local constant solutions so as to define a global constant solution. Here we take advantage of the local formulation at each interface in terms of the traces of w, say w(t, y k ) and w(t, y l ) on the D kand D l -side of H kl , respectively (relatively to its normal ν kl (y)):

θ k (w(t, y k )) = θ l (w(t, y l )), t > 0, y ∈ H kl .
The following augmented PDE formulation is based on a vector-valued color function that merges the description of the coupling problem. In this multi-domain approach, this function is based on the set of characteristic functions of each domain:

(2.25) v ∅ = χ D0 , v 1 = χ D1 , . . . , v L = χ D L ,
so that the change of variable (2.22) may be also rewritten

(2.26) u(t, x) = θ l (w(t, x)), x ∈ R d \ B such that v l (x) = 1.
Observe that since the (L+1) domains are a partition of the whole space R d , only L of the above characteristic functions are useful to complete the coupling description of the (L + 1) domains. Up to some relabeling we choose v 1 , . . . , v L , so that v ∅ is recovered thanks to

(2.27) v ∅ (x) = 1 - L l=1 v l (x), x ∈ R d \ B.
Multi-component coupling based on an augmented PDE's system. In the following we make use of the vector-valued color function v = (v 1 , . . . , v L ). At this stage, it takes values in the discrete set {0} ∪ {e 1 } ∪ . . . ∪ {e L } where e l stands for the l-th canonical vector of R L . This color function is intended to be regularized and to take values in the convex hull

B L + = v = (v 1 , . . . , v L ) ∈ R L v l ≥ 0, L l=1 v l ≤ 1 . The problem (2.21) is then understood in the augmented form (with t > 0, x ∈ R d ) (2.28) ∂ u C 0 (u, v)∂ t u + d i=1 ∂ u C i (u, v)∂ xi u = 0, ∂ t v = 0,
where the coupling functions C 0 and C i are assumed to restore the formulation (2.21) in terms of u in each open set D l , that is:

(2.29) lim v→0 C 0 (u, v) = γ 0 (u), lim v→e l C 0 (u, v) = γ l (u), lim v→0 C i (u, v) = a 0 i (γ 0 (u)), lim v→e l C i (u, v) = a l i (γ l (u)), 1 ≤ i ≤ d.
The following smoothness and monotonicity assumptions are required

(2.30) C 0 , C i ∈ C 2 (R × B L + ), (2.31) ∂ u C 0 (u, v) > 0, u ∈ R, v ∈ B L + .
This last property ensures the validity of the change of variable u → C 0 (u, v) for any fixed v, and the non-degenerate nature of the time-arrow in the augmented equations (2.28).

In this context, the augmented system in the main unknown u reads (2.32)

∂ t C 0 (u, v) + d i=1 ∂ xi C i (u, v) - d i=1 L l=1 ∂ v l C i (u, v)∂ xi v l = 0, ∂ t v = 0.
In the following, it will be useful to consider the same system written in the variable w = C 0 (u, v) (denoted by w(u, v), and with inverse u(w, v) for each fixed v). Equipped with such a change of unknown, (2.32) becomes (2.33)

∂ t w + d i=1 ∂ xi f i (w, v) - d i=1 L l=1 l i (w, v)∂ xi v l = 0, ∂ t v = 0,
where

f i (w, v) = C i (u(w, v), v) and l i (w, v) = ∂ v l C i|u (u(w, v), v
) with i ∈ {1, . . . , d} and l ∈ {1, . . . , L} (i.e. = ∇ v C). Hereafter and to shorten the notation, we write (2.34)

∂ t w + ∇ • f (w, v) -(w, v) : ∇v = 0, ∂ t v = 0,
with obvious notation.

Entropy stability and well-posedness. As already emphasized, in this work we propose a regularization mechanism based on thick interfaces that are modeled by any suitable regularized version of the discontinuous vector-valued color function v introduced in (2.25)-(2.27). For definiteness, we shall consider color functions

v in W 2,∞ (R + × R d , B L + )
. Obviously, it suffices to choose the initial data v 0 in W 2,∞ (R d , B L + ) so as to inherit from the required smoothness in the v solution of the augmented equations (2.34). In turn and arguing about this smoothness property, the equations under consideration reduce to an inhomogeneous scalar equation in w:

(2.35) ∂ t w + ∇ • f (w, v(x)) = (w, v(x)) : ∇v(x),
where the right-hand side just plays the role of a classical source term; namely this term does not contribute to the definition of the possible discontinuities of w. At a point of jump, (2.35) just resumes to the classical Rankine-Hugoniot condition

(2.36) -σ(w + -w -) + d i=1 f i (w + , v) -f i (w -, v) = 0.
A selection criterion of the admissible weak solutions w is of course needed, and we recast the balance law (2.35) in the main variable u:

(2.37)

∂ t C 0 (u, v) + d i=1 ∂ u C i (u, v)∂ xi u = 0
for all smooth solutions. For such solutions, additional equations are deduced and based on any (strictly) convex function → U( ), by multiplying (2.37) by U (C 0 (u, v)),

(2.38) ∂ t U(C 0 (u, v)) + d i=1 ∂ u Q i (u, v)∂ xi u = 0, where (2.39) Q i (u, v) = u U (C 0 (θ, v))∂ θ C i (θ, v)dθ, 1 ≤ i ≤ d.
We thus get from (2.38) the equivalent form for smooth solutions u:

(2.40)

∂ t U(C 0 (u, v)) + d i=1 ∂ xi Q i (u, v) = d i=1 L l=1 ∂ v l Q i (u, v)∂ xi v l .
Observe that the above right-hand side is nothing but a classical source term since we again emphasize that the color function v is smooth. As a consequence, the weak form of (2.40) for discontinuous solutions u reads:

(2.41)

∂ t U(C 0 (u, v)) + d i=1 ∂ xi Q i (u, v) ≤ d i=1 L l=1 ∂ v l Q i (u, v)∂ xi v l ,
which naturally plays the role of an (inhomogenous) entropy inequality for selecting the relevant weak solutions. Hereafter, we shall make use of the inequalities (2.41) for all convex entropy U. These will be alternatively invoked (essentially when the color function is locally constant) in the w variable:

(2.42)

∂ t U(w) + d i=1 ∂ xi F i (w, v) - d i=1 L l=1 L l i (w, v)∂ xi v l ≤ 0, with (2.43) F i (w, v) = Q i (u(w, v), v), L i (w, v) = ∂ v Q i| u (u(w, v), v), 1 ≤ i ≤ d.
To shorten the notation, equation (2.42) are written as

(2.44) ∂ t U(w) + ∇ • F(w, v) -L(w, v) : ∇v ≤ 0.
The inhomogeneous scalar conservation law (2.35) supplemented with all the entropy inequalities (2.42) naturally falls within Kruzkov's theory of entropy solutions, since the color function v belongs to W 2,∞ (R d , B L + ). Therefore, Kruzkov's uniqueness theorem for scalar conservation law with smooth inhomogeneities applies and asserts the uniqueness of the entropy weak solution of the Cauchy problem (2.35)-(2.42) with initial data

w 0 ∈ L 1 (R d ) ∩ L ∞ (R d ).
Hereafter, we shall prove existence and uniqueness of a solution to the coupled problem (2.35)-(2.42) thanks to a multidimensional well-balanced finite volume method formulated on general triangulations. Here, the well-balanced property means that the solutions in the u variable is kept constant in time and space as soon as the initial data u 0 is chosen constant whatever the definition of the (smoothly varying in space) color function v. This well-balanced property is obviously a constancy property of primary importance.

3.

A well-balanced finite volume scheme for coupling problems 3.1. Terminology and assumptions. Before stating our main result, we introduce some notation and motivate our formulation of the finite volume method under consideration. To meet the well-balancing property, the finite volume framework we develop uses two families of triangulations. The first triangulation, denoted by T h , is made of general polyhedra and will be referred to as the primal mesh. Then a closely related triangulation is of concern, the dual mesh T h , whose polyhedra are derived from the edges of the primal one. As we shall see, dual meshes may not uniquely defined from T h and it will turn that a given choice essentially affects the closed-form of expression of the CFL restriction in the (time explicit) finite volume method.

Equipped with these primal and dual meshes, approximate solutions u h and v h of the Cauchy problem (2.32) with initial data (u 0 , v 0 ), are sought as piecewise constant functions. In constrast with the usual approach, constant values for u h and v h will not be co-localized: u h (and v h , respectively) will assume constant values in each polyhedron of the primal mesh (and the dual mesh, resp.).

To facilitate the derivation of the proposed well-balanced scheme, we shall take advantage of the regularity of the color function v, which provides some room for the specific definition of the discrete approximation v h : it may range from a local averaged form to a point-wise evaluation. Without real loss of generality, we use an average value of v along each edge of the primal mesh. This choice allows to bypass the definition of the dual mesh from the edges of the primal one: a convex sequence of reals, in turn, provide sufficient information on the dual mesh. On the ground of this observation, we shall give a first brief but sustained mathematical presentation of the finite volume method under consideration. We shall then be in a position to state the main result of this paper. At last, we shall close this section with a comprehensive construction of the proposed finite volume approximation when deriving dual meshes from the primal one.

The primal mesh, T h , is a general (locally finite) triangulation of R d made of non-overlapping, non-empty, and open polyhedra :

∪ K∈T h K = R d .
We assume that for every pair of distinct polyhedra K, K ∈ T h the set K ∩ K is either an edge e of both K and K or a set with Haussdorf dimension less than or equal to d -2. The set of edges of a polyhedron K is denoted by ∂K; and for each e ∈ ∂K, ν K,e ∈ R d represents the outward unit normal vector to the edge e (see Figure 2). The volume of K and the (d -1)-measure of e are denoted |K| and |e|, respectively. Given an edge e in K, K e denotes the unique polyhedron in T h that shares the same edge e with K. We set h = sup K∈T h h K , where h K is the exterior perimeter of the polyhedron K, and assume that the triangulation T h satisfies the following non degeneracy condition (3.1) sup

K h K p K |K| ≤ C
for some constant C > 0. Here, p K denotes the perimeter of K defined by p K = e∈∂K |e|.

It is unnecessary, at this stage, to provide a comprehensive derivation of the dual mesh T h that one could define from the edges e in the primal mesh T h . Recall that, by design, a dual mesh is made of non-overlapping, non-empty, and open polyhedra denoted by K (e) with ∪ e∈T h K (e) = R d . By construction, both sets K (e)∩K and K (e)∩K e are non-empty for all pair (K, K e ) of adjacent polyhedra parametrized by the edges e in T h . Note that the set K (e) ∩ K is a subcell of K. Then, the only information about T h that is required in this section is a given convex sequence of reals prescribed in each polyhedron K in T h ; we denote by {α K,e } {e,e∈∂K} , that satisfies (for any

K in T h ) (3.2) 0 < α K,e < 1 (e ∈ ∂K), e∈∂K α K,e = 1.
We will see later that the coefficient α K,e is nothing but the ratio of the volume of K (e) ∩ K to the volume of K, where K (e) stands for the dual polyhedron of K attached to any edge e in ∂K:

(3.3) α K,e = |K (e) ∩ K| |K| , e ∈ ∂K.
At last, the time increment, denoted by τ , is assumed to satisfy τ h ≤ C and the primal mesh to be constrained by

(3.4) C 1 ≤ |e| h ≤ C 2
for some constants C, C 1 , C 2 > 0. Whereas the latter is probably not an optimal condition, it sufficies to ensure the non degeneracy of the mesh: all one-dimensional characteristic lengths are of order h. A key property for the forthcoming CFL condition, is that under these assumptions the area |K (e) ∩ K| is not smaller than O(h 2 ): there exists a positive constant c such that

(3.5) ch 2 ≤ |K (e) ∩ K|.
We use the notation t n = nτ . As already underlined, we will seek at each time level t n approximate solutions u h and v h of the Cauchy problem (2.32) with initial data (u 0 , v 0 ), under the form of piecewise constant functions with:

(3.6) u h (x, t n ) = u n K , x ∈ K, K ∈ T h , v h (x, t n ) = v h (x) = v e , x ∈ K (e), e ∈ T h .
Here and since the solution v in the Cauchy problem (2.32) does not depend on time, it seems natural to set v h (x, t n ) = v(x) = v h 0 (x) ∈ R L for all time level t n , for some discrete approximation v h 0 of the smooth function v 0 . We introduce

(3.7) v h (x) = v e = 1 |e| e v 0 (y)dy, x ∈ K (e), e ∈ T h ,
while the discrete version of the possibly discontinuous initial data u 0 is chosen according to the usual full averaging procedure over each polyhedron K:

(3.8) u 0 h (x) = u 0 K = 1 |K| K u 0 (y)dy, x ∈ K, K ∈ T h .
Remark 3.1. Since v 0 is regular, any other consistent definition for the constant value v e in K (e) would have been relevant. The interest in the particular choice (3.7) stems from the following Green formula, valid for each polygonal domain K:

X e∈∂K v el ν K,e |e| = K ∇ • (v l (x) X)dx = X K ∇v l (x)dx,
where X denotes any fixed vector in R d and v el (and v l , respectively) the l-th component of the vector v e ∈ R L (and v, resp.). Hence the proposed average value in (3.7) comes with the identity: K ∇v l (x)dx = e∈∂K v el ν K,e |e|. In a tensorial notation, we thus get K ∇v(x)dx = e∈∂K v e ⊗ ν K,e |e|.

The evolution in time of the discrete solution u h will rely on a family of numerical flux-functions, associated with each edge e of any polyhedron K in T h . Besides other properties, these numerical flux functions must meet some consistency property with the exact equation for governing u in (2.34), namely:

(3.9) ∂ t w(u, v)+∇•f (w(u, v), v)-(w(u, v), v) : ∇v = 0, x ∈ K, t ∈ (t n , t n+1 ).
Observe that in the neighborhood K (e) of each edge e, where v h reduces to a constant value v e , the above equation boils down to the scalar equation in the unknown w = w(u, v e ):

(3.10)

∂ t w + ∇ • f (w, v e ) = 0, x ∈ K (e) ∩ K, t ∈ (t n , t n+1 ).
This in turn leads us to define the required numerical flux function at each edge e in T h as a locally Lipschitz continuous two-point flux-function g e,K (., .; v e ) : R×R → R that satisfies the consistency property:

(3.11) g e,K (w, w; v e ) = f (w, v e ) • ν K,e ,
the conservation property:

(3.12) g e,K (w, w e ; v e ) = -g e,Ke (w e , w; v e ), for all reals w and w e , and the monotonicity property (3.13) ∂g(w, w e ; v e ) ∂w ≥ 0, ∂g(w, w e ; v e ) ∂w e ≤ 0.

In addition, we assume that the numerical flux depend (locally) Lipschitz continuously in the variable v e . Standard 3-point monotone schemes in the scalar framework obey (3.11)-(3.13) and that the main results in this paper are easily extended to all E-schemes (Osher [START_REF] Osher | Riemann solvers, the entropy condition, and difference approximations[END_REF]). For clarity, the dependence in the parameter v e appears explicitly in the numerical flux-function g e,K (., .; v e ). Remark 3.2. Since the function g(., .; .) : R 3 → R is locally Lipschitz continuous in its three arguments, for all compact K ⊂ R 3 , there exists some positive constant C K such that for all triple (w (1) , w

(1) e , v (1)
e ) and (w (2) , w

(2) e , v (2)
e ) in K, the following estimate holds true: |g e,K (w (2) , w (2) e ; v (2) e )g e,K (w (1) , w (1) e ; v (1) e )| ≤ C K |w (2) w (1) | + |w (2) ew (1) e | + |v (2) ev (1) e | . 3.2. Definition of the well-balanced scheme. We are now in a position to define the finite volume approximation of (3.9). Assuming that the approximate solution u h (., t n ) is known at time t n , we determine the evolution up to the next time level t n+1 as follows:

• Subcell reconstruction. At each time t n in each polyhedron K of T h , we consider any edge e ∈ ∂K and introduce the subcell state

(3.14) w n K,e = C 0 (u n K , v e
), e ∈ ∂K, as well the following average over all edges of K (3.15)

w n K = e∈∂K α K,e w n K,e .
• Evolution in time. In order to the discrete solution u h at time t n+1 , we define (in each polyhedron K) u n+1 K to be the unique solution of (3.16)

e∈∂K α K,e C 0 (u n+1 K , v e ) = w n+1 K ,
where the state w n+1 K is given by the finite volume scheme (3.17)

w n+1 K = w n K - τ |K| e∈∂K g e,K (w n K,e , w n Ke,e ; v e )|e| + τ |K| e∈∂K f (w n K,e , v e ) • ν K,e |e|.
This completes the description of our numerical method. The proposed finite volume method is explicit in time and, for the sake of stability, we need to impose a CFL (Courant, Friedrichs, Lewy) condition which reads, for all polyhedra K in T h and edges e ∈ ∂K,

|K| α K,e |e| = |K (e) ∩ K| |e| ≥ c C 2 h,
so that the CFL condition can not imply the degeneracy of the time step τ , that decreases at most as O(h). We will see in Section 4 how to build suitable primal and dual meshes. Several comments are in order. First observe that the constitutive assumptions (2.30)-(2.31) on the coupling function C 0 (., .) immediately yields existence and uniqueness of a solution to the nonlinear equation (3.14) so that the finite volume method (3.14)-(3.17) is well defined. The formulas (3.14) and (3.16) obviously express the same identity at the times t n and t n+1 , and are redundant: the finite volume method essentially reduces to (3.16)-(3.17). As they stand, they nevertheless ease the description of the method.

Next, it is worth observing that the consistency condition (3.11) allows in (3.17) to recast the flux balance e∈∂K f (w n K,e , v e )ν K,e |e| as e∈∂K g e,K (w n K,e , w n K,e ; v e )|e|. Here we stress that at each edge e in ∂K, both the numerical flux-function g e,K (w n K,e , w n Ke,e ; v e ) and its counterpart f (w n K,e , v e ) • ν K,e are evaluated thanks to the subcell values w n K,e (3.14) and not to their averaged form w n K in (3.15). The motivation is twofold. In a first hand, the two flux balances involved in (3.17), namely e∈∂K g e,K |e| and e∈∂K f (w n K,e , v e ) • ν K,e |e|, make the proposed formula to be a consistent finite volume approximation of the exact equation (3.9) for governing u: namely, the first one will be seen hereafter to be consistent with ∇ • f (w, v) while the second one actually provides a consistent approximation of the source term (w, v) : ∇v. In a second hand, the discretization of the source term is seen to be well-balanced. Proposition 3.3 (Well-balanced property). When the initial data u 0 for (3.9) is a constant function u 0 (x) = u (x ∈ R d ), then, for any choice of the color function v in (3.9), the discrete solution u h of (3.14)-(3.17) is also constant, with

(3.19) u h (x, t n ) = u 0 (x) = u , x ∈ R d
for all time level t n .

In other words, the finite volume method (3.14)-(3.17) is well-balanced with respect to all the natural equilibria of (3.9).

Proof. The discrete initial data (3.8) clearly reads u 0 h (x) = u for all x in R d so that at the first subcell reconstruction step, we get w 0 K,e = C 0 (u , v e ) = w 0 Ke,e for any edge e of an arbitrary polyhedron K in T h . Consequently, the numerical flux g e,K (w 0 K,e , w 0 Ke,e ; v e ) at any edge e boils down to f (w 0 K,e , v e ) • ν K,e in view of the consistency condition (3.11). Namely the two flux balances in the updating formula (3.17) cancel out and we end up with w 1 K = w 0 K = e∈∂K α K,e C 0 (u , v e ) thanks to the definition (3.15). Arguing about uniqueness, we thus get when solving (3.16) u 1 K = u for any polyhedron K of T h : namely u h (x, t 1 ) = u for all x in R d . An immediate recursion extends the result to the subsequent time levels.

To conclude this paragraph, it is worth illustrating that the last flux-balance entering the finite volume approximation (3.17) actually provides a consistent approximation of the source term (w, v) : ∇v. For the sake of simplicity, we temporarily adopt (cf. Remark 2.1):

C 0 (u, v) = (1 -v)γ -(u) + vγ + (u), C i (u, v) = (1 -v)a -(γ -(u)) + va + (γ + (u)), 1 ≤ i ≤ d, so that f (w, v) and (w, v) in (3.9) read f (w(u, v), v) = (1 -v)A -(γ -(u)) + vA + (γ + (u)), and (w(u, v), v) = A + (γ + (u)) -A -(γ -(u)) . It can be then readily computed: e∈∂K f (w n K,e , v e ) • ν K,e |e| = e∈∂K (1 -v e )A -(γ -(u n K )) + v e A + (γ + (u n K )) • ν K,e |e| = A + (γ + (u n K )) -A -(γ -(u n K )) • e∈∂K v e |e|ν K,e + A -(γ -(u n K )) • e∈∂K |e|ν K,e
which is nothing else a consistent discretization of (w(u, v)) : ∇v, in view of the representation formula in Remark 3.1 (for ∇v) and the identity e∈∂K |e|ν K,e = 0. These straightforward calculations allows to bridge the finite volume formula (3.17) to the governing equation (3.9) for w(u, v), expressed over K, namely where v h does achieve distinct values. The gap in between (3.9) and its reduced version (3.10) (i.e. with x ∈ K (e) ∩ K) will be definitely closed when revisiting the finite volume approximation (3.14)-(3.17) with primal-dual meshes (in Section 4).

Main convergence result.

We are now in a position to state the main result of this paper. The rest of this paper is devoted to a proof of this theorem. Then w n+1 K in (3.17) are recovered by the following averaging procedure:

(4.2) w n+1 K = e∈∂K α K,e w n+1,- K,e
.

Observe that the finite volume formula (4.1) for w n+1,- K,e is nothing but a consistent approximation of the one dimensional conservation law: ∂ t w + ∇ • f (w, v e ) = 0. The reason for calling w n+1,- K,e a subcell state will be explained in this paragraph and is at the core of the re-interpretation of the finite volume formula (3.17) with primal-dual meshes.

To further proceed, let us underline that the identity (4.2) just expresses that w n+1 K actually is a convex decomposition of the subcell states w n+1,- K,e . When understood in their quasi-one dimensional form (4.1), the latter can be recognized as extensions to the present inhomogenous setting of partial states entering similar convex decompositions that have proved well suited in the analysis of homogeneous multidimensional finite volume methods [START_REF] Coquel | Convergence of finite difference schemes for conservation laws in several space dimensions[END_REF][START_REF] Coquel | Convergence of finite difference schemes for conservation laws in several space dimensions: a general theory[END_REF]. Indeed, the interest in such a convex decomposition primary stems from the fact that many of the basic stability properties satisfied by the scheme (4.1) in one space variable are right away inherited in several space variables thanks to convexity under some CFL restriction. Observe that the relevant CFL condition for (4.1) reads .

(4.3) τ |K| |e| α K,e
In other words, all the steps involved in the method are locally conservative: this natural property will play a central role in the forthcoming analysis.

4.2.

A reformulation of the scheme. The derivation of a dual mesh T h from the edge of the primal one T h may be performed as follow. For any (open) polyhedron K, the idea is to pick an internal node x K in K which choice is left arbitrary at this stage. Such a procedure is given below a systematic definition independent of the mesh refinement h. Equipped with the node x K , we define for any edge e in K the convex hull of e and x K . The interior of this convex hull, we denote by E(x K , e), yields a non-empty open polyhedron made of (d + 1) edges. Observe that the following properties are met by construction: for any pair of distinct edges e, e in ∂K with K an arbitrary polyhedron in T h

(4.5) E(x K , e) ∩ K = E(x K , e), E(x K , e) ∩ E(x K , e ) = ∅,
while e∈∂K E(x K , e) = K modulo a null set. Then, the required definition of the polyhedron K (e) of the dual mesh T h , attached to any edge e in T h with adjacent polyhedron K and K e , follows from

(4.6) K (e) = E(x K , e) ∪ E(x Ke , e).
We refer the reader to Figure 2 for an illustration. The constructive procedure for defining the internal node x K independently of h relies on the set of vertices ϑ of the polyhedron K, together with a convex sequence of reals {β K,ϑ } {ϑ,ϑ∈K} satisfying:

0 < β K,ϑ < 1, ϑ ∈ K; ϑ∈K β K,ϑ = 1.
The required internal node x K in K is then defined by its coordinates in R d : x K = ϑ∈K β K,ϑ x ϑ , where x ϑ stands for the coordinates of the vertex ϑ. This construction ensures the correct behavior of the primal and dual meshes with the definition of the α K,e and with the previous non-degeneracy assumptions (3.3)-(3.5), the CFL condition (3.18) is then only modified according to the choice of the function v and its discrete representation.

To further proceed in the comprehensive derivation of the finite volume framework, some additional notation is in order. For any K in T h and e in ∂K, an edge of a dual polyhedron K (e) ∈ T h or of the subcell K (e) ∩ K of K will be indifferently denoted by e . Observe that with little abuse in the notation, an edge e of some cell K of the primal mesh T h is also a dual edge of the subcell K (e) ∩ K: see indeed Figure 2. At last ν K (e),e ∈ R d stands for the outward unit vector normal to the edge e .

Equipped with these notation, we are in a position to re-interpret the quasione dimensional state w n+1,- K,e introduced in (4.1) in term of a state in the subcell K (e) ∩ K of K, thanks to the following simple but key identity: We can clarify the origin of the definition g e ,K (e) = f (w n K,e , v e ) • ν K (e),e for edges e distinct from e. For such an edge e , it is worth introducing the adjacent subcell K (e )∩K to K (e)∩K in K: i.e. with e in ∂K such that K (e )∩K (e) = e . Note that e is of course distinct from e . We then successively rewrite the left-and right-hand numerical flux at e , say g e ,K (e) (respectively g e ,K (e ) ), as follows:

f (w(u n K , v e ), v e ) • ν K (e),e , respectively : -f (w(u n K , v e ), v e ) • ν K (e),e
, since by definition (3.14) w n K,e = w(u n K , v e ) and w K,e = w(u n K , v e ) and, equivalently, (4.10)

f (w(u, v), v) • ν K (e),e (ω(0-)), resp.: -f (w(u, v), v) • ν K (e),e (ω(0+)),
where ω(0 ∓ ) stands for the left and right traces at ξ = 0 of the self-similar function

ω : ξ ∈ R ξ → (u(ξ), v(ξ)) ∈ R × R L given by (4.11) ω(ξ) = (u n K , v e ), ξ < 0, (u n K , v e ), ξ > 0.
From Section 2, recall that the Riemann solution of (4.12)

∂ t w + ∂ x f (w(u, v), v) • ν K (e),e (u, v) -∂ v f (w(u, v), v) • ν K (e
),e : ∇v = 0,

∂ t v = 0
(with initial data ((u n K , v e ), x < 0, (u n K , v e ), x > 0)) consists in a standing wave separating (u n K , v e ) from (u n K , v e ), and thus coincides with ω(ξ) in (4.11). It is therefore clear that the flux-functions in (4.10) actually results from the Godunov method applied to the augmented system (4.12) at the edge e . In other terms, the finite volume formula (4.8)-(4.9) in each subcell K (e) ∩ K may be understood as an approximation of the balance law for governing u in (3.9):

∂ t w(u, v) + ∇ • f (w(u, v), v) -(w(u, v), v) : ∇v = 0, x ∈ K, t ∈ (t n , t n+1 ).
This interpretation closes the gap in between the governing equation (3.9) for u and its reduced form (3.10) expressed in w:

∂ t w + ∇ • f (w, v e ) = 0, x ∈ K (e) ∩ K, t ∈ (t n , t n+1 ).
4.3. Sup-norm estimates. Throughout the upcoming sections, the assumptions of Theorem 3.4 are tacitly assumed to be valid. Their formulations are thus skipped over in any forthcoming statements. The main result of this section ensures that the sequence of approximate solutions u h stays uniformly bounded in L ∞ (R + × R d ) as a consequence of the following result. Proposition 4.2 (Maximum principle). The finite volume method satisfies the following inequalities (in the variable u):

(4.13) min u n K , min e∈∂K u n Ke ≤ u n+1 K ≤ max u n K , max e∈∂K u n
Ke in each polyhedron K in T h and at all time level t n .

Since v 0 ∈ W 2,∞ immediately implies a sup-norm estimate for v h given by (3.7), we easily deduce, from the maximum principle (4.13), an additional uniform supnorm estimate but for w h = C 0 (u h , v h ) arguing about the regularity properties (2.30) of C 0 :

(4.14) ||w h || L ∞ (R + ×R d ) ≤ O(1).
Besides the monotonicity assumption (3.13) met by the numerical flux functions, we stress that the preservation of conservativity (4.4) in the subcell reconstruction procedure plays a central role in the validity of the reported maximum principle, as highlighted in the proof. The latter will be carried out using a recursion procedure based on subsequent partitions of the set of edges e in K. To fix the notation and up to some relabeling, {e 1 , . . . , e J K } represents the full set of edges e ∈ ∂K so that here the index J K is given by #{e, e ∈ ∂K}. Subsets of the form {e 1 , ..., e J }, with increasing index J ∈ {1, ..., J K }, will be of concern as follows. Being given J with 1 ≤ J ≤ K, let us attach to the subset {e 1 , ..., e J } the solution u n+1- K,{e1,...,e J } of the following nonlinear equation:

(4.15) 1≤j≤J α K,ej C 0 (u n+1- K,{e1,...,e J } , v ej ) = 1≤j≤J α K,ej w n+1- K,ej ,
where the subcell states w n+1- K,ej are defined in (4.1), Lemma 4.1. Again, the constitutive assumptions (2.30)-(2.31) ensure existence and uniqueness of a solution to (4.15).

Arguing about the conservation property (4.4) satisfied at the subcell reconstruction step, it is worth observing that u n+1- K,{e1,...,e J K } can be identified with the final state u n+1 K at time t n+1 in the finite volume approximation (3.14)-(3.17). Therefore, the recursion under consideration naturally ends up as soon as the index J To condense the notation, we introduce the two functions

u → Ψ J (u) = 1≤j≤J α K,ej C 0 (u, v ej ) and u → ψ (J+1) (u) = α K,e (J+1) C 0 (u, v e (J+1) )
so as to deduce:

Ψ J (u n+1- K,{e1,...,e (J+1) } ) -Ψ J (u n+1- K,{e1,...,e J } ) • ψ (J+1) (u n+1- K,{e1,...,e (J+1) } ) -ψ (J+1) (u n+1- K,e (J+1)
) ≤ 0. since by assumption (3.2) α K,e (J+1) > 0. But the monotonicity hypothesis (2.31) on C 0 together with again assumption (3.2) imply that both functions u → Ψ J (u) and u → ψ (J+1) (u) strictly increase with u so that the above inequality yields min(u n+1- K,{e1,...,e J } , u n+1- K,e (J+1) ) ≤ u n+1- K,{e1,...,e (J+1) } ≤ max(u n+1- K,{e1,...,e J } , u n+1- K,e (J+1) ). Lemma 4.3 implies

min(u n K , u n Ke (J+1) ) ≤ u n+1- K,e (J+1) ≤ max(u n K , u n Ke (J+1)
), and the proof is completed.

Entropy inequalities

5.1. Preliminaries. Proposition 4.2 asserts sup-norm boundedness for the sequence u h which in the absence of an a priori strong compactness argument, leads us to study the structure of the Young measure µ associated with {u h } h>0 . Recall that such a Young measure represents all the composite weak-star limits a(u h ) of u h with continuous functions a ∈ C 0 (R), namely for all continuous functions in a single variable

a(u h ) - < µ, a >= R a(λ)dµ(λ),
weakly-star in L ∞ . We propose to establish that the measure µ under consideration reduces to a Dirac measure, and hence to prove the strong convergence of u h , invoking DiPerna's uniqueness theorem [START_REF] Diperna | Measure-valued solutions to conservation laws[END_REF] for entropy measure-valued solutions.

In this section we derive the required discrete entropy inequalities together with the a priori estimates that are needed to handle the passage to the limit in the sense of measure valued solutions. In this respect, the main issue is to assess the relevance of the Young measure µ in such a limit. Indeed, discrete entropy inequalities generically involve numerical flux functions, that are continuous functions but of (at least) two arguments: the sequence u h (.) itself and its shift ∆ h u h = u h (. + h). Nonlinear superposition of possible discrete oscillations in u h and its shift ∆ h u h may prevent the usual Young measure µ to represent the composite weak-star limit of G(u h , ∆ h u h ). Counterexamples have been constructed in Coquel and LeFloch [START_REF] Coquel | Convergence of finite difference schemes for conservation laws in several space dimensions: a general theory[END_REF]. Some weak control over possible discrete oscillations is therefore mandatory in order to justify the applicability of µ in the limiting form of discrete entropy inequalities.

The requisite weak estimate corresponds to some estimate of the discrete entropy dissipation rate in the finite volume approximation. The derivation of several specific estimates with distinctive features have been the matter of a large literature following Coquel and LeFloch [START_REF] Coquel | Convergence of finite difference schemes for conservation laws in several space dimensions[END_REF]. (The reader is referred to the introduction where several subsequent contributions were quoted.) The estimates we derive now generalize the ones in Cockburn, Coquel, and LeFloch [START_REF] Cockburn | Convergence of finite volume methods for multidimensional conservation laws[END_REF]. The entropy dissipation estimate does not allow actually to pass to the weak limit in arbitrary numerical entropy-flux functions, but nevertheless turns out to be sufficient in order to handle discrete entropy inequalities. The main interest in such an estimate stems from the simplicity of its derivation. 5.2. Discrete entropy estimates. We first focus on the derivation of the discrete entropy inequalities and then the required weak estimate. The passage to the limit in the discrete inequalities is the subject of the following section. After Crandall and Majda [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF], assumptions (3.11)-(3.13) on the numerical flux functions g e,K are known to yield a full set of discrete entropy inequalities for scalar conservation laws.

Here and in the light of Section 2, the scalar conservation laws of concern have to be found locally at each edge e in T h , and take the generic form (5.1)

∂ t w + ∇ • f (w, v) = 0,
for a given v ∈ R. Associated entropy pairs were defined earlier in (2.39)-(2.43).

The inequalities stated below are naturally built from the subcell states w n+1,- K,e (4.1) of Lemma 4.1 and in this regard may be understood as subcell entropy inequalities. We refer the reader to [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF] for a proof of this classical result. As already claimed, the weak estimate will not allow to pass weakly to the limit in arbitrary numerical entropy flux-functions. We thus propose to merge inequalities (5.4) in such a way that solely exact entropy flux-functions F(w, v e ) • ν K,e enter the weak form. t n e φ(x, t)dxdt.

Then, the following discrete weak inequality holds The proof is postponed to the end of this section. We shall easily deduce from the discrete inequality (5.6) the following continuous (in space) inequality.

K∈T h e∈∂K α K,e U(w n+1,- K,e ) -U(w n K,e ) φ n e |K| - ]t n ,t n+1 [×R d Q(u n h , v(x)) • ∇φ(x, t) + φ(x, t)∂ v Q(u n h , v(x)) : ∇v(x)dxdt ≤ O(h)τ ||φ|| W 1,∞ (]t n ,t n+1 [×R d ) |supp(φ)|.
The proof, at the end of this section, essentially makes use of the uniform supnorm estimate (4.13) for the sequence u h together with the regularity assumption

v 0 ∈ W 2,∞ .
Clearly, the Young measure µ can tackle the weak limit of the space derivatives involved in inequality (5.7) extended to any time interval (0, T ), T > 0. Such a claim then naturally rises the question of passing to the weak limit in the discrete time derivative. The latter is conveniently decomposed as (5.8)

K∈T h e∈∂K α K,e U(w n+1,- K,e ) -U(w n K,e ) φ n e |K| = K∈T h U(w n+1 K ) -U(w n K ) φ n K |K| - K∈T h e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e
) φ n e |K| -K∈T h e∈∂K α K,e U(w n K,e ) -U(w n K ) φ n e |K|, where (5.9)

φ n K = e∈∂K α K,e φ n e .
The last two error terms entering the righ-hand side of (5.8) are devoted to sum up (5.10)

n≥0 K∈T h e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e ) φ n e |K| + n≥0 K∈T h e∈∂K α K,e U(w n K,e ) -U(w n K ) φ n e |K|,
with other error terms in the right-hand side of the discrete entropy inequalities (5.7). The former must therefore be proved to go to zero with h.

∂ v Q(u n K , v K ) : e∈∂K (v e -v K ) ⊗ ν K,e |e| = ∂ v Q(u n K , v K ) : 1 τ t n+1 t n K ∇v(x)dtdx ,
again thanks to the representation formula in Remark 3.1 (for ∇v). On the other hand, the latter error term in (5.17) is described by

(5.19) e∈∂K (φ n e -φ n K )∂ v Q(u n K , v K ) : ((v e -v K ) ⊗ ν K,e )|e| ≤ O(1) sup e∈∂K |(φ n e -φ n K )(v e -v K )|)p K ≤ O(h 2 K ) ∇φ L ∞ (]t n ,t n+1 [×K) p K ≤ O(h K ) ||∇φ|| L ∞ (]t n ,t n+1 [×K) |K|.
Here, we have successively used the sup-norm estimate (4.13) satisfied by u h , the definition of the perimeter p K of K, the estimate

∂ v Q(u n K , v K ) : e∈∂K φ n e (v e -v K ) ⊗ ν K,e |e| - 1 τ t n+1 t n K φ n K ∂ v Q(u n K , v K ) : ∇v(x)dtdx ≤ O(h K ) ||∇φ|| L ∞ (]t n ,t n+1 [×K) |K|.
For the final error term in the flux balance (5.16), we have the following bounds:

e∈∂K φ n e 1 0 (∂ v Q(u n K , v K + s(v e -v K )) -∂ v Q(u n K , v K ))ds : ((v e -v K ) ⊗ ν K,e )|e| ≤ O(1) sup e∈∂K |v e -v K | 2 p K ||φ|| L ∞ (]t n ,t n+1 [×K) ≤ O(h K )||φ|| L ∞ (]t n ,t n+1 [×K) |K|,
where we have used the regularity of the entropy flux Q, the sup-norm estimate (4.13), the estimate (5.20) satisfied by |v e -v K |, and the non-degeneracy assumption (3.1) on the triangulation T h .

To summarize, we have obtained the estimate for the flux balance on a single cell:

(5.22)

1 τ t n+1 t n K Q(u n K , v K ) • ∇φ(x, t) + φ n K ∂ v Q(u n K , v K ) : ∇v(x)dtdx - e∈∂K F(w n K,e , v e ) • ν K,e φ n e |e| ≤ O(h)||φ|| W 1,∞ (]t n ,t n+1 [×K) |K|.
We now derive the companion estimate (5.12), by starting from the decomposition

e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e ) φ n e = φ n K U(w n+1 K ) - e∈∂K α K,e U(w n+1,- K,e ) + e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e ) φ n e -φ n K ,
and observing, on one hand,

e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e ) φ n e -φ n K ≤ O(1) e∈∂K α K,e |φ n e -φ n K ||w n+1,- K,e -w n+1 K | ≤ O(h K ) e∈∂K α K,e |w n+1,- K,e -w n+1 K | ∇φ L ∞ (]t n ,t n+1 [×K)
and, on the other hand,

e∈∂K α K,e U(w n+1,- K,e ) -U(w n+1 K ) = U (w n+1 K ) e∈∂K α K,e w n+1,- K,e -w n+1 K + e∈∂K α K,e 1 0 U (w n+1,- K,e + s(w n+1 K -w n+1,- K,e ))ds(w n+1,- K,e -w n+1 K ) 2 .
Finally, in view of the convex decomposition (4.2) stating

w n+1 K = e∈∂K α K,e w n+1,- K,e , we get
U(w n+1 K ) - e∈∂K α K,e U(w n+1,- K,e ) ≤ -σ U e∈∂K α K,e |w n+1,- K,e -w n+1 K | 2 ,
where σ U denotes the convexity like-modulus of U introduced in Lemma 5.4. This concludes the proof.

5.3.

Entropy dissipation rate and strong convergence. The proposed estimates obtained in Lemma 5.4 deserve a few comments. Plugging first estimate (5.11) in (5.10) will be easily seen to yield the following upper-bound

n≥0 K∈T h e∈∂K α K,e U(w n K,e ) -U(w n K ) φ n e |K| ≤ O(h) φ W 1,∞ (R+×R d ) |supp(φ)|
that obviously suffices to conclude. By contrast and turning considering (5.12), a crude upper-bound based on the sup-norm estimate (4.14), say

e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e ) φ n e ≤ O(h) φ W 1,∞ (]t n ,t n+1 [×K)
would result in the useless estimate

n≥0 K∈T h e∈∂K α K,e U(w n+1 K )-U(w n+1,- K,e ) φ n e |K| ≤ O(1) φ W 1,∞ (R+×R d ) |supp(φ)|.
Proving that the error term of concern in (5.10) actually vanishes with h requires therefore in turn a sharper control in (5.12) of the oscillations of the w n+1,- K,e around their mean value w n+1 K . Such a control over these discrete oscillations results from a sharp evaluation of the discrete entropy rate of dissipation. (5.26)

N T n=0 K∈T h e∈∂K α K,e |w n+1 K -w n+1,- K,e | 2 ψ K |K| ≤ O(1),
where ψ K reads ψ K = e∈∂K α K,e ψ e , ψ e = 1 |e| e ψ(x)dx. Equipped with (5.26) we obtain the following entropy dissipation rate.

Corollary 5.6. The sequence u h satisfy the entropy like inequality (5.27)

R + ×R d U(C 0 (u h , v))∂ t φ(x, t) + Q(u h , v) • ∇φ + φ∂ v Q(u h , v) : ∇vdxdt ≥ O(h 1/2 ),
for any (smooth) convex entropy pair (U, Q) : R → R × R d introduced in (2.41) and (2.39). Equipped with the above inequality valid for any entropy pair (U, Q), we easily deduce that the Young measure µ = µ t,x associated with the sequence (u h ) h>0 is an entropy satisfying measure valued solution. In other words the uniformly bounded L ∞ sequence (u h ) h>0 , as announced at the beginning of this section, it is easy to check that the inequation (5.27) becomes as h tends to 0 the following inequation satisfied in the weak sense:

(5.28)

∂ t µ, U(C 0 (•, v)) + ∇ x µ, Q(•, v) -µ, ∂ v Q(•, v) : ∇v ≤ 0.
Relying on a direct extension of DiPerna's uniqueness theorem [START_REF] Boutin | Coupling techniques nonlinear hyperbolic equations. III. A regularization method based on thick interfaces[END_REF], we can deduce that the entropy measure-valued solution µ t,x reduces to a Dirac measure δ u(t,x) concentrated on a function u = u(t, x) since the initial data µ 0 coincides with the Dirac measure δ u0 (where u 0 is the initial data in the Cauchy problem (2.32)).

Proving that the inital data u 0 is correctly handled amounts to show that for every compact subset K of R we have (5.29) lim

t→0+ t 0 K µ s,x , |id -u 0 (x)| dxds = 0.
The condition (5.29)-(5.28) reduces to a Dirac measure concentrated at u(t, x), the Kruzkov entropy solution of (2.32)-(2.41) with same initial data u 0 . In other words, for all time T > 0 and for all compact K in R, the scheme converges strongly in L p loc ((0, T )×K) to the solution u. Theorem 3.4 of this paper is thus now established.

Proof of Proposition 5.5. We start from the discrete in time weak formulation (5.7) stated in Proposition 5.3:

K∈T h e∈∂K α K,e U(w n+1,- K,e ) -U(w n K,e ) φ n e |K| - ]t n ,t n+1 [×R d Q(u n h , v(x)) • ∇φ(x, t) + φ(x, t)∂ v Q(u n h , v(x)) : ∇v(x)dxdt ≤ O(h)τ ||φ|| W 1,∞ (]t n ,t n+1 [×R d ) |supp(φ)|,
in which we plug the decomposition (5.8)-(5.9). A discrete test function ψ K given for any given time-independent test function ψ ∈ D(R d ) is considered. We then get Then routine arguments give the conclusion from (5.30). where w ± (t, x) = lim θ→0+ w(t, x±θν x ) and ν x the exterior unit normal at x ∈ ∂D 1 .

The initial data plotted in Figure 3a is piecewise constant: w 0 (x, y) = 1, x < -0.8, 0, x ≥ -0.8.

The computations are performed on a Cartesian grid with 100 × 100 meshes, and the CFL number is chosen to be 0.5.

In an homogeneous domain with flux function f 0 , such an initial data would develop a shock front moving with the speed vector 0.5(1, 1) T . In the present heterogeneous domain, this shock front has the same behavior until it reaches the interface between both domains (see Figures 4a). The coupling condition at this interface is such that the value w = 2 arises then inside the domain D + . In this second domain, where the flux under consideration is f 1 , we observe then a (curved) shock wave connecting the states w = 2 and w = 0 and moving at the fixed speed given by the Rankine-Hugoniot relation, that is, 0.605(1, 1) T (see Figures 4b and4c). Finally, the shock front goes outside the whole domain [-1, 1] 2 (see Figure 4d). In Figures 4a, 4b, 4c, and4d) (right figures), we plot the u-variable, which is found to remain constant at each interface, as expected by the theory.

Figure 1 .

 1 Figure 1. Boundaries (in bold-face H 02 , circle points being excluded)

  sup u∈[m,M ] ∂f (w(u, v e ), v e ) ∂w ≤ 1, where m = inf x∈R d u 0 (x) and M = sup x∈R d u 0 (x). Due to the dimensional hypothesis (3.2)-(3.4)-(3.5) the ratio |K|α K,e /|e| satisfies

Theorem 3 . 4 .

 34 (Well-balanced finite volume method for multi-dimensional coupling problems). Consider the Cauchy problem (2.32)-(2.41) with initial data u 0 ∈ L ∞ (R d ) and v 0 ∈ W 2,∞ (R d) under the constitutive assumptions (2.30)-(2.31). Let u h be the sequence of approximate solutions defined by the finite volume method (3.7)-(3.8) and (3.14)-(3.17) with numerical flux-functions satisfying the conditions (3.11)-(3.13). Then under the CFL restriction (3.18), the sequence u h is uniformly bounded in L ∞ (R + × R d ) and converges (when h → 0) in the L p loc norm strongly (1 ≤ p < ∞) to the unique entropy solution u to the problem (2.32)-(2.41): namely for all time T > 0 and for all compact K in R d lim h→0 ||uu h || L p ((0,T )×K) = 0.

g

  e,K (w n K,e , w n Ke,e ; v e )g e,K (w n K,e , w n K,e ; v e ) w n Ke,ew n K,e ≤ 1, and hence the CFL restriction (3.18). At last and arguing about the definition (4.1), the subcell reconstruction step (3.15) at time t n+1 and the formula (3.16), we deduce the (seemingly trivial) identities (4.4) e∈∂K α K,e w n+1 K,e = w n+1 K = e∈∂K α K,e w n+1,- K,e

Figure 2 .

 2 Figure 2. Primal and dual meshes, edges and vertices.

e

  ∈K (e)∩K |e |ν K (e),e = 0, i.e. |e|ν K,e =e ∈K (e)∩K, e =e |e |ν K (e),e .It is then straightforward to recast w n+1,- e |K| g e,K (w n K,e , w n Ke,e ; v e )|e| + τ α K,e |K| f (w n Ke ) • ν K,e |e|, = w n K,e -τ |K (e) ∩ K| g e,K (w n K,e , w n Ke,e ; v e )|e| + e ∈K (e)∩K e =e f (w n K,e , v e ) • ν K (e),e |e | , where we have used the interpretation (3.3) of α K,e . Introducing the numerical flux formula: (4.8) g e ,K (e) = g e,K (w n K,e , w n Ke,e ; v e ), if e = e; f (w n K,e , v e ) • ν K (e),e ,

  K,e -τ |K (e) ∩ K| e ∈K (e)∩K g e ,K (e) |e |.

Lemma 5 . 1 ( 5

 515 Entropy inequalities per cell). Let (U, F) : R → R×R d be any convex entropy pair for the scalar conservation law(5.1), where e denotes any edge in ∂K for an arbitrarily K in T h . Then there exists a numerical entropy flux function G e,K : R 2 → R that satisfies the consistency property (5.2) G e,K (w, w; v e ) = F(w, v e ) • ν K,e , the conservation property (5.3) G e,K (w, w e ; v e ) = -G e,Ke (w e , w; v e )for all reals w and w e , so that the following discrete entropy inequality holds (G e,K (w n K,e , w n Ke,e ; v e ) -F(w n K,e , v e ) • ν K,e ≤ 0.

Lemma 5 . 2 .

 52 Let φ be any non-negative test function in D(R * + × R d ). Define for any edge e in T h , the average

  , v e ) • ν K,e φ n e |e| ≤ 0.

Proposition 5 . 3 .

 53 The finite volume approximation (3.14)-(3.17) obeys at each time level t n the following (discrete in time) entropy inequality(5.7)

Lemma 5 . 4 .

 54 For any polyhedron K in T h , one has(5.11) e∈∂K α K,e U(w n K,e ) -U(w n K ) φ n e ≤ O(h 2 ) φ W 1,∞ (]t n ,t n+1 [×K) ,Proof of Proposition 5.3. We begin with the discrete inequality (5.6) of Lemma 5.2 and specifically considerb the flux balance K∈T h e∈∂K F(w n K,e , v e ) • ν K,e φ n e |e|. Our purpose is to shift the mathematical expressions under consideration from the w to the u variable. Hence let us write F(w n K,e , v e ) = F(w(u n K , v e ), v e ) = Q(u n K , v e ) with Q(u, v) the exact entropy flux introduced in (2.39), which we repeat we get(5.18)

(5. 20 )

 20 |v ev K | ≤ e ∈∂K α K,e |v ev e | ≤ O(h K )from the definition of v K and the regularity propertyv 0 ∈ W 2,∞ , a similar estimate |φ n eφ n k | ≤ O(h K )and finally the non degeneracy assumption (3.1) on the triangulation T h . Involving (5.18)-(5.19), the identity (5.17) yields the following estimate(5.21)

Proposition 5 . 5 .

 55 Let T > 0 be any fixed time and let N T ∈ N be the floor of T /τ we denote [T /τ ]. Then, for any (time independent) non negative test function ψ ∈ D(R d ), the finite volume approximation (3.14)-(3.17) obeys the following estimate on the discrete oscillations:

  n+1 K ) -U(w n K ) ψ K |K| -]t n ,t n+1 [×R d Q(u n h , v(x)) • ∇ψ(x) + ψ(x)∂ v Q(u n h , v(x)) : ∇v(x)dxdt ≤ K∈T h e∈∂K α K,e U(w n+1 K ) -U(w n+1,- K,e) ψ e |K|+ K∈T h e∈∂K α K,e U(w n K,e ) -U(w n K )) ψ e |K| + O(h)τ ||ψ|| W 1,∞ (R d ) |supp(ψ)|.Invoquing estimates (5.11)-(5.12) then yieldsK∈T h U(w n+1 K ) -U(w n+1,- K,e ) ψ K |K| + σ U K∈T h e∈∂K α K,e |w n+1 K w n+1,- K,e | 2 ψ K |K| ≤ O(h)τ ||ψ|| W 1,∞ (R d ) |supp(ψ)| + O(h) K∈T h ||∇ψ|| L ∞ (K) |K| + O(h 2 ) K∈T h ||ψ|| W 1,∞ (K) |K| + ]t n ,t n+1 [×R d Q(u n h , v(x)) • ∇ψ(x) + ψ(x)∂ v Q(u n h , v(x)) : ∇v(x)dxdt.Observe that due to the estimate (4.13), the last contribution in the above righthand side can be given the following crude estimate O(τ) ψ W 1,∞ (R d ). Henceforth, we deduce thatK∈T h U(w n+1 K ) -U(w n K ) ψ K |K| + σ U K∈T h e∈∂K α K,e |w n+1 K w n+1,- K,e | 2 ψ K |K| ≤ O(h) ψ W 1,∞ (R d ) .Summing over time indices n ∈ [0, N T ] with N T = [T /τ ] for a fixed time T > 0, we getR d U(w h (x, T ))ψ h (x)dx + σ U 0 (x))ψ h (x)dx + O(1)T ψ W 1,∞ (R d ) ,which is the required result.Proof of Corollary 5.6. We start from (5.7)-(5.8)-(5.9) and consider the following discrete in time weak formulation for the time dependent test function φ ∈ D(R

| 2 χ φ |K|τ 1 / 2

 12 as a consequence of the convexity property of the α K,e -average. The estimate (5.26) then yields withψ = χ φ n≥0 K∈T h e∈∂K α K,e |w n+1 K w n+1,- K,e | χ φ ||∇φ|| L ∞ (]t n ,t n+1 [×K) |K|τ ≤ O(h 1/2 ).

6. Numerical experiments 6 . 1 . 1 .

 611 A two domain coupling problem. In this first test, we consider an heterogeneous medium which occupies the spatial domain [-1, 1] 2 and is constituted by a background domain D 0 and a ring-shaped inclusion D 1 centered at the origin (0, 0) with external radius √ 0.2 and with internal radius √ 0.1. In these two domains, the following respective flux functions are considered in term of the scalar unknown w = w(t, x):The regularized color function v plotted in Figure3bprovides us with a regularized version of the characteristic function of the domain D 1 . The coupling condition between D 0 and D 1 takes here the form 2 w -(t, x) = w + (t, x), x ∈ ∂D 1 ,

( a )Figure 3 .

 a3 Figure 3. Initial data for the multidimensional test.

( a )

 a Solution at t = 0.5 (b) Solution at t = 1.5 (c) Solution at t = 2.5 (d) Solution at t = 4.5

Figure 4 .

 4 Figure 4. Two domain test case. Evolution of the solution: w (left) and u (right).

6. 2 . 1 ,Figure 5 .

 215 Figure 5. Geometry of the three domains.

Figure 6 .

 6 Figure 6. Three domain test case. Solution w at different times.

4 .

 4 Finite volume approximations on primal-dual meshes 4.1. A convex combination. One of our objectives in this section is explaining how the coefficients α K,e should be determined. Arguing about the formuladefinitions (3.14)-(3.15) at time t n and the consistency condition (3.11), we obtain the following statement.Lemma 4.1 (Edge values and convex combination). For any polyhedron K of T h and edge e in ∂K, let us define the following subcell states:

	(4.1) w n+1,-K,e	= w n K,e -	|e| α K,e	τ |K|	g e,K (w n K,e , w n Ke,e ; v

e )g e,K (w n K,e , w n K,e ; v e ) .

  • ∇φ(x, t) + φ(x, t)∂ v Q(u n h , v(x)) : ∇v(x)dxdt ≤ O(h)τ ||φ|| W 1,∞ (]t n ,t n+1 [×R d ) |supp(φ)| + O(h 2 ) K∈T h ||φ|| W 1,∞ (]t n ,t n+1 [×K) |K| | ||∇φ|| L ∞ (]t n ,t n+1 [×K) |K|. | χ φ ||∇φ|| L ∞ (]t n ,t n+1 [×K) |K|τ ,making use of the characteristic function χ φ of 0<t<T supp(φ(•, t)), a compact subset of R d , where T is a finite time such that supp(φ(•, t)) = ∅ for t ≥ T . Cauchy-Schwarz's inequality then yields the following crude upper bound for the last term:| χ φ ||∇φ|| L ∞ (]t n ,t n+1 [×K) |K|τ

	R d ) and its discrete representation φ n K
	K∈T h	U(w n+1 K ) -U(w n K ) φ n K |K|
	-h , v(x)) + O(h) ]t n ,t n+1 [×R d Q(u n K∈T h e∈∂K α K,e |w n+1 K	-w n+1,-K,e
	where we have used estimates (5.11)-(5.12). Summing this inequality over time
	indices gives
	(5.30)
	-	n≥0 K∈T h	U(w n+1 K )	φ n+1 K τ -φ n K	τ |K|
	-		R + ×R d	Q(u n
				n≥0 K∈T h e∈∂K	α K,e |w n+1 K	-w n+1,-K,e
			n≥0 K∈T h e∈∂K	α K,e |w n+1 K	-w n+1,-K,e
			≤	n≥0 K∈T h e∈∂K	α K,e |w n+1 K	-w n+1,-K,e	| χ φ	2 |K|τ	1/2
				•	n≥0 K∈T h	||∇φ|| 2 L 1/2
			≤ O(1)
								+ * ×

h , v(x)) • ∇φ(x, t) + φ(x, t)∂ v Q(u n h , v(x)) : ∇v(x)dxdt ≤ O(h)||φ|| W 1,∞ (R + ×R d ) + O(1) ∞ (]t n ,t n+1 [×K) |K|τ

reaches the value J K . In order to initiate the recursion and propagate it, we need the following statement concerned with the values u n+1- K,{e J } , 1 ≤ J ≤ J K , solutions of C 0 (u n+1- K,{e J } , v e J ) = w n+1- K,e J .

Lemma 4.3 (Local maximum principle). The maximum principle holds true at any edge e J in ∂K:

), 1 ≤ J ≤ J K .

Then the maximum principle "propagates" to sets {e 1 , ..., e J }, as follows.

Lemma 4.4. The solution u n+1- K,{e1,...,e J } to (4.15) with J ∈ {1, ..., J K }, obeys the following maximum principle:

) .

The proposed lower and upper bounds for u n+1- K,{e1,...,e J K } , i.e. the estimate in the lemma with J = J K , just reads the expected local maximum principle (4.13) for u n+1 K , since again u n+1 K coincides with u n+1- K,{e1,...,e J K } by construction.

Proof of Lemma 4.3. To alleviate the notation we skip the index J and first point out an estimate valid under the CFL restriction (3.18) for any edge e in ∂K:

≤ max(w n K,e , w n Ke,e ) as a well-known consequence of the monotonicity assumptions (3.13) satisfied by the numerical flux function g e,K (., .; v e ). We then recall that the subcell reconstruction step (3.14) builds w n K,e = C 0 (u n K , v e ) while the identity w n+1- K,e = C 0 (u n+1- K,{e} , v e ) holds from our definition. We can thus recast (4.16) as:

, from which we immediately deduce the desired estimate, namely

), e ∈ ∂K since the function C 0 is by assumption (2.31) strictly increasing in its first argument.

Proof of Lemma 4.4. The desired lower-upper bounds with J = 1 are stated in Lemma 4.3. Then, assuming the validity of the maximum principle at rank J, 1 ≤ J < J K , this one is proved to hold at the rank (J + 1) starting from (4.15):

, v e (J+1) ).

We recast the above identity as follows:

, v e (J+1) ) .

while

(5.12)

where σ U denotes some convexity-like modulus of U: U (u) ≥ σ U > 0, for all u ∈ (m, M ) where the bounds m, M were introduced in (3.18) in agreement with the maximum principle (4.13).

Proof of Lemma 5.2. Let e be any edge in T h and K, K e the associated pair of adjacent polyhedra. Multiplying the subcell entropy inequality (5.4) valid for K by α K,e |K| and the companion inequality for K e by α K,e |K e |, we get

We then recast the flux balance as

where the average of the states v e is defined by

The treatment of the last remaining discrete term relies on the following identity:

which leads us to rewrite (5.15):

(5.16)

The matrix (v ev K ) ⊗ ν K,e |e| with size L × d appears as a discrete representation for the continuous function ∇v. The first term in the above right-hand side is rewritten as:

(5.17)

where the discrete flux function φ n K is obtained by averaging: φ n K = e∈∂K α K,e φ n e . On one hand, owing to the identity e∈∂K (v ev K ) ⊗ ν K,e |e| = e∈∂K v e ⊗ ν K,e |e| From the discrete weak entropy inequality (5.6) we recall that ) -U(w n K,e ) φ n e |K| -

Proof of Lemma 5.4. We first establish the estimate (5.11) and consider the following decomposition involving again the {α K,e } {e,e∈∂K} -average φ n K of the φ n e (5.9):

from which we deduce the following bound:

(5.23)

in view of the sup-norm estimate (4.13) satisfied by u h , the estimate |φ n eφ n K | ≤ O(h K ) and the convexity of the entropy U(w). The first error term in (5.23) is given the following bound: