Discontinuity Computing Using Physics-Informed Neural Networks

Abstract

Simulating discontinuities has been a long-standing challenge, especially when dealing with shock waves characterized by strong nonlinear features. Despite their promise, the recently developed physics-informed neural networks (PINNs) have not yet fully demonstrated their effectiveness in handling discontinuities when compared to traditional shock-capturing methods. In this study, we reveal a paradoxical phenomenon during the training of PINNs when computing problems with strong nonlinear discontinuities. To address this issue and enhance the PINNs’ ability to capture shocks, we propose PINNs-WE (Physics-Informed Neural Networks with Equation Weight) method by introducing three novel strategies. Firstly, we attenuate the neural network’s expression locally at ‘transition points’ within the shock waves by introducing a physics-dependent weight into the governing equations. Consequently, the neural network will concentrate on training the smoother parts of the solutions. As a result, due to the compressible property, sharp discontinuities emerge, with transition points being compressed into well-trained smooth regions akin to passive particles. Secondly, we also introduce the Rankine–Hugoniot (RH) relation, which is equivalent to the weak form of the conservation laws near the discontinuity, in order to improve the shock-capturing preformance. Lastly, we construct a global physical conservation constraint to enhance the conservation properties of PINNs which is key to resolve the right position of the discontinuity. To illustrate the impact of our novel approach, we investigate the behavior of the one-dimensional Burgers’ equation, as well as the one- and two-dimensional Euler equations. In our numerical experiments, we compare our proposed PINNs-WE method with a traditional high-order weighted essential non-oscillatory (WENO) approach. The results of our study highlight the significant enhancement in discontinuity computing by the PINNs-WE method when compared to traditional PINNs.

Appendices

Appendix A: Omissible Boundary Condition in PINNs

In traditional numerical methods used for solving initial boundary value problems, such as finite element and finite difference methods, boundary conditions play essential roles. They serve two primary purposes: mathematically defining the problem and numerically closing the discretization near the boundaries of the schemes. However, in PINNs, which do not rely on discretization or logical relations between sampling points, the need for boundary conditions to serve the latter purpose is eliminated. This simplifies the process of setting boundary conditions significantly.

In our experimentation with PINNs, we have observed that two types of boundary conditions may can be omitted when solving initial boundary value problems. The first type is constant boundary conditions, where the boundary values remain unchanged over time after being determined by the initial conditions. The second type is outflow boundary conditions, where information purely flows out of the domain.

Traditionally, setting outflow boundary conditions can be challenging because they require closed and discrete boundary conditions to exist while also ensuring that the set boundary conditions do not influence the internal flow. In PINNs, these boundary conditions can be omitted, yet they still ensure complete outflow characteristics at the boundary.

1.1 Appendix A.1: 1D Linear Transport Equation Problem

To illustrate the influence of omitting these two types of boundary conditions, we consider two problems solved with PINNs. The governing equation for both problems is chosen as:

$$\begin{aligned} \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0, \quad -1< x < 2, \end{aligned}$$

with the first initial condition as

