# Data-Driven Closure of Projection-Based Reduced Order Models for Unsteady
Compressible Flows

Victor Zucatti  William Wolf

###### Abstract

A data-driven closure modeling based on proper orthogonal decomposition (POD)
temporal modes is used to obtain stable and accurate reduced order models
(ROMs) of unsteady compressible flows. Model reduction is obtained via
Galerkin and Petrov-Galerkin projection of the non-conservative compressible
Navier-Stokes equations. The latter approach is implemented using the least-
squares Petrov-Galerkin (LSPG) technique and the present methodology allows
pre-computation of both Galerkin and LSPG coefficients. Closure is performed
by adding linear and non-linear coefficients to the original ROMs and
minimizing the error with respect to the POD temporal modes. In order to
further reduce the computational cost of the ROMs, an accelerated greedy
missing point estimation (MPE) hyper-reduction method is employed. A canonical
compressible cylinder flow is first analyzed and serves as a benchmark. The
second problem studied consists of the turbulent flow over a plunging airfoil
undergoing deep dynamic stall. For the first case, linear and non-linear
closure coefficients are both low in intrusiveness, capable of providing
results in excellent agreement with the full order model. Regularization of
calibrated models is also straightforward for this case. On the other hand,
the dynamic stall flow is significantly more challenging, specially when only
linear coefficients are used. Results show that non-linear calibration
coefficients outperform their linear counterparts when a POD basis with fewer
modes is used in the reconstruction. However, determining a correct level of
regularization is more complicated with non-linear coefficients. Hyper-reduced
models show good results when combined with non-linear calibration and an
appropriate sized POD basis.

###### keywords:

Reduced order model , Data-driven closure , Proper orthogonal decomposition ,
Galerkin projection , Petrov-Galerkin projection

[inst1]organization=University of Campinas, city=Campinas, postcode=13086-860,
country=Brazil

##  1 Introduction

Complex engineering flows can be accurately solved through numerical
simulations in current supercomputers. However, such flows can easily require
hundreds of millions of degrees of freedom to produce meaningful results [1]
and, consequently, the improvement in computational performance of the past
decades is still far from sufficient to produce real-time solutions to several
problems of interest. On the other hand, the large amount of data currently
being generated by high-fidelity simulations and experiments makes surrogate
modeling a viable alternative for the investigation of complex flows. These
methods are called reduced order models (ROMs) when data mining techniques
such as proper orthogonal decomposition (POD) [2, 3] are employed to reduce
the problem dimensionality. Recently, non-linear order reduction techniques
such as deep convolutional autoencoders [4] and kernel-POD [5] have become a
very active research area. Despite the advances in reduced order modeling, the
construction of accurate and long-term stable ROMs for complex non-linear
problems, such as those found in transport phenomena, still remains a
challenging task.

The implementation of techniques for order reduction is often combined with
physics-based or data-driven models. The first approach is usually performed
by projection of a reduced basis in the governing equations being solved.
Galerkin projection [6] is by far the most common physics-based technique, but
Petrov-Galerkin techniques such as the least-squares Petrov-Galerkin (LSPG)
[7, 8] have also been successfully employed. In the last few years, data-
driven ROMs based on regression [9, 10] have been also proposed and they
usually have no physical constrains from the governing equations. A comparison
of projection-based and data-driven ROMs for unsteady flows involving natural
convection can be found in [11].

The need for closure modeling is inherent to the simulation of unsteady
turbulent flows that are usually solved by formulations such as detached eddy
simulation (DES) and large eddy simulation (LES), which require modeling of
subgrid scales. Reduced order models suffer from similar problems; for
example, numerical errors arise from high-frequency filtering due to POD-basis
truncation. Moreover, numerical integration and derivation errors may arise
from the wavenumber interaction of the non-linear operators appearing in the
Navier-Stokes equations. These sources of errors lead to unstable and
inaccurate ROMs [6, 12, 13]. Closure modeling of ROMs is a very active
research field and, consequently, a plethora of methods have been proposed in
the last twenty years. However, these methods lack application to more
realistic problems involving a broad range of temporal and spatial scales. In
most applications, ROMs with closure methods are tested in canonical problems,
but some exceptions can be seen in [8, 14, 15, 7, 9].

Closure modeling can be divided in physics-based and data-driven. The first
category relies primarily on physical intuition such as the concept of an
energy cascade fundamental to turbulence. Artificial viscosity closure models
[16, 17, 18] fall into this category and the main idea of these methods is to
add the lacking dissipation effects of the truncated POD modes and subgrid
scales effects. Several approaches are proposed in literature and they can be
as simple as changing the viscosity coefficient [16, 17] or more sophisticated
such as Smagorinsky type models [18, 17] which vary in space and time. Another
physics-based closure model is based on the underlying streamline-upwind
Petrov–Galerkin (SUPG) applied to ROMs [19]. In this case, stabilization
parameters based on the underlying finite element discretization and the POD
basis truncation are derived with the appropriate theoretical support given by
finite element error analysis. Unfortunately, the performance of this type of
closure has been shown to be generally unsatisfactory when modeling complex
flows requiring several POD modes [18] and would probably be ineffective if
applied to chaotic turbulent flows.

Data-driven closure models, also known as calibration methods, try to model
all ROM inconsistencies as a whole (top-down approach) using the full order
model (FOM) data (e.g., simulation snapshots or POD temporal modes) as opposed
to attempting to model error sources individually (bottom-up approach). For
this category, the underlying principle is that a ROM should be capable of
recovering most of the information of the FOM from which it originated.
Therefore, it is only natural to exploit the data to its fullest extent before
attempting to extrapolate outside the training window. For example, this can
be achieved by adding corrective (calibration) coefficients to the original
ROM, capable of improving the model by taking the information provided by POD
temporal modes into consideration [20, 21, 22, 23, 15, 24]. Alternative
methods correct the ROM by directly using the snapshots [8, 25, 26]. Hybrid
closure models relying on both physical insight (e.g., energy conservation)
and data are also possible [27, 28]. Furthermore, it is relevant to point out
that minimization problems, commonly appearing when a data-driven approach is
employed, can be ill-conditioned and require some form of regularization. In
[29, 30], some regularization strategies for data-driven calibration of
projection-based ROMs are discussed.

In the present work, we use calibration based on POD temporal modes to improve
the accuracy and stability of Galerkin and Petrov-Galerkin ROMs. In fact, a
more accurate model can be obtained by adding corrective terms to the original
model and minimizing the error between the POD and ROM temporal modes. This
translates to solving a non-linear optimization problem. In [21], an
approximation transforming the original non-linear problem into a less
expensive linear optimization is proposed and it is also adopted here. A
comparison of the linear and non-linear optimization calibration methods is
provided in [31] for a low Reynolds cylinder wake flow and a separated flow
around an ONERA-D airfoil with snapshots obtained from PIV. Briefly, the
linear optimization approach was shown to be much faster while being almost as
accurate as its non-linear counterpart. Additionally, this calibration
methodology was used in reduced order modeling of unsteady flows around
airfoils at low Reynolds numbers [23, 22] and more recently at moderate
Reynolds numbers [32, 15].

Data-driven closure models have been mostly used with linear calibration
coefficients. However, more general (non-linear) coefficients can be
considered and have been tested to some extent [26, 27, 33, 34, 35]. In [34],
this calibration methodology not only is used to obtain linear calibration
terms but also quadratic and cubic operators. Methodology details are given
and results are presented for the one-dimensional Burgers equation.

In this work, calibration of projection-based ROMs is assessed using linear
and non-linear coefficients. For the latter approach, calibration is employed
for third and fourth-order tensors justified by the Galerkin and LSPG
projections, respectively. The model coefficients are pre-computed as result
of adopting the non-conservative form of the compressible Navier-Stokes
equations. A greedy hyper-reduction method is applied in the pre-computational
stage to obtain completely grid-independent ROMs. Methods are tested on two
compressible flow problems being a canonical compressible flow past a cylinder
and the turbulent flow over an airfoil undergoing deep dynamical stall. In
both cases, regularization is required and performed using a Tikhonov
methodology combined with an L-curve approach for appropriate parameter choice
[29, 30]. To the best of the authors’ knowledge, a pre-computed formulation of
the compressible Navier-Stokes equations has never been applied with the LSPG
method. Also, we believe that calibration of the fully-discrete ROMs based on
error minimization of the temporal modes has never been performed with non-
linear coefficients.

##  2 Reduced order modeling

###  2.1 Proper orthogonal decomposition

Proper orthogonal decomposition (POD) is a fundamental step in the calculation
of projection-based ROMs. Here, POD is applied to compute a low-dimensional
subspace of the primitive flow variables $\zeta$, $u$, $v$ and $p$, which
represent specific volume, velocity components and pressure, respectively.
These variables are chosen since the flows of interest in this work are solved
using a non-conservative compressible formulation as detailed in Section 2.4.
In the present notation, the unsteady flows can be decomposed as

$\mathbf{q}(\mathbf{x},t)=\mathbf{\bar{q}(x)}+\mathbf{q^{\prime}}(\mathbf{x},t)=\mathbf{\bar{q}(x)}+\sum_{i=1}^{M}\mathbf{\Phi_{i}(x)}a_{i}(t)\mbox{
,}$ (1)

where
$\mathbf{q}=\left\\{\zeta,u,v,p\right\\}^{\top}\in\mathbb{R}^{N}\times[0,T]$,
$\mathbf{\bar{q}(x)}\in\mathbb{R}^{N}$ is the mean flow,
$\mathbf{\Phi}\in\mathbb{R}^{N\,\times\,M}$ with
$\mathbf{\Phi_{i}}\in\mathbb{R}^{N}$ being the orthonormal spatial
eigenfunctions,
$\mathbf{a}=[a_{1},\ldots,a_{M}]^{\top}\in\mathbb{R}^{M}\times[0,T]$
represents the temporal modes and $\left\\{\cdot\right\\}^{\top}$ is the
transpose of $\left\\{\cdot\right\\}$. The parameter $M$ represents the number
of data sets extracted from the numerical simulation, $N$ is the number of
grid points multiplied by the number of flow variables, $T$ represents the
period of snapshot extraction and $i$ represents the mode index.

The POD consists of looking for the deterministic functions
$\\{\mathbf{\Phi_{i}}\\}_{i=1}^{M}$ that are most similar in an averaged sense
to the realizations $\mathbf{q}(\mathbf{x},t)$. Here, we employ the technique
introduced by [2] called “snapshot POD” where the following Fredholm integral
eigenvalue problem is solved as

$\int_{0}^{T}C_{ij}\,a_{i}(t^{\prime})\,dt^{\prime}=\lambda_{i}\,a_{i}(t)\mbox{
,}$ (2)

and where the temporal covariance matrix
$\mathbf{C}\in\mathbb{R}^{M\,\times\,M}$ is defined by

$C_{ij}=\frac{1}{T}\int_{\Omega}\mathbf{q}(\mathbf{x},t_{i})\,\mathbf{q}(\mathbf{x},t_{j})\,d\mathbf{x}\approx\frac{1}{T}\langle\mathbf{q}(\mathbf{x},t_{i}),\mathbf{q}(\mathbf{x},t_{j})\rangle_{\mathbf{\Pi}}\mbox{
.}$ (3)

