# Error analysis of a residual-based stabilization-motivated POD-ROM for
incompressible flows

Tomás Chacón Rebollo , Samuele Rubino , Mourad Oulghelou , Cyrille Allery
Departamento EDAN & IMUS, Universidad de Sevilla, Spain.
<EMAIL_ADDRESS>EDAN & IMUS, Universidad de Sevilla, Spain.
<EMAIL_ADDRESS>Arts et Métiers ParisTech, France.
<EMAIL_ADDRESS>Université de la Rochelle, France.
<EMAIL_ADDRESS>
###### Abstract

This article presents error bounds for a velocity-pressure segregated POD
reduced order model discretization of the Navier-Stokes equations. The
stability is proven in $L^{\infty}(L^{2})$ and energy norms for velocity, with
bounds that do not depend on the viscosity, while for pressure it is proven in
a semi-norm of the same asymptotic order as the $L^{2}$ norm with respect to
the mesh size. The proposed estimates are calculated for the two flow
problems, the flow past a cylinder and the lid-driven cavity flow. Their
quality is then assessed in terms of the predicted logarithmic slope with
respect to the velocity POD contribution ratio. We show that the proposed
error estimates allow a good approximation of the real errors slopes and thus
a good prediction of their rate of convergence.

Keywords: Navier-Stokes Equations, Residual-based Stabilization, Proper
Orthogonal Decomposition, Reduced Order Models, Incompressible Flows,
Numerical Analysis.

##  1 Introduction

Since the beginning of the XXI Century, there is an increasing number of
applications of Reduced Order Models (ROMs) addressed to the approximation of
incompressible Navier-Stokes flow solutions. By using carefully chosen basis
functions forming the projection subspace, ROMs collaborating with high
fidelity solutions, offer a computationally efficient tool to reduce the cost
of heavy numerical simulations. A very famous method that is extensively used
to construct projection subspaces is the POD (Proper Orthogonal Decomposition)
method [40, 8]. Starting from a set of ”snapshots” computed from the high-
fidelity model (or full-order model, FOM) in the off-line phase, the POD
allows trough solving an eigenvalue problem to obtain a truncated basis formed
by modes corresponding to the first larger eigenvalues. Once the POD basis is
built, the temporal dynamics of the targeted system (in our case, the Navier-
Stokes equations) is determined by solving a low dimension system obtained
from the Galerkin projection of the model equations onto this reduced subspace
[37]. The resulting model provides very fast solutions of the targeted
equations, with computational cost frequently several order of magnitude lower
than that of the FOM. The POD Galerkin approach has been applied in many
fields such as aeroelasticity [44, 33], structural dynamics [35, 2, 12], modal
analysis [24], flow prediction and control [4, 14, 22, 26, 41], particles
dispersion [3, 5, 13], etc. Other approaches to construct projection subspaces
relying on greedy algorithms can be found in [6, 36]. For a thorough review,
the interested reader can refer to [25].
In practice, the velocity snapshots are weakly divergence free up to the
computer accuracy. This allows to drop the pressure from projected equations
and form a reduced order model that only contains the velocity amplitudes as
unknown variables. However for many applications, the pressure is needed to be
recovered. Namely, to calculate lift and drag forces exerted on walls or
immersed boundaries. One option to recover the pressure is to build a ROM that
also accounts for its time trajectory in a reduced coupled velocity-pressure
setting. Specific techniques to ensure the discretisation of the pressure are
thus needed, for instance using velocity sumpremizers of the velocity-pressure
duality [7, 38] or stabilization techniques [39]. An alternative way is to
recover the pressure once the velocity is computed from the velocity ROM. The
governing reduced order model in this case is a pressure Poisson equation
derived by taking the divergence of the momentum conservation equation.
Although this sets the reduced problem for the pressure, the definition of
adequate boundary conditions remains an issue. Another way to construct the
pressure reduced model is by duality of the momentum conservation equation
with velocity supremizer test functions [29].
In this paper we afford the solution of the pressure equation by duality of
the momentum conservation equation with gradients of the reduced pressure
basis functions. This turns out to solving a Poisson equation for the
pressure, but with the advantage of setting the right boundary conditions and
thus ensuring an improved accuracy in its recovery. We shall call this method
SM (Stabilisation-Motivated) ROM, as it has been proposed in a framework of
stabilised solution of incompressible flows [15], see also [42, 32] where the
method is introduced from a minimum residual projection approach. We carry on
the numerical analysis of the SM-ROM and prove the stability in specific
norms. For velocity, in $L^{\infty}(L^{2})$ and energy norms, with bounds that
do not depend on the viscosity, and for pressure, in a semi-norm of the same
asymptotic order as the $L^{2}$ norm with respect to the mesh size. The error
estimates between the exact and the ROM-POD solutions are also obtained in the
same norms. These estimates are optimal with respect to the projection error
on the POD reduced spaces. To check the theoretical estimates, a numerical
study is performed on two flow examples, the flow past a cylinder at $Re=200$
and the lid-driven cavity flow at $Re=9500$. The quality of the estimate is
assessed by comparing the logarithmic slopes of the calculated versus the
predicted errors with with respect to the velocity POD contribution ratio. It
is shown on these examples that the error estimates allow a good approximation
of the real errors slopes and thus a good prediction of the their rate of
convergence.
The paper is organised as follows. In Section 2 we introduce the
discretization of the incompressible Navier-Stokes equations (NSE) that we
will consider, while in Section 3 we describe its full order discretization by
stable velocity-pressure finite elements. Section 4 deals with the POD-ROM
including the pressure recovery strategy that we are considering. Section 5 is
devoted to the numerical analysis, stability and error bounds. Numerical
experiments on the flow past a cylinder and the lid-driven cavity flow are
presented in Section 6. Finally conclusions are drawn in Section 7

##  2 Time-dependent NSE: model problem and variational formulation

We introduce an Initial-Boundary Value Problem (IBVP) for the incompressible
evolution NSE. For the sake of simplicity, we just impose the homogeneous
Dirichlet BC on the whole boundary.

Let $[0,T]$ be the time interval and $\Omega$ a bounded polyhedral domain in
$\mathbb{R}^{d}$, $d=2$ or $3$, with a Lipschitz-continuous boundary
$\Gamma=\partial\Omega$. The transient NSE for an incompressible fluid are
given by:

Find ${\bm{u}}:\Omega\times(0,T)\longrightarrow\mathbb{R}^{d}$ and
$p:\Omega\times(0,T)\longrightarrow\mathbb{R}$ such that:

$\left\\{\begin{array}[]{rcll}\partial_{t}{\bm{u}}+({\bm{u}}\cdot\nabla){\bm{u}}-\nu\Delta{\bm{u}}+\nabla
p&=&{\bm{f}}&\,\,\,\mbox{in}\,\,\,\Omega\times(0,T),\\\
\nabla\cdot{\bm{u}}&=&0&\,\,\,\mbox{in}\,\,\,\Omega\times(0,T),\\\
{\bm{u}}&=&\bf{0}&\,\,\,\mbox{on}\,\,\,\Gamma\times(0,T),\\\
{\bm{u}}({\bm{x}},0)&=&{\bm{u}}_{0}({\bm{x}})&\,\,\,\mbox{in}\,\,\,\Omega.\end{array}\right.$
(2.1)

The unknowns are the velocity ${\bm{u}}({\bm{x}},t)$ and the pressure
$p({\bm{x}},t)$ of the incompressible fluid. The data are the source term
${\bm{f}}({\bm{x}},t)$, which represents a body force per mass unit (typically
the gravity), the kinematic viscosity $\nu$ of the fluid, which is a positive
constant, and the initial velocity ${\bm{u}}_{0}({\bm{x}})$.

To define the weak formulation of problem (2.1), we need to introduce some
useful notations for functional spaces [10]. We consider the Sobolev spaces
$H^{s}(\Omega)$, $s\in\mathbb{R}$, $L^{p}(\Omega)$ and $W^{m,p}(\Omega)$,
$m\in\mathbb{N}$, $1\leq p\leq\infty$. We shall use the following notation for
vector-valued Sobolev spaces: ${\bf H}^{s}$, ${\bf L}^{p}$ and ${\bf W}^{m,p}$
respectively shall denote $[H^{s}(\Omega)]^{d}$, $[L^{p}(\Omega)]^{d}$ and
$[W^{m,p}(\Omega)]^{d}$ (similarly for tensor spaces of dimension $d\times
d$). Also, the parabolic Bochner function spaces $L^{p}(0,T;X)$ and
$L^{p}(0,T;{\bf X})$, where $X$ (${\bf X}$) stands for a scalar (vector-
valued) Sobolev space, shall be denoted by $L^{p}(X)$ and $L^{p}({\bf X})$,
respectively. In order to give a variational formulation of problem (2.1), let
us consider the velocity space:

${\bm{X}}={\bf
H}_{0}^{1}=[H_{0}^{1}(\Omega)]^{d}=\left\\{\bm{v}\in[H^{1}(\Omega)]^{d}:\bm{v}={\bf
0}\text{ on }\Gamma\right\\}.$

This is a closed linear subspace of ${\bf H}^{1}$ and thus a Hilbert space
endowed with the ${\bf H}^{1}$-norm. Thanks to Poincaré inequality, the ${\bf
H}^{1}$-norm is equivalent on ${\bf H}_{0}^{1}$ to the norm
$\|{\bm{v}}\|_{{\bf H}_{0}^{1}}=\|{\nabla\bm{v}}\|_{{\bf L}^{2}}$. Also, let
us consider the pressure space:

$Q=L_{0}^{2}(\Omega)=\left\\{q\in
L^{2}(\Omega):\int_{\Omega}q\,d{\bm{x}}=0\right\\}.$

We shall thus consider the following variational formulation of (2.1):

Given ${\bm{f}}\in L^{2}({\bf L}^{2})$, find
${\bm{u}}:(0,T)\longrightarrow{\bm{X}}$, $p:(0,T)\longrightarrow Q$ such that

$\left\\{\begin{array}[]{rcll}\displaystyle\frac{d}{dt}({\bm{u}},\bm{v})+b({\bm{u}},{\bm{u}},\bm{v})+\nu(\nabla{\bm{u}},\nabla\bm{v})-(p,\nabla\cdot\bm{v})&=&({\bm{f}},\bm{v})&\forall\bm{v}\in{\bm{X}},\quad\text{in
}\mathcal{D}^{\prime}(0,T),\\\ (\nabla\cdot{\bm{u}},q)&=&0&\forall q\in
Q,\quad\text{a.e. in }(0,T),\\\ {\bm{u}}(0)&=&{\bm{u}}_{0},\end{array}\right.$
(2.2)

where $(\cdot,\cdot)$ stands for the $L^{2}$-inner product in $\Omega$, and
$\mathcal{D}^{\prime}(0,T)$ is the space of distributions in $(0,T)$. The
trilinear form $b$ is given by: for ${\bm{u}},\,\bm{v},\,{\bm{w}}\in{\bm{X}}$

$b({\bm{u}},\bm{v},{\bm{w}})={1\over
2}\left[({\bm{u}}\cdot\nabla\,\bm{v},{\bm{w}})-({\bm{u}}\cdot\nabla\,{\bm{w}},\bm{v})\right].$
(2.3)

##  3 Finite element full order model

In order to give a FE approximation of (2.2), let $\\{{\cal T}_{h}\\}_{h>0}$
be a family of affine-equivalent, conforming (i.e., without hanging nodes) and
regular triangulations of $\overline{\Omega}$, formed by triangles or
quadrilaterals ($d=2$), tetrahedra or hexahedra ($d=3$). For any mesh cell
$K\in{\cal T}_{h}$, its diameter will be denoted by $h_{K}$ and
$h=\max_{K\in{\cal T}_{h}}h_{K}$. We consider ${\bm{X}}_{h}\subset{\bm{X}}$,
$Q_{h}\subset Q$ being suitable FE spaces for velocity and pressure,
respectively. The FE approximation of (2.2) can be written as follows:

Find $({\bm{u}}_{h},p_{h}):(0,T)\longrightarrow{\bm{X}}_{h}\times Q_{h}$ such
that

$\left\\{\begin{array}[]{rcll}\displaystyle\frac{d}{dt}({\bm{u}}_{h},{\bm{v}}_{h})+b({\bm{u}}_{h},{\bm{u}}_{h},{\bm{v}}_{h})+\nu(\nabla{\bm{u}}_{h},\nabla{\bm{v}}_{h})-(p_{h},\nabla\cdot{\bm{v}}_{h})&=&({\bm{f}},{\bm{v}}_{h})&\text{in
}\mathcal{D}^{\prime}(0,T),\\\ (\nabla\cdot{\bm{u}}_{h},q_{h})&=&0&\text{a.e.
in }(0,T),\\\ {\bm{u}}_{h}(0)&=&{\bm{u}}_{0h},\end{array}\right.$ (3.1)

for any $({\bm{v}}_{h},q_{h})\in{\bm{X}}_{h}\times Q_{h}$, and the initial
condition ${\bm{u}}_{0h}$ is some stable approximation to ${\bm{u}}_{0}$ in
$L^{2}$-norm belonging to ${\bm{X}}_{h}$.

Actually, as a direct method we consider a Galerkin method with grad-div
stabilization, since our aim is to get error bounds with constants independent
on inverse powers of $\nu$. To describe this approach, we define hereafter the
specific choice of FE spaces done both for the numerical analysis and
practical computations in the present work. Given an integer $l\geq 2$ and a
mesh cell $K\in{\cal T}_{h}$, denote by $\mathbb{R}^{l}(K)$ either
$\mathbb{P}^{l}(K)$ (i.e., the space of Lagrange polynomials of degree $\leq
l$, defined on $K$), if the grids are formed by triangles ($d=2$) or
tetrahedra ($d=3$), or $\mathbb{Q}^{l}(K)$ (i.e., the space of Lagrange
polynomials of degree $\leq l$ on each variable, defined on $K$), if the
family of triangulations is formed by quadrilaterals ($d=2$) or hexahedra
($d=3$). We consider the mixed FE pair known as Taylor–Hood elements
$({\bm{X}}_{h}^{l},Q_{h}^{l-1})$ [11, 43] for the velocity-pressure:

$\left\\{\begin{array}[]{lll}&Y_{h}^{l}=V_{h}^{l}(\Omega)=\\{v_{h}\in
C^{0}(\overline{\Omega}):v_{h}|_{K}\in\mathbb{R}^{l}(K),\,\forall K\in{\cal
T}_{h}\\},&\\\ \\\
&{\mathbf{Y}}_{h}^{l}=[Y_{h}^{l}]^{d}=\\{{\bm{v}}_{h}\in[C^{0}(\overline{\Omega})]^{d}:\bm{v}_{h}|_{K}\in[\mathbb{R}^{l}(K)]^{d},\,\forall
K\in{\cal T}_{h}\\},&\\\ \\\
&{\mathbf{X}}_{h}^{l}={\mathbf{Y}}_{h}^{l}\cap{\bf H}_{0}^{1},\quad
Q_{h}^{l-1}=Y_{h}^{l-1}\cap L_{0}^{2}.&\end{array}\right.$ (3.2)

Let us also consider the discrete space of divergence-free functions:

${\bm{V}}_{h}^{l}=\left\\{{\bm{v}}_{h}\in{\mathbf{X}}_{h}^{l}:(\nabla\cdot{\bm{v}}_{h},q_{h})=0\quad\forall
q_{h}\in Q_{h}^{l-1}\right\\}.$

The considered grad-div FOM is given by:

Find $({\bm{u}}_{h},p_{h}):(0,T)\longrightarrow{\mathbf{X}}_{h}^{l}\times
Q_{h}^{l-1}$ such that

$\left\\{\begin{array}[]{rcll}\displaystyle\frac{d}{dt}({\bm{u}}_{h},{\bm{v}}_{h})+b({\bm{u}}_{h},{\bm{u}}_{h},{\bm{v}}_{h})+\nu(\nabla{\bm{u}}_{h},\nabla{\bm{v}}_{h})&-&(p_{h},\nabla\cdot{\bm{v}}_{h})&\\\
+\mu(\nabla\cdot{\bm{u}}_{h},\nabla\cdot{\bm{v}}_{h})&=&({\bm{f}},{\bm{v}}_{h})&\text{in
}\mathcal{D}^{\prime}(0,T),\\\ \\\
(\nabla\cdot{\bm{u}}_{h},q_{h})&=&0&\text{a.e. in }(0,T),\\\
{\bm{u}}_{h}(0)&=&{\bm{u}}_{0h},\end{array}\right.$ (3.3)

for any $({\bm{v}}_{h},q_{h})\in{\mathbf{X}}_{h}^{l}\times Q_{h}^{l-1}$, where
$\mu(\nabla\cdot{\bm{u}}_{h},\nabla\cdot{\bm{v}}_{h})$ denotes the grad-div
stabilization term and $\mu$ is the positive grad-div stabilization parameter.

To state the full space-time discretization of the unsteady grad-div FOM
(3.3), consider a positive integer number $N$ and define $\Delta t=T/N$,
$t_{n}=n\Delta t$, $n=0,1,\ldots,N$. We compute the approximations
${\bm{u}}_{h}^{n}$, $p_{h}^{n}$ to ${\bm{u}}^{n}={\bm{u}}(\cdot,t_{n})$ and
$p^{n}=p(\cdot,t_{n})$ by using, for simplicity of the analysis, a backward
Euler scheme:

  * •

Initialization. Set: ${\bm{u}}_{h}^{0}={\bm{u}}_{0h}.$

  * •

Iteration. For $n=0,1,\ldots,N-1$: Given
${\bm{u}}_{h}^{n}\in{\mathbf{X}}_{h}^{l}$, find
$({\bm{u}}_{h}^{n+1},p_{h}^{n+1})\in{\mathbf{X}}_{h}^{l}\times Q_{h}^{l-1}$
such that:

$\left\\{\begin{array}[]{rcl}\left(\displaystyle\frac{{\bm{u}}_{h}^{n+1}-{\bm{u}}_{h}^{n}}{\Delta
t},{\bm{v}}_{h}\right)+b({\bm{u}}_{h}^{n+1},{\bm{u}}_{h}^{n+1},{\bm{v}}_{h})+\nu(\nabla{\bm{u}}_{h}^{n+1},\nabla{\bm{v}}_{h})&&\\\
-(p_{h}^{n+1},\nabla\cdot{\bm{v}}_{h})+\mu(\nabla\cdot{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\bm{u}}_{h}^{n+1}},\nabla\cdot{\bm{v}}_{h})&=&({\bm{f}}^{n+1},{\bm{v}}_{h}),\\\
\\\ (\nabla\cdot{\bm{u}}_{h}^{n+1},q_{h})&=&0,\end{array}\right.$ (3.4)

for any $({\bm{v}}_{h},q_{h})\in{\mathbf{X}}_{h}^{l}\times Q_{h}^{l-1}$.

It is well-known that considering the discrete divergence-free space
${\bm{V}}_{h}^{l}$ we can remove the pressure from (3.4) since
${\bm{u}}_{h}^{n+1}\in{\bm{V}}_{h}^{l}$ satisfies

$\displaystyle\left(\displaystyle\frac{{\bm{u}}_{h}^{n+1}-{\bm{u}}_{h}^{n}}{\Delta
t},{\bm{v}}_{h}\right)+b({\bm{u}}_{h}^{n+1},{\bm{u}}_{h}^{n+1},{\bm{v}}_{h})+\nu(\nabla{\bm{u}}_{h}^{n+1},\nabla{\bm{v}}_{h})$
$\displaystyle\quad+\mu(\nabla\cdot{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\bm{u}}_{h}^{n+1}},\nabla\cdot{\bm{v}}_{h})=({\bm{f}}^{n+1},{\bm{v}}_{h}),\quad\forall{\bm{v}}_{h}\in{\bm{V}}_{h}^{l}.$
(3.5)

For this method the following bounds hold (see [19, 20, 21]):

$\displaystyle\|{\bm{u}}^{n}-{\bm{u}}_{h}^{n}\|_{{\bf
L}^{2}}+\left(\nu\sum_{j=1}^{n}\Delta
t\|\nabla({\bm{u}}^{j}-{\bm{u}}_{h}^{j})\|_{{\bf L}^{2}}^{2}\right)^{1/2}\leq
C({\bm{u}},p,l+1)(h^{l}+\Delta t),\quad 1\leq n\leq N,$ (3.6)

and

$\displaystyle\left(\sum_{j=1}^{n}\Delta t\|p^{j}-p_{h}^{j}\|_{{\bf
L}^{2}}^{2}\right)^{1/2}\leq C({\bm{u}},p,l+1)h^{-1/2}(h^{l}+\Delta t),\quad
1\leq n\leq N,$ (3.7)

where the constant $C({\bm{u}},p,l+1)$ does not depend on inverse powers of
$\nu$.

##  4 Proper orthogonal decomposition reduced order model

We briefly describe the POD method, following [30], and apply it to the
projection-based stabilized FOM (3.4).

Let us consider the ensembles of velocity snapshots
$\chi^{v}=\text{span}\left\\{{\bm{u}}_{h}^{1},\ldots,{\bm{u}}_{h}^{N}\right\\}$
and pressure snapshots
$\chi^{p}=\text{span}\left\\{p_{h}^{1},\ldots,p_{h}^{N}\right\\}$, given by
the FE solutions to (3.4) at time $t_{n}$, $n=1,\ldots,N$. The POD method
seeks low-dimensional bases
$\left\\{\bm{\varphi}_{1},\ldots,\bm{\varphi}_{r_{v}}\right\\}$ and
$\left\\{\psi_{1},\ldots,\psi_{r_{p}}\right\\}$ in real Hilbert spaces
$\mathcal{H}_{v}$, $\mathcal{H}_{p}$ that optimally approximate the velocity
and pressure snapshots with respect to the discrete $L^{2}(\mathcal{H}_{v})$,
$L^{2}(\mathcal{H}_{p})$ norms, respectively (cf. [30]). It can be shown that
the following POD projection error formulas hold [27, 30]:

$\Delta
t\sum_{n=1}^{N}\left\|{\bm{u}}_{h}^{n}-\sum_{i=1}^{r_{v}}\left({\bm{u}}_{h}^{n},\bm{\varphi}_{i}\right)_{\mathcal{H}_{v}}\bm{\varphi}_{i}\right\|_{\mathcal{H}_{v}}^{2}=\sum_{i=r_{v}+1}^{M_{v}}\lambda_{i},$
(4.1)

and

$\Delta
t\sum_{n=1}^{N}\left\|p_{h}^{n}-\sum_{i=1}^{r_{p}}\left(p_{h}^{n},\psi_{i}\right)_{\mathcal{H}_{p}}\psi_{i}\right\|_{\mathcal{H}_{p}}^{2}=\sum_{i=r_{p}+1}^{M_{p}}\gamma_{i},$
(4.2)

where $M_{v},M_{p}$ are the rank of $\chi^{v}$ and $\chi^{p}$, respectively,
and $\lambda_{i},\gamma_{i}$ are the associated eigenvalues. Although
$\mathcal{H}_{v},\mathcal{H}_{p}$ can be any real Hilbert spaces, in what
follows we consider $\mathcal{H}_{v}={\bf L}^{2}$ and $\mathcal{H}_{p}=L^{2}$.
Also, we are going to take the same number of velocity and pressure POD bases
functions, i.e. $r_{v}=r_{p}=r$ in what follows.

We respectively consider the following velocity and pressure spaces for the
POD setting:

${\bm{X}}_{r}=\text{span}\left\\{\bm{\varphi}_{1},\ldots,\bm{\varphi}_{r}\right\\}\subset{\mathbf{X}}_{h}^{l},$

and

$Q_{r}=\text{span}\left\\{\psi_{1},\ldots,\psi_{r}\right\\}\subset
Q_{h}^{l-1}.$

######  Remark 4.1.

Since the POD velocity modes are linear combinations of the snapshots, they
satisfy the boundary conditions in (2.1) and are solenoidal. Thus, the POD
velocity modes belong to ${\bm{V}}_{h}^{l}$, which yields
${\bm{X}}_{r}\subset{\bm{V}}_{h}^{l}$.

The standard Galerkin projection-based POD-ROM uses both Galerkin truncation
and Galerkin projection. The former yields an approximation of the velocity
and pressure fields by a linear combination of the corresponding truncated POD
basis:

${\bm{u}}({\bm{x}},t)\approx{\bm{u}}_{r}({\bm{x}},t)=\sum_{i=1}^{r}a_{i}(t)\bm{\varphi}_{i}({\bm{x}}),$
(4.3)

and

$p({\bm{x}},t)\approx
p_{r}({\bm{x}},t)=\sum_{i=1}^{r}b_{i}(t)\psi_{i}({\bm{x}}),$ (4.4)

where $\left\\{a_{i}(t)\right\\}_{i=1}^{r}$ and
$\left\\{b_{i}(t)\right\\}_{i=1}^{r}$ are the sought time-varying coefficients
representing the POD-Galerkin velocity and pressure trajectories. Note that
$r<<{\cal N}$, where ${\cal N}$ denotes the number of degrees of freedom
(d.o.f.) in a full order simulation. Replacing the velocity-pressure FE pair
$({\bm{u}}_{h},p_{h})$ with $({\bm{u}}_{r},p_{r})$ in the FE approximation
(3.4) and projecting the resulted equations onto the POD product space
$({\bm{X}}_{r},Q_{r})$ using the POD basis
$\left(\left\\{\bm{\varphi}_{i}\right\\}_{i=1}^{r},\left\\{\psi_{i}\right\\}_{i=1}^{r}\right)$,
the full space-time discretization of the grad-div ROM reads as:

  * •

Initialization. Set:
${\bm{u}}_{r}^{0}=\displaystyle\sum_{i=1}^{r}({\bm{u}}_{0},\bm{\varphi}_{i})\bm{\varphi}_{i}.$

  * •

Iteration. For $n=0,1,\ldots,N-1$: Given ${\bm{u}}_{r}^{n}\in{\bm{X}}_{r}$,
find ${\bm{u}}_{r}^{n+1}\in{\bm{X}}_{r}$ such that:

$\displaystyle\left(\displaystyle\frac{{\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}}{\Delta
t},\bm{\varphi}\right)+b({\bm{u}}_{r}^{n+1},{\bm{u}}_{r}^{n+1},\bm{\varphi})+\nu(\nabla{\bm{u}}_{r}^{n+1},\nabla\bm{\varphi})+\mu(\nabla\cdot\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\bm{u}}_{r}^{n+1}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},\nabla\cdot\bm{v}_{r})$
$\displaystyle\quad=({\bm{f}}^{n+1},\bm{\varphi}),\quad\forall\bm{\varphi}\in{\bm{X}}_{r}.$
(4.5)

In (• ‣ 4), the pressure term vanishes due to the fact that ${\bm{u}}_{r}^{n}$
belongs to the discrete divergence-free space ${\bm{V}}_{h}^{l}$. In order to
recover the reduced order pressure, we propose to solve in a post-process
phase the residual-based equation, as follows: Find $p_{r}^{n+1}\in Q_{r}$
such that

$\displaystyle\sum_{K\in{\cal T}_{h}}\tau_{K}(\nabla
p_{r}^{n+1},\nabla\psi)_{K}$ (4.6) $\displaystyle\quad=-\sum_{K\in{\cal
T}_{h}}\tau_{K}(\frac{{\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}}{\Delta
t}+\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1}-\nu\,\Delta{\bm{u}}_{r}^{n+1}-{\bm{f}}^{n+1},\nabla\psi)_{K},$

for all $\psi\in Q_{r}$. Here, the $\tau_{K}$ are positive coefficients of
order $h_{K}^{2}$, that is, there exist two constants $C_{1}>0$, $C_{2}>0$,
such that

$C_{1}\,h_{K}^{2}\leq\tau_{K}\leq C_{2}\,h_{K}^{2},\quad\forall K\in{\cal
T}_{h},\,\,\forall h>0.$ (4.7)

Also, the symbol $\tilde{{\bm{u}}}_{r}$ stands for ${{\bm{u}}}_{r}$ when
$d=2$, and for some truncation of ${{\bm{u}}}_{r}$ such that

$\|\tilde{{\bm{u}}}_{r}\|_{L^{\infty}(K)}\leq C_{T}\,h_{K}^{-1},\quad\forall
K\in{\cal T}_{h},\,\,\forall h>0$ (4.8)

when $d=3$, for some constant $C_{T}>0$. We shall denote in this way all
generic constants that will appear along the paper, unless needed to specify
them by some index. Such a $\tilde{{\bm{u}}}_{r}$ may be, for instance,
defined as follows:

$\tilde{{\bm{u}}}_{r}=\Pi_{h}(T_{\alpha}({{\bm{u}}}_{r})),\quad\mbox{with
}\alpha=h_{K}^{-1}$

where $\Pi_{h}$ denotes the standard Lagrange interpolation operator on space
${\bm{V}}_{h}^{l}$, $T_{\alpha}$ is the truncation at height $\alpha>0$, and
$\alpha$ is the function defined on $\bar{\Omega}$ that assigns to any
$x\in\bar{\Omega}$ the inverse of the average of the grid elements that
contain $x$. It holds

$T_{\alpha}({\bm{u}})(x_{i})=\left\\{\begin{array}[]{ccl}{\bm{u}}(x_{i})&\mbox{if}&|{\bm{u}}(x_{i})|\leq\alpha(x_{i}),\\\
\alpha(x_{i})\,sign({\bm{u}}(x_{i}))&\mbox{if}&|{\bm{u}}(x_{i})|\geq\alpha(x_{i}).\end{array}\right.$

for any Lagrange interpolation node $x_{i}$. Then, as the grids are regular,
$|T_{\alpha}({\bm{u}})(x)|\leq Ch_{K}^{-1}$ for any $x\in\bar{\Omega}$, where
$K\in{\cal T}_{h}$ is the element containing $x$, for some constant $C>0$.

Due to inverse finite element error estimates and the local $L^{2}$ stability
of the Lagrange interpolation (see [17], Appendix),

$\|\tilde{{\bm{u}}}_{r}\|_{L^{\infty}(K)}\leq
C\,h_{K}^{-3/2}\|\tilde{{\bm{u}}}_{r}\|_{L^{2}(K)}\leq
C\,h_{K}^{-3/2}\,\|T_{\alpha}({{\bm{u}}}_{r})\|_{L^{2}(K)}\leq
C\,h_{K}^{-3/2}\,|K|^{1/2}\,\alpha\leq C\,h_{K}^{-1}.$

Note that the use of inverse finite element estimates combined with the
$L^{\infty}(L^{2})$ bound for the ${\bm{u}}_{r}$ yields

$\|{{\bm{u}}}_{r}\|_{L^{\infty}(K)}\leq
C\,h_{K}^{-3/2}\|{{\bm{u}}}_{r}\|_{L^{2}(K)}\leq C\,h_{K}^{-3/2}.$

Then, the truncation at height $h_{K}^{-1}$, although somewhat restrictive, is
rather mild, in particular it allows the blow-up of the $L^{\infty}(\Omega)$
norm of $\tilde{{\bm{u}}}_{r}$ as $h\rightarrow 0$ . Anyhow, this technical
hypothesis may be overcome when $d=2$ or a weaker norm of the pressure
gradient is bounded when $d=3$.

The structure of method (4.6) is inspired in that of stabilized methods for
fluid flows. Usually a weighted residual (the difference between the l.h.s.
and the r.h.s. of (4.6)) is added to the Galerkin formulation to stabilize
both convection dominance and pressure gradient [1]. The scaling of the
stabilizing coefficients $\tau_{K}$ as $h_{K}^{2}$ comes from dimensional
analysis. Here, we isolate this stabilizing term to recover the pressure, once
the velocity is known.

The presence of the factors $\tau_{K}$ in (4.6) only is needed whenever the
grids ${\cal T}_{h}$ are not uniformly regular. Otherwise all the $\tau_{K}$
may be taken equal to $h^{2}$ and (4.6) may be replaced by

$\displaystyle(\nabla
p_{r}^{n+1},\nabla\psi)_{\Omega}=-(\frac{{\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}}{\Delta
t}+\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1}-\nu\,\Delta{\bm{u}}_{r}^{n+1}-{\bm{f}}^{n+1},\nabla\psi)_{\Omega}.$
(4.9)

Note that equations (4.6) and (4.9) contain a discretization of the natural b.
c. for pressure, that is (assuming
$\tilde{{\bm{u}}}_{r}^{n+1}={\bm{u}}_{r}^{n+1}$)

$\nabla p_{r}^{n+1}\cdot{\bf
n}=\left(\frac{{\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}}{\Delta
t}+{\bm{u}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1}-\nu\,\Delta{\bm{u}}_{r}^{n+1}-{\bm{f}}^{n+1}\right)\cdot{\bf
n}\,\,\mbox{on }\partial\Omega.$

Thus, no artificial boundary conditions on the pressure are being imposed.

Hereafter, the grad-div ROM (• ‣ 4) together with (4.6) will be referred to as
Stabilization Motivated-ROM (SM-ROM). The SM-ROM (• ‣ 4)-(4.6), without grad-
div term, has been introduced and numerically investigated in [15].

Let us remark that an alternative time discretization could be given by the
semi-implicit Euler method, where the trilinear form in (3.4) is discretized
by $b({\bm{u}}_{h}^{n+1},{\bm{u}}_{h}^{n},{\bm{v}}_{h})$. This semi-implicit
time discretization is less costly from the computational point of view with
respect to a fully implicit one, which yields a nonlinear algebraic system of
equations to be solved. However, the numerical analysis will be performed in
Section for the more technical case of the fully implicit time discretization
given by (3.4). We shall present in Section 6 numerical results obtained with
the semi-implicit discretization.

##  5 Analysis of the residual-based SM-ROM

In this section, we perform the numerical analysis of the proposed unsteady
SM-ROM (• ‣ 4)-(4.6), dealing with stability and error estimates.

###  5.1 Technical background

This section provides some technical results that are required for the
numerical analysis. Throughout the paper, we shall denote by $C$ a positive
constant that may vary from a line to another, but which is always independent
of the viscosity $\nu$, the FE mesh size $h$, the FE velocity interpolation
order $l$, the time step $\Delta t$, and the velocity, pressure eigenvalues
$\lambda_{i}$, $\gamma_{i}$.

We shall denote by $\|\cdot\|_{m,p}$ the norm in $W^{m,p}(\Omega)$ and by
$\|\cdot\|_{m,p,K}$ the norm in $W^{m,p}(K)$ for some $K\in{\cal T}_{h}$.
Also, we shall denote by $\|\cdot\|_{m}$ the norm in $H^{m}(\Omega)$ and by
$\|\cdot\|_{m,K}$ the norm in $H^{m}(K)$ for some $K\in{\cal T}_{h}$.

We shall use some inverse estimates for finite element functions (cf. [18],
Theorem 3.2.6): there exists two constants $C_{I,1}>0$ and $C_{I,2}>0$ such
that

$\|v_{h}\|_{0,\infty,K}\leq
C_{I,1}\,h_{K}^{-\frac{d}{2}}\,\|v_{h}\|_{0,K},\,\,\|\Delta v_{h}\|_{0,K}\leq
C_{I,2}\,h_{K}^{-1}\,\|\nabla v_{h}\|_{0,K},\,\,\mbox{for all }v_{h}\in
V_{h}^{l}.$ (5.1)

We define the scalar product:

$(\cdot,\cdot)_{\tau}:{\bf L}^{2}(\Omega)\times{\bf
L}^{2}(\Omega)\to\mathbb{R},\quad(\bm{v},{\bm{w}})_{\tau}=\sum_{K\in{\cal
T}_{h}}\tau_{K}(\bm{v},{\bm{w}})_{K},$

and its associated norm:

$\|{\bm{v}}\|_{\tau}=(\bm{v},\bm{v})_{\tau}^{1/2}.$

We shall also use the following discrete Gronwall lemma:

######  Lemma 5.1.

Let $\alpha_{n}$, $\beta_{n}$, $\gamma_{n}$, $n=1,2,...$ be non-negative real
numbers satisfying

$(1-\sigma\,\Delta t)\,\alpha_{n+1}+\beta_{n+1}\leq(1+\tau\,\Delta
t)\,\alpha_{n}+\gamma_{n+1}$ (5.2)

for some $\sigma\geq 0$, $\tau\geq 0$. Assume that $\sigma\,\Delta t\leq
1-\delta$ for some $\delta>0$. Then it holds

$\alpha_{n}\leq
e^{\rho\,t_{n}}\,\alpha_{0}+\frac{1}{\delta}\,\sum_{l=1}^{n}e^{\rho\,(t_{n}-t_{l})}\,\gamma_{l};$
(5.3)
$\sum_{l=1}^{n}\beta_{l}\leq\left(1+\frac{\tau}{\sigma}+(\sigma+\tau)\,e^{\rho\,t_{n-1}}\,t_{n-1}\right)\,\alpha_{0}+\frac{1}{\delta}\,\left(1+(\sigma+\tau)\,e^{\rho\,t_{n-1}}\,t_{n-1}\right)\sum_{l=1}^{n}\,\gamma_{l},$
(5.4)

with $\displaystyle\rho=\frac{\sigma+\tau}{\delta}$.

###### Proof.

Dividing by $1-\sigma\,\Delta t$ in (5.2),

$\displaystyle\alpha_{n+1}+\beta_{n+1}$ $\displaystyle\leq$
$\displaystyle\alpha_{n+1}+\frac{1}{1-\sigma\Delta
t}\beta_{n+1}\leq(1+\frac{(\sigma+\tau)\Delta t}{1-\sigma\Delta
t})\alpha_{n}+\frac{1}{1-\sigma\Delta t}\,\gamma_{n+1}$ (5.5)
$\displaystyle\leq$ $\displaystyle(1+\rho\Delta
t)\alpha_{n}+\frac{1}{\delta}\,\gamma_{n+1}.$

Recursively using this inequality,

$\displaystyle\alpha_{n+1}+\beta_{n+1}$ $\displaystyle\leq$
$\displaystyle(1+\rho\Delta t)^{2}\,\alpha_{n-1}+(1+\rho\Delta
t)\,\,\frac{1}{\delta}\,\gamma_{n}+\frac{1}{\delta}\,\gamma_{n+1}\leq\cdots$
(5.6) $\displaystyle\leq$ $\displaystyle(1+\rho\Delta
t)^{n+1}\,\alpha_{0}+\frac{1}{\delta}\,\sum_{l=1}^{n+1}(1+\rho\Delta
t)^{n+1-l}\,\gamma_{l}$ $\displaystyle\leq$ $\displaystyle e^{\rho
t_{n+1}}\,\alpha_{0}+\frac{1}{\delta}\,\sum_{l=1}^{n+1}e^{\rho(t_{n+1}-t_{l})}\,\gamma_{l},$

hence (5.3). Summing now the inequalities (5.2) yields

$(1-\sigma\,\Delta
t)\,\sum_{l=1}^{n+1}\alpha_{i}+\sum_{l=1}^{n+1}\beta_{l}\leq(1+\tau\,\Delta
t)\,\sum_{l=0}^{n}\alpha_{i}+\frac{1}{\delta}\,\sum_{l=1}^{n+1}\gamma_{l}.$

Then,

$\displaystyle\sum_{l=1}^{n+1}\beta_{l}$ $\displaystyle\leq$
$\displaystyle(1+\tau\,\Delta t)\,\alpha_{0}+(\sigma+\tau)\Delta
t\,\sum_{l=1}^{n}\alpha_{l}+\frac{1}{\delta}\,\sum_{l=1}^{n+1}\gamma_{l}$
(5.7) $\displaystyle\leq$ $\displaystyle(1+\tau\,\Delta
t)\,\alpha_{0}+(\sigma+\tau)\,\Delta t\,\sum_{l=1}^{n}\left[e^{\rho
t_{l}}\,\alpha_{0}+\frac{1}{\delta}\,\sum_{r=1}^{l}e^{\rho(t_{l}-t_{r})}\,\gamma_{l}\right]+\frac{1}{\delta}\,\sum_{l=1}^{n+1}\gamma_{l}$
$\displaystyle\leq$
$\displaystyle\left(1+\frac{\tau}{\sigma}+(\sigma+\tau)\,e^{\rho
t_{n}}\,t_{n}\right)\,\alpha_{0}+\frac{1}{\delta}\,\left((\sigma+\tau)\,e^{\rho\,t_{n}}\,t_{n}+1\right)\,\sum_{l=1}^{n+1}\gamma_{l},$

hence (5.4). $\Box$

Note that the constant $\rho$ plays the role of an amplification factor for
the exponential dependence on the length of the time interval of the bounds
appearing in (5.3) and (5.4). It tends to 1 as $\Delta t$ tends to zero.

###  5.2 Existence and stability results for SM-ROM

We have the following existence and unconditional stability result for the SM-
ROM,(• ‣ 4)-(4.6):

######  Theorem 5.2.

Assume that $\Delta t\leq 1-\delta$ for some $\delta>0$, and that ${\bm{f}}\in
L^{2}(0,T;{\bf L}^{2}(\Omega))$. Assume also that the data ${\bm{f}}^{n}$ in
problems (• ‣ 4) and (4.6) satisfy

$\|{\bm{u}}^{0}_{r}\|_{0}\leq\|{\bm{u}}_{0}\|_{0},\quad\Delta
t\,\sum_{i=1}^{n}\|{\bm{f}}^{n}\|_{0}^{2}\leq\|{\bm{f}}\|^{2}_{L^{2}(0,t_{n};{\bf
L}^{2}(\Omega))},\quad\mbox{for }n=1,\cdots,N.$

Then the following statements hold.

  1. 1.

Equations (• ‣ 4) admits at least a solution ${\bm{u}}_{r}^{k+1}$, and (4.6)
admits a unique solution $p_{r}^{n+1}$.

  2. 2.

Any solution ${\bm{u}}_{r}^{k+1}$ of equation (• ‣ 4) satisfies the following
estimates

$\displaystyle\|{\bm{u}}_{r}^{k}\|_{0}^{2}\leq
e^{\rho\,t_{k}}\,\|{\bm{u}}_{0}\|_{0}^{2}+\frac{1}{\delta}\,e^{\rho\,t_{k}}\,\|{\bf
f}\|_{L^{2}(0,t_{k};{\bf L}^{2}(\Omega)^{2})}^{2},\quad k=1,\cdots,N;$ (5.8)
$\displaystyle\Delta
t\,\sum_{n=1}^{k}\left(\nu\,\|\nabla{\bm{u}}_{r}^{n}\|_{0}^{2}+\mu\,\|\nabla\cdot{\bm{u}}_{r}^{n}\|_{0}^{2}\right)\leq\frac{1}{2}\,\left(1+t_{k}\,e^{\rho\,t_{k}}\right)\|{\bm{u}}_{0}\|_{0}^{2}$
(5.9)
$\displaystyle+\frac{1}{2\delta}\,\left(1+t_{k}\,e^{\rho\,t_{k}}\right)\,\|{\bf
f}\|_{L^{2}(0,t_{k};{\bf L}^{2}(\Omega))}^{2},\quad k=1,\cdots,N$

with $\displaystyle\rho=\frac{1}{\delta}$.

  3. 3.

Assume in addition that $h\leq C_{S}\,\Delta t$ for some constant $C_{S}>0$,
$\Delta t\leq\displaystyle\frac{1-\delta}{s+1/5}$ with $s=5\,C_{2}\,C_{S}^{2}$
for some $0<\rho<1$, and that in problem (4.6) ,

     * 3a)

