Error analysis of proper orthogonal decomposition data assimilation schemes with grad–div stabilization for the Navier–Stokes equations

Journal of Computational and Applied Mathematics 411 (2022) 114246 Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam Error analysis of proper orthogonal decomposition data assimilation schemes with grad–div stabilization for the Navier–Stokes equations Bosco García-Archilla a,1, Julia Novo b,∗,2, Samuele Rubino c,3 aDepartamento de Matemática Aplicada II, Universidad de Sevilla, Spain bDepartamento de Matemáticas, Universidad Autónoma de Madrid, Spain cDepartment EDAN &IMUS, Universidad de Sevilla, Spain article info Article history: Received 11 February 2021 Received in revised form 8 February 2022 MSC: 35Q30 65M12 65M15 65M20 65M60 65M70 76B75 Keywords: Data assimilation Navie–Stokes equations Uniform-in-time error estimates Proper orthogonal decomposition Fully discrete schemes Mixed finite elements methods abstract The error analysis of a proper orthogonal decomposition (POD) data assimilation (DA) scheme for the Navier–Stokes equations is carried out. A grad–div stabilization term is added to the formulation of the POD method. Error bounds with constants independent on inverse powers of the viscosity parameter are derived for the POD algorithm. No upper bounds in the nudging parameter of the data assimilation method are required. Numerical experiments show that, for large values of the nudging parameter, the proposed method rapidly converges to the real solution, and greatly improves the overall accuracy of standard POD schemes up to low viscosities over predictive time intervals. ©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Reduced order models (ROM) are a fairly extensive technique applied in many different fields to reduce the computational cost of direct numerical simulations while keeping enough accurate numerical approximations. Proper Orthogonal Decomposition (POD) method provides the elements (modes) of the reduced basis from a given database (snapshots) which are computed by means of a direct or full order method. Data assimilation refers to a class of techniques that combine experimental data and simulations in order to obtain better predictions in a physical system. There is a vast literature on data assimilation methods (see e.g., [1–5], and the ∗Corresponding author. E-mail addresses: [email protected] (B. García-Archilla), [email protected] (J. Novo), [email protected] (S. Rubino). 1Research is supported by Spanish MCINYU, Spain under grants PGC2018-096265-B-I00 and PID2019-104141GB-I00. 2Research is supported by Spanish MINECO, Spain under grants PID2019-104141GB-I00 and VA169P20. 3Research is supported by Spanish MCINYU, Spain under grant RTI2018-093521-B-C31 and Spanish State Research Agency, Spain through the national programme Juan de la Cierva-Incorporación 2017. https://doi.org/10.1016/j.cam.2022.114246 0377-0427/©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons. org/licenses/by-nc-nd/4.0/). B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 references therein). One of these techniques is nudging in which a penalty term is added with the aim of driving the approximate solution towards coarse mesh observations of the data. In [6], a new approach, known as continuous data assimilation, is introduced for a large class of dissipative partial differential equations. In this paper we study the numerical approximation of the Navier–Stokes equations with a continuous data assimilation method defined over a reduced order space. The basis functions in the ROM are based only on velocity approximations at different times computed with a mixed finite element Galerkin method using inf–sup stable elements. Both the snapshots and the basis of the ROM satisfy a discrete divergence-free condition. We consider the Navier–Stokes equations (NSE) ∂tu−ν∆u+(u·∇)u+∇p=fin (0,T]×Ω, ∇ ·u=0 in (0,T]×Ω,(1) in a bounded domain Ω⊂Rd,d∈ {2,3}with initial condition u(0) =u0. In (1),uis the velocity field, pthe kinematic pressure, ν > 0 the kinematic viscosity coefficient, and frepresents the accelerations due to external body forces acting on the fluid. The Navier–Stokes equations (1) must be complemented with boundary conditions. For simplicity, we only consider homogeneous Dirichlet boundary conditions u=0on ∂Ω. As in [7] we consider given coarse spatial mesh measurements, corresponding to a solution uof (1), observed at a coarse spatial mesh. We assume that the measurements are continuous in time and error-free and we denote by IH(u) the operator used for interpolating these measurements, where Hdenotes the resolution of the coarse spatial mesh. Since no initial condition for uis available one cannot simulate Eq. (1) directly. To overcome this difficulty it was suggested in [6] to consider instead a solution vof the following system ∂tv−ν∆v+(v·∇)v+∇˜ p=f−β(IH(v)−IH(u)),in (0,T]×Ω, ∇ ·v=0,in (0,T]×Ω,(2) where βis the nudging parameter. In [7] a semidiscrete postprocessed Galerkin spectral method in considered and analyzed. A fully discrete method for the spatial discretization in [7] is analyzed in [8]. In [9] the continuous data assimilation algorithm is analyzed considering both a finite element Galerkin method and a Galerkin method with grad– div stabilization. The extension to the fully discrete case is carried out in [10]. For the Galerkin method with grad–div stabilization the constants in the error bounds in [9,10] are independent on inverse powers of the viscosity parameter. In [11] the authors consider also fully discrete approximations to (2) in which for the spatial discretization the Galerkin method with grad–div stabilization is considered. However, the constants in the error bounds in [11] are not independent on inverse powers of ν. Moreover, in [9,10] there is no need to impose an upper bound on the nudging parameter βas required in [7,8,11]. This fact is important because, on the one hand, there is numerical evidence that no upper bound is required in the numerical experiments and, on the other hand, better results are obtained in some experiments for values of βabove the upper bound assumed in Refs. [7,8,11]. In [12] a continuous data assimilation reduced order model (DA-ROM) method is introduced and analyzed. The idea is to consider a Galerkin approximation to (2) defined in a ROM space. The ROM space is based on a set of snapshots that are fully discrete Galerkin inf–sup stable mixed finite element approximations to (1) at different time steps. The DA-ROM method in [12] is a Galerkin method without any kind of stabilization. The implicit Euler method is used as time integrator and error bounds are proved that converge exponentially fast in time to the true solution. The constants in the error bounds in [12] depend on inverse powers of the viscosity parameter. In the present paper, we follow [12] and consider almost the same DA-ROM with the difference that we add grad–div stabilization. We will call the model grad–div-DA-ROM. We make some improvements compared with the error analysis in [12]. First of all, we prove error bounds in which the constants do not depend on inverse powers of the viscosity. This fact is important in many applications with large Reynolds numbers. A second difference with respect to [12] is the following. In [12] the correlation matrix is based on the inner products of the snapshots without dividing by the number of snapshots as it is standard (see [13]). The reason for not dividing by the number of snapshots is that proceeding in that way one can bound the maximum in time of the L2error between the true solution and the projection onto the ROM space instead of having a bound for a discrete primitive in time of the L2error (let say the mean error, see [13] again). Although an available bound for the maximum norm of the error in the projection simplifies the error analysis, one obtains for the correlation matrix not divided by the number of snapshots that the size of the eigenvalues scales exactly with the number of snapshots. This means that not dividing by the number of snapshots, say Mwhere Mis typically (∆t)−1,∆t being the time step, we get eigenvalues Mtimes larger than using the standard correlation matrix, which in practice implies that the error bounds are multiplied by M(say (∆t)−1). As a consequence, there is no gain using the correlation matrix considered in [12]. In the present paper, we use the standard correlation matrix as defined in [13] and we get error bounds for the error between the grad–div-DA-ROM and the orthogonal L2projection of the true solution onto the ROM space in which we apply the available bound for the mean error instead of requiring a bound for the maximum error. The last improvement respect to [12] is related to the nudging parameter. In the numerical experiments in [12] there is evidence that using a large value for β(say β=100,500) makes a significant difference between the DA-ROM and the 2 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 standard ROM, the first one being much more accurate. Although in [12, Remark 3.8] it is stated that with the analysis presented the usual upper bound on the nudging parameter can be relaxed or even eliminated this is not true. Actually, we found some inconsistencies in the statement of the main Theorem in [12], Theorem 3.5. More precisely, constants α1, α2are defined in the following way α1:= ν−2µ(β2−1)C2 IH2, α2:= 2µ−µC2 I 2β1−µ 2β2−6ν−1C2 b∥Sr∥2∥∇un+1∥2.(3) In (3), the value of µis β, i.e. µis the nudging parameter in (2),His the coarse mesh in (2),nis the time level, CIis a constant related to the interpolant operator IHand Cbis a constant related to a standard bound of the nonlinear term. In [12, Theorem 3.5] it is assumed that αi>0, βi>0, i=1,2. Following the error analysis in [12] we found that the correct value for the constant α2in (3) should be α2:= 2µ−µC2 I β1−2µ β2−6ν−1C2 b∥Sr∥2∥∇un+1∥2, while β2must be larger than 1. Then, in view of the assumption α1>0 we fall essentially into the upper bound ν−2µC2 IH2>0 assumed in Refs. [7,8,11], which means that the upper bound cannot be removed. On the other hand, if we want to relax condition ν−2µC2 IH2>0 we can take β2=1+ϵwith ϵ→0 but in that case in view of the correct value of α2we would need to take β1>(1 +ϵ)C2 I/(2ϵ), which increases as ϵgoes to zero. Since the factor β1µmultiplies the constant in the error bound of Theorem 3.5, relaxing the upper bound in the nudging parameter results in increasing the size of the constants in the error bounds. In the present paper, as in [9,10], we do not need to assume an upper bound on the nudging parameter. For the time integration we use the implicit Euler method although the error analysis for a second order time integrator as BDF2 can be carried out as in [10]. We prove error bounds for the method with constants independent on inverse powers of the viscosity. As in [12] and previous references the error in the initial condition goes to zero exponentially fast. The error in the grad–div-DA-ROM has three components, one coming from the time integrator used, one due to the error in the snapshots (finite element error) and a third one due to the POD method, measured in terms on the eigenvalues of the correlation matrix. Numerical experiments confirm that, for large values of the nudging parameter, the proposed grad– div-DA-ROM rapidly converges to the real solution, and greatly improves the overall accuracy of standard POD schemes up to low viscosities over predictive time intervals, similarly to the DA-ROM in [12]. We want to mention that there are other works in the literature in which reduced order models have been applied in the context of other data assimilation techniques different from the one considered here. See for example [14] for a reduced order approach for 3D variational data assimilation governed by parametrized partial differential equations, and [15,16] for reduced order modeling for four-dimensional variational (4D-Var) data assimilation problems. The outline of the paper is as follows. In Section 2we state some preliminaries and notation. In Section 3we recall the POD method and get some a priori bounds for the orthogonal projection of the true solution onto the POD space. In Section 4we describe the proposed grad–div-DA-ROM and bound the error. Section 5is devoted to show some numerical experiments. Finally, Section 6presents the main conclusions of this work. 2. Preliminaries and notation Let us denote by Q=L2 0(Ω)={q∈L2(Ω)|(q,1) =0}. Let Th=(τh j, φh j)j∈Jh,h>0 be a family of partitions of suitable domains Ωh, where hdenotes the maximum diameter of the elements τh j∈Th, and φh jare the mappings from the reference simplex τ0onto τh j. We shall assume that the partitions are shape-regular and quasi-uniform. Let r≥2, we consider the finite-element spaces Sh,r={χh∈C(Ωh)⏐⏐χh|τh j◦φh j∈Pr−1(τ0)}⊂H1(Ωh), S0 h,r=Sh,r∩H1 0(Ωh), where Pr−1(τ0) denotes the space of polynomials of degree at most r−1 on τ0. We shall denote by (Xh,r,Qh,r−1) the MFE pair known as Hood–Taylor elements [17,18] when r≥3, where Xh,r=(S0 h,r)d,Qh,r−1=Sh,r−1∩L2 0(Ωh),r≥3. To approximate the velocity we consider the discrete divergence-free space Vh,r=Xh,r∩{χh∈H1 0(Ωh)d|(qh,∇ ·χh)=0∀qh∈Qh,r−1}. 3 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 For n≥1 we define the fully discrete Galerkin approximation with the BDF2 time discretization (un h,pn h)∈Xh,r×Qh,r−1 satisfying for all (ϕh, ψh)∈Xh,r×Qh,r−1 (3un h−4un−1 h+un−2 h 2∆t,ϕh)+ν(∇un h,∇ϕh)+bh(un h,un h,ϕh)+(∇pn h,ϕh) =(fn,ϕh), (∇ ·un h, ψh)=0.(4) In (4) un his the Galerkin approximation at time tn,∆tis the time step and bh(·,·,·) is defined in the following way bh(uh,vh,ϕh)=((uh·∇)vh,ϕh)+1 2(∇ ·(uh)vh,ϕh),∀uh,vh,ϕh∈Xh,r. It is straightforward to verify that bhenjoys the skew-symmetry property bh(u,v,w)= −bh(u,w,v)∀u,v,w∈H1 0(Ω)d.(5) Let us fix T>0 and define M=T/∆t. For the fully discrete Galerkin approximation the following bounds hold, see for example [19,20]: ∥un−un h∥0≤C(u,p, ν, r)(hr+(∆t)2),1≤n≤M ∥un−un h∥1≤C(u,p, ν, r)(hr−1+(∆t)2),1≤n≤M.(6) Remark 2.1. The constant C(u,p, ν, r) in (6) depends explicitly on inverse powers of νand on norms of the true solution. In particular, it depends on ∥u∥r,∥p∥r−1. Let us observe that those norms of the solution could also depend on inverse powers of ν. As in [9,11,12], for the error analysis in the present paper we are assuming enough regularity for the solution so that the error bounds (6) hold. For the case in which non local compatibility conditions are not assumed one has to resort to the error bounds in [21]. If we use a stabilized method instead of the Galerkin one we can get bounds with constants independent on inverse powers of ν. For the error analysis we carry out in this paper we need to have velocity approximations with discrete divergence zero. Then, we could start from a Galerkin method with grad–div stabilization as proposed in [22]. A fully discrete version of the Galerkin method with grad–div stabilization and the implicit Euler method is analyzed in [22] resulting in the following bounds: ∥un−un h∥0+h∥un−un h∥1≤C(u,p,r)(hr−1+∆t),1≤n≤M,(7) where the constant C(u,p,r) depends on norms of the true solution but not directly on inverse powers of the viscosity parameter ν. Comparing the error bound (7) with (6) we can observe that instead of rate rin terms of ha rate of convergence r−1 is proved. The numerical experiments in [23] show that this rate is sharp for small values of the viscosity parameter ν. If the family of meshes is quasi-uniform then the following inverse inequality holds for each vh∈Sh,r, see e.g., [24, Theorem 3.2.6], ∥vh∥Wm,p(K)≤cinvhn−m−d(1 q−1 p) K∥vh∥Wn,q(K),(8) where 0 ≤n≤m≤1, 1 ≤q≤p≤ ∞, and hKis the diameter of K∈Th. We consider a modified Stokes projection that was introduced in [25] and that we denote by sm h:V→Vh,rsatisfying (∇sm h,∇ϕh)=(∇u,∇ϕh),∀ϕh∈Vh,r,(9) and the following error bound, see [25]: ∥u−sm h∥0+h∥u−sm h∥1≤C∥u∥jhj,1≤j≤r.(10) From [26], we also have ∥∇sm h∥∞≤C∥∇u∥∞,(11) where Cdoes not depend on νand [9, Lemma 3.8] ∥sm h∥∞≤C(∥u∥d−2∥u∥2)1/2,(12) ∥∇sm h∥L2d/(d−1) ≤C(∥u∥1∥u∥2)1/2,(13) where the constant Cis independent of ν. 4 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Let us denote by PQthe L2orthogonal projection onto Qh,r−1. It holds ∥p−PQp∥0≤Chr−1∥p∥r−1,p∈Q∩Hr−1(Ω).(14) We will also use the well-known property, see [27, Lemma 3.179] ∥∇ ·v∥0≤∥∇v∥0,v∈H1 0(Ω)d.(15) We will assume that the interpolation operator IHis stable in L2, that is, ∥IHu∥0≤c0∥u∥0,∀u∈L2(Ω)d,(16) and that it satisfies the following approximation property, ∥u−IHu∥0≤cIH∥∇u∥0,∀u∈H1 0(Ω)d.(17) The Bernardi–Girault [28], Girault–Lions [29], or the Scott–Zhang [30] interpolation operators satisfy (16) and (17). Notice that the interpolation can be on piecewise constants. 3. Proper orthogonal decomposition We will consider a proper orthogonal decomposition (POD) method. Let us fix T>0 and M>0 and take ∆t=T/M and let us consider the following space V= ⟨u1 h,...,uM h⟩. Let dpbe the dimension of the space V. Let Kbe the correlation matrix corresponding to the snapshots K=((ki,j)) ∈RM×Mwhere ki,j=1 M(ui h,uj h), and (·,·) is the inner product in L2(Ω)d. Following [13] we denote by λ1≥λ2≥ ··· ≥ λdp>0 the positive eigenvalues of Kand by v1,...,vdp∈RMthe associated eigenvectors. Then, the (orthonormal) POD basis is given by ψk=1 √M 1 √λk M ∑ j=1 vj kuh(·,tj),(18) where vj kis the jth component of the eigenvector vkand the following error formula holds, see [13, Proposition 1] 1 M M ∑ j=1uj h− l ∑ k=1 (uj h,ψk)ψk 2 0= dp ∑ k=l+1 λk,(19) where we have used the notation uj h=uh(·,tj). Denoting by Sthe stiffness matrix for the POD basis S=((si,j)) ∈Rdp×dpwith si,j=(∇ψi,∇ψj) then for any v∈Vthe following inverse inequality holds, see [13, Lemma 2] ∥∇v∥0≤√∥S∥2∥v∥0,(20) where ∥S∥2denotes the spectral norm of S. From this inverse inequality we get 1 M M ∑ j=1∇uj h− l ∑ k=1 (uj h,ψk)∇ψk 2 0 ≤∥S∥2 M M ∑ j=1uj h− l ∑ k=1 (uj h,ψk)ψk 2 0≤ ∥S∥2 dp ∑ k=l+1 λk.(21) Instead of (21) we can also apply the following result that is taken from [31, Lemma 3.2] 1 M M ∑ j=1∇uj h− l ∑ k=1 (uj h,ψk)∇ψk 2 0= dp ∑ k=l+1 λk∥∇ψk∥2 0.(22) In the sequel we will denote by Vl= ⟨ψ1,ψ2,...,ψl⟩, and by Plthe L2-orthogonal projection onto Vl. 5 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Although the proof of the following lemma can be found in [31, Lemma 3.3] we include it here for convenience of the readers. Lemma 3.1. Let ube the solution of (1) with initial condition u0and let us denote by uj=u(·,tj), then the following bounds hold, C(u,p, ν, r)being the constant in (6) 1 M M ∑ j=1∥uj−Pluj∥2 0≤C0,P:= 2C(u,p, ν, r)(h2r+(∆t)4)+2 dp ∑ k=l+1 λk, 1 M M ∑ j=1∥∇(uj−Pluj)∥2 0≤C1,P:= 3C(u,p, ν, r)(h2(r−1) +(∆t)4) (23) +3 dp ∑ k=l+1 λk∥∇ψk∥2 0+3C(u,p, ν, r)∥S∥2(h2r+(∆t)4). Proof. By definition of the Plprojection ∥uj−Pluj∥0≤ ∥uj−Pluj h∥0. Then ∥uj−Pluj∥2 0≤2∥uj−uj h∥2 0+2∥uj h−Pluj h∥2 0. Applying now (6) and (19) we prove the first inequality in (23). To prove the second one we write ∥∇(uj−Pluj)∥2 0≤3∥∇(uj−uj h)∥2 0+3∥∇(uj h−Pluj h)∥2 0 +3∥∇(Pluj h−Pluj)∥2 0. Taking into account that applying (20) we get ∥∇(Pluj h−Pluj)∥2 0≤ ∥S∥2∥Pl(uj h−uj)∥2 0≤ ∥S∥2∥uj h−uj∥2 0 we conclude by applying (6) and (22).□ 3.1. A priori bounds for the orthogonal projection onto Vl In this section we will prove some a priori bounds for the orthogonal projection Pluj,j=0,...,M, that are needed in the error analysis of the rest of the paper. To this end, in a first step, we get some a priori bounds for the Galerkin velocity approximation. Lemma 3.2. Let uj hbe the Galerkin velocity approximation at time tjdefined in (4). Let us assume ∆t≤Chd/4,(24) for any positive constant C. Then, the following bounds hold, C(u,p, ν, r), r =2,3, being the constant in (6) ∥uj h∥∞≤Cu,inf := C(C(u,p, ν, 2) +(∥uj∥d−2∥uj∥2)1/2),(25) ∥∇uj h∥∞≤Cu,1,inf := C(C(u,p, ν, 3) +∥∇u∥L∞(L∞)),(26) ∥∇uj h∥L2d/(d−1) ≤Cu,ld := C(C(u,p, ν, 2) +(∥u∥1∥u∥2)1/2).(27) Proof. To prove (25) we observe that using (8),(12),(6) and (10) we get ∥uj h∥∞≤ ∥uj h−sm h(·,tj)∥∞+∥sm h(·,tj)∥∞ ≤Ch−d/2∥uj h−sm h(·,tj)∥0+C(∥uj∥d−2∥uj∥2)1/2 ≤Ch−d/2C(u,p, ν, 2)(h2+∆t2)+C(∥uj∥d−2∥uj∥2)1/2≤Cu,inf, whenever we assume condition (24) holds. In the error bound (25) we have included the factor Ch2∥uj∥2coming from the error ∥uj−sm h(·,tj)∥0into the factor C(u,p, ν, 2)h2coming from the error of the Galerkin method since C(u,p, ν, 2) depends on ∥u∥L∞(H2). We refer to Remark 2.1 for the assumed regularity of the solution. 6 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 To prove (26) we apply (8),(11),(6) and (10) we get ∥∇uj h∥∞≤ ∥∇uj h−∇sm h(·,tj)∥∞+∥∇sm h(·,tj)∥∞ ≤Ch−d/2∥uj h−sm h(·,tj)∥1+C∥∇uj∥∞ ≤Ch−d/2C(u,p, ν, 3) (h2+(∆t)2)+C∥∇uj∥∞≤Cu,1,inf, whenever condition (24) holds. Finally, for the bound (27) we use (8),(13),(6) and (10) and assume again condition (24) holds (indeed the weaker condition ∆t≤Ch1/4would be enough) to get ∥∇uj h∥L2d/(d−1) ≤ ∥∇(uj h−sm h(·,tj))∥L2d/(d−1) +∥∇sm h(·,tj)∥L2d/(d−1) ≤Ch−1/2∥uj h−sm h(·,tj)∥1+C(∥u∥1∥u∥2)1/2 ≤Ch−1/2C(u,p, ν, 2)(h+(∆t)2)+C(∥u∥1∥u∥2)1/2≤Cu,ld.□ Now, we prove a priori bounds in the same norms for Pluj. Lemma 3.3. Under the same conditions of Lemma 3.2 the following bounds hold for the L2projection Pluj. ∥Pluj∥∞≤Cinf := Cu,inf +C+h−d/2√M√λl+1,(28) ∥∇Pluj∥∞≤C1,inf := Cu,1,inf +C+h−d/2√M∥S∥1/2 2⎛ ⎝ dp ∑ k=l+1 λk⎞ ⎠ 1/2 ,(29) ∥∇Pluj∥L2d/(d−1) ≤Cld := Cu,ld +C+h−1/2√M∥S∥1/2 2⎛ ⎝ dp ∑ k=l+1 λk⎞ ⎠ 1/2 .(30) Proof. Using inverse inequality (8),(25), the stability of the Plprojection, and (6) we get ∥Pluj∥∞≤ ∥uj h∥∞+∥Pluj−uj h∥∞≤Cu,inf +h−d/2∥Pluj−uj h∥0 ≤Cu,inf +h−d/2∥Pl(uj−uj h)∥0+h−d/2∥Pluj h−uj h∥0 ≤Cu,inf +h−d/2∥uj−uj h∥0+h−d/2∥Pluj h−uj h∥0 ≤Cu,inf +h−d/2C(u,p, ν, 2)(h2+(∆t)2)+h−d/2∥Pluj h−uj h∥0 ≤Cu,inf +C+h−d/2∥Pluj h−uj h∥0,(31) where in the last inequality we assume, as before, condition (24). In view of (19) we can write for the last term ∥Pluj h−uj h∥0≤M1/2(∑dp k=l+1λk)1/2 . Actually, this estimate can be slightly improved with the following argument. It is easy to see that uj h−Pluj h= dp ∑ k=l+1 (uj h,ψk)ψk. Using the definition of ψkit is also easy to observe that uj h−Pluj h= dp ∑ k=l+1 (uj h,ψk)ψk=√M dp ∑ k=l+1√λkvj kψk. And then ∥uj h−Pluj h∥0=√M⎛ ⎝ dp ∑ k=l+1 λk|vj k|2⎞ ⎠ 1/2 ≤√M√λl+1⎛ ⎝ dp ∑ k=l+1|vj k|2⎞ ⎠ 1/2 ≤√M√λl+1,(32) where in the last inequality we have used that (∑dp k=l+1|vj k|2)1/2 ≤1 since the matrix with columns the vectors vkcan be enlarged to an M×Morthogonal matrix. Inserting (32) into (31) we finally prove (28). 7 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Arguing similarly, using inverse inequality (8),(26),(20), the stability of the Plprojection, and (6) we get ∥∇Pluj∥∞≤ ∥∇uj h∥∞+∥∇(Pluj−uj h)∥∞ ≤Cu,1,inf +h−d/2∥∇(Pluj−uj h)∥0 ≤Cu,1,inf +h−d/2∥∇(Pl(uj−uj h))∥0+h−d/2∥∇(Pluj h−uj h)∥0 ≤Cu,inf +h−d/2∥S∥1/2 2∥uj−uj h∥0+h−d/2∥∇(Pluj h−uj h)∥0 ≤Cu,inf +h−d/2∥S∥1/2 2C(u,p, ν, 2)(h2+(∆t)2) +h−d/2∥∇(Pluj h−uj h)∥0 ≤Cu,inf +C+h−d/2∥∇(Pluj h−uj h)∥0.(33) Finally, from (22) we get ∥∇(Pluj h−uj h)∥0≤√M∥S∥1/2 2⎛ ⎝ dp ∑ k=l+1 λk⎞ ⎠ 1/2 , which inserted into (33) gives (29). Arguing exactly as before and applying (27) we also obtain (30).□ Remark 3.4. Let us observe that for the error analysis we assume that the constants Cinf,C1,inf and Cld in (28),(29) and (30), respectively, are bounded, which can always be obtained for llarge enough in the POD approximation (34). For the numerical experiments, however, we observed that the method worked equally well for all the different numbers of modes, l, we considered. Then, the restriction on the value of lcoming from the need of assuming bounded constants Cinf, C1,inf and Cld, seems not to be a problem for the proposed method in practice. 4. The POD data assimilation algorithm For any initial condition the POD data assimilation approximation using the implicit Euler method and grad–div stabilization is obtained by solving for n≥1: (un l−un−1 l ∆t,ϕl)+ν(∇un l,∇ϕl)+bh(un l,un l,ϕl)+µ(∇ ·un l,∇ ·ϕl) =(fn,ϕl)−β(IHun l−IHun,IHϕl),∀ϕl∈Vl,(34) where µis the grad–div stabilization parameter, βis the nudging parameter and IHis an interpolation operator over a coarse mesh. Theorem 4.1. Let un lbe the grad–div-DA-ROM approximation defined in (34), let unbe the velocity approximation of the Navier–Stokes equations (1) at time tnand let Plunbe its orthogonal projection over the POD space Vl. Assuming the solution (u,p)of (1) is smooth enough the following bound holds for a constant C independent on νand β ∥un l−Plun∥2 0≤1 (1+γ 2∆t)n∥e0 l∥2 0+TC1,P(ν+2µ+2 L(∥u∥2+Cld +Cinf)) +Tβc2 0C0,P+C µh2(r−1)∆t n ∑ j=1∥pj∥2 r−1+C(∆t)2 L∫tn 0∥utt (s)∥2 0ds,(35) where C0,P, C1,Pare the constants in (23), and Cld, Cinf are the constants in (28),(30). Proof. Following [12] we will compare un lwith Plun. It is easy to obtain (Plun−Plun−1 ∆t,ϕl)+ν(∇Plun,∇ϕl)+bh(Plun,Plun,ϕl) +µ(∇ ·Plun,∇ ·ϕl)=(fn,ϕl)+ν(∇τn 1,∇ϕl)+(τn 2,∇ ·ϕl) +(τn 3,ϕl)+(τn 4,ϕl),∀ϕl∈Vl,(36) where τn 1,τn 2,τn 3and τn 4are defined by: τn 1=(Plun−un), τn 2=(pn−PQ(pn))+µ(∇ ·(Plun−un)), 8 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 τn 3=1 ∆t(un−un−1)−un t,(37) (τn 4,ϕl)=bh(Plun,Plun,ϕl)−bh(un,un,ϕl), and we denote by PQthe L2orthogonal projection onto Qh,r−1. Let us denote by en l=un l−Plun. Subtracting (36) from (34) and taking ϕl=en lwe get 1 2∆t(∥en l∥2 0−∥en−1 l∥2 0)+ν∥∇en l∥2 0+µ∥∇ ·en l∥2 0+β∥IHen l∥2 0(38) ≤ −bh(un l,un l,en l)+bh(Plun,Plun,en l)+β(IH(un−Plun),IHen l) −ν(∇τn 1,∇en l)−(τn 2,∇ ·en l)−(τn 3,en l)−(τn 4,en l). We will argue as in [10]. For the first term on the right-hand side of (38) using the skew-symmetric property (5) we get ⏐⏐bh(un l,un l,en l)−bh(Plun,Plun,en l)⏐⏐=⏐⏐bh(en l,Plun,en l)⏐⏐ ≤ ∥∇Plun∥∞∥en l∥2 0+1 2∥∇ ·en l∥0∥Plun∥∞∥en l∥0 ≤L 2∥en l∥2 0+µ 4∥∇ ·en l∥2 0,(39) where L=2 max n≥0(∥∇Plun∥∞+1 4µ∥Plum∥2 ∞)≤2(C1,inf +C2 inf 4µ),(40) and we have applied (28) and (29) in the last inequality. As pointed out in Remark 3.4 Lis a bounded constant. For the second term on the right-hand side of (38), applying the L2-stability of the interpolation operator (16) we get β(IH(un−Plun),IHen l)≤βc0∥un−Plun∥0∥IHen l∥0 ≤β 2c2 0∥un−Plun∥2 0+β 2∥IHen l∥2 0.(41) For the truncation errors we write |ν(∇τn 1,∇en l)| ≤ ν 2∥∇τn 1∥2 0+ν 2∥∇en l∥2 0, |(τn 2,∇ ·en l)| ≤ ∥τn 2∥2 0 µ+µ 4∥∇ ·en l∥2 0,(42) |(τn 3+τn 4,en l)| ≤ 1 2L∥τn 3+τn 4∥2 0+L 2∥en l∥2 0. Inserting (39),(41) and (42) into (38) we get 1 2 1 ∆t(∥en l∥2 0−∥en−1 l∥2 0)+ν 2∥∇en l∥2 0+β 2∥IHen l∥2 0+µ 2∥∇ ·en l∥2 0(43) ≤L∥en l∥2 0+ν 2∥∇τn 1∥2 0+∥τn 2∥2 0 µ+1 2L∥τn 3+τn 4∥2 0+β 2c2 0∥un−Plun∥2 0. The following argument is taken from [9,10]. We first observe that L∥en l∥2 0≤2L∥IHen l∥2 0+2L∥(I−IH)en l∥2 0 so that assuming β≥8L(44) and multiplying (43) by 2 we obtain 1 ∆t(∥en l∥2 0−∥en−1 l∥2 0)+ν∥∇en l∥2 0+β 2∥IHen l∥2 0+µ∥∇ ·en l∥2 0−4L∥(I−IH)en l∥2 0 ≤ν∥∇τn 1∥2 0+2∥τn 2∥2 0 µ+1 L∥τn 3+τn 4∥2 0+βc2 0∥un−Plun∥2 0. Applying (17) we have ν∥∇en l∥2 0−4L∥(I−IH)en l∥2 0≥ν∥∇en l∥2 0−4Lc2 IH2∥∇en l∥2 0≥ν 2∥∇en l∥2 0, 9 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 6. Example 5.1 (Case Re =100): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes (166 snapshots used, which comprise one full period from t=5 s to t=5.332 s). but this falls outside the scope of the present work. Recently, the effect of adding the grad–div term as itself in the POD setting without the DA term has been numerically explored and extensively studied in [44] for the first time. Indeed, although the grad–div stabilization term has been already considered e.g. in [45,46] within a ROM framework, actually in [45] it has been embedded within a residual-based VMS [47,48] method, thus making difficult to understand its real contribution, while in [46] it has been neglected in the numerical studies. This term generally provides improvement of local discrete mass conservation [49,50], and thus it is particularly important in the present framework, in which mixed interpolations that satisfy the inf–sup condition but are not exactly divergence-free have been used to compute the snapshots. This allows to work with only velocity ROM, as in this case, since the POD velocity modes are solenoidal and the pressure term drops out, but could lead to a poor resolution, as the G-ROM results confirm. In [44], where we also discussed different techniques of pressure recovery for POD-ROM, we showed the benefits of adding the grad–div term to the G-ROM to improve in particular the energy prediction, thus we found convenient to add it to the G-ROM in the present numerical experiments too. Following Remark 4.3 and [44], to select the grad–div parameter µfor the following 16 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 7. Example 5.1 (Case Re =100): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes for DA-ROM with β=10,100,500 (166 snapshots used, which comprise one full period from t=5 s to t=5.332 s). numerical experiments, we have considered a constant value fixed minimizing the L∞error in time with respect to the snapshots energy computed in one period and then repeated in the rest of periods, thus being the snapshots data to construct the reduced basis sufficient to compute the constant µ, and no further information is needed. This allows in particular to obtain a reliable energy prediction. In [44], an adaptive in time algorithm for the grad–div parameter µhas been also proposed and numerically investigated. The adaptive in time strategy consists in adjusting µaround a constant value chosen as above, so that the contribution of the grad–div stabilization term removes dissipation if the ROM energy is too small, and adds dissipation if the energy is too large with respect to the FEM energy. This additionally allows to guarantee a very long time accuracy. For the time intervals considered in the present paper, we found that the constant µalready provides a reasonable accuracy. We emphasize, however, that when considering DA into the ROM, thus adding or not the grad–div stabilization term makes no significant difference and a reliable energy prediction for relatively long time integrations is similarly approached using just the DA term, as showed in the following numerical experiments. Nevertheless, introducing the grad–div term in the DA-POD setting improved the numerical analysis upon the existing 17 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 8. Example 5.1 (Case Re =100): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes for grad–div-DA-ROM with µ=0.15 and β=10,100,500 (166 snapshots used, which comprise one full period from t=5 s to t=5.332 s). DA-POD methods. Indeed, as showed in [44], obtaining error bounds with constants independent on inverse powers of the viscosity parameter can help to find good a priori error indicators that, at least for few modes (of interest in practice), almost match the computed errors over predictive time intervals. 5.1. Case Re =100 In this section, we discuss results for Re =100. In this case, we have used the computational grid represented in Fig. 1 on top to compute the snapshots, for which h=2.76 ×10−2m, resulting in 32 488 d.o.f. for velocities and 4 151 d.o.f. for pressure. Also, 166 snapshots were collected, which comprise one full period from t=5 s to t=5.332 s. All tested ROM have been run in the stable response time interval [5,7]s, corresponding to six periods for the lift coefficient. Thus, we are actually testing the ability of the considered ROM to predict/extrapolate in time, monitoring their performance over a six times larger time interval with respect to the one used to compute the snapshots and generate the POD modes. This 18 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 9. Example 5.1 (Case Re =100): First POD velocity modes (Euclidean norm) obtained with 166 snapshots (full period basis, left) and 106 snapshots (inaccurate basis corresponding to 64% of one full period, right). Table 2 Example 5.1 (Case Re =100): Errors levels with respect to DNS for G-ROM, grad–div-ROM (µ=0.15), DA-ROM (β=500), and grad–div-DA-ROM (µ=0.15, β =500) (166 snapshots used, which comprise one full period from t=5 s to t=5.332 s). Errors Re =100 G-ROM grad–div-ROM DA-ROM grad–div-DA-ROM Emax kin 4.56e−02 1.60e−05 8.20e−05 4.30e−05 cmax D3.84e−01 7.15e−02 2.72e−03 2.87e−03 cmax L6.78e−01 4.33e−02 4.18e−03 4.76e−03 ℓ2(L2)unorm 1.68e−01 9.73e−02 2.04e−02 2.04e−02 will show how the strategy to incorporate DA into the ROM can already provide long time stability and accuracy, thus proving its robustness. Numerical results for energy, drag and lift predictions using l=8 modes are shown in Figs. 6,7,8. In particular, Fig. 6 shows a comparison within DNS, G-ROM, grad–div-ROM with µ=0.15 (fixed as described above), DA-ROM with β=10, and grad–div-DA-ROM with µ=0.15 and β=10. From this figure, we observe that, whereas the G-ROM solution is totally inaccurate, the application of the grad–div stabilization term greatly improves the G-ROM solution, allowing to compute rather accurate quantities of interest. Indeed, the temporal evolution of the kinetic energy and lift coefficient is very close to that of the DNS, being the drag coefficient temporal evolution the most sensitive quantity presenting larger differences. A slight improvement is observed for using DA with β=10, being results for DA-ROM and grad–div-DA-ROM almost identical. Note that using DA, since we started from zero initial velocity conditions, the DNS results are approached around t=5.4 s with β=10. A significant improvement is observed by increasing the nudging parameter βfor DA reduced order methods. This is clearly displayed in Figs. 7,8, which respectively show the behavior of the DA-ROM and the grad–div-DA-ROM, varying the nudging parameter βfrom 10 to 500. Again, almost identical results are obtained with both DA reduced order methods, for which the best predictions are given by the largest values β=500 of the nudging parameter, although we observe a similar accuracy already for β=100. Note also that for large values of the nudging parameter (β=100,500), although we started from zero initial velocity conditions, the DNS results are approached with a rather accurate resolution just after very few iterations (around 20, i.e. 0.04 s, for β=100 and 5, i.e. 0.01 s, for β=500). All these results are also confirmed by Table 2, where we display the error levels with respect to DNS of maximum kinetic energy |Emax kin,l−Emax kin,DNS|, maximum drag coefficient |cmax D,l−cmax D,DNS|, maximum lift coefficient |cmax L,l−cmax L,DNS|, and velocity norm ∥ul−uDNS ∥ℓ2(L2)using l=8 modes for G-ROM, grad–div-ROM (µ=0.15), DA-ROM (β=500), and grad–div-DA-ROM (µ=0.15, β =500) in the time interval [5.01,7]s. Note how grad–div-ROM already reduces the error level in Emax kin of three orders of magnitude with respect to G-ROM, similarly to both DA reduced order methods, and in cmax Lof one order of magnitude, while both 19 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 10. Example 5.1 (Case Re =100): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes (106 snapshots used, which comprise 64% of one full period from t=5 s to t=5.212 s). DA reduced order methods of two orders of magnitude. However, for cmax D, while grad–div-ROM slightly reduces the error level with respect to G-ROM (five times), both DA reduced order methods guarantee again a reduction of two orders of magnitude. In terms of ℓ2(L2) velocity norm, both DA reduced order methods reduces the G-ROM error level eight times, while the grad–div-ROM is just slightly better accurate than G-ROM. We also investigate the considered ROM performances in predicting quantities of interest when inaccurate snapshots (64% of one full period) are used in their construction. Thus, we generate inaccurate snapshots using 64% of one full period of DNS data, which corresponds in this case to the first 106 DNS time step solutions from t=5 s to t=5.212 s. Fig. 9 displays the Euclidean norm of the first POD velocity modes obtained with the full set of snapshots (left) and the inaccurate set of snapshots (right). Results for the considered ROM using l=8 modes in this case are shown in Figs. 10, 11,12. Similar to the previous results, DA significantly improves the accuracy of the G-ROM, especially for large values of the nudging parameter, without the need to increase the number of reduced basis functions. While results for G-ROM becomes more and more inaccurate as time goes on, results for grad–div-ROM remain still acceptable if compared with 20 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 11. Example 5.1 (Case Re =100): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes for DA-ROM with β=10,100,500 (106 snapshots used, which comprise 64% of one full period from t=5 s to t=5.212 s). DA reduced order methods for a small value of the nudging parameter. Again, results for both DA-ROM (with and without grad–div term) are very close and almost approaches DNS results for large values of the nudging parameter. Actually, they are almost comparable to previous results for one full period of DNS data. All these considerations are also reflected by the error levels displayed in Table 3. These results suggest that, despite its simple implementation, DA can greatly improve the overall accuracy of the standard G-ROM in the computation of quantities of interest even when low-resolution data are available to construct the reduced basis, which is common in practice, whereas grad–div stabilization (without DA) continues providing reliable results. We notice, however, that as the Reynolds number is increased (see next section), results for grad–div-ROM (without DA) are less accurate, and maybe it should be combined with convection stabilization if one does not use DA in order to obtain more accurate results. 21 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 12. Example 5.1 (Case Re =100): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes for grad–div-DA-ROM with µ=0.15 and β=10,100,500 (106 snapshots used, which comprise 64% of one full period from t=5 s to t=5.212 s). 5.2. Case Re =1000 In this section, we discuss results for Re =1000. In this case, we have used a finer computational grid with respect to Re =100 to compute the snapshots (see Fig. 2 on top, for which h=1.46 ×10−2m, resulting in 101 820 d.o.f. for velocities and 12 885 d.o.f. for pressure). This has been necessary to obtain stable DNS results. However, the coarse mesh for DA in ROM is the same as for the previous case (see Fig. 2 on bottom). The full period length of the statistically steady state is now 0.22 s, so that 110 snapshots were collected, starting from t=5 s. Again, all tested ROM have been run in the stable response time interval [5,7]s, corresponding now to nine periods for the lift coefficient. This time range is thus nine times wider with respect to the time window used for the generation of the POD modes, so that at the higher Reynolds number we are performing the longer time integration with respect to the time interval used to compute the snapshots. 22 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 13. Example 5.2 (Case Re =1000): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes (110 snapshots used, which comprise one full period from t=5 s to t=5.22 s). Numerical results for energy, drag and lift predictions using l=8 modes are shown in Figs. 13,14,15. In particular, Fig. 13 shows a comparison within DNS, G-ROM, grad–div-ROM with µ=0.001 (fixed as described above), DA-ROM with β=10, and grad–div-DA-ROM with µ=0.001 and β=10. As already noticed in the previous case, from this figure we observe that, whereas the G-ROM solution is totally inaccurate, the application of the grad–div stabilization term helps to improve the G-ROM solution, although it shows larger error levels than the lower Reynolds number case Re =100 when compared to DNS results. A slight improvement is observed again for using DA with β=10, being results for DA-ROM and grad–div-DA-ROM almost identical. Looking at the temporal evolution of the kinetic energy (on top), we observe that also in this case the DA results almost stabilize around t=5.4 s with β=10, even if the reached values under-estimate the DNS results. Increasing the nudging parameter βfrom 10 to 500 for DA reduced order methods (see Figs. 14,15) already allows to almost approach DNS results, although we note a detachment in predicting cD,cLas time increases. Almost identical results are obtained with both DA reduced order methods, for which the best predictions are given by the largest values 23 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 14. Example 5.2 (Case Re =1000): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes for DA-ROM with β=10,100,500 (110 snapshots used, which comprise one full period from t=5 s to t=5.22 s). β=500 of the nudging parameter, although we observe a similar accuracy already for β=100. Note again that for large values of the nudging parameter (β=100,500), the DNS results are almost approached just after very few iterations (around 20, i.e. 0.04 s, for β=100 and 5, i.e. 0.01 s, for β=500). All these results are confirmed by Table 4. Note that grad–div-ROM now just slightly reduces the error levels with respect to G-ROM for all quantities, while both DA reduced order methods still guarantee a reduction of two orders of magnitude for Emax kin ,cmax D, and five times for cmax L. In terms of ℓ2(L2) velocity norm, both DA reduced order methods reduces the G-ROM error level by a factor of 6.5, while the grad–div-ROM is just slightly better accurate than G-ROM. Also in this case we finally investigate the considered ROM performances in predicting quantities of interest when inaccurate snapshots (64% of one full period) are used in their construction. Thus, we generate inaccurate snapshots using 64% of one full period of DNS data, which corresponds in this case to the first 70 DNS time step solutions from 24 B. García-Archilla, J. Novo and S. Rubino Journal of Computational and Applied Mathematics 411 (2022) 114246 Fig. 15. Example 5.2 (Case Re =1000): Temporal evolution of kinetic energy, drag coefficient and lift coefficient using l=8 modes for grad–div-DA-ROM with µ=0.001 and β=10,100,500 (110 snapshots used, which comprise one full period from t=5 s to t=5.22 s). Table 3 Example 5.1 (Case Re =100): Errors levels with respect to DNS for G-ROM, grad–div-ROM (µ=0.15), DA-ROM (β=500), and grad–div-DA-ROM (µ=0.15, β =500) (106 snapshots used, which comprise 64% of one full period from t=5 s to t=5.212 s). Errors Re =100 (Inaccurate snapshots) G-ROM grad–div-ROM DA-ROM grad–div-DA-ROM Emax kin 4.34e−02 4.13e−04 2.30e−05 6.10e−05 cmax D3.75e−01 6.40e−02 1.14e−02 1.16e−02 cmax L6.29e−01 1.89e−02 1.16e−02 1.11e−02 ℓ2(L2)unorm 1.76e−01 1.01e−01 2.99e−02 2.99e−02 25