5
SC’11 Special Issue

SC’11 Special Issue
Gordon Bell Prize

Special Achievements in Scalability and Time-to-Solution
Honorable Mention

Graph500 Ranking No.3
3 Technical Papers

http://www.gsic.titech.ac.jp/TSUBAME_ESJ

SC’11 Special Issue

03

ACM Gordon Bell Prize:
Special Achievements in Scalability and Time-to-Solution& SC'11 Technical Paper

Peta-scale Phase-Field Simulation for Dendritic
Solidification on the TSUBAME 2.0 Supercomputer
Takashi Shimokawabe
Takayuki Aoki
Tomohiro Takaki
Akinori Yamanaka
Akira Nukada
Toshio Endo
Naoya Maruyama
Satoshi Matsuoka

10

ACM Gordon Bell Prize: Honorable Mention

Large scale biofluidics simulations on TSUBAME2
Massimo Bernaschi
Massimiliano Fatica

Mauro Bisson
Toshio Endo
Satoshi Matsuoka
Simone Melchionna

Sauro Succi

17

21
26

Graph500 Ranking No.3

Graph500 Challenge on TSUBAME 2.0
Toyotaro Suzumura

Koji Ueno

SC'11 Technical Paper

Physis: A High-Level Stencil Framework for
Heterogeneous Supercomputers
Naoya Maruyama

Tatsuo Nomura

Kento Sato

Satoshi Matsuoka

SC'11 Technical Paper (Achieving a Perfect Score)

FTI: high performance Fault Tolerance
Interface for hybrid systems
Leonardo Bautista-Gomez
Dimitri Komatitch
Naoya Maruyama
Franck CappelloSatoshi Matsuoka
Takeshi Nakamura

Seiji Tsuboi

Peta-scale Phase-Field Simulation for Dendritic
Solidification on the TSUBAME 2.0 Supercomputer
Takashi Shimokawabe*
Akira Nukada**

Takayuki Aoki**
Toshio Endo**

Tomohiro Takaki***Naoya Maruyama**

Akinori Yamanaka****Satoshi Matsuoka*** Interdisciplinary Graduate School of Science and Engineering, Tokyo Institute of Technology** Global Scientific Information and Computing Center, Tokyo Institute of Technology*** Graduate School of Science and Technology, Kyoto Institute of Technology**** Graduate School of Engineering, Tokyo Institute of Technology

To realize a low-carbon society, it is indispensable to develop new strong and light materials. The material
microstructure controlling the mechanical property is determined in the solidification process. The phase-fieldmodel is a meso-scale physical theory describing phase transformations such as solidification. Realistic prediction ofmaterial solidification requires a millimeter-scale simulation heretofore impossible with supercomputers. The phasefield equation was discretized by using the finite-difference method, and a simulation code in CUDA was developed.
A large-scale solidification simulation of an aluminum-silicon alloy was carried out with 4,000 GPUs on TSUBAME 2.0
and achieved 2.0 petaflops in single precision. The ACM Gordon Bell Award - Special Achievement in Scalability andTime-to-Solution was awarded for this achievement in 2011.

Introduction

1

The development of new strong and light materials contributes
greatly to transportation systems with high fuel efficiency andrealization of a “low-carbon” society. The strength of a metal strongly
depends on its microstructure, which is, in turn, determined by thesolidification process of the metal (see Fig .1). However, determining
material properties (such as strength) requires millimeter-scale
macroscopic study.
The phase-field model [1] is a physical theory describing
the evolution of complicated material morphologies on the “meso-

Figure 1 Images of a material microstructure

scale,” namely, between the molecular scale and the macroscopic
scale of materials. It is a powerful numerical tool for studying the
dynamics of phase transformations such as solidification. The
derived equations are partial differential equations in time and2

Phase-field model

space and often discretized by FDM (finite-difference method) or
FEM (finite element method). The phase-field equations include
many complex nonlinear terms, and the amount of computations

T h e p ha s e - f i e l d m o d e l is d e r i ve d f ro m n o n - e q u il ib r ium

per mesh or element becomes large compared to normal stencil

statistical mechanics and can describe meso-scale phenomena

applications. Furthermore, the phase-field model assumes athat occur on a mid-scale between the molecular and macro

narrow-thickness interface between liquid phase and solid phase,scales. The phase-field model introduces a “continuous-orderso the time integration requires high spatial resolution and shortparameter,” i.e., a phase-field variable, to describe whether atime step, and phase-field simulations remain restricted to two-material is solid or liquid. The phase fieldor small three-dimensional computations due to their large

in the liquid phase and 1 in the solid phase. Interfaces between

computational cost.

solid and liquid phases are treated as diffuse ones described

In this study, we carry out a large-scale phase-fieldby localized regions, where the parameter changes smoothly

simulation for the aluminum-silicon dendritic growth during

bet ween the t wo fixed values. The inter face positions are

directional solidification on the GPU supercomputer TSUBAME

represented by

2.0. The simulation code is developed in CUDA, and almost all themethod can describe the locations of the interfaces without

capability of TSUBAME 2.0 is used for the simulation, which uses ausing a traditional interface-tracking method and carry out thedomain decomposition and inter-node communications with ansame calculations throughout a computational domain.

MPI library. Such large-scale phase-field simulation has not been
done before and will have a big impact on material science.

3

takes a value of 0

. Thanks to this approach, the phase-field

The dendritic solidification process in a binary alloy issimulated by solving the phase-field equation and the diffusion

SC’11 Special Issue

equation of the solute concentration. The time integration of thephase field is described by the following equation, which takes into
account the anisotropy of interfacial energy:

Single-GPU implementation

3

T h e t i m e i n te g r a t i o n o f t h e p h a s e f i e l d a n d t h e s o l u teconcentration given by Equations (1) and (2) are carried out by thesecond-order finite-difference scheme for space and the first-orderforward Euler-type finite-difference method for time on a threedimensional regular computational grid. The simulation code iswritten in CUDA to run on GPUs. All variables are allocated on GPU
where

,

,

,

, and

represent the gradientcoefficient that describes interface anisotropy, the mobility of thephase field, the height of the potential energy barrier, the entropyof fusion and undercooling, respectively. The two functions
andare given as

and ,video memory (also called global memory in CUDA), which virtuallyeliminates all the CPU-GPU data transfers through a PCI-Express
bus during simulation runs. The spatial patterns of the neighborpoints of phase fieldand solute concentrationthat arerequired to solve the discretized governing equations at thecenter point (i, j, k) of the grid are, respectively, shown in Figures 2

respectively.
The time integration of the solute concentration isderived by the following diffusion equation:(a) and (b). In one time step , 19 neighbor elements ofneighbor elements ofand sevenare used with the governing equations to

update the phase field value and concentration value at the center
point; thus, 26 elements must be read from the memory and two
updated values must be written back to the memory for each pointof the grid.

where

and

are the diffusion coefficients in the bulk solidand liquid phases , respectively. Solid concentration
concentration

satisfy

and liquid.

(a) Stencil of phase-field variable.

(b) Stencil of concentration variable.
Figure 2 Spatial access patterns of neighbor points
required for time integration.

4

ACM Gordon Bell Prize: Special Achievements in Scalability and Time-to-Solution & SC'11 Technical Paper

Peta-scale Phase-Field Simulation
for Dendritic Solidification on the TSUBAME 2.0 Supercomputer

The computation in the phase-field model involves alarge number of memory accesses. It is thus highly effective to
reduce access to the global memory. In GPU computing, a key

Optimization of multi-GPU computing

factor in performance improvement is how threads and blocks are

4

assigned to actual calculation. The elements of the computational
domain of size nx × ny × nz are calculated as follows. First, asWe describe our strategies of multi-GPU computation. Usingshown in Figure 3, the given domain is divided into pieces of size

multiple GPUs is necessary for high-performance computing of

64 × 4 × 32. Our kernel function is invoked on a GPU, and eachlarge-size problems. In addition to the intrinsic parallel-processing

thread block (which has 64 × 4 × 1 threads) handles each piece.

capabilities of GPUs, it is necessary to do parallel computing by

Each thread calculates 32 elements by marching in the z direction;using a large number of distributed GPUs. We decompose thea thread that corresponds to (i, j) updates grid points (i, j, k), where

whole computational domain in both the y- and z-directions (i.e.,k varies from k0 to k0 + 31, and k0 is a multiple of 32 (i.e., k0 = 0, 32,

2D decomposition) and allocate each subdomain to one GPU.

64, ..., nz − 32).

When a thread computes a point in the k + 1-th

We have chosen this method because 3D decomposition, whichplane, some elements to be referred to in this computation haveis better in terms of reducing communication amount, tends toalready been accessed by that thread in the previous k-th planedegrade GPU performance due to complicated memory accesscomputation. Such data are reused by holding them in registers,

patterns in the y-z plane for data exchanges between the GPU

thereby reducing global memory accesses. On the other hand, on-

and CPU. Similar to conventional multi-CPU implementations,chip shared memory is not used to store data shared by several

multi-GPU implementation requires boundary-data exchanges

neighbor threads; instead, the L1/L2 cache available on Fermi GPUs

between subdomains. Because a GPU cannot directly access theincluding M2050 is relied on.

global memory of other GPUs, host CPUs are used as bridges
for data exchange. For inter-node cases, this data exchange iscomposed of the following three steps: (1) data transfer from the
GPU to the CPU by using CUDA APIs, (2) data exchange between
nodes with the MPI library, and (3) data transfer back from the
CPU to the GPU by using CUDA APIs.
Unlike in conventional multi-CPU computation, in multiGPU computing, data-communication time with neighbor GPUs

Figure 3 Layouts of CUDA threads and blocks.

is not ignored in the total execution time, especially in large-scalecomputation, since the performance of the GPU is higher thanthat of the CPU. To achieve higher performance, it is important tohide communication overheads by overlapping communication
with computation. In this study, we implement the following threemethods for multi-GPU computing: (a) GPU-only method, (b)Hybrid-YZ method, and (c) Hybrid-Y method.
(a) GPU-only method: This basic method (which is used as areference) computes all subdomains without applying any specificoptimization. After that, the boundary regions of subdomains areexchanged between GPUs by using the above three steps.
(b) Hybrid-YZ method: By dividing each subdomain into five
regions, which are two y-boundaries, two z-boundaries and an
inside region, we can overlap communication for the boundary
data exchange with the computation of the inside region. A similar
overlapping technique has been described in previous reports [2].
In this technique, separate GPU kernels are invoked for separate
regions. Instead of dividing a GPU kernel, in this study, we let CPU
cores compute four y- and z-boundary regions, while the inside5

SC’11 Special Issue

region is computed by GPU. We observed that CPUs often become

Since each node of TSUBAME 2.0 has three GPUs and 12

bottlenecks by taking longer time to do computation of all

CPU cores (two 6-core Xeon CPUs), we assign four cores to each ofboundaries and communications than the GPU computation time.

the three GPUs. Each subdomain is thus cooperatively computed(c) Hybrid-Y method: This method is a slightly modified methodby a single GPU and four CPU cores. Multi cores are exploited byfrom the (b) method. It reduces CPU loads by assigning only theusing OpenMP.

y-boundaries to CPUs instead of both y- and z-boundaries. The
GPU can compute the z-boundaries effectively, since it accesses
consecutive addresses of the global memory.

Figure 4 Schematic diagramof the GPU-only method.

Figure 5 Schematic diagram of the Hybrid-YZ method.

Figure 6 Schematic diagram of the Hybrid-Y method.

6

