Model Order Reduction for the 1D Boltzmann-BGK Equation: Identifying
  Intrinsic Variables Using Neural Networks
Kinetic equations are crucial for modeling non-equilibrium phenomena, but
their computational complexity is a challenge. This paper presents a
data-driven approach using reduced order models (ROM) to efficiently model
non-equilibrium flows in kinetic equations by comparing two ROM approaches:
Proper Orthogonal Decomposition (POD) and autoencoder neural networks (AE).
While AE initially demonstrate higher accuracy, POD's precision improves as
more modes are considered. Notably, our work recognizes that the classical
POD-MOR approach, although capable of accurately representing the non-linear
solution manifold of the kinetic equation, may not provide a parsimonious model
of the data due to the inherently non-linear nature of the data manifold. We
demonstrate how AEs are used in finding the intrinsic dimension of a system and
to allow correlating the intrinsic quantities with macroscopic quantities that
have a physical interpretation.
Julian Koellermeier, Philipp Krah, Julius Reiss, Zachary Schellin
Model Order Reduction for the 1D Boltzmann-BGK Equation: Identifying Intrinsic Variables Using Neural Networks

###### Abstract

Kinetic equations are crucial for modeling non-equilibrium phenomena, but their computational complexity is a challenge. This paper presents a data-driven approach using reduced order models (ROM) to efficiently model non-equilibrium flows in kinetic equations by comparing two ROM approaches: Proper Orthogonal Decomposition (POD) and autoencoder neural networks (AE). While AE initially demonstrate higher accuracy, POD's precision improves as more modes are considered. Notably, our work recognizes that the classical POD-MOR approach, although capable of accurately representing the non-linear solution manifold of the kinetic equation, may not provide a parsimonious model of the data due to the inherently non-linear nature of the data manifold. We demonstrate how AEs are used in finding the intrinsic dimension of a system and to allow correlating the intrinsic quantities with macroscopic quantities that have a physical interpretation.

**Keywords**: Model order reduction, data-driven methods, kinetic equations, neural autoencoder networks, proper orthogonal decomposition, Sod shock tube, Boltzmann-BGK

## 1 Introduction

Kinetic equations are widely used in science and engineering [20, 25, 26, 33]. They allow the modeling of deviations from an equilibrium model which is given by an underlying macroscopic equation like the Euler equations, providing detailed insight into fundamental physical processes [35]. However, kinetic equations are often characterized by a large dimensional phase space, making them computationally expensive to solve and sometimes even unfeasible for realistic applications [35].

Investing in solving kinetic equations is only beneficial if large deviations from equilibrium are present [35]. Striking a balance between a fast but inaccurate equilibrium solver and a slow but accurate non-equilibrium solver remains an open challenge. Our work aims to address this challenge by providing a proof-of-concept for a data-driven solution to efficient modeling of flows in different non-equilibrium regimes.

In the field of non-equilibrium gas flows, several standard methods to discretize the high dimensional phase space exist. Particle-based Monte Carlo methods are only tractable in the free flight regime and out of scope in the transition regime of moderate non-equilibrium unless special techniques are used [5]. The straightforward Discrete Velocity Method (DVM) uses a pointwise discretization of the velocity space, potentially leading to a large number of equations [30]. Specially tailored moment models are based on the expansion of the particle distribution function and lead to a set of extended fluid dynamical equations [35]. However, it is by no means clear a-priori how many equations are sufficient and which variables are optimal [19, 34].
To tackle the computational complexity of kinetic equations, recently, reduced order models (ROM) have been introduced, enabling reductions in computational complexity by orders of magnitudes [1, 7, 8, 9]. Two different approaches have been followed in the literature. The classical offline-online decomposition as used by [1] involves a two-stage procedure. In the offline stage, the full order model (FOM) is assessed to create a database, which is then utilized to generate a data-dependent basis through proper orthogonal decomposition (POD). This basis allows an efficient description of the FOM on a low-dimensional linear subspace during the online phase. On the other hand, the online adaptive basis method called dynamic low-rank approximation [13] constructs the low dimensional linear basis during the online phase itself, eliminating the need to evaluate the expensive FOM. It has been successfully applied to kinetic equations in the works by Einkemmer et al. [7, 8, 9]. However, the additional complexity of updating the basis during the evolution makes it less online efficient than the classical offline-online approach shown in [15] for a shallow water moment model.

In this work, we adopt the same offline strategies as in [1]. Specifically, we sample data for a classical test case called Sod shock tube using a discrete velocity model as our FOM and compare the compression of the linear reduced subspace created by POD with a non-linear description provided by neural autoencoder networks. Neural networks, based on the universal approximation theorem [29], allow for the approximation of a wide range of function classes and appear promising in identifying the intrinsic dimension of a system. However, the non-linear relation between macroscopic model equations and the discrete velocity model hinders the determination of these dimensions using linear reduction methods like the POD.

This paper aims to utilize these data-driven model reduction techniques to reduce the number of describing variables and equations and determine how many and which variables are useful in specific test cases. For the non-vanishing Knudsen number, we expect to need more non-equilibrium variables with corresponding balance laws, while in the limit of vanishing Knudsen number, we expect to recover the Euler equations, given by conservation laws for mass, momentum, and energy. To the knowledge of the authors, this is the first paper aiming to bridge the gap between equilibrium and non-equilibrium flows using neural networks in this way.

The long-term objective of this line of work is to enable dynamically adapting the model by varying the number of variables during the online phase, paving the way for more efficient and accurate model adaptive simulations of kinetic equations.

The organization of the paper is as follows: In Section 2, we introduce the 1D model problem and the reference data used for model reduction. Section 3 describes the two model reduction techniques used in this study: Proper Orthogonal Decomposition (POD) and Autoencoder Networks. The results are presented in Section 4, and the paper concludes with a summary in Section 5.

