# On a certified VMS-Smagorinsky Reduced Basis model with LPS pressure
stabilization

Tomás Chacón Rebollo Departamento de Ecuaciones Diferenciales y Análisis
Numérico, and Instituto de Matemáticas de la Universidad de Sevilla (IMUS) ,
Apdo. de correos 1160, Universidad de Sevilla, 41080 Seville, Spain.
<EMAIL_ADDRESS>Enrique Delgado Ávila Departamento de Ecuaciones Diferenciales y
Análisis Numérico, Apdo. de correos 1160, Universidad de Sevilla, 41080
Seville, Spain<EMAIL_ADDRESS><EMAIL_ADDRESS>Macarena Gómez
Mármol22footnotemark: 2

###### Abstract

In this work we present the numerical analysis of a Reduced Basis VMS-
Smagorinsky model with local projection stabilization (LPS) on the pressure.
We construct the reduced velocity space by two different strategies, by
considering or not the enrichment of the reduced velocity space with the so-
called inner pressure supremizer. We present the development of an a
posteriori error estimator for the snapshot selection through a Greedy
algorithm, based on the Brezzi-Rappaz-Raviart (BRR) theory. Moreover, the
Empirical Interpolation Method (EIM) is considered for the approximation of
the non-linear terms. Finally, we present some numerical tests in which we
show an speedup on the computation of the reduced basis problem with the LPS
pressure stabilisation, with respect to the method using pressure supremizers.

Keywords: Reduced basis method, Emprical interpolation method, a posteriori
error estimation, VMS-Smagorinsky model, LPS pressure stabilization.

##  1 Introduction

In this work we present a reduced basis Variational Multi-Scale (VMS)
Smagorinsky turbulence model with local projection stabilisation (LPS) of the
pressure, and small-small VMS modelling of eddy viscosity. We address a
turbulence model to consider realistic situations, as turbulent flows
frequently appear in actual applications. The Smagorinsky model is a basic
Large Eddy Simulation (LES) turbulence model, that provides particularly
accurate solutions when the eddy viscosity acts on the resolved small scales
(the small-small VMS setting, cf. [21]).

Stabilised methods for finite element discretisations provide stable
discretisations of the pressure that allow to circumvent the Brezzi-Babuska
inf-sup condition. This occurs with both residual-based methods and penalty
methods. As residual-based stabilisation methods let us mention the Galerkin-
Least Squares (GALS) (cf. [20]) and modifications of this method as the
Streamline Upwind Galerkin (SUPG) method (cf. [8, 23]), or the Adjoint-
Stabilised method (cf. [17]). The stabilising effect of those methods is due
to the dissipative effect of the terms added to the Galerkin projection. In
particular, the pressure stabilisation is granted thanks to a weighted
pressure Laplacian term. The penalty methods are simplified versions that just
keep the dissipative terms, although loosing full consistency. We may mention
the penalty term-by-term stabilised method (cf. [10]), that is an extension of
the pure penalty method introduced in [6]. In the latter method, the penalty
term is just a scaled pressure Laplacian. Further, the Local Projection
Stabilisation (LPS) treatment of the pressure is a penalty method that just
stabilises the high-frequency components of the pressure gradient that cannot
be represented in the discrete velocity space. It allows to use non-stable
pairs of velocity-pressure finite elements, such as
$\mathbb{P}_{k}-\mathbb{P}_{k}$, $k\geq 2$, while keeping optimal order of
accuracy in both velocity and pressure (cf. [10, 12]).

There exist several well-established techniques to stabilise the
discretisation of the reduced pressure in reduced order models (ROM) of
incompressible flows. Let us mention to enrich the reduced velocity space with
reduced velocity divergence - reduced pressure supremizers (cf. [2, 24]), or
solving a Poisson equation for the reduced pressure –by testing the momentum
conservation equation with gradients of reduced pressures–(cf. [9, 26]), or
using SUPG stabilisation (cf. [1]). These techniques provide efficient
approximations of the reduced pressure, although are somewhat involved,
requiring either increasing the dimension of the reduced velocity space,
building an equation for the pressure with a number of terms, or increasing
the complexity of the discrete problem. The LPS stabilisation provides a
simplified alternative to these techniques. The POD-ROM solution of Navier-
Stokes equations with LPS stabilisation of the pressure has recently been
addressed in [22], which includes an a-priori error analysis of the method.

In this work, we address the a posteriori error analysis - based reduced order
modelling of incompressible flow equations with LPS stabilisation of the
pressure. This allows to construct “certified” solvers, with error level below
targeted values. We actually afford the reduced basis solution of the steady
Smagorinsky turbulence model for the sake of adressing realistic flows. We
include the small-small VMS formulation of the eddy diffusion which, actually,
is a three-level LES turbulence model (cf. [4]) that provides good accuracy in
the approximation of first and second order turbulence statistics (cf. [13]).
We prove the stability of the reduced problem with LPS stabilisation of the
pressure for velocity-pressure finite elements, whenever the pressure is
piecewise affine. We construct an a posteriori error bound estimator using the
Brezzi-Raviart-Rappaz (BRR) theory of approximation of regular branches of
solutions of non-linear equations (cf. [7]). This is the key to build a
certified RB solver. We use the Empirical Interpolation Method (EIM) to build
reduced approximations of the non-linear coefficients appearing in model: the
eddy diffusion and the statibilised coefficient appearing in the LPS
formulation. We present some numerical results were we compare the velocity
supremizers approach versus the LPS pressure discretization. We observe an
improved speedup for the LPS discretization with similar error levels.

The structure of this work is as follows. In section 2, we present the VMS-
Smagorinsky model with local projection stabilisation for the pressure. In
subsection 2.1, we present the FE model, while the reduced basis method is
considered in subsection 2.2. The numerical analysis for the well-posedness is
presented in section 3, and the sub-sequent development of the a posteriori
error bound estimator in section 4. Then, in section 5, we present the
approximation of the eddy-viscosity and the pressure stabilisation constant by
the EIM. Finally, we present the numerical tests in section 6, in which we
highlight the speedup rates for the computation of solution with the present
RB model presented in this work.

##  2 Pressure local projection stabilised VMS-Smagorinsky model

In this section we describe a local projection-stabilised VMS-Smagorinsky
model. In the turbulence LES (Large Eddy Simulation) high fidelity model that
we consider, the eddy-viscosity only acts on the small resolved scales. This
leads to a less diffusive model than the Smagorinsky one presented in [11],
where the eddy viscosity acts on both large and small resolved scales. In
[13], it has been shown that this model has a good precision for the
approximation of the mean flow and the second-order statistics of the
turbulent flow. Moreover, we use a LPS stabilisation treatment for the
pressure, that allows us to use non-stable pairs of finite elements.

###  2.1 Finite element problem

Let $\Omega$ be a bounded domain of $\mathbb{R}^{d}\,(d=2,3)$, with Lipschitz-
continuous boundary $\Gamma$. We assume that $\Gamma$ is split into
$\Gamma=\Gamma_{D}\cup\Gamma_{N}$ where $\Gamma_{D}$ and $\Gamma_{N}$ are two
connected measurable sets of positive $(d-1)$-dimensional measure, with
disjoint interiors. We intend to impose Dirichlet and Neumann boundary
conditions on $\Gamma_{D}$ and $\Gamma_{N}$, respectively.

We present a parametric steady Smagorinsky model, where we consider the
Reynolds number as a parameter, denoted by $\mu\in\mathcal{D}$, where
$\mathcal{D}$ is a compact sub-set of $\mathbb{R}$. Although the Smagorinsky
model is intrinsically discrete, we present it in a continuous form in order
to clarify its relationship with the Navier-Stokes equations: We search for a
velocity field $\textbf{u}\,:\Omega\mapsto\mathbb{R}^{d}$ and a pressure
function $p\,:\Omega\mapsto\mathbb{R}$ such that