ACM Gordon Bell Prize: Special Achievements in Scalability and Time-to-Solution & SC'11 Technical Paper

Peta-scale Phase-Field Simulation
for Dendritic Solidification on the TSUBAME 2.0 Supercomputer

Performance of multi-GPU computing

5

To measure the performance of the GPU computation,
we count the number of floating-point operations in the C/C++based phase-field simulation code by running it on a CPU witha performance counter provided by PAPI (Performance API). The
obtained count and GPU elapsed time are used for evaluating the

We show the performance of (a) the GPU-only method, (b) the

performance of the GPU computing. Figure 8 shows the results

Hybrid-YZ method, and (c) the Hybrid-Y method, all using multiple

on strong scaling, which show variation of the performance

M2050 GPUs of the GPU-rich supercomputer TSUBAME 2.0 at thewith the number of GPUs for a fixed total mesh size. This figure

Global Scientific Information and Computing Center of the Tokyo

compares the strong-scaling characteristics of the three proposed

Institute of Technology. In this study, we simulate solidification of

methods i.e., (a), (b), and (c). We perform simulations in single

aluminum–silicon (Al–Si) alloy.

precision for three different mesh sizes: 5123, 10243, and 20483. It is

An image of the solidification process of the alloy observed

clear from the figure that the overlapping methods, (b) Hhybrid-at the world’s largest synchrotron radiation facility, SPring- 8, by

YZ method and (c) Hhybrid-Y method, both work effectively in

Hideyuki Yasuda, a professor at Osaka University, and his researchhiding communication overhead as expected, resulting in angroup, is shown in Figure 7 (a) . Results of the phase-field simulation

improvement of overall performance when a small number ofusing a 4096 × 128 × 4096 mesh are shown in Figure 7 (b) . Although

GPUs are used. When the number of GPUs is larger, the volume

the alloy systems in the images differ, the growth process obtainedthat each GPU handles becomes smaller, and the percentage

by the simulation agrees well with the observed one.

of the boundar y region computed by the CPU in the whole
computational domain increases. Eventually, the GPU computation
is no longer able to hide the communication cost. In the case of theHybrid-YZ method, method (b), we observed that the CPUs often
become bottlenecks by taking a longer time to do their operations
than the GPU computation time especially when the number ofGPUs is larger. Since the Hybrid-Y method alleviates this bottleneck,it achieves a significantly improved performance compared withthat of the GPU-only method.

Figure 7(a) Solidification process of an alloy observed
at SPring-8 (courtesy of Professor Yasuda).

Figure 8 Results for strong scaling of multi-GPU
computation in single precision.
Next, we show the weak scaling results, which show how
the performance of the three methods varies with the number of
Figure 7(b) Solidification growth simulatedby the phase-field model using GPUs.

7

GPUs for a fixed mesh size per GPU. Each GPU handles a domainof 4096 × 160 × 128 for single precision. In demonstrating weakSC’11 Special Issue

scalability, the Hybrid-Y method has achieved extremely highFigure 11 demonstrates the dendritic growth duringperformance in single precision (Figure 9), namely, reaching

binary-alloy solidification using a mesh size of 4096 × 1024 × 4096

2.000 petaflops, specifically 1.975 petaflops by the GPU and 0.025

on TSUBAME 2.0. As an initial condition, 32 nuclei are put on the z

petaflops by CPU, for a 4096 × 6480 × 13000 mesh using 4,000= 0 plane. This simulation will allow us to clarify the mechanisms ofGPUs along with 16,000 CPU cores. This is the first peta-scale resultcompetitive growth between parallel dendrites. Evaluating theseas a real stencil application we know to date.

mechanisms is important to design and control the mechanicalproperties of solidified products.

Figure 9 Results on weak scaling of multi-GPU
computation in single precision.

Figure 10 shows that electric power consumption of
TSUBAME 2.0 when the simulation achieved 2 petaflops, namely,
44.5% of the efficiency to the peak performance, using 4,000 GPUs
along with 16,000 CPU cores. The 2-petaflops simulation consumed
electric power of 1.36 MW for all computational nodes and networks
on TSUBAME 2.0. The simulation thus achieved 1468 Mflops/W.
These results show that the simulation results were obtained bysmall electric power consumption.

Figure 11 Large-scale simulation of dendritic
solidification of aluminum–silicon alloy.

Figure 10 Electrical power consumption of
the 2-petaflops phase-field simulation
running on TSUBAME 2.0.

8

ACM Gordon Bell Prize: Special Achievements in Scalability and Time-to-Solution & SC'11 Technical Paper

Peta-scale Phase-Field Simulation
for Dendritic Solidification on the TSUBAME 2.0 Supercomputer

Summary

6[4]

T. Shimokawabe, T. Aoki, T. Takaki, A. Yamanaka, A. Nukada,
T. Endo, N. Maruyama, and S. Matsuoka : Peta-scale PhaseField Simulation for Dendritic Solidification on the TSUBAME
2 . 0 Supercomputer, in Proceedings of the 2 011 ACM/IEEEInternational Conference for High Performance Computing,A large-scale GPU simulation of the dendritic growth during binary-

Networking, Storage and Analysis, SC’11, IEEE Computer

alloy solidification was carried out on the basis of a phase-field

Society, Seattle, WA, USA, Nov. 2011.

model. An extremely high performance of 2.0 petaflops in single
precision was achieved on 4,000 GPUs of TSUBAME 2.0, in spite ofa stencil application, which is hard to extract high performance.
The computation reached 44.5 % of the peak performance andsimultaneously high efficiency from the viewpoint of electrical
power consumption. It means that we have the computational
result on TSUBAME 2.0 with smaller energy compare d toconventional supercomputing. It is concluded that a GPU
supercomputer is available for various practical stencil applications.
Acknowledgements
The peta-scale simulation was executed as a TSUBAME Grand
Challenge Program 2011 spring and we express our special
thanks to have a chance to use the whole TSUBAME resources.
This research was supported in part by Grant-in-Aid for Scientific
Research (B) 23360046 from The Ministry of Education, Culture,
Sports, Science and Technology (MEXT), two CREST projects, "ULPHPC: Ultra Low-Power, High Performance Computing via Modelingand Optimization of Next Generation HPC Technologies" and"Highly Productive, High Performance Application Frameworks for
Post Petascale Computing" from Japan Science and Technology
Agency (JST), and JSPS Global COE program "Computationism as aFoundation for the Sciences" from Japan Society for the Promotionof Science (JSPS).
References[1]

R . Ko b ay ashi: M o d e l in g a n d n u m e r i c a l sim u l at i o ns o fdendritic crystal growth. Physica D, Nonlinear Phenomena,63(3-4), 410 - 423 (1993)[2] Ta k a y u k i A o k i , S a t o i O g a w a , A k i n o r i Ya m a n a k a : G P UComputing for Dendritic Solidification based on Phase-Field
Model, TSUBAME e-Science Journal, Vol.1, pp. 5 - 8, Global
Scientific Information and Computing Center, Tokyo Tech (2010
September)[3] A . Ya m a n a k a , T. A o k i , S . O g a w a , a n d T. Ta k a k i : G P U accelerated phase-field simulation of dendritic solidification
in a binary alloy. Journal of Crystal Growth, 318(1):40 - 45 (2011).
The 16th International Conference on Crystal Growth (ICCG16)/
The 14th International Conference on Vapor Grow th andEpitaxy (ICVGE14)

9

Large scale biofluidics simulations
on TSUBAME2
Massimo Bernaschi*
Mauro Bisson*
Toshio Endo**
Massimiliano Fatica*** Satoshi Matsuoka** Simone Melchionna* Sauro Succi** CNR-IAC, Istituto Applicazioni Calcolo, Consiglio Nazionale delle Ricerche, Rome, Italy** Global Scientific Information and Computing Center, Tokyo Institute of Technology*** Nvidia Corp. Santa Clara, CA, USA

Multi-scale simulations of real-life biofluidic problems entail simulating suspensions composed by hundreds ofmillions of bodies interacting with each other and with a surrounding fluid in complex geometries. One such
example is the simulation of blood flow through the human coronary arteries, with spatial resolution comparablewith the size of red blood cells, and physiological levels of hematocrit (the red blood cell volume fraction).
We developed a methodology such that hemodynamic simulations exhibit excellent scalability on the TSUBAME 2
installation, achieving an aggregate sustained performance of hundreds of Teraflops. The result demonstrates
the capability of predicting the evolution of biofluidic phenomena of clinical significance, by using a suitablecombination of novel mathematical models, computational algorithms, hardware technology and code tuning.

1

Introduction

The simulations involve up to a billion fluid nodes, embedded in abounding space of about a three hundred billion voxels, with 10-450
million suspended bodies, as shown in Fig. 1. They are performed
with the (MUlti PHYsics/multiscale) MUPHY code, which couples
Lattice Boltzmann method for the fluid flow and a specialized

Transport phenomena are ubiquitous in living systems, underlying

version of Molecular Dynamics for the suspended bodies[4,5]. The

muscle contraction, digestion, the nourishment of cells in the body,simulation achieves an aggregate performance in excess of 600

blood circulation, to name but a few. Blood is a reference biofluid,

Teraflops, with a parallel efficiency of more than 90 percent on 4000-

being the carrier of those biological components that fuel the

GPU of the TSUBAME2 system.

most basic physiological functions, such as metabolism, immune
response and tissue repair.

Our work presents a number of unique features, both atthe level of high-performance computing technology and in terms

Building a detailed, realistic representation of bloodof physical/computational modeling. The extreme complicationand the vasculature represents a formidable challenge since theof the irregular geometries demands that the workload be evenly

computational model must combine the motion of the fluid within

distributed across the pool of as many as 4000 GPU of the TSUBAME2

an irregular geometry, subject to unsteady changes in flow andsupercomputer. The resulting domain-partitioning problem, evenpressure driven by the heartbeat, as coupled to the dynamics ofat the mere level of the fluid computation, poses a formidable

red and white blood cells and other suspended bodies of biologicalchallenge. On top of this, the application adds the further constraint

relevance.

of keeping a good workload balance also for the Molecular
Large-scale hemodynamic simulations have made
[1–3]

substantial progress in recent years

Dynamics (MD) component of the multiscale methodology.

, but until now the coupling

Complex and large-scale geometries, such as the oneof fluid dynamics with the motion of blood cells and otherconsidered here, are rare in the literature. By leveraging large-scalesuspended bodies in vessels with realistic shapes and sizes, hasparallel architectures, this work demonstrates the feasibility ofremained beyond reach. Owing to non-local correlations carriedcardiovascular simulations of unprecedented size that do not rely onby the flow pressure, the global geometry plays a significant role

any geometrical regularity of the computational domain.

on local circulation patterns, most notably on the shear stressat arterial walls. Wall shear stress is a recognized trigger for thecomplex biomechanical events that can lead to atherosclerotic
pathologies. Accurate and reliable hemodynamic simulations of thewall shear stress may provide a non-invasive tool for the predictionof the progression of cardiovascular diseases.
We illustrate here the first multiscale simulation ofcardiovascular flows in human coronary arteries reconstructed from

Figure 1 Geometry of the simulated
coronary arteries, with theunderlying level of red
blood cells embedded in
the Lattice Boltzmann mesh.

computed tomography angiography. The coronary arteries form
the network that supply blood to the heart muscle and span theentire heart extension. Spatial resolution extends from 5 cm down
to 10 μ m, where a red blood cell has a diameter of about 8 μ m.

10

ACM Gordon Bell Prize: Honorable Mention

Large scale biofluidics simulations
on TSUBAME2

Multiscale biofluidics