## 2 The Boltzmann-BGK Model and Data

This paper considers a proof-of-concept of using reduced models for the solution approximation of the 1D Boltzmann-BGK equation [2] for monoatomic, ideal gases

\[\partial_{t}f+c\partial_{x}f=\frac{1}{\tau}(f_{M}-f), \tag{1}\]

which is a potentially high-dimensional equation for the unknown probability density function \(f(t,x,c)\), where \(t\in\mathbb{R}^{+}\) is the time, \(x\in\mathbb{R}\) is the spatial variable, and \(c\in\mathbb{R}\) the microscopic particle velocity. For simplicity we consider the one-dimensional case in this paper, but the results can be extended to the multi-dimensional case.

Computing solutions and generating data of the Boltzmann-BGK model is essential for industrial and scientific applications, but often so computationally prohibitive that a large number of test cases is not feasible. To reduce time and cost during the data generating process, experiments or numerical simulations can be replaced by reduced-order models (ROMs).
For standard continuum flows the widely-used Euler equations can be applied, but more rarefied regimes require different extended fluid dynamical models. Rarefaction levels are distinguished with the help of the Knudsen number \(\mathrm{Kn}\) defined by the ratio of the mean free path length of the particles \(\lambda\) over a reference length \(l\) as

\[\mathrm{Kn}=\frac{\lambda}{l}. \tag{2}\]

The right-hand side of the BGK collision operator (1) models the relaxation with relaxation time \(\tau\in\mathbb{R}^{+}\) towards the equilibrium Maxwellian distribution \(f_{M}(t,x,c)\) given by

\[f_{M}(t,x,c)=\frac{\rho(t,x)}{(2\pi RT(t,x))^{\frac{3}{2}}}\exp\left(-\frac{(c -u(t,x))^{2}}{2RT(t,x)}\right), \tag{3}\]

where \(\rho(t,x)\), \(v(t,x)\) and \(T(t,x)\) are density, bulk velocity, and temperature of the flow, respectively. \(R\) is the universal gas constant. In this work, we consider the relaxation time \(\tau\) a parameter and set it equal to the Knudsen number, \(\tau=\mathrm{Kn}\), however, the relaxation time can also be changed, e.g., to depend on the gas density and temperature in addition.

For practical computations, we consider macroscopic moments of the distribution function, which are given by multiplying the distribution function with the co-called collision invariants \((1,c,\frac{1}{2}c^{2})\) and integrating in velocity space

\[\rho(t,x) = \int f(t,x,c)\,\mathrm{d}c, \tag{4}\] \[\rho(t,x)u(t,x) = \int cf(t,x,c)\,\mathrm{d}c,\] (5) \[E(t,x) = \int\frac{1}{2}c^{2}f(t,x,c)\,\mathrm{d}c, \tag{6}\]

where \(E\) denotes the total energy. The temperature \(T(t,x)\) and the pressure \(p(t,x)\) can be obtained by

\[T(t,x)=\frac{2E(t,x)}{3\rho(t,x)}-\frac{u(t,x)^{2}}{3}\quad\text{and}\quad p(t,x)=\rho(t,x)T(t,x). \tag{7}\]

Figure 1 illustrates the relation between the macroscopic moments and the distribution function \(f(t,x,c)\) at a certain position in time and space. The density \(\rho(t,x)\) is the integral of the distribution function, which is centered around the macroscopic velocity \(u(t,x)\), and the mean deviation is related to the temperature \(T(t,x)\).

The Boltzmann-BGK equation (1) is in equilibrium when \(f=f_{M}\). Multiplying the equilibrium solution with the collision invariants and integrating in velocity space, one finds the Euler equations of classical gas dynamics

\[\partial_{t}\rho+\partial_{x}(\rho u)=0, \tag{8}\] \[\partial_{t}(\rho u)+\partial_{x}(\rho u^{2}+p)=0,\] (9) \[\partial_{t}E+\partial_{x}(u(E+p))=0, \tag{10}\]

Figure 1: Illustration of the macroscopic moments corresponding to an example distribution function.
which are conservation laws for mass, momentum, and energy, respectively.

For distribution functions further away from equilibrium, for example due to a larger relaxation time \(\tau\) and a significantly large Knudsen number Kn, the Euler equations do not give accurate results. In this case, additional equations can be used, which are derived by the so-called method of moments [14, 35]. This effectively leads to an extended set of equations, called moment model. It is possible to preserve important properties like hyperbolicity with moment models [11, 17]. The additional equations (for example for the heat flux and higher-order moments) add complexity, but allow for more accurate solutions [18, 34]. However, it is often unclear a-priori, how many equations are needed for an efficiently accurate and computationally feasible solution. In this work, we aim to give a proof-of-concept for a data-based identification of the necessary number of variables, called the intrinsic physical dimension.

### Sod shock tube test case and reference data

Sod's shock tube is a well-established test case in the field of rarefied gases [19]. It uses discontinuous initial conditions based on equilibrium values

\[\begin{cases}(\rho_{L},u_{L},p_{L})=(1,0,1)&\text{if }x<0.5,\\ (\rho_{R},u_{R},p_{R})=(0.125,0,0.1)&\text{if }x>0.5,\end{cases} \tag{11}\]

corresponding to a jump in density and pressure at \(x=0.5\) due to a diaphragm at that position, which is removed at time \(t=0\).

The problem setup at \(t=0\) is shown in fig. 2, which is split into two regions left and right of the diaphragm.

For the generation of reference data we employ a discrete velocity method (DVM) [27], which uses a pointwise microscopic velocity space discretization

\[\partial_{t}f_{k}(t,x)=-(c_{k})\partial_{x}f_{k}(t,x)+\frac{1}{\tau}\left(M_{ fk}(t,x)-f_{k}(t,x)\right), \tag{12}\]