$\left\\{\begin{array}[]{ll}\textbf{u}\cdot\nabla\textbf{u}+\nabla
p-\nabla\cdot\left(\dfrac{1}{\mu}\nabla\textbf{u}+\nu_{T}(\textbf{u}^{\prime})\nabla\textbf{u}^{\prime}\right)=\textbf{f}&\mbox{
in }\Omega,\vspace{0.2cm}\\\ \nabla\cdot\textbf{u}={\bf 0}&\mbox{ in
}\Omega,\vspace{0.1cm}\\\ \textbf{u}=0&\mbox{ on }\Gamma_{D},\vspace{0.1cm}\\\
\textbf{n}\cdot\left(\dfrac{1}{\mu}\nabla\textbf{u}+\nu_{T}(\textbf{u}^{\prime})\nabla\textbf{u}^{\prime}\right)=0&\mbox{
on }\Gamma_{N},\vspace{0.2cm}\\\ \end{array}\right.$ (1)

where f is the kinetic momentum source, $\nu_{T}(\textbf{u}^{\prime})$ is the
eddy viscosity defined as

$\nu_{T}(\textbf{u}^{\prime})=C_{S}^{2}\displaystyle\sum_{K\in\mathcal{T}_{h}}h_{K}^{2}\big{|}\nabla\textbf{u}_{|_{K}}^{\prime}\big{|}\chi_{K},$
(2)

where $\big{|}\cdot\big{|}$ denotes the Frobenius norm in $\mathbb{R}^{d\times
d}$, $C_{S}$ is the Smagorinsky constant and $\textbf{u}^{\prime}$ stands for
the small-scales part of the resolved velocity field u, $\textbf{u}^{\prime}$
is assumed to lie in the inertial spectrum of the full velocity field. The
eddy viscosity setting in (1)-(2) corresponds to the small-small VMS setting
introduced in [21]. Other possibilities would correspond to replacing
$\textbf{u}^{\prime}$ by u in either (1) or (2). The eddy viscosity in (1)
only acts on the small scales of the resolved part of the velocity field. Let
$\\{\mathcal{T}_{h}\\}_{h>0}$ is uniformly regular mesh in the sense of
Ciarlet [Ciarlet1978].

To state the full order discretization that we consider for problem (1), let
us introduce the spaces

$Y=\\{\textbf{v}\in H^{1}(\Omega)^{d},\,\,\mbox{s.t.
}\textbf{v}_{|_{\Gamma_{D}}}={\bf 0}\,\\},\quad M=\\{q\in
L^{2}(\Omega),\,\,\mbox{s.t. }\,\int_{\Omega}q=0\,\\}.$

We assume $\textbf{f}\in Y^{\prime}$.

Given an integer $n\geq 1$ and an element $K\in\mathcal{T}_{h}$, we denote by
$\mathbb{P}_{n}(K)$ the Finite Element space given by Lagrange polynomials of
degree less than, or equal to, $n$ defined on $K$. We assume that $\Omega$ is
polygonal when $d=2$ or polyhedric when $d=3$, and consider a triangulation
$\mathcal{T}_{h}$ such that $\Gamma_{D}$ and $\Gamma_{N}$ are unions of full
simplices of the boundaries of elements of $\mathcal{T}_{h}$. Let us define
the finite element space

$V^{(n)}=\\{v_{h}\in C^{0}(\overline{\Omega})\text{ such that
}v_{h|_{K}}\in\mathbb{P}_{n}(K),\forall K\in\mathcal{T}_{h}\\}.$ (3)

Let $Y_{h}=Y\cap[V^{(n)}]^{d}$ and $M_{h}\subset M\cap V^{(m)}$ be two finite
element subspaces of $Y$ and $M$, respectively. We introduce the LPS-VMS
discretisation of the Smagorinsky model (1), with stabilisation on the
pressure,

$\left\\{\begin{array}[]{l}\mbox{Find
}(\textbf{u}_{h},p_{h})=(\textbf{u}_{h}(\mu),p_{h}(\mu))\in Y_{h}\times
M_{h}\mbox{ such that}\vspace{0.3cm}\\\
\begin{array}[]{ll}a(\textbf{u}_{h},\textbf{v}_{h};\mu)+b(\textbf{v}_{h},p_{h};\mu)+a_{S}^{\prime}(\textbf{u}_{h};\textbf{u}_{h},\textbf{v}_{h};\mu)\\\
+c(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h};\mu)=\left<\textbf{f},\textbf{v}\right>&\quad\forall\textbf{v}_{h}\in
Y_{h},\vspace{0.1cm}\\\
b(\textbf{u}_{h},q_{h};\mu)+s_{pres}(p_{h},q_{h};\mu)=0&\quad\forall q_{h}\in
M_{h},\end{array}\end{array}\right.$ (4)

where the bilinear forms $a(\cdot,\cdot;\mu)$ and $b(\cdot,\cdot;\mu)$ are
defined as

$a(\textbf{u},\textbf{v};\mu)=\frac{1}{\mu}\int_{\Omega}\nabla\textbf{u}:\nabla\textbf{v}\,d\Omega,\qquad
b(\textbf{v},q;\mu)=-\int_{\Omega}(\nabla\cdot\textbf{v})q\,d\Omega;$ (5)

while the trilinear form $c(\cdot,\cdot,\cdot;\mu)$ is defined as

$c(\textbf{z},\textbf{u},\textbf{v};\mu)=\frac{1}{2}\,\left[\int_{\Omega}(\textbf{z}\cdot\nabla\textbf{u})\textbf{v}\,d\Omega-\int_{\Omega}(\textbf{z}\cdot\nabla\textbf{v})\textbf{u}\,d\Omega\,\right].$
(6)

Moreover, the non-linear form $a^{\prime}_{S}(\cdot;\cdot,\cdot;\mu)$, is a
multi-scale Smagorinksy modelling for the eddy viscosity term, and it is given
by

$a_{S}^{\prime}(\textbf{z}_{h};\textbf{u}_{h},\textbf{v}_{h};\mu)=\displaystyle\int_{\Omega}\nu_{T}(\Pi_{h}^{*}\textbf{z}_{h})\nabla(\Pi_{h}^{*}\textbf{u}_{h}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\;d\Omega,$
(7)

where $\Pi_{h}^{*}=Id-\Pi_{h}$, $\Pi_{h}$ being a uniformly stable (in $L^{2}$
and $H^{1}$-norms) interpolation operator on $Y\cap[V_{h}^{(n-1)}]^{d}$. This
interpolation operator must satisfy optimal error estimates (cf. [5] Chap. IX)
and preserve the homogeneous Dirichlet boundary conditions when restricted to
$Y_{h}$. In this way the small-scale component of the velocity field that we
consider is defined by $\textbf{u}_{h}^{\prime}=\Pi_{h}^{*}\textbf{u}_{h}$.
The eddy viscosity only acts on $\Pi_{h}^{*}\textbf{u}_{h}$.

In VMS terminology, the Large-Small setting would corresponds to modeling the
turbulent viscosity by a function of the whole resolved velocity, taking

$a_{S}^{\prime}(\textbf{z}_{h};\textbf{u}_{h},\textbf{u}_{h})=\displaystyle\int_{\Omega}\nu_{T}(\textbf{z}_{h})\nabla(\Pi_{h}^{*}\textbf{u}_{h}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\;d\Omega.$
(8)

Further, the term $s_{pres}(\cdot,\cdot;\mu)$ is the projection-stabilisation
term for the pressure, defined as

$s_{pres}(p_{h},q_{h};\mu)=\displaystyle\sum_{K\in\mathcal{T}_{h}}\tau_{p,K}(\mu)(\sigma_{h}^{*}(\nabla
p_{h}),\sigma_{h}^{*}(\nabla q_{h}))_{K}.$ (9)

Here, $\sigma_{h}^{*}=Id-\sigma_{h}$, with $\sigma_{h}$ some locally stable
$L^{2}$ projection or interpolator operator from $L^{2}(\Omega)$ to
$[V_{h}^{(n)}]^{d}$. Moreover, $\tau_{p,K}(\mu)$ is the stabilisation
coefficient in (9), that must verify the following hypotesis:

######  Hypothesis 1.

There exists two positive constants $\alpha_{1},\alpha_{2}$, independent of
$h$, such that

$\alpha_{1}h_{K}^{2}\leq\tau_{p,K}(\mu)\leq\alpha_{2}h_{K}^{2},\quad\forall\mu\in\mathcal{D},\forall
K\in\mathcal{T}_{h},\forall h>0.$ (10)

In the following, we will use the stabilisation pressure coefficient
$\tau_{p,K}(\mu)$ proposed by Codina in [15], and used in [13, 25]:

$\tau_{p,K}(\mu)=\left[c_{1}\dfrac{1/\mu+\overline{\nu}_{T_{|K}}(\mu)}{h_{K}^{2}}+c_{2}\dfrac{U_{K}(\mu)}{h_{K}}\right]^{-1},$
(11)

where $\overline{\nu}_{T_{|K}}$ is some local eddy viscosity, $U_{K}$ is a
local velocity and $c_{1},c_{2}$ some positive experimental constants. Taking
$\tau_{p,K}(\mu)$ by this way, we are ensuring (10).

By denoting $X_{h}=Y_{h}\times M_{h}$, we rewrite problem (4), with compact
notation, as follows

$\left\\{\begin{array}[]{l}\mbox{Find }U_{h}(\mu)=(\textbf{u}_{h},p_{h})\in
X_{h}\mbox{ such that}\vspace{0.2cm}\\\
A(U_{h}(\mu),V_{h};\mu)=F(V_{h};\mu)\qquad\forall V_{h}\in
X_{h},\end{array}\right.$ (12)

where

$\begin{array}[]{ll}A(U_{h}(\mu),V_{h};\mu)&=\dfrac{1}{\mu}A_{Diff}(U_{h},V_{h})+A_{Div}(U_{h},V_{h})+A_{Conv}(U_{h};V_{h})\vspace{0.1cm}\\\
&+A_{Sma}(U_{h};V_{h})+A_{Pres}(U_{h},V_{h};\mu),\end{array}$ (13)

where, denoting $U=(\textbf{u},p^{u})$, $V=(\textbf{v},p^{v})$,

$\begin{array}[]{rll}A_{Diff}(U,V)&\\!\\!\\!\\!=\displaystyle\int_{\Omega}\nabla\textbf{u}:\nabla\textbf{v}\;d\Omega,\vspace{0.2cm}\\\
A_{Div}(U,V)&\\!\\!\\!\\!=\displaystyle\int_{\Omega}(\nabla\cdot\textbf{u})p^{v}\;d\Omega-\displaystyle\int_{\Omega}(\nabla\cdot\textbf{v})p^{u}\;d\Omega\vspace{0.1cm}\\\
A_{Conv}(U;V)&\\!\\!\\!\\!=\displaystyle\frac{1}{2}\,\displaystyle\int_{\Omega}\left((\textbf{u}\cdot\nabla\textbf{u})\textbf{v}-(\textbf{u}\cdot\nabla\textbf{v})\textbf{u}\,\right)\;d\Omega,\vspace{0.2cm}\\\
A_{Sma}(U;V)&\\!\\!\\!\\!=\displaystyle\int_{\Omega}\nu_{T}(\Pi_{h}^{*}\textbf{u})\,\nabla(\Pi_{h}^{*}\textbf{u}):\nabla(\Pi_{h}^{*}\textbf{v})\;d\Omega,\vspace{0.2cm}\\\
A_{Pres}(U,V;\mu)&\\!\\!\\!\\!=\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}\tau_{p,K}(\mu)\,\sigma_{h}^{*}(\nabla
p^{u})\cdot\sigma_{h}^{*}(\nabla p^{v})\;d\Omega,\vspace{0.2cm}\\\
F(V_{h};\mu)&\\!\\!\\!\\!=\langle\textbf{f},\textbf{v}_{h}\rangle.\end{array}$

The numerical analysis (existence of solutions, stability and error analysis)
of the full-order projection-based VMS-Smagorinsky model for problem (4) (or
equivalently (12)) with wall-law boundary conditions can be found in [13].

###  2.2 Reduced basis problem

In this section we introduce the Reduced Basis (RB) model. It is the
discretisation (12) with LPS stabilisation of the pressure built on the
reduced spaces. The RB problem reads

$\left\\{\begin{array}[]{l}\mbox{Find }(\textbf{u}_{N},p_{N})\in Y_{N}\times
M_{N}\mbox{ such that}\vspace{0.3cm}\\\
\begin{array}[]{ll}a(\textbf{u}_{N},\textbf{v}_{N};\mu)+b(\textbf{v}_{N},p_{N};\mu)+a_{S}^{\prime}(\textbf{u}_{N};\textbf{u}_{N},\textbf{v}_{N};\mu)\\\
+c(\textbf{u}_{N},\textbf{u}_{N},\textbf{v}_{N};\mu)=\langle\textbf{f},\textbf{v}_{N}\rangle&\quad\forall\textbf{v}_{N}\in
Y_{N},\vspace{0.1cm}\\\ b(\textbf{u}_{N},q_{N};\mu)=0&\quad\forall q_{N}\in
M_{N}.\end{array}\end{array}\right.$ (14)

Here, we denote by $Y_{N}$ the reduced velocity space, and by $M_{N}$ the
reduced pressure space. The computation of the reduced spaces is done through
the Greedy algorithm [16]. The Greedy algorithm summarizes:

  1. 1.

Set $\mu^{1}$, compute $(\textbf{u}_{h}(\mu^{1}),p_{h}(\mu^{1}))$ solving
problem (4), and define the reduced space
$Y_{1}=\mbox{span}\\{\textbf{u}_{h}(\mu^{1})\\}$,
$M_{1}=\mbox{span}\\{p_{h}(\mu^{1})\\}$. Orthonormalize the spaces
$Y_{1},M_{1}$.

  2. 2.

For $k\geq 2$, compute $\Delta_{k-1}(\mu),\forall\mu\in\mathcal{D}_{train}$
and set
$\mu^{k}=\arg\displaystyle\max_{\mu\in\mathcal{D}_{train}}\Delta_{k-1}(\mu).$
Here, we denote by $\Delta_{k}$ the a posteriori error estimator.

  3. 3.

Compute $(\textbf{u}_{h}(\mu^{k}),p_{h}(\mu^{k}))$, define the reduced spaces
$Y_{k},M_{K}$ by adding the new computed snapshot, and orhtonormalize them.

  4. 4.

Stop if
$\displaystyle\max_{\mu\in\mathcal{D}_{train}}\Delta_{k}(\mu)<\varepsilon_{RB}$.
If not, back to 2.

This yields reduced velocity and pressure spaces corresponding to the solution
of problem (14) for selected values of the Reynolds number:

$Y_{N}=\mbox{span}\\{\zeta^{\textbf{v}}_{i}:=\textbf{u}(\mu^{i}),\;i=1,\dots,N\\};$
(15) $M_{N}=\mbox{span}\\{\xi^{p}_{i}:=p(\mu^{i}),\;i=1,\dots,N\\},$ (16)

where the $\mu^{i}$ are selected by the Greedy algorithm. The stability of
this discretization is proved in the next section.

To test the numerical performances of this discretisation, we compare it with
the strategy consisting in enriching the reduced velocity space with the
supremizers of the reduced velocity divergence - reduced pressure duality. The
reduced velocity space in this case is given by

$Y_{N}=\mbox{span}\\{\zeta^{\textbf{v}}_{2i-1}:=\textbf{u}(\mu^{i}),\zeta^{\textbf{v}}_{2i}:=T_{p}^{\mu}(\xi^{p}_{i})),\;i=1,\dots,N\\},$
(17)

where here we are denoting $T_{p}^{\mu}:M_{h}\rightarrow Y_{h}$ the reduced
pressure supremizer, defined by

$\left(\nabla
T_{p}^{\mu}q_{h},\nabla\textbf{v}_{h}\right)_{\Omega}=b(q_{h},\textbf{v}_{h};\mu)\quad\forall\textbf{v}_{h}\in
Y_{h}.$ (18)

This strategy leads to an offline/online stabilisation where the dimension of
the reduced space $X_{N}=Y_{N}\times M_{N}$ in problem (14) is $3N$, while
using the LPS stabilisation the dimension of $X_{N}$ is $2N$.

When $Y_{N}$ is defined by (17), i. e., when it includes the reduced pressure
supremizers, then the discrete inf-sup condition between spaces $Y_{N}$ and
$M_{N}$ holds (cf. [2]). Consequently, for small enough data f, the reduced
problem (15) admits a solution by a finite-dimensional compactness argument,
similar to the one used in [13] to prove the existence of solutions of high
fidelity problem (4).

When the LPS approach is followed, this argument also holds when the
velocities are piecewise affine. Indeed, in this case we may bound a weaker
norm of the pressure,

$|||q_{N}|||_{N}=\sup_{\textbf{v}_{N}\in
Y_{N}}\dfrac{(q_{N},\nabla\cdot\textbf{v}_{N})_{\Omega}}{\|\nabla{\textbf{v}_{N}}\|_{0,2,\Omega}}+\|\sigma_{h}^{*}(\nabla
q_{N})\|_{\tau_{p}},\quad\forall q_{N}\in M_{N},$ (19)

where

$\|r\|_{\tau_{p}}=\left(\displaystyle\sum_{K\in\mathcal{T}_{h}}\tilde{\tau}_{p,K}\|r\|_{0,2,K}^{2}\right)^{1/2},\quad\mbox{for
}\,\,r\in L^{2}(\Omega),$

$\mbox{where
}\,\,\tilde{\tau}_{p,K}=\displaystyle\inf_{\mu\in\mathcal{D}}\tau_{p,K}(\mu)$,
and we denote by $\|\cdot\|_{l,p,S}$ the $W^{l,p}(S)$ norm on some open set
$S$. It holds

######  Proposition 1.

Assume that $M_{h}$ is formed by piecewise affine functions, and that the
space $Y_{N}$ defined by (15) contains some function $\textbf{w}_{N}$ such
that $\int_{\Omega}\textbf{w}_{N}\neq 0$. Then $|||\cdot|||_{N}$ is a norm on
$M_{N}$. Moreover, there exists a constant $C>0$ such that

$|||q_{N}|||_{N}\leq C\,\|q_{N}\|_{0,2,\Omega},\quad\forall q_{N}\in M_{N}.$
(20)

Proof. It is straightforward to prove that $|||q|||_{N}\geq 0$ for any $q\in
L^{2}(\Omega)$ and that $|||\cdot|||_{N}$ satisfies the triangular inequality.
Assume now that $|||q_{N}|||_{N}=0$ for some $q_{N}\in M_{N}$. Then
$\sigma_{h}^{*}(\nabla q_{N})=0$ and consequently $\nabla
q_{N}=\sigma_{h}(\nabla q_{N})\in[V_{h}^{n}]^{d}$. Then $\nabla q_{N}$ is a
continuous function and, as it is a constant on each element $K\in{\cal
T}_{h}$, then $\nabla q_{N}$ is constant on $\Omega$. Moreover, as
$\displaystyle\sup_{\textbf{v}_{N}\in
Y_{N}}\dfrac{(q_{N},\nabla\cdot\textbf{v}_{N})_{\Omega}}{\|\nabla{\textbf{v}_{N}}\|_{0,2,\Omega}}=0,$
then $(q_{N},\nabla\cdot\textbf{v}_{N})_{\Omega}=0$ for any $\textbf{v}_{N}\in
Y_{N}$. Then

$\nabla
q_{N}\,\int_{\Omega}\textbf{w}_{N}=-(q_{N},\nabla\cdot\textbf{w}_{N})_{\Omega}=0.$

Consequently, $\nabla q_{N}=0$ and as $\displaystyle\int_{\Omega}q_{N}=0$
because $q_{N}\in M_{h}$, it follows $q_{N}=0$.

Finally, estimate (20) readily follows from the inverse finite element
estimate (cf. [5])

$\|\nabla q_{N}\|_{0,2,K}\leq C\,h_{K}^{-2}\,\|q_{N}\|_{0,2,K}.$

∎

Then it follows

######  Corollary 1.

Under the hypotheses of Proposition 1, the reduced problem (14) admits at
least a solution, that satisfies the estimates

$\|\nabla\textbf{u}_{N}\|_{0,2,\Omega}+|||p_{N}|||_{N}\leq\Phi(\|\textbf{f}\|_{Y^{\prime}};\mu).$

for some positive continuous function $\Phi(\eta;\mu)$ increasing to $+\infty$
as $\eta\mapsto+\infty$ or $\mu\mapsto 0$.

Sketch of the proof. The proof of existence follows from a finite-dimensional
compactness argument using Brouwer’s fixed point theorem, similar to the one
used in [13] to prove the existence of solutions of the high fidelity problem
(4) with wall-law boundary conditions through a linearisation procedure. The
essential point is the bounds for velocity and pressure that on one hand imply
the uniqueness of solutions of the linearised problem and on another hand
guarantee the hypotheses of Brouwer’s fixed point theorem. Indeed, setting
$\textbf{v}_{N}=\textbf{u}_{N}$ and $q_{N}=-p_{N}$, and summing the two
equations appearing in (14), it follows

$\|\nabla\textbf{u}_{N}\|_{0,2,\Omega}+\|\sigma_{h}^{*}(\nabla
p_{N})\|_{\tau_{p}}\leq\Phi_{1}(\|\textbf{f}\|_{Y^{\prime}};\mu)$ (21)

for some positive continuous function $\Phi_{1}(\eta;\mu)$ increasing to
$+\infty$ as $\eta\mapsto+\infty$ or $\mu\mapsto 0$. Further, setting
$q_{N}=0$ in (14), it holds

$(\nabla\cdot\textbf{v}_{N},p_{N})_{\Omega}=a(\textbf{u}_{N},\textbf{u}_{N};\mu)+a_{S}^{\prime}(\textbf{u}_{N};\textbf{u}_{N},\textbf{v}_{N};\mu)+c(\textbf{u}_{N},\textbf{u}_{N},\textbf{v}_{N};\mu)-\langle\textbf{f},\textbf{v}_{N}\rangle$

for all $\textbf{v}_{N}\in Y_{N}$. Consequently,

$(\nabla\cdot\textbf{v}_{N},p_{N})_{\Omega}\leq\Phi_{2}(\|\textbf{f}\|_{Y^{\prime}};\mu)\,\|\nabla\textbf{v}_{N}\|_{0,2,\Omega}$

for some function $\Phi_{2}(\cdot;\cdot)$ with the same properties as
$\Phi_{1}$. Then,

$\sup_{\textbf{v}_{N}\in
Y_{N}}\dfrac{(\nabla\cdot\textbf{v}_{N},p_{N})_{\Omega}}{\|\nabla{\textbf{v}_{N}}\|_{0,2,\Omega}}\leq\Phi_{2}(\|\textbf{f}\|_{Y^{\prime}};\mu).$
(22)

Combining (21) and (22) it follows

$\|\nabla\textbf{u}_{N}\|_{0,2,\Omega}+|||p_{N}|||_{N}\leq\Phi_{1}(\|\textbf{f}\|_{Y^{\prime}};\mu)+\Phi_{2}(\|\textbf{f}\|_{Y^{\prime}};\mu).$

The same argument allows to prove the bounds for velocity and pressure in the
linearised problem and then to apply Brouwer’s fixed point theorem. We do not
describe the proof in detail for brevity. ∎

##  3 Well-posedness analysis

In this section we study the well-posedness of the full order (4) and of the
reduced problem (12), as well as of the linearised problem. This will be the
basis for the derivation of the a posteriori error estimation in the next
section.

We shall denote by $C>0$ along the paper the constants that may vary from line
to line but which are independent of $\nu$ and $h$. We need the following
basic results:

######  Lemma 1.

It holds

$|c(\textbf{z},\textbf{u},\textbf{v};\mu)|\leq
C_{S4}^{2}\,\|\textbf{z}\|_{0,2,\Omega}\,\|\textbf{u}\|_{0,2,\Omega}\,\|\textbf{v}\|_{0,2,\Omega}\,\forall\,\textbf{z},\,\textbf{u},\,\textbf{v}\in
Y,$ (23)

where $C_{S4}$ is the norm of the injection of $Y$ into $L^{4}(\Omega)^{d}$.

This is a standard result, we omit its proof.

######  Lemma 2.

For any $\textbf{u}_{h},\textbf{v}_{h},\textbf{z}_{h}\in Y_{h}$ and for any
$\mu\in\mathcal{D}$, the non-linear form
$a_{S}^{\prime}(\cdot;\cdot,\cdot;\mu)$ defined in (7) satisfies

$|a_{S}^{\prime}(\textbf{z}_{h};\textbf{u}_{h},\textbf{v}_{h};\mu)|\leq
C_{f}\,C_{S}^{2}h^{2-d/2}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega},$
(24)

for some constant $C_{f}>0$.

Proof. Using the local inverse finite element estimates (cf. [5])

$\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,\infty,K}\leq
C\,h_{K}^{-d/2}\,\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,K}$

