Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2978074

Formats

Article sections

Authors

Related links

Int j numer method biomed eng. Author manuscript; available in PMC 2010 November 10.

Published in final edited form as:

PMCID: PMC2978074

NIHMSID: NIHMS109648

Center for Computational Fluid Dynamics, George Mason University, Fairfax, VA 22030-4444. USA, Email: ude.umg@tumf

The publisher's final edited version of this article is available at Int j numer method biomed eng

See other articles in PMC that cite the published article.

The study of hemodynamics in arterial models constructed from patient-specific medical images requires the solution of the incompressible flow equations in geometries characterized by complex branching tubular structures. The main challenge with this kind of geometries is that the convergence rate of the pressure Poisson solver is dominated by the graph depth of the computational grid. This paper presents a deflated preconditioned conjugate gradients (DPCG) algorithm for accelerating the pressure Poisson solver. A subspace deflation technique is used to approximate the lowest eigenvalues along tubular domains. This methodology was tested with an idealized cylindrical model and three patient-specific models of cerebral arteries and aneurysms constructed from medical images. For these cases, the number of iterations decreased by up to a factor of 16, while the total CPU time was reduced by up to 4 times when compared with the standard PCG solver.

Detailed patient-specific hemodynamic information is important for better understanding vascular diseases and to improve and personalize their treatment. Image-based computational fluid dynamics models have been used in a wide variety of applications[1]. The geometry of these models is typically characterized by a complex network of interconnected tubular structures. These geometries make the numerical solution of the flow equations computationally expensive. This is one of the reasons why these methods have not yet been introduced into the routine clinical practice. The objective of this work was to improve the performance of incompressible flow solvers for patient-specific hemodynamics calculations, making them more practical for clinical purposes.

Several numerical schemes have been used to solve the incompressible Navier-Stokes equations [2–4]. Many of them are based on so-called projection techniques, where the advancement of the flowfield in time is split into the following three substeps [5]:

*advective-diffusive*prediction:**v**→^{n}$$\left[\frac{\rho}{\mathrm{\Delta}t}-\nabla \mu \nabla \right]\phantom{\rule{0.16667em}{0ex}}(\stackrel{\sim}{\mathbf{v}}-{\mathbf{v}}^{n})+\rho {\mathbf{v}}^{n}\xb7\nabla {\mathbf{v}}^{n}+\nabla {p}^{n}=\nabla \mu \nabla {\mathbf{v}}^{n}$$(1)*pressure*correction:*p*→^{n}*p*^{n}^{+1}$$\begin{array}{r}\nabla \xb7{\mathbf{v}}^{n+1}=0\\ \rho \frac{{\mathbf{v}}^{n+1}-\stackrel{\sim}{\mathbf{v}}}{\mathrm{\Delta}t}+\nabla ({p}^{n+1}-{p}^{n})=0\end{array}$$which, by taking the divergence, results in$${\nabla}^{2}({p}^{n+1}-{p}^{n})=\rho \frac{\nabla \xb7\stackrel{\sim}{\mathbf{v}}}{\mathrm{\Delta}t}$$(2)*velocity*correction: →**v**^{n}^{+1}$${\mathbf{v}}^{n+1}=\stackrel{\sim}{\mathbf{v}}-\mathrm{\Delta}t\nabla ({p}^{n+1}-{p}^{n})$$(3)

Stabilization is achieved by a consistant edge-based discretization of the numerical fluxes in space which lead to a fourth order damping for the divergence constraint as described in [5]. The discretization of the so-called pressure Poisson equation (2) leads to a very large but symmetric system of linear equations. Therefore a system of linear equations of the form *Ax* = *b* needs to be solved at each timestep. The solution of this system is typically carried out using a preconditioned conjugate gradients (PCG) technique. The classical a priori bound for the error of the conjugate gradients method is given by [6]:

$$\parallel {x}_{\ast}-{x}_{k}\mid {\mid}_{A}\phantom{\rule{0.16667em}{0ex}}\le 2{\left[\frac{\sqrt{\kappa}-1}{\sqrt{\kappa +1}}\right]}^{k}\parallel {x}_{\ast}-{x}_{1}\mid {\mid}_{A}$$

(4)

where *x*_{*} is the exact solution, *x _{k}* is the approximate solution at the

