15

TTX: A Direct Numerical Simulation Codefor Turbulent Reacting Flows
A High-level Framework for Parallel andEfficient AMR on GPUs
DS-CUDA: A Handy Tool to Use GPUs
in a Cloud Network

TTX : A Direct Numerical Simulation Codefor Turbulent Reacting Flows
Yuki Minamoto* Kozo Aoki** Mamoru Tanahashi** Department of Mechanical Engineering, Tokyo Institute of Technology**Mechanical and aerospace engineering, Tokyo Institute of Technology

Over eighty percent of world energy supply is estimated to be provided by burning fossil fuels for the next
30 years, and this imposes more and more stringent environmental regulations on the design of combustion
related devices. Direct numerical simulations (DNS) have played important roles in the research of turbulent
combustion. DNS data base provide key information for the development of turbulent combustion models,
which are to be used in computational fluid dynamics (CFD) of various combustion devices during design anddevelopment phases. This article summarises a DNS code and implemented numerical techniques that aredeveloped at Tokyo Institute of Technology.

Introduction

1

computing elementary reaction rates existing in chemicalmechanisms takes most of combustion DNS costs and amethod to improve this computation speed is necessary
to achieve such DNS. The present research implements
such a boundary technique and methods to accelerate

Over eighty percent of world energy supply is estimated tokinetic computation in a DNS code, TTX, and simulations are

be provided by burning fossil fuels for the next 30 years [1],performed to confirm its accuracy and performance for several

and this imposes more and more stringent environmental

laminar/turbulent reacting/non-reacting flow problems.

regulations on the design of combustion related devices.
Direct numerical simulations (DNS) have played importantroles in the research of turbulent combustion. DNS data baseprovide key information for the development of turbulent
combustion models, which are to be used in computational

Mathematical Backgroundfluid dynamics (CFD) of various combustion devices during

2

design and development phases.
For decades, DNS have been used for canonicalcombustion problems such as statistically 1D planar
propagating premixed flames and HCCI-type combustion

The governing equations consist of fully compressible

in a periodic rectangular domain. However, with the recent

conservation equations for mass, momentum, energy andadvancement in high-performance computing, DNS of slightly-

mass fractions of N -1 chemical species, where N is the number

more complicated and computationally costly combustionof chemical species involved in combustion. The equations

configurations such as V-flame, jet flame and flames stabilized

are discretised by a fourth order central finite differencein jet-in-cross-flow and swirl flow with using complex chemical

scheme and integrated in time by using a third ordermechanisms have been performed, and such simulations

explicit Runge-Kutta scheme on Cartesian mesh. Unphysicalwill further our understanding on the physics of turbulent

numerical oscillations resulted from the use of finite difference

combustion. Since these configurations include walls that do

scheme are removed by using either a sixth order explicit

not necessarily conform with the preferred structured meshspatial filtering or a fourth order compact filter. All open

coordinates for combustion DNS, most of these simulations

computational boundaries are described

use presumed profiles for inflow/near-wall flows as boundary

based on the Navier-Stokes characteristic boundary conditions

conditions. A portable high-order immersed boundary method(NSCBC) formulation. The detailed description can be found insuited for parallel computation is one way to improve theseprevious studies [2-6].

simulations. Also, the use of more practical and complicatedhydrocarbon fuels such as methane, n-heptane and gasoline

02

2.1 Governing equations and PDE solver

2.2 MTS and CODACT

surrogates with complex chemical kinetic models is important

One of the difficulties in conducting combustion DNS is theto reveal underlying flame — turbulence interaction that istime integration of the governing equations with chemicala key for the development of “green” combustor. However,reaction terms. The time stepping (Δt ) for fully explicit＊＊

Current position : Japan Atomic Energy Agency

schemes has to be smaller than the smallest chemicaltimescale which can become much smaller than the time

2.3 Immersed boundary method

stepping limited by a CFL condition. MTS (multi-timescale)

The present immersed boundary method (IBM) is based on

method [7] is an algorithm for integrating ODEs associated withwell-known ghost region/reverse profile approach [11] withchemical reaction using a large Δt . In MTS method, eachsome modification to achieve high accuracy and portability.

chemical species ( Y k ) belongs to one integration group based

The important change from the conventional ghost region

on its chemical timescale (τk ). The index (M ) of an integration

approach is that the ghost points are identified not only fromgroup which Y k belongs to is obtained as

the non-fluid points that face to the fluid region, but the

M = log10 ( Δt /τk ) + 1. Species in M th group are integrated

region has certain non-fluid mesh points in the wall normal

using a time stepping Δt M = Δt /10 ( -1). The integration

directions. The points width in wall normal direction is max ( n d ,Mprocedure starts with the smallest Δt M . After convergence ofn f /2 ) for uniform one-dimensional configuration, where

M th group, the mass fractions of the species in the group aren d and n f are the order of differentiation and explicit filter,

fixed and then the ODE system continues to be integrated

respectively. Figure 1 shows an example schematic of ghost

using Δt M -1. Thus, MTS enables us to use a large Δt , which isregion for a n d = 4 and n f = 6 case. Once the ghost points are

not limited by the smallest chemical timescale, as a global

identified, a reference point is identified for each ghost point:time stepping for DNS.
Another approach to reducing computational costis reduction of the reaction ODE system. CODAC (correlated
dynamic adaptive chemistry)[8]

is one of the methods forwhere the threevectors denote the locations of reference

