6

Turbulence Simulation Using 40963
Vortex Particles on 4096 GPUs
An Ultra-fast Computing Pipeline for
Metagenome Analysis
with Next-Generation DNA Sequencers
GPU-Accelerated Large-Scale
Simulation of Seismic-Wave Propagation

Turbulence Simulation Using 40963
Vortex Particles on 4096 GPUs
Rio Yokota* Lorena Barba** Tetsu Narumi*** Kenji Yasuoka*****King Abdullah University of Science and Technology**Department of Mechanical Engineering, Boston University***Department of Computer Science, The University of Electro-Communications****Department of Mechanical Engineering, Keio University

The simulation of high Reynolds number turbulence is one of the most computationally demanding applications in highperformance computing. In the present study, we perform the simulation of turbulence in a periodic box using 40963
vortex particles. The accuracy of the calculation is verified by comparing with a pseudo-spectral code. The relative parallel
efficiency is compared between the FMM based vortex particle method, and the FFT based pseudo-spectral method on
up to 4096 GPUs on the TSUBAME 2.0 system.

Introduction

1

despite the 300 fold increase in LINPACK performance between
the Earth Simulator and K computer. This can be either interpretedas; a) the difficultly to parallelize FFT on large distributed memory
systems with torus networks, or b) the inability of LINPACK toserve as an indicator for the actual performance that can be

Turbulence is ubiquitous in our daily lives and alsoachieved in real applications. The former interpretation prompts

appears in various engineering problems involving energy andthe development of alternative algorithms that have better

environmental applications. The non-linearity of the governing

scalability. The latter perception would most likely facilitate an

equation and the large degrees of freedom in the system, makeearly adaptation of an alternative benchmark, which could leadthe simulation of turbulence one of the most computationallyto better co-design between the software and hardware in highdemanding problems in computational physics. Recent advances

performance computing.

in the computational capability of supercomputers has made it

In the present study, we propose an alternative method

possible to perform direct numerical simulations of the Navier-for solving isotropic turbulence that can scale to the full node of

Stokes equation for reasonably high Reynolds numbers.

TSUBAME 2.0. We caution the reader that ``scalability” does not

Although typical engineering applications for turbulent

necessarily reflect the ``performance” of the code – slow codes

flows involve complex boundary conditions and require strong

usually scale better because it is easier to hide the communication.

coupling with other domains of physics, there exist cer tain

By comparing both the time-to-solution and scalability of two

canonical flow fields where the physics of turbulence itself can bedifferent algorithms for the same application, we hope to provide

investigated without the complications of these extra constraints.

data that supports decision making in the choice of algorithms

The homogeneous isotropic turbulence under periodic boundaryand co-design of next generation software and hardware in highconditions is one such case, where the fundamental characteristics

performance computing.

of high Reynolds number turbulence can be analyzed in isolationfrom mean shear and near wall effects.
The numerical method of preference for solving isotropic
turbulence with periodic boundary conditions has been pseudo-[1]

There have been reports of spectral method calculations for 81923 grid points,

but these are benchmarks that are not run for a significantly long time, and do
not provide any data on the actual physics of turbulence

spec tral methods [1] .   T he largest of such calculations was
performed with 40963 grid points at a Taylor-microscale Reynolds
number of                        on the Earth Simulator back in 2002[2]. Theircalculations revealed that the scaling of the high order turbulence
statistics is significantly different at high       , while the local scaling

Vortex Method2

exponents seemed universal [3] . Such conclusions would nothave been attainable without the use of supercomputers and an

There are a variety of methods that are referred to as

idealized flow configuration, and symbolizes the importance of``vortex methods”[4]. Some are Lagrangian while others are semi-such efforts.

01

Lagrangian. Some are based on the vorticity-streamfunction

We would like to point out however, that the record for theformulation while other use the vorticity-velocity formulation.

largest turbulence calculation has not been broken for a decade1,

There are also many different ways to treat the viscous diffusion

in the Lagrangian frame of reference. A common advantage ofthese methods is the treatment of convection in the Lagrangianframe. This allows the use of larger time step sizes while eliminating

Fast Multipole Method

numerical diffusion and dispersion effects. Furthermore, the3vorticity-based formulation allows the points to be placed only inregions with non-zero vorticity. This leads to substantial savings in

The N-body problem resulting from the Biot-Savar t

computational time and storage for external flows with unsteadyequation and the stretching term of the vorticity equation results

vortex dynamics, such as vortex rings and trailing edge vortices.

in O(N ) operations if solved directly. We use the fast multipole

In this work, we use the full Lagrangian vortex method with
[5]2

method (FMM)[6] to bring the number of operations down to O(N) .

vorticity-velocity formulation . The vorticity field is discretized bySince the Biot-Savart equation originates from a Poisson’s

a superposition of Lagrangian vortex particles with Gaussian basis

equation, the FMM for Laplace kernels can be used here. The curl

functions. In the vorticity-velocity formulation the velocity Poissonin the source term of the Biot-Savart equation and the gradient inequation.

the stretching term each require an extra a derivative to be applied
to the Laplace kernel. However, since the multipole expansion andlocal expansion already consist of these derivative terms, it is trivial

is solved along with the vorticity equation

to extract this information without any extra computation.
3.1 Load Balancing
When parallelizing FMMs, the global nature of the N-body

The velocity Poisson equation becomes the Biot-Savart equationinteraction and the dynamic nature of the particle distribution

when formulated as an integral equation with the Green’s functionpose a unique challenge. Unlike grid based methods, the use ofas the kernel. This becomes an N-body problem that can be solved

domain decomposition will not localize the data dependency

using the fast multipole method (FMM) described in the next

–the halo region will be the entire domain but with hierarchicallysection. In the vorticity equation, the convection term is calculatedcoarsened far regions. A common approach is to use a global tree

by advecting the Lagrangian particles. The stretching term is

structure to minimize the communication[7]. When the particle

calculated by substituting the Biot-Savart equation into the velocity,

distribution is irregular and the tree structure is adaptive, equallywhich results in another N-body problem. The diffusion term ispartitioning the domain will result in load imbalance. Furthermore,calculated by increasing the size of the Gaussian basis function

the dynamic nature of the N-body simulation necessitates frequentaccording to the analytical solution of the diffusion equation.

repartitioning, so the overhead must be small.

In order to maintain suf ficient overlap between the vor tex

A clever way to solve the load balancing issue is to record

