An efficient Chorin-Temam projection proper orthogonal decomposition based reduced-order model for nonstationary Stokes equations

Xi Li¹¹1School of Mathematics, Sichuan University, Chengdu, Sichuan 610064, China ([email protected]). The work of this author was supported by the National Natural Science Foundation of China(Grant No. 11971337)., Yan Luo²²2School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, Sichuan 611731, China ([email protected]). The work of this author was supported by the Young Scientists Fund of the National Natural Science Foundation of China(Grant No. 11901078) and by the Fundamental Research Funds for the Central Universities(Grant No. ZYGX2020J021). and Minfu Feng³³3Corresponding author. School of Mathematics, Sichuan University, Chengdu, Sichuan 610064, China ([email protected]). The work of this author was supported by the National Natural Science Foundation of China(Grant No. 11971337).

Abstract

In this paper, we propose an efficient proper orthogonal decomposition based reduced-order model(POD-ROM) for nonstationary Stokes equations, which combines the classical projection method with POD technique. This new scheme mainly owns two advantages: the first one is low computational costs since the classical projection method decouples the reduced-order velocity variable and reduced-order pressure variable, and POD technique further improves the computational efficiency; the second advantage consists of circumventing the verification of classical LBB/inf-sup condition for mixed POD spaces with the help of pressure stabilized Petrov-Galerkin(PSPG)-type projection method, where the pressure stabilization term is inherent which allows the use of non inf-sup stable elements without adding extra stabilization terms. We first obtain the convergence of PSPG-type finite element projection scheme, and then analyze the proposed projection POD-ROM’s stability and convergence. Numerical experiments validate out theoretical results.
Keywords: Projection method; PSPG stabilization; Proper orthogonal decomposition.

1 Introduction

Reduced Order Models(ROMs) have become widespread in the scientific community and owns constant attention in academic community because of its strength and efficiency in alleviating the huge computational costs needed in many modern engineering applications and academic purposes. As one of the most popular ROM approaches, the proper orthogonal decomposition(POD) strategy aims to use singular value decomposition(SVD) to extract the dominant modes(in the sense of energy) from a high-resolution numerical scheme(commonly referred to full order models(FOMs)) to form POD bases/modes. Onto these reduced-order bases, a Galerkin projection can be used to project into the space spanned by these bases/modes to obtain the reduced-order numerical solution, this pioneering work is done by Kunisch and Volkwein in [19, 20] which has laid the foundation of numerical analysis for Galerkin-POD. Other references about Galerkin-POD can be found in [30, 22, 21, 32] and so on. Applications of Galerkin-POD on incompressible flows, especially with the spatial variables discretized by finite element method(FE-POD), has a long history, from Luo et al.[22] obtaining FE-POD-ROM for nonstationary Navier-Stokes equations based on finite difference method and mixed finite element method to discretize temporal and spatial variables, Ravindran[26] analyzing FE-POD-ROM for nonstationary Boussinesq equations, to other finite element methods, like discontinuous Galerkin and hybrid discontinuous Galerkin, for spatial discrete are applied to get FE-POD-ROM[12].
Traditional mixed finite element POD-ROM had encountered a problem, that is, when considering the classical mixed finite element scheme as the snapshot source of the POD, the pressure variable is eliminated because of the classical LBB/inf-sup condition, so only the velocity finite element solution is used as the snapshot of the POD technique to construct the POD bases/modes, which makes the final obtained FE-POD-ROM does not contain the pressure variable. This makes it improper to use POD technique directly to reduce the model order in some engineering and academic models that need the pressure variables. Among that problem, the pressure field must be reconstructed a posterior using some pressure recovery techniques, like pressure Poisson equation[23], Supremizer stabilized technique[2] and so on. However, just as stated in [11] and [5], there are some deficiencies in these a posterior pressure recovery techniques, such as the usage of those techniques still need the fulfillment of LBB/inf-sup condition which restricts some flexible pairs of mixed finite element spaces, and the un-physical and unclear boundary condition of pressure Poisson equation poses a challenge for the utilization of this technique. With those in mind, recently, a new research direction related to FE-POD has aroused some researchers’ interest, we refer as stabilized FE-POD, which draws inspiration from the successful applications of stabilized techniques in finite element method. Caiazzo et al.[5] firstly proposed to combine stabilization technique, like residual-based streamline-upwind Petrov-Galerkin(SUPG) stabilization, with FE-POD to preserve the finite element pressure solution, so that the POD reduced-order pressure can be obtained from the snapshot of finite element pressure solution. From then on, other stabilized techniques were introduced into FE-POD, like local projection stabilization(LPS) FE-POD[27, 24], artificial compressible(AC) FE-POD[11], streamline diffusion(SD) FE-POD[1] and so on.
In this paper, our aim is to combine the classical projection method and POD technique to propose an efficient projection POD-ROM. As one of the fast decoupled numerical algorithms for solving incompressible fluid, projection method has become a classical one after it was proposed by Chorin and Temam [8, 31] and developed by Shen and Guermond [29, 15, 28, 14] etc. Apart from that, as mentioned by Rannacher in [25], another remarkable feature of this method is that the original Chorin-Temam projection method owns an inherent mechanism of pressure stability. This mechanism was proposed by Rannacher [25] and was realized by Frutos et al. [10]. These two advantages of the classical projection method are merged into the POD technique to get our projection POD-ROM. The main difference in our approach with respect to those stabilized FE-POD-ROM is that the pressure stabilized Petrov-Galerkin(PSPG)-like FE-POD-ROM has the inherent pressure stabilization, which means no additional pressure stabilization term is needed to circumvent the LBB/inf-sup condition of mixed POD spaces which is still an open problem to verify.
The main contribution of this paper is that a more efficient ROM numerical scheme for solving incompressible fluid is obtained by combining the advantageous computational efficiency of both projection method and POD technique. This POD-ROM scheme stems from replacing the finite element functions in the fully discrete scheme of the projection method with the POD reduced-order functions. To this end, we analyze the fully discrete scheme of classical projection method in the form of with inherent pressure stabilized term, and numerically verify some theoretical results concerning with this finite element scheme. Furthermore, we take the finite element solution as snapshots in order to use method of snapshots to construct POD modes and POD spaces, and obtain our projection POD scheme. The stability and convergence of this scheme are analyzed, and which are also validated by some numerical experiments.
The outline of this paper is as follows: In section 2, we introduce some essential notations. In section 3, the classical projection method will be presented, followed by its stability and convergence analysis. In section 4, we propose the projection POD numerical scheme and prove its stability and convergence. In section 5, some numerical experiments will be investigated to confirm the theoritical results analyzed before. Finally, in section 6, we draw conclusions to complete this paper.

2 Preliminaries and notations

we consider the nonstationary Stokes equations

