# Entropy-Stable Model Reduction of One-Dimensional Hyperbolic Systems using
Rational Quadratic Manifolds

R.B. Klein<EMAIL_ADDRESS>B. Sanderse<EMAIL_ADDRESS>P. Costa
<EMAIL_ADDRESS>R. Pecnik<EMAIL_ADDRESS>R.A.W.M. Henkes
<EMAIL_ADDRESS>Centrum Wiskunde & Informatica, Science Park 123,
Amsterdam, The Netherlands  Delft University of Technology, Process and
Energy, Leeghwaterstraat 39, Delft, The Netherlands  Eindhoven University of
Technology, Department of Mathematics and Computer Science, PO Box 513,
Eindhoven, The Netherlands

###### Abstract

In this work we propose a novel method to ensure important entropy
inequalities are satisfied semi-discretely when constructing reduced order
models (ROMs) on nonlinear reduced manifolds. We are in particular interested
in ROMs of systems of nonlinear hyperbolic conservation laws. The so-called
entropy stability property endows the semi-discrete ROMs with physically
admissible behaviour. The method generalizes earlier results on entropy-stable
ROMs constructed on linear spaces. The ROM works by evaluating the projected
system on a well-chosen approximation of the state that ensures entropy
stability. To ensure accuracy of the ROM after this approximation we locally
enrich the tangent space of the reduced manifold with important quantities.
Using numerical experiments on some well-known equations (the inviscid Burgers
equation, shallow water equations and compressible Euler equations) we show
the improved structure-preserving properties of our ROM compared to standard
approaches and that our approximations have minimal impact on the accuracy of
the ROM. We additionally generalize the recently proposed polynomial reduced
manifolds to rational polynomial manifolds and show that this leads to an
increase in accuracy for our experiments.

###### keywords:

Entropy stability , Manifold Galerkin method , Reduced order models , Rational
quadratic manifolds , Nonlinear conservation laws

###### Proof.

##  1 Introduction

Conservation laws are nearly universally present in any branch of physics and
engineering e.g. fluid dynamics, structural mechanics, plasma physics and
climate sciences; they express the conservation of some physical quantity of
interest. Often such conservation laws are described by hyperbolic equations
or systems thereof [34]. Physicists and engineers are increasingly reaching to
simulation tools for approximate solutions. With the increase in computational
power of recent decades, very large-scale problems have indeed been solved to
a satisfactory accuracy. Nonetheless, some applications of major engineering
interest like those of a multi-query (e.g. design optimization [108],
uncertainty quantification [35]) or real-time nature (e.g. model predictive
control [90], digital twin technology [55]) remain out of question for many
large scale systems. A way around these computational issues has been offered
by reduced order models (ROMs), which are low-dimensional surrogates of high-
fidelity models of interest, often referred to as full order models (FOM) in
the ROM community. ROMs rely on the offline-online decomposition paradigm [76]
for their efficiency, they are trained in an expensive offline phase and
subsequently evaluated at very low computational cost in new situations in the
online phase. Particularly, the low dimensionality of the ROM allows for fast
and cheap evaluation of its solution.

A popular class of ROMs are projection-based ROMs (pROMs) [11, 108], which
have traditionally been constructed by projecting equations of interest on
well-chosen linear subspaces. These subspaces are often found in a data-driven
manner using the proper orthogonal decomposition (POD) [61, 113] or greedy
methods [100] and the projections are carried out using a Galerkin [109, 73,
107, 69, 68, 10] or Petrov-Galerkin [24, 51, 23] approach. However, the
success of applications of linear subspace-based pROMs has been limited in the
field of hyperbolic equations. This is a result of an almost inherently slow
decay of the so-called Kolmogorov $n$-width (KnW) of the solution manifolds
$\mathcal{M}_{u}$ (i.e. the set of all solution trajectories for a range of
parameters and initial conditions of interest) of these
systems111$(\mathcal{X},||\cdot||_{\mathcal{X}})$ denoting a Banach space of
interest containing analytical PDE or approximate simulation solutions:

$d_{r}(\mathcal{M}_{u})=\inf_{\mathcal{V}\subset\mathcal{X};\text{dim}(\mathcal{V})=r}\sup_{\boldsymbol{u}\in\mathcal{M}_{u}}\inf_{\boldsymbol{v}\in\mathcal{V}}||\boldsymbol{u}-\boldsymbol{v}||_{\mathcal{X}},$
(1)

which measures the worst error that can be incurred when optimally
approximating $\mathcal{M}_{u}$ with an optimally chosen $r$-dimensional
linear subspace. The KnW decay of many hyperbolic systems is slow because
their solution trajectories are often not contained in low-dimensional linear
subspaces due to characteristically having moving features as part of their
solution. This has been shown analytically for some systems [99, 32, 49, 86,
6] and empirically for many others. In recent years a range of possible
solutions have been proposed, falling in roughly four categories. First, there
have been adaptive approaches that use linear subspaces that are changed
during the online phase to be better suited to new conditions [98, 96, 93, 56,
22, 39, 124, 45]; second, domain decomposition approaches that localize ROM
construction in time, space or parameter space [123, 30, 97, 63, 14, 1, 51,
33]; third, Lagrangian, registration and/or optimal transport based approaches
that track moving features improving linear data compressibility [91, 48, 64,
103, 88, 116, 38, 104, 89, 71, 122, 13, 87]; fourth, constructing ROMs on
nonlinear spaces (manifolds). Since it is not dependent on user-defined
localization or adaptation strategies and due to its high expressiveness we
will be interested in the latter category. In model reduction on manifolds the
linear reduced spaces of classical ROMs are replaced by nonlinear spaces.
Given sufficient expressiveness the nonlinear reduced spaces can potentially
approximate the solution manifolds of hyperbolic systems. These manifolds are
typically constructed from data. In [105, 106, 79, 18, 121, 31, 72, 37]
autoencoder neural networks are used to construct manifold ROMs. Another
popular nonlinear manifold construction approach are polynomial manifolds [47,
17, 111, 8, 66]. Finally, there has been an increasing interest in a more
physics-based approach, where ROMs are constructed on the invariant manifolds
of physical systems [5, 20, 92, 21].

Although, to the authors’ knowledge, there have not been any studies of
manifold ROMs on large scale under-resolved and shock-dominated cases within
fluid dynamical applications, it has long been known that their linear
counterparts can suffer from stability issues for such problems [109, 73, 41,
7, 2, 69, 68, 10]. It is quite reasonable to assume this will also be the case
for manifold ROMs. A promising approach to stabilization of such simulations
is the concept of entropy stability [118, 82], which has been widely used in
the context of obtaining stable FOMs. A numerical method is entropy stable if
it dissipates a convex functional associated to a conservation law, referred
to as entropy, given suitable boundary conditions. Entropy stability endows
numerical methods with physically admissible behaviour, which for fluid
dynamical applications manifests itself in satisfaction of the second law of
thermodynamics. Furthermore, it generalizes the $L^{2}$-stability properties
of linear hyperbolic systems to fully nonlinear systems [118]. Additionally,
stability in $L^{p}$ spaces can be shown formally [115, 44, 74].

However, in the projection step to construct pROMs an entropy-stable numerical
method generally loses its stability property. In recent years, some work has
been carried out in the ROM community to preserve the entropy stability
property. In Figure 1 several possible approaches have been visualized. In
[27] the ROM is evaluated at a corrected state that ensures entropy stability
is maintained (also known as ‘entropy projection’). Inspired by the approach
taken in the classical finite element work in [62], [69] writes the
conservation laws in symmetric form using an alternative set of variables and
projects the continuous equations leading to correct entropy estimates. In
[94] the Hilbert spaces in which the projections are carried out are defined
according to physical arguments. We also note the works in [68, 10, 107]. All
these approaches suffer a major drawback – they are built on linear reduced
spaces. As a result, relatively high-dimensional reduced spaces are required
to model many physical systems of interest, which comes at the cost of
computational efficiency of the ROM.

Our main contribution is to generalize these entropy-stable approaches to
nonlinear reduced spaces. This allows for lower-dimensional reduced spaces and
thus potentially more efficient ROMs. In particular, we will be interested in
generalizing the work in [27]. This approach is colored red in Figure 1. Our
main argument for not choosing [69, 94] stems from the argument given e.g. in
[25], namely that formulations in alternative variables are not consistent
with the famous Lax-Wendroff theorem [77] and can thus yield wrong shock
solutions. To our knowledge, our method is the first manifold ROM that is
provably entropy stable. At the same time, we note that preservation of other
mathematical structures on reduced manifolds has been successfully achieved in
the past: symplectic [111, 18, 16]; metriplectic [80]; conservative [78]. A
great overview of recent structure-preserving model reduction contributions in
general is given in [52].

A second contribution of this work is the development of a novel
generalization of polynomial manifolds [8, 17] to rational polynomial
manifolds. While these polynomial manifolds have shown successes in certain
applications, in our experience they are not sufficiently accurate for shock-
dominated problems. Rational polynomials are better at capturing
discontinuities, as we will show in this work.

Figure 1: Flowchart of possible entropy stable ROM approaches, our method is
indicated in red, $u$ are the conserved variables and $\eta$ are the alternate
variables used in [69]. FOM: full order model, ROM: reduced order model.

This article is organized as follows. In section 2 we introduce the theory of
nonlinear hyperbolic conservation laws and entropy analysis, we introduce a
baseline entropy stable FOM [44], and an entropy-stable linear ROM as proposed
in [27]. In section 3 we introduce our main contribution, the novel entropy
stable nonlinear manifold ROM. In section 4 we discuss our second
contribution, being rational polynomial manifolds. In section 5 we show the
effectiveness of our approach using several numerical experiments that are
based on a range of well-known conservation laws from fluid dynamics. We
conclude our work in section 6.

##  2 Preliminaries: Entropy inequality for conservation laws, entropy-stable
FOM, linear ROM

###  2.1 Introducing the entropy inequality

We give a short introduction to the concept of entropy, some related concepts
used in its analysis and its role in the theory of nonlinear conservation
laws. We consider conservation laws in one spatial dimension that can be
written as partial differential equations (PDE) of the form:

$\frac{\partial\boldsymbol{u}}{\partial
t}+\frac{\partial\boldsymbol{f}(\boldsymbol{u})}{\partial x}=0,$ (2)

where $\Omega$ is a spatial domain and $[0,T]$ is a temporal domain with
$T>0$. Furthermore $\boldsymbol{u}:\Omega\times(0,T]\rightarrow\mathbb{R}^{n}$
is the solution function,
$\boldsymbol{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the nonlinear flux
function, $n\in\mathbb{N}$ is the number of conserved quantities and
$x\in\Omega$ and $t\in[0,T]$ are the spatial coordinate and time,
respectively. To facilitate conservation statements and minimize the role of
boundary conditions we will focus on periodic spatial domains
$\Omega=\mathbb{T}([a,b])$ with $b>a$ and $b,a\in\mathbb{R}$ in this research,
here $\mathbb{T}$ is the torus. The equations are complemented by a set of
initial conditions $\boldsymbol{u}_{0}:\Omega\rightarrow\mathbb{R}^{n}$ so
that $\boldsymbol{u}_{0}(x)=\boldsymbol{u}(x,0)$. The conservation law (2) is
referred to as hyperbolic when the Jacobian matrix
$\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{u}}$ is diagonalizable with
real eigenvalues for all physically relevant $\boldsymbol{u}$ [82]. In many
cases, the solutions of physically relevant hyperbolic conservation laws also
satisfy additional conservation laws of the form:

$\frac{\partial s(\boldsymbol{u})}{\partial
t}+\frac{\partial\mathcal{F}(\boldsymbol{u})}{\partial x}=0,$ (3)

where the function $s:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is called the
entropy function which is defined to be convex222A function
$g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is convex if its Hessian,
$\frac{\partial^{2}g}{\partial\boldsymbol{u}^{2}}(\boldsymbol{u})$, is
positive definite for all $\boldsymbol{u}$. and
$\mathcal{F}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is called the entropy flux.
In particular, such an additional conservation law exists if the compatibility
relation:

$\boldsymbol{\eta}(\boldsymbol{u})^{T}\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{u}}=\frac{\partial\mathcal{F}}{\partial\boldsymbol{u}}^{T},$
(4)

is satisfied. Here,
$\boldsymbol{\eta}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n},\hbox{
}\boldsymbol{u}\mapsto\frac{\partial
s}{\partial\boldsymbol{u}}(\boldsymbol{u})$ is the gradient of the entropy
function $s$ with respect to $\boldsymbol{u}$. This mapping is injective due
to the convexity of $s$ and hence can be inverted. A pair $(s,\mathcal{F})$
satisfying the compatibility relation (4) is called an entropy pair of (2).

It is well-known that solutions $\boldsymbol{u}$ to (2) can develop
discontinuities in finite time for smooth $\boldsymbol{u}_{0}$ [40, 82, 81].
When this occurs the solution $\boldsymbol{u}$ is said to contain a shock or
contact discontinuity depending on the behaviour of the discontinuity [82]. In
this case both formulation (2) and the manipulations to obtain (3) are no
longer valid. We must therefore consider (2) in a weak sense to retain a
notion of solutions. This weak form of the conservation law is obtained by
integrating against a space of smooth test functions
$\boldsymbol{v}:\Omega\times[0,\infty)\rightarrow\mathbb{R}^{n}$ with compact
support i.e. $\boldsymbol{v}\in C_{0}^{\infty}(\Omega\times[0,\infty))$ and
transferring all derivatives to these test functions to obtain:

$\int_{0}^{\infty}\int_{\Omega}\boldsymbol{u}\cdot\frac{\partial\boldsymbol{v}}{\partial
t}+f(\boldsymbol{u})\cdot\frac{\partial\boldsymbol{v}}{\partial
x}dxdt+\int_{\Omega}\boldsymbol{u}_{0}\cdot\boldsymbol{v}(x,0)dx=0,$ (5)

on periodic $\Omega$. Note that this expression is valid even for
discontinuous $\boldsymbol{u}$ and that any smooth $\boldsymbol{u}$ satisfying
the strong form (2) also satisfies this weak form (5) [34, 40]. However,
weakening the notion of a solution like this comes at the cost that it does
not necessarily yield unique solutions (examples of such cases may be found in
[40, 81, 34, 82]). Out of all weak solutions, i.e. solutions to (5), those of
physical interest are the ones satisfying:

$\frac{\partial s(\boldsymbol{u})}{\partial
t}+\frac{\partial\mathcal{F}(\boldsymbol{u})}{\partial x}\leq 0,$ (6)

in the sense of distributions, with equality for smooth $\boldsymbol{u}$
following (3) and inequality for solutions containing shocks. This inequality
arises from considering limits of regularized conservation laws, more details
can be found in [82]. For scalar conservation laws Kruzkov [74] established
that weak solutions satisfying (6) are unique, but for systems uniqueness is
not yet completely established [12]. We can define the total entropy
functional as:

$\mathcal{S}[\boldsymbol{u}]:=\int_{\Omega}s(\boldsymbol{u})dx.$ (7)

Defining appropriate sequences of test functions and taking limits [83, 34],
it can be shown that the estimate:

$\frac{d\mathcal{S}[\boldsymbol{u}]}{dt}\leq 0,$ (8)

follows from (6) on periodic $\Omega$. This inequality will be the main
interest of this paper.

###  2.2 Entropy stable spatial discretization (FOM)

We will discretize the conservation law (2) with a finite volume method (FVM)
based on flux-differencing [82]. To introduce the general and frequently
recurring structure of our entropy stability proof we will provide some detail
on the full order model (FOM) discretization. The FVM discretization will be
constructed such that discrete analogues to (8) hold. Other discretization
methods that similarly mimic (8) are also possible, for example the split-form
discontinuous Galerkin (DG) methods described in [46], the summation-by-parts
schemes in [101] and the higher-order methods of [44, 42]. We choose the FVM
to keep the exposition simple, but note that our entropy-stable ROM framework
should work with other entropy-stable FOM discretizations as well.

The scheme is formulated as:

$\Delta
x_{i}\frac{d\boldsymbol{u}_{i}}{dt}+\boldsymbol{f}_{i+1/2}-\boldsymbol{f}_{i-1/2}=0,\quad
i\in\\{0,...,N-1\\},$ (9)

