=Paper= {{Paper |id=Vol-2964/article_187 |storemode=property |title=A Deep Learning Algorithm for Piecewise Linear Interface Construction (PLIC) |pdfUrl=https://ceur-ws.org/Vol-2964/article_187.pdf |volume=Vol-2964 |authors=Mohammadmehdi Ataei,Erfan Pirmorad,Franco Costa,Sejin Han,Chul Park,Markus Bussmann |dblpUrl=https://dblp.org/rec/conf/aaaiss/AtaeiPCHPB21 }} ==A Deep Learning Algorithm for Piecewise Linear Interface Construction (PLIC)== https://ceur-ws.org/Vol-2964/article_187.pdf

A Deep Learning Algorithm for Piecewise Linear Interface Construction (PLIC)
Mohammadmehdi Ataei,1,2 Erfan Pirmorad,2 Franco Costa,3 Sejin Han,4 Chul B Park,2 Markus
Bussmann2
1
Vector Institute for Artificial Intelligence, 661 University Ave Suite 710, Toronto, ON M5G 1M1, Canada
2
University of Toronto, 5 King’s College Rd, Toronto, ON M5S 3G8, Canada
3
Autodesk, Inc., 259-261 Colchester Rd., Kilsyth, VIC. 3137, Australia
4
Autodesk, Inc. 2353 North Triphammer Rd., Ithaca, NY 14850, USA
{ataei, pirmorad}@mie.utoronto.ca, {franco.costa, sejin.han}@autodesk.com, {park, bussmann}@mie.utoronto.ca

Abstract
Piecewise Linear Interface Construction (PLIC) is frequently
used to geometrically reconstruct fluid interfaces in Compu-
tational Fluid Dynamics (CFD) modeling of two-phase flows.
PLIC reconstructs interfaces from a scalar field that repre-
sents the volume fraction of each phase in each computational
cell. Given the volume fraction and interface normal, the lo-
cation of a linear interface is uniquely defined. For a cubic
computational cell (3D), the position of the planar interface
is determined by intersecting the cube with a plane, such that
the volume of the resulting truncated polyhedron cell is equal
to the volume fraction. Yet it is geometrically complex to find
the exact position of the plane, and it involves calculations
that can be a computational bottleneck of many CFD models.
However, while the forward problem of 3D PLIC is challeng-
ing, the inverse problem, of finding the volume of the trun-
cated polyhedron cell given a defined plane, is simple. In this
work, we propose a deep learning model for the solution to
the forward problem of PLIC by only making use of its in-
verse problem. The proposed model is up to several orders
of magnitude faster than traditional schemes, which signifi-
cantly reduces the computational bottleneck of PLIC in CFD
simulations.

Introduction
Multiphase flows (flows consisting of more than one com-
ponent) are ubiquitous in nature and are of interest in a Figure 1: a) An example of a VOF scalar field for a liquid-
variety of practical applications, and as such a significant gas system. α = 1 in liquid cells, α = 0 in gas cells, and
amount of work has been conducted on the computational 0 < α < 1 in interface cells. b) Interface reconstruction
modeling of such flows (e.g. Balachandar and Eaton 2010; using PLIC.
Rothman and Zaleski 1994; Wang et al. 2016). One of the
challenges of multiphase flow modeling is to track interfaces
between different phases. Among numerous Eulerian and and 0 < α < 1 in interface cells. An advection equation
Lagrangian interface tracking schemes, the Volume-of-Fluid is solved for α to track the location of interface through-
(VOF) Eulerian method is widely used (Hirt and Nichols out the simulation. There are two major categories of VOF
1981). schemes: algebraic and geometric (Lafaurie et al. 1994). Ge-
In the VOF method, as shown in figure 1a, the volume ometric VOF schemes provide more accurate fluid-fluid in-
fraction of one the two phases within each computational terface reconstruction compared to algebraic schemes, but
cell is represented by a scalar field α. For example, for a require more complicated geometrical calculations. Of these
liquid-gas system, α = 1 or α = 0 in each computational methods, Piecewise Linear Interface Construction (PLIC) is
cell that is occupied by the liquid or gas phase, respectively, frequently used to geometrically reconstruct an interface as
Copyright © 2021 for this paper by its authors, Use permitted under a line (2D) or plane (3D) (Youngs 1982). As shown in figure
Creative Commons License Attribution 4.0 International (CC BY 1b, for a cubic cell, PLIC linearly approximates an interface,
4.0). such that the volume of the resulting truncated polyhedron
Figure 2: All possible cases for intersection of a plane and a unit cube.

