# Coupling of the Finite Element Method with Physics Informed Neural Networks
for the Multi-Fluid Flow Problem

Michel Nohra  Département de mathématiques et de génie industriel,
Polytechnique Montréal, Montréal, Québec, Canada, H3T 1J4  Steven Dufour
Département de mathématiques et de génie industriel, Polytechnique Montréal,
Montréal, Québec, Canada, H3T 1J4

###### Abstract

Multi-fluid flows are found in various industrial processes, including metal
injection molding and 3D printing. The accuracy of multi-fluid flow modeling
is determined by how well interfaces and capillary forces are represented. In
this paper, the multi-fluid flow problem is discretized using a combination of
a Physics-Informed Neural Network (PINN) with a finite element discretization.
To determine the best PINN formulation, a comparative study is conducted using
a manufactured solution. We compare interface reinitialization methods to
determine the most suitable approach for our discretization strategy. We
devise a neural network architecture that better handles complex free surface
topologies. Finally, the coupled numerical strategy is used to model a rising
bubble problem.

##  1 Introduction

Multi-fluid flows play an important role in various industrial processes,
ranging from 3D printing to oil extraction processes, and metal injection
applications. Understanding and accurately modeling these complex flows is
important for optimizing process efficiency, and for ensuring reliable
outcomes. The numerical modeling of multi-fluid flows is challenging due to
the dynamics of the interfaces between the different fluids, and the influence
of capillary force. Numerical methods for the study of multi-fluid flows can
be broadly classified into Eulerian and Lagrangian techniques. The Lagrangian
approach splits the domain into subdomains and discretizes the interfaces. The
nodes on each interface are transported by the fluid flow. This process
necessitates updating physical quantities at every time step, and it requires
the remeshing of the entire domain if elements become too degenerated.
Geometric quantities, such as the curvature of the interface and the normal
vectors to the interfaces, are obtained from the relative position of these
nodes, which may become problematic for a distorted discretization of the
interface. The capillary force is treated as a boundary condition for each
subdomain.

On the other hand, Eulerian techniques employ a different strategy for
modeling multi-fluid flows. A marker variable is used to locate each fluid,
and the evolution of the interfaces is governed by the mass conservation law

$\partial_{t}F+\bm{u}\cdot\nabla F=0,$

where $\bm{u}$ denotes the velocity of each fluid, and $F$ is the marker.
Eulerian methods treat fluid properties as continuous fields, defined in the
entire computational domain. Geometric properties, such as the curvature and
normal vectors to an interface, are calculated using the marker, where normal
vectors become $\bm{n}=\frac{\nabla F}{\|\nabla F\|}$ and the curvature of the
interface is given by $\kappa=\nabla\cdot\bm{n}$, thus approximating the
capillary forces by a volumetric force. Eulerian techniques are more robust
and more suitable for distorted $3$D interfaces, but suffer from numerical
errors related to the dissipation of the interface marker. With the Eulerian
approach, three primary methods are commonly used to simulate multi-fluid
flows: the volume-of-fluid (VOF) method, the level-set method, and the pseudo-
concentration method. Each method offers unique features and advantages to
accurately represent and capture fluid interfaces.

With the volume-of-fluid method, the marker takes values between zero and one,
representing the fraction of a given fluid within a computational cell. If the
marker value in a cell is $F=1$, then the cell contains fluid A; if the value
is $F=0$, the cell contains only fluid B. This marker is then transported
using a finite volume method. The cells containing an interface are such that
$0<F<1$, but the position of the interface inside the cell is unknown. After
transporting the marker, algorithms are used to reconstruct an interface and
extract geometric quantities such as the normal vectors and the curvature of
an interface.

The level-set method, on the other hand, represents the location of the two
fluids using a signed distance function to an interface. The zero contour of
the marker determines the location of an interface. The normal vectors are
given by $\bm{n}=\frac{\nabla F}{\|\nabla F\|},$ and the curvature of
interfaces by $\kappa=\nabla\cdot\bm{n}.$ The pseudo-concentration method,
sometimes referred to as the conservative level-set method, combines
components of both the level-set and the volume-of-fluid methods. This method
uses a distance function to an interface to approximate a smooth Heaviside
function that varies between zero and one, which is used as a marker. The
marker’s contour level of $0.5$ determines each interface position, and the
geometric quantities are obtained using the expressions used with the level-
set method. Notably, the pseudo-concentration method exhibits a steeper
gradient at the interface while maintaining a constant value away from it. In
this work, we will use the level-set method due to the smoothness of the
marker.

Traditional numerical methods, like the finite difference, the finite element,
and the finite volume methods, are commonly used to numerically model multi-
fluid flow problems. The discretizations employed by these methods result in
high computational and memory requirements for high-dimensional problems. In
addition, the second-order derivatives required to calculate the curvature are
a source of additional numerical errors.

Neural networks are revolutionary mathematical approximators that have been
effectively applied in various fields, including speech recognition and
classification tasks [19, 2, 10]. With recent advancements in computer
hardware, the use of neural networks in scientific computing has recently
grown in popularity. These neural networks offer promising alternatives for
numerical modeling and simulation tasks, leveraging their ability to learn
complex patterns and relationships from large datasets. Neural networks have
been applied to multiple scientific problems, such as turbulence model design
from data [17, 4, 18, 27] and for simulating heat transfer problems [8, 25].

Among the methods that use deep learning for scientific computing, Physics-
Informed Neural Networks (PINNs)[22] are a popular approach for simulating
physical problems. PINNs, first introduced by Raissi et al.[22], leverage all
the information available about the problem, such as the underlying PDE,
boundary conditions, physical properties, and known measurement data, to train
a neural network that best approximates the model. These networks transform
boundary-value problems into optimization problems, providing flexibility to
include any constraint in the problem.

Since their introduction, PINNs have been applied to a variety of forward,
inverse, and parameterized problems [8, 13, 25, 16, 15, 29, 20], including
fluid dynamics problems. However, due to the second-order, non-linear and
saddle point nature of fluid dynamics problems, traditional numerical methods
still outperform PINNs for forward problems (well-posed boundary-value
problems). PINNs become interesting when experimental data are available [7].

In the case of hyperbolic equations, such as the transport equation,
traditional methods, such as the finite element method, often exhibit
instabilities that necessitate special treatment. This is where PINNs offer an
advantage, as they are not prone to the same instabilities, and they do not
require additional measures to ensure numerical stability. Therefore, PINNs
can be seen as a valuable alternative for the accurate discretization of
hyperbolic equations. Another advantage of PINNs over traditional methods is
that the numerical approximations are continuous and infinitely
differentiable, which is important for an accurate approximation of surface
tension. For the reasons listed above, Physics-Informed Neural Networks will
be used for the discretization of the transport equation of the level-set
function, and the finite element method will be used to discretize the Navier-
Stokes equation.

In this work, we will describe the multi-fluid flow problem and the governing
equations, the FEM for the Navier-Stokes equations, and the PINNs.
Additionally, we will propose three different time-integration schemes for the
advection equation with PINNs. We will present a new approach, based on the
multi-level PINN [3], to obtain a more accurate approximation of the level-set
function. We introduce three different methods to reinitialize the level-set
function. And finally, validate the proposed methodologies on various problems
with manufactured solutions, and on the rising bubble problems.

##  2 Governing equations

Let us consider a multi-fluid flow problem with two non-miscible and
incompressible fluids A and B. The density and viscosity of each fluid are
given by $\rho_{\text{A}},\mu_{\text{A}}$, $\rho_{\text{B}}$ and
$\mu_{\text{B}}$, and are defined over subdomains $\upOmega_{\text{A}}$ and
$\upOmega_{\text{B}}$, where
$\upOmega_{\text{A}}\cup\upOmega_{\text{B}}=\upOmega$, which is the whole
domain, and $\upOmega_{\text{A}}\cap\upOmega_{\text{B}}={\mathrm{\Gamma}},$
where $\mathrm{\Gamma}$ is the interface between the two fluids. With the
level-set method, we define a signed distance function to the interface
$\mathrm{\Gamma},$ denoted by $F$, where $F$ is negative in
$\upOmega_{\text{A}},$ positive in $\upOmega_{\text{B}}$, and zero at the
interface. We can define an approximate Heaviside function using the level-set
function $F$. Various approximations can be used, including the sigmoid
function $H(F)=\sigma(\alpha F)$, where $\alpha$ determines the sharpness of
$H$, which depends on the grid size. Using the level-set function, and the
approximate Heaviside function $H$, we can build a one-fluid model with
variable physical quantities:

$\displaystyle\rho(F)=\rho_{\text{A}}H(F)+\rho_{\text{B}}(1-H(F));$
$\displaystyle\mu(F)=\mu_{\text{A}}H(F)+\mu_{\text{B}}(1-H(F)).$

Using this one-fluid model, we obtain the set of equations that represent the
continuity or incompressibility equation, the conservation of momentum
equations, and the marker advection equation:

$\displaystyle\quad\nabla\cdot\bm{u}=0;$ (1)
$\displaystyle\quad\rho(F)\frac{\partial\bm{u}}{\partial
t}+\rho(F)(\bm{u}\cdot\nabla)\bm{u}+\nabla
p-\nabla\cdot\nabla(\mu(F)\bm{u})=\bm{F}_{\text{ext}};$
$\displaystyle\quad\partial_{t}F+\bm{u}\cdot\nabla F=0,$

where $\bm{u}$ represents the fluids velocity, $p$ represents their pressure,
$F$ is the level-set function, $\bm{F}_{\text{ext}}$ represents external
forces, $\mu$ is the viscosity and $\rho$ is the density of the fluids.

$\bm{F}_{\text{ext}}$ includes an approximation of the capillary force, which
is a surface force at the interface,

$\bm{f_{\text{s}}}=\delta\gamma\kappa\bm{n}_{\mathrm{\Gamma}},$

where $\delta$ is the Dirac function defined at the interface, $\gamma$ is the
surface tension coefficient, $\kappa$ is the interface curvature and
$\bm{n}_{\mathrm{\Gamma}}$ is a normal vector at the interface. Using the
continuum surface force approximation [6], we can approximate the normal
vectors using $\bm{n}_{\mathrm{\Gamma}}\ =\frac{\nabla F}{\|\nabla F\|}$, the
Dirac function by $\delta=\frac{\partial H}{\partial F}$, and the curvature of
the interface by $\kappa=\nabla\cdot\bm{n}_{\mathrm{\Gamma}},$ and obtain a
volumetric approximation of the capillary force.

###  2.1 Adimensional equations

Adimensionalizing, often referred to as non-dimensionalization, is a
fundamental technique in the fields of physics and engineering that helps
reduce the number of parameters in a problem, allowing for a clearer
understanding of the underlying dynamics and relationships. Results derived
from dimensionless equations are often more general, as they are not tied to
specific units of measurement. This universality facilitates easier comparison
across different systems and scales, leading to broader insights and more
robust solutions. Furthermore, adimensionalizing the problem means that the
calculated variables would tend to be in the order of unity, making the
problem easier from a numerical point of view. The quantities of interest
would be reformulated as follows:

$\rho=\hat{\rho}\bar{\rho}$ | $\mu=\hat{\mu}\bar{\mu}$ | $\bm{u}=\hat{u}\bar{\bm{u}}$ | $L=\hat{L}\bar{L}$
---|---|---|---
$t=\frac{\hat{L}}{\hat{u}}\bar{t}$ | $\partial_{t}=\frac{\hat{u}}{\hat{L}}\partial_{\bar{t}}$ | $\nabla=1/\hat{L}\bar{\nabla}$ | $p=\hat{\rho}\,\hat{u}^{2}\,\bar{p}$
Table 1: Adimensional properties.

where $\bar{A}$ represent the non-dimensional value of $A$ and $\hat{A}$ the
characteristic value of $A$. In the remainder of this work, all variables will
be presented in non-dimensional form, and the hat notation ($\hat{\cdot}$)
will be omitted for simplicity. After substituting in the differential
equation, we obtain the following non-dimensional numbers:

  * •

