13Multi-GPU Computation of Multi-phase Field
Simulation of the Evolution ofMetallic Polycrystalline Microstructure
A Large-scale Parallel Computation for VibrationalState
Analysis Based on Quantum Monte Carlo method
Large-scale DEM Simulations for Granular Dynamics

Multi-GPU Computation of Multi-phase
Field Simulation of the Evolution ofMetallic Polycrystalline Microstructure
Akinori Yamanaka* Masashi Okamoto** Takashi Shimokawabe*** Takayuki Aoki**** Division of Advanced Mechanical Systems Engineering, Institute of Engineering, Tokyo University of Agriculture and Technology** Department of Mechanical Systems Engineering, Graduate School of Engineering, Tokyo University of griculture and Technology*** Global Scientific Information and Computing Center, Tokyo Institute of Technology

The multi-phase field method is recognized as one of the most promising simulation tools for predicting the evolutionof polycrystalline microstructures in metals. We have recently developed a massively parallel computation techniquethat uses multiple graphics processing units (GPUs) for large-scale three-dimensional (3D) multi-phase field simulation.
In this article, we introduce this multiple GPU computation technique, including an overlapping method that enables
us to simultaneously perform computations concerning the GPU and data communication. We implemented our
computation technique on the TSUBAME 2 . 5 GPU supercomputer and evaluated its performance. Large-scale 3 D
simulations of polycrystalline grain growth performed on the TSUBAME2.5 exhibited high computing performance.

Introduction

1

overlapping method[3] that enables us to simultaneously carryout computations concerning the GPU and data communication.
By implementing our computation technique to the TSUBAME2.5
supercomputer, we conducted large-scale three-dimensional(3D) MPF simulations of polycrystalline grain growth. We present

In recent years, fuel-efficient vehicles such as hybrid cars havethe results of the performance evaluation of the large-scale MPF

been actively researched and developed to reduce thesimulations thus obtained[4].

environmental load of vehicular transportation. In order to further
improve the fuel efficiency of vehicles, it is necessary not only todevelop high performance engines and motors, but also to reduce
the weight of vehicles by reducing the thickness of the sheet
metal used in them. Meanwhile, sheet metal needs to be

2

Multi-phase field method

strengthened further to ensure the safety of vehicles’ occupantsin case of an accident. However, the strengthening of sheet metal
often causes low formability: when sheet metal is strengthened,

The MPF method used in this study was proposed by Steinbach

cracks and fractures occur in the press forming process. A large

and Pezzola in 1999[1]. The MPF method simulates microstructural

number of trial-and-error experiments have been carried out inevolution in materials by assuming that the total free energy

the past to control the strength and formability of the sheet metal.

of the material monotonically decreases with time. When we

At present, because of intensifying international
competition in the field of new materials, the infrastructure of

consider a system of N crystal grains, the total free energy of thesystem is represented by the following equation:

computational materials engineering has grown increasingly
important. Simulation-based material development is expected
to shorten the development period and help efficiently produce
new materials. In particular, the multi-phase field (MPF) method[1]where the first term of the right-hand side describes the potential

has attracted considerable attention as one of the most promisingenergy and the second term corresponds to the gradient energy.

simulation tools to predict microstructural evolution in metallic

The parameters W i j and a i j are functions of the interfacial energy

materials. In contrast to the phase field simulation of dendritic

and the interfacial thickness, andsolidification using the TSUBAME2.0 supercomputer of the Tokyo

parameter called the phase field variable.

Institute of Technology, which was awarded the ACM Gordon

local volume fraction of the ith crystal grain at coordinate r andBell Prize in 2011[2], an MPF simulation needs to solve multiple

time t . Therefore,

non-linear partial differential equations. Therefore, it requires

otherwise.

more memory and a longer computational time than needed forregion.

conventional phase-field simulations.
In this study, we develop and test a multiple GPU
computation technique for massively parallel computation oflarge-scale MPF simulations. This technique includes an

02(r, t ) is a continuous order(r, t ) describes the(r, t ) admits the value 1 for the i th grain and 0(r, t ) changes smoothly from 0 to 1 in the interfacial

The time evolution equation of the phase field variables(Allen–Cahn equation) can be derived by considering the monotonic
reduction of the total free energy expressed by Equation (1):of large-scale 3D MPF simulations. The details of the APT method
can be found in reference[5]. In this article, we categorize the APT
algorithm into APT1 and APT2. APT1 requires values of the phase
where n is the number of the phase field variables greater than 0,

field variables at each computational grid, whereas APT2 needs

andthe values of variables at neighboring computational grids.

is the mobility of the phase field variable.
In the MPF simulation, Equation (2) is solved by thesecond-order finite difference method for space and the first-order

3.2 Domain decomposition

forward Euler method for time on a regular 3D computational grid.

In order to perform the MPF simulation using multiple GPUs, we

The program code was written in CUDA Fortran.

decomposed an entire computational domain into subdomains.
A set of GPUs and a CPU is allocated for computation to each
subdomain. Figure1 shows the 3D domain composition as an example

Multiple GPU computation of
MPF method

3of domain decomposition. When the entire computational domain
divided by NX × NY × NZ finite difference grids is decomposed into
X , Y , and Z subdomains along X , Y , and z directions, the number ofcomputational grids in each subdomain is NX /X× NY /Y × NZ /Z .
The computation of the time evolution equation (Equation (2)) for
each subdomain requires data on the surfaces of the neighboring

3.1 Active parameter tracking method

subdomains. Therefore, we prepared boundary regions on the

Using the MPF method, we solve Equation (2) with respect not tosurfaces of each subdomain along x , y , and z directions. In the

N phase field variables but only n non-zero phase field variables.

