A parallel domain decomposition algorithm for fluid-structure interaction simulations of the left ventricle with patient-specific shape

1.   Introduction

The Biomechanics properties of the human heart have been attracting interests in both clinical and research communities, as they can reveal the fundamental physiological status of the heart, and are useful to predict and diagnose its diseases, such as myocardial ischemia and infarction, atrial and ventricular arrhythmias, systolic and diastolic dysfunctions, and so on ^[1,2]. Numerical simulation provides a practical way to study the cardiac mechanics ^[3,4], by which the kinetical and kinematical behaviors of the heart muscles and the enclosed blood are depicted via a system of partial differential equations. Within the framework, the interested mechanical features are summarized into corresponding model parameters, and can be quantified by solving an inverse problem. Following this path, a series of induced forward sub-problems are repeatedly solved, and the values of the interested parameters are approximated to match the prior knowledges in an optimal sense ^{[5,6,7,8,9,10,11,12,13,14,15,16]}.

In designing a numerical solver for such a system, a critical issue lays on the handling of the fluid-structure interaction ^[17,18,19] between the heart muscle and the enclosed blood. One popular approach is the Arbitrary Lagrange-Euler (ALE) method ^[20,21], in which the fluid sub-system is described in the Euler coordinate system, the structure sub-system is described in the Lagrange coordinate systems, and the relevant field variables are coupled through the fluid-structure interface continuation conditions. Since the fluid-structure interface is explicitly treated, the ALE method can achieve a good accuracy to approximate the physical quantities that are defined on the surface, such as the wall shear stress. Another type of methods, including the immersed boundary(IB) and its variants such as the immersed boundary/finite element (IB/FE) method ^{[22,23,24,25,26,27]}, works on a fixed mesh but requires a very fine mesh in the neighborhood of the interface to confine the interpolation error between the fluid mesh and structure mesh, which inevitably leads to a huge number of degrees of freedom. Watanabe et al. ^[28,29]. developed a program based on the finite element method in order to simulate the fluid-structure interaction of the human left ventricle during systole and diastole. Doyle et al. ^[30] used an ideal geometry to perform a fluid-structure interaction simulation of canine left ventricular blood flow, and analyzed in detail the parallel scalability and stability of their method. Nordsletten et al. ^[31] reconstructed the geometry of left ventricle based on Magnetic Resonance Imaging (MRI) and simulated the fluid-structure interaction of the heart to investigate the blood flow pattern and the pressure distribution inside the left ventricle. Gao et al. ^[32,33] used the IB/FE method to simulate the dynamic process of the left ventricle from end-diastole to end-systole in a heart beat, whose results show a very good match to the clinical observation.

In this paper, we propose an efficient parallel domain decomposition algorithm to solve the cardiac FSI system, and demonstrate the simulation results of large-scaled problems on a supercomputer. More specifically, we verify the computation efficiency, numerical stability and the parallel scalability of the proposed method, and demonstrate the fidelity of the simulated results by comparing them with clinical data. The rest of this paper is organized as follows: Section 2 briefly describes the reconstruction of the heart geometry from a medical image and the finite element mesh generation. The governing equations, boundary conditions, numerical solution algorithm, as well as the parameter settings are also introduced; in Section 3, we present the numerical experiments on the mesh size convergence, time step size convergence, and the parallel scalability of the solving algorithm, as well as give detailed analyses for a one-heart-beat simulation; the proposed method is summarized in Section 4, in which several possible improvements are also suggested for the future work.

2.   Model and method

2.1.   Geometry reconstruction and mesh generation

A CT image of a 38-year-old male is adopted to construct the 3D geometry model for our numerical simulations, which was taken at the end of a diastole phase with spatial resolution 0.5 mm. The axial, coronal and sagittal views of CT images are shown in Figure 1(a-c), respectively, in which the red region represents the blood contained in the left ventricle, while the myocardium is colored in cyan. Figure 1(d) shows the initial shape of the left ventricle obtained from the CT image reconstruction process, which is smoothed and repaired to get the final geometry as shown in Figure 2(a). In our tests, five meshes, denoted by ${\mathcal{L}}_{i}(i = \mathrm{1, 2}, \mathrm{3, 4}, 5)$, are generated for investigating the parallel efficiency and mesh convergence. The total number of mesh nodes and elements, the elements in the ventricle and cavity domains, and the element sizes are listed in Table 1. In particular, Figure 2(b) and (c) show the mesh structures of myocardium part and blood part of ${\mathcal{L}}_{1}$ in separate views.