particles, the coordinates are reinitialized every few steps, and thethe load imbalance from the previous time step, and use it tocorresponding vortex strength of each particle is calculated torepartition the tree for the present step. Such strategies werereproduce the vorticity field before the reinitialization. This is done

first introduced in the early 90’s on both shared memory [8] andusing radial basis function interpolation with a Gaussian basis.

distributed memory architectures[9]. This idea of using informationfrom the previous step to update the present step can be usedwith any existing partitioning scheme such as orthogonal recursive

02

Turbulence Simulation Using 40963
Vortex Particles on 4096 GPUs

bisection (ORB)[7], partitioning of Morton/Hilbert keys[9], and graphinteraction (P2P) kernel for the near field and the multipole-to-local

based partitioning [10]. The work/communication load from the(M2L) or multipole-to-particle (M2P) kernel for the far field.

previous step can be used as weights in any of these originallyThe ratio between these near field and far field calculations must

unweighted partitioning schemes. Our current FMM code supports

be optimized by properly selecting the number of particles perboth ORB and Morton key based partitioning schemes.

leaf cell during the tree construction. This is a non-trivial task on

Once the domain is partitioned, each process will own aheterogeneous architectures, since each kernel will accelerate

local portion of the global tree structure. It is then necessary toat different rates on GPUs [14]. There are also a wide variety ofcommunicate the local essential tree (LET), which is a significantly

mathematical formulations for each of these kernels, some of whichsmaller subset of the global tree that is required for the tree

are better suited for a particular architecture.

traversal on each process. The communication of the LET is

A simple strategy for selecting the optimal kernel on any

global in nature, and is a potential bottleneck for the scalability ofgiven architecture is to measure the time of dry kernel calls at the

FMMs. A common technique in both treecodes[11] and FMMs[12] isbeginning of the simulation, and to use those timings to select theto repeatedly exchange the LET on the other side of the bisector

fastest option during the tree traversal. By providing the optionof the equivalent binary tree structure. This N-D hypercube typeto select between P2P (direct), M2P (treecode),   and M2L (FMM)

communication is extended in our code to work with number ofkernels, the resulting algorithm will always run faster than any ofprocesses that are not a power of two. We have further extended

these methods alone, on any architecture, for any problem size. Our

this framework to work with periodic boundary conditions.

tests show that this is indeed the case[15].

3.2 Dual Tree Traversal
Both treecodes and FMMs build the same tree structure
and partition it in the same way. The difference between the two
methods lies in the way the tree is traversed. Treecodes loop over
the target leaf cells and traverse the entire source tree for eachtarget leaf cell. FMMs usually do not traverse the tree, but loop over
all target cells and construct an explicit source cell list per target
cell. For an adaptive tree, the construction of such source cell lists
becomes complicated, especially if one were to introduce theconcept of multipole acceptance criteria (MAC) into FMMs.
The dual tree traversal [13] can be used to simplif y theconstruction of source cell lists in FMMs. It can be thought of as a
generalization of U,V,W,X-lists[12] that also allows the use of MAC in
FMMs. A conceptual diagram of the dual tree traversal is shown in
Figure 1. Two trees are traversed simultaneously, one for the target
and one for the source. The MAC is applied for each pair during
the dual traversal. If the MAC is satisfied, the multipole-to-local
translation is calculated and the pair is removed from the traversal
queue. If the MAC is not satisfied, the traversal continues until both
trees reach their leaf cells, at which point the direct interaction
will be calculated between all particles in the cells. In other words,
the dual tree traversal provides a simple bookkeeping strategy
for constructing a MAC-based interaction stencil that is mutually
exclusive at each level, by letting the target cells inherit a uniquelist of source cells from their parents.
3.3 Auto-tuning on CPU/GPU
In hierarchical N-body methods such as treecodes andFMMs, the dominant part of the calculation is the particle-particle

03

Figure 1  Diagram of the dual tree traversal

Scalability Results

4GPUs. However, this part of the calculation scales perfectly, anddoes not affect the scalability of the FMM. The parts that do affectthe scalability are the tree construction and MPI communication.
Actually, the tree construction involves MPI communication forthe partitioning so it is the efficiency of MPI communication that

The present simulations are run on the TSUBAME 2.0 system,

determines the weak scalability.

which has 1408 nodes equipped with two Intel Xeon5670 CPUs,

It may be worth noting that the N-D   hypercube typethree NVIDIA M2050 GPUs, 54 GB of RAM, and 120 GB of local SSD

communication of the LET (described in section 3.1) turned out tostorage. The inter-node connection is a full bisection/non-blocking

be slower than a simple call to MPI_Alltoallv for sending the entire

fat-tree, where each node has two QDR Infiniband links.

LET at once. Therefore, the results shown in Figure 2 are with MPI_

The breakdown of the weak scaling runs for N=40962

Alltoallv and not the hypercube communication.

particles per GPU is shown in Figure 2. We assign one MPI process

We performed a similar weak scalability test for the spectral

per GPU so the number of processes is equal to the number

method, and compared the parallel efficiency with the FMM inof GPUs. The execution time is for the combined Biot-Savart +Figure 3. The spectral method was calculated on the CPU because

stretching kernel. In the legend, “P2P evaluation” is the evaluation

we did not have access to a GPU based FFT code that gave us

time of the P2P kernel in the near field, “FMM evaluation” issignificant acceleration over FFTW (when host-device data transfer

the sum of all other FMM kernels, “MPI communication” is theis included) at the time of the comparison. The parallel efficiency of

total time spent on all MPI communications except the ones

the FMM was 74 % at 4096 GPUs, while the parallel efficiency of thein the tree construction, “GPU buffering” is the time spent onspectral method was 14 % on the same number of CPUs. The time-

reordering/buffering the data and sending it to the GPU, andto-solution was approximately 100s for both the FMM and spectral“Tree construction” is the parallel construction/ partitioning of themethod.

global tree. The local portion of the P2P evaluation is overlappedwith the MPI communication. In Figure 2 the MPI communication
time is subtracted from the P2P evaluation time so that the totalheight of the bar reflects the actual wall clock time. The largest run
with N=40963 vortex particles on 4096 GPUs results in a sustained
performance of 1.01  Pflop/s.
We see that the GPU buffering is taking a significant amountof time, but we have found this to be necessary in order to achieve
high efficiency in the P2P evaluation and FMM evaluation on

Figure 2  Weak scaling of FMM on up to 4096 GPUs for                  N=40962  particles per GPU.

