超並列スカラ計算機における電磁流体コードの最適化と性能評価
荻野瀧樹1、深沢圭一郎2, 3, 4
1.
2.
3.
4.

名古屋大学太陽地球環境研究所
九州大学情報基盤研究開発センター
国際宇宙天気科学・教育センター
CREST, JST

April 26, 2012

ContextIntroduction
太陽地球惑星系科学の紹介
MHD方程式とVlasov方程式

Fujitsu PRIMERGY RX200S6における
電磁流体コードの性能評価
他計算機システムとの比較
並列ベクトル計算機用のMHDコードによる性能評価
Summary1

Introduction | 太陽地球惑星系科学とは2

Introduction | 宇宙天気3

Introduction | 地球磁気圏の形
磁気圏のグローバル構造と境界層の物理をself-consistent解くことが重要

Plasma Mantle
Bow ShockPlasma Sheet
Magnetic
Reconnection

MirrorFire-hose

Inner Boundaryof Plasma Sheet

Rayleigh-Taylor
Entry Layer
Magnetic
ReconnectionFilamentary

Ballooning
Interchange
Low Latitude
Boundary Layer
Kelvin-Helmholtz

Introduction | 地球磁気圏のシミュレーション結果例

5

Introduction | 地球磁気圏のシミュレーション結果例

6

Introduction | 目的
電磁流体コードの性能評価• 惑星磁気圏を解く電磁流体コードを用いて、九州大学のスパーコンピュータシステムでどれだけ速度が出るかを調べる。
• どのようなコードの組み方が最も効率的かを調べる。
• 最終的には上記のコードを利用して、高解像度土星磁気圏シミュレーションを行い、プラズマ対流の渦構造を詳しく調べることを目標とする。
↑午後に発表します7電磁流体コード |宇宙プラズマを取り扱う方程式（1）Vlasov方程式• 無衝突Boltzmann方程式とMaxwell方程式から成るプラズマの振る舞いを最も正確に表現できる方程式系。
速度分布関数f(x, v, t)を考えると、
f
f
f q
 v   (E  v  B ) 0tx mv• これとMaxwell方程式を連立して解く。
• しかし、位置(x, y, z)3次元、速度(VX, VY, VZ)3次元と時間から成る非線形方程式系で、解くことが困難。

8電磁流体コード |宇宙プラズマを取り扱う方程式（2）MHD (Magnetohydrodynamics)方程式• Vlasov方程式のn次モーメント取ることで、求められる。
• 0次(速度空間で積分)、1次(vかけて積分)、2次（v2かけて積分）より、

   ( v )tv
11  ( v   ) v  p  J  B
t


p
 ( v  ) p  p  v
tを得る。これらと磁場の誘導方程式
B   ( v  B) をまとめてMHD方程式という。
t9

MHDシミュレーションの性能評価条件
解法はLeap-frog+Lax Wendroff法の混合手法• 磁場を扱う方程式ではdiv・B=0を保証する手法を用いる必要があり、適用可能手法に限界がある。
• 今までのベンチマークとの比較、また、simpleな手法を用いることで、計算機の性能が出やすいと考えられる。

基本サイズ• 1coreあたりf(x, y, z, 8) = 16387064 (1273×8) grid を設定。
→weak scaling• 計算上この7倍の配列数を使用する。

MPI並列• のりしろは前後一つ。OpenMPは使用せず。
→グリッドを1/2メッシュ幅ずらした2個のシミュレーションボックスを利用し、1ループで2段計算するため、実際は2次精度。
10

MHDコードにおける性能評価手法
3種類の領域分割法
1次元

計算時間（TS）

TS 1  k1n 2 ( p  1)

通信時間（TC）

TC1  k 2 n 2 ( p  1)

2次元

3次元

[Fukazawa et al., 2010]

TS 2  k1n 2 ( p  1)
12

TC 2  2k 2 n 2 ( p  1)

TS 3  k1n 2 ( p  1)
13TC 3  3k 2 n 2 ( p  1)
11

Estimated Time

