A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations

Journal of Computational Physics 425 (2021) 109916 Contents lists available at ScienceDirect Journal of Computational Physics www.elsevier.com/locate/jcp A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations Mejdi Azaïez a, Tomás Chacón Rebollo b, Samuele Rubino b,∗ aInstitut Polytechnique de Bordeaux, Laboratoire I2M CNRS UMR5295, France bDepartamento EDAN & IMUS, Universidad de Sevilla, Spain a r t i c l e i n f o a b s t r a c t Article history: Received 6 May 2019 Received in revised form 3 June 2020 Accepted 7 September 2020 Available online 15 October 2020 Keywords: Finite element method Filtered advection stabilization A-posteriori stabilization Proper orthogonal decomposition Reduced order models Convection-dominated flows In this paper, we propose to improve the stabilized POD-ROM introduced in [48]to deal with the numerical simulation of advection-dominated advection-diffusion-reaction equations. In particular, we propose a three-stage stabilizing strategy that will be very useful when considering very low diffusion coefficients, i.e. in the strongly advectiondominated regime. This approach mainly consists in three ingredients: (1) the addition of a “streamline diffusion” stabilization term to the governing projected equations, (2) the modification of the correlation matrix defining the POD modes associated to the advection stabilization term, and (3) an a-posteriori stabilization scheme. Numerical studies are performed to discuss the accuracy and performance of the new method in handling strongly advection-dominated cases. ©2020 Elsevier Inc. All rights reserved. 1. Introduction Reduced Order Models (ROMs) applied to numerical design in modern engineering are a tool that is wide-spreading in the scientific community in the recent years in order to solve complex realistic multi-parameters, multi-physics and multiscale problems, where classical methods such as Finite Difference (FD), Finite Element (FE) or Finite Volume (FV) methods would require up to billions of unknowns. On the contrary, ROMs are based on a mathematically rigorous offline/online strategy, and the latter requires a reduced number of unknowns, which allows to face control, optimization, prediction and data analysis problems in almost real-time, that is, ultimately, a major goal for industrial applications. The reduced order modeling offline strategy relies on proper choices for data sampling and construction of the reduced basis (cf. [34]), which will be used then in the online phase, where a proper choice of the reduced model describing the dynamics of the system is needed. The key feature of ROMs is their capability to highly speedup computations, and thus drastically reduce the computational cost of numerical simulations, without compromising too much the physical accuracy of the solution from the engineering point of view. Among the most popular ROMs approaches, Proper Orthogonal Decomposition (POD) strategy provides optimal (from the energetic point of view) basis or modes to represent the dynamics from a given database (snapshots) obtained by a full-order system. Onto these reduced bases, a Galerkin projection of the governing equations can be employed to obtain a low-order dynamical system for the basis 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 *Corresponding author. E-mail addresses: [email protected] (M. Azaïez), [email protected] (T. Chacón Rebollo), [email protected] (S. Rubino). https://doi.org/10.1016/j.jcp.2020.109916 0021-9991/©2020 Elsevier Inc. All rights reserved. M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 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. This has led researchers to apply POD-ROMs to a variety of physical and engineering problems, including Computational Fluid Dynamics (CFD) problems in order to model advection-diffusion equations [25,26], see e.g. [27,28], and the Navier–Stokes Equations (NSE), see e.g. [11,12,17,29,42,46,57]. Once applied to the physical problem of interest, POD-ROMs can be used to solve engineering problems such as shape optimization [5,32] and flow control [6,15,31,55]. Although POD-ROMs can be very computationally efficient and relatively accurate in some flow configurations, they also present several drawbacks. For instance, for model reduction purposes, one only keeps few modes that are associated to the large eddies of the flow, which should be sufficient to give a good representation of the kinetic energy of the flow, due to the energetic optimality of the POD basis functions. However, the main amount of viscous dissipation takes place in the small eddies represented by basis functions that are not taken into account, and thus the leading reduced order system is not able to dissipate enough energy. So, although the disregarded modes do not contain a significant amount of kinetic energy, they have a significant role in the dynamics of the reduced order system. It is then necessary to close the POD-ROMs by modeling the interaction between the computed and the unresolved modes. This problem establishes a parallelism to Large Eddy Simulations (LES) [50]of turbulent flows, where the effect of the smallest flow structures on the largest ones is modeled. Since these are also in non-linear interactions, a proper non-linear efficient and accurate closure model should be proposed also in the POD context, considering that in this context the concepts of energy cascade and locality of energy transfer are still valid [23]. To prevent the loss of accuracy of POD-ROMs due to cutting out the POD modes corresponding to the viscous scales, various approaches have been proposed, both based on physical insights (cf., e.g., the survey in [57]), and numerical stabilization techniques (cf. [11–14,30,38]). We emphasize, however, that reduced order closure modeling and stabilization are two related, yet different issues. For example, if one considers a linear problem like the advection-diffusion-reaction problem investigated in the present manuscript, one could choose the solution norm and construct POD basis functions that are orthogonal in that norm. In that case, there would be no reduced order closure problem (in the corresponding inner product). Nevertheless, in the advection-dominated regime, the reduced order numerical stabilization would still be relevant, especially for low diffusion coefficients. On the other hand, the main goal of reduced order closure modeling is to increase the accuracy of ROMs, having some effect on their numerical stability too. Indeed, in order to increase the ROMs accuracy, reduced order closure models usually add numerical dissipation. This numerical dissipation aims at increasing the physical accuracy (i.e., matching the Kolmogorov energy cascade), and also allows to address numerical instabilities due to the truncation in Galerkin models [47]. To address this issue, in [48]a Streamline Derivative projection-based strategy for the numerical stabilization of POD- ROMs (SD-POD-ROM) has been introduced. The proposed model has been numerically analyzed for advection-diffusion- reaction equations, by mainly deriving the corresponding error estimates. Some preliminary numerical tests have been performed in [48]for a moderate Péclet number, showing the efficiency of the proposed method, as well as the increased accuracy over the standard POD-ROM that discovers its well-known limitations very soon in the numerical settings considered, i.e. for moderately low diffusion coefficients. In this paper, we aim to improve this approach by proposing a three-stage stabilizing strategy that will be very useful when considering very low diffusion coefficients, i.e. in the strongly advection-dominated regime. This approach mainly consists in three ingredients: (1) the addition of a “streamline diffusion” stabilization term to the governing projected equations, (2) the modification of the correlation matrix defining the POD modes associated to the advection stabilization term, and (3) an a-posteriori stabilization scheme. Parallel and independently to the current paper, a SUPG-POD-ROM combined with isogeometric analysis has been very recently proposed and analyzed in [45]to address, similarly to the present study, advection-dominated advection-diffusion-reaction problems. The latter and the present study independently perform a numerical investigation of two different stabilization POD-ROMs to address advection-dominance in POD solution to advection-diffusion-reaction equations. The rest of the paper is organized as follows: in section 2, we briefly describe the POD methodology and introduce the SD-POD-ROM for advection-diffusion-reaction problems. In section 3, we describe the process of a-posteriori stabilization in a general framework and how to apply it to the considered problems. Numerical studies are performed in section 4to discuss the accuracy and efficiency of our method in handling strongly advection-dominated cases, and also its robustness for long time integrations on periodic systems. Finally, section 5presents the main conclusions of this work and future research directions. 2. Streamline derivative projection-based POD-ROM In this paper, the proposed stabilization is preliminary analyzed and tested for the POD-ROM numerical approximation of advection-dominated advection-diffusion-reaction problems of the form: ⎧ ⎨ ⎩ ∂tu+b·∇u−νu+gu =fin ×(0,T), u=0on×(0,T), u(x,0)=u0(x)in , (2.1) 2 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 where bis the given advective field, ν<< 1the diffusion parameter, gthe reaction coefficient, fthe forcing term, the computational domain in Rd, d =2or 3, t∈[0, T], with Tthe final time, and u0the initial condition. For the sake of simplicity, we have imposed homogeneous Dirichlet boundary conditions on the whole boundary  =∂. To define the weak formulation of problem (2.1), let us consider the space: X=H1 0=v∈H1():v=0on, where H1is the usual Sobolev space [16]. We shall consider the following variational formulation of (2.1): Find u:(0, T) −→ Xsuch that d dt (u,v)+(b·∇u,v)+ν(∇u,∇v)+g(u,v)=(f,v)∀v∈X,(2.2) where (·, ·)stands for the L2-inner product in . In order to give a FE approximation of (2.2), let {Th}h>0be a family of affine-equivalent, conforming (i.e., without hanging nodes) and regular triangulations of , formed by triangles or quadrilaterals (d =2), tetrahedra or hexahedra (d =3). For any mesh cell K∈Th, its diameter will be denoted by hKand h =maxK∈ThhK. We consider Xh⊂Xa suitable FE space. The FE approximation of (2.2)can be written as follows: Find uh∈Xhsuch that d dt (uh,vh)+(b·∇uh,vh)+ν(∇uh,∇vh)+g(uh,vh)=(f,vh)∀vh∈Xh.(2.3) It is well-known that, in the case of low diffusion coefficient ν<< 1, the standard Galerkin method (2.3)is generally unstable and leads to globally polluted solutions presenting strong spurious oscillations. In this paper, we thus propose to first consider an offline stabilization procedure, which becomes necessary to deal with the numerical instabilities of the Galerkin method and to generate the snapshots for the online phase with a reasonable accuracy. In particular, we consider a simplification of the Streamline Derivative-based (SD-based) approach used by Knobloch and Lube (see [43]) in the FE context, which only acts on the high frequencies of the advective derivative. This approach consists in adding a filtered advection stabilization term by basically following the streamlines to prevent spurious instabilities due to dominant advection, but using a simple interpolation operator on a continuous buffer FE space instead of a local projection operator on a discontinuous enriched FE space (see [1]for more details). This stabilization term acts on the high frequency 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. This method falls into the class of Local Projection Stabilization (LPS) methods (cf. [2,4]). To briefly recall this approach, assume that the discrete space Xhis formed by piecewise polynomial functions of degree m ≥2, e.g. Xh=Pm∩X, where Pmdenotes the space of continuous functions whose restriction to each mesh cell K∈This the Lagrange polynomial of degree less than or equal to m. We define the scalar product: (·,·)τ:L2()×L2()→R,(v,w)τ= K∈Th τK(v,w)K, and its associated norm: vτ=(v,v)1/2 τ, where for any K∈Th, τKis in general a positive local stabilization parameter. The working expression for τKused in this context, designed by asymptotic scaling arguments, is: τK=c1 ν h2 K+c2b∞ hK+c3g−1 , where c1, c2and c3are positive algorithmic constants (see [48]for more details). The LPS method by interpolation applied to advection-diffusion-reaction equations is stated by: Find uh∈Xhsuch that ⎧ ⎪ ⎨ ⎪ ⎩ d dt (uh,vh)+(b·∇uh,vh)+(π h(b·∇uh), π h(b·∇vh))τ +ν(∇uh,∇vh)+g(uh,vh)=(f,vh)∀vh∈Xh, (2.4) 3 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 where π h=Id −πhis the “fluctuation operator”, with Id the identity operator and πha locally stable interpolation operator from L2()onto a projection space Dhdefined on the same mesh Thand formed by continuous FE (e.g., Dh=Pm−1), satisfying optimal error estimates (cf. [21]). In practical implementations, we choose πhas a Scott–Zhang-like [51]linear interpolation operator in the space P1(since we consider P2as FE solution space), implemented in the software FreeFem++ [33]. This interpolant may be defined as: ∀x∈,πh(v)(x)= a∈N Ih(v)(a)ψa(x), where Nis the set of Lagrange interpolation nodes of P1, ψaare the Lagrange basis functions associated to N, and Ih is the interpolation operator by local averaging of Scott–Zhang kind, which coincides with the standard nodal Lagrange interpolant when acting on continuous functions (cf. [21], section 4). 2.1. Proper orthogonal decomposition reduced order model For the report to be self-contained, this section briefly presents the computation of a basis for ROMs with POD. For more details, the reader is referred to [22,35,52,53,56]. 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 the best approximation: minu(x,t)− r  i=1 (u(x,t), ϕi)Hϕi 2 L2(0,T;H) ,(2.5) in a real Hilbert space H. Although Hcan be any real Hilbert space, in what follows we consider H=L2(), with induced norm · =(·, ·)1/2=⎛ ⎝  |·|2⎞ ⎠ 1/2 .Thus, the L2(0, T; L2()) norm is used, since it is directly related to the kinetic energy of the flow field. In the framework of the numerical solution of Partial Differential Equations (PDEs), 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 considered, 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 [44]: min 1 N+1 N  n=0u(·,tn)− r  i=1 (u(·,tn), ϕi)ϕi 2 ,(2.6) subject to the condition ϕj,ϕi=δij, 1 ≤i, j ≤r, where δij is the Kronecker delta. To solve the optimization problem (2.6), one can consider the eigenvalue problem: Kzi=λizi,for 1,...,r,(2.7) 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,(2.8) 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 (2.6), i.e. the POD basis, is given by: ϕi(·)=1 √λi 1 √N+1 N  n=0 (zi)nu(·,tn), 1≤i≤r,(2.9) where (zi)nis the n-th component of the eigenvector zi. It can also be shown that the following POD error formula holds [35,44]: 1 N+1 N  n=0u(·,tn)− r  i=1 (u(·,tn), ϕi)ϕi 2 = M  i=r+1 λi,(2.10) where Mis the rank of χ. 4 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 We consider the following space for the POD setting: Xr=span{ϕ1,...,ϕr}. Remark 2.1. Since, as shown in (2.9), the POD modes are linear combinations of the snapshots, the POD modes satisfy the boundary conditions in (2.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 homogeneous Dirichlet boundary conditions. This would lead to work with centered-trajectory method in the POD-ROMs setting [30]. One can also implement boundary conditions in ROMs constructed using continuous projection weakly, see [41]for more details on this issue. In the form it has been presented so far, POD seems to be a bivariate data compression or reduction technique, see e.g. [10]. Indeed, equation (2.6)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 reduced order system, 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 solution by a linear combination of the truncated POD basis: u(x,t)≈ur(x,t)= r  i=1 ai(t)ϕi(x), (2.11) where {ai(t)}r i=1are the sought time-varying coefficients representing the POD-Galerkin trajectories. Note that r<< Ndof , where Ndof denotes the number of degrees of freedom (d.o.f.) in a full order simulation (e.g., DNS). Replacing uwith ur in (2.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.(2.12) Despite its appealing computational efficiency, the standard POD-ROM (2.12) has generally been limited to diffusiondominated configurations. To overcome this restriction, we draw inspiration from the FE context, where stabilized formulations, such as (2.4)for instance, have been developed to deal with the numerical instabilities of the Galerkin method in advection-dominated configurations. 2.2. Streamline derivative projection-based method For ease of reading, we recall hereafter the approach leading to the SD-POD-ROM originally introduced and numerically analyzed in [48]. Let us introduce the POD space:  Xr=span{ ϕ1,..., ϕr}, where  ϕi, i =1, ..., r, are the POD modes associated to  K, defined as the snapshots correlation matrix with entries:  Kmn =1 N+1(b·∇u(·,tn), b·∇u(·,tm)),for m,n=0,...,N.(2.13) Note that for classical POD modes associated 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. [7–9]. With co-authors of the referred works, following the guidelines given there, we aim to derive a similar analysis for POD modes associated to the advection correlation matrix  Kmn defined in (2.13). We consider the L2-orthogonal projection on  Xr, Pr:L2() −→  Xr, defined by: (u−Pru, ϕr)=0,∀ ϕr∈ Xr.(2.14) Let P r=Id −Pr. We propose the Streamline Derivative projection-based POD-ROM (SD-POD-ROM) for (2.1): ⎧ ⎪ ⎨ ⎪ ⎩ d dt (ur,ϕr)+(b·∇ur,ϕr)+(P r(b·∇ur), P r(b·∇ϕr))τ +ν(∇ur,∇ϕr)+(gur,ϕr)=(f,ϕr)∀ϕr∈Xr. (2.15) We introduce the bilinear form A(u, v) =(b·∇u, v) +(P r(b·∇u), P r(b·∇v))τ+ν(∇u, ∇v) +(gu, v). The SD-POD-ROM (2.15)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.(2.16) 5 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Remark 2.2. In [48], we have proved that the solution of the fully discretized SD-POD-ROM (2.16)is stable and converges to the solution of the continuous problem (2.2). In particular, we have proved error estimates that are uniform with respect to the diffusion coefficient, which is extremely relevant when advection-dominated problems are considered, like in this work. Remark 2.3. When τK=0for any K∈Th, the SD-POD-ROM (2.15) coincides with the standard POD-ROM (2.12), since no numerical dissipation is introduced. Also, note that in this paper we directly consider the projection over the same number rof POD modes retained for the ROMs solution. Indeed, due to the slow convergence of the POD eigenvalues associated to the advection correlation matrix  Kmn in case of very low diffusion (see section 4) and the fact that error estimates for the SD-POD-ROM are directly proportional to them (cf. [48], Theorem 2.11), this improves results obtained by projecting over a number R <r, as initially proposed in [48]. Remark 2.4. Note that the SD-POD-ROM (2.15) rather differs from the VMS-POD-ROM introduced in [36]. Indeed, in [36], a gradient-based model for the standard POD-ROM is considered, which adds artificial viscosity by a term of the form: α(PR(∇ur), PR(∇ϕr)), being αa constant eddy viscosity coefficient, and PR=Id −PR, with PRthe L2-orthogonal projection on the POD space defined by span{∇ϕ1, ..., ∇ϕR}, R <r, making it applicable just to H1-POD basis (here, L2-POD basis is used), for which the decay of POD eigenvalues is rather slow in presence of strongly advection-dominated configurations (similar to the decay of POD eigenvalues associated to the advection correlation matrix (2.13), see, e.g., Figs. 12, 18), and this leads to higher POD errors [37]. On the contrary, in the present work, we are adding an advection stabilization term, by just following the streamlines, which seems to be a more natural choice when dealing especially with strongly advection-dominated regimes. We emphasize that the POD modes for the advection correlation matrix (2.13)are only used to construct the advection stabilization term through (2.14). This clearly differentiates the present work with respect to [36]. Also, the SD-POD-ROM (2.15)is different from the SUPG-POD-ROM introduced in [30], 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 frequency 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 (cf. [48], Lemma 2.7). We emphasize, however, that the SD-POD-ROM (2.15)is not fully consistent, but verifies optimal error estimates (cf. [48]). Instead, the SUPG-POD-ROM introduced in [30] retains numerical consistency, in the sense that the continuous solution exactly satisfies the discrete equations, whenever it is smooth enough. In terms of computational cost, the offline phase of the SUPG-POD-ROM is more expensive than the one of the SD-POD-ROM, since the former is fully residual-based, while the online phase is almost comparable. In terms of accuracy, both methods give similar reliable results (see section 4), especially when combined with a-posteriori stabilization described herein. 3. A-posteriori stabilization To describe the process of a-posteriori stabilization in a general framework, let us consider an elliptic variational problem: Find x∈Xsuch that b(x,w)=l(w)=f,w,∀w∈X,(3.1) where Xis a Hilbert space. The form bis defined on X×Xand l ∈X, being Xthe topological dual of X. Consider a family of sub-spaces of finite dimension of X, {Xi}i∈I, for some set of indices I. Let us assume that we solve problem (3.1)by the Galerkin method on Xi: Find xi∈Xisuch that b(xi,wi)=l(wi), ∀wi∈Xi.(3.2) Assume that the space Xiis decomposed into Xi=Yi⊕Zi, where Yiand Ziare subspaces of Xi. Let xi=yi+zibe the unique decomposition that xiadmits with yi∈Yiand zi∈Zi. Problem (3.2)may be recast as a variational problem for the only unknown yi, as follows. Denote by Athe operator from Xon Xdefined by the form b; that is for v ∈X, Avis the element of Xdefined by: Av,w=b(v,w), ∀w∈X. Denote by Ri:X→ Zithe “static condensation”operator on Zigenerated by the form b, defined for ϕ∈Xby: b(Ri(ϕ), wi)=ϕ,wi,∀wi∈Zi. Let us introduce the “condensed” variational formulation to problem (3.2). To do so, we consider the operators bcand lcas: bc(y,v)=b(y,v)−b(Ri(A∗v), Ri(Ay)), lc(v)=l(v)−b(Ri(A∗v), Ri(f)), ∀y,v∈X, 6 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 where A∗denotes the adjoint of the operator A. The “condensed” variational formulation to problem (3.2) reads: Find yi∈Xisuch that bc(yi,vi)=lc(vi), ∀vi∈Yi.(3.3) We next introduce the following definition: Definition 3.1. The family of finite-dimensional spaces {(Yi, Zi)}i∈I, where Iis a set if indices, is called to satisfy the saturation property if there exists a constant α>0such that yiX+ziX≤αxi+yiX,∀yi∈Yi,zi∈Zi,∀i∈I. The saturation property can be viewed as an inverse triangular inequality. It can be readily proved that this property is equivalent to the existence of some constant β>0such that |(yi,zi)X|≤(1−β)yiXziX,∀yi∈Yi,zi∈Zi;(3.4) actually we may take β=2 α2. Then, we can interpret the saturation property in the sense that the angle between spaces Yiand Zi, defined by arccossup yi∈Yi\{0},zi∈Zi\{0} (yi,zi)X yiXziX is uniformly bounded from below by a positive angle, with respect to i ∈I. Remark 3.2. Note that the argument of saturation property, applied here for the first time, up to our knowledge, to PODROM approximations to propose a cure for instabilities due to advection-dominance in POD solution to advection-diffusion- reaction equations, gave also a mathematical argument to perform the numerical analysis of recently proposed stabilization POD-ROMs [24,49] that take into account the pressure instability for incompressible flows governed by the NSE. Then, it holds (cf. [20]): Theorem 3.3. Assume that the spaces Yiand Zisatisfy Yi∩Zi=∅. Then: 1. Let xi=yi+zibe the unique decomposition that xiadmits with yi∈Yiand zi∈Zi. Then, xiis the solution of the Galerkin method (3.2)if and only if yiis the solution of the “condensed” variational formulation (3.3), and zi=Ri(l −A(yi)). 2. Assume, in addition, that the family of pairs of spaces {(Yi, Zi)}i∈Isatisfies the saturation property. Then, there exists a constant C>0such that yiX+ziX≤ClX,ciX≤ClX,(3.5) where ci=Ri(A(yi)). We may take advantage of this result to set up an a-posteriori stabilization procedure for the Galerkin solution of steady advection-reaction-diffusion equation. In this case, the framework Hilbert space is X=H1 0(). Assume that the space Yi contains in some sense the large scales (or low frequency) component of the space Xi. For instance, if Xiis a FE space constructed on a grid of a given diameter, Yicould be a FE subspace of Xiconstructed on a grid with a larger diameter, or with polynomials of lower degree. Also, if Xiis a POD space, then Yicould be a subspace formed by a truncated set of basis functions of low frequency. In both cases, Ziwill be a space containing the small scales (or high frequency) components of the space Xi. In this framework, ciis a representation on Zi(by means of the static condensation operator) of the small-scale components of the advection-diffusion-reaction operator Aacting on the large-scale component yiof the solution xi. Due to the second estimate in (3.5), ciis uniformly bounded in Xnorm. We interpret this bound as an a-posteriori stabilization effect. The stabilization effect largely depends on the actual choice for spaces Yiand Zi. For instance, for one-dimensional steady advection-diffusion equations with constant advection velocity, diffusion and forcing term, this choice may be made optimal when Xiis formed by piecewise affine finite elements, as follows. Assume that the space Xiis built on a grid of grid size h, Th. The subspace Yiis formed by piecewise affine finite elements on a grid with double grid size 2h, T2h. Then, there is a unique subspace Zisuch that the solution yiof the condensed variational formulation (3.3) coincides with the exact solution xof problem (3.1)at the nodes of the grid T2h. For some other choices of Zithere could be, however, an over-diffusive effect that yields a large damping of yi(cf. [20]). Note that to compute yifrom xiit is not necessary to build the space Yi. Indeed, it suffices to construct a projection operator i:xi∈Xi→ yi∈Yi. To each actual setting for ithere corresponds a space Zi, as Zi=(Id −i)(Xi). For Lagrange 7 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 FE spaces, in practice the simplest way to compute yiis to retain just the degrees of freedom of xithat correspond to the coarser grid on which Yiis built. Denote by {a1, a2, ···, ap}the Lagrange interpolation nodes of Yi, and by {ϕ1, ϕ2, ···, ϕp} the associated Lagrange basis functions of Yi. There exist a complementary set of interpolation nodes {ap+1, ap+2, ···, ar} and associated basis functions {ϕp+1, ϕp+2, ···, ϕr}such that {ϕ1, ϕ2, ···, ϕr}is a basis of Xi. Then, the operator iis defined, for any xi= r  k=1 αkϕk∈Xias: ir  k=1 αkϕk= p  k=1 αkϕk∈Yi.(3.6) The sub scale space Zifor this procedure is generated by the complementary basis functions {ϕp+1, ϕp+2, ···, ϕr}. In [20], it is proved that the pairs of spaces {(Yi, Zi)}i∈Iconstructed in this way indeed satisfy the saturation property. In this case the index imay be identified, as usual, with the diameter of the triangulation h. For POD approximations, the procedure is quite similar. The space Xiis generated by the basis functions {ϕ1, ϕ2, ···, ϕr}, then the operator iis defined by truncation of the POD series xi= r  k=1 αkϕk∈Xiright by (3.6), and again the spaces Yi and Ziare respectively spanned by {ϕ1, ϕ2, ···, ϕp}and {ϕp+1, ϕp+2, ···, ϕr}. In this case, the index imay be identified with the dimension rof the space Xi. In this paper we will apply the a-posteriori stabilization procedure in the offline stage, in which Xiis a FE space, and also in the online stage, in which Xiis a POD space. Hereafter, we detail the post-processing algorithm for the online stage: •Post-processing algorithm (online stage). (i) For n =0, 1, ..., N−1, given un r∈Xr, find un+1 rsuch that (2.16) holds. (ii) Represent the solution un+1 rusing R <rmodes. Although Rcould be estimated on the fly, minimizing for instance the error in a certain norm with respect to the snapshot solution (when available) at each ROM time step, for the considered numerical experiments choosing to truncate at R = r−10 gave the best balance between accuracy and suppression of spurious oscillations. 4. Numerical studies In this section, we present some numerical experiments to mainly assess accuracy and performance of the combination of the Streamline Derivative projection-based stabilization technique (2.15)with online stabilizing post-processing strategy. We consider the numerical computation of POD-ROM solutions to strongly advection-dominated advection-diffusion-reaction equations. As mentioned above, while for the Full Order Model (FOM) this strategy consists in interpolating the FOM solution on a coarser mesh (in practice, T2h), for the ROM the a-posteriori stabilization consists in truncating the ROM solution once obtained, as detailed in the algorithm above. This leads to a computationally efficient and mathematically founded offline/online algorithm (completely separated), implemented over the standard POD-Galerkin ROM. Actually, two applications (offline and online) of the stabilized post-processing technique are studied in this paper, where we will show the good performances of this technique to stabilize highly oscillatory FOMs and ROMs numerical solutions of strongly advectiondominated problems. From the following numerical results, we can observe that separately the two numerical stabilization strategies proposed (SD-POD-ROM and a-posteriori stabilization) already provide an improvement in general over the standard POD-Galerkin ROM. However, a further improvement is reached when we combine the two stabilization methods, which allows to obtain almost the same accuracy of more complex fully residual-based stabilization methods, such as SUPG-POD-ROM. The first numerical test 4.1 concerns an almost pure transient transport problem with a rotating cylinder. The second numerical test 4.2 concerns a 2D traveling wave displaying a sharp internal layer moving in time. In both cases, we employ P2(piecewise quadratic) FE on relatively coarse uniform spatial discretizations, and the backward Euler method for temporal discretization with time step t=10−3. In particular, FE meshes are significantly coarser than the width of the internal layers, which is common in practice. POD modes are represented using P2shape functions in order to perform the projection step of the ROM procedure, similar to the elements we use for the FOM discretization we compare to. The open-source FE software FreeFem++ [33] has been used to run all numerical experiments. In terms of computational cost, the CPU time of the tested ROMs is at least three orders of magnitude lower than the CPU time of the corresponding FOMs. Also, note that at online level adding the proposed stabilization techniques results in a not significant increase of the CPU time with respect to the standard POD-Galerkin ROM, thus proving the computational efficiency of the different ROMs stabilization strategies employed. 8 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 1. Example 4.1: Initial condition. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) 4.1. 2D rotating cylinder In this section, an almost pure transient transport problem with a rotating body will be considered. In particular, this problem is given in the unit disc  ={(x, y) ∈R2:x2+y2<1}by the advection-diffusion-reaction equation (2.1)with advection field b=(−y, x)T, reaction coefficient g=0, forcing term f=0, and a very small value for the diffusion parameter ν=10−20, as in [3]. The initial condition u0is given by: u0=0.5tanhe−10[(x−0.3)2+(y−0.3)2−0.5] 10−3+1,(4.1) which consists in a cylinder of height 1centered at (0.3, 0.3), as shown in Fig. 1. This condition is smooth, but has a sharp layer with thickness of order 10−3. The mesh is uniform with 256 triangles along the boundary of , which leads to mesh size h =4.26 ·10−2, thus the layer is under-resolved. The rotation is counter-clockwise and the solution after complete revolutions should be essentially the same as the initial condition, since the diffusion parameter ν=10−20 is very small. A pure transient transport problem with this data was considered in [18]. This example leads to astrongly advection-dominated problem, and therefore an offline stabilization procedure becomes necessary to deal with the numerical instabilities of the Galerkin method. As announced in section 2, in this work we preliminarily consider the LPS-FE by interpolation Method (LPS-FEM) given by (2.4), to which we further apply the a-posteriori stabilization described in section 3. 4.1.1. Short time behavior In first instance, we just compute one complete revolution of the cylinder being transported around the unit disc, i.e. the computational time interval is [0, T] =[0, 2π], and test the SD-POD-ROM in this interval where the snapshots are computed. Thus, we are evaluating the SD-POD-ROM in the reproductive (in time) regime. Note that the application of the a-posteriori stabilization described in the previous section further improves the accuracy provided by the LPS-FEM, as shown in Fig. 2, where we consider: varh(t)=max (x,y)∈ uh(x,y,t)−min (x,y)∈uh(x,y,t), as measure for under- and overshoots, as in [40]. Indeed, we observe that, even if both methods give similar error levels, LPS-FEM with post-processing is superior to LPS-FEM, for which the quantity varh(t)shows much larger oscillations. Note that the optimal value of varh(t)equals to 1for all t. As for the online phase, we perform a comparison between the SD-POD-ROM (2.15) and the SUPG-POD-ROM [30]by considering the application or not of the a-posteriori stabilization technique mentioned above, adapted to the POD-ROMs framework. The POD modes are generated in L2by the method of snapshots by storing every tenth FOM solution in the computational time interval [0, T] =[0, 2π], so that 629 snapshots were used. POD basis were constructed by using LPSFEM with stabilizing post-processing, to limit the influence of POD noisy data in the online phase. In Fig. 3, we show the decay of POD eigenvalues associated both to the snapshots correlation matrix (2.8) and the advection correlation matrix (2.13)in this case. To check the temporal behavior of the online spurious oscillations, we compute: varr(t)=max (x,y)∈ ur(x,y,t)−min (x,y)∈ur(x,y,t), 9 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 10. Example 4.1.2: Dominant POD modes for the correlation matrix (2.8) (left) and the advection correlation matrix (2.13)(right). 16 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 11. Example 4.2.1: Final solution profiles along the mean diagonal for different FOMs. Table 2 Example 4.2.1: L2-norm of the deviation from the final exact solution profile along the mean diagonal for different FOMs. Offline methods eFOM 0,ν=10−6 DNS-FEM 0.1828 DNS-FEM post-processing 0.1257 LPS-FEM 0.0576 LPS-FEM post-processing 0.0618 4.2.1. Case ν=10−6 In this case, we consider a uniform triangular mesh with mesh size h =1.41 ·10−2, which is relatively coarse with respect to the width of the internal layer. First, we tested different FOM: the Direct Numerical Simulation (2.3)(DNS-FEM), where no stabilization is introduced, a DNS with stabilized post-processing (DNS-FEM post-processing), the LPS (by interpolation)- FEM (2.4)(LPS-FEM), and the LPS-FEM with stabilized post-processing (LPS-FEM post-processing). In Fig. 11, we show for the different methods the final solution profiles along the mean diagonal of the computational domain compared with the corresponding exact solution profile. From this figure, it is evident that a DNS (i.e., no stabilization) gives oscillatory results, which are only in part corrected by applying the a-posteriori stabilization. Thus, since the problem is advection-dominated and the solution has already a steep internal layer, the use of a stabilized discretization is necessary when using relatively coarse meshes. For this purpose, we considered LPS by interpolation method, for which oscillations are rather reduced, and application or not of the a-posteriori stabilization almost gives similar results. A quantitative comparison between the different FOMs is given in Table 2, where the deviation eFOM 0from the final exact solution profile along the mean diagonal in a normalized discrete L2-norm subject to (4.5)is displayed. We may observe that, while for DNS methods errors are greater than 10%, for LPS-FEM methods are comparable and below 10%, being slightly better for the LPS-FEM method without a-posteriori stabilization. So, for this case, POD basis were constructed by using LPS-FEM method (2.4), and the studied ROMs thus used just slightly noisy POD data, which is unavoidable for strongly advection-dominated problems on realistic grids. In Fig. 12, we show the decay of POD eigenvalues associated both to the snapshots correlation matrix (2.8) and the advection correlation matrix (2.13). One can observe that the decay of the POD eigenvalues associated to the advection correlation matrix is rather slow, due to the low diffusion. Fig. 13 displays the dominant (i.e., most energetic) first five POD modes for the snapshots correlation matrix (2.8)(left) and the advection correlation matrix (2.13) (right). One can observe that the dominant POD modes for the advection correlation matrix (2.13) appear more oscillatory than the ones for the snapshots correlation matrix (2.8). However, adding the corresponding stabilization term in the online phase greatly improves the results over the standard POD-ROM, since allows to control the high frequency components of the advective derivative, main responsible for numerical oscillations. Fig. 14 presents results for all considered ROMs: the standard POD-Galerkin ROM (2.12)(G-ROM), the G-ROM with online stabilized post-processing (G-ROM post-processing), the SD-POD-ROM (2.15)(SD-ROM), and the SD-ROM with online stabilized post-processing (SD-ROM post-processing). In particular, we show for the different methods the final solution 17 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 12. Example 4.2.1:PODeigenvalues. Table 3 Example 4.2.1: Captured system’s kinetic energy and L2-norm of the deviation from the final exact solution profile along the mean diagonal for different ROMs at r=30, 60, 90. ν=10−6r=30 r=60 r=90 Captured system’s Ekin(%)99.76 99.98 >99.99 ν=10−6eROM 0 Online methods r=30 r=60 r=90 G-ROM 0.3743 0.1567 0.1067 G-ROM post-processing 0.3180 0.1389 0.0605 SD-ROM 0.3465 0.1435 0.0637 SD-ROM post-processing 0.2671 0.1383 0.0579 profiles along the mean diagonal of the computational domain compared with the corresponding exact solution profile, at r=30, 60, 90 (from top to bottom). One can observe that applying the online a-posteriori stabilization greatly improves results for the standard Galerkin-ROM (totally oscillatory), making it comparable with the stabilized SD-ROM, for which applying or not the online a-posteriori stabilization almost gives similar results. This is reflected by results depicted in Table 3, where the deviation eROM 0from the final exact solution profile along the mean diagonal in a normalized discrete L2-norm subject to (4.4)is displayed. One can see that, for r=90, SD-ROM post-processing method almost reaches the same accuracy of the offline phase by almost suppressing the influence of noisy modes. Also, note that although the first r=30 POD modes already capture more than 99% of the system’s kinetic energy, all ROMs yield poor quality results for which the peak of the front is not reached, and they display visible numerical oscillations, reflecting the complexity of the problem. Augmenting the number of POD modes allows to reach the peak of the front for all methods. However, whereas the solution of the G-ROM remains globally polluted with spurious oscillations, the application to it of the online a-posteriori stabilization already reduces to few oscillations and localizes them mainly near the steep layer, allowing to compute a rather accurate solution in this case, comparable with the one of the stabilized SD-ROM and of the offline phase. In Fig. 15, we show the numerical solution at T=1for the best performing SD-ROM with online a-posteriori stabilization for r=30, 60, 90 (from top to bottom). With this method, numerical unphysical oscillations are practically eliminated by gradually increasing the number of POD modes. 4.2.2. Case ν=10−8 In this case, we consider a uniform triangular mesh with mesh size h =9.43 ·10−3. Thus, a finer grid with respect to the previous case is used, which is necessary to maintain numerical diffusion within reasonable limits. Nevertheless, it remains relatively coarse with respect to the width of the internal layer. Again, we tested different FOMs: DNS-FEM, DNS-FEM postprocessing, LPS-FEM, and LPS-FEM post-processing. In Fig. 16, we show for the different methods the final solution profiles along the mean diagonal of the computational domain compared with the corresponding exact solution profile. Offline results proved again the necessity to consider LPS method to avoid globally spurious oscillations, but also that the application of the a-posteriori stabilization greatly improves the results of the LPS-FEM in this case. Indeed, error levels 18 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 13. Example 4.2.1: Dominant POD modes for the correlation matrix (2.8) (left) and the advection correlation matrix (2.13)(right). 19 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 14. Example 4.2.1: Final solution profiles along the mean diagonal for different ROMs at r=30,60,90 (from top to bottom). 20 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 15. Example 4.2.1: Numerical solution for SD-ROM with online stabilizing post-processing at T=1forr=30,60,90 (from top to bottom). Table 4 Example 4.2.2: L2-norm of the deviation from the final exact solution profile along the mean diagonal for different FOMs. Offline methods eFOM 0,ν=10−8 DNS-FEM 0.1816 DNS-FEM post-processing 0.1345 LPS-FEM 0.1247 LPS-FEM post-processing 0.0393 decrease from 12% to 4% when applying stabilizing post-processing to LPS-FEM, as shown in Table 4. Also, if we proceed by constructing POD basis from LPS-FEM (without stabilizing post-processing), being more influenced by spurious oscillations, it leads to online numerical solutions that are globally polluted with high spurious oscillations even for r=90, whatever it is the applied reduced order system, as shown in Fig. 17. Thus, we decided to proceed by constructing POD basis by using LPS-FEM with stabilizing post-processing, to limit the influence of POD noisy data in the online phase. In Fig. 18, we show the decay of POD eigenvalues associated both to the snapshots correlation matrix (2.8) and the advection correlation matrix (2.13)in this case. Again, one can observe that the decay of the POD eigenvalues associated to the advection correlation matrix is rather slow, due to the very low diffusion. However, adding the corresponding stabilization term in the online phase greatly improves the results over the standard POD-ROM also in this case. 21 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 16. Example 4.2.2: Final solution profiles along the mean diagonal for different FOMs. Fig. 17. Example 4.2.2: Final solution profiles along the mean diagonal for different ROMs at r=90 using noisy POD data from LPS-FEM. Fig. 19 presents results for all considered ROMs: G-ROM, G-ROM post-processing, SD-ROM, and SD-ROM post-processing. One can observe that results for G-ROM (with and without online a-posteriori stabilization) are globally quite oscillatory, even at r=90. However, applying SD-ROM already localizes oscillations just near the moving steep layer, and also SD-ROM with online stabilizing post-processing allows to further improve results, maintaining the amplitude of oscillations in a reasonable low range. This is reflected by results depicted in Table 5. One can see that, for r=90, SD-ROM post-processing method approaches the accuracy of the offline phase by considerably suppressing the influence of noisy modes. Comparing also to Table 3(Case ν=10−6), the SD-ROM performs well for the different values of νtested and displays a low sensitivity with respect to changes in the diffusion coefficient. This also provide a numerical support for the theoretical error estimate derived in [48], which is uniform with respect to ν(see Remark 2.2). Again, note that although the first r=30 POD modes already capture more than 99% of the system’s kinetic energy, all ROMs yield poor quality results for which the peak of the front is not reached, and they display globally spread numerical oscillations, reflecting the extreme complexity of the problem. Augmenting the number of POD modes allows to reach the peak of the front for all methods. However, whereas the solution of the G-ROM (with and without online a-posteriori stabilization in this case) remains globally polluted with spurious oscillations, the SD-ROM notably reduces the amplitude of oscillations, and its combination with online stabilizing post-processing allows to compute a rather accurate solution in this case, comparable with the one of the offline phase. In 22 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 18. Example 4.2.2:PODeigenvalues. Table 5 Example 4.2.2: Captured system’s kinetic energy and L2-norm of the deviation from the final exact solution profile along the mean diagonal for different ROMs at r=30, 60, 90. ν=10−8r=30 r=60 r=90 Captured system’s Ekin(%)99.71 99.96 >99.99 ν=10−8eROM 0 Online methods r=30 r=60 r=90 G-ROM 0.3733 0.1676 0.1224 G-ROM post-processing 0.3086 0.1493 0.0884 SD-ROM 0.3417 0.1463 0.0675 SD-ROM post-processing 0.2596 0.1449 0.0589 Fig. 20, we show the numerical solution at T=1for the best performing SD-ROM with online a-posteriori stabilization for r=30, 60, 90 (from top to bottom). Again, with this method, numerical unphysical oscillations are practically eliminated by gradually increasing the number of POD modes. 5. Summary and conclusions In this work, we have proposed to improve the stabilized POD-ROM introduced in [48]to deal with the numerical simulation of advection-dominated advection-diffusion-reaction equations. In particular, we have proposed a three-stage stabilizing strategy that has proved to be very useful when considering very low diffusion coefficients, i.e. in the strongly advection-dominated regime. This approach mainly consists in three ingredients: (1) the addition of a “streamline diffusion” stabilization term to the governing projected equations, (2) the modification of the correlation matrix defining the POD modes associated to the advection stabilization term, and (3) an a-posteriori stabilization scheme. The performed numerical studies have shown the potential of the new ROM in handling strongly advection-dominated cases, also tested for long time integrations on periodic systems, by extremely limiting spurious oscillations and thus obtaining rather accurate results in this framework. To remove the few remaining oscillations, one could think to apply more complex shock or discontinuity capturing methods (see [39]for a detailed review) and try to adapt them to the POD-ROM framework as future interesting research topic. Also, one could carry out a similar numerical investigation of the significantly more challenging Navier–Stokes equations in view of computing more complex convection-dominated and turbulent flows. Another interesting research direction could be to test the proposed method in the predictive regime for test cases not resembling a periodic behavior, such as the test case in section 4.2. In this case, one should endow the SD-ROM with a basis updating mechanism in order to get acceptable errors in this regime, using for instance a-posteriori error indicators. This study is in progress, following some hints given by the hybrid DNS/POD approach introduced in [14]. Apart from prediction in time considered in the present work, we are interested in extending the proposed method in order to make predictions across geometrical and/or physical parameters (see, e.g., [54]), of interest to solve engineering problems such as shape optimization and flow control. 23 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 19. Example 4.2.2: Final solution profiles along the mean diagonal for different ROMs at r=30,60,90 (from top to bottom). 24 M. Azaïez, T. Chacón Rebollo and S. Rubino Journal of Computational Physics 425 (2021) 109916 Fig. 20. Example 4.2.2: Numerical solution for SD-ROM with online stabilizing post-processing at T=1forr=30,60,90 (from top to bottom). CRediT authorship contribution statement All authors contribute to this paper as a whole. All authors read and approved the final manuscript. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work has been partially supported by the Spanish Government - FEDER EU grant RTI2018-093521-B-C31 and European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA). The first author has received financial support from the French State in the frame of the Investments for the future IdEx Bordeaux (Initiative d’Excellence de l’Université de Bordeaux) Programme, reference ANR-10-IDEX-03-02. The third author would gratefully acknowledge the financial support received from IdEx Bordeaux International Post-Doc Programme during his initial postdoctoral research involved in this article. The research of the third author has been also funded by the Spanish State Research Agency through the national programme Juan de la Cierva-incorporación 2017. References [1] N. Ahmed, T. Chacón Rebollo, V. John, S. Rubino, Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method, IMA J. Numer. Anal. 37 (3) (2017) 1437–1467. 25