for some constant $C>0$, it holds $|\nu_{T}(\Pi_{h}^{*}\textbf{z}_{h})({\bf
x})|\leq
C\,C_{S}^{2}\,h_{K}^{2-d/2}\,\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,K}$
for all ${\bf x}\in K$, for any $K\in{\cal T}_{h}$ . Then,

$\displaystyle|a_{S}^{\prime}(\textbf{z}_{h};\textbf{u}_{h},\textbf{v}_{h};\mu)|$
$\displaystyle\leq
C\,C_{S}^{2}\,h^{2-d/2}\,\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\,\sum_{K\in{\cal
T}_{h}}\|\nabla(\Pi_{h}^{*}\textbf{u}_{h})\|_{0,2,K}\,\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,K}$
$\displaystyle\leq
C\,C_{S}^{2}\,h^{2-d/2}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\,\|\nabla(\Pi_{h}^{*}\textbf{u}_{h})\|_{0,2,\Omega}\,\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}.$

Using the stability of operator $\Pi_{h}$ in $H^{1}(\Omega)^{d}$ norm,
estimate (24) follows.

∎

######  Lemma 3.

Assume Hypothesis 1 holds, and let $q_{h}\in V_{h}^{l}(\Omega)$. Then,

$\|\sigma_{h}^{*}(\nabla q_{h})\|_{\tau_{p}}\leq
C_{\tau}\|q_{h}\|_{0,2,\Omega}$ (25)

for some constant $C_{\tau}>0$.

Proof. Taking into account Hypothesis 1, and the local stability of operator
$\sigma$, it holds

$\tau_{K}\|{\sigma_{h}^{*}(\nabla q_{h})}\|_{0,2,K}^{2}\leq
Ch_{K}^{2}\|{\nabla q_{h}}\|_{0,2,K}^{2}$ (26)

Considering the local inverse inequalities for finite element functions (cf.
[5]) it follows

$\|{\nabla q_{h}}\|_{0,2,K}\leq C\,h_{K}^{-1}\,\|q_{h}\|_{0,2,K},\,\,\forall
K\in{\cal T}_{h}.$ (27)

Consequently,

$\|\sigma_{h}^{*}(\nabla q_{h})\|_{\tau_{p}}\leq C\,\|q_{h}\|_{0,2,\Omega}.$

∎

We shall use the following inf-sup generalised condition (cf. [12])

######  Proposition 2.

Assume that the family of triangulations $\\{\mathcal{T}_{h}\\}_{h>0}$ is
uniformly regular. Then it holds

$\|q_{h}\|_{0,2,\Omega}\leq\alpha\left(\sup_{\textbf{v}_{h}\in
Y_{h}}\dfrac{(q_{h},\nabla\cdot\textbf{v}_{h})_{\Omega}}{\|\nabla\textbf{v}_{h}\|_{0,2,\Omega}}+\|\sigma_{h}^{*}(\nabla
q_{h})\|_{\tau_{p}}\right)\quad\forall q_{h}\in M_{h},$ (28)

with $\alpha>0$, independent of $h$.

We will assume in the sequel for simplicity of derivations that the grids are
uniformly regular. The previous result can be extended to regular rather than
uniformly regular grids, with a more involved analysis. For further details,
see [12].

The following a priori estimates for the full order velocity and pressure
hold:

######  Theorem 1.

Let $(\textbf{u}_{h},p_{h})\in X_{h}$ a solution of problem (4). Let
$\textbf{f}\in Y_{h}^{\prime}$. Then, the following a priori estimates are
verified:

$\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}\leq\mu\,C_{P}\,\|\textbf{f}\,\|_{Y^{\prime}},$
(29) $\|\sigma_{h}^{*}(\nabla
p_{h})\|_{\tau_{p}}\leq\sqrt{\mu}\,C_{P}\,\|\textbf{f}\,\|_{Y^{\prime}},$ (30)
$\|p_{h}\|_{0,2,\Omega}\leq\alpha^{-1}\,\left(((1+\sqrt{\mu})\,C_{P}+1)\,\|\textbf{f}\,\|_{Y^{\prime}}+\frac{C}{\mu\,C_{P}}\,\|\textbf{f}\,\|_{Y^{\prime}}^{2}\right),$
(31)

where $C_{P}>0$ is the norm of the Poincaré injection of $Y$ into
$L^{2}(\Omega)$ and $C>0$ is a constant independent of $h$, $\mu$ and the data
f.

Proof. Estimates (29) and (30) are readily obtained taking
$\textbf{v}_{h}=\textbf{u}_{h}$ and $q_{h}=p_{h}$ in (4) and summing the two
identities appearing therein.

To prove (31), we take $\textbf{v}_{h}=\textbf{u}_{h}$ and $q_{h}=0$ in (4).
We have

$a(\textbf{u}_{h},\textbf{v}_{h};\mu)+a_{s}^{\prime}(\textbf{u}_{h};\textbf{u}_{h},\textbf{v}_{h};\mu)+c(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h};\mu)-(p_{h},\nabla\cdot\textbf{v}_{h})=\left<\textbf{f},\textbf{v}_{h}\right>.$

Thus, from Cauchy-Schwartz inequality, and Lemma 2,

$|(p_{h},\nabla\cdot\textbf{v}_{h})|=|a(\textbf{u}_{h},\textbf{v}_{h};\mu)+a_{s}^{\prime}(\textbf{u}_{h};\textbf{u}_{h},\textbf{v}_{h};\mu)+c(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h};\mu)-\left<\textbf{f},\textbf{v}_{h}\right>|$
$\leq\dfrac{1}{\mu}\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega}+C\,\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}^{2}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega}+\|\textbf{f}\,\|_{Y^{\prime}}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega},$

for some constant $C>0$. Therefore, from Proposition 2,

$\alpha\,\|p_{h}\|_{0,2,\Omega}\leq\sup_{\textbf{v}_{h}\in
Y_{h}}\dfrac{(p_{h},\nabla\cdot\textbf{v}_{h})_{\Omega}}{\|\nabla{\textbf{v}_{h}}\|_{0,2,\Omega}}+\|\sigma_{h}^{*}(\nabla
p_{h})\|_{\tau_{p}}$
$\leq\dfrac{1}{\mu}\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}+C\,\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}^{2}+(1+\sqrt{\mu}\,C_{P})\,\|\textbf{f}\,\|_{Y^{\prime}}$

Combining this estimate with (29), (31) follows. ∎

Let us now consider the tangent operator $\partial_{1}A(U,V;\mu)(Z)$, defined
as the Gateaux derivative of the operator $A(\cdot,\cdot;\mu)$, with respect
of the first variable, in the direction $Z\in X$. From (13) it holds

$\displaystyle\partial_{1}A(U,V;\mu)(Z)$ $\displaystyle=$
$\displaystyle\frac{1}{\mu}\,\partial_{1}A_{Diff}(U,V)(Z)+\partial_{1}A_{Div}(U,V)(Z)+\partial_{1}A_{Conv}(U,V)(Z)$
$\displaystyle+$
$\displaystyle\partial_{1}A_{Sma}(U;V)(Z)+\partial_{1}A_{Pres}(U,V;\mu)(Z),$

where

$\begin{array}[]{rl}\partial_{1}A_{Diff}(U,V)(Z)&=A_{Diff}(Z,V),\vspace{0.2cm}\\\
\partial_{1}A_{Div}(U,V)(Z)&=A_{Div}(Z,V),\vspace{0.2cm}\\\
\partial_{1}A_{Conv}(U;V)(Z)&=\displaystyle\frac{1}{2}\,\displaystyle\int_{\Omega}\left((\textbf{u}\cdot\nabla\textbf{z})\textbf{v}+(\textbf{z}\cdot\nabla\textbf{u})\textbf{v}-(\textbf{z}\cdot\nabla\textbf{v})\textbf{u}-(\textbf{u}\cdot\nabla\textbf{v})\textbf{z}\,\right)\;d\Omega,\vspace{0.2cm}\\\
\partial_{1}A_{Sma}(U;V)(Z)&=\displaystyle\int_{\Omega}\nu_{T}(\Pi_{h}^{*}\textbf{u})\,\nabla(\Pi_{h}^{*}\textbf{z}):\nabla(\Pi_{h}^{*}\textbf{v})\;d\Omega\vspace{0.1cm}\\\
&+\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}):\nabla(\Pi_{h}^{*}\textbf{z})}{|\nabla(\Pi_{h}^{*}\textbf{u})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}):\nabla(\Pi_{h}^{*}\textbf{v})\big{)}\;d\Omega,\vspace{0.2cm}\\\
\partial_{1}A_{Pres}(U,V;\mu)(Z)&=A_{Pres}(Z,V;\mu).\end{array}$

Note that if $\nabla(\Pi_{h}^{*}\textbf{u})=0$,
$\partial_{1}A_{Sma}(U;V)(Z)=0$, since $A_{Sma}(U;V)=0$. We further consider
that $\nabla(\Pi_{h}^{*}\textbf{u})\neq 0$. The same analysis can be
analogously performed when $\nabla(\Pi_{h}^{*}\textbf{u})=0$.

In the sequel we endow the product space $X=Y\times M$ with the hilbertian
norm

$\|U\|_{X}=\sqrt{\|\nabla\textbf{u}\|_{0,2,\Omega}^{2}+\|p\|_{0,2,\Omega}^{2}}\quad\forall\,U=(\textbf{u},p)\in
X.$ (32)

The derivation of the a posteriori error estimation is based upon the uniform
coerciveness and boundedness of the tangent operator
$\partial_{1}A(U,V;\mu)(Z)$ on the parametric solution of problem (4):

$\infty>\gamma_{0}\geq\gamma_{h}(\mu)\equiv\sup_{Z_{h}\in X_{h}}\sup_{V_{h}\in
X_{h}}\dfrac{\partial_{1}A(U_{h}(\mu),V_{h};\mu)(Z_{h})}{\|Z_{h}\|_{X}\|V_{h}\|_{X}}.$
(33) $0<\beta_{0}\leq\beta_{h}(\mu)\equiv\inf_{Z_{h}\in X_{h}}\sup_{V_{h}\in
X_{h}}\dfrac{\partial_{1}A(U_{h}(\mu),V_{h};\mu)(Z_{h})}{\|Z_{h}\|_{X}\|V_{h}\|_{X}}.$
(34)

for some constants $\gamma_{0}$, $\beta_{0}$. The coercivity (34) holds on
branches of non-singular solutions of problem (4), while the boundedness (33)
is a general property:

######  Proposition 3.

There exists $\gamma_{0}\in\mathbb{R}$ such that $\forall\mu\in\mathcal{D}$

$|\partial_{1}A(U_{h}(\mu),V_{h};\mu)(Z_{h})|\leq\gamma_{0}\|Z_{h}\|_{X}\|V_{h}\|_{X}\quad\forall
Z_{h},V_{h}\in X_{h}.$ (35)

Proof. Let us denote $U_{h}=(\textbf{u}_{h},p_{h}^{u})$,
$V_{h}=(\textbf{v}_{h},p_{h}^{v})$, $Z_{h}=(\textbf{z}_{h},p_{h}^{z})$. It
holds,

$|\partial_{1}A(U_{h}(\mu),V_{h};\mu)(Z_{h})|\leq\frac{1}{\mu}|\partial_{1}A_{Diff}(U_{h}(\mu),V_{h})(Z_{h})|+|\partial_{1}A_{Div}(U_{h}(\mu),V_{h})(Z_{h})|$
$+|\partial_{1}A_{Conv}(U_{h}(\mu),V_{h})(Z_{h})|+|\partial_{1}A_{Sma}(U_{h}(\mu);V_{h})(Z_{h})|+|\partial_{1}A_{Pres}(U_{h}(\mu),V_{h};\mu)(Z_{h})|.$

The boundedness of $|\partial_{1}A_{Diff}(U_{h}(\mu),V_{h};\mu)(Z_{h})|$,
$|\partial_{1}A_{Div}(U_{h}(\mu),V_{h})(Z_{h})|$, and
$|\partial_{1}A_{Conv}(U_{h}(\mu);V_{h})(Z_{h})|$ can be proved in a standard
way, the last one using estimate (23) (see [11]). To bound
$|\partial_{1}A_{Sma}(U_{h}(\mu);V_{h})(Z_{h})|$, we use the local inverse
finite element estimates and the boundedness of operator $\Pi_{h}$. It holds,

$\begin{array}[]{ll}|\partial_{1}A_{Sma}(U_{h}(\mu);V_{h})(Z_{h})|\leq&\displaystyle\int_{\Omega}(C_{S}h)^{2}|\nabla(\Pi_{h}^{*}\textbf{u}_{h})||\nabla(\Pi_{h}^{*}\textbf{z}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega,\vspace{0.1cm}\\\
&+\displaystyle\int_{\Omega}(C_{S}h)^{2}\dfrac{|\nabla(\Pi_{h}^{*}\textbf{u}_{h})||\nabla(\Pi_{h}^{*}\textbf{z}_{h})|}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h})|}|\nabla(\Pi_{h}^{*}\textbf{u}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega,\vspace{0.2cm}\\\
&\leq
2(C_{S}h)^{2}\|\nabla(\Pi_{h}^{*}\textbf{u}_{h})\|_{0,\infty,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}\vspace{0.1cm}\\\
&\leq
2C_{S}^{2}h^{2-d/2}\|\nabla(\Pi_{h}^{*}\textbf{u}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}\vspace{0.1cm}\\\
&\leq
2C_{1}C_{S}^{2}h^{2-d/2}\|\nabla(\Pi_{h}^{*}\textbf{u}_{h})\|_{0,2,\Omega}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega}\vspace{0.1cm}\\\
&\leq C_{3}\|Z_{h}\|_{X}\|V_{h}\|_{X},\end{array}$

