A STREAMLINE DERIVATIVE POD-ROM FOR ADVECTION-DIFFUSION-REACTION EQUATIONS

. We introduce a new streamline derivative projection-based closure modeling strategy for the numerical stabilization of Proper Orthogonal Decomposition-Reduced Order Models (POD-ROM). As a ﬁrst preliminary step, the proposed model is analyzed and tested for advection-dominated advection-diﬀusion-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 eﬃciency 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 diﬀusion coeﬃcients.


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 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: where b with b ∞ = O(1) is the given advective field, ε << 1 the diffusion parameter, g the reaction coefficient, f the forcing term, Ω the computational domain in R d (d = 2 or 3), t ∈ [0, T ], with T the final time, and u 0 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 ε << b ∞ ).
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.

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 the best approximation: in a real Hilbert space H. Although H can be any real Hilbert space, in what follows we consider H = L 2 (Ω), In the framework of the numerical solution of Partial Differential Equations (PDE), u is usually given at a finite number of times t 0 , . . . , t N , the so-called snapshots. Let us consider an ensemble of snapshots χ = span {u(·, t 0 ), . . . , u(·, t N )}, which is a collection of data from either numerical simulation results or experimental observations at time t n = n∆t, n = 0, 1, . . . , N and ∆t = T /N . Then, usually an approximation of the error in the square of the L 2 (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]: 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: where K ∈ R (N +1)×(N +1) is the snapshots correlation matrix with entries: (u(·, t n ), u(·, t m )) , for m, n = 0, . . . , N, z i is the i-th eigenvector, and λ i is 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: where (z i ) n is the n-th component of the eigenvector z i . It can also be shown that the following POD error formula holds [24,29]: where M is the rank of χ.
We consider the following space for the POD setting: 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 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 r of 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: where {a i (t)} r i=1 are the sought time-varying coefficients representing the POD-Galerkin trajectories. Note that r << N , where N denotes the number of degrees of freedom (d.o.f.) in a full order simulation (e.g., DNS). Replacing u with u r in (1), using the Galerkin method, and projecting the resulted equations onto the space X r , one obtains the standard POD-ROM: 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.

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 {T h } h>0 be a regular family of triangulations of Ω. For any mesh cell K ∈ T h , its diameter will be denoted by h K and h = max K∈T h h K . To describe our strategy, we define the scalar product: and its associated norm where for any K ∈ T h , τ K is a positive local stabilization parameter (to be determined later). Let us introduce the POD space: where ϕ i , i = 1, . . . , R, are the POD modes associated to K, defined as the snapshots correlation matrix with entries: Note that for classical POD modes asssociated to the standard correlation matrix K mn , there already exists a theory on convergence rates and error bounds for POD expansions of parameterized solutions of heat equations, see e.g. [6][7][8]. With co-authors of the referred works, we aim to derive a similar analysis for POD modes associated to the advection correlation matrix K mn defined in (9).
We consider the L 2 -orthogonal projection on X R , P R : Let P R = I − P R , where I is the identity operator. We propose the Streamline Derivative projection-based POD-ROM (SD-POD-ROM) for (1): Remark 1.2. When τ K = 0 for any K ∈ T h , 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: 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 M S − P OD − 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: α being a constant eddy viscosity coefficient, and P R = I − P R , with P R the L 2 -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.

Error estimates
In this section, we prove estimates for the average error: where the approximation u n is the weak solution of (1) and u n r is the solution of (11) discretized by FE, at time t n . 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.

