Reduced Numerical Approximation of Reduced Fluid-Structure Interaction Problems With Applications in Hemodynamics

1. Introduction

When modeling hemodynamics phenomena in big arteries, the resulting model is a complex unsteady fluid-dynamics system, usually coupled with a structural model for the vessel wall. In specific cases, suitable assumptions can be made to reduce the complexity of the model equations. In particular, when the displacement is small, the moving domain can be linearized around a reference steady configuration and the dynamics of the vessel motion can be embedded in the equations for the blood flow. In such way we obtain a Reduced Fluid-Structure Interaction (RFSI) formulation where a Navier-Stokes system in a fixed fluid domain is supplemented by a Robin boundary condition that represents a surrogate of the structure model.

Although the RFSI model is faster with respect to fully three-dimensional (3D) models where the structure is solved separately, the numerical computation of one heartbeat is still expensive: the resolution of an entire heartbeat, that typically lasts one physical second, takes orders of hours of computational time on a supercomputer. A big challenge in realistic applications is to achieve a real time resolution of fluid-structure interaction problems. In particular, in hemodynamics applications, this would grant the possibility to perform real time diagnosis. Nevertheless, the great variability of patient-specific data requires the parametrization of the model with respect to many physical and geometrical quantities. Moreover, as we have recalled above, the complexity of the hemodynamics phenomena requires a mathematical description with complex unsteady models that are difficult to be solved in real time. The RFSI model is already a simpler version of the fully 3D FSI system; a further reduction of the physical model would result in an inaccurate estimation of specific outputs, like e.g., the wall shear stress when using a rigid wall model [1, 2]. Thus, to further reduce numerical costs, in this work, we focus on the reduction of the discretization space. In realistic applications, the finite element space has order of 10⁶ degrees of freedom. The aim is to construct a discretization space such that the number of degrees of freedom is reduced to less than 100 and then to be able to solve one heartbeat in 1 s.

In the past few years, due to their relevance in realistic applications, a lot of interest has been devoted to discretization reduction techniques for parametrized Partial Differential Equation (PDE) problems (e.g., [3–6]). These techniques aim to define a suitable reduced order model which can be solved with marginal computational costs for different values of the model parameters. Reduced order models are then important in the many query context, when a parametrized model has to be solved for different values of the parameter, and in the real time problems, when the solution has to be computed with marginal computational costs. To obtain a suitable reduced order model, we typically start from a problem written in a high-fidelity approximation framework, e.g., using the finite element method. The dimension of the discretized system is then drastically reduced through suitable projection operators. The construction of these projection operators is the core of the reduced order technique. Another key concept in the reduction framework is the subdivision of the computational costs into two stages: an offline stage, expensive but performed once, and an online stage, real time and performed each time new values of the model parameters are considered. During the offline stage the projection space is generated by a reduced basis of functions of the high-fidelity approximation space.

Reduced order models applied to the Burgers equation parametrized with respect to the Péclet number is considered in Yano et al. [7] and Nguyen et al. [8]. Other applications of reduced basis techniques applied to fluid problems can be found (e.g., in [9–18] and in the recent volume [6]).

The aim of this work is indeed to propose a suitable discretization reduction algorithm that can be applied to a RFSI problem. The work is organized as follows. In section 2 we present the partial differential equations that we are interested in solving. We propose a possible parametrization of the unsteady equations with respect to temporal varying data and with respect to a perturbation of the boundary data. In section 3 we then present how the standard proper orthogonal decomposition algorithm can be applied to the problem at hand in order to generate a suitable reduced space. Moreover, we propose a way to improve the quality of the reduced approximation based on a greedy procedure. Finally, in section 4 we apply the reduction algorithms presented to a realistic hemodynamic problem. Conclusions follow.

2. Model Equation

Blood is in large vessels can be modeled as an incompressible viscous fluid the well-known Navier-Stokes equations (e.g., [19, 20]). Being [0, T] the temporal interval of interest the Navier-Stokes system reads as follows:

$\{\begin{array}{cc}{\rho }_f\frac{\partial u}{\partial t}+{\rho }_f(u·\nabla )u-\nabla ·{\sigma }^{n_S}=0& \text{in }\Omega \times [0,T],\\ \nabla ·u=0& \text{in }\Omega \times [0,T]\end{array}$

where u and p are the velocity and pressure of the blood, respectively, and ρ_f is its density. We denote by ${\sigma }^{n_S}$ the Cauchy stress tensor

${\sigma }^{n_S}={\mu }_f(\nabla u+{(\nabla u)}^T)-pI,$

with I being the identity tensor and μ_f the blood dynamic viscosity. Ω denote the domain of interest, in our case, the lumen of the vessels where we are interested in computing the dynamics of blood. Due to the compliant vessel wall, Ω should be time dependent. Considering that the wall displacement is relatively small with respect to the arterial diameter, we assume Ω as fixed which allows us to reduce the computational complexity otherwise generated from a moving domain. Nevertheless, to retrieve the physical effect of the wall compliance we introduce a non-rigid boundary condition on the lateral surface of the lumen. The condition is derived by a three dimensional linear isotropic elastic model condensed as a two dimensional membrane [1, 21]. Denoting Γ the lateral surface (i.e., the fluid-structure interface), Π_Γ(d_s), the stress-strain relation of this membrane, can be written as:

$\begin{array}{l}{\Pi }_{\Gamma }(d_s)=h_s\frac{E_s{\nu }_s}{(1-2{\nu }_s)(1+{\nu }_s)}\text{tr}\left(\frac{{\nabla }_{\Gamma }{d}_s+{({\nabla }_{\Gamma }{d}_s{\displaystyle )}}^T}{2}\right)\\ +h_s\frac{E_s}{2(1+{\nu }_s)}\left({\nabla }_{\Gamma }{d}_s+{({\nabla }_{\Gamma }{d}_s)}^T\right).\end{array}$

where ∇_Γd_s is the tangential gradient of d_s, E_s is the structural Young modulus, ν_s is the Poisson's ratio and h_s is the material thickness. All the physical parameters of the structure are assumed homogeneous in space.

Let us now suppose that the boundary ∂Ω is divided into three non intersecting parts such that ∂Ω = Γ ∪ Γ_D ∪ Γ_N. Γ_D is the Dirichlet boundary, typically the inflow of a vessel, Γ_N is the Neumann boundary, typically the outflows. We introduce the Hilbert space $V=H^1(\Omega ;\Gamma )\text{\hspace{0.22em}}=\{v\in H^1(\Omega )v{|}_{\Gamma }\in H^1(\Gamma )\}$ and the correspondent vectorial spaces V = [V]³ and V = [V]³. Moreover, we introduce a suitable couple standard finite element spaces V_h and Q_h such that V_h ⊂ V and $Q_h\subset L^2(\Omega )$ and they represent a stable coupled of finite element spaces for the Navier-Stokes equations. We set X_h: = V_h × Q_h. We define [t₀ T] a time interval of interest and we divide it into subintervals [t_n t_n+1] for n = 0, .., N_T − 1 such that t₀ < t₁ < t₂ < … < t_NT = T and t_n+1 − t_n = Δt; let us define ${{N}}_T=\{0,1,\dots ,N_T,N_T\}$ the collections of all the temporal indexes n. For a generic function ϕ(t) we use ${\varphi }^n:=\varphi (t_n)$. Finally, we define the operators D(·) and D_Γ(·) as follows:

$D(v)=\frac{\nabla v+{(\nabla v)}^T}{2}\text{ and }{D}_{\Gamma }(v)=\frac{{\nabla }_{\Gamma }v+{({\nabla }_{\Gamma }v)}^T}{2}\forall v\in V,$

where ∇(·) is the standard gradient operator and ∇_Γ(·) is the tangential component of the gradient with respect to the surface Γ.

The RFSI model as presented in Colciago et al. [1] is an unsteady Navier-Stokes model set on a fixed domain with generalized Robin boundary conditions (For similar models see e.g., [21–24]). Let us introduce the velocity and pressure unknowns [u_h, p_h] and the corresponding test functions [v_h, q_h]. Although the RFSI model lives in a fixed domain, it is necessary to define an auxiliary variable which stands for the displacement of the arterial wall d_{s, h}. Using a backward Euler finite difference method for the time derivatives, the fully discrete weak formulation of the RFSI problem is written as follows:

for each n = 0, .., N_T − 1, find $\left[{\text{u}}_h^{n+1},p_h^{n+1}\right]\in {\text{X}}_h$ such that ${\text{u}}_h^{n+1}={\text{g}}_D^{n+1}$ on Γ_D and

$\begin{array}{l}{a}_0([u_h^{n+1},p_h^{n+1}],[v_h,q_h])+a_1(u_h^n,u_h^{n+1},v_h)=F_0(v_h;h^{n+1})\\ F_u(v_h;u_h^n)+F_{d_s}(v_h;d_{s,h}^n)\forall [v_h,q_h]\in X_h,& (1)\end{array}$

where

$\begin{array}{l}{a}_0([u_h^{n+1},p_h^{n+1}],[v_h,q_h])={\displaystyle {\int }_{\Omega }(}{\rho }_f\frac{u_h^{n+1}}{\Delta t}·v_h+(2\mu D(u_h^{n+1})\\ -p_h^{n+1}I):\nabla v_h+q_h^{n+1}\nabla ·u_h)d\Omega \\ +{\displaystyle {\int }_{\Gamma }(}\frac{h_s{\rho }_s}{\Delta t}{u}_h^{n+1}·v_h\\ +h_s\Delta t{\Pi }_{\Gamma }(u_h^{n+1}):{\nabla }_{\Gamma }{v}_h)d\Gamma ,\\ a_1(u_h^n,u_h^{n+1},v_h)={\displaystyle {\int }_{\Omega }{\rho }_f}(u_h^n·\nabla )u_h^{n+1}·v_{h}d\Omega ,\\ F_0(v_h;h^{n+1})={\displaystyle {\int }_{{\Gamma }_N}{g}_N^{n+1}}·v_{h}d{\Gamma }_N,\\ F_u(v_h;p_h^n)\text{}\text{}=\text{}\text{}{\displaystyle {\int }_{\Omega }\frac{{\rho }_f}{\Delta t}}{u}_h^n·v_{h}d\Omega \text{}+\text{}\text{}{\displaystyle {\int }_{\Gamma }\frac{h_s{\rho }_s}{\Delta t}}{u}_h^n·v_{h}d\Gamma ,\\ F_{d_s}(v_h;d_{s,n}^n)=-{\displaystyle {\int }_{\Gamma }{h}_s}{\Pi }_{\Gamma }(d_{s,h}^n):{\nabla }_{\Gamma }{v}_{h}d\Gamma ,& (2)\end{array}$

with ${\text{d}}_{s,h}^{n+1}={\text{d}}_{s,h}^n+\Delta t{\text{u}}_h^{n+1}$ and ρ_s represents the density of the solid. The functions g_N(x, t) and g_D(x, t) are sufficiently regular functions that stand for the Dirichlet and Neumann boundary data, respectively. Finally the problem should be equipped with suitable initial condition that, without any loss of generality, we suppose to be equal to zero.

As said before, the RFSI problem (1) is indeed a linearized Navier-Stokes on a fixed domain with a non standard boundary condition on the interface Γ. In particular is a generalized Robin boundary condition that contains both a mass and a stiffness boundary term to mimic the presence of a compliant arterial wall surrounding the fluid domain (see [25] for more detailed on the analysis of partial differential equations with generalized Robin boundary condition). We remark that ${\text{d}}_{s,h}^n$ does not represents a problem unknown since it is indeed reconstructed as a weighted sum of the velocities at different time instants

2.1. Boundary Condition

Problem (1) is endowed with Dirichlet velocity boundary condition on the inlet surface Γ_D. Given the inlet velocity data g_D(x, t), at the time instant t_n+1 we impose:

$\begin{array}{ll}{u}_h^{n+1}=g_D^{n+1}\text{ on }{\Gamma }_D.& (3)\end{array}$

The Neumann boundary condition $\text{D}({\text{u}}_h^{n+1})\text{n}={\text{g}}_N$ is imposed weakly on Γ_N. The solution of problems (1)–(3) depends on the time variable t through the inlet and outlet data: g_N(x, t) and g_D(x, t). We suppose that

$\begin{array}{ll}{g}_D(x,t)={\sigma }_1(t){\tilde{g}}_D(x)\text{ and }{g}_N(x,t)={\sigma }_2(t){\tilde{g}}_N(x),& (4)\end{array}$

that is we separate the contribution of the space and temporal variables in the inlet and outlet data. In realistic applications, the separation of variables (4) often derives directly from modeling choices. If at the outlet we prescribe an average normal stress, no spatial variability is involved in the boundary condition data g_N. At the inlet, the Dirichlet data is imposed by means of a velocity profile; typically Poiseuille or Womersley profiles are chosen in hemodynamics applications [20]. The separation of variables in g_D(x, t) in this case is straightforward.

Assumption (4) allows to write an affine decomposition of the operators in problem (2) with respect to σ₁(t) and σ₂(t). With respect to the latter we have:

$F_0([v_h,q_h];{\sigma }_2^{n+1},g_N)={\sigma }_2^{n+1}{\displaystyle {\int }_{{\Gamma }_N}{\tilde{g}}_N}·v_{h}d{\Gamma }_{out}.$

The non homogeneous Dirichlet boundary condition (3) is not directly included in the variational form (1). In order to write the affine decomposition with respect to the parameter σ₁(t), a suitable choice to embed condition (3) into the variational formulation has to be made. In the literature two possible approaches are proposed: a strong imposition, using a lifting function or suitable Lagrange multipliers [17], and a weak imposition adding suitable penalty variational terms [2, 26]. Due to the fact that the Dirichlet data can be written in the form (4), a single time independent lifting function can be constructed and properly weighted by a scalar in order to represent the lifting at each temporal instant.

We explain how problem (1) is modified when a lifting function for the Dirichlet condition (3) is introduced. Let us directly consider the fully discretized formulation (1). We define the time independent lifting function $\text{R}\tilde{\text{g}}:{ℝ}^3\mapsto {ℝ}^3$ such that $\text{R}\tilde{\text{g}}\in {\text{V}}_h$ and

$R\tilde{g}(x)={\tilde{g}}_D(x)\text{ on }{\Gamma }_D\text{ and }R{\tilde{g}}_D(x)=0\text{ on }\partial \Omega \backslash {\Gamma }_D.$

At the time level t_n+1, the lifting function of the data ${\text{g}}_D^{n+1}={\sigma }_1^{n+1}{\tilde{\text{g}}}_D$ reads $\text{R}{\text{g}}_D^{n+1}={\sigma }_1^{n+1}\text{R}{\tilde{\text{g}}}_D$. Then, for each t_n+1, we introduce the following change of variable:

$\begin{array}{ll}{\tilde{u}}_h^{n+1}=u_h^{n+1}-R{g}_D^{n+1}.& (5)\end{array}$

We define the space X_h,ΓD as ${\text{X}}_{h,{\Gamma }_D}:={\text{V}}_h\cap {\left[H_{{\Gamma }_D}^1(\Omega )\right]}^d\times Q_h$ and we observe that ${\text{d}}_{s,h}^{n+1}=\sum _{s=0}^{n+1}\Delta t{\text{u}}_h^s=\sum _{s=0}^{n+1}\Delta t{\tilde{\text{u}}}_h^s$ on Γ.

2.1.1. Affine Decomposition

Using the definitions of the functionals as in (2), we are now ready to write the affine decomposition of problem (1) with respect to the temporal parameters σ₁(t) and σ₂(t). We remark that the lifting function $\text{R}\tilde{\text{g}}D$ does not depend on the time variable, thus the problem parameter at a fixed time level can be gathered in the following vector:

$\begin{array}{ll}{({\mu }^{\text{n+1}})}^T:\text{}\text{}=[{\mu }_0,{\mu }_1,{\mu }_2]:\text{}\text{}=[{\sigma }_1^{n+1},{\sigma }_2^{n+1},{\sigma }_1^n].& (6)\end{array}$

One single time step of finite element approximation of the RFSI problem can be written under the form:

for each n = 0, .., N_T − 1, find ${\tilde{\text{U}}}_h^{n+1}\in {\text{X}}_{h,{\Gamma }_D}$ such that

$\begin{array}{ll}a({\tilde{U}}_h^{n+1},W_h;{\tilde{U}}_h^n,{\mu }^{n+1})=F(W_h;\widetilde{U_h^n},D_h^n,{\mu }^{n+1})\forall W_h\in X_{h,{\Gamma }_D},\hfill & \left(7\right)\hfill \end{array}$

where

$\begin{array}{l}a({\tilde{U}}_h^{n+1},W_h;{\tilde{U}}_h^n,{\mu }^{\text{n+1}}):\text{}\text{}=a_0({\tilde{U}}_h^{n+1},W_h)\\ +{\mu }_2{a}_1(R{\tilde{g}}_D,{\tilde{U}}_h^{n+1},W_h)\\ +a_1({\tilde{U}}_h^n,{\tilde{U}}_h^{n+1},W_h),\\ F(W_h;{\tilde{U}}_h^n,d_{s,h}^n,{\mu }^{\text{n+1}})):\text{}\text{}\text{=}{\mu }_{\text{1}}{\text{F}}_{\text{0}}(W_h;\tilde{h})+F_u(W_h;{\tilde{U}}_h^n)\\ +{\mu }_2{F}_u(W_h;R{\tilde{g}}_D)+F_{d_s}(W_h;d_{s,h}^n)\\ -{\mu }_0{a}_0(R{\tilde{g}}_D,W_h)\\ -{\mu }_0{a}_1({\tilde{U}}_h^n,R{\tilde{g}}_D,W_h)\\ -{\mu }_2{\mu }_0{a}_1(R{\tilde{g}}_D,R{\tilde{g}}_D,W_h).& (8)\end{array}$

Due to the fact that we use a semi-implicit treatment of the convective term the formulation of the RFSI problem at one single time instant t_n+1 can be interpreted as a linear steady problem parametrized with respect to ${\mu }^{n+1},{\tilde{\text{U}}}_h^n$ and ${\text{d}}_{s,h}^n$.

Furthermore, we can introduce a parameter in the inlet flow rate function representing a small perturbation with respect to a reference value: the inlet flow rate function (4) is then modified as

$\begin{array}{ll}g(x,t;\alpha )=\theta (\alpha ,t){\sigma }_1(t){\tilde{g}}_D(x),& (9)\end{array}$

