Performance of preconditioned iterative linear solvers for cardiovascular simulations in rigid and deformable vessels

TL;DR

This paper evaluates various iterative linear solvers and preconditioners for cardiovascular simulations involving rigid and deformable vessels, highlighting the superior performance of incomplete factorization preconditioners and the BIPN method in specific scenarios.

Contribution

It provides a comprehensive comparison of solver and preconditioner performance in realistic cardiovascular models, introducing insights into their scalability and effectiveness.

Findings

Incomplete factorization preconditioners offer best scalability.

BIPN outperforms other methods in rigid wall models.

BICG with diagonal and LU preconditioners excels in deformable wall models.

Abstract

Computing the solution of linear systems of equations is invariably the most time consuming task in the numerical solutions of PDEs in many fields of computational science. In this study, we focus on the numerical simulation of cardiovascular hemodynamics with rigid and deformable walls, discretized in space and time through the variational multi-scale finite element method. We focus on three approaches: the problem agnostic generalized minimum residual (GMRES) and stabilized bi-conjugate gradient (BICGS) methods, and a recently proposed, problem specific, bi-partitioned (BIPN) method. We also perform a comparative analysis of several preconditioners, including diagonal, block-diagonal, incomplete factorization, multi-grid, and resistance based methods. Solver performance and matrix characteristics (diagonal dominance, symmetry, sparsity, bandwidth and spectral properties) are first examined for an idealized cylindrical geometry with physiologic boundary conditions and then successively tested on several patient-specific anatomies representative of realistic cardiovascular simulation problems. Incomplete factorization pre-conditioners provide the best performance and results in terms of both strong and weak scalability. The BIPN method was found to outperform other methods in patient-specific models with rigid walls. In models with deformable walls, BIPN was outperformed by BICG with diagonal and Incomplete LU preconditioners.

Tables

Table 1. Table 1: l 2 superscript 𝑙 2 l^{2} -norm of velocity errors from different residual norm tolerances. The solution error is obtained from nodal velocity errors after 10 time steps, using a reference simulation with a tolerance up to the machine precision, ϵ = 10 − 12 italic-ϵ superscript 10 12 \epsilon=10^{-12} .

	$ϵ={10}^{-1}$	$ϵ={10}^{-3}$	$ϵ={10}^{-6}$	$ϵ={10}^{-9}$
$Err$	0.2235	0.0028	1.618$\times {10}^{-7}$	1.107$\times {10}^{-9}$

Table 2. Table 2: The number of iterations and the portion of compute time consumed by the linear solver for the pipe benchmark problem with rigid walls. The number of linear solver iterations is counted for 10 time step calculations. t T T {}_{\text{T}} is the total compute time, t LS LS {}_{\text{LS}} is the compute time consumed by solving the linear system.

BIPN-RPC /$ϵ={10}^{-3}$
${\text{N}}_{\text{bipn}}$	${\text{t}}_{\text{LS}}/{\text{t}}_{\text{T}}$	${\text{N}}_{\text{gmres}}$	${\text{t}}_{\text{gmres}}/{\text{t}}_{\text{LS}}$	${\text{N}}_{\text{cg}}$	${\text{t}}_{\text{cg}}/{\text{t}}_{\text{LS}}$
129	89.5$\%$	1475	5.3$\%$	38582	77.6$\%$

Table 3. Table 3: The number of iterations and the portion of compute time consumed by the linear solver for the pipe benchmark problem with deformable walls. The number of linear solver iterations is counted for 10 time step calculations. t T T {}_{\text{T}} is the total compute time, t LS LS {}_{\text{LS}} is the compute time consumed by solving the linear system.

BIPN-RPC /$ϵ={10}^{-3}$
${\text{N}}_{\text{bipn}}$	${\text{t}}_{\text{LS}}/{\text{t}}_{\text{T}}$	${\text{N}}_{\text{gmres}}$	${\text{t}}_{\text{gmres}}/{\text{t}}_{\text{LS}}$	${\text{N}}_{\text{cg}}$	${\text{t}}_{\text{cg}}/{\text{t}}_{\text{LS}}$
608	78.6$\%$	117377	77$\%$	64865	15.6$\%$

Table 4. Table 4: Number of nodes, elements and non-zero entries in tangent stiffness matrix for selected meshes in scaling studies.

	Fluid
	$\#$ Nodes	$\#$ Tetrahedra	$\#$ Non-zeros
Mesh1	93,944	551,819	1,496,375
Mesh2	166,899	1,002,436	2,718,385
Mesh3	349,469	2,131,986	5,784,902
Mesh4	716,298	4,415,305	11,990,977
Mesh5	1,314,307	8,117,892	22,144,741

Table 5. Table 5: The 1-norm condition number estimates and five maximum and minimum eigenvalues ( λ i subscript 𝜆 𝑖 \lambda_{i} ) of a local matrix 𝐀 f subscript 𝐀 𝑓 \mathbf{A}_{f} , and blocks 𝐊 𝐊 \mathbf{K} , and 𝐋 𝐋 \mathbf{L} without resistance BC.

${𝐀}_f$: Condition number=1.293$\times {10}^6$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max$(\times {10}^{-3})$	3.304	3.301	3.298	3.282	3.270
Min$(\times {10}^{-8})$	0.748	1.042	1.348	1.524	1.672
$𝐊$: Condition number=163
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max$(\times {10}^{-3})$	3.304	3.301	3.298	3.282	3.270
Min$(\times {10}^{-5})$	4.615	4.621	5.958	6.381	6.614
$𝐋$: Condition number=$1.555\times {10}^{18}$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max$(\times {10}^{-7})$	9.053	8.180	8.085	8.083	7.983
Min$(\times {10}^{-9})$	0.000	0.456	1.224	1.765	2.214

Table 6. Table 6: The 1-norm condition number estimates and five maximum and minimum eigenvalues ( λ i subscript 𝜆 𝑖 \lambda_{i} ) of a local matrix 𝐀 ~ f subscript ~ 𝐀 𝑓 \mathbf{\tilde{A}}_{f} and 𝐊 ~ ~ 𝐊 \mathbf{\tilde{K}} with a resistance BC.

${\widetilde{𝐀}}_f$: Condition number=3.299$\times {10}^7$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max$(\times {10}^{-3})$	86.525	3.253	3.244	3.239	3.185
Min$(\times {10}^{-8})$	0.671	1.036	1.350	1.427	1.504
$\widetilde{𝐊}$: Condition number=5.883$\times {10}^3$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max$(\times {10}^{-3})$	86.525	3.253	3.244	3.239	3.185
Min$(\times {10}^{-5})$	2.878	2.886	3.457	3.981	4.002

Table 7. Table 7: The 1-norm condition number estimate and five maximum and minimum eigenvalues ( λ i subscript 𝜆 𝑖 \lambda_{i} ) of the raw local 𝐊 s subscript 𝐊 𝑠 \mathbf{K}_{s} matrix.

${𝐊}_s$: Condition number=3.169$\times {10}^4$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	7.706	7.469	7.389	7.325	7.294
Min$(\times {10}^{-3})$	1.374	1.404	1.468	1.510	1.529

Table 8. Table 8: The 1-norm condition number estimates and five maximum and minimum eigenvalues ( λ i subscript 𝜆 𝑖 \lambda_{i} ) of the preconditioned local matrices 𝐀 f ∗ subscript superscript 𝐀 𝑓 \mathbf{A}^{*}_{f} , 𝐊 ∗ superscript 𝐊 \mathbf{K}^{*} , 𝐋 ∗ superscript 𝐋 \mathbf{L}^{*} without resistance BC and 𝐊 s ∗ subscript superscript 𝐊 𝑠 \mathbf{K}^{*}_{s} . The 5th maximum eigenvalue of local 𝐀 f ∗ superscript subscript 𝐀 𝑓 \mathbf{A}_{f}^{*} is complex, with imaginary part of − 3.257 × 10 − 5 3.257 superscript 10 5 -3.257\times 10^{-5} .

${𝐀}_f^{\ast }$: Condition number=645
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	2.271	2.248	2.248	2.247	2.244
Min$(\times {10}^{-2})$	2.557	3.558	4.623	5.183	5.795
${𝐊}^{\ast }$: Condition number=14
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	2.249	2.248	2.248	2.244	2.244
Min	0.579	0.580	0.582	0.582	0.583
${𝐋}^{\ast }$: Condition number=$1.662\times {10}^{18}$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	2.290	2.224	2.216	2.213	2.201
Min$(\times {10}^{-2})$	0.000	0.200	0.555	0.798	1.019
${𝐊}_s^{\ast }$: Condition number=1.517$\times {10}^4$
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	3.567	3.556	3.214	3.181	3.161
Min$(\times {10}^{-3})$	1.173	1.193	1.205	1.222	1.238

Table 9. Table 9: The 1-norm condition number estimates and five maximum and minimum eigenvalues ( λ i subscript 𝜆 𝑖 \lambda_{i} ) of the preconditioned local matrix 𝐀 ~ f ∗ subscript superscript ~ 𝐀 𝑓 \mathbf{\tilde{A}}^{*}_{f} , 𝐊 ~ ∗ superscript ~ 𝐊 \mathbf{\tilde{K}}^{*} , and 𝐇 ​ 𝐊 ~ 𝐇 ~ 𝐊 \mathbf{H\tilde{K}} .

${\widetilde{𝐀}}_f^{\ast }$: Condition number=9243
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	55.15	2.262	2.247	2.247	2.243
Min$(\times {10}^{-1})$	0.232	0.358	0.467	0.492	3.99
${\widetilde{𝐊}}^{\ast }$: Condition number=1862
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	55.16	2.248	2.247	2.244	2.241
Min	0.217	0.384	0.392	0.396	0.399
$𝐇\widetilde{𝐊}$: Condition number=76
${\lambda }_i$	1st	2nd	3rd	4th	5th
Max	2.248	2.247	2.244	2.241	2.241
Min	0.217	0.384	0.392	0.396	0.399