where a uniform grid in velocity space is considered with \(c_{k}=k\Delta c\) to discretize the distribution function \(f_{k}(t,x)=f(t,x,c_{k})\), for some \(k\in\mathbb{Z}\). After a subsequent discretization in space, the DVM eq. (12) leads to a coupled ODE system in time than can be solved with standard methods.

For the numerical reference data, we use \(N_{x}=200\) spatial points \(x_{k}\in[0,1]\), \(N_{c}=40\) discrete velocities \(c_{j}\in[-10,10]\) and \(N_{t}=25\) time steps \(t^{n}\in[0,0.12]\) summarized in table 1. It is possible to choose another range for the discrete velocity points, but in typical applications the range of the bulk velocity is not known, such that one has to include a safety margin. We therefore chose the domain \([-10,10]\). The goal of the model order reduction is now to reduce the complexity of the computation using lower dimensional models. For that matter, it is not relevant what the actual error of the numerical reference data is or if is fully converged. It is fair to say that a full reference solution might easily take into account more spatial points, time steps, and discrete velocities, which makes it even more necessary to reduce the complexity.

Figure 2: Problem setup for the 1D Sod shock tube. A diaphragm at the center is initially separating the domain in two regions, where initial conditions for density \(\rho\), macroscopic velocity \(u\), and pressure \(p\) are indicated.
For the model order reduction later, we consider two different Knudsen numbers for Sod's shock-tube test case: \(\mathrm{Kn}=0.00001\) for a small Knudsen number in the hydrodynamic regime and \(\mathrm{Kn}=0.01\) for a relatively large Knudsen number in the rarefied regime.

To understand the behavior of the reference solutions for non-vanishing Knudsen numbers, we first describe the solution for vanishing Knudsen number in equilibrium, i.e., \(\mathrm{Kn}=0\), which can be obtained using the method of characteristics and the Rankine Hugoniot jump conditions connecting the states before and after the shocks [24].

Starting from the initial condition in fig. 3a, the solution evolves for \(t>0\) and five regions are formed that are depicted in fig. 3b [32]. A rarefaction wave is moving to the left between \(x_{1}\) and \(x_{2}\). The contact discontinuity is located at \(x_{3}\), where the macroscopic velocity \(u\) and the pressure \(p\) are continuous in contrast to the density \(\rho\) and the energy \(E\). \(x_{4}\) is the position of the shock wave.

In non-equilibrium, i.e., for solutions evolving with Knudsen numbers \(\mathrm{Kn}>0\), the solution does not have discontinuities due to the finite relaxation time \(\tau\). Figure 3c shows the reference solutions \(f(t,x,c)\) at \(t_{0}=0\), \(t_{1}=0.06\) and \(t_{3}=0.12\) for the two levels of rarefaction considered in this paper: \(\mathrm{Kn}=0.00001\) and \(\mathrm{Kn}=0.01\). Increasing the Knudsen number leads to a smoother transition from region 1 to region 5 with a less pronounced shock front.

## 3 Methods

In this section we present two common methods used for reducing the dimensionality of high dimensional data: (1) the proper orthogonal decomposition (POD) and (2) neural autoencoder networks (AE). The methods will be used to parameterise the high dimensional data stemming from the DVM simulation, using a linear mapping in case of the POD and a non-linear mapping in case of AE. Although the classical POD-MOR approach shows that linear mappings are sufficient to describe the non-linear solution manifold of the BGK equation to a good accuracy [1], it is in general not sufficient to determine a parsimonious model of the full model data, since the data manifold can be non-linear. Here, neural autoencoder networks can be used as they are capable to find the intrinsic dimension of a system.

### Proper orthogonal decomposition

The _proper orthogonal decomposition_[31] (POD) approximates the data with help of dyadic pairs:

\[f(x,t,c_{i})\approx\sum_{k=1}^{r}\hat{f}_{k}(x,t)\psi_{k}(c_{i})\qquad\text{ for }r\ll N_{c}\,. \tag{13}\]

The pairs \(\{(\hat{f}_{k}(x,t),\psi_{k}(c_{i}))\}_{k=1,\ldots,r}\) are the structures in the data that contain the most energy and they are chosen to minimize the gap between the data and the reconstruction eq. (13). In the following, \(\psi_{k}(c_{i})\) are termed POD-modes and \(\hat{f}_{k}(x,t)\) the corresponding reduced variables.

For notation, we define \(f^{(i)}(x,t)=f(x,t,c_{i})\) and the vectors \(f(x,t)=(f^{(1)}(x,t),\ldots,f^{(N_{c})}(x,t))\) and \(\psi_{k}=(\psi_{k}(c_{1}),\ldots,\psi_{k}(c_{N_{c}}))\). The proper orthogonal decomposition computes the solution of the minimization problem:

\[\min_{\psi_{k}}\|f(x,t)-\sum_{k=1}^{r}\langle f(x,t),\psi_{k}\rangle\psi_{k}\| _{2}^{2}\quad\text{such that}\quad\langle\psi_{k},\psi_{l}\rangle=\delta_{kl}. \tag{14}\]

\begin{table}
\begin{tabular}{c c c c} \hline Variable & Number of nodes \(i\) & Domain extension & Step size (uniform) \\ \hline \(x\) & 200 & [0, 1] & 0.005 \\ \(c\) & 40 & [-10,10] & \(\approx\) 0.51282051 \\ \(t\) & 25 & [0,0.12] & 0.005 \\ \hline \end{tabular}
\end{table}
Table 1: Problem setup for the Boltzmann-BGK model in Sod’s shock tube.
Figure 3: Sod shock tube and reference data. Initial conditions (a); equilibrium solution (b); Reference solutions in rarefied and hydrodynamic regime (c); Macroscopic quantities at \(t=0.12s\) (d).
Technically one can solve this optimization problem using a singular value decomposition (SVD) of the so-called snapshot matrix [22]:

\[\mathbf{F}=\begin{bmatrix}f^{(1)}(x_{1},t_{1})&\cdots&f^{(N_{c})}(x_{1},t_{1})\\ \vdots&\ddots&\vdots\\ f^{(1)}(x_{N_{a}},t_{1})&\cdots&f^{(N_{c})}(x_{N_{x}},t_{1})\\ f^{(1)}(x_{1},t_{2})&\cdots&f^{(N_{c})}(x_{1},t_{2})\\ \vdots&\ddots&\vdots\\ f^{(1)}(x_{N_{x}},t_{N_{t}})&\cdots&f^{(N_{c})}(x_{m},t_{N_{t}})\end{bmatrix} \in\mathbb{R}^{(N_{x}N_{t})\times N_{c}}\,. \tag{15}\]

The snapshot matrix collects the time and space discrete distribution function in its columns. Each column holds the time-spatial values for a different discrete velocity. Performing an SVD factorizes \(\mathbf{F}\) as

\[\mathbf{F}=\mathbf{\Phi}\mathbf{\Sigma}\mathbf{\Psi}^{\mathbf{T}}, \tag{16}\]

with diagonal matrix \(\mathbf{\Sigma}=\text{diag}(\sigma_{1},\ldots,\sigma_{m})\), \(m=\text{min}(N_{x}N_{t},N_{c})\), containing the singular values \(\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{m}\geq 0\) and \(\mathbf{\Phi}\in\mathbb{R}^{N_{x}N_{t}\times m},\mathbf{\Psi}\in\mathbb{R}^{N _{c}\times m}\) are orthogonal matrices containing the left and right singular vectors, respectively.

The first \(r\) columns of the truncated \(\mathbf{\Psi}_{r}=[\psi_{1},\ldots,\psi_{r}]\in\mathbb{R}^{N_{c}\times r}\) contain the POD modes in eq. (13). Together with \(\mathbf{\Sigma}_{r}=\text{diag}(\sigma_{1},\ldots,\sigma_{r})\) and the \(r\) leading left singular vectors \(\mathbf{\Phi}_{r}\in\mathbb{R}^{(N_{x}N_{t})\times r}\) they yield the rank \(r\)-term approximation \(\mathbf{F}_{r}\) of the snapshot matrix given by

\[\mathbf{F}_{r}:=\mathbf{\Phi}_{r}\mathbf{\Sigma}_{r}\mathbf{\Psi}_{r}^{T}. \tag{17}\]

According to the Eckart-Young-Mirsky theorem [6, 28]\(\mathbf{F}_{r}\) is the best rank \(r\) approximation and the resulting error in the Frobenius norm is rigorously computed from the trailing singular values

\[\|\mathbf{F}-\mathbf{F}_{r}\|_{\mathrm{F}}^{2}=\sum_{k=r+1}^{m}\sigma_{k}^{2}. \tag{18}\]

A common choice for \(r\) is to truncate after a certain energy percentage is reached in the reduced system compared to the full system:

\[E_{\text{cum}}=\frac{\|\mathbf{F}_{r}\|_{\mathrm{F}}}{\|\mathbf{F}\|_{ \mathrm{F}}}=\frac{\sum_{k=1}^{r}\sigma_{k}^{2}}{\sum_{k=1}^{m}\sigma_{k}^{2}}. \tag{19}\]

In this paper, POD is used to compare with the autoencoder from the next section.

### Autoencoders

Neural networks, particularly autoencoder networks, have become widely used tools for dimension reduction [36]. A comprehensive introduction to autoencoder networks can be found in [12]. Here we only summarize briefly the common idea of autoencoder networks and give the specific details of our implementation.

An autoencoder aims to reproduce the input data while compressing it through an information bottleneck. It consists of two main components:

* The encoder, denoted as \(g_{\text{enc}}\), maps the input data \(f\) to points \(\hat{f}\) in a lower-dimensional latent space: \(g_{\text{enc}}\colon\mathbb{R}^{M}\to\mathbb{R}^{r},f\mapsto\hat{f}=g_{\text{ enc}}(f)\), \(r\ll M\).
* The decoder, denoted as \(g_{\text{dec}}\), reconstructs the input space from the latent representation: \(g_{\text{dec}}\colon\mathbb{R}^{r}\to\mathbb{R}^{M},\hat{f}\mapsto g_{\text{ dec}}(\hat{f})=\tilde{f}\), \(r\ll M\).

Note that the dimension of the latent space is denoted by \(r\), to match the rank of the POD approximation.

The autoencoder is defined as the composition of both parts: \(\tilde{f}=g_{\text{dec}}\left(g_{\text{enc}}\left(f\right)\right)\). For our purpose, we identify the discrete velocity space as the input dimension \(M=N_{c}\). Thus, the autoencoder maps each time-spatial value of the distribution function \(f(x,t)\in\mathbb{R}^{N_{c}}\) onto a
smaller latent space \(\hat{f}(x,t)\in\mathbb{R}^{r}\), which parameterizes the necessary physical information of the system.

The goal of the optimization procedure is to determine \(g_{\text{dec}}\) and \(g_{\text{enc}}\) such that the reconstruction error over a set of training/testing data contained in \(\mathbf{F}\) is minimized. The reconstruction error is defined as:

\[\mathcal{L}=\frac{1}{N_{c}}\|f(x,t)-g_{\text{dec}}\circ g_{\text{enc}}\circ f(x,t)\|_{2}^{2}\,. \tag{20}\]