point, corresponding wall and ghost point, respectively. For

dynamically generating reduced kinetic mechanisms at each

each physical quantity,

spatial location and each time step. In CODAC method,equations, the reference value is obtained by interpolating thea phase space consisting of a few key parameters isfield onto the reference point. Using the reference point, theconstructed. Temperature and mass fractions of fuel, oxidizer

value at each ghost point is estimated depending on the walland important radicals are selected as the key parameters.

conditions as:, solved in the set of governing

Comparing the instantaneous phase parameters at a grid
point with those at another reference point (different point intime/space), correlation between the two points in the phasespace is examined. Threshold values used for examining thecorrelation are specified by a user. If the phase parameters at

Here, Eq. (2) is used for Dirichlet, and Eq. (3) is for Neumann

the two points are correlated, the reduced mechanism used atwall boundary conditions. In the present implementation,

the reference point is reused at the correlated grid point. If theghost values for velocity and temperature are obtained byphase parameters at a grid point are not correlated with any

using Eq. (2), while species mass fractions and pressure areother reference point, a new reduced mechanism needs to betreated with Eq. (3). The density values on ghost points are

generated by PFA (path flux analysis) [9].

obtained using a combination of Eqs. (2) and (3) to achieve

Further reduction of computational cost is achievedby using CODACT (CODAC and transport)[10]

mass conservation.

. The same

correlation procedure in CODAC is applied to reducing thecomputation oftransport properties, i.e. viscosity, diffusivity and thermal
conductivity. The phase parameters for reducing transport
calculation are temperature and molar fractions of the first
several abundant species. As with CODAC, transport properties
at a grid point are not calculated if the phase parameters at
the point are correlated with those at a different grid point.

Fig. 1 Schematic of DNS domain. Blue shadeindicates identified ghost region.

03

TTX : A Direct Numerical Simulation Code for Turbulent Reacting Flows

Performance and Accuracy

3

the use of CODACT does not influence the results. Figure 4
shows accumulated computational times for the cases (i) —(iii). Clearly, due to significantly small Δt , the computationaltime for the case (i) takes longest of all. Clearly, the use of MTS
and CODACT decreases computational by a factor of 7— 8

3.1 MTS and CODACT

times for complex chemistry. Similar trends are observed formore complex fuel combustion mechanism such as n-heptane.

A typical computational performance is shown in Fig. 2 interms of weak scaling. The used DNS configuration is a periodic
cuboid domain without using MTS, CODACT and IBM. The
domain was initially filled with a reactant mixture of CH 4 - air
preheated to 1000 K at 1 atm. The chemical reactions aredescribed by using GRI - 3.0 mechanism. The mixture is thenignited at the centre location by gradually imposing an energy
source term having a profile based on a Gaussian distribution.
Each CPU core are allotted 12 3 mesh points and simulations
were continued for the first 100 time steps. Clearly, TTX has
excellent parallel computation performance for a range of
number of cores.

Fig. 3 Comparison of temperature profiles
between the cases (i)—(iii).
Red: case (i),black dashed: case (ii),black solid: case (iii).

Fig. 2 Weak scaling plot of combustion DNS
in a typical configuration using TTX.
A comparison between explicit time integration
without CODACT (i), time integration using MTS without
CODACT (ii) and time integration using MTS and CODACT (iii)is shown in Fig. 3 for a homogeneous ignition simulation of
methane-air combustion. The time stepping for these cases

Fig. 4 Comparison of accumulated computationalcosts between the cases (i)—(iii).
Red: case (i), black dashed: case (ii),black solid: case (iii).
black dashed: case (ii), black solid: case (iii).

are set to be Δt =1ns (i), 5ns (ii) and 5ns (iii). Smaller Δt forthe case (i) is necessary as chemical time scales typicallydecreases rapidly with pressure rise. In Fig. 3, the temperature04

3.2 Immersed boundary method

profiles show identical variation except for the explicit caseTo access developed IBM, non-reacting laminar Taylor-Couette(i) which starts to show a deviation after t =0.5ms. This is due

flows are simulated. The inner and outer walls are describedto the lack of temporal resolution ( Δt ) since the chemical

based on the present IBM, and the DNS solutions at different

time scales are typically shortened as pressure increases withspatial resolutions are compared with the analytical solution.

ignition proceeds. Thus, without the aid of MTS, Δt needs to

The error convergence plot is shown in Fig. 5 in terms of

be set significantly smaller values than other cases. However,maximum and L 2 norm of relative errors. The convergencerate of L 2 norm decreases at a rate of between third andevolve from a current relatively small, canonical configuration

fourth orders, and this is comparable with the accuracy of

to a larger scale, more complex geometry that is more relevant

differentiation in the interior points in the present DNS.

to practical combustion devices.

Fig. 6 Configuration and streamwise velocity fieldof turbulent pipe flow DNS.
Fig. 5 Error convergence plot of the present IBM
obtained from a laminar Taylor-Couette DNS.
DNS of a turbulent pipe flow is performed. The
present pipe wall boundary conditions are specified as in thedeveloped IBM, whereas previous pipe flow DNS databaseapplied structured mesh that represent the pipe shape. The
initial nominal Reynolds number Reτ is set as 200, which isdefined based on the friction velocity and radius of the pipe.
The dimensions of the pipe are shown in Fig. 6 with the
DNS result for the streamwise velocity field. The simulation
was continued for over 40 mean-flow-through times toachieve quasi-steady state and samples are collected after 25
mean-flow-through times for statistics. Turbulent statistics
are compared with previous DNS with conventional wallFig. 7 Comparison of the mean streamwise velocity
profile of the pipe flow DNS.

