# RelaxNet: A structure-preserving neural network to approximate the Boltzmann
collision operator

Tianbai Xiao<EMAIL_ADDRESS>Martin Frank  Karlsruhe Institute of
Technology, Karlsruhe, Germany

###### Abstract

The extremely high dimensionality and nonlinearity in the Boltzmann equation
bring tremendous difficulties to the study of rarefied gas dynamics. This
paper addresses a neural network-based surrogate model that provides a
structure-preserving approximation for the fivefold collision integral. The
idea and notion originate from the similarity in structure between the BGK-
type relaxation model and residual neural network (ResNet) when a particle
distribution function is treated as the input to the neural network function.
Therefore, we extend the ResNet architecture and construct what we call the
relaxation neural network (RelaxNet). Specifically, two feed-forward neural
networks with physics-informed connections and activations are introduced as
building blocks in RelaxNet, which provide bounded and physically realizable
approximations of the equilibrium distribution and velocity-dependent
relaxation time respectively. The evaluation of the collision term is thus
significantly accelerated due to the fact that the convolution in the fivefold
integral is replaced by tensor multiplication in the neural network. We fuse
the mechanical advection operator and the RelaxNet-based collision operator
into a unified model named the universal Boltzmann equation (UBE). We prove
that UBE preserves the key structural properties in a many-particle system,
i.e., positivity, conservation, invariance, H-theorem, and correct fluid
dynamic limit. These properties promise that RelaxNet is superior to
strategies that naively approximate the right-hand side of the Boltzmann
equation using a machine learning model. The implementation of RelaxNet and
training strategy are demonstrated. The construction of the RelaxNet-based UBE
and its solution algorithm are presented in detail. Several numerical
experiments, where the ground-truth datasets for the supervised learning task
are produced by the Shakhov model, velocity-dependent $\nu$-BGK model, and the
full Boltzmann equation, are investigated. The capability of the current
approach for simulating non-equilibrium flow physics is validated through
excellent in- and out-of-distribution performance. The open-source codes to
reproduce the numerical results are available under the MIT license [1].

###### keywords:

kinetic theory, computational fluid dynamics, scientific machine learning,
artificial neural network

Table 1: Nomenclature. $t$, $\mathbf{x}$, $\mathbf{v}$ | time, space, and particle velocity variables
---|---
$f$ | particle distribution function
$\mathcal{Q}$, $\mathcal{A}$ | collision and advection operators of the Boltzmann equation
$\nu$ | relaxation frequency
$\bm{\psi}$ | collision invariants
$H$ | function used in the H-theorem
$\mathcal{M}$ | Maxwellian distribution function
$\rho$, $\mathbf{V}$, $T$, $E$ | density, bulk velocity, temperature, and energy
$\mathbf{c}$ | peculiar velocity
$m$ | molecular mass
$k$ | Boltzmann constant
$\mathcal{E}$ | equilibrium state in the relaxation model
$\sigma$ | activation function
$\mathcal{E}$-net | subnet in RelaxNet to approximate equilibrium state
$\tau$-net | subnet in RelaxNet to approximate relaxation time
$D$ | dimension
$\tau$ | relaxation time
$\bar{\nu}$, $\bar{\tau}$ | mean relaxation frequency and time
$\mu$ | viscosity
$p$ | pressure
$\mathcal{M}^{+}$ | correction term for equilibrium distribution
$\bm{\alpha}$, $\bm{\beta}$ | outputs of $\mathcal{E}$-net and $\tau$-net
ELU, ReLU | Exponential Linear Unit and Rectified Linear Unit functions
$\theta$ | trainable parameters
$\mathrm{NN}_{\theta}$ | RelaxNet function
$\theta_{1}$, $\theta_{2}$ | The parts of $\theta$ belonging to $\mathcal{E}$-net and $\tau$-net
$C$ | cost function
$\lambda$, $\xi$ | regularization strengths
$\varepsilon$ | minimum value of cost function
$\hat{\bm{\alpha}}$ | Coefficients of Maxwellian in summation form of collision invariants
$\tilde{\bm{\alpha}}$ | product of ${\bm{\alpha}}$ and $\hat{\bm{\alpha}}$
$\mathbf{P}$, $\mathbf{q}$ | pressure tensor and heat flux
$\mathbf{m}$ | moment variables
$\bm{\phi}$ | basis functions
$N_{m}$ | polynomial degree of $\bm{\phi}$
$\bm{\gamma}$ | Lagrange multipliers of the entropy closure problem
$\mathscr{H}$ | Heaviside step function
$N_{i}$ | number of samples
$f_{\mathrm{ref}}$, $\mathcal{Q}_{\mathrm{ref}}$ | referenced distribution function and collision term in the dataset
$N_{h}$ | number of hidden layers
CFL | Courant–Friedrichs–Lewy number

##  1 Introduction

The advances in machine learning are increasingly incorporated in
computational physics to enable versatile data-driven modeling and simulation.
Interestingly, not so long ago, many were convinced that there are significant
differences between scientific computing and machine learning tasks. For
example, scientific computing generally processes smaller but more
computationally intensive amounts of data. These data from different locations
are expected to be linked by physical principles, whereas classical machine
learning mainly deals with discrete and localized data. However, recent
progress, e.g., physics-informed neural networks [2], suggests that there is
ample room for building data-efficient and physics-enhanced machine learning
approaches for scientific applications at the intersection of classical
numerical methods.

Among possible applications, a key use of scientific machine learning (SciML)
is to help establish reliable physical models. This approach is especially
crucial when first-principle models are missing, e.g., fitting potential-
energy hypersurfaces whose topographical features are sufficiently close to
those of real but unknown surfaces such that the simulation results from
molecular dynamics (MD) or Monte-Carlo (MC) methods are experimentally
meaningful [3, 4]. Another promising direction is to seek simplified models
for high-dimensional and nonlinear terms. An appropriate surrogate model can
significantly reduce the dimensionality and complexity of differential and
integral equations, thereby improving computational efficiency while
maintaining satisfactory accuracy [5, 6].

The SciML-integrated model is expected to be as universally accurate as
possible to the ground-truth physics. It should state fundamental physical
laws, e.g., the conservation of mass in fluid dynamic equations. It is
desirable to preserve the structural properties under physical constraints,
e.g., symmetry, invariance, and non-decreasing entropy. Different from the
brute-force solution, this requires us to design physics-oriented
architectures, parameters, and connections in the machine learning model, or
seek efficient regularization approaches [7, 8, 9, 10]. Ideally, one would
like to perform a small set of experiments under idealized situations and
obtain a model that works under more general conditions [11]. In other words,
we expect the trained model to extrapolate well for both in- and out-of-
distribution cases [12] and the generalization gap to be as small as possible
[13]. This further requires that the model should have clear physical meaning
and interpretability rather than pure black-box nature.

In this paper, we focus on the kinetic theory of gases. Lying at its core, the
Boltzmann equation describes the evolution of a many-particle system in a
statistical manner, which serves as a cornerstone of many hydrodynamic and
thermodynamic theories [14, 15]. The Boltzmann equation is an integro-
differential equation defined in a high-dimensional phase space. The extremely
high dimensionality and nonlinearity pose a notorious challenge for
theoretical analysis and numerical simulation. Starting from the introduction
of the relaxation-type model [16], tremendous efforts have been made to seek a
sufficiently accurate simplified model of the Boltzmann equation. To further
improve the accuracy in highly dissipative regions and preserve the correct
fluid dynamic limit, without changing the form of the relaxation model,
explorations have been conducted to modify the objective equilibrium function
[17, 18] and the relaxation frequency [19, 20]. Although these improved models
can provide the correct Prandtl number, the accuracy in describing strongly
non-equilibrium flows is less satisfactory [21].

The burgeoning scientific machine learning offers an alternative for
establishing simplified models of the Boltzmann equation. The existence of the
multifold integral term precludes direct use of automatic differentiation-
based techniques [2, 22]. Preliminary attempts have been made to build
simplified surrogate models for the Boltzmann collision integral. The
underlying techniques of related work involve neural differential equations
[23], reduced order modeling [24], and spectral representations [25]. Despite
the good results obtained by these methods, very limited proof has been
provided showing that the numerical method preserves the key structures, e.g.,
the H-theorem [26], of a many-particle system. Therefore, the preservation of
such physical constraints must rely on the quality of training and the
reliability of the dataset. It is notoriously known that both tasks are
difficult to be perfect, and the generalization performance of machine
learning models based on this paradigm can be questionable.

In this paper, we exploit the structural similarity of the BGK-type relaxation
model and the residual neural network (ResNet) to construct the relaxation
neural network (RelaxNet), which plays as a surrogate model for the Boltzmann
collision integral with quantitative interpretability. Specifically, two feed-
forward neural networks with physics-informed connections and bounded
activation functions are introduced as building blocks to provide physically
realizable approximations of the parameterized equilibrium distribution and
velocity-dependent relaxation time. The evaluation of collision term is thus
significantly accelerated since the computation of convolutions in the
fivefold integral is replaced by the tensor multiplication in the neural
network. The advection operator (left-hand side) of the original Boltzmann
equation and the neural network-based collision term are fused into a
trainable framework, i.e., the so-called universal Boltzmann equation (UBE)
[6, 23]. We prove that the RelaxNet-based UBE preserves key structural
properties. The positivity, invariance, correct fluid dynamic limit, and
hyperbolicity of advection operator hold independent of trainable parameters.
The conservation and H-theorem are preserved up to the training error. The
implementation of the RelaxNet-based UBE, training strategy, and solution
algorithm are presented in detail.

The paper is organized as follows. In Section 2, we introduce some fundamental
concepts in the kinetic theory of gases. Section 3 presents the idea and
architecture of RelaxNet as well as the training method. Section 4 describes
the methodology for generating datasets for training and testing. Section 5
details the numerical solution algorithm for the RelaxNet-based universal
Boltzmann equation. Section 6 contains the numerical experiments for both
spatially homogeneous and inhomogeneous cases to validate the current
methodology. The last section is the conclusion. The nomenclature of this
paper can be found in Table 1.

##  2 Kinetic Theory

The kinetic theory of gases describes the evolution of a many-particle system
in terms of the particle distribution function $f(t,\mathbf{x},\mathbf{v})$,
where $t\in\Gamma\subseteq\mathbb{R}^{+}$ is time variable,
$\mathbf{x}\in\mathcal{D}\subseteq\mathbb{R}^{3}$ is space variable, and
$\mathbf{v}\in\mathbb{R}^{3}$ is particle velocity.. For dilute monatomic gas
in the absence of external force, the Boltzmann equation can be written as

$\partial_{t}f+\mathbf{v}\cdot\nabla_{\mathbf{x}}f=\mathcal{Q}(f,f)=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\left[f\left(\mathbf{v}^{\prime}\right)f\left(\mathbf{v}_{*}^{\prime}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]\mathcal{B}(\cos\theta,g)d\mathbf{\Omega}d\mathbf{v}_{*},$
(1)

where $\\{\mathbf{v},\mathbf{v}_{*}\\}$ denotes the pre-collision velocities
of two classes of colliding particles, and
$\\{\mathbf{v}^{\prime},\mathbf{v}_{*}^{\prime}\\}$ is the corresponding post-
collision velocities. The collision kernel $\mathcal{B}(\cos\theta,g)$
measures the probability of collisions in different directions, where $\theta$
is the deflection angle and $g=|\mathbf{g}|=|\mathbf{v}-\mathbf{v}_{*}|$ is
the magnitude of relative pre-collision velocity. The solid angle
$\mathbf{\Omega}$ denotes the unit vector along the relative post-collision
velocity $\mathbf{v}^{\prime}-\mathbf{v}_{*}^{\prime}$, and the deflection
angle satisfies the relation $\theta=\mathbf{\Omega}\cdot\mathbf{g}/g$. Note
that if we define the collision frequency as

$\nu(\mathbf{v})=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}f\left(\mathbf{v}_{*}\right)\mathcal{B}\left(\cos\theta,g\right)d\mathbf{\Omega}d\mathbf{v}_{*},$
(2)

the collision term in Eq.(1) can be written as a combination of gain and loss,
i.e.,

$\displaystyle\mathcal{Q}(f,f)=\mathcal{Q}^{+}(f,f)+\mathcal{Q}^{-}(f),$ (3)
$\displaystyle\mathcal{Q}^{+}=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}f\left(\mathbf{v}^{\prime}\right)f\left(\mathbf{v}_{*}^{\prime}\right)\mathcal{B}(\cos\theta,g)d\mathbf{\Omega}d\mathbf{v}_{*},\quad\mathcal{Q}^{-}=\nu(\mathbf{v})f(\mathbf{v}).$