Technical background
This section provides some technical results that are required for the numerical analysis. Throughout the paper, we shall denote by C, C 1 , C 2 , . . . 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 C from ε, 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)
For the FE discretization of the weak form of (1), we consider a family of finite dimensional subspaces X h of X = H 1 0 (Ω) such that, for all v ∈ H m+1 (Ω) ∩ X, the following assumption is satisfied.
where k is the order of accuracy of X h .
For the POD approximation, the following POD inverse estimate was proven in [29], Lemma 2: . . , r, be POD modes, M r be the POD mass matrix with entries [M r ] ij = (ϕ j , ϕ i ), and S r be the stiffness matrix with entries [S r ] ij = (∇ϕ j , ∇ϕ i ), i, j = 1, . . . , r. Let · 2 denote the matrix 2-norm. Then, for all ϕ r ∈ X r , the following estimate holds: Note that, since we have chosen H = L 2 in the POD method, M −1 r 2 = 1 in inequality (15). To prove optimal error estimates in time, we follow [29] and include the finite difference quotient∂u(t n ) = u n − u n−1 ∆t , for n = 1, . . . , N , in the set of snapshots χ = {u(t 0 ), . . . , u(t N ),∂u(t 1 ), . . . ,∂u(t N )}. As pointed out in [29], the POD error formula (6) becomes: For the subsequent numerical analysis, we need the following technical hypothesis on the stabilization parameters τ K : Hypothesis 2.4. The stabilization parameters τ K satisfy the following condition: for all K ∈ T h , and a positive constant C independent of h and ε.
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 ∈ L 2 (Ω), the following estimate is satisfied: Proof. By using (17) and the stability of P R in the L 2 -norm, it follows: Thus, the estimate (18) can be deduced.
Proof. Hypotheses 2.1 and 2.4 guarantee the well-posedness of (19). To prove estimate (20), we choose ϕ r = u n+1 r in (19), and decompose the bilinear form A into its symmetric and skew-symmetric parts: Using the identity: and Young's inequality, from (21) we get: 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 u n − u n r , we will first consider the Ritz projection w r ∈ X r of u ∈ X: The existence and uniqueness of w r follow from the Lax-Milgram lemma. We now state an estimate for u n −w n r , the error in the Ritz projection.
Lemma 2.9. The Ritz projection w n r of u n satisfies the following error estimate: 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 w n r of u n satisfies the following error estimate up tp O(∆t 2 ): 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(t n ) are included in the set of snapshots (see remark 1 in [29]).

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): Theorem 2.11. The solution of the SD-POD-ROM (19) satisfies the following error estimate: where the initial condition is given by the L 2 -orthogonal projection of u 0 on X r : and λ i , i = R+1, . . . , M , in the right-hand side of (26) are the eigenvalues associated to the snapshots correlation matrix K previously defined in (9).
Proof. We evaluate the weak form of (1) at t = t n+1 , let the function test v = ϕ r , and add and subtract the different quotient term u n+1 − u n ∆t : 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: (29) We now decompose the error as u n − u n r = (u n − w n r ) − (u n r − w n r ) = η n − φ n r , so that by triangle inequality we have: Note that η n has already been bounded in lemma 2.9. Thus, in order to estimate the error, we only need to estimate φ n r . The error equation (29) can be written as: (31) We consider ϕ r = φ n+1 r , and note that, since φ n+1 r ∈ X r , then A(η n+1 , φ n+1 r ) = 0, so that we get: where we have denoted r n = ∂ t u n+1 − u n+1 − u n ∆t . We proceed by estimating all the terms in (32). The terms on the left-hand side of (32) are estimated as follows: By using Cauchy-Schwarz and Young's inequalities, the terms on the right-hand side of (32) are estimated as follows: Using (33)- (36) and absorbing right-hand side terms into left-hand side terms, (32) becomes: By using Young's inequality, the first term on the left-hand side of (37) can be estimated as follows: Using (38) in (37) and multiplying by 2∆t, we get: Summing from n = 0 to k ≤ N − 1 in (39), we obtain: The second term on the right-hand side of (40) can be estimated as follows: which has been bounded in corollary 2.10. For the third term on the right-hand side of (40), we get: Finally, the last term on the right-hand side of (40) can be estimated as follows: where we have used assumption (17) in hypothesis 2.4 on the stabilization parameters, assumption (14) in hypothesis 2.2 on the FE approximation, (10), (6), and λ i , i = R + 1, . . . , M , are the eigenvalues associated to the snapshots correlation matrix K previously defined in (9). Using (41)-(43) in (40), (30) and the estimate (24) in lemma 2.9, we obtain the error estimate (26). This concludes the proof.
Remark 2.12. If one just consider the standard Galerkin POD-ROM (τ K = 0 for any K ∈ T h ), thus error estimate (26) can be recovered, without the appearance of the last term on the right-hand side of (26). In this case, any control on the high-frequency modes of the advective derivative is guaranteed. When R = 0, one has that the last term on the right-hand side of (26) is limited to √ h. This low convergence order appears linked to the diffusive nature of the penalty-stabilized POD-ROM (12), which extends the numerical diffusion to all the resolved modes.