When $d=2$, it holds

$C_{2}\leq\displaystyle\sqrt{\frac{\nu}{10}}\,\frac{1}{A\,C_{I,1}},\,\mbox{
with
}A=e^{\rho\,T}\,\|{\bm{u}}_{0}\|_{0}^{2}+\frac{1}{\delta}\,e^{\rho\,T}\,\|{\bf
f}\|_{L^{2}(0,T;{\bf L}^{2}(\Omega)^{2})}^{2},$ (5.10)

or

     * 3b)

When $d=3$, it holds $\|\tilde{{\bm{u}}}_{r}^{n}\|_{0,\infty,K}\leq
C_{T}\,h_{K}^{-1}$ with

$C_{T}\leq\nu\,\sqrt{C_{I,2}},\quad
C_{2}\leq\displaystyle\frac{1}{10}\,\frac{1}{\nu\,C_{I,2}}.$ (5.11)

Then the solution $p_{r}^{1},\cdots,p_{r}^{N}$ of equation (4.6) satisfies the
following estimates:

$\displaystyle\Delta t\,\sum_{n=1}^{N}\|\nabla p_{r}^{n}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle\left[2+\left(2s+1/5\right)\,e^{\rho
t}\,T\right]\|{\bm{u}}_{0}\|_{0}^{2}$
$\displaystyle+\frac{5}{\delta}\,\left[1+\left(2s+1/5\right)\,e^{\rho
t}\,T\right]\,(1+C_{2}\,h^{2})\,\|{\bm{f}}\|_{L^{2}(0,T;{\bf
L}^{2}(\Omega))}.$