where $\alpha \in {D}$, being ${D}$ the set of the admissible value of α. The same affine decomposition 8, with the following modification: the parameter becomes ${({\mu }^{n+1})}^T:=\left[{\mu }_0,{\mu }_1,{\mu }_2,{\mu }_3\right]:=\left[{\sigma }_1^{n+1},{\sigma }_2^{n+1},{\sigma }_1^n,{\theta }^n(\alpha )\right]$ and in (8) we substitute μ₀ with μ₃μ₀ and μ₂ with μ₂μ₃.

3. Numerical Reduction

In this section we briefly introduce some of the basic concepts of the reduced basis method that are useful to our purpose. For more details on the reduced basis theory we address the interested reader to e.g., Rozza et al. [5], Hesthaven et al. [27], and Quarteroni et al. [28]. We already introduced ${\tilde{\text{U}}}_h^{n+1}$ that, at each time instant is the a high-fidelity approximation of the exact solution and is computed as a finite element solution with a sufficiently fine mesh. The solutions ${\tilde{\text{U}}}_h^{n+1}$ of problem (7) are, in general, expensive to obtain from the computational point of view, since in realistic applications the finite element spaces has order of 10⁶ degrees of freedom and the complexity of the geometrical domain does not always allow for the generation of structured meshes. We conclude that due to the magnitude of the finite element problem a real time computation would be impossible to achieve.

As in the standard reduced basis theory, we state the following assumption: the family of solutions ${\tilde{\text{U}}}_h^{n+1}={\tilde{\text{U}}}_h^{n+1}({\mu }^{n+1})$ obtained for different realizations of the parameters belongs to a low dimensional manifold ${{M}}_h^{\mu }$. The aim of the reduction techniques is to find a suitable approximation of the manifold ${{M}}_h^{\mu }$ through the construction of a low dimensional space X_N ⊂ X_h,ΓD. The dimension of the reduced space N needs to be orders of magnitude lower that the dimension of the finite element space. The reduced approximation of RFSI problem reads:

given ${\tilde{\text{U}}}_N^0={\tilde{\text{U}}}_h^0$, for each n = 0, .., N_T − 1, find ${\tilde{\text{U}}}_N^{n+1}\in {\text{X}}_N$ such that

$\begin{array}{ll}a({\tilde{U}}_N^{n+1},W_N;{\tilde{U}}_N^n,{\mu }^{n+1})=F(W_N;{\tilde{U}}_N^n,D_N^n,{\mu }^{n+1})\forall W_N\in X_N,\hfill & \left(10\right)\hfill \end{array}$

where a(·, ·) and F(·) are defined as in (8).

3.1. Proper Orthogonal Decomposition

We apply a discretization reduction to the RFSI problem (7) and the Proper Orthogonal Decomposition (POD) method. In the context of this work we only detail the specific choices performed in relation to the problem at hand, for more details about POD applied to fluid problems we address the reader to e.g., Rowley [29] and Willcox and Peraire [30].

We define a subset of temporal indexes ${{N}}_S\subset {{N}}_T$ with cardinality N_S and consider the solutions of problem (1) at the time instants $t^{n_S}$ for $n_S\in {{N}}_S$. The solutions ${\tilde{\text{U}}}_h^{n_S}$, called snapshots, represent our starting point for the POD analysis. Since the RFSI problem (7) is a saddle point problem in two variables (velocity and pressure) with different characteristic order of magnitude, we split the POD into two eigenvalue decompositions: one for the velocity variable and another for the pressure one [31]. We measure the energy associated to the snapshots using the following scalar products: for the velocity, we set

$\begin{array}{l}{\left(u_h,v_h\right)}_V:\text{}={\left(u_h,v_h\right)}_{H^1(\Omega )}+{\left(u_h,v_h\right)}_{H^1(\Gamma )},\\ \forall u_h,v_h\in V_h\subset V(=H_{{\Gamma }_D}^1(\Omega ;\Gamma )),& (11)\end{array}$

and for the pressure,

$\begin{array}{ll}{(p_h,q_h)}_Q:\text{}\text{}={(p_h,q_h)}_{L^2(\Omega )},\forall p_h,q_h\in Q_h\subset Q:\text{}\text{}=L^2(\Omega ).& (12)\end{array}$

Then, we compute the two Gramian matrices

$\begin{array}{ll}{G}_{ij}^u={(u_h^i,u_h^j)}_V\text{ and }{G}_{ij}^p={(p_h^i,p_h^j)}_Q\forall j,i\in {{N}}_S,& (13)\end{array}$

and we perform the eigenvalue decomposition of G^u and the one of G^p, obtaining the pairs $({\lambda }_k^{\text{u}},{\zeta }_k^{\text{u}})$ and $({\lambda }_k^p,{\zeta }_k^p)$ where ${\lambda }_k^{\text{u}},{\lambda }_k^p\in ℝ$ and ${\zeta }_k^{\text{u}},{\zeta }_k^p\in {ℝ}^{N_S}$ are the k − th eigenvalues and eigenvectors of the velocity and pressure Gramian matrices, respectively, for $k\in {{N}}_S$. Fixing the same tolerance for both the velocity and pressure decompositions, we select the first N^u and N^p eigenpairs such that:

$\begin{array}{ll}\frac{{\displaystyle {\sum }_{j=1}^{N^u}{\lambda }_j^u}}{{\displaystyle {\sum }_{k=1}^{N_S}{\lambda }_k^u}}\ge 1-tol\text{ and }\frac{{\displaystyle {\sum }_{j=1}^{N^p}{\lambda }_j^p}}{{\displaystyle {\sum }_{k=1}^{N_S}{\lambda }_k^p}}\ge 1-tol,& (14)\end{array}$

respectively. The j−th velocity eigenfunction ${\varphi }_j^{\text{u}}\in {\text{V}}_h$ is reconstructed using the linear combination:

${\varphi }_j^u=\frac{1}{{\lambda }_j}{\displaystyle \sum _{n_S\in {{N}}_S}{[{\zeta }_j^u]}_{n_S}{u}_h^{n_S},\text{ for }j=1,\mathrm{..},N^u.}$

Similarly for ${\varphi }_j^p\in Q_h$ for j = 1, .., N^p. We remark that, since the velocity basis are linear combinations of solutions of problem (7), they all verify $\underset{\Omega }{\int }{q}_h\nabla ·{\varphi }_j^{\text{u}}=0$, ∀q_h ∈ Q_h for j = 1, .., N^u. Thus, the linear system induced by the bilinear form a(·, ·) as in (2) would be singular if we consider the functional spaces generated from the velocity functions ${\varphi }_j^{\text{u}}$ and the pressure modes ${\varphi }_j^p$. One of the possibilities often employed in the context of Navier-Stokes equations is to restrict the system and to solve the problem only for the velocity unknown (see e.g., [32]). Unfortunately, this is not possible when considering problem (2). The generalized boundary condition applied on Γ derives from a structural model which solution is driven by the pressure condition set on the external boundary in the structural model (see [1, 22]). If we solve the reduced system not taking into account the pressure variable, we cannot recover the velocity on the boundary Γ and the output functionals that depends on these values (e.g., wall shear stress). For these reasons, following Rozza and Veroy [33], for each selected pressure mode ${\varphi }_j^p$, we define the corresponding supremizer function σ_j ∈ V_h as the solution of the following problem:

$\begin{array}{ll}({\sigma }_j,v_h)={\displaystyle {\int }_{\Omega }{\varphi }_j^p}\nabla ·v_{h}d\Omega \forall v_h\in V_h,\text{ for }j=1,\mathrm{..},N^p.& (15)\end{array}$

We then add them to the POD basis. The POD reduced space ${\text{X}}_N^{POD}$ associated to the RFSI model is composed by the basis functions ${\left\{{\psi }_j\right\}}_{j=1}^{N^{\text{u}}+2\times N^p}$, ξ_j ∈ X_h defined as follows:

$\begin{array}{l}{\psi }_j={\left[{\varphi }_j^u,0\right]}^T\text{ for }j=1,\mathrm{..},N^u\\ {\psi }_{N^u+j}={\left[0,{\varphi }_j^p\right]}^T\text{ for }j=1,\mathrm{..},N^p\text{ and }\\ {\psi }_{N^p+N^u+j}={\left[{\varphi }_j^{\sigma },0\right]}^T\text{ for }j=1,\mathrm{..},N^p,& (16)\end{array}$

where ${\varphi }_j^{\sigma }$ for j = 1, .., N^p represent the orthonormalization of the supremizer functions σ_j, obtained with a Gram-Schmidt algorithm with respect to the scalar product (·, ·)_V.

3.2. Greedy Enrichment

The bottleneck of the POD procedure is the computation of the high-fidelity solutions ${\tilde{\text{U}}}_h^n$ necessary to build the correlation matrix: we have to solve a finite element problem N_T times. Moreover if we choose ${{N}}_S={{N}}_T$, the Gramian matrix becomes too large and its eigenvalue decomposition gets too much expensive. We can envision two situations where we would like to improve the quality of the approximation obtained with the POD reduced space without changing the snapshots sample. For example, if N_S is five times smaller than N_T, the information carried by the snapshots sample refers to only the 25% of the entire set of the truth solutions. Is it possible to improve the quality of the reduced approximation, without increasing the number of snapshots selected? In another scenario, suppose that a perturbation parameter α is introduced in the unsteady problem (7), as proposed in (9), and that the snapshots are computed for a specific value of α = α₁. We would like to generate a reduced space that suitably approximates also the truth solutions for other values of α without recomputing all the high-fidelity snapshots.