Numerical results
The mathematical model used for the numerical tests in this section is the advection-dominated advectiondiffusion-reaction equation (1)  sharp internal layer moving in time. As in the theoretical developments in section 2, in this section we employ the FE method for spatial discretization and the backward Euler method for temporal discretization. The open-source FE software FreeFem++ [23] has been used to run the numerical experiments. All computations are carried out on a MacBook Pro with a 3.1 GHz Intel Core i7 processor.
We perform a comparison between the SD-POD-ROM (11) and the standard POD-ROM (8). To generate the POD modes, we run offline a DNS with the following parameters: piecewise linear FE, uniform triangular mesh with mesh-size h = 10 −2 , and time step ∆t = 10 −3 . The POD modes are generated in L 2 by the method of snapshots by storing every tenth solution, so that 101 snapshots were used. The DNS average error  with r = 40 (middle). It is clear from this figure that, although the first 40 POD modes capture 99.96% of the system's kinetic energy, the standard POD-ROM yields poor quality results and displays visible numerical oscillations. This is confirmed by the standard POD-ROM high average error, which is almost one orders of magnitude higher than that of the performed SD-POD-ROM for r = 40. It happens that r = 40 in POD-ROM almost provides similar accuracy of r = 30 in SD-POD-ROM. It is thus clear that the standard POD-ROM, although computationally efficient (CPU times three orders of magnitude lower than DNS, at least), is rather inaccurate for advection-dominated configurations.   Now, we investigate the SD-POD-ROM (11). We make the following parameter choice: R = r/2. Also, the working expression of the stabilization coefficients is: by following the form proposed in [17,18], designed by asymptotic scaling arguments applied in the framework of stabilized methods aimed at taking into account the local balance between advection, diffusion and reaction.
In expression (44), c 1 , c 2 and c 3 are positive algorithmic constants, and U K is some local advection speed on the mesh cell K. The values of the constants are chosen to be c 1 = 4, c 2 = √ c 1 = 2, c 3 = 1, cf. [16].
These values are justified from the analysis of the one-dimensional advection-diffusion-reaction equation and from many numerical experiments, for which they are optimal, cf. [16]. In this case, h K = h, and we take U K = b ∞,K = b ∞ (= sin π/3), so that for all K ∈ T h : and assumption (17) figure 3 (right). We can observe from this figure that the SD-POD-ROM is more stable and accurate than the standard one, and numerical unphysical oscillations displayed by the latter are practically eliminated by adding numerical stabilization, while keeping the same level of computational efficiency. This is confirmed by figure 2 and table 2.
Also, if we slightly extend the time range used in the generation of the POD modes (just up to T = 1.25), we have that the final solution for the standard POD-ROM is already totally inaccurate and oscillatory, while the one for the new SD-POD-ROM is still rather acceptable, see figure 4. This suggests that the new SD-POD-ROM could be considered in principle a better predictive tool with respect to the standard POD-ROM. Nevertheless, for much larger time intervals than that used in the derivation of the input data, one should endow the new SD-POD-ROM with a basis updating mechanism, using for instance a posteriori error indicators. This study is today in progress, following some hints given by the hybrid DNS/POD approach introduced in [12].

Conclusions
In this work, we have proposed a new stabilized POD-ROM for the numerical simulation of advectiondominated advection-diffusion-reaction equations. This model, denoted SD-POD-ROM, is derived from highorder stabilized FE methods, and uses a streamline derivative projection-based operator to properly take into account the high frequencies advective derivative component of POD modes not included in the ROM.
We have performed a thorough numerical analysis of the arising fully discrete SD-POD-ROM applied to advection-dominated advection-diffusion-reaction problems. In particular, the numerical analysis makes apparent an extra-control on the high frequencies of the advective derivative, which is an extremely important feature in view of computing more complex convection-dominated and turbulent flows. We also emphasize that the theoretical error estimates are uniform with respect to the diffusion coefficient.
At a computational level, the new SD-POD-ROM has been tested on a representative problem displaying a sharp internal layer advected in time. We employed the theoretical error estimates to provide some guidance in choosing the stabilization parameters in practical computations. The numerical investigations yielded the following conclusions: the SD-POD-ROM is more stable and accurate than the standard POD-ROM, and numerical unphysical oscillations displayed by the latter are practically eliminated by the former, while keeping the same level of computational efficiency. Also, the new SD-POD-ROM showed better predictive features for short-time integrations of unsteady fields, which is rather promising in view of potentially becoming a powerful predictive tool in realistic physical and engineering applications.