# RiemannONets: Interpretable Neural Operators for Riemann Problems

Ahmad Peyvan  Vivek Oommen  Ameya D. Jagtap  George Em Karniadakis

###### Abstract

Developing the proper representations for simulating high-speed flows with
strong shock waves, rarefactions, and contact discontinuities has been a long-
standing question in numerical analysis. Herein, we employ neural operators to
solve Riemann problems encountered in compressible flows for extreme pressure
jumps (up to $10^{10}$ pressure ratio). In particular, we first consider the
DeepONet that we train in a two-stage process, following the recent work of
[1], wherein the first stage, a basis is extracted from the trunk net, which
is orthonormalized and subsequently is used in the second stage in training
the branch net. This simple modification of DeepONet has a profound effect on
its accuracy, efficiency, and robustness and leads to very accurate solutions
to Riemann problems compared to the vanilla version. It also enables us to
interpret the results physically as the hierarchical data-driven produced
basis reflects all the flow features that would otherwise be introduced using
ad hoc feature expansion layers. We also compare the results with another
neural operator based on the U-Net for low, intermediate, and very high-
pressure ratios that are very accurate for Riemann problems, especially for
large pressure ratios, due to their multiscale nature but computationally more
expensive. Overall, our study demonstrates that simple neural network
architectures, if properly pre-trained, can achieve very accurate solutions of
Riemann problems for real-time forecasting.

###### keywords:

Neural operator networks , Riemann problems , Compressible flows , DeepONet ,
U-Net , Data-driven basis

organization=Division of Applied Mathematics, 182 George Street, Brown
University, city=Providence,state=RI, postcode=02912,country=USA

##  1 Introduction

In recent years, data-driven modeling methods have shown great potential for
solving many challenging problems in the fields of computational science and
engineering. Some of these methods, in particular, can have a large impact on
large-scale computational problems and real-time forecasting as they are
efficient and, once trained, can be used for the same or similar tasks
repeatedly. Neural operators are a new paradigm for learning nonlinear
mappings between the input and output functions. In particular, there are
various neural operators available in the literature: The Deep Operator
Network (DeepONet), developed by Lu et al. [2] (first published in 2019 in
[3]), is the first proposed neural operator. DeepONet has been used to solve
many problems in computational science and engineering, such as stiff chemical
kinetics [4], multiscale bubble dynamics [5], brittle fracture analysis [6],
two-phase microstructural evolution [7] solar-thermal systems forecasting [8],
electroconvection [9], etc. In addition, several extensions of DeepONet have
been proposed in recent studies, including Partition-of-Unity (PoU) based
DeepONet [4], physics-informed DeepONet [10], DeepONet with proper orthogonal
decomposition (POD-DeepONet) [11, 12], multifidelity DeepONet [13, 14],
DeepONet with UQ [15, 16, 17] multiscale DeepONet [18], etc. Other neural
operator networks have also been proposed in the literature, such as the
Fourier neural operator (FNO) [19], the wavelet neural operator (WNO) [20],
the spectral neural operator (SNO) [21], the convolutional neural operator
[22], etc. The main difference between DeepONet and the aforementioned
operators is that DeepONet learns a new basis (through the trunk net) to
represent the operator whereas other operators, e.g., FNO, use a pre-specified
basis, e.g., Fourier expansions.

When implemented on a computer, not all models behave as operators, raising
doubts about what operator learning actually is. To address this, [23]
proposed a unifying mathematical framework for representation equivalent
neural operators (ReNO) to ensure that operations at the continuous and
discrete levels are equivalent. In their recent work, Lee and Shin [1]
proposed a novel two-step training procedure for DeepONet. The newly
introduced sequential two-step training approach begins with trunk network
training, which includes Gram-Schmidt orthonormalization via QR-factorization,
and then advances on to branch network training.

U-Net-based operators [24, 25, 26, 27] are another class of neural operators
that can be particularly effective for approximating mathematical operators
due to their inherent multi-scale nature. U-Net is a U-shaped fully
convolutional neural network [28], first proposed for biomedical image
segmentation tasks. The authors of [24, 25, 26] demonstrated that conditioning
the U-Net with respect to time can significantly improve the predictions of
the vanilla U-Net. In our current work, we extend this idea and condition the
U-Net with respect to the initial pressure and temperature states.

Neural operators have been successfully used to solve high-speed viscous
flows. In [29], Mao et al. employed DeepONet to solve hypersonic viscous flows
that give high-gradient solutions without exact shocks. In this work, we are
using DeepONet to solve the Riemann problem for the compressible Euler
equations of gas dynamics. To the best of our knowledge, this is the first
attempt to solve the Riemann problem that has discontinuous solutions using
neural operator networks. Specifically, the Riemann problem presents a
hyperbolic partial differential equation characterized by a discontinuous
initial solution. The inherent complexity of dealing with discontinuities
makes these problems particularly challenging to solve. Notably, the accuracy
of the solution tends to deteriorate in the vicinity of shocks and contact
waves.

In summary, the following are our main contributions:

  * •

Our study leverages the capabilities of deep neural operators to investigate
its efficacy in mapping input pressure ratios to the final solution at a
specified time. The final solution encompasses primitive field variables such
as density, pressure, and velocity.

  * •

We assess the impact of activation functions on prediction accuracy by
investigating both fixed and adaptive (Rowdy) activation functions [30].

  * •

We explore the performance of two training strategies for RiemannONet: the
traditional vanilla approach (one-step training) and the recently proposed
novel two-step training method [1].

  * •

We enforce positivity-preserving constraints during the training of
RiemannONet, a crucial consideration grounded in the governing physical
principles.

  * •

We systematically compare two types of neural operators for solving Riemann
problems: 1) a modified DeepONet architecture, and 2) a U-Net conditioned on
the pressure and temperature initial conditions.

  * •

We obtain interpretable basis functions for such discontinuous solutions. To
this end, we employ QR and SVD methods to investigate the solution spectrum
and diverse bases. Furthermore, we conduct a comprehensive investigation into
these basis functions, delving into the influence of network architecture on
their structure.

The structure of this paper unfolds as follows: Section 2 elucidates the
governing equations for the Riemann problem. In Section 3, we describe the two
neural operators employed in RiemannONets. Section 4 is dedicated to an in-
depth discussion of the results obtained across various test cases. Section 5
delves into the exploration of hierarchical and interpretable basis functions
for representing the two neural operators. Our conclusions and a summary of
our findings are encapsulated in Section 6.

##  2 Governing Equations

In this study, we consider the Riemann problem of the one-dimensional general
form of the hyperbolic Euler equations (1) as

