Xavier Blanc
email: blanc@ann.jussieu.fr

Emmanuel Labourasse

A positive scheme for diffusion problems on deformed meshes

Keywords: finite volume, diffusion, maximum principle, numerical analysis MSC

2010) 65M08,65M12,04A25

, and uses auxiliary unknowns at the nodes of the mesh. The flux is computed so as to be a two-point nonlinear flux, giving rise to a matrix which is the transpose of an M-matrix, which ensures that the scheme is positive. A particular attention is given to the computation of the auxiliary unknowns. We propose a new strategy, which aims at providing a scheme easy to implement in a parallel domain decomposition setting. An analysis of the scheme is provided: existence of a solution for the nonlinear system is proved, and the convergence of a fixed-point strategy is studied.

Introduction

We consider the unsteady diffusion equation

∂ t u -div (κ∇u) = f,
in Ω × (0, T ), γ (κ∇u) • n + δu = g on ∂Ω × (0, T ),

where Ω is a bounded open domain of R 2 , n is the outgoing normal to Ω. The data are f ∈ L 2 (Ω), g ∈ H 1/2 (∂Ω), and κ ∈ L ∞ (Ω) satisfies the ellipticity condition ∀x ∈ Ω, κ(x) ≥ κ 0 > 0.

The functions γ and δ are smooth functions such that ∀x ∈ ∂Ω, δ(x) ≥ 0, and γ(x) ≥ γ 0 > 0.

The system (1) is completed by an initial data

u(t = 0) = u 0 ∈ H 1 (Ω). (2)
For the sake of simplicity, we restrict our attention to the case of an isotropic diffusion, that is, κ(x) ∈ R. However, a simple adaptation allows to treat any bounded matrix κ(x) satisfying ellipticity conditions. We refer to [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF] for the details. Under the above conditions, the system (1)-( 2) has a unique solution in C 0 [0, T ], H 1 (Ω) ∩ C 1 [0, T ], L 2 (Ω) . See for instance [START_REF] Evans | Application of nonlinear semigroup theory to certain partial differential equations[END_REF] for the proof. Further, we have a maximum principle: if f ≥ 0, g ≥ 0 and u 0 ≥ 0, then u ≥ 0 in Ω×[0, T ]. The equation being linear, it also implies that, if m ≤ f ≤ M, mδ ≤ g ≤ M δ and m ≤ u 0 ≤ M , then m ≤ u ≤ M.

In the case when system (1)-( 2) is coupled with hydrodynamics, and if the hydrodynamic equations are solved using a Lagrangian or ALE method, the mesh may be highly distorted. Further, the unknown u is naturally discretized using a piece-wise constant approximation on this mesh. This is why we are interested in finite-volume approximations of (1)-( 2) on deformed meshes.

The scheme needs to be consistent, stable, and reproduce on the discrete level qualitative properties of the system (1)- [START_REF] Arbogast | Mixed finite elements for elliptic problems with tensor coefficients as cellcentered finite differences[END_REF], such as conservativity and maximum principle. Moreover, in the case when the unknown u is a temperature or a concentration, it is mandatory that u remains positive. Finally, it is desirable that the scheme produces a symmetric matrix, and is linear. However, as we will see below, these last two conditions will be dropped to the advantage of the previous ones. Further, we need the stencil of the scheme to be small, in order to make the parallelization of the scheme amenable. In summary, the ideal scheme should satisfy the following: 1. consistency; 2. stability; 3. conservativity; 4. small stencil (in link with parallelization); 5. maximum principle; 6. symmetry; 7. linearity.

The scheme we present here satisfies all these properties, except 5, 6 and 7. A weaker version of item 5 is proved, that is, the scheme is positive (see below). It should be noted that, in practice, we did not find any example in which the maximum principle is violated.

Many works have been devoted to the present subject. Let us mention those we know of, which might not be exhaustive (see the review article [START_REF] Droniou | Finite volume schemes for diffusion equations: Introduction to and review of modern methods[END_REF] for more details):

-To our knowledge, the first attempt to derive a consistent finite volume scheme for diffusion equations on a deformed mesh was the work of Kershaw [START_REF] Kershaw | Differencing of the diffusion equation in Lagrangian hydrodynamic codes[END_REF] (see also [START_REF] Pert | Physical constraints in numerical calculations of diffusion[END_REF] for a related scheme). The basic idea is to consider a transformation from a reference quadrilateral mesh to the actual deformed mesh, and write down a standard finite volume scheme in this reference configuration, then transform it in order to obtain a finite volume scheme on the deformed mesh. This scheme gives a symmetric matrix, but was proved to be consistent only on meshes consisting of parallelograms, and does not satisfy the maximum principle.

