Dan Chiorean
email: chiorean@cs.ubbcluj.ro

Benoit Combemale
email: benoit.combemale@irisa.fr

The 8th Educators' Symposium @ MODELS 2012: Software Modeling in Education

Keywords: A.0 [GENERAL]: Conference Proceedings, K.3.2 [Computer and Information Science Education] Computer science education, Curriculum Management, Documentation, Design, Human Factors Software Engineering Education, Introductory and Advanced Modeling Courses, successful teaching experiences

Modeling now plays an important role in software-intensive systems development and evolution. Modeling provides goaloriented abstractions in all phases of the system lifecycle, which requires deep knowledge on modeling techniques and broad experiences in applying these techniques. The Educators' Symposium is organized annually in the framework of the international conference MODELS for discussing teaching these technologies to software engineers at universities and software industries.

In this paper, we summarize the 8th edition of the Educators' Symposium: EduSymp 2012. First, we remind the general objectives and goals of the Symposium, and we describe the general organization of the 2012 edition (committees, program, statistics, and publication process). Then we introduce the papers selected by the program committee to be presented during the Symposium, and to be published into the proceedings. Finally, we provide a summary of presentations and discussions that occurred during the keynote as well as panel of EduSymp 2012.

MOTIVATION AND GOALS

"Modeling of software is becoming a pervasive technique to help software engineers understand, engineer, and communicate aspects of the software to appropriate stakeholders" [START_REF]Guide to the Software Engineering Body of Knowledge (SWEBOK) V3[END_REF]. However, "...MDE is still in the early adoption phase and to be successfully adopted by industry, it must prove its superiority over other development paradigms and be supported by a rich ecosystem of stable, compatible and standardized tools. It should also not introduce more complexity than it removes" [START_REF] Mohagheghi | MDE Adoption in Industry: Challenges and Success Criteria in Models in Software Engineering[END_REF].

As Bezivin remarked: "In never-ending Research-Development-Teaching cycles, MDE is now in the position where the teaching question is probably the more acute" [START_REF] Bezivin | Teaching Modeling: Why, When, What? in Models in Software Engineering[END_REF]. Many students have misconceptions or a superficial impression about modeling: "... they tend to view software modeling with great skepticism" and "... often feel that modeling adds accidental complexity to the software development process as they perceive it" [START_REF] France | Teaching Programming Students how to Model: Challenges & Opportunities, invited speak at EduSymp[END_REF].

The rate of success in case of companies that uses MDE techniques is very different from a company to another. Sensible differences in the results obtained in the teaching process were reported by academics in case of using different teaching methods and even when the order of topics introduced or associated was different. A first conclusion of the research presented in [START_REF] Whittle | What Factors Lead to Industrial Success or Failure with MDE[END_REF] is that "... the factors relating to whether new technologies succeed or fail are more often social or organizational rather than technical. This is true as well for MDE ... ".

In this context, offering teachers the opportunity to share both the more or less successful teaching experiences in the framework of MODELS conference became a tradition.

The paper is organized as follow. We give in Section 2 some insights on the organization of the edition 2012 of EduSymp. Then we present in Section 3 the accepted papers included in the post-proceedings as well as some notes in Section 4 on the invited talk and the final panel. Finally we conclude with personal thoughts on the symposium and current trends in the teaching of modeling.

ORGANIZATION

Dan Chiorean and Benoit Combemale organized the edition 2012. They were also co-chairs of the Program Committee. Publishing the proceedings by a well recognized international editor and announcing as earlier as possible this symposium edition were among the concerns of organizers. In previous editions the above-mentioned aspects were mentioned as potential factors diminishing the interest of submitting proposals at EduSymp. Organizers are very grateful to the conference general chairs for succeeding to enable the publication of all satellite events post-proceedings in the ACM DL. The symposium website 1 and the call for papers (CfP) was online, quite 10 months before the symposium date. Apart of the symposium site, the CfP was announced by different professional mailing lists (e.g., SEWORLD, pUML, planetmde). Despite the compliance with these "preconditions" the number of submissions received was in our opinion modest. We received 10 abstracts and 9 final papers. Finally 5 of these submissions were accepted (2 full papers and 3 short papers), resulting an acceptance rate of 55%. Each submission was reviewed by at least three members of the Program Committee. Based on the reviews, the decisions regarding papers acceptance were taken unanimously. Three papers were co-authored by one of the workshop PC members. The review process ensured that the authors had no influence on the acceptance/rejection decision for papers written by them. In addition to the workshop organizers, the PC consisted of: Colin Atkinson -University of Mannheim, Thomas Baar

ACCEPTED PAPERS

The first long paper "Model-Driven Paradigms -The Evolution of a University Course" -was written by László Lengyel and Gergely Mezei from Budapest University of Technology and Economics [START_REF] Lengyel | Model-Driven Paradigms -The Evolution of a University Course[END_REF]. It is about a master course entitled Model-Driven Paradigms, grouping topics from the followings domains: Domain-Specific Modeling, Model Processing and Software Methodologies. The topics were organized in 11 themes each having allocated a percent of the whole course. The main objectives of authors were to attain a higher satisfaction of students taking this course and better results at exams. In order to comply with these objectives the percent of different themes was changed, and the obtained results measured. The study was realized in the period 2009-2012, and the changes triggered by students' feedback were made in 2009, 2010 and 2012. The conclusion resulted highlight the importance of a correct balance between different themes. The second long paper, "Replacing Traditional Classroom Lectures with Lecture Videos -An Experience Report" was authored by Marion Brandsteidl, Tanja Mayerhofer, Martina Seidl and Christian Huemer from Institute of Software Technology and Interactive Systems Vienna and Vienna University of Technology [START_REF] Brandsteidl | Teaching Models @ BIG: Replacing Traditional Classroom Lectures with Lecture Videos -An Experience Report[END_REF]. As the title clearly suggest, the paper analyze the advantages obtained by replacing traditional classroom lectures with video lectures. The price to pay for obtaining the above mentioned advantages are also described. This time, it's about an introductory course in Object-Oriented Modeling offered by the Business Informatics Group (BIG) at the Vienna University of Technology. In our opinion, the large preference of students for video lectures compared with classical lectures is an explicit message to teachers that these modern techniques have to be considered in the future. Even if the price of producing a video lecture can be expensive for the first time, the R. O. I. is excellent. Students can consume each lecture anytime and at any place as often as they want. On the other hand, teachers are not forced to give the same lecture many times, and have the opportunity to use their time more efficient.

The three short papers were: "Improving Software Engineering Education by Modeling Real-World Implementations" [START_REF] Baar | Improving Software Engineering Education by Modeling Real-World Implementations[END_REF] by Thomas Baar -Hochschule für Technik und Wirtschaft (HTW) Berlin, "How Should Teaching Modeling and Programming Intertwine?" [START_REF] Demuth | How should Teaching Modelling and Programming Intertwine[END_REF] by Birgit Demuth -Technische Universität Dresden and "On Explaining Modeling Principles with Modeling Examples: A Classification Catalog" [START_REF] Gogolla | On Explaining Modeling Principles with Modeling[END_REF] written by Martin Gogolla from University of Bremen and Antonio Vallecillo from University of Malaga. Similar to proposals described in the long papers, Thomas Baar proposal was also triggered by students' feedback. If the connection between the model and the system is not enough clear, students become demotivated. This suggested Thomas to propose the use of real world models that can provide different opportunities like: to see the advantages of using models for medium and large systems and to have samples of efficient implementations for real models. In the section "Related work" the author analyze the well known experience of Repository for Model Driven Development (ReMoDD) and remarks that real models contained in ReMoDD are not enough detailed in order to be efficiently used in code generation. In our opinion, this is due to the fact that real models are proprietary models and, by consequence the owners try to protect them. In Birgit Demuth's paper, the author presents the conclusion of 15th year of experience in teaching modeling with UML and OO programming at Technische Universität Dresden. The main educator's question is: how to teach students in order to prepare them as better as possible both in programming and modeling. The conclusion confirms Robert France opinion, that modeling and programming must be teach intertwined. In this manner, students understand easily the utility of modeling in software engineering. Regarding this paper, it worth to mention that the results obtained by students represents the feedback and, in Dresden experiment were included two related courses: introduction in OO programming, analysis and design using Java and UML and a practical course of applying the knowledge acquired in analyzing, designing and implementing a medium size application working in a team. Finally, the paper written by Gogolla and Vallecillo proposes a topic different from all others: a catalogue in which examples are classified by rapport to their syntax, semantics, pragmatics, complexity or evolution. These criteria enable teachers to check if examples included in a curricula cover the concepts and principles introduced and analyzed in the course. The proposal will be validated and updated in future projects.

THE INVITED TALK AND THE PANEL

This year, the keynote talk entitled "Model-based Development in Large Project Courses with Real Customers" was given by professor Bernd Brügge from Technische Universität München. Based on the results obtained in an impressive number of real projects, Professor Brügge argued that one can teach industry relevant software engineering practices to students with a special emphasis on informal models the key of a successful communication between developers and users/clients. The speaker explained how to combine modeling with delivery and presentation activities, ranging from scenario-based design, requirements elicitation to the presentation of executable prototypes. The software lifecycle model is based on a hybrid process model using a combination of the Unified Process and Scrum. The presentation was very well received and has sparked a lot of questions and very interesting discussions. The keynote slides are posted on the symposium site http://edusymp2012.irisa.fr/. The paper entitled "Teaching the Tornado: From Communication Models to Releases" [START_REF] Bruegge | Teaching the Tornado -From Communication Models to Releases in Software Modeling[END_REF] describing the keynote was included in the EduSymp 2012 postproceedings.

