Least-squares pressure recovery in reduced order methods for incompressible flows

Journal of Computational Physics 519 (2024) 113397 Available online 17 September 2024 0021-9991/© 2024 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp Least-squares pressure recovery in reduced order methods for incompressible flows M. Azaïeza,∗, T. Chacón Rebollo b, M. Oulghelouc, I. Sánchez Muñoz d aBordeaux University, Bordeaux INP and I2M (UMR CNRS 5295), 33400 Talence, France bDepartamento EDAN & IMUS, Universidad de Sevilla, Spain cInstitut Jean Le rond d’Alembert, Sorbonne Universite, 4 Place Jussieu, Paris, France dDepartamento Matemática Aplicada I, Universidad de Sevilla, Spain A R T I C L E I N F O A B S T R A C T Keywords: Pressure recovery Inf-sup condition Reduced order methods Incompressible flow Navier-Stokes equations In this work, we introduce a method to recover the reduced pressure for Reduced Order Models (ROMs) of incompressible flows. The pressure is obtained as the least-squares minimum of the residual of the reduced velocity with respect to a dual norm. We prove that this procedure provides a unique solution whenever the full-order pair of velocity-pressure spaces is inf-sup stable. We also prove that the proposed method is equivalent to solving the reduced mixed problem with reduced velocity basis enriched with the supremizers of the reduced pressure gradients. Optimal error estimates for the reduced pressure are obtained for general incompressible flow equations and specifically, for the transient Navier-Stokes equations. We also perform some numerical tests for the flow past a cylinder and the lid-driven cavity flow which confirm the theoretical expectations, and show an improved convergence with respect to other pressure recovery methods. 1. Introduction Apart from addressing the non-linear convective term, the primary challenge in numerically solving the time-dependent, incompressible Navier-Stokes equations lies in the Stokes stage. Specifically, determining the pressure field acts as a Lagrange multiplier to ensure the continuity equation. This issue has led to a significant body of work and publications in both non-reduced and reducedorder methods. In the non-reduced framework, two main families of approaches can be considered: one based on time splitting and the projection method. For comprehensive coverage on this topic, refer to the following references [9,19]. The second family of methods operates without splitting and utilizes techniques such as Uzawa algorithms. For further insights, consult the following references [8,20,24]. The challenge of pressure calculation persists in the numerical solution of the Navier-Stokes equations using reduced bases. According to our review of the literature, solutions to this issue depend on the nature of the vectors in the bases of the reduced spaces, specifically on whether these vectors are divergence-free or not. For comments and justifications on this issue, refer, for example, to [3,7,13,16]and the references therein. Indeed, if the reduced base only weakly satisfies the constraint, it becomes impossible to directly calculate the pressure, which is precisely the Lagrange multiplier in the equations ensuring the constraint, using the Galerkin-POD-ROM application. * Corresponding author. E-mail addresses: [email protected] (M. Azaïez), [email protected] (T. Chacón Rebollo), [email protected] (M. Oulghelou), [email protected] (I. Sánchez Muñoz). https://doi.org/10.1016/j.jcp.2024.113397 Received 8 February 2024; Received in revised form 15 July 2024; Accepted 4 September 2024 Journal of Computational Physics 519 (2024) 113397 2 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. To overcome this challenge, various methods have been proposed to recover the pressure of the reduced model a posteriori. These include: i) solving the Poisson pressure equation (PPE) with appropriate boundary conditions. In this formula, the continuity equation of the reduced system is replaced by the Poisson pressure equation, which is obtained by taking the divergence of the momentum equation and exploiting the fact that the velocity field is divergence-free. For the Poisson equation to be consistent with the initial NSE, we need to impose appropriate boundary conditions. Some possibilities include enforcing a no-slip boundary on the velocity, or incorporating a Neumann boundary condition, which in some cases leads to an artificial boundary layer. Full details of these different approaches can be found in [2,11,15,16,23]. More recently, in [5], the authors proposed to solve a pressure equation by duality with the gradients of the reduced pressure basis functions. This amounts to solving a Poisson equation for the pressure, with the advantage of defining better adapted boundary conditions and ensuring better precision in their recovery. However, their approach does not seem to provide a clear answer for the case of non-ordinary boundary conditions, such as the non-stress or outflow boundary condition. Finally, it is important to note another inconsistency in the choice of the 𝐻1regularity of the pressure to solve the Poisson problem compared to that of NS, where the pressure is only 𝐿2. ii) solving the pressure via the momentum equation recovery formulation ∇𝑝𝑟=𝐟−𝜕𝐭𝒖𝐫…. In practice, additional velocity space modes, called supremizer modes, are introduced to compute the reduced pressure. They are chosen to satisfy the inf-sup condition [3,4,10,22]. This approach to pressure recovery was introduced in the context of the POD for parameterised NSEs (see [3]). A stability analysis in terms of the existence of an inf-sup condition was carried out in the same reference. The method was later extended in [7]to the case where a strongly divergence-free POD velocity basis is used. It has also been used in the context of Petrov-Galerkin methods in [1,6]. The way to ensure a stable enrichment of the reduced velocity basis is to go through each element of the reduced pressure basis and compute the Riesz representation of the linear, continuous functional associated with the gradient of it. It is obtained by solving a vectorial Poisson problem with as right hand side the gradient of the pressure element. More recently, in [12], data-driven closure/correction terms are given and studied to increase the pressure and velocity accuracy of reduced order models (ROMs) for fluid flows. These references propose some supremiser stabilisation techniques to ensure compatibility between the pressure and added velocity spaces. In [13], the authors proposed a numerical analysis of the two approaches. They proved stability and convergence results for the supremizer-stabilised approach and numerically investigated these properties, in addition to its performance against the pressure Poisson method. In this paper, we propose an approach falling into category (ii). We adopt a strategy that has the dual advantage of calculating pressure without resorting to additional boundary conditions and extending the choice of velocities in the supremizers procedure. We propose a recovery procedure using least-squares minimization of the residual dual norm. We prove that this procedure provides a unique solution whenever the full-order pair of velocity-pressure spaces is inf-sup stable. Additionally, we establish that this leastsquares recovery of the pressure is equivalent to solving the reduced mixed problem with the reduced velocity basis enriched with the supremizers of the reduced pressure gradients. In essence, the pressure computed by the supremizers procedure is treated as the least-squares pressure provided by the velocity solution of the reduced mixed problem. A similar result holds for the full-order mixed problem. As a consequence, the reduced pressure obtained by the least-squares method is equivalent to the pressure obtained from MER formulation with supremizers as test functions, which has been studied in [7,13]. We also prove error estimates for the least-squares recovered pressure for general incompressible flow equations. Finally, we present numerical results for transient Navier- Stokes equations, demonstrating good agreement with theoretical expectations. We obtain discretization errors in the approximation of the reduced pressure quite close to the POD projection, achieving error reductions of four orders of magnitude compared to the pressure Poisson Equation approach. The following sections of this paper are organized as follows: in Section 2, we introduce the parametrized incompressible Navier- Stokes equation and derive its reduced-order Galerkin formulation. Moving on to Section 3, we describe the reduced pressure recovery through the least-squares minimum residual, analyzing the existence and uniqueness of the solution. This section concludes with an algebraic description of the operator on the pressure. Section 4establishes the equivalence between the pressure gradient supremizers procedure and the least-squares recovery of the pressure, along with exploring the adoption of a new choice of supremizers procedure. Optimal error estimates for the pressure, in terms of the approximation, are proven in Section 5. In Section 6, we conduct a numerical investigation into the performance of these pressure recovery techniques. Numerical tests are performed to assess and quantify the error estimates derived in the previous section for the pressure, and a comparison is made with the results obtained by the method introduced by Chacon et al. (refer to [5]). 2. The need of pressure recovery for ROMs of incompressible flows We start from a ROM for an incompressible flow that is solved only for velocity, dropping the free divergence equation. To describe it, we consider a boundary value parametric model for the incompressible Navier-Stokes-Brinkman. Set Ω ⊂ℝ𝑑(𝑑=2or 3) be a bounded domain with a Lipschitz-continuous boundary Γ =𝜕Ωwhich is splitted into two disjoint parts Γ =Γ 𝐷∪Γ 𝑁. We consider the problem: Journal of Computational Physics 519 (2024) 113397 3 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Find a velocity field 𝒖(𝜇) ∶Ω →ℝ𝑑and a pressure 𝑝(𝜇) ∶Ω →ℝsuch that ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝛾𝒖(𝜇)+𝒖(𝜇)⋅∇𝒖(𝜇)−𝜇Δ𝒖(𝜇)+∇𝑝(𝜇)=𝒇in Ω, ∇⋅𝒖(𝜇)=0 in Ω, 𝒖(𝜇)=𝟎on Γ𝐷, −𝜇𝜕 𝒏𝒖(𝜇)+𝑝(𝜇)𝒏=𝟎on Γ𝑁, (1) where 𝛾is a positive constant, 𝒇∶Ω →ℝ𝑑is a given external body force field per unit mass and 𝒏is the external normal unit vector on Γ. Here 𝜇=1 𝑅𝑒 , where 𝑅𝑒 is the Reynolds number, is considered as a parameter. In order to present the basic ideas of the pressure recovery method that we are going to introduce in this paper, we will consider a steady-state problem that depends on a single parameter with homogeneous boundary conditions. It can easily be extended to the transient problem depending on different parameters, either physical or geometrical or both, and more general boundary conditions. In fact, the first term in the first equation of the model includes the case of the evolution problem, once it has been discretised in time. In order to formulate problem (1)in weak form we introduce for the velocity the space 𝐇1 𝐷(Ω) = {𝒗∈𝐻1(Ω)𝑑such that 𝒗=𝟎on Γ𝐷}, 𝐇−1(Ω) its dual space, while for the pressure we consider the space 𝑀that can be 𝐿2(Ω) or 𝐿2 0(Ω) ={𝑞∈𝐿2(Ω) such that ∫Ω𝑞𝑑𝒙=0 } when Γ𝑁=∅. We denote by (⋅, ⋅)0the 𝐿2inner product, either for scalar or vector functions, (⋅, ⋅)1the inner product in 𝐇1 𝐷(Ω). We also denote by ‖ ⋅‖0and ‖ ⋅‖1the norms defined by these inner products, respectively, and ⟨⋅, ⋅⟩the duality pairing between 𝐇−1(Ω) and 𝐇1 𝐷(Ω). The weak formulation of problem (1) consists in: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Find 𝒖(𝜇)∈𝐇1 𝐷(Ω) and 𝑝(𝜇)∈𝐿2 0(Ω) such that 𝑎(𝒖(𝜇),𝒖(𝜇),𝒗;𝜇)−(∇⋅𝒗,𝑝(𝜇))0=⟨𝒇,𝒗⟩,∀𝒗∈𝐇1 𝐷(Ω), (∇ ⋅𝒖(𝜇),𝑞)0=0 ∀𝑞∈𝐿2 0(Ω), (2) where 𝑎is the trilinear form given by 𝑎(𝒘,𝒖,𝒗;𝜇)=𝛾(𝒖,𝒗)0+(𝒘⋅∇𝒖,𝒗)0+𝜇(∇𝒖,∇𝒗)0,∀𝒘,𝒖,𝒗∈𝐇1 𝐷(Ω). To state the Galerkin discretisation of problem (2), let us consider a family of inf-sup stable discrete velocity-pressure spaces (𝑿ℎ, 𝑀ℎ) ⊂𝐇1 𝐷(Ω) ×𝐿2 0(Ω). That is, there exists a constant 𝛽>0such that 𝛽‖𝑞ℎ‖0≤sup 𝒗ℎ∈𝑿ℎ (∇ ⋅𝒗ℎ,𝑞ℎ)0 ‖∇𝒗ℎ‖0 ,∀𝑞ℎ∈𝑀ℎ,∀ℎ>0.(3) Note that for some discretisation techniques, in particular for spectral methods, the constant 𝛽may depend on the discretisation parameter ℎ. We consider for instance the Galerkin discretisation of problem (2) (the “Full Order Model”, FOM), that consists in: ⎧ ⎪ ⎨ ⎪ ⎩ Find 𝒖ℎ(𝜇)∈𝑿ℎand 𝑝ℎ(𝜇)∈𝑀ℎsuch that 𝑎(𝒖ℎ(𝜇),𝒖ℎ(𝜇),𝒗ℎ;𝜇)−(∇⋅𝒗ℎ,𝑝ℎ(𝜇))0=⟨𝒇ℎ,𝒗ℎ⟩,∀𝒗ℎ∈𝑿ℎ, (∇ ⋅𝒖ℎ(𝜇),𝑞 ℎ)0=0 ∀𝑞ℎ∈𝑀ℎ, (4) where 𝒇ℎis a given data coming from 𝒇. For example, for the transient case, approximating the time evolution by a first-order Euler scheme then 𝒇ℎ=𝒇+1 Δ𝑡𝒖𝑛−1 ℎ. To state the ROM approximation of problem (2), let us assume that we already have determined the reduced basis spaces for the discrete velocity and pressure (𝑿𝑟, 𝑀𝑟), such that: 𝑿𝑟⊂𝑿0ℎ={𝒗ℎ∈𝑿ℎsuch that (∇ ⋅𝒗ℎ,𝑞ℎ)=0,∀𝑞ℎ∈𝑀ℎ}.(5) We then consider the Galerkin projection of the momentum equation, as ROM for the velocity: {Find 𝒖𝑟(𝜇)∈𝑿𝑟such that 𝑎(𝒖𝑟(𝜇),𝒖𝑟(𝜇),𝒗𝑟;𝜇)=⟨𝒇𝑟,𝒗𝑟⟩,∀𝒗𝑟∈𝑿𝑟,(6) where 𝒇𝑟is a reduced data coming from 𝒇. Note that the reduced pressure 𝑝𝑟(𝜇) ∈𝑀𝑟verifies (∇ ⋅𝒗ℎ, 𝑝𝑟(𝜇)) =0, ∀ 𝒗ℎ∈𝑿0ℎ, because 𝑀𝑟⊂𝑀 ℎ, and thus the reduced pressure is removed from problem (6). However, this problem can be solved by itself, and only the reduced velocity is computed (as is done in practice for many engineering problems). Journal of Computational Physics 519 (2024) 113397 4 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Nevertheless, for many applications in aerodynamics, hemodynamics and hydrodynamics, knowing the pressure is of paramount interest, to compute for instance drag and lift coefficients, or the pressure exerted by an artery on its walls. We will therefore tackle the problem of recovering the pressure from the velocity solution 𝒖𝑟(𝜇)of the reduced problem (6). 3. Reduced pressure recovery by least-squares minimum residual In the following, we omit the dependency of the reduced solution on the parameters for brevity. The keystone of our method is to recover the reduced pressure by means of least-squares method. In order to describe it, let us introduce the following discrete operators: 𝐺∶𝑀ℎ↦𝐇−1(Ω) given by ⟨𝐺𝑞 ℎ,𝒗ℎ⟩=−(∇⋅𝒗ℎ,𝑞ℎ)0,∀𝑞ℎ∈𝑀ℎ,∀𝒗ℎ∈𝑿ℎ.(7) Note that ⟨𝐺𝑞 ℎ, 𝒗ℎ⟩ =⟨∇𝑞ℎ, 𝒗ℎ⟩ −∫ Γ𝑁 𝑞ℎ𝒗ℎ⋅𝒏and thus, if we consider Dirichlet boundary conditions on all the boundary Γ, the operator 𝐺would be the gradient operator. Π(𝑘) ℎ∶𝐇−1(Ω) ↦𝑿ℎdefined by (Π(𝑘) ℎ𝜑,𝒗ℎ)𝑘=⟨𝜑,𝒗ℎ⟩,∀𝜑∈𝐇−1(Ω),∀𝒗ℎ∈𝑿ℎ,(8) for 𝑘 =0or 𝑘 =1. Observe that from (8)and (7), (Π(𝑘) ℎ(𝐺𝑞 ℎ),𝒗ℎ)𝑘=−(∇⋅𝒗ℎ,𝑞ℎ)0,∀𝑞ℎ∈𝑀ℎ,∀𝒗ℎ∈𝑿ℎ.(9) We then propose to recover the reduced pressure as a solution of the following least-squares minimum residual problem, either with respect to the 𝐋2(Ω) or the 𝐇1 𝐷(Ω) norms, 𝑝𝑟=argmin 𝑞𝑟∈𝑀𝑟 𝐽(𝑞𝑟)∶=‖Π(𝑘) ℎ(𝐺𝑞 𝑟−𝑅(𝒖𝑟))‖2 𝑘,(10) where 𝑅(𝒖𝑟) ∈𝐇−1(Ω) given by ⟨𝑅(𝒖𝑟),𝒗⟩=⟨𝒇𝑟,𝒗⟩−𝑎(𝒖𝑟,𝒖𝑟,𝒗;𝜇),∀𝒗∈𝐇1 𝐷(Ω).(11) The optimality conditions of this problem are (Π(𝑘) ℎ(𝐺𝑝 𝑟),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘=(Π (𝑘) ℎ𝑅(𝒖𝑟),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘,∀𝑞𝑟∈𝑀𝑟.(12) The following results state that the method is well-defined. Lemma 1. Assume that the discrete inf-sup condition (3)holds. Then, for 𝑘 =0and 𝑘 =1, ‖𝑞ℎ‖ℎ,𝑘 =‖Π(𝑘) ℎ(𝐺𝑞 ℎ)‖𝑘,∀𝑞ℎ∈𝑀ℎ,(13) is a norm on 𝑀ℎwhich satisfies ‖𝑞ℎ‖0≤𝛼‖𝑞ℎ‖ℎ,𝑘,∀𝑞ℎ∈𝑀ℎ,(14) for some 𝛼>0. Moreover, for 𝑘 =1the norm (13)is uniformly equivalent to the 𝐿2(Ω) norm on 𝑀ℎ. Proof. It is straightforward to prove that ‖ ⋅‖ℎ,𝑘 is a semi-norm on 𝑀ℎ. To prove (14), let 𝑞ℎ∈𝑀ℎ, 𝑞ℎ≠0. Observe that from (9), the inf-sup condition (3) can be written as 𝛽‖𝑞ℎ‖0≤sup 𝒗ℎ∈𝑿ℎ (∇ ⋅𝒗ℎ,𝑞ℎ)0 ‖∇𝒗ℎ‖0 =sup 𝒗ℎ∈𝑿ℎ (Π(𝑘) ℎ(𝐺𝑞 ℎ),𝒗ℎ)𝑘 ‖∇𝒗ℎ‖0 =sup 𝒗ℎ∈𝑿ℎ ‖∇𝒗ℎ‖0=1 (Π(𝑘) ℎ(𝐺𝑞 ℎ),𝒗ℎ)𝑘=‖Π(𝑘) ℎ(𝐺𝑞 ℎ)‖2 𝑘 ‖∇Π(𝑘) ℎ(𝐺𝑞 ℎ)‖0 , (15) where in the last equality we use that the maximun of the inner product (Π(𝑘) ℎ(𝐺𝑞 ℎ), 𝒗ℎ)𝑘is reached when 𝒗ℎ=𝜆 Π(𝑘) ℎ(𝐺𝑞 ℎ)for any 𝜆 >0. Then, we obtain the equality with  𝒗ℎ=1 ‖∇Π(𝑘) ℎ(𝐺𝑞 ℎ)‖0 Π(𝑘) ℎ(𝐺𝑞 ℎ)such that ‖∇ 𝒗ℎ‖0=1. Journal of Computational Physics 519 (2024) 113397 5 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. When 𝑘 =0, ‖Π(𝑘) ℎ(𝐺𝑞 ℎ)‖2 𝑘 ‖∇Π(𝑘) ℎ(𝐺𝑞 ℎ)‖0 ≤𝑐𝑃‖Π(𝑘) ℎ(𝐺𝑞 ℎ)‖0=𝑐𝑃‖𝑞ℎ‖ℎ,0,(16) where 𝑐𝑃is the Poincaré’s constant yielding the inequality ‖𝒗‖0≤𝑐𝑃‖∇𝒗‖0,∀𝒗∈𝐇1 𝐷(Ω). When 𝑘 =1, ‖Π(𝑘) ℎ(𝐺𝑞 ℎ)‖2 1 ‖∇Π(𝑘) ℎ(𝐺𝑞 ℎ)‖0 =‖Π(𝑘) ℎ(𝐺𝑞 ℎ)‖1=‖𝑞ℎ‖ℎ,1.(17) Inequality (14) follows from (15), (16)and (17)with 𝛼=𝑀 𝛽being 𝑀=𝑐𝑃for 𝑘 =0or 𝑀=1for 𝑘 =1. Moreover, from (15)and (17), ‖𝑞ℎ‖ℎ,1=sup 𝒗ℎ∈𝑿ℎ (∇ ⋅𝒗ℎ,𝑞ℎ) ‖∇𝒗ℎ‖0 ≤√𝑑‖𝑞ℎ‖0, which, combined with (14)gives the equivalence of the norms ‖ ⋅‖0and ‖ ⋅‖ℎ,1on 𝑀ℎ.□ Remark 1. Actually, the equivalence of the norms ‖ ⋅‖0and ‖ ⋅‖ℎ,1on 𝑀ℎis a re-formulation of the discrete inf-sup property, that states that the 𝐿2(Ω) norm of the discrete pressures and the 𝐇−1(Ω) norm of their gradients are equivalent norms in 𝑀ℎ. Theorem 3.1. Assume that the discrete inf-sup condition (3)holds. Then problem (10)admits a unique solution, which is characterised as the solution of the normal equations (12). Proof. Let us consider the space 𝑺𝑟={Π (𝑘) ℎ(𝐺𝑞 𝑟)such that 𝑞𝑟∈𝑀𝑟}⊂𝑿ℎ.(18) Due to the inf-sup condition, the mapping =Π (𝑘) ℎ◦𝐺∶𝑀𝑟↦𝑺𝑟is bijective. Indeed, it is trivially surjective and from Lemma 1, the inf-sup condition implies that if Π(𝑘) ℎ(𝐺𝑞 𝑟) =0, then 𝑞𝑟=0. Therefore, the spaces 𝑺𝑟and 𝑀𝑟are isomorphic. Problem (10) can equivalently be expressed as 𝑝𝑟=−1(𝒔𝑟),with 𝒔𝑟=argmin 𝝃𝑟∈𝑺𝑟  𝐽(𝝃𝑟)∶=‖𝝃𝑟−Π (𝑘) ℎ(𝑅(𝒖𝑟))‖2 𝑘.(19) Due to the standard theory of least-squares approximation on finite-dimension sub-spaces, this problem admits a unique solution, characterised by the normal equations, that is (𝒔𝑟,𝝃𝑟)𝑘=(Π (𝑘) ℎ(𝑅(𝒖𝑟)),𝝃𝑟)𝑘,∀𝝃𝑟∈𝑺𝑟. These equations indeed are (12). □ 3.1. Matrix expression of pressure recovery problem Let {𝜙𝑖}𝑛ℎ 𝑖=1 be a basis of 𝑿ℎand {𝜓𝑖}𝑛𝑟 𝑖=1 a basis of 𝑀𝑟, then problem (12)is equivalent to Find 𝑝𝑟= 𝑛𝑟 ∑ 𝑖=1 𝑝𝑖𝜓𝑖∈𝑀𝑟such that 𝑝 ∈ℝ𝑛𝑟with (𝑝)𝑖=𝑝𝑖is the solution of the linear system: 𝑝 =𝑟 with =−1 𝑡and 𝑟 =−1  𝑅, (20) where ∈ℝ𝑛𝑟×𝑛ℎ∶𝑖𝑗 =−(∇⋅𝜙𝑗,𝜓𝑖)0, ∈ℝ𝑛ℎ×𝑛ℎ∶𝑖𝑗 =(𝜙𝑗,𝜙𝑖)𝑘  𝑅∈ℝ𝑛ℎ∶ 𝑅𝑖=⟨𝑅(𝒖𝑟),𝜙 𝑖⟩. Indeed, for all 𝑖 =1 … , 𝑛𝑟, Π(𝑘) ℎ(𝐺𝜓 𝑖)= 𝑛ℎ ∑ 𝑗=1 𝛼𝑖 𝑗𝜙𝑗,where 𝛼𝑖=−1  𝑏𝑖with ( 𝑏𝑖)𝑗=−(∇⋅𝜙𝑗,𝜓𝑖). Journal of Computational Physics 519 (2024) 113397 6 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Then, 𝑖𝑗 =(Π (𝑘) ℎ(𝐺𝜓 𝑗),Π(𝑘) ℎ(𝐺𝜓 𝑖))𝑘=(  𝑏𝑖)𝑡(−1)𝑡 𝑏𝑗. Note that  ∈ℝ𝑛𝑟×𝑛𝑟and thus, (20)is a low-dimensional system. Moreover, if we consider the Cholesky factorization of −1, −1 =𝑡, system (20) can be written as 𝑡𝑝 =𝑡 𝑅, with =.(21) In this case, 𝑝 ∈ℝ𝑛𝑟is the solution of problem 𝑝 =argmin 𝑞∈ℝ𝑛𝑟‖𝑡𝑞 −𝑡 𝑅‖2 ℝ𝑛𝑟.(22) That is, 𝑝 is the least-squares solution of the system 𝑡𝑝 =𝑡 𝑅. In particular, when the basis of 𝑿ℎis ortonormal with respect to the inner product (⋅, ⋅)𝑘, matrix  =𝐼𝑑 and 𝑝 is the least-squares solution of the system 𝑡𝑝 = 𝑅. From a purely algebraic point of view, the inf-sup condition (3)implies that system (20)admits a unique solution. This is stated as follows Proposition 1. Assume that the inf-sup condition (3)holds. Then the matrix of the linear system (20)is positive definite and thus, this system admits a unique solution. Proof. For simplicity of calculation, let us assume that the basis {𝜙𝑖}𝑛ℎ 𝑖=1 is orthogonal in 𝐿2(Ω). Then  =𝐼𝑑 and  = 𝑡. For any 𝑞 ∈ℝ𝑛𝑟, 𝑞𝑡𝑞 = 𝑛𝑟 ∑ 𝑖,𝑗=1 𝑞𝑖(𝑛ℎ ∑ 𝑘=1 (∇ ⋅𝜙𝑘,𝜓𝑖)0(∇ ⋅𝜙𝑘,𝜓𝑗)0)𝑞𝑗 = 𝑛ℎ ∑ 𝑘=1 (∇ ⋅𝜙𝑘, 𝑛𝑟 ∑ 𝑖=1 𝑞𝑖𝜓𝑖)0(∇ ⋅𝜙𝑘, 𝑛𝑟 ∑ 𝑗=1 𝑞𝑗𝜓𝑗)0= 𝑛ℎ ∑ 𝑘=1 (∇ ⋅𝜙𝑘,𝑞𝑟)2 0, with 𝑞𝑟= 𝑛𝑟 ∑ 𝑖=1 𝑞𝑖𝜓𝑖∈𝑀𝑟. Then 𝑞𝑡𝑞 ≥0, for all 𝑞 ∈ℝ𝑛𝑟. Moreover, 𝑞𝑡𝑞 =0if and only if (∇ ⋅𝜙𝑘, 𝑞𝑟)0, for all 𝑘 =1 … 𝑛ℎand thus, (∇ ⋅𝒗ℎ, 𝑞𝑟)0=0for all 𝒗ℎ∈𝑿ℎ. From (3), this implies that 𝑞𝑟=0and thus 𝑞 =0.□ 4. The pressure gradient supremisers revisited It is a common practice in ROM of incompressible flows to solve the reduced counterpart of the mixed problem (2)with enriched velocity spaces 𝑿𝑟in such a way that the couple (𝑿𝑟, 𝑀𝑟)satisfies an inf-sup condition. The reduced velocity space actually is built as 𝑿𝑟=𝑺𝑟⊕𝑿0𝑟,(23) where 𝑺𝑟is defined by (18)and 𝑿0𝑟is a reduced space formed by weakly divergence-free velocities, that is a subset of 𝑿0ℎgiven by (5). Observe that the elements of 𝑺𝑟are the pressure-gradient supremisers introduced in [21]and applied in several works to build reduced order models for incompressible flows, for instance [3,17,18,22]. In this way, the reduced problem: ⎧ ⎪ ⎨ ⎪ ⎩ Find 𝒖𝑟(𝜇)∈𝑿𝑟and 𝑝𝑟(𝜇)∈𝑀𝑟such that 𝑎(𝒖𝑟(𝜇),𝒖𝑟(𝜇),𝒗𝑟;𝜇)−(∇⋅𝒗𝑟,𝑝𝑟(𝜇))0=⟨𝒇𝑟,𝒗𝑟⟩,∀𝒗𝑟∈𝑿𝑟, (∇ ⋅𝒖𝑟(𝜇),𝑞𝑟)0=0,∀𝑞𝑟∈𝑀𝑟, (24) is well posed. The equivalence between this pressure gradient supremizers procedure and the least-squares recovery of the pressure introduced above is stated as follows. Proposition 2. Assume that (23)holds. Then the velocity 𝒖𝑟solution of problem (24)satisfies {Find 𝒖𝑟(𝜇)∈𝑿0𝑟such that 𝑎(𝒖𝑟(𝜇),𝒖𝑟(𝜇),𝒗0𝑟;𝜇)=⟨𝒇𝑟,𝒗0𝑟⟩,∀𝒗0𝑟∈𝑿0𝑟,(25) and the pressure 𝑝𝑟satisfies (12). Journal of Computational Physics 519 (2024) 113397 7 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Proof. From (23), the velocity 𝒖𝑟is decomposed into 𝒖𝑟=𝒔𝑟+𝒖0𝑟,with 𝒔𝑟∈𝑺𝑟and 𝒖0𝑟∈𝑿0𝑟. Let 𝜋𝑟=−1(𝒔𝑟) ∈𝑀𝑟, that is, 𝒔𝑟=Π (𝑘) ℎ(𝐺𝜋 𝑟). As 𝑀𝑟⊂𝑀 ℎ, (∇ ⋅𝒖0𝑟, 𝜋𝑟)0=0. From the second equation in (24), it follows 0=(∇⋅𝒖𝑟,𝜋 𝑟)0=(∇⋅𝒔𝑟,𝜋 𝑟)0=−(Π (𝑘) ℎ(𝐺𝜋 𝑟),𝒔𝑟)𝑘=−‖Π(𝑘) ℎ(𝐺𝜋 𝑟)‖2 𝑘=−‖𝜋𝑟‖2 ℎ,𝑘. Therefore, 𝜋𝑟=0and consequently 𝒔𝑟=0, what implies 𝒖𝑟=𝒖0𝑟∈𝑿0𝑟. Thus, the first equation in (24)gives (25). Moreover, the first equation in (24) can be written as −(𝑝𝑟,∇⋅𝒗𝑟)0=⟨𝑅(𝒖𝑟),𝒗𝑟⟩,∀𝒗𝑟∈𝑿𝑟, and, from the decomposition (23), this equation is equivalent to −(𝑝𝑟,∇⋅𝒔𝑟)0=⟨𝑅(𝒖𝑟),𝒔𝑟⟩,∀𝒔𝑟∈𝑺𝑟,(26) as (𝑝𝑟, ∇ ⋅𝒗0𝑟)0=0, for all 𝒗0𝑟∈𝑿0𝑟. Then, 𝑝𝑟satisfies (Π(𝑘) ℎ(𝐺𝑝 𝑟),𝒔𝑟)𝑘=(Π (𝑘) ℎ(𝑅(𝒖𝑟)),𝒔𝑟)𝑘, for all 𝒔𝑟=Π (𝑘) ℎ(𝐺𝑞𝑟) ∈𝑺𝑟with any 𝑞𝑟∈𝑀𝑟, that is, problem (12). □ Remark 2. From Proposition 2the pressure computed by the supremizers procedure is cast as the least-squares pressure provided by the velocity solution of the reduced mixed problem (24). In addition, the computation of this velocity can be partially decoupled from that of the pressure, as it solves problem (25). Therefore the (velocity, pressure) pair provided by the supremizers procedure can be computed by solving two coercive problems, sequentially in two steps. In the first step the velocity is computed by solving (25)and in the second step the pressure is computed by solving (12), that is: 𝒖𝑟∈𝑿0𝑟,𝑎(𝒖𝑟(𝜇),𝒖𝑟(𝜇),𝒗0𝑟;𝜇)=⟨𝒇𝑟,𝒗0𝑟⟩,∀𝒗0𝑟∈𝑿0𝑟; 𝑝𝑟∈𝑀𝑟,(Π(𝑘) ℎ(𝐺𝑝 𝑟),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘=(Π (𝑘) ℎ(𝑅(𝒖𝑟)),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘,∀𝑞𝑟∈𝑀𝑟. Remark 3. (26)is the formulation MER to recover the reduced pressure in ROMs (see [7,13]). Therefore, Proposition 2proves that the pressure obtained by this method is equivalent to the pressure given by least-squares method (10). Remark 4. All the preceding theory holds whenever the family of full-order velocity and reduced-order pressure spaces (𝑿ℎ, 𝑀𝑟) satisfies the discrete inf-sup condition. Then, 𝑿ℎneeds not to be the full finite element space appearing in the FOM problem (4), but it can be replaced by any smaller space that ensures the inf-sup condition. For instance, it can be a reduced space containing non-weakly divergence free velocities, or a finite element space in a coarse grid. This would decrease the time needed to compute in the off-line step the matrix −1 appearing in (20), without increasing the error in the computation of 𝑝𝑟(see Theorem 5.1 below). Remark 5. The previous analysis leads to casting general mixed methods in terms of pressure gradient supremisers. More specifically, if a pair of discrete velocity-pressure spaces (𝑿ℎ, 𝑀ℎ) ⊂𝐇1 𝐷(Ω) ×𝐿2 0(Ω) satisfies the inf-sup condition (3)then 𝑿ℎcan be expressed as 𝑿ℎ=𝑺ℎ⊕𝑿0ℎ,with 𝑺ℎ={Π (𝑘) ℎ(𝐺𝑞 ℎ)such that 𝑞ℎ∈𝑀ℎ},(27) for either 𝑘 =0or 𝑘 =1. Indeed, for any 𝒙ℎ∈𝑿ℎ, because of the inf-sup condition, there exists a unique (𝒙0ℎ, 𝜋ℎ) ∈𝑿ℎ×𝑀ℎsolution of problem {(𝒙0ℎ,𝒗ℎ)𝑘−(∇⋅𝒗ℎ,𝜋ℎ)0=(𝒙ℎ,𝒗ℎ)𝑘,∀𝒗ℎ∈𝑿ℎ, (∇ ⋅𝒙0ℎ,𝑞ℎ)0=0,∀𝑞ℎ∈𝑀ℎ, for 𝑘 =0or 𝑘 =1. Then, 𝒔ℎ=𝒙ℎ−𝒙0ℎverifies (𝒔ℎ,𝒗ℎ)𝑘=−(∇⋅𝒗ℎ,𝜋ℎ)0=(Π (𝑘) ℎ(𝐺𝜋 ℎ),𝒗ℎ)𝑘,∀𝒗ℎ∈𝑿ℎ. Therefore, 𝒔ℎ=Π (𝑘) ℎ(𝐺𝜋 ℎ) ∈𝑺ℎ. The sum is direct because if 𝒙ℎ∈𝑺ℎ∩𝑿0ℎthen 𝒙ℎ=0as occurs to 𝒔𝑟in Proposition 2. Moreover, 𝜋ℎis the solution of problem (Π(𝑘) ℎ(𝐺𝜋 ℎ),Π(𝑘) ℎ(𝐺𝑞 ℎ))𝑘=(Π (𝑘) ℎ(𝒙ℎ−𝒙0ℎ),Π(𝑘) ℎ(𝐺𝑞 ℎ))𝑘,∀𝑞ℎ∈𝑀ℎ, and thus, 𝜋ℎis the least-squares pressure recovered by (10)on 𝑀ℎinstead 𝑀𝑟with 𝑅(𝒙ℎ) =𝒙ℎ−𝒙0ℎ. Journal of Computational Physics 519 (2024) 113397 8 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. 5. Error estimates In this section we obtain optimal error estimates for the pressure in terms of the approximation properties of space 𝑀𝑟and the distance between the FOM and ROM velocities. For that, let us assume that an inverse inequality between the 𝐋2(Ω) and 𝐇1 0(Ω) semi-norms holds on space 𝑿ℎ: ‖∇𝒗ℎ‖0≤𝐶ℎ‖𝒗ℎ‖0,∀𝒗ℎ∈𝑿ℎ,(28) where 𝐶ℎ>0is unbounded as the discretisation parameter ℎdecreases to zero. This constant depends on the actual kind of discretisation used. For instance, for finite element discretisations this constant scales as ℎ−1 min, where ℎmin is the smallest diameter of the grid elements. In sequel, we will use the following notation: given a Banach space 𝑋endowed with the norm ‖ ⋅‖𝑋and a subset 𝑌⊂𝑋, we denote the distance from an element 𝑥 ∈𝑋to 𝑌by 𝑑𝑋(𝑥,𝑌 )=inf 𝜁∈𝑌‖𝑥−𝜁‖𝑋. Theorem 5.1. Assume that the discrete inf-sup condition (3)holds, then the solution of problem (12)satisfies the following error estimates: If 𝑘 =0, ‖𝑝ℎ−𝑝𝑟‖0≤𝐶ℎ(𝑑𝐿2(Ω)(𝑝ℎ,𝑀𝑟)+‖𝑅(𝒖ℎ)−𝑅(𝒖𝑟)‖𝐇−1(Ω)),(29) where 𝐶ℎis given by (28). If 𝑘 =1, ‖𝑝ℎ−𝑝𝑟‖0≤𝐶(𝑑𝐿2(Ω)(𝑝ℎ,𝑀𝑟)+‖𝑅(𝒖ℎ)−𝑅(𝒖𝑟)‖𝐇−1(Ω)),(30) where 𝐶is a positive constant independent of ℎ. Here (𝒖ℎ, 𝑝ℎ)is a solution of the FOM problem (4), 𝑅(𝒖𝑟)is defined by (11)and similarly, ⟨𝑅(𝒖ℎ), 𝒗⟩ =⟨𝒇ℎ, 𝒗⟩ −𝑎(𝒖ℎ, 𝒖ℎ, 𝒗; 𝜇), ∀𝒗∈𝐇1 𝐷(Ω). Proof. From (4)it holds (Π(𝑘) ℎ(𝐺𝑝 ℎ),𝒗ℎ)𝑘=(Π (𝑘) ℎ(𝑅(𝒖ℎ)),𝒗ℎ)𝑘,for any 𝒗ℎ∈𝑿ℎ. Let 𝑝𝑟be the 𝐿2(Ω) protection of 𝑝ℎon 𝑀𝑟. It holds (Π(𝑘) ℎ(𝐺𝑝𝑟),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘=(Π (𝑘) ℎ𝐺(𝑝𝑟−𝑝ℎ)+(Π (𝑘) ℎ(𝑅(𝒖ℎ)),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘,∀𝑞𝑟∈𝑀𝑟. Let us define the error 𝑒𝑟=𝑝𝑟−𝑝𝑟∈𝑀𝑟and subtract (12)from this last equality. We deduce the following equation for 𝑒𝑟: (Π(𝑘) ℎ(𝐺𝑒 𝑟),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘=(Π (𝑘) ℎ𝐺(𝑝𝑟−𝑝ℎ),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘 +(Π(𝑘) ℎ(𝑅(𝒖ℎ)−𝑅(𝒖𝑟)),Π(𝑘) ℎ(𝐺𝑞 𝑟))𝑘,∀𝑞𝑟∈𝑀𝑟. (31) Taking 𝑞𝑟=𝑒𝑟, ‖Π(𝑘) ℎ(𝐺𝑒 𝑟)‖2 𝑘 =(Π (𝑘) ℎ∇( 𝑝𝑟−𝑝ℎ),Π(𝑘) ℎ(𝐺𝑒 𝑟))𝑘+(Π (𝑘) ℎ(𝑅(𝒖ℎ)−𝑅(𝒖𝑟)),Π(𝑘) ℎ(𝐺𝑒 𝑟))𝑘. (32) We now bound the terms of the last line of (32): (Π(𝑘) ℎ∇( 𝑝𝑟−𝑝ℎ),Π(𝑘) ℎ(𝐺𝑒 𝑟))𝑘=−(∇⋅Π(𝑘) ℎ(𝐺𝑒 𝑟),𝑝𝑟−𝑝ℎ)0 ≤√𝑑‖∇Π(𝑘) ℎ(𝐺𝑒 𝑟)‖0‖𝑝𝑟−𝑝ℎ‖0. (33) (Π(𝑘) ℎ(𝑅(𝒖ℎ)−𝑅(𝒖𝑟)),Π(𝑘) ℎ(𝐺𝑒 𝑟))𝑘=⟨𝑅(𝒖ℎ)−𝑅(𝒖𝑟),Π(𝑘) ℎ(𝐺𝑒 𝑟)⟩ ≤‖𝑅(𝒖ℎ)−𝑅(𝒖𝑟)‖𝐇−1(Ω) ‖∇Π(𝑘) ℎ(𝐺𝑒 𝑟)‖0. (34) Then, from estimates (32), (33)and (34)we have: If 𝑘 =0, applying (28)to ‖∇ Π(𝑘) ℎ(𝐺𝑒 𝑟)‖0, ‖𝑒𝑟‖ℎ,0≤𝐶ℎ(√𝑑‖𝑝𝑟−𝑝ℎ‖0+‖𝑅(𝒖ℎ)−𝑅(𝒖𝑟)‖𝐇−1(Ω)),(35) Journal of Computational Physics 519 (2024) 113397 9 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. If 𝑘 =1, ‖𝑒𝑟‖ℎ,1≤√𝑑‖𝑝𝑟−𝑝ℎ‖0+‖𝑅(𝒖ℎ)−𝑅(𝒖𝑟)‖𝐇−1(Ω).(36) Estimate (29)and (30) follow from estimates (35), (36)and (14). □ 5.1. Error estimates for the evolutionary model In this section we apply the previous analysis to the evolutionary incompressible Navier-Stokes equation. That is, assuming that the flow takes place during a time interval [0, 𝑇], we consider the problem: Find a velocity field 𝒖(𝜇) ∶Ω ×(0, 𝑇) →ℝ𝑑and a pressure 𝑝(𝜇) ∶Ω ×(0, 𝑇) →ℝsuch that ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜕𝑡𝒖(𝜇)+𝒖(𝜇)⋅∇𝒖(𝜇)−𝜇Δ𝒖(𝜇)+∇𝑝(𝜇)=𝒇in Ω×(0,𝑇), ∇⋅𝒖(𝜇)=0 in Ω×(0,𝑇), 𝒖(𝜇)=𝟎on Γ𝐷×(0,𝑇), −𝜇∇𝒖(𝜇)⋅𝒏+𝑝(𝜇)𝒏=0 on Γ𝑁×(0,𝑇), 𝒖(𝜇,0) = 𝒖0,in Ω, (37) where 𝒖0is a initial field velocity given. Consider a uniform partition of the interval [0, 𝑇], {0 =𝑡0<𝑡 1<… 𝑡𝑁=𝑇}, with time-step size Δ𝑡 =𝑇∕𝑁. The time discretization of problem (37)by the Backward Euler scheme carries to the following family of stationary problems: Given the initialization 𝒖0(𝜇) =𝑢0, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Find 𝒖𝑛(𝜇)∈𝐇1 𝐷(Ω) and 𝑝𝑛(𝜇)∈𝐿2 0(Ω) such that 𝑎(𝒖𝑛(𝜇),𝒖𝑛(𝜇),𝒗;𝜇)−(∇⋅𝒗,𝑝𝑛(𝜇))0=⟨𝒇𝑛,𝒗⟩,∀𝒗∈𝐇1 𝐷(Ω), (∇ ⋅𝒖𝑛(𝜇),𝑞)0=0 ∀𝑞∈𝐿2 0(Ω), ∀𝑛=1,2,…𝑁. (38) For each time step, this problem fits into the formulation (2)with 𝛾=1 Δ𝑡and 𝒇𝑛=𝒇+1 Δ𝑡𝒖𝑛−1. Then, we obtain time-space approximations of a solution of problem (37) applying the FOM (4)to problem (38): (𝒖𝑛 ℎ(𝜇), 𝑝𝑛 ℎ(𝜇)) ∈𝑿ℎ×𝑀ℎsuch that {𝑎(𝒖𝑛 ℎ(𝜇),𝒖𝑛 ℎ(𝜇),𝒗ℎ;𝜇)−(∇⋅𝒗ℎ,𝑝𝑛 ℎ(𝜇))0=⟨𝒇𝑛−1 ℎ,𝒗ℎ⟩,∀𝒗ℎ∈𝑿ℎ, (∇ ⋅𝒖𝑛 ℎ(𝜇),𝑞 ℎ)0=0 ∀𝑞ℎ∈𝑀ℎ,(39) where 𝒇𝑛−1 ℎ=𝒇+1 Δ𝑡𝒖𝑛−1 ℎ. We also obtain ROM approximations: (𝒖𝑛 𝑟(𝜇), 𝑝𝑛 𝑟(𝜇)) ∈𝑿𝑟×𝑀𝑟, where the reduced velocity 𝒖𝑛 𝑟(𝜇)is computed by problem 𝑎(𝒖𝑛 𝑟(𝜇),𝒖𝑛 𝑟(𝜇),𝒗𝑟;𝜇)=⟨𝒇𝑛−1 𝑟,𝒗𝑟⟩,∀𝒗𝑟∈𝑿𝑟,(40) where 𝒇𝑛−1 𝑟=𝒇+1 Δ𝑡𝒖𝑛−1 𝑟, and the reduced pressure 𝑝𝑛 𝑟(𝜇)is recovered from the reduced velocity by the least-squares method (10). To obtain error estimates for the reduced pressure we introduce the following discrete functions. In the following we omit the dependency on the parameters for brevity. •𝒖ℎ∶[0, 𝑇] →𝑿ℎis the piecewise linear in time function such that 𝒖ℎ(𝑡𝑛) =𝒖𝑛 ℎ. •𝑝ℎ∶[0, 𝑇] →𝑀ℎis the piecewise constant in time function that takes the value 𝑝𝑛 ℎon (𝑡𝑛−1, 𝑡𝑛). •𝒖𝑟∶[0, 𝑇] →𝑿𝑟is the piecewise linear in time function such that 𝒖𝑟(𝑡𝑛) =𝒖𝑛 𝑟. •𝑝𝑟∶[0, 𝑇] →𝑀𝑟is the piecewise constant in time function that takes the value 𝑝𝑛 𝑟on (𝑡𝑛−1, 𝑡𝑛). Theorem 5.2. Assume that the discrete inf-sup condition (3)holds, then the least-squares reduced pressure for problem (37)satisfies the following error estimates: If 𝑘 =0, ‖𝑝ℎ−𝑝𝑟‖𝐿1(𝐿2(Ω)) ≤𝐶ℎ(𝑑𝐿1(𝐿2(Ω))(𝑝ℎ,𝑀𝑟)+‖𝒖ℎ−𝒖𝑟‖𝐿2(𝐇1 𝐷(Ω))),(41) where 𝐶ℎis given by (28). Journal of Computational Physics 519 (2024) 113397 16 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Fig. 6. Pressure mean absolute error of the cavity flow 𝑅𝑒 = 10000. As for the first case first, we propose to complete this quantitative study with a qualitative evaluation. We consider the pressures enriched by 40 POD modes and calculated with the different versions presented above. In Figs. 6we present respectively, the pressure mean absolute error and its isovalues. We observe that the solution obtained using LSpROM(2) and k=0 offers an error level of the same order as that of the POD projection. 6.2.3. Singular pressure test case In this final example, we push to the limit the pressure recovery by proposing to test the robustness of the LSpROM approaches on a specific complex fabricated case. We aim to recover the reduced pressure solution that is both highly singular and requires an impractically fine mesh for its approximation. We therefore propose to consider the velocity solution Journal of Computational Physics 519 (2024) 113397 17 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Fig. 7. Illustration of the analytical truncated Weirstrass pressure scalar fields with a horizontal cut at 𝑦=1 2, in three times instances 𝑡=4,𝑡=8and 𝑡=12. 𝐮= 20 ∑ 𝑘=1 10−5 2(20−𝑘)cos(𝑘𝑡)[sin2(𝑘𝜋𝑥)sin(2𝑘𝜋𝑦) −sin(2𝑘𝜋𝑥)sin 2(𝑘𝜋𝑦)](50) while the pressure is an approximation of the Weirstrass function, which is continuous and nowhere differentiable 𝑝= 20 ∑ 𝑘=1 10−2 (4 3)(20−𝑘) sin(𝑘𝑡)cos(𝑘𝑡3∕2) 35 ∑ 𝑛=1 (1 2)𝑛 (cos((𝑘+2) 𝑛𝑥)+cos((𝑘+2) 𝑛𝑦))(51) To illustrate the irregularity of the pressure solution, three snapshots with a horizontal cut at 𝑦 =1 2corresponding to 𝑡 =4, 𝑡 =8and 𝑡 =12are plotted in Fig. 7. Since pressure acts as a Lagrange multiplier to ensure the divergence free constraint, the challenge of developing pressure recovery schemes is equally important in the context of unsteady linear Stokes equations. We therefore propose to present our results both for the unsteady Stokes and for the Navier-Stokes equations. Moreover it will also allow us to appreciate the effect of nonlinear terms on the level of approximation of the reduced pressure. Before we start the discussion on the numerical results, it is important to note that the irregularity of the selected pressure excludes the Laplacian-type method of Chacon et al. [5]that suggest the pressure to be in 𝐻1. Therefore, the next results will only be achievable using the proposed LSpROM technique capable of recovering the reduced pressure in 𝐿2. Finally, regard to the previous numerical experiments, we chose to focus and show only the results using LSpROM(2). Results and discussions Using the given analytical solutions, the source term is calculated and the problem is solved in a domain Ω =]0, 1[×]0, 1[, with a regular triangular mesh of size 64 ×64and for in time window [0, 12]. A choice of the viscosity is made to 𝜈=10 −2 and a first order semi-implicit Euler scheme with step Δ𝑡 =10 −3 is used for time integration. The POD bases are afterwards built with snasphots uniformally retrieved from the simulations with a time jump of 100Δ𝑡. Unlike in the two previous cases, the convergence of the pressure POD eigenvalues towards zero is less evident. Fig. 8shows that for both unsteady Stokes and Navier Stokes snapshots, there is slow decay of pressure POD eigenvalues towards zero with a stagnation at a level around 10−6 achieved after a significant number of truncation modes. This indicates that a large number of modes is necessary to accurately capture the irregular structures present in the given Weirstrass pressure. Pressure approximation errors presented in Fig. 9exhibit a conforming decay with respect to the POD number of modes up to the truncation order 20. This behavior is observed consistently across LSpROM approaches for 𝑘 =0and for 𝑘 =1. For higher mode numbers, a stagnation in error reduction is noted for both unsteady Stokes and Navier Stokes cases. We can also note that the level of error approximation for the Navier Stokes is a thousand times larger than that obtained for Stokes. However for both problems, the obtained errors predictions align with the result predicted by Theorem 5.1. Furthermore, one can note a difference in the error magnitude between the LSpROM(2) for 𝑘 =0and 𝑘 =1. Ended, choosing 𝑘 =0outperforms 𝑘 =1by nearly an order of magnitude. This order gap underscores once again the superior performance of the choice 𝑘 =0approach in maintaining lower error levels. Journal of Computational Physics 519 (2024) 113397 18 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Fig. 8. POD eigenvalues for the unsteady Stokes (left) and unsteady Navier-Stokes (right). Fig. 9. Pressure errors using constrained LSpROM approach 2: unsteady Stokes (left) and unsteady Navier-Stokes (right). 7. Conclusions In this paper we have introduced a method to recover the reduced pressure for incompressible flows. We have given some fundamental theoretical results concerning the existence and uniqueness of the solution whenever the full-order pair of velocity-pressure spaces is inf-sup stable, and we have proved an optimal error estimate for the reduced pressure. An additional highlight lies in the fact that our method is equivalent to solving the reduced mixed problem with reduced velocity basis enriched with the supremizers of the reduced pressure gradients. We have introduced a constrained treatment to avoid instabilities for high number of modes, that provide nearly-optimal errors, close to the POD projection ones, that provide error reductions of four order of magnitudes with respect to the pressure Poisson equation procedure. All our theoretical results have been confirmed by numerical experiments on stiff test cases. We have specifically demonstrated, using a manufactured solution, the difference in pressure calculation levels for the Stokes and Navier-Stokes problems. This highlights the impact of the non-linear terms on the recovery procedure. CRediT authorship contribution statement M. Azaïez: Writing – original draft, Methodology, Conceptualization. T. Chacón Rebollo: Writing – original draft, Methodology, Formal analysis, Conceptualization. M. Oulghelou: Writing – original draft, Validation, Software, Methodology, Conceptualization. I. Sánchez Muñoz: Writing – original draft, Methodology, Investigation, Formal analysis. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Azaiez Majdi reports was provided by Bordeaux Polytechnic Institute. Azaiez Mejdi reports a relationship with Bordeaux Polytechnic Institute that includes: employment. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability No data was used for the research described in the article. Journal of Computational Physics 519 (2024) 113397 19 M. Azaïez, T. Chacón Rebollo, M. Oulghelou et al. Acknowledgements The research of M. Azaiez has been supported by the French National Research Agency ANR-23-ASTR-0003 (SINRAM). The research of T. Chacón and I. Sánchez has been partially supported by the Spanish Research Agency and EU Feder Fund Grant RTI2018- 093521-B-C31, as well as the Marie Skłodowska-Curie Actions 872442 MSCA-RISE-2019 “ARIA” Project. References [1] A. Abdulle, O. Budać, A Petrov-Galerkin reduced basis approximation of the Stokes equation in parameterized geometries, C. R. Math. 353 (2015) 641–645, https://doi .org /10 .1016 /j .crma .2015 .03 .019. [2] I. Akhtar, A. Nayfeh, C. Ribbens, On the stability and extension of reduced-order Galerkin models in incompressible flows, Theor. Comput. Fluid Dyn. 23 (2009) 213–237, https://doi .org /10 .1007 /s00162 -009 -0112 -y. [3] F. Ballarin, A. Manzoni, A. Quarteroni, G. Rozza, Supremizer stabilization of POD-Galerkin approximation of parametrized steady incompressible Navier-Stokes equations, Int. J. Numer. Methods Eng. 102 (5) (2015) 1136–1161. [4] A. Caiazzo, T. Iliescu, V. John, S. Schyschlowa, A numerical investigation of velocity-pressure reduced order models for incompressible flows, J. Comput. Phys. 259 (2014) 598–616. [5] T. Chacón Rebollo, S. Rubino, M. Oulghelou, C. Alléry, Error analysis of a residual-based stabilization-motivated POD-ROM for incompressible flows, Comput. Methods Appl. Mech. Eng. 401 (2022) 115627. [6] Y. Choi, K. Carlberg, Space-time least-squares Petrov-Galerkin projection for nonlinear model reduction, SIAM J. Sci. Comput. 41 (2019) A26–A58, https:// doi .org /10 .1137 /17M1120531. [7] E. Fonn, H. van Brummelenb, T. Kvamsdalac, A. Rasheeda, Fast divergence-conforming reduced basis methods for steady Navier-Stokes flow, Comput. Methods Appl. Mech. Eng. 346 (2019) 486–512. [8] V. Girault, P.A. Raviart, Finite Element Approximations of the Navier-Stokes Equations, Springer, New York, 1986. [9] J.L. Guermond, P. Minev, J. Shen, An Overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Engrg. 195 (2006) 6011–6045. [10] Saddam Hijazi, Giovanni Stabile, Andrea Mola, Gianluigi Rozza, Data-driven POD-Galerkin reduced order model for turbulent flows, J. Comput. Phys. 416 (2020) 109513. [11] G. Rozza, D. Huynh, A. Manzoni, Reduced basis approximation and a posteriori error estimation for Stokes flows in parametrized geometries: roles of the inf-sup stability constants, Numer. Math. 125 (1) (2013) 115–152. [12] A. Ivagnes, G. Stabile, A. Mola, Pressure data-driven variational multiscale reduced order models, J. Comput. Phys. 476 (2023) 111904, https://doi .org /10 . 1016 /j .jcp .2022 .111904. [13] K. Kean, M. Schneier, Error analysis of supremizer pressure recovery for POD based reducedorder models of the time-dependent Navier-Stokes equations, SIAM J. Numer. Anal. 58 (4) (2020) 2235–2264. [14] H.P. Langtangen, A. Logg, Solving PDEs in Python: The FEniCS Tutorial I, 1st edition, Springer Publishing Company, Incorporated, 2017. [15] M. Mohebujjaman, L.G. Rebholz, T. Iliescu, Physically-constrained data-driven correction for reduced order modeling of fluid flows, Int. J. Numer. Methods Fluids 89 (3) (2019) 103–122. [16] B.R. Noack, P. Papas, P.A. Monkewitz, The need for a pressure-term representation in empirical Galerkin models of incompressible shear flows, J. Fluid Mech. 523 (2005) 339–365. [17] N.C. Nguyen, K. Veroy, A.T. Patera, Certified real-time solution of parametrized partial differential equations, in: R. Catlow, H. Shercliff, S. Yip (Eds.), Handbook of Materials Modeling, Kluwer Academic, 2005. [18] A. Quarteroni, G. Rozza, Numerical solution of parametrized Navier-Stokes equations by reduced basis methods, Numer. Methods Partial Differ. Equ. 23 (2007) 923–948. [19] A. Quarteroni, F. Saleri, A. Veneziani, Factorization methods for the numerical approximation of Navier-Stokes equations, Comput. Methods Appl. Mech. Eng. 188 (1) (2000) 505–526. [20] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer, Berlin, 1994. [21] D.V. Rovas, Reduced-Basis output bound methods for parametrized partial differential equations, Ph D Thesis, Massachusetts Institute of Technology, 2003. [22] G. Rozza, K. Veroy, On the stability of the reduced basis method for Stokes equations in parametrized domains, Comput. Methods Appl. Mech. Eng. 196 (7) (2007) 1244–1260. [23] G. Stabile, S. Hijazi, A. Mola, S. Lorenzi, G. Rozza, POD-Galerkin reduced order methods for CFD using finite volume discretisation: vortex shedding around a circular cylinder, Commun. Appl. Ind. Math. 8 (2017), https://content .sciendo .com /view /journals /caim /8 /1 /article -p210 .xml. [24] R. Temam, Navier Stokes Equations, North Holland, Amsterdam, 1984.