Tommaso Taddei
email: tommaso.taddei@inria.fr

Lei Zhang
email: lei.a.zhang@inria.fr

A discretize-then-map approach for the treatment of parameterized geometries in model order reduction

Keywords: parameterized partial differential equations, model order reduction, parameterized geometries. 2010 MSC: 65N30, 41A45, 35J57

We present a general approach for the treatment of parameterized geometries in projection-based model order reduction. During the offline stage, given (i) a family of parameterized domains {Ω µ : µ ∈ P} ⊂ R D where µ ∈ P ⊂ R P denotes a vector of parameters, (ii) a parameterized mapping Φ µ between a reference domain Ω and the parameter-dependent domain Ω µ , and (iii) a finite element triangulation of Ω, we resort to an empirical quadrature procedure to select a subset of the elements of the grid. During the online stage, we first use the mapping to "move" the nodes of the selected elements and then we use standard element-wise residual evaluation routines to evaluate the residual and possibly its Jacobian. We discuss how to devise an online-efficient reduced-order model and we discuss the differences with the more standard "map-then-discretize" approach (e.g., Rozza, Huynh, Patera, ACME, 2007); in particular, we show how the discretize-then-map framework greatly simplifies the implementation of the reduced-order model. We apply our approach to a two-dimensional potential flow problem past a parameterized airfoil, and to the two-dimensional RANS simulations of the flow past the Ahmed body.

Introduction

Treatment of parameterized geometries in model order reduction

In science and engineering, we are often interested in performing parametric studies to assess the influence of geometry on the solution to a given partial differential equation (PDE) of interest. We denote by {Ω µ : µ ∈ P} a family of parameterized domains in R D , D = 2, 3; then, we introduce a reference domain Ω and the parameterized mapping Φ : Ω × P → R D such that Φ µ := Φ(•; µ) is a bijection in Ω and Φ µ (Ω) = Ω µ . In this paper, we present a general approach for the treatment of parameterized geometries in projection-based parameterized model order reduction (pMOR). The approach exploits the knowledge of the mapping Φ: automatic construction of rapidly-computable mappings for parameterized geometries is beyond the scope of the present work; we refer to the model reduction literature (e.g., [START_REF] Lassila | Model order reduction in fluid dynamics: challenges and perspectives[END_REF][START_REF] Rozza | Basic ideas and tools for projection-based model reduction of parametric partial differential equations[END_REF]) for a thorough discussion.

To make the setting mathematically precise, we shall introduce a model problem considered in the numerical results. Given the family of domains {Ω µ : µ ∈ P}, we define the Dirichlet datum h : P → H 1/2 (∂Ω µ ); then, we seek u : P → H 1 (Ω µ ) such that u µ ∂Ωµ = h µ and

Ωµ ∇u µ • ∇v dx = 0 ∀ v ∈ H 1 0 (Ω µ ). (1)
Since the domain Ω µ depends on the parameter, pMOR techniques cannot be directly applied to devise a reduced-order model (ROM) for (1) -pMOR techniques are indeed designed to approximate fields over a parameter-independent domain.

To address this issue, several authors (e.g., [START_REF] Lassila | Model order reduction in fluid dynamics: challenges and perspectives[END_REF][START_REF] Rozza | Basic ideas and tools for projection-based model reduction of parametric partial differential equations[END_REF][START_REF] Ballarin | Fast simulations of patient-specific haemodynamics of coronary artery bypass grafts based on a pod-galerkin method and a vascular shape parametrization[END_REF][START_REF] Ballarin | Numerical modeling of hemodynamics scenarios of patient-specific coronary artery bypass grafts[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]) have followed a map-then-discretize (MtD) approach. First, given the parameterized mapping Φ, we introduce the weak formulation of the mapped problem for

u µ = u µ •Φ µ ∈ H 1 (Ω), Ω K µ ∇ũ µ • ∇v dx = 0 ∀ v ∈ H 1 0 (Ω), ũµ | ∂Ω = hµ , (2)
where K µ = det(∇Φ µ ) ∇Φ -1 µ ∇Φ -T µ and hµ = h µ • Φ µ . Second, we discretize (2) using a standard high-fidelity (HF) solver -e.g., the finite element (FE) method -and we apply traditional pMOR techniques to construct a ROM.

Although the MtD approach has been successfully applied to a broad class of PDEs, it suffers from two major limitations. First, since we effectively modify the variational formulation of the problem, the MtD approach might require the implementation of new HF routines: for the problem at hand, the HF solver for the mapped problem should deal with non-constant coefficients and also with terms of the form Ω ∂ xi u ∂ xj v dx with i = j that are not present in the original problem. Although this does not represent a major limitation for this specific test case, the need for implementing additional functionalities in the HF code might prevent the application of projection-based pMOR techniques to more challenging problems: in section 2.3, we show through a simple example why the MtD approach might be considerably more involved to apply for stabilized FE formulations. Second, it might be convenient to choose parameterized mappings that are piecewise smooth in Ω (see the numerical examples in [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF] and also the two model problems of this paper): in this case, the underlying mesh used to discretize (2) should be conforming with the coarse-grained partition associated with the definition of the mapping to ensure convergence of the HF solver at optimal polynomial rate. Discontinuities of K µ inside elements of the mesh pose indeed major issues for the HF quadrature rule and ultimately affect performance of the solver (cf. section 2).

Objective of the paper and outline

In this paper, we present a discretize-then-map (DtM) approach to deal with parameterized geometries for projection-based pMOR. During the offline stage, we resort to an empirical quadrature (EQ) procedure to optimally choose a low-dimensional subset of the grid elements; during the online stage, we first use the mapping to "move" the nodes of the selected elements and then we resort to standard element-wise residual evaluation routines to assemble the reduced system. We show that our approach does not require the HF mesh to be conforming with the coarse-grained partition; furthermore, it does not modify the equations and thus does not require modifications to the HF code.

Our point of departure is a high-order isoparametric continuous FE discretization: we shall consider both Galerkin and Petrov-Galerkin FE formulations. As explained in section 2, given a mesh over the domain Ω, FE fields (and thus reduced-order spaces) should be interpreted as pairs of a FE vector and a mapping Φ (cf. Equation ( 6)): the vector contains the "coefficients" of the solution with respect to the FE basis; the mapping defines the deformation of the nodes of the reference mesh. We define inner products and norms in the reference configuration, while we define parameter-dependent variational forms in the parameter-dependent (physical) configuration. To make our discussion mathematically precise, we explicitly work with the algebraic representation of the FE fields: depending on the choice of Φ, the latter can be linked to the solution to [START_REF] Lassila | Model order reduction in fluid dynamics: challenges and perspectives[END_REF] and to the solution to [START_REF] Rozza | Basic ideas and tools for projection-based model reduction of parametric partial differential equations[END_REF].

We discuss in detail implementation aspects that are key to integrate the ROM with available HF codes; in particular, we show that the DtM framework greatly simplifies the implementation of ROMs, particularly for nonlinear PDEs. For completeness, we also discuss how to enforce Dirichlet boundary conditions (BCs): more in detail, we discuss the use of lifting functions -also known as "control function method" -to deal with inhomogeneous Dirichlet BCs. We refer to the recent work [START_REF] Star | Extension and comparison of techniques to enforce boundary conditions in Finite Volume POD-Galerkin reduced order models for fluid dynamic problems[END_REF] and to the references therein for a thorough survey on the imposition of Dirichlet BCs in (Petrov-)Galerkin ROMs: we envision that our DtM approach can cope with both strong and weak enforcement of Dirichlet BCs.

The paper is organized as follows. In section 2, we introduce the main features of the DtM framework; furthermore, we compare the approach with the more standard MtD framework and we offer a number of remarks. In section 3, we discuss the application of model reduction in the DtM framework through the vehicle of a two-dimensional model problem of the form (1): we present the HF FE discretization, and we introduce the lifting method to deal with Dirichlet BCs. Then, in section 4, we apply our technique to the incompressible Reynolds-averaged Navier-Stokes (RANS) simulations of the flow past a twodimensional Ahmed body to demonstrate the applicability of the approach to nonlinear PDEs: as in the previous section, we present in detail the integration of the ROM with a HF stabilized FE model. Finally, in section 5, we draw some conclusions and we outline a number of subjects of ongoing research.

Relationship with previous works