\left\{\begin{aligned} \partial_{t}\bm{u}-\nu\Delta\bm{u}+\nabla p&=\bm{f},\quad\text{in}\;\Omega,\\ \nabla\cdot\bm{u}&=0,\quad\text{in}\;\Omega,\\ \bm{u}&=\bm{0},\quad\text{on}\;\partial\Omega,\\ \bm{u}(0,x)&=\bm{u}_{0}(x),\quad\text{in}\;\partial\Omega,\end{aligned}\right.

(2.1)

where $\Omega$ is an open bounded domain in $\mathbb{R}^{2}$ with a sufficiently smooth boundary $\partial\Omega$ . The unknowns are the vector function $\bm{u}$ (velocity) and the scalar function $p$ (pressure).
In order to state the above problem properly, we need some settings. Let $L^{p}(\Omega),\;1\leq p\leq\infty$ denote the space of standard p-th power absolutely integrable functions with respect to the Lebesgue measure. In particular, $L^{2}(\Omega)$ is a Hilbert space endowed with the scalar product $(\cdot,\cdot)$ and its induced norm $\|\cdot\|_{0}$ . Moreover, we use $\|\cdot\|^{2}_{\ell^{2}(L^{2})}:=\Delta t\sum^{N}_{n=1}\|\cdot\|^{2}_{0}$ to denote the discrete integral in time with respect to $\|\cdot\|^{2}_{0}$ . The Sobolev space $W^{m,p}(\Omega)$ and its norm $\|\cdot\|_{m,p}$ are also standard in the sense of [9]. Furthermore, we use the abbreviation $H^{m}(\Omega):=W^{m,2}(\Omega)$ and $\|\cdot\|_{m}:=\|\cdot\|_{m,2}$ , it’s a Hilbert space with a scalar product. Moreover, we denote by $H^{1}_{0}(\Omega)$ the space of functions of $H^{1}(\Omega)$ with vanishing trace on $\partial\Omega$ and by $H^{-1}(\Omega)$ its dual space. The symbol $\langle\cdot,\cdot\rangle$ in general denotes the duality pairing between the space and it’s dual one. Vector analogues of the Sobolev spaces along with vector valued functions are denoted by boldface letters, for instance $\bm{H}^{m}(\Omega):=\left(H^{m}(\Omega)\right)^{2}$ . We all know that by Poincaré inequality[3, p.135], the seminorm $|u|_{1}=\|\nabla u\|_{0}$ in $H^{1}(\Omega)$ is a norm in $H^{1}_{0}(\Omega)$ .
Throughout this paper, we use $C$ to denote a positive constant independent of $\Delta t$ , $h$ , not necessarily the same at each occurrence.
Let us introduce some convenient spaces. The first one is

\bm{H}({\rm div};\Omega):=\{\bm{u}\in\bm{L}^{2}(\Omega)|\;\nabla\cdot\bm{u}\in\bm{L}^{2}({\Omega})\},

which is a Hilbert space with norm $\|\bm{u}\|_{{\rm div}}:=\|\bm{u}\|_{0}+\|\nabla\cdot\bm{u}\|_{0}$ , and

\bm{H}_{0}({\rm div};\Omega):=\{\bm{u}\in\bm{H}({\rm div};\Omega)|\;\bm{u}\cdot\bm{n}|_{\partial\Omega}=0\}.

For future use, we also need the following spaces, which play a key role in Helmholtz decomposition,

\bm{J}_{0}:=\left\{\bm{v}\in\bm{H}^{1}_{0}(\Omega):\;\nabla\cdot\bm{v}=0\right\},

and

\bm{J}_{1}:=\left\{\bm{v}\in\bm{L}^{2}(\Omega):\;\nabla\cdot\bm{v}=0,\;\text{and}\;\bm{v}\cdot\bm{n}|_{\partial\Omega}=0\right\}.

By means of spaces above, we give the classical Helmholtz decomposition[13, p.29] as

\bm{L}^{2}(\Omega)=\bm{J}_{1}(\Omega)\oplus[\bm{J}_{1}(\Omega)]^{\perp}=\bm{J}_{1}(\Omega)\oplus\{\nabla q:q\in H^{1}(\Omega)\}.

(2.2)

In order to give a variational formulation of problem (2.1), we consider the velocity space $\bm{V}:=\bm{H}^{1}_{0}$ and the pressure space $Q:=L^{2}_{0}=\{q\in L^{2}(\Omega)|\int_{\Omega}q\,{\rm d}x=0\}$ . For a given time constant $T$ , a weak formulation of the nonstationary Stokes equations are expressed as follows: for almost $\forall t\in(0,T]$ , find $(\bm{u},p):(0,T]\rightarrow\bm{V}\times Q$ , such that

\left\{\begin{aligned} (\partial_{t}\bm{u},\bm{v})+\nu(\nabla\bm{u},\nabla\bm{v})+(\nabla p,\bm{v})&=\langle\bm{f},\bm{v}\rangle,\quad\forall\bm{v}\in\bm{V},\\ (\nabla\cdot\bm{u},q)&=0,\quad\forall q\in Q.\\ \bm{u}(x,0)&=\bm{u}^{0}(x).\end{aligned}\right.

(2.3)

Let $\{\mathcal{T}_{h}\}$ be a uniformly regular family of triangulation of $\overline{\Omega}$ (see [9, p.111]) and $h:=\max_{K\in\mathcal{T}_{h}}\{h_{K}|\;h_{K}\\ :={\rm diam}(K)\}$ . Since we want to construct the POD space where snapshots are deriving from both the finite element approximated velocity and also pressure, contrary to the general choice which is to choose the inf-sup stable mixed finite element spaces, we choose equal-order mixed finite element spaces and then it’s obviously not stable in the sense of the classical inf-sup condition, as

Y^{l}_{h}:=\{v_{h}\in C^{0}(\overline{\Omega})|\;v_{h}|_{K}\in\mathbb{P}_{l}(K),\;\forall K\in\mathcal{T}_{h}\},\;l\geq 1,

where $\mathbb{P}_{l}(K)$ is the space of polynomials up to degree $l$ on $K$ . Since there are two types of velocities in the temporal semi-discrete projection scheme, i.e., intermediate velocity $\widetilde{\bm{u}}^{n+1}$ and end-of-step velocity $\bm{u}^{n+1}$ , we need the three finite element spaces to locate three full discrete variables of three unknowns. To this end, we define the equal-order finite element spaces for end-of-step velocity, intermediate velocity and pressure

\bm{V}_{h}:=\bm{Y}^{l}_{h}\cap\bm{V},\quad\bm{W}_{h}:=\bm{Y}^{l}_{h}\cap\bm{J}_{1},\quad Q_{h}:=Y^{l}_{h}\cap Q.

Since the triangulations are assumed to be shape regular, the following inverse inequality holds for each $\bm{v}_{h}\in\bm{V}_{h}$ on each mesh cell $K\in\{\mathcal{T}_{h}\}$ (see [3, Theorem 4.5.11]),

\|\bm{v}_{h}\|_{m,p,K}\leq C_{{\rm inv}}h^{l-m+2(\frac{1}{p}-\frac{1}{q})}_{K}\|\bm{v}_{h}\|_{l,q,K},

(2.4)

where $0\leq l\leq m\leq 1$ , $0\leq q\leq p\leq\infty$ and $h_{k}$ is the size of the mesh cell $K\in\{\mathcal{T}_{h}\}$ . Then, $I_{h}$ (resp. $J_{h}$ ) denotes a bounded linear interpolation operator $I_{h}:\;\bm{W}^{t,q}(K)\rightarrow\bm{V}_{h}$ (resp. $J_{h}:\;W^{t,q}(K)\rightarrow Q_{h}$ ), where $1\leq q\leq\infty$ , $0\leq s\leq t\leq l+1$ , that satisfies the following optimal bounds(see [3, Theorem 4.4.4])

	$\displaystyle\\|\bm{v}-I_{h}\bm{v}\\|_{s,q,K}$	$\displaystyle\leq Ch_{K}^{t-s}\\|\bm{v}\\|_{t,q,K},$		(2.5)
	$\displaystyle\\|q-J_{h}q\\|_{s,q,K}$	$\displaystyle\leq Ch_{K}^{t-s}\\|q\\|_{t,q,K},$		(2.5)

and also interpolation operator $J_{h}$ has following stability

\|\nabla J_{h}q\|_{0}\leq C\|\nabla q\|_{0}.

(2.6)

The Sobolev inclusion $L^{2}(\Omega)\hookrightarrow H^{-1}(\Omega)$ relation implies the following inequality holds:

\|v\|_{-1}\leq C\|v\|_{0},\quad\forall v\in L^{2}(\Omega).

(2.7)

3 Classical projection method

3.1 classical projection scheme

The original projection method consists of firstly finding an intermediate velocity $\widetilde{\bm{u}}^{n+1}$ from the momentum equation without the pressure term, and then utilizing the classical Helmholtz decomposition (2.2) to get an end-of-step velocity $\bm{u}^{n+1}$ which is solenoidal and its normal component $\bm{u}^{n+1}\cdot\bm{n}$ vanishes. Meanwhile, as a “by-product” of the above decomposition, a discrete pressure $p^{n+1}$ can be also get for which some investigations have been made to confirm that this “by-product” pressure is indeed the proper approximation of the exact pressure $p(t_{n+1})$ , see [25].
For a given timestep $\Delta t$ , we consider a uniform discretization of the time interval $[0,T]$ as $0=t_{0}<t_{1}<\cdots<t_{N}$ where $N=\lfloor T/\Delta t\rfloor$ is the rounding $T/\Delta t$ down. In mathematics, the above procedure can be interpreted as the following classical Chorin-Teman projection method[7, 31]:
(Scheme 1 - Semi-discrete of classical projection method) For $n=0,1,\ldots,N-1$ ,
Step 1- Given $\bm{u}^{n}\in\bm{J}_{0}$ , find $\widetilde{\bm{u}}^{n+1}\in\bm{V}$ that satisfies

\frac{\widetilde{\bm{u}}^{n+1}-\bm{u}^{n}}{\Delta t}-\nu\Delta\widetilde{\bm{u}}^{n+1}=\bm{f}^{n+1},\quad\text{in}\;\Omega.

Step 2- Given $\widetilde{\bm{u}}^{n+1}\in\bm{V}$ , compute $(\bm{u}^{n+1},p^{n+1})\in\bm{H}_{0}({\rm div},\Omega)\times\left[Q\cap H^{1}(\Omega)\right]$ from

\left\{\begin{aligned} \bm{u}^{n+1}+\Delta t\nabla p^{n+1}&=\widetilde{\bm{u}}^{n+1},\quad\text{in}\;\Omega,\\ \nabla\cdot\bm{u}^{n+1}&=0,\quad\text{in}\;\Omega,\\ \bm{u}^{n+1}\cdot\bm{n}&=0,\quad\text{on}\;\partial\Omega.\end{aligned}\right.

The usage of classical Helmholtz decomposition is reflected in decomposing the intermediate velocity $\widetilde{\bm{u}}^{n+1}$ from step 1 to get a solenoidal velocity $\bm{u}^{n+1}$ and a gradient of a $H^{1}$ -function, which can be stated as:
$\mathbf{Step\ 2^{\prime}}$ - Projecting $\widetilde{\bm{u}}^{n+1}$ onto the space $\bm{J}_{0}$

\bm{u}^{n+1}=P_{\bm{J}_{0}}(\widetilde{\bm{u}}^{n+1}).

Remark 3.1.

Let us explain the relation between $\bm{Step2}$ and $\bm{Step2^{\prime}}$ . Followed by (2.2), decomposing $\widetilde{\bm{u}}^{n+1}$ would directly gets:

\widetilde{\bm{u}}^{n+1}=\bm{u}^{n+1}+\nabla\phi^{n+1},

(3.8)

where $\bm{u}^{n+1}\in\bm{J}_{1}$ implies $\nabla\cdot\bm{u}^{n+1}=0$ and $\bm{u}^{n+1}\cdot\bm{n}|_{\partial\Omega}=0$ , which is exactly the second and third equations in $\bm{Step\ 2}$ . Multiplying the previous decomposition by $1/\Delta t$ and adding it to the first equation in $\bm{Step\ 1}$ , we get $(\bm{u}^{n+1}-\bm{u}^{n})/\Delta t-\nu\Delta\widetilde{\bm{u}}^{n+1}+\nabla(\phi^{n+1}/\Delta t)=\bm{f}^{n+1}$ , which means $\phi^{n+1}/\Delta t$ can be regarded as a proper approximation of $p(t_{n+1})$ , i.e., we can define $p^{n+1}:=\phi^{n+1}/\Delta t$ . Thus, on the one hand, we can compute $\phi^{n+1}$ in the former decomposition (3.8) and using $p^{n+1}:=\phi^{n+1}/\Delta t$ to get $p^{n+1}$ ; on the other hand, we can also substitute $\phi^{n+1}$ with $\Delta tp^{n+1}$ in (3.8) to directly compute $p^{n+1}$ , which is exactly what the first equation in step 2 does.

Or, Step 2 is equivalent to the following pressure Poisson equation using $\nabla\cdot\bm{u}^{n+1}=0$ :
$\mathbf{Step\ 2^{\prime\prime}}$ - Given $\widetilde{\bm{u}}^{n+1}\in\bm{V}$ , compute $p^{n+1}\in Q\cap H^{2}(\Omega)$ from

\left\{\begin{aligned} \Delta p^{n+1}&=\frac{1}{\Delta t}\nabla\cdot\widetilde{\bm{u}}^{n+1},\quad\text{in}\;\Omega,\\ \frac{\partial p^{n+1}}{\partial\bm{n}}&=0,\quad\text{on}\;\partial\Omega.\end{aligned}\right.

Based on the previous finite element spaces, we derive the fully discrete scheme of projection method.
(Scheme 2 - Fully-discrete of classical projection scheme) Let $\bm{u}^{0}_{h}$ be the Lagrangian interpolation or Ritz projection onto $\bm{V}_{h}$ of $\bm{u}_{0}$ , then for $n=0,1,\ldots,N-1$ , we can compute $(\widetilde{\bm{u}}^{n+1}_{h},\bm{u}^{n+1}_{h},p^{n+1}_{h})\in$ $\bm{V}_{h}\times\bm{W}_{h}\times Q_{h}$ iteratively by

\left\{\begin{aligned} \left(\frac{\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n}_{h}}{\Delta t},\bm{v}_{h}\right)+\nu(\nabla\widetilde{\bm{u}}^{n+1}_{h},\nabla\bm{v}_{h})&=\langle\bm{f}^{n+1},\bm{v}_{h}\rangle,\quad\forall\bm{v}_{h}\in\bm{V}_{h},\\ (\nabla\cdot\widetilde{\bm{u}}^{n+1}_{h},q_{h})+\Delta t(\nabla p^{n+1}_{h},\nabla q_{h})&=0,\quad\forall q_{h}\in Q_{h},\\ \bm{u}^{n+1}_{h}&=\widetilde{\bm{u}}^{n+1}_{h}-\Delta t\cdot\nabla p^{n+1}_{h}.\end{aligned}\right.

(3.9)

As pointed in [25], it is not necessary to compute the projected velocities $\{\bm{u}^{n}_{h}\}^{N}_{n=0}$ since these quantities can be eliminated. To see this, replacing $\bm{u}^{n}_{h}$ in the first equation with the third equation at $t_{n}$ in Scheme 3, we reach the following Proj-PSPG-FE scheme.
(Scheme 3 - Proj-PSPG-FE) Let $\bm{u}^{0}_{h}$ be the Lagrangian interpolation or Ritz projection onto $\bm{V}_{h}$ of $\bm{u}_{0}$ and $p^{0}_{h}$ denote some proper approximation of initial pressure $p(t_{0})$ , we set $\widetilde{\bm{u}}^{0}_{h}=\bm{u}^{0}_{h}+\Delta t\nabla p^{0}_{h}$ and replace $\bm{u}^{n}_{h}$ with the third equality in Scheme 2, we obtain the following equivalent scheme: for $n=0,1,\ldots,N-1$ , find $(\widetilde{\bm{u}}^{n+1}_{h},p^{n+1}_{h})\in\bm{V}_{h}\times Q_{h}$

\left\{\begin{aligned} \left(\frac{\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n}_{h}}{\Delta t},\bm{v}_{h}\right)+\nu(\nabla\widetilde{\bm{u}}^{n+1}_{h},\nabla\bm{v}_{h})+(\nabla p^{n}_{h},\bm{v}_{h})&=\langle\bm{f}^{n+1},\bm{v}_{h}\rangle,\quad\forall\bm{v}_{h}\in\bm{V}_{h},\\ (\nabla\cdot\widetilde{\bm{u}}^{n+1}_{h},q_{h})+\Delta t(\nabla p^{n+1}_{h},\nabla q_{h})&=0,\quad\forall q_{h}\in Q_{h}.\end{aligned}\right.

(3.10)

Remark 3.2.

For the initial discrete pressure $p^{0}_{h}$ , we simply take $p^{0}_{h}=0$ in later numerical experiments, the reasons behind this are twofold: firstly, it can simplify the solving of initial discrete pressure $p^{0}_{h}=0$ , compared to other techniques, like pressure Poisson equation used in [16, 17]; secondly, it can also directly get the initial intermediate discrete velocity $\widetilde{\bm{u}}^{0}_{h}$ since at this point $\widetilde{\bm{u}}^{0}_{h}=\bm{u}^{0}_{h}$ . We may also observe that this rude treatment would make the first few time steps discrete pressure appear oscillation. Indeed, numerical experiments have showed that the first few time steps discrete pressure values have to be abandoned because of big discrete errors; nevertheless, we will see in later numerical experiments that, only after a few time steps(for example, from $n=6$ on), the discrete pressure values can be adopted as the appropriate discrete approximation of the exact pressure.

Remark 3.3.

We note that although the obtained finite element velocity solution in scheme 3 is $\widetilde{\bm{u}}^{n+1}_{h}$ , rather than $\bm{u}^{n+1}_{h}$ , we will see in latter numerical experiment that when choosing $\bm{V}_{h}\times Q_{h}=P^{1}$ - $P^{1}$ , the solved $\widetilde{\bm{u}}^{n+1}_{h}$ can be used as a proper finite element approximation for the exact velocity $\bm{u}^{n+1}$ . If one wants to obtain $\bm{u}^{n+1}_{h}$ , we can multiply by the last equation in scheme 2 some test function $\bm{w}_{h}\in\bm{W}_{h}=\bm{Y}_{h}^{1}\cap\bm{J}_{1}$ to get

(\bm{u}^{n+1}_{h},\bm{w}_{h})=(\widetilde{\bm{u}}^{n+1}_{h},\bm{w}_{h})-\Delta t(\nabla p^{n+1}_{h},\bm{w}_{h}).

(3.11)

then, solving above equation can finally obtain $\bm{u}^{n+1}_{h}\in\bm{W}_{h}=\bm{Y}_{h}^{1}\cap\bm{J}_{1}$ .

We observe that, in the classical projection scheme (3.9) or (3.10), if we set $\Delta t=\mathcal{O}(h^{2})$ , the existence of PSPG-stabilized term $\Delta t(\nabla p^{n+1}_{h},\nabla q_{h})$ makes the spatial convergence order of the error can not achieve more than second order in $L^{2}$ norm of the velocity and first order of the error in the $L^{2}$ norm of the pressure, so $P^{1}$ - $P^{1}$ for finite element spaces and $(\bm{u},p)\in\bm{H}^{2}\times H^{1}$ seem to be the best choice concerning about the efficiency of computation and the regularity for the exact solutions. Thus, We set $l=1$ in the definition of $Y^{l}_{h}$ , thus the finite element spaces for velocity and pressure would be $\bm{V}_{h}=Y^{1}_{h}\cap\bm{V}$ and $Q_{h}=Y^{1}_{h}\cap Q$ , respectively.
In the sequel we assume $\Delta t=\mathcal{O}(h^{2})$ . Specifically, we assume there exist two positive constants $C_{1},C_{2}$ , such that

C_{1}h^{2}\leq\Delta t\leq C_{2}h^{2}.

(3.12)

The following modified inf-sup condition relaxes the classical inf-sup condition and makes many mixed finite element pairs “stable” in the sense of this modified one, whose detailed proof can be found in [4, Lemma 3] or [17, Lemma 2.1].

Lemma 3.1.

(modified inf-sup condition) Assuming (3.12) holds, then for $n=0,1,\ldots,N-1$ , we have the following pressure stability

\beta\|q_{h}\|_{0}\leq Ch\|\nabla q_{h}\|_{0}+\sup_{\bm{v}_{h}\in\bm{V}_{h}}\frac{(q_{h},\nabla\cdot\bm{v}_{h})}{\|\bm{v}_{h}\|_{1}},\quad\forall q_{h}\in Q_{h}.

(3.13)

3.2 convergence of classical projection scheme

In this subsection, we will cite some error estimation results of the finite element solution in the classical projection scheme and omit those proofs, since those estimates have been analyzed in [10] in detail. The necessity of citing is based on the fact that the following POD snapshots are the finite element solutions, which means that when we analyze the discretization error of the POD-based reduced-order solutions, apart from the ROM truncation error, the discretization error also involve spatial and temporal discretization errors.

Lemma 3.2.

(error estimates for Proj-PSPG-FE) For $n=1,\cdots,N$ , Let $(\bm{u}^{n},p^{n})\in\bm{X}\times Q$ be the solution of continuous variational form (2.3) at discrete time $t=t_{n}$ , $(\widetilde{\bm{u}}^{n}_{h},p^{n}_{h})\in\bm{X}_{h}\times Q_{h}$ is the solution obtained with the PSPG-like projection scheme (3.10), and assuming (3.12) holds, i.e. $\Delta t=\mathcal{O}(h^{2})$ , then we have

	$\displaystyle\\|\bm{u}^{n}-\widetilde{\bm{u}}^{n}_{h}\\|_{0}+h\\|p^{n}-p^{n}_{h}\\|_{0}+h\sqrt{\Delta t}\\|\nabla(p^{n}-p^{n}_{h})\\|_{0}$	$\displaystyle\leq C(h^{2}+\Delta t).$		(3.14)
	$\displaystyle\\|\nabla(\bm{u}-\bm{\widetilde{u}}_{h})\\|_{l^{2}(L^{2})}+\\|p-p_{h}\\|_{l^{2}(L^{2})}$	$\displaystyle\leq C(h+\sqrt{\Delta t}).$		(3.14)

In addition to the error estimates about the intermediate velocity $\widetilde{\bm{u}}^{n}_{h}$ analyzed in [10], we can also prove the $L^{2}$ error estimates about the end-of-step velocity $\bm{u}^{n}_{h}$ which is not made in [10]. In other words, we can obtain

Lemma 3.3.

For $n=1,\cdots,N$ , $\bm{u}^{n}_{h}$ is the velocity solution obtained in Scheme 2, and $\bm{u}^{n}$ is the velocity solution of continuous variational form (2.3) at discrete time $t=t_{n}$ , then

\|\bm{u}^{n}-\bm{u}^{n}_{h}\|_{0}\leq C(h^{2}+\Delta t).

(3.15)

Proof.

According to the convergent results given by Lemma 3.2, we can prove the stability of $\|\nabla p^{n}_{h}\|_{0}$ . That is, the first inequality in (3.14) implies

\|\nabla(p^{n}-p^{n}_{h})\|_{0}\leq C(h/\sqrt{\Delta t}+\sqrt{\Delta t}/h)\leq C,

where we have used (3.12). Then

\|\nabla p^{n}_{h}\|_{0}\leq\|\nabla(p^{n}-p^{n}_{h})\|_{0}+\|\nabla p^{n}\|_{0}\leq C,

where the constant $C$ can be chose to stay nearly fixed with the decrease of $h$ and $\Delta t$ . So, based on the result, testing $\bm{w}_{h}=\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n+1}_{h}$ in (3.11) and rearranging to get

\|\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n+1}_{h}\|_{0}\leq\Delta t\|\nabla p^{n+1}_{h}\|_{0}\leq C\Delta t.

Finally, utilizing triangle inequality and Lemma 3.2,

\|\bm{u}^{n+1}-\bm{u}^{n+1}_{h}\|_{0}\leq\|\bm{u}^{n+1}-\widetilde{\bm{u}}^{n+1}_{h}\|_{0}+\|\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n+1}_{h}\|_{0}\leq C(h^{2}+\Delta t).

∎

4 Projection POD

4.1 POD

In this subsection, we will briefly give some essential ingredients for constructing POD method, more details for constructing POD space and POD-based reduced-order models can be found in [21, 32]. As we remarked in Remark 3.2, the numerical pressure oscillation in the first few time steps caused by setting $p^{0}_{h}=0$ makes us have to discard the finite element solutions in the previous time steps, so when using finite element solution obtained by solving Scheme 4 to construct POD spaces, we also need to consider this situations. In view of the fact that the different choice of initial pressure values will affect the discrete error of finite element pressure in the previous time steps, we uniformly denote $n=n_{0},n_{0}\geq 1$ as the starting time steps of taking snapshots in Scheme 4; therefore, for a given positive integer $M$ , we choose the finite element solution $(\widetilde{\bm{u}}^{i}_{h},p^{i}_{h}),\;n_{0}\leq i\leq n_{0}+M-1$ in Scheme 4, and its difference quotients(DQs) $(\partial\widetilde{\bm{u}}^{i}_{h},\partial p^{i}_{h}),\;n_{0}+1\leq i\leq n_{0}+M-1$ , where the DQs $\partial f^{i}$ are defined by $\partial f^{i}:=(f^{i}-f^{i-1})/\Delta t$ for some function $f$ at discrete time $t=t_{i}$ , to formulate the snapshot spaces

	$\displaystyle\widetilde{\mathcal{U}}$	$\displaystyle:=\langle\widetilde{\bm{u}}^{n_{0}}_{h},\widetilde{\bm{u}}^{n_{0}+1}_{h},\cdots,\widetilde{\bm{u}}^{M+n_{0}-1}_{h},\partial\widetilde{\bm{u}}^{n_{0}+1}_{h},\partial\widetilde{\bm{u}}^{n_{0}+2}_{h},\cdots,\partial\widetilde{\bm{u}}^{M+n_{0}-1}_{h}\rangle,$
	$\displaystyle\mathcal{P}$	$\displaystyle:=\langle p^{n_{0}}_{h},p^{n_{0}+1}_{h},\cdots,p^{M+n_{0}-1}_{h},\partial p^{n_{0}+1}_{h},\partial p^{n_{0}+2}_{h},\cdots,\partial p^{M+n_{0}+1}_{h}\rangle.$

where $\langle S\rangle$ denotes the space spanned by the set $S$ and we denote by $N_{s}=2M-1$ the number of snapshots. Let $d_{\widetilde{u}},d_{p}$ be the dimensions of the spaces $\widetilde{\mathcal{U}}$ and $\mathcal{P}$ , respectively. The symbols $K_{\widetilde{u}},K_{p}$ are correlation matrices corresponding to the snapshots $K_{\widetilde{u}}=\left((K^{\widetilde{u}}_{i,j})\right)\in\mathbb{R}^{N_{s}\times N_{s}}$ and $K_{p}=\left((K^{p}_{i,j})\right)\in\mathbb{R}^{N_{s}\times N_{s}}$ where

K^{\widetilde{u}}_{i,j}:=\left(\widetilde{\bm{u}}^{i}_{h},\widetilde{\bm{u}}^{j}_{h}\right),\quad K^{p}_{i,j}:=\left(p^{i}_{h},p^{j}_{h}\right).

and $(\cdot,\cdot)$ is the $L^{2}$ inner product. Just as in [19], a singular value decomposition(SVD) is carried out and the leading generalized eigenfunctions are chosen as bases, referred to as the POD bases. We denote by $\lambda_{1}\geq\lambda_{2}\geq\cdots\lambda_{d_{\widetilde{u}}}>0$ the positive eigenvalues of $K_{\widetilde{u}}$ and by $\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{d_{\widetilde{u}}}\in\mathbb{R}^{N_{s}}$ the associated eigenvectors. Analogously, $\{\gamma_{i},\bm{y}_{i}\}^{d_{p}}_{i=1}$ are the eigen-pairs of $K_{p}$ . Then, the two POD bases can be written explicitly as

	$\displaystyle\bm{\varphi}_{i}$	$\displaystyle=\frac{1}{\sqrt{\lambda_{i}}}\left(\sum^{n_{0}+M-1}_{j=n_{0}}x^{j}_{i}\widetilde{\bm{u}}^{j}_{h}+\sum^{n_{0}+M-1}_{j=n_{0}+1}x^{j+M-1}_{i}\partial\widetilde{\bm{u}}^{j}_{h}\right),\;$	$\displaystyle\text{for}\;i=1,2,\cdots,d_{\widetilde{u}},$		(4.16)
	$\displaystyle\psi_{i}$	$\displaystyle=\frac{1}{\sqrt{\gamma_{i}}}\left(\sum^{n_{0}+M-1}_{j=n_{0}}y^{j}_{i}p^{j}_{h}+\sum^{n_{0}+M-1}_{j=n_{0}}y^{j+M-1}_{i}\partial p^{j}_{h}\right),\;$	$\displaystyle\text{for}\;i=1,2,\cdots,d_{p}.$		(4.16)

In what follows we will denote by

\displaystyle\bm{V}_{r}=\langle\bm{\varphi}_{1},\bm{\varphi}_{2},\cdots,\bm{\varphi}_{r}\rangle,\quad Q_{r}=\langle\psi_{1},\psi_{2},\cdots,\psi_{r}\rangle.

i.e., we choose the first $r$ POD bases to span the POD spaces, and define the traditional $L^{2}$ projection into the POD velocity space $\bm{V}_{r}$ and POD pressure space $Q_{r}$ as

Definition 4.1.

Let $\Pi^{v}_{r}:L^{2}(\Omega)\rightarrow\bm{V}_{r}$ and $\Pi^{q}_{r}:L^{2}(\Omega)\rightarrow Q_{r}$ such that

	$\displaystyle(\bm{u}-\Pi^{v}_{r}\bm{u},\bm{v}_{r})$	$\displaystyle=0\quad\forall\bm{v}_{r}\in\bm{V}_{r}\quad and$
	$\displaystyle(p-\Pi^{q}_{r}p,q_{r})$	$\displaystyle=0\quad\forall q_{r}\in Q_{r}.$

Then, we can get the following Proj-POD-ROM scheme:
(Scheme 4 - Proj-POD-ROM) For some $n_{0}\geq 1$ , let $(\widetilde{\bm{u}}^{n_{0}}_{r},p^{n_{0}}_{r})=(\Pi^{v}_{r}\widetilde{\bm{u}}^{n_{0}}_{h},\Pi^{q}_{r}p^{n_{0}}_{h})=(\sum^{r}_{i=1}(\widetilde{\bm{u}}^{n_{0}}_{h},\bm{\varphi}_{i})\bm{\varphi}_{i}$ , $\sum^{r}_{i=1}(p^{n_{0}}_{h},\psi_{i})\psi_{i})$ be the Galerkin projection of $(\widetilde{\bm{u}}^{n_{0}}_{h},p^{n_{0}}_{h})$ onto POD spaces $\bm{V}_{r}\times Q_{r}$ , then for $n=n_{0},n_{0}+1,n_{0}+2,\cdots,N-1$ , we can compute $(\widetilde{\bm{u}}^{n+1}_{r},p^{n+1}_{r})\in\bm{V}_{r}\times Q_{r}$ from

\left\{\begin{aligned} \left(\frac{\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}}{\Delta t},\bm{v}_{r}\right)+\nu(\nabla\widetilde{\bm{u}}^{n+1}_{r},\nabla\bm{v}_{r})+(\nabla p^{n}_{r},\bm{v}_{r})&=\langle\bm{f}^{n+1},\bm{v}_{r}\rangle,\quad\forall\bm{v}_{r}\in\bm{V}_{r},\\ (\nabla\cdot\widetilde{\bm{u}}^{n+1}_{r},q_{r})+\Delta t(\nabla p^{n+1}_{r},\nabla q_{r})&=0,\quad\forall q_{r}\in Q_{r}.\end{aligned}\right.

(4.17)

4.2 stability and convergence analysis of projection-POD scheme

We will carry on some numerical analysis about the newly proposed Proj-POD-ROM. To this end, we first recall the classical conclusions concerning about the projection error between the snapshots solution and POD-based reduced-order solution(see [19])

Lemma 4.1.

(time-averaged projection error estimates) We have the following time-averaged projection error estimates.

	$\displaystyle\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}}\left\\|\widetilde{\bm{u}}^{j}_{h}-\sum^{r}_{k=1}\left(\widetilde{\bm{u}}^{j}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\\|^{2}_{0}+\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}+1}\left\\|\partial\widetilde{\bm{u}}^{j}_{h}-\sum^{r}_{k=1}\left(\partial\widetilde{\bm{u}}^{j}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\\|^{2}_{0}$	$\displaystyle=\sum^{d_{\widetilde{u}}}_{k=r+1}\lambda_{k},$		(4.18)
	$\displaystyle\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}}\left\\|p^{j}_{h}-\sum^{r}_{k=1}\left(p^{j}_{h},\psi_{k}\right)\phi_{k}\right\\|^{2}_{0}+\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}+1}\left\\|\partial p^{j}_{h}-\sum^{r}_{k=1}\left(\partial p^{j}_{h},\psi_{k}\right)\phi_{k}\right\\|^{2}_{0}$	$\displaystyle=\sum^{d_{p}}_{k=r+1}\gamma_{k}.$		(4.18)