boundary method for mean streamwise velocity and Reynolds
shear stress (not shown) in Fig. 7. The present pipe flow DNS
results show excellent agreement with the previous DNS
results with Reτ = 180 and 181, ensuring that the present IBM
can capture the large velocity gradient and reproduce theturbulent boundary layer accurately. Also, the computational
cost per time step is 1020 ms with IBM and when IBM is
switched off, it decreases to 950 ms with 1296 cores at
Reedbush-U system at the University of Tokyo. Therefore, theadditional computational cost due to the IBM is less than 10%of the overall cost.
3.3 3D Turbulent combustion simulations
Typical DNS results of three-dimensional turbulent combustion
performed by using TTX is shown in Fig. 8. With the numerical
techniques described here, the turbulent combustion DNS will

Fig. 8 Typical combustion DNS results of swirl
flow stabilised turbulent flames of hydrogenair combustion. Iso-surfaces: heat release rate.
Color: normalised temperature. Left: low swirl
number (0.6) case.
Right: high swirl number (1.2) case05

TTX : A Direct Numerical Simulation Code for Turbulent Reacting Flows

Summary

timescale and correlated dynamic adaptive chemistry4

A direct numerical simulation code, TTX, is developed bymodeling of ignition and flame propagation using a real
jet fuel surrogate model, Combust. Flame, Vol. 162, pp.
1530-1539 (2015)[9]

W. Sun, Z. Chen, X. Gou and Y. Ju: A path flux analysis

implementing MTS, CODACT and IBM capability for bettermethod for the reduction of detailed chemical kinetic

computational speed and relatively flexible boundary

mechanisms, Combust. Flame, Vol. 157, pp. 1298-1307

configurations. The computational speed is increased by(2010)

a factor of 7 to 8 when DNS are performed with both MTS[10] W. Sun and Y. Ju: Multi-timescale and Correlated Dynamic

and CODACT for a complex chemical mechanism, while theAdaptive Chemistry and Transport Modeling of Flames

implemented IBM shows comparable accuracy to the accuracy

in n-Heptane/Air Mixtures, Proc. 53rd AIAA Aerospaceof present numerical schemes.

Sciences Meeting (2015).
[11] P. Parnaudeau, E. Lamballais, D. Heitz, J. H. Silverstrini:

Acknowledgementswith Compact Schemes for DNS of Flows in Complex

This research was supported in part by Grant-in-Aid for Young

Geometry, ERCOFTAC seriese, Direct and Large-Eddy

Scientists B (16K18026) of the Japan Science and Technology

Simulation V, Vol. 9, pp. 581-590.

Agency (JST).

References[1]

U.S. Energy Information Administration: Internationalenergy outlook 2010, DOE/EIA-0484 (2010)[2]

M. Tanahashi, M. Fujimura, and T. Miyauchi: Coherent
fine-scale eddies in turbulent premixed flames, Proc.
Combust. Inst., Vol. 28, pp. 529–535 (2000)[3]

Y. Minamoto, N. Fukushima, M. Tanahashi, T. Miyauchi, T.
D. Dunstan, and N. Swaminathan: Effect of flow-geometry
on turbulence-scalar interaction in premixed flames.
Phys. Fluids, Vol. 23, 125107 (2011)[4]

M. Shimura, K. Yamawaki, N. Fukushima, Y. S. Shim, Y.
Nada, M. Tanahashi, and T. Miyauchi: Flame and eddy
structures in hydrogenair turbulent jet premixed flame, J.
Turbulence, Vol. 13, pp. 1–17 (2010)[5]

B. Yenerdag, N. Fukushima, M. Shimura, M. Tanahashi,
and T. Miyauchi: Turbulence-flame interaction and fractal
characteristics of H2-air premixed flame under pressure
rising condition. Proc. Combust. Inst., Vol. 35, pp.1277–
1285 (2015)[6]

Sussman, P. Smereka and S. Osher: A Level Set Approachfor Computing Solutions to Incompressible Two-Phase
Flow, J. Comp. Phys., Vol. 114, pp.146-159 (1994)[7]

X. Gou, W. Sun, Z. Chen and Y. Ju: A dynamic multitimescale method for combustion modeling withdetailed and reduced chemical kinetic mechanisms,
Combust. Flame, Vol. 157, pp. 1111-1121 (2010)

[8]

06

Combination of the Immersed Boundary Method

W. Sun, X. Gou, H.A. El-Asrag, Z. Chen and Y. Ju: Multi-

A High-level Framework for
Parallel and Efficient AMR on GPUs
Mohamed Wahib* Naoya Maruyama* Takayuki Aoki*** RIKEN Advanced Institute for Computational Science**Global Scientific Information and Computing Center, Tokyo Institute of Technology

Adaptive Mesh Refinement methods reduce computational requirements of problems by increasing resolution
for only areas of interest. However, in practice, efficient AMR implementations are difficult considering thatthe mesh hierarchy management must be optimized for the underlying hardware. Architecture complexity ofGPUs can render efficient AMR to be particularity challenging in GPU-accelerated supercomputers.
This article presents a compiler-based high-level framework that can automatically transform serial uniform
mesh code annotated by the user into parallel adaptive mesh code optimized for GPU-accelerated
supercomputers. Experimental results on three applications show that the speedups of code generated by
our framework are comparable to hand-written AMR code, while achieving good and weak scaling up to 3,640
GPUs of the TSUBAME2.5 supercomputer.

