Discontinuity Computing Using Physics-Informed Neural Network

Li Liu1 ,Shengping Liu1 , Hui Xie1 , Fansheng Xiong1 ,Tengchao Yu1 , Mengjuan Xiao1 , Lufeng Liu1 , Heng Yong1,∗
1 Institute of Applied Physics and Computational Mathematics, Beijing 100094, China

Abstract

Simulating discontinuities is a long standing problem especially for shock waves with strong nonlinear feather. Despite
being a promising method, the recently developed physics-informed neural network (PINN) is still weak for calculating
discontinuities compared with traditional shock-capturing methods. In this paper, we intend to improve the shock-
capturing ability of the PINN. The primary strategy of this work is to weaken the expression of the network near
arXiv:2206.03864v2 [[Link]-dyn] 6 Aug 2022

discontinuities by adding a gradient-weight into the governing equations locally at each residual point. This strategy
allows the network to focus on training smooth parts of the solutions. Then, automatically affected by the compressible
property near shock waves, a sharp discontinuity appears with wrong inside shock transition-points compressed into well-
trained smooth regions as passive particles. We study the solutions of one-dimensional Burgers equation and one- and
two-dimensional Euler equations. Compared with the traditional high-order WENO-Z method in numerical examples,
the proposed method can substantially improve discontinuity computing.
Keywords: PINN, Shock capturing, Compressible flow, Euler equations, Discontinuity calculation

1. Introduction

Calculating shock waves and other discontinuities sharply and without oscillations is essential for solving hyperbolic
equations. The first study on shock-capturing dates back to the development by von Neumann and Richtmeyer [1],
who introduced artificial viscosity into a staggered Lagrangian scheme to solve compressible flows. Nowadays, various
advanced and high-order methods allow to simulate problems with shock waves. These methods include essential non-
oscillatory (ENO) [2], weighted ENO (WENO) [3], and discontinuous Galerkin [4] methods. More information about the
development of shock-capturing methods can be found in [5, 6].
Owing to their rapid development, machine learning and neural networks (NNs) have been used to solve partial
differential equations [7, 8, 9, 10, 11]. In particular, the physics-informed neural network (PINN) attracts much research
attention. PINN encodes partial differential equations (PDEs) or other model equations as one of its components. Given
its generality, PINN has been used for solving equations in many fields [12]. For hyperbolic equations, Patel et al. [13]
proposed a PINN that can discover thermodynamically consistent equations ensuring hyperbolicity for inverse problems
in shock hydrodynamics. Mao et al. [14] studied one- (1D) and two-dimensional (2D) Euler equations with shock waves
and used clustered training samples around a high gradient area to improve the solution accuracy in that area while
preventing error propagation to the entire domain. Jagtap et al. [15] proposed the conservative PINN that splits the
computing domain into several small subdomains using different NNs to solve Burgers and Euler equations. Jagtap et al.
[16] studied inverse problems in supersonic flows. The above mentioned studies have shown the effectiveness of PINN to
handle inverse problems with prior information about the development of the flow structures, such as density gradients.
However, to study forward problems, the applicability of the original PINN has been limited to simple problems such as
tracking moving shock waves. Based on traditional methods, Patel et al. [13] constructed a mesh-based control-volume
PINN and introduced entropy and total-variation-diminishing conditions into the network. Papados [17] simulated the
shock-tube problem with the computing domain extended, obtaining outstanding results without introducing non-physics
viscosity terms into the equations.

1
 We aim to enhance the shock-capturing ability of a PINN, especially for complex problems with shock generation.
