Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020

jh200036-MDHI

Highiresolutionisimulationioficardiacielectrophysiologyionirealisticiwhole-hearti
geometriesi
Kengo Nakajima（The University of Tokyo）
Abstract
This international project in JHPCN aims to combine the expertise of The University ofTokyo in HPC and the expertise of Simula Research Laboratory (Norway) in cardiac
modeling, with the objective of enhancing a 3D simulator of cardiac electrophysiology overrealistic whole-heart geometries. The enhanced 3D simulator is expected to efficiently use
modern supercomputers, such as Oakforest-PACS and Oakbridge-CX, to simulate realistic
scenarios of electrophysiology in the heart. These high-resolution and biologically-detailed
simulations are considered as an important tool for advancing the scientific understandingof the electrophysiology in the heart, and eventually for improving the medical treatment
and drug design of various heart diseases.
The start of this 2-year project is an existing simulator of cardiac electrophysiology. We
aim to carry out a series of improvements with respect to both the software implementationand the underlying numerical strategy. Specifically, SIMD vectorization will be enabled todeliver high single-core performance of the non-memory-traffic constrained computationaltasks. OpenMP parallelization will be combined with MPI-based parallelization to achieve
optimal single-node performance. Multi-thread-tasking will be investigated, together withsuitable data restructuring, for alleviating the MPI communication overhead in multi-node
simulations. Implicit integration in the time direction will be adopted and implemented toallow larger timesteps, thereby giving the potential of overall time saving.

1.

Basic Information

Masaharu Matsumoto (U Tokyo): Numericalalgorithms and parallel programming.

Glenn Terje Lines (Simula/Norway): Cardiac
electrophysiology, mathematical modeling.



Johannes Langguth (Simula/Norway): Codeparallelization, profiling and optimization.

Jonas van den Brink (Simula/Norway): Preparationof geometries and physiological parameters for
subcellular simulations.

(1) Collaborating JHPCN Centers
Information Technology Center, Univ. Tokyo(2) Research Areas

Very large-scale numerical computation

Very large-scale data processingVery large capacity network technology

Very large-scale information systems

(3) Roles of Project Members
 Kengo Nakajima (U Tokyo): Project administration,
numerical algorithms and parallel programming.

Kristian Gregorius Hustad (Simula/Norway): Codeparallelization, profiling and optimization.

Hermenegild Arevalo (Simula/Norway): Preparationof geometries and parameters, runningsimulations, result analysis.

2.

Purpose and Significance of ResearchXing Cai (Simula/Norway): Numerical algorithms,
code parallelization and optimization, as well as
project coordination together with Prof. Nakajima.

Akihiro Ida (U Tokyo): Numerical algorithms andparallel programming.

Coordinated electrical activities in the heart are vital for

Toshihiro Hanawa (U Tokyo): Code parallelization,
profiling and optimization.

electrophysiology require extremely fine spatial and

its function. Biophysically accurate simulations of cardiac

Masatoshi Kawai (U Tokyo): Numerical algorithms
and parallel programming.

temporal resolutions. In addition, unstructured

Tetsuya Hoshino (U Tokyo): Code parallelization,
profiling and optimization.

capture the realistic 3D geometry of the heart. The

computational meshes must be adopted to precisely

relevant mathematical model in this context is a nonlinear

1

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020

3D reaction-diffusion equation, with the transmembrane

Norwegian partner. Second, the Oakforest-PACS andelectrical potential being the primary unknown field,

Oakbridge-CX systems (plus the upcoming Wisteria

while calcium handling in each computational cell issystem) are of a suitable size for achieving the ambitious

described by a system of nonlinear ordinary differential

goal of this project, whereas access to world-leading

equations (ODEs). Such a whole-heart simulation, if

supercomputers has been very scarce for the Norwegianusing tens of millions of computational cells, can easily

partner. The two hardware systems also provide atake hours or even days of execution time by current

valuable testbed for the performance portability of themainstream simulation codes when running on a small or

simulation codes to be developed. The high-speed filemedium-scale cluster of multicore CPUs. The slow turn-cache systems available at UTokyo also provide goodaround times have so far limited the use of suchpossibilities of in-situ huge-scale data analysis.

ambitious simulations in a realistic setting where
computational cardiologists can for example experiment,

4.

“in-silico”, about how/when/where arrhythmia arises, or

jh20036 is is partially related to, but not a direct

the sensitivity of certain dysfunctions of the heart with

continuation of, JHPCN project jh180024/jh190040