The reconstruction error is the sum of the two-norm of the discrete velocities vector of the difference between the input data \(f\) and the reconstructed data \(\tilde{f}\) that has been squeezed through the informational bottleneck. The assumption is, that if the original data can be represented well while the information went through a smaller latent space, there exists a physical law in the latent space that describes the system sufficiently. The intrinsic latent dimension \(r=p^{*}\) which is sufficient to describe the data is then called the _intrinsic physical dimension_ similar to the intrinsic dimension defined in [23]. Such a reduced model is then termed parsimonious because it explains the data with a minimum number of variables.

In the training procedure, the functions \(g_{\text{enc}}\) and \(g_{\text{dec}}\) are determined by trainable parameters of the network, referred to as weights and biases. The networks are constructed using a composition of layers \(g_{\text{enc}}=L_{1}\circ L_{2}\circ\cdots\circ L_{N}\). Typically, each layer \(L_{n}\colon\mathbb{R}^{i}\to\mathbb{R}^{o}\) in the network consists of an affine linear mapping \(\mathbf{x}\mapsto h_{n}(\mathbf{W}_{n}x+b_{n})\), where \(\mathbf{W}_{n}\in\mathbb{R}^{o,i}\) represent the weights, \(\mathbf{b}_{n}\in\mathbb{R}^{o}\) denote the biases, and \(h_{n}\) are predefined non-linear functions. The configuration of the input and output dimensions \(i\) and \(o\) for each layer, the choice of activation function, and the number of layers collectively determine the architecture of the network. The choice of these so-called hyper-parameters is often difficult and a matter of trial and error.

ArchitectureIn our studies we have exploited _fully connected neural autoencoder networks_ (FCNN) and convolutional autoencoder networks (CNN). However, in this manuscript we restrict ourselves to the results of the fully connected network, since it gave structurally the best results. We have studied a variety of different activation functions, hidden layers, batch sizes and depths of the network. The best results concerning the validation error and acceptable training time where obtained by the network defined in table 2. A comprehensive study of the parameter optimization is attached to the manuscript in section 5.

TrainingBefore the training we initialize the weights of the network using the standard initialization implemented in pytorch. Thus the weights are randomly uniform distributed between \(m^{-1/2}\) and \(m^{1/2}\) with \(m\) being the number of input nodes in the layer. Our network is trained by splitting the data consisting of \(N_{x}\times N_{t}\) samples in a testing and training set with a 80/20 split over 3000 epochs using a batch size of 4. In each epoch the network is updated using the Adam optimizer with a learning rate of \(10^{-5}\). More information about hyperparameters and training of the network can be found in the appendix section 5.

\begin{table}
\begin{tabular}{l c c} \hline \hline
**Parameter** & hydrodynamic & rarefied \\ \hline Layer Sizes & [\(N_{c}\), 30, \(r\), 30, \(N_{c}\)] & [\(N_{c}\), 40, \(r\), 40, \(N_{c}\)] \\ Activation function & ELU & ReLU \\ Loss function & MSE eq. (20) & MSE eq. (20) \\ Optimizer & Adam & Adam \\ Learning rate & \(10^{-5}\) & \(10^{-5}\) \\ Epochs & 3000 & 3000 \\ Batch size & 4 & 4 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Hyper-parameters of Fully Connected Autoencoder Network (FCNN) for the hydrodynamic and rarefied regime.
Results

We reconstruct the full order model (FOM) solution with the help of POD and an autoencoder, the FCNN, for which the selection of hyperparameters and the training are described in the previous section. Note that we apply both model reduction techniques to reconstruct both the rarefied reference data and the hydrodynamic reference data. The goal is to later determine the intrinsic dimension of the data for both cases. We therefore compare the two dimension reduction techniques by means of different measures. The intrinsic variables obtained from POD and the FCNN will be referred to as \(\mathbf{h}\) and \(\mathbf{r}\), where the former describes the intrinsic variables when reducing the hydrodynamic data and the latter when reducing the rarefied data.

For the purpose of comparing the results we define the \(L_{2}\)-error

\[\mathcal{E}_{\text{rel}}=\frac{\|\mathbf{F}-\tilde{\mathbf{F}}\|_{2}}{\| \mathbf{F}\|_{2}}\,, \tag{21}\]

where \(\mathbf{F}\) the reference data is given in eq. (15) and the reconstructed data \(\tilde{\mathbf{F}}\) is either \(\mathbf{F}_{r}\) in case of the POD or the FCNN predictions with \(r\) latent variables for every \((x,t,c)\) in the data set.

### Singular value decay of reference data

As a first step, we perform a POD with the hydrodynamic data and with the rarefied data. The obtained singular values \(\sigma\), as well as the cumulative energy (cusum-e) defined in eq. (19), are shown in fig. 4. As expected, more modes are necessary in the rarefied regime compared to the hydrodynamic regime. With a total of \(p^{\text{POD}}=4\) intrinsic variables, a cumulative energy of over 99% can be achieved for the hydrodynamic regime. The cumulative energy of the singular values of the rarefied regime only reaches above 99% with \(p^{\text{POD}}=6\) singular values.

For the POD we define the _intrinsic dimension_\(p^{\text{POD}}\) as the smallest truncation rank \(r\) of the reduced system, at which the cumulative energy \(E_{\text{cum}}\) defined in eq. (19) reaches 99%. Although this choice is arbitrary it is a common practice in classical MOR to truncate eq. (13), whenever 99% of the cumulative energy is reached.

In fig. 4, we further see that the rate at which the singular values drop is approximately exponential in both regimes, which has been also observed by [1]. Consequently, a rapid decay of the Kolmogorov N-width is indicated. Note that the singular value decay is similar for both domains but not exactly the same, thus leading to an expected increase in the number of intrinsic variables in the rarefied regime necessary to achieve similar \(L_{2}\)-errors.

It is important to note that the parameter \(p^{\text{POD}}\) is not expected to precisely match the actual intrinsic dimension \(p^{*}\) of the solution manifold. The intrinsic dimension represents