The panel entitled: "How do we inspire students to model?" was focused on measures meant to support teachers in convincing students about the usefulness of understanding, learning and using modeling in Software Engineering. After a short introduction of panelists and of potential topics given by the moderator, the panelists presented their position. Finally, the participants addressed their questions to panelists and/or made comments related to different topics. Both the introductory slides and panelists' slides are available on the symposium website. Compared to the previous EduSymp panel organized in 2009 and entitled: "Teaching Modeling: Why, When, What?" this panel was focused on a narrow theme. Panelists positions were very appreciated. Some of the most important statements:

Colin Atkinson: Modeling is an integral, indispensible, essential part of effective SE. What we call models today are just views/parts of a subject. Code is just a view as well, and therefore also a model. All development is or should be modelaware.

Bernd Brügge: A problem solving method that works -break down recursive the problem into two or more sub-problems of the same type, until these become simple enough to be solved directly. The solutions to the sub-problems are then combined to give a solution to the original problem. The ability to understand and apply divide and conquer is a skill that takes time to master: Practice, practice, practice...

Jean Michel Bruel:

Considers that among the most appropriate manners for "inspiring" students to model is to ask them to do different activities like designing algorithms or coding patterns, (that are) natural based on modeling, without mentioning that this are modeling activities. Jean Michel explicitly stated about the difference of learning a modeling language, like UML and learning modeling.

Robert France:

Modeling is an intrinsic human capability. We all model. Often we use only implicitly held models to help understand and explain concepts, or to drive our decisionmaking process. The challenge for educators is how do we further develop/nurture this capability in software engineering students.

Jeff Gray: Focuses on "Big Ideas" of modeling that motivate its usage in a way that is more appealing than learning about the specific semantics of a modeling language.

Bran Selic: Convince students to use modeling by embedding MBE in other courses with a clearer value proposition e.g.: courses on software architecture, advanced programming, systems design. Teach essentials -not the syntax. Invite experts as guest lecturers to teach topics as design principles, metamodeling a. s. o. Use executable modeling tools and propose moderately complex team-based projects. Advice the use of tools easy to install and having effective tutorials.

CONCLUSION

This edition of EduSymp was a successful one both as concerning the topics debated, and the quality of presentations. The keynote focused on the original experience in teaching Software Engineering of Professor Brügge, at CMU and Technische Universität München. For sure, this experience deserves to be disseminated in other universities. The accepted papers were focused on: critic comparisons among the ponder that different topics have in the curriculum of introductory and advanced courses about modeling in software engineering, using new communication technologies in preparing and giving courses, the importance of using real world examples, classification criteria for examples -meant to support teachers in evaluating the degree in which the examples cover the concepts and principles discussed in the course. Finally, the panelists concluded that in software engineering, all development should be model-aware, that similar to teaching programming, teaching modeling have to be based on big ideas, not on secondary aspects and mentioned different tips and tricks meant to support teachers in succeeding to inspire students to model. Hoping that in the future the interest in presenting and sharing teaching experiences in modeling will grow offering so a better access to interesting examples and successful teaching experiences.

  -University of Applied Sciences Berlin, Marion Brandsteidl -Vienna University of Technology, Bernd Brügge -Technische Universität München, Peter J. Clarke -Florida International University, Birgit Demuth -TU Dresden, Robert B. France -Colorado State University, Martin Gogolla -University of Bremen, Jeff Gray -University of Alabama, Thomas Kühne -Victoria University of Wellington, Ludwik Kuzniarz -Blekinge Institute of Technology, Richard Paige -University of York, Alfonso Pierantonio -Università degli Studi dell'Aquila, Pascal Roques -PRFC France, Bran Selic -University of Toronto, Andreas Winter -Carl von Ossietzky University. The organizers are very grateful to PC members for accepting to do this job in favor of the scientific community and for the quality of their work. In order to increase the interest in this EduSymp, and in conformance with the tradition, just after opening remarks presented by Benoit Combemale, Bernd Bruegge presented an invited keynote talk. Similarly, in the last afternoon session, a panel discussion moderated by Dan Chiorean and having: Colin Atkinson -University of Manheim, Bernd Brügge -Technische Universität München, Jean Michel Bruel -University of Toulouse, Robert B. France -Colorado State University, Jeff Gray -University of Alabama and Bran Selic -Malina Software and University of Toronto as panelists was organized.


C Calgaro

E Creusé
email: creuse@math.univ-lille1.fr

T Goudon

Y Penel

Positivity-preserving schemes for Euler equations: sharp and practical CFL conditions

Keywords: Positivity-preserving MUSCL schemes, 2nd-order finite volume method, Euler equations, CFL condition, τ -limiters

When one solves PDEs modelling physical phenomena, it is of great importance to take physical constraints into account. More precisely, numerical schemes have to be designed such that discrete solutions satisfy the same constraints as exact solutions. For instance, the underlying physical assumptions for the Euler equations are the positivity of both density and pressure variables.

We consider in this paper an unstructured vertex-based tesselation in R 2 . Given a MUSCL finite volume scheme and given a reconstruction method (including a limiting process), the point is to determine whether the overall scheme ensures the positivity. The present work is issued from seminal papers from Perthame & Shu and Berthon. They proved in different frameworks that under assumptions on the corresponding one-dimensional numerical flux, a suitable CFL condition guarantees that density and pressure remain positive.

We first analyse Berthon's method by presenting the ins and outs. We then propose a more general approach adding non geometric degrees of freedom. This approach includes an optimization procedure in order to make the CFL condition explicit and as less restrictive as possible. The reconstruction method is handled independently by means of τ -limiters and of an additional damping parameter. An algorithm is provided in order to specify the adjustments to make in a preexisting code based on a certain numerical flux. Numerical simulations are carried out to prove the accuracy of the method and its ability to deal with low densities and pressures.

Introduction

Partial differential equations are widely used in physics modelling and are thus expected to reproduce actual situations. In particular, most models in fluid dynamics consist in systems of conservation laws whose solutions (must) satisfy some properties: from a physical point of view to make sense (for instance, we may think of the positivity of density or fluid height) and from the point of view of mathematical analysis for the problem to be well-posed.

The transition to industrial codes must take these elements into account in addition to other considerations such as accuracy and computational costs. The design of numerical schemes consists in a balance between these antagonistic criteria since the quest for a more and more accurate solution may increase the number of unknowns and thus the computational time. As for physical constraints which include maximum principle and positivity preservation, it is important to bear in mind that they may ensure discrete problems to be solvable at each time step.

Finite volume methods (FVm) seem to suit quite well to computational fluid dynamics since they stick to the derivation of equations from basic physical principles. In addition, they enable to simulate both classical and weak solutions as well as to handle general unstructured meshes. It is worth underlying that the term FVm may relate to very different approaches [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF]. FVm may differ from one another by the location of unknowns (cell-centered vs. vertex-based ); notice that the issue of which one is the "best" is still open [START_REF] Mavriplis | Unstructured Mesh Discretizations and Solvers for Computational Aerodynamics[END_REF], especially in the case of irregular meshes where major drawbacks may appear. For instance, in vertex-based approaches, the colocation point is generally not the centroid of the control volume which prevents from identifying mean values to pointwise values and may induce a large discrepency for this kind of meshes. Similarly in the framework of MUSCL techniques and multidirectional slope procedures, cell-centered methods can lead to severe CFL restrictions as well as to a loss of accuracy in the case of mesh distortions (see e.g. [START_REF] Clain | L ∞ -stability of the MUSCL methods[END_REF], particularly Remark 3 page 38).

Other distinctions between MUSCL procedures are for instance the reconstruction step or the numerical flux -see for example [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] for hyperbolic systems of conservation laws.

To illustrate what was introduced above, let us consider the following system of conservation laws:

∂ t W + ∇ • F(W) = 0, (t, x) ∈ [0, T ] × Ω, W(0, x) = W 0 (x), x ∈ Ω, (1)
for some time T > 0 and for a bounded domain Ω ⊂ R d together with suitable boundary conditions. The question of boundary conditions supplementing such hyperbolic systems is quite delicate and subtle. This is due to the fact that the number and the nature of boundary conditions depend on local properties of the solution itself and it might change with time (specializing to gas dynamics, it depends whether the flow is sub or supersonic). Hence, the difficulty is twofold: firstly one has to design physically relevant boundary conditions which guarantees well-posedness to the initial boundary value problem. About this point we refer to [START_REF] Aubert | Boundary conditions for Euler equations[END_REF] where practical criteria are detailed. Secondly the numerical scheme has to treat appropriately these conditions. Here we shall thus work with appropriate in-coming fields.

Then the solution W would be supposed to lie in a given set of admissible states W. For instance, for the Euler equations, W includes density, momentum and pressure; W is then the set of vectors W for which density and pressure are positive. From a theoretical point of view (well-posedness, relevance), it is particularly important to study the invariance of W with respect to Syst. (1) [START_REF] Serre | Systems of conservation laws. 2. Geometric structures, oscillations, and initial-boundary value problems[END_REF][START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF]: if W 0 ∈ W, do we have W(t, •) ∈ W for all t > 0?

This question being answered, it is relevant to study its discrete version. The present work deals with vertex-based strategies insofar as we analyze and modify ideas from works restricted to that framework [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF][START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF][START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF][START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]. We thus introduce a vertex-based tesselation of Ω made of triangles whose nodes are denoted by (M i ) 1≤i≤Nn . The control volume Ω i corresponding to M i is built joining barycenters of cells having M i as a vertex -see Fig. 1. V(i) is the set of neighbours' indices, Γ ij is the interface between Ω i and Ω j and n ij the unit normal vector to Γ ij oriented from Ω i to Ω j . Moreover, we assume that the mesh is smooth so that:

∀ i = j, Γ ij ∩ [M i M j ] = ∅. ( H0
)
If M i is located on the boundary, the associated control volume Ω i connects barycenters to the midpoints of the boundary edges. Likewise, the time interval is discretized by means of a sequence of time steps ∆t n . We thus set t 0 = 0 and t n+1 = t n + ∆t n .

To solve [START_REF] Aubert | Boundary conditions for Euler equations[END_REF], we can use a vertex-based finite-volume scheme:

W n+1 i = W n i -∆t n j∈V(i) |Γ ij | |Ω i | F W n i , W n j , n ij ,
where W n i is a cell-wise constant approximation at time t n of the solution over Ω i corresponding to an internal node M i and F is the numerical flux -see further for properties satisfied by F . For details about the numerical handling of the boundary conditions, the reader may refer to [START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF]. The question now reads: if

W n i ∈ W, do we have W n+1 i ∈ W?
The case of several first-order classical schemes (such as Lax-Friedrichs or Siliciu) is handled in Bouchut's review [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF].

With the quest for better accuracy by means of higher order schemes, the question has to be re-addressed. When one uses a MUSCL second-order (in space) scheme as designed by Van Leer [START_REF] Van Leer | Toward the ultimate conservative difference scheme V: A second order sequel of Godunov's methods[END_REF]: where W n ij and W n ji are reconstructed values of the solution in the vicinity of Γ ij , one has to take care of the reconstruction method in the answering process. A standard Runge-Kutta second-order (in time) extension would be straightforward since it uses a decomposition into two first-order schemes to which the present approach applies.

W n+1 i = W n i -∆t n j∈V(i) |Γ ij | |Ω i | F W n ij , W n ji , n ij , (2)
M j M k Mm Mp Mq Γ ij Ω i Ω l M i M l nij
For scalar conservation laws, discrete maximum principle has been extensively studied. For instance, we may refer to [START_REF] Barth | Finite volume methods: foundation and analysis[END_REF] where positive schemes are introduced as well as to [START_REF] Clain | L ∞ -stability of the MUSCL methods[END_REF] for cell-centered MUSCL schemes or to [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF] for vertexbased tesselations. The case of systems of conservation laws seems as usual to be more difficult. In the framework of gas dynamics (Euler equations for ideal gas), a partial answer has been yielded by Cournède et al. [START_REF] Cournède | A positive MUSCL scheme for triangulations[END_REF]. Khobalatte & Perthame [START_REF] Khobalatte | Maximum-principle on the entropy and second-order kinetic schemes[END_REF] and Estivalezes & Villedieu [START_REF] Estivalezes | High-order positivity-preserving kinetic schemes for the compressible Euler equations[END_REF] provided results for kinetic schemes. Then Perthame & Shu [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] (cell-centered), Linde & Roe [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF] (convex control volumes) and Berthon [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF][START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF] (vertex-based) proved that the second-order scheme (2) preserves positivity of density and pressure provided the 1D corresponding numerical flux does. Figuring out whether a 2D scheme preserves positivity thus comes down to studying properties of the 1D scheme. Besides, it is also the case for other kinds of numerical methods. We can for example mention the recent works of Parent applied to the resolution of the multidimensional Euler equations in generalized curvilinear coordinates using the so-called "rule of the positive coefficients" [START_REF] Parent | Positivity-preserving flux-limited method for compressible fluid flow[END_REF][START_REF] Parent | Positivity-preserving high-resolution schemes for systems of conservation laws[END_REF].

Moreover, while it is well known that CFL conditions ensure stability from a numerical point of view, they proved that CFL-like estimates for the time steps also seem to be necessary to preserve positivity. However, the latter conditions for the time steps are more restrictive than usual ones. As they result from sufficient conditions, the issue to know whether they are practically necessary remains open.

Noticing that the procedures in [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF][START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF][START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] do not seem optimal, we propose a slightly different approach in order to optimize the coefficients involved in their proofs and hence to make the CFL condition less restrictive and fully explicit. In Section 2, we recall Berthon's strategy showing that properties of a 1D numerical flux may imply the positivity-preserving property for a 2D finitevolume scheme. In Section 3, we revisit this approach but without any geometric interpretation in order to improve the computation of the time step. No matter what the approach which is used (geometric or algebraic), a central assumption is that the reconstruction step is "nicely" performed (in the sense that the reconstructed values are physically admissible and that a certain intermediate state to be specified exists in W). This will be the point in Section 4. Section 5 is a summary of § § 3 and 4: an analysis of the overall process is carried out leading to a practical procedure. The last part is devoted to numerical results for gas dynamics.

Berthon's approach for the Euler equations

Noting ρ, u = (u, v) and E the density, velocity and energy of a gas, the Euler equations read:

         ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • ρ|u| 2 + p = 0, ∂ t (ρE) + ∇ • (ρE + p)u = 0, (3)
where the pressure is given by the ideal gas equation of state:

p = ρ(γ -1) E - |u| 2 2 ,
γ ∈ [START_REF] Aubert | Boundary conditions for Euler equations[END_REF][START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF] being the ratio of specific heats. Eqs. (3) may be equivalently considered as a system of conservation laws and written like (1) with:

W =     ρ ρu ρv ρE     and F(W) =     ρu ρv ρu 2 + p ρuv ρuv ρv 2 + p (ρE + p)u (ρE + p)v     . (4)
From a physical point of view, the density and pressure variables must be positive which implies that the solution W of (1)-(4) must lie in the set:

W := W ∈ R 4 : ρ = W 1 > 0 and p = (γ -1) W 4 - W 2 2 + W 2 3 2W 1 > 0 .
A partial-order relation may be defined on R 4 to characterize W [START_REF] Estivalezes | High-order positivity-preserving kinetic schemes for the compressible Euler equations[END_REF]. It can also be shown that p = p(W) is homogeneous of degree 1 and a continuous concave function. Consequently, W is an open convex cone in R 4 . In addition, it has been proven that the Euler equations preserve the positivity of density and pressure [START_REF] Serre | Systems of conservation laws. 2. Geometric structures, oscillations, and initial-boundary value problems[END_REF][START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF]. Physically speaking, it is essential that the numerical scheme we wish to use to simulate Syst. ( 1) with (4) do preserve the invariance of W. Any procedure to investigate this property must rely on the convexity of the set of physical states W. A review has been carried out in [START_REF] Zhang | Maximum-principle-satisfying and positivitypreserving high order schemes for conservation laws: Survey and new developments[END_REF] where two main approaches are detailed. One approach is based on polynomial reconstructions (for a Discontinuous Galerkin method) and on quadrature integration formulae (which are convex procedures). Constraints then apply to the values of the polynomial at the quadrature nodes -see the series of papers by Zhang et al. [START_REF] Zhang | On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivitypreserving high order schemes for conservation laws: Survey and new developments[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF].

The other approach has been introduced in [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] and consists in expressing the updated value W n+1 i as a convex combination of states lying in W. Perthame and Shu use a cell-centered tesselation combined to a monoslope1 MUSCL procedure. More precisely, the second-order reconstructed values are computed by means of a linear function within the control volume. This results in the existence of positive coefficients ω ij such that:

W n i = j∈V(i) ω ij W n ij . (5)
The authors then proved that W n+1 i satisfies a similar relation with W n ij replaced by solutions of one-dimensional first order schemes (see below for more details). Assuming that reconstructed values belong to W and that the firstorder scheme preserves positivity leads to the conclusion. The same trick applies for the 2D MUSCL scheme on cartesian grids detailled in [START_REF] Spekreijse | Multigrid solution of monotone second-order discretizations of hyperbolic conservation laws[END_REF] and for convex control volumes in [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF].

This argument no longer holds when using a multislope procedure in the reconstruction step. Berthon adapted this idea to vertex-based meshes and multislope reconstructions in 1D [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF] and in 2D [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF]. We only focus on the 2D case in the sequel. This very general approach deals with the relevant CFL condition to apply given a reconstruction method and a (suitable) numerical flux.

As relations like (5) do not naturally occur in multislope methods, Berthon introduced a new state W * i ∈ W so that W n i is a convex combination of states in W -see [START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF]. But this new variable must be associated to a geometric element. That is why new points m ijp are added on the segment M i G ijp to build the sub-cell Ω * i where G ijp is the barycenter of triangle M i M j M p (see Fig. 2). For the sake of simplicity, we suppose that the ratios r ijp := M i m ijp /M i G ijp are equal to r i within the control volume Ω i , such that the quadrilaterals The convex combination reads in [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF]:

Ω ij := m ijp m ijl G ijl G ijp are trapezoidals. M j Mp M l G ijp G ijl m ijp m ijl M i Ω * i Ω ij
W n i = |Ω * i | |Ω i | W * i + j∈V(i) |Ω ij | |Ω i | W n ij . (6)
Nevertheless, given reconstructed values W n ij ∈ W (see Sect. 4) the existence of a state W * i verifying (6) and belonging to W is not ensured. Indeed, considering the equivalent formulation:

W * i = |Ω i | |Ω * i | W n i - j∈V(i) |Ω ij | |Ω * i | W n ij , (7)
we notice that the right hand side is known and does not necessarily lie in W.

The only way to satisfy this requirement is to modify the reconstruction step so that both "starred" density and pressure turn positive. A brief routine is presented in [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF]. We propose in Section 4 an algorithm that provides suitable reconstructed values. Once W * i ∈ W is computed, the next step consists in expressing W n+1 i as a combination similar to [START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF]. Denoting counter-clockwise Γ ij,k the edges of Ω ij (see Fig. 3) and Γ * ij the edges of Ω * i , Berthon sets: In accordance with the local indices of the edges, the notations are:

             W ij = W n ij -∆t n 4 k=1 |Γ ij,k | |Ω ij | F W n ij , W n ij,k , n ij,k , j ∈ V(i), (8a)