on a grid of $N$ grid cells so that $i\in\\{0,...,N-1\\}$. Here,
$\boldsymbol{u}_{i}:[0,\infty)\rightarrow\mathbb{R}^{n}$ is the numerical
solution vector in the $i$-th grid cell, $\Delta x_{i}:=x_{i+1/2}-x_{i-1/2}$
is the cell size of the $i$-th cell with $x_{i\pm 1/2}$ denoting respectively
the $x$ values of the right ($+$) and left ($-$) cell boundary and
$\boldsymbol{f}_{i+1/2}:=\boldsymbol{f}_{h}(\boldsymbol{u}_{i+1},\boldsymbol{u}_{i})$
is the numerical flux on the right cell boundary of the $i$-th cell with
$\boldsymbol{f}_{h}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$
being a two-point numerical flux function [117] approximating the flux
function $\boldsymbol{f}$ on a cell boundary based on two neighbouring
numerical solution values, similarly $\boldsymbol{f}_{i-1/2}$ approximates
$\boldsymbol{f}$ on the left boundary. We also define the total number of
unknowns as $N_{h}:=n\cdot N$. Periodic boundary conditions are enforced by
setting $\boldsymbol{u}_{N}:=\boldsymbol{u}_{0}$ and
$\boldsymbol{u}_{-1}:=\boldsymbol{u}_{N-1}.$ In the schemes we are considering
the numerical flux function is constructed from entropy-conservative flux
functions [118, 117, 119]. These are flux functions that assure discrete
analogues of (3) and (8) are satisfied with equality. This makes them suitable
starting points from which to construct flux functions that have appropriate
entropy-dissipative properties. We follow Tadmor’s framework [117] of entropy-
conservative fluxes, which are defined as follows:

######  Definition 1 (Entropy-conservative numerical flux).

An entropy-conservative two-point numerical flux
$\boldsymbol{f}^{*}_{h}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$
is a numerical two-point flux satisfying:

  1. 1.

consistency:
$\boldsymbol{f}^{*}_{h}(\boldsymbol{u},\boldsymbol{u})=\boldsymbol{f}(\boldsymbol{u})$;

  2. 2.

symmetry:
$\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{l},\boldsymbol{u}_{r})=\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{r},\boldsymbol{u}_{l})$;

  3. 3.

entropy conservation:
$(\boldsymbol{\eta}(\boldsymbol{u}_{l})-\boldsymbol{\eta}(\boldsymbol{u}_{r}))^{T}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{l},\boldsymbol{u}_{r})=\psi(\boldsymbol{u}_{l})-\psi(\boldsymbol{u}_{r})$.

The entropy-dissipative fluxes
$\boldsymbol{f}_{h}:\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{N_{h}}\rightarrow\mathbb{R}^{n}$
are now constructed from entropy-conservative fluxes $\boldsymbol{f}_{h}^{*}$
by adding (possibly solution dependent) entropy dissipation operators like:

$\boldsymbol{f}_{i+1/2}:=\boldsymbol{f}_{h}(\boldsymbol{u}_{i+1},\boldsymbol{u}_{i},\boldsymbol{u}_{h})=\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{i+1},\boldsymbol{u}_{i})-\boldsymbol{D}_{i+1/2}(\boldsymbol{u}_{h})\Delta\boldsymbol{\eta}_{i+1/2}$
(10)

with $\boldsymbol{D}_{i+1/2}:\mathbb{R}^{N_{h}}\rightarrow\mathbb{S}^{n}_{+}$
so that $\boldsymbol{D}_{i+1/2}(\boldsymbol{u}_{h})$ is symmetric positive
semi-definite (SPSD) for any $\boldsymbol{u}_{h}(t)\in\mathbb{R}^{N_{h}}$
($\mathbb{S}^{n}_{+}$ and $\mathbb{S}^{n}_{++}$ are the convex sets of
symmetric positive definite and symmetric positive semi-definite $n\times n$
matrices, respectively). Here,
$\boldsymbol{u}_{h}:[0,\infty)\rightarrow\mathbb{R}^{N_{h}}$ is the numerical
solution vector on the whole grid to be defined in what follows (this is
required for higher order reconstructions like in e.g. [44]). Additionally, we
have defined
$\Delta\boldsymbol{\eta}_{i+1/2}:=\boldsymbol{\eta}(\boldsymbol{u}_{i+1})-\boldsymbol{\eta}(\boldsymbol{u}_{i})$.
We will refrain from denoting explicitly the dependence on
$\boldsymbol{u}_{h}$ in the third argument of $\boldsymbol{f}_{h}$ and simply
write $\boldsymbol{f}_{h}(\boldsymbol{u}_{i+1},\boldsymbol{u}_{i})$ for (10).

For the purpose of model reduction in section 3 we rewrite discretization (9)
with flux (10) in a matrix-vector formulation. We will introduce the following
notations: volume-based quantities which live on cell centers and interface-
based quantities which live on cell interfaces. The numerical solution vector
is a volume-based quantity given by:

$\boldsymbol{u}_{h}(t):=[u^{1}_{0}(t)...,u^{k}_{i}(t),u^{k}_{i+1}(t),...,u^{n}_{i}(t),...]^{T}\in\mathbb{R}^{N_{h}},$

where $u^{k}_{i}(t)\in\mathbb{R}$ is the approximation of the $k$-th conserved
variable (the variable conserved by the $k$-th equation in (2)) in the $i$-th
cell evaluated at time $t$. The numerical flux vector is an interface-based
quantity, overloading the notation for the numerical flux function, it is
given by:

$\boldsymbol{f}_{h}(\boldsymbol{u}_{h}):=[f^{1}_{1/2},...,f^{1}_{N-1/2},f^{2}_{1/2},...f^{k}_{i-1/2},f^{k}_{i+1/2},...,f^{n}_{i+1/2},...,f^{n}_{N-1/2}]^{T}\in\mathbb{R}^{N_{h}},$

where $f^{k}_{i+1/2}$ is the numerical flux of the $k$-th conservation
equation evaluated at interface $i+1/2$ between cells $i+1$ and $i$. Periodic
boundary conditions are built into the flux vector by evaluating
$f^{k}_{N-1/2}(\boldsymbol{u}_{0},\boldsymbol{u}_{N-1})$ for all
$k\in\\{1,...,n\\}$.

To perform finite-difference operations as in (9) for volume-based and
interface-based quantities, respectively, the following matrices are defined:

$\bar{\Delta}_{v}=\begin{bmatrix}1&0&0&-1\\\ -1&1&0&0\\\ 0&\ddots&\ddots&0\\\
0&0&-1&1\\\ \end{bmatrix}\in\mathbb{R}^{N\times
N},\quad\bar{\Delta}_{i}=\begin{bmatrix}-1&1&0&0\\\ 0&\ddots&\ddots&0\\\
0&0&-1&1\\\ 1&0&0&-1\\\ \end{bmatrix}\in\mathbb{R}^{N\times N}.$

We note the skew-adjointness relation

$\bar{\Delta}_{v}=-\bar{\Delta}_{i}^{T},$ (11)

and that both have zero row and column sum. These properties will be used in
proving entropy stability of the scheme. For systems we define
$\Delta_{v}:=I\otimes\bar{\Delta}_{v}$ and
$\Delta_{i}:=I\otimes\bar{\Delta}_{i}$ with $I$ being the $n\times n$ identity
matrix and $\otimes$ the Kronecker product. Clearly, $\Delta_{v}$ and
$\Delta_{i}$ satisfy a similar skew-adjointness property (11). We will also
introduce the FVM mass matrices $\bar{\Omega}_{h}=\text{diag}(\Delta x_{i})$
with $i=0,1,...,N-1$ and $\Omega_{h}:=I\otimes\bar{\Omega}_{h}$. With these
operators we can write a compact form of the discretization (9) as follows:

$\Omega_{h}\frac{d\boldsymbol{u}_{h}}{dt}+\Delta_{v}\boldsymbol{f}_{h}(\boldsymbol{u}_{h})=0.$
(12)

To emphasize the role played by the dissipation operator
$\boldsymbol{D}_{i+1/2}$ in obtaining entropy-stable spatial discretizations
we will decompose $\Delta_{v}\boldsymbol{f}_{h}(\boldsymbol{u}_{h})$ in an
entropy-conserving part and an entropy-dissipating part, resulting in:

$\Omega_{h}\frac{d\boldsymbol{u}_{h}}{dt}+\Delta_{v}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{h})=\Delta_{v}\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\Delta_{i}\boldsymbol{\eta}_{h},$
(13)

here, $\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{h})$ is a vector of entropy
conservative numerical fluxes,
$\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\in\mathbb{S}_{+}^{N_{h}}$ is an SPSD
matrix containing the terms associated to the dissipation operators
$\boldsymbol{D}_{i+1/2}$ and $\boldsymbol{\eta}_{h}$ is a vector containing
the grid values of the entropy variables ordered similarly as
$\boldsymbol{u}_{h}$.

Having introduced an entropy-dissipative numerical flux, we evaluate the
discrete analogue to the continuous total entropy functional (7) which should
be suitably dissipated by the entropy stable discretization (13) or conserved
in the case of no dissipation. The discrete total entropy functional will be
defined as:

$S_{h}[\boldsymbol{u}_{h}]:=\boldsymbol{1}^{T}\bar{\Omega}_{h}\boldsymbol{s}_{h},$
(14)

where $\boldsymbol{1}$ is a vector of ones and the local entropy is defined
as:

$\boldsymbol{s}_{h}(\boldsymbol{u}_{h}(t)):=[s(\boldsymbol{u}_{0}(t)),...,s(\boldsymbol{u}_{i}(t)),...,s(\boldsymbol{u}_{N-1}(t))]^{T}\in\mathbb{R}^{N},$

which is a volume-based quantity. The time evolution of $S_{h}$ is given by

$\frac{dS_{h}[\boldsymbol{u}_{h}]}{dt}=\boldsymbol{1}^{T}\bar{\Omega}_{h}\frac{d\boldsymbol{s}_{h}}{dt}=\sum_{i}\Delta
x_{i}\boldsymbol{\eta}(\boldsymbol{u}_{i})^{T}\frac{d\boldsymbol{u}_{i}}{dt}=\boldsymbol{\eta}_{h}^{T}\Omega_{h}\frac{d\boldsymbol{u}_{h}}{dt}.$

We also define the entropy flux potential vector as:

$\boldsymbol{\psi}_{h}(\boldsymbol{u}_{h}(t)):=[\psi(\boldsymbol{u}_{0}(t)),...,\psi(\boldsymbol{u}_{i}(t)),...,\psi(\boldsymbol{u}_{N-1}(t))]^{T}\in\mathbb{R}^{N},$

which is a volume-based quantity, like the local entropy vector. To analyse
the entropy evolution we can substitute the spatial discretization (13) in the
previous expression to obtain:

$\displaystyle\frac{dS_{h}[\boldsymbol{u}_{h}]}{dt}$
$\displaystyle=\boldsymbol{\eta}_{h}^{T}\Omega_{h}\frac{d\boldsymbol{u}_{h}}{dt}$
$\displaystyle=-\boldsymbol{\eta}_{h}^{T}\Delta_{v}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{h})+\boldsymbol{\eta}_{h}^{T}\Delta_{v}\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\Delta_{i}\boldsymbol{\eta}_{h}$
$\displaystyle=(\Delta_{i}\boldsymbol{\eta}_{h})^{T}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}_{h})-\boldsymbol{\eta}_{h}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\Delta_{i}\boldsymbol{\eta}_{h}$
$\displaystyle=\boldsymbol{1}^{T}\bar{\Delta}_{i}\boldsymbol{\psi}_{h}-\boldsymbol{\eta}_{h}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\Delta_{i}\boldsymbol{\eta}_{h}$
$\displaystyle=0-\boldsymbol{\eta}_{h}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\Delta_{i}\boldsymbol{\eta}_{h}$
$\displaystyle\leq 0,$ (15)

where we used the skew-adjointness property (11), the entropy conservation
condition of the numerical flux, positive-definiteness of the dissipation
operator $\boldsymbol{D}_{h}(\boldsymbol{u}_{h})$ (and thus of
$\Delta_{i}^{T}\boldsymbol{D}_{h}(\boldsymbol{u}_{h})\Delta_{i}$) and the zero
column sum of $\bar{\Delta}_{i}$. Clearly, in case no entropy dissipation is
added in the numerical flux (10), equation (15) reduces to

$\frac{dS_{h}[\boldsymbol{u}_{h}]}{dt}=0.$

We note that the inequality (15) allows for formal $L^{p}$-stability
statements [115, 34, 44].

###  2.3 The entropy-stable linear ROM of Chan [27]

The main aim of this work is to propose reduced order models (ROMs) that are a
nonlinear generalization of the entropy stable ROM of Chan [27]. To highlight
key conceptual differences between the ROM in [27] and ours, and to introduce
the ROM methodology, we will briefly discuss the elements of Chan’s ROM
leading to its entropy stability. Classical reduced order models including
[27] make the assumption that the evolution of $\boldsymbol{u}_{h}(t)$ can be
accurately approximated with elements from a linear space
$\mathcal{V}\subset\mathbb{R}^{N_{h}}$ where $\text{dim}(\mathcal{V}):=r\ll
N_{h}$ so that $\mathcal{V}$ can be referred to as low-dimensional [73, 107,
69, 68, 10, 24, 51, 23]. Classically, the subspace $\mathcal{V}$ is
constructed using a truncated proper orthogonal decomposition (POD) based on
snapshot data collected in a matrix $X\in\mathbb{R}^{N_{h}\times n_{s}}$ with
$n_{s}\in\mathbb{N}$ the number of snapshots [109, 73]. The construction of
the ROM in [27] starts by defining the approximation
$\boldsymbol{u}_{h}\approx\boldsymbol{u}_{r}:=\Phi\boldsymbol{a}$ with
$\boldsymbol{a}\in\mathbb{R}^{r}$ being generalized coordinates in
$\mathcal{V}$ relative to the basis $\Phi$. Here, we assume $\Phi$ is
orthogonal in the $\Omega_{h}$-weighted inner product, i.e.
$\phi_{i}^{T}\Omega_{h}\phi_{j}=\delta_{ij}$ with $\phi_{i}$ the $i$-th column
of $\Phi$ and $\delta_{ij}$ the Kronecker delta function. Then, the
approximation is substituted in (13) introducing a semi-discrete residual
which is set orthogonal to $\mathcal{V}$ by solving the Galerkin projected
system:

$\frac{d\boldsymbol{a}}{dt}+\Phi^{T}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})=\Phi^{T}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}(\boldsymbol{a}),$
(16)

with
$\boldsymbol{f}_{r}^{*}(\boldsymbol{a}):=\boldsymbol{f}_{h}^{*}(\Phi\boldsymbol{a})$,
$\boldsymbol{D}_{r}(\boldsymbol{a}):=\boldsymbol{D}_{h}(\Phi\boldsymbol{a})$
and
$\boldsymbol{\eta}_{r}(\boldsymbol{a}):=\boldsymbol{\eta}_{h}(\Phi\boldsymbol{a})$.
Equation (16) defines a (linear) POD-Galerkin ROM [73, 107, 69, 68, 10].

Similar to the FOM case, see equation (14), we can evaluate the evolution of
the ROM total entropy. The reduced total entropy functional is defined as:

$S_{r}[\boldsymbol{a}]:=S_{h}[\Phi\boldsymbol{a}]=\boldsymbol{1}^{T}\bar{\Omega}_{h}\boldsymbol{s}_{r}(\boldsymbol{a}),$
(17)

with
$\boldsymbol{s}_{r}(\boldsymbol{a}):=\boldsymbol{s}_{h}(\Phi\boldsymbol{a})$.
The evolution of the reduced total entropy is:

$\frac{dS_{r}[\boldsymbol{a}]}{dt}=\boldsymbol{1}^{T}\bar{\Omega}_{h}\frac{d\boldsymbol{s}_{r}}{dt}=\sum_{i}\Delta
x_{i}\boldsymbol{\eta}(\Phi_{i}\boldsymbol{a})^{T}\Phi_{i}\frac{d\boldsymbol{a}}{dt}=\boldsymbol{\eta}_{r}^{T}\Omega_{h}\Phi\frac{d\boldsymbol{a}}{dt},$

with $\Phi_{i}\in\mathbb{R}^{n\times r}$ the rows of $\Phi$ approximating
values in cell $i$. The entropy evolution of (16) is:

$\displaystyle\frac{dS_{r}[\boldsymbol{a}]}{dt}$
$\displaystyle=\boldsymbol{\eta}_{r}^{T}\Omega_{h}\Phi\frac{d\boldsymbol{a}}{dt}$
$\displaystyle=-\boldsymbol{\eta}_{r}^{T}\Omega_{h}\Phi\Phi^{T}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})+\boldsymbol{\eta}_{r}^{T}\Omega_{h}\Phi\Phi^{T}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}$
$\displaystyle=-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})+\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r},$

where, since $\Phi\Phi^{T}\Omega_{h}$ defines an $\Omega_{h}$-orthogonal
projection operator [4],
$\tilde{\boldsymbol{\eta}}_{r}:=\Phi\Phi^{T}\Omega_{h}\boldsymbol{\eta}_{r}$
are the so-called projected entropy variables. It is unclear whether this
expression is bounded. To solve this, Chan [27] proposes a technique used
earlier in DG finite element literature [95, 26, 28]. Namely, the
discretization is not evaluated at $\Phi\boldsymbol{a}$ but at the entropy
projected state:

$\tilde{\boldsymbol{u}}_{r}:=\boldsymbol{u}(\Phi\Phi^{T}\Omega_{h}\boldsymbol{\eta}_{r})=\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}),$

where we have defined
$\boldsymbol{u}:\mathbb{R}^{N_{h}}\rightarrow\mathbb{R}^{N_{h}},\boldsymbol{\eta}_{h}\mapsto\boldsymbol{u}_{h}$
for notational convenience. Recall that the mapping $\boldsymbol{u}$ is indeed
available since $\boldsymbol{\eta}$ is injective. In this case we have:

$\displaystyle\frac{dS_{r}[\boldsymbol{a}]}{dt}$
$\displaystyle=-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}))+\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle=(\Delta_{i}\tilde{\boldsymbol{\eta}}_{r})^{T}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}))-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle=\boldsymbol{1}^{T}\bar{\Delta}_{i}\tilde{\boldsymbol{\psi}}_{r}-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle=0-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle\leq 0,$ (18)

so that we re-obtain an entropy estimate that mimics the FOM estimate (15).
Here, $\tilde{\boldsymbol{\psi}}_{r}$ is the entropy flux potential vector
evaluated at the entropy projected state. A key difference with the DG
literature [95, 26, 28] is that in the ROM case the basis $\Phi$ is only
constructed to resolve solutions present in the snapshot matrix $X$, whereas
in DG the trial basis is able to approximate a larger subspace of the relevant
PDE function spaces. As a result, the DG trial basis can resolve the entropy
variables well, but this may not be the case for the reduced basis $\Phi$. To
address this issue, [27] builds the basis $\Phi$ from a set of augmented
snapshots given (with some abuse of notation) by:

$\tilde{X}=[X,\boldsymbol{\eta}(X)],$

so that the projection of the entropy variables on the basis $\Phi$ is close
to the identity. As we will explain in section 3.3, for our proposed nonlinear
manifold ROMs, such a construction is not sufficient, and a new tangent space
enrichment technique will be proposed to ensure the accuracy of the entropy
projection.

##  3 An entropy-stable nonlinear manifold Galerkin ROM

The solution manifolds of many hyperbolic conservation laws (2) have slow
Kolmogorov $n$-width decay (1). Hence, approximations using linear subspaces
as in [27] may require very large reduced space dimensions $r$ before they
become accurate. This comes at the cost of their efficiency. For this reason
ROMs built on nonlinear spaces endowed (at least locally) with a manifold
structure have become a topic of interest [105, 106, 79, 111, 18, 19, 9, 8,
121, 31, 72, 37, 47, 66]. To address the shortcomings of the linear subspaces
employed in [27] we will generalize this method to nonlinear reduced spaces,
while keeping the entropy-stability property. We will be interested
specifically in nonlinear subsets of $\mathbb{R}^{N_{h}}$ endowed with some
inner product, instead of any abstract space. Therefore we will not be very
rigorous about our use of the word manifold, following predominantly the
treatise of [79] and standard multivariable calculus interpretations. For a
rigorous treatment we suggest consulting the recent preprint [19]. We will
give a brief description of nonlinear manifold ROMs and then propose our
generalization of [27].

###  3.1 Manifold Galerkin model reduction

In constructing ROMs on nonlinear manifolds we make the assumption that for
any $t\in[0,T]$ there are points $\boldsymbol{u}_{r}(t)$ on a low-dimensional
submanifold $\mathcal{M}\subset\mathbb{R}^{N_{h}}$ that accurately approximate
$\boldsymbol{u}_{h}(t)$. Here, we denote $r:=\text{dim}(\mathcal{M})$ and the
low-dimensionality of $\mathcal{M}$ implies that $r\ll N_{h}$. We will refer
to the submanifold $\mathcal{M}$ as the reduced manifold. Instead of the
classical affine reduced space parameterization seen in the previous section
we will use nonlinear manifold parameterizations given as:

$\boldsymbol{u}_{h}(t)\approx\boldsymbol{u}_{r}(t):=\boldsymbol{\varphi}(\boldsymbol{a}(t))\in\mathcal{M},$
(19)

where:

$\boldsymbol{\varphi}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{N_{h}},$ (20)

is assumed to be a smooth nonlinear injective function - at least when
restricted to some subset $\mathcal{A}\subseteq\mathbb{R}^{r}$ of interest
where the ROM will be well-defined. This means that
$\boldsymbol{\varphi}(\mathbb{R}^{r})=\mathcal{M}$ with a Jacobian
$\boldsymbol{J}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{N_{h}\times r},\hbox{
}\boldsymbol{a}\mapsto\frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}}(\boldsymbol{a})$
of full rank for any $\boldsymbol{a}\in\mathcal{A}\subseteq\mathbb{R}^{r}$,
where $\boldsymbol{a}:[0,T]\rightarrow\mathbb{R}^{r}$ are generalized
coordinates on the manifold $\mathcal{M}$. The function $\boldsymbol{\varphi}$
may be obtained in many ways: some examples are quadratic approximations [8,
111, 17, 47, 66] or neural networks [105, 106, 79, 72, 121, 37, 9, 31]. We
will propose a new method based on rational polynomials in section 4. In this
section, we develop an entropy-stable ROM which is agnostic of the choice for
$\boldsymbol{\varphi}$.

To construct a ROM we substitute (19) into the FOM discretization (13) so that
after applying the chain rule we find the residual:

$\boldsymbol{r}\left(\frac{d\boldsymbol{a}}{dt},\boldsymbol{a}\right):=\boldsymbol{J}\frac{d\boldsymbol{a}}{dt}+\Omega_{h}^{-1}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})-\Omega_{h}^{-1}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}(\boldsymbol{a}),$

where we define
$\boldsymbol{f}_{r}^{*}(\boldsymbol{a}):=\boldsymbol{f}_{h}^{*}(\boldsymbol{\varphi}(\boldsymbol{a}))$,
$\boldsymbol{D}_{r}(\boldsymbol{a}):=\boldsymbol{D}_{h}(\boldsymbol{\varphi}(\boldsymbol{a}))$
and
$\boldsymbol{\eta}_{r}(\boldsymbol{a}):=\boldsymbol{\eta}_{h}(\boldsymbol{\varphi}(\boldsymbol{a}))$.
The ROM is defined by minimizing this residual in the $\Omega_{h}$-norm for
$\frac{d\boldsymbol{a}}{dt}$ given some $\boldsymbol{a}$, this results in the
ROM:

$(\boldsymbol{J}^{T}\Omega_{h}\boldsymbol{J})\frac{d\boldsymbol{a}}{dt}+\boldsymbol{J}^{T}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})=\boldsymbol{J}^{T}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}(\boldsymbol{a}),$

which is indeed well-defined for
$\boldsymbol{a}\in\mathcal{A}\subseteq\mathbb{R}^{r}$ since the Jacobian
$\boldsymbol{J}$ is assumed to be of full-rank on the subset $\mathcal{A}$,
making the mass matrix $(\boldsymbol{J}^{T}\Omega_{h}\boldsymbol{J})$
invertible [4]. In this nonlinear case the ROM is given by the coefficients of
an orthogonal projection of the FOM on the tangent space of $\mathcal{M}$
defined by $T_{u_{r}}\mathcal{M}:=\text{span}(\boldsymbol{J}(\boldsymbol{a}))$
with $\boldsymbol{a}$ such that
$\boldsymbol{u}_{r}=\boldsymbol{\varphi}(\boldsymbol{a})$. This orthogonal
projection is carried out using the $\Omega_{h}$-weighted Moore-Penrose
pseudoinverse
$\boldsymbol{J}^{\dagger}:=(\boldsymbol{J}^{T}\Omega_{h}\boldsymbol{J})^{-1}\boldsymbol{J}^{T}\Omega_{h}$.
Constructing a ROM by projecting the FOM on the tangent space instead of the
reduced manifold itself will result in key differences in our approach
compared to the linear case outlined in [27]: in contrast to the linear case
where $\boldsymbol{J}=\Phi$, $T_{u_{r}}\mathcal{M}$ and $\mathcal{M}$ are no
longer the same space. We introduce
$\boldsymbol{J}^{+}=\boldsymbol{J}^{\dagger}\Omega_{h}^{-1}=(\boldsymbol{J}^{T}\Omega_{h}\boldsymbol{J})^{-1}\boldsymbol{J}^{T}$
and write the ROM in compact form:

$\frac{d\boldsymbol{a}}{dt}+\boldsymbol{J}^{+}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})=\boldsymbol{J}^{+}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}(\boldsymbol{a}).$
(21)

######  Remark 1.

The choice of inner-product spaces for ROMs of hyperbolic systems has recently
come into question [94]. Indeed when $n>1$ the norm
$||\boldsymbol{J}\boldsymbol{\alpha}||_{\Omega_{h}}$ for
$\boldsymbol{\alpha}\in\mathbb{R}^{r}$ is dimensionally inconsistent in
general. It is shown in [94, 69] that dimensionally consistent inner products
that are more appropriate in some sense can improve robustness of the ROMs.
For our approach however it will be important that the ROM is calculated with
the same inner product as used to calculate $\frac{dS_{h}}{dt}$. Therefore, we
will only deal with nondimensionalized conservation laws. Alternatively, our
results can also be applied at an equation-by-equation basis at the cost of
potentially introducing a larger number of generalized coordinates.

######  Remark 2.

A popular approach to construct nonlinear manifold ROMs is the least squares
Petrov-Galerkin (LSPG) method [106, 79]. Using this method a fully discrete
residual is minimized. We have chosen not to use this method because we want
to use the structure of our entropy stable FOM discretization in constructing
entropy stable ROMs. The fully discrete residual minimization approach of LSPG
makes it more difficult to apply this structure.

###  3.2 An entropy stable nonlinear manifold Galerkin ROM

The reduced total entropy functional of the nonlinear manifold ROM is now
defined similarly to the linear case (equation (17)) as:

$S_{r}[\boldsymbol{a}]:=S_{h}[\boldsymbol{\varphi}(\boldsymbol{a})]=\boldsymbol{1}^{T}\bar{\Omega}_{h}\boldsymbol{s}_{r}(\boldsymbol{a}),$

with
$\boldsymbol{s}_{r}(\boldsymbol{a}):=\boldsymbol{s}_{h}(\boldsymbol{\varphi}(\boldsymbol{a}))$.
The reduced total entropy evolution is given by:

$\frac{dS_{r}[\boldsymbol{a}]}{dt}=\boldsymbol{1}^{T}\bar{\Omega}_{h}\frac{d\boldsymbol{s}_{r}}{dt}=\sum_{i}\Delta
x_{i}\boldsymbol{\eta}(\boldsymbol{\varphi}_{i}(\boldsymbol{a}))^{T}\boldsymbol{J}_{i}\frac{d\boldsymbol{a}}{dt}=\boldsymbol{\eta}_{r}^{T}\Omega_{h}\boldsymbol{J}\frac{d\boldsymbol{a}}{dt},$
(22)

where $\boldsymbol{\varphi}_{i}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{n}$ is
the ROM approximation of the conserved variables in the $i$-th cell and
$\boldsymbol{J}_{i}\in\mathbb{R}^{n\times r}$ is the Jacobian matrix of
$\boldsymbol{\varphi}_{i}$ evaluated at $\boldsymbol{a}$. Using (21), the
entropy evolution of the nonlinear manifold Galerkin ROM is:

$\displaystyle\frac{dS_{r}[\boldsymbol{a}]}{dt}$
$\displaystyle=\boldsymbol{\eta}_{r}^{T}\Omega_{h}\boldsymbol{J}\frac{d\boldsymbol{a}}{dt}$
$\displaystyle=-\boldsymbol{\eta}_{r}^{T}\Omega_{h}\boldsymbol{J}\boldsymbol{J}^{+}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})+\boldsymbol{\eta}_{r}^{T}\Omega_{h}\boldsymbol{J}\boldsymbol{J}^{+}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}(\boldsymbol{a})$
$\displaystyle=-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{f}_{r}^{*}(\boldsymbol{a})+\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{D}_{r}(\boldsymbol{a})\Delta_{i}\boldsymbol{\eta}_{r}(\boldsymbol{a}),$
(23)

where the projected entropy variables are defined as
$\tilde{\boldsymbol{\eta}}_{r}=(\Omega_{h}\boldsymbol{J}\boldsymbol{J}^{+})^{T}\boldsymbol{\eta}_{r}=\boldsymbol{J}\boldsymbol{J}^{\dagger}\boldsymbol{\eta}_{r}$,
where $\boldsymbol{J}\boldsymbol{J}^{\dagger}$ is an $\Omega_{h}$-orthogonal
projection on $T_{u_{r}}\mathcal{M}$. It follows from (23) that the reduced
total entropy evolution satisfies an equation that is quite similar to the
total entropy evolution of the FOM. However, instead of the actual entropy
variables $\boldsymbol{\eta}_{r}$ evaluated at the point $\boldsymbol{u}_{r}$
on the reduced manifold $\mathcal{M}$, the inner product is taken with the
projected entropy variables $\tilde{\boldsymbol{\eta}}_{r}$. We would like to
use the entropy conservation condition at this point to show that the inner
product in (23) is zero, but this does not hold because
$\boldsymbol{\eta}(\boldsymbol{u}_{r})\neq\tilde{\boldsymbol{\eta}}_{r}$ in
general. To solve this, we can instead use the invertibility of the entropy
variables to find for what state
$\tilde{\boldsymbol{u}}_{r}\in\mathbb{R}^{N_{h}}$ we do have
$\boldsymbol{\eta}(\tilde{\boldsymbol{u}}_{r})=\tilde{\boldsymbol{\eta}}_{r}$.
If we evaluate our flux at this state instead we can invoke the entropy
conservation condition of the numerical flux to complete the proof of entropy
conservation (or stability). This is exactly what is done in the linear
setting in [27] and leads to our main novelty.

We now present the main novelty of our work. We introduce a novel nonlinear
manifold generalization of the linear entropy projection of [27]. It is given
by:

$\tilde{\boldsymbol{u}}_{r}=\boldsymbol{u}(\boldsymbol{J}\boldsymbol{J}^{\dagger}\boldsymbol{\eta}_{r})=\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}),$
(24)

where the entropy variables $\boldsymbol{\eta}_{r}$ evaluated at the ROM state
$\boldsymbol{u}_{r}$ are projected onto the tangent space
$T_{u_{r}}\mathcal{M}$ instead of the reduced space itself. We note that
projecting the entropy variables on the tangent space $T_{u_{r}}\mathcal{M}$
is a rather natural operation as the entropy variables can be interpreted as
the gradient vector field of the entropy functional $S_{h}$ and thus as
tangent vectors of $\mathbb{R}^{N_{h}}$. The vector
$\tilde{\boldsymbol{u}}_{r}$ is the entropy projected state of the ROM.
Carrying out this entropy projection on the tangent space is necessary as the
projected entropy variables $\tilde{\boldsymbol{\eta}}_{r}$ appearing in the
reduced total entropy evolution equation also are projected on the tangent
space.

A potential issue of our proposed form (24) is that the difference between the
projected entropy variables $\tilde{\boldsymbol{\eta}}_{r}$ and the original
entropy variables $\boldsymbol{\eta}_{r}$ can be very large. Namely, in
similar fashion to the naively constructed linear spaces spanned by $\Phi$
described in subsection 2.3 and [27], for arbitrary $\boldsymbol{\varphi}$ the
entropy variables $\boldsymbol{\eta}_{r}$ may not be well-resolved by the
columns of $\boldsymbol{J}(\boldsymbol{a})$ and thus by the
$\Omega_{h}$-orthogonal projection in (24). This will very likely cause
problems with accuracy of the ROM and the mapping $\boldsymbol{u}$ might not
even be well-defined at $\tilde{\boldsymbol{\eta}}_{r}$. We will propose a
novel method to assure the difference between $\boldsymbol{\eta}_{r}$ and
$\tilde{\boldsymbol{\eta}}_{r}$ remains small in the following section.

First, we continue constructing an entropy stable nonlinear manifold ROM and
perform a Galerkin projection of the FOM at the entropy projected state
$\tilde{\boldsymbol{u}}_{r}$. Doing so we obtain the following ROM:

$\frac{d\boldsymbol{a}}{dt}+\boldsymbol{J}^{+}\Delta_{v}\boldsymbol{f}_{h}^{*}(\tilde{\boldsymbol{u}}_{r})=\boldsymbol{J}^{+}\Delta_{v}D_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r},$
(25)

where we used
$\tilde{\boldsymbol{\eta}}_{r}=\boldsymbol{\eta}_{h}(\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}))$.
Here, the Jacobian matrix is still evaluated at the reduced coordinate
$\boldsymbol{a}\in\mathbb{R}^{r}$ such that
$\boldsymbol{u}_{r}=\boldsymbol{\varphi}(\boldsymbol{a})$ i.e. the non-
projected state. It can be seen that the reduced total entropy evolution is
bounded for this ROM since:

$\displaystyle\frac{dS_{r}[\boldsymbol{a}]}{dt}$
$\displaystyle=-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}))+\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle=(\Delta_{i}\tilde{\boldsymbol{\eta}}_{r})^{T}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r}))-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle=\boldsymbol{1}^{T}\bar{\Delta}_{i}\tilde{\boldsymbol{\psi}}_{r}-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle=0-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{i}^{T}\boldsymbol{D}_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r}$
$\displaystyle\leq 0,$ (26)

where
$\tilde{\boldsymbol{\psi}}_{r}=\boldsymbol{\psi}_{h}(\tilde{\boldsymbol{u}}_{r})$
is the entropy flux potential of the entropy projected state and we were
allowed to invoke the entropy conservation condition of the numerical fluxes.
Clearly, we have:

$\frac{dS_{r}[\boldsymbol{a}]}{dt}=0,$

when no entropy dissipation is present. Note that this approach exactly
recovers the linear approach of [27] when $\boldsymbol{J}=\Phi$, making it a
proper generalization. Thus, by changing the state at which the FOM is
evaluated and projected to the entropy projected state
$\tilde{\boldsymbol{u}}_{r}$, correct total entropy evolution estimates can be
recovered. We added a visualization of the ROM construction with an entropy
projection in Figure 2.

Figure 2: A visualization of the ROM construction. The entropy variables
$\boldsymbol{\eta}_{r}$ are projected on the tangent space
$T_{u_{r}}\mathcal{M}$ to obtain $\tilde{\boldsymbol{\eta}}_{r}$. A new state
$\boldsymbol{u}_{h}$ (not necessarily on the reduced manifold) is found with
entropy variables $\boldsymbol{\eta}_{h}$ such that
$\tilde{\boldsymbol{\eta}}_{r}=\boldsymbol{\eta}_{h}$. We set
$\tilde{\boldsymbol{u}}_{r}=\boldsymbol{u}_{h}$ and project
$-\Delta_{v}\boldsymbol{f}_{h}(\tilde{\boldsymbol{u}}_{r})$ orthogonally on
the tangent space to complete the ROM.

######  Remark 3.

To make relation (26) hold in a fully-discrete setting for Galerkin ROMs,
entropy stable time integration is necessary. This is not trivial as most
methods to satisfy entropy inequalities exactly during time integration
require convexity of the entropy for existence results [102, 118, 70].
Although this is the case for $S_{h}$, $S_{r}$ is not necessarily convex in
the generalized coordinates $\boldsymbol{a}$ unless
$\boldsymbol{\varphi}(\boldsymbol{a})$ is affine i.e.
$\boldsymbol{\varphi}(\boldsymbol{a})=\Phi\boldsymbol{a}+\boldsymbol{u}_{0}$
for some constant $\boldsymbol{u}_{0}\in\mathbb{R}^{N_{h}}$ and
$\Phi\in\mathbb{R}^{N_{h}\times r}$. In this work, we focus on semi-discrete
entropy stability and leave fully discrete entropy-stable ROMs as a suggestion
for future work. In numerical experiments we use sufficiently small time steps
to make entropy errors coming from the time integration negligible, for
details see subsection 5.2.

######  Remark 4.

The entropy conservative hyper-reduction method proposed in [27] does not
generalize to nonlinear spaces as it relies on precomputation of compositions
of linear operators. In the nonlinear case precomputation is not possible due
to the changing tangent space. An entropy conservative hyper-reduction method
suitable for nonlinear reduced spaces is also a suggestion for future work.

###  3.3 Tangent space enrichment

To arrive at the correct entropy estimate (26) we carried out the entropy
projection (24). Though the evolution of the entropy then satisfies a correct
estimate, it is not clear whether the ROM solution itself remains accurate.
Particularly, the difference between the entropy projected state (24) and the
original state (19) can be very large. To see this we consider the entropy
projection error:

$\varepsilon_{s}:=||\boldsymbol{u}_{r}-\tilde{\boldsymbol{u}}_{r}||_{\Omega_{h}}.$

Assuming the mapping $\boldsymbol{u}$ from entropy variables to conservative
variables is sufficiently smooth, using the mean-value theorem we can bound
this term as follows:

$\displaystyle||\boldsymbol{u}_{r}-\tilde{\boldsymbol{u}}_{r}||_{\Omega_{h}}$
$\displaystyle=||\boldsymbol{u}_{r}-\boldsymbol{u}(\boldsymbol{J}\boldsymbol{J}^{\dagger}\boldsymbol{\eta}_{r})||_{\Omega_{h}}$
$\displaystyle=||\boldsymbol{u}_{r}-\left[\boldsymbol{u}_{r}-\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{\eta}}(\boldsymbol{\theta})(I-\boldsymbol{J}\boldsymbol{J}^{\dagger})\boldsymbol{\eta}_{r}\right]||_{\Omega_{h}}\quad\theta_{i}\in[(\boldsymbol{u}_{r})_{i},(\tilde{\boldsymbol{u}}_{r})_{i}]\hbox{
}\forall i$
$\displaystyle\leq||\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{\eta}}(\boldsymbol{\theta})||_{\Omega_{h}}||(I-\boldsymbol{J}\boldsymbol{J}^{\dagger})\boldsymbol{\eta}_{r}||_{\Omega_{h}}$
$\displaystyle=||\left(\frac{\partial^{2}\boldsymbol{s}}{\partial\boldsymbol{u}^{2}}\right)^{-1}(\boldsymbol{\theta})||_{\Omega_{h}}||(I-\boldsymbol{J}\boldsymbol{J}^{\dagger})\boldsymbol{\eta}_{r}||_{\Omega_{h}},$

where we used the definition of the entropy variables in the last line and
used the induced operator $\Omega_{h}$-norm for the Hessian matrix of the
entropy function. It can be seen that there are two contributions to this
bound. There is one based on the model, specifically on the Hessian of the
entropy, and one given by the projection error of $\boldsymbol{\eta}_{r}$ on
the tangent space $T_{u_{r}}\mathcal{M}$. As long as the entropy is a convex
function at the mean value $\boldsymbol{\theta}$, the contribution of the
model-based term is bounded by a term involving the inverse of the smallest
eigenvalue of the entropy Hessian. We have little influence over this term.
Specifically, this term can be large when the entropy is close to being non-
convex. We do have control over the projection error. The magnitude of this
term is controlled by the choice of reduced space. For a general reduced space
constructed to contain solution snapshots this term can be very large, since
the columns of the Jacobian $\boldsymbol{J}$ can be close to orthogonal to
$\boldsymbol{\eta}_{r}$ while the standard nonlinear manifold ROM (21) works
fine. In the linear case, Chan [27] solved this problem by enriching the
snapshot data to construct $\Phi$ with snapshots of the entropy variables.
This lowered the projection error contribution to the bound on
$\varepsilon_{s}$ since $\boldsymbol{J}=\Phi$ in this case. However, for the
general nonlinear case, $\boldsymbol{J}$ is not the same as the reduced space
itself and we can no longer construct our reduced space to also contain the
entropy variables to keep the projection error low. Instead, we need a
different approach and therefore we propose a novel method to which we refer
as tangent space enrichment.

The key idea of tangent space enrichment is to construct an $r+1$-dimensional
manifold $\hat{\mathcal{M}}\subset\mathbb{R}^{N_{h}}$ from the original
$r$-dimensional manifold $\mathcal{M}$ by a ‘lifting’ operation. Consequently,
we use this new manifold for the ROM instead. This lifting operation is
defined so that the original manifold $\mathcal{M}$ is a subset of the lifted
manifold $\hat{\mathcal{M}}$, i.e. $\mathcal{M}\subset\hat{\mathcal{M}}$. Most
importantly however, for all points
$\boldsymbol{u}_{r}\in\mathcal{M}\subset\hat{\mathcal{M}}$ the lifting
operation is constructed such that $\boldsymbol{\eta}_{r}\in
T_{u_{r}}\hat{\mathcal{M}}$. This means that the entropy variable projection
error is precisely zero at the points contained in the old manifold when
projecting $\boldsymbol{\eta}_{r}$ on the tangent space of the new manifold.
This assures that the entropy projection is accurate for the points
$\boldsymbol{u}_{r}\in\mathcal{M}\subset\hat{\mathcal{M}}$ when using tangent
space enrichment.

We motivate this approach over a more straightforward generalization of Chan’s
snapshots enrichment [27] method by the following. Nonlinear reduced spaces
are often constructed iteratively by minimizing some loss function. A
nonlinear version of Chan’s enrichment method would require, for a given
$\boldsymbol{a}$, fitting $\boldsymbol{\varphi}(\boldsymbol{a})$ to a snapshot
$\boldsymbol{u}_{h}$ whilst the Jacobian $\boldsymbol{J}(\boldsymbol{a})$ has
a low entropy projection error $\varepsilon_{s}$. The construction of
$\mathcal{M}$ would therefore require including terms based on
$\boldsymbol{J}$ in the loss function. This can be very expensive and
generally will not exactly embed the entropy variables in the tangent space at
the appropriate points. As will be discussed, our method requires no extra
effort in constructing $\boldsymbol{\varphi}$ and contrary to the
straightforward approach exactly enriches the tangent spaces with the correct
entropy variables.

Our novel tangent space enrichment is defined by the following
parameterization for $\hat{\mathcal{M}}$:

$\hat{\boldsymbol{\varphi}}(\boldsymbol{a},\alpha)=\boldsymbol{\varphi}(\boldsymbol{a})+\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))\alpha,$
(27)

here $\boldsymbol{\varphi}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{N_{h}}$ is the
parameterization of the original manifold $\mathcal{M}$,
$\boldsymbol{a}\in\mathbb{R}^{r}$ are the $r$ reduced coordinates associated
to the original parameterization $\boldsymbol{\varphi}$ and
$\alpha\in\mathbb{R}$ is the $r+1$-th reduced coordinate associated to the
lifting operation. As in remark 1, here we see another reason for the
importance of non-dimensionalization. Namely, on dimensional grounds the
expression (27) does not make sense if $\boldsymbol{\varphi}$ and
$\boldsymbol{\eta}$ are not suitably normalized.

The new parameterization can be interpreted as follows. Given any point
$\boldsymbol{\varphi}(\boldsymbol{a})\in\mathcal{M}$ we generate new points
$\hat{\boldsymbol{u}}_{r}\in\hat{\mathcal{M}}$ by lifting the new points from
the point $\boldsymbol{\varphi}(\boldsymbol{a})$ in the direction of
$\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))$ by a distance
$||\hat{\boldsymbol{\varphi}}(\boldsymbol{a},\alpha)-\boldsymbol{\varphi}(\boldsymbol{a})||=||\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))||\cdot|\alpha|$.
Note that at $\alpha=0$ we do not lift the point at all, resulting in a point
at $\boldsymbol{\varphi}(\boldsymbol{a})$; in other words, the original
manifold $\mathcal{M}$ is the set
$\hat{\boldsymbol{\varphi}}(\mathbb{R}^{r},\alpha=0)$. The lifting operation
is visualized in Figure 3. The Jacobian matrix of the new parameterization
(27), whose columns span the tangent space $T_{\hat{u}_{r}}\hat{\mathcal{M}}$
of the lifted manifold $\hat{\mathcal{M}}$ at the point
$\hat{\boldsymbol{u}}_{r}\in\hat{\mathcal{M}}$, is given by:

$\displaystyle\hat{\boldsymbol{J}}(\boldsymbol{a},\alpha)$
$\displaystyle=\begin{bmatrix}\frac{\partial\hat{\boldsymbol{\varphi}}}{\partial\boldsymbol{a}}&\frac{\partial\hat{\boldsymbol{\varphi}}}{\partial\alpha}\end{bmatrix}$
(28) $\displaystyle=\begin{bmatrix}\left(I+\alpha\hbox{
}\frac{\partial\boldsymbol{\eta}}{\partial\boldsymbol{u}}(\boldsymbol{\varphi}(\boldsymbol{a}))\right)\boldsymbol{J}(\boldsymbol{a})&\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))\end{bmatrix},$

where $I\in\mathbb{R}^{N_{h}\times N_{h}}$ is the identity matrix and
$\frac{\partial\boldsymbol{\eta}}{\partial\boldsymbol{u}}(\boldsymbol{\varphi}(\boldsymbol{a}))=\frac{\partial^{2}\boldsymbol{s}_{h}}{\partial\boldsymbol{u}^{2}}(\boldsymbol{\varphi}(\boldsymbol{a}))\succeq
0$ is a sparse SPSD $2n-1$-diagonal matrix containing components of the
Hessian of the local entropy value with respect to the solution on each
diagonal. Note that the derivative with respect to the $r+1$-th tangent space
enrichment coordinate $\alpha$ is exactly
$\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))$. Furthermore, on the
old manifold $\mathcal{M}$ associated to $\alpha=0$ the matrix
$\frac{\partial\hat{\boldsymbol{\varphi}}}{\partial\boldsymbol{a}}$ is equal
the original Jacobian $\boldsymbol{J}(\boldsymbol{a})$. At the points
$\hat{\boldsymbol{\varphi}}(\boldsymbol{a},\alpha=0)=\boldsymbol{\varphi}(\boldsymbol{a})$
we have thus exactly enriched the tangent space with the entropy variables
$\boldsymbol{\eta}_{r}=\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))$
at those points. This is the direct result of the lifting operation. This can
be seen from the enriched parameterization (27). Lifting a point from
$\boldsymbol{\varphi}(\boldsymbol{a})$ by changing $\alpha$ while keeping
$\boldsymbol{a}$ constant, moves a point in the direction tangent to
$\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))$. As a consequence
$\boldsymbol{\eta}(\boldsymbol{\varphi}(\boldsymbol{a}))$ appears as a tangent
vector in the enriched Jacobian (28).

Figure 3: A visualization of tangent space enrichment. The dark purple curve
is the original manifold $\mathcal{M}$, the light purple region is a section
of the enriched manifold $\mathcal{M}^{*}$. $\mathcal{M}^{*}$ is constructed
by lifting points from $\varphi(\boldsymbol{a})$ along lines in the direction
of $\boldsymbol{\eta}(\varphi(\boldsymbol{a}))$.

Given the tangent space enrichment, the ROM is constructed in a similar
fashion as in the previous section, but using the enriched lifted manifold
$\hat{\mathcal{M}}$. The $r+1$-th tangent space enrichment coordinate $\alpha$
is simply treated as an additional reduced coordinate. The ROM thus takes the
form:

$\frac{d}{dt}\begin{bmatrix}\boldsymbol{a}\\\
\alpha\end{bmatrix}+\hat{\boldsymbol{J}}^{+}\Delta_{v}\boldsymbol{f}_{h}^{*}(\hat{\tilde{\boldsymbol{u}}}_{r})=\hat{\boldsymbol{J}}^{+}\Delta_{v}\boldsymbol{D}_{h}(\hat{\tilde{\boldsymbol{u}}}_{r})\Delta_{i}\hat{\tilde{\boldsymbol{\eta}}}_{r},$

where:

$\hat{\tilde{\boldsymbol{u}}}_{r}=\boldsymbol{u}(\hat{\boldsymbol{J}}\hat{\boldsymbol{J}}^{\dagger}\hat{\boldsymbol{\eta}}_{r}),$

and
$\hat{\tilde{\boldsymbol{\eta}}}_{r}=\hat{\boldsymbol{J}}\hat{\boldsymbol{J}}^{\dagger}\boldsymbol{\eta}(\hat{\boldsymbol{\varphi}}(\boldsymbol{a},\alpha))$.
Since the proof of entropy stability for our ROM is independent of the
manifold parameterization, this ROM is still entropy stable.

##  4 Rational polynomial manifolds

###  4.1 Background

In this article we are interested in systems that can exhibit strong spatial
gradients that are moving in time. Many existing data compression methods for
reduced manifold construction are not well suited for these types of systems.
Linear data compression methods like proper orthogonal decomposition (POD) and
general reduced basis methods fail because the moving gradients imply that the
data is often of very high rank. This indicates that the data is not well
represented in low-dimensional linear subspaces, which is the fundamental
assumption of linear approaches. Nonlinear data compression methods offer a
potential solution to this problem by instead compressing the highly nonlinear
data on nonlinear reduced manifolds. In model reduction different
parameterizations have become popular, in particular the decoder part of
autoencoder neural networks [105, 106, 79, 72, 121, 37, 9, 31] and
multivariate quadratic polynomials [8, 111, 17, 47, 66]. However, in the
vicinity of strong gradients these nonlinear methods can suffer from
oscillations [47]. These oscillations may be difficult or impossible to
remove. This is a problem for Galerkin projection-based ROMs which can be
sensitive to errors in the solution compression [23]. To accurately assess the
performance of our novel entropy stable manifold Galerkin ROM there is thus a
need for nonlinear data compression methods that are more capable of dealing
with large and moving spatial gradients, particularly without significant
oscillations. Recently, neural networks with discontinuous activation
functions have been proposed [36], however training these networks can be
cumbersome. Furthermore, registration based approaches [116, 89] have been
very effective, but have not yet been applied in the context of nonlinear
manifold ROMs similar to ours. In this research we will propose a novel
reduced manifold parameterization method based on rational polynomials, that
is far less oscillatory around strong spatial gradients than the previously
mentioned methods (neural networks and quadratic approaches), but is still
equally interpretable as the recently proposed quadratic manifolds.

