# Physics-informed neural networks for inverse problems in supersonic flows

Ameya D. Jagtap1,†, Zhiping Mao2,†, Nikolaus Adams3, and George Em
Karniadakis1 1 Division of Applied Mathematics, Brown University, 182 George
Street, Providence, RI, 02912, USA  2 School of Mathematical Sciences Xiamen
University, Xiamen, Fujian 361005, China.  3 Department of Mechanical
Engineering, Technical University of Munich, 85748 Garching, Germany.

###### Abstract

Accurate solutions to inverse supersonic compressible flow problems are often
required for designing specialized aerospace vehicles. In particular, we
consider the problem where we have data available for density gradients from
Schlieren photography as well as data at the inflow and part of wall
boundaries. These inverse problems are notoriously difficult and traditional
methods may not be adequate to solve such ill-posed inverse problems. To this
end, we employ the physics-informed neural networks (PINNs) and its extended
version, extended PINNs (XPINNs), where domain decomposition allows to deploy
locally powerful neural networks in each subdomain, which can provide
additional expressivity in subdomains, where a complex solution is expected.
Apart from the governing compressible Euler equations, we also enforce the
entropy conditions in order to obtain viscosity solutions. Moreover, we
enforce positivity conditions on density and pressure. We consider inverse
problems involving two-dimensional expansion waves, two-dimensional oblique
and bow shock waves. We compare solutions obtained by PINNs and XPINNs and
invoke some theoretical results that can be used to decide on the
generalization errors of the two methods.

###### keywords:

Extended physics-informed neural networks, Entropy conditions, supersonic
compressible flows, inverse problems

††journal: Journal

##  1 Introduction

In recent years physics-informed neural networks (PINNs) [1] emerged as an
alternative simple method to solve many problems in computational science and
engineering, see, for example [2, 3, 4, 5, 6, 7, 8, 9, 10, 4, 11, 12, 13, 14,
15]. In particular, PINNs do not require meshes and can efficiently solve
forward problems and even ill-posed inverse problems, which are otherwise
difficult or sometimes even impossible to solve using traditional numerical
methods. Moreover, PINN can easily handle noisy, sparse and multi-fidelity
data sets. The main advantage of the PINN methodology is that it can
seamlessly incorporate all the given information like governing equation,
experimental data, initial/boundary conditions into the loss function, thereby
recasting the original problem into an equivalent optimization problem.
Recently, Shin et al. [16] established the mathematical foundation of PINNs
for linear partial differential equations, whereas in [17], Mishra and
Molinaro presented estimate on the generalization error of the PINN
methodology.

In this work, we consider inverse problems of the shock wave problems in
supersonic compressible flows. The governing equations of such flows are
compressible Euler equations, which admit discontinuities or shocks, even
though the initial states are smooth. Such ill-posed inverse problems are
difficult or even sometimes impossible to solve using the traditional
numerical solvers. Moreover, for the shock wave problem, the traditional
numerical methods usually require boundary conditions (BCs) for all field
variables. However, for the simulations of high-speed flows, traditional
numerical methods are not only expensive, but also it is usually difficult to
determine the BCs, for example, as mentioned in [18], it is particularly vital
to couple the surface BCs into the flow-field calculation for the expansion
wave problem. Far-field BCs or symmetric BCs were employed for the
compressible flows [19, 18], but applying these kinds of BCs results in large
computational domains, and therefore becomes computationally expensive. Other
BCs like numerical boundary condition may cause numerical instability [20]. We
investigate in the present work the possibility of using the idea of PINN to
resolve this issue by extending the algorithm from the one-dimensional case
[3] to the two-dimensional shock wave problems. Another advantage of PINNs is
that PINNs can be employed in any arbitrarily shaped domain of physical
interest, while the traditional numerical methods need a specified large
computational domain. In the PINNs literature, Mao et al. [3] first solved
both forward and inverse one and two-dimensional Euler equations involving
shocks in the Cartesian domain. Recently, Patel et al. [21] introduced
thermodynamically consistent physics-informed neural networks, which enforce
the entropy conditions for the scalar conservation laws and the one-
dimensional Euler equations. In [22] Fuks and Tchelepi proposed a stable way
for PINNs to handle scalar conservation laws involving discontinuities.
Moreover, Jagtap et al. proposed the spatial domain decomposition based PINN
approach for conservative laws called as conservative PINN (cPINN) [6], and
more general space-time domain decomposition based PINN approach for any type
of differential equations, named eXtended PINNs (XPINNs) [5]. Apart from
parallel implementation [8], such domain decomposition approaches can allow
local powerful networks in subdomains where complex solutions can be expected.
In [23] Hu et al. provide a prior generalization bound via the complexity of
the target functions in the partial differential equation (PDE) problem, and a
posterior generalization bound via the posterior matrix norms of the networks
after optimization. Based on these bounds, one can analyze the conditions
under which XPINNs improve generalization. The main feature of XPINN, namely
the spatio-temporal domain decomposition, introduces a generalization trade-
off. The XPINNs decompose the complex solution of differential equations into
several simple parts, which decreases the complexity needed to learn by the
networks in the respective subdomains and thereby boosts generalization. In
contrast, such decomposition leads to less training data being available in
each subdomain, and hence such a model is typically vulnerable to over-
fitting, due to which it may become less generalizable. In the present study
we used bounds in [23] for the decomposition of computational domain.

Apart from PINNs, the other frameworks which employed deep neural networks for
tackling shock wave problems can be found in [24, 25, 26, 27, 28].

(a) Attached waves: oblique shock and expansion waves

(b) Detached waves: bow shock waves

Figure 1: Schematic of the three types of shock waves. Left: Oblique shock
(regions ① and ④) and expansion waves (regions ② and ③), right: bow shock
waves.

Inverse problems in supersonic compressible flows are often encountered in
designing specialized high-speed vehicles. Apart from the governing Euler
equations, which admit multiple weak solutions, we also enforce the entropy
conditions that are satisfied by the viscosity solutions as well as positivity
conditions on density and pressure fields. In this work, we shall infer all
the states (density, velocity and pressure) for the shock wave problems,
including an oblique shock wave, expansion wave and bow shock wave (see Fig.
1), in the two-dimensional case. In particular, we use the idea of hidden
fluid mechanics based on PINNs to solve steady state problems. More precisely,
we use the information of density gradients, a small data set on primitive
variables, as well as global constraints to predict the density, velocity, and
pressure for these inverse problems. In particular, we use both PINNs and
XPINNs to infer unknown states involving shocks and expansion waves.

This paper is arranged as follows: in Section 2 we present the governing
compressible Euler equations and the entropy conditions. In Section 3 we
describe the methodology, including the mathematical setup for the feed-
forward neural networks, with a brief introduction to PINNs and XPINNs
methodology, and the optimization methods. In section 4, various computational
examples such as expansion wave, oblique shock wave and the bow shock wave are
solved and various comparisons are made. Finally, we summarize our findings in
Section 5.

##  2 Governing Equations

The governing compressible Euler equations in the conservative form [29, 30,
31] can be written as:

$\partial_{t}U+\nabla\cdot
G(U)=0,~{}~{}~{}\bm{x}\coloneqq\\{x,y\\}\in\Omega\subset\mathbb{R}^{2},~{}~{}t\in(0,T],$
(1)

where

$U=\left[{\begin{array}[]{c}\rho\\\ \rho u\\\ \rho v\\\ \rho E\\\
\end{array}}\right],\;G=\\{G_{1},G_{2}\\},\;\text{with
}G_{1}(U)=\left[{\begin{array}[]{c}\rho u\\\ p+\rho u^{2}\\\ \rho uv\\\
pu+\rho uE\\\ \end{array}}\right],\;G_{2}(U)=\left[{\begin{array}[]{c}\rho
v\\\ \rho uv\\\ p+\rho v^{2}\\\ pv+\rho vE\\\ \end{array}}\right].$

Here, $\rho$ is the density, $p$ is the pressure, $(u,\,v)$ are the velocity
components, and $E$ is the total energy. We use the additional equation of
state, which describes the relation of the pressure and energy, to close the
above Euler equations. In particular, we consider the equation of state for a
polytropic gas given by

$p=(\gamma-1)\left(\rho E-\frac{1}{2}\rho\|{\bm{u}}\|^{2}\right),$ (2)

where $\gamma$ is the ratio of specific heats, whose value for air is 1.4, and
$\bm{u}=(u,v)$. Here, we solve the steady state Euler equations given as

$\nabla\cdot G(U)=0,\;\bm{x}\in\Omega.$ (3)

Remark: It is important to note that for neural network training we employ the
non-dimensional form of Euler equations, which look exactly the same as the
dimensional ones.

###  2.1 Entropy conditions

In general, equation (1) admits several weak solutions [32] and thus
additional conditions have to be imposed to select the ‘physically’ relevant
solution, among the others, the so-called Entropy condition. The viscosity
solutions [33] of conservation laws are associated with the zero-viscosity
limit given by $U=\lim_{\epsilon\rightarrow 0^{+}}U_{\epsilon}$, where
$\partial_{t}U_{{\epsilon}}+\nabla\cdot
G(U_{\epsilon})=\epsilon\nabla^{2}U_{\epsilon}$. These viscosity solutions
also satisfy entropy conditions.

A pair $(\eta,\phi)$ with a convex function $\eta(U)$ and an associated
entropy flux (a vector valued function)
$\bm{\phi}(U)=[\phi_{1}(U),\phi_{2}(U)]^{T}$ are called entropy-entropy flux
pair [34], if it satisfies the relations

$\eta^{\prime}(U)G_{1}^{\prime}(U)=\phi^{\prime}_{1}(U),\ \ \ \ \ \
\eta^{\prime}(U)G_{2}^{\prime}(U)=\phi^{\prime}_{2}(U).$

Thus, an additional conservation law can be added to the governing
compressible Euler equations

$\eta_{t}+\phi_{1_{x}}+\phi_{2_{y}}=0$ (4)

for smooth solutions. At discontinuities, the entropy is not always conserved
thus the admissible solution is chosen based on the following inequality.

$\eta_{t}+\phi_{1_{x}}+\phi_{2_{y}}\leq 0.$ (5)

All the admissible, bounded and non-smooth solutions satisfy the Entropy
condition in a distributional sense. It is important to note that the given
hyperbolic conservation laws may not have a unique entropy pair [34]. For
Euler equations, the most common choice of an entropy-entropy flux pair is

$\eta\triangleq\frac{-\rho s}{\gamma-1}\ \ \text{and}\ \
\bm{\phi}\triangleq\frac{-\rho\bm{u}s}{\gamma-1}.$ (6)

where the specific entropy $s=\log(p/\rho^{\gamma})$.Note that during the
network training, we also enforced the positivity conditions on density and
pressure. This will avoid the problem of handling negative values for pressure
and density. This can be done by choosing the following output from the neural
networks max$(\alpha,\rho)$ and max$(\alpha,p)$ for density and pressure,
respectively. Here $0<\alpha<<1$ is a small positive constant.

##  3 Methodology

###  3.1 Mathematical setup for feed-forward neural network (FFNN)