In the previous equation, $\mathbf{\Pi}\in\mathbb{R}^{N\,\times\,N}$ is a
symmetric positive definite matrix defining the inner-product
$\langle\mathbf{q}(\mathbf{x},t_{i}),\mathbf{q}(\mathbf{x},t_{j})\rangle_{\mathbf{\Pi}}=\mathbf{\langle
q_{i},q_{j}\rangle_{\Pi}=q_{i}^{T}\Pi q_{j}}$. In the present application, the
matrix $\mathbf{\Pi}$ is diagonal with non-zero elements defined as
$\mathbf{\Pi_{ii}=A_{i}}$, where $\mathbf{A_{i}}$ is the area associated to
the $i$-th vector element. The covariance matrix $\mathbf{C}$ is symmetric
positive semidefinite and, therefore, allows the use of singular value
decomposition to compute the singular values and singular vectors which are,
in turn, related to the eigenvalues $\lambda_{i}$ and eigenvectors (modes) of
the POD reconstruction. Such modes are calculated so that the reconstruction
is optimal in the sense of truncated mean quadratic error. The idea of writing
a temporal covariance matrix (snapshot POD) comes from the fact that the
solution cost grows rapidly for large computational grids. This is an issue
especially in multidimensional problems.

In Eq. 2, $a_{i}$ represent the $i-th$ time-dependent orthogonal POD
eigenfunctions, also called temporal modes. The reader is referred to [36, 6,
2, 3, 13] for more details on the calculation of these modes. The spatial
basis functions $\mathbf{\Phi_{i}}$ can be calculated from the realizations
$\mathbf{q}$ and the temporal modes $\mathbf{a}$ with

$\mathbf{\Phi_{i}}(\mathbf{x})=\frac{1}{T\lambda_{i}}\int_{0}^{T}q(\mathbf{x},t)\,a_{i}(t)\,dt\mbox{
.}$ (4)

In summary, the spatial modes are used in the Galerkin and LSPG projections to
reconstruct a system of ordinary differential equations that, in turn, will
determine the evolution of the temporal modes. With the calculation of both
terms, it is possible to reconstruct the fluctuation fields appearing in Eq. 1
as
$\mathbf{q}\approx\mathbf{\hat{q}}=\mathbf{\bar{q}}+\sum_{i=1}^{m}\mathbf{\Phi_{i}}a_{i}$.
In practical ROM applications, instead of employing the full set $M$ of POD
modes for the fluctuation field reconstruction, one seeks $m\ll M$. Therefore,
only the most energetic POD modes are used in the ROM.

###  2.2 Projection methods

Let us consider the system of non-linear partial differential equations
$\mathbf{F(q)}\in\mathbb{R}^{N}$ defined in a connected open region
$\Omega\subset\mathbb{R}^{N_{g}}$ with a boundary $\Gamma$

$\begin{cases}\mathbf{F(q})=\frac{d\mathbf{q}}{dt}-\mathbf{G(q)}=\mathbf{0}\quad
in\quad\Omega\\\ \mathbf{q}(\mathbf{x},0)=\mathbf{q_{0}}\\\
\mathbf{q=q_{b}}\quad on\quad\Gamma\mbox{ .}\end{cases}$ (5)

In the system above, $\mathbf{q}$ is a function of space and time, $N_{g}$ is
the number of grid points, and the non-linear operator $\mathbf{G(q)}$ is
given by the convective and diffusive operators appearing in the mass,
momentum and energy equations, herein referred to as Navier-Stokes equations.
Let $\mathbf{\Phi}\in\mathbb{R}^{N\,\times\,m}$ with
$\\{\mathbf{\Phi_{i}}\\}_{i=1}^{m}$ defining an orthonormal basis obtained by
POD. The state variable $\mathbf{q}$ is then approximated as the linear
combination of this basis vector as

$\mathbf{q}\approx\mathbf{\hat{q}}=\mathbf{\bar{q}}+\sum_{i=1}^{m}\mathbf{\Phi_{i}}\,a_{i}\mbox{
,}$ (6)

where the explicit dependencies on space and time are omitted for simplicity.

After approximation, the set of equations is written as
$\mathbf{F}(\mathbf{q})\approx\mathbf{R}(\mathbf{\hat{q}})\neq\mathbf{0}$ for
the physical problem being solved. Here,
$\mathbf{R(\hat{q})=R}(\dot{\mathbf{a}},\mathbf{a})$ is the residual after
order reduction and spatial discretization. A solution is sought by enforcing
the residual $\mathbf{R(\hat{q})}$ orthogonality as

$\mathbf{\langle\Psi_{i},R}(\mathbf{\bar{q}}+\sum_{j=1}^{m}\mathbf{\Phi_{j}}a_{j})\rangle_{\mathbf{\Pi}}=\mathbf{0}\mbox{
,}\quad\quad\mbox{i}=1,\ldots,m\mbox{ ,}$ (7)

where $\\{\mathbf{\Psi_{i}}\\}_{i=1}^{m}$ is the test basis. A projection
method is generally called Galerkin (Petrov-Galerkin) when the test and
solution bases are equal (different), i.e., $\mathbf{\Psi=\Phi}$
($\mathbf{\Psi\neq\Phi}$). The POD solution basis should satisfy the boundary
conditions of Eq. 5 and, in Ref. [36], the ability to obtain homogeneous
Dirichlet or Neumann boundary conditions from the snapshots is discussed.

####  2.2.1 Galerkin projection

Galerkin projection is the most popular alternative for reduced order modeling
of time dependent problems. This can be attributed to its implementation
simplicity and solid mathematical foundation. Applying the Galerkin projection
method ($\mathbf{\Psi=\Phi}$) to Eq. 5 we obtain

$\dot{a}_{i}=\langle\mathbf{\Phi_{i}},\mathbf{\hat{G}}(\mathbf{\bar{q}}+\sum_{j=1}^{m}\mathbf{\Phi_{j}}a_{j})\rangle_{\mathbf{\Pi}}\mbox{
,}\quad\quad\mbox{i}=1,\ldots,m\mbox{ ,}$ (8)

where $\mathbf{\hat{G}(\hat{q})}$ is the approximation of operator
$\mathbf{G(q)}$ after order reduction and spatial discretization. Moreover,
the initial conditions are determined by selection and projection of a single
snapshot $\mathbf{q_{0}}$ in the vector basis

$a_{i}(0)=\langle\mathbf{\Phi_{i}(x)},\mathbf{q_{0}}\rangle_{\mathbf{\Pi}}\mbox{
,}\quad\quad\mbox{i}=1,\ldots,m\mbox{ .}$ (9)

The previous system of ordinary differential equations represents the ROM
associated to the FOM and can be solved using a time-marching method. The
right-hand side of Eq. 8 should not scale with the full-order model so as to
achieve reduced computational cost. Following the POD-Galerkin approach, the
ROM obtained for the non-conservative form of the compressible Navier-Stokes
equations (see details in Section 2.4) can be written as

$\dot{a}_{i}+e_{i}+A_{ij}a_{j}+N_{ijk}a_{j}a_{k}=0_{i}\mbox{
,}\quad\quad\mbox{i, j, k}=1,\ldots,m\mbox{ ,}$ (10)

where the time-independent ODE coefficients $\mathbf{e}\in\mathbb{R}^{m}$,
$\mathbf{A}\in\mathbb{R}^{m\,\times\,m}$ and
$\mathbf{N}\in\mathbb{R}^{m\,\times\,m\,\times\,m}$ can be found in [15]. A
two-step (offline/online) computational approach is possible when the
governing equations are linear or have polynomial non-linearities. In the
offline stage, a reduced-order system of equations is obtained after
projecting the reduced-order basis in the governing equations. This stage is
performed only once but can be particularly costly depending on the number of
modes $m$ and the maximum tensor order of the model used in the reconstruction
given that the inner product $\mathbf{\Pi}$ is grid-size dependent. The online
step corresponds to solving the resulting reduced-order system of equations
such as the system of ODEs given by Eq. 10. This step no longer scales with
the FOM and, thus, is cheap in general.

In this work, for reasons that will become clearer later, it is convenient to
promptly discretize $\dot{a}_{i}$ in Eq. 10 by introducing a time-marching
method. As a result, the implicit Euler time-discretization of the previous
equation yields

$f_{i}^{G}(a_{i}^{n})=\frac{a_{i}^{n}-a_{i}^{n-1}}{\Delta
t}+e_{i}+A_{ij}a_{j}^{n}+N_{ijk}a_{j}^{n}a_{k}^{n}=0_{i}\mbox{
,}\quad\quad\mbox{i, j, k}=1,\ldots,m\mbox{ ,}$ (11)

where $\mathbf{f^{G}(a^{n})}$ is the non-linear algebraic system of equations
of the Galerkin ROM.

####  2.2.2 Least-squares Petrov-Galerkin

A Petrov-Galerkin approach uses different test and trial bases
($\mathbf{\Psi\neq\Phi}$). Particularly, the least-squares Petrov-Galerkin
(LSPG) is an alternative to the Galerkin method showing encouraging results
[7, 8, 14] when used to generate ROMs for non-linear time-dependent problems
such as those usually found in fluid flows. A thorough theoretical comparison
of both LSPG and Galerkin can be found in [37]. The LSPG method [7, 38] seeks
to minimize the fully discrete residual
$\mathbf{\hat{R}(a^{n})}\in\mathbb{R}^{N}$ (i.e., the residual after temporal
and spatial discretization) at each $n$-th time-step as

$\operatorname*{\textrm{minimize}}_{\mathbf{a^{n}}}\mathcal{F}\mathbf{(a^{n})}\mbox{
.}$ (12)

The objective function $\mathcal{F}\mathbf{(a^{n})}$ is defined in the
following special form

$\mathcal{F}\mathbf{(a^{n})}=\frac{1}{2}\mathbf{\|\hat{R}(a^{n})\|_{\Pi}^{2}}=\frac{1}{2}\mathbf{\langle\hat{R}(a^{n}),\hat{R}(a^{n})\rangle_{\Pi}}\mbox{
,}$ (13)

or equivalently

$\mathbf{a^{n}}=\operatorname*{\textrm{arg
min}}_{\mathbf{\alpha}^{n}\,\in\,\mathbb{R}^{m}}\|\mathbf{\hat{R}}(\mathbf{\alpha^{n})}\|_{\mathbf{\Pi}}^{2}\mbox{
.}$ (14)

In this work, the fully discrete residual $\mathbf{\hat{R}(a^{n})}$ represents
the $2D$ non-conservative compressible Navier-Stokes equations (Eqs. 22) after
the state variable $\mathbf{q}$ is approximated by the POD basis vector,
finite-difference spatial discretization and implicit time discretization. As
a result, the residual associated with its time discretization by the implicit
Euler scheme is

$\mathbf{\hat{R}(a^{n})}=\mathbf{\Phi}\frac{\mathbf{a^{n}-a^{n-1}}}{\Delta
t}+\mathbf{\hat{G}(a^{n})}\mbox{ ,}$ (15)

where $\Delta t$ is the time step, $\mathbf{\hat{G}}$ is the same operator as
in Eq. 8 and initial conditions are also given by Eq. 9.