However, the shock wave has zero thickness in compressible inviscid flow theoretically. Thus, it cannot be governed by
the strong form of the differential equations given its infinite gradient but instead controlled by a physical process from
the left and right regions (Rankine-Hugoniot conditions). Besides, there are also no theory that can guarantee NNs
approximate any first-order discontinuous functions. Hence, if residual points wrongly fall inside a shock region, large
equation losses exist in those points given their large gradients. These points are called transition-points [18] in this
paper. A NN may focus on handling transition-points as they carry most of the loss. However, a NN cannot increase
the gradient to decrease the thickness of the shock owing to the reason talked above. On the other hand, the physics
process compresses the region over time. As a result, transition-points fall into a paradoxical status, in which increasing
the gradient seriously increases the total loss because the points are not governed by the equations. However, decreasing
the gradient also increases the total loss as it conflicts with physical compression and away from the real solution. More
seriously, as the total loss is the sum or average of every point, the transition points attract training and influence the
convergence of other points in smooth regions. Here, we can compare the above process with a traditional high-order
method, such as the finite-volume WENO method, to illustrate the difficulty of PINN. In a WENO method, when a
cell is inside a discontinuous region, the order of the scheme automatically reduces to no more than the second-order.
Thus, a large dissipation is appended into the cell to obtain a non-oscillatory result. However, the given cell with large
dissipation and consequently large error does not influence the accuracy order of other cells beyond the discontinuity.
We introduce a ‘retreat to advance’ strategy into PINN. To break the paradoxical status of transition-points and
obtain a sharp discontinuity, PINN avoids training the shock waves and focuses on training other smooth regions by
weakening the network expression in strong compression regions. Then, the compressing property near shock waves
allows a sharp and exact shock to appear. This is done by multiplying a local positive and compression-related weight
to the governing equations to adjust the expressions of the NN in different residual points.
The remainder of this paper is organized as follows. In Section 2, we detail the weighted-equation (WE) method.
Then, various 1D and 2D forward examples are studied to show the effectiveness of the proposed method. Finally,
conclusions are drawn in Section 4.

2. Method

2.1. PINN for conservative hyperbolic PDE

We consider the following conservative hyperbolic PDE

∂U(x)
+ ∇ · F(U) = 0, x = (t, x1 , x2 , · · · ) ∈ Ω, (1)
∂t

with the initial and boundary conditions (IBs)

IB(U, x) = 0 on ∂Ω, (2)

and we treat the initial condition in the same way with Dirichlet boundary condition.
PINN mainly consists of two parts. The first part is a NN Û(x; θ) to approximate the relation of U(x) with trainable
parameters θ. The second part is informed with the governing equations and the IB conditions to train the NN. The
calculation of ∂/∂t and ∇· in the PDE is performed by automatic derivative evaluation. More details about the PINN
for concection PDE are available in[14, 17].
The loss function used to train Û(x; θ) contains at least two parts to define the problem. One is controlled by the
equations, and the other one is given by the IBs of the problem,

L = LPDE + ε1 LIBs . (3)

2
 Figure 1: Architecture of PINN-WE for solving conservative hyperbolic equations.

To define the loss, we choose a set of residual points inside the domain Ω and another set of points in ∂Ω as SPDE and
SIBs , respectively. Then
1 X 1 X
L= G2i + ε1 (Ûe,i − Ue0,i )2 , (4)
|SPDE | |SIBs |
xi ∈SPDE xi ∈SIBs

where Gi := ∂t Û(xi ) + ∇ · F(Û(xi )) and Gi = 0 is the governing equations at residual point xi ∈ SPDE and Ue0,i
represents the given initial IB at residual point xi ∈ SIBs .
ε1 is the weight to adjust the confinement strength of IBs[16, 19, 17]. A common condition is giving more weight to
the points on boundaries.
In each term of a loss function, averaging the residual across all residual points is common to obtain the total training
loss. Averaging may be suitable and convenient for problems with smooth solutions. However, when a discontinuity occurs,
the gradient becomes theoretically infinity and cannot be described directly by differential equations. Consequently,
transition-points defined inside the discontinuity may introduce large errors and function loss. We will show this in the
following case.

2.2. Analysis of transition-points based on Burgers equation

Here we consider inviscid Burgers equation which is a scalar nonlinear hyperbolic equation to illustrate the training
problem of transition-points.
The equation and IBs are given as

∂u ∂(u2 /2)
+ = 0, x ∈ [0, 2], t ∈ [0, 1],
∂t ∂x
u(0, x) = −sin(π(x − 1)), (5)

u(t, 0) = u(t, 2) = 0.

