EPINN-NSE: Enhanced Physics-Informed Neural Networks for Solving
  Navier-Stokes Equations
Fluid mechanics is a fundamental field in engineering and science. Solving
the Navier-Stokes equation (NSE) is critical for understanding the behavior of
fluids. However, the NSE is a complex partial differential equation that is
difficult to solve, and classical numerical methods can be computationally
expensive. In this paper, we present an innovative approach for solving the NSE
using Physics Informed Neural Networks (PINN) and several novel techniques that
improve their performance. The first model is based on an assumption that
involves approximating the velocity component by employing the derivative of a
stream function. This assumption serves to simplify the system and guarantees
that the velocity adheres to the divergence-free equation. We also developed a
second more flexible model that approximates the solution without any
assumptions. The proposed models can effectively solve two-dimensional NSE.
Moreover, we successfully applied the second model to solve the
three-dimensional NSE. The results show that the models can efficiently and
accurately solve the NSE in three dimensions. These approaches offer several
advantages, including high trainability, flexibility, and efficiency.
Ayoub Farkane, Mounir Ghogho, Mustapha Oudani, Mohamed Boutayeb
# Epinn-NSE: Enhanced Physics-Informed Neural Networks for Solving Navier-Stokes Equations1

###### Abstract

Fluid mechanics is a fundamental field in engineering and science. Solving the Navier-Stokes equation (NSE) is critical for understanding the behavior of fluids. However, the NSE is a complex partial differential equation that is difficult to solve, and classical numerical methods can be computationally expensive. In this paper, we present an innovative approach for solving the NSE using Physics Informed Neural Networks (PINN) and several novel techniques that improve their performance. The first model is based on an assumption that involves approximating the velocity component by employing the derivative of a stream function. This assumption serves to simplify the system and guarantees that the velocity adheres to the divergence-free equation. We also developed a second more flexible model that approximates the solution without any assumptions. The proposed models can effectively solve two-dimensional NSE. Moreover, we successfully applied the second model to solve the three-dimensional NSE. The results show that the models can efficiently and accurately solve the NSE in three dimensions. These approaches offer several advantages, including high trainability, flexibility, and efficiency.

N : 35Q35, 65M99, 68T05

## 1 Introduction

Fluid flow has a wide range of applications in mechanical and chemical engineering, biological systems, and astrophysics. Essentially, all moving bodies are connected to fluid flow, such as vehicles, airplanes, trucks, trains, birds, dolphins, and so on.

Navier-Stokes equations are Partial Differential Equations (PDEs) that are used to describe the flow of fluids. As with many PDEs, the NSE has no analytical solution, and even in three dimensions, it remains one of the Millennium Prize problems. Classical numerical methods, such as the finite difference method and finite element method (FEM), can produce an approximate solution for the NSE. However, deep learning techniques have recently become popular for solving PDEs. One such method is the physics informed neural network (PINN) [22, 19] which uses a deep neural network (DNN) based on optimization problems or residual loss functions to solve a PDE. Other deep learning techniques, such as the deep Galerkin method (DGM)[25] have also been proposed in the literature for solving PDEs. The DGM is particularly useful for solving high-dimensional free boundary PDEs such as the High-dimensional Hamilton-Jacobi-Bellman. However, DGM's numerical performance for other types of PDEs (elliptic, hyperbolic, and partial-integral differential equations, etc.) remains to be investigated. Different assumptions are considered on the operator of PDEs to guarantee the convergence of this neural network. Another approach to approximating PDEs using a variational form including DNNs is the deep Ritz method (DRM) [27], which reformulates the PDE as an optimization problem with a variational formulation. The DRM is naturally adaptive, less sensitive to the problem dimensionality, and works well with the stochastic gradient descent. Nonetheless, there are a variety
of drawbacks to DRM, such as the problem of convexity, and the issue of local minima and saddle points. Moreover, the DRM has no consistent conclusion about the convergence rate and the consideration of the necessary boundary condition is more complicated than in previous approaches. Additionally, there are some concerns with the network structure and activation function.

The work [5] solves a high dimensional semi-linear PDE with a condition on diffusion generator of the equation. To accomplish this, the method used the solution of a backward stochastic differential equation (BSDE) established in [16], and a time discretization scheme proposed in [2]. The approach connected the PDE and stochastic differential equation by using the deep BSDE. However, the numerical results showed that the method did not converge to a standard network of fully connected layers. Other works are based on the connection between PDEs and stochastic processes through the application of the Feynman-Kac formula [10]. This paper [1] solves a Kolmogorov PDE using deep learning. The Feynman-Kac formula may be extended to the full class of high linear parabolic equations.

Chebyshev neural networks are suggested in [14] to solve two-dimensional (2D) linear PDEs. The approach involves using a single hidden layer feedforward neural network and approximating the solution with the Kronecker product of two Chebyshev polynomials. However, this approach is still far from being suitable for solving nonlinear PDEs. In some relevant papers [20, 21], the nonlinear terms of PDEs were locally linearized in the temporal domain. However, this approach is only suitable for discrete-time domains and may reduce the accuracy of predictions in highly nonlinear systems. In another approach called MAgNet, proposed in [3], a graph neural network was used to forecast a spatially continuous solution of a PDE given a spatial position query. However, MAgNet uses a first-order explicit time-stepping scheme that can be unstable.

It is widely acknowledged that solving nonlinear-dynamic PDEs like the NSE is a challenging task, especially in high dimensions. The interesting paper [22] presented a solution to the 2D NSE using the PINN. In this study, we propose an improved model-based neural network for solving the NSE and compare it to the PINN [22] approach. We also conduct several experiments involving changes in data sizes and hyperparameters. Finally, we apply these enhanced models to solve the three-dimensional (3D) NSE using a test solution.

## 2 Governing Equations

### 2D NSE for Incompressible Flows

The NSE can be used to model various fluid flows [6] such as Couette-Poiseuille flows, Beltrami flows, pipe and cylinder flow, and oscillating plates. Consider the general form of the 2D NSE for an incompressible fluid:

\[\frac{\partial u}{\partial t}+\beta\left(u\frac{\partial u}{ \partial x}+v\frac{\partial u}{\partial y}\right) =-\frac{\partial p}{\partial x}+\nu\left(\frac{\partial^{2}u}{ \partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right) \tag{1}\] \[\frac{\partial v}{\partial t}+\beta\left(u\frac{\partial v}{ \partial x}+v\frac{\partial v}{\partial y}\right) =-\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^{2}v}{ \partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right), \tag{2}\]

where \(u\) is the x-component of velocity, \(v\) is the y-component of velocity, \(t\) is time, \(p\) is the flow pressure, \(\nu\) is the fluid viscosity1 and \(\beta\) is a fixed constant.
The solution of 2D NSE must satisfy the incompressibility equation, which requires the fluid velocity to be a divergence-free function:

\[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 \tag{1}\]

The previous equation describes the mass conservation of fluid.

When dealing with fluid dynamics, it is necessary to define both the initial and boundary conditions in order to fully specify a boundary value problem. The type of boundary condition used depends on the particular fluid phenomenon being modeled. In some cases, the boundary conditions may be Dirichlet boundary conditions, which specify the values of the solution at the boundaries. In other cases, the boundary conditions may be Neumann boundary conditions, which specify the derivatives of the solution at the borders. It is also possible to use mixed boundary conditions, which combine elements of Dirichlet and Neumann boundary conditions. The selection of the appropriate boundary condition is dependent on the specific problem under consideration.

### 3d Nse

The problems (1)-(3) can be extended mathematically to three dimensions using the following form:

\[\frac{\partial u}{\partial t}+\beta\left(u\frac{\partial u}{ \partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right) =-\frac{\partial p}{\partial x}+\nu\left(\frac{\partial^{2}u}{ \partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{ \partial z^{2}}\right) \tag{2}\] \[\frac{\partial v}{\partial t}+\beta\left(u\frac{\partial v}{ \partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right) =-\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^{2}v}{ \partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right)+\frac{\partial^{2 }v}{\partial z^{2}}\] (3) \[\frac{\partial w}{\partial t}+\beta\left(u\frac{\partial w}{ \partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right) =-\frac{\partial p}{\partial z}+\nu\left(\frac{\partial^{2}w}{ \partial x^{2}}+\frac{\partial^{2}w}{\partial y^{2}}\right)+\frac{\partial^{2 }w}{\partial z^{2}} \tag{4}\]

where the velocity components are \(u\), \(v\), and \(w\). The other parameters are defined as previously. Equation (4) introduced the conservation of momentum equation for the \(w\) velocity component. As previously stated, the solution can be sought within the divergence-free set:

\[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w} {\partial z}=0 \tag{5}\]

Solving this type of PDE can be challenging, but recent advancements in deep learning techniques have demonstrated potential in addressing these problems. One such approach is the PINN approach, which is elaborated on in the next section.

## 3 PINN model

The PINN is a deep learning algorithm that has been designed to solve PDEs. It combines the flexibility and expressiveness of neural networks with the ability to incorporate physical laws and constraints into the model.

### Methodology

The PINN approach consists of two main parts. The first part involves generating a candidate solution, referred to as "\(u\)", at predefined points using a multi-layer perception (MLP). The second part consists of calculating the terms associated with the PDE for the candidate solution and constructing a loss function that is optimized using a gradient descent algorithm. We will examine in greater detail how the methods work using the example below. Consider the following
one-dimensional PDE:

\[\frac{\partial u}{\partial t}+\frac{\partial^{2}u}{\partial x^{2}} =0,\ \ \forall x\in\Omega,\,\forall t\in[0,T]. \tag{1}\] \[u= 0\ \ \text{in}\ \Gamma_{\Omega},\,\forall t\in[0,T]\] (2) \[u(t=0) =u_{0},\ \ \forall x\in\Omega \tag{3}\]

In order to apply the PINN to solve the equation shown above, the neural network takes in \(x\) and \(t\) as inputs and produces a candidate \(u\) as the output. The solution is found through minimizing the following loss function:

\[\text{MSE}_{u}=\text{MSE}_{0}+\text{MSE}_{b}+\text{MSE}_{f}, \tag{4}\]

where:

\[\text{MSE}_{0} =\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\left|u\left(0,x_{0}^{i}\right) -u_{0}^{i}\right|^{2} \tag{5}\] \[\text{MSE}_{b} =\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\left(\left|u^{i}\left(t_{b}^{ i},.\right)-u_{b}^{i}\right|^{2}\right.\] (6) \[\text{MSE}_{f} =\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}\left|f\left(t_{f}^{i},x_{f}^{ i}\right)\right|^{2} \tag{7}\]

The set of data points that represent the initial conditions of the problem are given by \(\left\{x_{0}^{i},u_{0}^{i}\right\}_{i=1}^{N_{0}}\), where \(N_{0}\) is the number of data points; the collocation points on the boundary are represented by the set \(\left\{t_{b}^{i}\right\}_{i=1}^{N_{b}}\), where \(N_{b}\) is the number of points; and the set \(\left\{t_{f}^{i},x_{f}^{i}\right\}_{i=1}^{N_{f}}\) denotes the collocation points on \(f(t,x)\).

The mean squared error \(\text{MSE}_{u}\) is a measurement that calculates the average of the squared differences between the predicted and actual values. In the current context, \(\text{MSE}_{0}\) represents the MSE metric specifically associated with the initial conditions of the problem. Similarly, \(\text{MSE}_{b}\) is the MSE metric associated with the boundary conditions, and \(\text{MSE}_{f}\) is a loss function that ensures the solution satisfies the equation within the domain. The function \(f\) is exactly:

\[f(t,x)=\frac{\partial u}{\partial t}+\frac{\partial^{2}u}{\partial x^{2}} \tag{8}\]

The combination of these three loss functions plays a crucial role in obtaining an accurate solution for the problem at hand. This approach can be extended to different types of PDEs, with some exceptions for more complex problems where modifications to the loss function and neural network architectures may be necessary.

### PINN to Solve 2D NSE

A notable study [19] employed the PINN approach to solving the 2D NSE. The authors used an assumption to simplify the equation of incompressibility (3) by using a stream function \(\varphi\) to approximate the \(u\) and \(v\) components:

\[u=\frac{\partial\varphi}{\partial y},\ \ \ v=-\frac{\partial\varphi}{\partial x} \tag{9}\]

The stream function \(\varphi\) is calculated as the neural network's output, such as pressure \(p\). Consider the two functions \(f_{1}\) and \(f_{2}\) such that:
\[f_{1}(t,x,y) = \frac{\partial u}{\partial t}+\beta\left(u\frac{\partial u}{ \partial x}+v\frac{\partial u}{\partial y}\right)+\frac{\partial p}{\partial x }-\nu\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y ^{2}}\right) \tag{11}\] \[f_{2}(t,x,y) = \frac{\partial v}{\partial t}+\beta\left(u\frac{\partial v}{ \partial x}+v\frac{\partial v}{\partial y}\right)+\frac{\partial p}{\partial y }-\nu\left(\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y ^{2}}\right) \tag{10}\]

