# Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks

John Hanna, José Aguado, Sebastien Comas-Cardona, Ramzi Askri, Domenico Borzacchiello

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

## To cite this version:

John Hanna, José Aguado, Sebastien Comas-Cardona, Ramzi Askri, Domenico Borzacchiello. Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks. Computer Methods in Applied Mechanics and Engineering, 2022, 396, pp.115100. 10.1016/j.cma.2022.115100. hal-04586143

HAL Id: hal-04586143

https://hal.science/hal-04586143

Submitted on 28 May 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks

John Hanna, José Aguado, Sebastien Comas-Cardona, Ramzi Askri, Domenico Borzacchiello

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

## To cite this version:

John Hanna, José Aguado, Sebastien Comas-Cardona, Ramzi Askri, Domenico Borzacchiello. Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks. Computer Methods in Applied Mechanics and Engineering, 2022, 396, pp.115100. 10.1016/j.cma.2022.115100. hal-04586143

HAL Id: hal-04586143

https://hal.science/hal-04586143v1

Submitted on 28 May 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks

John M. Hanna<sup>a,b,*</sup>, José V. Aguado<sup>a</sup>, Sebastien Comas-Cardona<sup>a</sup>, Ramzi Askri<sup>b</sup>, Domenico Borzacchiello<sup>a</sup>

Received 27 September 2021; received in revised form 3 February 2022; accepted 30 April 2022

Available online xxxx

## Abstract

This paper aims to provide a machine learning framework to simulate two-phase flow in porous media. The proposed algorithm is based on Physics-informed neural networks (PINN). A novel residual-based adaptive PINN is developed and compared with the residual-based adaptive refinement (RAR) method and with PINN with fixed collocation points. The proposed algorithm is expected to have great potential to be applied to different fields where adaptivity is needed. In this paper, we focus on the two-phase flow in porous media problem. We provide two numerical examples to show the effectiveness of the new algorithm. It is found that adaptivity is essential to capture moving flow fronts. We show how the results obtained through this approach are more accurate than using RAR method or PINN with fixed collocation points, while having a comparable computational cost.

© 2022 Elsevier B.V. All rights reserved.

**Keywords:** Physics-informed neural networks; Adaptivity; Two-phase flow

## 1. Introduction

Multi-phase flow and transport phenomena in porous media appear in a variety of industrial applications, such as injection-based fiber-reinforced composites processes [1,2], exploitation of oil reservoirs [3,4] and water resources management [5,6]. That is why efficient stable modeling techniques are always required to understand the underlying physics of flow mechanisms.

Many numerical techniques have been developed to compute solutions to the partial differential equations governing these models. Grid-based classical methods such as finite element, finite volume, and finite difference are generally seen as the state of the art due to their numerical efficiency and stability [7]. To capture the flow front movement across the computational grid, these techniques are coupled with methods such as level-set [8], volume of fluids [9] or phase-field method [10]. On the other hand, meshless and particle-based methods [11] have been proved to provide more natural ways of tracking the interface although they are generally regarded as having

<footnote><sup>a</sup> Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000, 1 Rue de la Noë, 44300 Nantes, France</footnote>

<footnote><sup>b</sup> Nantes Université, IRT Jules Verne, F-44000 Nantes, 44340 Bouguenais, France</footnote>

<footnote>* Corresponding author at: Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000, 1 Rue de la Noë, 44300 Nantes, France.</footnote>

<footnote>E-mail address: john.hanna@ec-nantes.fr (J.M. Hanna).</footnote>
# Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks

John M. Hanna<sup>a,b,*</sup>, José V. Aguado<sup>a</sup>, Sebastien Comas-Cardona<sup>a</sup>, Ramzi Askri<sup>b</sup>, Domenico Borzacchiello<sup>a</sup>

Received 27 September 2021; received in revised form 3 February 2022; accepted 30 April 2022

Available online xxxx

## Abstract

This paper aims to provide a machine learning framework to simulate two-phase flow in porous media. The proposed algorithm is based on Physics-informed neural networks (PINN). A novel residual-based adaptive PINN is developed and compared with the residual-based adaptive refinement (RAR) method and with PINN with fixed collocation points. The proposed algorithm is expected to have great potential to be applied to different fields where adaptivity is needed. In this paper, we focus on the two-phase flow in porous media problem. We provide two numerical examples to show the effectiveness of the new algorithm. It is found that adaptivity is essential to capture moving flow fronts. We show how the results obtained through this approach are more accurate than using RAR method or PINN with fixed collocation points, while having a comparable computational cost.

© 2022 Elsevier B.V. All rights reserved.

**Keywords:** Physics-informed neural networks; Adaptivity; Two-phase flow

## 1. Introduction

Multi-phase flow and transport phenomena in porous media appear in a variety of industrial applications, such as injection-based fiber-reinforced composites processes [1,2], exploitation of oil reservoirs [3,4] and water resources management [5,6]. That is why efficient stable modeling techniques are always required to understand the underlying physics of flow mechanisms.

Many numerical techniques have been developed to compute solutions to the partial differential equations governing these models. Grid-based classical methods such as finite element, finite volume, and finite difference are generally seen as the state of the art due to their numerical efficiency and stability [7]. To capture the flow front movement across the computational grid, these techniques are coupled with methods such as level-set [8], volume of fluids [9] or phase-field method [10]. On the other hand, meshless and particle-based methods [11] have been proved to provide more natural ways of tracking the interface although they are generally regarded as having

<footnote><sup>a</sup> Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000, 1 Rue de la Noë, 44300 Nantes, France</footnote>

<footnote><sup>b</sup> Nantes Université, IRT Jules Verne, F-44000 Nantes, 44340 Bouguenais, France</footnote>

<footnote>* Corresponding author at: Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000, 1 Rue de la Noë, 44300 Nantes, France.</footnote>

<footnote>E-mail address: john.hanna@ec-nantes.fr (J.M. Hanna).</footnote>
a lower accuracy and reduced stability [12,13]. Practical applications include flow monitoring and control [14], identification of flow parameters such as porosity and permeability [15], uncertainty quantification [16] and flow optimization [17].

A successful approach for solving inverse problems relies on building surrogate parametric models for fast exploration of the state space of parameters [18,19]. These meta-models are specifically designed to be an optimal trade-off between accuracy and computational cost and can be effectively used for inverse problems when the minimization of the objective function requires multiple calls to the flow solver to perform parametric sweeps. Even so, when the dimension of the parametric space is high, the resulting complexity grows exponentially making some problems computationally intractable.

In recent years model order reduction techniques have been developed to tame this curse of dimensionality [20,21]. These rely on the fundamental assumption that the solution of the parametric problem lies in a low-dimensional manifold of the original subspace where the solution approximation is sought. Learning the structure of this manifold is done through an offline training procedure minimizing the L2 distance from existing data, usually collected from multiple runs of a high-fidelity solver, and produces a low-rank basis that can be reused to represent the solution of new problems for unseen choices of the parameters. In practice, the choice of reduced-bases representation is equivalent to assuming a tensor format for the problem solution. Among the different choices the Canonical Polyadic (CP), Tucker, and Tensor Train (TT) are the most commonly used, as they provide a compact representation of the parametric solution as well as a reduced complexity of the model [22–24].

The presence of a moving flow front in multi-phase flow introduces an additional difficulty to get effective model representations. Tensor formats are generally regarded as unfit to represent solutions exhibiting a moving discontinuity because, due to the dual-scale nature of the problem, there is a need for a large number of basis vectors to obtain a good approximation of the solution. For instance, applying a standard CP space-time decomposition to a simple 1D moving Heaviside function cannot provide an accurate approximation unless a high number of modes is used. In practice, the rank needed to obtain reasonably accurate results is not offering any computational advantage compared to full-order representations. This issue is well documented in the community of model order reduction and it affects not only hyperbolic equations that are likely to give rise to shocks or discontinuities but also to other problems in which the physics involved produces localized effects in the solutions. The issue was tackled by [25–27] who proposed ways to fix the problem, however, it is still a pressing difficulty that requires attention.

In recent years, with the ever-increasing number of machine learning applications that have emerged, neural networks (NN) have established themselves as a class of architectures offering superior modeling capabilities in problems such as regression and classification. When compared to the previously mentioned tensor formats NNs might provide a more flexible tool to represent functions with moving discontinuities resulting from the solution of transport equations.

It is the main aim of this paper to assess this opportunity and provide the basis for a proper framework to solve multi-phase flow in porous media through the use of NNs.

The proposed numerical technique is based on Physics-informed neural networks (PINN), a PDE solver based on the use of neural networks approximation as the search space for the PDE solution. A solution to a PDE can be obtained by minimizing the PDE and boundary conditions residuals over a finite set of points called collocation and training points, respectively. PINN was first introduced in [28] and has been used to solve many forward and inverse problems [29–39]. We show in this paper that the choice of collocations points is crucial to obtain accurate results.

The main contributions of this paper can be summarized as:

* Providing a basic framework to simulate two-phase flow in porous media using PINN.

* Developing a residual-based adaptive PINN for accurate flow front predictions.

* Extending the adaptivity algorithm for the training points to better capture the initial/boundary conditions and to have different collocation points for different PDEs in a coupled system.

The rest of the paper is organized as follows. In Section 2, the PINN method to solve two-phase flow in porous media is described. The novel residual-based adaptivity algorithm is detailed in Section 3. Numerical experiments were carried out and presented in Section 4, showing the advantages of using the adaptivity technique introduced in this work compared to using classical PINN with fixed collocation points or the RAR technique. Finally, a summary and a conclusion are given in Section 5.
a lower accuracy and reduced stability [12,13]. Practical applications include flow monitoring and control [14], identification of flow parameters such as porosity and permeability [15], uncertainty quantification [16] and flow optimization [17].

A successful approach for solving inverse problems relies on building surrogate parametric models for fast exploration of the state space of parameters [18,19]. These meta-models are specifically designed to be an optimal trade-off between accuracy and computational cost and can be effectively used for inverse problems when the minimization of the objective function requires multiple calls to the flow solver to perform parametric sweeps. Even so, when the dimension of the parametric space is high, the resulting complexity grows exponentially making some problems computationally intractable.