Equations

$$\text{centering}{\mathbf{A}}_f=\left[\begin{array}{cc}\mathbf{K}& \mathbf{G}\\ \mathbf{D}& \mathbf{L}\end{array}\right],\text{@add@centering}$$

$$\begin{array}{rl}\rho \frac{\partial \mathbf{u}}{\partial t}{|}_{\hat{\mathbf{x}}}+\rho \mathbf{v}\cdot \nabla \mathbf{u}& =\rho \mathbf{f}+\nabla \cdot {\sigma }_f\\ \nabla \cdot \mathbf{u}& =0\end{array}\text{in}{\Omega }_f,$$

$$\begin{array}{r}\mathbf{w}(\mathbf{x})=\sum _{a\in I_a}{N}^a(\mathbf{x}){\mathbf{w}}^a,\\ \mathbf{q}(\mathbf{x})=\sum _{a\in I_a}{N}^a(\mathbf{x}){\mathbf{q}}^a,\end{array}$$

$$\begin{array}{r}\mathbf{u}(\mathbf{x},t=t^n)={\mathbf{u}}^n(\mathbf{x})=\sum _{a\in I_a}{N}^a(\mathbf{x}){\mathbf{u}}^{a,n},\\ p(\mathbf{x},t=t^n)=p^n(\mathbf{x})=\sum _{a\in I_a}{N}^a(\mathbf{x})p^{a,n},\end{array}$$

$$\begin{gathered} {\mathbf{R}}_m^a(\dot{\mathbf{u}},\mathbf{u},\mathbf{v},p)= \\ \sum _{e\in {\mathbf{I}}_e}{\int }_{{\Omega }^e}\rho N^a\left(\dot{\mathbf{u}}-\mathbf{f}+(\mathbf{v}+{\mathbf{u}}_p)\cdot \nabla \mathbf{u}\right)d\Omega \\ +\sum _{e\in {\mathbf{I}}_e}{\int }_{{\Omega }^e}(\nabla N^a{)}^T(-p\mathbf{I}+\mu {\nabla }^s\mathbf{u}+\rho {\tau }_B{\mathbf{u}}_p\otimes ({\mathbf{u}}_p\cdot \nabla \mathbf{u}) \\ -\rho \mathbf{u}\otimes {\mathbf{u}}_p+\rho {\tau }_C\nabla \cdot \mathbf{u})d\Omega -{\int }_{{\Gamma }_h}{N}^a\mathbf{h}d\Gamma , \end{gathered}$$

$$R_c^a(\dot{\mathbf{u}},\mathbf{u},p)={\int }_{\Omega }\left(N^a\nabla \cdot \mathbf{u}-(\nabla N^a{)}^T{\mathbf{u}}_p\right)d\Omega ,$$

$$\begin{array}{r}{\mathbf{u}}_p=-{\tau }_M\left(\dot{\mathbf{u}}+\mathbf{v}\cdot \nabla \mathbf{u}+\frac{1}{\rho }\nabla p-\frac{\mu }{\rho }{\nabla }^2\mathbf{u}-\mathbf{f}\right),\\ {\tau }_M={\left(\frac{4}{\Delta t^2}+\mathbf{v}\cdot \mathbf{Gv}+C_I{\left(\frac{\mu }{\rho }\right)}^2\mathbf{G}:\mathbf{G}\right)}^{-1/2},\\ {\tau }_B=({\mathbf{u}}_p\cdot \mathbf{G}{\mathbf{u}}_p{)}^{-1/2},{\tau }_C=({\tau }_M\mathbf{g}\cdot \mathbf{g}{)}^{-1},\\ G_{ij}=\sum _{k=1}^3\frac{\partial {\xi }_k}{\partial x_i}\frac{\partial {\xi }_k}{\partial x_j},\\ \mathbf{g}\cdot \mathbf{g}=\sum _{i=1}^3{g}_i{g}_i,g_i=\sum _{k=1}^3\frac{\partial {\xi }_k}{\partial x_i},\end{array}$$

$${\dot{\mathbf{u}}}^{a,n+1}=\frac{\gamma -1}{\gamma }{\dot{\mathbf{u}}}^{a,n},{\mathbf{u}}^{a,n+1}={\mathbf{u}}^{a,n},p^{a,n+1}=p^{a,n},$$

$$\begin{array}{r}{\dot{\mathbf{u}}}^{a,n+{\alpha }_m}=(1-{\alpha }_m){\dot{\mathbf{u}}}^{a,n}+{\alpha }_m{\dot{\mathbf{u}}}^{a,n+1},\\ {\mathbf{u}}^{a,n+{\alpha }_f}=(1-{\alpha }_f){\mathbf{u}}^{a,n}+{\alpha }_f{\mathbf{u}}^{a,n+1}.\end{array}$$

$$\begin{array}{rl}\mathbf{K}\Delta \mathbf{u}+\mathbf{G}\Delta \mathbf{p}& =-{\mathbf{R}}_m\left({\dot{\mathbf{u}}}^{n+{\alpha }_m},{\mathbf{u}}^{n+{\alpha }_f},p^{n+1}\right),\\ \mathbf{D}\Delta \mathbf{u}+\mathbf{L}\Delta \mathbf{p}& =-{\mathbf{R}}_c\left({\dot{\mathbf{u}}}^{n+{\alpha }_m},{\mathbf{u}}^{n+{\alpha }_f},p^{n+1}\right),\end{array}$$

$$\begin{array}{r}{\mathbf{K}}^{ab}\approx \frac{\partial {\mathbf{R}}_m^a}{\partial \Delta {\mathbf{u}}^b},{\mathbf{G}}^{ab}\approx \frac{\partial {\mathbf{R}}_m^a}{\partial \Delta p^b},\\ {\mathbf{D}}^{ab}\approx \frac{\partial R_c^a}{\partial \Delta {\mathbf{u}}^b},{\mathbf{L}}^{ab}\approx \frac{\partial R_c^a}{\partial \Delta p^b}.\end{array}$$

$${\mathbf{A}}_f\mathbf{y}=-{\mathbf{R}}_f,$$

$${\mathbf{A}}_f=\left[\begin{array}{cc}\mathbf{K}& \mathbf{G}\\ \mathbf{D}& \mathbf{L}\end{array}\right],\mathbf{y}=\left[\begin{array}{c}\Delta \mathbf{u}\\ \Delta \mathbf{p}\end{array}\right],{\mathbf{R}}_f=\left[\begin{array}{c}{\mathbf{R}}_m\\ {\mathbf{R}}_c\end{array}\right],$$

$$\begin{array}{r}{\dot{\mathbf{u}}}^{a,n+1}\leftarrow {\dot{\mathbf{u}}}^{a,n+1}+\Delta {\mathbf{u}}^a,\\ {\mathbf{u}}^{a,n+1}\leftarrow {\mathbf{u}}^{a,n+1}+\gamma \Delta t\Delta {\mathbf{u}}^a,\\ p^{a,n+1}\leftarrow p^{a,n+1}+{\alpha }_f\gamma \Delta t\Delta p^a,\end{array}$$

$$\mathbf{K}^{ab}=\sum_{e\in\mathbf{I}_{e}}\int_{\Omega^{e}}\Big{[}\rho\alpha_{m}N^{a}N^{b}\mathbf{I}^{ab}+\rho\tilde{\alpha}_{f}N^{a}(\mathbf{v}+\mathbf{u}_{p})\cdot\nabla N^{b}\mathbf{I}^{ab}\\ +\mu\tilde{\alpha}_{f}(\nabla N^{a}\cdot\nabla N^{b}\mathbf{I}^{ab}+\nabla N^{b}\otimes\nabla N^{a})\\ +\rho\tilde{\alpha}_{f}{\tau}_{B}\mathbf{u}_{p}\cdot\nabla N^{a}\mathbf{u}_{p}\cdot\nabla N^{b}\mathbf{I}^{ab}\\ +\rho\tau_{M}\mathbf{u}\cdot\nabla N^{a}(\alpha_{m}N^{b}+\tilde{\alpha}_{f}\,\mathbf{u}\cdot\nabla N^{b})\mathbf{I}^{ab}\\ +\rho\tilde{\alpha}_{f}\tau_{C}\nabla N^{a}\otimes\nabla N^{b}\Big{]}d\Omega,$$

$$\mathbf{G}^{ab}=\sum_{e\in\mathbf{I}_{e}}\int_{\Omega^{e}}\Big{[}-\tilde{\alpha}_{f}\nabla N^{a}N^{b}+\tilde{\alpha}_{f}\tau_{M}\mathbf{u}\cdot\nabla N^{a}\nabla N^{b}\Big{]}d\Omega,$$

$$\mathbf{D}^{ab}=\sum_{e\in\mathbf{I}_{e}}\int_{\Omega^{e}}\Big{[}\tilde{\alpha}_{f}N^{a}\nabla N^{b}+\tilde{\alpha}_{f}\tau_{M}\bm{u}\cdot\nabla N^{a}\nabla N^{b}\\ +\tau_{M}\nabla N^{a}\alpha_{m}N^{b}\Big{]}d\Omega,$$

$$\mathbf{L}^{ab}=\sum_{e\in\mathbf{I}_{e}}\int_{\Omega^{e}}\Big{[}\frac{\tilde{\alpha}_{f}\tau_{M}}{\rho}\nabla N^{a}\cdot\nabla N^{b}\Big{]}d\Omega,$$

$$\begin{array}{r}\mathbf{F}=\nabla \mathbf{d}+\mathbf{I},\mathbf{C}={\mathbf{F}}^T\mathbf{F},\mathbf{E}=\frac{1}{2}(\mathbf{C}-\mathbf{I}),\end{array}$$

$$\mathbf{S}=\lambda \text{tr}(\mathbf{E})\mathbf{I}+2\mu \mathbf{E},$$