1

Introduction

own versions of the target-optimized solvers. Moreover, there
can be scalability limitations caused by the overhead of the
CPU-GPU communication schemes in those frameworks[2].
In this article, we present a high-level framework,
called Daino[3], which auto-generates efficient and scalable

In many scientific and engineering simulations, Partial

structured AMR solutions to scientific applications running on

Differential Equations (PDEs) are solved in a uniform meshGPU-accelerated systems.

arrangement by using finite difference schemes, referred toas iterative stencils. Typically, the resolution of the mesh is
uniformly set to the highest resolution to provide accurate
solutions. For meshes that require only high resolution forsome portions of the mesh, an alternative method, known

Backgroundas Adaptive Mesh Refinement (AMR), can be used instead of2

the uniform mesh. The AMR method solves the problem ona relatively coarse grid, and dynamically refines it in regions

Structured AMR methods use logically rectangular meshes inrequiring higher resolution. However, AMR codes tend to be

the implementation of the adaptive mesh. Structured AMR

far more complicated than their uniform mesh counterpartsutilizes a hierarchy of levels of spatial, and often temporal,

due to the software infrastructure necessary to dynamicallymesh spacing with each level being composed of a union

manage the hierarchical mesh framework. Despite thisof logically rectangular mesh regions. One way to represent

complexity, it is generally believed that future applications will

the structured AMR, namely tree-based AMR, divides theincreasingly rely on adaptive methods to study problems at

discretized domain into fixed blocks. If any cell within a block

unprecedented scale and resolution.

requires refinement, the whole block is refined.

Implementing efficient adaptive meshes in GPU-In the tree-based scheme, the mesh is organizedaccelerated systems is significantly hard in comparison

into a hierarchy of refinement levels. The mesh is usuallyto traditional CPU systems. More specifically, GPUs add

decomposed into relatively small fixed-sized octants ofcomplexity overhead for managing the mesh hierarchy andmesh cells. Each octant can be recursively refined into a set

optimization of data movement. This is made evident by theof octants of fine cells. The mesh configuration is managed

relatively wide use of AMR in CPU in comparison to GPU-basedusing a tree-based data structure that maintains explicit child-systems. For example, a mature AMR framework supporting

parent relationships between coarse and fine octants. SizeCPU, namely FLASH is reported to be in use by dozens ofrelations between neighboring octants are typically enforced[1]

production applications . On the contrary, a few number ofin structured AMR, which means neighboring octants can

individual applications adapted AMR solvers for GPU, withhave at most one level of refinement difference (referred tovarying levels of optimization and scaling. To summarize

as 2:1 balance). An important feature of octrees is that thethe problems with GPU-based AMR, only a few frameworks

traversal of an octree across its leaves corresponds to a Morton

enable automated AMR transformations for GPU, and theirz-shaped Space Filling Curve (SFC) in the geometric domain[4].

programming models require the programmer to write his

Accordingly, sorting the blocks by their Morton ID and equally

07

A High-level Framework for Parallel and Efficient AMR on GPUs

partitioning them leads to a uniform distribution of the blocks
of a mesh among different Processing Elements (PEs), while
benefiting from the locality provided by SFC affinity. Fig.1
illustrates how the domain and tree are represented in AMR,
and the use of SFC to divide the blocks among three PEs.

Fig. 1 Octree-based meshes. The blocks are equallypartitioned using a space filling curve.
(a) Adaptive mesh(b) Tree representation
		 (Note: 2D quadtree is used for illustration)
Fig. 2 Minimal example of using Daino directives

High-level Framework for
Efficient AMR

3

3.2 Optimizations
When an AMR code generated by Daino is executed on a GPUaccelerated cluster, the stencil and mesh adaptation kernels

08

We design a high-level programming framework that provides

run on the GPU, while managing the octree data structures

a highly productive programming environment for AMR. The

and load balancing is done on the CPU side. Since we pursue

framework is transparent and requires minimal involvement

efficiency and scalability, code on both the CPU and GPUfrom the programmer, while generating efficient and scalable

should be optimized. The stencil operations in the blocks

AMR code. The framework consists of a compiler and runtime

are themselves optimized for the GPU architecture[5]. Other

components. A set of directives allows the programmer to

optimizations, such as data layout in memory and using user-

identify stencils of a uniform mesh in an architecture- neutral

managed cache memory, are applied on the GPU kernels

way. The uniform mesh code is then translated to GPU-

responsible for adapting the mesh: error estimation, refinement

optimized parallel AMR code, which is then compiled to an(interpolation), and coarsening (extrapolation). Finally,

executable. The runtime component encapsulates the AMR

the generated code includes optimizations to reduce thehierarchy and provides an interface for the mesh management

communication between nodes, i.e. boundary data exchange,

operations.

and balancing the load, i.e. number of blocks per node.

3.1 Programming Model

3.3 Implementation

The framework provides directives to be used with standard

Our framework consists of a compiler and runtime

C (see [3] for details on directives). The programmer is requiredcomponents. We generate executables optimized for GPU

to add the directives to a serial uniform mesh code in order

execution by leveraging the LLVM compiler infrastructure.

to identify the operations and stencil data arrays that are the

The compiler builds on the LLVM compiler infrastructure[6].

