# Structure-Preserving Hyper-Reduction and Temporal Localization for Reduced
Order Models of Incompressible Flows

R.B. Klein<EMAIL_ADDRESS>B. Sanderse<EMAIL_ADDRESS>Centrum Wiskunde &
Informatica,Science Park 123, Amsterdam, The Netherlands  Delft University of
Technology, Process and Energy, Leeghwaterstraat 39, Delft, The Netherlands

###### Abstract

A novel hyper-reduction method is proposed that conserves kinetic energy and
momentum for reduced order models of the incompressible Navier-Stokes
equations. The main advantage of conservation of kinetic energy is that it
endows the hyper-reduced order model (hROM) with a nonlinear stability
property. The new method poses the discrete empirical interpolation method
(DEIM) as a minimization problem and subsequently imposes constraints to
conserve kinetic energy. Two methods are proposed to improve the robustness of
the new method against error accumulation: oversampling and Mahalanobis
regularization. Mahalanobis regularization has the benefit of not requiring
additional measurement points. Furthermore, a novel method is proposed to
perform structure-preserving temporal localization with the principle interval
decomposition: new interface conditions are derived such that energy and
momentum are conserved for a full time-integration instead of only during
separate intervals. The performance of the new structure-preserving hyper-
reduction methods and the structure-preserving temporal localization method is
analysed using two convection-dominated test cases; a shear-layer roll-up and
two-dimensional homogeneous isotropic turbulence. It is found that both
Mahalanobis regularization and oversampling allow hyper-reduction of these
test cases. Moreover, the Mahalanobis regularization provides comparable
robustness while being more efficient than oversampling.

###### keywords:

Discrete empirical interpolation method , Structure preservation ,
Incompressible Navier-Stokes equations , Reduced order models , Mahalanobis
regularization , Temporal localization

##  1 Introduction

Computational fluid dynamics (CFD) has become an integral part of many modern
engineering applications. The increase in computational power in recent
decades has allowed engineers to model increasingly larger fluid dynamical
systems. However, many modern applications are of a multi-query or real-time
nature e.g. design optimization [68] and uncertainty quantification [28] or
real-time control [67, 60] and digital twin technology [37]. These
applications still pose prohibitively large computational costs. Reduced order
models (ROMs) have been proposed as a solution to this problem. A ROM is a
type of surrogate model that approximates the high-dimensional full scale
model, or full order model (FOM), in a low-dimensional way by finding
approximate formulations of involved quantities or operators. The low-
dimensionality of the model consequently makes the ROM significantly cheaper
to evaluate than the FOM.

Traditionally in the CFD community, these ROMs have been constructed by
projecting the fluid dynamics equations of interest onto low-dimensional
linear spaces obtained from the proper orthogonal decomposition (POD)
algorithm [40, 79], using either Galerkin [73, 67, 81, 43, 53] or Petrov-
Galerkin [17, 35, 16] projection. More recently, alternatives have also been
explored. ROMs have been constructed without availability of a FOM by
inferring the ROM from data using operator inference methods [64, 78, 11, 52].
Machine learning methods like convolutional autoencoders have been used in a
projective sense [47, 66] and also inference methods [58, 51] have been
applied to obtain nonlinear low-dimensional approximations. Other nonlinear
dimensionality reduction methods like quadratic manifolds [9, 32] and
diffusion maps [80] have also been leveraged. However, the traditional POD-
Galerkin methods remain powerful and this article will primarily deal with
these methods and their natural extensions like the principal interval
decomposition (PID) [13, 42].

Nevertheless, the traditional methods possess several limitations. For
turbulent (and convection-dominated) flows that are of engineering interest it
is well-known that linear, projection-based ROMs suffer from stability and
accuracy issues [73, 31, 7, 5, 53]. Efforts have been made to solve this
issue, an overview is provided by [31, 73]. A promising solution is structure-
preservation, this entails constructing ROMs such that the underlying physics
of fluid flows are respected. Especially conservation of kinetic energy is an
important physical principle to uphold in a ROM with regards to stability as
it bounds the norm of the solution [73, 2, 20, 82].

Current attempts at constructing structure-preserving ROMs of fluid flow can
be classified into four categories. The first category is constrained
optimization projection [18, 77, 36, 12, 41]; here model reduction is cast
into an optimization problem and constrained to preserve structure. A second
category is formed by symmetry-preservation [73, 2, 20, 82] where symmetries
of the continuous analogues of ROM operators have been preserved resulting in
the conservation of energy or entropy and associated stability. In [34, 22,
65, 39, 1] Hamiltonian physics are preserved in a low-dimensional setting,
forming a third category. Finally, in [56, 85, 48, 57, 8, 4] physics-informed
data-driven approaches using inference and machine learning methods have been
adopted establishing a fourth category.

The idea of structure preservation is elegant but does not resolve directly
the well-known issue that reducing nonlinear models using POD requires
intermediate lifting of the reduced representation to the high-dimensional
(FOM) spaces. Thus, although the ROM is low-dimensional, the computational
effort to evaluate it is still high-dimensional, defeating its purpose.
Methods to overcome this problem are referred to as hyper-reduction methods
[15]. For sufficiently simple equations exact hyper-reduction methods exist
that eliminate computational dependence on the FOM dimensions [73, 3]. This
can be done when the underlying nonlinearity is of polynomial nature. Yet,
these methods can become prohibitively expensive for larger ROMs as noted in
[73]. Hence, approximate hyper-reduction methods may be considered since these
generally have better scaling properties. Examples of such approximate hyper-
reduction methods are the discrete empirical interpolation method (DEIM) [23,
25, 24, 10], Gauss-Newton with approximated tensors (GNAT) [19, 17] and
energy-conserving sampling and weighting (ECSW) [35, 30, 29].

Many existing hyper-reduction methods do not preserve the structure of the
operators to which they are applied and thus stability can be lost. The field
of model reduction of Hamiltonian systems has offered some solutions for this;
[22] proposes a structure-preserving DEIM variant for systems with Hamiltonian
functionals with non-quadratic terms; [83] improves this DEIM variant and [55]
proposes a structure-preserving DEIM variant for nonlinear Hamiltonian
operators that preserves skew-symmetry. However, the incompressible Navier-
Stokes equations considered in this work have a quadratic Hamiltonian
functional in the inviscid case [61, 6] (kinetic energy) making the first two
methods ([22, 83]) inapplicable. Furthermore, the approximate method proposed
in [55] scales computationally the same as the prohibitively expensive exact
method used in [73].

In this article we propose a novel structure-preserving DEIM that is capable
to robustly, accurately and efficiently deal with convection-dominated flows
described by the incompressible Navier-Stokes equations. The main idea will be
to relax the DEIM interpolation requirement and use the resulting available
degree of freedom to enforce energy conservation. Furthermore, we will propose
two methods to further enhance the robustness of the structure-preserving
DEIM. Namely, oversampling [63, 69] and the new approach of regularization
using the Mahalanobis distance (which, to our knowledge, has not yet been
applied to the DEIM). Finally, we will extend the feasibility of the hyper-
reduced order model (hROM) with a newly proposed structure-preserving
localization method based on the PID.

This article will be organized as follows. First, the governing equations will
be introduced in section 2, together with a structure-preserving FOM and a
structure-preserving ROM. Following this, the structure-preserving hyper-
reduction method will be proposed in section 3. In section 4 the new
structure-preserving interval decomposition approach for convection-dominated
flows will be discussed. Finally, in section 5 the performance of the proposed
structure-preserving hyper-reduction method and structure-preserving temporal
localization method will be analysed using two convection-dominated test
cases.

##  2 Preliminaries: Navier-Stokes equations, FOM and ROM, energy equations

###  2.1 Incompressible Navier-Stokes equations and kinetic energy
conservation

In this work we consider the incompressible Navier-Stokes equations:

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+\nabla\cdot\left(\boldsymbol{u}\otimes\boldsymbol{u}\right)=-\nabla
p+\nu\Delta\boldsymbol{u}+\boldsymbol{f}^{b},$ (1)
$\displaystyle\nabla\cdot\boldsymbol{u}=0,$ (2)

on a domain $\Omega\subset\mathbb{R}^{d}$ with periodic boundary conditions.
Several quantities related to the flow will also be considered, namely kinetic
energy $K(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$, defined as:

$K(t)=\frac{1}{2}\left|\left|\boldsymbol{u}\right|\right|_{L^{2}}^{2},$ (3)

and momentum $\boldsymbol{P}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{d}$,
defined as:

$\boldsymbol{P}(t)=\int_{\Omega}\boldsymbol{u}\hbox{ }d\Omega.$ (4)

Denoting
$\mathcal{C}(\boldsymbol{u},\boldsymbol{u}):=\nabla\cdot\left(\boldsymbol{u}\otimes\boldsymbol{u}\right)$,
$\mathcal{D}\boldsymbol{u}:=\Delta\boldsymbol{u}$, $\mathcal{G}p:=\nabla p$
and $\mathcal{M}\boldsymbol{u}:=\nabla\cdot\boldsymbol{u}$, an evolution
equation of the total kinetic energy can be found by differentiating (3) and
substituting (1):

$\displaystyle\frac{dK}{dt}$
$\displaystyle=-\left<\boldsymbol{u},\mathcal{C}(\boldsymbol{u},\boldsymbol{u})\right>_{L^{2}}-\left<\boldsymbol{u},\mathcal{G}p\right>_{L^{2}}+\nu\left<\boldsymbol{u},\mathcal{D}\boldsymbol{u}\right>_{L^{2}}.$

This evolution equation can be further simplified when the following is
considered. Since the convection operator is skew-adjoint given condition (2),
it holds that
$\left<\boldsymbol{u},\mathcal{C}(\boldsymbol{u},\boldsymbol{u})\right>_{L^{2}}=-\left<\mathcal{C}(\boldsymbol{u},\boldsymbol{u}),\boldsymbol{u}\right>_{L^{2}}=-\left<\boldsymbol{u},\mathcal{C}(\boldsymbol{u},\boldsymbol{u})\right>_{L^{2}}$,
meaning that
$\left<\boldsymbol{u},\mathcal{C}(\boldsymbol{u},\boldsymbol{u})\right>_{L^{2}}=0$.
Additionally, since the gradient operator and the divergence operator are each
other’s negated adjoints it can be stated that
$\left<\boldsymbol{u},\mathcal{G}p\right>_{L^{2}}=-\left<\mathcal{M}\boldsymbol{u},p\right>_{L^{2}}=0$
due to equation (2). Finally, the negative-definiteness of the diffusion
operator can be employed to simplify the evolution equation of the total
kinetic energy further to:

$\frac{dK}{dt}=-\nu\left|\left|\nabla\boldsymbol{u}\right|\right|_{L^{2}}^{2}\leq
0.$ (5)

Equation (5) implies that the kinetic energy, or the norm of the velocity
field $\boldsymbol{u}$, is a monotonically decreasing quantity for periodic or
homogeneous Dirichlet boundary conditions and is conserved in the inviscid
case ($\nu=0$). It will become clear soon that mimicking this property
discretely will be crucial to nonlinear stability of the FOM and ROMs.

Furthermore, it follows from integrating equation (1) over the domain $\Omega$
that:

$\frac{d\boldsymbol{P}}{dt}=0,$ (6)

for periodic boundary conditions and with $\boldsymbol{f}^{b}=0$. It is now
our goal to mimic these conservation properties at the discrete and especially
the reduced level.

###  2.2 Structure-preserving FOM

Equation (1) is discretized with a finite volume method (FVM) on a staggered
grid, resulting in a system of coupled ordinary differential equations (ODEs)
complemented by a set of linear constraints:

$\displaystyle\Omega_{h}\frac{d\boldsymbol{u}_{h}}{dt}+C_{h}(\boldsymbol{u}_{h})=-G_{h}\boldsymbol{p}_{h}+\nu
D_{h}\boldsymbol{u}_{h},$ (7) $\displaystyle M_{h}\boldsymbol{u}_{h}=0.$ (8)

Here $\boldsymbol{u}_{h}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{N}$ are the
numerical velocity values arranged in a vector,
$\Omega_{h}\in\mathbb{R}^{N\times N}$ is a diagonal matrix containing finite
volume sizes,
$C_{h}(\boldsymbol{u}_{h}):\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is the
spatial discretization of the nonlinear convection operator,
$G_{h}\in\mathbb{R}^{N\times N_{p}}$ is the spatial discretization of the
gradient operator,
$\boldsymbol{p}_{h}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{N_{p}}$ are the
numerical pressure values arranged in a vector, $D_{h}\in\mathbb{R}^{N\times
N}$ is the spatial discretization of the diffusion operator,
$M_{h}\in\mathbb{R}^{N_{p}\times N}$ is the spatial discretization of the
divergence operator, $N=N_{u}+N_{v}$ is the total number of velocity unknowns
and $N_{u}$ and $N_{v}$ are the numbers of velocity unknowns in the $x$ and
$y$ directions respectively (for the case of two-dimensional domains) and
$N_{p}$ is the number of pressure unknowns. Relations (7)-(8) are in turn
complemented by a vector of initial conditions
$\boldsymbol{u}_{0}=\boldsymbol{u}_{h}(0)$ and suitable boundary conditions.
The choice of boundary conditions for the flows considered in this work are
periodic.

We use the discretization as described in [71, 73] such that the following
three important properties of the continuous operators are inherited by their
spatial discretizations. First, the convection operator, which can be written
in quasi-linear form as
$C_{h}(\boldsymbol{u}_{h})=\widetilde{C}_{h}(\boldsymbol{u}_{h})\boldsymbol{u}_{h}$,
is skew-symmetric:

$\widetilde{C}_{h}(\boldsymbol{u}_{h})=-\widetilde{C}_{h}(\boldsymbol{u}_{h})^{T}.$
(9)

Second, the duality between gradient and divergence operators holds:

$G_{h}=-M_{h}^{T}.$ (10)

Third, the diffusive operator is symmetric negative-definite, allowing the
decomposition:

$D_{h}=-Q_{h}^{T}Q_{h}.$ (11)

We remark that the conservation properties discussed in this article are not
specific to the method described in [71, 73], any other finite volume
discretization of (1)-(2) satisfying the above properties can also be used.

Using these properties and the telescoping property of the FVM, discretization
(7)-(8) conserves discrete analogues of mass, total momentum and total kinetic
energy. To show this, we first define discrete total kinetic energy and
momentum on the basis of an inner product. To this end, the $\Omega_{h}$-inner
product is introduced. The $\Omega_{h}$-inner product and its induced norm are
defined by:

$\left<\boldsymbol{u},\boldsymbol{v}\right>_{\Omega_{h}}:=\left<\boldsymbol{u},\Omega_{h}\boldsymbol{v}\right>,\quad\left|\left|\boldsymbol{u}\right|\right|_{\Omega_{h}}:=\sqrt{\left<\boldsymbol{u},\boldsymbol{u}\right>_{\Omega_{h}}},$

where $\left<\cdot,\cdot\right>$ denotes the standard Euclidean inner product.
After deriving an evolution equation for the discrete total kinetic energy it
will be clear why this is a natural choice. Using the $\Omega_{h}$-inner
product, the discrete total kinetic energy
$K_{h}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ is defined as:

$K_{h}(t):=\frac{1}{2}\left|\left|\boldsymbol{u}_{h}\right|\right|_{\Omega_{h}}^{2}.$

The discrete total momentum
$\boldsymbol{P}_{h}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{d}$ is defined as:

$\left(\boldsymbol{P}_{h}(t)\right)_{i}=\left<\boldsymbol{e}_{i},\boldsymbol{u}_{h}\right>_{\Omega_{h}},$

where $\boldsymbol{e}_{i}\in\mathbb{R}^{N}$ is a vector of ones at the vector
indices where $\boldsymbol{u}_{h}$ has velocity components in direction
$i\in\\{1,...,d\\}$ and zeros elsewhere.

An evolution equation for discrete kinetic energy is found by temporal
differentiation of $K_{h}(t)$:

$\displaystyle\frac{dK_{h}}{dt}=\frac{1}{2}\frac{d}{dt}\left<\boldsymbol{u}_{h},\boldsymbol{u}_{h}\right>_{\Omega_{h}}$
$\displaystyle=-\left<\boldsymbol{u}_{h},\widetilde{C}_{h}(\boldsymbol{u}_{h})\boldsymbol{u}_{h}\right>-\left<\boldsymbol{u}_{h},G_{h}\boldsymbol{p}_{h}\right>+\nu\left<\boldsymbol{u}_{h},D_{h}\boldsymbol{u}_{h}\right>$
$\displaystyle=-\nu\left|\left|Q_{h}\boldsymbol{u}_{h}\right|\right|^{2}\leq
0,$ (12)

where $||\cdot||$ denotes the Euclidean norm. The simplification in the second
line follows from the application of properties (9)-(11). Indeed, evolution
equation (12) is a discrete analogue to equation (5) and shows that $K_{h}$,
and so the norm of $\boldsymbol{u}_{h}$, is a monotonically decreasing
quantity. The solution of the semi-discrete system of ODEs (7) can thus be
bounded from above by:

$\left|\left|\boldsymbol{u}(t)\right|\right|_{\Omega_{h}}\leq\left|\left|\boldsymbol{u}_{0}\right|\right|_{\Omega_{h}},$

guaranteeing stability of $\boldsymbol{u}_{h}(t)$.

It is also a well-known property of the FVM, referred to as the telescoping
property [84], that:

$\displaystyle\left(\frac{d\boldsymbol{P}_{h}}{dt}\right)_{i}=\frac{d}{dt}\left<\boldsymbol{e}_{i},\boldsymbol{u}_{h}\right>_{\Omega_{h}}$
$\displaystyle=-\left<\boldsymbol{e}_{i},\widetilde{C}_{h}(\boldsymbol{u}_{h})\boldsymbol{u}_{h}\right>-\left<\boldsymbol{e}_{i},G_{h}\boldsymbol{p}_{h}\right>+\nu\left<\boldsymbol{e}_{i},D_{h}\boldsymbol{u}_{h}\right>$
$\displaystyle=0,\quad\forall i\in\\{1,...,d\\}.$ (13)

This constitutes a discrete analogue to (6) and shows that discrete total
momentum is also conserved.

###  2.3 Structure-preserving ROM

To construct a reduced order model as in [73] we make the assumption that for
any $t\in[0,T]$ there are elements $\boldsymbol{u}_{r}(t)$ of a low-
dimensional linear subspace $\mathcal{V}\subset\mathbb{R}^{N}$ that accurately
approximate $\boldsymbol{u}_{h}(t)$. Here, we denote
$r:=\text{dim}(\mathcal{V})$ and the low dimensionality of $\mathcal{V}$
implies that $r\ll N$. Furthermore, we consider
$\mathcal{V}\subset\text{ker}(M_{h})$ such that all $\boldsymbol{u}_{r}(t)$
automatically satisfy condition (8). Let $\Phi\in\mathbb{R}^{N\times r}$ be a
POD basis for $\mathcal{V}$ that is orthonormal in the $\Omega_{h}$-inner
product, i.e. $\left<\Phi_{,i},\Phi_{,j}\right>_{\Omega_{h}}=\delta_{ij}$
111We use commas in the notation $(A_{,j})$ to the denote the $j^{\text{th}}$
column of matrix $A$, similarly $(A_{i,})$ denotes the $i^{\text{th}}$ row.,
then we write:

$\boldsymbol{u}_{h}(t)\approx\boldsymbol{u}_{r}(t):=\Phi\boldsymbol{a}(t)\in\mathcal{V}.$
(14)

Here $\boldsymbol{a}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{r}$ are the
generalized coordinates in $\mathcal{V}$ associated to the POD basis $\Phi$.
Substituting this approximation into the FOM (8) we obtain an overdetermined
system. To obtain a solvable system, we test the overdetermined system against
the POD modes in $\Phi$ and impose the Galerkin condition. This results in a
pressure-free222Note that since $\mathcal{V}\subset\text{ker}(M_{h})$, we can
write $\Phi^{T}G_{h}p_{h}=-(M_{h}^{T}\Phi)^{T}p_{h}=0$, making the ROM
velocity-only. ROM:

$\frac{d\boldsymbol{a}}{dt}=-\widetilde{C}_{r}(\boldsymbol{a})\boldsymbol{a}+\nu
D_{r}\boldsymbol{a},$ (15)

where
$\widetilde{C}_{r}(\boldsymbol{a}):=\Phi^{T}\widetilde{C}_{h}(\Phi\boldsymbol{a})\Phi$
and $D_{r}:=\Phi^{T}D_{h}\Phi$.

In order to consider structure-preservation we first define the following
quantities. A reduced order representation of kinetic energy is defined as:

$K_{r}(t):=\frac{1}{2}||\boldsymbol{u}_{r}||_{\Omega_{h}}^{2}=\frac{1}{2}||\boldsymbol{a}||^{2}.$
(16)

For total momentum a reduced order representation is:

$\left(\boldsymbol{P}_{r}(t)\right)_{i}:=\left<\boldsymbol{e}_{i},\boldsymbol{u}_{r}\right>_{\Omega_{h}}=\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},\boldsymbol{a}\right>.$

