What is GPU?
GP GPUGP GPU

Inexpensive

GPUスパコンTSUBAMEによるペタスケール物理シミュレーション

High Performance

Low Electronic Consumption
東京工業大学

学術国際情報センター

青 木尊 之Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

1Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

TOP 500 Ranking

TSUBAME Supercomputer

22013 November

MIC Xeon Phi

GP GPUGP GPUGPU K20X

100PF
K‐Computer
10.5 PF
Wire

10PF
４位

5位

TSUBAME 2.0

5位

17PF
Singleprecision

14位

2287.6TF

GPU K20X

1PF

Graph 500
No.3 (2011)

MIC Xeon Phi
41位 64位

100TF

25位 30位
7位

MEXT
MinisterAward(2012)9位

14位 109.7TF

49.5TF

10TF

TSUBAME

1.1

2.5

2.0

163.2TF

1.2
No.1 in Japan200720092011

GPU K20X

2013

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

GPU: Fermi core → Kepler core

Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

GP GPUTSUBAME 2.0

TSUBAME 2.5

NVIDIA Tesla M2050 ×4224

NVIDIA Tesla K20X ×4224

TSUBAME 2.5
Rack (30 nodes)
Performance: 122 TFLOPS
Memory: 2.28 TB

System (58 racks)
1442 nodes: 2952 CPU sockets,4264 GPUs

Memory: 116 TB(2 CPUs, 3 GPUs)
Performance: 4.08 TFLOPS
Memory: 58.0GB(CPU)+18 GB(GPU)

5GP GPU

Performance: 224.7 TFLOPS (CPU) ※ Turbo boost
5.562 PFLOPS (GPU)Total: 5.7 PFLOPS (double precision)
17.1 PFLOPS (single precision)Compute NodeCopyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

4Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

Green500 Ranking 2013 November
GP GPU

CUDA(Compute Unified Device Architecture)

GP GPU

CUDA is available for NVIDIA GPUs

■ GPU Computing (GPGPU) framework
C/C++ , FORTRAN Language + GPGPU extension(Compiler)
Runtime API
Libraries, Documentation,sample programs■ Currently CUDA 5.5 is released.
Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

8Domain decomposition

Bandwidth Bottlenecks

2D decomposition

GP GPUGP GPU■ Memory Hierarchy
Memorysub-domain #1

DDR3-1333

sub-domain #2x

Tesla K20X

32 GB/sec

CPUGPUGPU
Memory Interface 384bit

250 GB/sec(1) GPU→CPU

VRAM

1000BASE-T
0.125 GB/sec

PCI Express 2.0 x16
64 Gbps (8 GB/sec)

InfiniBand QDR
32 Gbps (4 GB/sec)

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

CPU
9Overlapping between
Computation and Communication

Sync

SendBuffer

MPI(3) CPU→GPU(2) CPU → CPUMPI

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Weather News

RecBuffer

Stream 1

Stream 2

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

10Japanese
GP GPU

Asynchronous Data transfer

Asynchronous Data transfer

PCI Express(cudaMemcpy)

11Copyright 気象庁（http://www.jma.go.jp/jma/index.html）© Global Scientific Information and Computing Center, Tokyo Institute of TechnologyGP GPUNext Generation

Weather Prediction

GP GPUCollaboration: Japan Meteorological Agency

Full 3-D Navior-Stokes Equation

u1
 u   u    P  2  u    (   r )  g  FtCloud Resolving Non-hydrostatic model
Compressible equation taking consideration of sound waves.

Meso-scalea few kmPhysical Process:

Tornado, Down burst
Heavy Rain

TyphoonGP GPU

Dynamical Process:Meso-scale Atmosphere Model:

2000 km

Atmosphere Model

Cloud Physics, Moist, Solar Radiation, Condensation,
Latent heat release, Chemical Process, Boundary Layer“Parameterization” including sin, cos, exp, …in empirical rules.

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

WRF GPU Computing

Full-GPU Implementation:
GP GPU WRF (Weather Research and Forecast)

Initial condition

WSM5 (WRF Single Moment 5-tracer) Microphysics*

Dynamics

GP GPU

Physics

CPURepresents condensation, precipitation and thermodynamic effects of latent heat release1 % of lines of code, 25 % of elapsed timeoverall improvement）