$$\begin{array}{r}{\rho }_s\frac{\partial \mathbf{u}}{\partial t}={\rho }_s\mathbf{f}+\nabla \cdot {\sigma }_s\end{array}\text{in}{\Omega }_s,$$

$${\int }_{{\Omega }_s^0}\left[{\rho }_s^0\mathbf{w}(\dot{\mathbf{u}}-\mathbf{f})+\nabla \mathbf{w}:\mathbf{P}\right]d\Omega =0,$$

$${\mathbf{R}}_m^a(\dot{\mathbf{u}},\mathbf{d})={\int }_{{\Omega }_s^0}\left[{\rho }_s^0{N}^a(\dot{\mathbf{u}}-\mathbf{f})+\mathbf{F}\mathbf{S}\nabla N^a\right]d\Omega .$$

$${\mathbf{d}}^{a,n+1}={\mathbf{d}}^{a,n}+{\mathbf{u}}^{a,n+1}\Delta t+\frac{0.5\gamma -\beta }{\gamma -1}{\dot{\mathbf{u}}}^{a,n+1}\Delta t^2,$$

$${\mathbf{d}}^{a,n+{\alpha }_f}=(1-{\alpha }_f){\mathbf{d}}^{a,n}+{\alpha }_f{\mathbf{d}}^{a,n+1}.$$

$${\mathbf{K}}_s\Delta \mathbf{d}=-{\mathbf{R}}_m({\dot{\mathbf{u}}}^{n+{\alpha }_m},{\mathbf{d}}^{n+{\alpha }_f}),$$

$$\mathbf{K}_{s}^{ab}=\int_{\Omega_{s}^{0}}\Big{[}\rho_{s}^{0}\alpha_{m}N^{a}N^{b}\mathbf{I}+\hat{\alpha}_{f}\big{(}\mathbf{S}\nabla N^{a}\big{)}\cdot\nabla N^{b}\mathbf{I}\\ +\lambda\,\hat{\alpha}_{f}(\mathbf{F}\nabla N^{a})\otimes(\mathbf{F}\nabla N^{b})+\mu\,\hat{\alpha}_{f}(\mathbf{F}\nabla N^{b})\otimes(\mathbf{F}\nabla N^{a})\\ +\mu\,\hat{\alpha}_{f}\mathbf{F}\mathbf{F}^{T}\nabla N^{a}\cdot\nabla N^{b}\Big{]}d\Omega,$$

$${\mathbf{d}}^{a,n+1}\leftarrow {\mathbf{d}}^{a,n+1}+\beta \Delta t^2\Delta {\mathbf{d}}^a,\forall a\in {\mathbf{I}}_a.$$

$$D(\mathbf{K})=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\Big{[}|K_{ii}|/\Big{(}\sum_{\begin{subarray}{c}j=1\\ i\neq j\end{subarray}}^{N_{c}}|K_{ij}|\Big{)}\Big{]},$$

$$S(\mathbf{A})=\frac{|{\tilde{\mathbf{A}}}_{sym}|-|{\tilde{\mathbf{A}}}_{skew}|}{|{\tilde{\mathbf{A}}}_{sym}|+|{\tilde{\mathbf{A}}}_{skew}|},$$

Full text

††journal: Computational Mechanics

1 Introduction

Cardiovascular simulations are increasingly used in clinical decision making, surgical planning, and medical device design. In this context, numerous modeling approaches have been proposed ranging from lumped parameter descriptions of the circulatory system to fully three-dimensional patient-specific representations. Patient-specific models are generated through a pipeline progressing from segmentation of medical image data, to branch lofting, boolean union, application of physiologic boundary conditions tuned to match patient data, and hemodynamics simulation. In diseased vessels, e.g., characterized by localized stenosis or aneurysms, computational fluid dynamics (CFD) has been widely used to diagnose important clinical indicators, such as pressure drop or flow reduction Taylor2013 . Measures of shear stress on the vessel lumen have also been correlated with the risk of endothelial damage and thrombus formation Segupta2012 ; Nesbitt2009 . These quantities are determined by discretization in space and time of the incompressible Navier Stokes equations. Multiscale models have been developed to simulate the local flow field in thee-dimensional patient-specific anatomies, while accounting for the presence of the peripheral circulation through closed-loop circuit models providing time dependent boundary conditions Lagana2002 ; Corsini2011 ; Marsden2015 ; Esmaily2012 . In addition, several approaches for fluid-structure interaction (FSI) have been suggested to account for vessel wall deformability Figueroa2006 ; Long2013 ; Ramachandra2017 . Recently, hemodynamic models have been used in the solution of complex problems in optimization Marsden2008 ; Marsden2014 and uncertainty quantification Sankaran2016 ; Schiavazzi2016 ; Schiavazzi2017a ; Schiavazzi2017 ; Tran2017 .

Efforts to improve realism in numerical simulations, however, often lead to an increase in the computational cost. Implementation of the coupled-momentum method (CMM) for FSI increases the simulation run time by roughly twice compared to rigid wall assumptions Figueroa2006 , and the cost can be substantially higher for Arbitrary Lagrangian and Eulerian (ALE) FSI. On top of this, optimization and uncertainty quantification studies often require the solution of a large number of simulations to obtain converged solutions. These requirements all point to a pressing need to reduce the computational cost to enable future integration of these tools in the clinical setting.

Preconditioned iterative approaches are widely used to solve linear systems, $\mathbf{Ay}=\mathbf{b}$, resulting from discretizations using variational multiscale finite element methods. However, few studies in the literature have examined in detail how linear solver performance depends on the properties of the coefficient matrix, to provide a concrete grounding for the choice and development of more efficient solvers. In addition, even fewer studies have carried out this analysis in the context of computational hemodynamics, i.e., the specific geometries, boundary conditions, mesh and material properties used to create numerical approximations of blood flow in rigid and deformable vessels. In this study, we investigate reductions in computational cost achievable through solving the discretized linear system efficiently, as this cost is well known to dominate the execution time. The objective of this study is to perform a systematic comparative analysis of existing linear solver strategies, linking their performance with the distributed coefficient matrix characteristics of cardiovascular modeling.

Krylov subspace based iterative solvers are typically preferred for the solution of large linear systems from CFD, due to their superior scalability and memory requirements compared to direct methods Benzi2002 . Popular Krylov subspace iterative solvers include the conjugate gradient method (CG) for symmetric positive definite (SPD) coefficient matrices, and the generalized minimum residual method (GMRES) or the bi-conjugate gradient stabilized method (BICGS) in the non-symmetric case. Alternatively, a recently proposed bi-partitioned linear solver (BIPN) Esmaily2015 leverages the block-structure in the coefficient matrix, separating contributions from the Navier-Stokes momentum and continuity equations. In BIPN, the coefficient matrix $\mathbf{A}_{f}$ arising from finite element spatial discretizations and the time discretization for the Navier-Stokes equations consists of four blocks,

[TABLE]

in which $\mathbf{K}$ and $\mathbf{G}$ stem from the momentum equation, $\mathbf{D}$ and $\mathbf{L}$ stem from the continuity equation and stabilization. BIPN solves the matrix block $\mathbf{K}$ using GMRES, while the rest is transformed to a Schur complement form as $\mathbf{L}-\mathbf{D}\mathbf{K}^{-1}\mathbf{G}$, approximated by a SPD matrix, $\mathbf{L}_{*}+\mathbf{G}_{*}^{T}\mathbf{G}_{*}$, in which the star subscript indicates the symmetric Jacobi scaling with diagonals of $\mathbf{K}$ and $\mathbf{L}$. The Schur complement form is solved with CG, and the solution time for CG takes more than 90$\%$ of the total computing time in benchmark testingEsmaily2015 .

It is also well known that preconditioning plays a key role in accelerating the convergence of Krylov subspace methods Trefethen ; Saad2003 by transforming the original linear system $\mathbf{A}\mathbf{y}=\mathbf{b}$ to $\mathbf{M}^{-1}\mathbf{A}\mathbf{y}=\mathbf{M}^{-1}\mathbf{b}$ (left preconditioning), $\mathbf{A}\mathbf{M}^{-1}\mathbf{z}=\mathbf{b},\mathbf{y}=\mathbf{M}^{-1}\mathbf{z}$ (right preconditioning), or $\mathbf{M}_{1}^{-1}\mathbf{A}\mathbf{M}_{2}^{-1}\mathbf{y}=\mathbf{M}_{1}^{-1}\mathbf{b}$, $\mathbf{x}=\mathbf{M}_{2}^{-1}\mathbf{y}$ (left and right preconditioning). In many cases, $\mathbf{M}$ is constructed in a way that $\mathbf{M}^{-1}$ approximates $\mathbf{A}^{-1}$. In general, an ideal preconditioner should be relatively cheap to apply and effective to reduce the overall solution time. In its simplest form, left, right or left-right Jacobi (diagonal) preconditioners are effective in shrinking the eigenvalue spectrum of diagonally dominant matrices. Preconditioners based instead on incomplete factorization (ILU) provide an approximate decomposition in the form $\mathbf{M}=\bar{\mathbf{L}}\bar{\mathbf{U}}$, where the sparsity pattern of $\mathbf{A}$ is preserved in the factors. ILUT preconditioners are slightly more general approaches allowing for adjustable inclusion of fill-ins, but require the user to specify an additional threshold parameter. We note that the efficiency of an ILU preconditioner results from a trade off between fewer Krylov iterations needed for convergence and the cost of incomplete factorization Saad2003 . Application-specific preconditioners have also been proposed in cardiovascular hemodynamics to improve performance when the model outlets are coupled through a resistance boundary condition, an RCR circuit, or for more general multi-domain configurations. In what follows, we refer to an in-house implementation of this class of preconditioners as resistance-based preconditioner (RPC) Esmaily2013a ; Esmaily2015 . Additional preconditioning techniques for cardiovascular simulations with FSI are suggested in Manguoglu2010 ; Deparis2016 . Finally, algebraic multigrid preconditioners have also received significant recent interest Santos2004 .

