フォーラム「粒⼦法による流体シミュレーションの並列化技術」九州⼤学 情報基盤研究開発センター 2017年8⽉21⽇

TSUBAME3.0
2,160 NVIDIA P100 GPUs

空間充填曲線を⽤いた動的負荷分散による12.15 PFLOPS (Peak)

⼤規模粒⼦法シミュレーション

東京⼯業⼤学 学術国際情報センター
⻘⽊ 尊之Pascal core GPU:

Tesla P100

3,584 cores

国際会議 Particles で発表のある粒⼦法
DEM, SPH, MPS, MPM, LBM

3Particle Distribution

DEM (Discrete Element Method)
Normal direction

Viscosity

Contact interaction

4

The load imbalance problem among subdomains
Many particles

Spring

No particle
Tangential direction
SpringFriction

ViscositySub-domains with High Aspect ratio
Slice Grid Method

The passive particles under vortex velocity field using 64 GPUs

Golf Bunker shot using 16.7 millions particles
with 64 GPUs

バンカーショット計算

Complex Shape Objects
Non-spherical model described by rigidly connected sphere particles

∆ /

∆∆

∆ / ∆∆ /
∆⁄

∆ /

∆)∆

∆ /

∆∆ /
∆∆

∆

∆ /

∆ /

∆∆∆
∆ /

∆∆
1112

Multiple GPU ScalabilityConditions
Particles: 2 × 106, 1.6 × 107, 1.29 × 108
Domain Decomposition : Dynamic load Balance using Slice Grid Method
Time-Integration: 2-stage
Runge-Kutta

13

Space Filling Curves

Space Filling Curves• Construction of a tree so that each leafs contain minimum number of particles• Each leafs is numbered by following the rules of space filling curves.

 The connections of leafs byPeano curve The domain decomposition
by Peano curve

A tree-base cell refinement

空間充填曲線による動的負荷分散
ペアノ、ヒルベルト、モートン空間充填曲線による 2 次元動的負荷分散
PeanoHilbert

Morton空間充填曲線を⽤いた領域分割
 粒⼦数 / 領域番号テーブルの利⽤
 ツリーの構築の際に各リーフ領域内の粒⼦数を計測するが、粒⼦配列に対するリーフ内の粒⼦の探索は探索時間が全粒⼦数に⽐例するためコストが⾼い．

 予め計算領域内に最⼩リーフと同サイズの空間格⼦を⽣成しておき，各格⼦に格⼦内の粒⼦数を登録した粒⼦数テーブルにより探索時間を削減する．

Scalability on TSUBAME2.5SPH （Smoothed Particle Hydrodynamics）

Comparison of 3 different SFC
with slice-grid method

Short-range Interaction （SPH: Smoothed Particle Hydrodynamics）
Each particle has its physical properties and interactions among particles in the range of kernel radius
are calculated at every time step.

Physical property

Kernel function W

Normalization
PeanoHilbert

MortonSlice-gridParticle

Discretization: mass: density

Kernel radius

A Dam Break Simulation

Description of the Object Shape

Arrangement & setting of parameters

Objects formed by particles generated from CAD data

12 m

CAD

4.8 mParticles

2.2m

Water
10 m

6m

Object

0.8 m

A Dam Break Simulation

Tsunami Debris Flow
Including a lot of floating objects
Computational conditions•••••••••

72 M particles with 80 GPUs

Simulation area : 180 m × 160 m × 20 m
Total particles: 117,561,285
Particles / object : 19 ~ 472
Objects: 10368
Physical time (s) : 10.0
Steps: 20,000
GPUs
: 256
Simulation (h): 100Initial height (m) : 10.0

Initial state of each objects

Representation of object shapes
Particles

CAD Data245

191944151

71472124309144410315

TSUBAME2.5を⽤いた⼤規模計算

Dynamic Domain Decomposition
Using the improved SPH with 87Million particles with 256 GPUs