Figure 3  Parallel efficiency of FMM and spectral method                  on up to 4096 processes.

04

Turbulence Simulation Using 40963
Vortex Particles on 4096 GPUs

Isotropic Turbulence5

the calculaiton. The duration of the run was severely limited by thewindow of access to the full node of TSUBAME 2.0. Therefore, theresults of the present turbulence simulation are only preliminary,
and are only interesting in terms of demonstrating the scalability
and validity of FMM based vortex methods on large scale GPU

The calculation of isotropic turbulence for             wasperformed by the FMM based vortex particle method and FFT
based pseudo-spectral method. The calculation domain was aperiodic box of [-π,π]3 and the number of particles/grid-points was
40963. For the FMM, we used 273 periodic images to approximate

systems, which is orthogonal to the issue of what new physics ofturbulence was obtained.
The kinetic energy spectra for the pseudo-spectral method
and vortex method at t/T=2 are compared in Figure 5. There is asmall deviation at the high end of the wave numbers, but otherwise
the results agree quantitatively. The improvement over similar

the periodic boundary condition, and the order of multipole

comparisons in the past[5][16] was achieved by using a higher order

expansion was set to p=14. The initial condition was generated

of multipole expansion (p=14), using more periodic images (273),

using the spectral code, by generating a velocity field in wave

more frequent reinitialization (every 5 steps), higher accuracy in thespace with a prescribed energy spectrum and random phases.

RBF iteration (relative |L|2=1e-5), and of course much more vortex

This initial velocity field is used directly as an input to the pseudo-

particles. We have confirmed that all the conditions mentioned

spectral method. For the vortex method, the velocity field is first

above must be satisfied in order to achieve the high accuracy inconverted to real space and the strength of the vortex particles arevortex methods that we show in Figure 5.

adjusted to reproduce this velocity field. The vortex particles are
initially positioned at the cell centers (same location as the pressure
on a staggered grid), and radial basis function interpolation is used
to determine the strength of each particle. The initial core radius of
the vortex particles are set to σ=Δx so that the overlap ratio is 1.
Figure 4 shows an isosurface of the second invariant of thevelocity gradient tensor at time t/T=2, where T is the large eddy
turn-over time. Many small vortex structres can be observed, but
no large scale coherent structures are seen at this early stage of

Figure 4  Isosurface of the second invariant of the velocity                  gradient tensor from the vortex method simulation

05

Figure 5  Energy spectra of the pseudo-spectral spectral                  method and vortex method at t/T=2

Conclusions

6

Comput. Phys., Vol. 226, pp. 1589-1606 (2007)[6]

H. Cheng, L. Greengard, and V. Rokhlin: A Fast Adaptive
Multipole Algorithm in Three Dimensions, Vol. 155, pp. 468498 (1999)[7] M . S. Warren and J. K . Salmon: A s trophysical N - b o dy
Simulation Using Hierarchical Tree Data Structures, Proc. ACM/

We compared the scalabilit y of a FMM based vor tex

IEEE Supercomput. Conf., pp. 570-576 (1992)

method against a FFT based pseudo-spectral method for thesimulation of isotropic turbulence with 40963 particles/grid points[8]

A d a p t i v e Fa s t M u l t i p o l e M e t h o d , P r o c . A C M / I E E E

on 4096 CPU/GPUs on TSUBAME 2.0. The FMM based calculation

Supercomput. Conf., pp.54-65 (1993)

achieved 74 % parallel efficiency, while the FFT based method had
14 % parallel efficiency on 4096 processes. The calculation using
40963 vortex particles reached 1.01 Pflop/s, which is a new record
for FMMs as far as the authors are aware. The comparison of theresults between the vortex method and pseudo-spectral method

J. P. Singh, C. Holt, J. L. Hennessy, and A. Gupta: A Parallel[9]M. S. Warren and J. K. Salmon: A Parallel Hashed Oct-tree
N-body Algorithm. Proc. ACM/IEEE Conf. Supercomput., pp.
12-21 (1993)[10] S.-H. Teng. Provably Good Partitioning and Load Balancing

showed quantitative agreement in the energy spectra. However,

Algorithms for Parallel Adaptive N-body Simulation. SIAM J.

due to the limited availability of the full node of TSUBAME 2.0, we

Sci. Comput., Vol. 19, pp. 635-656 (1998)

were not able to compare high order turbulence statistics between[11] T. Hamada, K. Nitadori, K. Benkrid, Y. Ohno, G. Morimoto, T.

the vortex method and spectral method. Future generation

Masada, Y. Shibata, K. Oguri, and M. Taiji: A Novel Multiple-

supercomputers may offer us the opportunity to make use of ourwalk Parallel Algorithm for the Barnes-Hut Treecode on GPUs –FMM based turbulence solver.

Towards Cost Effective, High Performance N-body Simulation,
Comput. Sci. – Res. Dev., Vol. 24, pp. 21-31 (2009)

Acknowledgements[12] I. Lashuk, A. Chandramowlishwaran, H. Langston, T.-A. Nguyen,
R. Sampath, A. Shringarpure, R. Vuduc, L. Ying, D. Zorin, andComputing time in the TSUBAME 2.0 system was made possible by

G. Biros: A Massively Parallel Adaptive Fast Multipole Method

the Grand Challenge Program of TSUBAME. The current work was

on Heterogeneous Architectures, Proc. Conf. High Perform.

partially supported by the Core Research for the Evolution Scienceand Technology (CREST) of the Japan Science and Technology
Corporation (JST). LAB acknowledges partial support from NSF
grant OCI-0946441, ONR award #N00014-11-1-0356, and Boston
University College of Engineering.

Comput. Netw. Stor. Anal., (2009)[13] W. Dehnen: A Hierarchical O(N) Force Calculation Algorithm, J.
Comput. Phys., Vol. 179, pp. 27-42 (2002)[14] R. Yokota and L. A. Barba: Treecode and Fast Multipole Method
for N-body Simulation with CUDA, GPU Computing Gems,
Morgan Kaufmann (2011)

References[15] R. Yokota and L. A. Barba: Hierarchical N-body Simulations with
Auto-tuning for Heterogenous Systems. Comput. Sci. Eng., Vol.

[1]

S. A. Orszag and G. S. J. Patterson: Numerical Simulation ofThree Dimensional Homogeneous Isotropic Turbulence. Phys.

[2]14, pp. 30-39 (2012)[16] R. Yokota, T. Narumi, R. Sakamaki, S. Kameoka, S. Obi, and K.

