A blog by Univ.-Prof. Dipl.-Ing. Dr. Gundolf Haase, University of Graz
One of the objectives of the Mont-Blanc project is to assess the behaviour of real applications on the different generations of Arm platforms made available by the project. As part of Mont-Blanc 3, the University of Graz worked on biophysical applications, and specifically on the particularly complex cardiovascular simulations.
The modelling and simulation of cardiovascular phenomena is a coupled multi-physics problem that consists of electrophysiology (excitation pattern in the heart), of non-linear elasticity (deformation of heart tissue) and of fluid dynamics (blood flow through heart chambers).
The electrophysiology is described via the bidomain equations that consist of two PDEs coupled with a system of ODEs (non-linear coupled PDE-ODE system):
Here φe and Vm are the extracellular potential and transmembrane voltage respectively; ~η represents the ionic current variables; σ¯i and σ¯e are conductivity tensors of intracellular and extracellular spaces respectively; σb is the isotropic conductivity of the fluid in which the heart is immersed (bath and cavities); Cm is the capacitance per unit area and β is surface to volume ratio; iion and g model ionic currents and specify the cell membrane model.
There are cheaper models as the monodomain equations or the Eikonal equation available in case only the arrival time of the excitation wave (=stimulation time of the heart muscles) is needed.
Cardiac electromechanics couples electrophysiology with the non-linear deformation, i.e., the bidomain equations are extended with the elasticity equation for large deformations and for non-linear material laws. The deformation is governed by the equilibrium equation
−divσ(u) = b,
with unknown displacement field u. The stress tensor σ = σp + σa consists of a passive contribution σp and an active contribution σa, while the latter depends on the normalized myocyte orientation. The passive contribution requires the right Cauchy-Green tensor taking into account the deformation gradient and a nonlinear function Ψ describing the strain energy density of the specific material.
The full cardiovascular model requires a Fluid-Structure interaction (FSI) and adds the Navier-Stokes equations to the coupled system from electromechanics including coupling terms at the interfaces between fluid domain and deformation domain. The fluid (blood) is assumed as incompressible.
The conditions on the interface between fluid domain and deformation domain have to guarantee continuity of pressure and deformation via the related physical quantities. The FSI has to be performed on changing computational domains which requires some caution in the numerical realization, especially a high quality mesh deformation for the moving discretized computational domains.
Discretization and numerical realization
The overall computational framework is embedded in the Cardiac Arrhythmia Research Package (CARP) (G. Plank, Graz and E. Vigmond, Bordeaux)  where the University of Graz contributes with parallel solvers (G. Haase, Graz; MontBlanc3) [2,3] and image processing [C. Bredies, Graz].
The computational domain is determined a priori by image processing (registration + segmentation) and physiological mapping. This continuous domain is discretized by finite elements (tetrahedrons) that define the computational domains. We have to distinguish between the fluid domain, the deformation domain and the electrical domain.
From images to the heart geometry
An electrical excitation wave on a discretized heart
The finite element discretization of the PDEs uses, linear elements, stable elements for saddle point problems, and higher order elements regarding the regularity requirements on the solution (that is, analysis of PDEs).
Independently from the model under investigation, we always have to solve a huge linear system of linear equations in each time step and in each step of the outer non-linear iteration. Solving the (block) systems takes advantage of the sparse matrix structure and the properties of the originating PDEs (symmetric?, positive definite?). The solution methods for the linear systems of equations are preconditioned Krylov subspace as the conjugate gradient (cg) or the generalized minimal residual (gmres) method with algebraic multigrid (AMG) as preconditioner.
Thanks to the Mont-Blanc 3 project we have been able to improve our solver package and parts of it (AMG, Jacobi smoother, Eikonal solver) with performance tools from other workgroups in MontBlanc3 on available and brand-new hardware as the Arm Clusters at BSC.
High Performance Computing
Solving the above-mentioned discretized (coupled) problems requires a lot of compute power, e.g., it takes 8 minutes on 16192 CPU cores (SuperMuc in Munich, Germany) to simulate a full electromechanic heart beat on a realistic discretization with 45 Million tetrahedrons. The implementation is usually done in C++ and uses shared memory parallelization (OpenMP, OpenSs) as well as distributed memory parallelization (MPI) and occasionally GPU-parallelization (CUDA, OpenACC) [2,3].
A 4×4 domain decomposition
Some strong scaling results for the electrophysiological subproblem on SuperMuc
Further topics (ongoing research)
In contrast to engineering applications, the material parameters (σ etc.) are only roughly known because dead tissue can be measured but behaves differently than living tissue, and experiments on beating human hearts are forbidden by themselves.
Therefore, simulation results have to be compared to measurable data (here, the ECG at 12 spots on the torso) and the material data will be calibrated such the computed ECG gets closer to the measured ECG. We follow currently several approaches for an automatic calibration (with reduced sets of equations) which contain
- optimization approaches,
- deep learning approaches,
- combinations of the two.
Here the highly optimized Eikonal solver  (based on several simplification of the bidomain equations) is the first choice to calculate a fast and sufficiently accurate electrical excitation wave, see figure. The Arm cluster outperforms the performance of other hardware platforms for our Eikonal solver in terms of power consumption.
 E. Vigmond, M. Hughes, G. Plank and L. Leonc: “Computational tools for modeling electrical activity in cardiac tissue”, Journal of Electrocardiology. 36. 2003. pp. 69-74
 G. Plank, M. Liebmann, R. Weber dos Santos, E.J. Vigmond and G. Haase: “Algebraic Multigrid Preconditioner for the Cardiac Bidomain Model“, IEEE Transactions on Biomedical Engineering, Vol 54, Issue 4, pp. 585-596 (2007).
 C.M. Augustin, A. Neic. M. Liebmann, A.J. Prassl, S.A. Niederer, G. Haase and G. Plank: “Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation“, Journal of Computational Physics. 305. 2016. pp. 622-646. doi: 10.1016/j.jcp.2015.10.045
 D. Ganellari, G. Haase and G. Zumbusch: “A massively parallel Eikonal solver on unstructured meshes”, Computing and Visualization in Science, pp. 1-16, 2018, doi: 10.1007/s00791-018-0288-z