Reynolds number $Re=\frac{\rho uL}{\mu}$ and it signifies the importance of
inertial forces over viscous forces. A higher Reynolds number means the flow
exhibits more turbulent behavior, and from a computational point of view, it
means that our problem has less diffusion, and necessitates a finer mesh for a
stable solution.

  * •

Weber number $We=\frac{\rho u^{2}L}{\sigma}$ and it represents the importance
of inertial forces over capillary forces. A higher Weber number means that the
interface has a greater tendency to break and contains higher curvatures.

  * •

Froud number $Fr=\frac{u}{\sqrt{gL}}$ and represent the importance of inertial
forces over gravitational forces.

The complete non-dimensional equations become:

$\displaystyle\nabla\cdot\bm{u}=0;$ (2)
$\displaystyle\rho(F)\frac{\partial\bm{u}}{\partial
t}+\rho(F)(\bm{u}\cdot\nabla)\bm{u}+\nabla
p-\frac{1}{Re}\mu(F)\nabla^{2}\bm{u}=\frac{1}{We}\frac{\nabla F}{||\nabla
F||}\nabla\cdot(\frac{\nabla F}{||\nabla F||})+\frac{1}{Fr^{2}}\rho(F);$
$\displaystyle\partial_{t}F+\bm{u}\cdot\nabla F=0.$

##  3 Numerical methods

In this work, we combine two numerical methods: the Finite Element Method
(FEM) for the discretization of the fluid dynamics problem, to determine the
fluid’s velocity and pressure, and the Physics-Informed Neural Network (PINN)
for the advection of the level-set function.

###  3.1 The finite element method for the Navier-Stokes equations

The level-set method for modeling multi-fluid flows involves the numerical
discretization of nonlinear, coupled partial differential equations. This
process typically progresses through discrete time steps. At each time step,
the Navier-Stokes equations are discretized based on the previously computed
level-set function. The level-set function is then advected using the computed
velocity field.

In this section, we explain the discretization of the Navier-Stokes equations
using the finite element method. Let us define the functional space of square-
integrable functions

$L^{2}(\upOmega)=\\{u:\upOmega\xrightarrow{}\mathbb{R}|\int_{\upOmega}u^{2}dx<\infty\\},$

with the scalar product

$(u,v)_{0,\upOmega}=\int_{\upOmega}u\,v\,dx.$

The vectorial and tensorial versions are

$L^{2}(\upOmega)^{n}=\\{\bm{u}:\upOmega\xrightarrow{}\mathbb{R}^{n}|\int_{\upOmega}\bm{u}\cdot\bm{u}\
dx<\infty\\},$

and

$L^{2}(\upOmega)^{(n\times
n)}=\\{\bm{\psi}:\upOmega\xrightarrow{}\mathbb{R}^{(n\times
n)}|\int_{\upOmega}\bm{\psi}:\bm{\psi}\ dx<\infty\\},$

where ’$:$’ is the double contraction operator. Their respective scalar
product are

$(\bm{u},\bm{v})_{0,\upOmega}=\int_{\upOmega}\bm{u}\cdot\bm{v}\ dx,$

and

$(\bm{\psi},\bm{\phi})_{0,\upOmega}=\int_{\upOmega}\bm{\psi}:\bm{\phi}\ dx.$

Next, we define the Sobolev spaces

$H^{1}(\upOmega)=\\{u\in L^{2}(\upOmega)\ |\ \partial_{x_{i}}u\in
L^{2}(\upOmega),i=1,2,...,n\\},$

and

$H^{1}(\upOmega)^{n}=\\{\bm{u}\in L^{2}(\upOmega)^{n}\ |\
\partial_{x_{i}}\bm{u}_{j}\in L^{2}(\upOmega),\ i,j=1,2,...,n\\},$

with the inner product

$(u,v)_{1,\upOmega}=(u,v)_{L^{2}(\upOmega)}+\sum_{i=1}^{n}(\partial_{x_{i}}u,\partial_{x_{i}}v)_{L^{2}(\upOmega)},$

and finally, for imposing boundary conditions, we define the Sobolev space

$H^{1}_{\text{D}\bm{u}}(\upOmega)^{n}=\\{\bm{u}\in H^{1}(\upOmega)^{n}\ |\
\bm{u}(\bm{x})=0\ \text{for}\ \bm{x}\in\mathrm{\Gamma}_{\text{D}\bm{u}}\\},$

and

$H^{1/2}(\mathrm{\Gamma}_{0})^{n}=\\{\gamma_{0}(\bm{u})\in
L^{2}({\rm\Gamma}_{0})^{n}|\ \bm{u}\in H^{1}(\upOmega)^{n}\\},$

where $\gamma_{0}(\bm{u})$ represents the trace of $\bm{u}$ on
${\rm\Gamma}_{0}$.

The weak formulation of the equations is obtained by multiplying these
equations by test functions belonging to appropriate functional spaces, and
integrating by parts over $\upOmega$. The mass conservation equation becomes

$(q,\nabla\cdot\bm{u})_{0,\upOmega}=0,\ \ \ \forall q\in L^{2}(\upOmega).$ (3)

Given $\bm{f}\in L^{2}(\upOmega)^{n}$, $\bm{u}_{\text{D}\bm{u}}\in
H^{1/2}(\mathrm{\Gamma}_{\text{D}\bm{u}})^{n}$ and $\bm{t}_{\text{N}\bm{u}}\in
H^{-1/2}(\mathrm{\Gamma}_{\text{N}\bm{u}})$, where $H^{-1/2}$ is the dual
space of $H^{1/2}$, we obtain the weak formulation of the momentum equation,

$\displaystyle\big{(}\rho\frac{\partial\bm{u}}{\partial
t},\bm{v}\big{)}_{0,\upOmega}+(\mu\nabla\bm{u},\nabla\bm{v})_{0,\upOmega}+$
$\displaystyle(\rho(\bm{u}\cdot\nabla)\bm{u},\bm{v})_{0,\upOmega}-(p,\nabla\cdot\bm{v})_{0,\upOmega}$
(4) $\displaystyle=$ $\displaystyle\ (\bm{f},\bm{v})_{0,\upOmega}+\langle
t_{\text{N}\bm{u}},v\rangle_{H^{-1/2},H^{1/2}},\quad\forall\>\bm{v}\in
H^{1}_{\text{D}\bm{u}}(\upOmega)^{n}.$

The bilinear form obtained for the weak formulations in equations (3) and (4)
can be written as:

$\displaystyle A(\bm{u},\bm{v})+B(\bm{v},p)$ $\displaystyle=G(\bm{v})\ ,\ \
\forall\bm{v}\in V;$ $\displaystyle B(\bm{u},q)$ $\displaystyle=0\ ,\ \
\forall q\in Q,$

where $V$ and $Q$ are Hilbert spaces, and $A$ and $B$ are bilinear forms given
by

$\displaystyle A(\bm{u},\bm{v})=\big{(}\rho\frac{\partial\bm{u}}{\partial
t},\bm{v}\big{)}_{0,\upOmega}+(\mu\nabla\bm{u},\nabla\bm{v})_{0,\upOmega}+(\rho(\bm{u}\cdot\nabla)\bm{u},\bm{v})_{0,\upOmega};$
$\displaystyle B(\bm{v},p)=-(p,\nabla\cdot\bm{v})_{0,\upOmega},$

and $G$ is a linear form,

$G(v)=(\bm{f},\bm{v})_{0,\upOmega}.$

A unique solution exists according to the Lax-Milgram theorem, and the
Ladyzhenskaya-Babuska-Brezzi (LBB) condition is satisfied, if the following
conditions are met:

  * •

Coercivity: This refers to the property that the bilinear form $A(u,v)$ is
coercive on the solution space V:

$\exists\alpha>0\text{ such that }\forall\>u\in V,\
A(u,u)\geq\alpha\|u\|_{V}^{2};$

  * •

Continuity (or boundedness): This condition ensures that the bilinear form is
continuous:

$|A(u,v)|\leq M\|u\|_{V}\|v\|_{V};$ $|B(u,q)|\leq N\|u\|_{V}\|q\|_{Q},$

for all $u,v\in V$, $q\in Q$, for positive constants $M$ and $N$.

  * •

LBB (inf-sup) condition: Ensuring that the chosen spaces for velocity and
pressure are compatible is crucial. This compatibility condition is often
expressed via the inf-sup condition,

$\inf_{v\in V}\sup_{q\in Q}\frac{|B(v,q)|}{\|v\|_{V}\|q\|_{Q}}>\beta,$

where $\beta$ is a positive constant.

In the finite element setting, the domain $\upOmega$ is discretized with
smaller subdomains, called elements. We take the discrete functional space
$V_{h}\subset H^{1}_{\text{D}\bm{u}}(\upOmega)^{n}$ and $Q_{h}\subset
L^{2}(\upOmega)$ of finite dimensions, with basis functions $\bm{\phi}_{i}$
and $\eta_{i}$ respectively, defined on each element, to express a discrete
approximation of the dependent variables as:

$\displaystyle\bm{u}_{h}$
$\displaystyle=\sum_{i=1}^{N_{u}}u_{i}\bm{\phi}_{i}(x,y);$ $\displaystyle
p_{h}$ $\displaystyle=\sum_{i=1}^{N_{p}}p_{i}\eta_{i}(x,y),$

where $\bm{u}_{h}$ is the discrete approximation of the velocity field, and
$p_{h}$ is the discrete approximation of the pressure. With the Galerkin
finite element method, we consider $\bm{v}_{h}\in V_{h}$ and $q_{h}\in Q_{h}$.
As for the choice of the subspaces $Q_{h}$ and $V_{h}$, they need to satisfy
the LBB condition defined above. An example of such spaces is the Taylor-Hood
element combination, where the $\bm{\phi}_{i}$ are the quadratic Lagrange
elements $P_{2}$, and the $\eta_{i}$ are the linear Lagrange elements $P_{1}$.

###  3.2 Physics-Informed Neural Networks

To gain a better understanding of neural networks, let us begin with a simple
example of a one-layer neural network. If we consider the input vector
$\bm{x}$ and the output vector $\bm{x}_{1}$ to have dimensions $n$ and $m$
respectively, the output can be expressed as

$\bm{x}_{1}=\sigma(W\bm{x}+\bm{b}),$

where $W$ is a weight matrix of dimension $m\times n$ and $\bm{b}$ is a bias
vector of dimension $m$. The activation function $\sigma(\bm{z})$ is user-
defined and can take various forms, such as a cosine or a hyperbolic tangent.
This function is applied to each component of a vector $\bm{z}$.

By stacking multiple layers, we can establish a feedforward neural network,
which involves passing the output of one layer to the input of the next layer,
as illustrated in figure 1. The output of the network is then expressed as

$N(\bm{x})=\sigma_{l}\big{(}W_{n}\sigma_{l-1}(W_{l-1}\sigma_{l-2}(\ldots(\sigma_{1}(W_{1}\bm{x}+\bm{b}_{1}))\ldots)+\bm{b}_{l-1})+\bm{b}_{l}\big{)}.$

Figure 1: Feedforward neural network

One can build a more complex neural network by changing the connections (the
architecture) between the layers.

Once the neural network has been constructed with the desired architecture,
the next step is to determine the optimal values for its parameters
$\bm{\zeta},$ a vector that contains all the weights and the biases of the
layers of the neural network, so that it produces the desired output. This
process of fine-tuning the weights is referred to as “training”. Training
involves minimizing the disparity between the actual output generated by the
network and the desired output. To formulate this as an optimization problem,
we aim to minimize a specific loss function, denoted as $L$, which is user-
defined. For example, for regression problems, the Mean Squared Error (MSE) is
commonly used as a loss function, while for classification tasks, the Cross-
Entropy loss function is often preferred. The optimization problem can then be
expressed as

$\min_{\bm{\zeta}}\ L(N_{\bm{\zeta}}(\bm{x}),D(\bm{x})),$ (5)

where $N_{\bm{\zeta}}(\bm{x})$ represents the output produced by the neural
network for a given ${\bm{\zeta}}$. $D(\bm{x})$ is the target, or desired
output. $L$ is the user-defined loss function that quantifies the
dissimilarity between the predicted and the actual outputs. While various
optimization methods could address this minimization problem, in the context
of this work, we will only consider gradient-based optimization methods.
Gradient-based methods leverage the property that the gradient of a function
points in the direction of the function’s steepest ascent. Therefore, to
minimize the function, one would generally move in the direction opposite to
the gradient, with a step size proportional to a user-defined learning rate.

Popular gradient-based methods include the Stochastic Gradient Descent (SGD),
the Limited-memory Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS) [5], and
the Adaptive Moment Estimation method (Adam) [14]. With neural networks, the
gradient is calculated using the back-propagation algorithm, which employs the
chain rule and automatic differentiation to compute the partial derivatives of
the network’s output with respect to the weights, the biases, and the
network’s input.

With the original PINN, introduced by Raissi et al. [22], which is usually
referred to as the vanilla PINN, we approximate the solution of a boundary-
value problem using a feedforward neural network. The derivatives of the
approximation are obtained by automatic differentiation. Let us consider a
domain $\upOmega$ with a boundary $\mathrm{\Gamma}=\partial\upOmega$. We want
to find an approximation of the solution $u(\bm{x})$ of a boundary-value
problem:

$\displaystyle R(\bm{x},u)=0\ ,\ \forall\bm{x}\in\upOmega;$ $\displaystyle
B(\bm{x},u)=0\ ,\ \forall\bm{x}\in\mathrm{\mathrm{\Gamma}},$

where $R$ is the residual of a PDE, and $B$ is the boundary condition on
$u(\bm{x})$ we want to impose on $\rm\Gamma.$ With the vanilla PINN, we
approximate the solution $u(\bm{x})$ by a feedforward neural network
$N_{\bm{\zeta}}(\bm{x})$, and we tune the weights and biases $\bm{\zeta}$ by
minimizing the loss function

$\displaystyle L$
$\displaystyle=\alpha_{1}\sum_{i=1}^{n_{\upOmega}}R(\bm{x}_{\upOmega}^{i},N_{\bm{\zeta}}(\bm{x}_{\upOmega}^{i}))^{2}+\alpha_{2}\sum_{j=1}^{n_{\mathrm{\Gamma}}}B(\bm{x}_{\mathrm{\Gamma}}^{j},N_{\bm{\zeta}}(\bm{x}_{\mathrm{\Gamma}}^{j}))^{2}+\alpha_{3}\sum_{k=1}^{n_{\rm
D}}(N_{\bm{\zeta}}(\bm{x}_{\rm D}^{k})-u(\bm{x}_{\rm D}^{k}))^{2}$
$\displaystyle=\alpha_{1}\|R(\bm{x}_{\upOmega},N_{\bm{\zeta}}(\bm{x}_{\upOmega}))\|+\alpha_{2}\|B(\bm{x}_{\mathrm{\Gamma}},N_{\bm{\zeta}}(\bm{x}_{\mathrm{\Gamma}}))\|+\alpha_{3}\|N_{\bm{\zeta}}(\bm{x}_{\rm
D})-u(\bm{x}_{\rm D})\|$

where $\\{\bm{x}_{\upOmega}^{i}\\}_{i=1}^{n_{\upOmega}}$ is a set of
$n_{\upOmega}$ points that belong to $\upOmega$,
$\\{\bm{x}_{\mathrm{\Gamma}}^{i}\\}_{i=1}^{n_{\mathrm{\Gamma}}}$ is a set of
$n_{\mathrm{\Gamma}}$ points that belong to $\mathrm{\Gamma}$, $\\{\bm{x}_{\rm
D}^{i}\\}_{i=1}^{n_{\rm D}}$ is a set of $n_{\rm D}$ points where the solution
is known (experimental measurements), $\alpha_{i}$ are weighting parameters
that are user-defined, and $\|\cdot\|$ is the discrete $L2$ norm.

###  3.3 Proposed PINN for the level-set advection problem

In this work, we will limit our use of PINN to approximate the level-set
function $F$, while the finite element method will give an approximation of
the velocity field $\bm{u}$ and pressure $p$.

####  3.3.1 Time integration schemes for the advection equation

Three approaches are considered in this work to handle time evolution. The
first two are continuous-in-time, and the third approach is discrete-in-time.
For the continuous-in-time approach, the time domain is decomposed into
smaller time intervals with a time step size ${\rm\Delta}t$. At each time step
$t_{n}$, we compute the FEM approximation of the velocity field $\bm{u}_{n}$.
This velocity field is then used to obtain an approximation of the level-set
function $F(\bm{x},t)$, defined over the entire time interval
$t\in[t_{n-1},t_{n}]$.

For a small time step size ${\rm\Delta}t$, the time interval becomes too
small. Neural networks perform better during training if the inputs are
normalized and have a magnitude of $O(1)$ [23], hence the need to normalize
the time input. Therefore, we define a normalized time
$t^{*}=\frac{t-t_{n-1}}{{\rm\Delta}t}$, which falls within the interval
$t^{*}\in[0,1]$ when $t\in[t_{n-1},t_{n}]$. The advection equation of the
level-set function becomes

$\partial_{t^{*}}F_{n}+{\rm\Delta}t\>\bm{\overline{u}}_{n}\cdot\nabla
F_{n}=0,\ \text{where}\ F_{n}(\bm{x},t^{*}=0)=F_{n-1}(\bm{x},t^{*}=1),$ (6)

where $F_{n-1}$ is the level-set function defined over the previous time
interval $[t_{n-2},t_{n-1}]$, $F_{n}$ is the discrete level-set function
defined over the current time interval $[t_{n-1},t_{n}]$, and
$\bm{\overline{u}}_{n}$ is the continuous advection velocity defined for all
$t\in[t_{n-1},t_{n}]$, obtained by interpolating the discrete velocities
calculated with the FEM $\bm{u}_{n},\bm{u}_{n-1},\dots\bm{u}_{0}$ in time. For
example, we can define
$\bm{\overline{u}}_{n}=\bm{u}_{n-1}(1-t^{*})+\bm{u}_{n}t^{*}$ as the linear
interpolation of the velocity fields. Figure 2 gives a graphical
representation of two time intervals, showing the discrete times and
velocities $t_{n}$ and $\bm{u}_{n}$, the level-set function $F_{n}$ defined on
each interval, the interpolated velocity $\bm{\overline{u}}_{n}$, and the
normalized time $t^{*}$ for each interval.

$t_{n-2}$$t_{n-1}$$t_{n}$$\bm{u}_{n-2}$$\bm{u}_{n-1}$$\bm{u}_{n}$$\bar{\bm{u}}_{n-1},F_{n-1}$$\bar{\bm{u}}_{n},F_{n}$$t^{*}\in[0,1]$$t^{*}\in[0,1]$
Figure 2: A representation of two time intervals, showing the discrete times
$t_{n}$, the discrete velocities $\bm{u}_{n}$, the continuous in time
interpolation of the velocities $\bar{\bm{u}}_{n}$, and the level-set function
$F_{n}$ defined over each time interval.

With the first continuous-in-time approach, which we will refer to as the weak
PINN, the initial conditions are weakly imposed. The neural network will take
the space-time coordinates $[x,y,z,t^{*}]$ as input, and the network is
trained using the loss function

$L=\|\partial_{t^{*}}F_{n}(\bm{x},t^{*})+{\rm\Delta}t\>\bm{\overline{u}}_{n}\cdot\nabla
F_{n}(\bm{x},t^{*})\|+\|F_{n}(\bm{x},t^{*}=0)-F_{n-1}(\bm{x},t^{*}=1)\|,$

where $\bm{\overline{u}}_{n}=\bm{u}_{n-1}(1-t^{*})+\bm{u}_{n}t^{*}$. We used a
linear interpolation for $\overline{\bm{u}},$ but higher-order interpolations
are possible.

With the second continuous-in-time approach, which we refer to as the strong
PINN, we strongly impose the initial conditions as proposed in [24], where the
output of the neural network $N_{\bm{\zeta}}(\bm{x},t^{*})$ is multiplied by a
distance function to the boundary $d(\bm{x},t^{*})$, and we add a lift
function $V(\bm{x},t^{*})$ that satisfies the desired boundary condition. In
our case, the distance function is $d(\bm{x},t^{*})=t^{*}$,
$V(\bm{x},t^{*})=F_{n-1}(\bm{x},t^{*}=1)\ (1-t^{*})$, and the level-set
function can be expressed as

$F_{n}(\bm{x},t^{*})=N_{\bm{\zeta}}(\bm{x},t^{*})\ t^{*}+F_{n-1}(\bm{x},1)\
(1-t^{*}),$ (7)

on $[t_{n-1},t_{n}].$ With this formulation, the level-set function at
$t^{*}=0$ satisfies the initial condition $F_{n}(\bm{x},0)=F_{n-1}(\bm{x},1),$
and the loss function becomes

$L=\|\partial_{t^{*}}F_{n}(\bm{x},t^{*})+{\rm\Delta}t\
\bm{\overline{u}}_{n}\cdot\nabla F_{n}(\bm{x},t^{*})\|.$

The reason why we define $V(\bm{x},t^{*})=F_{n-1}(\bm{x},1)\ (1-t^{*})$ is to
ensure that the level-set function $F_{n}$ only depends on the previous level-
set function $F_{n-1}$. If we would consider
$V(\bm{x},t^{*})=F_{n-1}(\bm{x},1),$ for example, we would obtain that

$\displaystyle F_{n}(\bm{x},t^{*})$
$\displaystyle=N_{\bm{\zeta},n}(\bm{x},t^{*})\ t^{*}+F_{n-1}(\bm{x},1)$
$\displaystyle=N_{\bm{\zeta},n}(\bm{x},t^{*})\
t^{*}+N_{\bm{\zeta},n-1}(\bm{x},1)+F_{n-2}(\bm{x},1)$
$\displaystyle=N_{\bm{\zeta},n}(\bm{x},t^{*})\
t^{*}+N_{\bm{\zeta},n-1}(\bm{x},1)+\ldots+N_{\bm{\zeta},1}(\bm{x},1)+F_{0}(\bm{x},1),$

where $N_{\bm{\zeta},n}$ is the neural network trained over the time interval
$[t_{n-1},t_{n}].$ On the other hand, for $V(\bm{x},t^{*})=F_{n-1}(\bm{x},1)\
(1-t^{*})$, we obtain that

$\displaystyle F_{n}(\bm{x},t^{*})$
$\displaystyle=N_{\bm{\zeta},n}(\bm{x},t^{*})\
t^{*}+F_{n-1}(\bm{x},1)(1-t^{*})$
$\displaystyle=N_{\bm{\zeta},n}(\bm{x},t^{*})\
t^{*}+(N_{\bm{\zeta},n-1}(\bm{x},1)\,(1)+F_{n-2}(\bm{x},1)\,(1-1))(1-t^{*})$
$\displaystyle=N_{\bm{\zeta},n}(\bm{x},t^{*})\
t^{*}+N_{\bm{\zeta},n-1}(\bm{x},1)(1-t^{*}).$

This results in less memory requirements, since we only need to keep the
parameters of two neural networks in memory, $N_{\bm{\zeta},n}$ and
$N_{\bm{\zeta},n-1}$, for each time interval, and the partial derivatives at
$t^{*}=1$, which are used to approximate the capillary forces, will depend
exclusively on the output of one neural network, ${\nabla
F_{n}(\bm{x},1)=\nabla N_{\bm{\zeta},n}(\bm{x},1)}$.

With the discrete-time approach, referred to as the finite difference PINN
(FDPINN), we use a finite difference approximation of the time derivative. The
neural network inputs are the spatial coordinates $[x,y,z]$. An implicit
scheme is used, and we consider the advection velocity to be
$\overline{\bm{u}}=\bm{u}_{n}.$ The loss function is then defined as

$L=\Big{\|}\frac{F_{n}(\bm{x})-F_{n-1}(\bm{x})}{{\rm\Delta}t}+\bm{u_{n}}\cdot\nabla
F_{n}(\bm{x})\Big{\|}.$

####  3.3.2 Neural network architecture

Neural networks exhibit a spectral bias, where the associated optimization
problem tends to converge faster toward the low-frequency components of the
solution [21, 9]. This spectral bias causes the problem to converge slowly
when the fluid interface undergoes significant topological changes, such as
bubble splitting, leading to prolonged computational times, and reduced
accuracy. Multiple methods have been proposed to address the spectral bias of
neural networks, most notably the Fourier Feature Neural Network [26], where
the presence of high-frequencies in the output is dealt with by mapping the
input space to a high-frequency component space. However, this method did not
produce satisfactory results for multi-fluid flow problems for two reasons.
First, these methods work better for cases where the solution has a uniform
frequency. This is not the case for multi-fluid flows, since the interface
could be smooth in one part of the domain, and the problem can become stiff in
another part of the domain. The second reason is that introducing high
frequencies in the input vector will introduce high frequencies in the output
vector, resulting in inaccurate free surface curvature approximations.

In this work, we use a multi-level neural network, inspired by the work of
Aldirany et al. [3]. We use two PINNs, giving $F_{n,1}$ and $F_{n,2}$ such
that

$F_{n}=F_{n,1}+\beta F_{n,2},$

where $\beta$ is a normalization parameter used to normalize the function
$F_{n,2}.$ The first function $F_{n,1}(\bm{x},t),$ approximates the level-set
function across the entire domain. The input of the neural network is the
space-time coordinates of points scattered over the entire domain. In this
work, we choose the nodes of the FEM mesh as collocation points and they are
assigned with random time coordinates $t^{*}$ within the interval
$t^{*}\in[0,1]$.

The PINN is trained for the boundary-value problem

$\partial_{t^{*}}F_{n,1}+{\rm\Delta}t\ \bm{\overline{u}}_{n}\cdot\nabla
F_{n,1}=0,\ \text{where}\ F_{n,1}(\bm{x},0)=F_{n-1}(\bm{x},1).$ (8)

The second function $F_{n,2}$ computes a higher-quality approximation near the
interface. Next, we define the input vector of the PINN that gives $F_{n,2},$
the boundary-value problem for $F_{n,2},$ and the normalization parameter
$\beta.$

Since interfaces evolve in their normal direction, adding the interface normal
vectors to the input of the PINN will help the network learn the evolution of
the free surface. The normal vectors to the interfaces are approximated using
$\nabla F_{n,1}$. Hence, the inputs of the second PINN, which produces
$F_{n,2},$ are the space-time coordinates, and the gradient of the approximate
solution $\nabla F_{n,1}$ at points around interfaces. These points are chosen
such that $|F_{n,1}(\bm{x},t^{*})|<\epsilon,$ where $\epsilon$ is a positive
constant that we choose to be $\epsilon=0.05.$ By replacing $F_{n}$ in
equation (6), we obtain

$\partial_{t}^{*}(F_{n,1}+\beta
F_{n,2})+{\rm\Delta}t\bar{\bm{u}}_{n}\cdot\nabla(F_{n,1}+\beta F_{n,2})=0,\
\text{where}\ F_{n,1}(\bm{x},0)+\beta F_{n,2}(\bm{x},0)=F_{n}(\bm{x},1),$

and by isolating $F_{n,2}$, we obtain the boundary-value problem

$\displaystyle\partial_{t}^{*}\beta
F_{n,2}+{\rm\Delta}t\bar{\bm{u}}_{n}\cdot\nabla\beta F_{n,2}$
$\displaystyle=-\partial_{t^{*}}F_{n,1}-{\rm\Delta}t\
\bar{\bm{u}}_{n}\cdot\nabla F_{n,1}$ (9) $\displaystyle=-R(F_{n,1}),\
\text{where}\ \beta F_{n,2}(\bm{x},0)=F_{n}(\bm{x},1)-F_{n,1}(\bm{x},0),$

where $R(F_{n,1})$ is the residual of $F_{n,1}$.

N.B: to accelerate the computation of the partial derivative of $F_{n,2}$ with
respect to $\bm{x}$, we can apply the chain rule

$\frac{\partial F_{n,2}(\bm{x})}{\partial x_{i}}=\frac{\partial
F_{n,2}(\bm{x},\nabla F_{n,1})}{\partial x_{i}}+\sum_{j=1}^{d}\frac{\partial
F_{n,2}(\bm{x},\nabla F_{n,1})}{\partial(\partial F_{n,1}(\bm{x})/\partial
x_{j})}\frac{\partial^{2}F_{n,1}(\bm{x})}{\partial x_{i}\partial x_{j}},$

and exclude $F_{n,1}$ from the computational graph.

As for the normalization of the function $F_{n,2}$ to obtain more accurate
results, we need to approximate the value of $\beta$, so that $F_{n,2}$ is of
the order of $O(1)$. We approximate $\beta$ using a finite difference scheme
for equation (9):

$\frac{\beta(F_{n,2}(\bm{x},1,\nabla F_{n,1})-F_{n,2}(\bm{x},0,\nabla
F_{n,1}))}{1-0}+{\rm\Delta}t\>\bar{\bm{u}}_{n}\cdot\nabla\beta
F_{n,2}(\bm{x},0,\nabla F_{n,1})=-\>R(F_{n,1}).$

Considering that $F_{n,1}(\bm{x},t)$ exactly satisfies the initial condition
$F_{n}(\bm{x},1)=F_{n,1}(\bm{x},0)$, which is the case for the “strong PINN”,
we obtain from (9) that $F_{n,2}(\bm{x},0,\nabla F_{n,1})=0$, leading to

$\max(\beta F_{n,2}(\bm{x},1,\nabla
F_{n,1}))\approx\beta\approx\max(|R(F_{n,1})|).$

The steps to obtain the level-set function approximation using the proposed
method are:

  * •

Compute the capillary forces using $F_{n-1}(\bm{x},1)$;

  * •

Compute the velocity $\bm{u}_{n}$ using the Finite Element Method;

  * •

Interpolate the velocities to obtain $\bar{\bm{u}}_{n}$ over
$t\in[t_{n-1},t_{n}]$;

  * •

Determine the first approximation of the level-set function $F_{n,1}$ over the
time interval $t^{*}\in[0,1]$ using the boundary value problem given by (8)
using PINN;

  * •

Compute $\beta=\max(R(F_{n,1}))$;

  * •

Compute a better approximation around the interface by approximating the
solution $F_{n,2}$ to the boundary value problem (9), using a PINN with the
input $[\bm{x},\nabla F_{n,1}]$;

  * •

Compute the final approximation of the level-set function $F_{n}=F_{n,1}+\beta
F_{n,2}.$

###  3.4 Re-initialization methods for the level-set function

The level-set function must be defined as a signed distance function to
achieve accurate approximations of the normal vectors to the interfaces, of
the curvature of the interfaces, and of the capillary force. This property is
lost when the interface is advected, and due to numerical diffusion. Hence,
the re-initialization of the level set function is necessary.

The typical approach for re-initialization involves finding an approximation
of the Eikonal boundary-value problem:

$\displaystyle\|\nabla\phi(\bm{x})\|=1,\quad$
$\displaystyle\forall\bm{x}\in\upOmega;$ $\displaystyle\phi(\bm{x})=0,\quad$
$\displaystyle\forall\bm{x}\in\mathrm{\Gamma}.$

The re-initialization methods fall mainly into two categories:

  * •

The Pseudo-Time Evolution technique, where the level-set function is advected
until steady-state is reached using:

$\displaystyle\partial_{\tau}\phi+\text{sgn}(F)(\|\nabla\phi\|-1)=0;$ (10)
$\displaystyle\phi(\bm{x},\tau=0)=F(\bm{x}),$ (11)

where $F(\bm{x})$ is the level-set function that needs to be re-initialized,
$\tau$ is a pseudo-time, $\text{sgn}(F)=1$ if $F$ is positive, and
$\text{sgn}(F)=-1$ if $F$ is negative.

However, re-initializing the level-set function discretely introduces errors
in the position of the interfaces, leading to mass conservation errors.
Several improvements have been introduced to alleviate this problem, such as
the constrained reinitialization method [11], where the position of the
interfaces is included in the formulation using the least-squares method.

  * •

The fast marching method is a technique that computes an approximate solution
to the Eikonal equation by marching through the grid points, starting from the
nearest points to the interface. It uses a priority queue to determine which
grid points are processed next, based on the current solution estimates. The
set of points is classified into three categories: the accepted, the
tentative, and the distant points. The key idea is to approximate the solution
for the tentative points using the approximate solution of the accepted points
with an upwind scheme, and to move the points with the lowest value of $\phi$
from the tentative to the accepted category.

In this section, we consider two re-initialization methods. The first one aims
at finding the distance from any point $\bm{x}_{0}$ in the domain to the
interface by minimizing: —l— x∥ x\- x_0 ∥; F(x)=0. Gradient-based optimization
techniques are now feasible due to the continuous representation of
$F(\bm{x})$, and by the fact that we have access to $\nabla F(\bm{x})$ by
back-propagation and automatic differentiation. The optimization is performed
based on two formulations of the problem. The first formulation, which we will
refer to as the penalty reinitialization, uses the penalty method —l— x∥
x-x_0∥ + λ∥F(x)∥, where $\lambda$ is a penalization parameter. The second
formulation, which we will refer to as the unconstrained reinitialization,
aims at removing the conflict between the terms $\|\bm{x}-\bm{x}_{0}\|$ and
$\|F(\bm{x})\|$, where we rather seek to align the vector $\bm{x}-\bm{x}_{0}$
with the normal vector $\nabla F$, and hence we minimize —l— x∥F(x)∥ + ∥1 -
∇F∥∇F∥ ⋅x-x0∥x-x0∥ ∥.

To accelerate the convergence, we can take advantage of the fact that the
position of the interfaces does not change significantly from one time step to
the next, which makes the points $\bm{x}$ obtained from the previous
reinitialization, a good initial guess for the current reinitialization.

$\bm{x}$neural network $N(\bm{x})$level-set
$F(\bm{x},1)$$\phi(\bm{x})=N(\bm{x})F(\bm{x},1)$ Figure 3: The architecture of
the PINN-R used to reinitialize the level-set function.

The third approach, which we will refer to as the PINN Reinitialization
(PINN-R), aims at approximating the solution to the Eikonal equation,

$\|\nabla\phi\|=1,$

using a PINN. The position of the interface is conserved by multiplying the
output of the PINN, which is designed to give strictly positive outputs, by
the level-set function $F(\bm{x},1)$. This is shown in figure 3, where
$F(\bm{x},1)$ is the level-set function that needs to be reinitialized, the
neural network $N(\bm{x})$ gives a strictly positive scalar output, and the
reinitialized level-set function $\phi(\bm{x})$.

##  4 Numerical results

###  4.1 Comparison of the time-integration schemes

We first investigate which of the previously described time-integration
schemes, the weak PINN, the FDPINN, and the strong PINN, is more accurate,
using a manufactured problem. The numerical test involves a cylinder with a
radius of $R=1$, with height $h=1$. Inside the cylinder, a bubble with a
radius of $r=0.25$, is positioned at the center of the cylinder, $\bm{x}_{\rm
c}=[0,0,0.5]$. The bubble is subjected to the velocity field
$\bm{u}=[-0.5x,-0.5y,z]$. With this velocity field, the exact position of any
particle, obtained by solving an ODE for each direction:

$\displaystyle\frac{dx}{dt}=-0.5x;$ $\displaystyle\frac{dy}{dt}=-0.5y;$
$\displaystyle\frac{dz}{dt}=z,$

is given by

$\bm{x}(t)=[x_{0}e^{-0.5t},y_{0}e^{-0.5t},z_{0}e^{t}],$

where $x_{0}$, $y_{0}$ and $z_{0}$ represent the particle’s initial position.
The accuracy of the solution is measured using
$E=\frac{1}{{n}}\sum_{i=1}^{n}|F(\bm{x}_{i})|,$ where $\bm{x}_{i}$ is a set of
${n}$ points that belong to the interface at time $t=0$. The chosen time step
size is ${\rm\Delta}t=0.1,$ for a total simulation time of $t=1$. The number
of collocation points is 70,000, and we train the PINN for 4000 iterations
using the Adam algorithm, with a learning rate of $\alpha=10^{-4}.$ We do not
use the multi-level approach, and use instead a feedforward neural network
with 5 layers and 50 neurons per layer.

As shown in figure 4, strongly imposing the initial conditions yields the most
accurate results in comparison with the other two formulations. This
formulation does not have a time discretization error, contrary to the FDPINN,
and it does not suffer from errors on the initial conditions, contrary to the
weak PINN.

(a)

(b)

(c)

Figure 4: (a) The interface shape obtained at $t=1$ for the comparison of the
three time-integration schemes, (b) a zoom-in on the interface shape at $t=1$,
(c) the error on the position of the interface, for the problem with a
manufactured solution using the three different time integration schemes.

###  4.2 Comparison between the multi-level PINN and the strong PINN

In this section, we compare the proposed multi-level PINN with a strong
imposition of the initial condition, with the single-level strong PINN. We use
the verification problem of subsection 3.4.1. The results are obtained with
$70,000$ collocation points, and ${\rm\Delta}t=0.1.$ We train the multi-level
PINN with ${n}=6000$ total iterations, ${n}_{1}=3000$ for the first function
$F_{n,1}$ and ${n}_{2}=3000$ for the second function $F_{n,2}.$ We train the
single-level PINN with ${n}=6000$ iterations.

(a)

(b)

(c)

Figure 5:  (a) The interface shape obtained at time $t=1$ with the strong
single-level PINN, and the proposed multi-level PINN, (b) a zoom-in on the
interface shapes, (c) the error $E$ on the interface for the proposed methods,
for the problem with a manufactured solution.

Figure 5 illustrates the capabilities of the proposed multi-level PINN, for
the same computational cost of the standard PINN formulation, where the error,
measured as ${E=\frac{1}{n}\sum_{i=1}^{n}|F(\bm{x}_{i})|}$, is equal to
${E=5\times 10^{-5}}$ for the multi-level PINN, and $E=10^{-3}$ for the
regular PINN.

###  4.3 Re-initialization methods

We apply the three reinitialization methods to a problem with a manufactured
solution to assess their effectiveness. The domain is a $2$D unit square. The
initial level-set function is given by

$F=(x+1)^{2}(r-0.25)^{3/2},$

where $r$ is the distance to the point $\bm{x}_{c}=[0.5,0.5].$ The exact
solution, after the reinitialization, is $\phi=r-0.25.$ We compare the
contours of the reinitialized level-set function obtained by the three
approaches, as shown in figure 6. We notice that the results obtained by the
penalty reinitialization are not accurate, when compared to the other
reinitialization methods. Next, we compare the curvature obtained by the
unconstrained reinitialization and the PINN-R. As seen in figure 7,
reinitializing the level-set function with the PINN-R results in a smoother
approximation, when compared to optimization-based methods.

The numerical experiment demonstrates the superiority of the PINN-R for
several key reasons. The optimization-based methods have more degrees of
freedom, equivalent to the number of collocation points. In contrast, the
PINN-R method involves a significantly lower number of degrees of freedom,
which are the weights and biases of the neural network. Consequently, the
PINN-R method is more computationally efficient. The PINN-R method conserves
the interface location during reinitialization, which leads to mass
conservation, which is not inherently guaranteed by the other methods. The
PINN-R approach leads to better curvature and capillary force approximations.

Figure 6:  The level-set function after the reinitialing
$F=(x+1)^{2}(r-0.25)^{3/2}$, using the penalty reinitialization, the
unconstrained reinitialization, and the PINN-R method.

(a)

(b)

Figure 7:  The error on the curvature of the level-set function after
reinitializing the level-set function $F=(x+1)^{2}(r-0.25)^{3/2},$ with (a)
the unconstrained reinitialization method, (b) the PINN-R method.

###  4.4 The steady-state Laplace problem

For this problem, a circular bubble with a radius of $r=0.2$ is placed in a
zero-velocity field without any external force. The viscosity, density, and
surface tension coefficient of the fluids are $\mu_{1}=10$, $\mu_{2}=0.1$,
$\rho_{1}=1000$, $\rho_{2}=1$, and $\gamma=1.96$. The bubble should remain
circular, and the magnitude of the velocity field should remain zero, with a
pressure jump across the interface. However, a velocity field appears due to
numerical errors, usually referred to as spurious currents [28, 1]. The
capillary force is balanced out by the pressure jump across the interface, as
pointed out by the Navier-Stokes equations when $\bm{u}=\bm{0}$,

$\nabla p=\kappa\gamma\delta\bm{n}.$

Hence, the pressure jump across the interface is given by
${\rm\Delta}p=\kappa\gamma$. The spurious currents appear when the numerical
pressure gradient cannot counteract the capillary force. Other factors are the
errors in the curvature of the interface approximation, errors in the normal
vectors to the interface, or a mismatch between the vector space of the
pressure gradient and capillary force in the FEM formulation [28].

In the first test, the elements size is $h=0.02,$ and the velocity field and
pressure will be discretized using the $\rm P2/\rm P1$ element combination.
The exact level-set function, given by

$F=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}-0.2,$ (12)

is used. The capillary force is modeled using the continuum surface force
model, and the partial derivatives of the level-set function are computed
using automatic differentiation. The capillary force is interpolated using
$\rm P0$ and $\rm P1$ elements to test the influence of the choice of the
vector space on the spurious currents.

(a)

(b)

Figure 8:  The spurious velocity field for the steady state Laplace problem
obtained by interpolating the capillary force (a) using the $\rm P_{1}$
element, and (b) using the $\rm P_{0}$ element.

The L2 norm of the spurious currents is $\|\bm{u}\|=10^{-4}$ for the $\rm
P_{0}$ approximation and $\|\bm{u}\|=5\times 10^{-4}$ for the $\rm P_{1}$
approximation. In both cases, the pressure jump across the interface
corresponds to the theoretical value with a relative error of $0.3\%$.

These results show that the choice of the elements used for the capillary
force interpolation is important for reducing the spurious currents magnitude,
even when the level-set function is exact.

(a)

(b)

(c)

(d)

Figure 9:  The spurious velocities in the steady-state Laplace problem for
different elements size, (a) with elements size $h=1/60,$ (b) $h=1/120,$ (c)
$h=1/180,$ and in (d) we show the $L2$ norm of the velocities.

Next, we study the effect of mesh refinement on the magnitude of the spurious
velocity. The velocity field and pressure are approximated using The $\rm
P_{2}/\rm P_{1}$ element combination. The exact level-set function (12) is
used. The capillary force is modeled using the continuum surface force model,
where the partial derivatives of the level-set function are computed using
automatic differentiation, and the surface tension is interpolated using the
$\rm P0$ element. As seen in figure 9, the element size plays an important
role in reducing the spurious currents.

Finally, we test the effect of the PINN approximation error on the spurious
currents. We train the neural network to approximate the initial condition
with the loss function

$L=\|F_{0}-(\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}-0.2)\|,$

using a different number of training iterations. The element size is
$h=1/180$, and the capillary force is interpolated using the $\rm P_{0}$
element. The number of training iterations are $n_{1}=1000,\ 2000,$ and $3000$
to obtain the first function $F_{0,1}$, and the number of training iterations
to obtain the second function $F_{0,2}$ is set to $n_{2}=2n_{1}.$

(a)

(b)

(c)

(d)

Figure 10:  The spurious velocities for the steady-state Laplace problem,
obtained using three different number of training iterations, where the first
function $F_{n,1}$ is obtained with (a) $n_{1}=1000$ training iterations, (b)
$n_{1}=2000$ training iterations,(c) $n_{1}=3000$ training iterations. In (d)
we show the $L2$ norm of the velocities for a different number of training
iterations. The second function $F_{n,2}$ is obtained by training for
$n_{2}=2n_{1}$ iterations.

The results shown in figure 10 indicate that the approximation error of the
PINN on the level-set function is an additional source for the creation of
spurious currents, that can be diminished by increasing the number of training
iterations.

###  4.5 Rising bubble problem

The next test is the 2D rising bubble problem, where we compare our results to
the benchmarked results performed in [12]. The domain geometry is rectangular
with a height of 2 and a width of 1. The initial conditions are $F_{0}=r-0.25$
for the level-set function, where $r$ is the distance to the point
$[x,y]=[0.5,0.5]$, and the velocity field $\bm{u}=\bm{0}$. The physical
parameters used in [12] can be found in Table 2. The non-dimensional numbers
are calculated with the following characteristic values: $\hat{L}=1$;
$\hat{u}=1;$ $\hat{\rho}=\rho_{1}$; $\hat{\mu}=\mu_{1}$. With these
characteristic values, we obtain that ${\rm Re}=100$, ${\rm We}=1000/24.5$,
and ${\rm Fr}=1/\sqrt{0.98}$ for the first test case, and ${\rm Re}=100$,
${\rm We}=1000/1.96$, and ${\rm Fr}=1/\sqrt{0.98}$ for the second test case.

The results documented in [12] are obtained using three different simulation
codes.

Test case | $\rho_{1}$ | $\rho_{2}$ | $\mu_{1}$ | $\mu_{2}$ | $\gamma$
---|---|---|---|---|---
1 | 1000 | 100 | 10 | 1 | 24.5
2 | 1000 | 1 | 10 | 0.1 | 1.96
Table 2: Physical parameters for the rising bubble test.

For the first case, the bubble undergoes little deformation, and the codes
tested by [12] gave the same bubble shape and rise velocity. The time step
chosen for our simulation is ${\rm\Delta}t=0.005$. We reinitialize the level-
set function at every time step, and we use 200 training iterations per time
step to obtain the first function $F_{n,1}$ and 600 training iterations per
time step to obtain the second function $F_{n,2}$, using the Adam optimizer.
From the results presented in figures 11 and 12, we notice good agreement with
the results obtained by [12] and the Comsol solver. The mass loss is around
$1.8\%,$ and the rise velocity of the bubble follows the same behavior as the
velocities found in [12].

Figure 11:  The bubble shape obtained at $t=3$s for the first rising bubble
test case problem, with the PINN-FEM coupling with elements size of $h=1/32$
and $h=1/64$, Comsol, and the results obtained in [12].

(a)

(b)

Figure 12:  (a) The bubble rise velocity, and (b) the mass loss, for the first
rising bubble test case problem, obtained with the PINN-FEM coupling with
elements size of $h=1/32$ and $h=1/64$, and the results obtained in [12].

The second case is more challenging because the bubble will break, leaving
behind thin filaments. The codes tested in [12] did not produce similar bubble
shapes, particularly in the filaments and the broken-off bubbles. In our
simulation, we use a time step of ${\rm\Delta}t=0.005$. We reinitialize the
level set function every ten time steps, and we use 200 training iterations
per time step to obtain the first function $F_{n,1}$, and 600 training
iterations per time step to obtain the second function $F_{n,2}$, using the
Adam optimizer. The results are portrayed in figure 13, where we can see that
the proposed methodology yields comparable shapes to the benchmarked results.

(a)

(b)

Figure 13: The bubble shape obtained at $t=3$s for the second rising bubble
test-case, obtained by the PINN-FEM coupling, (a) with elements size $h=1/80$,
and (b) with $h=1/160$, Comsol, and the results obtained in [12].

(a)

(b)

Figure 14:  (a)The bubble rise velocity, (b) the mass loss, for the second
rising bubble test case problem. The mass loss obtained is $2.5\%$ with
elements size $h=1/80$, and $1.2\%$ with $h=1/160$.

##  5 Conclusion

We developed a coupled PINN/FEM method for the simulation of multi-fluid
flows. The Navier-Stokes equations were discretized using the FEM, with the
P2/P1 element combination for the velocity and pressure. The level-set
function was approximated using a multi-level PINN, where the first PINN
approximates the level-set function in the entire domain, and the second PINN
refines the approximation near the interface. Three PINN formulations were
proposed to handle the imposition of the initial conditions: weak imposition,
strong imposition, and discrete-time steps. Numerical tests showed that the
strong imposition of initial conditions gives the most accurate results and
that the accuracy of the level-set function is mainly affected by the number
of training iterations per time step. Three reinitialization strategies were
proposed, that give a differentiable solution by automatic differentiation.
The first two strategies find the distance from any point in the domain to the
interface by constrained minimization. The third strategy reinitializes the
level-set function using a PINN for the Eikonal equation, with a strong
imposition of the interface position. Numerical tests showed that the
reinitialization by a PINN yields more accurate results for the curvature
calculations, and better mass-conservation properties. Finally, the proposed
numerical strategy was tested against the rising bubble problem, benchmarked
in [12]. It was found that the proposed strategy can produce results
comparable to traditional methods.

## References

  * [1] T. Abadie, J. Aubin, and D. Legendre.  On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks.  Journal of Computational Physics, 297:611–636, 2015.
  * [2] M. Alam, M. D. Samad, L. Vidyaratne, A. Glandon, and K. M. Iftekharuddin.  Survey on deep neural networks in speech and vision systems.  Neurocomputing, 417:302–321, 2020.
  * [3] Z. Aldirany, R. Cottereau, M. Laforest, and S. Prudhomme.  Multi-level neural networks for accurate solutions of boundary-value problems.  arXiv preprint arXiv:2308.11503, 2023.
  * [4] A. Beck, D. Flad, and C.-D. Munz.  Deep neural networks for data-driven les closure models.  Journal of Computational Physics, 398:108910, 2019.
  * [5] R. Bollapragada, J. Nocedal, D. Mudigere, H.-J. Shi, and P. T. P. Tang.  A progressive batching L-BFGS method for machine learning.  In International Conference on Machine Learning, pages 620–629. PMLR, 2018.
  * [6] J. U. Brackbill, D. B. Kothe, and C. Zemach.  A continuum method for modeling surface tension.  Journal of computational physics, 100(2):335–354, 1992.
  * [7] S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis.  Physics-informed neural networks (pinns) for fluid mechanics: A review.  Acta Mechanica Sinica, 37(12):1727–1738, 2021.
  * [8] S. Cai, Z. Wang, S. Wang, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks for heat transfer problems.  Journal of Heat Transfer, 143(6), 2021.
  * [9] Y. Cao, Z. Fang, Y. Wu, D.-X. Zhou, and Q. Gu.  Towards understanding the spectral bias of deep learning.  arXiv preprint arXiv:1912.01198, 2019.
  * [10] J. Chai, H. Zeng, A. Li, and E. W. Ngai.  Deep learning in computer vision: A critical review of emerging techniques and application scenarios.  Machine Learning with Applications, 6:100134, 2021.
  * [11] D. Hartmann, M. Meinke, and W. Schröder.  The constrained reinitialization equation for level set methods.  Journal of computational physics, 229(5):1514–1535, 2010.
  * [12] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska.  Quantitative benchmark computations of two-dimensional bubble dynamics.  International Journal for Numerical Methods in Fluids, 60(11):1259–1288, 2009.
  * [13] X. Jin, S. Cai, H. Li, and G. E. Karniadakis.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics, 426:109951, 2021.
  * [14] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [15] A. Kovacs, L. Exl, A. Kornell, J. Fischbacher, M. Hovorka, M. Gusenbauer, L. Breth, H. Oezelt, D. Praetorius, D. Suess, et al.  Magnetostatics and micromagnetics with physics informed neural networks.  Journal of Magnetism and Magnetic Materials, 548:168951, 2022.
  * [16] J. Lim and D. Psaltis.  Maxwellnet: Physics-driven deep neural network training based on maxwell’s equations.  Apl Photonics, 7(1):011301, 2022.
  * [17] J. Ling, A. Kurzawski, and J. Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 807:155–166, 2016.
  * [18] S. Miyazaki and Y. Hattori.  Improving accuracy of turbulence models by neural network.  arXiv preprint arXiv:2012.01723, 2020.
  * [19] A. B. Nassif, I. Shahin, I. Attili, M. Azzeh, and K. Shaalan.  Speech recognition using deep neural networks: A systematic review.  IEEE Access, 7:19143–19165, 2019.
  * [20] J. Oldenburg, F. Borowski, A. Öner, K.-P. Schmitz, and M. Stiehm.  Geometry aware physics informed neural network surrogate for solving navier–stokes equation (gapinn).  Advanced Modeling and Simulation in Engineering Sciences, 9(1):8, 2022.
  * [21] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville.  On the spectral bias of neural networks.  In International Conference on Machine Learning, pages 5301–5310. PMLR, 2019.
  * [22] M. Raissi.  Deep hidden physics models: Deep learning of nonlinear partial differential equations.  The Journal of Machine Learning Research, 19(1):932–955, 2018.
  * [23] J. Sola and J. Sevilla.  Importance of input data normalization for the application of neural networks to complex industrial problems.  IEEE Transactions on nuclear science, 44(3):1464–1468, 1997.
  * [24] N. Sukumar and A. Srivastava.  Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks.  Computer Methods in Applied Mechanics and Engineering, 389:114333, 2022.
  * [25] S. Tadeparti and V. V. Nandigana.  Convolutional neural networks for heat conduction.  Case Studies in Thermal Engineering, 38:102089, 2022.
  * [26] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. Barron, and R. Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  Advances in Neural Information Processing Systems, 33:7537–7547, 2020.
  * [27] C. Xie, X. Xiong, and J. Wang.  Artificial neural network approach for turbulence models: A local framework.  Physical Review Fluids, 6(8):084612, 2021.
  * [28] S. Zahedi, M. Kronbichler, and G. Kreiss.  Spurious currents in finite element based level set methods for two-phase flow.  International Journal for Numerical Methods in Fluids, 69(9):1433–1456, 2012.
  * [29] P. Zhang, Y. Hu, Y. Jin, S. Deng, X. Wu, and J. Chen.  A maxwell’s equations based deep learning method for time domain electromagnetic simulations.  IEEE Journal on Multiscale and Multiphysics Computational Techniques, 6:35–40, 2021.


11institutetext: Fontys University of Applied Science, Eindhoven, The
Netherlands 22institutetext: Eindhoven University of Technology, Eindhoven,
The Netherlands 33institutetext: IT University of Copenhagen, Denmark

# Dataset Distribution Impacts Model Fairness: Single vs. Multi-Task Learning

Ralf Raumanns 1122 Gerard Schouten 11 Josien P. W. Pluim 22
Veronika Cheplygina 33

###### Abstract

The influence of bias in datasets on the fairness of model predictions is a
topic of ongoing research in various fields. We evaluate the performance of
skin lesion classification using ResNet-based CNNs, focusing on patient sex
variations in training data and three different learning strategies. We
present a linear programming method for generating datasets with varying
patient sex and class labels, taking into account the correlations between
these variables. We evaluated the model performance using three different
learning strategies: a single-task model, a reinforcing multi-task model, and
an adversarial learning scheme. Our observations include: 1) sex-specific
training data yields better results, 2) single-task models exhibit sex bias,
3) the reinforcement approach does not remove sex bias, 4) the adversarial
model eliminates sex bias in cases involving only female patients, and 5)
datasets that include male patients enhance model performance for the male
subgroup, even when female patients are the majority. To generalise these
findings, in future research, we will examine more demographic attributes,
like age, and other possibly confounding factors, such as skin colour and
artefacts in the skin lesions. We make all code available on Github.

###### Keywords:

Skin lesions Bias Fairness Multi-task learning Adversarial learning

##  1 Introduction

Deep learning has shown many successes in medical image diagnosis [30, 12, 3],
but despite high overall performance, models can be biased against patients
from different demographic groups [1, 22, 15]. Bias and fairness are becoming
an active topic in medical imaging, with studies focusing for instance on skin
lesions [1, 17], chest x-rays [22] and brain MR scans [27]. Examples of
sensitive attributes include age, sex or race. For skin lesion classification,
the Fitzpatrick skin type is often studied [31, 4, 17, 34].

Fairness studies typically include baselines showing bias between groups,
and/or propose methods to improve fairness. The methods are based on sampling
or weighting strategies during training [17], and/or introducing training
strategies that try to debias the methods to rely on the sensitive attributes,
such as adversarial methods [1]. For instance, Yang and colleagues [35]
developed an adversarial debiasing framework to reduce biases in hospital
location and patient ethnicity. Similarly, Wu et al. introduced FairPrune
[34], a method for pruning parameters based on their significance to both
privileged and unprivileged groups, focusing on sex and skin tone. Moreover,
Bevan and Atapour-Abarghouei [5] employed various strategies to limit bias in
skin lesion images, specifically targeting discrepancies arising from medical
instruments, surgical markings, and rulers. Popular datasets include ISIC skin
lesion datasets [18, 10, 9, 32, 11, 29] and Fitzpatrick-17K [17, 16].
Researchers either use already provided data splits for evaluation, or split
the data ratios of patients with a specific demographic attribute, for example
male vs female patients.

Our current study builds on two crucial insights from other topics in medical
imaging: multi-task learning and shortcut learning [13, 25]. Firstly, some
studies use demographic attributes within multi-task learning settings; for
example,[23]. Here the attributes are _reinforcing_ the diagnosis during
optimization. This is at odds with the more recent adversarial strategies [2,
1] where models are encouraged to NOT predict the sensitive attribute.
Secondly, there are correlations between demographics and demographic
attributes and shortcut learning, including, for example, imaging devices and
surgical markers [33, 21, 15, 6, 5]. In such cases, simply splitting the data
according to a specific attribute can create imbalance in terms of the other
attributes, thus the observed (un)fairness could be due to the attributes that
were not considered.

Our contributions are as follows:

  1. 1.

We propose using the linear programming (LP) approach for skin lesions
retrieved from the ISIC archive [18, 10, 9, 32, 11, 29] via the gallery
browser [20]. This method gives more control over patient subset assignment.
It adjusts the proportions of selected dataset attributes while keeping others
constant.

  2. 2.

We systematically study two strategies that handle the demographic variable in
different ways: a reinforcing multi-task strategy [24, 28] and an adversarial
strategy [1, 2, 8]

  3. 3.

We evaluate our models using overall and subgroup Area Under the Curve (AUC)
based on sex, and show that:

     * •

Models perform better for male subgroups in the male-only and lightly skewed
male patient experiments. In the balanced dataset and lightly skewed female
patient experiments, there is no significant difference between the subgroups.
However, in the lightly skewed female patient scenario, the adversarial model
performs better for male patients.

     * •

Models trained exclusively on female patients exhibit a positive difference in
performance for female patients.

     * •

The base model reveals a significant sex bias, performing worse for female
patients, except when trained exclusively on female patients.

     * •

The reinforcement model has no significant effect on sex bias.

     * •

The adversarial model significantly reduces sex bias in scenarios involving
only female patients.

     * •

Eliminating model bias is challenging, with significant performance gaps
observed in datasets with skewed sex distributions.

  4. 4.

We make all data and models available on https://github.com/raumannsr/data-
fairness-impact.

##  2 Methods

Construction of datasets. We trained and validated our models on datasets with
different female (F) to male (M) patient ratios, equal numbers of malignant
and benign lesions, and equal number of patients below and above 60 (median
age) for each sex. We refer to the datasets as M100 (100% male patients),
F25M75 (25% female patients, 75% male patients), F50M50, F75M25, and F100 are
defined analogously. We evaluate the models using a balanced test set
mirroring F50M50.

We used the ISIC archive’s [18, 10, 9, 32, 11, 29] gallery browser [20], which
had 81,155 dermoscopic images of skin lesions, some with age and sex metadata.
We queried the archive for “dermoscopic” images diagnosed as “benign” or
“malignant” for all ages and both sexes. This gave us 71,035 images (62,439
benign, 8,596 malignant), which we processed using the steps in Fig. 1 (see
Appendix for more details).

pt1\. Remove all lesions with an undefined age attribute.pt2\. Remove
duplicates following Cassidy’s approach [7].pt3\. Remove multiplets by
selecting a random image from each patient.pt4\. Use linear programming to
find the optimal sampling accordingto sex, age (above or below median), and
class label.pt5\. Reserve 12.5% (1264 lesions) of the balanced dataset as the
test set.pt6\. Subsample the remaining data into subsets of 2206 malignant and
2206 benign instances each, see Table 1.pt7\. Use the subsets to create
training and validation sets with3528 and 884 lesions respectively. Figure 1:
Steps for filtering lesions and creating test, training and validation sets.
Steps 5 through 7 are repeated using 5 different seeds in a cross-validation
setup.

Linear programming for optimal dataset construction. We have developed a
method to create diverse dataset compositions using linear programming, a
common mathematical optimisation technique. The goal is to maximise the number
of instances of skin lesions within defined constraints, as we express below:

$\begin{array}[]{ll@{}ll}\text{Find a vector}&&x&\text{(decision
variables)}\\\ \text{that maximises}&&f=x_{1}&\text{(objective\ function)}\\\
\text{subject to}&&a_{i1}x_{1}+a_{i2}x_{2}+\dots+a_{in}x_{n}\leq
b_{i}&\text{(constraints)}\\\ &&\text{for }i=1,\dots,13&\\\
\text{and}&&x_{j}\geq 0.&\text{(non-negativity constraints)}\\\ &&\text{for
}j=1,\dots,14&\\\ \end{array}$

  * •

The decision variables ($x_{1},\dots,x_{14}$) denote specific categories, like
benign lesions in female patients aged 60 and above. See Appendix for more.

  * •

The objective function of the LP model is designed to maximise the number of
malignant instances $x_{1}$. There are fewer malignant instances than benign
ones in the ISIC archive, and the goal is to achieve a balance between the
two.

  * •

The constraints limit the solution by setting specific limits for each group
and maintaining ratios between these groups. Group examples include all benign
lesions, all female patients over 60, and all male patients under 60. The
primary constraint ensures an equal number of malignant and benign lesions
($x_{1}-x_{2}=0$). See Appendix for more.

  * •

Due to non-negativity constraints, decision variables cannot be negative.

Table 1: Datasets are distributed amongst malignant, benign, male patients
(M), and female patients (F) categories for both training and validation.

| M100 | F25M75 | F50M50 | F75M25 | F100
---|---|---|---|---|---
Malignant (M/F) | 2206 (2206/0) | 2941 (2206/735) | 4412 (2206/2206) | 3235 (809/2426) | 2426 (0/2426)
Benign (M/F) | 2206 (2206/0) | 2941 (2206/735) | 4412 (2206/2206) | 3235 (809/2426) | 2426 (0/2426)

Within set constraints, the optimal solution maximises malignant lesions and
assigns value to decision variables. To find this solution, we created a
unique LP model for each dataset. Table 1 shows the result of the LP model for
the different datasets.

Models. We used the ResNet50 model [19] in three ways, which include:

  * •

The single-task baseline model enhanced with two fully connected layers has a
sigmoid activation function and binary cross-entropy loss ($L_{c}$, see
Equation 1). We did not use class weights but rather the actual distribution
represented by the training dataset. We fine-tuned the model through a grid
search of three varying learning rates and batch sizes, selecting the
combination that yielded the highest performance across all experiments.

  * •

The multi-task reinforcing model, with three added layers to the convolutional
base, produces two outputs: one for classification and the other for the
binary sex attribute. We employed binary cross-entropy loss ($L_{c}$) and a
sigmoid activation function for both heads, giving equal weight to both
losses.

  * •

The multi-task adversarial model was implemented following the methodology [2,
1], using a network with a shared feature encoder and two classifier heads.
One classifier targeted skin cancer classification; the other predicted
confounding variables like sex or age. We used the ResNet architecture to
compare performance with baseline and reinforcing models equitably. We trained
the skin cancer classifier and encoder with cross-entropy loss ($L_{c}$) and
optimised the bias predictor with binary cross-entropy loss ($L_{c}$). To
diminish the confounder predictiveness of the encoded feature, we
adversarially adjusted the encoder using a third loss ($L_{br}$), setting λ as
the penalty for accurate demographic predictions. We used a lambda (λ, see
Equation 2) value of 5, as in [1], to assess subgroup performance and set the
penalty for correct predictions of the target demographic parameter.

To summarise, we use these loss functions:

$L_{c}=-\sum_{i=1}^{n}[y_{i}\log(\hat{y_{i}})+(1-y_{i})\log(1-\hat{y_{i}})]$
(1) $L_{br}=\lambda{L_{c}}$ (2)

where $n$ represents the number of lesions, and $y_{i}$ and $\hat{y_{i}}$
denote respectively the prediction and the expected outcome for lesion $i$.

We pre-trained all networks with ImageNet and resized images to 384 × 384 for
ResNet50’s input size. During model training, we used data augmentation
techniques, ran up to 40 epochs with a batch size of 20, and set a learning
rate of $2.0\mathrm{e}{-5}$. We stopped training if no significant improvement
occurred after 10 consecutive epochs to avoid overfitting. We implemented our
baseline and reinforcing models in Keras with the TensorFlow backend [14] and
our adversarial model in PyTorch [26].

Evaluation. For the purpose of an in-depth evaluation, we generated five
distinct instances of each dataset: M100, F25M75, F50M50, F75M25 and F100.
Each instance was created with a unique seed to ensure diversity and
robustness in our evaluation. Each seed corresponds to a balanced test set to
allow fair comparisons between dataset-model instances. Furthermore by using a
balanced test set, we ensure that the results are not skewed towards any
specific subgroup. We evaluate AUC overall and for male and female subgroups
for each model architecture and dataset combination.

Figure 2: The AUC score varies based on data splits ranging from only male
patients (M100) to only female patients (F100). We show base, reinforcing and
adversarial model performance for female and male patient subgroups.
Significance per Mann–Whitney U test (as used in [22]) is denoted by ****
$(P\leq 0.0001)$, *** $(0.0001<P\leq 0.001)$, ** $(0.001<P\leq 0.01)$, *
$(0.01<P\leq 0.1)$, and not significant (ns) $(P>0.1)$. < indicates lower
AUCs, > higher AUCs, and = comparable AUCs for female patients.

##  3 Results

Figure 2 shows the impact of dataset distributions on three learning
strategies, reporting AUC scores overall and for both sexes.

Sex-specific training data yields better results. Models perform better for
male subgroups in the male-only and lightly skewed male patient experiments.
In the balanced dataset and lightly skewed female patient experiments, there
is no significant difference between the subgroups. However, in the lightly
skewed female patient scenario, the adversarial model performs better for male
patients. An exception is observed when the training datasets consist only of
female patients; in such cases, there is a positive difference in AUC scores
favouring female patients. Thus, our models seem more attuned to male patients
in mixed-sex training sets, irrespective of the percentage of female patients.
The best results are achieved when both sexes are trained exclusively on their
respective data. Despite this, a male subgroup bias is apparent as the results
for female patients are significantly worse than for male patients when
trained exclusively on their data.

Base model reveals sex bias. The base model shows a significant sex bias in
performance. When only male patients are involved the base model reveals a
substantial performance gap between male and female patients. The results are
significantly worse for female patients in male-skewed scenarios, except for
the female-only experiments. Interestingly, the base model performs better for
female patients than the adversarial and reinforcing models in the F75M25
dataset scenario. It is worth noting that the performance gap between the
subgroups is not as apparent when the dataset includes male and female
patients, unlike in the experiment that only involved male patients.

Reinforcement model has no significant effect on sex bias. Training the model
only on male patients increases AUC score variability and reduces performance
differences between both sexes. The reinforcement model does not significantly
affect sex bias.

Adversarial model eliminates sex bias in cases involving only female patients.
The adversarial model reduces sex bias significantly in scenarios with only
female patients but is less effective in other scenarios, often favouring male
patients. Its performance varies across experiments and datasets.

##  4 Discussion and conclusions

We studied model and subgroup performance across datasets to identify the
influence of bias in datasets on fairness of model predictions. We used linear
programming (LP) to create various datasets with controlled male-female
ratios. This was done to systematically evaluate the performance of three
different learning strategies using ResNet-based CNNs. Other fairness and bias
studies that require a flexible method to create datasets under certain
constraints could benefit from this universal LP technique.

Our study shows that eliminating bias is challenging. The adversarial model
architecture is able to reduce sex bias in a female-only context but fails for
other datasets. Other model approaches do not show convincing results with
respect to bias reduction.

Skewed sex distributions still show a performance gap between male and female
patients. Our experiments demonstrate that the adversarial model better
corrects sex bias in female-only datasets and not in male-only datasets,
possibly due to other confounding and/or unidentified factors. Further
research is needed on this issue.

As expected the base model shows a sensitivity for sex bias, possibly due to
overfitting. The reinforcing and adversarial models both having a form of
regularisation (to counter overfitting), are potentially able to reduce sex
bias compared to the base model. However in our experiments we only see a bias
correction for adversarial models for female-only experiments.

Our outcomes show that sex-related information influences prediction tasks.
Future research should determine which specific sex-related factors are
essential to ensure fairness across different subgroups.

In contrast to categorical data like patient sex, where the groups are clearly
defined, this is not possible or only partially possible with continuous data
like age, which could lead to somewhat arbitrary subgroups. Therefore, we
started with the demographic attribute sex and will continue similar research
with the non-categorical age attribute.

Further, we have identified the following directions for future work:

  * •

Investigate whether using “early stopping” per task in a multi-task model
reduces subgroup bias.

  * •

Explore the impact of integrating segmentation with a classifier on sex-based
disparities in identifying skin lesions.

  * •

Study the roles of factors like skin colour and image artefacts in model
fairness for different subgroups.

  * •

Investigate the impact of shortcut learning on model fairness.

In conclusion, while we progress towards fairness, further advancements are
needed to ensure consistent and equitable performance across various data
distributions.

###  4.0.1 Acknowledgements

We gratefully acknowledge financial support from the Netherlands Organization
for Scientific Research (NWO), grant no. 023.014.010.

###  4.0.2

The authors declare that they have no known competing financial interests or
personal relationships that could have influenced the work reported in this
paper.