流動層の⼤規模シミュレーション
従来の固気・固液シミュレーション

直接相互作⽤による固気・固液シミュレーション

Fluid-solid Two Phase Flows-- Fluidized Bed -⾼精度
⾮球形粒⼦
乱流

 固体粒⼦は計算格⼦より⼩さいラグランジュ粒⼦
 球形粒⼦を仮定した相互作⽤モデル

 固体粒⼦は計算格⼦よりずっと⼤きい 流体−固体間の直接相互作⽤を計算

27

LBM-DEM coupled method

MRT型の格⼦ボルツマン法■ Solid fraction in each cell:■ Bounce back boundary condition∆ ,∆Δ ,,,,,,∆2,

ローカルな物理量で計算

̅

コヒーレント構造関数
1212：速度分布関数：変換マトリックス：緩和時間マトリックス

M(粘性＝分⼦粘性＋渦粘性)
̅:,乱流モデル（コヒーレント構造Smagorinsky）∗,,

隣接点に速度分布関数をコピー,/Δ/：速度モーメント: 速度分布の速度

緩和時間マトリックス を粘性 ∗ を⽤いて修正

渦粘性係数を局所的に決定・物体周りの流れに適⽤可能
・⼤規模計算に適⽤可能

Kobayashi, Hiromichi. "The subgrid-scale models based on coherent structures for rotatinghomogeneous turbulence and turbulent channel flow." Physics of Fluids 17.4 (2005): 045104.

D3Q27速度モデル

30

AMR (Adaptive Mesh Refinement)

格⼦解像度が異なる格⼦間の補間

Large-scale simulation flow around a sphere

1点空間補間各軸に線形補間 (1,2,4, または 8 点から)
時間補間
線形補間 (1つまたは2つの時刻から)レベル変換*
速度勾配テンソルが連続となる変換
Domain size
Mesh resolution (Δx/d)
Blockage ratio (d/D)
Total mesh numberGPU

D×D×2D
0.001 min.
0.512 max.
0.01529
2.7×108
NVIDIA K20X × 64

2点4点

空間補間
粗粗野化：

,

細分化：

,,2細,t*Filippova Olga and Hnel Dieter. Grid Refinement for Lattice-BGK Models.

時間補間

Journal of Computational Physics, Vol. 147, No. 1, pp. 219-228, 1998.

31

格⼦解像度ごとの時間ステップの設定

32

リーフの外殻格⼦のアクセスパターン

LBMの速度分布関数の移動速度 は⼀定

AMRのリーフ： 16×16×16等のパッチ⇒格⼦解像度に応じて時間刻み Δ を変化させる

橙：⾯上の外殻格⼦

格⼦は6⾯・12辺・8頂点の26通り
各カーネルの隣接点参照が27通り
全702通りの参照パターンが存在

緑：辺上の外殻格⼦

全⼿動の実装は⽣産性・保守性に問題あり

C++テンプレートによる実装
⻘：頂点上の外殻格⼦

空間充填曲線による領域分割

34

通信最適化① 通信する格⼦の最適化

袖領域リーフの全体

35

必要最⼩限通信する格⼦
2次元：8通り
3次元：26通り

36

通信最適化③ 格⼦ボルツマン法の分布関数

通信最適化② 時間ステップ
値が更新された格⼦のみ通信する2段の解像度の場合，• 偶数の時間ステップ：全格⼦を通信• 奇数の時間ステップ：細かい格⼦のみ通信

黒⾊：送る
⽩⾊：送らない

D2Q9LBM

C++テンプレートによる実装
偶数ステップ奇数ステップ
37

Drag Crisis

Interpolated Bounce-back移動境界
Re=2×105

Near Re = 105 〜 106,••38

Re=5×105
流体の格⼦点Drag coefficient drastically decreases
Vortexes translate past the sphere

格⼦点の間にある壁を2次精度で扱う境界条件

固体内部の格⼦点