Despite the availability of open-source implementations of iterative solvers and preconditioners, few studies have systematically compared the performance of these solvers for cardiovascuar models with rigid and deformable vessels. In addition, a thorough understanding of the factors affecting the performance of iterative linear solvers (e.g., diagonal dominance, condition number, sparsity pattern, symmetry, and positive-definiteness) is an important prerequisite for optimal choice of solver and for the development of new algorithms with improved performance.

In the current study, we first compare the performance of various iterative linear solvers and preconditioners for an idealized flow through a cylindrical vessel with a resistance outflow boundary condition. We then test our findings using three representative patient-specific cardiovascular models. The Trilinos software library Trilinos , developed at Sandia National Laboratory is coupled with the SimVascular svFSI open source finite element code to provide linear solvers such as GMRES and BICGS, as well as a variety of preconditioners: diagonal (Diag), block-diagonal (BlockD), incomplete LU (ILU), thresholded incomplete LU (ILUT), incomplete Cholesky (IC), and algebraic multigrid (ML). We use Kylov linear solvers and Diag, BlockD, ILU, ILUT preconditioners from the AztecOO package, and IC is implemented via the IFPACK package. The block-diagonal preconditioner scales the block matrix via the Trilinos Epetra Vbr class. The incomplete factorization methods use Additive Schwarz domain decomposition for the parallelization. The Trilinos ML package for algebraic multigrid is implemented on the AztecOO package. For detailed information on parallelization and preconditioning options, we refer readers to the Trilinos project Trilinos . BIPN and RPC are integrated and implemented directly in our flow solver with source code available through the SimVascular open source project (www.simvascular.org) Updergrove2016 .

This paper is organized as follows. In section 2, we review the formulation of the coefficient matrices resulting from finite element weak forms in fluid and solid mechanics, before discussing the performance of various linear solvers and preconditioners on a simple pipe benchmark in section 3. In section 4 we report the results of strong and weak scaling for BIPN with RPC and BICG with ILU, while in section 5 we examine the properties of the coefficient matrix. The effect of preconditioning on these properties is reported in section 6. In section 7 we compare performance of linear solvers in patient-specific models. We draw conclusions and discuss future work in section 8.

2 Linear systems in cardiovascular simulation

We begin by introducing the space-time discretization of the equations governing fluid and solid mechanics following an Arbitrary-Lagrangian-Eulerian (ALE) description of the interaction between fluid and structure Huges1981 ; Bazilevs2009 ; Esmaily2015 ; Vedula2017 . These equations are discretized with a variational multiscale finite element method, and are provided in the svFSI solver of the SimVascular open source project Updergrove2016 .

2.1 Linear system for fluid mechanics

Consider a domain $\Omega_{f}\in\mathbb{R}^{3}$, occupied by a Newtonian fluid whose evolution in space and time is modeled through the incompressible Navier-Stokes equations in ALE form,

[TABLE]

where $\rho,\mathbf{u}=\mathbf{u}(\mathbf{x},t)$, and $\mathbf{f}$ are fluid density, velocity vector, and body force, respectively. The fluid stress tensor is $\mathbf{\sigma}_{f}=-p\,{\mathbf{I}}+\mu\,(\nabla{\mathbf{u}}+\nabla{\mathbf{u}}^{T})=-p\,{\mathbf{I}}+\mu\,\nabla^{s}{\mathbf{u}}$, $\mu$ is the kinematic viscosity, $p=p(\mathbf{x},t)$ the pressure, $\mathbf{v}=\mathbf{u-\hat{u}}$ is the fluid velocity relative to the velocity of the domain boundary $\mathbf{\hat{u}}$. Variables are interpolated in space at time $t^{n}$ as

[TABLE]

[TABLE]

in which $\mathbf{I}_{a},N^{a},\mathbf{w}^{a},q^{a},\mathbf{u}^{a}$, and $p^{a}$ are the nodal connectivity set, interpolation functions at node $a$, test function weights, velocity and pressure at node $a$, respectively. In this study, we employ P1-P1 type (linear and continuous) spatial approximations of the fluid velocity and pressure. We consider a stabilized finite element discretization based on the variational multiscale method Bazilevs2008 ; Esmaily2015 , leading to the weak form of the Navier-Stokes momentum and continuity residuals

[TABLE]

[TABLE]

in which $\mathbf{R}_{m}^{a}$ and $R_{c}^{a}$ are momentum and continuity residuals at node $a$, and $\mathbf{h}$ is the surface traction on the Neumann boundary $\Gamma^{h}$. The stabilization parameters are defined as

[TABLE]

where $C_{I}$ is a constant set to 3, $\Delta t$ is the time step size, and $\xi$ represents natural coordinates. Integration in time is performed using the unconditionally stable, second order accurate generalized-$\alpha$ method jansen2000generalized , consisting of four steps: predictor, initiator, Newton-Raphson, and corrector. Initial values for accelerations, velocities and pressures at time $t^{n+1}$ are set in the prediction step as

[TABLE]

where $\gamma=0.5+\alpha_{m}-\alpha_{f}$, $\alpha_{m}=1/(1+\rho_{\infty})$, and $\alpha_{f}=(3-\rho_{\infty})/(2+2\rho_{\infty})$ are the generalized-$\alpha$ method coefficients, while $\rho_{\infty}$ is the spectral radius set to $\rho_{\infty}=0.2$ in this study. In the initiator step, accelerations and velocities are computed at an intermediate stage $n+\alpha_{m}$ and $n+\alpha_{f}$,

[TABLE]

A Newton-Raphson iteration is performed based on Equations (5) and (6), using $\dot{\mathbf{u}}^{n+\alpha_{m}}$, ${\mathbf{u}}^{n+\alpha_{f}}$, ${p}^{n+1}$ from (9) by solving a linear system of the form

[TABLE]

where the blocks $\mathbf{K}$, $\mathbf{G}$, $\mathbf{D}$, and $\mathbf{L}$ partition the tangent coefficient matrix with blocks for nodes $a$ and $b$ equal to

[TABLE]

We re-write this linear system in matrix form as

[TABLE]

where

[TABLE]

with blocks $\mathbf{K}$, $\mathbf{G}$, $\mathbf{D}$, $\mathbf{L}$ of size ($3N_{nd}\times 3N_{nd}$), $(N_{nd}\times 3N_{nd}$), ($3N_{nd}\times N_{nd}$), ($N_{nd}\times N_{nd}$), respectively. Here $N_{nd}$ is the total number of nodes, while $\Delta\mathbf{u}\in\mathbb{R}^{3\,N_{nd}}$ and $\Delta\mathbf{p}\in\mathbb{R}^{N_{nd}}$ contain nodal velocities and pressure increments. We note that the major focus of our study is on solving the linear system in equation (12). Once the momentum and continuity residual norms drop below a given tolerance, the unknowns at the next time step are determined through the corrections

[TABLE]

$\forall\,a\in\mathbf{I}_{a}$. Finally, detailed expressions for each block of the coefficient matrix from Eq. (5), Eq. (6), Eq. (11), and Eq. (14) are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

in which $\tilde{\alpha}_{f}=\gamma\,\Delta t\,\alpha_{f}$ and $\mathbf{I}_{e}$ is the list of elements containing nodes $a$ and $b$.

First we observe that the matrix $\mathbf{K}$ is diagonally dominant. Except for entries related to the stabilization terms, which are typically small, the most significant off-diagonal contribution is provided by the viscous term, which is also typically smaller than the acceleration and advection terms in cardiovascular flows. Second, the small magnitude of the stabilization terms suggests that $\mathbf{G}$ is similar to $-\mathbf{D}^{T}$. We also observe that the matrices $\mathbf{K}$ and $\mathbf{A}_{f}$ are non-symmetric, while $\mathbf{L}$ is symmetric and singular since it has an identical structure to the matrices arising from the discretization of generalized Laplace operators. We also note that $\mathbf{L}$ is characterized by small entries compared to the other blocks, since it only consists of stabilization terms.

2.2 Linear system for solid mechanics

In the solid domain, we start by introducing measures of deformation induced by a displacement field $\mathbf{d}=\mathbf{x}-\mathbf{X}$, i.e., the difference between the current and material configurations $\mathbf{x}\in\mathbb{R}^{3}$ and $\mathbf{X}\in\mathbb{R}^{3}$, respectively

[TABLE]

where $\mathbf{F}$, $\mathbf{C}$, $\mathbf{E}$, represent the deformation gradient, the Cauchy-Green deformation tensor and the Green strain tensor. The Jacobian is also defined as $J=\text{det}(\mathbf{F})$. We relate the second Piola-Kirchhoff stress tensor $\mathbf{S}$ with the Green strain tensor $\mathbf{E}$ through the Saint Venant-Kirchhoff hyperelastic constitutive model

[TABLE]

where $\lambda=\frac{\nu E_{s}}{(1+\nu)(1-2\nu)},\,\mu=\frac{E_{s}}{2(1+\nu)}$, $E_{s}$ and $\nu$ represent the Young’s modulus and Poisson’s ratio, respectively. The equilibrium equation is

[TABLE]

where $\rho_{s}$ and $\sigma_{s}$ denote the density and solid stress tensor, respectively. This leads to the weak form

[TABLE]

where $\mathbf{P}=\mathbf{F}\mathbf{S}$ is the first Piola-Kirchhoff stress tensor, $\mathbf{w}$ is a virtual displacement and $\Omega^{0}_{s}$ is the solid domain in the reference configuration. Discretization of (22) leads to the residual

[TABLE]

Using the generalized-$\alpha$ method, the displacements at time $t^{n+1}$ are predicted as

[TABLE]

in which $\beta=\frac{1}{4}(1+\alpha_{f}-\alpha_{m})^{2}$. In the initiator step, the intermediate displacements are provided by

[TABLE]

Solving (23) with the Newton-Raphson method, we obtain the linear system

