Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019jh190040

Physiologically realistic study of subcellular calcium dynamics withnanometer resolution
Kengo Nakajima（The University of Tokyo）
Abstract
This project has the aim of combining advanced mathematical modeling and verylarge-scale numerical simulations for studying subcellular calcium dynamics, which
are of vital importance for the function of the heart. By incorporating realistic
geometries of the calcium release units and their microanatomical structures, thenumerical experimentation that is enabled by this joint Japan-Norway JHPCN
project can improve the realism of state-of-the-art mathematical models of
subcellular calcium signaling, thereby promoting detailed studies of how diseasedriven structural changes affect the excitation-contraction coupling in the heart. At
the same time, the associated challenge of computational capacity and efficiency isaddressed by algorithmic improvements and optimization of the involved parallel
software on various coding levels. The obtained results with respect to cardiac
electrophysiology include a quantitative understanding of the impact of sarcoplasmic
reticulum load on the process of calcium induced calcium release, validating themathematical model adopted. Moreover, an old parallel simulator which existed prior
to this project has undergone a substantial restructuring and improvement process,including reduction of memory footprint and traffic, explicit code vectorization,
optimization of MPI communication and mixed OpenMP-MPI programming. Finally,effects of pipelined CG methods were evaluated on the Oakforest-PACS system.

1.

Basic Information

Masatoshi Kawai (U Tokyo): Numericalalgorithms and parallel programming.



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



Chad Jarvis (Simula/Norway): Codeparallelization, profiling and optimization, as
well as running simulations.



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

Jonas van den Brink (Simula/Norway):
Preparation of geometries and physiologicalparameters for subcellular simulations.

2.

Purpose and Significance of Research(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 andparallel programming.
With respect to HPC, this project aims to

Xing Cai (Simula/Norway): Numericalalgorithms, code parallelization andoptimization, as well as project coordination
together with Prof. Nakajima.

greatly improve an experimental subcellular
simulator previously cevelopec by the
Norwegian partner (Simula). The olc simulator

Glenn Terje Lines (Simula/Norway): Cardiac
electrophysiology, mathematical modelingand running simulations.

was inefficient cue to inefficient cata
structures, performance unfriencly loop nests

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

anc suboptimal parallelization. It thus urgently

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

With an efficient simulator (cevelopec in thisneecec coce restructuring anc optimization.
project) to execute very large-scale simulations,1

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019

this project will consolicate a multi-scalealso provice new possibilities of in-situ huge-mathematical mocel that gives ascale cata analysis.

physiologically accurate cescription of healthy
anc pathological calcium releases, thereby

4.

acvancing the scientific uncerstancing of

Outline of Research Achievements up to
FY2018subcellular calcium cynamics.
OpenMP parallelization
For this FY, the work relatec to further

The existing olc simulator hac only MPI

optimizing the subcellular simulator will also

parallelization implementec. The lack of sharec

procuce new knowlecge about efficiently cocing

memory-basec parallelization (via threacs) is a

anc parallelizing multiple inter-tanglec stencil

clear obstacle to fully utilizing the Oakforest-

computations for the Knights Lancing anc

PACS anc Oakbricge-CX systems. Our first coce

Skylake architectures. Consequently, large-

improvement was thus focusec on incorporating

scale simulations of subcellular calcium

OpenMPhancling, using realistic geometries anc

computational tasks. We also painstakingly

cistributions of calcium release units, can beremovec pitfalls such as race concitions anc

carriec out to further valicate thefalse sharing. The largest benefit of the OpenMP

mathematical mocel, as well as quantifying theversion is avoicing the overheac cue to explicit

impact of cisease-criven structural changes oninter-process MPI communication. For example,

the excitation-contraction coupling in the heart.

for a 3D solution comain of size (2µm)3, using a

Moreover, experience will be sought on using

168x168x168 computational mesh, the OpenMP-

high-speec file cache systems for in-situ cata

enablec version neecs 64.05 seconcs (256analytics.

threacs on one Knights Lancing noce), in

parallelization

into

the

variouscomparison with 79.70 seconcs neecec by the olc

3.

pure-MPI version (256 MPI processes).

Significance as JHPCN Joint Research
Project

The significance of this JHPCN joint research

Reducing memory footprintproject has two aspects. First, UTokyo has

The main computational kernel for simulating

worlc-leacing expertise in implementing anc

subcellular calcium signaling is a set of 3D

optimizing acvancec numerical coce for real-

ciffusion-reaction equations. To mocel theworlc applications to run on cutting-ecge

realistic compartments anc geometries insicesupercomputers. Such hancs-on experience inthe calcium release units, cifferent ciffusion

supercomputing is lacking for the Norwegian

constants are valic in cifferent physiologicalpartner. Seconc, the Oakforest-PACS anc

regions. A simple strategy, which was usec by

Oakbricge-CX systems are of a suitable size forthe olc simulator, is to acopt a variable ciffusion

achieving the ambitious goal of this project,coefficient throughout the comain. This strategy

whereas access to worlc-leacing

is easy to coce but has an extensive memory

supercomputers has tracitionally been very

overheac in that each ciffusive species will neec

scarce for the Norwegian partner. The high-three accitional 3D arrays. In orcer to recuce

speec file cache systems available at UTokyo

the memory footprint, a lookup table approach2

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019

was acoptec. The lookup table, which is of a tiny

Simulations involving many CRUs

size, can be usec to accurately incorporate theWe also stuciec the impact cue to cisease-criven

flux/no-flux interactions between each pair ofchanges of the microanatomical structuresneighboring computational voxels, therefore

insice the carciac cells. For example, chronic

completely eliminating the accitional 3D arrays.

heart failure can leac to a cisintegration of theTake for instance the mesh of size 168x168x168.

calcium release units (CRUs), in terms of aFour ciffusive species woulc have neececchange in the number, cistribution anc calcium

4x3x168x168x168x8=455MB

coublerelease channels. For this research topic, it isprecision) for extra memory usage, which canvery important to acopt simulations that involve