Using the previously defined functions (10) and (11), the model is trained with the following loss functions:

\[\text{MSE} = \frac{1}{N}\sum_{i=1}^{N}\left(\left|u\left(t^{i},x^{i},y^{i} \right)-u^{i}\right|^{2}+\left|v\left(t^{i},x^{i},y^{i}\right)-v^{i}\right|^{ 2}\right) \tag{12}\] \[\quad+\frac{1}{N}\sum_{i=1}^{N}\left(\left|f_{1}\left(t^{i},x^{ i},y^{i}\right)\right|^{2}+\left|f_{2}\left(t^{i},x^{i},y^{i}\right)\right|^{2} \right), \tag{13}\]

where the \(\left\{t^{i},x^{i},y^{i},u^{i},v^{i}\right\}\) represents accurate data obtained from a PDE simulation or real-world measurements.

The model is trained on data from experiments involving the flow of an incompressible fluid around a circular cylinder. The data is generated using a numerical solver called Nektar++, which discretizes the solution domain into triangular elements and approximates the solution as a linear combination of tenth-order hierarchical, semi-orthogonal Jacobi polynomial expansions. The method described in this section was used to solve the 2D NSE in this work [22]. While the PINN approach has proven advantageous in addressing this complex problem, current research is focused on refining the model and improving its effectiveness. The next section provides a comprehensive overview of the enhanced models, which aim to improve the methodology used to solve the NSE.

Figure 1: Method with the assumption for finding solutions to the 2D NSE using the PINN. The domain is represented by variables \(x\) and \(y\), while \(t\) represents time. The stream function is represented by \(\varphi\), and its partial derivatives with respect to \(x\) and \(y\) are represented by \(\varphi_{x}\) and \(\varphi_{y}\), respectively. “Velocity” involves extracting a candidate solution velocity using the assumption. “MLP” is a multi-layer perceptron. “\(L_{1}\), \(L_{2}\)” are the derivative operators associated with the equation.
## 4 EPINN-NSE Models

The PINN method is based on deep learning techniques and has shown promise in solving various PDEs. However, further research is necessary to improve its application to more complex problems. There are numerous papers in the literature that focus on improving PINN and its various applications. One of the significant challenges that PINN faces is the issue of training time. While the training duration is not a critical concern for simple cases, it can be a significant issue for complex equations, with training potentially taking days rather than hours. In this context, [24] proposed a technique called discretely trained PINNs (DTPINNs) to accelerate the approach by discretizing the PDE and replacing the exact spatial derivatives with high-order numerical discretizations using the finite differences method. Nonetheless, this approach still has a long way to go before it can be applied to problems that involve continuous time.

The PINN approach lacks a theoretical underpinning that ensures the convergence of the model and may be prone to failing to find the appropriate solution. For instance, the work [13] highlighted some instances where PINN fails to find an accurate solution for certain PDE cases, such as the reaction-diffusion and convection problems. The paper discussed various situations where these networks may fail, such as when the network fails to accurately capture the underlying physics of the system being simulated or when the network is not properly trained. Furthermore, the framework in [13] also provides potential solutions to these issues, such as incorporating more physics information into the network or using different training methods.

Solving the 2D NSE presents additional challenges for PINN. For instance, the model introduced in [22] required an extended training duration (more than 8 hours). The assumption (3.9) is based on the derivative relation between the stream function \(\varphi\) and the velocity components respectively \(u\), and \(v\). The factor of learning the pressure without training leads to a higher probability that the model produces inaccurate solutions. Another issue related to neural network learning is that the PINN approach uses full

Figure 2: Model of PINN without using any assumption for solving the 2D NSE.
batch learning, which may increase the probability of getting stuck on saddle points or local minima. Using full batch learning for training the model may present challenges when working with a large dataset, which can negatively impact the results.

This paper proposes improvements to the previously introduced PINN model for solving PDEs. Specifically, we enhance the model described in Section 3.2 by incorporating training on pressure and using advanced deep learning techniques. This involves introducing a pressure loss function to equations (3.12) and (3.13):

\[\frac{1}{N}\sum_{i=1}^{N}\left(\left|p\left(t^{i},x^{i},y^{i}\right)-p^{i} \right|^{2}\right), \tag{4.1}\]

where the \(p\left(t^{i},x^{i},y^{i}\right)\) represents the pressure candidate and \(p^{i}\) pressure obtained from data. The aim of incorporating training on pressure is to enhance the pace and accuracy of learning and to improve the results of approximating pressure.

In addition to those improvements, we develop a novel model that directly approximates the velocity components and pressure without using the assumption (3.9). Canceling this assumption based on the derivatives of the stream function can improve the training time and the overall results. The new model employs a neural network that takes space and time variables as inputs and produces solutions for velocity and pressure. The governing equations are computed using automatic differentiation to determine the terms of the loss function that the network will minimize.

