# A Reduced Basis Technique for Long-Time Unsteady Turbulent Flows

###### Abstract

We present a reduced basis technique for long-time integration of parametrized
incompressible turbulent flows. The new contributions are threefold. First, we
propose a constrained Galerkin formulation that corrects the standard Galerkin
statement by incorporating prior information about the long-time attractor.
For explicit and semi-implicit time discretizations, our statement reads as a
constrained quadratic programming problem where the objective function is the
Euclidean norm of the error in the reduced Galerkin (algebraic) formulation,
while the constraints correspond to bounds for the maximum and minimum value
of the coefficients of the $N$-term expansion. Second, we propose an _a
posteriori_ error indicator, which corresponds to the dual norm of the
residual associated with the time-averaged momentum equation. We demonstrate
that the error indicator is highly-correlated with the error in mean flow
prediction, and can be efficiently computed through an offline/online
strategy. Third, we propose a Greedy algorithm for the construction of an
approximation space/procedure valid over a range of parameters; the Greedy is
informed by the _a posteriori_ error indicator developed in this paper. We
illustrate our approach and we demonstrate its effectiveness by studying the
dependence of a two-dimensional turbulent lid-driven cavity flow on the
Reynolds number.

Lambert Fick1, Yvon Maday2,3, Anthony T Patera4, Tommaso Taddei2

1 Texas A&M University, Department of Nuclear Engineering, USA
<EMAIL_ADDRESS>
2 Sorbonne Universités, Laboratoire Jacques-Louis Lions, France
<EMAIL_ADDRESS><EMAIL_ADDRESS>
3 Brown University, Division of Applied Mathematics, USA
<EMAIL_ADDRESS>
4 MIT, Department of Mechanical Engineering, USA
<EMAIL_ADDRESS>

_Keywords:_ model order reduction, reduced basis method, CFD, _a posteriori_
error estimation

##  1 Introduction

For turbulent flows, estimation of the entire solution trajectory through a
low-dimensional Reduced Order Model (ROM) is infeasible due to the slow decay
of the Kolmogorov $N$-width, and due to the sensitivity of the dynamical
system to perturbations. Nevertheless, it might still be possible to estimate
various moments of the solution associated to a Direct Numerical Simulation
(DNS). The goal of this work is to develop a Reduced Basis (RB) technique for
long-time integration of turbulent flows. Our equations of interest are the
unsteady incompressible Navier-Stokes equations for high-Reynolds number flows
with no-slip boundary conditions:

$\left\\{\begin{array}[]{ll}\partial_{t}u+(u\cdot\nabla)u-\frac{1}{{\rm
Re}}\Delta u+\nabla p=f&{\rm in}\,\Omega\times\mathbb{R}_{+},\\\\[8.53581pt]
\nabla\cdot u=0&{\rm in}\,\Omega\times\mathbb{R}_{+},\\\\[8.53581pt] u=g&{\rm
on}\,\partial\Omega\times\mathbb{R}_{+},\\\\[8.53581pt] u=u_{0}&{\rm
on}\,\Omega\times\\{0\\},\\\\[8.53581pt] \end{array}\right.$ (1)

where $\Omega\subset\mathbb{R}^{d}$, and $f,g,u_{0}$ are suitable fields. We
denote by $\mu\in\mathcal{P}\subset\mathbb{R}^{P}$ the set of parameters
associated with the equations.

We consider two separate problems: the solution reproduction problem, and the
parametric problem. In the _solution reproduction problem_ , given the
velocity DNS data $\\{u(t_{\rm s}^{k};\mu)\\}_{k=1}^{K}$, and possibly the
pressure DNS data $\\{p(t_{\rm s}^{k};\mu)\\}_{k=1}^{K}$, at the sampling
times $\\{t_{\rm s}^{k}\\}_{k=1}^{K}\subset\mathbb{R}_{+}$, we wish to
construct a ROM that approximates — in a sense that will be defined soon — the
original DNS data for the same value of the parameter $\mu$. In the
_parametric problem_ , we wish to construct a ROM that approximates the DNS
data for all values of $\mu$ in a prescribed parameter range
$\mathcal{P}\subset\mathbb{R}^{P}$. For the parametric problem, we wish to
control the offline costs associated with the construction of the reduced
space: this implicates a Greedy (rather than POD) strategy in parameter.
Although the solution reproduction problem might be of limited interest in
practice, it represents the first step towards the development of a ROM for
the parametric problem.

Following [1], we quantify the accuracy of the ROM by computing the error in
the long-time average $\langle
u\rangle(x):=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}u(x,t)\,dt$, and the
error in the turbulent kinetic energy ${\rm
TKE}(t)=\frac{1}{2}\int_{\Omega}\|u(x,t)-\langle u\rangle(x)\|_{2}^{2}\,dx$
where $\|\cdot\|_{2}$ is the Euclidean norm. We remark that, at present, there
is no universally-accepted notion of ROM accuracy for turbulent flows. In [2],
the authors evaluate different ROMs based on five different criteria: the
kinetic energy spectrum, the mean velocity components, the Reynolds stresses,
the root mean square values of the velocity fluctuations, and the time
evaluations of the POD coefficients. In [3], the authors consider just the
time evaluations and the power spectra of selected POD coefficients. From an
engineering perspective, the definition of accuracy is entirely determined by
the particular quantity of interest we wish to predict: for this reason, we
envision that for several applications accurate estimates of the long-time
averages and possibly of the turbulent kinetic energy might suffice.

The most popular approach for the solution reproduction problem is the so-
called POD-Galerkin method [4, 5, 6, 7]: first, we generate a reduced space
$\mathcal{Z}^{\rm u}=\\{\zeta_{n}\\}_{n=1}^{N}$ for the velocity field by
applying the Proper Orthogonal Decomposition (POD, [8, 9, 10]) in the $L^{2}$
inner product; then, we estimate the velocity field for each time-step $t_{\rm
g}^{1},\ldots,t_{\rm g}^{J}$ as
$\hat{u}^{j}(\cdot)=\sum_{n=1}^{N}a_{n}^{j}\,\zeta_{n}(\cdot)$ where the
coefficients $\mathbf{a}^{j}=[a_{1}^{j},\ldots,a_{N}^{j}]$ are computed by
projecting the momentum equation onto the space $\mathcal{Z}^{\rm u}$. Since
all DNS data for the velocity are divergence-free, it is straighforward to
verify that POD modes $\zeta_{1},\ldots,\zeta_{N}$ are divergence-free, and so
is the ROM solution.

As observed by several authors, ROMs based on $L^{2}$ POD-Galerkin are prone
to instabilities [11, 12]. This can be explained through a physical argument.
In the limit of high-Reynolds numbers, large-scale flow features are broken
down into smaller and smaller scales until the scales are fine enough that
viscous forces can dissipate their energy ([13]). This implies that small-
scale modes have significant influence on the dynamics. POD modes based on the
$L^{2}$ inner product are biased toward large, high-energy scales: since large
scales are not endowed with the natural energy dissipation tendency of the
smaller lower-energy viscous scales, this leads to instabilities and/or large
errors in the estimate of the turbulent kinetic energy.

To address the issue of stability, several strategies have been proposed: (i)
including dissipation via a closure model, (ii) modifying the POD basis by
including functions that resolve a range of scales, (iii) employing a minimum
residual formulation, (iv) employing stabilizing inner products, (v)
calibration methods, and (vi) generating the reduced space through Dynamic
Mode Decomposition. Below, we briefly describe each strategy, and we provide
some references; we remark that most of the works presented below are
restricted to either laminar flows or short-time integration; therefore, they
do not directly address the problem of interest (the long-time integration of
fully-turbulent flows). We also recall that other topics are treated in the
literature: in particular, Noack et al. [14] proposed to incorporate pressure
in the ROM for cases with other than no-slip boundary conditions. For the
problem considered in this paper, the discussion of [14] is not relevant, and
is here omitted.

  * (i)

Starting with the pioneering work in [15], several authors have proposed to
include dissipation through the vehicle of a closure model. A first class of
models is designed and motivated by analogy with Large Eddy Simulation (see
[16] for an introduction to LES): in this respect, Couplet et al. ([17]) and
Noack et al. ([18, section 4.2]) provide theoretical and numerical evidence
that the energy transfer among $L^{2}$-based POD modes is similar to the
energy transfer among Fourier modes, and for this reason LES ideas based on
the energy cascade concept might be promising for POD-ROMs. We remark that in
[17, 18] the POD space is built for the fluctuating field $u_{\rm f}=u-\langle
u\rangle$ for a fixed value of the parameters (solution reproduction problem)
based on the $L^{2}$ inner product: for this reason, it appears difficult to
rigorously apply these ideas in the parametric setting in which we must
combine modes associated with different parameters. Another class of closure
models is based on the extension of stabilization techniques originally
introduced in the Finite Element or Spectral framework: two notable examples
are the Spectral Vanishing Viscosity Model (SVVM, [19], see also [20])
originally presented by Tadmor in [21] for spectral discretization of
nonlinear conservation laws for controlling high-wave number oscillations, and
the SUPG stabilization discussed in [22]. We refer to [2] for a numerical
comparison of four closure models for incompressible Navier-Stokes equations:
the mixing-length model, the Smagorinsky model, the variational multiscale
model, and the dynamic subgrid-scale model. We further refer to [23] for
another POD closure model based on approximate deconvolution, and we refer to
[24] for a numerical comparison of several closure models for the Burgers’
equation. Finally, we mention the nonlinear Galerkin method proposed by Marion
and Temam in [25], and applied to the simulation of turbulent flows in [26].
As the above-mentioned variational multiscale method, this approach corrects
the standard Galerkin model by exploiting the separation between large-scale
and small-scale modes. To our knowledge, the nonlinear Galerkin method has
never been applied in the model order reduction framework.

  * (ii)

Another approach is based on including in the POD basis functions that resolve
a range of scales. Bergmann et al. in [27] (see also [28]) proposed to augment
the original POD basis with a second POD performed on the residuals of the
momentum equation (and of the mass equation in case pressure is modelled by
the ROM); on the other hand, Balajewicz and Dowell proposed a Greedy technique
to include in the basis random linear combinations of low-energy POD modes —
associated with the $L^{2}$ inner product.

  * (iii)

Minimum residual formulation was first introduced in the reduced basis
framework in [29] for linear noncoercive problems, and then extended to fluid
problems in [30, 31, 11, 32]. Given the reduced space $\mathcal{Z}^{\rm u}$
for velocity (and possibly the reduced space $\mathcal{Z}^{\rm p}$ for
pressure), after having discretized the equation in time, the latter approach
computes the solution in $\mathcal{Z}^{\rm u}$ (or $\mathcal{Z}^{\rm
u}\times\mathcal{Z}^{\rm p}$) that minimizes a suitable dual residual at each
time-step. We remark that for problems with quadratic nonlinearities minimum
residual ROMs require $\mathcal{O}(N^{4})$ storage and the online cost for
each time-step is $\mathcal{O}(N^{4})$ for semi-implicit/explicit time-
discretizations — as opposed to $\mathcal{O}(N^{3})$ for standard POD-Galerkin
ROMs. For this reason, hyper-reduction techniques are employed to reduce the
online cost and the memory constraints ([30, 31, 11]).

  * (iv)

Iollo et al. in [33] proposed to employ the $H^{1}$ inner product — rather
than the more standard $L^{2}$ inner product — to generate the POD modes. This
choice is motivated by dynamic considerations: since small-scale modes have
relatively large $H^{1}$ norm compared to their $L^{2}$ norm, and recalling
that small scales are responsible for energy dissipation, the use of the
$H^{1}$ inner product leads to a more dissipative reduced order model. We
remark that several other authors proposed to not employ the standard $L^{2}$
inner product ([34, 35, 36]); however, their choices were not motivated by
long-time stability considerations.

  * (v)

If we denote by $\dot{\mathbf{a}}=\mathcal{F}(\mathbf{a})$ the ROM for the
coefficients of the POD expansion, in [37], Couplet et al. proposed to
calibrate the coefficients of $\mathcal{F}$ based on DNS data, under the
assumption that $\mathcal{F}$ is a polynomial of degree $2$ in $\mathbf{a}$.
We observe that the ROM $\mathcal{F}$ depends on the particular POD basis
selected; for large values of $N$ (dimension of the POD space), the
calibration procedure might require a substantial number of DNS snapshots.

  * (vi)

Dynamic Mode Decomposition (DMD) was first proposed by Schmid in [38]; as
shown by Rowley et al. in [39] DMD can be interpreted as an algorithm for
finding the Koopman modes associated with the nonlinear discrete dynamical
system obtained from the discretization of the Navier-Stokes equations.
Despite several authors have proven the effectiveness of DMD for the
extraction of physically-relevant time scales and their associated spatial
structures ([39, 38, 40]), the work by Alla and Kutz [41] represents one of
the few examples of application of DMD within the Galerkin framework.

Despite these advances, the solution reproduction problem remains an open
issue, particularly for turbulent flows. By performing a detailed analysis of
the performance of the POD-Galerkin approach, we empirically demonstrate that
in the case of turbulent flows POD-Galerkin ROMs might exhibit other spurious
effects such as false stable steady flows. This demonstrates the need for a
more fundamental correction to the POD-Galerkin formulation. We remark that a
similar issue has been observed in [19] by Sirisup and Karniadakis for long-
time integration of a POD-Galerkin ROM for a laminar flow past a cylinder, and
— in a different context — by Curry et al. in [42] for highly-truncated
spectral approximations to turbulent flows.

To our knowledge, there are very few works that systematically address the
parametric problem. Ma and Karniadakis ([6]), Galletti et al. ([7]), and
Stabile et al. ([43]) developed a reduced order model based on POD-Galerkin
for the flow past a cylinder for a wide range of Reynolds numbers in the
laminar regime. In these papers, the authors use DNS data for pre-selected
Reynolds numbers to generate reduced spaces for velocity ([6, 7]), and for
velocity and pressure ([43]). The choice of the parameters for which the DNS
data are computed is performed _a priori_. Non-adaptive explorations of the
parameter space typically require a large number of offline evaluations of the
Full Order Model (FOM); for this reason, they might not be practical in our
context.

The goal of this work is to develop a Model Order Reduction (MOR) procedure
for the parametrized incompressible Navier-Stokes equations. The three key
pieces of our MOR technique are (i) a reduced formulation for the computation
of the reduced-order solution, (ii) an _a posteriori_ indicator for the error
in the prediction of the mean flow, and (iii) a $H^{1}$-POD-$h$Greedy strategy
for the construction of the reduced space informed by the above-mentioned _a
posteriori_ indicator.

  * (i)

Our reduced formulation is based on a constrained Galerkin formulation. The
approach is designed to correct the standard Galerkin formulation, especially
for moderate values of $N$. For explicit and semi-implicit time
discretizations the formulation reads as a quadratic programming problem where
the objective function corresponds to the Euclidean norm of the error in the
reduced Galerkin (algebraic) formulation, while the constraints correspond to
bounds for the maximum and minimum value of the coefficients
$\\{a_{n}^{j}\\}_{n=1}^{N}\subset\mathbb{R}$ of the expansion. We discuss an
actionable procedure to estimate the lower and upper bounds associated with
each coefficient of the reduced expansion based on DNS data.

  * (ii)

Our error indicator corresponds to the dual norm of the residual associated
with the time-averaged momentum equation. Time-averaging is here motivated by
the chaotic behavior in time of the velocity field. We verify that the error
indicator can be efficiently computed through an offline/online strategy;
furthermore, we numerically demonstrate that the indicator is highly-
correlated with the error in the mean flow prediction: therefore, it is well-
suited to drive the Greedy procedure for the generation of the ROM.

  * (iii)

As in the seminal work by Haasdonk and Ohlberger [44], our POD- $h$Greedy
algorithm combines POD in time with Greedy in parameter. The procedure is a
simplified version of the $h$-type Greedy proposed in [45]. Given
$\mu^{1}\in\mathcal{P}$, we generate the DNS data for $\mu^{1}$, we apply POD
— based on the $H^{1}$ inner product — to generate the reduced space
$\mathcal{Z}_{1}^{\rm u}$, we build the POD-ROM, and we evaluate the error
indicator $\Delta_{1}^{\rm u}(\mu)$ for all $\mu\in\mathcal{P}_{\rm
train}\subset\mathcal{P}$. Then, we select $\mu^{2}$ that maximizes the error
estimate $\Delta_{1}^{\rm u}$ over the training set $\mathcal{P}_{\rm train}$.
During the second iteration, we perform the same steps as before for $\mu^{2}$
(generation of DNS data, POD, construction of the ROM, estimate of the error).
Finally, we select $\mu^{3}$ that maximizes $\Delta_{1,2}^{\rm
u}(\mu):=\min\\{\Delta_{1}^{\rm u}(\mu),\Delta_{2}^{\rm u}(\mu)\\}$ over
$\mathcal{P}_{\rm train}$. We then proceed to generate
$\mu^{4},\ldots,\mu^{L}$. At the end of the offline stage, the procedure
produces $L$ different ROMs; during the online stage, given a new value of
$\mu\in\mathcal{P}$, we first evaluate the ROMs associated with the $n_{\rm
cand}$ nearest anchor points, and then we select the ROM that minimizes the
error indicator.

We observe that in this work we restrict ourselves to linear approximation
spaces $\mathcal{Z}^{\rm u}$ that do not depend on time: this greatly
simplifies the implementation, and reduces the memory constraints for long-
time integration. We refer to [46, 47, 48] for MOR strategies based on
nonlinear approximation spaces for unsteady problems. On the other hand, we
refer to [49, 50] for space-time approximations of linear and nonlinear
parabolic problems.

The idea of employing a constrained formulation is new in the MOR framework.
We observe that a constrained formulation has been recently proposed in [51]
in the context of steady-state data assimilation: as in our work, the
constraints in [51] provide further information about the solution manifold;
however, while in our work the constraints are designed to compensate for the
effect of the unmodelled dynamics, in [51] the constraints are designed to
limit the effect of experimental noise. As opposed to calibration techniques
and also stabilized ROMs, the hyper-parameters of the ROM (the lower and upper
bounds for the coefficients of the expansion) are here tuned directly through
sparse DNS data, for an arbitrary reduced space $\mathcal{Z}^{\rm u}$, without
having to evaluate the ROM for several tentative candidates. This feature of
the approach greatly simplifies the implementation of the method, and in
practice reduces the offline costs.

The time-averaged error indicator is also new. In [44], the authors employ a
residual estimator that measures the error in the entire trajectory: for
turbulent flows, this metric is not appropriate due to the chaotic nature of
the dynamical system. This explains the importance of our new error indicator
for the problem at hand.

The POD-Greedy algorithm was first proposed in [44], and then analyzed in
[52]. The algorithm in [44] combines data from different parameters to
generate a single reduced space for the entire parameter space $\mathcal{P}$.
On the other hand, in our approach we build a reduced space for each of them.
Recalling the definitions of [45], the algorithm of [44] corresponds to a
POD-$p$Greedy, while our approach corresponds to a POD-$h$Greedy. For
turbulent flows, we empirically show in Appendix E that combining modes
associated with different values of the parameters might lead to poor
performance. On the other hand, $h$-refinement leads to more accurate and
stable ROMs.

The paper is organized as follows. In section 2, we introduce the model
problem considered in this work. In section 3, we consider the solution
reproduction problem. First, we consider the POD-Galerkin approach: we
introduce the formulation, and we assess the numerical performance. Then, we
present our constrained POD-Galerkin approach: as in the previous case, we
discuss the formulation, and then we numerically assess the performance. In
section 4, we consider the parametric problem: first, we present the
POD-$h$Greedy approach; second, we discuss how to adapt the constrained
Galerkin formulation to the parametric setting; third, we propose the time-
averaged error indicator; and fourth, we present the numerical assessment. In
section 5, we offer some concluding remarks, and we discuss potential
extensions of the current approach. A number of appendices provide further
analysis and numerical investigations: in Appendix A we provide an analysis of
the model problem considered; in Appendix B we discuss the selection of the
sampling times $\\{t_{\rm s}^{k}\\}_{k=1}^{K}$; in Appendix C we propose a
suitable definition of stability for ROMs; in Appendix D we investigate the
robustness of the constrained formulation proposed in this paper; in Appendix
E we illustrate the problem of $p$-refinement for the parametric case; and in
Appendix F we describe the offline/online strategy employed to compute the
error indicator.

##  2 A lid-driven cavity problem

We consider the following unsteady lid-driven cavity problem:

$\left\\{\begin{array}[]{ll}\partial_{t}u+(u\cdot\nabla)u-\nu({\rm Re})\Delta
u+\nabla p=0&{\rm in}\,\Omega\times\mathbb{R}_{+},\\\\[8.53581pt] \nabla\cdot
u=0&{\rm in}\,\Omega\times\mathbb{R}_{+},\\\\[8.53581pt] u=g(x)&{\rm
on}\,\Gamma_{\rm top}\times\mathbb{R}_{+},\\\\[8.53581pt] u=0&{\rm
on}\,\partial\Omega\setminus\Gamma_{\rm
top}\times\mathbb{R}_{+},\\\\[8.53581pt] u=0&{\rm
on}\,\Omega\times\\{0\\},\\\\[8.53581pt] \end{array}\right.$ (2a) where the
velocity $u:\Omega\times\mathbb{R}_{+}\to\mathbb{R}^{2}$ is a two-dimensional
vector field, the pressure $p:\Omega\times\mathbb{R}_{+}\to\mathbb{R}$ is a
scalar field, $\nu({\rm Re})=\frac{1}{{\rm Re}}$, $\Omega=(-1,1)^{2}$,
$\Gamma_{\rm top}=\\{x\in\bar{\Omega}:x_{2}=1\\}$, the Dirichlet datum is
given by $g(x)=\left[\begin{array}[]{c}(1-x_{1}^{2})^{2}\\\ 0\\\
\end{array}\right],$ (2b) and the Laplacian $\Delta$ should be interpreted as
component-wise. We remark that in (2) time is non-dimensionalized by the
convective scaling (i.e., dimensional boxside half-length divided by
dimensional maximum lid velocity).

The problem corresponds to a isothermal, incompressible, two-dimensional flow
inside a square cavity driven by a prescribed lid velocity. The problem is a
well-known prototypical example used to validate numerical schemes and reduced
order models ([1, 53, 54, 55, 56, 57]); unlike in the more standard lid-driven
cavity problem with $g(x)=[1,0]$, here we regularize the singularity near the
upper corners of the cavity.

In this paper, we study the dependence of the flow on the Reynolds number,
that is $\mu={\rm Re}$. It is well-known ([55]) that the flow exhibits a long-
time unsteady but stationary solution for ${\rm Re}>{\rm Re_{c}}$ ([55]); here
stationarity implies that all statistics are invariant under a shift in time
([58]). Since we are interested in long-time unsteady flows, we here consider
${\rm Re}\in\mathcal{P}=[{\rm Re}_{\rm LB},{\rm Re}_{\rm UB}]=[15000,25000]$:
for all values of ${\rm Re}$ in $\mathcal{P}$ the flow is asymptotically
statistically stationary. Balajewicz and Dowell considered the same problem —
for a single value of $\rm Re$ — in [1]; we remark that they define the
viscosity as $\nu({\rm Re})=\frac{2}{{\rm Re}}$, and they consider the case
${\rm Re}=30000$.

In view of the development of the ROM for (2) it is convenient to consider the
lifted equations. If we denote by $R_{g}$ the two-dimensional vector field
defined as the solution to the following Stokes problem:

$\left\\{\begin{array}[]{ll}-\Delta R_{g}+\nabla\lambda=0&{\rm
in}\,\Omega,\\\\[8.53581pt] \nabla\cdot R_{g}=0&{\rm
in}\,\Omega,\\\\[8.53581pt] R_{g}=g&{\rm on}\,\Gamma_{\rm top},\\\\[8.53581pt]
R_{g}=0&{\rm on}\,\partial\Omega\setminus\Gamma_{\rm top},\\\
\end{array}\right.$ (3a) we can define the lifted velocity solution
$\mathring{u}=u-R_{g}$ as the solution to:
$\left\\{\begin{array}[]{ll}\partial_{t}\mathring{u}+\left((\mathring{u}+R_{g})\cdot\nabla\right)(\mathring{u}+R_{g})-\frac{1}{\rm
Re}\Delta\left(\mathring{u}+R_{g}\right)+\nabla p=0&{\rm
in}\,\Omega\times\mathbb{R}_{+},\\\\[8.53581pt] \nabla\cdot\mathring{u}=0&{\rm
in}\,\Omega\times\mathbb{R}_{+},\\\\[8.53581pt] \mathring{u}=0&{\rm
on}\,\partial\Omega\times\mathbb{R}_{+},\\\\[8.53581pt]
\mathring{u}(t=0)=-R_{g}&{\rm on}\,\Omega\times\\{0\\}.\\\\[8.53581pt]
\end{array}\right.$ (3b)

Then, if we introduce the spaces $V:=[H_{0}^{1}(\Omega)]^{2}$, and $Q=\\{q\in
L^{2}(\Omega):\int_{\Omega}q\,dx=0\\}$, we can define the weak form associated
with (3): find $(\mathring{u},p)\in\mathcal{V}\times\mathcal{Q}$ such that for
a.e. $t>0$

$\left\\{\begin{array}[]{l}\langle\partial_{t}\mathring{u}(t),v\rangle_{\star}+\frac{1}{\rm
Re}(\mathring{u}(t)+R_{g},v)_{V}+c(\mathring{u}(t)+R_{g},\mathring{u}(t)+R_{g},v)+b(v,p(t))=0\\\\[5.69054pt]
\hfill\forall\,v\in V,\\\\[8.53581pt]
b(\mathring{u}(t),q)=0\hfill\forall\,q\in Q,\\\ \end{array}\right.$ (4a) where
$\mathcal{V}=\\{v\in L_{\rm loc}^{2}(\mathbb{R}_{+};V):\,\partial_{t}v\in
L_{\rm loc}^{2}(\mathbb{R}_{+};V^{\prime})\\}$, $\mathcal{Q}=L_{\rm
loc}^{2}(\mathbb{R}_{+};Q)$, $\langle\cdot,\cdot\rangle_{\star}$ denotes the
pairing between $V^{\prime}$ and $V$ which (for our smoothness assumptions and
numerical approximations) can be evaluated in terms of the pivot space
$L^{2}$, $(w,v)_{V}=\int_{\Omega}\,\nabla\,w:\nabla v\,dx$ is the inner
product associated with $V$, and $c(w,u,v)=\int_{\Omega}\,(w\cdot\nabla)u\cdot
v\,dx,\quad b(v,q)=-\int_{\Omega}\,(\nabla\cdot v)q\,dx.$ (4b)

We resort to a $\mathbb{Q}_{M}-\mathbb{Q}_{M-2}$ spectral element ([59])
discretization in space, and to an explicit three-step Adams-Bashforth (AB3)/
implicit two-step Adams-Moulton (AM2) discretization in time. DNS simulations
are performed using the open-source software nek5000 ([60]). We refer to the
spectral element literature (see, e.g., [61, 62, 63, 64]) for further details
about the spectral element method and its implementation for fluid dynamics
problems. More in detail, we consider a $16$ by $16$ structured quadrilateral
mesh, we consider $M=8$, and we resort to an equispaced time grid $\\{t_{\rm
g}^{j}=j\Delta t\\}_{j=0}^{J}$, with $\Delta t=5\cdot 10^{-3}$. We estimate
the long-time averaged velocity field as111In the current implementation,
$\langle u\rangle_{\rm g}$ is computed inside the time integration loop of the
Full Order Model.:

$\langle u\rangle_{\rm g}=\frac{\Delta t}{T-T_{0}}\sum_{j=J_{0}+1}^{J}u^{j},$
(5)

where $T_{0}=500$, $T=t_{\rm g}^{J}$, and $J_{0}$ is such that $t_{\rm
g}^{J_{0}}=T_{0}=500$. Consequently, we estimate the instantaneous turbulent
kinetic energy as

${\rm TKE}^{j}:=\frac{1}{2}\int_{\Omega}\,\|u^{j}-\langle u\rangle_{\rm
g}\|_{2}^{2}\,dx.$ (6)

In Appendix A, we provide a detailed analysis of the solution to the lid-
driven cavity problem (2).

In order to generate (and, later, assess) the ROM, we collect data at the
sampling times $\\{t_{\rm s}^{k}=T_{0}+\Delta t_{\rm s}k\\}_{k=1}^{K}$ with
$\Delta t_{\rm s}=1$. We observe that $\\{t_{\rm
s}^{k}\\}_{k=1}^{K}\subset\\{t_{\rm g}^{j}\\}_{j=J_{0}}^{J}$, and $K\ll J$:
this is dictated by memory constraints. We further observe that we do not
collect data in the transient region: this is motivated by the fact that we
are here ultimately interested in the long-time dynamics. In the remainder of
the paper, we use the subscript “s” to indicate the sampling times, and the
subscript “g” to indicate the time discretization. Furthermore, we use the
symbol $\langle\cdot\rangle_{\rm s}$ to indicate time averages performed based
on the sampling times, and the symbol $\langle\cdot\rangle_{\rm g}$ to
indicate time averages performed based on the time grid $\\{t_{\rm
g}^{j}\\}_{j=J_{0}}^{J}$. In Appendix B, we comment on the choice of $\Delta
t_{\rm s}$ and $K$.

##  3 The solution reproduction problem

In this section, we propose a MOR procedure for the solution reproduction
problem. As explained in the introduction, the solution reproduction problem
is of limited practical interest; however, it represents a key intermediate
step towards the development of a MOR procedure for the parametric problem.
Algorithm 1 outlines the general offline/online paradigm for the solution
reproduction problem. We recall that the offline stage is expected to be
expensive and is performed once, while the online stage should be inexpensive
and is performed many times — this distinction is of little relevance here,
but will be crucial in section 4 for the parametric problem.

Algorithm 1 Offline/online paradigm for the solution reproduction problem

Task: find an estimate of $\mathring{u}=\mathring{u}(x,t)$ of the form
$\hat{u}(x,t)=\sum_{n=1}^{N}a_{n}(t)\,\zeta_{n}(x)$.

Offline stage

1:Generate the DNS data $\\{\mathring{u}^{k}:=\mathring{u}(t_{\rm
s}^{k})\\}_{k=1}^{K}\subset V$.

2:Generate the reduced space $\mathcal{Z}^{\rm u}={\rm
span}\\{\zeta_{n}\\}_{n=1}^{N}$.

3:Formulate the Reduced Order Model.

Online stage

1:Estimate the coefficients $\\{a_{n}^{j}=a_{n}(t_{\rm g}^{j})\\}_{n=1}^{N}$
for $j=0,1,\ldots,J$.

2:Compute the QOIs (e.g., mean flow, TKE,…)

As anticipated in section 2, we here generate a ROM for the lifted velocity
field $\mathring{u}=u-R_{g}$, where $R_{g}$ is the solution to the Stokes
problem (3a). Reduction of the lifted equations is preferable from the MOR
perspective since it greatly simplifies the imposition of essential
(Dirichlet) inhomogenous boundary conditions. We observe that in the Fluid
Mechanics literature many authors consider $R_{g}=\langle u\rangle_{\rm g}$;
however, the latter choice of the lift cannot be extended to the parametric
case.

This section is organized as follows. In section 3.1, we present the POD-
Galerkin ROM. We first introduce the formulation, we review Proper Orthogonal
Decomposition for the generation of the reduced space, and then we present
numerical results that highlight the limitations of the approach. In section
3.2, we present the constrained POD-Galerkin ROM proposed in this paper. As
for POD-Galerkin, we first present and motivate the mathematical statement,
and then we present a number of numerical results to motivate the approach.

###  3.1 The POD-Galerkin ROM

####  3.1.1 The Galerkin formulation

Given the reduced space $\mathcal{Z}^{\rm u}={\rm
span}\\{\zeta_{n}\\}_{n=1}^{N}\subset V_{\rm div}=\\{v\in V:\nabla\cdot
v=0\\}$, we seek $\hat{u}\in\mathcal{V}_{N}:=H_{\rm
loc}^{1}(\mathbb{R}_{+};\mathcal{Z}^{\rm u})$ such that

$\left\\{\begin{array}[]{l}\displaystyle{\frac{d}{dt}(\hat{u}(t),v)_{L^{2}(\Omega)}+\frac{1}{\rm
Re}(\hat{u}(t)+R_{g},v)_{V}+c(\hat{u}(t)+R_{g},\hat{u}(t)+R_{g},v)=0}\\\
\hfill\forall\,v\in\mathcal{Z}^{\rm u},\\\\[8.53581pt]
\hat{u}(0)=-\Pi_{\mathcal{Z}^{\rm u}}^{L^{2}}R_{g},\\\ \end{array}\right.$ (7)

where $\Pi_{\mathcal{Z}^{\rm
u}}^{L^{2}}:[L^{2}(\Omega)]^{2}\to\mathcal{Z}^{\rm u}$ is the
$L^{2}(\Omega)$-projection operator on $\mathcal{Z}^{\rm u}$, and
$(\cdot,\cdot)_{L^{2}(\Omega)}$ is the $L^{2}(\Omega)$ inner product. If we
employ a semi-implicit time discretization, we obtain:

$\begin{array}[]{l}\displaystyle{\left(\frac{\hat{u}^{j+1}-\hat{u}^{j}}{\Delta
t},v\right)_{L^{2}(\Omega)}+\frac{1}{\rm
Re}(\hat{u}^{j+1}+R_{g},v)_{V}+c(\hat{u}^{j}+R_{g},\hat{u}^{j+1}+R_{g},v)=0}\\\
\hfill\forall\,v\in\mathcal{Z}^{\rm u},\;\;j=0,1,\ldots,\\\ \end{array}$ (8)

where $\Delta t=t_{\rm g}^{j+1}-t_{\rm g}^{j}$. We remark that the time scheme
is not the same used by nek5000. The Galerkin formulation (8) leads to the
following algebraic system for the coefficients
$\\{\mathbf{a}^{j}\\}_{j=0}^{J}$ of the $N$-term expansion:

$\mathbb{A}(\mathbf{a}^{j};{\rm
Re})\,\mathbf{a}^{j+1}=\mathbf{F}(\mathbf{a}^{j};{\rm Re}),\qquad
j=0,1,\ldots,$ (9a) where $\mathbb{A}(\mathbf{a}^{j};{\rm
Re}):=\mathbb{A}_{1}+\frac{1}{\rm
Re}\mathbb{A}_{2}+\mathbb{C}(\mathbf{a}^{j})$, $\mathbf{F}(\mathbf{a}^{j};{\rm
Re}):=\mathbb{E}\mathbf{a}^{j}-\frac{1}{\rm Re}\mathbf{G}$, with
$\begin{array}[]{l}\displaystyle{(\mathbb{A}_{1})_{m,n}=\frac{1}{\Delta
t}(\zeta_{n},\zeta_{m})_{L^{2}(\Omega)}+c(R_{g},\zeta_{n},\zeta_{m}),\quad(\mathbb{A}_{2})_{m,n}=(\zeta_{n},\zeta_{m})_{V},}\\\\[8.53581pt]
(\mathbb{C}(\mathbf{w}))_{m,n}=\sum_{i=1}^{N}\,w_{i}\,c(\zeta_{i},\zeta_{n},\zeta_{m}),\\\
\end{array}$ (9b) and
$\mathbf{G}_{m}=(R_{g},\zeta_{m})_{V},\quad(\mathbb{E})_{m,n}=\frac{1}{\Delta
t}(\zeta_{n},\zeta_{m})_{L^{2}(\Omega)}-c(\zeta_{n},R_{g},\zeta_{m}).$ (9c)

We observe that the Galerkin model for the velocity field does not contain the
pressure field. This follows from (i) the fact that the ROM is derived from
the weak form of the equations, (ii) the particular boundary conditions
prescribed, and (iii) the absence of parameters in the form $b(\cdot,\cdot)$.
We have indeed that for certain choices of the boundary conditions the ROM
should be obtained from the strong form of the Navier-Stokes equations: in
this respect, we recall that in [14] a Galerkin ROM is derived from the strong
form for a laminar flow problem with convective boundary condition ([65]) at
the outflow. In the parametric case it is possible to derive a ROM that does
not contain the pressure field if the form $b(\cdot,\cdot)$ in (4) is
parameter-independent; otherwise, it is not possible in general to generate a
space $\mathcal{Z}^{\rm u}$ such that $b(z,\cdot)\equiv 0$ for all
$z\in\mathcal{Z}^{\rm u}$ and for all values of the parameters. Since the
bilinear form $b(\cdot,\cdot)$ in (3) does not depend on the Reynolds number,
we will be able in section 4 to generate a ROM for the velocity only. We
remark that the case of parametrized $b$ form corresponds to the case of
geometric parametrizations, which is of particular interest for applications.
A potential strategy to handle this issue is to resort to the Piola’s
transform (see [66]). We refer to a future work for a detailed discussion of
this case. We also refer to the Reduced Basis literature ([67, 68, 69, 70])
for a thorough discussion about fluid problems in parametrized domains for
low-to-moderate Reynolds number flows.

The algebraic formulation (9) is the starting point for the development of the
offline/online decomposition. The matrices $\mathbb{A}_{1}$, $\mathbb{A}_{2}$,
$\mathbb{E}$, the third-order tensor $\mathbb{C}$ and the vector $\mathbf{G}$
can be pre-computed during the offline stage. Therefore, during the online
stage, the method only requires $\mathcal{O}(N^{3})$ storage, and the online
cost is $\mathcal{O}(N^{3}J)$. Provided that $N$ is much smaller than the
spatial mesh-size $\mathcal{N}$, the Galerkin ROM is significantly less
expensive and less memory-demanding than the Full Order Model. Other choices
of the time discretization lead to similar reduced systems that allow the same
offline/online decomposition.

####  3.1.2 Construction of the reduced space: Proper Orthogonal Decomposition

We employ Proper Orthogonal Decomposition (POD, [8, 9, 71]) to generate the
reduced space $\mathcal{Z}^{\rm u}$. Below we briefly review the numerical
strategy — known as _method of snapshots_ ([71]) — employed for the
computation of the POD modes. We refer to [72] for a review of the theoretical
results concerning the optimality properties of POD.

Given the snapshot set $\\{\mathring{u}^{k}\\}_{k=1}^{K}$, we assemble the
Gramian $\mathbb{U}\in\mathbb{R}^{K,K}$
$\mathbb{U}_{k,k^{\prime}}=(\mathring{u}^{k},\mathring{u}^{k^{\prime}})_{\star}$
where $(\cdot,\cdot)_{\star}$ is a suitable inner product that will be
introduced soon; then, we compute the first $N$ eigenmodes of the symmetric
matrix $\mathbb{U}$:

$\mathbb{U}\boldsymbol{\zeta}_{n}=\lambda_{n}\,\boldsymbol{\zeta}_{n},\qquad\lambda_{1}\geq\ldots\geq\lambda_{K}\geq
0;$ (10a) finally, we define the POD modes as
$\zeta_{n}:=\sum_{k=1}^{K}\,\left(\boldsymbol{\zeta}_{n}\right)_{k}\,\mathring{u}^{k},\qquad
n=1,\ldots,N.$ (10b)

It is easy to show that $\zeta_{1},\ldots,\zeta_{N}$ can be chosen to be
orthogonal in the $(\cdot,\cdot)_{\star}$-inner product; for stability
reasons, we also orthonormalize the POD modes so that
$(\zeta_{n},\zeta_{n^{\prime}})_{\star}=\delta_{n,n^{\prime}}$,
$n,n^{\prime}=1,\ldots,N$. In Appendix B, we discuss the choice of the
sampling times $\\{t_{\rm s}^{k}\\}_{k}$, and we propose a numerical technique
to assess the accuracy of the POD space for the full trajectory.

In this work, we employ the $H_{0}^{1}(\Omega)$ inner product:

$(w,v)_{\star}=(w,v)_{V}=\int_{\Omega}\,\nabla w:\nabla v\,dx.$ (11)

As explained in the introduction, this choice is motivated by dynamic
considerations. Since small-scale modes have relatively large $H^{1}$ norm
compared to their $L^{2}$ norm, and recalling that small scales are
responsible for energy dissipation, the use of the $H_{0}^{1}$ inner product
leads to a more dissipative reduced order model ([33]).

####  3.1.3 Performance of the POD-Galerkin ROM

We assess the numerical performance of the POD-Galerkin ROM presented in this
section. We here consider the lid-driven cavity problem (2) for ${\rm
Re}=15000$. We consider the time grid $\\{t_{\rm g}^{j}=\Delta
tj\\}_{j=0}^{J}$ with $\Delta t=5\cdot 10^{-3}$ and $J=2\cdot 10^{5}$
($T=t_{\rm g}^{J}=10^{3}$), and we acquire the snapshots
$\\{\mathring{u}^{k}=\mathring{u}(t_{\rm s}^{k})\\}_{k=1}^{K}$ where $t_{\rm
s}^{k}=500+k$ and $K=500$. The long-time averaged velocity field $\langle
u\rangle_{\rm g}$ is estimated through (5). On the other hand, we estimate the
mean TKE as follows:

$\langle{\rm TKE}\rangle_{\rm s}=\frac{1}{2K}\sum_{k=1}^{K}\,\|u^{k}-\langle
u\rangle_{\rm g}\|_{L^{2}(\Omega)}^{2}.$

Assembling and integration of the Reduced Order Model are performed in Matlab
[73].

Figure 1(a) shows the behavior of the eigenvalues
$\\{\lambda_{n}\\}_{n=1}^{K}$. The first eigenmode is roughly proportional to
$\langle u\rangle_{\rm g}-R_{g}$; provided that the estimate of the
coefficients is accurate, it does not contribute to the fluctuating field.
Therefore, we can identify the ratio

$r_{N}=\frac{\sum_{n=2}^{N}\lambda_{n}}{\sum_{n=2}^{K}\lambda_{n}}$

as the portion of $H_{0}^{1}$ energy of the fluctuating field associated with
the reduced POD space of dimension $N$. We find that $r_{N}=0.165$ for $N=2$,
$r_{N}=0.731$ for $N=20$, $r_{N}=0.797$ for $N=30$, and $r_{N}=0.87$ for
$N=50$. We observe that the decay with $N$ is rather slow; this suggests that
accurate estimates of the entire system dynamics are out of reach for fully
turbulent flows. Figure 1(b) shows the behavior with $N$ of the relative error
in the mean flow prediction:

$E_{N}^{0}=\frac{\|\langle u-\hat{u}\rangle_{\rm
g}\|_{L^{2}(\Omega)}}{\|\langle u\rangle_{\rm g}\|_{L^{2}(\Omega)}},\qquad
E_{N}^{1}=\frac{\|\langle u-\hat{u}\rangle_{\rm
g}\|_{H_{0}^{1}(\Omega)}}{\|\langle u\rangle_{\rm g}\|_{H_{0}^{1}(\Omega)}};$

while Figure 1(c) shows the behavior with $N$ of the mean predicted TKE:
$\langle\widehat{\rm TKE}\rangle_{\rm s}$. We observe that for small values of
$N$, we predict a false stable steady flow, while for moderate values of $N$
we substantially overestimate the TKE- Finally, for $N\gtrsim 50$ we observe a
slow convergence of the Galerkin ROM to the mean values predicted by the high-
fidelity model.

(a)

(b)

(c)

Figure 1: The solution reproduction problem; POD-Galerkin. Figure (a): POD
eigenvalues. Figure (b): behavior of the relative $L^{2}$ and $H^{1}$ errors
in mean flow prediction with $N$. Figure (c): behavior of the mean TKE with
$N$. (${\rm Re}=15000$).

Figure 2 shows the behavior for different values of $N$ of the sample mean and
sample variance of the coefficients $\\{a_{n}^{j}\\}_{j}$:

$\langle a_{n}\rangle_{\rm s}=\frac{1}{K}\sum_{k=1}^{K}\,a_{n}(t_{\rm
s}^{k}),\qquad V_{\rm
s}(a_{n})=\frac{1}{K-1}\sum_{k=1}^{K}\,\left(a_{n}(t_{\rm s}^{k})-\langle
a_{n}\rangle_{\rm s}\right)^{2},$

for the Full Order Model (FOM) and for the POD Galerkin ROM (POD-Gal). Figure
3 shows the behavior of the TKE as a function of time for three values of $N$;
predictions of first and second order moments — based on sampling times — are
reported in the caption of the Figure. Results are consistent with the results
in Figure 1. For small-to-moderate values of $N$, we observe several spurious
behaviors, namely convergence to false stable steady flows, and overly
unstable flows. As $N$ increases, the accuracy of the Galerkin ROM appears to
increase.

(a)  $N=20$

(b)  $N=40$

(c)  $N=60$

(d)  $N=20$

(e)  $N=40$

(f)  $N=60$

Figure 2: The solution reproduction problem; POD-Galerkin. Behavior of the
sample mean and sample variance of the coefficients $\\{a_{n}^{j}\\}_{j}$.
(${\rm Re}=15000$).

(a) $N=20$

(b) $N=40$

(c) $N=60$

Figure 3: The solution reproduction problem; POD-Galerkin. Behavior of the TKE
as a function of time for three values of $N$. $\langle\widehat{{\rm
TKE}}\rangle_{\rm s}=3.8\cdot 10^{-4}$ ($N=20$), $3.5\cdot 10^{-3}$ ($N=40$),
$1.1\cdot 10^{-3}$ ($N=60$). $V_{\rm s}(\widehat{{\rm TKE}})=8.8\cdot 10^{-8}$
($N=20$), $6.5\cdot 10^{-6}$ ($N=40$), $1.9\cdot 10^{-7}$ ($N=60$).
($\langle{\rm TKE}\rangle_{\rm s}=9.4\cdot 10^{-4}$, $V_{\rm s}({\rm
TKE})=8.5\cdot 10^{-8}$) (${\rm Re}=15000$).

Interestingly, the behavior of the ROM observed here is qualitatively similar
to the one observed for highly-truncated spectral approximations to turbulent
flows ([42]). We argue that the need for large reduced spaces might greatly
reduce the benefit of Model Reduction: if the value of $N$ required to obtain
sufficiently accurate results is too large, the resulting ROM might not lead
to significant computational speed-ups, and might also not be beneficial in
terms of memory. This observation motivates the correction to the Galerkin
formulation proposed in the next section. We finally remark that the results
shown in this section suggest the need for a pragmatic definition of long-time
stability: we address this issue in Appendix C.

###  3.2 The constrained POD-Galerkin formulation

####  3.2.1 Formulation

Given the reduced space $\mathcal{Z}^{\rm u}={\rm
span}\\{\zeta_{n}\\}_{n=1}^{N}\subset V_{\rm div}$, and the time grid
$\\{t_{\rm g}^{j}\\}_{j=0}^{J}$, we seek the coefficients
$\\{\mathbf{a}^{j}\\}_{j=0}^{J}\subset\mathbb{R}^{N}$ such that

$\mathbf{a}^{j+1}:={\rm
arg}\min_{\mathbf{a}\in\mathbb{R}^{N}}\,\|\mathbb{A}(\mathbf{a}^{j};{\rm
Re})\mathbf{a}-\mathbf{F}(\mathbf{a}^{j};{\rm Re})\|_{2}^{2},\quad{\rm
s.t.}\;\;\alpha_{n}\leq a_{n}\leq\beta_{n},\;\;n=1,\ldots,N;$ (12)

where $\\{\alpha_{n}\\}_{n=1}^{N}$ and $\\{\beta_{n}\\}_{n=1}^{N}$ are
suitable hyper-parameters that will be specified later, and
$\mathbb{A}(\cdot;{\rm Re})$, $\mathbf{F}(\cdot;{\rm Re})$ are defined in (9).
Formulation (12) reads as a constrained quadratic programming problem where
the objective function corresponds to the Euclidean norm of the error in the
reduced Galerkin formulation, while the constraints impose that each
coefficient of the $N$-term expansion remains in the interval
$[\alpha_{n},\beta_{n}]$, $n=1,\ldots,N$. We refer to (12) as _constrained
(POD-)Galerkin formulation_.

The hyper-parameters $\\{\alpha_{n}\\}_{n}$ and $\\{\beta_{n}\\}_{n}$ are
designed to embed in the ROM formulation information about the variation in
time of the coefficients $\\{a_{n}^{j}\\}_{j}$, for $n=1,\ldots,N$. For each
value of $n$, if we introduce the projection222We assume here that the POD
eigenmodes $\\{\zeta_{n}\\}_{n=1}^{N}$ are orthonormalized.  of the lifted
field on the $n$-th POD mode at time $t_{\rm g}^{j}$, $a_{n}^{\rm
FOM,j}:=(\mathring{u}^{j},\zeta_{n})_{V}$, we can interpret $\alpha_{n}$ and
$\beta_{n}$ as lower and upper bounds for the sequence $\\{a_{n}^{\rm
FOM,j}\\}_{j=J_{0}}^{J}$, where $J_{0}>0$ is introduced in (5) to discard the
transient dynamics. The hyper-parameters $\alpha_{n}$ and $\beta_{n}$ are not
directly related to the POD eigenvalues $\lambda_{n}$: the latter are — up to
a multiplicative constant — estimates of the squared $\ell^{2}$-norm of the
coefficients, $\lambda_{n}=\sum_{k}\left(a_{n}^{\rm
FOM,k}\right)^{2}\approx\frac{K}{J}\sum_{j}\left(a_{n}^{\rm
FOM,j}\right)^{2}$.

Based on the interpretation of the hyper-parameters, we propose to estimate
$\\{\alpha_{n}\\}_{n=1}^{N}$ and $\\{\beta_{n}\\}_{n=1}^{N}$ based on the
sample minima and the sample maxima associated with the snapshots
$\\{\mathring{u}^{k}\\}_{k=1}^{K}$:

$\alpha_{n}:=m_{n}^{\rm u}-\epsilon\Delta_{n}^{\rm
u},\qquad\beta_{n}:=M_{n}^{\rm u}+\epsilon\Delta_{n}^{\rm u},$ (13a) where
$m_{n}^{\rm u}$ and $M_{n}^{\rm u}$ are sample minimum and sample maximum
associated with the projection of the lifted field on the $n$-th POD mode,
$m_{n}^{\rm u}=\min_{k}\,a_{n}^{\rm
FOM,k}:=(\mathring{u}^{k},\zeta_{n})_{V},\quad M_{n}^{\rm
u}=\max_{k}\,a_{n}^{\rm FOM,k};$ (13b) $\Delta_{n}^{\rm u}$ is the sample
estimate of the difference between maximum and minimum, $\Delta_{n}^{\rm
u}:=M_{n}^{\rm u}-m_{n}^{\rm u};$ (13c) and the constant $\epsilon>0$ takes
into account the fact that sample minima and sample maxima in (13b) are upper
and lower bounds for the true minima and true maxima, respectively. We
emphasize that in our framework $K\ll J$; therefore, $\\{a_{n}^{\rm
FOM,k}\\}_{k=1}^{K}$ should be interpreted as a (deterministic) sample from
the population $\\{a_{n}^{\rm FOM,j}\\}_{j=J_{0}}^{J}$. Given the special
features of the learning task at hand — the estimation of minima and maxima of
a population — we expect that we can estimate the hyper-parameters based on
_sparse_ DNS data (i.e., data that are not dense in any specific region of the
time interval).

Accurate estimates of the hyper-parameters of the formulation based on sparse
DNS data represent an important feature of our constrained formulation. As
observed by many authors, low-frequency features of the turbulent flow — which
largely contribute to long-time flow averages — are well-represented by the
snapshots $\\{\mathring{u}^{k}\\}_{k=1}^{K}$ and consequently by the POD space
only if the sampling times $\\{t_{\rm s}^{k}\\}_{k=1}^{K}$ are not clustered
in any specific region of the time interval. This implies that both the
ingredients of the ROM — the space $\mathcal{Z}^{\rm u}$ and the hyper-
parameters $\\{\alpha_{n}\\}_{n}$ and $\\{\beta_{n}\\}_{n}$ — require the same
sampling strategy for the construction of the snapshot set. Therefore, the
same dataset used to generate the POD space is well-suited to estimate the
hyper-parameters of the ROM. This observation allows us to limit the size $K$
of the snapshot set, and ultimately leads to a reduction of the offline memory
cost.

Unlike the standard POD-Galerkin ROM, we here use DNS data twice: first, to
build the space $\mathcal{Z}^{\rm u}$; second, to estimate the hyper-
parameters $\\{\alpha_{n}\\}_{n=1}^{N}$ and $\\{\beta_{n}\\}_{n=1}^{N}$.
Furthermore, while POD-Galerkin is independent of the particular basis
$\zeta_{1},\ldots,\zeta_{N}$ chosen for $\mathcal{Z}^{\rm u}$, the box
constraints in (9) depend on the choice of the basis. We emphasize that by
choosing $\\{\zeta_{n}\\}_{n}$ as basis of $\mathcal{Z}^{\rm u}$ we explicitly
incorporate (prior) information about the decay of the POD coefficients
directly in the formulation.

We observe that if the solution to Galerkin ROM (9) — $\mathbf{a}_{\rm
Gal}^{j+1}=\mathbb{A}(\mathbf{a}^{j};{\rm Re})^{-1}$
$\mathbf{F}(\mathbf{a}^{j};{\rm Re})$ — satisfies the box constraints in (12),
then $\mathbf{a}^{j+1}=\mathbf{a}_{\rm Gal}^{j+1}$. Therefore, our constrained
formulation _corrects_ the Galerkin formulation only if $\mathbf{a}_{\rm
Gal}^{j+1}$ does not satisfy the prescribed bounds. This represents the main
difference between our approach and the other stabilized ROMs proposed in the
literature and briefly mentioned in the introduction: rather than introducing
artificial dissipation in the Galerkin model, we exploit prior information
about the attractor to correct the ROM.

We finally comment on time discretization. In this work, we employ the first-
order semi-implicit time-discretization introduced in (8). However, the
approach can be trivially extended to other time discretizations: first, we
derive the discrete Galerkin ROM from (7), then we substitute the resulting
algebraic formulation in the objective function of (12). For explicit and
semi-implicit single-step time integrators, the resulting constrained
formulation corresponds to a quadratic programming problem, which can be
solved using interior-point methods (see, e.g., [74]). For fully-implicit
single-step methods, the constrained formulation reads as a nonlinear
constrained optimization problem, which again can be solved using interior-
point methods or sequential quadratic programming. We envision that the
extension to multistep methods might require some additional care since the
solution is not guaranteed to be smooth in time when the constraints are
active. A thorough analysis of different time integrators is beyond the scope
of this paper.

####  3.2.2 Performance of the constrained POD-Galerkin ROM

We present numerical results for ${\rm Re}=15000$. Time grid $\\{t_{\rm
g}^{j}\\}_{j=0}^{J}$ and sampling times $\\{t_{\rm s}^{k}\\}_{k=1}^{K}$ are
the same considered for POD-Galerkin, if not specified otherwise. As for the
previous test the ROM is implemented in Matlab; the quadratic programming
problem is solved using the routine quadprog based on an interior-point
algorithm. We here set $\epsilon=0.01$ in (13).

Figure 4(a) shows the behavior of the relative $L^{2}$ and $H^{1}$ errors in
the mean flow prediction with $N$, while Figure 4(b) shows the behavior of the
mean TKE with $N$. We observe that the constrained formulation leads to a
substantial improvement in performance compared to the standard POD-Galerkin
method (cf. Figures 1(b) and 1(c)): for $N\gtrsim 40$ the relative error in
the mean is less than $2\%$, while the predicted mean TKE is bounded from
above by $\langle{\rm TKE}\rangle_{\rm s}$ for all values of $N$. Furthermore,
we observe that the TKE of our constrained Galerkin formulation is larger than
the one predicted by the Galerkin ROM for certain values of $N$, and is
smaller for other values of $N$: this empirically proves that our approach
does not necessarily add dissipation to the Galerkin ROM. In Figure 5, we
repeat the tests of Figure 4 for $\Delta t^{\prime}=0.5\Delta t=2.5\cdot
10^{-3}$. We observe that results are consistent with the results shown in
Figure 4: this provides empirical evidence for the stability of our
constrained formulation under time-step refinement. Figure 6 shows the
behavior of the sample mean and sample variance of the coefficients
$\\{a_{n}^{j}\\}_{j}$ for three different values of $N$. Also in this case, we
observe a substantial improvement in performance compared to POD-Galerkin,
particularly for high modes. Finally, Figure 7 shows the behavior of the TKE
as a function of time for three values of $N$. We observe that for $N=40$ and
$N=60$ the predicted TKE is in good qualitative agreement with the truth; in
addition, predictions of first- and second-order moments (reported in the
caption) are accurate. In Appendix D, we present additional results to
demonstrate the efficiency of the constrained formulation, and also the
robustness with respect to the choice of $\epsilon$.

(a)

(b)

Figure 4: The solution reproduction problem; constrained POD-Galerkin. Figure
(a): behavior of the relative $L^{2}$ and $H^{1}$ errors in mean flow
prediction with $N$. Figure (b): behavior of the mean TKE with $N$. (${\rm
Re}=15000$, $\epsilon=0.01$).

(a)

(b)

Figure 5: The solution reproduction problem; constrained POD-Galerkin for a
finer time grid. Figure (a): behavior of the relative $L^{2}$ and $H^{1}$
errors in mean flow prediction with $N$. Figure (b): behavior of the mean TKE
with $N$. (${\rm Re}=15000$, $\epsilon=0.01$, $\Delta t=2.5\cdot 10^{-3}$).

(a)  $N=20$

(b)  $N=40$

(c)  $N=60$

(d)  $N=20$

(e)  $N=40$

(f)  $N=60$

Figure 6: The solution reproduction problem; constrained POD-Galerkin.
Behavior of the sample mean and sample variance of the coefficients
$\\{a_{n}^{j}\\}_{j}$. (${\rm Re}=15000$, $\epsilon=0.01$).

(a) $N=20$

(b) $N=40$

(c) $N=60$

Figure 7: The solution reproduction problem; constrained POD-Galerkin.
Behavior of the TKE as a function of time for three values of $N$.
$\langle\widehat{{\rm TKE}}\rangle_{\rm s}=8.6\cdot 10^{-4}$ ($N=20$),
$9.4\cdot 10^{-4}$ ($N=40$), $7.7\cdot 10^{-4}$ ($N=60$). $V_{\rm
s}(\widehat{{\rm TKE}})=5.5\cdot 10^{-9}$ ($N=20$), $1.7\cdot 10^{-7}$
($N=40$), $5.8\cdot 10^{-8}$ ($N=60$). ($\langle{\rm TKE}\rangle_{\rm
s}=9.4\cdot 10^{-4}$, $V_{\rm s}({\rm TKE})=8.5\cdot 10^{-8}$) (${\rm
Re}=15000$, $\epsilon=0.01$).

##  4 The parametric problem

We consider the extension of our MOR approach to the parametric context. For
the purpose of exposition, we focus our discussion on the lid-driven cavity
problem presented in section 2: we wish to estimate the solution to (2) for
${\rm Re}\in\mathcal{P}=[15000,25000]$. In view of the $h$-refinement, we
introduce the partition of $\mathcal{P}$,
$\mathcal{I}_{1},\ldots,\mathcal{I}_{M}$ such that
$\bigcup_{m=1}^{M}\mathcal{I}_{m}=\mathcal{P}$,
$\mathcal{I}_{m}\cap\mathcal{I}_{m^{\prime}}=\emptyset$. We seek an estimate
of the lifted velocity field $\mathring{u}=\mathring{u}(x,t;{\rm Re})$ of the
form $\hat{u}(x,t;{\rm Re})=\sum_{n=1}^{N}a_{n}^{m}(t;{\rm
Re})\zeta_{n}^{m}(x)$ for all ${\rm Re}\in\mathcal{I}_{m}$. The approach can
be trivially extended to other parametric problems that do not involve
geometric variations; as already mentioned in section 3, the extension to the
latter case is beyond the scope of the present work. Algorithm 2 summarizes
the general offline/online paradigm for the parametric problem. We highlight
that, for the sake of generality, in Algorithm 2 we distinguish between $L$
(number of offline solves) and $M$ (number of reduced spaces). However, in
this work we consider the case $L=M$.

In order to tackle the parametric problem outlined above, we should address
two challenges: first, we should extend the constrained formulation to the
parametric case; second, we should develop a Greedy strategy for the proper
selection of the parameters ${\rm Re}_{1}^{\star},\ldots,{\rm
Re}_{L}^{\star}$, and the partition $\\{\mathcal{I}_{m}\\}_{m}$. We emphasize
that a proper selection of the parameters reduces the number of offline full
order solves, and is thus crucial for the feasibility of the approach. In
order to address the first challenge, we propose an actionable procedure for
the selection of the hyper-parameters $\\{\alpha_{n}\\}_{n}$ and
$\\{\beta_{n}\\}_{n}$ associated with the constrained formulation (12) in the
parametric case. On the other hand, the Greedy approach relies on an
inexpensive error indicator, which corresponds to the dual norm of the
residual associated with the time-averaged momentum equation.

Algorithm 2 Offline/online paradigm for the parametric problem

Task: find an estimate of the lifted velocity field
$\mathring{u}=\mathring{u}(x,t;{\rm Re})$ of the form $\hat{u}(x,t;{\rm
Re})=\sum_{n=1}^{N}a_{n}^{m}(t;{\rm Re})\zeta_{n}^{m}(x)$ for all ${\rm
Re}\in\mathcal{I}_{m}$, $m=1,\ldots,M$.

Offline stage

1:Generate the DNS data $\\{\mathring{u}^{k}({\rm
Re_{\ell}}^{\star}):=\mathring{u}(t_{\rm s}^{k};{\rm
Re_{\ell}^{\star}})\\}_{k=1}^{K}\subset V$, and ${\rm
Re}_{1}^{\star},\ldots,{\rm Re}_{L}^{\star}\in\mathcal{P}$.

2:Generate the partition $\\{\mathcal{I}_{m}\\}_{m}$ of $\mathcal{P}$, and the
reduced spaces $\mathcal{Z}_{m}^{\rm u}={\rm
span}\\{\zeta_{n}^{m}\\}_{n=1}^{N}$, $m=1,\ldots,M$.

3:Formulate the Reduced Order Models for each subregion.

Online stage

1:Given ${\rm Re}\in\mathcal{P}$, find $m\in\\{1,\ldots,M\\}$ such that ${\rm
Re}\in\mathcal{I}_{m}$.

2:Estimate the coefficients $\\{a_{n}^{m,j}({\rm Re})=a_{n}^{m}(t_{\rm
g}^{j};{\rm Re})\\}_{n=1}^{N}$ for $j=0,1,\ldots,J$.

3:Compute the QOIs (e.g., mean flow, TKE,…).

The section is organized as follows. In section 4.1 we present the
POD-$h$Greedy approach, in section 4.2 we present the ROM formulation, and we
discuss the choice of the hyper-parameters $\\{\alpha_{n}\\}_{n}$ and
$\\{\beta_{n}\\}_{n}$. Then, in section 4.3 we propose the time-averaged error
indicator. Finally, in section 4.4, we present the numerical results for the
lid-driven cavity problem.

###  4.1 POD-Greedy algorithm

We first present the POD-$h$Greedy algorithm for the construction of the
reduced spaces $\\{\mathcal{Z}_{m}^{\rm u}\\}_{m=1}^{M}$, and the partition
$\\{\mathcal{I}_{m}\\}_{m=1}^{M}$ of $\mathcal{P}$, based on the results of
$L$ FOM simulations associated with the parameters ${\rm
Re}_{1}^{\star},\ldots,{\rm Re}_{L}^{\star}$. The approach is a simplified
version of the $h$-refinement procedure for parabolic problems proposed in
[45]. In particular, we here consider $L=M$: this implies that each reduced
space $\mathcal{Z}_{\ell}^{\rm u}$ is based on the POD of a single full-order
solve.

In view of the presentation of the algorithm, we introduce the discretized
parameter space $\mathcal{P}_{\rm train}=\\{{\rm Re}_{i}\\}_{i=1}^{n_{\rm
train}}$, ${\rm Re}_{1}\leq\ldots\leq{\rm Re}_{n_{\rm train}}$, the integers
$L$ and $N$, which fix the maximum number of offline solves and the size of
the reduced space $\mathcal{Z}^{\rm u}$, the integer $n_{\rm cand}<L$, which
is the number of ROM evaluations performed online for a given value of the
parameters, and the _a posteriori_ error indicator $\Delta^{\rm
u}:\bigotimes_{j=0}^{J}V\times\mathcal{P}\to\mathbb{R}_{+}$. The error
indicator takes as input a sequence $\\{w^{j}\\}_{j=0}^{J}\subset V$ and the
value of the parameter, and returns an estimate of the error in the prediction
of the mean flow; we formally present the indicator in section 4.3. We further
introduce the functions

$\begin{array}[]{l}\displaystyle{[\\{\zeta_{n}\\}_{n=1}^{N}]=\texttt{POD}_{V}\left(\mathcal{S},N\right);\quad[\\{\mathring{u}^{k}({\rm
Re})\\}_{k=1}^{K}]=\texttt{DNS-solver}\left({\rm Re},\\{t_{\rm
s}^{k}\\}_{k=1}^{K}\right);}\\\\[11.38109pt]
\displaystyle{[\\{\hat{u}^{j}({\rm Re})\\}_{j=0}^{J}]=\texttt{ROM-
solver}\left({\rm Re},\mathcal{Z}^{\rm u}\right).}\\\ \end{array}$

$\texttt{POD}_{V}$ takes as input the set of snapshots
$\mathcal{S}=\\{w^{i}\\}_{i=1}^{|\mathcal{S}|}$ and an integer $N>0$, and
returns the orthonormalized first $N$ POD eigenmodes (see section 3.1.2); on
the other hand, DNS-solver takes as input the value of the Reynolds number and
the sampling times $\\{t_{\rm s}^{k}\\}_{k=1}^{K}$, and returns the
instantaneous velocity at times $\\{t_{\rm s}^{k}\\}_{k=1}^{K}$; finally, ROM-
solver takes as input the value of the Reynolds number and the reduced space
$\mathcal{Z}^{\rm u}$, and returns the ROM solution
$\hat{u}^{j}=\sum_{n=1}^{N}\,a_{n}^{j}\zeta_{n}$ for each time step of the
grid $\\{t_{\rm g}^{j}\\}_{j=0}^{J}$. Algorithm 3 presents the computational
procedure for both offline and online stage. With some abuse of notation, we
use $\Delta_{\ell}^{\rm u}(\cdot)$ to refer to the error estimate associated
with the $\ell$-th model.

Algorithm 3 POD-$h$Greedy algorithm for the construction of
$\\{\mathcal{Z}_{\ell}^{\rm u},\mathcal{I}_{\ell}\\}_{\ell}$

Offline stage: $[\\{\mathcal{Z}_{\ell}^{\rm u}\\}_{\ell=1}^{L}]=$ Offline
$(\mathcal{P}_{\rm train},N,L,\Delta^{\rm u},\\{t_{\rm s}^{k}\\}_{k=1}^{K})$.

_Inputs:_ $\mathcal{P}_{\rm train}=\\{{\rm Re}_{i}\\}_{i=1}^{n_{\rm train}}=$
discretized parameter space, $N=$ dimension of each reduced space, $L=$
maximum number of offline solves, $\Delta^{\rm u}=$ error indicator,
$\\{t_{\rm s}^{k}\\}_{k=1}^{K}=$ sampling times.

_Output:_ $\\{(\mathcal{Z}_{\ell}^{\rm u},{\rm
Re}_{\ell}^{\star}\\}_{\ell=1}^{L}=$ reduced space/anchor point pairs.

1:${\rm Re}_{1}^{\star}={\rm rand}(\mathcal{P}_{\rm train})$

2:for $\ell=1,\ldots,L$ do

3: $[\\{\mathring{u}^{k}({\rm Re}_{\ell}^{\star})\\}_{k=1}^{K}]=\texttt{DNS-
solver}\left({\rm Re}_{\ell}^{\star},\\{t_{\rm s}^{k}\\}_{k=1}^{K}\right)$

4:
$[\\{\zeta_{n}^{\ell}\\}_{n=1}^{N}]=\texttt{POD}_{V}\left(\\{\mathring{u}^{k}({\rm
Re}_{\ell}^{\star})\\}_{k=1}^{K},N\right)$

5: Define $\mathcal{Z}_{\ell}^{\rm u}={\rm
span}\\{\zeta_{n}^{\ell}\\}_{n=1}^{N}$, build the ROM structures (cf. sections
3.1.1 and 4.3).

6: for $i=1,\ldots,n_{\rm train}$ do

7: $[\\{\hat{u}_{\ell}^{j}({\rm Re}_{i})\\}_{j=0}^{J}]=\texttt{ROM-
solver}\left({\rm Re}_{i},\mathcal{Z}_{\ell}^{\rm u}\right)$;

8: Compute the error estimate $\Delta_{\ell}^{\rm u}({\rm Re}_{i})$

9: end for

10: ${\rm Re}_{\ell+1}^{\star}={\rm arg}\max_{{\rm Re}\in\mathcal{P}_{\rm
train}}\,\min_{\ell^{\prime}=1,\ldots,\ell}\Delta_{\ell^{\prime}}^{\rm u}({\rm
Re})$.

11:end for

Online stage: $[\\{\hat{u}^{j}\\}_{j}]=$ Online $(\\{(\mathcal{Z}_{\ell}^{\rm
u},{\rm Re}_{\ell}^{\star}\\}_{\ell=1}^{L},\Delta^{\rm u},n_{\rm cand},{\rm
Re})$.

_Inputs:_ $\\{(\mathcal{Z}_{\ell}^{\rm u},{\rm
Re}_{\ell}^{\star}\\}_{\ell=1}^{L}=$ reduced space/anchor point pairs, $n_{\rm
cand}=$ online ROM evaluations, ${\rm Re}=$ input parameter.

_Output:_ $\\{\hat{u}^{j}\\}_{j}=$ solution estimate.

1:Find the $n_{\rm cand}$ nearest anchors to ${\rm Re}$: ${\rm
Re}_{(1)}^{\star},\ldots,{\rm Re}_{(n_{\rm cand})}^{\star}$

2:for $i=1,\ldots,n_{\rm cand}$ do

3: $[\\{\hat{u}_{(i)}^{j}\\}_{j=0}^{J}]=\texttt{ROM-solver}\left({\rm
Re},\mathcal{Z}_{(i)}^{\rm u}\right)$

4: Compute the error estimate $\Delta_{(i)}^{\rm u}({\rm Re})$

5:end for

6:Return $\\{\hat{u}^{j}=\hat{u}_{(i^{\star})}^{j}\\}_{j}$, where $i^{\star}$
is the minimizer of $\\{\Delta_{(i)}^{\rm u}({\rm Re})\\}_{i}$.

Algorithm 3 combines a POD in time with a Greedy in parameter. As explained in
the introduction, Greedy techniques are crucial to allow efficient parameter
explorations at an affordable offline computational cost. We emphasize that
our approach is different from the POD-Greedy strategy proposed in [44]:
rather than building $L$ different reduced spaces, the authors of [44] combine
data from different parameters to generate a single reduced space.

We also remark that our definition of the partition — see Algorithm 3, Online
stage — can be formally expressed as follows:

$\mathcal{I}_{\ell}=\left\\{{\rm Re}\in\mathcal{P}:\,\ell\in I({\rm
Re}),\;\Delta_{\ell}^{\rm u}({\rm Re})<\Delta_{\ell^{\prime}}^{\rm u}({\rm
Re}),\ell^{\prime}\in I({\rm Re}),\ell^{\prime}\neq\ell\right\\},$ (14)

where $I({\rm Re})\subset\\{1,\ldots,L\\}$ is the set of indices associated
with the $n_{\rm cand}$ nearest anchor points to ${\rm Re}$.

###  4.2 Constrained Galerkin formulation

Given the reduced spaces $\\{\mathcal{Z}_{\ell}^{\rm u}={\rm
span}\\{\zeta_{n}^{\ell}\\}_{n=1}^{N}\\}_{\ell=1}^{L}$, we consider the
constrained Galerkin formulation proposed in section 3: given ${\rm
Re}\in\mathcal{P}$, and the time grid $\\{t_{\rm g}^{j}\\}_{j=0}^{J}$, the
$\ell$-th ROM seeks the coefficients
$\\{\mathbf{a}_{\ell}^{j}\\}_{j=0}^{J}\subset\mathbb{R}^{N}$ such that

$\begin{array}[]{l}\displaystyle{\mathbf{a}_{\ell}^{j+1}:={\rm
arg}\min_{\mathbf{a}\in\mathbb{R}^{N}}\,\|\mathbb{A}_{\ell}(\mathbf{a}^{j};{\rm
Re})\mathbf{a}-\mathbf{F}_{\ell}(\mathbf{a}^{j};{\rm Re})\|_{2}^{2},\;\;\;{\rm
s.t.}\;\;\alpha_{n}^{\ell}({\rm Re})\leq a_{n}\leq\beta_{n}^{\ell}({\rm
Re}),}\\\ \hfill n=1,\ldots,N;\\\ \end{array}$ (15)

where $\mathbb{A}_{\ell}$ and $\mathbf{F}_{\ell}$ can be computed by
exploiting (9) for the reduced space $\mathcal{Z}_{\ell}^{\rm u}$, and the
constraints $\\{\alpha_{n}^{\ell}\\}_{n,\ell}$ and
$\\{\beta_{n}^{\ell}\\}_{n,\ell}$ are based on the DNS data for the anchor
point ${\rm Re}_{\ell}^{\star}$. In greater detail, given
$\ell\in\\{1,\ldots,L\\}$, assuming that
$\zeta_{1}^{\ell},\ldots,\zeta_{N}^{\ell}$ are orthonormal in $V$, we define
$\\{\alpha_{n}^{\ell}\\}_{n}$ and $\\{\beta_{n}^{\ell}\\}_{n}$ such that

$\alpha_{n}^{\ell}:=m_{n,\ell}^{\rm u}-\epsilon\Delta_{n,\ell}^{\rm
u},\quad\beta_{n}^{\ell}:=M_{n,\ell}^{\rm u}+\epsilon\Delta_{n,\ell}^{\rm u},$
(16a) where $m_{n,\ell}^{\rm u}:=\min_{k}\,a_{n,\ell}^{\rm FOM,k},\quad
M_{n,\ell}^{\rm u}:=\max_{k}\,a_{n,\ell}^{\rm FOM,k},\quad\Delta_{n,\ell}^{\rm
u}:=M_{n,\ell}^{\rm u}-m_{n,\ell}^{\rm u},$ (16b) and $a_{n,\ell}^{\rm
FOM,k}:=(\mathring{u}^{k}({\rm Re}_{\ell}^{\star}),\zeta_{n}^{\ell})_{V}$.

The offline/online decomposition is equivalent to the one described in section
3.1.1. We omit the details.

We observe that our choices of $\alpha_{n}^{\ell}$ and $\beta_{n}^{\ell}$
correspond to a constant approximation of the functions

$m_{n,\ell}^{\rm FOM,u}({\rm
Re}):=\min_{j=J_{0},\ldots,J}\,\left(\mathring{u}^{j}({\rm
Re}),\zeta_{n}^{\ell}\right)_{V},\;\;M_{n,\ell}^{\rm FOM,u}({\rm
Re}):=\max_{j=J_{0},\ldots,J}\,\left(\mathring{u}^{j}({\rm
Re}),\zeta_{n}^{\ell}\right)_{V};$ (17)

where $\mathcal{I}_{\ell}\subset\mathcal{P}$ is defined in (14). We observe
that the piece-wise constant approximations of $m_{n,\ell}^{\rm FOM,u}$ and
$M_{n,\ell}^{\rm FOM,u}$ are justified by our Greedy algorithm, which
adaptively determines the partition of $\mathcal{P}$ based on the error
indicator. For practical parametrizations, and practical values of $L$ (i.e.,
number of offline solves) we expect that accurate estimates of
$m_{n,\ell}^{\rm FOM,u}$ and $M_{n,\ell}^{\rm FOM,u}$ over $\mathcal{P}$ might
be out of reach. Therefore, we here effectively rely on (i) the robustness of
our constrained approach to perturbations in the value of the hyper-
parameters, and (ii) the weak sensitivity of the functions $m_{n,\ell}^{\rm
ROM,u}$ and $M_{n,\ell}^{\rm ROM,u}$ with respect to the parameter. For the
lid-driven cavity problem considered in this work, we provide numerical
evidence to support these two assumptions in Appendix D.

###  4.3 A time-averaged error indicator

Given the sequence $\\{w^{j}\\}_{j=0}^{J}\subset V$ and ${\rm
Re}\in\mathcal{P}$, we define the discrete time-averaged residual $\langle
R\rangle:\bigotimes_{j=0}^{J}V\times V_{\rm
div}\times\mathcal{P}\to\mathbb{R}$ associated with (8):

$\langle R\rangle\left(\\{w^{j}\\}_{j=0}^{J},v;\,{\rm Re}\right)=\frac{\Delta
t}{T-T_{0}}\sum_{j=J_{0}}^{J-1}\,e(\hat{u}^{j},\hat{u}^{j+1},{\rm Re})$ (18a)
where $T=t_{\rm g}^{J}$, $T_{0}=t_{\rm g}^{J_{0}}$, and
$\begin{array}[]{l}\displaystyle{e(\hat{u}^{j},\hat{u}^{j+1},{\rm
Re}):=\left(\frac{\hat{u}^{j+1}-\hat{u}^{j}}{\Delta
t},v\right)_{L^{2}(\Omega)}\,+\frac{1}{\rm
Re}(\hat{u}^{j+1}+R_{g},v)_{V}}\\\\[8.53581pt]
\hfill+c(\hat{u}^{j}+R_{g},\hat{u}^{j+1}+R_{g},v),\,j=J_{0},\ldots,J-1\\\
\end{array}$ (18b)

Then, we define the error indicator $\Delta^{\rm
u}:\bigotimes_{j=0}^{J}V\times\mathcal{P}\to\mathbb{R}_{+}$ as follows:

$\Delta^{\rm u}\left(\\{w^{j}\\}_{j=0}^{J};\,{\rm Re}\right):=\big{\|}\langle
R\rangle\left(\\{w^{j}\\}_{j=0}^{J},\cdot;\,{\rm Re}\right)\big{\|}_{V_{\rm
div}^{\prime}}$ (19)

where $\|\cdot\|_{V_{\rm div}^{\prime}}$ denotes the norm of the dual space
$V_{\rm div}^{\prime}$.

In our numerical tests, as in (5), we consider $J_{0}$ such that $t_{\rm
g}^{J_{0}}=T_{0}=500$: this choice is designed to limit the effect of the
transient dynamics. It is easy to verify that the solution to the FOM for any
initial condition — provided that the same time discretization is employed —
satisfies $\Delta^{\rm u}\equiv 0$. This implies that two sequences
$\\{w^{j}\\}_{j=0}^{J},\\{\tilde{w}^{j}\\}_{j=0}^{J}\subset V$ satisfying
$\Delta^{\rm u}\equiv 0$ might be far from each other at each time step (i.e.,
$\|\tilde{w}^{j}-w^{j}\|_{V}$ is large for any $j\geq 0$). However, for
sufficiently large values of $J$, we expect $\Delta^{\rm u}$ to be highly-
correlated with the error in the mean flow prediction; for this reason, we can
exploit $\Delta^{\rm u}$ to guide the Greedy algorithm presented in section
4.1. We empirically investigate the correlation between $\Delta^{\rm u}$ and
the error in the mean flow prediction in the numerical experiments at the end
of the section. A theoretical justification of the error indicator is beyond
the scope of the present work.

The error indicator $\Delta^{\rm u}$ can be computed efficiently for sequences
in $\mathcal{Z}^{\rm u}$ exploiting an offline/online computational
decomposition; the procedure is standard in the Reduced Basis literature, and
is reported in Appendix F.

###### Remark 4.1

We do not expect that the residual indicator (19) is in good quantitative
agreement with the error in mean flow prediction $\|\langle
u-\hat{u}\rangle_{\rm g}\|_{V}$. More precisely, if we define the effectivity
$\eta:=\frac{\Delta^{\rm u}}{\|\langle u-\hat{u}\rangle_{\rm g}\|_{V}}$ of the
residual error indicator, we do not expect that $\eta$ is close to one.

In order to obtain a quantitative estimate of the error of the ROM anchored in
${\rm Re}^{\star}$, we can consider the corrected estimator

$\Delta^{\rm u,corr}({\rm Re};{\rm Re}^{\star}):=\frac{1}{\eta({\rm
Re}^{\star})}\,\Delta^{\rm u}({\rm Re}),$ (20)

where $\Delta^{\rm u}({\rm Re})$ is the error indicator associated with the
ROM anchored in ${\rm Re}^{\star}$, and $\eta({\rm Re}^{\star})$ is the
effectivity evaluated at $\rm Re={\rm Re}^{\star}$. Note that the computation
of $\eta({\rm Re}^{\star})$ does not require any additional call to the DNS
solver.

###  4.4 Numerical results

Figure 8 shows the results of the application of Algorithm 3 for the
construction of the ROM for the parametric problem. In order to assess
performance, we generate DNS data for $t_{\rm g}^{j}\in\\{0,\ldots,1500\\}$,
$\\{t_{\rm s}^{k}=500+k\\}_{k=1}^{K=1000}$ for ${\rm
Re}=15000,16000,\ldots,25000$ ($n_{\rm train}=11$ datapoints). Then, we apply
Algorithm 3 with ${\rm Re}_{1}^{\star}=15000$, $N=80$, and $\epsilon=0.05$. We
perform $L=3$ iterations of the Greedy procedure. Figure 8(a) shows the
behavior of $\Delta^{\rm u}$ with ${\rm Re}$ for the three iterations, while
Figure 8(b) shows the behavior of the relative $H^{1}$ error in mean flow
prediction with ${\rm Re}$. The black continuous line denotes the performance
of the reduced model which minimizes the error indicator, and thus is selected
by the Greedy procedure (cf. Algorithm 3, $n_{\rm cand}=2$). We observe that
the maximum relative error decreases at each iteration, and it is roughly
$13\%$ after the third iteration.

(a)

(b)

Figure 8: The parametric problem; performance of POD-$h$Greedy. Figure (a):
behavior of $\Delta^{\rm u}$ with ${\rm Re}$ for three iterations. Figure (b):
behavior of the relative $H^{1}$ error in mean flow prediction with ${\rm Re}$
for three iterations; the black line (est) shows the performance of the
reduced model which minimizes the error indicator (and thus is selected by the
Greedy procedure). ($\epsilon=0.05$, $N=80$,${\rm Re_{1}}=15000,{\rm
Re_{2}}=25000,{\rm Re_{3}}=17000$).

Results of Figure 8 show the importance of the error indicator $\Delta^{\rm
u}$ in (19) to select the parameters ${\rm Re}_{2}$ and ${\rm Re}_{3}$, and
also motivate the choice of the partition $\\{\mathcal{I}_{\ell}\\}_{\ell}$ in
(14): after the third iteration, for $10$ out of $11$ values of the Reynolds
number, the reduced model that minimizes the error indicator (over all models)
is the same that minimizes the true error. On the other hand, we observe that
the indicator is in poor quantitative agreement with the true error: Figure
9(a) shows that the effectivity $\eta$ of the error indicator is
$\mathcal{O}(10^{-3})$ for all three reduced order models and for all values
of the Reynolds numbers considered. However, Figure 9(b) shows that the
correction proposed in Remark 4.1 leads to an indicator that is in reasonable
quantitative agreement with the error in mean flow prediction.

(a)

(b)

Figure 9: The parametric problem; error estimator. Figure (a): behavior of the
effectivity $\eta$ of the error indicator $\Delta^{\rm u}$ (19) for the three
ROMs. Figure (b): behavior of the effectivity $\eta^{\rm corr}$ of the
corrected error indicator $\Delta^{\rm u,corr}$ (20) for the three ROMs.

Figure 10 shows the behavior of the TKE with time for three values of the
Reynolds number, ${\rm Re=}16000,20000,23000$, which have not been selected by
the Greedy procedure. Here, predictions are based on the ROM after three
iterations of Algorithm 2.

(a)  ${\rm Re}=16000$

(b) ${\rm Re}=20000$

(c) ${\rm Re}=23000$

Figure 10: The parametric problem; behavior of the TKE with time for three
values of the Reynolds number. ($\epsilon=0.05$, $N=80$, ${\rm
Re_{1}}=15000,{\rm Re_{2}}=25000,{\rm Re_{3}}=17000$). Predictions for ${\rm
Re}_{1}$ and ${\rm Re}_{2}$ rely on the ROM anchored in ${\rm
Re}^{\star}=17000$, while predictions for ${\rm Re}_{3}$ rely on the ROM
anchored in ${\rm Re}^{\star}=25000$.

##  5 Conclusions

In this paper, we present a Reduced Basis technique for long-time integration
of turbulent flows. The three contributions of this work are (i) a constrained
Galerkin formulation that correct the Galerkin statement by incorporating
prior information about the long-time attactor, (ii) an inexpensive time-
averaged indicator for the error in mean flow prediction, and (iii) a
POD-$h$Greedy technique for the construction of the ROM. In order to assess
performance, we apply our approach to a lid-driven cavity problem parametrized
with respect to the Reynolds number: first, we consider the solution
reproduction problem (non-predictive case) to demonstrate the effectivity of
our new constrained formulation; second, we consider the parametric problem
(predictive case) to validate our error indicator, and more broadly the POD-
Greedy procedure.

Our constrained Galerkin formulation is able to accurately predict mean flow
and also the TKE. The error indicator, despite it is not corroborated by a
firm theoretical analysis, is found to be highly-correlated with the error in
the prediction of the mean flow; hence, it is naturally suited to drive the
offline Greedy.

In this paper, we also highlight a number of challenges, which are
particularly relevant for turbulent flows, and that should be taken into
consideration in the design of MOR strategies for turbulent flows: first, the
slow convergence of the Kolmogorov $N$-width suggested by Figure 1(a) which
prevents us from accurately representing the full dynamics; second, the
difficulty to combine modes associated with different parameters (cf. Appendix
E); third, the large offline costs both in terms of computational time and
required storage. In this paper, we propose to address the first challenge by
_reducing our goal_ : rather than trying to estimate the full trajectory, we
develop a ROM uniquely for the prediction of first and second moments of the
long-time dynamics. Furthermore, we propose to address the second challenge by
resorting to an $h$-refinement in parameter. On the other hand, we here
postulate that the snapshot set $\\{u^{k}\\}_{k=1}^{K}$ is rich enough to
accurately estimate the first $N$ POD modes associated with the full
trajectory $\\{\mathring{u}^{j}\\}_{j=J_{0}}^{J}$, and also that it is
possible to compute and store the Riesz representers
$\tilde{\xi}_{1},\ldots,\tilde{\xi}_{M}$, $M=N^{2}+3N+2$, for residual
calculations. In Appendix B, we review a computational strategy to assess _a
posteriori_ the representativity of our snapshot set; on the other hand, we
refer to a future work for the development of computational strategies to
reduce the offline costs related to residual calculations.

We finally outline a number of potential next steps that we wish to pursue in
the future.

  * •

_Constrained formulation_ Our constrained formulation minimizes the $\ell^{2}$
error in the reduced Galerkin statement subject to lower and upper bounds for
the coefficients of the $N$-term expansion. We wish to consider other choices
both for the objective functions and for the constraints. In particular, we
wish to minimize the residual in a suitable dual norm, and we wish to design
other constraints to take into account the properties of the attractor.
Furthermore, we also wish to consider the post-processing rectification method
proposed in [75] to improve the accuracy of the mean flow. Finally, we wish to
consider alternative strategies for writing the nonlinear term in the momentum
equation, and also for imposing strong boundary conditions.

  * •

_hp-Greedy_ In Appendix E, we discuss the limitation of the traditional
POD-($p$)Greedy algorithm. However, in order to tackle complex
parametrizations, we envision that the $h$-refinement strategy proposed in
this paper might require an unfeasible number of offline simulations. This is
why we wish to consider more advanced sampling strategies that combine
$h$-refinement and $p$-refinement.

  * •

_Extension to more challenging problems_ We wish to consider geometry
variations, which are particularly relevant for applications. As explained in
the body of the paper, this might be accomplished by resorting to the Piola
transform, or by considering a two-field (velocity and pressure) formulation.
Furthermore, we wish to apply our approach to transient problems: in order to
face this task, we envision that time-dependent constraints should be
considered, and also the time-averaged error indicator should be modified
based on the particular quantity of interest we wish to predict. Finally, we
wish to apply our approach to the reduction of LES/URANS flow simulations.
This would substantially increase the range of engineering applications we can
tackle with our method.

##  Appendix A Analysis of the solution to the lid-driven cavity problem

Figure 11 shows the velocity streamlines for three different times for ${\rm
Re}=15000$, while Figure 12 shows the velocity streamlines for several times
for ${\rm Re}=20000$: for the latter value of the Reynolds number, we observe
the presence of vortices along the edges and in the center of the cavity. We
remark that the same behavior has been observed by Cazemier et al. in [54]
(cf. Figure 3, page 1687). This rare behavior makes estimates of long-time
averages particularly difficult, especially for the TKE. Figure 13 shows the
behavior of the turbulent kinetic energy ${\rm TKE}$ with time for three
values of the Reynolds number. We observe that for sufficiently large values
of the Reynolds number we have significant peaks in the TKE. These peaks
correspond to eddies that are ejected into the core region and cross the
cavity.

(a) $t=501$

(b) $t=600$

(c) $t=700$

Figure 11: A lid-driven cavity problem. Velocity streamlines for ${\rm
Re}=15000$.

(a) $t=1252$

(b) $t=1266$

(c) $t=1276$

(d) $t=1286$

(e) $t=1320$

(f) $t=1344$

Figure 12: A lid-driven cavity problem. Velocity streamlines for ${\rm
Re}=20000$ for several time steps.

(a) ${\rm Re}=15000$

(b) ${\rm Re}=20000$

(c) ${\rm Re}=25000$

Figure 13: A lid-driven cavity problem. Behavior of the turbulent kinetic
energy ${\rm TKE}$ with time for three values of the Reynolds number.

Figure 14 shows the behavior with time of the first and second components of
the velocity field at three spatial locations for ${\rm Re}=15000$ and ${\rm
Re}=25000$. We observe that for $t\gtrsim T_{0}=500$ the effects of the
transient dynamics are negligible. Figure 15 shows the autocorrelation factors
associated with the time series $\\{u_{i}(x_{\ell}^{\rm probe},t^{j},{\rm
Re})\\}_{j=J_{0}}^{J}$ for $i,\ell=1,2$ and for ${\rm Re}=15000,25000$. We
here define the autocorrelation factors for a time sequence
$\\{y^{j}\\}_{j=0}^{J}$ as follows:

$\rho_{\rm g}(\tau=\kappa\Delta t)=\frac{\displaystyle{\frac{1}{J-\kappa-
J_{0}+1}\sum_{j=J_{0}}^{J-\kappa}\left(y^{j}-\langle y\rangle_{\rm
g}\right)\left(y^{j+\kappa}-\langle y\rangle_{\rm
g}\right)}}{\displaystyle{\langle\left(y-\langle y\rangle_{\rm
g}\right)^{2}\rangle_{\rm g}}}.$

We observe that the autocorrelation factor decreases as $\tau$ increases, and
is roughly $0.8$ for $\tau=1$, for all probes considered.

(a)

(b) ${\rm Re}=15000$, $x_{1}^{\rm probe}$

(c) ${\rm Re}=15000$, $x_{2}^{\rm probe}$

(d) ${\rm Re}=15000$, $x_{3}^{\rm probe}$

(e)

(f) ${\rm Re}=25000$, $x_{1}^{\rm probe}$

(g) ${\rm Re}=25000$, $x_{2}^{\rm probe}$

(h) ${\rm Re}=25000$, $x_{3}^{\rm probe}$

Figure 14: A lid-driven cavity problem. Behavior of the velocity components at
three spatial locations, for two values of ${\rm Re}$ ($x_{1}^{\rm
probe}=[0,-0.8]$, $x_{2}^{\rm probe}=[0,-0.95]$, $x_{3}^{\rm probe}=[0.8,0]$).

(a) ${\rm Re}=1510^{3}$, $x_{1}^{\rm probe}$

(b) ${\rm Re}=1510^{3}$, $x_{2}^{\rm probe}$

(c) ${\rm Re}=2510^{3}$, $x_{1}^{\rm probe}$

(d) ${\rm Re}=2510^{3}$, $x_{2}^{\rm probe}$

Figure 15: A lid-driven cavity problem. Behavior of the autocorrelation for
the velocity components at two spatial locations, for two values of ${\rm Re}$
($x_{1}^{\rm probe}=[0,-0.8]$, $x_{2}^{\rm probe}=[0,-0.95]$).

##  Appendix B _A posteriori_ assessment of the POD accuracy

As explained in the main body of the paper, POD relies on a snapshot set
$\\{\mathring{u}^{k}\\}_{k=1}^{K}$ to generate a $N$-dimensional approximation
space for the full trajectory $\\{\mathring{u}^{j}\\}_{j=J_{0}}^{J}$ in the
limit $J\to\infty$. The snapshot set is associated with the sampling times
$\\{t_{\rm s}^{k}:=T_{0}+\Delta t_{\rm s}k\\}_{k=1}^{K}$, where $T_{0}=t_{\rm
g}^{J_{0}}$, $K$ is the cardinality of the snapshot set, and $\Delta t_{\rm
s}$ is the sampling period.

The choices of $\Delta t_{\rm s}$ and $K$ are a trade-off between (i)
information content of the snapshot set, and (ii) computational resources. The
snapshot set should be rich enough to accurately estimate the first $N$ POD
modes associated with the full trajectory
$\\{\mathring{u}^{j}\\}_{j=J_{0}}^{J}$ in the limit $J\to\infty$. On the other
hand, it is well-known that POD suffers from (i) the quadratic growth in $K$
in computational complexity for computing the Gramian, and for computing the
symmetric eigen-decomposition; and (ii) the memory requirements related to the
storage of the snapshots, which scale linearly with $K$. Furthermore, by
increasing $\Delta t_{\rm s}$ and $K$, we ultimately increase the number of
time steps performed by the spectral element solver — which is given by
$J=\left(T_{0}+\Delta t_{\rm s}K\right)/\Delta t$.

In this Appendix, we propose a cross-validation (CV, see, e.g., [76] and [77,
Chapter 7.10]) strategy to estimate the $\ell^{2}$-averaged projection error
associated with the POD reduced space over the full-trajectory,

$\mathcal{E}(\\{\mathring{u}^{j}\\}_{j=J_{0}}^{J},\mathcal{Z}^{\rm
u})=\frac{1}{J+1-J_{0}}\,\sum_{j=J_{0}}^{J}\,\|\mathring{u}^{j}-\Pi_{\mathcal{Z}^{\rm
u}}^{V}\mathring{u}^{j}\|_{V}^{2}.$

Estimates of this quantity might be employed to decide whether or not to
acquire new snapshots and/or to increase the dimension $N$ of the reduced
space. On the other hand, evaluations of the autocorrelation factor introduced
in Appendix A can be used to assess _a posteriori_ the amount of redundancy in
the snapshot set.

Since in our setting the snapshots are correlated in time (cf. Appendix A), we
here rely on the $h$-block variant proposed in [78] (see also [79]). The
approach relies on the assumption that the snapshot set is associated with a
stationary process: under this assumption, the covariance matrix between
$\mathring{u}^{j}$ and $\mathring{u}^{j+\kappa}$ is only a function of
$\kappa$, and approaches $0$ as $\kappa\to\infty$. The key idea of $h$-block
CV is to reduce the training set by removing the $h$ observations preceding
and following the observation in the test set. In section B.1, we adapt the
computational procedure discussed in [78] to the particular learning task of
interest; then, in section B.2, we apply the procedure to the case ${\rm
Re}=15000$ to support our choice $K=500$.

Before proceeding with the presentation of the methodology, we remark that, in
the statistics literature, several authors have proposed validation techniques
to assess the accuracy of POD (or, equivalently, PCA and Karhunen-Loéve)
spaces. We refer to [80, Chapter 6] and to the references therein for a number
of different proposals. We further recall the work by Chowdhary and Najm [81]
that relies on a Bayesian framework to account for inaccuracies due to limited
sample size. The approaches presented in [80, 81] aim at generating confidence
(credible) regions for the estimate of the POD modes; on the other hand, we
are here primarily interested in assessing the _out-of-sample_ accuracy of the
$N$-dimensional POD reduced space $\mathcal{Z}^{\rm u}$ for the full
trajectory $\\{\mathring{u}^{j}\\}_{j=J_{0}}^{J}$. For completeness, we also
recall that several authors ([82, 83]) have proposed and analyzed hierarchical
POD approaches to reduce the size $K$ of the snapshot set, without
significantly compromising the accuracy of the POD space.

###  B.1 $h$-block Cross-Validation

Algorithm 4 summarizes the computational procedure for the estimation of the
$\ell^{2}$-averaged projection error
$\mathcal{E}(\\{\mathring{u}^{j}\\}_{j=J_{0}}^{J},\mathcal{Z}^{\rm u})$. We
observe that for $h=0$ the procedure reduces to Leave-One-Out-Cross-Validation
(LOOCV). We further observe that the approach requires the assembling of the
Gramian matrix $\mathbb{U}$, and then the solution to $K$ dense eigenvalue
problems of size $K-2h-1$: for the particular problem at hand, the
computational cost associated with the procedure is negligible compared to the
computational cost associated with the solution to the FOM. Finally, we
emphasize that the procedure relies on the input parameter $h$. We here
propose to choose $h$ based on the analysis of the autocorrelation factor:
recalling the results presented in Appendix A, we consider $h=\frac{4}{\Delta
t_{s}}=4$.

Algorithm 4 $h$-block Cross-Validation

$[\widehat{\mathcal{E}}]=$ hblock-CV $(\\{\mathring{u}^{k}\\}_{k=1}^{K},h,N)$

_Inputs:_ $\\{\mathring{u}^{k}\\}_{k=1}^{K}$ = snapshot set, $h$ = correlation
parameter, $N$ = size of the POD space.

_Output:_ $\widehat{\mathcal{E}}=$ CV estimate of
$\mathcal{E}(\\{\mathring{u}\\}_{j=J_{0}}^{J},\mathcal{Z}^{\rm u})$.

1:for $k=1,\ldots,K$ do

2: $[\mathcal{Z}^{\rm u,(k)}:={\rm
span}\\{\zeta_{n}^{(k)}\\}_{n=1}^{N}]=\texttt{POD}_{V}\left(\\{\mathring{u}^{1},\ldots,\mathring{u}^{k-h-1},\mathring{u}^{k+h+1},\ldots,\mathring{u}^{K}\\}_{k=1}^{K},N\right)$

3:end for

4:Compute the CV estimate $\widehat{\mathcal{E}}$ as
$\widehat{\mathcal{E}}=\frac{1}{K}\sum_{k=1}^{K}\,\|\mathring{u}^{k}-\Pi_{\mathcal{Z}^{\rm
u,(k)}}^{V}\mathring{u}^{k}\|_{V}^{2}.$

###  B.2 Results for ${\rm Re}=15000$

Figure 16 shows the behavior of $\widehat{\mathcal{E}}$ for different values
of $N$ for the snapshot set $\\{\mathring{u}^{k}\\}_{k=1}^{K}$ associated with
$\\{t_{\rm s}^{k}=500+k\\}_{k=1}^{K=500}$, and ${\rm Re}=15000$. We compare
results with the in-sample estimate

$\mathcal{E}^{\rm
in}=\frac{1}{K}\sum_{k=1}^{K}\,\|\mathring{u}^{k}-\Pi_{\mathcal{Z}^{\rm
u}}^{V}\mathring{u}^{k}\|_{V}^{2},$

and the out-of-sample estimate

$\mathcal{E}^{\rm
out}=\frac{1}{K}\sum_{k=K+1}^{2K}\,\|\mathring{u}^{k}-\Pi_{\mathcal{Z}^{\rm
u}}^{V}\mathring{u}^{k}\|_{V}^{2},$

where $\\{\mathring{u}^{k}\\}_{k=K+1}^{2K}$ are associated with the sampling
times $\\{t_{\rm s}^{k}=500+k\\}_{k=K+1}^{2K}$. For visualization purposes, we
normalize $\widehat{\mathcal{E}}$, $\mathcal{E}^{\rm in}$, and
$\mathcal{E}^{\rm out}$ by $\mathcal{E}^{\rm in}(N=1)$: for $N=60$,
$\widehat{\mathcal{E}}\approx 15\%\times\mathcal{E}^{\rm in}(N=1)$,
$\mathcal{E}^{\rm in}\approx 10\%\times\mathcal{E}^{\rm in}(N=1)$, and
$\mathcal{E}^{\rm out}\approx 17\%\times\mathcal{E}^{\rm in}(N=1)$. We observe
that $\widehat{\mathcal{E}}$ is a more accurate estimate of $\mathcal{E}^{\rm
out}$ compared to the in-sample estimate $\mathcal{E}^{\rm in}$.

Figure 16: A Cross-Validation procedure for the _a posteriori_ assessment of
the POD accuracy. Behavior of $\mathcal{E}^{\rm in}$, $\widehat{\mathcal{E}}$,
and $\mathcal{E}^{\rm out}$ with $N$. All quantities are normalized by
$\mathcal{E}^{\rm in}(N=1)$ ($K=500$, $h=2$, ${\rm Re}=15000$).

##  Appendix C On the definition of ROM stability

Based on the results of section 3, we could take a pragmatic view of the long-
time ROM stability. Given the reduced space $\mathcal{Z}^{\rm u}\subset V_{\rm
div}$, we define the best-fit errors associated with mean flow and TKE:

$e_{1}^{\rm opt}:=\frac{\|\langle\mathring{u}\rangle_{\rm
g}-\Pi_{\mathcal{Z}^{\rm u}}^{V}\langle\mathring{u}\rangle_{\rm
g}\|_{V}}{\|\langle u\rangle_{\rm g}\|_{V}},\quad e_{2}^{\rm
opt}:=\frac{|\langle{\rm TKE}\rangle_{\rm s}-\langle{\rm TKE}^{\rm
opt}(\cdot;\mathcal{Z}^{\rm u})\rangle_{\rm s}|}{\langle{\rm TKE}\rangle_{\rm
s}},$

where $\Pi_{\mathcal{Z}^{\rm u}}^{V}$ denotes the projection operator
associated with the $V$ inner product on the subspace $\mathcal{Z}^{\rm u}$,
and ${\rm TKE}^{\rm opt}(t;\mathcal{Z}^{\rm
u})=\frac{1}{2}\int_{\Omega}\|\Pi_{\mathcal{Z}^{\rm
u}}^{L^{2}}(\mathring{u}(t)-\langle\mathring{u}\rangle_{\rm
g})\|_{2}^{2}\,dx$. Then, we define the effective stability constants as the
ratios between the optimal mean error and the actual error:

${\rm m}(\mathcal{Z}^{\rm u}):=\frac{\|\langle\mathring{u}\rangle_{\rm
g}-\langle\hat{u}\rangle_{\rm g}\|_{V}}{\|\langle u\rangle_{\rm
g}\|_{V}e_{1}^{\rm opt}},\qquad{\rm\sigma}(\mathcal{Z}^{\rm
u}):=\frac{|\langle{\rm TKE}\rangle_{\rm s}-\langle\widehat{{\rm
TKE}}\rangle_{\rm s}|}{\langle{\rm TKE}\rangle_{\rm s}e_{2}^{\rm opt}}.$ (21)

The stability constants ${\rm m}$ and ${\rm\sigma}$ can be used to
quantitatively measure the stability of the ROM. We observe that, in the limit
$T\to\infty$, our definition of long-time stability for ROMs is independent of
transient dynamics. We further observe that a ROM of dimension $N$ is stable
if and only if mean and variance of the time coefficients
$\\{a_{n}^{j}\\}_{j}$ are correctly estimated for $n=1,\ldots,N$. We finally
remark that our definition of stability is close to the one proposed in [1];
however, while the definition in [1] is tailored to $L^{2}$ POD spaces and
$R_{g}=\langle u\rangle_{\rm g}$, our definition applies to any reduced space
and to any choice of the lift.

Figure 17 shows the behavior of ${\rm m}(\mathcal{Z}^{\rm u})$ and
${\rm\sigma}(\mathcal{Z}^{\rm u})$ defined in (21) for POD-Galerkin and
constrained POD-Galerkin for ${\rm Re}=15000$: our constrained POD-Galerkin
ROM is more stable — according to the definition given in this Appendix — than
the standard Galerkin ROM.

(a)

(b)

Figure 17: The solution reproduction problem; behavior of ${\rm
m}(\mathcal{Z}^{\rm u})$ and ${\rm\sigma}(\mathcal{Z}^{\rm u})$ (21) for POD-
Galerkin and constrained POD-Galerkin (${\rm Re}=15000$, $\epsilon=0.01$).

##  Appendix D Robustness of the constrained formulation

In this Appendix, we present a number of numerical results that provide
further insights about the constrained formulation proposed in this paper. In
greater detail, we study the activation rate of the box constraints, the
dependence of the solution to the choice of $\epsilon$, and the behavior of
$m_{n}^{\rm FOM,u}$ and $M_{n}^{\rm FOM,u}$ defined in (17) with respect to
the Reynolds number ${\rm Re}$.

In Figure 18, we study the behavior of the activation rate of each box
constraint for two values of $N$ for ${\rm Re}=15000$, $\\{t_{\rm
s}^{k}=500+k\\}_{k=1}^{K=500}$. In greater detail, we count how many times the
$n$-th component of the solution to Galerkin satisfies the prescribed
constraints:

$\mbox{\\#}{\rm
Gal}_{n}:=\frac{1}{J-J_{0}}\sum_{j=J_{0}+1}^{J}\,\mathbbm{1}\left((\mathbf{a}_{\rm
Gal}^{j})_{n}\in[\alpha_{n},\beta_{n}]\right),\quad n=1,\ldots,N.$

We observe that the behavior with $n$ of $\mbox{\\#}{\rm Gal}_{n}$ is
irregular, and strongly depends on $N$. This suggests that selecting _a
priori_ the active constraints might be impractical. In Figure 19, we study
the behavior of the relative error in the mean flow prediction, the behavior
of the mean TKE, and the behavior of the activation rate of the box
constraints

$\mbox{\\#}{\rm
Gal}:=\frac{1}{J-J_{0}}\sum_{j=J_{0}+1}^{J}\,\mathbbm{1}\left((\mathbf{a}_{\rm
Gal}^{j})_{n}\in[\alpha_{n},\beta_{n}],\,n=1,\ldots,N\right),$

with respect to $\epsilon$, for two values of $N$. We observe that for
$\epsilon\lesssim\bar{\epsilon}=0.1$ results do not seem to depend on the
value of $\epsilon$. This provides evidence that the current approach is
robust with respect to the choice of $\epsilon$. We further observe that for
all values of $\epsilon$ considered $\mbox{\\#}{\rm Gal}(N=40)\gtrsim 0.85$
and $\mbox{\\#}{\rm Gal}(N=60)\gtrsim 0.90$. Therefore, our constrained
formulation corrects the original formulation only for $10-15\%$ time steps.
For this reason, we envision that efficient implementations of the constrained
ROM might be nearly as inexpensive as the Galerkin ROM. We further observe
that $\mbox{\\#}{\rm Gal}$ increases as $N$ increases: this can be explained
by observing that the POD-Galerkin ROM becomes more and more accurate as $N$
increases, and thus requires less corrections.

(a) $N=40$

(b) $N=60$

Figure 18: The solution reproduction problem; activity of the box constraints
for two values of $N$. (${\rm Re}=15000$, $\epsilon=0.01$).

(a)  $N=40$

(b)  $N=40$

(c)  $N=40$

(d)  $N=60$

(e)  $N=60$

(f)  $N=60$

Figure 19: The solution reproduction problem; sensitivity with respect to
$\epsilon$ for constrained POD-Galerkin for two values of $N$. Figures (a) and
(d): behavior of the relative $L^{2}$ and $H^{1}$ errors. Figures (b) and (e):
behavior of the mean TKE. Figures (c) and (f): percentage of pure Galerkin
solves. (${\rm Re}=15000$).

Figure 20 investigates the behavior of $m_{n}^{\rm FOM,u}$ and $M_{n}^{\rm
FOM,u}$ defined in (17) with respect to the Reynolds number ${\rm Re}$. For
this test, we consider the POD space associated with ${\rm Re}=20000$ and the
sampling times $\\{t_{\rm s}^{k}=500+k\\}_{k=1}^{K=1000}$, and we show results
for $n=1,\ldots,12$. Results suggest that the sensitivity of $m_{n}^{\rm
FOM,u}$ and $M_{n}^{\rm FOM,u}$ with ${\rm Re}$ are relatively modest if
compared to $M_{n}^{\rm FOM,u}-m_{n}^{\rm FOM,u}$.

Figure 20: The parametric problem; behavior of $m_{n}^{\rm FOM,u}({\rm Re})$
and $M_{n}^{\rm FOM,u}({\rm Re})$ with ${\rm Re}$. The POD space is generated
from the DNS data for ${\rm Re}=20000$.

##  Appendix E On the problem of $p$-refinement

We here illustrate the major issue associated with the combination of POD
modes associated with different values of the parameter. We here simulate the
application of the first two iterations of the POD-$p$Greedy algorithm as
proposed in [44]. In more detail, we consider the following test.

  1. 1.

Generate DNS data for ${\rm Re}=20000$, and use them to build the
$N_{1}=60$-dimensional POD space.

  2. 2.

Generate DNS data for ${\rm Re}=15000$, and build the $N_{2}=60$-dimensional
POD space for the set of snapshots $\\{(\mathbb{I}-\Pi_{\mathcal{Z}^{\rm
u}}^{V})\mathring{u}^{k}({\rm Re}=15000)\\}_{k=1}^{K}$.

  3. 3.

Define $\mathcal{Z}^{\rm u}={\rm span}\\{\zeta_{n}\\}_{n=1}^{N_{1}+N_{2}}$
where $\zeta_{1},\ldots,\zeta_{N_{1}}$ are associated with Step 1 and
$\zeta_{N_{1}+1},\ldots,\zeta_{N_{2}}$ are associated with Step 2.

  4. 4.

Perform a convergence study in $N$ for ${\rm Re}=15000$ and ${\rm Re}=20000$
for both pure Galerkin and constrained Galerkin ($\epsilon=0.01$).

We consider here $T_{0}=500$ and $T=1500$, $\\{t_{\rm
s}^{k}=500+k\\}_{k=1}^{K=1000}$. We recall that for ${\rm Re}=15000$ (cf.
section 3) we were able to obtain accurate ROMs for $N\gtrsim 40$ both in
terms of mean flow prediction and TKE.

Figures 21 and 22 show the results of this test for the constrained-Galerkin
ROM. Figures 21(a) and (c) show the behavior of the relative error in mean
flow prediction for the constrained formulation, for ${\rm Re}=15000$ and
${\rm Re}=20000$, respectively. We here compute lower and upper bounds
$\\{\alpha_{n}\\}_{n}$ and $\\{\beta_{n}\\}_{n}$ using (13) with
$\epsilon=0.01$. Figures 21(b) and (d) show the behavior of the mean TKE for
the same values of the Reynolds number. Similarly, Figures 22(a) and (b) show
the behavior of the TKE in time for $N=120$. Finally, Figures 23 (a) and (b)
show the behavior of the TKE in time for $N=120$ for the unconstrained
formulation. Results — especially for ${\rm Re}=15000$ — show the key issue of
combining modes associated with different parameters. For $N_{2}=60$ (and
$N_{1}+N_{2}=120$), the error in mean flow prediction is roughly $10\%$, and
we also significantly overestimate the mean and the peaks of the TKE. As
expected, these issues are even more severe for the unconstrained formulation:
the behavior with time of the TKE predicted by the unconstrained ROM is
roughly the same for the two values of the Reynolds number considered.

We offer a physical explanation for the poor performance of POD-$p$Greedy. As
observed in Appendix A, for sufficiently large values of ${\rm Re}$ eddies are
ejected into the core region of the cavity. This instability is observed for
${\rm Re}=20000$, but is not observed for ${\rm Re}=15000$. As a result, the
ejection of the eddies into the core region of the cavity is well-represented
by the POD space associated with ${\rm Re}=20000$, and then, by construction,
by the final reduced space $\mathcal{Z}^{\rm u}$. The presence of modes
associated with the core eddies makes the ROM more prone to show this
instability even for values of the Reynolds number at which the full-order
solution does not show it.

(a) ${\rm Re}=15000$

(b) ${\rm Re}=15000$

(c) ${\rm Re}=20000$

(d) ${\rm Re}=20000$

Figure 21: The parametric problem; on the problems of $p$-refinement.
Performance of constrained Galerkin ($\epsilon=0.01$, $N_{1}=60$, $N_{2}=60$).

(a) ${\rm Re}=15000$, $N=120$

(b) ${\rm Re}=20000$, $N=120$

Figure 22:  The parametric problem; on the problems of $p$-refinement.
Behavior of the TKE with time for constrained Galerkin ($\epsilon=0.01$,
$N_{1}=60$, $N_{2}=60$).

(a) ${\rm Re}=15000$, $N=120$

(b) ${\rm Re}=20000$, $N=120$

Figure 23:  The parametric problem; on the problems of $p$-refinement.
Behavior of the TKE with time for unconstrained Galerkin ($N_{1}=60$,
$N_{2}=60$).

##  Appendix F Offline/online computational decomposition for the residual
indicator

We here describe the offline/online computational decomposition for the
computation of the residual error indicator introduced in this paper. We omit
the subscript $\ell$ associated with the partition of the parameter domain to
simplify notation. We first introduce the Riesz representers:

$\begin{array}[]{l}(\xi_{n}^{\rm
m},v)_{V}=(\zeta_{n},v)_{L^{2}(\Omega)},\quad(\xi_{n}^{\rm
a},v)_{V}=(\zeta_{n},v)_{V},\quad(\xi_{n}^{\rm
cg},v)_{V}=c(R_{g},\zeta_{n},v),\\\\[8.53581pt] (\xi_{m,n}^{\rm
c},v)_{V}=c(\zeta_{n},\zeta_{m},v),\quad(\xi_{n}^{\rm
mg},v)_{V}=c(\zeta_{n},R_{g},v),\quad(\xi_{1}^{\rm
f},v)_{V}=(R_{g},v)_{V},\\\\[8.53581pt] (\xi_{2}^{\rm
f},v)_{V}=c(R_{g},R_{g},v),\\\ \end{array}$ (22)

for $n=1,\ldots,N$ and for all $v\in V_{\rm div}$. Then, it is easy to verify
that, if $w^{j}=\sum_{n=1}^{N}a_{n}^{j}\zeta_{n}$ for $j=0,\ldots,J$, we can
rewrite $\langle R\rangle$ as follows:

$\begin{array}[]{ll}\displaystyle{\langle
R\rangle\left(\\{w^{j}\\}_{j=0}^{J},v;\,{\rm
Re}\right)=}&\displaystyle{\Big{(}\sum_{n=1}^{N}\Big{(}\xi_{n}^{\rm
m}\left(\frac{a_{n}^{J}-a_{n}^{J_{0}}}{T-T_{0}}\right)+\xi_{n}^{\rm
a}\left(\frac{1}{\rm Re}\bar{a}_{n}^{+}\right)+\xi_{n}^{\rm
cg}\bar{a}_{n}^{+}}\\\\[11.38109pt]
&\displaystyle{+\sum_{m=1}^{N}\xi_{m,n}^{\rm c}\bar{c}_{m,n}+\xi_{n}^{\rm
mg}\bar{a}_{n}^{-}\Big{)},\,v\Big{)}_{V}}\\\\[11.38109pt]
&\displaystyle{+\frac{1}{\rm Re}(\xi_{1}^{\rm f},v)_{V}+(\xi_{2}^{\rm
f},v)_{V};}\\\ \end{array}$ (23a) where $\bar{a}_{n}^{+}=\frac{\Delta
t}{T-T_{0}}\sum_{j=J_{0}+1}^{J}\,a_{n}^{j},\quad\bar{a}_{n}^{-}=\frac{\Delta
t}{T-T_{0}}\sum_{j=J_{0}}^{J-1}\,a_{n}^{j},\quad\bar{c}_{m,n}=\frac{\Delta
t}{T-T_{0}}\,\sum_{j=J_{0}}^{J-1}a_{m}^{j+1}a_{n}^{j}.$ (23b)

Equation (23) can be rewritten as $\langle
R\rangle\left(\\{w^{j}\\}_{j=0}^{J},v;\,{\rm
Re}\right)=\sum_{m=1}^{M}\,\Theta_{i}(\\{\mathbf{a}^{j}\\}_{j};{\rm
Re})\tilde{\xi}_{m}$ where $M=N^{2}+3N+2$ and

$\begin{array}[]{l}\displaystyle{[\tilde{\xi}_{1},\ldots,\tilde{\xi}_{M}]=[\xi_{1}^{\rm
m},\ldots,\xi_{N}^{\rm m},\xi_{1}^{\rm a},\ldots,\xi_{N}^{\rm a},\xi_{1}^{\rm
cg},\ldots,\xi_{N}^{\rm cg},}\\\\[5.69054pt]
\hfill\displaystyle{\xi_{1,1}^{\rm c},\ldots,\xi_{N,N}^{\rm c},\xi_{1}^{\rm
mg},\ldots,\xi_{N}^{\rm mg},\xi_{1}^{\rm f},\xi_{2}^{\rm f}]}\\\\[11.38109pt]
\displaystyle{[\Theta_{1},\ldots,\Theta_{M}]=\Big{[}\frac{a_{1}^{J}-a_{1}^{J_{0}}}{T-T_{0}},\ldots,\frac{a_{N}^{J}-a_{N}^{J_{0}}}{T-T_{0}},\frac{\bar{a}_{1}^{+}}{\rm
Re},\ldots,\frac{\bar{a}_{N}^{+}}{\rm
Re},\bar{a}_{1}^{+},\ldots,\bar{a}_{N}^{+},}\\\\[8.53581pt]
\hfill\displaystyle{\bar{c}_{1,1},\ldots,\bar{c}_{N,N},\bar{a}_{1}^{-},\ldots,\bar{a}_{N}^{-},\frac{1}{\rm
Re},1\Big{]}.}\end{array}$

Therefore, recalling the Riesz representation theorem, we find

$\Delta^{\rm u}(\\{w^{j}\\}_{j};{\rm
Re})=\sqrt{\boldsymbol{\Theta}^{T}\Sigma\boldsymbol{\Theta}},\qquad\boldsymbol{\Theta}=\boldsymbol{\Theta}(\\{\mathbf{a}^{j}\\}_{j};{\rm
Re}),$ (24)

where $\Sigma_{i,i^{\prime}}=(\tilde{\xi}_{i},\tilde{\xi}_{i^{\prime}})_{V}$.
Equation (24) clarifies the offline/online decomposition: during the offline
stage, we compute the Riesz representers (22) — this corresponds to the
solution to $M$ Stokes problems — and we assemble the matrix $\Sigma$; during
the online stage, we compute the vector $\boldsymbol{\Theta}$ and we exploit
(24) to compute the error estimator $\Delta^{\rm u}$.

## Acknowledgements

The authors thank Prof. Paul Fischer (UIUC), and Dr. Elia Merzari (Argonne
National Lab) for their support with the software Nek5000.

## References

  * [1] M Balajewicz and EH Dowell.  Stabilization of projection-based reduced order models of the navier–stokes.  Nonlinear Dynamics, 70(2):1619–1632, 2012.
  * [2] Z Wang, I Akhtar, J Borggaard, and T Iliescu.  Proper orthogonal decomposition closure models for turbulent flows: a numerical comparison.  Computer Methods in Applied Mechanics and Engineering, 237:10–26, 2012.
  * [3] L Cordier, B R Noack, G Tissot, G Lehnasch, J Delville, M Balajewicz, G Daviller, and R K Niven.  Identification strategies for model-based control.  Experiments in fluids, 54(8):1–21, 2013.
  * [4] AE Deane, IG Kevrekidis, G E Karniadakis, and SA Orszag.  Low-dimensional models for complex geometry flows: Application to grooved channels and circular cylinders.  Physics of Fluids A: Fluid Dynamics (1989-1993), 3(10):2337–2354, 1991.
  * [5] X Ma, G-S Karamanos, and GE Karniadakis.  Dynamics and low-dimensionality of a turbulent near wake.  Journal of Fluid Mechanics, 410:29–65, 2000.
  * [6] X Ma and G E Karniadakis.  A low-dimensional model for simulating three-dimensional cylinder flow.  Journal of Fluid Mechanics, 458:181–190, 2002.
  * [7] B Galletti, CH Bruneau, L Zannetti, and A Iollo.  Low-order modelling of laminar flow regimes past a confined square cylinder.  Journal of Fluid Mechanics, 503:161–170, 2004.
  * [8] G Berkooz, P Holmes, and J L Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics, 25(1):539–575, 1993.
  * [9] K Kunisch and S Volkwein.  Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics.  SIAM Journal on Numerical analysis, 40(2):492–515, 2002.
  * [10] M Kahlbacher and S Volkwein.  Galerkin proper orthogonal decomposition methods for parameter dependent elliptic systems.  Discussiones Mathematicae, Differential Inclusions, Control and Optimization, 27(1):95–117, 2007.
  * [11] 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.
  * [12] D Rempfer.  On low-dimensional galerkin models for fluid flow.  Theoretical and Computational Fluid Dynamics, 14(2):75–88, 2000\.
  * [13] P Moin and K Mahesh.  Direct numerical simulation: a tool in turbulence research.  Annual review of fluid mechanics, 30(1):539–578, 1998.
  * [14] B R Noack, P Papas, and P A Monkewitz.  The need for a pressure-term representation in empirical galerkin models of incompressible shear flows.  Journal of Fluid Mechanics, 523:339–365, 2005.
  * [15] N Aubry, P Holmes, J L Lumley, and E Stone.  The dynamics of coherent structures in the wall region of a turbulent boundary layer.  Journal of Fluid Mechanics, 192:115–173, 1988.
  * [16] P Sagaut.  Large eddy simulation for incompressible flows: an introduction.  Springer Science & Business Media, 2006.
  * [17] M Couplet, P Sagaut, and C Basdevant.  Intermodal energy transfers in a proper orthogonal decomposition–galerkin representation of a turbulent separated flow.  Journal of Fluid Mechanics, 491:275–284, 2003.
  * [18] B R Noack, K Afanasiev, M Morzynski, G Tadmor, and F Thiele.  A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.  Journal of Fluid Mechanics, 497:335–363, 2003.
  * [19] S Sirisup and GE Karniadakis.  A spectral viscosity method for correcting the long-term behavior of pod models.  Journal of Computational Physics, 194(1):92–116, 2004.
  * [20] Y Maday, A Manzoni, and A Quarteroni.  An online intrinsic stabilization strategy for the reduced basis approximation of parametrized advection-dominated problems.  Comptes Rendus Mathematique, 354(12):1188–1194, 2016.
  * [21] E Tadmor.  Convergence of spectral methods for nonlinear conservation laws.  SIAM Journal on Numerical Analysis, 26(1):30–44, 1989.
  * [22] S Giere, T Iliescu, J Volker, and D Wells.  Supg reduced order models for convection-dominated convection–diffusion–reaction equations.  Computer Methods in Applied Mechanics and Engineering, 289:454–474, 2015.
  * [23] X Xie, D Wells, Z Wang, and T Iliescu.  Approximate deconvolution reduced order modeling.  Computer Methods in Applied Mechanics and Engineering, 313:512–534, 2017.
  * [24] O San and T Iliescu.  Proper orthogonal decomposition closure models for fluid flows: Burgers equation.  International Journal of Numerical Analysis & Modeling - Series B, 5(3):217–237, 2014.
  * [25] M Marion and R Temam.  Nonlinear galerkin methods: the finite elements case.  Numerische Mathematik, 57(1):205–226, 1990.
  * [26] A Debussche, T Dubois, and R Temam.  The nonlinear galerkin method: a multiscale method applied to the simulation of homogeneous turbulent flows.  Theoretical and Computational Fluid Dynamics, 7(4):279–315, 1995\.
  * [27] M Bergmann, C-H Bruneau, and A Iollo.  Enablers for robust pod models.  Journal of Computational Physics, 228(2):516–538, 2009.
  * [28] M Bergmann.  Optimisation aérodynamique par réduction de modèle POD et contrôle optimal: application au sillage laminaire d’un cylindre circulaire.  PhD thesis, Vandoeuvre-les-Nancy, INPL, 2004.
  * [29] Y Maday, A T Patera, and D V Rovas.  A blackbox reduced-basis output bound method for noncoercive linear problems.  in Studies in Mathematics and its Applications, D. Cioranescu and J. L. Lions, eds., Elsevier Science B. V, pages 533–569, 2001.
  * [30] 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(2):155–181, 2011.
  * [31] 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.
  * [32] A Tallet, C Allery, C Leblond, and E Liberge.  A minimum residual projection to build coupled velocity–pressure pod–rom for incompressible navier–stokes equations.  Communications in Nonlinear Science and Numerical Simulation, 22(1):909–932, 2015.
  * [33] A Iollo, S Lanteri, and J A Désidéri.  Stability properties of pod–galerkin approximations for the compressible navier–stokes equations.  Theoretical and Computational Fluid Dynamics, 13(6):377–396, 2000\.
  * [34] N Aubry, R Guyonnet, and R Lima.  Spatio-temporal symmetries and bifurcations via bi-orthogonal decompositions.  Journal of Nonlinear Science, 2(2):183–215, 1992.
  * [35] M F Barone, I Kalashnikova, D J Segalman, and H K Thornquist.  Stable galerkin reduced order models for linearized compressible flow.  Journal of Computational Physics, 228(6):1932–1946, 2009.
  * [36] 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.
  * [37] M Couplet, C Basdevant, and P Sagaut.  Calibrated reduced-order pod-galerkin system for fluid flow modelling.  Journal of Computational Physics, 207(1):192–220, 2005.
  * [38] P J Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of fluid mechanics, 656:5–28, 2010.
  * [39] C W Rowley, I Mezić, S Bagheri, P Schlatter, and D S Henningson.  Spectral analysis of nonlinear flows.  Journal of fluid mechanics, 641:115–127, 2009.
  * [40] P J Schmid, L Li, M P Juniper, and O Pust.  Applications of the dynamic mode decomposition.  Theoretical and Computational Fluid Dynamics, 25(1-4):249–259, 2011\.
  * [41] A Alla and J N Kutz.  Nonlinear model order reduction via dynamic mode decomposition.  arXiv preprint arXiv:1602.05080, 2016.
  * [42] J H Curry, J R Herring, J Loncaric, and S A Orszag.  Order and disorder in two-and three-dimensional benard convection.  Journal of Fluid Mechanics, 147:1–38, 1984.
  * [43] G Stabile, S Hijazi, A Mola, S Lorenzi, and G Rozza.  Advances in reduced order modelling for cfd: vortex shedding around a circular cylinder using a pod-galerkin method.  arXiv preprint arXiv:1701.03424, 2017.
  * [44] B Haasdonk and M Ohlberger.  Reduced basis method for finite volume approximations of parametrized linear evolution equations.  ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique, 42(2):277–302, 2008.
  * [45] J L Eftang, D J Knezevic, and A T Patera.  An hp certified reduced basis method for parametrized parabolic partial differential equations.  Mathematical and Computer Modelling of Dynamical Systems, 17(4):395–422, 2011.
  * [46] T Taddei, S Perotto, and A Quarteroni.  Reduced basis techniques for nonlinear conservation laws.  ESAIM: Mathematical Modelling and Numerical Analysis, 49(3):787–814, 2015.
  * [47] N Cagniart, Y Maday, and B Stamm.  Model order reduction for problems with large convection effects.  2016\.
  * [48] M Ohlberger and S Rave.  Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing.  Comptes Rendus Mathematique, 351(23):901–906, 2013.
  * [49] K Urban and A T Patera.  A new error bound for reduced basis approximation of parabolic partial differential equations.  Comptes Rendus Mathematique, 350(3):203–207, 2012.
  * [50] M Yano.  A space-time petrov–galerkin certified reduced basis method: Application to the boussinesq equations.  SIAM Journal on Scientific Computing, 36(1):A232–A266, 2014.
  * [51] JP Argaud, B Bouriquet, H Gong, Y Maday, and O Mula.  Stabilization of (g) eim in presence of measurement noise: application to nuclear reactor physics.  arXiv preprint arXiv:1611.02219, 2016.
  * [52] B Haasdonk.  Convergence rates of the pod–greedy method.  ESAIM: Mathematical Modelling and Numerical Analysis, 47(03):859–873, 2013.
  * [53] M Balajewicz.  A New Approach to Model Order Reduction of the Navier-Stokes Equations.  PhD thesis, Duke University, 2012.
  * [54] W Cazemier, RWCP Verstappen, and AEP Veldman.  Proper orthogonal decomposition and low-dimensional models for driven cavity flows.  Physics of Fluids (1994-present), 10(7):1685–1699, 1998.
  * [55] PN Shankar and MD Deshpande.  Fluid mechanics in the driven cavity.  Annual Review of Fluid Mechanics, 32(1):93–136, 2000.
  * [56] F Terragni, E Valero, and J M Vega.  Local pod plus galerkin projection in the unsteady lid-driven cavity problem.  SIAM Journal on Scientific Computing, 33(6):3538–3561, 2011.
  * [57] S Lorenzi, A Cammi, L Luzzi, and G Rozza.  Pod-galerkin method for finite volume approximation of navier–stokes and rans equations.  Computer Methods in Applied Mechanics and Engineering, 311:151–179, 2016.
  * [58] S B Pope.  Turbulent flows.  Cambridge University Press, 2000.
  * [59] A T Patera.  A spectral element method for fluid dynamics: laminar flow in a channel expansion.  Journal of Computational Physics, 54(3):468–488, 1984.
  * [60] JW Lottes PF Fischer and SG Kerkemeier.  nek5000 Web page, 2008.  http://nek5000.mcs.anl.gov.
  * [61] C Bernardi and Y Maday.  Spectral methods.  Handbook of numerical analysis, 5:209–485, 1997.
  * [62] C Canuto, M Y Hussaini, A Quarteroni, and T A Zang.  Spectral methods in fluid dynamics.  Springer Science & Business Media, 2012.
  * [63] G Karniadakis and S Sherwin.  Spectral/hp element methods for computational fluid dynamics.  Oxford University Press, 2013.
  * [64] M O Deville, P F Fischer, and E H Mund.  High-order methods for incompressible fluid flow, volume 9.  Cambridge University Press, 2002.
  * [65] A Sohankar, C Norberg, and L Davidson.  Low-reynolds-number flow around a square cylinder at incidence: study of blockage, onset of vortex shedding and outlet boundary condition.  International journal for numerical methods in fluids, 26(1):39–56, 1998.
  * [66] A E Løvgren, Y Maday, and E M Rønquist.  A reduced basis element method for the steady stokes problem.  ESAIM: Mathematical Modelling and Numerical Analysis, 40(3):529–552, 2006.
  * [67] G Rozza and K Veroy.  On the stability of the reduced basis method for stokes equations in parametrized domains.  Computer methods in applied mechanics and engineering, 196(7):1244–1260, 2007.
  * [68] G Rozza, DB P Huynh, and A Manzoni.  Reduced basis approximation and a posteriori error estimation for stokes flows in parametrized geometries: roles of the inf-sup stability constants.  Numerische Mathematik, 125(1):115–152, 2013.
  * [69] T Lassila, A Manzoni, A Quarteroni, and G Rozza.  Model order reduction in fluid dynamics: challenges and perspectives.  In Reduced Order Methods for modeling and computational reduction, pages 235–273. Springer, 2014.
  * [70] Francesco Ballarin, Gianluigi Rozza, and Yvon Maday.  Reduced-order semi-implicit schemes for fluid-structure interaction problems.  In Model Reduction of Parametrized Systems, pages 149–167. Springer, 2017.
  * [71] L Sirovich.  Turbulence and the dynamics of coherent structures. part i: Coherent structures.  Quarterly of applied mathematics, 45(3):561–571, 1987.
  * [72] S Volkwein.  Model reduction using proper orthogonal decomposition.  Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz. see http://www. uni-graz. at/imawww/volkwein/POD. pdf, 2011.
  * [73] MATLAB.  version 9.0 (R2016a).  The MathWorks Inc., Natick, Massachusetts, 2016.
  * [74] J Nocedal and S Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [75] Y Maday R Chakir.  A two-grid finite-element/reduced basis scheme for the approximation of the solution of parametric dependent p.d.e.  Comptes Rendus Mathématique, 1513(7):343–462, 2009.
  * [76] R Kohavi et al.  A study of cross-validation and bootstrap for accuracy estimation and model selection.  In Ijcai, volume 14, pages 1137–1145, 1995.
  * [77] T. Hastie, R. Tibshirani, and J. Friedman.  The elements of statistical learning, volume 2.  Springer, 2009.
  * [78] P Burman, E Chow, and D Nolan.  A cross-validatory method for dependent data.  Biometrika, 81(2):351–358, 1994.
  * [79] J Racine.  Consistent cross-validatory model-selection for dependent data: hv-block cross-validation.  Journal of econometrics, 99(1):39–61, 2000.
  * [80] Ian Jolliffe.  Principal component analysis.  Wiley Online Library, 2002.
  * [81] K Chowdhary and H N Najm.  Bayesian estimation of karhunen–loève expansions; a random subspace approach.  Journal of Computational Physics, 319:280–293, 2016.
  * [82] A Paul-Dubois-Taine and D Amsallem.  An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models.  International Journal for Numerical Methods in Engineering, 102(5):1262–1292, 2015.
  * [83] C Himpe, T Leibner, and S Rave.  Hierarchical approximate proper orthogonal decomposition.  arXiv preprint arXiv:160705210, 2016.