With these two scenarios in mind, we propose a strategy to improve the quality of the reduced approximation based on a greedy algorithm. For references to standard greedy algorithms applied to parametrized PDEs see e.g., Hesthaven et al. [27] and Quarteroni et al. [28].

We introduce another solution ${\text{U}}_{N,h}^n$ that belongs to an intermediate problem between (7) and (10): find ${\text{U}}_{N,h}^n\in {\text{X}}_h$ such that

$\begin{array}{ll}a(U_{N,h}^{n+1},W_h;{\tilde{U}}_N^n,{\mathit{\mu }}^{n+1})=F(W_h;{\tilde{U}}_N^n,D_N^n,{\mathit{\mu }}^{n+1}) \forall W_h\in X_{h,{\Gamma }_D},\hfill & \left(17\right)\hfill \end{array}$

We notice that problem (17) is set in the high-fidelity approximation framework but the right hand side and the advection field are defined by (10). In fact, in (7), these terms are evaluated using the truth solution ${\tilde{\text{U}}}_h^n$, while in (17) it is evaluated using the reduced solution ${\tilde{\text{U}}}_N^n$, as in problem (10). Considering the error between ${\tilde{\text{U}}}_N^n$ and ${\tilde{\text{U}}}_h^n$ in a generic norm ||·||_*, the following triangular inequality holds:

$\begin{array}{l}\Vert {\tilde{U}}_N^{n+1}-{\tilde{U}}_h^{n+1}{\Vert }_{*}=\Vert {\tilde{U}}_N^{n+1}-U_{N,h}^{n+1}+U_{N,h}^{n+1}-{\tilde{U}}_h^{n+1}{\Vert }_{*}\\ \le \Vert {\tilde{U}}_N^{n+1}-U_{N,h}^{n+1}{\Vert }_{*}+\Vert U_{N,h}^{n+1}-{\tilde{U}}_h^{n+1}{\Vert }_{*}\end{array}$

The greedy procedure that we propose focuses on the first contribution $||{\tilde{\text{U}}}_N^{n+1}-{\text{U}}_{N,h}^{n+1}|{|}_{*}$. Subtracting problem (10) from (17) allows to state a result of Galerkin orthogonality:

$a(U_{N,h}^{n+1}-{\tilde{U}}_N^{n+1},W_h;{\tilde{U}}_N^n,{\mu }^{\text{n+1}})=0.$

We assume that the dual norm of the residual can be used as an indicator of the error $||{\text{U}}_N^{n+1}-{\text{U}}_h^{n+1}|{|}_{\text{X}}$. In particular, at each time level t_n+1, we consider

$\begin{array}{ll}{r}_N^{n+1}(W_h):\text{}\text{}=F(W_h;{\tilde{U}}_N^n,D_N^n,{\mathit{\mu }}^{n+1})-a({\tilde{U}}_N^{n+1},W_h;{\tilde{U}}_N^n,{\mathit{\mu }}^{n+1})\hfill & \left(18\right)\hfill \end{array}$

and its associated dual norm $||r_N^{n+1}({\text{W}}_h)|{|}_{{\text{X}}^{\prime }}$.

We now have defined all the necessary quantities, we can proceed presenting the steps to be performed when we want to enrich the POD basis with a greedy algorithm. First, perform a POD on the snapshots ${\tilde{\text{U}}}_h^{n_S}$, for $n_S\in {{N}}_S$ and we construct the reduced space ${\text{X}}_N^{POD}$. Then, we start the greedy enrichment setting ${\text{X}}_N={\text{X}}_N^{POD}$:

1. Generate the reduced basis solutions ${\tilde{\text{U}}}_N^n$, $n\in {{N}}_T$, by solving the reduced order problem (10).

2. Compute the dual norms of the residuals $||r_N^n({\text{W}}_h)|{|}_{{\text{X}}^{\prime }}$, $n\in {{N}}_T$, which are used as error indicators.

3. Select n* such that

$n^{*}=\arg \underset{n\in {{N}}_T}{\max }\Vert r_N^n(W_h){\Vert }_{X^{\prime }}.$

4. Compute the ${\text{U}}_{N,h}^{n^{*}}$ by solving the reduced order problem (17).

5. Split ${\text{U}}_{N,h}^{n^{*}}$ into its velocity and pressure components, ${\text{u}}_h^{n^{*}}$ and $p_h^{n^{*}}$, respectively. Compute the supremizer ${\sigma }_{n^{*}}$ associated with the pressure component.

6. Compute ϕ^u representing the orthonormalization of the velocity function ${\text{u}}_h^{n^{*}}$ with respect to the reduced space X_N, obtained with a Gram-Schmidt algorithm considering the scalar product (·, ·)_V; similarly for ϕ^p and $p_h^{n^{*}}$.

7. Build ${\text{X}}_{N+2}={\text{X}}_N\oplus \left\{{\psi }^{\text{u}},{\psi }^p\right\}$ defined as is (16).

8. Compute ϕ^σ representing the orthonormalization of the velocity function ${\sigma }^{n*}$ with respect to the reduced space X_N+2, obtained with a Gram-Schmidt algorithm considering the scalar product (·, ·)_V.

9. Build ${\text{X}}_{N+3}={\text{X}}_{N+2}\oplus \left\{{\psi }^{\sigma }\right\}$ defined as is (16).

10. Update the structures for the online computation of the reduced solutions and the dual norms of the residuals.

11. Set N = N + 3 and X_N = X_N+3. Repeat until a predefined stopping criterion is satisfied.

Remark. We remark that the functions that are added to the space X_N in step 5 are derived from ${\text{U}}_{N,h}^{n^{*}}$ and not the truth solution ${\tilde{\text{U}}}_h^{n^{*}}$. We have no guarantee that ${\text{U}}_{N,h}^{n^{*}}$ is close to ${\tilde{\text{U}}}_h^{n^{*}}$ or that it belongs to the low dimensional manifold ${{M}}_h^{\mu }$ of the truth solutions. We would like also to remark that even if we are trying to reduce the error ${\tilde{\text{U}}}_N^{n+1}-{\text{U}}_{N,h}^{n+1}$, to date we have no proof that the algorithm converges. In fact, we cannot theoretically prove that

$\begin{array}{ll}\Vert |{\tilde{U}}_N^{n+1}-{\tilde{U}}_h^{n+1}\Vert {|}_{*}\le \Vert |{\tilde{U}}_{N-1}^{n+1}-{\tilde{U}}_h^{n+1}\Vert {|}_{*}.& (19)\end{array}$

For this lack of theoretical convergence results, to stop the greedy enrichment procedure, we rather opt for a fixed number of solutions N_max chosen a priori, instead of using a certain tolerance on the a posteriori error estimator. Nevertheless, in the next section we will show some numerical evidence that the greedy enrichment is able to improve the quality of the approximation space.

4. Application to a Femoropopliteal Bypass

4.1. Application and Motivation

Atherosclerotic plaques often occur in the femoral arteries. The obstruction of the blood flow results in a lower perfusion of the lower limbs and the most common symptom of this disease is an intermittent claudication, which affects the 4% of people over the age of 55 years [34]. In order to restore the physiological blood circulation, different medical treatments are possible. In critical cases, the stenosis is treated with surgical intervention that bypasses the obstruction using a graft and providing an alternative way where blood can flow. The bypass creates a side-to-end anastomosis between the graft and the upstream artery (before the occlusion) and an end-to-side anastomosis with the distal downstream part. In particular, the design of end-to-side anastomosis affects the flow downstream the bypass and provokes remodeling phenomena inside the arterial wall. The arteries adapt their size in order to maintain a certain level of shear stress, which results in a thickening of the intima layer and in an increasing risk of thrombi formation. The arterial wall remodeling is in fact linked with hemodynamic factors such as the wall shear stress magnitude and direction. Moreover, velocity profiles and separation of flows have been investigated when studying the bypass end-to-side anastomosis [35, 36]. Studies with idealized geometrical models have been proposed in order to define an optimal design for the anastomosis [37]. Nevertheless, the geometry of the vessel is one of the most important factors that affect the pattern of the wall shear stress. Further, patient specific data would be required in order to analyse each particular case.

We focus our attention on the patient-specific femoropopliteal bypass performed with a venous graft bridging the circulation from the femoral artery to the popliteal one. As a domain of interest we select the end-to-side anastomosis (see Figure 1). The geometry was reconstructed by CT-scan images as it is detailed in Marchandise et al. [38] and inlet and outlet flow rates are provided from the experimental data in Marchandise et al. [38]. We compute the Reynolds number as $Re=\frac{4{\rho }_f{Q}_{in}}{\pi D\mu }$, being ρ_f the blood density, Q_in the inlet flow rate, D the vessel diameter and μ the blood viscosity. The average Reynolds number ranges from 144 and 380 (at the systolic peak), in agreement with the values provided in Loth et al. [36].

FIGURE 1

Figure 1. Realistic geometry of the end-to-side anastomosis. Graph of the inlet and outlet boundary conditions. (A) End-to-side geometry. (B) Inlet flow rate. (C) Outlet pressure.

4.2. Test Case

4.2.1. Application of the POD Algorithm