The idea of deforming the reduced mesh as opposed to restate the equations in the reference domains has been previously considered in several works. Washabaugh et al. [START_REF] Washabaugh | On the use of discrete nonlinear reduced-order models for the prediction of steady-state flows past parametrically deformed complex geometries[END_REF] have resorted to a similar discretize-then-map approach for steady aerodynamics flows in the finite volume framework: the authors resort to minimum residual projection and rely on Gappy POD, [START_REF] Carlberg | The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF], for hyperreduction. Similarly, Dal Santo and Manzoni [START_REF] Santo | Hyper-reduced order models for parametrized unsteady navier-stokes equations on domains with variable shape[END_REF] have resorted to a Galerkin ROM in combination with the (matrix) discrete empirical interpolation method ( [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Negri | Efficient model reduction of parametrized systems by matrix discrete empirical interpolation[END_REF]) for hyper-reduction. However, in these works, the authors did not extensively compare the DtM approach with the MtD framework, and also did not discuss in detail the issue of dual residual norm estimation, which is key to assess the accuracy of the ROM during the online stage.

EQ procedures are a class of hyper-reduction techniques, which aim to devise an efficient quadrature rule to approximate a suitable set of parameterized integrals associated with the reduced problem. We here rely on the EQ procedure in [START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF] that combines the approaches in [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF] and [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF]. As explained in Remark 3 in the appendix, EQ procedures are related to the hyper-reduction procedure in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF]: a thorough introduction to EQ procedures for hyper-reduction of projection-based ROMs is beyond the scope of the present work. We refer to the recent reviews [START_REF] Farhat | Computational bottlenecks for PROMs: precomputation and hyperreduction[END_REF] and [START_REF] Yano | Model reduction in computational aerodynamics[END_REF] for a thorough discussion on EQ (also dubbed "mesh sampling and weighting'); we further refer to [START_REF] Chan | Entropy stable reduced order modeling of nonlinear conservation laws[END_REF] for a discussion of the stability properties of the hyper-reduced system for hyperbolic conservation laws.

In section 4, we resort to Least-Squares Petrov-Galerkin (LSPG, [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF]) projection to devise an efficient ROM for the RANS equations. As observed in [START_REF] Grimberg | On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows[END_REF][START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF], Petrov-Galerkin/minimum residual ROMs exhibit superior stability properties for non-coercive nonlinear PDEs. We remark that several authors have considered stabilized FE HF discretizations for ROM generation and construction, [START_REF] Giere | SUPG reduced order models for convection-dominated convection-diffusion-reaction equations[END_REF][START_REF] Pacciarini | Stabilized reduced basis method for parametrized advection-diffusion PDEs[END_REF][START_REF] Stabile | A reduced order variational multiscale approach for turbulent flows[END_REF] : the above-mentioned papers, however, do not consider geometric parameterizations.

Discretize-then-map treatment of parameterized geometries

Finite element spaces

For simplicity, we omit dependence on parameter. Given the domain Ω ⊂ R D , we define the FE mesh T = {D k } Ne k=1 , where D k ⊂ Ω denotes the k-th element of the mesh. We define the reference element D = {X ∈ [0, 1] D : D d=1 X d < 1} and the bijection Ψ k from D to D k for k = 1, . . . , N e ; we denote the Jacobian and Jacobian determinant of Ψ k by G k := ∇Ψ k and g k := det(G k ). Finally, we define the FE space of order p associated with the mesh T :

X T := v ∈ C(Ω) : v • Ψ k ∈ P p ( D), k = 1, . . . , N e . (3)
We consider a FE isoparametric discretization. We define the Lagrangian basis { i } n lp i=1 of the polynomial space P p ( D) associated with the nodes {X i } n lp i=1 ; then, we define the mappings

{Ψ k } k such that Ψ k (X) = n lp i=1 x hf i,k i (X), (4)
where {x hf i,k := Ψ k (X i ) : i = 1, . . . , n lp , k = 1, . . . , N e } are the nodes of the mesh, and we define the basis functions i,k

:= i • Ψ -1 k : D k → R. Note that Ψ k is completely characterized by the nodes in the k-th element X hf k := {x hf i,k } n lp
i=1 , k = 1, . . . , N e . We further introduce the nodes of the mesh {x hf j } N hf j=1 taken without repetitions and the connectivity matrix T ∈ N n lp ,Ne , such that x hf i,k = x hf T i,k . We observe that we might use a different polynomial order for the mappings {Ψ k } k : in view of mesh deformation (cf. ( 5)), it might be preferable for problems with straight boundaries to consider subparametric discretizations based on linear elemental maps (cf. section 2.2.2).

Given the mesh T over Ω and the bijection Φ : Ω → Φ(Ω), we introduce the mapped mesh Φ(T ) that shares with T the same connectivity matrix T and has nodes {Φ(x hf j )} N hf j=1 . We denote by Ψ k,Φ the elemental mapping associated with the k-th element D k,Φ of Φ(T ); Ψ k,Φ is given by

Ψ k,Φ (X) = n lp i=1 Φ(x hf i,k ) i (X). (5)
We also introduce the FE space X Φ(T ) .

Given the vector u ∈ R N hf , we define the corresponding FE field u Φ ∈ X Φ(T ) such that u Φ ∈ C ( k D k,Φ ) and

u Φ D k,Φ = n lp i=1 (u) T i,k i,k,Φ , k = 1, . . . , N e , (6a)
where i,k,Φ = i • Ψ -1 k,Φ . If Φ is the identity map, id(x) ≡ x, we use notation u = u Φ=id ∈ X T . (6b)
Note that by changing Φ we effectively change the function field associated with the vector u: we should thus interpret FE fields (and thus reduced-order spaces) as pairs of a FE vector and a mapping.

MtD and DtM formulation of the Laplace problem

In the DtM approach, we directly discretize (1): given the mesh T , the parameterized mapping Φ and the parameter µ ∈ P, we seek

u hf µ,Φµ ∈ X µ := X Φµ(T ) such that Ωµ ∇u hf µ,Φµ • ∇v dx = 0 ∀ v ∈ X µ ∩ H 1 0 (Ω µ ), u hf µ,Φµ | ∂Ωµ = h µ , (7)
with Ω µ = Φ µ (Ω). On the other hand, in the MtD approach, we discretize the mapped problem (2): we seek u hf µ ∈ X := X T such that

Ω K µ ∇u hf µ • ∇v dx = 0 ∀ v ∈ X ∩ H 1 0 (Ω), u hf µ | ∂Ω = hµ . (8)
In the remainder of this section, we compare the two formulations for the Laplace problem.

Equivalence for piecewise-polynomial maps

Exploiting the change-of-variable formula, we have

Ne k=1 D k K µ ∇w • ∇v dx = Ne k=1 Φµ(D k ) ∇w • ∇v dx. Note that Ψ k,Φµ (5) coincides with Φ • Ψ k in the nodes {X i } n lp i=1 . However, Ψ k,Φ and Φ • Ψ k coincide in D -and thus D k,Φ = Φ(D k ) -only if Φ • Ψ k ∈ P p for k = 1, . . . , N e .
In conclusion, we obtain the result contained in Lemma 1: the proof is a tedious but straightforward application of the chain rule and is here omitted.

Lemma 1. Consider an isoparametric FE discretization of order p. Then, [START_REF] Washabaugh | On the use of discrete nonlinear reduced-order models for the prediction of steady-state flows past parametrically deformed complex geometries[END_REF] and [START_REF] Carlberg | The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF]

are equivalent if Φ • Ψ k ∈ P p for k = 1, . . . , N e .

Optimal convergence and discrete bijectivity

In order to highlight relevant features of the MtD and DtM approaches, we consider the problem

-∂ xx u = sin(πx) x ∈ Ω = (0, 1), u(0) = u(1) = 0, (9a)
which admits the unique solution u(x) = 1 π 2 sin(πx). We further define the mapping

Φ(x) =      1 2 x x ≤ x 0 , 1 2 x 0 + 2 -x 0 1 -x 0 (x -x 0 ) x ≥ x 0 ; x 0 = 1 √ 2 . ( 9b
)
Clearly, Φ is a Lipschitz map from Ω in itself. We denote by T an uniform FE grid of Ω with N e elements. Problem ( 9) is intended to study the performance of the two strategies when combined with piecewise smooth maps, which are broadly used in MOR.

In Figure 1(a), we show the behavior of the L 2 error associated with three FE discretizations: a P3 isoparametric FE discretization of the problem associated with Φ(T ) (P3 DtM); a P1-P3 subparametric FE discretization of the problem associated with Φ(T ) (P1-P3 DtM ); a P3 FE discretization of the mapped problem associated with T (P3 MtD). We observe that the MtD approach fails to recover the optimal convergence rate (r = 4) due to quadrature error and due to the fact that the MtD approach approximates the mapped field u • Φ, which is considerably less smooth than u. The isoparametric discretization suffers due to stability issues: even if Φ is bijective, the local map Ψ k,Φ that contains x 0 is not necessarily bijective. On the other hand, the subparametric DtM approach recovers the optimal convergence rate as expected from standard FE theory (e.g., [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]Proposition 3.4.1]). Note that for P1 discretizations in one-dimensional domains Φ(T ) is guaranteed to not have inverted elements; however, the same result does not hold in higher dimensions. Consider the element depicted in Figure 1(b), the dashed line denotes the boundary of the domain: we can construct a mapping Φ that maps the point A in the point B in the Figure . The deformed element thus reduces to a straight line regardless of the polynomial degree -which implies that the corresponding elemental mapping Ψ k,Φ is not bijective. Note that for the MtD approach well-posedness depends on the fact that K µ in ( 8) is positive definite: the latter is independent of the underlying mesh and thus follows from the bijectivity of Φ. We conclude that, unlike MtD, the DtM approach might fail for smooth but large deformations, particularly in the presence of anisotropic elements.

Application of hyper-reduction methods

Hyper-reduction refers to a broad class of methods that aims to reduce the online assembling cost for projection-based ROMs; hyper-reduction techniques might be applied either at the continuous or at the discrete level -continuousbased and discrete-based hyper-reduction. Continuous-based hyper-reduction methods (e.g., [START_REF] Ballarin | Fast simulations of patient-specific haemodynamics of coronary artery bypass grafts based on a pod-galerkin method and a vascular shape parametrization[END_REF][START_REF] Ballarin | Numerical modeling of hemodynamics scenarios of patient-specific coronary artery bypass grafts[END_REF][START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF]) aim to devise an affine approximation of the integrand to ultimately obtain a parametrically-affine (cf. [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]) approximation of the problem:

R aff (u hf µ , v) = Qa q=1 Θ q µ R q (u hf µ , v) = 0,
where R 1 , . . . , R Qa are parameter-independent forms that are polynomial in the first argument and linear in the second argument, and Θ 1 µ , . . . , Θ Qa µ are parameter-dependent coefficients that can be rapidly evaluated. On the other hand, discrete-based hyper-reduction methods (e.g., [START_REF] Carlberg | The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF][START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF]) aim to identify a "reduced integration domain".

The MtD approach can cope with both continuous and discrete methods; furthermore, if the mapping Φ is piecewise-linear and discontinuities of Φ align with the elements' facets of the mesh, ( 8) admits a parametrically-affine decomposition without having to apply hyper-reduction. On the other hand, the DtM approach can only cope with discrete-based hyper-reduction techniques, which require integration over a portion of the domain.

It is important to observe that parametrically-affine problems have been extensively studied in the model reduction literature and there exist specialized techniques to tackle them. In particular, rigorous a posteriori error estimation techniques are largely restricted to parametrically-affine problems and thus fit exclusively with the MtD framework.

