Pressure stabilization strategies for a LES filtering Reduced Order Model

We present a stabilized POD-Galerkin reduced order method (ROM) for a Leray model. For the implementation of the model, we combine a two-step algorithm called Evolve-Filter (EF) with a computationally efficient finite volume method. In both steps of the EF algorithm, velocity and pressure fields are approximated using different POD basis and coefficients. To achieve pressure stabilization, we consider and compare two strategies: the pressure Poisson equation and the supremizer enrichment of the velocity space. We show that the evolve and filtered velocity spaces have to be enriched with the supremizer solutions related to both evolve and filter pressure fields in order to obtain stable and accurate solutions with the supremizer enrichment method. We test our ROM approach on 2D unsteady flow past a cylinder at Reynolds number 0<= Re<= 100. We find that both stabilization strategies produce comparable errors in the reconstruction of the lift and drag coefficients, with the pressure Poisson equation method being more computationally efficient.


Introduction
For about a couple of decades, reduced order models (ROMs) have emerged as an efficient tool for the approximation of problems governed by parametrized partial differential equations. This success is owed to the fact that ROMs can significantly reduce the computational cost required by classical full order models (FOMs), e.g. Finite Element Method or Finite Volume Method, when several solutions associated to different parameter values are needed. The basic ROM framework consists of two steps. During a first phase (called offline), a database of several solutions is collected by solving a FOM of choice for different parameter values. Then, during a second phase (called online) the database collected in the offline phase is used to compute the solution for newly specified values of the parameters in a short amount of time. For a comprehensive review on ROMs, we refer the reader to, e.g., [26,40,6,5,3,7].
The particular ROM we consider in this paper is based on a POD-Galerkin approach, which consists in extracting of the most energetic modes representing the system dynamics and projecting the governing equations onto the space spanned by these modes. For the specific application that we target, i.e. incompressible fluid flow at moderately high Reynolds number, it is well known that standard POD-Galerkin models lead to instabilities [14,37]. A successful way to cure these instabilities in advection dominated flows is to adopt subgridscale closure models. See, e.g., [51,2]. Thus, we choose to work with a Large Eddy Simulation (LES) approach.
We focus on a variant of the so-called Leray model [34], where the small-scale effects are described by a set of equations to be added to the discrete Navier-Stokes equations. This extra problem acts as a differential low-pass filter [11]. For its actual implementation, we use the Evolve-Filter (EF) algorithm [12,19,18,33]. One of the novelties of our approach is that we use a Finite Volume (FV) method [22,21], while the vast majority of the works on Leraytype models use a Finite Element framework. In the context of regularized ROMs, the Leray model and the EF algorithm have been thoroughly investigated. See, e.g., [25,55,53,24,54]. In all of these works, the filtering approach is only employed at reduced order level. In [23], for the first time the LES filtering is used also at the FOM level, i.e. to generate the snapshot data. Such an approach provides a ROM which is fully consistent with the FOM as the same mathematical framework is used during both the offline and online stages.
This paper could be seen as an extension of [23]. Therein, we used a pressure Poisson equation (PPE) in the online stage [1,46] as a pressure stabilization technique. Here, we compare the PPE method with the supremizer enrichment of the velocity space [44,4,46,20], i.e. another techniques that provides pressure stability. The main objective of this work is to test the accuracy and efficiency of these two methods within our LES filtering approach. We show that adapting the supremizer enrichment to the EF algorithm is not a trivial exercise. Indeed, the supremizer method becomes accurate only when the evolve and filtered velocity spaces are enriched by the supremizer solutions associated to both evolve and filter pressure fields. Since we are interested in the reconstruction of the pressure field at reduced order level, time is the only parameter we consider. We vary no physical and/or geometrical parameters. We test our framework on 2D flow past a cylinder with time-dependent Reynolds number 0 ≤ Re(t) ≤ 100 [49,28].
The work is organized as follows. Sec. 2 describes the full order model and the numerical method we use for it. In Sec. 3, we describe the reduced order model. The numerical examples are reported in Sec. 4. Finally, Sec. 5 provides conclusions and perspectives.

