1979

自己紹介
西浦 泰介（ニシウラ ダイスケ）（国研）海洋研究開発機構（JAMSTEC）
数理科学・先端技術研究分野
役職：技術研究員
専門：粉体工学，計算科学
地球シミュレータGPU

液体と固体粒子の混相流シミュレーション
 干渉沈降 懸濁液の乾燥収縮
 凝集粒子の分散

GPUで個別要素法（DEM）ソフト開発
 SPHによる流体計算との連成
 QDEMによる鉄道軌道解析

1999

同志社大学入学
物質化学工学科

2003

同志社大学
大学院進学

2008

同志社大学卒業
博士(工学)取得
JAMSTEC入社
ポスドク研究員

2011

2013

「京」コンピュータ

大規模並列化プログラムの実用化
 実スケール砂箱シミュレーション
 SPH-DEMによる津波堆積物の輸送
 地盤の遠心載荷シミュレーション

PEZY

福島県に生まれる

現在

JAMSTEC
研究員

JAMSTEC
技術研究員

• 粒子間相互作用力に基づいて，個々の粒子（または計算点）の運動方程式を解きLagrange的に追跡することで，全体の運動を表す方法
開始DEM

初期配置

BD・SD・DPD

物性値

固体粒子
離散体
原子MD相互作用候補ペア作製

SPH

粒子法
流体・固体相互作用力計算
粒子情報更新

MPS

MPI通信
連続体
電磁場

LBM

更新判定

PIC

データ出力
終了

• SPH, DEM等の近接粒子相互作用による粒子法では、粒子が動き回り計算コストの空間的不均質性が時間的に変化する。
それゆえ空間領域分割による分散メモリ環境での効率的なノード間並列は最高ランクの技術的難易度の課題
スパコンを用いた大規模粒子法計算を実現するための動的負荷分散アルゴリズムを開発
粒子法コード並列化の概要

共有メモリ並列・最適化
プロセッサーレベル
ループアンローリングによるベクトル化促進（情報基盤センター並びにNECの協力）

我々の既往研究で最高速な粒子間計算アルゴリズムの採用（Nishiura, Furuichi andSakaguchi, CPC, 2015）

分散メモリ並列・最適化
スライスグリッド領域分割と擬ニュートン法を用いた動的負荷分散（Furuichi and Nishiura, CPC, 2017） 共有メモリ書き込み競合

3

 セルに粒子番号を格納する• セル内の粒子数と粒子番号を記憶• 近傍のセルを探索して相互作用相手粒子を見つける
 粒子に働く力を足し合わせる• 作用反作用の法則を用いる• 相互作用ペア間の力を計算してペア粒子に力を渡す

165・・・ 85

 相互作用相手の粒子数が不均一

6 ロードバランスが悪い

9 10
3

13

 メモリのランダムアクセス
 メモリアクセス速度が低下
 キャッシュが効かない1.セルへの粒子番号登録

2.相互作用力の足し合わせ① Atomic_add関数② 作用反作用の法則を使わない③ 参照表を使う

粒子番号をソート

粒子番号の最小値と最大値を登録

D. Nishiura and H. Sakaguchi, J. Comput. Phys., 230, 1923-1938, (2011)

～Type2～
for i = 1 to Np do
for l = 1 to pair_num[i] doj = Listj[i][l]if distance is less than h then
FV = {forces or relative velocities}
Code1：Type1のペアリスト配列
Code2：並列Scanを使用
Code3：Type2のペアリスト配列
Code4：Atomic_add関数を使用
Code5：作用反作用を使わない

0.20

Code 1
Code 2
Code 3

Code 4
Code 5(a) Xeon CPU

0.15

0.10

0.05

0.00

Interaction

Update particle

Pair list

0.8

Computation time per iteration [s]～Type1～for l = 1 to Nl doi = Listi[l]j = Listj[l]if distance is less than h then
FV = {forces or relative velocities}

Computation time per iteration [s]0.250.6

Code 1
Code 2
Code 3

Code 4
Code 5(b) FX100.40.20.0

Interaction

Update particle

Pair list