Rev. Lett. Vol. 28, pp. 76-79 (1972)

Yasuoka: Fast Multipole Methods on a Cluster of GPUs for the

M. Yokokawa, K. Itakura, A. Uno, T. Ishihara, and Y. Kaneda:

Meshless Simulation of Turbulence. Comput. Phys. Comm., Vol.

16.4-Tflops Direct Numerical Simulation of Turbulence by a180, pp. 2066-2078 (2009)

Fourier Spectral Method on the Earth Simulator, Proc. ACM/IEEE Supercomput. Conf., Washington DC, pp. 50-66 (2002)[3]

T. Ishihara, T., Gotoh, and Y. Kaneda: Study of High-Reynolds
Number Isotropic Turbulence by Direct Numerical Simulation,
Annu. Rev. Fluid Mech., Vol. 41, pp. 165-180 (2009)[4]

G. H. Cottet and P. Koumoutsakos: Vortex Methods, Cambridge
University Press (2000)[5]

R. Yokota, T. K. Sheel, and S. Obi: Calculation of Isotropic
Turbulence Using a Pure Lagrangian Vor tex Method, J.

06

An Ultra-fast Computing Pipeline for
Metagenome Analysis with
Next-Generation DNA Sequencers
Takashi Ishida* Shuji Suzuki* Yutaka Akiyama**Department of Computer Science, Tokyo Institute of Technology

Metagenome analysis is useful for not only understanding symbiotic systems but also watching environment
pollutions. However, metagenome analysis requires sensitive sequence homology searches which require large
computation time and it is thus a bottleneck in current metagenome analysis based on the data from the latest DNA
sequencers generally called a next-generation sequencer. To solve the problem, we developed a faster homology
search program based on GPU-calculation and a large-scale computing pipeline for metagenome analysis on
TSUBAME2.

sequencer in a few hours. We believe our pipeline will accelerate

1

Introduction

metagenome analyses with next-generation sequencers.

Metagenome analysis is the study of the genomesof unculture d microb es obtaine d dire c tly f rom microbial

Metagenome mapping

communities in their natural habitats such as soils, seas, and2human bodies. The analysis is useful for not only understanding
symbiotic systems but also watching environment pollutions[1].

Current sequencers produce only information in short

These days, the latest DNA sequencers, generally called next-

fragments, whose lengths range between 50 and 700 base-pair

generation se quencer, b e came to pro duce huge amount(bp). Thus, mapping or assembling processes are required foro f g en o mi c dat a in a sh o r t tim e an d i t is e x p e c te d thatobtaining meaningful results from the output.

metagenomic researches are promoted based on such huge

DNA sequence fragment on known genome information. In single-

genomic data.
However, metagenome analysis requires comparisonsorganism genomics, it is easy because the reference genome hasof sequence data obtained from a sequencer with sequencealready been obtained and we do not have to mention about many

data of remote homologues in databases. Because, current

mismatches and gaps. Thus, the process is similar to general text

databases do not include sequence data for most of microbes

search and does not require large computation. Therefore, by using

in the sample. Therefore, sensitive sequence homology searchefficient short-read mapping programs such as BWA[5] and Bowtie[6],processes are required in metagenome analysis. Unfortunately,

we can process the output of a next generation sequencer on athis pro cess ne e ds large computation time and is thus a

few workstations.

bottleneck in current metagenome analysis based on the data

In contrast, in metagenomic analysis, the DNA sequence
fragments obtained from environmental samples frequently[2]from a next-generation sequencer .

07

Mapping is a process to identify the location of each

Here, we developed a large -scale automatedinclude DNA sequences from many different species, and closelycomputing pipeline for analyzing huge amount of metagenomic

related reference genome sequences are often unavailable. Thus,

data obtained from a next-generation sequencer. This pipeline

more sensitive search, generally called a homology search, isenables us to analyze metagenomic data from a next generation

required for the identification of novel genes. In addition, in thesequencer in real time by utilizing huge computational resources

typical metagenomic analyses, sequenced DNA fragments areon TSUBAME2. Also, we developed a fast homology searchtranslated into protein coding sequences and then further assigned[3]

program based on GPU-calculation, which is named GHOSTM .

to protein families to improve search sensitivity. This also increases

This program has sufficient search sensitivity for metagenomic

the computational cost of homology search.

analysis and is much faster than BLASTX program [4] , which is

The BLAST algorithm is sufficiently sensitive for searching

generally used for previous metagenomic researches. By using

protein families and is much faster than the classical dynamic

these new program and pipeline, we can process metagenome

programming method. Thus, BLASTX program has been used in

information obtained from a single run of a next generation

previous metagenomic researches[7]. However, its performance isinsufficient for analyzing the large quantities of data produced by alocal alignments around the candidate positions using thenext-generation sequencer. In practice, by using BLASTX program,Smith-Waterman algorithm [8] for calculating the alignment

approximately 25,000 CPU days are needed for querying 600 GBscores. Finally, the third component sorted the alignment scoresof short reads generated from one run of the latest sequencer

and output the search results.

Illumina’s Hiseq 2000. Thus, the metagenomic mapping based

B ot h t h e c a n d i date s ea rc h a n d l o c a l a li g nm e nt

on a homology search is a bottleneck in current analysis and its

components required a large amount of computing time.

acceleration is highly required.

Therefore, we processed queries on both components in paralleland mapped them onto GPUs. Thus, multiple queries were
simultaneously processed on different GPU cores.

GHOSTM: a GPU-accelerated homology
search tool for metagenomics

3

3.2 Search sensitivity and speed
Because metagenome analyses require highly
sensitive searches, it is difficult to use homology search programwith high speed but low sensitivity, such as BLAT [9]. In contrast,

To a d d re s s t h e is su e, w e d e ve l o p e d a n e w a n d

GHOSTM has suf ficient search sensitivit y for metagenomic

efficient homology search algorithm suitable for GPU calculation

analysis. Figure 2 shows the comparison of search sensitivity

and implemented the system on GPUs. The system accepts

for each homology search program. To evaluate the searcha large number of short DNA fragment sequences produced

sensitivity, we used the search results obtained with the Smith-by a nex t- generation sequencer as the input and, like theWaterman local alignment method implemented in SSEARCH[10],BLASTX program, performs DNA sequence homology searches

and these results were assumed to be the correct answers.

against a protein sequence database. We used NVIDIA CUDA to

We analyzed the performance of a particular method in terms

implement the GPU-calculation. Thus, the system requires CUDAof the fraction of its results that corresponded to the correct