The first paper to consider a deflation method for the CG solver is perhaps due to Nicolaides[8]. The main idea was to remove certain components of the initial residual that may impede convergence. It also introduced the idea of subdomain deflation, where the slow varying components of the solution are approximated by a coarser discretization of the domain. However, no numerical experiments were reported. The subdomain deflation approach was successfully applied to the bending of a cantilever beam and to the stationary Stokes problem, as reported in Mansfield[9]. Deflation was also used for an augmented conjugate gradients [10] where the Krylov subspace generated by a previous system is recycled for further solves in subsequent systems. An approximation of the eigenvectors is used to augment the space of subsequent systems to deflate the lowest eigenmodes [11]. In this case the eigenvectors are computed explicitly and not approximated through a coarser discretization of the domain. In this paper we present the application of the DPCG method to incompressible flows in tubular domains. A subspace deflation technique is used to approximate the lowest eigenmodes (along the vessels) in order to accelerate the convergence.

The Deflated Conjugate Gradients (DCG) technique starts by selecting a subspace * ^{n}* which is called the

Given an initial guess *x*_{0}, the initial error projection step is carried out by writing:

$${e}_{0}=W{c}_{0}+{e}_{1}$$

(5)

where *e*_{0} = *x*_{*} − *x*_{0} is the initial error, *W* is a matrix whose columns are a basis of and *c*_{0} is an *m*-dimensional vector (*m* = *dim*()). The vector *c*_{0} can be computed by enforcing *e*_{1} to be in ^{}* ^{A}*, i.e.

$${W}^{T}AW{c}_{0}={W}^{T}{r}_{0}$$

(6)

where *Ae*_{0} = *r*_{0} is the initial residual. The initial guess is then updated as follows:

$${x}_{1}={x}_{0}+W{c}_{0}$$

(7)

The remaining components of the error are eliminated by applying a CG algorithm whose search space is restricted to ^{}* ^{A}*. Recall that from the standard CG algorithm, each subsequent search direction

$${d}_{k+1}={r}_{k+1}+{\beta}_{k}{d}_{k};\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\beta}_{k}=\frac{{r}_{k+1}^{T}{r}_{k+1}}{{r}_{k}^{T}{r}_{k}}$$

(8)