Optimality conditions are derived from Taylor’s theorem and can be determined
by examining the gradient $\nabla\mathcal{F}\mathbf{(a^{n}})$ and Hessian
$\Delta\mathcal{F}\mathbf{(a^{n}})$ matrices [39]. The derivatives of
$\mathcal{F}\mathbf{(a^{n}})$ can be expressed in terms of the Jacobian
$\mathbf{J(a^{n})=\frac{\partial\hat{R}(a^{n})}{\partial a^{n}}}$. Applying
the first-order necessary condition
$\mathbf{\nabla}\mathcal{F}(\mathbf{a^{n}})=0$ yields

$\nabla\mathcal{F}(\mathbf{a^{n}})=\mathbf{\langle
J(a^{n}),\hat{R}(a^{n})\rangle_{\Pi}}=\mathbf{\Bigg{\langle}\frac{\partial\hat{R}(a^{n})}{\partial
a^{n}},\hat{R}(a^{n})\Bigg{\rangle}_{\Pi}=0}\mbox{ ,}$ (16)

which is equivalent to Eq. 12, provided that convergence to the same local
minimum occurs [8]. Therefore, in the LSPG method, the discrete test basis
$\mathbf{\Psi(a^{n})}\in\mathbb{R}^{N\,\times\,m}$ is given by

$\mathbf{\Psi(a^{n})}=\mathbf{J(a^{n})}=\mathbf{\frac{\partial\mathbf{\hat{R}(a^{n}})}{\partial\mathbf{a^{n}}}}\mbox{
.}$ (17)

The offline/online approach discussed previously is also possible when using
the LSPG method [40, 41]. However, pre-calculating the LSPG coefficients can
quickly become a cumbersome task, specially when compared to the Galerkin
method, given the complexity introduced by the test basis
$\mathbf{\Psi(a^{n})}$. In fact, the LSPG ROM not only has additional terms
but also comes with the downside of higher-orders tensors that can lead to
prohibitive computational costs when several POD modes are required in the
model construction. This cost is associated with the multiple spatial terms
arising from the inner product from Eq. 16 that contains the Jacobian of the
residual, which is extensive for the discrete compressible Navier-Stokes
equations. Similarly, decomposing the state variable $\mathbf{q}$ following
Eq. 6 contributes to further intensify this issue. For the implicit Euler time
marching scheme, the LSPG ROM of the non-conservative compressible Navier-
Stokes equations can be written as

$\begin{split}f_{i}^{PG}(a_{i}^{n})=\frac{a_{i}^{n}-a_{i}^{n-1}}{\Delta
t}+(e_{i}^{\prime}+A_{ij}^{\prime}a_{j}^{n}+B_{ij}^{\prime}a_{j}^{n-1}+N_{ijk}^{\prime}a_{j}^{n}a_{k}^{n}+L_{ijk}^{\prime}a_{j}^{n}a_{k}^{n-1})\,+&\\\
\Delta
t(e_{i}^{\prime\prime}+A_{ij}^{\prime\prime}a_{j}^{n}+N_{ijk}^{\prime\prime}a_{j}^{n}a_{k}^{n}+Q_{ijkl}^{\prime\prime}a_{j}^{n}a_{k}^{n}a_{l}^{n})=0\mbox{
,}\quad\quad\mbox{i, j, k, l}=1,\ldots,m\mbox{ ,}&\end{split}$ (18)

where $\mathbf{f^{PG}}$ describes the non-linear algebraic system of equations
of the LSPG ROM and $\mathbf{e}^{\prime}\in\mathbb{R}^{m}$,
$\mathbf{A}^{\prime}\in\mathbb{R}^{m\,\times\,m}$,
$\mathbf{B}^{\prime}\in\mathbb{R}^{m\,\times\,m}$,
$\mathbf{N}^{\prime}\in\mathbb{R}^{m\,\times\,m\,\times\,m}$,
$\mathbf{L}^{\prime}\in\mathbb{R}^{m\,\times\,m\,\times\,m}$,
$\mathbf{e}^{\prime\prime}\in\mathbb{R}^{m}$,
$\mathbf{A}^{\prime\prime}\in\mathbb{R}^{m\,\times\,m}$,
$\mathbf{N}^{\prime\prime}\in\mathbb{R}^{m\,\times\,m\,\times\,m}$ and
$\mathbf{Q}^{\prime\prime}\in\mathbb{R}^{m\,\times\,m\,\times\,m\,\times\,m}$
are the pre-computable coefficients defining the algebraic system of
equations. In particular, the non-conservative form of the compressible
Navier-Stokes equations has quadratic non-linearities that lead to the
emergence of the fourth-order tensor $Q_{ijkl}^{\prime\prime}$ in Eq. 18.
Consequently, not only is the LSPG ROM more computationally expensive during
the offline stage but also it has a significant overhead after pre-computation
(online stage) compared to the Galerkin ROM. While the governing equations
produce third-order tensors in the Galerkin model (Eq. 10), they lead to
fourth-order tensors in the LSPG.

###  2.3 Hyper-reduction

Although an order reduction method is applied, models generated using
projection-based techniques can fail to produce relevant computational time
gains. Pre-computation of Galerkin/Petrov-Galerkin coefficients
(offline/online decoupling procedure) is impossible for problems containing
strong non-linearities (i.e. non-polynomial) or non-affine parameter
dependence. On the other hand, even when pre-computation is possible, the
offline stage can be very costly depending on the number of modes used in the
reconstruction and complexity of the coefficients. In particular, Petrov-
Galerkin methods can be very computationally demanding during the offline
stage. Non-linear methods [4, 5] recently being applied in a reduced order
modeling context also require further approximation and techniques are being
developed to specifically address this issue. In all cases, ROMs could benefit
from an additional layer of approximation provided by hyper-reduction methods.
This is typically performed by sampling the spatial modes using a greedy
method before projection.

Many studies [42, 7, 15] adopt a gappy POD framework [43] when using hyper-
reduction. In this approach, given a subset of indices
$J=\\{j_{1},\dots,j_{s}\\}\subset\\{1,\dots,N\\}$, a solution vector
$\mathbf{\hat{q}}$ is approximated by sampling the spatial modes using a mask
projection matrix
$\mathbf{P}=[\mathbf{e_{j_{1}},\dots,e_{j_{s}}}]^{T}\in\mathbb{R}^{N\times s}$
to construct the estimated solution
$\tilde{\mathbf{q}}\approx\mathbf{P^{T}\Phi a}$. Here, $s$ is the number of
indices retained from the original vector of size $N$ and $\mathbf{e_{j_{k}}}$
denotes the vector with a $1$ in the $j_{k}$-th coordinate and $0$’s
elsewhere. Ideally, the number of sampled entries should be kept to a minimum
but produce an optimal approximation of the temporal modes of the modified
mask problem. In other words, the error $\epsilon$ between the gapless
solution $\mathbf{\hat{q}}$ and $\tilde{\mathbf{q}}$, defined as

$\epsilon=\|\mathbf{P^{T}\hat{q}}-\mathbf{P^{T}\Phi
a}\|_{\mathbf{\Pi}}^{2}\,\,\mbox{ ,}$ (19)

should be minimal for a given mask matrix $\mathbf{P}$ sampling $s$ vector
entries. This is equivalent to the minimization problem

$\operatorname*{\textrm{minimize}}_{\mathbf{P}}\|\mathbf{I}-\langle\mathbf{\tilde{\Phi},\tilde{\Phi}}\rangle_{\mathbf{\Pi}}\|_{2}^{2}\,\,\mbox{
,}$ (20)

where $\mathbf{\tilde{\Phi}}=\mathbf{P^{T}\Phi}$. In this work, hyper-
reduction is performed using the POD basis $\mathbf{\Phi}$ for the
reconstruction of the solution vector $\mathbf{\hat{q}}$. However, a basis
built from the FOM residual or fully discrete residual vector
$\mathbf{\bar{R}}$ could also be adopted.

The optimal solution of Eq. 20 for a set of given size can be an intractable
combinatorial optimization problem even for relatively small problems.
Consequently, this problem is usually solved for a sub-optimal set of indices
using a greedy algorithm such as the missing point estimation (MPE) method
[44]. In an iteration of this method, the condition number

$c\left(\mathbf{\langle\tilde{\Phi},\tilde{\Phi}\rangle_{\mathbf{\Pi}}}\right)\equiv\frac{\lambda_{max}\left(\mathbf{\langle\tilde{\Phi},\tilde{\Phi}\rangle_{\mathbf{\Pi}}}\right)}{\lambda_{min}\left(\mathbf{\langle\tilde{\Phi},\tilde{\Phi}\rangle_{\mathbf{\Pi}}}\right)}$
(21)

of the approximated identity matrix
$\langle\mathbf{\tilde{\Phi},\tilde{\Phi}}\rangle_{\mathbf{\Pi}}\approx\mathbf{I}$
(i.e., inner product of the mask POD basis) is evaluated for all the remaining
candidate points $J_{c}=\\{1,\ldots,N\\}-J_{e}$ and eliminates the index with
the highest condition number. Here, $J_{e}$ corresponds to the set of
eliminated indices and is initially an empty set. This procedure is repeated
until a user specified near-optimal set of points is left. Computational cost
of the MPE method may quickly become prohibitive, although the method is
cheaper than the optimal solution. Furthermore, the condition number can serve
as an a priori non-sharp error indicator of the hyper-reduced problem.

In the present work, we adopt the accelerated greedy MPE technique for hyper-
reduction. This procedure [45] allows for substantial time gains compared to
the naive MPE while still using the same underlying principles. Similarly to
the condition number problem previously discussed, the greedy point selection
is equivalent to picking the spatial mode entry responsible for the largest
growth in the condition number of the modified eigenvalue problem. In the
accelerated version of the MPE, the costly modified symmetric eigenvalue
problem is not actually solved, but the index selection occurs by analyzing
properties of the candidate indices.

###  2.4 Non-conservative compressible Navier-Stokes equations

Consider the $2D$ non-conservative compressible Navier-Stokes equations
presented by [46] as

$\displaystyle\begin{split}\zeta_{t}=&\zeta(u_{x}+v_{y})-u\zeta_{x}-v\zeta_{y}\mbox{
,}\end{split}$ (22a) $\displaystyle\begin{split}u_{t}=&-uu_{x}-vu_{y}-\zeta
p_{x}+\frac{Ma}{Re}\zeta\bigg{[}\bigg{(}\frac{4}{3}u_{x}-\frac{2}{3}v_{y}\bigg{)}_{x}+(v_{x}+u_{y})_{y}\bigg{]}\mbox{
,}\end{split}$ (22b) $\displaystyle\begin{split}v_{t}=&-uv_{x}-vv_{y}-\zeta
p_{y}+\frac{Ma}{Re}\zeta\bigg{[}\bigg{(}\frac{4}{3}v_{y}-\frac{2}{3}u_{x}\bigg{)}_{y}+(v_{x}+u_{y})_{x}\bigg{]}\mbox{
, and}\end{split}$ (22c)
$\displaystyle\begin{split}p_{t}=&-up_{x}-vp_{y}-\gamma
p(u_{x}+v_{y})+\frac{\gamma Ma}{RePr}[(p\zeta)_{xx}+(p\zeta)_{yy}]\\\
&+\frac{(\gamma-1)Ma}{Re}\bigg{[}u_{x}\bigg{(}\frac{4}{3}u_{x}-\frac{2}{3}v_{y}\bigg{)}+v_{y}\bigg{(}\frac{4}{3}v_{y}-\frac{2}{3}u_{x}\bigg{)}+(v_{x}+u_{y})^{2}\bigg{]}\mbox{
,}\end{split}$ (22d)

where $\zeta=1/\rho$ is the specific volume, $u$ and $v$ are the $x$ and $y$
velocity components, respectively, and $p$ is the pressure. In the set of Eqs.
22, $Pr$, $Re$ and $Ma$ denote the reference Prandtl, Reynolds and Mach
numbers, respectively. Subscripts denote partial derivatives and $\gamma$ is
the specific heat ratio.

The quadratic non-linearities of the non-conservative compressible Navier-
Stokes equations allow the pre-computation of Galerkin and LSPG coefficients.
This set of equations was previously used for construction of ROMs in Refs.
[46, 24, 15, 32]. However, to the authors knowledge, this work is the first to
perform pre-computation of LSPG coefficients for Eq. 22. In CFD applications,
the compressible Navier-Stokes equations are usually expressed in conservative
form but this would add tensors of even higher orders than those computed in
this work. This issue has been typically remedied by using an additional
approximation layer such as hyper-reduction [47, 48]. Further approximation
can be detrimental to the models and it has motivated corrective techniques
such as the empirical quadrature procedure (EQP) [49], stabilization by
entropy preservation [50] and the energy-conserving sampling and weighting
(ECSW) methods [51, 8].

###  2.5 Linear least-squares calibration of projection-based ROMs

In general, calibration methods attempt to improve stability and accuracy of
time-dependent reduced order models by adding corrective terms to the non-
linear system of algebraic (or differential) equations obtained by projection-
based or data-driven models. Considering the tensors appearing in the LSPG
formulation, the corrective terms would be given by

$f_{i}^{n}+e_{i}^{c}+A_{ij}^{c}a_{j}^{n}+N_{ijk}^{c}a_{j}^{n}a_{k}^{n}+Q_{ijkl}^{c}a_{j}^{n}a_{k}^{n}a_{l}^{n}=0\mbox{
,}\quad\quad\mbox{i, j, k, l}=1,\ldots,m\mbox{ ,}$ (23)

where $\mathbf{f}$ could be either the Galerkin ROM $\mathbf{f^{G}}$ or the
LSPG ROM $\mathbf{f^{PG}}$. The terms $\mathbf{e^{c}}\in\mathbb{R}^{m}$,
$\mathbf{A^{c}}\in\mathbb{R}^{m\,\times\,m}$,
$\mathbf{N^{c}}\in\mathbb{R}^{m\,\times\,m\,\times\,m}$ and
$\mathbf{Q^{c}}\in\mathbb{R}^{m\,\times\,m\,\times\,m\,\times\,m}$ are the
calibration operators added to the initial ROM. In this section, calibration
is presented for a specific set of calibration terms (first to fourth-order
tensors), but can be applied without loss of generality for a different set of
calibration tensors. For example, the same reasoning is valid when only linear
terms (i.e., $\mathbf{e^{c}}$ and $\mathbf{A^{c}}$) are considered. In fact,
the choice for linear calibration terms is prevalent in the literature.
However, calibration by non-linear operators have been tested to a smaller
degree with higher-order tensors [27, 26, 33, 34] and deep residual neural
networks [35].

In a linear least-squares calibration, the goal may be to minimize the error
$E_{1}^{ROM}$ between temporal modes obtained by solving the system of non-
linear equations of the ROM $\mathbf{a^{ROM}}(t)$ and the original POD
temporal modes $\mathbf{a}(t)$ in a user specified training window $0\leq
t\leq T$ as

$E_{1}^{ROM}=\sum_{i=1}^{m}\int_{0}^{T}\left(a_{i}(t)-a_{i}^{ROM}(t)\right)^{2}dt\mbox{
.}$ (24)

In the above equation, $\mathbf{a}(t)$ refers to the temporal modes obtained
directly by the POD snapshot method in Eq. 2, while the term
$\mathbf{a^{ROM}}(t)$ represents the temporal modes obtained by the solutions
of Eqs. 11 and 18 for the Galerkin and LSPG models, respectively. A second
error norm $E_{2}^{ROM}$ can be also defined by the ODE system and is
frequently preferred when performing calibration [15, 27, 31, 33, 52, 26]

$E_{2}^{ROM}=\sum_{i=1}^{m}\int_{0}^{T}\left(\dot{a}_{i}(t)-\dot{a}^{ROM}_{i}(t)\right)^{2}dt\mbox{
.}$ (25)

However, given that the LSPG model is naturally an algebraic system of
equations, error norm $E_{1}^{ROM}$ is used in this work and it was also
applied in [23, 22] in the context of Galerkin-ROM calibration. The non-linear
optimization problem that arises when minimizing the error norm $E_{1}^{ROM}$
from Eq. 24 can be transformed to a linear one by the approximation
$\mathbf{a^{ROM}}(t)\approx\tilde{\mathbf{a}}(t)$, as discussed in [21]. The
new error norm $E_{1}$ is expressed as

$E_{1}=\sum_{i=1}^{m}\int_{0}^{T}\left(a_{i}(t)-\tilde{a}_{i}(t)\right)^{2}dt\mbox{
,}$ (26)

where $\tilde{a}_{i}(t_{n})$ is given by

$\frac{\tilde{a}_{i}^{n}-a_{i}^{n-1}}{\Delta
t_{s}}-(e_{i}+A_{ij}\tilde{a}_{j}^{n}+N_{ijk}\tilde{a}_{j}^{n}\tilde{a}_{k}^{n})=0\mbox{
,}\quad\quad\mbox{i, j, k}=1,\ldots,m\mbox{, \,\,\,\, or}$ (27a)
$\begin{split}&\frac{\tilde{a}_{i}^{n}-a_{i}^{n-1}}{\Delta
t_{s}}+(e_{i}^{\prime}+A_{ij}^{\prime}\tilde{a}_{j}^{n}+B_{ij}^{\prime}a_{j}^{n-1}+N_{ijk}^{\prime}\tilde{a}_{j}^{n}\tilde{a}_{k}^{n}+L_{ijk}^{\prime}\tilde{a}_{j}^{n}a_{k}^{n-1})\,+\\\
&\Delta
t_{s}(e_{i}^{\prime\prime}+A_{ij}^{\prime\prime}\tilde{a}_{j}^{n}+N_{ijk}^{\prime\prime}\tilde{a}_{j}^{n}\tilde{a}_{k}^{n}+Q_{ijkl}^{\prime\prime}\tilde{a}_{j}^{n}\tilde{a}_{k}^{n}\tilde{a}_{l}^{n})=0\mbox{
,}\quad\quad\mbox{i, j, k, l}=1,\ldots,m\mbox{ .}\end{split}$ (27b)

In the equations above, $\Delta t_{s}$ is the time step between snapshots and
$\mathbf{a^{n-1}}$ is the POD temporal mode at $t_{n-1}$. A pre-computable
approximation of the ROM temporal modes can be obtained for all time instants
using Eq. 27a (Galerkin) or 27b (LSPG). Also, in both cases, time
discretization is performed using the implicit Euler method and for the first
time step $\tilde{a}_{i}^{1}=a_{i}^{1}$. Suppression of the non-linear
constraint to obtain an affine operator may seem extreme at first, but this
can be understood as a measure of how the non-calibrated ROM alters the
solution in each time step. If successful, calibration should systematically
compensate for any major deviation from the POD temporal modes.

Here, $E_{1}$ will be used to indicate the non-calibrated affine approximation
error while $E_{1}^{c}$ will indicate the same error with calibration terms
$\mathbf{e^{c}}\in\mathbb{R}^{m}$,
$\mathbf{A^{c}}\in\mathbb{R}^{m\,\times\,m}$,
$\mathbf{N^{c}}\in\mathbb{R}^{m\,\times\,m\,\times\,m}$ and
$\mathbf{Q^{c}}\in\mathbb{R}^{m\,\times\,m\,\times\,m\,\times\,m}$ taken into
consideration

$\begin{split}&E_{1}^{c}(\mathbf{e^{c},A^{c},N^{c},Q^{c}})=\sum_{i=1}^{m}\int_{0}^{T}\Big{[}a_{i}(t)-\tilde{a}_{i}(t)\\\
&-\int_{t}^{t+\Delta
t_{s}}\Big{(}e_{i}^{c}+A_{ij}^{c}a_{j}(t^{\prime})+N_{ijk}^{c}a_{j}(t^{\prime})a_{k}(t^{\prime})+Q_{ijkl}^{c}a_{j}(t^{\prime})a_{k}(t^{\prime})a_{l}(t^{\prime})\Big{)}dt^{\prime}\Big{]}^{2}dt\mbox{
.}\end{split}$ (28)

Adding a regularization term when minimizing error $E_{1}^{c}$ can prevent
overfitting of the solution which in turn could lead to a better model. Taking
this into account, a functional
$\mathcal{L}(\mathbf{e^{c},A^{c},N^{c},Q^{c}})$ is defined as

$\mathcal{L}(\mathbf{e^{c},A^{c},N^{c},Q^{c}})=E_{1}^{c}(\mathbf{e^{c},A^{c},N^{c},Q^{c}})+\tilde{\theta}(\mathbf{\|e^{c}\|_{F}^{2}+\|A^{c}\|_{F}^{2}+\|N^{c}\|_{F}^{2}+\|Q^{c}\|_{F}^{2}})\mbox{
,}$ (29)

with

$\tilde{\theta}=\frac{1-\theta}{\theta}\frac{E_{1}}{\mathbf{\|e\|_{F}^{2}+\|A\|_{F}^{2}+\|N\|_{F}^{2}+\|Q\|_{F}^{2}}}\mbox{
,}$ (30)

where the second term on the right hand side of Eq. 29 corresponds to the
regularization term and $\|\cdot\|_{\mathbf{F}}$ represents the Frobenius
norm. The parameter $\theta\in(0,1]$, and values of $\theta$ close to $0$ add
more importance to the original ROM while $\theta$ close to $1$ adds a higher
weight to the prediction quality of the calibrated model along the training
window. It is also relevant to point out that the calibration parameter
$\theta$ can be regarded as a regularization parameter in the Tikhonov
regularization methodology [23].

Minimizing the functional $\mathcal{L}$ gives rise to the linear system

$\mathcal{L}(\mathbf{K^{c}})=\sum_{i=1}^{m}\int_{0}^{T}\left(a_{i}(t)-\tilde{a}_{i}(t)-\int_{t}^{t+\Delta
t_{s}}\sum_{j=1}^{m^{*}}K_{ij}^{c}a_{j}^{*}(t^{\prime})dt^{\prime}\right)^{2}dt+\tilde{\theta}\mathbf{\|K^{c}\|_{F}^{2}}\mbox{
,}$ (31)

where
$\mathbf{a^{*}}=[1\;\;\mathbf{a}\;\;\mathbf{a}\hat{\otimes}\mathbf{a}\;\;\mathbf{a}\hat{\otimes}\mathbf{a}\hat{\otimes}\mathbf{a}]\in\mathbb{R}^{m^{*}}\times[0,T]$
is the enriched temporal mode vector. Here, $\otimes$ symbolize the column-
wise Kronecker product, which for a column vector
$\mathbf{a}=[a_{1},\ldots,a_{m}]^{\top}\in\mathbb{R}^{m}$ is given by

$\mathbf{a}\otimes\mathbf{a}=[a_{1}^{2}\;\;a_{1}a_{2}\;\ldots\;a_{1}a_{m}\;\;a_{2}a_{1}\;\ldots\;a_{2}^{2}\;\ldots\;a_{2}a_{m}\;\ldots\;a_{m}^{2}]\in\mathbb{R}^{m^{2}}\mbox{
,}$ (32)

and $\hat{\otimes}$ corresponds to the (modified) Kronecker product without
duplicate terms (for instance, $a_{1}a_{2}=a_{2}a_{1}$ and
$a_{1}a_{1}a_{2}=a_{1}a_{2}a_{1}=a_{2}a_{1}a_{1}$). Similarly, the different
calibration coefficients are grouped such that
$\mathbf{K^{c}=[e^{c},\;A^{c},\;\hat{N}^{c},\;\hat{Q}^{c}]}\in\mathbb{R}^{m\,\times\,m^{*}}$.
Moreover, $\hat{\mathbf{N}}^{c}\in\mathbb{R}^{m\,\times\,m(m+1)/2}$ and
$\hat{\mathbf{Q}}^{c}\in\mathbb{R}^{m\,\times\,m(m+1)(m+2)/6}$ are,
respectively, the matricized equivalents of tensors $\mathbf{N}^{c}$ and
$\mathbf{Q}^{c}$ without redundant coefficients. Failure to eliminate
duplicate terms leads to a dependent (ill-posed) linear system of equations.
For the calibration terms being considered, the dimension $m^{*}$ is given by

$m^{*}=1+m+\frac{m(m+1)}{2}+\frac{m(m+1)(m+2)}{6}\mbox{ .}$ (33)

For a given parameter $\theta$, the optimality condition
$\partial\mathcal{L}/\partial K_{ij}^{c}=0$ leads to solving $m$ linear
systems of size $m^{*}$ or, for $i=1,\dots,m$

$\mathbf{D^{T}K_{i}^{c}=b^{i}}\mbox{ ,}$ (34)

where $\mathbf{K_{i}^{c}}$ is the $i^{th}$ row of $\mathbf{K^{c}}$,

$D_{ij}=\int_{0}^{T}\left(\int_{t}^{t+\Delta
t_{s}}a_{i}^{*}(t^{\prime})dt^{\prime}\right)\left(\int_{t}^{t+\Delta
t_{s}}a_{j}^{*}(t^{\prime})dt^{\prime}\right)dt+\tilde{\theta}\delta_{ij}$
(35)

and

$b_{i}^{j}=\int_{0}^{T}\bigg{(}a_{j}(t)-\tilde{a}_{j}(t)\bigg{)}\bigg{(}\int_{t}^{t+\Delta
t_{s}}a_{i}^{*}(t^{\prime})dt^{\prime}\bigg{)}dt\mbox{ .}$ (36)

Matrix $\mathbf{D}\in\mathbb{R}^{m^{*}\,\times\,m^{*}}$ is calculated only
once while vector $\mathbf{b^{i}}$ must be evaluated for each mode during the
calibration phase. Further details and a thorough discussion can be found in
[53].

An appropriate value for the calibration (regularization) parameter $\theta$
needs to be chosen. This can be done systematically with the L-curve method
and is discussed by [29, 30] in a reduced order modeling context. This method
is based on the examination of the curve representing the norm of the
regularized solution $\mathbf{\|K^{c}\|_{F}}$ versus a corresponding residual
norm ($E_{1}^{c}$ in this paper). Typically, this curve exhibits an “L” shape
such as the blue curves in Fig. 1 if the problem is ill-conditioned or
contains noise. The corner of the L-curve usually represents a fair trade-off
between the solution norm $\mathbf{\|K^{c}\|_{F}}$ (ordinate) and minimization
of the error norm $E_{1}^{c}$ (abscissa). For a typical well-behaved L-curve,
the identification of the corner can be performed by curvature maximization
and corresponds to the maximum curvature point [30].

Usually the L-curve is obtained after minimizing the functional
$\mathcal{L(\mathbf{K^{c}})}$ in Eq. 31 for different values of
$\theta\in(0,1]$ [15, 32, 29]. In this work, an iterative heuristic approach
is preferred. To begin with, a first set of calibration terms is determined
for a suitably chosen calibration parameter $\theta_{1}$ and the non-
calibrated ROM. Afterwards, these calibration coefficients are added to the
original ROM coefficients. Next, a second set of calibration terms is obtained
using the modified reduced order coefficients as a replacement for the non-
calibrated coefficients combined with a second regularization parameter
$\theta_{2}$. At each iteration, the error norm $E_{1}^{c}$ and the norm of
the calibration coefficients (recovered after subtracting the original terms)
are computed. The iterations only come to a halt when an appropriate L-curve
is obtained. In practice, the regularization parameters are chosen
heuristically. If too big, the calibration coefficients could already be to
the left side of the L-curve corner. On the other hand, too small
regularization parameters may lead to a substantial number of iterations.
Moreover, this approach allows to systematically take into consideration the
errors associated with the time integration of the calibration coefficients.
As a consequence, calibration coefficients obtained by this procedure could
perform better compared to their counterparts obtained by solving a single
optimization problem.

##  3 Results and Discussion

This section presents results of Galerkin and LSPG ROMs for two test cases: a
compressible cylinder flow at low Reynolds number, and a turbulent flow over a
plunging airfoil under deep dynamic stall. An assessment of linear and non-
linear calibration is presented for the ROMs computed with and without hyper-
reduction.

###  3.1 Cylinder flow

A compressible cylinder flow is studied to perform an initial assessment of
the Galerkin and Petrov-Galerkin ROMs. The effects of linear and non-linear
calibration are analyzed for both models. The flow has a freestream Mach
number $M_{\infty}=0.4$ and the Reynolds number based on the cylinder diameter
is set as $Re_{D}=150$. Results obtained from the FOM are sampled for
$1{,}120$ snapshots with constant non-dimensional time step $\Delta
t_{s}=0{,}125$. The first $560$ snapshots, which represent $70$ non-
dimensional time units based on freestream velocity, are used to construct a
reduced order basis by the snapshot POD method introduced in section 2.1. The
computational O-grid employed in the simulation has $420\times 599$ points in
the streamwise and wall-normal directions, respectively.

The optimality property of the POD method is expected to produce a basis where
only a small number of modes is necessary to reconstruct the input data and,
thus, the benefit of adding further modes is expected to rapidly decay.
Usually, the number of POD modes used in the ROM is chosen according to the
relative information content
$\mbox{RIC}(m)=\sum_{i=1}^{m}\lambda_{i}/\sum_{i=1}^{M}\lambda_{i}$ and should
satisfy a predefined threshold. For the problem at hand, only four modes
($98{,}87\%$ of the model RIC) are used in the ROM reconstructions. The
implicit Euler method is the chosen time-marching scheme and spatial
derivatives of the POD basis are computed using a $10th$-order explicit
central finite difference scheme for both Galerkin and LSPG methods. Although
Zucatti et al. [11] show that high-resolution methods can be beneficial,
calibration should be able to compensate for this source of error. Some tests
were also made using a $2nd$-order central finite difference scheme and
resulted in ROMs producing the same outcome but requiring slightly more
invasive calibration operators (i.e., with higher $\|\cdot\|_{\mathbf{F}}$
norms) compared to the $10th$-order method. Also, higher-order spatial
derivatives imply in little pre-computation overhead to the models. Finally,
reduced order models obtained either by Galerkin or LSPG schemes are
expressed, after spatial and temporal discretization, as a set of coupled non-
linear algebraic system of equations which are solved using the Levenberg-
Marquardt algorithm.

The performance of both linear and non-linear calibration coefficients is
evaluated. Here, linear calibration (using only constant and linear
coefficients) translates to minimizing the functional $\mathcal{L}$ of Eq. 31
for $\mathbf{e^{c}}$ and $\mathbf{A^{c}}$ ($\mathbf{K^{c}=[e^{c},A^{c}]}$).
Similarly, non-linear calibration (considering linear and non-linear
coefficients) alludes to the same minimization problem but with higher order
tensors ($\mathbf{N^{c}}$ and $\mathbf{Q^{c}}$) also being taken into
consideration. Specifically, non-linear calibration of the LSPG ROMs also adds
operators $\mathbf{N^{c}}$ and $\mathbf{Q^{c}}$
($\mathbf{K^{c}=[e^{c},A^{c},\hat{N}^{c},\hat{Q}^{c}]}$) but non-linear
calibration of the Galerkin ROMs only adds $\mathbf{N^{c}}$
($\mathbf{K^{c}=[e^{c},A^{c},\hat{N}^{c}]}$).

As discussed in section 2.5, calibration may require regularization to avoid
adding coefficients resulting from an ill-posed minimization problem. A
suitable set of calibration coefficients for each ROM is determined after
analyzing the L-curves from Fig. 1. Here, $\mathbf{K_{G}}$ and
$\mathbf{K_{PG}}$ correspond to the grouped calibration terms of the original
non-calibrated Galerkin and LSPG ROMs, respectively. For a similar
approximation error $E_{1}^{c}$, the sum of the Frobenius norms
$\|\cdot\|_{F}$ of the linear calibration coefficients is larger than that
computed for the non-linear calibration. Therefore, the non-linear calibration
operators are less intrusive until the L-curve reaches a corner. As can be
also observed in Fig. 1, the coefficients from calibration are nearly an order
of magnitude smaller compared to those from the original ROMs for both the
Galerkin and LSPG models. Moreover, the latter method exhibits calibration
coefficients which are less intrusive compared to their counterparts from the
Galerkin method. It should be pointed out that $E_{1}^{c}$ serves only as an
error indicator and a better model performance with respect to other
calibrated ROMs is uncertain.

When calibration with non-linear coefficients is performed, a rapid growth of
the relative norm is observed without further reduction of $E_{1}^{c}$ when a
threshold is met. In these cases, the coefficients corresponding to the corner
of the L-curve are selected as the suitable calibration operators. On the
other hand, linear calibration operators converge in norm for both LSPG and
Galerkin methods and coefficients are selected after convergence.

(a) Galerkin L-curves.

(b) LSPG L-curves.

Figure 1: Ratio of Frobenius norms computed for calibrated and original ROM
coefficients as a function of the approximation error using the iterative
procedure described in section 2.5.