In recent years model order reduction techniques have been developed to tame this curse of dimensionality [20,21]. These rely on the fundamental assumption that the solution of the parametric problem lies in a low-dimensional manifold of the original subspace where the solution approximation is sought. Learning the structure of this manifold is done through an offline training procedure minimizing the L2 distance from existing data, usually collected from multiple runs of a high-fidelity solver, and produces a low-rank basis that can be reused to represent the solution of new problems for unseen choices of the parameters. In practice, the choice of reduced-bases representation is equivalent to assuming a tensor format for the problem solution. Among the different choices the Canonical Polyadic (CP), Tucker, and Tensor Train (TT) are the most commonly used, as they provide a compact representation of the parametric solution as well as a reduced complexity of the model [22–24].

The presence of a moving flow front in multi-phase flow introduces an additional difficulty to get effective model representations. Tensor formats are generally regarded as unfit to represent solutions exhibiting a moving discontinuity because, due to the dual-scale nature of the problem, there is a need for a large number of basis vectors to obtain a good approximation of the solution. For instance, applying a standard CP space-time decomposition to a simple 1D moving Heaviside function cannot provide an accurate approximation unless a high number of modes is used. In practice, the rank needed to obtain reasonably accurate results is not offering any computational advantage compared to full-order representations. This issue is well documented in the community of model order reduction and it affects not only hyperbolic equations that are likely to give rise to shocks or discontinuities but also to other problems in which the physics involved produces localized effects in the solutions. The issue was tackled by [25–27] who proposed ways to fix the problem, however, it is still a pressing difficulty that requires attention.

In recent years, with the ever-increasing number of machine learning applications that have emerged, neural networks (NN) have established themselves as a class of architectures offering superior modeling capabilities in problems such as regression and classification. When compared to the previously mentioned tensor formats NNs might provide a more flexible tool to represent functions with moving discontinuities resulting from the solution of transport equations.

It is the main aim of this paper to assess this opportunity and provide the basis for a proper framework to solve multi-phase flow in porous media through the use of NNs.

The proposed numerical technique is based on Physics-informed neural networks (PINN), a PDE solver based on the use of neural networks approximation as the search space for the PDE solution. A solution to a PDE can be obtained by minimizing the PDE and boundary conditions residuals over a finite set of points called collocation and training points, respectively. PINN was first introduced in [28] and has been used to solve many forward and inverse problems [29–39]. We show in this paper that the choice of collocations points is crucial to obtain accurate results.

The main contributions of this paper can be summarized as:

* Providing a basic framework to simulate two-phase flow in porous media using PINN.

* Developing a residual-based adaptive PINN for accurate flow front predictions.

* Extending the adaptivity algorithm for the training points to better capture the initial/boundary conditions and to have different collocation points for different PDEs in a coupled system.

The rest of the paper is organized as follows. In Section 2, the PINN method to solve two-phase flow in porous media is described. The novel residual-based adaptivity algorithm is detailed in Section 3. Numerical experiments were carried out and presented in Section 4, showing the advantages of using the adaptivity technique introduced in this work compared to using classical PINN with fixed collocation points or the RAR technique. Finally, a summary and a conclusion are given in Section 5.
## 2. Model problem and PINN formulation

In this section, we first introduce the model problem for two-phase flow in porous media. The basics for PINN are then recalled. Finally, PINN is applied to solve the model problem.

### 2.1. Model problem for two-phase flow in porous media

Flow in porous media can be described by Darcy's law, which reads as

$$ \mathbf{v} = -\frac{1}{\mu}\mathbf{K} \cdot \nabla p \qquad (1) $$

where $\mathbf{v}$ is the volume average Darcy's velocity, $\mathbf{K}$ the permeability tensor, $\mu$ the viscosity, and $\nabla p$ the pressure gradient. Both fluids are assumed to be incompressible, therefore, the mass conservation equation reduces to

$$ \nabla \cdot \mathbf{v} = 0 \qquad (2) $$

To differentiate between the two fluid phases, a fraction function $c$ is introduced which takes a value 1 for one fluid and 0 for the other one. The viscosity $\mu$ is redefined as

$$ \mu = c\mu_2 + (1-c)\mu_1 \qquad (3) $$

where $\mu_2$ and $\mu_1$ are the viscosities of the two fluids. $c$ evolves with time according to the following advection equation

$$ c_t + \mathbf{v} \cdot \nabla c = 0 \qquad (4) $$

where $c_t$ is the time derivative of the fraction function $c$.

The problem formulation is completed by assigning boundary conditions for $\mathbf{v}$, $p$ and $c$, as well as initial conditions for $c$:

$$ c(\mathbf{x}, t = 0) = c_0(\mathbf{x}). \qquad (5) $$

Pressure Dirichlet boundary conditions are prescribed on inlet and outlet boundaries:

$$ p(\mathbf{x}_{\text{inlet}}, t) = p_{\text{in}} \quad (\text{Inlet pressure}) \qquad (6) $$

and

$$ p(\mathbf{x}_{\text{outlet}}, t) = p_{\text{out}} \quad (\text{Outlet pressure}) \qquad (7) $$

Impermeable boundaries are characterized by zero normal velocity:

$$ \mathbf{v} \cdot \mathbf{n} = 0 \quad (\text{Impermeable wall}) \qquad (8) $$

Inlet flow also requires the assignment of boundary conditions for $c$:

$$ c(\mathbf{x}_{\text{inlet}}, t) = 1 \quad (\text{Inlet}) \qquad (9) $$

### 2.2. PINN

The first basic idea in PINN is the choice of the search space as a fully connected neural network approximation. The approximation has the form of

$$ \mathbf{z}^i = \sigma^i(\mathbf{W}^i \mathbf{z}^{i-1} + \mathbf{b}^i), \quad i = 1, \dots, L \qquad (10) $$

where $\mathbf{z}^0$ and $\mathbf{z}^L$ are the input (temporal or space dimensions) and output (an approximation to the solution) of the network, respectively. While, $\mathbf{W}^i$ and $\mathbf{b}^i$ are the parameters of each layer, known as the weights and biases, respectively. $\sigma^i$ is a nonlinear function used to add nonlinearity in the representation and is called the activation function. $L$ is the number of layers in the network.

Secondly, the residual of the PDE, to be solved, is obtained by differentiating the neural network output with respect to the inputs using automatic differentiation. A solution is reached by finding the weights and biases that minimize a loss function, composed of the PDEs and boundary conditions residual over a set of finite points namely, collocation and training points, respectively.
## 2. Model problem and PINN formulation

In this section, we first introduce the model problem for two-phase flow in porous media. The basics for PINN are then recalled. Finally, PINN is applied to solve the model problem.

### 2.1. Model problem for two-phase flow in porous media

Flow in porous media can be described by Darcy's law, which reads as

$$ \mathbf{v} = -\frac{1}{\mu}\mathbf{K} \cdot \nabla p \qquad (1) $$

where $\mathbf{v}$ is the volume average Darcy's velocity, $\mathbf{K}$ the permeability tensor, $\mu$ the viscosity, and $\nabla p$ the pressure gradient. Both fluids are assumed to be incompressible, therefore, the mass conservation equation reduces to

$$ \nabla \cdot \mathbf{v} = 0 \qquad (2) $$

To differentiate between the two fluid phases, a fraction function $c$ is introduced which takes a value 1 for one fluid and 0 for the other one. The viscosity $\mu$ is redefined as

$$ \mu = c\mu_2 + (1-c)\mu_1 \qquad (3) $$

where $\mu_2$ and $\mu_1$ are the viscosities of the two fluids. $c$ evolves with time according to the following advection equation

$$ c_t + \mathbf{v} \cdot \nabla c = 0 \qquad (4) $$

where $c_t$ is the time derivative of the fraction function $c$.

The problem formulation is completed by assigning boundary conditions for $\mathbf{v}$, $p$ and $c$, as well as initial conditions for $c$:

$$ c(\mathbf{x}, t=0) = c_0(\mathbf{x}). \qquad (5) $$

Pressure Dirichlet boundary conditions are prescribed on inlet and outlet boundaries:

$$ p(\mathbf{x}_{\text{inlet}}, t) = p_{\text{in}} \quad (\text{Inlet pressure}) \qquad (6) $$

and

$$ p(\mathbf{x}_{\text{outlet}}, t) = p_{\text{out}} \quad (\text{Outlet pressure}) \qquad (7) $$

Impermeable boundaries are characterized by zero normal velocity:

$$ \mathbf{v} \cdot \mathbf{n} = 0 \quad (\text{Impermeable wall}) \qquad (8) $$

Inlet flow also requires the assignment of boundary conditions for $c$:

$$ c(\mathbf{x}_{\text{inlet}}, t) = 1 \quad (\text{Inlet}) \qquad (9) $$

### 2.2. PINN

The first basic idea in PINN is the choice of the search space as a fully connected neural network approximation. The approximation has the form of

$$ \mathbf{z}^i = \sigma^i(\mathbf{W}^i \mathbf{z}^{i-1} + \mathbf{b}^i), \quad i = 1, \dots, L \qquad (10) $$

where $\mathbf{z}^0$ and $\mathbf{z}^L$ are the input (temporal or space dimensions) and output (an approximation to the solution) of the network, respectively. While, $\mathbf{W}^i$ and $\mathbf{b}^i$ are the parameters of each layer, known as the weights and biases, respectively. $\sigma^i$ is a nonlinear function used to add nonlinearity in the representation and is called the activation function. $L$ is the number of layers in the network.

Secondly, the residual of the PDE, to be solved, is obtained by differentiating the neural network output with respect to the inputs using automatic differentiation. A solution is reached by finding the weights and biases that minimize a loss function, composed of the PDEs and boundary conditions residual over a set of finite points namely, collocation and training points, respectively.
<figure><img src="image_3.png" /><figcaption>Fig. 1. PINN structure for a general two-phase flow in porous media problem.</figcaption></figure>

## 2.3. PINN structure and loss function definition for the model problem

To solve the problem using PINN, three distinct neural networks are used: one for the velocity, another for pressure, and another for the fraction function. Each of these networks has space x and time t inputs. The outputs of these neural networks (v, p and c) are differentiated with respect to the inputs, using automatic differentiation, forming the residuals of the three differential equations (4), (1) and (2). It should be noted that the network of v has 1, 2, or 3 outputs, according to the problem dimension, corresponding to the velocity components. Another possible structure is to have different networks for each velocity component. The PINN structure for a general two-phase flow in porous media problem is summarized in Fig. 1.