for some constants $C_{1}>0$, $C_{2}>0$, $C_{3}>0$, where in the last
inequality we have used estimate (29). The last term can be bounded using
Lemma 3, and Cauchy-Schwartz inequality:

$\begin{array}[]{ll}|\partial_{1}A_{Pres}(U_{h}(\mu);V_{h};\mu)(Z_{h})|\leq&\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}\tau_{p,K}(\mu)|\sigma_{h}^{*}(\nabla
p_{h}^{z})||\sigma_{h}^{*}(\nabla p_{h}^{v})|\;d\Omega\vspace{0.1cm}\\\
&\leq\|\sigma_{h}^{*}(\nabla p_{h}^{z})\|_{\tau_{p}}\|\sigma_{h}^{*}(\nabla
p_{h}^{v})\|_{\tau_{p}}\leq
C_{\tau}^{2}\|p_{h}^{z}\|_{0,2,\Omega}\|p^{v}_{h}\|_{0,2,\Omega}\vspace{0.1cm}\\\
&\leq C_{4}\|Z_{h}\|_{X}\|V_{h}\|_{X}.\end{array}$

Thus, we have proved (35). ∎

The coercivity (34) may be proved for small enough data f:

######  Proposition 4.

Suppose that

$\|\textbf{f}\|_{Y^{\prime}}<\inf_{\mu\in\cal{D}}{\dfrac{1}{2\,\mu^{2}\,C_{P}\,C_{4S}^{2}}},$

where $C_{4S}$ is the norm of the injection of $H^{1}_{0}(\Omega)^{d}$ in
$L^{4}(\Omega)^{d}$. Then there exists $\tilde{\beta}>0$ such that,

$\tilde{\beta}\,\|Z_{h}\|_{X}\leq\sup_{V_{h}\in
X_{h}}\dfrac{\partial_{1}A(U_{h}(\mu),V_{h};\mu)(Z_{h})}{\|V_{h}\|_{X}}\quad\forall\,Z_{h}\in
X_{h},\quad\forall\,\mu\in{\cal D}.$ (36)

Proof. Let us denote, $U_{h}(\mu)=(\textbf{u}_{h},p_{h}^{u})$,
$Z_{h}=(\textbf{z}_{h},p_{h}^{z})$. It holds

$\partial_{1}A_{Pres}(U_{h}(\mu),Z_{h};\mu)(Z_{h})=\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}^{2}.$

Moreover,

$\partial_{1}A_{Sma}(U_{h}(\mu);Z_{h})(Z_{h})=\displaystyle\int_{\Omega}\nu_{T}(\nabla(\Pi_{h}^{*}\textbf{u}_{h}))|\nabla(\Pi_{h}^{*}\textbf{z}_{h})|^{2}\;d\Omega$
$+\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})|^{2}}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h})|}\;d\Omega\geq
0.$

Also,

$\partial_{1}A_{Diff}(U_{h}(\mu),Z_{h};\mu)(Z_{h})+\partial_{1}A_{Conv}(U_{h}(\mu),Z_{h})(Z_{h})$
$\geq\dfrac{1}{\mu}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2}+\partial_{1}A_{Conv}(U_{h}(\mu),Z_{h})(Z_{h})$
$\geq\displaystyle\dfrac{1}{\mu}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2}+\frac{1}{2}\,\displaystyle\int_{\Omega}\left((\textbf{z}_{h}\cdot\nabla\textbf{u}_{h})\textbf{z}_{h}-(\textbf{z}_{h}\cdot\nabla\textbf{z}_{h})\textbf{u}_{h}\,\right)\;d\Omega$
$\geq\frac{1}{\mu}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2}-\|\textbf{z}_{h}\|_{0,4,\Omega}^{2}\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}-\|\textbf{u}_{h}\|_{0,4,\Omega}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}\|\textbf{z}_{h}\|_{0,4,\Omega}$
$\geq\left(\frac{1}{\mu}-2\,C_{S4}^{2}\,\,\|\nabla\textbf{u}_{h}\|_{0,2,\Omega}\right)\,\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2}\geq\left(\frac{1}{\mu}-2\,C_{S4}^{2}\,\mu\,C_{P}\,\|\textbf{f}\|_{Y^{\prime}}\right)\,\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2}$
$\geq\rho\,\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2},\quad\forall\,\mu\in\mathcal{D},\hfill$

for
$\rho=\displaystyle\inf_{\mu\in\mathcal{D}}\left(\frac{1}{\mu}-2\,C_{S4}^{2}\,\mu\,C_{P}\,\|\textbf{f}\|_{Y^{\prime}}\right)>0$.
Consequently,

$\partial_{1}A(U_{h}(\mu),Z_{h};\mu)(Z_{h})\geq\rho\,\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}^{2}+\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}^{2}.$

Let

$S(\mu)=\displaystyle\sup_{V_{h}\in
X_{h}}\dfrac{\partial_{1}A(U_{h}(\mu),V_{h};\mu)(Z_{h})}{\|V_{h}\|_{X}}.$

Then, in abridged notation for brevity,

$\rho\,\|\textbf{z}_{h}\|_{Y}^{2}+\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}^{2}\leq
S(\mu)\,\|Z_{h}\|_{X}=S(\mu)\,\left(\|\textbf{z}_{h}\|_{Y}^{2}+\|p_{h}^{z}\|_{M}^{2}\right)^{1/2}.$
(37)

Also, letting $V_{h}=(\textbf{v}_{h},0)\in X_{h}$, estimates similar to those
carried on in the proof of Proposition 3 yield

$(\nabla\cdot\textbf{v}_{h},p_{h}^{z})\leq\left(S(\mu)+(\frac{1}{\mu}+C)\,\|\textbf{z}_{h}\|_{Y}\right)\,\|\textbf{v}_{h}\|_{Y}$

for some constant $C>0$ independent of $\nu$ and $h$. Then, letting
$\displaystyle z=\sup_{\textbf{v}_{h}\in
Y_{h}}\dfrac{(p_{h}^{z},\nabla\cdot\textbf{v}_{h})_{\Omega}}{\|\textbf{v}_{h}\|_{Y}}$,

$z\leq S(\mu)+C^{\prime}(\mu)\,\|\textbf{z}_{h}\|_{Y}\,\quad\mbox{with
}C^{\prime}(\mu)=\frac{1}{\mu}+C.$ (38)

Combining (37) and (38), and using (28),

$\displaystyle\rho\,\|\textbf{z}_{h}\|_{Y}^{2}+\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}^{2}\leq
S(\mu)\,\left(\|\textbf{z}_{h}\|_{Y}+z+\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}\right)$ $\displaystyle\leq
S(\mu)\,\left(\alpha\,S(\mu)+(1+\alpha\,C^{\prime}(\mu))\,\|\textbf{z}_{h}\|_{Y}^{2}+\alpha\,\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}\right)$

Using Young’s inequality, it follows

$\displaystyle\rho\,\|\textbf{z}_{h}\|_{Y}^{2}+\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}^{2}$ $\displaystyle\leq$ $\displaystyle
C\,\left(1+\frac{1+C^{\prime}(\mu)}{\rho}\right)\,S(\mu)^{2}$ (39)
$\displaystyle\leq$ $\displaystyle
C\,\left(1+\frac{1}{\mu\,\rho}\right)\,S(\mu)^{2}.$

Combining this inequality with (38),

$z\leq\left(1+C\,C^{\prime}(\mu)\,\sqrt{1+\frac{1}{\mu\,\rho}}\right)\,S(\mu).$
(40)

Taking the square rooth of (39) and summing it with (40) finally yields

$\displaystyle\displaystyle\min\\{\sqrt{\rho},\,\alpha^{-1}\\}\,\|Z_{h}\|_{X}\leq\sqrt{\rho}\,\|\textbf{z}_{h}\|_{Y}+z+\|\sigma_{h}^{*}(\nabla
p_{h}^{z})\|_{\tau_{p}}$ $\displaystyle\leq\displaystyle
C\,\left(1+\frac{1}{\mu}\right)\,\sqrt{1+\frac{1}{\mu\,\rho}}\,S(\mu)\leq
M\,S(\mu),\,\,$

with
$M=\displaystyle\sup\left\\{C\,\left(1+\frac{1}{\mu}\right)\,\sqrt{1+\frac{1}{\mu\,\rho}},\,\,\,\mu\in\mathcal{D}\right\\}$.
Then, (42) holds with $\tilde{\beta}=\min\\{\sqrt{\rho},\,\alpha^{-1}\\}/M$. ∎

The same demonstration proves the coercivity of the reduced problem for small
data, when the reduced space $X_{N}$ is endowed with the norm

$\|U_{N}\|_{N}=\sqrt{\|\nabla\textbf{u}_{N}\|_{0,2,\Omega}^{2}+|||p_{N}|||_{N}^{2}}\quad\forall\,U_{N}=(\textbf{u}_{N},p_{N})\in
X_{N},$ (41)

where $|||\cdot|||_{N}$ is defined in (19).

######  Proposition 5.

Assume that the pressure finite element space $M_{h}$ is formed by piecewise
affine functions. Suppose that

$\|\textbf{f}\|_{Y^{\prime}}\leq\inf_{\mu\in\cal{D}}{\dfrac{1}{2\,\mu^{2}\,C_{P}\,C_{4S}^{2}}}.$

Then there exists a constant $\tilde{\beta}>0$ such that,

$\tilde{\beta}\,\|Z_{N}\|_{N}\leq\sup_{V_{N}\in
X_{N}}\dfrac{\partial_{1}A(U_{N}(\mu),V_{N};\mu)(Z_{N})}{\|V_{N}\|_{N}}\quad\forall\,Z_{N}\in
X_{N},\quad\forall\,\mu\in{\cal D}.$ (42)

In this case, the $\|\cdot\|_{N}$ norm is bounded from above by a constant
times the $\|\cdot\|_{0,2,\Omega}$ norm but, to the best of our knowledge,
they are no longer equivalent.

##  4 A posteriori error bound estimator

In this section we develop the a posteriori error bound estimator used in the
hierarchical construction of the reduced space during the Greedy algorithm. We
use the Brezzi-Rappaz-Raviart (BRR) theory [7] for approximation on a regular
branch of solutions of non-linear problems. We start by proving the local
Lipschitz continuity of the tangent Smagorinsky operator, which is a key
result to apply this theory:

######  Lemma 4.

It holds

$\displaystyle\left|\partial_{1}A(U_{h}^{1},V_{h};\mu)(Z_{h})-\partial_{1}A(U_{h}^{2},V_{h};\mu)(Z_{h})\right|\leq$
(43)
$\displaystyle\rho_{T}\|U_{h}^{1}-U_{h}^{2}\|_{X}\|Z_{h}\|_{X}\|V_{h}\|_{X},$

for all $U_{h}^{1},U_{h}^{2},Z_{h},V_{h}\in X_{h}$, with
$\rho_{T}=2\,C_{S4}^{2}+C\,h^{2-d/2}$, for some constant $C>0$.

Proof. We have

$\displaystyle\partial_{1}A(U_{h}^{1},V_{h};\mu)(Z_{h})-\partial_{1}A(U_{h}^{2},V_{h};\mu)(Z_{h})=$
$\displaystyle\displaystyle\frac{1}{2}\left(\displaystyle\int_{\Omega}((\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\cdot\nabla\textbf{z}_{h})\textbf{v}_{h}\;d\Omega+\displaystyle\int_{\Omega}(\textbf{z}_{h}\cdot\nabla(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2}))\textbf{v}_{h}\;d\Omega\right.$
$\displaystyle\left.-\displaystyle\int_{\Omega}(\textbf{z}_{h}\cdot\nabla\textbf{v}_{h})(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\;d\Omega-\displaystyle\int_{\Omega}((\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\cdot\nabla\textbf{v}_{h})\textbf{z}_{h}\;d\Omega\right)$
$\displaystyle+\displaystyle\int_{\Omega}\left[\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}^{1})-\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}^{2})\right]\nabla(\Pi_{h}^{*}\textbf{z}_{h}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\;d\Omega$
$\displaystyle+\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\;d\Omega$
$\displaystyle-\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\;d\Omega.$

Thanks to the triangular inequality, it holds

$\begin{array}[]{lc}&\left|\partial_{1}A(U_{h}^{1},V_{h};\mu)(Z_{h})-\partial_{1}A(U_{h}^{2},V_{h};\mu)(Z_{h})\right|\leq\\\
&\displaystyle\frac{1}{2}\left|\displaystyle\int_{\Omega}((\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\cdot\nabla\textbf{z}_{h})\textbf{v}_{h}\;d\Omega\right|+\displaystyle\frac{1}{2}\left|\displaystyle\int_{\Omega}(\textbf{z}_{h}\cdot\nabla(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2}))\textbf{v}_{h}\;d\Omega\right|\\\
&+\displaystyle\frac{1}{2}\left|\displaystyle\int_{\Omega}(\textbf{z}_{h}\cdot\nabla\textbf{v}_{h})(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\;d\Omega\right|+\displaystyle\frac{1}{2}\left|\displaystyle\int_{\Omega}((\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\cdot\nabla\textbf{v}_{h})\textbf{z}_{h}\;d\Omega\right|\\\
&+\left|\displaystyle\int_{\Omega}\left[\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}^{1})-\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}^{2})\right]\nabla(\Pi_{h}^{*}\textbf{z}_{h}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\;d\Omega\right|\\\
&+\left|\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\;d\Omega\right.\\\
&\left.-\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\;d\Omega\right|.\\\
\end{array}$ (44)

We bound each term in (44) separately. The four first terms can be bounded in
a standard way using estimate (23). For the third term, we use the local
inverse inequalities of finite element functions (cf. [5]) and the stability
of operator $\Pi_{h}$,

$\left|\displaystyle\int_{\Omega}\left(\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}^{1})-\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}^{2})\right)\nabla(\Pi_{h}^{*}\textbf{z}_{h}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\;d\Omega\right|$
$\leq\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\big{|}|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})|-|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2})|\big{|}|\nabla(\Pi_{h}^{*}\textbf{z}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega$
$\leq(C_{S}h)^{2}\displaystyle\int_{\Omega}|\nabla\big{(}\Pi_{h}^{*}(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\big{)}||\nabla(\Pi_{h}^{*}\textbf{z}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega$
$\leq(C_{S}h)^{2}\|\nabla\big{(}\Pi_{h}^{*}(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\big{)}\|_{0,\infty,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}$
$\leq
C_{S}^{2}h^{2-d/2}C\|\nabla\big{(}\Pi_{h}^{*}(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\big{)}\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}$
$\leq
C_{S}^{2}h^{2-d/2}CC_{f}\,\|\nabla(\textbf{u}_{h}^{1}-\textbf{u}_{h}^{2})\|_{0,2,\Omega}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega}$
$\leq
C_{S}^{2}h^{2-d/2}CC_{f}\,\|U_{h}^{1}-U_{h}^{2}\|_{X}\|Z_{h}\|_{X}\|V_{h}\|_{X}$

The last term in (44) is bounded as follows:

$\left|\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\;d\Omega\right.$
$\left.-\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\dfrac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\;d\Omega\right|$
$=\left|\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\left[\frac{\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})|}\left(\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h}-\Pi_{h}^{*}\textbf{u}^{2}_{h}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\right)\right.\right.$
$+\frac{\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h}-\Pi_{h}^{*}\textbf{u}^{2}_{h}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}$
$+\left.\left.\frac{(|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2})|-|\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h})|)\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1}):\nabla(\Pi_{h}^{*}\textbf{z}_{h})}{|\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})||\nabla(\Pi_{h}^{*}\textbf{u}^{2}_{h})|}\big{(}\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{2}):\nabla(\Pi_{h}^{*}\textbf{v}_{h})\big{)}\right]\;d\Omega\right|$
$\leq\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}|\nabla(\Pi_{h}^{*}\textbf{z}_{h})||\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h}-\Pi_{h}^{*}\textbf{u}^{2}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega$
$+\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}|\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h}-\Pi_{h}^{*}\textbf{u}^{2}_{h})||\nabla(\Pi_{h}^{*}\textbf{z}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega$
$+\displaystyle\sum_{K\in\mathcal{T}_{h}}\displaystyle\int_{K}(C_{S}h_{K})^{2}\left||\nabla(\Pi_{h}^{*}\textbf{u}_{h}^{1})|-|\nabla(\Pi_{h}^{*}\textbf{u}^{2}_{h})|\right||\nabla(\Pi_{h}^{*}\textbf{z}_{h})||\nabla(\Pi_{h}^{*}\textbf{v}_{h})|\;d\Omega$
$\leq
3(C_{S}h)^{2}\|\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h}-\Pi_{h}^{*}\textbf{u}^{2}_{h})\|_{0,\infty,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}$
$\leq
3C_{S}^{2}h^{2-d/2}C\|\nabla(\Pi_{h}^{*}\textbf{u}^{1}_{h}-\Pi_{h}^{*}\textbf{u}^{2}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{z}_{h})\|_{0,2,\Omega}\|\nabla(\Pi_{h}^{*}\textbf{v}_{h})\|_{0,2,\Omega}$
$\leq
3C_{S}^{2}h^{2-d/2}CC_{f}\,\|\nabla(\textbf{u}^{1}_{h}-\textbf{u}^{2}_{h})\|_{0,2,\Omega}\|\nabla\textbf{z}_{h}\|_{0,2,\Omega}\|\nabla\textbf{v}_{h}\|_{0,2,\Omega}$
$\leq
3C_{S}^{2}h^{2-d/2}CC_{f}\,\|U^{1}_{h}-U^{2}_{h}\|_{X}\|Z_{h}\|_{X}\|V_{h}\|_{X}.$