Furthermore, our models incorporate technical improvements in neural network architecture and optimize hyperparameters. Several studies [11, 9, 17] have highlighted the advantages of using mini-batch gradient descent. To take advantage of these benefits, we use mini-batch learning in our approach. This enables us to handle large datasets and make use of parallel processing, leading to faster convergence and lower memory requirements. We employed the Adam optimizer [12], which is less sensitive to the initialization of the model parameters. Compared to other optimization algorithms, Adam is less likely to get stuck in local minima or saddle points. Adam adjusts the learning rate for each parameter based on its past gradients, resulting in faster convergence and improved optimization. To effectively train the models and find the most-effective one, we use the validation loss function to evaluate the models at each epoch and compare the results with previous loss values. Afterward, we save the model that performs best. Moreover, we integrated an early stopping criterion or "patience." This technique seeks to terminate the model training process when it fails to exhibit significant improvement after a predetermined number of iterations.

### Solving the 2D NSE

For this case, we aim to employ both models to solve the problem. The first model based on assumption (3.9) is explained in Figure 1 with adding the techniques introduced previously and the loss function associated with pressure (4.1). In the second approach (Figure 2), the neural network generates a candidate solution \(u\), \(v\), and \(p\). This solution is then evaluated using a loss function. Removing the assumption (3.9) requires the addition of the loss function linked to the incompressibility equation (2.7). The neural network is trained to find the solution and learn the functions \(f_{1}(t,x,y)\) and \(f_{2}(t,x,y)\) by minimizing the following loss
function:

\[\text{MSE} = \text{MSE}_{e}+\text{MSE}_{f}+\text{MSE}_{d} \tag{10}\] \[\text{MSE}_{e} = \frac{1}{N}\sum_{i=1}^{N}\left(\left|u\left(t^{i},x^{i},y^{i} \right)-u^{i}\right|^{2}+\left|v\left(t^{i},x^{i},y^{i}\right)-v^{i}\right|^{2}\right)\] (11) \[\qquad+\frac{1}{N}\sum_{i=1}^{N}\left(\left|p\left(t^{i},x^{i},y ^{i}\right)-p^{i}\right|^{2}\right)\] \[\text{MSE}_{f} = \frac{1}{N}\sum_{i=1}^{N}\left(\left|f_{1}\left(t^{i},x^{i},y^{i} \right)\right|^{2}+\left|f_{2}\left(t^{i},x^{i},y^{i}\right)\right|^{2}\right)\] (12) \[\text{MSE}_{d} = \frac{1}{N}\sum_{i=1}^{N}\left(\left|u_{x}\left(t^{i},x^{i},y^{i} \right)+v_{y}\left(t^{i},x^{i},y^{i}\right)\right|^{2}\right),\]

where \(N\) is the number of collocation points denoted as \(\left\{t^{i},\ x^{i},\ y^{i}\right\}\) and \(u(t^{i},x^{i},y^{i})\), \(v(t^{i},x^{i},y^{i})\) and \(p(t^{i},x^{i},y^{i})\) are the candidate solutions being evaluated. \(u^{i}\), \(v^{i}\), and \(p^{i}\) are the real solutions obtained from a PDE solver or real-world observations.

The loss function specified in equation (10) is employed to measure the difference between the exact and approximate solutions. Additionally, the loss function defined in equation (11) ensures that the solution adheres to the equation, while the loss function outlined in equation (12) represents the divergence-free equation.

### Solving the 3D NSE

The 3D NSE is one of the most difficult PDEs to solve. Only a few studies have used deep learning techniques to tackle this equation, particularly in a 3D problem. Here we cite some papers [15, 23] that have been conducted using numerical methods such as the FEM and finite volume method (FVM) integrated with DNNs. Other works [7, 26, 18] addressed specific cases of the NSE in fluid mechanics problems using PINNs. In this section, we adapted the \(I_{2}\) model to approximate the solution to the 3D problem.

Consider the 3D NSE:

\[\frac{\partial u}{\partial t}+\beta\left(u\frac{\partial u}{ \partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z} \right)+\frac{\partial p}{\partial x}-\nu\left(\frac{\partial^{2}u}{ \partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{ \partial z^{2}}\right) = h(t,x,y,z) \tag{13}\] \[\frac{\partial v}{\partial t}+\beta\left(u\frac{\partial v}{ \partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z} \right)+\frac{\partial p}{\partial y}-\nu\left(\frac{\partial^{2}v}{\partial x ^{2}}+\frac{\partial^{2}v}{\partial y^{2}}+\frac{\partial^{2}v}{\partial z^{2} }\right) = g(t,x,y,z)\] (14) \[\frac{\partial w}{\partial t}+\beta\left(u\frac{\partial w}{ \partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z} \right)+\frac{\partial p}{\partial z}-\nu\left(\frac{\partial^{2}w}{\partial x ^{2}}+\frac{\partial^{2}w}{\partial y^{2}}+\frac{\partial^{2}w}{\partial z^{2} }\right) = k(t,x,y,z),\]

where \((t,x,y,z)\in[0,T]\times\Omega\), such that \(\Omega\subset\mathbb{R}^{3}\) (we assume \(\Omega=[0,1]^{3}\), \(T=20\)). The functions \(h\), \(g\), and \(k\) represent the force components applied to the system under consideration, While the other variables are as defined in Section 2. The solution sought satisfies:

\[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+ \frac{\partial w}{\partial z} = I(t,x,y,z) \tag{15}\]

The homogeneous Dirichlet boundary condition is considered for this problem:

\[u(t,x,y,z)=v(t,x,y,z)=w(t,x,y,z)=0,\ \ \forall(x,y,z)\in\Gamma_{\Omega},\forall t \in[0,T], \tag{16}\]
where \(\Gamma_{\Omega}\) is the boundary of domain \(\Omega\).

For 3D NSE, as illustrated in Figure 3, the neural network inputs consist of the spatial variables \(x,y,z\) and the temporal variable \(t\). The neural network outputs comprise the three velocity components, represented by \(u\), \(v\), and \(w\), as well as the pressure \(p\). Additionally, a similar approach is used in constructing the loss function as the 2D case. In this case, the loss functions are represented as:

\[\text{MSE} =\text{MSE}_{e}+\text{MSE}_{f}+\text{MSE}_{d} \tag{4.11}\] \[\text{MSE}_{e} =\frac{1}{N}\sum_{i=1}^{N}\left(\left|u\left(t^{i},x^{i},y^{i},z^ {i}\right)-u^{i}\right|^{2}+\left|v\left(t^{i},x^{i},y^{i},z^{i}\right)-v^{i} \right|^{2}\right.\] (4.12) \[\qquad+\left.\left|w\left(t^{i},x^{i},y^{i},z^{i}\right)-w^{i} \right|^{2}\right)+\frac{1}{N}\sum_{i=1}^{N}\left(\left|p\left(t^{i},x^{i},y^ {i},z^{i}\right)-p^{i}\right|^{2}\right)\] (4.13) \[\text{MSE}_{f} =\frac{1}{N}\sum_{i=1}^{N}\left(\left|f_{1}\left(t^{i},x^{i},y^{ i},z^{i}\right)\right|^{2}+\left|f_{2}\left(t^{i},x^{i},y^{i},z^{i}\right) \right|^{2}+\left|f_{3}\left(t^{i},x^{i},y^{i},z^{i}\right)\right|^{2}\right)\] (4.14) \[\text{MSE}_{d} =\frac{1}{N}\sum_{i=1}^{N}\left(\left|u_{x}(t^{i},x^{i},y^{i},z^ {i})+v_{y}(t^{i},x^{i},y^{i},z^{i})+w_{z}(t^{i},x^{i},y^{i},z^{i})\right.\] \[\qquad\left.-I(t^{i},x^{i},y^{i},z^{i})\right|^{2}\right)\]

The neural network generates outputs identified as \(u\left(t^{i},x^{i},y^{i},z^{i}\right)\), \(v\left(t^{i},x^{i},y^{i},z^{i}\right)\), \(w\left(t^{i},x^{i},y^{i},z^{i}\right)\), and \(p\left(t^{i},x^{i},y^{i},z^{i}\right)\). On the other hand, \(u^{i}\), \(v^{i}\), \(w^{i}\), and \(p^{i}\) represent the true solution obtained from the data. Additionally, \(f_{1}\), \(f_{2}\), and \(f_{3}\) are functions defined as follows:

\[f_{1}(t,x,y,z) =\frac{\partial u}{\partial t}+\beta\left(u\frac{\partial u}{ \partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z} \right)+\frac{\partial p}{\partial x}-\nu\left(\frac{\partial^{2}u}{\partial x ^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}}\right) \tag{4.15}\] \[\qquad-h(t,x,y,z)\] (4.16) \[f_{2}(t,x,y,z) =\frac{\partial v}{\partial t}+\beta\left(u\frac{\partial v}{ \partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z} \right)+\frac{\partial p}{\partial y}-\nu\left(\frac{\partial^{2}v}{\partial x ^{2}}+\frac{\partial^{2}v}{\partial y^{2}}+\frac{\partial^{2}v}{\partial z^{2}}\right)\] (4.17) \[\qquad-g(t,x,y,z)\] \[f_{3}(t,x,y,z) =\frac{\partial w}{\partial t}+\beta\left(u\frac{\partial w}{ \partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z} \right)+\frac{\partial p}{\partial z}-\nu\left(\frac{\partial^{2}w}{\partial x ^{2}}+\frac{\partial^{2}w}{\partial y^{2}}+\frac{\partial^{2}w}{\partial z^{2}}\right)\] \[\qquad-k(t,x,y,z)\]

The dataset used to train and evaluate the model's performance was generated
using a test solution, which is given by:

\[u :=(t,x,y,z)\mapsto\mathrm{e}^{-t}\sin(\pi x)\sin(\pi y)\sin(\pi z) \tag{4.18}\] \[v :=(t,x,y,z)\mapsto\mathrm{e}^{-t}\left(x^{2}-x\right)\left(y^{2}-y \right)\left(z^{2}-z\right)\] (4.19) \[w :=(t,x,y,z)\mapsto\mathrm{e}^{-t}\sin(\pi x)\sin(\pi y)\left(z^{2} -z\right)\] (4.20) \[p :=(t,x,y,z)\mapsto\mathrm{e}^{-t}xyz \tag{4.21}\]

## 5 Experimental Results and Analysis

In this section, we assess the performance of the improved models in solving the NSE in both 2D and 3D spaces.

### 2D Problem

#### 5.1.1 Dataset

In this study, the data used to train and validate the model is obtained from a computational fluid dynamics simulation for incompressible fluid flow around a circular cylinder. The simulation was performed using the Nektar++ [4] PDE solver, as described in the paper [22]. The simulation was carried out in a rectangular domain \(\Omega=[-15,25]\times[-8,8]\), with a constant velocity applied to the left boundary, zero pressure assumed at the right boundary, and periodic conditions applied to the top and bottom boundaries. The simulation is based on the Reynolds number, which is calculated as the ratio of the free stream velocity, cylinder diameter, and kinematic viscosity (\(R=\frac{U_{\infty}L}{\nu}\), assumed \(U_{\infty}=1\), \(L=1\), \(\nu=0.01\) ). There are 5000 points in the data set, with each point representing 200 time moments, and it includes values for the velocity components \(u\) and \(v\), as well as the pressure \(p\). The data is presented as a collection of \(\left\{t^{i},x^{i},y^{i},u^{i},v^{i},p^{i}\right\}\) values.

Figure 3: Approach to solve 3D NSE using the \(I_{2}\) model. \(n\) represents the number of neurons per layer and \(k\) represents the number of hidden layers.
#### 5.1.2 Implementations Details

We conducted several experiments2 to evaluate the model's performance. These experiments involved analyzing various factors that influence the model's ability by examining their loss functions, validation loss function, relative errors of predicted solutions, and the training duration required to learn and predict solutions. The objective of these experiments is to extract subsets of data from the global dataset.