Let $\mathcal{N}^{L}:\mathbb{R}^{N_{0}}\rightarrow\mathbb{R}^{N_{L}}$ be a
FFNN of $L$ layers and $N_{k}$ neurons in the $k^{th}$ layer, where the input
layer has $N_{0}$ and the output layer $N_{L}$ neurons. The weights and biases
in the $k^{th}$ layer ($1\leq k\leq L$) are given as
$\bm{W}^{k}\in\mathbb{R}^{N_{k}\times N_{k-1}}$ and
$\bm{b}^{k}\in\mathbb{R}^{N_{k}}$, respectively. If
$\bm{z}\in\mathbb{R}^{N_{0}}$ is the input vector at the $k^{th}$ layer, then
the output vector at $k^{th}$ layer is denoted by $\mathcal{N}^{k}(\bm{z})$.
Thus, the input vector can be expressed as $\mathcal{N}^{0}(\bm{z})=\bm{z}$.
We employ the layer-wise adaptive activation function, denoted by
$\sigma(na^{k})$, where $a^{k}$ are the trainable parameters, $n$ is the
predefined scaling factor and $\sigma$ is the activation function, see Jagtap
et al. [35, 36]. The introduction of the additional activation slope
parameters $a^{k}$ dynamically changes the slope of activation function during
training, thereby increases the training speed of the neural networks, see
[35, 36]. Recently, Jagtap et al. [37] proposed the Rowdy activation function
that is designed to eliminate any saturation region by injecting sinusoidal
noise like effects, which further increases the training speed. In this paper,
we employ layer-wise locally adaptive activation functions with scaling factor
$n=10$ for all hidden-layers and initialize $na^{k}=1,~{}\forall k$. The
$(L-1)$-hidden layer FFNN is defined as

$\displaystyle\mathcal{N}^{1}(\bm{z})=\bm{W}^{1}\bm{z}+\bm{b}^{1},\quad\text{
for }\quad k=1,$ (7)
$\displaystyle\mathcal{N}^{k}(\bm{z})=\bm{W}^{k}\sigma(a^{k-1}\mathcal{N}^{k-1}(\bm{z}))+\bm{b}^{k}\in\mathbb{R}^{N_{k}},\quad\text{
for }\quad 2\leq k\leq L\ ,$

where in the activation function is identity in the last layer. Let
$\tilde{\bm{\Theta}}=\\{\bm{W}^{k},\bm{b}^{k},a^{k}\\}_{k=1}^{L}\in\mathcal{V}$
be the collection of all trainable parameters including weights, biases, and
activation slopes, and $\mathcal{V}$ is the parameter space, then, the neural
network output can be written as

$\xi_{\tilde{\bm{\Theta}}}(\bm{z})=\mathcal{N}^{L}(\bm{z};\tilde{\bm{\Theta}})\
,$

where $\mathcal{N}^{L}(\bm{z};\tilde{\bm{\Theta}})$ emphasizes the dependence
of the neural network output $\mathcal{N}^{L}(\bm{z})$ on the trainable
parameters $\tilde{\bm{\Theta}}$. In general, weights and biases are
initialized from known probability distributions. In particular, we used the
Xavier initialization [38] for all the trainable parameters.

###  3.2 Overview of Physics-Informed Neural Networks (PINNs)

For a steady compressible inviscid flow governed by the Euler equations (3),
we assume that apart from the inflow and wall boundary data on partial or full
primitive variables, there is a set of experimental data (Schlieren data)
available on density gradient. By using such data and combining them with the
Euler equations (3), we would like to design PINNs to infer all the states of
interest, i.e., the density $\rho(x,y)$, the velocity $\bm{u}(x,y)$ and the
pressure $p(x,y)$. Fig. 2 shows the schematic representation of the PINNs for
the Euler equations. The left part is the feed-forward neural networks with
input $x$ and $y$ and the primitive variables as the output. The right part is
the physics-informed part where network output is forced to satisfy
conservation laws. Both parts contribute to loss function. In the present
work, we shall learn the density, velocity, and pressure fields by using a
feed-forward neural network. In particular, we shall approximate the density,
velocity and pressure by setting the neural network prior on these primitive
variables denoted jointly as $\xi_{\tilde{\bm{\Theta}}}(x,t)$.

Figure 2:  Schematic representation of the PINNs for the Euler equations. The
left part is the 2 hidden-layer deep feed-forward neural networks with input
$x$ and $y$, and with the primitive variables as the output. Here,
$\sigma(\cdot)$ is the nonlinear activation function, whereas $a$ is the
tunable activation slope. The right part is the physics-informed part, where
the network output is forced to satisfy conservation laws as well as entropy
conditions. Both parts contribute to loss function. Note that during the
network training the positivity condition on the density and the pressure is
enforced by choosing the maximum value between them and the small positive
constant $\alpha$.

Let $\\{\bm{x}^{(i)}_{\xi}\\}_{i=1}^{N_{\xi}}$ and
$\\{\bm{x}^{(i)}_{F}\\}_{i=1}^{N_{F}}$ be the set of randomly selected
training data and residual points, respectively. These points are drawn from a
distribution, which is usually not known a priori and often needs to be chosen
from the given input training data. The PINN algorithm aims to learn a
surrogate $\xi\approx\xi_{\tilde{\bm{\Theta}}}$ for predicting the solution
$\xi$ of the given PDE. The PINN loss function is given as

$\mathcal{J}(\tilde{\bm{\Theta}})=W_{\xi}~{}[\text{MSE}_{\xi}(\tilde{\bm{\Theta}};\\{\bm{x}^{(i)}_{\xi}\\}_{i=1}^{N_{\xi}})+\text{MSE}_{\nabla\rho}(\tilde{\bm{\Theta}};\\{\bm{x}^{(i)}_{\nabla\rho}\\}_{i=1}^{N_{\nabla\rho}})]+W_{\mathcal{F}}~{}\text{MSE}_{\mathcal{F}}(\tilde{\bm{\Theta}};\\{\bm{x}^{(i)}_{F}\\}_{i=1}^{N_{F}})\
,$ (8)

where $W_{\xi}$ and $W_{\mathcal{F}}$ are the weights for the data mismatch
and residual terms, respectively. The mean squared errors (MSEs) are given by

$\displaystyle\text{MSE}_{\xi}(\tilde{\bm{\Theta}};\\{\bm{x}^{(i)}_{\xi}\\}_{i=1}^{N_{\xi}})$
$\displaystyle=\frac{1}{N_{\xi}}\sum_{i=1}^{N_{\xi}}\left|\xi^{(i)}-\xi_{{\tilde{\bm{\Theta}}}}(\bm{x}^{(i)}_{\xi})\right|^{2}\
,\ $
$\displaystyle\text{MSE}_{\mathcal{F}}(\tilde{\bm{\Theta}};\\{\bm{x}^{(i)}_{F}\\}_{i=1}^{N_{F}})$
$\displaystyle=\frac{1}{N_{F}}\sum_{i=1}^{N_{F}}\left|\mathcal{F}_{{\tilde{\bm{\Theta}}}}(\bm{x}^{(i)}_{F})\right|^{2},$
$\displaystyle\text{MSE}_{\nabla\rho}(\tilde{\bm{\Theta}};\\{\bm{x}^{(i)}_{\nabla\rho}\\}_{i=1}^{N_{\nabla\rho}})$
$\displaystyle=\frac{1}{N_{\nabla}}\sum_{i=1}^{N_{\nabla}}|\nabla\rho^{i}-\nabla\rho_{{\tilde{\bm{\Theta}}}}(\bm{x}^{i}_{\nabla})|^{2}.$

The term $\text{MSE}_{\xi}$ is the MSE for data mismatch term for primitive
variables, which enforces the given experimental and/or synthetic data.
$\text{MSE}_{\nabla\rho}$ represent the data mismatch term for density
gradient. The term $\text{MSE}_{\mathcal{F}}$ is the MSE for the governing PDE
residue where
$\mathcal{F}_{{\tilde{\bm{\Theta}}}}(\bm{x}_{F})\coloneqq\mathcal{F}(\xi_{{\tilde{\bm{\Theta}}}}(\bm{x}_{F}))$
represents the residual of the governing system of conservation laws given by
equation (3) and the entropy conditions. The notations $N_{\xi}$, $N_{F}$ and
$N_{\nabla}$ are the number of training data, residual points and number of
data points for density gradient, respectively. The trainable parameters of
the neural networks $\xi_{\tilde{\bm{\Theta}}}$ can be estimated by minimizing
the loss function given by Eq. (8).

The PDE residual $\mathcal{F}_{{\tilde{\bm{\Theta}}}}$ in the loss function is
constructed using automatic differentiation [39], which is a graph-based
differentiation method for calculating derivatives of field variables
accurately in a computational graph. Moreover, compared to traditional
numerical differentiation, automatic differentiation does not suffer from
truncation and round-off errors. Highly popular machine learning libraries
such as TensorFlow [40] and Pytorch [41] provide automatic differentiation
package.

###  3.3 Extended Physics-Informed Neural Networks (XPINNs)

The extended physics-informed neural networks (XPINNs) methodology [5] is a
recently developed generalized space-time domain decomposition approach for
deep learning of PDEs. It overcomes many limitations of the vanilla PINN
method such as parallel implementation capacity, deployment of locally
powerful networks, large representation capacity, etc. Here, we shall briefly
discuss the XPINN method.

The computational domain $\Omega$ is divided into $N_{sd}$ number of
arbitrary, non-overlapping subdomains. Let
$\\{\bm{x}^{(i)}_{\xi_{q}}\\}_{i=1}^{N_{\xi_{q}}}$,
$\\{\bm{x}^{(i)}_{F_{q}}\\}_{i=1}^{N_{F_{q}}}$ and
$\\{\bm{x}^{(i)}_{I_{q}}\\}_{i=1}^{N_{I_{q}}}$ be the set of randomly selected
training, residual, and the common interface points, respectively in the
$q^{th}$ subdomain. The $N_{\xi_{q}},N_{F_{q}}$, and $N_{I_{q}}$ represent the
number of training data points, the number of residual points, and the number
of points on the common interface in the $q^{th}$ subdomain, respectively.
Similar to PINN, the XPINN algorithm aims to learn a surrogate
$\xi_{q}\approx\xi_{\tilde{\bm{\Theta}}_{q}}$, $q=1,2,\cdots,N_{sd}$ for
predicting the solution
$\xi\approx\xi_{\tilde{\bm{\Theta}}}=\sum_{q=1}^{N_{sd}}\xi_{\tilde{\bm{\Theta}}_{q}}\cdot
1_{q}$ of the given PDE over the entire computational domain. Here $1_{q}$
denotes the indicator function; see [5] for more details. The loss function of
XPINN is defined subdomain-wise, which has a similar structure as the PINN
loss function in each subdomain but is endowed with the interface conditions
for stitching the subdomains together. For the forward problem, the loss
function in the $q^{th}$ subdomain is defined as

$\displaystyle\mathcal{J}(\tilde{\bm{\Theta}}_{q})=~{}$ $\displaystyle
W_{\mathcal{F}_{q}}~{}(\text{MSE}_{G}^{Mass}+\text{MSE}_{G}^{MomX}+\text{MSE}_{G}^{MomY}+\text{MSE}_{G}^{Ene}+\text{MSE}_{G}^{Entropy})$
$\displaystyle+$ $\displaystyle
W_{\xi_{q}}~{}(\text{MSE}_{\xi}^{D}+MSE_{\nabla\rho})+W_{I_{q}}~{}\text{MSE}_{\xi_{avg}}+W_{I_{\mathcal{F}_{q}}}~{}\text{MSE}_{\mathcal{R}},~{}~{}q=1,2,.\cdots,N_{sd},$
(9)

where the MSE for mass, momentum ($x$ and $y$ directions), energy conservation
laws, and entropy condition are given as