target for transformation. Note that the directives are notFirst, we use the front end to analyze and translate the stencil

changing the uniform mesh implementation; the programmer

source code into GPU-optimized code in the form of LLVM

can still use the uniform mesh implementation if the directives

Intermediate Representation (IR). Next, compiler passes areare ignored by compiler. A sample example of using directives

applied on the IR to add the AMR management code, which

to annotate a C kernel in Daino is shown in Fig. 2.

in turn make API calls to the runtime API and GPU-optimized

code generated by the Nvidia backend code generator. Finally,

First, more than 1.7x speedup is achieved using Daino using

LLVM IR is compiled and linked with the runtime libraries to

up to 3640 GPUs of TSUBAME for the phase-field simulation.

generate and executable.

This is a considerable improvement considering that the

The runtime includes two libraries: the first library

uniform mesh implementation is a Gordon Bell prize winnerencapsulates the AMR hierarchy management software andfor time-to-solution[7]. Second, Daino achieves good scaling

the second is a communication library that wraps the MPI

that comparable to the scalability of the hand-written AMR

runtime library to simplify data movement operations for thecode.

AMR driver. The stages of compilation and layout of Diano areshown in Fig. 3.

Fig. 3 The framework overview

Evaluation

4

We demonstrate the scalability of auto-generated AMR codeusing three production applications. We compare the speedup
and scalability with hand-written AMR of all three applications
using the TSUBAME2.5 supercomputer at Tokyo Institute ofTechnology.
4.1 Applications
Phase-field Simulation : This application simulates 3D dendritic

growth during binary ally solidification[7].
Hydrodynamics Solver : This solver models a 2nd orderdirectionally split hyperbolic schemes to solve Euler
equations[8].
Shallow-water : Modelling shallow water by depth-averaging

Navier-Stokes equations[9].
4.1 Results
In a weak scaling experiment, shown in Fig. 4, the run- time foruniform mesh, hand-written AMR, and auto-generated AMR
are compared. The following points are important to note.

Fig. 4 Weak scaling of uniform mesh,hand-written and automated AMR

09

A High-level Framework for Parallel and Efficient AMR on GPUs

Fig. 5 shows a strong scaling comparison for handwritten and auto-generated AMR against uniform mesh
implementation. The code generated by Daino achieves

Summary

speedups and scalability comparable to hand-written

5

implementations. However, when using more GPUs, reductionin speedup starts to occur as the management of AMR starts

Producing efficient AMR code is a challenge, especially forto dominate the simulation runtime.

GPUs. We introduce a framework for producing efficientand distributed AMR code for GPU-accelerated systems. We
demonstrated the efficacy and scalability of three applications
using the full TSUBAME supercomputer. To the authors
knowledge, this is the first study to scale auto-generated AMR
code to O(1,000s) of GPUs. However, there are still problems tobe solved, for example, handling of customized error functions
and boundary conditions.

Acknowledgements
This project is partially supported by JST, CREST through its
research program: "Highly Productive, High Performance
Application Frameworks for Post Petascale Computing.".
This research is also partly supported by KAKENHI, Grant-inAid for Scientific Research (S) 26220002 from the Ministry of
Education, Culture, Sports, Science and Technology (MEXT) of
Japan and "Joint Usage/Research Center for Interdisciplinary
Large-scale Information Infrastructures (JHPCN)" and "High
Performance Computing Infrastructure (HPCI)" in Japan.

References[1]

A. Dubey, K. Antypas, M. K. Ganapathy, L. B. Reid, K.
Riley, D. Sheeler, A. Siegel, and K. Weide, “Extensible
Component-based Architecture for FLASH, a Massively
Parallel, Multiphysics Simulation Code,” Parallel Comp.,
vol. 35, no. 10-11, pp. 512–522, (2009)[2]

M. Wahib, N. Maruyama, Data-centric GPU-based
Adaptive Mesh Refinement, IA^3 2015, Workshop on
Irregular Applications: Architectures and Algorithms, colocated with SC’15, Austin, TX (2015)[3]

M. Wahib, N. Maruyama, T. Aoki, Daino: A High-level
Framework for Parallel and Efficient AMR on GPUs, SC'16,
ACM/IEEE Proceedings of the International Conference
for High Performance Computing, Networking, Storage
and Analysis (2016)[4]

H. Tropf and H. Herzog, “Multidimensional Range Searchin Dynamically Balanced Trees,” Angewandte Informatik,

Fig. 5 Strong scaling of uniform mesh,hand-written and automated AMR

101981 (1981)[5]

N. Maruyama and T. Aoki, Optimizing Stencil

Computations for NVIDIA Kepler GPUs, 1st International
Workshop on High-Performance Stencil Computations
HiStencils'14 (2014)[6]

http://www.llvm.org[7]

T. Shimokawabe, T. Aoki, T. Takaki, T. Endo, A. Yamanaka,
N. Maruyama, A. Nukada, and S. Matsuoka, “Peta-scale
Phase-field Simulation for Dendritic Solidification on the
TSUBAME 2.0 Super- computer,” ser. SC ’11 (2011)[8]

H.-Y. Schive, U.-H. Zhang, and T. Chiueh, “DirectionallyUnsplit Hydrodynamic Schemes with Hybrid MPI/
OpenMP/GPU Parallelization in AMR,” Int. J. High Perform.
Comput. Appl., vol. 26, no. 4, pp. 367–377 (2012)[9]