Thus, we have just proved that

$\left|\partial_{1}A(U_{h}^{1},V_{h};\mu)(Z_{h})-\partial_{1}A(U_{h}^{2},V_{h};\mu)(Z_{h})\right|\leq\rho_{T}\|U^{1}_{h}-U^{2}_{h}\|_{X}\|Z_{h}\|_{X}\|V_{h}\|_{X},$

with, $\rho_{T}=2C_{T}+C\,h^{2-d/2}$ for some constant $C>0$. ∎

We assume that the RB problem (14) is well posed, and define the following
inf-sup and continuity constants:

$0<\beta_{N}(\mu)\equiv\inf_{Z_{h}\in X_{h}}\sup_{V_{h}\in
X_{h}}\dfrac{\partial_{1}A(U_{N}(\mu),V_{h};\mu)(Z_{h})}{\|Z_{h}\|_{X}\|V_{h}\|_{X}},$
(45) $\infty>\gamma_{N}(\mu)\equiv\sup_{Z_{h}\in X_{h}}\sup_{V_{h}\in
X_{h}}\dfrac{\partial_{1}A(U_{N}(\mu),V_{h};\mu)(Z_{h})}{\|Z_{h}\|_{X}\|V_{h}\|_{X}},$
(46)

Let us introduce the a posteriori error bound estimator:

$\Delta_{N}(\mu)=\frac{\beta_{N}(\mu)}{2\rho_{T}}\left[1-\sqrt{1-\tau_{N}(\mu)}\right],$
(47)

where $\tau_{N}(\mu)$ is given by

$\displaystyle\tau_{N}(\mu)$
$\displaystyle=\frac{4\epsilon_{N}(\mu)\rho_{T}}{\beta_{N}^{2}(\mu)},$ (48)

with $\epsilon_{N}(\mu)$ the dual norm, on $X_{h}^{\prime}$, of the residual
${\cal R}(U_{N}(\mu);\mu)$, defined as

${\cal R}(U_{N}(\mu);\mu)(V_{h})=A(U_{N}(\mu),V_{h};\mu)-F(V_{h};\mu);$

that is,

$\epsilon_{N}(\mu)=\sup_{V_{h}\in X_{h}}\frac{{\cal
R}(U_{N}(\mu);\mu)(V_{h})}{\|V_{h}\|_{X}}.$

The a posteriori error bound estimator is stated by the following result

######  Theorem 2.

Let $\mu\in\mathcal{D}$, and assume that $\beta_{N}(\mu)>0$. If problem (12)
admits a solution $U_{h}(\mu)$ such that

$\|U_{h}(\mu)-U_{N}(\mu)\|_{X}\leq\frac{\beta_{N}(\mu)}{\rho_{T}},$

then this solution is unique in the ball
$B_{X}\left(U_{N}(\mu),\dfrac{\beta_{N}(\mu)}{\rho_{T}}\right)$.

Moreover, assume that $\tau_{N}(\mu)\leq 1$ for all $\mu\in\mathcal{D}$. Then
there exists a unique solution $U_{h}(\mu)$ of (12) such that the error with
respect to $U_{N}(\mu)$, solution of (14), is bounded by

$\|U_{h}(\mu)-U_{N}(\mu)\|_{X}\leq\Delta_{N}(\mu)\leq\left[\frac{2\gamma_{N}(\mu)}{\beta_{N}(\mu)}+\tau_{N}(\mu)\right]\|U_{h}(\mu)-U_{N}(\mu)\|_{X}.$
(49)

The proof of this theorem follows that of Theorem 5.3 in [11], we omit it for
brevity.

######  Remark 1.

Note that $\Delta_{N}(\mu)$ is defined only if $\tau_{N}(\mu)\leq 1$. In
practice, at the first iterations of the Greedy algorithm, $\tau_{N}(\mu)>1$.
For these cases, we consider a Taylor approximation of the a posteriori error
estimator, that is given by

$\Delta_{N}(\mu)\simeq\dfrac{\beta_{N}(\mu)}{4\rho_{T}}\tau_{N}(\mu).$ (50)

In Section 6 we show that this Taylor approximation is more accurate than
$\tau_{N}(\mu)$ as estimator.

##  5 Approximation of eddy viscosity term and pressure stabilizing
coefficient

In this section, we present the approximation of the non-linear terms with
respect to the parameter, throughout the Empirical Interpolation Method [18,
3].

Both the Small-Small setting of the Smagorinsky eddy-diffusion term defined in
(13), $\nu_{T}(\Pi_{h}^{*}(\nabla\textbf{u})):=\nu_{T}(\mu)$, and the pressure
stabilisation coefficient, $\tau_{p,K}(\mu)$, defined in (9) are non-linear
functions of the parameter, and consequently need to be linearised with the
EIM to build the RB model to reduce the on-line computation times.

For this purpose, we build two reduced-basis spaces
$W_{M_{1}}^{S}=\\{q_{1}^{S}(\mu),\dots,q_{M_{1}}^{S}(\mu)\\}$ and
$W_{M_{2}}^{P}=\\{q_{1}^{P}(\mu),\dots,q_{M_{2}}^{P}(\mu)\\}$ by a greedy
procedure selection, with $W_{M_{1}}^{S}$ and $W_{M_{2}}^{P}$ the EIM reduced
spaces associated to the eddy-viscosity term and the pressure-stabilisation
term, respectively.

Thus, we approximate them by the following trilinear forms:

$\begin{array}[]{rl}a_{S}^{\prime}(\textbf{u}_{N};\textbf{u}_{N},\textbf{v}_{N};\mu)&\approx\hat{a}_{S}^{\prime}(\textbf{u}_{N},\textbf{v}_{N};\mu),\vspace{0.1cm}\\\
s_{pres}(p_{N},q_{N};\mu)&\approx\hat{s}_{pres}(p_{N},q_{N};\mu),\end{array}$
(51)

where,

$\begin{array}[]{rl}\hat{a}_{S}^{\prime}(\textbf{u}_{N},\textbf{v}_{N};\mu)&=\displaystyle\sum_{k=1}^{M_{1}}\sigma^{S}_{k}(\mu)s(q_{k}^{S},\textbf{u}_{N},\textbf{v}_{N}),\vspace{0.1cm}\\\
\hat{s}_{pres}(p_{N},q_{N};\mu)&=\displaystyle\sum_{k=1}^{M_{2}}\sigma^{P}_{k}(\mu)r(q_{k}^{P},p_{N},q_{N}),\end{array}$
(52)

with,

$\begin{array}[]{rl}s(q_{k}^{S},\textbf{u}_{N},\textbf{v}_{N})&=\displaystyle\sum_{K\in\mathcal{T}_{h}}\big{(}q^{S}_{k}\,\nabla(\Pi_{h}^{*}\textbf{u}_{N}),\nabla(\Pi_{h}^{*}\textbf{v}_{N})\big{)}_{K},\vspace{0.1cm}\\\
r(q_{k}^{P},p_{N},q_{N})&=\displaystyle\sum_{K\in\mathcal{T}_{h}}\big{(}q^{P}_{k}\,\sigma_{h}^{*}(\nabla
p_{N}),\sigma_{h}^{*}(\nabla q_{N})\big{)}_{K}.\end{array}$ (53)

Here we are considering that the approximations given by the EIM for
$\nu_{T}(\mu)$ and $\tau_{K,p}(\mu)$ are respectively

$\begin{array}[]{rl}\mathcal{I}_{M_{1}}[\nu_{T}(\mu)]&=\displaystyle\sum_{k=1}^{M_{1}}\sigma^{S}_{k}(\mu)q_{k}^{S},\vspace{0.1cm}\\\
\mathcal{I}_{M_{2}}[\tau_{K,p}(\mu)]&=\displaystyle\sum_{k=1}^{M_{2}}\sigma^{P}_{k}(\mu)q_{k}^{P}.\end{array}$
(54)

In practice, to deal with non-homogeneous Dirichlet boundary conditions, we
split the solution as $\textbf{u}_{N}+\textbf{u}_{D}$, where
$\textbf{u}_{N}\in X_{h}$ and $\textbf{u}_{D}\in H^{1}(\Omega)^{d}$ is a lift
of the Dirichlet boundary condition, compatible with the remaining boundary
conditions (if any). We are thus led to a “practical” RB problem, given by

