Uniform in Time Error Estimates for a Finite Element Method Applied to a Downscaling Data Assimilation Algorithm for the Navier--Stokes Equations

UNIFORM IN TIME ERROR ESTIMATES FOR A FINITE ELEMENT METHOD APPLIED TO A DOWNSCALING DATA ASSIMILATION ALGORITHM FOR THE NAVIER STOKES EQUATIONS BOSCO GARC´ IA-ARCHILLA∗, JULIA NOVO†,AND EDRISS S. TITI‡ Abstract. In this paper we analyze a finite element method applied to a continuous downscaling data assimilation algorithm for the numerical approximation of the two and three dimensional Navier- Stokes equations corresponding to given measurements on a coarse spatial scale. For representing the coarse mesh measurements we consider different types of interpolation operators including a Lagrange interpolant. We obtain uniform-in-time estimates for the error between a finite element approximation and the reference solution corresponding to the coarse mesh measurements. We consider both the case of a plain Galerkin method and a Galerkin method with grad-div stabilization. For the stabilized method we prove error bounds in which the constants do not depend on inverse powers of the viscosity. Some numerical experiments illustrate the theoretical results. Key words. data assimilation, downscaling, Navier-Stokes equations, uniform-in-time error estimates, mixed finite elements AMS subject classifications. 35Q30, 65M12, 65M15, 65M20, 65M60, 65M70, 76B75 1. Introduction. Data assimilation refers to a class of techniques that combine experimental data and simulation in order to obtain better predictions in a physical system. There is a vast literature on data assimilation methods, specially in the recent years (see e.g., [4], [33], [35], [41], and the references therein). One of these techniques is nudging, where a penalty term is added in order to drive the approximate solution towards coarse mesh or large scale spatial observations of the data. In a recent work [7], a new approach, known as continuous data assimilation, is introduced for a large class of dissipative partial differential equations, including Rayleigh-B´enard convection [19], the planetary geostrophic ocean dynamics model [20], etc. (see also references therein). Continuous data assimilation has also been used in numerical studies, for example, with the Chafee-Infante reaction-diffusion equation the Kuramoto-Sivashinsky equation (in the context of feedback control) [36], Rayleigh-B´enard convection equations [3], [18], and the Navier-Stokes equations [25], [28]. However, there is much less numerical analysis of this technique. The present work concerns with the numerical analysis of continuous data assimilation for the Navier-Stokes equations when discretized with mixed finite element methods (MFE). To be more precise, we consider the Navier-Stokes equations (NSE) ∂tu−ν∆u+ (u·∇)u+∇p=fin (0, T]×Ω, ∇·u= 0 in (0, T]×Ω,(1.1) ∗Departamento de Matem´atica Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research is supported by Spanish MINECO under grant MTM2015-65608-P ([email protected]). †Departamento de Matem´aticas, Universidad Aut´onoma de Madrid, Spain. Research is supported by Spanish MINECO under grant MTM2016-78995-P (AEI/FEDER, UE) and VA024P17 (Junta de Castilla y Leon, ES) cofinanced by FEDER funds (julia.nov[email protected]). ‡Department of Mathematics, Texas A&M University, College Station, TX 77843, USA. Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK. Department of Computer Science and Applied Mathematics, Weizmann Institute of Science, Rehovot 76100, Israel. Research is supported in part by the ONR grant N00014-15-1-2333, the Einstein Stiftung/Foundation - Berlin, through the Einstein Visiting Fellow Program, and by the John Simon Guggenheim Memorial Foundation ([email protected] [email protected]). 1 Esta es la versión de autor del artículo publicado en = This is an author produced version of a paper published in: SIAM Journal on Numerical Analysis, 58.1 (2020): 410-429 DOI: https://doi.org/10.1137/19M1246845 Copyright: © 2020 by SIAM. Unauthorized reproduction of this article is prohibited El acceso a la versión del editor puede requerir la suscripción del recurso = Access to the published version may require subscription 2B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI in a bounded domain Ω ⊂Rd,d∈ {2,3}. In (1.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.1) must be complemented with boundary conditions. For simplicity, we only consider homogeneous Dirichlet boundary conditions u=0on ∂Ω. Following [37] we consider given coarse spatial scale measurements, corresponding to a solution uof (1.1), observed at a coarse spatial mesh. The measurements are assumed to be continuous in time and error-free. We denote by IH(u) the operator used for interpolating these measurements, where Hdenotes the resolution of the coarse spatial mesh. Since the initial condition for uis missing one cannot compute uby simulating equation (1.1) directly. To overcome this difficulty it was suggested in [7] to consider instead a solution vof the following approximating system ∂tv−ν∆v+ (v·∇)v+∇˜p=f−β(IH(v)−IH(u)),in (0, T]×Ω, ∇·v= 0,in (0, T]×Ω,(1.2) where βis the relaxation (nudging) parameter. In the case of the Navier-Stokes equations (and indeed, of many other nonlinear dissipative systems), it is well-known that for relatively not so small Reynolds numbers, solutions are unstable and even chaotic. For this reason, it is expected that any small error in the initial data could lead to exponentially growing error in the solutions. Notably, the instabilities in the NSE occur at the large spatial scales, while the fine scales are stabilized by the viscosity. For this reason once the large spatial scales are stabilized, as it is done in the proposed downscaling data assimilation approximation, equation (1.2), the corresponding solution are stable and converge to the same solution uthat is corresponding to IH(u). This is the very reason that small errors are not magnified in time and allows to obtain uniform in time error bounds. In this paper we consider a semidiscretization in space with inf-sup stable mixed finite elements for equation (1.2) and analyze two different methods: the Galerkin method and the Galerkin method and grad-div stabilization. Grad-div stabilization was originally proposed in [21] to improve the conservation of mass in finite element methods (see also [38], [39]) However, it has been observed in the simulation of turbulent flows, [32], [42], that using only grad-div stabilization produced stable (nonoscillating) simulations. We prove uniform-in-time error estimates for approximating the unknown reference solution, u, that corresponds to the coarse spatial scale measurement IH(u). For the Galerkin method without grad-div stabilization, the spatial error bounds we prove are optimal, in the sense that the rate of convergence is that of the best interpolant. In the case we add grad-div stabilization, as in [15], [16], we get error bounds in which the error constants do not depend on inverse powers of the viscosity parameter ν. This fact is of importance in many applications where viscosity is orders of magnitude smaller than the velocity (i.e., large Rynolds number). The convergence rates we prove in our error bounds are sharp and confirmed by numerical experiments. We now comment on the analysis of numerical methods for (1.2). In [37], a semidiscrete postprocessed Galerkin spectral method for the two-dimensional Navier- Stokes equations is studied. Under suitable conditions on the nudging parameter β and the coarse mesh resolution H, uniform-in-time error estimates are obtained for the difference between the numerical approximation to vand u. Furthermore, the use of a postprocessing technique introduced in [23] [24], allows for higher convergence rates than a standard spectral Galerkin method. A fully-discrete method for the spatial ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 3 discretization in [37] is analyzed in [31], where the backward Euler method is used for time discretization. Fully implicit and semi-implicit methods are considered, and optimal uniform-in-time error estimates are obtained with the same convergence rate in space as in [37]. More closely related to the present work are [34] and [40]. In [40] they only analyze linear problems and, for the proof of the results on the Navier-Stokes equations they present, they refer to [34] with some differences that they point out. They also present a wide collection of numerical experiments. In [34], the authors consider fully discrete approximations to equation (1.2) where the spatial discretization is performed with a MFE Galerkin method plus grad-div stabilization. A second order IMEX in time scheme is analyzed in [34], and, as in [31], [37] and the present paper, uniform-in- time error bounds are obtained. Compared with [34], for the same convergence rate, the error bounds in the present paper have constants that do not depend on inverse powers of the viscosity parameter ν(Theorem 3.3) or, for similar error constants, error bounds in the present paper have an order of convergence one unit larger (Theorem 3.2 below). Also, the analysis in [34] is restricted to IHubeing an interpolant for non smooth functions (Cl´ement, Scott-Zhang, etc), since it makes explicit use of bound (2.19), which is not valid for nodal (Lagrange) interpolation (neither it is (2.20)). In the present paper, we prove error bounds for the case in which (2.19) holds, but also for the case in which IHuis a standard Lagrange interpolant (Theorem 3.12 below). To our knowledge, this is the first time in the literature where such kind of bounds are proved. Also, compared with [34] and [40], we remove the upper bound assumed on the nudging parameter β. The authors of [34] had observed (see [34, Remark 3.8]) that the upper bound they required in the analysis does not hold in the numerical experiments and they state that a different approach to the analysis should be used to remove the upper bound on β. An analogous upper bound on βappears also in [31] and [37], where the value of Hdepends on the inverse of the nudging parameter βwhich means that increasing the value of βwould require a smaller value of H. Although the analysis of the present paper could be extended to fully discrete methods following for example the techniques in [15], [16] we believe that the new ideas introduced in the present paper are easier to understand in the framework of the semidiscrete methods. The extension of the analysis of the present paper to the fully discrete case will be subject of future work. The rest of the paper is as follows. Section 2is devoted to preliminary material, in Section 3we introduce and analyze the finite element method for equation (1.2) with and without grad-div stabilization. In Subsection 3.1 we analyze the case in which IHuis the standard Lagrange interpolant. Finally, in Section 4some numerical experiments are shown to illustrate the theoretical results. 2. Preliminaries and Notation. Throughout the paper, Ws,p(D) will denote the Sobolev space of real-valued functions defined on the domain D⊂Rdwith distributional derivatives of order up to sin Lp(D). We denote by |·|s,p,D standard seminorm, and, following [14], for Ws,p(D) we will use the norm k·ks,p,D defined by kfkp s,p,D = s X j=0 |D|p(j−s) d|f|p j,p,D , where |D|stands for the Lebesgue measure of Dso that kfkm,p,D |D|m d−1 pis scale invariant. If sis not a positive integer, Ws,p(D) is defined by interpolation [1]. In the case s= 0 one has W0,p(D) = Lp(D). As it is standard, Ws,p(D)dwill be endowed 4B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI with the product norm and, since no confusion can arise, it will be denoted again by k·kWs,p(D). The case p= 2 will be distinguished by using Hs(D) to denote the space Ws,2(D). The space H1 0(D) is the closure in H1(D) of the set of infinitely differentiable functions with compact support in D. For simplicity, k·ks(resp. |·|s) is used to denote the norm (resp. semi norm) both in Hs(Ω) or Hs(Ω)d. The exact meaning will be clear by the context. The inner product of L2(Ω) or L2(Ω)dwill be denoted by (·,·) and the corresponding norm by k·k0in general Dis skipped in the notation for the norm when D= Ω. For vector-valued functions, the same conventions will be used as before. The norm of the dual space H−1(Ω) of H1 0(Ω) is denoted by k · k−1. As usual, L2(Ω) is always identified with its dual, so one has H1 0(Ω) ⊂L2(Ω) ⊂H−1(Ω) with compact injection. The following Sobolev’s embedding [1] will be used in the analysis: For s > 0, let 1 ≤p < d/s and qbe such that 1 q=1 p−s d. Then, there exists a positive scale invariant constant cssuch that (2.1) kvkLq0(Ω) ≤cs|Ω|s d−1 p+1 q0kvkWs,p(Ω),1 q0≥1 q,∀v∈Ws,p(Ω). If p > d/s the above relation is valid for q0=∞. A similar embedding inequality holds for vector-valued functions. We will also use the following interpolation inequality (see, e.g., [14, formula (6.7)] and [22, Exercise II.2.9]) (2.2) kvkL2d/(d−1)(Ω) ≤c1kvk1/2 0kvk1/2 1,∀v∈H1(Ω), (where, for simplicity, by enlarging the constants if necessary, we may take the constant c1in (2.2) equal to csin (2.1) for s= 1) and Agmon’s inequality (2.3) kvk∞≤cAkvk1/2 d−2kvk1/2 2, d = 2,3,∀v∈H2(Ω). The case d= 2 is a direct consequence of [2, Theorem 3.9]. For d= 3, a proof for domains of class C2can be found in [14, Lemma 4.10]. By means of the Calder´on extension theorem (see e.g., [1, Theorem 4.32] the proof is also valid for bounded Lipschitz domains. Finally, we will use Poincar´e’s inequality, (2.4) kvk0≤cP|Ω|1/dk∇vk0,∀v∈H1 0(Ω), where the constant cPcan be taken cP≤√2/2. Denoting by (2.5) ˆcP= 1 + c2 P, observe that from (2.4) it follows that (2.6) kvk1≤(ˆcP)1/2k∇vk0,∀v∈H1 0(Ω). In all previous inequalities, the constants cs,c1,cAand cPare scale-invariant, as it will be the case of all constants in the present paper unless explicitly stated otherwise. Let Hand Vbe the Hilbert spaces H={u∈L2(Ω))d|div(u)=0,u·n|∂Ω= 0}, V={u∈H1 0(Ω))d|div(u) = 0}, endowed with the inner product of L2(Ω)dand H1 0(Ω)d,respectively. Let Th= (τh j, φh j)j∈Jh,h > 0 be a family of partitions of suitable domains Ωh, where his the maximum diameter of the elements τh j∈ Th, and φh jare the mappings ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 5 from the reference simplex τ0onto τh j. We shall assume that the partitions are shaperegular and quasi-uniform. Let r≥2, we consider the finite-element spaces Sh,r =nχh∈ CΩhχh|τh j◦φh j∈Pr−1(τ0)o⊂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. For r= 1, Sh,1stands for the space of piecewise constants. When Ω has polygonal or polyhedral boundary Ωh= Ω and mappings φh jfrom the reference simplex are affine. When Ω has a smooth boundary, for the purpose of analysis we will assume that Ωhexactly matches Ω, as it is done for example in [12], [43], although at a price of a more complex analysis discrepancies between Ωhand Ω can also be taken into account (see, e.g., [6], [44]). We shall denote by (Xh,r, Qh,r−1) the MFE pair known as Hood–Taylor elements [10,46], when r≥3, where Xh,r =S0 h,rd, Qh,r−1=Sh,r−1∩L2(Ωh)/R, r ≥3, and, when r= 2, the MFE pair known as the mini-element [11] where Qh,1=Sh,2∩ L2(Ωh)/R, and Xh,2= (S0 h,2)d⊕Bh. Here, Bhis spanned by the bubble functions bτ,τ∈ Th, defined by bτ(x) = (d+ 1)d+1λ1(x)···λd+1(x), if x∈τand 0 elsewhere, where λ1(x), . . . , λd+1(x) denote the barycentric coordinates of x. For these elements a uniform inf-sup condition is satisfied (see [10]), that is, there exists a constant βis >0 independent of the mesh grid size hsuch that (2.7) inf qh∈Qh,r−1 sup vh∈Xh,r (qh,∇·vh) kvhk1kqhkL2/R≥βis. The velocity will be approximated by elements of the discrete divergence-free space Vh,r =Xh,r ∩χh∈H1 0(Ωh)d|(qh,∇·χh)=0 ∀qh∈Qh,r−1. For each fixed time t∈[0, T] the solution (u, p) of (1.1) is also the solution of a Stokes problem with right-hand side f−ut−(u·∇)u. We will denote by (sh, qh)∈ (Xh,r, Qh,r−1),its MFE approximation satisfying ν(∇sh,∇ϕh)−(qh,∇·ϕh) = ν(∇u, ∇ϕh)−(p, ∇·ϕh) = (f−ut−(u·∇u),ϕh)∀ϕh∈Xh,r,(2.8) (∇·sh, ψh) = 0 ∀ψh∈Qh,r−1. We observe that sh=Sh(u) : V→Vh,r is the discrete Stokes projection of the solution (u, p) of (1.1) (see [29]) and satisfies ν(∇Sh(u),∇ϕh) = ν(∇u,∇ϕh)−(p, ∇·ϕh)=(f−ut−(u·∇)u,ϕh),∀ϕh∈Vh,r. The following bound holds: (2.9) ku−shk0+hku−shk1≤CNj(u, p)hj,1≤j≤r, where here and in the sequel, for v∈V∩Hj(Ω)dand q∈L2 0(Ω)∩Hj−1(Ω) we denote (2.10) Nj(v, q) = kvkj+ν−1kqkHj−1/R, j ≥1. 6B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI The proof of (2.9) when Ω = Ωhcan be found in [30, Lemma 5.3]. For the general case, superparametric approximation at the boundary is assumed, see [5]. Under the same conditions, the bound for the pressure is (cf. [27]) (2.11) kp−qhkL2/R≤Cβis νNj(u, p)hj−1,1≤j≤r, where the constant Cβis depends on the constant βis in (2.7). Assuming that Ω is of class Cm, with m≥3, and using standard duality arguments and (2.9), one obtains (2.12) ku−shk−s≤CNr(u, p)hr+s,0≤s≤min(r−2,1). We also consider a modified Stokes projection that was introduced in [15] and that we denote by sm h:V→Vh,r satisfying ν(∇sm h,∇ϕh)=(f−ut−(u·∇)u−∇p, ϕh),∀ϕh∈Vh,r.(2.13) The following bound holds, see [15]: (2.14) ku−sm hk0+hku−sm hk1≤Ckukjhj,1≤j≤r. Following [12], one can also obtain the following bound k∇(u−sm h)k∞≤Ck∇uk∞,(2.15) where Cdoes not depend on ν. We will denote by πhpthe L2projection of the pressure ponto Qh,r−1. It holds (2.16) kp−πhpk0≤Chj−1kpkHj−1/R,1≤j≤r. If the family of meshes is quasi-uniform then the following inverse inequality holds for each vh∈Sh,r, see e.g., [13, Theorem 3.2.6], (2.17) kvhkWm,p(K)≤cinvhn−m−d(1 q−1 p) KkvhkWn,q(K), where 0 ≤n≤m≤1, 1 ≤q≤p≤ ∞, and hKis the diameter of K∈ Th. In the sequel ILa hu∈Xh,r will denote the Lagrange interpolant of a continuous function u. The following bound can be found in [9, Theorem 4.4.4] (2.18) |u−ILa hu|Wm,p(K)≤cinthn−m|u|Wn,p(K),0≤m≤n≤k+ 1, where n > d/p when 1 < p ≤ ∞ and n≥dwhen p= 1. We will assume that the interpolation operator IHis stable in L2, that is, kIHuk0≤c0kuk0,∀u∈L2(Ω)d,(2.19) and that it satisfies the following approximation property, ku−IHuk0≤cIHk∇uk0,∀u∈H1 0(Ω)d.(2.20) The Bernardi–Girault [8], Girault–Lions [26], or the Scott–Zhang [45] interpolation operators satisfy (2.20) and (2.19). Notice that the interpolation can be on piecewise constants, as we use in the numerical experiments in Section 4. We remark that, for the error analysis, we do not need condition (3.105) in [37], i.e., we do not assume that ku−IH(u)k−1≤c−1Hkuk0, for u∈L2(Ω)d. ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 7 3. The finite element method. We consider the following method to approximate (1.2). Find (uh, ph)∈Xh,r ×Qh,r−1satisfying for all (ϕh, ψh)∈Xh,r ×Qh,r−1 (˙ uh,ϕh) + ν(∇uh,∇ϕh) + bh(uh,uh,ϕh) + µ(∇·uh,∇·ϕh)+(∇ph,ϕh) = (f,ϕh)−β(IH(uh)−IH(u), IHϕh),(3.1) (∇·uh, ψh) = 0, where µis a stabilization parameter that can be zero in case we do not stabilize the divergence or different from zero in case we add grad-div stabilization 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. Hereafter, we denote by (·,·) both the inner product in L2and the duality action between H−1and H1 0, depending on the context. It is straightforward to verify that bhenjoys the skew-symmetry property (3.2) bh(u,v,w) = −bh(u,w,v)∀u,v,w∈H1 0(Ω)d. Let us observe that taking ϕh∈Vh,r from (3.1) we get (˙ uh,ϕh) + ν(∇uh,∇ϕh) + bh(uh,uh,ϕh) + µ(∇·uh,∇·ϕh) =(3.3) (f,ϕh)−β(IH(uh)−IH(u), IHϕh). For the analysis below, we need to introduce the values µand k, defined as follows (3.4) µ=0,if µ= 0, 1,otherwise,k=0,if µ= 0, 1/µ, otherwise. The following lemma will be used for proving the main results of the section. Lemma 3.1. Let uhbe the finite element approximation defined in (3.3)and let wh,τ1 h,τ2 h: [0, T]→Vh,r be functions satisying (˙ wh,ϕh) + ν(∇wh,∇ϕh) + bh(wh,wh,ϕh) + µ(∇·wh,∇·ϕh) =(3.5) (f,ϕh)+(τ1 h,ϕh) + µ(τ2 h,∇·ϕh), Assume that the quantity Ldefined in (3.14), below, when µ= 0, and in (3.15), below, when µ > 0is bounded. Then, if β≥8Land Hsatisfies condition (3.21), below, the following bounds hold for eh=uh−wh, keh(t)k2 0≤e−γt/2keh(0)k2 0+Zt 0 e−γ(t−s)/2(1 −µ)2ˆcP ν+µ Lkτ1 hk2 −1+µds +Zt 0 e−γ(t−s)/2βc2 0ku(s)−wh(s)k2 0+ 2kkτ2 hk2 0ds,(3.6) where, µand kare defined in (3.4), and γis defined in (3.24)below. Proof. Subtracting (3.5) from (3.3) we get the error equation (˙ eh,ϕh) + ν(∇eh,∇ϕh) + β(IHeh, IHϕh) + bh(uh,uh,ϕh)−bh(wh,wh,ϕh)(3.7) +µ(∇·eh,∇·ϕh) = β(IHu−IHwh, IHϕh)+(τ1 h,ϕh) + µ(τ2 h,∇·ϕh),∀ϕh∈Vh,r. 8B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI Taking ϕh=ehin (3.7) we get 1 2 d dtkehk2 0+νk∇ehk2 0+βkIHehk2 0+µk∇·ehk2 0≤ |bh(uh,uh,eh)−bh(wh,wh,eh)| +β|(IHu−IHwh, IHeh)|+|(τ1 h,eh)|+|µ(τ2 h,∇·eh)|.(3.8) We will bound the terms on the right-hand side of (3.8). For the nonlinear term and the truncation errors we argue differently depending on whether µ= 0 or µ > 0. If µ= 0, using the skew-symmetry property (3.2), (2.1) and (2.6), and when d= 3, we have |bh(uh,uh,eh)−bh(wh,wh,eh)|=|bh(eh,wh,eh)| ≤ k∇whkL2d/(d−1) kehkL2dkehk0 +1 2|(∇·eh)wh,eh)| ≤ ˆc1k∇whkL2d/(d−1) k∇ehk0kehk0+1 2k∇ehk0kwhk∞kehk0 ≤2ˆc2 1k∇whk2 L2d/(d−1) ν+kwhk2 ∞ νkehk2 0+ν 4k∇ehk2 0, (3.9) where (3.10) ˆc1= (ˆcP)1/2c1, c1being the constant in (2.1) for s= 1. In the case d= 2, and noticing that 2d= 2d/(d−1), the first term on the right-hand side above, using (2.2), (2.6) and Young’s inequality is bounded as follows k∇whkL2d/(d−1) kehkL2dkehk0≤(ˆcP)1/4c1k∇whkL2d/(d−1) (k∇ehk0kehk0)1/2kehk0 ≤3(ˆcP)1/3c4/3 1k∇whk4/3 L2d/(d−1) (4ν)1/3kehk2 0+ν 4k∇ehk2 0.(3.11) For the truncation error when µ= 0 using (2.6) we get (3.12) |(τ1 h,eh)| ≤ kτ1 hk−1kehk1≤(ˆcP)1/2kτ1 hk−1k∇ehk0≤ˆcP νkτ1 hk2 −1+ν 4k∇ehk2 0. When µ6= 0, we bound the nonlinear term in the following way. Using again the skew-symmetry property (3.2) we get |bh(uh,uh,eh)−bh(wh,wh,eh)|=|bh(eh,wh,eh)| ≤ k∇whk∞kehk2 0 +1 2k∇·ehk0kwhk∞kehk0≤ k∇whk∞kehk2 0+µ 4k∇·ehk2 0+kwhk2 ∞ 4µkehk2 0.(3.13) In the sequel we denote L= max t≥0 2ˆc2d/3 1k∇wh(t)k2d/3 L2d/(d−1) ν(2d−3)/3+kwh(t)k2 ∞ ν!,if µ= 0,(3.14) L= 2 max t≥0k∇wh(t)k∞+kwh(t)k2 ∞ 4µ,if µ > 0,(3.15) Observe that in the case µ= 0, bounding the factor 3(ˆcP)1/3/41/3in (3.11) by 2(ˆcP)2/3 we have the left-hand side of (3.9) can be bounded by Lkehk2 0+(ν/2)k∇ehk2 0, and, in the case µ > 0 the left-hand side of (3.13) is bounded by (L/2)kehk2 0+(µ/4)k∇·ehk2 0. ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 9 Next, we bound the truncation error when µ > 0, |(τ1 h,eh)|+|µ(τ2 h,∇·eh)| ≤ 1 2Lkτ1 hk2 0+L 2kehk2 0+kkτ2 hk2 0+µ 4k∇·ehk2 0,(3.16) where kis defined in (3.4). For the second term on the right-hand side of (3.8) applying (2.19) we get β|(IHu−IHwh, IHeh)| ≤ βc0ku−whk0kIHehk0 ≤β 2c2 0ku−whk2 0+β 2kIHehk2 0.(3.17) Inserting (3.9), (3.11), (3.12), (3.13), (3.16) and (3.17) into (3.8) we get 1 2 d dtkehk2 0+(1 + µ)ν 2k∇ehk2 0+β 2kIHehk2 0+µ 2k∇·ehk2 0≤Lkehk2 0 (3.18) +kkτ2 hk2 0+β 2c2 0ku−whk2 0+(1 −µ)ˆcP ν+µ 2Lkτ1 hk2 −1+µ. Now we bound Lkehk2 0≤2LkIHehk2 0+ 2Lk(I−IH)ehk2 0. Since we are assuming that β≥8Lwe have that β/2−2L≥β/4, so that taking into account that 1 + µ≥1 and (µ/2)k∇·ehk ≥ 0 we get d dtkehk2 0+νk∇ehk2 0−4Lk(I−IH)ehk2 0+β 2kIHehk2 0≤(3.19) 2kkτ2 hk2 0+βc2 0ku−whk2 0+(1 −µ)2ˆcP ν+µ Lkτ1 hk2 −1+µ. For the second and third terms on the left-hand side above, applying (2.20) to the latter, we write νk∇ehk2 0−4Lk(I−IH)ehk2 0≥νk∇ehk2 0−4Lc2 IH2k∇ehk2 0≥ν 2k∇ehk2 0,(3.20) whenever (3.21) H≤ν1/2 (8L)1/2cI . Therefore, for the last three terms on the left-hand side of (3.19) we have (3.22) νk∇ehk2 0+β 2kIHehk2 0−4Lk(I−IH)ehk2 0≥ν 2k∇ehk2 0+β 2kIHehk2 0. Now, applying (2.20) again to bound below the right-hand side above we have that ν 2k∇ehk2 0+β 2kIHehk2 0≥ν 2c−2 IH−2k(I−IH)ehk2 0+β 2kIHehk2 0 (3.23) ≥γ(kIHehk2 0+k(I−IH)ehk2 0), where (3.24) γ= min ν 2c−2 IH−2,β 2. 16 B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI spatial discretization we used P2/P1elements on a regular triangulation with SW-NE diagonals, with the same number of subdivisions on each coordinate direction. For for coarse mesh interpolation we take piecewise constants. The time integration was done with an implicit/explicit (IMEX) method based on the second order backward differentiation formula (BDF), where, to avoid solving nonlinear steady problems at each step, linear extrapolation of the form bh(2uh(t−∆t)−uh(t−2∆t),uh(t),ϕh) was used in the convection term, where ∆tis the time step, except in the first step where bh(uh(t−∆t),uh(t),ϕh) was used. The time step was chosen so that the error arising from the spatial discretization was dominant. To check that this was the case, we made sure that results were not essentially altered if recomputed with a smaller ∆t. For h= 1/12,1/18,1/24,1/36 and 1/48, the values of ∆tused were, respectively, ∆t= 1/160,1/160,1/320,1/320 and 1/640. Unless stated otherwise, in what follows the initial condition was set to uh=0and p= 0, so that there is an O(1) error at time t= 0. We first check that there is no upper bound on the nudging parameter β. The left plot in Fig. 1shows the velocity errors in L2vs time for different values of β for ν= 10−6, including β= 100. It can be seen a clear difference between β= 0, where the initial errors do not decay with time, and β > 0 where they do, and for the four largest values of βshown, they do so exponentially in time, until an asymptotic regime is reached. We also notice that the results are little altered for β≥10. In 0 10 20 30 40 50 60 10-3 10-2 10-1 100 101 0 10 20 30 40 50 60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fig. 1.Velocity errors vs time view of (3.24), one may be temped to question the advantage of taking β≥ν(cIH)−2, since the rate of decay of the initial errors, γ, is unaltered for larger values of β. If we assume that cIis of order one, then, the value of ν(cIH)−2in the present example is unlikely to be larger than 10−5, so that Fig. 1(and more examples in [34] and [40]) seems to suggest that there is some advantage in taking β≥ν(cIH)−2if we want a faster decay of the initial errors. Since this is in apparent contradiction with the analysis in the previous section, we now propose an alternative explanation. Let us consider for some integer k≥2 the value r=kIH(u(0))k0/(kku(0)k0). If we take the intial condtion uh= 0, then, by continuity there exist t0>0 such that (4.3) kIH(eh(t)k0≥rkeh(t)k0, t ∈[0, t0]. ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 17 Consequently, from (3.18) it follows that d dtkehk2 0+(βr2−2L)kehk2 0≤2kkτ2 hk2 0+βc2 0ku−whk2 0+(1 −µ)2ˆcP ν+µL kτ1 hk2 −1+µ, for t∈[0, t0], which, for β > 2L/r2, could explain that initial errors decay faster when larger values of βare taken. In Fig. 1we also show the ratios kIH(eh)k0/kehk0. It can be seen that although they became smaller as βis increased, they are sufficiently away from zero to suggest that the analysis in the present section may explain the faster rates of decay of the initial errors when larger values of βare taken. In the reminder of this section we take β= 1. Since, as shown in Fig. 1, after an initial decay, the errors show an oscillatory behaviour, in the examples below by L2errors we mean the maximum of errors kuh(t)−u(t)k0for values of tafter the asymptotic regime has shown itself. We now check the rates of convergence proved in the present paper. In Fig. 2we show errors vs hfor different values of the diffusion parameter νand compare the cases of positive µ(µ= 0.05) and µ= 0. The value of His H= 3hand βis set to β= 1. Results corresponding to the smallest value of νare represented with discontinuous lines in both plots so that they can be seen superimposed to those corresponding to larger values of ν. Slopes of least squares fits to the results corresponding to each value of νare shown, so that the order of convergence can be checked. In both cases, µ= 0.05 and µ= 0, O(h3) errors are obtained for large values of ν, which is what Theorem 3.2 and Remark 3.5 predict. However, for smaller values of ν, while the errors with positive µbecome O(h2) and independent of ν, as Theorem 3.3 predicts, for µ= 0 the method does not have convergent behaviour for the values of hshown (presumably, the method will show convergence for h≤ν). 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 10-4 10-3 10-2 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 10-4 10-3 10-2 10-1 100 101 Fig. 2.Velocity errors. Left, µ= 0.05 Right, µ= 0. Finally, we check that the requirement H/h bounded is required for convergence if Lagrange interpolants are used. In Fig. 3we show velocity errors when h→0 in two different scenarios: h=H/3 (left) and Hfixed to H= 0.25. We see that while H= 3h, the method converge as predicted by Theorem 3.12 (the value of ν= 10−6 and that of µ= 0.05). If His kept fixed, however, the method using the Lagrange Interpolant does not exhibit convergent behaviour. We remark, however, that with larger values of βor ν, convergence is not altered as much as in Fig. 3when H/h grows. Nevertheless, this example shows the risks of not keeping (H/h) bounded with Lagrange interpolants. 18 B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI 0.03 0.04 0.05 0.06 0.07 0.08 10-3 10-2 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 10-2 Fig. 3.Velocity errors. Left H= 3h. Right H= 0.25. 5. Conclusions. We have analyzed a semidiscretization in space by inf-sup stable mixed finite elements of a continuous downscaling data assimilation method for the two and three-dimensional Navier-Stokes equations. The data assimilation method, introduced in [7], combines observational data (measurements) on large spatial scales or coarse mesh, IHu, with simulations in order to improve predictions of the physical phenomenon being studied. We have considered the Galerkin method with and without grad-div stabilization. Uniform error bounds in time have been obtained for the approximation to velocity field, under standard assumptions in finite element analysis. The order of convergence proved for the method without stabilization is optimal, in the sense that it is the best that can be obtained with the finite element space being used (i.e., errors of the same order as interpolation). For the Galerkin method with grad-div stabilization error bounds in which the constants are independent on inverse powers of the viscosity are proved. Convergence rates and dependence or independence of νare corroborated in numerical experiments. As opposed to previous works in the literature, our analysis also covers the case in which IHuis the standard Lagrange interpolant, where we show that H/h must be kept bounded in order to get convergence. Also, the upper bound on the nudging parameter assumed in previous references is removed. The techniques of analysis used in the present paper allow to improve the available error bounds for a closely-related finite element method in [34]. REFERENCES [1] R. A. Adams and J. J. F. Fournier,Sobolev spaces, vol. 140 of Pure and Applied Mathematics (Amsterdam), Elsevier/Academic Press, Amsterdam, second ed., 2003. [2] S. Agmon,Lectures on elliptic boundary value problems, AMS Chelsea Publishing, Providence, RI, 2010, https://doi.org/10.1090/chel/369,https://doi.org/10.1090/chel/369. Prepared for publication by B. Frank Jones, Jr. with the assistance of George W. Batten, Jr., Revised edition of the 1965 original. [3] M. U. Altaf, E. S. Titi, T. Gebrael, O. M. Knio, L. Zhao, M. F. McCabe, and I. Hoteit,Downscaling the 2d b´enard convection equations using continuous data assimilation, Computational Geosciences, 21 (2017), pp. 393–410, https://doi.org/10.1007/ s10596-017-9619-2. [4] M. Asch, M. Bocquet, and M. Nodet,Data assimilation, vol. 11 of Fundamentals of Algorithms, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2016, https://doi.org/10.1137/1.9781611974546.pt1,https://doi.org/10.1137/1.9781611974546. pt1. Methods, algorithms, and applications. ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 19 [5] B. Ayuso and B. Garc´ ıa-Archilla,Regularity constants of the Stokes problem. Application to finite-element methods on curved domains, Math. Models Methods Appl. Sci., 15 (2005), pp. 437–470, https://doi.org/10.1142/S021820250500042X,https://doi.org/10. 1142/S021820250500042X. [6] B. Ayuso, B. Garc´ ıa-Archilla, and J. Novo,The postprocessed mixed finite-element method for the Navier-Stokes equations, SIAM J. Numer. Anal., 43 (2005), pp. 1091–1111, https: //doi.org/10.1137/040602821,https://doi.org/10.1137/040602821. [7] A. Azouani, E. Olson, and E. S. Titi,Continuous data assimilation using general interpolant observables, J. Nonlinear Sci., 24 (2014), pp. 277–304, https://doi.org/10.1007/ s00332-013-9189-y,https://doi.org/10.1007/s00332-013-9189-y. [8] C. Bernardi and V. Girault,A local regularization operator for triangular and quadrilateral finite elements, SIAM J. Numer. Anal., 35 (1998), pp. 1893–1916, https://doi.org/10.1137/ S0036142995293766,https://doi.org/10.1137/S0036142995293766. [9] S. C. Brenner and L. R. Scott,The mathematical theory of finite element methods, vol. 15 of Texts in Applied Mathematics, Springer-Verlag, New York, second ed., 2002, https: //doi.org/10.1007/978-1-4757-3658-8,https://doi.org/10.1007/978-1-4757-3658-8. [10] F. Brezzi and R. S. Falk,Stability of higher-order Hood-Taylor methods, SIAM J. Numer. Anal., 28 (1991), pp. 581–590, https://doi.org/10.1137/0728032,https://doi.org/10.1137/ 0728032. [11] F. Brezzi and M. Fortin,Mixed and hybrid finite element methods, vol. 15 of Springer Series in Computational Mathematics, Springer-Verlag, New York, 1991, https://doi.org/10.1007/ 978-1-4612-3172-1,https://doi.org/10.1007/978-1-4612-3172-1. [12] H. Chen,Pointwise error estimates for finite element solutions of the Stokes problem, SIAM J. Numer. Anal., 44 (2006), pp. 1–28, https://doi.org/10.1137/S0036142903438100,https: //doi.org/10.1137/S0036142903438100. [13] P. G. Ciarlet,The finite element method for elliptic problems, vol. 40 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002, https://doi.org/10.1137/1.9780898719208,https://doi.org/10.1137/1. 9780898719208. Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58 #25001)]. [14] P. Constantin and C. Foias,Navier-Stokes equations, Chicago Lectures in Mathematics, University of Chicago Press, Chicago, IL, 1988. [15] J. de Frutos, B. Garc´ ıa-Archilla, V. John, and J. Novo,Grad-div stabilization for the evolutionary Oseen problem with inf-sup stable finite elements, J. Sci. Comput., 66 (2016), pp. 991–1024, https://doi.org/10.1007/s10915-015-0052-1,https://doi.org/10. 1007/s10915-015-0052-1. [16] J. de Frutos, B. Garc´ ıa-Archilla, V. John, and J. Novo,Analysis of the grad-div stabilization for the time-dependent Navier-Stokes equations with inf-sup stable finite elements, Adv. Comput. Math., 44 (2018), pp. 195–225, https://doi.org/10.1007/s10444-017-9540-1, https://doi.org/10.1007/s10444-017-9540-1. [17] J. de Frutos, B. Garc´ ıa-Archilla, and J. Novo,Error analysis of projection methods for non inf-sup stable mixed finite elements: the Navier-Stokes equations, J. Sci. Comput., 74 (2018), pp. 426–455, https://doi.org/10.1007/s10915-017-0446-3,https://doi.org/10.1007/ s10915-017-0446-3. [18] A. Farhat, H. Johnston, M. Jolly, and E. S. Titi,Assimilation of nearly turbulent rayleigh– b´enard flow through vorticity or local circulation measurements: A computational study, Journal of Scientific Computing, 77(3) (2018), pp. 1519–1533, https://doi.org/10.1007/ s10915-018-0686-x. [19] A. Farhat, M. S. Jolly, and E. S. Titi,Continuous data assimilation for the 2D B´enard convection through velocity measurements alone, Phys. D, 303 (2015), pp. 59–66, https: //doi.org/10.1016/j.physd.2015.03.011,https://doi.org/10.1016/j.physd.2015.03.011. [20] A. Farhat, E. S. Lunasin, and E. S. Titi,On the charney conjecture of data assimilation employing temperature measurements alone: the paradigm of 3d planetary geostrophic model,, Math. Clim. Weather Forecast., 2 (2016), pp. 59–66, https://doi.org/10.1515/ mcwf-2016-0004,https://doi.org/10.1515/mcwf-2016-0004. [21] L. P. Franca and T. J. R. Hughes,Two classes of mixed finite element methods, Comput. Methods Appl. Mech. Engrg., 69 (1988), pp. 89–129, https://doi.org/10.1016/ 0045-7825(88)90168-5,https://doi.org/10.1016/0045-7825(88)90168-5. [22] G. P. Galdi,An introduction to the mathematical theory of the Navier-Stokes equations. Vol. I, vol. 38 of Springer Tracts in Natural Philosophy, Springer-Verlag, New York, 1994, https: //doi.org/10.1007/978-1-4612-5364-8,https://doi.org/10.1007/978-1-4612-5364-8. Linearized steady problems. 20 B. GARC´ IA-ARCHILLA, J. NOVO, AND E. S. TITI [23] B. Garc´ ıa-Archilla, J. Novo, and E. S. Titi,Postprocessing the Galerkin method: a novel approach to approximate inertial manifolds, SIAM J. Numer. Anal., 35 (1998), pp. 941–972, https://doi.org/10.1137/S0036142995296096,https://doi.org/10. 1137/S0036142995296096. [24] B. Garc´ ıa-Archilla, J. Novo, and E. S. Titi,An approximate inertial manifolds approach to postprocessing the Galerkin method for the Navier-Stokes equations, Math. Comp., 68 (1999), pp. 893–911, https://doi.org/10.1090/S0025-5718-99-01057-1,https://doi.org/10. 1090/S0025-5718-99-01057-1. [25] M. Gesho, E. Olson, and E. S. Titi,A computational study of a data assimilation algorithm for the two-dimensional Navier-Stokes equations, Commun. Comput. Phys., 19 (2016), pp. 1094–1110. [26] V. Girault and J.-L. Lions,Two-grid finite-element schemes for the transient Navier-Stokes problem, M2AN Math. Model. Numer. Anal., 35 (2001), pp. 945–980, https://doi.org/10. 1051/m2an:2001145,https://doi.org/10.1051/m2an:2001145. [27] V. Girault and P.-A. Raviart,Finite element methods for Navier-Stokes equations, vol. 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1986, https://doi. org/10.1007/978-3-642-61623-5,https://doi.org/10.1007/978-3-642-61623-5. Theory and algorithms. [28] K. Hayden, E. Olson, and E. S. Titi,Discrete data assimilation in the Lorenz and 2D Navier-Stokes equations, Phys. D, 240 (2011), pp. 1416–1425, https://doi.org/10.1016/j. physd.2011.04.021,https://doi.org/10.1016/j.physd.2011.04.021. [29] J. G. Heywood and R. 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 (1982), pp. 275–311, https://doi.org/10. 1137/0719018,https://doi.org/10.1137/0719018. [30] J. G. Heywood and R. Rannacher,Finite element approximation of the nonstationary Navier-Stokes problem. III. Smoothing property and higher order error estimates for spatial discretization, SIAM J. Numer. Anal., 25 (1988), pp. 489–512, https://doi.org/10. 1137/0725032,https://doi.org/10.1137/0725032. [31] H. A. Ibdah, C. F. Mondaini, and E. S. Titi,Uniform in time error estimates for fully discrete numerical schemes of a data assimilation algorithm, IMA Journal of Numerical Analysis, to appear (2019), https://arxiv.org/abs/1805.01595v1. [32] V. John and A. Kindl,Numerical studies of finite element variational multiscale methods for turbulent flow simulations, Comput. Methods Appl. Mech. Engrg., 199 (2010), pp. 841–852, https://doi.org/10.1016/j.cma.2009.01.010,https://doi.org/10.1016/j.cma.2009.01.010. [33] E. Kalnay,Atmospheric Modeling, Data Assimilation and Predictability, Cambridge University Press, 2002, https://doi.org/10.1017/CBO9780511802270. [34] A. Larios, L. G. Rebholz, and C. Zerfas,Global in time statbility and accuracy of imex-fem data assimilation schemes for the navier-stokes equations, arXiv:1805.04090v1, (2018), https://arxiv.org/abs/1805.04090. [35] K. Law, A. Stuart, and K. Zygalakis,Data assimilation, vol. 62 of Texts in Applied Mathematics, Springer, Cham, 2015, https://doi.org/10.1007/978-3-319-20325-6,https: //doi.org/10.1007/978-3-319-20325-6. A mathematical introduction. [36] E. Lunasin and E. S. Titi,Finite determining parameters feedback control for distributed nonlinear dissipative systems—a computational study, Evol. Equ. Control Theory, 6 (2017), pp. 535–557, https://doi.org/10.3934/eect.2017027,https://doi.org/10.3934/eect.2017027. [37] C. F. Mondaini and E. S. Titi,Uniform-in-time error estimates for the postprocessing Galerkin method applied to a data assimilation algorithm, SIAM J. Numer. Anal., 56 (2018), pp. 78–110, https://doi.org/10.1137/16M110962X,https://doi.org/10.1137/ 16M110962X. [38] M. A. Olshanskii,A low order Galerkin finite element method for the Navier-Stokes equations of steady incompressible flow: a stabilization issue and iterative methods, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 5515–5536, https://doi.org/10.1016/ S0045-7825(02)00513-3,https://doi.org/10.1016/S0045-7825(02)00513-3. [39] M. A. Olshanskii and A. Reusken,Grad-div stabilization for Stokes equations, Math. Comp., 73 (2004), pp. 1699–1718, https://doi.org/10.1090/S0025-5718-03-01629-6,https: //doi.org/10.1090/S0025-5718-03-01629-6. [40] L. G. Rebholz and C. Zerfas,Simple and efficient continuous data assimilation of evolution equations via algebraic nudging, arXiv e-prints, (2018), arXiv:1810.03512, p. arXiv:1810.03512, https://arxiv.org/abs/1810.03512. [41] S. Reich and C. Cotter,Probabilistic forecasting and Bayesian data assimilation, Cambridge University Press, New York, 2015, https://doi.org/10.1017/CBO9781107706804,https:// ERROR ESTIMATES FOR A FEM APPLIED TO A DOWNSCALING DA ALGORITHM 21 doi.org/10.1017/CBO9781107706804. [42] L. R¨ ohe and G. Lube,Analysis of a variational multiscale method for large-eddy simulation and its application to homogeneous isotropic turbulence, Comput. Methods Appl. Mech. Engrg., 199 (2010), pp. 2331–2342, https://doi.org/10.1016/j.cma.2010.02.015,https:// doi.org/10.1016/j.cma.2010.02.015. [43] A. H. Schatz,Pointwise error estimates and asymptotic error expansion inequalities for the finite element method on irregular grids. I. Global estimates, Math. Comp., 67 (1998), pp. 877–899, https://doi.org/10.1090/S0025-5718-98-00959-4,https://doi.org/10. 1090/S0025-5718-98-00959-4. [44] A. H. Schatz and L. B. Wahlbin,On the quasi-optimality in L∞of the ˙ H1-projection into finite element spaces, Math. Comp., 38 (1982), pp. 1–22, https://doi.org/10.2307/2007461, https://doi.org/10.2307/2007461. [45] L. R. Scott and S. Zhang,Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp., 54 (1990), pp. 483–493, https://doi.org/10.2307/ 2008497,https://doi.org/10.2307/2008497. [46] C. Taylor and P. Hood,A numerical solution of the Navier-Stokes equations using the finite element technique, Internat. J. Comput. & Fluids, 1 (1973), pp. 73–100, https://doi.org/ 10.1016/0045-7930(73)90027-3,https://doi.org/10.1016/0045-7930(73)90027-3.