parallel computing, the data in the boundary region is transferred

Therefore, we do not need to save the values of all N phaseby using the Message Passing Interface (MPI) library. In thefield variables. In this study, we employed the active parameter

assessment of the performance of parallel computation,[5]

tracking (APT) method to reduce memory consumption. The APT

we attempted one-, two- and three-dimensional domain

method is an essential algorithm for the efficient implementation

decompositions and selected the most-effective manner of these.

Fig. 1 3D domain decomposition.

03

Multi-GPU Computation of Multi-phase Field Simulation of the Evolution ofMetallic Polycrystalline Microstructure

3. 3 Assignment of threads and blocks

and 4, data in the boundary regions is asynchronously transferredfrom the global memory on the GPU to the host (CPU) memory

It is well-known that the assignment of CUDA threads and blocks

using the cudaMemCpyAsync function in the CUDA applicationis crucial to obtaining high performance from GPU computing.

programming interfaces (APIs).
The CPU transfers data in the boundary regions, calculated

Figure 2 shows the assignment of CUDA threads and blocks toeach subdomain employed in this study. We divided subdomainsby the GPU in the previous time step, to the host memory. Followingof size nx × ny × nz into X ’ × Y ’ small domains along the x and ydata communication between the GPU and the CPU in Streams 2, 3,

directions. Thus, the size of each small domain was nx /X ’ × ny /Y ’

and 4, the CPU initiates the procedures for APT1 and APT2 for the× nz . A thread block of size nx /X ’ × ny /Y ’ × 1 handled each small

boundary regions and carries out data communication among CPUs

domain by marching in the z direction[2]. The optimum number

with the MPI library. Data updates by APT1 and APT2 are transferredof CUDA threads and block were determined by trial and error

back from the CPU to the GPU while Stream 1 is executed on the

because performance depends on the size of the computationalGPU.

domain and the number of GPUs used in the simulation.

Once all streams on the GPU are synchronized withthe computation and data communication on the CPU, the GPU
handles the procedure for APT2 for the boundary regions andupdates the calculated data. The calculated data on the CPU canbe updated while all streams are executed on the GPU.

Fig. 2 Assignment of CUDA threads and blocks.
Fig. 3 The overlapping method.
3. 4 Overlapping method
In order to perform parallel computation using multiple GPUs andCPUs, data communication is needed not only among CPUs, butalso among the GPUs and between GPUs and CPUs. Therefore,
the time required for data communication can degrade parallel
efficiency. Thus, we propose an overlapping method that masks
the time needed for data communication with the time required
to perform computation on the GPUs[4].
Figure 3 shows the computational diagram used for
the MPF simulation using the overlapping method. All operations(kernel execution and data communication) on the GPU are
concurrently run in four CUDA streams (Stream 1~4). In Stream 1,
we compute the time evolution equation for the internal regionof the subdomain and handle the procedure for APT1. The
computation of the time evolution equation for the boundary
regions is carried out in Streams 2, 3, and 4 simultaneously.
Following the completion of the computations in Streams 2, 3,

04

4

Performance evaluation of multiple
GPU computing on TSUBAME2.5

different colors. Polycrystalline grain growth can be observedin the heat treatment of metallic materials, e.g., annealing. The
simulation result successfully reproduced the coarsening andthe shrinking of crystal grains. The statistical evaluation of thepolycrystalline microstructure, e.g., the distribution of crystal
grain size and average grain size, is only possible when a large-

4. 1 Polycrystalline grain growth simulation

scale MPF simulation taking into account a large number of

We evaluated the performance of our proposed multiple GPU
computation by simulating large-scale 3D polycrystalline grain

crystal grains is performed.
In order to evaluate the performance of multiple GPU

growth on the TSUBAME2.5 GPU supercomputer at the Global

computing, we measured variation in the value of Floating-

Scientific Information and Computing Center of the Tokyo

point Operations Per Second (FLOPS) by changing the number

Institute of Technology. All simulations were performed usingof GPUs used for the simulation. Figure 5 shows the results ofsingle-precision floating-point calculation.

the performance evaluation for weak scaling. Here, we evaluatedthe variation in the value of FLOPS with the number of GPUs fora fixed computational grid size per GPU. In this case, each GPU
handled 2563 computational grids and 512 crystal grains. The
results show that the performance improved in proportion to
the number of GPUs, and that the overlapping method improved
performance. As a result, we attained 1.9 TFLOPS using 729 GPUs
for 23043 computational grids and 373,248 crystal grains.

Fig. 4 Polycrystalline grain growth behavior simulatedusing the TSUBAME2.5 supercomputer.
The simulation was performed using 256 GPUs
on 10243 finite difference grids. Crystal grains
were visualized by different colors.

Fig. 5 Weak scaling of multiple GPU computing.

Furthermore, we evaluated the performance of oursystem for strong scaling. Here, the variation in the value ofAs an example of the simulation results, we show

FLOPS with the number of GPUs for a fixed computational

snapshots of the polycrystalline grain growth behavior simulated

domain size was measured. The number of crystal grains was set

using 256 GPUs on TSUBAME2.5 in Figure 4. In this simulation,

to 512 per 2563 computational grids. Figure 6 shows the results

33

the size of the computational domain was 0.512 mm and theof the strong scaling for three different computational domain

number of computational grids was 10243. The initial number ofsizes: 2563, 5123, and 10243. Similar to the results for weak scaling

crystal grains was 32,768. The crystal grains were visualized byshown in Figure 4, the performance improved with increasing

05

Multi-GPU Computation of Multi-phase Field Simulation of the Evolution ofMetallic Polycrystalline Microstructure

number of GPUs. However, in case a small computational domain
size was used, the performance degraded when we used a large