W * i = W * i -∆t n j∈V(i) |Γ * ij | |Ω * i | F W * i , W n ij , n * ij . ( 8b
) G ij,4 = G ijp G ij,2 = G ijl m ij,4 = m ijp m ijl = m ij,2 Z ij Γ ij,4 T ij,4 Γ ij,3 = Γ ij T ij,3 Γ ij,2 T ij,2 Γ ij,1 = Γ * ij T ij,1 nij,3 = nij nij,2 nij,4 (a) Ω ij = 1≤k≤4 T ij,k M i m ij,4 = m ijp m ijl = m ij,2 Γ ij,1 = Γ * ij n * ij = -nij,1 = nij (b) Ω * i = j∈V(i) T * ij
W n ij,1 = W * i , W n ij,2 = W n il , W n ij,3 = W n ji and W n ij,4 = W n ip .
Due to ( 6) and ( 8), a straightforward calculation shows that (2) also reads:

W n+1 i = |Ω * i | |Ω i | W * i + j∈V(i) |Ω ij | |Ω i | W ij . (9)
Relations (8a) and (8b) correspond to 1st-order schemes applied to local variables W n ij and W * i . To prove that W n+1 i belongs to W by means of a convexity argument, it suffices to show that W * i and W ij lie in W. These requirements imply assumptions on the numerical flux F as it will be shown in the sequel. A numerical flux must satisfy the following standard properties [START_REF] Eymard | Finite volume methods[END_REF]:

• Consistency with the physical flux:

∀ W ∈ W, ∀ n ∈ S 2 , F (W, W, n) = F(W) • n;
(H1)

• Conservativity:

∀ (W l , W r , n) ∈ W 2 × S 2 , F (W l , W r , n) = -F (W r , W l , -n); (H2) • Continuity: F is locally Lipschitz-continuous. ( H3
)
In addition to the previous requirements, the proof is based on the invariance of W wrt the flux. More precisely, Berthon proved that properties of the two-dimensional flux F may reduce to properties of the corresponding onedimensional flux in the direction orthogonal to the interface. In the present case, (8a) and (8b) can be decomposed respectively as:

               W ij = 4 k=1 |T ij,k | |Ω ij | W ij,k , (10a)
W ij,k = W n ij -∆t n |Γ ij,k | |T ij,k | F (W n ij , W n ij,k , n ij,k ) -F (W n ij , W n ij , n ij,k ) , (10b)
and:

           W * i = j∈V(i) |T * ij | |Ω * i | W * ij , (11a)
W * ij = W * i -∆t n |Γ * ij | |T * ij | F (W * i , W n ij , n * ij ) -F (W * i , W * i , n * ij ) . (11b)
To derive one-dimensional problems, additional geometric elements were required. A point Z ij ∈ Ω ij is introduced in order to build a splitting of Ω ij in four triangles -see Fig. 3a. Its position is left free. Likewise, Ω * i is split into #V(i) triangles -see Fig. 3b.

The equivalence between (8a) and ( 10) relies on the Green formula and on the consistency assumption (H1):

4 k=1 |Γ ij,k | |Ω ij | F (W n ij , W n ij , n ij,k ) = 4 k=1 |Γ ij,k | |Ω ij | F(W n ij ) • n ij,k = F(W n ij ) |Ω ij | • 4 k=1 |Γ ij,k |n ij,k = F(W n ij ) |Ω ij | • ∂Ωij n ij dσ = 1 |Ω ij | Ωij ∇ • F(W n ij ) =constant dx = 0.
What is interesting in formulations [START_REF] Cournède | A positive MUSCL scheme for triangulations[END_REF] and ( 11) is that W ij,k and W * ij are now solutions of first-order one-dimensional schemes. To lead to the conclusion, it suffices to have:

∀ (W 1 , W 2 , n) ∈ W 2 × S 2 , ∃ ∆t > 0, W 1 - ∆t [F (W 1 , W 2 , n) -F (W 1 , W 1 , n)] ∈ W. ( 12
)
This is obviously satisfied for ∆t small enough since W is open and W 1 ∈ W.

The issue reduces to determining how small ∆t should be which provides the CFL condition. It is possible to go further for most schemes by taking into account the rotational invariance of the physical flux in the Euler equations [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]

, § 3.2.1]. Setting: R n =     1 0 0 0 0 cos θ n sin θ n 0 0 -sin θ n cos θ n 0 0 0 0 1     , θ n = arccos(n • e 1 ),
and:

F n (W) = F(W) • n,
we have:

F n (W) = R -1 n F e1 (R n W).
Notice that W is also invariant by rotation. It is thus natural to assume that the numerical flux has the same property, which is the case for most classical fluxes such as Lax-Friedrichs and Godunov schemes:

2 F n (W l , W r ) := F (W l , W r , n) = R -1 n F e1 (R n W l , R n W r ). Hence we assume that, given (W 1 , W 2 ) ∈ W 2 , if ∆t > 0 satisfies ∆t λ(W 1 , W 2 ) ≤ α 0 , (H4.CFL) then W 1 - ∆t [F e1 (W 1 , W 2 ) -F e1 (W 1 , W 1 )] ∈ W. ( H4
)
Here α 0 > 0 is specific to the flux F e1 . For instance, for the Lax-Friedrichs scheme, α 0 = 1 since in the present situation -see (10b) and (11b) -there is no wave arising on the left side. λ(W l , W r ) denotes the largest eigenvalue of the Riemann problem with numerical flux F e1 and initial data (W l , W r ). In the Euler case (4), we have:

λ(W l , W r ) = max{|u l | + c l , |u r | + c r }, c = γp ρ .
Then (H4) implies ( 12) under the CFL condition (H4.CFL) with λ replaced by λ n :

λ n (W l , W r ) = max{|u l • n| + c l , |u r • n| + c r }.
We are now able to conclude. Given a numerical flux F satisfying (H1-H2-H3-H4) and given reconstructed values W n ij ∈ W such that W * i ∈ W, we deduce from (10b) and (11b) that W ij,k and W * ij belong to W under the following CFL conditions:

∆t n max j∈V(i) 1≤k≤4 |Γ ij,k | |T ij,k | λ n ij,k (W n ij , W n ij,k ) ≤ α 0 , (13a)
∆t n max j∈V(i) |Γ * ij | |T * ij | λ n * ij (W * i , W n ij ) ≤ α 0 . (13b)
Then, W ij and W * i are in W as (10a) and (11a) are convex combinations. For the same reason, [START_REF] Clain | A multislope MUSCL method on unstructured meshes applied to compressible Euler equations for axisymmetric swirling flows[END_REF] shows that W n+1 i ∈ W. To compute explicitly the time step, we have to combine the finite collection of conditions [START_REF] Estivalezes | High-order positivity-preserving kinetic schemes for the compressible Euler equations[END_REF] for each node M i .

In a nutshell, MUSCL scheme (2) guarantees that density and pressure variables remain positive as long as the time step satisfies inequalities [START_REF] Estivalezes | High-order positivity-preserving kinetic schemes for the compressible Euler equations[END_REF]. In addition to the preservation of positivity, this procedure also guarantees that solutions satisfy entropy inequalities [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF]. Indeed, the modifications we present in this paper do not affect Berthon's proof about entropy-satisfying solutions.

An important feature of this approach is that the scheme is actually implemented under formulation (2): intermediate variables such as W ij or W * i are never computed. They are only useful from a theoretical point of view. Nevertheless, one must emphasize that CFL conditions (13) may lead to severe computational costs. Indeed, two parameters are not specified in [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF]: the ratio r i corresponding to the localization of nodes m ijp and the position of nodes Z ij which influences the values of |T ij,k |. Arbitrary choices for these parameters may significantly decrease ∆t n as it will be shown in § 6. For Z ij , it is easy to prove that its optimal position corresponds either to a segment or to a single node depending on r i . Indeed, the ratio

|Γ ij,k |/|T ij,k | is equal to 2/h ij,k
where h ij,k is the height associated to Z ij in the triangle T ij,k . The optimization process comes down to maximizing the smallest height in Ω ij . As for r i , it influences all geometric values as well as the eigenvalues

λ n * ij (W * i , W n ij ) (the state W *
i is computed from [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF] where the areas |Ω ij | and |Ω * i | appear). To avoid this dependance on geometric aspects, we present in the next section a more abstract approach based on Berthon's procedure: we keep the idea of expressing the second-order two-dimensional scheme as a convex combination of first-order one-dimensional schemes but we introduce non geometric coefficients in order to optimize CFL conditions.

Algebraic approach

Deriving a new CFL condition

The core of the method is still the convexity of the physical set W. We keep the same notations for the splitting of Ω i as the union of Ω * i and Ω ij even if these elements will not be used a posteriori. Rather than considering the geometry-dependent mean (6), we introduce some coefficients η ij and η * i such that:

W n i = η * i W * i + j∈V(i) η ij W n ij . (14)
We must underline that this very section based on Eq. ( 14) is not restricted to vertex-based approaches. It can adapt to the cell-centered framework directly, the set V(i) consisting of exactly three indices.

From now on, we suppose that coefficients η = {η * i , η ij } are known (see § 5 for more details) and verify:

η ij ≥ 0, η * i ≥ 0, η * i + j∈V(i) η ij = 1. ( 15
)
Instead of (8), we update in time W n ij and W * i by means of first-order schemes:

           W ij = W n ij -∆t n 4 k=1 ζ ij,k F W n ij , W n ij,k , n ij,k , j ∈ V(i), W * i = W * i -∆t n j∈V(i) ζ * ij F W * i , W n ij , n * ij , (16)
where we naturally assume that:

ζ ij,k ≥ 0, ( 17a
) 4 k=1 ζ ij,k n ij,k = 0, (17b)
ζ * ij ≥ 0, ( 17c
)
j∈V(i) ζ * ij n * ij = 0. (17d)
Equalities (17b) and (17d) correspond to the preservation of steady states specific to Eqs. [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. To obtain from ( 2), ( 14) and ( 16) the convex combination for

W n+1 i : W n+1 i = η * i W * i + j∈V(i) η ij W ij ,
a straightforward computation shows that it is necessary and sufficient to have (see Fig. 3):

           η ij ζ ij,1 = η * i ζ * ij , (18a)
η ij ζ ij,2 = η il ζ il,4 , (18b)
η ij ζ ij,3 = |Γ ij | |Ω i | . (18c)
Next step consists in introducing one-dimensional schemes similarly to [START_REF] Cournède | A positive MUSCL scheme for triangulations[END_REF]:

               W ij = 4 k=1 ν ij,k W ij,k , W ij,k = W n ij -∆t n µ ij,k F (W n ij , W n ij,k , n ij,k ) -F (W n ij , W n ij , n ij,k ) ,
and:

     W * i = j∈V(i) ν * ij W * ij , W * ij = W * i -∆t n µ * ij F (W * i , W n ij , n * ij ) -F (W * i , W * i , n * ij ) .
Given (17b) and (17d), the equivalence with ( 16) holds iff:

           4 k=1 ν ij,k = 1, ν ij,k µ ij,k = ζ ij,k , j∈V(i) ν * ij = 1, ν * ij µ * ij = ζ * ij .
Unknowns ν = {ν * ij , ν ij,k } can be easily eliminated:

ν ij,k = ζ ij,k µ ij,k , ν * ij = ζ * ij µ * ij ,
so that hypotheses on ζ and µ reduce to:

             4 k=1 ζ ij,k µ ij,k = 1, ( 19a
)
j∈V(i) ζ * ij µ * ij = 1. (19b)
Consequently, this yields the fact that W n+1 i ∈ W under the assumptions (H1-H2-H3-H4) and the CFL conditions:

       ∆t n max j∈V(i) 1≤k≤4 µ ij,k λ n ij,k (W n ij , W n ij,k ) ≤ α 0 , ∆t n max j∈V(i) µ * ij λ n * ij (W * i , W n ij ) ≤ α 0 . Given that n * ij = -n ij,1 = n ij,3 = n ij , it is convenient to consider the weaker formulation: ∆t n max j∈V(i) µ * ij , max 1≤k≤4 µ ij,k × λn i ≤ α 0 , (20)
λn i := max j∈V(i) 1≤k≤4 |u n ij • n ij,k | + c n ij , |u n ij,k • n ij,k | + c n ij,k .
This CFL condition consists of classical MUSCL eigenvalues together with the ones (for k = 1) associated to the new states W * i and which will be investigated in Sect. 5.

To make this approach legitimate, we aim at minimizing the coefficients µ = {µ * ij , µ ij,k } in order to maximize ∆t n . As we shall see, there is a balance between the order of the method and the CFL condition through the coefficients η ij . To better understand the dependency wrt to these coefficients, we solve the system of constraints.

Solving the constraints

We first notice that constraint (17d) is redundant. Indeed:

j∈V(i) ζ * ij n * ij = - j∈V(i) ζ * ij n ij,1 (18a)
= -

1 η * i η ij ζ ij,1 n ij,1 (17b)
= 1 η * i j∈V(i) η ij [ζ ij,2 n ij,2 + ζ ij,3 n ij,3 + ζ ij,4 n ij,4 ] (18c) = 1 η * i j∈V(i) |Γ ij | |Ω i | n ij =0 (Green) + j∈V(i) η ij ζ ij,2 n ij,2 + j∈V(i) η ij ζ ij,4 n ij,4 (18b)
=

p∈V(i) ηipζip,2(-nip,2) = 0.
Another remark is that (17c) is a consequence of (17a) due to (18a). Moreover, from a linear algebra argument, there are only two independent vectors in the set {n ij,k } 1≤k≤4 . In our case, we choose n ij,3 and n ij,4 as linearly independent vectors. Others are connected by the relations:

4 k=1 |Γ ij,k |n ij,k = 0, n ij,1 = -n ij,3 . (21a)
Due to the fact that m ij,2 G ij,2 G ij,4 m ij,4 is a trapezoidal, the combination of these equalities yields (see Fig. 3 for notations):

n ij,2 = - |Γ ij | M i G ij,2 n ij,3 - M i G ij,4 M i G ij,2 n ij,4 . (21b)
Taking ( 21) into account, (17b) now reads:

0 = -ζ ij,1 - |Γ ij | M i G ij,2 ζ ij,2 + ζ ij,3 n ij,3 + - M i G ij,4 M i G ij,2 ζ ij,2 + ζ ij,4 n ij,4 .
As the vectors are independent, we obtain:

         ζ ij,1 = ζ ij,3 - |Γ ij | M i G ij,2 ζ ij,2 , ζ ij,4 = M i G ij,4 M i G ij,2 ζ ij,2 ,
or equivalently, due to (18c):

η ij ζ ij,2 M i G ij,2 = η ij ζ ij,4 M i G ij,4 = 1 |Ω i | - η ij ζ ij,1 |Γ ij | . ( 22
)
The system of constraints (17a-18a-18b-18c-19a-19b-22) is now solvable. Let

j 0 ∈ V(i) be a fixed index. Set X = ζ ij0,1 . As M i G ij,2 = M i G il,4 = M i G ijl and
due to (18b), we induce that:

1 |Ω i | - η ij ζ ij,1 |Γ ij | = 1 |Ω i | - η il ζ il,1 |Γ il | .
We then iteratively verify that for all j ∈ V(i):

ζ ij,1 = |Γ ij | |Γ ij0 | η ij0 η ij X, ζ * ij = |Γ ij | |Γ ij0 | η ij0 η * i X, ζ ij,3 = |Γ ij | η ij |Ω i | , ζ ij,2 = M i G ijl η ij 1 |Ω i | - η ij0 |Γ ij0 | X , ζ ij,4 = M i G ijk η ij 1 |Ω i | - η ij0 |Γ ij0 | X , with: 0 ≤ X ≤ X max := |Γ ij0 | η ij0 |Ω i | .
The upper bound for X comes from the positivity of ζ ij,2 and ζ ij,4 .

Optimizing the CFL condition

Minimizing the largest µ ij,k satisfying (19a) comes down to taking µ ij,k equal to each other for all k (cf. Lemma 1 in the Appendix), i.e.:

µ ij,k ≡ µ ij (X) = 4 k=1 ζ ij,k = |∂T ij | η ij |Ω i | - X |Γ ij0 | η ij0 η ij |∂T ij | -2|Γ ij | .
Here T ij denotes the triangle M i G ijl G ijk . As for the condition (19b) associated to Ω * i , we similarly take (Lemma 1):

µ * ij ≡ µ * i (X) = j∈V(i) ζ * ij = |∂Ω i | |Γ ij0 | η ij0 η * i X.
The optimization problem thus reduces to:

µ opt i := min 0≤X≤Xmax max µ * i (X), max j∈V(i) µ ij (X) .
As each constraint is linear wrt X and monotone (increasing for µ * i and decreasing for µ ij ), its solution is (see Fig. 4): Given the expressions of µ * i and µ ij with respect to X, this solution reads:

(a). If µ * i (X max ) ≤ max j∈V(i) µ ij (X max ), µ opt i = max j∈V(i) µ ij (X max ); (b). Otherwise, µ opt i = µ * i ( X) where X = max j∈V(i) {X j : µ * i (X j ) = µ ij (X j )}. X µ µ opt i Xmax (a) µ * i (Xmax) ≤ max j∈V(i) µ ij (Xmax) X µ µ opt i X Xmax (b) µ * i (Xmax) > max j∈V(i) µ ij (Xmax)
µ opt i =          2 |Γ ij1 | η ij1 |Ω i | , if |∂Ω i | η * i ≤ 2 |Γ ij1 | η ij1 , |∂Ω i | |Ω i | |∂T ij2 | η ij2 |∂Ω i | + η * i (|∂T ij2 | -2|Γ ij2 |) , otherwise, (23)
where

j 1 = argmax j∈V(i) |Γ ij | η ij and j 2 = argmax j∈V(i) {X j : µ * i (X j ) = µ ij (X j )}.
The optimal choice µ opt i therefore highly depends on the coefficients η ij . The final step would consist in optimizing µ opt i wrt all (η ij ) such that (15) holds. However, these coefficients also influence the reconstruction procedure. That is why it is necessary first to get interested in the computation of reconstructed values.

Reconstruction step

This section is largely independent from the previous one. The overall process in the design of positivity-preserving schemes relies on the convex combination ( 14) and thus on the existence of an intermediate state such that:

W * i = 1 η * i   W n i - j∈V(i) η ij W n ij   ∈ W. (24)
We assumed previously that the reconstructed values W n ij belong to W. In this section, we first explain how to compute these values and secondly how to ensure that W * i ∈ W. Let ∆W n ij be the gradients on each interface of the control volume for the conservative variables, i.e. such that W n ij = W n i +∆W n ij . There are three issues to address when computing ∆W n ij : guaranteeing that W n ij ∈ W, improving the accuracy of the method and avoiding numerical instabilities. It has been proven in the literature that the use of flux/slope limiters may help achieve these purposes [START_REF] Van Leer | Toward the ultimate conservative difference scheme V: A second order sequel of Godunov's methods[END_REF][START_REF] Sweby | High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws[END_REF]. A limiter enables to combine a high-order scheme (where the solution is smooth) and a first-order scheme (elsewhere). Most limiters are functions of consecutive gradients (in one dimension) or upstream/downstream gradients (in higher dimensions) -see Fig. 5.

More precisely, we use the approach described in [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]. Let M i and M l be two nodes (with M l an internal node). We aim at approaching values of a reconstructed variable ξ at the point Q il by means of the formulae:

ξ il = ξ i + α il ϕ(r il )∆ξ il = ξ i + α il φ(r il )∆ξ up il + 1 -φ(r il ) ∆ξ il , (25a)
ξ li = ξ l -α li ϕ(r li )∆ξ il = ξ l -α li φ(r li )∆ξ dn il + 1 -φ(r li ) ∆ξ il . (25b)
The notations are -see Fig. 5:

• α il = M i Q il M i M l , α li = 1 -α il ; • ∆ξ il = ξ l -ξ i = ∇ ξijl • ---→ M i M l where ξijl is a linear approximation of ξ within M i M j M l ; • ∆ξ up il = ∇ ξikm • ---→ M i M l
where ξikm is a linear approximation of ξ within the upstream triangle M i M k M m ; if this triangle does not exist (when M i is located on the boundary), we set ∆ξ up il = 0;

• ∆ξ dn il = ∇ ξlqp • ---→ M i M l
where ξlqp is a linear approximation of ξ within the downstream triangle M l M q M p ;

• r il = ∆ξ up il ∆ξ il and r li = ∆ξ dn il ∆ξ il ; if ∆ξ il = 0, we set r il = r li = 0 so that

ξ il = ξ li = ξ i = ξ l ;
• ϕ is the limiter.

For scalar equations in one dimension, the TVD property requires a slope limiter ϕ to satisfy estimates [START_REF] Sweby | High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws[END_REF] like: 0 ≤ ϕ(r) ≤ min(2r, 2).

To reach higher orders, it is necessary to have ϕ(1) = 1 and ϕ as smooth as possible in the neighbourhood of 1 [START_REF] Spekreijse | Multigrid solution of monotone second-order discretizations of hyperbolic conservation laws[END_REF], which is equivalent to assuming smoothness for φ such that ϕ(r) = 1 -(1 -r)φ(r). The readers may refer to [START_REF] Piperno | Criteria for the design of limiters yielding efficient high resolution TVD schemes[END_REF] and [START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF] about procedures for designing new limiters. For scalar equations in two dimensions, the same inequality turns out to be sufficient for reconstruction at the middle of the edges even if the TVD concept no longer applies [START_REF] Goodman | On the accuracy of stable schemes for 2 D scalar conservation laws[END_REF]. However, for control volumes as on Fig. 1, interfaces do not cross edges at their middles. That is why limiters have to be adapted. Q-limiters have been introduced for cell-centered meshes [START_REF] Clain | L ∞ -stability of the MUSCL methods[END_REF] and τ -limiters for vertex-based triangulations [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]. These limiters satisfy:

0 ≤ ϕ(r) ≤ min(τ r, τ ) ( 26
)
where τ is a geometric parameter to be specified. To guarantee that ξ il > 0 and ξ li > 0 for a limiter ϕ satisfying [START_REF] Piperno | Criteria for the design of limiters yielding efficient high resolution TVD schemes[END_REF], it is necessary to take:

τ = min i,l 1 α il . (27)
Due to Hyp. H0, Q il lies between M i and M l which ensures that α il ∈ [0, 1] and it can be shown [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF] that τ ∈ [START_REF] Aubert | Boundary conditions for Euler equations[END_REF][START_REF] Barth | Finite volume methods: foundation and analysis[END_REF].

To extend this approach to systems of equations, the first question that arises is about the set of variables to which to apply the limitation procedure. Indeed for the Euler equations, it is possible to limit the conservative variables (ρ, ρu, ρv, ρE) [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF], the physical variables (ρ, u, v, p) [START_REF] Clain | A multislope MUSCL method on unstructured meshes applied to compressible Euler equations for axisymmetric swirling flows[END_REF], the characteristic variables [START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF] or the entropic variables [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF]. Berthon carried out numerical simulations in [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF] to highlight the influence of this choice. The fact still remains that there is no compelling argument to make the decision. However as far as this study is concerned, it seems more relevant to limit the physical variables insofar as we are designing a method to ensure that two physical variables (ρ, p) remain positive. Denoting by U = (ρ, u, v, p), the reconstruction is computed like this: Algorithm 1 1. Perform the change of variable U i = κ(W i ); 2. Compute "physical" gradients ∆U ij by applying [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] to each component of

U i ; 3. Compute "conservative" gradients ∆W ij = κ -1 U i + ∆U ij -W i .
By construction, W ij := W i + ∆W ij lies in W.

Construction of W * i

We then investigate the existence of W * i belonging to W. Due to the facts that W is open and that W n i ∈ W, there exist η ij small enough such that (24) holds. However, according to § 3.3 (see [START_REF] Parent | Positivity-preserving high-resolution schemes for systems of conservation laws[END_REF] for example), the smaller η ij , the worse the CFL condition. That is why another strategy has to be found out. Let us introduce β i ∈ [0, 1] and set:

W ij := W i + β i ∆W ij . (28)
The equality 15) and ( 28), Eq. ( 24) then reads:

W ij = (1 -β i )W i + β i W ij shows that W ij ∈ W. Due to (
W * i = W n i -β n i j∈V(i) η ij η * i ∆W n ij = W n i + 1 -η * i η * i β n i ∆W * i , ∆W * i := -1 k∈V(i) η ik j∈V(i) η ij ∆W ij . ( 29
)
∆W * i is (up to the minus sign) an average of the gradients over the cell Ω i . As above, the facts that W n i ∈ W and that W is open imply that there exists β n i small enough so that W * i ∈ W. The case β n i = 1 corresponds to the second order (except for extrema, where the gradient is zero and the scheme degenerates to order 1 owing to the limiter) while β n i = 0 corresponds to the first order. β n i < 1 thus implies a loss of accuracy compared to the pure secondorder MUSCL scheme but turns out to be necessary for guaranteeing positivity. As it will be highlighted in the numerical simulations, β n i = 1 only locally. The requirements ρ * i > 0 and p * i > 0 corresponding to W * i ∈ W reduce to: 3

P 1 1 -η * i η * i β n i := 1 + D n i 1 -η * i η * i β n i > 0, ( 30a
)
P 2 1 -η * i η * i β n i := 1 + B n i 1 -η * i η * i β n i + A n i 1 -η * i η * i β n i 2 > 0, ( 30b
)
3 P 1 and P 2 are such that

ρ * i = ρ n i P 1 1-η * i η * i β n i and p * i = p n i P 2 P 1 1-η * i η * i β n i .
with:

A n i = γ -1 ρ n i p n i ∆ρ * i ∆(ρE) * i - |∆(ρu) * i | 2
2 ;

B n i = γ -1 p n i E n i ∆ρ * i + ∆(ρE) * i -u n i • ∆(ρu) * i ; D n i = ∆ρ * i ρ n i .
Consider the roots of the linear inequality (30a) and of the (at most) quadratic constraint (30b). We must bear in mind that our aim is to take β n i as large as possible in [0, 1] not to damage accuracy. For the first inequality, we set:

ϑ (ρ) i :=      +∞, if D n i ≥ 0, - 1 D n i , otherwise, and
β (ρ) i := min 1, η * i 1 -η * i ϑ (ρ) i . (31a)
As for the other inequality, we set

δ n i = (B n i ) 2 -4A n i . This term is clearly positive: 4 δ n i = γ -1 p n i 2 ∆(ρE) * i -u n i • ∆(ρu) * i + ∆ρ * i |u n i | 2 -E n i 2 + 2 γ -1 p n i ρ n i ∆ρ * i |u n i | - u n i • ∆(ρu) * i |u n i | 2 + ∆(ρu) * i 2 - [u n i • ∆(ρu) * i ] 2 |u n i | 2 .
Hence we set:

ϑ (p) i :=                +∞, if (A n i = 0, B n i ≥ 0) or (A n i > 0, B n i > 0), - 1 B n i , if (A n i = 0, B n i < 0), - B n i + δ n i 2A n i , otherwise, (31b) and β
(p) i := min 1, η * i 1 -η * i ϑ (p) i
. To satisfy (30a) and (30b) simultaneously, we then have to take:

β n i ≤ min β (ρ) i , β (p) i . ( 32
)
The analysis of this result and the actual choice for β n i are specified in the next section. 4 This expression holds if u n i = 0.

Otherwise,

p n i = (γ -1)ρ n i E n i and δ n i = E n i ∆ρ * i -∆(ρE) * i 2 + 2E n i ∆(ρu) * i 2 /(ρ n i E n i ) 2 > 0. 20

Practical issues

The final step of our approach consists in choosing coefficients η ij so as to obtain to a balance between the CFL condition (characterized by λn i and µ opt i ) and the order of the scheme (characterized by β n i ). A first comment is that it is always possible to introduce ηij such that:

η ij = (1 -η * i )η ij and j∈V(i) ηij = 1 =⇒ ∆W * i = - j∈V(i) ηij ∆W ij .
Note that ∆W * i is independent from η * i . The optimal choice would consist in maximizing ∆t n to avoid extra computations (equivalent to minimizing λn i and µ opt i ) and maximizing β n i to improve accuracy. However, as it will be demonstrated in the sequel, these two goals are antagonistic. Moreover, due to numerous nonlinearities, it does not seem achievable to state an optimal choice for the three aforementioned coefficients. That is why we choose to mainly focus on µ opt i . Indeed, the worst case for β n i corresponds to β n i = 0 and the scheme locally degenerates to order 1. An unrelevant choice for η * i and ηij would make µ opt i go to ∞ and ∆t n to 0. Hence the priority given to µ opt i involved in the CFL condition even if we first study the influence of η * i on the two other coefficients (parameters ηij are assumed to be given).

In Sect. 4, we detailed x how to obtain reconstructed values in W by means of τ -limiters and y how to ensure W * i ∈ W thanks to the use of a damping coefficient β n i . The latter parameter must satisfy the constraint [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]. A first comment is that if

β (ρ) i < 1 (resp. β (p) i < 1), we cannot take β n i = β (ρ) i (resp. β n i = β (p)
i ) since this would imply ρ * i = 0 (resp. p * i = 0) and W * i ∈ W. Secondly, it is important to mention where W * i is involved from a practical point of view. Although this additional variable is used to prove positivity, it is not part of the numerical scheme [START_REF] Barth | Finite volume methods: foundation and analysis[END_REF]. It only determines the CFL condition (20) through the

term |u n ij,1 • n ij,1 | + c n ij,1 = (ρu) * i ρ * i • n ij + c * i . To prevent c * i = γp * i /ρ * i (

and so λn

i ) from growing drastically, it seems better to have

β n i close to β (p) i than to β (ρ)
i . More precisely, we finally take:

β n i := min 1, η * i 1 -η * i σ ρ ϑ (ρ) i , η * i 1 -η * i σ p ϑ (p) i , (33)
with σ ρ < σ p < 1. Fig. 6 shows how β n i evolves depending on η * i . Hence it is tempting to take η * i large enough as it would ensure β n i = 1. But as we shall see later, large η * i must be avoided. To go back to λn i , we sketch out the evolution of the "starred" eigenvalue

η * i β n i 1 1 1+ϑ i 1 (a) ϑ i < +∞ η * i β n i 1 (b) ϑ i = +∞ Figure 6: Plot of β n i (η * i ) depending on ϑ i := min{σ ρ ϑ (ρ) i , σ p ϑ (p) i } η * i µ opt i η * i 1 (a) ε i < 0 η * i µ opt i η * i 1 (b) ε i ≥ 0 Figure 7: Plot of µ opt i (η * i ) depending on ε i wrt η * i : |u * i • n ij | + c * i = (ρu) * i ρ * i • n ij + c * i = 1 P 1 1-η * i η * i β n i u n i • n ij + 1 -η * i η * i β n i ∆(ρu) * i ρ n i • n ij + c n i P 2 1 -η * i η * i β n i . (34)
This eigenvalue may be predominant in λn i in [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF]. We first notice that other eigenvalues

|u n ij,k • n ij,k | + c n ij,k , k = 1, remain bounded no matter what η * i (they only depend on β n i ∈ [0, 1]
). We then deal with two cases: if

ϑ i := min{σ ρ ϑ (ρ) i , σ p ϑ (p) i } < ∞, then the term 1-η * i η * i β n i belongs to [0, ϑ i ].
The rhs in (34) is thus bounded.

However, when ϑ

(ρ) i = ϑ (p) i = +∞, then β n i = 1, 1/P 1 ( 1-η * i η * i
) ≤ 1 and the rhs in [START_REF] Woodward | The numerical simulations of two-dimensional fluid flow with strong shocks[END_REF] goes to +∞ when η * i → 0. Once more, large η * i seem more suitable. We now investigate the tuning of µ opt i in the CFL condition [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF]. The resulting time step is such that the scheme does preserve positivity of density and pressure. As mentioned above, we expect µ opt i to be as small as possible. The problem reads:

min (ηij )∈R #V(i) + j ηij =1 min 0≤η * i ≤1 µ opt i (η * i , ηij ). (P)
With the η * i -parameterization, (23) reads:

µ opt i (η * i , ηij ) =            2 (1 -η * i )|Ω i | max j∈V(i) |Γ ij | ηij , if η * i ≥ η * i , |∂Ω i | |Ω i | min j∈V(i) η * i 1 - 2|Γ ij | |∂T ij | - ηij |∂Ω i | |∂T ij | + ηij |∂Ω i | |∂T ij | -1 , (35) where:
η * i (η ij ) := |∂Ω i | |∂Ω i | + 2 max j∈V(i) |Γ ij | ηij .
We still denote by j 2 the index for which the minimum is reached in the second case in (35) and we set:

ε i := 1 - 2|Γ ij2 | |∂T ij2 | - ηij2 |∂Ω i | |∂T ij2 |
which is the coefficient in front of η * i in [START_REF] Zhang | On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF]. The profile of µ opt i is pictured on Fig. 7 depending on the sign of ε i .