2.2 or higher. The search system, which we named GHOSTM,

answers obtained by SSEARCH. The search accuracy of GHOSTM

demonstrated a calculation speed that was 130 times faster withwas clearly higher than BLAT. Low-scoring hits (e.g., <50) areone GPU than BLAST on a CPU.
3.1 Algorithm
The GHOSTM is mainly composed of threecomponents, as shown in Figure 1. The first component searched

generally not used in practice because such hits can occur
by chance. With the exception of the low-score hits, GHOSTM
successfully identified more than 90% of the hits identified

by SSEARCH. This result suggests that GHOSTM is sufficiently
accurate for general usage.

the candidate alignment positions for a sequence from theTable 1 shows the computational times of BLAST, BLAT,

database using the indexes. The second component calculated

and GHOSTM for 100 thousand reads. Each query read has the

Figure 1  Data flow and processing within GHOSTM.

Figure 2  Comparison of search accuracy08

An Ultra-fast Computing Pipeline for
Metagenome Analysis with
Next-Generation DNA Sequencers

length from 60 to 75 bp and the search target is KEGG Genes (data has 224 million DNA reads with 75 bp, and after excluding“genes.pep”) database with approximately 2.5GB. The GHOSTM

low-quality data, the dataset became 71 million DNA reads.

program achieved a calculation speed approximately 130 andThus, homology searches were performed with 71 million DNA

400 times faster than the BLAST program using 1 thread and 4

reads for NCBI nr amino-acid sequence DB (4.2GB). We evaluatedthreads, respectively. Moreover, GHOSTM was approximately 3.4

the effective performance of the pipeline with both BLASTX on

times faster than BLAT despite of its higher search sensitivity.

CPUs and GHOSTM on GPUs.

GHOSTM achieves both high search speed and high search

As results, the pipeline shows almost linear speedup

sensitivity compared with previous homology search programs.

to the number of computing cores. When we used BLASTX asa homology search program, the pipeline achieved to process
about 24 million reads per an hour with 16,008 CPU cores (1,334
computing nodes) (Figure. 3). When we used GHOSTM as ahomology search program, the pipeline achieved to process
about 60 million reads per an hour with 2,520 GPUs (8 4 0
computing nodes) (Figure. 4). When we used 2,520 GPUs, thespeed up is only 20% compared with the case that we used 1,260
GPUs. This was because the dataset was too small for 2,560 GPUs
and it might cause the failure of the load balancing. Therefore,
the speedup will be linear even with more than 2,520 GPUs byprocessing much more metagenomic data.

Table 1 Comparison of search speed

Automated computing pipeline on
TSUBAME2

4

Even using GHOSTM, it is dif ficult to analyze hugeamount of metagenomic data on a single workstation. Thus,
we developed a large-scale automated computing pipeline foranalyzing huge amount of metagenomic data obtained from anext-generation sequencer. This pipeline enables us to analyze
metagenomic data from a next generation sequencer in real

Figure 3  Speedup of the BLASTX-based system                  for the number of CPU corestime by utilizing huge computational resources on TSUBAME2.
File I/O processes including a database copy and writing
search results caused a contention problem when we used
many computation nodes. Thus, we changed to employ asophisticated file transfer manner where data are simultaneously
copied from local disk of a node to another in a binary-tree
manner.
4.1 Large-scale experiments
By the assistance of TSUBAME Grand Challenge (Large
Scale GPU application) program, we performed a large-scale
metagenome analysis by using our pipeline on whole TSUBAME2
system. We used data sampled from polluted soils and obtainedby using a next-generation sequencer. Original metagenomic

09

Figure 4  Speedup of the GHOSTM-based                  system for the number of GPUs

Summary5

microbiome with increased capacity for energy harvest.,”Nature, vol. 444, no. 7122, pp. 1027-31, Dec. 2006.
[8]

T. F. Smith and M. S. Waterman, “Identification of common
molecular subsequences.,” Journal of molecular biology,
vol. 147, no. 1, pp. 195-7, Mar. 1981.

For analyzing metagenomic data obtained from a next
generation sequencer in real time, we developed a new efficient[9]

W. J. Kent, “BLAT--the BLAST-like alignment tool.,” Genome
research, vol. 12, no. 4, pp. 656-64, Apr. 2002.

GPU-based homology search program GHOSTM and a large-[10] W. R . Pearson, “ Searching protein se quence libraries:

scale automated computing pipeline which enable us to utilize

comparison of the sensitivity and selectivity of the Smith-huge computational resources on TSUBAME2.

Waterman and FASTA algorithms.,” Genomics, vol. 11, no. 3,

The results of the experiment with whole TSUBMAME2 system

pp. 635-50, Nov. 1991.

indicate the pipeline can process genome information obtained
from a single run of next generation sequencers in a few hours.
We believe our pipeline will accelerate metagenome analysis
with next generation sequencers.
Acknowledgements
This research was supported in part by HPCI Strategic
Program Computational Life Science and Application in Drug
Discovery and Medical Development by MEXT of Japan, Cancer
Research Development funding by National Cancer Center,
Japan, and the CUDA COE Program by NVIDIA.  
References[1]

M . Ar umugam et al., “Enterot y p es of the human gut
microbiome.,” Nature, vol. 473, no. 7346, pp. 174-80, May
2011.

[2]

J. C. Wooley, A. Godzik , and I. Friedberg, “A primer on
metagenomics.,” PLoS computational biology, vol. 6, no. 2,
p. e1000667, Jan. 2010.

[3]

S. Suzuki, T. Ishida, K. Kurokawa, and Y. Akiyama, “GHOSTM: A
GPU-Accelerated Homology Search Tool for Metagenomics,”
PLoS ONE, vol. 7, no. 5, p. e36060, May 2012.

[4]

S. F. Altschul, W. Gish, W. Miller, E. W. Myers, and D. J. Lipman,“Basic local alignment search tool.,” Journal of molecularbiology, vol. 215, no. 3, pp. 403-10, Oct. 1990.

[5] H. Li and R. Durbin, “Fast and accurate short read alignment
with Burrows-Wheeler transform.,” Bioinformatics (Oxford,
England), vol. 25, no. 14, pp. 1754-60, Jul. 2009.
[6]

B. Langmead, C. Trapnell, M. Pop, and S. L. Salzberg, “Ultrafast
and memory-efficient alignment of short DNA sequences
to the human genome.,” Genome biology, vol. 10, no. 3, p.
R25, Jan. 2009.