where the parameter *β _{k}* is obtained by enforcing the conjugacy between subsequent search directions ((

$${d}_{k+1}={r}_{k+1}+{\beta}_{k}{d}_{k}-W{c}_{k+1}$$

(9)

where the parameter *β _{k}* is obtained as before and the

$${W}^{T}AW{c}_{k+1}={W}^{T}A{r}_{k+1}$$

(10)

Preconditioning can be easily included in the DCG algorithm by simply replacing the usual Euclidean inner product with the M-inner product, where M is an SPD preconditioner matrix [6].

The deflated preconditioned CG algorithm reads as follows:

1: | Given A, b, W, M and an initial guess x_{0} |

2: | Compute r_{0}:= b − Ax_{0} |

3: | Solve W^{T} AWc_{0} = W ^{T}r_{0} |

4: | Set x_{1}:= x_{0} + Wc_{0} |

5: | Compute r_{1}:= b − Ax_{1} |

6: | Solve Mz_{1} = r_{1} |

7: | Solve W^{T} AWc_{1} = W ^{T}Az_{1} |

8: | Set d_{1}:= z_{1} − Wc_{1} |

9: | do until convergence |

10: | α := (_{k}r)/(_{k}, z_{k}d)_{k}, Ad_{k} |

11: | x_{k}_{+1}:= x + _{k}α_{k}d_{k} |

12: | r_{k}_{+1}:=r − _{k}α_{k}Ad_{k} |

13: | Mz_{k}_{+1} = r_{k}_{+1} |

14: | β:= (_{k}r_{k}_{+1}, z_{k}_{+1})/(r, _{k}z)_{k} |

15: | W^{T}AWc_{k}_{+1} = W ^{T}Az_{k}_{+1} |

16: | d_{k}_{+1}:= z_{k}_{+1} + β − _{k}d_{k}Wc_{k}_{+1} |

17: | end do |

The DCG technique requires the definition of a deflation subspace . For matrices arising from a discretization of the domain (e.g. finite elements), the subdomain deflation approach described in [8] defines a deflation subspace based on a coarser discretization of the domain. This technique attempts to spatially approximate the smallest eigenmodes of the underlying system so that they can be removed from the search space of the conjugate gradients method.

Since the lowest eigenmodes of the Laplace operator in elongated or tubular domains are oriented along the domain longest direction, a good strategy for approximating these eigenmodes would be a coarser discretization of the domain along this longest direction. One of the simplest ways of defining a coarser discretization from a given unstructured mesh is to agglomerate the nodes of the mesh into subdomains. Two alternative methods have been developed following this strategy:

- Seedpoints Alternative (SA): For this (manual) technique, the user defines an arbitrary set of points, called
*seedpoints*. Given a mesh, the closest mesh points to the seedpoints are found, and a region number is assigned accordingly. Points not assigned to any region are then added one layer at a time until all points have been assigned a region number. - Advancing Layers Method (ALM): Starting from a prescribed point, neighboring points are added one layer at a time, until a specified number of points per region is exceeded. The last set of points added is then used as a starting point for the next group. The procedure is repeated until all points have been assigned a region number.

Figure 1 shows an schematics of the seedpoints alternative (top) and the advancing layers method (bottom) applied to a 2D pipe. The contours in this figure represent layers of elements being added to each deflation group. For the seedpoints alternative the deflation groups are grown in parallel starting from the seed points (black dots), while for the advancing layers method the deflation groups are build sequentially from left to right, where the start points (black dots) are selected automatically from the last layer of the last built group.

Schematics of the seedpoints alternative (top) and the advancing layers method (bottom) applied to a 2D pipe.

The last step is to define an interpolant polynomial over the deflation regions. The simplest choice for this is a constant function. Using this approach, each column of the matrix *W* will represent one deflation subdomain. For a given column, a unity value will be assigned if a point belongs to that region, and zero otherwise.

$${W}_{ij}=\{\begin{array}{cc}1& {p}_{i}\in {g}_{j}\\ 0& {p}_{i}\notin {g}_{j}\end{array}$$

(11)

The Deflated Preconditioned Conjugate Gradients (DPCG) was tested on an idealized pipe flow model and three patient-specific models of cerebral arteries and aneurysms constructed from medical images using a previously developed methodology[12]. For all the cases a simple diagonal preconditioned was used. It was also verified that the results obtained using the deflated PCG solver coincide with those obtained by the standard PCG solver.

For the hemodynamic simulations, the blood was considered as a Newtonian fluid, which is modeled by the *unsteady* incompressible Navier-Stokes equations. Pulsatile flow boundary conditions were computed using a Womersley model with measured flow rates of a normal subject as described in [12]. The simulations were performed using implicit timestepping, solving a pseudo-steady problem at each timestep [13]. Within each pseudo-timestep the advective terms were integrated implicitly using 5 LU-SGS passes (local Courant number *C* = 5.0), followed by the pressure projection. The timestep was set to Δ*t* = 0.01 sec. The material properties of the fluid were taken to be *ρ* = 1.0 g/cm3 and *μ*= 0.04 Poise. Two cardiac cycles were computed with 100 timesteps per cycle.

The aim of these examples was the comparison of the number of iterations of the pressure Poisson solver and the speed of the flow calculations. An Intel Xeon processor (E5345/2.33MHz) computer with 16GB of RAM was used in all the examples.

The first example is the classic Poiseuille pipe flow, a steady flow of a viscous Newtonian fluid in a straight circular domain. A uniform velocity profile is prescribed at the inflow, while a constant pressure is prescribed at the outflow. Since the pressure field has to establish itself along the pipe, the number of iterations required increases with the graph depth of the finite element mesh. The physical dimensions and parameters were set as follows:

- pipe radius:
*r*= 1 - pipe length:
*l*= 20, 40, 80 - density:
*ρ*= 1 - inflow velocity:
*v*= 1 - viscosity:
*μ*= 0.01

The element size was set to *h* = 0.1, implying approximately a graph depth of 200, 400 and 800 for the cases considered. This resulted in grids of 129 Kels, 260 Kels and 516 Kels respectively. All cases were run for 100 timesteps using explicit timestepping. The number of groups were chosen to be 15, 30 and 60. The deflation subdomains were generated by the advancing layers technique, starting from the exit. Figure 2–3 show the surface mesh, deflation subdomain boundaries for 15 groups, pressure and absolute value of the velocity for *l* = 20. As expected the velocity profile at the outlet is parabolic, and the pressure field is linear along the cylinder axis after an initial development distance.

Figure 4 shows the number of iterations required for the PCG solver. The sudden ‘dips’ in the number of iterations at some timesteps are due to the fact that a projective prediction of pressure increments with 2 Krylov vectors [14] was used. Note the dramatic decrease in the number of iterations achieved by the deflated PCG solver. This decrease may also be seen in Figure 5, which depicts the average number of iterations for the first 20 steps for the different options chosen.

Pipe Flow: Average number of iterations required for the PCG solver as a function of the pipe length.

While the number of iterations increases linearly with the pipe length for the conventional PCG, the performance of the deflated PCG solver seems to be insensitive to the pipe length and the number of groups chosen. Figure 6 shows the total CPU time required for the simulation, highlighting the importance of a fast pressure Poisson solver. Note that for the case *l* = 80, the deflated PCG case performs seven times faster.

The second example is a patient-specific model of a cerebral aneurysm located at the Internal Carotid Artery (ICA). This case is particularly interesting since it was possible to reconstruct the entire ICA from the carotid bifurcation to the carotid terminus. The vascular model was reconstructed from 3DRA images.

The volume mesh obtained for this model had 646 Kpts and 3.6 Mels. The deflation groups were generated using both methods. For the seedpoints alternative 45 points were manually selected from the surface model. For the advancing layers method 20, 50, 100 and 150 groups were automatically generated. Figure 7 shows the overall domain, as well as the group boundaries for the 45 (top left) and 50 (top right) groups cases.

ICA Aneurysm: Top: deflation boundary subdomains for 50 groups using advancing layers method (left) and 45 groups using the seedpoints alternative (right). Bottom: Pressure (left) and Wall Shear Stress distribution (right) at peek systole.

The pressure was prescribed (homogeneous bc) at the outflow boundary at the ICA terminal (top part), while a time-dependent velocity profile was prescribed at the inflow boundary at the origin of the ICA (bottom part). No-slip boundary conditions were applied at the vessel wall.

Figure 7 depicts the pressure drop (bottom left) and the wall shear stress distribution (bottom right) at the peak systole (the inflow rate peak) of the second cardiac cycle. As it was expected, the pressure gradient is aligned with the parent vessel, while the pressure is almost constant inside the aneurysm.

Figure 8 shows the average number of iterations per time step (top) for the pressure Poisson solver and the total CPU time (bottom) needed for completing the 200 timesteps. In this case, the number of iterations were reduced by up to 8.5 times while the CPU time was reduced by up to 3.2 times. This result corresponds to the case with 150 groups using the advancing layers method.

The third example is a patient-specific model of a cerebral aneurysm located at the tip of the basilar artery constructed, as in the previous example, from 3DRA images. The model included in the upstream direction both left and right vertebral arteries. The volume mesh generated for this model had 535 Kpts and 2.9 Mels. The deflation subdomains were generated using both methods. For the seedpoints alternative 90 points were manually selected from the surface model. For the advancing layers method 20, 50, 100 and 150 groups were automatically generated. Figure 9 shows the overall domain, as well as the deflation group boundaries for the 90 and 100 groups cases.

BT Aneurysm: Top: deflation boundary subdomains for 100 groups using the advancing layers method (left) and 90 groups using the seedpoints alternative (right). Bottom: Pressure (left) and Wall Shear Stress distribution (right) at peek systole.

The pressure was prescribed (homogeneous bc) at the three outflows (top part), while a time-dependent velocity profile was prescribed at the two in-flows (bottom part). No-slip boundary conditions were applied at the vessel wall.

Figure 9 depicts the pressure (left) and the wall shear stress distribution (right) at the peak systole of the second cardiac cycle. As it was expected, the pressure changes along the vessels, while it is almost constant inside the aneurysm.

Figure 10 shows the average number of iterations per time step (top) for the pressure Poisson solver and the total CPU time (bottom) required to compute 200 timesteps. The highest reduction in the number of iterations was about 16 times while the total CPU time was reduced by up to a factor of 4.2 times. Notice that in this case, the best result was obtained with 90 groups using the seedpoints alternative.

The last example is a subject-specific model of the circle of Willis of a normal volunteer constructed from magnetic resonance images. The volume mesh generated for this model had approximately 836 Kpts and 4.6 Mels. The deflation groups were generated using both methods. For the seedpoints alternative 131 points were manually selected from the surface of the model, while for the advancing layers method 20, 50, 100 and 150 groups were automatically generated. Figure 11 shows the overall domain, as well as the deflation subdomain boundaries for the 131 and 150 groups cases.

Circle of Willis: Top: deflation boundary subdomains for 100 groups using the advancing layers method (left) and 131 groups using the seedpoints alternative (right). Bottom: Pressure (left) and Wall Shear Stress distribution (right) at peek systole.

The pressure was prescribed (homogeneous bc) at the six outflows, while a time-dependent velocity profile was prescribed at the three inflows (bottom part). No-slip boundary conditions were applied at the vessel wall.

Figure 11 depicts the pressure (left) and the wall shear stress distribution (right) at the peak systole of the second cardiac cycle. As it was expected, the pressure changes slowly on the large vessels when compared with the small vessels.

Figure 12 shows the average number of iterations per time step (top) for the pressure Poisson solver and the total CPU time (bottom). The highest reduction in the number of iterations was about 10 times while the total CPU time was reduced by up to a factor of 3.9 times. Notice that as in the previous case, the best result was obtained using the seedpoints alternative.

A Deflated Preconditioned Conjugate Gradients (DPCG) technique has been developed for the pressure Poisson equation within an incompressible flow solver. A simple idealized pipe flow model and three patient-specific image-based blood flow computations were carried out to assess the performance of the deflated PCG solver. The number of iterations was significantly reduced (up to a factor of 16) while speedups of about 4 times were obtained in the solution of the incompressible flow equations.

Automatic (advancing layers) and manual (seedpoints) methods for constructing the deflation subdomains in elongated or tubular domains have been developed. These methods resulted in similar performances. However, in some cases, the manual method outperformed the automatic one, possibly because it produced more regular groups with more uniform pressure distributions.

In conclusion, this methodology has improved substantially the performance of the incompressible flow solver especially for vascular domains, making larger and more complex hemodynamic models practical. Although the methodology was developed in the context of blood flow simulations, it is general and can be used for any incompressible flow calculation.

We thank Philips Medical Systems for financial support.

1. Steinman DA, Taylor CA. Flow imaging and computing: Large artery hemodynamics. Ann Biomed Eng. 2005;33(12):1704–1709. [PubMed]

2. Thomaset F. Implementation of Finite Element Methods for Navier-Stokes Equations. Springer-Verlag; 1981.

3. Gunsburger MD, Nicolaides R. Incompressible Computational Fluid Dynamics: Trends and Advances. Cambridge University Press; 1993.

4. Hafez MM. Numerical Simulation of Incompressible Flows. World Scientific; 2003.

5. Lohner R. Applied CFD techniques: An Introduction based on Finite Element Methods. 2. Wiley; 2008.

6. Saad Y. Iterative methods for sparse linear systems. 2. SIAM; 2003.

7. van der Sluis A, van der Vorst HA. The rate of convergence of conjugate gradients. Numerische Mathematik. 1986;48:543–560.

8. Nicolaides RA. Deflation of conjugate gradients with applications to boundary value problems. SIAM J Numer Anal. 1987;24(2):355–365.

9. Mansfield L. On the use of deflation to improve the convergence of the conjugate gradient iteration. Comm Appl Num Meth. 1988;4:151–156.

10. Erhel J, Burrage K, Pohl B. Restarted gmres preconditioned by deflation. J Comput Appl Math. 1996;69:303–318.

11. Saad Y, Yeung M, Erhel J, Guyomarc’h F. A deflated version of the conjugate gradient algorithm. SIAM J Sci Comput. 2000;21(5):1909– 1926.

12. Cebral JR, Castro MA, Appanaboyina S, Putman CM, Millán D, Frangi AF. Efficient pipeline for image-based patient-specific analysis of cerebral aneurysm hemodynamics: Technique and sensitivity. IEEE Transactions in Medical Imaging. 2005;24(1):457–467. [PubMed]

13. Löhner Rainald, Yang Chi, Juan Raul Cebral, Camelli Fernando, Soto Orlando, Waltz J. Improving the speed and accuracy of projection-type incompressible flow solvers. Comp Meth Appl Mech Eng. 2006;195(23–24):3087–3109.

14. Löhner Rainald. Projective prediction of pressure increments. Comm Num Meth Eng. 2005;21(4):201–207.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |