A non-intrusive data-driven ROM framework for hemodynamics problems

Reduced order modeling (ROM) techniques are numerical methods that approximate the solution of parametric partial differential equation (PDE) by properly combining the high-fidelity solutions of the problem obtained for several configurations, i.e. for several properly chosen values of the physical/geometrical parameters characterizing the problem. In this contribution, we propose an efficient non-intrusive data-driven framework involving ROM techniques in computational fluid dynamics (CFD) for hemodynamics applications. By starting from a database of high-fidelity solutions related to a certain values of the parameters, we apply the proper orthogonal decomposition with interpolation (PODI) and then reconstruct the variables of interest for new values of the parameters, i.e. different values from the ones included in the database. Furthermore, we present a preliminary web application through which one can run the ROM with a very user-friendly approach, without the need of having expertise in the numerical analysis and scientific computing field. The case study we have chosen to test the efficiency of our algorithm is represented by the aortic blood flow pattern in presence of a Left Ventricular Assist Device (LVAD) when varying the pump flow rate.


Introduction
Reduced order modeling (ROM) (see, e.g., [8]) is a well-spread technique used both in academia and in industry. It has been introduced as an efficient tool to approximate full order systems by significantly reducing the computational cost required to obtain numerical solutions in a parametric setting. ROM consists in two main stages: an offline phase that can be carried out on high performance computing facilities, and an online one that hinges on a system of reduced dimensionality to perform the parametric computation on portable devices. In the offline phase, the reduced order space is built starting from full order complex simulations computed for certain values of the physical and/or geometrical parameters. In this work, we employ the proper orthogonal decomposition (POD) for the detection of the reduced basis functions that span this new reduced space. After the creation of such a space, in the online phase a new parametric solution is obtained as a linear combination of the precomputed reduced basis functions, by means of an interpolation carried out by using RBF functions [22]. The resulting ROM is thus called proper orthogonal decomposition with interpolation (PODI) [9].
The aim of this work is the development of an efficient non-intrusive data-driven reduced order model to be used within hemodynamics framework. The reader can find examples of the ROM application in the hemodynamics field in [4,5,30,34,15]. We highlight that the online evaluation of the data-driven approach used here is based only on data and does not require knowledge about the governing equations that describe the system. It is also non-intrusive, i.e. no modification of the simulation software is carried out. For this reason it is particularly versatile thanks to its capability to be coupled with commercial solvers as well. It should be noted that many efforts are making in order to integrate ROM and technological innovation. From this viewpoint, a crucial step is the web server ARGOS [2], developed by mathLab group at SISSA that will make possible the exploitation of reduced order models to a wide category of people working in industrial and biomedical contexts. Through specific web applications the user will be able to solve many complex problems without the need of being an expert in numerical analysis and scientific computing. In particular, it is expected that the ATLAS project [3] will collect all cardiovascular applications. In this framework, we present a preliminary web application through which one can run the ROM by using a very user-friendly GUI interface. The benchmark we have chosen to test the efficiency of our algorithm is represented by the aortic blood flow pattern in presence of a Left Ventricular Assist Device (LVAD) (see, e.g., [1,29,21,11,10,14,6]) when varying the pump flow rate (see, e.g., [15,7,25]).
The work is organized as follows. In Sec. 2 we present the general parametric full order model governing hemodynamics problems, over which we apply the proposed numerical methodology. In Sec. 3 we present the PODI method, whilst in Sec. 4 we show the numerical setting of the problem and the achieved results, as well as provide a brief description of the web application developed. Finally, in Sec. 5 conclusions and perspectives are provided.

The full order model
In this work we consider the blood as modeled by the unsteady incompressible Navier-Stokes equations described in an Eulerian framework. We consider a fixed domain Ω ⊂ R d with d = 2, 3 over a time interval of interest (t 0 , T ) ⊂ R + . Let π ∈ P ⊂ R P be a parameter vector in a P -dimensional parameter space P. We have ρ ∂ t u(x, t; π) + ρ ∇ · (u(x, t; π) ⊗ u(x, t; π)) − 2µ∆u(x, t; π) + ∇p(x, t; π) = 0, in Ω × [t 0 , T ], endowed with proper boundary conditions. In (1)-(2), ∂ t denotes the time derivative, ρ = 1060 kg/m 3 is the blood density, µ = 0.004 Pa · s is the blood dynamic viscosity, u is the blood velocity and p is the pressure. We impose a no slip boundary condition on the wall of the domain. At the inflow, we prescribe a known flow-rate and ∂p/∂n = 0 where n is the outward normal. On the other hand, in order to enforce realistic outflow boundary conditions at each outlet of the domain, we consider the Windkessel model based on the electric-hydraulic analogy [28]. By representing the blood pressure and flow rate with voltage and current, respectively, and by describing the effects of friction and inertia in blood flow and of vessel elasticity with resistance R, inductance L and capacitance C, respectively, the methods for analysis of electric circuits can be borrowed and applied to the investigation of cardiovascular dynamics. In this work, we consider a three-element Windkessel RCR model [33]. It consists of a proximal resistance R p,k , a compliance C k , and a distal resistance R d,k , for each outlet k (Fig. 1).
The downstream pressure, p k , is expressed through the following differential algebraic equations (DAE) system: where Q k is the flow rate, and p p,k and p d,k are the proximal and the distal pressure, respectively. For the space discretization of problems (1)-(2), we adopt the Finite Volume (FV) approximation. A partitioned approach has been used to deal with the pressure-velocity coupling. In particular a Poisson equation for pressure has been used. This is obtained by taking the divergence of the momentum equation (1) and exploiting the divergence free constraint (2), We have used the PISO algorithm [19] employed in the finite volume C++ library OpenFOAM ® [32]. For more details, we refer the reader to [15].

The reduced order model
The reduced order model we propose is the so-called proper orthogonal decomposition with interpolation. It is a technique widely used within the reduced order modeling community in the study of parametric problems. POD allows to extract, from a set of high-dimensional snapshots, the optimal basis which minimizes the error between the original snapshots and their orthogonal projection. The data-driven approach used in this work is based only on data and does not require knowledge about the governing equations that describe the system (and which generated the snapshots). It is also non-intrusive, i.e., no modification of the simulation software is carried out. Still, there are works that use non-intrusive methods that are not data-driven (see, e.g., [35]). The original snapshots are projected onto the POD space in order to reduce their dimensionality. Then the solution manifold is approximated using an interpolation technique. In this work, we will use a radial basis function (RBF) interpolation [22]. Several examples of applications based on this so-called POD with interpolation (PODI) [9] technique can be found in literature, in a wide range of contexts (see, e.g., [13]).
We are going to briefly describe the method that consists in two phases: -Offline: let N denote the number of degrees of freedom, e.g. associated to the FV discretization introduced in the previous section. Let ϕ i , with i = 1, . . . , N s , be the snapshots related to a generic variable of interest collected by solving the high-fidelity problem, with different values of the input parameters π i , resulting in N s input-output pairs (π i , ϕ i ). The snapshots matrix S is assembled by arranging the snapshots as columns, i.e. S = [ϕ 1 , ϕ 2 , . . . , ϕ Ns ]. By applying the singular value decomposition to this matrix, we have: where U ∈ A N ×Ns is the unitary matrix containing the left-singular vectors, Σ ∈ A Ns×Ns is the diagonal matrix containing the singular values λ i , and V ∈ A Ns×Ns , with the symbol * denoting the conjugate transpose. The left-singular vectors, namely the columns of U = [φ 1 , φ 2 , . . . , φ Ns ], are the so-called POD modes. In order to reduce the dimensionality of the problem, we can keep the first k modes to span the optimal space with dimension k to represent the snapshots. By considering that the singular values are returned in decreasing order, we could truncate the number of modes simply selecting the first k columns of U .
Therefore, the matrices U k ∈ A N ×k , Σ k ∈ A k×k , V k ∈ A Ns×k in Eq. 5 are the truncated matrices with rank k. After constructing the POD space, we can project the original snapshots onto this space. We compute C ∈ R k×Ns as C = U k T S, where the columns of C are the so-called modal coefficients. We express the input snapshots as a linear combination of the modes using such coefficients. Then, we have: where α j (π i ) are the elements of C. At a given mode φ j , the (π i , α j (π i )) pairs sample the solution manifold in the parametric space. The interpolation of the modal coefficients α j (π i ) in the parameter space is carried by using RBF functions. It is based on the following formula: where w j,m are proper weights and ζ j,m are the RBF functions which are chosen to be Gaussian functions, centered in π m . For the computation of the weights w j,m , the following property has to be used: The last equation can be rewritten in form of a linear system: Thus, one could solve the latter linear system to obtain the weights w j related to the mode φ j which will be stored to be then used in the online stage.
-Online: we are able, for any new parameter value π to calculate the new coefficients α j (π ), which are given simply by: Then, we compute the high-dimensional solution by projecting back the (approximated) modal coefficients to the original space: We remark that the procedure can be repeated for several variables of interests. Furthermore, it is not necessary for such a variable to be an unknown of the original system (such as velocity and pressure); indeed, we will use the PODI technique not only for primal quantities, but also for derived quantities such as wall shear stress (WSS).
Regarding the technical implementation of the PODI method, we use the Python package called EZyRB [12].

Numerical results and discussion
In order to test the performance of the presented computational pipeline, we investigate the aortic blood flow pattern in presence of a Left Ventricular Assist Device (LVAD) when varying the pump flow rate. This case study has been thoroughly discussed both at high-fidelity (or full order model, FOM) and ROM level in [15]. Here, after summarizing some relevant computational details related to the clinical data, geometrical model and full order simulations (Sec  Fig. 2a), the ascending aorta (with inlet in the center of Fig. 2a, which is observed from below in Fig. 2b), brachiocephalic artery, right subclavian artery, right common carotid artery, left common carotid artery, left subclavian artery and descending aorta, as shown in Fig. 2. We consider a tetrahedral computational grid with h min = 5.83e − 4 and h max = 3.2e − 3 for a total of 200k cells. The quality of this mesh is suitable for a FV solver: it features very low values of average non-orthogonality (30 • ) and skewness (around 1). Fig. 2 shows the mesh. It should be noted that such a computational grid has been used in [15] where a mesh convergence analysis is carried out.
By making reference to the equations (1)-(4), the convective term is discretized by using a first order upwind scheme. On the other hand, for the diffusive term, a central differencing interpolation scheme with non-orthogonality correction is preferred. Regarding the pressure gradient, we use a linear interpolation scheme. For more details about such schemes, the reader can make reference to [32,20]. Finally, to discretize in time the equations (1) and (3), we adopt Backward Differentiation Formula of order 1 (BDF1), see e.g. [27].
Coefficients values of the Windkessel models used for the enforcement of the outlet boundary conditions shown in Table 1 are based on [15].
Although the high-order simulations are obtained through time-stepping by solving the governing equations (1)-(4), we are interested in computing and collecting the steady states, i.e. solutions where ∂ t u vanishes.  Table 2. Cumulative energy of the eigenvalues for pressure p, wall shear stress WSS, and velocity components, u x , u y and u z .

p
WSS u x u y u z E X 0.5% 7.7% 8.5% 12.2% 11.4% Table 3. L 2 norm relative errors for pressure p, wall shear stress WSS, and velocity components, u x , u y and u z , for P F = 4 l/min.

ROM results.
To train the ROM, we consider the range of pump flow rate (i.e., the flow rate at the inlet of outflow cannula) P F ∈ [3,5] that covers typical clinical values. Thus, we consider as parameter π the pump flow rate P F . In particular, we choose equispaced distributions inside the ranges P F ∈ [3, 3.8] and P F ∈ [4.2, 5]. The sampling frequency is 0.2 for both ranges, so that we have a database including 10 snapshots related to the high-fidelity steady state solutions. It should be noted that in [15] we have showed that, for this benchmark, the number of snapshots does not affect significantly the accuracy of the ROM. One new value of P F in which the ROM has not been trained but which belongs to the range of the training space, P F = 4, is used to evaluate the performance of the parametrized ROM. In [15] a equispaced distribution of 11 snapshots inside the range P F ∈ [3,5] was used, i.e. the snapshot related P F = 4 was included in the FOM database, and ROM was performed for P F = 3.45 and P F = 4.35. We note that, in [15], the distance between the parameter values for which the ROM is performed and the nearest snapshot is 0.05. On the other hand, here such a distance is larger, 0.2. Table 2 shows the cumulative energy of the eigenvalues for pressure p, wall shear stress WSS, and velocity components, u x , u y and u z . In order to retain the 99% of the system's energy, 1 mode for p, 1 for WSS, 2 for u x , u y and u z are selected. It has been verified that considering a larger number of POD modes does not increase the accuracy of the ROM. To provide some quantitative results, the relative error in the L 2 -norm, calculated as where X F OM is the value of a particular field in the FOM model, and X ROM the one that is calculated using the ROM, is considered. In Table 3, the relative error for all the variables of interest is reported. Fig. 3 displays a comparison between FOM and ROM for p and WSS, and for the velocity related to a section of the ascending aorta next to the anastomosis location. The comparison indicates that the ROM is able to provide a good reconstruction for all the variables. Finally, we comment on the computational costs. The CPU time required by a FOM simulation is 9600s and the one of the ROM, that is related to the computation of the modal coefficients and reconstruction of the fields, is 40s. This corresponds to a speed-up of ≈ 240, that demonstrates the fact that it is possible to use the ROM in the place of the FOM in order to obtain accurate simulations with a significant reduction of the computational cost.

Web application.
Due to the aforementioned speedup, research activities based on techniques (e.g., ROMs) that leads to technological innovation for real time calculation is acquiring considerable relevance and popularity also in the biomedical field. The combination of ROMs with technological development through a web interface would allow real time data to be accessed in hospitals and operating rooms on portable devices. In this scenario, the web server ARGOS [2] has been created, which has the task of proposing a platform to favor a more widespread exploitation of real time computing through a simple "click". It is a very intuitive and smooth web platform, which does not need a strong experience in numerical analysis, fluid dynamics or scientific computing field to be used. ARGOS offers a wide variety of applications related to several problems and in particular it contains the section ATLAS [3] focused on the cardiovascular field. Here, we are going to provide a brief description of the web application under development which has the aim to support the user (which in this case could either be a scientist involved in the manufacturing of the pump, or a medical doctor interested in evaluation the hemodynamics in different operating scenarios) to set the LVAD device according to the need of the patient. Figure 4 a) displays a screenshot of the application under development [24]. On the left side, the user can set up the pump configuration. We have two control panels related to two different settings denoted as Panel 1 (Figure 4 b)) and Panel 2 (Figure 4 c)). In Panel 1, required input data are the pressure head ∆P and pump speed ω, and the corresponding output is the pump flow rate P F . The target user for this panel is a scientist involved in the design of the pump (i.e., for instance, by changing the pressure head ∆P ). In Panel 2, the user provides a measured pump speed ω and the corresponding measured pump flow rate P F , and the application returns the corresponding pressure head ∆P . Then, for this value of ∆P , it is possible to vary the value of ω to obtain a different value of P F . The target user for this panel is a medical doctor, who reads measured values of pump speed ω and flow rate P F during an LVAD ramp test, and is then interested in predicting possible hemodynamics outcomes when changing the pump configuration set during the ramp test. The relationship between ω, ∆P and P F used in the application is given by the following analytical relationship where K A , K B and K C are constants which depend on pump design (Table 4), that provides an acceptable fit as showed in Fig. 5. The ∆P − P F analytical curve also is displayed in the application, below the control panels, and updates in real time when the pump speed ω is changed by the user. Currently, the application can be used only for P F ∈ [3,5] because the FOM snapshots are collected for such a range of values. Then, if the values of ω and ∆P selected are such that the corresponding P F is outside the range [3,5], an error message will appear. On the right part of the screen, we report a brief description of the web application, including relevant references, logos and acknowledgements, and we visualize in real time the solutions provided by the ROM related to the P F value corresponding to the pump setting for the variables of interest: pressure, velocity and wall shear stress.

Conclusion and perspectives
In this work, an efficient non-intrusive data-driven reduced order modelling to be used within hemodynamics framework is presented. The FOM is represented by the incompressible Navier-Stokes equations discretized by using a FV technique. Furthermore, the development of the ROM   is carried out by using the proper orthogonal decomposition with interpolation (PODI). The online phase of the ROM results in a data-driven approach which is based only on data and does not require knowledge about the governing equations that describe the system. It is also non-intrusive, i.e. no modification of the simulation software is required. For this reason it is particularly versatile thanks to its capability to be coupled also with commercial solvers. Moreover, we have presented a preliminary web application through which one can run the ROM by using a very user-friendly interface, without the need of having a specific numerical expertise, and thus possibly widening the use of numerical tools to practitioners. The benchmark we have chosen to test the efficiency of our algorithm is represented by the aortic blood flow pattern in presence of a Left Ventricular Assist Device (LVAD) when varying the pump flow rate. We show that the ROM provides accurate solutions with a significant reduction of the computational cost, up to at least two orders of magnitudes. As a follow-up of the present work, we are going to make further efforts in order to improve the ease of use of the web application. We are also moving towards geometrical parametrization in the context of patient-specific geometries, extending e.g. the work carried out in [5] to different problems and different model reduction techniques. Finally, we are interested in improving the full order model, by considering turbulence effects (see, e.g., [16,18]), as well as coupling the fluid model with an elasticity model to simulate fluid-structure interaction (FSI) (see, e.g., [17,7,23,26]). This would make the presented pipeline more complete and versatile.