Figure 2 shows the time histories of u-velocity fluctuations for the FOM,
Galerkin and LSPG models with and without calibration. Solutions are presented
for the training region ($0\leq t\leq 70$) and for an extrapolation period for
which FOM solutions are available ($70<t\leq 140$). In this case, calibration
of both linear and non-linear terms is evaluated. All calibrated models are
visually indistinguishable from the FOM solution in Fig. 2. In comparison,
very dissipative ROM solutions are obtained when the implicit Euler method is
employed without calibration. Specifically, the uncalibrated LSPG model is
less dissipative than the uncalibrated Galerkin. The present first-order
implicit time-marching scheme is generally inappropriate for long-time
integration in spite of the good stability properties. It is important to
remind that each calibrated ROM is computed for a specific time step $\Delta
t_{s}$ and altering this parameter after the model is already generated can be
harmful. For example, time-step refinement diminishes the dissipative nature
of the implicit Euler method which may lead to unstable solutions. As a
positive effect, calibration allows for improved solution accuracy even with
simple time-marching schemes.

(a) Galerkin models.

(b) LSPG models.

Figure 2: Time histories of u-velocity fluctuations computed by the implicit
Euler time-marching scheme at a probe located on the cylinder wake.

Figure 3 shows the u-velocity fluctuation contours for the FOM and ROMs
considered. Solutions are computed at $t=140$ and the contour levels are kept
equal for all cases to provide a fair comparison. All calibrated ROM solutions
present a good comparison to the FOM with only very small discrepancies. On
the other hand, the uncalibrated ROMs perform very poorly. Despite the small
phase error, one can clearly see the effects of the diminishing magnitudes of
the fluctuation field accused by the probes in Fig. 2. In particular, the
magnitudes of the fluctuation field of the uncalibrated Galerkin ROM is
generally smaller compared to that of the uncalibrated LSPG ROM.

(a) Galerkin ROM with calibrated linear terms.

(b) LSPG ROM with calibrated linear terms.

(c) Galerkin ROM with calibrated non-linear terms.

(d) LSPG ROM with calibrated non-linear terms.

(e) Uncalibrated Galerkin ROM.

(f) Uncalibrated LSPG ROM.

(g) Full-order model.

Figure 3: Velocity fluctuations $u^{\prime}$ obtained for different ROMs and
the FOM at $t=140$. The color limits are kept the same for all figures
$-0.20<u^{\prime}<0.20$.

###  3.2 Deep dynamic stall of a plunging airfoil

This second test case investigates a periodically plunging SD7003 airfoil
undergoing deep dynamic stall. The airfoil is subject to a sinusoidal motion
with reduced frequency $k=\pi fL/U_{\infty}=0.5$, where $f$ in the plunging
frequency, $L$ is the chord of the airfoil and $U_{\infty}$ in the reference
free stream velocity. The motion amplitude is set as $h_{0}/L=0.5$ with a
static angle of attack $\alpha_{0}=8^{\circ}$. The chord Reynolds number based
on the freestream velocity is $Re=\text{60,000}$ and the freestream Mach
number is $M=0.1$. This flow configuration is interesting for micro air
vehicle applications and further details of the flow setup and the LES
methodology employed in the present case can be found in [54].

Figure 4: Iso-surfaces of Q-criterion colored by pressure at two instants of
the plunge motion.

The mesh configuration is a body-fitted O-grid with $481\times 351\times 96$
points in the streamwise, wall-normal and spanwise directions, respectively.
The iso-surfaces of Q-criterion colored by pressure are shown in Fig. 4. In
the figure, some of the main flow features can be observed, for example the
massive leading edge vortex which is transported over the airfoil suction
side, and a trailing edge vortex that develops during the dynamic stall
process. A detailed analysis of the present flow can be found in [54].

Figure 5: Spectrum of singular values (left) and relative information content
(right).

This problem was studied previously in the context of both data-driven [9] and
projection-based reduced order modeling [32]. In [9], a ROM is constructed via
deep neural networks (DNNs) and a comparison with sparse regression [10] is
provided. The DNN method was capable of learning and reproducing the transient
features of the flow with great accuracy in both training and testing windows.
On the other hand, a calibration approach with linear coefficients based in
minimizing the error norm given by Eq. 25, and without pre-computation of the
LSPG coefficients, was performed in [32]. Results showed that linear
calibration significantly improves stability and accuracy for both Galerkin
and LSPG ROMs. The impact of hyper-reduction was also evaluated. In
particular, the work developed in [32] was fundamental to the current article
and will be further discussed throughout the text when appropriate.

Reduced order modeling is tested using 2,244 spanwise-averaged snapshot
solutions (i.e., 2D solutions) of the flow with a constant non-dimensional
time step $\Delta t=0.036$. Half of these snapshots (three plunging cycles)
are used to construct a reduced order basis by the snapshot POD method. An
accurate reconstruction of the dataset from a small POD basis is expected and,
as shown by the present results, the benefits of including additional modes
rapidly vanish. The fast magnitude decay of the singular values for the first
30 most energetic modes (out of 1,122) is presented in Fig. 5. The magnitudes
are similar for mode pairs, what indicates that such modes contain similar
frequency information, differing almost exclusively with respect to phase.
Additionally, the evolution of the RIC for these singular values is also
presented in Fig. 5 and correspond to a RIC of $99.23\%$.

Some POD spatial eigenfunctions are shown in Fig. 6. These modes are presented
for u and v velocity fluctuations, and pressure ($\zeta$ has similar spatial
distribution as the latter variable). The first mode pair contains the large-
scale structure related to the dynamic stall vortex. The other modes also
contain information related to advection of the dynamic stall vortex combined
with the trailing edge vortex and some higher frequency content from the
turbulent flow. The higher POD modes are composed of high-frequency
information and their convergence is not guaranteed due to limitations in the
data acquisition frequency. The corresponding POD temporal modes are shown in
Fig. 7. The first mode has a quasi-sinusoidal behavior related to the plunging
motion. Higher POD temporal modes display a more complex dynamics including
multiple frequencies in a wavepacket-like content and some inherent noise from
the turbulent simulation.

Figure 6: Contours of POD spatial modes for modes 1, 4, 8 and 12 (top to
bottom) for u-velocity (left), v-velocity (center) and pressure fluctuations
(right).

Figure 7: Temporal dynamics of POD modes 1, 4, 8 and 12.

For this flow problem, the ROMs are constructed using 4, 8 and 12 POD modes.
Similarly to the prevous test case, the spatial derivatives of the POD modes
are computed using a $10th$-order central finite difference scheme and the
equations are integrated in time with the first-order implicit Euler method
for both Galerkin and LSPG techniques.

Coefficients arising from the projection are pre-computed for two conditions:
without further approximation (gapless) and with application of hyper-
reduction. In the previous case, modes are sampled by the accelerated greedy
MPE for 1,000 points ($0.594\%$ of total points) to render pre-computation
grid size independent and, hence, to evaluate the impact of hyper-reduction on
the calibration techniques. Furthermore, attempts to use additional modes led
to ROMs being either unstable or inaccurate despite calibration. This can be
attributed, at least in part, to the spatio-temporal resolution of the data.
For this turbulent flow, the temporal and spatial resolutions are lower than
those considered for the previous cylinder flow or in the test cases discussed
in [11, 15]. In [32], ROMs with 16 modes and similar linear calibration terms
showed good results in both training and testing windows, but were unstable
for long-term prediction.

Figures 8 and 9 show the impact of the calibration coefficients relative to
the ROMs’ original coefficients without and with hyper-reduction,
respectively. Similarly to the cylinder flow problem, the curves only exhibit
a typical L-shape when non-linear calibration coefficients are considered.
Differently from the cylinder flow, the present calibrated ROMs require
significantly more invasive coefficients (at least the same order of
magnitude) to provide substantial reduction of the error indicator
$E_{1}^{c}$. This could be attributed to the lower quality of the original
reduced order model coefficients, attested in Fig. 10. This observation is
even more evident when hyper-reduction is taken into consideration. In fact,
hyper-reduced models produce very invasive calibration coefficients despite
being similar in the shape of the L-curves compared to their gapless
counterparts. Another observation from Figs. 8 and 9 is that calibration of
linear and non-linear terms present similar trajectories at first. Eventually,
the linear coefficients converge in norm while the non-linear ones are capable
of achieving further error reductions at the cost of more intrusive
coefficients. As can be seen from the figures, without calibration, the models
with 4 modes are intriguingly more precise compared to models with 8 and 12
modes. This shows that these additional modes are, at least initially,
negatively impacting the models. On the other hand, the extra modes are
capable of further reducing the error indicator with less intrusive
calibration coefficients.

(a) Galerkin curves.

(b) LSPG curves.

Figure 8: Ratio of Frobenius norms computed for calibrated and non-calibrated
coefficients without hyper-reduction as a function of the approximation error
using the iterative procedure described in section 2.5.

(a) Galerkin curves.

(b) LSPG curves.

Figure 9: Ratio of Frobenius norms computed for calibrated and non-calibrated
coefficients hyper-reduced to 1,000 points as a function of the approximation
error using the iterative procedure described in section 2.5.

Time histories of u-velocity fluctuations computed by the FOM and ROMs are
shown in Figs. 10 – 15 for probes placed near the airfoil leading and trailing
edges. Additionally, Table 1 contains the relative Frobenius norms $\rho$ and
approximation errors $E_{1}^{c}$ for the ROMs. As can be see from the figures,
uncalibrated ROMs are stable but completely inaccurate for both the Galerkin
and LSPG methods, independently of the number of POD modes and hyper-
reduction. In Figs. 10 and 11 ROMs are constructed using 4 POD modes (RIC =
$74.65\%$). For this number of modes, both linear and non-linear calibration
models adopt fully modeled errors (i.e., no regularization required for stable
and accurate ROMs). In particular, non-linear coefficients are capable of
substantially reducing the error indicator but the resulting calibration
coefficients are very invasive (specially when hyper-reduction is considered).
This manifests itself visually by the gradual increase of the relative
Frobenius norms shown in Table 1. For hyper-reduced models, purely linear
coefficients have slightly higher errors when compared to gapless models. As
shown in the u-velocity fluctuations of Figs. 10 and 11, solutions obtained
with linear coefficients recover some of the lower frequency oscillations but
suffer to learn finer flow features. On the other hand, the non-linear
calibration solutions are capable of not only better recover the lower
frequency fluctuations but also a significant part of the finer scales. For
the gapless models, Galerkin and LSPG results are almost visually identical
when non-linear calibration is employed. On the other hand, both linear and
non-linear calibrated hyper-reduced models are visually identical.

Models constructed with 8 modes can be seen in Figs. 12 and 13 (RIC =
$92.42\%$). Results for linear calibration coefficients are shown for models
with norm convergence and still suffer to reproduce the finer flow scales.
Compared to the models with four POD modes, the gapless models show a small
improvement of the error indicator after calibration while the error worsens
when hyper-reduction is considered (see table 1). For the non-linear
coefficients, the L-curves are steep and it is possible to obtain less
intrusive coefficients with smaller errors. On the other hand, the additional
modes used in the reconstruction makes regularization a requirement to obtain
stable and accurate models. Stable models with good accuracy could be obtained
using the guidance of the L-curve and running a few models to check for
stability. Probes show results for stable models with less intrusive
calibration coefficients but marginally larger error indicators (see Table 1).