to solve PLIC in 2D for a square cell. However, analyti-
cal solutions in 3D for a cubic cell include multiple if-else
conditions, and require more sophisticated operations (e.g.
dealing with the inverse of third degree polynomials and
complex numbers) (Scardovelli and Zaleski 2000; Yang and
James 2006). For example, Scardovelli and Zaleski (Scar-
dovelli and Zaleski 2000) found an analytical solution to the
forward problem by inverting the inverse problem of PLIC
(i.e. finding the volume fraction of a polyhedron formed
given a plane with ~n and d), which is a simpler geomet-
ric calculation. This analytical solution includes a number
of conditional statements to select a solution from a set of
equations, corresponding to different cases of the intersec-
tion of a plane with a cube. This selection can slow down
a CFD solver, and is not easily parallelizable. To reduce
the computational cost, alternative iterative and approxima-
tion techniques have been developed to calculate d (Skarysz,
Garmory, and Dianat 2018; Rider and Kothe 1998; Kawano
2016).
Machine Learning (ML) has attracted considerable atten-
tion in the computational sciences, including CFD. Over
Figure 3: DNN architecture. the last decade, ML techniques have been used extensively
in fluid mechanics to solve a variety of problems, includ-
ing model discovery, modeling of reduced order, process-
ing of experimental data, shape optimization, and flow con-
cell is equal to the volume fraction α, where the plane nor- trol (Brunton, Noack, and Koumoutsakos 2020; Raissi, Yaz-
mal unit vector ~n is calculated as the gradient of the volume dani, and Karniadakis 2020; Brenner, Eldredge, and Freund
∇α
fraction scalar field, as ~n = − ||∇α|| (Kothe and Mjolsness 2019; Singh, Medida, and Duraisamy 2017; Pfaff et al. 2020;
1992). Eliminating symmetrical cases, the number of possi- Nikolaev, Richter, and Sadowski 2020; D’Elia et al. 2020;
ble cube-plane intersections can be reduced to the five sce- Yoon, Melander, and Verzi 2020). Previously, we have pub-
narios shown in figure 2, where α ≤ 0.5 (for α > 0.5, α is lished an alternative neural network based technique to per-
changed to 1 − α to match one of the five cases). The con- form PLIC calculations, by training a neural network with a
stant of the plane d is the smallest distance from the plane large synthetic dataset of PLIC solutions (Ataei et al. 2021a).
to the origin of the Cartesian coordinate that conforms with In this work, we propose a direct solution to the forward
the cubic cell edges, which must be calculated. The problem problem of PLIC using a fully connected Deep Neural Net-
of computing d is known as the forward problem of PLIC work (DNN). Unlike the previous work, we demonstrate that
(Scardovelli and Zaleski 2000). a DNN can solve the forward problem given only the inverse
To date, many studies have used PLIC for interface re- problem of PLIC, without requiring a synthetic dataset. We
construction on structured and unstructured grids, and also show that, in addition to its simplicity, this approach signifi-
for related calculations such as calculating the distance be- cantly outperforms other traditional PLIC implementations.
tween two interfaces, and have proposed improvements to
reduce the implementation complexity and increase its effi- Methodology
ciency (López et al. 2005; Ataei et al. 2021b; Dai and Tong
2019; Ataei et al. 2021a; Nahed and Dgheim 2020; Scheufler We solve the forward problem of PLIC by converting it into
and Roenby 2019). A simple analytical solution can be used an optimization problem that can be solved using a DNN. In
the forward problem, we must find d such that the volume This DNN is implemented in the Pytorch deep learning
V of a polyhedron cell is equal to the volume fraction α in a library (Paszke et al. 2019). We use ReLU activation func-
unit cubic cell: tions for all layers except for the output layer, which is a
linear activation function. The gradient of the custom loss
V =α (1) function in Eq. 2 is calculated using the Pytorch automatic
This equation can be converted into an optimization problem differentiation autograd package.
for a DNN by defining a loss function L as: For the inputs of the DNN m ~ is uniformly sampled where
all three coordinates are positive. Similarly, α is sampled
uniformly between 0 and 0.5. We used 5000 randomly shuf-
X
L= ~ − α)2
(V (d, m) (2)
batch
fled sample points in total, split into three datasets: 70%
training, 20% test, and 10% validation. The model was
where m~ = (m1 , m2 , m3 ) are the components of ~n along trained on the training dataset, the validation dataset was
the axes shown in figure 2 (Jafari, Shirani, and Ashgriz used to ensure that the model was not overfitting during the
2007), calculated as: training, and the test dataset was used to evaluate the final
m1 := min (|nx | , |ny | , |nz |) ≥ 0 accuracy of the model.
m3 := max (|nx | , |ny | , |nz |) > 0 (3) The network weights and biases were randomly initial-
m2 := |nx | + |ny | + |nz | − m1 − m3 ≥ 0 ized at the beginning using Xavier initialization (Glorot and
Bengio 2010). The results presented here are the average of
As shown in figure 3, we design a multilayer perceptron six runs with six different random seeds for network initial-
(MLP) neural network with fully connected layers, com- ization. The weights and biases were tuned using the Adam
prised of an input layer, H hidden layers each with N neu- Optimizer (Kingma and Ba 2015) with a batch size of 64 for
rons, and an output layer, to solve this optimization problem. 512 epochs.
The network outputs d based on m ~ and α as inputs. If the dif-
ference between V (d, m),~ calculated based on the predicted
d (the output of the DNN), and the known α is sufficiently
Results
minimized, then the forward PLIC problem is solved by the Evaluation metrics
DNN. One of the advantages of this technique is that (un- The training was performed on a PC running Ubuntu 20.04
like conventional methods for training a DNN), there is no (Intel i7-8700K, NVIDIA 1080 Ti graphics card (CUDA
need for synthetic data to train the DNN, as the DNN auto- version 11.0), 32GBs of RAM). It takes about 20 minutes
matically finds the solution to the forward PLIC problem by to train each network.
minimizing this loss function. The training progress of a DNN with H = 1, N = 24
Fortunately, formulating V (i.e. the inverse problem of is shown in figure 4. The Mean Squared Error (MSE) of
PLIC) in terms of the plane parameters (d, m) ~ is a straight- both training and validation datasets decreased during train-
forward geometrical problem. The volume under the plane ing. The model performance was evaluated against the test
can be found by calculating the volume of the pyramid dataset. Figure 5 shows the MSE and the maximum abso-
formed between the intersection of the plane with the three lute error for DNNs with different H and N . Notice the net-
coordinates (S1 , S2 , S3 ), and then subtracting the volume of work with H = 2, N = 48 is the most accurate; however,
the parts that extend beyond the cube (see figure 2). Such other networks with fewer neurons and hidden layers also
calculations result in a set of equations for V for each of the perform on-par or better when compared to PLIC models
five cases, as follows (Scardovelli and Zaleski 2000): that approximate PLIC (e.g. Kawano 2016).
To visualize the accuracy of the DNN model, figure 6
0 ≤ d ≤ m1 plots the value of α based on the value of d predicted by

 f1 (d, m)
~
1
 f (d, m)

2 ~ m ≤d≤m 1 2 the DNN (H = 2, N = 48) against the real value of α,
V = · f3 (d, m)
~ m2 ≤ d ≤ min (m1 + m2 , m3 ) (4) for 200 random data points from the test dataset. It can be
λ   f (d, m)
~ m3 ≤ d
 4 seen that there is excellent agreement between the predicted
f5 (d, m)
~ min (m1 + m2 , m3 ) ≤ d ≤ m3 α and real α, close to the ideal where all points would lie on
where: the diagonal line.

f1 = d3 Speedup
f2 = d3 − δ13 The advantages of using a DNN for PLIC calculations is that
f3 = d3 − δ13 − δ23 (5) it performs much faster than an analytical solution. More-
f4 = d3 − δ13 − δ23 − δ33 over, the DNN model can be easily parallelized on both
= 6m1 m2 d − 12 (m1 + m2 )

f5 CPUs and GPUs.
where: Figure 7 shows a comparison of the DNN (H = 2, N =
48) versus the Scardovelli and Zaleski (Scardovelli and Za-
δ1 = d − m1 leski 2000) analytical solution. To compare the models, the
δ2 = d − m2 analytical solution was implemented in C++, and compiled
(6) using optimization flag -O3 using the GCC compiler ver-
δ3 = d − m3
λ = 6 m1 m2 m3 sion 9.3. The Pytorch DNN model was converted into Torch-
Figure 4: Training progress over 512 epochs with batch size
of 64 for a DNN with one hidden layer with 24 neurons with
seed = 100.

Figure 6: Plots of the DNN predictions vs real α by a DNN
with two hidden layers each with 48 neurons (seed = 400).
To avoid clutter, only 200 random points are selected from
the test dataset.

that, while traditional algorithms such as Scardovelli and
Zaleski (Scardovelli and Zaleski 2000) require fewer float-
ing operations than a neural network approach, the neural
network model better utilizes the hardware and operates at
higher floating operations per second (Ataei et al. 2021a).

Figure 5: Mean Squared Error (MSE) and maximum abso-
lute error on the test dataset, for different DNN models with
different H and N .

script and imported into C++ as well. The std::chrono
library and Pytorch’s Profiler were used to measure the ex-
ecution time of the Scardovelli and Zaleski PLIC solution
and the neural model, respectively.
First, we compare the execution time of both models on
one CPU core. Since the DNN can perform PLIC calcula-
tions in batches, we ran the models for different numbers
of cells with randomly and uniformly generated α and m. ~ Figure 7: Speedups of the DNN model compared to the an-
On one CPU, the DNN model is 3.4 times faster than the alytical solution of Scardovelli and Zaleski (running on one
analytical model for 1000 cells. However, as the number of CPU core as the baseline).
cells increases to 10 million, the speedup gain is only 1.3.
This is expected, as the overhead associated with moving the
tensors to RAM becomes substantial for large batch sizes. Implementation in a CFD Solver
The DNN model accelerates significantly better on a GPU LBfoam (Ataei et al. 2021b) is a free surface lattice Boltz-
(even taking into account the overhead of moving data to mann solver for foaming simulations. In LBfoam, a disjoin-
the GPU), and with 10 million cells we achieve more than ing pressure Π between two adjacent bubble interfaces is
41 times speedup. In our previous work, we demonstrated responsible for stabilizing liquid lamellae (the thin layer of
liquid separating two bubbles). The disjoining pressure is ac- Dai, D.; and Tong, A. Y. 2019. Analytical interface recon-
tive up to a distance δmax , and is a function of the distance struction algorithms in the PLIC-VOF method for 3D poly-
between bubble interfaces δ as: hedral unstructured meshes. International Journal for Nu-
merical Methods in Fluids 91(5): 213–227.
0 δ > δmax D’Elia, M.; Karniadakis, G.; Pang, G.; and Parks, M. 2020.
Π= δ (7)
kΠ 1 − δmax δ < δmax Nonlocal physics-informed neural networks - A unified the-
where kΠ is a constant. In order to calculate δ, LBfoam re- oretical and computational framework for nonlocal models.
constructs adjacent bubble interfaces with PLIC. In AAAI Spring Symposium: MLPS.
To verify that the DNN model works as expected in a CFD Glorot, X.; and Bengio, Y. 2010. Understanding the diffi-
solver, we replaced the PLIC calculation model in LBfoam culty of training deep feedforward neural networks. In Pro-
with the trained DNN model with H = 2 and N = 48. ceedings of the Thirteenth International Conference on Ar-
Figure 8 shows simulation results based on the same pa- tificial Intelligence and Statistics, volume 9 of Proceedings
rameters and at the same timestep. The DNN was used for of Machine Learning Research, 249–256. JMLR Workshop
the simulation on the left, and the analytical PLIC was used and Conference Proceedings.
on the right. As can be seen, the simulation results are iden- Hirt, C.; and Nichols, B. 1981. Volume of fluid (VOF)
tical. Depending on the number of cells that require interface method for the dynamics of free boundaries. Journal of
reconstruction (i.e. the interface cells closer than δmax ), the Computational Physics 39(1): 201–225.
speedup gain for interface construction was between 1.5-3
times. Jafari, A.; Shirani, E.; and Ashgriz, N. 2007. An improved
three-dimensional model for interface pressure calculations
in free-surface flows. International Journal of Computa-
Conclusions tional Fluid Dynamics 21(2): 87–97.
In this paper, we present a DNN model to solve the forward Kawano, A. 2016. A simple volume-of-fluid reconstruction
problem of PLIC. We show that a DNN can learn to solve the method for three-dimensional two-phase flows. Computers
PLIC problem without a synthetic dataset, using its inverse & Fluids 134-135: 130–145.
problem. The DNN implementation of PLIC is up to sev-
eral orders of magnitude faster than analytical approaches, Kingma, D. P.; and Ba, J. 2015. Adam: A method for
for the same accuracy. The DNN implementation is easily stochastic optimization. In 3rd International Conference on
portable to various CFD solvers, and can be efficiently ac- Learning Representations, ICLR 2015, San Diego, CA, USA,
celerated by CPUs and GPUs in parallel. We implemented May 7-9, 2015, Conference Track Proceedings.
the DNN model in the LBfoam CFD solver to verify that the Kothe, D. B.; and Mjolsness, R. C. 1992. RIPPLE - A new
model can perform well in a practical setting. model for incompressible flows with free surfaces. AIAA
Journal 30(11): 2694–2700.
Acknowledgements Lafaurie, B.; Nardone, C.; Scardovelli, R.; Zaleski, S.; and
We thank the Natural Sciences and Engineering Research Zanetti, G. 1994. Modelling merging and fragmentation in
Council of Canada (NSERC) and Autodesk Inc. for their fi- multiphase flows with SURFER. Journal of Computational
nancial support. Physics 113(1): 134–147.
López, J.; Hernández, J.; Gómez, P.; and Faura, F. 2005. An
References improved PLIC-VOF method for tracking thin fluid struc-
Ataei, M.; Bussmann, M.; Shaayegan, V.; Costa, F.; Han, tures in incompressible two-phase flows. Journal of Com-
S.; and Park, C. B. 2021a. NPLIC: A Machine Learning putational Physics 208(1): 51–74.
Approach to Piecewise Linear Interface Construction. Nahed, J.; and Dgheim, J. 2020. Estimation curvature in
Ataei, M.; Shaayegan, V.; Costa, F.; Han, S.; Park, C. B.; PLIC-VOF method for interface advection. Heat and Mass
and Bussmann, M. 2021b. LBfoam: An open-source soft- Transfer 56(3): 773–787.
ware package for the simulation of foaming using the Lat- Nikolaev, A.; Richter, I.; and Sadowski, P. 2020. Deep learn-
tice Boltzmann Method. Computer Physics Communica- ing for climate models of the Atlantic ocean. In AAAI Spring
tions 259: 107698. Symposium: MLPS.
Balachandar, S.; and Eaton, J. K. 2010. Turbulent dispersed Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.;
multiphase flow. Annual Review of Fluid Mechanics 42(1): Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga,
111–133. L.; Desmaison, A.; Kopf, A.; Yang, E.; DeVito, Z.; Raison,
M.; Tejani, A.; Chilamkurthy, S.; Steiner, B.; Fang, L.; Bai,
Brenner, M. P.; Eldredge, J. D.; and Freund, J. B. 2019. Per- J.; and Chintala, S. 2019. PyTorch: An imperative style,
spective on machine learning for advancing fluid mechanics. high-performance deep learning library. In Wallach, H.;
Phys. Rev. Fluids 4: 100501. Larochelle, H.; Beygelzimer, A.; d'Alché-Buc, F.; Fox, E.;
Brunton, S. L.; Noack, B. R.; and Koumoutsakos, P. 2020. and Garnett, R., eds., Advances in Neural Information Pro-
Machine learning for fluid mechanics. Annual Review of cessing Systems, volume 32, 8026–8037. Curran Associates,
Fluid Mechanics 52(1): 477–508. Inc.
Figure 8: 3D foaming simulation using the LBfoam software. The disjoining pressure between bubbles is calculated using the
DNN (left) and the analytical PLIC (right).

Pfaff, T.; Fortunato, M.; Sanchez-Gonzalez, A.; and Yang, X.; and James, A. J. 2006. Analytic relations for
Battaglia, P. W. 2020. Learning mesh-based simulation with reconstructing piecewise linear interfaces in triangular and
graph networks. tetrahedral grids. Journal of Computational Physics 214(1):
41–54.
Raissi, M.; Yazdani, A.; and Karniadakis, G. E. 2020. Hid-
den fluid mechanics: Learning velocity and pressure fields Yoon, H.; Melander, D.; and Verzi, S. J. 2020. Permeability
from flow visualizations. Science 367(6481): 1026–1030. prediction of porous media using convolutional neural net-
works with physical properties. In AAAI Spring Symposium:
Rider, W. J.; and Kothe, D. B. 1998. Reconstructing volume MLPS.
tracking. Journal of Computational Physics 141(2): 112–
152. Youngs, D. L. 1982. Time-dependent multi-material flow
with large fluid distortion. Numerical Methods for Fluid Dy-
Rothman, D. H.; and Zaleski, S. 1994. Lattice-gas models namics .
of phase separation: interfaces, phase transitions, and multi-
phase flow. Rev. Mod. Phys. 66: 1417–1479.
Scardovelli, R.; and Zaleski, S. 2000. Analytical relations
connecting linear interfaces and volume fractions in rectan-
gular grids. Journal of Computational Physics 164(1): 228–
237.
Scheufler, H.; and Roenby, J. 2019. Accurate and efficient
surface reconstruction from volume fraction data on general
meshes. Journal of Computational Physics 383: 1–23.
Singh, A. P.; Medida, S.; and Duraisamy, K. 2017. Machine-
learning-augmented predictive modeling of turbulent sepa-
rated flows over airfoils. AIAA Journal 55(7): 2215–2227.
Skarysz, M.; Garmory, A.; and Dianat, M. 2018. An iterative
interface reconstruction method for PLIC in general convex
grids as part of a Coupled Level Set Volume of Fluid solver.
Journal of Computational Physics 368: 254–276.
Wang, Z.-B.; Chen, R.; Wang, H.; Liao, Q.; Zhu, X.; and Li,
S.-Z. 2016. An overview of smoothed particle hydrodynam-
ics for simulating multiphase flow. Applied Mathematical
Modelling 40(23): 9625–9655.