[7] P. J. Turnbaugh, R. E. Ley, M. a Mahowald, V. Magrini, E.
R. Mardis, and J. I. Gordon, “An obesit y-associated gut

10GPU-Accelerated Large-Scale
Simulation of Seismic-Wave Propagation
Taro Okamoto* Hiroshi Takenaka** Takeshi Nakamura*** Takayuki Aoki*****Department of Earth and Planetary Sciences, Tokyo Institute of Technology **Department of Earth and Planetary Sciences, Kyushu University***Earthquake and Tsunami Research Project for Disaster Prevention, Japan Agency for Marine-Earth Science and Technology****Global Scientific Information and Computing Center, Tokyo Institute of Technology

Simulating seismic wave propagation is important for the study of earthquake sources, the generation of strong
ground motions and the excitation of large tsunamis. We describe methods for accelerating large-scale finitedifference time-domain simulation of the seismic wave propagation by the use of graphics processing units (GPUs).
We then present examples of the wave-filed from the 2 011 Tohoku-Oki earthquake simulated by using the GPUsof TSUBAME supercomputer. The simulated wave-field exhibits strongly complex pattern reflecting the sourcecomplexity and the heterogeneities around the source region.

Introduction

1

The Tohoku-Oki earthquake on March 11, 2011 (M JMA
9.0; Fig. 1) generated strong shaking reaching the maximum
intensity (seven) on the JMA’s scale and caused devastatingtsunamis with run-up heights exceeding 30 m. 15,845 people
were lost, and 3,368 people are still missing (January 27, 2012).
Such magnitude 9 earthquake was not expected to occur along
the plate interface off the northeastern Japan. Thus it is very
important to study this event for understanding the geophysicalcondition of the generation of mega-thrust earthquake, thecharacteristics of the induced strong ground motions, and themechanism of the excitation of the large tsunamis.
The ground motion records of this event are quite
important data for the quantitative studies on the earthquake
source and the induced damages. However, modeling of theground motions is not a simple task because of the strong lateralheterogeneity in and around the Japan trench: steeply varying
topography, oceanic water layer, thick sediments, crust withvarying thickness and subducting oceanic plate can all affectFigure 1  
Aftershock distribution of the 2011 Tohoku-Oki earthquake.
The epicenters of earthquakes (MJMA ≧ 4.0) from
Mar. 09, 2011 to Apr. 11, 2011 (JST) are shown.
Earthquakes with magnitude larger than or equal to 7 areshown by red stars. Hypocentral data determined by
Japan Meteorological Agency are used.

the seismic waves radiated from suboceanic earthquakes.
The modeling of the ground motion induced by thisevent is a computational challenge: large memor y size andfast computing devices are required because the huge source
size of the earthquake imposes a very large domain size for thesimulation.

Finite-diference time domain method

2

In this paper we review our 3-D finite - dif ference

11

time domain (FDTD) method [1-3]. In order to simulate the wave

We apply the time domain, staggered-grid three-propagation with a large grid size, we adopt the GPU (graphics

dimensional finite-difference scheme [4].  The components of theprocessing unit) computing to our finite-difference program.

particle velocity (v i; i=x, y, z) and the stress tensor (τij; i, j=x, y, z)

We then present the results of the simulation of the wave

are the field variables. We use three material parameters (Lamépropagation based on a preliminary source model of the 2011

coef ficients and density) as we assume isotropic and elastic

Tohoku-Oki earthquake.

material for the simulation. Thus twelve variables are assigned

to a unit cell (i.e., a heterogeneous formulation: Fig. 2). The

with 1408 nodes and 4224 of NVIDIA M2050 GPUs, and has aaccuracies of finite-differences are fourth-order in space andpeak performance of 2.4 PFlops (peta-flops). It is ranked as thesecond-order in time, respectively. In order to incorporate theworld fifth fastest supercomputer in the November 2011 TOP-

effects due to heterogeneities and irregularities in the structure,

500 list (www.top500.org).

we have adopted an approach to model structures with both ofthe land and ocean topography and heterogeneity in 3D seismic
modeling with the finite-difference method [5]. The approachunifies the implementation for irregular free-surface and that for
irregular liquid-solid interface.

Figure 3  A simplified diagram of a                  GPU (NVIDIA M2050) used in this study.
                  A single GPU has 14 processors (multiprocessors)                  and a single multiprocessor has 32 processor cores.
                  Thus 448 processor cores are integrated in a GPU.
                  The size of the global memory is 3 GB.

Memory optimization
Figure 2  A single unit cell of the staggered grid.
                 Open symbols denote the variables that belong                 to the neighboring cells [2].

The GPU is characterized by its hundreds of cores forarithmetic instructions (Fig. 3): we need to provide sufficientamount of data to these cores. The bandwidth of the main
memory (called “global memory” in NVIDIA GPUs) connected

GPU computing

3

to the GPU is much faster than that of the memor y system
connected to the conventional CPUs. However, 400 to 600 clock
cycles of memory latency still occurs in transferring the data
between the global memory and the GPU. It is thus critical for
performance to optimize the memory manipulation in memory-

The GPU (graphics processing unit) is a remarkableintensive applications such as the simulation of seismic wave

d ev i ce fo r i t s many - co re archi te c t ure an d hi gh m em o r ypropagation. Therefore, we use the fast (but small) memories

bandwidth (Fig. 3). The GPU delivers extremely high computing

in the GPU, the registers and the shared memor y , as sof tware

p er fo r man ce at a re du ce d p ower an d cos t co mp are d to

managed cache memories to reuse the data and hence toconventional central processing units (CPUs). The simulation

reduce the data transfer from the global memory [1, 6-9] (see Fig. 4of seismic wave propagation is a t ypical memor y-intensive

for the method).

problem which involves a large amount of data transfer between

In order to reduce the access to the global memory, we

the memor y and the arithmetic units, while the number ofdefine the material parameters only at the center of the unit cell.

arithmetic operations is relatively small. This is the reason we

The material parameters at the grid points for particle velocities

adopt the GPU computing to the simulation of seismic wave

and shear stresses are computed by using the values definedpropagation: it can benefit from the high memory-bandwidth ofat the center of the unit cells at every time steps. The unified

the GPUs.

implementation method [5] allows us to compute the appropriate
We use the TSUBAME 2.0 supercomputer in Global

values of the material parameters. We further use the look-up