now be savec.

many CRUs with realistic geometries anc(inmicroanatomical cetails. Simulations involving
Explicit code vectorization with AVX-512

93 CRUs that were cone on Oakforest-PACS

The “lookup table” approach, however, requires

showec that a calcium wave can inceec occur

elaborate cocing. It incuces more instructions to

uncer pathological concitions.

be carriec out on the processor level. To offset
this potential performance cisacvantage, we

5.

explicitly vectorizec the ciffusion anc reaction

Detailsof

FY2019Research

Achievements

computations using AVX-512 intrinsics. (This is

Overview

becausecompiler-

New efforts have been investec in FY2019 toenablec auto vectorization is not sufficiently

further optimize the parallel simulator ofeffective.) For example, on a 672x672x168 mesh,subcellular calcium hancling. The concrete

using 256 MPI processes, the time neecec by the

topics that have been stuciec incluce

experiments

show

that

manually vectorizec “lookup table” version is26.1 seconcs (for 100 time steps), versus 41.7



seconcs neecec by the manually vectorizec“coefficient array” counterpart.

Recucing OpenMP overheac;
Appropriate choice of number of
OpenMP threacs per MPI process.


Quantitative study of the impact of SR load

Improving overlap of MPI
communication with computation;

Many factors can affect the calcium-incucec

calcium release process, inclucing geometrical

Preliminary Investigation towarcs
Implicit Time-Marching Scheme

configuration, RyR sensitivity anc initial SR
concentration. As an example, we investigatec

Also, in preparation for running simulations on

the impact of varying the SR concentration (alsothe Oakbricge-CX system, the cocecallec SR loac). Thanks to the new simulator

optimization work cone in FY2018 has been

anc the Oakforest-PACS system, we were able

tunec anc measurec on a 4-socket server withto carry out thousancs of simulations neecec to

Xeon Skylake Golc 16-core CPUs, which is

procuce a quantitative analysis about how cells

available at Simula. Table 1 confirms that the

behave cifferently uncer cifferent SR loacs.

coce optimization solutions (e.g. using AVX-512

These simulations are necessary for valicating

intrinsics) acoptec for the Oakforest-PACS

the mathematical mocel.

system also work on the Xeon architecture.
3

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 20194

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019

In Table 1, CA cenotes the original “coefficient

Scalability of Krylov subspace methocs, such asarray” approach, whereas LUT cenotes the newCG“lookup table” approach. The subscripts autosynchronization steps that arise in cot-procucts

anc man cenote, respectively, compiler-enable

anc norm calculations on parallel machines. In

vectorization anc manual vectorization using

this work, a mocifiec preconcitionec CG methoc

AVX-512 intrinsics. Moreover, some snapshots

is presentec that removes the costly globalof a multi-CRU simulation are inclucec.

synchronization steps from the stancarc CG

methoc,suffersfromcostlyglobalalgorithm by only performing a single nonKrylov Iterative Solvers
Becauseweareblocking recuction per iteration basec on thesolving

nonlinearwork by P. Ghysels et al in 2014 [P. Ghysels, W.

phenomena, we acoptec fully explicit time-

Vanroose, Hicing global synchronization latency

marching as the scheme for integration in timein the preconcitionec Conjugate Gracient

cirection. Generally, explicit time-marchingalgorithm, Parallel Computing 40-7, 224-238,

requires huge number of time steps cue to2014

constraint of CFL concitions. We are also

communication phase can be overlappec by theconsicering

schemematrix–vector procuct, which typically only

where longer time steps are allowec. We neec torequires local communication. The resulting

solve large scale linear equations if we acopt

algorithm will be referrec to as pipelinec CG.

implicit time-marching scheme. In FY.2019, we

These methocs are comparec to a recently

evaluatec performance of preconcitionec Krylov

proposec asynchronous CG algorithm by Gropp

Iterative solvers on the Oakforest-PACS (OFP)[Ghysels 2014].

implicit

very

time-marchingsystem. In the present work, we solvec largematrices cerivec from 3D FEM applications for
solicmechanics[5].

Linearequations are solvec by Conjugate Gracientmethoc (CG) preconcitionec by SGS (Symmetric
Gauss-Seicel).

This

globalFigure 2: Algorithm of Preconcitionec CG
Methoc (PCG) [Ghysels 2014]z
Uniform Distributed Force inz-direction @ z=Zmax

In the original CG, results of cot-procucts are

Uy=0 @ y=Ymin

usec just after they are calculatec, as shown inlines 4-5 (), anc lines 8-9 () in Figure 2.

Ux=0 @ x=Xmin

Therefore,(Nz-1) elements
Nz nodes2014]].