We first solve this problem with a traditional PINN introduced in above. Fig.2 gives the loss history and u and the
residual as functions of x at t = 1 at different training epochs (1000,3000,5000,8000 and 11500). It is clearly that PINN
first tends to reach a global smooth solution to fit the IBs to get a small LIBs (Slice 1). And then the training attends
to further deduce the residual to reach the exact solution of the problem. However, when the transition-points reach a
relatively high gradient (Slice 2), they nearly carry all the function loss LPDE (Fig,2.c). Then the training will fall into
a paradoxical status: no matter increase (Slice 3) or decrease the gradient (Slice 4) of transition-points will increase the
function loss. This is due to those points can not be controlled by the equation directly, then decrease the gradient will
remove from the exact solution but increase the gradient also will increase the residual. As shown in the loss history,
after 3000 epochs, the training struggles in those transition-points and can not decrease the total loss effectively. More
seriously, as a global method, the total loss also decide the error size of other regions.

3
 Loss_PDE
101 Loss_IBs 0.6
Loss_Total
Slice1 Slice2 Slice3 Slice4 Slice5
0
0.4
10

0.2
-1
10
Loss

u
10-2 -0.2 Slice1
Slice2
Slice3
-0.4 Slice4
-3
10 Slice5
Reference
-0.6
-4
10
5000 10000 15000 0 0.5 1 1.5
Epochs x

(a) Loss history with PINN (b) u at t = 1 at different epochs

6
Slice1
Slice2
Slice3
5 Slice4
Slice5

4
Residual

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x
(c) Residual of x at t = 1 at different epochs

Figure 2: Results of Burgers equation with PINN. NN with 4 hidden layers and 30 neurons per layer. The PDE residual points are set with
uniform 100 × 100 grids on X × T space and the IBs points are set with uniform 100 [Link] optimizer is Adam with a learn rate 0.001.
The reference result is given with WENO-Z on a refine mesh with 10000 spatial grids.

4
2.3. PINN-WE

Based on above analysis of Burgers equation, the training will fall into straggling due to the transition-points if we
take an average of each point residual into the total loss. We propose an weighted equation method to weaken the effect
of points in highly compressible regions by assigning a local positive weight λ to the governing equations. This method
is based on the fact that strong discontinuity is formed by the convergence of the characteristic lines. As we weaken
the expression of transition-points, the NN will focus on training the smooth regions and get a high accuracy in those
region, then automatically effected by the compression property of the strong discontinuity, the transition-points will be
compressed into smooth region. Then a sharp and exact discontinuity solution will appear with the training.
For a general conservative equation (Eq. 1),

∂U
Gnew := λ( + ∇ · F), (6)
∂t

and Gnew = 0 is WE. Then, Gnew = 0 has the same solutions as G = 0 if λ is always positive. In addition, we can adjust
the NN expression in different points by the design of gradient-dependent weight λ. Correspondingly, the new loss is
defined as
1 X 1 X
L= G2new,i + ε1 (Ûe,i − Ue0,i )2 , (7)
|SPDE | |SIBs |
xi ∈SPDE xi ∈SIBs

The architecture of PINN-WE for solving conservative hyperbolic equations is shown in Fig. 1.
We define the gradient-dependent weight as

1
λ= . (8)
ε2 (|∇ · ~u| − ∇ · ~u) + 1

where ~u is the velocity field. In [1], artificial viscosity is added into the governing equations according to the velocity
divergence ∇ · ~u considering that the field is compressed when ∇ · ~u < 0. Here, we only use the velocity divergence to
detect shocks and apply λ into the equations without adding any numerical dissipation to the equations. As λ is constant
positive, so it will not influence the exact solution of the equations.
Fig.1 shows the architecture of the proposed method to solve the given PDE (1). We evaluate our proposal in the
following numerical examples.

3. Numerical examples

3.1. Burgers equation

We first resolve Burgers equation (5) using the new PINN-WE method with a totally same setting with section 2.2.
In Burgers equation the velocity is u itself. So the weight in Gnew is

1
λ= . (9)
ε2 (| ∂u
∂x | − ∂u
∂x ) +1

Fig.3 gives the loss history of epochs and u, the residual and λ as functions of x at t = 1 at different training epochs
(1000,3000,5000,8000,11500 and 14931 (which has the minimum total loss)). Compared to the result in Fig.2, in the first
1000 epochs (Slice 1), PINN-WE has a similar result with traditional PINN as the gradient of transition-point is not
large enough to cause problems. However, after 3000 epochs, PINN-WE can effectively deduce the total loss. All the
transition-points are compressed into smooth region besides the one u = 0 which has zero velocity. As shown in Fig.3.d,
after 3000 epochs, all the weights change to nearly 1 besides the central one which has large gradient.

1 Code can be found in https : //[Link]/bf ly123/P IN N W E

5
 Loss_PDE
Loss_IBs Slice1
0.6
10
0 Loss_Total Slice2
Slice1 Slice2 Slice3
Slice3
Slice4 Slice5 Slice6
Slice4
0.4 Slice5
10
-1 Slice6
Reference
0.2

10-2
Loss

u
0

-3
10 -0.2

-0.4
10-4

-0.6
10-5
5000 10000 15000 0 0.5 1 1.5 2
Epochs x

(a) Loss history (b) u at t = 1 at different epochs

x
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
6 Slice1 1
Slice2
Slice3
Slice4
5 Slice5
Slice6 0.8

4
Slice1
Residual

Slice2
0.6 Slice3
Weight

3 Slice4
Slice5
Slice6
2 0.4

1
0.2

0
0 0.5 1 1.5 2
x

(c) Residual of x at t = 1 at different epochs (d) Weight λ of x at t = 1 at different epochs

Figure 3: Results of Burgers equation with PINN-WE. Same computational setting with Fig.2

6
3.2. Euler equations

Then we consider the following 1D and 2D Euler equations in conservative forms:

∂U
+ ∇ · F = 0. (10)
∂t

For the 1D case,    

ρ ρu
   
   
2
U = ρu , F = ρu + p  . (11)
  
   
E u(E + p)
For 2D cases, F = (F1 , F2 ), where
     
ρ ρu ρv
     
   2   
ρu ρu + p  ρuv 
U =   , F1 =   , F2 =  . (12)
     
 2
ρv  ρuv ρv + p 
    

     
E u(E + p) v(E + p)

Above, ρ is the density, u is the velocity, (u, v) is the 2D velocity vector, p is the pressure, and E is the total energy. To
close the equations, we use the following equation of state for ideal gas:

1 2 p
E= ρu + , (13)
2 γ−1

where γ = 1.4 is the specific heat ratio.

We consider forward problems with PINN to solve weighted Euler equations Gnew = 0 directly without more data
besides the IBs. In 1D cases, we consider the classical Sod and Lax problems. Then, we consider a 2D Riemann problem
and a more complex moving shock problem. We compare all the results with those obtained from traditional high-order
finite differential methods, namely, fifth-order WENO-Z method in spatial discretization [20] and third-order Runge-Kutta
method [21] for time integration.

3.2.1. Sod problem

The Sod problem has been extensively studied. It is a 1D Riemann problem with the initial constant states in a tube
with unit length formulated as 
(1, 0, 1), if 0 ≤ x ≤ 0.5,
(ρ, u, p) = (14)
(0.1, 0, 0.125), if 0.5 < x ≤ 1.


We used an NN with 7 hidden layers and 50 neurons per layer. Then, we randomly selected 10,000 residual points
from a uniform mesh of 100 × 200 in the X × T space. The number of IBs points was 1000. After training, we constructed
a test set with 100 uniform points in x ∈ [0, 1] at final time t = 0.2.
We first compared the results with the traditional high-order WENO-Z method for 100 cells in space. Fig. 4 shows
that PINN-WE achieves similar and even better results compared with the WENO-Z method. Especially in capturing
shock waves, no transition-points occur inside the shock because no dissipation is introduced into the equations for
PINN-WE.
Fig. 3 shows the training loss of PINN and PINN-WE, and Fig. 4 shows results at different training epochs to evaluate
the training evolution. In the first 1000 epochs, there is a small difference between PINN and PINN-WE because no
strong compression occurs. Then, PINN straggles when training a shock, while PINN-WE can easily avoid this region.
With the loss decreasing at smooth regions, a sharp shock appears with the transition-points compressed by the left and
right smooth regions.

7
 1 Velocity 1
Velocity

0.8 0.8

0.6 0.6

Density Density
0.4 0.4
Pressure Pressure

0.2 0.2
PINN-WE WENO-Z
Exact Exact

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x

Figure 4: Evaluation of Sod problem at t = 0.2.

101
ADAM

100
Total loss

-1
10

-2 LBFGS
10
PINNs
PINNs-WE

-3
10

0 2000 4000 6000 8000 10000

Epochs

Figure 5: Training loss evolution for Sod problem.

8
 1000 Epochs 2000 Epochs 4000 Epochs

PINNs-WE
PINNs PINNs-WE PINNs-WE
PINNs PINNs
0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x x

1 6000 Epochs 8000 Epochs 10000 Epochs

0.8

0.6

0.4

0.2
PINNs-WE PINNs-WE
PINNs-WE PINNs
PINNs PINNs
0
0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x x

Figure 6: Results at different training epochs for Sod problem.

3.2.2. Lax problem

The Lax problem is also a Riemann problem that contains a strong shock and strong contact. Its initial condition is
given by 
(0.445, 0.698, 3.528), if 0 ≤ x ≤ 0.5,
(ρ, u, p) = (15)
(0.5, 0, 0.571), if 0.5 < x ≤ 1.


We used the same NN as for the Sod problem and randomly selected 50,000 interior points from a uniform 1000 × 5000
mesh in the X × T space. We considered 1000 initial points. This problem was more difficult to train than the Sod
problem, and the total loss was 0.01 after training for 20,000 epochs. We compared the results with those obtained from
the high-order WENO-Z method. Fig. 5 shows that PINN-WE correctly captures shocks.

3.2.3. 2D Riemann problem

Then, we considered a 2D problem with strong contact discontinuities. The basic settings were retrieved from case 8
in [22]. The initial condition is given by

(1, −0.75, 0.5, 1)

 if 0 ≤ x ≤ 0.5, 0 ≤ y ≤ 0.5,



(2, 0.75, 0.5, 1) if 0 ≤ x ≤ 0.5, 0.5 < y ≤ 1,

(ρ, u, v, p) = (16)
(3, −0.75, −0.5, 1) if 0.5 < x ≤ 1, 0 ≤ y ≤ 0.5,







(1, 0.75, −0.5, 1) if 0.5 ≤ x ≤ 0.5, 0.5 < y ≤ 1.


We used an NN with 6 hidden layers and 60 neurons per layer. The training points were obtained by Latin hypercube
sampling, with 200,000 interior points in the T × X × Y space and 10,000 initial points in the X × Y space. The final
training loss was 0.009. We evaluated the model with meshes of 100 × 100 and 400 × 400 in the X × Y space at instants
0.2 and 0.4, respectively. Then, we compared the results with those obtained from the WENO-Z scheme computed in a
100 × 100 mesh in the X × Y space. The corresponding results are shown in Fig. 6.
We provide test results for 100 × 100 mesh points, which outnumber the training points (approximately 60) along

9
 PINNs-WE WENO-Z
Exact Exact

Pressure Pressure

Velocity Velocity

Density Density

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

x x

Figure 7: Results and loss for Lax problem.

each dimension. Comparison results are shown in Figs. 8 and 9 with the WENO-Z method in the same 100 × 100
mesh. For clarity, we show the original data without using the smoothing effect of the plotting software (Tecplot). The
proposed method can capture contact discontinuities more sharply and nearly without transport and smoothing points,
which are unavoidable by the traditional high-order method. Then, we increased the test points to a 400 × 400 mesh,
which is substantially larger than training data along each dimension. Then, we performed an unfair comparison with
the WENO-Z method in a 400 × 400 mesh. The computation of the detailed structure, especially in the middle region,
is weaker in PINN-WE owing to the available training data, but discontinuities are still sharper than those computed by
WENO-Z. These results illustrate the advantages of PINN-WE in high dimensions given its meshless feature.

3.2.4. Transonic flow through circular cylinder