Scientific Information and Computing Center, Tokyo Institute oftable method for the material parameters. In addition to theseTechnology, for the seismic wave simulation employed in thistechniques we optimize the number of grid points (i.e., the block

paper. The TSUBAME 2.0 is a large-scale GPU cluster equipped

size ) assigned to a single execution unit (called a multiprocessor ),12

GPU-Accelerated Large-Scale
Simulation of Seismic-Wave Propagation

because the memory transfer rate is better for grouped memory

In the parallel computing we need to exchange thetransaction using blocks of proper size (typically 16 to 64) thandata in the ghost zones between the neighboring subdomains.

that for serialized, one-by-one memory transaction.

We adopt the three-dimensional (3D) domain decomposition

Fig. 5 shows the per formance of the (multi- GPU-

(Fig. 6). Here we note that, in GPU computing, the data withinc ap ab l e) FDT D p ro gr am e xe cu te d o n a sin gl e - GPU. A s a

the ghost zones in the GPU memory must be once copied tocomparison we show the per formance of a FDTD program

the host memory in order to exchange them with other GPUs

executed on the host CPUs. Note that this example for CPU does

installed on the other nodes (see Fig. 3). Thus in decomposing

not involve MPI inter-node communications: the program isthe domain, we use dedicated memory buffers for the ghost

parallelized with OpenMP and executed on a single node. For

zones to improve the data transfer rate bet ween the GPU

these programs, the performance of a single GPU is roughly

memory and the host memory [1,3].

three-fold faster than that of a single node (2 CPUs, 12 cores).

We also overlap the procedures for the communicationand the computation to reduce the total processing time. In

Parallel computing with multi-gpus

our program we first compute for the grid points at and near
the sides of the subdomains for ghost point exchange. Then

For parallel computing with multi-GPUs, we decompose thewe start the exchange procedure for the ghost points and the

FDTD domain into subdomains and allocate a single subdomain

computations for the internal grid points simultaneously.

to a single GPU. This is necessary for large-scale simulations

With these optimizations, the performance of our 3D

because the size of the global memory of a GPU is not large (e.g.,

FDTD program scales well with the increasing number of GPUs

3 GB in the case of NVIDIA M2050). We here use the MPI library(Fig. 7): by using 400 nodes of the TSUBAME supercomputer andfor our parallel FDTD program.

under the condition of 2 GPUs activated per a node, a very highsingle-precision performance of about 50 TFlops is achievedin the case of 800 GPUs (85 % of the complete scalability). Even
higher performance of about 61 TFlops (69 % ) is achieved in thecase of 1200 GPUs under the condition of 3 GPUs per a node [3].

Figure 4  
Illustrations showing the usage of the shared memory andregisters. (a) All the variables in the whole subdomain is stored in
the global memory in the GPU. We further copy the variables on achosen level into the small, two-dimensional (2D) blocks assignedin the shared memory. We assign a single thread to each unit cell.
Thick blue arrow denotes the marching direction of thecomputation. (b) Illustration of grid points for the finite-difference.
Values on red points are stored in shared memory and thoseon gray point in registers. To the variables not on the 2D plane,
only “vertical” finite-differences are operated so that they are notrequired to be shared among different threads.
When the computation proceeds to the next level, the variables inthe shared memory are copied to the registers and vice versa.
This reduces the data transfer from the global memory.

13

Figure 5
Performance of GPU and CPU in GFlops.
A single GPU is used to execute the (multi-GPU-capable)
FDTD program with a domain size of 320 × 320 × 320.
On CPUs we parallelized a Fortran source code by using
OpenMP with 1D domain decomposition, compiled it by PGIFortran compiler with "-fastsse" option, and executed theprogram on a single host node with 12 cores (2 CPUs) in total.

stations that are far from the source region and hence the model
represents long-wavelength feature of the source. The source
model consists of point sources distributed on the grid points(23 × 11) on the fault plane with a spacing of about 20 km. We
compile several detailed structure models [11-14] in and aroundthe northeastern Japan region to construct the 3D model for thesimulation. The 3D model incorporates the irregular topographyof the land and ocean area, ocean water layer, sediments, crustand subducting plates. For the example shown here we wererequired to use 1000 GPUs (334 nodes) of the TSUBAME 2.0
supercomputer. This is about 24% of the entire resources of the

Figure 6
Schematic illustration of the 3D domain decomposition.
Blue regions indicate the ghost zones.

TSUBAME 2.0. The FDTD parameters are as follows: the grid
spacing 0.15 km, the time increment 0.005 s, the grid size 6400
× 3200 × 1600, time steps 44000, and the maximum frequency
0.61 Hz.
Fig. 8 and 9 show the snapshots of the ground motion(particle velocity) computed for the preliminary source modelof the 2011 Tohoku-Oki earthquake [2]. The ground motions at
the free-surface in the land area and at the solid side of theocean bottom in the oceanic area are visualized in color scale.
Although not visualized, the propagation of the pressure waves
within the ocean water layer is also modeled in the simulation.
We are able to observe strong interference due to the complex
source model and the heterogeneous structure. The complexity
is especially significant in the later part of the simulation. It isalso remarkable that, from about 130 s to 170 s, regions withstrong positive motion (red) and negative motion (blue) emerge

Figure 7
Weak scaling curve of Multi-GPU case on TSUBAME–2.0.
The subdomain size was fixed to 320 × 320 × 320.
The total number of subdomains is equal to the number of used
GPUs. The numbers of unit cells are also shown for several points.
The experiments were performed with 2 GPUs per a node, except
in the case of 1200 GPUs that was performed with 3 GPUs per anode.

near the coast of the Fukushima prefecture (Fig. 9, yellow circle).
Motions in the regions overwhelm the wave-field propagatedfrom around the epicenter, and largely distort the pattern of
the wave propagation. The source that provides this feature is
the relatively large slip near the Fukushima prefecture retrievedin the source model [2] . Such a feature was also observed in
the strong ground motion records and the strong-motiongeneration-areas (SMGAs) were identified below the coast of
Miyagi to Fukushima prefectures [15]. These characteristics must

Simulation of Tohoku-Oki earthquake4

be studied further in the future analysis, with improvements inthe source model as well as in the 3D structure model.

The aftershock distribution of the 2011 Tohoku-Oki
earthquake indicates the approximate source area (i.e., size ofthe fault) of about 500 km long and 200 km wide. We use apreliminary source model determined by ourselves [10] for the
results shown in this paper. The preliminary source model was
determined by using waveform data recorded at world-wide

14

