7th Symposium

学際大規模情報基盤共同利用・共同研究拠点公募型共同研究 平成２７年度採択課題

jh150048-NA26

青木尊之 (東京工業大学 学術国際情報センター 教授)

動的負荷分散によるGPUスパコンを用いた粒子法の大規模シミュレーション手法の開発 II − 流体・構造連成計算への適用 −改良型粒子法(井元、田上 2014)

研究背景と目的

流体と構造物が共存する現象は実際に非常に多く、さまざまな分野で流体構造連成解析の重要性が広く認識されている。粒子法は複雑形状や流体構造連成計算を容易に行うことができる利点を持っているが、粒子法の並列計算には不均一な粒子分布が発生するためにプロセッサ間の負荷分散が著しく低下するなど大規模計算化が難しい問題がある。複数の球形粒子を連結させ、多数の球形粒子で表現する複雑形状の構造物とSPH(Smoothed Particle Hydrodynamics)を用いた流体計算との大規模流体−構造連成計算を、演算性能、メモリバンド幅や電力効率の点で有利なGPU (Graphics Processing Unit)を用いて高効率に実行するためのアルゴリズムの開発と実問題の大規模シミュレーションを行う。

＜ 改良型粒子法の陽解法の場合の計算手順 ＞－ 有限差分法の一般化としてとらえる新しい粒子法 －1. 予測子を計算:
2. Birch-Murnaghan’s の方程式から計算される仮の圧力を計算:

 2 次のスパイキー型のカーネル関数を利用

3. 仮の圧力をもとに粒子座標を修正:
4. 補間により圧力を再計算:

DEM (Discrete Element Method)と流体シミュレーションのため研究計画 紛体シミュレーションのためのの SPH (Smoothed Particle Hydrodynamics) 法に対し、GPUスパコンにおいて大規模シミュ: カーネル半径一階導関数の例

レーションを可能にするための計算手法の開発を行う。

5. 修正子により流速を計算:

 複数の球形粒子を剛体連結させた非球形DEMのGPU計算の効率化と通信負荷の低減化手法を開発する。
 第一段階として、スライスグリッド法による動的負荷分散を流体構造連成の大規模GPU計算にも適用し、演算や通信コストの問題を定量的に評価する。

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

 第二段階は空間充填曲線による動的負荷分散について検討し、細分化の段数の制御、データ構造の管理コスト、領域間通信のコスト、リーフのマイグレーションのコスト等を検証する。
 1000GPU以上の大規模アプリケーションを実行し、高い並列実行性能を達成できる動的負荷分散法を開発する。

粒子法による流体-構造連成計算手法

 2次元領域分割までが限界（分割手続きが複雑） 一つの分割領域が複数の分割領域と接触
 分割領域が高アスペクト比（粒子が分割領域境界を横切り易い：ノード間データ通信量の増加）
スライスグリッド法による動的負荷分散によるDEM計算は、弱スケーリングですら 512GPUで50%程度まで並列化性能が低下。

非球形粒子による個別要素法

32

 粒子間で 4 種類の相互作用計算を行う。

41

DEMに対する強・弱スケーリング

Spring

1. 流体粒子-流体粒子間 (流体計算)

Friction

粉体を有限の大きさを持つ粒子で表現し、周囲の接触している粒子から力を受Viscosity
けて個々の粒子が運動する計算手法である。粒子間の接触相互作用はバネとダンパでモデル化され、接線方向に作用する力は粒子間の摩擦に相当する。
Viscosity

計算領域を8分木または27分木で再帰的に分割し、全てのリーフの中の粒子数が所定の数以下になるまで細分化を続ける。計算領域内のリーフを空間充填曲線でたどり、通過したリーフ内の粒子数の積算により領域を分割。

改良型粒子法(上述参照)を用いる。

2. 流体粒子-壁粒子・物体構成粒子間 (境界条件)
壁粒子は静止した境界条件、物体構成粒子は移動境界条件として扱う。

Spring 非球形粒子モデルの導入

3. 物体構成粒子-流体粒子間の相互作用計算

球形粒子モデルでは、粉体同士が噛みあうような複雑な相互作用を表現できない。複数の球形粒子を剛体連結させて、非球形粒子をモデル化する。以下のことから球形粒子に比べて計算コストが高い。

粒子 i への圧力勾配による流体作用力を与える：(v : 現在の速度 v* : 仮の速度)

4. 物体構成粒子-物体構成粒子間の相互作用計算
個別要素法(DEM)による衝突力を与える。

• 計算に必要な粒子数の増加• クォータニオンによる回転管理• 重心に作用する合力とトルクの総和計算

自身と隣接セル内の粒子とのみ相互作用O(N2)からO(N)へと計算量を削減

 リンクリスト法によるメモリ使用量の削減
空間格子には粒子を一つだけ登録し、セル内で粒子同士がリンクリスト構造を生成
小領域
66