The last set of ROMs analysed is obtained using a 12-mode basis (RIC =
$97.34\%$) and probed results can be seen in Figs. 14 and 15. Again, linear
calibration coefficients are fully modeled while non-linear calibration
coefficients require regularization. Discrepancies between gapless and hyper-
reduced models are considerably larger in this case. Gapless models obtained
with linear calibration coefficients benefit the most from having a more
complete basis of POD modes. In fact, the error indicator $E_{1}^{c}$ is
almost halved compared to the previous cases as shown in Table 1. Furthermore,
regularization becomes an even bigger issue for models with non-linear
calibration operators. Again, stable solutions are obtained after studying the
L-curve and running a few models for a stability check. As can be observed in
Table 1, the ratio of the Frobenius norms of these models is smaller compared
to those obtained with fewer modes and this leads to an increase in
$E_{1}^{c}$. On the other hand, hyper-reduced models perform poorly for both
linear and non-linear calibration coefficients. The performance of the hyper-
reduced linear models is very different compared to the gapless linear models
despite the small difference between error proxies. As shown by the figures,
the performance of the former deteriorates with 12 modes. In a similar
fashion, the hyper-reduced non-linear calibrations models present almost the
same errors as their gapless counterparts. However, hyper-reduced solutions
are unstable and attempting to obtain more accurate ROMs using less intrusive
non-linear coefficients was unsuccessful. In fact, this led to non-linear
models performing similarly to the linear ones.

Hyper-reduction Modes Galerkin Petrov-Galerkin (Grid points) $\rho_{G}^{L}$
$E_{1}^{c}\cdot 10^{7}$ $\rho_{G}^{NL}$ $E_{1}^{c}\cdot 10^{7}$
$\rho_{PG}^{L}$ $E_{1}^{c}\cdot 10^{7}$ $\rho_{PG}^{NL}$ $E_{1}^{c}\cdot
10^{7}$ Yes (1,000) 4 14.058 273.29 228.18 9.7357 8.9025 273.35 144.46 9.7198
8 7.3549 281.78 29.586 17.354 3.9919 281.89 15.255 19.811 12 6.1504 154.67
16.431 32.661 3.4044 154.91 9.2299 32.329 No (168,350) 4 1.3182 260.35 21.7529
9.7341 0.4519 259.34 7.1836 9.7600 8 1.2191 254.04 4.9123 17.143 0.5018 253.86
2.2382 13.860 12 0.9751 139.13 2.5158 32.734 0.3025 138.67 0.7660 32.966

Table 1: Ratio of the Frobenius norms $\rho$ and approximation error
$E_{1}^{c}$ of the ROMs. Here, the superscript specifies if calibration
coefficients are linear (L) or non-linear (NL) and the subscript specifies if
the projection method is Galerkin (G) or Petrov-Galerkin (PG).

(a) Probe in proximity of the leading edge.

(b) Probe in proximity of the trailing edge.

Figure 10: Fluctuation time histories computed by the FOM and gapless ROMs
with 4 POD modes.

(a) Probe in proximity of the leading edge.

(b) Probe in proximity of the trailing edge.

Figure 11: Fluctuation time histories computed by the FOM and hyper-reduced
ROMs with 4 POD modes.

(a) Probe in proximity of the leading edge.

(b) Probe in proximity of the trailing edge.

Figure 12: Fluctuation time histories computed by the FOM and gapless ROMs
with 8 POD modes.

(a) Probe in proximity of the leading edge.

(b) Probe in proximity of the trailing edge.

Figure 13: Fluctuation time histories computed by the FOM and hyper-reduced
ROMs with 8 POD modes.

(a) Probe in proximity of the leading edge.

(b) Probe in proximity of the trailing edge.

Figure 14: Fluctuation time histories computed by the FOM and gapless ROMs
with 12 POD modes.

(a) Probe in proximity of the leading edge.

(b) Probe in proximity of the trailing edge.

Figure 15: Fluctuation time histories computed by the FOM and hyper-reduced
ROMs with 12 POD modes.

Snapshots of u-velocity fluctuations at $t=67.5$ are presented in Figs. 16 and
17. The first figure allows a comparison of results between the FOM and the
different gapless calibrated ROMs with an increasing number of POD modes used
in the reconstruction. Similarly, the second figure can be used to observe
some of the effects of hyper-reduction on the solutions. For this time instant
and, without employing hyper-reduction, the main features of the flow are
recovered when using non-linear calibration coefficients for all cases
considered. The further addition of POD modes to these ROMs improves accuracy
to some degree but not significantly. On the other hand, solutions obtained
with linear calibration coefficients display a poorer performance when only 4
modes are used in the reconstruction. However, results show that they benefit
from a more complete basis of POD modes. In particular, they are capable of
delivering results of similar accuracy compared to those of the non-linear
method. For example, linear and non-linear calibration solutions are similar
when 12 modes are used in the reconstruction with the advantage of linear
models not requiring regularization. On the other hand, performance
deteriorates for most cases when hyper-reduction is applied. Fortunately,
cheap hyper-reduced models with non-linear calibration coefficients show good
results for smaller POD bases (4 and 8 modes). However, adding further modes
to these models is inefficient since they become inaccurate and unstable.
Furthermore, linear methods perform poorly for all hyper-reduced cases
considered.

(a) Linear 4-mode G-ROM.

(b) Linear 8-mode G-ROM.

(c) Linear 12-mode G-ROM.

(d) Non-linear 4-mode G-ROM.

(e) Non-linear 8-mode G-ROM.

(f) Non-linear 12-mode G-ROM.

(g) Linear 4-modes PG-ROM.

(h) Linear 8-modes PG-ROM.

(i) Linear 12-modes PG-ROM.

(j) Non-linear 4-modes PG-ROM.

(k) Non-linear 8-modes PG-ROM.

(l) Non-linear 12-modes PG-ROM.

(m) Full order model.

Figure 16: Contours of u-velocity fluctuations computed at $t=67.5$ using
calibrated gapless models with implicit Euler integration. The terms “G” and
“PG” stand for Galerkin and LSPG models, respectively.

(a) Linear 4-mode G-ROM.

(b) Linear 8-mode G-ROM.

(c) Linear 12-mode G-ROM.

(d) Non-linear 4-mode G-ROM.

(e) Non-linear 8-mode G-ROM.

(f) Non-linear 12-mode G-ROM.

(g) Linear 4-modes PG-ROM.

(h) Linear 8-modes PG-ROM.

(i) Linear 12-modes PG-ROM.

(j) Non-linear 4-modes PG-ROM.

(k) Non-linear 8-modes PG-ROM.

(l) Non-linear 12-modes PG-ROM.

(m) Full order model.

Figure 17: Contours of u-velocity fluctuations computed at $t=67.5$ using
calibrated hyper-reduced models with implicit Euler integration. The terms “G”
and “PG” stand for Galerkin and LSPG models, respectively.

##  4 Conclusions

In this work, a data-driven closure methodology is applied to projection-based
ROMs. The closure approach is based on calibration of model solutions using
the POD temporal modes. The performance of calibration is assessed for both
the Galerkin and least-squares Petrov-Galerkin techniques for two unsteady
compressible flows. The first is the canonical compressible flow past a
cylinder and it serves for testing and verification purposes. The second test
case involves the turbulent flow over a plunging airfoil undergoing deep
dynamic stall. This problem is particularly challenging because of the wide
range of spatial and temporal flow scales.

The ROMs are obtained by projection of the non-conservative form of the
Navier-Stokes equations in order to reduce the polynomial complexity of the
non-linear ODEs arising in the models. This form of the Navier-Stokes
equations facilitates an offline/online procedure where all grid-dependent
calculations are performed only once during the offline step. In order to
achieve further cost reductions in the ROMs’ construction, hyper-reduction is
also applied. The consequences of an additional layer of approximation by an
accelerated greedy MPE method is implemented and evaluated for the deep
dynamic stall problem. Results show that cheap and accurate ROMs can be
obtained when hyper-reduction is combined with non-linear calibration
coefficients and a smaller POD basis. However, for these cases, the
calibration operators are an order of magnitude more intrusive than those
obtained for gapless ROMs.

In the absence of closure, the cylinder flow ROMs display a high amplitude
error. The model solutions are very dissipative, with the LSPG method being
more accurate. This can be explained by the choice of a first-order implicit
Euler time-marching scheme and can be diminished by time step refinement.
Linear and non-linear calibration coefficients are both capable of correcting
the amplitude error without changing the time step, and the resulting
calibrated model solutions are visually indistinguishable from the full order
model. In all calibration cases, resulting operators have relatively low
intrusiveness. Also, non-linear calibration leads to L-curves with clear
corners, resulting in an easy choice of adequate closure operators.

For the dynamic stall problem, non-calibrated ROMs fail completely in
recovering the flow dynamics despite the stable solutions achieved by the
implicit Euler time-marching scheme. Differently from the first problem, non-
linear calibration is far superior to that performed exclusively with linear
coefficients. This is specially true when a smaller POD basis is considered.
The differences between linear and non-linear calibration solutions are
reduced as more POD modes are added to the model bases. This is particularly
observed for gapless ROMs since hyper-reduced calibrated models show stability
issues when higher POD modes are used in their reconstruction. Eventually,
further adding modes makes all models unstable as this type of closure does
not guarantee stability or accuracy for long-term temporal integration. One
should also remind that the present dynamic stall problem is obtained by a
large eddy simulation of a turbulent flow and, in this case, the limited
sampling frequency resolution of snapshots directly impacts the convergence of
higher POD modes, impairing the calibration procedure. For this case,
determining an adequate level of regularization is also more complicated given
that the L-curves do not display clear corners. Despite these issues, present
results show that non-linear calibration with few POD modes can still deliver
stable and accurate solutions for this complex unsteady flow problem.

## Acknowledgments

The authors of this work would like to acknowledge Fundação de Amparo à
Pesquisa do Estado de São Paulo, FAPESP, for supporting the present work under
research grants No. 2013/08293-7, 2018/11410-9 and 2019/18809-7 and Conselho
Nacional de Desenvolvimento Científico e Tecnológico, CNPq, for supporting
this research under grants No. 407842/2018-7 and 304335/2018-5. We also thank
SDUMONT-LNCC (Project SimTurb) and CENAPAD-SP (Project 551) for providing the
computational resources used in this work.