2

enabling scale-specific hydrodynamic interactions at sub-meshspacing resolution; iv) avoiding the detailed representation of thefluid-body dividing interface.
The last two points are crucial to boost the performancesof the simulations and to guarantee numerical stability of the dual

Living matter is typically composed of two main components: amethod, in particular as related to the stiff forces exerted on the fluid

generic liquid, given by an aqueous solution, plasma, cytosol, etc.,

by the immersed bodies. The framework employs either a single

and suspended bodies, such as bioparticles, cells, proteins, DNA,or a multiple timestep algorithm to handle stiff forces. It should be

etc. , that move within the embedding solvent. We have developed

stressed that the LB method is largely tolerant towards the presence

the MUPHY[4,5] software to simulate such generic conditions.

of rapidly varying forces and allows simulating dense suspensions

As compared to the previous versions of MUPHY[4,5], thefrom creeping flow conditions up to Reynolds number in the ordercurrent simulation framework handles the numerical solution ofof 1000. This strategic asset allows covering a wide range of physicalthe motion of a generic fluid and suspended bodies. Within thisphenomena in real-world physiological conditions, enabling tovision, MUPHY represents a major effort to design a computational

reproduce highly non-local rheological response that cannot be

infrastructure for biofluidics research based on the use of a libraryassimilated to a continuum governed by constitutive relations.

of solvents and solutes. A non-exhaustive list of functionalities
includes: selection of Newtonian versus non-Newtonian rheologicalresponse; selection of dif ferent kernels for the collision, to

The LB method[10] is based on the evolution of the singlet distribution

reproduce from molecular hydrodynamics to stochastic frictional

representing the probability of finding, at mesh location x and atdynamics; inclusion of stochastic fluctuations; selection of solutes

time t, a “fluid particle” traveling with a given discrete speed. “Fluid

as disconnected particles or forming polymer and molecules;

particles” represent the collective motion of a group of physicalsuspended bodies with polymorphic hydrodynamic shapes;

particles (often referred to as populations). The hydrodynamic fluid-

selection of different body-body forces; scale-adaptive body-

body coupling is based on specific roto-translational kernels thatfluid coupling mechanisms; handling of irregular confining media

represent either rigid entities moving in the fluid as impenetrable

and tissues. Within this range of options, biofluids are modeled

bodies, soft vesicles, or a combination thereof.

in multiple ways and the multiscale/multiphysics behavior of
biological or physiological systems can be studied accordingly.

11

2-1 Lattice Boltzmann/Molecular Dynamics

The fluid-body hydrodynamic interaction is constructedthaccording to the transfer function centered on the i particle and

MUPHY leverages two computational engines. The firsthaving spherical or ellipsoidal symmetry and compact support,

one handles the motion of the generic fluid within the hydro-with a hydrodynamic shape that can be smaller than the meshkinetic formulation embodied by the Lattice Boltzmann (LB)

spacing[6]. The fluid-particle coupling requires the computation of

method: collision-driven mechanisms reproduce the dynamics

convolutions over the mesh points {x} and for each configuration ofof fluids in the continuum (in contrast to the direct macroscopic

the N suspended bodies.

description of fluid motion via the Navier-Stokes equations). The

While, in principle, the calculation of hydrodynamic

second engine handles the motion of Lagrangian bodies by using

interactions grows cubically with the system size, the concurrent

a technique that shares several technical aspects with Molecular

evolution of fluid and bodies provides a strategic algorithmic

Dynamics (MD). However, the nature of the suspended bodies

advantage. In fact, the LB method features a computational

is non- conventional and requires a substantial extension of thecost which scales linearly with the number of mesh points M,

basic MD technique. Finally, the coupling between fluid andthat is, O(M )=O(N/c ), for a fixed solute concentration c = N/M .

moving particles takes place via specifically designed kernels based

By employing the link-cell method to compute the direct forceson kinetic modeling, again significantly distinct from stresslet,

among bodies, the particle dynamics also shows O(N ) complexity.

boundary-integral and other methods based on macroscopic

Thanks to the fact that the solvent-mediated particle-particle

hydrodynamics. The end result of the computational environment

interactions are localized and explicit, the LB-MD coupling scales

is a fully time-explicit simulation technique that offers strategic

linearly with the number of mesh elements and suspended

advantages for the study of biofluids under realistic conditions, inbodies. However, hydrodynamic coupling represents the mostparticular: i) employing geometries with irregular boundaries, suchtime-consuming component of the methodology due to theas in the case of the blood vasculature; ii) describing non-trivial

large overhead arising from the O(100 ) number of mesh points

rheology, as emerging from the underlying particle dynamics; iii)

enclosed in the particle support and scattered access to memory.

SC’11 Special Issue

3

MUPHY

CPU (like MPI primitives) may be concurrent with the execution ofkernels on the GPU. As a result, the CPU should be seen as a MPI coprocessor of the GPU.
3-1 Domain decomposition

The MUPHY (MUlti PHYsics/multiscale) code was originally written

The geometry used in our simulations is highly irregular, as shown

in Fortran 90 by using MPI for the parallelization[4]. For a flexible

in Figure 1, and partitioning in subdomains handled by the available

and efficient handling of complex geometries, MUPHY makes

computing resources represents a major challenge in itself. In

use of an indirect addressing scheme that not only limits thea previous attempt to solve this issue, we used the third-party

memory requirements to the bare minimum, but also facilitates

software PT-SCOTCH, the parallel version of the SCOTCH graph/the achievement of a good load balancing and the selection of themesh partitioning tool[12], in order to distribute the computational

most suitable communication pattern depending on the platform

load in an even manner. The graph-based procedure utilizes a graphin use. MUPHY was originally developed for the IBM BlueGenebisection algorithm and it is completely unaware of the geometryarchitecture[7], a system characterized by thousands of relativelyof the computational domain. We then realized that the lack ofslow PowerPC processors connected by a high performance custom

geometrical information degrades the quality of the partitioning as

network. On BlueGene/P we achieved an excellent scalability up

the number of partitions increases, in which case the subdomains[11]

to the largest configuration available to us (294,912 cores) . The

reduce to highly irregular shapes with large contact areas each otherobtained total performance was, in that case, in the range of tens

that increase the communication overhead. An optimal solutionof TeraFlops, with limited room for further improvements since

was found by combining the graph-based partitioning with aour code cannot exploit the SIMD-like operations of the PowerPC

flooding-based approach (also known as graph-growing method)architecture (those operations require stride-one access whereas

according to the following procedure: the mesh is first partitioned,

both the LB and the MD components of MUPHY have a “scattered”

by using PT-SCOTCH, in a fixed number of subdomains (256). Thendata access pattern). In the meantime, the steadily increasing

each subdomain is further divided by using a flooding algorithm.

computing power of modern GPUs motivated us to develop a[7]version of MUPHY targeted to clusters of GPUs .

In MUPHY, the communication pattern is set up at run
time. Each task determines the neighboring tasks owning mesh

A multi-GPU code resorts to parallelism at two levels:

points that need to be accessed during the simulation, regarding

intra- and inter- GPU. Unlike the case of codes developed for clusters

the non-local streaming of populations in the LB algorithm, theof traditional multi-core systems that can be implemented either

migration of particles and calculation of inter-domain forces forin a hybrid (e.g., OpenMP+MPI) or in a simple distributed memory

the particle dynamics. During this pre-processing step, we employ

layout (counting on the avail- ability of highly efficient MPI libraries

mainly MPI collective communication primitives whereas, in thefor shared memory systems), for a multi-GPU platform, a hybrid

remainder of the execution, most communications are point-to-

paradigm is the only choice. This is not the only difference thatpoint and make use of the following scheme: the receive operations

needs to be taken into account: on a traditional multi-core platform

are always posted in advance by using corresponding non-blocking

the parallelism is limited to a few threads, at most tens on high-end

MPI primitives, then the send operations are carried out. Finally,systems. On a GPU hundreds of threads are required to keep theeach task waits for the completion of its receive operations, by using

hardware busy, so that a much finer-grain parallelization is required.

the MPI wait primitives.

Finally, albeit the combination of latest generation Nvidia GPUs and

For the Lattice Boltzmann component, only the evaluation

CUDA drivers offers the chance to copy data from/to the global (i.e.,of global quantities (e.g. , the momentum along the x, y, z directions)main) memory of GPU to/from the global memory of another GPU

is carried out by using MPI collective (reduction) primitives. Also forsitting on the same system, GPUs can not, in general, exchange data

the Molecular Dynamics component most of the communications

with each other without using the CPU.

are point-to-point but, here, the irregularity of the geometry and, asAt first glance, the passage through the CPU introduces
an overhead in the communication among the GPUs but, by using

a consequence, of the corresponding domain decomposition resultsin other issues, as detailed out in the next subsection.

the CUDA concepts of stream and asynchronous memory copies, it
is possible to overlap data transfers between GPU and CPU memory
with the execution of kernels (functions in the CUDA jargon) on
the GPU. Moreover, also the execution of functions running on the

12

ACM Gordon Bell Prize: Honorable Mention

Large scale biofluidics simulations
on TSUBAME2

3-2 Parallel Molecular Dynamics

inside the frontier cells of neighboring domains. For this reason,
the coupling is done by exploiting the cell tiling, analogously

In most parallel Molecular Dynamics applications the geometry ofto the computation of inter-domain forces. In order to computethe spatial domain is regular and simple Cartesian decompositions

the forces/torques acting from the fluid to the particles, eachare applied such that each task has (approximately) the sameprocessor exchanges particles inside the frontier cells with its

number of particles. In an irregular domain, this strategy wouldneighbors. Finally, forces and torques associated to the external

produce two distinct domain decompositions: one for the LB

particles are exchanged back with neighboring processors and,(as described above) and one for the particle dynamics. As aon the receiving side, the external contributions are distributed

consequence, the same subdomain might belong to two or moreto the frontier particles. This approach is particularly efficient asprocessors for the LB and for the MD components and consequently

compared to other strategies. For example, by exchanging thethe interaction between particles and fluid would become highly

mesh points located inside frontier cells in place of particles,

non-local with a complex and expensive communication pat- tern.

one could lower the number of exchanged data to a single

We then decided to resort to a domain decomposition strategy

communication step. However, such alternative would impose

where the MD parallel domains and the LB parallel domains

a large communication overhead since all mesh nodes near thecoincide. In this way, each computational task performs both the LB

interfaces among domains should be exchanged.

and MD calculations and the inter- actions of the particles with the

Finally, regarding the computation of the momentum

fluid are quasi-local. The underlying LB mesh serves the purposetransfers exerted from the particles to the fluid, it is carried out byof identifying particles that belong to the domain via a test of

exchanging the frontier particles as in the particle-on-fluid casemembership: a particle with position R belongs to the domain if

but without a second data exchange.

the vector of nearest integers coincides with a mesh point of thedomain. Since an even number of bodies is expected to populate
the domains, a pretty good load balancing of the MD parts is
granted by the mesh partitioning.
We have developed a novel parallelization strategy
suitable for domains with irregular geometry. Among others,a major issue regards the identification of particles that residenext to the subdomain frontiers and that interact with intra- andinterdomain particles. Our solution[8] relies on the notion of cells,
cubes with side greater or equal to the interaction cutoff, thattile the whole irregular domain handled by a computational task(see Fig.2). This representation allows the processors to perform
an efficient search of both interdomain and intradomain pairs ofparticles and to reduce data transfers by exchanging a limited
superset of both the particles actually involved in the evaluation offorces among interdomain pairs and the particles moving acrossdomains.
The final component of our multiscale application
deals with the fluid-particle coupling. Each suspended particle
experiences hydrodynamic forces and torques arising from thefluid macroscopic velocity and vorticity, smeared over a domain
made of 4×4×4 mesh points. Analogously, mesh points experiencea momentum transfer arising from the surrounding particles.
These non-local operations require multiple communication steps
such that each processor owning a given particle, exchanges
hydrodynamic quantities with the surrounding domains.
Par ti cl e - f lui d co up lin g is a n o n - l o c al o p er atio n
involving the interaction between particles and mesh nodes

