Learning Thermodynamically Stable and Galilean Invariant Partial Differential
Equations for Non-equilibrium Flows

Juntao Huang 111Department of Mathematics, Michigan State University, East
Lansing, MI 48824, USA. Zhiting Ma 222Department of Mathematical Sciences,
Tsinghua University, Beijing, China. Yizhou Zhou 22footnotemark: 2 Wen-An Yong
22footnotemark: 2, ‡‡‡E-mail<EMAIL_ADDRESS>Corresponding author

Abstract

In this work, we develop a method for learning interpretable,
thermodynamically stable and Galilean invariant partial differential equations
(PDEs) based on the Conservation-dissipation Formalism of irreversible
thermodynamics. As governing equations for non-equilibrium flows in one
dimension, the learned PDEs are parameterized by fully-connected neural
networks and satisfy the conservation-dissipation principle automatically. In
particular, they are hyperbolic balance laws and Galilean invariant. The
training data are generated from a kinetic model with smooth initial data.
Numerical results indicate that the learned PDEs can achieve good accuracy in
a wide range of Knudsen numbers. Remarkably, the learned dynamics can give
satisfactory results with randomly sampled discontinuous initial data and
Sod’s shock tube problem although it is trained only with smooth initial data.

Key Words: Machine learning; Non-equilibrium thermodynamics; Conservation-
dissipation formalism; Hyperbolic balance laws; Kinetic equation; Galilean
invariance.

##  1 Introduction

Emergence of machine learning and data science has brought new vitality to the
development of various scientific disciplines. Mostly, such disciplines are
based on certain fundamental principles which can be described with
mathematical frameworks. The frameworks usually contain problem-dependent
freedoms. Besides the traditional approaches like experiments, machine
learning is a promising alternative in fixing the freedoms. In this sense, the
booming data-based technique can display its power by a proper combination
with reliable scientific theories.

Thermodynamics is a classical theory studying relationships among various
apparently unrelated variables or parameters characterizing a thermodynamic
system. For non-equilibrium systems, it was explained in [41] that the most
natural form of the expected relations are evolution partial differential
equations (PDEs). Thus, the core task of irreversible thermodynamics is
choosing suitable variables and determining their evolution equations.
However, no well-accepted rules have been given so far and consequently there
are many different schools, such as Classical Irreversible Thermodynamics
[10], Rational Extended Thermodynamics [27], Extended Irreversible
Thermodynamics [20, 22], General Equation for Non-equilibrium
Reversible–irreversible Coupling [29, 30], Conservation-dissipation Formalism
(CDF) [47], Energetic Variational Approach [17], and so on. An exposition of
different versions of thermodynamics can be found in [28] and references cited
therein.

This project focuses on a combination of machine learning with the CDF theory
mentioned above. CDF is based on the first and second laws of thermodynamics.
It assumes that certain conservation laws are known a priori. By postulating
the existence of an entropy function and a dissipative matrix, we follow a
Conservation-dissipation Principle distilled in [40] to choose the non-
equilibrium variables and derive the corresponding evolution equations. The
PDEs thus obtained are hyperbolic balance laws (first-order PDEs) which
naturally meet four fundamental physical requirements proposed in [41]: the
observability of physical phenomena, time-irreversibility, long-time tendency
to equilibrium, and compatibility with existing classical theories. The
importance of these requirements were expounded in [41]. In CDF, the entropy
function and dissipative matrix are problem-dependent and constitute its
freedoms. These freedoms can be fixed by resorting to machine learning. Once
the freedoms are determined, we obtain interpretable and thermodynamically
stable PDEs governing the thermodynamic system under consideration.

As a start of this project, the present paper is devoted to finding the
governing macroscopic equations for one-dimensional non-equilibrium flows. To
be specific, we assume that the flow obeys the mesoscopic BGK (Bhatnagar-
Gross-Krook) equation in kinetic theory. Thus, the classical conservation laws
of mass, momentum and energy are known a priori. By adding one more (non-
equilibrium) variable, we obtain a hyperbolic system of first-order PDEs with
four unknown variables, which is Galilean invariant. Then we use fully
connected neural networks to represent the entropy function and dissipation
matrix. The training data are generated by numerically solving the BGK model
with smooth initial data. The neural network representations for the hidden
dynamics are trained by minimizing a loss function based on a discrete version
of the CDF-based PDEs acting on the training data.

Numerical results show that our CDF-based machine learning model indeed
satisfies the conservation-dissipation principle. This particularly implies
that the learned PDEs are symmetrizable hyperbolic. Moreover, the model can
achieve good accuracy in a wide range of Knudsen numbers for both randomly
sampled smooth and discontinuous initial data. In particular, its performance
on Sod’s shock tube problem [35] is better than the Euler equations in both
hydrodynamic and kinetic regimes. Remarkably, the learned dynamics can give
satisfactory results with discontinuous initial data although it is trained
only with smooth initial data.

The macroscopic governing equations obtained above can also be understood as a
moment closure system of the BGK model. However, it is important to point out
that the CDF does not necessarily rely on the kinetic theory. Here we use the
kinetic model just because we have no other means to obtain the training data.
We are expecting to have such data from experiments or other means so that the
CDF framework can exhibit its potential in a wide range of problems. In
addition, note that the CDF theory itself allows more non-equilibrium
variables. Here we choose only one non-equilibrium variable merely for the
sake of simplicity. It is possible to achieve better results by adding more
non-equilibrium variables.

As to the moment closure system, we mention the recent work [14], where the
authors introduced a framework to construct interpretable moment closure
models for kinetic problems. They firstly learned a set of generalized moments
using the auto-encoder to optimally represent the underlying velocity
distribution, and then learned the moment closure model for the generalized
moments using neural networks with the aim of best capturing the associated
dynamics of the kinetic equation. In this way, some novel moment closure
systems are obtained and give good prediction results. Nevertheless, it is not
clear whether the moment closure systems obtained in [14] are hyperbolic or
preserve the dissipation property of the kinetic equation characterized by the
H-theorem. In contrast, the hyperbolicity and dissipation property are
enforced in our CDF-based PDEs automatically.

Other recent works on learning PDEs have the following. In [24], the authors
designed a neural network with Galilean invariance to discover the Reynolds
stress anisotropy tensor in turbulence models from high-fidelity simulation
data. In [23], a neural network which preserves the rotational symmetry was
constructed to approximate the constitutive laws in non-Newtonian fluid
models. In [31, 32], by resorting to the Gaussian process regression or neural
networks named Physics-Informed Neural Networks, some classical PDEs with
unknown parameters or functions are explored. This was further applied in many
realistic problems [8, 43, 39]. In [4, 34], a data-based method called SINDy
(Sparse Identification of Nonlinear Dynamics) was developed to discover
dynamical systems from the perspective of sparse regression and compressed
sensing. In [26, 25], the authors used a symbolic neural network called PDE-
Net to uncover the analytic form of PDEs.