[TABLE]

where $\mathbf{K}_{s}^{ab}\approx\frac{\partial\mathbf{R}_{m}^{a}}{\partial\Delta d^{b}}$, with tangent stiffness matrix

[TABLE]

and $\hat{\alpha}_{f}=\alpha_{f}\,\beta\,\Delta t^{2}$. At this point, $\mathbf{d}$ is corrected using

[TABLE]

Again, we note that the focus of this work is on solving the linear system in equation (26) together with solving equation (12) when including FSI. Finally, we observe that most terms in (27) (i.e., the third, fourth, and fifth terms) contribute to the off-diagonals of $\mathbf{K}_{s}$ and that the matrix $\mathbf{K}_{s}$ is symmetric.

3 Linear solver performance on pipe flow benchmark

All tests discussed in this study are performed using the svFSI finite element solver from the SimVascular open source project, leveraging the message passing interface (MPI) library and optimized to run efficiently on large computational clusters Esmaily2015a . We note this implementation assigns fluid and solid elements to separate cores and assumes unique interface nodes, i.e., common nodes on the fluid and solid mesh are expected to match. An FSI mesh with matching interface nodes is generated through the freely-available Meshmixer software, with details reported in Vedula2017 .

We consider a simple pipe benchmark whose size and boundary conditions are chosen to represent the ascending aorta of a healthy subject, having $4$ cm diameter, $30$ cm length and $0.2$ cm thickness, assuming a radius/thickness ratio of 10%. The inflow is steady with a parabolic velocity profile, having a mean flow rate of $Q=83$ mL/s ($5$ L/min). A resistance boundary condition is used at the outlet equal to $R=1600$ g/cm4/s, which produces a mean pressure of approximately $\Delta P=100$ mmHg, typical of the systemic circulation of a healthy subject. Simulations are performed with rigid and deformable walls (see Figure 1), using 1,002,436 tetrahedral elements for the fluid domain and 192,668 tetrahedral elements for the wall, generated in SimVascular with the TetGen mesh generator plug-in Updergrove2016 ; Hongzhi2017 . We measure the wall clock time on the XSEDE Comet cluster for simulations consisting of 10 time steps of 1 millisecond each, using 38 and 48 cores for rigid and deformable simulations, respectively. The XSEDE Comet cluster has 1,944 compute nodes with the Intel Xeon E5-2680v3 cores, 24 cores/node, 2.5GHz clock speed, 960 GFlop/s flop speed, and 120GB/s memory bandwidth. For more information about the Comet cluster, please refer to the XSEDE user portal. We use the restarting scheme for GMRES. We set the restart number of 200 which showed superior performance against smaller restart numbers (See Appendix A). We set the ILUT drop tolerance to ${10^{-2}}$ and the fill-in level to 2. In our test, changing the drop tolerance to ${10^{-4}}$ and ${10^{-6}}$ did not significantly change the performance of the linear solver reported here. For the multigrid preconditioner we selected a maximum level of four, a Gauss-Seidel smoother, and the symmetric Gauss-Seidel for the subsmoother. We confirmed that this setting was superior to other choices for smoother and subsmoother (See Appendix B). Finally, as the node ordering affects the amount of fill-in produced by an ILU decomposition, we apply Reverse-Cuthill-McKee (RCM) reordering prior to the incomplete factorization. RCM has been shown to be effective among many reordering schemes in the solution of non-symmetric sparse linear systems (see, e.g., Benzi1999 ). From our testing, ILUT with RCM reordering provided superior performance against ILUT without reordering and ILUT with METIS reordering (See Appendix C).

3.1 Rigid wall benchmark

We plot the linear solver performance measured by wall clock time in Figure 2. In Table 2, we report the number of iterations and the portion of total compute time consumed by solving the linear system. Three iterative solver tolerances, $\epsilon=10^{-3},10^{-6},10^{-9}$, are tested and compared. Table 1 shows the effect of tolerance on the velocities at $t=10$ ms, suggesting the velocity error norm is of the same order of the selected tolerance.

Figure 2 shows that incomplete factorization preconditioners are fast and exhibit robust performance across all tolerances, irrespective of the underlying iterative linear solver.

Despite similar performance for ILU, ILUT and IC preconditioners across all cases, the slightly worse performance of ILUT with respect to ILU seems to suggest that constructing a more accurate factor may lead to smaller run times than the savings in the factorization cost acheived by dropping additional fill-ins. GMRES with diagonal preconditioners (either diagonal or block-diagonal) appears to be significantly slower than other schemes, particularly as the tolerance $\epsilon$ becomes smaller. This degrading performance of standard GMRES for cardiovascular modeling is consistent with previous studies, showing that resistance boundary conditions are responsible for an increase in the condition number Esmaily2013a ; Esmaily2015 . While BICGS seems to perform better than GMRES with diagonal preconditioners, its performance becomes increasingly unstable with smaller tolerance $\epsilon$. Furthermore, while algebraic multigrid preconditioners appear to be superior to diagonal preconditioning, they are inferior to BIPN or GMRES/BICGS with ILU. Finally, the performance of BIPN with RPC preconditioning is comparable to ILU for small tolerances (i.e., $\epsilon=10^{-3}$) but significantly degrades as the tolerance value decreases.

From Table 2, we see that the time consumed by the linear system solvers constitutes the majority of the total compute time. In BIPN, the compute time of GMRES is significantly smaller than the compute time of CG. As the tolerance value decreases, the relative percentage of compute time for GMRES solve becomes larger. All Trilinos preconditioners show larger iteration numbers with decreasing tolerance. The relatively small percentage of linear solver compute time against the total compute time with IC and ML implies that building such preconditioners is expensive. This suggests that storing and reusing a preconditioner for several time steps could increase efficiency.

3.2 Deformable wall benchmark (FSI)

We illustrate the compute times for the deformable wall case in Figure 3 and summarize the number of iterations and percentage compute time of the linear solvers in Table 3. For FSI simulations with a tolerance $\epsilon=10^{-3}$, BIPN with RPC shows more than an 8 fold increase in compute time compared to the rigid wall case, which also suggests the limitations of this approach for the deformable wall case. The increase of iteration number for GMRES as tolerance values decrease is notably higher than in the rigid wall case, and the percentage of compute time for the GMRES solve in BIPN is significantly higher, about 80$\%$ of the total compute time while the CG part continues to make up a smaller percentage. This suggests directions for future improvement of BIPN in the GMRES part for FSI simulations. Conversely, incomplete factorization preconditioners (both for GMRES and BICGS), exhibit good performance across all tolerances and a limited increase in compute time compared to the rigid case. Diagonal preconditioners show a comparable performance to incomplete factorization schemes at large tolerances, but the performance degrades for smaller tolerances. Among all algorithms implemented in Trilinos, the algebraic multigrid preconditioners is the slowest, while BICG with ILU appears to be the best solution scheme overall.

4 Parallel scalability

Parallel scalability is investigated for two preconditioned linear solvers, BIPN-RPC and BICG-ILU, in terms of speedup (see, e.g., Chen2006 ), defined as the computing speed of multiple cores compared to a single core calculation, i.e., $S_{p}=T_{1}/T_{p}$, where $T_{p}$ is the compute time on $N_{p}$ cores. The ideal strong scalability performance between $S_{p}$ and $N_{p}$ is linear, however, sublinear scaling is expected due to communication cost.

4.1 Strong scaling

In this section, we monitor the compute time for a model with a fixed number of degrees of freedom, while progressively increasing the number of cores. We first test the strong scalability by varying the number of cores, testing 1, 2, 4, 8, 24, 48, and 96 cores for an $\approx$1 million element mesh (Mesh1 in Table 4) (Figure 4). We chose the number of cores as multiples of 24, to use all cores in any given node. We note, however, that one should use only between 2/3 and 3/4 of the total number of cores in a given node since the local memory bandwidth is often a bottleneck resulting in higher overhead and reduced speed improvement. In the rigid wall model, BIPN-RPC and BICG-ILU show similar performance across all core numbers as shown in Figure 2. The parallel speedup shows that both methods scale well up to 24 cores, leading to about 40,000 elements per core, while their parallel performance is reduced when running on more than 48 cores. In Figure  2, the zig-zag behavior of BIPN-RPC (FSI) reveals excessive inter-core communications and memory references. In the FSI problem, the total compute time of BIPN-RPC is larger, and the speed up is worse than BICG-ILU. BICG-ILU shows consistently good scalability for both rigid and FSI models whereas the scalability of BIPN-RPC degrades significantly in the FSI case.

We conduct an additional scaling study on a refined mesh with 8M elements (Mesh5 in Table 4). We use 24, 48, 96, 192, and 384 cores and we use the $N_{p}=24$ case as a reference for $S_{p}$. As shown in Figure 5, the BICG-ILU speedup scales almost linearly up to 192 cores, i.e., $\sim$40,000 elements per core. In the 8M mesh, the compute time of BIPN-RPC becomes significantly larger due to a poor weak scalability, as discussed in the next section.

4.2 Weak scaling

In this section, we solve models of increasing size and report the simulation time by keeping approximately the same number of elements per core. The number of nodes, elements and non-zeros in the coefficient matrix is summarized for each mesh in Table 4. We use 20, 40, 80, 160, 320 cores for meshes 1 to 5, resulting in 27,000 elements and about 4,500 nodes per a core.

BICG-ILU shows increased compute time but does not show any significant changes in scaling when increasing the number of cores. For the rigid model, BIPN-RPC shows poor scalability as the number of cores exceeds 160. This explains the loss of superior performance of BIPN-RPC against BICG-ILU in the strong scalability study with 8M mesh. Similar to the rigid case, BIPN-RPC shows performance loss after 160 cores for the FSI model.

5 Characteristics of unpreconditioned matrices