Noting that symmetry properties (9) and (11) are still satisfied by the
reduced operators in (15), the evolution equation for $K_{r}(t)$ follows by
differentiating (16):

$\displaystyle\frac{dK_{r}}{dt}=\frac{1}{2}\frac{d}{dt}\left<\boldsymbol{a},\boldsymbol{a}\right>$
$\displaystyle=-\left<\boldsymbol{a},\widetilde{C}_{r}(\boldsymbol{a})\boldsymbol{a}\right>+\nu\left<\boldsymbol{a},D_{r}\boldsymbol{a}\right>$
$\displaystyle=-\nu\left|\left|Q_{r}\boldsymbol{a}\right|\right|^{2}\leq 0,$
(17)

where $Q_{r}:=Q_{h}\Phi$. This equation is the ROM analogue of equation (5)
and shows that the norm of the generalized coordinates satisfies:

$\left|\left|\boldsymbol{a}(t)\right|\right|\leq\left|\left|\boldsymbol{a}_{0}\right|\right|,$

with equality for $\nu=0$.

To conserve total momentum at the reduced level we require
$\boldsymbol{e}_{i}\in\mathcal{V}\hbox{ }\forall i\in\\{1,...,d\\}$. This is
done by explicitly selecting
$\boldsymbol{e}_{i}/||\boldsymbol{e}_{i}||_{\Omega_{h}}$ as the first $d$ POD
modes in $\Phi$. An evolution equation for $\boldsymbol{P}_{r}(t)$ is then
found as:

$\displaystyle\left(\frac{d\boldsymbol{P}_{r}}{dt}\right)_{i}=\frac{d}{dt}\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},\boldsymbol{a}\right>$
$\displaystyle=-\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},\widetilde{C}_{r}(\boldsymbol{a})\boldsymbol{a}\right>+\nu\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},D_{r}\boldsymbol{a}\right>$
$\displaystyle=-\left<\boldsymbol{e}_{i},\widetilde{C}_{h}(\Phi\boldsymbol{a})\Phi\boldsymbol{a}\right>+\nu\left<\boldsymbol{e}_{i},D_{h}\Phi\boldsymbol{a}\right>$
(18) $\displaystyle=0.$

In the second line the telescoping property of the FOM operators was evoked.

Finally, a POD basis $\Phi$ satisfying all required criteria
($\text{span}(\Phi)\subset\text{ker}(M_{h})$, $\Omega_{h}$-orthonormality and
containing all $\boldsymbol{e}_{i}/||\boldsymbol{e}_{i}||_{\Omega_{h}}$ in its
columns) is constructed from snapshot data using an altered POD method as
described in [73]. The FOM and the ROM can both be integrated in time with
energy-conserving Runge-Kutta methods to conserve energy fully discretely (see
[73, 72] for further details).

##  3 A structure-preserving discrete empirical interpolation method

###  3.1 Introduction to DEIM

The ROM (15) as constructed in the previous section is $r$-dimensional, but
evaluating the convective terms in a naive fashion still requires a
computational effort that scales with the FOM dimension $N$. One approach to
alleviate this issue is the exact tensor decomposition. This method results in
an exact low-dimensional representation of the convection operator but scales
with $\mathcal{O}(r^{3})$, which becomes computationally expensive for
realistic values of $r$ [73]. A more efficient alternative approach, which
will be considered here, is the discrete empirical interpolation method (DEIM)
[23]. Using the DEIM we will construct a so-called hyper-reduced order model
(hROM).

The DEIM approximates $C_{h}(\boldsymbol{u})$ with elements from a low-
dimensional linear subspace $\mathcal{M}_{d}\subset\mathbb{R}^{N}$. Here, we
denote $m:=\text{dim}(\mathcal{M}_{d})$ and the low dimensionality of
$\mathcal{M}_{d}$ implies that $m\ll N$. Let $M\in\mathbb{R}^{N\times m}$ be
an orthonormal basis for $\mathcal{M}_{d}$ constructed from snapshot data of
the operator using the POD algorithm, where the snapshot data is gathered in a
matrix $\Xi\in\mathbb{R}^{N\times n_{s}}$ given as:

$\Xi=\left[C_{h}(\boldsymbol{u}_{h}(t^{0})),C_{h}(\boldsymbol{u}_{h}(t^{0}+\Delta
t)),...,C_{h}(\boldsymbol{u}_{h}(t^{0}+(n_{s}-1)\Delta t))\right].$ (19)

The basis $M$ thus satisfies $\left<M_{,i},M_{,j}\right>=\delta_{ij}$. Note
the distinction between the basis $M$ and the discrete divergence operator
$M_{h}$. The DEIM approximation then takes the form:

$C_{h}(\Phi\boldsymbol{a})\approx
M\boldsymbol{c}(\boldsymbol{a})\in\mathcal{M}_{d},$ (20)

where $\boldsymbol{c}(\boldsymbol{a}):\mathbb{R}^{r}\rightarrow\mathbb{R}^{m}$
are the generalized coordinates in $\mathcal{M}_{d}$ associated to basis $M$
and are referred to as DEIM coordinates. Using the DEIM approximation the ROM
takes the form:

$\frac{d\boldsymbol{a}}{dt}=-\Phi^{T}M\boldsymbol{c}(\boldsymbol{a})+\nu
D_{r}\boldsymbol{a}.$ (21)

The linear nature of the approximation allows the Galerkin projection
$\Phi^{T}M$ to be precomputed in an offline step. Hence, with the values of
the DEIM coordinates available, the ROM can be evaluated in a completely low-
dimensional fashion.

However, the problem (20) to calculate the DEIM coordinates is overdetermined.
Furthermore, the least-squares solution to system (20),
$\boldsymbol{c}=M^{T}C_{h}(\Phi\boldsymbol{a})$, requires an expensive
projection operation. To solve this, the DEIM bases its approximation on only
an $m$-dimensional subset of evaluations of $C_{h}(\Phi\boldsymbol{a})$ at
points on the grid, referred to as the measurement points, and denoted by
$\mathcal{P}$. More specifically, a measurement matrix $P\in\\{0,1\\}^{N\times
m}$ is defined that consists of selected columns of the $N\times N$ identity
matrix corresponding to the measurement points. To calculate the DEIM
coordinates the following system is then solved:

$P^{T}C_{h}(\Phi\boldsymbol{a})=P^{T}M\boldsymbol{c}.$ (22)

The formal solution to this system is
$\boldsymbol{c}=(P^{T}M)^{-1}P^{T}C_{h}(\Phi\boldsymbol{a})$. A greedy
measurement point selection procedure that ensures that $(P^{T}M)^{-1}$ exists
and that provides near-optimal approximations of $C_{h}(\Phi\boldsymbol{a})$
was proposed in [23]. Using this method, the DEIM coordinates can be
determined using only $m$ evaluations of $C_{h}(\Phi\boldsymbol{a})$ and
$(P^{T}M)^{-1}$ can be precomputed. This procedure results in a completely
low-dimensional ROM.

We will analyse the conservation properties of this ‘conventional’ DEIM
approach. Because the telescoping property of operators discretized using the
FVM is a linear property, it is conserved under application of the POD
algorithm. This may be shown by analysing the eigenvectors of the correlation
matrix $\Xi\Xi^{T}$. Using this property it can be shown that momentum is
conserved for a hROM using the DEIM:

$\displaystyle\left(\frac{d\boldsymbol{P}_{r}}{dt}\right)_{i}=\frac{d}{dt}\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},\boldsymbol{a}\right>$
$\displaystyle=-\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},\Phi^{T}M\boldsymbol{c}(\boldsymbol{a})\right>+\nu\left<\Phi^{T}\Omega_{h}\boldsymbol{e}_{i},D_{r}\boldsymbol{a}\right>$
$\displaystyle=-\left<\boldsymbol{e}_{i},M\boldsymbol{c}(\boldsymbol{a})\right>+\nu\left<\boldsymbol{e}_{i},D_{h}\Phi\boldsymbol{a}\right>$
(23) $\displaystyle=0.$

We now arrive at the key issue addressed in this article. This issue is that
the hROM (21) does not conserve reduced total kinetic energy, i.e. it can
occur that:

$\displaystyle\frac{dK_{r}}{dt}=\frac{1}{2}\frac{d}{dt}\left<\boldsymbol{a},\boldsymbol{a}\right>$
$\displaystyle=-\left<\boldsymbol{a},\Phi^{T}M\boldsymbol{c}(\boldsymbol{a})\right>-\nu\left|\left|Q_{r}\boldsymbol{a}\right|\right|^{2}>0,$
(24)

because generally
$\left<\boldsymbol{a},\Phi^{T}M\boldsymbol{c}(\boldsymbol{a})\right>$ does not
equal zero. This is because the DEIM-interpolated convection operator:

$\Phi^{T}M(P^{T}M)^{-1}P^{T}\widetilde{C}_{h}(\Phi\boldsymbol{a})\Phi,$ (25)

is no longer skew-symmetric. This can lead to instability as the norm of the
generalized coordinates $\boldsymbol{a}$ is not bounded from above. An
approach to retain the skew-symmetry of a DEIM approximation and therefore the
energy-conservation property of the ROM has been proposed in [55]. However,
this method scales equivalently to the exact tensor decomposition and is
therefore considered infeasible for realistic fluid dynamics applications. We
will present a new approach to make the hROM energy-conserving, while also
retaining the momentum conservation property.

###  3.2 A novel structure-preserving DEIM

The method we propose will be named the constrained least-squares discrete
empirical interpolation method (CLSDEIM). It is based on the realization that
skew-symmetry of the interpolated quasi-linear convection operator (25) is a
sufficient but not a necessary condition for energy-conservation. Instead, a
necessary condition to conserve reduced total kinetic energy is:

$\left<\boldsymbol{a},\Phi^{T}M\boldsymbol{c}(\boldsymbol{a})\right>=0.$ (26)

Condition (26) can be satisfied even if the Galerkin-projected DEIM operator
$\Phi^{T}M\boldsymbol{c}(\boldsymbol{a}):\mathbb{R}^{r}\rightarrow\mathbb{R}^{r}$
is not skew-symmetric and, when satisfied, results in the correct reduced
kinetic energy evolution equation (17). The new idea of the CLSDEIM is to
enforce this condition by posing the DEIM as a constrained minimization
problem.

The conventional DEIM finds the DEIM coordinates $\boldsymbol{c}$ by
minimizing the Euclidean norm between the nonlinearity and the DEIM
approximation in the measurement space. The Euclidean norm can be considered
as minimized since the difference in $\mathcal{P}$ i.e.
$P^{T}C_{h}(\Phi\boldsymbol{a})-P^{T}M\boldsymbol{c}$, is zero. The idea of
the proposed CLSDEIM is to employ this view of DEIM as a minimization problem
and constrain it to take place over the set $\mathcal{F}(\boldsymbol{a})$ of
DEIM approximations satisfying condition (26), defined using the DEIM
coordinates as:

$\mathcal{F}(\boldsymbol{a}):=\\{\boldsymbol{c}\in\mathbb{R}^{m}\hbox{
}|\hbox{ }\boldsymbol{a}^{T}\Phi^{T}M\boldsymbol{c}=0\\}.$

The set $\mathcal{F}(\boldsymbol{a})$ is referred to as the feasible set. As
DEIM approximations with
$\boldsymbol{c}(\boldsymbol{a})\in\mathcal{F}(\boldsymbol{a})$ satisfy
condition (26), the CLSDEIM produces approximations that conserve reduced
total kinetic energy. The constrained minimization problem to find the DEIM
coordinates $\boldsymbol{c}$ will be posed as the following linearly
constrained least-squares problem:

$\boxed{\boldsymbol{c}(\boldsymbol{a})=\operatorname*{arg\,min}_{\widetilde{\boldsymbol{c}}\in\mathbb{R}^{m}}\left|\left|P^{T}C_{h}(\Phi\boldsymbol{a})-P^{T}M\widetilde{\boldsymbol{c}}\right|\right|^{2}\quad\text{s.t.}\quad\boldsymbol{a}^{T}\Phi^{T}M\widetilde{\boldsymbol{c}}=0.}$
(27)

This means the CLSDEIM relaxes condition (22) of exact correspondence in
$\mathcal{P}$ between the FOM’s convection operator and the DEIM approximation
imposed by the conventional DEIM. Rather, the CLSDEIM minimizes the difference
between the approximation and the FOM operator in $\mathcal{P}$ while
simultaneously constraining the approximation to be energy-conserving.

Considering the objective function and the geometric interpretation of the
DEIM [25, 23], the CLSDEIM can be interpreted geometrically as an oblique
projection of $C_{h}(\Phi\boldsymbol{a})$ on the subspace
$\mathcal{M}_{\mathcal{F}}\subset\mathbb{R}^{N}$ of all DEIM approximations
$M\boldsymbol{c}$ with $\boldsymbol{c}\in\mathcal{F}(\boldsymbol{a})$. The
subspace $\mathcal{M}_{\mathcal{F}}$ is defined as:

$\mathcal{M}_{\mathcal{F}}=\mathcal{M}_{d}\cap\text{ker}(\left<\Phi\boldsymbol{a},\cdot\right>),$

where
$\left<\Phi\boldsymbol{a},\cdot\right>:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is
the functional taking the Euclidean inner product with $\Phi\boldsymbol{a}$.
Particularly, $\mathcal{M}_{\mathcal{F}}$ are all vectors in $\mathbb{R}^{N}$
that can be written as a linear combination of the columns of $M$ and satisfy
condition (26). Since $\mathcal{M}_{\mathcal{F}}$ is the intersection between
two linear subspaces of $\mathbb{R}^{N}$ it is also a linear subspace of
$\mathbb{R}^{N}$. Contrary to the conventional DEIM, the CLSDEIM also projects
obliquely through the measurement space $\mathcal{P}$. Hence, the orthogonal
projection of the DEIM residual $\boldsymbol{r}(t)\in\mathbb{R}^{N}$ on
$\mathcal{P}$, as given by:

$P^{T}\boldsymbol{r}(t)=P^{T}\left[C_{h}(\Phi\boldsymbol{a}(t))-M\boldsymbol{c}(t)\right],$

is generally not zero. However, the CLSDEIM projector
$\Pi_{\mathcal{M}_{\mathcal{F}}}:\mathbb{R}^{N}\rightarrow\mathcal{M}_{\mathcal{F}}$
does project $C_{h}(\Phi\boldsymbol{a})$ onto $\mathcal{M}_{\mathcal{F}}$ such
that the orthogonally projected residual $P^{T}\boldsymbol{r}(t)$ is minimal
in the Euclidean norm. This statement follows exactly from the formulation of
the CLSDEIM minimization problem (27).

Furthermore, the CLSDEIM basis $M$ and the measurement space $\mathcal{P}$ can
simply be found following the procedures of the conventional DEIM algorithm.
Indeed, as we are using the conventional DEIM basis the reduced total momentum
will also remain a conserved quantity for the CLSDEIM as was shown in equation
(23).

The constrained minimization problem (27) is solved using the method of
Lagrange multipliers [14]. The Lagrangian
$\mathcal{L}(\boldsymbol{c},\lambda):\mathbb{R}^{m}\times\mathbb{R}\rightarrow\mathbb{R}$
is defined as:

$\mathcal{L}(\boldsymbol{c},\lambda)=\left|\left|P^{T}C_{h}(\Phi\boldsymbol{a})-P^{T}M\boldsymbol{c}\right|\right|^{2}+\lambda\boldsymbol{a}^{T}\Phi^{T}M\boldsymbol{c},$

where $\lambda\in\mathbb{R}$ is a Lagrange multiplier. Taking partial
derivatives of the Lagrangian and setting them to zero leads to the following
system for the optima:

$\begin{bmatrix}2(P^{T}M)^{T}P^{T}M&(\boldsymbol{a}^{T}\Phi^{T}M)^{T}\\\
\boldsymbol{a}^{T}\Phi^{T}M&0\end{bmatrix}\begin{bmatrix}\boldsymbol{c}_{o}\\\
\lambda_{o}\end{bmatrix}=\begin{bmatrix}2(P^{T}M)^{T}P^{T}C_{h}(\Phi\boldsymbol{a})\\\
0\end{bmatrix}.$ (28)

Note that this matrix is symmetric and has constant coefficients with the
exception of the last row and column which depend on $\boldsymbol{a}$. Solving
(33) can be done explicitly using the inverse of a symmetric $2\times 2$ block
matrix [49]:

$\begin{bmatrix}A&B\\\
B^{T}&0\end{bmatrix}^{-1}=\begin{bmatrix}A^{-1}-A^{-1}B(B^{T}A^{-1}B)^{-1}B^{T}A^{-1}&A^{-1}B(B^{T}A^{-1}B)^{-1}\\\
(B^{T}A^{-1}B)^{-1}B^{T}A^{-1}&-(B^{T}A^{-1}B)^{-1}\end{bmatrix},$

where $A$ is symmetric. Taking:

$\displaystyle A:=2(P^{T}M)^{T}P^{T}M$ $\displaystyle\in\mathbb{R}^{m\times
m}$ $\displaystyle
B:=(\boldsymbol{a}^{T}\Phi^{T}M)^{T}=M^{T}\Phi\boldsymbol{a}$
$\displaystyle\in\mathbb{R}^{m\times 1},$

the DEIM coordinates $\boldsymbol{c}_{o}$ solving (27) are determined by:

$\boldsymbol{c}(\boldsymbol{a})=\left(A^{-1}-A^{-1}B(B^{T}A^{-1}B)^{-1}B^{T}A^{-1}\right)2(P^{T}M)^{T}P^{T}C_{h}(\Phi\boldsymbol{a}).$

We use the notation $B:=\boldsymbol{b}(\boldsymbol{a})$ to explicitly denote
the dependence of $B$ on $\boldsymbol{a}$ and the fact that it can be treated
as a vector. Since
$B^{T}A^{-1}B=\boldsymbol{b}(\boldsymbol{a})^{T}A^{-1}\boldsymbol{b}(\boldsymbol{a})$
and
$B^{T}A^{-1}2(P^{T}M)^{T}P^{T}C_{h}(\Phi\boldsymbol{a})=\boldsymbol{b}(\boldsymbol{a})^{T}(P^{T}M)^{-1}P^{T}C_{h}(\Phi\boldsymbol{a})$
are simply scalars, the final expression for $\boldsymbol{c}_{o}$ is:

$\boxed{\text{CLSDEIM:}\qquad\boldsymbol{c}(\boldsymbol{a})=(P^{T}M)^{-1}P^{T}C_{h}(\Phi\boldsymbol{a})-\frac{\boldsymbol{b}(\boldsymbol{a})^{T}(P^{T}M)^{-1}P^{T}C_{h}(\Phi\boldsymbol{a})}{(\boldsymbol{b}(\boldsymbol{a})^{T}A^{-1}\boldsymbol{b}(\boldsymbol{a}))}A^{-1}\boldsymbol{b}(\boldsymbol{a}).}$
(29)

The DEIM coordinates found from equation (29) result in a DEIM approximation
$M\boldsymbol{c}$ satisfying condition (26) and consequently also equation
(17). All terms in (29) are either inner products or matrix-vector products
between low-dimensional vectors or matrices and vectors respectively, hence
the CLSDEIM can be solved in a low-dimensional fashion. Note that the first
term on the right hand side of equation (29) is the solution to equation (22).
Hence, equation (29) can be interpreted as the conventional DEIM with a
correction term. The correction term projects the conventional DEIM
coordinates on $\text{ker}(\left<\boldsymbol{b}(\boldsymbol{a}),\cdot\right>)$
such that condition (26) is satisfied and minimization problem (27) is solved.

###### Remark 1

Both implicit and explicit Runge-Kutta (RK) methods will be considered to
integrate the proposed hROMs. To achieve exact energy conservation (for
inviscid flows), implicit energy-conserving RK methods need to be used [72,
73]. Using such methods and structure-preserving hyper-reduction, the change
in reduced total kinetic energy $K_{r}$ is given by:

$\boxed{K_{r}^{n+1}-K_{r}^{n}=-\Delta
t\sum_{i=1}^{s}b_{i}\nu\left|\left|Q_{r}\boldsymbol{A}_{i}\right|\right|^{2},}$
(30)

which is the time-discrete equivalent of (17). Here, $s\in\mathbb{N}$ and
$b_{i}\in\mathbb{R}$ are parameters of the given RK method and
$\boldsymbol{A}_{i}\in\mathbb{R}^{r}$ are the generalized coordinates at the
$i^{\text{th}}$ RK stage. In the inviscid case it can be seen from (30) that
the hROMs conserve reduced kinetic energy. In the viscous case the resulting
hROMs will be nonlinearly stable when $b_{i}\geq 0\hbox{ }\forall\hbox{
}i\in\\{1,...,s\\}$, which holds for the Gauss-Legendre RK methods (see [73,
72] for further details).

###  3.3 Preventing overfitting errors

Interpolation-based methods like the DEIM and CLSDEIM can be prone to
overfitting when noise is present [38]. It was shown in [63] that adding
Gaussian noise $\epsilon\in\mathbb{R}^{m}$ to the measurements, i.e.
$P^{T}C_{h}(\boldsymbol{u})+\epsilon$, results in error bounds on the DEIM
approximation of $C_{h}(\boldsymbol{u})$ that grow as $\mathcal{O}(\sqrt{m})$.
Similar overfitting errors may occur in the presence of numerical errors in
$\boldsymbol{a}$ as we will show using experiments later. As a result of such
overfitting errors, numerical errors in the solution may accumulate, further
harming the capability of the DEIM to reconstruct the correct
$C_{h}(\Phi\boldsymbol{a})$ of the following time steps. This is especially
notable when the interpolation procedure is carried out many times, which is
the case when solving time-dependent problems. We remark that this can happen
even in the case of fully-discrete energy-conserving schemes since energy
conservation is a restriction on the norm of the solution and not on its
accuracy. In this case, overfitting errors typically manifest themselves in
non-convergence of the nonlinear solver. To address overfitting we therefore
propose two further stabilization approaches. The first one is oversampling,
which has recently gained significant attention [63, 69]. The second approach
is regularization, which is a well-known technique from optimization and has
thus far, to the author’s knowledge, not been applied to the DEIM. This latter
aproach has the benefit that additional measurement points are not necessary
and therefore reduces the computational burden compared to oversampling.

####  3.3.1 Oversampling

When using oversampling we replace the interpolation of DEIM (22) and CLSDEIM
(27) with a regression procedure. This is achieved by considering a larger
measurement space than the dimension of the DEIM space $\mathcal{M}_{d}$, thus
$m_{p}:=\text{dim}\left(\mathcal{P}\right)>m$. As a result, the measurement
matrix becomes $P\in\\{0,1\\}^{N\times m_{p}}$. The formulation of the
underlying optimization problem remains the same as optimization problem (27),
but now the larger measurement matrix is used. Due to the fact that we are
trying to fit a linear combination of $m$ DEIM modes to $m_{p}$ measurements
with $m_{p}>m$, the resulting DEIM approximation $M\boldsymbol{c}$ is less
sensitive to noise in measurements of $C_{h}(\Phi\boldsymbol{a})$ in
$\mathcal{P}$. This was also analysed in [63], where it was shown that
oversampling using randomized measurement points resulted in a contribution of
noise to the error bound of the approximation of $C_{h}(\boldsymbol{u})$ that
scales with $\mathcal{O}\left(\sqrt{m/m_{p}}\right)$. Hence, overfitting
errors as a result of noise can be mitigated by choosing $m_{p}$ sufficiently
large compared to $m$. We will refer to CLSDEIM being extended with
oversampling as OCLSDEIM. The DEIM coordinates as calculated using the
OCLSDEIM require solving optimization problem (27) with the new measurement
matrix $P$. The solution procedure of this problem is completely analogous to
that of the CLSDEIM (27) and the final expression for $\boldsymbol{c}_{o}$ is:

$\boxed{\text{OCLSDEIM:}\qquad\boldsymbol{c}(\boldsymbol{a})=(P^{T}M)^{\dagger}P^{T}C_{h}(\Phi\boldsymbol{a})-\frac{\boldsymbol{b}(\boldsymbol{a})^{T}(P^{T}M)^{\dagger}P^{T}C_{h}(\Phi\boldsymbol{a})}{(\boldsymbol{b}(\boldsymbol{a})^{T}A^{-1}\boldsymbol{b}(\boldsymbol{a}))}A^{-1}\boldsymbol{b}(\boldsymbol{a}),}$
(31)

where instead of $(P^{T}M)^{-1}$ the Moore-Penrose pseudoinverse
$(P^{T}M)^{\dagger}$ appears.

The measurement point selection procedure will be as follows. The first $m$
measurement points are found using the classic DEIM procedure described in
[23]. The next $m_{p}-m$ points are generated from a uniform random
distribution where care is taken to not select repeating indices. This
procedure strikes a good balance between efficiency and accuracy and is simple
to implement. Alternatives are possible, see for example the following
articles for deterministic approaches to measurement point selection: [63, 69,
86, 50].

####  3.3.2 Generalized Tikhonov regularization

When overfitting occurs, DEIM modes interpolate measurements containing noise
or numerical errors resulting from time integration, negatively effecting
accuracy in subsequent time steps. Although the DEIM measurement points are
chosen to minimize $||(P^{T}M)^{-1}||$ thereby mitigating the effect of
overfitting [23], accumulation of errors can still take place as we will see
during the numerical experiments. This accumulation can eventually result in
severe errors and even numerical instability. To prevent this issue we propose
regularizing the least-squares problem using a generalized Tikhonov
regularization [38]. Using such a regularization procedure we can penalize
undesired properties of solutions to the DEIM optimization problem by adding
extra terms to its associated cost function. In this article we will propose
to use regularization terms of the form:

$\boxed{\boldsymbol{c}(\boldsymbol{a})=\operatorname*{arg\,min}_{\widetilde{\boldsymbol{c}}\in\mathbb{R}^{m}}\left|\left|P^{T}C_{h}(\Phi\boldsymbol{a})-P^{T}M\widetilde{\boldsymbol{c}}\right|\right|^{2}+\alpha\left|\left|\widetilde{\boldsymbol{c}}-\boldsymbol{\mu}\right|\right|_{\Gamma}^{2}\quad\text{s.t.}\quad\boldsymbol{a}^{T}\Phi^{T}M\widetilde{\boldsymbol{c}}=0.}$
(32)

Here
$||\boldsymbol{c}||_{\Gamma}:=\sqrt{\left<\boldsymbol{c},\Gamma\boldsymbol{c}\right>}$
with $\Gamma\in\mathbb{R}^{m\times m}$ symmetric positive definite (SPD),
$\boldsymbol{\mu}\in\mathbb{R}^{m}$ is some desired or reference state and
$\alpha\in\mathbb{R}^{+}$ is a hyperparameter. The minimization problem can
again be solved using the method of Lagrange multipliers. The Lagrangian
$\mathcal{L}(\boldsymbol{c},\lambda;\alpha):\mathbb{R}^{m}\times\mathbb{R}\rightarrow\mathbb{R}$
of this minimization is defined as:

$\mathcal{L}(\boldsymbol{c},\lambda;\alpha)=\left|\left|P^{T}C_{h}(\Phi\boldsymbol{a})-P^{T}M\boldsymbol{c}\right|\right|^{2}+\alpha\left|\left|\boldsymbol{c}-\boldsymbol{\mu}\right|\right|_{\Gamma}^{2}+\lambda\boldsymbol{a}^{T}\Phi^{T}M\boldsymbol{c},$

where $\lambda\in\mathbb{R}$ is a Lagrange multiplier. Taking partial
derivatives of the Lagrangian and setting them to zero leads to the following
system for the optima:

$\begin{bmatrix}2\left((P^{T}M)^{T}P^{T}M+\alpha\Gamma\right)&(\boldsymbol{a}^{T}\Phi^{T}M)^{T}\\\
\boldsymbol{a}^{T}\Phi^{T}M&0\end{bmatrix}\begin{bmatrix}\boldsymbol{c}_{o}\\\
\lambda_{o}\end{bmatrix}=\begin{bmatrix}2(P^{T}M)^{T}P^{T}C_{h}(\Phi\boldsymbol{a})+2\alpha\Gamma\boldsymbol{\mu}\\\
0\end{bmatrix}.$ (33)

Again, this matrix is symmetric and has constant coefficients with the
exception of the last row and column which depend on $\boldsymbol{a}$.
Denoting:

$A_{\Gamma}:=(P^{T}M)^{T}P^{T}M+\alpha\Gamma\in\mathbb{R}^{m\times m},$

the DEIM coordinates as calculated by the regularized CLSDEIM are given as:

$\boxed{\begin{gathered}\text{Regularized CLSDEIM: }\\\
\boldsymbol{c}(\boldsymbol{a})=A_{\Gamma}^{-1}\left((P^{T}M)^{T}P^{T}C_{h}(\Phi\boldsymbol{a})+\alpha\Gamma\boldsymbol{\mu}\right)-\frac{\boldsymbol{b}(\boldsymbol{a})^{T}A_{\Gamma}^{-1}\left((P^{T}M)^{T}P^{T}C_{h}(\Phi\boldsymbol{a})+\alpha\Gamma\boldsymbol{\mu}\right)}{(\boldsymbol{b}(\boldsymbol{a})^{T}A_{\Gamma}^{-1}\boldsymbol{b}(\boldsymbol{a}))}A_{\Gamma}^{-1}\boldsymbol{b}(\boldsymbol{a}).\end{gathered}}$
(34)

Note that this equation also applies to a combination of oversampling and
regularization. Taking $\alpha=0$ in equation (34) we obtain the normal
CLSDEIM since $A_{\Gamma}^{-1}(P^{T}M)^{T}=(P^{T}M)^{\dagger}=(P^{T}M)^{-1}$
in this case. A data-driven method to calculate $\Gamma$ and
$\boldsymbol{\mu}$ that uses the knowledge we have on the DEIM coordinates of
our snapshot data $\Xi$ will be proposed in the following section. The
hyperparameter tuning of $\alpha$ was carried out manually for simplicity. A
more advanced approach to find well-performing hyperparameters for regularized
minimization problems like (32) is generalized cross-validation [33], which we
suggest for future work.

####  3.3.3 Mahalanobis distance

The method we propose to determine the regularization term is to use the
Mahalanobis distance
$d_{M}(\boldsymbol{c}):\mathbb{R}^{m}\rightarrow\mathbb{R}^{+}$ [38] of the
solution $\boldsymbol{c}$ with respect to the data matrix
$C=M^{T}\Xi\in\mathbb{R}^{m\times n_{s}}$, note the difference with the
discrete convection operator $C_{h}$. The Mahalanobis distance is computed as:

$d_{M}(\boldsymbol{c})=\left|\left|\boldsymbol{c}-\boldsymbol{\mu}_{C}\right|\right|_{S^{-1}},$

and thus corresponds to taking $\Gamma=S^{-1}$ and
$\boldsymbol{\mu}=\boldsymbol{\mu}_{C}$ in optimization problem (32). Here,
the matrix $S\in\mathbb{R}^{m\times m}$ and vector
$\boldsymbol{\mu}_{C}\in\mathbb{R}^{m}$ are the covariance matrix and the mean
associated to the data in $C$, respectively. This method is based on the idea
of promoting solutions $\boldsymbol{c}$ that are similar to those observed in
the data matrix $C$ (which has the DEIM coordinates of the snapshots in $\Xi$
as its columns). Here, similar is in the sense of originating from the same
distribution. In fact, regularization using the Mahalanobis distance $d_{M}$
corresponds to assuming a prior distribution on the optimal $\boldsymbol{c}$
that can be accurately characterized by a mean and a covariance matrix $S$.
The Mahalanobis distance is then the multidimensional generalization of the
number of standard deviations an observation is removed from the
distribution’s mean.

In this paper we will concern ourselves with solution reproduction problems
for simplicity, and not with parameterized problems. We will see that, as
numerical errors accumulate during time integration, DEIM coordinates take
increasingly different values from their reference values in $C$ when using
the DEIM or CLSDEIM. Hence, such values will contribute to a large
$d_{M}(\boldsymbol{c})$ and as a result be less optimal than a more well-
behaved solution given the regularized CLSDEIM loss-function (32). We assume
that this idea can be extended to the parametric case, i.e. that the DEIM
coordinates of a parametric problem can also be treated as generated by a
single prior distribution. We will refer to the CLSDEIM regularized by the
Mahalanobis distance as the MCLSDEIM.

The diagonal entries $(S)_{ii}$ of $S$ can be very small for $i\gg 1$, due to
the fact that the DEIM coordinates for these $i$ are generally very small.
Furthermore $S$ may be nearly rank deficient. Both these aspects can result in
an ill-conditioned $S^{-1}$ having large entries. Future research will involve
improving the conditioning of the regularization term using e.g. a penalized
Mahalanobis distance [38]. In the current work we solve this issue by taking
very small hyperparameter values $\alpha$.

##  4 A structure-preserving temporal localization method

Having constructed new energy-conserving and robust DEIM methods, we now
extend the feasibility of our approach with a new structure-preserving
temporal localization method. Namely, it is well known that as flows become
more convection dominated ($\nu\rightarrow 0$) the Kolmogorov N-width decay
becomes increasingly slower. This makes it more difficult for conventional POD
to generate a low-dimensional basis set that accurately captures a significant
amount of energy contained in both $\Xi$ and the solution snapshots
$X\in\mathbb{R}^{N\times n_{s}}$. To deal with this problem researchers have
proposed several solutions, one of which is the application of the principal
interval decomposition (PID) for both the construction of DEIM and POD spaces,
$\mathcal{M}_{d}$ and $\mathcal{V}$ respectively [21, 62, 42, 13, 3, 35, 26].
However, as far as we know structure-preservation has not been taken into
account in existing PID approaches. In what follows, we will propose a method
to preserve structure throughout a full time integration using a hROM with
temporal localization.

###  4.1 Temporal localization using PID

The PID decomposes the snapshot sets over $n$ intervals in time
$[t_{i},t_{i+1}]$ and applies the POD algorithm to the individual intervals.
The idea is that by calculating modes tailored to specific intervals, the
local timescales within the respective intervals are captured significantly
better than by a set of modes calculated from the full set of snapshots [13,
21]. Based on snapshot sets:

$X=[X_{1},X_{2},...,X_{n}],\quad\Xi=[\Xi_{1},\Xi_{2},...,\Xi_{n}],$

the PID provides sets of POD modes:

$\Phi_{i},M_{i},\hbox{ }i\in\\{1,...,n\\},$ (35)

applicable to use at times $t\in[t_{i},t_{i+1}]$ within their respective
intervals. The local POD modes $\Phi_{i}$ follow from an SVD of the local
snapshot matrix $X_{i}\in\mathbb{R}^{N\times n_{s}^{i}}$, where
$n_{s}^{i}\in\mathbb{N}$ is the number of snapshots in the $i^{\text{th}}$
interval. The temporally localized DEIM measurement space $\mathcal{P}_{i}$
for the $i^{\text{th}}$ interval is determined solely based on the DEIM modes
in $M_{i}$ using the algorithm as in [23] specifically for $M_{i}$.
Furthermore, we have
$\mathcal{V}^{i}:=\text{span}(\Phi_{i})\subset\mathbb{R}^{N}$ and
$\mathcal{M}_{d}^{i}:=\text{span}(M_{i})\subset\mathbb{R}^{N}$.

Setting up the hROM using the PID, the dynamical system (21) takes the form:

$\frac{d\boldsymbol{a}_{i}}{dt}=-\Phi_{i}^{T}M_{i}\boldsymbol{c}_{i}(\boldsymbol{a}_{i})+\nu
D_{r}^{i}\boldsymbol{a}_{i}\quad t\in[t_{i},t_{i+1}],$ (36)

where $D_{r}^{i}:=\Phi_{i}^{T}D_{h}\Phi_{i}$. Note that we also introduce
subscripts for the generalized and DEIM coordinates $\boldsymbol{a}_{i}$ and
$\boldsymbol{c}_{i}$ respectively, as they are only valid during the interval
$[t_{i},t_{i+1}]$. The DEIM coordinates are calculated using any of the
previously discussed methods using the appropriate measurement spaces and DEIM
and POD bases. It should be noted that all the conservation properties of the
structure-preserving DEIM methods hold within intervals. This is trivial as
nothing changes in the way the DEIM coordinates are calculated. Instead, the
art is to develop interface conditions such that during the transition from
one interval to the next, these conservation properties are still satisfied.

Note that in the subsequent sections a hROM using the PID and a certain DEIM
algorithm will be referred to as PID-CLSDEIM hROM (when for example the
CLSDEIM is used without regularization).

###  4.2 Interface conditions

The bases in (35) are used only during the interval $t\in[t_{i},t_{i+1}]$. At
$t=t_{i+1}$ an interface is reached and it is required to transition from a
solution described in terms of
$\boldsymbol{u}_{r}^{-}:=\Phi_{i}\boldsymbol{a}_{i}(t_{i+1})\in\mathcal{V}^{i}$
to one described in terms of
$\boldsymbol{u}_{r}^{+}:=\Phi_{i+1}\boldsymbol{a}_{i+1}(t_{i+1})\in\mathcal{V}^{i+1}$,
where $\boldsymbol{a}_{i}\in\mathbb{R}^{r_{i}}$,
$\boldsymbol{a}_{i+1}\in\mathbb{R}^{r_{i+1}}$ and in general
$\text{dim}(\mathcal{V}^{i})=r_{i}\neq r_{i+1}=\text{dim}(\mathcal{V}^{i+1})$.
During the numerical experiments in this article we will consider the case
$r_{i}=r\hbox{ }\forall\hbox{ }i\in\\{1,...,n\\}$. To determine the interval
parameters, one could use heuristic methods [21] or a clustering algorithm
[35], but this is outside the scope of this article.

We will define how to transition between intervals based on an interface
condition. Using the conventional PID [3] the following interface condition is
defined:

$\text{find}\quad\boldsymbol{a}_{i+1}(t_{i+1})\quad\text{s.t.}\quad\left<\boldsymbol{u}_{r}^{+}-\boldsymbol{u}_{r}^{-},\boldsymbol{v}\right>_{\Omega_{h}}=0\quad\forall\hbox{
}\boldsymbol{v}\in\mathcal{V}^{i+1},$ (37)

which is stable in the sense that it does not lead to an energy increase, as
we will show. This condition is satisfied by the $\Omega_{h}$-orthogonal
projection of $\boldsymbol{u}_{r}^{-}$ onto $\mathcal{V}^{i+1}$:

$\boldsymbol{a}_{i+1}(t_{i+1})=\Phi_{i+1}^{T}\Omega_{h}\boldsymbol{u}_{r}^{-}=\Phi_{i+1}^{T}\Omega_{h}\Phi_{i}\boldsymbol{a}_{i}(t_{i+1}).$
(38)