• Estimation of residual dual norm is straightforward for affine problems (cf. [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]); furthermore, for certain classes of linear and nonlinear PDEs, the successive constraint method (SCM, [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF][START_REF] Huynh | A naturalnorm successive constraint method for inf-sup lower bounds[END_REF][START_REF] Yano | A space-time Petrov-Galerkin certified reduced basis method: application to the Boussinesq equations[END_REF]) can be applied to obtain a rigorous error bound for the prediction error. In addition, the strategy in [START_REF] Schmidt | Rigorous and effective a-posteriori error bounds for nonlinear problems-application to rb methods[END_REF] can later be applied to significantly sharpen the error bound.

• For a class of parametrically-affine PDEs, exact estimates can be derived to bound the error with respect to the exact solution ( [START_REF] Ali | Reduced basis methods with adaptive snapshot computations[END_REF][START_REF] Yano | A minimum-residual mixed reduced basis method: Exact residual certification and simultaneous finite-element reduced-basis refinement[END_REF]): this class of bounds is extremely important to devise spatio-parameter adaptive strategies for parametric problems, [START_REF] Yano | A reduced basis method for coercive equations with an exact solution certificate and spatio-parameter adaptivity: energy-norm and output error bounds[END_REF].

In this work, we adapt the approach in [START_REF] Taddei | An offline/online procedure for dual norm calculations of parameterized functionals: empirical quadrature and empirical test spaces[END_REF] to devise a rapid and reliable estimator of the residual dual norm; furthermore, we show that the dual residual norm is highly correlated with the actual error: we can thus apply the approach in [START_REF] Drohmann | The ROMES method for statistical modeling of reduced-order-model error[END_REF] to devise effective probabilistic error bounds.

Implementation for stabilized FE discretizations

As stated in the introduction, practical implementation of the MtD approach might be considerably more involved, particularly for stabilized FE formulations and nonlinear problems. To motivate our claim, consider the streamlineupwind/ Petrov-Galerkin, (SUPG, [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF]) discretization of the advection-diffusion equation -∆u

+ b • ∇u = f in Ω µ , u| ∂Ωµ = 0;
which can be stated as follows:

Ne k=1 D k,Φµ (∇u • ∇v + b • ∇uv -f v) dx + αh k,Φµ D k,Φµ R(u) b b 2 • ∇v dx = 0, where R(u) := -∆u + b • ∇u -f is the strong residual, α > 0 is a user-defined constant, h k,Φµ = |D k,Φµ | 1/D and |D k,Φµ | is the volume of the element.
In order to apply the MtD approach, we should rewrite -∆u in the reference configuration: if we denote by ∇ the gradient with respect to the reference coordinates, and we define ũ := u • Φ µ , ṽ := u • Φ µ , we find

D k,Φµ ∆u b b 2 • ∇v dx = D k ∇Φ -T µ ∇ • ∇Φ -T µ ∇ ũ ∇Φ -1 µ b b 2
• ∇ṽ dx Note that the expression at right-hand side requires the implementation of new local integration routines to deal with the mapped terms; in addition, it requires the evaluation of the mapping Φ µ , the gradient ∇Φ µ and the Hessian HΦ µ in the quadrature points. We envision that computation of second-order derivatives might be involved and computationally expensive; furthermore, second-order derivatives are not well-defined for Lipschitz maps. On the other hand, the DtM approach only requires evaluation of Φ µ in the nodes of the mesh (cf. ( 5)); furthermore, it allows to reuse local integration routines that are used in the HF solver, cf. Algorithms 3 and 4 in section 4.

Theory and results for the linear case

We here present a pMOR strategy to estimate elements of the manifold M := {u hf µ : µ ∈ P} ⊂ U, where u hf µ denotes the vector associated with the FE approximation to (1) for a given µ ∈ P. Given an estimate u µ of u hf µ for some µ ∈ P, recalling the definitions in (6), we might define the FE fields u µ ∈ X T and u µ,Φµ ∈ X Φµ(T ) : provided that the mapping is well-defined, we expect that

u hf µ -u µ,Φµ H 1 (Ωµ) ≈ ũµ -u µ H 1 (Ω) .
In section 3.1, we introduce useful notation and definitions. In section 3.2, we present the HF formulation. Then, in section 3.3, we derive a Galerkin ROM based on lifting: we discuss the offline/online computational decomposition, empirical interpolation of the Dirichlet BCs, the introduction of an empirical quadrature rule for the reduction of online costs to assemble the reduced matrix and vector, and a posteriori error estimation. In section 3.4, we present the model problem considered in the numerical experiments, a two-dimensional potential flow past a parameterized airfoil. Finally, in section 3.5, we present the numerical results to prove the effectiveness of the approach.

Preliminary definitions

We omit dependence on parameter. Given the Gaussian quadrature rule in D, {(ω g q , X g q )} nq q=1 , we define the high-fidelity quadrature (hfq) rule in Ω as

ω hfq q,k := ω g q g k (X g q ), x hfq q,k = Ψ k (X g q ), q = 1, . . . , n q , k = 1, . . . , N e . ( 10
)
Furthermore, we define L ∈ R nq,n lp and L ∇ ∈ R nq,n lp ,Ne,D such that

(L) q,i = i (X g q ), (L ∇ ) q,i,k,d = ∂ x d i,k (x hfq q,k ), (11)
with q = 1, . . . , n q , i = 1, . . . , n lp , k = 1, . . . , N e , d = 1, . . . , D. Note that i (X g q ) = i,k (x hfq q,k ) for all k = 1, . . . , N e . Quadrature rule and shape functions are used for FE assembling; the HF quadrature rule also provides the baseline for the subsequent hyper-reduction procedure. Given the mapping Φ, we can introduce the counterparts L Φ , L ∇,Φ of L, L ∇ defined in [START_REF] Negri | Efficient model reduction of parametrized systems by matrix discrete empirical interpolation[END_REF]: explicit expressions are obtained by replacing 11): since L Φ = L for any choice of Φ, we shall omit the subscript Φ below.