13

Figure 2 Tiling of an irregular domain in external,
frontier and internal cells, and selective
exchange with neighboring domains.

SC’11 Special Issue

4

Resultsof hematocrit. The efficiency features superlinear scaling for anumber of GPUs up to 1024, whereas at the largest number ofGPUs available the efficiency slightly lowers to ≃ 80% .
The excellent scalabilit y should be ascribed to the optimalload balancing obtained with our hybrid graph partitioning/

In our benchmark s, we measure the total runtime for theflooding scheme for the multi-branched arteries. Although thesimulation together with the breakdown between computation

distribution of red blood cells over a test case of 1200 domains

and communication. All simulations include about 1 billion latticeshows a broad distribution of values, the execution time for bothsites for the fluid, within a bounding box having a total of almost

the aggregate MD and the LB component is compact. The result

300 billion nodes, and 450 million Red Blood Cells (RBC). We choose

is a direct consequence of the minimal extension of contact

the fluid-body coupling parameters such that a single RBC carries

regions among domains that optimizes the share between

a hydrodynamic shape with linear size of 4 μ m and 8 μ m for thecomputations and communications.

smallest and largest principal directions of the globule, respectively,
and corresponding to a hematocrit level of 58% .

Fig.4 shows the distribution of the time spent incommunication for a run employing 4000 GPUs. On average, most

Most computations are performed in single floatingof the tasks (∼3000 out of 4000) spend about 50% of the total time

point precision with only few reduction operations carried out inin communication. This result confirms that the asynchronous

double precision. Each GPU thread is responsible for the update, in

communication scheme in use works at sustained rate, and thethe LB phase, of a number of mesh nodes that depends on the total

CPU actually plays the role of MPI co-processor of the GPU, where

number of available GPU (the thread configuration is fixed on eachmost of the computations take place.

GPU). For instance, with 512 GPU, each thread is in charge of 8 mesh

For 1334 nodes of the TSUBAME2 (4000 MPI tasks), we

nodes. For the MD phase, interactions among the particles are

estimate a (weighted) average performance of slightly less thanprocessed on a per-particle basis. In this case, the grid of threads is

600 TeraFlops (see Fig. 5 for the corresponding breakdown). Just

directly mapped onto the arrays of particles. Threads are assigned

to convey the flavor of the practical impact of this application,

to particles according to a global id and the search for interacting

the above performance corresponds to simulating a complete

pairs proceeds for each thread in an independent fashion. Eachheartbeat at microsecond resolution and fully inclusive of red

thread scans the cell neighbors and, for each interacting pair, it

blood cells, in 48 hours time on the full TSUBAME2 system.

computes the contribution to the total force. The obtained results
provide a fundamental check of the reliability of the code up tophysiological levels of hematocrit.
Fig.3 (Upp er p anel) shows the elapse d time p er
simulation-step, as well as the breakdown for the LB and MD
components separately. A few comments are in order. At first, theperformance of the Lattice Boltzmann component of MUPHY on
a single GPU is in line with other, highly tuned, CUDA LB kernels.
Secondly, the elapsed time decreases significantly with the numberof cores, with a speed-up of 12.5 between the 256 and 4000 GPU
configurations, corresponding to a parallel efficiency around 80% .
The calculation of the MD direct forces features an efficiency ashigh as 95% . The LB component performs in an optimal way, witha work share that stays always below 4% .
The results show that the combination of asynchronous
communication and overlap bet ween communication andcomputation, significantly alleviates the lack of support for the
direct exchange of data among GPUs.
Fig.3 (Lower panel) shows the total parallel efficiency
referred to either 256 GPUs, if the available memor y allows

Figure 3 Upper Panel: Elapsed time per timestep.
Lower Panel: Parallel efficiency.

simulating the hematocrit level, or 512 GPUs for higher levels

14ACM Gordon Bell Prize: Honorable Mention

Large scale biofluidics simulations
on TSUBAME2

Summary5

In the context of computational biofluidics, we have performed
the first large-scale simulation of the entire heart-sized coronary
system on a world-class cluster of GPUs. The hemodynamic
system consists of a realistic representation of the complex humanarterial geometry at the spatial resolution of red-blood cells. The
computational environment involves one-billion fluid nodes,
embedded in a bounding space of one trillion voxels and coupledwith the concurrent motion of hundreds of millions of red-blood
cells. We achieved a close to 600 Teraflops performance on the4000 GPU configuration of the TSUBAME2 supercomputer, with a
Figure 4 Distribution of the percentage of total timespent in MPI running with 4000 tasks.

parallel efficiency in excess of 90 percent, performing about 2000
billion lattice updates per second concurrently with the simulationof the dynamics of up to 450 million red blood cells.
The above accomplishment results from thedevelopment of several unique features, in terms of both highperformance technology and physical/computational modeling.
We are not aware of any previous implementation dealing with
non-idealized geometries. The present work represents a major
progress in the predictive capabilities of computer simulations forreal-life biofluidic research, with special, yet not exclusive, focus on
cardiovascular clinical practice.
Acknowledgements
This research was carried out as a TSUBAME Grand Challenge
Program.
We thank E. Kaxiras, C.L. Feldman, A.U. Coskun, F.J. Rybicki, A.G.
Whitmore, G. Amati, F. Pozzati and F. Schifano for numerous
discussions.
References[1]

D.A. Vorp, D.A. Steinman, C.R. Ethier, Comput. Sci. Eng., pp. 51[2]

A. Quarteroni, A. Veneziani, P. Zunino, SIAM J. Num. Analysis,[3]

L. Grinberg, T. Anor, E. Cheever, et al., Phil. Trans. Royal Soc. A,(2001).
Figure 5 Breakdown of computational components
and relative performances.

39, 1488 (2002).
367 1896 2371 (2009).
[4]

M. Bernaschi, S. Melchionna, S. Succi et al., Comp.Phys. Comm.,
180, 1495, (2009).

[5] S. Melchionna, M. Bernaschi, S. Succi et al, Comp.Phys. Comm.,
181, 462, (2010).
[6]

S. M el chio nna, Ma cro m o l. T h e o r y Sim ., D O I: 1 0 .1 0 0 2 /mats.201100012 (2011).

15

SC’11 Special Issue[7] M. Bernaschi, M. Fatica, S. Melchionna, S. Succi and E. Kaxiras,
Concurrency and Computation: Practice and Experience, DOI:
10.1002/cpe.1466 (2009).
[8]

M. Bisson, M.Bernaschi, S.Melchionna, Commun. Comput.
Phys., 10, 1071 (2011).

[9]

S. Melchionna, J. Comput. Phys. 230, 3966 (2011).