wall

11J. Almedij, 2008
Morrison, 2013present method

壁が近い場合2,

Δ ,’CD,

0.1
1.E+05

Re,

Δ

2,

12Δ

Δ ,

2·

1.E+0639

固体粒⼦に作⽤する流体⼒

検証：移動する球に作⽤する抗⼒
wall

壁で反射される速度分布関数 の運動量変化から評価,,

1

球を⼀定速度で流体中を移動させ，抗⼒係数を測定

2,

Galilean invariant momentum exchange method*40

Δ ,0.02,

Δ

20 Δ
102

壁を含む格⼦点におけるの総和から，物体に作⽤する流体⼒とトルク,Δを計算

抗⼒係数2,,∈∈*Wen, Binghai, et al. "Galilean invariant fluid–solid interfacial dynamics in lattice Boltzmann simulations.“, Journal of Computational Physics 266(2014): 161-170.
41/4: 抗⼒: 流体の密度: 粒⼦の移動速度: 粒⼦直径MRT
MRT+LES
Experiment

101

CD’

100

10-1 -1
10

100

101

102

103

Re104

105

10642

複数GPU計算

検証：球形粒⼦の沈降計算

velocity [m/s]015 mm 15Δ
120 mm
100 mm 100 mm 160 mm
Δ
1.0 mm1120 kg/m
970, 965, 962, 960 kg/m

-0.04

-0.08

-0.12
01.5, 4.1, 11.6, 31.9全体の計算領域

Re=1.5Re=4.1
Re=11.6
Re=31.9
12time [s]3

ハロー領域の通信
全ノード間総和計算• 各⼩領域を1台のGPUに割り当て，LBMの計算を実⾏

4• 各GPUがすべての粒⼦を計算• MPI_Allreduceを⽤いて各GPUで計算された流体⼒の総和

粒⼦の沈降速度の時間変化

Ten Cate, A., Nieuwstad, C., Derksen, J. and Van den Akker, H.: "Particle imaging velocimetry experiments and lattice-Boltzmann
simulations on a single sphere settling under gravity'", Physics of Fluids, Vol. 14, No. 11, pp. 4012-4025 (2002).

⼆次元領域分割

4443

Spouted Bed

DEM
68,921 particles

FluidLBM

48 GPUs
512 512 1680
1.2 kg/m

288x288x576 mesh

1.5ΔΔ8 GPUs Tesla K80
512 512 1024 mesh1000 kg/m

81,920 sphere particles
Diameter 0.5 mm
100 kg/m
1.67 sΔ

計算結果

たわら型粒⼦を含む流動層

1.510m /s

5512 particles
85 kg/m
1.33 s
Δ0.03,
0.9

80 mm

0.133 mm
1.33 s
0.5 m/s

480 mm

円柱形粒⼦

流体（空気）

Δ
Δ

0.05 mm1.67 sSolid

10,300 cubic particles
185,800 particles
2500 kg/m

225 750 3600 mesh
54 GPU
1.2 kg/m

10 m /s
1.5 m/s

2.6 mm
6.0 mm
半径1.3 mmの球を9個並べた粒⼦モデル

固体粒⼦の
初期配置

0.5 m/s
30 mm
47

100 mm

流体の鉛直⽅向の速度48

Falling Gingko leavesModeling of Gingko leaf

RealPolygons

Level SetParticles(Collision detection)
49

Summary 空間充填曲線を⽤いた動的負荷分散による⼤規模粒⼦計算の実現
 格⼦ボルツマン法による⾮球形固体粒⼦との直接相互作⽤計算
 数10GPUの⼤規模計算• ⾮球形固体粒⼦の沈降計算• 噴流計算• 舞い落ちる多数のイチョウの葉の計算
 格⼦ボルツマン法のAMR実装• リーフの外殻格⼦の隣接アクセスの最適化• 空間充填曲線による領域分割• 分割領域間通信の最適化