There are more works about learning ordinary differential equations (ODEs).
For Hamiltonian systems in classical mechanics, different neural networks such
as SympNets and HénonNet have been designed in [19, 5]. More works on learning
Hamiltonian systems can be found in [13, 37, 46] and references cited therein.
In [21], ODEs with dissipation structure are considered and the dissipation
structure was built into the neural network by jointly learning the dynamic
model and a Lyapunov function. In a very recent work [42], a generalized
Onsager principle was proposed and a so-called OnsagerNet was designed. In a
series of works [15, 44, 38, 45], deep neural networks with symmetry-
preserving property are explored to present the potential energy surface (PES)
that describes the interaction between all the nuclei in the system in the
molecular dynamics. For more details of the development of the physics-based
modeling with machine learning, we refer readers to the review paper [16].

The paper is organized as follows. The general CDF theory is briefly reviewed
in Section 2. In Section 3, we construct a simple CDF-based model for non-
equilibrium flows. Section 4 is devoted to a detailed procedure to learn the
freedoms in the model. The effectiveness of the procedure is demonstrated
through numerical results in Section 5. Some concluding remarks are given in
Section 6.

##  2 Preliminaries

In this section, we briefly review the CDF theory [47]. It is concerned with
an irreversible process which obeys some conservation laws of the form

$\partial_{t}u+\sum_{j=1}^{d}\partial_{x_{j}}f_{j}=0.$ (2.1)

Here $u=u(t,x)\in\mathbb{R}^{n}$ represents conserved variables depending on
the time and spatial coordinates $(t,x)$, $x\in\mathbb{R}^{d}$ with $d$ the
dimension of space, and $f_{j}$ is the corresponding flux along the
$x_{j}$-direction. If each $f_{j}$ is given in terms of the conserved
variables, the system (2.1) becomes closed. In this case, the system is
considered to be in local equilibrium and $u$ is also referred to as
equilibrium variables. However, very often $f_{j}$ depends on some extra
variables in addition to the conserved ones. The extra variables characterize
non-equilibrium features of the system under consideration, called non-
equilibrium or dissipative variables, and their choice is not unique.

Motivated by the Extended Thermodynamics [27, 20], in CDF we choose a set of
dissipative variables $v\in\mathbb{R}^{m}$ so that the flux $f_{j}$ in (2.1)
can be expressed as $f_{j}=f_{j}(u,v)$ and seek evolution equations of $v$ in
the form

$\partial_{t}v+\sum_{j=1}^{d}\partial_{x_{j}}g_{j}(u,v)=q(u,v).$ (2.2)

This is our constitutive equation. Here $g_{j}(u,v)$ is the corresponding flux
and $q=q(u,v)$ is the source vanishing at equilibrium. Together with the
conservation laws (2.1), the evolution of the non-equilibrium system is
governed by first-order PDEs

$\partial_{t}U+\sum_{j=1}^{d}\partial_{x_{j}}F_{j}(U)=Q(U),$ (2.3)

where

$\displaystyle U=\begin{pmatrix}u\\\\[5.69054pt] v\end{pmatrix},\qquad
F_{j}(U)=\begin{pmatrix}f_{j}(U)\\\\[5.69054pt] g_{j}(U)\end{pmatrix},\qquad
Q(U)=\begin{pmatrix}0\\\\[5.69054pt] q(U)\end{pmatrix}.$

It was observed in [40] that many classical PDEs of the form (2.3) respect the
following so-called Conservation-dissipation Principle:

  1. (i)

There is a strictly concave smooth function $\eta=\eta(U)$, called entropy
(density), such that the matrix product $\eta_{UU}F_{jU}$ is symmetric for
each $j$ and for all $U$ under consideration.

  2. (ii)

There is a positive definite matrix $M=M(U)$, called dissipation matrix, such
that the non-zero source can be written as $q(U)=M(U)\eta_{v}(U)$.

Here the subscript stands for the corresponding partial derivative, for
instance $\eta_{v}=\partial\eta/\partial v$ and
$\eta_{UU}=\partial^{2}\eta/\partial U^{2}$, and $\eta_{v}(U)$ should be
understood as a column vector. Note that the dissipation matrix is not assumed
to be symmetric, and its positive definiteness means that of the symmetric
part $(M+M^{T})/2$.

Inspired by the above observation, we try to determine the constitutive
relation (2.2) so that the resultant balance laws (2.3) also respect the
conservation-dissipation principle (see an example in the next section). This
is our Conservation-dissipation Formalism (CDF) proposed in [47].

This formalism has two freedoms: the entropy function $\eta=\eta(U)$ and the
dissipation matrix $M=M(U)$. They are both functions of the state variable
$U$. The former is strictly concave and the latter is positive definite.
Except these, no further restriction is imposed on the freedoms. Specific
expressions of $\eta(U)$ and $M(U)$ should be problem-dependent. With these
freedoms, CDF provides a tailored platform for machine learning to display its
power.

On the conservation-dissipation principle, we make the following simple
comments. Condition (i) is the well-known entropy condition for hyperbolic
conservation laws [11, 12]. It ensures that the first-order system (2.3) is
globally symmetrizable hyperbolic and thereby well-posed [9]. Condition (ii)
is a nonlinearization of the celebrated Onsager reciprocal relation for scalar
irreversible processes [10]. Together with the first condition, it can be
regarded as a stability criterion for non-equilibrium thermodynamics. As shown
in [41], this criterion guarantees that the CDF-based PDEs reflect the
observability of physical phenomena, time-irreversibility, long-time tendency
to equilibrium, and compatibility with existing classical theories. Further
comments can be founded in [41].

##  3 A CDF-based model

Consider an electrically neutral fluid without external forces. Its motion
obeys the conservation laws of mass, momentum and energy [20]:

$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho v)$ $\displaystyle=0,$ (3.1)
$\displaystyle\partial_{t}(\rho v)+\nabla\cdot(\rho v\otimes v+{\bf P})$
$\displaystyle=0,$ $\displaystyle\partial_{t}E+\nabla\cdot(Ev+q+{\bf P}v)$
$\displaystyle=0.$