1: r0:= b-Ax0; u0:= M-1r0; p0:= u0
2: for i=0 … do3:
s:= Api4:a:= (ri,ui)/(s,pi)5:xi+1:= xi + api
6:ri+1:= ri- as7:
ui+1:= M-1ri+18:
b:= (ri+1,ui+1)/(ri,ui)
9:
pi+1:= ui+1 + bpi
10: end do

scale linear equations with sparse coefficient
static-linear[Ghyselsalgorithm

issignificantly

affectec by the communication overheac in the(Ny-1) elements
Ny nodescot-procuctsyUz=0 @ z=Zmin

this

using

MPI_Allreduceonmassively parallel supercomputers.

(Nx-1) elements
Nx nodesx

In [Ghysels 2014], they proposec Pipelinec CG

Figure 1: 3D FEM Appl. for Solic Mechanics [5]methoc (Fig.3), where procecures of expensive
5Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019

computing (e.g. SpMV, Preconcitioning) are

In the previous work [塙敏博，中島研吾，大島聡

insertec just after cot-procucts. In Fig.3, cot-

史，星野哲也，伊田明弘，パイプライン型共役勾

procucts are calculatec on line-3 anc line-4, anc

配法の性能評価，情報処理学会研究報告（2016-

before the results are usec in line-7 anc line-8,

HPC-157-6 ） ， 2016 [Hanawa 2016]], we

preconcitioning (line-5) anc SpMV (line-6) are

implementec these algorism to Reecbush-U in

concuctec. In [Ghysels 2014], they introcucec

ITC/U.Tokyo using up to 384 noces (12,288

MPI_Iallreduce for Asynchronous Collective

cores). Figure 5 shows the results for strong

Communication which is supportec in MPI 3 orscaling.

later,

insteaccommunications

ofby

MPI_Allreduce.

Thus,

MPI_Iallrecuce

are

10000

8000

overlappec with computation on line-5 anc lineSpeed-Up

6 in Figure 3. Figure 4 cescribes Gropp’s
algorithm [Ghysels 2014], where asynchronous
collective

communication

coulc

Alg.1
Alg.2
Alg.3
Alg.4
Ideal

overlap

communications anc computations for line 3/4

6000

4000

2000

anc line 9/10.
00

1: r0:= b-Ax0; u0:= M-1r0; w0:= Au0
2: for i=0 … do
3:gi:= (ri,ui)
4:d:= (wi,ui)
5:
mi:= M-1wi
6:ni:= Ami
7:
if i>0 then8:bi:= gi/ gi-1; ai:= gi/ ( d - bi gi/ai-1)
9:
else
10:
bi:= 0; ai:= gi/d
11:
end if12:
zi:= ni + bi zi-113:
qi:= mi + bi qi-1
14:
si:= wi + bi si-115:
pi:= ui + bi pi-1
16:
xi+1:= xi + ai pi
17:
ri+1:= ri - ai si18:
ui+1:= ui - ai qi
19:
wi+1:= wi - ai zi
20: end for

2048

4096

6144

8192 10240 12288

CORE#

Figure 5: Results for Strong Scaling on
Reecbush-U System using up to 384 noces(12,288 cores), Alg.1: Original CG, Alg.2:
Chronopoulos/Gear, Alg.3: Pipelinec CG, Alg.4:
Gropp’s CG [Hanawa 2016]
While performance of original algorithm (Alg.1)
is saturatec at 8,192 cores anc at 12,288 cores,
Alg,3 (Pipelinec CG (Fig.3)) anc Alg.4 (Gropp’s
CG (Fig.4)) provice excellent scalability.
These results in Fig.5 were obtainec in July

Figure 3: Pipelinec CG [Ghysels 2014]2016 using Intel MPI Library 2015 or 2016. Intel1: r0:= b-Ax0; u0:= M-1r0; p0:= u0; s0:= Ap0;
g0:= (r0,u0)
2: for i=0 … do
3:d:= (pi,si)4:
qi:= M-1si
5:ai:= gi/d6:
xi+1:= xi + ai pi
7:ri+1:= ri - ai si8:
ui+1:= ui - ai qi
9:gi+1:= (ri+1,ui+1)10:
wi+1:= Aui+1
11:
bi+1:= gi+1/ gi
12:pi+1:= ui+1 + bi+1 pi
13:si+1:= wi+1 + bi+1 si
14: end forCompiler anc Intel MPI Library have been
improvec significantly in recent years, therefore
performance of MPI_Allreduce is much better
now. If the performance is evaluatec using IntelMPI Library 2017 or later, cifference among
Alg.1, Alg.3 anc Alg.4 is small. Moreover, Alg.3
anc Alg.4 are rather slower because they neec
more computation comparec to Alg.1. Generally,
it is cifficult to observe effects of asynchronous

Figure 4: Gropps’a CG [Ghysels 2014]

collective communications without harcware
6

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019support, such as communication assist cores on

6.

Fujitsu’s FX100 system.

Progress during FY2019 and Future
Prospects

In Intel MPI Library 2019, performance ofIn FY.2019, we have evaluated performance ofAsynchronous

Progress

Threadbeen

pipelined CG methods using up to 4,096 nodes ofsignificantly

improved,

whereefficient

OFP system with capability of Asynchronous

asynchronous collective communications are

Progress Thread in Intel MPI Library 2019.

supported

Therefore,

Pipelined CG based on Gropp’s algorithm was

Communication-Computation Overlapping will

40% faster than the original CG method at 4,096

be improved by Asynchronous Progress.

nodes (262,144 cores) of OFP.

Figure 6 and Table 2 show the performance of

Due to the coronavirus epidemic, the research

Alg.1, Alg.3 and Alg.4 in Fig.5 and [Hanawa

activities of the project were seriously disrupted

2016] using up to 4,096 nodes (262,144 cores) ofduring the last two months of FY2019.

the Oakforest-PACS system [5]. Performance ofSpecifically, the in-situ data analytics part of theAlg.4-IAR (Gropp’s CG with MPI_Iallreduce)

work didn’t get started. Moreover, the originallyattained 40% speed-up compared to Alg.1-AR

planned(Original CG) at 4,096 nodes of OFP.

simulations was only partially initiated.

bysoftware.

has

work

on

runninghuge-scale

It is envisioned that our research effort will
restart, as soon as the coronavirus situation
alleviates to allow an appropriate work mode,for continuing the work on running large-scale
simulations. The simulation results will beanalyzed so that they will constitute the main
physiological content of a journal submission.
The HPC findings associated with the large-

Figure 6: Results for Strong Scaling on
Oakforest-PACS (OFP) System using up to
4,096 nodes (262,144 cores), Alg.1-AR: Original
CG, Alg.3-AR/3-IAR: Pipelined CG, Alg.4-AR/4IAR: Gropp’s CG, AR: MPI_Allreduce, IAR:MPI_Iallreduce [5]scale simulations will be collected and drafted asanother journal submission.

Table 2: Results for Strong Scaling on OakforestPACS (OFP) System using up to 4,096 nodes(262,144 cores), Speed-up of Alg.3-IAR and
Alg.4-IAR against Alg.1-AR at each core/node
number
Node # Core #Alg.3-IAR
Alg.4-IAR
128
8,192
1.03
1.06
256
16,384
1.13
1.21
512
32,768
1.05
1.20
1,024
65,5361.06
1.25
2,048 131,072
1.08
1.34
4,096 262,144
1.19
1.38

7.

List of Publications and Presentations(1)

Journal Papers (Refereed)(2)

Proceedings of International Conferences(Refereed)[1] Chad Jarvis, Glenn Terje Lines, Johannes
Langguth, Kengo Nakajima, Xing Cai.
Combiningalgorithmic

rethinking

andAVX-512 intrinsics for efficient simulationof7

subcellular

calciumsignaling.

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019(4) Presentations at domestic conference (Non-Proceecings of ICCS 2019 Conference, Junerefereed)

12-14, 2019, Faro, Portugal.

[4] Xing Cai. Heterogeneous computing for[2] Johannes Langguth, Hermenegilc Arevalo,
Kristian

Hustac,

Xing

Cai.

cardiac electrophysiology. Invitec keynote

Towards

talk at PREAPP workshop on Efficientdetailed real-time simulations of cardiac
arrhythmia.

ProceecingsofFrameworks

theSeptember

8-11,

anc

Data-

26, 2019, Tromsø, Norway.

2019,[5] Masashi

Singapore.
(3) International

Compute-intensive Computing (EFFECT), April 25-

International Conference in Computing in
Carciology,

forHorikoshi,

Kengo

Nakajima,

Balazs Gerofi, Yutaka Ishikawa, Parallel
conference

PapersPreconcitionec(Non-

Iterative

Solverson

Oakforest-PACS, 2019 年並列／分散／協調

refereed)[3] Nakajima, K., Parallel Multigric with

処理に関する『北見』サマー・ワークショッ

Acaptive Multilevel hCGA on Manycore

プ（SWoPP 北見 2018），日本応用数理学会

Clusters,

Extreme-Scale/Exascale「行列・固有値問題の解法とその応用」研究

Applications China, Japan, Worlc, ISC

部会（MEPA）（北見，2019 年 7 月 25 日）High

Performance

2019,

Frankfurt,(5) Other (patents, press releases, books and so

Germany, 2019 (Invitec Talk)on)
N/A

8