$\frac{\partial\mathbf{U}(x,t)}{\partial
t}+\frac{\partial\mathbf{F}(\mathbf{U})}{\partial x}=\mathbf{S}(x,t),\quad
t\in[0,t_{f}],\quad\textrm{and}\quad\mathbf{U}(x,0)=\Bigg{\\{}\begin{array}[]{lc}\mathbf{U}_{L}&x\leq
x_{c}\\\ \mathbf{U}_{R}&x>x_{c}\end{array}$ (1)

where $t$ is time, $x$ is the spatial coordinate, $\mathbf{U}$ indicates the
vector of conservative variables, $\mathbf{F}$ denotes the corresponding
advective fluxes, and $\mathbf{S}$ is a vector of source terms. The hyperbolic
Euler equation is subjected to a discontinuous initial condition composed of
two constant states of the conservative variable vectors as $\mathbf{U}_{L}$
and $\mathbf{U}_{R}$ separated at $x=x_{c}$. We are interested in determining
the solution at a final time, as $t=t_{f}$. The definitions of the
conservative variable, flux, and source term vectors are explained below. The
Euler equations takes the form of Eq. (1) with $\mathbf{S}=\mathbf{0}$, the
conservative variables vector and flux vectors defined as

$\mathbf{U}=\begin{pmatrix}\rho\\\ \rho u\\\ \rho
E\end{pmatrix}^{T},\quad\mathbf{F}=\begin{pmatrix}\rho u\\\ \rho u^{2}+p\\\
u(\rho E+p)\end{pmatrix}^{T},$ (2)

respectively. In Eq. (2), the total energy is described as

$\rho E=\frac{p}{\gamma-1}+\frac{1}{2}u^{2},$ (3)

where $\gamma=1.4$ for all the test cases without chemistry. In Eqs. (2) and
(3), $\rho$ is density, $p$ is pressure, and $u$ indicates the velocity.

##  3 Methodology

Within this section, we introduce two distinct neural operators: DeepONet and
U-Net. We provide a comprehensive overview of the necessary modifications
essential for their effective training within the framework of RiemannONets.

###  3.1 Deep Operator Networks

The development of modern machine learning models has made it possible to
create rapid simulators for solving parametric PDEs. Instead of dealing with
the PDE explicitly, a substitute model of the PDE solution operator that can
regularly simulate PDE solutions for varied initial and boundary conditions is
frequently required, e.g. in inverse problems, in design, in uncertainty
quantification, etc. This promise is fulfilled by the novel idea of neural
operators, which was first proposed in 2019 in [3] in the form of DeepONet [2]
and is inspired by the universal approximation theorem of operators. In this
section, we will discuss the architecture of DeepONet for the Riemann problem
named RiemannONet.

Figure 1: Schematic representation of RiemannONet for the Sod’s problem. The
input to the trunk net is the sensor location in spatial dimensions, while the
input to the branch net is different realizations of left side pressure
(keeping right side pressure fixed). The output of RiemannONet is the
primitive variables consisting of velocity, density, and pressure at the final
time.

####  3.1.1 Training procedure

##### Single step

Our RiemannONet consists of two networks, namely, the trunk net and the branch
net; see Fig. 1. In the branch net, the encoded input is the different
realization of input pressure $p_{l}$ (only the left side, as the right side
has a fixed value), whereas the trunk net input consists of the spatial
dimensions $x$. The aim is to learn the solution (primitive variables) at a
later time. The output of the network is given by

$U_{r,s,v}^{NN}=\sum_{l}B_{r,l,v}\cdot T_{s,l},$ (4)

where $U_{r,s,v}^{NN}$ is the profile of the $v^{th}$ primitive variable,
$B_{r,l,v}$ is the output of the branch net, and $T_{s,l}$ is the output of
the trunk net. Here $l$ is the index of the latent dimension, which is a
hyperparameter. The RiemannONet is being trained on a data set consisting of
input-output data, which is divided into training and testing datasets. The
loss function consists of a data mismatch term:

$\mathcal{J}(\Theta)=\left|\left|\mathbf{U}^{\text{Exact}}-\mathbf{U}^{NN}(\Theta)\right|\right|_{L_{2}},$
(5)

where $\Theta$ represents all the trainable parameters of the network.

##### Two-step

The training of DeepOnet can be split into two stages employing the approach
of Lee and Shin [1]. First, the trunk net is trained by minimizing the loss
function

$\mathcal{L}^{t}(\Theta)=\left|\left|\sum_{l}T_{l}(\Theta)\mathcal{A}_{l}-\mathbf{U}^{\text{Exact}}\right|\right|_{L_{2}},$
(6)

where $\Theta$ represents all the trainable parameters of the trunk network
and $\mathcal{A}_{l}$ denotes the elements of a trainable matrix defined as
$\mathbf{A}\in\mathbb{R}^{r,vK}$ with $K$ the number of neurons of the output
layer of the trunk net. In the equation (6), the index $l$ iterates over the
neurons in the final layer of the trunk network. Let $\Theta^{*}$ and
$\mathbf{A}^{*}$ be the optimal values of the trainable parameters. The second
step consists of training the branch net by first performing a
$QR$-factorization of the trunk net represented by the matrix
$\mathbf{T}(\Theta^{*})$ using the following formula:

$\mathbf{Q}^{*}\mathbf{R}^{*}=\mathbf{qr}\left(\mathbf{T}(\Theta^{*})\right),$
(7)

where $\mathbf{qr}$ represents the QR-factorization function. In this study,
we replaced the QR-factorization with the SVD decomposition since the SVD
approach provides a unique solution and generates a hierarchical set of
orthonormal basis functions. We can construct the equivalent Q, and R matrices
for the SVD method as

$\mathbf{U}\mathbf{S}\mathbf{V}=\mathbf{svd}\left(\mathbf{T}(\Theta^{*})\right),\quad\mathbf{Q}^{*}=\mathbf{U},\textrm{
and}\quad\mathbf{R}^{*}=\mathbf{S}\mathbf{V}.$ (8)

Therefore, we can replace the QR with the SVD decomposition without any other
modification. The column of matrix $\mathbf{Q}^{*}$ forms a set of orthonormal
basis functions.

Next, the branch net is trained to match $\mathbf{R}^{*}\mathbf{\hat{A}}^{*}$
where $\mathbf{\hat{A}}$ is the reshape of matrix $\mathbf{A}$ with
$\mathbf{\hat{A}}\in\mathbb{R}^{r,K,v}$. The matrix multiplication occurs
along the second axis with a dimension of $K$. Hence, we minimize the loss
function for the training of the branch net as

$\mathcal{L}^{b}(\mu)=\left|\left|\mathbf{B}(\mu)-\mathcal{R}^{*}\mathcal{\hat{A}}^{*}\right|\right|_{L_{2}},$
(9)

where $\mathbf{B}$ represents the branch network with
$\mathbf{B}\in\mathbb{R}^{r,K,v}$. We can now construct a trained DeepONet
model by multiplying a trunk network defined as
$\mathbf{\hat{T}}=\mathbf{T}(\Theta^{*})\left(\mathbf{R}^{*}\right)^{-1}$ into
the trained branch network defined as $\mathbf{B}(\mu^{*})$. The trained
branch network is also described as $\mathbf{B}(\mu^{*})$, where $\mu^{*}$ is
the optimum value of the branch network parameters that minimize the loss
function in Eq. (9).

####  3.1.2 Rowdy activation functions

Adaptive activation functions are state-of-the-art activation functions that
can give superior performance with respect to their fixed counterparts [31].
There are various adaptive activation functions proposed for deep and physics-
informed neural networks [32]; see, for example, the work of Jagtap et al.
[33, 34]. In this work, we employed the Rowdy activation function [30]. These
adaptive activation functions have been successfully employed to solve
problems involving high-frequency complex structures. In the Rowdy activation
function, the base activation function $\phi_{1}$ is any standard activation
function, such as hyperbolic tangent. The rest of the activation functions
$\phi_{k},k=2,\cdots,K$ are defined as sine functions as follows:

$\phi_{k}(x)=n\sin((k-1)nx),$

where $n$ is the scaling factor, and here we select $n=10$. In this work, we
modify the Rowdy activation function by adding a shiftable parameter to the
base and sine functions. We choose $K=2$ and define the activation function as

$g(x)=h{\left(10a\left(x\right)+c\right)}+10a_{1}\sin{\left(10F_{1}\left(x\right)+c_{1}\right)},$
(10)

where $x$ is the input to the activation function and
$h(\left(10a\left(x\right)+c\right))$ denotes the base function of the Rowdy
activation functions that could be $\cos(\left(10a\left(x\right)+c\right))$ or
$\tanh(\left(10a\left(x\right)+c\right))$ functions for various Riemann
problems in this study. The Rowdy adaptive activation functions incorporate
five trainable parameters, including $a$, $c$, $a_{1}$, $F_{1}$, and $c_{1}$.
The $a$ and $c$ trainable coefficients are initialized using a constant value
of $0.1$. The $a_{1}$ coefficient is the amplitude hyper-parameter that is
initialized as zero. The phase shift parameter $c_{1}$ is also initialized
with the zero value. The $F_{1}$ is the frequency hyper-parameter that is
initialized with a constant value of $0.1$.

####  3.1.3 Positivity preservation for density and pressure

Ensuring the positivity of density and pressure helps maintain the physical
integrity of these conservation laws. Neglecting these constraints can lead to
physically unrealistic results and numerical instabilities in simulations and
calculations. In this work, we enforce the positivity conditions of pressure
and density during the training of the neural network. For each epoch, the
predicted density and pressure are forced to be above a small positive value,
$\epsilon=10^{-10}$. The same positivity conditions have been enforced in the
work of Jagtap et al. [35], where they employed physics-informed neural
networks to solve inverse problems in supersonic flows; see also [36, 37], and
in the work of Peyvan et al. [38] for hypersonic flow simulations.

###  3.2 U-Net with parameter conditioning

The U-Net is a multiscale convolutional network architecture that is
extensively used for solving image segmentation [28] problems. Gupta _et al._
[24] introduced the temporal conditioning mechanism that enables the U-Net to
learn time-dependent systems. The efficiency of this approach was further
extended and demonstrated in [25, 26]. In this work, we extend the same idea
and condition the U-Net with respect to the pressure states on the left side
of the initial condition profile.

The model consists of two networks. 1) a U-shaped fully convolutional neural
network; and 2) a Multi-Layered Perceptron (MLP). The MLP takes $p_{l}$ as the
input and learns a collection of functions $\vec{f}(p_{l})$ as the output. The
MLP used in this work consists of 2 hidden layers with 128 neurons and learns
non-linear functions of $P_{l}$ as follows:

$\vec{f}(p_{l})=w^{MLP}_{2}\sin(w^{MLP}_{1}\sin(w^{MLP}_{0}p_{l}+b^{MLP}_{0})+b^{MLP}_{1})+b^{MLP}_{2},$
(11)

where $\\{w^{MLP}_{i},b^{MLP}_{i}\\}_{i=0}^{2}$ represents the parameters
learned during the training of the operator.

The U-Net takes a subset of reference fields from the training dataset
($U^{ref}$) as the input. $U^{ref}$ is a collection of discrete
representations of density ($\rho$), velocity ($v$) and pressure ($p$) fields
sub-sampled from the already available training dataset. In this context,
$U^{ref}$ is a two-dimensional tensor with $C$ channels and a spatial
dimension of $W$. The U-Net learns to project $U^{ref}$ to multiple basis
functions of the same spatial resolution using the convolutional block that
comprises 1) a 1D convolutional layer [39], 2) a group normalization layer
[40], and 3) a non-linear activation layer.

The one-dimensional convolution can be performed on a two-dimensional tensor
$u$, using a three-dimensional weight tensor $w^{conv}$ with an input channel
size of $C$, output channel size of $C^{\prime}$, and a dimensionality of
$W^{\prime}$, in the following manner:

$\text{conv}(u)_{k^{\prime},i}=\sum_{k=0}^{C-1}\sum_{m=0}^{W^{\prime}-1}u_{k,i+m}.w^{conv}_{k,k^{\prime},m}$
(12)

where $w^{conv}$ is the weight tensor learned during the training process.

For the group normalization operation on a two-dimensional tensor $u$, we
separate $C$ channels to $G$ groups of $\tilde{C}$ channels
($C=G\times\tilde{C}$), compute $G$ means and standard deviations separately
as,

$\mu_{g}=\frac{1}{\tilde{C}.W}\sum_{\tilde{k}=0}^{\tilde{C}-1}\sum_{m=0}^{W-1}u_{g,\tilde{k},m}\quad\text{,
$g$=0,1,...,$G$-1},$ (13)
$\sigma_{g}^{2}=\frac{1}{\tilde{C}.W}\sum_{\tilde{k}=0}^{\tilde{C}-1}\sum_{m=0}^{W-1}(u_{g,\tilde{k},m}-\mu_{g})^{2}\quad\text{,
$g$=0,1,...,$G$-1}.$ (14)

We then normalize $u$ as,

$\hat{u}_{g,\tilde{k},i}=\frac{u_{g,\tilde{k},i}-\mu_{g}}{\sqrt{\sigma_{g}^{2}+\epsilon}}\quad\text{,
$g$=0,1,...,$G$-1}$ (15)

where $\epsilon$ is a small positive constant. The output of the group
normalization operation is,

$\text{GN}(u)_{k,i}=\gamma_{k}.\hat{u}_{k,i}+\beta_{k}\quad\text{,
$k$=0,1,...,$C$-1}.$ (16)

Here, $\gamma_{k}$ and $\beta_{k}$ are trainable $C$ dimensional parameters to
learn the ideal shift and scale operation.

The conv($u$) and GN($u$) are linear transformations of $u$. We introduce non-
linearity using Gaussian Error Linear Unit (GELU) activation function [41],
which can be approximately represented as,

$\text{GELU}(u)=0.5u(1+\tanh(\sqrt{\frac{2}{\pi}}(u+0.044715u^{3}))).$ (17)

The convolutional block that non-linearly transforms $u$ can be represented
as,

$\text{conv\\_block}(u)=\text{GELU}(\text{GN}(\text{conv}(u))).$ (18)

The downsampling operation is performed by the one-dimensional max-pooling
layer [42]. This operation is expressed as,

$\text{down}(u)_{k,i}=\max_{0\leq m<W^{\prime}}u_{k,i.S_{w}+m},$ (19)

where $S_{w}$ represents strides along width and height. In order to obtain a
scale-down factor of 2 during the downsampling operation, $S_{w}=W^{\prime}=2$
during maxpooling. The up-sampling is performed by the one-dimensional
transpose convolutional operation [43] in the following manner:

$\text{up}(u)_{k,i}=\sum_{m=0}^{W^{\prime}-1}u_{k,i.S_{w}+m}.w^{tconv}_{m},$
(20)

where $w^{tconv}_{m,n}$ is a trainable two-dimensional tensor with
dimnesionality $W^{\prime}$. To achieve a scale-up factor of 2,
$S_{w}=W^{\prime}=2$.

Borrowing the notations from [25], the latent representations of $U^{ref}$ can
be expressed as,

$\vec{z}^{\mathcal{L}_{p}}=\begin{cases}\text{conv\\_block}(U^{ref})&\text{if
$p$=1}\\\
\text{conv\\_block}(\text{down}(\vec{z}^{\mathcal{L}_{p-1}}))&\text{if
$p$=2,3,4}\end{cases}$ (21)

The $\vec{z}^{\mathcal{L}_{p}}$ can be interpreted as the learned basis
functions at 4 different scales. We further visualize the eigen spectrum and
analyze the eigen modes of the multiscale basis functions in section 5.2.

Next, we condition the latent basis functions $\vec{z}^{\mathcal{L}_{p}}$ with
the pressure initialization $P_{l}$ using an element-wise product operation as
shown below.

$\vec{z}^{\mathcal{L}_{p}}_{p_{l}}=\vec{z}^{\mathcal{L}_{p}}\odot
w^{\mathcal{L}_{p}}\vec{f}(p_{l})\\\
\implies(z^{\mathcal{L}_{p}}_{p_{l}})_{k,i}=(z^{\mathcal{L}_{p}})_{k,i}.(w^{\mathcal{L}_{p}}\vec{f}(p_{l}))_{k}\quad\forall
p,$ (22)

where $w^{\mathcal{L}_{p}}$ linearly projects $\vec{f}(p_{l})$ to number of
channels at the $p^{th}$ latent level before conditioning the latent
representation. We note that the conditioned latent variable
$\vec{z}^{\mathcal{L}_{p}}_{p_{l}}$ is discrete with respect to space and
continuous with respect to the parameter $P_{l}$. Next, we upsample the
conditioned latent vectors in the following manner,

$\vec{d}^{\mathcal{L}_{p}}_{p_{l}}=\begin{cases}\text{up}(\text{conv\\_block}(\vec{z}^{\mathcal{L}_{p+1}}_{p_{l}}))&\text{if
$p$=3}\\\
\text{up}(\text{conv\\_block}(\vec{z}^{\mathcal{L}_{p+1}}_{p_{l}}\mathrel{\ooalign{$\bigcirc$\cr$+$\cr}}\vec{d}^{\mathcal{L}_{p}}_{p_{l}}))&\text{if
$p$=1,2},\end{cases}$ (23)

where $\mathrel{\ooalign{$\bigcirc$\cr$+$\cr}}$ concatenates the tensors along
the dimension of channels.

The output of the model is computed as,

$U(P_{l})=\text{GN}(\text{conv}(\text{conv\\_block}(\vec{z}^{\mathcal{L}_{1}}_{p_{l}}\mathrel{\ooalign{$\bigcirc$\cr$+$\cr}}\vec{d}^{\mathcal{L}_{1}}_{p_{l}}))).$
(24)

A schematic of the architecture used in this study is presented in Fig. 2.

Figure 2: U-Net conditioned on pressure initialization ($p_{l}$).
$U^{ref}\subseteq U^{\text{train}}$ is provided as input to a U-Net which
behaves like a multi-scale neural operator. The output of each encoder block,
$\vec{z}^{\mathcal{L}_{p}}$, is conditioned on the parameter $p_{l}$ through
an element-wise product operation. The corresponding representation is
concatenated with the previous decoder block’s output
($\vec{d}^{\mathcal{L}_{p-1}}_{p_{l}}$), and subsequently projected back as
the output of the model.

##  4 Results

In the subsequent test cases, we utilize the RiemannONets to establish a
mapping from the initial pressure value on the left side of the membrane to
the density, velocity, and pressure profiles at a final time. This mapping is
performed across 200 equidistant coordinate points situated within the one-
dimensional physical domain. The RiemannONets refer to either DeepONet or
U-Net architectures. Table 1 presents an overview of the results of this study
and provides a comprehensive comparison of several deep neural operators
trained and applied for the inference of three Riemann problems. Here, we
investigate low-pressure ratio (LPR), intermediate-pressure ratio (IPR), and
high-pressure ratio (HPR) Sod problems. For each test case in Table 1, we
train the neural operators with ten different initializations of weights and
biases and present the mean and standard deviation of the $L_{2}$ norm of the
error across the ten ensembles. The mean training time of 10 ensembles is also
shown in Table 1. The training time is computed by deploying the code on the
NVIDIA RTX 3090 GPUs with Ampere architecture and the on-card memory of 24
GBs. Both branch and trunk nets consist of five hidden layers, each with a
width of 150 neurons for the 2-step training DeepONet and 50 for the vanilla
DeepONet. Using a single shared trunk neural net, we set up the branch net to
infer three variables, including density, velocity, and pressure. The relative
$L_{2}$ norm of the error is computed over the entire test data set. Each
sample in the test data set contains 200 density, velocity, and pressure
values corresponding to 200 coordinate points. Let $\mathcal{D}$ be a third-
order tensor indicating the inferred solution of the test data set with the
shape of $\left(N_{s},N_{p},3\right)$, where $N_{s}=100$ is the number of
samples, $N_{p}=200$ is the number of points, and 3 referring to the density,
velocity, and pressure. The relative $L_{2}$ norm for each quantity of
interest is computed as:

$L_{2}\left(E_{k}\right)=\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\frac{\sqrt{\sum_{j=1}^{N_{p}}\left(D_{i,j,k}-G_{i,j,k}\right)^{2}}}{\sqrt{(\sum_{j=1}^{N_{p}}G_{i,j,k}^{2})}},\quad
k=1,2,3,$ (25)

where $E_{k}$ denotes the point-wise error of predicted density ($k=1$),
velocity ($k=2$), and pressure $k=1$). In Eq. (25), $G_{i,j,k}$ refers to the
ground truth for $i^{\textrm{th}}$ test sample at $j^{\textrm{th}}$ coordinate
for density, velocity, or pressure solution. The total relative $L_{2}$ norm
of the error is the mean of the $L_{2}$ norm of error in density, velocity,
and pressure.

Table 1: Relative L2 norms mean and standard deviation obtained using 10 runs. The $L_{2}$ norm of the error is calculated over the entire testing dataset for density, velocity, and pressure profiles. The time reported is the training time; the inference time is negligible. Cases |  $L_{2}(\rho)$ % |  $L_{2}(u)$ % |  $L_{2}(p)$ % |  total $L_{2}$ norm % | Time (Min)
---|---|---|---|---|---
LPR(1 step Tanh) | $0.96\pm 0.064$ | $3.92\pm 0.251$ | $0.74\pm 0.044$ | $1.88\pm 0.120$ | $12.95$
LPR(1 step Rowdy) | $\mathbf{0.70\pm 0.103}$ | $\mathbf{2.57\pm 0.580}$ | $\mathbf{0.53\pm 0.081}$ | $\mathbf{1.27\pm 0.254}$ | $38.07$
LPR(2 step Rowdy) | $\mathbf{0.41\pm 0.017}$ | $1.28\pm 0.119$ | $\mathbf{0.33\pm 0.047}$ | $0.67\pm 0.061$ | $33.45$
LPR(U-Net) | $0.49\pm 0.070$ | $\mathbf{0.86\pm 0.129}$ | $0.36\pm 0.094$ | $\mathbf{0.57\pm 0.098}$ | 1388.58
IPR(2 step Rowdy) | $\mathbf{0.33\pm 0.027}$ | $\mathbf{0.86\pm 0.071}$ | $\mathbf{0.20\pm 0.030}$ | $\mathbf{0.46\pm 0.043}$ | 33.06
IPR(U-Net) | $0.48\pm 0.069$ | $0.98\pm 0.117$ | $0.50\pm 0.156$ | $0.66\pm 0.114$ | 1257.57
HPR(2 step Rowdy(QR)) | $0.71\pm 0.113$ | $3.50\pm 0.248$ | $3.70\pm 3.79$ | $2.64\pm 1.38$ | $34.60$
HPR(2 step Rowdy(SVD)) | $\mathbf{0.66\pm 0.093}$ | $3.39\pm 0.104$ | $2.86\pm 1.680$ | $2.31\pm 0.626$ | $26.98$
HPR(U-Net) | $1.00\pm 0.208$ | $\mathbf{2.65\pm 0.115}$ | $\mathbf{2.27\pm 0.457}$ | $\mathbf{1.97\pm 0.260}$ | 1235.67

In Table 1, we first compare the accuracy of the vanilla DeepONet framework
with adaptive and normal activation functions. We employed the Rowdy adaptive
activation function described in Eq. (10) with a $\tanh$ base activation
function. For the first comparison, the (LPR) test problem is used. For the
vanilla DeepONet, we use weight ($L_{2}$) regularization and apply density and
pressure constraint positivity during training. For all the variables, the
vanilla DeepONet with Rowdy activation function performs superior to the
DeepONet with $\tanh$ activation function. The training time for the Rowdy
activation function is longer than $\tanh$ since the optimizer must minimize
the loss function for five additional coefficients of the adaptive activation
function besides the weights and biases.

We have also compared the accuracy of one-step Rowdy and two-step Rowdy
DeepONet approaches for the LPR case. We can conclude that the two-step
training significantly improves the inference accuracy while requiring less
training time than the vanilla DeepONet. We can observe that for the LPR case,
the U-Net accuracy in inferring velocity and pressure is better than the two-
step Rowdy approach. In comparison, the two-step Rowdy provides greater
accuracy for density than the U-Net. The training time for the two-step Rowdy
approach is significantly less than the U-Net, while the accuracy is
comparable. For the IPR case, the two-step Rowdy method obtains better
accuracy for all the variables. The two-step Rowdy method provides a small
standard deviation for the IPR cases, indicating that the two-step DeepONet is
independent of the initialization point. For the HPR test case, we have
compared the accuracy of the two-step training approach with SVD and QR-
factorization of the trunk net. The SVD factorization demonstrates higher
accuracy compared to the QR factorization. Here, we also compare the accuracy
of the two-step approach with U-Net. The U-Net provides higher accuracy for
velocity and pressure and a higher error for density than the two-step
training with the SVD approach. For all the test cases, the computational cost
of U-Net is much higher than the two-step training DeepONet approach.

|
---|---
(a) $\rho_{\textrm{test}}$, $\tanh$ | (b) $\rho_{\textrm{test}}$, Rowdy
Figure 3: Low pressure ratio test case: Comparison of $\tanh$ versus Rowdy
$\tanh$ adaptive activation functions The density of four samples is inferred
from the testing data set. The predictive accuracy of the Rowdy $\tanh$
adaptive activation functions is better than that of its fixed counterpart.

###  4.1 Low pressure ratio Sod problem

In this case, we solve Eq. (1) on a spatial domain defined such as
$x\in[0,1]$. The initial conditions for primitive variables are imposed as

$\left(\rho,u,p\right)=\begin{cases}\left(1.0,0.0,p_{l}\right)&x\leq 0.5\\\
\left(0.125,0.0,0.1\right)&x>0.5,\end{cases}$ (26)

where $p_{l}\in[1.0,5.0]$. We employed the analytical method described in [44]
to obtain the results for primitive variables at $t_{f}=0.1$. From the 500
cases we randomly choose 400 trajectories for training and 100 trajectories
for testing data-sets.

|  |
---|---|---
(a) $\rho_{\textrm{test}}$, DeepONet | (b) $u_{\textrm{test}}$, DeepONet | (c) $p_{\textrm{test}}$, DeepONet
|  |
(d) $\rho_{\textrm{test}}$, U-Net | (e) $u_{\textrm{test}}$, U-Net | (f) $p_{\textrm{test}}$, U-Net
Figure 4: Low Pressure Ratio Sod Problem: Comparison of DeepONet and U-Net
results. The first row shows the DeepONet results for density, velocity, and
pressure, whereas the second row shows the corresponding results for U-Net.
The density of four samples is inferred from the testing data set.

We explore the utilization of adaptive activation functions in the DeepONet
framework in contrast to conventional $\tanh$ activation functions. Figure 3
illustrates the testing density results for a low-pressure ratio. A standard
activation function is employed in the left plot, while the adaptive Rowdy
activation function is used in the right plot. It is evident that, compared to
the Rowdy activation functions, the fixed activation function introduces more
oscillations, particularly in the vicinity of discontinuous solutions such as
contact and shock waves. Despite the increased computational cost associated
with Rowdy activation, relying solely on fixed activation functions for
training is unlikely to enhance the accuracy of predictions. Henceforth, we
have employed the Rowdy activation function for all subsequent test cases.

The results of the DeepONet and U-Net models are illustrated in Figure 4.
Specifically, the first row displays the outcomes of DeepONet for density,
velocity, and pressure, while the second row presents the corresponding
results for U-Net. The density for four distinct samples is inferred from the
testing dataset. Notably, both DeepONet and U-Net exhibit good performance,
particularly in addressing the challenges posed by the low-pressure ratio
problem. In particular, while capturing the shock, contact, and expansion
waves, both architectures exhibit negligible overshoots and undershoots.
Focusing on the density profiles, we can observe that the U-Net provides some
low-amplitude oscillations in the blue curve compared to the DeepONet
counterpart. However, DeepONet induces an overshoot in the velocity profile of
$pl=1.15$ at the shock wave location. The visual analysis from the figure
matches the $L_{2}$ norm results shown in Table 1.

|  |
---|---|---
(a) $\rho_{\textrm{test}}$, DeepONet | (b) $u_{\textrm{test}}$, DeepONet | (c) $p_{\textrm{test}}$, DeepONet
|  |
(d) $\rho_{\textrm{test}}$, U-Net | (e) $u_{\textrm{test}}$, U-Net | (f) $p_{\textrm{test}}$, U-Net
Figure 5: Intermediate Pressure Ratio Sod Problem: Comparison of DeepONet and
U-Net results. The first row shows the DeepONet results for density, velocity,
and pressure, whereas the second row shows the corresponding results for
U-Net. The density of four samples is inferred from the testing data set.

###  4.2 Medium pressure ratio Sod problem

Next, we select a more challenging Sod problem for which we increase the
pressure and density ratio. We solve Eq. (1) on a spatial domain defined as
$x\in[-1,2]$. For this case, the initial conditions are selected as

$\left(\rho,u,p\right)=\begin{cases}\left(2.0,0.0,p_{l}\right)&x\leq 0.5\\\
\left(0.125,0.0,0.1\right)&x>0.5,\end{cases}$ (27)

where $p_{l}\in[50.0,100.0]$. We again use the exact method to obtain the
results at $t_{f}=0.1$. We choose 500 equispaced various $p_{l}$ values and
randomly assign 400 cases for training and 100 for testing.

Figure 5 depicts the outcomes obtained from DeepONet and U-Net for the IPR
problem, showcasing the testing results for density, velocity, and pressure.
The initial row presents the results of DeepONet for density, velocity, and
pressure, while the subsequent row displays the corresponding outcomes for
U-Net. The density values for four distinct samples are extrapolated from the
testing data set. Notably, in this instance, the accuracy of DeepONet is
slightly better than that of U-Net, particularly evident in the absence of
oscillations near the shock wave location at the right-most discontinuities in
the curves.

|  |
---|---|---
(a) $\rho_{\textrm{test}}$, DeepONet | (b) $u_{\textrm{test}}$, DeepONet | (c) $p_{\textrm{test}}$, DeepONet
|  |
(d) $\rho_{\textrm{test}}$, U-Net | (e) $u_{\textrm{test}}$, U-Net | (f) $p_{\textrm{test}}$, U-Net
Figure 6: High Pressure Ratio Sod Problem: Comparison of DeepONet and U-Net
results. The first row shows the DeepONet results for density, velocity, and
pressure, whereas the second row shows the corresponding results for U-Net.
The density of four samples is inferred from the testing data set.

###  4.3 High pressure ratio Sod problem (LeBlanc problem)

For the last test case, we increase the pressure ratio to the extreme values
employed in the so-called LeBlanc-Sod problem. We solve Eq. (1) on a spatial
domain defined as $x\in[-20,20]$. For this case, the initial conditions are
selected as

$\left(\rho,u,p\right)=\begin{cases}\left(2.0,0.0,p_{l}\right)&x\leq-10\\\
\left(0.001,0.0,1.0\right)&x>-10,\end{cases}$ (28)

where $p_{l}\in[10^{9},10^{10}]$. We use the exact method to obtain the
results at $t_{f}=0.0001$. We choose 500 equispaced various $p_{l}$ values and
randomly assign 400 cases for training and 100 for testing. For training the
neural nets on this problem, we employ the logarithm of density and pressure
values to construct the loss function. For the inference stage, we apply the
exponential function on the density and pressure to convert the predicted
values to physical values of pressure and density.

The LeBlanc problem poses a significant challenge due to a very high-pressure
ratio. Figure 6 presents the inferred results of DeepONet and U-Net for
density, velocity, and pressure. Notably, these results highlight the testing
outcomes. U-Net adeptly captures discontinuous features such as shocks and
contacts with higher accuracy (without oscillations) than the DeepONet.

##  5 Operator Representation

This section is dedicated to the analysis of data-driven basis functions,
forming the foundation for obtaining representations of the two neural
operators. It’s worth highlighting that the basis for DeepOnet is continuous,
while for U-Net, it takes a discrete form.

|  |
---|---|---
(a) Spectrum | (b) Mode$=0$ basis | (c) Mode$=10$ basis
|  |
(d) Mode$=30$ basis | (e) Mode$=80$ basis | (f) Mode$=120$ basis
Figure 7: QR vs. SVD spectrum and basis comparison: The top-left plot shows
the spectrum comparison between QR and SVD. The rest of the plots compare the
different modes of QR and SVD.

###  5.1 Basis functions - DeepONet

In Fig. 7 we compare the spectrum of singular values obtained using QR and SVD
decomposition of the trunk net output using the optimized trunk parameters for
the training of the HPR problem. The QR eigenvalues are sorted after creating
a positive set of eigenvalues by adding the minimum eigenvalue’s absolute
value, which is negative to all the modes eigenvalues. Some of the basis
functions for five different modes are depicted for comparison. We can observe
that the SVD basis functions exhibit a hierarchical structure from the lowest
to highest modes, where we observe high-frequency oscillations in the basis
function profile. In contrast, the basis functions computed using the QR
factorization are not hierarchical since mode zero exhibits high-frequency
features like the higher modes. Moreover, the QR spectrum is flat across all
the modes, inferring that the contribution of all the modes is equal. Both
sets of 150 eigenfunctions obtained from QR and SVD decomposition are
orthonormal; however, the spectrum of SVD approaches zero at the highest mode,
indicating a smaller amount of oscillations than the QR spectrum. Therefore,
the SVD factorization is more robust than the QR factorization and can be used
for feature extraction since the low-mode eigenfunctions mimic the physical
shape of the solution profiles.

|  |
---|---|---
(a) LPR | (b) IPR | (c) HPR
Figure 8: Effect of jump in pressure ratio: We compare the third basis for low
(a), intermediate (b), and high pressure (c) ratios.

In Fig. 8, we compare the third mode basis functions for three problems from
low to high-pressure ratios. The pressure ratio range varies such that
$p_{l}\in[10,50]$, $p_{l}\in[500,1000]$, and $p_{l}\in[10^{9},10^{10}]$ for
LPR, IPR, and HPR problems, respectively. We can observe that the shape of the
low modes remains approximately the same from LPR to IPR cases, indicating
that the data-driven basis function obtained using the low-pressure training
data set can train a branch net for the intermediate pressure data set for at
least the low modes.