Finally, we considered a 2D example with strong shock. A transonic flow with Mach number 0.728 passed through a
stationary circular cylinder. A similar problem and its detailed analysis can be found in [23]. The initial condition was
given by a uniform flow with (ρ, u, v, p) = (2, 112, 1.028, 0, 3.011). The computational domain was T = [0, 0.4],X = [0, 1.5]
and Y = [0, 2]. The center of the cylinder of radius 0.25 was located at (1, 1). We used an NN with 7 hidden layers
and 90 neurons per layer. The collocations were obtained from Latin hypercube sampling in 3D domain T × X × Y for
300, 000 points. In addition, 15, 000 boundary points were randomly sampled on the cylinder, and 15, 000 initial points
were obtained from Latin hypercube sampling. The total loss was 0.028 after 2000 steps of optimization with the L-BFGS
algorithm.
Figs. 10–14 show the comparison results with the WENO-Z method. Sharp and similar but smoother results are
obtained from PINN-WE compared with WENO-Z. Moreover, our proposal can clearly capture the vertex behind the
cylinder. However, limited by the number of the collocation samples owing to the hardware we used, we could not
further decrease the loss to obtain a more exact result of small structures. Nevertheless, PINN-WE seems promising for
computing complex transonic and supersonic flows.
The results at time 0.4 in Fig. 6 show that the proposed PINN-WE is also effective for relatively complex high-
dimensional problems. Compared with the results from the high-order WENO-Z method with 100 × 100 grids, more
sharp shocks are obtained using our proposal.

10
 (a) Density PINN-WE (b) Density WENO-Z

(c) U PINN-WE (d) U WENO-Z

Figure 8: Results of 2D Riemann problem (part 1) with 100 × 100 test points for PINN-WE and the same number of mesh grids for WENO-Z.

11
 (a) V PINN-WE (b) V WENO-Z

(c) Pressure PINN-WE (d) Pressure WENO-Z

Figure 9: Results of 2D Riemann problem (part 2) with 100 × 100 test points for PINN-WE and the same number of mesh grids for WENO-Z.

12
 (a) Density PINN-WE (b) Density WENO-Z

(c) U PINN-WE (d) U WENO-Z

Figure 10: Results of 2D Riemann problem (part 3) with 400×400 test points for PINN-WE and the same number of mesh grids for WENO-Z.

13
 (a) V PINN-WE (b) V WENO-Z

(c) Pressure PINN-WE (d) Pressure WENO-Z

Figure 11: Results of 2D Riemann problem (part 4) with 400 × 400 test points for PINN-WE and the same number of mesh grids for WENO-Z.

14
Figure 12: Resulting pressure for transonic flow through circular cylinder using proposed PINN-WE (left) and WENO-Z method (right).

Figure 13: Resulting density for transonic flow through circular cylinder using proposed PINN-WE (left) and WENO-Z method (right).

15
Figure 14: Resulting velocity u for transonic flow through circular cylinder using proposed PINN-WE (left) and WENO-Z method (right).

Figure 15: Resulting velocity v for transonic flow through circular cylinder using proposed PINN-WE (left) and WENO-Z method (right).

16
Figure 16: Resulting streamline for transonic flow through circular cylinder using proposed PINN-WE (left) and WENO-Z method (right).

4. Conclusions

We propose a concept to capture discontinuities, especially shock waves, using a PINN for solving Euler equations.
Unlike the idea of enhancing the NN expression in a large gradient domain, we consider the existence of a paradoxical
problem within wrong inside shock points (transition-points). Regardless of increasing or decreasing gradients, these
points likely increase the total loss. Thus, NN training may fall in conflict. Accordingly, we introduce a positive gradient-
dependent weight into the governing equations to adjust the expression of a PINN in regions with different physical
features. Then, for solving the Euler equations, we construct a weight inverse to the local physics compression by
measuring the velocity divergence.
By solving the WEs using PINN, training focuses on smooth regions, while shock regions have very small weights.
By relying on the real physics compression from the trained smooth regions, discontinuities automatically appear as the
transition-points move out into smooth regions like passive particles.
We only focused on capturing discontinuities in this study, but many questions for solving Euler equations remain
open. In fact, we found various problems that have not been addressed when simulating complex discontinuous problems.
First, convergence to a real physical weak solution is difficult. Using WE, PINN can converge fast to a discontinuous
result, but a nonphysical weak solution may be obtained for the Euler equations. Thus, effective entropy conditions and
other physical limitations should be integrated into PINN. Second, the shock position may be inaccurate without the
exact conservation of mass, momentum, and total energy in the PINN model. Third, as a global method, PINN can
suitably simulate big structures but may inadequately reflect detailed structures (Fig. 5). Thus, the resolution of PINN
should be improved to increase its applicability to diverse problems.