Here $\rho$ is the fluid density, $v\in\mathbb{R}^{d}$ is the velocity,
$E=\rho e+\frac{1}{2}\rho|v|^{2}$ is the total energy with $e$ the specific
internal energy, ${\bf P}$ is the pressure tensor, and $q$ represents the heat
flux. This system gives the evolution of variables $(\rho,\rho v,E)$. It is
closed only if ${\bf P}$ and $q$ are specified. In [47], the CDF was used to
obtain a set of evolution equations for these variables with the two freedoms—
the corresponding entropy function and the dissipation matrix.

The goal of this paper is to fix the freedoms with machine learning. For
simplicity, we only consider the one-dimensional problem. Moreover, in order
to have reasonable training data, we turn to the one-dimensional BGK model

$\frac{\partial f}{\partial t}+\xi\frac{\partial f}{\partial
x}=\frac{1}{\varepsilon}(f_{M}-f),$ (3.2)

meaning that the fluid flow is governed by this kinetic model. Here
$f=f(x,t,\xi)$ is a distribution function with $\xi\in\mathbb{R}$ the particle
velocity, $\varepsilon$ is the dimensionless Knudsen number, and $f_{M}$ is
the Maxwellian

$f_{M}=f_{M}(\xi;\rho,v,T)=\frac{\rho}{(2\pi
T)^{1/2}}\exp\left(-\frac{(\xi-v)^{2}}{2T}\right)$

with

$\displaystyle\rho=\int_{\mathbb{R}}fd\xi,\quad\rho v=\int_{\mathbb{R}}\xi
fd\xi,$ $\displaystyle\rho T=\int_{\mathbb{R}}(\xi-v)^{2}fd\xi.$ (3.3)

It is well-known that the one-dimensional BGK model implies the conservation
laws (3.1) with the pressure tensor ${\bf P}$ becoming a scalar $p=\rho T$,
$e=T/2$,

$\displaystyle E=\frac{1}{2}\int_{\mathbb{R}}\xi^{2}fd\xi=\rho
e+\frac{1}{2}\rho v^{2},\qquad
q=\frac{1}{2}\int_{\mathbb{R}}(\xi-v)^{2}(\xi-v)fd\xi.$ (3.4)

In fact, multiplying (3.2) with $1,\xi,\frac{1}{2}\xi^{2}$ and integrating the
resultant equations over $\xi\in\mathbb{R}$ will immediately lead to the
conservation laws (3.1) (see e.g. [36]). Notice that for the one-dimensional
kinetic model, we have the following equation of state $p=2\rho e$. Thus, only
the heat flux $q$ in (3.1) is an extra variable.

To close (3.1) with CDF, we firstly recall that the system in equilibrium has
a specific entropy [36]

$s^{\textrm{(eq)}}=s^{\textrm{(eq)}}(\nu,e)=-k_{b}\nu\int_{\mathbb{R}}f_{M}\ln
f_{M}d\xi=k_{b}\left(\frac{1}{2}\ln e+\ln\nu\right)+C,$

where $k_{b}$ is the Boltzmann constant, $\nu=1/\rho$ is the specific volume,
and $C$ is a constant. Clearly, this entropy function is concave with respect
to $(e,\nu)$. In what follows, we take $k_{b}=1$ for simplicity.

Furthermore, we introduce a new dissipative variable $w$ and postulate that a
specific entropy of the form

$s=s(\nu,e,w;\varepsilon)=s^{\textrm{(eq)}}(\nu,e)+s^{\textrm{(neq)}}(w;\varepsilon)$

exists for the non-equilibrium system under consideration. Here for simplicity
we have assumed that there is no cross-terms between $w$ and $(\nu,e)$.
According to the conservation-dissipation principle, the undetermined and non-
equilibrium part $s^{\textrm{(neq)}}(w;\varepsilon)$ should be a concave
function of $w$ with the Knudsen number $\varepsilon$ as a parameter. For such
a concave function $s=s(\nu,e,w;\varepsilon)$, it is well-known that the
function $\eta$ defined as

$\eta=\eta(\rho,\rho v,E,\rho w;\varepsilon)=\rho s(\nu,e,w;\varepsilon)$

can act as the entropy in the conservation-dissipation principle. Namely, it
is strictly concave with respect to $(\rho,\rho v,E,\rho w)$.

Next we follow [47] and introduce the differential operator $D$ acting on a
function $h=h(x,t)$ as $Dh:=\partial_{t}(\rho h)+\partial_{x}(\rho
vh)=\rho(\partial_{t}h+v\partial_{x}h)$. Using the conservation laws (3.1), we
calculate

$\displaystyle\eta_{t}+\partial_{x}(v\eta)\equiv Ds$
$\displaystyle=s_{\nu}D\nu+s_{e}De+s_{w}Dw$
$\displaystyle=\rho\partial_{x}v+\theta^{-1}(-\partial_{x}q-\rho\theta\partial_{x}v)+s_{w}Dw$
$\displaystyle=-\theta^{-1}\partial_{x}q+s_{w}Dw$
$\displaystyle=-\partial_{x}(\theta^{-1}q)+s_{w}Dw+q\partial_{x}\theta^{-1}$

with $\theta^{-1}:=s_{e}$. Here we have used $s_{\nu}=\rho$ and
$s_{e}=1/(2e)$. This calculation suggests that $\theta^{-1}q$ is the entropy
flux and $s_{w}Dw+q\partial_{x}\theta^{-1}$ is the entropy production.

According to the calculation above, we choose the heat flux

$q=s_{w}=s^{\textrm{(neq)}}_{w}(w;\varepsilon)$

and the evolution equation for $w$ as

$\partial_{t}(\rho w)+\partial_{x}(\rho vw)+\partial_{x}\theta^{-1}=Mq$

with $M=M(\rho,e,w;\varepsilon)$ a positive function. Here $M$ is independent
of the velocity $v$ such that the model is Galilean invariant. By such a
choice, we have the entropy production

$s_{w}Dw+q\partial_{x}\theta^{-1}=qMq\geq 0.$

Moreover, we notice that the non-equilibrium variable $w$ can be globally
expressed in terms of $q$ and $\varepsilon$ due to the strict concavity of
$s^{\textrm{(neq)}}(w;\varepsilon)$ with respect to $w$.

Consequently, we arrive at the following balance laws

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)$ $\displaystyle=0,$ (3.5)
$\displaystyle\partial_{t}(\rho v)+\partial_{x}(\rho v^{2}+\rho\theta)$
$\displaystyle=0,$
$\displaystyle\partial_{t}E+\partial_{x}(Ev+s^{\textrm{(neq)}}_{w}+\rho\theta
v)$ $\displaystyle=0,$ $\displaystyle\partial_{t}(\rho w)+\partial_{x}(\rho
vw)+\partial_{x}\theta^{-1}$ $\displaystyle=Ms^{\textrm{(neq)}}_{w}$