respect to the calcium channel mechanism and/or the(Physiologically realistic study of subcellular calcium

conductivity properties of the heart.

dynamics with nanometer resolution), which has obtained

In this project, we aim to increase the scale of state-

Outline of Research Achievements up to FY2019

good results about optimizing cardiac simulations onof-the-art whole-heart simulations by tenfold, i.e., using

regular computational meshes.

hundreds of millions of computational cells instead oftens of millions. At the same time, we have the ambition

5.

of reducing the turn-around time from hours to minutes.

(1) Overview

We will achieve these ambitious goals by combining

Overall, the 3D simulator of cardiac electrophysiology

performance-enhancing numerical strategies, hardware-carries out a time integration where each iteration

aware code optimization and parallelization-overhead

consists of first individually solving a system of nonlinear

reduction. It is also remarked that the core computations,

ODEs (the so-called cell model that models thei.e., numerically solving partial differential equationstransmembrane ionic currents and cellular calcium(PDEs) on unstructured meshes and millions of ODE

handling) on each computational cell, and then solving asystems, are present in many other computational science

3D diffusion equation (PDE) that couples all theapplications. This means that the advances achieved incomputational cells in a realistic whole-heart geometry.

this proposed project will extend well beyond the domain

The research achievements during FY2020 center aroundof cardiac electrophysiology simulations.

improving the performance of the ODE and PDE

Details of FY2020 Research Achievements

computations, which are detailed below. Generally,3.

problems are solved explicitly in time direction, where

Significance as JHPCN Joint Research Project

The necessity of implementing this JHPCN joint researchinvserins of large-scale matrices are not required.

project is due to two aspects. First, UTokyo has world-

Because implicit integration in the time direction allows

leading expertise in implementing and optimizing

larger timesteps, thereby giving the potential of overall

advanced numerical code. This expertise has been built

time saving, we plan to introduce and implement implicit

up via developing real-world applications for running on

integration in time direction. In FY.2020, preliminary

cutting-edge supercomputers at UTokyo. Such hands-on

studies for implicit methods using multigrid methods

experience on supercomputing is lacking for thewere also conducted.

2

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020

(2) Enabling SIMD vectorization of ODE computation

512 intrinsics are automatically enabled). Therefore, theThe numerical solution of the ODEs per computational

automatic vectorization approach is favored, due to its

cell typically adopts explicit time integrators such as thesimplicity and minimum code restructuring required.

forward Euler (FE) method or the general first-order

In the following, we report the achieved ODE

Rush-Larsen (GRL1) method. The associatedcomputation performance in various tests. These

computational intensity (number of floating-point

experiments cover both Oakforest-PACS and Oakbridge-

operations per memory load/store) is relatively high. ThisCX, and three biologically detailed cell ODE models aremeans that enabling SIMD vectorization on modern CPU

investigated: the 19-state Ten-Tusscher-Panfilov modelarchitectures is important for achieving goodfrom 2006 (TT06), the 24-state Jæger-Tveito model (JT),performance of the ODE computation, which was lacking

and the 39-state Grandi-Pasqualini-Bers model (GPB).

on the existing simulator before the project start.

The used ODE solvers are the forward-Euler (FE) method

There are in general two strategies for SIMD

and the generalized first-order Rush-Larsen (GRL1)

vectorization. The first is to rely on the automatic

method. For all the experiments, instead of reporting thevectorization capability of compilers, whereas the second

detailed time results, we use the performance metric as

is to directly insert SIMD vectorization intrinsics (such as

millions of cell steps executed per second (the higher the

AVX2 or AVX-512) or use a portable vectorization

better). We can see from Figures 1-8 that compiler-

library such as VCL. During FY2020 we have tested both

enabled automatic vectorization provides substantial

strategies of SIMD vectorization.

performance improvement over a naïve version without

To help trigger automatic vectorization by the compiler,

SIMD vectorization. The Oakbridge-CX system provides

we have found that it is important to restructure thebetter ODE performance than the Oakforest-PACS

overall data structure in the following way. Instead ofsystem, in both single-thread and single-node scenarios,

storing together the various ODE state variables for eachindependent of the cell model and the ODE solver.

cell (such that the data structure for all the cells is “anmethod.

array of structs”), which is adopted by the old simulator,
we have grouped each state variable from all the cells
together (such that the overall data structure is now “astruct of arrays”). Another useful hint given to thecompiler is the omp simd directive associated with
OpenMP parallelization of the loop that traverses all thecells. (The ODE computation is embarrassingly parallelover the cells.) On the Oakbrige-CX system, in order toencourage the compiler to automatically adopt AVX-512
intrinsics instead of AVX2, we found it necessary to add
the clause of simdlen(8) in addition, which already seems

