Simulations of cardiac function at the organ scale — electrophysiology, deformation, hemodynamics — are computationally demanding. Electrophysiology needs a fine spatio-temporal mesh to resolve the fast, steep depolarization wavefront (under 1 ms, under 1 mm); the non-linear equations of tissue deformation need less resolution, but are correspondingly more expensive to solve.

Clinical use compounds the cost: personalizing a model for diagnosis [Niederer et al, 2010], disease stratification [Arevalo et al, 2016], or therapy planning requires many forward simulations inside an optimization loop that re-parameterizes the model until it matches clinical data. We address this with efficient numerical methods that exploit high performance computing hardware, and with model formulations that need less resolution in the first place.

Finite Element Discretization

Patient-specific anatomical models are discretized with the finite element method (FEM) on unstructured, boundary-fitted grids (see image-based mesh generation). The same mesh serves both the electrophysiology and the nonlinear elasticity equations, so there is no need to map data between separate electrical and mechanical grids.

Figure 1: Finite element discretization

An unstructured triangular mesh, illustrating the general idea behind finite element discretization used throughout our models — the same principle applies to the tetrahedral meshes generated from patient imaging data. Public domain image by Oleg Alexandrov, via Wikimedia Commons.

Domain Decomposition

The sparse linear systems produced by FEM discretization are partitioned across compute nodes for parallel solution. We compare our own non-overlapping partitioning scheme, pt, against the widely used overlapping partitioning of PETSc — trading fragmented interface rows for a leaner memory footprint and less inter-node communication [Neic et al, 2012].

Figure 2: Domain-decomposition parallelization strategies

Parallelization strategies of the linear algebra package PETSc and the domain decomposition solver pt. Left panel: FE mesh of a domain consisting of two subdomains p₀ and p₁ and sparsity pattern of the corresponding global accumulated FE stiffness matrix A. In the FE mesh, black and white circles are inner nodes of p₀ and p₁, respectively, while squares are boundary nodes. Edges in the FE mesh correspond to nonzero entries of the stiffness matrix A. In the depicted sparsity pattern of A, white and black circles represent FE edges between inner nodes of p₀ and p₁, respectively, and squares represent FE edges from or to a boundary node. Right upper panel: Grid partitioning in the PETSc case. Submatrices A₀ and A₁, corresponding to p₀ and p₁, contain entire rows for all nodal indices assigned to a partition. Vector b₀ holds the entire local index range but for an operation equivalent to A₀·b₀ the node values corresponding to the triangles on p₀ need to be communicated, respectively, b₁ and p₁. Hatched areas in the FE mesh indicate domain overlap between p₀ and p₁. Right lower panel: Grid partitioning in the pt case. Domains p₀ and p₁ do not overlap. Therefore, rows corresponding to interface nodes in submatrices A₀ and A₁ are fragmented, but vectors b₀ and b₁ are complete, holding any entry corresponding to matrix entries in A₀ and A₁. From Neic et al, IEEE Trans. Biomed. Eng. 2012.

Scalable Solvers for Nonlinear Elasticity

Solving Cauchy’s equation of motion at high spatial resolution produces very large, ill-conditioned nonlinear systems. We developed a novel algebraic multigrid (AMG) preconditioner for an iterative Krylov solver, designed specifically for favorable strong scaling of both the setup and the solution phase [Augustin et al, 2016]. On whole heart electromechanics models of up to 184.6 million tetrahedral elements, it scales efficiently to 8192 compute cores — see strong scaling below.

Strong Scaling

Increasing the number of compute cores to reduce execution time is strong scaling. It is limited by the ratio of local compute work (volume) to communication cost (surface): as cores are added, the volume-to-surface ratio shrinks until communication starts to dominate — the limit of strong scalability — beyond which adding cores no longer helps, and can even hurt, execution time. At very high core counts, I/O cost grows for the same reason.

Optimizing strong scaling therefore means distributing workload and communication evenly, decomposing the domain into volumes with minimal surface area, and hiding the cost of I/O. The two figures below show this in practice, in two of our simulation codes.

Figure 3: Strong scaling from improved domain decomposition and I/O

Improved domain decomposition (DD) to reduce communication costs combined with asynchronuous I/O techniques to hide the otherwise significant cost of I/O due to blocking of computations led to significantly more favorable strong scaling properties. Shown are improvements over an implementation using a naive grid partitioning based on linear decomposition based on nodal indices and a standard sequential blocking I/O implementation. Efficient strong scaling could be improved from 1024 cores up to 16k cores. For details see Niederer et al, Front. Physiol 2011.

Figure 4: Strong scaling of the AMG-preconditioned EP/electromechanics solver

Strong scaling profile of an electromechanical model using a monodomain (left) or bidomain (right) EP model. Red, blue and green traces refer to CM, MM and FM simulation runs respectively, over a time frame of 500 ms. On the left hand side the assembly time for mechanics and the solve time (i.e. preconditioner setup and solving) for all linearized systems of equations is given. In addition, the total solving time of a monodomain electromechanical simulation is shown. Depicted in the right plot is the computational time for all elliptic systems, the computational time for all ODE systems and the computational time for all parabolic systems. In these three cases the computational time includes assembling, preconditioner setup and solving times. Solid lines correspond to the total solving time of a bidomain electromechanical simulation. From Augustin et al, J. Comput. Phys. 2016.

Use of Accelerators

GPUs are well suited to the large, sparse linear systems that arise from bidomain simulations. Porting our simulator to multi-GPU execution required only minor changes to the existing CPU code base: benchmarks on a state-of-the-art rabbit ventricle model showed bidomain simulations sped up by a factor of 11.8–16.3× on 6–20 GPUs versus the same number of CPU cores — matching the fastest 20-GPU run took 476 CPU cores on a national supercomputing facility [Neic et al, 2012].

Figure 5: GPU accelerator hardware

An opened GPU compute blade with glycol liquid cooling, illustrating the dense packaging of GPU and CPU modules typical of modern accelerated HPC systems. Photo by Steve Jurvetson, CC BY 4.0, via Wikimedia Commons.