2.2.   Governing Equations

Define the computational region at time t as $\Omega ^{t} = \Omega _{f}^{t}\cup \Omega _{s}^{t} \subset {R}^{3}$, where $\Omega _{f}^{t}$ is fluid domain, $\Omega _{s}^{t}$ is structure domain, and ${\Gamma }_{w}^{t} = {\partial \Omega }_{f}^{t}\cap {\partial \Omega }_{s}^{t}$ is the interface between fluid and structure. In addition, ${\Gamma }_{f, I}^{t}$ and ${\Gamma }_{f, O}^{t}$ are the inlet and outlet boundaries of the fluid domain at time 𝑡, and ${\Gamma }_{s, I}^{t}$ and ${\Gamma }_{s, O}^{t}$ are the inlet and outlet boundaries of the structure domain at time 𝑡. In particular, the initial configuration of each of the above is denoted by letting $t = 0$, as shown in Figure 3.

The myocardium of the left ventricle is modeled as a linear elastic material. Namely, the movement of each material point in the myocardium region is described by the following equation:

where ${\bf{x}}_{s}$ represents the displacement of a material point in ${\Omega }_{s}^{0}$, ${\rho }_{s}$ is density of myocardium, ${\boldsymbol{\sigma }}_{s} = {\mathrm{\lambda }}_{s}\left(\nabla · {\bf{x}}_{s}\right)\boldsymbol{I}+{\mu }_{s}\left(\nabla {\bf{x}}_{s}+\nabla {{\bf{x}}_{s}}^{T}\right)$ is the Cauchy stress tensor in which $\boldsymbol{I}$ denoting the rank two identity tensor, $\alpha$ is a damping parameter, ${\mathrm{\lambda }}_{s}$ and ${\mu }_{s}$ are two Lamé coefficients, which are generally calculated from provided material Young's modulus $E$ and Poisson's ratio ${v}_{s}$ : ${\mathrm{\lambda }}_{s} = {v}_{s}E/\left(\left(1+{v}_{s}\right)\left(1-2{v}_{s}\right)\right)$, ${\mu }_{s} = E/\left(2\left(1+{v}_{s}\right)\right)$.

On the other hand, we use an unsteady incompressible Navier-Stokes equation to depict the blood flow. Note that, because of the deformation of the myocardium, the fluid domain occupied by the blood changes in time. To deal with this issue, an ALE method ^[34] is used to describe the fluid domain $\Omega _{f}^{t}$ at time $t$ :

Where $\mathrm{Y}$ represents a reference point in the initial configuration of the fluid domain, and ${\bf{x}}_{f}$ is its displacement. In practice, the reference points match the mesh vertices, and ${\bf{x}}_{f}$ are assumed to satisfy the following Laplace equation:

Under the above settings, the ALE form of the fluid equation reads:

where ${\rho }_{f}$ is fluid density, ${\bf{u}}_{f}$ is fluid velocity, and ${\boldsymbol{\sigma }}_{f} = -{p}_{f}\boldsymbol{I}+{\mu }_{f}\left(\nabla {\bf{u}}_{f}+{\nabla {\bf{u}}_{f}}^{T}\right)$ is Cauchy stress tensor, ${p}_{f}$ is fluid pressure, ${\mu }_{f}$ is fluid viscosity. ${\boldsymbol{\omega }}_{g} = \partial {\bf{x}}_{f}/\partial t$ represents mesh moving velocity, and the first term in Eq (2) ${\left.(\partial {\bf{u}}_{f}/\partial t)\right|}_{Y}$ denotes the time derivative of velocity is taken as the ALE coordinates are unchanged.

Besides, the following conditions are satisfied at the fluid-structure interface:

where ${\bf{n}}_{s}$ and ${\bf{n}}_{f}$ are the unit outer normal vectors in structure and fluid domains, respectively.

2.3.   Material parameters and boundary conditions settings