According to these figures, the minimum is thus reached either for η * i = 0 or

η * i = η * i .
In view of what was shown above, the former case must not be chosen as λn i may tend to +∞. Therefore, a reasonable choice seems to be η * i = η * i . In that case, µ opt i reads:

µ opt i = |∂Ω i | + 2 max j∈V(i) |Γ ij | ηij |Ω i | . (36)
Hence: (P) ≤ min

(ηij )∈R #V(i) + j ηij =1 µ opt i η * i (η ij ), ηij .
Applying Lemma 1, we obtain the exact solution for the right optimization problem which is:

ηij = |Γ ij | |∂Ω i | , η * i = 1 3 , µ opt i = 3 |∂Ω i | |Ω i | , (37a)
β n i = min 1, σ ρ 2 ϑ (ρ) i , σ p 2 ϑ (p) i . (37b)
The most striking feature of this approach is that η * i is independent from the cell Ω i and this enables to save computational time. Moreover, it is far away from 0 to avoid the issues mentioned previously about β n i and λn i . Eventually, numerical simulations go to show that η * i is large enough to provide order 2 almost everywhere (see § 6).

To conclude, we should underline that coefficients (37a) are optimized from the point of view of the CFL conditions. Other choices are possible depending on the feature one wants to focus on (accuracy, computational time, . . . ): these coefficients remain user-tuned.

We see on Fig. 8 the influence of this choice on µ opt i and thus on the computation of the time step through [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF] for a single control volume. Our choice (37a) pictured in dashed red line on Fig. 8 provides a significantly small lower bound for the CFL condition compared to the other coefficients implemented:

Case 1. ηij = 1 #V(i) ; Case 2. ηij = |∂T ij | |Ω i | ; Case 3. ηij = |Γ ij | |∂Ω i | for various η * i .
Case 2. corresponds to coefficients used by Berthon [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF] (η * i = r 2 i ). Although η * i is finally tuned by the user, we notice that an arbitrary choice may have dramatic consequences on the computational time.

We end up this section with a practical algorithm. Given a numerical code solving the Euler equations, here are the few modifications to make: More precisely, there is an initial step that replaces the computation of the smallest height. The very difference with classical codes consists in computing steps 2.(b) and 2.(c). They correspond to an additional loop but nevertheless do not produce prohibitive extra computational cost. It is the price to pay to ensure that no matter how close to vaccum the physical solution may be, the numerical solution remains admissible.

Numerical results

This section deals with numerical simulations on benchmarks in order to highlight the effects of the method. Our approach is very general and adapts to many numerical fluxes provided they satisfy (H1)-(H2)-(H3) which are very classical and (H4). This is the latter requirement which may fail. The most Figure 9: Profile of the pressure at time t famous (and maybe the most commonly used) scheme that does not work is the Roe scheme [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF]. This defect has been corrected in the HLL family (see for instance the HLLE scheme [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF]). Entropy fixes have also been derived to avoid non-physical (i.e. not entropy-satisfying) solutions. Some of them may be reinterpreted as positive schemes -see e.g. [START_REF] Pelanti | A review of entropy fixes as applied to Roe's linearization[END_REF].

Other schemes that do satisfy (H4) are for instance the Godunov method, the Lax-Friedrichs flux (and consequently the Rusanov flux) -see [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF]Appendix] or [5, § 2.4.2] -and the Siliciu relaxation scheme -see [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF] or [5, § 2.4.4].

We must underline that in many situations (H4) is not required if the numerical state is far enough from the critical region (ρ = 0, p = 0) and if gradients are not too steep. That is why Roe scheme provides admissible results in many cases. What we state in this paper is that our approach including reconstruction (with τ -limiters and β n i ) and relevant CFL conditions guarantees the positivity of density and pression no matter what the data.

We present in the sequel some numerical simulations with the Rusanov flux. Despite diffusive effects, it yields accurate results. As far as the limiter is concerned, we use a modified Van Leer limiter [START_REF] Clain | L ∞ -stability of the MUSCL methods[END_REF] so that (26) holds:

ϕ(r) =              0, if r < 0, τ r 1 + (τ -1)r , if 0 ≤ r ≤ 1, τ r τ -1 + r , if r > 1,
with τ defined by [START_REF] Serre | Systems of conservation laws. 2. Geometric structures, oscillations, and initial-boundary value problems[END_REF]. See Tab. 1 for concrete values.

123 problem

The first problem is a one-dimensional Riemann problem for the Euler equations (3) which is known to be one of the most suitable tests to assess the robustness of a scheme. Indeed, the so-called 123 problem consists of two rarefaction waves where the intermediate state U # is close to vacuum (ρ 1, p 1). It is presented in [32, § 6.4]. The left and right initial physical states are U l = (ρ 0 , -u 0 , 0, p 0 ) and U r = (ρ 0 , u 0 , 0, p 0 ). The profile of the pressure at time t is pictured on Fig. 9 (the density has a similar profile).

Mesh

ρ # = ρ 0 1 - γ -1 2
u 0 c 0 2 γ-1 , u # = 0, v # = 0, p # = p 0 1 - γ -1 2
u 0 c 0 2γ γ-1
.

For instance, for γ = 1.4 and the original data (ρ 0 = 1, u 0 = 2, p 0 = 0.4), then the initial value in conservative variables is W r = (1, 2, 0, 3) and:

ρ # ≈ 0.0219 and p # ≈ 0.0019.

But it is possible to tune the initial data to get closer to vacuum. In particular, the critical set of data is given by u 0 = 2c0 γ-1 for which ρ # = 0 and p # = 0. To prove the robustness of the procedure, we take:

ρ 0 = 1, p 0 = 0.4, u 0 = 3.74 < 2c 0 γ -1 ≈ 3.742.
The corresponding intermediate states are about the machine epsilon. The exact density remains positive and so must do the numerical solution under the CFL condition [START_REF] Linde | Robust Euler codes, 13th AIAA CFD Conf[END_REF]. This 1D problem is simulated with a 2D code. We consider the 2D domain Ω = [-0.5, 0.5] × [0, 0.25] for which two meshes are considered: a cartesian grid (Fig. 10a -8180 nodes, 15920 triangles) and an unstructured grid with a straight interface generated by Triangle Mesh Generator [START_REF] Shewchuk | Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator[END_REF] (Fig. 10b -8064 nodes, 15799 triangles).