In this section, we investigate the properties of the coefficient matrices discussed in section 2.1 to better understand the performance results. With reference to the pipe benchmark, we visualize the matrix sparsity pattern and investigate its properties including bandwidth, symmetry, diagonal dominance and spectrum. We also investigate the characteristics of both global and local matrices. The global matrix contains all mesh nodes from all cores, while a local matrix contains only a subset of the nodes in the global matrix assigned to a single core upon partitioning. We report single-core, local matrix characteristics with $\text{N}_{lnd}\sim 5000$ nodes to represent the typical case of distributed discretizations consisting of $\sim$ 25000 tetrahedral elements per core. This way we focus on detailed local information in a region of specific interest (e.g. resistance boundary), and also calculate properties of the matrix such as eigenvalues and condition numbers in a cost-effective way. We discuss properties for two groups of coefficient matrices, i.e., matrices associated with fluid flow and matrices associated with solid mechanics.

5.1 Matrix properties for fluid flows in rigid vessels

Sparsity pattern. The structure of $\mathbf{A}_{f}$ is determined by element nodal connectivity, i.e., the non-zero columns for node $a$ are from nodes belonging to the element star associated to $a$. The star of elements associated to a given node $a$ is the set of all elements connected to $a$. An example of global sparsity pattern for the pipe benchmark unstructured mesh is reported in Figure 7. The node ordering starts from the outer surface and goes to the interior of the pipe, as seen from the reduced connectivity between the upper-left block and the remaining inner nodes. The density is less than 1 percent, and sparsity is more than 99.99 percent (Table 4).

A Reverse-Cuthill-McKee (RCM) bandwidth minimizing permutation of the same raw matrix shows a banded sparse structure illustrated in Figure 7. We report quantitative estimates for bandwidth and number of non-zeros, showing that bandwidth is approximately 1000 with a maximum of about 1500, and most nodes are connected to 15-16 other nodes.

Diagonal dominance. A closer look at the magnitudes of entries in the global matrix $\mathbf{A}_{f}$ reveals a clear block-structure (Figure 8). The 4-by-4 block structure corresponds to the matrix blocks in equation (13). This shows that diagonal entries are larger than the off-diagonals in each block of $\mathbf{K}$, suggesting that magnitudes of entries in the blocks of $\mathbf{G},\mathbf{D},\mathbf{L}$ are small compared to diagonals in $\mathbf{K}$. This relates to the dominant contribution of the acceleration and advection terms, over the stabilization and viscous terms in (15) (see, e.g., Esmaily2015 ). The matrix $\mathbf{A}_{f}$ is, however, not diagonally dominant (i.e. $|A_{ii}|<\sum_{j=1}|A_{ij}|$). To show this, we quantified the relative magnitudes of off-diagonal and diagonal entries (Figure 9). Specifically, we counted the number of elements that have a certain percentage of absolute magnitudes compared to the diagonal values, showing that most off-diagonal values are less than $20\%$ of the associated diagonal. To report quantitative estimates of diagonal dominance, we measure the mean of the relative magnitude of the sum of off-diagonal values to the diagonal,

[TABLE]

in which $N_{r}$ and $N_{c}$ is the number of rows and columns in $\mathbf{K}$ respectively. $D$ increases with the diagonal dominance. D($\mathbf{A}_{f}$) for the global matrix is 0.514 and D($\mathbf{K}$) for the global matrix is 0.678.

Symmetry. We use an index, $S$, to quantify how close a matrix is to symmetric. We first obtain off-diagonal elements of $\mathbf{A}$ by subtracting its diagonal as $\tilde{\mathbf{A}}=(\mathbf{A}-\text{diag}(\mathbf{A}))$. We then decompose the $\tilde{\mathbf{A}}$ into the symmetric part, $\tilde{\mathbf{A}}_{sym}$=$(\tilde{\mathbf{A}}+\tilde{\mathbf{A}}^{T})/2$, and the skew-symmetric part, $\tilde{\mathbf{A}}_{skew}=(\tilde{\mathbf{A}}-\tilde{\mathbf{A}}^{T})/2$. The index $S$ is defined as

[TABLE]

in which we use the 2-norm for $|\mathbf{A}|$. The index equals $-1$ for a perfectly skew-symmetric matrix and $1$ for a perfectly symmetric matrix. As shown in Section 2.1, the matrix is nonsymmetric in the $\mathbf{K}$, $\mathbf{G}$, $\mathbf{D}$ blocks due to stabilization and convective terms, with $S(\mathbf{A}_{f})$ equal to 0.9859. That is, $S(\mathbf{A}_{f})$ is a nearly symmetric matrix in the analyzed regime (ideal aortic flow). Finally, $\mathbf{L}$ is a symmetric and semi-positive definite matrix with $S(\mathbf{L})=1$.

Eigenvalues. Spectral properties are widely used to characterize the convergence and robustness of iterative solvers. It is well known, for example, that the rate of convergence of CG depends on the spectral radius of a left-hand-side SPD matrix. Despite eigenvalues clustered around 1 leading to rapid convergence of iterative solvers for well-conditioned SPD matrices, the eigenvalues may not be solely responsible for the convergence rate of these solvers and other matrix characteristics may play a role Saad2003 . Calculation of all eigenvalues ($\lambda_{i}$) of the global matrix for a typical cardiovascular model with order 1 million mesh elements is prohibitively expensive. In this paper we therefore report the spectrum of local matrices instead of the global matrix. For a smaller size system, we also demonstrate that the distribution of eigenvalues from local matrices is a good approximation to the distribution of eigenvalue of the global matrix (See appendix D).

In Figure 10, we plotted the spectrum of local $\mathbf{A}_{f}$ matrices from the pipe benchmark with rigid walls. The eigenvalues of $\mathbf{A}_{f}$ are complex with small magnitudes of the imaginary part up to $O(10^{-8})$, while the magnitudes of the real part ranges from $O(10^{-9})$ to $O(10^{-1})$. In Figure 10, there are three distinct groups of eigenvalues with different ranges of the real part. The first group contains those with the real part less than $10^{-5}$. The second group contains those with real part ranging between $10^{-5}$ and $10^{-2}$. The spectrums of $\mathbf{K}$ and $\mathbf{L}$ in Figure 10 show that the third group, with eigenvalues larger than $10^{-5}$ in $\mathbf{A}_{f}$, is attributed to the $\mathbf{K}$ block, while the eigenvalues of $\mathbf{L}$, $\mathbf{G}$ and $\mathbf{D}$ are responsible for the group of smallest real eigenvalues in $\mathbf{A}_{f}$. We list several minimum and maximum eigenvalues of a local matrix without resistance boundary condition (BC) in Table 3. The maximum eigenvalues of $\mathbf{K}$ and $\mathbf{A}_{f}$ appear to be the same, suggesting that block $\mathbf{K}$ dominates the high portion of the spectrum, with the smallest eigenvalues provided by the blocks $\mathbf{L}$, $\mathbf{G}$ and $\mathbf{D}$. This suggests that the large condition number of $\mathbf{A}_{f}\sim O(10^{6})$, obtained from MATLAB condest, relates to the inhomogeneous eigenvalue spectrum observed in the momentum, continuity and coupling blocks. Additionally, the small condition number of $\mathbf{K}$ justifies the idea behind the BIPN approach, i.e., to solve the $\mathbf{K}$ block separately, while expressing the other blocks in Schur complement form Esmaily2015 . $\mathbf{L}$ is singular, and thus has zero eigenvalue and an extremely large condition number. Lastly, the resistance boundary condition is responsible for the few largest eigenvalues order of $O(10^{-1})$ in Figure 10, as we discuss in the next section.

Effects of the resistance boundary condition. A resistance boundary condition perturbs the condition number of the coefficient matrix $\mathbf{A}_{f}$, and may be responsible for a significant increase in the solution time for the tangent linear system Esmaily2013a . The boundary traction $\mathbf{h}$, is given, in this case by

[TABLE]

in which $P_{i}$ is the pressure at surface $i$, evaluated as $P_{i}=R_{i}\,Q_{i}$, i.e., proportional to the flow rate $Q_{i}$ across the surface

[TABLE]

through the prescribed resistance $R_{i}$. Thus, the contribution of the resistance boundary condition to the coefficient matrix is

[TABLE]

where $n^{bc}$ is the number of resistance boundaries, and $\tilde{R}_{i}=\gamma\,\Delta t\,R_{i}$. $\mathbf{K}^{bc}$ is finally added to the $\mathbf{K}$ sub-matrix, resulting in the coefficient matrix $\mathbf{\tilde{K}}=\mathbf{K}+\mathbf{K}^{bc}$. Generalization from an outlet resistance to a coupled lumped parameter network model is accomplished using a slightly more general expression for $\mathbf{K}^{bc}$, i.e.

[TABLE]

where the resistance matrix $M_{kl}$ is obtained by coupling pressures and flow rates at different outlets Esmaily2012 .

Addition of a resistance boundary condition alters the topology of the coefficient matrix due to the rank one contribution $\mathbf{S}_{i}\otimes\mathbf{S}_{i}$. This, in practice, couples all velocity degrees of freedom on a given outlet, significantly affecting the performance of matrix multiplication for the $\mathbf{K}$ block and the fill-in generated by its LU decomposition. Thus, the vector $\mathbf{S}_{i}$ is stored separately, to improve the efficiency of matrix multiplication and for RPC preconditioning. Figure 8 shows how the global matrix entries are affected by the presence of a resistance boundary condition, i.e., large magnitude components in the z-directional velocity block (lower-right block in $\mathbf{K}$) arise. To better highlight this effect, we show two local matrices with and without a resistance boundary condition in Figure 11. Addition of $\mathbf{K}^{bc}$ increases the contribution of off-diagonal entries moving the matrix further away from diagonal dominance.