For all experiments performed in this paper, the duration of a cardiac cycle was set to 0.80 s, with the diastolic phase lasts for 0.48 s and the systolic phase for 0.32 s, if not mentioned otherwise. Fluid density is $1.05g/{cm}^{3}$, and viscosity is $0.03g/\left(cm · s\right)$, structure density is $1.37g/{cm}^{3}$, Young's modulus is $6\times {10}^{6}g/\left(cm · {s}^{2}\right)$, and Poisson's ratio is $0.49$. Each numerical simulation starts from a diastolic phase, during this stage, the boundary conditions are given in Eq (5). That is, the inlet of left ventricle, ${\Gamma }_{f, I}^{0}$ keeps open, and the blood flow passing this boundary according to a prescribed rate, as the time-variant curve shown in Figure 4. Note that, this plot indicates two peaks of the inflow, the first is the early diastolic filling wave (E wave) and the second is the late diastolic filling wave (A wave). In the meantime, the outlet of left ventricle is closed, thus the blood flow rate on ${\Gamma }_{f, O}^{0}$ keeps at zero. Boundary conditions of the systolic phase are prescribed by Eq (6): during this stage, the inlet of left ventricle is closed (no blood inflow), while its outlet is open and the contained blood flows out freely in a pressure-free state.

2.4.   Numerical algorithms

We used a finite element method to spatially discretize the model equations, and solve the resulting system in a monolithic way. More specifically, the solid equation is discretized by the conventional P1 finite, and the fluid domain equations are discretized by a P1-P1 finite element that is enhanced with stabilization terms ^[35]. Further, a one-step backward-difference method ^[36] is used to get a fully discretized system, of which the solution vector is denoted by $x$. To solve the system, an inexact Newton method is invoked at each time step: Denote the nonlinear algebraic system as $\mathcal{F}\left(x\right) = 0$, and given the initial value ${x}^{\left(0\right)}$, the following iterations are performed:

Here, ${\alpha }^{\left(k\right)}$ is the step size obtained by a line search step, and the search direction ${s}^{\left(k\right)}$ is obtained by using Krylov subspace method to approximately solve the Jacobian system as follows:

where ${J}_{k} = \nabla \mathcal{F}\left({x}^{\left(k\right)}\right)$ is the Jacobian matrix of the nonlinear function $\mathcal{F}\left(x\right)$ at point ${x}^{\left(k\right)}$. Due to the intrinsic properties of Navier-Stokes equations, ${J}_{k}$ is non-symmetric, so the GMRES method is used to approximately solve the linear system (8). ${\eta }_{k}$ is the termination condition of the linear iteration, and ${10}^{-6}$ is chosen in this paper.

In above, one key argument is on the construction of ${M}_{k}$, which represents a precondition operator for accelerating the convergence of Krylov iterations. In this paper, we use the Restricted Additive Schwarz (RAS) precondition method ^[37] to construct ${M}_{k}$. To do so, the whole mesh $\Omega$ is divided into N non-overlapping sub-domains:

then each sub-domain $\Omega _{l}$ is extended by $\delta$ layers to obtain overlapped sub-domain ${\Omega }_{l}^{\delta }$. The preconditioner ${M}_{k}$ is constructed by:

where ${R}_{l}^{\delta }$ and ${R}_{0}^{\delta }$ denote the restriction operators of $\Omega$ to ${\Omega }_{l}^{\delta }$ and ${\Omega }_{l}$ :

for $v$ represents a vector defined on the global mesh domain. Beside, their transposes represent the corresponding prolongation operators, ${B}_{l}^{-1}$ denotes the subproblem solver defined on ${\Omega }_{l}^{\delta }$, which is set as an in-complete LU decomposition operator of the subproblem Jacobian matrix in this paper. More details on the implementation of the preconditioner can be found in ^[32].

For clarity, we summarize Newton-Krylov-Schwarz(NKS) algorithm for solving the system of nonlinear algebraic equations of the FSI problem in each time step as follows.

3.   Numerical experiments

3.1.   Mesh convergence tests

It is well-known that, when solving the above fluid-structure interaction problem, the smaller size of mesh element, the higher accuracy of the solution will achieve, which, however, consumes more computational resources and time. In this section, we discuss the mesh convergence of the proposed algorithm and analyze the effect of the size of the mesh element on the solution accuracy. For this purpose, four meshes, ${\mathcal{L}}_{1}, {\mathcal{L}}_{2}, {\mathcal{L}}_{3}, {\mathcal{L}}_{4}$ are used to repeat a same simulation for one cardiac cycle. The time step is set to 0.0005 s, and 20 nodes of Tianhe-2A supercomputer (a total of $20\times 24 = 480$ processor cores) are used for each case.