The full order model
In this section, we describe our Full Order Model (FOM). We consider a fixed domain Ω ⊂ R D with D = 2, 3 over a time interval of interest (t 0 , T ) ⊂ R + . The so-called Leray model couples the Navier-Stokes equations (NSE) with a differential filter as follows: where ρ is the fluid density, µ is the dynamic viscosity, u is the fluid velocity, p is the fluid pressure, u is the filtered velocity, and variable λ is a Lagrange multiplier to enforce the incompressibility constraint for u. Parameter α in eq. (3) can be interpreted as a filtering radius. In this paper, α will be constant in space and time. More sophisticated choices are possible but will not be considered here. Problem (1)-(4) is endowed with suitable boundary conditions (2α 2 ∇u − λI)n = 0 on ∂Ω N × (t 0 , T ). (8) and the initial data In addition, u D and u 0 are given. Note that we restrict our attention to homogeneous Neumann boundary conditions. Of course, the methodology we propose can be extended to non-homogeneous Neumann conditions, as well as other boundary conditions (e.g., Robin conditions). [9] 2.1 The Evolve-Filter algorithm We start with the time discretization of model (1)-(4). Let ∆t ∈ R, t n = t 0 + n∆t, with n = 0, ..., N T and T = t 0 + N T ∆t. Moreover, we denote by y n the approximation of a generic quantity y at the time t n .
-Filter : find (u n+1 , λ n+1 ) such that with boundary conditions We consider u n+1 and q n+1 the approximation of the fluid velocity and pressure at the time t n+1 , respectively.
Remark 2.1. Filter problem (13)- (14) can be considered a generalized Stokes problem. In fact, if we multiply eq. (13) by ρ/∆t and rearrange the terms we obtain: where q n+1 = ρλ n+1 /∆t can be seen as a filtered pressure. Problem (17), (14) can be seen as a time dependent Stokes problem discretized by the Backward Euler (or BDF1) scheme and with viscosity µ. A solver for problem (17), (14) can be obtained by adapting a standard linearized Navier-Stokes solver.

Space discrete problem by a Finite Volume method
For the space discretization of problems (9)- (10) and (17), (14), we adopt a FV method. We partition the computational domain Ω into cells or control volumes Ω i , with i = 1, . . . , N c , where N c is the total number of cells in the mesh. Let A j be the surface vector of each face of the control volume, with j = 1, . . . , M . The fully discretized form of problem (9)-(10) is given by In (18)- (20), v n+1 i and b n+1 i denote the average velocity and source term in control volume Ω i , respectively. Moreover, we denote with v n+1 i,j and q n+1 i,j the velocity and pressure associated to the centroid of face j normalized by the volume of Ω i .
The fully discrete problem associated to the filter problem (17), (14) is given by In (21)-(23), we denoted with u n+1 i the average filtered velocity in control volume Ω i , while q n+1 i,j is the auxiliary pressure at the centroid of face j normalized by the volume of Ω i . For more details on the full discretization of both problems (18)- (20) and (21)- (23), we refer the reader to [22].
We have implemented the EF algorithm within the C++ finite volume library OpenFOAM ® [52]. For the solution of the linear system associated with (18)- (19) we used the PISO algorithm [27], while for problem (21)-(22) we chose a slightly modified version of the SIMPLE algorithm [39], called SIMPLEC algorithm [50]. Both PISO and SIMPLEC are partitioned algorithms that decouple the computation of the pressure from the computation of the velocity.

The reduced order model
The Reduced Order Model (ROM) we propose is based on a framework introduced in [23]. Here in Sec 3.1 we describe the procedure we use to construct a POD-Galerkin ROM, while in Sec. 3.2 we present two different strategies for pressure stabilization at reduced order level. Finally, Sec. 3.3 describes the method we apply to enforce non-homogeneous Dirichlet boundary conditions (5), (7) at the reduced order level. The ROM computations are carried out using ITHACA-FV [45], an in-house open source C++ library.