Footnote 2: The implementations were executed using the Tesla T4 GPU.

We consider the following notation:

* \(I_{0}\): represents the initial implementation of PINN established in [22].
* \(I_{1}\): represents the implementation associated with the developed model based on assumption (11).
* \(I_{2}\): represents the implementation of the model without using any assumption.

Initially, we compared the training performance of one of our proposed models with the approach described in [22].

Table 1 presents a comparison between the initial model \(I_{0}\), introduced in paper [22], and one of the models suggested in this paper \(I_{2}\). The model \(I_{2}\) is trained using Pytorch, while the other model uses the TensorFlow \(1.x\). Both models are trained with the same GPU. For training, \(I_{0}\) randomly selected 5000 points from the global dataset. To compare with \(I_{2}\), we choose 8334 points from the global data and take 60% (\(\simeq 5000\)) of this data for training and 20% (\(\simeq 1667\)) for validation. As Table 1 shows, \(I_{2}\) outperforms \(I_{0}\) using just 2000 epochs, with lower training loss and shorter training time; \(I_{0}\) required 200,000 iterations.

Using a full batch for training model \(I_{0}\) can create challenges when trying to conduct additional experiments. While more data can be added to allow for longer training times, this also raises the risk of the optimization problem getting trapped at a saddle point or local minimum. To address this, we decided to implement the same model using mini-batch learning, validation loss, and early stopping criteria to alleviate these concerns. Additionally, we further improved the model, which we refer to as model \(I_{1}\), by including training on pressure.

Table 2 illustrates the common parameters between the two models \(I_{1}\) and \(I_{2}\). These parameters include the neural network's hyperparameters and the optimization parameters. It should be noted that both models were built from scratch using the

\begin{table}
\begin{tabular}{l c c} \hline  & \(I_{0}\) & \(I_{2}\) \\ \hline Software library & TF1.x & Pytorch \\ \hline Dataset & 1000000 & 8334 \\ \hline Training data & 5000 & 5000 \\ \hline Validation data & - & 1667 \\ \hline Iterations-Epochs & 200000 & **2000** \\ \hline Batch learning & full-batch & **mini-batch** \\ \hline Batch size & 5000 & 64 \\ \hline Training loss & 2.825\(\times 10^{-1}\) & **1.2\(\times 10^{-2}\)** \\ \hline Validation loss & - & **6.03e-5** \\ \hline Training duration & \(>\)8h & **1h40min** \\ \hline GPU & Yes & Yes \\ \hline \end{tabular}
\end{table}
Table 1: The initial training comparison of \(I_{0}\) and \(I_{2}\).
PyTorch software library.

In order to evaluate the precision of the models, we calculate the Relative Error (RE) between the predicted solution and the actual solution. RE is calculated using the following formula:

\[RE_{u}=\frac{||predicted\ value(u)-true\ value(u)||}{||true\ value(u)||} \tag{5.1}\]

To conduct additional experiments with larger datasets, we increased the batch size for both implementations, as shown in Table 4, to prevent excessive training times. We also have further experiments planned in Tables 5 and 6 to evaluate the impact of major hyperparameters, such as the number of neurons and layers, on the models. For training and evaluating the models, we used a batch size of 256 and an

\begin{table}
\begin{tabular}{l c} \hline
**Parameters** & **Details** \\ \hline Global data size & 5000*200 \\ \hline Hidden layer & 8 \\ \hline Neurons per layer & 20 \\ \hline Activation function & Tanh \\ \hline Optimizer & Adam \\ \hline Initial learning rate & 0.001 \\ \hline Early stopping criteria & 100 \\ \hline Training data size & 60\% \\ \hline Validation data size & 20\% \\ \hline Test data size & 20\% \\ \hline GPU & Yes \\ \hline \end{tabular}
\end{table}
Table 2: Parameters shared in both models \(I_{1}\) and \(I_{2}\).

\begin{table}
\begin{tabular}{l|c|c|c|c|c|c} \hline Experiment & \multicolumn{2}{c|}{\(E_{1}\)} & \multicolumn{2}{c|}{\(E_{2}\)} & \multicolumn{2}{c}{\(E_{3}\)} \\ \hline Model & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) \\ \hline Data size & \multicolumn{2}{c|}{9000} & \multicolumn{2}{c|}{10000} & \multicolumn{2}{c}{20000} \\ \hline Early stopping criteria & \multicolumn{2}{c|}{} & \multicolumn{2}{c}{100} & \multicolumn{2}{c}{} \\ \hline Batch size & \multicolumn{2}{c|}{} & \multicolumn{2}{c}{256} \\ \hline Training Loss & 1.62e-3 & 1.65e-3 & 1.06e-3 & 1.53e-3 & 6.98e-4 & 1.38e-3 \\ \hline Validation loss & 1.66e-3 & 1.86e-3 & 1.33e-3 & 1.78e-3 & 6.81e-4 & 1.18e-3 \\ \hline \(RE_{u}\) & 2.19e-2 & 2.26e-2 & 2.03e-2 & 2.44e-2 & 1.44e-2 & 1.84e-2 \\ \hline \(RE_{v}\) & 7.07e-2 & 6.14e-2 & 5.59e-2 & 6.45e-2 & 3.99e-2 & 5.12e-2 \\ \hline \(RE_{p}\) & 7.30e-2 & 8.15e-2 & 7.30e-2 & 8.26e-2 & 4.81e-2 & 5.87e-2 \\ \hline Training duration & 1h15min & 28min & 1h45min & 32min & 2h16min & 46min \\ \hline \end{tabular}
\end{table}
Table 3: Numerical results of experiments \(E_{1},E_{2}\) and \(E_{3}\) for solving the 2D NSE using models \(I_{1}\)and \(I_{2}\). Both models are executed with the same early stopping criteria and batch size.
early stopping criterion of 100 epochs as patience, along with the selection of 20,000 collocation points. Based on the numerical results presented in Tables 1, 3, 4 5, and 6, it is evident that the suggested models \(I_{1}\) and \(I_{2}\) outperform other models in solving the 2D NSE. While the \(I_{1}\) model has been considerably enhanced by integrating deep learning techniques, the training process still requires a considerable amount of time. In contrast, the \(I_{2}\) model demonstrates its capability to compete with the \(I_{1}\) model and achieve good results, while requiring a training time that is at most half that of the \(I_{1}\). Furthermore, the results indicate that modifying the architecture and hyperparameters can have a discernible impact on the final results. In fact, as the number of neurons and hidden layers increase, the predicted solution associated with model \(I_{2}\) achieves a higher level of convergence.

