# A 3D Machine Learning based Volume Of Fluid scheme without explicit interface reconstruction

Moreno Pintore, Bruno Després

<figure><img src="image_1.png" /></figure>

## To cite this version:

Moreno Pintore, Bruno Després. A 3D Machine Learning based Volume Of Fluid scheme without explicit interface reconstruction. 2025. hal-05149322

HAL Id: hal-05149322

https://hal.sorbonne-universite.fr/hal-05149322v1

Preprint submitted on 7 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A 3D MACHINE LEARNING BASED VOLUME OF FLUID SCHEME WITHOUT EXPLICIT INTERFACE RECONSTRUCTION

MORENO PINTORE\* AND BRUNO DESPRÉS†

**Abstract.** We present a machine-learning based Volume Of Fluid method to simulate multi-material flows on three-dimensional domains. One of the novelties of the method is that the flux fraction is computed by evaluating a previously trained neural network and without explicitly reconstructing any local interface approximating the exact one. The network is trained on a purely synthetic dataset generated by randomly sampling numerous local interfaces and which can be adapted to improve the scheme on less regular interfaces when needed. Several strategies to ensure the efficiency of the method and the satisfaction of physical constraints and properties are suggested and formalized. Numerical results on the advection equation are provided to show the performance of the method. We observe numerical convergence as the size of the mesh tends to zero $h = 1/N_h \searrow 0$, with a better rate than two reference schemes.

**Key words.** Volume Of Fluid, Machine Learning, CFD, 3D, Interface

**MSC codes.** 35Q35, 68T07, 76-10, 76M12.

**1. Introduction.** Multi-material or multi-phase numerical simulations are crucial in numerous engineering applications, a non-exhaustive list of examples includes water waves, nuclear reactors, jets and sprays, oil or gas pipelines, thermal power stations, and mold filling applications. Due to the presence of sharp interfaces between the involved materials, performing accurate simulations is complex and requires ad-hoc techniques. Indeed, a numerical scheme with too much diffusion would not be able to capture the sharp and discontinuous nature of the interface, whereas a scheme with less diffusion might not be stable enough to properly simulate the evolution of the phenomenon of interest. In this context, the Volume Of Fluid (VOF) method is one of the most commonly used techniques when exact conservation of one or more physical quantities is needed. Even though the first papers on the VOF method are not recent, see [17] and later references, the research field is still very active because of the importance of the applications and the complexity of the problem, which further increases for three-dimensional problems involving more than two materials. See, for example, [25, 34, 16, 39] and references therein. VOF methods are largely based on geometrical heuristic algorithms, even though a recent work establishes error estimates in some cases [7]. The essence of most VOF methods we know about (see our reference list) is to reconstruct some interface parameters from cell averages, and then to use this information to design numerical fluxes in a Finite Volume scheme.

Our work in three dimensions is based on the coupled VOF and Machine Learning methodology (coined as VOFML) for two-dimensional interfaces proposed in [9, 2]. Its main feature, which is also its main originality, is that we do not reconstruct interface parameters. Instead, we reconstruct the flux directly from the cell averages. The contribution of this work is to provide evidence that such a 3D VOFML is efficient. Since we do not explicitly reconstruct three-dimensional interface parameters, it results in a tremendous conceptual simplification with respect to most of the literature on VOF schemes with or without Machine Learning, for which we refer the reader to [17, 13, 30, 35, 12, 18, 29, 6] and references therein. In our opinion, this conceptual simplification is an example of the kind of new algorithms that Machine Learning can contribute to generate.

The general context of VOF methods is easily explained on our model problem. It is a practical simple scenario with two incompressible immiscible fluids, fluid A and fluid B, in the space-time domain $\Omega \times [0, T]$,

<footnote>\*Laboratoire Jacques-Louis Lions, LJLL, Sorbonne Université; MEGAVOLT Team, INRIA; F-75005 Paris, France, (moreno.pintore@sorbonne-universite.fr).</footnote>

<footnote>†Laboratoire Jacques-Louis Lions, LJLL, Sorbonne Université; Université Paris Cité; CNRS; MEGAVOLT Team, INRIA; F-75005 Paris, France, (bruno.despres@sorbonne-universite.fr).</footnote>
with $\Omega \in \mathbb{R}^3$ and $T > 0$. Let $\alpha : \Omega \times [0, T] \to \mathbb{R}$ be the characteristic function of the fluid A, i.e. $\alpha(x, y, z, t) = 1$ if the position $(x, y, z) \in \Omega$ is occupied by the fluid A at time $t \in [0, T]$, and $\alpha(x, y, z, t) = 0$ otherwise. Because of the immiscible and incompressible nature of the fluids, the conservation of mass implies the following conservation equation for $\alpha$:

$$ (1.1) \qquad \partial_t \alpha + \nabla \cdot (\alpha u) = \partial_t \alpha + u \cdot \nabla \alpha u = 0, $$

where $u$ is a divergence-free velocity field $\nabla \cdot u = 0$. In order to numerically solve (1.1), one introduces a mesh on $\Omega$ and associates with each cell $C$ of such a mesh the volume fraction

$$ \alpha_C = \frac{1}{|C|} \int_C \alpha \, dx \, dy \, dz, $$

where $|C|$ is the volume of the cell $C$. Note that $\alpha_C = 1$ if $C$ contains only fluid A, $\alpha_C = 0$ if it contains only fluid B, and $\alpha_C \in (0, 1)$ if it contains both fluids. The idea of the VOF method is to use the volume fractions available at time $t$ to compute the volume fractions at time $t + \Delta t$ without explicitly knowing the distribution of the two fluids, and to iterate this procedure to cover the entire time interval $[0, T]$ to mimic the evolution of the unknown function $\alpha$. The incompressibility hypothesis can be removed for real applications.

In order to compute the fluxes required to compute the volume fractions at the subsequent time step, most VOF schemes explicitly reconstruct, for each cell $C$, a simple interface that approximates the exact unknown interface between the two fluids. This simple interface is computed using the volume fractions of the cell $C$ and its neighbor cells (such a set of cells is named stencil), and it is often represented by affine or quadratic functions of $\mathbf{x}$. In particular, the most popular class of VOF methods adopts the Piecewise Linear Interface Calculation (PLIC) approach, where the interface is represented by a planar polygon in each cell. The most famous PLIC representatives are the Youngs method [37], where the plane normal $\mathbf{n}_\alpha$ is computed by using the relation $\mathbf{n}_\alpha = \frac{\nabla \alpha}{|\nabla \alpha|}$ and its location is computed by matching the volume fraction inside $C$, the LVIRA method (Least-square Volume-of-fluid Interface Reconstruction Algorithm) and its efficient variant ELVIRA (Efficient LVIRA) [30], where the planar interface is computed by solving a least squares minimization problem in each cell, and the height function method [35], where the interface normal is efficiently computed using a height-function technique. Because of the longevity of VOF methods, several alternatives exist. For example, it is possible to avoid the use of a stencil by introducing additional variables in each cell that represent the position of the center of mass of the fluid inside the cell, as in the Moment Of Fluid method [12], or to bypass the interface reconstruction by using explicit algebraic formulas [10] to directly compute the required fluxes.

When considering three-dimensional simulations, however, even the implementation of standard PLIC methods on Cartesian meshes is complex and it is often advisable to introduce further approximations [18] to simplify the implementation and reduce the associated computational cost. The situation is even more complex when more than two fluids are present in the simulation, because simple interfaces are not expressive enough to properly represent the complex distributions of more than two fluids in a single cell. For example, in [5] an additional minimization procedure in each cell with 3 fluids is considered to localize the triple point, and in [29] a different minimization strategy is adopted, assuming that the fluid ordering is known a priori, to localize each interface. Such minimization strategies often rely on iterative clipping of polygons, which is computationally expensive in 3D. Therefore, it is often necessary to consider alternative and more complex implementations like the one proposed in [21] to reduce the computational overhead. Because of all these sources of complexity and thanks to the rapidly increasing popularity of Scientific Machine Learning, in the last few years, numerous techniques based on neural networks [15] have been proposed. For example, in [6] and in [33] two PLIC variants in which the interface normal is computed through a neural network are proposed, in [27] a graph neural network is used to describe the local interface in unstructured meshes, in [32] and in [38] the neural network is involved only in the curvature computation, whereas in [24] the authors propose a
neural network to directly predict the volume fractions at the subsequent time step without computing and using the fluxes.

In this manuscript, we focus on extending to three-dimensional simulations the Volume Of Fluid - Machine Learning (VOF-ML) method proposed in [9], which is an alternative VOF method relying on a neural network. As discussed in [9] for two-dimensional problems, the main idea of VOF-ML is to train a neural network to predict, given the volume fractions in a suitable stencil, the flux that has to be used in a standard Finite Volume (FV) scheme [22]. Moreover, no simulations or physical measurements are required to train the network, because the training dataset is cheaply assembled from geometrical considerations. In fact, the training dataset contains several distributions of pairs of fluids in a stencil and the corresponding fluxes, and this information does not depend on the equation of interest, the type or number of materials, or the domain $\Omega$. All this additional information is only used by the FV scheme that uses the flux predicted by the neural network. For example, even if the network is trained on pairs of fluids, simulations with more than two fluids can be performed as discussed in [2] without retraining the neural network on distributions including more materials. The manuscript contains the description of the method and a few implementation details that can be used to enforce physical properties and more accurately compute the fluxes.

To summarize, we provide a list of features that characterize the proposed method:

* it addresses 3D configurations;

* it exploits the capabilities and nonlinearity of deep neural networks;

* in the training phase it exploits the interface information, but then in the inference phase it computes the fluxes without explicitly reconstructing the interface;

* it easily handles distributions where the interface is not regular, which is particularly important in presence of more than two materials;

* its implementation is easier than standard three-dimensional VOF methods and its efficiency relies on highly efficient deep learning libraries (e.g. Tensorflow [1], PyTorch [28] or JAX [4]).

The manuscript is organized as follows. In Section 2 we describe the basic concepts of the method, with particular focus on the neural network architecture, on the training dataset generation, on the network training and on its usage inside a FV solver. More advanced strategies, introduced to improve the efficiency and accuracy of the scheme and to enforce physical constraints, are discussed in Section 3. Numerical results are shown in Section 4 to empirically validate the method, whereas the final conclusions are presented in Section 5.

## 2. Method description.
As discussed in the Introduction, let us consider equation (1.1) on the space-time domain $\Omega \times [0, T]$. We consider a uniform Cartesian mesh $\mathcal{T}_h$ on the cubic domain $\Omega = [0, L]^3$ comprising $N_h$ cells in each direction. Generalization to Cartesian meshes with elongated cells or with a different number of cells in each direction is straightforward. Using a standard notation, we denote by $C_{i,j,k}$, $i, j, k = 1, \dots, N_h$, a generic cell of $\mathcal{T}_h$

$$C_{i,j,k} = \{(x, y, z) \in \Omega : ih < x < (i+1)h, jh < y < (j+1)h, kh < z < (k+1)h.\}$$

Moreover, we also introduce the stencil $S_{i,j,k}^m$ of margin $m \in \mathbb{N}$ and centered in $C_{i,j,k}$ as the indices of the $m^3$ cells around $C_{i,j,k}$