Remark 4.1.

We remark that choosing $L^{2}$ inner product to construct POD bases and by virtue of Parseval’s identity makes the above three identities actually represent $L^{2}$ POD projection errors.

Apart from the above time-averaged projection error estimates, we also need the following up-to-data optimal in time error estimates(see, [18]), which is presented in the form of hypothesis in the previous papers and is proved to be valid in the presence of DQs in [18].

Lemma 4.2.

(POD Optimal pointwise in time error estimate) We have the following optimal in time projection errors

	$\displaystyle\max_{0\leq k\leq N}\\|\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\bm{\widetilde{u}}^{k}_{h}\\|^{2}_{0}$	$\displaystyle\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i},$
	$\displaystyle\max_{0\leq k\leq N}\\|p^{k}_{h}-\Pi^{q}_{r}p^{k}_{h}\\|^{2}_{0}$	$\displaystyle\leq C\sum^{d_{p}}_{i=r+1}\gamma_{i},$

where $C=6\max\{1,T^{2}\}$ , and

	$\displaystyle\max_{0\leq k\leq N}\\|\nabla(\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\bm{\widetilde{u}}^{k}_{h})\\|^{2}_{0}$	$\displaystyle\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i}\\|\nabla\bm{\varphi}_{i}\\|^{2}_{0},$
	$\displaystyle\max_{0\leq k\leq N}\\|\nabla(p^{k}_{h}-\Pi^{q}_{r}p^{k}_{h})\\|^{2}_{0}$	$\displaystyle\leq C\sum^{d_{p}}_{i=r+1}\gamma_{i}\\|\nabla\psi_{i}\\|^{2}_{0},$