i,k = i • Ψ -1 k with i,k,Φ = i • Ψ -1 k,Φ in (
We define the Hilbert space

U := (R N hf , (•, •)) where the inner product (•, •) is given by (w, v) = Ω ∇w • ∇v + w • v dx. ( 12
)
We further define the induced norm

• := (•, •). If we define the H 1 norm in U ⊂ R D as w H 1 (U ) := w 2 L 2 (U ) + ∇w 2 L 2 (U ) , we obtain that u = u H 1 (Ω) ; furthermore, if Φ is a diffeomorphism, there exist c, C > 0 such that c u ≤ u Φ H 1 ( k D k,Φ ) ≤ C u for all u ∈ U.
We define the "elemental" restrictions such that

E k u = [(u) T 1,k , . . . , (u) T n lp ,k ]
T , for k = 1, . . . , N e : the restriction E k depends exclusively on the connectivity matrix T and is thus independent of the mapping Φ. Given u ∈ U, we define u un ∈ R n lp ,Ne such that (u un ) i,k = (E k u) i . Similarly, given the reduced-order basis Z = [ζ 1 , . . . , ζ N ] ∈ R N hf ,N , we define the tensor Z un ∈ R n lp ,Ne,N such that

(Z un ) i,k,n = (ζ un n ) i,k = (E k ζ n ) i ,    i = 1, . . . , n lp , k = 1, . . . , N e , n = 1, . . . , N. (13)
We denote by Γ dir ⊂ ∂Ω the Dirichlet boundary; in this work, we assume that each facet of the mesh is either strictly contained in Γ dir or has empty intersection with Γ dir . We introduce the set of indices of the mesh I dir ⊂ {1, . . . , N hf } belonging to Γ dir ; we denote by M hf := |I dir | the cardinality of the set. Then, we define the space U 0 := {v ∈ U : v(I dir ) = 0} endowed with the norm • , and the space

H = (R M hf , • dir = (•, •) dir ) where the norm
• dir is introduced below. Furthermore, we introduce the extension operator

H : H → U such that for all h ∈ H Hh(I dir ) = h, (Hh, v) = 0 ∀ v ∈ U 0 . (14)
H is used in section 3.3 to derive the Galerkin ROM for the lifted solution field. Finally, we define the norm • dir such that for all w ∈ U

w(I dir ) dir = Ne k=1 ∂D k ∩Γ dir (w(x)) 2 dx. (15)

High-fidelity formulation of the model problem

As in the previous section, we omit dependence on parameter. We discretize (1) using a CG FE method: given the mesh T , the mapping Φ, and the indices I dir ⊂ {1, . . . , N hf } associated with Γ dir = ∂Ω, we denote by h ∈ R M hf the vector of boundary conditions such that (h

) i = h Φ(x dir i ) with x dir i = x hf (I dir )i , i = 1, . . . , M hf := |I dir |. Then, we seek u hf ∈ U such that u hf (I dir ) = h and R hf (u hf , v, Φ) = Ne k=1 r E k u hf , E k v, Φ X hf k = 0 ∀ v ∈ U 0 , (16a)
where the local residual r is given by

r w, v, Φ X hf k = nq q=1 ω hfq q,k,Φ • (∇w Φ • ∇v Φ ) x=x hfq q,k,Φ = nq q=1 ω hfq q,k,Φ • n lp j=1 (v) j • n lp i=1 (w) i • D d=1 (L ∇,Φ ) q,i,k,d • (L ∇,Φ ) q,j,k,d = v T A un k,Φ w , ∀ w, v ∈ R n lp .
(16b) Here, w Φ , v Φ are defined as in [START_REF] Star | Extension and comparison of techniques to enforce boundary conditions in Finite Volume POD-Galerkin reduced order models for fluid dynamic problems[END_REF], and A un k,Φ is given by

A un k,Φ i,j = nq q=1 ω hfq q,k,Φ • D d=1
(L ∇,Φ ) q,i,k,d • (L ∇,Φ ) q,j,k,d , i, j = 1, . . . , n lp .

(16c) Given the solution u hf to (16), we might define the corresponding FE fields u hf and u hf Φ using (6): the former is an approximation of the solution ũ to (2); the latter is an approximation of the solution to [START_REF] Lassila | Model order reduction in fluid dynamics: challenges and perspectives[END_REF].

Computation of (L ∇,Φ ) •,•,k,• ∈ R nq,n lp ,D for a given k is dominated by the cost of computing ∇Ψ k in all quadrature points: the latter involves the multiplication of a D × n lp matrix (associated with the nodes) with a n lp × n q D matrix (associated with the gradient of the shape functions). Assuming n q = O(n lp ), computation of (L ∇,Φ ) •,•,k,• requires O(n2 lp ) flops1 ; similarly, computation of quadrature weights {ω hfq q,k,Φ } q and assembling of A un k,Φ using (16c) requires O(n3 lp ) flops. We remark that, thanks to a large body of research in the past few decades in the FE community, efficient fully-vectorized routines are now available to compute local tensors of the form (16c) for structured and unstructured grids, for a broad range of PDEs; furthermore, the loop over the elements to compute {A un k,Φ } k can be trivially parallelized.

Galerkin reduced-order model for the lifted solution field

Given the reduced-order basis Z = [ζ 1 , . . . , ζ N ] ∈ R N hf ,N such that Z(I dir , : ) = 0, we aim to approximate the solution u hf µ for a given µ ∈ P as

u µ = Zα µ + Hh µ , where R hf µ (u µ , ζ n ; Φ µ ) = 0, n = 1, . . . , N. (17)
Note that computation of u µ is highly inefficient: first, evaluation of Hh µ requires the solution to a problem of size N hf ; second, the assembling of the N ×N system associated with [START_REF] Yano | Model reduction in computational aerodynamics[END_REF] requires integration over the whole FE mesh and thus scales with the total number of elements N e . We shall resort to the empirical interpolation method (EIM, [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF]) to obtain a low-dimensional approximation of the Dirichlet datum of the form

h µ = H h µ (I ei ), where H ∈ R M hf ,M , I ei ⊂ {1, . . . , M hf }, |I ei | = M M hf . ( 18
)
Since H is parameter-independent, we can precompute the surrogate basis W = [ψ 1 , . . . , ψ M ] = HH ∈ R N hf ,M ; finally, we define the approximate lift e µ = Wh µ (I ei ).

Exploiting (16c), it is easy to verify that

R hf µ Z α hf µ + e µ , ζ n ; Φ µ = 0 ⇔ A hf µ α hf µ = F hf µ (19a)
where

           A hf µ := Ne k=1 A k,µ , A k,µ n,n = (E k ζ n ) T A un k,Φµ E k ζ n F hf µ := Ne k=1 F k,µ , F k,µ n = -(E k ζ n ) T A un k,Φµ E k e µ . (19b)
Assuming that {E k Z} k and {E k W} k are precomputed during the offline stage, the cost of assembling the reduced system scales with O(N e n lp ((N + n lp ) 2 + M )). To reduce assembling costs, we replace the HF residual with the weighted residual R eq µ ,

R eq µ (w, v, Φ µ ) = Ne k=1 ρ eq k r E k w, E k v, Φ µ X hf k , (20)
where ρ eq = [ρ eq 1 , . . . , ρ eq Neq ] T is a sparse vector of positive weights such that ρ eq 0 = Q N e . We denote by I eq = {k ∈ {1, . . . , N e } : ρ eq k > 0} the set of indices associated with the sampled elements.

In conclusion, given µ ∈ P, we define the estimate u µ = Z α µ + e µ of the solution u hf µ to [START_REF] Farhat | Computational bottlenecks for PROMs: precomputation and hyperreduction[END_REF], such that

e µ = W h µ (I ei ), R eq µ (Z α µ + e µ , ζ n , Φ µ ) = 0, n = 1, . . . , N. (21a)
Exploiting previous identities, we find that α µ satisfies

A eq µ α µ = F eq µ , where A eq µ = k∈Ieq ρ eq k A k,µ , F eq µ = k∈Ieq ρ eq k F k,µ .
(21b) Online assembling of the reduced system requires the computation of the Q • n lp nodes of the mapped mesh, {Φ µ (x hf i,k ) : i = 1, . . . , n lp , k ∈ I eq }, the assembling of the elemental matrices {A un k,Φµ } k∈Ieq in (16c), the computation of { A k,µ } k in (19b), and then the summation over the sampled elements. Provided that evaluation of Φ µ in a given point requires O(1) flops, the overall assembling cost scales with O(Qn lp ((N + n lp ) 2 + M )).

Algorithm 1 Offline-online computational decomposition

Offline stage 1: Compute samples of the HF solution manifold {u hf µ k } ntrain k=1 and compute the lifted fields {ů hf

µ k = u hf µ k -Hh µ k } k 2: Compute the POD space Z = POD {ů hf µ k } k , (•, •), tol pod .
3: Compute H and I ei in (18) using EIM and then set W = HH (cf. section 3.3.1).

4: Find the sparse vector of weights ρ eq in ( 20) (cf. section 3.3.2).

5: Construct the offline structures for dual residual norm calculation (cf. section 3.3.3).

Online stage (for any given µ ∈ P)

1: Compute the deformed nodes {Φ µ (x hf i,k ) : i = 1, . . . , n lp , k ∈ I eq }.

2: Assemble A eq µ and F eq µ in (21b).

3: Solve A eq µ α µ = F eq µ and return u µ = Z α µ + Wh µ (I ei ). 4: Estimate the dual residual norm (26) (cf. section 3.3.3).

In Algorithm 1, we outline the offline/online computational decomposition considered in the numerical experiments; the procedure includes the offline construction and the online evaluation of the dual residual estimator that is used to assess the accuracy of the ROM. Computation of ( 21) requires the storage of the reduced-order bases Z, W and of the nodes {x hf i,k : i = 1, . . . , n lp , k ∈ I eq } in the selected elements (n lp Q(M + N + D) doubles), and the storage of quantities (shape functions, quadrature rule) associated with the master element D (O(n lp n q ) doubles ) -computational and storage costs do not consider dual residual norm estimation, which is addressed in section 3.3.3.

We here resort to a non-adaptive (deterministic or random) sampling of the parameter domain and we resort to proper orthogonal decomposition (POD, [START_REF] Volkwein | Model reduction using proper orthogonal decomposition[END_REF]) to generate the reduced-order basis for the ROM. The size of the spaces is chosen according to the criterion

N := min    N : N n=1 λ n ≥ (1 -tol pod ) ntrain i=1 λ i    , (22)
where {λ i } i are the eigenvalues in descending order of the POD covariance matrix, and tol pod is a suitable tolerance. Given snapshots

{w k } k ⊂ R N pod , notation W = POD {w k } k , (•, •) pod , tol pod signifies that W ∈ R N pod
,N is constructed using POD based on the inner product (•, •) pod and tolerance tol pod in [START_REF] Pacciarini | Stabilized reduced basis method for parametrized advection-diffusion PDEs[END_REF].

In the next three sections, we review how to construct the approximate lift e, the weights ρ eq ∈ R Ne in [START_REF] Grimberg | On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows[END_REF], and how to estimate the dual residual error indicator; we state upfront that a thorough discussion and analysis of EIM and EQ procedures is beyond the scope of the present work: we refer to [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF] (see also [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF]) and to the reviews [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF] for a detailed presentation of EIM; on the other hand, we refer to the previously mentioned papers [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Farhat | Computational bottlenecks for PROMs: precomputation and hyperreduction[END_REF][START_REF] Yano | Model reduction in computational aerodynamics[END_REF] for a discussion on EQ procedures for linear and nonlinear problems.

Affine approximation of the Dirichlet datum

Given the parameters {µ k } ntrain k=1 ⊂ P, we compute {h k = h µ k } k and we apply POD based on the (•, •) dir inner product to obtain a M -dimensional approximation space H M = span{ξ m } M m=1 ⊂ H dir where the integer M is chosen according to [START_REF] Pacciarini | Stabilized reduced basis method for parametrized advection-diffusion PDEs[END_REF]. Then, we pass Ξ := [ξ 1 , . . . , ξ M ] to the EIM Greedy algorithm (see, e.g., [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF]Algorithm 10.3]) to find the set of indices I ei ⊂ {1, . . . , M hf }, |I ei | = M . Finally, we define the matrix W in (21a) as W = HΞ (Ξ(I ei , :)) -1 .

We might avoid the use of EIM to approximate the Dirichlet datum by resorting to a weak imposition of the Dirichlet condition -see, e.g., [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF] and also [START_REF] Babuška | The finite element method with penalty[END_REF][START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF] for a discussion in the continuous setting. In this case, we can use EQ to simultaneously deal with hyper-reduction and BC approximation. Since we do not pursue this strategy in this paper, we do not further comment on this issue.

Empirical quadrature procedure

Following [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF], we seek ρ eq ∈ R Ne + such that (i) the number of nonzero entries in ρ eq , ρ eq 0 , is as small as possible; (ii, constant function constraint) the constant function is approximated correctly in Ω (i.e., Φ = id),

Ne k=1 ρ eq k |D k | -|Ω| 1; (23)
(iii, manifold accuracy constraint) for all µ ∈ P train = {µ k } ntrain k=1 , the empirical residual satisfies

A hf µ -1 A eq µ α hf µ -F eq µ 2 1. ( 24
)
We refer to [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF] for a detailed justification of ( 23) and ( 24) for linear and nonlinear PDEs.

It is easy to verify that ( 23) and ( 24) could be rewritten in matrix form as

Ne k=1 ρ eq k |D k | -|Ω| = |G const ρ eq -|Ω| 1, with G const = [|D 1 |, . . . , |D Ne |],
and

A hf µ -1 A eq µ α hf µ -F eq µ 2 = G µ ρ eq 2 1,
with

G µ = A hf µ -1 A 1,µ α hf µ -F 1,µ , . . . , A Ne,µ α hf µ -F Ne,µ .
The problem of finding ρ eq can thus be reformulated as a sparse-representation (or bestsubset selection) problem:

min ρ∈R Ne ρ 0 , s.t Gρ -b ≤ δ; ρ ≥ 0; (25a)
for suitable choices of the vector norm • and the tolerance δ > 0, and

G =       G const G µ 1 . . . G µ n train       , b = |Ω| 0 (25b)
A large body of research focuses on the development of robust algorithms to find approximate solutions to [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF]. Following [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF], we resort to the nonnegative linear least-squares method (see [START_REF] Lawson | Solving least squares problems[END_REF]) -more in detail, we rely on the Matlab function lsqnonneg, which takes as input the pair (G, b) and a tolerance tol eq > 0 and returns the sparse vector ρ eq , [ρ eq ] = lsqnonneg (G, b, tol eq ) .

We remark that for large-scale problems computation of the solution to the nonnegative linear least-squares problem might be prohibitively expensive: to address this issue, efficient partitioned approaches have been developed in [START_REF] Grimberg | Mesh sampling and weighting for the hyperreduction of nonlinear petrov-galerkin reduced-order models with local reduced-order bases[END_REF].

A posteriori error estimation

Due to the absence of reliable a priori error estimates, it is of paramount importance to assess the accuracy of the ROM during the online stage; for this reason, a large body of work has been devoted to the development of rapid and reliable a posteriori error indicators -see [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF] and the references therein. Most proposals rely on the computation (estimation) of the dual residual norm

R hf µ ( u µ ) := sup v∈U0 R hf µ ( u µ , v; Φ µ ) v . ( 26
)
for all µ ∈ P train,r . Proceeding as in section 3.3.2, we find that R eq,r µ -R hf,r

µ 2 = G r µ ρ eq,r -R hf,r µ 2 with G r µ = [ R r 1,µ , . . . , R r Ne,µ ]: min ρ∈R Ne ρ 0 , s.t G r ρ -b r ≤ δ; ρ ≥ 0; G r =     G r µ 1 r . . . G r µ n train,r r     b r =     R hf,r µ 1 r . . . R hf,r µ n train,r r     .
(32) Similarly to [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF], approximate solutions to (32) can be obtained using the nonnegative least-squares method.

Algorithm 2 summarizes the offline procedure to compute the empirical test space Y r Jr and the weights ρ eq,r , and the online procedure to rapidly compute R eq µ . Note that the offline procedure takes as input the ROM and the tolerances tol es (for Y r Jr ) and tol eq,r (for ρ eq,r ). In Appendix B (cf. Remark 6), we present a rigorous justification of our approach.

Algorithm 2 Offline-online estimation of the dual residual norm Offline stage

1: Given P train,r = {µ k r } ntrain,r k=1
, estimate u µ by solving [START_REF] Giere | SUPG reduced order models for convection-dominated convection-diffusion-reaction equations[END_REF] for all µ ∈ P train,r .

2: Compute the Riesz representer ψ r µ using [START_REF] Ali | Reduced basis methods with adaptive snapshot computations[END_REF] for all µ ∈ P train,r . 5: [ρ eq,r ] = lsqnonneg (G r , b r , tol eq,r ) .

Online stage (for any given µ ∈ P and estimate α µ ) 1: Compute the local residuals { R r k,µ : k ∈ I eq,r } using (29b).

2: Compute R eq,r µ using (29a) and return R eq µ = R eq,r µ 2 .

Model problem: potential flow past a parameterized airfoil

We introduce the domain

Ω µ = Ω box \ Ω naca,µ ⊂ R 2 such that Ω box = (x min , x max ) × (-H, H), x min = -2, x max = 6, H = 4, and Ω naca,µ = x ∈ (0, 1) 2 : -f naca (x 1 , µ 1 ) < x 2 < f naca (x 1 , µ 2 ) , (33)
where f naca (s, th) = 5th 0.2969 √ s -0.1260s -0.3516s 2 + 0.2843s 3 -0.1036s 4 , and µ 1 , µ 2 are positive parameters. Then, we introduce problem (1) with Dirich-let datum h given by

h µ (x) =                            0 on Γ btm = (x min , x max ) × {-H} 1 on Γ top = (x min , x max ) × {H}, x 2 + H 2H on Γ out = {x max } × (-H, H), hµ ( x2+H 2H ) -hµ (0) hµ (1) -hµ (0) on Γ in = {x min } × (-H, H), 1 on ∂Ω naca,µ , (34a)
with

hµ (t) = 1 2 1 + 1 π arctan (10(t -µ 3 )) + 1 π arctan (10(t -µ 4 )) . (34b)
Finally, we shall define the vector of parameters µ = [µ 1 , µ 2 , µ 3 , µ 4 ]: in the following we assume that µ belongs to the compact region

P = [0.09, 0.15] 2 × [0.1, 0.3] × [0.6, 0.8]. (34c)
In order to define the mapping Φ, (i) we partition Ω µ into the four regions {Ω i,µ } 4 i=1 depicted in Figure 2(a) and we define Ω := Ω μ and {Ω i = Ω i,μ } 4 i=1 with μ = [0.12, 0.12, 0, 0]; (ii) we define the parameterized Gordon-Hall maps (see [START_REF] Gordon | Construction of curvilinear co-ordinate systems and applications to mesh generation[END_REF]) Ψ gh,i : (0, 1) 2 × P → Ω i,µ for i = 1, . . . , 4; and (iii) we define Φ such that Φ µ (x) = Ψ gh,i Ψ -1 gh,i (x, μ) , µ , if x ∈ Ω i , i = 1, 2, 3, 4; [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF] note that Φ µ coincides with the identity map id for µ = μ and Φ µ ≡ id in Ω 1 ∪Ω 4 for all µ ∈ P. We remark that several other approaches might be considered to construct exact or approximate mappings for parameterized geometries: we refer to the pMOR literature for a thorough overview (see in particular [START_REF] Lassila | Model order reduction in fluid dynamics: challenges and perspectives[END_REF][START_REF] Rozza | Basic ideas and tools for projection-based model reduction of parametric partial differential equations[END_REF]). Figure 2(b) shows the mesh used in our implementation. Note that the mesh is not conforming with the coarse-grained partition. To efficiently evaluate [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF], during the offline stage, we compute I Φ ∈ {1, 2, 3, 4} N hf and {x hf,ref j } N hf j=1 ⊂ [0, 1] 2 such that (I Φ ) j denotes the label of the region to which the j-th node of the mesh x hf j belongs, and x hf,ref j := Ψ -1 gh,(IΦ)j x hf j , μ . Then, given a new value of the parameter µ ∈ P, we simply compute the mapped nodes using the identity

Φ µ (x hf j ) = Ψ gh,(IΦ)j x hf,ref j , µ , j = 1, . . . , N hf . (36)

Numerical results

We consider a p = 3 FE discretization with N e = 2816 elements (N hf = 12840). To train our model, we consider n train = 200 snapshots associated with randomly-sampled parameters {µ k } ntrain k=1 . We compute the EIM approximation using a tolerance tol eim = 10 -14 : for this choice of the tolerance, we find an expansion with M = 21 terms. In order to construct the dual-residual estimator, we consider n train,r = 100 randomly-selected parameters, and we set tol es = 10 -4 and tol eq,r = 10 -10 in Algorithm 2. We assess performance based on n test = 20 out-of-sample parameters {µ i t } ntest i=1 . Figure 3(a) shows the average out-of-sample relative error

E avg := 1 n test ntest i=1 u hf µ i t -u µ i t u hf µ i t (37)
for various choices of the size N of the reduced-order basis Z and for two different tolerances tol eq for ρ eq . We observe that the ROM reaches a plateau that depends on the error max µ∈P H(h µ -h µ ) associated with EIM. Figure 3(b) shows the percentage of sampled elements Q/N e • 100 for several choices of N and for the two tolerances considered: note that Q grows linearly with N . Speedups of the hyper-reduced ROM compared to the ROM with HF quadrature scale with N e /Q and thus range from O(30) to O(100) for N = 1, . . . , 6 and tol eq = 10 -10 : this can be explained by observing that, for the range of N and Q considered and for serial implementations, the online cost is dominated by the cost of assembling the reduced system at each Newton iteration. Jr and of the percentage Q r /N e • 100 of sampled elements for several ROMs associated with different choices of N and two EQ tolerances tol eq . Note that J r and Q increase rapidly for moderate values of N and then reach a plateau: interestingly, this is in qualitative agreement with the behavior of the error shown in Figure 3(a). Figure 4(c) shows the relationship between the (estimated) dual residual norm and the relative error; to generate the results, we consider the ROM with N = 5 modes and two different choices of tol eq . We observe that the relative error is approximately proportional to the dual residual norm; furthermore, the prediction obtained using our approach is extremely accurate for all tests considered.

Extension to nonlinear problems

We discuss the application of the DtM framework to vector-valued nonlinear steady PDEs. Towards this end, we denote by D eq > 1 the number of equations/state variables and we introduce the FE space . (c) exact and estimated residual vs relative H 1 error. (toleq,r = 10 -10 , toles = 10 -4 )

X v T = Deq i=1 X T . Given u = [u T 1 , . . . , u T Deq ] T ∈ R N hf,v , N hf,v = N hf • D eq ,
vector-valued restriction operator

E v k : R N hf,v → R n lp •Deq such that E v k w =    E k w 1 . . . E k w Deq    , w =    w 1 . . . w Deq    , k = 1, . . . , N e .
We further introduce the vector-valued extension operator H : H → U, which is a generalization of the extension operator ( 14): we provide further comments about its definition in section 4.3. We consider problems of the form: given µ ∈ P, find u hf µ ∈ U such that

u hf µ (I dir ) = h µ , R hf µ u hf µ , v = 0 ∀ v ∈ U 0 , (38)
and

J un µ (w) i,i ,k,d,d = ∂ ∂w un i ,k,d r µ,k E v k w, e i+(d-1)n lp , Φ µ (X hf k ) ,    i, i = 1, . . . , n lp d, d = 1, . . . , D eq k = 1, . . . , N e (42b
) and then on a sparse assembly routine to define the global vector and matrix R hf µ (w), J hf µ (w) in (40a) and (40b). Algorithm 3 outlines the procedure implemented in a typical FE assembler. The function local assembler computes the tensors defined in [START_REF] Lawson | Solving least squares problems[END_REF]: note that the function takes as input (i) the restriction of the FE field to the k-th element (w un ) •,k,• , (ii) the nodes of the element Φ µ (X hf k ), (iii) a structure ref that contains quantities associated with the reference domain -quadrature rule {(ω g q , X g q )} nq q=1 , evaluations of the shape functions and their gradient in the quadrature points {X g q } nq q=1 -and (iv) the structure input data that contains inputs associated with the particular problem at hand -e.g., functional expression of the source term, coefficients. Given the local tensors R un µ and J un µ , we define the corresponding indices associated with the global enumeration (cf. Lines 4 and 6 of the Algorithm) and then we assemble the global vector and matrix by accumulation: the functions accumarray and sparse are the Matlab routines used to create vectors and sparse matrices by accumulation -analogous routines can be found in any software for scientific computing.

: J hf µ (w) = sparse ii j , jj j , J un µ (w), N hf,v , N hf,v .
To reduce the level of intrusiveness of the pMOR procedure, and ultimately simplify the implementation, we shall reuse the local integration routine local assembler to assemble ROM matrices and vectors: in the next section, we illustrate how to achieve this goal. We observe that for several problems in computational mechanics, the naive Newton's method in [START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF] might not converge: in section 4.3, we illustrate a possible remedy for the RANS equations, which is used in the numerical simulations.

Remark 1. Although we here focus on continuous FE discretizations, Algorithm 3 -and the subsequent pMOR procedure -might be generalized to finite volume and discontinuous Galerkin (DG) discretizations. More in detail, for secondorder problems, DG discretizations rely on local assembling routines that take as inputs the solution in the k-th element and in its neighbors. As a result, computation of the ROM residual requires to store the basis functions in the sampled elements and in their neighbors. Algorithm 3 might also be extended to Taylor-Hood discretizations: the implementation should take into account the fact that the number of degrees of freedom per element differs among the different variables (e.g., pressure and velocity).

Least-squares Petrov-Galerkin reduced-order model

For simplicity, we assume that h µ = ΞΘ µ with Θ : P → R M , M = O(1); then, given the extension operator H : H → U, we obtain the affine-in-parameter lifting e µ = Hh µ = WΘ µ with W = HΞ; the extension to the non-affine case can be addressed as in section 3 using EIM. Given the reduced-order bases

Z = [ζ 1 , . . . , ζ N ] ∈ R N hf,v ,N and Y = [ξ 1 , . . . , ξ J ] ∈ R N hf,v
,J with Z(I dir , :) = 0, Y(I dir , :) = 0, we seek approximations of the form

u µ = Z α µ + e µ , such that α µ ∈ arg min α∈R N sup ξ∈col(Y) R eq µ (Zα + e µ , ξ) |||ξ||| , (43a)
where R eq µ is the empirical weighted residual

R eq µ (w, v) = Ne k=1 ρ eq k r µ,k E v k w, E v k v, Φ µ (X hf k ) ; (43b)
with ρ eq ∈ R Ne + , ρ eq 0 = Q N e . Provided that ξ 1 , . . . , ξ J are orthonormal in U, we can rewrite (43a) as

α µ ∈ arg min α∈R N R N,J (α, µ) 2 , R N,J (α, µ) j = R eq µ (Zα + e µ , ξ j ), (43c)
for j = 1, . . . , J. Problem (43c) is a nonlinear least-squares problem: we can thus resort to the Gauss-Newton algorithm to efficiently compute the solution.

We rely on the dataset {(µ k , α train µ k )} ntrain k=1 (cf. (45b)) to train a regressor µ ∈ P → α (0) µ that is used to initialize the iterative Gauss-Newton scheme: we refer to [12, section 4.1] for a thorough description of the multi-target regression algorithm employed in this work.

Similarly to the linear case, we should discuss the choice of (i) the trial reduced-order basis Z, (ii) the empirical weights ρ eq , (iii) the test reducedorder basis Y, and the construction of the approximate dual residual estimator. Before addressing these tasks, we derive explicit expressions for R N,J (α, µ) and the Jacobian J N,J (α, µ) = ∇ α R N,J (α, µ).

We denote by Z un ∈ R n lp ,Ne,Deq,N and Y un ∈ R n lp ,Ne,Deq,J the unassembled trial and test reduced spaces, and we define I eq = {k ∈ {1, . . . , N e } : ρ eq k > 0}.

Then, it is easy to verify that

               R N,J (α, µ) j = k∈Ieq ρ eq k   n lp i=1 Deq d=1 (Y un ) i,k,d,j (R un ( u µ )) i,k,d   J N,J (α, µ) j,n = k∈Ieq ρ eq k   n lp i,i =1 Deq d,d =1 (Y un ) i,k,d,j (Z un ) i ,k,d ,n (J un ( u µ )) i,i ,k,d,d   (44
) for j = 1, . . . , J, n = 1, . . . , N .

Algorithm 4 summarizes the assembly procedure for the computation of the reduced residual. Note that Algorithm 4 exploits the same local integration routines employed in Algorithm 3; as in the linear case, evaluation of R N,J (α, µ) , J N,J (α, µ) requires the storage of the reduced-order basis in Q = |I eq | elements of the mesh and can thus be performed in O(Q) operations. We highlight that Algorithm 4 depends on the underlying FE model exclusively through the local assembly routine: provided that the latter is available, implementation of the ROM assembly routine is straightforward.

Algorithm 4 ROM assembly (µ, α) → R N,J (α, µ), J N,J (α, µ) . R un

µ ( u µ ) •,k,• , J un µ ( u µ ) •,•,k,•,• = local assembler u un µ •,k,• , Φ µ (X hf k )
, ref, input data . 5: end for 6: Compute R N,J (α, µ) and J N,J (α, µ) using (44).

Construction of trial and test reduced-order spaces

Given the snapshots {u hf µ k } ntrain k=1 , we apply POD to generate the reducedorder basis Z:

Z = POD {ů hf µ k } ntrain k=1 , •, • , tol pod , ůhf µ := u hf µ -e µ µ ∈ P.
On the other hand, we resort to the procedure proposed in [START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF] to construct the test space. For completeness, Algorithm 5 summarizes the data compression algorithm for the construction of the test space: the FE vector ψ k,n ∈ U 0 on Line 2 of the Algorithm is the Riesz representer of the Fréchet derivative of the residual at u hf µ k applied to the n-th element of the reduced-order basis; we refer to [START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF] for further details and for the mathematical justification of the approach for linear problems. where λ v , λ p , λ ν > 0 are suitable positive constraints. In this work, given the lifted snapshots { Ůk } ntrain k=1 , we define λ v , λ p , λ ν as the largest eigenvalues of the Gramian matrices C v , C p , C ν such that

C v i,j = d=1,2 ( Ůj d , Ůi d ), C p i,j = Ω ůj 3 ůi 3 dx, C ν i,j = (ů j 4 , ůi 4
), i, j = 1, . . . , n train .

Our choice is motivated by the need to properly "weight" variables (velocity, pressure, artificial viscosity) characterized by different magnitudes (and different units!). Due to our construction of the empirical test space, our choice of the inner product has also the effect of weighting the contributions of the residuals associated with the different equations: similarly, in [7, section IV] the authors consider a weighted residual to deal with the difference in magnitude between the different equations of the model.

Numerical results

We train the reduced-order model based on n train = 40 snapshots associated with equispaced parameters in P. In order to generate the EQ weights ρ eq , we consider additional n train,eq = 10 randomly-selected parameters (cf. comment after (45b)). The dual residual estimator is built using tol es = 10 -4 and tol eq,r = 10 -10 . We assess performance based on n test = 10 randomly-selected out-ofsample parameters.

Figure 7(a) shows the average out-of-sample error E avg (cf. ( 37)) for various choices of the size N of the reduced-order basis Z and for three different tolerances tol eq for ρ eq , the test space Y is chosen according to Algorithm 5 with J = 2N . As for the linear problem, we observe that the ROM guarantees near-optimal performance with respect to the projection error, for sufficiently tight tolerances. Figure 7(b) shows the percentage Q/N e • 100 of sampled elements with respect to N for various tolerances: for all test cases considered, Q is less than 2% of the size N e of the complete mesh. As for the linear case, speedups of the hyper-reduced ROM compared to the ROM with HF quadrature scale with N e /Q and thus range from O(50) to O(300) for N = 1, . . . , 10 and tol eq = 10 -10 . In Figure 8, we replicate the results shown in Figure 4 for the linear problem. For the test in Figure 8(c), we consider various ROMs associated with N = 3, 6, 9 and tol eq = 10 -6 , 10 -10 , 10 -14 . We observe that the dual residual is highly-correlated with the H 1 relative error: this empirical finding motivates the use of the (approximate) dual residual estimator to drive greedy sampling methods (cf. [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]). , toleq = 10 -10 , toleq = 10 -14 (b) percentage of sampled elements Q/Ne • 100 for three choices of toleq (see (a)).

Conclusions

In this work, we presented and numerically assessed a general approach for the treatment of parameterized geometries in projection-based pMOR. The approach relies on recently-developed EQ procedures to sample the HF mesh and ultimately reduce online costs. We discussed in detail the comparison with the more standard map-then-discretize approach; we presented the application to a linear diffusion problem to illustrate the key features of the methodology and analyze the offline and online computational and memory costs; then, we considered the extension to steady nonlinear PDEs and we discussed the application to the two-dimensional RANS equations.

We here resorted to continuous FE discretizations; the extension to DG and finite volume discretizations is part of ongoing research. Furthermore, we wish to combine the approach with adaptive sampling strategies to reduce offline costs: in this respect, the results in Figure 8 further motivate the use of residual-based error indicators to drive sampling. We also wish to investigate if the conservation properties of the hyper-reduced ROM derived in [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Chan | Entropy stable reduced order modeling of nonlinear conservation laws[END_REF] can be extended to problems with varying geometry, within the DtM framework. Finally, we shall combine the approach with solution-aware registration procedures, [START_REF] Taddei | A registration method for model order reduction: data compression and geometry reduction[END_REF][START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF], for the construction of the parameterized mapping Φ. and toleq = 10 -14 . (c) exact and estimated residual vs relative H 1 error. (toleq,r = 10 -10 , toles = 10 -4 , J = 2N ) grant agreement 872442 (ARIA). Tommaso Taddei also acknowledges the support of IdEx Bordeaux (projet EMERGENCE 2019). The authors thank Dr. Andrea Ferrero (Politecnico di Torino) and Professor Angelo Iollo (Inria Bordeaux) for fruitful discussions on the numerical example of section 4.

Figure 1 :

 1 Figure 1: optimal convergence and discrete bijectivity. (a) Performance of three FE discretizations for a Laplace problem. Isoparametric P3 DtM ; subparametric P1-P3 DtM , P3 MtD . (b) Failure of the discretize-then-map approach for large deformations. The mapping maps the point A in B: this makes the resulting deformed mesh singular.

3 :

 3 [η 1 , . . . , η Jr ] = POD ({ψ r µ k r } k , (•, •), tol es ). 4: Compute G r , b r in (32).

Figure 2 :

 2 Figure 2: Model problem. (a) coarse-grained partition associated with Φ, (b) computational mesh.

Figure 3 :

 3 Figure 3: Potential flow past an airfoil. Performance with respect to N and toleq. (a) average out-of-sample error Eavg: projection , Galerkin ROM with HF quadrature , toleq = 10 -6, toleq = 10 -10 . (b) number of sampled elements Q for two choices of toleq (see (a)).

Figure 4

 4 Figure 4 investigates the performance of the error estimation procedure. Figures 4(a) and (b) show the behavior of the size J r of the empirical test space Y rJr and of the percentage Q r /N e • 100 of sampled elements for several ROMs associated with different choices of N and two EQ tolerances tol eq . Note that J r and Q increase rapidly for moderate values of N and then reach a plateau: interestingly, this is in qualitative agreement with the behavior of the error shown in Figure3(a). Figure4(c)shows the relationship between the (estimated) dual residual norm and the relative error; to generate the results, we consider the ROM with N = 5 modes and two different choices of tol eq . We observe that the relative error is approximately proportional to the dual residual norm; furthermore, the prediction obtained using our approach is extremely accurate for all tests considered.

Figure 4 :

 4 Figure 4: Potential flow past an airfoil. Dual residual norm estimation. (a) size of the empirical test space Jr with respect to N , for toleq = 10 -10 and toleq = 10 -12 . (b) percentage of sampled elements Qr with respect to N , for toleq = 10 -10 toleq = 10 -12. (c) exact and estimated residual vs relative H 1 error. (toleq,r = 10 -10 , toles = 10 -4 )

Algorithm 3

 3 High-fidelity assembly. (µ, w) → R hf µ (w), J hf µ (w) .

1 : 4 : 6 :

 146 for k = 1, . . . , N e do 2:R un µ (w) •,k,• , J un µ (w) •,•,k,•,• = local assembler (w un ) •,k,• , Φ µ (X hf k ), ref, input data 3: end for Define ii r i,k,d = T i,k + N hf (d -1) for i = 1, . . . , n lp , k = 1, . . . , N e , d = 1, . . . , D eq . 5: R hf µ (w) = accumarray ii r , R un µ (w), N hf,v , 1 . Define ii j i,i ,k,d,d = T i,k + N hf (d -1) and jj j i,i ,k,d,d = T i ,k + N hf (d-1) for i, i = 1, . . . , n lp , k = 1, . . . , N e , d, d = 1, . . . , D eq .

7

 7

1 :

 1 Compute the deformed nodes {Φ µ (x hf i,k ) : i = 1, . . . , n lp , k ∈ I eq }.

2 :

 2 Compute u un µ i,k,d = n (Z un ) i,k,d,n (α) n , i = 1, . . . , n lp , k ∈ I eq , d = 1, . . . , D eq . 3: for k ∈ I eq do 4:

Figure 5 :Figure 6 :

 56 Figure 5: Ahmed's body. (a) description of the body geometry. (b) coarse-grained partition associated with Φ. (c) computational mesh. L = 1.044, R = 0.1, h b = 0.05, Hc = 0.288, δ = 0.222.

Figure 7 :

 7 Figure 7: Ahmed's body problem. Performance with respect to N and toleq (J = 2N ). (a) average out-of-sample error Eavg: projection , AMR ROM with HF quadrature , toleq = 10 -6, toleq = 10 -10 , toleq = 10 -14 (b) percentage of sampled elements Q/Ne • 100 for three choices of toleq (see (a)).

Figure 8 :

 8 Figure 8: Ahmed's body problem. Dual residual norm estimation. (a) size of the empirical test space Jr with respect to N , for toleq = 10 -10 and toleq = 10 -14 . (b) percentage of sampled elements Qr/Ne •100 with respect to N , for toleq = 10 -10and toleq = 10 -14 . (c) exact and estimated residual vs relative H 1 error. (toleq,r = 10 -10 , toles = 10 -4 , J = 2N )

We here assume that quadrature weights and quadrature points in D and shape functions are computed in advance. Exploiting fast matrix matrix multiplication algorithms, we might be able to assemble the tensor L ∇,Φ in O(n β lp ) flops with

< β <

Computation of {A k,Φµ } k : O(n

lp Qr); computation of E k uµ (given αµ): O(n lp Qr(N + M )); computation of R r k,µ : O(Qr(n 2 lp + Jrn lp )); computation of R r µ 2 : O((Qr + 1)Jr).

The same result holds for the nonlinear case.

Acknowledgements

The authors acknowledge the support by European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Actions,

Here, we propose to adapt the strategy in [START_REF] Taddei | An offline/online procedure for dual norm calculations of parameterized functionals: empirical quadrature and empirical test spaces[END_REF] to devise a rapid and reliable estimator of the dual residual norm [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF]: the approach involves the definition of a low-dimensional empirical test space and of an empirical quadrature rule.

We approximate the dual residual norm R hf µ ( u µ ) for any given µ as

where Y r Jr = span{η j } Jr j=1 ⊂ U 0 is referred to as empirical test space and R eq,r

(28) Similarly to [START_REF] Grimberg | On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows[END_REF], ρ eq,r = [ρ eq 1 , . . . , ρ eq Neq ] T is a sparse vector of positive weights such that ρ eq,r 0 = Q r N e ; we denote by I eq,r = {k ∈ {1, . . . , N e } : ρ eq,r k > 0} the set of indices associated with the sampled elements. Provided that {η j } Jr j=1 is an orthonormal basis of Y r Jr , we find that

where { R r k,µ } k ⊂ R Jr are given by R r k,µ j

) for j = 1, . . . , J r and k = 1, . . . , N e . Given µ ∈ P and the estimate α µ , the cost of estimating R eq µ is dominated by the cost of assembling the local residuals { R r k,µ } k and scales with 2 O Q r n 3 lp + n lp (N + M ) + J r n lp . In the remainder of this section, we illustrate how to construct the test space Y r Jr and the quadrature weights ρ eq,r .

We shall choose Y r Jr to approximate the manifold of Riesz representers ⊂ P, we compute the solution to the ROM u µ for all µ ∈ P train,r and then the Riesz representation of the residual ψ r µ using [START_REF] Ali | Reduced basis methods with adaptive snapshot computations[END_REF]; finally, we use POD to construct the test space Y r Jr . On the other hand, we choose the weights ρ eq,r so that (i) ρ eq,r is as sparse as possible, (ii) all weights are non-negative, and (iii)

where R hf : U × U 0 × P → R is nonlinear in the first (and possibly the third) argument and is associated with the system of D eq > 1 PDEs of interest. Note that, for continuous discretizations of second-order elliptic PDEs, we might write the residual R hf as the sum of local residuals associated with integration over the elements of the HF discretization:

As in section 3, we denote by u hf µ ∈ X T (and u hf µ,Φµ ∈ X Φµ(T ) ) the HF FE solution field associated with the parameter µ ∈ P.

In section 4.1, we present the HF solver: we introduce relevant notation and we highlight the building blocks of the full-order-model implementation. In section 4.2, we introduce the projection-based LSPG ROM. In section 4.3 we introduce the model problem and the geometric parameterization; finally, in section 4.4, we present results of the numerical experiments.

Implementation of the high-fidelity solver

We introduce the algebraic counterpart of R hf ,

where e 1 , . . . , e N hf,v are the vectors of the canonical basis in R N hf,v , and its Jacobian J hf : R N hf,v × P → R N hf,v ,N hf,v such that

) Exploiting the previous definitions, we can rewrite [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] as