D. Nishiura, et. al., Comp. Phys. Comm., 194, 18-32, (2015) 2次元のスライスグリッド法を採用（i）各列の計算コストが一様になるように計算領域を垂直に分割（ii）各列を水平方向に分割し、計算コストが一様な空間要素を作成（iii）各要素にMPIプロセスを割り当てる 特徴１：他の空間充填曲線、グラフ、曲線、階層、3次元分割などの領域分割法と比べて、実装が比較的単純であり、格子の更新に必要な計算コストが少ない。
 特徴２：スライスグリッドは、長方形のジオメトリの制約が性能劣化を引き起こす場合がある。しかし粒子分布がそれほど離散的でなく2次元的な広がりを持つ場合には、性能劣化は限定的。

 動的負荷分散は各MPIプロセスの計算コストを均質にするための前処理 既往研究では“粒子数＝計算負荷”の例が多数。しかし“粒子数＝計算負荷”では、粒子毎に異なる計算コストや境界条件処理、またヘテロな計算機環境に負荷分散を適応できないといった問題がある。
 各MPIプロセスで費やされた“計算経過時間＝計算負荷“とする。
Elapsed Timerank 0𝑡𝑎𝑟𝑖, 0rank 1

Arithmetic

Non blocking MPI
functionssend

Arithmetic

recv

 経過時間𝑡 = 𝑡𝑎𝑟𝑖 +𝑡𝑐𝑜𝑚𝑚 は算術（arithmetic）に用いられた時間
𝑡𝑎𝑟𝑖 とMPI通信に用いられた時間𝑡𝑐𝑜𝑚𝑚 を測定して求める。

id
ol𝑡𝑐𝑜𝑚𝑚, 0
Non blocking MPI
functions 残差（評価）関数𝑓を領域境界𝑥𝐼 𝑛 で定義し𝑓 𝑛 𝑥𝐼 𝑛 = 𝑡𝐼 𝑛 -𝑡𝐼−1 𝑛 ,負荷分散問題をノルム 𝑓 𝑛 𝑥𝐼 𝑛 を最小化する非線形問題と置き換える。
 擬ニュートン法とマルチグリッド法により境界位置を移動
𝑥𝐼 𝑛+1 = 𝑥𝐼 𝑛 + 𝛿𝑥𝐼 ,演算時間基準の負荷バランサー

Line 1: 𝑛 =0: given 𝑥𝐼 0 , 𝐽 ,𝛿𝑥𝐼 = 0
Line 2: Solve SPH kernels and observe 𝑡𝑎𝑟𝑖 and𝑡𝑐𝑜𝑚𝑚 .
Line 3: Compute residual 𝑓 𝑛 𝑥𝐼 𝑛Line 4: if 𝛿𝑥𝐼 ≠ 0 then update 𝐽 𝐼𝐼Line 5: Solve linearized problem 𝐽 𝛿𝑥 = −𝛼𝑓 𝑛Line 6: Update 𝑥𝐼 𝑛+1 = 𝑥𝐼 𝑛 + 𝛿𝑥𝐼Line 7: If 𝛿𝑥𝐼 = 0 then change 𝛼 with line search.
Line 8: 𝑛 = 𝑛 + 1: goto Line 2

粒子の空間分布と種類によって計算負荷不均質性が現れる水塊問題（Water cone test）をサンプル問題とし、並列化性能を測定
計算例：異なる色は各MPIプロセスの担当領域
Water cone test

Wall particles

LyLz

Water particles
Lx• 水粒子の分布形状＝計算負荷の偏り• 水粒子と壁粒子で異なる計算負荷＝計算負荷の偏り
動的負荷分散アルゴリズムの効果を検証

100

 Case1：総粒子数 2.3億個(うち水粒子 1.8億粒子)(Lx, Ly, Lz)=(64, 4, 64)

ESbetter

1/Maximum elapsed time [𝑠 −1 ]並列化率：99.947% on ES ｖｓ 99.961% on K computer

K10 京コンピュータはESと比較して高い並列化性能
 同じノード数の経過時間ではESの方が優れている。