The Boltzmann equation possesses structural properties that reflect the
physical constraints of a many-particle system. Such structure is manifested
in the properties of the collision operator $\mathcal{Q}$ and the advection
operator $\mathcal{A}=\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}$, and it
plays a key role in design of the analytical and numerical solution algorithms
of the Boltzmann equation. These structural properties can be highlighted as
follows [27].

  1. (1)

Invariant range: There exists a set $B\subseteq[0,\infty)$ consistent with the
physical bounds of $f$ such that $\mathrm{Range}(f(t,\cdot,\cdot))\subseteq B$
if $\mathrm{Range}(f(0,\cdot,\cdot))\subseteq B$. A minimum requirement for
the invariant range is $f\geq 0$ since $f$ represents the density of
particles.

  2. (2)

Conservation: There exists the collision invariants
$\bm{\psi}=(1,\mathbf{v},\mathbf{v}^{2}/2)^{T}$ such that

$\langle\bm{\psi}\mathcal{Q}(g,g)\rangle=0,\quad\forall
g\in\mathrm{Dom}(\mathcal{Q}),$ (4)

where $\langle\cdot\rangle=\int_{\mathbb{R}^{3}}\cdot d\mathbf{v}$ denotes the
integral over particle velocity space. This property leads to the conservation
laws

$\partial_{t}\langle\bm{\psi}f\rangle+\nabla_{\mathbf{x}}\cdot\langle\mathbf{v}\bm{\psi}f\rangle=0,$
(5)

where $\mathbf{W}(t,\mathbf{x})=(\rho,\rho\mathbf{V},\rho
E)^{T}=\langle\bm{\psi}f\rangle$ are named as conservative variables in fluid
mechanics. These velocity-space integrals are also commonly referred as
moments of particle distribution function $f$.

  3. (3)

Hyperbolicity: For each fixed $\mathbf{v}$, the advection operator
$\mathcal{A}$ is hyperbolic over $(t,\mathbf{x})\in\Gamma\times\mathcal{D}$.

  4. (4)

H-theorem: There exists a differentiable function $H(f)=-f\log f$ such that

$\langle H^{\prime}(g)\mathcal{Q}(g,g)\rangle\leq 0,\quad\forall
g\in\mathrm{Dom}(\mathcal{Q})\ \mathrm{s.t.}\
\mathrm{Range}(g)\subseteq\mathcal{D},$ (6)

which further results in the dissipation law

$\partial_{t}\langle
H(f)\rangle+\nabla_{\mathbf{x}}\cdot\langle\mathbf{v}H(f)\rangle\leq 0.$ (7)

The H-theorem implies that $H$ is a Lyapunov function for the Boltzmann
equation, and the logarithm of its minimizer $\mathcal{M}$ must be a linear
combination of the collision invariants $\bm{\psi}$, which can be further
proved to have the following form, i.e., the Maxwellian distribution function
[28],

$\mathcal{M}:=\rho\left(\frac{m}{2\pi
kT}\right)^{3/2}\exp\left(-\frac{m}{2kT}\mathbf{c}^{2}\right),$ (8)

where $m$ is molecular mass, $k$ is the Boltzmann constant, $T$ is
temperature, and $\mathbf{c}=\mathbf{v}-\mathbf{V}$ is the peculiar velocity.

  5. (5)

Galilean invariance: There exists the Galilean transformation defined by

$\mathcal{G}(g)(t,\mathbf{x},\mathbf{v})=g(t,\mathbf{x}-\bm{\omega}t,\mathbf{v}-\bm{\omega}),$
(9)

where $\bm{\omega}\in\mathbb{R}^{3}$ is a translational velocity, that commute
with the advection and collision operators, i.e.,

$\displaystyle\mathcal{A}(\mathcal{G}(g))=\mathcal{G}(\mathcal{A}(g)),\quad\forall
g\in\mathrm{Dom}(\mathcal{A}),$ (10)
$\displaystyle\mathcal{Q}(\mathcal{G}(g))=\mathcal{G}(\mathcal{Q}(g)),\quad\forall
g\in\mathrm{Dom}(\mathcal{Q}).$

As a consequence, the transformed particle distribution function
$g=\mathcal{G}(f)$ still satisfies Eq.(1) if $f$ is a solution of it.

Following the structure in Eq.(3) and the H-theorem, several relaxation models
have been developed to simplify the computation of the fivefold integral on
the right-hand side of the Boltzmann equation. These simplified models take
the form

$\mathcal{Q}(f)=\nu(\mathcal{E}-f),$ (11)

where $\mathcal{E}$ denotes the target equilibrium state and $\nu$ is the
relaxation frequency. In the Bhatnagar–Gross–Krook (BGK) model [16], the
equilibrium distribution $\mathcal{E}$ is set as Maxwellian, and $\nu$ is
independent of $\mathbf{v}$, which implies that particles with different
velocities share the same probability of collision. Although the BGK model can
satisfy the structural properties shown in Section 2, the underlying
assumption becomes particularly unreasonable when the magnitude of $\nu$ in
rarefied gases is moderate, which makes the BGK model inaccurate in describing
highly dissipative flows. One direction to improve the BGK model is to
introduce $\mathbf{v}$-dependent relaxation frequency in Eq.(11) [19, 20].
Following the principle that high-speed particles should have a greater
probability of collision, these models attempt to recover the true collision
frequency in Eq.(2). Another idea is to supplement correction terms in
$\mathcal{E}$ based on high-order moments of $f$, e.g., stress tensor [17] and
heat flux [18], which builds a multi-level relaxation mechanism towards
Maxwellian. It is shown that both strategies above can provide the correct
Prandtl number. However, no model has yet provided as satisfactory the
solution as the full Boltzmann equation in describing the evolution of gases
far away from equilibrium.

##  3 Relaxation Neural Network

###  3.1 Deep residual learning

Our idea is to build a surrogate model of the nonlinear function described by
the Boltzmann equation based on a neural network. To preserve as many
structural properties as shown in Section 2, we prefer a physics-oriented
architecture rather than a naive feed-forward neural network.

We choose to construct the model on the basis of residual neural network
(ResNet). To clarify this, we briefly introduce here the basic principles of
ResNets. The concept of residual learning was proposed to mitigate the
degradation problem in deep neural networks, i.e., as the network depth
increases, the accuracy gets saturated and then degrades rapidly. This
phenomenon implies that it is non-trivial to learn a deeper model based on its
shallower counterpart in which the added layers only perform identity mapping.
Therefore, instead of building each few stacked layers to fit a desired
underlying mapping directly, one can explicitly let these layers fit a
residual mapping. It has been shown that it is easier to optimize the residual
mapping than the original unreferenced mapping [29].

Figure 1 shows the schematic diagram of the building block in a ResNet. We
denote $\mathcal{H}(\mathbf{u})$ as the underlying mapping to be fit by a few
stacked layers and $\mathbf{u}$ as the input to the first layer. If a
continuous function can be approximated asymptotically by multiple nonlinear
layers, then we can equivalently conclude that they can asymptotically
approximate the residual function, i.e.,
$\mathcal{F}(\mathbf{u})=\mathcal{H}(\mathbf{u})-\mathbf{u}$, provided that
the input and output are of the same dimensions. The original function thus
becomes $\mathcal{H}(\mathbf{u})=\mathcal{F}(\mathbf{u})+\mathbf{u}$, and its
formulation is realized by adding shortcut connections [30] to feedforward or
convolutional neural networks. The shortcut can skip one or several hidden
layers and performs identity mapping. The nonlinear activation function inside
the hidden layers is denoted as $\sigma_{h}$. The addition operation can be
regarded as the connection function which combines the contributions from
$\mathcal{F}$ and identity mapping of $\mathbf{u}$. The residual block is
further activated with $\sigma$, and the final output can be written as

$\mathcal{N}(\mathbf{u})=\sigma(\mathcal{F}(\mathbf{u})+\mathbf{u}).$ (12)

Figure 1: Schematic diagram of a building block in the residual neural network
(ResNet).

It is noticeable that the relaxation model of the Boltzmann equation in
Eq.(11) has a similar structure to the ResNet. Specifically, if the particle
distribution function is considered as the input to the solver of relaxation
term, then the equilibrium distribution function $\mathcal{E}=\mathcal{E}(f)$
plays an equivalent role to that of $\mathcal{H}$. Next, unlike in neural
networks where addition is employed as the connection function, subtraction is
used in the relaxation model. Finally, the relaxation frequency can be viewed
as an activation function acting on the equilibrium function $\mathcal{E}$ and
the shortcut connection from $f$. The building block of the solver of
relaxation term can thus be written as

$\mathcal{Q}(f)=\nu(\mathcal{E}(f)-f).$ (13)

The structural similarity between the relaxation model of Boltzmann equation
in Eq.(13) and ResNet in Eq.(12) can be easily noticed. This provide us an
opportunity to build a physics-oriented neural network based on the relaxation
model and preserve the structural properties of the Boltzmann equation.

###  3.2 Relaxation neural network

Following the spirit of the relaxation model and ResNet, we build a novel
neural network model which is named as relaxation neural network (RelaxNet).
The neural network is built on top of the discrete velocity formulation of the
Boltzmann equation, i.e., the particle distribution function is discretized in
the velocity space and the macroscopic variables as its moments can evaluated
by numerical quadrature. The schematic diagram of the model with forward pass
is shown in Figure 2. As can be seen, RelaxNet consists of two feedforward
subnets and the corresponding activation and connection functions. We describe
these components in detail below.

Figure 2: Schematic diagram of forward pass in the relaxation neural network
(RelaxNet).

  1. (1)

Neural networks

     * (a)

$\mathcal{E}$-net: feedforward network
$\mathcal{N}_{\mathcal{E}}:\mathbb{R}^{N_{v}}\rightarrow\mathbb{R}^{D+2}$,
where $N_{v}$ is the number of quadrature points to discretize $f$ in the
velocity space and $D\in\\{1,2,3\\}$ is the dimension of interest. The neural
network takes the difference between Maxwellian and distribution function,
i.e., $\mathcal{M}-f$, as input. The output of neural network is a vector
$\bm{\alpha}=(\alpha_{1},\bm{\alpha}_{2},\alpha_{3})^{T}$, where
$\bm{\alpha}_{2}$ is a vector of dimension $D$. The hyperbolic tangent (tanh)
function is employed as the activation function for the hidden layers. We
require the $\mathcal{E}$-net to have only weights without biases.

     * (b)

$\tau$-net: feedforward network
$\mathcal{N}_{\tau}:\mathbb{R}^{N_{v}+1}\rightarrow\mathbb{R}^{N_{v}}$. The
input to the neural network $(\mathcal{M}-f,\bar{\tau})^{T}$ is a vector of
dimension $N_{v}+1$. Here $\bar{\tau}$ is the mean relaxation time calculated
by

$\bar{\tau}=\frac{\mu}{p},$ (14)

where $\mu$ is the viscosity that can be decided from a molecular model, and
$p$ is the pressure. The relationship between mean relaxation time and
frequency is $\bar{\nu}\bar{\tau}=1$. The output of neural network
$\bm{\beta}$ is a vector of the same dimension as $f$. The hyperbolic tangent
(tanh) function is used as the activation function for the hidden layers. We
require the $\tau$-net to have only weights without biases.

  2. (2)

Structure functions

     * (a)

$\mathcal{M}(f)$: Maxwellian function described in Eq.(8), where the
macroscopic variables can be evaluated as

$\rho=\langle
f\rangle,\quad\mathbf{V}=\frac{\langle\mathbf{v}f\rangle}{\langle
f\rangle},\quad T=\frac{m}{3k\rho}\langle\mathbf{c}^{2}f\rangle.$ (15)

     * (b)

$\mathcal{M}^{+}(\bm{\alpha})$: Constructor of the correction term for
equilibrium distribution via $\bm{\alpha}$, i.e.,

$\mathcal{M}^{+}=\exp\left(\alpha_{1}+\bm{\alpha}_{2}\cdot\mathbf{v}+\frac{1}{2}\alpha_{3}\mathbf{v}^{2}\right),$
(16)

  3. (3)

Connection and activation functions

     * (a)

ELU: Exponential Linear Unit function defined as

