平成 24 年度先端的計算科学研究プロジェクト
第一原理無衝突プラズマシミュレーション用
5 次元ブラソフコードの性能評価
梅田隆行（名古屋大学太陽地球環境研究所）
深沢圭一郎（九州大学情報基盤研究開発センター）

1. 研究の目的と意義
太陽から地球に至るジオスペース環境の変動を理解することは、人類の活動が宇宙へと拡大しつつある今日、極めて重要な課題である。本研究の大きな目的は、ジオスペースで生起するプラズマの非線形現象を解明し、宇宙環境変造の因果関係を理解するのと共に、数値宇宙天気予報に適用することである。人類の活動に影響を与えるジオスペースの変動現象としては、突発的な磁気嵐やオーロラの爆発現象、放射線高エネルギー粒子生成、高エネルギー粒子線による人工衛星の誤作動などが挙げられる。これらの現象は、電磁気圏プラズマのグローバルな対流循環、境界層で生起するメソスケールでの突発的な不安定性（平衡状態の破れ）及び、電子・イオンが粒子として振舞うミクロスケール現象（粒子加速・加熱）が複雑に結びついており、マルチスケール結合過程である。これらの多様な時空間スケールの非線形現象を解明するために、様々なグループによって、グローバル現象を扱う磁気流体力学（MHD）／多流体モデル、ミクロ現象を扱う運動論（粒子／ブラソフ）モデル及び、両者の中間（メソ）スケール現象を扱う流体と運動論のハイブリッドモデルが、それぞれ個別の現象に対して適用されてきた。本研究ではその中でも特に、（現行の MHD・粒子モデル、次世代の多流体・ハイブリッドモデルに対して）次々世代の宇宙プラズマシミュレーション手法であるブラソフモデルに着目しており、粒子運動論スケールからメソ・マクロスケールまでをシームレスに扱う大規模シミュレーションにより、無衝突宇宙プラズマのマルチスケール結合過程の解明を目指している。本プロジェクトでは特に、最新のスーパーコンピュータの能力を最大限に活用できるように、ブラソフコードの並列化・最適化を行うことを目的とする。
ブラソフモデルは、無衝突宇宙プラズマの運動論シミュレーション手法の 1 つであり、位置－速度位相空間に定義されたプラズマ粒子の分布関数の時間発展をブラソフ方程式により直接解き進めている。もう１つの手法である粒子モデルでは、プラズマ粒子であるイオンや電子などの個々の荷電粒子の運動をニュートン－ローレンツ方程式により解き進めている。格子（Cell）上に定義された電磁場中を粒子が動きまわることから、PIC（Particle-In-Cell）法と呼ばれている。宇宙空間に存在する膨大な数の荷電粒子を有限の計算機資源で扱うことは不可能であるため、ある程度まとまった数の荷電粒子の集団を 1 つの“超”粒子として扱う。PIC 法はその数値解法の完成度が高く、プラズマ科学分野では広く用いられている。しかし、超粒子として扱うことにより熱雑音が大きくなり、電荷密度や電流密度などの荷電粒子の運動に起因する場の量を格子上に割り振る際に生じる高波数モードが数値誤差として蓄積することが欠点として挙げられる。一方、格子上に定義された分布関数は熱雑音を持たないため、PIC 法の欠点を克服できる。また、メソ・マクロスケールの現象においてプラズマ密度の空間的非一様性は本質的であり、PIC 法において並列化の際に負荷のバランス（各プロセス内の粒子数の均一性）を保つためには特殊なデータ分割が必要になるが、格子上に定義された分布関数を用いるブラソフモデルにおいては、流体モデルと同様に領域分割による並列計算も容易である。
1しかし、ブラソフ方程式は実空間 3 次元及び速度空間 3 次元の計 6 次元を扱う方程式であり、コンピュータで解くには膨大なリソースを必要とする。このため、その手法の開発はあまり進んでいない。実際、ここ数年の HPC プロジェクトによる計算機環境の飛躍的な向上によって手法の開発が進み、実空間 2 次元及び速度空間 3 次元の 5 次元テストシミュレーションがようやく実行できるようになった段階である。本研究グループではこれまで、磁気リコネクション（テアリング不安定性）やケルビン－ヘルムホルツ不安定性などのメソスケール現象のみならず、イオンジ
ャイロ半径スケールの半径を持つ非磁化小天体と太陽風との相互作用などのグローバルシミュレーションにも世界で唯一成功している[Umeda et al., 2010a,b; Umeda et al., 2011; Umeda, 2012;
Umeda et al., 2012a]。
本プロジェクトで用いるブラソフコードは、東京大学・名古屋大学・九州大学等のスーパーコンピュータシステムにおいて全ノード計測を行っており、HPCI「京」若手人材育成利用にも採択されている。これまでの弱いスケーリングテストの実績では、Infiniband 規格以上で全ノードフルバイセクションバンド幅構成の高速ネットワーク環境において、非常に高いスケーリング性能を達成しており、FX10 の Tofu ネットワークでは 90%以上のスケーラビリティが期待できる。一方
CX400 はネットワーク構成がブロック化されており（図 1 参照）、回線の細いネットワークが実効性能にどのように影響するのかが興味深い点である。また大規模ベンチマークの意義としては、コードの性能評価のみならず超並列計算機システム自体の性能向上が挙げられる。本研究グルー
プはこれまでに様々なシステムにおいて OS やネットワークパラメータの設定不備や障害検知に貢献してきたことを追記する。

