GPUを用いた大規模シミュレーションと可視化

計算科学技術部門セッション
「シミュレーション可視化技術の最前線」
GPUを用いた大規模シミュレーションと可視化
小野寺直幸、青木尊之、下川辺隆史、杉原健太
東京工業大学学術国際情報センター
発表内容
•
研究背景と研究目的
•
大規模計算に適した計算手法と可視化
•
スーパーコンピュータ TSUBAME を用いた計算コードの性能測定
•
計算結果
•
•
東京都心部気流計算
•
車体周りのラージエディ・シミュレーション
•
卓球競技におけるピンポン球の飛翔計算
•
気液二相流計算
まとめ
研究概要
・Large-eddy simulationによる乱
流解析
・Graphics-processing unit(GPU)
を用いた大規模並列計算
複雑物体周りの Large-eddy
simulation
2009年 ∼
GTX285(2GBMem)
160x160x256, 4 GPU
http://www.top500.org/lists/2014/06/
GPUとCPUの比較
NVIDIA Tesla K20X
演算性能
Intel Xeon E5-2697 v2
1.31 TFlops (倍精度)
演算性能
0.52 TFlops (倍精度)
3.95 TFlops (単精度)
メモリバンド幅
250 GB/sec
メモリバンド幅
59.7 GB/sec
メモリ量
6 GB (GDDR5)
メモリ量
128 GB (DDR3)
CUDA コア数
2688 個
コア数
2 個 (24コア)
CUDAコア
CPUコア
差分法による非圧縮流れの解析
•
非圧縮性 Navier-Stokes 方程式
•
速度・圧力の時間発展方程式を計算
u
¯i =
u
¯ni
p¯
x2i
2 n+1
u
¯n+1
i
+
u
¯ni u
¯nj
1
+
xj
Re
u
¯i
=
xi
=u
¯i
p¯n+1
xi
t
2 n
u
¯i
x2j
yj+1
ij
xj
t
yj
ui− 12 ,j pi,j
vi,j− 12
yj−1
xi−1
xi
xi+1
スタッガード格子配置
メモリのRead/Writeが頻発し計算効率が悪化
格子ボルツマン法による流体解析
•
ボルツマン方程式
•
有限個の方向を持つ速度分布関数の時間発展方程式を計算
1
fi (x + ci t, t + t) = fi (x, t)
(fi (x, t) fieq (x, t))
方程式がシンプルで計算効率が良い
東京工業大学青木先生の講義資料より抜粋
D2Q9モデル
D3Q19モデル
差分法と格子ボルツマン法の演算密度の比較
•
差分法による非圧縮性流れの解析(FDM)
u
¯i =
u
¯ni
p¯
x2i
2 n+1
u
¯n+1
i
•
u
¯ni u
¯nj
1
+
xj
Re
+
=
2 n
u
¯i
x2j
ij
xj
NVIDIA M2050
t
能
理
u
¯i
xi
=u
¯i
e
p¯n+1
xi
性
論
in
oﬂ
t
v
ro
p
Im
ed
の
で
el
od
M
o
R
FDM≒1
LBM≒2
格子ボルツマン法による流体解析(LBM)
fi (x + ci t, t +
t) = fi (x, t)
1
(fi (x, t)
fieq (x, t))
Rooﬂineモデルを用いた理論性能
一般的に、格子ボルツマン法は差分法と比較して演算密度が高いため、
高い演算性能を得ることが可能
行列の反復計算を含む問題での計算速度の悪化
8
N=1283 / GPU
Time (s/step)
7
6億7千万格子
6
(320GPUs)
5
⇒4億格子
4
⇒2億格子
3
(96GPUs)
2
⇒1億格子
1
気液界面（VOF = 0.5）
気液界面の撹拌計算
0
(192GPUs)
0
50
100
150
200
250
300
350
Number of GPUs
１ステップあたりの計算時間
問題サイズが大きくなるに従って、Poisson方程式の反復解法の収束性
が悪化し、計算速度が低下
乱流の計算手法の分類
モデル
DNS
LES
特徴
すべての変動を計算
格子解像度以下の成分
をモデル化
計算負荷
大
中
解析精度・解像度
高い
・非常に多くの格子点数
が必要
中∼高い
・格子解像度を増やすこ
とでDNSに近づく
低い
RANS
変動成分をモデル化
小
・幅広い現象が計算可能
である反面、精度が低い
・モデルの影響が大きい
ラージエディ・シミュレーション
•
フィルター化されたNavier-Stokes方程式
物理量を格子で解像できる成分(GS)と格子解像度以下の成分(SGS)に分割し、SGS成分
をモデル化
u=u
¯ + u0
@u
¯i
@u
¯i u
¯j
+
=
@t
@xj
•
@ p¯
1 @2u
¯i
+
@xi
Re @x2j
@⌧ij
@xj
エネルギー
散逸
渦粘性モデルによるSGS応力
SGSのエネルギー散逸を分子粘性と同様にモデル化
⌧ij = ui uj
=
u
¯i u
¯j
2⌫SGS Sij
1
Sij =
2
✓
@uj
@ui
+
@xi
@xj
◆
GS
SGS
コルモゴロフの相似則
ラージエディ・シミュレーションの乱流モデル
Smagorinsky モデル、Dynamic Smagorinsky モデル
•
Smagorinsky モデル
幅広く使われる代表的なSGSモデル。モデル係数の値を一定値とするため、経験が必要
○ 計算がシンプル
⌫SGS = C4 |S|
2
•
：定数
Dynamic Smagorinsky モデル
△ モデル係数が一定(定数値)なため、
壁面近傍などで精度が低下
△ モデル係数の決定に経験則が必要
物理量に二種類のフィルターをかけることで、動的にモデル係数が決定可能
⌫SGS = C42 |S|
< Lij Lij >
C=
< Mij Mij >
○ 壁などの物体に対しても適用可能
△ 計算が煩雑
¯d
¯j u
¯ˆi u
¯ˆj
Lij = u
iu
ˆ¯ S¯ˆ △ 全計算領域での平均化が必要
¯ 2 |S|
¯ˆ2 |S|
¯d
Mij = 24
S¯ij 24
ij
: 計算領域全体での平均操作
→複雑物体周りの流れに適用できない
→大規模計算に不向き
ラージエディ・シミュレーションの乱流モデル
Coherent-structure Smagorinsky モデル
•
Coherent-structure Smagorinsky モデル
速度勾配テンソルの第二不変量を局所的に求めることで、渦粘性係数を決定
⌫SGS = C42 |S|
C = C1 |FCS |
3/2
FCS =
Q
E
○ 渦粘性係数を局所的に決定可能
→複雑物体周りの流れに適用可能
→局所的なメモリアクセスのため大規模計算に適している
速度勾配テンソルの第二不変量(Q)と
エネルギー散逸(ε)
格子ボルツマン法によるラージエディ・シミュレーション
•
Boltzmann 方程式
fi (x + ci t, t +
t) = fi (x, t)
streaming(並進移動)
•
1
(fi (x, t)
fieq (x, t)) + Fi
collision (衝突・減衰)
forcing term (外力)
サブグリッドスケールモデルの適用
渦粘性モデルでは、SGSの効果を分子粘性と同様に表現
=
3
1
+
c2 t 2
=
0
+
(LBMの緩和時間)
SGS
(粘性 = 分子粘性＋渦粘性)
平行平板間乱流の
速度勾配テンソルの第二不変量
格子ボルツマン法での物体境界面の表現手法
•
Interpolated bounce-back法
壁面上で鏡面反射条件を与え、逆方向に跳ね返す境界条件
fluid
(鏡面反射条件)
fluid
壁面上での
運動量交換
（
bounce-back rule
wall
wall
P.H.Kao, R.J.Yang, An investigation into curved and moving boundary treatments in the
lattice Boltzmann method, (2008)
Pierre Lallemand, Li-Shi Luo, Lattice Boltzmann method for moving boundaries, (2003)
Postprocessing
•
Paraview (http://www.paraview.org)
•
Povray (http://www.povray.org)
共に Open-source multi-platform
TSUBAME 2.0 の概要
System (58racks)
TSUBAME 2.0
Peak performance 2.4 PFlops
GPU :
NVIDIA Tesla M2050 x 4224
(1408Node 3GPU/Node)
Network：
QDR InﬁniBand x 2 (80Gbps)
Rack
(30nodes)
Compute Node
GPU×3
CPU 2
TSUBAME 2.5 の概要
System (58racks)
TSUBAME 2.5
Peak performance 17 PFlops(SP)
GPU :
NVIDIA Tesla K20X (6GB)
(1408Node 3GPU/Node)
Network：
QDR InﬁniBand x 2 (80Gbps)
Rack
(30nodes)
Compute Node
GPU×3
CPU 2
GPUによる並列計算
•
MPIライブラリを用いて、各計算領域の端部分を隣接ランクと交換
GPU
GPU
cudaMemcpy
cudaMemcpy
GPU CPU
CPU GPU
PCI Express
PCI Express
MPIによる領域分割
CPU
MPI通信
CPU CPU
QDR Inﬁniband
CPU
大規模並列計算での通信のオーバーヘッド
•
GPUによる並列計算での計算プロセスの流れ
１．GPUによる計算
Kernel function
2．CPUによるMPI通信
D2H
MPI
send & recv
H2D
ノード内の計算速度が速いGPUでは、通信時間の割合が大きくなる
MPIによるデータ転送の時間がオーバーヘッド
オーバーラップ計算による通信の隠
stream[0]
GPUのカーネル計算とCPUを用いたMPIの通信をオーバーラップ
stream[1]
•
Kernel function[0]
GPU
D2H
CPU
MPI
send & recv
H2D
MPI
CPU
Kernel function[1]
GPU
MPIの通信時間の隠が可能となり、良いスケーリングが期待
Strong scalability on TSUBAME 2.0
•
全計算量域の格子点数を固定：192x512x512, 192x2048x2048, 192x4096x4096
・オーバーラップ手法の導入により、
30%の速度向上
・良い強スケーリングが得られたこ
とより、より高速な解析が可能
Weak scalability on TSUBAME 2.0 and 2.5
•
プロセスあたりの格子点数を固定：196
256 256 / GPU
(N = 192 x 256 x 256)
Performance [TFlops]
1250
TSUBAME 2.5 (overlap)
TSUBAME 2.0 (overlap)
1000
TSUBAME 2.5
1142 TFlops (3968 GPUs)
288 GFlops / GPU
750
x 1.93
500
TSUBAME 2.0
149 TFlops (1000 GPUs)
149 GFlops / GPU
250
0
0
1000
2000
3000
Number of GPUs
4000
・TSUBAMEのアップデートによ
り、GPUあたり 1.93倍の高速化
・3968 GPUを用いた計算で、
1.14 PFlops(単精度)の非常に高い
実効性能を達成
東京都心部の大規模気流計算
計算条件
計算領域サイズ
格子解像度
格子点数
GPUあたりの格子点
10km
10km
500m
1 m立方格子
新宿駅
528 億格子
東京駅
10,080 10,240 512
160 160 512
とデータサイズ
52 MB/変数 (単精度)
GPU数
4032 台
境界条件
上空100mで10m/sの北風
渋谷駅
品川駅
東京都心部の大規模気流計算
計算データサイズ
計算条件
計算領域サイズ
格子解像度
格子点数
GPUあたりの格子点
10km
10km
500m
1 m立方格子
新宿駅
528 億格子
東京駅
10,080 10,240 512
160 160 512
とデータサイズ
52 MB/変数 (単精度)
GPU数
4032 台
渋谷駅
境界条件
上空100mで10m/sの北風211GB/変数と非常に大きいため、
データサイズが4032並列では
品川駅
各計算ノード内のローカルSSDを使用した並列 I/O 処理が必要
東京都心部の大規模気流計算
建物データ
計算条件
読み込んだ建物データ (株)パスコ
TDM3D 簡易版
計算領域サイズ
格子解像度
格子点数
GPUあたりの格子点
10km
10km
500m
1 m立方格子
新宿駅
528 億格子
東京駅
10,080 10,240 512
160 160 512
とデータサイズ
52 MB/変数 (単精度)
GPU数
4032 台
境界条件
上空100mで10m/sの北風
渋谷駅
品川駅
東京都心部1m解像度10ｋｍ四方気流計算 10,080x10,240x512(4,032 GPUs)
東京都心部1m解像度10ｋｍ四方気流計算 10,080x10,240x512(4,032 GPUs)
東京都心部1m解像度10ｋｍ四方気流計算 10,080x10,240x512(4,032 GPUs)
東京都心部1m解像度10ｋｍ四方気流計算 10,080x10,240x512(4,032 GPUs)
東京都心部1m解像度10ｋｍ四方気流計算 10,080x10,240x512(4,032 GPUs)
東京都庁前
東京都庁前の速度分布高さ200m
640 m
風向き(北風)
960 m
東京都庁前の速度分布高さ100m
640 m
風向き(北風)
960 m
東京都庁前の速度分布垂直断面
風向き(北風)
風向き
垂直断面図
(北風)
車体周りの流体計算
物体表面データ形式
•
STLデータと距離を用いた物体表現
STLデータ(三角形の頂点と法線ベクトルからな
るデータ)から、物体表面の距離関数を作成
車体周りの流体計算
計算条件
格子点数
3,072 1,536 768
4.2 mm
格子解像度
計算領域
60km/h
13m
6.5m
3.25m
GPU数
288 台
計算時間
約20時間(20万ステップ)
ピンポン球の飛翔計算
計算条件
格子点数
1,152 576 576
格子解像度
0.4 mm
速度
20 m/s
回転の種類
(i) 無回転
(ii) 横軸回転
(ii) Z-axis
y
z
20 m/s
x
(i) non-rotation
rotation
40mm
ピンポン球周りのラージエディ・シミュレーション(物体座標固定)
(i) non-rotation
(ii) z-axis rotation
ピンポン球周りのラージエディ・シミュレーション(物体移動)
> no spin
ピンポン球周りのラージエディ・シミュレーション(物体移動)
> top spin
ピンポン球周りのラージエディ・シミュレーション(物体移動)
> side spin
Two-phase ﬂow simulation
•
差分法による非圧縮性二相流解析
r·u=0
@u
+ (u · r) u =
@t
@
+ (u · r)
@t
=0
Positive area
1
1
2
rp + ⌫r u + F
⇢
⇢
@
= sgn( ) (1
@⌧
@
+ r · (u ) = 0
@t
⇢ = ⇢g + (⇢l ⇢g ) H( )
⌫ = ⌫g + (⌫l
⌫g ) H( )
|r |)
φ>0
Liquid
Interface
φ=0
Gas
Negative area
φ<0
符号付き距離関数を用いた
気液界面捕獲
Dam break simulation
P.K.Stanby, A.Chegini and T.C.D.Barnes (1998)
15 cm
⇢g
⇢l
⌫g
36 cm
⌫l
45
st
1.8 cm
Interface
thickness
1.18[kg/m^3]
997 [kg/m^3]
15.4 e-06[m^2/s]
0.89 e-06[m^2/s]
0.72 [N/m]
2 cells
72 cm
乾いた床へのダム崩壊では水面の先端が床をうように安定に浸水するが、
濡れた床の場合には水柱の先端が激しく乱れ、砕波する
Dam break
simulation
⇢g
⇢l
⌫g
⌫l
st
Interface
thickness
1.18[kg/m^3]
997 [kg/m^3]
15.4 e-06[m^2/s]
0.89 e-06[m^2/s]
0.72 [N/m]
2 cells
576x96x288
Bubble rising
unsteady simulation
⇢g
⇢l
⌫g
⌫l
24 cm
22 cm
st
Interface
thickness
1.18[kg/m^3]
997 [kg/m^3]
15.4 e-06[m^2/s]
0.89 e-06[m^2/s]
0.72 [N/m]
3 cells
・Symmetric wall boundary
!
Substitute air value and
vertical velocity 0.5 m/s
1cm
6 cm
まとめと今後の計画
•
GPUを用いた流体解析コードを開発し、TSUBAME 2.5において3968台
のGPUを用いた計算で、1.19 PFlopsと非常に高い演算性能を達成した。
•
格子ボルツマン法を用いた東京都心部1m解像度10km四方の気流計算に
おいて、10,080 x 10,240 x 512の大規模計算を実施し、高層ビル群に
より気流が乱される様子や幹線道路に沿って流れる「風の道」等の現象を
再現した。
•
気液二相流計算を実施し、細かな気泡が巻き込まれる様子を再現した。
•
今後の計画として、気液二相流計算において百億格子点以上の大規模計算
を実施する。

Download Report