[Fukazawa et al., 2010]12高効率の並列プログラムの作成
１．先ず、プログラムの構造の確定が最も重要・計算での変数参照を元にプログラム構造を決定
配列変数定義の最小化
計算の明確なブロック化
計算と通信の完全な分離定数と変数は意識して明確に区別
・領域分割法の導入（２、３次元領域分割）
・Dependencyを可能な限り取り除く
計算のメイン部分ではIf文や条件文は絶対に使用不可
定数定義の利用
Mask、List vector、Gather and scatterの最適な使用
・通信では変数の一括変換と一括転送を実施
これを条件文やdo loopを導入せずに実行
13

MPIを用いたMHDコードの構造
計算環境設定2-stepmodifiedLeap-frog法

初期値
前データの read

境界条件
領域分割間データの転送 (send+recv)
First Step の計算
領域分割間データの転送 (send+recv)
Second Step の計算

CPU間のデータ転送を２箇所に集めて一括転送
その場合、条件文とdo loop文を使用せず内部境界の修正
上流境界の更新
サンプリングデータの write
計算終了
14

スカラチューニング | キャッシュチューニング
キャッシュの基本的働き• CPU/core内部に設けられた高速な記憶装置。ここに使用頻度の高いデータを蓄積しておくことにより、低速なメインメモリへのアクセスを減らすことができ、処理を高速化することができる。
• メインメモリにアクセスしたときに、その周辺データをキャッシュに格納する。なので、連続データアクセスが多いプログラムは効率的にキャッシュを利用可能。

性能評価• スカラ機ではキャッシュヒット向上を見込んで、配列の並び替えを行うと性能が上がる機種がある（非x86系）ので、その効果も評価する。
• Type A：f(i, j, k, m)

m: MHD方程式の成分• Type B：f(m, i, j, k)

15

PRIMERGY
RX200S6 (sugoka)における性能評価
16

PRIMERGY RX200 S6
2011年5月に稼働した九州大学の計算機システム
PRIMERGY RX200 S6CPU

型式

6 core Xeon X5670

周波数2.93 GHz (70.32GFlops)

Cache

L1: 64KB/core
L2: 256MB/core
L3: 12MB/CPU

Memory Band幅
B/F値

32/70.32

0.45

Node

Core数

12

メモリ

48GB

Node数

392 (4704core)

理論性能55.13TFlops

富士通製PCクラスタ型計算機
Westmereコア

32GB/s /CPU

SystemNode間通信 InfiniBand 4x QDR (4GB/s)
17

PRIMERGY RX200 S6（2）利用計算機資源• このプロジェクトで利用できた全ノード4704coreで、55.13TFlopsの理論性能

コンパイラオプション• 富士通製とインテル製のFortranコンパイラが使用可能• 富士通製（mpifrt）:-Am –Kfast,fsimple,fuse,prefetch_indirect,prefetch_sequential,SSE4,loop,nomfunc,noparallel,tiling• Intel製（ifort）： MPIが使えない?

18

MHDコードの性能 | 4種類の結果
sugoka（PRIMERGY RX200 S6）での測定
Performance

20,000

3D decomposition Type A
3D decomposition Type B
2D decomposition
1D decomposition

GFlops

16,00012,000
8,000
4,000
00

1,000

2,000
3,000
Number of cores

4,000

5,000
19

MHDコードの性能評価 | 3次元領域分割
Scalability

Performance
20,00017,627

3D decomposition Type A

GFlops

0.96

0.95

0.94
0.92
0.91

3D decomposition Type B

16,0000.99

1 0.970.93
0.9 0.92

12,0000.940.94

0.91

9,023

0.89

12,031

0.87

8,000
4,737
3,506

4,000
0

1,495
767
1,009

0 501

0.8

3D decomposition Type B

3,158
2,265

1,000

3D decomposition Type A

6,113

0.7
2,000

3,000

Number of cores

4,000

5,000

0

1,000

2,000

3,000

4,000

5,000