17
References

[1] J. VonNeumann, R. D. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, Journal of
applied physics 21 (3) (1950) 232–237.

[2] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order accurate essentially non-oscillatory
schemes, III, in: Upwind and high-resolution schemes, Springer, 1987, pp. 218–290.

[3] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted eno schemes, Journal of computational physics 126 (1)
(1996) 202–228.

[4] B. Cockburn, C.-W. Shu, The local discontinuous galerkin method for time-dependent convection-diffusion systems,
SIAM Journal on Numerical Analysis 35 (6) (1998) 2440–2463.

[5] S. Pirozzoli, Numerical methods for high-speed flows, Annual review of fluid mechanics 43 (2011) 163–194.

[6] D. Zhang, C. Jiang, D. Liang, L. Cheng, A review on TVD schemes and a refined flux-limiter for steady-state
calculations, Journal of Computational Physics 302 (2015) 114–154.

[7] 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, Journal of Computational physics
378 (2019) 686–707.

[8] G. Pang, L. Yang, G. E. Karniadakis, Neural-net-induced gaussian process regression for function approximation
and PDE solution, Journal of Computational Physics 384 (2019) 270–288.

[9] K. O. Lye, S. Mishra, D. Ray, Deep learning observables in computational fluid dynamics, Journal of Computational
Physics 410 (2020) 109339.

[10] J. Magiera, D. Ray, J. S. Hesthaven, C. Rohde, Constraint-aware neural networks for riemann problems, Journal of
Computational Physics 409 (2020) 109345.

[11] H. Huang, Y. Liu, V. Yang, Neural networks with local converging inputs (NNLCI) for solving conservation laws,
part II: 2D problems, arXiv preprint arXiv:2204.10424.

[12] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific machine learning through
physics-informed neural networks: Where we are and what’s next, arXiv preprint arXiv:2201.05624.

[13] R. G. Patel, I. Manickam, N. A. Trask, M. A. Wood, M. Lee, I. Tomas, E. C. Cyr, Thermodynamically consistent
physics-informed neural networks for hyperbolic systems, Journal of Computational Physics 449 (2022) 110754.

[14] Z. Mao, A. D. Jagtap, G. E. Karniadakis, Physics-informed neural networks for high-speed flows, Computer Methods
in Applied Mechanics and Engineering 360 (2020) 112789.

[15] 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, Computer Methods in Applied Mechanics and
Engineering 365 (2020) 113028.

[16] A. D. Jagtap, Z. Mao, N. Adams, G. E. Karniadakis, Physics-informed neural networks for inverse problems in
supersonic flows, arXiv preprint arXiv:2202.11821.

[17] A. Papados, Solving hydrodynamic shock-tube problems using weighted physics-informed neural networks with
domain extension, DOI: 10.13140/RG.2.2.29724.00642/1.

18
[18] Y. Shen, L. Liu, Y. Yang, Multistep weighted essentially non-oscillatory scheme, International Journal for Numerical
Methods in Fluids 75 (4) (2014) 231–249.

[19] J. Yu, L. Lu, X. Meng, G. E. Karniadakis, Gradient-enhanced physics-informed neural networks for forward and
inverse pde problems, Computer Methods in Applied Mechanics and Engineering 393 (2022) 114823.

[20] R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic
conservation laws, Journal of Computational Physics 227 (6) (2008) 3191–3211.

[21] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of
computational physics 77 (2) (1988) 439–471.

[22] A. Kurganov, E. Tadmor, Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem
solvers, Numerical Methods for Partial Differential Equations: An International Journal 18 (5) (2002) 584–608.

[23] H. Mo, F.-S. Lien, F. Zhang, D. S. Cronin, An immersed boundary method for solving compressible flow with
arbitrarily irregular and moving geometry, International Journal for Numerical Methods in Fluids 88 (5) (2018)
239–263.

19