where $C=6\max\{1,T^{2}\}$ .

We also need in the later analysis the time-averaged optimal projection error estimate about DQs, which is derived in [19].

Lemma 4.3.

(POD Optimal time-averaged error estimate about DQs) For the difference quotients we have the estimate:

\frac{1}{M}\sum^{n_{0}+M-1}_{k=n_{0}}\|\partial\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\partial\bm{\widetilde{u}}^{k}_{h}\|^{2}_{0}\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i}.

The following analogues of finite element inverse inequalities will be used in later analysis(see, [19])

Lemma 4.4.

(POD inverse estimates) Denote by $S_{\widetilde{u}}=\left((s^{\widetilde{u}}_{i,j})\right)\in\mathbb{R}^{{N_{s}}\times{N_{s}}}$ and $S_{p}=\left((s^{p}_{i,j})\right)\in\mathbb{R}^{{N_{s}}\times{N_{s}}}$ are stiffness matrices for POD spaces $\bm{V}_{r},Q_{r}$ respectively, with entries $s^{\widetilde{u}}_{i,j}=(\nabla\bm{\varphi}_{i},\nabla\bm{\varphi}_{j}),s^{p}_{i,j}=(\nabla\psi_{i},\nabla\psi_{j})$ , and by ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$ the spectral norm of the matrix, then

	$\displaystyle\\|\nabla\widetilde{\bm{v}}\\|_{0}$	$\displaystyle\leq\sqrt{{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{S}_{\widetilde{u}}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}_{2}}\\|\widetilde{\bm{v}}\\|_{0},\ \text{for}\ \forall\widetilde{\bm{v}}\in\bm{V}_{r},$		(4.19)
	$\displaystyle\\|\nabla q\\|_{0}$	$\displaystyle\leq\sqrt{{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|S_{p}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}_{2}}\\|q\\|_{0},\ \text{for}\;\forall q\in Q_{r}.$		(4.19)

