{'fraction_non_alphanumeric': 0.037698320202632356, 'fraction_numerical': 0.016393952107904838, 'mean_word_length': 4.43622233485386, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, '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}}
Computational reduction strategies for the detection of steady bifurcations in incompressible fluid-dynamics: applications to Coanda effect in cardiology

We focus on reducing the computational costs associated with the hydrodynamic stability of solutions of the incompressible Navier-Stokes equations for a Newtonian and viscous fluid in contraction-expansion channels. In particular, we are interested in studying steady bifurcations, occurring when non-unique stable solutions appear as physical and/or geometric control parameters are varied. The formulation of the stability problem requires solving an eigenvalue problem for a partial differential operator. An alternative to this approach is the direct simulation of the flow to characterize the asymptotic behavior of the solution. Both approaches can be extremely expensive in terms of computational time. We propose to apply Reduced Order Modeling (ROM) techniques to reduce the demanding computational costs associated with the detection of a type of steady bifurcations in fluid dynamics. The application that motivated the present study is the onset of asymmetries (i.e., symmetry breaking bifurcation) in blood flow through a regurgitant mitral valve, depending on the Reynolds number and the regurgitant mitral valve orifice shape.


Introduction
We focus on the hydrodynamic stability of solutions of the incompressible Navier-Stokes equations for a Newtonian and viscous fluid in contraction-expansion channels, with a particular concern on steady bifurcations. Steady bifurcations occur when new, non-unique solution branches of the Navier-Stokes equations appear as physical and/or geometric control parameters are varied. When the fluid domain is characterized by two or three dimensions with non-periodic boundary conditions, the formulation of the stability problem requires solving an eigenvalue problem for a partial differential operator. See [23] for a review on numerical methods for stability analysis based on linearized eigenvalue problems. An alternative to the eigenvalue problem approach is the direct simulation of the flow to characterize the asymptotic behavior of the solution; see, e.g., [33,3,54]. Both approaches can be extremely expensive in terms of computational time. In this paper, we propose to apply Reduced Order Modeling (ROM) techniques to reduce the demanding computational costs associated with flow stability analysis.
Practical applications of contraction-expansion channel flows include equipments such as heat exchangers, combustion chambers, and mixing vessel. An application that motivated the present study is the onset of asymmetries (i.e., symmetry breaking bifurcation) in blood flow through  a regurgitant mitral valve, depending on the Reynolds number and the regurgitant mitral valve orifice shape. Mitral regurgitation is a valvular disease characterized by abnormal leaking of blood through the mitral valve from the left ventricle into the left atrium of the heart. See Figure 1.
In certain cases the regurgitant jet "hugs" the wall of the heart's atrium as shown in Figure 1 (c). These eccentric, wall-hugging, non-symmetric regurgitant jets have been observed at low Reynolds numbers [67,1] and are said to undergo the Coanda effect [65,70]. This effect, described as the tendency of a fluid jet to be attracted to a nearby surface, owes its name to Romanian aerodynamics pioneer Henri Coanda. The primary tool to assess the severity of mitral regurgitation is echocardiography [73]. One of the biggest challenges in echocardiographic assessment of mitral regurgitation is the Coanda effect: the wall-hugging jets appear smaller in the color Doppler image of regurgitant flow, leading to a gross under-estimation of regurgitant volume by inexperienced observers [31,16]. As a result, patients requiring treatment may not be recognized. Despite the large cardiovascular and bioengineering literature reporting on the Coanda effect in echocardiographic assessment of mitral regurgitation, there is very little connection with the fluid dynamics literature that could help identify and understand the main features of the corresponding flow conditions. In this paper, our goal is to understand what triggers the Coanda effect in a simplified setting. A contraction-expansion channel is a simplified setting which has the same geometric features of mitral regurgitation. In fact, a mitral regurgitant jet flows from the left ventricle through the contraction between the mitral leaflet, called regurgitant orifice, into the left atrium. First, we focus on planar contraction-expansion channels (see Fig. 2) and investigate the influence of the Reynolds number and the contraction width w c (i.e., the orifice height) on the flow. Then, we consider the 3D geometry reported in Fig. 3 to understand the role played by the channel depth h (i.e., the orifice length). Eccentric regurgitant jets typically occur in prolapsed mitral valves, i.e. when two valve flaps of the mitral valve do not close evenly. Thus, another parameter of interest, although not considered in this work, could be the orifice depth. Moreover, for a more realistic setting one would have to account for the pulsatility of the flow and include the Strouhal number among the parameters.
We remark that the focus of this paper is to investigate the cause of the Coanda effect in simplified settings. Nonetheless, it is thanks to the results reported here that our medical collaborators at the Houston Methodist DeBakey Heart & Vascular Center were able to reproduce the Coanda effect in a mock heart chamber (see Sec. 4.4). A comparison between the experiments in vitro and corresponding 3D simulations is presented in [69].
The incompressible fluid dynamics in a planar contraction-expansion channel has been widely studied from both theoretical and practical perspectives; see, e.g., [49,24,63,26,36,48] and references therein. In the two-dimensional geometry reported in Fig. 2, the wall-hugging effect happens only above a critical Reynolds number (12), which depends on the expansion ratio λ defined in (11). Compare Fig. 5(e) and Fig. 5 (b), which correspond to a Reynolds number above and below the critical value, respectively. The asymmetric, wall-hugging solution remains stable for a certain range of Reynolds number and asymmetries become stronger with the increasing Reynolds number, as shown in [48]. The formation of stable asymmetric vortices in 2D planar expansion is attributed to an increase in velocity near one wall that leads to a decrease in pressure near that wall [70]. Once a pressure difference is established across the channel, it will maintain the asymmetry of the flow. The critical value of the Reynolds number has been identified for different expansion ratios λ. In particular, it was found that such critical value decreases with increasing value of λ (see [24,57]).
In the three-dimensional geometry reported in Fig. 3, the critical Reynolds number for the symmetry-breaking (i.e., the wall-hugging) varies with the expansion ratio and the aspect ratio defined in (13), as shown in [17,18,52]. When the expansion ratio is fixed and the aspect ratio decreases, the endwall influence becomes more important: the critical Reynolds number increases [18,52]. For moderate aspect ratios, the flow is steady in time but highly three-dimensional, and complex spiraling structures are observed, which are not closed recirculating cells as in the case of 2D flows. See Fig. 15. The numerical studies in [66] found that the flow only resembles a 2D flow for very large aspect ratios. The theoretical study of Lauga et al. [41] shows that for low aspect ratios the flow is highly three-dimensional. The numerical and experimental studies in [52] show that the strong three-dimensional effects appearing for low aspect ratios inhibit the wall-hugging effect observed in geometries with high aspect ratios at the same Reynolds number. This suggests that the eccentric regurgitant jets, such as the one in Figure 1 (c), occur when the regurgitant orifice is long (large aspect ratio) and narrow (large expansion ratio).
Given the relatively fast decay of energy spectrum for flows at sufficiently low Reynolds numbers, a ROM technique is expected to be an efficient tool to reduce the prohibitive computational costs associated to identifying the flow conditions and geometries that trigger asymmetries. Recent developments of ROM techniques have focused on the reduction of computational time for a wide range of differential problems [19,55], while maintaining a prescribed tolerance on error bounds [59,38,56]. Terragni and Vega [64] showed that a Proper Orthogonal Decomposition (POD) approach allows for considerable computational time savings for the analysis of bifurcations in some nonlinear dissipative systems. Herrero, Maday and Pla [37] have used a Reduced Basis (RB) method to speed up the computations of different solution branches of a two-dimensional natural convection problem (Rayleigh-Bénard), achieving a good accuracy but without investigating the approximation of bifurcation points. For each fixed aspect ratio, multiple steady solutions for the Rayleigh-Bénard problem can be found for different Rayleigh numbers and stable solutions coexist at the same values of external physical parameters. In [37], it is shown that stable and unstable solutions are correctly identified by the RB method. Yano and Patera [72] introduced a RB method for the stability of flows under perturbations in the forcing term or in the boundary conditions, which is based on a space-time framework that allows for particularly sharp error estimates. Furthermore, in [72] it is shown how a space-time inf-sup constant approaches zero as the computed solutions get close to a bifurcating value.
In a previous work [53] we have investigated steady and Hopf bifurcations in a natural convection problem dealing with a geometrical parameter (the cavity length), and a physical parameter (the Grashof number). This work is an extension of that study and provides a proof of concept of the applicability of reduced order methods to investigate stability and bifurcations in complex fluid dynamic problems at a reasonable computational cost. The proposed framework allows the use of a black-box input-output toolbox to be managed also by non-expert scientists in computational sciences. The offline-online splitting of the computational procedure is crucial in view of the use of a High Performance Computing (HPC) infrastructure for the offline computational step (expensive and time consuming) and a light modern device, such as tablet or smart phone, for online calculations. The idea is to use different platforms (and methodologies) for a strategic computational collaboration between high order and reduced order methods, with competitive computational costs for complex simulations. This could still be considered a research frontier in computational fluid dynamics, especially in view of real life applications.
The outline of the paper is as follows. The general problem setting is described in Section 2, while the numerical methods are treated in Section 3. Numerical results in 2D and 3D contractionexpansion channels under different parametrizations are reported in Section 4. Conclusions and perspectives follow in Section 5.

Problem setting
Let Ω ⊂ R d , d = 2, 3, be the computational domain. The motion of an incompressible, viscous fluid in a spatial domain Ω over a time interval of interest (0, T ) is governed by the incompressible Navier-Stokes equations where u is the velocity and ß is the Cauchy stress tensor. In large arteries and inside the heart, it is widely accepted to model blood as a Newtonian fluid. See, e.g., [29] and references therein. For such fluids, σ(u, p) = −pI + 2ν (u), where p is the pressure, ν is the fluid kinematic viscosity, and (u) = (∇u + (∇u) T )/2 is the strain rate tensor. Eq. (1) represents the conservation of the linear momentum, while eq. (2) represents the conservation of the mass. Eq. (1)-(2) need to be endowed with boundary and initial conditions, e.g.: where Γ D ∪ Γ N = ∂Ω and Γ D ∩ Γ N = ∅. Here, d, g, and u 0 are given. We assume that there is no external body force acting on the fluid (see eq. (1)) and the motion is driven by the boundary conditions (3)-(4).
When the fluid acceleration is negligible (i.e., the system has evolved towards a steady state), eq. (1)-(2) can be replaced by: To characterize the flow regime under consideration, we define the Reynolds number as where U and L are characteristic macroscopic velocity and length respectively. We will characterize U and L for the specific cases under consideration in Sec. 2.1 and 2.2.
x y L c w c Figure 2: Scheme of the planar contraction-expansion channel considered in this work. The notation is the same used in [52].
If eq. (5)-(6) are used to model the fluid dynamics, then the first term in eq. (8) is disregarded.
The nonlinearity in problem (1)-(2) can produce a loss of uniqueness for the solution, with multiple solutions branching from a known solution at a bifurcation point. As will be explained in Sec. 2.1, the Coanda effect is associated with a steady bifurcation point. To detect numerically the presence of a steady bifurcation point, we will rely on the spectrum analysis of a linearized operator. See, e.g., [2] for a theoretical introduction to bifurcation theory, and [20] and [23] for applications to numerical analysis. Following [20, par. 7.3], we introduce the linearization L : V D × V → V of the convection operator in eq. (5), obtained by Fréchet differentiation about a base point u * of the term u · ∇u: At a symmetry breaking bifurcation point, a simple eigenvalue of L changes sign. We remark that for the kind of bifurcations we are interested in (namely supercritical pitchfork bifurcations) it is sufficient to study the spectrum of the antisymmetric part of the linearized operator [20]. Other techniques would have to be used for different kind of singular points (e.g., Hopf bifurcations, fold points) occurring at higher Reynolds numbers or in different settings.

2D case
The domain under consideration for the two dimensional case is shown in figure 2. The following boundary conditions are imposed in this case: homogeneous Dirichlet (no-slip) boundary condition on the sides drawn with a continuous line, non-homogeneous Dirichlet boundary condition (timeindependent parabolic velocity profile) on the red dashed line which corresponds to the inlet, and homogeneous Neumann (stress-free) boundary condition on the blue dashed line which corresponds to the outlet. For space limitation, the channel depicted in Fig. 2 is shorter than the actual one. The actual domain length past the expansion is 6 times the channel height L c . The length of the contraction channel is equal to (L c − w c )/2, and the distance between the inflow section and the contraction is equal to L c .
For the characterization of the flow in the 2D case, we introduce the following quantities: expansion ratio: average horizontal velocity: v x = Q w c , where Q is the flow rate; -Reynolds number: Notice that the above definition of Re 2D does not coincide with taking L = w c and U = v x in (7). The reason for the extra factor 2 will be explained in the next subsection. In the numerical simulations, we vary Re 2D by changing the value of the viscosity ν. In the geometry considered by [52] the expansion ratio λ is equal to 15.4. In order to reproduce the results in [52] for validation purposes, in Sec -Symmetric jet: for slightly larger Reynolds numbers there is a breaking of the vertical symmetry. The flow is still symmetric with respect to the horizontal axis, but the two vortices downstream of the expansion are larger than the vortices upstream. See Fig. 5 (b). At a further increase of the Reynolds number, the vertical asymmetry of the flow becomes increasingly evident, yet the horizontal symmetry is maintained.
-Asymmetric jets: when the Reynolds number is sufficiently large, the configuration with a symmetric jet is still possible, but unstable [63]. In fact, small perturbations 1 expand one recirculation zone and shrink the other, causing a drastic variation in the flow. See Fig. 5 (c) and (d). This horizontally asymmetric solution remains stable for a certain range of Re 2D and asymmetries become stronger with the increasing Reynolds number, as shown in [48]. At the minimum value of the Reynolds number for which the asymmetric jet configuration exists there is a symmetry breaking bifurcation point or steady-state bifurcation point. In Sec. 4.3 we will show how the critical value of the Reynolds number for the symmetry breaking bifurcation changes as λ varies.

3D case
A 3D channel is obtained by extruding the 2D geometry considered in the previous section in the direction orthogonal to the flow plane. The goal of this test case is to study the influence of the channel depth on the flow pattern, and in particular on the symmetry breaking bifurcation point. Once λ is fixed, the 3D problem depends on two parameters: the Reynolds number and the channel depth h. See Figure 3. Note that also in this case the geometry reported in the figure has been cropped due to space limitation. In reality, we considered a channel length past the expansion equal to 6 times the channel height L c .
The boundary conditions are the same as in the 2D case: we impose a parabolic velocity profile at the inlet, a stress-free condition at the outlet, and a no-slip condition everywhere else.
We introduce the following quantities, which are useful in the characterization of the numerical simulation [52]: aspect ratio: normalized channel depth: Note that H = 1 is the limit case of infinite channel depth, which corresponds to the 2D configuration. The definition of Re 3D has been obtained by setting U = v x and L equal to the hydraulic diameter of the contraction channel in (7). We remark that the 2D case can be seen as a limit of the 3D case for AR → ∞. This justifies the factor 2 in the definition of Re 2D in (12). Another reason to define Re 2D in (12)

Numerical method
We are interested in adopting a Reduced Order Model (ROM) for the numerical solution of the problems presented in the previous section. Reduced Order Models have been introduced for parametrized problems requiring real-time capabilities due to a many-query setting. The goal is to compute reliable results at a fraction of the cost of a conventional (full order) method. A practical way to realize this is to organize the computation in two steps: -An offline phase: full order approximation solutions corresponding to selected representative parameters values/system configurations are computed and stored, together with other information concerning the parametrized problem. This is a computationally expensive step usually performed on high performance computing facilities.
-An online phase: the information obtained during the offline phase is used to compute the solution for a newly specified value of the parameters in a short amount of time (ideally in real time), even on a relatively low power device such as a laptop or a smartphone.
These split computational procedures are built in such a way that new parameter dependent quantities are easily and quickly computed online, while representative basis functions for selected parameter values and more demanding quantities are pre-computed offline. We refer to [40] for a review of ROM in Fluid Mechanics.
The problem under consideration might depend on several parameters, with each parameter varying in a certain range. We introduce a parameter vector µ that contains all parameters. If the problem depends on two parameters, we have µ = (µ 1 , We consider both physical parameters (e.g., Re) and geometric parameters (e.g., λ and H). To stress the solution dependence on the parameter(s), we will use the notation u = u(µ) and p = p(µ), without implying that there is a one-to-one correspondence between µ and u or p. The maps µ → u(µ) and µ → p(µ) are one-to-one only on the region of D where there exists a unique solution of problem (1)- (2).
The treatment of the geometric parametrization deserves further explanation. Let ξ ∈ D be a geometric parameter the problem depends on. We select a reference domain Ω that is mapped to the parametrized domain Ω (ξ) through a one-to-one, orientation preserving transformation T : D × Ω → Ω (ξ). Using this map, we can cast eq. (8)-(9) into the reference domain. For instance, eq. (9) becomes: where F(ξ) is the Jacobian matrix of transformation T (ξ) and J(ξ) its determinant. Among many Reduced Order Models available in the literature, we choose a Reduced Basis (RB) method. We will briefly recall the main features of RB methods in the following sections. For a general review on the RB method we refer to, e.g., [59,38,56].

Full order approximation
As full order approximation scheme for eq. (8)- (9) to be used in the offline phase, we choose the Spectral Element Method (SEM). See, e.g., [22,14,13] for a general review of SEM and application to fluid mechanics. We adopt the SEM implementation available in open source software Nek5000 [27], where the basis functions for each element are the Lagrange interpolants on a Gauss-Lobatto-Legendre tensor product grid. We refer to [22] for an introduction to efficient SEM implementation.
Let V N , Q N , and V N 0 be the Spectral Element spaces, which are finite dimensional approximations of the infinite dimensional spaces V D , Q, and V , respectively. The full order approximation problem reads: for a given µ ∈ D, Notice that in eq. (15)- (16) we have already accounted for the fact that g = 0 in (4) for both the 2D and 3D case. For the computations in Sec. 4, we choose the stable P 11 − P 9 couple for velocity and pressure approximation. In the Nek5000 solver, the aliasing errors associated with the choice of high order polynomials for approximating the nonlinearity are dealt with the 3/2 rule (also called zero-padding rule, see [10]). This rule consists in evaluating the integrals of the nonlinear term (to be liearized) using a quadrature formula with 3/2 times the quadrature points of the other terms, so that the aliasing errors contribute only for those wavelengths that are filtered out by the grid size. For the time discretization of eq. (15)-(16) we adopt a Backward Differentiation Formula of order 3 (BDF3; see, e.g., [22]). The convective term is treated explicitly, with a third order extrapolation formula as explained in Nek5000 documentation [50]. Such a treatment of the convective term does not guarantee the unconditional stability in time of the linearized numerical scheme. A CFL condition has to be verified at every collocation point.
Given an initial solution, we consider the system to be close enough to the steady state when the following stopping condition is satisfied: with tolerance tol = 10 −8 . When the stopping criterion (17) is met, the simulation is interrupted.

Sampling
The sampling process consist in selecting N parameters {µ i }, with i = 1, . . . , N , in the parameter space D, whose corresponding solutions {u N (µ i )} ⊂ V N and {p N (µ i )} ⊂ Q N will be used to construct the Reduced Basis spaces for velocity and pressure, respectively. Solutions u N (µ i ) and p N (µ i ), with i = 1, . . . , N , are called snapshots. In order to simplify the notation, we will denote u N (µ i ) by u(µ i ) and p N (µ i ) by p(µ i ). In this section, we are going to explain how to select the velocity snapshots. The same procedure can be applied to obtain the pressure snapshots. Let µ k be the k-th component of parameter vector µ and let D k be the interval of interest for such component. The sampling procedure described below will select N k values of µ k in D k and the total number of sample parameter vectors is N = k N k : 8 For each component µ k of the parameter vector µ we choose as µ j k the Chebyshev points: where µ k,min = min This procedure to sample sample points in D is called Gauss-Lobatto-Chebyshev (GLC) tensor product collocation strategy [71]. For example, Fig. 4   The GLC collocation points are a practical choice, since in this case many sampling points share the same geometric parameter, requiring to start the continuation method very few times [71]. Other sampling methods such as Greedy [59] or CVT [53] would require to start the continuation method for each new sampling point. A disadvantage of the GLC collocation strategy is that the approximation spaces are not hierarchical, meaning that the RB spaces for a certain value of N are not in general subspaces of the RB spaces obtained for a higher value of N . This could increase the offline computational cost in case we need to enrich the RB spaces.
Our conjecture is that for bifurcation problems it may be useful to cluster the sampling points close to the bifurcation points. However, to the best of our knowledge, there are no error estimates for different sampling methods for steady state Navier-Stokes equations involving bifurcation points. In any case we are enriching our investigation with an eigenvalue analysis to detect the bifurcation point at the reduced level as well, as we will introduce in Sec. 3.5.

Construction of the RB spaces
For every selected sample µ i , we solve the full order approximation problem (15)-(16) until stopping criterion (17) is satisfied to get u(µ i ) and p(µ i ). After the sampling is complete, we have two which generate the finite dimensional subspaces V N and Q N , called Reduced Basis spaces. The key feature of a correct ROM is that the dimension of the reduced order space is much lower that the dimension of the full order space: In this way, all the computations required by the online phase (see Sec. 3.4) will be much less expensive that the computations required during the offline phase.
In this section, we are going to focus on how to construct the Reduced Basis The same procedure can be applied to obtain the Reduced Basis To actually build the space V N , it is usually preferred not to evaluate directly the inner products in eq. (15)- (16) for all the snapshots {u(µ i )}, for two reasons: -the snapshots may contain some redundant information, due to the sampling procedure, thus leading to linear dependence and ill-conditioning during the matrix assembling.
-in general the snapshots will not be orthogonal to each other, and consequently they will generate a full mass matrix, increasing the storage requirement and operations count.
Thus, it is preferred to compute an orthonormal generating set {φ i } N i=1 for V N so that the resulting linearized problem is well conditioned. Two of the most popular techniques to compute othogonal basis functions are the Proper Orthogonal decomposition (POD) [68] and the Gram-Schmidt orthogonalization (GS) with its variants [60].
One of the ways to build a POD is to compute the correlation matrix for the set of snapshots, defined as as The eigenvalues λ i and eigenvectors ψ i of C are computed, and each POD basis vector is defined as where ψ i,k denotes the k-th component of the i-th eigenvector. The eigenvalue λ i associated to each POD mode is related to the fraction of energy stored in the corresponding mode. The POD modes are automatically orthogonal in the L 2 inner product, but not normal in general. It can be shown [68,35] that the space generated by the POD, denoted by V N POD , minimizes the projection error in the L 2 norm: where Π V N : V → V N is the projection operator on the space V N generated by the POD modes.
For the results reported in Sec. 4, we have used POD to compute the basis functions. However, we would like to remark that POD has one major drawback when the number of snapshots is large: the number of operations required to compute the correlation matrix (22) and its eigenpairs becomes prohibitive. If this is the case, a Gram-Schmidt orthogonalization is usually preferred over computing the basis functions by means of a POD [60]: where the normalization step is not always adopted, since it may lead to an ill-conditioned linear system. As mentioned above, the same orthonomalization method can be applied to obtain the pressure basis {σ i } N i=1 . The reduced basis spaces V N and Q N are defined as:

Online phase computation
After the construction of the Reduced Basis spaces, an online approximation of the solution can be computed by applying the Galerkin projection to the spaces V N and Q N . Namely, given a target parameter µ ∈ D we search for (u N (µ), For convenience, for the rest of this section we are going to assume that the first term in eq. (27) is negligible, as if the flow was modeled by eq. (5)- (6).
The convective term is linearized with a fixed point scheme. Suppose that an initial tentative solution u N 0 (µ) is known. Given u N k−1 (µ), at the k-th iteration of the fixed point method we solve problem: The iterative scheme can be stopped for instance when an increment-based residual: is below a given tolerance. The solution scheme described so far requires the V N − Q N pair to satisfy a stability condition called inf-sup condition or Ladyzhenskaya-Brezzi-Babuška (LBB) condition: See, e.g., [39,4,11,9,25]. Spaces V N and Q N in (26) computed using the POD or GS modes as explained in Sec. 3.3 are not guaranteed to fulfill condition (32). There are two options for circumventing this issue: casting the problem into a divergence-free space (see, e.g., [53]) and enforcing approximation stability properties for the V N − Q N pair (see, e.g., [60,58,6]). Here, we choose the former approach. This means that we require V N to be a subset of H 1 div (Ω ): which is a subspace of H 1 (Ω ). If the basis functions for V N are divergence-free, eq. (30) is no longer needed. Thus, the pressure disappears from the variational formulation and we do not need to build the space Q N . See, e.g., [28]. With a divergence-free basis set for V N , at every fixed-point iteration we have to solve the following linear system: where u k is the vector containing the projection coefficients of u N k−1 onto the space V N , b k ∈ R N depends from the specified boundary conditions, and A k (µ) ∈ R N ×N is given by: Once the velocity u N (µ) ∈ H 1 div (Ω ) has been computed, the pressure can be recovered, for example, by solving a Poisson problem online: We refer to,e.g., [12] for an analysis of velocity-pressure reduced order models.
In equation (35), we wrote explicitly the dependence of matrix A on the parameter vector µ. Such dependence is more or less evident for the different type of parameters. For instance, if the Reynolds number is the only parameter, i.e. µ = µ 1 = Re, from (7) we have ν = U L/Re and matrix A can be written as: with a linear dependence on Re −1 . On the other hand, if H is the only parameter, i.e. µ = µ 1 = H, the dependence of A on it is hidden in the inner products and differential operators. This holds true in general for geometric parameters. Let ξ ∈ D be a geometric parameter. If the geometric transformation T : D × Ω → Ω is affine, it is possible to express the inner products as a linear combination of the inner products on the reference domain: Only functions Θ i depend on ξ and need to be evaluated online. Matrices A i are assembled offline since they do not depend on ξ. Thus, the affine decomposition (36) allows for important computational time savings. In this work, we will consider only affine decompositions. If A depends nonlinearly on ξ, it has to be computed from scratch for each value of ξ. The efficient assembling of A when the geometric transformation is non-affine is still an active research area, one of the most popular techniques being the Empirical Interpolation Method [7]. The construction of a divergence-free basis set for V N when geometric parameters are considered is less trivial than in the case of physical parameters only. Thus, it requires further explanation. The Piola transformation P can be seen as the composition of the map T in eq. (14) with any function f defined on the image (or preimage) of T . For example, if f : Ω → R, a new function g : Ω → R can be obtained by considering g(x) = f (T −1 x) for x ∈ Ω . The Piola transformation P acts as a map between finite dimensional Hilbert spaces D×V N ( Ω ) and V N (Ω ): Its use in an offline-online setting is as follows: 1. The snapshots {u N (µ i )} N i=1 are divergence-free on the original domain Ω . By pulling back the divergence operator to the reference domain Ω through the Piola map, we obtain a set of snapshots that are divergence free on the reference domain.
2. Perform the POD or GS orthogonalization for the divergence-free snapshots on the reference domain to obtain a basis for V N ( Ω ). These basis functions are divergence free on Ω , but not on Ω unless mapped with the Piola transformation.
3. Compute the matrices A i in (36) on the reference domain and with the orthogonal divergencefree basis set.
4. During the online phase, apply the Piola transformation to the matrices A i computed at step 3 so that their entries coincides with the Piola-transformed divergence-free basis functions computed on Ω .
We refer to [9] for details on the Piola transformation, and to [46] for an application to RB methods in incompressible fluid mechanics in laminar regime. For the application of RB methods to moderately turbulent flows we refer for example to [45], and references therein.
Regarding the boundary conditions, the global support of the RB modes does not allow to impose pointwise values for the non-homogeneous Dirichlet condition. An equivalent way to impose the desired flow conditions is to impose the mass flow rate, instead of the inflow velocity profile. The physically correct definition of mass flow rate is: where u x is the x-component of the velocity. Notice that due to incompressibility and the prescribed boundary conditions, (37) is equivalent to: From the implementation point of view, it is more convenient to impose the integrated mass flow rateẇ x , that for a given inlet velocity profile is defined as: We impose the average mass flow condition for the RB simulation through a Lagrange multiplier approach as follows. We compute the integrated mass flow rate for each of the RB functions: and collect all theċ i in a vector C ∈ R N . Let α k ∈ R be the Lagrange multiplier associated with the mass flow rate constraint at the k-th iteration of the fixed point method described above. Notice that this is a new unknown in the problem. Then, instead of solving system (34), at each fixed-point iteration we solve the following linear system: We remark that imposing a constrainted condition by a Lagrange multiplier is fairly common in the Reduced Basis context, see, e.g., [51].

Bifurcation detection
In the configuration described in Section 2, the first pitchfork bifurcation point is determined by a classical modal stability analysis, that can be set up as follows. Let us consider the 3D case, for which the parameters are Re 3D and H. Suppose that an initial RB solution u N (µ i ) of the steady state problem is known for a given value of the parameter µ i = (Re 3D,i , H i ), characterized by a sufficiently small Reynolds number Re 3D,i so that the solution is surely unique. We proceed as follows: set s = 1 and Re s = Re 3D,i , then: 1. Keeping fixed the value of the geometric parameter H i , increase the value of the Reynolds number by a sufficiently small increment ∆Re (i.e., small enough so that the corrector step will converge to a solution in the desired branch) and set µ s+1 = (Re s + ∆Re, H i ).
2. Compute the RB solution u N (µ s+1 ) of the steady state problem for the new parameter value µ s+1 .
3. Compute the Galerkin projection of operator L defined in (10) on the RB space V N to form the matrix L(µ s+1 ): 4. Compute the eigenvalues of L(µ s+1 ) and check if there is one eigenvalue that has changed sign with respect to the previous iteration. If not, set s = s + 1 and go back to step 1.
We remark that the above algorithm may be unstable in the sense that in a neighborhood of the bifurcation point it may abruptly switch the approximated solution branch, or fail to converge. To make sure that the approximation is always laying on the correct branch a continuation method may be used.
Continuation methods rely on a predictor-corrector iteration to compute solutions lying on the same branch. Suppose that the u(µ s ) is a solution of equations (8)- (9), and is known to lie on a certain branch of interest. The predictor step consists in the computation of an initial guess for the velocity increment ∆u due to an increase of a single parameter, denoted by ∆µ, by solving the linearized Navier-Stokes equations. Then, starting from the predictionũ = u(µ s )+∆u, a Newton iteration is set up to impose that the new solution u(µ s + ∆µ) solves the original problem (8)- (9), under the constraint that the solution be orthogonal to the tangent plane at the point (µ s , u(µ s )) in the parameter-solution space. We refer to e.g. [23] for an introduction to continuation methods in fluid mechanics.
The continuation method is computationally quite expensive. Thanks to the fact the the GLC collocation strategy keeps the number of sample values small, the number of times the continuation method has to be restarted is reduced, allowing for important computational time savings.
Note that the matrix L is dense but has rank equal to N , with N of the order of a few tenths at most. Hence all the eigenvalues can be computed inexpensively with QR iterations [32], for instance. If the spectrum analysis had to be carried out on the full-order model, only a few of the eigenvalues closer to zero could be computed. Moreover, the computations would be much more expensive, requiring Krylov subspace methods [61] and most likely a supercomputer.
For the 2D case, we use an analogous algorithm, the only difference being that the parameters are Re 2D and λ.
Lately, increasing attention has been devoted to eigenvalue calculation (as bifurcation detector tool) at the reduced order level [47,34]. We refer to [30] for a theoretical analysis of bifurcation detection techniques in Navier-Stokes equations and to [21] for a bifurcation detection method in a similar geometry.

Results
In this section the method described in section 3 will be validated against benchmark problems reported in [52,24]. We start with the test cases in two dimensions and then consider problems in three dimensions. We show that our RB method successfully captures the bifurcation points reported in [52,24]. We compare our results with full order solutions and provide an estimate of the computational savings. Moreover, we carry out an extensive set of simulations that will allow us to confirm that the eccentric mitral regurgitant jets occur when the regurgitant orifice is long (large aspect ratio) and narrow (large expansion ratio).

2D case: one parameter study
We start with the validation of the bifurcation detection method presented in section 3.5 for the 2D test case with the Reynolds number as the only varying parameter. For the moment, the geometry is kept fixed. We set the expansion ratio λ to 15.4 in order to compare our results with those reported in reference [52]. In this case we choose a mesh with 308 spectral elements of order 11, with careful refinement near the re-entrant corners of the domain, where we can expect a loss of regularity for the solution.
As shown in table 1, we sample nine values for the Reynolds number in the interval Re 2D ∈ [0.01, 90]. For the first four values of Re 2D in table 1 the offline solver returned only the symmetric solution, as expected. For the remaining five values, the solver returned two snapshots: one for the symmetric solution (unstable [63]) and one for the asymmetric solution (stable). As mentioned in Sec. 2.1, at a Reynolds number higher than the critical value for the symmetry breaking two stable solutions co-exist, which are one the mirrored image of the other with respect to the horizontal axis (see, e.g., [8]). Bifurcation theory allows to clarify the nature of the multiplicity of possible flows, whereas a (numerical or laboratory) experiment will give one or the other of the stable symmetric solutions. Thus, for the multi-parameter case we will disregard the symmetric unstable solution and retain only the stable solutions.
The online phase for the 2D problem is performed with a RB space of dimension N = 9. In Fig. 5 we report representative snapshots for the 2D case, corresponding to Reynolds numbers Re 2D = 0.01, 13.2, 27.7, 62.2. For very low Reynolds number the solution is characterized by symmetry about the horizontal axis and a vertical axis, with a couples of vortices both upstream and downstream of the contraction called Moffatt eddies [49]. See Fig. 5(a). As the inertial effects of fluid become more important (i.e., as Re 2D increases), the Moffatt eddies upstream of the contraction gradually diminish in size and two recirculation regions of equal size develop downstream of the expansion. See Fig. 5 (b). Symmetry about the vertical axis is lost, but the solution is still symmetric about the horizontal axis. Past the bifurcation point we can see two solutions: a symmetric one (unstable) and a slightly asymmetric one (stable). See Fig. 5(c) and (d). The formation of stable asymmetric vortices in 2D planar expansion is attributed to the Coanda effect (see [70]): an increase in velocity near one wall will lead to a decrease in pressure near that wall and once a pressure difference is established across the channel it will maintain the asymmetry of the flow. This asymmetric solution remains stable for a certain range of Re 2D and asymmetries become stronger with the increasing Reynolds number. See Fig. 5(e).  To test the bifurcation detection method described in section 3.5, we run the online solver parametrized using the 2D basis set with N = 9 snapshots. In Fig. 6, we plot the real part of the eigenvalue of matrix L in (41) responsible for the symmetry breaking. We see that the curve crosses the horizontal axis at a Reynolds number of about Re 2D,sb = 26. This is in good agreement with the critical values for the symmetry breaking reported by [52,54] (Re 2D,sb = 28) and [48] (Re 2D,sb = 27).   [20,55]. The arrows indicate the direction of the increasing Reynolds numbers. Fig. 7(b) is a zoomed-in view of Fig. 7(a), and Fig. 7(c) is in turn a zoomed-in view of Fig. 7(b). In Fig. 7 (c) we see the eigenvalue responsible for the bifurcation: it is the simple eigenvalue colored in blue that changes sign as the Reynolds number increases.
For this 2D case, the computational savings are significant. The detection of the bifurcation point using the continuation method required about 80 runs, with a total computational time of around 5 minutes (0.08h) on a common desktop computer, which means 3.75s per online single run. The same computations using the full order model described in sec. 3.1 would have required about 10 CPU-hours per run. Hence, adding to the online cost the time required for the RB spaces generation (i.e., the 2 CPU hours required by the POD computations), we can estimate that the computational cost for the reduced model is around 11.5% of the computational cost for the full order model, considering all the operations needed for the bifurcation detection and computation (N = 9).
Time to build the RB spaces + Online time to detect the bifurcation point Time of the equivalent full order computation = 9 · 10h + 2h + 0.08h 80 · 10h 11.5%.
More generally, if only the online runtimes are considered, the computational savings become much more relevant compared with the offline runtimes per single query: RB online query time Equivalent full order single computation = 3.75s 10h 0.01%.
An important quantity to be used as indication if a reduced computational model is competitive is the break-even, comparing all the offline computational times needed to prepare the reduced basis problem (N = 9) and an online query with full order model: All full order computations for RB prep.
suggesting that the use of Reduced Order Methods becomes more and more competitive as the number of queries increases (with 10 or more queries this approach brings already important computational advantages). Also, this confirms that for the one-parameter scenario, this method  could be efficiently adapted to a real-time query tool to be used, e.g., on smartphones or other mobile devices with appropriate apps. We conclude this section by showing that there is no visible qualitative difference between the solutions obtained with the full order method and with the RB method for values of the Reynolds number not associated with the snapshots. See the comparison in Fig. 8

Unstable solution branch
As already mentioned, for a given expansion ratio λ and given Re 2D the symmetric flow configuration exists regardless of whether Re 2D is smaller or grater than the critical value Re 2D,sb for the bifurcation. Indeed, the symmetric branch is the only solution branch existing for Reynolds numbers below Re 2D,sb , but for Reynolds numbers above Re 2D,sb it becomes unstable. See, for example, the unstable symmetric flow configuration for Re 2D = 27.7 in Fig. 5(c) and the corresponding stable asymmetric configuration in Fig. 5(d).
The numerical tests have shown that the RB approximation of the unstable branch can be achieved, but some care is required with the choice of the trial and test RB spaces. One way to reconstruct the unstable branch is to use only the basis functions coming from the sampling of the unstable branch itself for both trial and test spaces. If this strategy is adopted, all the flow configurations of the unstable branch will be correctly approximated, but the bifurcation point will not be detected. On the other hand, if basis functions coming from both the symmetric and the asymmetric branch are employed, the bifurcation point can be successfully detected but the fixed point scheme fails to converge after the bifurcation point, oscillating without damping between the two solution branches. In this case, convergence to the symmetric or asymmetric branch after the bifurcation point can be achieved through e.g. a predictor-corrector or a pseudo-arclength continuation method (see [23]) during the online phase, with a further programming effort.
On the other hand, if one is interested only in the approximation of the stable solution branches, there is no need for basis functions coming from the unstable branch and no need for a continuation method in the online phase. The reduced basis for the velocity is constructed only with basis functions arising from stable branches. This will allow to detect the bifurcation point and compute the stable solution for every parameter value.

2D case: two-parameter study
In this section, we still consider a slightly modified 2D channel: the part of the channel upstream of the sudden expansion in Fig. 2 is removed, since we focus now on the flow downstream of the contraction. The new geometry is thus a rectangle. We let vary both the Reynolds number and the contraction width, so the parameter vector has now two components: µ = (Re 2D , λ).
In this very simple case, the geometry can be parametrized in two different ways: -Geometric parametrization: the contraction width is treated as explained in section 3 and the incompressibility constraint can be imposed through the Piola transformation as explained in section 3.4; -Boundary condition parametrization: the different aspect ratio of the contraction is imposed by parametrizing the boundary conditions. Indeed, a channel with a contraction of width λ will produce in our model a parabolic inner velocity profile dependent on λ: where the y coordinate has origin on the symmetry axis of the contraction.
One advantage of the second strategy is that the RB functions are automatically divergence-free and the relatively complex procedure of the Piola transformation does not need to be performed. Thus, we choose the boundary condition parametrization. However, we need to be careful in imposing the inlet velocity profile because the mass flow rate constraint as expressed in equations (39) and (40) is not sufficient to ensure uniqueness of the RB solution. One possible workaround for this issue is to split the boundary integral (38) used for the mass flow rate constraint in two parts: where Γ in is the part of ∂Ω where the inlet velocity profile is imposed, Γ 0 the part of Γ in where v x = 0 and Γ λ the part of Γ in where v x = 0. Notice that Γ 0 ∪ Γ λ = Γ in and Γ 0 ∩ Γ λ = ∅. We introduce two Lagrange multipliers α 0 for Γ 0 and α λ for Γ λ , in order to enforce (42) in integral form as: Finally, the two Lagrange multipliers α 0 and α λ are treated as additional unknowns, and a linear system analogous to that in equation (40) is solved. The GLC collocation sampling has been carried out on the kinematic viscosity set ν ∈ [1.5, 5] · 10 −3 and on the contraction width set w c ∈ [1/10, 1/2]. We obtained 6 values for the kinematic viscosity ν = 1.5 · 10 −3 , 1.73446 · 10 −3 , 2.375 · 10 −3 , 4.125 · 10 −3 , 4.76554 · 10 −3 , 5 · 10 −3 and 7 values for the expansion ratio λ = 2, 3, 4, 5, 6, 8, 10, so N = 42. Note that the sampling has not been performed directly on the Reynolds number due to its dependence on the contraction width. In table 3, we report the critical Reynolds numbers for the symmetry breaking computed with the RB method for different values of the contraction width. The same values in table 3 are plotted in figure 9. We clearly see that as the aspect ratio λ decreases, the critical Reynolds number for the symmetry breaking increases, as observed also in [24]. We see that Re 2D,sb decreases fast for small values of λ, while it decreases mildly for λ ≥ 6 (also recall that for λ = 15.4 we found Re 2D,sb ≈ 26). We remark that also in this case the results match closely the ones in [24]. Let us analyze the computational time savings allowed by our RB method. Since we have N = 42 with a 10h cputime needed per single run, the computational time analysis is given by: Time to build the RB spaces + Online time to detect the bifurcation point Time of the equivalent full order computation = 42 · 10h + 2h + 7 · 0.08h 7 · 10h · 80 7.5%.
The ratio between a single online reduced order run and a single full order one are the same as the one considered in the single parameter case (order 10 −4 ). The break-even, comparing all the offline computational times needed to prepare the reduced basis problem (N = 42) and an online query with full order model is: All full order computations for RB prep. Full order one query comp. time = 42 · 10h + 2h 10h 42.2.
After 43 queries a reduced order computational model brings savings. For λ = 6, which is one value among those listed in table 3, we plot in figure 10 the vertical component of the velocity is taken on the horizontal axis, at distance 1 from the inlet, versus the Reynolds number. This bifurcation diagram with both the stable and unstable solution branches compares very well with the one in [24], but it has been obtained at a fraction of the computational time as explained above. Keeping λ = 6, we check how the flow evolves as Re 2D is pushed to a higher value, well beyond the parameter range considered in this work. Fig. 11 reports the streamlines of both the unstable and stable solution at Re 2D = 600. The stable solution in Fig. 11 (b) shows that the flow structure becomes more complex, with existing recirculations changing shape and growing in size. This is consistent with the results presented in [24,54].

3D case
The three-dimensional channel in Fig. 3 has been obtained by extruding the two-dimensional geometry in Fig. 2 along the z-axis. Thus, for the 3D case, we would have three parameters: the Reynold number Re 3D , the contraction width, and the channel depth. However, since we have already investigated in Sec. 4.3 the influence of the expansion ratio (i.e., the contraction width) on the critical Reynolds number for the symmetry breaking, we fix the contraction width and consider the Reynolds number and the channel depth as the only parameters.
We set the expansion ratio λ to 15.4, due to the richness of flow patterns described in Sec. 4.1 and reference [52]. Of course, we expect the vortex structure to be much more complex than in the 2D case. We are interested in understanding how varying the Reynolds number and the aspect ratio AR (and thus H) affects the flow in the expansion channel. The goal of this section is to evaluate the effect of the walls on the bifurcating phenomenon. Intuitively, when the walls are very far apart (large values of H), their influence on the central region of the channel will be quite small, and the flow pattern can be expected to be close to the 2D case. On the other hand, when the walls are very close with respect to the channel height (small values of H), a relatively large fraction of the sectional area will be occupied by low velocity fluid. As a result, we can expect that the bifurcation will take place at higher Reynolds numbers. As reported in tables 1 and 2, we sample sample nine values for the Reynolds number in the interval Re 2D ∈ [0.01, 90] and eight value of H. Notice that the eighth "value" of H in table 2 corresponds to the 2D case.
In order to show the sequence of events as the Reynolds number is increased when the aspect ratio is fixed, we set it to 1.6398 which corresponds to H = 0.6210. In Fig. 12, we display the streamlines on the xy-plane for different values of the Reynolds number Re 3D . At Re 3D = 0.01, the 3D flow looks similar to the 2D flow: (compare Fig. 12 (a) with Fig. 5(a)) but it features smaller Moffatt eddies. As the Reynolds number increases, "lip vortices" form, as shown in Fig. 12 (b). This is in agreement with the observations in [52] and references therein. The size of the lip vortices increases as Re 3D increases and once they reach the corner, the vortices continue to grow in the downstream direction, i.e. along the x-axis. See Fig. 12 (c), (d), and (e). By convention, once they expand in the downstream direction they are called "corner vortices". Notice that the flow downstream of the expansion is symmetric about the xz-plane up to Re 3D = 76.821, while asymmetries in 2D (i.e., for H = 1) arise around Re 2D = 26.
Let us consider the geometry with H = 0.9517, which corresponds to the largest aspect ratio among those in Table 2 for which we have an actual 3D geometry. We proceed with the computation of the symmetry breaking bifurcation point using the bifurcation detection method described in section 3.5. Since H is fixed, we consider a total of 9 basis functions for the online computation, corresponding to the different values of Reynolds number reported in table 2. In Fig. 13, we plot the real part of the eigenvalue of matrix L in (41) responsible for the symmetry breaking. We see that the curve crosses the horizontal axis at a Reynolds number of about 35. This coincides with the critical value for the symmetry breaking reported by [52]. For the sake of completeness, in figure 14 we report the path of all the eigenvalues in the complex plane.
For λ = 15.4, the critical Reynolds number for the symmetry breaking in the 2D geometry (i.e., H = 1) found in Sec. 4.1 is Re 2D,sb = 26. See Fig. 6. When H is decreased to 0.9517, the critical Reynolds number for the symmetry breaking increases to Re 3D,sb = 35, as shown in Fig. 13. If H is further decreased to 0.6210, we saw in Fig. 12 that the flow remains symmetric up to Re 3D = 76.821. As expected, at low values of H the proximity of vertical walls make the flow fully three-dimensional (instead of quasi-2D) inhibiting the symmetry breaking. Thus, as H decreases Re 3D,sb becomes larger and larger.
Next, we let both the geometric parameter H and the Reynolds number vary. We display in figure 15 the streamlines on the xy-plane (left) and yz-plane for representative values of the two parameters. For low values of H and Re 3D , the flow develops without forming vortices, with the streamlines deviating only slightly out of plane. See Fig. 15(a) and (b), which have been obtained for Re 3D = 27.79 and H = 0.8165. As the channel increases in width, the streamlines gradually become fully three-dimensional, especially in the vortex region. See Fig. 15(c) through (f). Notice how Fig. 15 (c) and (d), obtained for H = 0.9517 and Re 3D = 27.79, differ from Fig. 15(a) and (b), obtained for the same Reynolds number but a smaller H. The corner vortices in Fig. 15(e) looks similar to the recirculations observed in 2D (see, e.g., Fig. 5(b)). However, in a 3D geometry the presence of a top and bottom bounding wall leads to complex 3D spiraling recirculation structures [18,66], as shown in Fig. 15(f). See also Fig. 15(g) and (h). Fig. 16 shows that streamlines on the xy-plane and yz-plane of both the unstable (symmetric) and the stable (asymmetric) solution for a value of the Reynolds number (Re 3D = 76.82) past the bifurcation point, the usual expansion ratio λ = 15.4, and H = 0.9517. The vortex pattern becomes even more intricate after the bifurcation point, with the vortices promoting the mixing    Each 3D full order computation requires about 240h of CPU time, and the preprocessing time is about 40h. The computational time savings estimate for the two parameter (Re and H) case is given by: time to build the RB spaces + online time to detect the bifurcation point time of the equivalent full order computation = 56 · 240h + 40h + 0.05h 7 · 7 · 10 · 240h 11.4%.
where, based on the experience acquired with the 2D case, we suppose that 7 runs per each parameter are required to have a reasonable tracking of the bifurcation points in the parameter space. With 2 parameters this amounts to 49 runs, each run requiring on average 10 full simulations. Thus, in the 3D case the break-even is given by: All full order computations for RB prep. Full order one query comp. time = 56 · 240h + 40h 240h 56.2.
The interpretation of this result is that a reduced order model can be expected to bring savings if more than 56 runs are planned.
To test our method, we select a geometric aspect ratio not considered in the sampling phase, and we try to recover some characterizing flow features as a function of the Reynolds number. We where Π c is any plane crossing the contraction section and orthogonal to the channel axis and |Ω ∩ Π c | is the measure of the intersection between the plane Π c and the domain Ω . We also (a)   consider the normalized axial velocity gradient: We plot the normalized axial velocity (45) and normalized axial gradient (46) along the center line for different values of the Reynolds number in figure 17 (a) and (b), respectively. The results are in good qualitative agreement with those reported in [52]. Concerning the normalized axial velocity, for small Reynolds numbers the curve is almost a symmetric step function, since the viscosity is sufficiently high to avoid large velocity gradients both inside the cross-section and along the channel length. As the Reynolds number is increased, the curve becomes more and more asymmetric, and the averaging effect of the viscosity takes longer to smooth out the velocity gradients. This is visible from the long tail of the curves with higher Reynolds number. The viscosity has also a clear effect on the normalized axial gradient in Fig. 17 (b): the two spikes show that the velocity gradients in proximity of the variations in channel width increase as the Reynolds number increases. We remark that the graphs in figure 17 can be easily drawn by saving the normalized axial velocity and normalized axial gradient for the RB functions and using these as to interpolate the desired output in real time. This feature is particularly interesting in the real-time query case, since it does not need to search a large database during the postprocessing phase.
We conclude the section with the streamlines for the flow associated to H = 0.2085 (in Fig. 18), H = 0.6210 (in Fig. 19), and H = 0.9517 (in Fig. 20) for a small value, a medium value, and a large of Re 3D ∈ [0.01, 90]. In particular, compare the solutions for Re 3D = 90 (leftmost panel in Fig. 18, 19, and 20). They clearly show that at low values of H the symmetry breaking bifurcation is pushed to higher values of Re 3D due the vertical walls.
Based on the results presented in Sec. 4.3 and 4.4, we conclude that eccentric mitral regurgitant jets are produced by long (large H) and narrow (large λ) orifices. In fact, such slender orifices associated with eccentric jets, seem to resemble the coaptation geometry of the mitral valve. Coaptation is the region where the two leaflets of the mitral valve meet Our hypothesis is    that Coanda effect occurs in mitral valves in which the leakage, i.e., regurgitation, occurs along a large section of the coaptation zone, rather than at an isolated point, leading to a possibly significant regurgitant volume. This is corroborated by clinical observations indicating that eccentric regurgitant jets are, indeed, prevalent in patients with severe MR [43,15,62].
Before the study presented in this manuscript, our collaborators at the Houston Methodist DeBakey Heart & Vascular Center had never succeeded in reproducing the Coanda effect in vitro. Following our results, they designed a long and narrow orifice in a divider plate that mimics a closed leaky mitral valve. A close-up view of the orifice is in Fig. 21(a). The divider plate was mounted on an anatomically correct mock (left) heart chamber developed to study the use of 2D and 3D color Doppler techniques in imaging the clinically relevant intra-cardiac flow events associated with regurgitant jets [44,42]. See Figure 21 (b). The chamber is connected to a pulsatile flow loop. The fluid in the mock heart chamber is water with 30% glycerin added to mimic blood viscosity. Notice that this is consistent with modeling blood as a Newtonian fluid in Sec. 2. From the 2D Doppler echocardiographic image in Fig. 21 (c) we see that indeed the slender orifice in Fig. 21(a) generates a regurgitant jet that hugs the wall. See also [69]. We expect also that these studies could enhance in the near future in vivo studies and applications.

Conclusions and perspectives
The symmetry breaking bifurcation (Coanda effect) has been studied in parametric flows, representing a simplified test case for regurgitant mitral valve flows. Our preliminary work shows that standard reduced order methods (e.g., Reduced Basis and /or Proper Orthogonal Decomposition) allow to capture complex physical and mathematical phenomena, such as bifurcations in the parametrized Navier-Stokes equations, at a fraction of the computational cost required by full order order methods. In order to detect the bifurcation points, the reduced parametric Navier-Stokes equations have been supplemented with a generalized eigenvalue problem, also cast into the reduced order setting. This work is also an example of computational collaboration between high performance computing and reduced order methods: thanks to the computational gains with the same resources we can treat more complex problems. This computational collaboration has demonstrated the ability to provide reliable and accurate results with significant reduction of computational times. Results have been validated both with the full-order model and by comparison with parametric studies available in literature for both 2D and 3D cases.
Research perspectives in this field include the development of proper error bounds for the detection of the bifurcation points and the verification of the accuracy. At the state of the art this aspect is carried out by supplementing the state equation with a generalized eigenvalue problem, solved with the same reduced order method proposed for the state equation. Moreover, we plan on taking into account the interaction of the fluid with elastic walls (i.e., elastic valve leaflet) [5].
This would lead to important improvements in the study of this complex multiphysics nonlinear problem and a better understanding of how the Coanda effect is influenced by the valve elasticity.