1101001000
Number of Nodes

10000

計算コスト＝経過時間を算術部分と通信部分に分離して解析

Elapsed time [s]ave comm on ES

0.1

0.1

0.01

0.01

ave comm on K

ave on ES
ave on K0.001

0.001
101001000
Number of nodes

10000

10K-computer (K)Earth Simulator (ES)Peak performance

128 GFLOPS

256 GFLOPS

Memory size

16 GB64 GB

Memory bandwidth

64 GB/s

256 GB/s

Interconnect

5GB/s×24GB/s×2100

1000
Number of nodes10000

 ESとKのハードウエア性能差によって、経過時間と並列化率の差が説明可能
 ノード間通信による律速

M. Furuichi and D. Nishiura, Comp. Phys. Comm., 219, 135-148, (2017)• 1ｓｔLoop

ペア粒子間の相互作用力の計算– 境界セル領域の粒子のみ，粒子情報を更新する• 2ndLoop– OpenMPのThread#0が境界セル領域の粒子情報をMPI通信する– 他のThreadは境界セル領域以外の粒子情報を更新する– OpenMPのスケジューリングはDynamicにし，ループの分割サイズを適切に設定する– Thread#0はMPI通信が終わると，粒子情報の更新作業に回る力の足し合わせ計算と粒子情報の更新
通信と計算のオーバーラップを適用
1st loop: 境界領域更新
2nd loop: 残りの更新と通信従来手法（Type 1）

本手法（Type 2）6

本手法（Type 3）

65442

3

123
1

5

オーバーラップのためのループ構造
12halo send34

2nd loop

5send halo
1st loop

1st loop
Type 1612

2nd loop

34

56

send send halo halo
1st loop
Type 2，3

ESにおけるプログラム全体の1Step当りの計算時間（9MPIプロセス，2600万粒子）
Nonoverlap
120.0
演算[ms]
5.5
通信 [ms]Total [ms]125.5
ADB hit [%] 76.24

Type 1 Type 2 Type 3
264.4
264.4
61.07

123.1
123.174.24

120.9
120.9
75.50

Performance（粒子数/計算時間）

Type 3

 ループ分割による性能劣化は小さい 通信コストをほぼ隠蔽することに成功している．
 Type 1：袖領域の粒子へのメモリアクセス効率が悪く，大きく性能が悪化
 Type 2：袖領域のメモリアクセスパターンの効率化により性能が改善
 Type 3：キャッシュ（ADB）効率の向上により，さらに性能が改善
 「京」では，粒子データの再配列作業に必要な整数演算がSIMD化されないため，オーバーヘッドが大きくなりすぎ，性能の改善が見られなかった．
M. Furuichi et. al., SC17, Poster, (submitted)

• 津波シミュレーションにおける動的負荷分散の計算例• 地形データは、計算が進むにつれ必要な所を呼び出す• 粒子数は水部分で約2000万粒子• 色は担当するMPIプロセス
48CPU（MPIプロセス）:384Coreで12時間計算した結果

RadiusRadiusofofinfluence
influence

SPH
Water
Water
Particles
ParticlesW(r-r’,h)W(r-r’,h)

DEM

rr
2h2h
Compact
Compact support
supportofkernelof kernel

粘弾性力と摩擦力と回転トルク

D. Valdez-Balderas et. al., J. Parallel Dist. Comput., 73, 1483-1493, 2013

粘性力と圧力接触力

 2   r  r   aWab pdvap  g   mb  a2  b2   aWab   mb  f a b
  va  vb dt
bbab  a b 
DEMは 摩擦による接線力が有る 回転の自由度が有る
 接線力は過去の値を引き継がないといけない
 粒子間の力をMPI通信する必要が有るvdem
 g   Fcol mdem
t demtI dem  r  Fcol
t