ASUCA⇒ 20 x boost in microphysics (1.2 - 1.3 x

Full GPU Approach

outputGPUWRF-Chem**

provides the capability to simulate chemistry and aerosols from cloudscales to regional ⇒ x 8.5 increase

Initial condition

Accelerator Approach

Dynamics

Physics

output

CPU*J. Michalakes, and M. Vachharajani: GPU Acceleration of Numerical
Weather Prediction. Parallel Processing Letters Vol. 18 No. 4. World
Scientific. Dec. 2008. pp. 531—548

GPUJ. Ishida, C. Muroi, K. Kawano, Y. Kitamura, Development of a new nonhydrostatic
model “ASUCA” at JMA, CAS/JSC WGNE Reserch Activities in Atomospheric andOceanic Modelling.

**John C. Linford, John Michalakes, Manish Vachharijani, and Adrian Sandu. Multi-core acceleration ofchemical kinetics for simulation and prediction, proceedings of the 2009 ACM/IEEE conference onsupercomputing (SC'09), ACM, 2009.
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

15

Entire Porting Fortran to CUDA

GP GPU

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

TSUBAME 2.0 (1 GPU)

16GP GPU Rewrite from Scratch

FortranC/C++Original code
at JMAChangingarray order

z,x,y (k,i,j)-ordering

x,z,y (i,k,j)-ordering

CUDA
12x speedup
50x speedup
to compared
comparedto Xeon
toX5670 1 socket
1 CPU coreGPU code
x,z,y (i,k,j)-ordering 1 Year by one Ph.D student
Introducing many optimizations, overlapping the computation with thecommunication, kernel fuse, re-ordering kernel, . . .
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

17

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

18

ASUCA Typhoon Simulation
5km-horizontal resolution 479×466×48

TSUBAME 2.0 Weak Scaling

GP GPUGP GPU

145.0 Tflops
Single precision

76.1 Tflops
Doublele precision

Fermi core Tesla
M2050

3990 GPU■ SC’10 Technical Paper
Best Student Award finalist

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

19Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

20

ASUCA Typhoon Simulation
500m-horizontal resolution 4792×4696×48
Using 437 GPUs
GP GPUGP GPUAir flow in a wide area of Tokyo
by using Lattice Boltzmann Method

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology21

22

LES (Large-Eddy Simulation)

Lattice Boltzmann Method
f i
1eq
 e i  f i   f i  f i
tCopyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

GP GPUGP GPU

15
1216

3
93

2f i eq  wi 1  2 ei  u   4 ei  u   2 u  u 2c2c
 c


11

5

3

8

7

2

Energy spectrum

Relaxation timefor LES model

1

10

94

Strongly Memory Bound Problem:

17Collision step:

Streaming step:

14

613

Molecular viscosity andEddy viscosity

18relaxation time

GS: kinetic viscosity (SGS)
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

23

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

SGS

SGS Model

Computational Area
GP GPU

Smagorinsky model○ Simpleᇞ inaccurate for the flow with wall boundaryᇞ emperical tuning for the constant model coefficient

Dynamic Smagorinsky model○ applicable to wall boundary
ᇞ complicated calculation
ᇞ average process over the wide area
→ not available for complex shaped body→ not suitable for large-scale problem

Coherent-Structure Smagorinsky model*H.Kobayashi, Phys. Fluids.17, (2005).

→ model coefficient determined by the second invariant of
the velocity gradient tensor

GP GPUMajor part of Tokyo
Including Shnjuku-ku,
Chiyoda-ku, Minato-ku,
Meguro-ku, Chuou-ku,

ShinjyukuTokyo

10km×10km
Building Data:
Pasco Co. Ltd.
TDM 3D

Shibuya

Shinagawaᇞ model coefficient○ applicable to wall boundary○ model coefficient is locallydetermined.

Map ©2012 Google, ZENRIN

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology26

Weak Scalability
GP GPU

600 TFLOPS
on 4000 GPUs

27

Performance [TFLOPS]TSUBAME 2.0 → 2.5▲ TSUBAME 2.5 (overlap)● TSUBAME 2.0 (overlap)

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

28

DriVer: BMW-Audi
GP GPU

Lehrstuhl für Aerodynamik und Strömungsmechanik
Technische Universität München

GP GPUTSUBAME 2.5
1142 TFLOPS(3968 GPUs)
288 GFLOPS/GPUTSUBAME 2.0
149 Tflops(1000 GPUs)
149 GFlops / GPU×1.93
Number of GPUs

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

30

Development New Materials
GP GPUGP GPU

Mechanical Structure

Material Microstructure

Low-carbon society
Improvement of fuel efficiency by
reducing the weight of transportation

Developing lightweight strengthening
material by controlling microstructure
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

31

Phase-Field Model
GP GPUThe phase-field model is derived from non-equilibrium
statistical physics and f = 0 represents the phase A and f = 1for phase B.

Dendritic Growth

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Al-Si: Binary Alloy

GP GPUTime evolution of the phase-field （Allen-Cahn equation)

Phase B
Phase A
Phase-field 1Time evolution of the condensation: cdiffusive interface
with finite thickness
0Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Finite Difference Method
GP GPU

Requirement for
Peta-scale Computing
Previous works

Phase Field : Condensation : c

19 points to solve i , j , kGP GPU

×1000 large-scale computation

2D computation

on TSUBAME 2.0

7 points to solve ci , j , kz  k 1

3D computationzkz  k 1

Single dendrite
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

35～ mm scale

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

36

Weak scaling

Domain decomposition
2D decomposition

4096 x 6400 x 12800
4000 (40 x 100) GPUs
16,000 CPU cores□ GPU-Only

GP GPU（No overlapping）◯ Hybrid-YZ
sub-domain #1（y,z boundary by CPU）

sub-domain #2x▲ Hybrid-Y（y boundary by CPU）

GPUHybrid-Y method2.0000045 PFlops(1) GPU→CPUCPU

GP GPUPCI Express(cudaMemcpy)

GPU: 1.975 PFlops
CPU: 24.69 TFlops(3) CPU→GPU

single precision(2) CPU → CPU

MPIEfficiency 44.5%(2.000 PFlops / 4.497 PFlops)•

Mesh size: 4096 x160x128/GPU
NVIDIA Tesla M2050 card / Intel Xeon X5670 2.93 GHz on TSUBAME
Copyright ©2.0
Global Scientific Information and Computing Center, Tokyo Institute of Technology
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

3738TSUBAME 2.0 → 2.5 性能向上
Mesh size of a subdomain(1GPU + 4 CPU cores): 4096 x 162 x 130

4,096 x 5,022 x 16,640

GP GPU(3,968 GPUs+15,872 CPU cores)

TSUBAME 2.5
3.406 PFLOPS

×1.75
TSUBAME 2.0
2.000 PFLOPS
4,096 x 6,480 x 13,000(4,000 GPUs+16,000 CPU cores)

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Dendrite Analysis

Simulation for Granular Materials
GP GPUGP GPUDEM (Distinct Element Method)

Normal direction

Contact interaction
Viscosity

SpringTangential direction
SpringFriction

Fij   kxij  xij
Viscosity
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

41

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of TechnologyParticles Across Boundaries

DEM Simulation on Multi-GPUGP GPU

 Load balance among GPUs to keep parallel efficiency

GP GPU Particle search on GPU to move boundary De-fragmentation of GPU Memory
Many particles

New boundary line

2th search

Old boundary line

1th search

Few particles

Node1

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

De-fragmentation

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Optimization of De-fragmentation frequency
GP GPU

Node111

2GP GPU

Computational time increase with linear by re-numbering.
111

Without sortingSome particles are sent to neighbors.
1Node213Time variation

1Every 2000
Every 20
Every 1000
Every 100

Some particles newly come.
1213111Re-numbering is executed on CPU.
1111

1Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Benchmark Test

Neighbor Particle List
GP GPU

Computational field are divided into 64 sub-domains.
Single GPU is assigned to each sub-domain.

GP GPU

Linked-list methodLocal domain
03606NULL
87 percent of memory usage is reduced
compared to regular neighbor list.
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

3Granular Convey

粉体の撹拌シミュレーション

64 GPUs, 5 million particles

64 GPUs, 5 million particles

Gas-Liquid Two-Phase Flow

Level-Set method (LSM)
GP GPU

Mesh Method (Surface Capture)■ Navier-Stokes solver：Fractional Step
■ Time integration：3rd TVD Runge-Kutta■ Advection term：5th WENO
■ Diffusion term：4th FD■ Poisson：MG-BiCGstab■ Surface tension：CSF model■ Surface capture：CLSVOF(THINC + Level-Set)

Particle Method
ex. SPH
Low accuracy< 106-7 particles

GP GPUThe Level-Set methods (LSM) use the signed distancefunction to capture the interface. The interface isrepresented by the zero-level set (zero-contour).

: Level-Set function(distance function): Heaviside function

High accuracy > 108-9 mesh points

Re-initialization for Level-Set function

notsplash

Advantage : Curvature calculation, Interface boundary
Drawback : Volume conservation

Numerical noise and unphysical oscillation
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of TechnologyContinuous Surface Force(CSF) model by Brackbill, Kothe and Zemach (1991)

Anti-diffusive Interface Capture
GP GPUThe interfacial surface force is transformed to a volumeforce in the region near the interface via a delta function

THINC (tangent of hyperbola for interface capturing) Scheme

GP GPU[ Xiao, etal, Int. J. Numer. Meth. Fluid. 48(2005)1023 ]・VOF(volume of fluid) type interface capturing method
・Flux from tangent of hyperbola function
・Semi-Lagrangian time integration

Curvature
Surface tension
forceFig. Takehiro Himeno, et. Al., JSME,65-635,B(1999),pp2333-2340

Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

Normalvector

Surface tension represented
by volume force
Approximate delta function

・1D implementation can be applied to 2D & 3D → Simple

・Finite Volume like usage＊ THINC is the method how to compute flux
→ 3 krenel (x, y, z) can be fused to 1 kernel. Merit in memory R/WCopyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of TechnologySparse Matrix Solver1Ax = buSet k  0 r0  p0  M 1 b  Ax0 for    p  tPre-conditioner:
Incomplete Cholesky,ILU, MG, AMG,
Block Diagonal Jacobi

Non-zero Packing:
CRS → ELL, JDL

MG V-Cycle
L0 f 0  S 0
Poisosn Eq.

Restriction

GL1v1  R1

AF ： Fine Matrix
AC ： Corse Matrix
R : Restriction
P : Prolongation

Correction Eq.

Gq , M Aq M Aq , M Aq 1kk

11kkkx k 1  x k   k pk  k qk
rk 1  qk  k M 1 Aqkifrk 1 , rk 1    2 b, b exit;

r0 , rk 1 k r0 , M 1 Apk
pk 1  rk 1   k  pk  k M 1 Apk k 

Lv RCopyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

56

Milk Crown

n+1 Step

L v  R1L2 v 2  R 2
Prolongation

L3v 3  R 33

Restriction

4k Prolongation

Restriction

G30qk  rk   k M 1 Apk

ProlongationAC = RAFP

L2 v 2  R 2
3 3

r , M Ap 
11 1Restriction

G255

GP GPU

Smoother : Red & Black ILUn Step1

r0 , rk 

loop end

Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

G0k GP GPUCollaboration withMizuho Information & Research Institutefor k = 0; k < N; k++;××××××× ×× ×××× ××××× ××××× ×××××× ×××××× ×××××× ××××××× ×××××××××××× ××××× ××××× ×××××× ××××××× ××××× ××××× ××××× ×××××× ×××××××Krylov sub-space methods:
CG, BiCGStab, GMRes, , ,BiCGStab + MG

GP GPUProlongation

L v  R4
4 4Copyright © Takayuki Aoki / Global Scientific Information and Computing Center, Tokyo Institute of Technology

57

Drop on dry floor

60Industrial Appl. Steering Oil

Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

SUMMARY■ Peta-scale GPU applications successfully
running on TSUBAME2.0/2.5
・ Meso-scale weather prediction
・ Lattice Boltzmann Method
・ Gas-Liquid Two-phase flow … on going・ Phase-Field model for developing newmaterials■ GPU Green computing : less electrical power
to get application results.
Copyright © Global Scientific Information and Computing Center, Tokyo Institute of Technology

63

GP GPU62