The loss function can, then, be defined as follows

$$ \text{Loss} = \lambda_v \text{loss}_v + \lambda_c \text{loss}_c + \lambda_p \text{loss}_p + \lambda_1 \text{loss}_{f1} + \lambda_2 \text{loss}_{f2} + \lambda_3 \text{loss}_{f3} \quad (11) $$

where

$$ \text{loss}_v = \frac{1}{N_v} \sum_{i=1}^{N_v} r_v^2(\mathbf{x}_v^i, t_v^i) = \frac{1}{N_v} \sum_{i=1}^{N_v} \|\mathbf{v}(\mathbf{x}_v^i, t_v^i) \cdot \mathbf{n}(\mathbf{x}_v^i, t_v^i)\|^2 \quad (12) $$

$$ \text{loss}_c = \frac{1}{N_c} \sum_{i=1}^{N_c} r_c^2(\mathbf{x}_c^i, t_c^i) = \frac{1}{N_c} \sum_{i=1}^{N_c} \|c(\mathbf{x}_c^i, t_c^i) - c_b^i\|^2 \quad (13) $$

$$ \text{loss}_p = \frac{1}{N_p} \sum_{i=1}^{N_p} r_p^2(\mathbf{x}_p^i, t_p^i) = \frac{1}{N_p} \sum_{i=1}^{N_p} \|p(\mathbf{x}_p^i, t_p^i) - p_b^i\|^2 \quad (14) $$

$$ \text{loss}_{f1} = \frac{1}{N_{f1}} \sum_{i=1}^{N_{f1}} \|f_1(\mathbf{x}_{f1}^i, t_{f1}^i)\|^2 = \frac{1}{N_{f1}} \sum_{i=1}^{N_{f1}} \|c_t + \mathbf{v} \cdot \nabla c\|_{(\mathbf{x}_{f1}^i, t_{f1}^i)}^2 \quad (15) $$

$$ \text{loss}_{f2} = \frac{1}{N_{f2}} \sum_{i=1}^{N_{f2}} \|f_2(\mathbf{x}_{f2}^i, t_{f2}^i)\|^2 = \frac{1}{N_{f2}} \sum_{i=1}^{N_{f2}} \|\mathbf{v} + \frac{1}{\mu} \mathbf{K} \cdot \nabla p\|_{(\mathbf{x}_{f2}^i, t_{f2}^i)}^2 \quad (16) $$

$$ \text{loss}_{f3} = \frac{1}{N_{f3}} \sum_{i=1}^{N_{f3}} \|f_3(\mathbf{x}_{f3}^i, t_{f3}^i)\|^2 = \frac{1}{N_{f3}} \sum_{i=1}^{N_{f3}} \|\nabla \cdot \mathbf{v}\|_{(\mathbf{x}_{f3}^i, t_{f3}^i)}^2 \quad (17) $$

where {$\mathbf{x}_v^i, t_v^i, v_b^i$}$_{i=1}^{N_v}$, {$\mathbf{x}_c^i, t_c^i, c_b^i$}$_{i=1}^{N_c}$ and {$\mathbf{x}_p^i, t_p^i, p_b^i$}$_{i=1}^{N_p}$ are the points where the initial/boundary conditions are defined for v, c and p, respectively. While {$\mathbf{x}_{f1}^i, t_{f1}^i$}$_{i=1}^{N_{f1}}$, {$\mathbf{x}_{f2}^i, t_{f2}^i$}$_{i=1}^{N_{f2}}$ and {$\mathbf{x}_{f3}^i, t_{f3}^i$}$_{i=1}^{N_{f3}}$ are the collocation points in space and time for the three residuals, respectively, where the physics is enforced and $\lambda_i$ are the weights of
<figure><img src="image_3.png" /><figcaption>Fig. 1. PINN structure for a general two-phase flow in porous media problem.</figcaption></figure>

## 2.3. PINN structure and loss function definition for the model problem

To solve the problem using PINN, three distinct neural networks are used: one for the velocity, another for pressure, and another for the fraction function. Each of these networks has space x and time t inputs. The outputs of these neural networks (v, p and c) are differentiated with respect to the inputs, using automatic differentiation, forming the residuals of the three differential equations (4), (1) and (2). It should be noted that the network of v has 1, 2, or 3 outputs, according to the problem dimension, corresponding to the velocity components. Another possible structure is to have different networks for each velocity component. The PINN structure for a general two-phase flow in porous media problem is summarized in Fig. 1.

The loss function can, then, be defined as follows

$$ \text{Loss} = \lambda_v \text{loss}_v + \lambda_c \text{loss}_c + \lambda_p \text{loss}_p + \lambda_1 \text{loss}_{f_1} + \lambda_2 \text{loss}_{f_2} + \lambda_3 \text{loss}_{f_3} \quad (11) $$

where

$$ \text{loss}_v = \frac{1}{N_v} \sum_{i=1}^{N_v} r_v^2(\mathbf{x}_v^i, t_v^i) = \frac{1}{N_v} \sum_{i=1}^{N_v} \|\mathbf{v}(\mathbf{x}_v^i, t_v^i) \cdot \mathbf{n}(\mathbf{x}_v^i, t_v^i)\|^2 \quad (12) $$

$$ \text{loss}_c = \frac{1}{N_c} \sum_{i=1}^{N_c} r_c^2(\mathbf{x}_c^i, t_c^i) = \frac{1}{N_c} \sum_{i=1}^{N_c} \|c(\mathbf{x}_c^i, t_c^i) - c_b^i\|^2 \quad (13) $$

$$ \text{loss}_p = \frac{1}{N_p} \sum_{i=1}^{N_p} r_p^2(\mathbf{x}_p^i, t_p^i) = \frac{1}{N_p} \sum_{i=1}^{N_p} \|p(\mathbf{x}_p^i, t_p^i) - p_b^i\|^2 \quad (14) $$

$$ \text{loss}_{f1} = \frac{1}{N_{f1}} \sum_{i=1}^{N_{f1}} \|f_1(\mathbf{x}_{f1}^i, t_{f1}^i)\|^2 = \frac{1}{N_{f1}} \sum_{i=1}^{N_{f1}} \|c_i + \mathbf{v} \cdot \nabla c\|_{(\mathbf{x}_{f1}^i, t_{f1}^i)}^2 \quad (15) $$

$$ \text{loss}_{f2} = \frac{1}{N_{f2}} \sum_{i=1}^{N_{f2}} \|f_2(\mathbf{x}_{f2}^i, t_{f2}^i)\|^2 = \frac{1}{N_{f2}} \sum_{i=1}^{N_{f2}} \|\mathbf{v} + \frac{1}{\mu} \mathbf{K} \cdot \nabla p\|_{(\mathbf{x}_{f2}^i, t_{f2}^i)}^2 \quad (16) $$

$$ \text{loss}_{f3} = \frac{1}{N_{f3}} \sum_{i=1}^{N_{f3}} \|f_3(\mathbf{x}_{f3}^i, t_{f3}^i)\|^2 = \frac{1}{N_{f3}} \sum_{i=1}^{N_{f3}} \|\nabla \cdot \mathbf{v}\|_{(\mathbf{x}_{f3}^i, t_{f3}^i)}^2 \quad (17) $$

where {$\mathbf{x}_v^i, t_v^i, v_b^i$}$_{i=1}^{N_v}$, {$\mathbf{x}_c^i, t_c^i, c_b^i$}$_{i=1}^{N_c}$ and {$\mathbf{x}_p^i, t_p^i, p_b^i$}$_{i=1}^{N_p}$ are the points where the initial/boundary conditions are defined for v, c and p, respectively. While {$\mathbf{x}_{f1}^i, t_{f1}^i$}$_{i=1}^{N_{f1}}$, {$\mathbf{x}_{f2}^i, t_{f2}^i$}$_{i=1}^{N_{f2}}$ and {$\mathbf{x}_{f3}^i, t_{f3}^i$}$_{i=1}^{N_{f3}}$ are the collocation points in space and time for the three residuals, respectively, where the physics is enforced and $\lambda_i$ are the weights of
each term in the loss function. The $\lambda_i$ weights are important to make each contribution to the loss function has comparable magnitude, thus, aiding the optimization process. Finally, a solution for the fields, p, v and c, is obtained by minimizing the loss function with respect to the neural networks' parameters.

### 3. Collocation points adaptivity

Mesh Refinement is a basic idea in classical numerical techniques such as finite element (FEM) and finite volume (FVM) methods [40]. There are three basic techniques to mesh adaptation: h-adaptivity [41], r-adaptivity [42] and p-adaptivity [43]. H-adaptivity adds more nodes, thus increasing the degrees of freedom and the mesh connectivity. While r-adaptivity keeps the same number of nodes and degrees of freedom however the nodes are relocated while keeping the same connectivity. Finally, p-adaptivity increases the polynomial degrees of elements while keeping the mesh fixed. Other adaptivity methods exist that combines some of the basic techniques together such as: hp-adaptivity [44] and rh-adaptivity [43].

There are mainly three drivers for mesh adaptation: error [45], PDE residual [46] and solution features [47].

Error-based adaptation adds more degrees of freedom where the solution error is high. This technique ignores the fact that the error is transported in the domain. Therefore, adapting where the error is high might ignore the region of the error source itself, where adaption is more useful [48].

Residual-based adaptation refines the mesh where the discretized PDE residual is high. The residual can be seen as the source of error in the solution [49]. Thus, refining where the residual is high is seen as a way of refining where the error source is. Therefore, this technique usually performs better than error-based adaptation.

Solution-based adaptation utilizes the solution features such as gradients or discontinuities for adaptation. The philosophy behind this technique is that by using more points in these locations, these features can be resolved, thus, leading to improving the overall accuracy of the solution. However, if multiple features are present in the same problem, the adaptation results in over-refining some features while others are ignored. An example of this adaptivity failure can be found in [50].

The adaptation process is usually computationally expensive since certain requirements have to be satisfied and the mesh connectivity needs to be updated. Moreover, parallelization becomes complex for unstructured grids. However, in the case of PINN, changing the collocation points locations or adding more points are cheap processes. The main reason is that PINN is a meshless method, thus, there are no element volumes to take care of or mesh connectivity to update. Moreover, the approximation of the derivatives is independent of the collocation point position. Therefore, there is no discretization error resulting from the distribution of the collocation points. The only thing to do is identify the location where more points are needed.

