Bruno Després

Machine Learning design of Volume of Fluid schemes for compressible flows

Keywords: VOF, CFD, ML

Our aim is to establish the feasibility of Machine-Learning-designed Volume of Fluid algorithms for compressible flows. We detail the incremental steps of the construction of a new family of Volume of Fluid-Machine Learning (VOF-ML) schemes adapted to bi-material compressible Euler calculations on Cartesian grids. An additivity principle is formulated for the Machine Learning datasets. We explain a key feature of this approach which is how to adapt the compressible solver to the preservation of natural symmetries. The VOF-ML schemes show good accuracy for advection of a variety of interfaces, including regular interfaces (straight lines and arcs of circle), Lipschitz interfaces (corners) and non Lipschitz triple point (the Trifolium test problem). Basic comparisons with a SLIC/Downwind scheme are presented together with elementary bi-material calculations with shocks.

Introduction

Machine Learning (ML) [START_REF] Goodfellow | Deep Learning[END_REF] for the construction of numerical fluxes adapted to Finite Volume (FV) discretizations is becoming a research subject of its own, see recent contributions for the discretization of hyperbolic equations by Hesthaven et al [START_REF] Ray | An artificial neural network as a troubled-cell indicator[END_REF][START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF]. For viscous incompressible flows, like bubble flows where the curvature of the interface controls the dynamics, it seems that ML techniques are established techniques now [37]: we refer to Zaleski et al. [START_REF] Qi | Computing curvature for volume of fluid methods using machine learning[END_REF][START_REF] Aniszewski | PArallel, Robust, Interface Simulator[END_REF] for the reconstruction of the curvature of interfaces and to [START_REF] Fuster | An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF] for an extension to compressible effects. The exact curvature is function of the second derivative of the function that defines the interface, and a comprehensive review centered on incompressible flows with surface tension is [START_REF] Gibou | Sharp interface approaches and deep learning techniques for multiphase flows[END_REF]. More general references can be found in [START_REF] Raissi | Hidden physics models: Machine learning of nonlinear partial differential equations[END_REF]. On the contrary, for compressible non viscous flows, the interfaces are more related to contact discontinuities and material discontinuities: the dynamics of interfaces is more similar to the one of passive scalars; so it is needed to address interfaces with low regularity. In this context, we refer to [START_REF] Debar | Fundamentals of the KRAKEN code[END_REF][START_REF] Youngs | Time-dependent multi-material flow with large fluid distortion, Numerical methods for fluid dynamics[END_REF] for historical references on VOF methods (KRAKEN code, YOUNGS method) for compressible flows and [START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF] for a more comprehensive presentation of the topic (algorithms LVIRA, ELVIRA and GRAD). Another name in the field is the PLIC method (Piecewise linear interface calculation) [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF], the Youngs algorithm being an example. Recent developments on the MOF method which adresses high order extensions with different ideas are described in [START_REF] Ahn | Multi-material interface reconstruction on generalized polyhedral meshes[END_REF]. The objective of this work is to explain that an avenue completely different from PLIC, YOUNGS, LVIRA, ELVIRA, GRAD or MOF can be walked through for compressible solvers, which is the development of a ML strategy. ML techniques have their own philosophy and techniques [START_REF] Goodfellow | Deep Learning[END_REF][START_REF] Chollet | Deep learning with Python[END_REF] since they are not based on analytical formulas, but on the construction of large datasets and on algorithmic learning of essential features encoded in these datasets. This class of methods has proved to be efficient for image identification and image comparison, so we believe it makes sense to consider that interface reconstruction and VOF procedures can also be addressed within the ML paradigm. Also it is reasonable to state that the performances of standard VOF methods for triple point problems often encountered for multi-material problems (3 phases and more, 2 phases near a wall, 2 materials sliding on a 3rd one, . . . ) suffer restrictions, with the notable exception of the new family of MOF methods [START_REF] Ahn | Multi-material interface reconstruction on generalized polyhedral meshes[END_REF]. This class of problems, which is our longterm objective, is another reason why the desire to establish the feasibility of ML techniques in the context of VOF. In this work, we do not make an extensive comparison of the quality of our new solvers with respect to the literature, even if some simple comparisons will be proposed with the SLIC method of Noh and Woodward [START_REF] Noh | SLIC (Simple Line Interface Calculation[END_REF] implemented as the anti-diffusive scheme of Lagoutière [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF] (it will be referred to as the SLIC/Downwind scheme in the core of this article). We concentrate hereafter on the feasibility of this new family of methods.

The kind of ML algorithms that we use corresponds to supervised training. We refer to [START_REF] Goodfellow | Deep Learning[END_REF][START_REF] Chollet | Deep learning with Python[END_REF] for a state of art description of this approach. In more numerical words, it corresponds to the construction of an approximate/interpolated function defined through given points, the ensemble of these points is the dataset. The two main steps are the construction of one (or more) dedicated dataset(s), then the construction of the approximate function (this is called the training session). Both steps are crucial for the quality of the final results. This methodology is very classical, however the recent progresses of publicly available dedicated softwares make this task easy and powerful in terms of the size of the dataset which can be large and in terms of the quality of the training which is based in our case on dense Convolution Neural Network with many layers. A recent observation is that Neural Networks with high number of layers are quite efficient (we reach the same conclusion) and that the ReLU activation function provides enough accuracy (this is confirmed by the recent theoretical works [START_REF] Daubechies | Nonlinear Approximation and (Deep) ReLU Networks[END_REF]). It explains why we do not base our training session on any high order activation functions within a 2-layer structure (contrary to [START_REF] Qi | Computing curvature for volume of fluid methods using machine learning[END_REF]), but as already mentioned on the ReLU function within a higher number of layers (up to 5 in our case).

Another key ingredient in our approach is the use of the Lagrange+remap strategy [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF] on a Cartesian grid for the compressible flow solver (our implementation is based on [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF][START_REF] Després | Numerical methods for Eulerian and Lagrangian conservation laws[END_REF]). The main virtue of Lagrange+remap schemes is the natural decoupling between a Lagrange step where the acoustic part is treated and the interfaces are considered as fixed, and a remap (or projection) step where only the transport part of the equations is treated and the interfaces move. Lagrange+remap schemes are usually deployed as Finite Volume schemes, so there are naturally conservative in masses, total momentum and total energy: this is a valuable property for an accurate calculation of shocks. Also it has been proved in [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF] that the stability of the first order scheme (in terms of the preservation of the bounds for the volume fractions and of the entropy inequality) is independent of the flux for the volume fractions in the remap stage. The Cartesian grid provides a simplification of the data structure which is amenable to reduce the computational burden for the description of the local geometry of the mesh. In practice the quality of the numerical treatment of the interfaces depends of the interface reconstruction technique in the remap stage, ultimately only the flux of mass fractions or volume fractions matters. These features are very specific to Lagrange+remap compressible schemes.

This brief tour of ML features and Lagrange+remap features motivates the development of a ML flux function which aims at an accurate transport/remap of a reconstructed interface. Following Hesthaven's approach [START_REF] Ray | An artificial neural network as a troubled-cell indicator[END_REF][START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF], we focus on the numerical construction of a flux function where the inputs contain the local volume fractions and output is the flux. That is the reconstruction of the interface features (in terms of curvature, angle at the corners, . . . ) is not the main goal, only the remap Finite Volume flux. One advantage of this approach is that it is versatile with respect to the type of different interfaces: we will consider straight lines, arcs of circle, corners with different angles (right angles, acute angles, . . . ) and even a non Lipschitz profile to emulate a simple triple point geometry. The new VOF-ML schemes are restricted to the dimension d " 2 on uniform Cartesian grids, however the proposed methodology can be used a priori for the reconstruction of internal boundaries in many fields, like in references [START_REF] Basting | Extended ALE Method for fluid-structure interaction problems with large structural displacements[END_REF][START_REF] Boilevin-Kayl | A loosely coupled scheme for fictitious domain approximations of fluidstructure interaction problems with immersed thin-walled structures[END_REF] and references therein.

In summary, the main steps of the development of our VOF-ML schemes are firstly the design of good datasets and accurate training, and secondly the incorporation of new fluxes in Lagrange+remap schemes. We will adopt an incremental presentation of our methods. The contributions of this work can be summarized as follows. ' We construct new VOF-ML schemes adapted to bi-material compressible Euler calculations and describe the calculation chain which is made of (1) construction of good datasets, (2) training session and (3) modification of a Lagrange+remap solver. The accuracy of the critical steps is controlled uniformly over the chain (here around 1%). In the training, we used up to 5 layers of neurons (4 dense hidden layers of neurons). ' We show that VOF-ML has a good ability to recover regular interfaces, details of Lipschitz interfaces, such as the corners of the Zalezak test problem [START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF], and even non Lipschitz interfaces (a new Trifolium test problem is proposed). Since the methodology is quite general in terms of the features of the interface, it is an improvement with respect to curvature reconstruction only [START_REF] Qi | Computing curvature for volume of fluid methods using machine learning[END_REF] and to straight line interface reconstruction only [START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF]. ' The implementation proposed in this work preserves the natural symmetries on a Cartesian mesh. ' The cost of the new VOF-ML schemes scales as

T total " T solver `CML interface N 1D (1)
where T solver is the cost of the Finite Volume solver, C ML interface is a constant which depends of the geometry of the interface and of the new VOF-ML schemes and N 1D is the number of cells in one direction. This feature was expected and it is not an original one with respect to the literature. A proof is provided in 2D at the end of this work, the same scaling holds in 3D. Numerical measurements show that C ML interface can be high with respect to the unit cost of the FV solver. However, asymptotically for large N 1D , the cost of the new VOF-ML schemes is negligible. Our contribution is the numerical observation that, for Lagrange+remap calculations, C ML interface is bounded uniformly with respect to the mesh size.

The organization is as follows. Section 2 explains the algorithmic structure. We introduce the key ideas of our approach for straight line interfaces in Section 3. In next Section 4, we generalize the material to other types of interfaces and construct a dimensionless flux. Section 5 is dedicated to technical details about the inference in a C++ code and Section 6 presents an implementation of dimensionless flux which respects natural symmetries and the maximum principle for the mass and volume fractions. The bi-material compressible model which serves in the numerical section is given in Section 7. The numerical results are presented and discussed in Section 8. We end with a conclusion and perspectives. More mathematical comments are in the appendix section, with the description of a new Trifolium test problem and a test problem with vorticity beyond pure solid body rotation.

Notation. Depending on the context, the volume fraction (which is between 0 and 1) will be denoted as α or a 1 . For similar reasons, the interpretation of the indices will be provided by the context. This abuse of notation has the advantage of never using the heavy notation pa 1 q i j which is the value of the volume fraction a 1 in cell pi, jq.

Description of the calculation chain

As stated in the introduction, this work is based on an incremental calculation chain where one step is performed after the other. Since this incremental algorithmic structure has a great influence on the organization of the calculations, we describe below some details and the software stack. A software may change without changing the chain structure. ' First step is the description of the geometry and construction of datasets. These datasets contain vectors of inputs of size less or equal to 26 in our case in 2D. In some of our tests, the number of vectors for the training dataset runs up to 1.4 10 6 and up to 3.5 10 5 for the validation dataset. The numerical values which determine the vectors are obtained by computing simple integrals which correspond to areas. We use Python code to perform these calculations. ' Second step is the training session which is performed with the Keras-Tensorflow suite in this work. It consists in the construction with a stochastic gradient algorithm of a function which interpolates the training dataset at best. The structure of this function depends on the number of layers and neurons per layer. The validation dataset is used to assess that the level of overfitting or underfitting is under control (which is the case in our calculations). ' Third step consists in passing the function to a C++ code for the inference (which means calling the function). Various possibilities exist so far. One can built an API but it has been considered as too heavy so far. Fortunately simple softwares exist, for example frugally-deep [START_REF] Dobias | Frugally-deep github page[END_REF], kerasify [21] and keras2cpp [START_REF]Keras2cpp[END_REF], which all three can be used to call the function from C++. Preliminary tests show that frugally-deep is a possibility, but the unitary cost of calling the function is too high (« 0.03 s). Since it is needed to call the function in all cells near the interface and at all time steps, the CPU cost in a C++ code has been considered too heavy with frugally-deep. The software keras2cpp is used here. ' Fourth step. It concerns the implementation of the new VOF-ML scheme in our finite volume Lagrange+remap solver. In particular, attention is paid to the preservation of natural symmetries since it is not immediate with a VOF-ML scheme (here with accuracy around 1%).

Straight lines: a case of study

Straight-line interfaces serve as a case of study to present the basic features of the construction of the datasets and of the training session. It is generalized to more general types of interfaces in the next sections.

Parametrization of the geometry

For our CFD calculations, we use a Two Dimensional (2D) Cartesian mesh with a mesh size ∆x ą 0. As required in Learning textbooks [START_REF] Goodfellow | Deep Learning[END_REF][START_REF] Chollet | Deep learning with Python[END_REF], it is better to normalize the data. Normalization is also a standard technique in CFD, so we believe it is a good idea to systematically use normalized data in our ML procedure. This is why we consider a normalized Cartesian mesh

C i j " " x " px 1 , x 2 q P R 2 | ´1 2 `i ă x 1 ă 1 2 `i and ´1 2 `j ă x 2 ă 1 2 `j* .
By construction the area of all square cells is 1. The numeration is such that the central point is the center of mass of the central square C 00 . To have notations compatible with the ones needed for the description of ML algorithms, the description of the local geometry needed for VOF methods starts from the data of two functions D : R par ÝÑ R in and E : R par ÝÑ R out , where ' par P N ˚is the dimension of the space of parameters which describe the interfaces between the fluids, ' in P N ˚is the dimension of the space of inputs, ' out P N ˚is the dimension of the space of outputs. The goal is to construct with ML algorithms a third function

F : R in ÝÑ R out (2)
such that FpDpzqq « Epzq for all z P R par .

If the function D has a left inverse D ´1 : R in ÝÑ R par , then the best solution is to take F " EoD ´1 (the situation can be thought as similar to the one of auto-encoders, in ML language [START_REF] Goodfellow | Deep Learning[END_REF]). In what follows, we will use a ML software to construct a function F which realizes (3) to the best, without even questioning about the invertibility of D. We detail hereafter the functions D and E used in this work.

Straight lines

The first and main example directly comes from the VOF literature [START_REF] Debar | Fundamentals of the KRAKEN code[END_REF][START_REF] Youngs | Time-dependent multi-material flow with large fluid distortion, Numerical methods for fluid dynamics[END_REF][START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF]. For small mesh sizes, regular interfaces between fluids are asymptotically straight lines: that is why we consider only the latter limit case in this section. This example has a pedagogical virtue, and other cases will be variations around this theme. Straight lines

I θ,r px 1 , x 2 q " cos θx 1 `sin θx 2 ´r, (4)
are described by 2 parameters pθ, rq P r0, 2πq ˆR`, as illustrated in Figure 1. The straight line delimits two half planes P ´" tI θ,r px 1 , x 2 q ă 0u and P `" tI θ,r px 1 , x 2 q ą 0u .

The intersection of square cells C i, j with the first half plane tI θ,r px 1 , x 2 qu ă 0 yields the volume fractions

α i j " ż C i j XP ´d x 1 dx 2
which depend on the parameters θ and r. By definition the volume fractions are normalized between 0 and 1, that is α i j P r0, 1s. Considering the other half plane P `would result in calculating 1 ´αi j instead of α i j . To construct the inputs, we decide of a certain block p2N `1q ˆp2N `1q of cells, symmetric with respect to the central cell, and gather the volume fractions in one vector α " pα i, j q ´Nďi, jďN P r0, 1s p2N`1q 2 Ă R p2N`1q 2 .

With these notations, the function D is Dpθ, rq " α.

(

) 5
Then one has to decide the output which might be the periodic angle θ P r0, 2πq in this example. However forcing periodicity is not in the default implementation in Keras-TensorFlow. Therefore it is better for implementation purposes to work with 2 outputs which are the components pcos θ, sin θq of the direction, because the 2π periodicity of the angle is naturally taken into account and this procedure is easy to generalize in 3D. So our function E is Epθ, rq " pcos θ, sin θq .

By construction the outputs are also normalized. For this example ppar, in, outq " p2, p2N `1q 2 , 2q.

Basic properties of ML methods

To have a self-contained presentation, we summarize hereafter the basic mathematical features which are used to construct the function F defined in [START_REF] Aniszewski | PArallel, Robust, Interface Simulator[END_REF][START_REF] Basting | Extended ALE Method for fluid-structure interaction problems with large structural displacements[END_REF]. More theoretical material is proposed in the appendix to explain the good properties of ML algorithms for our purposes.

Functions

ML frameworks make use of high level recursive numerical functions, and three basic functions are sufficient to construct the function F. Linear algebra. Functions like y " WX `b are implemented for X P R n , W P R mˆn and b P R m .

Non linearity. Some non linear functions are implemented for x P R, such as the sigmoid σpxq " p1 `expp´xqq ´1 or the rectified linear unit (ReLU) function Rpxq " 1 2 px `|x|q " maxp0, xq.

Recursivity. Recursive calculations like W 2 RpW 1 X `b1 q `b2 or R `W2 RpW 1 X `b1 q `b2 ˘are implemented.

These functions can be used to construct other functions like maxpa, bq " a `Rpb ´aq and minpa, bq " ´maxp´a, ´bq " a ´Rpa ´bq. Many more linear and non linear functions (convolution neural network, maxpool, . . . ) are also implemented in Keras-TensorFlow [START_REF] Chollet | Deep learning with Python[END_REF]33]. Once a combination is chosen, the structure of the function F is known. For example a function with one hidden layer (parameters W 1 , b 1 ) is FpXq " W 2 RpW 1 X `b1 q `b2 , a function with two hidden layers is FpXq " W 3 pR `W2 RpW 1 X `b1 q `b2 ˘`b 3 , and so on and so forth. A key feature of modern ML software is the calculation of the coefficients W 1 , b 1 , . . . of a function F which is represented or approximated within the above structure. This is performed by optimization of a functional which encodes [START_REF] Aniszewski | PArallel, Robust, Interface Simulator[END_REF][START_REF] Basting | Extended ALE Method for fluid-structure interaction problems with large structural displacements[END_REF]. This optimization uses stochastic gradient descent with automatic differentiation for the calculation of the gradient. The ReLU function R is piecewise differentiable, its derivative is the Heaviside function almost everywhere: it allows symbolic differentiation with the chain rule to calculate the first gradient of functions defined recursively. Comprehensive references can be found in [START_REF] Chollet | Deep learning with Python[END_REF]33].

Elementary VOF advection and ML

We detail in this Section two simple reasons why the non linear ReLU function R is one of the best, better than a sigmoid σ, for application to VOF reconstruction on Cartesian grids.

The first reason is that many Finite Volume schemes with second-order accuracy use limiters as minmodpa, bq " " 0 for ab ď 0, minp|a|, |b|qsignpaq for ab ą 0.

One has the formula minmodpa, bq " Rpa ´Rpa ´bqq ´Rp´a ´Rpb ´aqq which shows that limiting techniques are also prone to be rewritten recursively with the ReLU function R. In view of the mathematical theory of hyperbolic equations [START_REF] Godlevski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF], it is also instructing to remark that the ReLU function R is deeply connected to Kruzkov entropies η k pxq " |x ´k| " Rpx ´kq `Rpk ´xq.

The second reason is based on the following geometric observation. Consider Figure 2 where, on the left and central parts, an interface splits a 2D cell with area ∆x 2 in 2 pieces. In terms of a reconstructed volume fraction, it corresponds to a function αpx, yq " " 1 for 0 ă x ă γ∆x, 0 ď γ ď 1, 0 for γ∆x ă x ă ∆x.

Consider that a velocity u ˚ą 0 is given and compute the flux at the right boundary during a time 0 ă ∆t ď ∆x u f p∆tq " In gray the value of the indicatrix function is 1, in white the value is 0. On the left part, the time is t 1 " ∆x 4u ˚: the flux is f 1 " 0 and the mean volume fraction is f 1 ∆x 2 {4 " 0. On the center part, t 2 " 3∆x 4u ˚: the flux is f 2 " ∆x 2 4 and the mean volume fraction is f 2 3∆x 2 {4 " 1 3 .

ż ∆x y"0 ż ∆x x"p1´βq∆x αpx, yqdxdy, β " u ˚∆t ∆x .
An exact calculation yields f p∆tq " ∆x 2 Rpβ ´γq that is f p∆tq " ∆xR pu ˚∆t ´γ∆xq .

One recognizes the SLIC method [START_REF] Noh | SLIC (Simple Line Interface Calculation[END_REF][START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF] or the antidiffusive scheme [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF]. The formula [START_REF] Cockburn | A priori error estimates for numerical methods for scalar conservation laws. III. Multidimensional flux-splitting monotone schemes on non-Cartesian grids[END_REF] shows that the rectified linear unit function has the ability to be exact for some elementary VOF procedure. On the right part of the Figure, the interface has an angle and a reasonable approximation is f p∆tq « ∆xα 1 R pu ˚∆t ´γ1 ∆xq ´∆xα 2 R pu ˚∆t ´γ2 ∆xq with α 1 , α 2 , γ 1 and γ 1 conveniently chosen.

More arguments in favor of using the function R are in the recent works [START_REF] Daubechies | Nonlinear Approximation and (Deep) ReLU Networks[END_REF][START_REF] Ward | AdaGrad stepsizes: sharp convergence over nonconvex landscapes, from any initialization[END_REF]. These basic observations are the reason why we use only the ReLU function R in our tests, because we believe it is more adapted.

Application to ML angle reconstruction for straight lines

We address the interpolation features of Keras-TensorFlow and measure its ability to construct the function F for the reconstruction of the angle. A reference, but with a completely different method, is the ELVIRA algorithm [START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF] which is exact for the reconstruction of the angle for straight line interfaces.

This section is a continuation of the pedagogical example Section 3.2. The vectors in the ML dataset B are constructed as follows. The distance to the origin, parameter r, is sampled uniformly in the range p´0.5 ? 2, 0.5 ? 2q (here 41 values): so the central cell is always crossed by the interface (even if at a corner of the cell in the extreme case). We sample the angle θ uniformly in r0, 2πq every degree (that is 360 values). It yields n θ ˆnr vectors in R 9 for 3 ˆ3 blocks, or in R 25 for 5 ˆ5 blocks. The volume fractions are calculated with numerical integration in 2D (N quadra ˆNquadra points with N quadra " 100). The parameter N quadra is the number of integration point per dimension and it yields a quadrature error

ǫ quadra « 1 N quadra " 1%. ( 8
)
It produces cardpBq vectors z i P R par . Then we apply the function D in (5) which calculates volume fractions: it yields cardpBq vectors pX i q 1ďiďcardpBq , X i " Dpz i q P R in , in " p2N `1q 2 .

We also apply the function E in (6) which calculates the normal to the interface: it yields cardpBq vectors

py i q 1ďiďcardpBq , y i " Epz i q P R out , par " 2.
Systematically, 80% of the data belongs to the training dataset while 20% of the data belongs to the test dataset: this choice is made at random. At the end of this stage we have a dataset for training

B train " " pX i , y i q, 1 ď i À 4 cardpBq 5 * 6
and a dataset for testing

B test " " pX i , y i q, 1 ď i À cardpBq 5 * .
We use a dense two-layer ML reconstruction: the number of neurons of the dense hidden layer is systematically equal to the size of the inputs (parameter in), for example in=9 neurons for the first line of Table 1. Theoretically we postulate that an approximation formula (see Appendix A) holds for the function

F : R in ÝÑ R par X ÝÑ y " FpXq
and we let Keras-TensorFlow find weights such that the least square error is quasi-minimized, in the context of implementation with batches, 1 cardpB train q cardpB train q ÿ

i"1

|y i ´FpX i q| 2 . (9)
In all our tests, we use the default parameters. The stochastic gradient descent runs with the Adam-optimizer [START_REF] Kingma | ADAM: a method for stochastic optimization[END_REF]. Here the batch size is 128. The number of epochs is 200. It yields Table 1 where the third column is the number of training/validation data. The loss error defined by ( 9) is in column 5. The validation error (equal to (9) but with the dataset for testing instead of the dataset for training) is in column 6. We observe that they are of the same order. Actually in all training sessions, we have observed that the loss error and the validation error are systematically of the same order, which is hopefully the sign that no overfitting or underfitting is attached to the ML treatment of our data. We will not report anymore on this feature in the rest of the paper, since it is not a restriction.

The last column is the mean value of the L 2 error on the test data (equal to the square root of the error-validation). Clearly, the direction vector is recovered with an error around 2%. The results are promising. Many tests have been performed with the 2D data 1 : it has been observed that the L 2 error can be decreased by increasing the number of dense hidden layers (that is by increasing the quality of the interpolator). But here the gain is marginal. Indeed making the additional hypothesis that the predicted direction has norm equal to one, the mean L 2 error, when it is small, scales like

dim
mean L 2 error « d 1 cardpB test q ř ˆˇˇˇc os θ true sin θ true ´ˇˇˇc os θ predic sin θ predic ˙2 « c 1 cardpB test q ř `θtrue ´θpredic ˘2 " ǫ θ .
In 2D, it means ǫ θ is around 1% in relative unit ǫ θ « 1%. [START_REF] Debar | Fundamentals of the KRAKEN code[END_REF] Considering that the numerical integration (8) has itself a maximal error ǫ quadra of the same order, it is probably meaningless to increase the accuracy, without increasing the accuracy of the ML dataset (by increasing the accuracy of the numerical integration).

1 A similar test has been performed in 3D. It is a confirmation in our framework of the 3D results in [START_REF] Qi | Computing curvature for volume of fluid methods using machine learning[END_REF]. Now blocks are 27 cells or 125 cells, and the angles pϕ, θq P r0, 2πq ˆr0, πs which determine the planes are sampled quasi-uniformly to avoid over sampling near the poles. The numerical integration is uniform (N quadra ˆNquadra ˆNquadra points with N quadra " 50). The results for the recovery of the direction vector psin θ cos ϕ, sin θ sin ϕ, cos θq for planes in 3D are displayed below and they are quantitatively the same as in 2D.

dim

More general geometries and definition of a dimensionless VOF-ML flux

In the following, we parametrize the geometry by arcs of circle and corners which are assembled to describe interfaces more general than just straight lines, even if the principles are exactly the same. We also explain how to define and train the normalized original VOF-ML flux which will be used in the remap stage of the compressible solver.

Arcs of circle

We refer in [START_REF] Qi | Computing curvature for volume of fluid methods using machine learning[END_REF] for the ML modeling of interfaces between incompressible fluids. Our implementation is slightly different because we ask that the modeling of curved lines degenerates to the one of straight lines if applied to a mesh with smaller and smaller mesh size (it allows natural mesh convergence tests). Consider the circle of Figure 3: the center is

A " rpcos θ, sin θq `dp´sin θ, cos θq ´Lpcos θ, sin θq, the radius is L ě 0 and the offset in the direction perpendicular to pcos θ, sin θq is d P R. Take a point x " px 1 , x 2 q P R2 inside the disk of radius L, that is |x ´A| 2 ă L 2 . The expansion is

px 1 ´r cos θ `d sin θ `L cos θq 2 `px 2 ´r sin θ ´d cos θ `L sin θq 2 ă L 2 .
Cancellation of some terms yields 2L pcos θx 1 `sin θx 2 ´rq `px 1 ´r cos θ `d sin θq 2 `px 2 ´r sin θ ´d cos θq 2 ă 0, that is

I θ,r px 1 , x 2 q `1 2L " px 1 ´r cos θ `d sin θq 2 `px 2 ´r sin θ ´d cos θq 2 ı ă 0.
Using a rescaling of all lengths x 1 , x 2 , r, d, L by a factor which corresponds to the mesh size ∆x, one gets

I θ,r px 1 , x 2 q `1 2 ν " px 1 ´r cos θ `d sin θq 2 `px 2 ´r sin θ ´d cos θq 2 ı ă 0, ( 11
)
where the factor is ν " ∆x L ą 0. For small mesh size ∆x, the continuous limit ν Ñ 0 recovers the equation (4) of straight lines 2 . Taking par " 4 and in " p2N `1q 2 as before, the function D : R par Ñ R in is defined by Dpθ, r, d, Rq " α. The function E will be described in Section 4.3 in formula [START_REF] Després | An explicit a priori estimate for a finite volume approximation of linear advection on non-Cartesian grids[END_REF].

Corners

The interfaces described above are smooth and so are not satisfactory for the local modeling of interfaces which are not C 1 but only Lipschitz, such as the cross test problem and the Zalezak test problem which will be considered in the numerical section. For such test problems, normal vectors at the interfaces have points of discontinuities: these points are called corners or corner points in this work. Consider the illustration in Figure 4. Straight lines correspond to µ " π (it is a degenerate case), right angles correspond to µ " π 2 and µ " 3 π 2 and acute angles correspond to the other values of µ. The corner is denoted as A " px A 1 , x A 2 q. Actually this modeling of corners allows also to model straight line interfaces, it is sufficient either to take µ " π or to sample A outside of the block with an additional condition on the aperture of the corner.

A " px A 1 , x A 2 q µ θ O Figure 4: Description of the parameters pθ, µ, x A 1 , x A 2 q for corners.
With these notations, one has par " 4 and in " p2N `1q 2 . The function D : R par Ñ R in is defined by Dpθ, µ, x A 1 , x A 2 q " α. The function E is described in the next Section 4.3 in formula [START_REF] Després | An explicit a priori estimate for a finite volume approximation of linear advection on non-Cartesian grids[END_REF].

Definition of a dimensionless VOF flux and adapted inputs

The function f in (7) used for the pedagogical case of straight lines is not normalized. For normalization reasons, it is better to consider the average flux which is the flux divided by the area of the swept region

gp∆tq " f p∆tq p1 ´βq∆x 2 . ( 12
)
This quantity is normalized in the sense that, on the one hand gp∆tq can be evaluated in function of the normalized quantities β P r0, 1s and γ P r0, 1s, and on the other hand gp∆tq P r0, 1s under the CFL condition u ˚∆t ∆x P r0, 1s. Other tests performed for the ML reconstruction of f p∆tq show a loss of accuracy for small velocity or for small Courant number: it was expected because f p∆tq P r0, ∆xu ˚∆ts is systematically small for small ∆t; on the contrary gp∆tq is much less sensitive to small ∆t. For all these reasons, we advocate using g instead of f as output for VOF flux reconstruction. This gives the function E " gp∆tq that will be used from now on with out " 1.

The reconstruction of the parameters of the different types of interfaces is not sufficient for the remap stage of CFD calculations. Indeed, even with an accurate reconstruction of these values, one must still construct the flux at the cell face. Therefore, to have a satisfactory ML modeling of the flux, we enlarge by one the size of the inputs. The new input is the rescaled velocity which scales like a non dimensional Courant number β " u ˚∆t{∆x P r0, 1s.

We add it to the volume fractions, so the function D maps the parameters to a vector made of the volume fractions plus β Dpinterface parametersq " pα, βq P R p2N`1q 2 `1.

(

) 13
For example for a 3 ˆ3 block, the size of the inputs is now in " 1 `9 " 10. The output size is out " 1: the output is the dimensionless outgoing flux

g " 1 β ż 1 2 1 2 ´β ż 1 2 ´1 2 αpx, yqdxdy P r0, 1s (14)
where α is the indicatrix function, as described in [START_REF] Després | An explicit a priori estimate for a finite volume approximation of linear advection on non-Cartesian grids[END_REF]. It yields the function E Epinterface parametersq " g P R.

(15)

Validation of the VOF-ML flux for straight lines

Redoing the tests of Table 1 for N " 3 and N " 5, the dimensions are now par " 2, in " 1 `p2N `1q 2 and out " 1. The size of the ML dataset is 100 ˆ360 ˆ41. The batch size is 16 ˆ1028. The number of epochs is still 200. We use 2 dense hidden layers: a first one with 4 ˆin neurons, a second one with 2 ˆin neurons. The results are in Table 2. One observes that the results in the last column are even slightly better than in Table 1 At the end of this procedure, all the geometry (interface reconstruction, numerical integration of indicatrix functions in normalized cells, numerical integration of the flux) is offline. That is, it is considered only in the Python-Keras-TensorFlow stage of the calculation chain.

Construction of ML datasets and accuracy of the VOF-ML flux for general geometries

The construction of the ML datasets relies so far on uniform sampling of all parameters. This was possible because the number of parameters was small, equal to 2 in this example.

But we need more parameters to describe complex interfaces like arcs of circles and corners, so uniform sampling is no more possible because the number of configurations would blow up. That is why we decided, quite arbitrarily, to keep the uniform sampling of the angle θ and of the rescaled velocity β, and to take random values for all other parameters. Now the size of the ML dataset is n θ ˆnβ ˆ2 ˆP: the factor 2 is because we systematically take the contrasted values (r α " 1 ´α) and the extra factor P guarantees that the same value of the parameters θ and β are used P times while the other parameters are randomized. The sampling of the rescaled velocity β is equal to the sampling for the numerical integration, because it is convenient for implementation purposes, that is n β " N quadra .

The result for the reconstruction of the flux for arcs of circle (Section 4.1) are in Finally we consider the corner configuration of Section 4.2. We take n θ " 24 ˆ360 uniformly distributed angles. We take n β " 100 values of β which are uniformly distributed. We double the results (taking the negative value α Ð 1´α and β Ð 1´β). The parameters x A , y A are taken random in the square r´3{4, 3{4s for in " 10 and in the square r´5{4, 5{4s for in " 26 . The parameter µ takes two different values π 2 and 3π 2 at random, it models only right angles. The results are displayed in Table 4. It seems the 3 ˆ3 block has more difficulty to recover the solution, even if the accuracy is still of the same order as previous accuracies. It is probably related to some lack of invertibility of the function D for 3 ˆ3 blocks (a topic already evoked at the end of Section 3.1). It is also probably related to the lower regularity of this type of Lipschitz interfaces.

Inference of the model in a C++ code

Now that ML datasets are built with various parameters, and that the training and testing are performed, it remains to pass the model to a C++ code for the inference, i.e. the F function call. The software keras2cpp offers the ability of calling the function many times on tensors (which are arrays)

T input P R N tensor ˆin . ( 16
)
The size of the rectangular tensor T is N tensor ˆin where N tensor is the number of fluxes that one desires to calculate and in is the size of the inputs. The function F then returns a tensor

T output " F pT input q P R N tensor ˆout ( 17
)
where out is the size of the outputs (actually equal to 1 in our case). This procedure allows a natural acceleration which is highly profitable for performance: a typical mean value of a call to F is between 10 ´5s and 10 ´4s, see the end of Section 8.

At the end of this stage, one has a function F callable in C++ which realizes the task (2-3) for the considered dataset: the output is the dimensionless flux (12).

Implementation in a Finite Volume solver

We provide some important details of the implementation of the method within a Finite Volume solver.

Preservation of natural symmetries

Calculations on a Cartesian mesh should in principle respect natural symmetry principles, when they are available. Moreover the accuracy of the ML interface reconstruction is achieved with an accuracy « 1% for the reason explained above. So any loss of symmetry can induce important discrepancies. The situation is different for incompressible fluid flow where the curvature of bubbles is recovered within an error of 10 ´6 with ML algorithms [37]. We particularize 3 different symmetries, and the consequences for the implementation of the VOF-ML scheme.

Rotational symmetry on a Cartesian mesh. The scheme must be invariant when a rotation of angle π{2 is applied. This is easy to obtain with a convenient ordering of the volume fractions. This principle is already used in standard Finite Volume implementations [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][Chapter 24], so there is no need to describe further. See Figure 5 for an illustration. Figure 5: A rotational symmetry is described for a 3 ˆ3 block. On the left the call is Fpβ, α 1 , α 2 , α 3 , α 4 , α 5 , α 6 , α 7 , α 8 , α 9 q. On the right it is Fpβ 1 , α 3 , α 6 , α 9 , α 2 , α 5 , α 8 , α 1 , α 4 , α 7 q where the ordering of the volume fractions is the same up to a rotation. A similar procedure is used for calculation of the flux across the right boundary and the top boundary.

Mirror symmetry. Consider transport in the left direction, as in Figure 5. Let Y " pα 1 , α 2 , α 3 , α 4 , α 5 , α 6 , α 7 , α 8 , α 9 q and consider its rearrangement ΠY " pα 7 , α 8 , α 9 , α 4 , α 5 , α 6 , α 1 , α 2 , α 3 q, where the bottom line components α 1 , α 2 , α 3 are exchanged with the top line components α 7 , α 8 , α 9 . The function should satisfy the mirror symmetry Fpβ, Yq " Fpβ, ΠYq, but such a symmetry is not guaranteed by our construction. So we force it by implementing

F mir pβ, Yq " 1 2 pFpβ, Yq `Fpβ, ΠYqq .
A similar treatment is applied for the left-right mirror symmetry.

Numbering symmetry. It means that the result should be the same if we exchange the role of material 1 and material 2 in the model, that is if we exchange a 1 and a 2 . Mathematically, it can be described by considering the vector e " p1, . . . , 1q and the function ΛY " e ´Y. We note the commutation ΠΛ " ΛΠ. The numbering symmetry corresponds to the property Fpβ, ΛYq " 1 ´Fpβ, Yq.

Since the construction does not guarantee this property, we force it by implementing

F numb pβ, Yq " 1 2 pFpβ, Yq `1 ´Fpβ, ΛYqq .
The scheme which is implemented satisfies all symmetries by using Fpβ, Yq ML VOF "

1 4 pFpβ, Yq `Fpβ, ΠYq `1 ´Fpβ, ΛYq `1 ´Fpβ, ΛΠYqq . ( 18
)
This flux is called the VOF-ML flux in the rest of this article3 .

Preservation of natural bounds

It is expected that the VOF-ML flux ( 18) can violate natural bounds like 0 ď a 1 ď 1, since no such considerations have been used so far. Of course, any kind of a posteriori method can be used to recover the maximum principle. In this work we use a simple a priori bound which is well adapted to directional splitting implementations. Essentially we constraint the flux so that no negative volume or mass is created.

Consider a square cell p0, ∆xq 2 with initial volumes 0 ď V 1 " a 1 ∆x 2 and 0 ď V 2 " a 2 ∆x 2 with a 1 `a2 " 1. Assume a flux a 1 is applied on one side for a velocity u ˚ą 0, and disregard the flux on the other side. The remaining volumes at the end of the time step will be V 1 " V 1 ´∆xpa 1 ∆xu ˚q and V 2 " V 2 ´∆xpa 2 ∆xu ˚q where a 2 " 1 ´a1 . Note that V 1 `V2 ă ∆x 2 , but this is not the important point. The conditions V 1 ě 0 and V 2 ě 0 yield two inequalities a 1 ď ∆x u ˚∆t a 1 and 1 ´a1 ď ∆x u ˚∆t p1 ´a1 q which are summarized in

1 ´∆x u ˚∆t p1 ´a1 q ď a 1 ď ∆x u ˚∆t a 1 .
This double bound is systematically added in our implementation of the VOF-ML flux. It yields a non negativity principle, which transfers into a maximum principle for the volume fractions (and mass fractions as well). In practice we check if the VOF-ML flux [START_REF] Godlevski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] belongs to the interval " 1 ´∆x u ˚∆t p1 ´a1 q, ∆x u ˚∆t a 1 ‰ : if it is in the interval we do nothing; if it is outside the interval, we replace the value of the flux with the nearest extreme value in the interval, which is equal either to the lowest value 1 ´∆x u ˚∆t p1 ´a1 q or to the highest one ∆x u ˚∆t a 1 .

Hybridation of the VOF flux

VOF schemes are often used in combination with other numerical techniques [START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF][START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF]. For example a VOF interface reconstruction technique can be used near the interface, while a more traditional Finite Volume scheme may be used away from the interface. This is why it is probably better to conceive that VOF-ML schemes must be hybrid ones. Hybrid schemes have also the virtue that, in principle, the higher numerical cost of sophisticated VOF schemes contributes to the total cost only for cells near the interface.

The schemes developed in this article for test purposes are hybrid schemes. Various sensors have been tested, but the simplest one is retained. A small threshold value ǫ flag ą 0 is taken and one checks at the beginning of every iteration and for all cells if a 1 is in between ǫ flag and 1 ´ǫflag : if ǫ flag ď a 1 ď 1 ´ǫflag ùñ the cell is flagged.

If the cell is flagged, one uses the VOF-ML flux [START_REF] Godlevski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]: in the other case, we use another much cheaper method; in the Euler2D code, the anti-dissipative SLIC/Downwind flux is available, it is cheap and with strong anti-dissipative properties [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF]. Considering that the quadrature error and average accuracy at the end of the training session are of the order of 1%, we take a slightly greater value ǫ flag " 2 10 ´2.

With this value, all thresholds (8), ( 10) and ( 19) have the same magnitude. It reflects a kind of compatibility of the whole chain of calculation.

An additivity principle for the ML datasets

The additivity principle states that the ML dataset used for computing a given profile must be the union (or the addition) of ML datasets dedicated to different local features of the profile. In mathematical words, it means the ML dataset uses different functions D 1 , D 2 , . . . , but the same function E, and that the function F is the best that ML can construct such that FpD i pzqq « Epzq for all z P R par and for all i " 1, 2, . . . [START_REF]Keras2cpp[END_REF] This additivity principle is more a "good practice" requirement, and no theory is available so far. However its efficiency is illustrated by the following simple numerical experiment. Consider the advection of a square. Locally the square is made of straight lines and corners. In the top left part of Figure 6, a simulation is displayed where the ML training is done with a dataset with straight line profiles on 5 ˆ5 blocks. It is clear that the numerical procedure is efficient for the capturing of straight lines, but the accuracy is very poor at corners. On the contrary, on the top right part of the Figure, the ML dataset is the addition of two ML datasets dedicated to straight lines and corners. Clearly, the overall accuracy is correct. On the bottom of Figure 6, the number of flagged cells is plotted in function of the number of iterations in time. On the left part, the number of flagged cells increases linearly. It renders the computation inefficient and inaccurate. On the right part of the Figure, with the good resolution ML dataset (that is with straight lines and corners), the number of flagged cells is bounded uniformly with respect to the iterations.

Not only the ML algorithm is able to find good parameters for the advection of straight lines and corners separately, but also it is able to interpolate smoothly from case to the other. More generally we have observed in all our calculations that this principle is true:

Additivity principle ùñ Flagged cells number is bounded ùñ Computational efficiency.

A simple bi-material hydrodynamics model

Numerical tests of the VOF-ML flux have been performed in the Euler2D code developed by F. Lagoutière [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF] that is described below.

In its basic setting, it uses a first-order accurate Lagrange+remap Finite Volume scheme with directional splitting [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF] on a 2D Cartesian grid. The global scheme is conservative for masses, total momentum and total energy, and is endowed with strong stability properties based on entropy inequalities [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF]. It solves a bi-material compressible Euler model $ ' ' & ' ' % B t pρc 1 q `∇ ¨pρc 1 uq " 0, B t pρc 2 q `∇ ¨pρc 2 uq " 0, B t pρuq `∇ ¨pρu b uq `∇p " 0, B t pρeq `∇ ¨pρue `puq " 0.

The additivity of mass fractions, specific volumes (1{ρ " τ) and internal energies is assumed

$ & % 1 " c 1 `c2 , τ " c 1 τ 1 `c2 τ 2 , e " c 1 ǫ 1 `c2 ǫ 2 `1 2 |u| 2 .
We use perfect gas pressure laws p 1 " pγ 1 ´1qρ 1 ǫ 1 and p 2 " pγ 2 ´1qρ 2 ǫ 2 .

For γ 1 " γ 2 , the model degenerates to a mono-fluid or mono-material model. Truly bi-fluid or bi-material configurations correspond to γ 1 ‰ γ 2 . The closure (that is the global pressure p) is calculated in accordance with, either the iso-pressure/isotemperature model (p 1 " p 2 , T 1 " T 2 , ε 1 " C v1 T 1 and ε 2 " C v2 T 2 ), or the iso-pressure/iso-δQ model

p 1 " p 2 , D t ǫ 1 `p1 D t τ 1 " D t ǫ 2 `p2 D t τ 2 .
The iso-pressure/iso-δQ model is our preference because the quality of the results is often better, but its implementation is non trivial because this closure model is non conservative. On the other hand the mathematical analysis of the iso-pressure/isotemperature model is simpler, moreover its implementation is easier and the code is more robust: these are valuable features for comparisons. We will show results with both closures. If the initial data is constant pressure ppt " 0, xq " p ref and constant velocity up0, xq " u ref , then the exact solution of the model (whatever the closure) is provided by the solution of pure advection at velocity u ref .

The volume fractions are related to the mass fractions and specific volumes

a 1 " c 1 τ 1 τ and a 2 " c 2 τ 2 τ " 1 ´a1 .
The boundary conditions can be periodic, neutral or wall conditions. Periodic boundary conditions are required for pure advection test cases. Two advection schemes are available in the Euler 2D for the remap stage. The first one is the first-order accurate upwind scheme which has poor accuracy for interfaces. The second one is the downwind scheme under upwind constraints [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF]: it is similar to the SLIC algorithm [START_REF] Noh | SLIC (Simple Line Interface Calculation[END_REF] and will be referred to as the SLIC/Downwind scheme. The other details of such a code are standard, we refer to [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF][START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF].

Numerical results

Here we report on the accuracy of the new VOF-ML procedure, by considering basic test problems of the literature calculated within Euler2D. We distinguish pure advection and shock problems.

In order to explain what is needed for the calculation of a given type of interface, we consider VOF-ML fluxes issued from different datasets. For the Zalezak notched circle problem shown in Figure 9, the dataset is complete in the sense that it contains arcs of circles, right corners and corners with angle. In the appendix, a dedicated and somewhat artificial dataset is specially designed for the Trifolium test problem. An accurate model callable with keras2cpp for a 5 ˆ5 block trained with a ML dataset with arcs of circles, right and acute angles can be downloaded at [13]. It can be used to reproduce the results.

Advection test cases

The initial data is pressure equilibrium and the vectorial velocities are the same in all cells u " u ref " pa, aq with a " 1. We use periodic boundary conditions. The domain of calculation is a square with characteristic length equal to 1. The exact solution at time t " 1 is also the initial solution. The schemes use directional splitting. With an appropriate choice of the initial data, the sound speed c is smaller than the material velocity in both dimension, that is c ă 1. The indicated Courant number is CFL " a ∆t ∆x . In our tests, the Courant number varies between 0.1 to 0.6 in order to explore the sensibility of the algorithm. It must be noticed that Courant numbers too close to 1 are of less interests since our advection scheme becomes perfect at the limit. Also the Courant number is one parameter in the vector of inputs (parameter β): since the training session can be interpreted as an interpolation procedure, it is expected that the quality of the interpolation is smaller at the boundary of the domain of inputs (in phase space), so for Courant numbers close to 0 or to 1. For Courant numbers close to 1, this phenomenon probably mitigates the fact that an upwind advection scheme is perfect at the limit.

The cross

This is an extension of the square test problem displayed in Figure 6. The Courant number is 0.4. The ML dataset is made of curved lines and corners. Actually when ν " 0, the curved lines (arcs of circle) dataset contains the straight lines dataset. So the accuracy for straight lines description is a priori as good as it is with only straight lines in the dataset. In the ML dataset we oversample the corners (1163196 items, obtained with 4 ˆ18 ˆ20 angles, only right angle corners that is µ " π{2 or µ " 3π{2 only, and 10 different values of the corners positions) with respect to arcs of circle (348754 items, obtained with 4 ˆ360 ˆ2 angles, other parameters are randomized). In this case, the oversampling of corners is a good way to obtain a satisfactory accuracy. The interpretation is that it counteracts the lower accuracy of corner reconstruction with respect to straight lines reconstruction (just compare the last column of Table 2 and4). Here we use 5 dense layers in the neural network (meaning we use 5 levels of recursivity as described in Section 3.3.1).

The results are displayed in Figure 7. The results are good, in particular for the 5 ˆ5 blocks. Looking in more detail at the top right part of the Figure where the result is computed with the SLIC/Downwind scheme available in the code [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF], one sees large errors because the branches of the cross are not aligned. This is not observed with the VOF-ML scheme.

The disk

The disk is discretized on a 40 ˆ40 cells grid. The Courant number is 0.4. The results with 3 ˆ3 blocks and 5 ˆ5 blocks are displayed in Figure 8. We also display the result obtained with the anti-dissipative scheme [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF] which is over-compressive: it quickly captures an octagon instead of a disk. We use of course the ML dataset with arcs of circle. The results with the VOF-ML schemes are qualitatively accurate.

The Zalezak notched circle

We consider the advection of the Zalezak notched circle [START_REF] Pilliod | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF] tilted π 7 with respect to the mesh axis. The elevations are displayed in Figure 9. The results with the SLIC/Downwind scheme are not accurate since the disk is transformed into a polygon (as for the disk). The results calculated with a 5 ˆ5 block are better but a discrepancy is visible at one of the acute angles. But, following the additivity principle, with right angles (µ " π 2 or 3π 2 ) and acute angles (parameter 0 ď µ ď π{8) in the dataset, the results are very accurate once again. Below in Figure 10 we reproduce the contour 0.5 of the Zalesak profile at time t " 1 and t " 3. The same setting is used (directional splitting, 15 cells in one radius). One observes that the VOF-ML scheme provides a good numerical rendering of the corners and acute angles at t " 1 even if a slight degradation is visible at t " 3.

Convergence tests

It is instructive to perform convergence tests to get some guarantee of the good quality of these schemes with very fine meshes. It must be mentioned that the schemes being written in a Finite Volume framework, they are conservative for local masses, total impulse and total energy. We consider the advection of 2 smooth interfaces which are a disk and an ellipse. The disk has a radius L " 0.4. The equation of the ellipse is 2x 2 `y2 " L 2 . In the dataset the parameter r is carefully interpolated in the interval r0, 1.5 ? 2s. The results are displayed in Table 5 where we measure the difference between the initial data and the final data at time t " 1 in the L 1 norm. One observes a low rate of convergence, however it is remarkable that, even on a coarse mesh, the level of error is quite small in L 1 norm. Then, by refining the grid, the L 1 norm diminishes at a rate which is approximatively 1/2: that is doubling the number of cells in every direction decreases the error by approximatively a factor ? 2. This feature is clear on Figure 11. Therefore one can consider that the new VOF-ML schemes have the attributes of a Finite Volume scheme: they are conservative and converge in L 1 norm at rate 1/2 for bounded variation (BV) data, that is the error scales like « C VOF´ML ∆x 1 2 . We refer to [START_REF] Cockburn | A priori error estimates for numerical methods for scalar conservation laws. III. Multidimensional flux-splitting monotone schemes on non-Cartesian grids[END_REF][START_REF] Després | An explicit a priori estimate for a finite volume approximation of linear advection on non-Cartesian grids[END_REF][START_REF] Merlet | L 8 -and L 2 -error estimates for a finite volume approximation of linear advection[END_REF] for estimates obtained in the context of standard numerical analysis and to [START_REF] Delarue | Probabilistic analysis of the upwind scheme for transport equations[END_REF] for an optimal estimate obtained with probabilistic tools which shows convergence of Finite Volume schemes at a rate h 1 2 in L 1 for initial data in BV. For the completeness of this work, Table 6 provides the L 1 error norm computed with the Upwind scheme and with the SLIC/Downwind scheme. For the Upwind scheme, one observes convergence rate « C Upwind ∆x 1 2 , with a constant C Upwind « 600ˆC VOF´ML which is much larger than the one of the VOF-ML scheme. For the SLIC/Downwind scheme, one observes a first order convergence rate « C SLIC ∆x: what is remarkable is that the constant C VOF´ML is so small so that the error of the SLIC/Downwind scheme for the 480 ˆ480 very fine mesh is still greater than the error of the VOF-ML scheme on the coarse 30 ˆ30 mesh.

In these tests, a rigorous protocol has been adopted. That is the parameters for the construction of the dataset for the 5 ˆ5 blocks are exactly the same as the parameters for the construction of the dataset for the 3 ˆ3 blocks. Also the training is done with 5 000 epochs. As a result, on the Table and the Figure, the accuracy is very similar for 3 ˆ3 blocks and 5 ˆ5 blocks. For the finer mesh, the accuracy is slightly in favor of 5 ˆ5 blocks. With less rigorous protocols for the construction of datasets, the comparison between 3 ˆ3 blocks and 5 ˆ5 blocks is less clear. Here we see that the results are of equal quality between 3 ˆ3 and 5 ˆ5 blocks. A simple explanation is that the boundaries are smooth and smooth curves tend to straight lines asymptotically for small mesh size. Since straight lines are already well captured with 3 ˆ3 blocks it is not surprising that the gain offered by 5 ˆ5 is marginal here. On the contrary, if the interfaces are less regular as in Table 4, it must be reminded that 5 ˆ5 blocks are more accurate than 3 ˆ3 blocks. 5 and L 1 convergence at a rate « 1{2 is satisfied asymptotically.

Shock test problems

Test problems with shocks are often simpler than pure advection, since the interfaces coincide with contact discontinuities which are constrained by the dynamics of the equations (in particular they are often less singular). . The 3 ˆ3 blocks and 5 ˆ5 blocks are run with a VOF-ML scheme trained with 5 dense layers, with a large number of epochs (5 000 here) to reach accuracy of the training. For reasons explained in Section 8.3, the finer the mesh, the smaller the relative cost of VOF-ML with respect to the Finite Volume solver. In column 5, the L 1 numerical order of convergence is o n " log e n´1 ´log en log 2

. In column 6, it is o n " log e n´2 ´log en log 4

. 30 60 120 240 480 order Upwind 0.6633 0.4895 0.3565 0.2551 0.1818 « 0.5 SLIC/Downwind 0.0527 0.0352 0.0195 0.0108 0.0051 « 1 Table 6: Reference L 1 norm of the error for the ellipse test case with the Upwind scheme and with SLIC/Downwind scheme. The best error with the SLIC/Downwind scheme on the finest mesh is still much greater than the error of the VOF-ML scheme (in Table 5) on the coarsest mesh.

such that maxp|u x |, |u y |, cq ∆t ∆x ď 1. For the last test case where the sound speed takes high values, then the Courant number for advection is small. For such simple shock problems, the SLIC/Downwind scheme already provides a satisfactory accuracy. Any dataset made only with arcs of circles and with fine training is sufficient. Here the interfaces are smooth curves. We used 5 ˆ5 blocks. These problems also contribute to the validation of the VOF-ML flux for a non-trivial velocity field that contains some degree of vorticity beyond pure translation or pure solid body rotation.

1D shock tube problem

We consider a symmetric 1D shock tube problem computed with the 2D code. The initial data come from the Sod shock tube problem but with

γ 1 ‰ γ 2 0.3 ă x ă 0.7 : a 1 " 1, a 2 " 0, ρ 1 " 1, p 1 " 1, u x " 0, u y " 0, γ 1 " 1.6,
x ă 0.3 or 0.7 ă x : a 1 " 0, a 2 " 1, ρ 2 " 0.125, p 2 " 0.1, u x " 0, u y " 0, γ 2 " 1.4.

To have good resolution even with a first-order accurate acoustic Lagrangian Riemann solver, we use 1000 ˆ5 cells. The Courant number is 0.6. The results calculated with the iso-pressure/iso δQ model are displayed in Figure 12. The steepness of the contact discontinuity is perfectly captured both on the density profile and on the volume fraction profile. In particular the perfect right-left symmetry of the profile is visible on the profile for a 1 .

2D convergent shock problem

We consider a convergent mild shock bi-material test problem where γ 1 " 1.4, γ 2 " 1.6. The isobar-isothermal closure model is used with C v1 " C v2 " 1. The initial data are p 1 " 0.1, p 2 " 1, ρ 1 " ρ 2 " 1. The shock is driven by the difference of pressure only. With ρ 2 ą ρ 1 the shock would be more violent, while with ρ 2 ă ρ 1 the shock would be weaker. The equality of density also simplifies our implementation of the initial condition. We take wall boundary conditions. The Courant number is 0.4. The interface is initially at R " 0.3. The mesh is made of 60 ˆ60 cells. We plot in Figure 13 the volume fraction at time t " 0, 0.1, 0.2 and 0.3. We observe a good rendering of the circular nature of the interface which evolves dynamically.

2D divergent shock problem

We exchange the pressure p 1 " 1, p 2 " 0.1 with respect to the previous test problem. We plot in Figure 14 the volume fraction at time t " 0, 0.1, 0.2 and 0.3. It also shows a good rendering of the steepness of the interface/contact discontinuity.

A Richmyer-Meshkov instability

We finally consider a 2D test problem with the data from [START_REF] Lagoutière | Numerical resolution of a two-component compressible fluid model with interfaces[END_REF]. A strong shock hits a one-mode sinusoidal interface which becomes unstable. The domain of calculation is r0, 3.6s ˆr0, 14s with an interface ypxq " 12 `0.5 cosp2πx{3.6q. The data on both sides of the interface are a 1 " 1, a 2 " 0, ρ 1 " 2.95, p 1 " 50000, u x " 0, u y " 453, γ 1 " 5{3, a 1 " 0, a 2 " 1, ρ 2 " 1.87, p 2 " 50000, u x " 0, u y " 453, γ 2 " 5{3.

The same coefficient γ 1 " γ 2 " 5{3 is used, which allows simple comparisons with purely mono-fluid implementations. The shock is created by a discontinuity at y " 7 The problem can be run with the iso-pressure/iso-δQ closure model, but it is simpler to run it with the iso-pressure/iso-temperature closure model by taking C v1 " 2 and C v2 " 3. With these values, the discontinuity of the density is very well recovered at the interface at the end of the simulation, with the advantage of a simpler and more robust closure model which can be used for comparisons (these important considerations about the choice of a closure model are however not directly related to the core of this work). The physics of this problem is interesting for our algorithms, since the length of the interface increases dynamically during the non linear stage of the instability. The number of square cells is 72 ˆ280 to mitigate the low accuracy of the first-order accurate Lagrangian phase. We represent the volume fraction a 1 at time t " 0.115 (CFL=0.5) in Figure 15: on the left part it is with the VOF-ML flux, on the right part it is with the SLIC/Downwind flux. The dynamics is the same even if they are some differences. The structure is narrower with VOF-ML. Actually the quality of results is also very sensitive to the CFL number: if it is too small, extra-diffusion occurs; if it is too large, some oscillations appear. Nevertheless these results shows that the VOF-ML flux has the ability to capture the dynamics of the interface even at the transition between the linear regime and the non linear regime. The symmetry of VOF-ML result is perfectly preserved by our implementation. The number of flagged cells displayed on the bottom-right part of the Figure is clearly related to the length of the interface. Indeed one recognizes the classical strong compression at iteration « 150 followed by a growth of the size of the interface in the linear range. Until the strong compression, the number of flagged cells is approximatively constant as in the right part of Figure 6. Then the number of flagged cells increases linearly, following the physics of the interface. We also plot the velocity in the y direction in the bottom-left part of the figure: since u y takes different signs, it naturally explains the stretching of the details of the volume fraction.

Computational costs

Finally we give measurements of the CPU costs for the VOF-ML schemes. We remind the reader it is in our tests implemented in C++ using the keras2cpp library. Clearly the CPU cost of one inference (that is one call of the function F) depends on our specific implementation and on the performance of, on the one hand the Euler2D code, and on the other hand the keras2cpp library: the Euler2D code and the keras2cpp library seem to have been implemented with the best practice for such pieces of C++ code.

We begin with a theoretical formula Time " `Ncell ˆTsolver `Nflagged ˆTtensor ˘ˆN cycle . Here T solver is the cost per cell per cycle of Euler2D, T tensor is the average cost per flagged cell per cycle of the VOF-ML schemes [START_REF] Godlevski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] when applied to rectangular tensors [START_REF] Fuster | An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF][START_REF] Gibou | Sharp interface approaches and deep learning techniques for multiphase flows[END_REF], N cell the total number of cells and N flagged the number of flagged cells. The total cost per cell per cycle is T total " T solver `σT tensor , σ " N flagged N cell P r0, 1s.

To determine the values of the various quantities on the right hand side of the equality, we make measurements (in µs per cell per cycle) on a standard Macbook Pro 2.9 GHz Intel Core i7. This is reported in Table 7.

The number of dense layers of the neural network used for measurements is indicated in the second column. For 3 ˆ3 blocks, the number of neurons per layer l ě 2 is decreasing function of the layer number, equal to 2 l´2 ˆin: the number of neurons of the last layer is par. For 5 ˆ5 blocks, the number of neurons per layer is the same (that is we use in " 10 and not in " 26), except for the first one for which it is equal to 2 l´2 ˆ26. The performance of the solver implemented in the Euler2D code is in the third column of the Table: one observes that the cost of the global Finite Volume Lagrange + remap scheme is approximatively 1 µs. It is representative of a reasonably optimized code. The mesh is 60 ˆ60 cells and the size of the tensors is N tensor " 100. What we observe by comparing the third and fourth columns of Table is that the cost of the VOF-ML scheme implemented with keras2cpp is higher than the cost of the Finite Volume scheme itself. For small test problems such as the ones considered in this work, it is not a real restriction. For more refined grids or for 3D problems, it would become a heavy burden, in particular if one uses a model with higher number of layers. But fortunately, interface calculations have the benefit of a good scaling. The good control of the number of cells in the neighborhood of the interface, like the ones we show in the numerical results produced in this work, has the effect that the extra cost of the VOF-ML scheme is asymptotically negligible for large number of cells. Indeed in dimension d " 2, one has N cell " N 2 1D and N flagged « λ interface N 1D where λ interface ą 0 is a constant which depends on the type of interface and on the numerical parameters, but which is independent of N 1D . Therefore one can rewrite the theoretical formula as T total " T solver `λinterface ˆTtensor N 1D .

Provided λ interface ˆTtensor is bounded uniformly with respect to N 1D , one obtains that T total « T solver for large N 1D .

It is the result (1) announced in the introduction.

Conclusion, perspectives and open problems

In this work we detailed a calculation chain for development of a family of Volume of Fluid-Machine Learning algorithms suitable for bi-material compressible fluid flow calculations on Cartesian grids. We explained a key feature which is the adaptation of the compressible fluid flow solver in order to respect the natural symmetries on a Cartesian grid. The last column gives the relative cost of the VOF-ML scheme, after multiplication of T tensor by an arbitrary scaling parameter σ " 100{60 2 which models the relative size of flagged cells near an interface relatively to the total number of cells.

The results show that, using the additivity principle for the construction of the ML dataset, it is possible to incorporate different features of curves (in terms of regularity), and that the final schemes interpolate between the various types of profiles with no noticeable difficulties. A variety of basic test problems common to the compressible community has been computed with satisfactory accuracy. For non smooth interfaces like corners, 5 ˆ5 blocks seem to be more accurate than 3 ˆ3 blocks, but for smooth interfaces, 5 ˆ5 and 3 ˆ3 blocks exhibit the same accuracy. To reproduce the results, an accurate model callable with keras2cpp for a 5 ˆ5 block trained with a ML dataset with arcs of circles, right and acute angles can be downloaded at [13].

A main issue of the training method is its global accuracy. To compute the numerical values of the volume fractions, one needs a numerical integration procedure adapted to different kind of interfaces (regular, right angle corner, acute angle, non Lipschitz triple point). That is why the Riemann integration with first-order 100 ˆ100 uniform sampling is used throughout this work. It yields an accuracy 10 ´2 which is an acceptable scale of error for a compressible calculation: all the rest of the simulation chain follows the same level of accuracy. In principle, it is possible to start with a better accuracy (10 ´3 for example). But it would result in a huge overhead for numerical integration of profiles with only Lipschitz regularity (not to speak of the Trifolium triple point and alike). In this direction, acceleration techniques available for smooth profiles would be a help [START_REF] Aniszewski | PArallel, Robust, Interface Simulator[END_REF] but with limitations. However it is known that some level of noise is sometimes better than no noise, because it yields robustness. This is why it is not guaranteed that a better global accuracy (10 ´3 for example) would result in any gain for compressible CFD with real data. New developments are necessary to investigate these delicate issues.

The main limitation so far is the overhead of the Volume of Fluid-Machine Learning calculation with respect to the flow solver. However the theoretical scaling shows the overhead is negligible for small mesh size. With dedicated hardware support, it is expected that advanced C++ implementations of Machine Learning functions will reduce this burden. So this limitation is expected to be less stringent in a near future.

The potential of the method for the Volume of Fluid-Machine Learning calculation of interfaces with low regularity seems rich. We have in mind to apply the method to multiphase and multi-material problems (3 phases and more, 2 phases near a wall, 2 materials sliding on a 3rd one, . . . ) and to tridimensional configurations for which traditional techniques based on regular approximation of interfaces suffer huge limitations.

Among open theoretical problems raised by our numerical experiments, we mention two. A first one is to understand for what reason the dynamical system made of the Volume of Fluid-Machine Learning flux coupled with a more traditional Finite Volume bi-material solver is able to get a good control of the number of flagged cells near the interface. As shown in Section 6.4, this possibility depends crucially of the dataset. Therefore, to get a theoretical answer to this question, it will be necessary to understand the algorithm globally. A second theoretical question is to prove the 1/2 convergence rate for bounded variation data, as it is observed in Section 8.1.4. So far our work is restricted to Cartesian meshes. An open problem of great practical interest is to go beyond such meshes of simple structure. The main difficulty is that adding the local parameters of the mesh in the vector of inputs might result in a very expensive method.

Appendix A.1. Interpolation of functions in dimension one

Take a smooth function f : R Ñ R with compact support. Note the Heaviside function Hpxq " 1 for x ą 0 and Hpxq " 0 otherwise. One has the formula f pxq " ş x ´8 f 1 pyqdy " ş `8 ´8 f 1 pyqHpx ´yqdy. Since R 1 " H almost everywhere, an integration by part yields f pxq " ş `8 ´8 f 2 pyqRpx ´yqdy. A standard discretization procedure (mesh size h ą 0) yields the approximation f h pxq " ř iPZ f 2 pihqRpx ´ihq. This formula shows that linear combinations of offsets of the function R can approximate the function f as accurately as desired. If the function f is not sufficiently smooth (for example if is piece-wise constant), then it is sufficient to regularize it by a factor ǫ ą 0 (the function f ǫ is now smooth with enough derivatives)) to get another approximation

f ǫ h pxq " ÿ iPZ f 2 ǫ pihqRpx ´ihqh. (A.1)
This formula can be recovered with one dense hidden layer of neural network, using the rectified linear unit function.

Appendix A.2. Interpolation of functions in higher dimensions

Similar identities hold in any dimension, and it provides a simple approximate version of the Cybenko theorem [START_REF] Cybenko | Approximation by Superpositions of a Sigmoidal Function[END_REF]. Let f : R d Ñ R be a smooth function with enough derivatives in L 2 pR d q. If it is not the case, it is always possible to regularize it. One has the general formula

Appendix A.3. Recursive construction of "hat" elements

The previous formulas did not make use of the recursive features. Now we state a known result [START_REF] Daubechies | Nonlinear Approximation and (Deep) ReLU Networks[END_REF] which shows that recursivity helps to reconstruct basic finite elements shape functions. In dimension one, the Finite Element hat function ϕ j is shows the interest of having a recursive ML construction. Generalization of this principle in higher dimensions is left to the reader. in Figure 15) in Section 8.2.4. Inspired by many works in the literature, see [START_REF] Ahn | Multi-material interface reconstruction on generalized polyhedral meshes[END_REF] and references therein, we add one more test problem with a vorticity field beyond pure solid body rotation and with simple analytical solution. The domain is r0, 1s ˆr0, 1s. We note r " a |x ´0.75| 2 `|y ´0.5| 2 . The initial data is a 1 " 1 in the disk tr ă 0.15u and a 1 " 0 outside the disk. The free-divergence velocity field u " 8π a x 2 `y2 p´y, xq is constant in time for 0 ă t ă 1 2 t end , and is reversed for 1 2 t end ă t ă t end . Due to the term r, the velocity field can be seen as the composition of a pure body rotation with an additional distortion which has a "slicing effect" on the solutions. The reversal of the velocity field has the effect that the exact solution at t " t end is equal to the initial one. We take t end " 0.6 and CFL " 0.2. The blocks are 5 ˆ5 cells. The VOF-ML flux is the one that was carefully trained for the convergence tests of Section 8.1.4. We remind the reader that the code is run with the splitting method (x, y, x, y,. . . ) which a priori is not the best one for rotation dominated problems. We first take a 60 ˆ60 cells discretization.

The results are displayed in Figure C.17. One observes that the SLIC/Downwind scheme performs already quite well, however the VOF-ML flux provides smoother results, it is visible on the isolines. The L 1 norm of the error at final time is smaller with the VOF-ML flux than it is with the SLIC/Downwind flux, that is err VOF´ML 2 t end , even if the VOF-ML flux is perhaps better at keeping the connection of the plateau. At final time t " t end , the results are of bad quality with the SLIC/Downwind flux, because the connectivity is lost. On the other hand the VOF-ML flux shows some ability to recover the connectivity of the correct profile, even if a corner is visible: our interpretation is that the intermediate profiles are smoother,

Figure 1 :

 1 Figure 1: Description of the parameters pθ, rq for straight lines.

1 0Figure 2 :

 12 Figure2: Description in a square cell of the swept region delimited by a moving interface. The swept region (in dashed) depends on the time t. In gray the value of the indicatrix function is 1, in white the value is 0. On the left part, the time is t 1 " ∆x 4u ˚: the flux is f 1 " 0 and the mean volume fraction is f 1 ∆x 2 {4 " 0. On the center part, t 2 " 3∆x 4u ˚: the flux is f 2 " ∆x 2 4 and the mean volume fraction is f 2 3∆x 2 {4 " 1 3 .

Figure 3 :

 3 Figure 3: Description of the parameters pθ, r, d, Lq for arcs of circle.

Figure 6 :

 6 Figure 6: Value of volume fraction a 1 at time t " 0.1. Top left: the ML dataset has only straight lines. Top right: the ML dataset is made of straight lines and corners, following the addition principle. Bottom: the number of flagged cells versus number of time steps.

Figure 7 :

 7 Figure 7: Value of volume fraction a 1 at time t " 1. Top left initial data (also equal to final data). Top right, SLIC/Downwind scheme. Bottom left, blocks are 3 ˆ3. Bottom right, blocks are 5 ˆ5.

Figure 8 :

 8 Figure 8: Value of volume fraction a 1 at time t " 1 after one round turn. Top left, initial and final data. Top right, SLIC/Downwind scheme. Bottom left, blocks are 3 ˆ3. Bottom right, blocks are 5 ˆ5.

Figure 9 :

 9 Figure 9: Value of volume fraction a 1 at time t " 1 after one round turn. Top left, initial and final data. Top right, SLIC/Downwind scheme: the notched circle is transformed in a polygonal line as shown by the contour lines. Bottom left, blocks are 5 ˆ5, but the ML dataset does not have acute angles. Bottom right, blocks are 5 ˆ5, now the ML dataset has acute angles. The level set of the bottom right result is presented in Figure 10.

Figure 10 :

 10 Figure10: Level set 0.5 of the advected Zalezak profile. Value of a 1 at time t " 1 on the left and t " 3 on the right. A small degradation of the quality of the isoline at one internal corner is visible at t " 3.

Figure 11 :

 11 Figure 11: Logarithmic plot of the L 1 error versus the number of cells. The data are from Table5and L 1 convergence at a rate « 1{2 is satisfied asymptotically.

Figure 12 :

 12 Figure 12: T " 0.1. Top left: pressure. Top right velocity. Bottom left: density. Bottom right: volume fraction.

ρ shock 2 " 6 .01, p shock 2 "Figure 13 :

 26213 Figure 13: Interface/contact discontinuity in a bi-material convergent test problem. From top left to bottom right, the time is T " 0, 0.1, 0.2 and 0.3.

Figure 14 :

 14 Figure 14: Interface/contact discontinuity in a bi-material divergent test problem. From top left to bottom right, the time is T " 0, 0.1, 0.2 and 0.3.

Figure 15 :

 15 Figure 15: Top: volume fraction a 1 for the Richtmyer-Meshkov test problem in a smaller part of the whole domain r0, 3.6s ˆr6.2, 10.2s Ă r0, 3.6s ˆr0, 14s; top left with the VOF-ML flux; top right with the SLIC/Downwind flux. Bottom left: the velocity field u y . Bottom right: number of flagged cells versus the number of time steps.

  ¨x ´bqdΩdb where S d " t|Ω| " 1u Ă R d and w : S d ˆR Ñ R is a function to identify. To construct w, differentiate f pxq " ş ΩPS d ş băΩ¨x wpΩ, bqdΩdb. It yields ∇ f pxq " ´şΩPS d wpΩ, Ω ¨xqΩdΩ. Take w as the half Fourier transform of another function z, that is wpΩ, aq " ş µą0 e iµa zpΩ, µqdµ. One gets ∇ f pxq " ´żΩPS d ż µą0 e iµΩ¨x zpΩ, µqΩdµdΩ. Make the change of variable k " Ωµ that is µ " |k| and Ω " k |k| . The equivalence of the measures writes dk " µ d´1 dµdΩ. One gets ∇ f pxq " ż kPR d ˆ´zpΩ, µq µ d ˙ke ik¨x dk. By identification with the formula ∇ f pxq " i 2π ş kPR d p f pkqke ik¨x dk, it shows that zpΩ, µq " ´i 2π µ d p f pµΩq. With this formula, one can get w, and finally represent f .

for p j ´1qh ď x ď jh, x´jh h for jh ď x ď p j ` 1qh .

 1qh One has ϕ j pxq " R `1 h Rpx ´p j ´1qhq ´2 h Rp jh ´xq ˘. So one can also approximate a given function with the well known Finite Element interpolation formula f hh pxq "

1 Figure C. 17 :

 117 Figure C.17:The volume fraction at mid-time1 2 t end and final time t end (where the initial profile is reconstructed). The velocity field is u " ˘8rπp´y, xq . Top row with the VOF-ML flux. Bottom row with the SLIC/Downwind flux.

1 " 2 " 1 "

 121 0.055. It is similar in quadratic norm, that is err VOF´ML L 0.140. The enhanced accuracy of the VOF-ML flux is a confirmation of the previous results. Taking a more stringent velocity field u " 8πpx 2 `y2 qp´y, xq and a coarser mesh made of 40 ˆ40 square cells results in a more severe test problem. It also illustrates the natural limits of VOF methods. The results in Figure C.18 show that the VOF-ML flux and the SLIC/Downwind flux have difficulties to capture the tail of the filament at mid-time t " 1

Table 1 :

 1

		par	in out	B train /B test	loss	val	mean L 2 error
	2	2	9	2	11807/2953 5.4 e-4 5.5 e-4	2.3 e-2
	2	2	25	2	11869/2891 1.8 e-5 1.7 e-5	4.2 e-3

Accuracy of the direction vector pcos θ, sin θq for straight lines in 2D.

Table 2 :

 2 . It validates the use of the mean flux as the output. Accuracy of the VOF-ML flux for straight lines.

	dim par	in out	B train /B test	loss	val	mean L 2 error
	2	2	10	1	1192451/298309	7e-5	7.1 e-5	9.8e-3
	2	2	26	1	1193738/297022 8.8e-5 9.4e-5	6.5e-3

Table 3

 3

	dim par	in out	B train /B test	loss	val	mean L 2 error
	2	5	10	1	1395790/349490 4.2e-5	4e-5	6.3e-3
	2	5	26	1	1395422/349858 3.1e-5 4.2e-5	6.5e-3

, where n θ " 8 ˆ360, n β " N quadra " 100, P " 3 and we randomize the parameters r ď 2 ? 2, d ď b N `1 2 and L ď 2 ? 2. The error (last column) is still considered small enough.

Table 3 :

 3 Accuracy of the VOF-ML flux for arcs of circle.

Table 4 :

 4 Accuracy of the VOF-ML flux for corners (right angles).

	dim par	in out	B train /B test	loss	val	mean L 2 error
	2	5	10	1	1396604/348676 1.9e-3 2.1e-3	4.5e-2
	2	5	26	1	1396744/348536 1.6e-4 2.2e-4	1.5e-2

Table 5 :

 5 Convergence table

Table 7 :

 7 Performance measurements in µs.

	block layers	T solver T tensor rescaled T tensor
	3 ˆ3 3 ˆ3 3 ˆ3 3 ˆ3 5 ˆ5 5 ˆ5 5 ˆ5 5 ˆ5	2 3 4 5 2 3 4 5	0.958 0.958 0.9807 90.8 9.78 31.9 0.977 255 0.946 15.5 0.928 44.7 0.962 113 1.007 304	0.271 0.886 2.52 7.08 0.432 1.24 3.134 8.44

If one desires to capture instead the small radius limit, that is L Ñ 0 `, it is better to rescale the equations. Then take σ " 1 ν " L ∆x and the parametrization 2σI θ,r px 1 , x 2 q `px 1 ´r cos θ `d sin θq 2 `px 2 ´r sin θ ´d cos θq 2 ă 0.

The cost of this procedure is

times the cost of one individual call to the function F. On a Cartesian grid in dimension d with m different materials, the cost of a similar procedure would be equal to 2 d´1 m! times the cost of one individual call.

Acknowledgements. The first author thanks Olivier Pironneau for many fruitful discussions during the elaboration stage of this work and acknowledges the support of CEA. The latest version of Euler2D was provided by Frédéric Lagoutière who is kindly thanked.

Appendix A. Interpolation features of ML algorithms

For the completeness of this work, we detail hereafter basic interpolation properties of ML algorithms which justify some of our choices.

Appendix B. The Trifolium pure advection test problem

This test problem illustrates the ability of VOF-ML schemes to treat some non Lipschitz profiles which are difficult to conceptualize with standard methods. We consider the pure advection u " p1, 1q of an interface defined in polar coordinates as r 3 " 1 27 cos 3θ. It is a modification with more pronounced lobes of the Trifolium rosace r " cos 3θ. The Trifolium has a triple point where the interface is non Lipshitz (even if it is piecewise Lipshitz). It also has a symmetry with angle 2π 3 which is not tackled by the natural symmetries of a Cartesian grid exposed in Section 6.1.

The results computed with a ML dataset made only of arcs of circle are displayed on the left of Figure B.16, calculated on a mesh with 100ˆ100 cells. They show that the Trifolium is already accurately and robustly recovered. However the structure of the triple point is lost, in particular its connectivity because two lobes stay connected and the third one is disconnected. The 2π 3 perfect symmetry near the triple point is not respected. So one can think of a new somewhat artificial ML dataset with a specific treatment of the triple point of the Trifolium. Following the additivity principle, a new ML dataset is created with Trifolium profiles 0 ď ǫr 3 ď cosp3θq where ǫ is a scaling parameter which is drawn randomly (0 ď ǫ ď 1{20 in our tests). In the right part of Figure B.16, we represent the new results. We observe that the connectivity of the 3 lobes and the 2π 3 symmetry of the Trifolium are now preserved with the VOF-ML scheme constructed from the full dataset. The Courant number is 0.1.

Appendix C. A test problem with vorticity field beyond pure solid body rotation

The design principle of the VOF-ML flux is based on the incorporation of local pure translations/advections in datasets. However it must be noticed that the method has been validated for test problems which display velocity fields which are quite different from global pure translation. Indeed the dynamically created velocity vector field for the compressible Euler equations is convergent in Section 8.2.2, divergent in Section 8.2.3 and more complex with a vortex (which can be seen as a double shear which is ultimately what is important for the quality at final time. With filament profiles in the dataset, it is possible that the results would be better (this is a completely open question).