M. Sætra, A. Brodtkorb, K. Lie, “Efficient GPUimplementation of Adaptive Mesh Refinement for the
Shallow-Water Equations,” J. Sci. Comput., vol. 63, no. 1,
pp. 23–48 (2015)

11

DS-CUDA : A Handy Tool to Use GPUs
in a Cloud Network
Tetsu Narumi** Faculty of Informatics and Engineering, University of Electro-Communications

Cloud services have recently added GPU support for gaming and HPC purposes.
Our DS-CUDA offloads only the calculation on GPUs to a cloud system using CUDA.
The merit of DS-CUDA is that users can freely modify their applications on the client side,
as well as supporting tablet sensors and fault tolerance of networks.
We report four topics on DS-CUDA: 1) Overview of DS-CUDA, 2) 880 fold acceleration of replica-exchange
molecular dynamics simulation using 1,024 GPUs on the TSUBAME 2.5 supercomputer,3) Overhead of fault tolerant function using overclocked GPUs,
and 4) Energy efficiency using a tablet and GPU notebook.

Introduction

1

Role sharing of a client device
and cloud servers

2

Several services and tools have been released to reduce theDoing everything on one computer often results in higher cost

cost of preparing high-performance GPUs (Graphics Processing

or lower energy-efficiency compared to role sharing of client-

Units) on a client device. Amazon EC2 [1] provides a cloudand cloud-side computers. Fig. 1 shows three ways of role

server with 16 GPUs for GPGPU (General-Purpose computing

sharing between a client device and cloud servers: A) Most of

on GPUs) applications. Nvidia GRID [2] enables us to run high-the calculation and rendering are performed on a cloud server,

end GUI applications on a low-performance client device byB) Only rendering is performed on a cloud server, and C) Only

utilizing GPU-enabled calculation and rendering on the cloud

CUDA calculation is performed on a cloud server.

server and transferring the result via a low-latency network.
In this report, we consider an interactive molecular
dynamics simulation [3], which calculates on GPUs with CUDA
(Compute Unified Device Architecture) and visualizes atthe same time. Seamless methods to use high-end GPUs
on a cloud device to achieve interactive simulations andvisualizations would be required from mobile devices liketablets which supports many sensors.
DS-CUDA [4, 5] is a tool to virtualize cloud-side GPUs as
if they are attached to a client side. It supports fault toleranceof GPUs and networks. In this report, the overview of DS-CUDA
and their recent results are presented.

Fig. 1 Role sharing of a client device and cloud servers

In the case of A), the client works as a zero-clientcomputer, which means that the keyboard input or touch
information of the client device is transferred to the cloudserver, and the rendered video or sound is transferred fromthe cloud server. However, cloud service like Nvidia GRID does12

not support the recompiling and execution of user’s own

calculation errors or a network is unstable. The reliability of

applications on the cloud server. They only provide popular

GPUs are usually lower than those of CPUs even when they are

ISV applications. Amazon EC2 supports the compilation andspecially designed for GPGPU [10]. The fault tolerance should beexecution of any application since users can customize theuseful for reliable simulations, especially for unstable network

development environment on the virtual machine, while theenvironments in cloud computing.

video transfer via a network might become the bottleneck.

In the method of A), implementing fault tolerance

Special mechanism will be needed to support tablet sensors

is a little difficult. Users can customize their applications tolike accelerometers on a client device for both services.

support fault tolerance if Amazon EC2 is used. However, theIn the case of B), a tool can hook rendering APIs
like OpenGL on a client device and transfer them to a cloud[6]

server . However, CUDA is not supported, which is a big

cost of such modification is usually very high. In the method ofC), users have potentially no cost for fault tolerance since thetool can virtually provide reliable CUDA APIs to them.

disadvantage for HPC applications.
In the case of C), rendering and light processing
including sensors are performed on a client device. Users cancompile their applications on a familiar environment withoutchanging their source codes. DS-CUDA or rCUDA [7] is suitable
for such situations. Users can enjoy the benefit of highend GPUs on a cloud server by hooking CUDA APIs in theirapplications.

Comparison against other tools

3

Table 1 Comparison of tools for using cloud GPUs

Table 1 compares the cloud tools mentioned in the previous
section. DS-CUDA and rCUDA in method C) can easily support
specific functions of applications since most parts run on theclient side. Only the information called from CUDA APIs are

Overview of DS-CUDA

transferred via a network, and attackers cannot see readable

4

information like files on the cloud server, which should
increase security. However, they currently do not support

Fig. 2 and 3 show the procedure of compilation and the blockinteroperability functions of OpenGL and Direct3D, which isdiagram of DS-CUDA system, respectively.

necessary to render the result of CUDA calculation without
moving data from the GPU memory.
The recompilation of CUDA code is not needed in
the case of rCUDA since the change of the path of dynamic
link library is enough, while DS-CUDA needs recompilation.
rCUDA requires a CUDA environment on the client side,
while DS-CUDA supports client devices without CUDA when
specially prepared [8, 9]. Special environment like CUDA is
not needed on a client side for Nvidia GRID or Amazon EC2,
but users have to prepare the same software-development
environment on the cloud server.
The merit of DS-CUDA is that it supports fault

Fig. 2 Compilation with DS-CUDA compiler

tolerant mechanism for the cases when GPUs caused

13DS- CUDA : A Handy Tool to Use GPUs in a Cloud Network

Using many GPUs
from a single thread

5

CUDA provides APIs to use multiple GPUs in a node even froma single thread program. However, parallel programming like
MPI is needed to use many GPUs in a cluster by controllingthem node by node. Beginners in GPGPU might feel difficulty
in programming since a distributed memory paradigm must
be considered. Programming of GPUs with DS-CUDA is simple
since the client sees all the GPUs in the server side as if they
are attached to the client. In the following, result of using
1,024 GPUs from a single thread program is explained as anextreme case [11].
Fig. 3 Block diagram of a DS-CUDA system

In this case, replica-exchange Molecular Dynamics(MD) simulation with 14,336 replicas was executed with 1,024

First, a CUDA source code is compiled with DS-

GPUs in the TSUBAME 2.5 supercomputer. Replica-exchange

CUDA compiler, which is actually a Ruby script, and two

MD simulation can search energy minimum state efficiently

executables are generated for both the client and cloudby exchanging temperature between parallel MD simulations(or server) sides (Fig. 2). On the client side, CUDA APIs arewith different temperatures every constant steps (100 in thisreplaced with DS-CUDA APIs, and compiled with the DS-CUDA

result). One replica holds 256 or 2,048 Argon particles. Since

library to generate the executable. On the server side, thethe amount and frequency of the communication between

main function is removed from the source code, and compiled

different replicas are small, a single thread is enough to control

using the nvcc compiler with the DS-CUDA server library,

all the GPUs in the cluster.

which includes the body of a server as well as the CUDA libraryto generate a server executable.
Before running the user application, a DS-CUDA

Fig. 4 shows the strong scaling of the relative
calculation speed. 880 and 340 times acceleration is achieved
with 1,024 GPUs for 256 and 2,048 particles, respectively.

daemon should be ready for each node (Fig. 3). When theclient executable is started, the server executable is transferredfrom the client side to the server side and started via the DSCUDA daemon. Every time the client executable accesses thevirtual GPU in the client device, its data is transferred to theserver side via a network through the DS-CUDA server, and theserver accesses the physical GPU in the node.
When the function of redundant calculations is
activated for fault tolerance, each CUDA API is automaticallycopied to two identical APIs and executed in two servers.
Then the results from cudaMemcpy(DeviceToHost) APIs arecompared. If the results do not match, previous APIs areautomatically re-executed. When the migration functionis activated, data in GPU memory which is allocated bycudaMalloc() API are periodically backed up to the client side.
If unrecoverable error occurs on GPUs or network between
client and servers is disconnected, backed-up data aretransferred to a new GPU to continue the calculation.

14

Fig. 4 Strong scaling of performanceagainst the number of GPUs

Fig. 5 shows the time ratio of calculation andcommunication with 2,048 particles for 256, 1,024, and 4,096
GPUs. We made a performance model which fits up to 1,024

Overhead to support fault toleranceGPUs and applied it to 4,096 GPUs.

6

Users have only to set environment variables to activate
the fault tolerant function, but the overhead in calculation
time exists. Backing-up GPU memory to the client side every
few steps, known as checkpointing, causes performance
degregation, especially when the frequency of checkpoinging
is not optimal.
Fig. 6 shows the relative calculation time when the
period of checkpointing was changed [12]. The same applicationin the previous section was used for this experiment. Overclocked GPUs were used to artificially generate calculation
errors, and we chose data when the period is close to 900
seconds per one error, which is shown as crosses in Fig. 6.
A special DS-CUDA server, which is designed to generate an
Fig. 5 Ratio of calculation and communication

error every 900 seconds, is also used, and its data was plotted
as circles in Fig. 6.

In the performance model, the total calculation time
is divided into three parts: GPU calculation, CPU calculation,and communication. The communication is further
divided into two parts: Latency-oriented communication
time and bandwidth-oriented communication time. The
average amount of the communication is only 1.9 kbytes,
which is relatively small. Therefore, the latency-oriented
communication time is 57.5 msec for 1,024 GPUs, while
the bandwidth-oriented one is only 6.7 msec. In total, GPUcalculation time is still dominant, which is 471 msec. However
for 4,096 GPUs, GPU calculation time would be reduced to
118 msec, while latency-oriented communication time beincreased to 230 msec. The bottleneck would be caused by

Fig. 6 Overhead of calculation withfault tolerant function

the latency-oriented communication for this case. To avoid thebottleneck, one idea is to use multi thread programming like
OpenMP on the client side.

On the one hand, the calculation time becomeslarger when the period of checkpointing is small, since theoperation of backing-up GPU memory takes more time. On
the other hand, the calculation time also becomes larger when
the period of checkpointing is large, since the re-executionof previous CUDA APIs takes time as the possibility of errors
occurring becomes larger. In this configuration, the optimal
checkpointing interval was 100 seconds, and the overhead
was roughly 10% . Non-over-clocked GPUs are more reliable [10],and the overhead would be much less than this result.

15

DS- CUDA : A Handy Tool to Use GPUs in a Cloud Network

Combining tablet and GPU notebook

7

Concluding remarks

8

Fig. 7 shows a block diagram of a DS-CUDA system, which

The effectiveness of the idea of DS-CUDA to use multiple GPUs

composed of a non-GPU tablet on the client side, a GPU

was discussed as a cloud and visualization tool.

notebook on the cloud side, and a WiFi router between

DS-CUDA enables the use of many GPUs via athem. Table 2 shows the energy efficiency (Gflops/Watt)

