Stéphane Del Pino
email: stephane.delpino@cea.fr

Isabelle Marmajou
email: isabelle.marmajou@cea.fr

Triangular metric-based mesh adaptation for compressible multi-material flows in semi-Lagrangian coordinates

Keywords: Gas dynamics, Lagrange, finite-volume, adaptive mesh refinement, Maximum Principle. mixtures

In this paper, we propose an adaptive mesh refinement method for 2D multi-material compressible nonviscous flows in semi-Lagrangian coordinates.

The mesh adaptation procedure is local and relies on discrete metric field evaluation. The remapping method is second-order accurate and we prove its stability. We propose a multi-material treatment using two ingredients: the local remeshing is performed in a way to reduce as much as possible mixture creation and an interface reconstruction method is used to avoid material diffusion in mixing cells.

The obtained method is almost Lagrangian and can be implemented in a parallel framework. We provide some numerical tests which attest the validity of the method and its robustness.

Introduction

The aim of this paper is to approximate multi-component compressible gas dynamics in semi-Lagrangian coordinates with mesh adaptation. This work is an extension of [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF] to the second-order of accuracy and to the multi-material case, including mixing.

In this paper, we consider the simple case of a mono-fluid model for multi-material problems: all constituents share the same velocity and temperature. This kind of model is suitable to deal with gas or plasma at the fluid velocity u,

d dt ω(t) 1 = ∂ω(t) u • n, (1)
∀1 ≤ i ≤ N, d dt ω(t) ρ i = 0, (2)
ρ = N i=1 ρ i , (3)
d dt ω(t) ρu + ∂ω(t) pn = 0, (4)
d dt ω(t) ρE + ∂ω(t) pu • n = 0. (5)
In ( 1)-( 5), u denotes the fluid velocity, ∀1 ≤ i ≤ N, ρ i is the partial density of material i and ρ refers to the density of the fluid, E is the total energy density: E = e +1 2 u2 , where e denotes the internal energy density of the mixture. Finally p is the pressure of the fluid, which is given by an equation of state: p = p(ρ 1 , . . . , ρN , T ). ∀1 ≤ i ≤ N , ρi is the density of material i and denoting f i its volume fraction: one has f i ρi = ρ i .

We do not detail the mixture equation of state itself, one generally use a Dalton law [START_REF] Gillespie | The Gibbs-Dalton Law of Partial Pressures[END_REF]. In the following we simply assume that the pressure law is such that the system is hyperbolic. Also, we do not describe the evolution of the (f i ) 1≤i≤N , since it is not the purpose of this paper.

Injecting the mass conservation of each material (2) into (3), one trivially retrieves the mass conservation of the fluid d dt ω(t) ρ = 0.

The mass fraction of material i, c i := ρi ρ is just transported by the fluid

∀1 ≤ i ≤ N, D t c i = 0, (7)
where D t := ∂ t + u • ∇ is the Lagrangian derivative.

Thus to solve (1)-( 5) one can equivalently consider equations ( 2)-( 3) or ( 6)- [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF]. In the later case, one retrieves the classical integral formulation of compressible Euler equations, where the mass fractions are just transported 1 .

In this paper, we use a finite-volume semi-Lagrangian scheme (Glace [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF] or Eucclhyd [START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF]) to approach the solutions of (1)- [START_REF] Anderson | An arbitrary Lagrangian-Eulerian method with adaptive mesh refinement for the solution of the Euler equations[END_REF].

Classically, the use of semi-Lagrangian methods 2 to approximate fluid flows can lead to distorted meshes:

especially in case of flows containing vortexes or shears. If the mesh is too distorted, the numerical approximation of the solutions of (1)-( 5) can be of very bad quality, leading eventually to the calculation failure.

To overcome this classical issue, a standard approach consists into replacing the mesh by a new one before it becomes improper to calculation. After such a change, the numerical solution is remapped on the new mesh and one can use the semi-Lagrangian hydro-scheme to proceed with the calculation 3 .

Generally, it is enough to define a new mesh by simply displacing the vertices of the cells to improve their quality. In that case, the new mesh and the old mesh share the same connectivity. This approach is usually referred as indirect ALE (Arbitrary Lagrangian Eulerian), see [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF] for instance.

The indirect ALE approach is very successful and is used in many hydro-codes. However, ALE has some intrinsic limitations since the connectivity of the mesh does not change. For instance, in the case of dominant vortexes, ALE can degenerate to a quasi-Eulerian method (see the introduction of [START_REF] Loubère | ReALE: A Reconnection-based Arbitrary-Lagrangian-Eulerian Method[END_REF] where this phenomenon is fully described). If this behavior does not invalidate the indirect ALE approach, it reduces greatly the advantages of the Lagrange plus remap strategy.

Obviously allowing changes of connectivity during the remeshing step has more potential than just smoothing the mesh, but it is more complex by nature:

• one has to provide a systematic way to compute the new mesh,

• remapping onto the new mesh can be much more difficult to handle, especially in the case of parallel computations.

This kind of methods (see for instance [START_REF] Fritts | The Free-Lagrange Method[END_REF][START_REF]Advances in the Free-Lagrange Method[END_REF]) have regain popularity recently. For instance, in [START_REF] Anderson | An arbitrary Lagrangian-Eulerian method with adaptive mesh refinement for the solution of the Euler equations[END_REF] a patch-based adaptive mesh refinement method was proposed in the indirect ALE context.

In [START_REF] Loubère | ReALE: A Reconnection-based Arbitrary-Lagrangian-Eulerian Method[END_REF], the ReALE method is proposed. In a few words, it consists in maintaining a Voronoï tessellation all along the calculation. In [START_REF] Loubère | ReALE: A Reconnection-based Arbitrary-Lagrangian-Eulerian Method[END_REF], the number of cells is constant, but it was made variable in [START_REF] Bo | Adaptive reconnection-based arbitrary lagrangian eulerian method[END_REF] to allow adaptivity.

Another approach has been proposed in [START_REF] Hoch | Semi-conformal polygonal mesh adaptation seen as grid velocity formulation for ale simulations[END_REF]. The idea is to define a way to deal with polygonal mesh adaptation, the first results are promising but it is still a work in progress.

Obviously, alongside with Voronoï based methods are their dual approaches: triangular-based (Delaunay) remeshing methods. An example of their use in the Lagrangian hydrodynamics context is given in [START_REF] Lin | A local rezoning and remapping method for unstructured mesh[END_REF]. In this study, Authors replace large portions of the mesh by new triangulations when the quality of the mesh deteriorates. The main drawback of this method is that one must perform a non-local remapping from one mesh to the other. Even in the case of triangles, it remains a complex task, especially for parallel computations.

Finally, in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF] we proposed a triangular-based local remeshing method. It relies on a geometric error estimate to define adaptation criteria (see for instance [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF][START_REF] Alauzet | Estimation d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage[END_REF][START_REF] Alauzet | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Dobrzynski | Fast and accurate simulations of air-cooled structures[END_REF][START_REF] Dobrzynski | Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments[END_REF]). The aim of this method is not only to maintain good quality of the mesh to use Lagrangian schemes but also to adapt the mesh to improve the quality of the numerical approximation. In other words this method could be defined as a conforming AMR (Adaptive Mesh Refinement) technique in the Lagrangian hydrodynamics context. The adaptation criterion follows [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF] using the very simple local adaptation patterns depicted in [START_REF] Dobrzynski | Fast and accurate simulations of air-cooled structures[END_REF]. It has been successfully extended to 3D in [START_REF] Hansen | Alexa: MPI+X Shock Hydrodynamics on Dynamically Adaptative Tetrahedral Meshes[END_REF][START_REF] Ibanez | multi-material dynamic domain topoplogy changes in the Lagrangian Grid Reconnection (LGR) code[END_REF] in the case of Euclidean metrics.

Let us emphasize that mesh adaptation is completely local: the new quantities and discrete metric fields are updated after each local mesh modification. The proposed mesh adaptation procedure is convergent so that the mesh remains adapted all along the simulation. As a consequence, very few mesh modifications are required from one time step to the other.

In this paper, extending [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF], the multi-material treatment is improved in two ways. In the one hand, the mesh adaptation strategy is changed to take into account mixing flows and in the other hand an interface reconstruction procedure is proposed to reduce numerical material diffusion.

Also, the remapping procedure is extended to second-order while ensuring stability by means of Maximum Principles [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF].

Finally, let us recall that the use of local mesh adaptation techniques used in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF] has two more advantages:

local adaptation makes the method cheap (the adaptation criteria do not change a lot from one time step to the other) and makes it possible to handle parallelism 4 .

The remaining of the paper articulates as follows. In the next part, we describe the mesh adaptation procedure itself. Then, we discuss a second-order conservative remapping strategy which satisfies Maximum Principles: a key to stability. Next, we present some numerical tests to assess the efficiency of the method.

Finally, we conclude and draw some potential extensions of the present work. Some proofs are given in the Appendix to maintain a fluent reading of the manuscript.

Notations

We define the following notations:

• M denotes a mesh,

• J is the set of all cells of M,

• L is the set of all edges of M,

• and R is the set of all vertices of M.

Also, we define the following subsets:

• ∀j ∈ J , R j is the set of vertices of the cell j,

• ∀l ∈ L, R l is the set of vertices of the edge l.

Triangular mesh adaptation

In this section, we present the triangular mesh adaptation strategy that we developed to deal with multimaterial flows. One of the key aspects of the method resides in the Lagrangian context itself. Actually, the method is built to take advantages of the finite-volume semi-Lagrangian scheme used to approximate the fluid flow [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique Lagrangienne bidimensionelle[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF].

In our view, if the semi-Lagrangian context adds some constraints, it gives more coherence to the whole method (Lagrangian-AMR). Actually, semi-Lagrangian schemes provide some kind of "automatic mesh adaptation" with regard to the flow: contact discontinuities are preserved exactly, the mesh is naturally refined in shocks, the numerical solution is just displaced in pure transport regions, etc. . . As a consequence, the AMR approach must try not to break these features, but take advantage of them.

These constrains were successfully satisfied in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF]: for instance, mesh was still refined into shocks, mesh was not modified in transport regions and interface were preserved since no mixing was allowed.