Figure 1 Single-thread ODE performance on Oakbridge

to be the default behavior of the Intel compiler on
Oakforest-PACS. For comparison, we have also used the
VCL library to explicitly vectorize the ODE computation.
Experiments have shown that this manual approach does
not bring any noticeable performance benefits over the
automatic compiler vectorization (at least when AVX-

3

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020

Figure 2 Single-thread ODE performance on Oakbridge

Figure 5 Single-node ODE performance on Oakbridge

Figure 3 Single-thread ODE performance on Oakforest

Figure 6 Single-node ODE performance on Oakbridge

Figure 7 Single-node ODE peformance on Oakforest
Figure 4 Single-thread ODE performance on Oakforest

4

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020

Figure 8 Single-node ODE performance on Oakforest

Figure 9 Single-node PDE performance on Oakbridge(3) Mixing OpenMP and MPI parallelization for PDE

On the Oakforest-PACS system, good single-node PDE

computation

performance is also associated with using a single MPI

Each compute node of the Oakbridge-CX and Oakforst-

process that spawns at least 68 OpenMP threads. In

PACS systems has, respectively, 56 and 68 cores. The

Figure 10, we further study the effect of using more thanlatter is also well known to support 4 threads per core.

one thread per core. It shows that the best single-node

The old simulator focused on using MPI parallelization

PDE performance is associated with using one MPI

both inter-node and intra-node, thus having the risk ofprocess that spawns in total 272 OpenMP threads.

extensive MPI overhead when using a large number ofcompute nodes. We have therefore enabled
OpenMP+MPI parallelization for the PDE computation(the same as for the ODE computation).
In Figure 9, we study the obtained PDE performance on asingle compute node of Oakbridge-CX, where we change
the number of MPI processes used while keeping the total
number of OpenMP threads at 56. The reported PDE
performance uses the metric of “effective memory
bandwidth achieved”, which is calculated as the
Figure 10 Single-node PDE performance with respect toOpenMP threads per coreminimum incurred amount of memory traffic (assuming
perfect caching) divided by the time used. We can seethat using one MPI process that spawns 56 OpenMP(4) Multi-node strong scaling study

threads gives the best single-node PDE performance. It

We have also measured the strong scaling results of

can also be observed that the effective memory

running a parallel simulation (ODE+PDE) of cardiac

bandwidth achieved is quite close to the STREAM

electrophysiology, which uses a realistic whole-heartbenchmark measured memory bandwidth, meaning thatgeometry that involves over 14 million computational

the single-node PDE performance is approaching its

cells (each being a tetrahedron). As shown in Figure 11,

realistic upper limit on Oakbridge-CX.

up to 256 compute nodes have been used on OakbridgeCX, where each node always uses a single MPI processthat spawns 56 OpenMP threads. The performance metric5

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020

used is number of cell steps computed per second. It can

(a)

(b)

(c)

(d)

observed a gap between the actual parallel performanceCsand the perfect scaling, mostly indicating the impact of
MPI communication overhead.

Figure 12 Formats of sparse matrix storage. (a) Compressedrow storage (CRS); (b) Ellpack-Itpack (ELL), (c) Sliced ELL,(d) SELL-C- (SELL-2-8)
6.00

(a)

Lev.1

Lev.2

Lev.3

Lev.4

Rest

5.00

Figure 11 Strong scaling measurements on Oakbridge

4.00
sec.

(5) Multigrid Method for Impicit Time Integration
The parallel multigrid method is expected to play animportant role in scientific computing on exa-scale

3.00

2.00

supercomputer systems for solving large-scale linear

1.00

equations with sparse coefficient matrices. Because
0.00

solving sparse linear systems is a very memory-bound

CRS

ELL

6.00

process, efficient method for storage of coefficient(b)matrices is a crucial issue. Nakajima et al. implemented

Lev.1

SCS-a:
C=8

Lev.2

SCS-b:
C=8

Lev.3

SCS-b:
C=128

Lev.4

Rest

5.00

sliced ELL method to parallel conjugate gradient solvers

4.00

sec.

with multigrid preconditioning (MGCG) for theapplication on 3D groundwater flow through3.00

heterogeneous porous media (pGW3D-FVM), and

2.00

excellent performance has been obtained on large-scale

1.00

multicore/manycore clusters [Nakajima, K., IEEE
0.00
CRS

ICPADS 2014, 2014]. In the present work, we introduced

ELL

SCS-a:
C=8

SCS-b:
C=8

SCS-b:
C=128

Figure 13 Elapsed Computation Time for MGCG Solver: 8 nodesof OFP (HB 416), 33,554,432 DOF, (a) Double Precision(FP64), (b) Single Precision (FP32)

SELL-C- with double/single precision computing to the
MGCG solver, and evaluated the performance of thesolver with OpenMP/MPI hybrid parallel programing

Figuremodels on the Oakforest-PACS (OFP) system using up to13(a) and Figure 13(b) show computation

2,048 nodes of Intel Xeon Phi. Because SELL-C- istime of MGCG for HB 416 with 8 nodes of OFP, where

suitable for wide-SIMD architecture, such as Xeon Phi,

total problem size is 33,554,432 DOF. Computation time

improvement of the performance over the sliced ELL wasof MGCG for HB 416 with 8 nodes of OFP. Each of Lev.h

more than 35% for double precision and more than 45%(h=1 (finest level)~4 (coarser level)) shows total time forfor single precision on OFP.

smoothing at the h-th level of multigrid computing, while
Rest includes time for communications, coarse grid solver
and conjugate gradient solver except multigrid. SELL-8-8

6

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020and SELL-128-128 are applied to SCS-b. Improvement of(C= =128) with HB 416. In both of double precision andperformance of MGCG solver over CRS is shown in Table

single precision, performance of SCS-b (C==8) and SCS-

1. Generally, improvement of performance by SCS-b isb (C= =128) is competitive, but SCS-b (C= =128) isslightly faster for single precision.

excellent, and the effects are more significant in FP32

10.00

cases. Generally, performance improvement by single
8.00

precision for CRS is lower compared to Sliced ELL and

6.00

sec.

SELL-C-.

CRS-d
ELL-d
SCS-a(C=8)-d
SCS-b(C=8)-d
SCS-b(C=128)-d
SCS-b(C=128)-s

Table 1. Summary of Results in Fig.9 for OFP

4.00

a. Improvement over CRS (MGCG solver, Level-1 of Smoother)

Sliced ELL

Double Precision(FP64)

Single Precision(FP32)

36.6%, 46.6%

15.6%, 44.7%

2.00

0.00
1.E+01

SCS-a (C==8)

79.4%, 101.3%

58.7%, 81.0%

SCS-b (C==8)

84.9%, 108.7%

90.9 %, 152.2%

SCS-b (C==128)

90.1%, 107.8%

137.5%, 174.9%

1.E+02

1.E+03
1.E+04
Core #

1.E+05

1.E+06

Figure 15 Elapsed Computation Time for MGCG Solver, Weak
Scaling up to 2,048 nodes, Max. Problem Size: 8,589,934,592
DOF, HB 416, “-d”: Double Precision (FP64), “-s”: Single
Precision (FP32)

b. Improvement over Sliced ELL (MGCG solver, Level-1 of
Smoother)

Figure 16 compares HB 416 and HB 88 for SCS-b

SCS-a (C==8)

Double Precision(FP64)

Single Precision(FP32)

31.4%, 39.2%

20.3%, 25.1%

SCS-b (C==8)

35.4%, 44.4%

49.6%, 74.3%

SCS-b (C==128)

39.2%, 46.4%

59.4%, 89.9%(C= =8) with double precision and for SCS-b (C==128)
with single precision. Generally speaking, HB 416 isfaster than HB 88, but performance is similar at 1,024 and
2,048 nodes. Problem of CGA (Coarse Grid Aggregation)

Performance of weak scaling has been evaluated using

adopted in this work (Figure 17) is that the coarse grid

up to 2,048 nodes of OFP (131,072 cores). Maximum

solver works on a single MPI process and, problem size ofproblem size is 8,589,934,592 DOF. Figure 14 shows

the coarse grid solver is proportional to total number ofresults of weak scaling up to 2,048 nodes of OFP using HB

MPI processes. Therefore, cost of the coarse grid solver is

416, where “-d” denotes double precision (FP64), and “-

more significant with many nodes. Moreover, this effect is

s” is for single precision (FP32).

more significant, if number of threads for each MPI

5.50
5.00
4.50

SCS-b(C=8)-d
SCS-b(C=128)-d
SCS-b(C=8)-s
SCS-b(C=128)-s

process is smaller, such as HB 416 in this work.
5.00

sec.

4.00

4.00

3.50

sec.

3.00
2.50
2.00
1.50
1.E+01

3.00
HB 4x16-d
HB 4x16-s
HB 8x8-d
HB 8x8-s
2.00
1.E+02

1.E+03
1.E+04
Core #

1.E+05

1.E+061.00
1.E+01Figure 14 Elapsed Computation Time for MGCG Solver, Weak
Scaling up to 2,048 nodes, Max. Problem Size: 8,589,934,592
DOF, HB 416, “-d”: Double Precision (FP64), “-s”: Single
Precision (FP32)

1.E+02

1.E+03
1.E+04
Core #

1.E+05

1.E+06

Figure 16 Elapsed Computation Time for MGCG Solver, Weak
Scaling up to 2,048 nodes, Max. Problem Size: 8,589,934,592
DOF, Comparison of HB 416 and HB 88, “-d”: Double
Precision (FP64), SCS-b (C==8), “-s”: Single Precision(FP32), SCS-b (C==128)

Figure 15 compares double precision (FP64) and single
precision (FP32) for both of SCS-b (C==8) and SCS-b

7

Joino Usage/Research Cenoer for Inoerdisciplinary Large-scale Informaoion Infrasorucoures
Final Reporo for JHPCN Joino Research of FY 2020
Finecollection of cardiac cells. Frontiers in Physics, 8:

Level=1

579461, 2021, DOI: 10.3389/fphy.2020.579461
Level=2[2] Proceedings of International Conferences
Level=m-3(Refereed)

Level=m-2

Coarse• Communication overhead
could be reduced• Coarse grid solver is more
expensive than original
approach.
• If process number is larger,
this effect might besignificant[2] K.H. Jæger, K.G. Hustad, X. Cai, A. Tveito.
Operator Splitting and Finite Difference Schemes

Coarse grid solver ona single MPI Process(multi-threaded,further MG)

for Solving the EMI Model. Book chapter in“Modeling Excitable Tissue: The EMI Framework”,
pages 44-55, 2021, Springer, DOI: 10.1007/978-3030-61157-6_4

Figure 17 Procedures of coarse grid aggregation (CGA), where
information of each MPI process is gathered in a single MPI
process for computation at level=m-2 [14][3] K.G. Hustad, X. Cai, J. Langguth, H. Arevalo,
Efficient simulations of patient-specific electricalheart activity on the DGX-2. Poster presented at the6.

Progress during FY2020 and Future

GTC-2020 Conference.

Prospects[4] J. Langguth, N. Gaur, H. Arevalo, C. Jarvis, N.

Following the original 2-year project plan, we have

Altanaite, Q. Lan, X. Cai. Towards detailed Organ-

carried out two activities in FY2020: (1) Single-node

Scale Simulations in Cardiac Electrophysiology.

optimization of an explicit-method based simulator; (2)Poster presented at the GTC-2020 Conference.

Scale-out optimization of the explicit-method based[5] Nakajima, K., Gerofi, B., Ishikawa, Y., Horikoshi,

simulator, and (3) Preliminary Works on Mulgirgrid

M., Efficient Parallel Multigrid Solver on Intel Xeon

Solver. In (1) and (2), we have achieved close-to-

Phi Cluster, IXPUG (Intel Extreme Performancemaximum single-node performance, thanks to SIMD

Users Group) HPC Asia 2021, 2021

vectorization of the ODE computation and an effective[6] Nakajima, K., Ogita, T., Kawai, M., Efficient Parallel

OpenMP+MPI parallelization of the PDE computation.

Multigrid Methods on Manycore Clusters with

The multi-node scale-out optimization has had a good

Double/Single

start, but will require more detailed profiling to identify

PrecisionComputing,IEEE

Proceedings of the 16th International Workshop on

potential optimization opportunities. This activity (1) and

Automatic Performance Tuning (iWAPT 2021) in(2) will be continued in FY2021, while we will also carry

conjunction with 35th IEEE International Parallelout the other two planned research activities: (4)and Distributed Processing Symposium (IPDPS

Development of a new, implicit-method based simulator

2021), 2021（in press）(based on (3) in FY.2020); (5) Large-scale, realistic
simulations of cardiac electrophysiology.

[3] International

conference

Papers(Non-refereed)7.

List of Publications and Presentations[4] Presentations at domestic conference (Non-[1] Journal Papers (Refereed)refereed)[1] K.H. Jæger, K.G. Hustad (+), X. Cai (+), A. Tveito.

[5] Published library and relating data

Efficient numerical solution of the EMI model[6] Other (patents, press releases, books and so

representing the extracellular space (E), cellon)membrane (M) and intracellular space (I) of a8