The resistance BC perturbs the spectrum and increases the condition number of $\mathbf{A}_{f}$. The largest few eigenvalues in the spectrum of $\mathbf{A}_{f}$ and $\mathbf{K}$ in Figure 10 are calculated from local matrices in partitioned domains interfacing the resistance boundary. The change of spectral properties of $\mathbf{\tilde{K}}$ due to rank one contribution $\mathbf{K}^{bc}$ (see, e.g., Esmaily2015 ) is measured with maximum and minimum eigenvalues reported in Table 6. The maximum eigenvalue of $\mathbf{\tilde{K}}$ is significantly larger than $\mathbf{K}$, leading to a $\sim O(10)$-fold increase in the condition number of $\mathbf{\tilde{K}}$, thus increasing the spectral radius of the whole spectrum as shown in Figure 10. Additionally, our tests confirm that the first maximum eigenvalue $\mathbf{A}_{f}$ increases linearly with the assigned resistance. Thus, in a more general case, we expect several large eigenvalues to be added to $\mathbf{A}_{f}$ for models with multiple outlet resistances.

5.2 Matrix properties for fluid flow in deformable vessels

Sparsity pattern. The global sparsity pattern for the FSI mesh is illustrated in Figure 12, where nodes in the solid-fluid interface are ordered first, followed by nodes in the fluid region next to the interface, and nodes in the solid domain. As shown in Figure 7 when the connectivity matrix is reordered by RCM, the global sparsity pattern has a banded sparse matrix structure with the larger maximum bandwidth of $\approx 1750$ compared to the rigid case, while the number of non-zeros in a row is mostly clustered around 15 to 16, similar to the rigid case. Magnitudes of entries associated with solid nodes appear to be significantly larger than those in the fluid domain. For example, the magnitude of $\mathbf{K}_{s}$ is order one (Figure 13), whereas the magnitude of $\mathbf{K}$ is order $10^{-3}$ (Figure 11). In what follows, we focus on the local matrix $\mathbf{K}_{s}$ from Eq. (27) since the characteristics of $\mathbf{A}_{f}$ in FSI are similar to the rigid wall case.

Diagonal dominance. The metric introduced above to quantify diagonal dominance drops to $D(\mathbf{A}_{\text{FSI}})=0.4775$ for the global matrix $\mathbf{A}_{\text{FSI}}$, i.e., the additional FSI terms reduce the diagonal dominance of the system. The magnitudes of a local $\mathbf{K}_{s}$ matrix before and after diagonal scaling is compared, in Figure 13, to a local $\mathbf{K}$ from the fluid domain. The diagonally scaled block $\mathbf{K}_{s}$ qualitatively shows how the off-diagonals in $\mathbf{K}_{s}$ are larger than those in $\mathbf{K}$. The diagonal dominance metric for the local blocks $\mathbf{K}_{s}$ and $\mathbf{K}$ are found to be $D(\mathbf{K}_{s})=0.379$ and $D(\mathbf{K})=0.678$, respectively.

Symmetry and positive definiteness. The symmetry metric for matrix $\mathbf{K}_{s}$ is one, i.e., $S(\mathbf{K_{s}})=1$ and positive definite as expected from its properties and confirmed numerically.

Eigenvalues. We calculated and plotted the eigenvalue spectrum of local matrices from the FSI benchmark in Figure 14. All eigenvalues of $\mathbf{K}_{s}$ are real due to the symmetry of $\mathbf{K}_{s}$. As listed in Table 7, the magnitudes of eigenvalues in $\mathbf{K}_{s}$ is significantly larger than the magnitudes of eigenvalues from $\mathbf{A}_{f}$.

5.3 Discussion

Results from the previous sections suggest the following conclusions. First, the condition number of both the fluid and solid tangent matrices $\mathbf{A}_{f}$ and $\mathbf{K}_{s}$ appears to be large, and therefore preconditioning is necessary. Second, the fluid matrix $\mathbf{A}_{f}$ is more diagonally dominant than the solid matrix $\mathbf{K}_{s}$. This suggests that diagonal preconditioning is expected to be more effective for rigid wall simulations, but incomplete factorization preconditioners are expected to work better under fluid-structure interaction, consistent with the results obtained in the pipe benchmark. Third, resistance and coupled multidomain boundary conditions need a special treatment for preconditioning, due to their effect on the maximum eigenvalue and condition number.

6 Effect of preconditioning

In this section we investigate how application of various preconditioners affects the spectral properties of the coefficient matrix in both the rigid and deformable case by explicitly computing the preconditioned matrix $\mathbf{M}_{l}^{-1}\mathbf{A}\mathbf{M}_{r}^{-1}$, where $\mathbf{M}_{l}^{-1}$ is a left preconditioner and $\mathbf{M}_{r}^{-1}$ is a right preconditioner.

Consider a left and right Jacobi preconditioning for $\mathbf{A}_{f}$:

[TABLE]

resulting in the linear system $\mathbf{A}_{f}^{*}\mathbf{y}^{*}=-\mathbf{R}^{*}$. Spectral properties of the preconditioned matrix $\mathbf{A}_{f}^{*}$ and $\mathbf{K}_{s}^{*}$ are reported in Figure 15 and Table 8.

Table 8 shows how diagonal preconditioning is effective in improving the conditioning of $\mathbf{A}_{f}$, particularly for $\mathbf{K}$, without resistance BC. The condition number of $\mathbf{K}^{*}$ is reduced to $\sim$10, and only a few linear solver iterations are expected to be sufficient to substantially reduce the approximated residual. This again justifies the approach followed by the BIPN solver, where the linear system involving the momentum block $\mathbf{K}$ is solved separately, thus shifting the computational cost to the iterative solution of its Schur-complement. The condition number of the approximated Schur complement block $\mathbf{L}^{*}+{\mathbf{G}^{*}}^{T}{\mathbf{G}^{*}}$ is, in this example, equal to 120, consistent with previous findings in Esmaily2015 . This also justifies the fact that $\sim$80 percent of total compute time in BIPN is dedicated to the solution of the Schur complement linear system. Conversely, a symmetric Jacobi preconditioning does not seem to significantly reduce the condition number of the solid block $\mathbf{K}_{s}$. This is attributed to the presence of large off-diagonal values in $\mathbf{K}_{s}$ that are only marginally affected by diagonal preconditioning. As a result, the eigenvalues of $\mathbf{K}_{s}^{*}$ range from $O(10^{-3})$ to $O(1)$ while the eigenvalues of $\mathbf{A}_{f}^{*}$ range $O(10^{-2})$ to $O(1)$, shown in Figure 15, with the exception of a few eigenvalues from the resistance BC. This, in turn, justifies the superiority of incomplete factorization preconditioners for FSI simulations.

Application of a symmetric Jacobi preconditioner to a local coefficient with resistance boundary condition leads to the eigenvalues reported in Table 9. Despite a reduction by three orders of magnitude, the condition number is still one order of magnitude larger than in the case without resistance BC. Note that the maximum eigenvalue of the preconditioned matrix is one order of magnitude larger than the second maximum eigenvalue, consistent with previous observations. Thus, the RPC preconditioning proposed in Esmaily2013a seeks a preconditioning matrix $\mathbf{H}$ such that $\mathbf{H}\simeq{(\mathbf{K}^{total})}^{-1}$. The idea is to construct $\mathbf{H}$ by combining the diagonal components of $\mathbf{K}$ with the resistance contributions stored in $\mathbf{S}_{j}$ as

[TABLE]

where $\mathbf{K}^{d}$ is $\text{diag}(\mathbf{K})$. The preconditioned matrix $\mathbf{H\tilde{K}}$ has a small condition number (Table 7) and smaller off-diagonal entries (Figure 16).

7 Performance of linear solvers in patient-specific models

In this section we tested the performance of linear solvers on patient-specific cardiovascular models, in an effort to extrapolate the results obtained for the pipe benchmark to more realistic problems. We use three different models with a wide range of boundary conditions (i.e., resistance, RCR, coronary BC, closed-loop multidomain), with and without wall deformability and covering various patient-specific anatomies. All anatomic models were constructed from medical image data using SimVascular.

7.1 Pulmonary hypertension

The first model represents the left and right pulmonary arteries with associated branches and is used to investigate the effects of pulmonary hypertension (PH). The finite element mesh contains 3,223,072 tetrahedral elements to represent the pulmonary lumen, has rigid walls and 88 outlets with Windkessel (RCR) boundary conditions, prescribed through a coupled 0-D multi-domain approach Esmaily2012 ; Yang2018 (Figure 17). A pulsatile inflow waveform extracted from PC-MRI was imposed at the pulmonary artery inlet. This model is solved using a time step of 0.46 milliseconds and 120 cores ($\approx 25,000$ elements per core). The tolerance on the Newton-Raphson residual is set to $\epsilon=10^{-4}$.

We briefly report global matrix characteristics in the PH model. The diagonal dominance metric is $D(\mathbf{A}_{f})=0.5598$, and $D(\mathbf{K})=0.7397$. The metric value for $D(\mathbf{A}_{f})$ is similar to values from the pipe model with rigid walls. The matrix is nearly-symmetric, $S(\mathbf{A}_{f})=0.9903$.

Results in Figure 17(c) compare the performance of diagonal, block-diagonal and ILU preconditioning. BIPN with RPC shows the best performance, followed by ILU-BICG. This is expected due to the large number of resistance boundary conditions (i.e., 88) at the model outlets. Diagonal preconditioners with GMRES instead perform poorly, consistent with our observations in the pipe benchmark.

7.2 Coronary artery bypass graft

Second, we consider a model of coronary artery bypass graft surgery (see, e.g., Ramachandra2017 ; Tran2017 ) with 4,199,945 tetrahedral elements and rigid walls, coupled with a closed-loop 0D lumped parameter network (LPN), including coupled heart, coronary and systemic circulation models (Figure 18). Simulations were performed using 168 cores ($\sim 24,000$ elements per core), with a time step of 0.87 millisecond and a non linear iteration tolerance of $\epsilon=10^{-3}$.