[10] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics
and Beyond, Oxford University Press, USA (2001).
[11] A. Peters et al., Proceedings of Supercomputing 2010, New
Orleans, 2010.
[12] http://www.labri.fr/perso/pelegrin/scotch16Graph500 Challenge on TSUBAME 2.0
Toyotaro Suzumura*/** Koji Ueno**Graduate School of Information Science and Engineering, Tokyo Institute of Technology**IBM Research – Tokyo

Graph500 is a new benchmark that ranks supercomputers by executing a large-scale graph search problem.
Our early study reveals that the provided reference implementations by the Graph500 benchmark are not scalable
in a large-scale distributed environment. For the Graph500 benchmark on TSUBAME 2.0, we ran our highly scalable
BFS method that divides adjacent matrix - representing large-scale graph - with 2D partitioning and distributes the
portion to all the processors. In contrast to traditional 1D partitioning, this can greatly reduce the communicationexchange. With the optimization method, we successfully solved BFS (Breadth First Search) for large-scale graph with
236（68.7 billion）vertices and 240（1.1 trillion）edges for 10.955 seconds with 1366 nodes and 16392 CPU cores. This
record corresponds to 100.366 GE/s, which was the 3rd-ranked score in the latest ranking announced in SC2011.

Introduction

1

undirected graph from the graph generator in a format usable byKernel 2. The first kernel transforms the edge tuples (pairs of startand end vertices) to efficient data structures with sparse formats,such as CSR (Compressed Sparse Row) or CSC (Compressed Sparse
Column). Then Kernel 2 does a breadth-first search of the graph

Large-scale graph analysis is a hot topic for various fields of study,

from a randomly chosen source vertex in the graph.

such as social networks, micro-blogs, cyber security, protein-protein

The benchmark uses the elapsed times for both kernels,

interactions, and the connectivity of the Web. The numbers ofbut the rankings for Graph500 are determined by how large thevertices in the analyzed graph networks have grown from billions

problem is and by the throughput in TEPS (Traversed Edges Per

to tens of billions and the edges have grown from tens of billions toSecond). This means that the ranking results basically depend on

hundreds of billions. Since 1994, the best known de facto ranking

the time used by the second kernel.

of the world’s fastest computers is TOP500, which is based on a

After both kernels have finished, there is a validationhigh performance Linpack benchmark for linear equations. As anphase to check if the result is correct. When the amount of dataalternative to Linpack, Graph500 [1] was recently developed. We

is extremely large, it becomes difficult to show that the resulting

con-ducted a thorough study of the algorithms of the reference

breadth-first tree matches the reference result. Therefore theimplementations and their performance in an earlier paper [2]. Basedvalidation phase uses 5 validation rules. For example, the first ruleon that work, we implemented a scalable and high-performanceis that the BFS graph is a tree and does not contain any cycles.

implementation of an optimized Graph500 benchmark for large

There are six problem classes: toy, mini, small, medium,

distributed environments. In this paper, we give an overview of thelarge, and huge. Each problem solves a different size graph definedGraph500 benchmark in Section 2 and its basic parallel algorithm

by a Scale parameter, which is the base 2 logarithm of the number

called level-synchronized BFS in Section 3. Our proposed scalableof vertices. For example, the level Scale 26 for toy means 226 and

BFS method is described in Section 4 and the per formancecorresponds to 1010 bytes occupying 17 GB of memory. The six

evaluation is shown in Section 5, and then give a conclusion in

Scale values are 26, 29, 32, 36, 39, and 42 for the six classes. The

Section 6.

largest problem, huge (Scale 42), needs to handle around 1.1 PBof memory. As of this writing, Scale 38 is the largest that has been

GRAPH500 Benchmark

2

In this section, we give an overview of the Graph500 benchmark [1].
In contrast to the computation-intensive benchmark used by
TOP500, Graph500 is a data-intensive benchmark.
It does breadth-first searches in undirected large graphs
generated by a scalable data generator based on a Kronecker
graph [16]. The benchmark has two kernels: Kernel 1 constructs an

17solved by a top-ranked supercomputer.

SC’11 Special Issue

neighbors (the neighboring vertices that have not been visited
at levels 0 or 1). The algorithm also maintains NQ (Next Queue),containing the vertices that should be visited at the next level.
After visiting all of the nodes at each level, the queues CQ and NQ
are swapped at line 16.
VISITED is a bitmap that represents each vertex with one
bit. Each bit of VISITED is 1 if the corresponding vertex has been
already visited and 0 if not. PRED has a predecessor vertex for eachvertex. If an unvisited vertex v is found at line 12, the vertex u is thepredecessor vertex of the vertex v at line 14. When we complete
BFS, PRED forms a BFS tree, the output of kernel2 in the Graph500
benchmark.
At each level, the set of all vertices v is the NQ nominee,
now called “NQ-N”. NQ-N has all the adjacent vertices of the

Figure 1 Kronecker Graph [4]Level-Synchronized BFS

vertices in CQ that would be potentially stored in NQ. NQ-N is3

All of the MPI reference implementation algorithms of theGraph500 benchmark use a “level-synchronized breadth-firstobtained from line 9 to line 11 in the algorithm.
The Graph500 benchmark provides 4 different reference
implementations based on this level-synchronized BFS method.
Basically all the reference implementations has the scalability
issue and are mostly saturated with 32 nodes on TSUBAME 2.0.
Please refer our work [2] on their details and algorithms and theirperformance analysis on TSUBMAE 2.0.

search”, which means that all of the vertices at a given level ofthe BFS tree will be processed (potentially in parallel) before any
vertices from a lower level in the tree are processed.

U-BFS: Our Scalable BFS Method

4

To solve the scalability issue of reference implementations, ourscalable BFS method named U-BFS is based upon the techniqueproposed in [3]. Their approach is based on the level-synchronized
BFS and 2D (two dimensions) partitioning technique. This isscalable technique to reduce the communication cost unlike the1D partitioning including vertical partitioning and horizontal
partitioning of reference implementations. Our proposed method
also optimizes this 2D partitioning technique and also other some
optimization techniques. Therefore we give a brief overview of the2D partitioning technique here.
Assume that we have P processors in total, P=R*C
processors are logically deployed in two dimensional mesh of
Algorithm I is the abstract pseudocode for the algorithm

R (processor row) * C (processor column). An adjacent matrix isthat implements level-synchronized BFS. Each MPI process has two

divided as shown in Figure 2 and the processor (i, j) is responsible

queues, CQ and NQ, and two arrays, PRED for a predecessor arrayfor handling the C blocks fromand VISITED to track whether or not each vertex has been visited.

R*C blocks and the processor (i, j) handles the k th block where k is

At any given time, CQ (Current Queue) is the set of ver-. Vertices are divided into

computed by (j-1)*R + i.

tices that must be visited at the current level. At level 1, CQ will

Each level of level-synchronized BFS method with 2D

contain the neighbors of r, so at level 2, it will contain their pending

partitioning is performed by 2 phases called “expand” and “fold”.

18Graph500 Ranking No.3

Graph500 Challenge on TSUBAME 2.0

In the expand phase, every processor copies its CQ to all the otherprocessors that exists in the same column just as the vertical 1D

5.1 Evaluation Environment

partitioning. In the fold phase, each row of processors integrates

Each TSUBAME 2.0 node has two Intel Westmere EP 2.93 GHz

the compute results after the expand phase. Therefore, this isprocessors (Xeon X5670, 256-KB L2 cache, 12-MB L3) and 50 GB

equivalent to a method of combining two types of 1D partitioning.

of local memory. As the software environment we used gcc 4.3.4

If C is 1, this corresponds to the vertical 1D partitioning and if R is 1,(OpenMP 2.5), MVAPICH2 version 1.6 and we used at maximum

it corresponds to the horizontal 1D partitioning.

1366 nodes. TSUBAME 2.0 is also characterized as a supercomputer

In the fold phase, it firstly searches all the adjacentwith heterogeneous processors including tremendous amountvertices against each vertex, v of CQ obtained by the expandof GPU devices, but we do not use the environment. One node

phase, and then sends a tuple of (v, u) - where u is one of the foundof TSUBAME 2.0 has physical 12 CPU cores and virtually 24 cores

adjacent vertices – to the corresponding processor where v iswith SMT (Simultaneous Multithreading). Our implementation

located.

treats 24 cores for one single node and the same number of
The advantage of 2D par titioning is to reduce theprocessors are allocated to each MPI processes. Each computing

number of processors that needs communication among them.

node is connected to two QDR Infiniband network links, so the

The two types of 1D partitioning requires all-to-all communication.

communication bandwidth for the node is about 80 times larger

However, the 2D p ar titioning c an re duce the numb er ofthan a fast LAN (1 Gbps). Not only the link speed at the end-point

communication processors and most importantly becomes highly

nodes, but the network topology of the entire system heavily

scalable in large computing environments since the expand phaseaffects the performance such an I/O intensive application as

only requires the communication among the nodes in the same

Graph500. TSUBAME 2.0 uses a full-bisection fat-tree topology,

column and the fold phase only requires the communication

which accommodates applications that need more bandwidthamong the processors in the same row.

than provided by such topologies as a torus or mesh.

5

Performance Evaluation

5.2 Performance Result
We compare U-BFS with the latest version (2.1.4) of reference
implementations. Figure 3 compares our optimized
implementation, U-BFS with reference implementations. This
experiment is conducted in a weak-scaling fashion, so the problem

We conducted the Graph500 benchmark on TSUBAME 2.0. Here issize for one node is SCALE 26. The horizontal axis is the number ofthe evaluation environment and performance result.

nodes and the vertical axis is TEPS (GE/s).
U-BFS and two reference implementations, replicatedcsr and replicated-csc (called R-CSR and R-CSC hereafter) use
2 MPI processes for 1 node. The reference implementation,
simple (called SIM hereafter), uses 16 MPI processes for one node
since the implementation has not been implemented with themultithreading parallelism. As shown in the graph, there exists
some result that cannot be measured due to the errors suchas validation error, segmentation fault, and memory error withreference implementations. Figure 3 shows that U-BFS outperforms
R-CSC and SIM. With smaller number of nodes less than 32
nodes, R-CSR shows better performance, but our method shows
performance advantage with more than 32 nodes. For instance, our
optimized method is 2.8 times faster than R-CSR with 128 nodes
and SCALE 26 for one node (All the problem size is SCALE 33).
The scalability of our approach is shown in Figure 4. As
shown in the figure, the throughput in TEPS is linearly increasing
with larger number of nodes and finally achieves 99.0 GE per

Figure 2 2D Partitioning Based BFS

19

[3]

second (TEPS) with 1024 nodes and SCALE 36. We successfully

SC’11 Special Issue

solved BFS (Breadth First Search) for large-scale graph with 236（68.7 billion）vertices and 240（1.1 trillion）edges for 10.955
seconds with 1366 nodes and 16392 CPU cores. This record

Concluding Remarks

corresponds to 100.366 GE/s, which was the 3rd-ranked score in6

the latest ranking announced in SC2011.
In this paper we described our challenge and our basic algorithm for
the Graph500 benchmark on TSUBAME 2.0. Our proposed approachbased on 2D partitioning greatly outperforms the benchmark
reference implementations and also shows great scalability with up
to 1366 nodes. Overall our score was 3rd score announced in 2011/11.
Our achievement could not be done with only 2D partitioning but
also more optimization techniques. This is the first challenge forGraph500, but this continuous challenge and research towards
more optimization and scalable algorithms waits us.
Acknowledgements
This research was supported in part by TSUBAME Grand Challenge
Program and a CREST project "Highly Productive, High Performance
Figure 3 Performance Comparison with
Reference Implementations

Application Frameworks for Post Petascale Computing" from Japan
Science and Technology Agency (JST).
References[1]

Graph500 : http://www.graph500.org/.

[2]

Toyotaro Suzumura, Koji Ueno, Hitoshi Sato, Katsuki Fujisawa
and Satoshi Matsuoka, "Performance Evaluation of Graph500
on Large-Scale Distributed Environment", IEEE IISWC 2011
( IEEE International Symposium on Workload Characterization) ,2011/11, Austin, TX, US[3]

Andy Yoo, et al, A Scalable Distributed Parallel Breadth-First
Search Algorithm on BlueGene/L. SC 2005.

[4]

J. Leskovec, D. Chakrabarti, J. Kleinberg, and C. Faloutsos, “Realistic,mathematically tractable graph generation and evolution, using kronecker multiplication,” in Conf. on Principles and Prac-

Figure 4 Performance of Our Optimized
Implementation with
Scale 26 per 1 node

tice of Knowledge Discovery in Databases, 2005.

20

Physis: A High-Level Stencil Framework for
Heterogeneous Supercomputers
Naoya Maruyama* Tatsuo Nomura** Kento Sato*** Satoshi Matsuoka**Global Scientific Information and Computing Center, Tokyo Institute of Technology**Google, Inc.

***Graduate School of Information Science and Engineering, Tokyo Institute of Technology

We propose a compiler-based programming framework that automatically translates user-written structured gridcode into scalable parallel implementation code for GPU-equipped clusters such as TSUBAME2.0.
Our framework automatically translates user-written stencil functions to GPU execution code as well as message
passing parallel code for inter-node parallelism. It also includes several optimizations for better scalability with alarge number of GPUs, such as compute and communication overlapping. We present an overview of our framework
and report performance results using TSUBAME2.0, which demonstrate good scalability up to 256 GPUs.

Introduction

1

To s o lve the prob lem and improve pro gr amm er
productivit y, we envision a high-level programming modelthat provides a uniform application programming inter face
for heterogeneous systems. While low-level inter faces areindeed essential when the maximum programming flexibility is

21

Heterogeneous computing with both conventional CPUs andrequired, such a case should not be common but exceptional,

vector-oriented GPU accelerators is becoming common because

and simplifying programming even with limited flexibilities andof superior performance as well as power efficiency. The peaksmall performance cost should be highly important to allow theperformance of latest NVIDIA GPUs can be as high as 515 GFLOPS per

adoption of heterogeneous systems for a wider range of application

chip, which is faster than the latest CPUs by several factors, allowing

programmers.

significant performance boost in compute-bound applications such

This ar ticle presents our high-level programmingas N-body problems. GPU’s memory bandwidth is also much greater

framework called Physis that is specialized to stencil computations

than conventional CPU memory, reaching over 100 GB/s in properlywith re gular multidimensional Car tesian grid. [1] In s tencilaligned memory accesses, making it possible to achieve significant

computations, each grid point is repeatedly updated by only

speedups in memory-bound applications such as computational

using neighbor points, exhibiting regular spatial locality. Such a

fluid dynamics.[2]computation pattern, called “ structured grids ”, frequently appears in

Programming such heterogeneous systems, however, isnumerical simulation codes for solving partial differential equations.

a notoriously difficult task. The reason is two-fold. First, most of the

The performance of stencil applications is often determined byexisting programming models for such systems only provide low-memory system performance since the typical byte-per-flop ratio inlevel platform-specific abstractions. The lack of high-level unified

such code is higher than the ratio of today’s processor and memory

programming models forces the programmer to learn multiple

systems, including GPUs. Therefore, optimizing data movements isdistinctive models for parallel computing, e.g., message passing

the most important to improve performance of such applications.

for distributed memory machines and GPU-centric models for

Typical such optimizations for the GPU include latency hiding by

accelerators, often resulting in ad hoc hybrid programming models.

scheduling a large number of concurrent threads, data alignment

Since parallel programming even with a single model is known toto allow coalesced memory accesses, and locality optimizations by

be difficult and error prone, exploiting the potential performancethread blocking.[3] In addition to these optimizations, more coding

advantage with hybrid models is thus a highly difficult task. Second,

effort is required to scale well with a large number of GPUs, such as

in the current heterogeneous architecture, data movements often

communication and computation overlapping. GPU performance

involve complex performance considerations such as locality

scalability is especially important for applications using a large

optimizations for keeping data close to processor cores andamount of data, since a single GPU is equipped only with a few giga

overlapping of communications and computations. While thesebytes of memory.

techniques have long been well known and studied on parallel

In the Physis framework, we design its programming

platforms, realizing them on complex heterogeneous systems

model such that architecture neutralit y can be realized on

further increases the programmer burden. As a result, further scaling

various parallel plat forms, with a par ticular focus on GPU-

performance with a large number of GPUs remains to be challenging

based heterogeneous supercomputers. It provides portable andeven for highly skilled experts.

declarative constructs for describing stencil computations, such as

SC’11 Special Issue

creating multidimensional grids, data copying to and from them,and applying stencils over them. Global view memory model

High-Level Framework For Stencil

and implicit parallelism are adopted to realize high productivity

Computations

as well as architecture neutrality. The declarative programming

2

interface at the same time allows for static compilation techniques
to automatically parallelize stencil computations over distributed

We design a high-level programming framework that provides

memory environments with optimizations such as compute-a highly produc tive programming environment for stencil

communication overlapping. While it is beyond the scope of thiscomputations. The framework consists of a domain-specificpaper, the framework is designed to allow for further advanced

language and platform-specific runtimes. The DSL allows forsoftware techniques to be applied transparently for the application

declarative and flexible descriptions of stencils in an architecture-

programmer, such as model-based and experimental performanceneutral way, which is then translated to architecture-specific code bytuning, resiliency through error checking and scalable and fastsource-to-source translators. The framework runtime encapsulates

checkpointing.

architecture-specific data management tasks and provides a

We describe an implementation of the framework based

uniform interface of virtual shared memory for multidimensionalon the standard C language. We introduce a small set of custom data

grids. The rest of this section discusses our major design goals of thetypes and intrinsics for stencil computations into C as an embedded

framework.

DSL .

[4]

Those custom extensions are translated to platform native

code, such as CUDA for GPU and MPI for message passing. Programs

2-1 Design Goals

written in the Physis DSL can also be automatically translated toAutomatic parallelization: We design the Physis DSL amenable toparallel code using MPI with the overlapping optimization for better

compiler-based automatic parallelization on distributed-memory

scalability.

parallel environments. Although automatic parallelization has been
To evaluate our framework, we implement several

an active research topic for the past decades, it has not been widely

stencil applications in the Physis DSL and evaluate its performance

successful in practice for general-purpose languages, especiallyusing the TSUBAME2.0 supercomputer at Tokyo Tech, which is theon distributed memory environments, since effectively exploiting

fifth fastest machine at the Nov 2011 list of Top500. We present

data localities available in applications is a complex and difficult

results of performance studies using up to 256 NVIDIA Fermi GPUs,task. In contrast, our framework is limited to a small set of domain-and demonstrate that our framework can achieve performance

specific computations, but by doing so we eliminate the difficulties

comparable to hand-written versions with good strong and weakof automatic parallelization in conventional general-purpose

scalability. For more complete presentation, refer to our SC11 paper .[1]languages, and realize implicit parallelism on a variety of parallel
platforms.

Embedded DSL rather than external DSL: Inventing a completely
new language for a given problem domain, i.e., an external DSL
approach, potentially allows for a maximally optimized language
design. In practice, however, being dissimilar to existing familiar
languages may hinder adoption by a wide body of application
programmers. We design our DSL as a small set of extensions
on existing general-purpose languages, i.e., an embedded DSL.
We choose C as the base language in our current design andimplementation since it is one of the most commonly used languages
in high performance computing.

Declarative and expressive programming model: In order toimprove productivity, we maximize programming abstraction by
adopting a declarative programming model that allows for lessmanual programming than imperative models. For example, in
the Physis DSL the programmer expresses how each grid element
is computed, but it is determined by the framework how thewhole grid is processed with the user-specified computation; the22

SC’11 Technical Paper

Physis: A High-Level Stencil Framework for
Heterogeneous Supercomputers

overall computation may be performed sequentially or in parallel

The set of size_t parameters specify the indices of adepending on target environments of the framework. Having too

point within the given grid, so the number of index parametersmuch abstraction, however, can be too restricted to implement

depend on the dimensionality of the grid (e.g., three for 3-D grids).

real-world scientific simulations. For example, some stencils may

The return type of PSGridGet and the v parameter of PSGridSet andbe applied only to a part of a whole grid, such as boundary regions.

PSGridEmit have the same type as the element type of the grid,

Although we attempt to keep the language extensions minimal,

which is either float or double.

we adopt additional domain-specific constructs if they can further

PSGridGet returns the value of the specified point, while

improve productivity of programmers and performance of final

PSGridSet writes a new value to the specified point. PSGridEmit

implementation codes.

performs similarly to PSGridSet, but does not accept the index
parameters, and is solely used in stencil functions.

Programming Model

3

3-1 Example
Fig. 1 shows a function of 9-point stencil on 2-D grids. Such afunction can be applied to grids by using two declarative intrinsics:
PSStencilMap and PSStencilRun. Fig. 2 illustrates how these intrinsics

The Physis DSL extends the standard C with several new data types

can be used to invoke the diffusion stencil of Fig. 1 on 2-D grids.

and intrinsics for stencil computations. The user is required to use
the extensions to express stencil-based applications, which are thentranslated to actual implementation code by the Physis translator.
Physis suppor ts multidimensional Car tesian gridsof floating-point values (either float or double). To represent
multidimensional grids, we introduce several new data types named
based on its dimensionality and element type, e.g., PSGrid3DFloat
for 3-D grids of float values and PSGrid2DDouble for 2-D grids ofdouble values. The type does not expose its internal structure, butrather works as an opaque handle to actual implementation, which
may differ depending on translation targets.
Since many of the Physis intrinsics are overloaded withrespect to the grid types, below we simply use PSGrid to specify
different grid types when not ambiguous.
Grids of type PSGridFloat3D can be created and destructedwith intrinsics PSGridFloat3DNew and PSGridFree, as defined as follows:

Figure 1 Example 9-point stencil
PSStencilMap creates an object of PSStencil, which
encapsulates a given stencil function with its parameters bound toactual arguments. It is analogous to closures in functional
programming. PSStencilRun executes PSStencil objects
in a batch v. Each stencil function may be executed in parallel,
exhibiting implicit parallelism.

PSGrid3DFloat PSGrid3DFloatNew(
size_t dimx, size_t dimy,size_t dimz,
enum PS_GRID_ATTRIBUTE attr)
void PSGridFree(PSGrid g)
They can be accessed both in bulk and point-wise ways
using the following intrinsics:
void PSGridCopyin(PSGrid g,const void *src)
void PSGridCopyout(PSGrid g, void *dst)
PSGridGet(PSGrid g, size_t i,size_t j, size_t k)void PSGridSet(PSGrid g,
size_t i, size_t j, size_t k, T v)
void PSGridEmit(PSGrid g, T v)

23

Figure 2 Example code to apply stencils to grids

SC’11 Special Issue

Experimental Evaluation4

performance of seismic benchmark exhibits significant drop at 64
GPUs and relatively low scalability afterward, which remains to be a
subject of more detailed performance analysis.

To evaluate our proposed framework, we have implemented
a protot ype for parallel machines using the ROSE compiler
framework . [5] It consists of a source-to-source translator andruntime components for each target platform. As translation
targets, we currently generate C for CPU execution and CUDA
for GPU execution. In addition, for platforms involving multiple
distributed compute resources, we generate message-passing
parallel code using MPI.
As performance benchmark programs, we implemented
two stencil codes in Physis:
・Diffusion: 3-D 7-point stencil.

Figure 3 Weak scaling performance of Diffusion
benchmark with up to 256 GPUs.

・Seismic: 3-D seismic wave simulation with 27 stencil functions.
Diffusion is relatively small scale, consisting of only one
stencil function, whereas the seismic code consists of 27 stencils withstaggered grids. Among the 27 stencils, six are for computing the 2-D
surfaces of the 3-D domain, which are implemented with PSDomain
objects. All benchmarks use single-precision floating-point data.
We use the TSUBAME2.0 supercomputer at Tok yo
Institute of Technology, which consists of 1408 compute nodes.
Each node has two Intel Xeon Westmere-EP 2.9GHz CPUs andthree NVIDIA M2050 GPUs with 52GB and 3GB of system and GPU
memory, running SUSE Linux Enterprise Server 11 SP1. The computenodes are interconnected by dual QDR Infiniband networks with afull bisection-bandwidth fat-tree topology network. We use CUDA
v3.2 for GPU code and gcc/g++ v4.1.2 for CPU code.

4-2 Strong Scaling Evaluation
Fig. 5 shows the strong scaling performance of the diffusion stencil
with the problem size fixed at 512 x 512 x 4096. We evaluated 1-D,
2-D, and 3-D decompositions using up to 128 GPUs. In the 1-D
decomposition, we uniformly decomposed the z-direction by the
number of GPUs. In 2-D, we also decomposed the y-direction by
two GPUs and again uniformly decomposed the z-direction with therest of GPUs. Similarly, in the 3-D decomposition, we decomposed
the x-direction into two GPUs in addition to the uniform y- andz- direction decompositions. As expected, the 1-D and 2-D cases
performed better with a smaller number of GPUs, but as the number
increases, the 3-D version outperformed the other two versions.
While it is well known that 3-D decomposition often performs

4-1 Weak Scaling Evaluation

better in large-scale settings, our contribution is to allow application

Fig. 3 shows the results of weak scaling evaluation with thescalability without investing significant amount of efforts.

diffusion code. The red and blue lines are the cases where eachscientists to transparently enjoy such better performance andGPU is assigned a subdomain of 256x128x128 and 512x256x256,
respectively. In both cases, the 3-D domain is decomposed only
over y and z dimensions. As expected, the larger problem size
allowed for better performance and scaling, which is almost linear
scaling up to 256 GPUs, although even the smaller case achieved 11
times speedup with 256 GPUs compared to the 16-GPU case.
Fig. 4 shows the results of weak scaling evaluations with theseismic code, where each GPU computes a subdomain of 256 3 region.
Unlike the diffusion case, the problem domain is decomposed over x
and y dimensions; in other words, the domain is expanded in the x-yplane with the problem size of each GPU fixed. The decomposition
implies that boundary exchanges involve non-unit stride data
accesses, thus resulting lower scalability than the diffusion code. The

Figure 4 Weak scaling performance of Seismic benchmark
with up to 144 GPUs. The problem size of each GPU
is fixed at 256x256x256

24

SC’11 Technical Paper

Physis: A High-Level Stencil Framework for
Heterogeneous Supercomputers

Related Work5further optimized code and to evaluate its effectiveness using awider variety of stencil applications.
Acknowledgements
This project was par tially suppor ted by JST, CREST throughRecent developments of GPU accelerators for scientific computing haveits research program: ``Highly Productive, High Performanceenabled low-cost power-efficient approaches to increase compute

Application Frameworks for Post Petascale Computing.'' We also

performances. However, one of the side effects of this trend is thethank MEXT Grant-in-Aid for Young Scientists (22700047), JST-ANR

significant decrease of programmer productivity due to the complexities

FP3C, and NVIDIA CUDA Center of Excellence.

involved in programming heterogeneous architecture. Although little
work has been done for large-scale heterogeneous supercomputers,
several recent projects attempted to solve the problem.

References[1]

N . M a r u y a m a , T. N o m u r a , K . S a t o , a n d S . M a t s u o k a ,Mint is a high-level directive-based framework for stencil“ Phy sis: an imp licitl y p ar all el p ro gr ammin g m o d el fo r

computations .[6] It allows for regular loop-based stencil programss ten cil co mp u t atio ns o n larg e - s c al e GPU -a ccel er ate dto be annotated with its custom directives so that stencil loopss u p e r c o m p u t e r s , ” I n P r o c e e d i n g s o f 2 0 11 AC M / I E E E

can be executed on GPUs. Ypnos is a Haskell-based DSL for stencil

International Conference for High Performance Computing,
Networking, Storage and Analysis (SC'11), Seattle, WA, USA (2011)

computations that is designed so that compiler-based automaticparallelization is possible .[7] Both of them share common objectives[2]

T.Shimokawabe, T.Aoki, C.Muroi, J.Ishida, K.Kawano, T.Endo,with ours, such as automatic parallelization, but so far they are limited

A.Nukada, N.Maruyama, S.Matsuoka, “An80-fold speedup, 15.0

to single-GPU platforms, whereas our primary focus is to realize scalable

TFlops full GPU acceleration of non-hydrostatic weather model

multi-GPU implementations.

ASUCA production code” in Proceedings of the 2010 ACM/IEEE

Listz is a DSL for unstructured mesh-based simulations .[8]International Conference for High Performance Computing,
Networking, Storage and Analysis, (SC’10), New Orleans, LA, USA (2010)

As in our framework, actual implementations of the mesh interface arehidden from the programmer, which allows Listz to perform aggressive[3]

S. Ryoo, C. I. Rodrigues, S. S. Baghsorkhi, S. S. Stone, D. B. Kirk, and W. M.

domain-specific optimizations. The implementation of the Listz is based

W. Hwu, “Optimization principles and application performanceon Scala’s extensive language features, which facilitate developments

evaluation of a multi- threaded GPU using CUDA,” In Proceedingsof DSLs. Listz could be use to implement structured-grid stencilof the 13th ACM SIGPLAN Symposium on Principles and practiceof parallel programming, PPoPP ’ 0 8 , zpages 73 – 82, (20 0 8)

applications; however, since it targets unstructured meshes, it may not
be able to fully exploit the optimization opportunities of structured data.

[4]

M. Fowler. Domain-Specific Languages. Addison- Wesley
Professional (2010)

Conclusion

6[5] M. Schordan and D. Quinlan, “A source-to-source architecture
for user-defined optimizations,” In Modular Programming
Languages, volume 2 7 8 9 of Lec ture Notes in Computer
Science, pages 214–223. Springer Berlin / Heidelberg (2003)[6]

In order to improve programmer productivity on large-scale

per formance in 3 D stencil methods with annotatedheterogeneous GPU clusters, we have designed and implemented

c ,” In Pro ce e din gs o f th e Inter natio nal Co nferen ce o n

the Physis framework that supports portable programming of stencil

Supercomputing (ICS’11), ICS ’11, pages 214–224 (2011)

computations with structured grids. The C-based DSL represents a[7] D. A. Orchard, M. Bolingbroke, and A. Mycroft. “Ypnos: declarative,

high-level declarative programming model for stencil computations.

parallel structured grid programming,” In DAMP ’10: Proceedings

The DSL translator and runtime together realize an ef ficientof the 5th ACM SIGPLAN workshop on Declarative aspects ofmulticore programming, pages 15–24 (2010)

implementation of the programming model with optimizations suchas automatic overlapping of computations and communications.

25

D. Unat , X . C a i, a n d S . B . B a d e n . “M int : rea li z in g CUDA

[8]

[H. Chafi, Z. DeVito, A. Moors, T. Rompf, A. K. Su- jeeth, P. Hanrahan,This paper presented our current framework implementation and

M. Odersky, and K. Olukotun, “Language virtualization forevaluations of its productivity and performance. We have shown thatheterogeneous parallel computing,” In Proceedings of the ACM

our framework successfully generates scalable code for up to 256

international conference on Object oriented programming systems

GPUs. We plan to extend the presented framework for generating

languages and applications, OOPSLA ’10, pages 835–847 (2010)

FTI: high performance Fault Tolerance
Interface for hybrid systems
Leonardo Bautista-Gomez* Dimitri Komatitch** Naoya Maruyama*
Franck Cappello † Satoshi Matsuoka ‡ Takeshi Nakamura**** Tokyo Institute of Technology, INRIA*** JAMSTEC

Seiji Tsuboi***** Observatoire Midi-Pyrénées, University of Toulouse† INRIA, University of Illinois‡ Tokyo Institute of Technology

Large scientific applications deployed on current petascale systems expend a significant amount of their execution
time dumping checkpoint files to remote storage. New fault tolerant techniques will be critical to efficiently exploit
post-petascale systems. In this work, we propose a low-overhead high-frequency multi-level checkpoint techniquein which we integrate a highly-reliable topology-aware Reed-Solomon encoding in a three-level checkpoint scheme.
We efficiently hide the encoding time using one Fault-Tolerance dedicated thread per node. We implement our
technique in the Fault Tolerance Interface FTI. We evaluate the correctness of our performance model and conduct
a study of the reliability of our library. To demonstrate the performance of FTI, we present a case study of the Mw9.0
Tohoku Japan earthquake simulation with SPECFEM3D on TSUBAME2.0. We demonstrate a checkpoint overhead as
low as 8% on sustained 0.1 petaflops runs (1152 GPUs) while checkpointing at high frequency.

1

Introduction

hybrid checkpoint using new technologies such as phase change
memories (PCM) [2]. Moreover, some hard failures can be toleratedusing solid-state-drives (SSD) [8] and cross-node redundancy
schemes, such as checkpoint replication or XOR encoding [9] whichallows to leverage multi-level checkpointing, as proposed by

In high performance computing (HPC), systems are built from

Moody et al. [1]. Furthermore, Cheng et al. demonstrated that more

highly reliable components. However, the overall failure rate ofcomplex erasure codes such as Reed-Solomon (RS) encoding cansupercomputers increases with component count. Nowadays,

be used to further increase the percentage of hard failures tolerated

petascale machines have a mean time between failures (MTBF)

without stressing the PFS [3]. Our work goes in the same direction as[14]

measured in hours or days

and fault tolerance (FT) is a well-known

issue. Long running large applications rely on FT techniques tosuccessfully finish their long executions. Checkpoint/Restart (CR)these three studies and partially leverages some of those results.
1-1 Contributions

is a popular technique in which the applications save their state in

We propose a model of a highly reliable erasure code scheme

stable storage, frequently a parallel file system (PFS); upon a failure,based on our topology-aware RS encoding published in previous

the application restarts from the last saved checkpoint. CR is awork [4]. We extend our previous research by studying not only the

relatively inexpensive technique in comparison with the process-scalability of the encoding algorithm but also the impact of thereplication scheme that imposes over 100 % of overhead.

checkpoint size per node and the group size, on encoding andHowever, when a large application is checkpointed, tens

decoding performance. We evaluate and prove the accuracy ofof thousands of processes will each write several GBs of data andour performance model and show that our topology-aware RS

the total checkpoint size will be in the order of several tens of TBs.

encoding scheme is several orders of magnitude more reliable

Since the I/O bandwidth of supercomputers does not increase atthan XOR encoding. We apply our FT-dedicated thread scheme

the same speed as computational capabilities, large checkpoints

presented in previous work [5] in order to further decrease thecan lead to an I/O bottleneck, which causes up to 25% of overhead

checkpoint overhead and we integrate it for the first time in a multi-

in current petascale systems[12]. Post-petascale systems will havelevel checkpoint technique that we implement in our FTI library. Oura significantly larger number of components and an importantevaluation shows that by using FT-dedicated threads in the nodes

amount of memory. This will have an impact on the system’

we can efficiently hide the encoding time, making its overheads reliability. With a shorter MTBF, those systems may require anegligible in comparison with a simple local write checkpoint. We

higher checkpoint frequency and at the same time they will haveextend our evaluation with a functional test in a real case study bysignificantly larger amounts of data to save.

simulating the March 11th

Mw9.0 Tohoku, Japan earthquake andAlthough the overall failure rate of future post-petascale

we present synthetic seismograms for the Hirono seismic recoding

systems is a common fac tor to study when designing F T-

station in Fukushima prefecture. We perform a large scalability andtechniques, another important point to take into account is theoverhead evaluation of our library with SPECFEM3D and we show

pattern of the failures. Indeed, when moving from 90nm to 16nm

that FTI can successfully scale to more than 1000 GPUs and reachtechnology, the soft error rate (SER) is likely to increase significantly,

over 100 TFlops while checkpointing at high frequency and causing

as shown in a recent study from Intel [13, 6]. A recent study by Dong

only about 8 % overhead in comparison with a not checkpointed

et al. explains how this provides an opportunity for local/global

execution.

26

SC’11 Technical Paper (Achieving a Perfect Score)

FTI: high performance Fault Tolerance Interface forhybrid systems

Japan earthquake in the same manner as Nakamura et al[22]. We

2

EVALUATION

use broadband seismograms of IRIS GSN and IFREE OHP seismic
stations with epicentral distance between 30 and 100 degrees. The
broadband original data are integrated into ground displacement
and band-pass filtered in the frequency band 0.002-1 Hz. We use

In this section we evaluate the performance, scalability andthe velocity structure model of the earth IASP91 [18] to calculate theefficiency of our FTI library using SPECFEM3D. All our experiments

wavefield near the source and stations. We assume that the strike

were done on TSUBAME2.0 with the configuration given in table 1.

of the fault plane is 201 degree and the dip angle is 9 degree, basedth

2-1 Simulating the March 11

Mw9.0 Tohoku Japan earthquake

source. The length of a subfault is 20 km along strike. The assumed

In an effort to extend our evaluation with a functional test offault length is 440 km in total, consistent with the aftershock

FTI in a real case simulation with a production level application,

distribution.

we decided to simulate the devastating Mw9.0 Tohoku JapanThe nonnegative least-squares method [21] is employed forearthquake that struck the northeast part of the island on March

constraining the rake angle in the waveform inversion. The resultsth

11 , 2011. The simulation is done with SPECFEM3D on TSUBAME2.0

of the inversion show the bilateral rupture to the northeast and theusing an input model that describes the source fault.

southwest with two main asperities along the fault; maximum slip

SPECFEM3D is used by more than 300 research groups inis of around 40 m with the reverse fault mechanism approximately

the world for a large number of applications, for example to model

100 km northeast of the epicenter and another large slip withthe propagation of seismic waves resulting from earthquakes,

reverse fault mechanism at 100 km southeast of the epicenter. The

seismic acquisition experiments carried out in the oil industry,

total amount of released seismic moment corresponds to moment

or laborator y experiments with ultrasounds in cr ystals. Thismagnitude Mw = 9.1. We calculate synthetic seismograms withapplication won the Gordon Bell SuperComputing award for Best

this source propagation model for a realistic 3D Earth model using

Performance [11] for a calculation of seismograms in the whole 3D

the spectral-element method [20, 23]. As we can observe in figure 1,

Earth down to periods of approximately 5 seconds, carried out

the synthetic seismogram for the Hirono seismic recording station

at 5 teraflops (sustained) on 1944 processors using 14.6 billion

located in the Fukushima prefecture, the East component showsdegrees of freedom stored in 2.5 terabytes of memory on the Earthabout 2m static displacement to the East, which seems to be

Simulator, the fastest computer in the world at that time (2002).

consistent with the observed crustal deformation caused by thisFor the source model, we apply waveform inversion[19]

to obtain slip distribution in the source fault at the 2011 Tohoku,

Table 1 : TSUBME 2.0 Architecture

27

on the Global Centroid Moment Tensor model of the earthquake

earthquake. Figures 1 also shows the other two components at thesame station.

SC’11 Special Issue

finite difference or finite element code; this is intrinsically related
to the fact that in such numerical methods few operations are
performed per grid point, and thus the cost comes mostly frommemory accesses [15, 16].
First, we start with a strong scalability test in which
we evaluate the performance of SPECFEM3D in three cases: nocheckpointing, checkpointing with FTI (L2) and checkpointing
on Lustre [7]. Since the problem size is fixed, the memory used(and therefore the checkpoint size) per GPU decreases when thenumber of GPUs increases. In this experiment, for the FTI tests all
the checkpoints were done with the L2 of FTI, thus we do not takeadvantage of the multi-level scheme of FTI at this point. Sincecheckpoint size decreases, we also decrease the checkpoint intervalin order to decrease the recovery cost in case of failure. All thecheckpoints are done at the application level and we checkpoint
only the strictly necessary data in order to restart the execution; thiscorresponds to about 20% of the memory used by the application.

Figure 1 Synthetic seismograms for the Hirono station
2-2 FTI scalability study with SPECFEM3D
In order to demonstrate the efficiency and scalability of FTI we
decided to evaluate it at large scale with SPECFEM3D. Recently,
SPECFEM3D was ported to GPU clusters using CUDA [15, 16], so it canbe used in hybrid systems such as TSUBAME2.0. It is important to
notice that such seismic simulations do not need double precisionand perform their runs in single precision [15]. Also, we want tohighlight that SPECFEM3D is a memory-bound application, as any

Figure 2 FTI evaluation on TSUBAME2.028

SC’11 Technical Paper (Achieving a Perfect Score)

FTI: high performance Fault Tolerance Interface forhybrid systems

As we can see in figure 7a, SPECFEM3D strong-scales
well from 5 TFlops on 48GPUs to almost 23 TFlops on 384GPUs

as TSUBAME2.0 and checkpointing with our library FTI every 6
minutes (high frequency checkpointing).

without checkpointing. FTI follows closely this progression bycausing only about 4 % of overhead for 384GPUs. In contrast,
checkpointing to Lustre becomes prohibitively costly at highfrequencies. Performance was measured using PAPI to measure

CONCLUSIONS

the floating point operations [17] and each dot in the figure is the3average of 5 runs.
To evaluate the overhead of FTI at large scale, we stressed

In this work we have proposed a highly reliable technique basedeven more our library by weak-scaling to more than 1000 GPUs.

on a topology-aware RS encoding. Also, we have exploited someIn this second experiment, we populate the GPUs memory withcharacteristics of GPU computing through which many GPU

2.1GBs of data, out of 2.6GBs available for the user (12.5 % is usedapplications are capable of spawning one extra FT-dedicatedfor ECC in Fermi GPUs ) and we keep the checkpoint interval fixed

thread per node in order to improve checkpoint performance.

to 6 minutes, which is the Young’s optimal checkpoint interval for

We have integrated both techniques in a multi-level checkpoint

a MTBF of 12 hours and a L1 checkpoint of 2 seconds. Then, we

model that we have implemented in our FTI library and we haverun SPECFEM3D for several configurations: The first one is without

conducted an exhaustive study of correctness of our performancecheckpoint (No ckpt.); the second one is checkpointing to the local

model and the reliability of our library.

SSDs without any encoding (L1); the third one is using FTI, thus in

Moreover, we have conducted for the first time a large

addition to the local checkpoint, every 2 checkpoints FTI will use

scale evaluation of such a multi-level technique with a production

the RS encoding proposed in our model (FTI-L1,L2); the fourth one

level scientific application on an hybrid platform. Our evaluation

is similar to the previous one but in addition every 6 checkpoints

with SPECFEM3D on TSUBAME2.0 shows that FTI imposes only 8%

the latest checkpoint files are flushed to Lustre (FTI-L1,L2,L3); andof checkpoint overhead while running at over 0.1 petaflops andfinally checkpointing with BLCR on Lustre (BLCR+Lustre). Although

checkpointing every 6 minutes.

there are some ongoing works [10] to make it possible, BLCR cannot
currently checkpoint GPU-accelerated systems. Hence, we emulate
it by writing 2.1GBs of data per process (therefore per GPU) to Lustre

This work was supported in part by the JSPS, the ANR/JST FP3C

in the same way BLCR would do it. It is important to highlight thatproject, and by the INRIA-Illinois Joint Laboratory for Petascale

BLCR, as any other kernel-level checkpoint, will save the complete

Computing.

memory of every process, creating a 5 times larger checkpoint.
In figure 7b, we can see that SPECFEM3D has an almost

References

perfect weak scaling, from 43TFlops to 117TFlops on 1152GPUs[ 1] A . Mo o dy, G . Bronevet sk y, K . Mohror, B. R . de Supinsk i,for the No ckpt. test. Also, in the figure the L1 results are actually

Design, Modeling, and Evaluation of a Scalable Multi-level

hidden by the FTI-L1,L2 results. Indeed, both scenarios achieve

Checkpointing System. In ACM/IEEE International Conference

almost identical results causing about 8 % overhead in comparison

for High Performance Computing, Networking, Storage andwith the No ckpt. case. This means, that the RS encoding done at

Analysis, New Orleans, 2010

L2 is completely hidden thanks to the FT-dedicated threads. The[2] X. Dong, N. Muralimanohar, N. Jouppi, R. Kaufmann, Y. Xie.

L1 checkpoints, capable of tolerating transient failures, are done

Leveraging 3D PCRAM Technologies to Reduce Checkpoint

between two L2 checkpoints while the FT-threads are still encodingO v e r h e a d f o r Fu t u r e E x a s c a l e Sy s t e m s . I n AC M / I E E E

the previous, more reliable, checkpoint. The FTI-L1,L2,L3 scheme

International Conference for High Performance Computing,adds an extra 3 % overhead due to Lustre writing performance.

Networking, Storage and Analysis, Portland, 2009.

Finally, the BLCR+Lustre scheme imposes an always larger and[ 3 ] Z. Cheng, J. Dongarra, A scalable Checkpoint Encoding

prohibitive overhead as the size of the problem increases. For each

Algorithm for Diskless Checkpointing. Proceedings of the 11th

run we let the application run between 30 and 40 minutes andIEEE High Assurance Systems Engineering Symposium, HASE

every point in the figure is the average of 3 runs.

29

Acknowledgements

2008, Nanjing, China, December, 2008.

At this p o int , we have achieve d over 10 0T Fl o ps[4] L. Bautista-Gomez, N. Maruyama, A. Nukada, F. Cappello, S.

of sustained performance with a production-level scientific

Matsuoka, “Low-overhead diskless checkpoint for hybrid

application such as SPECFEM3D, on an hybrid supercomputer suchcomputing systems”, International Conference on HighSC’11 Special Issue

Performance Computing, Goa, India, December 2010.
[5] L. Bautista-Gomez, N. Maruyama, F. Cappello, S. Matsuoka,“Distributed Diskless Checkpoint for large scale systems”, IEEE/
ACM International Symposium on Cluster, Cloud and Grid
computing (CCGrid2010), Melbourne, Australia, May 2010.

[20 ] D. Komatitsch, J. Ritsema, J. Tromp, The spectral-element
method, Beowulf computing, and global seismology, Science
298, 1737-1742, 2002.
[ 21] C . Lawson , R . Hanson, Solving Least Squares Problems,
Prentice-Hall, New Jersey, 340 pp, 1974.

[ 6 ] B. Schroeder, E. Pinheiro, W. Web er. DR AM errors in the[22 ] T. Nakamura, S. Tsuboi, Y. Kaneda, Y. Yamanaka, Rupture

wild: A Large-Scale Field Study. In Proceedings of the 11th

process of the 2008 Wenchuan, China earthquake inferredinternational joint conference on Measurement and modelingfrom teleseismic waveform inversion and forward modelingof computer systems (SIGMETRICS), ACM, New York, NY, USA,of broadband seismic waves, Tectonophysics, vol. 491, 72-84,

2009.
[7]

2010.

S. Microsystems. Lustre file system, October 2008[8] A. Moody, G.

[23] S. Tsuboi, D. Komatitsch, C. Ji, J. Tromp, Broadband modelling

Bronevetsky, Scalable I/O Systems via Node-Local Storage:of the 2002 Denali fault earthquake on the Earth Simulator,

Approaching 1 TB/sec File I/O. DOE technical report, 2009

Phys. Earth Planet. Inter. 139, 305-312, 2003.

[ 9 ] W. D. Gropp, R. Ross, and N. Miller. Providing ef ficient I/O
redundancy in MPI environments. Lecture Notes in Computer
Science, 3241:7786, September 2004.
[10] A. Nukada, S. Matsuoka, NVCR : A Transparent CheckpointRes tar t Libr ar y for NV IDIA CUDA in Pro ce e dings at theInternational Heterogeneity in Computing Workshop, Alaska,
2011. (To appear)[11] D. Komatitsch, S. Tsuboi, C. Ji and J. Tromp, A 14 . 6 billion
degrees of freedom, 5 teraflops, 2 . 5 teraby te earthquake
simulation on the Earth Simulator, Proceedings of the ACM /IEEE Supercomputing SC’2003 conference, November 2003.
[12 ] G. Grider, J. Loncaric, and D. Limpart, Roadrunner System
Management Report, Los Alamos National Laboratory, Tech.
Rep. LA-UR-07-7405, 2007.
[13] S. Y. Borkar, Designing Reliable Systems from Unreliable
Components: The Challenges of Transistor Variability and
Degradation, IEEE Micro, vol. 25, no. 6, pp. 10-16, 2005.
[14 ] D. Reed, High-End Computing : The Challenge of Scale ,Director’s Colloquium, LANL, May 2004.
[15] D. Komatitsch, D. Michéa, G. Erlebacher, Porting a high-orderfinite-element earthquake modeling application to NVIDIA
graphics cards using CUDA, Journal of Parallel and Distributed
Computing, vol. 69(5), p. 451-460, doi: 10.1016/j.jpdc.2009.01.006,
2009.
[16] D. Komatitsch, G. Erlebacher, D. Göddeke, D. Michéa, Highorder finite-element seismic wave propagation modeling with

● TSUBAME e-Science Journal No.5

MPI on a large GPU cluster, Journal of Computational Physics,

Published 02/28/2012 by GSIC, Tokyo Institute of Technology ©
ISSN 2185-6028
Design & Layout: Kick and Punch
Editor: TSUBAME e-Science Journal - Editorial room
Takayuki AOKI, Thirapong PIPATPONGSA,
Toshio WATANABE, Atsushi SASAKI, Eri Nakagawa
Address: 2-12-1-E2-1 O-okayama, Meguro-ku, Tokyo 152-8550
Tel:+81-3-5734-2087 Fax: +81-3-5734-3198
E-mail: tsubame_j@sim.gsic.titech.ac.jp
URL:
http://www.gsic.titech.ac.jp/

vol. 229(20), p. 7692-7714, doi: 10.1016/j.jcp.2010.06.024, 2010.
[17] http://icl.cs.utk.edu/papi/[18 ] B. Kennet, E. Engdahl, Traveltimes for global ear thquake
location and phase identification. Geophys. J. Int., 105, 429465, 1991.
[19] M. Kikuchi, H. Kanamori, Note on Teleseismic Body-Wave Inversion
Program, 2003. http://www.eri.u-tokyo.ac.jp/ETAL/KIKUCHI/

30