In this study, we develop a residual-based algorithm by enriching the locations where the residual is high with more collocation points. We build the algorithm based on the work of [51] in which the authors developed the residual-based adaptive refinement method (RAR). In their work, the authors used a dense set of randomly drawn points in the space-time domain, where residuals are evaluated. Points corresponding to the largest residual values are then added to the training set of collocation points. The progressive refinement of the training set allows for residual control. However when residual is showing high values in very narrow regions, this sampling strategy tends to produce excessively clustered points ignoring other solution features if existing and leading to unnecessary over-refinement. This behavior appears to be related to solutions exhibiting moving sharp fronts or discontinuities, as in the model considered in this work. An example of that is seen in Fig. 2 where the added points are focused in a very small region.

To avoid this point clusters that might lead to over-fitting of the model, we designed a probability density function based on the residual field to control the spread of the added points. The density function is, then, used to draw points from the dense set and these points are added to the training set. By doing this, collocation points will be more evenly spread in the domain; more points will be added where the probability is high (high residual) and fewer points where the probability is lower (low residual). Fig. 3 explains the procedure used for the adaptivity technique developed in this paper.

The algorithm is extended for coupled differential equations so that each PDE residual will have different collocation points. Moreover, data points are also enriched using a similar strategy to capture the initial/boundary conditions better. The algorithm is as follows:
each term in the loss function. The $\lambda_i$ weights are important to make each contribution to the loss function has comparable magnitude, thus, aiding the optimization process. Finally, a solution for the fields, p, v and c, is obtained by minimizing the loss function with respect to the neural networks' parameters.

### 3. Collocation points adaptivity

Mesh Refinement is a basic idea in classical numerical techniques such as finite element (FEM) and finite volume (FVM) methods [40]. There are three basic techniques to mesh adaptation: h-adaptivity [41], r-adaptivity [42] and p-adaptivity [43]. H-adaptivity adds more nodes, thus increasing the degrees of freedom and the mesh connectivity. While r-adaptivity keeps the same number of nodes and degrees of freedom however the nodes are relocated while keeping the same connectivity. Finally, p-adaptivity increases the polynomial degrees of elements while keeping the mesh fixed. Other adaptivity methods exist that combines some of the basic techniques together such as: hp-adaptivity [44] and rh-adaptivity [43].

There are mainly three drivers for mesh adaptation: error [45], PDE residual [46] and solution features [47].

Error-based adaptation adds more degrees of freedom where the solution error is high. This technique ignores the fact that the error is transported in the domain. Therefore, adapting where the error is high might ignore the region of the error source itself, where adaption is more useful [48].

Residual-based adaptation refines the mesh where the discretized PDE residual is high. The residual can be seen as the source of error in the solution [49]. Thus, refining where the residual is high is seen as a way of refining where the error source is. Therefore, this technique usually performs better than error-based adaptation.

Solution-based adaptation utilizes the solution features such as gradients or discontinuities for adaptation. The philosophy behind this technique is that by using more points in these locations, these features can be resolved, thus, leading to improving the overall accuracy of the solution. However, if multiple features are present in the same problem, the adaptation results in over-refining some features while others are ignored. An example of this adaptivity failure can be found in [50].

The adaptation process is usually computationally expensive since certain requirements have to be satisfied and the mesh connectivity needs to be updated. Moreover, parallelization becomes complex for unstructured grids. However, in the case of PINN, changing the collocation points locations or adding more points are cheap processes. The main reason is that PINN is a meshless method, thus, there are no element volumes to take care of or mesh connectivity to update. Moreover, the approximation of the derivatives is independent of the collocation point position. Therefore, there is no discretization error resulting from the distribution of the collocation points. The only thing to do is identify the location where more points are needed.

In this study, we develop a residual-based algorithm by enriching the locations where the residual is high with more collocation points. We build the algorithm based on the work of [51] in which the authors developed the residual-based adaptive refinement method (RAR). In their work, the authors used a dense set of randomly drawn points in the space-time domain, where residuals are evaluated. Points corresponding to the largest residual values are then added to the training set of collocation points. The progressive refinement of the training set allows for residual control. However when residual is showing high values in very narrow regions, this sampling strategy tends to produce excessively clustered points ignoring other solution features if existing and leading to unnecessary over-refinement. This behavior appears to be related to solutions exhibiting moving sharp fronts or discontinuities, as in the model considered in this work. An example of that is seen in Fig. 2 where the added points are focused in a very small region.

To avoid this point clusters that might lead to over-fitting of the model, we designed a probability density function based on the residual field to control the spread of the added points. The density function is, then, used to draw points from the dense set and these points are added to the training set. By doing this, collocation points will be more evenly spread in the domain; more points will be added where the probability is high (high residual) and fewer points where the probability is lower (low residual). Fig. 3 explains the procedure used for the adaptivity technique developed in this paper.

The algorithm is extended for coupled differential equations so that each PDE residual will have different collocation points. Moreover, data points are also enriched using a similar strategy to capture the initial/boundary conditions better. The algorithm is as follows:
**Algorithm 1: Residual-based adaptivity**

Inputs: Number of adaptivity steps *M*, number of iterations *n*, tolerances $\epsilon_i$, $\epsilon_v$, $\epsilon_c$ and $\epsilon_p$;

while $m < M$ & ($\mu_i > \epsilon_i$ || $\mu_v > \epsilon_v$ || $\mu_c > \epsilon_c$ || $\mu_p > \epsilon_p$) do

- Calculate $f_i$ (3 PDE residuals), $r_v$, $r_c$ and $r_p$ using dense sets;

- Build the probability functions $p_i$, $p_v$, $p_c$ and $p_p$ using dense sets Eq: (18);

- Draw points from the dense sets using the probabilities and add it to the training sets;

- Use minimization algorithm for *n* fixed iterations;

- Calculate the residuals' mean using the dense sets;

$$ \mu_i = \frac{1}{N_i} \sum |f_i| $$

$$ \mu_v = \frac{1}{N_v} \sum |\mathbf{v}(\mathbf{x}_v^i, t_v^i) \cdot \mathbf{n}(\mathbf{x}_v^i, t_v^i)| $$

$$ \mu_c = \frac{1}{N_c} \sum |c(x_c^i, t_c^i) - c_b^i| $$

$$ \mu_p = \frac{1}{N_p} \sum |p(x_p^i, t_p^i) - p_b^i| $$

end

The probability functions used have the form of

$$ p(\mathbf{X}) = \frac{\max(\log|r(\mathbf{X})/\epsilon|, 0)}{\int_{\Omega} \max(\log|r(\mathbf{X})/\epsilon|, 0) d\mathbf{X}} \quad (18) $$

where $\mathbf{X}$ is the random vector $[x, t]^T$, $r$ the considered residual, $\Omega$ the spatio-temporal domain and $\epsilon$ a small tolerance to filter small residual values. In practice, the choice of the value of $\epsilon$ is chosen to control the spread of the point distributions. The function is designed in a way to ensure that its integral over the space-time domain is 1, hence, the presence of the term in the denominator which is calculated using Monte Carlo integration over the dense set of points.

# 4. Numerical experiments

## 4.1. One dimensional injection

We consider a one-dimensional problem shown in Fig. 4. At $t = 0$, the domain is filled with one fluid (fluid 1). Another fluid (fluid 2) is being injected from the left end at constant pressure $p_{in}$, while the pressure at the other end is fixed to $p_{out}$.

Eq. (1) will reduce to

$$ v = -\frac{k}{\mu} p_x \quad (19) $$

where $v$ is the volume average Darcy's velocity, $k$ the permeability, $\mu$ the viscosity and $p_x$ the pressure gradient with respect to $x$. Eq. (2) reduces to

$$ v_x = 0 \quad (20) $$

Boundary conditions are given for the pressure inlet and outlet at all times.

$$ p(0, t) = p_{in}, \quad p(l, t) = p_{out} \quad (21) $$

where *l* is the domain length.

To differentiate between the two phases, a fraction function *c* is defined. *c* takes a value of 1 in the domain that is filled with fluid 2 and 0 when fluid 1 is present. The viscosity $\mu$ is redefined according to (3). Eq. (4) reduces
**Algorithm 1: Residual-based adaptivity**

Inputs: Number of adaptivity steps M, number of iterations n, tolerances $\epsilon_i$, $\epsilon_v$, $\epsilon_c$ and $\epsilon_p$;

while $m < M$ & ($\mu_i > \epsilon_i$ || $\mu_v > \epsilon_v$ || $\mu_c > \epsilon_c$ || $\mu_p > \epsilon_p$) do

- Calculate $f_i$ (3 PDE residuals), $r_v$, $r_c$ and $r_p$ using dense sets;

- Build the probability functions $p_i$, $p_v$, $p_c$ and $p_p$ using dense sets Eq: (18);

- Draw points from the dense sets using the probabilities and add it to the training sets;

- Use minimization algorithm for n fixed iterations;

- Calculate the residuals' mean using the dense sets;

$$ \mu_i = \frac{1}{N_i} \sum |f_i| $$

$$ \mu_v = \frac{1}{N_v} \sum |\mathbf{v}(\mathbf{x}_v^i, t_v^i) \cdot \mathbf{n}(\mathbf{x}_v^i, t_v^i)| $$

$$ \mu_c = \frac{1}{N_c} \sum |c(x_c^i, t_c^i) - c_b^i| $$

$$ \mu_p = \frac{1}{N_p} \sum |p(x_p^i, t_p^i) - p_b^i| $$

end

The probability functions used have the form of

$$ p(\mathbf{X}) = \frac{\max(\log|r(\mathbf{X})/\epsilon|, 0)}{\int_{\Omega} \max(\log|r(\mathbf{X})/\epsilon|, 0) d\mathbf{X}} \quad (18) $$

where $\mathbf{X}$ is the random vector $[x, t]^T$, $r$ the considered residual, $\Omega$ the spatio-temporal domain and $\epsilon$ a small tolerance to filter small residual values. In practice, the choice of the value of $\epsilon$ is chosen to control the spread of the point distributions. The function is designed in a way to ensure that its integral over the space-time domain is 1, hence, the presence of the term in the denominator which is calculated using Monte Carlo integration over the dense set of points.

# 4. Numerical experiments

## 4.1. One dimensional injection

We consider a one-dimensional problem shown in Fig. 4. At $t = 0$, the domain is filled with one fluid (fluid 1). Another fluid (fluid 2) is being injected from the left end at constant pressure $p_{in}$, while the pressure at the other end is fixed to $p_{out}$.

Eq. (1) will reduce to

$$ v = -\frac{k}{\mu} p_x \quad (19) $$

where $v$ is the volume average Darcy's velocity, $k$ the permeability, $\mu$ the viscosity and $p_x$ the pressure gradient with respect to $x$. Eq. (2) reduces to

$$ v_x = 0 \quad (20) $$

Boundary conditions are given for the pressure inlet and outlet at all times.

$$ p(0, t) = p_{in}, \quad p(l, t) = p_{out} \quad (21) $$

where $l$ is the domain length.

To differentiate between the two phases, a fraction function $c$ is defined. $c$ takes a value of 1 in the domain that is filled with fluid 2 and 0 when fluid 1 is present. The viscosity $\mu$ is redefined according to (3). Eq. (4) reduces
<figure><img src="image_4.png" /><figcaption>Fig. 2. Left: log of the absolute residual field $f_1$, right: the chosen points according to RAR.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Left: log of the absolute residual field $f_1$, middle: the probability density function built from the residual, right: the chosen points drawn from the density function.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. One-dimensional domain (filling problem).</figcaption></figure>

to

$$c_t + v c_x = 0 \qquad (22)$$

where $c_t$ and $c_x$ are time and position gradients of $c$, respectively.

To solve this equation, initial and boundary conditions need to be enforced on $c$.

$$c(0, t) = 1, \quad c(x, 0) = 0 \qquad (23)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Left: log of the absolute residual field $f_1$, right: the chosen points according to RAR.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Left: log of the absolute residual field $f_1$, middle: the probability density function built from the residual, right: the chosen points drawn from the density function.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. One-dimensional domain (filling problem).</figcaption></figure>

to

$$c_t + v c_x = 0 \qquad (22)$$

where $c_t$ and $c_x$ are time and position gradients of $c$, respectively.

To solve this equation, initial and boundary conditions need to be enforced on $c$.

$$c(0, t) = 1, \quad c(x, 0) = 0 \qquad (23)$$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Front (fluid 1/fluid 2 interface) position with time for the fixed, RAR and new adaptivity cases along with analytical front position for the 1D example and a zoomed image to differentiate.</figcaption></figure>

To sum up, the strong form of the problem can be written as:

$$
\begin{align}
c_t + v c_x &= 0, \quad x \in [0, l], \quad t \in [0, T], \nonumber \\
v &= -k/\mu p_x, \quad x \in [0, l], \quad t \in [0, T], \nonumber \\
v_x &= 0, \quad x \in [0, l], \quad t \in [0, T], \nonumber \\
\mu &= c\mu_2 + (1-c)\mu_1 \nonumber \\
p(0, t) &= p_{in}, \quad p(l, t) = p_{out} \nonumber \\
c(0, t) &= 1, \quad c(x, 0) = 0 \tag{24}
\end{align}
$$

where T is the end of the time domain where the problem is solved.

The training is performed using a neural network of 5 hidden layers and 20 neurons in each layer. Hyperbolic tangent activation function is used in all the hidden layers; this choice proved to provide good results [35,36,38]. However, for the output layer sigmoid function is used for the pressure and fraction function networks since their values go from 0 to 1, while, linear activation function is used for the output layer of the velocity network. Adam optimizer is used for the initial training phase with a fixed learning rate of 0.001 and 1000 iterations. Afterwards, BFGS optimization algorithm is used. This optimization strategy is followed by [31,51,52]. The reason is that second-order methods like BFGS are prone to fall into local minima, therefore, Adam is firstly used to reach the zone of the global minimum, afterwards, BFGS is used to reach the minimum easier since it iteratively builds an approximation of the Hessian matrix. The different terms of the loss function are weighted similarly, thus, $\lambda_c = \lambda_p = \lambda_1 = \lambda_2 = \lambda_3 = 1$. The parameters used to solve the problem are given in Table 1.

Three numerical experiments are performed and compared. The first one using fixed number and position of collocation points (2500 points organized in 50 × 50 grid) for the whole training phase. The second experiment using the RAR technique starting with 1600 points organized as 40 × 40 grid points during the Adam training phase. Afterwards, point enrichments are performed every 50 BFGS iterations till the stopping criteria are satisfied (at 2500 points as well). The final experiment using the provided adaptivity algorithm starting with 1600 points organized as 40 × 40 grid points during the Adam training phase. Similar enrichments are done as in the RAR experiment. All cases took near 200 s to converge using laptop Intel core i7-6700HQ CPU @ 2.60 GHz 2.59 GHz with 8 Go RAM.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Front (fluid 1/fluid 2 interface) position with time for the fixed, RAR and new adaptivity cases along with analytical front position for the 1D example and a zoomed image to differentiate.</figcaption></figure>

To sum up, the strong form of the problem can be written as:

$$
\begin{align}
c_t + v c_x &= 0, \quad x \in [0, l], \quad t \in [0, T], \nonumber \\
v &= -k/\mu p_x, \quad x \in [0, l], \quad t \in [0, T], \nonumber \\
v_x &= 0, \quad x \in [0, l], \quad t \in [0, T], \nonumber \\
\mu &= c\mu_2 + (1-c)\mu_1 \nonumber \\
p(0, t) &= p_{in}, \quad p(l, t) = p_{out} \nonumber \\
c(0, t) &= 1, \quad c(x, 0) = 0 \tag{24}
\end{align}
$$

where *T* is the end of the time domain where the problem is solved.

The training is performed using a neural network of 5 hidden layers and 20 neurons in each layer. Hyperbolic tangent activation function is used in all the hidden layers; this choice proved to provide good results [35,36,38]. However, for the output layer sigmoid function is used for the pressure and fraction function networks since their values go from 0 to 1, while, linear activation function is used for the output layer of the velocity network. Adam optimizer is used for the initial training phase with a fixed learning rate of 0.001 and 1000 iterations. Afterwards, BFGS optimization algorithm is used. This optimization strategy is followed by [31,51,52]. The reason is that second-order methods like BFGS are prone to fall into local minima, therefore, Adam is firstly used to reach the zone of the global minimum, afterwards, BFGS is used to reach the minimum easier since it iteratively builds an approximation of the Hessian matrix. The different terms of the loss function are weighted similarly, thus, $\lambda_c = \lambda_p = \lambda_1 = \lambda_2 = \lambda_3 = 1$. The parameters used to solve the problem are given in Table 1.

Three numerical experiments are performed and compared. The first one using fixed number and position of collocation points (2500 points organized in 50 × 50 grid) for the whole training phase. The second experiment using the RAR technique starting with 1600 points organized as 40 × 40 grid points during the Adam training phase. Afterwards, point enrichments are performed every 50 BFGS iterations till the stopping criteria are satisfied (at 2500 points as well). The final experiment using the provided adaptivity algorithm starting with 1600 points organized as 40 × 40 grid points during the Adam training phase. Similar enrichments are done as in the RAR experiment. All cases took near 200 s to converge using laptop Intel core i7-6700HQ CPU @ 2.60 GHz 2.59 GHz with 8 Go RAM.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Pressure profiles at different times ($t = 0, t = 0.15, t = 0.25$ and $t = 0.5$) for the adaptive and fixed collocations cases for the 1D example.</figcaption></figure>

<table><caption>**Table 1**
Parameters used for the one-dimensional two-phase flow.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>$l$</td><td>1</td></tr><tr><td>$T$</td><td>0.5</td></tr><tr><td>$k$</td><td>1</td></tr><tr><td>$\mu_2$</td><td>1</td></tr><tr><td>$\mu_1$</td><td>10<sup>-5</sup></td></tr><tr><td>$p_{in}$</td><td>1</td></tr><tr><td>$p_{out}$</td><td>0</td></tr></tbody></table>

A fixed 1000 points randomly distributed over the domain are used as a test set to have a sense of the generalization error committed during the training phase. It should be noted that these points are only used for testing but are not used in the training phase.

Firstly, the flow front (fluid 1/fluid 2 interface) is extracted as the 0.5 level set of $c$ and plotted for both cases along with the analytical front solution in Fig. 5. The analytical solution for the flow front $x_f$ is obtained as

$$x_f = \frac{-\mu_1 l + \sqrt{\mu_1^2 l^2 + 2(\mu_2 - \mu_1)k(p_{in} - p_{out})t}}{\mu_2 - \mu_1} \quad (25)$$

The new adaptivity technique provided the best results among the three cases. PINN with RAR did not provide a significant improvement to the approximation; that is probably due to focusing of the enrichments in a small spatio-temporal region leading to harder optimization and higher generalization error.

The pressure profiles at different times for the three cases are shown in Fig. 6 and compared to the analytical solution. The analytical solution for pressure can be written as

$$p(x, t) = \begin{cases} \frac{-\mu_2(p_{in} - p_{out})}{(\mu_2 - \mu_1)x_f(t) + \mu_1 l} & x + p_{in} & x < x_f(t) \\ \frac{-\mu_1(p_{in} - p_{out})}{(\mu_2 - \mu_1)x_f(t) + \mu_1 l} & x + \frac{\mu_1(p_{in} - p_{out})}{(\mu_2 - \mu_1)x_f(t) + \mu_1 l} & l + p_{out} & x \ge x_f(t) \end{cases} \quad (26)$$
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Pressure profiles at different times ($t = 0, t = 0.15, t = 0.25$ and $t = 0.5$) for the adaptive and fixed collocations cases for the 1D example.</figcaption></figure>

<table><caption>**Table 1**
Parameters used for the one-dimensional two-phase flow.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>$l$</td><td>1</td></tr><tr><td>$T$</td><td>0.5</td></tr><tr><td>$k$</td><td>1</td></tr><tr><td>$\mu_2$</td><td>1</td></tr><tr><td>$\mu_1$</td><td>$10^{-5}$</td></tr><tr><td>$p_{in}$</td><td>1</td></tr><tr><td>$p_{out}$</td><td>0</td></tr></tbody></table>

A fixed 1000 points randomly distributed over the domain are used as a test set to have a sense of the generalization error committed during the training phase. It should be noted that these points are only used for testing but are not used in the training phase.

Firstly, the flow front (fluid 1/fluid 2 interface) is extracted as the 0.5 level set of $c$ and plotted for both cases along with the analytical front solution in Fig. 5. The analytical solution for the flow front $x_f$ is obtained as

$$x_f = \frac{-\mu_1 l + \sqrt{\mu_1^2 l^2 + 2(\mu_2 - \mu_1)k(p_{in} - p_{out})t}}{\mu_2 - \mu_1} \quad (25)$$

The new adaptivity technique provided the best results among the three cases. PINN with RAR did not provide a significant improvement to the approximation; that is probably due to focusing of the enrichments in a small spatio-temporal region leading to harder optimization and higher generalization error.

The pressure profiles at different times for the three cases are shown in Fig. 6 and compared to the analytical solution. The analytical solution for pressure can be written as

$$p(x, t) = \begin{cases} \frac{-\mu_2(p_{in} - p_{out})}{(\mu_2 - \mu_1)x_f(t) + \mu_1 l} & x + p_{in} & x < x_f(t) \\ \frac{-\mu_1(p_{in} - p_{out})}{(\mu_2 - \mu_1)x_f(t) + \mu_1 l} & x + \frac{\mu_1(p_{in} - p_{out})}{(\mu_2 - \mu_1)x_f(t) + \mu_1 l} & l + p_{out} & x \ge x_f(t) \end{cases} \quad (26)$$
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Evolution of collocation points with new adaptivity (from left to right: collocations for $f_1$, collocations for $f_2$ and collocations for $f_3$) for the 1D example. New collocation points after each adaptation step are shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The new adaptive case provided a pressure solution closer to the analytical solution than both of the fixed PINN and PINN with RAR cases. PINN with RAR shows a prediction far from the analytical solution near $t = 0$.

The distribution of the collocation points for the different PDEs is shown in the case of new adaptive PINN in Fig. 7. The figure shows the evolution of these distributions at different stages of using the adaptivity algorithm.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Evolution of collocation points with new adaptivity (from left to right: collocations for $f_1$, collocations for $f_2$ and collocations for $f_3$) for the 1D example. New collocation points after each adaptation step are shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The new adaptive case provided a pressure solution closer to the analytical solution than both of the fixed PINN and PINN with RAR cases. PINN with RAR shows a prediction far from the analytical solution near $t = 0$.

The distribution of the collocation points for the different PDEs is shown in the case of new adaptive PINN in Fig. 7. The figure shows the evolution of these distributions at different stages of using the adaptivity algorithm.
<figure><img src="image_10.png" /><figcaption>Fig. 8. Loss vs. Iteration graphs for fixed PINN (left plot), new adaptive PINN (middle plot), and PINN with RAR (right plot) for the 1D example.</figcaption></figure>

From Fig. 7, it can be seen for the collocation points of $f_1$ that the points seem to be dense near the location of the front, where the residual is higher. This helped in capturing the interface location accurately. While for the other collocation points ($f_2$ and $f_3$), they were distributed almost randomly in the domain. That is because the residual field is spread all over the domain since there are no sharp solution features to capture.

The loss function is compared in the three cases by plotting the loss vs. iteration graph for both cases. The loss using the training set is compared to that using the testing set for all cases as shown in Fig. 8. It should be noted that there is a deviation between the training and testing loss in the fixed collocations case, meaning that the generalization error is high. This deviation is marginally decreased in the PINN with RAR case. For the new adaptive PINN case, the deviation significantly decreased meaning that less generalization error is committed using the newly-developed adaptive technique. From a deep learning perspective, using adaptive collocation points can be seen as a form of regularization of the neural network solution; adaptivity prevents overfitting, thus making the solution more accurate for unseen points (points not used in the training process).

The different terms in the loss function are plotted for the new adaptive case against the number of iterations in Fig. 9. This is done to assure the convergence of all the terms.

All the loss terms have more or less a similar trend in convergence. We can note that the loss corresponding to the initial and boundary conditions of the fraction function $c$ has the highest values (harder convergence) that is due to the discontinuity in the initial/boundary condition values at $t = 0$.

## 4.2. Two dimensional central injection

The next example is a two-dimensional problem that is encountered in structural composites manufacturing processes such as resin transfer molding. The domain is a square of a unity area with an elliptic injection port placed at the center with constant pressure ($p = 1$). The four outer sides are outlets where the pressure is set to zero. The problem domain is plotted in Fig. 10. The analytical solution of this problem exists in [53].

The material properties are the same as in the first example except for the permeability of the domain which is uniform orthotropic and can written as:

$$ \mathbf{K} = \begin{bmatrix} 1.5 & 0 \\ 0 & 1 \end{bmatrix} \qquad (27) $$

All training parameters, that are used in the first example, are the same in this example except for the network architecture. In this case, 4 neural networks are used to approximate the fraction function, pressure and two velocity components; each network is composed of 5 hidden layers and 20 neurons.

The fraction function and pressure fields are plotted at different times for the adaptive case to visualize the evolution of the flow front and pressure with time in Fig. 11.

The flow front position in x and y directions is plotted in Figs. 12 and 13 for the 3 cases along with the analytical solution to compare.

It can be seen that the new adaptive algorithm provided better prediction of the flow front position.
<figure><img src="image_10.png" /><figcaption>Fig. 8. Loss vs. Iteration graphs for fixed PINN (left plot), new adaptive PINN (middle plot), and PINN with RAR (right plot) for the 1D example.</figcaption></figure>

From Fig. 7, it can be seen for the collocation points of $f_1$ that the points seem to be dense near the location of the front, where the residual is higher. This helped in capturing the interface location accurately. While for the other collocation points ($f_2$ and $f_3$), they were distributed almost randomly in the domain. That is because the residual field is spread all over the domain since there are no sharp solution features to capture.

The loss function is compared in the three cases by plotting the loss vs. iteration graph for both cases. The loss using the training set is compared to that using the testing set for all cases as shown in Fig. 8. It should be noted that there is a deviation between the training and testing loss in the fixed collocations case, meaning that the generalization error is high. This deviation is marginally decreased in the PINN with RAR case. For the new adaptive PINN case, the deviation significantly decreased meaning that less generalization error is committed using the newly-developed adaptive technique. From a deep learning perspective, using adaptive collocation points can be seen as a form of regularization of the neural network solution; adaptivity prevents overfitting, thus making the solution more accurate for unseen points (points not used in the training process).

The different terms in the loss function are plotted for the new adaptive case against the number of iterations in Fig. 9. This is done to assure the convergence of all the terms.

All the loss terms have more or less a similar trend in convergence. We can note that the loss corresponding to the initial and boundary conditions of the fraction function $c$ has the highest values (harder convergence) that is due to the discontinuity in the initial/boundary condition values at $t = 0$.

## 4.2. Two dimensional central injection

The next example is a two-dimensional problem that is encountered in structural composites manufacturing processes such as resin transfer molding. The domain is a square of a unity area with an elliptic injection port placed at the center with constant pressure ($p = 1$). The four outer sides are outlets where the pressure is set to zero. The problem domain is plotted in Fig. 10. The analytical solution of this problem exists in [53].

The material properties are the same as in the first example except for the permeability of the domain which is uniform orthotropic and can written as:

$$ K = \begin{bmatrix} 1.5 & 0 \\ 0 & 1 \end{bmatrix} \qquad (27) $$

All training parameters, that are used in the first example, are the same in this example except for the network architecture. In this case, 4 neural networks are used to approximate the fraction function, pressure and two velocity components; each network is composed of 5 hidden layers and 20 neurons.

The fraction function and pressure fields are plotted at different times for the adaptive case to visualize the evolution of the flow front and pressure with time in Fig. 11.

The flow front position in x and y directions is plotted in Figs. 12 and 13 for the 3 cases along with the analytical solution to compare.

It can be seen that the new adaptive algorithm provided better prediction of the flow front position.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Different loss terms versus iterations for the new adaptivity case for the 1D example.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Domain of the 2D central injection problem. $\Gamma_{out}$ is the outlet boundary where the pressure is set to 0, while $\Gamma_{in}$ is the boundary of the inlet where pressure is set to 1.</figcaption></figure>

The loss function is plotted for the 3 cases in Fig. 14. 10,000 points in the space–time domain are chosen randomly to assess the testing loss, while they are not used in the training.

Using the new adaptivity algorithm, the discrepancy between the training and testing loss is greatly reduced. That means that it offers a mean to reduce the generalization error thus provides better accuracy to the approximated fields.

**5. Conclusion**

In this study, Physics-Informed Neural Networks (PINN) are applied to solve two-phase flow problems in porous media. A novel residual-based adaptive algorithm is developed. The key is to utilize the PDE residual to build
<figure><img src="image_11.png" /><figcaption>Fig. 9. Different loss terms versus iterations for the new adaptivity case for the 1D example.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Domain of the 2D central injection problem. $\Gamma_{out}$ is the outlet boundary where the pressure is set to 0, while $\Gamma_{in}$ is the boundary of the inlet where pressure is set to 1.</figcaption></figure>

The loss function is plotted for the 3 cases in Fig. 14. 10,000 points in the space-time domain are chosen randomly to assess the testing loss, while they are not used in the training.

Using the new adaptivity algorithm, the discrepancy between the training and testing loss is greatly reduced. That means that it offers a mean to reduce the generalization error thus provides better accuracy to the approximated fields.

## 5. Conclusion

In this study, Physics-Informed Neural Networks (PINN) are applied to solve two-phase flow problems in porous media. A novel residual-based adaptive algorithm is developed. The key is to utilize the PDE residual to build
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Left: Front function field. Right: pressure field. The plots are shown at times 0, 0.039 and 0.078 from top to bottom.</figcaption></figure>

a probability density function, from which new collocation points are drawn and added to the training set. The technique is applied to the different PDEs in the coupled system independently, thus, different collocations points are used for the different PDEs. Moreover, the technique is applied to enrich the points used to capture the initial and boundary conditions.
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Left: Front function field. Right: pressure field. The plots are shown at times 0, 0.039 and 0.078 from top to bottom.</figcaption></figure>