-Another approach is to use a diamond scheme, which was analyzed in [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem. M2AN[END_REF]. The idea is to introduce auxiliary unknowns at the nodes of the mesh. These additional unknowns are computed using an interpolation method as functions of the cell unknowns. The node unknowns are then used to compute a second-order approximation of the fluxes. This method converges at order two if the interpolation method is sufficiently precise. However, this scheme is not positive, and the associated matrix is in general not symmetric.

-The method of mixed finite element [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] may be recast in a finite volume formulation. Here, we use additional degrees of freedom at the edges of the mesh. A hybrid formulation of (1) then allows to eliminate the cell unknowns, so that the system we need to solve gives the edge unknowns, and the cell unknowns are then computed accordingly. This scheme is convergent of order two, the corresponding matrix is symmetric, but it is not positive. A numerical analysis of such schemes was given in [START_REF] Arbogast | Mixed finite elements for elliptic problems with tensor coefficients as cellcentered finite differences[END_REF].

-More recently, the discrete duality finite volume method (DDFV) was proposed by F. Hermeline [START_REF] Hermeline | A finite volume method for second-order elliptic equations. (Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF][START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF][START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF][START_REF] Hermeline | Approximation of 2D and 3D diffusion operators with variable full tensor coefficients on arbitrary meshes[END_REF][START_REF] Hermeline | Nouvelles méthodes de volumes finis pour approcher des équations aux dérivées partielles sur des maillages quelconques[END_REF][START_REF] Hermeline | A finite volume method for approximating 3D diffusion operators on general meshes[END_REF]. In this method, one also uses additional unknowns localized at the nodes of the mesh. These unknowns are not related to the cell unknowns by an interpolation procedure, but are computed by writing down, on a dual mesh, a diffusion scheme approximating [START_REF] Aavatsmark | Convergence of a symmetric MPFA method on quadrilateral grids[END_REF]. Doing so, one solves two diffusion problems instead of one. This scheme is convergent of order two, even if the mesh contains non-convex cells, and the associated matrices are symmetric, but not positive.

-The mimetic finite difference method is described in the papers [START_REF] Brezzi | Mimetic finite differences for elliptic problems[END_REF][START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF][START_REF] Kuznetsov | The mimetic finite difference method on polygonal meshes for diffusion-type problems[END_REF][START_REF] Lipnikov | The mimetic finite difference discretization of diffusion problem on unstructured polyhedral meshes[END_REF] (see also the review paper [START_REF] Lipnikov | Mimetic finite difference method[END_REF]). In this method, the fluxes are considered as additional unknowns. Moreover, the discrete system is designed so that some properties of the continuous system are reproduced, such as, for instance, the Green formula. The scheme is convergent of order two, and the corresponding matrix is symmetric. However, it does not satisfy the maximum principle (see [START_REF] Lipnikov | Monotonicity conditions in the mimetic finite difference method[END_REF]), and the number of degrees of freedom is much larger than in other schemes.

-The scheme using stabilization and harmonic interfaces (SUSHI) was proposed by R. Eymard, T. Gallouët and R. Herbin [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: A scheme using stabilization and hybrid interfaces[END_REF]. It is based on the same ideas as the standard diamond scheme. However, a stabilization term is added to the discrete gradient, which preserves consistency, while improving robustness. This scheme is convergent of order two, the corresponding matrix being symmetric, but the scheme is not positive. It is proved in [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF] that, actually, SUSHI and mimetic finite difference methods are part of the same family of methods.

-The multi-point flux approximation (MPFA) [START_REF] Aavatsmark | Convergence of a symmetric MPFA method on quadrilateral grids[END_REF][START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF][START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] uses additional unknowns at the edges of the mesh, with possibly several of them on each face. These additional unknowns are used to compute a consistent approximation of the flux, and are then eliminated by imposing that the flux is continuous across each edge. These scheme are convergent, except on random meshes, and they give rise to non-symmetric matrices. They do not satisfy the maximum principle (see [START_REF] Michael | A quasi-positive family of continuous Darcy-flux finite-volume schemes with full pressure support[END_REF][START_REF] Eigestad | Symmetry and M -matrix issues for the O-method on an unstructured grid[END_REF][START_REF] Helmer | A family of MPFA finite-volume schemes with full pressure support for the general tensor pressure equation on cell-centered triangular grids[END_REF]).

-In [START_REF] Le | A linear scheme satisfying a maximum principle for anisotropic diffusion operators on distorted grids. (Un schéma linéaire vérifiant le principe du maximum pour des opérateurs de diffusion très anisotropes sur des maillages déformés[END_REF], a finite difference type scheme has been proposed by C. Le Potier. This scheme is convergent of order one, and satisfies the maximum principle. However, in order to manage this, the stencil is enlarged to a larger distance. Thus, this scheme is difficult to use in a decomposition domain strategy.

-In [START_REF] Siess | A linear and accurate diffusion scheme respecting the maximum principle on distorted meshes[END_REF], V. Siess studies a scheme which is convergent of order one, and satisfies the maximum principle. In order to do so, the mesh is replaced by the Voronoï mesh associated with the cell centers of the original mesh. On such a mesh, it is possible (and simple) to write down a second order consistent approximation of the fluxes using the neighbouring cell unknowns. Therefore, this scheme is consistent and satisfies the maximum principle. However, the change of mesh induces an error at order one, so that convergence is only of order one. Further, the stencil may be larger than for standard finite volume methods, and extension to anisotropic diffusion is not obvious.

-In [START_REF] Droniou | Construction and convergence study of schemes preserving the elliptic local maximum principle[END_REF], a nonlinear scheme is designed with additional unknowns which are not at the nodes nor at the faces of the mesh. These additional unknowns are computed by interpolation. The scheme is convergent of order two, satisfies the maximum principle, and gives a non-symmetric matrix in general. Its main drawback is that it is nonlinear.

-A scheme which is very similar to the preceding one was proposed in [START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF], in which the additional unknowns are located at the edges of the mesh. The implementation of this scheme is slightly different from [START_REF] Droniou | Construction and convergence study of schemes preserving the elliptic local maximum principle[END_REF], but the properties are very similar: it converges with the same rate, satisfies the maximum principle, produces an non-symmetric matrix, and is nonlinear.

-The scheme proposed in [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF][START_REF] Wang | A monotone finite volume scheme for advection-diffusion equations on distorted meshes[END_REF][START_REF] Yuan | Monotone finite volume schemes for diffusion equations on polygonal meshes[END_REF] is based on the same idea, that is, defining a nonlinear scheme so that it is positive. However, our implementation is simpler and more natural, since it uses only two-point fluxes, in the spirit of [START_REF] Bertolazzi | A second-order maximum principle preserving finite volume method for steady convection-diffusion problems[END_REF][START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF]. This scheme is convergent of order two, and positive, that is, if f ≥ 0, g ≥ 0 and u 0 ≥ 0, then u ≥ 0. However, since the scheme is nonlinear, this does not imply the maximum principle.

-Finally, monotone nonlinear schemes based on the idea of building consistent two-point flux approximations, but without interpolation of additional unknowns, was proposed in [START_REF] Danilov | A monotone nonlinear finite volume method for diffusion equations on conformal polyhedral meshes[END_REF][START_REF] Lipnikov | Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes[END_REF][START_REF] Lipnikov | A monotone finite volume method for advection-diffusion equations on unstructured polygon meshes[END_REF][START_REF] Lipnikov | Minimal stencil finite volume scheme with the discrete maximum principle[END_REF][START_REF] Nikitin | A monotone nonlinear finite volume method for advection-diffusion equations on unstructured polyhedral meshes in 3D. Russian[END_REF]. It avoids the problem of interpolation of additional unknowns, satisfies the maximum principle, and is precise at order two. However, in the case of highly deformed meshes, it may be necessary to enlarge the stencil in order to preserve this properties.

The scheme described in [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF][START_REF] Wang | A monotone finite volume scheme for advection-diffusion equations on distorted meshes[END_REF][START_REF] Yuan | Monotone finite volume schemes for diffusion equations on polygonal meshes[END_REF] is the one we have tested. It satisfies all the points mentioned above, except items 5, 6 and 7. The maximum principle is replaced by positivity of the scheme, which is weaker since the scheme is nonlinear. Since it uses only two-point fluxes, its stencil is reduced to edge-neighbouring cells. This implies that it is very simple to use it in a decomposition domain strategy. Note that we are not able to prove that the scheme is coercive (uniformly with respect to the mesh size). However, the numerical results indicate that it is the case. Apart from studying test cases which are different from [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF], our contribution is that we give proof of the existence of a solution to the scheme. We also provide a convergence analysis of the fixed point iteration needed to deal with the nonlinearity of the scheme.

Let us end this introduction by a practical remark: in order to have a positive scheme, one needs either to enlarge the stencil, or to use a nonlinear scheme (see [START_REF] Keilegavlen | Sufficient criteria are necessary for monotone control volume methods[END_REF][START_REF] Kershaw | Differencing of the diffusion equation in Lagrangian hydrodynamic codes[END_REF]). We use here the second strategy, which leads to a nonlinear scheme even if the problem (1)-( 2) is linear. In principle, this induces a higher computational cost, since the nonlinear equation is in general solved using a fixed point strategy. However, for the applications we have in mind, the original problem is nonlinear. Therefore, the model problem (1)-( 2) is only an intermediate step to solve a nonlinear problem. In practice, the numerical tests show that it is not necessary to reach convergence of the fixed point strategy. Only one step of this algorithm already gives satisfactory results. The outer loop (for instance a Newton algorithm) ensures consistency when it converges. Hence, the additional cost related to the nonlinearity of the scheme is relatively small, as far as the mesh distortion is not too important. For highly distorted mesh, the convergence of the Newton algorithm may be relatively slow, thereby inducing a more important additional computational cost. However, when applied to a linear problem, the present scheme is much more costly than linear schemes.

The article is organized as follows: in Section 2, we describe the method in detail, and prove that the method is positive and consistent. In Section 3, we give the implementation details. We prove in particular that, if the time step is small enough, the fixed point algorithm converges. Finally, Section 4 gives numerical tests asserting the accuracy and positivity of the scheme.

2 The method of Sheng, Yuan and Yue

Notation

We give here the notation used throughout the article: we assume that a mesh is given, and denote by -K the set of all cells of the mesh; -E the set of edges of the mesh; -N the set of nodes of the mesh; For any K ∈ K, we still denote by K the center of this cell. For any L ∈ K sharing an edge with K, we denote by e = K|L ∈ E this common edge.

Given an edge e of a cell K, we will denote by n K,e the outer normal to K on e. Finally, we define the global mesh size by ∆x = max {|e|, e ∈ E} .

(3)

Discrete fluxes

In order to write down a finite volume scheme for (1)-( 2), we integrate it over a cell K:

K ∂u ∂t dx - K ∇ • (κ∇u)dx = K f dx,
and use the divergence formula:

K ∂u ∂t dx - ∂K κ∇u • n K dΓ = K f dx,
where n K is the outer unit normal to cell K, and ∂K is the boundary of K. We define an edge e as the intersection of two neighbouring cells K and L. We note e = ∂K ∩ ∂L = K|L. We have

K ∂u ∂t dx + e∈∂K - e κ∇u • n Ke dΓ = K f dx. (4)
We set

F K,e = - e κ(x)∇u(x, t) • n Ke dΓ,
the flux going out of K through e. The point is to write down an approximation of F K,e for each edge e, as a function of the cell unknowns.

For this purpose, we introduce additional unknowns at the nodes of the mesh. In a second step, we will express these additional unknowns by an interpolation method. Let us first define the notation: (recall that K and L denote both the cells and the corresponding centers.) The points M 1 and M 2 are nodes of K such that the basis (

---→ KM 1 , ---→ KM 2
) is direct, and such that the decomposition of n K,e in this basis gives non-negative coefficients. See figure 1. The points M 3 and M 4 are defined in a similar way, with the constraint that the basis (

--→ LM 3 , --→ LM 4
) is direct, and that the decomposition of the normal vector n Le in this basis gives non-negative coefficients. We define O 1 (respectively O 2 ) the intersection between the half line starting at K (respectively L) with direction n K,e (respectively n L,e ) and the boundary of K (respectively L). We also define the angles

θ K1 = ( ---→ KM 1 , ---→ KO 1 ), θ K2 = ( ---→ KO 1 , ---→ KM 2 ), θ L1 = ( --→ LM 3 , --→ LO 2 ), θ L2 = ( --→ LO 2 , --→ LM 4 ), θ K = θ K1 + θ K2 , θ L = θ L1 + θ L2 .
Note that the nodes M i may be different from the endpoints of the edge e (see Figure 2). We then write the vector n K,e as a linear combination of ---→ KM 1 and ---→ KM 2 . We thus compute α ∈ R and β ∈ R such that 1 The cells K and L, their centers, and the points Mi for

n Ke = α ---→ KM 1 ---→ KM 1 + β ---→ KM 2 ---→ KM 2 . ( 5
) M 2 M 3 M 4 M 1 n Ke n Le O 2 O 1 L K Fig.
1 ≤ i ≤ 4. M 1 M 3 M 4 M 2 O 1 O 2 n Le n Ke K L Fig. 2
The cells K, L and the points Mi for 1 ≤ i ≤ 4 in the case when these points are not the endpoints of the edge e = K|L.

We compute the outer product of ( 5) with ---→ KM 1 , and find

n Ke ∧ ---→ KM 1 = ---→ KM 1 | sin θ K1 | n Ke ∧ ---→ KM 1 = |β| ---→ KM 1 | sin θ K | =⇒ |β| = sin θ K1 sin θ K . Since θ K ∈ (0, π] and θ K1 ∈ (0, π], we have 0 < sin θ K ≤ 1 and 0 < sin θ K1 ≤ 1 whence β = sin θ K1 sin θ K Likewise, α = sin θ K2 sin θ K . We thus have n Ke = sin θ K2 sin θ K ---→ KM 1 ---→ KM 1 + sin θ K1 sin θ K ---→ KM 2 ---→ KM 2 . (6)
A similar argument gives

n Le = sin θ L2 sin θ L --→ LM 3 --→ LM 3 + sin θ L1 sin θ L --→ LM 4 --→ LM 4 . (7)
We infer

F K,e = - e sin θ K2 sin θ K ∇u • ---→ KM 1 ---→ KM 1 + sin θ K1 sin θ K ∇u • ---→ KM 2 ---→ KM 2 κ(x)dΓ F L,e = - e sin θ L2 sin θ L ∇u • --→ LM 3 --→ LM 3 + sin θ L1 sin θ L ∇u • --→ LM 4 --→ LM 4 κ(x)dΓ
We use a finite difference approximation to compute the integrand in the above formula:

∇u • ---→ KM i ---→ KM i = u(M i ) -u(K) ---→ KM i + O(∆x).
Thus,

F K,e = -|e| sin θ K2 sin θ K u M1 -u K ---→ KM 1 + sin θ K1 sin θ K u M2 -u K ---→ KM 2 κ e + O(∆x 2 ), F L,e = -|e| sin θ L2 sin θ L u M3 -u L --→ LM 3 + sin θ L1 sin θ L u M4 -u L --→ LM 4 κ e + O(∆x 2 ),
where κ e is the value of κ(x) at the center of the edge e. We define

F 1 = -|e|κ e sin θ K2 sin θ K u M1 -u K ---→ KM 1 + sin θ K1 sin θ K u M2 -u K ---→ KM 2 , (8)
F 2 = -|e|κ e sin θ L2 sin θ L u M3 -u L --→ LM 3 + sin θ L1 sin θ L u M4 -u L --→ LM 4 . (9)
In order to use formulas ( 8) and ( 9), we need to compute the values u Mi at the nodes M i , and the value of κ on the edge e. This value is in general computed as a linear combination of κ K and κ L , as for instance

κ e = λ K,e + λ L,e λ L,e κ L + λ K,e κ K ,
where λ K,e and λ L,e are, respectively, the distance between the center of K (resp. L) and the edge e.

The unknowns u Mi are computed via an interpolation method, for which we refer to Section 3.1 below.

Formulas ( 8) and ( 9) give a consistent approximation of the flux through the edge e (up to a change of sign). We are going to combine them in order to obtain a positive scheme. For this purpose, we write

F K,e = µ 1 (u)F 1 -µ 2 (u)F 2 , (10)
F L,e = µ 2 (u)F 2 -µ 1 (u)F 1 , (11)
where the coefficients µ 1 and µ 2 are chosen below. We need these formulas to be consistent, so we impose

µ 1 + µ 2 = 1.
The idea of [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF] is to compute µ 1 and µ 2 in such a way that [START_REF] Danilov | A monotone nonlinear finite volume method for diffusion equations on conformal polyhedral meshes[END_REF] and [START_REF] Després | Lagrangian gas dynamics in 2D and lagrangian systems[END_REF] are two-point approximations of the flux. Going back to (8) and ( 9), we insert them into [START_REF] Danilov | A monotone nonlinear finite volume method for diffusion equations on conformal polyhedral meshes[END_REF] and [START_REF] Després | Lagrangian gas dynamics in 2D and lagrangian systems[END_REF]:

F K,e = µ 1 F 1 -µ 2 F 2 = -µ 1 |e|κ e sin θ K2 sin θ K u M1 -u K ---→ KM 1 + sin θ K1 sin θ K u M2 -u K ---→ KM 2 +µ 2 |e|κ e sin θ L2 sin θ L u M3 -u L --→ LM 3 + sin θ L1 sin θ L u M4 -u L --→ LM 4 .
Thus,

F K,e = µ 1 |e|κ e sin θ K1 sin θ K 1 ---→ KM 2 + sin θ K2 sin θ K 1 ---→ KM 1 u K -µ 2 |e|κ e sin θ L1 sin θ L 1 --→ LM 4 + sin θ L2 sin θ L 1 --→ LM 3 u L -µ 1 |e|κ e sin θ K1 sin θ K u M2 ---→ KM 2 + sin θ K2 sin θ K u M1 ---→ KM 1 =a1 +µ 2 |e|κ e sin θ L1 sin θ L u M4 --→ LM 4 + sin θ L2 sin θ L u M3 --→ LM 3 =a2 . ( 12
)
Note that, due to our choice or M 1 , M 2 , M 3 , M 4 , (see figures 1 and 2), we have a 1 ≥ 0 and a 2 ≥ 0 (provided u ≥ 0, which is the case, as we will see below). We are going to impose that a 1 µ 1 = a 2 µ 2 in order to have two-point fluxes. We thus need to solve the following system:

µ 1 + µ 2 = 1, a 1 µ 1 -a 2 µ 2 = 0.
As far as a 1 + a 2 = 0, this system has a unique solution, which reads

µ 1 = a 2 a 1 + a 2 , µ 2 = a 1 a 1 + a 2 .
When a 1 + a 2 = 0, in [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF][START_REF] Yuan | Monotone finite volume schemes for diffusion equations on polygonal meshes[END_REF], it is proposed to arbitrarily choose µ 1 = µ 2 = 1 2 . However, this may cause the flux to be a discontinuous function of u. This is not desirable, so we use slightly different values for a 1 and a 2 , namely

ã1 = a 1 + ∆x 2 , ã2 = a 2 + ∆x 2 ,
where ∆x is defined by (3). Then we solve the system defining µ 1 and µ 2 with these values

μ1 = a 2 + ∆x 2 a 1 + a 2 + 2∆x 2 , μ2 = a 1 + ∆x 2 a 1 + a 2 + 2∆x 2 .
Note that, if the coefficient a i are non-negative, which is the case in the interpolation procedures mentioned below (see Section 3.1), this definition of μ1 and μ2 is a continuous function of u (see the proof of Proposition 3.1 below). Finally, we re-define the fluxes as

F K,e = -F L,e = μ1 |e|κ e sin θ K1 sin θ K 1 ---→ KM 2 + sin θ K2 sin θ K 1 ---→ KM 1 u K -μ2 |e|κ e sin θ L1 sin θ L 1 --→ LM 4 + sin θ L2 sin θ L 1 --→ LM 3 u L . ( 13
)
Since the definition (12) of the fluxes is consistent with the exact fluxes of order 2 in ∆x, so is [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF].

Remark 1 In the definition of ã1 and ã2 , we have chosen to use a global mesh size ∆x, but it is also possible to use a local one, for instance

∆x LM = max {|e|, e ∈ E ∩ (K ∪ L)} .
This would allow smaller values of the regularizing terms where the mesh is fine, which numerically is relevant.

Since we have imposed F L,e = -F K,e , the scheme is conservative (see Proposition 2.4 below). Moreover, if the values u Mi are non-negative, which is the case for the interpolation procedures we use, we have μ1 > 0 and μ2 > 0, hence

F K,e = -F L,e = A K,e u K -A L,e u L ,
with A K,e ≥ 0 and A L,e ≥ 0. This implies that the matrix is the transpose of an M-matrix (see Definition 2.1 below). This in turn implies that the scheme is well-defined (and positive), as we will see below in Proposition 2.3.

S 2 S 1 K e n H
Fig. 3 The boundary cell K and its boundary edge e

Boundary conditions

The boundary conditions of ( 1) are taken into account as follows: let K be a boundary cell, with an edge e on the boundary (see figure 3). We denote by S 1 and S 2 the nodes of e, and we impose the value

-γ e F K,e + δ e u e = g e , (14)
where g e is the value of g on e, and similarly for γ e and δ e . The value u e is an approximation of the unknown u on e, yet to be determined, together with the flux F K,e through e. Now, let us define the point H, which is the orthogonal projection of K onto the line (S 1 S 2 ). Then, a consistent approximation of the flux is given by

F K,e = κ e u K -u H --→ KH , hence, u H = u K - --→ KH κ e F K,e (15)
Moreover, given the values of u at the nodes S 1 and S 2 , the following approximation of u e is of second order:

u e = 1 2 (u S1 + u S2 ) .
It is then possible to express the value u H by an affine approximation:

u e -u H = 1 2 + --→ S 2 K • ---→ S 1 S 2 ---→ S 1 S 2 2 u S1 + 1 2 - --→ S 1 K • ---→ S 1 S 2 ---→ S 1 S 2 2 u S2 . (16)
We now re-write ( 14) as -γ e F K,e + δ e u H + δ e (u e -u H ) = g e .

Inserting the value of u e -u H given by ( 16) into this equation, and then using the value of u H given by ( 15), we infer

-γ e + δ e --→ KH κ e F K,e + δ e u K = g e -δ e 1 2 + --→ S 2 K • ---→ S 1 S 2 ---→ S 1 S 2 2 u S1 + 1 2 - --→ S 1 K • ---→ S 1 S 2 ---→ S 1 S 2 2 u S2 .
Hence, we have the following value for the flux: [START_REF] Eigestad | Symmetry and M -matrix issues for the O-method on an unstructured grid[END_REF] Copyright line will be provided by the publisher The values u S1 and u S2 being computed by the interpolation method described below (Section 3.1), we use this value of the flux in the finite volume formulation.

F K,e = δ e 1 2 + --→ S2K• ---→ S1S2 ---→ S1S2 2 u S1 + 1 2 - --→ S1K• ---→ S1S2 ---→ S1S2 2 u S2 -g e + δ e u K γ e + δ e --→ KH κe
The advantages of this method (comparing with the one described in [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF]) are that it avoids the introduction of additional degrees of freedom on the boundary. Further, in the case of a Cartesian mesh, we recover the approximation by standard finite difference methods.

Note that the node values u S1 and u S2 are necessary to compute F K,e . As we will see below, in practice, we use an explicit value for these unknowns: u Si are the values computed at the preceding step of the fixed point algorithm.

A summary of the scheme

Let us summarize the scheme we have just derived: a consistent approximation of the flux from cell K to cell L through edge e = K|L is given by F K,e = A K,e u K -A L,e u L , and F L,e = -F K,e , with

A K,e = μ1 |e|κ e sin θ K1 sin θ K 1 ---→ KM 2 + sin θ K2 sin θ K 1 ---→ KM 1 , (18)
A L,e = μ2 |e|κ e sin θ L1 sin θ L 1 --→ LM 4 + sin θ L2 sin θ L 1 --→ LM 3 . ( 19
)
The coefficients μ1 and μ2 are given by

μ1 = a 2 + ∆x 2 a 1 + a 2 + 2∆x 2 , μ2 = a 1 + ∆x 2 a 1 + a 2 + 2∆x 2 , with a 1 = |e|κ e sin θ K1 sin θ K u M2 ---→ KM 2 + sin θ K2 sin θ K u M1 ---→ KM 1 , a 2 = |e|κ e sin θ L1 sin θ L u M4 --→ LM 4 + sin θ L2 sin θ L u M3 --→ LM 3
.

The node values u Mi are given by the interpolation procedure described in Section 3.1 below. The important fact is that it should satisfy u Mi ≥ 0.

Finally, we use an implicit discretization in time, so that the scheme for system (1)-( 2) is given by

   1 ∆t u n+1 -u n + M (u n+1 )u n+1 = f + g, u n+1 ≥ 0, (20)
where we have denoted by u n the vector (u n K ) K∈K . Here, K is the set of all cells of the mesh. The vector f is (f K ) K∈K , where f K = K f, so that u K is an approximation of K u. The vector g corresponds to the boundary conditions:

g K = e∈K∩∂Ω g K,e , with g K,e = g e -δ e 1 2 + --→ S2K• ---→ S1S2 ---→ S1S2 2 u n+1 S1 + 1 2 - --→ S1K• ---→ S1S2 ---→ S1S2 2 u n+1 S2 γ e + δ e --→ KH κe
The matrix M (u) is defined as follows:

[M (u)] KK = e∈K A K,e + e∈K∩∂Ω   δ e γ e + δ e --→ KH κe   , (21)
[M (u)] KL = -A L,e if K = L. (22)
Implicitly, we set A L,e = 0 if e is not an edge of L. Hence, if K and L do not share an edge, then [M (u)] KL = 0.

Properties of the scheme

First, we prove that the scheme is well-posed, that is, (20) has a solution. For this purpose, we first need some definitions: Definition 2.1 We say that a matrix M ∈ R N ×N is an M-matrix if it satisfies the following inequalities:

∀i = j, m ij ≤ 0, (23) ∀i,
N j=1 m ij ≥ 0. (24)
Moreover, if [START_REF] Hermeline | Approximation of 2D and 3D diffusion operators with variable full tensor coefficients on arbitrary meshes[END_REF] is strict for all i, we say that M is a strict M-matrix. Remark 2 Definition 2.1 is not the standard definition of an M-matrix. Indeed, in the literature, an M-matrix is a matrix M such that M = cId -A, where A has non negative entries, and c ≥ ρ(A), where ρ(A) is the spectral radius of A (see [START_REF] Plemmons | M-matrix characterizations.i -nonsingular m-matrices[END_REF] for instance). It is an easy exercise to prove that Definition 2.1 is a special case of this latter definition.

We have the following standard fact: Lemma 1 Assume that M is a strict M-matrix. Then, M satisfies the following properties:

if X ∈ R N is such that ∀i, (M X) i ≥ 0, then ∀i, X i ≥ 0. ( 25
) if X ∈ R N is such that ∀i, (M X) i > 0, then ∀i, X i > 0. ( 26
)
These properties are equivalent to the fact that M is invertible, with M -1 having non-negative coefficients.

The proof of this result can be found in many linear algebra textbooks, as for instance [START_REF] Plemmons | M-matrix characterizations.i -nonsingular m-matrices[END_REF][START_REF] Serre | Matrices[END_REF]. Definition 2.2 In the sequel, for any X ∈ R N , we will denote by X ≥ 0 the fact that all the components of X are non-negative:

X ≥ 0 ⇐⇒ ∀1 ≤ i ≤ N, X i ≥ 0.
Similarly, X > 0 means that all components of X are positive:

X > 0 ⇐⇒ ∀1 ≤ i ≤ N, X i > 0.
We are now in position to prove that the scheme is well-posed: Proposition 2.3 If f + g ≥ 0 and u n ≥ 0, then system (20) has a solution u n+1 .

Recall that, according to Definition 2.2, u ≥ 0 means that each component of the vector u are non-negative. Proof: we reproduce here the proof given in [START_REF] Droniou | Construction and convergence study of schemes preserving the elliptic local maximum principle[END_REF]. Equation ( 20) may be written φ(u n+1 ) = u n+1 , where

φ(u) = (Id + ∆tM (u)) -1 (u n + ∆t(f + g)) . (27)
Hence, we need to prove that φ has a fixed point u ≥ 0 in order to prove that (20) has a solution.

To this end, we first note that, in view of ( 21) and ( 22), M (u n+1 ) T is always an M-matrix. Hence, Id + ∆tM (u n+1 )

T is a strict M-matrix. Applying Lemma 1, we infer that the matrix Id + ∆tM (u n+1 ) -T has positive coefficients. Hence,

since u n + ∆t(f + g) ≥ 0, we have ∀u ∈ R N , φ(u) ≥ 0.
Moreover, multiplying (27) by the constant vector (1, . . . , 1) on the left, we find that

K∈K [φ(u)] K = K∈K (u n K + f K + g K ) := C 0 ,
which is a constant independent of u. We thus define the set

C = u ∈ R N , u ≥ 0, K∈K u K = C 0 . (28)
It is a convex compact subset of R N , and the application φ maps C into itself. Moreover, φ is continuous. Indeed, each coefficient of M (u) is a continuous function of u, and M → (I + M ) -1 is a continuous application from the set of Mmatrices to the set of matrices. Hence, we may apply Brouwer's theorem [START_REF] Rogers | A less strange version of Milnor's proof of Brouwer's fixed-point theorem[END_REF], which implies that φ has a fixed point in C, hence (20) has a solution.

Remark 3 Here again, the positivity property is a consequence of the structure of the matrix M , and does not depend on the fact that we use M = M u n+1 in the scheme [START_REF] Helmer | A family of MPFA finite-volume schemes with full pressure support for the general tensor pressure equation on cell-centered triangular grids[END_REF]. In particular, as far as u n ≥ 0, we still have positivity for the scheme [START_REF] Le | A linear scheme satisfying a maximum principle for anisotropic diffusion operators on distorted grids. (Un schéma linéaire vérifiant le principe du maximum pour des opérateurs de diffusion très anisotropes sur des maillages déformés[END_REF]. However, only the exact solution to the linear system is non-negative. In particular, if the linear system is solved with a poor precision, it may still exhibit negative values.

Proposition 2. [START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2d anisotropic diffusion problems[END_REF] The scheme defined by (20) is conservative.

Proof: the only thing we need to check is that if f = 0, δ = 0 and g = 0, then the scheme satisfies a discrete version of the equality d dt Ω u = 0, that is,

∀n ≥ 0, K∈K u n+1 K = K∈K u n K . (29)
For this purpose, we write (20):

∀K ∈ K, 1 ∆t u n+1 K -u n K + L∈K M u n+1 KL u n+1 L = 0,
and we sum up with respect to K. This gives

K∈K u n+1 K - K∈K u n K + ∆t K∈K L∈K M u n+1 KL u n+1 L = 0.
Using formulas ( 21) and ( 22), we find

K∈K M KL = M LL + K =L M KL = e∈L A L,e - e∈L A L,e = 0.
Remark 4 Here again, the above property is a consequence of the structure of the matrix M . In particular, M need not be exactly equal to M (u n+1 ) in order to have this property. In other words, we do not need to reach convergence of the fixed point strategy in order to have conservativity. In particular, if we replace (20) by

1 ∆t u n+1 -u n + M (u n )u n+1 = f + g, ( 30
)
the scheme is still conservative.

Implementation

Interpolation

Here we explain how we compute the node unknowns as functions of the cell unknowns, in order to compute the matrix M (u) of the scheme. For this purpose, a natural way to proceed is the one proposed in [START_REF] Sheng | A nine point scheme for the approximation of diffusion operators on distorted quadrilateral meshes[END_REF] or [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF]: take a node P of the mesh, and consider the neighbouring cells (K i ) 1≤i≤p , that is, all the cells having P as a node. Then, write down a linear approximation of u P as a function of the values u Ki :

u P = P i=1 ω i u Ki . (31)
In order to have a second-order approximation, so that the scheme may still be second-order convergent, we require formula [START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF] to be exact for any affine function. This gives the following system:

                     p i=1 ω i = 1, p i=1 ω i x Ki = x P , p i=1 ω i y Ki = y P . (32)
P Fig. 4 An example in which the point P is not in the convex envelop of the neighbouring cell centers.

In this system, the unknowns are the weights ω i , and for each point M , we denote its components by (x M , y M ). Hence, we have three equations and p unknowns. In general, p = 3. We solve the system using a least squares procedure [START_REF] Bertolazzi | On vertex reconstructions for cell-centered finite volume approximations of 2d anisotropic diffusion problems[END_REF]. In general, p ≥ 3 (for a structured mesh of quadrangles, p = 4 when P is an interior node), and we thus find a solution to this system. However, as we have seen, it is important that u P ≥ 0 in order for the scheme to be well-defined (see Proposition 2.3). A way to ensure this is to have only non-negative weights in the interpolation procedure. This is always the case if P is in the convex envelop of the points K i . When the mesh is highly deformed, this may not be the case (see Figure 4). It is possible to circumvent this difficulty by extending the interpolation to non-neighbouring cells. However, by doing so, we extend the stencil of the scheme, which we want to avoid do to parallelization concerns. At this point, we see two main possible strategies:

-first compute the orthogonal projection Q of P onto the convex hull of the points K i . Then use this projection Q in system (32) instead of P itself. Doing so, you are sure that ω i ≥ 0 for all i, so that u P ≥ 0. However, the projection may result in poor accuracy, since Q might be rather far from P .

-Another possibility is to first solve system [START_REF] Lipnikov | Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes[END_REF], and compute u P , but then truncate the result:

u P = max 0, P i=1 ω i u Ki .
Doing so, we are sure that u P ≥ 0, and in principle we do not affect the accuracy, since we know, a posteriori, that u P ≥ 0. Moreover, contrary to what happens when truncating cell values in the same way, the scheme is still conservative. Indeed, this is a property of the matrix which is not modified by the present truncation.

This last method is the one we have retained, since it seems to us the best compromise between robustness and precision, given the fact that we do not want to enlarge the stencil of the scheme as for the strategy proposed in [START_REF] Sheng | Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes[END_REF].

Fixed-point iteration

As we have seen above, the system is nonlinear (the matrix M in (20) depends on u n+1 ). Hence in order to solve [START_REF] Helmer | A family of MPFA finite-volume schemes with full pressure support for the general tensor pressure equation on cell-centered triangular grids[END_REF], one needs to use a fixed point strategy. For instance, the following algorithm is the most simple and most natural one:

v 0 = u n , ∀k ≥ 0, 1 ∆t v k+1 -u n + M (v k )v k+1 = f + g, ( 33
)
with a stopping criterion of the form v k+1 -v k ≤ ε v k . The value of u n+1 is then given by v k+1 .

We have the following convergence result: Proposition 3.1 Assume that f + g ≥ 0 and u n > 0. Then, if ∆t is small enough, the sequence (v k ) k∈N defined by (33) converges to u n+1 , which is a solution to [START_REF] Helmer | A family of MPFA finite-volume schemes with full pressure support for the general tensor pressure equation on cell-centered triangular grids[END_REF].

Proof: Here again, we use the mapping φ defined in the proof of Proposition 2.3 :

φ : C -→ C u -→ (Id + ∆tM (u)) -1 (u n + ∆t(f + g)) .
The set C is defined by [START_REF] Kershaw | Differencing of the diffusion equation in Lagrangian hydrodynamic codes[END_REF]. We are going to prove that this map is a contraction, so that the iteration v k+1 = φ(v k ) converges to a fixed point of φ.

In order to do so, we first prove the following fact: if ∆t is small enough, then we have

∀u ∈ C, ∀v ∈ C, ∀x ∈ R N , Id + ∆t 2 (M (u) + M (v)) x ≥ 1 2 x . (34)
To prove our claim, we note that, in view of ( 18) and ( 19), we have

|A K,e | ≤ 2 |e|κ e min ---→ KM 1 , ---→ KM 2 ,
where M 1 and M 2 are the endpoints of e. This upper bound, in particular, does not depend on u. Inserting this into ( 21) and ( 22), we have

|[M (u)] KK | ≤ e∈K   2 |e|κ e min ---→ KM 1 , ---→ KM 2   + e∈K∩∂Ω δ e γ e ,
and

|[M (u)] KL | ≤ 2 |e|κ e min --→ LM 1 , --→ LM 2
, All these bounds are independent of u, hence, we have proved that there is a constant C 1 , which depends only on the mesh, and on δ, γ and κ, such that

M (u) ≤ C 1 . Hence, if 2C 1 ∆t ≤ 1, we have Id + ∆t 2 (M (u) + M (v)) x ≥ x - ∆t 2 (M (u) + M (v))x ≥ x - ∆t 2 2C 1 x ≥ 1 2 x ,
which proves [START_REF] Lipnikov | Anderson acceleration for nonlinear finite volume scheme for advection-diffusion problems[END_REF]. Next, we prove the following fact: there exists a constant C 2 depending only on the mesh and on δ, γ and κ, such that,

∀u ∈ C, ∀v ∈ C, M (u) -M (v) ≤ C 2 ∆x 2 u -v . (35)
It is clear from the above argument (and from the formulas ( 21) and ( 22)), that it is sufficient to prove this inequality for A K,e . Now, A K,e depends on u only through μ1 and is a linear function of μ1 , so it is sufficient to prove (35) for µ 1 . We thus write :

μ1 (u) = ψ(a 1 (u), a 2 (u)), where ψ(α, β) := β + ∆x 2 α + β + 2∆x 2 .
The function (α, β) → ψ(α, β) is differentiable in the set {α ≥ 0, β ≥ 0}, and

∂ψ ∂α = - β + ∆x 2 (α + β + 2∆x 2 ) 2 , ∂ψ ∂β = α (α + β + 2∆x 2 ) 2 .
Hence, |∇ψ| ≤ 2∆x 2 -1 , thus

|μ 1 (u) -μ1 (v)| = |ψ(a 1 (u), a 1 (u)) -ψ(a 1 (v), a 2 (v))| ≤ 1 2∆x 2 (|a 1 (u) -a 1 (v)| + |a 2 (u) -a 2 (v)|) .
Finally, in view of the definition of a 1 and a 2 , we have

|a i (u) -a i (v)| ≤ C 2 u -v ,
where C 2 is a constant depending only on κ and on the mesh. Hence, we infer

|μ 1 (u) -μ1 (v)| ≤ C 2 2∆x 2 u -v .
This proves [START_REF] Lipnikov | Minimal stencil finite volume scheme with the discrete maximum principle[END_REF] for μ1 , hence for M . We are now in position to prove that φ is a contraction: indeed, we have, from the definition of φ,

Id + ∆t 2 (M (u) + M (v)) (φ(u) -φ(v)) = 1 2 ∆t (M (v) -M (u)) (φ(u) + φ(v)) .
Applying [START_REF] Lipnikov | Anderson acceleration for nonlinear finite volume scheme for advection-diffusion problems[END_REF] on the one hand, and (35) on the other hand, we have

1 2 φ(u) -φ(v) ≤ C 2 ∆t ∆x 2 u -v φ(u) + φ(v) . (36)
Due to the definition of C, we know that φ(u) and φ(v) are bounded independently of u and v, respectively. Hence, choosing ∆t even smaller if necessary, we have proved that φ is a contraction on C. This proves the result.

Remark 5 This result is of limited practical interest for several reasons. First, a careful examination of the constants in the proof indicate that, at best, the condition would read

∆t ≤ min 1 2C 1 , C 2 ∆x 2 , (37)
for some C 1 and C 2 that in principle still depend on the mesh. Hence, we end up with a condition which is the stability condition for an explicit scheme, loosing the benefit of an implicit scheme. Second, the constant C 2 is not explicit and could in principle be large for highly deformed meshes. However, it indicates that, if convergence of the fixed point should fail, decreasing the time step would help. Actually, it should be noted that in the numerical tests we have done, we did not find any situation in which the sequence defined by [START_REF] Lipnikov | A monotone finite volume method for advection-diffusion equations on unstructured polygon meshes[END_REF] does not converge (even with rather large values of ∆t). An alternative strategy would be to use acceleration techniques, such as the one presented in [START_REF] Lipnikov | Anderson acceleration for nonlinear finite volume scheme for advection-diffusion problems[END_REF].

Remark 6

It is also possible to prove a similar result when assuming that u ≥ m > 0. In such a case, roughly speaking, ∆x 2 in (36) would be replaced by m + ∆x 2 , therefore making condition (37) nicer. However, in such a case, it is necessary that the interpolation procedure giving the node unknowns satisfies the following property:

if ∀K ∈ K, u K ≥ m, then ∀N ∈ N , u N ≥ m. ( 38
)
This property is not satisfied by the interpolation method using truncation (see section 3.1).

If the problem we need to solve is really (1)-( 2), then the above fixed-point strategy is necessary to obtain good results. However, the problems we deal with in practice are nonlinear. Hence, problem (1)-( 2) is only an intermediate step in solving the master equation using, for instance, a Newton algorithm. When this algorithm reaches convergence, u n is close to u n+1 , so the convergence of the fixed point algorithm may not be needed. In numerical simulations, we have compared two strategies: the first one is to reach convergence in the fixed point algorithm. The second one consists in dropping the fixed points strategy, that is, use only one iteration. It happens that, in most cases, it is less costly to use the second strategy. This amounts to using the scheme [START_REF] Le | A linear scheme satisfying a maximum principle for anisotropic diffusion operators on distorted grids. (Un schéma linéaire vérifiant le principe du maximum pour des opérateurs de diffusion très anisotropes sur des maillages déformés[END_REF]. The precision of the two methods is comparable. Here, Remarks 4 and 3 are of particular interest: this strategy does not change the main properties of the scheme.

Numerical tests 4.1 Steady analytic problem

We first assess the scheme on a steady problem with an analytically solution. Since the positivity of the solution is a requirement for our scheme, we choose a solution positive everywhere:

u(x, y) = 1 + cos(πx) cos(πy) ( 39
)
This solution is obtained in solving the following equation in Ω = (0, 1) 2

-div (∇u) = 2π 2 cos(πx) cos(πy), in Ω, (∇u)

• n = 0 on ∂Ω, (40)
with the constraint Ω u = 1. We perform this problem on several meshes, which are classical for diffusion operator. These meshes are obtained by perturbing an initial uniform mesh made of n 2 squares of size 1/n 2 :

-Random mesh -all vertex are translated by the vector 0.4(r x /n, r y /n), where r x and r y are two random numbers uniformly distributed on [0, 1]. This mesh includes non-convex zones, see left part of fig. 5.

-Kershaw mesh [START_REF] Kershaw | Differencing of the diffusion equation in Lagrangian hydrodynamic codes[END_REF] -mesh is highly skewed, see middle part of fig. 5. -Strongly non-convex mesh, constructed by moving the central vertex of every the four-cell patterns of the Cartesian mesh of a vector (0.8 √ 2/n, 0.8 √ 2/n), see right part of fig. 5.

We perform a convergence study on this problem and for these meshes, starting with the 8 2 mesh to a 512 2 mesh. We also compute the ratio

c r = e∈E F K,e (u K -u L ) e∈E κe de |e| (u K -u L ) 2 .
to measure numerically the coercivity of the scheme. Results are summarized in Table The scheme achieves almost second order accuracy even on these highly deformed meshes. In the case of non-convex mesh (last column), it seems that the interpolation procedure we use degrades the order. The coefficient c r , which accounts for the coercivity, converges for all the grid patterns.

Planar nonlinear heat wave in a cold wall

A well-known analytic solution of the system (1) is referred in the literature as Marshak waves. This self-similar solution has been at first exhibited by Marshak [START_REF] Marshak | Effect of radiation on shock wave behavior[END_REF]. It consists in considering a diffusion coefficient κ = κ 0 u k , where κ 0 is a constant and k ∈ N + , a boundary condition u(0, t) = u 0 and a initial condition u(x, t) = 0, x > 0. This generates a self-similar wave propagating into the half space x > 0.

Introducing dimensionless variables

µ = µ(ξ) = u u 0 , ξ = k + 1 2κ 0 u k 0 1/2 x √ u , (41)
the first equation of system (1) reduces to an ordinary differential equation

d 2 µ k+1 dξ 2 + ξ dµ dξ = 0. (42)
The boundary conditions associated with equation ( 42) are µ(0) = 1 and µ = µ k (dµ/dξ) = 0 for ξ = ξ 0 . The position ξ 0 is the solution of the boundary value problem. For our application, we choose k = 3, κ 0 = 1, u 0 = 1 and run the code until t = 1. We perform this calculation on 4 Kershaw meshes. Domain size is 1 × 1, and the number of zones goes from 32 2 zones for the coarsest mesh to 256 2 zones for the finest one. Numerical results are displayed on figures 6 and 7. On these figures, we have plotted the values of the temperature for all the zones of the mesh, as a function of x. The oscillatory aspect of the plot is due to the fact that the mesh is not 1D. We illustrate this point on figure 8 for the coarsest grid. We see on this plot that the solution has no oscillation.

We observe convergence to the analytic solution. Moreover, the solution remains positive all along the simulation. To have an idea of the convergence rate, we compute the L ∞ norm of the function u 4 , for which the spatial derivative is bounded. We obtain a convergence rate of ≈ 1.2.

Radiative shock

Our last test problem is more challenging. The goal is to assess the robustness and accuracy of the method, in a configuration which is relevant for our applications. In theses cases, radiation is strongly coupled with the Euler equations for hydrodynamics. Since we use Arbitrary Lagrangian-Eulerian methods to discretize the Euler equations, the mesh moves at each time-step of the calculation, and the distortion of the zones is induced by the fluid motion.

In the following, we consider the following set of equations corresponding to the Euler equation for the hydrodynamics coupled with an equilibrium diffusion model for the radiation. We solve the equations in the Lagrangian frame:

          ρd t τ = div u, ρd t u = -∇(p + p r ), ρd t (e + e r ) = -div ((p + p r )u) + div c 3σ R ∇e r , ρd t (ε + e r ) = -div (p r u) + div c 3σ R ∇e r , (43)
where d t ϕ ≡ (∂ t + u • ∇) ϕ is the Lagrangian derivative, ρ the density, τ = 1/ρ the specific volume, u the fluid velocity, p the material pressure, p r = 4/3aT 4 the radiative pressure, T the temperature, e = ε + 1/2u 2 the total material energy, e r = aT 4 the radiative energy, c the speed of light, σ R (T ) the Rosseland opacity and ε(ρ, T ) the internal energy.

To solve system (43), we split it into a hydrodynamic part and a radiative part. The hydrodynamic part is solved using the Godunov-type Lagrangian scheme Glace [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Després | Lagrangian gas dynamics in 2D and lagrangian systems[END_REF]. The radiative part consists in solving the following equation:

ρ∂ t (ε + e r ) = div c 3σ R ∇e r . (44)
Note that this equation is nonlinear even if σ R does not depend on T . We use a Newton procedure to solve the equation [START_REF] Plemmons | M-matrix characterizations.i -nonsingular m-matrices[END_REF]. As for the Marshak wave, we use the procedure described in section 3.2 in order to converge at the same time the non-linearity of the system and of the scheme.

In order to properly assess the capability of our code to compute such flows, we compare our results to the semi-analytic solution calculated by Lowrie and Rauenzahn [START_REF] Lowrie | Radiative shock solutions in the equilibrium diffusion limit[END_REF]. We refer to this work for all details concerning the semi-analytic solution. We choose the configuration so that the shock is super-critical. The structure of the solution consists in a radiation precursor and then an isothermal, hydrodynamic shock.

We use the following values ahead of the shock: ρ 0 = 1, u 0 = (1.58551e8, 0., 0.) T , T 0 = 1.53783e6, σ = 485.8. Moreover, we use perfect gas equation of state, with γ = 5/3. This gives the following values for the main parameters of the simulation

P 0 = aT 4 0 ρa 2 0 = 10 -4 , M 0 = |u 0 | -s shock a 0 = 10 and K = aT 4 0 c 3σρ 0 a 3 0 = 10 -4 ,
with a 0 = γp 0 /ρ 0 the sound speed ahead of the shock and s shock the shock speed. Lowrie and Rauenzahn show that these three numbers are sufficient to determine entirely the self-similar solution of the problem. M 0 is the shock Mach number, K controls the amount of thermal diffusion and consequently the extent of the radiation precursor in front of the To assess the scheme for this kind of problem, we artificially extend the calculation in the y-direction: y ∈ [0, 0.1], and perform a convergence study. Number of zones goes from 64 × 10 to 512 × 80 and we apply a Kershaw like pattern on the initial mesh. Resulting coarsest mesh (64 × 10) is depicted in figure 9.

During the calculation, the Lagrangian motion of the mesh deforms the zones giving rise of a more perturbed mesh (see figure 10).

We display on figures 11 to 14 the results of the simulations on the different meshes in term of velocity and temperature, compared to the analytic solution. All the values of the different layers are displayed on these figures.

It shows the convergence of the numerical simulation to the semi-analytic solution, assessing the global hydro-radiative scheme to capture such radiative shock configuration. Moreover, we don't observe any undershoots on the temperature profiles (undershoots on velocity are due to wall-heating effect [START_REF] Noh | Errors for calculations of strong shocks using artificial viscosity and an artificial heat flux[END_REF]). On these meshes, the numerical cost of the nonlinearity of the scheme is weak, because the non-linearity of the problem is as hard to solve than the non-linearity of the scheme.

Conclusion

We have presented in this article a finite-volume scheme for diffusion equation which has a two-point stencil, which is consistent, and is positive, even on deformed meshes. Therefore, assuming the solution is positive, we propose a method to ensure the positivity of the scheme, while keeping a very compact stencil. Consequently, the scheme is easy to use on a parallel computing architecture.

We have conducted numerical tests which assess convergence on highly deformed mesh for smooth solutions. If the deformation is not too strong, second-order accuracy is achieved. Further, we have observed convergence for nonlinear Marshak waves on the one hand, and for hydro-radiative shock solutions on the other hand. In the first case, analytic solutions are available, while in the second case, semi-analytic solutions are used for comparison.

The main drawback is that the scheme is non-linear. Therefore, a fixed-point loop is in principle necessary for computing the solution. We prove that if ∆t is small enough, this algorithm converges.

The fact that the scheme is non-linear indicates that it is not efficient for solving simple linear problems. However, in the case of nonlinear problems, the outer loop (in our case a Newton algorithm) allows for dropping the fixed-point strategy. Indeed, at convergence of the Newton algorithm, two successive values of u, namely u n and u n+1 , are very close, so only one iteration of the fixed-point strategy is necessary. Hence, for non-linear problems, the cost of the present scheme remains comparable to that of linear schemes.

Fig. 5 8 2

 2 Fig. 5 8 2 zones meshes: random (left), Kershaw (middle), non-convex (right).

Fig. 6

 6 Fig. 6 Plain line: analytic solution, symbols: numerical solution. Left 32 2 Kershaw mesh, right 64 2 Kershaw mesh.

Fig. 7

 7 Fig. 7 Plain line: analytic solution, symbols: numerical solution. Left 128 2 Kershaw mesh, right 256 2 Kershaw mesh.

Fig. 8

 8 Fig. 8 Black line: analytic solution, color lines: numerical solutions on the 32 2 Kershaw mesh. Each color corresponds to the colored line on the Kershaw mesh on the right part of the figure, where the numerical solution is measured.



Fig. 9

 9 Fig. 9 Initial coarsest mesh for the radiative shock problem.

Fig. 10

 10 Fig. 10 Final coarsest mesh for the radiative shock problem.

Fig. 11

 11 Fig. 11 1D-plot of the numerical solution (symbols) compared to analytic (plain line) for the 64 × 10 zones mesh. Left: velocity, right: temperature.

Fig. 12

 12 Fig. 12 1D-plot of the numerical solution (symbols) compared to analytic (plain line) for the 128 × 20 zones mesh. Left: velocity, right: temperature.

Fig. 13

 13 Fig. 13 1D-plot of the numerical solution (symbols) compared to analytic (plain line) for the 256 × 40 zones mesh. Left: velocity, right: temperature.

Fig. 14

 14 Fig. 14 1D-plot of the numerical solution (symbols) compared to analytic (plain line) for the 512 × 80 zones mesh. Left: velocity, right: temperature.

Table 1

 1 10 -3 2.27 0.81 4.12 × 10 -2 0.83 0.26 1.03 × 10 -2 1.92 0.15 32 5.18 × 10 -4 2.10 0.80 1.91 × 10 -2 1.11 0.22 2.61 × 10 -3 1.98 0.13 64 1.11 × 10 -4 2.22 0.79 7.35 × 10 -3 1.38 0.19 7.14 × 10 -4 1.87 0.12 128 2.69 × 10 -5 2.04 0.79 2.27 × 10 -3 1.70 0.18 1.99 × 10 -4 1.84 0.12 256 9.57 × 10 -6 1.49 0.79 6.18 × 10 -4 1.88 0.18 5.82 × 10 -5 1.77 0.12 512 2.69 × 10 -6 1.83 0.79 1.60 × 10 -4 1.96 0.18 1.84 × 10 -5 1.66 0.12 Convergence table for the analytic problem

	1.

Copyright line will be provided by the publisher

Acknowledgments

The authors address many thanks to G. Samba and F. Hermeline for useful advises and fruitful discussions about this work.