with I in = {1, . . . , N hf,v } \ I dir . Problem (40c) can be tackled using variants of the Newton's method: for k = 1, . . . until convergence,

with J hf,(k)

. This implies that a typical HF solver for [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] involves (i) an efficient sparse assembly routine to compute R hf,(k) µ , J hf,(k) µ at each iteration, and (ii) a direct or iterative sparse linear solver to solve (41) at each iteration k.

Exploiting [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF], given µ ∈ P and w ∈ U, the assembly routine relies on local integration to build the tensors R un µ , J un µ such that

Algorithm 5 Construction of the test space Y J = (col(Y), •, • ).

Empirical quadrature

We apply a variant of the empirical quadrature procedure considered in [START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF]. More precisely, given the snapshots {u hf µ k } ntrain k=1 and the reduced-order bases Z and Y, we seek ρ eq ∈ R Ne

+ such that (i) ρ eq 0 , is as small as possible, (ii) the constant function is accurately approximated in Ω (see [START_REF] Stabile | A reduced order variational multiscale approach for turbulent flows[END_REF]), and (iii) for all µ ∈ P train,eq = P train ∪ {µ eq,i } ntrain,eq i=1 ⊂ P,

where R hf N,J : R N × P → R J is the reduced residual associated with the choice ρ eq k = 1 for k = 1, . . . , N e (reduced residual with HF quadrature), while α train µ is chosen such that