a probability density function, from which new collocation points are drawn and added to the training set. The technique is applied to the different PDEs in the coupled system independently, thus, different collocations points are used for the different PDEs. Moreover, the technique is applied to enrich the points used to capture the initial and boundary conditions.
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Front (fluid 1/fluid 2 interface) position in the *x*-direction with time for the fixed and adaptive cases along with analytical front position for the 2D case.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Front (fluid 1/fluid 2 interface) position in the *y*-direction with time for the fixed and adaptive cases along with analytical front position for the 2D case.</figcaption></figure>

The adaptivity algorithm is used to solve 1D and 2D two-phase flow in porous media. Using the new technique provided better results than using the classical PINN with fixed collocation points and also showed an improvement over the RAR technique. The adaptive technique can be seen as a form of regularization of the neural networks, thus reducing the generalization error.
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Front (fluid 1/fluid 2 interface) position in the *x*-direction with time for the fixed and adaptive cases along with analytical front position for the 2D case.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Front (fluid 1/fluid 2 interface) position in the *y*-direction with time for the fixed and adaptive cases along with analytical front position for the 2D case.</figcaption></figure>

The adaptivity algorithm is used to solve 1D and 2D two-phase flow in porous media. Using the new technique provided better results than using the classical PINN with fixed collocation points and also showed an improvement over the RAR technique. The adaptive technique can be seen as a form of regularization of the neural networks, thus reducing the generalization error.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Loss vs. Iteration graphs for fixed PINN (left plot), new adaptive PINN (middle plot), and PINN with RAR (right plot).</figcaption></figure>

## Acknowledgment

This study was funded under the PERFORM Thesis program of IRT Jules Verne, Bouguenais, France.

## References

[1] C.H. Park, L. Woo, Modeling void formation and unsaturated flow in liquid composite molding processes: a survey and review, J. Reinf. Plast. Compos. 30 (11) (2011) 957-977.

[2] K.M. Pillai, Modeling the unsaturated flow in liquid composite molding processes: a review and some thoughts, J. Compos. Mater. 38 (23) (2004) 2097-2118.

[3] L.A. Dillard, H.I. Essaid, W.N. Herkelrath, Multiphase flow modeling of a crude-oil spill site with a bimodal permeability distribution, Water Resour. Res. 33 (7) (1997) 1617-1632.

[4] M.G. Gerritsen, L.J. Durlofsky, Modeling fluid flow in oil reservoirs, Annu. Rev. Fluid Mech. 37 (2005) 211-238.

[5] L.M. Abriola, G.F. Pinder, A multiphase approach to the modeling of porous media contamination by organic compounds: 1. Equation development, Water Resour. Res. 21 (1) (1985) 11-18.

[6] M.B. Allen, Why upwinding is reasonable, in: Finite Elements in Water Resources, Springer, 1984, pp. 13-23.

[7] Z. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, 2006.

[8] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146-159.

[9] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201-225.

[10] L.-Q. Chen, Phase-field models for microstructure evolution, Annu. Rev. Mater. Res. 32 (1) (2002) 113-140.

[11] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375-389.

[12] Z.-B. Wang, R. Chen, H. Wang, Q. Liao, X. Zhu, S.-Z. Li, An overview of smoothed particle hydrodynamics for simulating multiphase flow, Appl. Math. Model. 40 (23-24) (2016) 9625-9655.

[13] M. Sawley, P. Cleary, J. Ha, Modelling of flow in porous media and resin transfer moulding using smoothed particle hydrodynamics, 1999.

[14] B. Liu, S. Bickerton, S.G. Advani, Modelling and simulation of resin transfer moulding (RTM)—gate control, venting and dry spot prediction, Composites A 27 (2) (1996) 135-141.

[15] M. Piller, G. Schena, M. Nolich, S. Favretto, F. Radaelli, E. Rossi, Analysis of hydraulic permeability in porous media: from high resolution X-ray tomography to direct numerical simulation, Transp. Porous Media 80 (1) (2009) 57-78.

[16] A.R. Valdez, B.M. Rocha, G. Chapiro, R.W. dos Santos, Uncertainty quantification and sensitivity analysis for relative permeability models of two-phase flow in porous media, J. Pet. Sci. Eng. 192 (2020) 107297.

[17] B.-Y. Kim, G.-J. Nam, H.-S. Ryu, J.-W. Lee, Optimization of filling process in RTM using genetic algorithm, Korea-Australia Rheol. J. 12 (1) (2000) 83-92.

[18] H. Fang, C. Gong, C. Li, X. Li, H. Su, L. Gu, A surrogate model based nested optimization framework for inverse problem considering interval uncertainty, Struct. Multidiscip. Optim. 58 (3) (2018) 869-883.

[19] B.Y. Mirghani, E.M. Zechman, R.S. Ranjithan, G. Mahinthakumar, Enhanced simulation-optimization approach using surrogate modeling for solving inverse problems, Environ. Forensics 13 (4) (2012) 348-363.

[20] S. Krogstad, A sparse basis POD for model reduction of multiphase compressible flow, in: SPE Reservoir Simulation Symposium, OnePetro, 2011.

[21] H. Bansal, S. Rave, L. Iapichino, W.H.A. Schilders, N. van de Wouw, Model order reduction framework for problems with moving discontinuities, in: ENUMATH, 2019, pp. 83-91.

[22] T.G. Kolda, B.W. Bader, Tensor decompositions and applications, SIAM Rev. 51 (3) (2009) 455-500.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Loss vs. Iteration graphs for fixed PINN (left plot), new adaptive PINN (middle plot), and PINN with RAR (right plot).</figcaption></figure>

## Acknowledgment

This study was funded under the PERFORM Thesis program of IRT Jules Verne, Bouguenais, France.

## References

[1] C.H. Park, L. Woo, Modeling void formation and unsaturated flow in liquid composite molding processes: a survey and review, J. Reinf. Plast. Compos. 30 (11) (2011) 957-977.

[2] K.M. Pillai, Modeling the unsaturated flow in liquid composite molding processes: a review and some thoughts, J. Compos. Mater. 38 (23) (2004) 2097-2118.

[3] L.A. Dillard, H.I. Essaid, W.N. Herkelrath, Multiphase flow modeling of a crude-oil spill site with a bimodal permeability distribution, Water Resour. Res. 33 (7) (1997) 1617-1632.

[4] M.G. Gerritsen, L.J. Durlofsky, Modeling fluid flow in oil reservoirs, Annu. Rev. Fluid Mech. 37 (2005) 211-238.

[5] L.M. Abriola, G.F. Pinder, A multiphase approach to the modeling of porous media contamination by organic compounds: 1. Equation development, Water Resour. Res. 21 (1) (1985) 11-18.

[6] M.B. Allen, Why upwinding is reasonable, in: Finite Elements in Water Resources, Springer, 1984, pp. 13-23.

[7] Z. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, 2006.

[8] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146-159.

[9] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201-225.

[10] L.-Q. Chen, Phase-field models for microstructure evolution, Annu. Rev. Mater. Res. 32 (1) (2002) 113-140.

[11] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375-389.

[12] Z.-B. Wang, R. Chen, H. Wang, Q. Liao, X. Zhu, S.-Z. Li, An overview of smoothed particle hydrodynamics for simulating multiphase flow, Appl. Math. Model. 40 (23-24) (2016) 9625-9655.

[13] M. Sawley, P. Cleary, J. Ha, Modelling of flow in porous media and resin transfer moulding using smoothed particle hydrodynamics, 1999.

[14] B. Liu, S. Bickerton, S.G. Advani, Modelling and simulation of resin transfer moulding (RTM)—gate control, venting and dry spot prediction, Composites A 27 (2) (1996) 135-141.

[15] M. Piller, G. Schena, M. Nolich, S. Favretto, F. Radaelli, E. Rossi, Analysis of hydraulic permeability in porous media: from high resolution X-ray tomography to direct numerical simulation, Transp. Porous Media 80 (1) (2009) 57-78.

[16] A.R. Valdez, B.M. Rocha, G. Chapiro, R.W. dos Santos, Uncertainty quantification and sensitivity analysis for relative permeability models of two-phase flow in porous media, J. Pet. Sci. Eng. 192 (2020) 107297.

[17] B.-Y. Kim, G.-J. Nam, H.-S. Ryu, J.-W. Lee, Optimization of filling process in RTM using genetic algorithm, Korea-Australia Rheol. J. 12 (1) (2000) 83-92.

[18] H. Fang, C. Gong, C. Li, X. Li, H. Su, L. Gu, A surrogate model based nested optimization framework for inverse problem considering interval uncertainty, Struct. Multidiscip. Optim. 58 (3) (2018) 869-883.

[19] B.Y. Mirghani, E.M. Zechman, R.S. Ranjithan, G. Mahinthakumar, Enhanced simulation-optimization approach using surrogate modeling for solving inverse problems, Environ. Forensics 13 (4) (2012) 348-363.

[20] S. Krogstad, A sparse basis POD for model reduction of multiphase compressible flow, in: SPE Reservoir Simulation Symposium, OnePetro, 2011.

[21] H. Bansal, S. Rave, L. Iapichino, W.H.A. Schilders, N. van de Wouw, Model order reduction framework for problems with moving discontinuities, in: ENUMATH, 2019, pp. 83-91.

[22] T.G. Kolda, B.W. Bader, Tensor decompositions and applications, SIAM Rev. 51 (3) (2009) 455-500.
[23] I.V. Oseledets, D. Savostianov, E.E. Tyrtyshnikov, Tucker dimensionality reduction of three-dimensional arrays in linear time, SIAM J. Matrix Anal. Appl. 30 (3) (2008) 939–956.

[24] J.V. Aguado, D. Borzacchiello, K.S. Kollepara, F. Chinesta, A. Huerta, Tensor representation of non-linear models using cross approximations, J. Sci. Comput. 81 (1) (2019) 22–47.

[25] D. Neron, P. Ladeveze, Idelsohns benchmark, Technical report, 2013.

[26] L. Peng, K. Mohseni, Nonlinear model reduction via a locally weighted POD method, Int. J. Numer. Metho. Eng. 106 (5) (2016) 372–396.

[27] A. Ferrero, A. Iollo, F. Larocca, Global and local POD models for the prediction of compressible flows with DG methods, Internat. J. Numer. Methods Engrg. 116 (5) (2018) 332–357.

[28] M. Raissi, P. Perdikaris, G.E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686–707.

