{'fraction_non_alphanumeric': 0.04259298618490967, 'fraction_numerical': 0.01390010626992561, 'mean_word_length': 4.40816091954023, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical Model Reduction Techniques for Flow Modeling in a Parametrized Setting

In this work we focus on two different methods to deal with parametrized partial differential equations in an efficient and accurate way. Starting from high fidelity approximations built via the hierarchical model reduction discretization, we consider two approaches, both based on a projection model reduction technique. The two methods differ for the algorithm employed during the construction of the reduced basis. In particular, the former employs the proper orthogonal decomposition, while the latter relies on a greedy algorithm according to the certified reduced basis technique. The two approaches are preliminarily compared on two-dimensional scalar and vector test cases.


Introduction
The interest for fluid dynamics simulations is growing more and more in the scientific community and in the nowadays society, both in terms of spread and of relevance. This is, most of all, due to practical issues and time reasons. Physical experiments are often very expensive and time demanding so that, for some specific applications (e.g., in naval or aeronautic applications as well as in medical surgery planning), they are not well-suited, and numerical simulations become the actual tool for modeling reliable scenarios in such contexts.
Although the computational power is continuously growing, standard methods in Computational Fluid Dynamics (CFD), such as direct numerical simulations based on finite elements, may be very demanding in terms of computational time and numerical sources, especially when interested in simulating challenging phenomena in complex domains with a certain accuracy, or, even more, when dealing with multi-query or parametric frameworks [24].
For these reasons, many different methods have been proposed in the scientific panorama with the aim of offering a compromise between modeling accuracy and computational efficiency. Model reduction techniques represent a relevant solution in such a direction [24]. Some of them are strictly intertwined with the model of the phenomenon at hand, while others perform the reduction only under specific physical assumptions on the described configuration [36,9,8].
In this work, we focus on the Hierarchical Model (HiMod) reduction technique [17,31,30,34]. This procedure has been devised to describe CFD configurations where a principal dynamics overwhelms the transverse ones, with a strong interest for hemodynamic configurations [5,20,32]. The leading dynamics is aligned with the main stream of the flow, while transverse dynamics are generally induced by geometric irregularities in the computational domain and play a role only in localized areas. In practice, the idea is to discretize the different dynamics by resorting to different numerical methods, in the spirit of a separation of variable. For instance, in the original proposal of HiMod reduction, the main direction of the flux is discretized with 1D finite elements, while the transverse dynamics are reconstructed by using few degrees of freedom via a suitable modal basis. This separate discretization, independently of the dimension of the (full) problem at hand, leads to solve a system of coupled 1D problems, whose coefficients include the effect of the transverse dynamics. This ensures to the HiMod reduction a reliability which is considerably higher compared with standard 1D reduced models, and at a computational cost which remains absolutely affordable. The computational advantages provided by a HiMod discretization have been also verified for simulations in real geometries [3,20,27]. In particular, HiMod reduction guarantees a linear dependence of the computational cost on the number of degrees of freedom in contrast to a standard full finite element model which demands a suitable power of such a number.
The interest of this paper is on a parametric setting, where the reference model, coinciding with a parametric Partial Differential Equation, has to be solved several times, for many different values of the parameter. The goal we pursue is to approximate, for each value of the parameter, the HiMod discretization by a modeling procedure which turns out to be computationally cheaper than HiMod reduction itself. A first effort in such a direction is proposed in [6,26]. The authors apply a Proper Orthogonal Decomposition procedure to HiMod approximations, to extract a reduced basis which allows us to predict the HiMod discretization associated with any value of the parameter. The new procedure, named HiPOD, is numerically investigated on scalar elliptic problems and on the Stokes equations in [6]. In this paper, we investigate a new procedure, alternative to HiPOD, to pursue the same goal of managing, in a cheap way, a parametric framework. In particular, we aim at exploiting the computational advantages provided by a greedy algorithm in the construction of a reduced basis [16,22]. For this purpose, we combine HiMod reduction with the Reduced Basis (RB) approach [37,21], into the new technique called HiRB. HiPOD or HiRB approximations considerably decrease the computational effort due to the lower dimension of the high fidelity problems. According to an offline/online paradigm, the offline stage remains the bottleneck from a practical viewpoint. However, the employment of HiMod discretizations as high fidelity solutions significantly reduces the computational effort of this phase. Finally, a system of very small order is solved during the online phase and yields a reliable approximation for the parametric problem at hand.
The paper is organized as follows. Section 2 introduces the HiMod setting, and particularizes such a procedure both to a scalar advection-diffusion-reaction problem and to the Stokes equations. Sections 3 and 4 exemplify the HiPOD and the HiRB procedures, respectively on the test problems introduced in the previous section. Particular care is devoted to the inf-sup condition characterizing the discretization of the Stokes problem. Actually, although the high fidelity solutions satisfy the Ladyzhenskaya-Brezzi-Babuška (LBB) condition, this is not ensured either by the POD and the reduced basis formulations. To overcome this issue, we propose here to resort to the supremizer enrichment stabilization technique [4]. Section 5 performs a preliminary comparison between HiPOD and HiRB, starting from the (two-dimensional) test cases considered throughout the paper. Finally, some conclusions are drawn in Section 6 and future developments are summarily itemized.

The HiMod setting
We summarize here the main features of the Hierarchical Model (HiMod) reduction technique, following the original setting in [17,31,34,30]. We assume that the d-dimensional domain, Ω, with d = 2, 3, coincides with the fiber bundle where Ω 1D is the supporting fiber aligned with the main flow, while γ x denotes the (d − 1)dimensional transverse fiber at point x ∈ Ω 1D , parallel to the secondary dynamics. In practice, computations are performed in a reference domain,Ω, so that Ψ(Ω) =Ω, Ψ : Ω →Ω being a sufficiently regular map (see Figure 1). In general, domainΩ coincides with a rectangle (d = 2) or with a right circular cylinder (d = 3). For simplicity, we consider a rectilinear axis Ω 1D = (0, L) with L > 0, so that Ψ(x, y) = (x, ψ x (y)), for any (x, y) ∈ Ω. Map Ψ preserves the supporting fiber and deforms only the transverse shape of the domain via the map ψ x : γ x →γ between the generic, γ x , and the reference,γ, transverse section. This induces a decomposition similar to (1) on the reference domain as well, beingΩ = x∈Ω 1D {x} ×γ. We refer to [29,32,5] for the more general case of a curvilinear fiber Ω 1D .
In the next sections, we apply the HiMod discretization to a scalar and to a vector problem, in order to detail the involved procedures.
HiMod reduction applies to the weak form of the full problem, where, to simplify notation, we assume g = 0 in (2) and we drop the dependence on (x, y). The assumptions above on the problem data ensure the well-posedness of (3) [15]. Thus, the HiMod formulation for problem (2) can be stated as for a certain m ∈ N + , and with k=1 is a modal basis of functions defined onγ, orthonormal with respect to the L 2 (γ)-scalar product, and m is the modal index, i.e., the number of modes employed to model the transverse dynamics. In the sequel, we identify V h 1D with the standard space of the (continuous) finite elements (FE) [15], while referring to [29,32,5] for different discretizations. As far as the choice of m is concerned, it can be fixed a priori, thanks to heuristic considerations or to a partial knowledge of the full problem, or a posteriori, driven by a modeling error analysis as in [33,35].
The HiMod space has to be endowed with a conformity and a spectral approximability assumption to ensure the well-posedness of formulation (5), as well as a standard density hypothesis has to be advanced on the discrete space V h 1D to guarantee the convergence of the HiMod approximation, u m , to u (we refer to [31] for the details).
Concerning the boundary conditions completing problem (2), we have to distinguish between data assigned on the inflow/outflow boundaries and on the lateral surface of Ω. In the first case, we employ a modal expansion of the data to be imposed in an essential way. With reference to lateral boundary conditions, we resort to the approach proposed in [3], where the authors set a general way to incorporate, essentially, the lateral boundary data by defining a customized basis referred to as an educated modal basis. The effectiveness of such a procedure is successfully investigated both from a theoretical and a numerical point of view in the same work. In the numerical assessment below, we resort to an educated modal basis to manage the lateral boundary conditions. From a computational viewpoint, discretization (5) turns the full model (2) into a system of m coupled 1D problems defined on Ω 1D . This represents the strength-point of a HiMod formulation due to the expected saving in terms of computational effort, for m reasonably small. Actually, we are led to solve the HiMod linear system where A m ∈ R mN h ×mN h and f m ∈ R mN h are the HiMod stiffness matrix and right-hand side associated with the bilinear and linear forms in (4), with N h = dim(V h 1D ), and where u m ∈ R mN h collects the (unknown) coefficients of the HiMod expansion For a full characterization of system (7), we refer to [17,31]. To qualitatively investigate the reliability of the HiMod reduction, we solve problem (2) by means of a FE solver and of a HiMod discretization on the two-dimensional (2D) domain, Ω, identified by the map ψ x (y) = y − 0.  Γ D and Γ N are identified with the inflow and outflow boundary, respectively while Γ R coincides with the lateral surface, being g = h = l = 0 and ρ = 1. The FE solver employs affine finite elements on a 2D unstructured mesh consisting of 12800 triangles. The HiMod reduction discretizes the main stream with linear finite elements associated with a uniform subdivision of the supporting fiber into 80 intervals, while resorting to m = 8 educated modal basis functions in the transverse direction. Figure 2 compares the FE with the HiMod approximation and highlights the good qualitative matching between the two solutions. A quantitative investigation of the HiMod procedure is beyond the goal of this paper and can be found, e.g., in [31,3,20], together with a modeling convergence analysis both with respect to the modal expansion and the FE discretization.

2.2.
HiMod reduction for the Stokes equations. In this section we generalize the HiMod procedure to a vector problem, namely, to the Stokes equations, find u : Ω → R d and p : Ω → R such that in Ω, where u = (u x , u y ) T and p denote the velocity and the pressure of the flow, ν > 0 is the kinematic viscosity, D(u) = 0.5 ∇u + (∇u) T is the strain rate, f = [f x , f y ] T is the force per unit mass, and C in and C out are the inflow and outflow data, respectively. The boundary ∂Ω is partitioned so that the inflow and outflow sections, Γ in and Γ out , coincide with the fibers γ 0 and γ L , respectively while Γ w denotes the lateral walls x∈Ω 1D {x} × ∂γ x . On Γ in and Γ out we impose a non-homogeneous tangential Neumann condition, with C in and C out constant values, and we assume the transverse component of the velocity to be null. Finally, a no-slip boundary condition is enforced on the velocity along the wall surface. By introducing the bilinear forms a(·, ·) : and the functional F (·) : V → R, given by , the weak form of (9) can be stated as find u ∈ V , p ∈ Q such that where the notation is simplified by removing the dependence on (x, y) and where the natural boundary conditions still have to be properly included in F . The generalization of a HiMod reduction to the Stokes equations deserves particular attention, especially with reference to the two-field formulation involved by saddle point problems [11,12]. While the search for inf-sup (or LBB) compatible spaces for velocity and pressure is largely investigated for standard finite element and spectral discretizations [11,12,14], we are not aware of any theoretical result for hybrid methods involving both the techniques. In [20,27], empirical criteria to select the HiMod velocity and pressure are provided and numerically checked. A first theoretical assessment of these criteria is currently under investigation [7]. The reduced spaces involved in the HiMod discretization of the Stokes equations are for the velocity and the pressure, respectively where space V mu is the scalar space defined as in (6). Notice that we employ the same (educated) modal basis, {ϕ k } mu k=1 , for all the components of the HiMod velocity, while we resort to the modal basis {η k } mp k=1 to discretize the pressure. For what concerns the compatibility of the velocity with the pressure HiMod spaces, in this work we adopt the empirical criterion in [3,20], so that we set m u = m p + 2 and we choose the 1D FE pair (V h 1D , Q h 1D ) as the Taylor-Hood P2/P1 elements [11]. We denote by N h,u and N h,p the dimension of V h 1D and Q h 1D , so that the dimension of V mu and Q mp becomes dm u N h,u and m p N h,p , respectively.
Let us now describe the algebraic formulation for the HiMod discretization of the Stokes problem. After assembling the matrices A mu ∈ R dmuN h,u ×dmuN h,u , B mp,mu ∈ R mpN h,p ×dmuN h,u and the vector f mu ∈ R dmuN h,u associated with the HiMod discretization of the forms a(u, v), b(u, q) and F (v) in (10), the linear system has to be solved, where u mu ∈ R dmuN h,u and p mp ∈ R mpN h,p collect the unknown coefficients of the HiMod expansion for the velocity, u mu , and the pressure, p mp , respectively and with 0 the null vector in R mpN h,p . We conclude this section by exemplifying the HiMod procedure on a benchmark Stokes test case. The reference domain isΩ = (0, L) × (−0.5, 0.5) ⊂ R 2 , while the map ψ x is given by , so that Ω coincides with a sinusoidal domain. In particular, we select L = 6 and H = 1. Furthermore, we assign ν = 5, f = [3, 0] T , C in = 10, C out = 0.
Concerning the HiMod discretization, we enrich the Taylor-Hood P2/P1 discretization of the mainstream by resorting to m p = 5 and m u = 7 educated modes to discretize the transverse components of the pressure and the velocity, respectively. In particular, both the finite element approximations rely on a uniform subdivision of the supporting fiber into 80 subintervals.    Remark 2.1. As investigated more into the details in [42], the conservative form (9) of the Stokes problem allows one to obtain a more accurate HiMod approximation with respect to a nonconservative formulation. This is due to the coupling between the velocity components (namely, between the off-diagonal blocks of the HiMod matrix) ensured by the conservative form.
Remark 2.2. Cylindrical domains demand a careful selection of the modal basis as investigated in [20], where a polar coordinate system is employed to model the transverse dynamics in circular and elliptic pipes. As a first alternative, one can resort to the Transversally Enriched Pipe Element Method (TEPEM) [27]. In such a case, the physical domain is mapped to a reference slab, so that the modal basis coincides with the tensor product of the one-dimensional modal functions. In [20] TEPEM is compared with the HiMod reduction based on polar coordinates to highlight pros and cons of the two approaches. As expected, TEPEM turns out to be easier to implement but less accurate than HiMod, especially in the presence of highly oscillatory flows. The higher reliability characterizing the HiMod scheme can be ascribed to the tight correspondence between the domain geometry and the modal basis. Another alternative to a polar coordinate system is represented by the isogeometric version of the HiMod approach, as recently investigated in [5] for patient-specific geometries.

The HiPOD approach
HiMod reduction allows one to recast a d-dimensional problem as a system of 1D problems. Even though this leads to a computational benefit, the overall computational cost might be still not negligible when dealing with multi-query or inverse problems or, more in general, with parametrized settings. In this context, to further lighten the computational effort, we rely on projection-based model reduction techniques, by properly combining the HiMod discretization with the Proper Orthogonal Decomposition (POD). In particular, we adopt the standard offline/online paradigm [21]. The offline phase is meant to build the POD basis, starting from the hierarchical reduction of a certain number of full problems associated with a sampling of the parameter domain. The online phase approximates the HiMod discretization for any new value of the parameter by employing the POD basis. The combination between HiMod and POD justifies the name of this method, HiPOD [6,26,28].
3.1. HiPOD reduction for ADR problems. We generalize problem (3) to a parameter dependent setting, so that the new problem is (12) find where µ ∈ D ⊂ R P denotes a vector of P real numbers collecting the problem parameters, and D is the parameter domain.
For each value µ (j) ∈ S, we approximate the corresponding solution, u(µ j ), to (12) by computing the HiMod discretization, u m (µ (j) ), for a certain value, m, of the modal index. According to the modal expansion in (8), this yields the M vectors collecting, by mode, the HiMod coefficients. These vectors are employed to assemble the response matrix which will be used to extract the POD basis. To this aim, we define the correlation matrix associated with U , with X m,u ∈ R mN h ×mN h the HiMod matrix associated to the inner product in V m . Then, we consider the spectral decomposition of matrix C, so that The POD basis is thus identified by the vectors with N ≤ M . Integer N can be selected driven by heuristic considerations (e.g., by studying the trend of the spectrum of C) or by an energy criterion, for instance, we pick N such that and ε a user-defined tolerance. Independently of the adopted criterion, we denote the reduced POD space by V m,N = span{ϕ 1 , . . . , ϕ N }, and the matrix collecting, by column, the POD basis Remark 3.1. As an alternative to the approach based on the correlation matrix, one can exploit directly the spectral properties of the response matrix U to extract the reduced POD basis, by setting X m,u = I in (13), with I ∈ R mN h ×mN h the identity matrix [6]. This two-fold possibility is justified by the relation between the singular vectors of U and the eigenvectors of C [19].
3.1.2. The online phase. Goal of the online phase is to build a HiMod approximation to problem (12) for any value µ ∈ D of the parameter, by skipping the solution of the associated HiMod system (7), where the dependence on the parameter µ has been highlighted. This task is accomplished by means of a projection step, i.e., by solving the system Notice that the order of system (17) is significantly smaller compared with the HiMod system (16), being, in general, N mN h . Successively, u m,N (µ) is projected back to the original HiMod space, thus obtaining the approximation In the sequel, we will denote by u m,N,POD (µ) the HiPOD approximation for the HiMod solution u m (µ) associated with vector u m,N,POD (µ). As known, the bottleneck of the projection approach lies in the assembly of A m (µ) and f m (µ). An efficient assembly can be obtained under an affine parameter dependence hypothesis. This requirement will be accomplished in the considered test cases. Alternative procedures, such as the empirical interpolation method, are adopted in more complex cases [21].
3.1.3. Numerical assessment. We apply the HiPOD procedure to the test case in Section 2.1. We remind that the high fidelity solution is provided, in such a case, by a HiMod approximation.
During the offline phase, we hierarchically reduce the corresponding 100 ADR problems, by employing the same HiMod discretization as in Figure 2, right. The number N of POD basis functions is picked by analyzing the spectrum of the correlation matrix C (see Figure 6). The eigenvalues quickly decrease. For the sake of comparison, we select N = 20 for both ranges. The corresponding eigenvalue, normalized to the maximum one, is O(10 −6 ) and O(10 −7 ) for D 1 and D 2 , respectively.
Then, we run the online phase to approximate the HiMod solution associated with the parameter µ = [5,20,75,25] T , i.e., the solution in Figure 2, right.
By comparing the contour plot of the HiPOD approximation in Figure 7 (left for D 1 , right for D 2 ) with the HiMod discretization in Figure 2, right, we recognize that the global trend of the HiMod solution is correctly detected by both the HiPOD solutions. As shown in Figure 8, a more quantitative investigation based on the distribution of the error u 8 (µ) − u 8,20,POD (µ) shows that the solution associated with the smallest parameter domain is, as expected, more accurate (about of two orders of magnitude) with respect to the approximation obtained when dealing with D 1 . The highest accuracy is particularly evident in correspondence with the outflow boundary. In Section 5 we provide a further error analysis for this test case, based on a random sampling.
3.2.1. The offline phase. As in Section 3.1.1, we assume to deal with a training set S consisting of M different parameters. According to a segregated approach, we generate the POD basis for the velocity and for the pressure independently. In [6] it has been shown that a segregated procedure is more effective compared with a monolithic approach, where a unique POD basis (for both velocity and pressure) is built. Moreover, we resort to a vector-valued POD basis for the velocity, in contrast to what has been done in Section 2.2, where the same (scalar) modal basis is adopted for each component of the velocity. This choice is consistent with standard reduced order modeling techniques in a finite element framework [4,40]. Thus, we assemble two distinct response matrices, U u ∈ R dmuN h,u ×M for the velocity and U p ∈ R mpN h,p ×M for the pressure. Then, by mimicking the scalar case, we compute the correlation matrices associated with U u and U p , and we retain the first N u and N p eigenvectors for the velocity and for the pressure, respectively. This leads us to identify the reduced order spaces V mu,Nu and Q mp,Np , together with the corresponding matrices Υ mu,Nu and Π mp,Np collecting, by column, the POD basis functions for velocity and pressure, respectively. The basic HiPOD procedure is here modified to take into account the stability issue characterizing the approximation provided by a projection of the Stokes equations. Actually, it turns out that, even though the solutions involved in the offline phase are inf-sup stable, this does not guarantee a priori the inf-sup condition to the reduced space, with the possible generation of spurious pressure modes. Following [4,40], to recover the inf-sup property for the POD approximation, we enrich the velocity space V mu,Nu with the so-called supremizer solutions. In particular, to preserve the offline/online paradigm, we properly modify the procedure proposed in [4]. Let p mp (µ (i) ) be the i-th column of matrix U p , for i = 1, . . . , M . We solve the additional HiMod systems (20) X mu,u s mu (µ (i) ) = B T mp,mu (µ (i) )p mp (µ (i) ), for = 1, . . . , M , thus obtaining the M supremizer solutions s mu (µ (i) ) ∈ R dmuN h,u . Here X mu,u ∈ R dmuN h,u ×dmuN h,u denotes the HiMod matrix associated with the inner product in V mu (so that we use the same modal basis for both velocity and supremizers), while B mp,mu (µ) encodes the HiMod discretization of the bilinear form b(v, p(µ); µ). Successively, we assemble the response matrix U s ∈ R dmuN h,u ×M associated with the supremizers together with the corresponding correlation matrix, and we build the matrix Ξ mu,Ns collecting the first N s POD supremizer basis functions, with N s < M . Finally, we define the matrix Φ mu,Nu+Ns = [Υ mu,Nu , Ξ mu,Ns ] ∈ R dmuN h,u ×(Nu+Ns) and the enriched velocity space V mu,Nu+Ns spanned by the columns of Φ mu,Nu+Ns . Throughout the paper, we will refer to V mu,Nu+Ns simply as to the reduced velocity space. Furthermore, we will always assume N u = N s = N p = N .

3.2.2.
The online phase. We extend here the procedure introduced in Section 3.1.2. For any µ ∈ D, with µ = µ (i) and i = 1, . . . , M , rather than solving the corresponding HiMod system (21) A mu (µ) B T mp,mu (µ) B mp,mu (µ) 0 we rely on the reduced system where u mu,2N (µ) ∈ R 2N and p mp,N (µ) ∈ R N denote the POD reduced approximations for the velocity and the pressure, 0 is the null vector in R N , and where we assume that matrices A mu,N (µ) = Φ T mu,2N A mu (µ)Φ mu,2N ∈ R 2N ×2N , B mp,mu,N (µ) = Π T mp,N B mp,mu (µ)Φ mu,2N ∈ R N ×2N and the vector f mu,N (µ) = Φ T mu,2N f mu (µ) ∈ R 2N can be efficiently assembled owing to affine parameter dependence. Finally, POD solutions u mu,2N (µ) and p mp,N (µ) are projected back to the HiMod space, to yield the approximations for the HiMod velocity and pressure in (21). The HiPOD approximations for the HiMod solutions u mu (µ) and p mp (µ) will be denoted in the sequel by u mu,2N,POD and p mp,N,POD , respectively.
Concerning the stability of the POD reduced problem, for any µ ∈ D, one can numerically compare the inf-sup constant associated with the HiMod discretization, (23) β mu,mp (µ) = inf qm p ∈Qm p , with the corresponding constant resulting from the POD reduction procedure, where, with an abuse of notation, we have adopted the same symbol for the continuous HiPOD spaces as for the corresponding discrete counterparts. Practical computations for these constants rely on generalized eigenvalue problems (see, e.g., [11]). In particular, we resort to the formulas mu,mp,N , denote the minimum eigenvalue of the generalized problems respectively, where X mu,u is defined as in (20), X mp,p ∈ R mpN h,p ×mpN h,p denotes the HiMod matrix associated with the inner product in Q mp , while X mu,u,2N ∈ R 2N ×2N , X mp,p,N ∈ R N ×N are the corresponding reduced order matrices, given by X mu,u,2N = Φ T mu,2N X mu,u Φ mu,2N , X mp,p,N = Π T mp,N X mp,p Π mp,N , and where, to simplify the notation, we have removed the subscript POD to the HiPOD approximations. We consider a sampling set S consisting of 100 randomly selected values, {µ (1) , µ (2) , . . . , µ (100) }. Then, we hierarchically reduce problem (10) for each parameter in S, by preserving the same HiMod discretization as the one adopted in Figures 3-5, right. Figure 9 shows the trend of the eigenvalues of the correlation matrix associated with the HiMod velocity, pressure and supremizers. The drop to the numerical precision occurring at the fourth eigenvalue in all the plots suggests us to set N = 4. Indeed, due to the superposition property and since we deal with a linear problem, four independent basis functions are enough to span the space of the solutions to this parametrized Stokes problem. In the online phase we yield an approximation for the HiMod discretization corresponding to the parameter µ = [5, 10, 0, 3, 0] T , i.e., for the solution provided in Figure 3-5, right. Figures 10-12, left show the contour plot of the HiPOD approximation for the two components of the velocity and for the pressure. Comparing the three plots with the corresponding ones in Figures 3-5, right we recognize that the selected POD basis suffices to provide a reliable approximation, at least qualitatively. Nevertheless, if we investigate more into the details the distribution of the HiPOD error, we remark that while the two components of the HiMod velocity are approximated up to the machine precision, this is not the case for the pressure as shown by the contour plots in Figure 12, right. This lack of accuracy for the pressure is consistent with what already observed in [6], and deserves a further investigation, for instance, by considering different choices for the supremizers [38,2].
Finally, in Figure 13 we compare the trend of the HiPOD inf-sup constant β mu,mp,N,POD (µ) with β mu,mp (µ) in (23). For this purpose, we set N u = N p = 4 and we make N s varying between 0 and 4. To simplify notation, we preserve the same notation as in (24) although hypothesis N u = N p = N s is here removed. It is evident that the system is unstable for N s = 0 so that    supremizers are required to recover a reliable pressure, while β mu,mp,N,POD (µ) reaches a value comparable with β mu,mp (µ) when N s = 4.

The HiRB approach
As an alternative to the HiPOD reduction, we introduce a new technique to deal with a parametrized setting. We aim at applying a Reduced Basis (RB) approach to HiMod solutions, by relying on a greedy algorithm during the offline phase (see, e.g., [16,22]). The combination of HiMod with RB justifies the name HiRB adopted to denote the new procedure.

HiRB reduction for ADR problems.
Analogously to what done in Section 3, we first exemplify the HiRB reduction on problem (12). In particular, we focus on the offline stage since HiPOD and HiRB essentially resort to the same projection procedure during the online phase. The key idea of the RB method is to generate the reduced space,Ṽ m,N , by a greedy algorithm, i.e., by adding a single function at a time to the reduced basis [21]. LetṼ m,k denote the reduced space known at the k-th iteration. Additionally, we assume to have an error estimator, η m,k (µ), for the modeling error associated with the reduced solution u m,k,RB (µ) ∈Ṽ m,k , such that with u m (µ) the high fidelity HiMod discretization. Algorithm 1 itemizes the steps constituting the offline phase of HiRB method. First, the greedy algorithm identifies as new parameter the value which corresponds to the most informative HiMod solution not yet included in the reduced space, since it maximizes the discrepancy between the HiMod space, V m , and the reduced space,Ṽ m,k (step 1.). The search performed by the greedy algorithm starts from a random choice, µ g , for the parameter and goes on until when a value k ∈ N is found such that max µ∈S η m,k (µ) < η, or the reduced space dimension reaches N , with η a user-defined threshold. Then, the HiMod approximation, u m (µ (k+1) g ) ∈ R mN h , is computed (step 2.) and orthonormalized with respect to the functions already included in the RB basis, collected in the matrix Φ m,k = [φ 1 , . . . ,φ k ] ∈ R mN h ×k (step 3.). This is justified by the fact that solutions at step 2. can be linearly dependent so that they cannot constitute a basis. We denote the new basis function yielded at step 3. byφ k+1 ∈ R mN h . Finally, the RB matrix is extended to include the new basis function, so that we havẽ ). This allows us to define the spaceṼ m,k+1 of the HiRB approximations associated with Φ m,k+1 .
Throughout the paper, the HiRB reduced space eventually yielded by Algorithm 1 is denoted byṼ m,N , with N possibly equal to k, for k < N , if the check at step 1. succeeds at the k-th iteration.
Finally, the HiRB online phase follows, by mimicking exactly what performed in Section 3.1.2, with matrix Φ m,N replaced byΦ m,N . In particular, we denote the HiRB approximation associated with the vector u m,N,RB (µ) :=Φ m,Nũm,N (µ) by u m,N,RB (µ), withũ m,N (µ) the solution of the reduced system corresponding to (17).
The choice of the error estimator, η m,k (µ), represents a key issue of the RB approach, in particular to ensure the convergence of the greedy algorithm as well as the reliability of the reduced order model. In general, η m,k (µ) demands the computation of the reduced solutioñ Φ m,kũm,k (µ), so that, at each iteration of Algorithm 1, an online phase of dimension k has to be carried out. For additional details in a standard RB setting, we refer the interested reader, e.g., to [10,13,16,21,38,39]. The most common choice for the modeling error estimator relies on the ratio between the dual norm of the weak residual associated with the reduced model and a lower bound for the coercivity constant of the reduced problem [21]. This is the choice pursued to define η m,k (µ) in (25), where we set the coercivity constant to one since it is not trivial to compute exactly this value. The reliability of this estimator is numerically verified in Section 5. The extension of the successive constraint method proposed in [23] to the HiMod setting represents a topic for a possible future investigation. 4.1.2. Numerical assessment. We adopt exactly the same setting as in Section 3.1.3, so that the parameter µ coincides with the vector [ν, b x , b y , σ] T and varies over the ranges, D 1 and D 2 , in (18). Algorithm 1 is run over a training set S consisting of 100 samples. Nevertheless, we omit to set the threshold η at step 1. and we fix a priori the dimension N of the reduced space to 20, also with a view to the comparison performed in Section 5. This leads to hierarchically reduce only 20 ADR problems, in contrast to 100 ADR problems with the HiPOD procedure. Finally, we run the HiRB online phase for µ = [5,20,75,25] T to approximate the HiMod solution to problem (2). Figure 14 shows the contour plot of the HiRB approximation for the two ranges of the parameter µ. The qualitative matching between these solutions and the HiMod approximation in Figure 2, right is good. By analyzing the distribution of the modeling error, u 8 (µ) − u 8,20,RB (µ), in Figure 15, it is confirmed that the smaller the parameter range the higher the accuracy of the HiRB approximation. In particular, the maximum error reduces more than one order when sampling µ in D 2 . A cross-comparison with the HiPOD approximations in Figures 7-8 highlights a slightly higher reliability for the HiRB approach for this test case. A more thorough investigation in such a direction will be performed in Section 5, together with an error analysis over a random testing.