計算領域
144 m × 160 m × 60 m
物体構成粒子数1000 / 物体
物理時間 (秒) 6.8
GPU 数
256

 3次元領域分割が可能
 粒子分布の変化があったところだけ変更可能
 毎ステップの領域再分割が容易 隣接ノードとの接続関係が明確
 通信と計算のオーバーラップが可能 領域形状が複雑（隣接領域との通信が複雑）
 リーフの隣接ノードへのマイグレーションの負荷が不明

最大収容数を想定して静的メモリを確保
メモリ消費量が著しく増加

0 サスペンションフロー計算

 空間充填曲線による動的領域分割の利点／欠点

 近傍粒子探索リストの利用

3

流体構造連成の大規模GPU計算

Morton Curve

Hilbert Curve

GPUにおける高速化手法

0

壁粒子流体粒子物体構成粒子

 壁、流体、物体のすべてを粒子により表現する。

Performance [s-1]（東京工業大学）（海洋研究開発機構）（東京工業大学）（東京工業大学）（東北大学）（東北大学）（九州大学） Discrete Element Method

粒子法としては、渦を非常に良く捉えている。

x 方向に粒子数が一定になるように領域を分割し、分割された各小領域をy方向に独立に分割する。

研究グループ
青木 尊之
西浦 泰介
都築 怜理
渡辺 勢也
森口 周二
今井 陽介
井元 佑介

Re = 1,000 スライスグリッド法の問題点

 分割領域の境界周辺リーフと内側リーフに分け、境界周辺リーフを先に計算した後に内部リーフの計算とノード間通信をオーバーラップさせ、その有効性を検証する。

代表者
副代表者
課題参加者
課題参加者
課題参加者
課題参加者
課題参加者

2次元キャビティ・フロー問題による精度検証 (Ghia らとの比較)

粒子数
物体数
ステップ数
計算時間

87,430,000
2304 個
13600
43.1 hPeano Curve 不均一形状の物体を多数含む大規模津波シミュレーション

3

符号付距離関数による物体表現
 計算開始前にCADデータから符号付き距離場を生成し、格子データに格納し、計算時にはデータを参照するだけで距離判定を行う。
 符号付距離関数を用いて、流体計算用の粒子群データも生成できる。

NULL

メモリ消費量を大幅に削減 (DEMなら1/8、SPHなら1/64)

 セル番号をキーとした粒子ソートによる高速化
 同一セル内の粒子がメモリ上に連続する様ソートを実行 し 、 GPU におけるコアレッシングを実現

Particle

4,500

3,500
Elapsed time [msec]

 140万個の粒子を用いたSPH法によるダム崩壊を、単体GPU 上で実行した場合の計算時間の比較

3,000

••
•
•

2,500
2,000
1,500
1,000

約 8.6 倍に高速化

500001000

20003000Time step

CADデータ

符号付き距離関数

2次元スライスグリッド法を用いた場合の弱スケーリング

非球形DEMの大規模撹拌計算

セル番号による並び替えあり
セル番号による並び替えなし

4,000

Positive area Negative area>0<040005000

計算粒子数 : 6,205,952 個
非球形粒子数 : 299,584 個
ステップ数: 800,000
計算時間: 13.7 h

1,600

サイズ 構成粒子数 立方体数
3 mm8262,144

5 mm

56

32,768

9 mm

296

4,096

17 mm

1,352

512

33 mm

5,768

64

52.8 %
27.3 %

動的負荷分散法と物体管理法の改善が必要と判明

1,000
800
600

4GPU

blockIdx.x=2
ペアノ曲線

ヒルベルト曲線32GPU

256GPU

モートン曲線

スライスグリッド法

空間充填曲線による動的負荷分散をSPHの複数GPU計算に適用(10万粒子、24GPUの検証例)

200
Memory bandwidth[GB/sec]

1,2000 シェアードメモリを用いたリダクションで総和計算
2048 objects

150

1thread1block100

構成粒子数が多いほどリダクションの効果大

ノード内のReduction

1,400200 １つの非球形粒子の総和計算を1blockが担当

２０４８個の非球形粒子の総和計算でバンド幅測定

ノード間のReduction400 力とトルクのリダクションによる総和計算
blockIdx.x=1

総和計算の内訳
1,800

Elapsed time [msec]

A piece ofCAD polygon

謝辞

50
0100101102103
104
Number of particles / object105

本研究の一部は、科学研究費補助金・基盤研究(S) 課題番号26220002, 日本学術振興会特別研究員奨励費・課題番号
14J12004 , 科学技術振興機構CREST「ポストペタスケール高性能計算に資するシステムソフトウェア技術の創出」、革新的ハイパフォーマンス・コンピューティング・インフラから支援を頂いた。また、本研究は東京工業大学・学術国際情報センター
TSUBAME2.5 を利用して計算を行ったものであり、記して謝意を表す。

学際大規模情報基盤共同利用・共同研究拠点 第７回シンポジウム
2015年 7月 9日，10日

THE GRAND HALL （品川）