Time and memory cost by using each of the meshes are listed in Table 2, where "time" represents the total compute time for all time steps and "memory" is the maximum memory allocated by all processor cores during the simulation. Effect of the element size on the compute time and memory cost is given in Figure 5. As it is expected, as mesh element size decreases, both time and memory required for simulation get increased.

To verify our results, the calculated left ventricular cavity volume and the blood flow rate at its outlet by using different meshes are compared: Let ${V}_{i}\left(t\right)$ be the left ventricular cavity volume at time $t$ obtained by using mesh ${\mathcal{L}}_{i}(i = \mathrm{1, 2}, \mathrm{3, 4})$, and ${Q}_{i}\left(t\right)$ the flow rate through the outlet. We regard the results of the finest mesh ${\mathcal{L}}_{4}$ as a reference, and calculate the deviations of the cavity volume and outlet flow rate of other coarser meshes. Deviations at several representative moments are listed in Table 3 and 4, and the overall pattern are plotted in Figures 6 and 7 for the whole cardiac cycle. All of these clearly show that the deviation becomes smaller as the element size decreases, which indicates that the solution accuracy gets improved.

3.2.   Time-step convergence

In order to verify the convergence of time step size, we use ${\mathcal{L}}_{1}$ mesh and set the size of the time step to be ${\Delta t}_{1} = 0.005, {\Delta t}_{2} = 0.002, {\Delta t}_{3} = 0.001, {\Delta t}_{4} = 0.0005$ (unit: second) to perform simulations, respectively. ${V}^{i}\left(t\right)$ denotes the volume occupied by the blood contained in left ventricle at time $t$ by using time step ${\Delta t}_{i}$, and ${Q}^{i}\left(t\right)$ denotes the rate of the blood flow that passes through the left ventricular outlet. Taking the results computed by using the smallest ${\Delta t}_{4}$ as the reference, the deviations of left ventricle blood volume and outlet flow rate obtained by using other time steps at several representative moments are presented in Tables 5 and 6, and the overall patterns during the cardiac cycle are presented in Figures 8 and 9. Based on these results, it can be seen that the difference of different computing results becomes smaller as the time step decreases, as desired.

3.3.   Parallel performance

To investigate the parallel scalability of the proposed algorithm, two computational meshes, ${\mathcal{L}}_{4}$ ($4.06\times {10}^{6}$ mesh elements) and ${\mathcal{L}}_{5}$ ($7.92\times {10}^{6}$ mesh elements), are chosen to repeat simulations on the Tianhe-2A supercomputer, by using 72, 144, 288, 576, 1152, and 2304 processor cores, respectively. For each case, the average compute time of the first 10 time steps is counted to calculate the speedup. More specifically, let ${T}_{ini}$ be the average computation time per time step when using 72 processors, and ${T}_{p}$ the average computation time per time step when using $p$ processors, the speedup and parallel efficiency are defined as follows:

Table 7, Figures 10 and 11 present the parallel scalability performance of the numerical algorithm for solving the left ventricle fluid-structure interaction problem, where "np" denotes the number of processor cores, "Newton" the average number of Newton iteration steps per time step, "GMRES" the average number of GMRES iteration steps in each Newton step, "Time" the average compute time per time step, "Speedup" the speedup of the algorithm, and "Efficiency" the parallel efficiency of the algorithm. It can be seen that when the number of processor cores gets increased while the problem size keeps unchanged, the number of Newton iteration steps remains almost the same and the number of GMRES iterations increases just few steps, and the total compute time decreases significantly. Moreover, as the number of mesh elements is increased from $4.06\times {10}^{6}$ to $7.92\times {10}^{6}$, both speedup and parallel efficiency are improved, indicating that the algorithms can achieve higher parallel efficiency when the scale of the problem matches appropriately the number of processor cores. Nevertheless, it can be found that, when for case of $7.92\times {10}^{6}$ elements, the algorithm still has 40% parallel efficiency when the number of processor cores is expanded to 2304 cores. These results verify that the proposed algorithm has good parallel scalability for solving the left ventricle fluid-structure interaction problem.