$\displaystyle\text{MSE}_{G}^{Mass}$
$\displaystyle=\frac{1}{N_{F_{q}}}\sum_{i=1}^{N_{F_{q}}}|(G_{1}^{Mass}(\bm{x}^{i}_{F_{q}}))_{x}+(G_{2}^{Mass}(\bm{x}^{i}_{F_{q}}))_{y}|^{2},$
$\displaystyle\text{MSE}_{G}^{MomX}$
$\displaystyle=\frac{1}{N_{F_{q}}}\sum_{i=1}^{N_{F_{q}}}|(G_{1}^{MomX}(\bm{x}^{i}_{F_{q}}))_{x}+(G_{2}^{MomX}(\bm{x}^{i}_{F_{q}}))_{y}|^{2},$
$\displaystyle\text{MSE}_{G}^{MomY}$
$\displaystyle=\frac{1}{N_{F_{q}}}\sum_{i=1}^{N_{F_{q}}}|(G_{1}^{MomY}(\bm{x}^{i}_{F_{q}}))_{x}+(G_{2}^{MomY}(\bm{x}^{i}_{F_{q}}))_{y}|^{2},$
$\displaystyle\text{MSE}_{G}^{Ene}$
$\displaystyle=\frac{1}{N_{F_{q}}}\sum_{i=1}^{N_{F_{q}}}|(G_{1}^{Ene}(\bm{x}^{i}_{F_{q}}))_{x}+(G_{2}^{Ene}(\bm{x}^{i}_{F_{q}}))_{y}|^{2},$
$\displaystyle\text{MSE}_{G}^{Entropy}$
$\displaystyle=\frac{1}{N_{F_{q}}}\sum_{i=1}^{N_{F_{q}}}|(-\phi_{1}^{Ene}(\bm{x}^{i}_{F_{q}}))_{x}+(-\phi_{2}^{Ene}(\bm{x}^{i}_{F_{q}}))_{y}+\epsilon|^{2},$

where $\epsilon$ is a small positive constant. The MSE for non-zero Dirichlet
boundary condition is given as

$\displaystyle\text{MSE}_{\xi}^{D}$
$\displaystyle=\frac{1}{N_{\xi_{q}}}\sum_{i=1}^{N_{\xi_{q}}}|\xi^{i}-\xi_{{\tilde{\bm{\Theta}}}_{q}}(\bm{x}^{i}_{\xi_{q}})|^{2}.$

Along the common interface, the MSE average term for all primitive variables
represented by a common symbol $\xi$ and residual continuity condition
($\text{MSE}_{\mathcal{R}}$) are given as

$\displaystyle\text{MSE}_{\xi_{avg}}(\tilde{\bm{\Theta}}_{q};\\{\bm{x}^{(i)}_{I_{q}}\\}_{i=1}^{N_{I_{q}}})$
$\displaystyle=\sum_{\forall
q^{+}}\left(\frac{1}{N_{I_{q}}}\sum_{i=1}^{N_{I_{q}}}\left|\xi_{{\tilde{\bm{\Theta}}}_{q}}(\bm{x}^{(i)}_{I_{q}})-\left\\{\\!\left\\{\xi_{{\tilde{\bm{\Theta}}}_{q}}(\bm{x}^{(i)}_{I_{q}})\right\\}\\!\right\\}\right|^{2}\right),$
$\displaystyle\text{MSE}_{\mathcal{R}}(\tilde{\bm{\Theta}}_{q};\\{\bm{x}^{(i)}_{I_{q}}\\}_{i=1}^{N_{I_{q}}})$
$\displaystyle=\sum_{\forall
q^{+}}\left(\frac{1}{N_{I_{q}}}\sum_{i=1}^{N_{I_{q}}}\left|\mathcal{F}_{{\tilde{\bm{\Theta}}}_{q}}(\bm{x}^{(i)}_{I_{q}})-\mathcal{F}_{{\tilde{\bm{\Theta}}}_{q^{+}}}(\bm{x}^{(i)}_{I_{q}})\right|^{2}\right),$

respectively. Notation $\\{\\!\\{\cdot\\}\\!\\}$ represent the average value
of primitive variables along the common interface. The term
$\text{MSE}_{\nabla\rho}=\frac{1}{N_{\nabla_{q}}}\sum_{i=1}^{N_{\nabla_{q}}}|\nabla\rho^{i}-\nabla\rho_{{\tilde{\bm{\Theta}}}_{q}}(\bm{x}^{i}_{\nabla_{q}})|^{2}$
represents the data mismatch term for density gradient obtained from the
Schlieren photography experimental technique traditionally used in the study
of high-speed aerodynamics.

As discussed in [5], the XPINN can reduce the generalization error by
carefully selecting residual points in each subdomain. Moreover, with the
XPINN methodology, the computational domain can be decomposed in any arbitrary
way according to the required number and the location of residual points. This
can be achieved with some prior knowledge about the solution behavior such
that the localized powerful neural network can be employed in that subdomain.
The XPINN can also reduce the approximation error by selecting the network
size, and other hyperparameters in each subdomain depending on the complexity
of the solution. Furthermore, XPINN can also reduce optimization error by
using the different network size for each subdomain.

###  3.4 Optimization Methods

By minimizing the loss function, we seek the optimal parameters of the network
$\tilde{\bm{\Theta}}^{*}$ from the parameter space. For this purpose there are
several optimization algorithms available that can be used to minimize loss
function. The gradient-based methods are one of the most popular class of
optimization methods that can be employed for the minimization of loss
functions. In the basic form, given an initial value of tunable parameters
$\tilde{\bm{\Theta}}$, these parameters are updated as
$\tilde{\bm{\Theta}}^{m+1}=\tilde{\bm{\Theta}}^{m}-\eta_{l}\left.\frac{\partial\mathcal{J}(\tilde{\bm{\Theta}})}{\partial\tilde{\bm{\Theta}}}\right|_{\tilde{\bm{\Theta}}=\tilde{\bm{\Theta}}^{m}}$,
where the learning rate is denoted by $\eta_{l}$. In particular, we shall use
the Adam optimizer [42] followed by limited-memory
Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimizer [43], which is a quasi-
Newton optimization method. For smooth and regular solutions the L-BFGS
optimizer can find a better solution with a small number of iterations
compared to the Adam optimizer, due to second-order accuracy as opposed to
Adam, which is first-order accurate but in general more robust. Hence, we
often start with the Adam optimizer and after we reach a small value of the
loss function we switch to L-BFGS.

##  4 Computational results

In this section we shall solve several examples including expansion and shock
waves problems. In appendix A we present a pedagogical example involving the
inverse problem for the two-dimensional Euler equations with smooth solutions.
We provide all the details for this simple example so that the interested
reader can repeat our results.

In this work, the comparison study has been performed between the PINN and
XPINN methods based on the following points.

  * 1.

In both PINNs and XPINNs, we used unity weights for the different MSE terms in
the loss function for all the problems, except for the case where dynamics
weights are employed.

  * 2.

It is important to note that in this study we do not aim to compare the
computational cost associated with PINNs and XPINNs, because such study has
already been performed in [8]. The main objective of this work is to solve the
ill-posed inverse problems in supersonic flow, which has many real-world
applications in designing aerospace vehicles.

  * 3.

We are solving an inverse problem with given data inside the domain as well as
on the boundary of the domain. Such data availability allows us to decompose
the computational domain into several smaller subdomains as per the
compleixity of the data, especially for the XPINN method.

  * 4.

For all the test cases, the primitive variable data, which is available in the
dimensional form is converted into the non-dimensional one for the network
training.

###  4.1 Expansion wave problem

A supersonic expansion fan is a centered expansion process that occurs when a
supersonic flow turns around a convex corner. The expansion fan consists of an
infinite number of Mach waves, diverging from a sharp corner; see the left
plot of Figure 3 for the schematic. To overcome the aforementioned boundary
conditions issue, we are going to use PINNs and XPINNs to solve the expansion
wave problem with the inlet flow conditions as well as the density gradient
information from Schlieren photography; see Figure 3 (right).

We solve the inverse problem of the expansion wave problem by using PINN
without using the outflow boundary conditions. Instead, we use the information
of density gradient $\nabla\rho$ and a small number of pressure data points on
the surface of the wedge, following the one-dimensional problem [3]. In
addition, we use the inlet conditions here. However, unlike the one-
dimensional case where we use the density gradient in the entire computational
domain, here we only use the density gradients in a small region of the
computational domain. By combining the mathematical model, we have the
weighted loss function of PINN given by

$\mathcal{J}=\omega_{1}~{}\text{MSE}_{\mathcal{F}}+\omega_{2}~{}\text{MSE}_{\nabla\rho|_{D}}+\omega_{3}~{}\text{MSE}_{Inflow}+\omega_{4}~{}\text{MSE}_{p^{*}},$

where $\text{MSE}_{\mathcal{F}}$ corresponds to the Euler equation (3)
including the entropy conditions, $\text{MSE}_{Inflow}$ corresponds to the
inflow boundary conditions on primitive variables,
$\text{MSE}_{\nabla\rho|_{D}}$ corresponds to the density gradient and
$D\subset\Omega$ is a small region of the computational domain,
$\text{MSE}_{p^{*}}$ corresponds to the pressure data on the wall surface. We
also employ dynamic weights [44] denoted as $\omega_{i},i=1,2,\cdots$ in front
of all the MSE terms in the loss function. Hereafter, for brevity, the loss
function $\mathcal{J}$ is denoted without showing its dependence on the
network parameters $\tilde{\bm{\Theta}}$.

Consider the following example: Let $\gamma=1.4$. Assume the inlet flow for a
tuning corner with $\theta=10^{\circ}$ is given by

$M_{\infty}=2,\;\rho_{\infty}=1.23kg/m^{3},\;u_{\infty}=678.1m/s,\;v_{\infty}=0,\;p_{\infty}=1.01\times
10^{5}Pa,\;T_{\infty}=286.1K.$

The exact solutions can be obtained by using the isentropic relations [18].
The distributions of the residual points and the data points of the density
gradient are shown in Fig. 3 (right). For the pressure, we use 50 randomly
distributed data points on the surface of the wedge, i.e.,
$y=-tanh(\theta)x,\,x\in(0,1).$

Figure 3: Expansion wave problem: (Left) schematic of the expansion waves, and
(right) distributions of the 1200 residual points (blue $\ast$ points, where
we compute the residuals), the 340 data points for the density gradient (red
$+$ points) and the 100 data points for the inflow conditions (magenta $\circ$
points).

Besides using the typical hyperbolic tangent activation function and fixed
weights, we additionally investigate the effectiveness of using the adaptive
‘tanh’ activation function [45]. For PINNs, we use 6 hidden-layers with 40
neurons in each layer for all cases. The results are shown in Figs. 5-8. We
observe that all the predictions of the density obtained with or without
adaptive activation function/dynamic weights are in good agreement with the
exact solution. However, for the predictions of the pressure and the velocity,
the PINN solutions obtained without activation functions or dynamic weights
are not as good as the ones obtained with activation functions or dynamic
weights, especially for the velocity (see Figs. 7-8). This means that we
predict all the states by using PINNs with adaptive activation function or
dynamic weights along with the above setup, and the use of adaptive activate
functions or dynamical weights helps to improve the accuracy of PINNs. This is
very important when we consider the oblique shock or bow shock wave problems
later. Additionally, we observe that the results obtained by using dynamical
weights have supreme accuracy among all the results of PINNs. This is also
verified by presenting the relative $L_{2}$ error (see Table 1). Specifically,
by using the above setup, we obtain good predictions with PINN without using
the outflow boundary condition and the wall boundary condition. Moreover, we
verify the condition

$\tilde{n}\cdot{\bm{u}}=0$ (10)

on the wall, where $\tilde{n}$ is the unit outer normal vector. This will be
used as additional constraint for the bow shock problem considered later.

We also solved this problem using the XPINN methodology, and compared it with
PINN results. For that, we divide the domain into two subdomains as shown in
Fig. 4, keeping the network architecture same as PINNs in these subdomains.
The number of evenly distributed interface points on both interfaces is 300.

Figure 4: Expansion wave problem: Schematic presentation of subdomains and
interfaces for the XPINN method.

The adaptive activation functions are used for the XPINN method. It can be
seen that the XPINN results are slightly better than PINN. The relative
$L_{2}$ errors in all the primitive variables for all cases are given in Table
1. In the case of XPINNs, we present the total relative error from all the
subdomains. We also plot the loss function versus iterations for PINNs; see
Fig. 9 (left). Fig. 9 (right plot) shows the evolution of the dynamic weights
($\omega_{2},\omega_{3},\omega_{4}$) with the number of iterations for PINNs.
Fig. 10 shows the total loss function plot for XPINNs (left) and the combined
entropy loss (right) from all the subdomains.

(a) Exact solution

(b) PINN solution with DW (c) Relative error with DW (d) Relative error with
AA (e) Relative error without DW/AA (f) Relative error with XPINN.

Figure 5: Expansion wave problem: Comparison of the _density_ between PINN
solutions with/without adaptive activation function (AA) or dynamic weights
(DW). Here we use ‘tanh’ as the basic activation function. (a) Exact solution.
(b) PINN solution using dynamic weights. (c) Relative point-wise error with
dynamic weights. (d) Relative point-wise error with adaptive activation
function. (e) Relative point-wise error without dynamic weights or adaptive
activation function. (f) Relative point-wise error with XPINN using adaptive
activations.

(a) Exact solution

(b) PINN solution with DW (c) Relative error with DW (d) Relative error with
AA (e) Relative error without DW/AA (f) Relative error with XPINN.

Figure 6: Expansion wave problem: Comparison of the _pressure_ between PINN
solutions with/without adaptive activation function (AA) or dynamic weights
(DW). Here we use ‘tanh’ as the basic activation function. (a) Exact solution.
(b) PINN solution using dynamic weights. (c) Relative point-wise error with
dynamic weights. (d) Relative point-wise error with adaptive activation
function. (e) Relative point-wise error without dynamic weights or adaptive
activation function. (f) Relative point-wise error with XPINN using adaptive
activations.

(a) Exact solution

(b) PINN solution with DW (c) Relative error with DW (d) Relative error with
AA (e) Relative error without DW/AA (f) Relative error with XPINN.

Figure 7: Expansion wave problem: Comparison of the _velocity $u$_ between
PINN solutions with/without adaptive activation function (AA) or dynamic
weights (DW). Here we use ‘tanh’ as the basic activation function. (a) Exact
solution. (b) PINN solutions using dynamic weights. (c) Relative point-wise
error with dynamic weights. (d) Relative point-wise error with adaptive
activation function. (e) Relative point-wise error without dynamic weights or
adaptive activation function. (f) Relative point-wise error with XPINN using
adaptive activations.

(a) Exact solution

(b) PINN solution with DW (c) Relative error with DW (d) Relative error with
AA (e) Relative error without DW/AA (f) Relative error with XPINN.

Figure 8: Expansion wave problem: Comparison of the _velocity $v$_ between PINN solutions with/without adaptive activation function (AA) or dynamic weights (DW). Here we use ‘tanh’ as the basic activation function. (a) Exact solution. (b) PINN solutions using dynamic weights. (c) Relative point-wise error with dynamic weights. (d) Relative point-wise error with adaptive activation function. (e) Relative point-wise error without dynamic weights or adaptive activation function. (f) Relative point-wise error with XPINN using adaptive activations. | PINNs without AA/DW | PINNs with AA | PINNs with DW | XPINNs with AA
---|---|---|---|---
$\rho$ | 1.3494e-02 | 6.5783e-03 | 4.0532e-03 | 3.3864e-03
$u$ | 7.4464e-02 | 2.7493e-02 | 1.5532e-02 | 1.1842e-02
$v$ | 6.9054e-02 | 2.8942e-02 | 1.6578e-02 | 1.3864e-02
$p$ | 2.2646e-02 | 5.5783e-03 | 4.3532e-03 | 3.0864e-03
Table 1: Expansion wave problem: Relative $L_{2}$ errors in all the primitive
variables for all cases. For XPINNs, we present the total relative error from
all the subdomains.

Figure 9: Expansion wave problem: Loss function variation (left) and the
variation of dynamic weights (right) using PINNs.

Figure 10: Expansion wave problem: Loss function variation in two subdomains
(left) and the combined entropy loss function for both subdomains (right)
using XPINNs.

Next, we show the theoretical generalization bound in XPINN and PINN models
given in [23].

Method | Rel. $L_{2}$ error in $\rho$ | Norms | Theoretical Bounds
---|---|---|---
PINN using AA | 6.5783e-03 | 100.00% | 100.00%
XPINN-1 | 3.3864e-03 | 125.31% | 83.23%
XPINN-2 | 32.59%
Table 2: The Rademacher complexity and the theoretical bounds for the
Expansion wave problem. These are computed using the relative $L_{2}$ error in
density.

These results are calculated using the relative $L_{2}$ error in density. We
present the results including the relative $L_{2}$ error in density, as well
as the calculated theoretical generalization bound in XPINN and PINN models in
Table 2. The results show that XPINN generalizes better than PINN (see the
theoretical bound column). In addition, the Rademacher complexity is presented
in the column ‘Norms’, where that of PINN is set 100% for comparison. The
complexity of PINN is the largest in subdomain 1, since PINN has to fit the
entire target function into one network, which is very complex, whereas in
subdomain 2, the XPINN complexity is slightly larger due to the reduced
training data set.

###  4.2 Oblique shock problem

We now consider another type of attached shock wave problem, i.e., the oblique
shock problem (see regions ① and ④ of the left plot of Fig. 1). Let
$\gamma=1.4$, consider the oblique shock problem with the exact solution given
by

$(\rho,u,v,p)=\left\\{{\begin{array}[]{cc}(0.06688~{}kg/m^{3},738.2~{}m/s,0,~{}9485~{}Pa),~{}~{}~{}~{}~{}~{}~{}\text{
pre-shock,}\\\
(0.09515~{}kg/m^{3},635.9\sin\theta~{}m/s,635.9\cos\theta~{}m/s,1.5\times
10^{4}~{}Pa),\text{ post-shock,}\\\ \end{array}}\right.$ (11)

where $\theta=10^{\circ}$. The loss function of PINN is given as follows:

$\displaystyle\mathcal{J}=$
$\displaystyle~{}\omega_{1}~{}\text{MSE}_{\mathcal{F}}+\omega_{2}~{}\text{MSE}_{\nabla\rho|_{D}}+\omega_{3}~{}\text{MSE}_{Inflow}+\omega_{4}~{}\text{MSE}_{p^{*}}$
(12)
$\displaystyle+\omega_{5}~{}(\text{MSE}_{Mass}+\text{MSE}_{Momentum}+\text{MSE}_{Energy}),$

where $\text{MSE}_{\mathcal{F}}$ corresponds to the Euler equation and the
entropy condition, $\text{MSE}_{\nabla\rho|_{D}}$ corresponds to the density
gradient and $D\subset\Omega:=[0,1]^{2}$ is a region shown in the Fig. 11
(left) (domain with red $+$ points). Note that due to discontinuous solution
we employed finite different formula to enforced the density gradient; see [3]
for more details. $\text{MSE}_{Inflow}$ corresponds to the inflow boundary
condition along the inflow boundary points shown by magenta circles in Fig. 11
(left), $\text{MSE}_{p^{*}}$ corresponds to the pressure values at the point
$(0.4,0)$, and $\text{MSE}_{Mass}$, $\text{MSE}_{Momentum}$ and
$\text{MSE}_{Energy}$ correspond to the total mass, the total momentum and the
total energy, respectively. For PINNs, we use 7 hidden-layers with 30 neurons
in each layer for all cases. In the case of XPINN, the domain is divided into
3 subdomains, and a separate neural network with same network architecture as
PINNs is employed in each subdomain. The Fig. 11 (right) shows the schematic
of three XPINN subdomains along with common interfaces. For stitching these
subdomains, the conservative flux continuity condition [6] is also used along
with residual continuity condition. Also, as discussed in [6, 5], an average
solution is enforced along the common interface for faster convergence. The
number of evenly distributed interface points on both interfaces is 200.

Figure 11: Oblique shock problem: (Left) distributions of the 3900 residual
points (blue $\ast$ points), the 1200 data points for the density gradient
(red $+$ points), and the 120 data points for the inflow conditions (magenta
$\circ$ points). (Right) schematic presentation of subdomains and interfaces
for the XPINN method.

The results are shown in Fig. 12 at locations $x=0.3$ and 0.7. The PINN is
used with or without adaptive activations and dynamic weights, whereas XPINN
is employed with adaptive activations. We note that the vanilla PINN solutions
gives large errors. This implies that the vanilla PINNs cannot deal well with
the problems with discontinuous solutions. On the other hand, it can be
observed that the predicted density, velocity and pressure match well with the
exact solution using adaptive activations for both PINNs and XPINNs as well as
the case of using dynamical weights for PINNs. The relative $L_{2}$ errors in
all the primitive variabls for all the cases are given in Table 3. In the case
of XPINNs, the total relative error in all the subdomains is given together.
This means that the use of adaptive activations or dynamical weights
significantly improves the accuracy of the predictions. We also compare the
results of XPINN with and without entropy condition as shown in Fig. 13. In
both cases, the oblique shock wave is captured exactly.

Figure 12: Oblique shock problem: Primitive variables predicted by PINNs and
XPINNs at $x=0.3$ (top row) and $x=0.7$ (bottom row).

Figure 13: Oblique shock problem: Primitive variables predicted by XPINNs with and without entropy condition at $x=0.3$ (top row) and $x=0.7$ (bottom row). | PINNs without AA/DW | PINNs with AA | PINNs with DW | XPINNs with AA
---|---|---|---|---
$\rho$ | 3.5373e-02 | 6.5533e-03 | 7.1642e-03 | 5.3864e-03
$u$ | 9.9734e-02 | 7.7493e-03 | 1.8442e-02 | 8.1842e-03
$v$ | 8.1637e-02 | 7.8424e-03 | 9.4577e-03 | 7.3864e-03
$p$ | 6.8367e-02 | 4.5363e-03 | 6.7382e-03 | 4.0864e-03
Table 3: Oblique shock problem: Relative $L_{2}$ errors in all the primitive
variables for all cases. In the case of XPINNs, the total relative error from
all the subdomains is given together.

Figure 14: Oblique shock problem: Loss versus iterations for PINNs and XPINNs
given by the left and right figures, respectively. Middle figure shows the
variation of PINN dynamic weights ($\omega_{2}$\- $\omega_{5}$) with the
number of iterations.

Fig. 14 shows the loss function variations with resect to iterations for PINNs
and XPINNs given by the left and right figures, respectively. The middle
figure shows the variation of the PINN dynamic weights ($\omega_{2}$\-
$\omega_{5}$) with the number of iterations.

Now, we show the theoretical generalization bound in XPINN and PINN models for
the oblique shock wave problem.

Method | Rel. $L_{2}$ error in $\rho$ | Norms | Theoretical Bounds
---|---|---|---
PINN using AA | 6.5533e-03 | 100.00% | 100.00%
XPINN-1 | 5.3864e-03 | 41.34% | 89.43%
XPINN-2 | 121.69%
XPINN-3 | 49.72%
Table 4: The Rademacher complexity and the theoretical bounds for the oblique
shock wave problems. These are computed using the relative $L_{2}$ error in
density.

Table 4 shows the results including relative $L_{2}$ error in density and the
calculated theoretical generalization bound in XPINN (for subdomains 1,2, and
3) and PINN models. Again, XPINN generalizes better than PINN (see the
theoretical bound). After observing the Rademacher complexity given in the
column ‘Norms’, the complexity of PINN is higher compared to XPINNs in
subdomains 1 and 3, since PINN has to fit the entire discontinuous target
function into one network. In subdomain 2, the complexity of XPINN is slightly
higher than PINN due to the presence of oblique shock as well as the reduced
training data set.

###  4.3 Detached bow shock problem

We now consider a detached shock problem, i.e., the bow shock problem; see
Fig. 1 (b). Again, we would like to solve the inverse problem of the bow shock
problem by using PINNs and XPINNs. We now consider the bow shock problems with
the following inlet flow conditions

$M_{\infty}=4,\;p_{\infty}=101253.6Pa,\;\rho_{\infty}=1.225kg/m^{3},\;u_{\infty}=1360.6963m/s,\;v_{\infty}=0,\;T_{\infty}=288K.$

The distributions of the residual points, the data points for the density
gradient and the inflow conditions are shown in Fig. 15 (left). The data
points for the pressure are located on the surface of the body. The CFD data
were obtained by using the ‘Openbowshock’ code [46].

Figure 15: Bow shock problem: (Left) distributions of the 5600 residual points
(blue $\ast$ points), the 700 data points for the density gradient (red $+$
points) and the 300 data points for the inflow conditions (magenta $\circ$
points). (Right) schematic presentation of subdomains and curved interface for
the XPINN method.

In addition to the loss function (12), we also enforce the condition (10)
along the surface of the cylinder, namely, we have the following loss function

$\displaystyle\mathcal{J}=$
$\displaystyle~{}\omega_{1}~{}\text{MSE}_{\mathcal{F}}+\omega_{2}~{}\text{MSE}_{\nabla\rho|_{D}}+\omega_{3}~{}\text{MSE}_{Inflow}+\omega_{4}~{}\text{MSE}_{p^{*}}$
(13)
$\displaystyle+\omega_{5}~{}(\text{MSE}_{Mass}+\text{MSE}_{Momentum}+\text{MSE}_{Energy})+\omega_{6}~{}\text{MSE}_{\tilde{n}\cdot{\bm{u}}},$

where
$\text{MSE}_{\tilde{n}\cdot{\bm{u}}}=\frac{1}{N_{\bm{u}}}\sum_{k=1}^{N_{\bm{u}}}(\tilde{n}\cdot{\bm{u}})^{2}(x_{k},y_{k})$
and $(x_{k},y_{k}),\;k=1,\ldots,N_{\bm{u}}$ are points located on the surface
of the body. Moreover, according to previous investigation, PINNs with dynamic
weights and adaptive activation functions usually perform better than the
vanilla PINNs. Therefore, in this case, we use PINNs with dynamic weights.
Further, we used 5 hidden-layers with 160 neurons in each layer. We also
employ XPINN and compare the results with the PINN results. To this end, we
divide the domain into two subdomains as shown in figure 15 (right) keeping
the network architecture same in these two subdomains. The number of points on
the interface is 200, which are uniformly distributed along the interface
line.

Figure 16: Bow shock problem: Reference solutions for all primitive variables.

Figure 17: Bow shock problem: Point-wise errors for PINNs (top row), and
XPINNs (bottom row). The point-wise error in the density, pressure and
velocity components are shown from left to right columns, respectively.
Figure 18: Bow shock problem: Variation of dynamic weights for PINNs.

Figure 19: Bow shock problem: Loss function variation for XPINNs in two
subdomains (left) and the combined entropy loss in the entire computational
domain (right).

We present in Fig. 16 the reference solutions for all primitive variables, and
in Figure 17 the point-wise errors. The point-wise error in the density,
pressure and velocity components are shown from left to right columns. It can
be observed that the predictive accuracy of XPINN is better than that of
PINNs. The dynamic weights for PINNs are shown in Fig. 18. The dynamic weights
$\omega_{4}$ and $\omega_{5}$ converge around $10^{3}$, whereas the weights
$\omega_{2}$ and $\omega_{3}$ converge around $10^{1}$ and $10^{2}$. The loss
functions versus number of iterations are plotted in Figure 19 (left), whereas
the right side figure shows the combined entropy loss in the entire domain.

##  5 Summary

In this work we solved the inverse supersonic flow problems involving
expansion and compression (shock) waves. Accurate solutions to such inverse
problems are often required for designing aerospace vehicles, which are
otherwise difficult or even impossible to solve using the traditional
numerical solvers. For this purpose, we employ physics-informed neural
networks (PINNs) as well as it’s extended space-time domain decomposition-
based version, named as extended PINNs (XPINNs). For the computation, we used
the density gradient data obtained from the Schlieren technique, inflow
boundary data and the pressure wall boundary data. Apart from the governing
compressible Euler equations, the entropy conditions satisfied by the more
relevant viscosity solutions are also enforced. Moreover, we also enforced the
positivity conditions on density and pressure. We compared the results of
PINNs and XPINNs in terms of predictive accuracy. We used the theoretical
generalization bounds in XPINN and PINN models, which show the complexity of
PINN and XPINNs methods as well as how well these methods generalize,
particularly for the expansion wave and the oblique shock wave test cases. We
employed adaptive activation functions for both PINNs and XPINNs, which showed
good predictive accuracy compared to fixed activation functions. Moreover, the
use of dynamic weights shows accuracy improvements over fixed weights. For the
problems with discontinuities, additional techniques, such as adaptive
activate functions and dynamical weights, are important and can significantly
improve the accuracy of the PINN solutions. XPINNs with locally tuned neural
networks also gives more accurate results suggesting that XPINN is also a
promising and scalable method for tackling problems with complex solutions
involving local discontinuous features such as shocks.

##  Appendix A Pedagogical example: Inverse problem for the two-dimensional
Euler equations with smooth solutions

(a) Data for $\nabla\rho$

(b) Exact solution: density

(c) PINN solution: density

(d) PINN solution: pressure

(e) PINN solution: velocity $u$

(f) PINN solution: velocity $v$

Figure 20: (a) We only use the data of density gradient along the line $y=x$
shown in the plot. (b) Exact solution of the density. (c) PINN solution of the
density. (d) PINN solution of the pressure. The exact solution of the pressure
is $p(x,y)=1$. (e) PINN solution of the velocity $u$. The exact solution of
the velocity $u$ is $u(x,y)=0.7$. (f) PINN solution of the velocity $v$. The
exact solution of the velocity $v$ is $v(x,y)=0.3$.

_Consider the Euler equation ( 1) in the two-dimensional case with the
periodic boundary conditions and initial conditions_

$U_{0}=(\rho_{0},~{}~{}u_{0},~{}~{}v_{0},~{}~{}p_{0})=\left(1.0+0.2\sin(\pi(x+y)),~{}~{}0.7,~{}~{}0.3,~{}~{}1.0\right),$

_in which case we have the exact solutions_

$(\rho,~{}~{}u,~{}~{}p)=\left(1.0+0.2\sin(\pi(x+y-(u+v)t)),~{}~{}0.7,~{}~{}0.3,~{}~{}1.0\right),$

_where $(x,y)\in(-1,1)^{2}$. _

In the one-dimensional case, we choose $D$ to be the whole computational
domain. However, in the two-dimensional case, we use density gradient only in
a subdomain. In this example, we use the density gradient along the line
$y=x$, see Figure 20(a). For the pressure, we only use one single point value
of the pressure at different time, here we use $p(-1,-1,t)$. We also use the
inflow conditions. Therefore, the loss function is given by:

$\mathcal{J}=\text{MSE}_{\mathcal{F}}+\text{MSE}_{\nabla\rho|_{y=x}}+\text{MSE}_{p^{*}}+\text{MSE}_{Inflow}.$
(14)

The results are shown in Figs. 20 (c)-(f). Observe that we now have good
predictions for all the states, i.e., the density, velocity and the pressure.

## Acknowledgements

This work was supported by the OSD/AFOSR MURI grant FA9550- 20-1-0358 and the
Alexander von Humboldt fellowship to G.E. Karniadakis. The work of Z. Mao was
conducted while he was a postdoc at Brown University and also during his one
month-visit to the Technical University of Munich working with Prof. N. Adams.

## References

  * [1] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686–707.
  * [2] M. Raissi, Z. Wang, M. S. Triantafyllou, G. E. Karniadakis, Deep learning of vortex induced vibrations, arXiv preprint arXiv:1808.08952.
  * [3] 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.
  * [4] M. Raissi, A. Yazdani, G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367 (6481) (2020) 1026–1030.
  * [5] A. D. Jagtap, G. E. Karniadakis, Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations, Communications in Computational Physics 28 (5) (2020) 2002–2041.
  * [6] 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.
  * [7] K. Shukla, A. D. Jagtap, J. L. Blackshire, D. Sparkman, G. E. Karniadakis, A physics-informed neural network for quantifying the microstructural properties of polycrystalline nickel using ultrasound data: A promising approach for solving inverse problems, IEEE Signal Processing Magazine 39 (1) (2022) 68–77.
  * [8] K. Shukla, A. D. Jagtap, G. E. Karniadakis, Parallel physics-informed neural networks via domain decomposition, Journal of Computational Physics 447 (2021) 110683.
  * [9] Y. Chen, L. Lu, G. E. Karniadakis, L. Dal Negro, Physics-informed neural networks for inverse problems in nano-optics and metamaterials, Optics express 28 (8) (2020) 11618–11633.
  * [10] U. bin Waheed, E. Haghighat, T. Alkhalifah, C. Song, Q. Hao, Pinneik: Eikonal solution using physics-informed neural networks, Computers & Geosciences (2021) 104833.
  * [11] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes, A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics, Computer Methods in Applied Mechanics and Engineering 379 (2021) 113741.
  * [12] Y. Yang, P. Perdikaris, Adversarial uncertainty quantification in physics-informed neural networks, Journal of Computational Physics 394 (2019) 136–152.
  * [13] S. Cai, Z. Wang, F. Fuest, Y. J. Jeon, C. Gray, G. E. Karniadakis, Flow over an espresso cup: inferring 3-d velocity and pressure fields from tomographic background oriented schlieren via physics-informed neural networks, Journal of Fluid Mechanics 915.
  * [14] E. Kharazmi, Z. Zhang, G. E. Karniadakis, hp-vpinns: Variational physics-informed neural networks with domain decomposition, Computer Methods in Applied Mechanics and Engineering 374 (2021) 113547.
  * [15] A. D. Jagtap, D. Mitsotakis, G. E. Karniadakis, Deep learning of inverse water waves problems using multi-fidelity data: Application to serre-green-naghdi equations, arXiv preprint arXiv:2202.02899 (2022).
  * [16] Y. Shin, J. Darbon, G. E. Karniadakis, On the convergence and generalization of physics informed neural networks, arXiv e-prints (2020) arXiv–2004.
  * [17] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes, IMA Journal of Numerical Analysis.
  * [18] J. D. Anderson, Computational fluid dynamics: the basics with applications, International edition: McGraw-Hill.
  * [19] A. Bayliss, E. Turkel, Far field boundary conditions for compressible flows, Journal of Computational Physics 48 (2) (1982) 182–199.
  * [20] D. Hall, S. Yen, Implementation of nonreflective boundary condition at the outflow boundary, Numerical Boundary Condition Procedures (1981) 45.
  * [21] 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.
  * [22] O. Fuks, H. A. Tchelepi, Limitations of physics informed machine learning for nonlinear two-phase transport in porous media, Journal of Machine Learning for Modeling and Computing 1 (1).
  * [23] Z. Hu, A. D. Jagtap, G. E. Karniadakis, K. Kawaguchi, When do extended physics-informed neural networks (xpinns) improve generalization?, arXiv preprint arXiv:2109.09444 (2021).
  * [24] J. Magiera, D. Ray, J. S. Hesthaven, C. Rohde, Constraint-aware neural networks for riemann problems, Journal of Computational Physics 409 (2020) 109345.
  * [25] D. A. Bezgin, S. J. Schmidt, N. A. Adams, A data-driven physics-informed finite-volume scheme for nonclassical undercompressive shocks, Journal of Computational Physics 437 (2021) 110324.
  * [26] Z. Mao, L. Lu, O. Marxen, T. A. Zaki, G. E. Karniadakis, Deepm&mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators, Journal of Computational Physics 447 (2021) 110698.
  * [27] M. Monfort, T. Luciani, J. Komperda, B. Ziebart, F. Mashayek, G. E. Marai, A deep learning approach to identifying shock locations in turbulent combustion tensor fields, in: Modeling, Analysis, and Visualization of Anisotropy, Springer, 2017, pp. 375–392.
  * [28] Z. He, E. Guo, X. Huang, C. Mo, W. Kang, F. Zhang, C. Wang, H. Zhang, X. Chu, G. Jia, et al., Inverting shock-wave temperatures via artificial neural networks, Journal of Applied Physics 127 (12) (2020) 125901.
  * [29] R. Courant, K. O. Friedrichs, Supersonic flow and shock waves, Vol. 21, Springer Science & Business Media, 1999.
  * [30] H. W. Liepmann, A. Roshko, Elements of gasdynamics, Courier Corporation, 2001.
  * [31] R. D. Zucker, O. Biblarz, Fundamentals of gas dynamics, John Wiley & Sons, 2002\.
  * [32] C. M. Dafermos, Hyperbolic conservation laws in continuum physics, 4th Edition, Vol. 325 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, Berlin, 2016.
  * [33] S. Bianchini, A. Bressan, Vanishing viscosity solutions of nonlinear hyperbolic systems, Annals of mathematics (2005) 223–342.
  * [34] E. Godlewski, P.-A. Raviart, Numerical approximation of hyperbolic systems of conservation laws, Vol. 118, Springer, 1996.
  * [35] A. D. Jagtap, K. Kawaguchi, G. E. Karniadakis, Adaptive activation functions accelerate convergence in deep and physics-informed neural networks, Journal of Computational Physics 404 (2020) 109136.
  * [36] A. D. Jagtap, K. Kawaguchi, G. E. Em Karniadakis, Locally adaptive activation functions with slope recovery for deep and physics-informed neural networks, Proceedings of the Royal Society A 476 (2239) (2020) 20200334.
  * [37] A. D. Jagtap, Y. Shin, K. Kawaguchi, G. E. Karniadakis, Deep kronecker neural networks: A general framework for neural networks with adaptive activation functions, Neurocomputing 468 (2022) 165–180.
  * [38] X. Glorot, Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, 2010, pp. 249–256.
  * [39] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of machine learning research 18.
  * [40] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., Tensorflow: A system for large-scale machine learning, in: 12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16), 2016, pp. 265–283.
  * [41] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, A. Lerer, Automatic differentiation in pytorch.
  * [42] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980.
  * [43] R. H. Byrd, P. Lu, J. Nocedal, C. Zhu, A limited memory algorithm for bound constrained optimization, SIAM Journal on scientific computing 16 (5) (1995) 1190–1208.
  * [44] S. Wang, Y. Teng, P. Perdikaris, Understanding and mitigating gradient flow pathologies in physics-informed neural networks, SIAM Journal on Scientific Computing 43 (5) (2021) A3055–A3081.
  * [45] A. D. Jagtap, G. E. Karniadakis, Adaptive activation functions accelerate convergence in deep and physics-informed neural networks, arXiv preprint arXiv:1906.01170 (2019).
  * [46] S. Karpuk, Openbowshock 2d (https://www.mathworks.com/matlabcentral/fileexchange/69993-openbowshock-2d), MATLAB Central File Exchange.


# Record statistics of emitted energies $-$ prediction of an upcoming failure

Subhadeep Roy<EMAIL_ADDRESS>Dept. of Physics, Birla
Institute of Technology & Science Pilani, Hyderabad Campus, Secunderabad,
Telangana 500078.

###### Abstract

The article reports a numerical investigation of the breakdown of a disordered
system considering the effect of local stress concentration under the action
of an external tensile force. The statistics of the record-breaking magnitudes
of emitted energies during the failure process as well as the waiting time to
achieve those record events show rich behavior. The latter includes
information about the acceleration and subsequent catastrophic failure through
its non-monotonic behavior. The maximum waiting time is also correlated with
the maximum change in elastic energy as the model evolves, a different way of
predicting an upcoming failure, which is consistent with our hypothesis as
well. At a moderate disorder, such a prediction can be done with higher
accuracy while at a low disorder, due to the abrupt nature of the failure
process our hypothesis does not hold well.

###### keywords:

Disordered systems, Fiber bundle model, Avalanche statistics, Critical stress,
Record statistics.

##  1 Introduction

The fracture process in a disordered system has been explored widely from the
point of view of statistical mechanics due to its practical application as
well as fundamental theoretical problems [1]. One interesting aspect in this
regard is the fact that heterogeneous materials do not break at once when
applied to external mechanical stress.

The micro and mesoscale heterogeneity of materials has the consequence that
their fracture process evolves through local breaking events, known as
intermittent bursts, that produce crackling noise which can be captured
through an acoustic emission (AE) experiment [2]. The crackling noises, when
translated into bursts and subsequently emitted energies, can be exploited to
predict an upcoming catastrophic event [3]. Such acoustic signals (also known
as precursors) become more and more populated closer to the global failure and
approach the catastrophic event in an accelerating manner indicated by the
increasing rate of deformation [4]. During the initiation and propagation of
cracks through a disordered media, complexity arises due to material disorders
such as the distribution of flaws, dislocations, and micro-cracks present. As
the external stress is increased, micro-cracks generate randomly and the
stress ﬁelds of such micro-cracks interact with each other making this process
correlated.

The intensity and density of the precursory activity and its accuracy for an
upcoming failure process depend on the strength of the disorder. In one
limiting case, where the disorder strength is close to zero, the catastrophic
failure occurs suddenly without any prior warning. In this limit, no
precursors are observed [5]. On the other hand, when the disorder strength is
sufficiently high, a gradual accumulation of damage is observed showing bursts
and avalanches with a rate growing with time [6].

Examples of such avalanches are observed on different scales, starting from
laboratory scales like fracturing of wood [7], cellular glass [8], concrete
[4], etc., to a much larger scale of seismic events where this phenomenon is
reminiscent of the Gutenberg-Richter law [9]. Apart from failure, there are
many other statistical models that incorporate avalanches such as magnetic
materials showing the Barkhausen effect [10] and flux lines in superconductors
[11], friction [12], Heines jumps during the flow of fluid through a porous
media [13], and so on. Therefore, understanding the physical mechanisms of
precursor activities and avalanche dynamics goes well beyond the study of
fracture and breakdown.

One robust behavior that is observed during a failure process is how the
emitted energies, captured during an acoustic emission (AE) experiment, are
distributed statistically. In most cases, the distribution of emitted energy
shows a scale-free behavior with a non-universal exponent that depends on the
exact nature of the material. For example, a value 1.41 for human teeth [14],
1.25 in paper [15], 1.3 in synthetic plaster [4], 1.51 in wood and 2.0 in
fiber-glass [7], 1.9 in fiberglass [16] and 1.5 in cellular glass [8] etc.

Apart from experiments, a considerable amount of effort has been given to
large-scale simulations of lattice models [1] for a better understanding of an
evolving disordered system and the subsequent failure associated with the
motion. Such a lattice model can be an electrical analog like a random
resistor network [17], a mechanical analog like a random spring network [18],
or other more complex systems like continuum models for damage mechanics [19].
The lattice models usually evolve by solving a number of coupled equations
(depending on the complexity of the model). A long-range interaction is
introduced in these models in terms of disorder incorporated through random
threshold values or bond dilution. A scaling concept [20] can also be
extracted from the lattice models through their description of topological
properties during the failure process.

A relatively simpler lattice model like the fiber bundle model (FBM) shows a
scale-free avalanche size distribution in the mean-field limit with a
universal exponent of 5/2 [21]. Also, a relatively lower exponent (= 3/2) is
found in the mean-field limit when the model is close to the failure point
[22]. Though, as the FBM deviates from the mean-field limit and the local
stress concentration starts playing a major role in the failure process, the
distribution becomes exponential instead of scale-free [23]. A recent study on
emitted energies during the failure process of an FBM shows that the
distribution of energies is scale-free irrespective of what the distribution
of avalanche is. The energy distribution shows a slope (= 5/2) similar to the
avalanche distribution in a mean-field limit and a higher slope (= 3.5) in
presence of local stress concentration [24]. Moreover, with such local stress
concentration, the model shows different modes [25, 26, 27, 28] of failure as
well as fluctuation in failure time [29].

In this paper, I will lead with a previous work where the record statistics of
the avalanche sizes are explored in the context of the mean-field FBM [30, 31]
and later in presence of heterogeneous stress field [32]. Apart from this,
record statistics have been explored during compressive failure of porous
materials [33] as well as in the context of mining collapse [34]. In these
articles, a warning of a catastrophic event is reported which is not complete
without correlating the failure point with an observable system parameter.
Moreover, I have used emitted energies as the recorded event, and not the
avalanches, since it’s rather realistic as well as convenient to capture
energies through an Acoustic Emission (AE) experiment. A couple of recent
papers [24, 35] have discussed the correlation between avalanche and energy
emitted during the failure of a fiber bundle model and how considering energy
is more advantageous from a statistical point of view especially when the
stress field is heterogeneous in space. These modifications make the present
work much more relevant to real-life scenarios making it available for future
experimental validation.

##  2 Description of Fiber Bundle Model

A conventional fiber bundle model [36, 37] consists of a number of parallel
Hookean fibers attached between two bars (in case of 1d) or plate (in case of
2d) which are pulled apart. If there are $L$ fibers and the externally applied
stress is $F$ then the stress a single fiber experience is $\sigma=F/L$. The
disorder is introduced in the model in terms of fluctuation in local strength
values (say $h$) from fiber to fiber. In the present work, I have considered
the following distribution

$p(h)\sim\begin{cases}h^{-1},&(10^{-\beta}\leq h\leq 10^{\beta})\\\ 0.&({\rm
otherwise})\end{cases}$ (1)

Here $\beta$ is related to the span as well as the disorder strength - the
higher the $\beta$ value higher the disorder strength. Such a power-law
distribution has been observed earlier in local material strength values [38].
If the local stress of a certain fiber exceeds its threshold value, that fiber
breaks irreversible and the stress carried by the broken fiber is
redistributed.

We have adopted here a local load-sharing scheme [39, 40] for the
redistribution where the extra stress due to the rupture is redistributed
among the nearest neighbors only. Since the model we consider is 1d, there are
only two nearest neighbors (the redistribution scheme is elaborated in the
‘Method’ section towards the end of this article). The other extreme is the
mean-field limit where the stress of the broken fiber is redistributed equally
among all surviving fibers. Apart from these two extreme limits, an
intermediate stress redistribution scheme [41, 42] has also been adopted that
allows a scale-free stress relaxation making the amount of stress to be
redistributed a function of distance from the broken fiber. Further rupture
can happen at this point without increasing the applied stress as a result of
the stress redistribution, which makes the local stress values elevate. This
starts an avalanche which might lead to a global failure or to the next steady
state with a fraction of fibers broken and the other fraction intact. In the
latter situation, we increase the applied stress just enough to break the next
weakest fiber only, starting a new avalanche. Such a process goes on until all
fibers are broken. The last stress at which the model attains global failure
is known as the critical stress of the system.

##  3 Numerical Results

Numerically, I have studied a fiber bundle model with system size varying
between $10^{3}$ and $10^{7}$. The disorder is varied by tuning the parameter
$\beta$ from 0.5 to 4.0. A low $\beta$ suggests a very low disorder where the
failure process is observed to be brittle-like abrupt [43, 44, 45]. On the
other hand, for a high disorder, the bundle breaks in a number of avalanches
[25] with energies emitted in every avalanche [24]. This letter deals with the
‘record events’ associated with the energies emitted. We will discuss a unique
way to connect the statistics of the record energy bursts with the critical
value of the applied stress at which the bundle hits the global failure.

Figure 1: (a) Energy emitted during the failure process vs the number of
steady states. The red dots correspond to the record events with magnitude
$r_{k}$ corresponding to the $k$th event. $\tau_{k}$ is the waiting time to
achieve the $k$th record event. The figure shows $r_{k}$ for $k=10$ and
$\tau_{k}$ for $k=14$. (b) Variation of the average number $\langle
N_{r}\rangle$ of record events for a certain number $n$ of total avalanches
for $\beta$ ranging between 0.5 and 4.0. Initially, $\langle
N_{r}\rangle\approx n$ as a certain event has a higher probability to be
greater than the last one. Later the process slows down and speed-up again
closer to the onset of acceleration.

We start by observing the nature of the energy emitted as the model evolves
with increasing applied stress. Figure 1(a) shows this variation of emitted
energy $E$ as a function of the number $n$ of the steady states. Each vertical
line on the diagram represents the magnitude of the energy emitted during the
avalanche corresponding to that particular steady state. As already mentioned
in the model description, a steady state is a point during the evolution of
the model where external stress is required to initiate further ruptures. From
the definition of an ‘avalanche’, $n$ can stand for the number of avalanches
as well. The size $s$ of an avalanche is the number of fibers broken in
between two stress increments (and hence in between two steady states). The
corresponding emitted energy will then be

$\displaystyle
E(s)=\displaystyle\sum_{j=1}^{s}\displaystyle\frac{1}{2}h_{j}^{2}$ (2)

where $j$ runs from $1$ to $s$ (size of the avalanche) and $h_{j}$’s are the
threshold values of the fibers broken during that avalanche. The term
$1/2h^{2}_{j}$ is the damage energy for the rupture of the jth fiber. A recent
article [24] shows the lack of correlation between $s$ and $E$ in presence of
local stress concentration in FBM, suggesting that a high energy can be
achieved through a small burst and vice versa. We have chosen energy over
avalanche for our present study since it is a quantity that can be captured
easily during an acoustic emission (AE) experiment. The red dots in the same
figure denotes the record events: an event that is larger than the last one.
The magnitude $k$ of a record event is given by $r_{k}$ while $\tau_{k}$ is
the waiting time for that event. $\tau_{k}$ can also be defined as the
lifetime of the $(k-1)th$ event. One main aim of the article is to study the
behavior of $\tau_{k}$ and this can be used to predict the failure point.

Figure 1(b) shows the variation of average number of record events $\langle
N_{r}\rangle$ increases with the number $n$ of the total avalanches.
Initially, $\langle N_{r}\rangle$ increases linearly with $n$. This happens as
almost every avalanche emits record energy due to the fact that the reference
for record events is very low at the beginning. Later $\langle N_{r}\rangle$
slows down and increases again close to the failure point with a faster rate
but slower than the initial linear increment. This final increased rate of
$\langle N_{r}\rangle$ suggests the onset of acceleration when the model is
closer to global failure. Inset of figure 1(b) shows the same behavior but for
a single configuration. The vertical jumps take place whenever there is a new
record event and the horizontal lines corresponds to no activity zone and
hence contributes to the waiting time. The figure explicitly shows that there
are a lot of activities with a number of vertical jumps in the beginning as
well as at the end while somewhere in the middle, mostly horizontal lines are
observed leading to a peak in the waiting time. Such a slowdown can also be
observed during a seismic event, a phenomenon widely known as the seismic
quiescence. The acceleration after the quiescence and just before the main
shock has been measured quantitatively in the literature by the decrease in
slope in the distribution of fore-shock magnitudes and famously known as the
b-value decrease [46]. Such sequential record events can be characterized as
IID sequences which universally shows a logarithmic increment [47]. Such an
increment has also been observed in the context of fiber bundle model recently
[30]. In this present article a scale-free increase is observed rather a
logarithmic one, most probably due to the lack of correlation between
avalanches and energies in the presence of local stress concentration [24,
35].

Figure 2: (a) Average (over $10^{4}$ realizations) waiting time
$\langle\tau_{k}\rangle$ as a function of the order $k$ of the record event.
$\langle\tau_{k}\rangle$ increases initially with $k$, and decreases again
prior to the failure point. $\langle\tau_{k}\rangle$ reaches its maximum at
$k=k^{\ast}$. (b) Average critical stress $\langle\sigma_{c}\rangle$ plotted
against $\langle\sigma^{\ast}\rangle(=\langle\sigma\rangle\ \text{at}\
k=k^{\ast})$ for $L=10^{5}$. The color gradient shows the values of $\beta$.
At low $\beta$, $\langle\sigma_{c}\rangle=\langle\sigma^{\ast}\rangle$,
suggesting no scope of prediction. For moderate $\beta$, a prediction can be
done since we observe
$\langle\sigma_{c}\rangle=a\langle\sigma^{\ast}\rangle-b$, where $a=5.3\pm
0.18$, and $b=1.05\pm 0.81$. Inset (I) and (II) shows how $k_{m}$, the maximum
number of record events registered, is related to disorder strength $\beta$
and critical stress $\sigma_{c}$ respectively.

Next, we turn to study the average waiting time $\langle\tau_{k}\rangle$ as a
function of the index $k$ of record events. The averaging is done for each $k$
value over $10^{4}$ realizations. Figure 2(a) shows $\langle\tau_{k}\rangle$
as a function of $k$ for $\beta=0.5$, 1.0, 1.5 and 2.0.
$\langle\tau_{k}\rangle$ shows a non-monotonic behavior with a peak at a
certain value of $k$, say $k^{\ast}$. On both sides of $k^{\ast}$, $\tau_{k}$
decreases suggesting that we will see closely spaced record events as we move
away from $k^{\ast}$. For $k<k^{\ast}$, this happens since the reference frame
that decides an event to be the record one is itself small. On the other side
($k>k^{\ast}$), though the reference is quite large, a decrease in $\tau_{k}$
suggests we get many bigger events within that limit. This $k^{\ast}$, as we
will see later in this paper, can be connected to the failure point predicting
an upcoming catastrophic event.

Figure 3: (a) Change in average elastic energy $\Delta E_{e}$ as a function of
the external stress for $\beta=0.8$, 1.0 and 1.2. $\Delta E_{e}$ shows a non-
monotonic behavior with a peak at $\sigma_{m}$. For a certain $\beta$, such
peak takes place before $\langle\sigma_{c}\rangle$ and after
$\langle\sigma^{\ast}\rangle$, as it is shown in the inset. (b) The figure
shows how $\langle\sigma_{m}\rangle$ varies with $\langle\sigma^{\ast}\rangle$
and $\langle\sigma_{c}\rangle$ as $\beta$ is continuously varied. We observe:
$\langle\sigma_{m}\rangle=p\langle\sigma^{\ast}\rangle-q$ and
$\langle\sigma_{c}\rangle=c\langle\sigma_{m}\rangle-d$, where $p=4.1\pm 0.09$,
$q=0.75\pm 0.04$, $c=1.3\pm 0.04$, and $d=0.1\pm 0.05$.

Figure 2(b) shows the variation in $\langle\sigma_{c}\rangle$ and
$\langle\sigma^{\ast}\rangle$, where $\langle\rangle$ denotes an average over
$10^{4}$ realizations. $\langle\sigma_{c}\rangle$ is the average critical
stress which is applied externally on the bundle just before the global
failure. $\langle\sigma^{\ast}\rangle$ is the stress at which $\tau_{k}$
reaches its maximum value. So,
$\langle\sigma^{\ast}\rangle=\langle\sigma\rangle_{k=k^{\ast}}$.

We study the correlation between $\langle\sigma^{\ast}\rangle$ and
$\langle\sigma_{c}\rangle$. Each point, within the figure, is for a certain
value of $\beta$ ranging between 0 and 1.6 (see the color gradient). We see
two branches of data. The first branch corresponds to low disorder strength
when the failure process is brittle-like abrupt [44, 45]. Here,
$\langle\sigma^{\ast}\rangle=\langle\sigma_{c}\rangle$, which means the peak
in $\langle\tau_{k}\rangle$ takes place at the same time the model reaches the
global failure. A prediction will not be possible in this case. However, for
high disorder, we observe $\langle\sigma_{c}\rangle$ to increase linearly with
$\langle\sigma^{\ast}\rangle$:

$\displaystyle\langle\sigma_{c}\rangle=a\langle\sigma^{\ast}\rangle-b$ (3)

where $a=5.3\pm 0.18$, and $b=1.05\pm 0.81$. This means,
$\langle\sigma_{c}\rangle$ has a value much higher than
$\langle\sigma^{\ast}\rangle$ and we can have information about the failure
point much before the failure occurs. The inset of figure 2(b) shows how
$k_{m}$, the maximum number of record events registered varies with disorder
strength $\beta$ and the critical stress $\sigma_{c}$ respectively. A higher
$k_{m}$ suggests the bundle is relatively stable as the final acceleration
starts at a higher time. This agrees with how $k_{m}$ behaves with
$\sigma_{c}$ \- a higher critical stress when more record events are captured,
supporting the higher stability at higher disorder strength.

We also want to take this opportunity to explore a different technique to
predict the critical stress which was suggested recently [48] (but in the
context of the mean-field limit) and see whether it is consistent with our
hypothesis. For this, we will study the elastic energy associated with system
evolution. The elastic energy at a certain applied stress is defined as:

$\displaystyle
E_{e}=\displaystyle\sum_{j}\displaystyle\frac{1}{2}[\sigma(j)]^{2}$ (4)

where $j$ runs over all intact fibers and $\sigma(j)$ is the local stress of
fiber $j$. $E_{e}$ will then be a non-monotonic function since there are many
fibers in the beginning that are less strained while at the end there is a
very less number of fibers that are highly strained. Unfortunately, such a
peak in $E_{e}$ takes place after $\langle\sigma_{c}\rangle$ when the model
has already become unstable [48]. Instead, if we calculate the change in
elastic energy $\Delta E_{e}$ as we go from one steady state to the next one,
we observe a maximum value of $\Delta E_{e}$ at $\langle\sigma_{m}\rangle$,
which takes place before $\langle\sigma_{c}\rangle$.

Figure 4: (a) $\tau_{k}$ as a function of $k$ for a 10 different realizations
for a bundle of size $10^{7}$. $\tau_{k}$ shows a peak in all the cases making
our hypothesis valid not only for average parameters but for a single
realization as well. (b) Prediction of $\sigma_{c}$ from $\sigma^{\ast}$ for
100 different configurations for $\beta=1.0$. We observe:
$\sigma_{c}=a^{\prime}\sigma^{\ast}-b^{\prime}$ from least square fit, where
$a^{\prime}=6.18\pm 1.37$ and $b^{\prime}=1.96\pm 0.94$.

Figure 3(a) shows $\Delta E_{e}$ to increase as $\langle\sigma\rangle^{2}$
prior to the peak independent of $\beta$ and decreases sharply afterward. The
inset of the same figure shows schematically where
$\langle\sigma^{\ast}\rangle$, $\langle\sigma_{m}\rangle$ and
$\langle\sigma_{c}\rangle$ lie. The difference between
$\langle\sigma_{m}\rangle$ and $\langle\sigma_{c}\rangle$ is very less while
$\langle\sigma^{\ast}\rangle$ takes place much early (see the inset). So,
$\langle\sigma^{\ast}\rangle$ serves as a better prediction point rather than
$\langle\sigma_{m}\rangle$ but we can treat it to be the second warning point
prior to global failure. Figure 3(b) shows the following correlation

$\displaystyle\langle\sigma_{m}\rangle=p\langle\sigma^{\ast}\rangle-q\ \
\text{and}\ \ \langle\sigma_{c}\rangle=c\langle\sigma_{m}\rangle-d$ (5)

where where $p=4.1\pm 0.09$, $q=0.75\pm 0.04$, $c=1.3\pm 0.04$, and $d=0.1\pm
0.05$. Above correlation shows that not only $\langle\sigma_{c}\rangle$ but
$\langle\sigma_{m}\rangle$ also can be predicted if we know
$\langle\sigma^{\ast}\rangle$.

Finally, we will discuss our hypothesis in the context of a practical
application. So far, we have only discussed the average quantities which are
not useful when comparing with experiments since a single experiment is
equivalent to a single realization. Here, we chose a $\beta$ value large
enough ($=1.0$) so that the prediction is possible. Figure 4(a) shows the
variation in $\tau_{k}$ with index $k$ for 10 different realizations for
$L=10^{7}$ and $\beta=1.0$. For each realization, $\tau_{k}$ consistently
shows a peak before the global failure and we can extract both $\sigma^{\ast}$
and $\sigma_{c}$ for each realization. Figure 4(b) shows how $\sigma_{c}$
varies with $\sigma^{\ast}$ for 100 such realizations. The dotted line shows
the best fit according to the least square method:

$\displaystyle\sigma_{c}=a^{\prime}\sigma^{\ast}-b^{\prime}$ (6)

where where $a^{\prime}=6.18\pm 1.37$ and $b^{\prime}=1.96\pm 0.94$. The above
correlation shows that the prediction of critical stress $\sigma_{c}$ is
possible for individual configurations as well.

## Discussion

When a disordered system is acted upon a constant stress (creep [49, 50]) or a
slowly varying external load, the process leads to a jerky evolution of the
fracture process which is often detected externally by a sequence of acoustic
signals. The present work is mainly based on such a sequence of acoustic
signals where the record events (most prominent ones) are the center of
interest. The paper explores an efficient method through studying the
statistics of the record signals to predict the onset of an accelerating
nature and eventual catastrophic failure.

As an important outcome of the work, I have represented that during quasi-
brittle fracture, by combining the record statistics with avalanche size
distribution during a failure process it is possible to predict the failure
point with high accuracy. The technique is built on the behavior of waiting
times between consecutive record bursts. Such waiting-time shows a non-
monotonic behavior where the value of applied stress at which the waiting time
is maximum shows a high correlation with critical stress. Such a technique has
a two-fold purpose: first, a warning of an upcoming failure through the
maximum value of waiting time as well as the maximum change in elastic energy,
and second, correlating the critical point with the stress/strain at this
maxima. The average applied stress $\langle\sigma^{\ast}\rangle$ at maximum
waiting time, as well as the average applied stress $\langle\sigma_{m}\rangle$
at the maximum change in elastic energy, are observed to be correlated
linearly with the average critical stress $\langle\sigma_{c}\rangle$. One
challenge is such an early warning of imminent failure requires to be done
much before the model enters the unstable state. Since
$\langle\sigma^{\ast}\rangle$ takes place much before
$\langle\sigma_{m}\rangle$, the process with the record statistics gives us a
larger time for damage control than the elastic energy treatment. As it was
already observed in the literature, a disordered system like the fiber bundle
model shows brittle-like abrupt failure [43, 44, 45, 25] (in a single
avalanche) when the disorder strength is low. Clearly, this hypothesis fails
in this limit of low disorder, where the failure process is brittle-like
abrupt [43, 44, 45, 25], due to a lack of avalanche statistics. Particularly,
it was shown that at low disorder $\langle\sigma^{\ast}\rangle$ is almost
equal to $\langle\sigma_{c}\rangle$, and a prediction is not possible. That’s
why this treatment is only fruitful in the quasi-brittle region only.

This paper is a step closer toward failure prediction from a statistical point
of view. It is imperative to extend the study to realistic systems and
engineering materials towards understanding the validity of this novel process
starting from laboratory scales like glassy materials, and damage mechanics in
continuum media to geological scales like Heines jumps during multi-phase
flows in water and oil industries as well as large scale seismic events.

## Methods

The methods used for the work are numerical. Codes written by the author in
C-language are used for the simulations. The fibers are broken following the
weakest-link-of-the-chain theory. If $h_{i}$’s are the threshold values of the
fibers and $\sigma_{i}$ are corresponding local stress values, then the fiber
that breaks has the minimum value of $(h_{i}-\sigma_{i})$. The stress of the
broken fiber is redistributed among the neighboring fibers. In one dimension,
which I have used for this study, the two neighbors (right and left) have the
following local stress:
$\sigma_{r(l)}\rightarrow\sigma_{r(l)}+\displaystyle\frac{\sigma_{b}d_{l(r)}}{d_{r}+d_{l}}$,
where $\sigma_{b}$ is the stress of the broken fiber, $d_{l}$ and $d_{r}$ are
the distance of the neighboring fiber from the broken one on the left and
right sides respectively. After the redistribution, some other fibers may
break due to the elevated local stress profile starting an avalanche. After an
avalanche, the next weakest fiber is broken by calculating the minimum of
$(h_{i}-\sigma_{i})$ again. For sufficiently nice statistics, I have simulated
a bundle of size ranging between $10^{3}$ and $10^{7}$ with $10^{4}$
realization.

## Acknowledgements

SR acknowledges Birla Institute of Technology & Science, Hyderabad campus for
the support during this work. A special thank goes to Prof. Rajesh Ravindran,
Prof. Anuradha Banerjee, Prof. Purusattam Ray, and Subrat Senapati for
discussions and valuable comments from time to time.

## Author contributions statement

The author confirms being the sole contributor of this work and has approved
it for publication.

## Accession codes

The codes can be accessed through the link
https://github.com/sroy2807/failurePredictionFBM.git.

## Competing interests

The author declares no competing financial interests.

## References

  * [1] H. J. Herrmann and S. Roux, Statistical Models for the Fracture of Disordered Media, North-Holland, Amsterdam (1990).
  * [2] J. Baro, A. Corral, X. Illa, A. Planes, E. K. H. Salje, W. Schranz, D. E. Soto-Parra, and E. Vives, Phys. Rev. Lett. 110, 088702 (2013).
  * [3] Guillaume F Nataf et. al, J. Physics: Condens. Matter 26, 275401 (2014). Vasseur, J. et. al, Sci. Rep. 5, 13259 (2015).
  * [4] A. Petri, G. Paparo, A. Vespignani, A. Alippi, and M. Costantini, Phys. Rev. Lett. 73, 3423 (1994).
  * [5] Zapperi, S., Ray, P., Stanley, H. E. and Vespignani, Phys. Rev. Lett. 78, 1408 (1997).
  * [6] Ramos, O., Cortet, P.-P., Ciliberto, S. and Vanel, Phys. Rev. Lett. 110, 165506 (2013).
  * [7] A. Garcimartin, A. Guarino, L. Bellon, and S. Ciliberto, Phys. Rev. Lett. 79, 3202 (1997);
  * [8] C. Maes, A. Van Moffaert, H. Frederix, and H. Strauven, Phys. Rev. B 57, 4987 (1998).
  * [9] B. Gutenberg and C. F. Richter, Bull. Seismol. Soc. Am. 34, 185 (1994).
  * [10] S. Zapperi, P. Cizeau, G. Durin, and H. E. Stanley, Phys. Rev. B 58, 6353 (1998).
  * [11] S. Field, J. Witt, F. Nori, and X. Ling, Phys. Rev. Lett. 74, 1206 (1995).
  * [12] S. Ciliberto and C. Laroche, J. Phys. I 4, 223 (1994).
  * [13] A. H. Thompson, A. J. Katz, and R. A. Raschke, Phys. Rev. Lett. 58, 29 (1987).
  * [14] L. Wang, S. Cao, X. Jiang and E. K.H. Salje, J. Mech. Bio. Mat. 122 104666 (2021).
  * [15] L. I. Saliminen, A. I. Tolvanen, and M. J. Alava, Phys. Rev. Lett 89, 185503 (2003).
  * [16] A. Guarino, S. Ciliberto, A. Garciamartin, M. Zei and R. Scorretti, Eur. Phys. J. B 26, 141 (2002).
  * [17] A. B. Harris, and R. Fisch, Phys. Rev. Lett. 38, 796 (1977).
  * [18] Ashwij Mayya, Anuradha Banerjee, and R. Rajesh, Phys. Rev. E 96, 053001 (2017).
  * [19] E. Berthier, A. Mayya, and L. Ponson, J. Mech. Phys. Sol. 162, 104826 (2022).
  * [20] A. Hansen and P. C. Hemmer, Trends Stat. Phys. 1, 213 (1994).
  * [21] P. C. Hemmer and A. Hansen, J. Appl. Mech. 59 (4), 909 (1992).
  * [22] S. Pradhan, A. Hansen, and P. C. Hemmer, Phys. Rev. Lett. 95, 125501 (2005).
  * [23] M. Kloster, A. Hansen, and P. C. Hemmer, Phys. Rev. E 56, 2615 (1997).
  * [24] S. Roy and S. Biswas, Front. Phys. 9, 643602 (2021).
  * [25] S. Roy, S. Biswas, and P. Ray, Phys. Rev. E 96, 063003 (2017).
  * [26] S. Roy, and T. Hatano, P. Ray, Phys. Rev. E 105, 055003 (2022).
  * [27] S. Sinha, S. Roy, and A. Hansen, Physica A 569, 125782 (2021).
  * [28] S. Sinha, S. Roy, and A. Hansen, Phys. Rev. Research 2, 043108 (2020).
  * [29] S. Roy, S. Biswas, and P. Ray, Phys. Rev. Research 1, 033047 (2019).
  * [30] V. Kadar, G. Pal, and F. Kun, Sci. Rep. 10, 2508 (2020).
  * [31] Z. Danku, and F. Kun, Frontiers in Physics 2, 8 (2014).
  * [32] V. Kadar, Z. Danku, G. Pal, and F. Kun, Physica A 594, 127015 (2022).
  * [33] G. Pal, F. Raischel, S. Lennartz-Sassinek, F. Kun, and I. G. Main, Phys. Rev. E 93, 033006 (2016).
  * [34] X. Jiang, H. Liu, I. G. Main, and E. K. H. Salje, Phys. Rev. E 96, 023004 (2017).
  * [35] N. K. Bodaballa, S. Biswas, S. Roy, Front. Phys. 10, 768853 (2022)
  * [36] F. T. Pierce, J. Text. Inst. 17, 355 (1926).
  * [37] A. Hansen, P. C. Hemmer, and S. Pradhan, The Fibre Bundle Model, Wiley-Vch, Germany (2015).
  * [38] R. A. Fisher, and L. H. C. Tippett, Proc. Cambridge Philos. Soc. 24, 180 (1928).
  * [39] S. L. Phoenix, Adv. Appl. Probab. 11, 153 (1979).
  * [40] R. L. Smith and S. L. Phoenix, J. Appl. Mech. 48, 75 (1981).
  * [41] S. Biswas, S. Roy, and P. Ray, Phys. Rev. E 91, 050105(R) (2015).
  * [42] R. C. Hidalgo, Y. Moreno, F. Kun, and H. J. Herrmann, Phys. Rev. E 65, 046148 (2002).
  * [43] S. Roy and P. Ray, Europhysics Letters, 112(2), 26004 (2015).
  * [44] S. Roy, Phys. Rev. E 96, 042142 (2017).
  * [45] S. Roy, Front. Phys. 9, 752086 (2021).
  * [46] Y-M. Wu, S. K. Chen, T-C. Huang, H-H. Huang, W.-An Chao, and I. Koulakov, Geo. Rev. Lett. 45(4), 1832 (2018).
  * [47] G. Wergen, J. Phys. A: Math. Theor. 46, 223001 (2013).
  * [48] S. Pradhan, J. T. Kjellstadli and A. Hansen, Front. Phys. 7, 106 (2019).
  * [49] S. Roy, and T. Hatano, Phys. Rev. E 97, 062149 (2018).
  * [50] Subhadeep Roy, and Takahiro Hatano, Phys. Rev. Research 2, 023104 (2020).