However, in order to prevent the creation of mixing, we arbitrarily canceled material mass fluxes when coarsening the mesh at material interfaces. This has the unpleasant effect of reducing the stability of the method since the new density violates the Maximum Principle in favor of the conservation of the mass of each material. To be more precise, as soon as mixing creation is forbidden, one has to make the choice between mass conservation and stability (Maximum Principle). Also, preventing mixing may lead to very small cells which can impose arbitrary small time steps.

In this paper we explore the second alternative (which is probably more natural): we allow creation of mixing in some circumstances and then satisfy mass conservation for each material. We shall see in Section 2 that doing so, the stability of the method can be proven.

We now present the remeshing strategy that we developed in this context. One shall see in Section 1.5 how the mesh is changed to limit the creation of mixing.

Remeshing constrains

As stated before, the mesh modifications must satisfy a few constrains.

First of all, we want the whole method to be as Lagrangian as possible. In other words if the mesh M n is adapted to the approximation of the solution at some time t n , we want the mesh M n+1 to be as close as possible to M n after the Lagrangian motion at time t n+1 . Being "as Lagrangian as possible" has two direct consequences. In the one hand, mesh adaptation does not deteriorate the good properties of the semi-Lagrangian scheme. In the other hand, in practice, the less the mesh will be modified, the less expensive the method will be.

A second practical constrain is that this method must integrate a massively parallel multi-physics code [START_REF] Grospellier | The arcane development framework[END_REF][START_REF] Lefebvre | Development and validation of the TROLL radiationhydrodynamics code for 3D hohlraum calculations[END_REF]. In this context, it is highly desirable to build parallel algorithms that give exactly the same result regardless of the number of processors.

All these requirements indicate that the mesh changes should be as local as possible. We rely on [START_REF] Dobrzynski | Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments[END_REF][START_REF] Dobrzynski | Fast and accurate simulations of air-cooled structures[END_REF] where the Authors define a sufficient set of local mesh modifications patterns to achieve adaptivity. These local modifications patterns are recalled in the following paragraph.

1.2. Mesh modification patterns: how to modify the mesh?

In the sake of completeness, we recall the patterns (that we already described in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF]). All three patterns are edge-oriented: they are defined by means of edge modifications. Each pattern is dedicated to a specific treatment. The swapping pattern is defined to improve the local quality of the mesh. The splitting pattern allows to refine the mesh and the collapsing pattern permits local coarsening of the mesh.

The patterns are defined in such a way that remeshing conserves locally the volume of the mesh. As a consequence, Geometrical Conservation Law (GCL) is satisfied by construction.

We define a cavity (or a patch) as the smallest set of triangles involved in a local mesh modification.

edge swapping For this operation, the cavity is the set of the two adjacent triangles to a given edge.

Swapping consists in considering the quadrangle formed by the two initial triangles and to replace them by the two new triangles built on its other diagonal, see Figure 1. Obviously, the swapping pattern is only defined on inner edges. edge splitting This pattern corresponds to the local refinement of the mesh. It consists in subdividing the old cavity triangles by joining their opposite vertex to the middle of the considered edge, see Figure 2.

edge collapsing This pattern is used to provide an edge-based mesh coarsening. As depicted on Figure 3, removing an edge can be viewed as collapsing it by displacing one of its vertices to the other. The corresponding initial cavity is thus the set of all the cells connected to the "displaced" vertex. Figure 3: Edge collapsing. On the left, the old cavity associated to the "displaced" vertex. On the right, the new cavity.

Remark 1 (Edge collapsing on boundaries). The edge collapsing pattern can be applied to any edge of the mesh (including boundaries). However, removing edges on a non-straight boundary would change the geometry of the domain and thus lead to a violation of the GCL. In this work, we only allow collapsing of boundary edges if the mesh modification does not change the geometry.

Now that the question "How to modify the mesh?" has been addressed, it remains to set when should a given pattern be applied?

Edge splitting and collapsing control the mesh size and edge swapping controls the quality of the mesh.

In order to decide if an edge has to be split or has to be collapsed, one must consider a relevant length function. This function relies on the definition of a discrete metric field [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF]. The adaptation procedure is described in Section 1.4.

In the next paragraph, we explain when an edge has to be swapped.

Mesh quality: edge swapping criterion

Classically, the Delaunay criterion is used to trigger edge swapping, but it is known that it is subject to numerical instabilities and that it does not extend straightfully to 3D. So, following [START_REF] Dobrzynski | Fast and accurate simulations of air-cooled structures[END_REF][START_REF] Alauzet | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Alauzet | Estimation d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage[END_REF] we use another quality criterion.

Let us denote by c ∈ {1, 2} the two possible configurations (with and without swapping). For each triangle T c j , j ∈ {1, 2} of each configuration we call r c j the radius of the inscribed circle of T c j and R c j the radius of its circumscribed circle. One can now define the quality of each triangle as the ratio q c j :=

r c j R c j .
The quality of the configuration c is then defined as the minimum quality of the two triangles

q c := min j∈{1,2} q c j . (8)
The best configuration c is the one that provides the greater quality q c c := arg max

c∈{1,2} q c . (9)
Swapping is required if the current mesh does not provide the best configuration5 (see Figure 4). Observe that this criterion often coincides with the Delaunay criterion.

Mesh adaptation: refining and coarsening of the mesh

The mesh adaptation procedure consists in controlling the size of the edges of the mesh.

As we already stated in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF], the adaptation is achieved following the concept of unit mesh. The idea is to build a mesh where all edges have a length of 1, for some length function. In that case, the mesh adaptation problem reduces to finding a relevant length function. Obviously this concept is not new, see for instance [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF].

Unit mesh and metric field

Let us recall that in dimension d, a metric-field is a field of symmetric positive-definite matrices of R d×d .

The Riemannian distance between two points a and b of R d , associated to the metric field M , is defined by

d M (a, b) = inf γ 1 0 γ (t), M γ(t) γ (t) 1/2 dt,
where γ is a C 1 -path joining a to b (γ(0) = a et γ(1) = b).

Therefore the length of an edge joining the points a and b, is given by

l M (a, b) = 1 0 (b -a), M a+t(b-a) (b -a) 1/2 dt. (10)
Obviously, as soon as M is a non-constant field, then l M (•, •) is not a distance: l M (•, •) does not satisfy triangle inequality. We refer to [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF] for a discussion on the consequences of the loss of the triangle inequality.

Actually, even if the metric field is uniform, it is generally not possible to build an exact unit mesh, so instead of imposing l M (a, b) = 1 for each edge (a, b), the constrain is relaxed as

η 2 ≤ l M (a, b) ≤ 2 η , ( 11
)
where η is a correction constant used to bypass the triangle inequality loss which could produced nonconvergence of the remeshing algorithm. In practice, we always take the value η = 1.2.

In this paper, we only treat the case of isotropic metrics. The reason for that is that our experiments dealing with anisotropic meshes were not successful: it requires a precise calculation of the anisotropic metric itself, which is not an easy task starting from finite-volume data. Moreover, using very distorted triangles may deteriorate the quality of the numerical approximation in the case of finite-volume semi-Lagrangian schemes.

In the isotropic case, the metric field is nothing else but a field of the form

M x = λ(x)I, with λ(x) > 0. Actually, h(x) := 1 √ λ(x)
is the desired edge length at the position x in the euclidean space (the physical space).

In practice, one has to evaluate a discrete metric field defined at each node of the mesh. Then to compute the length of the edge e = [a, b], one must interpolate the metric field along e to apply formula [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF]. There are many ways to do it. In this paper, we use

∀t ∈ [0, 1], M t := 1 1-t λa + t λ b I, where M a = λ a I, M b = λ b I,
With this choice, one can compute (10) analytically6 .

Metric field evaluation

An optimal way to adapt a mesh to improve the approximation would rely on a posteriori error estimates.

However, such estimators remain out of range for gas dynamics. Thus, one has to consider more empirical approaches. Following [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF] a good estimator is the geometric estimator. It consists in correlating locally the prescribed edge length to the curvature of the graph of a scalar quantity.

Let us denote by φ a finite-volume scalar quantity (piecewise constant by cell). For all vertex r, one estimates H φ r , an approximation of the Hessian ∇ 2 φ xr . In this work, we use a least squares method. One can then define the anisotropic discrete metric field associated to H φ r by

∀r ∈ R, M φ r aniso := P φ r Λ φ r P φ r -1 ,
where

H φ r = P φ r Λ φ r P φ r -1
and Λ φ r is the diagonal matrix of eigenvalues of H φ r : λ φ r 1 and λ φ r 2 . Finally, since we only consider here the isotropic case, we choose λ φ r := max(|λ φ r 1 |, |λ φ r 2 |), which is equivalent to imposing the most constraining mesh size in all directions at vertex r.

The isotropic discrete metric field associated to φ at node r is defined by

∀r ∈ R, M φ r iso := λ φ r I. ( 12
)
Building a unit mesh in the metric field M φ r iso provides an adapted mesh to the scalar quantity φ in the physical space.

Observe, that it is also possible to define a discrete metric field that prescribes a mesh size according to other considerations. For instance one may impose the mesh to be as refined as possible in parts of the computational domain.

Metric intersection

In practice, one generally wants to adapt the mesh according to various criteria (different physical quantities or geometrical constrains).

Let us assume that each of these criteria i is represented by a discrete isotropic metric field ∀r ∈ R,

M i r = λ i r I, with λ i r > 0. Then we define ∀r ∈ R, λ ∩ r := max i λ i r ,
so that the intersection metric field reads

∀r ∈ R, M ∩ r := λ ∩ r I. Since λ ∩ r = max i λ i r ⇐⇒ h ∩ r = min i h i r
, then the intersection metric imposes at each node r the smallest mesh size prescribed by all the criteria.

Metric filtering

Actually min r∈R h ∩ r can be arbitrary small and h ∩ r can be very oscillatory. So, M ∩ cannot be used in practice to define a mesh. One must filter the metric field.

The first filter is set to bound the size of the mesh edges. Thus, we define