with the freedoms $M=M(\rho,e,w;\varepsilon)$ and
$s^{\textrm{(neq)}}=s^{\textrm{(neq)}}(w;\varepsilon)$. Here we notice that
the temperature $T=\theta$ for the one-dimensional kinetic model. This system
satisfies the conservation-dissipation principle and thereby is globally
symmetrizable hyperbolic. Moreover, it is Galilean invariant.

##  4 Learning the freedoms

In this section, we present a procedure to learn the two undetermined
functions $s^{\textrm{(neq)}}=s^{\textrm{(neq)}}(w;\varepsilon)$ and
$M=M(\rho,e,w;\varepsilon)$ in (3.5). Note that an unlimited amount of
training data $(\rho,\rho v,E,q)$ can be obtained by numerically solving the
BGK model and then by computing the moments in (3.3)-(3.4). The details are
given in Section 5.

Although the training data $(\rho,\rho v,E,q)$ are available, the values of
$w$ cannot be obtained directly since
$s^{\textrm{(neq)}}=s^{\textrm{(neq)}}(w;\varepsilon)$ has not been fixed yet.
Therefore, we use the variables $(\rho,\rho v,E,q)$ with
$q=s^{\textrm{(neq)}}_{w}(w;\varepsilon)$ and rewrite the balance laws (3.5)
as

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)$ $\displaystyle=0,$ (4.1)
$\displaystyle\partial_{t}(\rho v)+\partial_{x}(\rho v^{2}+\rho\theta)$
$\displaystyle=0,$ $\displaystyle\partial_{t}E+\partial_{x}(Ev+q+\rho\theta
v)$ $\displaystyle=0,$
$\displaystyle\partial_{t}q+v\partial_{x}q+\frac{g}{\rho}\partial_{x}\theta^{-1}$
$\displaystyle=\frac{gMq}{\rho}$

with

$g=s^{\textrm{(neq)}}_{ww}(w;\varepsilon)<0.$

For smooth solutions, the system (4.1) is equivalent to the balance laws
(3.5). Recall that the non-equilibrium variable $w$ can be globally expressed
in terms of $q$ and $\varepsilon$ due to the strict concavity of
$s^{\textrm{(neq)}}(w;\varepsilon)$ with respect to $w$. Then $g$ can be
viewed as a function of $(q;\varepsilon)$. Therefore, our task becomes to
learn the negative $g=g(q;\varepsilon)$ and the positive
$M=M(\rho,e,q;\varepsilon)$.

Because the training data are known only at discrete space-time points, it is
a common practice to replace or approximate the PDEs with their discrete
versions [14]. Here we discretize the last equation in (4.1) as

$q_{j}^{n+1}=q_{j}^{n}-\frac{\Delta t}{2\Delta
x}v_{j}^{n}(q_{j+1}^{n}-q_{j-1}^{n})-\frac{\Delta t}{2\Delta
x}\frac{g_{j}^{n}}{\rho_{j}^{n}}\left((\theta_{j+1}^{n})^{-1}-(\theta_{j-1}^{n})^{-1}\right)+\Delta
t\left(\frac{gMq}{\rho}\right)_{j}^{n},$

where the indices $j$ and $n$ together denote the space-time point $(n\Delta
t,j\Delta x)$. By writing the last equation in the abstract form

$q_{j}^{n+1}=\mathcal{S}[g,M](V_{j-1}^{n},V_{j}^{n},V_{j+1}^{n};\Delta
t,\Delta x)$ (4.2)

with $V=(\rho,v,E,q)$, we define our loss function as the mean squared error
(MSE):

$\mathcal{L}=\sum_{\text{training
data}}|q_{j}^{n+1}-\mathcal{S}[g,M](V_{j-1}^{n},V_{j}^{n},V_{j+1}^{n};\Delta
t,\Delta x)|^{2}.$ (4.3)

Note that the first three equations in (4.1) do not involve $g$ and $M$ and
thereby they are irrelevant at this point. In addition, the equation (4.2) is
just a discrete version of the fourth equation in (4.1). It should not be
understood as a numerical scheme thereof.

With the loss function defined in (4.3), we use fully-connected neural
networks to approximate $g$ and $M$ in (4.1). To ensure the positivity of $M$
and $-g$, the softplus function is added in the output layer. Moreover, only
smooth solutions are taken as training data since it is not clear how to
define weak solutions for the non-conservative system (4.1) [9].

After having $g$ and $M$, we need another neural network to get the balance
laws (3.5) where
$s^{\textrm{(neq)}}_{w}=s^{\textrm{(neq)}}_{w}(w;\varepsilon)$ is to be
determined. To do this, we recall from
$q=s^{\textrm{(neq)}}_{w}(w;\varepsilon)$ that $w$ can be solved as a function
$w=F(q;\varepsilon)$ satisfying

$q=s^{\textrm{(neq)}}_{w}(F(q;\varepsilon);\varepsilon).$

From this we have the relation

$g=s^{\textrm{(neq)}}_{ww}(w;\varepsilon)=s^{\textrm{(neq)}}_{ww}(F(q,\varepsilon);\varepsilon)=\frac{1}{F_{q}(q;\varepsilon)}.$

Namely, $F_{q}(q;\varepsilon)$ has been trained.

Then we design a fully-connected neural network $\mathcal{F}(q;\varepsilon)$
to approximate $F(q;\varepsilon)$ which should be strictly decreasing. With
the automatic differentiation technique [3] to compute the derivative
$\mathcal{F}_{q}(q;\varepsilon)$ with respect to $q$, we define the loss
function as

$\mathcal{L}_{F}=\sum_{q}\left|\frac{1}{g(q;\varepsilon)}-\mathcal{F}_{q}(q;\varepsilon)\right|^{2}.$
(4.4)

Here the summation is over some discrete $q$ in a reasonable interval. Having
$F(q;\varepsilon)$, the function $s^{\textrm{(neq)}}_{w}(w;\varepsilon)$ can
be uniquely determined by solving the nonlinear algebraic equation
$w=F(q;\varepsilon)$ with the bisection method. In this way, the balance laws
(3.5) are fixed.

We conclude this section with the following remarks.

######  Remark 4.1.

The learned balance laws (3.5) can be used to predict non-smooth solutions
although they are learned through the equivalent version (4.1) only with
smooth training data. In the prediction, we will use the Lax-Friedrichs scheme
which is irrelevant to (4.2).

######  Remark 4.2.