Remark 4.2.

Mass matrix $M$ , for both $(\bm{\varphi}_{i},\bm{\varphi}_{j})$ and $(\psi_{i},\psi_{j})$ , equals identity $I$ since the standard orthogonality of POD base in the sense of $L^{2}$ inner product.

We can prove the unconditional stability of the Proj-POD-ROM scheme.

Theorem 4.1.

(Unconditional stability of Proj-POD-ROM) We have the following unconditional stability

	$\displaystyle\left\\|\widetilde{\bm{u}}_{r}^{N}\right\\|_{0}^{2}+\sum_{n=n_{0}}^{N-1}\left\\|\widetilde{\bm{u}}_{r}^{n+1}-\widetilde{\bm{u}}_{r}^{n}\right\\|_{0}^{2}$	$\displaystyle+\Delta t\sum_{n=n_{0}}^{N-1}\left(\nu\left\\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}+\Delta t\left\\|\nabla p_{r}^{n+1}\right\\|_{0}^{2},\right)$		(4.20)
	$\displaystyle\leq$	$\displaystyle C\left(\left\\|\widetilde{\bm{u}}_{r}^{0}\right\\|_{0}^{2}+\frac{\Delta t}{\nu}\sum_{n=n_{0}}^{N-1}\\|\bm{f}^{n+1}\\|_{-1}\right).$		(4.20)

Proof.

Taking $(\bm{v}_{r},q_{r})=(\widetilde{\bm{u}}^{n+1}_{r},p^{n}_{r})$ in (4.17) and multiply by $2\Delta t$ , adding both equations and integrating by parts, adding and subtracting $2\Delta t^{2}\|\nabla p_{r}^{n+1}\|^{2}_{0}$ , we get

	$\displaystyle\left\\|\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}-\left\\|\widetilde{\bm{u}}_{r}^{n}\right\\|_{0}^{2}$	$\displaystyle+\left\\|\widetilde{\bm{u}}_{r}^{n+1}-\widetilde{\bm{u}}_{r}^{n}\right\\|_{0}^{2}+2\nu\Delta t\left\\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}+2\Delta t^{2}\\|\nabla p_{r}^{n+1}\\|^{2}_{0},$		(4.21)
		$\displaystyle\leq C\frac{\Delta t}{\nu}\\|\bm{f}^{n+1}\\|_{-1}+\nu\Delta t\left\\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}+2\Delta t^{2}\left(\nabla p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right).$		(4.21)

For the last term above, we utilize the second equation in (4.17) twice to obtain

\Delta t\left(\nabla p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right)=-(\nabla\cdot(\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}),p_{r}^{n+1})=(\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r},\nabla p_{r}^{n+1}).

This equality implies

2\Delta t^{2}\left(\nabla p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right)\leq\frac{2}{3}\|\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}\|^{2}_{0}+\frac{3}{2}\Delta t^{2}\|\nabla p_{r}^{n+1}\|^{2}_{0}.

Inserting the above equality into (4.21) and adding $n$ from $0$ to $N-1$ we finally obtain (4.20). ∎

We have the following truncation error estimate between the POD reduced-order solution and the snapshot data.

Lemma 4.5.

(POD truncation error) For $n=1,2,\cdots,N$ , $(\widetilde{\bm{u}}^{n}_{h},p^{n}_{h})$ is the solution of finite element projection scheme (3.10) and $(\widetilde{\bm{u}}^{n}_{r},p^{n}_{r})$ is the solution of POD projection scheme (4.17), then we have the following error estimate:

	$\displaystyle\\|\widetilde{\bm{u}}^{N}_{h}-\widetilde{\bm{u}}^{N}_{r}\\|^{2}_{0}$	$\displaystyle+\nu\Delta t\sum^{N-1}_{n=n_{0}}\\|\nabla(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r})\\|^{2}_{0}+\Delta t^{2}\sum^{N-1}_{n=n_{0}}\\|\nabla(p^{n+1}_{h}-p^{n+1}_{r})\\|^{2}_{0},$		(4.22)
		$\displaystyle\leq C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\\|\nabla\varphi_{i}\\|^{2}_{0}\Big{)}+(\Delta t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+\Delta t^{2}\right).$		(4.22)

Proof.

Finite element projection scheme (3.10) subtract POD projection scheme (4.17) to get:

\left\{\begin{aligned} \left(\frac{\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n}_{h}}{\Delta t}-\frac{\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}}{\Delta t},\bm{v}_{r}\right)+\nu(\nabla(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}),\nabla\bm{v}_{r})+(\nabla(p^{n}_{h}-p^{n}_{r}),\bm{v}_{r})&=0,\\ (\nabla\cdot(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}),q_{r})+\Delta t(\nabla(p^{n+1}_{h}-p^{n+1}_{r}),\nabla q_{r})&=0.\end{aligned}\right.

(4.23)

Denoting

	$\displaystyle\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}=\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\widetilde{\bm{u}}^{n+1}_{h}+\Pi^{v}_{r}\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}$	$\displaystyle:=\widetilde{\bm{\eta}}^{n+1}_{u,r}+\widetilde{\bm{w}}^{n+1}_{u,r},$
	$\displaystyle p^{n+1}_{h}-p^{n+1}_{r}=p^{n+1}_{h}-\Pi^{q}_{r}p^{n+1}_{h}+\Pi^{q}_{r}p^{n+1}_{h}-p^{n+1}_{r}$	$\displaystyle:=\eta^{n+1}_{p,r}+w^{n+1}_{p,r}.$

Testing $(\bm{v}_{r},q_{r})=(\widetilde{\bm{w}}^{n+1}_{u,r},w^{n+1}_{p,r})$ , we obtain

\left\{\begin{aligned} \left(\frac{\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}}{\Delta t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)&+\nu\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}+(\nabla w^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r})=\left(\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)+\nu(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\nabla\widetilde{\bm{w}}^{n+1}_{u,r}),\\ &\quad+(\nabla(p^{n+1}_{h}-p^{n}_{h}),\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla(p^{n}_{r}-p^{n+1}_{r}),\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla\eta^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r}),\\ (\nabla\cdot\widetilde{\bm{w}}^{n+1}_{u,r},w^{n+1}_{p,r})&+\Delta t\|\nabla w^{n+1}_{p,r}\|^{2}_{0}=(\nabla\cdot\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})+\Delta t(\nabla\eta^{n+1}_{p,r},\nabla w^{n+1}_{p,r}).\end{aligned}\right.

Adding both equations we have the following error equation

$\displaystyle\left(\frac{\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}}{\Delta t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)$	$\displaystyle+\nu\\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}+\Delta t\\|\nabla w^{n+1}_{p,r}\\|^{2}_{0}=\left(\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta t},\widetilde{\bm{w}}^{n+1}_{u,r}\right),$	(4.24)
	$\displaystyle+\nu(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\nabla\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla(p^{n+1}_{h}-p^{n}_{h}),\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla(p^{n}_{r}-p^{n+1}_{r}),\widetilde{\bm{w}}^{n+1}_{u,r}),$
	$\displaystyle+(\nabla\eta^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla\cdot\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})+\Delta t(\nabla\eta^{n+1}_{p,r},\nabla w^{n+1}_{p,r}),$
	$\displaystyle:=A_{1}+A_{2}+A_{3}+A_{4}+A_{5}+A_{6}+A_{7}.$