Geometric parameters involved in CFL conditions (smallest height over the whole tesselation for order 1, (µ opt i ) -1 for order 2) are specified in Tab. 1. To ensure positivity of density and pressure, the present study requires (according to the figures) to divide the time step by about 10 (wrt the 1st-order CFL condition) which is equivalent to processing 10 times more iterations. This is in accordance with a similar study devoted to the scalar case [7, Eq. ( 40)].

The Riemann problem corresponds to the initial condition:

W 0 (x, y) = W l , if x < 0, W r , if x > 0.
All the results in the sequel are presented over the domain [0, 0.5] × {0.125} insofar as the solution is symmetric (wrt to x = 0) and invariant wrt y.

The two specific issues related to this test are about the positivity of density and pressure and the profile of energy (cf. [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]Fig. 6.14] or [START_REF] Clain | A multislope MUSCL method on unstructured meshes applied to compressible Euler equations for axisymmetric swirling flows[END_REF]Fig. 14] where there is a noticeable overshoot in the vicinity of 0). Typically, the 2nd-order scheme (2) without the β-procedure fails to preserve the positivity of p which becomes negative at iteration 13 on both kinds of grid no matter what the time step. The use of the coefficient β n i (and of the additional state W * i ) is necessary to ensure that W n+1 i ∈ W as expected. In many non critical cases, the classical version of the FV scheme (2) does not suffer from positivity defects. However, this case proves that the additional numerical tools we present in this paper help to cure this issue. Coefficient β n i may be interpreted as a new class of limiters devoted to the preservation of positivity.

More precisely, this coefficient is "activated" (in the sense that β n i < 1) only for a few iterations (up to iteration 19 for the structured mesh, up to iteration 110 for the unstructured grid) and in a narrow area (see Fig. 11 for the structured mesh) consisted of the vertices which belong to the interface. It can be accounted for by the steep gradients on the density variable -this can be inferred from (31a) that the greater ∆ρ * i , the lower β n i . A similar analysis may be carried out for the pressure variable but expressions for β (p) i is much more intricate.

Fig. 12 shows the comparison (for the density variable) between the 1storder scheme together with the standard CFL condition

∆t n max i (|u n i | + c n i ) ≤ α 0 min i h i
and the modified 2nd-order scheme (i.e. including the β n i coefficient) -referred to as "order β" in the legend -processed on the structured grid. The latter scheme provides more accurate results as expected.

Results at a larger time (β n i is no longer activated) are depicted on Fig. 13 for the energy and on Figs. 14 for the density. On Fig. 13, we notice that energy is correctly handled (no overshoot) in the vicinity of 0 no matter what the order of the scheme. However, results are not accurate due to the fact that computing E # involves the ratio p # /ρ # while these two variables are close to 0.

Fig. 14a shows that 1st-and β-order versions of the scheme satisfy the positivity constraint (unlike the 2nd-order scheme which fails as mentioned above). The structured solutions obviously seem better than their unstructured counterparts. This clear advantage seems in accordance with the fact that this test is one-dimensional in essence. Indeed, if we pay attention to two specific areas (namely where the solution is less smooth at the origin of the rarefaction wave -Fig. 14b and where it gets close to 0 -Fig. 14c), we remark different qualitative evolutions. In particular, in the vicinity of 0, order β (unstr) yields less accurate results than order 1 (unstr). The conclusion is reversed in the structured case. Moreover, although β n i = 1 everywhere (corresponding to order 2), solutions are close to order 1. The first iterations are thus predominant in this test.

Mach 3 wind tunnel case

The second problem is a genuine 2D case presented in [START_REF] Woodward | The numerical simulations of two-dimensional fluid flow with strong shocks[END_REF]. It consists in an inward Mach 3 flow within a singular domain. More precisely, the domain is a rectangle which is 3 meter long and 1 meter high with a step (2.4m × [START_REF] Ch | A finite volume method with efficient limiters for solving conservation laws[END_REF][START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF][START_REF] Woodward | The numerical simulations of two-dimensional fluid flow with strong shocks[END_REF] for instance) the solution may consist of a front wave before the step and then two reflections. From a numerical point of view, two main difficulties are about the handling of the singularity (the edge of the step) and the location of reflections.

Two unstructured meshes are considered (resp. made of 20291 -Mesh 1 -and 80441 nodes -Mesh 2). Using the 1st-order scheme (Fig. 15a) does not enable to locate accurately the first reflection (the second one is not even detectable) and the front wave is blurred. The 2nd-order scheme provides better results (Fig. 15b). Indeed, we recover the two reflections and the front wave is clearer. A striking point is that in this simulation, β n i is constantly equal to 1 which means that the pure 2nd-order MUSCL scheme is able to tackle this problem. Decreasing the value of η * i (which is the scalar parameter in § 5 and set to 1/3 in all computations) to 0.1 leads to the activation of β n i in the vicinity of the front wave and close to the corner. This remark goes to show that the choice η * i = 1/3 not only improves the computational time, but also preserves accuracy. Furthermore, we see on Fig. 15c that the addition of numerical tools (coefficient β n i in the reconstruction, modification of the CFL condition) does not damage the convergence of the scheme: when refining the mesh, we obtain more accurate results with clearer waves.

Conclusion

We carried out a theoretical study of whether MUSCL schemes for the Euler equations ensure positivity of density and pressure. This is a generalization of seminal works from Perthame & Shu and Berthon. On the one hand, we classicaly rewrote a standard MUSCL scheme as a convex combination of 1st-order 1D schemes and we introduced abstract coefficients allowing for an optimization of the CFL condition (which is thus explicit). This process assumes properties of the numerical flux which are satisfied by many classical schemes but not Roe scheme. It heavily relies on the 2D geometry (and could be extended to 3D) since it consists of one more step compared to dimension 1. On the other hand, to take advantage of an additional state as detailled in Berthon [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF], we compute a damping coefficient β n i whose role consists in maintaining updated values in the set of physically admissible states. We gave simple and practical directions to modify existing codes in order to guarantee the robustness of the scheme.

Choices have been made to tune the parameters we introduced in this study in order to find out a balance between accuracy and computational costs. Although we favoured the latter one through the optimization of the CFL condition, our choice turned out to be relevant about accuracy as well. To be more precise about the CFL condition, it is practically about 10 times more constraining than standard ones in accordance with a similar scalar study [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]. We should however bear in mind that our results come from sufficient conditions. Simulations have been performed with the Rusanov flux which is known to be diffusive. As our approach can adapt to any scheme, we may obtain more striking results with less diffusive fluxes. Numerical simulations also showed that the coefficient β n i is particularly useful when the solution is close to the boundary of the physical set together with steep gradients.

This study was devoted to the Euler equations for ideal gas. A first extension would include the handling of other equations of state. This would only modify the computation of β n i . Future work might also deal with other physical models. The process first requires to identify the physical constraints (provided the resulting set is convex) and then to adapt the computation of β n i and µ opt i .

Appendix Lemma 1

Let ϑ ∈ R d with positive components, ϑ = 0. The optimization problem:

X ϑ := max χ∈E ϑ min 1≤j≤d χ j with E ϑ = {χ ∈ R d , χ j ≥ 0, χ • ϑ = 1} is solvable with X ϑ = |ϑ| -1
1 and a maximizer is given by χ j = |ϑ| -1

1 for all j.

Proof. For any χ ∈ E ϑ , we observe that:

1 = χ • ϑ ≥ |ϑ| 1 min 1≤j≤d χ j .
Thus X ϑ ≤ |ϑ| 1 and the maximum is reached fo the constant vector χ.

Figure 1 :

 1 Figure 1: Vertex-based mesh: control volumes Ω i and Ω l

Figure 2 :

 2 Figure 2: Sub-cells within a control volume

Figure 3 :

 3 Figure 3: Splitting of sub-cells

Figure 4 :

 4 Figure 4: Graphs of µ * i (dashed blue line) and µ ij (plain red lines)

Figure 5 :

 5 Figure 5: Reconstruction step: upstream and downstream gradients

Figure 8 : 2 1.

 82 Figure 8: Plot of µ opt i (η * i ) (right) for a single control volume (left)

Figure 10 :

 10 Figure 10: Grids for 123 problem

Figure 11 :

 11 Figure 11: Location of vertices where β n i = 1

Figure 12 :

 12 Figure 12: Density plot at time 0.0125: 1st-order and modified 2nd-order schemes on a structured grid

Figure 13 :

 13 Figure 13: Energy plot at time 0.075

Figure 14 :

 14 Figure 14: Density plot at time 0.075

Figure 15 :

 15 Figure 15: Density plot at time 4.0: 30 contours from 0.7 to 6.5

- 1 1. 1 1. 1 1.

 111 Let χ be defined by χ j = |ϑ| -The vector χ obviously satisfies χ ∈ E ϑ and min 1≤j≤d χ j = |ϑ| -Hence we conclude that X ϑ = |ϑ| -1

Table 1 :

 1 Coefficients involved in CFL conditions and in reconstruction stepThe intermediate state U # is given by:

		min h k	1/µ opt i	τ
	Structured	0.0025 0.00026 1.50
	Unstructured 0.0022 0.00017 1.24

Monoslope means that a unique gradient is computed within the control volume. This gradient is used to reconstruct values at all interfaces. Multisplope methods use one gradient per interface.

For a numerical flux which would not satisfy the rotational invariance, (H4.CFL) is replaced by a genuinely multidimensional version.

Acknowledgements

The authors are grateful to the referees for their helpful comments and suggestions which made the reading of technical parts more pleasant.