Instead of resorting to the non-conservative form (4.1), we have also explored
the possibility in learning the balance laws (3.5) directly. However, the
computational cost is much larger. The reason is that the algebraic equation
$q=s^{\textrm{(neq)}}_{w}(w;\varepsilon)$ has to be solved to obtain $w$
whenever the neural network defining the equation is updated.

######  Remark 4.3.

In learning the function $F(q;\varepsilon)$, we have also tried to use the
convex neural network [1] to enforce its monotonicity with respect to $q$. But
the performance does not seem to be better.

##  5 Numerical results

This section contains three subsections, where we present some results of the
learning procedure given in the previous section. The domain of the balance
laws (3.5) is taken to be $(x,t)\in[-\pi,\pi]\times[0,0.5]$ with periodic
boundary conditions.

###  5.1 Generating Data

In this part, we show how to generate the training data. To do this, we solve
the BGK model (3.2) in the domain $(x,t)\in[-\pi,\pi]\times[0,0.5]$ with
periodic boundary conditions and initial data constructed below. Once the
distribution $f=f(x,t;\xi)$ is solved, the training data $(\rho,\rho v,E,q)$
can be obtained by using (3.3) and (3.4).

For the purpose of computation, the velocity space $\xi\in\mathbb{R}$ is
truncated into a sufficiently large interval
$\xi\in[-\left|\xi\right|_{\max},\left|\xi\right|_{\max}]$, which is
discretized into $N_{\xi}$ grid points with
$\Delta\xi=2\left|\xi\right|_{\max}/N_{\xi}$. The spatial interval
$[-\pi,\pi]$ is discretized into $N_{x}$ grid points with $\Delta
x=2\pi/N_{x}$. To ensure that the numerical error is sufficiently small, we
use high-order numerical methods: the fifth-order WENO method [18] in the
space variable and the third-order implicit-explicit (IMEX) method [2] in
time.

The time step is taken to be $\Delta t=0.1\Delta x$. We fix the parameters
$\left|\xi\right|_{\max}=10$ and $N_{\xi}=100$ such that the discretization
error in velocity is much smaller than that in space and time. The number of
grid points in space is taken to be $N_{x}=80$. Our numerical experience shows
that the numerical error is generally less than $10^{-6}$ with these
parameters for initial data below.

Following [14], we consider two types of initial conditions. The first one is
the smooth function $f_{\textrm{smooth}}$. It is a convex combination of two
Maxwellians $f_{M}(\xi;U_{1})$ and $f_{M}(\xi;U_{2})$:

$f_{\textrm{smooth}}=\alpha f_{M}(\xi;U_{1})+(1-\alpha)f_{M}(\xi;U_{2})$

with $\alpha$ sampled from $[0,1]$. Here the macroscopic variables
$U_{i}=(\rho_{i},v_{i},T_{i})$ for $i=1,2$ are taken to be the sine waves