\begin{table}
\begin{tabular}{l|c|c|c|c|c|c} \hline Experiment & \multicolumn{2}{c|}{\(E_{7}\)} & \multicolumn{2}{c|}{\(E_{8}\)} & \multicolumn{2}{c}{\(E_{9}\)} \\ \hline Model & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) \\ \hline Hidden layers & \multicolumn{5}{c}{10} \\ \hline Neurons per & \multicolumn{5}{c|}{10} \\ layer & \multicolumn{5}{c|}{10} \\ \hline Training loss & 3.05e-3 & 5.21e-3 & 9.60e-4 & 9.04e-4 & 2.96e-4 & 2.67e-4 \\ \hline Validation loss & 2.93e-3 & 5.00e-3 & 8.95e-4 & 9.00e-4 & 3.19e-4 & 2.62e-4 \\ \hline \(RE_{u}\) & 2.82e-2 & 4.36e-2 & 1.65e-2 & 1.56e-2 & 8.91e-3 & 7.11e-3 \\ \hline \(RE_{v}\) & 9.72e-2 & 1.14e-1 & 4.54e-2 & 4.37e-2 & 2.88e-2 & 2.34e-2 \\ \hline \(RE_{p}\) & 1.19e-1 & 1.41e-1 & 5.60e-2 & 4.84e-2 & 3.08e-2 & 2.78e-2 \\ \hline Training duration & 4h40min & 1h58min & 4h40min & 1h58min & 2h28min & 1h20min \\ \hline \end{tabular}
\end{table}
Table 4: Numerical outcomes obtained from experiments \(E_{4},E_{5}\), and \(E_{6}\) after augmenting the data size and elevating the batch size.

\begin{table}
\begin{tabular}{l|c|c|c|c|c|c} \hline Experiment & \multicolumn{2}{c|}{\(E_{4}\)} & \multicolumn{2}{c|}{\(E_{5}\)} & \multicolumn{2}{c}{\(E_{6}\)} \\ \hline Model & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) \\ \hline Data size & \multicolumn{5}{c|}{30000} & \multicolumn{5}{c|}{40000} & \multicolumn{5}{c}{60000} \\ \hline Early stopping criteria & \multicolumn{5}{c}{100} \\ \hline Batch size & \multicolumn{5}{c}{500} \\ \hline Training Loss & 7.26e-4 & 7.11e-4 & 4.02e-4 & 7.33e-4 & 6.62e-4 & 7.00e-4 \\ \hline Validation loss & 6.85e-4 & 6.82e-4 & 3.69e-4 & 6.41e-4 & 5.44e-4 & 6.09e-4 \\ \hline \(RE_{u}\) & 1.37e-2 & 1.30e-2 & 9.20e-3 & 1.30e-2 & 1.19e-2 & 1.30e-2 \\ \hline \(RE_{v}\) & 3.93e-2 & 3.82e-2 & 2.85e-2 & 3.61e-2 & 3.66e-2 & 3.48e-2 \\ \hline \(RE_{p}\) & 4.53e-2 & 4.30e-2 & 4.29e-2 & 3.99e-2 & 4.23e-2 & 4.10e-2 \\ \hline Training duration & 2h30min & 1h12min & 3h30min & 1h40min & 2h55min & 1h30min \\ \hline \end{tabular}
\end{table}
Table 5: Experiments with fixing the number of hidden layers and changing the number of neurons per layer.
In the next section, we try to solve the non-stationary NSE in a 3D domain using model \(I_{2}\). We selected this approach due to its perceived efficacy in effectively addressing the complexity of the problem at hand.

### 3D Problem

#### 5.2.1 Dataset

In order to construct a model for the 3D problem, a dataset comprising variables \(\left\{x,y,z,t,u,v,w,p\right\}\) is required. In this study, a total of 7318 collocation points were generated in the spatial domain using a popular and freely available mesh generator known as Gmsh [8]. The time interval was discretized into 200 moments. The solution associated with the points was calculated using the equations (4.18)-(4.21).

Afterward, the resulting dataset was rearranged into the form \(\left\{x^{i},y^{i},z^{i},t^{i},u^{i},v^{i},w^{i},p^{i}\right\}\) for all \(1\leq i\leq(7318\times 200)\), which allowed for more in-depth analysis and interpretation of the model's results.

#### 5.2.2 Numerical Results

Due to the complexity of the problem, a substantial dataset was utilized for this experiment3.

Footnote 3: The Tesla A40 GPU was employed to run the implementations.

Figure 4 displays the training and validation losses with respect to the data size. Based on the conducted experiments, we have successfully run the model using 35% of the entire dataset, which comprises over one million examples. Notably, the model's training can be completed in just 3 hours and 30 minutes using 500,000 examples, demonstrating its excellent efficiency. All experiments were conducted using the same parameters listed in Table 2 except the global data size in this case is \(7318\times 200\).

The results presented in Table 7 were obtained using the number of epochs as a stopping criterion, with 2000 epochs chosen as the optimal value for training the model. This number was considered large enough to provide sufficient time for the model to optimize and investigate its performance capabilities. By examining the results shown in Figure 5, we can conclude that the relative errors exhibit a decreasing trend as the data size increases. This observation aligns with the expected behavior of relative errors with increasing data size, as larger datasets enable more precise calculations by providing a greater amount of information.