We will bound the seven terms above separately.

\displaystyle|A_{1}|

\displaystyle=\left|\left(\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)\right|\leq C\left\|\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta t}\right\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}.

For the first term above, we utilize the expression of $\widetilde{\bm{u}}^{n}_{h}$ and $\Pi^{v}_{r}\widetilde{\bm{u}}^{n}_{h}$ in the form of POD bases $\bm{\varphi}_{i}$ to get

\displaystyle\widetilde{\bm{\eta}}^{n+1}_{u,r}=\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\widetilde{\bm{u}}^{n+1}_{h}=\sum^{d_{\widetilde{u}}}_{k=1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}-\sum^{r}_{k=1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}=\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}.

Then

	$\displaystyle\left\\|\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta t}\right\\|^{2}_{0}$	$\displaystyle=\left\\|\frac{1}{\Delta t}\left(\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}-\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\widetilde{\bm{u}}^{n}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right)\right\\|^{2}_{0},$
		$\displaystyle=\left\\|\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\frac{\widetilde{\bm{u}}^{n}_{h}-\widetilde{\bm{u}}^{n}_{h}}{\Delta t},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\\|^{2}_{0},$
		$\displaystyle=\left\\|\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\partial\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\\|^{2}_{0},$
		$\displaystyle=\left\\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\\|^{2}_{0}.$

For other terms on the right-side hand of (4.24),

	$\displaystyle\|A_{2}\|$	$\displaystyle=\left\|\nu(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\nabla\widetilde{\bm{w}}^{n+1}_{u,r})\right\|\leq C\nu\\|\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r}\\|^{2}_{0}+\frac{1}{2}\nu\\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}\leq C\nu\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\\|\nabla\varphi_{i}\\|^{2}_{0}+\frac{1}{2}\nu\\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0},$
	$\displaystyle\|A_{3}\|$	$\displaystyle=\left\|(\nabla(p^{n+1}_{h}-p^{n}_{h}),\widetilde{\bm{w}}^{n+1}_{u,r})\right\|\leq C\\|\nabla(p^{n+1}_{h}-p^{n}_{h})\\|^{2}_{0}+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}\leq C\Delta t\int^{t_{n+1}}_{t_{n}}\\|\nabla(\partial_{t}p_{h})\\|^{2}_{0}\,{\rm d}t+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0},$
	$\displaystyle\|A_{4}\|$	$\displaystyle=\left\|(\nabla(p^{n}_{r}-p^{n+1}_{r}),\widetilde{\bm{w}}^{n+1}_{u,r})\right\|\leq C\\|\nabla(p^{n+1}_{r}-p^{n}_{r})\\|^{2}_{0}+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}\leq C\Delta t\int^{t_{n+1}}_{t_{n}}\\|\nabla(\partial_{t}p_{r})\\|^{2}_{0}\,{\rm d}t+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0},$
	$\displaystyle\|A_{5}\|$	$\displaystyle=\left\|(\nabla\eta^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r})\right\|\leq C\\|\nabla\eta^{n+1}_{p,r}\\|^{2}_{0}+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}\leq C\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0},$
	$\displaystyle\|A_{6}\|$	$\displaystyle=\left\|(\nabla\cdot\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})\right\|\leq C\|(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})\|\leq C\\|\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r}\\|^{2}_{0}+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\\|\nabla\varphi\\|^{2}_{0}+\frac{1}{10}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0},$
	$\displaystyle\|A_{7}\|$	$\displaystyle=\left\|\Delta t(\nabla\eta^{n+1}_{p,r},\nabla w^{n+1}_{p,r})\right\|\leq\frac{1}{2}\Delta t\\|\nabla\eta^{n+1}_{p,r}\\|^{2}_{0}+\frac{1}{2}\Delta t\\|\nabla w^{n+1}_{p,r}\\|^{2}_{0}\leq C\Delta t\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+\frac{1}{2}\Delta t\\|\nabla w^{n+1}_{p,r}\\|^{2}_{0}.$

Combining with all the seven terms’ results, we have

	$\displaystyle\left(\frac{\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}}{\Delta t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)$	$\displaystyle+\frac{1}{2}\nu\\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}+\frac{1}{2}\Delta t\\|\nabla w^{n+1}_{p,r}\\|^{2}_{0}\leq\left\\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\\|^{2}_{0},$
		$\displaystyle\quad+C\Delta t\int^{t_{n+1}}_{t_{n}}\Big{(}\\|\partial_{t}p_{h}\\|^{2}_{0}+\\|\partial_{t}p_{r}\\|^{2}_{0}\Big{)}\,{\rm d}t+C(\nu+1)\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\\|\nabla\varphi_{i}\\|^{2}_{0},$
		$\displaystyle\quad+C(\Delta t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+\frac{1}{2}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}.$

Multiplying by $\Delta t$ , adding with $n$ from $n_{0}$ to $N-1$ and rearranging, we get

	$\displaystyle\\|\widetilde{\bm{w}}^{N}_{u,r}\\|^{2}_{0}$	$\displaystyle+\Delta t\sum^{N-1}_{n=n_{0}}\\|\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}\\|^{2}_{0}+\nu\Delta t\sum^{N-1}_{n=n_{0}}\\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}+\Delta t^{2}\sum^{N-1}_{n=n_{0}}\\|\nabla w^{n+1}_{p,r}\\|^{2}_{0},$
		$\displaystyle\quad\leq\Delta t\sum^{N-1}_{n=n_{0}}\left\\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\\|^{2}_{0}+C(\nu+1)\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\\|\nabla\varphi_{i}\\|^{2}_{0}+\frac{1}{2}\Delta t\sum^{N-1}_{n=n_{0}}\\|\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0},$
		$\displaystyle\quad+C(\Delta t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+C\Delta t^{2}\int^{T}_{t_{n_{0}}}\Big{(}\\|\partial_{t}p_{h}\\|^{2}_{0}+\\|\partial_{t}p_{r}\\|^{2}_{0}\Big{)}\,{\rm d}t.$

For the first term on the right-hand side above, we use Lemma 4.3 to have

\displaystyle\Delta t\sum^{N-1}_{n=n_{0}}\left\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\|^{2}_{0}\leq C\frac{1}{M}\sum^{M}_{k=1}\|\partial\widetilde{\bm{u}}^{k}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{k}_{h}\|^{2}_{0}\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}.

By discrete Gronwall inequality and using stability result to bound the last term above, we get

	$\displaystyle\\|\widetilde{\bm{w}}^{N}_{u,r}\\|^{2}_{0}$	$\displaystyle+\Delta t\sum^{N-1}_{n=n_{0}}\\|\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}\\|^{2}_{0}+\nu\Delta t\sum^{N-1}_{n=n_{0}}\\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\\|^{2}_{0}+\Delta t^{2}\sum^{N-1}_{n=n_{0}}\\|\nabla w^{n+1}_{p,r}\\|^{2}_{0},$		(4.25)
		$\displaystyle\quad\leq C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\\|\nabla\varphi_{i}\\|^{2}_{0}\Big{)}+(\Delta t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+\Delta t^{2}\right).$		(4.25)

∎

Finally, by triangle inequality, Lemma 3.2 and Lemma 4.5, we obtain the following theorem which states the convergence between continuous variational form and POD projection scheme.

Theorem 4.2.

(error estimate for Proj-POD-ROM) For $n=n_{0},n_{0}+1,\cdots,N$ , let $(\bm{u}^{n},p^{n})$ is the solution of (2.3) at $t=t_{n}$ , $(\widetilde{\bm{u}}^{n}_{r},p^{n}_{r})$ denotes the POD-based reduced-order solutions obtained in (4.17), then we have the following convergence estimate:

	$\displaystyle\left\\|\bm{u}^{N}-\widetilde{\bm{u}}^{N}_{r}\right\\|^{2}_{0}$	$\displaystyle+\nu\Delta t\sum^{N-1}_{n=n_{0}}\left\\|\nabla(\widetilde{\bm{u}}^{n}-\widetilde{\bm{u}}^{n}_{r})\right\\|^{2}_{L^{2}}+\Delta t^{2}\sum^{N-1}_{n=n_{0}}\left\\|\nabla(p-p_{r})\right\\|^{2}_{L^{2}},$		(4.26)
		$\displaystyle\leq C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\left\\|\nabla\varphi_{i}\right\\|^{2}_{0}\Big{)}+(\Delta t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\left\\|\nabla\psi_{i}\right\\|^{2}_{0}+\Delta t^{2}+h^{2}\right).$		(4.26)

5 Numerical Tests

In this section, we present some numerical experiments to confirm the $a\,priori$ error estimates derived in Theorem 4.2 for the POD projection scheme 4 - Proj-POD-ROM. To this end, we first determine some necessary parameters used before and verify the convergence in Lemma 3.2 and Lemma 3.3 of the finite element projection scheme, since with which we formed our POD snapshot spaces, and also the convergence result in Theorem 4.2 is closely related to the one of the finite element scheme. The open-source finite element package iFEM [6] has been used to run the numerical experiments.

5.1 Problem setting

We take $\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$ and the time interval $(0,1]$ with the viscosity coefficient $\nu=1$ and the prescribed solution

	$\displaystyle\bm{u}$	$\displaystyle=\cos(t)\cdot\left(\begin{array}[]{c}\pi\sin(\pi x)^{2}\sin(2\pi y)\\ -\pi\sin(2\pi x)\sin(\pi y)^{2}\\ \end{array}\right),$
	$\displaystyle p$	$\displaystyle=\cos(t)\cdot 10\cos(\pi x)\cos(\pi y).$