現在の粒子の位置関係からだけでは，法線方向のバネの縮みは解っても，接線方向のバネの縮みは解らない• DEMでは安定ソートで粒子番号を並べ替えないといけない（ペアリストの実装方法に依る）
作用反作用の法則を用い～Type1～る場合，例えば i > j の場
for l = 1 to Nl do 合のみペアリストにする
i = Listi[l]j = Listj[l]～Type2～for i = 1 to Np do
for l = 1 to pair_num[i] do
j = Listj[i][l]4 Type1は検索でListを引き継ぐ
 安定ソートが必要 Type2はListの引き継ぎにソートが必要 ソートの負荷が大きい

768

安定ソート

正しい
6が8から受ける力として引き継ぐ33が5から受けてた力を
5ソート

8

766が7から受ける力として引き継ぐ
符号が逆になる，
或いは検索不能• 砂粒１個レベルに注目した研究は無い• ２次元的な研究や解釈が多い T.Hori(2011), Y.Yamada(2014), T.Dotare(2016)• 粘弾塑性のレオロジーモデルを使った計算が多い
L.Wenk(2013), J.B.Ruh(2014)• 不均質を入れて，付加体の構造を意図的に再現– 摩擦に空間分布を付ける J.B.Ruh(2014)– 海山を沈み込ませる S.Dominguez(1998), P.Zweigel(1998)

破壊がスポット状に発生 → 伝搬 → 応力鎖が交差して固定化→ 応力鎖のアーチ形成
水平方向に伸びる応力鎖のアーチが付加体の形状を決めているDEMの運動方程式

DEM

 sphVdem nvdemmdem
 mdem g   Fss f sf  f ff
tnsphDEM間相互作用項
SPH-DEM相互作用項

SPHSPH-DEM（流体が浸透する粒子群）相互作用
SPH粒子周りの局所的な粒子濃度から粒子-流体間の相互作用力を求める
DEMとSPH間の相対速度

f fs  1  vnn sph dem

1f sf  f fs
1 sph vdem  1    1      1.75  v  v 150

sphdemsph 2r
2r
demdem  3  sph vdem  vsph 1   CD 1.7
8rdem  0.8    0.8

β：粒子層の空隙率を考慮した流体抵抗係数（ Darcy則に基づいた経験則）

SPH-DEM（流体が浸透しない固体連結粒子）相互作用
DEM粒子をSPH粒子として（具体的にはDEM粒子のところにSPH粒子も配置している）密度，粘性項，圧力項を計算し，DEM粒子にそれらの力を与える．
DEM粒子のところに配置したSPH粒子はDEMの座標と速度に強制する．
mgdvg
dtdg

 mg g   Fss   Fsf
  Fss  rgs  Fsf

DEMSPHIg

SPH

SPH

F sf  sphVdem f ff

dt

流体の体積保存（間隙水圧）
状態方程式に空隙率が考慮されているので、空隙率が低下すればTait方程式より圧力が増加する．
t 1 7a t 1t
 a   a  t  mb  va  vb  Wabpa  B   1b
 f • 捨て石（ピンク色粒子）は流体が浸透する相互作用モデル，防波堤（緑色粒子）は流体が浸透しない相互作用モデルを適用～引き波による堤防の転倒～

１GPUで１日計算(独)港湾空港技術研究所 有川太郎上席研究員 提供

計算負荷が粒子の種類（DEM orSPH or 壁粒子）によって違う
計算経過時間基準の負荷バランサーが有効

SPH: １億粒子
DEM: １千万粒子「京」512ノード72時間

Sheet pileΔh : 水位差• SPHやDEMを対象とした粒子法の大規模並列化– 相互作用ペア粒子の探索や力の足し合わせ，DEMでは接線力の引き継ぎなど，近接相互作用系の粒子法には特有の難易度が有る– ２次元スライスグリッド法を採用– 経過時間基準で負荷バランスを均等化– 擬ニュートン法とマルチグリッド法による分割領域サイズの変更– 粒子番号の並び替えとループ分割による通信と演算のオーバーラップ– 空間充填曲線によるキャッシュ効率の改善• 並列化粒子法アプリケーション– 津波の都市部への遡上シミュレーション– 24億粒子による実スケール砂箱シミュレーション– DEMとSPHの連成による混相流シミュレーション– 津波による堤防の倒壊，土砂移動，噴砂