In this section we investigate the behavior of the POD and the greedy enriched POD algorithms on a case representing the femoropopliteal bypass application where the finite element resolution is performed on a coarse mesh. The usage of a coarse grid allows us to lower the offline computational costs and, thus, to test and compare several reduced basis approximations. Since we are interested in the realistic application of the femoropopliteal bypass, the physical parameters and boundary data are patient-specific. The coarse mesh is composed by 5,823 tetrahedra and 1.309 vertices. To obtain the high fidelity solutions of the RFSI model we use standard ℙ¹+Bubble-ℙ¹ finite elements for a total of 22,702 degrees of freedom. The boundary conditions are periodic with period of 0.8 s (one heartbeat). We set the solutions at time t = 0 equal to zero. To get rid of the dependence of these initial condition we perform the simulation of an entire heartbeat and we focus on the solutions obtained for the subsequent heartbeat. Thus, to test the POD reduction algorithm, we compute the high fidelity numerical solutions for a time lapse corresponding to the second heartbeat, from t₀ = 0.8 s to t_NT = 1.6 s with a time step Δt = 0.001 for a number of time intervals N_T = 800. We denote with the superscript $n\in {{N}}_T$ varying from 0 to N_T the sequence of computed solutions:

$\begin{array}{l}{U}_h^n\approx U_h(t_n)\text{ where }{t}_0=0.8,t_1=t_0+\Delta t,t_2=t_0+2\Delta t,\\ t_3=t_0+3\Delta t,\mathrm{..},t_{N_T}=1.6\text{s}.\end{array}$

We save the finite element solutions every five time steps and we use the apex $n_S\in {{N}}_S$, n_S = 5k, with k = 0, .., N_S (N_S = 160) to address the stored functions, that represent the snapshot sample:

$\begin{array}{l}{U}_h^{n_S}\approx U_h(t_{n_S})\text{ where }{t}_5=0.805,t_{10}=0.810,t_{15}=0.815,\\ t_{20}=0.820,\mathrm{..},t_{n_{N_S}}=1.6\text{s}.\end{array}$

Indeed, we compute the POD starting from the 160 snapshots ${\text{U}}_h^{n_S}$, $n_S\in {{N}}_S$, which represent the 25% of the finite element solutions computed for the second heartbeat. To check the quality of the reduced space approximations, we monitor the following errors:

• relative error of the velocity at time t_nS and correspondent space-time error:

$\begin{array}{l}{\epsilon }_N(u^{n_S}):\text{}\text{}=\frac{\Vert u_N^{n_S}-u_h^{n_S}{\Vert }_V}{\Vert u_h^{n_S}{\Vert }_V}\text{and}\\ E_N(u):\text{}\text{}=\frac{{\left({\displaystyle {\sum }_{n_S\in {{N}}_S}{(\Vert u_N^{n_S}-u_h^{n_S}{\Vert }_V)}^2}\right)}^{1/2}}{{\left({\displaystyle {\sum }_{n_S\in {{N}}_S}{(\Vert u_h^{n_S}{\Vert }_V)}^2}\right)}^{1/2}};& (20)\end{array}$

• relative error of the pressure at time t_nS and correspondent space-time error:

$\begin{array}{l}{\epsilon }_N(p^{n_S}):\text{}\text{}=\frac{\Vert p_N^{n_S}-p_h^{n_S}{\Vert }_{L^2(\Omega )}}{\Vert p_h^{n_S}{\Vert }_{L^2(\Omega )}}\text{and}\\ E_N(p):\text{}\text{}=\frac{{\left({{\displaystyle {\sum }_{n_S\in {{N}}_S}(\Vert p_N^{n_S}-p_h^{n_S}{\Vert }_{L^2(\Omega )})}}^2\right)}^{1/2}}{{\left({\displaystyle {\sum }_{n_S\in {{N}}_S}{(\Vert p_h^{n_S}{\Vert }_{L^2(\Omega )})}^2}\right)}^{1/2}};& (21)\end{array}$

• space-time dual norm of the residual scaled with respect to the space-time norm of the global solution

$\begin{array}{ll}{R}_N(U):\text{}\text{}={\left(\frac{N_S}{N_T}\right)}^{1/2}\frac{{\left({\displaystyle {\sum }_{n\in {{N}}_T}\Vert r_N^n(W_h){\Vert }_{X^{\prime }}^2}\right)}^{1/2}}{{\left({\displaystyle {\sum }_{n_S\in {{N}}_S}\Vert U_h^{n_S}{\Vert }_X{)}^2}\right)}^{1/2}};& (22)\end{array}$

