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 ofl 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)) Rooflineモデルを用いた理論性能 一般的に、格子ボルツマン法は差分法と比較して演算密度が高いため、 高い演算性能を得ることが可能 行列の反復計算を含む問題での計算速度の悪化 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 InfiniBand 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 InfiniBand 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 Infiniband 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 flow 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の大規模計算を実施し、高層ビル群に より気流が乱される様子や幹線道路に沿って流れる「風の道」等の現象を 再現した。 • 気液二相流計算を実施し、細かな気泡が巻き込まれる様子を再現した。 • 今後の計画として、気液二相流計算において百億格子点以上の大規模計算 を実施する。
© Copyright 2024