Katia Ait-Ameur
email: aitameur.katia@gmail.com

Yvon Maday
email: maday@ann.jussieu.fr

Marc Tajchman
email: marc.tajchman@cea.fr

Multi-step variant of parareal algorithm

come   L'archive ouverte pluridisciplinaire

Introduction

In the field of nuclear energy, computations of complex two-phase flows are required for the design and safety studies of nuclear reactors. System codes are dedicated to the thermal-hydraulic analysis of nuclear reactors at system scale by simulating the whole reactor. We are here interested in the Cathare code developed by CEA, [START_REF] Bestion | The physical closure laws in the cathare code[END_REF]. Typical cases involve up to a million of numerical time iterations, computing the approximate solution during long physical simulation times. A space domain decomposition method has already been implemented. To improve the response time, we will consider a strategy of time domain decomposition, based on the parareal method [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF]. The Cathare time discretization is based on a multi-step time scheme (see [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]). In this paper, we derive a strategy to adapt the parareal algorithm to multistep schemes that is not implemented in the code. The paper is organized as follows.

In Section 2, we recall the classical version of the parareal algorithm, and then detail the variant that allows us to use multi-step time schemes for the fine solvers. A couple of remarks on the algorithm will be discussed. The numerical convergence for a simplified test case is shown in Section 3 on a Dahlquist test equation followed by numerical results on an advection-diffusion equation and on an industrial test case with an application on the Cathare code.

Parareal algorithm and the multi-step variant

After the discretization of a PDE in space, we obtain an ODE system of the form:

∂u ∂t + A(t, u) = 0, t ∈ [0, T], u(t = 0) = u 0 ( 1
)
where A : R × R N → R N , and N denotes the number of degrees of freedom. We here recall the classical parareal algorithm as initially proposed in [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF], [START_REF] Baffico | Parallel-in-time moleculardynamics simulations[END_REF], [START_REF] Bal | A "parareal" time discretization for non-linear pde's with application to the pricing of an american put[END_REF].

Let G and F be two propagators such that, for any given t ∈ [0, T], s ∈ [0, Tt] and any function w in a Banach space, G(t, s, w) (respectively F(t, s, w)) takes w as an initial value at time t and propagates it at time t + s. The full time interval is divided into N c sub-intervals [T n , T n+1 ] of size ∆T that will each be assigned to a processor. The algorithm is defined using two propagation operators:

• G(T n , ∆T, u n ) computes a coarse approximation of u(T n+1 ) with initial condition u(T n ) u n (low computational cost) • F(T n , ∆T, u n ) computes a more accurate approximation of u(T n+1 ) with initial condition u(T n ) u n (high computational cost)

Starting from a coarse approximation u n 0 at times T 0 , T 1 , • • • , T N c , obtained using G, the parareal algorithm performs for k = 0, 1, • • • the following iteration:

u n+1 k+1 = G(T n , ∆T, u n k+1 ) + F(T n , ∆T, u n k ) -G(T n , ∆T, u n k )
In the parareal algorithm, the value u(T n ) is approximated by u n k at each iteration k with an accuracy that tends rapidly to the one achieved by the fine solver, when k increases. The coarse approximation G can be chosen much less expensive than the fine solver F by the use of a scheme with a much larger time step (even δT = ∆T) δT δt (time step of the fine solver) or by using a reduced model. All the fine propagations are made in parallel over the time windows and the coarse propagations are computed in a sequential way but have a low computational cost. We refer to [START_REF] Maday | An adaptive parareal algorithm[END_REF] about the parallel efficiency of parareal and a recent work offering a new formulation of the algorithm to improve the parallel efficiency of the original one. The main convergence properties were studied in [START_REF] Gander | Analysis of the parareal time-parallel time-integration method[END_REF] and stability analysis was made in [START_REF] Staff | Stability of the parareal algorithm[END_REF], [START_REF]Parallelization in time of (stochastic) ordinary differential equations[END_REF]. In the sequel, we will consider the case that the coarse solver is based on a one-step time scheme and the fine solver on a two-step time scheme. Hence we will use the following notation for the fine solver that takes two initial values: F(t, s, x, y), for t ∈ [0, T], s ∈ [0, Tt[ and x, y in a Banach space.

Example

If one solves (1) with a multi-step time scheme as fine propagator F like the second-order BDF method:

3 2 u j+1 -2u j + 1 2 u j-1 = -δt A(u j+1 , t j+1 ), j = 1, • • • , N f , t j+1 -t j = δt
Here the fine solver reads as: u j+1 = F(t j , δt, u j-1 , u j ). Now, we apply the parareal algorithm with a coarse grid: T 0 , • • • , T N c where:

T n+1 -T n = ∆T = N f δt.
Then we can write:

u(T n + jδt) u n, j , j = 1, • • • , N f , n = 1, • • • , N c .
In order to perform the fine propagation, in a given time window [T n , T n+1 ], we only need the local initial condition u n k and a consistent approximation of u(T nδt).

In [START_REF] Audouze | Symplectic multi-time step parareal algorithms applied to molecular dynamics[END_REF], the authors propose a consistent approximation in the context of the simulation of molecular dynamics. The proposed method was linked to the nature of the model and the symplectic character of their algorithm is shown, which is an important property to verify for molecular dynamics. In the context of our application to the thermalhydraulic code Cathare, we want to derive a multi-step variant of parareal that will not be intrusive in the software. We seek a consistent approximation of u(T nδt). The only fine trajectory at our disposal is

F(T n-1 , ∆T, u n-2, N f -1 k , u n-1 k ).
Its final value at T n is:

F(T n-1 , ∆T, u n-2, N f -1 k , u n-1 k )(T n
) from which we compute u n k+1 by the parareal correction. Hence, we translate the solution:

F(T n-1 , ∆T -δt, u n-2, N f -1 k , u n-1 k )(T n -δt) by the same correction: u n k+1 -F(T n-1 , ∆T, u n-2, N f -1 k , u n-1 k
) and obtain the so called consistent approximation u n-1, N f -1 k+1 to initialize the fine propagation in [T n , T n+1 ]. We now detail our algorithm:

                     u n+1 0 = G(T n , ∆T, u n 0 ), 0 ≤ n ≤ N -1 u n+1 k+1 = G(T n , ∆T, u n k+1 ) + F(T n , ∆T, u n-1, N f -1 k , u n k ) -G(T n , ∆T, u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n, N f -1 k+1 = F(T n , ∆T -δt, u n-1, N f -1 k , u n k ) + u n+1 k+1 -F(T n , ∆T, u n-1, N f -1 k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 (2)
Another option to treat this issue is to use a one-step time scheme to initialize the fine computation or to make one iteration with a second-order Runge Kutta method. We will see from the numerical results that these choices modify the fine scheme and prevent the parareal algorithm to converge to the fine solution: even after N c iterations (where N c is the number of time windows), the parareal algorithm does not converge to the fine solution with the machine precision but the parareal error stagnates around 10 -6 . This method adds consistency with the fine scheme. Also, this strategy can be applied to multi-step time schemes involving several fine time steps preceding the time T n by applying the same correction to terms taking the form:

u n, N f -i k+1 , i = 1, • • • , I.
The convergence analysis will be shown in a forthcoming work.

Numerical results

We now show some numerical experiments, first for an ordinary differential equation, the test equation of Dahlquist, then for a partial differential equation, namely the advection-diffusion equation and finally on an industrial application with the Cathare code.

Dahlquist equation

We first use the Dahlquist test equation :

u (t) = λu(t), t ∈ (0, T), with u(0) = 1,
discretized by a second-order BDF method. With λ = -1, T = 5, ∆T = T/50, which correspond to 50 processors, and δt = T/5000, we obtain the convergence curves shown in 3.1. Here, the fine solver is based on a two-step time scheme where the computation of the solution u n, j+1 at time T n + ( j + 1)δt needs the knowledge of the solutions u n, j and u n, j-1 at times T n + jδt and T n + ( j -1)δt, respectively. We use the multi-step variant of parareal (2) to initialize the fine solver in each time window, starting from the parareal iteration k ≥ 2. At the parareal iteration k = 1, we use a Backward Euler method to initialize the fine solver since we don't have the fine solution yet. The coarse solver is based on a one-step time scheme, namely the Backward Euler method. We plot the relative error in L ∞ (0, T) between the fine solution computed in a sequential way and the parareal solution as a function of iteration number for the classical parareal algorithm where the Backward Euler method is used at each iteration for the initialization of the fine solver (circles), and for the multi-step variant of the parareal algorithm that we introduced in the previous section (squares). We see in the Fig. 1 (without multi-step) that starting from the fourth parareal iteration the error stagnates around 10 -6 without recovering the fine solution at the machine precision, even after 50 iterations. On the other hand, we see in Fig. 1 (with multi-step) that the error continues to decrease after the fourth parareal iteration until reaching the machine precision at the eleventh iteration. In this case, if we don't use the multi-step variant of parareal, we loose one of the well known property of parareal: to recover the fine solution at the machine precision after N c iterations.

Advection-diffusion equation

We now study the behavior of the multi-step parareal (2) applied to the advectiondiffusion equation:

         ∂u ∂t -∂ 2 u ∂x 2 + ∂u ∂x = 0, (x, t) ∈ (0, 2π) × (0, T) u(x, 0) = u 0 (x), x ∈ (0, 2π) u(x, t) = u(x + 2π, t), t ∈ (0, T) (3)
We have chosen a spectral Fourier approximation in space (truncated series with N = 16) and a second-order BDF method in time for a propagation over [0, 2]. The parareal in time algorithm is implemented with ∆T = 0.1 and δt = 10 -3 . We have chosen the following initial condition:

u 0 (x) = l=N /2 l=-N /2+1
ûl e ilx , with ûl = sgn(l) |l | p .

We choose p = 4, hence the initial condition is sufficiently regular.

The coarse and fine solvers are the same as in the previous numerical example and we use the same initialization of the fine solver at the first parareal iteration. We plot the following error as a function of iteration number:

E n k = max n | |U s eq (T n )-U n k | | L 2 ((0,2π)) max n | |U s eq (T n ) | | L 2 ((0,2π))
, where U seq is the fine solution computed in a sequential way. In Fig. 2 (circles), we observe a similar behavior as in the previous case, the error stagnates around 10 -6 after the fifth parareal iteration without reaching the machine precision when the number of iterations is equal to the number of time windows (20 in this case). In Fig. 2 (squares), using the multi-step variant, the error continues to decrease after the fifth iteration until reaching the machine precision around 10 -16 at iteration number 20. this study show that this variant converges numerically on different examples: the two simple test equations allow to see clearly the advantage of our strategy. The application on an industrial code shows its efficiency on a more realistic test case without being intrusive in the software. The convergence analysis of this algorithm will be the subject of a forthcoming paper. The extension of this method to the use of a multi-step scheme in the coarse solver will be also investigated. This point was treated in [START_REF] Audouze | Symplectic multi-time step parareal algorithms applied to molecular dynamics[END_REF] and a similar strategy will be studied.

Fig. 1

 1 Fig. 1 Convergence of the multi-step parareal for the Dahlquist test equation

Fig. 2

 2 Fig. 2 Convergence of the parareal for the advection-diffusion equation

Fig. 3

 3 Fig. 3 Multi-step parareal for an industrial application with Cathare code

Acknowledgements

This research is partially supported by ANR project CINE-PARA (ANR-15-CE23-0019).

Application to the Cathare code

The Cathare code simulates two-phase flows at a macroscopic scale and the model used is the six-equation two-fluid model ( [START_REF] Drew | Theory of multicomponent fluids[END_REF], [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF], [START_REF] Ndjinga | Influence of interfacial pressure on the hyperbolicity of the two-fluid model[END_REF]) that considers a set of balance laws (mass, momentum and energy) for each phase liquid and vapor. The unknowns are the volume fraction α k ∈ [0, 1], the enthalpies H k , the velocity u k of each phase and the pressure p. The density ρ k ≥ 0 is computed with equations of state (k = l, g). The Cathare scheme is based on a finite volume method on a staggered grid (MAC scheme) and on a two-step time scheme. Here, we write the time discretization of the Cathare scheme:

This choice of time discretization was made for stability purpose.

Here we apply the parareal algorithm to the solution of an oscillating manometer. This test case is proposed in [START_REF] Hewitt | [END_REF] for system codes to test the ability of each numerical scheme to preserve system mass and to retain the gas-liquid interface.

We have used coarse and fine solvers such that δt = 10 -5 for F and ∆T = 10δt for G. All calculations have been evaluated with a stopping criterion where the tolerance is fixed to the precision of the numerical scheme, = 5 • 10 -2 . In order to perform the fine propagation, in a given time window [T n , T n+1 ], at the first parareal iteration we need to choose a different consistent approximation of u(T nδt), since we have not used the fine solver yet. In the context of the application to the Cathare code, we choose a non intrusive initialization by imposing

In Fig. 3, we plot the evolution of the relative error in L 2 norm between the parareal solution and the fine one accross the time. This result illustrates that the multistep variant (2) of the parareal algorithm effectively converges when applied to the problem of the oscillating manometer.

Conclusion

We have built a new variant of parareal algorithm allowing to overcome the issue of initializing the fine solver when the time scheme involve the numerical solution at times preceding the local initial condition in a given time window. The results of


David Cohen-Steiner

André Lieutier

Julien Vuillamy

Regular triangulations as lexicographic optimal chains

We introduce a total order on n-simplices in the n-Euclidean space for which the support of the lexicographic-minimal chain with the convex hull boundary as boundary constraint is precisely the n-dimensional Delaunay triangulation, or in a more general setting, the regular triangulation of a set of weighted points. This new characterization of regular and Delaunay triangulations is motivated by its possible generalization to submanifold triangulations as well as the recent development of polynomial-time triangulation algorithms taking advantage of this order.

Introduction

Algorithms computing the Delaunay triangulation and its variants [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF] can be described as a set of combinatorial operations, evaluating geometric predicates and acting on algebraic quantities such as coordinates or equations. However, the Delaunay triangulation, or more generally regular triangulations, also have a variational formulation, introduced in [START_REF] Chen | Optimal delaunay triangulations[END_REF] and further studied in [START_REF] Alliez | Variational tetrahedral meshing[END_REF][START_REF] Chen | Efficient mesh optimization schemes based on optimal delaunay triangulations[END_REF], characterizing them as solutions of a linear programming problem on triangulations. We show that this formulation can in fact be extended to the space of simplicial chains (Proposition 4.3).

The main part of this work consists in proving this variational formulation can be transposed in terms of lexicographic minimum (Theorem 3.1), for which polynomial time algorithms have recently been devised [START_REF] Cohen-Steiner | Lexicographic optimal homologous chains and applications to point cloud triangulations[END_REF]. The proof of the theorem may be of interest in itself by offering an unusual point of view on regular triangulations.

One motivation for this work is that the same lexicographic order may be applied to the meshing of submanifold of Euclidean space. To be usable in this context, a such order must be isometry invariant, so that it can be evaluated from edge lengths only. It has been shown indeed that minimum chains in their homology class, for this lexicographic order, provide triangulations of well sampled (with respect to the reach) smooth 2-submanifolds of Euclidean n-space [START_REF] Cohen-Steiner | Lexicographic optimal chains and manifold triangulations[END_REF], akin to the tangential Delaunay complex [START_REF] Boissonnat | Manifold reconstruction using tangential delaunay complexes[END_REF]. In practice, this order enables efficient algorithms providing minimal solutions that inherit the optimality properties of 2-dimensional Delaunay triangulations and create pertinent and convincing meshes for surface reconstruction, in particular in the case of noisy point clouds with non uniform densities and outliers [START_REF] Cohen-Steiner | Lexicographic optimal homologous chains and applications to point cloud triangulations[END_REF].

Conventions and notations

A k-simplex σ being a set of k + 1 vertices in ambient space R n , we allow ourselves to use set theoretic operators on simplices. For example, τ ⊂ σ means that τ is a face of σ and σ 1 ∪ σ 2 is the join of σ 1 and σ 2 . |σ| denotes the underlying space of the simplex σ, i.e. the convex hull of its vertices which, thanks to the generic condition 3.1, are affinely independent.

Homology coefficients are implicitly in Z 2 = Z/2Z, i.e. integers modulo 2 so that 1 = -1. The vector space of k-chains over K is denoted C k (K) in place of C k (K, Z 2 ). In this context we can allow ourselves to see chains as sets of simplices. For example, for Γ ∈ C k (K) and σ is a k-simplex in K, we write interchangeably σ ∈ Γ and Γ(σ) = 1. Similarly, for Γ 1 , Γ 2 ∈ C k (K), we use interchangeably vector and set theoretic operators:

Γ 1 + Γ 2 = Γ 1 -Γ 2 = (Γ 1 ∪ Γ 2 ) \ (Γ 1 ∩ Γ 2 ).
For a chain Γ in C k (K), we denote by |Γ| the support of Γ, which is the sub-complex of K made of all k-simplices in Γ together with all their faces.

A total order ≤ on k-simplices induces a lexicographic order on k-chains with coefficients in Z 2 as follows: Definition 2.1 (Lexicographic Order on chains). Assume there is a total order ≤ on the set of k-simplices of K, defining the max on sets of k-simplices. For Γ 1 , Γ 2 ∈ C k (K):

Γ 1 lex Γ 2 ⇐⇒ def.      Γ 1 = Γ 2 or max {σ ∈ Γ 1 + Γ 2 } ∈ Γ 2

Main result

We consider a set P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n × R of weighted points in n-dimensional Euclidean space. A weighted point (P, 0) is seen as a usual point P ∈ R n , while, when µ > 0, it is associated to the sphere centered at P with radius r = √ µ. The convex hull CH(P) denotes the convex hull in R n of the set of points: CH(P) = CH ({P 1 , . . . , P N }).

The n-dimensional full complex over P, denoted K P , is the simplicial complex made of all possible simplices up to dimension n with vertices in P. The aim of this paper is to prove the following: When all the weights are zero (∀i, µ i = 0), the regular triangulation (Definition 4.1) is the Delaunay triangulation.

Theorem 3.1. Let P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n × R, with N ≥ n +
The relation lex among n-chains is the lexicographic order defined according to Definition 2.1, where the total order on n-simplices is given at the end of Section 3.3. The general position assumption is a generic condition formalized in Condition 3.2.

Obviously, replacing in Theorem 3.1 the full complex K P by a complex containing the regular triangulation would again give this triangulation as a minimum.

Outline of the proof.

The main argument of the proof is given in Section 6.

Statement of Proposition 4.3 is the same as Theorem 3.1 except for the order along which the minimum is taken. In Proposition 4.3, the minimization is meant for the preorder (which is in fact generically an order) p , induced by a weighted L 1 norm • (p) on chains. All the proof consists then in showing that, while the two orders differ, they share the same minimum under the theorem constraints. The bounding weight of a simplex is a generalization, for weighted points, of the radius of the smallest ball enclosing the simplex. It is the dominant quantity in the definition of the order on simplices which induces the lexicographic order lex on chains. The proof then proceeds in 3 main steps:

1) For p large enough, the weight w p (σ) p of any single simplex σ in the • (p) norm is larger than the sum of all weights of all simplices with smaller bounding weight than σ (Lemma 4.5 in Section 4.1). This fact allows to focus on the link of a single simplex, the one with the largest bounding weight µ = B defined in (30), for which some simplices in the minimum of lex and p would differ as explained in Section 6.

2) We introduce (Sections 4.2 and 4.3) an extension of the classical lifted paraboloid construction that allows to see the simplicial structure of the link of a simplex τ in the Delaunay triangulation as a convex cone. In this representation, we study the subcomplex of this link corresponding to cofaces of τ with same bounding weight as τ . We show that this subcomplex is isomorphic to a set of bounded and visible faces of a convex polytope. It is therefore a simplicial ball (Lemma 4.13 in Section 4.3).

3) By induction on the dimension of convex cones and convex polytopes, one shows that this bounded subcomplex of the boundary of a convex polytope visible from the origin can be expressed as a the minimum under boundary constraint (Lemma 5.4 in Section 5) for another lexicographic order.

This second lexicographic order corresponds, on the restriction of the link of τ to simplices with the same bounding weight, to the first lexicographic order on corresponding full dimensional simplices in the star of τ (Lemmas 6.1 and 6.2 in Section 6).

Weighted points and weighted distances

We follow the terminology, notations and conventions from Section 4.4 of [START_REF] Boissonnat | Geometric and topological inference[END_REF]. Definition 3.2 (Section 4.4 in [START_REF] Boissonnat | Geometric and topological inference[END_REF]). Given two weighted points (P 1 , µ 1 ), (P 2 , µ 2 ) ∈ R n × R their weighted distance is defined as:

D ((P 1 , µ 1 ), (P 2 , µ 2 )) = def. (P 1 -P 2 ) 2 -µ 1 -µ 2
We say (P 1 , µ 1 ) and (P 2 , µ 2 ) are orthogonal if D ((P 1 , µ 1 ), (P 2 , µ 2 )) = 0.

We introduce a first generic condition: Condition 3.1. We say that P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n × R satisfies the first generic condition if no (n + 1) points {P i1 , . . . , P in+1 } lie on a same (n -1)-dimensional affine space.

We now define the generalization of the circumsphere and the smallest enclosing ball for sets of weighted points. Definition 3.3. Assume that P satisfies Condition 3.1. Given a k-simplex σ ⊂ P with 0 ≤ k ≤ n, the generalized circumsphere and smallest enclosing ball of σ are the weighted points (P C , µ C )(σ) and (P B , µ B )(σ) respectively defined as:

µ C (σ) = def. min {µ ∈ R, ∃P ∈ R n , ∀(P i , µ i ) ∈ σ, D ((P, µ), (P i , µ i )) = 0}
(1)

µ B (σ) = def. min {µ ∈ R, ∃P ∈ R n , ∀(P i , µ i ) ∈ σ, D ((P, µ), (P i , µ i )) ≤ 0} (2)
P C (σ) and P B (σ) are respectively the unique points P that realize the minimum in Equations (1) and [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF]. The weights µ C (σ) and µ B (σ) are called respectively circumweight and bounding weight of σ.When ∀i, µ i = 0, they correspond respectively to the square of the circumradius and the square of the radius of the smallest ball enclosing σ.

The set {(P, µ), ∀(P i , µ i ) ∈ σ, D ((P, µ), (P i , µ i )) = 0}, on which the first arg min is taken, is not empty, thanks to the generic condition 3.1.

Lemma 3.4 (Proof in Appendix B).

For any k-simplex σ ∈ K P , one has P B (σ) ∈ |σ|. Condition 3.2. We say that P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n × R is in general position if it satisfies the first generic condition 3.1 and if, for a k-simplex σ and a k -simplex σ in K P with 2 ≤ k, k ≤ n, one has:

µ C (σ) = µ C (σ ) ⇒ σ = σ (3)
From now on, we assume P to be in general position i.e. it satisfies Condition 3.2.

Lemma 3.5 (Proof in Appendix C

). Under generic condition 3.2 on P, for any simplex σ, there exists a unique inclusion minimal face Θ(σ) of σ such that (P B , µ B )(σ) = (P C , µ C )(Θ(σ)).

Moreover one has (P C , µ C )(Θ(σ)) = (P B , µ B )(Θ(σ)).

Figure 1 illustrates the possibilities for Θ(σ) in the case n = 3 and zero weights.

Regular triangulation order on simplices

For a n-simplex σ, we define a (k + dim(Θ(σ)))-dimensional face Θ k (σ) as follows. For k = 0, Θ 0 (σ) = Θ(σ), where Θ(σ) is defined in Lemma 3.5.

For k > 0, Θ k (σ) is the (dim(Θ k-1 (σ)) + 1)-dimensional coface of Θ k-1 (σ) with minimal circumradius: Θ k (σ) = arg min Θ k-1 (σ)⊂τ ⊂σ dim(τ )=dim(Θ k-1 (σ))+1 µ C (τ ) (4)
and µ k (σ) is the circumweight of Θ k (σ):

µ k (σ) = µ C (Θ k (σ)). In particular, µ 0 (σ) = µ C (Θ(σ)) = µ B (σ) (by Lemma 3.5) and if k = dim(σ) -dim(Θ(σ)) then µ k (σ) = µ C (σ).
Observe that, thanks to generic condition 3.2 and Lemma 3.5, one has, for two n-simplices σ 1 , σ 2 : µ B (σ 1 ) = µ B (σ 2 ) ⇒ Θ(σ 1 ) = Θ(σ 2 ) and therefore, if µ B (σ 1 ) = µ B (σ 2 ), µ k (σ 1 ) and µ k (σ 2 ) are defined for the same range of values of k.

We define the following order relation on n-simplices (recall that µ 0 (σ) = µ B (σ)):

σ 1 ≤ σ 2 ⇐⇒ def. σ 1 = σ 2 or          µ 0 (σ 1 ) < µ 0 (σ 2 ) or ∃k ≥ 1, µ k (σ 1 ) > µ k (σ 2 ) and ∀j, 0 ≤ j < k, µ j (σ 1 ) = µ j (σ 2 ) (5)
One can check that when P is in general position, the relation ≤ is a total order. For example, when n = 2 and the weights are zero, this order on triangles consists in first comparing the radii of the smallest circles enclosing the triangles T i , i = 1, 2, whose squares are R B (T i ) 2 = µ B (T i ) = µ 0 (T i ). This is generically enough for acute triangles, but not for obtuse triangles that could generically share their longest edge. In this case the tie is broken by comparing in reverse order the circumradii, whose squares are R

C (T i ) 2 = µ C (T i ) = µ 1 (T i ).
Following Definition 2.1 the order ≤ on n-simplices induces a lexicographic order lex on the n-chains of K P . Observation 3.6. From the definition of the order on simplices, the lexicographic minimum is invariant under a global translation of the weights by a common shift s: ∀i, µ i ← µ i + s. (as explicitly explained in observation 7.1 in Appendix A). The same holds for regular triangulations. Therefore proving Theorem 3.1 for non positive weights is enough to extend it to any weights.

Regular triangulations of weighted points

We recall now the definition of a regular triangulation over a set of weighted points. Regular triangulation can alternatively be defined as the dual of the Laguerre (or power) diagram of a set of weighted points. We use here a generalization of the empty sphere property of Delaunay triangulations. Definition 4.1 (Lemma 4.5 in [START_REF] Boissonnat | Geometric and topological inference[END_REF]). A regular triangulation T of the set of weighted points P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n ×R, N ≥ n+1, is a triangulation of the convex hull of {P 1 , . . . , P N } taking its vertices in {P 1 , . . . , P N } such that for any simplex σ ∈ T , if (P C (σ), µ C (σ)) is the generalized circumsphere of σ, then:

(P i , µ i ) ∈ P \ σ ⇒ D ((P C (σ), µ C (σ)), (P i , µ i )) > 0

Lift of weighted points and p-norms

Given a weighted point (P, µ) ∈ R n × R, its lift with respect to an implicit origin O ∈ R n , denoted by lift(P, µ), is a point in R n × R given by: lift(P, µ) = def.

((P -O), (P -O) 2 -µ)
Similarly to Delaunay triangulations, it is a well known fact that simplices of the regular triangulation of P are in one-to-one correspondence with the lower convex hull of lift(P):

Proposition 4.2. A simplex σ is in the regular triangulation of P if and only if lift(σ) is a simplex on the lower convex hull of lift(P).

Based on this lifted paraboloid formulation, the idea of variational formulation for Delaunay triangulations has emerged [START_REF] Chen | Optimal delaunay triangulations[END_REF]. This idea has been exploited further in order to optimize triangulations in [START_REF] Chen | Efficient mesh optimization schemes based on optimal delaunay triangulations[END_REF][START_REF] Alliez | Variational tetrahedral meshing[END_REF]. We follow here the same idea but the variational formulation, while using the same criterion, is applied on the linear space of chains, which can be seen as a superset of the space of triangulations.

We define a function on the convex hull of a k-simplex f σ : |σ| → R where σ = {(P 0 , µ 0 ), . . . , (P k , µ k )} as the difference between the linear interpolation of the height of the lifted vertices and the function x → (x -O) 2 . More precisely, for a point x ∈ |σ| with barycentric coordinates λ i ≥ 0, i λ i = 1, we have x = i λ i P i and:

f σ : x → f σ (x) = def. i λ i ((P i -O) 2 -µ i ) -(x -O) 2 (6)
A short computation shows that the function f σ , expressed in terms of barycentric coordinates, is invariant by isometry (translation, rotation or symmetry on σ). In particular f σ (x) does not depend on the origin O of the lift. It follows from Proposition 4.2 that, if σ reg is a simplex containing x in the regular triangulation of P, for any other simplex σ containing x with vertices in P:

f σreg (x) ≤ f σ (x) (7)
In the particular case where all weights are non positive, that is ∀i, µ i ≤ 0, the convexity of x → x 2 says that the expression of f σ (x) in ( 6) is never negative and in this case [START_REF] Cohen-Steiner | Lexicographic optimal chains and manifold triangulations[END_REF] implies that defining the weight w p of a n-simplex σ as:

w p (σ) = def. f σ p = |σ| f σ (x) p dx 1 p (8)
allows to characterize the regular triangulation as the one induced by the chain Γ reg that, among all chains with boundary β P , minimizes:

Γ → Γ (p) = def. σ |Γ(σ)|w p (σ) p (9)
In this last equation, the notation |Γ(σ)| instead of Γ(σ) is there since Γ(σ) ∈ Z 2 and the sum is in R: Γ(σ) → |Γ(σ)| convert coefficients in Z 2 into binary real numbers in {0, 1}. Formally, we have the following Proposition 4.3 that characterizes regular triangulations as a linear programming problem over Z 2 . Proposition 4.3 (Proof in Appendix D). Let P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n × R, with N ≥ n + 1 be in general position with non positive weights, and let K P be the n-dimensional full complex over P. Denote by β P ∈ C n-1 (K P ) the (n -1)-chain made of simplices belonging to the boundary of CH(P). For any p ∈ [1, ∞), if:

Γ reg = arg min Γ∈Cn (K P ) ∂Γ=β P Γ (p)
then the simplicial complex |Γ reg | support of Γ reg is the regular triangulation of P.

Lemma 4.4 (Proof in Appendix E). One has:

sup x∈|σ| f σ (x) = µ B (σ) ( 10
)
The following is an immediate consequence used in the proof of Lemma 4.5:

lim p→∞ w p (σ) = w ∞ (σ) = f σ ∞ = sup x∈|σ| f σ (x) = µ B (σ) (11) Lemma 4.5 (Proof in Appendix F). Let P = {(P 1 , µ 1 ), . . . , (P N , µ N )} ⊂ R n × R, with N ≥ n + 1,
be in general position with non positive weights. Let K P be the corresponding n-dimensional full complex. For p large enough, the weight w p (σ) p of any n-simplex σ ∈ K P , is larger than the sum of all n-simplices in K P with smaller bounding weight µ B . In other words, if

K [n] P is the set of n-simplices in K P : ∃p , ∀p ≥ p , ∀σ ∈ K [n] P , w p (σ) p > τ ∈K [n] P ,µ B (τ )<µ B (σ) w p (τ ) p
As explained in the proof of the main theorem of Section 6, Lemma 4.5 allows us to focus on the link of a single simplex τ . However, before that, we need to introduce geometrical constructions that give an explicit representation of this link (Sections 4.2 and 4.3).

Projection on the bisector of a simplex

We denote by bis τ the (n -k)-dimensional affine space bisector of τ , formally defined as:

bis τ = def. {x ∈ R n , ∀v 1 , v 2 ∈ τ, D ((x, 0), v 1 ) = D ((x, 0), v 2 )} (12)
In the particular case where dim(τ ) = 0, one has bis τ = R n . Let x → π bisτ (x) and x → d(x, bis τ ) = d(x, π bisτ (x)) denote respectively the orthogonal projection on and the minimal distance to bis τ . We define a projection π τ : P → bis τ × R as follows:

π τ (P, µ) = def.
π bisτ (P ), µ -d(P, bis τ ) 2 (13) Figure 2 illustrates π τ for ambient dimension 3 and dim(τ ) = 1.

Let

o τ = P C (τ ) ∈ bis τ denote the (generalized) circumcenter of τ . If (P i , µ i ) ∈ τ , then D ((o τ , µ C (τ )), (P i , µ i )) = 0. Since o τ = π bisτ (P i ) we have (P i -o τ ) 2 -µ i -µ C (τ ) = d(P, bis τ ) 2 -µ i -µ C (τ ) = 0. It follows that if we denote µ(π τ (P i , µ i )) the weight of π τ (P i , µ i ), one has µ(π τ (P i , µ i )) = µ i -d(P, bis τ ) 2 = -µ C (τ )
. Therefore, π τ sends all vertices of τ to a single weighted point:

(P i , µ i ) ∈ τ ⇒ π τ (P i , µ i ) = (o τ , -µ C (τ )) (14)

Polytope and shadow associated to a link in a regular triangulation

In this section, we study the link of a k-simplex τ in the regular triangulation of P that satisfies:

(P B , µ B )(τ ) = (P C , µ C )(τ ) (15)
We know from Proposition 4.2 that the link of τ in the regular triangulation is isomorphic to the link of lift(π τ (τ )) on the boundary of the lower convex hull of lift(π τ (P)). We consider the lift with the origin at o τ , in other words, the image of a vertex (P, µ) ∈ P is:

Φ τ (P, µ) = def.
lift (π τ (P, µ))

= π bisτ (P ) -o τ , (π bisτ (P ) -o τ ) 2 -µ + d(P, bis τ ) 2
Observe that: Φ τ (τ ) = {(0, µ C (τ )} Call P τ the set of weighted points in (P i , µ i ) ∈ P such that τ ∪ (P i , µ i ) has same bounding weight as τ . In the case of the Delaunay triangulation, P τ corresponds to the set of points in P inside the smallest ball enclosing τ . Formally, using the assumption (15):

P τ = def. {(P i , µ i ) ∈ P \ τ, D ((P C (τ ), µ C (τ )), (P i , µ i )) < 0} (16)
Denote by K τ the (n-k-1)-dimensional simplicial complex made of all up to dimension (n-k-1) simplices over vertices in P τ . Observe that:

D ((P C (τ ), µ C (τ )), (P i , µ i )) < 0 ⇐⇒ (π bisτ (P i ) -o τ ) 2 -µ i + d (P i , bis τ ) 2 < µ C (τ ) (17)
Denote by Height(lift((P, µ)) the height of the lift of a point (P, µ), defined as the last coordinate of the lift, so that:

Height(Φ τ (P, µ)) = (π bisτ (P ) -o τ ) 2 -µ + d(P, bis τ ) 2 (18)
Since under our generic conditions we have π bisτ (P ) -o τ = 0, (16), ( 17) and (18) imply that ∃v ∈ P τ ⇒ Height(Φ τ (v, µ)) > 0 and (17) can be rephrased as: Observation 4.8. A vertex belongs to P τ if and only if the height of its image by Φ τ is strictly less than µ C (τ ) > 0:

v ∈ P τ ⇐⇒ 0 < Height(Φ τ (v)) < µ C (τ )
This observation allows to define the following conical projection: Definition 4.9 (Shadow). Let v be a vertex in P τ . The shadow Sh τ (v) of v is a point in the (n -k)-dimensional Euclidean space bis τ defined as the intersection of the half-line starting at (0, µ C (τ )) and going through Φ τ (v) with the space bis τ .

Sh τ (P, µ) = def. µ C (τ ) µ C (τ ) -Height(Φ τ (P, µ)) (π bisτ (P ) -o τ )
The shadow of a simplex σ ∈ K τ is a simplex in bis τ whose vertices are the shadows of vertices of σ.

Let Γ reg be the n-chain containing the n-simplices of the regular triangulation of P and |Γ reg | the corresponding simplicial complex. Denote by X(τ )

∈ C n-k-1 (K τ ) the (n -k -1)-chain made of simplices σ ∈ K τ such that τ ∪ σ ∈ |Γ reg |: X(τ ) = def. {σ ∈ K τ , dim(σ) = n -k -1, τ ∪ σ ∈ Γ reg } (19)
In the following, we call polytope a finite intersection of closed half spaces ∩ i H i . The convex cone of a polytope at a point p is the intersection of all such H i whose boundary contain p. Definition 4.10 (Polytope facet visible from the point 0). We say that a facet f of a polytope is visible from the point 0, or visible for short, if the closed half-space H containing the polytope and whose boundary is the supporting plane of f does not contains 0. For each upper half-space H j contributing to the convex cone of the lower convex hull of Φ τ (P) at (0, µ C (τ )), the intersection H j ∩ bis τ is a (n -k)-dimensional half-space in bis τ . The shadow polytope is precisely defined as the intersection of all such half-spaces H j ∩ bis τ . Since each H j is a upper half-space and since by observation 4.8 one has µ C (τ ) > 0, H j does not contains (0, 0), which implies that H j ∩ bis τ does not contain the point o τ in bis τ . It follows that: Observation 4.12. All facets of the shadow polytope are visible from 0. Lemma 4.13 (Proof in Appendix H). Let P be in general position with non positive weights.

1. τ is in the regular triangulation of P if and only if Φ τ (τ ) = (0, µ C (τ )) is an extremal point of the convex hull of Φ τ (P).

2. When τ is in the regular triangulation of P, its link is isomorphic to the link of the vertex Φ τ (τ ) = (0, µ C (τ )) in the simplicial complex corresponding to the boundary of the lower convex hull of Φ τ (P).

3. When τ is in the regular triangulation of P, Sh τ induces a bijection between the simplices in X(τ ) and the set of bounded facets of the boundary of the shadow polytope.

Convex cone as lexicographic minimal chain

This section is self-contained and does not relies on previous constructions. Lemma 5.4 is a key ingredient of the proof of Theorem 3.1 but is also a result of independent interest: visible convex hulls can be defined as minimal lexicographic chains. Lemmas 5.2 and 5.3 (Proof in Appendix I and J) will be instrumental in the proof of Lemma 5.4.

Definition 5.1 (Trace of a chain in a link

). Given a k-simplex τ in a simplical complex K and a n-chain on K for n > k, we call trace of Γ on the link of τ the (n -k -1)-chain Tr τ (Γ) defined in the link of τ by:

Tr τ (Γ)(σ) = def.
Γ(τ ∪ σ) Lemma 5.2 (Proof in Appendix I). Given a k-simplex τ in a simplical complex K and a n-chain Γ on K for n > k, one has:

∂ Tr τ (Γ) = Tr τ (∂Γ) Lemma 5.3 (Proof in Appendix J). Let C ⊂ R n be a polytope and O ∈ R n \ C. Let X ⊂ ∂C be a compact set union of facets of C visible from O. If x ∈ X maximizes the distance to O, then x is in the closure of ∂C \ X.
We need to define another order on simplices together with its induced order on chains, respectively denoted ≤ Sh and Sh . We associate to a (n -1)-simplex σ in R n that does not contain 0 a dimension increasing sequence of faces ∅ = τ -1 (σ) ⊂ τ 0 (σ) ⊂ . . . ⊂ τ n-1 (σ) = σ with dim(τ i ) = i. Under a simple generic condition, it is defined as follows.

τ -1 (σ) = ∅ and τ 0 (σ) is the vertex of σ farthest from 0. More generally, define the distance from a flat (an affine space) F to 0 as d 0 (F ) = inf p∈F d(p, 0). If ζ is a non degenerate i-simplex for i ≥ 0, defines d 0 (ζ) = d 0 (F (ζ)) where F (ζ) is the i-dimensional flat support of ζ. For i ≥ 0, τ i (σ) is the coface of dimension i of τ i-1 (σ) whose supporting i-flat is farthest from 0:

τ i (σ) = def. arg max ζ⊃τ i-1 (σ) dim(ζ)=i d 0 (ζ) (20)
For i = 0, . . . , n-1 we set δ i (σ) = d 0 (τ i (σ)) and the comparison < Sh between two (n-1)-simplices σ 1 and σ 2 is a lexicographic order on the sequences (δ i (σ 1 )) i=0,...n-1 and (δ i (σ 2 )) i=0,...n-1 :

σ 1 < Sh σ 2 ⇐⇒ def. ∃k ≥ 0, δ k (σ 1 ) < δ k (σ 2 ) and ∀j, 0 ≤ j < k, δ j (σ 1 ) = δ j (σ 2 ) ( 21
)
which defines an order relation:

σ 1 ≤ Sh σ 2 ⇐⇒ def. σ 1 = σ 2 or σ 1 < Sh σ 2 (22)
Condition 5.1. Let K be a (n -1)-dimensional simplicial complex. For any pair of simplices

σ 1 , σ 2 ∈ K: dim(σ 1 ) = dim(σ 2 ) = k and d 0 (σ 1 ) = d 0 (σ 2 ) ⇒ σ 1 = σ 2 .
Under Condition 5.1, ≤ Sh is a total order on simplices and, following Definition 2.1, the order ≤ Sh on simplices induces a lexicographic order Sh on k-chains of K. Lemma 5.4. Let P be a set of points in R n such that 0 ∈ R n is not in the convex hull of P . Let K be the complete (n -1)-dimensional simplicial complex over P , i.e. the simplicial complex made of all (n -1)-simplices whose vertices are points in P with all their faces. Assume that K satisfies the generic condition 5.1. Let X be a (n -1)-chain in K whose (n -1)-simplices are on the boundary of the convex hull of P and are all visible from 0 ∈ R n . Then:

P X

X = min Sh {Γ ∈ C n-1 (K), ∂Γ = ∂X} ( 23
)
where when n = 1 the boundary operator in (23) is meant as the boundary operator of reduced homology, i.e the linear operator ∂0 : C n-1 (K) → Z 2 that counts the parity.

The lemma is illustrated for n = 2, 3 on figures 5 and 6.

Proof. We first claim that the lemma holds for n = 1. In this case the fact that 0 is not in the convex hull of P means that the 1-dimensional points in P are either all positive, either all negative. The single simplex in the convex hull boundary visible from 0 is the point in P closest to 0, i.e. the one with the smallest absolute value, which corresponds to the minimum chain with odd parity in the Sh order, which proves the claim.

We assume then the theorem to be true for the dimension n -1 and proceed by induction. This recursion is illustrated on figure 5 for n = 2 and figure 6 for n = 3.

Consider the minimum:

Γ min = min Sh {Γ ∈ C n-1 (K), ∂Γ = ∂X} (24)
We need to prove that Γ min = X.

Let v be the (unique) vertex in the simplices of ∂X which is farthest from 0. Since v is a vertex in at least one simplex in ∂X = ∂Γ min , it must be a vertex in some simplex in Γ min .

Thanks to Lemma 5.3, if a point x is a local maximum in X of the distance to 0 one has x ∈ ∂X. It follows that v is also the vertex in the simplices of X which is farthest from 0.

Since v is the vertex in X farthest from 0 and since by definition Γ min Sh X, we know that Γ min does not contain any vertex farther from the origin than v, therefore v is also the vertex in the simplices of Γ min farthest from 0.

Since ∂Γ min = ∂X, Lemma 5.2 implies that:

∂ Tr v (Γ min ) = Tr v (∂Γ min ) = Tr v (∂X) (25)
In order to define a lexicographic order on chains on the link of v in K, we consider the hyperplane Π containing 0 and orthogonal to the line 0v. We associate to any (n -2)-simplex η ∈ Lk K (v) the (n -2) simplex π vΠ (η) conical projection of η on Π with center v. In other words, if u is a vertex of η: where uv denote the line going through u and v. The map π vΠ is a conical projection on vertices but it extends to a bijection on simplices and an isomorphism on chains that trivially commutes with the boundary operator. definition of the lexicographic order Sh , the comparison of two chains whose farthest vertex is v starts by comparing their restrictions to the star of v. Therefore, since v is the farthest vertex in Γ min , the restriction of Γ min to the star of v must be minimum under the constraint ∂Γ = ∂X. The constraint ∂Γ = ∂X for the restriction of Γ min to the star of v is equivalent to the constraint given by equation (25) or equivalently by:

{π vΠ (u)} = Π ∩ uv
∂π vΠ (Tr v (Γ min )) = π vΠ (Tr v (∂X))
and the minimization on the restriction of the (n -1)-chain Γ min to the star of v can equivalently be expressed as the minimization of the (n -2)-chain γ min = π vΠ (Tr v (Γ min )) under the constraint ∂γ min = π vΠ (Tr v (∂X)), we have then:

γ min = π vΠ (Tr v (Γ min )) = π vΠ Tr v min Sh {Γ ∈ C n-1 (K), ∂Γ = ∂X} (26) = min Sh {γ ∈ C n-2 (π vΠ (Lk K (v))) , ∂γ = π vΠ (Tr v (∂X))} = min Sh {γ ∈ C n-2 (π vΠ (Lk K (v))) , ∂γ = ∂π vΠ (Tr v (X))}
In the third equality of (26) we have used the fact that the orders on (n -1)-simplices in the star of v in K and the order on corresponding (n -2)-simplices in the image by π vΠ of the link of v are compatible.

Indeed, if F is a k-flat in R n going through v, we have (see Figure 7):

d 0 (π vΠ (F )) = d 0 (F ∩ Π) = d 0 (F ) v -0 (v -0) 2 -d 0 (F ) 2 (27)
with the convention d 0 (π vΠ (F )) = +∞ in the non generic case where F ∩ Π = ∅ (the denominator vanishes in this case while since v ∈ F one has

F ∩ Π = ∅ ⇒ d 0 (F ) < v -0 ). As seen on (27) d 0 (F ) → d 0 (π vΠ (F )
) is an increasing function and the orders are therefore consistent along the induction. We claim that the minimization problem in the last member of (26) satisfies the condition of the theorem for n = n -1 which is assumed true by induction.

The recursion is as follows: hyperplane Π corresponds to R n with n = n -1 and:

• n ← n -1

• P ← π vΠ (P \ {v}) • K ← π vΠ (Lk K (v)) • X ← π vΠ (Tr v (X))
Since (n -1)-simplices in X are in convex positions, hyperplanes supporting these simplices, in particular the simplices in the star of v, separate all points of P from 0. It follows that the intersection of theses hyperplanes with the horizontal hyperplane, i.e. the images by π vΠ of the hyperplanes, separates P = π vΠ (P \ {v}) from 0. It follows that the (n -1)-simplices in X = π vΠ (Tr v (X)) are in convex position and are visible from 0. Therefore one can apply our lemma recursively, which gives us, using (26):

γ min = π vΠ (Tr v (Γ min )) = π vΠ (Tr v (X))
It follows that the faces in the star of v corresponding to Tr v (Γ min ) belong to X. Call Y the (n -1)-chain made of these simplices in the star of v. We have both ⊂ X and Y ⊂ Γ min . Since v is the vertex farthest from 0 in both X and in Γ min one has by definition of the lexicographic order:

Γ min = min lex {Γ ∈ C n-1 (K), ∂Γ = ∂X} = Y + min lex {Γ ∈ C n-1 (K), ∂Γ = ∂(X -Y )}
So, by considering the new problem X ← (X -Y ) and iterating as long as X is not empty, we get our final result Γ min = X.

Proof of Theorem 3.1

The next two lemmas establish the between the orders lex and Sh .

In the remaining of this proof we use the lexicographic order Sh on shadows of (n -k -1)chains in K τ , defined at the beginning of Section 5. This order is equivalent to the order lex on corresponding n-chains restricted to the set of simplices with bounding weight µ = B (Lemmas 6.1 and 6.2). This correspondence allows to conclude the proof by applying Lemma 5.4 that says the chain defined by visible faces of a polytope minimises the lexicographic order Sh among chains with same boundary. More formally, thanks to Lemma 6.2 we have:

↓ Γ 1 lex ↓ Γ 2 ⇒ Sh τ (Tr τ (↓ Γ 1 )) Sh Sh τ (Tr τ (↓ Γ 2 ))
It follows that Sh τ (Tr τ (↓ Γ min )) is, among all chains in the complex Sh τ (K τ ), the one that minimises Sh under the constraint (32), or equivalently:

∂ Sh τ (Tr τ (↓ Γ min )) = ∂ Sh τ (Tr τ (↓ Γ reg )) Lemma 5.4 applied with: n ← n -k X ← Sh τ (Tr τ (↓ Γ reg )) = Sh τ (X(τ )) P ← Sh τ (P τ ) K ← Sh τ (K τ ) implies: Sh τ (Tr τ (↓ Γ min )) = Sh τ (Tr τ (↓ Γ reg ))
In other words, Γ min and Γ reg coincide on simplices with bounding weight µ = B , a contradiction with the definition of µ = B .

Appendices

A Invariance by global weight translation

Observation 7.1. Let ψ λ : R n × R → R n × R, be the transformation that shifts the weight by λ:

ψ λ (P, µ) = (P, µ + λ)
Let σ ∈ K P be a simplex, from definitions 3.3 and 3.2 we have:

P C (ψ λ (σ)) = P C (σ) and µ C (ψ λ (σ)) = µ C (σ) -λ P B (ψ λ (σ)) = P B (σ) and µ B (ψ λ (σ)) = µ B (σ) -λ
It follows that a global shift by a constant value λ results in an opposite shift on the weights of generalized circumcenters. It therefore preserves the relative order between simplices weights µ C and µ B :

µ C (σ 1 ) ≤ µ C (σ 2 ) ⇐⇒ µ C (ψ λ (σ 1 )) ≤ µ C (ψ λ (σ 2 ))
and the same relation holds for µ B . Since the order ≤ between simplices (5) defined in section 3.3 relies entirely on comparisons on µ C and µ B , this total order is preserved by a global weight translation.

B Proof of Lemma 3.4

Proof. If σ = {(P 0 , µ 0 ), . . . , (P k , µ k )}, |σ| is the convex hull of {P 0 , . . . , P k }. If P / ∈ |σ| the projection of P on |σ| decreases the weighted distance from P to all the vertices of σ, which shows that (P, µ) cannot realize the arg min in [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF].

C Proof of Lemma 3.5

Figure 1 illustrates the possibilities for Θ(σ) in the case n = 3 and zero weights.

Proof. In (2), denote by τ ⊂ σ the set of vertex (P i , µ i ) ∈ σ for which: D ((P, µ), (P i , µ i )) = 0 This set cannot be empty as, if all inequalities in (2) where strict, a strictly smaller value of µ would still match the inequality, which would contradict with the arg min in [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF]. One has then τ = ∅ and:

∀(P i , µ i ) ∈ τ, D ((P, µ), (P i , µ i )) = 0 (33) and of course:

∀(P i , µ i ) ∈ τ, D ((P, µ), (P i , µ i )) ≤ 0 (34)
No other (P, µ) with a smaller value of µ can satisfies (33) nor (34) as it would again similarly contradict the arg min in [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF]. It follows that:

(P B , µ B )(σ) = (P B , µ B )(τ ) = (P C , µ C )(τ )

D Proof of Proposition 4.3

Proof. Note that Γ reg is unique under the assumed general position. Since, in the regular triangulation, all (n -1)-simplices that are not on the boundary of CH(P) are shared by exactly two n-simplices, while only those in β P have a single n-coface, we have:

∂Γ reg = β P
We claim now that :

∂Γ = β P ⇒ Γ reg (p) ≤ Γ (p) (35)
Indeed, ( 8) and ( 9 From the equivalence between regular triangulations and convex hull on lifted points we know that if σ reg ∈ Γ reg , then for any n-simplex in σ ∈ K :

x ∈ |σ| ∩ |σ reg | ⇒ δ σreg (x) ≤ δ σ (x) (37)
According to Lemma 7.2, in (36), there is an odd number of, and therefore at least one, simplex σ ∈ Γ satisfying |σ| x in the condition on the sum. Therefore (37) gives:

x ∈ |σ reg | ⇒ δ σreg (x) p ≤ σ∈Γ |σ| x δ σ (x) p (38)
And, since, from definition of triangulation, for x ∈ CH(P) \ |K

[x x] ∩ |K n-2 P | = ∅ (40)
where [x x] denote the line segment in R n between x and x. Indeed, we consider moving a point x t from x 0 to some x 1 , picking x 0 , x 1 far away enough to have [x 0 x 1 ] ∩ CH(P) = ∅ and in such a way that [x 0 x 1 ] belongs to none of the affine hyper-planes spanned by (n -1)-simplices in K P ), which occurs generically. Then, the negation of condition (40) with x = x t , occurs only as isolated values of t. We can pick a value t for which is does not occur: set x = x t and the claim is proved. | Next, we navigate a point y(t) = (1 -t)x + tx along segment [x x]. This segment intersect transversally the (n -1)-faces |τ |, for τ ∈ K P . At each intersection point we can keep track of the change in the number of covering n-simplices, where, by covering simplices we name the n-simplices σ ∈ Γ such that x ∈ |σ|.

We know that this number is zero at x since x / ∈ CH(P). Since CH(P) ∩ [x x] is convex there is a single intersection point y(t b ) between [x x] ans the boundary of CH(P). This point y(t b ) ∈ (x , x) hits a (n -1)-simplex τ b ∈ |β P |, face of the convex hull boundary. Since ∂Γ = β P , we know that τ b is shared by an odd number n b of n-simplices in |Γ|. By definition of the convex hull, and since P is in general position, for each n-simplex σ coface of τ b , |σ| is on the inner side of the convex hull supporting half plane. It follows that the number of covering simplices become the odd number n b just after the first crossing.

Then when crossing any other (n -1)-simplex τ i / ∈ |β P |, at some point y(t i ), the condition ∂Γ = β P requires the number n i of n-simplices in |Γ| coface of τ i to be even. When crossing |τ i |, along [x x], point y(t i ) exits k -, and enters k + n-simplices in |Γ|, with k -+ k + = n i . The current number of covering n-simplices value is incremented by k + -k -. Since n i is even and:

k + -k -= k + + k --2k -= n i -2k - k + -k -is

E Proof of Lemma 4.4

Proof. Since the expression (6) does not depend on the origin O, let us choose this origin to be O = P B (σ). With barycentric coordinates based on the vertices P i , i = 0, k of σ, i.e. λ i ≥ 0, i λ i = 1 such that x = i λ i P i , the expression of f σ is:

f σ (x) = i λ i (P i -P B (σ)) 2 -µ i -(x -P B (σ)) 2 one has (P i -P B (σ)) 2 -µ i -µ B (σ) = D ((P B (σ), µ B (σ)), (P i , µ i )) ≤ 0 so that (P i -P B (σ)) 2 -µ i ≤ µ B (σ) and it follows that: ∀x, f σ (x) ≤ µ B (σ) -(x -P B (σ)) 2 (41)
We have from Lemma 3.4 that P B (σ) ∈ |Θ(σ)| ⊂ |σ| so that, in the expression of P B (σ) as a barycenter of vertices of σ, only coefficients λ i corresponding to vertices of Θ(σ) ⊂ σ are non zero:

P B (σ) = (Pi,µi)∈Θ(σ) λ i P i
One has by definition of (P C , µ C ):

(P i , µ i ) ∈ Θ(σ) ⇒ (P i -P C (Θ(σ))) 2 -µ i = µ C (Θ(σ))
and since we know from Lemma 3.5 that:

(P B , µ B )(σ) = (P B , µ B )(Θ(σ)) = (P C , µ C )(Θ(σ))
one gets:

(P i , µ i ) ∈ Θ(σ) ⇒ (P i -P B (σ)) 2 -µ i = µ B (σ)
and:

f σ (P B (σ)) =   (Pi,µi)∈Θ(σ) λ i ((P i -P B (σ)) 2 -µ i )   -0 2 = µ B (σ)
This with (41) ends the proof.

F Proof of Lemma 4.5

Proof. Consider the smallest ratio between two bounding of n-simplices in K P :

ι = inf σ 1 ,σ 2 ∈K [n] P µ B (σ 1 )<µ B (σ 2 ) µ B (σ 2 ) µ B (σ 1 )
From the finiteness of simplices, one has ι > 1 and since for a finite number of n-simplices σ ∈ K

[n] P we have from (11): lim

p→∞ w p (σ) = µ B (σ)
It follows that there is some p 0 large enough such that for any p > p 0 one has:

∀σ 1 , σ 2 ∈ K [n] P , µ B (σ 1 ) < µ B (σ 2 ) ⇒ µ B (σ 2 ) µ B (σ 1 ) > 1 + ι 2
Let N be total number of n-simplices in K P . Taking:

p = max p 0 , log N log 1+ι 2
realizes the statement of the lemma.

G Proof of Lemma 4.6

Proof of Lemma 4.6. Under generic condition, σ is in the regular triangulation of P if and only if there is a weighted point (P C (σ), µ C (σ)) such that:

∀(P i , µ i ) ∈ σ, D ((P C (σ), µ C (σ)), (P i , µ i )) = 0 (42) ∀(P i , µ i ) ∈ P \ σ, D ((P C (σ), µ C (σ)), (P i , µ i )) > 0 (43)
Observe that, since τ ⊂ σ, (42) implies:

∀(P i , µ i ) ∈ τ, D ((P C (σ), 0), (P i , µ i )) = µ C (σ)
This and the definition (12) of bisector show that P C (σ) must be on the bisector of τ :

H Proof of Lemma 4.13

Proof. For 3. consider a simplex σ ∈ X(τ ). Φ τ (σ) is a simplex in the link of Φ τ (τ ) in the lower convex hull of Φ τ (P) and therefore the convex cone CC σ with apex Φ τ (τ ) and going through Φ τ (|σ|) is on the boundary of the convex cone CC P with apex Φ τ (τ ) and going through the convex hull of Φ τ (P).

Therefore Sh τ (σ), intersection of CC σ with bis τ is on the boundary of the shadow polytope, intersection of CC P with bis τ .

Sh τ (σ) is bounded as being the convex hull of the shadow of its vertices. Thanks to observation 4.12, it is a facets of the boundary of the shadow polytope visible from 0. In the reverse direction, a bounded facet of the boundary of the shadow polytope is precisely the shadow of a simplex Φ τ (σ) in the link of Φ τ (τ ) in the lower convex hull of Φ τ (P) with all its vertices in P τ . Therefore σ ∈ X(τ ).

I Proof of Lemma 5.2

Proof. We need to prove that for any simplex σ in the link of τ , one has:

∂ Tr τ (Γ)(σ) = Tr τ (∂Γ)(σ)
We have by definition that a (n -k -2)-simplex σ is in Tr τ (∂Γ)(σ) if and only if τ ∩ σ = ∅ and τ ∪ σ is in ∂Γ. In other words, τ ∪ σ has an odd number of n-cofaces in the chain Γ. This in turn means that σ has an odd number of (n -k)-cofaces in the trace of Γ in the link of τ , i.e. σ ∈ ∂ Tr τ (Γ).

J Proof of Lemma 5.3

Proof. Assume for a contradiction that x is in the relative interior of X, that is there is some ρ > 0 such that B(x, ρ) ∩ X = B(x, ρ) ∩ ∂C. Then all facets containing x are visible from O. if x is not a vertex pf C, it belongs then to the relative interior of a convex face f in ∂C with dim f ≥ 1. Then we have a contradiction since the function y → d(P, y) is convex on f and cannot have an interior local maximum at x. We assume now that x is a vertex of C.

Following for example [START_REF] Federer | Curvature measures[END_REF][START_REF] Tyrrell | Convex analysis[END_REF], denote by Tan x C and Nor x C respectively the Tangent and Normal cone to C at x. In case of a closed polytope they can be expressed as: Since C is a convex polytope, Tan x C is a convex closed cone and one has [START_REF] Tyrrell | Convex analysis[END_REF]:

Tan x C = (Nor x C) ⊥ = def.
{v, ∀u ∈ Nor x C, u, v ≤ 0} (47)

Each facet F i of C containing x is supported by a half-space H i = {y, y -x, n i ≤ 0} and one has:

Nor x C = i λ i n i , ∀i, λ i ≥ 0 (48) ∀i, O -x, n i > 0
This with (48) gives: ∀u ∈ Nor x C, O -x, u > 0 (49) using (47) we get that O -x is in the interior of -Tan x C i.e.:

x -O ∈ (Tan x C)

• ( 50
)
Since x is a vertex of C, Tan x C is the convex hull of Tan x ∂C, therefore (50) implies at there are t 1 , . . . t n+1 ∈ Tan x ∂C and λ 1 , . . . λ n+1 ≥ 0 such that:

x -O = i λ i t i

1

 1 be in general position and let K P be the n-dimensional full complex over P. Denote by β P ∈ C n-1 (K P ) the (n-1)-chain made of simplices belonging to the boundary of CH(P). If Γ min = min lex {Γ ∈ C n (K P ), ∂Γ = β P } the simplicial complex |Γ min | support of Γ min is the regular triangulation of P.

Figure 1 :

 1 Figure 1: Illustration of the definition of Θ(σ) for a tetrahedron σ = abcd in the case of zero weights.

Figure 2 :

 2 Figure 2: Illustration for the definition of bis τ , π τ (top left and right) and Φ τ (bottom).

Lemma

Figure 3 :

 3 Figure 3: Illustration of the shadow polytope (grey area) of Definition 4.11 corresponding to the example of Figure 2.

Definition 4 .

 4 11 (Shadow Polytope). The (possibly empty) intersection of the convex cone of the lower convex hull of Φτ (P) at Φ τ (τ ) = (0, µ C (τ )) with bis τ is called shadow polytope of τ .

Figure 3

 3 Figure3depicts the shadow polytope corresponding to the example of Figure2as the hatched area and the bounded cells of its boundary by the two edges in blue.For each upper half-space H j contributing to the convex cone of the lower convex hull of Φ τ (P) at (0, µ C (τ )), the intersection H j ∩ bis τ is a (n -k)-dimensional half-space in bis τ . The shadow polytope is precisely defined as the intersection of all such half-spaces H j ∩ bis τ . Since each H j is a upper half-space and since by observation 4.8 one has µ C (τ ) > 0, H j does not contains (0, 0), which implies that H j ∩ bis τ does not contain the point o τ in bis τ . It follows that:

Figure 4 :

 4 Figure 4: Illustration of Lemma 5.4 for n = 1.

Figure 5 :

 5 Figure 5: Illustration of the recursion in the proof of Lemma 5.4 for n = 2.

Figure 6 :

 6 Figure 6: Illustration of the recursion in the proof of Lemma 5.4 for n = 3.

Figure 7 :

 7 Figure 7: Illustration for Equation (27).

  even and the number of covering simplices remains odd all along the path [x x].

P

  C (σ) ∈ bis τ(44) We get, for (P i , µ i ) ∈ P:D ((P C (σ), µ C (σ)), (P i , µ i )) = (P C (σ) -P i ) 2 -µ C (σ) -µ i = (π bisτ (P i ) -P C (σ)) 2 + d(P i , bis τ ) 2 -µ C (σ) -µ i = D bisτ (P C (σ), µ C (σ)), (π bisτ (P i ), µ i -d(P i , bis τ ) 2 ) = D bisτ ((P C (σ), µ C (σ)), π τ (P i , µ i ))(45)In the last two lines, the weighted distance is denoted D bisτ instead of D in order to stress that, thanks to (44) , it occurs on weighted point of bis τ × R rather than R n × R. It result that (42) and (43) are equivalent to:∀(P i , µ i ) ∈ σ, D bisτ ((P C (σ), µ C (σ)), π τ (P i , µ i )) = 0 ∀(P i , µ i ) ∈ P \ σ, D bisτ ((P C (σ), µ C (σ)), π τ (P i , µ i )) > 0which precisely means that π τ (σ) is a coface of the vertex π τ (τ ) = {(o τ , -µ C (τ ))} in the regular triangulation of π τ (P).

  Tan x C = ρ>0 {λ(c -x), λ ≥ 0, c ∈ B(x, ρ)} and Nor x C = (Tan x C) ⊥ = def.{u, ∀v ∈ Tan x C, u, v ≤ 0} (46)

  is in the regular triangulation of P, π τ induces a simplicial isomorphism between the link of τ in the regular triangulation of P and the link of vertex π τ (τ ) in the regular triangulation of π τ (P).

4.6 (Proof in Appendix G). Let P be in general position, τ ∈ K P a k-simplex and σ ∈ K P a coface of τ . Then σ is in the regular triangulation of P if and only if π τ (σ) is a coface of the vertex {(o τ , -µ C (τ ))} = π τ (τ ) in the regular triangulation of π τ (P). An immediate consequence of Lemma 4.6 is: Corollary 4.7. The projection π τ preserves the structure of the regular triangulation around τ , more precisely: 1. the simplex τ is in the regular triangulation of P if and only if the vertex π τ

(τ ) = {(o τ , -µ C (τ ))}

is a vertex of the regular triangulation of π τ (P),

2. if τ

  n-1 P | there is exactly one simplex σ reg such that |σ reg | x , (38) can be rewritten as: and since the function is continuous, this implies Γ (p) > Γ reg (p) . Given P = {P 1 , . . . , P N } ⊂ R n , with N ≥ n + 1, in general position, denote by β P ∈ C n-1 (K P ) the (n -1)-chain made of simplices belonging to the boundary of CH(P). Let Γ ∈ C n (K P ) be such that:∂Γ = β P If x ∈ CH(P) \ |K n-1 P |then there is an odd number of n-simplices σ ∈ Γ such that x ∈ |σ|. Proof. We claim that since x ∈ CH(P) \ |K n-1 P | there is x ∈ R n \ CH(P) such that:

	δ σ (x) p ≤	δ σ (x) p	(39)
	σ∈Γreg	σ∈Γ
	|σ| x	|σ| x
	which, together with (36) gives the claim (35).
	Now, if some n-simplex σ ∈ Γ with |σ| x, for some x ∈ CH(P) \ |K n-1 P |, is not Delaunay,
	then
	δ σ (x) p = δ σreg (x) p <	δ σ (x) p
	σ∈Γreg	|σ| x
	|σ| x
	Lemma 7.2.

  1. follows from Corollary 4.7 item 1. together with Proposition 4.2. 2. follows from Corollary 4.7 item 2. together with Proposition 4.2.

Lemma 6.1 (Proof in Appendix K). For σ 1 , σ 2 ∈ K τ , one has:

For a n-chain Γ denote by ↓ ρ Γ the chain obtained by removing from Γ all simplices with bounding weight strictly greater than ρ.

{σ ∈ Γ, µ B (σ) ≤ ρ} (29) Lemma 6.2 (Proof in Appendix L). For two n-chains Γ 1 , Γ 2 ∈ C n (K P ) if Γ 1 = Γ 2 one has:

Proof of Theorem 3.1. We prove Theorem 3.1 in the case of non positive weights which then extends to any weights thanks to Observations 3.6 and 7.1. As in Proposition 4.3, denote by Γ reg the chain that defines the regular triangulation of CH(P). As in Theorem 3.1 denote by β P ∈ C n-1 (K P ) the (n -1)-chain made of simplices belonging to the boundary of CH(P).

According to Proposition 4.3, Γ reg minimizes Γ → Γ (p) among the chains with boundary β P for any p ≥ 1. In particular, Γ reg minimizes Γ → Γ p for the value p of Lemma 4.5.

Proposition 4.3 and Theorem 3.1 consider a minimum with respect to the same boundary condition while their objective differ. In order to prove Theorem 3.1, we have to show that both minimum still agree. By contradiction, we assume now that they differ, which means that, Γ min = Γ reg where Γ min is the minimal chain of Theorem 3.1. Consider µ = B to be the largest bounding weight for which some simplex in Γ min and Γ reg differ:

There must be at least one simplex with bounding weight µ = B in Γ reg as otherwise, by definition of µ = B , there would be a simplex with radius µ = B in Γ min and this would give Γ reg lex Γ min with Γ reg = Γ min and since ∂Γ reg = ∂Γ min = β P this contradicts the definition of Γ min .

Similarly, it follows from Lemma 4.5 that if there was no simplex with bounding weight µ = B in Γ min , one would have Γ min p < Γ reg p and ∂Γ reg = ∂Γ min : a contradiction with the minimality of Γ reg for norm • p (Proposition 4.3). We have shown that if they differ, both Γ reg and Γ min must contain at least one simplex with bounding weight µ = B . We know from Lemma 3.5 and the generic condition 3.2 that the set of simplices with bounding weight µ = B are all cofaces of some unique dimension minimal simplex

is the unique simplex in K P whose bounding weight is µ = B . But then Γ reg and Γ min coincide on simplices with bounding weight µ = B , a contradiction with the definition of

B one has:

where

is defined in (29). In order to spare our eyes, we allow ourselves to replace for the rest of the section ↓ µ = B by ↓ and τ µ = B by τ . It follows from (31) that:

↓ Γ reg and ↓ Γ min have therefore the same boundary and by Lemma 5.2, their trace also have the same boundary:

Observe that Tr τ (↓ Γ reg ) coincide with the definition of X(τ ) in (19). We know then from Lemma 4.13 that the shadows of simplices in X(τ ) = Tr τ (↓ Γ reg ), that is Sh τ (X(τ )) = Sh τ (Tr τ (↓ Γ reg )) is a chain in Sh τ (K τ ) made of the faces of the convex hull of Sh τ (P τ ) visible from the origin 0. which gives:

Since all the λ i are not negative, there must be at least one i for which:

This precisely means that y → d(P, y) is increasing in the direction t i in a neighborhood of x in X, a contradiction since x is assumed to be a local maximum of y → d(P, y) in X.

K Proof of Lemma 6.1

Proof. Using the definition 3.3, one has:

Looking at Definition 3.3 in the light of (45) in the proof of Lemma 4.6 in Appendix G, we get that:

Since both terms of (28) are invariant by a global translation, we can assume without loss of generality and in order to make the computations simpler that o τ = 0.

In this case, as seen in ( 14) the coordinates of π τ (τ ) are (0, -µ C (π τ (τ )) = (0, -µ C (τ )) by ( 52). So that D ((P, µ), π τ (τ )) = 0 gives us:

It follows that among the weighted points (P, µ) that satisfy D ((P, µ), π τ (τ )) = 0, minimizing µ is equivalent to minimizing P 2 and one can reformulate the characterization (51) of (P C , µ C ) (π τ (τ ∪ σ)) = (P C , µ C ) (π τ (τ ) ∪ π τ (σ)) as:

(P,µ)∈R n ×R ∀(P i ,µ i )∈πτ (τ )∪πτ (σ),D((P,µ),(P i ,µ i ))=0 P 2 (54)

Observe that: D ((P, µ),

So that the definition of (P C , µ C ) (π τ (τ ) ∪ π τ (σ)) given in (54) can be equivalently formulated as Π (P C ,µ C ) being the hyperplane in R n-k × R that minimizes P 2 C among all hyperplanes containing both lift(π τ (τ )) and all the points in lift(π τ (σ)).

But, as seen on (55), 2 P C is the slope of the hyperplane Π (P C ,µ C ) so that Π (P,µ) is the hyperplane with minimal slope going through Φ τ (τ ) = lift(π τ (τ )) and all the points in Φ τ (σ) = lift(π τ (σ)). This slope 2 P C is also the slope of the unique (dim(σ) + 1)-dimensional affine space F going through Φ τ (τ ) = lift(π τ (τ )) = (0, µ C (τ )) and all the points in Φ τ (σ) = lift(π τ (σ)). Since F ∩ R n-k × {0} is the affine space supporting Sh τ (σ), one has:

so that, using (52) for the second equality:

It follows that the map:

is decreasing.

L Proof of Lemma 6.2

Proof. If we denote by → ρ Γ the set of simplex in Γ with bounding weight equal to ρ:

We claim that:

It remains to show that the order lex restricted to simplices τ ∪ σ with µ B (τ ∪ σ) = ρ corresponds to the order Sh on the shadow of σ.

By definition of lex , since in (5) one has always µ 0 (τ ∪ σ 1 ) = µ 0 (τ ∪ σ 2 ) = ρ, it goes like this:

Observe that this expression is similar to (21). For a 0-simplex {v} ∈ K P , the circumweight µ C (τ ∪ {v}) is, according to Lemma 6.1, a decreasing function of the distance d 0 (Sh τ (η)) between its shadow and the origin. It follows that for a (n -k -1)-simplex σ ∈ K P , the vertex v for which the circumweight µ C (τ ∪ {v}) is minimal has its shadow Sh τ (v) maximizing the distance to the origin. This minimal circumweight is µ 1 (τ ∪ {v}) while this maximal distance is δ 0 (σ).

More generally, looking at (20) and (4), Lemma 6.1 allows to check that the simplex Θ k (σ) of (4) in the star of τ in K P corresponds to the simplex τ k-1 (σ) in (20) in the link of τ :

So that for σ 1 , σ 2 ∈ K P referring to (21) and (57):

It follows that, for Γ 1 , Γ 2 ∈ C n (K P ) and τ = τ ρ :

which, with claim (56), ends the proof.


Katia Ait-Ameur
email: aitameur.katia@gmail.com

Yvon Maday
email: maday@ann.jussieu.fr

Marc Tajchman
email: marc.tajchman@cea.fr

Time-parallel algorithm for two phase flows simulation

published or not. The documents may come

Introduction

In the nuclear energy domain, computations of complex two-phase flows are required for the design and safety studies of nuclear reactors. System codes are dedicated to the thermalhydraulic analysis of nuclear reactors at the system scale by simulating the whole reactor. We are here interested in the Cathare code developed by the CEA. Like all system codes, Cathare essentially simulates assemblies of onedimensional elements (pipes) and 3D elements (vessels). The discretization level of each element is kept intentionnaly at a coarse level to be able to handle whole systems simulations. Typical meshes used for the simulations are about 10 2 to 10 3 cells with a 3D element. Typical cases involve up to a million of numerical time iterations, computing the approximate solution during long physical simulation times. A space domain decomposition method has already been implemented and to improve the response time, we will consider a strategy of time domain decomposition, based on the parareal method ( [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF]). The paper is organized as follows: after the presentation of the six equation two-fluid model and the Cathare numerical scheme in section 2, the main aspects of the parareal methods will be recalled in section 3. The numerical convergence observed in our example is shown in section 4 followed by the performances we obtain by applying the multistep variant of the parareal algorithm.

Model

At the system scale, the finest details of the flow (description of the liquid-vapor interfaces for example) are not absolutely necessary to obtain a satisfactory macroscale description of the dynamics. For this reason macroscopic models have been developped that focus on the evolution of averaged quantities (see [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF] and [START_REF] Drew | Theory of multicomponent fluids[END_REF]). There are many different averaged models depending on the simplifying assumptions. The model used in Cathare is the 6 equation two-fluid model that considers a set of balance laws (mass, momentum and energy) for each phase, liquid and vapor. It assumes independent velocities and a pressure equality. The unknowns are the volume fraction α k ∈ [0, 1], the pressure p ≥ 0, the velocity u k and the enthalpy H k of each phase. The subscript k stands for l if it is the liquid phase and g for the gas phase. For the sake of simplicity, we write the terms of the model involved in our test case, studied in section 4.

           ∂ t (α k ρ k ) + ∂ x (α k ρ k u k ) = 0 α k ρ k ∂ t u k + α k ρ k u k ∂ x u k + α k ∂ x p = α k ρ k g + F int k ∂ t α k ρ k H k + u 2 k 2 + ∂ x α k ρ k u k H k + u 2 k 2 = α k ∂ t p + α k ρ k u k g (1)
with α g + α l = 1 and the two equations of state :

ρ k = ρ k (p, H k ). The interfacial forces F int k
are of 2 types. The first ensures hyperbolicity of the system (see [START_REF] Ndjinga | Influence of interfacial pressure on the hyperbolicity of the two-fluid model[END_REF] for the well-posedness of the 6 equation model). The second is the interfacial friction term which will be important in the sequel for our test case. In this configuration, the phases are separated which means that one of the two phases vanishes in some parts of the domain. It is numerically challenging to compute the velocity of the ghost phase (see [START_REF] Ndjinga | A 2x2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF]). For this reason, the Cathare scheme forces the two velocities to be equal with this damping term.

Numerical method

The Cathare scheme is based on a finite volume method on a staggered grid (MAC scheme) and on a two step time scheme. In a staggered scheme the i-th component of the velocity is located at the center of the edge orthogonal to the i-th unit vector. Pressures, void fractions and enthalpies are cell-centered. Given a time discretization T 0 , T 1 , T 2 , • • • of the full time interval [0, T ), we use the following notations: (α k ρ k ) n is an approximation of (α k ρ k ) at time T n . Here, we write the time discretization of the Cathare scheme:

                       (α k ρ k ) n+1 -(α k ρ k ) n ∆t + ∂ x (α k ρ k u k ) n+1 = 0 (α k ρ k ) n+1 u n+1 k -u n k ∆t + (α k ρ k u k ) n+1 ∂ x u n+1 k + α n+1 k ∂ x p n+1 = (α k ρ k ) n+1 g + F n,n+1 k 1 ∆t (α k ρ k ) n+1 H k + u 2 k 2 n,n+1 -(α k ρ k ) n H k + u 2 k 2 n-1,n +∂ x α k ρ k u k H k + u 2 k 2 n+1 = α n+1 k p n+1 -p n ∆t + (α k ρ k u k ) n+1 g (2)
Where the notation F n,n+1 k means that the discretization of F k is a function of the approximate solution at times T n and T n+1 . After discretization, the non linear system is solved by a Newton method. Here, we highlight some characteristics of the Cathare scheme, some advantages and limitations:

• The scheme must be accurate enough at the incompressible limit to be able to capture the correct streamlines and pressure fields. This characteristic was mainly studied in the monophasic case:

-Staggered schemes enjoy good precision at the incompressible limit -However, Riemann solvers have poor precision in the incompressible limit. Corrections are proposed in [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF] to overcome this issue

• The vanishing phase is a numerical challenge of two phase flows simulation. It is important to capture well the volume fraction since it governs the composition of the mixture and two-phase/single phase transition . An important issue is to guarantee the positivity of the volume fraction. Many schemes were designed to ensure this property (like [START_REF] Ndjinga | A 2x2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF] for two incompressible phases). Cathare uses a high interfacial friction to deal numerically with these transitions.

The Parareal algorithm

Several approaches have been proposed over the years to decompose the time direction when solving a partial differential equation (see [START_REF] Gander | 50 years of time parallel time integration[END_REF] for an overview). Of these, the parareal in time algorithm, which performances we explore in this work, has received an increasing amount of attention in the last twenty years with many applications (see [START_REF] Baffico | Parallel-in-time molecular dynamics simulations[END_REF], [START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm[END_REF] among many others). In the sequel, we recall the classical parareal algorithm as initially proposed in [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF], [START_REF] Baffico | Parallel-in-time molecular dynamics simulations[END_REF], [START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF] and the principle of the multistep variant we will apply in section 4.

Original parareal algorithm

After the discretization in space of a PDE with N the number of degrees of freedom:

∂ u ∂t + A(t, u) = 0, t ∈ [0, T ], u(t = 0) = u 0 (3) A : R × R N → R N , u ∈ R N (4)
We recall here the classical parareal algorithm as initially proposed in [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF], [START_REF] Baffico | Parallel-in-time molecular dynamics simulations[END_REF], [START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF].

Let G and F be two propagators such that, for any given t ∈ [0, T ], s ∈ [0, Tt] and any function w in a Banach space, G(t, s, w) (respectively F(t, s, w)) takes w as an initial value at time t and propagates it at time t + s. The full time interval is divided into N c sub-intervals [T n , T n+1 ] of size ∆ T that will each be assigned to a processor. The algorithm is defined using two propagation operators:

• G(T n , ∆ T, u n ) computes a coarse approximation of u(T n+1 ) with initial condition u(T n ) u n (low computational cost) • F(T n , ∆ T, u n ) computes a more accurate approximation of u(T n+1 ) with initial condition u(T n ) u n (high computational cost)

Starting from a coarse approximation u n 0 at times T 0 , T 1 , • • • , T N c , obtained using G, the parareal algorithm performs for k = 0, 1, • • • the following iteration:

u n+1 k+1 = G(T n , ∆ T, u n k+1 ) + F(T n , ∆ T, u n k ) -G(T n , ∆ T, u n k )
In the parareal algorithm, the value u(T n ) is approximated by u n k at each iteration k with an accuracy that tends rapidly to the one achieved by the fine solver, when k increases. The coarse approximation G can be chosen much less expensive than the fine solver F by the use of a scheme with a much larger time step (even δ T = ∆ T ) δ T δt (time step of the fine solver) or by using a reduced model. All the fine propagations are made in parallel over the time windows and the coarse propagations are computed in a sequential way but have a low computational cost. The main convergence properties were studied in [START_REF] Gander | Analysis of the Parareal Time-Parallel Time-Integration Method[END_REF] and stability analysis was made in [START_REF] Staff | Stability of the Parareal algorithm[END_REF], [START_REF] Bal | Parallelization in time of (stochastic) ordinary differential equations[END_REF]. We refer to [START_REF] Maday | An adaptive parareal algorithm[END_REF] about the parallel efficiency of parareal and a recent work offering a new formulation of the algorithm to improve the parallel efficiency of the original one.

Multistep variant of parareal

This variant of the parareal algorithm was proposed in [START_REF] Ait-Ameur | Multi-step variant of the parareal algorithm[END_REF]. In the sequel, we will consider that the fine solver is based on a two step time scheme like the Cathare time scheme. Hence we will use the following notation for the fine solver that takes two initial values: F(t, s, x, y), for t ∈ [0, T ], s ∈ [0, Tt[ and x, y in a Banach space.

Example

If one solves (3) with a multistep time scheme as fine propagator F like the order 2 BDF method:

3 2 u j+1 -2u j + 1 2 u j-1 = -δtA(u j+1 ,t j+1 ), j = 1, • • • , N f ,t j+1 -t j = δt
Here the fine solver reads: u j+1 = F(t j , δt, u j-1 , u j ). Now, we apply the parareal algorithm with a coarse grid:

T 0 , • • • , T N c where: T n+1 -T n = ∆ T = N f δt.
Then we can write:

u(T n + jδt) u n, j , j = 1, • • • , N f , n = 1, • • • , N c .
In order to perform the fine propagation, in a given time window [T n , T n+1 ], we only need the local initial condition u n k and a consistent approximation of u(T nδt).

In [2], the authors propose a consistent approximation in the context of the simulation of molecular dynamics. The proposed method was linked to the nature of the model and the symplectic character of their algorithm is shown, which is an important property to verify for molecular dynamics. In the context of our application to the thermalhydraulic code Cathare, we want to derive a multistep variant of parareal that will not be intrusive in the software. We seek a consistent approximation of u(T nδt). The only fine trajectory at our dis-

posal is F(T n-1 , ∆ T, u n-2,N f -1 k , u n-1 k ).
Its final value at T n is:

F(T n-1 , ∆ T, u n-2,N f -1 k , u n-1 k )(T n
) from which we compute u n k+1 by the parareal correction. Hence, we translate the solution:

F(T n-1 , ∆ T -δt, u n-2,N f -1 k , u n-1 k )(T n -δt) by the same correction: u n k+1 -F(T n-1 , ∆ T, u n-2,N f -1 k , u n-1 k
) and obtain the so called consistent approximation u n-1,N f -1 k+1 to initialize the fine propagation in [T n , T n+1 ]. We now detail our algorithm:

                 u n+1 0 = G(T n , ∆ T, u n 0 ), 0 ≤ n ≤ N -1 u n+1 k+1 = G(T n , ∆ T, u n k+1 ) + F(T n , ∆ T, u n-1,N f -1 k , u n k ) -G(T n , ∆ T, u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n,N f -1 k+1 = F(T n , ∆ T -δt, u n-1,N f -1 k , u n k ) + u n+1 k+1 -F(T n , ∆ T, u n-1,N f -1 k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 (5)
Remark: In order to perform the fine propagation, in a given time window [T n , T n+1 ], at the first parareal iteration we need to choose a different consistent approximation of u(T nδt), since we have not used the fine solver yet. To treat this, we could make one iteration with a Backward Euler method or one iteration with a second order Runge Kutta method. In the context of the application to the Cathare code, we choose a non intrusive initialization by imposing u n-1,N f -1 0 = u n 0 .

Test case

Here we apply the multistep parareal algorithm to the resolution of an oscillating manometer. This test case is proposed in [11] for system codes to test the ability of each numerical scheme to preserve system mass and to retain the gas-liquid interface. In this test case, the phases are separated and the interfacial friction term will be important in this configuration.

Note that here we have used the same physical model and the same mesh (110 cells) for both the coarse and the fine solvers: the only difference is the size of the time steps, δt for F and ∆ T for G. All calculations have been evaluated with a stopping criteria where the tolerance is fixed to the precision of the numerical scheme, ε = 5 • 10 -2 . With this threshold, parareal convergence is achieved after 2 or 3 iterations.

In the following subsections, after giving a numerical proof of the convergence of the parareal algorithm in our test case, some results about measured speed-up will be presented.

About the convergence

Figure 1 illustrates that the multistep parareal algorithm effectively converges when applied to the problem of the oscillating manometer. For a given time step T n and parareal iteration k, the relative error in L 2 norm between the parareal solution and the sequential fine solver decreases beyond our given convergence threshold ε. In the figure, the test case has been solved with the multistep parareal algorithm when δt = 10 -5 and ∆ T = 10δt.

These results are obtained on 16 time windows and we will use this configuration in the sequel to study the performances.

Speed-up performances

In the following strong scaling tests, the same setting is used for the multistep parareal algorithm. The test case has been solved on an increasing number N proc of processes N proc = 5, 10, 15, • • • , 70. In figure 2, with 25 processes, we obtain a speed up of 3.4 and of 3.7 with 50 processes. Here, we observe two global trends:

• For N proc = {5, 10, 15, 20, 25, 40, 50}, the speed up first monotonically increases until reaching 25 processes and then increase again with 40 and 50 processes. This is due to the number of parareal iterations that is equal to 2 in this case • For N proc = {30, 35, 45, 55, 60, 65, 70}, the speed up is drastically reduced because the parareal algorithm converges in 3 iterations in this case

In the sequel, we highlight the well-known dependance of the computational cost of the parareal algorithm on the number of iterations. Let T f ine be the cpu time to run the fine solver in a sequential way on the whole time interval [0, T ). Since the coarse time step is ten times greater than the fine time step we suppose that the cpu time of the coarse solver T coarse = T f ine 10 . This ratio between coarse and fine solvers should be as high as possible to minimize the computational cost of the coarse solver which is launched in a sequential way. This aspect will be studied in a forthcoming work to obtain better speed-up performances by coarsening more the solver G. When the algorithm converges in N it iterations, the coarse solver is launched N it times and the fine solver N it -1 times in parallel over the number of processes N proc . Hence, we can write the cpu time in parallel T para in terms of T f ine :

T para = (N it -1) T f ine N proc + N it T coarse + τ = N it -1 N proc + N it 10 T f ine + τ
where τ contains the time of communication between processes and the cpu time for the computation of the parareal corrections and of the error. Now, we can deduce an upper bound of the speed up S when the parareal algorithm converges in 2 or 3 iterations by neglecting τ:

S = T f ine T para ≤ 1 N it -1 N proc + N it 10
Example: On 25 processes, the algorithm converges in 2 iterations: S ≈ 4 when the measured speed up is 3.4. On 35 processes, the algorithm converges in 3 iterations: S ≈ 2.8 when the measured speed up is 2.3.

Conclusion

The results of this study show that the parareal algorithm can effectively speedup two-phase flows simulations. Here we have tested the multistep variant of the parareal algorithm on a test case that is representative of the numerical challenges for two phase flows. These results can certainly be improved by using a new version of parareal (see [START_REF] Maday | An adaptive parareal algorithm[END_REF]) that improves the parallel efficiency of the original parareal. This new algorithm proposes to solve the subproblems at increasing accuracy accross the parareal iterations and should allow us to increase the speed up performances obtained on our test case.

Fig. 1

 1 Fig. 1 Convergence of the multistep parareal algorithm when δt = 10 -5 and ∆ T = 10δt

Fig. 2

 2 Fig. 2 Strong scaling results with the multistep variant of the parareal algorithm

Acknowledgements This research is partially by ANR project CINE-PARA (ANR-15-CE23-0019). We acknowledge the travel support provided by SMAI to attend the Spanish-French School Jacques-Louis Lions.