###### Proof.

1.- Once ${\bm{u}}_{r}^{n}$ and ${\bm{u}}_{r}^{n+1}$ are known, problem (4.6)
is a standard elliptic problem with a form $a_{pres}(p,q)=(\nabla p,\nabla
q)_{\Omega}$ which is coercive on sub-spaces of $H^{1}(\Omega)$ with vanishing
mean in $\Omega$. Then, problem (4.6) admits a unique solution.

The existence of solution of equation (• ‣ 4) follows from a standard finite-
dimensional compactness argument, that we omit for brevity.

2.- Estimate (5.9) follows from a standard argument for finite element
approximation of incompressible flow problems using the discrete Gronwall
Lemma 5.1 that can be found, for instance, in [17], Chapter 10. We omit it for
brevity

3.- Let us take $\bm{v}_{r}={\bm{u}}_{r}^{n+1}$ in problem (• ‣ 4). Using the
identities

$({\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n},{\bm{u}}_{r}^{n+1})_{\Omega}=\frac{1}{2}\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}-\frac{1}{2}\,\|{\bm{u}}_{r}^{n}\|_{0}^{2}+\frac{1}{2}\,\|{\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}\|_{0}^{2},$

and

$b({\bm{u}}_{r}^{n};{\bm{u}}_{r}^{n+1},{\bm{u}}_{r}^{n+1})=0$

it holds

$\displaystyle\frac{1}{2}\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}-\frac{1}{2}\,\|{\bm{u}}_{r}^{n}\|_{0}^{2}+\frac{1}{2}\,\|{\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}\|_{0}^{2}+\nu\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2}+2\,\mu\,\Delta
t\,\|\nabla\cdot{\bm{u}}_{r}^{n+1}\|_{0}^{2}$ $\displaystyle=\Delta
t\,({\bm{f}}^{n+1},{\bm{u}}_{r}^{n+1})_{\Omega}\leq\frac{\Delta
t}{2}\,\|{\bm{f}}^{n+1}\|_{0}^{2}+\frac{\Delta
t}{2}\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}.$ (5.13)

Taking $\psi=p_{r}^{n+1}$ and multiplying by $\Delta t$ in problem (4.6), and
summing up the first identity in (5.13) yields

$\displaystyle\frac{1}{2}\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}-\frac{1}{2}\,\|{\bm{u}}_{r}^{n}\|_{0}^{2}+\nu\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\Delta t\,\|\nabla
p_{r}^{n+1}\|_{\tau}^{2}$ (5.14) $\displaystyle\leq-\sum_{K\in{\cal
T}_{h}}\tau_{K}\,({\bm{u}}_{r}^{n+1}-{\bm{u}}_{r}^{n}+\Delta
t\,(\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1}-\nu\,\Delta{\bm{u}}_{r}^{n+1}-{\bm{f}}^{n+1}),\nabla
p_{r}^{n+1})_{K}$ $\displaystyle+\Delta
t\,({\bm{f}}^{n+1},{\bm{u}}_{r}^{n+1})_{\Omega}$

We bound each term on the r.h.s. of this expression. Observe that

$\sum_{K\in{\cal
T}_{h}}\tau_{K}\,(a,b)_{K}\leq\|a\|_{\tau}\,\|b\|_{\tau}\leq\frac{1}{2}\beta^{-1}\,\|a\|_{\tau}^{2}+\frac{1}{2}\beta\,\|b\|_{\tau}^{2}\quad$

for any functions $a,\,b\in L^{2}(\Omega)$ or $a,\,b\in{\bf L}^{2}(\Omega)$
and any $\beta>0$. Using $\tau_{K}\leq C_{2}\,h_{K}^{2}$, $h_{K}\leq
C_{S}\,\Delta t$ and estimate (5.8) it follows

$\displaystyle\sum_{K\in{\cal T}_{h}}\tau_{K}\,({\bm{u}}_{r}^{n+1},\nabla
p_{r}^{n+1})_{K}\leq\frac{1}{2\varepsilon\Delta
t}\,\|{\bm{u}}_{r}^{n+1}\|_{\tau}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}$
$\displaystyle\leq\frac{1}{2\varepsilon}\,\Delta
t\,C_{2}\,C_{S}^{2}\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2};$ (5.15)

where $\varepsilon>0$ is a constant to be determined later. Similarly,

$\displaystyle\sum_{n=0}^{N-1}\sum_{K\in{\cal
T}_{h}}\tau_{K}\,({\bm{u}}_{r}^{n},\nabla
p_{r}^{n+1})_{K}\leq\frac{1}{2\varepsilon}\,C_{2}\,C_{S}^{2}\,\|{\bm{u}}_{r}^{n}\|_{0}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\sum_{n=1}^{N}\|\nabla p_{r}^{n}\|_{\tau}^{2};$ (5.16)

Further,

$\displaystyle-\Delta t\,\sum_{K\in{\cal
T}_{h}}\tau_{K}\,(\nu\,\Delta{\bm{u}}_{r}^{n+1},\nabla
p_{r}^{n+1})_{K}\leq\frac{\nu^{2}}{2\varepsilon}\,\Delta
t\,\|\Delta{\bm{u}}_{r}^{n+1}\|_{\tau}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}$
$\displaystyle\leq\frac{\nu^{2}}{2\varepsilon}\,C_{I,2}\,\max_{K\in{\cal
T}_{h}}\tau_{K}\,h_{K}^{-2}\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}$
$\displaystyle\leq\frac{\nu^{2}}{2\varepsilon}\,C_{I,2}\,C_{2}\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}.$ (5.17)

The convection-pressure gradient interaction term will be specifically bounded
in each of the situations 3a) and 3b), as follows.

3a) In this case $d=2$ and $\tilde{{\bm{u}}}_{r}^{n+1}={\bm{u}}_{r}^{n+1}$. It
holds

$\displaystyle\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,(\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1},\nabla
p_{r}^{n+1})_{K}\leq\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,\|{\bm{u}}_{r}^{n+1}\|_{0,\infty,K}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}\,\|\nabla
p_{r}^{n+1}\|_{0,K}$ $\displaystyle\displaystyle\leq C_{I,1}\,\,\Delta
t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,h_{K}^{-d/2}\,\|{\bm{u}}_{r}^{n+1}\|_{0,K}\,\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}\,\|\nabla
p_{r}^{n+1}\|_{0,K}$ $\displaystyle\displaystyle\leq A\,C_{I,1}\,\,\Delta
t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,h_{K}^{-d/2}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}\,\|\nabla
p_{r}^{n+1}\|_{0,K},$
$\displaystyle\displaystyle\leq\frac{1}{2\varepsilon}\,A^{2}\,C_{I,1}^{2}\,\Delta
t\,\sum_{K\in{\cal
T}_{h}}\tau_{K}\,h_{K}^{-d}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}$
$\displaystyle\displaystyle\leq\frac{1}{2\varepsilon}\,A^{2}\,C_{I,1}^{2}\,C_{2}\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}.$ (5.18)

where $A$ is defined in (5.10).

3b) In this case $d=3$ and $\|\tilde{{\bm{u}}}_{r}^{n+1}\|_{0,\infty,K}\leq
C_{T}\,h_{K}^{-1}$. It holds

$\displaystyle\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,(\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1},\nabla
p_{r}^{n+1})_{K}\leq\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,\|\tilde{{\bm{u}}}_{r}^{n+1}\|_{0,\infty,K}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}\,\|\nabla
p_{r}^{n+1}\|_{0,K}$ $\displaystyle\displaystyle\leq C_{T}\,\Delta
t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,h_{K}^{-1}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}\,\|\nabla
p_{r}^{n+1}\|_{0,K}$
$\displaystyle\displaystyle\leq\frac{1}{2\varepsilon}\,C_{T}^{2}\,\Delta
t\,\sum_{K\in{\cal
T}_{h}}\tau_{K}\,h_{K}^{-2}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0,K}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}$
$\displaystyle\displaystyle\leq\frac{1}{2\varepsilon}\,C_{T}^{2}\,C_{2}\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}.$ (5.19)

Further,

$\displaystyle\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\,({\bm{f}}^{n+1},\nabla
p_{r}^{n+1})_{K}\leq\displaystyle\frac{1}{2\varepsilon}\,\Delta
t\,\|{\bm{f}}^{n+1}\|_{\tau}^{2}+\frac{\varepsilon}{2}\,\Delta t\,\|\nabla
p_{r}^{n+1}\|_{\tau}^{2};$ (5.20)

and

$\displaystyle\displaystyle\Delta
t\,({\bm{f}}^{n+1},{\bm{u}}_{r}^{n+1})\leq\displaystyle\frac{1}{2\varepsilon}\,\Delta
t\,\|{\bm{f}}^{n+1}\|_{0}^{2}+\frac{\varepsilon}{2}\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2};$ (5.21)

Taking now $\varepsilon=1/5$ and summing up (5.15), (5.16), (5.17), (5.20) and
(5.21) with (5.18) (case 3a)), or (5.19) (case 3b)) yields

$\displaystyle(1-(s+1/5)\Delta t)\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\mu\,\Delta
t\,\|\nabla{\bm{u}}_{r}^{n+1}\|_{0}^{2}+\Delta t\,\|\nabla
p_{r}^{n+1}\|_{\tau}^{2}$ $\displaystyle\leq(1+s\Delta
t)\,\|{\bm{u}}_{r}^{n}\|_{0}^{2}+5\,\Delta
t\,\left(\|{\bm{f}}^{n+1}\|_{0}^{2}+\|{\bm{f}}^{n+1}\|_{\tau}^{2}\right),$
$\displaystyle\leq(1+s\Delta t)\,\|{\bm{u}}_{r}^{n}\|_{0}^{2}+5\,\Delta
t\,(1+C_{2}\,h^{2})\,\|{\bm{f}}^{n+1}\|_{0}^{2}$ (5.22)

where

$\mu=2\,\left(\nu-\frac{1}{2\varepsilon}\,A^{2}\,C_{I,1}^{2}\,C_{2}-\frac{\nu^{2}}{2\varepsilon}\,C_{2}\,C_{I,2}\right)\,\quad\mbox{when
}d=2,$

or

$\mu=2\,\left(\nu-\frac{1}{2\varepsilon}\,C_{2}\,C_{T}^{2}-\frac{\nu^{2}}{2\varepsilon}\,C_{2}\,C_{I,2}\right)\,\quad\mbox{when
}d=3.$

In any case, in view of conditions (5.10) and (5.11), it holds $\mu\geq\nu$.
Applying Lemma 5.1 to (5.22) with $\sigma=s+1/5$, $\tau=s$,
$\alpha_{n}=\mu\,\Delta t\,\|{\bm{u}}_{r}^{n+1}\|_{0}^{2}$, $\beta_{n}=\Delta
t\,\|\nabla p_{r}^{n+1}\|_{\tau}^{2}$ and $\gamma_{n}=5\,\Delta
t\,(1+h^{2})\,\|{\bm{f}}^{n+1}\|_{0}^{2}$ yields estimate (3) . $\Box$

######  Remark 5.3.

Estimates for velocity (5.8), (5.9) and for pressure (3) do not degenerate as
$\nu\to 0$. This is obtained using the enhanced regularity ${\bm{f}}\in
L^{2}(0,T;{\bf L}^{2}(\Omega))$. Also, due to the conditions satisfied by
constant $C_{2}$ in either (5.10) or (5.11), (3) provides either an estimate
for the norm

$\nu\,\Delta t\,\sum_{i=1}^{N-1}\|p_{r}^{n+1}\|^{2}_{h}\,\,\,\mbox{when
}d=2,\,\,\mbox{or }\Delta
t\,\sum_{i=1}^{N-1}\|p_{r}^{n+1}\|^{2}_{h},\,\mbox{when }d=3,\,\,\mbox{with }$
$\|p\|_{h}=\left(\sum_{i=1}^{N-1}h_{K}^{2}\,\|\nabla
p\|^{2}_{0,K}\right)^{1/2}.$