\begin{table}
\begin{tabular}{l|c|c|c|c|c|c} \hline Experiment & \multicolumn{2}{c|}{\(E_{10}\)} & \multicolumn{2}{c|}{\(E_{11}\)} & \multicolumn{2}{c}{\(E_{12}\)} \\ \hline Model & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) & \(I_{1}\) & \(I_{2}\) \\ \hline Hidden layers & \multicolumn{2}{c|}{12} & \multicolumn{2}{c|}{15} & \multicolumn{2}{c}{18} \\ \hline Neurons per & \multicolumn{2}{c|}{20} \\ layer & \multicolumn{2}{c|}{} & \multicolumn{2}{c|}{} & \multicolumn{2}{c|}{} \\ \hline Training loss & 1.12e-3 & 1.49e-3 & 1.16e-3 & 8.86e-4 & 1.27e-3 & 9.40e-4 \\ \hline Validation loss & 1.07e-3 & 1.38e-3 & 1.10e-3 & 8.50e-4 & 1.11e-3 & 8.69e-4 \\ \hline \(RE_{u}\) & 1.75e-2 & 2.05e-2 & 1.86e-2 & 1.54e-2 & 1.76e-2 & 1.68e-2 \\ \hline \(RE_{v}\) & 5.11e-2 & 5.42e-2 & 5.21e-2 & 4.25e-2 & 5.59e-2 & 4.23e-2 \\ \hline \(RE_{p}\) & 6.47e-2 & 5.85e-2 & 6.18e-2 & 5.08e-2 & 5.90e-2 & 4.96e-2 \\ \hline Training & \multicolumn{2}{c|}{} & \multicolumn{2}{c|}{} & \multicolumn{2}{c|}{} & \multicolumn{2}{c|}{} \\ duration & \multicolumn{2}{c|}{2h20min} & \multicolumn{2}{c|}{1h} & \multicolumn{2}{c|}{3h15min} & \multicolumn{2}{c|}{1h43min} & \multicolumn{2}{c|}{4h10min} & \multicolumn{2}{c|}{2h} \\ \hline \end{tabular}
\end{table}
Table 6: Outcomes of fixing the number of neurons per layer and changing the number of layers.
## 6 Conclusions and Perspectives

The NSE represent a highly intricate challenge to solve, given their dynamic and nonlinear nature as PDEs. Recent research has shown that physics-informed neural networks (PINN) offer a promising solution to this challenge. Our study has demonstrated the effectiveness of improved models based on PINN for solving the NSE in 2D and 3D cases. We have shown that these models offer significant advantages over the original PINN approach, including enhanced computational efficiency, more stable convergence, and reduced relative errors. The approach employed in this paper involves utilizing a test solution in the three-dimensional problem for training and evaluating the efficacy of models aimed at solving complex PDEs. This methodology serves as a viable alternative to conducting real experiments or employing PDE solvers, which may prove more costly. These findings suggest that further refinement of this method could open up new opportunities for its application to real-world data from experiments in three dimensions. Overall,

\begin{table}
\begin{tabular}{l|c|c|c|c|c|c} \hline Experiment & \(E_{1}\) & \(E_{2}\) & \(E_{3}\) & \(E_{4}\) & \(E_{5}\) & \(E_{6}\) \\ \hline Data size & 40000 & 100000 & 200000 & 300000 & 400000 & 500000 \\ \hline Loss function & 2.13e-5 & 7.23e-6 & 3.46e-6 & 2.56e-5 & 1.49e-5 & 5.37e-6 \\ \hline Validation loss & 1.99e-5 & 5.84e-6 & 3.23e-6 & 2.75e-6 & 1.98e-6 & 2.31e-6 \\ \hline \(RE_{u}\) & 2.13e-2 & 1.20e-2 & 1.06e-2 & 1.13e-2 & 6.97e-3 & 5.69e-3 \\ \hline \(RE_{v}\) & 3.89e-1 & 2.50e-1 & 4.56e-1 & 1.94e-1 & 9.65e-2 & 8.19e-2 \\ \hline \(RE_{w}\) & 6.23e-2 & 2.25e-2 & 2.42e-2 & 2.48e-2 & 2.37e-2 & 1.35e-2 \\ \hline \(RE_{p}\) & 1.58e-2 & 8.47e-3 & 4.24e-3 & 4.98e-3 & 3.99e-3 & 2.91e-3 \\ \hline Training duration & 1h5min & 2h36min & 5h27min & 8h & 11h & 13h41min \\ \hline \end{tabular}
\end{table}
Table 7: The outcomes obtained through the implementation of the \(I_{2}\) model for solving the 3D NSE, wherein the number of epochs was employed as a stopping criterion.

Figure 4: Effects of varying dataset size on training and validation loss.
our work highlights the potential of PINN as a powerful tool for tackling complex physical modeling problems and paves the way for exciting future developments in this area.

In addition to the promising results presented in this study, there are several exciting perspectives for future research. One possibility is to explore the potential of PINN in solving more complex problems, such as those involving multiphase flows or turbulence. Another avenue for investigation is to extend the methodology developed in this paper to other types of PDEs, including those with non-linear boundary conditions or time-dependent coefficients. Moreover, the use of PINN could be extended to other areas of physics and engineering, such as solid mechanics, electromagnetics, and quantum mechanics, to name a few.

Another perspective to consider is the potential for combining PINN with other machine learning techniques, such as deep learning or reinforcement learning, to further enhance their capabilities. The integration of PINN with other methods may allow for more efficient and accurate solutions to complex physical problems, while also providing insights into the underlying physical phenomena.

Finally, it is essential to recognize the potential impact of solving PDEs in various fields of science and engineering. The ability to accurately model and predict complex physical phenomena have significant implications for the design of new materials, optimization of manufacturing processes, and the development of new technologies. As such, the continued development and refinement of PDE solvers have the potential to transform how we approach modeling and simulation in various domains, leading to exciting discoveries and advancements.

Figure 5: Relative errors of velocity and pressure components in various experiments of 3D NSE solved using the \(I_{2}\) model.