The right-hand side, the Dirichlet boundary condition and the initial condition are chosen in accordance to the above prescribed solution.
We set $\Delta t=0.1h^{2}$ and all grids are regular $N\times N$ triangular grids with SWNE diagonals for different $N$ (i.e., diagonals coming from connecting the southwest and the northeast vertexes on all rectangles), and we take $N=4,8,16,32,64$ sequentially. For simplicity, we use $P^{1}$ - $P^{1}$ element pair for spatial discretization, and we take the snapshots on the finest computational mesh, i.e., h=1/64, and show this mesh in Figure 1.
When constructing the snapshots spaces, we take into account the fact that the finite element projection scheme involves the initial pressure $p^{0}_{h}$ which is not the part of the definition of the problem, so we take $p^{0}_{h}=0$ and it might be good to think of the numerical discrete errors of pressure $\|p^{n}-p^{n}_{h}\|_{L^{2}}$ in the previous steps $n$ are so large that it is not suitable to take finite element solution on those steps into snapshot spaces, thus we choose the finite element solution $(\bm{u}^{n}_{h},p^{n}_{h})$ and its difference quotients $(\partial\bm{u}^{n}_{h},\partial p^{n}_{h})$ from $n=n_{0}=6$ to $n=25$ to formulate snapshots spaces, which means we take $M=20$ and thus the number of snapshots is $N_{s}=2M-1=39$ .
The finite element solution $(\bm{u}^{n}_{h},\widetilde{\bm{u}}^{n}_{j},p^{n}_{h})$ and exact solution $(\bm{u}^{n},p^{n})$ at discrete termination time $t=\Delta t\cdot N$ on finest mesh are shown in Figure 2, and following POD bases are formed from snapshots via $L^{2}$ inner product.

Refer to caption — Figure 1: Computational mesh with h=1/64.

5.2 Convergence of the finite element projection scheme

In this subsection, numerical tests will be used to numerically verify two theoretical results, that is Lemma 3.2 and Lemma 3.3, with which we utilize to get convergence analysis Theorem 4.2 of POD scheme 4. Although Lemma 3.2 has been analyzed in [10], it lacks the corresponding numerical verifications, so we supplement the numerical tests about Lemma 3.2 and also provide that of the Lemma 3.3, which is newly analyzed in this paper. Since $\Delta t=\mathcal{O}(h^{2})$ , then from the Table 1 and Table 2 we can see

	$\displaystyle\\|\nabla(\bm{u}-\bm{\widetilde{u}}_{h})\\|_{l^{2}(L^{2})}+\\|p-p_{h}\\|_{l^{2}(L^{2})}$	$\displaystyle=\mathcal{O}(h),$
	$\displaystyle\\|\bm{u}^{n}-\widetilde{\bm{u}}^{n}_{h}\\|_{0}+h\\|p^{n}-p^{n}_{h}\\|_{0}+h\sqrt{\Delta t}\\|\nabla(p^{n}-p^{n}_{h})\\|_{0}$	$\displaystyle=\mathcal{O}(h^{2}).$

which are what the theoretical results in Lemma 3.2 shows. Especially, we can see from Table 1, $\left\|\bm{u}^{n}-\bm{u}^{n}_{h}\right\|_{L^{2}}=\mathcal{O}(h^{2})$ , which is in accordance of Lemma 3.3.

Table 1: Errors and convergence orders of finite element solution

\widetilde{\bm{u}}^{n}_{h},\bm{u}^{n}_{h}

with

P^{1}

P^{1}

pair.

1/h	$\max\limits_{n}\left\\|\bm{u}^{n}-\widetilde{\bm{u}}^{n}_{h}\right\\|_{L^{2}}$		$\max\limits_{n}\left\\|\bm{u}^{n}-\bm{u}^{n}_{h}\right\\|_{L^{2}}$		$\left\\|\nabla(\bm{u}-\widetilde{\bm{u}}_{h})\right\\|_{\ell^{2}(L^{2})}$
1/h	error	rate	error	rate	error	rate
4	5.3013e $-$ 01	-	4.6509e $-$ 01	-	4.8187e $+$ 00	-
8	1.6490e $-$ 01	1.7078	1.4931e $-$ 01	1.6392	2.6626e $+$ 00	0.85582
16	4.3368e $-$ 02	1.9259	4.0108e $-$ 02	1.8963	1.3785e $+$ 00	0.94976
32	1.0969e $-$ 02	1.9828	1.0527e $-$ 02	1.9298	7.1098e $-$ 01	0.95516
64	2.7499e $-$ 03	1.9960	2.8273e $-$ 03	1.8966	3.7409e $-$ 01	0.92642

Table 2: Errors and convergence orders of finite element pressure

p^{n}_{h}

with

P^{1}

P^{1}

pair.

1/h	$\max\limits_{n}\\|p^{n}-p^{n}_{h}\\|_{L^{2}}$		$\left\\|p-p_{h}\right\\|_{\ell^{2}(L^{2})}$		$\sqrt{\Delta t}\left\\|\nabla(p-p_{h})\right\\|_{\ell^{2}(L^{2})}$
1/h	error	rate	error	rate	error	rate
4	2.7636e $+$ 00	-	2.2987e $+$ 00	-	1.4286e $+$ 00	-
8	1.1144e $+$ 00	1.3103	8.8892e $-$ 01	1.3707	4.6814e $-$ 01	1.6096
16	3.6664e $-$ 01	1.6038	2.7275e $-$ 01	1.7045	1.3827e $-$ 01	1.7595
32	1.2463e $-$ 01	1.5567	8.1260e $-$ 02	1.7470	4.5158e $-$ 01	1.6144
64	4.6335e $-$ 02	1.4275	2.5152e $-$ 02	1.6919	1.5553e $-$ 01	1.5378

5.3 Convergence of projection POD

After having the snapshots spaces, we determine $d_{\widetilde{u}}=\text{rank}(\widetilde{\mathcal{U}})=20$ (for which $\gamma_{i}<10^{-15},i>d_{\widetilde{u}}$ ), $d_{p}=\text{rank}(\mathcal{P})=20$ (for which $\epsilon_{i}<10^{-13},i>d_{p}$ ). Figure 3 shows the decay of POD eigenvalues(left) of velocity $\gamma_{i},i=1,\cdots,d_{\widetilde{u}}$ and pressure $\epsilon_{i},i=1,\cdots,d_{p}$ , together with the corresponding captured system’s energy(right) in the form of $100\sum^{r}_{i=1}\gamma_{i}/\sum^{d_{\widetilde{u}}}_{i=1}\gamma_{i}$ for velocity and $100\sum^{r}_{i=1}\epsilon_{i}/\sum^{d_{p}}_{i=1}\epsilon_{i}$ for pressure. We note that the first $r=4$ POD modes already capture more than $99.99\%$ of the system’s velocity-pressure energy.
In order to confirm the computational efficiency advantage of the POD-ROM scheme over the finite element FOM scheme, we take the number of POD modes is $r=4$ , whereas the number of degree of freedoms(DOFs) for velocity and pressure in finite element scheme are 8450 and 4225. We show in Table 3 the comparison of cumulative time spent by the two schemes to run to some certain time in the process of time-stepping iterative, and the $L^{2}$ numerical spatial discrete error of velocity and pressure in two schemes at that time. We can see from Table 3 that, compared with the finite element FOM scheme, POD-ROM scheme can effectively improve computational efficiency in the context of only consuming approximately 1/6 times to get the numerical solution with even higher accuracy.
Figure 4 plots the temporal evolution of the discrete $L^{2}$ relative error(in semilogarithmic scale) of the reduced-order velocity and pressure with respect to the full order ones: $\|\widetilde{\bm{u}}^{n}_{h}-\widetilde{\bm{u}}^{n}_{r}\|_{L^{2}}/\|\widetilde{\bm{u}}^{n}_{h}\|_{L^{2}}$ and $\|p^{n}_{h}–p^{n}_{r}\|_{L^{2}}/\|p^{n}_{h}\|_{L^{2}}$ for different number of POD modes. As expected in Theorem 4.2, the errors would decrease as the number of POD modes $r$ increasing. For velocity, the left figure in Figure 4 numerically demonstrate this result, and for pressure, the error would increase slightly when $r$ from $6$ to $8$ , but when we continue to increase $r$ , the error will decrease, which is consistent with the Theorem 4.2 as a whole.

Table 3: Comparison of errors and time consumed by time-stepping iterative between finite element scheme with

P^{1}

P^{1}

pair and POD scheme with

r=4

on the mesh of

h=1/64

and

\Delta t=0.1h^{2}

n	finite element scheme			POD scheme
n	$\\|\bm{u}^{n}\!-\!\widetilde{\bm{u}}^{n}_{h}\\|_{L^{2}}$	$\\|p^{n}\!-\!p^{n}_{h}\\|_{L^{2}}$	CPU run time(s)	$\\|\bm{u}^{n}\!-\!\widetilde{\bm{u}}^{n}_{r}\\|_{L^{2}}$	$\\|p^{n}\!-\!p^{n}_{r}\\|_{L^{2}}$	CPU run time(s)
2500	2.3789e $-$ 03	2.9458e $-$ 02	689	2.2860e $-$ 03	2.7823e $-$ 02	113
5000	2.3929e $-$ 03	2.9253e $-$ 02	1338	2.3007e $-$ 03	2.7642e $-$ 02	227
7500	2.3740e $-$ 03	2.8975e $-$ 02	2005	2.2826e $-$ 03	2.7379e $-$ 02	339
10000	2.3452e $-$ 03	2.8591e $-$ 02	2668	2.2549e $-$ 03	2.7015e $-$ 02	452
20000	2.1443e $-$ 03	2.6010e $-$ 02	5339	2.0618e $-$ 03	2.4573e $-$ 02	914
30000	1.8163e $-$ 03	2.1886e $-$ 02	8060	1.7464e $-$ 03	2.0673e $-$ 02	1399
40000	1.3805e $-$ 03	1.6464e $-$ 02	10746	1.3274e $-$ 03	1.5547e $-$ 02	1900

6 Conclusions