3.4.   Analysis on simulation results

In the following, we use 480 processors on Tianhe-2A supercomputer to perform a cardiac cycle simulation based on ${\mathcal{L}}_{4}$, and the time step size is set to 0.0005 s. Due to the simplification of the model, we mainly discuss the velocity field distribution in the blood region contained left ventricle, as well as the deformation and stress status of the myocardium.

Figure 12 shows the rate of blood flow at the inlet and outlet of left ventricle as a time-dependent function during a cardiac cycle, and Figure 13 shows instantaneous blood flow velocity streamline in left ventricle at several time moments. It can be seen that the blood flow velocity increased rapidly at the beginning stage, and the first local maximal flow rate appears as the early diastolic filling wave (E wave) arrives. Besides, an obvious vortex formed at 0.08s can be observed. After that, the inlet flow rate drops, but the previously observed vortex area gets expanded, with some new vortexes generated. At 0.25 s, the inlet flow velocity drops to a low level, the velocity field also gradually gets stabilized, and the vortex area decrease to almost vanished. The late diastolic filling wave (A wave) arrives around 0.4 s, which causes the inlet flow rate increased again. It can be seen that the velocity field in the left ventricle fluctuated, and a vortex appears again. The diastole phase terminates at 0.48 s, and the systole begins right away, during which, the left ventricle inlet keeps closed and the outlet is open. The outlet blood flow velocity quickly reaches its peak, and then falls gradually. During the systolic phase, a vortex can also be observed at 0.70 s.

In addition, Figure 14 presents the history of the cavity volume of the left ventricle during one cardiac cycle, and Figure 15 demonstrates the instantaneous velocity field of the blood as well as the displacement field of the myocardium. According to these results, it can be seen that the velocity field inside the left ventricle fluctuates with the inlet flow rate throughout diastole, and the volume of the left ventricular cavity and the displacement of the myocardium undergo two drastic changes after the arrival of the E and A waves, as large amount of blood rapidly flushes into the left ventricle. Besides, during the systolic phase, the outlet blood flow reached its peak at the beginning of contraction and then decrease gradually, and the velocity near the outlet region was significantly higher than elsewhere. As expected, after blood flowed out of left ventricle, the volume of the left ventricular cavity and the displacement of the myocardium around the left ventricle get decreased.

Furthermore, we cut the left ventricle on three representative cross section along the axial axis, coronal axis and sagittal axis, respectively, and compare the simulation results with the CT image data at the end of the systole, as shown in Figure 16. It can be seen that the results obtained by simulation are within the normal physiological range and visually match the CT images. But at the same time, due to the simplicity of the underlying model and experimental settings (e.g., the left ventricle is linearly elastic and the active contraction is ignored), some artifacts and flaws can be found in the current simulation results. For example, as shown in Figure 14, the volume of the left ventricle cavity at the end of the drainage stage is larger than the initial one, and we found that the blood flowing into the left ventricle during diastole is larger than that flowing out from during the drainage stage (74 mL vs 56 mL), indicating an un-normal situation. In addition, backflow has already appeared at the outlet at the end of contraction, and the backflow speed is relatively high. At this moment, the blood flow rate is about 20cm/s, but the speed of backflow has reached about 80 cm/s. If the simulation continues, the backflow speed may reach up to 1000 cm/s, which does not lay in a physiological range. In the future, we will revise our model, algorithm and experimental settings, in aiming to get more realistic simulation results.

4.   Conclusions

Due to the complexity of human heart, fluid-structure interaction simulation on the hemodynamic is a very challenging task. In this paper, we proposed a high-performance parallel solver to perform efficient FSI simulations of the human heart left ventricle on Tianhe-2A supercomputer. Through experiments, we verify that the algorithm is highly efficient, and has good numerical stability and parallel scalability up to more than 2000 processor cores. By using the solver, we further provide a preliminary demonstration of its possible application scenarios in clinical research.

Acknowledgements

This work was partially supported by the National Key R & D Program of China 2018YFE0198400, NSFC grants under 12001520, 12071461 and 62161160312, and Shenzhen Fund RCYX20200714114735074.

Conflict of interest

The authors declare that there is no conflict of interest.