Acknowledgements

number of GPUs. This is because the size of subdomain handledThis work was partially supported by the Joint Usage/Researchby each GPU decreased by increasing the number of GPUs andCenter for Interdisciplinary Large-scale Information Infrastructures

the time for the data communication cannot be masked within Japan. This work was also supported by a Grant-in-Aid forthat for computation. On the other hand, when we used large

Scientific Research (KAKENHI 25630322) from the Japan Society

computational domain sizes and more than 100 GPUs, thefor the Promotion of Science (JSPS).

performance clearly exhibited high scalability.

References[1]

I. Steinbach, F. Pezzola: A Generalized field method for
multiphase transformations using interface fields, Physica
D, Vol. 45, pp. 385-393 (1999)[2]

T. Shimokawabe, T. Aoki, T. Takaki, A. Yamanaka, A. Nukada,
T. Endo, N. Maruyama, S. Matsuoka: Peta-scale phase-field
simulation for dendritic solidification on the TSUBAME2.0
supercomputer, Proceedings of the 2011 ACM/IEEEInternational Conference for High Performance Computing,
Networking, Storage and Analysis, SC’11, IEEE Computer
Society, Seattle, WA, USA, (2011)[3]

M. Okamoto, A. Yamanaka, T. Shimokawabe, T. Aoki:
Multiple GPU Computing of Polycrystalline Grain Growth
Simulation using Multi-phase-Field Method, Transactionof the Japan Society for Computational Engineering and

Fig. 6 Strong scaling of multiple GPU computing.

Science, Vol. 2013, p. 20130018 (2013) (in Japanese)[4]

A. Yamanaka, M. Okamoto, T. Shimokawabe, T. Aoki: Largescale 3D multi-phase field simulation of microstructure
evolution using TSUBAME2.5 GPU supercomputer,

Summary

5

Proceedings of 2nd International Congress on 3D Materials
Science, The Minerals, Metals & Materials Society, pp. 59-64,(2014)[5]

S. G. Kim, D. I. Kim, W. T. Kim, Y. B. Park: Computer simulationof two-dimensional and three-dimensional ideal grain

The multiple GPU computation technique has been developed
for massively parallel computation of large-scale 3D MPF
simulations. In this paper, we successfully applied this technique
was successfully to the TSUBAME2.5 supercomputer. We showedthat extremely large-scale MPF simulations of polycrystalline
grain growth can be efficiently performed using our technique.
The performance evaluation showed that high scalability wasachieved by using the overlapping method.
It remains time-consuming to optimize the size anddistribution of the polycrystalline microstructure that produces
desirable mechanical properties in the material only by trial-anderror experiments. Therefore, we expect that efficient large-scale
MPF simulations using the multiple GPU computing techniqueproposed in this article will contribute to the acceleration ofmaterials development and the reduction of experimental cost.

06

growth, Phys. Rev. E, Vol. 74, p. 061605 (2006)

A Large-scale Parallel Computation for
Vibrational State Analysis Based on
Quantum Monte Carlo method
Ryota Nakayama* Osamu Fujioka** Yukiumi Kita** Masanori Tachikawa*** International College of Arts and Sciences, Yokohama City University** Quantum Chemistry Division, Yokohama City University

We show theoretical outline of vibrational quantum Monte Carlo (vibQMC) method that we have recently developed,
and its parallel implementation method toward a large-scale parallel computation on a super-computer system.
The parallelization ratio of our vibQMC program code is 99.9981% , and the practical parallel efficiency using 5376
cores on TSUBAME 2.5 super-computer system is about 91% .
We also show theoretical results of vibrational state analysis of the monohydrated negative core ion, H 3 O 2 - , which
is a precursor ion in forming aerosols in the atmosphere.

Introduction

1

2

MethodQuantum Monte Carlo (QMC) method is one of the most accurate

In this study, we used two types of QMC technique, variational

ab initio methods to solve the Schrödinger equations of atoms,

Monte Carlo (VMC) and RMC methods, for anharmonic vibrational

molecules, and solids [1]. In particular, diffusion Monte Carlo (DMC)

state analyses of polyatomic molecules. We here briefly describe

method gives a very accurate variational energy of the systems.

theoretical outline of the both methods.

However, the scope of application of the DMC method ispractically limited in the estimation of a variational energy due

2.1. Variational Monte Carlo (VMC) method

to some serious problem known as the population control bias,

We consider the following expectation value

the difficulty to compute the expectation values of the non-operator

commutative operators with Hamiltonian, etc .