We build a sequence of POD reduced spaces with decreasing values of the tolerance tol and we compute the aforementioned indicators for each one of the reduced spaces generated. The space-time errors are reported in Table 1. In particular, we show: the number of selected velocity modes (#u basis); the number of selected pressure modes (#p basis); the total number of basis functions composing the reduced space (# basis = #u basis + 2 × #p basis); the space-time errors and residuals as defined above. Since the problem at hand is unsteady and the solution at a time instant t_n depends on the solutions at the previous instants, the POD model errors E_N(u) and E_N(p) are bounded from above by the fixed tolerance but they are however of the same order of magnitude (see Table 1). We notice that, even if $||r_N^n({\text{W}}_h)|{|}_{{\text{X}}^{\prime }}$ does not represent an upper bound for the error, nevertheless, from experimental results, we can use it as an indicator of $||{\text{U}}_N^{n_S}-{\text{U}}_h^{n_S}|{|}_{\text{X}}$ (see Figure 2). The apparent strong correlation between the dual norm of the residual and the global error norm is probably due to the strong contribution of the mass term in the unsteady formulation. Indeed, if we choose a time step of 0.001, the mass matrix is multiplied for a factor of 10³. We remark that the magnitudes of the absolute errors for the velocity span from 10⁻¹ to 10² and the associated velocity solutions norms are of order of 10² − 10³. For the pressure, we have absolute errors of order 10⁰ − 10³, while their solutions norms are of order 10³⁻⁵. Indeed, in absolute terms the global error is mostly related with the pressure one.

TABLE 1

Table 1. Number of basis functions and space-time errors for the velocity and pressure.

FIGURE 2

Figure 2. Dual norms of the residuals and norms of the global errors with respect to time for different choices of the POD tolerance tol. (A) tol = 1e-5. (B) tol = 1e-6. (C) tol = 1e-7.

4.2.2. Application of the Greedy Enriched POD Algorithm

The POD algorithm provided satisfactory results and we were able to reduce the approximation space dimension from 10⁵ to 10 − 100. In this section we aim at comparing the greedy enrichment algorithm with the POD one, in order to understand if using different basis functions than POD modes provides the same quality of reduced approximations. Thus, we compare the magnitude of the reduced approximation errors obtained using a reduced space generated through a standard POD algorithm with the ones obtained using the POD coupled with the greedy enrichment as introduced in section 3.2. We recall that the snapshots sample represents a subset of the time instants we solve in the unsteady simulation: indeed we store only the 25% of the time instants solutions computed. As there is no error bound available, we use the dual norm of the residual as surrogate. This is a rough approximation, also because the real error includes time integration, while the dual norm of the residual can only represent a space error. Of course we do not expect the greedy enriched POD to perform better, on the contrary it can have (and actually has) limitations.

Remark We are interested to simulate a fluid-dynamics phenomena with cyclic inputs. Typically in hemodynamics applications, we are interested in several heartbeats. Thus, instead of performing the greedy research only on one single heartbeat, we exploit as much as we can the information on the truth solutions coming from the snapshots. For each single snapshot ${\tilde{\text{U}}}_h^{n_s}$, we perform a simulation that starts from the initial time t_ns and ends at t_ns+NT = t_ns+0.8. We define a vector index n = (n_T, n_S) with n_T = n_S+n such that ${\tilde{\text{U}}}_h^{\text{n}}={\tilde{\text{U}}}_h^{n_S,n_T}$ being the approximate solution at time $t^{n_T}$ obtained starting from the initial condition ${\tilde{\text{U}}}_h^{n_s}$. We define the set of indexes ${N}=\left\{(n_T,n_S):n_T=n_S+n,n\in {N}\mathrm{\text{ and }}{n}_S\in {{N}}_S\right\}$. The generalization of the greedy enrichment presented in Section 3.2 is straightforwards substituting n with n. In particular, the selection of the worst approximated index n* in the greedy enrichment can be generalized as follows:

$n^{*}=\arg \underset{n\in {N}}{\max }\Vert r_N^n(W_h){\Vert }_{X^{\prime }}.$

In order to initialize the greedy enrichment algorithm we compute a POD basis fixing the tolerance tol = 1e − 5 (50 velocity modes, 3 pressure modes, 3 supremizers). To compare the POD approximation with the greedy enriched one, we augment the initial reduced space with two strategies. On one side, we apply the greedy enrichment and, at each iteration, we add the triplet of functions selected by the largest dual norm of the residual in space. On the other side, we augment the basis by adding, at each algorithm iteration, one POD mode for the velocity and one POD mode for the pressure with its associated supremizer. In both cases, at each iteration, we increase the reduced space dimensions of three units. The results obtained using only POD modes are displayed in black and addressed with the label POD, while the results obtained with the greedy enrichment are shown in red and addressed with the label Greedy enriched POD (see Figure 3).

FIGURE 3

Figure 3. Space-time errors vs. the number of basis, comparison between the standard POD and the greedy enriched POD, starting from 50 POD basis functions. FEM solutions obtained on a coarse mesh, bypass application. (A) Space-time velocity error. (B) Space-time pressure error. (C) Space-time dual norm of the residuals.

From Figure 3, we note that the decrements of the errors in the greedy enrichment algorithm are slower than when adding POD modes. Nevertheless, we notice that both the space-time pressure error and residuals are comparable when adding POD modes or greedy basis functions (see Figures 3B,C). On the contrary the decrements of the velocity is much slower when we use the greedy enrichment with respect to adding POD modes. We recall, however, than the residual is mostly related to the pressure error component.

4.3. Realistic Case

4.3.1. Application of the POD Algorithm

In this section we perform a discretization reduction of the RFSI model applied to the femoropopliteal bypass case, where the high fidelity approximations are computed using a fine mesh. As before, a parabolic velocity profile is imposed at the inlet section and a mean pressure condition at the outlet. The ℙ¹+Bubble-ℙ¹ discretization yields 1,410,475 degrees of freedom on the fine mesh. We first test the discretization reduction using a standard POD procedure: we compute the high fidelity numerical solutions for two heartbeats with a time step Δt = 0.001 and we store the ones related to the second heartbeat every five time steps. Thus, ${{N}}_T=1,2,3,\dots ,800$ and ${{N}}_S=5,10,15,..,800$. We compute the Gramian matrices associated to the 160 snapshots ${\text{U}}_h^{n_k}$, separating the velocity and pressure components. We denote ${\lambda }_k^{\text{u}}$ and ${\lambda }_k^p$ for k = 1, ..N_S the eigenvalues associated to the decomposition of the correlation matrices of the velocity and pressure, respectively (see section 3.1). In both cases they decrease exponentially fast. The eigenvalues ${\lambda }_k^p$ associated to the pressure snapshots (Figure 4B) decrease faster than the ones associated to the velocity (Figure 4A). Thus, by fixing the same tolerance, we expect that a fewer number of pressure modes will be selected with respect to the velocity ones.

FIGURE 4

Figure 4. Velocity and pressure eigenvalues computed with 160 snapshots sampled every 0.005 s. (A) Velocity eigenvalues, ${\lambda }_k^{\text{u}}$. (B) Pressure eigenvalues, ${\lambda }_k^p$.

We compute the POD reduced spaces using three different values for the tolerance: tol ∈ {1e − 4, 1e − 5, 1e − 6}. As it was done in section 4.2, in Table 2, we record the selected number of modes and we compute the space-time errors E_N(u) and E_N(p) of the velocity and pressure, respectively. By taking advantage of the generated reduced space, at each time iteration we solve the reduced system and we compute a linear functional of the approximate solution that evaluates the outlet flow rate. We record the computational time associated with the offline and online computations in Table 3. We can appreciate that the resolution of the reduced problem combined with the evaluation of a linear output functional is performed in almost real time: using 35 basis functions we solve 1.6 physical seconds in 0.8 computational seconds, while a ten heartbeats simulation (8 physical seconds) selecting the POD space with 54 basis functions takes 12.6 s on a notebook. Performing the same simulation with the high fidelity model would have taken around 40 h on 256 processors of a supercomputer. The offline costs of the POD reduction (without the snapshots generation) are reported in the third column of Table 3. We remark that most of this time is spent in the generation of the structures for the residual evaluation.

TABLE 2

Table 2. Number of basis functions and space-time errors for the velocity and pressure.

TABLE 3

Table 3. CPU time $X_N^{POD}$: offline computations costs for the generation of the POD reduced spaces (without the finite element computations) on 512 processors; CPU time 2HB - RB: online computational time corresponding to the simulation of 2 heartbeats (2HB) on a personal laptop; CPU time 2HB - FE: finite element computational time corresponding to the simulation of 2 heartbeats (2HB) on 256 processes on a supercomputer.

Note that the POD model errors E_N(u) and E_N(p) decrease significantly when increasing the number of basis functions, as it is reported from both the values of Table 2. Once again, we notice that the dual norm of the residual $||r_N^{n+1}({\text{W}}_h)|{|}_{{\text{X}}^{\prime }}$ is a good indicator of the approximation error $||{\text{U}}_N^{n_k}-{\text{U}}_h^{n_k}|{|}_{\text{X}}$ (see Figure 5).

FIGURE 5

Figure 5. Dual norm of the residuals and X-norm of the errors with respect to the time instant for different choices of the POD tolerance tol. FEM solutions obtained on a fine mesh, bypass application. (A) tol = 1e-4. (B) tol = 1e-5. (C) tol = 1e-6.

In the femoropopliteal bypass application, we are interested in measuring also the errors on the output of interests. Being ${\sigma }^{n_S}$ the stress tensor and n the normal vector to the surface Γ, we compute the wall shear stress as ${\tau }^{n_S}:={\sigma }^{n_S}n-({\sigma }^{n_S}n·n)n$ and we also consider the averaged wall shear stresses on a generic area A: ${\tau }_A^{n_S}=1/A\underset{A}{\int }\left|{\tau }^{n_S}\right|dA$. We remark that to properly estimate the selected output of interest we need accurate high fidelity solutions with a mesh refined at the wall, as shown by Marchandise et al. [38]; the fine mesh used in this work is similar to the fine one used in that paper. Reducing the dimension of the finite element space does not lead to the same results that we obtain reducing the degrees of freedom using the POD decomposition. In fact, the wall shear stress values computed using a coarse finite element space underestimate considerably the values obtained with the fine grid (see Figure 6), while the results obtained with the POD reduced approximation mostly overlapped with the ones computed with the finite element discretization.

FIGURE 6

Figure 6. Averaged wall shear stress (in [g/cm s]) computed in different location of the interface. Comparisons between the values obtained with the finite element discretization on the fine mesh (FEM Fine), the coarse one (FEM Coarse) and the reduced basis approximation (POD tol = 1e-5). (A) Areas location. (B) Wss absolute values - Area 1. (C) Wss absolute values - Area 2.

4.3.2. Percentage of Flow Coming From the Occluded Artery

In many cases, the artery is not completely occluded but a residual flow is still provided by the original vessel. Studying the distribution of fluid-dynamic quantities could be important to identify, for example, whether it would be better or not to surgically close the original artery. Our approach allows to study the variation of the flow and wall shear stress for different values of flow percentage coming from the occluded artery with low computational costs.

Here the percentage of flow coming from the occluded artery is an additional parameter. We solve the high fidelity model for two extreme cases: μ₁ = 0% (full occlusion) and μ₂ = 50% of flow. Using a POD algorithm with tol = 10⁻⁵, we construct the reduced space ${\text{X}}_{N_1}^{POD,{\mu }_1}$ of dimension N₁, associated with the snapshots computed with μ₁. We obtain N₁ = 54, where 48 basis functions are velocity modes, 3 pressure modes and 3 supremizers. Then, we consider the snapshots associated with μ₂, we build a second POD reduced space ${\text{X}}_{N_2}^{POD,{\mu }_2}$ of dimension N₂, also fixing tol = 10⁻⁵. We obtain N₂ = 55 where 47 basis functions are velocity modes, 4 pressure ones and 4 supremizers. Finally, we check that the basis functions are linear independent and orthonormalize the basis of ${\text{X}}_{N_2}^{POD,{\mu }_2}$ with respect to ${\text{X}}_{N_1}^{POD,{\mu }_1}$. Indeed, none of the basis functions obtained for μ₂ = 0.5 is as a linear combination of the basis related to μ₁ = 0 and the final reduced space X_N has dimension N = N₁ + N₂ = 109.

We then choose a parameter μ₃ = 25% of flow coming from the occluded artery. We focus on the velocity profiles near the systolic peak and in the early diastole (see Figure 7) and the averaged wall shear stress (see Figure 8). We compare the target outputs obtained using the finite element approximation (label: FEM 25%) and the reduced ones (label: POD 25%). Moreover, we display the selected quantities also for the high fidelity solutions obtained with μ₁ = 0 (label: FEM 0%), μ₂ = 0.5 (label: FEM 50%) in order to clarify how the system dynamics changes when different values of the percentage are considered.

FIGURE 7

Figure 7. Velocity profiles for two time instants of the diastolic phase. Comparisons between finite element and reduced approximations. Femoropopliteal bypass application in which the high fidelity solutions are obtained using a finite element approximation on a fine mesh. (A) Systole, 0.9 s. (B) Peak systole, 1 s. (C) Diastole, 1.1 s. (D) Diastole, 1.2 s.

FIGURE 8

Figure 8. Averaged wall shear stress computed in the location A₁, cf. Figure 6. Comparisons between finite element (corse and fine meshe) and reduced approximations.

During the systolic phase the reduced solutions well reproduce the high fidelity ones, while during the diastole, the differences are more visible, in particular in some locations as, for example, near the anastomosis between the arterial vessel and the bypass graft (see Figure 7). We remark that the velocity profiles for different values of the parameter are significantly different between themselves. Our reduced solution approximates well the value of the wall shear stress associated to μ₃ (see Figure 8). We can appreciate the good agreement between the reduced and the high fidelity results, even when the values of wall shear stress associated to μ₃ are not berween those associated to μ₁ and μ₂ (see Figure 8B). Moreover, the approximation of the wall shear stress using a finite element approximation with a coarse grid leads to a consistent underestimation of their values (see Figure 8C).

4.3.3. Application of the Greedy Enriched Algorithm With Perturbed Data

In this section we apply the greedy enrichment in the case of perturbed boundary data. In particular, as in (9), we introduce a parameter in the inlet flow rate function representing a small perturbation with respect to a reference value. The perturbation function θ(α, t) is define as follows:

$\theta (\alpha ,t)=1+\alpha \sin \left(\frac{2\pi t}{0.8}\right)$

where α is supposed to vary between 0 and 0.2. Thus, the maximum relative difference with the original flow rate is equal to the 20%. We denote with ${\tilde{\text{U}}}_{*}^{\text{n}}(\alpha )$ with * = {h}, {N} or {N, h} the numerical solutions at the time instant tⁿ that depend on the parameter α and with $r_N^{\text{n}}({\text{W}}_h;\alpha )$ the residuals. In the perturbed case, the algorithm steps in section 3.2 are modified as follows. First we perform a POD algorithm fixing α = α₁ = 0.0; the resulting reduced space is addressed with ${\text{X}}_N^{{\alpha }_1,POD}$, where the apex α₁ denotes the choice of the α parameter. Then, we set α₂ = 0.2:

1. Generate the reduced basis solutions ${\tilde{\text{U}}}_N^{\text{n}}({\alpha }_2)$, $\text{n}\in {N}$ by solving the reduced order problem (10).

2. Compute the dual norms of the residuals $||r_N^{\text{n}}({\text{W}}_h;{\alpha }_2)|{|}_{{\text{X}}^{\prime }}$, $\text{n}\in {N}$, which are used as error indicators.

3. Select n* such that ${\text{n}}^{*}=arg\underset{\text{n}\in {N}}{max}||r_N^{\text{n}}({\text{W}}_h;{\alpha }_2)|{|}_{{\text{X}}^{\prime }}$.

4. Compute the ${\text{U}}_{N,h}^{{\text{n}}^{*}}({\alpha }_2)$ by solving the reduced order problem (17).

5. Split ${\text{U}}_{N,h}^{{\text{n}}^{*}}({\alpha }_2)$ into its velocity and pressure components, ${\text{u}}_h^{{\text{n}}^{*}}$ and $p_h^{{\text{n}}^{*}}$, respectively. Compute the supremizer ${\sigma }_{{\text{n}}^{*}}$ associated with the pressure component.

6. Compute ϕ^u representing the orthonormalization of the velocity function ${\sigma }^{n*}$ with respect to the reduced space X_N, obtained with a Gram-Schmidt algorithm considering the scalar product (·, ·)_V; similarly for ϕ^p and $p_h^{{\text{n}}^{*}}$.

7. Build ${\text{X}}_{N+2}={\text{X}}_N\oplus \left\{{\psi }^{\text{u}},{\psi }^p\right\}$ built as is (16).

8. Compute ϕ^σ representing the orthonormalization of the velocity function ${\text{u}}_h^{n^{*}}$ with respect to the reduced space X_N+2, obtained with a Gram-Schmidt algorithm considering the scalar product (·, ·)_V.

9. Build ${\text{X}}_{N+3}={\text{X}}_{N+2}\oplus \left\{{\psi }^{\sigma }\right\}$ built as is (16).

10. Update the structures for the online computation of the reduced solutions and the dual norms of the residuals.

11. Set N = N + 3 and X_N = X_N+3. Repeat until a predefined stopping criterion is satisfied.

The real modification is indeed related to the fact that the initial POD is computed for α = α₁ = 0, while the greedy enrichment is performed fixing α = α₂ = 0.2. The resulting reduced space X_N aims to represent a suitable space of approximation for both values of α. In the parametrized case, by using greedy enrichment we aim at saving a part of the offline computational costs: indeed, in a standard POD-Greedy procedure (see [39]), each new evaluation of the parameter α requires the computation of the associated finite element solutions for each time instant $n\in {{N}}_T$. In our application, this would require about 8 h on 256 processors. Instead, during the greedy enrichment, we perform only one finite element resolution for a single time step, while the remaining computations are dedicated to reduced basis structures.

We test the greedy enrichment algorithm by initializing it with two different starting POD reduced spaces: in one case we consider the modes selected with tol = 1e − 4 (35 POD basis functions) and in the other one we consider the POD modes corresponding to tol = 1e − 5 (54 POD basis functions). In the first case, we enrich the space ${\text{X}}_{35}^{{\alpha }_1,POD}$ by adding 8 triplets selected by the greedy algorithm; we obtain the reduced space X₅₉. In the second case, starting from ${\text{X}}_{54}^{{\alpha }_1,POD}$, we enrich the space adding 12 triplets, obtaining X₉₀. All the errors and residuals computed and shown below are referred to the solutions obtained with α₂ = 0.2. In particular, in Table 4 we report the velocity and pressure errors generated by the greedy enriched reduced spaces as well as the ones obtained with the standard POD ones. Moreover, we compute the space-time dual norm of the residual, scaled by the solution norm (sixth column of Table 4).

TABLE 4

Table 4. Number of basis functions and space-time errors for the velocity and pressure.

We note that the space-time velocity error does not decrease significantly neither when adding greedy basis functions nor when augmenting the number of selected POD modes. If we look at the pressure, using the greedy enrichment we manage to decrease its error more than if we use POD modes. Also the space-time dual norm of the residual is smaller when considering the greedy enriched space than the POD ones.

Regarding the offline costs, to generate the space X₅₉ starting form the ${\text{X}}_{35}^{0,POD}$, we perform 8 iterations of the greedy enrichment algorithm: this takes 82 min on 512 processors where the most of the time is devoted to the generation of the reduced structures for the residual evaluation. We remark that computing a standard POD reduced space for the parameter evaluation corresponding to α = 0.2 would require about 8 h on 256 processors for the finite element computation of two periods, plus about 1 h on 512 processors for the generation of the reduced space itself.

To explain why we obtain better results for the pressure than for the velocity, we investigate the absolute values of velocity, pressure and global solutions errors and we compare them to the dual norms of the residuals (see Figure 9). Since the velocity and pressure norms have two different magnitudes (10 − 10² for the velocity and 10³ − 10⁵ for the pressure), the corresponding absolute values of the pressure errors are bigger than the velocity ones, even if the relative errors are lower. The greedy procedure selects the worst approximated time instant based on the dual norms of the residuals and these quantities are indicators of the global absolute errors. Since the latter is mostly due to the pressure error, this can explain why the greedy enrichment provides better results for the pressure than for the velocity.

FIGURE 9

Figure 9. Dual norms of the residuals and norms of the global errors with respect to time for different choices of the POD tolerance tol. Femoropopliteal bypass application in which the high fidelity solutions are obtained using a finite element approximation on a fine mesh. (A) 54 Basis - POD. (B) 78 Basis - POD. (C) 90 Basis - Greedy enriched POD.

5. Conclusions

In this work we presented an application of reduced order modeling to a RFSI problem that is indeed an unsteady Navier-Stokes problem with generalized Robin boundary conditions. We detailed how an affine decomposition with respect to boundary data varying in time can be obtained under suitable hypothesis. Moreover, we presented and detailed how the POD can be applied to the RFSI problem in order to take into account the different order of magnitudes of the variables. We discussed the introduction of the supremizer functions inside the reduced basis, necessary to include the pressure in the reduced system. Afterwards, we proposed an enrichment of the POD reduced basis based on a greedy algorithm. All the algorithms presented were then numerically tested on a realistic hemodynamics problem. We tested the POD and greedy enrichment algorithm on two cases: a test case, where the finite element solution is obtained with a coarse grid, and a fine case, where the finite element space has order of 10⁶ degrees of freedom. The results showed the good performances of the POD reduction algorithm on the RFSI problem, also with respect to the evaluation of specific hemodynamics target output (wall shear stress). Moreover we provided numerical evidence of how the reduced approximation can be improved using the greedy enrichment algorithm, in particular regarding the pressure error. The different behavior of the velocity and pressure errors is due to the use of the dual norm of the residual as an indicator of the global solution error. Indeed, since we do not have suitable a-posteriori error estimators, one for the velocity and one for the pressure variables, we measure the dual norm of the residual as a surrogate estimator. Being the pressure variable and the correspondent error two order of magnitudes grater that the velocity ones, the residual is indeed an indicator of the pressure errors. Nevertheless, even in lack of theoretical results, numerical experiments showed that the greedy enrichment is able to improve the quality of the reduced approximation allowing us to save computational time. The development of suitable a-posteriori error estimators for the pressure and velocity in the case of RFSI problem would be required to improve the performances of the greedy enrichment.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer, JC, and handling Editor declared their shared affiliation.

Acknowledgments

The authors acknowledge support of the European Research Council under the Advanced Grant ERC-2008-AdG 227058, Mathcard, Mathematical Modeling and Simulation of the Cardiovascular System, and of the FNS project 200021E-168311, Domain-Decomposition-Based Fluid Structure Interaction Algorithms for Highly Nonlinear and Anisotropic Elastic Arterial Wall Models in 3D. Moreover, we thank Prof. Bernard Haasdonk and Prof. Alfio Quarteroni for the insightful discussion. We gratefully acknowledge the Swiss National Supercomputing Center (CSCS) for providing the CPU resources for the numerical simulations under the projects IDs s475 and s796.

References

1. Colciago CM, Deparis S, Quarteroni A. Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics. J Comput Appl Math. (2014) 265:120–38. doi: 10.1016/j.cam.2013.09.049

CrossRef Full Text | Google Scholar

2. Bazilevs Y, Michler C, Calo VM, Hughes TJR. Weak Dirichlet boundary conditions for wall-bounded turbulent flows. Comput Methods Appl Mech Eng. (2007) 196:4853–62. doi: 10.1016/j.cma.2007.06.026

CrossRef Full Text | Google Scholar

3. Manzoni A, Negri F. Rigorous and Heuristic Strategies for the Approximation of Stability Factors in Nonlinear Parametrized PDEs. Preprint available as Mathicse Report Nr. 08.2014 (2014).

4. Lassila T, Manzoni A, Quarteroni A, Rozza G. Model order reduction in fluid dynamics: challenges and perspectives. In: Quarteroni A, Rozza G, editors. Reduced Order Methods for Modeling and Computational Reduction. Vol. 9 of MS&A Series. Milano: Springer (2014).

5. Rozza G, Huynh DBP, Patera AT. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch Comput Methods Eng. (2008) 15:229–75. doi: 10.1007/s11831-008-9019-9

CrossRef Full Text | Google Scholar

6. Quarteroni A, Rozza G, editors. Reduced Order Methods for Modeling and Computational Reduction. MS&A, Modeling, Simulation and Applications. Milano: Springer (2013). Available online at: http://www.springer.com/mathematics/computational+science+

Google Scholar

7. Yano M, Patera AT, Urban K. A space-time hp-interpolation-based certified reduced basis method for Burgers' equation. Math Models Methods Appl Sci. (2014) 24:1903–35. doi: 10.1142/S0218202514500110

CrossRef Full Text | Google Scholar

8. Nguyen NC, Rozza G, Patera AT. Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers equation. Calcolo (2009) 46:157–85. doi: 10.1007/s10092-009-0005-x

CrossRef Full Text | Google Scholar

9. Wirtz D, Sorensen DC, Haasdonk B. A posteriori error estimation for DEIM reduced nonlinear dynamical systems. SIAM J Sci Comput. (2014) 36:A311–38. doi: 10.1137/120899042

CrossRef Full Text | Google Scholar

10. Deparis S, Lovgren AE. Stabilized reduced basis approximation of incompressible three-dimensional Navier-Stokes equations in parametrized deformed domains. J Sci Comput. (2012) 50:198–212. doi: 10.1007/s10915-011-9478-2

CrossRef Full Text | Google Scholar

11. Iapichino L, Quarteroni A, Rozza G. A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks. Comput Methods Appl Mech Eng. (2012) 221–222:63–82. doi: 10.1016/j.cma.2012.02.005

CrossRef Full Text | Google Scholar

12. Amsallem D, Farhat C. An online method for interpolating linear parametric reduced-order models. SIAM J Sci Comput. (2011) 33:2169–98. doi: 10.1137/100813051

CrossRef Full Text | Google Scholar

13. Gerner AL, Veroy K. Reduced basis a posteriori error bounds for the Stokes equations in parametrized domains: a penalty approach. Math Models Methods Appl Sci. (2011) 21:2103–34. doi: 10.1142/S0218202511005672

CrossRef Full Text | Google Scholar

14. Deparis S, Rozza G. Reduced basis method for multi-parameter-dependent steady NavierStokes equations: applications to natural convection in a cavity. J Comput Phys. (2009) 228:4359–78. doi: 10.1016/j.jcp.2009.03.008

CrossRef Full Text | Google Scholar

15. Grinberg L, Yakhot A, Karniadakis G. Analyzing transient turbulence in a stenosed carotid artery by proper orthogonal decomposition. Ann Biomed Eng. (2009) 37:2200–17. doi: 10.1007/s10439-009-9769-z

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Deparis S. Reduced basis error bound computation of parameter-dependent Navier-Stokes equations by the natural norm approach. SIAM J Numer Anal. (2008) 46:2039–67. doi: 10.1137/060674181

CrossRef Full Text | Google Scholar

17. Gunzburger MD, Peterson JS, Shadid JN. Reduced-order modeling of time-dependent {PDEs} with multiple parameters in the boundary data. Comput Methods Appl Mech Eng. (2007) 196:1030–47. doi: 10.1016/j.cma.2006.08.004

CrossRef Full Text | Google Scholar

18. Sen S, Veroy K, Huynh DBP, Deparis S, Nguyen NC, Patera AT. Natural norm a posteriori error estimators for reduced basis approximations. J Comput Phys. (2006) 217:37–62. doi: 10.1016/j.jcp.2006.02.012

CrossRef Full Text | Google Scholar

19. Nichols WW, O'Rourke MF. McDonald's Blood Flow in Arteries - Theoretical, Experimental and Clinical Principles. London: Arnold (1998).

Google Scholar

20. Formaggia L, Quarteroni A, Veneziani A. Cardiovascular Mathematics. Heidelberg: Springer (2009).

Google Scholar

21. Figueroa CA, Vignon-Clementel IE, Jansen KE, Hughes TJR, Taylor CA. A coupled momentum method for modeling blood flow in three-dimensional deformable arteries. Comput Methods Appl Mech Eng. (2006) 195:5685–706. doi: 10.1016/j.cma.2005.11.011

CrossRef Full Text | Google Scholar

22. Moireau P, Xiao N, Astorino M, Figueroa CA, Chapelle D, Taylor CA, et al. External tissue support and fluid-structure simulation in blood flows. Biomech Model Mechanobiol. (2012) 11:1–18. doi: 10.1007/s10237-011-0289-z

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Figueroa CA, Baek S, Taylor CA, Humphrey JD. A computational framework for fluid-solid-growth modeling in cardiovascular simulations. Comput Methods Appl Mech Eng. (2009) 198:3583–601. doi: 10.1016/j.cma.2008.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Nobile F, Vergara C. An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions. SIAM J Sci Comput. (2008) 30:731–63. doi: 10.1137/060678439

CrossRef Full Text | Google Scholar

25. Kashiwabara T, Colciago CM, Dedé L, Quarteroni A. Well-posedness, regularity, and convergence analysis of the finite element approximation of a generalized Robin boundary value problem. SIAM J Numer Anal. (2015) 53:105–26. doi: 10.1137/140954477

CrossRef Full Text | Google Scholar

26. Sirisup S, Karniadakis GE. Stability and accuracy of periodic flow solutions obtained by a POD-penalty method. Physica D (2005) 202:218–37. doi: 10.1016/j.physd.2005.02.006

CrossRef Full Text | Google Scholar

27. Hesthaven JS, Rozza G, Stamm B. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. Cham: Springer (2016).

Google Scholar

28. Quarteroni A, Manzoni A, Negri F. Reduced Basis Methods for Partial Differential Equations: An Introduction, Vol. 92. Heildelberg: Springer (2015).

Google Scholar

29. Rowley CW. Model reduction for fluids, using balanced proper orthogonal decomposition. Int J Bifurcat Chaos. (2005) 15:997–1013. doi: 10.1142/S0218127405012429

CrossRef Full Text | Google Scholar

30. Willcox K, Peraire J. Balanced model reduction via the proper orthogonal decomposition. AIAA J. (2002) 40:2323–30. doi: 10.2514/2.1570

CrossRef Full Text | Google Scholar

31. Gerner AL, Reusken A, Veroy K. Reduced Basis a Posteriori Error Bounds for the Instationary Stokes Equations. Singapore: World Scientific (2012).

Google Scholar

32. Burkardt J, Gunzburger M, Lee HC. {POD} and CVT-based reduced-order modeling of NavierStokes flows. Comput Methods Appl Mech Eng. (2006) 196:337–55. doi: 10.1016/j.cma.2006.04.004

CrossRef Full Text | Google Scholar

33. Rozza G, Veroy K. On the stability of the reduced basis method for Stokes equations in parametrized domains. Comput Methods Appl Mech Eng. (2007) 196:1244–60. doi: 10.1016/j.cma.2006.09.005

CrossRef Full Text | Google Scholar

34. Giordana S, Sherwin SJ, Peir J, Doorly DJ, Crane JS, Lee KE, et al. Local and global geometric influence on steady flow in distal anastomoses of peripheral bypass grafts. J Biomech Eng. (2005) 127:1087–98. doi: 10.1115/1.2073507

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Loth F, Jones SA, Zarins CK, Giddens DP, Nassar RF, Glagov S, et al. Relative contribution of wall shear stress and injury in experimental intimal thickening at PTFE end-to-side arterial anastomoses. J Biomech Eng. (2002) 124:44–51. doi: 10.1115/1.1428554

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Loth F, Fischer PF, Bassiouny HS. Blood flow in end-to-side anastomoses. Annu Rev Fluid Mech. (2008) 40:367–93. doi: 10.1146/annurev.fluid.40.111406.102119

CrossRef Full Text | Google Scholar

37. Lassila T, Manzoni A, Quarteroni A, Rozza G. Boundary control and shape optimization for the robust design of bypass anastomoses under uncertainty. ESAIM Math Model Numer Anal. (2013) 47:1107–31. doi: 10.1051/m2an/2012059

CrossRef Full Text | Google Scholar

38. Marchandise E, Crosetto P, Geuzaine C, Remacle JF, Sauvage E. Quality open source mesh generation for cardiovascula flow simulation. In: Ambrosi D, Quarteroni A, Rozza G, editors. Modeling of Physiological Flows. Milano: Springer (2011) 395–414.

39. Haasdonk B, Ohlberger M. Reduced basis method for finite volume approximations of parametrized linear evolution equations. ESAIM Math Model Numer Anal. (2008) 42:277–302. doi: 10.1051/m2an:2008001

CrossRef Full Text | Google Scholar