A streamline derivative POD-ROM for advection-diffusion-reaction equations

ESAIM: PROCEEDINGS AND SURVEYS, October 2018, Vol. 64, p. 121-136 Laurence CARASSUS, Marion DARBAS, Ghislaine GAYRAUD, Olivier GOUBET, St´ephanie SALMON A STREAMLINE DERIVATIVE POD-ROM FOR ADVECTION-DIFFUSION-REACTION EQUATIONS ∗ Samuele Rubino1 Abstract. We introduce a new streamline derivative projection-based closure modeling strategy for the numerical stabilization of Proper Orthogonal Decomposition-Reduced Order Models (PODROM). As a first preliminary step, the proposed model is analyzed and tested for advection-dominated advection-diffusion-reaction equations. In this framework, the numerical analysis for the Finite Element (FE) discretization of the proposed new POD-ROM is presented, by mainly deriving the corresponding error estimates. Numerical tests for advection-dominated regime show the efficiency of the proposed method, as well the increased accuracy over the standard POD-ROM that discovers its well-known limitations very soon in the numerical settings considered, i.e. for low diffusion coefficients. R´esum´e. Nous introduisons une nouvelle strat´egie de mod´elisation de type streamline derivative bas´ee sur projection pour la stabilisation num´erique de mod`eles d’ordre r´eduit de type POD (PODROM). Comme premi`ere ´etape pr´eliminaire, le mod`ele propos´e est analys´e et test´e pour les ´equations d’advection-diffusion-r´eaction domin´ees par l’advection. Dans ce cadre, l’analyse num´erique de la discr´etisation par ´el´ements finis (FE) du nouveau POD-ROM propos´e est pr´esent´ee, en d´erivant principalement les estimations d’erreur correspondantes. Des tests num´eriques pour le r´egime domin´e par l’advection montrent l’efficacit´e de la m´ethode propos´ee, ainsi que la pr´ecision accrue par rapport `a la m´ethode POD-ROM standard qui d´ecouvre tr`es rapidement ses limites bien connues dans le cas des param`etres num´eriques consid´er´es, c’est-`a-dire pour de faibles coefficients de diffusion. Introduction Among the most popular Reduced Order Models (ROM) approaches, Proper Orthogonal Decomposition (POD) strategy provides optimal (from the energetic point of view) modes to represent the dynamics from a given database (snapshots) obtained by a full-order system. Onto these POD modes, a Galerkin projection of the governing equations can be employed to obtain a low-order dynamical system for the modes coefficients. The resulting low-order model is named standard POD-ROM, which thus consists in the projection of high-fidelity (full-order) representations of physical problems onto low-dimensional spaces of solutions, with a dramatically reduced dimension. These low-dimensional spaces are capable of capturing the dominant characteristics of the solution, their main advantage being that the computations in the low-dimensional space can be done at a reduced computational cost. In particular, the computational cost in a Direct Numerical Simulation (DNS) of a complex problem could be reduced by several orders of magnitude when POD-ROM is employed. This has led researchers to apply POD-ROM to a variety of physical and engineering problems, including Computational Fluid Dynamics (CFD) problems, see e.g. [10, 14, 19, 27, 30, 37]. Once applied to the physical problem of ∗The author would gratefully acknowledge the financial support received from IdEx (Initiative d’Excellence de l’Universit´e de Bordeaux) International Post-Doc Program during his postdoctoral research involved in this article. 1Departmento EDAN & IMUS, Universidad of Sevilla, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain; e-mail: [email protected] c EDP Sciences, SMAI 2018 Article published online by EDP Sciences and available at https://www.esaim-proc.org or https://doi.org/10.1051/proc/201864121 122 ESAIM: PROCEEDINGS AND SURVEYS interest, POD-ROM can be used to solve engineering problems such as shape optimization [2, 22] and flow control [3,13,21,35], making their effective resolution increasingly affordable in almost real-time. Although POD-ROM can be very computationally efficient and relatively accurate in some configurations, they also present several drawbacks. In this report, we address one of them, namely the numerical instability of a straightforward POD-Galerkin procedure applied to advection-dominated problems. To address this issue, we draw inspiration from the Finite Element (FE) context, where stabilized formulations have been developed to deal with the numerical instabilities of the Galerkin method. In particular, we consider a Streamline Derivativebased (SD-based) approach used by Knobloch and Lube in [28] in the FE context, which only acts on the high frequencies of the advective derivative (see also [1] for its extension to Navier–Stokes Equations (NSE)). This approach consists in adding a filtered advection stabilization term by basically following the streamlines to prevent spurious instabilities due do dominant advection. This stabilization term acts on the high frequencies component (main responsible for numerical oscillations) of the advection/streamline derivative, which seems to be a natural choice when dealing especially with strongly advection-dominated configurations. Although applications of stabilized methods can already be found in the ROM literature (see [10–12, 20, 25, 26] for the POD context, and also [31,32] for the Reduced-Basis (RB) context), to the authors’ knowledge this is the first time that the SD-based formulation in [28] has been applied in a POD setting. A different strategy used in the recent literature to obtain surrogate ROM for nonlinear dynamical systems is the Dynamic Mode Decomposition (DMD) method. In particular, a DMD-Galerkin method has been applied in [4] to advection-diffusion problems. In this report, the proposed SD-based POD-ROM (SD-POD-ROM) is preliminary analyzed and tested for the numerical approximation of advection-dominated advection-diffusion-reaction problems of the form:    ∂tu+b·∇u−ε∆u+gu =fin Ω ×(0, T), u= 0 on Γ ×(0, T), u(x,0) = u0(x) in Ω, (1) where bwith kbk∞=O(1) is the given advective field, ε << 1 the diffusion parameter, gthe reaction coefficient, fthe forcing term, Ω the computational domain in Rd(d= 2 or 3), t∈[0, T], with Tthe final time, and u0 the initial condition. For the sake of simplicity, we have imposed homogeneous Dirichlet boundary condition on the whole boundary Γ = ∂Ω. Although the new SD-POD-ROM is being developed to derive a low-order approximation of convectiondominated and turbulent flows described by the NSE [9], as a first preliminary step we have decided to analyze it for the mathematical setting in (1), which is simpler to work out, yet relevant to our ultimate goal (since ε << kbk∞). The rest of the paper is organized as follows. In section 1, we briefly describe the POD methodology and introduce the new SD-POD-ROM. The error analysis for the full discretization (FE in space and backward Euler in time) of the new model is presented in section 2. The new method is tested numerically in section 3 for a 2D traveling wave problem, presenting a sharp internal layer moving in time. Finally, section 4 presents the main conclusions of this work and future research directions. 1. Streamline derivative projection-based POD-ROM 1.1. Proper orthogonal decomposition reduced order model For the report to be self-contained, this section briefly presents the computation of a basis for ROM with POD. For more details, the reader is referred to [15,24,33,34,36]. We first present the continuous version of POD method. Consider a function u(x, t):Ω×[0, T]→R, and r∈N. Then, the goal of POD consists in finding the set of orthonormal POD basis {ϕ1, . . . , ϕr}that deliver ESAIM: PROCEEDINGS AND SURVEYS 123 the best approximation: min     u(x, t)− r X i=1 (u(x, t), ϕi)Hϕi     2 L2(0,T ;H) ,(2) in a real Hilbert space H. Although Hcan be any real Hilbert space, in what follows we consider H=L2(Ω), with induced norm k·k = (·,·)1/2=ZΩ|·|21/2 . In the framework of the numerical solution of Partial Differential Equations (PDE), uis usually given at a finite number of times t0, . . . , tN, the so-called snapshots. Let us consider an ensemble of snapshots χ= span {u(·, t0), . . . , u(·, tN)}, which is a collection of data from either numerical simulation results or experimental observations at time tn=n∆t,n= 0,1, . . . , N and ∆t=T/N. Then, usually an approximation of the error in the square of the L2(0, T) norm is consideed, e.g., by a modification of the composite trapezoidal rule. Thus, in its discrete version (method of snapshots), the POD method seeks a low-dimensional basis {ϕ1, . . . , ϕr}that optimally approximates the snapshots in the following sense, see for instance [29]: min 1 N+ 1 N X n=0     u(·, tn)− r X i=1 (u(·, tn), ϕi)ϕi     2 ,(3) subject to the condition (ϕj, ϕi) = δij, 1 ≤i, j ≤r, where δij is the Kronecker delta. To solve the optimization problem (3), one can consider the eigenvalue problem: Kzi=λizi,for 1, . . . , r, (4) where K∈R(N+1)×(N+1) is the snapshots correlation matrix with entries: Kmn =1 N+ 1 (u(·, tn), u(·, tm)) ,for m, n = 0, . . . , N, ziis the i-th eigenvector, and λiis the associated eigenvalue. The eigenvalues are positive and sorted in descending order λ1≥. . . ≥λr>0. It can be shown that the solution of (3), i.e. the POD basis, is given by: ϕi(·) = 1 √λi N X n=0 (zi)nu(·, tn),1≤i≤r, (5) where (zi)nis the n-th component of the eigenvector zi. It can also be shown that the following POD error formula holds [24,29]: 1 N+ 1 N X n=0     u(·, tn)− r X i=1 (u(·, tn), ϕi)ϕi     2 = M X i=r+1 λi,(6) where Mis the rank of χ. We consider the following space for the POD setting: Xr= span {ϕ1, . . . , ϕr}. Remark 1.1. Since, as shown in (5), the POD modes are linear combinations of the snapshots, the POD modes satisfy the boundary conditions in (1). This is because of the particular choice we have made at the beginning to work with homogeneous Dirichlet boundary conditions. In general, one has to manipulate the snapshots set. This is the case, for instance, of steady-state non-homogeneous Dirichlet boundary conditions, for which is preferable to consider a proper lift in order to generate POD modes for the lifted snapshots, satisfying 124 ESAIM: PROCEEDINGS AND SURVEYS homogeneous Dirichlet boundary conditions. This would lead to work with centered-trajectory method in the POD-ROM setting [20]. In the form it has been presented so far, POD seems to be only a bivariate data compression or reduction technique. Indeed, equation (3) simply says that the POD basis is the best possible approximation of order rof the given data set. In order to make POD a predictive tool, one couples the POD with the Galerkin procedure. This, in turn, yields a ROM, i.e., a dynamical system that represents the evolution in time of the Galerkin truncation. Thus, the Galerkin POD-ROM uses both Galerkin truncation and Galerkin projection. The former yields an approximation of the velocity field by a linear combination of the truncated POD basis: u(x, t)≈ur(x, t) = r X i=1 ai(t)ϕi(x),(7) where {ai(t)}r i=1 are the sought time-varying coefficients representing the POD-Galerkin trajectories. Note that r << N, where Ndenotes the number of degrees of freedom (d.o.f.) in a full order simulation (e.g., DNS). Replacing uwith urin (1), using the Galerkin method, and projecting the resulted equations onto the space Xr, one obtains the standard POD-ROM: d dt(ur, ϕr)+(b·∇ur, ϕr) + ε(∇ur,∇ϕr)+(gur, ϕr)=(f, ϕr)∀ϕr∈Xr.(8) Despite its appealing computational efficiency, the standard POD-ROM (8) has generally been limited to diffusion-dominated configurations. To overcome this restriction, we draw inspiration from the FE context, where stabilized formulations have been developed to deal with the numerical instabilities of the Galerkin method in advection-dominated configurations. 1.2. Streamline derivative projection-based method It is well known that a simple Galerkin truncation of POD basis leads to unstable results for advectiondominated configurations [5], and although the disregarded modes do not contain a significant amount of the system’s kinetic energy, they have a significant role in the dynamics of the reduced-order system. To model the effect of the discarded POD modes, various approaches have been proposed, both based on physical insights (cf., e.g., the survey in [37]), or on numerical stabilization techniques (cf. [10,12,20,26]). In this paper, we develop an approach that enters in the second group (no ad-hoc eddy viscosity is required, as it is in [37]), and aims to improve the previous works, because on one side a projection-stabilized structure is used (contrary to strategies in [10, 12,20]), which allows to act only on the high frequencies components of the advective derivative, and to control them, aspect of extreme importance. On the other side, a SD-based model is considered here, which is more adequate (with respect to a gradient-based model used in [26]) when dealing especially with advection-dominated configurations. This would allow to improve numerical stability and physical accuracy of the standard Galerkin POD-ROM also for convection-dominated and turbulent flows, with a rather simple driven structure, both for practical implementations such as to perform the numerical analysis. Let {Th}h>0be a regular family of triangulations of Ω. For any mesh cell K∈ Th, its diameter will be denoted by hKand h= maxK∈ThhK. To describe our strategy, we define the scalar product: (·,·)τ:L2(Ω) ×L2(Ω) →R,(v, w)τ=X K∈Th τK(v, w)K, and its associated norm: kvkτ= (v, v)1/2 τ, ESAIM: PROCEEDINGS AND SURVEYS 125 where for any K∈ Th,τKis a positive local stabilization parameter (to be determined later). Let us introduce the POD space: b XR= span {bϕ1,...,bϕR}, R ≤r, where bϕi,i= 1, . . . , R, are the POD modes associated to b K, defined as the snapshots correlation matrix with entries: b Kmn =1 N(b·∇u(·, tn+1),b·∇u(·, tm+1)) ,for m, n = 0, . . . , N −1.(9) Note that for classical POD modes asssociated to the standard correlation matrix Kmn, there already exists a theory on convergence rates and error bounds for POD expansions of parameterized solutions of heat equations, see e.g. [6–8]. With co-authors of the referred works, we aim to derive a similar analysis for POD modes associated to the advection correlation matrix b Kmn defined in (9). We consider the L2-orthogonal projection on b XR,PR:L2(Ω) −→ b XR, defined by: (u−PRu, bϕR) = 0,∀bϕR∈b XR.(10) Let P0 R=I−PR, where Iis the identity operator. We propose the Streamline Derivative projection-based POD-ROM (SD-POD-ROM) for (1): d dt(ur, ϕr)+(b·∇ur, ϕr) + ε(∇ur,∇ϕr)+(P0 R(b·∇ur), P0 R(b·∇ϕr))τ+ (gur, ϕr) = (f, ϕr)∀ϕr∈Xr.(11) Remark 1.2. When τK= 0 for any K∈ Th, the SD-POD-ROM (11) coincides with the standard POD-ROM, since no numerical dissipation is introduced. When R= 0, since numerical diffusion is extended to all the resolved modes {ϕ1, . . . , ϕr}, the SD-POD-ROM (11) becomes a penalty-stabilized method of the form: d dt(ur, ϕr)+(b·∇ur, ϕr) + ε(∇ur,∇ϕr)+(b·∇ur,b·∇ϕr)τ+ (gur, ϕr) = (f, ϕr)∀ϕr∈Xr,(12) which provides less accuracy with respect to the SD-POD-ROM (11), see remark 2.12 in section 2.2. Remark 1.3. Note that the new SD-POD-ROM (11) proposed in the present work rather differs from the V MS −POD −ROM introduced in [25]. Indeed, in [25], a gradient-based model for the standard POD-ROM is considered, which adds artificial viscosity by a term of the form: α(P0 R(∇ur), P0 R(∇ϕr)), αbeing a constant eddy viscosity coefficient, and P0 R=I−PR, with PRthe L2-orthogonal projection on the POD space defined by span{∇ϕ1,...,∇ϕR}. On the contrary, in the present work, we are adding an advection stabilization term, by basically following the streamlines, which seems to be a more natural choice when dealing especially with strongly advection-dominated configurations. This clearly differentiate the present work with respect to [25]. Also, the new SD-POD-ROM (11) is different from the SUPG-POD-ROM introduced in [20], since the former does not involve the full residual (only a streamline derivative stabilization term is introduced), thus presenting a simpler and cheaper structure for practical implementations such as to perform the numerical analysis, and also uses a projection-stabilized structure, which allows to act only on the high frequencies components of the advective derivative: this guarantees an extra-control on them that prevents high-frequency oscillations without polluting the large scale components of the approximation for advection-dominated problems, see remark 2.8 in section 2.1. 126 ESAIM: PROCEEDINGS AND SURVEYS 2. Error estimates In this section, we prove estimates for the average error: 1 N+ 1 N X n=0 kun−un rk, where the approximation unis the weak solution of (1) and un ris the solution of (11) discretized by FE, at time tn. We derive the error estimate in two steps. First, we gather some necessary assumptions and preliminary results in section 2.1. Then, we present the main result in section 2.2. 2.1. Technical background This section provides some technical results that are required for the numerical analysis. Throughout the paper, we shall denote by C,C1,C2,. . . positive constants that may vary from a line to another, but which are always independent of the FE mesh size h, the eigenvalues λi, and the diffusion coefficient ε. Of particular interest is the independence of the generic constant Cfrom ε, that is we will prove estimates that are uniform with respect to the diffusion coefficient, which is extremely relevant when advection-dominated problems are considered. We first make the following assumption, which is needed to prove the well-posedness of the weak formulation of (1). Hypothesis 2.1. (Coercivity and continuity) g−1 2∇·b≥α > 0; max{kgk∞,kbk∞}=β > 0.(13) For the FE discretization of the weak form of (1), we consider a family of finite dimensional subspaces Xh of X=H1 0(Ω) such that, for all v∈Hm+1(Ω) ∩X, the following assumption is satisfied. Hypothesis 2.2. (Approximability) inf vh∈Xh{kv−vhk+hk∇(v−vh)k} ≤ C hm+1kvkHm+1 ,1≤m≤k, (14) where kis the order of accuracy of Xh. For the POD approximation, the following POD inverse estimate was proven in [29], Lemma 2: Lemma 2.3. Let ϕi,i= 1, . . . , r, be POD modes, Mrbe the POD mass matrix with entries [Mr]ij = (ϕj, ϕi), and Srbe the stiffness matrix with entries [Sr]ij = (∇ϕj,∇ϕi),i, j = 1, . . . , r. Let k·k2denote the matrix 2-norm. Then, for all ϕr∈Xr, the following estimate holds: k∇ϕrk ≤ qkSrk2kM−1 rk2kϕrk.(15) Note that, since we have chosen H=L2in the POD method, kM−1 rk2= 1 in inequality (15). To prove optimal error estimates in time, we follow [29] and include the finite difference quotient ¯ ∂u(tn) = un−un−1 ∆t, for n= 1, . . . , N, in the set of snapshots χ={u(t0), . . . , u(tN),¯ ∂u(t1),...,¯ ∂u(tN)}. As pointed out in [29], the POD error formula (6) becomes: 1 2N+ 1 N X n=0     u(tn)− r X i=1 (u(tn), ϕi)ϕi     2 +1 2N+ 1 N X n=1     ¯ ∂u(tn)− r X i=1 ¯ ∂u(tn), ϕiϕi     2 = M X i=r+1 λi,(16) ESAIM: PROCEEDINGS AND SURVEYS 127 For the subsequent numerical analysis, we need the following technical hypothesis on the stabilization parameters τK: Hypothesis 2.4. The stabilization parameters τKsatisfy the following condition: 0< τK≤C hK,(17) for all K∈ Th, and a positive constant Cindependent of hand ε. Remark 2.5. The question whether the stabilization parameters should depend on the spatial resolution of the underlying FE space, or on the number of POD modes used has been addressed in [20], by means of numerical analysis arguments. In that work, numerical investigations using both definitions suggested that the one based on estimates from the underlying FE discretization provides a better suppression of numerical oscillations, and thus guarantees a more effective numerical stabilization. For this reason, we make here assumption 2.4 on the stabilization parameters, which is also essential for the subsequent numerical analysis. Lemma 2.6. Assume that Hypothesis 2.4 holds. Then, for all g∈L2(Ω), the following estimate is satisfied: kP0 R(g)kτ≤C√hkgk.(18) Proof. By using (17) and the stability of PRin the L2-norm, it follows: kP0 R(g)k2 τ≤C hkP0 R(g)k2≤C hkgk2. Thus, the estimate (18) can be deduced.  We introduce the bilinear form A(u, v)=(b·∇u, v) + ε(∇u, ∇v) + (P0 R(b·∇u), P0 R(b·∇v))τ+ (gu, v). The SD-POD-ROM (11) with a backward Euler time discretization reads: 1 ∆t(un+1 r−un r, ϕr) + A(un+1 r, ϕr)=(fn+1, ϕr)∀ϕr∈Xr.(19) We have the following existence and stability result for the fully discretized SD-POD-ROM (19): Lemma 2.7. Problem (19) admits a unique solution that satisfies the following bound: kuk rk2+ ∆t N−1 X n=0 εk∇un+1 rk2+kP0 R(b·∇un+1 r)k2 τ≤ ku0 rk2+∆t α N−1 X n=0 kfn+1k2,(20) for k= 0, . . . , N, with αthe coercivity constant in (13). Proof. Hypotheses 2.1 and 2.4 guarantee the well-posedness of (19). To prove estimate (20), we choose ϕr=un+1 rin (19), and decompose the bilinear form Ainto its symmetric and skew-symmetric parts: A= As+Ass, where As(u, v) = ε(∇u, ∇v)+(P0 R(b· ∇u), P0 R(b· ∇v))τ+(g−1 2∇·b)u, vand Ass(u, v) = b·∇u+1 2(∇·b)u, v, nothing that Ass(u, u) = 0. From (19), we obtain: (un+1 r−un r, un+1 r)+∆tAs(un+1 r, un+1 r)=∆t(fn+1, un+1 r).(21) Using the identity: (a−b)a=1 2(|a|2−|b|2+|a−b|2),∀a, b ∈R, 128 ESAIM: PROCEEDINGS AND SURVEYS and Young’s inequality, from (21) we get: kun+1 rk2−kun rk2+ε∆tk∇un+1 rk2+ ∆tkP0 R(b·∇un+1 r)k2 τ≤∆t αkfn+1k2.(22) Then, the stability estimate (20) follows by summing (22) from n= 0 to k≤N−1.  Remark 2.8. The stability estimate (20), which makes apparent the estimate of the advective stabilization term, guarantees an extra-control on the high frequencies of the advective derivative, which is not obtained by the standard Galerkin POD-ROM. This is an aspect of extreme importance, especially when dealing with advection-dominated configurations. In order to prove an estimate for kun−un rk, we will first consider the Ritz projection wr∈Xrof u∈X: A(u−wr, ϕr)=0,∀ϕr∈Xr.(23) The existence and uniqueness of wrfollow from the Lax–Milgram lemma. We now state an estimate for un−wn r, the error in the Ritz projection. Lemma 2.9. The Ritz projection wn rof unsatisfies the following error estimate: 1 N N X n=1 kun−wn rk ≤ C (1 + kSrk2)1/2 hm+1 1 N N X n=1 kunkHm+1 +v u u t M X i=r+1 λi  + (1 + ε)1/2 hm1 N N X n=1 kunkHm+1 +v u u t M X i=r+1 λi  .(24) The proof of this lemma can be directly derived by the one performed for the VMS-POD-ROM in [25]. Corollary 2.10. The Ritz projection wn rof unsatisfies the following error estimate up tp O(∆t2): 1 N N X n=1 k∂t(un−wn r)k ≤ C (1 + kSrk2)1/2 hm+1k∂tukL2(Hm+1)+v u u t M X i=r+1 λi  + (1 + ε)1/2 hmk∂tukL2(Hm+1)+v u u t M X i=r+1 λi  .(25) The proof of this corollary follows along the same lines as the proof of lemma 2.9. Note that it is exactly at this point that we use the fact that the finite difference quotients ¯ ∂u(tn) are included in the set of snapshots (see remark 1 in [29]). 2.2. Error estimate for the SD-POD-ROM We are now in position to prove the following error estimate result for the SD-POD-ROM defined by (19): ESAIM: PROCEEDINGS AND SURVEYS 129 Theorem 2.11. The solution of the SD-POD-ROM (19) satisfies the following error estimate: 1 N+ 1 N X n=0 kun−un rk ≤ C (1 + kSrk2)1/2 hm+1 1 N N X n=1 kunkHm+1 +k∂tukL2(Hm+1)+v u u t M X i=r+1 λi  + (1 + ε)1/2 hm1 N N X n=1 kunkHm+1 +k∂tukL2(Hm+1)+v u u t M X i=r+1 λi  +ku0−u0 rk+ ∆tk∂ttukL2(L2)+v u u th M X i=R+1 b λi ,(26) where the initial condition is given by the L2-orthogonal projection of u0on Xr: u0 r= r X i=1 (u0, ϕi)ϕi,(27) and b λi,i=R+1, . . . , M, in the right-hand side of (26) are the eigenvalues associated to the snapshots correlation matrix b Kpreviously defined in (9). Proof. We evaluate the weak form of (1) at t=tn+1, let the function test v=ϕr, and add and subtract the different quotient term un+1 −un ∆t: ∂tun+1 −un+1 −un ∆t, ϕr+un+1 −un ∆t, ϕr+a(un+1, ϕr)=(fn+1, ϕr),(28) where we have considered the bilinear form a(u, v)=(b·∇u, v) + ε(∇u, ∇v)+(gu, v). Subtracting (19) from (28), we obtain the error equation: ∂tun+1 −un+1 −un ∆t, ϕr+un+1 −un+1 r ∆t, ϕr−un−un r ∆t, ϕr+A(un+1−un+1 r, ϕr)+(a−A)(un+1, ϕr)=0. (29) We now decompose the error as un−un r= (un−wn r)−(un r−wn r) = ηn−φn r, so that by triangle inequality we have: kun−un rk ≤ kηnk+kφn rk.(30) Note that kηnkhas already been bounded in lemma 2.9. Thus, in order to estimate the error, we only need to estimate kφn rk. The error equation (29) can be written as: ∂tun+1 −un+1 −un ∆t, ϕr+ηn+1 −ηn ∆t, ϕr−φn+1 r−φn r ∆t, ϕr+A(ηn+1−φn+1 r, ϕr)+(a−A)(un+1, ϕr)=0. (31) We consider ϕr=φn+1 r, and note that, since φn+1 r∈Xr, then A(ηn+1, φn+1 r) = 0, so that we get: A(φn+1 r, φn+1 r) + 1 ∆tφn+1 r−φn r, φn+1 r=1 ∆tηn+1 −ηn, φn+1 r+rn, φn+1 r+ (a−A)(un+1, φn+1 r),(32) 136 ESAIM: PROCEEDINGS AND SURVEYS [25] T. Iliescu and Z. Wang,Variational multiscale proper orthogonal decomposition: convection-dominated convection-diffusion- reaction equations, Math. Comp., 82 (2013), pp. 1357–1378. [26] , Variational multiscale proper orthogonal decomposition: Navier-Stokes equations, Numer. Methods Partial Differential Equations, 30 (2014), pp. 641–663. [27] I. Kalashnikova and M. F. Barone,Efficient non-linear proper orthogonal decomposition/Galerkin reduced order models with stable penalty enforcement of boundary conditions, Internat. J. Numer. Methods Engrg., 90 (2012), pp. 1337–1362. [28] P. Knobloch and G. Lube,Local projection stabilization for advection-diffusion-reaction problems: one-level vs. two-level approach, Appl. Numer. Math., 59 (2009), pp. 2891–2907. [29] K. Kunisch and S. Volkwein,Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math., 90 (2001), pp. 117–148. [30] D. J. Lucia and P. S. Beran,Projection methods for reduced order models of compressible flows, J. Comput. Phys., 188 (2003), pp. 252–280. [31] Y. Maday, A. Manzoni, and A. Quarteroni,An online intrinsic stabilization strategy for the reduced basis approximation of parametrized advection-dominated problems, C. R. Math. Acad. Sci. Paris, 354 (2016), pp. 1188–1194. [32] P. Pacciarini and G. Rozza,Stabilized reduced basis method for parametrized advection-diffusion PDEs, Comput. Methods Appl. Mech. Engrg., 274 (2014), pp. 1–18. [33] J. R. Singler,New POD error expressions, error bounds, and asymptotic results for reduced order models of parabolic PDEs, SIAM J. Numer. Anal., 52 (2014), pp. 852–876. [34] L. Sirovich,Turbulence and the dynamics of coherent structures. I. Coherent structures, Quart. Appl. Math., 45 (1987), pp. 561–571. [35] G. Tadmor, O. Lehmann, B. R. Noack, L. Cordier, J. Delville, J.-P. Bonnet, and M. Morzy´ nski,Reduced-order models for closed-loop wake control, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 369 (2011), pp. 1513–1524. [36] S. Volkwein,Model reduction using proper orthogonal decomposition, tech. rep., University of Konstanz, Available at: http://www.math.uni-konstanz.de/numerik/personen/volkwein/teaching/POD-Vorlesung.pdf, 2011. [37] Z. Wang, I. Akhtar, J. Borggaard, and T. Iliescu,Proper orthogonal decomposition closure models for turbulent flows: a numerical comparison, Comput. Methods Appl. Mech. Engrg., 237/240 (2012), pp. 10–26.