Number of cores• 実効性能はTypeAで4700コア時に最大17TFlops（実行効率33%）出ている。
• スケーラビリティは低並列時に不安定で、緩やかに下がっており、TypeBのほうが10%近く高い。
20

MHDコードの性能評価 | 他システムとの比較
Core数/CPU数

実効性能[GFlops]

理論性能[GFlops]

実行効率[%]

領域分割CPU種類

SX-964/64

2188

6553

33

2次元ベクトル

SX-8R8/88028228

1次元ベクトル

HA8000 8192/1024

10038

75366

13

3次元A

Opteron

HX600

1024/256

2166

10240

21

3次元AOpteron

FX11024/256

2081

10240

21

3次元B

SPARC64VII

SR16000 1344/672

5375

25267

21

3次元B

POWER6

PG S64704/144

18267

55130

33

3次元AXeonPG S31536/768

2536

18432

14

3次元AXeon

21

MHDコードの性能評価 | CPU当たりの性能
SX-9を1としたときの各種スパコンCPU、coreの性能
@core

@CPU

CPU種類

SX-9

1.00001.0000

ベクトル

SX-8R

0.29250.2925ベクトル

HA80000.03580.1432

Opteron

HX600

0.06180.2474

Opteron

FX1

0.05940.2377SPARC64VII

SR16000

0.11700.2340POWER6

PG S60.1187

0.7122

Xeon

PG S3

0.0483

0.0966

Xeon

実アプリではFX1、SR16000やHX600が4CPUでSX-9の1CPUと同等の性能になっている。
PRIMEGY RX200 S6では1CPUでSX-9のCPUの7割近くの性能

22並列ベクトル計算機用の１次元領域分割MHDコード|
他システムとの比較
１次元分割数

Thread数

格子点数

実効性能[GFlops]

実行効率[%]

CPU種類

PG S6

72

12902*802*802

2259

22.3

Xeon

PG S6

144

6902*802*8022390

23.6

Xeon

PG S6216

4902*802*802

2142

21.2

Xeon

PG S639212902*802*1602

5538

10.0Xeon

FX1644902*802*802

325

12.7

SPARC64VII

HX600

644902*802*802

360

14.1

Opteron

DELL
PowerEdge

240

1902*802*802

27212.9

ADM Opteron• この表での実行性能と実行効率は、データI/Oの処理が入っているために前の２つの表の値と比べた場合、約78％と小さくなっている。
• Xeonはベクトル化機能が有効であることが分かる。4704コアは要検討。
23

今後の課題１．Core数を1-10万個と更に増やした場合もscalabilityは保たれるか？（どうしたらに予測できるか？）
２．通信は本当にボトルネックになるのか？この限界を明確に示すにはどうすればよいのか？
３．Cacheのヒット率を更に上げるにはどうするのか？
４．Flat MPIとHybrid（プロセス並列＋自動並列）はCore数（5000以上）が増えた時、どちらが高効率なのか？
５．スパコンの種類によって高効率のプログラムは異なるのか？
x86系：ベクトル機能が有効、それ以外：cacheヒットが有効？
６．最善のテストは最大数のCPUを用いて、粒度を意識して絶対効率を実測すること
24Summary
MHDシミュレーションの性能評価をおこなった
 PRIMERGY RX200S6では3次元領域分割TypeAが最も性能が高く、キャッシュヒットを狙ったType Bは一番性能が出なかった。
 最高で4700コア時に17TFlops、実行効率33%が出ている。
 他の計算システムと比較すると、スカラではx86系は3次元領域分割TypeA、それ以外では3次元領域分割TypeBが最適であった。
 おそらくSSEの効果によりベクトル化が効いており、キャッシュヒット率を上げることよりも効果が出ている。
 ほぼSX-9の1CPUにXeonが追いついた。
 並列ベクトル機用に開発したほとんどの１次元領域分割コードが864コアまでは実行効率21-24％の高性能を得た。4704コアでは10％程度。
 プログラム構造を確定することが重要（通信量の最小化と一括転送）
 プログラムサイズの最小化（変数の１次元化と領域分割の両立性？）
25