GPU-Accelerated Large-Scale
Simulation of Seismic-Wave Propagation

Figure 9
Same as Fig. 8 but for snapshots at 130 s, 150 s, 170 s and 190 s.

Conclusion

5

In order to simulate the wave propagation from the2011 Tohoku-Oki ear thquake, we constructed a preliminar y
3D struc ture model for the nor theastern Japan region byco mp ilin g th e m o d els fo r to p o gr ap hy, s e dim ent s, cr us t ,Figure 8
Snapshots of the simulated wave-field. The vertical componentof the particle velocity is visualized with color scale. Red and blue
colors denote upward and downward motions, respectively.
From left to right, top to bottom the snapshot at 10 s after the
onset of the rupture, at 30 s, 50 s, 70 s, 90 s, and 110 s, respectively,
are shown.

15

and subducting plates. We adopted the GPU computing toaccelerate the large-scale finite-dif ference simulation of theseismic wave propagation. By using about 24% of the resourcesof the GPU supercomputer, TSUBAME 2.0, at Tokyo Institute ofTechnology, Japan, we are now able to simulate the seismic
wave propagation from the source of the huge fault plane up to0.61 Hz within a tolerable computation time. The snapshots of

Performance Computing & Simulation, pp. 36-43 (2009).

simulated wave-field exhibit strongly complex pattern because[7] Micikevicius, P. “3D finite-difference computation on GPUsof the complex source model, the irregular topography andusing CUDA”, GPGPU-2: Proc. 2nd Workshop on General

the shallow heterogeneities. Such effect must be considered

Purpose Processing on Graphics Processing Units, pp. 79-

in improving the source model, the 3D structure model, and in

84, Washington DC, USA (2009).

the future quantitative study on the observed ground motion[8]

Michéa, D. and Komatit sch, D. “Acceler ating a thre e dimensional finite-difference wave propagation code using

records.

GPU graphics cards”, Geophys. J. Int., doi: 10.1111/j.1365246X.2010.04616.x (2010).

Acknowledgements[9]

Okamoto, T., Takenaka, H. and Nakamura, T. “Simulation

We are grateful to the researchers who provided the structuralo f s e is m i c w av e p r o p a g a t i o n b y G PU ”, P r o ce e d i n g smodels. We are also grateful to Global Scientific Informationof s y mp osium on a d v an ce d computin g s ys tems andand Computing Center (GSIC), Tokyo Institute of Technology, forinfrastructures, pp. 141-142 (2010).

providing opportunity (Grand Challenge) to use large number[10] O k amoto, T., Takenak a, H., Har a, T., Nak amur a, T. andof nodes of TSUBAME. This research was partially supported by

Aoki, T. “Rupture Process And Waveform Modeling of The

Grant-in-Aid for Scientific Research (KAKENHI 23310122) from2011 Tohoku- Oki, Magnitude -9 Ear thquake”, American

Japan Society for the Promotion of Science (JSPS).

Geophysical Union, Fall Meeting, U51B-0038, San Francisco,
USA (2011).

References[11] Kisimoto, K. “Combined bathymetric and topographic mesh
data: Japan250m.grd”, Geological Survey of Japan, Open-

[1]

[2]

O k am oto, T., Takenak a, H . , Nak amur a, T. an d Ao k i, T.

[12] Fujiwara, H., Kawai, S., Aoi, S., Morikawa, N., Senna, S.,

propagation by multi-GPUs and three-dimensional domain

Kudo, N., Ooi, M., Hao, K. X.-S., Hayakawa, Y., Toyama, N.,

decomposition”, Earth, Planets and Space, vol. 62, no. 12, pp.

Matsuyama, H., Iwamoto, K., Suzuki, H. and Liu, Y. “A study

939-942 (2010).

on subsurface structure model for deep sedimentary layers

O k amoto, T., Takenak a, H., Nak amur a, T., and Aok i, T.

of Japan for strong-motion evaluation”, Technical Note

"Large-scale simulation of seismic-wave propagation ofof the National Research Institute for Earth Science andthe 2011 Tohoku-Oki M9 earthquake", Proceedings of the

[3]

[4]

file Report, No. 353 (1999).

“Acceler ating large -sc ale simulation of seismic wave

Disaster Prevention, No.337 (2009).

International Symposium on Engineering Lessons Learned[13] Baba, T., Ito, A., Kaneda, Y., Hayakawa, T. and Furumura, T. “3-Dfrom the 2011 Great East Japan Earthquake, pp. 349-360

seismic wave velocity structures in the Nankai and Japan(2012).

Trench subduc tion zones derived from marine seismic

Okamoto, T., Takenaka, H., Nakamura, T. and Aoki, T. “GPU-surveys”, Abstr. Japan Geoscience Union Meet., S111-006,

accelerated simulation of seismic wave propagation”, in

Makuhari, Japan (2006).

GPU Solutions to Multi-scale Problems in Science and[14] Nakamura, T., Okamoto, T., Sugioka, H., Ishihara, Y., Ito, A.,

Engineering, Yuen, D., Wang, J., Johnsson, L., Chi, C.-H., Shi,

Obana, K., Kodaira, S., Suetsugu, D., Kinoshita, M., Fukao, Y.

Y. (Eds.), 250 pp., Springer, in press.

and Kaneda, Y. “3D FDM simulation for very-low-frequency

Graves, R. W. “Simulating seismic wave propagation in 3D

earthquakes off Kii Peninsula”, Abstr. Seism. Soc. Japan, P1-

elastic media using staggered-grid finite differences”, Bull.

06, Hiroshima, Japan (2010).

Seism. Soc. Am., vol. 86, pp. 1091-1106 (1996).

[15] Kurahashi, S. and Irikura, K. “Source model for generating[5] Takenaka, H., Nakamura, T., Okamoto, T. and Kaneda, Y. “A

strong ground motions during the 2011 off the Pacific coast

unified approach implementing land and ocean-bottomof Tohoku earthquake”, Earth Planets Space, Vol. 63, 571-576

top o graphies in the s taggere d - grid f inite - dif ference(2011).

method for seismic wave modeling”, Proc. 9th SEGJ Int.
Symp., CD-ROM Paper No.37 (2009).
[6]

Abdelkhalek, R., Calandra, H., Coulaud, O., Roman, J. andLatu, G. “Fast seismic modeling and reverse time migration
on a GPU clus ter ”, Inter national Conference on High16

● TSUBAME e-Science Journal No.6
Published 07/31/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/

17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