The diagonal dominance metric for the global matrix in the CABG model is $D(\mathbf{A}_{f})=0.5200$, and $D(\mathbf{K})=0.6914$. The matrix in the CABG model is less diagonally dominant than the previous cylinder or pulmonary hypertension model. The matrix is also near symmetric, $S(\mathbf{A}_{f})=0.9938$.

As expected, due to the presence of coupled multi-domain boundary conditions Esmaily2013a , BIPN results in the best performance, followed closely by BICG with ILU, while GMRES with a diagonal preconditioner performs poorly. The relative performance of ILU against the diagonal preconditioner is better than seen in previous models. The smaller diagonal dominance metric in CABG model confirms superiority of ILU over the diagonal preconditioner.

7.3 Left coronary

Next, we tested performance of linear solvers in a left coronary model. The left coronary model was extracted from a full coronary artery model used in Tran et al.Tran2017 . The pulsatile flow waveform at the inlet of the left coronary branch was extracted from the full model simulation and imposed at the inlet of the model. The model has six outlets, each with applied open-loop coronary outlet boundary conditions Kim2009 . All resistance and compliance values as well as inflow pulsatile waveforms are determined to produce a normal physiologic response of the left coronary artery following our prior work. We ran simulations with rigid and deformable walls with a lumen mesh containing 486,066 tetrahedral elements and a vessel wall mesh with 206,369 tetrahedral elements, time step of 1 millisecond and tolerance of $\epsilon=10^{-4}$. We used 20 cores for the rigid wall simulation and 24 cores for the deformable wall simulation.

The diagonal dominance metric for the left coronary model with rigid wall is $D(\mathbf{A}_{f})=0.5238$, and $D(\mathbf{K})=0.6972$, which are similar to the CABG model. The matrix is near symmetric, $S(\mathbf{A}_{f})=0.9855$, however, this model is furthest from symmetric among all models considered.

For the ALE FSI model, the diagonal dominance is reduced as $D(\mathbf{A}_{FSI})=0.4323$, $D(\mathbf{K})=0.5320$. Note that the number of elements in a wall mesh is 40 percent of the fluid mesh, so the effect of adding solid mechanics in the linear system is more significant than the pipe case where only 20 percent of elements were in the in wall mesh. As a result, the symmetry metric is very close to 1, $S(\mathbf{A}_{FSI})=0.99998$, since $\mathbf{K}_{s}$ is symmetric.

As shown in Figure 20, performance test results are consistent with previous findings. RPC-BIPN is the fastest method for the rigid wall simulation. In FSI, the performance of BIPN is poor, while diagonal and ILU preconditioners with GMRES perform better.

7.4 Discussion

From the performance results in the patient-specific models, we find that RPC-BIPN is the fastest method for rigid wall simulations with many resistance BCs in agreement with the pipe model. ILU-BICG is only slightly slower than RPC-BIPN while the standard diagonal scaled GMRES fails. The performance degradation of RPC-BIPN in FSI models is consistent with the pipe model and suggests the need for future improvements of BIPN for ALE FSI.

8 Summary and conclusions

In this paper we study the performance of preconditioned iterative linear solvers for cardiovascular simulation in rigid and deformable walls. To this end, we implement various iterative linear solvers - GMRES, BICGS, and BIPN - and preconditioners - diagonal, block-diagonal, ILU, ILUT, ML, and RPC in a single flow solver. Standard iterative solvers and preconditioners are employed from the Trilinos library and compared against RPC-BIPN, implemented in our in house solver. Simulation wall clock time is measured and compared in a benchmark pipe flow with a resistance BC.

ILU preconditioned BICG provides the best overall performance in both rigid and deformable wall simulations. RPC-BIPN in the FSI simulation shows $\sim 8$ fold increase in compute time compared to the rigid wall case. Strong and weak scalings of ILU-BICG and RPC-BIPN are reported.

To better understand the observed performance, characteristics of the left-hand-side matrix in the linear system are examined. We report sparsity patterns, diagonal dominance, symmetry, eigenvalues and condition numbers in global and local matrices. Results show that the sparse matrix structure has a narrow banded structure after Reverse-Cuthill-McKee reordering. The matrix from the fluid domain has larger diagonal values than off-diagonals and is nearly symmetric. Eigenvalues and the condition number of the matrix from the fluid domain show that the $\mathbf{K}$ block has a significantly smaller condition number compared to the $\mathbf{A}_{f}$ matrix, supporting the main premise of BIPN. Effects of preconditioning on matrix characteristics are investigated by explicitly forming the preconditioned matrix. A diagonal preconditioner is shown to be effective to reduce the range of eigenvalues in the fluid domain, especially for the $\mathbf{K}$ matrix.

Adding wall deformability to the fluid simulation increases the bandwidth of the matrix and decreases the relative magnitudes of the diagonal values compared to the off-diagonal values. Due to the reduction of diagonal dominance, a diagonal preconditioner does not significantly reduce the condition number of the original matrix.

The resistance boundary condition disturbs the sparsity and diagonal dominance of the original fluid matrix, and causes an ill-conditioned system by adding an eigenvalue which is larger than the maximum eigenvalue of matrix without resistance BC. The resistance based preconditioner successfully reduces the condition number of the system with a resistance boundary condition by four orders of magnitude, while a diagonal preconditioner only reduces the condition number by two orders of magnitudes.

The performance of various preconditioned linear solvers is evaluated in four patient-specific models. In these models, RPC-BIPN is best for rigid wall models with multiple resistance or coupled LPN outlet boundary conditions. In the deformable wall simulation, RPC-BIPN shows significant performance degradation and diagonal preconditioners or ILU with BICG achieve the best performance.

This study motivates several new research directions to develop new preconditioned linear solver strategies. The effectiveness of BIPN for the solution of fluids problems with rigid walls has been proven in the current study. Currently our in-house code (RPC-BIPN) uses the bi-partitioned approach for ALE FSI, forming a linear system from the momentum equations for the fluid and the solid domains together, and another system for the continuity equation for the fluid domain. However, since the characteristics of the matrix in the solid domain is different from the fluid domain, most notably diagonal dominance, linear systems from these two domains should be separately solved (i.e. Tri-partitioning). The inefficiency of solving FSI in BIPN stems from adding off-diagonal dominance to the left-hand-side matrix block $\mathbf{K}$. Since RPC is based on a simple diagonal preconditioner, solving the system $\mathbf{K}$ becomes less efficient. We suggest solving $\mathbf{K}_{s}$ separately with an Incomplete Cholesky preconditioner, exploiting its symmetric property, rather than a simple diagonal preconditioning. Exploration of this idea is the subject of future work.

Additionally, we point out that the Schur complement block in BIPN is not preconditioned. Since a major portion of the computational cost in BIPN is consumed when solving the Schur complement block Esmaily2015 , acceleration of the linear solver performance by a proper preconditioning technique could significantly reduce the compute time. To form a preconditioner for the Schur complement block, one would need an efficient sparse matrix-matrix multiplication scheme as well as explicit formation of the Schur complement block. The open-source Trilinos library provides this option as well as various preconditioners so combining a partitioning approach with Trilinos is expected to provide consistent performance in both rigid and deformable wall simulations of cardiovascular hemodynamics. Implementation and testing of this approach is left for future investigation. Testing of linear solver performance on more complex patient-specific disease with large wall deformations (e.g. aortic dissection) are warranted, and would likely lead to further insights. Future studies are also warranted to further assess solver performance and matrix characteristics, towards development of new solver and preconditioner strategies.

9 Appendix

Appendix A GMRES restart

We tested different GMRES restart numbers in the pipe benchmark problems. In Figure 20, we plot compute times for preconditioned GMRES using a pipe model with rigid wall and a pipe with deformable wall. Our test shows that decreasing restart number increases the compute time of linear solver in the rigid pipe model. The FSI model does not show a notable difference.

Appendix B Choice of smoother and subsmoother for ML

In the ML package, multiple options are available for the smoother and the subsmoother. As shown in Figure 21, the Gauss-Seidel smoother works best among Chebyshev, symmetric Gauss-Seidel, and ILUT. For the subsmoother, the symmetric Gauss-Seidel is the best among Chebyshev and MLS.

Appendix C Effect of reordering in ILU

We evaluated and compared compute times of linear solvers with different reordering methods. RCM and METIS reordering for ILUT via the Trilinos IFPACK are implemented. We use 2 level fill-in and $10^{-2}$ dropping tolerance for this test. Figure 23 shows performance differences between ILUT with different reordering schemes. From the testing, we confirm that the RCM is the fastest method against METIS and no reordering. The superior performance of RCM is notable when GMRES is used with ILUT.

Appendix D Eigenvalue spectrums of the local and global matrix.

In this section we compare the spectrum of eigenvalues in the local and global matrices and investigate how our analysis on local eigenvalues can be generalized to the global matrix. We use a pipe model in the same dimension shown in Figure 1 meshed with 24,450 elements with $N_{nd}=5462$. We use one core to extract the global matrix, and four cores to examine local matrices. As shown in Figure 24, the eigenvalue distributions of the global and local matrices are similar. Although the eigenvalues from the global and local matrices are not exactly the same, the distribution of eigenvalues of local matrices is a good approximations to the distribution of eigenvalues in the global matrix.

Acknowledgements.

This work was supported by NIH grant (NIH R01-EB018302), NSF SSI grants 1663671 and 1339824, and NSF CDSE CBET 1508794. This work used the Extreme Science and Engineering Discovery Environment (XSEDE)XSEDE , which is supported by National Science Foundation grant number ACI-1548562. We thank Mahidhar Tatineni for assisting on building Trilinos on Comet cluster, which was made possible through the XSEDE Extended Collaborative Support Service (ECSS) program EECS . The authors also thank Michael Saunders, Michael Heroux, Mahdi Esmaily, Ju Liu, and Vijay Vedula, for fruitful discussions that helped in the preparation of this paper. The authors would like to thank the two anonymous reviewers whose comments greatly contributed to improve the completeness of the present study. We also acknowledge support from the open source SimVascular project at www.simvascular.org.