## References

  * [1] T. R. Ricciardi, W. R. Wolf, N. J. Moffitt, J. R. Kreitzman, P. Bent, Numerical noise prediction and source identification of a realistic landing gear, Journal of Sound and Vibration 496 (2021) 115933.  doi:10.1016/j.jsv.2021.115933.
  * [2] L. Sirovich, Turbulence and the dynamics of coherent structures. Part 1: Coherent structures, Quarterly of Applied Mathematics 45 (1987) 561–571.  doi:10.1090/qam/910463.
  * [3] L. Cordier, M. Bergmann, Proper orthogonal decomposition: an overview, in: Lecture series 2002-04, 2003-03 and 2008-01 on post-processing of experimental and numerical data, Von Karman Institute for Fluid Dynamics, 2008., VKI, 2008, p. 46 pages.
  * [4] K. Lee, K. T. Carlberg, Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders, Journal of Computational Physics 404 (2020) 108973.  doi:https://doi.org/10.1016/j.jcp.2019.108973.
  * [5] R. P. Marcondes, T. R. Ricciardi, W. R. Wolf, Spatio-temporal data reconstruction analysis via kernel-based proper orthogonal decomposition, 2021\.  arXiv:https://arc.aiaa.org/doi/pdf/10.2514/6.2021-1096, doi:10.2514/6.2021-1096.
  * [6] C. W. Rowley, T. Colonius, R. M. Murray, Model reduction for compressible flows using POD and Galerkin projection, Physica D: Nonlinear Phenomena 189 (1–2) (2004) 115 – 129.
  * [7] K. Carlberg, C. Farhat, D. Amsallem, The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and tubulent flows, Journal of Computational Physics 243 (2013) 623–647.  doi:10.1016/j.jcp.2013.02.028.
  * [8] S. J. Grimberg, C. Farhat, N. Youkilis, On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows, Journal of Computational Physics 419 (2020) 1–28.  doi:10.1016/j.jcp.2020.109681.
  * [9] H. F. S. Lui, W. R. Wolf, Construction of reduced-order models for fluid flows using deep feedforward neural networks, Journal of Fluid Mechanics 872 (2019) 963–994.  doi:10.1017/jfm.2019.358.
  * [10] S. L. Brunton, J. L. Proctor, N. J. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937.  arXiv:http://www.pnas.org/content/113/15/3932.full.pdf, doi:10.1073/pnas.1517384113.
  * [11] V. Zucatti, H. F. S. Lui, D. B. Pitz, W. R. Wolf, Assessment of reduced-order modeling strategies for convective heat transfer, Numerical Heat Transfer, Part A: Applications 77 (2020) 702–729.
  * [12] W. Cazemier, R. W. C. P. Verstappen, A. E. P. Veldman, Proper orthogonal decomposition and low-dimensional models for driven cavity flows, Physics of Fluids 10 (7) (1998) 1685–1699.  doi:10.1063/1.869686.
  * [13] B. R. Noack, P. Papas, P. A. Monkewitz, The need for a pressure-term representation in empirical Galerkin models of incompressible shear flows, Journal of Fluid Mechanics 523 (2005) 339–365.  doi:10.1017/S0022112004002149.
  * [14] S. Grimberg, C. Farhat, R. Tezaur, 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 (2021) 1–29arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6603, doi:10.1002/nme.6603.
  * [15] V. Zucatti, W. Wolf, M. Bergmann, Calibration of projection-based reduced-order models for unsteady compressible flows, Journal of Computational Physics (2021) 110196doi:10.1016/j.jcp.2021.110196.
  * [16] M. Bergmann, C. H. Bruneau, A. Iollo, Enablers for robust POD models, Journal of Computational Physics 228 (2009) 516–538.  doi:10.1016/j.jcp.2008.09.024.
  * [17] O. San, T. Iliescu, Proper orthogonal decomposition closure models for fluid flows: Burgers equation, International Journal of Numerical Analysis and Modeling, Series B 5 (2014) 217–237.
  * [18] Z. Wang, I. Akhtar, J. Borggaard, T. Iliescu, Proper orthogonal decomposition closure models for turbulent flows: A numerical comparison, Comput. Methods Appl. Mech. Eng. 237 (2012) 10–26.
  * [19] S. Giere, T. Iliescu, V. John, D. Wells, Supg reduced order models for convection-dominated convection–diffusion–reaction equations, Computer Methods in Applied Mechanics and Engineering 289 (2015) 454–474.  doi:10.1016/j.cma.2015.01.020.
  * [20] B. Galletti, A. Bottaro, C. H. Bruneau, A. Iollo, Accurate model reduction of transient and forced wakes, European Journal of Mechanics B/Fluids 26 (2006) 354–366.  doi:10.1016/j.euromechflu.2006.09.004.
  * [21] M. Couplet, C. Basdevant, P. Sagaut, Calibrated reduced-order POD-Galerkin system for fluid flow modelling, Journal of Computational Physics 207 (2005) 192–220.  doi:0.1016/j.jcp.2005.01.008.
  * [22] R. Bourguet, M. Braza, A. Dervieux, Reduced-order modeling for unsteady transonic flows around an airfoil, Physics of Fluids 19 (2007) 111701.  doi:10.1063/1.2800042.
  * [23] R. Bourguet, M. Braza, A. Sévrain, A. Bouhadji, Capturing transition features around a wing by reduced-order modeling based on compressible Navier–Stokes equations, Physics of Fluids 21 (2009) 094104.  doi:10.1063/1.3234398.
  * [24] M. Balajewicz, I. Tezaur, E. Dowell, Minimal subspace rotation on the Stiefel manifold for stabilization and enhancement of projection-based reduced order models for the compressible Navier-Stokes equations, Journal of Computational Physics 321 (2016) 224–241.  doi:10.1016/j.jcp.2016.05.037.
  * [25] S. J. Grimberg, C. Farhat, Hyperreduction of cfd models of turbulent flows using a machine learning approach, in: AIAA Scitech Forum, 2020.  doi:10.2514/6.2020-0363.
  * [26] C. Mou, B. Koc, O. San, L. G. Rebholz, T. Iliescu, Data-driven variational multiscale reduced order models, Computer Methods in Applied Mechanics and Engineering 373 (2021) 113470.  doi:https://doi.org/10.1016/j.cma.2020.113470.
  * [27] X. Xie, M. Mohebujjaman, L. G. Rebholz, T. Iliescu, Data-driven filtered reduced order modeling of fluid flows, SIAM Journal on Scientific Computing 40 (3) (2018) 834–857.  arXiv:https://doi.org/10.1137/17M1145136, doi:10.1137/17M1145136.
  * [28] M. Mohebujjaman, L. Rebholz, T. Iliescu, Physically constrained data-driven correction for reduced-order modeling of fluid flows, International Journal for Numerical Methods in Fluids 89 (3) (2019) 103–122.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.4684, doi:10.1002/fld.4684.
  * [29] L. Cordier, B. Abou El Majd, J. Favier, Calibration of pod reduced-order models using tikhonov regularization, International Journal for Numerical Methods in Fluids 63 (2) (2010) 269–296.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.2074, doi:10.1002/fld.2074.
  * [30] B. Abou El Majd, L. Cordier, New regularization method for calibrated POD reduced-order models, Mathematical Modelling and Analysis 21 (2016) 47–62.  doi:10.3846/13926292.2016.1132486.
  * [31] J. Favier, L. Cordier, A. Kourta, A. Iollo, Calibrated POD reduced-order models of massively separated flows in the perspective of their control, Fluids Engineering Division Summer Meeting, 2006, pp. 743–748.  doi:10.1115/FEDSM2006-98495.
  * [32] V. Zucatti, W. Wolf, Assessment of projection-based reduced-order modeling strategies for unsteady flows, in: AIAA Scitech Forum, 2021.  arXiv:https://arc.aiaa.org/doi/pdf/10.2514/6.2021-0241, doi:10.2514/6.2021-0241.
  * [33] C. Mou, Z. Wang, D. R. Wells, X. Xie, T. Iliescu, Reduced order models for the quasi-geostrophic equations: A brief survey, Fluids 6 (1) (2021).  doi:10.3390/fluids6010016.
  * [34] C. Mou, Cross-validation of data-driven correction reduced order modeling, Ph.D. thesis, Virginia Polytechnic Institute and State University (2018).
  * [35] X. Xie, C. Webster, T. Iliescu, Closure learning for nonlinear model reduction using deep residual neural network, Fluids 5 (1) (2020).  doi:10.3390/fluids5010039.
  * [36] M. Bergmann, L. Cordier, J.-P. Brancher, Optimal rotary control of the cylinder wake using proper orthogonal decomposition reduced-order model, Physics of Fluids 17 (2005) 097101.  doi:10.1063/1.2033624.
  * [37] K. Carlberg, M. Barone, H. Antil, Galerkin v. least-squares Petrov-Galerkin projection in nonlinear model reduction, Journal of Computational Physics 330 (2017) 693–734.  doi:10.1016/j.jcp.2016.10.033.
  * [38] A. M. A. Quarteroni, F. Negri, Reduced Basis Methods for Partial Differential Equations: An Introduction, 1st Edition, Springer, Switzerland, 2016\.
  * [39] J. Nocedal, S. J. Wright, Numerical Optimization, Springer, 2006.  doi:10.1007/978-0-387-40065-5.
  * [40] S. Grimberg, Projection-based model order reduction and hyperreduction of turbulent flow models, Ph.D. thesis, Stanford University (2020).
  * [41] M. Zahr, Adaptive model reduction to accelerate optimization problems governed by partial differential equations, Ph.D. thesis, Stanford University (2016).
  * [42] K. Willcox, Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition, Computers & Fluids 35 (2006) 208–226.  doi:10.1016/j.compfluid.2004.11.006.
  * [43] R. Everson, L. Sirovich, Karhunun-Loève procedure for gappy data, Optical Society of America 12 (1995) 1657–1664.  doi:10.1364/JOSAA.12.001657.
  * [44] P. Astrid, S. Weiland, K. Willcox, T. Backx, Missing point estimation in models described by proper orthogonal decomposition, IEEE Transactions on Automatic Control 53 (2008) 2237–2251.  doi:10.1109/TAC.2008.2006102.
  * [45] R. Zimmermann, K. Willcox, An accelerated greedy missing point estimation procedure, SIAM Journal on Scientific Computing 38 (2016) A2827–A2850.  doi:10.1137/15M1042899.
  * [46] A. Iollo, S. Lanteri, J.-A. Désidéri, Stability properties of POD-Galerkin approximations for the compressible Navier-Stokes equations, Theoretical and Computational Fluid Dynamics 13 (2000) 377–396.
  * [47] S. Chaturantabut, D. C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal on Scientific Computing 32 (2010) 2737–3764.  doi:10.1137/090766498.
  * [48] B. Peherstorfer, K. Willcox, Online adaptive model reduction for nonlinear systems via low-rank updates, SIAM J. Sci. Comput. 37 (2015) A2123–A2150.  doi:10.1137/140989169.
  * [49] M. Yano, Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws, Adv. Comput. Math. 45 (2019) 2287–2320.
  * [50] J. Chan, Entropy stable reduced order modeling of nonlinear conservation laws, Journal of Computational Physics 423 (2020) 109789.  doi:10.1016/j.jcp.2020.109789.
  * [51] C. Farhat, P. Avery, T. Chapman, 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 (9) (2014) 625–662.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.4668, doi:10.1002/nme.4668.
  * [52] C. Mou, H. Liu, D. R. Wells, T. Iliescu, Data-driven correction reduced order models for the quasi-geostrophic equations: a numerical investigation, International Journal of Computational Fluid Dynamics 34 (2) (2020) 147–159.  arXiv:https://doi.org/10.1080/10618562.2020.1723556, doi:10.1080/10618562.2020.1723556.
  * [53] R. Bourguet, Analyse physique et modélisation d’écoulements turbulents instationnaires compressibles autour de surfaces portantes par approches statistiques haute-fidélité et de dimension réduite dans le contexte de l’interaction fluide-structure, Ph.D. thesis, l’Institut National Polytechnique de Toulouse (2008).
  * [54] B. L. O. Ramos, W. R. Wolf, C.-A. Yeh, K. Taira, Active flow control for drag reduction of a plunging airfoil under deep dynamic stall, Phys. Rev. Fluids 4 (2019) 074603.  doi:10.1103/PhysRevFluids.4.074603.