###  4.2 Pole-free rational quadratic manifolds

We now give a description of rational polynomial manifolds. A rational
polynomial manifold is the element-wise ratio of two polynomial manifolds:

$\boldsymbol{\varphi}(\boldsymbol{a})=\frac{\sum_{i=1}^{p_{\text{num}}}\boldsymbol{H}^{i}:\boldsymbol{a}^{\otimes
i}+\boldsymbol{u}_{\text{ref}}}{\sum_{i=1}^{p_{\text{den}}}\boldsymbol{G}^{i}:\boldsymbol{a}^{\otimes
i}+\boldsymbol{1}},$ (29)

here, $\boldsymbol{H}^{i},\boldsymbol{G}^{i}\in\mathbb{R}^{N_{h}\times
r\times...\times r}$ are $(i+1)^{\text{th}}$-order tensors with the first axis
of size $N_{h}$ and $i$ axes of length $r$, $\boldsymbol{a}^{\otimes i}$ is
the $i$-fold outer product such that for example $(\boldsymbol{a}^{\otimes
3})_{ijk}=a_{i}a_{j}a_{k}$, $\boldsymbol{H}^{i}:\boldsymbol{a}^{\otimes i}$
denotes summation of the components of $\boldsymbol{a}^{\otimes i}$ and the
components of slices along the first axis of $\boldsymbol{H}^{i}$, again as
example $\left(\boldsymbol{H}^{3}:\boldsymbol{a}^{\otimes
3}\right)_{i}=\sum_{j,k,l=0}^{r-1}\left(\boldsymbol{H}^{3}\right)_{ijkl}a_{j}a_{k}a_{l}$.
Furthermore, we have $\boldsymbol{u}_{\text{ref}}\in\mathbb{R}^{N_{h}}$ and we
consider division of two vectors element-wise. The constant vector in the
denominator has been set to one without loss of generality. The expression
(29) generalizes polynomial manifolds [47, 8], in the sense that those are
recovered by setting $p_{\text{den}}=0$. This shows that rational polynomial
manifold encapsulate a larger class of functions than polynomial manifolds. By
introducing a polynomial in the denominator and allowing it to rapidly and
smoothly approach zero for small changes in $\boldsymbol{a}$ we can obtain
very fast and smooth increases in the function value of $\varphi$ without
oscillations. This gives us the opportunity to model steep gradients in the
snapshot data which may be present in the form of advected shocks. Since
higher-order tensors can become quite expensive to deal with, we restrict our
attention to rational quadratic manifolds by setting
$p_{\text{num}}=p_{\text{den}}=2$.

To compromise between efficiency and expressiveness we will take
$p_{\text{num}}=p_{\text{den}}=2$. The $i$-th component of the vector-valued
function output $\varphi_{i}(\boldsymbol{a})$ can then be written as:

$\varphi_{i}(\boldsymbol{a})=\frac{\boldsymbol{a}^{T}\boldsymbol{H}_{i}^{2}\boldsymbol{a}+\boldsymbol{H}_{i}^{1}\boldsymbol{a}+(\boldsymbol{u}_{\text{ref}})_{i}}{\boldsymbol{a}^{T}\boldsymbol{G}_{i}^{2}\boldsymbol{a}+\boldsymbol{G}_{i}^{1}\boldsymbol{a}+1},$
(30)

where $\boldsymbol{H}_{i}^{2},\boldsymbol{G}_{i}^{2}\in\mathbb{R}^{r\times r}$
are the $i$-th slices along the first axes of $\boldsymbol{H}^{2}$ and
$\boldsymbol{G}^{2}$, respectively and
$\boldsymbol{H}_{i}^{1},\boldsymbol{G}_{i}^{1}\in\mathbb{R}^{1\times r}$ are
the $i$-th rows of $\boldsymbol{H}^{1}$ and $\boldsymbol{G}^{1}$,
respectively. Since the matrices are only used in quadratic forms, we can,
without loss of generality, assume
$\boldsymbol{H}_{i}^{2},\boldsymbol{G}_{i}^{2}$ to be symmetric. An obvious
concern with expressions of this form is the occurrence of spurious poles,
i.e. unwanted division by zero. We avoid this issue for the case
$p_{\text{num}}=p_{\text{den}}=2$ and all $\boldsymbol{a}\in\mathbb{R}^{r}$ by
constraining the quadratic form in the denominator to be positive semi-
definite and setting the linear term to zero:

$\boldsymbol{G}_{i}^{2}\succeq 0,\quad\boldsymbol{G}_{i}^{1}=0,\quad\forall
i\in\\{0,...,N_{h}\\}.$

It is easily seen that, in this case, the denominator is never less than 1.
Consequently, spurious poles cannot occur for any real $\boldsymbol{a}$.
Removing the linear term has not resulted in a significant loss in accuracy in
our numerical experiments. The full $\boldsymbol{\varphi}$ is then given as:

$\boldsymbol{\varphi}(\boldsymbol{a})=\frac{\boldsymbol{H}^{2}:[\boldsymbol{a}\otimes\boldsymbol{a}]+\boldsymbol{H}^{1}\boldsymbol{a}+\boldsymbol{u}_{\text{ref}}}{\boldsymbol{G}:[\boldsymbol{a}\otimes\boldsymbol{a}]+\boldsymbol{1}},\quad\boldsymbol{G}_{i}\succeq
0\hbox{ }\forall i,$ (31)

since there is no linear term in the denominator we write $\boldsymbol{G}$
instead of $\boldsymbol{G}^{2}$. To construct manifold Galerkin ROMs we will
require the Jacobian matrix of this expression - see equations (21) and (28).
The Jacobian matrix is given by the following:

$\frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}}=\frac{2\boldsymbol{H}^{2}\cdot\boldsymbol{a}+\boldsymbol{H}^{1}}{\left(\boldsymbol{G}:[\boldsymbol{a}\otimes\boldsymbol{a}]+\boldsymbol{1}\right)\otimes\boldsymbol{1}}-\left[\frac{\boldsymbol{H}^{2}:[\boldsymbol{a}\otimes\boldsymbol{a}]+\boldsymbol{H}^{1}\boldsymbol{a}+\boldsymbol{u}_{\text{ref}}}{\left(\boldsymbol{G}:[\boldsymbol{a}\otimes\boldsymbol{a}]+\boldsymbol{1}\right)^{2}}\otimes\boldsymbol{1}\right]\circ(2\boldsymbol{G}\cdot\boldsymbol{a}),$

where the operation
$2\boldsymbol{H}^{2}\cdot\boldsymbol{a}\in\mathbb{R}^{N_{h}\times r}$
indicates slice-wise matrix multiplication, i.e. for the $i$-th row it holds
that
$(2\boldsymbol{H}^{2}\cdot\boldsymbol{a})_{i}=2\boldsymbol{H}_{i}^{2}\boldsymbol{a}$,
due to symmetry of $\boldsymbol{H}_{i}^{2},\boldsymbol{G}_{i}$ the order of
axes is not relevant, division of matrices is understood element-wise and
$\circ$ is the Hadamard matrix multiplication operator.

###  4.3 Manifold construction

We will determine the coefficient tensors in (31) from data. Like for
quadratic manifolds [8, 47], it holds for rational quadratic manifolds that
the coefficients in two different slices of the coefficient tensors are
independent, as can be seen in (30). Consequently, the coefficients can be
determined purely from the parametric and temporal behaviour of the data in
the specific cell and solution variable associated to the $i$-th component
$\varphi_{i}(\boldsymbol{a})$ of $\boldsymbol{\varphi}(\boldsymbol{a})$. The
task of fitting a rational manifold thus reduces to fitting an expression (30)
to data for each cell and solution variable.

In the spirit of quadratic manifolds we will compress the snapshot data
defined by:

$X=[\boldsymbol{u}_{h}(t^{0}),\boldsymbol{u}_{h}(t^{1}),...,\boldsymbol{u}_{h}(t^{n_{s}-1})]\in\mathbb{R}^{N_{h}\times
n_{s}},$

where $n_{s}\in\mathbb{N}$ is the number of snapshots, as the coefficients of
their projection on the first $r$ left singular vectors of $X$ given in
$\Phi\in\mathbb{R}^{N_{h}\times r}$. Following this we define
$A:=(\Phi^{T}X)^{T}\in\mathbb{R}^{n_{s}\times r}$. We aim to find the
coefficients
$\boldsymbol{G}_{i},\boldsymbol{H}_{i}^{2},\boldsymbol{H}_{i}^{1},(\boldsymbol{u}_{\text{ref}})_{i}$
such that for all $j=0,...,n_{s}-1$ we have
$\varphi_{i}(\boldsymbol{a}_{j})\approx X_{ij}$, where $\boldsymbol{a}_{j}$ is
the $j$-th row of $A$ written as a column vector. Defining the $i$-th row of
$X$ as $\boldsymbol{y}^{i}\in\mathbb{R}^{n_{s}}$, the optimization procedure
for the coefficients
$\boldsymbol{G}_{i},\boldsymbol{H}_{i}^{2},\boldsymbol{H}_{i}^{1},(\boldsymbol{u}_{\text{ref}})_{i}$
is formulated as:

$\boldsymbol{G}_{i},\boldsymbol{H}_{i}^{2},\boldsymbol{H}_{i}^{1},(\boldsymbol{u}_{\text{ref}})_{i}=\operatorname*{arg\,min}_{\boldsymbol{G}\in\mathbb{S}^{r}_{+},\boldsymbol{H}\in\mathbb{S}^{r},\boldsymbol{h}\in\mathbb{R}^{r},u\in\mathbb{R}}\left|\left|\boldsymbol{y}^{i}-\frac{(A\boldsymbol{H}\circ
A)\boldsymbol{1}+A\boldsymbol{h}+u}{(A\boldsymbol{G}\circ
A)\boldsymbol{1}+1}\right|\right|^{2},$ (32)

though the sets $\mathbb{S}^{r}_{+}$ and $\mathbb{S}^{r}$ are convex subsets
of $\mathbb{R}^{r\times r}$ the objective function is non-convex due to the
division operation.

A popular approach is to linearize this nonlinear optimization problem [59, 3]
by multiplying with the denominator. This results in a convex semi-definite
program which can be solved very efficiently with e.g. interior point methods.
However, the optimum value of this linearized problem is generally not the
same as the nonlinear problem (32). Iterative approaches that attempt to
somehow refine the solution of the linearized problem and that can potentially
be supplemented with our semi-definite constraint exist [59, 3, 53]. However,
convergence to (local) minima of the nonlinear problem is generally not
guaranteed, nor is convergence in general [112, 60]. For this reason, we will
stick to the fully nonlinear and expensive optimization problem (32).
Nonetheless, we believe that the linearized approach holds significant promise
and that it will be crucial for future work in order to scale the approach to
larger meshes and datasets. Finally, we will implement the semi-definite
constraint in this fully nonlinear setting by optimizing for the Cholesky
decomposition of $\boldsymbol{G}_{i}=L_{i}L_{i}^{T}$ [4]:

$L_{i},\boldsymbol{H}_{i}^{2},\boldsymbol{H}_{i}^{1},(\boldsymbol{u}_{\text{ref}})_{i}=\operatorname*{arg\,min}_{L\in\mathbb{L}^{r},\boldsymbol{H}\in\mathbb{S}^{r},\boldsymbol{h}\in\mathbb{R}^{r},u\in\mathbb{R}}\left|\left|\boldsymbol{y}^{i}-\frac{(A\boldsymbol{H}\circ
A)\boldsymbol{1}+A\boldsymbol{h}+u}{(AL\circ
AL)\boldsymbol{1}+1}\right|\right|^{2},$ (33)

where $\mathbb{L}^{r}\subset\mathbb{R}^{r\times r}$ is the vector subspace of
lower triangular $r\times r$ matrices. As initial guess we can either use
$L_{i-1},\boldsymbol{H}_{i-1}^{2},\boldsymbol{H}_{i-1}^{1},(\boldsymbol{u}_{\text{ref}})_{i-1}$
if available and corresponding to the same solution variable, or otherwise
simply vectors or tensors consisting of “ones”. We will carry out the fitting
procedure using the JAXFit package [58] for GPU-accelerated nonlinear least
squares solutions.

##  5 Numerical experiments

To show that our entropy stable manifold Galerkin ROMs satisfy appropriate
semi-discrete entropy inequalities we will perform numerical experiments on a
range of one-dimensional nonlinear conservation laws. We will carry out the
experiments using the rational quadratic manifolds proposed in section 4. We
will also compare the ability of the rational quadratic manifolds to compress
convection dominated data to that of linear POD-based methods and quadratic
manifolds [8, 66, 47]. We do not compare against neural network based
approaches [106, 79, 9] since in our experience they struggle with
approximating discontinuities, and require careful hyperparameter tuning to
give reasonable results. The underlying FOM discretizations will be the
existing TeCNO schemes of [44], so that we will only mention some aspects of
the discretizations but for details we will refer to [44]. We will start with
the inviscid Burgers equation in subsection 5.1, then we will treat the
shallow water equations in subsection 5.2 and finally we will treat the
compressible Euler equations with ideal thermodynamics in subsection 5.3. We
test different aspects of the ROM with the different test cases, an overview
of the different test purposes has been provided in Table 1.

Experiment | Purpose
---|---
Inviscid Burgers | Manifold accuracy
Shallow water | Entropy conservation properties
Compressible Euler | Impact of entropy projection and tangent space enrichment
Table 1: Overview of numerical experiments.

The experiments have been implemented using the JAX library [15] in Python,
which allows for automatic differentiation to compute Jacobian matrices and,
where possible, accelerated computing using an Nvidia A2000 laptop GPU.

###  5.1 Inviscid Burgers equation

We will use this experiment mainly to highlight our proposed rational
quadratic manifolds when compared to existing (‘standard’) Galerkin ROMs on
different types of manifolds. We will already include the entropy stable ROM
(25) here, but the focus on the role of entropy stability will be in the next
test cases.

The inviscid Burgers equation is given by:

$\frac{\partial u}{\partial t}+\frac{\partial}{\partial
x}\left(\frac{u^{2}}{2}\right)=0,$ (34)

with conserved variable $u:\Omega\times[0,T]\rightarrow\mathbb{R}$. As
continuous and discrete entropy we take [82]:

$\mathcal{S}[u]=\frac{1}{2}\int_{\Omega}u^{2}dx,$

and

$S_{h}[\boldsymbol{u}_{h}]=\frac{1}{2}||\boldsymbol{u}_{h}||_{\Omega_{h}}^{2}.$

The reduced entropy functional follows in a straightforward fashion from the
discrete entropy functional. Using these specific entropies we have for the
entropy variables:

$\eta(u)=u.$

Consequently, the manifold parameterization with TSE is given in the
particularly simple form:

$\boldsymbol{\varphi}^{*}(\boldsymbol{a},\alpha)=(1+\alpha)\cdot\boldsymbol{\varphi}(\boldsymbol{a}).$

An entropy conservative flux is given by [44]:

$f_{i+1/2}=\frac{u_{i+1}^{2}+u_{i+1}u_{i}+u_{i}^{2}}{6},$

and we use a local Lax-Friedrichs-type of entropy dissipation operator [82]:

$D_{i+1/2}(\boldsymbol{u}_{h})=\max(|u_{i+1}|,|u_{i}|).$

We will discretize (34) on a domain $\Omega=\mathbb{T}([0,L])$ of length $L=1$
using a numerical grid consisting of 300 cells. Discretization in time will be
done using the classical Runge-Kutta 4 (RK4) method [50] with a time step of
size $\Delta t=0.001$. We will integrate in time until $T=1$. Solution
snapshots are captured after every 5 timesteps resulting in $n_{s}=201$. We
perform two tests: in the first, we compress the data to a reduced dimension
of $r=15$ for all manifolds, and in the second we compress the data to a
reduced dimension such that the reconstruction errors (to be defined later)
are of the same order as the rational quadratic manifold with $r=15$. The
initial condition is a simple offset sine wave:

$u_{0}(x)=\sin(2\pi x)+1.$

Figure 4: Normalized singular values of the inviscid Burgers equation data.

The KnW decay (1) for this system is very slow as is evident from the
normalized singular values depicted in Figure 4. Defining the relative
information content (RIC) as in [76] we have $\text{RIC}\approx 99.5\%$ for
$r=15$.

Figure 5: Original data and reconstructions for $r=15$ in space-time plots.
Figure 6: Space–time errors $\varepsilon_{xt}$ of the reconstructions.

We will first compare the reconstruction accuracy of our proposed rational
quadratic manifold to existing quadratic manifold [8] and POD linear [109]
manifold approaches for the solution data of the FOM with $r=15$. To do this
we will save the matrix of generalized coordinates
$A=\Phi^{T}X\in\mathbb{R}^{r\times n_{s}}$ associated to the snapshots in $X$.
Note that these coordinates form the reduced representation for all manifolds
since all of the tested manifolds are constructed from the POD compression of
the data. In turn we will try to reconstruct the snapshots in $X$ from their
reduced representations in $A$. We will construct the quadratic manifold as in
[8] with a manually determined regularization coefficient $\lambda=0.5$
($\alpha$ in (27) of [8]) and the rational manifold using the fully nonlinear
curve-fitting approach outlined in the previous section. In Figure 5 we
display the reconstruction in addition to the original data using an $x-t$
plot. Furthermore, in Figure 6 we plot the local error in space–time defined,
with some abuse of notation, as:

$\varepsilon_{xt}=\boldsymbol{\varphi}(A)-X.$

It can be seen that for $r=15$, the reconstruction accuracy of both the
quadratic as POD linear manifold is poor, whereas the reconstruction of the
rational quadratic manifold is visually nearly identical to the data. Indeed,
the largest local error of the rational quadratic manifold is at most
approximately $0.003$ which is two orders of magnitude lower than the largest
errors of the quadratic and linear manifolds (approximately $0.4$ and $0.6$
respectively). The sources of error for the linear and quadratic manifolds are
predominantly oscillations around the moving shock discontinuity as can be
seen in Figure 6. This shows the poor performance of these methods for such
problems. The rational quadratic manifold also oscillates around the shock,
but with a much smaller amplitude, indicating that it is better-suited for
shock-dominated problems. The accuracy of the rational quadratic manifold is
much higher than the quadratic and POD linear manifolds.

We note that the increased accuracy comes at the cost of a slow fitting
procedure. The precise fitting times and maximum absolute space-time errors,
$\varepsilon_{xt}^{\text{max}}:=\max_{i,j}|\varepsilon_{xt}|_{i,j}$, have been
displayed in Table 2.

| $t_{\text{fit}}$ | $t_{\text{online}}$ | $\varepsilon_{xt}^{\text{max}}$
---|---|---|---
entropy stable rational quadratic ($r=15)$ | 1851 s | 7.01 s | $3.36\cdot 10^{-3}$
rational quadratic ($r=15)$ | 1851 s | 6.74 s | $3.36\cdot 10^{-3}$
POD linear ($r=15$) | 0.107 s | 0.543 s | 0.675
POD linear ($r=160$) | 0.107 s | 0.622 s | $2.61\cdot 10^{-3}$
quadratic ($r=15$) | 0.964 s | 6.14 s | 0.456
quadratic ($r=150$) | 122 s | 121 s | $4.70\cdot 10^{-3}$
Table 2: Accuracy and computational cost for fitting Burgers solution with
different manifold parameterizations

When we construct the POD linear manifold and quadratic manifold to an
accuracy of $\varepsilon_{xt}^{\text{max}}\approx 3\cdot 10^{-3}$ we see that
we require approximately $r=160$ and $r=150$, respectively. The
reconstructions are given in Figure 7 and the largest space–time errors
$\varepsilon_{xt}$ are $\mathcal{O}(4\cdot 10^{-3})$. The changes in fitting
time have also been denoted in Table 2. A large increase can be observed for
quadratic manifolds, while the linear POD stays constant as the implementation
calculates all $n_{s}$ singular vectors at once.

Figure 7: Reconstructions for $\varepsilon_{xt}^{\text{max}}\approx 3\cdot
10^{-3}$ in space-time plots.

We continue to consider the ROM performance and accuracy in more detail. We
compare the rational quadratic manifold ROM in entropy stable (25) (ES-ROM)
and generic (21) (RQ-ROM) form to a linear manifold POD-Galerkin ROM (L-ROM)
and a quadratic manifold Galerkin ROM (Q-ROM). We will make one comparison of
the ROMs using the previously obtained manifolds with $r=15$ and another
comparison where $r$ is chosen such that each manifold has approximately the
same space–time reconstruction error $\varepsilon_{xt}^{\text{max}}\approx
3\cdot 10^{-3}$. The initial conditions for the simulations with $r=15$ will
be taken as the first column of the matrix $A=\Phi^{T}X$ and the entropy
stable form of the rational quadratic ROM will have $\alpha=0$ at $t=0$. For
the manifolds that have approximately the same accuracy we will take the first
columns of the matrices $A_{r}=\Phi_{r}^{T}X\in\mathbb{R}^{r\times n_{s}}$
with $r=160$ and $r=150$ for the linear and quadratic manifolds, respectively.
Here, $\Phi_{r}\in\mathbb{R}^{N_{h}\times r}$ are the $r$ first singular
vectors of $X$. We plot the temporal evolution of the total error norm:

$\varepsilon_{u}(t):=||\boldsymbol{u}_{h}(t)-\boldsymbol{u}_{r}(t)||_{\Omega_{h}},$

and the ideal linear projection error (L-ideal):

$\varepsilon_{\text{proj}}(t):=||(I-\Pi_{r_{\text{lin}}})\boldsymbol{u}_{h}(t)||_{\Omega_{h}},$

where $\Pi_{r_{\text{lin}}}:\mathbb{R}^{N_{h}}\rightarrow\mathcal{V}$ projects
on the $r_{\text{lin}}$-dimensional reduced space of the respective linear
ROMs in the different experiments. The ideal projection error forms a lower
bound for the linear POD-Galerkin ROM error. To measure computational
performance we will track the runtimes of the online phases which we denote
$t_{\text{online}}$. The results for the simulations with constant $r$ are
shown in Figure 9 and the results for approximately constant space–time error
are given in Figure 9. We also show the spatial profile of the solution at
$t=T$ as predicted by the different ROMs and different manifolds in Figure 11
for constant $r$ and Figure 11 for constant space–time error, respectively.
The rational manifold ROMs clearly outperform the others in case of constant
$r$ and the differences between the results of the rational manifolds
themselves are nearly zero. Steep increases in errors occur for all ROMs upon
formation of the shock which indicates this is indeed a difficult instant of
the flow for the reduced manifolds to fit to the data. For the simulation with
manifolds with constant error the performance of the ROMs is more equal, with
the linear and quadratic ROMs only suffering of some oscillations in the
spatial profile. The oscillations occur at moments when the ideal projections
error $\varepsilon_{\text{proj}}$ is also oscillatory in time. The rational
quadratic manifolds do not suffer from oscillations. Because of the large
reduced spaces required to obtain approximately the same reconstruction errors
as the rational manifolds, the linear and quadratic manifold ROMs have more
expensive online phases than when tested at constant $r$. This is especially
notable for quadratic manifolds where computing the Jacobian and its Moore-
Penrose pseudoinverse contribute heavily to the increase in cost. At constant
$r$ the quadratic manifold ROM and the rational manifold ROMs are nearly
equally fast showing that the Jacobian of the quadratic manifold
parameterization and of the rational manifolds with and without enrichment are
nearly equally expensive to compute. From this experiment we conclude that at
the cost of a slower fitting process rational quadratic manifolds can
significantly outperform quadratic and linear POD manifolds in terms of
reconstruction accuracy for the same number reduced space dimensions.

Figure 8: Total error evolution for the ROMs $\varepsilon_{u}(t)$ and the
ideal linear projection error $\varepsilon_{\text{proj}}(t)$ for reduced
manifolds with constant $r=15$.

Figure 9: Total error evolution for the ROMs $\varepsilon_{u}(t)$ and the
ideal linear projection error $\varepsilon_{\text{proj}}(t)$ for reduced
manifolds with approximately constant $\varepsilon_{xt}^{\text{max}}$.

Figure 10: Spatial profile $\boldsymbol{u}_{r}(t)$ for the ROMs at $t=T$ for
reduced manifolds with constant $r=15$.

Figure 11: Spatial profile $\boldsymbol{u}_{r}(t)$ for the ROMs at $t=T$ for
reduced manifolds with approximately constant $\varepsilon_{xt}^{\text{max}}$.

###  5.2 Shallow water equations

This experiment will mainly focus on the entropy conservation and stability
aspect of our proposed ROM on rational manifolds. We show that our novel
entropy stable ROM satisfies the reduced total entropy estimate (26). To this
end we will carry out experiments with the shallow water equations. We use the
shallow water equations due to their nontrivial entropy function, which will
be defined later, as compared to the Burgers equation. We will perform one
experiment where discontinuities appear in the solution and one where the
solution remains smooth during the time interval of interest. In the smooth
case we can run the FOM and ROMs without entropy dissipation operators. As a
result we can analyse the entropy conservation properties of the ROM. In the
discontinuous case we will analyse the behaviour of the reduced entropy as
compared to the FOM entropy.

In the following we briefly introduce the shallow water equations and the
entropy stable numerical scheme used to obtain the FOM. The shallow water
equations are given by:

$\frac{\partial}{\partial t}\begin{bmatrix}h\\\
hu\end{bmatrix}+\frac{\partial}{\partial x}\begin{bmatrix}hu\\\
hu^{2}+\frac{1}{2}gh^{2}\end{bmatrix}=0,$ (35)

with conserved variables $h:\Omega\times[0,T]\rightarrow\mathbb{R}$ and
$hu:\Omega\times[0,T]\rightarrow\mathbb{R}$ representing the local water
column height and the momentum per unit mass, respectively. We collect the
conserved variables in a vector $\boldsymbol{u}:=\left[h,hu\right]^{T}$. The
constant $g\in\mathbb{R}^{+}$ is the positive real gravitational acceleration,
we assume the equations are normalized such that $g=3$, which gave challenging
test cases for our spatial domain size and initial conditions. As continuous
and discrete total entropy we take the common choice [43]:

$\mathcal{S}[\boldsymbol{u}]=\int_{\Omega}\frac{1}{2}\left(\frac{u_{2}^{2}}{u_{1}}+gu_{1}^{2}\right)dx,$

leading to:

$S_{h}[\boldsymbol{u}_{h}]=\sum_{i}\Delta
x_{i}\frac{h_{i}u_{i}^{2}+gh_{i}^{2}}{2},$

from which the reduced total entropy follows. This choice leads to the
following entropy variables:

$\boldsymbol{\eta}(\boldsymbol{u})=\begin{bmatrix}gu_{1}-\frac{1}{2}\left(\frac{u_{2}}{u_{1}}\right)^{2}\\\
\frac{u_{2}}{u_{1}}\end{bmatrix},$

where $u_{1}=h$ and $u_{2}=hu$, and inverse function:

$\boldsymbol{u}(\boldsymbol{\eta})=\frac{2\eta_{1}+\eta_{2}^{2}}{2g}\begin{bmatrix}1\\\
\eta_{2}\end{bmatrix}.$

An entropy conservative flux is given by [43, 44]:

$\boldsymbol{f}_{i+1/2}=\begin{bmatrix}\overline{h}_{i+1/2}\overline{u}_{i+1/2}\\\
\overline{h}_{i+1/2}\cdot\overline{u}_{i+1/2}^{2}+\frac{1}{2}g\overline{h^{2}}_{i+1/2}\end{bmatrix},$

where $\overline{a}_{i+1/2}=\frac{1}{2}(a_{i+1}+a_{i})$ indicates taking the
average of neighbouring volume based quantities. As entropy dissipation
operator we take the diffusion operators of Roe type (see [44]) with the
eigenvalues and eigenvectors of the flux Jacobian evaluated at the arithmetic
average of neighbouring values. We obtain a second accurate entropy
dissipation operator using the entropy stable total variation diminishing
(TVD) reconstruction based on the minmod limiter (see [44]).

We will discretize (35) for both experiments on a domain
$\Omega=\mathbb{T}([-L,L])$ with $L=1$ using a numerical grid consisting of
300 cells. Discretization in time will be done using the RK4 method with a
time step of size $\Delta t=0.0005$. We will integrate the discontinuous
experiment in time until $T=1$ and the smooth experiment until $T=0.5$.
Solution snapshots are captured every 5 timesteps resulting in $n_{s}=401$ and
$n_{s}=201$ for the discontinuous and smooth experiment, respectively. For the
discontinuous case we will be interested in a dam break problem, this means we
will take as initial condition:

$h_{0}(x)=\begin{cases}1.5&|x|<0.2,\\\ 1&|x|\geq
0.2,\end{cases}\quad(hu)_{0}(x)=0.$

The smooth case will consist of a quiescent water level with a small
perturbation, so that the initial condition is given by:

$h_{0}(x)=1+0.1\cdot\exp\left(-100\cdot x^{2}\right),\quad(hu)_{0}(x)=0.$

The reduced space dimension is taken at $r=15$ for both experiments.

Although it is important that entropy is appropriately conserved or
dissipated, accuracy with respect to the FOM is also required for an effective
ROM. Hence, before we analyse the conservation properties of our proposed ROM,
we compare the FOM solution approximation quality of our entropy stable ROM
and the generic ROM. We will provide space–time plots of both the
discontinuous and smooth experiments. The discontinuous case is given in
Figure 12 and the smooth case is given in Figure 13. Visually, both ROMs
closely resemble the FOM in both cases. Furthermore, it can be seen from the
sharp color gradients that the dam break problem develops shocks. To emphasize
these are difficult cases for linear model reduction approaches we also plot
the normalized singular value decay in Figure 15 and Figure 15 for the dam
break and water height perturbation problems, respectively. For the dam beak
problem the decay is very slow and the water height perturbation decays
moderately slow, indicating slow and moderately slow Kolmogorov $n$-width
decay (1).

Figure 12: $x-t$ plot of solution approximation by the ROMs and FOM for the
dam break problem. Figure 13: $x-t$ plot of solution approximation by the ROMs
and FOM for the water height perturbation problem.

Figure 14: Normalized singular values of the shallow water equations data for
the dam break problem.

Figure 15: Normalized singular values of the shallow water equations data for
the water height perturbation problem.

We will analyse the entropy conservation and stability properties of the
entropy stable (ES-ROM) and generic (RQ-ROM) rational manifold ROMs. To this
end, we define the entropy error:

$\varepsilon_{\mathcal{S}}(t):=\left|S_{h}[\boldsymbol{u}_{h}(t)]-S_{r}[\boldsymbol{a}(t)]\right|,$

giving the absolute instantaneous deviation of the entropy of the ROM from the
entropy of the FOM. Similarly we will define the entropy conservation error:

$\varepsilon_{\mathcal{S}_{0}}(t):=\left|S_{r}[\boldsymbol{a}(0)]-S_{r}[\boldsymbol{a}(t)]\right|,$

which measures the departure from the initial entropy and thus the error in
exact conservation of the entropy in time. Since our models are semi-
discretely entropy stable we have to monitor the instantaneous time rate of
change of reduced total entropy (26) to verify that our proposed theoretical
framework works. Hence, we will analyse the contribution to the total entropy
production (22) of two separate parts of the ROM (25), namely the entropy
conserving part:

$\left(\frac{dS_{r}}{dt}\right)_{\text{cons}}:=-\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}\boldsymbol{f}_{h}^{*}(\boldsymbol{u}(\tilde{\boldsymbol{\eta}}_{r})),$

which should equal zero to machine precision, and the entropy dissipative
part:

$\left(\frac{dS_{r}}{dt}\right)_{\text{diss}}:=\tilde{\boldsymbol{\eta}}_{r}^{T}\Delta_{v}D_{h}(\tilde{\boldsymbol{u}}_{r})\Delta_{i}\tilde{\boldsymbol{\eta}}_{r},$

which should always be negative or zero. Similar quantities can be defined in
an obvious manner for the generic ROM without entropy projection. In the
results given by Figure 17 and Figure 19 we have used symmetric log plots
which are linear around zero so that negative values can also be plotted. This
allows us to see when a ROM is unphysically producing entropy i.e.
$\left(\frac{dS_{r}}{dt}\right)_{\text{cons}},\left(\frac{dS_{r}}{dt}\right)_{\text{diss}}>0$.

Figure 16: The evolution of the reduced total entropy $\mathcal{S}_{r}$ of the
entropy stable and generic ROM compared to the entropy of the FOM and the
evolution of the entropy error $\varepsilon_{\mathcal{S}}$ of the entropy
stable and generic ROM for the dam break problem. Figure 17: The entropy
production $\left(\frac{dS_{r}}{dt}\right)_{\text{cons}}$ of the conservative
part of the ROMs and the entropy production
$\left(\frac{dS_{r}}{dt}\right)_{\text{diss}}$ of the dissipative part of the
ROMs for the dam break problem.

The results of the discontinuous dam break experiment are displayed in Figure
16 and Figure 17. The results confirm that the proposed entropy stable
framework works as expected. This is the case since the entropy production of
the conservative part is zero to machine precision for the entropy stable ROM
while the entropy dissipative part does not change sign and is indeed
negative. The entropy production of the conservative part of the generic ROM
is orders of magnitude larger than that of the entropy stable ROM. An
important point is that large portions in time of the entropy production by
the entropy conservative part are positive (instead of zero). This indicates
physically incorrect behaviour as entropy is being produced. The contribution
of the entropy dissipation operator from the generic ROM is erratic. Moreover,
it is also occasionally positive showing that this part of the ROM is also
sometimes producing physically incorrect results. The temporal evolution of
the entropy error, $\varepsilon_{\mathcal{S}}$, is also given in Figure 16.
The temporal evolution is approximately constant in time for the entropy
stable ROM and $\varepsilon_{\mathcal{S}}$ is small. This indicates that the
evolution of the entropy behaves roughly the same as the FOM and is off mainly
due to an error in representation of the initial condition and of subsequent
FOM solutions. The general behaviour of the entropy error of the generic ROM
is erratic and shows that the entropy of the generic ROM oscillates around the
values predicted by the FOM. This can also be seen in the temporal evolution
of the reduced entropy as in the top panel of Figure 16.

Figure 18: The evolution of the reduced total entropy $\mathcal{S}_{r}$ of the
entropy conserving and generic ROM compared to the entropy of the FOM and the
evolution of the entropy error $\varepsilon_{\mathcal{S}_{0}}$ of the entropy
conserving and generic ROM for the water height perturbation problem. Figure
19: The entropy production $\left(\frac{dS_{r}}{dt}\right)_{\text{cons}}$ of
the conservative part of the ROMs for the water height perturbation problem.

The results of the smooth experiment are shown in Figure 18 and Figure 19. As
there is no entropy dissipation present in the FOM or ROM the entropy should
remain approximately constant (exact conservation is difficult since RK4 is
not an entropy conservative time-integrator [84, 85, 118]). For our entropy
stable ROM this is indeed the case as can be seen from the bottom panel of
Figure 18, the entropy conservation error, $\varepsilon_{\mathcal{S}_{0}}$,
does not exceed $\mathcal{O}(10^{-9})$. The entropy of our entropy stable ROM
stays almost exactly constant. The error in entropy with respect to the FOM is
almost entirely dictated by the initial error
$\varepsilon_{\mathcal{S}}(0)=\left|S_{h}[\boldsymbol{u}_{h}(t)]-S_{r}[\boldsymbol{a}(0)]\right|$.
The generic ROM does not conserve entropy and its entropy conservation error
behaves erratically. This manifests itself in clear deviations from the FOM
entropy which can be observed in the top panel of Figure 18. For completeness
we also plot the entropy production of the conservative part of the spatial
discretization of the ROM $\left(\frac{dS_{r}}{dt}\right)_{\text{cons}}$ in
Figure 19. Again, it can be seen that entropy is conserved up to machine
precision by the spatial discretization of our entropy stable ROM where this
is not the case for the generic ROM. Additionally, the generic ROM produces
entropy during several intervals of the simulation and is therefore not
physically correct. From both experiments, we conclude that our novel entropy
stable ROM ensures physically correct behaviour, whereas this cannot be
assumed for the generic manifold Galerkin ROM.

###  5.3 Compressible Euler equations

The focus of this experiment is on the effect of the entropy projection and
tangent space enrichment on the accuracy of the ROM. We will be interested in
particular in the benefit of tangent space enrichment in the reconstruction
accuracy of the entropy projection. In addition, we will analyse the error
that can be incurred with respect to the FOM by the introduction of an entropy
projection step in the ROM as we propose. A good case to study for this
experiment are the compressible Euler equations. Due to their nontrivial
entropy functional and corresponding entropy variables it is not expected that
without extra measures, like TSE, the entropy projection will be accurate.

A short introduction to the compressible Euler equations now follows. The
compressible Euler equations are given by:

$\frac{\partial}{\partial t}\begin{bmatrix}\rho\\\ \rho u\\\
E\end{bmatrix}+\frac{\partial}{\partial x}\begin{bmatrix}\rho u\\\ \rho
u^{2}+p\\\ (E+p)u\end{bmatrix}=0,$ (36)

where $\rho:\Omega\times[0,T]\rightarrow\mathbb{R}$ is the density, $\rho
u:\Omega\times[0,T]\rightarrow\mathbb{R}$ is the momentum and
$E:\Omega\times[0,T]\rightarrow\mathbb{R}$ is the total energy. We gather
these conserved variables in a vector $\boldsymbol{u}=\left[\rho,\rho
u,E\right]^{T}$. Furthermore, we assume the equations are suitably normalized
so that they are dimensionless. The pressure
$p:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is related to the conserved variables
through an equation of state, representing the thermodynamics at hand:

$p(\boldsymbol{u})=\left(u_{3}-\frac{1}{2}\frac{u_{2}^{2}}{u_{1}}\right)(\gamma-1),$
(37)

where $\gamma\in\mathbb{R}^{+}$ is the specific heat ratio, which we take at
the standard choice $\gamma=1.4$. The thermodynamic quantities, i.e. pressure,
density and total energy are necessarily nonnegative. Throughout the
experiments we will assume our FOM and ROMs respect this condition, assuring
this condition mathematically may be the subject of future work. The entropy
functional of interest will be taken as:

$\mathcal{S}[\boldsymbol{u}]=\int_{\Omega}\frac{-u_{1}\sigma}{\gamma-1}dx,$

where $\sigma:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the specific entropy
defined as a function of the conserved variables like:

$\sigma(\boldsymbol{u})=\ln\left(\frac{p}{u_{1}^{\gamma}}\right),$

where $p$ is evaluated using (37). Different entropies are also possible, see
for instance [54]. In turn, this gives rise to the discrete total entropy
functional:

$S_{h}[\boldsymbol{u}_{h}]=\sum_{i}\Delta
x_{i}\frac{-\rho_{i}\sigma_{i}}{\gamma-1}.$

The associated entropy variables are given by:

$\boldsymbol{\eta}(\boldsymbol{u})=\begin{bmatrix}\frac{\gamma-\sigma}{\gamma-1}-\frac{u_{2}^{2}}{2u_{1}p}\\\
u_{2}/p\\\ -u_{1}/p\end{bmatrix},$

and consequently the inverse of the entropy variables is:

$\boldsymbol{u}(\boldsymbol{\eta})=\exp\left(\frac{\gamma}{1-\gamma}-\left[\eta_{1}-\frac{1}{2}\frac{\eta_{2}^{2}}{\eta_{3}}\right]\right)\begin{bmatrix}(-\eta_{3})^{\frac{1}{1-\gamma}}\\\
-\eta_{2}(-\eta_{3})^{\frac{\gamma}{1-\gamma}}\\\
\left(\frac{1}{2}\eta_{2}^{2}-\frac{\eta_{3}}{\gamma-1}\right)\cdot(-\eta_{3})^{\frac{2\gamma-1}{1-\gamma}}\end{bmatrix}.$

Considering $-\eta_{3}=\rho/p$ is generally exponentiated to some noninteger
power we see the importance of positivity of the thermodynamic variables. As
in [44], we will use the entropy conserving numerical flux by Ismail and Roe
[65] for which we define the following variables:

$\boldsymbol{z}=\begin{bmatrix}z^{1}\\\ z^{2}\\\
z^{3}\end{bmatrix}=\sqrt{\frac{\rho}{p}}\begin{bmatrix}1\\\ u\\\
p\end{bmatrix},$

finally an entropy conservative flux is given by
$\boldsymbol{f}_{i+1/2}=\begin{bmatrix}\boldsymbol{F}_{i+1/2}^{1}&\boldsymbol{F}_{i+1/2}^{2}&\boldsymbol{F}_{i+1/2}^{3}\end{bmatrix}^{T}$:

$\displaystyle\boldsymbol{F}_{i+1/2}^{1}$
$\displaystyle=\overline{z^{2}}_{i+1/2}\cdot(z^{3})_{i+1/2}^{\ln},$
$\displaystyle\boldsymbol{F}_{i+1/2}^{2}$
$\displaystyle=\frac{\overline{z^{3}}_{i+1/2}}{\overline{z^{1}}_{i+1/2}}+\frac{\overline{z^{2}}_{i+1/2}}{\overline{z^{1}}_{i+1/2}}\boldsymbol{F}_{i+1/2}^{1},$
$\displaystyle\boldsymbol{F}_{i+1/2}^{3}$
$\displaystyle=\frac{1}{2}\frac{\overline{z^{2}}_{i+1/2}}{\overline{z^{1}}_{i+1/2}}\left(\frac{\gamma+1}{\gamma-1}\frac{(z^{3})_{i+1/2}^{\ln}}{(z^{1})_{i+1/2}^{\ln}}+\boldsymbol{F}_{i+1/2}^{2}\right),$

where $a^{\ln}$ denotes the logarithmic mean, which is defined as:

$a_{i+1/2}^{\ln}:=\frac{a_{i+1}-a_{i}}{\ln{a_{i+1}}-\ln{a_{i}}}.$

Computation of the logarithmic mean is generally not numerically stable when
$a_{i+1}\approx a_{i}$, but a popular algorithm which we will use to deal with
this is also given in [65]. There is an abundance of alternative entropy
conservative numerical fluxes that can be used [101, 29, 57, 117, 75, 120]
some of which also conserve kinetic energy in the sense of Jameson [67]. As an
entropy dissipation operator we take the Roe-type diffusion operator [44]
where the eigenvalues and vectors of the flux Jacobian are evaluated at the
arithmetic average of the neighbouring conserved variables. As with the
shallow water equations, we obtain a second accurate entropy dissipation
operator using the entropy stable total variation diminishing (TVD)
reconstruction based on the minmod limiter [44].

For the experiment we will consider a periodic modification of the famous
Sod’s shock tube [114], which avoids the need to implement entropy stable
boundary conditions. We will discretize (36) on a domain
$\Omega=\mathbb{T}([0,L])$ with $L=1$ on a numerical grid of $250$ cells. The
number of cells is relatively small to facilitate a relatively short manifold
learning process. Integration of the ROM in time will be carried out using the
RK4 time integrator with a time step size $\Delta t=0.0001$. We will integrate
in time until $T=0.5$ (beyond the typical time used for this experiment),
resulting in interesting shock-rarefaction interactions. Again, we will
capture snapshots after every 5 timesteps so that we have $n_{s}=1001$. Our
periodic modification of Sod’s shock tube experiment has an initial condition
given by:

$\rho_{0}(x)=\begin{cases}1&0.25<x<0.75\\\
0.125&\text{elsewhere}\end{cases}\quad u_{0}(x)=0,\quad
p_{0}(x)=\begin{cases}1&0.25<x<0.75\\\ 0.1&\text{elsewhere}\end{cases},$

where the conserved variables $(\rho,\rho u,E)$ are calculated from these
primitive variables using the equation of state (37) and the definition of
momentum. We take a reduced space dimension $r=15$.

We will be primarily interested in the entropy projection and tangent space
enrichment during this experiment, but for completeness we also plot the ROM
approximations to the FOM solution and the singular values. The singular
values are displayed in Figure 20 and the ROM approximations are shown using
$x-t$ plots in Figure 21. A relatively slow decay of singular values can be
observed in Figure 20, hence linear model reduction methods are likely to not
perform well for this problem. The FOM solution is approximated well by our
novel entropy stable manifold Galerkin ROM. The solution approximations of the
ROM as displayed in Figure 21 are nearly identical to the FOM.

Figure 20: Normalized singular values of the compressible Euler equations data
for the modified Sod’s shock tube problem. Figure 21: $x-t$ plot of solution
approximation by the entropy stable ROM compared to the FOM for the modified
Sod’s shock tube problem.

We are interested in the accuracy of the entropy projection with and without
tangent space enrichment. Accordingly, we introduce a metric to measure this
accuracy. Since we are not only interested in comparing errors, but also to
get an idea of the absolute size of the error we specifically introduce the
relative entropy projection error:

$\varepsilon_{\Pi}(t):=\frac{||\boldsymbol{u}_{r}(t)-\boldsymbol{u}(\Pi_{T\mathcal{M}}\boldsymbol{\eta}_{r}(t))||_{\Omega_{h}}}{||\boldsymbol{u}_{r}(t)||_{\Omega_{h}}},$

measuring not only how far the entropy projection is from the identity mapping
as with the entropy projection error $\varepsilon_{s}$ of subsection 3.3, but
also the size of the error $\varepsilon_{s}$ relative to the approximated
value $\boldsymbol{u}_{r}$. We will plot this value for two ROM simulations of
the compressible Euler equations with an entropy projection, where one has an
enriched tangent space and the other not. The results are provided in Figure
22. The entropy projection error with TSE can be seen to be a very small
fraction of the norm $||\boldsymbol{u}_{r}(t)||_{\Omega_{h}}$, indicating
minimal impact on the accuracy of the ROM given it is well-conditioned. In
contrast, the ROM without TSE instantly produces NaN values and could
therefore not be included in Figure 22. To obtain a further comparison we plot
the spatial profiles of the projected entropy variables at two selected
moments $t_{p}\in\mathbb{R}^{+}$ in time, namely $t_{p}\in\\{0.1,0.5\\}$. To
have a meaningful comparison, i.e. one where we are not projecting NaN values
to start with, we calculate the entropy variables from the stable ROM with
enriched tangent space. Furthermore, we use the the generalized coordinates
$\boldsymbol{a}_{p}=\Phi^{T}X_{p}$ to compute the tangent space basis for the
ROM without TSE. The results are shown in Figure 23 and Figure 24. Very poor
reconstruction of entropy variables can be observed for the ROM without TSE,
whereas with TSE the reconstruction is accurate at both moments. From Figure
23 and Figure 24 the NaN values in Figure 22 can be explained by the
projection of the entropy variables taking unphysical values (positive
$\eta_{3}$). From this we conclude that tangent space enrichment or any other
manner of assuring the accuracy of the entropy projection is vital for a
properly functioning ROM when using an entropy projection.

Figure 22: The relative entropy projection error for a ROM with TSE. The ROM
without TSE resulted in all NaN values and is therefore not included in the
figure. Figure 23: Entropy variable $\boldsymbol{\eta}_{r}$ and conserved
variable $\boldsymbol{u}_{r}$ approximation by entropy projection with and
without TSE at $t_{p}=0.1$. Figure 24: Entropy variable
$\boldsymbol{\eta}_{r}$ and conserved variable $\boldsymbol{u}_{r}$
approximation by entropy projection with and without TSE at $t_{p}=0.5$.

In applying the tangent space enrichment framework, we rely on the artificial
TSE coordinate $\alpha$ staying small ($\alpha\ll 1$) during simulations. If
this is not the case we cannot assure that the reduced space can accurately
represent the solution nor that the local tangent space can accurately
represent the FOM dynamical system $\frac{d\boldsymbol{u}_{h}}{dt}$ at that
point. The reason for this being that the enriched manifold parameterization
$\hat{\boldsymbol{\varphi}}$ and its Jacobian matrix $\hat{\boldsymbol{J}}$
are no longer close to the original parameterization $\boldsymbol{\varphi}$
and Jacobian $\boldsymbol{J}$ which are accurate by assumption. To verify this
is indeed not the case we will monitor the value of $\alpha$ throughout a
simulation of the compressible Euler equations. The results are shown in
Figure 25. We have also plotted the error $\varepsilon_{u}$ of the ROM with
entropy projection and tangent space enrichment and of a generic manifold
Galerkin ROM for reference in Figure 26. It can be seen in Figure 25 that the
value of $\alpha$ remains small around $\mathcal{O}(10^{-5})$. Consequently,
the original manifold parameterization $\boldsymbol{\varphi}$ and Jacobian
$\boldsymbol{J}$ are well-approximated by their enriched counterparts
$\hat{\boldsymbol{\varphi}}$ and $\hat{\boldsymbol{J}}$. It can be seen in
Figure 26 that this is, in fact, the case, as the errors of the ROM differ by
at most $\mathcal{O}(10^{-5})$ and evolve very similarly. Hence, we conclude
that the entropy projection step, introduced to obtain an entropy stable
framework, is not detrimental to the accuracy of the ROM, provided that an
enriched tangent space is used.

Figure 25: Value of the artificial TSE coordinate $\alpha$ throughout a
simulation of the modified Sod’s shock tube experiment. Figure 26: Value of
the ROM error $\varepsilon_{u}$ throughout a simulation of the modified Sod’s
shock tube experiment.

##  6 Conclusion

In this article we have proposed a method to construct nonlinear manifold
Galerkin reduced order models (ROMs) in such a way that the important total
entropy functional of the ROM approximation is appropriately conserved or
dissipated. This is a crucial concept in obtaining stable and physically
admissible ROM solutions. In particular, we have focused on systems of one-
dimensional nonlinear conservation laws. Correct semi-discrete entropy
estimates upon orthogonal projection were obtained for these systems by
evaluating the projected system not at the current state, but at its entropy
projection. This was proposed earlier for linear ROMs and extended in this
work to nonlinear manifold ROMs.

The entropy projection of the state is obtained by transforming conservative
variables to entropy variables, consequently projecting these on the tangent
space of the reduced manifold and finally transforming back to conserved
variables. To assure accuracy, it is important that the entropy projection is
as close as possible to the identity mapping. This is generally not the case
for general nonlinear reduced spaces and hence we have proposed the method of
tangent space enrichment (TSE). With TSE the manifold is lifted along an
additional dimension parameterized by a new coordinate. This coordinate
direction is constructed to linearly extend in the direction of the local
entropy variables, so that the tangent space spans the entropy variable at
least approximately given the absolute value of the TSE coordinate.
Accordingly, the entropy projection error will stay small.

We have tested our proposed framework on several nonlinear conservation laws
from fluid dynamics. We verified that the entropy estimates are satisfied
(semi-discretely): the projection of entropy-conserving flux differences
produces no total entropy and the projection of entropy dissipative terms
dissipates total entropy. This is in contrast to the generic manifold Galerkin
framework which leads to production of entropy in our numerical experiments,
which is physically incorrect. We have also shown that the introduction of the
artificial TSE coordinate is vital for the accuracy of the entropy projection
and leads to minimal decreases in accuracy.

We have also for the first time generalized the recently proposed quadratic
manifolds to rational quadratic manifolds. We have suggested a framework to
find the coefficients of the rational quadratic manifolds based on a nonlinear
curve fitting approach. We have also formulated the rational quadratic
polynomials such that they do no not have real poles. This was achieved
through semi-definite constraints, avoiding division by zero for any point in
the reduced space. Numerical experiments on the inviscid Burgers equation have
shown the increased performance of these rational quadratic manifold
parameterizations as compared to existing quadratic manifold parameterizations
and linear methods.

In future work, two challenges need to be tackled to make the approach
computationally efficient: (i) we need a faster way to fit the rational
quadratic manifolds and (ii) we need an entropy-stable hyperreduction
approach. The former can possibly be achieved through linearization and
iterative techniques combined with better choices of generalized coordinates
[110], whereas the latter could be achieved by adapting the constrained
optimization formulation that we proposed for energy-conserving systems in
[73]. In addition, the framework would benefit from extension with an entropy-
stable time integration method and entropy-stable treatment of boundary
conditions.

## CRediT authorship contribution statement

R.B. Klein: conceptualization, methodology, software, validation, formal
analysis, investigation, writing - original draft
B. Sanderse: writing - review & editing, supervision
P. Costa: writing - review & editing, supervision
R. Pecnik: writing - review & editing, supervision
R.A.W.M. Henkes: writing - review & editing, supervision, project
administration, funding acquisition

## Declaration of competing interests

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.

## Data availability

Data will be made available on request.

## Acknowledgements

R.B. Klein gratefully acknowledges the funding for this project obtained from
Delft University of Technology.

## References

  * [1] S.E. Ahmed and O. San.  Breaking the Kolmogorov barrier in model reduction of fluid flows.  Fluids, 5(1), 2020.
  * [2] D. Amsallem and C. Farhat.  Stabilization of projection-based reduced-order models.  International Journal for Numerical Methods in Engineering, 91(4):358–377, 2012.
  * [3] A.P. Austin, M. Krishnamoorthy, S. Leyffer, S. Mrenna, J. Müller, and H. Schulz.  Practical algorithms for multivariate rational approximation.  Computer Physics Communications, 261:107663, 2021.
  * [4] S. Axler.  Linear Algebra Done Right.  Undergraduate Texts in Mathematics. Springer International Publishing, 2014.
  * [5] J. Axås and G. Haller.  Model reduction for nonlinearizable dynamics via delay-embedded spectral submanifolds, 2023.
  * [6] M. Bachmayr and A. Cohen.  Kolmogorov widths and low-rank approximations of parametric ellipitic PDEs.  Mathematics of Computation, 86(304):pp. 701–724, 2017.
  * [7] M. Balajewicz and E.H. Dowell.  Stabilization of projection-based reduced order models of the Navier-Stokes equations.  Nonlinear Dynamics, 70:1619–1632, 10 2012.
  * [8] J. Barnett and C. Farhat.  Quadratic approximation manifold for mitigating the Kolmogorov barrier in nonlinear projection-based model order reduction.  Journal of Computational Physics, 464:111348, 2022.
  * [9] J. Barnett, C. Farhat, and Y. Maday.  Neural-network-augmented projection-based model order reduction for mitigating the Kolmogorov barrier to reducibility.  Journal of Computational Physics, 492:112420, 2023.
  * [10] M.F. Barone, I. Kalashnikova, D.J. Segalman, and H.K. Thornquist.  Stable Galerkin reduced order models for linearized compressible flow.  Journal of Computational Physics, 228(6):1932–1946, 2009.
  * [11] P. Benner, S. Gugercin, and K. Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Review, 57(4):483–531, 2015.
  * [12] S. Bianchini and A. Bressan.  Vanishing viscosity solutions of nonlinear hyperbolic systems.  Annals of Mathematics, 161(1):223–342, 2005.
  * [13] T. Blickhan.  A registration method for reduced basis problems using linear optimal transport, 2023.
  * [14] J. Borggaard, A. Hay, and D. Pelletier.  Interval-based reduced-order models for unsteady fluid flow.  International Journal of Numerical Analysis and Modeling, 4:353–367, 01 2007.
  * [15] J. Bradbury, R. Frostig, P. Hawkins, M.J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang.  JAX: composable transformations of Python+NumPy programs, 2018.
  * [16] B. Brantner, G. de Romemont, M. Kraus, and Z. Li.  Structure-preserving transformers for learning parametrized Hamiltonian systems, 2023.
  * [17] P. Buchfink, S. Glas, and B. Haasdonk.  Approximation bounds for model reduction on polynomially mapped manifolds, 2023.
  * [18] P. Buchfink, S. Glas, and B. Haasdonk.  Symplectic model reduction of Hamiltonian systems on nonlinear manifolds and approximation with weakly symplectic autoencoder.  SIAM Journal on Scientific Computing, 45:A289–A311, 03 2023.
  * [19] P. Buchfink, S. Glas, B. Haasdonk, and B. Unger.  Model reduction on manifolds: A differential geometric framework, 2024.
  * [20] G. Buza.  Spectral submanifolds of the Navier-Stokes equations, 2023.
  * [21] G. Buza, S. Jain, and G. Haller.  Using spectral submanifolds for optimal mode selection in nonlinear model reduction.  Proceedings of the Royal Society A, 477(2246):20200725, 2021.
  * [22] K.T. Carlberg.  Adaptive h-refinement for reduced-order models.  International Journal for Numerical Methods in Engineering, 102(5):1192–1210, 2015.
  * [23] K.T. Carlberg, M. Barone, and H. Antil.  Galerkin v. least-squares Petrov–Galerkin projection in nonlinear model reduction.  Journal of Computational Physics, 330:693–734, 2017.
  * [24] K.T. Carlberg, C. Bou-Mosleh, and C. Farhat.  Efficient non-linear model reduction via a least-squares Petrov–Galerkin projection and compressive tensor approximations.  International Journal for Numerical Methods in Engineering, 86(2):155–181, 2011.
  * [25] M.H. Carpenter, T.C. Fisher, E.J. Nielsen, and S.H. Frankel.  Entropy stable spectral collocation schemes for the Navier-Stokes equations: Discontinuous interfaces.  SIAM Journal on Scientific Computing, 36(5):B835–B867, 2014.
  * [26] J. Chan.  On discretely entropy conservative and entropy stable discontinuous Galerkin methods.  Journal of Computational Physics, 362:346–374, 2018.
  * [27] J. Chan.  Entropy stable reduced order modeling of nonlinear conservation laws.  Journal of Computational Physics, 423:109789, 2020.
  * [28] J. Chan, H. Ranocha, A.M. Rueda-Ramírez, G. Gassner, and T. Warburton.  On the entropy projection and the robustness of high order entropy stable discontinuous Galerkin schemes for under-resolved flows.  Frontiers in Physics, 10, 2022.
  * [29] P. Chandrashekar.  Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations.  Communications in Computational Physics, 14(5):1252–1286, 2013.
  * [30] S. Chaturantabut.  Temporal localized nonlinear model reduction with a priori error estimate.  Applied Numerical Mathematics, 119:225–238, 2017.
  * [31] J. Cocola, J. Tencer, F. Rizzi, E. Parish, and P. Blonigan.  Hyper-reduced autoencoders for efficient and accurate nonlinear model reductions, 2023.
  * [32] A. Cohen and R. DeVore.  Kolmogorov widths under holomorphic mappings.  IMA Journal of Numerical Analysis, 36(1):1–12, 03 2015.
  * [33] D.M. Copeland, S.W. Cheung, K. Huynh, and Y. Choi.  Reduced order models for Lagrangian hydrodynamics.  Computer Methods in Applied Mechanics and Engineering, 388:114259, 2022.
  * [34] C.M. Dafermos.  Hyperbolic Conservation Laws in Continuum Physics.  Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 2006.
  * [35] M. D’Elia, M. Gunzburger, and G. Rozza.  Quantification of Uncertainty: Improving Efficiency and Technology QUIET selected contributions: QUIET selected contributions.  01 2020.
  * [36] F. Della Santa and S. Pieraccini.  Discontinuous neural networks and discontinuity learning.  Journal of Computational and Applied Mathematics, 419:114678, 2023.
  * [37] A.N. Diaz, Y. Choi, and M. Heinkenschloss.  A fast and accurate domain-decomposition nonlinear manifold reduced order model, 2024.
  * [38] V. Ehrlacher, D. Lombardi, O. Mula, and F. Vialard.  Nonlinear model reduction on metric spaces. Application to one-dimensional conservative PDEs in Wasserstein spaces.  ESAIM: Mathematical Modelling and Numerical Analysis, 54, 03 2020.
  * [39] P.A. Etter and K.T. Carlberg.  Online adaptive basis refinement and compression for reduced-order models via vector-space sieving.  Computer Methods in Applied Mechanics and Engineering, 364:112931, 2020.
  * [40] L.C. Evans.  Partial Differential Equations.  Graduate studies in mathematics. American Mathematical Society, 2010.
  * [41] L. Fick, Y. Maday, A.T. Patera, and T. Taddei.  A stabilized POD model for turbulent flows over a range of Reynolds numbers: Optimal parameter sampling and constrained projection.  Journal of Computational Physics, 371:214–243, 2018.
  * [42] T.C. Fisher and M.H. Carpenter.  High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains.  Journal of Computational Physics, 252:518–557, 2013.
  * [43] U. Fjordholm, S. Mishra, and E. Tadmor.  Energy Preserving and Energy Stable Schemes for the Shallow Water Equations, page 93–139.  London Mathematical Society Lecture Note Series. Cambridge University Press, 2009.
  * [44] U.S. Fjordholm, S. Mishra, and E. Tadmor.  Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws.  SIAM Journal on Numerical Analysis, 50(2):544–573, 2012.
  * [45] H. Gao and M.J. Zahr.  An adaptive model reduction method leveraging locally supported basis functions, 2023.
  * [46] G.J. Gassner, A.R. Winters, and D.A. Kopriva.  Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations.  Journal of Computational Physics, 327:39–66, 2016.
  * [47] R. Geelen, S. Wright, and K. Willcox.  Operator inference for non-intrusive model reduction with quadratic manifolds.  Computer Methods in Applied Mechanics and Engineering, 403:115717, 2023.
  * [48] J. Gerbeau and D. Lombardi.  Approximated Lax pairs for the reduced order integration of nonlinear evolution equations.  Journal of Computational Physics, 265:246–269, 2014.
  * [49] C. Greif and K. Urban.  Decay of the Kolmogorov n-width for wave problems.  Applied Mathematics Letters, 96:216–222, 2019.
  * [50] D.F. Griffiths and D.J. Higham.  Numerical Methods for Ordinary Differential Equations: Initial Value Problems.  Springer Undergraduate Mathematics Series. Springer London, 2010.
  * [51] S. Grimberg, C. Farhat, R. Tezaur, and C. Bou-Mosleh.  Mesh sampling and weighting for the hyperreduction of nonlinear Petrov–Galerkin reduced-order models with local reduced-order bases.  International Journal for Numerical Methods in Engineering, 122(7):1846–1874, 2021.
  * [52] A. Gruber and I. Tezaur.  Canonical and noncanonical Hamiltonian operator inference.  Computer Methods in Applied Mechanics and Engineering, 416:116334, 2023.
  * [53] B. Gustavsen.  Improving the pole relocating properties of vector fitting.  IEEE Transactions on Power Delivery, 21(3):1587–1592, 2006.
  * [54] A. Harten.  On the symmetric form of systems of conservation laws with entropy.  Journal of Computational Physics, 49(1):151–164, 1983.
  * [55] D. Hartmann, M. Herz, and U. Wever.  Model Order Reduction a Key Technology for Digital Twins, pages 167–179.  04 2018.
  * [56] J.S. Hesthaven, C. Pagliantini, and N. Ripamonti.  Rank-adaptive structure-preserving model order reduction of Hamiltonian systems.  ESAIM: Mathematical Modelling and Numerical Analysis, 56(2):617–650, 2022.
  * [57] J.E. Hicken and J. Crean.  A family of entropy-conservative flux functions for the Euler equations, 2019.
  * [58] L.R. Hofer, M. Krstajić, and R.P. Smith.  JAXFit: Trust region method for nonlinear least-squares curve fitting on the GPU, 2022.
  * [59] J.M. Hokanson.  Multivariate rational approximation using a stabilized Sanathanan-Koerner iteration, 2020.
  * [60] J.M. Hokanson and C.C. Magruder.  Least squares rational approximation, 2018.
  * [61] P. Holmes, J.L. Lumley, G. Berkooz, and C.W. Rowley.  Turbulence, Coherent Structures, Dynamical Systems and Symmetry.  Cambridge Monographs on Mechanics. Cambridge University Press, 2012.
  * [62] T.J.R. Hughes and M. Mallet.  A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective-diffusive systems.  Computer Methods in Applied Mechanics and Engineering, 58(3):305–328, 1986.
  * [63] W.L. IJzerman.  Signal Representation and Modeling of Spatial Structures in Fluids.  PhD thesis, University of Twente, April 2000.
  * [64] A. Iollo and D. Lombardi.  Advection modes by optimal mass transfer.  Phys. Rev. E, 89:022923, Feb 2014.
  * [65] F. Ismail and P.L. Roe.  Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks.  Journal of Computational Physics, 228(15):5410–5436, 2009.
  * [66] S. Jain, P. Tiso, J.B. Rutzmoser, and D.J. Rixen.  A quadratic manifold for model order reduction of nonlinear structural dynamics.  Computers & Structures, 188:80–94, 2017.
  * [67] A. Jameson.  Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes.  Journal of Scientific Computing, 34:188–208, 02 2008.
  * [68] I. Kalashnikova and S. Arunajatesan.  A stable Galerkin reduced order model (ROM) for compressible flow.  Blucher Mechanical Engineering Proceedings, 1(1):1399 – 1423, 2014.
  * [69] I. Kalashnikova and M. Barone.  Stable and Efficient Galerkin Reduced Order Models for Non-Linear Fluid Flow.
  * [70] D.I. Ketcheson.  Relaxation Runge-Kutta methods: Conservation and stability for inner-product norms.  SIAM Journal on Numerical Analysis, 57(6):2850–2870, 2019.
  * [71] M. Khamlich, F. Pichi, and G. Rozza.  Optimal transport-inspired deep learning framework for slow-decaying problems: Exploiting Sinkhorn loss and Wasserstein kernel, 2023.
  * [72] Y. Kim, Y. Choi, D. Widemann, and T. Zohdi.  A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder.  Journal of Computational Physics, 451:110841, 2022.
  * [73] R.B. Klein and B. Sanderse.  Energy-conserving hyper-reduction and temporal localization for reduced order models of the incompressible Navier-Stokes equations.  Journal of Computational Physics, 499:112697, 2024.
  * [74] S.N. Kružkov.  First order quasilinear equations in several independent variables.  Mathematics of the USSR-Sbornik, 10(2):217, feb 1970.
  * [75] Y. Kuya and S. Kawai.  High-order accurate kinetic-energy and entropy preserving (KEEP) schemes on curvilinear grids.  Journal of Computational Physics, 442:110482, 2021.
  * [76] T. Lassila, A. Manzoni, A. Quarteroni, and G. Rozza.  Model Order Reduction in Fluid Dynamics: Challenges and Perspectives, pages 235–273.  Springer International Publishing, Cham, 2014.
  * [77] P.D. Lax.  Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves.  Society for Industrial and Applied Mathematics, 1973.