Note that the interface matrix
$T_{i}:=\Phi_{i+1}^{T}\Omega_{h}\Phi_{i}\in\mathbb{R}^{r_{i+1}\times r_{i}}$
can be precomputed in an offline stage to obtain a cheap online evaluation.
The kinetic energy
$K_{r}^{+}:=\frac{1}{2}\left|\left|\boldsymbol{u}_{r}^{+}\right|\right|_{\Omega_{h}}^{2}=\sum_{k}\frac{1}{2}\left<\boldsymbol{u}_{r}^{-},(\Phi_{i+1})_{,k}\right>_{\Omega_{h}}^{2}$
cannot increase beyond
$K_{r}^{-}:=\frac{1}{2}\left|\left|\boldsymbol{u}_{r}^{-}\right|\right|_{\Omega_{h}}^{2}$
due to the Bessel inequality [45], which states that for the (complete) inner
product space defined as
$\mathcal{H}_{\Omega_{h}}:=\left(\mathbb{R}^{N},\left<\cdot,\cdot\right>_{\Omega_{h}}\right)$
and orthonormal sequence $\left((\Phi_{i+1})_{,k}\right)$ with
$(\Phi_{i+1})_{,k}\in\mathcal{H}_{\Omega_{h}},\hbox{ }k\in\\{1,...,r_{i+1}\\}$
we have the following:

$\text{Bessel Inequality in
$\mathcal{H}_{\Omega_{h}}$:}\quad\sum_{k=1}^{r_{i+1}}\left<\boldsymbol{u}_{r}^{-},(\Phi_{i+1})_{,k}\right>_{\Omega_{h}}^{2}\leq\left|\left|\boldsymbol{u}_{r}^{-}\right|\right|_{\Omega_{h}}^{2}.$

Hence, application of the conventional PID using interface condition (37)
cannot lead to an increase in kinetic energy for any number of intervals $n$,
making it nonlinearly stable. Furthermore, the solution is also optimal in the
sense that (38) is also the solution to the minimization problem:

$\boldsymbol{a}_{i+1}(t_{i+1})=\operatorname*{arg\,min}_{\boldsymbol{a}\in\mathbb{R}^{r_{i+1}}}\left|\left|\boldsymbol{u}_{r}^{-}-\Phi_{i+1}\boldsymbol{a}\right|\right|_{\Omega_{h}}^{2}.$
(39)

In summary, the conventional PID is nonlinearly stable. However, kinetic
energy and momentum are in general not exactly conserved when transitioning
between intervals using condition (37), which can lead to artificial
(numerical) dissipation of the solution. To solve this issue of lack of
conservation, we propose the following new structure-preserving interface
condition.

###  4.3 Structure-preserving interface conditions

To consider structure preservation over interfaces we define
$\left(\boldsymbol{P}_{r}^{+}\right)_{k}:=\left<\boldsymbol{e}_{k},\boldsymbol{u}_{r}^{+}\right>_{\Omega_{h}}$
and similarly
$\left(\boldsymbol{P}_{r}^{-}\right)_{k}:=\left<\boldsymbol{e}_{k},\boldsymbol{u}_{r}^{-}\right>_{\Omega_{h}}$.
In order to conserve energy and momentum over interfaces we benefit from the
optimization formulation (39) underlying the non-conserving interface
condition (37). We propose to constrain the optimization problem (39) to
conserve $K_{r}$ and $\boldsymbol{P}_{r}$, leading to:

$\boxed{\boldsymbol{a}_{i+1}(t_{i+1})=\operatorname*{arg\,min}_{\boldsymbol{a}\in\mathbb{R}^{r_{i+1}}}\left|\left|\boldsymbol{u}_{r}^{-}-\Phi_{i+1}\boldsymbol{a}\right|\right|_{\Omega_{h}}^{2}\quad\text{s.t.}\quad\frac{1}{2}\left|\left|\boldsymbol{a}\right|\right|^{2}=K_{r}^{-},\quad\left<\boldsymbol{e}_{k},\Phi_{i+1}\boldsymbol{a}\right>_{\Omega_{h}}=\left(\boldsymbol{P}_{r}^{-}\right)_{k}\forall
k\in\\{1,...,d\\}.}$ (40)

The constraints on this problem can be simplified by realizing that the
momentum component $(\boldsymbol{P}_{r})_{k}$ in direction $k$ is fully
determined by the $k^{\text{th}}$ component $a_{k}$ of the generalized
coordinates, namely:

$\displaystyle\left(\boldsymbol{P}_{r}\right)_{k}=\left<\boldsymbol{e}_{k},\sum_{j=1}^{r}\Phi_{,j}a_{j}\right>_{\Omega_{h}}=\left<\boldsymbol{e}_{k},\frac{\boldsymbol{e}_{k}}{\left|\left|\boldsymbol{e}_{k}\right|\right|_{\Omega_{h}}}a_{k}\right>_{\Omega_{h}}=\left|\left|\boldsymbol{e}_{k}\right|\right|_{\Omega_{h}}a_{k},$

due to $\Omega_{h}$-orthogonality. Furthermore, the first $d$ POD modes for
every interval are the same and given by the momentum-conserving modes
$\boldsymbol{e}_{k}/\left|\left|\boldsymbol{e}_{k}\right|\right|_{\Omega_{h}}$.
As a result, the momentum-conserving interface constraint can only be
satisfied by setting
$(\boldsymbol{a}_{i+1})_{k}=(\boldsymbol{a}_{i})_{k}\hbox{ }\forall
k\in\\{1,...,d\\}$. The following matrices and vectors are now defined:

$\widetilde{\Phi}_{i}:=\left[(\Phi_{i})_{,d+1},\hbox{ }...\hbox{
},\hbox{}(\Phi_{i})_{,r_{i}}\right],\quad\widetilde{\boldsymbol{a}}_{i}^{T}:=\left[(\boldsymbol{a}_{i})_{d+1},\hbox{
}...\hbox{ },\hbox{ }(\boldsymbol{a}_{i})_{r_{i}}\right].$

Since the momentum conservation constraint fully determines the first $d$
components of the generalized coordinates $\boldsymbol{a}_{i+1}$ we can
simplify the optimization problem (40) by formulating it in terms of the new
$\widetilde{(\cdot)}$ variables. Note that at any interface:

$\displaystyle\boldsymbol{u}_{r}^{+}-\boldsymbol{u}_{r}^{-}$
$\displaystyle=\sum_{k=d+1}^{r_{i+1}}(\Phi_{i+1})_{,k}(\boldsymbol{a}_{i+1})_{k}-\sum_{l=d+1}^{r_{i}}(\Phi_{i})_{,l}(\boldsymbol{a}_{i})_{l}+\sum_{j=1}^{d}\frac{\boldsymbol{e}_{j}}{||\boldsymbol{e}_{j}||_{\Omega_{h}}}((\boldsymbol{a}_{i+1})_{j}-(\boldsymbol{a}_{i})_{j})$
$\displaystyle=\widetilde{\Phi}_{i+1}\widetilde{\boldsymbol{a}}_{i+1}-\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}$

due to the momentum-conserving interface conditions. The kinetic energy
constraint can also be simplified to
$||\widetilde{\boldsymbol{a}}_{i+1}||^{2}=||\widetilde{\boldsymbol{a}}_{i}||^{2}$
since the first $d$ components of the generalized coordinates
$\boldsymbol{a}_{i+1}$ and $\boldsymbol{a}_{i}$ are the same. The reduced
minimization problem for the remaining $r_{i+1}-d$ components of the
generalized coordinates $\boldsymbol{a}_{i+1}(t_{i+1})$ then becomes:

$\widetilde{\boldsymbol{a}}_{i+1}(t_{i+1})=\operatorname*{arg\,min}_{\boldsymbol{a}\in\mathbb{R}^{r_{i+1}-d}}\left|\left|\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}-\widetilde{\Phi}_{i+1}\boldsymbol{a}\right|\right|_{\Omega_{h}}^{2}\quad\text{s.t.}\quad\left|\left|\boldsymbol{a}\right|\right|^{2}=\left|\left|\widetilde{\boldsymbol{a}}_{i}\right|\right|^{2}.$
(41)

As before, we solve this constrained minimization problem using the method of
Lagrange multipliers [14]. The Lagrangian
$\mathcal{L}(\boldsymbol{a},\lambda):\mathbb{R}^{r_{i+1}-d}\times\mathbb{R}\rightarrow\mathbb{R}$
of this minimization is as follows:

$\mathcal{L}(\boldsymbol{a},\lambda)=\boldsymbol{a}^{T}\boldsymbol{a}-2\boldsymbol{a}^{T}\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}+\lambda(\boldsymbol{a}^{T}\boldsymbol{a}-\widetilde{\boldsymbol{a}}_{i}^{T}\widetilde{\boldsymbol{a}}_{i}),$

where
$\widetilde{T}_{i}:=\widetilde{\Phi}_{i+1}^{T}\Omega_{h}\widetilde{\Phi}_{i}$
and $\lambda\in\mathbb{R}$ is a Lagrange multiplier. Taking partial
derivatives of the Lagrangian and setting them to zero leads to the following
nonlinear system for the optima ($\boldsymbol{a}_{o},\lambda_{o}$):

$\displaystyle\nabla_{a}\mathcal{L}=(1+\lambda_{o})\boldsymbol{a}_{o}-\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}=0$
$\displaystyle\nabla_{\lambda}\mathcal{L}=\boldsymbol{a}_{o}^{T}\boldsymbol{a}_{o}-\widetilde{\boldsymbol{a}}_{i}^{T}\widetilde{\boldsymbol{a}}_{i}=0.$

Rewriting the equation in the first line, an expression for the generalized
coordinates at the optimum $\boldsymbol{a}_{o}$ is obtained. This expression
can be substituted into the second optimality condition to solve for
$\frac{1}{1+\lambda_{o}}$. Finally, this results in:

$\widetilde{\boldsymbol{a}}_{i+1}(t_{i+1})=\boldsymbol{a}_{o}=\frac{||\widetilde{\boldsymbol{a}}_{i}||}{\left|\left|\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}\right|\right|}\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}.$
(42)

To interpret this solution, we consider the subset
$\mathcal{S}(\delta):=\\{\boldsymbol{u}\in\mathcal{H}_{\Omega_{h}}\hbox{
}|\hbox{ }||\boldsymbol{u}||_{\Omega_{h}}=\delta\\}$ of all vectors in
$\mathcal{H}_{\Omega_{h}}$ with the same norm of $\delta\in\mathbb{R}^{+}$
(and thus the same kinetic energy) and define the feasible set of minimization
problem (41):

$\mathcal{F}_{a}:=\\{\boldsymbol{a}\in\mathbb{R}^{r_{i+1}-d}\hbox{ }|\hbox{
}\widetilde{\Phi}_{i+1}\boldsymbol{a}\in\mathcal{S}(||\widetilde{\boldsymbol{a}_{i}}||)\\}$

and the set
$\mathcal{F}:=\text{span}(\widetilde{\Phi}_{i+1})\cap\mathcal{S}(||\widetilde{\boldsymbol{a}_{i}}||)=\\{\widetilde{\Phi}_{i+1}\boldsymbol{a}\hbox{
}|\hbox{ }\boldsymbol{a}\in\mathcal{F}_{a}\\}$. The solution then corresponds
to an orthogonal projection (in $\mathcal{H}_{\Omega_{h}}$) of
$\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}$ on
$\text{span}(\widetilde{\Phi}_{i+1})$ and a subsequent rescaling to obtain an
element of $\mathcal{F}$. Solution (42) then produces the generalized
coordinates with respect to $\widetilde{\Phi}_{i+1}$ of this scaled
projection. A visual interpretation for $N=3$ and $\Omega_{h}=I$ is provided
in Figure 2.

Using equation (42) and the momentum conservation constraint the generalized
coordinates of $\boldsymbol{u}_{r}^{+}$ as determined by the fully constrained
minimization problem (40) are given by:

$\boxed{\boldsymbol{a}_{i+1}(t_{i+1})^{T}=\left[(\boldsymbol{a}_{i})_{1}\hbox{
},\hbox{ }...\hbox{ },\hbox{ }(\boldsymbol{a}_{i})_{d}\hbox{
},\left(\frac{||\widetilde{\boldsymbol{a}}_{i}||}{||\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}||}\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}\right)^{T}\right].}$

We will refer to using the PID with this structure-preserving interface
condition as SP-PID.

It may happen that $||\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}||=0$, in
this case we cannot scale the result to have the same norm, and thus kinetic
energy, as $\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}$. This
situation occurs when
$\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}\perp_{\Omega_{h}}\text{span}(\widetilde{\Phi}_{i+1})$
and it can be shown that, if
$\text{dim}(\text{span}(\widetilde{\Phi}_{i+1})),N>1$, minimization problem
(41) has an infinite number of solutions. Namely, writing out the loss
function of minimization problem (41), taking into account that
$\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}\perp_{\Omega_{h}}\text{span}(\widetilde{\Phi}_{i+1})$,
we have
$||\widetilde{\Phi}_{i}\widetilde{\boldsymbol{a}}_{i}-\widetilde{\Phi}_{i+1}\boldsymbol{a}||_{\Omega_{h}}=\sqrt{2}||\widetilde{\boldsymbol{a}}_{i}||\hbox{
}\forall\hbox{ }\boldsymbol{a}\in\mathcal{F}_{a}$ and since
$|\mathcal{F}_{a}|$ is infinite in this case we have an infinite number of
solutions. A visual representation of this case for $N=3$ and $\Omega_{h}=I$
is given in Figure 2. In principle any element of $\mathcal{F}_{a}$ can be
used to determine $\boldsymbol{u}_{r}^{+}$ in this situation. We note that in
general $||\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}||=0$ is an unlikely
scenario since the snapshots at the end of $X_{i}$ tend to resemble those at
the beginning of $X_{i+1}$. As a result, $\mathcal{V}^{i}$ and
$\mathcal{V}^{i+1}$ are typically capable of resolving snapshots close to
either side of their shared interface. However, to make sure that we prevent
situations where $||\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}||=0$ (or
close to zero) we suggest to use overlapping snapshots. This entails
concatenating $\gamma_{l}\in\mathbb{N}$ of the rightmost snapshots in
$X_{i-1},\Xi_{i-1}$ and $\gamma_{r}\in\mathbb{N}$ of the leftmost snapshots in
$X_{i+1},\Xi_{i+1}$ with $X_{i}$ and $\Xi_{i}$ respectively for the
construction of $\Phi_{i}$ and $M_{i}$, as illustrated in Figure 3. This will
promote non-empty and larger intersections between the reduced spaces and as a
result discourage $||\widetilde{T}_{i}\widetilde{\boldsymbol{a}}_{i}||=0$,
because both reduced spaces $\mathcal{V}^{i}$ and $\mathcal{V}^{i+1}$ are
constructed to resolve $\boldsymbol{u}_{r}$ on either side of the interface.
In addition, the overlapping of snapshot sets tends to improve the accuracy of
the transition, as will be shown in the results section.

Figure 1: Graphical intuition for optimization solution with $N=3$ and
$\Omega_{h}=I$.

Figure 2: Graphical intuition for infinite number of optimization solutions
with $N=3$ and $\Omega_{h}=I$

Figure 3: Constructing $\mathcal{V}^{i}$ using overlapping snapshots.

##  5 Results

In this section we will discuss the results of two test cases to be introduced
in what follows. The source code for these experiments is a custom finite
volume incompressible flow solver that has been written in the C++ programming
language using the Armadillo linear algebra library [75, 76] and the LIS
library of iterative solvers for linear systems [59]. All experiments in this
paper have been carried out on a simple workstation with a single Intel(R)
Core(TM) i7-6700HQ CPU.

###  5.1 Shear layer roll-up

####  5.1.1 Problem set-up

The first test case concerns the shear layer roll-up (SLR) [54, 44], which is
a flow on a double-periodic domain $[0,2\pi]\times[0,2\pi]$ with a central
band of flow in a positive coordinate direction and neighbouring bands of flow
in the opposite direction. The bands are joined with a thin region of strong
velocity gradients and hence strong shear forces. The initial conditions are:

$u(x,y,0)=\begin{cases}\text{tanh}\left(\frac{y-\pi/2}{\delta}\right),&y\leq\pi\\\
\text{tanh}\left(\frac{3\pi/2-y}{\delta}\right),&y>\pi\end{cases},\quad
v(x,y,0)=\epsilon\sin(x),$

where $\delta=\frac{\pi}{15}$ determines the initial thickness of the shear
layers and the parameter $\epsilon=0.05$ determines the initial amplitude of
an unstable perturbation in the second coordinate direction to trigger the so-
called roll-up.

Using the SLR the proposed hROMs will be tested for their structure-preserving
capabilities, accuracy and computational performance. Structure preservation
will be tested by considering the temporal evolution of errors in conserved
quantities (momentum and kinetic energy) for the first four seconds of the
inviscid SLR flow. Accuracy will be tested by considering two different
errors, namely the error with respect to the FOM, which is defined as

$\epsilon_{u}(t)=\left|\left|\boldsymbol{u}_{h}(t)-\boldsymbol{u}_{r}(t)\right|\right|_{\Omega_{h}},$
(43)

and the best-approximation error

$\epsilon_{b}(t)=\left|\left|(I-\Phi\Phi^{T}\Omega_{h})\boldsymbol{u}_{h}(t)\right|\right|_{\Omega_{h}},$
(44)

which forms a lower bound to the accuracy that can be obtained using $r$ POD
modes. Comparing $\epsilon_{u}$ and $\epsilon_{b}$ will give an indication on
how close to optimal the hROM is. We will also analyse whether overfitting
takes place by considering the temporal evolution of individual DEIM
coordinates and their Mahalanobis distance $d_{M}(\boldsymbol{c})$. Finally,
the computational performance will be tested by measuring execution times
associated to the offline and online phases of the hROMs and comparing these
to the FOM.

Based on a grid convergence study it was concluded that a $256\times 256$ grid
and time step $\Delta t=0.01$ were sufficiently fine to have a reliable FOM
reference result, and these settings will be used in the numerical experiments
unless otherwise mentioned. The same time step is also used to build the
snapshot matrix and to simulate the ROMs. The hyper-parameter of the MCLSDEIM
will be set to $\alpha=1\cdot 10^{-14}$ and the hyper-parameter of the
OCLSDEIM will be set to $m_{p}=2m$.

Figure 4: The normalized POD and DEIM singular value decay of snapshot data
from a FOM simulation with $\nu=0$ and $t=4$.

####  5.1.2 Conservation properties

The proposed hyper-reduction methods will be tested for their structure-
preserving properties. We expect $K_{r}$ to be conserved for the inviscid case
only when using the proposed hROMs in conjunction with energy-conserving
Runge-Kutta (RK) methods (see remark 1). Thus, the proposed methods and the
conventional DEIM will be tested for the case $\nu=0$ using the (implicit)
energy-conserving Gauss-Legendre 4 (GL4) method [72], and compared also to the
classical (explicit) Runge-Kutta 4 (RK4) method [74]. Time integration will
take place until $t=4$, since for longer time intervals numerical oscillations
develop in the inviscid FOM simulation. These numerical oscillations do not
destabilize the FOM solution due to the nonlinear stability property, but they
render the simulation results inaccurate and corrupt the snapshot data sets
and hence the ROM. Until $t=4$ the FOM solution has a quickly decaying
Kolmogorov N-width, and 8 POD and 8 DEIM modes suffice to accurately capture
most of the energy in the snapshots as can be seen in Figure 4. Note that the
POD and DEIM modes are obtained from snapshots taken at every individual time
step of the FOM simulation. Since the initial conditions of the SLR have a net
total momentum of zero in all coordinate direction we will plot the absolute
value of its components $|(\boldsymbol{P}_{r})_{i}|$. We expect this value to
remain zero up to machine precision for all hROMs. For a momentum-conserving
semi-discretization any RK method trivially conserves momentum, hence we only
use RK4 (and not GL4) in order to not clutter the results. For the kinetic
energy we will plot the value:

$\epsilon_{K}(t):=|K_{r}(t)-K_{r}(0)|,$

which denotes the absolute deviation of $K_{r}$ with respect to its initial
value. As kinetic energy is conserved for $\nu=0$ by the structure-preserving
hROMs, this value is also expected to remain zero up to machine precision
using energy-conserving RK methods.

Figure 5: The error $\epsilon_{K}$ in kinetic energy with respect to its
initial value for the inviscid SLR using the proposed hROMs and the DEIM-hROM,
the solid lines are for RK4 and the dashed lines for GL4.

Figure 6: Conservation of momentum component $(\boldsymbol{P}_{r})_{i}$ in
absolute value for the inviscid SLR using the proposed hROMs and the DEIM-hROM
and RK4, the solid lines are for $i=1$ and the dashed lines are for $i=2$.

In Figure 6 the error in the kinetic energy is displayed. Colors are used to
distinguish the different hROM methods and solid lines and dashed lines are
used to distinguish RK4 from GL4 time integration. The key observation is that
when using the standard DEIM algorithm (with either RK4 or GL4), $K_{r}$
deviates significantly from its initial value, while our proposed structure-
preserving hROMs exactly conserve kinetic energy up to machine precision
(using GL4 time integration). When using the RK4 time integrator the
structure-preserving methods deviate slightly from their initial value, due to
the non-energy-conserving nature of RK4. These results are in line with [73]
in which similar small energy errors were observed when using high-order
explicit RK methods to integrate systems with structure-preserving semi-
discretizations.

In Figure 6 the error in the momentum components is displayed. In this figure
colors are used to distinguish which hROM method is applied and solid and
dashed lines are used to indicate the first and second component of
$\boldsymbol{P}_{r}$, respectively. It can be observed that the theoretical
results on momentum conservation hold for all hROMs including the DEIM-hROM.
The deviation from zero of all momentum components
$|(\boldsymbol{P}_{r})_{i}|$ can be seen to increase minutely. We assume this
can be attributed to the finite precision of the SVD implementation of the
Armadillo C++ library, which causes the modes in the DEIM basis $M$ to be
momentum-conserving only up to the SVD precision.

####  5.1.3 Convergence behaviour

To assess the behaviour of $\epsilon_{u}$ as $r$ and $m$ are varied we perform
a convergence study. This will be done for the case of viscous flow, with
$\nu=0.001$, allowing integration of the SLR flow until $t=8$. Given the small
differences between RK4 and GL4 in the study of conservation properties, we
only consider RK4 here for efficiency. Our primary interest is the behaviour
of $\epsilon_{u}$ as a function of $m$, hence we will sweep this parameter for
some values of $r$. To determine what values of $m$ and $r$ are necessary such
that the corresponding reduced spaces $\mathcal{V}$ and $\mathcal{M}_{d}$ can
accurately resolve the snapshot data we study again the singular value decay
of $X$ and $\Xi$. The singular value decays of the snapshot data for a FOM
simulation with $\nu=0.001$ and until $t=8$ have been plotted in Figure 7 for
the data sets $X$ and $\Xi$. From this figure we determine a sweep will be
performed over all values $m\in\\{1,2,...,60\\}$ corresponding to a relative
information content (RIC) [46] ranging from $89.86\%$ to $100\%$. The reduced
space dimension will be taken as $r\in\\{30,40\\}$ both corresponding to
$\text{RIC}\approx 100\%$.

Figure 7: The normalized POD and DEIM singular value decay of snapshot data
from a FOM simulation with $\nu=0.001$ and $t=8$.

The results of the convergence study have been plotted in Figure 9 and Figure
9 for $r=30$ and $r=40$ respectively. In these figures we also plot
$\epsilon_{u}$ for the non-hyper-reduced order model (dashed lines). This
value will be referred to as the ‘ROM error’, obtained by using the exact
tensor decomposition for the convective terms. Clearly, the ROM error is
independent of the dimension $m$ of $\mathcal{M}_{d}$. Although the ROM error
is not strictly a lower bound for $\epsilon_{u}$ of hROMs, the methods should
converge to it if the convection operator is approximated accurately enough.
All methods show erratic error behaviour as $m$ increases for $m<r$. This
shows that the convection operator cannot be expected to be approximated well
throughout a long time-integration procedure when $m$ is not sufficiently
large compared to $r$. When $m\geq r$ better accuracy is obtained with the
OCLSDEIM and MCLSDEIM, and the ROM error is closely approximated. This means
that the convection operator is approximated well enough throughout the time
integration such that $\Phi^{T}C_{h}(\Phi\boldsymbol{a})$ and
$\Phi^{T}M\boldsymbol{c}$ are nearly indistinguishable. It can be seen that
$\epsilon_{u}$ as obtained using the DEIM and CLSDEIM remains behaving
erratically when $m\geq r$ and generally the ROM error is not obtained. This
shows that exactly (DEIM) or near exactly (CLSDEIM) fitting the approximation
$M\boldsymbol{c}$ to $C_{h}(\Phi\boldsymbol{a})$ in $\mathcal{P}$ when
$m_{p}=m$ is not sufficient to guarantee accuracy in a long time integration
process. Indeed, using more measurements (OCLSDEIM) or smarter choices of DEIM
coordinates (MCLSDEIM) is required to obtain the ROM error.

Figure 8: Convergence study of $\epsilon_{u}$ evaluated at $t=8$ for different
hROMs and different values of $m$ for $r=30$, including the non-hyper-reduced
model results (dashed).

Figure 9: Convergence study of $\epsilon_{u}$ evaluated at $t=8$ for different
hROMs and different values of $m$ for $r=40$, including the non-hyper-reduced
model results (dashed).

Figure 10: The velocity and vorticity fields as predicted by the different
hyper-reduction methods at $r=30,m=50,t=8$ ($t=5$ for the unstable DEIM).

####  5.1.4 Temporal error evolution

To obtain a better understanding of the performance of the methods we will
analyse the full temporal evolution of several errors. We will examine a
particularly bad case for the DEIM and CLSDEIM such that the errors are
further exaggerated, specifically we will look at $r=30,m=50$, for which
Figure 9 shows high errors for both DEIM and CLSDEIM. For illustration, the
velocity and vorticity fields as obtained using $r=30,m=50$ are shown in
Figure 10. For this case the error due to the DEIM becomes unstable and due to
the CLSDEIM becomes very large. Besides $\epsilon_{u}$, we will also analyse
the temporal evolution of the operator error:

$\epsilon_{M}:=\frac{||C_{h}(\Phi\boldsymbol{a})-M\boldsymbol{c}(\boldsymbol{a})||}{||C_{h}(\Phi\boldsymbol{a})||}.$

This error provides insight in how well a hyper-reduction method is
approximating the convection operator.

In Figure 11 the error evolution of both $\epsilon_{u}$ and $\epsilon_{M}$
have been displayed for $r=30,m=50$. The ideal projection error $\epsilon_{b}$
has also been plotted as function of time. We once more use colors to
distinguish different methods, solid lines are associated to the evolution of
$\epsilon_{u}$ using the RK4 integrator and dashed lines to $\epsilon_{M}$. We
can see that the errors due to DEIM and CLSDEIM grow gradually as time
progresses. Furthermore, as the solution becomes more dominated by errors it
becomes more difficult for the DEIM and CLSDEIM to accurately approximate the
correct convection operator output. The OCLSDEIM and MCLSDEIM are much more
accurate throughout the full time integration process. It can be seen that
their respective $\epsilon_{u}$ values stay close to the lower bound
$\epsilon_{b}$ and that the operator output is approximated very accurately.

Figure 11: The evolution of $\epsilon_{u}$ and $\epsilon_{M}$ (dashed lines)
for $r=30$ and $m=50$ using different hROMs and for both RK4 (solid lines)
compared to the ideal projection error $\epsilon_{b}$.

Another experiment that gives insight in the error is to plot the temporal
evolution of individual DEIM coordinates. We have carried out this experiment
again for the ‘bad’ case $r=30,m=50$. The results are displayed in Figure 12,
where the evolution of a subset of the first 26 DEIM coordinates is given.
Colors are used to indicate different DEIM coordinates and different linetypes
are used to indicate what hyper-reduction method was used. We also plot the
reference values for the DEIM coordinates as given by the columns of
$M^{T}\Xi$. In this reproduction case the DEIM coordinates as given by the
different hyper-reduction methods should be as close as possible to these
ideal reference values. We can see the DEIM coordinates as determined by the
DEIM gradually becoming unstable around 5 seconds into the simulation. The
DEIM coordinates as determined by the CLSDEIM differ significantly from the
reference values. The OCLSDEIM and MCLSDEIM perform much better and almost
coincide with the reference values throughout the full simulation. Plotting
the Mahalanobis distance $d_{M}(\boldsymbol{c})$ in Figure 13, the poor
accuracy of the DEIM coordinates as determined by the DEIM and CLSDEIM is
accentuated. As time progresses the values of the DEIM coordinates becomes
increasingly less similar to the reference values in $M^{T}\Xi$, whereas the
DEIM coordinates as determined by the OCLSDEIM and MCLSDEIM remain close to
the reference distribution. Furthermore, Figure 13 shows that the Mahalanobis
distance, at least in the reproduction case, is a good indicator of a poor
choice in DEIM coordinates as it becomes increasingly larger when
$\epsilon_{M}$ becomes larger. Combining these results and the previous
results on error evolution, we see that by exactly (DEIM) or near exactly
(CLSDEIM) fitting the approximation $M\boldsymbol{c}$ to
$C_{h}(\Phi\boldsymbol{a})$ in $\mathcal{P}$ with $m_{p}=m$, DEIM coordinates
are determined that cause gradual accumulation of errors in the time
integration process. Instead, more robustness is required to deal with
demanding situation like time integration where accumulation of errors can
cause problems. As shown in our experiments, robustness of the DEIM and
CLSDEIM can be significantly improved by:

  1. 1.

performing a regression through a larger set of measurements than free
parameters (OCLSDEIM);

  2. 2.

using knowledge on the prior distribution of DEIM coordinates (MCLSDEIM).

Figure 12: Temporal evolution of DEIM coordinates $c_{i}$,
$i\in(1,2,11,16,26)$ using $r=30,m=50$ compared to the reference values as
given by the columns of $M^{T}\Xi$. Figure 13: Temporal evolution of the
Mahalanobis distance $d_{M}(\boldsymbol{c})$ of the DEIM coordinates using
$r=30,m=50$ for the different hyper-reduction methods.

####  5.1.5 Performance

To end our discussion on the SLR experiment we will show that the hyper-
reduction methods exhibit similar computational performance as the
conventional DEIM. This will be done by measuring the on- and offline phase
wall-clock times and the associated speed-ups with respect to the FOM. Here,
the offline phase is constituted by the necessary SVDs to construct
$\mathcal{M}_{d}$ and $\mathcal{V}$, the construction of $\mathcal{P}$ and the
precomputation of the necessary matrices and LU-decompositions. The online
phase is constituted by the numerical integration of the hROMs. We will
evaluate this for $r=30,m=40$ as all methods provide accurate results for this
configuration of reduced space dimensions. To see the effect of increasing
$m_{p}$ for the OCLSDEIM we will test for $m_{p}\in\\{2m,4m\\}$, where we
expect to see a decline in performance as $m_{p}$ is increased.

The results are tabulated in Table 1 for the online and offline phase. It can
be seen that the computational performance is similar for all methods, both in
terms of online and offline performance. The DEIM, CLSDEIM and MCLSDEIM are
especially close due to the fact that for all these cases $m_{p}=m$. The
effect of the larger number of measurement evaluations required can be seen
for the OCLSDEIM(2m) and OCLSDEIM(4m), which are slower than the other
methods. Furthermore, using $m_{p}=2m$ of the OCLSDEIM indeed provides better
computational performance than $m_{p}=4m$, as expected.

Method | Wall-clock time online (s) | Speed-up | Wall-clock time offline (s)
---|---|---|---
DEIM | 0.236 | 808 | 45.105
CLSDEIM | 0.300 | 636 | 45.104
OCLSDEIM(2m) | 0.494 | 386 | 44.972
OCLSDEIM(4m) | 1.099 | 174 | 45.123
MCLSDEIM | 0.302 | 632 | 44.215
Table 1: Computational performance (offline and online) of the different
hyper-reduction methods for $r=30,m=40$.

###  5.2 Two-dimensional homogeneous isotropic turbulence

####  5.2.1 Problem set-up

The second test case is two-dimensional homogeneous isotropic turbulence
(2DT). 2DT is a flow on a double-periodic domain $[0,2\pi]\times[0,2\pi]$. We
will generate initial conditions following the procedure outlined in [70]. In
short, we construct an initial energy spectrum:

$E(k)=\frac{a_{s}}{2}\frac{1}{k_{p}}\left(\frac{k}{k_{p}}\right)^{2s+1}\exp{\left[-\left(s+\frac{1}{2}\right)\left(\frac{k}{k_{p}}\right)^{2}\right]},$

where $k:=||\boldsymbol{k}||=\sqrt{k_{x}^{2}+k_{y}^{2}}$ is the wave-vector
magnitude. The maximum value of the initial spectrum takes place at
$k_{p}=12$, $a_{s}\in\mathbb{R}$ is a normalization factor given by
$a_{s}:=(2s+1)^{s+1}/(2^{s}s!)$ and $s\in\mathbb{R}$ is a shape parameter
taken as $s=3$. From this spectrum we follow [70] and calculate vorticity
Fourier-coefficients like:

$|\hat{\omega}(\boldsymbol{k})|=\sqrt{\frac{k}{\pi}E(k)}.$

Consequently, we generate a vorticity spectrum
$\hat{\omega}(\boldsymbol{k}):\mathbb{R}^{2}\rightarrow\mathbb{C}$ by
providing random phases to the Fourier-coefficients using a process for which
we refer to [70]. Finally, we obtain a divergence-free velocity field by
solving the relation:

$\Delta\psi=-\omega,\quad\boldsymbol{u}=\begin{bmatrix}\frac{\partial\psi}{\partial
y},-\frac{\partial\psi}{\partial x}\end{bmatrix}^{T},$

numerically. Note that this procedure is strictly for two-dimensional flows.

2DT is a fluid flow that has been studied extensively [3, 27, 70].
Furthermore, it is a highly convection-dominated flow with many convected
spatial structures and as a result it has a slowly decaying Kolmogorov
N-width. For this reason the purpose of this test case is two-fold. Namely,
the complexity of the flow makes it an adequate test case to see how well our
proposed hyper-reduction methods will perform for more turbulent flows of
engineering interest. Additionally, the slow Kolmogorov N-width decay will
allow us to compare the PID and our newly proposed structure-preserving
alternative, the SP-PID, in a meaningful setting.

The FOM simulations are performed with a $1024\times 1024$ grid and a time
step of $\Delta t=0.0002$. These settings are based on the observation that
these also led to converged results for similar spatially second-order energy-
conserving schemes in [70]. Saving snapshots every time step for this test
case is not feasible given the simple workstation used to perform all
calculations. Hence, the simulation will be sampled such that the Nyquist-
Shannon criterion is met [3]; this criterion is satisfied by sampling every
$\Delta t_{s}=0.005$.

####  5.2.2 Hyper-reduction methods

Due to the previous results for the SLR flow, in this section we will only
test the two best-performing methods, MCLSDEIM and the OCLSDEIM. We simulate
the FOM until $t=4$ with $\nu=0.001$ and compare how well the resulting energy
spectra $E(k)$ of the FOM and hROMs align at $t\in\\{1,2,3,4\\}$. We will also
provide a visual comparison of the vorticity fields at these time instances.
To deal with the slow Kolmogorov N-width decay we will apply the SP-PID with
$n=8$ evenly sized intervals of 100 snapshots and for all intervals we will
take an overlap parameters $\gamma_{l}=\gamma_{r}=\gamma=20$. The singular
values of all intervals and the full snapshot sets have been displayed in
Figure 15 and Figure 15 for the POD and DEIM, respectively. As discussed in
[3], it can be seen that the application of temporal localization indeed leads
to a faster decay of singular values within intervals. This effectively
increases the Kolmogorov N-width decay. For the reduced space dimensions we
will use $r=30$ and $m=40$ for all intervals. Like [70], we will calculate the
energy spectra using the following equations:

$E(k,t)=\sum_{k\leq||\widetilde{\boldsymbol{k}}||<{k+1}}\hat{E}\left(\widetilde{\boldsymbol{k}},t\right),\quad\hat{E}(\boldsymbol{k},t)=\frac{1}{2}k^{2}|\hat{\psi}(\boldsymbol{k},t)|^{2},$

where
$\hat{\psi}(\boldsymbol{k},t):\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{C}$
are the Fourier-coefficients of the numerically calculated stream-function
associated to a discrete velocity field $\boldsymbol{u}_{h}(t)$.

Figure 14: Singular value decay of solution snapshots of intervals versus the
full snapshot set.

Figure 15: Singular value decay of operator snapshots of intervals versus the
full snapshot set.

We have provided the energy spectra for the different hyper-reduction methods
in Figure 17 and Figure 17. Excellent agreement at all points in time can be
observed for both the hyper-reduction methods. Furthermore, the results are in
line with those observed in [70] for $Re=1000$ in terms of spectra, where an
accurate high-order pseudospectral method was used. We can see clearly that
the spectrum is shifting to the left as time progresses, implying the presence
of an inverse energy cascade. Furthermore, the inertial range nearly scales as
$\mathcal{O}(k^{-3})$ which should be the case as $\nu\rightarrow 0$ [27]. We
have also shown the vorticity fields as calculated using the FOM and the
different hROMs in Figure 18. We can see that, besides the energy spectra
matching between the FOM and hROMs, the vorticity fields are also almost
identical. From this we conclude that the hROMs using the MCLSDEIM and
OCLSDEIM and the SP-PID:

  1. 1.

can accurately reconstruct the FOM results;

  2. 2.

adhere to the physics underlying the 2DT flow.

Figure 16: Energy spectra at $t\in\\{1,2,3,4\\}$ as calculated by the FOM
(solid line) and MCLSDEIM-SP-PID hROM (dashed line).

Figure 17: Energy spectra at $t\in\\{1,2,3,4\\}$ as calculated by the FOM
(solid line) and OCLSDEIM-SP-PID hROM (dashed line).

Figure 18: 2DT vorticity fields as calculated using the FOM and hROMs for
$t\in\\{1,2,3,4\\}$.

####  5.2.3 Temporal localization methods

We will now study in more detail the conservation properties of SP-PID at
interfaces and compare its transition accuracy against the orthogonal
projection of the standard PID as in equation (38). Furthermore, we will
briefly investigate the impact of the overlap parameter $\gamma$. We will do
this by simulating the FOM until $t=1$ with $\nu=0.001$ and only considering
$n=2$ evenly sized intervals. We will then construct the reduced spaces
associated to these intervals with overlap parameters $\gamma\in\\{0,20\\}$
and analyse the differences this makes on the accuracy of the transition at
$t=0.5$. We will also consider the accuracy of the conserved quantities energy
and momentum, which should be exactly conserved by the SP-PID and not by the
PID. However, for the PID, we expect these quantities to be better conserved
for larger $\gamma$. We also expect the accuracy of the transition to be
improved by taking larger $\gamma$, for both SP-PID and PID. During this
experiment we will use the MCLSDEIM with $\alpha=1\cdot 10^{-14}$ and $m=40$.
Additionally, we will use a reduced space of $r=30$. Finally, we will
determine the accuracy of the transition using the transition error:

$\epsilon_{T}:=||\boldsymbol{u}_{r}^{+}-\boldsymbol{u}_{r}^{-}||_{\Omega_{h}}.$

| SP-PID ($\gamma=0$) | PID ($\gamma=0$) | SP-PID ($\gamma=20$) | PID ($\gamma=20$)
---|---|---|---|---
$\epsilon_{T}$ | 0.0629 | 0.0629 | 0.00515 | 0.00515
$|K_{r}^{+}-K_{r}^{-}|$ | 0 | 0.00198 | 0 | $1.33\cdot 10^{-5}$
$|(\boldsymbol{P}_{r})_{1}^{+}-(\boldsymbol{P}_{r})_{1}^{-}|$ | 0 | $2.83\cdot 10^{-14}$ | 0 | $3.32\cdot 10^{-15}$
$|(\boldsymbol{P}_{r})_{2}^{+}-(\boldsymbol{P}_{r})_{2}^{-}|$ | 0 | $8.49\cdot 10^{-17}$ | 0 | $2.16\cdot 10^{-16}$
Table 2: Transition errors in velocity field and conserved quantities at the
interface on $t=0.5$.

The results of the experiment have been displayed in Table 2. It can be seen
that the transition accuracy indeed increases as the overlap parameter is
increased, since the transition error $\epsilon_{T}$ decreases significantly.
Furthermore, $\epsilon_{T}$ using the SP-PID and the PID is nearly identical
for both overlap parameter values. This shows that even in the unconstrained
setting the SP-PID interface condition (40) is nearly optimal. As expected,
the conserved quantities are exactly conserved over interfaces using the SP-
PID. Increasingly smaller errors are made using the standard PID as the
overlap parameter increases. We have also shown the energy spectra of the
velocity fields before ($\boldsymbol{u}_{r}^{-}$) and after
($\boldsymbol{u}_{r}^{+}$) the interface in Figure 20 and Figure 20 for
$\gamma=0$ and $\gamma=20$, respectively. We see good agreement between the
new and old spectra for both methods, however for larger overlap parameters
the smaller scales in the flow are captured more accurately. This can be seen
in the zoomed Figure 22 and Figure 22 for $\gamma=0$ and $\gamma=20$,
respectively. From this we conclude that:

  1. 1.

the SP-PID exactly conserves both energy and momentum over interfaces, whereas
the PID does not;

  2. 2.

the overall accuracy of the transition (conserved quantities and
$\epsilon_{T}$) improves with greater overlap parameters $\gamma$ for the SP-
PID and PID.

We suspect there will be an optimal overlap parameter, this will be the
subject of future research.

Figure 19: Spectra before and after transition at $t=0.5$ using PID and SP-PID
for $\gamma=0$.

Figure 20: Spectra before and after transition at $t=0.5$ using PID and SP-PID
for $\gamma=20$.

Figure 21: Spectra before and after transition at $t=0.5$ using PID and SP-PID
for $\gamma=0$ (zoomed).

Figure 22: Spectra before and after transition at $t=0.5$ using PID and SP-PID
for $\gamma=20$ (zoomed).

##  6 Conclusion

In this article we have proposed new structure-preserving hyper-reduction
methods referred to as the CLSDEIM, OCLSDEIM and MCLSDEIM. In these new hyper-
reduction methods the conventional DEIM is posed as a minimization problem and
subsequently constrained to conserve kinetic energy, in addition to conserving
momentum. The OCLSDEIM and MCLSDEIM enhance the robustness of the CLSDEIM by
oversampling and by Mahalanobis regularization, respectively.

The second novelty in this work is that we have proposed a method to conserve
energy and momentum over interfaces in order to use the new hyper-reduction
methods in combination with the PID temporal localization method. The method
is based on solving a constrained minimization problem at every interface, and
its solution can be interpreted as a projection followed by a scaling to
conserve energy. This new SP-PID allows the construction of fully structure-
preserving and temporally localized hROMs of convection-dominated flows.

We have performed two numerical experiments to show the performance of the new
structure-preserving hyper-reduction methods and the structure-preserving
temporal localization method. We used the shear layer roll-up (SLR) to test
the structure-preserving properties, accuracy and computational performance of
the proposed hyper-reduction methods. We confirmed the theoretical results
that the CLSDEIM, OCLSDEIM and MCLSDEIM conserve energy and momentum in a
fully discrete setting when the Gauss-Legendre 4 (GL4) Runge-Kutta (RK) time
integration method is used. It was also shown that the standard DEIM caused a
significant error in kinetic energy in this case. Convergence experiments
showed that both the DEIM and CLSDEIM exhibit erratic behaviour as a function
of the DEIM space dimension $m$, due to error accumulation during the time
integration process. This erratic behaviour and error accumulation is
alleviated by the robust OCLSDEIM and MCLSDEIM. Finally, the computational
performance of the DEIM, CLSDEIM and MCLSDEIM was comparable, whereas the
OCLSDEIM was slower than the other methods. This shows the benefit of using
regularization instead of oversampling to increase robustness of the CLSDEIM.

The second numerical experiment was the two-dimensional homogeneous isotropic
turbulence (2DT). The results of this test case further highlighted the
accuracy and robustness of the new hyper-reduction methods. Furthermore, the
2DT experiment showed that using the SP-PID energy and momentum can be
conserved over interfaces, in contrast to the classical PID. Finally, it was
shown that overlapping snapshots in the construction of local reduced spaces
can increase the transition accuracy.

There are several directions of future research that we believe to be valuable
or interesting. First of all is the generalization of the MCLSDEIM and SP-PID
to parametric problems. For the MCLSDEIM this entails confirming that the DEIM
coordinates of a parametric problem can indeed be thought of as originating
from a single prior distribution. For the SP-PID this entails finding methods
to efficiently and accurately distribute the full parametric snapshot set over
multiple subsets. Besides these methods, the extension to arbitrary boundary
conditions (instead of periodic conditions) will be very valuable. Finally,
the methods proposed in this article could be merged with ideas from e.g. [22,
83] to obtain structure-preserving hROMs for problems with general first
integrals and general nonlinear skew-symmetric operators (like the convection
operator in this article or those in [55]). We believe this could be done by
suitably altering condition (26) to conserve a general hyper-reduced first
integral.

## CRediT authorship contribution statement

R.B. Klein: Conceptualization, Methodology, Software, Formal analysis,
Investigation, Writing - Original Draft, Visualization
B. Sanderse: Conceptualization, Writing - Review & Editing, Supervision,
Project administration

## 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

The authors thank R.A.W.M. Henkes for his supervision during the master thesis
behind this work and thank C. Pagliantini for the stimulating discussions on
discrete empirical interpolation during the MORe conference at TU Berlin and
while visiting the CWI institute. This publication is part of the project
“Discretize first, reduce next” (with project number VI.Vidi.193.105 of the
research programme NWO Talent Programme Vidi which is (partly) financed by the
Dutch Research Council (NWO).

## References

  * [1] B.M. Afkham and J.S. Hesthaven.  Structure-preserving model-reduction of dissipative hamiltonian systems.  Journal of Scientific Computing, 81(1):3–21, 2019.
  * [2] B.M. Afkham, N. Ripamonti, Q. Wang, and J.S. Hesthaven.  Conservative model order reduction for fluid flow, pages 67–99.  Springer International Publishing, Cham, 2020.
  * [3] S.E. Ahmed and O. San.  Breaking the kolmogorov barrier in model reduction of fluid flows.  Fluids, 5(1), 2020.
  * [4] N. Akkari, F. Casenave, and V. Moureau.  Time stable reduced order modeling by an enhanced reduced order basis of the turbulent and incompressible 3d navier–stokes equations.  Mathematical and Computational Applications, 24:45, 04 2019.
  * [5] D. Amsallem and C. Farhat.  Stabilization of projection-based reduced-order models.  International Journal for Numerical Methods in Engineering, 91(4):358 – 377, 2012.
  * [6] V.I. Arnol’d.  The hamiltonian nature of the euler equations in the dynamics of a rigid body and of an ideal fluid.  1969\.
  * [7] M. Balajewicz and E.H. Dowell.  Stabilization of projection-based reduced order models of the navier-stokes.  Nonlinear Dynamics, 70(2):1619 – 1632, 2012.
  * [8] M. Balajewicz, E.H. Dowell, and B. Noack.  Low-dimensional modelling of high reynolds number shear flows incorporating constraints from the navier-stokes equation.  Journal of Fluid Mechanics, 729:285–308, 08 2013.
  * [9] 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.
  * [10] M. Barrault, Y. Maday, N.C. Nguyen, and A.T. Patera.  An ‘empirical interpolation’ method: Application to efficient reduced-basis discretization of partial differential equations.  Comptes Rendus Mathematique, 339(9):667–672, 2004.
  * [11] P. Benner, P. Goyal, J. Heiland, and I. Pontes Duff.  Operator inference and physics-informed learning of low-dimensional models for incompressible flows.  ArXiv, 2020.
  * [12] P. Blonigan, K. Carlberg, F. Rizzi, M. Howard, and J. Fike.  Model reduction for hypersonic aerodynamics via conservative lspg projection and hyper-reduction.  AIAA Paper 2020-0104, AIAA Scitech 2020 Forum, Orlando, FL, January 6–10, 2020.
  * [13] 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.
  * [14] S. Boyd, S.P. Boyd, and L. Vandenberghe.  Convex Optimization.  Number pt. 1 in Berichte über verteilte messysteme. Cambridge University Press, 2004.
  * [15] S.L. Brunton and J.N. Kutz.  Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press, 2022.
  * [16] K. 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.
  * [17] K. 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:155 – 181, 04 2011.
  * [18] K. Carlberg, Y. Choi, and S. Sargsyan.  Conservative model reduction for finite-volume models.  Journal of Computational Physics, 371, 11 2017.
  * [19] K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem.  The gnat method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows.  Journal of Computational Physics, 242:623–647, 2013.
  * [20] J. Chan.  Entropy stable reduced order modeling of nonlinear conservation laws.  Journal of Computational Physics, 423:109789, dec 2020.
  * [21] S. Chaturantabut.  Temporal localized nonlinear model reduction with a priori error estimate.  Applied Numerical Mathematics, 119:225–238, 2017.
  * [22] S. Chaturantabut, C. Beattie, and S. Gugercin.  Structure-preserving model reduction for nonlinear port-hamiltonian systems.  SIAM Journal on Scientific Computing, 38(5):B837–B865, 2016.
  * [23] S. Chaturantabut and D.C. Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  SIAM J. Sci. Comput., 32:2737–2764, 2010.
  * [24] S. Chaturantabut and D.C. Sorensen.  A state space error estimate for pod-deim nonlinear model reduction.  SIAM J. Numer. Anal., 50:46–63, 2012.
  * [25] Saifon Chaturantabut.  Nonlinear Model Reduction via Discrete Empirical Interpolation.  PhD thesis, Rice University, 2011.
  * [26] 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.
  * [27] P.A. Davidson.  Turbulence: An Introduction for Scientists and Engineers.  OUP Oxford, 2004.
  * [28] M. D’Elia, M. Gunzburger, and G. Rozza.  Quantification of Uncertainty: Improving Efficiency and Technology QUIET selected contributions: QUIET selected contributions.  01 2020.
  * [29] C. Farhat, P. Avery, T. Chapman, and J. Cortial.  Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency.  International Journal for Numerical Methods in Engineering, 98, 06 2014.
  * [30] C. Farhat, T. Chapman, and P. Avery.  Structure‐preserving, stability, and accuracy properties of the energy‐conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models.  International Journal for Numerical Methods in Engineering, 102:1077–1110, 2015.
  * [31] 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.
  * [32] 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.
  * [33] G. Golub, M. Heath, and G. Wahba.  Generalized cross-validation as a method for choosing a good ridge parameter.  Technometrics, 21:215–223, 05 1979.
  * [34] Y. Gong, Q. Wang, and Z. Wang.  Structure-preserving galerkin pod reduced-order modeling of hamiltonian systems.  Computer Methods in Applied Mechanics and Engineering, 315:780–798, 2017.
  * [35] 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, 12 2020.
  * [36] S. Grimberg, C. Farhat, and N. Youkilis.  On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows.  J. Comput. Phys., 419:109681, 2020.
  * [37] D. Hartmann, M. Herz, and U. Wever.  Model order reduction a key technology for digital twins, pages 167–179.  04 2018.
  * [38] T. Hastie, R. Tibshirani, and J.H. Friedman.  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer series in statistics. Springer, 2001.
  * [39] J.S. Hesthaven, C. Pagliantini, and N. Ripamonti.  Structure-preserving model order reduction of hamiltonian systems.  ArXiv, 2021.
  * [40] 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.
  * [41] C. Huang, C. Wentland, K. Duraisamy, and C. Merkle.  Model reduction for multi-scale transport problems using structure-preserving least-squares projections with variable transformation.  ArXiv, 11 2020.
  * [42] W.L. IJzerman.  Signal Representation and Modeling of Spatial Structures in Fluids.  PhD thesis, University of Twente, April 2000.
  * [43] I. Kalashnikova and S. Arunajatesan.  A stable galerkin reduced order model for compressible flow.  Blucher Mechanical Engineering Proceedings, 1(1):1399 – 1423, 2014\.
  * [44] R. Knikker.  Study of a staggered fourth‐order compact scheme for unsteady incompressible viscous flows.  International Journal for Numerical Methods in Fluids, 59, 2009\.
  * [45] E. Kreyszig.  Introductory Functional Analysis with Applications.  Wiley Classics Library. Wiley, 1991.
  * [46] 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.
  * [47] K. Lee and K.T. Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [48] J. Loiseau and S.L. Brunton.  Constrained sparse galerkin regression.  Journal of Fluid Mechanics, 838:42–67, 2018.
  * [49] T. Lu and S. Shiou.  Inverses of 2 × 2 block matrices.  Computers & Mathematics with Applications, 43(1):119–129, 2002\.
  * [50] K. Manohar, B.W. Brunton, J.N. Kutz, and S.L. Brunton.  Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns.  IEEE Control Systems Magazine, 38(3):63–86, 2018.
  * [51] R. Maulik, B. Lusch, and P. Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  Physics of Fluids, 33(3):037106, Mar 2021.
  * [52] S.A. McQuarrie, S. Huang, and K.E. Willcox.  Data-driven reduced-order models via regularised operator inference for a single-injector combustion process.  Journal of the Royal Society of New Zealand, 51(2):194–211, 2021\.
  * [53] Barone M.F., Kalashnikova I., Segalman D.J., and Thornquist H.K.  Stable galerkin reduced order models for linearized compressible flow.  Journal of Computational Physics, 228(6):1932 – 1946, 2009.
  * [54] M.L. Minion and D.L. Brown.  Performance of under-resolved two-dimensional incompressible flow simulations, ii.  Journal of Computational Physics, 138(2):734–765, 1997.
  * [55] Y. Miyatake.  Structure-preserving model reduction for dynamical systems with a first integral.  Japan Journal of Industrial and Applied Mathematics, 2019.
  * [56] M. Mohebujjaman, L. Rebholz, and T. Iliescu.  Physically‐constrained data‐driven correction for reduced order modeling of fluid flows.  International Journal for Numerical Methods in Fluids, 89, 09 2018\.
  * [57] M. Mohebujjaman, L.G. Rebholz, X. Xie, and T. Iliescu.  Energy balance and mass conservation in reduced order models of fluid flows.  Journal of Computational Physics, 346:262–277, October 2017.
  * [58] N.T. Mücke, S.M. Bohté, and C.W. Oosterlee.  Reduced order modeling for parameterized time-dependent pdes using spatially and memory aware deep learning.  J. Comput. Sci., 53:101408, 2021.
  * [59] A. Nishida.  Experience in developing an open source scalable software infrastructure in japan.  volume 6017, pages 448–462, 03 2010.
  * [60] B.R. Noack, M. Morzynski, and G. Tadmor.  Reduced-Order Modelling for Flow Control.  CISM International Centre for Mechanical Sciences. Springer Vienna, 2011\.
  * [61] P.J. Olver.  A nonlinear hamiltonian structure for the euler equations.  Journal of Mathematical Analysis and Applications, 1982.
  * [62] B. Peherstorfer, D. Butnaru, K. Willcox, and H. Bungartz.  Localized discrete empirical interpolation method.  SIAM Journal on Scientific Computing, 36(1):A168–A192, 2014.
  * [63] B. Peherstorfer, Z. Drmač, and S. Gugercin.  Stability of discrete empirical interpolation and gappy proper orthogonal decomposition with randomized and deterministic sampling points.  SIAM Journal on Scientific Computing, 42(5):A2837–A2864, 2020.
  * [64] B. Peherstorfer and K. Willcox.  Data-driven operator inference for nonintrusive projection-based model reduction.  Computer Methods in Applied Mechanics and Engineering, 306:196–215, 2016.
  * [65] R. Polyuga and A. Schaft.  Effort and flow-constraint reduction methods for structure preserving model reduction of port-hamiltonian systems.  Systems & Control Letters, 61, 03 2012.
  * [66] F. Romor, G. Stabile, and G. Rozza.  Non-linear manifold rom with convolutional autoencoders and reduced over-collocation method.  ArXiv, 2022.
  * [67] C.W. Rowley, T. Colonius, and R.M. Murray.  Model reduction for compressible flows using pod and galerkin projection.  Physica D: Nonlinear Phenomena, 189(1):115–129, 2004.
  * [68] G. Rozza, M. Hess, G. Stabile, M. Tezzele, and F. Ballarin.  1 Basic ideas and tools for projection-based model reduction of parametric partial differential equations, pages 1–47.  De Gruyter, 2020.
  * [69] A.K. Saibaba.  Randomized discrete empirical interpolation method for nonlinear model reduction.  SIAM Journal on Scientific Computing, 42(3):A1582–A1608, 2020.
  * [70] O. San and A. Staples.  High-order methods for decaying two-dimensional homogeneous isotropic turbulence.  Computers & Fluids, 63, 12 2012.
  * [71] B. Sanderse.  Energy-Conserving Discretization Methods for the Incompressible Navier-Stokes Equations : Application to the Simulation of Wind-Turbine Wakes.  PhD thesis, Centrum voor Wiskunde en Informatica, 2013.
  * [72] B. Sanderse.  Energy-conserving runge–kutta methods for the incompressible navier–stokes equations.  Journal of Computational Physics, 233:100–131, 2013.
  * [73] B. Sanderse.  Non-linearly stable reduced-order models for incompressible flow with energy-conserving finite volume methods.  Journal of Computational Physics, 421:109736, 2020.
  * [74] B. Sanderse and B. Koren.  Accuracy analysis of explicit runge–kutta methods applied to the incompressible navier–stokes equations.  Journal of Computational Physics, 231(8):3041–3063, 2012.
  * [75] C. Sanderson and R. Curtin.  Armadillo: A template-based c++ library for linear algebra.  Journal of Open Source Software, 1:26, 07 2016.
  * [76] C. Sanderson and R. Curtin.  A user-friendly hybrid sparse matrix class in c++.  pages 422–430, 07 2018.
  * [77] A.M. Schein, K.T. Carlberg, and M.J. Zahr.  Preserving general physical properties in model reduction of dynamical systems via constrained-optimization projection.  ArXiv, abs/2011.13998, 2020.
  * [78] H. Sharma, Z. Wang, and B. Kramer.  Hamiltonian operator inference: Physics-preserving learning of reduced-order models for canonical hamiltonian systems.  Physica D: Nonlinear Phenomena, 431:133122, 2022.
  * [79] L. Sirovich.  Turbulence and the dynamics of coherent structures. i - coherent structures. ii - symmetries and transformations. iii - dynamics and scaling.  Quarterly of Applied Mathematics - QUART APPL MATH, 45, 10 1987\.
  * [80] B.E. Sonday, A. Singer, C.W. Gear, and I.G. Kevrekidis.  Manifold learning techniques and model reduction applied to dissipative pdes.  ArXiv, 2010.
  * [81] I. Tezaur and M. Barone.  Stable and efficient galerkin reduced order models for non-linear fluid flow.  6th AIAA Theoretical Fluid Mechanics Conference, 06 2011.
  * [82] A. Treuille, A. Lewis, and Z. Popovic.  Model reduction for real-time fluids.  ACM Trans. Graph., 25:826–834, 07 2006.
  * [83] Z. Wang.  Structure-preserving galerkin pod-deim reduced-order modeling of hamiltonian systems.  ArXiv, abs/2103.00388, 2021.
  * [84] P. Wesseling.  Principles of Computational Fluid Dynamics.  Springer Series in Computational Mathematics. Springer, 2001.
  * [85] X. Xie, M. Mohebujjaman, L. Rebholz, and T. Iliescu.  Data-driven filtered reduced order modeling of fluid flows.  SIAM Journal on Scientific Computing, 40, 09 2017.
  * [86] R. Zimmermann and K. Willcox.  An accelerated greedy missing point estimation procedure.  SIAM Journal on Scientific Computing, 38(5):A2827–A2850, 2016.