$$S_{i,j,k}^m = \{(i', j', k') : |i - i'| \le m, |j - j'| \le m, |k - k'| \le m\}.$$

The stencil has a central role in any VOF scheme, because the volume fractions of its cells are used to compute the flux needed to compute the volume fractions at the next time instant. A two-dimensional schematic representation of the VOF method is shown in Figure 1.

The basic idea of the VOF-ML method is to exploit the nonlinear approximation capabilities of neural networks to capture the complex relationship between the volume fractions in each cell of a given stencil $S_{i,j,k}^m$,
<figure><img src="image_3.png" /><figcaption>(1.a) Example of volume fractions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(1.b) Representation of the unknown flux.</figcaption></figure>

**Fig. 1** Example of volume fractions and flux for a given distribution. The fluid A occupies the red region. The velocity in the central cell is $u^*$, the time-step length is $\Delta t$.

the velocity of the fluid inside the cell $C_{i,j,k}$, and the flux exiting from the cell $C_{i,j,k}$. Such a relationship is the central part of a VOF method, since the volume fractions and the velocity fields are known, but the flux is unknown. The flux reconstruction is then used in a standard FV scheme to obtain the new volume fractions.

The VOF-ML method is characterized by three main components:

* the generation of the dataset, used to train a neural network;

* the training of the involved neural network;

* the use of the trained neural network inside a standard FV scheme.

Each of these components is described in detail in the next subsections.

**2.1. A generic neural network architecture.** We only consider fully-connected feed-forward neural networks. This is the simplest existing architecture, but it is very flexible and allows to accurately approximate very nonlinear and non-smooth operators.

Let $L$ be the number of layers, let $A_\ell \in \mathbb{R}^{a_\ell \times a_{\ell-1}}$ and $b_\ell \in \mathbb{R}^{a_\ell}$, $\ell = 1, ..., L$, be $L$ matrices and vectors of known dimension and let $\rho: \mathbb{R} \to \mathbb{R}$ be a fixed nonlinear function, which can be applied element-wise to define the operator $\rho: \mathbb{R}^{a_\ell} \to \mathbb{R}^{a_\ell}$, $\ell = 1, ..., L$. We define the parametric function $\alpha^{NN}: \mathbb{R}^{a_0} \to \mathbb{R}^{a_L}$ associated with the proposed neural network as

$$ (2.1) \qquad \begin{aligned} & x_0 = \mathbf{x}, \\ & x_\ell = \rho(A_\ell x_{\ell-1} + b_\ell), \quad \ell = 1, ..., L-1, \\ & \alpha^{NN}(\mathbf{x}) = A_L x_{L-1} + b_L. \end{aligned} $$

The input of the neural network is a vector containing the volume fractions in a stencil and the normalized velocity in a specific direction, whereas its output represents the approximation of the flux in the chosen direction. Therefore, the input dimension $a_0$ and the output dimension $a_L$ have to be $a_0 = m^3 + 1$ and $a_L = 1$.

Moreover, the nonlinear function $\rho$, also named activation function [11], can be arbitrarily chosen. Some common choices are $\rho(x) = \text{ReLU}(x) = \max(0, x)$, $\rho(x) = \text{RePU}(x) = \max\{0, x^p\}$ for $1 < p \in \mathbb{N}$, $\rho(x) = 1/(1 + e^{-x})$, $\rho(x) = \tanh(x)$ and $\rho(x) = \log(1 + e^x)$. In this work, we only consider the Rectified
Linear Unit ($\rho(x) = ReLU(x) = max(0, x)$) because many FV flux limiters can be exactly expressed by using such a function. For example, the operator

$$ \minmod(a, b) = \begin{cases} 0 & ab \le 0 \\ \min(|a|, |b|) \operatorname{sign}(a) & ab > 0 \end{cases} $$

can be expressed as $\minmod(a, b) = R(a - R(a - b)) - R(-a - R(b - a))$, when R is the ReLU function. The entries of all matrices and vectors in (2.1) are named trainable weights of the neural network. The training of the neural network is the optimization procedure through which such weights are optimized to minimize a suitable cost function, named loss function.

## 2.2. Dataset generation.

In this section, we describe the way in which we generate the training dataset. It consists of numerous input-output pairs, where each data pair is associated with a suitable distribution of the fluids A and B on a single stencil.

* The input is a vector of size $m^3 + 1$. It contains the $m^3$ volume fractions associated with the given fluids distribution and the local Courant number in the positive x-direction. The local Courant number $\beta$ is sampled in an interval $(C_{\min}, C_{\max}) \in [0, 1]$.

* The output is the numerical flux in the x-direction for the model equation (1.1).

Every input-output pair is obtained by purely geometrical considerations. In most papers combining neural networks and VOF schemes [35, 18, 6], the training dataset is generated using regular shapes or functions. Here, instead, we are interested in training a neural network that can accurately predict the flux in the presence of both smooth and non-smooth interfaces, because in many applications both cases are present in the same simulation, in different regions of the domain. A generic case with more fluids can be solved by coupling 2-fluids sub-problems as discussed in [2].

A generic 2D example is in Figure 1, where we subdivide the reference stencil in two regions A and B. However, in 3D the situation is more complicated as sketched in the next Figures. In this work, we consider the following synthetic distributions, where the region A is:

* either the intersection between a cube and one, two, or three parametric half-spaces,

* or the intersection between a cube and the interior of a parametric ellipsoid.

All technical formulas are provided in Appendix A.

In Figure 2 we show examples of a region cut by one plane, in Figure 3 we show examples of a region cut by two planes, in Figure 4 we show examples of a region cut by three planes, and, in Figure 5, we show examples of portions of the surface of the cutting ellipsoids, discretized by a finite number of points.

<figure><img src="image_5.png" /><figcaption>Fig. 2 Examples of regions defined by a single half-space.</figcaption></figure>

These regions are convex, but concave regions can be obtained by switching the roles of the two fluids. Note that this is a continuous set of geometrical configurations, in the sense that the intersection of two or three half-spaces and the interior part of an ellipsoid can tend to a single half-space, which is the simplest
<figure><img src="image_6.png" /><figcaption>Fig. 3 Examples of regions defined by two half-spaces.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4 Examples of regions defined by three half-spaces.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5 Examples of smooth regions cut by an ellipsoid generated as in Figure 6.</figcaption></figure>

considered case. Considering in our dataset the intersection between two or three half-spaces allows us to train
the neural network on less regular distributions. We hope it results in improving the prediction capabilities
for specific applications, but without changing the nature of the VOF-ML approach.

REMARK 2.1. *The results of our numerical tests show that this ensemble of configurations is sufficient to reach good prediction capabilities of our VOFML scheme in 3D.*

For each cell of the stencil, we have to compute the volume fraction associated with one of the two regions, say region A. Performing such a computation by numerical quadrature would be very expensive or inaccurate, due to the discontinuous nature of the target function. Instead, when considering one or more half-spaces, we explicitly construct the polytopes defined by the intersection between the region A and each cell of the stencil and we exactly calculate the volumes. In the case of an ellipsoid, we start by approximating it with a polytope with $N_{pts-el} = 10000$ vertices (generated as in Remark 2.2), then we construct the polytopes defined by the intersection of the approximated ellipsoid and each cell. This way, it is always possible to use known and efficient algorithms to evaluate the volume of a convex hull [3] and obtain the exact volume fractions in the first case, and a good approximation in the second one.
REMARK 2.2 (Discretization of the ellipsoid surface). In this work, all the integrals involved in the dataset construction are computed through the convex hull algorithm. It is therefore necessary to discretize each ellipsoid with a polytope, which has to be an accurate approximation of the ellipsoid to avoid introducing large errors in the training dataset.

We construct the polytope with $N_{pts-el}$ vertices using the Fibonacci spiral [19]. To define such a set of points $\{x_i, y_i, z_i\}_{i=0}^{N_{pts-el}-1}$, let us consider the corresponding spherical coordinates $\{\phi_i, \theta_i\}_{i=0}^{N_{pts-el}-1}$. Denoting by $G_R$ the golden ratio $G_R = (\sqrt{5}+1)/2$, $\phi_i$ and $\theta_i$ can be computed as:

$$ \phi_i = 2 \frac{\pi i}{G_R}, \quad \theta_i = \cos^{-1} \left( 1 - \frac{2i + 1}{n} \right), \quad i = 0, \dots, N_{pts-el} - 1. $$

Then, the Cartesian coordinates $x_i$, $y_i$ and $z_i$ can be retrieved as:

$$ [x_i, y_i, z_i] = [\cos(\phi) \sin(\theta), \sin(\phi) \sin(\theta), \cos(\theta)], \quad i = 0, \dots, N_{pts-el} - 1. $$

The obtained distribution is then mapped to the surface of the ellipsoid through a simple affine transformation. An example of the obtained distribution on the unit sphere, with $N_{pts-el} = 1500$, is shown in Figure 6. Note that the points have a quasi-uniform distribution on the entire surface.

<figure><img src="image_9.png" /><figcaption>Fig. 6 Fibonacci spiral with $N_{pts-el} = 1500$.</figcaption></figure>

Due to the high dimension of the corresponding parameter spaces (see Appendix A), we sample the training distributions using a Latin Hypercube Sampling technique in order to weaken the effect of the curse of dimensionality. At the end of this stage, one has designed the finite dataset of inputs

$$ (2.2) \qquad \mathcal{D}_{\text{inp}} = \{ (\mathbf{x}, \beta) \in \mathbb{R}^{m^3} \times \mathbb{R} \}, \quad \#\mathcal{D}_{\text{inp}} < \infty. $$

The exact numerical flux is represented as a geometrical function $\alpha$ calculated as the volume fractions, i.e. using the convex hull algorithm

$$ \alpha : \mathbb{R}^{m^3+1} \times \mathbb{R} \to \mathbb{R}. $$

The dataset of outputs is

$$ \mathcal{D}_{\text{out}} = \{ \gamma \in \mathbb{R} \mid \gamma = \alpha(\mathbf{x}, \beta) \text{ for } (\mathbf{x}, \beta) \in \mathcal{D}_{\text{inp}} \} \subset \mathbb{R}. $$

The total or training dataset is $\mathcal{D} = \{ (\mathbf{x}, \beta, \gamma) \in \mathbb{R}^{m^3+2} \mid (\mathbf{x}, \beta) \in \mathcal{D}_{\text{inp}}, \gamma = \alpha(\mathbf{x}, \beta) \}$.
**2.3. Training.** After the generation of the training dataset $\mathcal{D}$, a large set of $\#\mathcal{D}$ input-output pairs is available, and a supervised training procedure is performed. Let us consider the neural network described in Section 2.1 and fix the architecture free hyper-parameters $L, a_1, a_2, \dots, a_{L-1}$. The objective of training is to find a set of weights of the neural network such that $\alpha^{\mathcal{NN}}(\mathbf{x}, \beta) \approx \alpha(\mathbf{x}, \beta)$ for any possible pair $(\mathbf{x}, \beta) \in \mathcal{D}$.

To do so, we construct a loss function that penalizes large differences between $\alpha^{\mathcal{NN}}$ and $\alpha$. We consider the standard Mean Squared Error (MSE) loss function $\mathcal{L}$ defined as:

$$ (2.3) \qquad \mathcal{L} = \frac{1}{\#\mathcal{D}_{\text{inp}}} \sum_{(\mathbf{x}, \beta) \in \mathcal{D}_{\text{inp}}} (\alpha^{\mathcal{NN}}(\mathbf{x}, \beta) - \alpha(\mathbf{x}, \beta))^2. $$

Note that, given a neural network architecture, the loss function $\mathcal{L}$ depends only on the dataset $\mathcal{D}$ and on the neural network's trainable weights. In particular, the neural network and its training do not depend on the PDE that has to be solved (in this case equation (1.1)) or on any previously computed expensive simulation. All our training data are synthetic.

Because of the intrinsic nonlinearity of the neural network representation, the cost function $\mathcal{L}$ is highly non-convex, and it is not possible to find global minima of such a cost function. Therefore, we content ourselves with finding local minima through an iterative descent algorithm. Numerous algorithms are available, in this work we use a combination of the ADAM optimizer [20] and the BFGS optimizer [36]. At the end of the training, we have a function $\alpha^{\mathcal{NN}}$ that approximately represents the relation between the input vector and the target for all pair $(\mathbf{x}, \beta) \in \mathcal{D}$, that is

$$ (2.4) \qquad \alpha^{\mathcal{NN}}(\mathbf{x}, \beta) \approx \alpha(\mathbf{x}, \beta) \quad \text{for all } (\mathbf{x}, \beta) \in \mathcal{D}_{\text{inp}}. $$

**2.4. Definition of the 3D VOFML scheme.** Let us consider Problem (1.1) on the domain $\Omega$, discretized by the mesh $\mathcal{T}_h$ with mesh size $\Delta x$. Let $\Delta t$ be the time step size. The initial condition is evolved in time one time step at a time until the entire time interval $[0, T]$ is covered.

Since we adopt a FV scheme, the volume fraction $\alpha_{i,j,k}$ at a given time $t$ is updated at the time $t + \Delta t$ as $\overline{\alpha_{i,j,k}}$ using the generic formula

$$ (2.5) \qquad \overline{\alpha_{i,j,k}} = \alpha_{i,j,k} - \sum_{f \in \mathcal{F}(C_{i,j,k})} \tilde{g}(f, n) $$

where the sum is extended over the set of faces $\mathcal{F}(C_{i,j,k})$

$$ \mathcal{F}(C_{i,j,k}) = F_{i,j,k}^x \cup F_{i,j,k}^y \cup F_{i,j,k}^z, \quad F_{i,j,k}^x = \{f_{i \pm \frac{1}{2}, j, k}\}, \quad F_{i,j,k}^y = \{f_{i, j \pm \frac{1}{2}, k}\}, \quad F_{i,j,k}^z = \{f_{i, j, k \pm \frac{1}{2}}\}, $$

and $\tilde{g}(f, n)$ is the signed flux across the face $f$ in the direction normal $n$ in the time interval $(t, t + \Delta t)$, and $f_{i \pm \frac{1}{2}, j, k}$, $f_{i, j \pm \frac{1}{2}, k}$ or $f_{i, j, k \pm \frac{1}{2}}$ refer to the faces of the structured mesh. In practice, we use a directional splitting scheme to reduce the 3-dimensional problem (2.5) to a triplet of 1-dimensional problems. For the simplicity of the presentation only, we assume that the velocity field is constant with non-negative components in all directions

$$ u = (u_x, u_y, u_z) \quad \text{with } u_x, u_y, u_z \ge 0 $$

Then the splitting method is made of the following 3 steps

$$ (2.6) \qquad \begin{aligned} \alpha_{i,j,k}^x &= \alpha_{i,j,k} - u_x \frac{\Delta t}{\Delta x} \left( \alpha_{i+\frac{1}{2},j,k} - \alpha_{i-\frac{1}{2},j,k} \right), \\ \alpha_{i,j,k}^y &= \alpha_{i,j,k}^x - u_y \frac{\Delta t}{\Delta x} \left( \alpha_{i,j+\frac{1}{2},k} - \alpha_{i,j-\frac{1}{2},k} \right), \\ \overline{\alpha_{i,j,k}} &= \alpha_{i,j,k}^y - u_z \frac{\Delta t}{\Delta x} \left( \alpha_{i,j,k+\frac{1}{2}} - \alpha_{i,j,k-\frac{1}{2}} \right). \end{aligned} $$
For stability reasons, we assume a natural CFL constraint max(u<sub>x</sub>, u<sub>y</sub>, u<sub>z</sub>) Δt Δx ≤ 1.

The fluxes α<sub>i+1/2,j,k</sub> and α<sub>i-1/2,j,k</sub> are the fluxes going outside and inside on the two faces in F<sup>x</sup><sub>i,j,k</sub>. Similar reasonings hold for the other directions y and z. The iterations of equation (2.6) allows to simulate the entire flow from time t = 0 to the final time t = T.

DEFINITION 2.3 (Flux of the VOFML scheme). The originality of the VOFML scheme in 3D relies on the following choice of the numerical flux, which is hereafter given for the flux in the x direction for all i, j, k

$$
(2.7) \qquad \alpha_{i+\frac{1}{2},j,k} = \alpha^{NN}(\mathbf{x}, \beta_x)
$$

where

• the volume fractions are taken locally x = (α<sub>i',j',k'</sub>)(i',j',k')∈S<sup>m</sup><sub>i',j',k'</sub>,

• β<sub>x</sub> = u<sub>x</sub> Δt Δx ≥ 0 is the Courant number in the x direction,

• α<sup>NN</sup> is the function trained in (2.4) for the dataset D (2.2).

The same definition holds with evident modification in the y and z directions. If some components of the velocity are negative, the stencil S<sup>m</sup><sub>i,j,k</sub> is upwinded accordingly, and the Courant number is kept non-negative.

Note that (2.6) is the general update rule used in many VOF schemes, the only difference is in the way the fluxes are computed. In more classical VOF schemes such as [17, 13, 30, 35, 12, 18, 29, 6], the local stencil is used to reconstruct interface parameters, which are post-processed to evaluate the numerical flux. In this VOFML scheme the numerical flux is directly computed by evaluating a previously trained neural network. It saves one step by comparison with the general VOF procedure which results in an important algorithmic simplification.

**3. Physical constraints.** The VOFML scheme defined in the previous Section is a purely data-driven approach which does not take into account the natural symmetries attached to the model transport equation 1.1. This is a common problem for most data driven techniques used in combination with numerical solvers. It is therefore advisable to enforce additional physical constraints or symmetries in the network or in its training, as described in [14, 26, 9] in a different context. This enforcement can be performed a priori or a posteriori. What we call an a priori approach in the context of ML schemes is a ML architecture which naturally provides a quantity satisfying the symmetries. What we call an a posteriori approach is based on a post-processing. In this Section, we propose and analyze a few possible a priori strategies to increase the symmetry properties of the VOF-ML scheme. We pay particular attention to the symmetry group of a cubic structure in 3D. An a posteriori strategy to avoid creation of negative mass is also proposed in Section 3.2.

**3.1. 3D symmetries.** Several natural symmetries exist on Cartesian meshes. An important fact [23] is that a cube has a symmetry group of order 48 (the octahedral group), that is left invariant by 47 non-trivial symmetries/permutations generically denoted as σ (the trivial symmetry is the identity). The octohedral group is isomorphic to S<sub>4</sub> × Z<sub>2</sub>, that is once again 24 × 2 elements. It is clearly more symmetries than in 2D.

It explains that some care is needed to implement a symmetry-preserving method in our 3D neural network approach. In the following, we provide an arbitrary list of basic symmetries that are actually sub-groups of cubic symmetries. The notation "(1+c)" refers to the fact that the sub-group is made of 1 + c elements, where the first one is the identity and the other ones are non-trivial. We consider the following basic transformations, which are ordered with respect to one particular direction, arbitrarily chosen to the x-direction.

• T<sub>1</sub>: rotations leading to a reordering of the axes: the first axis, i.e. the one in the direction of the flux, can be the x-axis, the y-axis or the z-axis (1+2);

• T<sub>2</sub>: reflection with respect to the plane orthogonal to the first axis (1+1);

• T<sub>3</sub>: rotations of multiples of π/2 around the first axis (1+3);

• T<sub>4</sub>: reflection with respect to the plane orthogonal to the second axis (1+1).
<figure><img src="image_10.png" /><figcaption>Fig. 7 Comparison with the flux of Figure 1 (on a 2D configuration).</figcaption></figure>

All 48 symmetries of a cube can be obtained by combining these basic transformations. For example, a reflection around the third axis can be obtained by applying the transformations $T_3$ (with $\pi$ as the rotation angle) and $T_4$. Note that each of the considered transformations maps the region occupied by the stencil to itself but permutes the cells composing it. Therefore, since in a VOF scheme one is only interested in the volume fractions, we always work with the corresponding permutations of the volume fractions, and we do not apply the geometrical transformations to the cube associated with the stencil.

Let us now consider an arbitrary distribution of two fluids inside a cubic stencil, and let us focus on the flux in the x-direction for a constant velocity of the form $(v, 0, 0)$. Applying any transformation of type $T_1$ or $T_2$, then the flux may vary as illustrated on the left of Figure 7. On the other hand, the flux remains the same whenever any transformation of type $T_3$ or $T_4$ is applied, as illustrated on the left of Figure 7. For this reason, we distinguish two sub-groups of the symmetry group

$$S_{\neq} = \{\text{symmetries obtained by applying transformations } T_1 \text{ and } T_2\}, \quad \text{(exact flux may vary)},$$

$$S_{=} = \{\text{symmetries obtained by applying transformations } T_3 \text{ and } T_4\}, \quad \text{(exact flux is the same).}$$

Due to the different effects of these sets of transformations on the flux, we use the set $S_{=}^-$ to impose physical properties to the neural network (see Section 3.1.1) and the set $S_{=}^+$ to cheaply enrich the training dataset. In particular, note that to generate the triplet $(\mathbf{x}, \beta, \gamma) \in \mathcal{D}$ the computation of $m^3 + 1$ integrals is required (the $m^3$ volume fractions included in the vector $\mathbf{x}$ and the flux $\gamma = \alpha(\mathbf{x}, \beta)$). On the other hand, given the untransformed triplet $(\mathbf{x}, \beta, \gamma) \in \mathcal{D}$ and for any $\sigma \in S_{=}^+$, the new triplet $(\sigma\mathbf{x}, \beta, \gamma_{\sigma})$ can be more cheaply computed by permuting the entries of the vector $\mathbf{x}$ according to the transformation $\sigma$ and by computing just the integral $\gamma_{\sigma} = \alpha(\sigma\mathbf{x}, \beta)$. We thus consider the following enriched input dataset:

$$ (3.1) \qquad \mathcal{D}_{\text{inp}}^{\sigma} = \{ (\sigma \mathbf{x}, \beta) \in \mathbb{R}^{m^3} \times \mathbb{R}, \sigma \in S_{=} \}, \qquad \#\mathcal{D}_{\text{inp}}^{\sigma} < \infty, $$
and the corresponding total dataset

$$ \mathcal{D}^{\sigma} = \{ (\mathbf{x}, \beta, \gamma) \in \mathbb{R}^{m^3+2} \mid (\mathbf{x}, \beta) \in \mathcal{D}_{\text{inp}}^{\sigma}, \gamma = \alpha(\mathbf{x}, \beta) \}. $$

Note that $\mathcal{D}_{\text{inp}} \subset \mathcal{D}_{\text{inp}}^{\sigma}$ and $\mathcal{D} \subset \mathcal{D}^{\sigma}$ because the identity belongs to $S_{\neq}$. In the following, for the sake of simplicity, we drop the superscript $\sigma$ and we always refer to $\mathcal{D}_{\text{inp}}^{\sigma}$ and $\mathcal{D}^{\sigma}$ as $\mathcal{D}_{\text{inp}}$ and $\mathcal{D}$ because the method description is not influenced by the particular choice of the training dataset.

**3.1.1. Enforcing symmetry constraints.** Let us consider a set of volume fractions in a stencil and a given velocity $u^*$ in its central cell. As discussed in Section 3.1, if a transformation in $S_{=}$ is applied to the stencil, then the exact flux in the $x$-direction does not vary. As a consequence, the flux computed with any VOF scheme should not change with respect to the flux computed on the original stencil. However, even with a perfectly devised training, the neural network would not satisfy such an important physical constraint. Therefore, we need to add such a constraint in our neural network architecture to obtain networks that satisfy it by construction. We thus want to modify the flux computation such that, if a transformation in $S_{=}$ is applied to the stencil, the output of the neural network does not vary. There are 8 different permutations that preserve the flux

$$ S_{=} = \{\sigma_1, \sigma_2, \dots, \sigma_8\}. $$

Given a neural network $\alpha^{NN}$ as in (2.4), we define a new neural network $\hat{\alpha}^{NN}$ by

$$ (3.2) \qquad \hat{\alpha}^{NN}(\mathbf{x}, \beta) = \frac{1}{8} \sum_{j=1}^{8} \alpha^{NN}(\sigma_j \mathbf{x}, \beta) \quad \text{for all } (\mathbf{x}, \beta) \in \mathbb{R}^{m^3} \times \mathbb{R}. $$

This operation is performed by a linear combination, so it is an admissible neural network. The cost of evaluating $\hat{\alpha}^{NN}$ is a priori 8 times the cost of evaluating $\alpha^{NN}$, in practice less because some saving is possible with vectorization of the calculations.

**LEMMA 3.1.** *The new neural network preserves symmetries in $S_{=}$, that is*

$$ \hat{\alpha}^{NN}(\sigma_j \mathbf{x}, \beta) = \hat{\alpha}^{NN}(\mathbf{x}, \beta) \quad \text{for all } (\mathbf{x}, \beta) \in \mathbb{R}^{m^3} \times \mathbb{R} \text{ and all } \sigma_j \in S_{=}. $$

*Proof.* This is due to the group structure of $S_{=}$. Indeed

$$ \hat{\alpha}^{NN}(\sigma_j \mathbf{x}, \beta) = \frac{1}{8} \sum_{k=1}^{8} \alpha^{NN}(\sigma_k \sigma_j \mathbf{x}, \beta) = \frac{1}{8} \sum_{l=1}^{8} \alpha^{NN}(\sigma_l \mathbf{x}, \beta) $$

because, for any $j = 1, \dots, 8$, the index $l$ in $\sigma_l = \sigma_k \sigma_j$ covers $l = 1, \dots, 8$ as $k = 1, \dots, 8$. $\square$

The same post-processing technique can be used to enforce additional symmetries. For example, we can consider that the numbering of the fluid of interest leads to an additional physical symmetry. If the flux of the fluid A is $f_A$, then the flux of the fluid B is $f_B = 1 - f_A$. It can be represented with the operator $M : \mathbb{R}^{m^3} \to \mathbb{R}^{m^3}$ be such that

$$ M(\mathbf{x}) = (1 - x_1, \dots, 1 - x_{m^3}) \quad \text{for all } \mathbf{x} = (x_1, \dots, x_{m^3}) \in \mathbb{R}^{m^3}. $$

When the roles of the fluids are switched (i.e. when $M$ is applied to the input vector), the output $\alpha^{NN}(M\mathbf{x})$ of the network should coincide with $1 - \alpha^{NN}(\mathbf{x})$. Since it is not the case for a generic neural network such as (2.4), we define the new neural network $\bar{\alpha}^{NN}$

$$ (3.3) \qquad \bar{\alpha}^{NN}(\mathbf{x}, \beta) = \frac{1}{2} \alpha^{NN}(\mathbf{x}, \beta) + \frac{1}{2} (1 - \alpha^{NN}(M\mathbf{x}, \beta)) \quad \text{for all } (\mathbf{x}, \beta) \in \mathbb{R}^{m^3} \times \mathbb{R}. $$
LEMMA 3.2. The neural network $\bar{\alpha}^{NN}$ is additive in the sense

$$ \bar{\alpha}^{NN}(\mathbf{x}, \beta) + \bar{\alpha}^{NN}(M\mathbf{x}, \beta) = 1 \text{ for all } (\mathbf{x}, \beta) \in \mathbb{R}^{m^3} \times \mathbb{R}. $$

*Proof.* One has

$$ \begin{aligned} \bar{\alpha}^{NN}(M\mathbf{x}, \beta) &= \frac{1}{2}\alpha^{NN}(M\mathbf{x}, \beta) + \frac{1}{2}(1 - \alpha^{NN}(M(M\mathbf{x}, \beta))) \\ &= \frac{1}{2}\alpha^{NN}(M\mathbf{x}, \beta) + \frac{1}{2}(1 - \alpha^{NN}(\mathbf{x}, \beta)) \\ &= -\frac{1}{2}\alpha^{NN}(\mathbf{x}, \beta) + \frac{1}{2}(1 + \alpha^{NN}(M\mathbf{x}, \beta)) \\ &= 1 - \left[\frac{1}{2}\alpha^{NN}(\mathbf{x}, \beta) + \frac{1}{2}(1 - \alpha^{NN}(M\mathbf{x}, \beta))\right] = 1 - \bar{\alpha}^{NN}(\mathbf{x}, \beta), \end{aligned} $$

because $M \circ M(\mathbf{x}) = M(M\mathbf{x}) = \mathbf{x}$.

<figure><img src="image_11.png" /></figure>

To obtain a neural network satisfying both types of constraints, it is sufficient to combine the two operations by linear combination. It yields the new neural network

$$ (3.4) \qquad \tilde{\alpha}^{NN}(\mathbf{x}, \beta) = \frac{1}{2}\hat{\alpha}^{NN}(\mathbf{x}, \beta) + \frac{1}{2}(1 - \hat{\alpha}^{NN}(M\mathbf{x}, \beta)). $$

By construction, the neural network $\tilde{\alpha}^{NN}$ also satisfies the symmetry constraints of Lemma 3.1 and Lemma 3.2.

**3.2. Non-negativity constraints.** Since we are interested in using the neural network's predictions in (2.6) to update the volume fractions for the subsequent time step, we need to enforce the preservation of natural bounds. In our case, such bounds come from the fact that a volume fraction is necessarily between 0 and 1 since $\alpha$ is an indicator function

$$ 0 \le \alpha_C = \frac{1}{|C|} \int_C \alpha \, dx \, dy \, dz \le 1 $$

which shows up at three different places, where the first two items pose no problem. The real issue is the third one.

* For all $(\mathbf{x}, \beta) \in D_{\text{inp}}$, one has $\mathbf{x} \in [0, 1]^{m^3}$. Note that a CFL number is also between 0 and 1, so that $D_{\text{inp}} \subset [0, 1]^{m^3+1}$. This constraint is satisfied by construction of the dataset.

* Similarly, the exact flux is in bounds, that is $\alpha(\mathbf{x}, \beta) \in [0, 1]$ for all $(\mathbf{x}, \beta)$. This constraint is also satisfied by construction of the dataset. So one has by construction $D_{\text{out}} \subset [0, 1]$ and finally $D \subset [0, 1]^{m^3+2}$.

* More technical is the satisfaction of the maximum principle attached to the transport of the indicator function with the transport equation (1.1). The practical question is the following. Assume the initial discrete solution of the discrete scheme (2.6) is in bounds, that is $\alpha_{i,j,k} \in [0, 1]$ for all $(i, j, k)$. Then one must make sure that the final discrete solution is in bounds as well, that is $\bar{\alpha}_{i,j,k} \in [0, 1]$ for all $(i, j, k)$.

Our solution, inspired by [10, 9], is to analyze the last of these three items, starting from the consideration of the first step of the discrete scheme (2.6). Let us consider a cubic cell with edge length $\Delta x$ and with volume fraction $\alpha_{i,j,k} \in [0, 1]$. The volume of the cell is $\Delta x^3$, the volume occupied by the fluid A is $0 \le V_A = \alpha_{i,j,k}\Delta x^3$. The volume occupied by the complementary fluid B is $0 \le V_B = (1 - \alpha_{i,j,k})\Delta x^3$.
The first step of the splitting scheme (2.6) is

$$ (3.5) \qquad \alpha_{i,j,k}^x = \alpha_{i,j,k} - \beta_x \left( \alpha_{i+\frac{1}{2},j,k} - \alpha_{i-\frac{1}{2},j,k} \right), \qquad \beta_x = u_x \frac{\Delta t}{\Delta x} \ge 0. $$

For convenience, we use a complementary variable $\mu$ such that $\alpha + \mu = 1$ everywhere. That is $\alpha_{i,j,k} + \mu_{i,j,k} = 1$, $\alpha_{i \pm \frac{1}{2}, j,k} + \mu_{i \pm \frac{1}{2}, j,k} = 1$ and $\alpha_{i,j,k}^x + \mu_{i,j,k}^x = 1$. Since we use a flux which satisfies the additivity condition of Lemma 3.2, one can write

$$ (3.6) \qquad \mu_{i,j,k}^x = \mu_{i,j,k} - \beta_x \left( \mu_{i+\frac{1}{2},j,k} - \mu_{i-\frac{1}{2},j,k} \right), \qquad \beta_x = u_x \frac{\Delta t}{\Delta x} \ge 0. $$

Therefore the double inequality $0 \le \alpha_{i,j,k}^x \le 1$ is equivalent to $0 \le \alpha_{i,j,k}^x$ and $0 \le \mu_{i,j,k}^x$.

LEMMA 3.3. Assume the initial data is in bounds $0 \le \alpha_{i,j,k} \le 1$ for all $(i, j, k)$. Assume the fluxes satisfy $0 \le \alpha_{i,j,k}^x \le 1$ for all $(i, j, k)$. Assume the numerical fluxes satisfy

$$ (3.7) \qquad 0 \le \alpha_{i+\frac{1}{2},j,k} \le \frac{1}{\beta_x} \alpha_{i,j,k}^x \quad \text{and} \quad 0 \le 1 - \alpha_{i+\frac{1}{2},j,k} \le \frac{1}{\beta_x} (1 - \alpha_{i,j,k}^x) \quad \text{for all } (i,j,k). $$

Then the final data is in bounds $0 \le \overline{\alpha_{i,j,k}} \le 1$ for all $(i, j, k)$.

The second set of inequalities can be rewritten for the complementary variable $0 \le \mu_{i+\frac{1}{2},j,k} \le \frac{1}{\beta_x} \mu_{i,j,k}^x$.

*Proof.* Let us make the hypothesis that only the bounds from below are satisfied $0 \le \alpha_{i+\frac{1}{2},j,k}$ and $0 \le \mu_{i+\frac{1}{2},j,k}$. An easy interpretation is that non-negative volumes of fluids are exchanged between cells. Then we observe that (3.5-3.6) implies

$$ \left\{ \begin{array}{l} \alpha_{i,j,k}^x \ge \alpha_{i,j,k} - \beta_x \alpha_{i+\frac{1}{2},j,k}, \\ \mu_{i,j,k}^x \ge \mu_{i,j,k} - \beta_x \mu_{i+\frac{1}{2},j,k}. \end{array} \right. $$

One observes that $\beta_x \alpha_{i+\frac{1}{2},j,k}$ (resp. $\beta_x \mu_{i+\frac{1}{2},j,k}$) is some non-negative quantity subtracted from the budget in the cell. Of course, it is highly recommended to have non-negative cell fractions at the end of the process. Then it is natural to enforce $\alpha_{i,j,k} - \beta_x \alpha_{i+\frac{1}{2},j,k} \ge 0$ and $\mu_{i,j,k} - \beta_x \mu_{i+\frac{1}{2},j,k} \ge 0$ which are nothing than the upper bounds in the hypothesis (3.7) of the Lemma. One obtains $\alpha_{i,j,k}^x \ge 0$ and $\mu_{i,j,k}^x \ge 0$ which are equivalent to the claim. $\square$

REMARK 3.4. Using the method proposed in [10], all four inequalities in (3.7) can be rewritten as

$$ m_{i,j,k}^x \le \alpha_{i+\frac{1}{2},j,k} \le M_{i,j,k}^x $$

where

$$ m_{i,j,k}^x = \max \left( 0, 1 - \frac{1}{\beta_x} (1 - \alpha_{i,j,k}^x) \right) \quad \text{and} \quad M_{i,j,k}^x = \min \left( 1, \frac{1}{\beta_x} \alpha_{i,j,k}^x \right). $$

These inequalities are compatible under CFL $\beta_x \le 1$, because it is evident to check that the choice $\alpha_{i+\frac{1}{2},j,k} = \alpha_{i,j,k}$ satisfies all inequalities.

To enforce them, we adopt a simple projection strategy to ensure that the predicted flux $\alpha^{NN}(\mathbf{x})$ is always included in the admissibility interval. We predict the flux with the neural network and, if its output is outside the interval, we just select the extremum closer to the prediction. The formal formula that incorporates all issues previously discussed is in the next definition.
DEFINITION 3.5. The VOFML flux in 3D is given by the formula

$$ (3.8) \qquad \alpha_{i+\frac{1}{2},j,k}^{x,\text{final}} = P_{[m_{i,j,k}^x, M_{i,j,k}^x]} \tilde{\alpha}^{NN}(\mathbf{x}, \beta_x), \quad \mathbf{x} = (\alpha_{i',j',k'})(i',j',k') \in S_{i',j',k'}^m. $$

The same procedure is repeated in the directional splitting procedure for the directions y and z. If some components of the velocity are negative, the stencil $S_{i,j,k}^m$ is upwinded accordingly, and the Courant number is kept non-negative.

4. Numerical results. In this section, we show and discuss numerical experiments empirically proving the better performance of the VOF-ML scheme. Note that the implementation of the VOF-ML is done in the context of hybridization, see Appendix B. Hybridization simply means that the VOF-ML scheme is activated close to the interface and that another standard and much less costly scheme is activated away from the interface. It results in some reduction of the implementation cost.

For a fair comparison, we only compare the VOF-ML scheme with other schemes exploiting the directional splitting and that do not explicitly reconstruct the interface. In particular, we test the commonly used upwind scheme (UW) and the limited downwind scheme (LD). The former, despite its popularity, is known to be excessively diffusive for sharp interface advection problems, whereas the latter is very accurate and can even exactly reconstruct interfaces that are aligned with the grid.

The correctness of our implementation has been verified on a simple test (not reported here) with constant velocity $u = [1, 1, 1]^T$ and where the region A is the cube $[0.2, 0.8]^3$ inside a cubic domain $\Omega = (0, 1)^3$. In particular, we have checked that our implementation satisfies two properties of the upwind scheme: firstly the upwind scheme is exact for a Courant number equal to 1, and secondly the theoretical asymptotic convergence rate is 0.5 in the $L^1$ norm for a Courant number less than 1 [8]. This rate is the asymptotic rate of convergence. However, in our numerical tests we observe a pre-asymptotic convergence rate < 0.5 which is interpreted as the effect of the excessive numerical diffusion due to the specific velocity and initial conditions.

In [9, 2], the authors observed that the VOF-ML performances were significantly better than the upwind ones and slightly better than the limited downwind ones for two-dimensional problems. We now show that, for three-dimensional problems, where more complex interfaces may arise, the advantage of the VOF-ML scheme is much greater.

In all the numerical tests, we use the same neural network (without any retraining or fine-tuning). Using the notation of Section 2.1, the chosen architecture is characterized by the following hyperparameters:

$$ \rho(x) = ReLU(x), \quad L = 5, \quad a_0 = 3^3 + 1, \quad a_L = 1, \quad a_i = 50, \forall i = 1, 2, 3, 4. $$

Such a network contains 9151 trainable coefficients, which are optimized by minimizing (2.3) through 5000 epochs performed with the ADAM optimizer [20], followed by 5000 steps of the BFGS optimizer [36]. To be coherent with the fact that the distributions used to create the dataset have different parametrizations, we generate the training dataset using the set of synthetic distributions in Table 1, where a random Courant number

$$ \beta \in [0, 0.6] $$

is associated with each distribution.

As discussed in Section 2.2, to each of these distributions we associate the vector $\mathbf{x}$ of the corresponding volume fractions. Moreover, to enrich the dataset, for each of its element, we add all the vectors of the form $(\sigma\mathbf{x}, \beta, \gamma_\sigma)$, for all $\sigma \in S_{\neq}$ as discussed in Section 3.1. The dimension of the dataset is thus $(3000 + 6000 + 9000 + 6000) \cdot 6 = 144000$.

As usually done, the entire dataset is split into three smaller dataset: the training dataset (80% of the data) used to optimize the weights, the validation dataset (10% of the data) used to avoid overfitting, and
<table><caption>**Table 1** Distributions used to generate the training dataset $\mathcal{D}$.</caption><thead><tr><th>Number of samples</th><th>Generating geometrical objects</th><th>Number of associated parameters</th></tr></thead><tbody><tr><td>3000</td><td>1 plane</td><td>3</td></tr><tr><td>6000</td><td>2 planes</td><td>6</td></tr><tr><td>9000</td><td>3 planes</td><td>9</td></tr><tr><td>6000</td><td>1 ellipsoid</td><td>6</td></tr></tbody></table>

the test dataset (10% of the data) used to test the model. In particular, we compare the flux reconstruction capability of the three considered methods on the test dataset and show the results in Table 2. For each method, we compute the Mean Squared Error (MSE) and the Mean Absolute Error (MAE) defined as:

$$ \text{MSE} = \frac{1}{\#\mathcal{D}_{\text{test}}} \sum_{(\mathbf{x},\beta) \in \mathcal{D}_{\text{test}}} (\alpha^*(\mathbf{x},\beta) - \alpha(\mathbf{x},\beta))^2, \quad \text{MAE} = \frac{1}{\#\mathcal{D}_{\text{test}}} \sum_{(\mathbf{x},\beta) \in \mathcal{D}_{\text{test}}} |\alpha^*(\mathbf{x},\beta) - \alpha(\mathbf{x},\beta)|, $$

where $\mathcal{D}_{\text{test}}$ is the test dataset, $\#\mathcal{D}_{\text{test}}$ is its size, and $\alpha^*$ denotes the flux computed with one of the three methods. In Table 2, we highlight that the flux reconstruction is much more accurate, using both types of errors, when the VOFML method is used. In particular, when considering the MSE metric, VOFML is 43.7 times more accurate than UW and 13.2 times more accurate than LD, whereas when the MAE metric is involved, VOFML is 8.57 times more accurate than UW and 3.70 times more accurate than LD. Such an error is propagated at each time step during the simulation, it is thus crucial to reconstruct the flux in the most accurate way possible.

<table><caption>**Table 2** Flux reconstruction error on the test dataset with the different metrics and methods.</caption><thead><tr><th colspan="3">MSE</th><th colspan="3">MAE</th></tr><tr><th>UW</th><th>LD</th><th>VOFML</th><th>UW</th><th>LD</th><th>VOFML</th></tr></thead><tbody><tr><td>1.339e-02</td><td>4.048e-03</td><td>3.071e-04</td><td>3.995e-02</td><td>1.724e-02</td><td>4.660e-03</td></tr></tbody></table>

**4.1. Test 1: constant velocity.** Let us consider the domain $\Omega = (-1, 1)^3$ and the problem

$$ (4.1) \qquad \partial_t \alpha + \nabla \cdot (\alpha u) = 0, $$

with $u = [1, 2, 3]^T$, periodic boundary conditions on the entire boundary $\partial\Omega$ and initial condition $\alpha_0$. To mimic the Zalezak's initial condition in this three-dimensional domain, let us consider a sphere, from which we subtract a rectangular prism. Let us denote the two three-dimensional regions $R_1$ and $R_2$ as:

$$ R_1 = \{(x, y, z) \in \Omega : x^2 + y^2 + z^2 < 0.4^2\}, $$

$$ R_2 = \{(x, y, z) \in \Omega : |x| < 0.2, |y| < 0.2, z < 0\}. $$

Then, the region $R_{\text{zal}}$ is defined as $R_{\text{zal}} = R_1 \setminus R_2$. The initial condition $\alpha_0$ is the indicator function of the set $\hat{R}_{\text{zal}}$, obtained from $R_{\text{zal}}$ through a rotation $R_0$ around the center of the domain to avoid mesh-related artifacts during the simulation. The rotation $R_0$ can be obtained combining a anti-clockwise rotations of $\pi/5$, $\pi/7$ and $\pi/9$ radians around the $x$, $y$ and $z$-axis, respectively.
The domain is discretized with a uniform Cartesian grid, with $N_h$ cells in each direction and a total of $N_h^3$ cells. The considered values of $N_h$ are 10, 14, 20, 27, 38, 54, 75, and 105. The values are chosen to obtain almost equally spaced points in the plots of the error, which are in logarithmic scale. The initial volume fractions in each cell are computed with tensorial quadrature rules, higher order quadrature rules would be equally accurate because of the discontinuous nature of the initial condition.

We are interested in simulating the dynamic up to time $T = 2$. This way, the exact distribution $\alpha(\cdot, T)$ at time $t = T$ coincides with the initial distribution $\alpha(\cdot, 0)$ at time $t = 0$. Indeed, due to the periodic boundary conditions, the initial distribution $\alpha_0$ is advected through the entire domain exactly one, two, and three times in the $x$, $y$, and $z$ direction, respectively. This is useful to evaluate the accuracy of the method by computing the error with respect to the exact solution.

We solve problem (4.1) with the three proposed methods (upwind, limited downwind, and VOF-ML) on different meshes to understand the convergence with respect to the mesh size. In order to fix the CFL condition, in this case equal to 0.1, 0.2, and 0.3 in the three directions, we fix $\Delta t = 0.1\Delta x$. At the end of each simulation, we measure the accuracy by computing the three relative $L^1$ errors

$$E^{\text{UW}} = \frac{\|\alpha^{\text{UW}}(T) - \alpha_0\|_1}{\|\alpha_0\|_1}, \quad E^{\text{LD}} = \frac{\|\alpha^{\text{LD}}(T) - \alpha_0\|_1}{\|\alpha_0\|_1}, \quad E^{\text{VOF-ML}} = \frac{\|\alpha^{\text{VOF-ML}}(T) - \alpha_0\|_1}{\|\alpha_0\|_1}.$$

Here $\alpha^{\text{UW}}(T)$ and $\alpha^{\text{LD}}(T)$ are the solutions at time $t = T$, obtained with the upwind scheme and the limited downwind scheme, respectively. Instead, $\alpha^{\text{VOF-ML}}(T)$ is the solution, at time $t = T$, obtained with the hybrid method described in Appendix B. The fluxes close to the interface are computed with the neural network, whereas the other ones are computed with the limited downwind scheme.

The results are shown in Fig. 8. It is evident that the upwind scheme is not able to correctly simulate the dynamic. The excessive diffusion leads to a very large error, which decreases slowly when the mesh is refined. On the other hand, the limited downwind scheme and the VOF-ML perform very similarly when the mesh is coarse enough ($N_h \le 20$ for this specific test), but for finer meshes the VOF-ML method converges faster with respect to the number of cells. The empirical convergence rates with respect to the number of cells, highlighted inside the brackets in the legend of the figure, are 0.16, 0.58, and 1.01 for the upwind, limited downwind, and VOF-ML method, respectively.

In Figure 9, instead, we analyze the evolution of the ratio $R_{\text{mix}}$ between the number of mixed cells and the total number of cells in the mesh for three different refinements. Note that, in this first test problem, the initial condition is simply translated. Therefore, the size of the exact interface between the two fluids remains the same at any time instant $t \in [0, T]$, and $R_{\text{mix}}$ should remain approximately constant during the entire simulation. In the subfigures, we report results for the coarsest mesh ($N_h = 10$), an intermediate mesh ($N_h = 27$), and the finest mesh ($N_h = 105$). We can observe that the noise on the curves decreases when the mesh is refined.

Unsurprisingly, $R_{\text{mix}}$ grows when the upwind scheme is adopted. This happens because the upwind scheme is extremely diffusive and is not able to preserve the sharp interfaces present in the initial condition. On the other hand, $R_{\text{mix}}$ remains approximately constant with both the limited downwind scheme and the VOF-ML scheme, pretty close to the value associated with the initial condition. The fact that the lowest values of $R_{\text{mix}}$ are obtained with the limited downwind scheme, even if its overall performance is worse than the VOF-ML one, can be explained by the fact that the scheme is developed to exactly advect one-dimensional interfaces. The diffusion in each update of the directional splitting is thus minimized when the limited downwind scheme is adopted. On the other hand, no information about the diffusion is used to train the neural network, and the constant behaviour of the associated $R_{\text{mix}}$ is only a consequence of its high accuracy.

We also show, for the same meshes, the values of the ratio $R_{\text{mix}}^{T/0}$ between $R_{\text{mix}}$ at the end of the simulation and the initial condition. A value larger than 1 means that the interface between the two fluids expanded
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Errors for Test 1. The number in the brackets is the convergence rate and the slope of the corresponding regression dashed line.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9** Test 1. Ratio of mixed cells over the total number of cells. Evolution (in logarithmic scale) during the simulation for a coarse, intermediate, and fine mesh.</figcaption></figure>

because more mixed cells are present, whereas a value smaller than 1 means that the interface shrank, for the opposite reason. In Table 3, we can observe that $R_{\text{mix}}^{T/0}$ decreases only with the limited downwind scheme when the mesh is refined. However, with both the limited downwind scheme and the VOF-ML scheme, $R_{\text{mix}}^{T/0}$ is of order 1, whereas it increases dramatically with the upwind scheme, highlighting once more the excessive
diffusion of the method.

<table><caption>**Table 3** Test 1. Ratios $R_{\text{mix}}^{T/0}$ for different meshes.</caption><thead><tr><th>N<sub>h</sub></th><th>UW</th><th>LD</th><th>VOF-ML</th></tr></thead><tbody><tr><td>10</td><td>35.71</td><td>1.68</td><td>1.89</td></tr><tr><td>27</td><td>71.12</td><td>1.41</td><td>2.23</td></tr><tr><td>105</td><td>103.3</td><td>1.14</td><td>2.33</td></tr></tbody></table>

**4.2. Test 2: velocity fields** $u = [u_1, u_2, u_3]$ **with** $\frac{\partial u_1}{\partial x} = \frac{\partial u_2}{\partial y} = \frac{\partial u_3}{\partial z} = 0$. Let us consider a slightly more general test case. We consider equation (4.1) on the domain $\Omega = (0, 1)^3$, before the time limit $T = 1$, with periodic boundary conditions and with a velocity field $u = [u_1, u_2, u_3]$ such that:

$$ (4.2) \qquad \begin{aligned} u_1 &= 25 \sin(2\pi y)^2 \sin(2\pi z) y(y-1)z(z-1) \cos(\pi t/T), \\ u_2 &= 25 \sin(2\pi z)^2 \sin(2\pi x) z(z-1)x(x-1) \cos(\pi t/T), \\ u_3 &= 25 \sin(2\pi x)^2 \sin(2\pi y) x(x-1)y(y-1) \cos(\pi t/T). \end{aligned} $$

Note that $\frac{\partial u_1}{\partial x} = \frac{\partial u_2}{\partial y} = \frac{\partial u_3}{\partial z} = 0$. It is thus evident that $u$ is divergence-free, but the stricter property on the three derivatives ensures that, during the directional splitting, the intermediate updates are perfectly conservative (up to machine precision). This is sufficient to ensure global conservation, even if a directional splitting approach is adopted.

Inspired by an initial distribution considered in [2], let us also consider a more complex initial distribution. Let us define the following regions:

$$ \begin{aligned} R_{\text{sphere}} &= \{(x, y, z) \in \Omega : (x - 0.5)^2 + (y - 0.5)^2 + (z - 0.5)^2 < 0.2^2\}, \\ R_{3,x} &= \{(x, y, z) \in \Omega : |x - 0.5| < 0.3, |y - 0.5| < 0.075, |z - 0.5| < 0.075\}, \\ R_{3,y} &= \{(x, y, z) \in \Omega : |x - 0.5| < 0.075, |y - 0.5| < 0.3, |z - 0.5| < 0.075\}, \\ R_{3,z} &= \{(x, y, z) \in \Omega : |x - 0.5| < 0.075, |y - 0.5| < 0.075, |z - 0.5| < 0.3\}, \end{aligned} $$

We define $R_{\text{init}}$ as $R_{\text{init}} = R_{\text{sphere}} \cup R_{3,x} \cup R_{3,y} \cup R_{3,z}$. As in the previous case, the initial condition $\alpha_0$ is the indicator function of the set $\tilde{R}_{\text{init}}$ obtained from $R_{\text{init}}$ through the rotation $R_0$ defined before.

As in the previous case, we fix the time step size $\Delta t$ as $\Delta t = 0.1\Delta x$ and perform three simulations for each mesh, one with the upwind scheme, one with the limited downwind scheme, and one with the VOF-ML scheme. The results are reported in Figure 10. Once more, the upwind scheme is excessively diffusive and not competitive with the other schemes. The limited downwind scheme and the VOF-ML scheme reach similar accuracies on coarse meshes, but the latter is characterized by a faster convergence rate with respect to $N_h$ and rapidly becomes more accurate than the former when the mesh is refined.

In Figure 11, we show the evolution of $R_{\text{mix}}$ on this second test. Once more, it is evident that the upwind scheme is excessively diffusive, whereas the limited downwind scheme and the VOF-ML scheme are able to capture the sharp interfaces. Note that, in this case, the initial distribution is deformed by the velocity field (4.2) and therefore $R_{\text{mix}}$ should not remain constant. However, because of the factor $\cos(\pi t/T)$ in (4.2), the velocity field is symmetrical, in time, around the instant $t = T/2$. This implies that the initial distribution, deformed up to time $t = T/2$, is deformed back to the initial condition exactly in the same way, implying, as
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Errors for Test 2. The number in the brackets is the convergence rate and the slope of the corresponding regression dashed line.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 11** Test 2. Ratio of mixed cells over the total number of cells. Evolution (in logarithmic scale) during the simulation for a coarse, intermediate, and fine mesh.</figcaption></figure>

a consequence, a symmetrical behaviour in the expected evolution of $R_{\text{mix}}$. This symmetry is more evident in Figure 11.c because the fine mesh minimizes the additional noise that may ruin it.

The ratios $R_{\text{mix}}^{T/0}$ are shown in Table 4. Once more, only the upwind scheme's values grow excessively. It is interesting to remark that, in this case, $R_{\text{mix}}^{T/0}$ decreases below 1 for the upwind scheme, highlighting a slight
antidiffusivity property of the method. In particular, the values of the VOF-ML scheme and of the limited downwind scheme for the last mesh are almost equidistant from the target value 1.

<table><caption>**Table 4** Test 2. Ratios $R_{\text{mix}}^{T/0}$ for different meshes.</caption><thead><tr><th>N<sub>h</sub></th><th>UW</th><th>LD</th><th>VOF-ML</th></tr></thead><tbody><tr><td>10</td><td>7.61</td><td>1.29</td><td>1.34</td></tr><tr><td>27</td><td>10.14</td><td>0.97</td><td>1.30</td></tr><tr><td>105</td><td>20.45</td><td>0.81</td><td>1.26</td></tr></tbody></table>

**4.3. Test 3: general divergence-free velocity fields.** Lastly, inspired by a test in [13], we consider a velocity field $u = [u_1, u_2, u_3]$ where

$$ (4.3) \qquad \begin{aligned} u_1 &= 2 \sin(\pi x)^2 \sin(2\pi y) \sin(2\pi z) \cos(\pi t/T), \\ u_2 &= - \sin(\pi y)^2 \sin(2\pi x) \sin(2\pi z) \cos(\pi t/T), \\ u_3 &= - \sin(\pi z)^2 \sin(2\pi x) \sin(2\pi y) \cos(\pi t/T). \end{aligned} $$

Such velocity field is divergence-free, but the derivatives $\frac{\partial u_1}{\partial x}$, $\frac{\partial u_2}{\partial y}$ and $\frac{\partial u_3}{\partial z}$ are different than 0. This implies that the velocity used in the update in a single direction is not divergence-free, and therefore the mass conservation may be violated.

Let us denote by $\alpha_A$ and $\alpha_B$ the volume fractions of fluid A and B, respectively. Our goal is to ensure that $0 \le \alpha_A \le 1$, $0 \le \alpha_B \le 1$, $\alpha_A + \alpha_B = 1$ in each cell and at each time step. As discussed in Section 3.2, a simple flux projection is sufficient to avoid negative values of $\alpha_A$ and $\alpha_B$, but a more elaborate strategy is needed to ensure the other constraints.

Let us assume that, at a given instant, all volume fractions $\alpha^n$ are in the interval $[0, 1]$, and let us define $\alpha_A^n = \alpha^n$ and $\alpha_B^n = 1 - \alpha^n$. Because of the interchangeability of the fluids, both $\alpha_A$ and $\alpha_B$ satisfy the equation:

$$ \partial_n \alpha_{\uparrow} + \nabla \cdot (\alpha_{\uparrow} u) = 0, \quad \uparrow \in \{A, B\}, $$

that can be discretized as discussed in Section 2.4 via directional splitting. After a single directional update for both $\alpha_A$ and $\alpha_B$, we obtain the intermediate volume fractions $\tilde{\alpha}_A^{n+1}$ and $\tilde{\alpha}_B^{n+1}$ which may be greater than 1 and whose sum may be different than 1. We thus compute the normalized volume fractions $\alpha_A^{n+1}$ and $\alpha_B^{n+1}$ as

$$ \alpha_A^{n+1} = \frac{\tilde{\alpha}_A^{n+1}}{\tilde{\alpha}_A^{n+1} + \tilde{\alpha}_B^{n+1}}, \qquad \alpha_B^{n+1} = \frac{\tilde{\alpha}_B^{n+1}}{\tilde{\alpha}_A^{n+1} + \tilde{\alpha}_B^{n+1}}. $$

Note that $\alpha_A^{n+1}$ and $\alpha_B^{n+1}$ satisfy all the required constraints. It is thus possible to define $\alpha^{n+1}$ as $\alpha^{n+1} = \alpha_A^{n+1}$ to obtain an update that does not violate anymore the mass conservation constraints. Therefore, the described procedure is able to preserve the mass conservation of the initial volume fractions for a single directional update. Since the initial condition of our problem at time $t = 0$ does not violate the mass conservation constraints, the described procedure can be iteratively applied to advect the initial condition while preserving the mass of both fluids.

For efficiency reasons, we highlight that, after having computed a flux $\alpha_A^*$ to update $\alpha_A^n$, the flux $\alpha_B^* = 1 - \alpha_A^*$ is already available and does not require any further (possibly expensive) computations. For this reason, the computational cost of the procedure is almost the same as the one of the simple update used in the previous sections.
We can finally consider the problem

$$ \partial_t \alpha + \nabla \cdot (\alpha u) = 0, $$

with the velocity defined in (4.3), on the domain $\Omega = (0, 1)^3$, up to time $T = 2$, and with the initial condition as the indicator function of a sphere centered in $[0.35, 0.35, 0.35]$ and radius 0.15.

The convergence rates are shown in Figure 12. We can observe that, due to the more complex nature of the problem, a more evident pre-asymptotic phase is present for meshes with less than $N_h = 20$ cells in each direction. However, the error decays are similar to the ones shown in the previous cases, with the VOF-ML that significantly outperforms the limited downwind scheme and the upwind scheme that is excessively diffusive and converges very slowly.

<figure><img src="image_16.png" /><figcaption>**Fig. 12** Errors for Test 3. The number in the brackets is the convergence rate and the slope of the corresponding regression dashed line.</figcaption></figure>

Also in this case, we can analyze the behaviour of $R_{\text{mix}}$, plotted in Figure 13. As in the previous test problem, the upwind scheme is not able to accurately preserve the initial discontinuous interface, and the evolution of $R_{\text{mix}}$ is symmetrical when the limited downwind or the VOF-ML schemes are adopted. Because of the similarity between Figure 11 and Figure 13, we can infer that the additional normalization introduced in this last test case does not interfere with the previously observed properties of the methods.

In Table 5 we show the values of $R_{\text{mix}}^{T/0}$. As in the previous cases, only the limited downwind and the VOF-ML display values close to 1, with the VOF-ML scheme that is particularly close to such a target value.

**5. Conclusion.** In this paper, we have presented a machine learning based three-dimensional Volume Of Fluid (VOF) scheme on Cartesian grids, where a fully-connected feed-forward neural network is used to
<figure><img src="image_17.png" /><figcaption>Fig. 13 Test 3. Ratio of mixed cells over the total number of cells. Evolution (in logarithmic scale) during the simulation for a coarse, intermediate, and fine mesh.</figcaption></figure>

<table><caption>Table 5 Test 3. Ratios $R_{\text{mix}}^{T/0}$ for different meshes.</caption><thead><tr><th>N<sub>h</sub></th><th>UW</th><th>LD</th><th>VOF-ML</th></tr></thead><tbody><tr><td>10</td><td>20.95</td><td>4.10</td><td>2.35</td></tr><tr><td>27</td><td>20.75</td><td>0.87</td><td>1.09</td></tr><tr><td>105</td><td>41.13</td><td>0.75</td><td>1.07</td></tr></tbody></table>

approximate the unknown function mapping the volume fractions and the Courant number to the target flux used in the Finite Volume update. Once the network is trained, it is therefore possible to directly evaluate the flux, without the need to reconstruct the interface in advance, which is the current bottleneck of most VOF schemes. Moreover, the absence of the interface reconstruction step allows a greater flexibility because it allows the method to accurately handle less regular interfaces, commonly present in the presence of more than two materials, where classical linear interfaces are not enough to provide good approximations.

The implementation of the method is also easier than most three-dimensional VOF schemes, because the flux is directly computed by evaluating a trained neural network, which can be efficiently done thanks to any of the several existing machine-learning libraries. In order to train the network, we propose a standard fitting procedure on a synthetic dataset, which can be generated by using only geometrical information. This allows the user to train a single network and to use it to solve different equations, on different domains, and with different initial and boundary conditions without any retraining. Nevertheless, the training dataset can be adapted, for example, by considering more or less non-regular interfaces, if some additional knowledge is available. A great focus is also devoted to the enforcement of physical constraints, which is achieved by suitably modifying the neural network.

The method has been validated on advection problems with accurate results for initial condition problems. Future perspectives include the analysis of the scheme on more physically relevant problems and the comparison with more complex or machine-learning based VOF schemes. One of our next goals of research will be to couple the tridimensional VOF-ML scheme to a more general hyperbolic scheme for solving compressible Euler equations with internal free boundaries in 3D.

**Acknowledgements.** This work is funded by PEPR/IA (https://www.pepr-ia.fr/). Both authors thank Stéphane Jaouen for fruitful discussions during the elaboration stage of this work.
REFERENCES

[1] M. ABADI ET AL., *TensorFlow: Large-scale machine learning on heterogeneous systems*, 2015. Software available from `tensorflow.org`.

[2] M. ANCELLIN, B. DESPRÉS, AND S. JAOUEN, *Extension of generic two-component vof interface advection schemes to an arbitrary number of components*, Journal of Computational Physics, 473 (2023), p. 111721.

[3] C. B. BARBER, D. P. DOBKIN, AND H. HUHDANPAA, *The quickhull algorithm for convex hulls*, ACM Transactions on Mathematical Software (TOMS), 22 (1996), pp. 469–483.

[4] J. BRADBURY, R. FROSTIG, P. HAWKINS, M. J. JOHNSON, C. LEARY, D. MACLAURIN, G. NECULA, A. PASZKE, J. VANDERPLAS, S. WANDERMAN-MILNE, AND Q. ZHANG, *JAX: composable transformations of Python+NumPy programs*, http://github.com/google/jax, (2018).

[5] A. CABOUSSAT, M. FRANCOIS, R. GLOWINSKI, D. KOTHE, AND J. SICILIAN, *A numerical method for interface reconstruction of triple points within a volume tracking algorithm*, Mathematical and Computer Modelling, 48 (2008), pp. 1957–1971.

[6] A. CAHALY, F. EVRARD, AND O. DESJARDINS, *Plic-net: A machine learning approach for 3d interface reconstruction in volume of fluid methods*, International Journal of Multiphase Flow, 178 (2024), p. 104888.

[7] A. COHEN, O. MULA, AND A. SOMACAL, *High order recovery of geometric interfaces from cell-average data*, ESAIM: Mathematical Modelling and Numerical Analysis, 59 (2025), pp. 693–727.

[8] F. DELARUE AND F. LAGOUTIÈRE, *Probabilistic analysis of the upwind scheme for transport equations*, Archive for rational mechanics and analysis, 199 (2011), pp. 229–268.

[9] B. DESPRÉS AND H. JOURDREN, *Machine learning design of volume of fluid schemes for compressible flows*, Journal of Computational Physics, 408 (2020), p. 109275.

[10] B. DESPRÉS AND F. LAGOUTIÈRE, *Contact discontinuity capturing schemes for linear advection and compressible gas dynamics*, Journal of Scientific Computing, 16 (2001), pp. 479–524.

[11] S. R. DUBEY, S. K. SINGH, AND B. B. CHAUDHURI, *Activation functions in deep learning: A comprehensive survey and benchmark*, Neurocomputing, 503 (2022), pp. 92–108.

[12] V. DYADECHKO AND M. SHASHKOV, *Reconstruction of multi-material interfaces from moment data*, Journal of Computational Physics, 227 (2008), pp. 5361–5384.

[13] A. FAKHREDDINE, K. ALAMÉ, AND K. MAHESH, *Directionally-split volume-of-fluid technique for front propagation under curvature flow*, International Journal for Numerical Methods in Fluids, 96 (2024), pp. 1517–1554.

[14] K. GARANGER, J. KRAUS, AND J. J. RIMOLI, *Symmetry-enforcing neural networks with applications to constitutive modeling*, Extreme Mechanics Letters, 71 (2024), p. 102188.

[15] I. GOODFELLOW, Y. BENGIO, AND A. COURVILLE, *Deep Learning*, MIT Press, 2016. http://www.deeplearningbook.org.

[16] A. HAN, R. CHIODI, AND O. DESJARDINS, *Capturing thin structures in vof simulations with two-plane reconstruction*, Journal of Computational Physics, 519 (2024), p. 113453.

[17] C. W. HIRT AND B. D. NICHOLS, *Volume of fluid (VOF) method for the dynamics of free boundaries*, Journal of computational physics, 39 (1981), pp. 201–225.

[18] A. KAWANO, *A simple volume-of-fluid reconstruction method for three-dimensional two-phase flows*, Computers & Fluids, 134-135 (2016), pp. 130–145.

[19] B. KEINERT, M. INNMANN, M. SÄNGER, AND M. STAMMINGER, *Spherical fibonacci mapping*, ACM Transactions on Graphics (TOG), 34 (2015), pp. 1–7.

[20] D. P. KINGMA AND J. BA, *Adam: a method for stochastic optimization*, arXiv preprint arXiv:1412.6980, (2014).

[21] J. KROMER, J. POTYKA, K. SCHULTE, AND D. BOTHE, *Efficient sequential plic interface positioning for enhanced performance of the three-phase vof method*, Computers & Fluids, 266 (2023), p. 106051.

[22] R. J. LEVEQUE, *Finite volume methods for hyperbolic problems*, vol. 31, Cambridge university press, 2002.

[23] J. S. LOMONT, *Applications of finite groups*, Academic Press, 2014.

[24] A. MAK AND M. RAESSI, *A machine learning approach to volume tracking in multiphase flow simulations*, Fluids, 10 (2025).

[25] W. MENG, C. HAO YU, J. LI, Z. WU, AND R. AN, *Modeling of interfacial two-phase flows via vof-based ls method with weno scheme in the finite volume method*, Computers & Fluids, 291 (2025), p. 106580.

[26] B. K. MILLER, M. GEIGER, T. E. SMIDT, AND F. NOË, *Relevance of rotationally equivariant convolutions for predicting molecular properties*, ArXiv, abs/2008.08461 (2020).

[27] T. NAKANO, M. A. BUCCI, J.-M. GRATIEN, AND T. FANEY, *Machine learning model for gas-liquid interface reconstruction in cfd numerical simulations*, Fluids, 10 (2025).

[28] A. PASZKE ET AL., *Pytorch: An imperative style, high-performance deep learning library*, in Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035.

[29] A. PATHAK AND M. RAESSI, *A three-dimensional volume-of-fluid method for reconstructing and advecting three-material interfaces forming contact lines*, Journal of Computational Physics, 307 (2016), pp. 550–573.

[30] J. E. PILLIOD AND E. G. PUCKETT, *Second-order accurate volume-of-fluid algorithms for tracking material interfaces*, Journal of Computational Physics, 199 (2004), pp. 465–502.
[31] J. PITMAN AND N. ROSS, Archimedes, Gauss, and Stein, Notices AMS, 59 (2012), pp. 1416–1421.

[32] Y. QI, J. LU, R. SCARDOVELLI, S. ZALESKI, AND G. TRYGGVASON, Computing curvature for volume of fluid methods using machine learning, Journal of Computational Physics, 377 (2019), pp. 155–161.

[33] M. A. RUSHDI, S. YOSHIDA, C. HU, T. N. DIEF, A. E. SALEM, AND M. M. K. AND, Efficient and accurate machine learning models for volume of fluid (vof) simulation on uniform cartesian mesh, Engineering Applications of Computational Fluid Mechanics, 19 (2025), p. 2451774.

[34] M. TAVARES, C. JOSSERAND, A. LIMARE, J. M. LOPEZ-HERRERA, AND S. POPINET, A coupled vof/embedded boundary method to model two-phase flows on arbitrary solid surfaces, Computers & Fluids, 278 (2024), p. 106317.

[35] G. WEYMOUTH AND D. K.-P. YUE, Conservative volume-of-fluid method for free-surface simulations on cartesian-grids, Journal of Computational Physics, 229 (2010), pp. 2853–2865.

[36] S. WRIGHT, J. NOCEDAL, ET AL., Numerical Optimization, vol. 35, Springer, 1999.

[37] D. L. YOUNGS, Time-dependent multi-material flow with large fluid distortion, Numerical methods for fluid dynamics, (1982).

[38] J. F. Z. KRAUS, F. KÖHLER AND M. SCHÄFER, A deep learning approach for curvature prediction in algebraic volume of fluid methodsg, International Conference on Computational Methods for Coupled Problems in Science and Engineering, 9 (2021).

[39] C. ZHAO, J. MAAREK, S. M. TALEGHANI, AND S. ZALESKI, A hybrid continuum surface tension force for the three-phase vof method, Journal of Computational Physics, 504 (2024), p. 112872.

**Appendix A. Parametrization of the training configurations.** The considered training configurations are parametrized in the following way, in order for them to be described by parameters that vary in suitable hyper-cubes. This is important to efficiently sample random configurations in these high-dimensional parameter spaces. We always consider a normalized stencil, i.e. a stencil centered in the origin and comprising cubic cells of edge 1.

We start by introducing an important sampling strategy to uniformly sample points on the unit sphere, because it will be widely used in the following discussion.

REMARK A.1. A strategy to uniformly sample points $[a, b, c] \in \mathbb{R}^3$ on the unit sphere is to project a uniform distribution of points $[\alpha, \beta, \gamma] \in \mathbb{R}^3$ on a cylinder to the sphere. Such a simple method is based on Archimedes' theorem, see [31] for more details. Given $\theta_1 \in [0, 2\pi]$ and $\theta_2 \in [-1, 1]$, the distribution on a cylinder of height 2 and radius 1 is obtained as $a = \cos(\theta_1)$, $b = \sin(\theta_1)$ and $c = \theta_2$. It is then mapped to the unit sphere as $\alpha = \sqrt{1 - c^2} \cdot a$, $\beta = \sqrt{1 - c^2} \cdot b$ and $\gamma = c$. Only two parameters are needed to parameterize such a distribution.

We now describe the parametrization of the selected distributions.

• The region A is a half-space. The plane separating the region A and the region B is parametrized as:

$$p_1(x, y, z; \theta) = a_1(\theta)x + b_1(\theta)y + c_1(\theta)z - d_1(\theta),$$

where $\theta = [\theta_1, \theta_2, \theta_3] \in \mathbb{R}^3$ is a vector of free parameters and the region A is defined as the set of points $(x, y, z) \in \mathbb{R}^3$ such that $p_1(x, y, z; \theta) < 0$.

In order to obtain a uniform distribution of all the existing planes, we want $n_1(\theta) = [a_1(\theta), b_1(\theta), c_1(\theta)]$ to be uniformly distributed on the unit sphere (i.e. $\|n_1(\theta)\|_2 = 1$) and $d_1(\theta)$ to vary between in $[0, \frac{\|n_1(\theta)\|_1}{2}]$. This ensures that the region A and the region B coexist in the central cell of the stencil. The first two parameters $\theta_1$ and $\theta_2$ are involved in the sampling of all possible vectors $n_1(\theta)$ as in Remark A.1, whereas $d_1$ is computed as $d_1(\theta) = \theta_3 \frac{\|n_1(\theta)\|_1}{2}$, with $\theta_3 \in [0, 1]$. Therefore, the parameter $\theta$ must assume values in $[0, 2\pi] \times [-1, 1] \times [0, 1]$.

Examples of distributions that can be generated with this approach are represented in Figure 2.

• The region A is the intersection between two half-spaces. The interfaces of the two half-spaces can be defined as:

$$p_i(x, y, z; \theta) = a_i(\theta)x + b_i(\theta)y + c_i(\theta)z - d_i(\theta), \quad i = 1, 2,$$
with $\theta \in \mathbb{R}^6$. Note that it is always possible to sample the two planes independently with the strategy described in the previous case, but this would lead to two main issues. It might happen that the region A, if non-empty, is outside the stencil, or otherwise it might coincide with a half-space, at least in the central cell. To avoid these issues and ensure that the two planes intersect inside the central cell, we have to modify the sampling strategy. The main idea is to generate two parametric planes intersecting on the y-axis and to apply parametric rotations and translations to move them such that the origin is mapped into a random point inside the central cell.

Given the first parameter $\theta_1 \in [0, 2\pi]$, we define the two initial planes intersecting on the y-axis as

$$ \tilde{p}_1(x, y, z) = z, \quad \tilde{p}_2(x, y, z; \theta_1) = -\sin(\theta_1)x + \cos(\theta_1)z, $$

with normal vectors $\tilde{n}_1 = [0, 0, 1]^T$ and $\tilde{n}_2 = [-\sin(\theta_1), 0, \cos(\theta_1)]^T$.

Let $\theta_2, \theta_3, \theta_4 \in [0, 2\pi]$ be the angles of rotation around the three Cartesian axes. The three corresponding rotation matrices are:

$$ R_x(\theta_2) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\theta_2) & -\sin(\theta_2) \\ 0 & \sin(\theta_2) & \cos(\theta_2) \end{bmatrix}, \quad R_y(\theta_3) = \begin{bmatrix} \cos(\theta_3) & 0 & \sin(\theta_3) \\ 0 & 1 & 0 \\ -\sin(\theta_3) & 0 & \cos(\theta_3) \end{bmatrix}, $$

$$ R_z(\theta_4) = \begin{bmatrix} \cos(\theta_4) & -\sin(\theta_4) & 0 \\ \sin(\theta_4) & \cos(\theta_4) & 0 \\ 0 & 0 & 1 \end{bmatrix}, $$

whereas the general three-dimensional rotation matrix can be defined as $R = R_x(\theta_2)R_y(\theta_3)R_z(\theta_4)$. Such a matrix is used to rotate $\tilde{p}_1$ and $\tilde{p}_2$ by multiplying their normal vectors. We thus define $n_1$ and $n_2$

as $n_1 = R\tilde{n}_1$ and $n_2 = R\tilde{n}_2$. Note that the origin belongs to both rotated planes $\hat{p}_i(x, y, z) = n_i^T \begin{bmatrix} x \\ y \\ z \end{bmatrix}$, $i = 1, 2$. To obtain the sought planes $p_1(x, y, z; \theta)$ and $p_2(x, y, z; \theta)$, we apply a last translation map. Since translations in the direction of the intersection line between the two planes would not affect the planes position, the direction of the applied translation lies in the plane orthogonal to such direction. The admissible translation vectors are thus parametrized in polar coordinates in such plane, where $\theta_5 \in [0, 2\pi]$ represents the angle and $\theta_6 \in [0, 1]$ the scaling factor of the translation distance (i.e. when $\theta_6 = 0$ no translation is applied, when $\theta_6 = 1$ the translation is such that a single point of the line belonging to both planes lies in the central cell).

Therefore, the parameter $\theta \in \mathbb{R}^6$ describing the two final planes can vary in $[0, 2\pi]^5 \times [0, 1]$. Once the two planes $p_1(x, y, z; \theta)$ and $p_2(x, y, z; \theta)$ are known, the region A can be easily defined as the set of points $(x, y, z) \in \mathbb{R}^3$ such that $p_1(x, y, z; \theta) < 0$ and $p_2(x, y, z; \theta) < 0$.

Examples of distributions that can be generated with this approach are represented in Figure 3.

* The region A is the intersection of three half-spaces. The idea is similar to the previous case, we define three planes

$$ p_i(x, y, z; \theta) = a_i(\theta)x + b_i(\theta)y + c_i(\theta)z - d_i(\theta), \quad i = 1, 2, 3, $$

with $\theta \in \mathbb{R}^9$. The first 6 parameters are chosen as in the previous case for $p_1$ and $p_2$, whereas the remaining three are chosen as in the case with one plane, but adapting the definition of $d_3$ to ensure the intersection of the three planes exists and lies in the central cell. In particular, we change the sign of $d_3$ if the region $p_3(x, y, z; \theta) < 0$ does not intersects the region where $p_1(x, y, z; \theta) < 0$ and
$p_2(x, y, z; \theta) < 0$. Therefore, $\theta \in \mathbb{R}^0$ can vary in $[0, 2\pi]^5 \times [0, 1] \times [0, 2\pi] \times [-1, 1] \times [0, 1]$ and the region A is the set of points $(x, y, z) \in \mathbb{R}^3$ such that $p_1(x, y, z; \theta) < 0$, $p_2(x, y, z; \theta) < 0$ and $p_3(x, y, z; \theta) < 0$. Examples of distributions that can be generated with this approach are represented in Figure 4.

* The region A is the interior part of an ellipsoid. Let us consider the following parametric ellipsoid:

$$s(x, y, z; \theta) = \frac{(x - x_0(\theta))^2}{a(\theta)^2} + \frac{(y - y_0(\theta))^2}{b(\theta)^2} + \frac{(z - z_0(\theta))^2}{c(\theta)^2} - 1,$$

where $\theta \in \mathbb{R}^6$. Once more, a naive parametrization may lead to ellipsoids outside the central cell, we thus propose the following parametrization.

The center $[x_0(\theta), y_0(\theta), z_0(\theta)]$ is chosen as

$$[x_0(\theta), y_0(\theta), z_0(\theta)] = [\theta_1, \theta_2, \theta_3] \in [-\sqrt{3(m + 0.5)}, \sqrt{3(m + 0.5)}]^3$$

where $[-m - 0.5, m + 0.5]^3$ is the region covered by the stencil (we remind that $m$ is the stencil margin).

Let us now consider the values $a(\theta)$, $b(\theta)$ and $c(\theta)$. Three temporary values $[\tilde{a}(\theta), \tilde{b}(\theta), \tilde{c}(\theta)]$ are selected on the unit sphere as in Remark A.1 sampling the parameters $\theta_4 \in [0, 2\pi]$ and $\theta_5 \in [-1, 1]$. Moreover, the last parameter $\theta_6$ can vary in $[0, 1]$ and is used as a nonlinear scaling to compute $[a(\theta), b(\theta), c(\theta)]$ from the values $[\tilde{a}(\theta), \tilde{b}(\theta), \tilde{c}(\theta)]$. The value $\theta_6 = 0$ is associated with a scaling such that a single point (or more generally, a finite number of points) of the ellipsoid surface intersects the border of the central cell. The value $\theta_6 = 1$ is associated with a scaling such that the entire central cell is inside the ellipsoid and the ellipsoid volume is minimal.

Therefore, $\theta \in \mathbb{R}^6$ can vary in $[-\sqrt{3(m + 0.5)}, \sqrt{3(m + 0.5)}]^3 \times [0, 2\pi] \times [-1, 1] \times [0, 1]$ and the region A is defined as the set of points $(x, y, z) \in \mathbb{R}^3$ such that $s(x, y, z; \theta) < 0$.

Examples of distributions that can be generated with this approach are represented in Figure 5. Here, the surface of the ellipsoid is represented by the vertices of the polyhedron used to approximately compute its volume. Note that, to obtain a similar integration precision for any ellipsoid, more points are needed for surfaces with larger curvature. This is automatically handled by the fact that we approximate a sphere using the points of the Fibonacci spiral [19], and then we map the surface of the sphere to the surface of the ellipsoid through an affine transformation (see Section 2.2).

## Appendix B. Hybridization of the VOF flux.

One of the advantages of neural networks is that, after an usually expensive training phase, evaluating the trained neural network is extremely efficient. It is thus possible to train a model on a large dataset just once, and then use the trained network multiple times, without the need to retrain. However, evaluating a neural network is still more expensive than some simple and less accurate VOF schemes. Therefore, we propose a hybrid strategy in which the VOF-ML scheme is used close to the interface, where high accuracy is required, and a less expensive scheme is adopted far from it, where the computation of the flux is trivial.

Observe that, if a volume fraction is exactly equal to 1, it means that inside the corresponding cell there is only fluid A. Therefore, the entirety of the flux going out from such a cell will be composed of fluid A, i.e. the flux will be 1 too. Analogously, if a volume fraction is exactly equal to 0, also the flux fraction going out from the cell will be 0. We thus fix a small threshold $\epsilon_{\text{mark}} > 0$ to mark the cells that are far enough from these two trivial cases. In Section 4, we set $\epsilon_{\text{mark}} = 0.01$, but this value may be chosen according to the used mesh size.

At the beginning of each iteration, we mark all the cells where we want to compute the flux going out from the cell with the VOF-ML scheme. All the remaining fluxes are computed with the Limited Downwind
scheme, which is cheaper but less precise. For each time step, the marking rule is

$$ \text{if } \varepsilon_{\text{mark}} \le \alpha_{i,j,k}^n \le 1 - \varepsilon_{\text{mark}} \Rightarrow \text{the cell } C_{i,j,k} \text{ is marked.} $$

Such a marking procedure can be applied to each intermediate step of (2.6) to further enhance the efficiency of the method by changing the adopted flux computation scheme more often. The marked cells are named *mixed cells*, whereas the unmarked cells are named *pure cells*.