A POD-Galerkin projection method
We approximate velocity fields v and u and pressure fields q and q as linear combinations of the dominant modes (basis functions), which are assumed to be dependent on space variables only, multiplied by scalar coefficients that depend on the time: In (24)- (25), N Φr denotes the cardinality of a reduced basis for the space field Φ belongs to. Note that Φ could be either a scalar or a vector field. Using (24) to approximate v n+1 and q n+1 in (9)-(10), we obtain where u * r = 2u n r − u n−1 r and b n+1 r = (4u n r − u n−1 r )/(2∆t). Then, using (25) to approximate u n+1 and q n+1 in (17), (14) we get: Several techniques to generate the reduced basis spaces are available in the literature, e.g. the Proper Orthogonal Decomposition (POD), the Proper Generalized Decomposition (PGD) and the Reduced Basis (RB) with a greedy sampling strategy. See, e.g., [43,15,30,40,16,17,48,7]. We find the reduced basis by using the method of snapshots. To this purpose, we solve the FOM described in Sec. 2 for each time t k ∈ {t 1 , . . . , t Ns } ⊂ (t 0 , T ]. The snapshots matrices are obtained from the full-order snapshots: where the subscript h denotes a solution computed with the FOM and N Φ h is the dimension of the space field Φ belong to in the FOM. The POD problem consists in finding, for each value of the dimension of the POD space N P OD = 1, . . . , N s , the scalar coefficients a 1 1 , . . . , a Ns 1 , . . . , a 1 Ns , . . . , a Ns Ns and functions ζ 1 , . . . , ζ Ns that minimize the error between the snapshots and their projection onto the POD basis. In the L 2 -norm, we have It can be shown [31] that eq. (31) is equivalent to the following eigenvalue problem where C Φ is the correlation matrix computed from the snapshot matrix S Φ , Q Φ is the matrix of eigenvectors and Λ Φ is a diagonal matrix whose diagonal entries are the eigenvalues of C Φ . Then, the basis functions are obtained as follows: The POD modes resulting from the aforementioned methodology are: where N Φr < N s are chosen according to the eigenvalue decay. The reduced order model can be obtained through a Galerkin projection of the governing equations onto the POD spaces. Let where ϕ i and ψ i are the basis functions in (24). At time t n+1 , the reduced algebraic system for problem (26)- (27) is: where vectors β n+1 and γ n+1 contain the values of coefficients β i and γ i in (24) at time t n+1 . The term G r (β n , β n−1 )β n+1 in (38) is related to the non-linear convective term: where G r is a third-order tensor defined as follows See [41,42] for more details. Next, let where ϕ i and ψ i are the basis functions in (25). At time t n+1 , the reduced algebraic system for problem (28)-(29) is where vectors β n+1 and γ n+1 contain the values of coefficients β i and γ i in (25) at time t n+1 .
Finally, the initial conditions for the ROM algebraic system (38)-(39) are obtained performing a Galerkin projection of the initial full order condition onto the POD basis spaces: . The complete reduced algebraic system at time t n+1 is given by (38)

Pressure fields reconstruction and pressure stability
It is well known that the reduced problem (38), (39), (44), (45) presents stability issues because the approximation spaces need to satisfy the inf-sup (Ladyzhenskaya-Brezzi-Babuska) condition [13,10]. In a standard finite element (FE) NSE framework, the inf-sup condition reads: inf where Q h is the FE space for the pressure approximation, V h is the FE space for the approximation of the velocity field, and γ is a constant that does not depend on the mesh size h.
In order to obtain a stable and accurate reconstruction of the pressure field at the reduced level, different approaches have been proposed. One option is to use a global POD basis for both pressure and velocity field and same temporal coefficients [8,36]. Another option is represented by the supremizer enrichment method; see, e.g., [44,4,46,20]. Finally, one can take the divergence of the momentum equation to obtain a Poisson equation for the pressure that is projected onto a POD basis; see, e.g., [1,46]. This third method is called Poisson pressure equation (PPE).
In this work, we test and compare the performances of two methods: the PPE method (already combined with the EF algorithm in [23]) and the supremizer enrichment method (not yet tested for the EF algorithm). As we will show, the extension of the supremizer enrichment method to the EF algorithm is not straightforward.