Figure 4: Singular value decay \(\sigma\) and cumulative energy increase for the number of singular variables \(k\) in the hydrodynamic regime (a) and in the rarefied regime (b). A black cross marker corresponds to over 99% cumulative energy.
the minimum number of variables required to accurately describe the system's exact solution manifold. This discrepancy arises because the solution manifold is fundamentally nonlinear, making it challenging to adequately capture with a parsimonious linear approximation.

For the FCNN, the intrinsic dimension is defined as the smallest number of intrinsic variables that minimizes the error. In well-trained models, the FCNN's intrinsic dimension should ideally align with \(p^{*}\).

From a fluid mechanics perspective, the hydrodynamic regime theoretically requires only \(p^{*}=3\) intrinsic variables. This is because near-equilibrium flows in this regime can be effectively characterized by three macroscopic quantities: density \(\rho\), macroscopic velocity \(u\), and total energy \(E\), as outlined in eq. (8)-eq. (10), see also [1, 19].

Conversely, the rarefied regime demands a larger intrinsic dimension, denoted as \(p^{*}\). This is due to the need for more than only the equilibrium Maxwellian distribution function to describe the microscopic velocities adequately. Therefore, we initially set \(p^{*}=3\) intrinsic variables (\(\mathbf{h}\)) for the FCNN in the hydrodynamic case and choose \(p^{*}=5\) intrinsic variables (\(\mathbf{r}\)) for the rarefied regime. This choice aligns with extended fluid dynamic models as described in [19, 35].

We note that each FCNN with different latent space dimension \(r\) needs to be trained separately. This is different from the POD, where the decomposition is only performed once. Thereafter the approximation quality is given by the truncation rank \(r\).

### Variation of the number of intrinsic variables

The variation of the number of intrinsic variables \(r\) in fig. 5 sheds light on the performance of the autoencoder with different bottleneck layer sizes. In the case of the POD \(r\) is the truncation rank of the decomposition eq. (13) and the latent space dimension in case of the FCNN. To this end, \(r\) is varied for both the POD and the FCNN over \(r\in\{1,2,3,4,8,16,32\}\) for the hydrodynamic case and over \(r\in\{1,2,4,5,8,16,32\}\) for the rarefied case. We note that the loss of information when applying POD goes exponentially to zero with increasing \(r\), which is not surprising when consulting the _Eckard-Young Theorem_[6]. Note that the FCNN is retrained for each different \(r\). By changing \(r\), i.e. widening the bottleneck layer, a gain or loss of capacity occurs that can be connected to stability during training.

Both for the hydrodynamic and the rarefied regime, POD initially yields a larger error than the FCNN for small number of intrinsic variables \(r\). Not surprisingly, the POD accuracy increases with the number of singular values taken into account until the error reaches machine precision. The FCNN error decreases as well and then reaches a plateau, with a

Figure 5: The \(L_{2}\)-Error over the variation of the latent space dimension/truncation rank \(r\) using FCNN/POD for the hydrodynamic regime (left) and the rarefied regime (right).
typical remaining error due to the network architecture and training. For the previously identified values \(p^{*}=3\) in the hydrodynamic case and \(p^{*}=5\) in the rarefied case, the FCNN results in a more accurate approximation than the POD.

We note that when testing the FCNN against POD and fixing \(r\) the FCNN is limited by the estimation error of the training and performs under its abilities. However, POD uses five to six times more parameters than the FCNN while the deterministic character enables POD to achieve any possible accuracy, which was not observed with the neural network.

In the following, we consider the intrinsic variables with constant values \(p^{*}=3\) in the hydrodynamic case and \(p^{*}=5\) in the rarefied case.

This leads to the number of trainable parameters of the POD and the FCNN shown in table 3. We can see that the FCNN achieves a relatively small error with a small number of parameters in comparison with the POD for this choice of the number of intrinsic variables.

Next, a qualitative analysis with the actual reconstructions is presented. From computations of the \(L_{2}\)-error over time \(t\), which are not shown due to conciseness, it became clear that the time step that contributes most to the error is the last time step in case of POD, while the FCNN distributes the error more evenly over all time steps.

### Reconstruction quality

The reconstructed solutions compared to the full order model (FOM) at the final time step \(t=0.12s\) are given in fig. 6. Because of the small overall errors indicated in table 3 both the POD and the FCNN reproduce the FOM solution without any visible differences at first sight.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Algorithm & \multicolumn{2}{c}{Parameters \(\theta\)} & \multicolumn{2}{c}{Int. variables \(p\)} & \multicolumn{2}{c}{\(L_{2}\)-error} \\  & **H** & **R** & **H** & **R** & **H** & **R** \\ \hline POD & 15129 & 25225 & 3 & 5 & 0.0205 & 0.0087 \\ FCNN & 2683 & 3725 & 3 & 5 & 0.0008 & 0.0009 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Amount of parameters used to reconstruct \(f\), the number of intrinsic variables \(p\) and the corresponding \(L_{2}\)-Error for POD and FCNN, both for hydrodynamic (**H**) and rarefied (**R**) regimes.

Figure 6: Comparison of the FOM solutions \(f\) (left column) with reconstructions \(\hat{f}\) obtained from POD (middle column) and FCNN (right column) at end time \(t=0.12s\) for \(x\in[0.375,0.75]\) in the hydrodynamic regime (top row) and in the rarefied regime (bottom row).
As shown in fig. 7, the more subtle information loss from the model reduction can unfold in actual differences in the macroscopic quantities \(\rho\), \(\rho u\) and \(E\). Overall, fig. 7 shows that the errors are larger for the hydrodynamic regime (top row), most notably for the momentum and energy of the POD model close to the contact discontinuity. However, the position of the shock is well approximated. In contrast, the FCNN model yields a very good agreement in the hydrodynamic case. For the rarefied regime, both models approximate the FOM solution very well. The lack of sharp shock structures in the full model and the increased intrinsic dimension \(p^{*}=5\) combined seem to notably influence the accuracy.

### Conservation properties

The physical consistency of the reduced \(\tilde{f}\), in terms of conservation of mass, momentum, and energy, is a critical criterion for its validity. Hence, conservation properties are analyzed in the following. We note that conservation of mass, momentum, and energy is not directly built-in by means of a specifically tailored loss function. Even though we can expect to recover some conservation properties as they are implicitly built into the numerical reference data. The investigation of different loss functions to improve upon this is left for future work.

We investigate the conservation properties by means of the derivative of the cell-averaged conserved quantities mass, momentum, and total energy, defined exemplary as

\[\frac{\mathrm{d}}{\mathrm{d}t}\int\rho(x,t)\,\mathrm{d}x\Delta t=\overline{ \rho}, \tag{22}\]

for the mass. Note that a derivative \(\overline{\rho}=0\) denotes conservation of mass, for example. We expect the conservation of mass and total energy, while the momentum increases with a constant rate, due to the boundary conditions of the test case, featuring a larger pressure on the left-hand side of the domain.

Figure 8 shows the evolution of the derivatives of mass, momentum, and total energy as a function of time for the hydrodynamic regime (top row) and the rarefied regime (bottom

Figure 7: Comparison of FOM macroscopic quantities \(\rho\), \(\rho u\) and \(E\) with POD and the FCNN for the hydrodynamic regime (top row) and the rarefied regime (bottom row) at time \(t=0.12s\).
row). Indeed, fig. 8 indicates that conservation of mass and energy are achieved with

reasonable accuracy for the FCNN, while the error is larger for the POD reconstruction for both regimes. Also the increase in momentum of the full order model (FOM) is accurately described by the FCNN with a larger error for the POD method for both regimes. Overall, the errors are slightly smaller for the rarefied case compared to the hydrodynamic case, which might be due to the higher capacity of the neural network and more modes for the POD using five intrinsic variables in the rarefied regime.

### Physical interpretability

An important question in the context of model reduction with neural networks is the interpretability of the results, because of the usual black-box nature of neural networks [10]. Especially when benchmarking neural networks for model order reduction with POD, evaluating the interpretability of the intrinsic variables is important, since POD by construction achieves a so-called physically interpretable decomposition of the input data [3] as outlined previously.

Following the assumption that the hydrodynamic case can be fully described in terms of three macroscopic quantities and that the rarefied case is reasonably describable in a similar way with an extended set of five variables, we test the intrinsic variables \(\mathbf{h}\) and \(\mathbf{r}\) for similarities and investigate if they match any macroscopic quantity. Two related macroscopic quantities, namely the temperature \(T\) and macroscopic velocity \(u\), are added to the three macroscopic variables \(\rho\), \(\rho u\), and \(E\). In fig. 9 and fig. 12, these are plotted first over the whole domain of \(x\) and \(t\) and for the end time \(t=0.12s\) for both regimes. Similarly, we plot both the FCNN's and the POD's first 3 intrinsic variables \(\mathbf{h}\) of the hydrodynamic

Figure 8: Comparison of the conservative properties of reconstructions obtained from POD and the FCNNs with the conservative properties of the FOM solution using the mass, momentum, and energy derivative. A value of 0 indicates conservation. FCNN approximately conserves mass and energy, while the momentum increases with the correct rate of the test case.
case and 5 intrinsic variables \(\mathsf{r}\) of the rarefied test case

\[\mathsf{h} =[h_{0}(x,t),h_{1}(x,t),h_{2}(x,t)],\] \[\mathsf{r} =[r_{0}(x,t),r_{1}(x,t),r_{2}(x,t),r_{3}(x,t),r_{4}(x,t)],\]

depicted in fig. 10 and fig. 13 respectively.

Strikingly, most intrinsic variables of the FCNN in fig. 10 and fig. 13 and the POD in fig. 11 and fig. 14 appear to be a combination of the five intrinsic variables shown in fig. 9 and fig. 12, respectively. In particular consider FCNN in the hydrodynamic case by comparing fig. 9 and fig. 10. The rarefaction wave, shock wave, and contact discontinuity, which can be identified in \(h_{0}\) reflect a combination of those found in the density \(\rho\) and the total energy \(E\). Furthermore, \(h_{1}\) seems to model the negative momentum \(\rho u\), with different boundary values. The temperature \(T\) appears linked to \(h_{2}\), where the same fluctuation appears. Similar results hold true for the POD variables in fig. 11, where especially the first two

Figure 11: First three intrinsic variables \(h_{0}(x,t)\), \(h_{1}(x,t)\) and \(h_{2}(x,t)\) of hydrodynamic case obtained by the POD. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

Figure 10: Intrinsic variables \(h_{0}(x,t)\), \(h_{1}(x,t)\) and \(h_{2}(x,t)\) of hydrodynamic case obtained by the FCNN. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

Figure 9: Macroscopic quantities of the hydrodynamic case obtained by the FOM. Density \(\rho\), momentum \(\rho u\), total energy \(E\), temperature \(T\), and velocity \(u\) over time \(t\) and space \(x\) in the top row and at \(t=0.12\) in the bottom row.
Figure 14: First five intrinsic variables \(r_{0}(x,t)\), \(r_{1}(x,t)\), \(r_{2}(x,t)\), \(r_{3}(x,t)\) and \(r_{4}(x,t)\) of rarefied case obtained by the POD. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

Figure 12: Macroscopic quantities of the rarefied case obtained by the FOM. Density \(\rho\), momentum \(\rho u\), total energy \(E\), temperature \(T\), and velocity \(u\) over time \(t\) and space \(x\) in the top row and at \(t=0.12\) in bottom row.

Figure 13: Intrinsic variables \(r_{0}(x,t)\), \(r_{1}(x,t)\), \(r_{2}(x,t)\), \(r_{3}(x,t)\) and \(r_{4}(x,t)\) of rarefied case obtained by the FCNN. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).
variables closely resemble the density \(\rho\) and the momentum \(\rho u\). Interestingly, this does not relate to very good conservation properties of the POD in comparison with FCNN as shown in the previous section and fig. 8. Considering the FCNN in the rarefied case, we compare fig. 12 and fig. 13. Here \(r_{3}\) clearly reflects the shape of the density \(\rho\). Moreover, the peak of the velocity \(u\) can be observed in \(r_{0}\). For the other intrinsic variables of \(r_{1}\), \(r_{2}\) and \(r_{4}\) a clear discernability of macroscopic quantities is difficult to observe and might require linear or nonlinear combinations. Additionally, it is possible that those intrinsic variables resemble non-equilibrium variables not present among the macroscopic variables. Considering the POD results in fig. 14, we again see a very good agreement of the first intrinsic variables with density and momentum.

For more physical insight into the relation between macroscopic variables and intrinsic variables, the Pearson correlation between all variable combinations is computed in [12] for the FCNN. Note that we expect similar results for the POD based on the previous results, but do not present them here for conciseness. The Pearson correlation coefficient \(r_{X,Y}=r_{Y,X}\) is a measure of linear correlation between two sets of data, here represented by a macroscopic variable \(X\in\{\rho,\rho u,E,T,u\}\) and an intrinsic variable \(Y\in\{r_{0},r_{1},r_{2},\}\) for the hydrodynamic case and \(Y\in\{r_{0},r_{1},r_{2},r_{3},r_{4}\}\) for the rarefied case. It is commonly defined as

\[r_{X,Y}=\frac{\sum\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sum \sqrt{(x_{i}-\bar{x})}\sqrt{\sum\left(y_{i}-\bar{y}\right)}}. \tag{23}\]

Note that \(r_{X,Y}\in[-1,1]\), with \(r_{X,Y}=0\) meaning that there is no correlation between both data sets, \(r_{X,Y}=1\) indicating a perfect correlation, and \(r_{X,Y}=-1\) indicating a perfect anti-correlation.

The Pearson correlation coefficients for the hydrodynamic test case are presented in fig. 15. As predicted by the previous analysis, there appears to be an almost perfect correlation of the first intrinsic variable \(r_{0}\) with the density \(\rho\). This means that the FCNN succeeds at identifying the density \(\rho\) precisely as an internal variable. Note that \(r_{0}\) also correlates almost perfectly with the energy \(E\), as the energy depends linearly on \(\rho\). Additionally, there is a relatively strong linear correlation between \(r_{2}\) and both the momentum \(\rho u\) as well as the temperature \(T\). The intrinsic variable \(r_{2}\) on the other hand seems to be correlated with all variables.

The Pearson correlation coefficients for the rarefied test case are presented in fig. 16. In agreement with the previous results, there is no clear correlation of most of the intrinsic variables. An exception is \(r_{0}\), which is anti-correlated with \(u\) and \(r_{3}\), which correlates almost perfectly with the energy \(E\). Note that only linear correlations are tested here. For further analysis on the disentanglement of the intrinsic variables, it might be more suitable to consider nonlinear correlations beyond the Pearson correlation coefficients. One option that could be explored in future work would be to train the functional relation between the intrinsic variables and macroscopic variables, for example by means of symbolic regression [4].

## 5 Conclusion

This paper marks the first comparison of velocity space model reduction techniques for rarefied flows using proper orthogonal decomposition (POD) and a fully connected neural network (FCNN).

As physically expected, the rarefied regime needs more modes than the hydrodynamic regime. Choosing three and five intrinsic variables for the hydrodynamic and rarefied case, respectively, leads to less than one percent error. The FCNN is initially more accurate but has a remaining error even for larger number of intrinsic variables, while POD achieves subsequently higher accuracy with increasing the number of parameters. The resulting errors of the macroscopic variables are small, especially in the smoother rarefied case.

Even though not strictly enforced, the FCNN approximately exhibits the correct conservation of mass, momentum, and energy, while POD has a slightly larger error.
Figure 16: Pearson correlation between macroscopic quantities and intrinsic variables for the rarefied case.

Figure 15: Pearson correlation between macroscopic quantities and intrinsic variables for the hydrodynamic case.
The correlation of intrinsic variables and macroscopic variables is investigated by means of the evolution of the reconstructed values and the pairwise correlation for the FCNN. The density is directly included in the latent space while the relation with other macroscopic variables seems to be more complex.

In addition to optimizing the neural network's performance, our research endeavors to enhance its predictive capabilities by incorporating fundamental physical properties such as conservation and interpretability. This objective includes the exploration of different loss functions and the regularization techniques. Furthermore, future work can utilize the model reduction approach in real-world scenarios, including simulations and parameter predictions. Lastly, the presented concepts could be applied in related fields, such as shallow flows [16] or fusion plasmas [21].

## Author Contribution Statement (CRediT)

All authors contributed to this publication. The authors contributions differ in the following points:

**Julian Koellermeier:**: initial idea, methodology, setup of test cases, reference data, writing
**Philipp Krah:**: initial idea, methodology, implementation, visualization, writing
**Julius Reiss:**: initial idea, supervision
**Zachary Schellin:**: implementation, training, visualization

## Conflict of Interest

The authors declare that they have no conflict of interest.

## Code and Data Availability

All data and scripts to reproduce the results can be obtained from the authors on reasonable request.

## Acknowledgement

The authors would like to acknowledge the financial support of the CogniGron research center and the Ubbo Emmius Funds (University of Groningen). Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. P. Krah acknowledges partial funding from the Agence Nationale de la Recherche (ANR), project CM2E, grant ANR-20-CE46-0010-01.