Note that [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] involves the solution to n train,eq ROMs with HF quadrature. By combining [START_REF] Stabile | A reduced order variational multiscale approach for turbulent flows[END_REF] with [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF], we obtain a sparse representation problem of the form [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF], which can be tackled using the non-negative least-squares method (cf. section 3.3.2). Some comments are in order. First, in [START_REF] Taddei | Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs[END_REF], we consider a slightly different accuracy constraint in which we pre-multiply R N,J (α train µ , µ) -R hf N,J (α train µ , µ) by the transposed reduced Jacobian: in our experience, the constraint (45) performs better for moderate values of n train . Second, to improve performance of EQ, similarly to [13, Algorithm 1], we add training points to P train by solving the ROM with HF quadrature: we empirically observe that this choice improves performance of the hyper-reduced ROM, particularly for moderate values of n train . In addition, since the solution to (45b) 2 is significantly less expensive than an HF solve, increasing the size of P train,eq does not significantly increase offline costs.

A posteriori error estimation

As for the linear case, we rely on an approximate dual residual estimator:

where the space Y r Jr ⊂ U 0 and the weights ρ eq,r ∈ R Ne + are computed by extending the approach presented in section 3.3.3 to nonlinear problems -since the extension is straightforward, we omit the details. As discussed in the pMOR literature, rigorous and rapidly-computable a posteriori error estimators are rarely available for nonlinear PDEs: in the numerical experiments, we verify numerically that the residual estimator is highly-correlated with the reconstruction error and can thus be used to drive an adaptive sampling strategy for the optimal selection of the training parameters {µ k } ntrain k=1 .