## References

  * [1] Abbasi-Sureshjani, S., Raumanns, R., Michels, B.E., Schouten, G., Cheplygina, V.: Risk of training diagnostic algorithms on data with demographic bias. In: MICCAI LABELS workshop, Lecture Notes in Computer Science. vol. 12446, pp. 183–192. Springer (2020)
  * [2] Adeli, E., Zhao, Q., Pfefferbaum, A., Sullivan, E.V., Fei-Fei, L., Niebles, J.C., Pohl, K.M.: Representation learning with statistical independence to mitigate bias. IEEE Winter Conf Appl Comput Vis 2021, 2512–2522 (Jan 2021)
  * [3] Bejnordi, B.E., Veta, M., van Diest, P.J., van Ginneken, B., Karssemeijer, N., Litjens, G., van der Laak, J.A., Hermsen, M., Manson, Q.F., Balkenhol, M., et al.: Diagnostic assessment of deep learning algorithms for detection of lymph node metastases in women with breast cancer. JAMA 318(22), 2199–2210 (2017)
  * [4] Benčević, M., Habijan, M., Galić, I., Babin, D., Pižurica, A.: Understanding skin color bias in deep learning-based skin lesion segmentation. Comput. Methods Programs Biomed. 245, 108044 (Mar 2024)
  * [5] Bevan, P.J., Atapour-Abarghouei, A.: Skin deep unlearning: Artefact and instrument debiasing in the context of melanoma classification. arXiv preprint arXiv:2109.09818 (Apr 2023)
  * [6] Bissoto, A., Valle, E., Avila, S.: Debiasing skin lesion datasets and models? not so fast (2020)
  * [7] Cassidy, B., Kendrick, C., Brodzicki, A., Jaworek-Korjakowska, J., Yap, M.H.: Analysis of the ISIC image datasets: Usage, benchmarks and recommendations. Med. Image Anal. 75, 102305 (Jan 2022)
  * [8] Chu, Z., Rathbun, S.L., Li, S.: Multi-Task adversarial learning for treatment effect estimation in basket trials. In: Flores, G., Chen, G.H., Pollard, T., Ho, J.C., Naumann, T. (eds.) Proceedings of the Conference on Health, Inference, and Learning. Proceedings of Machine Learning Research, vol. 174, pp. 79–91. PMLR (2022)
  * [9] Codella, N., Rotemberg, V., Tschandl, P., Celebi, M.E., Dusza, S., Gutman, D., Helba, B., Kalloo, A., Liopyris, K., Marchetti, M., Kittler, H., Halpern, A.: Skin lesion analysis toward melanoma detection 2018: A challenge hosted by the international skin imaging collaboration (isic) (2019)
  * [10] Codella, N.C.F., Gutman, D., Celebi, M.E., Helba, B., Marchetti, M.A., Dusza, S.W., Kalloo, A., Liopyris, K., Mishra, N., Kittler, H., Halpern, A.: Skin lesion analysis toward melanoma detection: A challenge at the 2017 international symposium on biomedical imaging (isbi), hosted by the international skin imaging collaboration (isic) (2018)
  * [11] Combalia, M., Codella, N.C.F., Rotemberg, V., Helba, B., Vilaplana, V., Reiter, O., Carrera, C., Barreiro, A., Halpern, A.C., Puig, S., Malvehy, J.: Bcn20000: Dermoscopic lesions in the wild (2019)
  * [12] Esteva, A., Kuprel, B., Novoa, R.A., Ko, J., Swetter, S.M., Blau, H.M., Thrun, S.: Dermatologist-level classification of skin cancer with deep neural networks. Nature 542(7639), 115–118 (Feb 2017)
  * [13] Geirhos, R., Jacobsen, J.H., Michaelis, C., Zemel, R., Brendel, W., Bethge, M., Wichmann, F.A.: Shortcut learning in deep neural networks. Nature Machine Intelligence 2(11), 665–673 (Nov 2020)
  * [14] Géron, A.: Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow. “O’Reilly Media, Inc.” (Oct 2022)
  * [15] Gichoya, J.W., Banerjee, I., Bhimireddy, A.R., Burns, J.L., Celi, L.A., Chen, L.C., Correa, R., Dullerud, N., Ghassemi, M., Huang, S.C., et al.: AI recognition of patient race in medical imaging: a modelling study. The Lancet Digital Health 4(6), e406–e414 (2022)
  * [16] Groh, M., Harris, C., Daneshjou, R., Badri, O., Koochek, A.: Towards transparency in dermatology image datasets with skin tone annotations by experts, crowds, and an algorithm. Proc. ACM Hum.-Comput. Interact. 6(CSCW2), 1–26 (Nov 2022)
  * [17] Groh, M., Harris, C., Soenksen, L., Lau, F., Han, R., Kim, A., Koochek, A., Badri, O.: Evaluating deep neural networks trained on clinical images in dermatology with the fitzpatrick 17k dataset. 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW) pp. 1820–1828 (Apr 2021)
  * [18] Gutman, D., Codella, N.C.F., Celebi, E., Helba, B., Marchetti, M., Mishra, N., Halpern, A.: Skin lesion analysis toward melanoma detection: A challenge at the international symposium on biomedical imaging (isbi) 2016, hosted by the international skin imaging collaboration (isic) (2016)
  * [19] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778 (2016)
  * [20] ISIC archive. https://gallery.isic-archive.com, accessed: 2024-06-07
  * [21] Jiménez-Sánchez, A., Juodelyte, D., Chamberlain, B., Cheplygina, V.: Detecting shortcuts in medical images-a case study in chest x-rays. In: 2023 IEEE 20th International Symposium on Biomedical Imaging (ISBI). pp. 1–5. IEEE (2023)
  * [22] Larrazabal, A.J., Nieto, N., Peterson, V., Milone, D.H., Ferrante, E.: Gender imbalance in medical imaging datasets produces biased classifiers for computer-aided diagnosis. Proceedings of the National Academy of Sciences 117(23), 12592–12594 (2020)
  * [23] Liu, X., Shi, J., Zhou, S., Lu, M.: An iterated Laplacian based semi-supervised dimensionality reduction for classification of breast cancer on ultrasound images. In: International Conference of the IEEE Engineering in Medicine and Biology Society. vol. 2014, pp. 4679–4682 (2014). https://doi.org/10.1109/EMBC.2014.6944668
  * [24] Marques, S., Schiavo, F., Ferreira, C.A., Pedrosa, J., Cunha, A., Campilho, A.: A multi-task CNN approach for lung nodule malignancy classification and characterization. Expert Syst. Appl. 184, 115469 (Dec 2021)
  * [25] Nauta, M., Walsh, R., Dubowski, A., Seifert, C.: Uncovering and correcting shortcut learning in machine learning models for skin cancer diagnosis. Diagnostics (Basel) 12(1) (Dec 2021)
  * [26] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J.P., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Others: An imperative style, high-performance deep learning library. Adv. Neural Inf. Process. Syst. 32, 8026–8037 (2019)
  * [27] Petersen, E., Feragen, A., da Costa Zemsch, M.L., Henriksen, A., Wiese Christensen, O.E., Ganz, M., Initiative, A.D.N.: Feature robustness and sex differences in medical imaging: a case study in mri-based alzheimer’s disease detection. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 88–98. Springer (2022)
  * [28] Raumanns, R., Schouten, G., Joosten, M., Pluim, J.P.W., Cheplygina, V.: Enhance (enriching health data by annotations of crowd and experts): A case study for skin lesion classification. Machine Learning for Biomedical Imaging 1, 1–26 (2021). https://doi.org/https://doi.org/10.59275/j.melba.2021-geb9, https://melba-journal.org/2021:020
  * [29] Rotemberg, V., Kurtansky, N., Betz-Stablein, B., Caffery, L., Chousakos, E., Codella, N., Combalia, M., Dusza, S., Guitera, P., Gutman, D., Halpern, A., Helba, B., Kittler, H., Kose, K., Langer, S., Lioprys, K., Malvehy, J., Musthaq, S., Nanda, J., Reiter, O., Shih, G., Stratigos, A., Tschandl, P., Weber, J., Soyer, H.: A patient-centric dataset of images and metadata for identifying melanomas using clinical context. Scientific Data; London 8(1), s41597–021 (2021)
  * [30] Saha, A., Bosma, J., Twilt, J., van Ginneken, B., Bjartell, A., Padhani, A., Bonekamp, D., Villeirs, G., Salomon, G., Giannarini, G., Kalpathy-Cramer, J., Barentsz, J., Maier-Hein, K., Rusu, M., Rouvière, O., van den Bergh, R., Panebianco, V., Kasivisvanathan, V., Obuchowski, N., Yakar, D., et al.: Artificial intelligence and radiologists in prostate cancer detection on mri (pi-cai): an international, paired, non-inferiority, confirmatory study. Lancet Oncol. (Jun 2024)
  * [31] Seth, P., Pai, A.K.: Does the fairness of your Pre-Training hold up? examining the influence of Pre-Training techniques on skin tone bias in skin lesion classification. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 570–577 (2024)
  * [32] Tschandl, P., Rosendahl, C., Kittler, H.: The HAM10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions (2018)
  * [33] Willemink, M.J., Koszek, W.A., Hardell, C., Wu, J., Fleischmann, D., Harvey, H., Folio, L.R., Summers, R.M., Rubin, D.L., Lungren, M.P.: Preparing medical imaging data for machine learning. Radiology p. 192224 (2020)
  * [34] Wu, Y., Zeng, D., Xu, X., Shi, Y., Hu, J.: FairPrune: Achieving fairness through pruning for dermatological disease diagnosis. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2022. pp. 743–753. Springer Nature Switzerland (2022)
  * [35] Yang, J., Soltan, A.A.S., Eyre, D.W., Yang, Y., Clifton, D.A.: An adversarial training framework for mitigating algorithmic biases in clinical machine learning. NPJ Digit Med 6(1), 55 (Mar 2023)

Upper Boundaries LP Model

Table 2: The distribution of instances after the final step, which involves removing multiplets, establishes the upper boundaries of the decision variables in the LP model. | Aged | Female subgroup | Male subgroup
---|---|---|---
Malignant | $<60$ | 1,371 | 1,261
| $\geq 60$ | 1,641 | 2,801
Benign | $<60$ | 10,810 | 12,239
| $\geq 60$ | 2,397 | 3,364

Decision variables LP Model

$x_{1}$:# malignant instances
$x_{2}$:# benign instances
$x_{3}$:# male patients (M) with malignant lesions
$x_{4}$:# female patients (F) with malignant lesions
$x_{5}$:# benign M
$x_{6}$:# benign F
$x_{7}$:# malignant lesions of M (age $<60$)
$x_{8}$:# malignant lesions of M (age $\geq 60$)
$x_{9}$:# malignant lesions of F (age $<60$)
$x_{10}$:# malignant lesions of F (age $\geq 60$)
$x_{11}$:# benign M (age $<60$)
$x_{12}$:# benign M (age $\geq 60$)
$x_{13}$:# benign F (age $<60$)
$x_{14}$:# benign F (age $\geq 60$)

Constraints LP Model

$x_{1}-x_{2}=0$ (1) $rx_{4}-x_{3}=0$ (2) $sx_{8}-x_{7}=0$ (3)
$tx_{10}-x_{9}=0$ (4) $x_{1}-x_{3}-x_{4}=0$ (5) $x_{3}-x_{7}-x_{8}=0$ (6)
$x_{4}-x_{9}-x_{10}=0$ (7) $ux_{12}-x_{11}=0$ (8) $vx_{14}-x_{13}=0$ (9)
$x_{2}-x_{5}-x_{6}=0$ (10) $x_{6}-x_{13}-x_{14}=0$ (11)
$x_{5}-x_{11}-x_{12}=0$ (12) $wx_{6}-x_{5}=0$ (13)

Eq. 1: # of malignant records equals # benign records.
Eq. 2: The ratio $r$ of malignant male patients (M) to female patients (F).
Eq. 3: The ratio $s$ of malignant M (age $<60$) to M (age $\geq 60$).
Eq. 4: The ratio $t$ of malignant F (age $<60$) to F (age $\geq 60$).
Eq. 5: # of malignant lesions equals the # M and F with malignant lesions.
Eq. 6: # M with malignant lesions is equal to that of all ages.
Eq. 7: # F with malignant lesions is equal to that of all ages.
Eq. 8: The ratio $u$ of benign M (age $<60$) to M (age $\geq 60$).
Eq. 9: The ratio $v$ of benign F (age $<60$) to F (age $\geq 60$).
Eq. 10: # benign lesions equals the # M and F with benign lesions.
Eq. 11: # F with benign lesions is equal to # all ages.
Eq. 12: # M with benign lesions is equal to # all ages.
Eq. 13: The ratio $w$ of benign M to F.