The $\|\cdot\|_{h}$ norm plays a role in the discrete inf-sup condition for
non-stable pairs of finite element velocity-pressure spaces (such as
${\mathbf{X}}_{h}^{l}$ and $Q_{h}^{l}$, for instance, with equal interpolation
in velocity and pressure, cf. [16]):

$\|p_{h}\|_{0}\leq\sup_{{\bm{v}}_{h}\in{\mathbf{X}}_{h}^{l}}\frac{(\nabla\cdot{\bm{v}}_{h},p_{h})}{\|{\bm{v}}_{h}\|_{1,\Omega}}+C\,\|p_{h}\|_{h}\,\mbox{for
any }\,\,p_{h}\in Q_{h}^{l}.$

Also, a straightforward argument using inverse inequalities proves that

$\|p_{h}\|_{h}\leq C\,\|p_{h}\|_{0}\,\,\mbox{for any }\,\,p_{h}\in
Q_{h}^{l-1}.$

###  5.3 Error estimates for SM-ROM

We are now in position to state the following error estimate results for the
SM-ROM defined by (• ‣ 4)-(4.6):

######  Theorem 5.1 (Velocity error estimate).

Let ${\bm{u}}$ be the velocity in the NSE (2.1), let ${\bm{u}}_{r}$ be the
grad-div ROM velocity defined in (• ‣ 4), and assume that the solution
$({\bm{u}},p)$ of (2.1) is regular enough. Then, the following bound holds:

$\Delta t\sum_{n=0}^{N-1}\|{\bm{u}}_{r}^{n+1}-{\bm{u}}^{n+1}\|_{0}^{2}\leq
C\left(\sum_{i=r+1}^{M_{v}}\lambda_{i}\left(1+\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}\right)+h^{2l}+\Delta
t^{2}\right).$ (5.23)

The proof of this result follows along the same lines as the proof of Theorem
5.3 in [34], we omit it for brevity.

######  Theorem 5.2 (Pressure error estimate).

Let $p$ be the pressure in the NSE (2.1), let $p_{r}$ be the SM-ROM pressure
defined in (4.6), and assume that the solution $({\bm{u}},p)$ of (2.1) is
regular enough. Then, the following bound holds:

$\Delta t\sum_{n=0}^{N-1}\|\nabla(p_{r}^{n+1}-p^{n+1})\|_{\tau}^{2}\leq
CE^{*},$ (5.24)

where:

$\displaystyle E^{*}$ $\displaystyle=$ $\displaystyle{h^{2}}+(h^{2l}+\Delta
t^{2})(1+\nu^{-1}+\|{S_{r}^{v}}\|_{2}+h^{-1})$ (5.25)
$\displaystyle+\sum_{i=r+1}^{M_{v}}\lambda_{i}\left[1+\left(1+\|{S_{r}^{v}}\|_{2}\right)\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}\right]+h^{2}\sum_{i=r+1}^{M_{p}}\gamma_{i}\|{\nabla\psi_{i}}\|_{0}^{2},$

and the penalty term $h^{2}$ in (5.25) actually vanishes when $d=2$.

Proof. We start deriving the pressure error bound by splitting the error for
the velocity and the pressure into two terms:

${\bm{u}}^{n+1}-{\bm{u}}_{r}^{n+1}=({\bm{u}}^{n+1}-P_{r}^{v}{\bm{u}}^{n+1})-({\bm{u}}_{r}^{n+1}-P_{r}^{v}{\bm{u}}^{n+1})=\bm{\eta}^{n+1}-\bm{\phi}_{r}^{n+1},$
(5.26)
$p^{n+1}-p_{r}^{n+1}=(p^{n+1}-P_{r}^{p}p^{n+1})-(p_{r}^{n+1}-P_{r}^{p}p^{n+1})=\rho^{n+1}-s_{r}^{n+1}.$
(5.27)

In (5.26), the first term,
$\bm{\eta}^{n+1}={\bm{u}}^{n+1}-P_{r}^{v}{\bm{u}}^{n+1}$, represents the
difference between ${\bm{u}}^{n+1}$ and its $L^{2}$-orthogonal projection on
${\bm{X}}_{r}$. The second term,
$\bm{\phi}_{r}^{n+1}={\bm{u}}_{r}^{n+1}-P_{r}^{v}{\bm{u}}^{n+1}$, is the
remainder. Similarly, in (5.27), the first term,
$\rho^{n+1}=p^{n+1}-P_{r}^{p}p^{n+1}$, represents the difference between
$p^{n+1}$ and its $L^{2}$-orthogonal projection on $Q_{r}$. The second term,
$s_{r}^{n+1}=p_{r}^{n+1}-P_{r}^{p}p^{n+1}$, is the remainder.

Next, we construct the error equation. We first consider the solution
$({\bm{u}},p)$ of (2.1) at $t=t_{n+1}$ and:

$\sum_{K\in{\cal
T}_{h}}\tau_{K}(\partial_{t}{\bm{u}}^{n+1}+{\bm{u}}^{n+1}\cdot\nabla{\bm{u}}^{n+1}-\nu\,\Delta{\bm{u}}^{n+1}+\nabla
p^{n+1}-{\bm{f}}^{n+1},\nabla s_{r}^{n+1})_{K}=0,$ (5.28)

then subtract (4.6) with $\psi=s_{r}^{n+1}$ from it. By adding and subtracting
the difference quotient term $\displaystyle\frac{1}{\Delta
t}({\bm{u}}^{n+1}-{\bm{u}}^{n})$, and applying the decompositions
(5.26)-(5.27), we get:

$\begin{array}[]{lll}&&\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(\partial_{t}{\bm{u}}^{n+1}-\displaystyle\frac{{\bm{u}}^{n+1}-{\bm{u}}^{n}}{\Delta
t}+\displaystyle\frac{1}{\Delta
t}(\bm{\eta}^{n+1}-\bm{\phi}_{r}^{n+1})-\displaystyle\frac{1}{\Delta
t}(\bm{\eta}^{n}-\bm{\phi}_{r}^{n}),\nabla s_{r}^{n+1}\right)_{K}\\\ \\\
&+&\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(({\bm{u}}^{n+1}\cdot\nabla{\bm{u}}^{n+1}-\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1})-\nu\Delta(\bm{\eta}^{n+1}-\bm{\phi}_{r}^{n+1}),\nabla
s_{r}^{n+1}\right)_{K}\\\ \\\ &+&\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(\nabla(\rho^{n+1}-s_{r}^{n+1}),\nabla
s_{r}^{n+1}\right)_{K}=0.\end{array}$ (5.29)

Letting
$\bm{c}^{n}=\partial_{t}{\bm{u}}^{n+1}-\displaystyle\frac{{\bm{u}}^{n+1}-{\bm{u}}^{n}}{\Delta
t}$, from (5.29) we obtain:

$\begin{array}[]{lll}&&\Delta t\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\|{\nabla s_{r}^{n+1}}\|_{0,K}^{2}=\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(\Delta
t\bm{c}^{n}+(\bm{\eta}^{n+1}-\bm{\eta}^{n})-(\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n}),\nabla
s_{r}^{n+1}\right)_{K}\\\ \\\ &+&\Delta t\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(({\bm{u}}^{n+1}\cdot\nabla{\bm{u}}^{n+1}-\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1})-\nu\Delta(\bm{\eta}^{n+1}-\bm{\phi}_{r}^{n+1}),\nabla
s_{r}^{n+1}\right)_{K}\\\ \\\ &+&\Delta t\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(\nabla\rho^{n+1},\nabla
s_{r}^{n+1}\right)_{K}=I+II+III+IV+V+VI.\end{array}$ (5.30)

We estimate the terms on the r.h.s. of (5.30) one by one. By Cauchy-Schwarz
and Young’ s inequalities, we bound the first term on the r.h.s. of (5.30):

$I=\Delta t\sum_{K\in{\cal T}_{h}}\tau_{K}(\bm{c}^{n},\nabla
s_{r}^{n+1})_{K}\leq\Delta t\sum_{K\in{\cal
T}_{h}}\frac{\tau_{K}}{2\varepsilon}\|{\bm{c}^{n}}\|_{0,K}^{2}+\Delta
t\sum_{K\in{\cal T}_{h}}\tau_{K}\frac{\varepsilon}{2}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2},$ (5.31)

for some small positive constant $\varepsilon$ (to be determined later).

For the second term on the r.h.s. of (5.30), assuming $h\sim\Delta t$, we
have:

$\displaystyle II$ $\displaystyle=$ $\displaystyle-\sum_{K\in{\cal
T}_{h}}\tau_{K}(\bm{\eta}^{n+1}-\bm{\eta}^{n},\nabla s_{r}^{n+1})_{K}$ (5.32)
$\displaystyle\leq$ $\displaystyle\sum_{K\in{\cal
T}_{h}}\frac{\tau_{K}}{2\varepsilon\Delta
t}\|{\bm{\eta}^{n+1}-\bm{\eta}^{n}}\|_{0,K}^{2}+\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\frac{\varepsilon}{2}\|{\nabla s_{r}^{n+1}}\|_{0,K}^{2}$
$\displaystyle\leq$ $\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\|{\bm{\eta}^{n+1}-\bm{\eta}^{n}}\|_{0,K}^{2}+\Delta
t\sum_{K\in{\cal T}_{h}}\tau_{K}\frac{\varepsilon}{2}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}.$

Note that this term actually vanishes for uniformly regular grids (take
$\tau_{K}=h^{2}$ and integrate by parts).

Similarly, for the third term on the r.h.s. of (5.30), we have:

$\displaystyle III$ $\displaystyle=$ $\displaystyle\sum_{K\in{\cal
T}_{h}}\tau_{K}(\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n},\nabla s_{r}^{n+1})_{K}$
(5.33) $\displaystyle\leq$ $\displaystyle\sum_{K\in{\cal
T}_{h}}\frac{\tau_{K}}{2\varepsilon\Delta
t}\|{\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n}}\|_{0,K}^{2}+\Delta
t\sum_{K\in{\cal T}_{h}}\tau_{K}\frac{\varepsilon}{2}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}$ $\displaystyle\leq$ $\displaystyle\Delta
t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\|{\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n}}\|_{0,K}^{2}+\Delta
t\sum_{K\in{\cal T}_{h}}\tau_{K}\frac{\varepsilon}{2}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}.$

The fourth term on the r.h.s. of (5.30) can be bounded as follows:

$\displaystyle IV$ $\displaystyle=$ $\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}({\bm{u}}^{n+1}\cdot\nabla{\bm{u}}^{n+1}-\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla{\bm{u}}_{r}^{n+1},\nabla
s_{r}^{n+1})_{K}$ (5.34) $\displaystyle=$ $\displaystyle\Delta
t\sum_{K\in{\cal
T}_{h}}\tau_{K}(\bm{\eta}^{n+1}\cdot\nabla{\bm{u}}^{n+1}-\bm{\phi}_{r}^{n+1}\cdot\nabla{\bm{u}}^{n+1},\nabla
s_{r}^{n+1})_{K}$ $\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}(\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla\bm{\eta}^{n+1}-\tilde{{\bm{u}}}_{r}^{n+1}\cdot\nabla\bm{\phi}_{r},\nabla
s_{r}^{n+1})_{K}$ $\displaystyle+{\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(({\bm{u}}_{r}^{n+1}-\tilde{{\bm{u}}}_{r}^{n+1})\cdot\nabla{\bm{u}}^{n+1},\nabla
s_{r}^{n+1}\right)_{K}}$ $\displaystyle\leq$ $\displaystyle\Delta
t\sum_{K\in{\cal
T}_{h}}\tau_{K}\|{\bm{\eta}^{n+1}}\|_{0,K}\|{\nabla{\bm{u}}^{n+1}}\|_{0,\infty,K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}$ $\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\|{\bm{\phi}_{r}^{n+1}}\|_{0,K}\|{\nabla{\bm{u}}^{n+1}}\|_{0,\infty,K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}$ $\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\|{\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,\infty,K}\|{\nabla\bm{\eta}^{n+1}}\|_{0,K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}$ $\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\|{\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,\infty,K}\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0,K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}$ $\displaystyle+{\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\|{{\bm{u}}_{r}^{n+1}-\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,K}\|{\nabla{\bm{u}}^{n+1}}\|_{0,\infty,K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}}$ $\displaystyle\leq$ $\displaystyle\Delta
t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\tau_{K}(\|{\bm{\eta}^{n+1}}\|_{0,K}^{2}+\|{\bm{\phi}_{r}^{n+1}}\|_{0,K}^{2})+5\Delta
t\sum_{K\in{\cal T}_{h}}\frac{\varepsilon}{2}\tau_{K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}$ $\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\frac{\tau_{K}}{h_{K}^{2}}\left(\|{\nabla\bm{\eta}^{n+1}}\|_{0,K}^{2}+\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0,K}^{2}\right)$
$\displaystyle+{\Delta t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\tau_{K}(\|{{\bm{u}}_{r}^{n+1}}\|_{0,K}^{2}+\|{\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,K}^{2})}$
$\displaystyle\leq$ $\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\tau_{K}(\|{\bm{\eta}^{n+1}}\|_{0,K}^{2}+\|{\bm{\phi}_{r}^{n+1}}\|_{0,K}^{2})+5\Delta
t\sum_{K\in{\cal T}_{h}}\frac{\varepsilon}{2}\tau_{K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}$ $\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\frac{C}{\varepsilon}\frac{\tau_{K}}{h_{K}^{2}}\left(\|{\nabla\bm{\eta}^{n+1}}\|_{0,K}^{2}+\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0,K}^{2}\right)+{\Delta
t\frac{C}{\varepsilon}h^{2}},$

where we used Hölder’s, Young’s and triangle inequalities, the velocity
stability estimate $\eqref{estvel1}$, and (4.8) to bound
$\tilde{{\bm{u}}}_{r}^{n+1}$ when $d=3$. Note that the last penalty term in
(5.34) vanishes when $d=2$, since
$\tilde{{\bm{u}}}_{r}^{n+1}={\bm{u}}_{r}^{n+1}$.

Using local inverse estimates [9], we can bound the fifth term on the r.h.s.
of (5.30) as:

$\displaystyle V$ $\displaystyle=$ $\displaystyle-\nu\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\left(\Delta(\bm{\eta}^{n+1}-\bm{\phi}_{r}^{n+1}),\nabla
s_{r}^{n+1}\right)_{K}$ (5.35) $\displaystyle\leq$ $\displaystyle\nu^{2}\Delta
t\sum_{K\in{\cal
T}_{h}}\frac{\tau_{K}}{2\varepsilon}\left(\|{\Delta\bm{\eta}^{n+1}}\|_{0,K}^{2}+\|{\Delta\bm{\phi}_{r}^{n+1}}\|_{0,K}^{2}\right)+2\Delta
t\sum_{K\in{\cal T}_{h}}\frac{\varepsilon}{2}\tau_{K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}$ $\displaystyle\leq$
$\displaystyle\frac{C}{\varepsilon}\nu^{2}\,\Delta
t(\|{\Delta\bm{\eta}^{n+1}}\|_{\tau}^{2}+\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0}^{2})+2\Delta
t\sum_{K\in{\cal T}_{h}}\frac{\varepsilon}{2}\tau_{K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2},$

Finally, by Cauchy-Schwarz and Young’s inequalities, we bound the sixth term
on the r.h.s. of (5.30):

$\displaystyle VI=\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}(\nabla\rho^{n+1},\nabla s_{r}^{n+1})_{K}$ $\displaystyle\leq$
$\displaystyle\Delta t\sum_{K\in{\cal
T}_{h}}\frac{\tau_{K}}{2\varepsilon}\|{\nabla\rho^{n+1}}\|_{0,K}^{2}$ (5.36)
$\displaystyle+\Delta t\sum_{K\in{\cal
T}_{h}}\tau_{K}\frac{\varepsilon}{2}\|{\nabla s_{r}^{n+1}}\|_{0,K}^{2}.$

Summarizing (5.31)-(5.36) and taking $\varepsilon=1/11$ from (5.30), we get:

$\displaystyle\Delta t\sum_{K\in{\cal T}_{h}}\tau_{K}\|{\nabla
s_{r}^{n+1}}\|_{0,K}^{2}$ $\displaystyle\leq$ $\displaystyle C\Delta
t\left(h^{2}\|{\bm{c}^{n}}\|_{0}^{2}+\|{\bm{\eta}^{n+1}-\bm{\eta}^{n}}\|_{0}^{2}+\|{\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n}}\|_{0}^{2}\right)$
(5.37) $\displaystyle+C\Delta
th^{2}(\|{\bm{\eta}^{n+1}}\|_{0}^{2}+\|{\bm{\phi}_{r}^{n+1}}\|_{0}^{2}+{1})$
$\displaystyle+C\Delta
t(\|{\nabla\bm{\eta}^{n+1}}\|_{0}^{2}+\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0}^{2})$
$\displaystyle+C\nu^{2}\Delta t\|{\Delta\bm{\eta}^{n+1}}\|_{\tau}^{2}+C\Delta
t\|{\nabla\rho^{n+1}}\|_{\tau}^{2}.$

Summing (5.37) from $n=0$ to $N-1$, we have:

$\displaystyle\Delta t\sum_{n=0}^{N-1}\|{\nabla s_{r}^{n+1}}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle C\Delta
t\sum_{n=0}^{N-1}\left(h^{2}\|{\bm{c}^{n}}\|_{0}^{2}+\|{\bm{\eta}^{n+1}-\bm{\eta}^{n}}\|_{0}^{2}+\|{\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n}}\|_{0}^{2}\right)$
(5.38) $\displaystyle+C\Delta
t\sum_{n=0}^{N-1}h^{2}(\|{\bm{\eta}^{n+1}}\|_{0}^{2}+\|{\bm{\phi}_{r}^{n+1}}\|_{0}^{2})+{Ch^{2}}$
$\displaystyle+C\Delta
t\sum_{n=0}^{N-1}(\|{\nabla\bm{\eta}^{n+1}}\|_{0}^{2}+\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0}^{2})$
$\displaystyle+C\Delta
t\sum_{n=0}^{N-1}\left(\nu^{2}\|{\Delta\bm{\eta}^{n+1}}\|_{\tau}^{2}+\|{\nabla\rho^{n+1}}\|_{\tau}^{2}\right).$

Next, we estimate each term on the r.h.s. of (5.38).

By Taylor’s theorem, the first term on the r.h.s. of (5.38) can be estimated
as follows:

$\Delta t\sum_{n=0}^{N-1}h^{2}\|{\bm{c}^{n}}\|_{{\bf L}^{2}}^{2}\leq C\Delta
t^{2}h^{2}\|{\partial_{tt}{\bm{u}}}\|_{L^{2}({\bf L}^{2})}^{2}\leq C\Delta
t^{2}h^{2}.$ (5.39)

To estimate the second term on the r.h.s. of (5.38), we use triangle
inequality and the $L^{2}$ error estimate for $P_{r}^{v}{\bm{u}}$ (see [28],
Lemma 3.3):

$\Delta t\sum_{n=0}^{N-1}\|{\bm{\eta}^{n+1}-\bm{\eta}^{n}}\|_{0}^{2}\leq
C\left(h^{2l}+\Delta t^{2}+\sum_{i=r+1}^{M_{v}}\lambda_{i}\right).$ (5.40)

By using triangle and Cauchy-Schwarz inequalities, the bound (3.6), and
estimate (73) in [34], the third term on the r.h.s. of (5.38) can be estimated
as follows:

$\displaystyle\Delta
t\sum_{n=0}^{N-1}\|{\bm{\phi}_{r}^{n+1}-\bm{\phi}_{r}^{n}}\|_{0}^{2}$
$\displaystyle\leq$ $\displaystyle 2\Delta
t\sum_{n=0}^{N-1}\left(\|{({\bm{u}}_{r}^{n+1}-P_{r}^{v}{\bm{u}}_{h}^{n+1})-({\bm{u}}_{r}^{n}-P_{r}^{v}{\bm{u}}_{h}^{n})}\|_{0}^{2}\right)$
(5.41) $\displaystyle+2\Delta
t\sum_{n=0}^{N-1}\left(\|{P_{r}^{v}\left[({\bm{u}}_{h}^{n+1}-{\bm{u}}^{n+1})-({\bm{u}}_{h}^{n}-{\bm{u}}^{n})\right]}\|_{0}^{2}\right)$
$\displaystyle\leq$ $\displaystyle
C\sum_{i=r+1}^{M_{v}}\lambda_{i}\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}$
$\displaystyle+C\Delta
t\sum_{n=0}^{N-1}\left(\|{({\bm{u}}_{h}^{n+1}-{\bm{u}}^{n+1})-({\bm{u}}_{h}^{n}-{\bm{u}}^{n})}\|_{0}^{2}\right)$
$\displaystyle\leq$ $\displaystyle
C\left(\sum_{i=r+1}^{M_{v}}\lambda_{i}\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}+h^{2l}+\Delta
t^{2}\right).$

Similarly, for the fourth and fifth term on the r.h.s. of (5.38) we have:

$\Delta t\sum_{n=0}^{N-1}h^{2}\|{\bm{\eta}^{n+1}}\|_{0}^{2}\leq
C\,h^{2}\left(h^{2l}+\Delta t^{2}+\sum_{i=r+1}^{M_{v}}\lambda_{i}\right),$
(5.42)

and

$\displaystyle\Delta t\sum_{n=0}^{N-1}h^{2}\|{\bm{\phi}_{r}^{n+1}}\|_{0}^{2}$
$\displaystyle\leq$ $\displaystyle 2\Delta
t\sum_{n=0}^{N-1}h^{2}\left(\|{{\bm{u}}_{r}^{n+1}-P_{r}^{v}{\bm{u}}_{h}^{n+1}}\|_{0}^{2}+\|{P_{r}^{v}({\bm{u}}_{h}^{n+1}-{\bm{u}}^{n+1})}\|_{0}^{2}\right)$
(5.43) $\displaystyle\leq$ $\displaystyle
C\,h^{2}\left(\sum_{i=r+1}^{M_{v}}\lambda_{i}\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}+h^{2l}+\Delta
t^{2}\right).$

For the sixth term on the r.h.s. of (5.38), we have:

$\Delta t\sum_{n=0}^{N-1}\|{\nabla\bm{\eta}^{n+1}}\|_{0}^{2}\leq
C\left[(h^{2l}+\Delta
t^{2})(\nu^{-1}+\|{S_{r}^{v}}\|_{2})+\sum_{i=r+1}^{M_{v}}\lambda_{i}\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}\right],$
(5.44)

where we used the $H_{0}^{1}$ error estimate for $P_{r}^{v}{\bm{u}}$ (see
[28], Lemma 3.3), and $\|{S_{r}^{v}}\|_{2}$ denotes the $2$-norm of the
stiffness velocity matrix with entries
$[S_{r}^{v}]_{ij}=(\nabla\bm{\varphi}_{i},\nabla\bm{\varphi}_{j})$,
$i,j=1,\ldots,r$ that comes from the use of POD inverse estimates (see [30],
Lemma 2).

For the seventh term on the r.h.s. of (5.38), we have:

$\displaystyle\Delta t\sum_{n=0}^{N-1}\|{\nabla\bm{\phi}_{r}^{n+1}}\|_{0}^{2}$
$\displaystyle\leq$ $\displaystyle 2\Delta
t\sum_{n=0}^{N-1}\left(\|{\nabla({\bm{u}}_{r}^{n+1}-P_{r}^{v}{\bm{u}}_{h}^{n+1})}\|_{0}^{2}+\|{\nabla\left[P_{r}^{v}({\bm{u}}_{h}^{n+1}-{\bm{u}}^{n+1})\right]}\|_{0}^{2}\right)$
(5.45) $\displaystyle\leq$ $\displaystyle
C\|{S_{r}^{v}}\|_{2}\left(\sum_{i=r+1}^{M_{v}}\lambda_{i}\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}+h^{2l}+\Delta
t^{2}\right),$

where we used the triangle inequality, POD inverse estimates (see [30], Lemma
2), estimate (73) in [34], and bound (3.6).

The eight term on the r.h.s. of (5.38) can be bounded as follows:

$\displaystyle\Delta
t\sum_{n=0}^{N-1}\nu^{2}\|{\Delta\bm{\eta}^{n+1}}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle C\nu^{2}\Delta
t\sum_{n=0}^{N-1}\left[\|{\Delta({\bm{u}}^{n+1}-{\bm{u}}_{h}^{n+1})}\|_{\tau}^{2}+\|{\Delta({\bm{u}}_{h}^{n+1}-P_{r}^{v}{\bm{u}}_{h}^{n+1})}\|_{\tau}^{2}\right]$
(5.46) $\displaystyle+C\nu^{2}\Delta
t\sum_{n=0}^{N-1}\|{\Delta[P_{r}^{v}({\bm{u}}_{h}^{n+1}-{\bm{u}}^{n+1})]}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle C\nu^{2}\Delta
t\sum_{n=0}^{N-1}\left(h^{2l}+\|{\nabla({\bm{u}}^{n+1}-{\bm{u}}_{h}^{n+1})}\|_{0}^{2}+\|{\nabla({\bm{u}}_{h}^{n+1}-P_{r}^{v}{\bm{u}}_{h}^{n+1})}\|_{0}^{2}\right)$
$\displaystyle+C\nu^{2}\|{S_{r}^{v}}\|_{2}\,\Delta
t\sum_{n=0}^{N-1}\|{{\bm{u}}^{n+1}-{\bm{u}}_{h}^{n+1}}\|_{0}^{2}$
$\displaystyle\leq$ $\displaystyle
C\left[\nu^{2}\,h^{2l}+\nu(1+\nu\|{S_{r}^{v}}\|_{2})(h^{2l}+\Delta
t^{2})+\nu^{2}\sum_{i=r+1}^{M_{v}}\lambda_{i}\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}\right],$

where we used the triangle inequality, optimal and local inverse estimates
[9], POD inverse estimates (see [30], Lemma 2), bound (3.6), and the POD
projection error in $H^{1}$-seminorm for ${\bm{u}}_{h}$ (see [23], Lemma 3.2).

Finally, for the last term on the r.h.s. of (5.38), similarly to (5.44) we
obtain:

$\displaystyle\Delta t\sum_{n=0}^{N-1}\|{\nabla\rho^{n+1}}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle C\Delta
t\sum_{n=0}^{N-1}\left[\|{\nabla(p^{n+1}-p_{h}^{n+1})}\|_{\tau}^{2}+\|{\nabla(p_{h}^{n+1}-P_{r}^{p}p_{h}^{n+1})}\|_{\tau}^{2}\right]$
(5.47) $\displaystyle+C\Delta
t\sum_{n=0}^{N-1}\|{\nabla[P_{r}^{p}(p_{h}^{n+1}-p^{n+1})]}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle C\Delta
t\sum_{n=0}^{N-1}\left(h^{2l}+\|{p^{n+1}-p_{h}^{n+1}}\|_{0}^{2}+h^{2}\|{\nabla(p_{h}^{n+1}-P_{r}^{p}p_{h}^{n+1})}\|_{0}^{2}\right)$
$\displaystyle\leq$ $\displaystyle C\left[h^{2l}+h^{-1}(h^{2l}+\Delta
t^{2})+h^{2}\sum_{i=r+1}^{M_{p}}\gamma_{i}\|{\nabla\psi_{i}}\|_{0}^{2}\right],$

where we used the triangle inequality, optimal and local inverse estimates
[9], bound (3.7), and the POD projection error in $H^{1}$-seminorm for $p_{h}$
(see [23], Lemma 3.2).

Collecting (5.39)-(5.47), estimate (5.38) becomes:

$\displaystyle\Delta t\sum_{n=0}^{N-1}\|{\nabla s_{r}^{n+1}}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle C\left[{h^{2}}+(h^{2l}+\Delta
t^{2})(1+\nu^{-1}+\|{S_{r}^{v}}\|_{2}+h^{-1})\right]$ (5.48)
$\displaystyle+C\sum_{i=r+1}^{M_{v}}\lambda_{i}\left[1+\left(1+\|{S_{r}^{v}}\|_{2}\right)\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}\right]$
$\displaystyle+C\,h^{2}\sum_{i=r+1}^{M_{p}}\gamma_{i}\|{\nabla\psi_{i}}\|_{0}^{2}.$

Finally, using the triangle inequality and estimate (5.47), we conclude:

$\displaystyle\Delta
t\sum_{n=0}^{N-1}\|{\nabla(p_{r}^{n+1}-p^{n+1})}\|_{\tau}^{2}$
$\displaystyle\leq$ $\displaystyle 2\Delta
t\sum_{n=0}^{N-1}\left(\|{\nabla\rho^{n+1}}\|_{\tau}^{2}+\|{\nabla
s_{r}^{n+1}}\|_{\tau}^{2}\right)$ (5.49) $\displaystyle\leq$ $\displaystyle
C\left[{h^{2}}+(h^{2l}+\Delta
t^{2})(1+\nu^{-1}+\|{S_{r}^{v}}\|_{2}+h^{-1})\right]$
$\displaystyle+C\sum_{i=r+1}^{M_{v}}\lambda_{i}\left[1+\left(1+\|{S_{r}^{v}}\|_{2}\right)\|{\nabla\bm{\varphi}_{i}}\|_{0}^{2}\right]$
$\displaystyle+C\,h^{2}\sum_{i=r+1}^{M_{p}}\gamma_{i}\|{\nabla\psi_{i}}\|_{0}^{2},$

where we recall that the penalty term $h^{2}$ in (5.49) vanishes when $d=2$.
This concludes the proof. $\Box$

######  Remark 5.4.

An alternative pressure error estimate can be obtained by using a different
bound for the term $IV$ in (5.34) when $d=3$, taking advantage of the fact
that $\tilde{{\bm{u}}}_{r}^{n+1}={{\bm{u}}}_{r}^{n+1}$ when
$\|{{{\bm{u}}}_{r}^{n+1}}\|_{0,\infty,K}\leq Ch_{K}^{-1}$. Indeed, we can
replace in (5.34) the last term by:

$\displaystyle{\Delta t\sum_{K\in{\cal
T}_{h}^{*}}\frac{C}{\varepsilon}\tau_{K}(\|{{\bm{u}}_{r}^{n+1}}\|_{0,K}^{2}+\|{\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,K}^{2})}$
$\displaystyle\leq$ $\displaystyle{\Delta t\sum_{K\in{\cal
T}_{h}^{*}}\frac{C}{\varepsilon}\tau_{K}|K|(\|{{\bm{u}}_{r}^{n+1}}\|_{0,\infty,K}^{2}+\|{\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,\infty,K}^{2})}$
(5.50) $\displaystyle\leq$ $\displaystyle{\Delta t\sum_{K\in{\cal
T}_{h}^{*}}\frac{C}{\varepsilon}h_{K}^{2}|K|\left(h_{K}^{-3}\|{{\bm{u}}_{r}^{n+1}}\|_{0,K}^{2}+h_{K}^{-2}\right)}$
$\displaystyle\leq$ $\displaystyle{\Delta
t\,\frac{C}{\varepsilon}\,h^{-1}|\Omega^{*}|},$

where ${\cal T}_{h}^{*}=\\{K\in{\cal
T}_{h}:\|{{\bm{u}}_{r}^{n+1}}\|_{0,\infty,K}\geq Ch_{K}^{-1}\\}$,
$|\Omega^{*}|$ is the measure of ${\cal T}_{h}^{*}$ (i.e.,
$|\Omega^{*}|=\displaystyle\sum_{K\in{\cal T}_{h}^{*}}|K|$), and we used local
inverse estimates [9], the velocity stability estimate $\eqref{estvel1}$, and
(4.8) to bound $\tilde{{\bm{u}}}_{r}^{n+1}$, assuming also uniformly regular
grids.

Now, we have:

$\displaystyle|\Omega^{*}|$ $\displaystyle=$ $\displaystyle\sum_{K\in{\cal
T}_{h}^{*}}\int_{K}\frac{h_{K}}{h_{K}}\,d{\bm{x}}\leq C\sum_{K\in{\cal
T}_{h}^{*}}\int_{K}h_{K}\|{{\bm{u}}_{r}^{n+1}}\|_{0,\infty,K}d{\bm{x}}$ (5.51)
$\displaystyle\leq$ $\displaystyle C\sum_{K\in{\cal
T}_{h}^{*}}\int_{K}h_{K}^{1/2}\|{{\bm{u}}_{r}^{n+1}}\|_{0,6,K}d{\bm{x}}\leq
C\sum_{K\in{\cal
T}_{h}^{*}}\|{\nabla{\bm{u}}_{r}^{n+1}}\|_{0,K}\,|K|\,h_{K}^{1/2}$
$\displaystyle\leq$ $\displaystyle C\left(\sum_{K\in{\cal
T}_{h}^{*}}\|{\nabla{\bm{u}}_{r}^{n+1}}\|_{0,K}^{2}\right)^{1/2}\left(\sum_{K\in{\cal
T}_{h}^{*}}|K|^{2}h_{K}\right)^{1/2}$ $\displaystyle\leq$ $\displaystyle
C\|{\nabla{\bm{u}}_{r}^{n+1}}\|_{0,\Omega}\left(\sum_{K\in{\cal
T}_{h}^{*}}|K|h_{K}^{4}\right)^{1/2}\leq
C\|{\nabla{\bm{u}}_{r}^{n+1}}\|_{0,\Omega}h^{2}|\Omega^{*}|^{1/2}.$

where we used the definition of ${\cal T}_{h}^{*}$, local inverse estimates
[9], Hölder’s inequality, the Sobolev embedding $H^{1}\hookrightarrow L^{6}$,
and Cauchy-Schwarz inequality. Thus, from (5.51) we get:

$|\Omega^{*}|\leq C\|{\nabla{\bm{u}}_{r}^{n+1}}\|_{0,\Omega}^{2}\,h^{4},$
(5.52)

and substituting (5.52) in (5.50) we obtain:

${\Delta t\sum_{K\in{\cal
T}_{h}^{*}}\frac{C}{\varepsilon}\tau_{K}(\|{{\bm{u}}_{r}^{n+1}}\|_{0,K}^{2}+\|{\tilde{{\bm{u}}}_{r}^{n+1}}\|_{0,K}^{2})}\leq{\Delta
t\,\frac{C}{\varepsilon}\,h^{3}\|{\nabla{\bm{u}}_{r}^{n+1}}\|_{0,\Omega}^{2}}.$
(5.53)

Thus, by applying the velocity stability estimate $\eqref{estvel2}$, in this
case we can replace the penalty term $h^{2}$ in (5.25) by $h^{3}\nu^{-1}$.

##  6 Numerical tests

The data used for the construction of the POD subspaces are obtained by using
Fenics [31] with Taylor-Hood element $\mathbb{P}_{2}/\mathbb{P}_{1}$. To
quantify the error estimates derived in the previous section, we consider two
classical benchmark numerical examples, the flow past a cylinder and the lid-
driven cavity flow. Since the errors depend strongly on the velocity POD
truncation order, the error analyses are presented in terms of the velocity
POD contribution ratio given as

$\mathcal{R}_{u}=1-\displaystyle{\overset{r}{\underset{i=1}{\sum}}}\lambda_{i}/\displaystyle{\overset{\mathcal{N}}{\underset{i=1}{\sum}}}\lambda_{i}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}$
(6.1)

where $r$ and $\mathcal{N}$ are respectively the truncation order and the
total number of snapshots. This contribution ratio can be seen as a refinement
parameter of the error. It can be used along with the predicted error-slope to
determine a priori the accuracy of the SM-ROM.

###  6.1 Flow past a cylinder Re = 200

Consider the two dimensional flow past a cylinder of diameter $D$, placed in a
rectangular domain of length $H=30D$ and width = $45D$. The center of the
cylinder is situated at $L_{1}=10D$ from the inlet and $H/2$ from horizontal
walls. The flow is driven by an inlet velocity $U$ and is allowed to flow
through the outlet. Free slip boundary conditions are applied at the
horizontal edges while no slip boundary condition are considered on the
cylinder’s wall. Illustration of boundary condition is given in Figure 1.

For the numerical simulation, the semi-implicit Euler method is used with time
step $\Delta t=10^{-2}$, while a non-uniform triangular mesh made of $21174$
cells is used for the spatial discretization. At $Re=200$, the flow creates
alternating low-high pressure vortices downstream the cylinder, triggering the
generation of the Von Karman vortex pattern in the wake of the flow past the
cylinder. These structures are illustrated in Figure 2.
In order to construct the velocity and pressure POD bases, $400$ uniformly
distributed snapshots covering $8$ periods of the periodic regime of the flow
are considered. The evolution of the contribution ratio of POD modes is given
as a function of the modes numbers in Figure 3. A first inspection of the
contribution ratio reveals that the first three modes carry the biggest share
of the flow energy. After that, the contribution evolves by couple of modes.
This behavior is essentially due the double multiplicity of POD eigenvalues
revealing the presence of eigenspaces of dimension two. In that case, it’s
worth noting that the SM-ROM error is most likely to oscillate if the POD
basis is truncated at a mode forming the first vector of a two dimensional POD
eigenspace. To prevent that from happening, it is preferred to construct the
SM-ROM gradually with respect to eigenspaces rather than to eigenvectors.

The POD contribution ratio in Figure 3 tends rapidly to the numerical zero
meaning by that the few first modes are sufficiently enough to catch the
quasi-totality of the flow energy. The prediction quality of the SM-ROM
improves with decreasing values of $\mathcal{R}_{u}$. This can be seen from
the hydrodynamics coefficients $C_{D}$ and $C_{L}$ given by

$C_{D}=\displaystyle{\frac{F_{D}}{\frac{1}{2}\rho U^{2}D}}\ \ ,\ \
C_{L}=\displaystyle{\frac{F_{L}}{\frac{1}{2}\rho
U^{2}D}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}$

where $\rho$ and $U$ are respectively the density of the fluid and the
characteristic velocity, and $F_{D}$ and $F_{L}$ are the drag and lift forces
exerted by the fluid on the cylinder’s boundary

$\left(\begin{matrix}F_{D}\\\
F_{L}\end{matrix}\right)=\int_{{{}_{\textnormal{cylinder}}}}\left(\mu(\nabla{\bm{u}}+\nabla^{T}{\bm{u}})-pI\right)\cdot{\bm{n}}\,d\sigma{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}$

Here $\mu$ refers to the dynamic viscosity. Figure 4 gives a confrontation of
the hydrodynamic coefficients for variable magnitudes of the ratio
$\mathcal{R}_{u}$. It shows that the predicted coefficients by SM-ROM converge
towards the reference coefficients as $\mathcal{R}_{u}$ drops down to zero.
Precisely, when the ratio attains an order of $10^{-4}$ or less, the solutions
are well recovered and the SM-ROM coefficients are in good agreement with the
underlying high fidelity model.
In the following, we investigate the errors evolution with respect to the POD
contribution ratio. In theory, the error by the SM-ROM and the POD
contribution ratio must decrease simultaneously. However, this is not always
the case. It can be noticed from Figure 5 that while the SM-ROM velocity
errors decrease as expected, the SM-ROM pressure error experiences a
stagnation for ratios smaller than $10^{-6}$, meaning that the SM-ROM has
attained its maximum precision at this stage. This can be explained as the
high fidelity weak formulation (3.4) is not exactly verified for gradients of
pressure POD modes, because these do not lie in the discrete velocity space
$\bm{X}_{h}$. This effect does not appear in the error estimates (5.24), as
these apply to the error between POD-ROM and exact pressures. This limit
accuracy for the reduced pressure computation should decrease as the
discretization parameters $h$ and $\Delta t$ decrease.
Reduced order model errors with respect to the rate of convergence are
reported in Figure 5. For the sake of clarity, and since we are interested
only in the error slopes, the constant in the error estimation is adjusted in
such a way that the estimated and SM-ROM errors start from the same point.
Moreover, a power law regression with solid lines is overlayed on the computed
errors and a curve of logarithmic slope equal to one is plotted alongside for
comparison. Given the likely stagnation behavior of the errors at some values
of the contribution ratio, the regression is only applied to the first points
where variations are noticed.

It can be seen that the logarithmic slops for the SM-ROM errors are close to
one meaning that the rate of convergence of the errors with respect to
$\mathcal{R}_{u}$ is close to $1$. In other words, a jump of one order of
magnitude in the ratio $\mathcal{R}_{u}$ results in a jump of the same order
in the error. This observation is in good agreement with Figure 4 which shows
that between POD ratios $10^{-2}$ and $10^{-4}$, the SM-ROM jumps two orders
of magnitudes and attains $4\times 10^{-4}$ for velocity and $2\times 10^{-2}$
for pressure. These errors represent a sufficient accuracy to reproduce the
dynamics of the flow and thus meet the high fidelity model hydrodynamic
coefficients.

Now by inspecting the estimator, one can see that even though the slopes are
slightly smaller, it provides a good indicator of the convergence rate of the
SM-ROM. If we look at the regression of the estimated errors at ratio
$10^{-4}$, it indicates a value of $\frac{3}{2}\times 10^{-3}$ for velocity
and $6\times 10^{-2}$ for the pressure. These are sufficiently good
approximations of the actual SM-ROM error that confirm the robustness of the
proposed estimator.

Figure 1: Two-dimensional domain and boundary conditions for the problem of
flow past a circular cylinder.

Figure 2: Illustration of Von Karman vortices arising in the flow past a
cylinder for $Re=200$ at three different time instants. Figure 3: Evolution of
the velocity and pressure POD contribution ratio with respect to modes numbers
for the flow past a cylinder at $Re=200$. Figure 4: Hydrodynamics coefficients
Figure 5: Log-Log scale representation of the estimator versus SM-ROM errors
for the flow past a cylinder at $Re=200$. The dots refer to the calculated
values while the solid lines are their power regression. The green solid line
corresponds to the curve of logarithmic slope equal to $1$. Only the errors
corresponding to $\mathcal{R}_{u}$ larger than $10^{-6}$ have been taken into
account to calculate the regression lines.

###  6.2 Test case 2 : Lid driven cavity Re = 9500