network from a usual client computer. Acceleration ofand rendering speed (Frame/sec) of an MD simulation andcalculation and visualization can be achieved while utilizing

visualization using 5,832 particles[13, 14]. Power consumption ofboth of the tablet and the notebook are summed for the last

touch and accelerometer sensors on a tablet client device.
Moreover, a fault tolerant capability can be easily added.
Although mobile GPUs might have better energy

row.

performance and high-end GPUs might have better absoluteperformance, combining a client with cloud computers canachieve balanced performance in calculation speed, rendering
speed, and power efficiency.
Recent advancement of tablets makes them an ideal
tool for interactive simulation and visualization. Cloud tools
Fig. 7 Using a tablet and notebook with DS-CUDA

like DS-CUDA would be an interesting approach to realize suchrequirements.

Acknowledgements
This research was supported in part by the Japan Science andTechnology Agency’s Core Research of Evolutional Scienceand Technology research program:“Highly Productive, High
Table 2 Calculation speed and energy efficiency

Performance Application Frameworks for Post PetascaleComputing.” TSUBAME 2.5 supercomputer in Tokyo Institute ofTechnology was used in part for this research. Author thanks

The DS-CUDA system achieves 9.0 Gflops/Watt,

Atsushi Kawai in K&F Computing Research, Minoru Oikawa in

which is good enough compared with Green500 [15], though

Chiba University, Edgar Josafat Martinez-Noriega in University

the number is a little lower than that of a Notebook of 12.8

of Electro-Communications, and Kentaro Nomura and Kenji

Gflops/Watt. Note that it should not be compared with

Yasuoka in Keio University for the development and evaluation

Green500 numbers precisely, since it was calculated by aof DS-CUDA.

single-precision flops count and complicated function anddivision is counted as several operations.
Though the tablet itself is slow both in calculation
and rendering, combining a tablet and a notebook with DS-

References[1]

speed, energy efficiency, and rendering speed. Especially,

[2]
[3]

N. Luehr, A. G. Jin, and T. J. Mart´ınez, “Ab Initio Interactive
Molecular Dynamics on Graphical Processing Units

Although not shown in this table, a SHIELD tablet,which supports native CUDA, was found to achieve lower(GPUs),” Journal of chemical theory and computation,

numbers in calculation speed, energy efficiency, and rendering

Vol. 11, No. 10, pp. 4536–4544, (2015)

speed compared with the DS-CUDA system.

16

Nvidia GRID: http://www.nvidia.com/object/gridtechnology.html (Accessed in Jan. 2017)

rendering speed of 20 Frame/sec is smooth enough forinteractive simulation.

Amazon EC2: https://aws.amazon.com/jp/ec2/instancetypes/p2/ (Accessed in Jan. 2017)

CUDA achieved several tens of acceleration in calculation[4]

A. Kawai, K. Yasuoka, K. Yoshikawa, and T. Narumi,

“Distributed-Shared CUDA: Virtualization of large-scale
GPU systems for programmability and reliability,” Future
Computing 2012, Nice, France, (2012)[5]

DS-CUDA: http://narumi.cs.uec.ac.jp/dscuda/ (Accessedin Jan. 2017)[6]

S. Shi, C.H. Hsu, “A Survey of Interactive Remote
Rendering Systems,” ACM Computing Surveys (CSUR),
Vol. 47, No. 4, Article 57, pp. 1-29, (2015)[7]

J. Duato, A. J. Pena, F. Silla, R. Mayo, and E. S. QuintanaOrti, “Performance of CUDA virtualized remote GPUs inhigh performance clusters,” International Conference on
Parallel Processing, IEEE, pp. 365–374, (2011)[8]

T. Shimada, “PluginGPU Box : Virtual CUDA environment
for Windows PC,” Graduation thesis in Department ofCommunication Engineering and Informatics, University
of Electro-Communications, (2014)[9]

E. J. Martinez-Noriega, “Running CUDA on Android
Through GPU Virtualization,” GPU Technology
Conference, San Jose, P4160, (2014)[10] I. S. Haque, V. S. Pande, “Hard Data on Soft Errors: A LargeScale Assessment of Real-World Error Rate in GPGPU,”
10th IEEE/ACM International Conference on Cluster,
Cloud and Grid Computing, pp.691-696, (2010)[11] M. Oikawa, K. Nomura, K. Yasuoka, T. Narumi, “ Parallel
Molecular Dynamics Simulation with Replica-exchange
Method Using 1,024 GPUs,” IPSJ Transactions on
Advanced Computing Systems, Vol. 7, No. 4, pp. 1-14,(2014)[12] T. Narumi, M. Oikawa，E. J. Martinez-Noriega，K.
Yasuoka, “DS-CUDA: GPU Virtualization Middleware to
Support Migration Functionality,” IPSJ SIG Technical
Report, 2016-HPC-153(17), pp.1-6, (2016)[13] E. J. Martinez-Noriega, “High Performance Computing on
Mobile Devices through Distributed Shared CUDA,” GPU
Technology Conference, San Jose, S5290, (2015)[14] E. Martinez-Noriega, S. Yazaki, T. Narumi, “Performance
Evaluation of Remote CUDA Offloading on Mobile
Devices for Energy Efficient Systems,” submitted[15]

Green500 List: https://www.top500.org/green500/● TSUBAME e-Science Journal vol.15(Accessed in Jan. 2017)

Published 3/8/2017 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, Yuki ITAKURA
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/

17

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

15