$R(z)=\left\\{\begin{array}[]{cc}z,&z>0,\\\ e^{z}-1,&z<=0.\\\
\end{array}\right.$ (17)

     * (b)

ReLU: Rectified Linear Unit function defined as

$R(z)=\left\\{\begin{array}[]{cc}z,&z>0,\\\ 0,&z<=0.\\\ \end{array}\right.$
(18)

     * (c)

$\oplus$: Addition.

     * (d)

$\ominus$: Subtraction.

     * (e)

$\otimes$: Multiplication.

Based on the proposed RelaxNet, a neural network enhanced Boltzmann equation
can be formulated as

$\partial_{t}f+\mathbf{v}\cdot\nabla_{\mathbf{x}}f=\mathrm{NN}_{\theta}(f,\bar{\tau}),$
(19)

where $\theta$ denotes trainable parameters in the model. Such an equation is
consistent of the idea of universal differential equation [6], where the
mechanical advection operator
$\mathcal{A}=\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}$ and the neural
collision operator $\mathrm{NN}_{\theta}$ constitute a differentiable
surrogate model. Therefore, we can call Eq.(19) as the universal Boltzmann
equation (UBE).

###  3.3 Training strategy

UBE in Eq.(19) is a data-driven model and it forms a typical supervised
training task. Given a dataset consisting of reference solutions to the
Boltzmann equation, an optimization problem can be set up which minimizes the
difference between the current prediction of UBE model and ground-truth
solution. We define the cost function as

$C(\theta)=\sum_{n,i,j}\Arrowvert\mathrm{NN}_{\theta}-\mathcal{Q}_{\mathrm{ref}}\Arrowvert(t^{n},\mathbf{x}_{i},\mathbf{v}_{j})+\lambda\sum_{n,i}\Arrowvert\langle\bm{\psi}\mathrm{NN}_{\theta}\rangle\Arrowvert(t^{n},\mathbf{x}_{i})+\xi\sum_{l}^{L}\Arrowvert\theta^{l}\Arrowvert,$
(20)

where $\Arrowvert\cdot\Arrowvert$ denotes the Euclidean distance, and
$\\{n,i,j\\}$ are the indices of the solution collocation points at
$\\{t^{n},\mathbf{x}_{i},\mathbf{v}_{j}\\}$. The generation of datasets used
for the above function will be detailed in Section 4.

The latter two terms in Eq.(20) serve as regularization terms to mitigate
overfitting. The second term in Eq.(20) comes naturally from the physical
constraint (2) shown in Section 2, i.e., the collision term should conserve
mass, momentum and energy in the Boltzmann equation. This physics-informed
regularization [31] leads to more accurate preservation of conservation laws,
and improves the physical interpretability of the model. The last term in
Eq.(20) sums over the squared weight parameters of the network, where
$\theta^{l}$ denotes the weight parameters of $l$-th layers and $L$ is the
total number of layers. The parameters $\\{\lambda,\xi\\}$ are the
regularization strengths. We require

$\lambda N_{t}N_{x}>1,\quad\xi N_{t}N_{x}N_{v}>1,$ (21)

where $\\{N_{t},N_{x},N_{v}\\}$ are the number of time steps, physical cells,
and velocity points respectively in the dataset.

The cost function in Eq.(20) can then be minimized by different optimization
algorithms. In this paper, we adopt the Adam [32] method which is an
improvement of the stochastic gradient descent method with adaptive moment
estimation. The gradient of Eq.(20) with respect to $\theta$ is computed with
the help of reverse-mode automatic differentiation (AD). Consider a smooth
function $y=\mathscr{F}(x)$, the reverse-mode AD computes the dual (conjugate-
transpose) matrix of Jacobian $\mathcal{J}=\nabla\mathscr{F}$ at $x=x_{0}$
with the chain rule

$\left(\mathcal{J}(\mathscr{F})\left(x_{0}\right)\right)^{*}=\left(\mathcal{J}\left(G_{1}\right)\left(x_{0}\right)\right)^{*}\times\cdots\times\left(\mathcal{J}\left(G_{k}\right)\left(x_{k-1}\right)\right)^{*},$
(22)

with $x_{i}:=G_{i}\left(x_{i-1}\right)\text{ for }i=1,\ldots,k-1$. The
implementation of the reverse-mode AD can be found in the open-source library
Zygote [33].

###  3.4 Structure of RelaxNet

The universal Boltzmann equation based on RelaxNet preserve the structural
properties illustrated in Section 2. We prove these properties in detail
below.

######  Theorem 1.

Let $f$ be the solution of the RelaxNet-based universal Boltzmann equation
given in Eq.(19), then the range of $f$ satisfy

$\mathrm{Range}(f(t,\mathbf{x},\mathbf{v}))\subseteq B=\left\\{f:0\leq
f<\infty\right\\},$ (23)

given $\mathrm{Range}(f(0,\mathbf{x},\mathbf{v}))\subseteq B$.

###### Proof.

We define shifted $f$ along characteristics as

$f^{\\#}(t,\mathbf{x},\mathbf{v})=f(t,\mathbf{x}+\mathbf{v}t,\mathbf{v}).$
(24)

Then UBE in Eq.(19) can be written as

$\frac{d}{dt}f^{\\#}+\nu f^{\\#}=\nu\mathcal{E}^{\\#},$ (25)

where $d/dt$ denotes the full derivative. Given the initial value

$f(0,\mathbf{x},\mathbf{v})=f_{0}(\mathbf{x},\mathbf{v}),$ (26)

then the solution of the initial value problem consisting of Eq.(25) and (26)
can be written as [34]

$f^{\\#}(t,\mathbf{x},\mathbf{v})=f_{0}(\mathbf{x},\mathbf{v})e^{-\nu
t}+\int_{0}^{t}e^{-\nu(t-s)}\nu\mathcal{E}^{\\#}(s,\mathbf{x},\mathbf{v})ds.$
(27)

Proving Theorem 1 is equivalent to proving that
$\mathrm{Range}(f^{\\#})\subseteq B$. As shown in Section 3.2, $\mathcal{E}$
and $\nu$ can be expressed as

$\displaystyle\mathcal{E}=\mathrm{ReLU}(\mathcal{M}\mathcal{M}^{+}),\quad\mathcal{M}^{+}=\exp(\bm{\alpha}_{\theta_{1}}\cdot\bm{\psi}),$
(28)
$\displaystyle\nu=\frac{1}{\tau},\quad\tau=\bar{\tau}(1+\mathrm{ELU}(\bm{\beta}_{\theta_{2}})),$

where $\theta=(\theta_{1},\theta_{2})^{T}$ denotes all the trainable
parameters. With the help of the physics-informed activation functions, i.e.,
ELU and ReLU, it is clear that $\mathcal{E}$ and $\nu$ in Eq.(28) are non-
negative independent of $\theta$. Besides, since the hidden layers inside
$\mathcal{E}$-net and $\tau$-net are activated by the tanh function, the
following relation holds,

$\mathrm{Range}(\mathcal{E})\subseteq B,\quad\mathrm{Range}(\nu)\subseteq B$
(29)

given any finite layers and parameters. So far, Theorem 1 is proved. ∎

######  Theorem 2.

Let $\mathrm{NN}_{\theta}$ be the optimized RelaxNet, then the moments of
$\mathrm{NN}_{\theta}$ with respect to the collision variants $\bm{\psi}$
satisfy

$\Arrowvert\langle\bm{\psi}\mathrm{NN}_{\theta}(g,\bar{\tau})\rangle\Arrowvert\leq\varepsilon,\quad\forall(g,\bar{\tau})\in\mathrm{Dom}(\mathrm{NN}_{\theta}),$
(30)

where $\varepsilon$ is the minimum value of the cost function.

###### Proof.

The optimized cost function in Eq.(20) satisfy

$C(\theta)=\sum_{n,i,j}\Arrowvert\mathrm{NN}_{\theta}-\mathcal{Q}\Arrowvert(t^{n},\mathbf{x}_{i},\mathbf{v}_{j})+\lambda\sum_{n,i}\Arrowvert\langle\bm{\psi}\mathrm{NN}_{\theta}\rangle\Arrowvert(t^{n},\mathbf{x}_{i})+\xi\sum_{l}^{L}\Arrowvert\theta^{l}\Arrowvert=\varepsilon,$
(31)

which results in

$\lambda\sum_{n,i}\Arrowvert\langle\bm{\psi}\mathrm{NN}_{\theta}\rangle\Arrowvert(t^{n},\mathbf{x}_{i})\leq\varepsilon.$
(32)

Due to the requirement in Eq.(21), we have

$\Arrowvert\langle\bm{\psi}\mathrm{NN}_{\theta}\rangle\Arrowvert\leq\varepsilon,$
(33)

at any solution point in the dataset. ∎

######  Theorem 3.

For each fixed $\mathbf{v}$, the advection operator $\mathcal{A}$ in the
RelaxNet-based universal Boltzmann equation is hyperbolic over
$(t,\mathbf{x})\in\Gamma\times\mathcal{D}$.

###### Proof.

As RelaxNet surrogates only the right-hand side of the Boltzmann equation,
this theorem naturally holds according to Property (3) in Section 2. ∎

######  Theorem 4.

The function $H(f)=-f\log f$ satisty

$\partial_{t}\langle
H(f)\rangle+\nabla_{\mathbf{x}}\cdot\langle\mathbf{v}H(f)\rangle\leq\sum_{i}\chi_{i}\varepsilon,$
(34)

where $\bm{\chi}=\sum_{i}\chi_{i}$ is a vector of the same dimension of
$\bm{\alpha}$. The H-theorem holds asymptotically as $\varepsilon\rightarrow
0$.

###### Proof.

Substituting $H^{\prime}(f)=-(1+\log f)$ into the RelaxNet-based UBE in
Eq.(19) and considering Eq.(28) yield

$\displaystyle\langle(\partial_{t}f+\mathbf{v}\cdot\nabla_{\mathbf{x}}f)(1+\log
f)\rangle=\partial_{t}\langle
H(f)\rangle+\nabla_{\mathbf{x}}\cdot\langle\mathbf{v}H(f)\rangle=\langle\nu(\mathcal{E}-f)(1+\log
f)\rangle.$ (35)

Notice that the Maxwellian distribution in Eq.(8) can be written as

$\mathcal{M}=\exp(\hat{\bm{\alpha}}\cdot\bm{\psi}),$ (36)

where

$\displaystyle\hat{\bm{\alpha}}=(\hat{\alpha}_{1},\hat{\bm{\alpha}}_{2},\hat{\alpha}_{3})^{T},$
(37) $\displaystyle\hat{\alpha}_{1}=\log\left(\rho\left(\frac{m}{2\pi
kT}\right)^{D/2}-\frac{m\mathbf{V}^{2}}{2kT}\right),\quad\hat{\bm{\alpha}}_{2}=\frac{m\mathbf{V}}{kT},\quad\hat{\alpha}_{3}=-\frac{m}{kT}.$

Given the structure function $\mathcal{M}^{+}$ defined in Eq.(16), the
equilibrium state can be written as

$\mathcal{E}=\mathcal{M}\mathcal{M}^{+}=\exp(\tilde{\bm{\alpha}}\cdot\bm{\psi})$
(38)

where $\tilde{\alpha}_{i}=\alpha_{i}\hat{\alpha}_{i}$. From Theorem 2, we
obtain

$\langle\nu(\mathcal{E}-f)\log\mathcal{E}\rangle=\left\langle\nu(\mathcal{E}-f)\left(\tilde{\alpha}_{1}+\tilde{\bm{\alpha}}_{2}\cdot\mathbf{v}+\frac{1}{2}\tilde{\alpha}_{3}\mathbf{v}^{2}\right)\right\rangle\leq\sum_{i}\tilde{\alpha}_{i}\varepsilon.$
(39)

The right-hand side of Eq.(35) can be written as

$\langle\nu(\mathcal{E}-f)(1+\log
f)\rangle=\langle\nu(\mathcal{E}-f)\rangle+\langle\nu(\mathcal{E}-f)\log
f\rangle\leq\varepsilon+\langle\nu(\mathcal{E}-f)\log f\rangle.$ (40)

We rewrite the last term of the above equation as

$\langle\nu(\mathcal{E}-f)\log f\rangle=\langle\nu(\mathcal{E}-f)(\log
f-\log\mathcal{E})\rangle+\langle\nu(\mathcal{E}-f)\log\mathcal{E}\rangle.$
(41)

Since the logarithm is a monotonically increasing function, the first in the
above equation is non-positive. Considering Eq.(35), (39), and (41), we end up
with

$\partial_{t}\langle
H(f)\rangle+\nabla_{\mathbf{x}}\cdot\langle\mathbf{v}H(f)\rangle\leq\left(1+\sum_{i}\tilde{\alpha}_{i}\right)\varepsilon.$
(42)

So far, Theorem 4 is proved. ∎

######  Theorem 5.

The RelaxNet-based universal Boltzmann equation is Galilean invariant.

###### Proof.

We assume $f$ is a solution to UBE in Eq.(19) and consider the Galilean
transformation defined in Eq.(9), i.e.,

$g(t,\mathbf{x},\mathbf{v})=\mathcal{G}(f)=f(t,\mathbf{x}-\bm{\omega}t,\mathbf{v}-\bm{\omega}).$
(43)

As the advection operator in UBE is the same as the original Boltzmann
equation, the commute naturally holds, i.e.,

$\mathcal{A}(\mathcal{G}(f))=\mathcal{G}(\mathcal{A}(f)).$ (44)

Besides, the inputs to RelaxNet, i.e., $\\{f,\bar{\tau}\\}$, are scalar and
Galilean invariant, and thus we have

$\mathcal{Q}(\mathcal{G}(f),\bar{\tau})=\mathcal{G}(\mathcal{Q}(f,\bar{\tau})).$
(45)

Therefore, $g$ is also a soslution to Eq.(19). ∎

######  Theorem 6.

The RelaxNet-based universal Boltzmann equation preserves the correct
continuum limit.

###### Proof.

Based on Theorem 2, we obtain the transport equations for conservative
variables

$\partial_{t}\langle\bm{\psi}f\rangle+\nabla_{\mathbf{x}}\cdot\langle\mathbf{v}\bm{\psi}f\rangle\leq\varepsilon,$
(46)

or specifically,

$\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\\ \rho\mathbf{V}\\\
\rho
E\end{array}\right)+\nabla_{\mathbf{x}}\cdot\left(\begin{array}[]{c}\rho\mathbf{V}\\\
\rho\mathbf{V}\mathbf{V}\\\ \rho
E\mathbf{V}\end{array}\right)\leq-\nabla_{\mathbf{x}}\cdot\left(\begin{array}[]{c}0\\\
\mathbf{P}\\\
\mathbf{P}\cdot\mathbf{V}+\mathbf{q}\end{array}\right)+\left(\begin{array}[]{c}\varepsilon\\\
\varepsilon\\\ \varepsilon\end{array}\right),$ (47)

where

$\mathbf{P}=\langle\mathbf{c}\mathbf{c}f\rangle,\quad\mathbf{q}=\frac{1}{2}\langle\mathbf{c}\mathbf{c}^{2}f\rangle,$
(48)

denote the pressure tensor and heat flux respectively. UBE in Eq.(19) can be
rewritten as

$f=\mathcal{E}-\tau\mathrm{NN}_{\theta}(f,\bar{\tau}),$ (49)

where $\tau=1/\nu$ is the relaxation time. In the continuum limit, infinitely
many intermolecular interactions happen in a unit of time, i.e.,
$\nu\rightarrow\infty$. From Eq.(28), we have $\tau\rightarrow 0$ and thus
$f\rightarrow\mathcal{E}$. We next prove that $\mathcal{E}$ converges to
$\mathcal{M}$ in the continuum limit. Note that the difference between
Maxwellian and distribution function, i.e., $\mathcal{M}-f$, are used as input
to the zero-bias $\mathcal{E}$-net in RelaxNet. As a result, the
$\mathcal{E}$-net produces $\bm{\alpha}\rightarrow\mathbf{0}$ and
$\mathcal{M}^{+}\rightarrow 1$, and thus we have

$\mathcal{E}=\mathcal{M}\mathcal{M}^{+}\sim\mathcal{M}\quad(\mathrm{as}\
f\rightarrow\mathcal{M}).$ (50)

In addition, given $\langle\bm{\psi}(\mathcal{M}-f)\rangle=\mathbf{0}$, the
residuals in Eq.(46) converge to zero,

$\varepsilon\rightarrow 0\quad(\mathrm{as}\ f\rightarrow\mathcal{M}),$ (51)

and we recover the exact Euler equations independent of $\theta$, i.e.,

$\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\\ \rho\mathbf{V}\\\
\rho
E\end{array}\right)+\nabla_{\mathbf{x}}\cdot\left(\begin{array}[]{c}\rho\mathbf{V}\\\
\rho\mathbf{V}\mathbf{V}\\\ \rho
E\mathbf{V}\end{array}\right)=-\nabla_{\mathbf{x}}\cdot\left(\begin{array}[]{c}0\\\
p\mathbf{I}\\\ p\mathbf{V}\end{array}\right),$ (52)

given the following relationship as $f=\mathcal{M}$,

$\mathbf{P}=p\mathbf{I},\quad\mathbf{q}=\mathbf{0}.$ (53)

∎

##  4 Data Sampling

As presented in Section 3.3, a dataset with ground-truth distribution
functions $f_{\mathrm{ref}}$ and collision terms $\mathcal{Q}_{\mathrm{ref}}$
is needed in Eq.(20) to perform the supervised learning task. We consider the
space of $f$ as a sample space under the constraint,

$f(t,\mathbf{x},\mathbf{v})\in B=\left\\{f:0\leq f<\infty\right\\},$ (54)

where $B$ can be referred as the physically realizable set. A data-
distribution $p_{f}$ can be incorporated to generate $f$ from $B$.

A common strategy for data-driven modeling is to first perform deterministic
simulations and extract solutions in the intermediate steps as sample data.
One disadvantage of this approach is that $p_{f}$ is implicitly determined by
the chosen test cases and can thus be heavily biased. It may not be possible
to cover enough different particle distributions and flow regimes. The high
expense of CFD simulations makes it challenging to establish an all-round
database.

In this paper, we provide an alternative sampling strategy to generate
particle distribution functions and corresponding collision terms. The
theoretical basis of this approach is the closure hierarchies of moment system
of the Boltzmann equation [35]. In the following we briefly explain the
principle and implementation of the approach.

###  4.1 Generation of distribution function

A closure strategy of the Boltzmann equation aims to reconstruct the particle
distribution function $f$ from moments

$\mathbf{m}=\langle\bm{\phi}f\rangle,$ (55)

under the constraint

$f(t,\mathbf{x},\mathbf{v})\in B_{M}=\left\\{f:f\geq 0,\
\Arrowvert\langle\bm{\phi}f\rangle\Arrowvert<\infty\right\\},$ (56)

where $\bm{\phi}(\mathbf{v})\in\mathbb{R}^{N_{m}}$ is a set of basis
functions. Here we choose the basis in such a way that the first moments in
$\mathbf{m}$ are consistent with the conservative variables
$\mathbf{W}=(\rho,\rho\mathbf{V},\rho E)^{T}$. Therefore, we can write
$\bm{\phi}$ as

$\bm{\phi}(\mathbf{v})=(\bm{\psi},\tilde{\bm{\phi}}(\mathbf{v}))^{T}$ (57)

where $\tilde{\bm{\phi}}(\mathbf{v})$ are chosen as monomials and mixed
polynomials of $\mathbf{v}$ up to degree $N_{m}$.

The entropy closure seeks for a unique map $\mathbf{m}\mapsto f$ by forming an
optimization problem. The objective function of the optimization problem is
built with the help of a convex entropy function. For classical Maxwell-
Statistics statistics, the minimal entropy closure problem writes

$\min_{g\in B_{m}}\langle g\log(g)-g\rangle,\quad\text{ s.t.
}\mathbf{m}=\langle\bm{\phi}g\rangle.$ (58)

The solution of the above optimization problem can be expressed as

$f_{u}(\mathbf{v})=\exp(\bm{\gamma}_{u}\cdot\bm{\phi}),$ (59)

where $\bm{\gamma}_{u}\in\mathbb{R}^{N_{m}}$ is a vector of Lagrange
multipliers of the dual formulation of the optimization problem, which reads

$\bm{\gamma}_{u}=\underset{\bm{\gamma}\in\mathbb{R}^{N_{m}}}{\operatorname{argmax}}\\{\bm{\gamma}\cdot\mathbf{m}-\langle\exp(\bm{\gamma}\cdot\bm{\phi})\rangle\\}.$
(60)

The idea is to sample particle distribution functions from Eq.(59). Note that
the condition number of the minimal entropy closure at a moment $\mathbf{u}$
can be computed via the positive semi-definite Hessian of the dual problem,
i.e.,

$\mathbf{H}_{u}=\langle\bm{\phi}\bm{\phi}\exp(\bm{\gamma}_{u}\cdot\bm{\phi})\rangle.$
(61)

where $\bm{\phi}\bm{\phi}$ denotes the union of $\bm{\phi}$. As analyzed in
[36, 37], generally the lower the condition number, the less the reconstructed
particle distribution differs from the Maxwellian. Therefore, by controlling
the condition number of $\mathbf{H}_{u}$, we can generate distribution
functions that are close and far from equilibrium. We let the first three
elements of $\bm{\gamma}$ be consistent with $\hat{\bm{\alpha}}$ in Eq.(37),
and sample $\gamma_{n}$ for $n>3$ with a prescribed probability density.

It is worthwhile to mention that the condition number cannot be arbitrarily
large due to the constraint of realizability. The realizable set of the
minimal entropy problem in Eq.(58) writes

$\mathcal{R}=\\{\mathbf{m}:\mathbf{m}=\langle\bm{\phi}g\rangle,g\in B_{m}\\}.$
(62)

It is known that the minimal entropy problem near the boundary of
$\mathcal{R}$, denoted as $\partial\mathcal{R}$, may become ill-conditioned
[38]. Therefore, we keep the condition number within a certain range, and
employ an upwind reconstruction to enable the capability to generate highly
non-equilibrium distribution functions that may exist in highly dissipative
flows. The distribution functions produced from Eq.(59) are divided randomly
into two parts, and the referenced $f$ is reconstructed as

$f_{\mathrm{ref}}=f_{L}\mathscr{H}(\mathbf{n}\cdot\mathbf{v})+f_{R}(1-\mathscr{H}(\mathbf{n}\cdot\mathbf{v})),$
(63)

where $\\{f_{L},f_{R}\\}$ are the distribution functions from these two parts,
$\mathbf{n}$ is a randomly generated unit vector, and $\mathscr{H}(x)$ is the
heaviside step function. Figure 3 shows typical near-equilibrium and non-
equilibrium particle distribution functions reconstructed from Eq.(63) using
the degree $N_{m}=4$. As the polynomial degree increases, the non-equilibrium
effect is expected to be stronger. The detailed sampling strategy is
summarized in Algorithm 1. Note that the sampling strategy here can work
together with classical simulation-based data samplers, where the data
produced by a specific simulation can be employed as reinforcement data on top
of the pre-generated data. We will explain this further in Section 5.

(a) Near-equilibrium reconstructed distribution function (b) Non-equilibrium
reconstruction distribution function

Figure 3: Typical one-dimensional near-equilibrium and non-equilibrium
particle distribution functions generated by the sampling strategy in Section
4.1 with the polynomial degree $N_{m}=4$. The values of velocity $u$ and
distribution function $f$ are normalized by $(2RT_{0})^{1/2}$ and
$\rho_{0}/(2RT_{0})^{1/2}$, respectively.

Input : Number of dataset entries $N_{i}$, maximum moment order $N_{m}$, and
spatial dimension $D$
Range of density, velocity, and temperature to formulate $\hat{\bm{\alpha}}$
Standard deviation $\sigma$ to sample $\bm{\gamma}$
Condition number threshold $c$

Result: Dataset of reference distribution functions $f$

for  _$i=1,\dots,2N_{i}$_ do

Sample $\rho\in[\rho_{min},\rho_{max}]$

Sample $\mathbf{V}\in[\mathbf{V}_{min},\mathbf{V}_{max}]$

Sample $T\in[T_{min},T_{max}]$

Compute mean for
$(\gamma_{u,1},\bm{\gamma}_{u,2},\gamma_{u,3})^{T}=\hat{\bm{\alpha}}$

do

Sample $\gamma_{u,n}$, $n=4,\dots,N_{m}$

Compute $\mathbf{H}_{u}(\bm{\gamma}_{u})$

while  _$\sigma_{\mathbf{H}_{u}(\bm{\gamma}_{u})} <c$_

Compute $f_{\text{ent},i}=\exp(\bm{\gamma}_{u}\cdot\bm{\phi})$

end for

Divide $f_{\mathrm{ent}}$ into $f_{L}$ and $f_{R}$ equally

for  _$i=1,\dots,N_{i}$_ do

Sample unit vector $\mathbf{n}$

Reconstruct
$f_{\mathrm{ref},i}=f_{L,i}\mathscr{H}(\mathbf{n}\cdot\mathbf{v})+f_{R,i}(1-\mathscr{H}(\mathbf{n}\cdot\mathbf{v}))$

end for

Algorithm 1 Sampling of reference particle distribution functions

###  4.2 Computation of collision term

After the particle distribution functions are collected, the collision terms
need to be computed and used in the cost function Eq.(20). In this paper, we
consider three different collision terms for the Boltzmann and its model
equations, which are briefly introduced below.

  1. (1)

Full Boltzmann collision integral with fast spectral method: The full
Boltzmann term can be transformed into Carleman-type [39],

$\displaystyle\mathcal{Q}(f,f)$
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\Theta
g\left[f\left(\mathbf{v}+\frac{g\Omega-\mathbf{g}}{2}\right)f\left(\mathbf{v}_{*}-\frac{g\Omega-\mathbf{g}}{2}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]d\Omega
d\mathbf{v}_{*}$ (64)
$\displaystyle=2\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\Theta\delta\left(2\mathbf{y}\cdot\mathbf{g}+\mathbf{y}^{2}\right)\left[f\left(\mathbf{v}+\frac{\mathbf{y}}{2}\right)f\left(\mathbf{v}_{1}-\frac{\mathbf{y}}{2}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]d\mathbf{y}d\mathbf{v}_{*}$
$\displaystyle=4\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\Theta\delta(\mathbf{y}\cdot\mathbf{z})[f(\mathbf{v}+\mathbf{y})f(\mathbf{v}+\mathbf{z})-f(\mathbf{v})f(\mathbf{v}+\mathbf{y}+\mathbf{z})]d\mathbf{y}d\mathbf{z},$

which can be solved with a discrete Fourier transform-based spectral method
with a computational cost of $O(N^{3}\log N)$ [40]. Specifically, the particle
distribution function and collision term discretized with $N_{v}$ quadrature
points are expanded into the Fourier series,

$\displaystyle
f(t,\mathbf{x},\mathbf{v})=\sum_{k=-N_{v}/2}^{N_{v}/2-1}{f}_{k}(t,\mathbf{x})\exp\left(i\xi_{k}\cdot\mathbf{v}\right),$
(65) $\displaystyle{f}_{k}=\frac{1}{L^{D}}\langle
f(t,\mathbf{x},\mathbf{v})\exp\left(-i\xi_{k}\cdot\mathbf{v}\right)\rangle,$
$\displaystyle\mathcal{Q}(t,\mathbf{x},\mathbf{v})=\sum_{k=-N_{v}/2}^{N_{v}/2-1}\mathcal{Q}_{k}(t,\mathbf{x})\exp\left(i\xi_{k}\cdot\mathbf{v}\right),$
$\displaystyle\mathcal{Q}_{k}=\sum_{l,m=-N_{v}/2,(l+m=k)}^{N_{v}/2-1}{f}_{l}{f}_{m}[\beta(l,m)-\beta(m,m)],$

where $L$ is the span of velocity space in each dimension, and $\beta$ is the
kernel mode, and the convolutions defined in Eq.(64) can be solved in the
frequency domain. The details of this method can be found in [40], and its
numerical implementation is available in [1].

  2. (2)

Shakhov relaxation model: The Shakhov relaxation model builds a heat flux-
based correction term on top of the BGK model to ensure that the Chapman-
Enskog expansion of the model has the correct Prandtl number [18]. The
equilibrium state in the Shakhov model writes

$\mathcal{E}=\mathcal{M}\left[1+(1-\mathrm{Pr})\mathbf{c}\cdot\mathbf{q}\left(\frac{\mathbf{c}^{2}}{RT}-5\right)/(5pRT)\right],$
(66)

where Pr is the Prandtl number, $\mathbf{c}$ is the peculiar velocity,
$\mathbf{q}$ is the heat flux, and $R=k/m$ is the gas constant. It has been
shown that the Shakhov model is more accurate in predicting highly non-
equilibrium flows in normal shock wave compared to the BGK model [41]. With
the inclusion of heat flux, the Shakhov model loses some of the structural
properties shown in Section 2. For example, non-negative solution of particle
distribution function can appear, and the H-theorem is no longer strictly
valid in the Shakhov model. Using the Shakhov model to construct the dataset
is to validate the ability of RelaxNet to correctly approximate $\mathcal{E}$.

  3. (3)

$\nu$-BGK model: Another strategy to improve the BGK model is to introduce
velocity-dependent relaxation frequency. Here we choose the $\nu$-BGK model
[19], where the relaxation term writes

$\mathcal{Q}=\nu(\mathbf{v})(\mathcal{M}-f).$ (67)

Different curves of $\nu$ as a function of $\mathbf{v}$ can be fitted to
provide the correct Prandtl number of monatomic gas. Here we adopt the
relaxation frequency defined in [42], which is an improvement over the
frequencies provided in [19] and can provide more accurate solution in non-
equilibrium flows. In this model, the relaxation frequency is defined as

$\displaystyle\nu(|\mathbf{v}|)=A\frac{p}{\mu}\left[\nu_{eq}^{0}(\xi)+2\nu_{eq}^{0}(\xi)\right],$
(68)
$\displaystyle\nu_{eq}^{0}(\xi)=\frac{3}{2}\left[\exp(-\xi^{2})+\frac{\sqrt{\pi}}{2}\left(\frac{1}{\xi}+2\xi\right)\mathrm{erf}(\xi)\right],$

where $\xi=|\mathbf{c}|/\sqrt{2kT/m}$, and $A$ is a parameter depending on the
viscosity index. Using the $\nu$-BGK model to construct the dataset is to
validate the ability of RelaxNet to correctly approximate velocity-dependent
$\nu$.

##  5 Solution Algorithm

###  5.1 Update algorithm

The finite volume method is employed to build the solution algorithm of the
RelaxNet-based universal Boltzmann equation. We denote the mean value of
particle distribution function in a control volume as

$f(t^{n},\mathbf{x}_{i},\mathbf{v}_{j})=f_{i,j}^{n}=\frac{1}{\Omega_{i}(\mathbf{x})\Omega_{j}(\mathbf{v})}\int_{\Omega_{i}}\int_{\Omega_{j}}f(t^{n},\mathbf{x},\mathbf{v})d\mathbf{x}d\mathbf{v},$
(69)

where $\Omega_{i}$ and $\Omega_{j}$ are cell areas in the discrete physical
and velocity space. The update algorithm of distribution function in each
control volume can be written as

$f_{i,j}^{n+1}=f_{i,j}^{n}+\frac{1}{\Omega_{i}}\int_{t^{n}}^{t^{n+1}}\sum_{r=1}^{n_{f}}F_{r}\Delta
S_{r}dt+\int_{t^{n}}^{t^{n+1}}\mathrm{NN}_{\theta}(f_{i},\bar{\tau}_{i})_{j}dt,$
(70)

where $F_{r}$ is the time-dependent flux function of distribution function at
$r$-th cell interface, $\Delta S_{r}$ is the interface area, and $n_{f}$ is
the number of interfaces. Details of the update algorithm can be found in
[43].

###  5.2 Numerical flux

The numerical flux in Eq.(70) is evaluated from the reconstructed particle
distribution function at the cell interface. At the cell interface
$\mathbf{x}_{i+1/2}$, the distribution function is constructed in an upwind
fashion,

$f_{i+1/2,j}=f_{i+1/2,j}^{L}\mathscr{H}\left(\mathbf{n}_{i+1/2}\cdot\mathbf{v}_{j}\right)+f_{i+1/2,j}^{R}(1-\mathscr{H}\left(\mathbf{n}_{i+1/2}\cdot\mathbf{v}_{j}\right)),$
(71)

where $\mathbf{n}_{i+1/2}$ is the unit normal vector of the interface, and
$\mathscr{H}(x)$ is the heaviside step function. The left and right value of
distribution function $f_{i+1/2,j}^{L,R}$ can be obtained through
reconstruction and interpolation. For example, a second-order interpolation
results in

$\displaystyle
f_{i+1/2,j}^{L}=f_{i,j}+\nabla_{\mathbf{x}}f_{i,j}\cdot(\mathbf{x}_{i+1/2}-\mathbf{x}_{i}),$
(72) $\displaystyle
f_{i+1/2,j}^{R}=f_{i+1,j}+\nabla_{\mathbf{x}}f_{i+1,j}\cdot(\mathbf{x}_{i+1/2}-\mathbf{x}_{i+1}),$

where $\nabla_{\mathbf{x}}f$ is the reconstructed gradient with slope
limiters. The principle of higher-order interpolation can be found in [44].
The numerical flux of particle distribution function can then be evaluated as

$F_{i+1/2,j}=f_{i+1/2,j}\mathbf{n}_{i+1/2}\cdot\mathbf{v}_{j}.$ (73)

###  5.3 Collision term

The RelaxNet-based collision operator in each cell reads

$\mathcal{Q}_{i,j}=\mathrm{NN}_{\theta}(f_{i},\bar{\tau}_{i})_{j}.$ (74)

According to the kinetic theory of gases, the mean relaxation time is defined
as $\bar{\tau}=\mu/{p}$ where $p$ is pressure. In this paper, we adopt the
variational hard-sphere (VHS) molecule model to compute the viscosity
coefficient, i.e.,

$\mu=\mu_{\mathrm{ref}}\left(\frac{T}{T_{\mathrm{ref}}}\right)^{\omega},$ (75)

where $\mu_{\mathrm{ref}}$ and $T_{\mathrm{ref}}$ are the viscosity and
temperature in the reference state, and $\omega$ is the viscosity index.

Once the parameters in RelaxNet get optimized, the collision term can be
integrated with the help of time-integral algorithms. While the forward Euler
method is the most straightforward integrator for solving Eq.(70), it is
feasible to adopt other higher-order methods, e.g., Tsitouras’ 5/4 Runge-Kutta
method [45]. When the time step is much larger than mean relaxation time, it
is preferable to employ implicit [46] and implicit-explicit (IMEX) methods
[47] to improve numerical stability of the scheme.

Note that the update procedure in Eq.(70) produces case-specific distribution
functions that may be beyond the dataset used for training and testing. The
newly generated distribution functions can be supplemented to the dataset on-
the-fly by computing the collision term
$\mathcal{Q}_{i,j}=\mathcal{Q}(f_{i,j})$ with the methods shown in Section
4.2. We use the rise of residuals as a criterion to generate supplementary
data. A schematic diagram of the solution algorithm of the RelaxNet-based UBE
is presented in Figure 4.

Figure 4: Schematic diagram of solution algorithm of the RelaxNet-based
universal Boltzmann equation.

##  6 Numerical Experiments

In this section, we present numerical experiments to validate the RelaxNet-
based universal Boltzmann equation and its solution algorithm. The
dimensionless variables are introduced in the simulation, i.e.,

$\displaystyle\tilde{t}=\frac{t}{L_{0}/(2RT_{0})^{1/2}},\
\tilde{\mathbf{x}}=\frac{\mathbf{x}}{L_{0}},\
\tilde{\rho}=\frac{\rho}{\rho_{0}},\
\tilde{\mathbf{v}}=\frac{\mathbf{v}}{(2RT_{0})^{1/2}},\
\tilde{\mathbf{V}}=\frac{\mathbf{V}}{(2RT_{0})^{1/2}},\
\tilde{T}=\frac{T}{T_{0}},$ $\displaystyle\tilde{E}=\frac{E}{2RT_{0}},\
\tilde{\mathbf{P}}=\frac{\mathbf{P}}{2\rho_{0}RT_{0}},\
\tilde{\mathbf{q}}=\frac{\mathbf{q}}{\rho_{0}(2RT_{0})^{3/2}},\
\tilde{\mu}=\frac{\mu}{\rho_{0}L_{0}(2RT_{0})^{1/2}},\
\tilde{f}=\frac{f}{\rho_{0}/(2RT_{0})^{D/2}}.$

The subscript zero represents the reference state. For brevity, the tilde
notation for dimensionless variables will be removed henceforth. In all
simulations we employ monatomic hard-sphere gas, where the viscosity
coefficient can be determined by Eq.(75).

###  6.1 Data generation

To further explain the principle of the data sampling strategy described in
Section 4 and illustrate its superiority, we first compare the distribution of
the data generated by Algorithm 1 with the data collected directly from
numerical simulation results of the Boltzmann equation. The computational
setup for Algorithm 1 is presented in Table 2, where
$\\{\rho_{s},U_{s},T_{s}\\}$ is the macroscopic variables used for sampling
distribution functions, $\bm{\gamma}$ is the Lagrange multipliers in Eq.(58),
and $\mathcal{U}(a,b)$, $\mathcal{N}(\mu,\sigma)$, and
$\mathcal{T}(\mu,\sigma,a,b)$ denote uniform, normal, and truncated normal
distributions, respectively.

Table 2: Computational setup for sampling particle distribution functions in the data generation problem. $N_{i}$ | $N_{m}$ | $D$ | $\rho_{s}$ | $U_{s}$ | $T_{s}$ | $\gamma_{4}$ | $\gamma_{5}$
---|---|---|---|---|---|---|---
$3000$ | 4 | 1 | $\mathcal{U}(0,2)$ | $\mathcal{U}(-1.5,1.5)$ | $\mathcal{U}(0,4)$ | $\mathcal{N}(0,0.005)$ | $\mathcal{T}(0,0.005,-\infty,0)$

Two problems are considered below as illustrative examples to produce
simulation-based datasets.

  1. (1)

Wave propagation problem. The initial particle distribution function is set to
Maxwellian everywhere in correspondence with the following flow field,

$\left(\begin{array}[]{c}\rho\\\ U\\\ p\\\
\end{array}\right)_{{t=0}}=\left(\begin{array}[]{c}A(1+a\sin(2\pi x))\\\ b\\\
c\\\ \end{array}\right).$

The detailed computational setup is presented in Table 3.

  2. (2)

Sod shock tube problem. The initial particle distribution function is set to
Maxwellian everywhere in correspondence with the one-dimensional Riemann
problem,

$\left(\begin{array}[]{c}\rho\\\ U\\\ p\\\
\end{array}\right)_{{t=0,x_{L}}}=\left(\begin{array}[]{c}\rho_{L}\\\ 0\\\
p_{L}\\\ \end{array}\right),\ \left(\begin{array}[]{c}\rho\\\ U\\\ p\\\
\end{array}\right)_{{t=0,x_{R}}}=\left(\begin{array}[]{c}\rho_{R}\\\ 0\\\
p_{R}\\\ \end{array}\right).$

The detailed computational setup is presented in Table 4.

Table 3: Computational setup for wave propagation problem. $t$ | $x$ | $N_{x}$ | $u$ | $N_{u}$ | Quadrature | Kn | $A_{s}$
---|---|---|---|---|---|---|---
$(0,5]$ | $[0,1]$ | $200$ | $[-8,8]$ | 80 | Rectangular | $10^{-3}$ | $\mathcal{U}(0.5,1.5)$
|  |  |  |  |  |  |
$a_{s}$ | $b_{s}$ | $c_{s}$ | Collision | Integrator | Boundary | CFL |
$\mathcal{U}(0,0.2)$ | $\mathcal{U}(-1,1)$ | $\mathcal{U}(0.5,1.5)$ | Shakhov | RK4 | Periodic | 0.5 |
Table 4: Computational setup for Sod shock tube problem. $t$ | $x$ | $N_{x}$ | $u$ | $N_{u}$ | Quadrature | Kn | $\rho_{Ls}$
---|---|---|---|---|---|---|---
$(0,0.2]$ | $[0,1]$ | $200$ | $[-8,8]$ | 80 | Rectangular | $10^{-4}$ | $\mathcal{U}(0,2)$
|  |  |  |  |  |  |
$\rho_{Rs}$ | $T_{Ls}$ | $T_{Rs}$ | Collision | Integrator | Boundary | CFL |
$\mathcal{U}(0,0.25)$ | $\mathcal{U}(1,3)$ | $\mathcal{U}(0,1.25)$ | Shakhov | IMEX | Dirichlet | 0.5 |

We compare the datasets generated from the algorithmic sampler described in
Section 4 and the above two numerical examples with the same sample size
$N_{i}=3000$. Two cases are considered. In the first case, the input
parameters of two examples are fixed as $\\{A,a,b,c\\}=\\{1,0.1,1,1\\}$ and
$\\{\rho_{L},\rho_{R},T_{L},T_{R}\\}=\\{1,0.125,2,0.625\\}$ to analyze the
dataset produced by a single numerical simulation. In the second case, we
sample the parameters ten times according to the uniform distribution of the
intervals given in Table 3 and 4 to characterize the datasets collected from
multiple numerical simulations with various inputs.

Figure 5 and 6 present the phase diagrams between macroscopic variables
generated by the algorithmic sampler using Algorithm 1 and the simulation
results. As shown in Figure 5, the samples from a single numerical simulation
can suffer from a strong bias. Such a phenomenon is caused by the physical
constraints in a particular flow problem. In the wave propagation problem, the
density and temperature is linearly coupled due to the isobaric pressure. In
the Sod shock tube problem, the density, velocity and temperature are
correlated around wave structures, i.e., the rarefaction wave, contact
discontinuity, and shock wave. The discontinuous flow field in the Sod problem
makes the samples show a even sparser distribution in the phase diagram. As
particle distribution functions generated from multiple sampling of parameters
and numerical simulations are added to the dataset, as presented in Figure 6,
a wider range of flow variables can be obtained. From this trend, it may be
feasible to generate reliable datasets through numerical simulations if the
sampling ranges of parameters are further extended, but this obviously
requires the use of appropriate physical problems and adequate sample size. To
accomplish the goal, the use of trial-and-error processes is inevitable.

In contrast, the algorithmic sampler generates a dataset of particle
distribution function with a much wider range of densities, velocities and
temperatures. The obtained samples are homogeneously distributed without a
bias towards a preferred direction. Therefore, the predefined domain of
$\\{\rho,U,T\\}$ is fully covered with a relatively small sample size. Note
that this also implies a broader and unbiased distribution of particle
distribution functions, since there is a one-to-one correspondence between
distribution functions and macroscopic quantities. Compared with the
simulation-based sampling, the broader dataset provided by the algorithmic
sampler reduces the risk of overfitting and helps improve the generalization
of the neural network model. Table 5 provides the computation time and
allocations for the three methods (including the garbage collection with Julia
programming language [48]). Since the presented sampling strategy does not
require a full simulation of the Boltzmann equation with multiple cases or
initial conditions, the computational cost is significantly reduced by three
to five orders of magnitude.

Table 5: Computational cost of the data generation problem. | Time | Allocations | Memory |
---|---|---|---|---
Wave |  $197.12$ s | $1.49\times 10^{9}$ |  $485.94$ GB |
Sod |  $41.46$ s | $4.18\times 10^{6}$ |  $13.15$ GB |
Sampler |  $0.03$ s | $7.27\times 10^{3}$ |  $13.33$ MB |

###  6.2 Homogeneous relaxation

In this numerical experiment, we study the homogeneous relaxation problem of
non-equilibrium particle distribution functions. To fully validate the ability
of RelaxNet to approximate the solution of the Boltzmann and its model
equations, three collision models, i.e., the Shakhov model, the velocity-
dependent $\nu$-BGK model, and the full Boltzmann integral model described in
Section 4.2, are employed to provide ground-truth datasets for training the
neural network.

####  6.2.1 Shakhov model

First we consider the non-equilibrium distribution function in the initial
state,

$f(t=0,u)=\frac{1}{2\sqrt{\pi}}\left(\exp\left(-\left(u-1.5\right)^{2}\right)+\frac{7}{10}\exp\left(-(u+1.5)^{2}\right)\right).$
(76)

We construct RelaxNet in such a way that there are two hidden layers of the
same dimension as the input layer in both $\mathcal{E}$-net and $\tau$-net. A
dataset is generated by Algorithm 1 to train RelaxNet, which is then used to
solve the universal Boltzmann equation. The setups of the neural network,
training and computation are provided in Table 6, where $N_{h}=2$ is the
number of hidden layers and Tsit5 denotes Tsitouras’ 5/4 Runge-Kutta method
[45].

Table 6: Computational setup for the homogeneous relaxation problem based on the Shakhov model. $D$ | $u$ | $N_{u}$ | $N_{i}$ | $\rho_{s}$ | $U_{s}$
---|---|---|---|---|---
1 | $[-8,8]$ | 80 | 10000 | $\mathcal{U}(0.1,10)$ | $\mathcal{U}(-1,1)$
|  |  |  |  |
$V_{s}$ | $W_{s}$ | $T_{s}$ | $\gamma_{4}$ | $\gamma_{5}$ | $\mathrm{Kn}_{s}$
$0$ | $0$ | $\mathcal{U}(0.1,5)$ | $\mathcal{N}(0,0.005)$ | $\mathcal{T}(0,0.005,-\infty,0)$ | $\mathcal{U}(0.001,1)$
|  |  |  |  |
$N_{h}$ | $N_{m}$ | Reference | Optimizer | $t$ | Kn
2 | 4 | Boltzmann | Adam | $(0,8]$ | 1
|  |  |  |  |
Integrator | Quadrature | CFL |  |  |
Tsit5 | Rectangular | 0.3 |  |  |

The values of cost function in Eq.(20) for both training and test datasets
during the training process are plotted in Figure 7. After the parameters of
RelaxNet are optimized, the initial value problem of the UBE model is solved.
Note that the bimodal distribution in Eq.(76) can fall outside the training
set, so this example actually acts as a validation set. Figure 8 presents the
solutions of particle distribution functions at different time instants
computed by the BGK, Shakhov, and UBE model. It can be seen that the numerical
solution provided by UBE is closer to the Shakhov reference solution than the
BGK model. This is attributed to a more accurate approximation of the
collision term of the Boltzmann equation, which is shown in Figure 9. Figure
10 presents the evolution of macroscopic density and energy with time. It is
clear that the conservation of mass is exactly satisfied, which provides a
corroboration of Theorem 2. The error in energy slightly exceeds the range of
the Shakhov reference solution. Due to the interpretability of RelaxNet, we
can open the black box of forward pass in the UBE model and explain the issue.
Figure 11 presents the equilibrium state $\mathcal{E}$ approximated by
RelaxNet at different time instants, which intuitively explains the mechanism
of UBE. Note that the Shakhov model does not satisfy $H$-theorem and the
logarithm of equilibrium distribution corrected by heat flux cannot be written
as a sum of collision invariants. As a result, the approximation of RelaxNet
deviates a bit from the reference Shakhov solution.

####  6.2.2 $\nu$-BGK model

Next we consider the same setup as in Section 6.2.1, except that the velocity-
dependent $\nu$-BGK model shown in Section 4.2 is employed to produce
reference solution. Figure 12 presents the solutions of particle distribution
functions at different time instants computed by the BGK, the $\nu$-BGK, and
the UBE model. It is clear that there is a more significant difference between
the BGK and $\nu$-BGK solutions, while the solution provided by the UBE model
is in excellent agreement with the reference $\nu$-BGK solution. As we
explained $\mathcal{E}$ in Section 6.2.1, RelaxNet allows us to quantitatively
interpret the relaxation frequency $\nu$ in the UBE model. Figure 13 and 14
show the collision terms and relaxation frequencies used in RelaxNet, and
compare them with the BGK and $\nu$-BGK models. It is clear that the
dependence between relaxation frequencies and particle velocities are
precisely recovered, bringing a more accurate approximation of the collision
term. Figure 15 presents the evolution of macroscopic density and energy with
time. It can be seen that although the parameters in RelaxNet are optimized
based on the solution of $\nu$-BGK model in a supervised learning manner, the
physics-informed regularization in the cost function in Eq.(20) allows UBE to
better preserve the conservation property, which provides a corroboration of
Theorem 2.

####  6.2.3 Full Boltzmann model

We then turn to the full Boltzmann equation. In this case, the initial
particle distribution function is set as

$f(t=0,u,v,w)=\frac{1}{\sqrt{\pi}}\left(\exp\left(-\left(u-1\right)^{2}\right)+\frac{7}{10}\exp\left(-(u+1)^{2}\right)\right)\exp\left(-v^{2}\right)\exp\left(-w^{2}\right).$

RelaxNet is constructed in the same way as in Section 6.2.1 and 6.2.2.
Following Algorithm 1, the fast spectral method introduced in Section 4.2 is
employed to generate the dataset used for training RelaxNet. The computational
setup is presented in Table 7.

Table 7: Computational setup for the homogeneous relaxation problem based on the full Boltzmann model. $D$ | $\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | $N_{i}$ |
---|---|---|---|---|---|---
3 | $[-8,8]^{3}$ | 80 | 28 | 28 | 10000 |
|  |  |  |  |  |
$\rho_{s}$ | $U_{s}$ | $V_{s}$ | $W_{s}$ | $T_{s}$ | $\gamma_{4}$ |
$\mathcal{U}(0.1,10)$ | $\mathcal{U}(-1,1)$ | $0$ | $0$ | $\mathcal{U}(0.1,5)$ | $\mathcal{N}(0,0.005)$ |
|  |  |  |  |  |
$\gamma_{5}$ | $\mathrm{Kn}_{s}$ | $N_{h}$ | $N_{m}$ | Reference | Optimizer |
$\mathcal{T}(0,0.005,-\infty,0)$ | $\mathcal{U}(0.001,1)$ | 2 | 4 | Boltzmann | Adam |
|  |  |  |  |  |
$t$ | Kn | Integrator | Quadrature | CFL |  |
$(0,8]$ | 1 | Tsit5 | Rectangular | 0.3 |  |

Figure 16 provides the evolution of particle distribution functions with time.
The variable used for display is the reduced distribution function on $u$
axis, i.e.,

$h(t,u)=\int_{\mathbb{R}^{2}}f(t,u,v,w)dvdw.$

As shown, UBE based on RelaxNet provides a solution consistent with the full
Boltzmann equation, while the solution of the BGK model shows significant
deviations. Figure 17 presents the collision terms computed by the three
approaches at different time instants. Figure 18 shows the equilibrium
distributions and relaxation frequencies approximated by RelaxNet. It can be
seen that the difference between the equilibrium distribution in RelaxNet and
the Maxwellian distribution is slight, and thus the correction of the BGK
solution is mainly achieved by the velocity-dependent relaxation frequency.
Different from predefined frequencies in the $\nu$-BGK model, e.g., in
Eq.(68), RelaxNet provides a variable relaxation frequency that depends on the
solution of local particle distribution function, and thus provides an
approximation that better fits the real collision frequency in Eq.(2). Figure
19 shows the evolution of entropy density defined as $\eta=\langle
f\log(f)\rangle$, from which the principle of increase of entropy and the
H-theorem demonstrated in Section 2 adn Theorem 4 are validated.

Table 8 presents the computational costs of the three methods to simulate the
time-series solution of the homogeneous relaxation problem. Since a large
number of convolution operations for solving the Boltzmann collision integral
are replaced by tensor multiplication, the computational cost of RelaxNet is
significantly reduced compared. Under the current numerical setup, the
simulation time is reduced to $1/50$ of the fast spectral method for the full
Boltzmann equation, and the memory usage is reduced to $1/17$. This numerical
experiment verifies the ability of the RelaxNet-based UBE to solve the
Boltzmann equation accurately and efficiently.

Table 8: Computational cost of computing time-series solutions in the homogeneous relaxation problem. | Time | Allocations | Memory |
---|---|---|---|---
Boltzmann |  $2.167$ s | $90322$ |  $6.43$ GB |
BGK |  $1.076$ ms | $2939$ |  $3.38$ MB |
UBE |  $43.472$ ms | $14344$ |  $375.66$ MB |

###  6.3 Normal shock structure

We then study flow transport in inhomogeneous gases. Due to the existence of
non-Maxwellian particle distribution functions [49], the normal shock wave
problem is an ideal case to validate the RelaxNet-based UBE in solving highly
non-equilibrium flows. The initial flow field is initialized according to the
Rankine-Hugoniot jump conditions, i.e.,

$\left(\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ T\\\
\end{array}\right)_{{t=0,x_{L}}}=\left(\begin{array}[]{c}\rho_{-}\\\ U_{-}\\\
0\\\ 0\\\ T_{-}\\\ \end{array}\right),\ \left(\begin{array}[]{c}\rho\\\ U\\\
V\\\ W\\\ T\\\
\end{array}\right)_{{t=0,x_{R}}}=\left(\begin{array}[]{c}\rho_{+}\\\ U_{+}\\\
0\\\ 0\\\ T_{+}\\\ \end{array}\right),$

where the upstream and downstream density, velocity, and temperature are
denoted by $\\{\rho_{-},U_{-},T_{-}\\}$ and $\\{\rho_{+},U_{+},T_{+}\\}$,
respectively. The initial particle distribution functions are set to
Maxwellian according to the local flow variables,

$\displaystyle\frac{\rho_{+}}{\rho_{-}}=\frac{(\varpi+1)\rm{Ma}^{2}}{(\varpi-1)\rm{Ma}^{2}+2},$
$\displaystyle\frac{U_{+}}{U_{-}}=\frac{(\varpi-1)\rm{Ma}^{2}+2}{(\varpi+1)\rm{Ma}^{2}},$
$\displaystyle\frac{T_{+}}{T_{-}}=\frac{((\varpi-1)\rm{Ma}^{2}+2)(2\varpi\rm{Ma}^{2}-\varpi+1)}{(\varpi+1)^{2}\rm{Ma}^{2}},$

where $\rm Ma$ is the Mach number, and $\varpi$ is the ratio of specific heat.
The upstream flow variables are used as references to dimensionaless the
system. The detailed computational setup is provided in Table 9. Note that the
dataset used to optimize the parameters of RelaxNet consists of two parts. In
addition to the data generated by Algorithm 5 with predefined parameters, as
the numerical simulation progresses, the case-specific data can be
supplemented on the fly, as shown in Figure 4.

Table 9: Computational setup for the normal shock structure problem. $D$ | $\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | $N_{i}$ |
---|---|---|---|---|---|---
3 | $[-10,10]^{3}$ | 80 | 28 | 28 | 10000 |
|  |  |  |  |  |
$\rho_{s}$ | $U_{s}$ | $V_{s}$ | $W_{s}$ | $T_{s}$ | $\gamma_{4}$ |
$\mathcal{U}(0.5,5)$ | $\mathcal{U}(-3,3)$ | $0$ | $0$ | $\mathcal{U}(0.5,5)$ | $\mathcal{N}(0,0.005)$ |
|  |  |  |  |  |
$\gamma_{5}$ | $\mathrm{Kn}_{s}$ | $N_{h}$ | $N_{m}$ | Reference | Optimizer |
$\mathcal{T}(0,0.005,-\infty,0)$ | $\mathcal{U}(0.01,5)$ | 2 | 4 | Boltzmann | Adam |
|  |  |  |  |  |
$t$ | $x$ | $N_{x}$ | Kn | Ma | $\varpi$ |
$(0,250]$ | $[-35,35]$ | 100 | 1 | $[2,3]$ | $5/3$ |
|  |  |  |  |  |
Integrator | Boundary | Quadrature | CFL |  |  |
Tsit5 | Dirichlet | Rectangular | 0.5 |  |  |

Figure 20 shows the profiles of density, velocity, temperature, viscous
stress, heat flux, and entropy density in the shock structure at $\rm Ma=2$.
The definition of viscous stress is the difference between the first element
in the pressure tensor and the isotropic pressure, i.e.,

$\Delta=\mathbf{P}_{xx}-p.$

When the Mach number is small, the particle distribution functions inside the
shock wave deviate moderately away from Maxwellian, and thus the Boltzmann,
BGK, and UBE solutions are roughly equivalent. The discrepancy of the viscous
stress and heat flux provided by the Boltzmann and BGK models can be observed
due to the fact that high-order moments are more sensitive to the nuances of
distribution function than density and bulk velocity. Figure 21 presents the
macroscopic flow field at $\rm Ma=3$. It is obvious that the difference
between the Boltzmann and BGK solutions become more significant with the
increasing Mach number. The unit Prandtl number in the BGK model leads to an
incorrect heat transfer rate and an early rise in the temperature profile.
Thanks to the solution-dependent equilibrium and relaxation frequency, the
current UBE is able to provide solutions equivalent to the reference Boltzmann
results in all cases.

Figure 22 and 23 provide the reduced distribution functions and collision
terms in the $x-u$ phase space at $\rm Ma=2$ and $\rm Ma=3$. An one-to-one
correspondence can be clearly observed between the evolution of macroscopic
quantities and particle distribution functions. When $\rm Ma=3$, significant
difference can be obtained between the BGK and UBE collision terms, resulting
in different distributions of particle distribution functions and macroscopic
flow variables. Table 10 shows the computational cost of computing the right-
hand side operator once in the Boltzmann, BGK, and UBE models. The benchmark
verifies that the RelaxNet-based UBE increases the computational efficiency by
a factor of 108 and reduces the memory load to $1/43$ compared to the fast
spectral method of the Boltzmann equation.

Table 10: Computational cost of computing the right-hand side of the kinetic model in the normal shock structure problem. | Time | Allocations | Memory |
---|---|---|---|---
Boltzmann |  $32.114$ ms | $1446$ |  $102.98$ MB |
BGK |  $5.893$ $\rm\mu s$ | $36$ |  $12.92$ KB |
UBE |  $295.947$ $\rm\mu s$ | $93$ |  $2.38$ MB |

###  6.4 Lid-driven cavity

In the last numerical experiment, we study the lid-driven cavity to validate
the ability of current approach to solve non-equilibrium flows under multi-
dimensional geometry. The initial flow field inside the square cavity is set
as

$\left(\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ T\\\
\end{array}\right)_{{t=0}}=\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\ 1\\\
\end{array}\right),$

and the particle distribution function is set to Maxwellian everywhere. The
flow domain is enclosed by four isothermal solid walls. The upper wall moves
in the tangent direction with $\mathbf{V}_{w}^{u}=(0.15,0,0)^{T}$, with the
rest three walls kept still with $\mathbf{V}_{w}^{d,l,r}=(0,0,0)^{T}$.
Maxwell’s diffusive boundary is considered at all wall boundaries. The
detailed computational setup is provided in Table 11.

Table 11: Computational setup for the lid-driven cavity problem. $D$ | $\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | $N_{i}$ |
---|---|---|---|---|---|---
3 | $[-5,5]^{3}$ | 28 | 28 | 28 | 10000 |
|  |  |  |  |  |
$\rho_{s}$ | $U_{s}$ | $V_{s}$ | $W_{s}$ | $T_{s}$ | $\gamma_{4}$ |
$\mathcal{U}(0.5,5)$ | $\mathcal{U}(-1,1)$ | $\mathcal{U}(-1,1)$ | $0$ | $\mathcal{U}(0.5,5)$ | $\mathcal{N}(0,0.005)$ |
|  |  |  |  |  |
$\gamma_{5}$ | $\mathrm{Kn}_{s}$ | $N_{h}$ | $N_{m}$ | Reference | Optimizer |
$\mathcal{T}(0,0.005,-\infty,0)$ | $\mathcal{U}(0.01,1)$ | 2 | 4 | Boltzmann | Adam |
|  |  |  |  |  |
$\mathbf{x}$ | $N_{x}$ | $N_{y}$ | $\mathbf{V}_{w}^{u}$ | $\mathbf{V}_{w}^{d,l,r}$ | $T_{w}$ |
$[0,1]^{2}$ | 45 | 45 | $(0.15,0,0)^{T}$ | $(0,0,0)^{T}$ | 1 |
|  |  |  |  |  |
$\varpi$ | Kn | Integrator | Boundary | Quadrature | CFL |
$5/3$ | 0.075 | Tsit5 | Maxwell | Rectangular | 0.5 |

Figure 24 shows the contours of density with streamlines and temperature with
heat flux vectors inside the cavity. Consistent with the report in [50], the
inverse-Fourier heat flux driven by the viscous stress in the non-equilibrium
regime is accurately identified. Figure 25 presents the velocity profiles
along the horizontal and vertical central lines of the cavity. The validity of
UBE and the corresponding numerical method is fully validated by a
quantitative comparison with the reference DSMC solution. Figure 26 and 27
further provide the particle distribution functions and collision terms along
the two central lines. For low-speed microflows, the non-equilibrium effect is
weak and the difference between the Boltzmann and BGK solutions is mild. The
correction provided by the RelaxNet-based UBE for the BGK model can be
observed, which ensures the correct heat transfer rate and Prandtl number.

## Conclusion

Scientific machine learning is increasingly showing its power and versatility
in modeling and simulation in computational physics. In this paper, a novel
relaxation neural network (RelaxNet) is proposed as a surrogate model of the
collision operator in the Boltzmann equation. Based on the architecture of
residual neural networks, RelaxNet builds a parameterized relaxation model
with solution-dependent equilibrium state and frequency. The quantitative
interpretability of RelaxNet makes it the first machine learning model that
strictly preserves all the key structural properties of the Boltzmann
equation, including invariance, conservation, H-theorem, and the correct
continuum limit. Since the convolution in the original Boltzmann equation is
simplified by the tensor multiplication in RelaxNet, the dimensionality and
computational complexity of the collision term are greatly reduced. Based on
RelaxNet, a universal Boltzmann equation (UBE) model is developed, which fuses
mechanical and neural models into a single differentiable framework that is
compatible with source-to-source automatic differentiation. Based on the
entropy closure of the Boltzmann moment system, an algorithmic sampler is
designed to generate datasets for model training and testing, and the
superiority of this strategy over simulation-based samplers is demonstrated
through numerical experiments. The solution algorithm for solving the
RelaxNet-based UBE model is described in detail. Both spatially homogeneous
and inhomogeneous test cases are presented to validate UBE and its solution
algorithm. The current approach provides an accurate and efficient tool for
the study of non-equilibrium flow physics. Like the network chain that
contains multiple ResNet blocks, it is promising to extend RelaxNet to build
multi-level relaxation models [51], which remains to be investigated in future
work.

## CRediT authorship contribution statement

Tianbai Xiao: : Conceptualization, Formal analysis, Investigation,
Methodology, Project administration, Software, Visualization, Writing –
original draft, Writing – review & editing. Martin Frank: Formal analysis,
Project administration, Resources, Supervision, Writing – review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgement

The current research is funded by the German Research Foundation in the frame
of the priority programm SPP 2298 ”Theoretical Foundations of Deep Learning”
with the project number 441826958.

## References

  * [1] Tianbai Xiao.  Kinetic.jl: A portable finite volume toolbox for scientific and neural computing.  Journal of Open Source Software, 6(62):3060, 2021.
  * [2] Maziar Raissi, Paris Perdikaris, and George 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.
  * [3] Lionel Raff, Ranga Komanduri, Martin Hagan, and Satish Bukkapatnam.  Neural networks in chemical reaction dynamics.  OUP USA, 2012.
  * [4] Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car, and Weinan E.  Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics.  Physical review letters, 120(14):143001, 2018.
  * [5] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [6] Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman.  Universal differential equations for scientific machine learning.  arXiv preprint arXiv:2001.04385, 2020.
  * [7] Samuel Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  Advances in neural information processing systems, 32, 2019.
  * [8] Erik J Bekkers.  B-spline cnns on lie groups.  arXiv preprint arXiv:1909.12057, 2019.
  * [9] Muhammed Shuaibi, Adeesh Kolluru, Abhishek Das, Aditya Grover, Anuroop Sriram, Zachary Ulissi, and C Lawrence Zitnick.  Rotation invariant graph neural networks using spin convolutions.  arXiv preprint arXiv:2106.09575, 2021.
  * [10] Steffen Schotthöfer, Tianbai Xiao, Martin Frank, and Cory D Hauck.  Structure preserving neural networks: A case study in the entropy closure of the boltzmann equation.  In Proceedings of the International Conference on Machine Learning, PMLR, Baltimore, MD, USA, pages 17–23, 2022.
  * [11] Weinan E, Jiequn Han, Linfeng Zhang, et al.  Integrating machine learning with physics-based modeling.  arXiv preprint arXiv:2006.02619, 2020.
  * [12] Keyulu Xu, Mozhi Zhang, Jingling Li, Simon S Du, Ken-ichi Kawarabayashi, and Stefanie Jegelka.  How neural networks extrapolate: From feedforward to graph neural networks.  arXiv preprint arXiv:2009.11848, 2020.
  * [13] Zheyan Shen, Jiashuo Liu, Yue He, Xingxuan Zhang, Renzhe Xu, Han Yu, and Peng Cui.  Towards out-of-distribution generalization: A survey.  arXiv preprint arXiv:2108.13624, 2021.
  * [14] Manuel Torrilhon.  Modeling nonequilibrium gas flow based on moment equations.  Annual review of fluid mechanics, 48:429–458, 2016.
  * [15] Byung C Eu.  Kinetic theory and irreversible thermodynamics.  NASA STI/Recon Technical Report A, 93:24498, 1992.
  * [16] Prabhu Lal Bhatnagar, Eugene P Gross, and Max Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  Physical review, 94(3):511, 1954.
  * [17] Lowell H Holway Jr.  New statistical models for kinetic theory: methods of construction.  Physics of Fluids (1958-1988), 9(9):1658–1673, 1966.
  * [18] EM Shakhov.  Generalization of the krook kinetic relaxation equation.  Fluid Dynamics, 3(5):95–96, 1968.
  * [19] Luc Mieussens and Henning Struchtrup.  Numerical comparison of bhatnagar–gross–krook models with proper prandtl number.  Physics of Fluids, 16(8):2797–2813, 2004.
  * [20] Jeff Haack, Cory Hauck, Christian Klingenberg, Marlies Pirner, and Sandra Warnecke.  A consistent bgk model with velocity-dependent collision frequency for gas mixtures.  Journal of Statistical Physics, 184(3):1–17, 2021.
  * [21] Ching Shen.  Rarefied gas dynamics: fundamentals, simulations and micro flows.  Springer Science & Business Media, 2006.
  * [22] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [23] Tianbai Xiao and Martin Frank.  Using neural networks to accelerate the solution of the boltzmann equation.  Journal of Computational Physics, 443:110521, 2021.
  * [24] Alexander Alekseenko, Robert Martin, and Aihua Wood.  Fast evaluation of the boltzmann collision operator using data driven reduced order models.  Journal of Computational Physics, page 111526, 2022.
  * [25] Sean T Miller, Nathan V Roberts, Stephen D Bond, and Eric C Cyr.  Neural-network based collision operators for the boltzmann equation.  Journal of Computational Physics, page 111541, 2022.
  * [26] Sydney Chapman and Thomas George Cowling.  The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases.  Cambridge University Press, 1970.
  * [27] Graham W Alldredge, Martin Frank, and Cory D Hauck.  A regularized entropy-based moment method for kinetic equations.  SIAM Journal on Applied Mathematics, 79(5):1627–1653, 2019.
  * [28] François Bouchut, François Golse, and Mario Pulvirenti.  Kinetic equations and asymptotic theory.  Elsevier, 2000.
  * [29] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [30] Christopher M Bishop et al.  Neural networks for pattern recognition.  Oxford university press, 1995.
  * [31] Mohammad Amin Nabian and Hadi Meidani.  Physics-informed regularization of deep neural networks.  arXiv preprint arXiv:1810.05547, 2018.
  * [32] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [33] Mike Innes, Alan Edelman, Keno Fischer, Chris Rackauckas, Elliot Saba, Viral B Shah, and Will Tebbutt.  A differentiable programming system to bridge machine learning and scientific computing.  arXiv preprint arXiv:1907.07587, 2019.
  * [34] Hans Babovsky et al.  Die Boltzmann-Gleichung: Modellbildung-Numerik-Anwendungen.  Springer, 1998.
  * [35] C David Levermore.  Moment closure hierarchies for kinetic theories.  Journal of statistical Physics, 83(5):1021–1065, 1996.
  * [36] Raúl E Curto.  Recursiveness, positivity and truncated moment problems.  Houston Journal of Mathematics, 17:603–635, 1991.
  * [37] Michael Junk.  Maximum entropy for reduced moment problems.  Mathematical Models and Methods in Applied Sciences, 10(07):1001–1025, 2000.
  * [38] Tianbai Xiao, Steffen Schotthöfer, and Martin Frank.  Predicting continuum breakdown with deep neural networks.  arXiv preprint arXiv:2203.02933, 2022.
  * [39] T Carleman.  L’intégrale de fourier et questions que s’y rattachent, volume 1 of publications scientifiques de l’institut mittag-leffler.  Almqvist & Wiksells, Uppsala, 1944.
  * [40] Clément Mouhot and Lorenzo Pareschi.  Fast algorithms for computing the boltzmann collision operator.  Mathematics of computation, 75(256):1833–1852, 2006.
  * [41] Kun Xu and Juan-Chen Huang.  An improved unified gas-kinetic scheme and the study of shock structures.  IMA Journal of Applied Mathematics, 76(5):698–711, 2011.
  * [42] Ruifeng Yuan and Lei Wu.  Capturing the influence of intermolecular potential in rarefied gas flows by a kinetic model with velocity-dependent collision frequency.  Journal of Fluid Mechanics, 942, 2022.
  * [43] Tianbai Xiao and Martin Frank.  A stochastic kinetic scheme for multi-scale flow transport with uncertainty quantification.  Journal of Computational Physics, 437:110337, 2021.
  * [44] Tianbai Xiao.  A flux reconstruction kinetic scheme for the boltzmann equation.  Journal of Computational Physics, 447:110689, 2021.
  * [45] Ch Tsitouras, I Th Famelis, and TE Simos.  On modified runge–kutta trees and methods.  Computers & Mathematics with Applications, 62(4):2101–2111, 2011\.
  * [46] Lawrence F Shampine.  Implementation of rosenbrock methods.  ACM Transactions on Mathematical Software (TOMS), 8(2):93–113, 1982\.
  * [47] Uri M Ascher, Steven J Ruuth, and Brian TR Wetton.  Implicit-explicit methods for time-dependent partial differential equations.  SIAM Journal on Numerical Analysis, 32(3):797–823, 1995.
  * [48] Harold Abelson and Gerald Jay Sussman.  Structure and interpretation of computer programs.  The MIT Press, 1996.
  * [49] Harold Meade Mott-Smith.  The solution of the boltzmann equation for a shock wave.  Physical Review, 82(6):885, 1951.
  * [50] Benzi John, Xiao-Jun Gu, and David R Emerson.  Investigation of heat and mass transfer in a lid-driven cavity under nonequilibrium flow conditions.  Numerical Heat Transfer, Part B: Fundamentals, 58(5):287–303, 2010\.
  * [51] Dominique d’Humières.  Multiple–relaxation–time lattice boltzmann models in three dimensions.  Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 360(1792):437–451, 2002.

(a) $\rho-U$ (b) $\rho-T$

Figure 5: phase diagrams of the data generated by the algorithmic sampler and
a single simulation.

(a) $\rho-U$ (b) $\rho-T$

Figure 6: phase diagrams of the data generated by the algorithmic sampler and
multiple simulations. Figure 7: Values of cost function versus optimization
time in the training and test datasets for the homogeneous relaxation problem.

(a) $t=0.5$ (b) $t=1$ (c) $t=2$

Figure 8: Particle distribution functions at different time instants in the
homogeneous relaxation problem (Shakhov model as reference solution).

(a) $t=0$ (b) $t=0.5$ (c) $t=1$

Figure 9: Collision terms at different time instants in the homogeneous
relaxation problem (Shakhov model as reference solution).

(a) Density (b) Energy

Figure 10: Evolution of total density and energy with time in the homogeneous
relaxation problem (Shakhov model as reference solution).

(a) $t=0$ (b) $t=0.5$ (c) $t=1$

Figure 11: Equilibrium functions at different time instants in the homogeneous
relaxation problem (Shakhov model as reference solution).

(a) $t=0.5$ (b) $t=1$ (c) $t=2$

Figure 12: Particle distribution functions at different time instants in the
homogeneous relaxation problem ($\nu$-BGK model as reference solution).

(a) $t=0$ (b) $t=0.5$ (c) $t=1$

Figure 13: Collision terms at different time instants in the homogeneous
relaxation problem ($\nu$-BGK model as reference solution).

(a) $t=0$ (b) $t=0.5$ (c) $t=1$

Figure 14: Relaxation frequencies at different time instants in the
homogeneous relaxation problem ($\nu$-BGK model as reference solution).

(a) Density (b) Energy

Figure 15: Evolution of total density and energy with time in the homogeneous
relaxation problem ($\nu$-BGK model as reference solution).

(a) $t=0.2$ (b) $t=0.5$ (c) $t=1$ (d) $t=2$

(e) $t=4$ (f) $t=6$

Figure 16: Particle distribution functions at different time instants in the
homogeneous relaxation problem (full Boltzmann model as reference solution).

(a) $t=0$ (b) $t=1$ (c) $t=2$

Figure 17: Collision terms at different time instants in the homogeneous
relaxation problem (full Boltzmann model as reference solution).

(a) $t=0$ (b) $t=1$ (c) $t=2$

Figure 18: Equilibrium (dark colors) and relaxation frequencies (light colors)
at different time instants in the homogeneous relaxation problem (full
Boltzmann model as reference solution). Figure 19: Evolution of entropy with
time in the homogeneous relaxation problem (full Boltzmann model as reference
solution).

(a) Density (b) Velocity (c) Temperature (d) Viscous stress

(e) Heat flux (f) Entropy

Figure 20: Profiles of macroscopic variables at $\rm Ma=2$ in the normal shock
structure.

(a) Density (b) Velocity (c) Temperature (d) Viscous stress

(e) Heat flux (f) Entropy

Figure 21: Profiles of macroscopic variables at $\rm Ma=3$ in the normal shock
structure.

(a) UBE (b) BGK (c) UBE (d) BGK

Figure 22: Contours of particle distribution function (first row) and
collision term (second row) at $\rm Ma=2$ in the normal shock structure.

(a) UBE (b) BGK (c) UBE (d) BGK

Figure 23: Contours of particle distribution function (first row) and
collision term (second row) at $\rm Ma=3$ in the normal shock structure.

(a) Density with streamlines (b) Temperature with heat flux vectors

Figure 24: Contours of density with streamlines and temperature with heat flux
vectors inside the cavity.

(a) $U$-velocity along vertical central line (b) $V$-velocity along horizontal
central line

Figure 25: Velocity profiles along vertical and horizontal central lines of
the cavity. The results are normalized with the wall speed $U_{w}$.

(a) UBE (b) BGK (c) UBE (d) BGK

Figure 26: Particle distribution functions along vertical (first row) and
horizontal (second row) central lines of the cavity.

(a) UBE (b) BGK (c) UBE (d) BGK

Figure 27: Collision terms along vertical (first row) and horizontal (second
row) central lines of the cavity.