4.2.
HiRB reduction for the Stokes equations. We detail the offline step of the HiRB reduction procedure on problem (19), while referring to Section 3.2.2 for the online phase. 4.2.1. The offline phase. In this context, we assume to have an error estimator both for the velocity and the pressure [18,38], such that with u mu (µ), p mp (µ) the high fidelity HiMod solution pair and u mu,2k,RB (µ), p mp,k,RB (µ) the HiRB approximation belonging to the RB spaceṼ mu,2k ×Q mp,k .
Algorithm 2 details the operations characterizing the HiRB offline phase when applied to the Stokes problem. Two are the main differences with respect to Algorithm 1, namely i) we pursue a segregated approach to build the reduced spaces for the velocity and the pressure, ii) the space for the velocity is enriched via the supremizers. After the greedy selection on the training set S (step 1.), we solve both the HiMod problem (21) and the HiMod supremizer equation (20) by setting µ = µ (k+1) g and µ (i) = µ (k+1) g , respectively thus obtaining the HiMod velocity and pressure pair, u mu (µ . compute the HiMod supremizer s mu (µ (k+1) g ); 3a. compute the new element,υ k+1 , of the RB basis for the velocity by orthonormalizing u mu (µ (k+1) g ) with respect toΥ mu,k ; 3b. compute the new element,π k+1 , of the RB basis for the pressure by orthonormalizing p mp (µ (k+1) g ) with respect toΠ mp,k ; 3c. compute the new element,ξ k+1 , of the RB supremizer basis for the velocity by orthonormalizing s mu (µ (k+1) g ) with respect toΞ mu,k ; 4a. build the updated RB matrixΥ mu,k+1 by includingυ k+1 as the (k + 1)-th column; 4b. build the updated RB matrix for the pressureΠ mp,k+1 by includingπ k+1 as the (k+1)-th column; 4c. build the updated RB matrixΞ mu,k+1 by includingξ k+1 as the (k + 1)-th column; 4d. build the updated RB matrix for the velocityΦ mu,2(k+1) , and then go back to 1.

4.2.2.
Numerical assessment. We adopt the same numerical setting as in Section 3.2.3 with the goal of approximating the HiMod solution in Figures 3-5, right with a reduced basis approach. Analogously to Section 4.1.2, we waive the opportunity to employ the threshold η in Algorithm 2, and we fix the dimension of the reduced spaces to N = 4 to match the choice in Section 3.2.3. Figures 16-18, left show the contour plot of the HiRB approximation for the two components of the velocity and for the pressure. The qualitative agreement both with the HiMod solution and with the HiPOD approximation in Figures 3-5, right and 10-12, left respectively confirms the reliability of the proposed procedure. The distribution of the HiRB modeling error in Ω is provided in Figures 16-18, right. The pressure is the quantity characterized by the worst accuracy, analogously to what obtained with the HiPOD approach. Nevertheless, we remark that HiRB technique furnishes an approximation of lower quality also for the y-component of the velocity when compared with the HiPOD approximation (notice the difference in terms of order of magnitude for the corresponding modeling errors in Figure 17, right and Figure 11, right, respectively). Finally, we recognize a more uniform distribution of the modeling error in Figure 16, right with respect to the corresponding trend of Figure 10, right. In the former case, the error is spread across the whole domain, whereas in the latter the error is mostly confined to the outflow boundary.

HiPOD versus HiRB
This section is meant to compare HiPOD and HiRB techniques. The comparison is carried out in terms of three main issues, which are investigated hereafter, separately.
For this purpose, we remind that both HiPOD and HiRB approaches actually carry out a twofold reduction. The first one is obtained with the HiMod discretization [3,20,27], while the second reduction is performed via a projection step during the online phase. The final expectation is the capability to have a reliable approximation for the (full) problem at hand, with a very contained computational effort. Figure 19. ADR test case, parameter range D 1 : H 1 (Ω)-norm of the modeling error associated with the HiPOD (left) and with the HiRB (right) reduction as a function of N . 5.1. Accuracy of the reduced problems. To compare HiPOD and HiRB in terms of accuracy, we plot the average of the associated error over a testing set of 100 randomly selected parameters for both the test cases in Sections 2.1 and 2.2.
In particular, in Figures 19 and 20 we show the trend of the H 1 (Ω)-norm of the modeling error characterizing the ADR test case and for the two choices of the parameter range in (18). For HiRB approximations, we provide also the trend of the error estimator. HiPOD exhibits a slightly higher accuracy with respect to HiRB (about half an order of magnitude), for both D 1 and D 2 . This is likely related to the adopted estimator which under-estimates the exact error of about one order of magnitude (see Figures 19 and 20, right). Actually, η m,k (µ) is an error indicator rather than an error estimator, since we have set the coercivity constant to one. This might result in a sub-optimal greedy selection, although the (monotonic) decreasing trend of the exact error is correctly captured by η m,k (µ). Finally, as expected, the reduced order approximation associated with D 2 is more accurate for both the procedures.
The discrepancy between HiPOD and HiRB in terms of accuracy is less evident when considering the Stokes test case. For N = 4, the velocity is approximated almost at the machine precision, while the pressure is characterized by a modeling error of the order of 10 −2 with respect to the L 2 (Ω)-norm. The lower accuracy of the pressure is consistent with what remarked in Figures 12 and 18. Possible improvements in such a direction are suggested in Section 3.2.3. Moreover, the computation of separate error bounds for the velocity and the pressure would be helpful in improving the accuracy, despite requiring further evaluation of stability factors [18].

5.2.
Speedup of the reduced problems. We investigate here the computational effort demanded by the online stage of the HiPOD and HiRB methods. We quantify such an effort in terms of CPU time 1 . In particular, we quantify the speedup characterizing the two approaches with the ratio τ m (µ)/τ m,N (µ), where we denote the elapsed time associated with the standard HiMod approximation by τ m (µ), and the time required to solve the corresponding HiPOD or HiRB system by τ m,N (µ). A value of speedup greater than one results in a computational gain. Table 1 gathers the values of this investigation. In order to filter out any dependence on µ, we compute the speedup index over a testing set of 100 randomly selected parameters. We observe a large speedup for the Stokes test case and a very mild sensitivity with respect to N for the  ADR problem, independently of the adopted technique. More in detail, we point out a general lower speedup (of about one third) for the HiRB procedure when compared with HiPOD, in particular for the ADR test case. This is due to the fact that quantity τ m,N (µ) includes also the time elapsed for the evaluation of the error estimator in the HiRB case, whereas this is not required by the HiPOD procedure. Actually, the HiPOD and the HiRB speedups become very similar for the Stokes test case.

5.3.
Cost of the offline phase. We focus now on the offline phase, by comparing the total time required by the HiPOD and HiRB procedures to build the reduced basis.
It is reasonable that, for a fixed dimension, N , of the reduced space, the HiRB reduction requires less offline time than HiPOD. Actually, to extract the reduced basis, the HiPOD approach computes the HiMod discretization for each of the M parameters in the sampling set, S, and, only a posteriori, compresses such information into a reduced basis of dimension N , with N < M . On the contrary, the HiRB method iteratively generates the reduced basis by adding  a new basis function at each iteration of the greedy algorithm. Thus, we compute exactly N HiMod solutions, out of the M possible approximations, associated with the parameters in S. Nevertheless, the offline stage of the HiRB method includes the evaluation of the error estimator during the greedy selection. A key requirement is that this evaluation is computationally cheap. However, the construction of the data structures (e.g., higher order tensors [21]) required for this purpose usually entails an additional computational cost, which might dominate the overall offline cost if M is small. Finally, further less relevant differences between the two methods can be pointed out, such as the CPU time required to solve the eigenvalue problem associated with the HiPOD reduction. Figure 23 compares the trend of the total CPU time demanded by the offline stages of the HiPOD and HiRB procedures, as a function of the size M of the sampling set, when applied to the test case in Section 2.1 and for N set to 20. In agreement with what remarked above, it follows that the HiRB training is more expensive than the HiPOD one for small values of M . For instance, for M = 50, HiPOD takes 25s, whereas HiRB requires more than 90s, most of the time being spent in the setup of the error estimator. Vice versa, for large values of M , HiRB demands less time than HiPOD. For example, when M = 300, HiPOD is more time-consuming than HiRB, by requiring 140s to be compared with 110s.
Finally, we remark the different slope characterizing the two plots in Figure 23. The mild slope of the HiRB curve confirms that the computational effort required by the evaluation of the error estimator is essentially independent of the size M . On the contrary, the considerable slope of the HiPOD curve highlights that the computation of the HiMod approximations is not negligible and, in general, heavier than the evaluation of η m,k (µ).

Conclusions
This work has to be meant as a first attempt to compare the new reduction techniques HiPOD and HiRB, for the modeling of parametrized problems. HiPOD has been introduced in [26,6], whereas the HiRB approach is proposed here for the first time. The two methods are then compared on a benchmark ADR and Stokes problem. Starting from this comparison, we can state that HiRB is better performing than HiPOD for large training sets, thus turning out to be the ideal tool to tackle, for instance, time demanding fluid dynamics problems. As expected, the greedy algorithm allows us to reduce the offline time. The weak point of the HiRB approach remains the availability of a reliable error estimator. So far, to simplify the introduction of the new method, we have adopted an error indicator coinciding with the residual associated with the reduced solution, by completely neglecting the coercivity constant. This rough choice actually leads to under-estimate the exact error, with a consequent performance loss in terms of speedup and slightly of accuracy with respect to the HiPOD procedure. However, these conclusions have to be considered as preliminary since we have limited our analysis only to two test cases and, clearly, a more thorough investigation is deserved.
An important issue related to both HiPOD and HiRB has concerned the inf-sup stability which is not necessarily guaranteed for the reduced formulations. To tackle this matter in both the cases, supremizer enrichment has been employed with significative improvements. The pressure approximation for the Stokes equations still demands some amendment for both the techniques. The proposal of different supremizers likely represents a viable remedy in such a direction. Among other future developments of possible interest, we cite the generalization to threedimensional and to nonlinear problems, as well as to an unsteady framework. Finally, to certify the reliability of the two methods, a more rigorous investigation of the accuracy characterizing HiPOD and HiRB procedures is desirable, by properly combining HiMod estimates in [31,3] with the well-established accuracy results on POD and RB [41,21].