$\left\\{\begin{aligned}
\rho_{i}(x,0)&=a_{\rho,i}\sin({k_{\rho,i}}x+\psi_{\rho,i})+b_{\rho,i},\\\
v_{i}(x,0)&=0,\\\
T_{i}(x,0)&=a_{T,i}\sin({k_{T,i}}x+\psi_{T,i})+b_{T,i}.\end{aligned}\right.$
(5.1)

In (5.1), the parameters $a_{\rho,i}$ and $a_{T,i}$ are sampled from
$[0.2,0.3]$, $k_{\rho,i}$ and $k_{T,i}$ are from $\\{1,2\\}$, $\psi_{\rho,i}$
and $\psi_{T,i}$ are from $[0,2\pi]$, and $b_{\rho,i}$ and $b_{T,i}$ from
$[0.5,0.7]$. Note that with such parameters, $\rho_{i}$ and $T_{i}$ are always
positive.

The second type of initial conditions are of the form

$f_{\textrm{shock}}=\alpha
f_{M}(\xi;U_{\textrm{smooth}})+(1-\alpha)f_{M}(\xi;U_{\textrm{shock}})$

with $\alpha$ sampled from $[0,1]$. Here $U_{\textrm{smooth}}$ is taken to be
$U_{1}$ or $U_{2}$ above, and

$U_{\textrm{shock}}(x,0)=\left\\{\begin{aligned} &(\rho_{1},0,T_{1}),\quad
x\in[-\pi,x_{1}]\cup[x_{2},\pi],\\\ &(\rho_{2},0,T_{2}),\quad
x\in(x_{1},x_{2}).\\\ \end{aligned}\right.$ (5.2)

In (5.2), $\rho_{1}$ and $T_{1}$ are two constants sampled from $[1,1.1]$;
$\rho_{2}$ and $T_{2}$ from $[0.55,0.65]$; and $x_{1}$ and $x_{2}$ are the
discontinuity locations sampled from $[-2,-1.8]$ and $[1.5,1.7]$,
respectively.

###  5.2 Training

To train the functions $g=g(q;\varepsilon)$ and $M=M(\rho,e,q;\varepsilon)$ in
(4.1), we design two fully-connected neural networks based on the loss
function (4.3). Each of them has three hidden layers with 30 nodes in each
layer. As the activation function, the softplus function is used in the output
layer to ensure the positivity of $M$ and $-g$, while the $\tanh$ function is
used in other layers. The networks are trained by using back propagation with
the stochastic gradient descent (SGD) algorithm [33]. The learning rate is
taken to be 0.05, the momentum is 0.9, and the batch size is 50.

In our training, we only use $50$ different initial values of the first type
(5.1) with $k_{\rho,i}=k_{T,i}=1$. The training data are the solutions,
initiated at these initial values, evaluated at the space-time points
$(x_{j},t^{n})=(-\pi+j\Delta x,n\Delta t)$ with $j=1,2,\cdots,80$ and
$n=0,1,\cdots,10$. Here and in the loss function (4.3) $\Delta t=0.1$ and
$\Delta x=2\pi/80$. Note that only the smooth initial data
$f_{\textrm{smooth}}$ (the first type) are used.

To train the inverse function $w=F(q;\varepsilon)$ of
$q=s^{\textrm{(neq)}}_{w}(w;\varepsilon)$ in (3.5), we use the loss function
(4.4) and design another fully-connected neural network
$\mathcal{F}(q;\varepsilon)$ of three hidden layers with 20 nodes in each
layer. For computing the derivative $\mathcal{F}_{q}(q;\varepsilon)$ in (4.4),
the automatic differentiation technique [3] is employed. In addition, we
approximate $w=F(q;\varepsilon)$ with
$\mathcal{F}(q;\varepsilon)-\mathcal{F}(0;\varepsilon)$, which ensures
$F(0;\varepsilon)=0$ and thereby the uniqueness.

The learned $w=F(q;\varepsilon)$ is illustrated in Figure 5.1. Clearly, it is
a strictly decreasing function of $q$ for each Knudsen number $\varepsilon$.
This monotonicity guarantees the strictly concavity of the entropy function
$s^{\textrm{(neq)}}=s^{\textrm{(neq)}}(w;\varepsilon)$ and thereby the
conservation-dissipation principle.

Figure 5.1: Profiles of function $w=F(q;\varepsilon)$ for different Knudsen
numbers $\varepsilon=10^{-3},10^{-2},10^{-1},1,10$. The intersection point
reflects $F(0;\varepsilon)=0$.

###  5.3 Prediction

To show the performance of the CDF-based machine learning model (3.5), we
solve it numerically and compare the numerical solutions with those computed
from the kinetic equation (3.2). To compute the solutions, we discretize the
learned balance laws (3.5) by using the Lax-Friedrichs scheme with a
relatively fine mesh. The scheme is first-order, monotone and can produce
entropy solutions in the presence of shocks.

Let $\hat{U}=(\hat{\rho},\hat{\rho}\hat{v},\hat{E})$ be the solution computed
by the scheme and $U=(\rho,\rho v,E)$ be the solution generated by the BGK
model (3.2) with the same initial data. We define their relative $L^{1}$ and
$L^{2}$ errors as

$E_{1}:=\frac{\sum_{j}\left|U_{j}-\hat{U}_{j}\right|}{\sum_{j}\left|U_{j}\right|},$

and

$E_{2}:=\sqrt{\frac{\sum_{j}(U_{j}-\hat{U}_{j})^{2}}{\sum_{j}(U_{j})^{2}}},$

where $U_{j}$ and $\hat{U}_{j}$ are the values of $U$ and $\hat{U}$ at the
$j$-th spatial grid point $x_{j}$ and the final time $t=0.5$. We note that the
prediction solution $\hat{U}$ is computed by the Lax-Friedrich scheme with a
refined mesh grid $N_{x}=400$.

The relative errors for smooth and discontinuous solutions are given in Table
5.1 and Table 5.2, respectively. The mean and the standard derivation (std)
are taken from 10 randomly different runs. From Table 5.1 and Table 5.2, we
see that the magnitude of relative errors is less than 3% for smooth solutions
and 6% for discontinuous solutions. Moreover, the errors become larger for
larger Knudsen number $\varepsilon$, which is reasonable because it is more
difficult to model non-equilibrium flows.

It is remarkable that although only the smooth training data are used to learn
the balance laws (3.5), the learned model performs well also with the
discontinuous initial data.

| relative $L^{1}$ error | relative $L^{2}$ error
---|---|---
Knudsen number $\varepsilon$ | mean | std | mean | std
$10^{-3}$ | 2.46e-03 | 4.13e-04 | 2.62e-03 | 4.32e-04 [t]
$10^{-2}$ | 1.72e-03 | 1.57e-04 | 1.89e-03 | 1.35e-04
$10^{-1}$ | 7.83e-03 | 2.62e-03 | 8.39e-03 | 2.71e-03
$10^{0}$ | 2.17e-02 | 6.68e-03 | 2.33e-02 | 7.14e-03
$10^{1}$ | 2.18e-02 | 5.37e-03 | 2.35e-02 | 5.73e-03 [b]
Table 5.1: Relative $L^{1}$ and $L^{2}$ errors for smooth solutions at $t=0.5$. The mean and the standard derivation (std) are taken from 10 randomly different runs. | relative $L^{1}$ error | relative $L^{2}$ error
---|---|---
Knudsen number $\varepsilon$ | mean | std | mean | std
$10^{-3}$ | 1.03e-02 | 1.12e-03 | 1.49e-02 | 2.24e-03 [t]
$10^{-2}$ | 6.30e-03 | 1.06e-03 | 8.60e-03 | 1.81e-03
$10^{-1}$ | 2.25e-02 | 2.65e-03 | 2.66e-02 | 3.02e-03
$10^{0}$ | 4.54e-02 | 6.32e-03 | 5.16e-02 | 6.62e-03
$10^{1}$ | 5.37e-02 | 1.06e-02 | 6.05e-02 | 1.20e-02 [b]
Table 5.2: Relative $L^{1}$ and $L^{2}$ errors for discontinuous solutions at
$t=0.5$. The mean and the standard derivation (std) are taken from 10 randomly
different runs.

Furthermore, we plot the solution profiles of density, momentum, and energy
obtained from the kinetic equation and the learned model at $t=0$ and $t=0.5$.
Figure 5.2 is for smooth initial data with $k_{\rho,i}=k_{T,i}=2$ in (5.1),
which are not used in the training data, and with Knudsen number
$\varepsilon=10$. Here and below, the ‘exact’ means that obtained from the
kinetic equation while ‘predict’ is for the learned model. It is observed that
the solution profiles of our model agree well with those of the BGK model. The
performances with other Knudsen numbers are similar and are omitted here.

[5pt](i) density at $t=0$ [5pt](ii) momentum at $t=0$ [5pt](iii) energy at
$t=0$
[5pt](i’) density at $t=0.5$ [5pt](ii’) momentum at $t=0.5$ [5pt](iii’) energy
at $t=0.5$

Figure 5.2: Solution profiles of density, momentum, and energy (from left to
right) at $t=0$ and $t=0.5$ (from top to bottom) with $\varepsilon=10$ and
with smooth initial data.

Below is the solution profiles with discontinuous initial data. Figures
5.3-5.7 are for different Knudsen numbers
$\varepsilon=10^{-3},10^{-2},10^{-1},1$ and $10$. For small Knudsen numbers
($\varepsilon=10^{-3},10^{-2},1$), our learned model agrees quite well with
the BGK model, see Figures 5.3-5.5. However, the deviation becomes large for
$\varepsilon=1$ and 10, see Figures 5.6 and 5.7.

[5pt](i) density at $t=0$ [5pt](ii) momentum at $t=0$ [5pt](iii) energy at
$t=0$
[5pt](i’) density at $t=0.5$ [5pt](ii’) momentum at $t=0.5$ [5pt](iii’) energy
at $t=0.5$

Figure 5.3: Solution profiles of density, momentum, and energy (from left to
right) at $t=0$ and $t=0.5$ (from top to bottom) with $\varepsilon=10^{-3}$
and with discontinuous initial data.

[5pt](i) density at $t=0$ [5pt](ii) momentum at $t=0$ [5pt](iii) energy at
$t=0$
[5pt](i’) density at $t=0.5$ [5pt](ii’) momentum at $t=0.5$ [5pt](iii’) energy
at $t=0.5$

Figure 5.4: Solution profiles of density, momentum, and energy (from left to
right) at $t=0$ and $t=0.5$ (from top to bottom) with $\varepsilon=10^{-2}$
and with discontinuous initial data.

[5pt](i) density at $t=0$ [5pt](ii) momentum at $t=0$ [5pt](iii) energy at
$t=0$
[5pt](i’) density at $t=0.5$ [5pt](ii’) momentum at $t=0.5$ [5pt](iii’) energy
at $t=0.5$

Figure 5.5: Solution profiles of density, momentum, and energy (from left to
right) at $t=0$ and $t=0.5$ (from top to bottom) with $\varepsilon=10^{-1}$
and with discontinuous initial data.

[5pt](i) density at $t=0$ [5pt](ii) momentum at $t=0$ [5pt](iii) energy at
$t=0$
[5pt](i’) density at $t=0.5$ [5pt](ii’) momentum at $t=0.5$ [5pt](iii’) energy
at $t=0.5$

Figure 5.6: Solution profiles of density, momentum, and energy (from left to
right) at $t=0$ and $t=0.5$ (from top to bottom) with $\varepsilon=1$ and with
discontinuous initial data.

[5pt](i) density at $t=0$ [5pt](ii) momentum at $t=0$ [5pt](iii) energy at
$t=0$
[5pt](i’) density at $t=0.5$ [5pt](ii’) momentum at $t=0.5$ [5pt](iii’) energy
at $t=0.5$

Figure 5.7: Solution profiles of density, momentum, and energy (from left to
right) at $t=0$ and $t=0.5$ (from top to bottom) with $\varepsilon=10$ and
with discontinuous initial data.

We further examine our model on Sod’s shock tube problem [35] with initial
data

$(\rho,v,p)_{\textrm{}}(x,0)=\left\\{\begin{aligned} &(1,0,1),\quad&x<0,\\\
&(0.125,0,0.1),\quad&x>0.\\\ \end{aligned}\right.$ (5.3)

With these initial data, we numerically solve our model, the Euler equations
and the kinetic equation up to $t=0.3$. In the computations of our learned
equations and the Euler equations, we use the Lax-Friedrich scheme with
$N_{x}=1600$. The solution profiles of density, momentum, and energy obtained
from the three models are shown in Figure 5.8 and Figure 5.9 with Knudsen
number $\varepsilon=10^{-2}$ and $\varepsilon=1$.

From Figure 5.8 with $\varepsilon=10^{-2}$, it is observed that the solution
profiles of our model agree quite well with those of the BGK model. In
contrast, the Euler equations generate sharp discontinuity, which deviates
from the BGK model. On the other hand, Figure 5.9 is for larger Knudsen number
$\varepsilon=1$. In this case, it is far away from equilibrium and difficult
to build a good macroscopic model. We observe that the solutions of the Euler
equations are quite different from those of the kinetic model, while our
learned model also has a slight discrepancy but can well capture the tendency
of the solution. The discrepancy may be due to the fact that only one non-
equilibrium variable $w$ is introduced in our model. The poor behavior was
also observed for the moment closure systems with few number of moments and
improved with larger number of moments [6, 7]. Thus, we expect that the
discrepancy could be reduced by adding more non-equilibrium variables under
the CDF framework and this is left for future work.

[5pt](i) density at $t=0.3$ [5pt](ii) momentum at $t=0.3$ [5pt](iii) energy at
$t=0.3$

Figure 5.8: Sod shock tube problem. Solution profiles of density, momentum,
and energy (from left to right) at $t=0.3$ with $\varepsilon=10^{-2}$.

[5pt](i) density at $t=0.3$ [5pt](ii) momentum at $t=0.3$ [5pt](iii) energy at
$t=0.3$

Figure 5.9: Sod shock tube problem. Solution profiles of density, momentum,
and energy (from left to right) at $t=0.3$ with $\varepsilon=1$.

##  6 Concluding remarks

In this work, we advocate combining machine learning with non-equilibrium
thermodynamics. Precisely, we develop a method for learning interpretable,
thermodynamically stable and Galilean invariant PDEs based on CDF. As
governing equations for one-dimensional non-equilibrium flows, the learned
PDEs are parameterized by fully-connected neural networks and satisfy the
conservation-dissipation principle automatically. In particular, they are
hyperbolic balance laws.

The training data are generated from the BGK model with smooth initial data.
Numerical results indicate that our CDF-based machine learning model achieves
good accuracy in a wide range of Knudsen numbers. It is remarkable that the
learned dynamics can give satisfactory results with randomly sampled
discontinuous initial data although it is trained only with smooth initial
data. Particularly, for the classical Sod’s shock tube problem, our model
behaves much better than the Euler equations.

This work could be improved in the following aspects: (1) Our model has only
one non-equilibrium variable which may not be the best. Better numerical
results are expected by introducing more non-equilibrium variables. (2) Only
one-dimensional problems are considered here. It is more challenging to
generalize our method to problems in multi-dimension.

## Acknowledgment

JH would like to thank Qi Tang in Los Alamos National Laboratory for helpful
discussions in training of the neural networks.

## References

  * [1] B. Amos, L. Xu, and J. Z. Kolter.  Input convex neural networks.  In International Conference on Machine Learning, pages 146–155, 2017.
  * [2] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri.  Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2-3):151–167, 1997.
  * [3] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind.  Automatic differentiation in machine learning: a survey.  The Journal of Machine Learning Research, 18(1):5595–5637, 2017\.
  * [4] S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the national academy of sciences, 113(15):3932–3937, 2016.
  * [5] J. W. Burby, Q. Tang, and R. Maulik.  Fast neural Poincaré maps for toroidal magnetic fields.  arXiv preprint arXiv:2007.04496, 2020.
  * [6] Z. Cai, Y. Fan, R. Li, and Z. Qiao.  Dimension-reduced hyperbolic moment method for the Boltzmann equation with BGK-type collision.  Communications in Computational Physics, 15(5):1368–1406, 2014\.
  * [7] Z. Cai and R. Li.  Numerical regularized moment method of arbitrary order for Boltzmann-BGK equation.  SIAM Journal on Scientific Computing, 32(5):2875–2907, 2010.
  * [8] X. Chen, L. Yang, J. Duan, and G. E. Karniadakis.  Solving inverse stochastic problems from discrete particle observations using the fokker-planck equation and physics-informed neural networks.  arXiv preprint arXiv:2008.10653, 2020.
  * [9] C. M. Dafermos.  Hyperbolic Conservation Laws in Continuum Physics, volume 3.  Springer, 2005.
  * [10] S. R. De Groot and P. Mazur.  Non-equilibrium Thermodynamics.  Courier Corporation, 2013.
  * [11] K. O. Friedrichs and P. D. Lax.  Systems of conservation equations with a convex extension.  Proceedings of the National Academy of Sciences, 68(8):1686–1688, 1971.
  * [12] S. K. Godunov.  An interesting class of quasilinear systems.  In Dokl. Acad. Nauk SSSR, volume 139, pages 521–523, 1961.
  * [13] S. Greydanus, M. Dzamba, and J. Yosinski.  Hamiltonian neural networks.  In Advances in Neural Information Processing Systems, pages 15379–15389, 2019.
  * [14] J. Han, C. Ma, Z. Ma, and W. E.  Uniformly accurate machine learning-based hydrodynamic models for kinetic equations.  Proceedings of the National Academy of Sciences, 116(44):21983–21991, 2019.
  * [15] J. Han, L. Zhang, R. Car, and W. E.  Deep potential: A general representation of a many-body potential energy surface.  arXiv preprint arXiv:1707.01478, 2017.
  * [16] J. Han, L. Zhang, and W. E.  Integrating machine learning with physics-based modeling.  arXiv preprint arXiv:2006.02619, 2020.
  * [17] Y. Hyon, D. Y. Kwak, and C. Liu.  Energetic variational approach in complex fluids: maximum dissipation principle.  Discrete & Continuous Dynamical Systems-A, 26(4):1291, 2010.
  * [18] G.-S. Jiang and C.-W. Shu.  Efficient implementation of weighted ENO schemes.  Journal of Computational Physics, 126(1):202–228, 1996.
  * [19] P. Jin, A. Zhu, G. E. Karniadakis, and Y. Tang.  Symplectic networks: Intrinsic structure-preserving networks for identifying Hamiltonian systems.  arXiv preprint arXiv:2001.03750, 2020.
  * [20] D. Jou, J. Casas-Vázquez, and G. Lebon.  Extended Irreversible Thermodynamics.  In Extended Irreversible Thermodynamics, pages 41–74. Springer, 1996.
  * [21] J. Z. Kolter and G. Manek.  Learning stable deep dynamics models.  In Advances in Neural Information Processing Systems, pages 11128–11136, 2019.
  * [22] G. Lebon, D. Jou, and J. Casas-Vázquez.  Understanding Non-equilibrium Thermodynamics, volume 295.  Springer, 2008.
  * [23] H. Lei, L. Wu, and W. E.  Machine learning based non-Newtonian fluid model with molecular fidelity.  arXiv preprint arXiv:2003.03672, 2020.
  * [24] J. Ling, A. Kurzawski, and J. Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 807:155–166, 2016.
  * [25] Z. Long, Y. Lu, and B. Dong.  PDE-Net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network.  Journal of Computational Physics, 399:108925, 2019.
  * [26] Z. Long, Y. Lu, X. Ma, and B. Dong.  PDE-Net: Learning PDEs from data.  In International Conference on Machine Learning, pages 3208–3216, 2018.
  * [27] I. Müller and T. Ruggeri.  Rational Extended Thermodynamics.  1998\.
  * [28] W. Muschik.  Contact temperature and internal variables: A glance back, 20 years later.  Journal of Non-Equilibrium Thermodynamics, 39(3):113–121, 2014\.
  * [29] H. C. Öttinger.  Beyond equilibrium thermodynamics.  John Wiley & Sons, 2005.
  * [30] M. Pavelka, V. Klika, and M. Grmela.  Multiscale Thermo-dynamics: Introduction to GENERIC.  Walter de Gruyter GmbH & Co KG, 2018.
  * [31] M. Raissi and G. E. Karniadakis.  Hidden physics models: Machine learning of nonlinear partial differential equations.  Journal of Computational Physics, 357:125–141, 2018.
  * [32] M. Raissi, P. Perdikaris, and 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:686–707, 2019.
  * [33] S. Ruder.  An overview of gradient descent optimization algorithms.  arXiv preprint arXiv:1609.04747, 2016.
  * [34] S. Rudy, A. Alla, S. L. Brunton, and J. N. Kutz.  Data-driven identification of parametric partial differential equations.  SIAM Journal on Applied Dynamical Systems, 18(2):643–660, 2019\.
  * [35] G. A. Sod.  A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.  Journal of Computational Physics, 27(1):1–31, 1978.
  * [36] H. Struchtrup.  Macroscopic transport equations for rarefied gas flows.  In Macroscopic Transport Equations for Rarefied Gas Flows, pages 145–160. Springer, 2005.
  * [37] P. Toth, D. J. Rezende, A. Jaegle, S. Racanière, A. Botev, and I. Higgins.  Hamiltonian generative networks.  arXiv preprint arXiv:1909.13789, 2019.
  * [38] H. Wang, L. Zhang, J. Han, and W. E.  Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics.  Computer Physics Communications, 228:178–184, 2018.
  * [39] M. Yin, X. Zheng, J. D. Humphrey, and G. E. Karniadakis.  Non-invasive inference of thrombus material properties with physics-informed neural networks.  arXiv preprint arXiv:2005.11380, 2020.
  * [40] W.-A. Yong.  An interesting class of partial differential equations.  Journal of Mathematical Physics, 49(3):033503, 2008.
  * [41] W.-A. Yong.  Intrinsic properties of conservation-dissipation formalism of irreversible thermodynamics.  Philosophical Transactions of the Royal Society A, 378(2170):20190177, 2020.
  * [42] H. Yu, X. Tian, W. E, and Q. Li.  OnsagerNet: Learning stable and interpretable dynamics using a generalized Onsager principle.  arXiv preprint arXiv:2009.02327, 2020.
  * [43] E. Zhang, M. Yin, and G. E. Karniadakis.  Physics-informed neural networks for nonhomogeneous material identification in elasticity imaging.  arXiv preprint arXiv:2009.04525, 2020.
  * [44] L. Zhang, J. Han, H. Wang, R. Car, and W. E.  Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics.  Physical Review Letters, 120(14):143001, 2018.
  * [45] L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, and W. E.  End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems.  In Advances in Neural Information Processing Systems, pages 4436–4446, 2018.
  * [46] Y. D. Zhong, B. Dey, and A. Chakraborty.  Symplectic ODE-net: Learning Hamiltonian dynamics with control.  arXiv preprint arXiv:1909.12077, 2019.
  * [47] Y. Zhu, L. Hong, Z. Yang, and W.-A. Yong.  Conservation-dissipation formalism of irreversible thermodynamics.  Journal of Non-Equilibrium Thermodynamics, 40(2):67–74, 2015.