Model problem: RANS simulations of the flow past the Ahmed body

We consider a two-dimensional flow past a parameterized Ahmed body, [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF][START_REF] Lienhart | Flow and turbulence structures in the wake of a simplified car model (Ahmed model)[END_REF], at moderate Reynolds number. The geometry of the body is prescribed in Figure 5(a). The problem is parameterized with respect to the slant angle µ: in the numerical simulations, we consider angles between 5 o and 50 o degrees, P = [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Tezduyar | Finite element stabilization parameters computed from element matrices and vectors[END_REF]. Reynolds number is defined as the ratio between horizontal inflow velocity times the height H c of the body divided by kinematic viscosity; in our tests, we consider Re = 3 • 10 3 .

We consider the incompressible RANS equations with Spalart-Allmaras closure model (see [START_REF] Spalart | One-equation turbulence model for aerodynamic flows[END_REF][START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF]). We denote by U the set of D eq = 4 state variables which comprises the two velocity components, pressure and turbulent viscosity. We refer to Appendix A for further details concerning the mathematical model and boundary conditions.

We discretize the problem using a SUPG FE method, [START_REF] Tezduyar | Stabilized finite element formulations for incompressible flow computations[END_REF][START_REF] Tezduyar | Finite element stabilization parameters computed from element matrices and vectors[END_REF]. We add the Pressure-Stabilized Petrov-Galerkin (PSPG) term to eliminate spurious modes in the pressure solution when considering the same polynomial order for pressure and velocity; furthermore, we consider the least-squares incompressibility constraint (LSIC) stabilization term to improve accuracy and conditioning of the discrete problem, [START_REF] Franca | Stabilized finite element methods: II. The incompressible Navier-Stokes equations[END_REF][START_REF] Gelhard | Stabilized finite element schemes with LBB-stable elements for incompressible flows[END_REF][START_REF] Braack | Stabilized finite element methods for the generalized oseen problem[END_REF]. In all our tests, we consider the mesh depicted in Figure 5(c) that consists of N e = 26602 P2 elements (N hf,v = 221960). To compute the HF solution, we resort to a pseudo-time continuation (see, e.g., [START_REF] Kelley | Convergence analysis of pseudo-transient continuation[END_REF]) procedure to determine an accurate initial condition for the Newton solver; then, we apply the Newton scheme [START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF] with 1D line search.

Similarly to section 3.4, we partition the domain Ω µ into the seven regions {Ω i,µ } 7 i=1 depicted in Figure 5(b), we define parameterized Gordon-Hall maps in each of the seven regions, and we define Φ as in [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF]. Note that Φ is the identity in the regions 1, 2, 6, 7.

Figure 6(a) shows the behavior of U 1,µ,Φµ (horizontal velocity) for µ = 25 o , while Figure 6(b) shows the profile of U 1,Φµ for x 1 = 1.1 and x 2 ∈ [0, 0.4] for three values of µ: we observe the presence of recirculating flow in the proximity of the body that becomes more prominent for larger slant angles.

Since the boundary conditions are parameter-independent, we bypass the definition of the extension operator and we simply define e µ := U hf μ with μ = 25 o . We equip the space U with the inner product

Appendix A. Explicit expressions of the RANS equations

We denote by u the velocity field, by p the fluid pressure, and by ν the turbulent viscosity; we further denote by U = [u, p, ν] the vector of state variables; 31 then, we introduce the RANS equations as

where

where ν > 0 is the kinematic viscosity, P, D are the wall production and destruction terms which are explicit functions of vorticity and turbulent viscosity (see [START_REF] Spalart | One-equation turbulence model for aerodynamic flows[END_REF][START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF]). We are here interested in the limit solution for t → ∞. FE discretizations of (A.1) can be obtained by projecting the equations onto finite dimensional FE spaces. For stability reasons, we add SUPG and LSIC stabilization forms to the momentum equation

where τ supg and τ lsic are defined as in [START_REF] Peterson | Overview of the incompressible Navier-Stokes simulation capabilities in the MOOSE framework[END_REF] and [START_REF] Braack | Stabilized finite element methods for the generalized oseen problem[END_REF], respectively. We also add the PSPG term to the continuity equation

with τ pspg = τ supg . Finally, we add the SUPG form to the transport equation for the turbulent viscosity:

In the implementation, we resort to pseudo-time integration to solve the HF system and we consider the strategy in [START_REF] Crivellini | A Spalart-Allmaras turbulence model implementation in a discontinuous Galerkin solver for incompressible flows[END_REF] to treat negative viscosities that might be generated during the iterations. Finally, boundary conditions are set as follows:

• Top and outlet: homogeneous Neumann conditions.

• Bottom and body: u = 0, ν = 0.

Appendix B. Additional remarks

We provide below further interpretations and remarks.

Remark 2. Choice of the boundary norm in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF]. Due to our choice of • dir , the space H is a discrete counterpart of L 2 (Γ dir ); for second-order elliptic PDEs, this choice is not consistent with the infinite-dimensional formulation, which involves the fractional Sobolev space H 1/2 (Γ dir ). A consistent choice of • dir would be • dir = H • : clearly, evaluations of the latter require the solution to a HF problem of size N hf . Our choice is thus a compromise between mathematical rigor and simplicity of implementation.