[29] A.D. Jagtap, G.E. Karniadakis, Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations, Commun. Comput. Phys. 28 (5) (2020) 2002–2041.

[30] A.D. Jagtap, E. Kharazmi, G.E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Comput. Methods Appl. Mech. Engrg. 365 (2020) 113028.

[31] Z. Mao, A.D. Jagtap, G.E. Karniadakis, Physics-informed neural networks for high-speed flows, Comput. Methods Appl. Mech. Engrg. 360 (2020) 112789.

[32] G. Pang, L. Lu, G.E. Karniadakis, Fpinns: Fractional physics-informed neural networks, SIAM J. Sci. Comput. 41 (4) (2019) A2603–A2626.

[33] K. Shukla, P.C. Di Leoni, J. Blackshire, D. Sparkman, G.E. Karniadakis, Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks, J. Nondestruct. Eval. 39 (3) (2020) 1–20.

[34] Q. Zhu, Z. Liu, J. Yan, Machine learning for metal additive manufacturing: predicting temperature and melt pool fluid dynamics using physics-informed neural networks, Comput. Mech. 67 (2) (2021) 619–635.

[35] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes, A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics, Comput. Methods Appl. Mech. Engrg. 379 (2021) 113741.

[36] S.A. Niaki, E. Haghighat, T. Campbell, A. Poursartip, R. Vaziri, Physics-informed neural network for modelling the thermochemical curing process of composite-tool systems during manufacture, Comput. Methods Appl. Mech. Engrg. 384 (2021) 113959.

[37] L. Sun, H. Gao, S. Pan, J.-X. Wang, Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data, Comput. Methods Appl. Mech. Engrg. 361 (2020) 112732.

[38] A.M. Tartakovsky, C.O. Marrero, P. Perdikaris, G.D. Tartakovsky, D. Barajas-Solano, Physics-informed deep neural networks for learning parameters and constitutive relationships in subsurface flow problems, Water Resour. Res. 56 (5) (2020) e2019WR026731.

[39] S. Goswami, C. Anitescu, S. Chakraborty, T. Rabczuk, Transfer learning enhanced physics informed neural network for phase-field modeling of fracture, Theor. Appl. Mech. 106 (2020) 102447.

[40] M.J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53 (3) (1984) 484–512.

[41] P. Díez, A. Huerta, A unified approach to remeshing strategies for finite element h-adaptivity, Comput. Methods Appl. Mech. Engrg. 176 (1–4) (1999) 215–229.

[42] H. Askes, A. Rodríguez-Ferran, A combined rh-adaptive scheme based on domain subdivision. Formulation and linear examples, Internat. J. Numer. Methods Engrg. 51 (3) (2001) 253–273.

[43] H. Askes, A. Rodríguez-Ferran, A combined rh-adaptive scheme based on domain subdivision. Formulation and linear examples, Internat. J. Numer. Methods Engrg. 51 (3) (2001) 253–273.

[44] I. Babuška, M. Suri, The p and h-p versions of the finite element method, basic principles and properties, SIAM Rev. 36 (4) (1994) 578–632.

[45] X.D. Zhang, J.-Y. Trépanier, R. Camarero, A posteriori error estimation for finite-volume solutions of hyperbolic conservation laws, Comput. Methods Appl. Mech. Engrg. 185 (1) (2000) 1–19.

[46] D. Lee, Y. Tsuei, A formula for estimation of truncation errors of convection terms in a curvilinear coordinate system, J. Comput. Phys. 98 (1) (1992) 90–100.

[47] D.S. McRae, R-refinement grid adaptation algorithms and issues, Comput. Methods Appl. Mech. Engrg. 189 (4) (2000) 1161–1182.

[48] A. Choudhary, C. Roy, Efficient residual-based mesh adaptation for 1D and 2D CFD applications, in: 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2011, p. 214.

[49] C. Roy, Strategies for driving mesh adaptation in CFD, in: 47th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2009, p. 1302.

[50] R.P. Dwight, Heuristic a posteriori estimation of error due to dissipation in finite volume schemes and application to mesh adaptation, J. Comput. Phys. 227 (5) (2008) 2845–2863.

[51] L. Lu, X. Meng, Z. Mao, G.E. Karniadakis, Deepxde: A deep learning library for solving differential equations, SIAM Rev. 63 (1) (2021) 208–228.

[52] T. Kadeethum, T.M. Jørgensen, H.M. Nick, Physics-informed neural networks for solving nonlinear diffusivity and Biot's equations, PLoS One 15 (5) (2020) e0232683.

[53] J. Weitzenböck, R. Shenoi, P. Wilson, Radial flow permeability measurement. Part a: Theory, Composites A 30 (6) (1999) 781–796.
[23] I.V. Oseledets, D. Savostianov, E.E. Tyrtyshnikov, Tucker dimensionality reduction of three-dimensional arrays in linear time, SIAM J. Matrix Anal. Appl. 30 (3) (2008) 939-956.

[24] J.V. Aguado, D. Borzacchiello, K.S. Kollepara, F. Chinesta, A. Huerta, Tensor representation of non-linear models using cross approximations, J. Sci. Comput. 81 (1) (2019) 22-47.

[25] D. Neron, P. Ladeveze, Idelsohns benchmark, Technical report, 2013.

[26] L. Peng, K. Mohseni, Nonlinear model reduction via a locally weighted POD method, Int. J. Numer. Metho. Eng. 106 (5) (2016) 372-396.

[27] A. Ferrero, A. Iollo, F. Larocca, Global and local POD models for the prediction of compressible flows with DG methods, Internat. J. Numer. Methods Engrg. 116 (5) (2018) 332-357.

[28] M. Raissi, P. Perdikaris, G.E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686-707.

[29] A.D. Jagtap, G.E. Karniadakis, Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations, Commun. Comput. Phys. 28 (5) (2020) 2002-2041.

[30] A.D. Jagtap, E. Kharazmi, G.E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Comput. Methods Appl. Mech. Engrg. 365 (2020) 113028.

[31] Z. Mao, A.D. Jagtap, G.E. Karniadakis, Physics-informed neural networks for high-speed flows, Comput. Methods Appl. Mech. Engrg. 360 (2020) 112789.

[32] G. Pang, L. Lu, G.E. Karniadakis, Fpinns: Fractional physics-informed neural networks, SIAM J. Sci. Comput. 41 (4) (2019) A2603-A2626.

[33] K. Shukla, P.C. Di Leoni, J. Blackshire, D. Sparkman, G.E. Karniadakis, Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks, J. Nondestruct. Eval. 39 (3) (2020) 1-20.

[34] Q. Zhu, Z. Liu, J. Yan, Machine learning for metal additive manufacturing: predicting temperature and melt pool fluid dynamics using physics-informed neural networks, Comput. Mech. 67 (2) (2021) 619-635.

[35] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes, A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics, Comput. Methods Appl. Mech. Engrg. 379 (2021) 113741.

[36] S.A. Niaki, E. Haghighat, T. Campbell, A. Poursartip, R. Vaziri, Physics-informed neural network for modelling the thermochemical curing process of composite-tool systems during manufacture, Comput. Methods Appl. Mech. Engrg. 384 (2021) 113959.

[37] L. Sun, H. Gao, S. Pan, J.-X. Wang, Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data, Comput. Methods Appl. Mech. Engrg. 361 (2020) 112732.

[38] A.M. Tartakovsky, C.O. Marrero, P. Perdikaris, G.D. Tartakovsky, D. Barajas-Solano, Physics-informed deep neural networks for learning parameters and constitutive relationships in subsurface flow problems, Water Resour. Res. 56 (5) (2020) e2019WR026731.

[39] S. Goswami, C. Anitescu, S. Chakraborty, T. Rabczuk, Transfer learning enhanced physics informed neural network for phase-field modeling of fracture, Theor. Appl. Mech. 106 (2020) 102447.

[40] M.J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53 (3) (1984) 484-512.

[41] P. Díez, A. Huerta, A unified approach to remeshing strategies for finite element h-adaptivity, Comput. Methods Appl. Mech. Engrg. 176 (1-4) (1999) 215-229.

[42] H. Askes, A. Rodríguez-Ferran, A combined rh-adaptive scheme based on domain subdivision. Formulation and linear examples, Internat. J. Numer. Methods Engrg. 51 (3) (2001) 253-273.

[43] H. Askes, A. Rodríguez-Ferran, A combined rh-adaptive scheme based on domain subdivision. Formulation and linear examples, Internat. J. Numer. Methods Engrg. 51 (3) (2001) 253-273.

[44] I. Babuška, M. Suri, The p and h-p versions of the finite element method, basic principles and properties, SIAM Rev. 36 (4) (1994) 578-632.

[45] X.D. Zhang, J.-Y. Trépanier, R. Camarero, A posteriori error estimation for finite-volume solutions of hyperbolic conservation laws, Comput. Methods Appl. Mech. Engrg. 185 (1) (2000) 1-19.

[46] D. Lee, Y. Tsuei, A formula for estimation of truncation errors of convection terms in a curvilinear coordinate system, J. Comput. Phys. 98 (1) (1992) 90-100.

[47] D.S. McRae, R-refinement grid adaptation algorithms and issues, Comput. Methods Appl. Mech. Engrg. 189 (4) (2000) 1161-1182.

[48] A. Choudhary, C. Roy, Efficient residual-based mesh adaptation for 1D and 2D CFD applications, in: 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2011, p. 214.

[49] C. Roy, Strategies for driving mesh adaptation in CFD, in: 47th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2009, p. 1302.

[50] R.P. Dwight, Heuristic a posteriori estimation of error due to dissipation in finite volume schemes and application to mesh adaptation, J. Comput. Phys. 227 (5) (2008) 2845-2863.

[51] L. Lu, X. Meng, Z. Mao, G.E. Karniadakis, Deepxde: A deep learning library for solving differential equations, SIAM Rev. 63 (1) (2021) 208-228.

[52] T. Kadeethum, T.M. Jørgensen, H.M. Nick, Physics-informed neural networks for solving nonlinear diffusivity and Biot's equations, PLoS One 15 (5) (2020) e0232683.

[53] J. Weitzenböck, R. Shenoi, P. Wilson, Radial flow permeability measurement. Part a: Theory, Composites A 30 (6) (1999) 781-796.