In this test case we assess the error estimates on the two dimensional lid
driven cavity flow problem with $Re=9500$. The problem domain consists of a
square cavity $]0,D[\times]0,D[$ where the fluid is driven by a tangential
velocity $U$ applied at the top wall. The remaining walls are defined as no-
slip conditions. The high fidelity problem is solved in a triangular mesh
composed of $32928$ cells by using the semi-implicit Euler method with a fixed
time step $\Delta t=10^{-3}$.

The dynamics of the flow solutions is shown in figure 6. We observe that in
the lower left corner, the secondary vortex separates into two small vortexes
that reincorporate in a periodic manner. In the same way the secondary vortex
in the upper left corner has a similar behavior. In order to construct the POD
bases, $200$ snapshots uniformly selected from the periodic regime with $10$
periods are used. In Figure 7, we show the evolution of the contribution ratio
of POD modes with respect to modes numbers. The effect of the POD contribution
ratio on the SM-ROM solution is illustrated in Figure 8. It shows the phase
portraits111The evolution of the horizontal velocity in terms of the vertical
velocity at points located at the bottom left $(2/16,2/16)$, top left
$(2/16,13/16)$ and top right $(0.95,0.95)$ of the cavity. Due to the periodic
nature of the flow at these points, the curves of the phase portraits are
closed and follow the same path over the periods of the flow.

Starting from a POD contribution ratio of order $10^{-2}$, it can be seen that
the SM-ROM phase portraits are not periodic and do not match the high fidelity
solution. By decreasing the ratio to $10^{-4}$, the SM-ROM phase portraits
tend to converge to the high fidelity solution except for the the bottom left
point of the cavity where small differences can still be noticed. These
differences disappear with ratios of order $10^{-5}$ or less. This shows that
the predictions by the SM-ROM converge to the groundtruth solution with
decreasing values of the POD contribution ratio.

In terms of error estimates, we observe from Figure 9, that the estimator
succeeds to predict the SM-ROM error slopes with a sufficiently good accuracy.
The logarithmic slopes for the SM-ROM errors are close to one which gives an
approximate rate of convergence of $10^{-1}$. This is in good agreement with
Figure 8 since between POD ratios $10^{-2}$ and $10^{-5}$, the SM-ROM jumps
three orders of magnitudes and attains $10^{-5}$ for velocity and $4\times
10^{-6}$ for pressure. These errors represent a sufficient accuracy to
reproduce well the dynamics of the flow and thus meet the high fidelity model
hydrodynamics coefficients. As for the flow past a cylinder, the slopes of the
estimators are smaller than the SM-ROM slopes. Nevertheless, they provide a
good indicator of the convergence rate of the SM-ROM. Typically, at ratio
$10^{-5}$, the estimator indicates an error value of $\frac{3}{2}\times
10^{-6}$ for velocity and $2\times 10^{-5}$ for the pressure. These are
sufficiently good approximations of the SM-ROM velocity and pressure errors of
values $\frac{3}{2}\times 10^{-7}$ and $4\times 10^{-6}$ respectively. Which
confirms once again the robustness of the proposed estimator.

Figure 6: Streamlines in three different instants of the periodic regime of
the flow in a lid driven cavity for $Re=9500$. Figure 7: Evolution of the
velocity and pressure POD contribution ratio with respect to modes numbers for
the Lid driven cavity problem at $Re=9500$. Figure 8: Velocity phase portrait
at the corner of the cavity domain for $Re=9500$. $u$ and $v$ are respectively
the horizontal and vertical components of the velocity vector $\bm{u}$. Figure
9: Log-Log scale representation of the estimator versus SM-ROM errors for the
flow in a lid driven cavity at $Re=9500$. The dots refer to the calculated
values while the solid lines are their power regression. The green solid line
corresponds to the curve of logarithmic slop equal to $1$. Only the errors
corresponding to $\mathcal{R}_{u}$ larger than $10^{-5}$ have been taken into
account to calculate the regression lines.

##  7 Summary and conclusions

In this paper we have carried on the numerical analysis of the stabilized
motivated reduced order (SM-ROM) solution of the pressure in incompressible
flows. We treat the pressure equation by duality of the momentum conservation
equation with gradients of the reduced pressure space, as function tests. We
have used the Proper Orthogonal decomposition (POD) method with $L^{2}$ norms
to build the reduced spaces for both velocity and pressure, and proven the
stability of the velocity-pressure segregated SM-ROM discretization in adapted
norms. In particular for pressure, the stability is proven in a discrete norm
that has the same asymptotic order as the $L^{2}$ norm with respect to the
mesh size. Both estimates for velocity and pressure are optimal with respect
to the projection error on the POD reduced spaces. To validate the theoretical
results, we have performed two two-dimensional numerical tests (lid-driven
cavity flow and for flow past a cylinder), for which we nearly recover the
optimality of the error estimates in both problems. We observe an error
stagnation effect for high ROM-POD dimensions, that possibly arises because
the weak high-fidelity model is not exact for gradients of discrete pressures.
Besides previous numerical investigation of the SM-ROM, the present analysis
supports the use of this method to recover the pressure from ROMs of
incompressible fluids that only retain weakly divergence-free velocities, with
accuracy enough for practical applications of interest.

Acknowledgments: This work has been partially supported by the Spanish
Government-EU Feder grant RTI2018-093521-B-C31. T. Chacón deeply thanks the
University of La Rochelle for its invitation in October of 2019 that allowed
to start this research. The research of Samuele Rubino has been also funded by
the Spanish State Research Agency through the national program Juan de la
Cierva-Incorporación 2017.

## References

  * [1] N. Ahmed, T. Chacón Rebollo, V. John, and S. Rubino.  Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method.  IMA J. Numer. Anal., 37(3):1437–1467, 2017.
  * [2] A. AL-Dmour and K. Mohammad.  Active control of flexible structures using principal component analysis in the time domain.  Journal of Sound and Vibration, 253, 2002.
  * [3] C. Allery, C. Beghein, and A. Hamdouni.  On investigation of particle dispersion by a POD approach.  International Applied Mechanics, 44(1):110–119, 2008,  https$:$//doi.org/10.1007/s10778-008-0025-2.
  * [4] C. Allery, C. Béghein, and A. Hamdouni.  Applying Proper orthogonal decomposition to the computation of particle dispersion in a two-dimensional ventilated cavity.  Communications in Nonlinear Science and Numerical Simulation, 10:907–920, 2005.
  * [5] C. Allery, C. Béghein, and A. Hamdouni.  Applying Proper Orthogonal Decomposition to the computation of particle dispersion in a two-dimensional ventilated cavity.  Communications in Nonlinear Science and Numerical Simulation, 10(8):907–920, 2005,  https$:$//doi.org/10.1016/j.cnsns.2004.05.005.
  * [6] F. Ballarin, T. Chacón Rebollo, E. Delgado Ávila, M. Gómez Mármol, and G. Rozza.  Certified reduced basis VMS-Smagorinsky model for natural convection flow in a cavity with variable height.  Comput. Math. Appl., 80(5):973–989, 2020.
  * [7] F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza.  Supremizer stabilization of podgalerkin approximation of parametrized steady incompressible navier-stokes equations.  Internat. J. Numer. Methods Engrg., 102:1136–1161, 2015.
  * [8] G. Berkooz, P. Holmes, and J. L. Lumley.  On the relation between low-dimensional models and the dynamics of coherent structures in the turbulent wall layer.  Theor. Comput. Fluid Dyn., 4:255–269, 1993.
  * [9] C. Bernardi, Y. Maday, and F. Rapetti.  Discrétisations variationnelles de problèmes aux limites elliptiques, volume 45 of Mathématiques & Applications.  Springer-Verlag, 2004.
  * [10] H. Brezis.  Functional analysis, Sobolev spaces and partial differential equations.  Universitext. Springer, New York, 2011.
  * [11] F. Brezzi and R. S. Falk.  Stability of higher-order Hood-Taylor methods.  SIAM J. Numer. Anal., 28(3):581–590, 1991.
  * [12] V. Buljak and G. Maier.  Proper orthogonal decomposition and radial basis functions in material characterization based on instrumented indentation.  Engineering Structures, 33, 2011.
  * [13] C. Béghein, C. Allery, J. Pozorski, and M. Waclawczyk.  Application of POD-based dynamical systems to dispersion and deposition of particles in turbulent channel flow.  International Journal of Multiphase Flow, 58:97 – 113, 2014,  https$:$//doi.org/10.1016/j.ijmultiphaseflow.2013.09.001.
  * [14] C. Béghein, C. Allery, M. Wacławczyk, and J. Pozorski.  Application of POD-based dynamical systems to dispersion and deposition of particles in turbulent channel flow.  International Journal of Multiphase Flow, 58:97–113, 1 2014.
  * [15] A. Caiazzo, T. Iliescu, V. John, and S. Schyschlowa.  A numerical investigation of velocity-pressure reduced order models for incompressible flows.  J. Comput. Phys., 259:598–616, 2014.
  * [16] T. Chacón Rebollo, M. Gómez Mármol, V. Girault, and I. Sánchez Muñoz.  A high order term-by-term stabilization solver for incompressible flow problems.  IMA J. Numer. Anal., 33(3):974–1007, 2013.
  * [17] T. Chacón Rebollo and R. Lewandowski.  Mathematical and numerical foundations of turbulence models and applications.  Modeling and Simulation in Science, Engineering and Technology. Birkhäuser/Springer, New York, 2014.
  * [18] P. Ciarlet.  The finite element method for elliptic problems. Studies in Mathematics and its Applications, Vol. 4. North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.
  * [19] J. de Frutos, B. García-Archilla, V. John, and J. Novo.  Analysis of the grad-div stabilization for the time-dependent Navier-Stokes equations with inf-sup stable finite elements.  Adv. Comput. Math., 44(1):195–225, 2018.
  * [20] J. de Frutos, B. García-Archilla, and J. Novo.  Fully discrete approximations to the time-dependent Navier-Stokes equations with a projection method in time and grad-div stabilization.  J. Sci. Comput., 80(2):1330–1368, 2019.
  * [21] J. de Frutos, B. García-Archilla, and J. Novo.  Corrigenda: Fully discrete approximations to the time-dependent Navier-Stokes equations with a projection method in time and grad-div stabilization.  J. Sci. Comput., 88(2):Paper No. 40, 3, 2021.
  * [22] L. Fick, Y. Maday, A. T. Patera, and T. Taddei.  A stabilized POD model for turbulent flows over a range of Reynolds numbers: optimal parameter sampling and constrained projection.  J. Comput. Phys., 371:214–243, 2018.
  * [23] S. Giere, T. Iliescu, V. John, and D. Wells.  SUPG reduced order models for convection-dominated convection-diffusion-reaction equations.  Comput. Methods Appl. Mech. Engrg., 289:454–474, 2015.
  * [24] S. Han and B. Feeny.  Application of Proper orthogonal decomposition to structural vibration analysis.  Mechanical Systems and Signal Processing, 17, 2003.
  * [25] J. S. Hesthaven, G. Rozza, and B. Stamm.  Certified reduced basis methods for parametrized partial differential equations.  SpringerBriefs in Mathematics. Springer, Cham; BCAM Basque Center for Applied Mathematics, Bilbao, 2016.  BCAM SpringerBriefs.
  * [26] S. Hijazi, G. Stabile, A. Mola, and G. Rozza.  Data-driven POD-Galerkin reduced order model for turbulent flows.  J. Comput. Phys., 416:109513, 30, 2020.
  * [27] P. Holmes, J. L. Lumley, and G. Berkooz.  Turbulence, coherent structures, dynamical systems and symmetry.  Cambridge Monographs on Mechanics. Cambridge University Press, Cambridge, 1996.
  * [28] T. Iliescu and Z. Wang.  Variational multiscale proper orthogonal decomposition: Navier-Stokes equations.  Numer. Methods Partial Differential Equations, 30(2):641–663, 2014\.
  * [29] K. Kean and M. Schneier.  Error analysis of supremizer pressure recovery for pod based reduced-order models of the time dependent navier-stokes equations.  SIAM J. Numer. Anal., 58 (4):2235–2264, 2020.
  * [30] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for parabolic problems.  Numer. Math., 90(1):117–148, 2001.
  * [31] H. P. Langtangen and A. Logg.  Solving PDEs in Python: The FEniCS Tutorial I.  Springer Publishing Company, Incorporated, 1st edition, 2017.
  * [32] C. Leblond, C. Allery, and C. Inard.  An optimal projection method for the reduced-order modeling of incompressible flows.  Computer Methods in App. Mechanics and Engineering, 200(33-36):2507–2527, 2011,  https$:$//doi.org/10.1016/j.cma.2011.04.020.
  * [33] D. J. Lucia, P. S. Beran, and W. A. Silva.  Reduced-order modeling: new approaches for computational physics.  Progress in Aerospace Sciences, 40, 2004.
  * [34] J. Novo and S. Rubino.  Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows.  SIAM J. Numer. Anal., 59(1):334–369, 2021.
  * [35] S. Park, J.-J. Lee, C.-B. Yun, and D. J. Inman.  Electro-mechanical impedance-based wireless structural health monitoring using PCA-Data compression and k-means Clustering Algorithms.  Journal of Intelligent Material Systems and Structures, 19, 05 2007\.
  * [36] A. Quarteroni and G. Rozza.  Numerical solution of parametrized navier–stokes equations by reduced basis methods.  Numerical Methods for Partial Differential Equations, 23(4):923–948, 2007.
  * [37] C. W. Rowley, T. Colonius, and R. M. Murray.  Model reduction for compressible flows using pod and galerkin projection.  Phys. D: Nonlinear Phenom., 189 (1):115–129, 2004.
  * [38] G. Rozza and K. Veroy.  On the stability of the reduced basis method for stokes equations in parametrized domains.  Comput. Methods Appl. Mech. Engrg., 196:1244–1260, 2007.
  * [39] S. Rubino.  Numerical analysis of a projection-based stabilized POD-ROM for incompressible flows.  SIAM J. Numer. Anal., 58(4):2019–2058, 2020.
  * [40] L. Sirovich.  Turbulence and the dynamics of coherent structures : Part I, II and III.  Quarterly of Applied Mathematics, pages 461–590, 1987.
  * [41] A. Tallet, C. Allery, and C. Leblond.  Optimal flow control using a POD-based reduced-order model.  Numerical Heat Transfer, Part B: Fundamentals, 70(1):1–24, 2016,  https$:$//doi.org/10.1080/10407790.2016.1173472.
  * [42] 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-3):909–932, 2015,  https$:$//doi.org/10.1016/j.cnsns.2014.09.009.
  * [43] C. Taylor and P. Hood.  A numerical solution of the Navier-Stokes equations using the finite element technique.  Internat. J. Comput. & Fluids, 1(1):73–100, 1973.
  * [44] J. P. Thomas, E. H. Dowell, and K. C. Hall.  Three-dimensional transonic aeroelasticity using Proper orthogonal decomposition-based reduced-order models.  Journal of Aircraft, 40, 05 2003.