Remark 3. Reduced integration domain. Following [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF][START_REF] Fritzen | An algorithmic comparison of the hyper-reduction and the discrete empirical interpolation method for a nonlinear thermal problem[END_REF], we might refer to the region k∈Ieq D k as to reduced integration domain (RID). Note that, while the EQ procedure employed in this paper directly reduces the integration domain without modifying the test functions, the hyper-reduction approach in [START_REF] Fritzen | An algorithmic comparison of the hyper-reduction and the discrete empirical interpolation method for a nonlinear thermal problem[END_REF] similarly to [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF] -ultimately performs projection of the residual on the nodes of the mesh: as a result, for the linear problem 3 considered in section 3, the latter enforces that if R µ (Zα + e µ , e i , Φ µ ) = 0 ∀ i = 1, . . . , N hf then α = α µ , which is in general false for the EQ approach considered here. We further remark that the strategy employed in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF][START_REF] Fritzen | An algorithmic comparison of the hyper-reduction and the discrete empirical interpolation method for a nonlinear thermal problem[END_REF] to construct the RID is different, and the authors do not explicitly discuss the problem of parameterized geometries. Remark 4. Connection between linear and nonlinear procedures. From the discussion in section 3.2, we deduce that the FE solver for (16) shall include a routine of the form:

where F un k,Φ , A un k,Φ are the contributions to the global system associated with the k-th element, Φ(X hf k ) are the nodes of the deformed mesh associated with the k-th element, ref and input data are data structures introduced in section 4.1 -note that F un k,Φ ≡ 0 in the example of section 3. Then, recalling (13), we find that (21b) could be rewritten as

The latter identity is the linear (and scalar) counterpart of the general relation provided in [START_REF] Gordon | Construction of curvilinear co-ordinate systems and applications to mesh generation[END_REF].

Remark 5. Positivity of the empirical weights. Positivity of the weights ρ eq in ( 20) is justified by the analysis in [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF]; on the other hand, the constraint ρ eq,r ≥ 0 is not strictly necessary and might be omitted. We might thus construct the weights ρ eq,r by approximating the solution to the sparse representation problem min

To apply nonnegative least-squares to (B.1), we first observe that (see the discussion in [33, section 2.3.1]), given any solution (ρ eq,r + , ρ eq,r -) to min ρ + ,ρ -∈R Ne ρ + 0 + ρ -0 , s.t. G r (ρ + -ρ -) -b r ≤ δ, ρ + , ρ -≥ 0, (B.2) ρ eq,r = ρ eq,r + -ρ eq,r -solves (B.1). Problem (B.2) can be tackled using nonnegative least-squares; note, however, that (B.2) has twice as many degrees of freedom as [START_REF] Yano | A reduced basis method for coercive equations with an exact solution certificate and spatio-parameter adaptivity: energy-norm and output error bounds[END_REF] and is thus more challenging to be approximated.

Remark 6. Mathematical justification of the choice of Y r

Jr , ρ eq,r in section 3.3.3. We denote by Π Y r Jr : U → Y r Jr the orthogonal projection operator onto Y r Jr . We can then exploit the same argument as in [START_REF] Taddei | An offline/online procedure for dual norm calculations of parameterized functionals: empirical quadrature and empirical test spaces[END_REF]Proposition 2] to obtain

.

(B.3) We observe that (I) is associated with the approximation properties of the empirical test space for the manifold M r test , while (II) -which is (31) -is related to the accuracy of the quadrature rule.

Graphical Abstract

A discretize-then-map approach for the treatment of parameterized geometries in model order reduction Tommaso Taddei, Lei Zhang