図１：九州大学 CX400 のネットワーク構成。256 ノードをフルバイセクションバンド幅（FBB）接続した構成が 5 セットと、192 ノードを FBB 接続した 1 セットが、総バンド幅 2.6TB/sec で接続されている。

2. 計算手法の概要
無衝突プラズマの振る舞いは、以下のブラソフ（無衝突ボルツマン）方程式によって記述される。

f s  f s q s    f sv  ( E  v  B)   0(1)tr msv
   
ここで E 、 B 、 r 及び v はそれぞれ電場、磁場、位置、速度を表す。また、 f s (r , v , t ) は位置－速度位相空間におけるプラズマ粒子の分布関数であり、 s はイオンや電子など粒子種を示す。 q s と ms はそれぞれ電荷と質量を表す。
プラズマ粒子の分布関数は、式(1)で示される通り電磁場によって変形する。電磁場の時空間発展は以下のマックスウェル方程式によって記述される。
2 1 E  B  0 J  2(2.1)
tc


B
 E  (2.2)
t (2.3)
E 0
(2.4)
 B  0
ここで J は電流密度、 は電荷密度、 0 は真空中の透磁率、 0 は真空中の誘電率、c は光速を示す。
ブラソフ方程式(1)を速度空間で積分すると、以下の電荷保存則が得られる。

 J  0t(3 )

電磁場を変化させるための電流密度 J はプラズマの運動によって生じるが、これはブラソフ方程式(1)の第二項にある実空間の流束 v f s を速度空間で積分することによって求まり、電流密度 J が電荷保存則(3)を満足する限り、ポアソン方程式(2.3)は自動的に満たされる。以上の方程式は、ブラソフコードにおいて解いているプラズマ粒子の運動論方程式であり、無衝突プラズマの第一原理と呼ぶ。
ブラソフ方程式は 4 次元以上の「超多次元」を扱う方程式であり、そのままの形で多次元数値積分を行うのは非常に困難であるため、演算子分離（operator splitting）法が古くから用いられてきた[e.g.,Cheng & Knorr, 1976]。過去の研究では、各次元それぞれを 1 次元移流方程式に分解する方法が採用されていたが、本研究では、以下のように実空間移流、速度空間移流、速度空間回転の 3 つの物理的な演算子に分離する手法を開発した[Umeda et al., 2009]。

f s  f sv   0
tr
f s qs  f sE  0t ms v
f s qs   f s
(v  B )   0
t ms
v

(4.1)
(4.2)
(4.3)

この演算子分離は、粒子（PIC: Particle-In-Cell）法においてニュートン－ローレンツ式（荷電粒子の運動方程式）を時間 2 次精度で解く手法として広く用いられている Boris アルゴリズム[Boris, 1970]に基づいている。

式(4.1)及び(4.2)は多次元の移流方程式であり、 v 及び E がそれぞれ r 及び v に依存しない線形方程式である。ブラソフ方程式の数値積分では、古くからセミラグランジュ法が採用されてきた。即ち、格子点上に与えられている分布関数の値を更新する際に、次の時刻にその格子点上に来るべき分布関数の値の座標をバックトレースし、隣接する格子点より値を内挿する。2000 年代前半までは、スプライン補間や CIP（Constraint Interpolation Profile）法などの高次補間法を用いたセミラグランジュ法が用いられた、実空間 1 次元、速度空間 1 次元の、2 次元位相空間シミュレーションが主流であった。本
研究をはじめとする近年の研究では保存型セミラグランジュ法を採用し、超多次元シミュレーションに成功している。これは、保存則であるブラソフ方程式において、プラズマの数密度の保存を数値的に保証し、数値誤差を軽減しようとするものである。PIC シミュレーションでは、電荷保存則(3)を厳密に満たすことによって数値誤差の軽減を行っており、ブラソフシミュレーションでは保存型解法を用いることにより、式(3)が自動的に満たされるため、保存型解法を採用するメリットは大きい。また
式(4.1)及び(4.2)の多次元の線形移流方程式を解く上で、計算精度を上げることは非常に重要であるが、2 次精度以上の手法を用いた際に現れる数値振動を抑制するためのリミッタを保存型解法に導入するのは比較的容易である。

3

本研究では、演算子分離による数値拡散を抑制するために、多次元の線形移流方程式に対する演算子非分離（unspliting）法を新たに開発している[Umeda et al., 2009]。また本研究では、無振動性及び正値性を保証するリミッタを新たに開発し、数値振動の抑制を行っている[Umeda, 2008; Umeda et al.,2012b]。ここで無振動スキームとは、ある区間において新たな極値（極大、極小）を生じず、既に存在する極地を（できるだけ）減衰させないスキームである。本研究で行ったベンチマークテストでは、開発途上の 5 次多項式を用いた補間をテスト的に採用している。
式(4.3)は荷電粒子の速度が磁力線により運動エネルギーを保ったまま変化する回転方程式を表す。
直交座標系における回転方程式は剛体回転問題と等価であり、線形移流問題と同様に、数値計算において最も基本的であるが、計算精度が重要となる問題である。本研究で採用している back-substitution法[Schmitz & Grauer, 2006]では、Boris アルゴリズム[Boris, 1970]に基づいて速度空間での粒子の軌道をバックトレースし、vx、 vy、 vz 方向それぞれの演算子を分離して回転運動を解いている。剛体回転問題では、系の外側（速度空間では速度が速い成分）になればなるほど移動量が大きくなり、クーラン条件の影響を受けやすくなる点に注意が必要であり、今後、陰解法や演算子非分離法の開発が必要である。
以上のように、ブラソフ方程式の数値解法は未だ発展途上である。この大きな原因は、ブラソフコードで扱う次元が多いためであり、開発やデバッグにために大容量の共有メモリ環境が必要となるからである。
一方、マックスウェル方程式(2.1)及び(2.2)は、FDTD（Finite Difference Time Domain）法と呼ばれる電磁場解析法を用いて解く。FDTD 法では、Yee 格子[Yee, 1966]と呼ばれる staggered 格子を用いており、式(2.4)が自動的に満たされるように物理量が配置されている。また leap-frog アルゴリズムに基づいて電場と磁場の時刻も半タイムステップずらしており、空間精度と時間精度は共に 2 次である。
ブラソフシミュレーションでは非常に多くのメモリを必要とするため、並列計算が必須となる。ブラソフコードで使用する物理量は全て格子点上で与えられており、並列化においては領域分割法が有効である。図 2 は実空間 2 次元及び速度空間 3 次元を使用するブラソフコードにおける並列化の概念を示す。我々の目は 4 次元以上の空間を認識できないが、2 次元実空間の各格子上に 3 次元速度空間（速度分布関数）が定義されていると考えると分かりやすい。本研究では図２のように実空間（x－y平面）においてのみ領域分割を行い、速度空間の領域分割は行わない[Umeda et al., 2012a]。これは、電荷密度や電流密度などのモーメント量を計算する際に必要な速度空間の積分において、各実空間でのリダクション処理を行わないようにするためである。またオプションとして、プロセス毎に OpenMPによるスレッド並列化も許可している。

図２：5 次元ブラソフコードにおける空間領域分割[Umeda et al., 2012a]。
43. 性能測定
本プロジェクトでは、FX10（768 ノード）及び CX400（1476 ノード）を用いたブラソフコードの性能評価を行った。それぞれのシステムにおいて用いたコンパイラオプションは以下のとおりである。
FX10: mpifrt –Kfast,visimpact,preex,openmp,simd=2 –x250-W0,-zmpa=Ailfunc-for-floor
CX400(富士通コンパイラ): mpifrt –Kfast,openmp,parallel –x250
CX400(インテルコンパイラ): mpiifort –ip –ipo –O3 –xAVX -opemmp
こ こ で 、 -x250 は ユ ー ザ 定 義 関 数 の イ ン ラ イ ン 化 用 の オ プ シ ョ ン で あ る 。 -Ksimd=2は SIMD 化 強 化 用 の オ プ シ ョ ン で あ る が 、必 ず -Kvisimpact よ り も 後 ろ に つ け る 必 要 があ る 。ま た 、-W0 以 下 は floor 関 数 強 化 用 の 非 公 式 オ プ シ ョ ン で あ る が 、次 の コ ン パ イ ラ
バ ー ジ ョ ン で は オ プ シ ョ ン な し で 関 数 の 強 化 が 行 わ れ る 予 定 で あ る 。 な お 、 FX10 とCX400 で は コ ー デ ィ ン グ の 異 な る コ ー ド を 用 い て 性 能 測 定 行 っ た 。 こ れ は 、 東 京 大 学 の大 規 模 HPC チ ャ レ ン ジ に よ っ て 、 FX10 と CX400 で は チ ュ ー ニ ン グ 方 法 が 異 な る こ と が事 前 に 分 か っ て お り 、特 に FX10 で は 、整 数 の 組 み 込 み 関 数 を 用 い ず に 倍 精 度 に 型 変 換 して 組 み 込 み 関 数 を 呼 び 出 し た ほ う が 高 速 だ か ら で あ る 。 [梅 田 &深 沢 , 2013]。
まず参考値として、単コアを用いた場合の各システム／コンパイラの性能を表 1 に示す。CX400 とインテルコンパイラの組み合わせは、CX400 と富士通コンパイラの組み合わせの約 1.5 倍、FX10 の約 2.2 倍高速であることが分かる。なお、FX10 の経過時間に対する実行性能と実行効率が CX400 よりも高い理由は、上記の型変換により FX10 用コードの浮動小数点演算数が若干多いためである。
本プロジェクトの全ノード計測では、“コアあたり”の格子数を実空間では 40×20、速度空間では30×30×30 に設定する。これはコアあたりのメモリ使用量を約 1GB に固定してノード数を増やしていく弱いスケーリングに該当する。ブラソフモデルは４次元以上の超次元を扱うためにメモリ使用量が非常に多いため、速度空間の格子点を 303－603 に固定してノードのメモリ使用量 1-4GB に設定しつつ、使用ノード数を増やして計算領域を拡張していくのが実際の超並列計算機の利用方法である。
1 コアあたり 1GB にメモリ使用量を固定したときの、5 次元ブラソフコードの弱いスケーリング性能を図 3 及び図 4 に示す。また参考値として、FX1 での計測結果も示す。図 3 はコア数に対する実効速度（GFlops）を表し、図 4 はコア数に対する実効並列効率（スケーラビリティ：N コアの実効効率を 1 コアの実効効率で割った量）を表す。FX10 においてはスレッド数４のハイブリッド並列が最速であり、CX400 と FX1 ではフラット MPI 並列が最速であったため、主にそれらの結果を用いているが、CX400 の一部のデータはスレッド数４のハイブリッド並列のものである。
表１：ブラソフコードの単コア性能。
格子点数 40x20x30x30x30/core、時間ステップ 5 での結果。

システム（コンパイラ）

経過時間

実効性能

実効効率

FX10

200.0sec

2.54GFlops

17.2%

CX400（富士通）

153.2sec

3.11GFlops

14.4%

CX400（インテル）

103.6sec

4.59GFlops

21.2%

5

FX10 においては、スムーズに 768 ノード（12,288 コア）までの計測が行えた。FX10 の結果は東 京
大 学 の 大 規 模 HPC チ ャ レ ン ジ [梅 田 &深 沢 , 2013]と ほ ぼ 同 じ で あ り 、九 大 と 東 大 の シ ス テ
ム に お い て 環 境 設 定 に 差 は ほ と ん ど 無 い と 言 え る 。 FX1 で は 、全ノードを用いて計測した場
合に極端に性能が劣化する現象が見られた（名大の 3072 コア：図 4 の FX1 の実線データ、JAXA の12,032 コア：図 4 の FX1 の点線データ）が、FX10 では並列性能に多少のブレがあるものの、全ノードを用いた場合でも FX1 で見られたような性能劣化は全く起こらなかった。
また 256 コアまでは 100%を超えており、これはスレッド化によって単コアよりも性能が向上しているためである。FX10 上では 1000 コア以上でも 97%以上の実効並列化率でスケールしており、Tofu ネットワークが非常に優れた性能であることを示している。
一方 CX400 においては、最大で 1024 ノード（16,384 コア）までの計測しか行えず、全ノードでの計測は実現できなかった。1024 ノード以上を用いた場合、ジョブからの出力がないままタイムアウトする現象が多発した。図 3 及び 4 の結果は、唯一成功したスレッド数４のハイブリッド並列の結果を示している。なお、この現象に原因は未だ不明である。
また 16 ノードから 256 ノードの計測において、性能が激しく劣化する現象（経過時間にして約 20秒の遅延、実効効率にして約 3%の低下）が見られた。使用ノード数を増やすとこの現象の発生頻度も増えるため、図 1 に示すネットワーク構成に原因があると考えられる。即ち、使用ノード数が多い場合に、ジョブが１つの FBB グループ内に割り当てられるとは限らず、FBB グループを跨ぐ遅いネットワークを利用した場合に性能劣化が生じる。