operators, respectively) with a given trial wave function(

andof Hamiltonian

are kinetic and potential energy:

Reptation Monte Carlo (RMC) method [2] is one of QMC
methods using the imaginary-time propagation of timedependent Schrödinger equation as well as the DMC method. In
principle, RMC method has the equivalent theoretical accuracy to
the DMC method, but allows the population control bias freealgorithm and the exact estimation of non-commutative operatorswith Hamiltonian. The computational cost in RMC calculations is,wherehowever, extremely expensive compared to that in the DMC

and is referred to as a configuration or walker in QMC methods.

calculations. Thus, the scope of the application of the RMC method

The variable

is a local energy. We assume a real trial

is limited in small molecular systems such as H2O molecule [3].

wave function (). In Eq. (1), the expectation value

In this report, we show the parallel implementation

is a generalized multi-dimension coordinate,

be evaluated as the average ofcanover whole configuration

method of our RMC program developed toward RMC calculations

spaceof large-scale molecular systems. The RMC program has been

simply calculate the value ofdeveloped as a part of vibrational quantum Monte Carlo (vibQMC)

with Metropolis’s method [4], where M is the number of sampling

method that we have recently proposed. In the following sections,points.

we describe the theoretical outline of our vibQMC method, thewith the statistical weight. In VMC method, thus, we

by generating a set of

The VMC method enable us to analyze not only abenchmark calculations of parallel efficiency of vibQMC program

total energy of a system, but also a expectation value of a given

on TSUBAME 2.5 super-computer system, and theoretical results

physical properties, although its theoretical accuracy stronglyof vibrational state analysis of monohydrated negative coredepends on the quality of trial wave function. In QMC calculations-

ion, H3O2 , which is a precursor ion in forming aerosols in thefor many electron system such as atoms, molecules, and solids,

atmosphere, as an application of a large-scale parallel RMC

Slater-Jastrow type trial wave function, which consists of a single

computation.

(or multi) Slater determinant(s) obtained with molecular orbital07

A Large-scale Parallel Computation for
Vibrational State Analysis Based on Quantum Monte Carlo method

calculations and Jastrow factor involving many body effects,

imaginary-time path or reptile . In Eqs.(5)-(7), we employed theare generally employed. In the present study in which we focus

following second order decomposition:

on quantum many body problems of molecular vibrations, weemployed the following Vibrational Self- Consistent Field (VSCF)
type wave function [5]:
where

,

,

and
where q i andof i th normal vibrational mode, respectively. Each modal function.
In Eq.(5),are a vibrational coordinate and modal functionis the propagator to generate

reptiles whose distributions are according to, andis the statistical weight of the generated reptile.

is expanded by a set of eigenfunctions of harmonic oscillator.
Variational parameters in the trial wave function such as centers,

In our RMC program, reptiles are generated with the Langevinexponents, and expansion coefficients of basis functions areequation defined as

optimized with the linear optimization method [6] proposed by
Umrigar et al . Optimized VSCF functions are used in the initial trial
wherewave function in RMC calculations described in the next section.
2.2. Reptation Monte Carlo (RMC) method

is the reduced mass of a vibrational mode,, anddrift velocity. In reject/accept procedures in RMC

Here, we consider a one-dimension case for the simplification ofimaginary-timetime evolution of time-dependent Schrödinger equation, thedistribution

are obtained by acting an imaginary-time propagator on a given initial trial wave function

calculations, we employed Metropolis’s method.
After a lot of Monte Carlo samplings with a long timemathematical expression. In QMC method using an imaginaryexact wave function,

the Gaussian distribution with a variance(0), the distributions of q and q
×density distribution:, and that of q

(N /2)(N)

converge to theconverge to the exact. To evaluate the expectation value ofHamiltonian operator, we used the following estimator (mixed

whereestimator) with the distribution

×:The exact density distribution

is used to calculate theis an imaginary-time. In RMC method, we consider thefollowing pseudo-partition function Z 0：

Decomposing the imaginary-time propagator

into N shorttime propagators with an imaginary-time step

expectation value of non-commutative operators with

, we

Hamiltonian such as potential energy operator, and to analyze

obtain the following partition function under the second order[7]

the geometry of molecules.

approximation (Suzuki-Torotter decomposition ):where

The variable q (i) is the configuration at the imaginary-timeand a set of configurations

08

is referred to as an,

Parallelization of QMC algorithm

3

As described in the previous sections, we generate multiple
Markov chains with Metropolis’s method for walkers and reptiles

Results and discussion

4

4.1. Parallel efficiency

in VMC and RMC calculations, respectively. Thus, the most

We performed a benchmark calculation of the parallel efficiency

straightforward and useful approach toward a parallel computingof our vibQMC program on TSUBAME 2.5 super-computer system

of QMC calculations is simultaneous samplings in configurationin Tokyo Institute of Technology: the number of cores is up to

space using statistical independent multiple walkers/ reptiles

5376 cores (1 process/core × 12 cores/node × 448 nodes).

(see Fig.1). In our program, multiple Markov chain generations

The vibrational ground state of a monohydrated positive core

are simply performed with multiple processes. Figure 2 showion, H5O2+, are used in benchmark calculations. We used thethe schematic illustration of our parallelized QMC program code.

analytical potential function proposed by Huang et al .[8] which

Data I/O from/to a external storage are performed in Master

well reproduces potential energy surface at CCSD(T)/aug-ccpVTZ

process. Markov chain generations for walkers/reptiles arelevel of ab initio calculations.

performed on both Master and Slave processes. We implement
such parallelized algorithm with MPI (Message Passing Interface),and use OpenMPI library version 1.4.2 on TSUBAME 2.5 supercomputer system.

Fig. 1 Schematic diagram of decomposed Markov chain.
Fig. 3 Relative speed up to the computational timewith 12 parallels (1 node) in VMC and RMC
calculations on TSUBAME 2.5 super-computer
system.
Figure 3 shows the relative speed up (=Time[12process] /Time [ N processes]) to the computational time with 12
parallels (1 node) in VMC and RMC calculations. The speed up with
5376 cores is 4868 (=405.70 × 12) in VMC calculations, and 2307(=192.25 × 12) in RMC calculations: the parallel efficiencies in
VMC and RMC calculations are about 91% and 43% , respectively.
Assuming that the size of problems is constant in each calculation(Amdahl’s law [9] ), we can estimate the ratio of sequential
Fig. 2 Schematic illustration of parallelized
QMC program code.

execution process in 12 parallels calculation as r （12）VMC =0.023% andr RMC =0.3% in VMC and RMC calculations, respectively.
（12）

09

A Large-scale Parallel Computation for
Vibrational State Analysis Based on Quantum Monte Carlo method

These ratio denote that 99.9981% and 99.975% of total
processing in a non-parallel (1 core) calculation are parallelized in
VMC and RMC calculations, respectively.
4.2. Monohydorated negative ion core: H3O2We performed anharmonic vibrational state analyses of amonohydrated negative core ion, H 3O2-, which is one of precursor
ions in forming aerosols in the atmosphere, as an applicationof large-scale parallel computation of our QMC program. Witha conventional ab inito calculation which does not include thenuclear quantum effect (NQE) or the effect of quantum molecularvibrations, the most stable equilibrium structure of H3O2- system
becomes asymmetric O—H*···O or O···H*—O structures as shownin the Fig. 4 (a), because the central hydrogen atom (H*) in H3O2has the double-well potential along the H* transfer coordinate.
On the other hand, with theoretical calculations including the
NQE [10,11], the effective potential energy curve along the H*

Fig. 4 Schematic illustrations of (a) equilibrium and(b) vibrational averaged geometries of H 3 O 2 - .

transfer coordinate changes to the single-well due to the small
potential barrier height (0.88 kJ/mol). Then, the H* is located at
the center between two oxygen atoms as O···H*···O as shownin Fig. 4(b) which corresponds to the transition state (TS) in
conventional ab initio calculations.
We analyzed the vibrational ground state andfundamental tone state of vibrational mode associated with the
H* transfer between two oxygen atoms (here we call bridgevibration) for H3O2- and its deuterium (D) and tritium (T) species.
The VSCF type wave function was used as the trial wave function
in VMC and RMC calculations. The analytical potential function
proposed by Huang et al . [8] which well reproduces potential
energy surface at CCSD(T)/aug-ccpVTZ level of ab initio calculation
are used.
Table 1 shows the zero-point vibrational energy (ZPE)
and the fundamental frequency of the bridge-vibration for all

Table 1. Zero-point vibrationalal energy (ZPE) andfundamental frequency (ω) of the bridge
vibration mode of H 3 O 2 - system. Unit in cm -1

species. The ZPEs, which are the variational energy, obtainedwith RMC method are lower than that with VMC method for all
species. The theoretical accuracy is, thus, improved with the RMC

In order to analyze structural properties of H*, D*, and T*
at each vibrational state, we focus on the parameter

calculations compared to the VMC calculations. In addition, theδ OH* = R O1H* − R O2H*, where the two oxygen-hydrogen distances,

RMC calculation well reproduce the experimental fundamental

R O1H*, O2H*, are defined in Fig. 4(b). The H* (or D*, T*) is located at

frequency (697cm-1) of the bridge- vibration of H-species within

the central position between two oxygen atoms at δ OH* = 0.

the error of 12 cm-1.

Figure 5 shows one dimensional distribution of δ OH* obtainedwith RMC calculations. In the vibrational ground states (Fig. 5(a)),
the H-species has a single peak in the distribution at δ OH* = 0.
Similar results have been reported with DMC calculations by
McCoy et al .[10] and path-integral molecular dynamics (PIMD)
calculation by Suzuki et al .[11] In their PIMD calculations, D- andT-species have almost the same distributions of δOD* and δ OT*
each other at 50 K. Our RMC calculation, however, shows that the10

distribution of T-species is more localized than that of D-species

equilibrium geometries. Each peak position of H-, D-, and T-species

at 0 K (vibrational ground state).

in one side of the distribution is shifted to a large |δOH* | region as
the nuclear mass becomes lighter due to the anharmonicity ofpotential energy curve along the H* (or D*, T*) transfer coordinate.
It is interesting that the probability density at δOH* = 0 ( TS withrespect to H* transfer) increases as the nuclear mass becomes
heavier at both the vibrational ground and fundamental tone
states.
Figure 5(c) shows the distributions of δOH* at 600K,
where we assumed Boltzmann statistics with two vibrational
states. The H-species has a single peak in distribution, but its
distribution at 600 K is more delocalized than that at 0 K. On theother hand, the distributions of the D- and T-species slightly split
as reported in the previous PIMD calculations at 600 K [11]. The twostate model with RMC calculations well reproduces PIMD results
at 600 K. The consistency between both methods indicates thatthe splits of the distribution of D- and T species at hightemperature region are mainly due the vibrational excitation of
the bridge-vibration.

Conclusion

5

In this report, we describe theoretical outline of vibrational
quantum Monte Carlo (vibQMC) method, variational Monte
Carlo (VMC) and reptation Monte Carlo (RMC) methods, and its
parallel implementation method toward a large-scale parallel
computation on a super-computer system. The parallelization
ratios of our vibQMC program code are 99.9981 % and 99.975 %in VMC and RMC programs, respectively. The parallel efficiencies
using 5376 cores on TSUBAME 2.5 super-computer system isabout 91% and 43 % in VMC and RMC calculations, respectively.
We also show theoretical results of vibrational state analysis of
the monohydrated negative core ion, H3O2-, which is a precursor
ion in forming aerosols in the atmosphere as an application
Fig. 5 One dimensional distribution of δOH* in H 3 O 2 and its D- and T-species. (a) the vibrational
ground state (v=0), (b) the fundamental tone
state of the bridge vibrational mode (v=1), (c)at 600K. Unit in Bohr.

of a large-scale parallel computation of our vibQMC program.
We used less time-consuming analytical potential functions invibrational state analyses from a computational cost point ofview. A further improvement of parallel efficiency is strongly
expected by a combination of vibQMC and on-the-fly ab initio
method in potential energy calculations. Such approach must be

The distribution of δ OH* at the fundamental tone state

quite useful for analyzing a larger molecular system, because anof the bridge-vibration is shown in Fig. 5(b). All species haveaccurate analytical potential function is no longer available fordouble peaks in the distribution at around two equivalentsuch systems.

11

A Large-scale Parallel Computation for
Vibrational State Analysis Based on Quantum Monte Carlo method

Acknowledgements
The present calculations were carried out on the TSUBAME 2.5
supercomputer in the Tokyo Institute of Technology as a subjectof the TSUBAME grand challenge (category B). Financial support
was provided by Grant-in-Aid for Scientific Research and for thepriority area by Ministry of Education, Culture, Sports, Scienceand Technology, Japan, for Y.K. and M.T.
References[1]

B.L. Hammond, W.A. Lester Jr. and P.J. Reynolds, “Monte
Carlo Methods in Ab Initio Quantum Chemistry” (World
Scientific, 1994).

[2]

S. Baroni and S. Moroni, Phys. Rev. Lett., 82, 4745 (1999).

[3]

D.G. Oblinsky, W.K. Yuen, S.M. Rothstein, J. Mol. Struct.
(THEOCHEM) 961, 219 (2010).

[4]

N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller
and E. Teller, J. Chem. Phys., 21, 1087 (1953).

[5]

J. M. Bowman, J. Chem. Phys. 68, 608 (1978).

[6]

J.Toulouse, C.J. Umrigar, J. Chem. Phys. 126, 084102 (2007).

[7]

M. Suzuki, Proc. Japan Acad., 69, 161 (1993).

[8]Huang et al. J. Am. Chem. Soc. 126, 5042 (2004) .

[9]

Gene M. Amdahl. Validity of the single processor approachto achieving large scale computing capabilities. In AFIPS
Conference Proceedings, pp. 483-485 (1967) .

[10] A. B. McCoy, X. Huang, S. Carter, and J. M. Bowman, J. Chem.
Phys. 123, 064317 (2005)[11] K.Suzuki, M. Shiga, and M. Tachikawa, J. Chem. Phys. 129,
144310 (2008).
[12] E.G. Diken, J.M. Headrick, J.R. Roscioli, J.C. Bopp, and M.A.
Johnson, A.B. McCoy J. Phys. Chem. A, 109, 8 (2005)

12

Large-scale DEM Simulations
for Granular Dynamics
Satori Tsuzuki Seiya Watanabe Takayuki Aoki
Global Scientific Information and Computing Center (GSIC) at Tokyo Institute of Technology

Granular materials such as sands and powders sometimes play important roles in science and engineering. In the numericalmethod DEM (discrete element method), the collisions between granular particles are described as classical spring force andfriction force models. Particle simulations using the DEM have been commonly applied to study granular materials. Evenin a spoon of sugar, there are more than 1 million grains. A large number of particles have to be used to simulate granularmaterials. Although coarse-graining models have been developed to reduce the particle number, it is quite meaningful tocarry out DEM simulations with realistic-sized particles for quantitative granular analyses. Recent supercomputers perform
well enough to carry out such granular simulations. Since granular material changes its particle distribution in time and space,
it is necessary to introduce a dynamic load balance in spite of the large computation overhead. After developing several
numerical methods for GPU implementation, we have succeeded in carrying out practical DEM applications using over 10
million particles. A golf bunker shot simulation demonstrates quite realistic results.

Introduction

1

The computational cost of the DEM is proportional
to the number of particles. We divide the particles into groups
with the same number of particles. If we divide the particles
by the numbers initially assigned to them, unacceptable data
communication among nodes will occur. It is natural to apply

Beside natural phenomena, granular materials often appear in

spatial domain decomposition to the DEM simulation sincemachines as toners in laser printers, pharmaceutical tablets inthe particles interact by being in contact with each other.

formulation processes, suspension flows in chemical plants,

Unfortunately, since granular material changes its particle

and so on. Granular materials sometimes behave as a liquid anddistribution in time and space, we do not keep the same particle

sometimes behave as a solid.

number in each subdomain with static domain decomposition.

There is a lot of demand for studying granularmaterials by computer simulation. In the commonly used DEM

Therefore, it is necessary to introduce a dynamic load balance.
The GPU has the advantages of high performance(Discrete Element Method), spring and friction forces workfor floating-point operation and wide memory bandwidth.

only among particles that are in contact. Since the number

We have to use the device memory on the GPU board, so theof contact interactions with a particle is small, the cost of thecommunication cost among the device memory becomes large.

memory access is greater than that of floating-point operation

This makes large-scale DEM simulations more difficult on GPU

in interaction calculation. On the other hand, MD (molecularsupercomputers.

dynamics) or astrophysics N-Body problems are similar types ofparticle simulations; however, the cost of floating-point operationis quite dominant due to their long-range interactions.
There are few studies on large-scale granularsimulations because the computational cost proportionally

DEM computation on GPU

increases with increasing number of particles. When we want to2

simulate granular materials in a spoon of sugar, there are morethan 1 million sugar grains, and the computational cost becomes

Recent GPUs have more than 2,000 processing cores (CUDA

quite large. Although coarse-graining models have beencores) on a single chip. Fine-grain parallelization and multi-

developed to reduce the particle number, it is very meaningful

thread programming are required to have high performance,

to carry out DEM simulations with realistic-sized particles forand we have to take into consideration the hierarchical structure

quantitative analyses. Recent supercomputers perform welland the limited size of the memory. In the GPU implementation

enough to carry out granular simulations using the real particleof our code, we use the CUDA programming framework given

number. Almost all supercomputers consist of multiple nodes,

by NVIDIA.

each node has a few processors with or without accelerators such

In the DEM, the particle interaction is modeled as aas the GPU (graphics processing unit), and are interconnected byspring and a dumping force proportional to the penetration

high-speed networks. Large-scale DEM simulations have to run

depth and the relative velocity of the two particles in contact

efficiently on the memory-distributed system of supercomputers.

with each other in the normal direction. In the tangential
direction, friction is also taken into account as shown in Fig. 1.

13

Large-scale DEM Simulations for Granular Dynamics

The linked-list method is a candidate to reduce thememory use as shown in Fig. 2 [1][2]. Each particle has a memory
pointer referring to the next particle in the same cell. Using chain
access we can reduce the memory usage to 1/8.

Fig. 1 DEM computational model.

The equation of translational motion is describedin classical mechanics in Eq. (1), and the right-hand side is thecontribution from all the particles in contact.

Multi-GPU DEM simulation using
Dynamic Load Balance

3

In large-scale DEM simulations requiring a lot of GPUs, thecomputational domain is decomposed into subdomains. A GPU
is assigned to each subdomain and computes particles locatedin the subdomain. Since particle distributions change in timeand space, static domain decomposition does not keep the same
number of particles in each subdomain. The slice-grid method [3]where x i is the position of the i -th particle and x ij is the penetration
depth between the i -th and j -th particles. The notations k and γ
are the spring constant and dumping coefficient, respectively.
In rotational motion, similar equations for the angular velocities
with torques are solved, and the particle positions and velocities
are updated by time integration with the leap-frog scheme oris introduced to maintain equal numbers of particles to keep
the memory usage equal and the computational load balanceamong GPUs. Figure 3 illustrates that the vertical boundaries ofthe subdomains move first to keep the vertical load balance ofthe horizontal subdomain group, and the horizontal boundariesof the subdomains move individually next.

the Runge-Kutta method.
All the dependent variables of particles are allocatedto the device memory (so called “global memory” in CUDA
programming). In the thread assignment for CUDA cores, one
thread computes one particle motion by solving Eq. (1).
It is quite inefficient to make judgments on whether
particles are in contact for all the particles. Neighbor-particle
lists are commonly used to reduce the cost to find the particles
in contact; however, the amount of memory needed to save thelist in the cell often becomes a severe problem in large-scale
simulations.

Fig. 3 Dynamic load balance based on
two-dimensional slice-grid method.
To determine the moving distance of the previous
subdomain boundary, we have to count the particles locatednear the boundary. We propose an efficient way to find nearboundary particles on “global memory” without copying theparticle data to the host CPU memory. The subdomain is divided
with a proper space, as shown in Fig. 4, and we count theparticle number within thespace by means of the Thrust

library. The particles in the neighbor subdomain after moving
Fig. 2 Neighbor particle search using linked-list
method on GPU.

the boundary are transferred there through the PCI-Express bus.
When the neighbor subdomain is allocated to a different node,
the data transfer includes the node-to-node communication by
the MPI library.

14

Strong and Weak scalabilities on
TSUBAME 2.5

4

We examined the performance of our DEM simulation on the
TSUBAME 2.5 supercomputer with GPUs (NVIDIA Tesla K20X) withthe dynamic load balance based on the 2-dimensional slice-grid
method. The scalabilities are studied for the agitation simulation
shown in Fig. 6 in the range of 2 M (million) to 1.0 billion particles.

Fig. 4 Particle counting and boundary
moving on GPU.

Frequent data transfer of particles among subdomains
causes fragmentation of the GPU memory, which degrades
the access performance and memory usage. In Fig. 5, adefragmentation should be executed with a proper frequency inspite of the overhead of data movement to the host memory.
Fig. 6 Agitation simulation using 4,120,000
particles on 64 GPUs.
In Fig. 7, the vertical axis indicates the performance
defined as the particle number divided by the elapsed time. The
solid lines indicate strong scalability with 2M (2 million), 16M,
and 129M particles, respectively. The square symbols show theresults with 2M particles using 4 to 64 GPUs, the triangles denote
16M particles using 32 to 256 GPUs, and the cross marks are with
129M particles using 256 to 512 GPUs. According to Fig. 7, theperformances keep improving in proportion to the number ofFig. 5 Defragmentation of particle data
on GPU device memory.

GPUs with 8 to 16-fold and become sluggish with more than
16 GPUs when using 2M particles.
We study the weak scalability by comparing theperformances for 2M, 16M, and 129M particles using 4 GPUs,
32 GPUs, and 256 GPUs, respectively. It was found that the
weak scalability was degraded from the ideal dashed line withincreasing GPU number. Some subdomains have shapes withhigh-aspect ratio and particles move across the boundary easily,so the amount of data communication increases and the total
performance becomes worse. It is meaningful that we succeeded
in a large-scale DEM simulation with 129M particles on 512 GPUs
regardless of low parallel efficiency.

15

Large-scale DEM Simulations for Granular Dynamics

Figure 9 shows the simulation for a granular conveyor
with 4,330,000 particles on 64 GPUs. The sand sliding down on
the spiral structure is demonstrated with 4,160,000 particles on
32 GPUs, as shown in Fig. 10.

Fig. 7 Scalabilities of DEM simulation on
TSUBAME 2.5.
Fig.9 Conveyor simulation using 4,330,000
particles on 64 GPUs.

Application to practical problems

5

To solve the problems when interacting with complex shapesof objects, the contact judgment is easily done by introducing
the signed-distance function from the object surface. When
we have the CAD data of the object, we have the distance fromthe object surface at the particle positions instead of having tocalculate the minimum distance from all the polygons of the
CAD data, as shown in Fig. 8[4].
Fig.10 Sand simulation on spiral slider using
4,160,000 particles on 32 GPUs.

As a typical simulation of granular structureinteraction, we apply our code to a simulation for a golf bunker
shot that had been studied by 2-dimensional simulation inprevious studies [5] due to the computational cost. We havesuccessfully performed a 3-dimensional DEM simulation byusing 16.7M particles with a realistic particle size on 64 GPUs of
TSUBAME 2.5, which has never been done before. To have theFig. 8 Representation of object shape using
Level Set method.

initial condition called “eye-ball”, we carried out a simulation of agolf ball falling on the sand in advance. We determine the swing
path of the sand wedge by using a model based on a rotationalor double pendulum. The swinging speed of the sand wedge
is 5.0 m/s at the head edge. Figure 11 shows a snapshot of thesimulation with 0.4-mm particles representing typical silica sand.
The golf ball is pushed by the sand without the iron head of thesand wedge hitting it.

16

Fig.11 Golf bunker shot with 16.7M particles on 64 GPUs.

DEM simulations using Non-sphericalparticles

6

Real granular particles are not spherical. Using the model ofnon-spherical particles, we can conduct more realistic granularsimulations [6]. It requires much higher computational cost andlarger memory even if we use a simple model in which several
spherical particles are rigidly connected. We carried out a largescale simulation for a foot stamp using 405,000 tetrapod-shaped
particles that consist of the same four spherical particles locatedat the vertex of a tetrahedron as shown in Fig. 12.
The numerical results are shown in Fig. 13 inFig.12 Tetrahedral non-spherical particle.

comparison with using spherical particles. The footprint withthe tetrahedral particles is clearly visible, and the shear friction isenhanced due to the inter-locking among tetrahedral particles.

17

Large-scale DEM Simulations for Granular Dynamics

Acknowledgements
This research was partly supported by KAKENHI, Grant-in-Aid for
Scientific Research (S) 26220002 from the Ministry of Education,
Culture, Sports, Science and Technology (MEXT) of Japan,
partly by the Japan Science and Technology Agency (JST) Core
Research of Evolutional Science and Technology (CREST) research
program "Highly Productive, High Performance Application
Frameworks for Post Peta-scale Computing", and partly by"Joint Usage/Research Center for Interdisciplinary Large-scale（ a ） Using spherical particles

Information Infrastructures (JHPCN)" and "High Performance
Computing Infrastructure (HPCI)" in Japan. The authors thankthe Global Scientific Information and Computing Center, Tokyo
Institute of Technology for use of the resources of the TSUBAME
2.5 supercomputer.
References[1]

G. S. Grest, B. D¨unweg, and K. Kremer, “Vectorized link cell
Fortran code for molecular dynamics simulations for a large
number of particles,” Computer Physics Communications,vol. 55, pp. 269–285, Oct. 1989.

[2]（ b） Using tetrahedral non-spherical particles

Gomez-Gesteira, M., Crespo, A., Rogers, B., Dalrymple, R.,
Dominguez, J., and Barreiro, A.: fSPHysicsg development ofa free-surface fluid solver Part 2: Efficiency and test cases,

Fig.13 Foot stamp simulations.

Computers and Geosciences, Vol. 48, No. 0, pp. 300-307(2012).
[3]

Conclusion

S. Tsuzuki and T. Aoki: Large-scale granular simulations
using Dynamic load balance on a GPU supercomputer, in

7

Poster at the 26th IEEE/ACM International Conference on
High Performance Computing, Networking, Storage andAnalysis (SC) 2014, New Orleans, US-LA, USA (2014).
[4]

J. A. Bærentzen and H. Aanæs, “Computing discrete signed

By introducing a dynamic load balance, we successfully

distance fields from triangle meshes,” Informatics andperformed large-scale DEM simulations with a maximum of

Mathematical Modeling, Technical University of Denmark,129M particles on a GPU supercomputer, TSUBAME 2.5. The

DTU, Richard Petersons Plads, Building 321, DK-2800 Kgs.
Lyngby, Tech. Rep., 2002.

two-dimensional slice-grid method works well as a dynamic
domain decomposition to keep equal memory consumption and[5]

5th International Conference on Engineering of Sport,

We applied the simulation code to several practical

September (2004)

problems including a golf bunker shot, and the scalabilities werealso examined. In addition, it is found that we can conduct more
realistic granular simulations by using non-spherical particles.

H. Horii, T. Koizumi, N. Tsujiuchi, M. Miki, and J. Hidaka:Computational Simulation Modeling of Bunker Shots,computational load balance.

[6]

Ikuya Ono, Hiroshi Nakashima, Hiroshi Shimizu, Juro
Miyasaka, and Katsuaki Ohdoi, Investigation of elemental
shape for 3D DEM modeling of interaction between soil anda narrow cutting tool, Journal of Terramechanics, Volume
50, Issue 4, August 2013, Pages 265-276, ISSN0022-4898.

18

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

19

vol.

International Research Collaboration
The high performance of supercomputer TSUBAME has been extended to theinternational arena. We promote international research collaborations using
TSUBAME between researchers of Tokyo Institute of Technology and overseas
research institutions as well as research groups worldwide.

Recent research collaborations using TSUBAME
1. Simulation of Tsunamis Generated by Earthquakes using Parallel
Computing Technique
2. Numerical Simulation of Energy Conversion with MHD Plasma-fluid Flow
3. GPU computing for Computational Fluid Dynamics

Application Guidance
Candidates to initiate research collaborations are expected to conclude
MOU (Memorandum of Understanding) with the partner organizations/departments. Committee reviews the “Agreement for Collaboration” for joint
research to ensure that the proposed research meet academic qualifications
and contributions to international society. Overseas users must observe
rules and regulations on using TSUBAME. User fees are paid by Tokyo Tech’s
researcher as part of research collaboration. The results of joint research are
expected to be released for academic publication.

Inquiry
Please see the following website for more details.
http://www.gsic.titech.ac.jp/en/InternationalCollaboration

13