$\left\\{\begin{array}[]{l}\mbox{Find
}(\textbf{u}_{N},p_{N})=(\textbf{u}_{N}(\mu),p_{N}(\mu))\in Y_{N}\times
M_{N}\mbox{ such that}\vspace{0.3cm}\\\
\begin{array}[]{ll}a(\textbf{u}_{N},\textbf{v}_{N};\mu)+b(\textbf{v}_{N},p_{N};\mu)+\hat{a}_{S}^{\prime}(\textbf{u}_{N};\textbf{v}_{N};\mu)\\\
+c(\textbf{u}_{N},\textbf{u}_{N},\textbf{v}_{N};\mu)+c(\textbf{u}_{D},\textbf{u}_{N},\textbf{v}_{N};\mu)\\\
+c(\textbf{u}_{N},\textbf{u}_{D},\textbf{v}_{N};\mu)=F(\textbf{v}_{N};\mu)&\quad\forall\textbf{v}_{N}\in
Y_{N}\vspace{0.1cm}\\\
b(\textbf{u}_{N},q_{N};\mu)+\hat{s}_{pres}(p_{N},q_{N};\mu)=0&\quad\forall
q_{N}\in M_{N}.\end{array}\end{array}\right.$ (55)

We can express the solution $(\textbf{u}_{N}(\mu),p_{N}(\mu))\in X_{N}$ of
(55) as a linear combination of the basis functions:

$\textbf{u}_{N}(\mu)=\sum_{j=1}^{N}u_{j}^{N}(\mu)\zeta^{\textbf{v}}_{j},\quad
p_{N}(\mu)=\sum_{j=1}^{N}p_{j}^{N}(\mu)\xi^{p}_{j}.$

The tensor representation associated to the eddy-viscosity term and the
pressure stabilisation coefficient are define as follow:

$(\mathbb{S}^{\prime}_{N}(q^{S}_{s}))_{ij}=\sum_{K\in\mathcal{T}_{h}}(q^{S}_{s}\nabla(\Pi_{h}^{*}\zeta^{\textbf{v}}_{j}),\nabla(\Pi_{h}^{*}\zeta^{\textbf{v}}_{i}))_{K},\quad
i,j=1\dots,N,s=1,\dots,M_{1},$
$(\mathbb{P}_{N}(q^{P}_{s}))_{ij}=\sum_{K\in\mathcal{T}_{h}}(q^{P}_{s}\sigma_{h}^{*}(\nabla\xi^{p}_{j}),\sigma_{h}^{*}(\nabla\xi^{p}_{i}))_{K},\quad
i,j=1,\dots,N,s=1,\dots,M_{2}.$

With this tensor representation for the non-linear terms in (14) done in the
offline phase, it holds that

$\hat{a}_{S}^{\prime}(\zeta^{\textbf{v}}_{j};\zeta^{\textbf{v}}_{i};\mu)=\sum_{s=1}^{M_{1}}\sigma_{s}^{S}(\mu)\mathbb{S}^{\prime}_{N}(q^{S}_{s})\quad\mbox{
and
}\quad\hat{s}_{pres}(\xi^{p}_{j},\xi^{p}_{i};\mu)=\sum_{s=1}^{M_{2}}\sigma_{s}^{P}(\mu)\mathbb{P}_{N}(q^{P}_{s}).$

Problem (55) is solved by a semi-implicit evolution approach. Thanks to these
linearised representations, we are able to efficiently solve the RB model in
the online phase.

##  6 Numerical Results

In this section, we present some numerical results to test the LPS pressure
stabilisation for the Smagorinsky reduced basis model. We consider the lid-
driven cavity problem for two different finite element discretisations. To let
this problem fit into the preceding theory we consider a lifting
$\textbf{u}_{D}\in[H^{1}(\Omega)]^{2}$ of the non-homogeneous boundary
conditions, and search for the solution as
$\textbf{u}=\textbf{u}_{0}+\textbf{u}_{D}$, where the new unknown
$\textbf{u}_{0}$ satisfies homogeneous Dirichlet boundary conditions. The
change in the convection term makes it appear a linear term with the structure
$\textbf{u}_{D}\nabla\cdot\textbf{u}_{0}+\textbf{u}_{0}\nabla\cdot\textbf{u}_{D}$,
the remaining operator terms in the Smagorinsky model have the same structure.
An involved, but straightforward, extension of the preceding analysis, proves
that it also applies to this situation whenever $\textbf{u}_{D}$ is small
enough with respect to the viscosity coefficient. We have preferred not to
include all this analysis for brevity. In our tests we use as
$\textbf{u}_{D}\in[H^{1}(\Omega)]^{2}$ the standard boundary condition on the
cavity lid $\textbf{u}=(1,0)$, extended by zero to the remaining of
$\partial\Omega$.

We include LPS stabilisation of the pressure for both the full order and the
reduced methods. On one hand the Taylor-Hood ($\mathbb{P}_{2}-\mathbb{P}_{1}$)
finite element, that fits into the theory developed above, then both the full
order finite element and the reduced methods are stable. On another hand, the
$\mathbb{P}_{2}-\mathbb{P}_{2}$ finite element for velocity-pressure, for
which the full order method is stable, but -to the best of our knowledge- not
the reduced basis problem without pressure supremizers. For both cases, we
consider the two strategies of either including pressure supremizers, or not,
to build the reduced pressure spaces, to compare the performances of both
strategies. We use FreeFem++ to perform the computations (cf. [19]).

We consider the Reynolds number as a parameter, ranging in
$\mathcal{D}=[1000,5100]$. In this range of Reynolds number, the solution
reaches a steady state regime. We consider a regular mesh with 5000 triangles
and 2601 nodes. In both cases, we implement the EIM for the eddy viscosity
term and the pressure stabilisation constant.

We are aware that in the turbulent regime the lid-driven cavity flow is 3D
(Cf. [14], Chapter 13). We consider these 2D tests to assess the theoretical
derivations within the paper, and to compare the two strategies of including -
not including pressure supremizers in the reduced velocity spaces, when the
reduced model includes the LPS stabilisation of the pressure.

For the numerical tests, we consider a weighted the scalar product for the
velocity space, that lets us improve the efficiency of the a posteriori error
estimations as shown in [11]. This norm is defined as

$(\textbf{u}_{h},\textbf{v}_{h})_{T}=\displaystyle\int_{\Omega}\left[\frac{1}{\overline{\mu}}+\nu_{T}^{\prime*}\right]\nabla\textbf{u}_{h}:\nabla\textbf{v}_{h}\;d\Omega\qquad\forall\textbf{u}_{h},\textbf{v}_{h}\in
Y_{h},$ (56)

where $\nu_{T}^{\prime*}=\nu_{T}(\Pi_{h}^{*}\textbf{u}_{h}(\overline{\mu}))$,
and

$\overline{\mu}=\arg\displaystyle\min_{\mu\in\mathcal{D}}\displaystyle\sum_{K\in\mathcal{T}_{h}}(C_{S}h_{K})^{2}\displaystyle\min_{x\in
K}|\nabla(\Pi_{h}^{*}\textbf{u}_{h})(\mu)|(x)\chi_{K}(x),$

with $\textbf{u}_{h}(\mu)$ the velocity solution of (4)

###  6.1 Test 1. Taylor-Hood finite element

For this test, we need $M_{1}=49$ basis functions $q_{k}^{S}$ for the eddy
viscosity $\nu_{T}$ EIM approximation and $M_{2}=29$ basis functions
$q_{l}^{S}$ for the stabilisation pressure coefficient $\tau_{p,K}$ EIM
approximation to reach a prescribed tolerance of $\varepsilon_{EIM}=5\cdot
10^{-5}$ for both terms. In Fig. 1 we show the evolution of the error for the
both EIM approximations.

Figure 1: Test 1. Evolution of the parametric $L^{\infty}$ error in the EIM
approximation of eddy viscosity and stabilisation pressure coefficient.

For the Greedy algorithm that selects the reduced basis functions for velocity
and pressure, we prescribe a tolerance of $\varepsilon_{RB}=7\cdot 10^{-5}$,
which is reached for $N=16$ basis functions in both cases. In Fig. 2 we show
the evolution of the a posteriori error bound estimator during the Greedy
algorithm for both strategies (including pressure supremizers or not in the
reduced velocity space). We can observe that the evolution of the maximum
value for the error estimator either considering supremizers (left) or not
(right) are equivalent, giving the same number of basis functions until
reaching the tolerance.

Figure 2: Test 1. Evolution of the error in the Greedy algorithm with (left)
and without pressure supremizers (right).

We also have measured, for some randomly selected parameter values, the
maximum exact error for velocity and pressure snapshot solutions with respect
to the reduced solution computed, during the Greedy algorithm. We can observe
in Fig. 3 that for both the velocity and pressure the maximum error values in
each iteration of the Greedy algorithm, for both velocity and pressure, are
quite close.

Figure 3: Test 1. Velocity (left) and pressure (right) maximum error during
the Greedy algorithm.

In Fig. 4, we show the comparison of the exact error for some snapshots and
the a posteriori error bound estimator considering either including the
supremizers in the velocity spaces (left) or not (right). We can observe that
the effectivity in this case is about one order of magnitude. The estimator is
somewhat less smooth as function of the parameter $\mu$ for the LPS-VMS
method.

Figure 4: Error and a posteriori error bound estimator for $N=16$, with (left)
and without pressure supremizers (right).

Finally, in Table 1 we summarize results obtained for several values of
Reynold numbers, where we can observe that the absolute error between FE and
RB solution is $10^{-6}$ for the velocity and $10^{-8}$ for pressure.
Moreover, we can observe that the speedup of the solution is several thousands
of magnitude, with higher speedups (increase of nearly 50%) and somewhat
larger errors for the strategy of not including the pressure supremizers in
the reduced pressure spaces.

Without supremizer
---
Data | $\mu=1610$ | $\mu=2751$ | $\mu=3886$ | $\mu=4521$
$T_{FE}$ | 1372.65s | 1513.53s | 3184.17s | 3608.32s
$T_{online}$ | 0.57s | 0.65s | 0.7s | 0.77s
speedup | 2407 | 2327 | 4548 | 4685
$\|\textbf{u}_{h}-\textbf{u}_{N}\|_{T}$ | $2.01\cdot 10^{-6}$ | $1.98\cdot 10^{-6}$ | $2.04\cdot 10^{-6}$ | $2.64\cdot 10^{-6}$
$\|p_{h}-p_{N}\|_{0,2,\Omega}$ | $1.51\cdot 10^{-8}$ | $2.03\cdot 10^{-8}$ | $9.73\cdot 10^{-9}$ | $2.05\cdot 10^{-8}$

With supremizer
---
Data | $\mu=1610$ | $\mu=2751$ | $\mu=3886$ | $\mu=4521$
$T_{FE}$ | 1372.65s | 1513.53s | 3184.17s | 3608.32s
$T_{online}$ | 0.85s | 0.97s | 1.08s | 1.11s
speedup | 1596 | 1544 | 2947 | 3224
$\|\textbf{u}_{h}-\textbf{u}_{N}\|_{T}$ | $2.21\cdot 10^{-6}$ | $2.85\cdot 10^{-6}$ | $3.83\cdot 10^{-7}$ | $1.98\cdot 10^{-6}$
$\|p_{h}-p_{N}\|_{0,2,\Omega}$ | $1.09\cdot 10^{-8}$ | $1.37\cdot 10^{-8}$ | $7.08\cdot 10^{-9}$ | $1.10\cdot 10^{-8}$

Table 1: Computational time for FE solution and RB online phase, with the
speedup and the error.

###  6.2 Test 2: $\mathbb{P}_{2}-\mathbb{P}_{2}$ velocity-pressure finite
element

For this test, we need $M_{1}=52$ basis for the eddy viscosity approximation
and $M_{2}=28$ basis for the pressure stabilisation constant approximation to
reach a prescribe tolerance of $\varepsilon_{EIM}=5\cdot 10^{-5}$ for both
terms. In Fig. 5 we show the evolution of the error for both non-linear terms.

Figure 5: Evolution of the error in the EIM.

For the Greedy algorithm to determine the reduced velocity and pressure
spaces, in this case we prescribe the same tolerance of
$\varepsilon_{RB}=7\cdot 10^{-5}$. This tolerance is reached for $N=16$ basis
functions when we do not consider supremizers, while in the case of
considering the supremizer, the tolerance is reached for $N=17$. In Fig. 6 we
show the evolution of the a posteriori error bound estimator during the Greedy
algorithm for both cases. We can observe, again, that the evolution of the
maximum value for the error estimator either considering supremizers (left) or
not (right) are equivalent. Although we do not know wether the reduced problem
without supremizers is stable, all computations take place as if it was.

Figure 6: Test 2. Evolution of the error in the Greedy algorithm with (left)
and without pressure supremizers (right).

We can observe in Fig. 7 that for both the velocity and pressure the maximum
error value in each iteration of the Greedy algorithm, for both velocity and
pressure, are quite close, as in the Taylor-Hood finite element test.

Figure 7: Test 2. Velocity (left) and pressure (right) maximum error during
the Greedy algorithm.

In Fig. 8 we show the comparison of the exact error for some snapshots and the
a posteriori error bound estimator considering either including the
supremizers in the velocity spaces (left) or not (right).

Figure 8: Test 2. Error and a posteriori error bound estimator for $N=16$,
with (left) and without pressure supremizers (right).

Finally, in Table 2 we summarize results obtained for several values of
Reynold numbers, where we can observe that the absolute error between FE and
RB solution is around order of $10^{-6}$ for the velocity and $10^{-8}$ for
pressure. Moreover, we can observe that the speedup of the solution is of
several thousands of magnitude, with increase of nearly 50%, and somewhat
larger errors, when no supremizers are used.

Without supremizer
---
Data | $\mu=1610$ | $\mu=2751$ | $\mu=3886$ | $\mu=4521$
$T_{FE}$ | 2259.04s | 3008.81s | 5574.5s | 6171.41s
$T_{online}$ | 0.70s | 0.75s | 0.78s | 0.81s
speedup | 3227 | 4011 | 7146 | 7619
$\|\textbf{u}_{h}-\textbf{u}_{N}\|_{T}$ | $5.4\cdot 10^{-7}$ | $1.44\cdot 10^{-6}$ | $1.49\cdot 10^{-6}$ | $5.44\cdot 10^{-6}$
$\|p_{h}-p_{N}\|_{0,2,\Omega}$ | $1.34\cdot 10^{-8}$ | $3.35\cdot 10^{-8}$ | $1.55\cdot 10^{-7}$ | $2.25\cdot 10^{-8}$

With supremizer
---
Data | $\mu=1610$ | $\mu=2751$ | $\mu=3886$ | $\mu=4521$
$T_{FE}$ | 2259.04s | 3008.81s | 5574.5s | 6171.41s
$T_{online}$ | 0.99s | 1.03s | 1.26s | 1.23s
speedup | 2267 | 2895 | 4391 | 5016
$\|\textbf{u}_{h}-\textbf{u}_{N}\|_{T}$ | $4.94\cdot 10^{-7}$ | $4.51\cdot 10^{-6}$ | $6.52\cdot 10^{-7}$ | $3.52\cdot 10^{-6}$
$\|p_{h}-p_{N}\|_{0,2,\Omega}$ | $5.23\cdot 10^{-8}$ | $3.22\cdot 10^{-8}$ | $6.38\cdot 10^{-8}$ | $8.45\cdot 10^{-8}$

Table 2: Computational time for FE solution and RB online phase, with the
speedup and the error.

##  7 Conclusions

In this paper we have developed a VMS-Smagorinsky reduced basis model with
local projection stabilisation (LPS) to stabilise the pressure discretization.
We have proved the stability of the LPS stabilised reduced problem for
piecewise affine pressure discretisations. We use the EIM to approximate the
eddy viscosity and the pressure stabilisation coefficient in order to reduce
the non-linearities that appear in both terms, letting us to efficiently store
the tensors associated to both the VMS-Smagorinsky term and the LPS pressure
term. Also, we have developed an a posteriori error estimator for this problem
via the BRR theory.

Finally, we have presented some numerical results of the RB method constructed
by means of a Greedy algorithm, based upon this error estimator, applied to
the 2D lid-driven cavity problem. We use both Taylor-Hood and
$\mathbb{P}_{2}-\mathbb{P}_{2}$ pairs of velocity-pressure discretizations of
the full order model. Our theory proves that the Taylor-Hood discretisation
reduced problem including LPS stabilisation of the pressure is stable. We have
observed that in both strategies (including or not pressure supremizers in the
reduced velocity spaces), the errors obtained for the velocity and pressure
reduced solutions are quite close. Further, the computational speedup is
higher, with increases of nearly 50%, when we do not consider the enrichment
with pressure supremizers. We thus have built a RB method with enhanced
reduction and speedup, while keeping error levels close to using pressure
supremizers.

## Acknowledgements

This work has been partially supported by the Andalusian Government - Feder
Fund Project US-1254587, that provided the funding to hire Dr. Delgado for
some time, and the Spanish Government - Feder Fund Project
RTI2018-093521-B-C31.

## References

  * [1] S. Ali, F. Ballarin, and G. Rozza, Stabilized reduced basis methods for parametrized steady Stokes and Navier–Stokes equations, Computers & Mathematics with Applications, 80 (2020), pp. 2399–2416.
  * [2] F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza, Supremizer stabilization of POD-Galerkin approximation of parametrized steady incompressible Navier-Stokes equations., Int J Numer Methods Eng, (2014), pp. 1136–1161.
  * [3] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera, An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations, C.R. Acad. Sci. Paris Sér. I Math., 339 (2004), pp. 667–672.
  * [4] C. Bernardi, T. Chacón, and M. Gómez, Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier-Stokes equations, SEMA Journal, (2012), pp. 51–74.
  * [5] C. Bernardi, Y. Maday, and F. Rapetti, Discrétisations variationnelles de problèmes aux limites elliptiques, Springer Berlin Heidelberg, 2004.
  * [6] F. Brezzi and J. Pitkäranta, On the stabilization of finite element approximations of the stokes equations, in Efficient Solutions of Elliptic Systems, Vieweg+Teubner Verlag, 1984, pp. 11–19.
  * [7] F. Brezzi, J. Rappaz, and P. Raviart, Finite dimensional approximation of nonlinear problems, Numer. Maht., 36 (1980), pp. 1–25.
  * [8] A. N. Brooks and T. J. Hughes, Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations, Computer Methods in Applied Mechanics and Engineering, 32 (1982), pp. 199–259.
  * [9] A. Caiazzo, T. Iliescu, V. John, and S. Schyschlowa, A numerical investigation of velocity–pressure reduced order models for incompressible flows, Journal of Computational Physics, 259 (2014), pp. 598–616.
  * [10] T. Chacón Rebollo, A term by term stabilization algorithm for finite element solution of incompressible flow problems, Numerische Mathematik, 79 (1998), pp. 283–319.
  * [11] T. Chacón Rebollo, E. Delgado Ávila, M. Gómez Mármol, F. Ballarin, and G. Rozza, On a certified smagorinsky reduced basis turbulence model, SIAM Journal on Numerical Analysis, 55 (2017), pp. 3047–3067.
  * [12] 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 Journal of Numerical Analysis, 33 (2012), pp. 974–1007.
  * [13] T. Chacón Rebollo, M. Gómez Mármol, and S. Rubino, Numerical analysis of a finite element projection-based VMS turbulence model with wall laws, Computer Methods in Applied Mechanics and Engineering, 285 (2015), pp. 379–405.
  * [14] T. Chacón Rebollo and R. Lewandowski, Mathematical and Numerical Foundations of Turbulence Models and Applications, Springer New York, 2014.
  * [15] R. Codina, Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods, Computer Methods in Applied Mechanics and Engineering, 190 (2000), pp. 1579–1599.
  * [16] J. Edmonds, Matroids and the greedy algorithm, Math. Program., 1 (1971), pp. 127–136.
  * [17] L. P. Franca and C. Farhat, Bubble functions prompt unusual stabilized finite element methods, Computer Methods in Applied Mechanics and Engineering, 123 (1995), pp. 299–308.
  * [18] M. A. Grepl, Y. Maday, N. C. Nguyen, and A. T. Patera, Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations, ESAIM: Mathematical Modelling and Numerical Analysis, 41 (2007), pp. 575–605.
  * [19] F. Hecht, New development in FreeFem++, J. Numer. Math., 20 (2012), pp. 251–265.
  * [20] T. J. Hughes, L. P. Franca, and G. M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. the galerkin/least-squares method for advective-diffusive equations, Computer Methods in Applied Mechanics and Engineering, 73 (1989), pp. 173–189.
  * [21] T. J. R. Hughes, A. A. Oberai, and L. Mazzei, Large eddy simulation of turbulent channel flows by the variational multiscale method, Physics of Fluids, 13 (2001), pp. 1784–1799.
  * [22] J. Novo and S. Rubino, Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows, SIAM Journal on Numerical Analysis, 59 (2021), pp. 334–369.
  * [23] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, vol. 23, Springer Berlin Heidelberg, 1994.
  * [24] G. Rozza and K. Veroy, On the stability of the reduced basis method for Stokes equations in parametrized domains, Comput. Meth. Appl. Mech. Engrg., (2007), pp. 1244–1260.
  * [25] S. Rubino, Numerical modeling of turbulence by Richardson number-based and VMS models, PhD thesis, University of Seville, 2014.
  * [26] K. Veroy and A. T. Patera, Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds, Int. J. Numer. Methods Fluids, 47 (2005), pp. 773–788.


# On entire $f$-maximal graphs in the Lorentzian product
$\mathbb{G}^{n}\times\mathbb{R}_{1}$

H. V. Q. An Duy Tan University, Danang, Vietnam D. V. Cuong Duy Tan
University, Danang, Vietnam N. T. M. Duyen College of Education, Hue
University, Hue, Vietnam D. T. Hieu T. L. Nam Dong Thap University, Dong Thap,
Vietnam

###### Abstract

In the Lorentzian product $\mathbb{G}^{n}\times\mathbb{R}_{1},$ we give a
comparison theorem between the $f$-volume of an entire $f$-maximal graph and
the $f$-volume of the hyperbolic $H_{r}^{+}$ under the condition that the
gradient of the function defining the graph is bounded away from 1. This
condition comes from an example of non-planar entire $f$-maximal graph in
$\mathbb{G}^{n}\times\mathbb{R}_{1}$ and is equivalent to the hyperbolic angle
function of the graph being bounded. As a consequence, we obtain a Calabi-
Bernstein type theorem for $f$-maximal graphs in
$\mathbb{G}^{n}\times\mathbb{R}_{1}.$

AMS Subject Classification (2000): Primary 53C42; Secondary 53C50; 53C25
Keywords: Lorentzian product, Calabi-Bernstein’s Theorem, Gauss space,
$f$-maximal graphs

##  1 Introduction

Suppose that $\Sigma$ is the graph of a $C^{2}$ function $f$ over a domain
$U\subset\mathbb{R}^{n}.$ The graph $\Sigma$ is minimal if the function $f$
satisfies the Minimal Surface Equation

$\operatorname{div}\left(\frac{\nabla f}{\sqrt{1+|\nabla f|^{2}}}\right)=0.$

S. N. Bernstein (1915-1917) proved a surprising theorem in the case $n=2$,
named after him, that such a function defined over $\mathbb{R}^{2}$ must be
affine. This means that a minimal graph over $\mathbb{R}^{2}$ must be a plane.
He conjectured that the theorem holds true for $n>2.$

Bernstein’s conjecture has been a longstanding problem and it was proved to be
true for $n\leq 7$ by the works of De Giorgi, F. Almgren and J. Simons (see
[6], [11] and [21]). However, Bombieri, De Giorgi and Giusti [7] showed the
existence of some entire minimal graphs other than the hyperplanes for $n\geq
8.$

The Lorentzian version of the Bernstein’s Theorem for maximal graphs in the
Lorentz-Minkowski spaces $\mathbb{R}^{n+1}_{1}$ is called the Calabi-
Bernstein’s Theorem. Quite different from the Euclidean case, the Calabi-
Bernstein’s Theorem holds true for any dimension. The theorem was first proved
by Calabi [8] for the case $n=2$ and later by Cheng and Yau [10] for the case
$n>2.$

The Calabi-Bernstein’s Theorem has been generalized to Lorentzian product
spaces (see [2], [3] and [4], for instance) as well as to other ambient spaces
such as warped products, Robertson-Walker spacetimes, manifolds with
density…(see [5], [9], [14], [20] and [22], for instance). Albujer and Alías
[3] proved a Calabi-Bernstein type result, that any entire maximal graph in
$M^{2}\times\mathbb{R}_{1},$ where $M$ is a Riemannian manifold with non-
negative Gaussian curvature, must be totally geodesic. They also proved that
if $M$ is non-flat, then the graph must be a slice
$M\times\\{t_{0}\\},t_{0}\in\mathbb{R}.$ In [1], some non-planar entire
maximal graphs in $\mathbb{H}^{2}\times\mathbb{R}_{1}$ were constructed. This
shows that without the assumption of nonnegative Gauss curvature of $M,$ the
theorem is no longer true.

The Calabi-Bernstein’s Theorem has been also generalized to manifolds with
density. L. Wang (see [22]) proved a Bernstein type theorem for self-shrinkers
in $\mathbb{R}^{n},$ i.e., for $f$-minimal hypersurfaces in Gauss space
$\mathbb{G}^{n}.$ By a calibration argument, without using second order
differential equations, the last two authors of this paper (see [14]) proved a
Bernstein type theorem for entire $f$-minimal graphs in
$\mathbb{G}^{n}\times\mathbb{R}.$ In these cases the theorem holds true
without any condition.

The purpose of this paper is to establish a similar Calabi-Bernstein type
theorem for entire $f$-maximal graphs in the Lorentzian product
$\mathbb{G}^{n}\times\mathbb{R}_{1}$ by using Lorentzian calibration
arguments.

The results are quite different from the case of the Riemannian product. The
first main result of the paper is that an example of a non-planar entire
$f$-maximal graph in $\mathbb{G}^{n}\times\mathbb{R}_{1}$ was found. We see
that along the first coordinate axis, the gradient of the function determining
the graph converges to 1. So it seems the condition for an entire $f$-maximal
graph in $\mathbb{G}^{n}\times\mathbb{R}_{1}$ being a hyperplane is that the
gradient of the function determining the graph is bounded away from 1, or
equivalently, the hyperbolic angle function of the graph is bounded. In fact,
with this condition, a comparison theorem between the $f$-volume of an entire
$f$-maximal graph and the $f$-volume of the hyperbolobic $H_{r}^{+}$ is given.
As a consequence, we obtain the second main result, a Calabi-Bernstein type
theorem for $f$-maximal graphs in $\mathbb{G}^{n}\times\mathbb{R}_{1}.$

For more details about Lorentz-Minkowski spaces, maximal surface,
calibrations, manifolds with density and Gauss space, we refer the reader to
[12], [13], [14], [15], [16], [17], [18] and [19].

##  2 $f$-maximal graphs

Let $\mathbb{R}^{n+1}_{1}$ be the $(n+1)$-dimensional Lorentz-Minkowski space,
i.e., $\mathbb{R}^{n+1}$ endowed with the Lorentzian scalar product

$\langle\ ,\ \rangle=dx_{1}^{2}+dx_{2}^{2}+\ldots+dx_{n}^{2}-dx_{n+1}^{2}.$

Since $\langle,\rangle$ is non-positive definite, for
$\textbf{x}\in\mathbb{R}^{n+1}_{1},$ $\langle\textbf{x},\textbf{x}\rangle$ may
be zero or negative. A nonzero vector $\textbf{x}\in\mathbb{R}^{n+1}_{1}$ is
called spacelike, lightlike or timelike if
$\langle\textbf{x},\textbf{x}\rangle>0$,
$\langle\textbf{x},\textbf{x}\rangle=0$ or
$\langle\textbf{x},\textbf{x}\rangle<0$, respectively. Two vectors
$\textbf{x},\textbf{y}\in\mathbb{R}^{n+1}_{1}$ are said to be Lorentzian
orthogonal if $\langle\textbf{x},\textbf{y}\rangle=0.$ The norm of a tangent
vector x is defined by
$\|\textbf{x}\|=\sqrt{|\langle\textbf{x},\textbf{x}\rangle|}.$

We say that two timelike tangent vectors x and y lie in the same timelike cone
if $\langle\textbf{x},\textbf{y}\rangle<0.$ If x and $\textbf{e}_{n+1},$ where
$\textbf{e}_{n+1}$ denotes the timelike coordinate vector field, lie in the
same timelike cone, we say that x is future-directed. For timelike tangent
vectors x and $\textbf{y},$ we have the Cauchy-Schwarz inequality

$|\langle\textbf{x},\textbf{y}\rangle|\geq\sqrt{-\langle\textbf{x},\textbf{x}\rangle}\sqrt{-\langle\textbf{y},\textbf{y}\rangle},$

and the equality holds if and only if x and y are proportional. If x and y lie
in the same timelike cone, then there exists a unique number $\varphi\geq 0$
such that

$\langle\textbf{x},\textbf{y}\rangle=-|\textbf{x}||\textbf{y}|\cosh\varphi.$

The number $\varphi$ is called the hyperbolic angle between x and
$\textbf{y}.$

A hypersurface in $\mathbb{R}^{n+1}_{1}$ is called spacelike if its induced
metric from $\mathbb{R}^{n+1}_{1}$ is Riemannian. Equivalently, every tangent
vector of the hypersurface is spacelike. A spacelike hypersurface is
orientable, i.e., there exists a unit timelike normal vector field $N$ on it.
Moreover, we can suppose that $N$ is future-directed and then we can define
the hyperbolic angle function $\theta$ of the hypersurface. The value of
$\theta$ at a point $p$ on the hypersurface is the hyperbolic angle between
$N(p)$ and $\textbf{e}_{n+1}.$

Let $\Sigma$ be a spacelike hypersurface in $\mathbb{R}^{n+1}_{1}$ and $N$ be
the unit timelike normal vector field. Denote by $\overline{D}$ and $D$ the
Levi-Civita connections of $\mathbb{R}^{n+1}_{1}$ and $\Sigma,$ respectively.
For $X,Y\in\mathfrak{X}(\Sigma),$ we have

$\overline{D}_{X}Y=D_{X}Y-\langle A(X),Y\rangle N,$

and

$A(X)=-\overline{D}_{X}N,$

where $A:\mathfrak{X}(\Sigma)\rightarrow\mathfrak{X}(\Sigma)$ is the shape
operator associated to $N.$ The mean curvature of $\Sigma,$ denoted by $H,$ is
defined by

$H=-\frac{1}{n}\operatorname{trace}(A).$

It is well-known that

$nH=-\operatorname{div}N.$

A spacelike hypersurface is called _maximal_ if its mean curvature $H$ is zero
everywhere. When $\Sigma$ is the graph of a function
$u:U\rightarrow\mathbb{R},$ where $U\subset\mathbb{R}^{n},$ then it is easy to
see that

$N=\frac{1}{\sqrt{1-|\nabla u|^{2}}}(\nabla u+\textbf{e}_{n+1}),$

is the future-directed unit timelike normal vector field on $\Sigma,$ where
$\nabla u$ denotes the gradient of the function $u$ on $U.$ Note that since
$\Sigma$ is spacelike, $|\nabla u|<1.$ Hence the mean curvature of $\Sigma$ is
computed on $\mathbb{R}^{n}$ to be

$H=-\frac{1}{n}{\rm div}_{{}_{\mathbb{R}^{n}}}\left(\frac{\nabla
u}{\sqrt{1-|\nabla u|^{2}}}\right),$

where ${\rm div}_{{}_{\mathbb{R}^{n}}}$ is the divergence operator on
$\mathbb{R}^{n}.$

Therefore, $\Sigma$ is maximal if and only if $u$ satisfies the maximal
equation

${\rm div}_{{}_{\mathbb{R}^{n}}}\left(\frac{\nabla u}{\sqrt{1-|\nabla
u|^{2}}}\right)=0.$

A density on $\mathbb{R}^{n+1}_{1}$ is a positive function $e^{-f},$ and it is
used to weight volumes of $k$-dimensional spacelike submanifolds. In
$\mathbb{R}^{n+1}_{1}$ with density $e^{-f},$ the weighted volume, or
$f$-volume, of a spacelike hypersurface $\Sigma,$ denoted by
$\operatorname{Vol}_{f}(\Sigma),$ is defined by

$\operatorname{Vol}_{f}(\Sigma)=\int_{\Sigma}e^{-f}dV_{\Sigma},$

where $dV_{\Sigma}$ is the volume element of the hypersurface. The weighted
mean curvature or $f$-mean curvature of $\Sigma,$ denoted by $H_{f},$ is
defined by

$H_{f}=H+\frac{1}{n}\langle\nabla f,N\rangle.$

Now, $\Sigma$ is called _$f$ -maximal_ provided that $H_{f}=0$ everywhere,
i.e., $H=-\frac{1}{n}\langle\nabla f,N\rangle.$

If $\Sigma$ is the graph of a function $u:U\rightarrow\mathbb{R},$ then
$\Sigma$ is $f$-maximal if and only if $u$ satisfies the $f$-maximal equation

${\rm div}_{{}_{\mathbb{R}^{n}}}\left(\frac{\nabla u}{\sqrt{1-|\nabla
u|^{2}}}\right)-\langle\nabla f,N\rangle=0.$

The Gauss space $\mathbb{G}^{n},$ a typical example of manifolds with density,
is just $\mathbb{R}^{n}$ with Gaussian probability density

$e^{-f(\textbf{x})}=e^{c-\frac{|{\bf x}|^{2}}{2}},$

where $c=\log[(2\pi)^{-\frac{n}{2}}].$ Gauss space has many applications to
probability and statistics (see [17] and [18], for instance).

The space $\mathbb{G}^{n}\times\mathbb{R}_{1}$ is
$\mathbb{R}^{n+1}_{1}=\mathbb{R}^{n}\times\mathbb{R}_{1}$ endowed with the
Gaussian-Euclidean density

$e^{-f(\textbf{x},t)}=e^{c-\frac{|{\bf x}|^{2}}{2}},$

where $(\textbf{x},t)\in\mathbb{G}^{n}\times\mathbb{R}_{1},\
\textbf{x}\in\mathbb{G}^{n},\ t\in\mathbb{R}_{1}.$

Below is the example about a non-planar entire $f$-maximal graph in
$\mathbb{G}^{n}\times\mathbb{R}_{1}$ mentioned in the introduction section.

Consider the entire graph $\Sigma_{0}$ of the function
$u:\mathbb{G}^{n}\longrightarrow\mathbb{R},\
\textbf{x}\longmapsto\displaystyle\int_{0}^{x_{1}}\sqrt{\dfrac{e^{\tau^{2}}}{1+e^{\tau^{2}}}}d\tau,$
where $\textbf{x}=(x_{1},x_{2},\ldots,x_{n}).$ The graph $\Sigma_{0}$ has a
parametrization as follows

$\displaystyle X(\textbf{x},t)=\left(\textbf{x},u(\textbf{x})\right).$

Since

$1-|\nabla u|^{2}=\frac{1}{1+e^{x_{1}^{2}}}>0,$

$\Sigma_{0}$ is spacelike.

It follows that the future-directed unit timelike normal vector field is

${N}=\left(\dfrac{u_{x_{1}}}{\sqrt{1-u_{x_{1}}^{2}}},0,\ldots,0,\dfrac{1}{\sqrt{1-u_{x_{1}}^{2}}}\right).$

Therefore, the mean curvature of $\Sigma_{0}$ is

$H=-\frac{1}{n}{\rm
div}N=-\frac{1}{n}\frac{u_{x_{1}x_{1}}}{(1-u_{x_{1}}^{2})^{3/2}}=-\frac{x_{1}e^{x_{1}^{2}/2}}{n}.$

Since

$\langle\nabla
f,N\rangle=\frac{x_{1}u_{x_{1}}}{\sqrt{1-u_{x_{1}}^{2}}}=x_{1}e^{x_{1}^{2}/2},$
$H_{f}=H+\frac{1}{n}\langle\nabla f,N\rangle=0,$

i.e., $\Sigma_{0}$ is $f$-maximal.

This example shows that we need an assumption in order to obtain a uniqueness
result.

##  3 Calabi-Bernstein result

Let $\Sigma\subset\mathbb{G}^{n}\times\mathbb{R}_{1}$ be the graph of a
function $u({\textbf{x}})=t$ over $\mathbb{G}^{n},$ that is $f$-maximal, and
$N$ be the future-directed unit timelike normal vector field of $\Sigma.$
Consider the smooth extension of $N$ by translations along $t$-axis, that we
will also denote by $N,$ and the $n$-differential form $w$ defined by

$w=i_{N}dV.$

Similarly as in the Riemannian case, we obtain the following,

######  Lemma 1.

The form $w$ is a Lorentzian calibration that calibrates $\Sigma,$ i.e.,

  1. 1.

$d(e^{-f}w)=0;$

  2. 2.

For orthonormal spacelike vector fields $X_{i},\ i=1,2,\ldots,n,$
$|w(X_{1},X_{2},\ldots,X_{n})|\geq 1,$ and the equality holds at a point if
and only if $X_{1},X_{2},\ldots,X_{n}$ are tangent to $\Sigma$ at such a
point.

###### Proof.

  1. 1.

Because $\Sigma$ is $f$-maximal,

$\displaystyle d(e^{-f}w)$ $\displaystyle=d(e^{-f}i_{N}dV)={\rm
div}(e^{-f}{N})dV$ $\displaystyle=(e^{-f}{\rm div}N-e^{-f}\langle\nabla
f,N\rangle)dV$ $\displaystyle=-e^{-f}(nH+\langle\nabla f,N\rangle)dV=0.$

  2. 2.

At any point $p$ of the Lorentz Minkowski spacetime, we can find a tangent
vector $\overline{N}(p)$ such that
$\\{X_{1}(p),\ldots,X_{n}(p),\overline{N}(p)\\}$ is a frame [19, page 84].
Then we can write

$N=\pm\cosh\varphi\overline{N}+\sinh\varphi u_{0},$

for $\varphi\in\mathbb{R},$ where $\langle N,u_{0}\rangle=0,\ \langle
u_{0},u_{0}\rangle=1,$ and the sign depends on the timelike orientation.
Hence,

$\displaystyle|w(X_{1},\ldots,X_{n})(p)|=|dV(X_{1},\ldots,X_{n},N)(p)|=|dV(X_{1},\ldots,X_{n},\pm\cosh\varphi\overline{N})|=\cosh\varphi.$

The previous equation clearly proves this part of the lemma.

∎

By Stokes’ Theorem, as in the case of $f$-minimal submanifolds [13, Theorem
2.1], it can be proved that $\Sigma$ is $f$-area-maximizing, i.e., any compact
portion of $\Sigma$ has largest area among all spacelike hypersurfaces in its
homology class, the Lorentzian version follows analogously.

As a next step, we establish a comparison between the $f$-volume of $\Sigma$
and the $f$-volume of the spacelike Lorentzian hypersphere
$H_{r}^{+}=\\{({\textbf{x}},t)\in\mathbb{G}^{n}\times\mathbb{R}_{1}:\
\langle\textbf{x},\textbf{x}\rangle-\langle t,t\rangle=-r^{2},\ t>0\\}.$

######  Lemma 2.

If $|\nabla u|$ is bounded away form 1, then for every $r>0$,

$\operatorname{Vol}_{f}(\Sigma)\geq\operatorname{Vol}_{f}(H^{+}_{r}).$

######  Remark 3.

Because $\langle
N,\textbf{e}_{n+1}\rangle=\displaystyle\frac{1}{\sqrt{1-|\nabla u|^{2}}},$ it
follows that $|\nabla u|$ is bounded from 1 is equivalent to the hyperbolic
angle function $\theta$ of $\Sigma$ being bounded. This condition appeared in
some Calabi-Bernstein type results for constant mean curvature spacelike
hypersurfaces in Lorentzian warped products (see [5], for instance).

###### Proof.

Figure 1: For the proof of Lemma 2

Since the density does not depend on the last coordinate, the $f$-volume and
$f$-mean curvature are invariant under translations along the $x_{n+1}$-axis.
Therefore, we can assume that the origin $O\in\Sigma.$ Let ${B^{n}_{R}}$ be
the $n$-ball in $\mathbb{G}^{n}$ with center $O$ and radius $R$ and ${\cal
C}=B^{n}_{R}\times\mathbb{R}$ be the $(n+1)$-cylinder over ${B^{n}_{R}}.$
Since $\Sigma$ is spacelike, $\widetilde{\Sigma}:=\Sigma\cap{\cal C}$ is
bounded and lies in the slab between parallel hyperplanes $t=\pm R.$ Let
$\Omega$ be the region bounded by ${\cal C},\widetilde{\Sigma}$ and
$\widetilde{H^{+}_{r}}:=H^{+}_{r}\cap{\cal C}$ with inward-pointing
orientation; $\widetilde{\partial\cal C}$ be the part of $\partial{\cal C}$
between $\widetilde{\Sigma}$ and $\widetilde{H^{+}_{r}}$ and $w$ be the
Lorentzian calibration constructed as above. By Stokes’ Theorem (note that,
the directions of both $\Sigma$ and $H_{r}^{+}$ are future-directed), we
obtain

$0=\int_{\Omega}d(e^{-f}w)=\int_{\partial\Omega}e^{-f}w=\int_{\widetilde{\Sigma}}e^{-f}w-\int_{\widetilde{H_{r}^{+}}}e^{-f}w+\int_{\widetilde{\partial\cal
C}}e^{-f}w,$

or

$\int_{\widetilde{H_{r}^{+}}}e^{-f}w=\int_{\widetilde{\Sigma}}e^{-f}w+\int_{\widetilde{\partial\cal
C}}e^{-f}w.$

By Lemma 1,
$\int_{\widetilde{H_{r}^{+}}}e^{-f}w\geq\operatorname{Vol}_{f}({\widetilde{H_{r}^{+}}})$
and
$\int_{\widetilde{\Sigma}}e^{-f}w=\operatorname{Vol}_{f}(\widetilde{\Sigma}).$
Thus, we obtain the following estimate

$\operatorname{Vol}_{f}(\widetilde{H_{r}^{+}})\leq\operatorname{Vol}_{f}(\widetilde{\Sigma})+\int_{\widetilde{\partial\cal
C}}e^{-f}w.$

Since
$\displaystyle\lim_{R\rightarrow\infty}\operatorname{Vol}_{f}(\widetilde{H_{r}^{+}})=\operatorname{Vol}_{f}(H_{r}^{+})$
and
$\displaystyle\lim_{R\rightarrow\infty}\operatorname{Vol}_{f}(\widetilde{\Sigma_{R}})=\operatorname{Vol}_{f}(\Sigma),$
so the remain we have to prove is
$\displaystyle\lim_{R\rightarrow\infty}\int_{\widetilde{\partial\cal
C}}e^{-f}w=0.$

Let $\overline{\partial{\cal
C}}=S_{R}^{n-1}\times\left[-~{}\sqrt{R^{2}+r^{2}},\sqrt{R^{2}+r^{2}}\;\right]$
be the part of $\partial{\cal C}$ between two parallel hyperplanes
$t=\pm\sqrt{R^{2}+r^{2}}.$ Because $\widetilde{\partial{\cal
C}}\subset\overline{\partial\cal C}$ and the density on
$\overline{\partial{\cal C}}$ is $e^{c-\frac{R^{2}}{2}}$ (a constant), we get

$\int_{\widetilde{\partial\cal C}}e^{-f}w\leq
e^{c-\frac{R^{2}}{2}}\int_{\overline{\partial\cal C}}w.$

Note that a unit normal vector field of ${\cal C}$ is of the form
$(\textbf{x},0),$ where $|\textbf{x}|=1.$ Therefore

$\int_{\overline{\partial\cal C}}w=\int_{\overline{\partial{\cal
C}}}\frac{\langle\nabla u,\textbf{x}\rangle}{\sqrt{1-|\nabla u|^{2}}}.$

By the assumption $|\nabla u|$ is bounded away form 1, we have

$\int_{\overline{\partial{\cal C}}}\frac{\langle\nabla
u,\textbf{x}\rangle}{\sqrt{1-|\nabla u|^{2}}}\leq
K\operatorname{Vol}(\overline{\partial{\cal
C}})=2K\sqrt{R^{2}+r^{2}}\operatorname{Vol}(S^{n-1}_{R}),$

where $K>0$ is a constant. Finally, we get the following inequality

$\int_{\widetilde{\partial\cal C}}e^{-f}w\leq
2e^{c-\frac{R^{2}}{2}}K\sqrt{R^{2}+r^{2}}\operatorname{Vol}(S^{n-1}_{R}).$

It is easy to see that, the right-hand side of the inequality goes to zero
when $R$ approaches infinity and therefore the proof of the lemma is complete.
∎

As a consequence, we obtain the following Calabi-Bernstein type theorem.

######  Theorem 4.

In $\mathbb{G}^{n}\times\mathbb{R}_{1},$ the graph $\Sigma$ of a function
$u(\text{\bf x})=t$ over $\mathbb{G}^{n},$ with $|\nabla u|$ is bounded away
from 1, is $f$-maximal if and only if $u$ is constant, i.e., $\Sigma$ is a
spacelike hyperplane.

###### Proof.

It is clear that, if $u$ is constant, then $\Sigma$ is $f$-maximal.

Note that $H_{r}^{+}$ is the graph of the function
$g(\textbf{x})=\sqrt{{\textbf{x}}^{2}+r^{2}},$ so

$\operatorname{Vol}_{f}(H_{r}^{+})=\int_{\mathbb{R}^{n}}e^{-f}\sqrt{\frac{r^{2}}{\textbf{x}^{2}+r^{2}}}dV_{{}_{\mathbb{R}^{n}}}.$

It is easy to see that

$\lim_{r\rightarrow\infty}\operatorname{Vol}_{f}(H_{r}^{+})=1.$

Since

$\displaystyle 1$
$\displaystyle=\operatorname{Vol}_{f}\mathbb{G}^{n}=\int_{\mathbb{R}^{n}}e^{-f}dV_{{}_{\mathbb{R}^{n}}}$
$\displaystyle\geq\int_{\mathbb{R}^{n}}e^{-f}\sqrt{1-|\nabla
u|^{2}}dV_{{}_{\mathbb{R}^{n}}}=\operatorname{Vol}_{f}(\Sigma)$
$\displaystyle\geq\lim_{r\rightarrow\infty}\operatorname{Vol}_{f}(H_{r}^{+})=1,$

we conclude that $|\nabla u|^{2}=0,$ i.e., $u$ is constant. ∎

Acknowledgements. The authors would like to thank the referee for carefully
reading our paper and for his/her valuable comments and suggestions which
helped to improve the manuscript. This research is funded by Vietnam National
Foundation for Science and Technology Development (NAFOSTED) under grant
number 101.04.2014.26.

## References

  * [1] A. L. Albujer, New examples of entire maximal graphs in $\mathbb{H}^{2}\times\mathbb{R}_{1}$, Differential Geom. Appl. 26 (2008), no. 4, 456-462.
  * [2] A. L. Albujer, L. J. Alías, Calabi-Bernstein results and parabolicity of maximal surfaces in Lorentzian product spaces, Recent trends in Lorentzian geometry, 49-85, Springer Proc. Math. Stat., 26, Springer, New York, 2013.
  * [3] A. L. Albujer, L. J. Alías, Calabi-Bernstein results for maximal surfaces in Lorentzian product spaces, J. Geom. Phys. 59 (2009), no. 5, 620-631.
  * [4] J. A. Aledo, A. Romero, R. M. Rubio, The classical Calabi-Bernstein Theorem revisited, J. Math. Anal. Appl. 431 (2015), no. 2, 1172-1177.
  * [5] J. A. Aledo, A. Romero, R. M. Rubio, Constant mean curvature spacelike hypersurfaces in Lorentzian warped products and Calabi-Bernstein type problems, Nonlinear Anal. 106 (2014), 57-69.
  * [6] F. J. Almgren, _Some interior regularity theorems for minimal surfaces and an extension of Bernstein’s theorem_ , Ann. of Math. 84 (1966), 277-292.
  * [7] E. Bombieri, E. De Giorgi, E. Giusti, _Minimal cones and the Bernstein theorem_ , Inventiones Math. 7 (1969), 243-269.
  * [8] E. Calabi, Examples of Bernstein problems for some nonlinear equations, in Global Analysis (Proc. Sympos. Pure Math., Vol. XV, Berkeley, Calif., 1968), Amer. Math. Soc., Providence, R.I., 1970 pp. 223-230.
  * [9] M. P. Cavalcante, H. F. de Lima, M. S. Santos, New Calabi-Bernstein type results in weighted generalized Robertson-Walker spacetimes, Acta Math. Hungar. 145 (2015), no. 2, 440-454.
  * [10] S. Y. Cheng, S. T. Yau, Differential equations on Riemannian manifolds and their geometric applications, Comm. Pure Appl. Math. 28 (1975), no. 3, 333-354.
  * [11] E. de Giorgi, _Una extensione del teorema di Bernstein_ , Ann. Scuola Norm. Sup. Pisa 19 (1965), 79-85.
  * [12] R. Harvey and H. B. Lawson, Jr., Calibrated geometries, Acta Math. 148 (1982) 47-157.
  * [13] D. T. Hieu, Some calibrated surfaces in manifolds with density, J. Geom. Phys. 61 (2011), no. 8, 1625-1629.
  * [14] D. T. Hieu, T. L. Nam, Bernstein type theorem for entire weighted minimal graphs in $\mathbb{G}^{n}\times\mathbb{R}$, J. Geom. Phys. 81 (2014), 87-91.
  * [15] R. López, Differential geometry of curves and surfaces in Lorentz-Minkowski space, Int. Electron. J. Geom. 7 (2014), no. 1, 44-107.
  * [16] J. Mealy, Calibrations on semi-Riemannian manifolds, Thesis (Ph.D.)-Rice University. ProQuest LLC, Ann Arbor, MI, 1989. 72 pp.
  * [17] F. Morgan, Manifolds with density, Notices Amer. Math. Soc. 52 (2005), 853-858.
  * [18] F. Morgan, Geometric Measure Theory: a Beginner’s Guide, 4th ed. Academic Press, London, 2009.
  * [19] B. O’Neil, Semi-Riemannian Geometry with Applications to Relativity, Academic Press, London (1983).
  * [20] A. Romero, Maximal graphs and Calabi-Bernstein’s type problems in some Robertson-Walker spacetimes, Bull. Transilv. Univ. Brasov Ser. III 1(50) (2008), 309-316.
  * [21] J. Simons, _Minimal varieties in Riemannian manifolds_ , Ann. of Math. 88 (1968), 62-105.
  * [22] L. Wang, _A Bernstein type theorem for self-similar shrinkers_ , Geom. Dedicata 151 (2011), 297-303.