In this paper, we proposed an efficient projection POD-ROM, which combined the advantages of classical projection method and POD technique.
The main contribution of the present paper consisted of two aspects: the first one was high computational efficiency. Through auxiliary intermediate velocity variable, the classical projection method decoupled the velocity variable and pressure variable, meanwhile decoupled the saddle-point system arose from Stokes equations, so one strength of projection method lied in its high computational efficiency, or low computational costs; Furthermore, POD technique was utilized to get the ROM, which have made the newly proposed projection POD-ROM had high computational efficiency. The second contribution was based on the fact that, in the fully discrete scheme of the classical projection method, the original scheme could be rewritten into a PSPG-type pressure stabilization scheme by eliminating the end-of-step velocity, where the pressure stabilized term $\Delta t(\nabla p^{n+1}_{h},\nabla q^{n+1}_{h})$ , $\Delta t=\mathcal{O}(h^{2})$ , was inherent, so that some flexible mixed finite element spaces pairs(for example, $P^{1}$ - $P^{1}$ pair) could be used without considering the classical LBB/inf-sup condition for mixed POD spaces, which was different from other stabilized FE-POD-ROM to overcome LBB/inf-sup condition by adding extra stabilization terms.
Numerical experiments have been conducted to confirm the convergence for both projection finite element scheme and projection POD scheme. We first numerically confirmed the PSPG-type classical projection owned the desired convergence orders consistent with theoretical analysis, and then, finite element solutions were taken as the snapshots to formulate POD bases/modes which are used to get reduced-order velocity and pressure variables, so in theory, the discrete error between exact solutions and reduced-order solutions should converge to zero as increasing the number of POD modes, we numerically verify this result. Apart from the projection POD-ROM convergence, we also conduct the experiment to compare the error and computational time between projection finite element FOM and projection POD-ROM, and the results revealed the fact that projection POD-ROM not only had less discretization error, but also less computational costs, compared with projection finite element FOM.
One future research direction will be the applied of projection POD on nonlinear nonstationary Navier-Stokes equations. Another investigation is validation of fulfillment of LBB/inf-sup condition for mixed POD spaces.

Acknowledgment

The authors thank the anonymous referees for their constructive comments and suggestions which improved the manuscript.

References

[1] Mejdi Azaïez, Tomás Chacón Rebollo, and Samuele Rubino. A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations. J. Comput. Phys., 425:Paper No. 109916, 27, 2021.
[2] Francesco Ballarin, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza. Supremizer stabilization of POD-Galerkin approximation of parametrized steady incompressible Navier-Stokes equations. Internat. J. Numer. Methods Engrg., 102(5):1136–1161, 2015.
[3] Susanne C. Brenner and L. Ridgway Scott. The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics. Springer, New York, third edition, 2008.
[4] Erik Burman and Miguel A. Fernández. Analysis of the PSPG method for the transient Stokes’ problem. Comput. Methods Appl. Mech. Engrg., 200(41-44):2882–2890, 2011.
[5] Alfonso Caiazzo, Traian Iliescu, Volker John, and Swetlana Schyschlowa. A numerical investigation of velocity-pressure reduced order models for incompressible flows. J. Comput. Phys., 259:598–616, 2014.
[6] Long Chen. ifem: an innovative finite element methods package in matlab. Preprint, University of Maryland, 2008.
[7] Alexandre Joel Chorin. Numerical solution of the Navier-Stokes equations. Math. Comp., 22:745–762, 1968.
[8] Alexandre Joel Chorin. On the convergence of discrete approximations to the Navier-Stokes equations. Math. Comp., 23:341–353, 1969.
[9] Philippe G. Ciarlet. The finite element method for elliptic problems, volume 40 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.
[10] Javier de Frutos, Bosco García-Archilla, and Julia Novo. Error analysis of projection methods for non inf-sup stable mixed finite elements. The transient Stokes problem. Appl. Math. Comput., 322:154–173, 2018.
[11] Victor DeCaria, Traian Iliescu, William Layton, Michael McLaughlin, and Michael Schneier. An artificial compression reduced order model. SIAM J. Numer. Anal., 58(1):565–589, 2020.
[12] Guosheng Fu and Zhu Wang. POD-(H)DG method for incompressible flow simulations. J. Sci. Comput., 85(2):Paper No. 24, 20, 2020.
[13] Vivette Girault and Pierre-Arnaud Raviart. Finite element methods for Navier-Stokes equations - Theory and algorithms, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986.
[14] J. L. Guermond, P. Minev, and Jie Shen. An overview of projection methods for incompressible flows. Comput. Methods Appl. Mech. Engrg., 195(44-47):6011–6045, 2006.
[15] J.-L. Guermond and L. Quartapelle. On the approximation of the unsteady Navier-Stokes equations by finite element projection methods. Numer. Math., 80(2):207–238, 1998.
[16] John G. Heywood and Rolf Rannacher. Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization. SIAM J. Numer. Anal., 19(2):275–311, 1982.
[17] Volker John and Julia Novo. Analysis of the pressure stabilized Petrov-Galerkin method for the evolutionary Stokes equations avoiding time step restrictions. SIAM J. Numer. Anal., 53(2):1005–1031, 2015.
[18] Birgul Koc, Samuele Rubino, Michael Schneier, John Singler, and Traian Iliescu. On optimal pointwise in time error bounds and difference quotients for the proper orthogonal decomposition. SIAM J. Numer. Anal., 59(4):2163–2196, 2021.
[19] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition methods for parabolic problems. Numer. Math., 90(1):117–148, 2001.
[20] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J. Numer. Anal., 40(2):492–515, 2002.
[21] Zhendong Luo and Goong Chen. Proper orthogonal decomposition methods for partial differential equations. Mathematics in Science and Engineering. Elsevier/Academic Press, London, 2019.
[22] Zhendong Luo, Jing Chen, I. M. Navon, and Xiaozhong Yang. Mixed finite element formulation and error estimates based on proper orthogonal decomposition for the nonstationary Navier-Stokes equations. SIAM J. Numer. Anal., 47(1):1–19, 2008/09.
[23] Bernd R. Noack, Paul Papas, and Peter A. Monkewitz. The need for a pressure-term representation in empirical galerkin models of incompressible shear flows. J. Fluid Mech., 523:339–365, 2005.
[24] Julia Novo and Samuele Rubino. Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows. SIAM J. Numer. Anal., 59(1):334–369, 2021.
[25] Rolf Rannacher. On Chorin’s projection method for the incompressible Navier-Stokes equations. In The Navier-Stokes equations II—theory and numerical methods (Oberwolfach, 1991), volume 1530 of Lecture Notes in Math., pages 167–183. Springer, Berlin, 1992.
[26] S. S. Ravindran. Error analysis for Galerkin POD approximation of the nonstationary Boussinesq equations. Numer. Methods Partial Differential Equations, 27(6):1639–1665, 2011.
[27] Samuele Rubino. Numerical analysis of a projection-based stabilized POD-ROM for incompressible flows. SIAM J. Numer. Anal., 58(4):2019–2058, 2020.
[28] Jie Shen. On error estimates of projection methods for Navier-Stokes equations: first-order schemes. SIAM J. Numer. Anal., 29(1):57–77, 1992.
[29] Jie Shen. On error estimates of some higher order projection and penalty-projection methods for Navier-Stokes equations. Numer. Math., 62(1):49–73, 1992.
[30] John R. Singler. New POD error expressions, error bounds, and asymptotic results for reduced order models of parabolic PDEs. SIAM J. Numer. Anal., 52(2):852–876, 2014.
[31] Roger Temam. Une méthode d’approximation de la solution des équations de Navier-Stokes. Bull. Soc. Math. France, 96:115–152, 1968.
[32] S. Volkwein. Model reduction using proper orthogonal decomposition. Lecture Notes, Faculty of Mathematics and Statistics, University of Konstanz, 2011.

	$\displaystyle\\|\bm{v}-I_{h}\bm{v}\\|_{s,q,K}$	$\displaystyle\leq Ch_{K}^{t-s}\\|\bm{v}\\|_{t,q,K},$		(2.5)
	$\displaystyle\\|q-J_{h}q\\|_{s,q,K}$	$\displaystyle\leq Ch_{K}^{t-s}\\|q\\|_{t,q,K},$		(2.5)

	$\displaystyle\max_{0\leq k\leq N}\\|\nabla(\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\bm{\widetilde{u}}^{k}_{h})\\|^{2}_{0}$	$\displaystyle\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i}\\|\nabla\bm{\varphi}_{i}\\|^{2}_{0},$
	$\displaystyle\max_{0\leq k\leq N}\\|\nabla(p^{k}_{h}-\Pi^{q}_{r}p^{k}_{h})\\|^{2}_{0}$	$\displaystyle\leq C\sum^{d_{p}}_{i=r+1}\gamma_{i}\\|\nabla\psi_{i}\\|^{2}_{0},$

	$\displaystyle\\|\nabla\widetilde{\bm{v}}\\|_{0}$	$\displaystyle\leq\sqrt{{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{S}_{\widetilde{u}}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}_{2}}\\|\widetilde{\bm{v}}\\|_{0},\ \text{for}\ \forall\widetilde{\bm{v}}\in\bm{V}_{r},$		(4.19)
	$\displaystyle\\|\nabla q\\|_{0}$	$\displaystyle\leq\sqrt{{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|S_{p}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}_{2}}\\|q\\|_{0},\ \text{for}\;\forall q\in Q_{r}.$		(4.19)

	$\displaystyle\left\\|\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}-\left\\|\widetilde{\bm{u}}_{r}^{n}\right\\|_{0}^{2}$	$\displaystyle+\left\\|\widetilde{\bm{u}}_{r}^{n+1}-\widetilde{\bm{u}}_{r}^{n}\right\\|_{0}^{2}+2\nu\Delta t\left\\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}+2\Delta t^{2}\\|\nabla p_{r}^{n+1}\\|^{2}_{0},$		(4.21)
		$\displaystyle\leq C\frac{\Delta t}{\nu}\\|\bm{f}^{n+1}\\|_{-1}+\nu\Delta t\left\\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\\|_{0}^{2}+2\Delta t^{2}\left(\nabla p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right).$		(4.21)

	$\displaystyle\\|\widetilde{\bm{u}}^{N}_{h}-\widetilde{\bm{u}}^{N}_{r}\\|^{2}_{0}$	$\displaystyle+\nu\Delta t\sum^{N-1}_{n=n_{0}}\\|\nabla(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r})\\|^{2}_{0}+\Delta t^{2}\sum^{N-1}_{n=n_{0}}\\|\nabla(p^{n+1}_{h}-p^{n+1}_{r})\\|^{2}_{0},$		(4.22)
		$\displaystyle\leq C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\\|\nabla\varphi_{i}\\|^{2}_{0}\Big{)}+(\Delta t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\\|\nabla\psi_{i}\\|^{2}_{0}+\Delta t^{2}\right).$		(4.22)