図 3：1GB/core 使用時の 5 次元ブラソフコードの

図 4：1GB/core 使用時の 5 次元ブラソフコードの弱いスケーリング性能(1)。コア数に対する実効速弱いスケーリング性能(2)。コア数に対する実効並度を示している。

列効率を示している。

4. まとめブラソフコードの全ノード計測において、FX10 の 768 ノード（12,288 コア）を用いた場合の最高性能は 57.5TFlops（実効性能 16.7%）であり、CX400 の 1024 ノード（16,384 コア）を用いた場合の最高性能は 57.2TFlops（実効性能 16.2%）であった。Sandy Bridge プロセッサの単コアでは 20%以上の実効性能を達成しており、並列性能の劣化はネットワーク構成に依るところが大きい。また CX400は FX10 に比べて障害発生率が高いように思われ、1024 ノード以上の計測が実現できなかったのはノードまたはネットワークの障害によるものと推測している。
6

参 考 文献

Boris, J. P.: Relativistic plasma simulation-optimization of a hybrid code, in: Proc. Fourth Int. Conf. Num. Sim.
Plasmas, ed. by J. P. Boris and R. A. Shanny, pp.3—67, Naval Research Laboratory, Washington D. C. (Nov.
1970).
Cheng, C. Z., Knorr, G. : The integration of the Vlasov equation in configuration space, J. Comput. Phys., Vol.22,No.3, 330—351 (1976).
Harris, E. G., On a plasma sheath separating regions of oppositely directed magnetic fields, Nuovo Cimento,Vol.23, No.1, 115—121 (1962).
Schmitz, H., R. Grauer, R.: Comparison of time splitting and backsubstitution methods for integrating Vlasov's
equation with magnetic fields, Comput.Phys. Commun., Vol.175, No.2, 86—92 (2006).
Umeda, T.: A conservative and non-oscillatory scheme for Vlasov code simulations, Earth Planets Space, Vol.60,No.7, 773—779 (2008).
Umeda, T.: Effect of ion cyclotron motion on the structure of wakes: A Vlasov simulation, Earth Planets Space,Vol.64, No.2, 231—236 (2012).
梅田隆行, 深沢圭一郎: 5 次元ブラソフコードによる無衝突プラズマシミュレーションの性能評価, 東京
大学情報基盤センター スーパーコンピューティングニュース, Vol.15, No.1 47—55 (2013).
Umeda, T., Fukazawa, K., Nariyuki, Y., Ogino, T.: A scalable full electromagnetic Vlasov solver for cross-scale
coupling in space plasma, IEEE Trans. Plasma Sci., Vol.40, No.5, 1421—1428 (2012a).
Umeda, T., Kimura, T., Togano, K., Fukazawa, K., Matsumoto, Y., Miyoshi, T., Terada, N., Nakamura, T. K. M.,
Ogino, T.: Vlasov simulation of the interaction between the solar wind and a dielectric body, Phys. Plasmas,Vol.18, No.1, 012908(7pp.) (2011).
Umeda, T., Miwa, J., Matsumoto, Y., Nakamura, T. K. M., Togano, K., Fukazawa, K., Shinohara, I.: Fullelectromagnetic Vlasov code simulation of the Kelvin-Helmholtz instability, Phys. Plasmas, Vol.17, No.5,052311(10pp.) (2010a).
Umeda, T., Nariyuki, Y., Kariya, D.: A non-oscillatory and conservative semi-Lagrangian scheme with
fourth-degree polynomial interpolation for solving the Vlasov equation, Comput. Phys. Commun., Vol.183,No.5, 1094—1100 (2012b).
Umeda, T., Togano, K., Ogino, T.: Two-dimensional full-electromagnetic Vlasov code with conservative scheme
and its application to magnetic reconnection, Comput. Phys. Commun., Vol.180, No.3, 365—374 (2009).
Umeda, T., Togano, K., Ogino, T.: Structures of diffusion regions in collisionless magnetic reconnection, Phys.
Plasmas, Vol.17, No.5, 052103(6pp.) (2010b).
Yee, K. S., Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic
media, IEEE Trans. Antenn. Propagat., Nol.AP-14, No.3, 302—307 (1966).

7