∀r ∈ R, λ # = min max C ε λ ∩ r , 1 h 2 max , 1 h 2 min , (13)
where 0 < h min ≤ h max are the minimal and maximal edge lengths and where C ε is a constant coefficient used to control the interpolation error. One should refer to [START_REF] Alauzet | Estimation d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage[END_REF][START_REF] Alauzet | Anisotropic mesh adaptation for CFD computations[END_REF] for details, where in dimension d, C := d 2 2(d+1) 2 and ε is a prescribed interpolation error. This allows to define the sharp discrete metric field M # = λ # I.

However, this metric field can still be very oscillatory, i.e. it could lead to meshes with very fast variations of the edges size. This phenomenon is called h-shocks (see [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF] for an illustration). These kind of rapid mesh size variations can deteriorate the quality of the numerical approximation. So, it is common to replace this sharp metric field by a smooth one.

Definition 1 (Smooth isotropic metric). Let M be a discrete metric field defined at nodes

∀r ∈ R, M r = λ r I.
The metric field is said to be smooth if

∀l ∈ L, ∀r, s ∈ R l , h r ≤ h s + σ x r -x s , (14)
where σ ≥ 0 is a fixed constant, and where ∀r ∈ R, h r = 1 √ λr .

In practice we usually set the smoothing factor σ = 0.25.

The smoothing procedure consists in replacing the imposed edge lengths in order to reduce the size growth by satisfying [START_REF] Després | Machine learning design of volume of fluid schemes for compressible flows[END_REF], see Algorithm 1.

Algorithm 1: Metric smoothing.

1 i ← 0 2 M 0 ← M # 3 while ∀l ∈ L, ∃r, s ∈ R l , h i r > h i s + σ x r -x s do 4 if h i r > h i s + σ x r -x s then 5 h i+1 r ← h i s + σ x r -x s 6 else 7 h i+1 r ← h i r 8 i ← i + 1 9 M ← M i
After the last smoothing step, one finally gets the discrete metric field M that is used to adapt the mesh.

Adaptation constrains in multi-material case

As discussed before, we would like to preserve material interfaces as much as possible. This is achieved by imposing some constrains at material interfaces that reduce as much as possible the creation of mixing cells.

Let us first observe that the splitting pattern cannot create mixing cells from pure cells only. In this case, new triangles are defined as subdivision of old ones (see Figure 2). So, no additional constrain is imposed to edge splitting.

Obviously, if a swapping is applied between two pure cells containing two different materials, then the new ones will be both mixing cells. Thus we do not allow swapping of material interface edges.

The collapsing may or may not create mixing according to the geometry of interface itself. This is illustrated on Figure 5. In the case of coarsening of a cavity that contains only pure cells, there are two situations. Firstly, if the collapsed edge is an interface edge (see Figure 5a) then generally, only one of the new cells will be a mixing cell. Secondly, if the collapsed edge is not an interface edge (see Figure 5b) then generally most of the new cells will be mixing cells. Following this simple analysis, we adopt some rules to reduce mixing creation. These rules take place in cavities which contain inner interface edges and are applied only to swapping and collapsing patterns.

Swapping We forbid swapping of material interface edges.

Collapsing We allow the collapsing of a material interface edge, if it does not modify any other material interface. This collapsing of interface edge can generate mixing.

The Figure 6 Cleaning Due to the interfaces remeshing constrains defined above, it may occur that some triangles cannot be remeshed. This can lead to arbitrary small or very distorted triangles.

These non-remeshed triangles may have bad consequences. The quality of the numerical approximation in these triangles can be poor. Moreover, the local time-step can be very small.

To avoid these issues, if an edge becomes too small (for instance ten times smaller than the imposed minimal edge length h min ), then the interfaces constrains applied to the collapsing pattern are ignored.

We will refer this "strategy" as a cleaning pattern, while it is just the collapsing pattern with relaxed interfaces constrains.

Observe that the cleaning pattern is applied at the very end of the whole mesh adaptation procedure.

In practice, it is rarely triggered.

Remapping