|  |
---|---|---
(a) Spectrum | (b) Mode$=0$ basis | (c) Mode$=10$ basis
|  |
(d) Mode$=30$ basis | (e) Mode$=80$ basis | (f) Mode$=120$ basis
Figure 9: The top-left plot shows the SVD spectrum for the intermediate
pressure ratio, and rest of the plots shows different SVD basis functions.

In Fig. 9 we present the eigenvalue spectrum of the basis functions for the
IPR problem across 150 modes corresponding to the number of neurons in the
last layer of the trunk net. According to Fig. 9(a), the eigenvalues decay
when the mode number increases. The descending trend of the spectrum describes
the shapes of the basis functions for modes $0$, $10$, $30$, $80$, and $120$.
From mode $0$, it is evident that the expansion waves are the first feature
that are extracted from the data set. The highest-mode basis functions capture
the wide range of frequencies exhibited at discontinuities such as contact and
shock waves. Considering Fig. 9(c)-(e), we can realize the fact that these
modes are hierarchically contributing to resolving the shock wave feature. The
order of the SVD orthonormal basis functions can be employed to construct the
inferred solution manually using a specific number of high modes to avoid
oscillations at the discontinuities of the solution. In the following section,
we explain an approach that modifies the number of basis functions to infer
the solution accurately.

|  |
---|---|---
(a) Spectrum | (b) Mode$=0$ basis | (c) Mode$=1$ basis
|  |
(d) Mode$=10$ basis | (e) Mode$=20$ basis | (f) Mode$=49$ basis
Figure 10:  Comparison of the SVD spectrum and the basis functions for
different number of neurons. It is evident that as the number of neurons
increases, more high frequencies are captured, which is useful for accurately
resolving shock structure present in the solution.

Fig. 10(a)-(f) depicts the spectrum of eigenvalues for a trunk net constructed
by five hidden layers with a constant width of 50, 100, and 150 neurons
corresponding to 50, 100, and 150 eigenvalues and basis functions. By
comparing the spectrum of eigenvalues, we observe that using a higher number
of modes can recover a higher content of information at the higher modes.
Employing a higher number of neurons in the trunk net layers decreases the
$L_{2}$ norm of the error, but it could induce oscillations from the highest
modes. According to Fig. 10(b) and (c), the first feature extracted from the
data is the expansion wave, while the second one is the shape of the velocity
profile. Considering Fig. 10(e) and (f), we can deduce that by increasing the
number of neurons, we are increasing the capacity of capturing the shock wave
more accurately since we can observe high-frequency jumps in the basis
function of 50 neuron cases close to the location of the shock wave near the
right boundary exits for Mode$=49$. At the same time, there is no high-
frequency oscillation for the 100 and 150-neuron cases. The shock wave feature
is captured using modes higher than $49$ in the 100 and 150-neuron cases.

|
---|---
(a) Spectrum | (b) Mode$=0$ basis
|
(c) Mode$=1$ basis | (d) Mode$=2$ basis
Figure 11: Layer-wise contribution to the spectrum and the basis: The figure
shows the contribution of each hidden layer to the final basis. We have
plotted the contribution of the first five layers. We see that with each
layer, more high-frequencies are captured, and the final layer (layer 5)
output (excluding the linear layer) is almost similar to the basis function.

In Fig. 11, we attempt to reveal the contribution of each hidden layer to the
shape of the final basis functions of the HPR problem. For this purpose, we
performed an SVD decomposition of the output of hidden layers from the first
to the last linear layer. Using the eigenvalues of each layer output, we can
observe the flow of information learned by the trunk net. Considering Fig.
11(a), the red line shows the spectrum of eigenvalues of the output of the
last linear layer of the trunk net. The label “$m$ layers” refers to the
output of the $m^{th}$ hidden layer. The spectrum shows that the first layers
always contribute to learning the lowest-mode basis functions. According to
Fig. 11(a), the flow of information does not reach the highest modes until the
$4^{th}$ hidden layer, where we can see there is no plateau region in the
spectrum of the output of the first $4$, first $5$ and last linear layer. In
Fig. 11(b)-(d), the shape of the first three modes out of the layers is shown.
From Fig. 11(b), we see that the hidden layers keep adding more information as
we move from the trunk input to the trunk net output.

|
---|---
(a) Spectrum rectangle architecture | (b) Spectrum cone architecture
|
(c) Mode$=0$ basis rectangle architecture | (d) Mode$=0$ basis cone architecture
Figure 12: Layer-wise contribution to the spectrum and the basis for
rectangular and decoder architectures: The decoder architecture refers to a
trunk shape such as [30, 70, 108, 150, 150] and rectangle network refers to
[150, 150, 150, 150, 150] hidden layers arrangement.

According to Fig. 11(a), we can use a smaller number of neurons for the first
hidden layer of the neural net as the eigenvalues of most of the neurons are
approximately zero. So, we can generate a cone-shaped configuration of the
hidden layers according to the non-zero eigenvalues and use the downsized
trunk net for the training procedure. To this end, we chose
$[30,70,108,150,150]$ number of neurons for the hidden layers of the trunk
net. We computed the optimal number of neurons for each layer based on the
points on the spectrum that reach the plateau state. A comparison of the
original rectangle shape versus cone shape architectures is depicted in Fig.
12. We are comparing the rectangle with a cone shape (Fig. 12(a) and (b); we
can see that when we use the cone shape type, the flow of information to the
last layers gets delayed for higher modes. This fact can be seen by comparing
the shape of the basis functions of the first $5$ hidden layers. We can see
that the shape of the basis functions changes for higher modes in the
rectangle architecture.

|  |
---|---|---
(a) Branch coefficients $\rho$ train | (b) Branch coefficients $u$ train | (c) Branch coefficients $p$ train
|  |
(d) Branch coefficients $\rho$ test | (e) Branch coefficients $u$ test | (f) Branch coefficients $p$ test
Figure 13: Contribution of coefficients of basis functions in the final
inferred solution of the IPR problem: We plotted the value of branch
coefficients for the intermediate pressure ratio case. The top row shows the
basis coefficients for training, and the bottom row shows the basis
coefficients for testing. The first, second, and third columns give the basis
coefficients for density, velocity, and pressure. The insets show
representative basis spatial modes.

Figure 13 depicts the value of the coefficients of the basis functions that
are predicted by the branch net. The bar plots show the values on a log scale
for the vertical scale. The values are sorted from lowest to highest modes for
density, velocity, and pressure. The first row depicts the branch coefficients
for the training data set, and the second row shows the branch coefficients
for the test data set. The inset plots show the shape of the basis function at
a particular mode. Here, we are trying to investigate the contribution of
various features to the inferred solution of the IPR problem. According to
Fig.13(a), the highest contribution to the density profile comes from mode $0$
while the highest contribution for the velocity is at mode $1$ (see Fig. 13
(b) and (e)). Another interesting fact is the large contribution of the
highest modes in the inferred profile of density, velocity, and pressure. This
large contribution results from predicting a discontinuous solution consisting
of shock and contact waves.

|  |
---|---|---
(a) Spectrum | (b) Mode$=0$ basis | (c) Mode$=10$ basis
|  |
(d) Mode$=20$ basis | (e) Mode$=30$ basis | (f) Mode$=40$ basis
Figure 14: (a) Eigenspectra of U-Net at $\mathcal{L}_{1}$, $\mathcal{L}_{2}$,
$\mathcal{L}_{3}$, $\mathcal{L}_{4}$. (b)-(f) Orthogonally decomposed
representation of the basis learned by the U-Net at different latent levels.

###  5.2 Basis functions - U-Net

In this section, we visualize and analyze the basis functions learned by the
U-Net conditioned on $p_{l}$. Specifically, we consider the IPR test case and
train the U-Net to learn the dependence of the density field ($\rho$) on
$p_{l}$. From Eq. (21), $\vec{z}^{\mathcal{L}_{p}}$ represents the discrete
basis functions learned by the U-Net at the $p^{th}$ latent level spanning
$\mathcal{L}_{p}$. We perform SVD on $\vec{z}^{\mathcal{L}_{p}}$ separately
for $p=1,2,3,4$ and visualize the eigenvalues and eigenmodes.

$\vec{z}^{\mathcal{L}_{p}}=U_{p}\Sigma_{p}V_{p}^{T}$ (29)

where $\Sigma_{p}$ represents the ordered list of eigenvalues and $V_{p}$
represents the corresponding eigenmodes at the latent space $\mathcal{L}_{p}$.
We visualize the decay of eigen values $\Sigma_{p}$ for $p=1,2,3,4$ in Fig. 14
(a). We also plot the scaled basis functions for eigen modes = 0, 10, 20, 30,
40 for all $p$ in Fig. 14 (b)-(f).

In Fig. 14(a) we observe that the rate of decay of eigenvalues is greater for
$\mathcal{L}_{4}$ than $\mathcal{L}_{3}$ than $\mathcal{L}_{2}$ than
$\mathcal{L}_{1}$. Therefore, the basis functions $\vec{z}^{\mathcal{L}_{p}}$
with a larger $p$ only learn the high energy modes with lower frequencies, and
subsequently, the basis functions with a lower $p$ become increasingly
responsible for learning the lower energy modes that carry higher frequencies.
Furthermore, from Fig. 14 (e) and (f), the basis functions corresponding to
$\mathcal{L}_{4}$ are constant, which again indicates its inability to capture
high frequencies. The 0th modes of $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$
shown in Fig. 14 (b) have a striking similarity with the pattern of the
density fields learned by the U-Net, as well as the 0th mode learned by
DeepONet shown in Fig. 9(b).

###  5.3 Constructing optimal sets of basis functions for higher accuracy

We investigate the effect of the trunk net width on the accuracy of the two-
step approach to learning the solution to the HPR problem. According to Table
2, we can conclude that increasing the number of neurons in the trunk network
can lead to higher accuracy. Figure 10 shows that the accuracy improvement is
achieved by incorporating higher modes basis functions to construct the
discontinuous solution of the Riemann problems. However, by increasing the
number of neurons, we are adding high oscillatory basis functions to the set,
which can give rise to unwanted oscillations in the inferred solution.
Therefore, at the inference stage, we can take advantage of the hierarchical
features of the SVD basis functions and remove higher modes one by one from
the basis set. We can monitor the corresponding accuracy by computing the
$L_{2}$ norm of the solution predicted by the truncated set of basis
functions. We then use the set of basis functions that provides the lowest
$L_{2}$ norm of the error. In brief, we can design a procedure to use the two-
step DeepONet for the best-inferred solution of the Riemann problems as
follows:

Table 2: Relative L2 norms mean and standard deviation obtained using 10 runs. The $L_{2}$ norm of the error is calculated over the entire testing dataset for density, velocity, and pressure profiles for the HPR test cases using 50, 100, and 150 neurons for the width of the trunk net. The time reported is the training time; the inference time is negligible. Cases |  $L_{2}(\rho)$ % |  $L_{2}(u)$ % |  $L_{2}(p)$ % |  total $L_{2}$ norm % | Time (Min)
---|---|---|---|---|---
HPR(2 step Rowdy(50)) | $2.96\pm 0.050$ | $6.51\pm 0.046$ | $8.47\pm 0.147$ | $5.98\pm 0.080$ | $33.23$
HPR(2 step Rowdy(100)) | $0.92\pm 0.043$ | $3.83\pm 0.147$ | $3.28\pm 0.839$ | $2.67\pm 0.343$ | $34.51$
HPR(2 step Rowdy(150)) | $\mathbf{0.66\pm 0.093}$ | $\mathbf{3.39\pm 0.104}$ | $\mathbf{2.86\pm 1.680}$ | $\mathbf{2.31\pm 0.626}$ | $26.98$

  1. 1.

Increase the number of neurons at the last layer of the trunk net to obtain
the lowest value of loss during training.

  2. 2.

If oscillations are near the discontinuous solution, we remove the highest
basis functions from the set and compute the optimal error value.

  3. 3.

Use the optimum number of basis functions for the trained model for future
inference.

Figure 15:  The $L_{2}$ norm of the error in density, velocity, and pressure
for the HPR problem by using optimal number of basis functions to avoid
numerical artifact due to the existence of high mode basis functions in the
orthonormal set.

We have followed this procedure on the HPR test case with ten ensembles to
ascertain our claim. First, we took the best-trained DeepONet from the set of
10 ensembles, removed the highest mode basis functions from the set, and used
the remaining basis functions to infer the solution. We then computed the
$L_{2}$ norm of the variables using the first 60 to 150 basis functions from
the orthonormal basis set. The results of the errors are depicted in Fig. 15.
We can observe that the optimal number of basis functions to use here is 114
first basis functions. We applied the approach to the ten ensembles of the two
step-trained DeepONets and recalculated the $\%L_{2}(\rho)=0.573\pm 0.091$
,$\%L_{2}(u)=3.39\pm 0.105$, and $\%L_{2}(p)=2.29\pm 1.78$. We improved the
accuracy of $\rho$ by $13.6$ percent and of $p$ by $20.0$ percent. The
accuracy of the velocity remained intact, as we may need to use an even higher
number than 150 neurons to see a similar trend as in pressure and density
inference results.

##  6 Summary

Enhancing the prediction of solutions for high-speed flows governed by
compressible Euler equations holds significant implications for the design of
aerospace vehicles, including airplanes, re-entry vehicles, missiles, and
more. In this study, we leveraged the properties of deep neural operators –
RiemannONets – to tackle Riemann problems, crucial for simulating high-speed
flows. These problems entail discontinuous solutions, such as shocks and
contact discontinuities, representing some of the most challenging aspects in
the realm of scientific computing.

We devised the RiemannONets by incorporating two distinct neural operators.
The first operator is built upon DeepONet, which underwent modifications for a
two-stage training approach, enhancing prediction accuracy. The second
operator is a U-Net, tailored to be conditioned on pressure initialization.
Our training and testing of RiemannONets focused on the Sod shock tube
problem, encompassing pressure ratios ranging from 10 to $10^{10}$.
Specifically, RiemannONets were trained on input-output data sets,
establishing a mapping from the initial solution to the final time step. Once
trained, RiemannONets can predict solutions for unseen datasets in real-time,
requiring no additional optimization. The predictions from both neural
operators exhibited remarkable accuracy across low, intermediate, and very
high pressure ratios, with an error margin below $2\%$. However, it is
noteworthy that U-Net’s computational speed lags behind DeepONet by orders of
magnitude. These results correspond to a single mapping from the initial
condition to the final time. We have also shown that the use of adaptive
activation function in the structure of the DeepONet increases the accuracy of
the prediction comparing to the fixed activation functions. Notably, we
achieved similarly accurate results for the time-evolving Sod problem using
the DeepONet trained in two stages.

We systematically explored the basis functions generated by the trunk net of
DeepONet for representing the operator in a continuous manner. Using the
orthonormalization process, we constructed a set of orthonormal basis
functions. For orthonormalization, we performed QR-factorization and singular
value decomposition (SVD) on the output of trunk net. Our in-depth
investigation of the data-driven basis functions has led us to the following
conclusions:

  * •

The SVD decomposition resulted in a hierarchical orthonormal basis with
distinctive features, which can be used to remove Gibbs oscillations at the
discontinuities during the inference.

  * •

A comparison of the eigenvalue spectra of QR and SVD revealed that the QR
eigenvalues are equally contributing to the final solution, while the SVD
eigenvalues show a descending trend in contribution from low to high modes.

  * •

The unique basis of the SVD decomposition proves particularly advantageous for
high-speed flows.

  * •

The SVD basis functions exhibit a similar shape for low, intermediate, and
high-pressure ratio Sod problems. Consequently, we can utilize a trunk network
trained on a fixed range of pressure ratios to train the branch network for
higher pressure ratios ranges.

  * •

Using a larger width of the trunk network in the two-step training improves
the accuracy of the inference by learning a higher amount of information at
high modes.

  * •

The first hidden layers of the trunk network are responsible for learning low
modes while the later hidden layers contribute to learning the high-frequency
features.

  * •

The specific contribution of each basis function in constructing density,
velocity, and pressure fields are explored, revealing the distinct hierarchy
of feature learning for density, velocity, and pressure profiles.

  * •

Employing the hierarchy of the SVD basis functions, we constructed a procedure
that can effectively remove high frequency artifacts near the discontinuous
region of the Riemann problem solutions.

At present, there are no numerical methods in existence that can attain the
SVD data-driven orthonormal basis functions. This requires the incorporation
of ad hoc features in the basis functions, such as the utilization of
enrichment methods in finite elements or the integration of a feature layer in
a neural network to address discontinuous, singular, or multiscale solutions.
Additionally, we visualize the basis functions learned by different latent
levels of the U-Net. This rvealed that the first mode of U-Net is similar to
the first mode of DeepOnet. Moreover, the lower latent layers focus on
learning the low modes whereas the higher levels capture the highest modes.

Our ongoing research is focused on examining basis functions tailored for
addressing two- and three-dimensional high-speed flow problems characterized
by discontinuous solutions.

## Acknowledgments

This work was supported by the U.S. Army Research Laboratory W911NF-22-2-0047
and by the MURI-AFOSR FA9550-20-1-0358.

## References

  * [1] S. Lee, Y. Shin, On the training and generalization of deep operator networks, arXiv preprint arXiv:2309.01020 (2023).
  * [2] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature machine intelligence 3 (3) (2021) 218–229.
  * [3] L. Lu, P. Jin, G. E. Karniadakis, Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators, arXiv preprint arXiv:1910.03193 (2019).
  * [4] S. Goswami, A. D. Jagtap, H. Babaee, B. T. Susi, G. E. Karniadakis, Learning stiff chemical kinetics using extended deep neural operators, Computer Methods in Applied Mechanics and Engineering 419 (2024) 116674.
  * [5] C. Lin, Z. Li, L. Lu, S. Cai, M. Maxey, G. E. Karniadakis, Operator learning for predicting multiscale bubble growth dynamics, The Journal of Chemical Physics 154 (10) (2021).
  * [6] S. Goswami, M. Yin, Y. Yu, G. E. Karniadakis, A physics-informed variational deeponet for predicting crack path in quasi-brittle materials, Computer Methods in Applied Mechanics and Engineering 391 (2022) 114587.
  * [7] V. Oommen, K. Shukla, S. Goswami, R. Dingreville, G. E. Karniadakis, Learning two-phase microstructure evolution using neural operators and autoencoder architectures, npj Computational Materials 8 (1) (2022) 190.
  * [8] J. D. Osorio, Z. Wang, G. Karniadakis, S. Cai, C. Chryssostomidis, M. Panwar, R. Hovsapian, Forecasting solar-thermal systems performance under transient operation using a data-driven machine learning approach based on the deep operator network architecture, Energy Conversion and Management 252 (2022) 115063\.
  * [9] S. Cai, Z. Wang, L. Lu, T. A. Zaki, G. E. Karniadakis, DeepM&Mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks, Journal of Computational Physics 436 (2021) 110296.
  * [10] S. Wang, H. Wang, P. Perdikaris, Learning the solution operator of parametric partial differential equations with physics-informed DeepONets, Science advances 7 (40) (2021) eabi8605.
  * [11] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, G. E. Karniadakis, A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data, Computer Methods in Applied Mechanics and Engineering 393 (2022) 114778.
  * [12] S. Venturi, T. Casey, SVD perspectives for augmenting DeepONet flexibility and interpretability, Computer Methods in Applied Mechanics and Engineering 403 (2023) 115718.
  * [13] A. A. Howard, M. Perego, G. E. Karniadakis, P. Stinis, Multifidelity deep operator networks, arXiv preprint arXiv:2204.09157 (2022).
  * [14] L. Lu, R. Pestourie, S. G. Johnson, G. Romano, Multifidelity deep neural operators for efficient learning of partial differential equations with application to fast inverse design of nanoscale heat transport, Physical Review Research 4 (2) (2022) 023210.
  * [15] Y. Yang, G. Kissas, P. Perdikaris, Scalable uncertainty quantification for deep operator networks using randomized priors, Computer Methods in Applied Mechanics and Engineering 399 (2022) 115399.
  * [16] C. Moya, S. Zhang, G. Lin, M. Yue, Deeponet-grid-uq: A trustworthy deep operator framework for predicting the power grid’s post-fault trajectories, Neurocomputing 535 (2023) 166–182.
  * [17] G. Lin, C. Moya, Z. Zhang, B-DeepONet: An enhanced Bayesian DeepONet for solving noisy parametric PDEs using accelerated replica exchange SGLD, Journal of Computational Physics 473 (2023) 111713.
  * [18] L. Liu, W. Cai, Multiscale DeepONet for nonlinear operators in oscillatory function spaces for building seismic wave responses, arXiv preprint arXiv:2111.04860 (2021).
  * [19] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895 (2020).
  * [20] T. Tripura, S. Chakraborty, Wavelet neural operator: a neural operator for parametric partial differential equations, arXiv preprint arXiv:2205.02191 (2022).
  * [21] V. Fanaskov, I. Oseledets, Spectral neural operators, arXiv preprint arXiv:2205.10573 (2022).
  * [22] B. Raonić, R. Molinaro, T. Rohner, S. Mishra, E. de Bezenac, Convolutional neural operators, arXiv preprint arXiv:2302.01178 (2023).
  * [23] F. Bartolucci, E. de Bézenac, B. Raonić, R. Molinaro, S. Mishra, R. Alaifari, Are neural operators really neural operators? frame theory meets operator learning, arXiv preprint arXiv:2305.19913 (2023).
  * [24] J. K. Gupta, J. Brandstetter, Towards multi-spatiotemporal-scale generalized PDE modeling, arXiv preprint arXiv:2209.15616 (2022).
  * [25] V. Oommen, K. Shukla, S. Desai, R. Dingreville, G. E. Karniadakis, Rethinking materials simulations: Blending direct numerical simulations with neural operators (2023).  arXiv:2312.05410.
  * [26] O. Ovadia, V. Oommen, A. Kahana, A. Peyvan, E. Turkel, G. E. Karniadakis, Real-time inference and extrapolation via a diffusion-inspired temporal transformer operator (DiTTO) (Dec 2023).
  * [27] M. A. Rahman, Z. E. Ross, K. Azizzadenesheli, U-NO: U-shaped neural operators, arXiv preprint arXiv:2204.11127 (2022).
  * [28] O. Ronneberger, P. Fischer, T. Brox, U-Net: Convolutional networks for biomedical image segmentation, in: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, Springer, 2015, pp. 234–241.
  * [29] 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.
  * [30] 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.
  * [31] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378 (2019) 686–707.
  * [32] A. D. Jagtap, G. E. Karniadakis, How important are activation functions in regression and classification? a survey, performance comparison, and future directions, Journal of Machine Learning for Modeling and Computing 4 (1) (2023).
  * [33] 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.
  * [34] A. D. Jagtap, K. Kawaguchi, G. 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.
  * [35] A. D. Jagtap, Z. Mao, N. Adams, G. E. Karniadakis, Physics-informed neural networks for inverse problems in supersonic flows, Journal of Computational Physics 466 (2022) 111402.
  * [36] 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.
  * [37] 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.
  * [38] A. Peyvan, K. Shukla, J. Chan, G. Karniadakis, High-order methods for hypersonic flows with strong shocks and real chemistry, Journal of Computational Physics 490 (2023) 112310.
  * [39] A. Krizhevsky, I. Sutskever, G. E. Hinton, Imagenet classification with deep convolutional neural networks, Advances in Neural Information Processing Systems 25 (2012).
  * [40] Y. Wu, K. He, Group normalization, in: Proceedings of the European conference on computer vision (ECCV), 2018, pp. 3–19.
  * [41] D. Hendrycks, K. Gimpel, Gaussian error linear units (gelus), arXiv preprint arXiv:1606.08415 (2016).
  * [42] K. Yamaguchi, K. Sakamoto, T. Akabane, Y. Fujimoto, A neural network for speaker-independent isolated word recognition., in: First International Conference on Spoken Language Processing, ICSLP, 1990, pp. 1077–1080.
  * [43] V. Dumoulin, F. Visin, A guide to convolution arithmetic for deep learning, arXiv preprint arXiv:1603.07285 (2016).
  * [44] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: a Practical Introduction, Springer Science & Business Media, 2013.