Pressure Poisson Equation method
We take the divergence of eq. (9) and account for conditions (10) to obtain the Poisson pressure equation for the Evolve step: So, at the Evolve step instead of solving (9)-(10), we solve the modified systems (9), (47) with boundary condition (11) and where ∂ n denotes the derivative with respect to outgoing normal n. We proceed similarly for the Filter step. We take the divergence of eq. (17) and account for condition (14) to obtain the Poisson pressure equation: The Filter step becomes solving (17), (49) with boundary condition (15) and The reader interested in enforcing non-homogeneous Neumann conditions for the pressure fields is referred to [38,29]. Remark 3.1. Systems (9)- (10) and (17), (14) are not equivalent to systems (9), (47) and (17), (49) for steady flows [35,38,29]. As discussed in Remark 2.1, the filter problem can be seen as a time-dependent Stokes problem.
Using (24) to approximate v n+1 and q n+1 in (47), we obtain ∆q n+1 After space discretization, the matrix form of eq. (51) reads: where The residual associated with the non-linear term in the equation (52) is evaluated with the same strategy used for eq. (38). We have where J r is a third-order tensor defined as follows Using (25) to approximate q n+1 in eq. (49), we get ∆q n+1 Once discretized in space, eq. (57) can be written in the matrix form as: where With the PPE method, at every time step the ROM algebraic system that has to be solved is (38), (52), (44), (58).

Supremizer enrichment method
For a given pressure basis function, the supremizer is the solution that permits the realization of the inf-sup condition (46). One has to find the supremizer for each pressure basis function. Here, we use an approximated supremizer enrichment procedure: instead of pressure basis functions we use pressure snapshots. This procedure allows to drastically reduce the online computational cost. In fact, the supremizer basis functions do not depend on the particular pressure basis functions but are computed directly from the pressure snapshots during the offline phase. The downside of this approximated procedure is that it is not possible to rigorously show that the inf-sup condition is satisfied. One only relies on heuristic criteria or checks during a post-processing stage.
For NSE, the supremizer s(t i ) relative to pressure snapshot q(t i ) is found by solving the following problem: in Ω, (60) for i = 1, . . . , N s . For more details, the reader is referred to [44,4,46,20]. For the EF algorithm, in addition to (60)-(61) one has to solve an analogous problem for each auxiliary pressure snapshot q(t i ): A POD procedure is applied to the matrices above in order to obtain the supremizer POD basis functions: where N sr and N sr conform to the notation introduced in Sec. 3.1. These velocity supremizer basis functions are added to the reduced velocity spaces, which become: Let us call SUP 1 the supremizer enrichment method described above. In Sec. 4, we will show that SUP 1 does not lead to an accurate reconstruction of the pressure fields. To increase the accuracy, we add the supremizer solutions associated to both pressure fields to the evolve and filtered velocity spaces, i.e.: We call SUP 2 this realization of the supremizer enrichment. With either SUP 1 or SUP 2 , the ROM algebraic system that has to be solved at every time step is (38), (39), (44), (45).

Treatment of the Dirichlet boundary conditions: the lifting function method
In order to homogeneize the velocity fields snapshots and make them independent of the boundary conditions, we use the lifting function method [46]. The lifting functions are problem-dependent: they have to be divergence free in order to retain the divergence-free property of the basis functions and they have to satisfy the boundary conditions of the FOM.
The velocity snapshots are modified as follows: where N BC is the number of non-homogeneous Dirichlet boundary conditions, χ j (x) are the lifting functions, and v BC j and u BC j are suitable temporal coefficients. The POD is applied to the snapshots satisfying the homogeneous boundary conditions. Then, the boundary value is added back in this way:

Numerical results
We consider 2D flow past a cylinder [28,49]. This well-know benchmark will allow us to assess the ability of the ROM approaches presented in Sec. 3 to reconstruct the time evolution of the velocity and pressure fields. We have thoroughly investigated the 2D flow past a cylinder with a finite volume FOM in [22] and with a ROM using a PPE approach in [23]. Here, we aim at comparing the supremizer enrichment method and the PPE method in terms of errors and speed-up. The computational domain is a 2.2 × 0.41 rectangular channel with a cylinder of radius 0.05 centered at (0.2, 0.2), when taking the bottom left corner of the channel as the origin of the axes. The channel is filled with fluid with density ρ = 1 and viscosity µ = 10 −3 . We impose a no slip boundary condition on the upper and lower walls and on the cylinder. At the inflow, we prescribe the following velocity profile: and ∂q/∂n = ∂q/∂n = 0. At the outflow, we prescribe ∇v · n = 0 and q = q = 0. We start the simulations from fluid at rest. Note that the Reynolds number is time dependent, with 0 ≤ Re(t) ≤ 100 [49]. We consider a hexaedral computational grid with h min = 4.2e − 3 and h max = 1.1e − 2 for a total of 1.59e4 cells. The quality of this mesh is high: it features very low values of maximum non-orthogonality (36 • ), average non-orthogonality (4 • ), skewness (0.7), and maximum aspect ratio (2). Fig. 1 (left) shows a part of the mesh. We chose this mesh because it is the coarsest among all the meshes considered in [22] and thus the most challenging for our LES approach. We fix the time step to ∆t = 4e − 4 for both FOM and ROM simulations. For the convective term, we use a second-order accurate Central Differencing scheme [32]. In this way, we avoid introducing stabilization and are thus able to assess the effect of the filter. We need to enforce a time-dependent non-uniform Dirichlet boundary condition at the inlet. For this purpose, we consider a divergence free function with the following non-uniform velocity distribution: and uniform null values on the rest of the boundary. See Figure 1 (right). We will compare our findings with those in [46]. The choice of Ref. [46] is due to the fact that therein the authors develop a NSE-ROM finite volume framework both with PPE and supremizer enrichment methods for the reconstruction of the pressure field.
We set α = 0.0032 and refer to [22] for details on this choice. The snapshots are collected every 0.1 s (i.e., N s = 80) using an equispaced grid method in time. Therefore, the dimension of the correlation matrix C Φ in (32) is 80 × 80. For a convergence test as the number of snapshots increases, see [23]. Table 1 reports the first 4 cumulative eigenvalues for the velocity, pressure, and supremizer fields. In order to retain 99% of the energy for the ROM, we need 2 modes for v, 2 modes for q, 2 modes for u, and 1 mode for q. As for SUP 1 , we consider a number of supremizer modes greater than pressure modes as suggested in [46]: 4 modes for s and 3 modes for s. On the other hand, for SUP 2 we take into account an equal number of pressure and supremizer modes: 2 modes for s and 1 mode for s.  We calculate the L 2 relative error: where Φ h and Φ r are the FOM approximation of a given field (i.e., v h , u h , q h or q h ) and the corresponding ROM approximation (i.e., v r , u r , q r or q r ), respectively. Fig. 2 shows error (70) for the two velocity and pressure fields over time for the three different ROM techniques under investigation: PPE, SUP 1 and SUP 2 . We also present the errors related to the ROM computations with no stabilization technique, referred to as NOS. As one would expect, Fig. 2 shows that the model with no stabilization is completely unreliable. From Fig. 2, we see that the SUP 1 -ROM is a big improvement over NOS-ROM. However, while the errors for the velocity fields are acceptable, the pressure errors remain large and far above the values obtained with NSE in [46]. This shows that SUP 1 -ROM is not a reliable stabilization of the ROM for the EF algorithm. The SUP 2 -ROM produces much better results in terms of the pressure fields, with errors for the velocity fields that are comparable with those given by the SUP 1 -ROM. We speculate that the better performance of the SUP 2 model (which adds the supremizer solutions related to both pressure fields to the evolve and filtered velocity spaces) with respect to the SUP 1 model could be due to the strong coupling between the evolve velocity v and the filter velocity u, However, the stability of the inf-sup ROM formulation for the EF algorithm needs to be investigate in more depth and will be object of future work.
Finally, from Fig. 2 we observe that the PPE-ROM provides the lowest errors for the velocity fields, while the pressure errors are comparable (for q) or worse (for q) than the errors given by the SUP 2 -ROM. Our findings for the EF algorithm are in agreement with what observed in [46] for NSE. Indeed, therein it is shown that the PPE model produces better results for the velocity field but worse results for the pressure field when compared to the supremizer enrichment model. For insights on the behavior of the errors at the first and last time steps of the simulation, we refer to [23]. We report minimum, average, and maximum relative errors for PPE and SUP 2 in Table  2. The average errors for u and q are comparable to the values obtained in [46].  For a visual comparison, we report in Fig. 3 velocity and pressure fields at t = 5 computed by FOM, PPE-ROM, and SUP 2 -ROM. We observe that both PPE-ROM and SUP 2 -ROM can capture well the main flow features.
For a quantitative comparison of FOM, PPE-ROM, and SUP 2 -ROM, we consider the quantities of interest for this benchmark, i.e. the drag and lift coefficients [28,49]: where U r = 1 is the maximum velocity at the inlet/outlet, L r = 0.1 is the cylinder diameter, S is the cylinder surface, and t and n are the tangential and normal unit vectors to the cylinder, respectively. See Fig. 4 for the coefficients in (71) computed by the three approaches. We observe that the amplitude of both coefficients is slightly underestimated (resp., overestimated) by the PPE-ROM (resp., SUP 2 -ROM) over the entire time interval. The ROM reconstruction of the lift coefficient appears to be more critical, especially for the SUP 2 -ROM and around the center of the time interval. For a further quantitative assessment of the reconstruction of the coefficients in (71), we computed the following errors We notice that errors defined in (72) are different from the ones in [23], which are related to the maximum values only. Table 3 reports the errors (72) for PPE-ROM and SUP 2 -ROM. We see that the two ROM strategies provide comparable results: slightly lower than 9% relative error for the the drag coefficient and around 14% relative error for the lift coefficient.
E c d E c l PPE-ROM 0.088 0.144 SUP 2 -ROM 0.086 0.139 Table 3: Relative errors (72) for lift and drag coefficients for PPE-ROM and SUP 2 -ROM.
Finally, we provide some information on the computational cost. The total CPU time required by a FOM simulation is about 1900 s. The solution of the PPE-ROM algebraic systems (38), (52), (44), and (58) takes 1.61 s, while the solution of the SUP 2 -ROM systems (38), (39), (44), and (45) takes 2.56 s. So, the SUP 2 -ROM is less efficient than the PPE-ROM. The additional cost of the SUP 2 -ROM comes from the supremizer modes, which increases the size of the reduced dynamical system. However, it is possible to assert that both ROMs allow to obtain a considerable speed-up.