In this paragraph, we shall present the remapping procedure we use. Actually, it follows a classical van Leer procedure [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme. V. a Second-Order Secquel to Godunov's Method[END_REF]. We will show its properties and present how to enforce its stability by means of Maximum Principles. Finally we will adapt a simple interface reconstruction [START_REF] Youngs | Time dependent Multi-Material Flow with Large Fluid Distortion[END_REF] method to our context.

In this paper, we focus on the resolution of the simple multi-constituent model ( 1)-( 5) to keep the presentation simple, but the method applies for more complex (realistic) models. Actually even for more complex problems, the remapping strategy remains the same: one remaps conservative variables using a conservative scheme and others quantities are deduced from the remapped ones. Note that density ρ is not remapped but simply deduced from remapped partial densities using (3).

In the case of ( 1)-( 5), the remapped quantities are (ρ i ) 1≤i≤N , ρE and ρu, and one deduces

ρ = N i=1 ρ i , u = ρu ρ , E = ρE ρ , . . .
The remapped quantities are piecewise constant by cell since we use a semi-Lagrangian finite-volume scheme [START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique Lagrangienne bidimensionelle[END_REF][START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF] to solve (1)-( 5).

Second-order remapping

As showed in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF], the mesh adaptation leads to very few mesh changes from one time step to the other.

Therefore the method is intrinsically cheap and the whole numerical scheme has low numerical diffusion compared to ALE methods (it remains almost semi-Lagrangian). Thus, one can wonder if a second-order remapping is really important. We will illustrate its benefits through numerical tests in Section 4, however one can already anticipate improvements: since the mesh is generally modified in regions of interest (shocks, shears, rarefaction waves,. . . ), the gain of precision will not be marginal.

It is enough to detail the remapping algorithms and their properties considering just the mesh changes inside a cavity (see Section 1.2). So, let us denote by O (respectively by N ) the set of old cells (respectively the set of new cells) used when applying one remeshing pattern to a given edge of the mesh.

Since these mesh modifications do not change the computational domain, one has

j∈O T O j = k∈N T N k ,
where (T O j ) j∈O and (T N k ) k∈N are respectively the triangles composing the old and the new cavities. In the same way, φ O j (respectively φ N k ) denotes the value of a quantity φ defined in the triangle T O j (respectively T N k ). Also, let us define the intersection of two triangles by

∀(j, k) ∈ O × N , I k j := T O j ∩ T N k .
Obviously, one has by construction

∀j ∈ O, k∈N I k j = T O j , and ∀k ∈ N , j∈O I k j = T N k ,
and

∀j 1 , j 2 ∈ O, ∀k 1 , k 2 ∈ N , I k1 j1 ∩ I k2 j2 = ∅ ⇐⇒ (j 1 , k 1 ) = (j 2 , k 2 ),
which just means that two different intersections do not overlap.

Finally, let us define the subsets of old cells of O that intersect a given new cell k ∈ N .

∀k ∈ N , O k = j ∈ O such that T N k ∩ T O j = ∅ .
The values of a remapped quantity φ are given by

∀k ∈ N , |T N k |φ N k = j∈O k I k j φO j (x), (15)
where ( φO j ) j∈O are conservative reconstructions of φ,

∀j ∈ O, |T O j |φ O j = T O j φO j (x). ( 16
)
Property 1 (Conservation). For any volumic quantity φ, the remapping scheme ( 15) is locally conservative

k∈N φ N k |T N k | = j∈O φ O j |T O j |. ( 17
)
Proof. The proof is straightforward. Starting from [START_REF] Dobrzynski | Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments[END_REF], Therefore

k∈N φ N k |T N k | = k∈N j∈O k I k j φO j (x).
k∈N φ N k |T N k | = j∈O k∈N I k j φO j (x). So, since k∈N I k j = T O j , k∈N φ N k |T N k | = j∈O T O j φO j (x),
and finally, since the reconstruction is conservative, using ( 16), one gets [START_REF] Dyadechko | Moment-of-fluid interface reconstruction[END_REF].

In the remaining of the paper we consider conservative affine reconstructions

7 ∀j ∈ O, φO j (x) :=    φ O j + g O j • (x -x j ), if x ∈ T O j , 0, else,
where g O j ∈ R 2 and where x j := 1

|T O j | T O j
x is the mass center of the cell j. 7 For the numerical tests presented in Section 4, we use a least-squares method to compute the conservative reconstructions.

Property 2. The remapping scheme ( 15) is exact for affine functions that are exactly reconstructed.

Proof. Assuming that φ is affine: φ(x) = a + g • x, with a ∈ R and g ∈ R 2 . Thus, if the reconstruction is exact,

∀j ∈ O, φ(x) = φO j (x) = φ O j + g • (x -x j ), with φ O j = a + g • x j .
Applying the remapping scheme [START_REF] Dobrzynski | Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments[END_REF], one has

∀k ∈ N , |T N k |φ N k = j∈O k I k j φO j (x) = j∈O k I k j a + g • x = T N k a + g • x.
So that φ N k = a + g • x k , which is the exact mean value of φ on T N k .

In practice, conservative reconstructions ( φO j ) j∈O are computed using a least squares method which is exact for affine data.

Stability

Another important property that one desires is stability. In this paragraph, we study the stability of the remapping scheme [START_REF] Dobrzynski | Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments[END_REF]. In a first time, we discuss its properties at first-order (setting ∀j ∈ O, g O j = 0 in ( 15)). Then, we define the Maximum Principles that should be satisfied for higher-order remapping.

Finally, we recall briefly the APITALI procedure [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF] adapted to our method to ensure stability of the second-order scheme.

First-order stability

Setting g O j = 0 in [START_REF] Dobrzynski | Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments[END_REF], one gets the first-order remapping scheme

∀k ∈ N , |T N k |φ N k = j∈O k |I k j |φ O j . (18)
We shall now show that (φ N k ) k∈N given by the scheme (18) satisfy local Maximum Principles.

Property 3 (Maximum Principle for a volumic quantity). Let O and N , be the sets of old and new cells of a cavity. Let φ be a scalar volumic quantity remapped using [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF]. One has ∀k ∈ N , min

j∈O k φ O j ≤ φ N k ≤ max j∈O k φ O j . (19)
Proof. According to [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF], one has

∀k ∈ N , φ N k = j∈O k |I k j | |T N k | φ O j .
Observing that

∀k ∈ N , ∀j ∈ O k , |I k j | |T N k | ≥ 0, and ∀k ∈ N , j∈O k |I k j | |T N k | = 1, then φ N k is a convex combination of (φ O j ) j∈O k
, which ends the proof.

According to Property 3, using [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF] to remap conservative quantities gives a first stability result: partial densities (ρ i ) 1≤i≤N , momentum components (ρu i ) 1≤i≤2 and total energy ρE satisfy Maximum Principles.

Actually one can establish more stability properties. One can show Maximum Principles for specific quantities at first-order. To do so, it remains to proof that density ρ satisfies also a Maximum Principle.

Property 4 (Maximum Principle for density). If partial densities (ρ i ) 1≤i≤N are remapped using [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF], then new densities defined by

∀k ∈ N , ρ N k = N i=1 ρ i N k , satisfy ∀k ∈ N , min j∈O k ρ O j ≤ ρ N k ≤ max j∈O k ρ O j .
Proof. The proof is a direct consequence of the linearity of the scheme [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF].

∀k ∈ N , |T N k |ρ N k = N i=1 |T N k |ρ i N k = N i=1 j∈O k |I k j |ρ i O j = j∈O k |I k j | N i=1 ρ i O j , ∀k ∈ N , |T N k |ρ N k = j∈O k |I k j |ρ O j .
Actually at first-order, computing the new density as the sum of new partial densities is equivalent to remapping the density using [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF]. So, it satisfies a Maximum Principle according to Property 3.

Property 5 (Maximum Principle for a specific quantity). Let ψ be a specific quantity. If φ = ρψ is remapped using [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF], then the new specific quantities obtained by

∀k ∈ N , ψ N k = φ N k ρ N k , satisfy ∀k ∈ N , min j∈O k ψ O j ≤ ψ N k ≤ max j∈O k ψ O j . (20)
Proof. One has ∀j ∈ O, φ O j = ρ O j ψ O j . So, [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF] gives

∀k ∈ N , |T N k |ρ N k ψ N k = j∈O k |I k j |ρ O j ψ O j , or ∀k ∈ N , ψ N k = j∈O k |I k j |ρ O j |T N k |ρ N k ψ O j .
One has

∀k ∈ N , ∀j ∈ O k , |I k j |ρ O j |T N k |ρ N k ≥ 0, and ∀k ∈ N , j∈O k |I k j |ρ O j |T N k |ρ N k = 1.
The left inequality is a direct consequence of Maximum Principle for the density (since ρ O j > 0). The right equation is due to the fact that at first-order, density satisfies [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF].

So, ψ N

k is a convex combination of (ψ O j ) j∈O so it satisfies [START_REF] Gillespie | The Gibbs-Dalton Law of Partial Pressures[END_REF]. This establishes that (u i ) 1≤i≤2 and E satisfy Maximum Principles.

Actually in order to ensure the stability of the whole method, it remains to prove that the internal energy also satisfies a Maximum Principle.

Property 6 (Maximum Principle for specific internal energy). If (ρ i ) 1≤i≤N , (ρu i ) 1≤i≤2 and ρE are remapped using [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF], then the new internal energies defined by

∀k ∈ N , e N k = E N k - 1 2 u N k 2 , satisfy ∀k ∈ N , min j∈O k e O j ≤ e N k ≤ max j∈O k e O j + 1 4 max i,j∈O k u O i -u O j 2 . ( 21
)
Proof. The proof is a bit more technical. It is given in Appendix A.

Maximum Principles for second-order

In the previous paragraph, we established the stability of the first-order remapping scheme through a set of Maximum Principles. However, as usual, the obtained bounds in ( 19), ( 20) and ( 21) are too strict for remapping of orders greater than 1. For instance, in the case of the refinement of a single cell j, to satisfy Maximum Principle ( 19) is equivalent to set ∀k ∈ N , φ N k = φ O j , so that g φ j = 0. In that case, the first-order scheme is the only one that satisfies [START_REF] Fritts | The Free-Lagrange Method[END_REF].

Therefore, to achieve higher-order of accuracy, one must consider a weaker set of Maximum Principles.

As usual, this consists simply in extending the set of admissible values by taking into account cell values in a vicinity of the modified cells. Thus, one defines the set of cells (before adaptation) that share a vertex with cell j

V j := {i ∈ J /R j ∩ R i = ∅},
One has obviously j ∈ V j .

We now present the Maximum Principles that we consider for higher-order remapping schemes. We distinguish three cases, the scalar case (for volumic remapped quantities, for deduced specific quantities and for the density as the sum of partial densities), the vector case (for the velocity) and the case of the internal energy.

Definition 2 (Maximum Principle for a scalar quantity). After projection, a scalar quantity φ is said to satisfy Maximum Principle if one has

∀k ∈ N , φ N k ∈ min j∈O k min i∈Vj φ O i , max j∈O k max i∈Vj φ O i .
Definition 3 (Local Convex Hull Principle for a vector quantity [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF]). After projection, a vector quantity φ is said to satisfy Maximum Principle if, one has

∀k ∈ N , φ N k ∈ C φ k ,
where C φ k is the convex hull defined by

(φ i ) i∈ j∈O k Vj .
Observe that this definition (following [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF][START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF]) is a generalization of the scalar Maximum Principle to vectors: the convex hull for scalar quantities is nothing else but the interval formed by the minimal and maximal values of φ on the set V j .

In the case of the velocity, one could also use a scalar component-wise approach to satisfy the Maximum Principle, but it would break the Galilean invariance of the method.

Finally let us define the Maximum Principle for internal energy in the case of approximating the system ( 1)-( 5)

Definition 4 (Maximum Principle for specific internal energy). After remapping, specific internal energy is said to satisfy Maximum Principle, if one has

∀k ∈ N , e N k ∈ min j∈O k min l∈Vj e O l , max j∈O k max l∈Vj e O l + 1 4 max i,j∈O k max l∈Vj u O i -u O l 2 .
It is worth noticing the following obvious result.

Property 7. The first-order remapping scheme [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF] satisfies the Maximum Principles in the sense of Definitions 2, 3 and 4, for partial densities and density, velocity, specific total energy and internal energy.

Proof. The proof is obvious since ∀k ∈ N , ∀j ∈ O k , j ∈ V j . For the vector case, it is enough to remark that according to [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF], u N k is a convex combination of (u O j ) j∈O k and thus belongs to C u k .

A priori limitation

For second-order schemes, we use a priori slope limiters. See [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme. V. a Second-Order Secquel to Godunov's Method[END_REF] for scalar quantities and see [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] for vectors: we use Vector Image Polygon (VIP) in order to preserve Galilean invariance.

It allows to define the limited reconstructions ∀j ∈ O,

ρ i L j (x) := ρ ij + α ρi j g ρi j • (x -x j ), ∀1 ≤ i ≤ N, ρu i L j (x) := ρu ij + α ρu j g ρui j • (x -x j ), ∀1 ≤ i ≤ 2, ρE L j (x) := ρE j + α ρE j g ρE j • (x -x j ),
where 0 ≤ α ρi j , α ρu j , α ρE j ≤ 1 are the limitation factors. Note that α ρu j is common to all momentum's components.

This standard a priori limitation procedure may not be enough to ensure stability in the sense of the Maximum Principles defined in the previous paragraph. So, one can additionally use an a posteriori limitation procedure to ensure stability as it is described in the following.

A posteriori limitation

One can actually satisfy the Maximum Principles defined in Section 2.2.2, using a posteriori limitation. In this work, we use the APITALI 8 [START_REF] Hoch | An Arbitrary Lagrangian-Eulerian strategy to solve compressible fluid flows[END_REF][START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF] procedure to achieve it. The idea of the method can be summarized as follows.

Let S ≥ 2 be the maximum number of trials. Let φ be a remapped scalar quantity, then one defines a finite sequence of reconstructions ∀1 ≤ s ≤ S,

∀j ∈ O, φs j (x) := φ O j + β s α φ j g φ j • (x -x j ), (22)
where (β s ) 1≤s≤S , is a finite non-negative and strictly decreasing real sequence such that β 1 = 1 and β S = 0.

One recalls that α φ j is the a priori limitation coefficient. The method consists in an iterative limitation procedure which locally degrades the reconstruction to first-order until Maximum Principles are satisfied.

More precisely, it consists in setting s j = 1 in each cell j ∈ O for each remapped quantity (ρ i ) 1≤i≤N , ρu and ρE using [START_REF] Hansen | Alexa: MPI+X Shock Hydrodynamics on Dynamically Adaptative Tetrahedral Meshes[END_REF]. In other words, in each old cell j, the first trial set of reconstructions is ( ρi

sj =1 j ) 1≤i≤N , ρu sj =1 j and ρE sj =1 j . Then if one of the computed quantities (ρ i N k ) 1≤i≤N , ρ N k , u N k , E N k or e N
k in the new cell k does not satisfy its associate Maximum Principle, one uses the next reconstruction (s j + 1) of all the quantities it depends on, in all the old cells j ∈ O k . This procedure is repeated until all Maximum Principles are satisfied.

APITALI converges in at most S steps since the Maximum Principles are satisfied at least for the first-order scheme, see Property 7.

Remark 2. Setting S = 2, that is (β s ) 1≤s≤2 = (1, 0), one retrieves the MOOD method [START_REF] Clain | A high-order finite volume method for systems of conservation laws -multi-dimensional optimal order detection (mood)[END_REF] for scalar quantities.

Interfaces Reconstruction

As discussed in Section 1.5, starting from pure cells, numerical mixing can only be produced when the collapsing pattern 9 is used (coarsening). But once non-pure cells are created, mixing can be propagated by any of the remeshing operations. To reduce the numerical diffusion of the interfaces, it is common to use interface reconstruction methods [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF]. Interface reconstruction itself is a complex task and remains active field of research [START_REF] Benson | Volume of fluid interface reconstruction methods for multi-material problems[END_REF][START_REF] Dyadechko | Moment-of-fluid interface reconstruction[END_REF][START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF][START_REF] Ahn | Multi-material interface reconstruction on generalized polyhedral meshes[END_REF][START_REF] Després | Machine learning design of volume of fluid schemes for compressible flows[END_REF]. The aim of this paper is not to compare methods nor to propose innovative ones to compute material interface positions. In the tests presented in Section 4, we use a VoF method [START_REF] Youngs | Time dependent Multi-Material Flow with Large Fluid Distortion[END_REF], but one can use his favorite method instead. 8 A Posteriori ITerAtive LImiter 9 Let us recall that we choose not to apply swapping for material interface edges.

20

Which ever the interface reconstruction method is, in the end, any mixing cell

T O j is subdivided into N pure sub-cells (T O j,i ) 1≤i≤N , such that N i=1 T O j,i = T O j and i 1 = i 2 ⇐⇒ T O j,i1 ∩ T O j,i2 = ∅.
This allows to define the volume fraction occupied by each material

∀1 ≤ i ≤ N, f ij := |T O j,i | |T O j |
.

In each old mixing cell j, one has

∀1 ≤ i ≤ N, f ij ≥ 0 and N i=1 f ij = 1
From mass conservation of each material, one defines the density of each material i in the old mixing cell j by

∀1 ≤ i ≤ N, ρi O j :=      |T O j | |T O j,i | ρ i O j if |T O j,i | > 0, 0 else, =      ρi O j fi j if f ij > 0, 0 else.

Mixing cells Interface reconstruction Material mesh

∩

New cavity With this reconstruction, one gets the new partial densities by means of exact intersection of the new triangles and the old sub-cells -see Figure 7. It reads

∀k ∈ N , ∀1 ≤ i ≤ N, ρ i N k = j∈O k |T N k ∩ T O j,i | |T N k | ρi O j . ( 23
)
It is quite obvious that according to [START_REF] Hoch | An Arbitrary Lagrangian-Eulerian strategy to solve compressible fluid flows[END_REF], when using interface reconstruction, partial densities do not satisfy Maximum Principle. However to obtain stability, one requires a Maximum Principle on the density which would give stability to other quantities: E, u and e, as discussed in Section 2.2.1.

Actually one can show the following result.

Property 8. ∀j ∈ O, ∀1 ≤ i ≤ N , denoting respectively by

f i O j := |T O j,i | |T O j | and c i O j := ρ i O j ρ O j ,
the old volume fractions and old mass fractions. If the interface reconstruction method is such that volume fractions are equal to mass fractions

∀j ∈ O, ∀1 ≤ i ≤ N, f i O j = c i O j , (24)
and if new densities are defined as the sum of the new partial densities given by ( 23), then the new densities satisfy the following Maximum Principle ∀k ∈ N , min

j∈O k ρ O j ≤ ρ N k ≤ max j∈O k ρ O j .
Proof. The proof is given in Appendix B.

In the general case, ( 24) is only a sufficient condition, but it is easy to check that in many configurations it is actually a necessary condition -see for instance the example given in Appendix C. Therefore in order to ensure the Maximum Principle, one should use an interface reconstruction method that satisfies [START_REF] Hoch | Semi-conformal polygonal mesh adaptation seen as grid velocity formulation for ale simulations[END_REF]. Moreover, since partial densities are discontinuous at material interfaces, we enforce first-order for all remapped quantities in cavities where interface reconstruction takes place. The purpose is to avoid unnecessary iterations of the APITALI procedure in these cavities, since generally, only first-order satisfies Maximum Principle at material interfaces.

Remark 3. In practice, when dealing with interface reconstruction, the method itself is used in cavities that contain mixing cells and second-order elsewhere.

One can wonder if a Maximum Principle for the remapped density is relevant in the multi-material context. However, if ρ satisfies a Maximum Principle, then one obtains a Maximum Principle for the internal energy (see Property 6) which is important for stability.

Implementation considerations

Let us emphasize that the method we describe in this paper does not decouple the mesh adaptation phase from the remapping phase. Remapping is performed after each step of the mesh modification. This has two main consequences. First, one performs remapping locally (from one cavity to another) which simplifies calculations and parallelism. Second, since the mesh is adapted up to convergence to the numerical solution at current time, we expect just a few modifications to adapt the mesh to numerical solution at next time.

Finally, the adaptation procedure is summarized in Algorithm 2.

Algorithm 2: Adaptation procedure.

1 Function swap()

2
while some edges must be swapped do

Numerical tests

We shall now present a few numerical tests to illustrate the good behavior of the method. The tests are organized as follow. Starting from academic reference tests, we will go through result comparisons of more complex flows and finally provide an example of a much more challenging multi-material flow which has no reference solution.

In all the following tests, an adaptation loop is performed at time t = 0 in order to prepare the mesh for the calculation. For each test, the adaptation criteria are based on ρ, (c i ) 1≤i≤N and E10 . In the case of multi-material tests, the use of (c i ) 1≤i≤N allows the mesh to be automatically refined at material interfaces.

The mesh adaptation options are σ = 0.25 in ( 14) and η = 1.2 in [START_REF] Clain | A high-order finite volume method for systems of conservation laws -multi-dimensional optimal order detection (mood)[END_REF]. The cleaning threshold is set to 10 -1 h min .

Finally, when APITALI is used to ensure Maximum Principles, the decreasing sequence is set to

     ∀1 ≤ s < 10, β s = 0.75 s-1 ,
and

β 10 = 0.
Order 1 is forced when reaching the 10 th iteration.

Remark 4. Interpolation error ε remains difficult to predetermine. However in practice, it is quite easy to define by testing few values. In this work, we prefer to underestimate its value, which has the consequence of slightly over refining the mesh.

The following results were obtained using a second-order version of the Glace [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF] scheme to solve (1)-( 5), similar results are obtained using Eucclhyd [START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF] scheme.

Two-dimensional Riemann problems

In order to validate and analyze the behavior of the method, we first use three 2D Riemann problems that were defined by P. Lax and X.-D. Liu in [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF]: following [START_REF] Liska | Comparison of several difference schemes on 1D and 2D test problems for the Euler equations[END_REF] we consider the cases 3, 12 and 15.

In [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF] all these tests are performed in the domain ]0, 1[ 2 and discontinuity lines are set along the axes x = 1 2 and y = 1 2 . Since our method does not allow mass fluxes naturally, the initial domain is modified in order to allow comparisons with the reference solutions at final time. To do so, we always impose the appropriate velocity boundary conditions. The fluid follows a perfect law with γ = 1.4 in all these test cases.

We set the minimal mesh size to h min = 2.5 × 10 -3 (following [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF]) and the interpolation error in (13) to ε = 2.5 × 10 -3 .

Configuration 3

Observe that since the solutions of these tests present shears or vortexes, a completely Lagrangian method cannot succeed in their calculation. Initial data is given on Figure 8. The final time is set to t = 0.3. In this case, the number of cells grows from 19524 to 37324, compared to the 159201 cells of the reference [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF].

The solution we obtain is presented on Figure 9. The four 1D Riemann problems yield four single shock waves. One observes the trace of the wall-heating which is exactly preserved since once created. These numerical artifacts would disappear with few numerical dissipation as it is observed in ALE simulations. This shows that in these regions the scheme is almost Lagrangian and that the mesh is just transported.

The solution is very close to the reference one.

As it was already observed in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF], the whole method has very little numerical dissipation. In one hand, the second-order Lagrangian scheme itself produces small numerical dissipation. In the other hand, the convergence of the adaptation loop at each iteration provides very few mesh changes from one iteration to the other. Therefore the remapping procedure is somehow a rare event. However one can observe on Figure 10, the benefits of the second-order remapping: the instability is more dissipated at first-order (the mushroom is more developed and shocks are sharper at second-order). This could be explained by the fact that the mesh is modified in regions of interest and so, that a more accurate remapping is relevant. However one should note that the use of APITALI procedure, while ensuring stability, has very low impact on the quality of the solution.

Remark 5 (Spurious mesh refinement). One observes on Figures 9 and 10 spurious mesh refinements.

These are due to some kind of wall-heating: the pressure is at equilibrium around these points keeping density and energy variations unchanged. These density or energy variations produce locally a non null hessian which imposes to keep the mesh refined. Since the pressure is constant, the solution and the metric field are just transported in these regions: the mesh has no reason to be coarsened. Again a few numerical dissipation could cure these artifacts. Initial data is given on Figure 11. The gas follows a perfect law with γ = 1.4. The final time is set to t = 0.2. The number of cells grows from 15272 to 36573, compared to the 159201 cells of the reference [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF].

On Figure 12, we compare our solution to the reference one. The bottom and left 1D Riemann problems result in a sliding contact discontinuity. Observe that without any special treatment, the sliding is quite well captured due to the mesh adaptation. Nevertheless, as expected, one observes Kelvin-Helmoltz instabilities due to the sliding.

The agreement to the reference solution is quite good and the mesh is correctly adapted. Again, one observes the remaining of the wall-heating which in some way reflects the very low dissipation of the method. Initial data is given on Figure 13. The gas follows a perfect law with γ = 1.4. The final time is set to

t = 0.2.
The solution we obtain is close to the reference one as depicted on Figure 14. In this case, the number of cells grows from 14950 to 39487, compared to the 159201 cells of the reference [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF].

As previously the bottom and the left Riemann problems yield single contact discontinuity waves. But this time these slide lines have a normal velocity which is somehow more complex to deal with. The top Riemann problem generates a single rarefaction wave and the right one, a shock.

One can see that, without special treatment, the sliding is again well approximated. Also, one still observes the imprint of the wall-heating in the mesh adaptation.

Noh problem

The aim of this test case is to illustrate the benefits of ensuring Maximum Principles during the remapping phase. In this regard, it is enough to consider a Noh problem [START_REF] Noh | Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux[END_REF] in Cartesian coordinates.

This test is performed in 2D, the computational domain is ]0, 1[ × ]0, 0.1[ and initial state is given by ρ = 1, p = 10 -9 and u = (-1, 0). The pressure follows a perfect gas law with γ = 5 3 . The final computation time is 0.6, we set h min = 10 -3 and the interpolation error to ε = 2 × 10 -2 . Without Maximum Principles, the second-order remapping fails.

The results of the second-order remapping satisfying Maximum Principles are presented on Figures 15 and16. The numerical solution is in accordance with the theoretical one. One can observe the very good agreement of the numerical solutions. This is quite remarkable since the flow itself is quite complex and since there is no stabilization term in the approximated model -we consider Euler equations (no viscosity nor conductivity) without surface tension. Moreover the involved numerical methods are very different.

Mixer test

We designed this test to illustrate the robustness and the flexibility of the method. The test consists in two rigid bars rotating in a fluid domain. The fluid itself is made of three constituents of various densities and initially at rest. Dealing with Euler equations it is natural to impose perfect sliding at the bars boundaries.

Actually, with classic numerical approaches this test is very challenging. Indeed, due to the large displacements, ALE method cannot handle the mesh deformations, which makes reconnection-based methods mandatory in the case of conforming approaches. One could use fictitious-domain-like methods, but in that case the treatment of perfect sliding imposed at the boundaries of the bars is very difficult to take into account.

Obviously there is no reference solution for this test. We illustrate here the effect of the interface reconstruction compared to a second-order remapping of partial densities.

The initial configuration is given on Figure 20. The boundary conditions on the bounding box are set to u • n = 0. On the rotating bars, one imposes normal velocities such that they rotate at the angular velocities 1 and 7 4 . The final time is set to t = 2π so that the left bar performs a complete rotation, whereas the right one almost performs 2. We impose h min = 10 -2 and, in [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF], the interpolation error is set to 5 × 10 -3 . Even if these calculations could be performed without it, we use APITALI to enforce Maximum Principles of the remapping step.

We compare the numerical solutions we obtain using interface reconstruction or not (see Figure 21).

One observes that at the end of the calculations the fluids are extremely mixed. Nevertheless, the good Right: Eulerian AMR reference from [START_REF] Banks | A high-resolution method for compressible multimarterial flow on overlapping grids[END_REF].

accordance of the big structures is quite noticeable for such a complex flow. One can see the benefits (Table 2) of the interface reconstruction method since the number of mixing cells remains quite small all along the calculation (44.7% at final time), while 82.7% of the cells are already mixing cells at time π 2 without interface reconstruction. One can also observe that the number of cells at final time is greater when using interface reconstruction. Actually, we interpret this result as a consequence of a better description of Right: Eulerian AMR reference from [START_REF] Banks | A high-resolution method for compressible multimarterial flow on overlapping grids[END_REF].

the interfaces all along the simulation. In fact, the numbers of cells are very close until time t = 5.2.

Conclusions

In this paper we presented improvements of our previous work [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF]: second-order remapping and multimaterial management.

Concerning the remapping step, we proved a set of Maximum Principles satisfied at first-order which allowed to derive an APITALI-like procedure to ensure the stability of the second-order scheme.

The mesh adaptation now allows creation of mixing during the local coarsening of the mesh while preserving Lagrangian interfaces as much as possible. Remapping in cavities that contain mixing cells uses a simple interface reconstruction method to reduce the numerical diffusion.

We performed various numerical tests that demonstrate the quality of the method as well as its robustness.

As already observed in [START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF], the method achieves to be almost Lagrangian: less that 1% of the cells are remeshed at each time step. Moreover, since the mesh is modified through local patterns, the method is One represents the mass fraction of the middle material and the mesh evolution.

Left: without interface reconstruction. Right: using Youngs' method.

quite cheap and permits a parallel implementation 11 .

11 The numerical results do not depend on the number of processors: they are exactly the same for any number of processors.

Future works could concern the 3D extension. First results have been published in [START_REF] Hansen | Alexa: MPI+X Shock Hydrodynamics on Dynamically Adaptative Tetrahedral Meshes[END_REF]. Nonetheless in 3D, the difficulty lies in the quality pattern: one of the key points of the method is the convergence of the three mesh modification steps (quality improvement, coarsening and refinement of the mesh) to obtain a quasi-Lagrangian method. In 3D, it is more difficult and more expensive to achieve convergence of the edges swapping procedure.

One could also try to improve the metric calculation. For instance, in the case of very different materials, it could be relevant to compute discrete metric fields separately in each material to be able to detect waves that can be neglected when considering the whole data.

Finally, the use of better interface reconstruction algorithms could improve the accuracy of the method at interfaces.

∀k ∈ N , d N k = j∈O k α k j u O j 2 + 1 2 j,l∈O k α k j α k l u O j -u O l 2 -u O j 2 -u O l 2 ,
which rewrites

∀k ∈ N , d N k = j∈O k α k j u O j 2 + 1 2 j,l∈O k α k j α k l u O j -u O l 2 - j∈O k l∈O k α k l α k j u O j 2 .
Observing that since l∈O k α k l = 1, one has This ends the proof.

Appendix B. Proof of Property 8

We recall that Property 8 establishes a Maximum Principle for density when interface reconstruction is such that volume fractions are set equal to the mass fractions.

Let us denote by k ∈ N one of the new cells and by i ∈ {1, . . . , N } one of the material. According to [START_REF] Hoch | An Arbitrary Lagrangian-Eulerian strategy to solve compressible fluid flows[END_REF], one has

ρ i N k = j∈O k |T N k ∩ T O j,i | |T N k | ρi O j , where ρi O j =      ρi O j fi j if f ij > 0, 0 else.
Let us denote by O k,i := j ∈ O k , s.t. f i O j > 0 , the set of old cells that intersect k and that contain material i. Then one has

ρ i N k = l∈O k,i |T N k ∩ T O j,i | |T N k | ρi O j = l∈O k,i |T N k ∩ T O j,i | |T N k | 1 f i O j ρ i O j .

Since by hypothesis c

i O j = f i O j , one has ρ i N k = l∈O k,i |T N k ∩ T O j,i | |T N k | 1 c i O j ρ i O j .
By definition of the mass fractions, c i

ρ i N k = j∈O k |T N k ∩ T O j,i | |T N k | ρ O j .
So, the new density in cell k reads

ρ N k = N i=1 ρ i N k = j∈O k N i=1 |T N k ∩ T O j,i | |T N k | ρ O j = j∈O k |T N k ∩ T O j | |T N k | ρ O j .
Therefore, the new density ρ N k is a convex combination of (ρ O j ) j∈O k .

This ends the proof.

Appendix C. Refinement of a cell containing two materials

If the new cell k only contains material i, then

ρ N k = ρ i N k = |T N k ∩ T O j,i | |T N k | c i O j f i O j ρ O j ,
according to [START_REF] Hoch | An Arbitrary Lagrangian-Eulerian strategy to solve compressible fluid flows[END_REF]. Since in that case, T N k ⊂ T O j,i , one has

ρ N k = c i O j f i O j ρ O j .
However at first-order if the scheme satisfies Maximum Principle, one has

ρ O j ≤ ρ N k ≤ ρ O j ,
which implies

f i O j = c i O j .
So in this case, ( 24) is a necessary condition to satisfy Maximum Principle at first-order.

Figure 1 :

 1 Figure 1: Edge swapping. On the left, the initial cavity. On the right, the new cavity.

Figure 2 :

 2 Figure 2: Edge splitting. On the left, the initial cavity (case of an inner edge). On the right, the new cavity. The new vertex is positioned at the middle of the initial edge according to some prescribed metric, see Section 1.4.1.

Figure 4 :

 4 Figure 4: Quality criterion for swapping. The right configuration as the best quality. In this case, the right configuration also satisfies the Delaunay criterion.

  Collapsing an interface edge. Three of the four new cells are kept pure. Actually, if nodes A, R and D are aligned, no mixing cell is created. Collapsing a non interface edge from an interface node.

Figure 5 :

 5 Figure 5: Collapsing pattern: eventual mixing creation. Pure cells are colored in blue and red, mixing cells are colored in green.

Figure 6 :

 6 Figure 6: Rules to reduce mixing creation. Cells are pure cells colored by material. The blue edges cannot be swapped. The blue edges can be collapsed as long as the node C is untouched. Red edges (AD and DE) cannot be collapsed but can be swapped. Black edges can be collapsed as soon as the displaced node is not an interface node: for instance to collapse edges BA or BC, one has to move node B to A or C.

  and since ∀j ∈ O \ O k , I k j = ∅, one has

Figure 7 :

 7 Figure 7: From left to right. A cavity containing mixing cells. A reconstruction of material interface. A local conforming meshthat respects the material interfaces (each sub-cell contains only one material). This later mesh is suitable to compute new partial densities using exact intersections.

3 Get a set of edges than can be swapped simultaneously 4 Swap these edges 5 Remap the reliable quantities locally 6 Function collapse() 7 while some edges must be removed do 8 Get a set of edges than can be removed simultaneously 9 Remove these edges 10 Remap the reliable quantities locally 11 Update metric according to new values 12 Call swap() 13 Function split() 14 while some edges must be split do 15 Get a set of edges than can be split simultaneously 16 Split edges 17 Remap the reliable quantities locally 18 Update metric according to new values 19 Call swap() 20 Function clean() 21 Call

 3456789101112131415161718192021 collapse() for edges smaller than 1 10 in the metric and ignoring interfaces constrains 22 Function adapt()

Figure 8 :

 8 Figure 8: Configuration 3: initial domain for the Lagrangian-AMR calculation.

Figure 9 :

 9 Figure 9: Configuration 3: comparison of the numerical solution with second-order remapping with Maximum Principles to the reference at time t = 0.3. Left: the final mesh colored with the density. Middle: density isolines. Right: density isolines from [28].

Figure 10 :

 10 Figure 10: Comparison of the instability at time t = 0.3, according to the remapping strategy. Left: first-order. Middle: second-order with Maximum Principles. Right: second-order.

Figure 11 :

 11 Figure 11: Configuration 12: initial domain for the Lagrangian-AMR calculation.

Figure 12 :

 12 Figure 12: Configuration 12: comparison of the numerical solution with second-order remapping and Maximum Principles to the reference. Left: the final mesh colored with the density. Middle: density isolines. Right: density isolines from [28].

3 Figure 13 :

 313 Figure 13: Configuration 15: initial domain for the Lagrangian-AMR calculation. A security margin is taken to avoid tangling of cells at the boundary.

Figure 14 :

 14 Figure 14: Configuration 15: comparison of the numerical solution with second-order remapping and Maximum Principles to the reference. Left: the final mesh colored with the density. Middle: density isolines. Right: density isolines from [28].

Figure 15 :Figure 17 :

 1517 Figure 15: The density at time t = 0.6 in all cells compared to the exact solution.

  (a) Comparison at time t = 2. (b) Comparison at time t = 4. (c) Comparison at time t = 6.

Figure 18 :

 18 Figure 18: Shock-bubbles interaction. Comparison of Schlieren in the case of helium bubbles. Left: semi-Lagrangian AMR.

  (a) Comparison at time t = 2. (b) Comparison at time t = 4. (c) Comparison at time t = 6.

Figure 19 :

 19 Figure 19: Shock-bubbles interaction. Comparison of Schlieren in the case of refrigerant bubbles. Left: semi-Lagrangian AMR.

Figure 20 :

 20 Figure 20: Mixer test: initial configuration. The fluid is at rest: u = 0 and p = 1. Each fluid follows a perfect gas law (γ = 1.4).

  (a) Time t = 0. (b) Time t = π 2 . (c) Time t = π. (d) Time t = 3π 2 .(e) Time t = 2π.

Figure 21 :

 21 Figure 21: Comparison of the numerical solutions with and without interface reconstruction. Both simulations are computed at second-order and satisfy Maximum Principles.

Sol

  Since ∀v, w ∈ R n , one has vw 2 = v 2 -2(v, w) + w 2 , it comes that

j∈O k l∈O k α k l α k j u O j 2

 2 left part of[START_REF] Grospellier | The arcane development framework[END_REF]. Finally, one has ∀k ∈ N ,

  since j ∈ O k \ O k,i ⇐⇒ |T N k ∩ T O j,i | = 0,

Table 2 :

 2 Evolution of the number of mixed cells for both approaches according to time.

From the discrete point of view, it means that c i will remain constant in each cell as long as the calculation is Lagrangian.

Semi-Lagrangian methods are also referred as Lagrange plus update or updated Lagrangian in the literature.

The calculation is no more purely lagrangian: mass fluxes occur while remapping the numerical solution onto the new mesh.

Similarly to[START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF], the numerical results presented in this paper are totally independent of the number of processors.

If both configurations quality are equal, the edge is not swapped.

There are other possibilities (see[START_REF] Alauzet | Estimation d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage[END_REF][START_REF] Alauzet | Anisotropic mesh adaptation for CFD computations[END_REF] or[START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF] for instance) for which analytical calculation of the edge length is also possible.

One can additionally add u as a criterion, but it breaks the Galilean invariance of the method. For the tests presented in this paper, the effect of a criterion based on u is marginal.

Acknowledgements

The Authors are grateful to Philippe Hoch and Emmanuel Labourasse for the discussions we had on this topic, especially concerning Maximum Principle results at first-order and the adaptation of APITALI procedure to this context.

Shock interaction with cylindrical bubbles

This is a multi-material test proposed in [START_REF] Banks | A high-resolution method for compressible multimarterial flow on overlapping grids[END_REF]. It consists in computing the interaction of a shock wave propagating in the air with six cylindrical bubbles of a given material. Actually two different materials are considered for the bubbles. First, a set of six helium bubbles and then, using the same geometrical configuration, a set of six bubbles of a refrigerant material (a much more stiff and heavy material). The computational domain is Ω = ]0, 10[ 2 and the initial geometry of the bubbles is given in surrounding the bubbles is initially at rest: ρ = 1, u = 0 and p = 1, and follows a perfect gas law with γ = 1.4. The helium bubbles are initially at equilibrium ρ = 0.138, u = 0 and p = 1. The helium follows a perfect gas law γ = 1.67. In the case of the refrigerant material, one sets ρ = 3.15, u = 0 and p = 1. It also follows a perfect gas law with γ = 1.25. In both configurations a shock is initialized at position x = 0.4.

The shock Mach number is 1.22. See Figure 17.

Symmetry boundary conditions are imposed at the top and the bottom of the computational domain.

Since in our method we do not impose mass fluxes, it is modeled by the use of a piston whose speed is given by the post-shock state on the left. Also to avoid reflections on the right boundary, we extend the computational domain and impose u = 0 on the right.

For the calculations, the minimal mesh size is set to h min = 2.5 × 10 -3 following [START_REF] Banks | A high-resolution method for compressible multimarterial flow on overlapping grids[END_REF]. The interpolation error is set to ε = 5 × 10 -4 . Although it is not necessary to run this test, we activate APITALI in order to satisfy Maximum Principles after remapping. Finally, interface reconstruction is used to reduce numerical diffusion.

The number of cells grows from 317575 to 1471717 in the case of the helium and from 317575 to 1393532 in the case of the refrigerant material.

We compare our results to the ones obtained in [START_REF] Banks | A high-resolution method for compressible multimarterial flow on overlapping grids[END_REF] at times t = 2, t = 4 and t = 6. Figure 18 compares the numerical solutions in the case of helium bubbles, and Figure 19 in the case of refrigerant bubbles.

Appendix A. Proof of Property 6

We recall that Property 6 establishes Maximum Principle for specific internal energy.

To simplify notations, let us define

At first-order, assuming that ∀j ∈ O, ρ O k ≥ 0 and since ∀k ∈ N , ρ N k satisfies [START_REF] Dyadechko | Multi-material interface reconstruction from the moment data[END_REF], one has

Using this notation and ( 18), at first-order new total energies and new velocities are given by

So the new internal energies write Let j be a mixing cell that is refined. Let us assume a configuration such that one of the two new obtained cells is pure.


Rodrick Wallace
email: rodrick.wallace@gmail.com

Preventing Black Plague 2.0: Stochastic sterilization of agribusiness-driven pandemics

Keywords: agricultural reform, control theory, Data Rate Theorem, infectious disease, noise induced stabilization, phase transition

Current 'economies of scale' in global agribusiness -across space, time, genetics and organizational structure -constitute an insurmountable barrier to public health efforts at control of pandemic disease. A relatively elementary model of noise-induced stabilization suggests that sufficient and properly designed and monitored disruption of global agricartels -breakup -would act as a de-facto sterilizing temperature to prevent pathogen emergence.

Introduction

Much recent work has highlighted the role of agribusiness 'economies of scale' in grossly accelerating the rate of pathogenic challenge to human populations, a dynamic underlying repeated outbreaks of pandemic disease. See R.G. [START_REF] Wallace | Big Farms Make Big Flu: Dispatches on infectious disease, agribusiness and the nature of science[END_REF][START_REF] Wallace | Dead Epidemimiologists: On the origins of COVID-19[END_REF] for detailed expositions and comprehensive overviews of the literature. Similar matters are more abstractly studied in [START_REF] Wallace | Big Farms Make Big Flu: Dispatches on infectious disease, agribusiness and the nature of science[END_REF] and R. [START_REF] Wallace | Clear-Cutting Disease Control: Capital-led deforestation, public health austerity, and vector-borne infection[END_REF]. More recently, R. Wallace et al. (2022, Ch. 5) provide mathematical models showing how such 'economies', aimed at a draconian -Procrustean -smoothing of workflows in time, space, and genetic structure, are, from the perspectives of public health, actually diseconomies of scale, raising the probability of pandemic pathogen escape to virtual certainty.

Here, we reconsider these matters, focusing on interactions between such 'economies' and the failures of public health systems constrained by the rate at which disease control actions can be delivered in the face of rising rates of disease spread. The stochastic theory of differential equations (Appleby et al. 2008;[START_REF] Mao | Stochastic stabilisation and destabilisation[END_REF]) details how obviating those 'economies' -fragmenting and roughening them -can reimpose public health control of pandemic disease via a kind of sterilization. The argument -at least mathematically -is surprisingly straightforward.

We begin with a reconsideration of 'economies of scale', examine public health entities as control systems, and construct a simple mathematical model based on the Data Rate Theorem linking information and control theory illustrating how carefully-designed policies of 'stochasticization' can drive emerging pathogens to extinction.

'Economies of scale' redux

We suppose there are n basic inverse measures of corporate agribusiness 'economies of scale' as studied in Wallace et al. (2022, Ch.5). These are, in addition to temporal, spatial, and genetic, also across financial, organizational, and 'social' matters. Examples of the latter forms include 'simple' verticalization, illegal cartels suppressing competition and/or governmental regulation, or in some polities, private armies interfacing with local police to intimidate smallholders and indigenous organizations, and so on. The larger index, the 'smoother' the workflow, and the more concentrated the 'corporate socialism' privatizing profits and socializing costs.

In reality, these measures of such 'corporate efficiency' must interact, so that one is confronted with a square n × n matrix, say R, having elements R i,j . A square matrix of order n will have n scalar invariants, that is, n real numbers that remain constant under certain symmetry operations. The first invariant is usually taken as the trace, and the last as ± the determinant.

Given the complexity of modern agribusiness enterprises, n may not be small, particularly across powerful cartel/political networks.

The invariants can be expressed as the coefficients of the characteristic polynomial P(λ)

P(λ) = det(R -λI) = λ n + r 1 λ n-1 + ... + r n-1 λ + r n (1)
where λ is a parameter that is an element of a ring, det is the determinant, and I is the n × n identity matrix. λ may even be taken as the matrix R itself, since square matrices form a ring. We will assume it possible, as in Principal Component Analysis, to project most of the 'variance' -in some appropriate sense -down on to a single scalar function of the n matrix invariants, writing a scalar function as

ρ = ρ(r 1 , r 2 , ...r n ) (2)
Figure 1: A linearized control system near the nonequilbirium steady state. X t+1 is the system's response to the control signal U t and the earlier state X t . U t is the output of the control information source, and W t a 'noise' vector.

3 Public health as control of an inherently unstable system

As Wallace et al. (2022, Ch.5) note at some length, modern nation-states confronted by emerging infection are not simply natural populations of infective, susceptible and 'removed' individuals. Modern states are active cognitive entities that, when confronted by challenge or opportunity, must and do choose one or more actions from the much larger set of what is available to them. Such choice implies a reduction in uncertainty, and any reduction of uncertainty necessarily implies the existence of an information source 'dual' to the cognitive process studied. The argument is direct, unambiguous, and compelling across many modalities and scales (e.g., [START_REF] Wallace | Consciousness: A Mathematical Treatment of the Global Neuronal Workspace Model[END_REF][START_REF] Wallace | Consciousness, crosstalk, and the mereological fallacy: an evolutionary perspective[END_REF][START_REF] Wallace | Computational Psychiatry: A systems biology approach to the epigenetics of mental disorders[END_REF][START_REF] Wallace | Consciousness, Cognition and Crosstalk: The evolutionary exaptation of nonergodic groupoid symmetry-breaking[END_REF].

In the context of an emerging pathogen, a public health entity attempts to impose control on an inherently unstable system -an explosively growing infection. Such matters are the purview of the Data Rate Theorem linking control and information theories [START_REF] Nair | Feedback control under data rate constraints: an overview[END_REF].

Figure 1 outlines a polity's public health system in the context of an emerging pathogen. At time t, the system receives a multidimensional state vector X t representing both disease and institutional crossectional data.

It then produces an output vector at time t + 1 as X t+1 . At time t the system is also affected by a 'noise' vector W t including uncontrolled inputs and a control vector U t representing the cognitive activities of the state. The standard approach [START_REF] Nair | Feedback control under data rate constraints: an overview[END_REF]) examines the dynamics of a linear approximation at the nonequilibrium steady state (nss) written as

X t+1 = AX t + BU t + W t (3)
A and B are fixed matrices characteristic of the system at nss. In the first period of an emerging infection, this system must be seen as inherently -indeed, explosively -unstable. The matrix A can be factored by a similarity transformation into two diagonal submatrices A U and A S , with two zero off-diagonal matrices. A U has eigenvalues ≥ 1 and A S has eigenvalues < 1. The portion of the system under A U is growing explosively, and must be brought under control by the imposition of the control signal vector U t . The Data Rate Theorem states that the rate at which control information is provided, H, must be greater than the rate at which the system generates its own 'topological information' according to the relation

H > log || det(A U )|| ≡ H 0 ( 4
)
where det is the determinant of the unstable subcomponent of A.

In the context of a public health system under siege from concentrated agribusiness actions, this relation becomes more general, and we write it as

H(ρ) > f (ρ)H 0 (5)
H 0 now represents the basic topology of the contact patterns at national, municipal, and social network scales and levels of organization, i.e., the basic backcloth on which infection spreads. What are H(ρ) can, somewhat surprisingly, be approximated using a standard Black-Scholes treatment abducted from financial engineering [START_REF] Black | The pricing of options and corporate liabilities[END_REF], as described in the Mathematical Appendix. In first order

H ≈ κ 1 ρ + κ 2 (6)
Making the same level of approximation for f (ρ) lead to the condition

T ≡ κ 1 ρ + κ 2 κ 3 ρ + κ 4 > H 0 (7)
T will be characterized as the 'control temperature' of the system. At low values of ρ the stability condition is κ 2 /κ 4 > H 0 , and at high values it is κ 1 /κ 3 > H 0 . If κ 2 /κ 4 > κ 1 /κ 3 , then, at some intermediate value of agribusiness 'smoothing' ρ, the stability criterion is violated, triggering uncontrolled pandemic spread, as in figure 2,where κ

1 = 1/4, κ 2 = 8, κ 3 = κ 4 = 1¿

Stochastic stabilization

As Appleby et al. (2008) and [START_REF] Mao | Stochastic stabilisation and destabilisation[END_REF] show at great length, stochastic differential equations of the form Figure 2: The horizontal line is the critical value H 0 . If κ 2 /κ 4 > κ 1 /κ 3 , then at some value of the 'environmental insult' ρ, the system falls below the critical control level H 0 from Eq.( 7) and agribusiness-driven pandemic disease cannot be controlled. Here,

κ 1 = 1/4, κ 2 = 8, κ 3 = κ 4 = 1. dX t = f (X t )dt + g(X t )dB t (8)
where dB t is Brownian white noise, have a peculiar characteristic, namely that any stable/unstable relation dX/dt = f (X(t)) can be destabilized/stabilized by proper choice of g(X t ) in Eq.( 8).

In a fundamental sense, this is the whole ballgame, but it interesting to carry out a more detailed example using Eq.( 7). That is, we will now explore the stochastic stability of Eq.( 7) originally set in the unstable mode as a function of applied 'noise' using the basic model

dρ t = f (ρ t )dt + σh(ρ t )dB t (9)
taking dρ/dt = f (ρ). In the examples below, we will assume an exponential model, i.e., f (ρ) = β -αρ, but proceed here in general.

We calculate the expression for dT t from Eq.( 7) via the Ito Chain Rule (Protter 2005), and then evaluate < dT t >= 0 for the nonequilibrium steady state.

This procedure produces a surprisingly simple result, The first relation can be solved for ρ t in terms of σ, and then plugged in to the expression for T from Eq.( 7), giving T (σ) at nss. The second and third relations illustrate the results of Appleby et al. (2008) and [START_REF] Mao | Stochastic stabilisation and destabilisation[END_REF] in that proper choice of the function h(ρ) will always make ρ = V -1 (σ 2 ) a monotonic decreasing function of σ.

h(ρ t ) 2 - f (ρ t ) σ 2 (ρ t + κ 4 /κ 3 ) = 0 σ 2 = f (ρ t ) h(ρ t ) 2 (ρ t + κ 4 /κ 3 ) ≡ V (ρ t ) ρ t = V -1 (σ 2 ) (10)
As an example, we again assume an exponential model for f (ρ t ) = β -αρ t , and T as in figure 2, with β = 2, α = 1. For h(ρ) = ρ, T (σ) is shown in figure 3, where increasing 'noise' eventually stabilizes the system, raising T above the critical level H 0 .

For h(ρ) = 1, a more complicated pattern emerges in figure 4: above σ = 3/2 the system assumes a non-zero imaginary component, implying a phase transition to instability. Clearly, just how stochasticity is imposed -the precise nature of g(X) in Eq.( 8) -is indeed critical to system dynamics. Some exploration of Eq.( 10) suggests that if, in the second term in Eq.( 9), h(ρ) is positive and monotonic increasing, the system will follow the pattern of figure 3 rather than figure 4. That is, the 'noise' level should be linked to the observed stress index ρ to avoid disruptive phase transitions.

Discussion

Although, from an epidemiological modeling perspective, Eq.( 8) most cleanly encompasses the canonical intervention -use 'noise' as a temperature analog to Figure 4: T (ρ(σ)) for h(ρ) = 1, again based on the Ito Chain Rule expression for < dT t >= 0 as in figure 3. Above σ = 3/2, however, the system assumes a nonzero imaginary component, implying a phase transition to instability. Evidently -and consonant with the arguments surrounding Eq.( 8) -just how stochasticity is imposed is critical to system dynamics under rising stress. carefully disperse and 'sterilize' current agribusiness cartels -Eqs.( 9) and ( 10) and figures 3 and 4 provide an interesting detailed example. Our approach is based on a Data Rate Theorem perspective on public health as an infectious disease control system. That is, as have many others, we argue that prevention and control of mass-fatal pandemic disease is predicated on the imposition of a sterilizing dispersal across current large-scale agribusiness enterprise. Dispersal and reconfiguration of large-scale global agribusiness in time, space, and genetic structure -along with escape from the excesses of 'vertical integration' and associated agricartel political corruption and violence -will be required to avert repeated pandemic outbreaks likely to rapidly reduce human populations to a 'global carrying capacity' of perhaps 500 million stunned survivors.

Figures 3 and4, however, also suggest that the breakup of current agricartels, even beyond the obvious political difficulties, is likely to be a subtle matter, and that the exact form in which sterilizing stochasticity is imposed can be of singular importance. The critical transition shown in figure 4 implies that policies of agricartel dispersal must be carefully designed, continuously monitored and revised as data indicate to ensure proper public health outcomes.

Mathematical Appendix

A Black-Scholes model

We look at H(ρ) as the control information rate 'cost' of stability at the integrated environmental insult ρ. To determine the mathematical form of H(ρ) under conditions of volatility i.e., variability proportional to a signal, we must first model the variability of ρ, most simply taken as

dρ t = g(t, ρ t )dt + bρ t dB t (11)
Here, dB t is white noise and -counterintuitively -the function g(t, ρ) will fall out of the calculation on the assumption of certain regularities.

H(ρ t , t) is the minimum needed incoming rate of control information under the Data Rate Theorem. Expand H in ρ using the Ito chain rule [START_REF] Protter | Stochastic Integration and Differential Equations[END_REF]: It is now possible to define a Legendre transform, L, of the rate H, by convention having the form

dH t =
L = -H + ρ∂H/∂ρ (13)
H is an information index, a free energy measure in the sense of Feynman (2000), so that L is a classic entropy measure.

We make an approximation, replacing dX with ∆X and applying Eq.( 12), so that ∆L = (-∂H/∂t -1 2 b 2 ρ 2 ∂ 2 H/∂ρ 2 )∆t ( 14)

According to the classical Black-Scholes model [START_REF] Black | The pricing of options and corporate liabilities[END_REF], the terms in g and dB t 'cancel out', and white noise has been subsumed into the Ito correction factor, a regularity assumption making this an exactly solvable but highly approximate model.

The conventional Black-Scholes calculation takes ∆L/∆T ∝ L. At nonequilibrium steady state, by some contrast, we can assume ∆L/∆t = ∂H/∂t = 0, giving

- 1 2 b 2 ρ 2 ∂ 2 H/∂ρ 2 = 0 (15) so that H = κ 1 ρ + κ 2 (16)
The κ i will be nonnegative constants.

Figure 3 :

 3 Figure 3: T (ρ(σ)) for h(ρ) = ρ, based on the relation < dT t >= 0 via application of the Ito Chain Rule from Eq.(10), using the κ-values of figure 2, taking β = 2, α = 1. At low σ the stability criterion is violated, but becomes enabled at sufficiently high 'noise' σ, a classic example of noise-induced stabilization.

  [∂H/∂t + g(ρ t , t)∂H/∂ρ + 1 2 b 2 ρ 2 t ∂ 2 H/∂ρ 2 ]dt +[bρ t ∂H/∂ρ]dB t(12)

Data availability

There are no external data sets associated with this paper.