$$\begin{aligned} u(x,0)= \left\{ \begin{aligned}&\textrm{sin} (2\pi x) + 1 \quad \textrm{if}&0\le x \le 1, \\&1&\textrm{else}. \end{aligned}\right. \end{aligned}$$

The computational time is set as $0 \le t \le 0.5$. Both left and right boundary conditions are kept constant throughout the computation time. We compare the results with and without the inclusion of boundary conditions in Table 2.

Our results indicate that when boundary conditions are constant, there is minimal influence when omitting them from the loss function. The number of residual points sampled within the computational domain is denoted as $N_f = 10000$, while $N_{\textrm{IC}} = 100$ represents the number of residual points sampled on the initial condition. Additionally, $N_{\textrm{BC}} = 100$ is the number of residual points sampled on the left and right boundary conditions. All these residual points are uniformly distributed.

Table 2 Accuracy comparison with and without constant boundary constraint

For the second initial condition, it is defined as:

$$\begin{aligned} u(x,0)= \textrm{sin} (2\pi x) + 1, \end{aligned}$$

and the computational time is extended to $0 \le t \le 5$. The left boundary condition is set as an inflow boundary condition:

$$\begin{aligned} u(-1,t) = \textrm{sin}(2\pi (-t -1)) + 1. \end{aligned}$$

On the other hand, the right boundary condition is chosen as an outflow condition. In traditional methods with outflow boundaries, there are various approaches to closing the discretization. One common method is to extrapolate the value at $u(2 + \Delta x, t)$ using u(2, t). Here, $\Delta x$ represents the mesh size and can be calculated as $\Delta x = L_x/N_x = 0.03$, where $L_x$ denotes the length in space, and $N_x$ is the number of points in the x direction. Alternatively, more complex zero-gradient outflow conditions, such as setting $u_x(2+\Delta x, t) = u_x(2, t)$, can be employed to minimize the impact of boundary conditions. In this context, as you mentioned, using uniform points simplifies the choice of $\Delta x$.

Table 3 Accuracy comparison with different set of outflow boundary constraint

In our analysis, we have tested three cases: extrapolated boundary conditions, zero-gradient boundary conditions, and omitting boundary conditions altogether with PINNs. To minimize the impact of network randomness, the results reported in Table 3 represent the averages of ten separate runs with different random seeds.

The results demonstrate that even when using zero-gradient boundary conditions, the boundary conditions still influence the accuracy of internal flows. Conversely, omitting the outflow boundary condition appears to be a better fit for representing the true physical behavior.

1.2 Appendix A.2: 2D Vortex Evaluation Problem

Then we test a 2D vortex problem governed by 2D Euler equations (16). The initial condition is considered as

$$\begin{aligned} \left( \begin{array}{c} \rho \\ u \\ v \\ p \end{array}\right) =\left( \begin{array}{c} (1+\delta T)^{1 /(\gamma -1)} \\ 1+\left( l_{y} / 2-y\right) \frac{\sigma }{2 \pi } e^{0.5\left( 1-r^{2}\right) } \\ 1+\left( x-l_{x} / 2\right) \frac{\sigma }{2 \pi } e^{0.5\left( 1-r^{2}\right) } \\ (1+\delta T)^{\gamma /(\gamma -1)} \end{array}\right) . \end{aligned}$$

(A.1)

Here $\delta T$ is the perturbation in the temperature and is given by

$$\begin{aligned} \delta T=-\frac{(\gamma -1) \sigma ^{2}}{8 \gamma \pi ^{2}} e^{\left( 1-r^{2}\right) }, \end{aligned}$$

Fig. 21

The result and point-wise error of the 2D vortex evolution problem at $t=1$

where $r^{2}=\left( x-x_c\right) ^{2}+\left( y-y_c / 2\right) ^{2}$ and the vortex strength $\sigma =5$, $x_c = y_c = 2$. The computational domain is given as $[0,5] \times [0,5]$. Here $\gamma$ is 1.4. The initial conditions lead advection of a non-linear vortex at an angle of $45^{\circ }$ with the x-axis and the numerical solutions are obtained after $t=1$. The left and bottom boundaries are constant inflow while right and top boundaries are outflows. The $N_{f}=400000$ and $N_{\textrm{IBCs}}= 10000$. Result at $t=1$ is compared with the exact solution in Fig. 21. And convergence of $L_2$ and $L_\infty$ relative errors are presented in Table 4. We show that omitting the boundary conditions in this case can obviously improve the accuracy.

In summary, the PINNs method offers the significant advantage of simplifying the handling of boundary conditions. By not relying on a discrete grid and having the capability to capture temporal behavior, it becomes feasible to omit constant boundary conditions and simplify the setting of outflow boundary conditions. This flexibility and ease of handling boundary conditions make the PINNs method a practical choice for solving hyperbolic equations.

Appendix B: Classical PINNs in Computing Linear and Weak Discontinuities

The effectiveness of Physics-Informed Neural Networks (PINNs) in solving problems with smooth solutions has been extensively studied and demonstrated. However, research on solving problems involving discontinuities is still relatively limited. In this section, we aim to first validate the performance of the classical PINNs approach in solving linear discontinuities and weak discontinuities. Here, we consider several test cases, include linear transport equation problems involving linear discontinuities to verify the solution with linear discontinuities, and one-dimensional Riemann problems with rarefaction waves to verify the solution with derivatives discontinuities.

Table 4 Accuracy comparison with and without constant boundary constraint, 2D vortex evolution problem

1.1 Appendix B.1: 1D Linear Transport Equation with Complex Waveforms

The equation is given as

$$\begin{aligned} \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0, \quad -1< x < 2, \end{aligned}$$

with the initial condition

$$\begin{aligned} u(x,0)=\left\{ \begin{array}{ll} \frac{1}{6}(G(x, \beta , z-\delta )+G(x, \beta , z+\delta )+4 G(x, \beta , z)), &{} -0.8 \le x \le -0.6 \\ 1, &{} -0.4 \le x \le -0.2 \\ 1-|10(x-0.1)|, &{} 0 \le x \le 0.2 \\ \frac{1}{6}(F(x, \alpha , \alpha -\delta )+F(x, \alpha , \alpha +\delta )+4 F(x, \alpha , a)), &{} 0.4 \le x \le 0.6 \\ 0, &{} \text{ otherwise } \end{array}\right. \end{aligned}$$

As in Ref. [19], the constants have been assigned specific values: $a = 0.5$, $z = -0.7$, $\delta = 0.005$, $\alpha = 10$, and $\beta = \log {2}/36\delta ^2$. The solution exhibits a diverse set of discontinuities, including a smooth combination of Gaussians, a square wave, a sharp triangle wave, and a half ellipse. The equation possesses an exact solution of the form $u(x,t) = u(x-t,0)$.

Here we use the classical PINNs with the loss function of the mean square error from two part

$$\begin{aligned} \textrm{Loss} = \textrm{MSE}_f + 10 \textrm{MSE}_{\textrm{IC}} \end{aligned}$$

and ignoring the boundary conditions as we talked in Appendix A.

The function residual points is taken as $N_f = 10000$ and the initial residual points is $N_{\textrm{IC}}=1000$, and all they are sampled uniformly.

We use the ADAM optimizer with a learning rate of 0.01 then follows L-BFGS optimizer with a learning rate of 1 until the loss converges.

The solution and the error solved by PINNs are present in Fig.. Then we test the trained network with 100 uniform sampling points at time $t=0.5$ and compare it we the exact solution. And the L2 relative error is 0.005 ( the average of ten cases with different random seeds to eliminate the influence randomness).

Fig. 22

The result and point-wise absolute error of the 1D linear transport equation with complex waveforms

Fig. 23

Test result at $t=0.5$ of the 1D linear transport equation with complex waveforms

1.2 Appendix B.2: 2D Linear Transport Equation with Initial Interface Evolution

Interface tracking is a widely encountered scientific and engineering problem that often requires solving the linear transport equation with a given velocity field:

$$\begin{aligned} \frac{\partial u}{\partial t} + a_x(t,x,y) \frac{\partial u}{\partial x} + a_y(t,x,y) \frac{\partial u}{\partial y} = 0. \end{aligned}$$

We consider a complex case introduced by [28]. The velocity field is given as

$$\begin{aligned} \left\{ \begin{aligned}&a_x = \sin ^2(\pi x)\sin (2\pi y), \\&a_y = -\sin ^2(\pi y)\sin (2\pi x), \\ \end{aligned} \right. \end{aligned}$$

when $t<T/2$. Then we take an opposite velocity field to rotate it back. and the initial conditions are

$$\begin{aligned} u_0(x,y) = \left\{ \begin{aligned}&0,&\sqrt{(x-x_r)^2+(y-y_r)^2}> 0.15,\\&1,&\textrm{else}. \end{aligned} \right. \end{aligned}$$

The computation domain is $0 \le x,y \le 1$.

The numbers of residual points sampled in the computational domain and on the initial condition are $N_f = 10000$ and $N_\textrm{IC}=1000$, respectively. They are sampled with the Latin hypercube sampling (LHS) method.

Fig. 24

The result of the 2D linear transport equation with vortex stretching, $t_1 = 0.5$, $t_2 = 1$ and $t_3 = 2$, black line is the initial position of the interface

We have test three networks with different t. After the training, we take two test sets with 10000 uniform points at $T=t/2$ and $T=t/2$ in each case, respectively. The result are shown in Fig. 24. We can see PINNs capture the interface sharply with little dissipation as there is no necessary dissipation introduced to capture the discontinuity. So the possible dissipation is the error from the approximation that can be controlled by the convergence of the loss function.

Fig. 25

The result of the double rarefaction problem

1.3 Appendix B.3: 1D Riemann Problem with Double Rarefaction Waves

After we test two linear discontinuity cases, then we will test the performance of PINNs in solving weak discontinuities. Rarefaction wave is typically one kind of weak discontinuities that they have only derivative discontinuity and $C^0$ smoothness.

We choose a classical double rarefaction problem described by Toro [29], Chapter 6. In this test, the center of the domain is evacuated as two rarefaction waves propagate in each direction, outward from the center.

The initial conditions are: The governing equation is the 1D Euler equation that is given in Sect. 2 with the initial condition as

The initial condition is given as

$$\begin{aligned} (\rho ,u,p) = \left\{ \begin{aligned}&(1,-2,0.4), \quad&\textrm{if} \quad 0\le x \le 0.5,\\&(1,2,0.4),&\quad \textrm{if} \quad 0.5< x \le 1.\\ \end{aligned}\right. \end{aligned}$$

(B.1)

And the computational time is $t = 0.1$. This problem is hard to solve as their are vacuum created at the center after the expansion. Especially the internal energy $e=\frac{p}{/}{\rho (\gamma -1)}$ may have large error inside the low pressure/density region.

We use a two step training strategy, the initial residual points $N_{IC} = 100$ and we pre-train the network with $N_{f} = 2000$ first to convergence and then refine it with more residual points as $N_f = 10000$. The results are shown in Fig. 25. It shows that PINNs performs good in solving rarefaction waves even with vacuum regions.

Then we give a conclusion of appendix B. We test three classical interface tracking cases with the classical PINNs without any optimization. The results show that PINNs is very good at solving linear discontinuous problem, the discontinuous keeps sharp without obviously dissipation. So PINNs have much potential to be powerful methods for solving multi-materials problems.