Conclusions and future perspectives
We presented a POD-Galerkin based reduced order method for a Leray model implemented through the Evolve-Filter (EF) algorithm. Unlike the large majority of the works on Leraytype models, we choose a Finite Volume method for the space discretizaion because of its computational efficiency. The novelty of this work is the investigation and comparison of two techniques for the stabilization of the pressure fields: (i) Poisson Pressure equation, and (ii) supremizer enrichment. We showed that the standard supremizer enrichment, which works well for the Navier-Stokes equations with no filter, needs to be modified in order to obtain stable and accurate solutions with the EF algorithm. The modification consists in adding to the evolve and filter velocity spaces the supremizer solutions related to both evolve and filter pressure fields. We assessed our ROM through the classical 2D flow past a cylinder benchmark. We found that our ROM with both Poisson Pressure equation and modified supremizer enrichment captures the flow features with an accuracy comparable to ROMs applied to the Navier-Stokes equations with no filter [46]. Moreover, we quantified the relative error of the drag and lift coefficients computed by ROM and FOM and found that both stabilization approaches produce comparable errors.
In the future, we would like to investigate in more depth the inf-sup stability of a ROM formulation with supremizer enrichment for the EF algorithm. In addition, we would like to extend to the EF algorithm other efficient stabilization techniques, such as the one proposed in [47].

Disclosure statement
No potential conflict of interest was reported by the author(s).