平成 25 年度先端的計算科学研究プロジェクト（ベンチマーク課題）
マルチモーメント移流法を用いた超多次元ブラソフシミュレーションの性能評価
簑島松本敬（海洋研究開発機構）

洋介（千葉大学）、 天野 孝伸（東京大学）

概要
無衝突プラズマの第一原理運動論シミュレーション手法の一つであるブラソフシミュレーションのためのアルゴリズムを開発し、5 次元電磁ブラソフシミュレーションコードに適用した。シミュレーションコードは MPI と OpenMP のハイブリッド並列を施し、集団通信の除去や MPI 通信の隠蔽を行った。高性能演算サーバ Fujitsu PRIMERGY CX400 を用いてシミュレーションコードの性能評価を行い、128 − 16384 並列で並列化率 99.998% 以上、実行効率 14% を得た。

1 はじめに地球を取り巻く宇宙空間は希薄なプラズマで満ちている。宇宙プラズマは 2 体間の衝突が殆ど起こらない無衝突状態のため、流体スケールから粒子スケールまで、様々なスケールの現象が混在している。これらの現象は極めて非線形のため、その定量的理解にコンピュータシミュレーションが果たしてきた役割は非常に大きい。磁気流体力学モデルや第一原理運動論モデルなど、注目する現象のスケールに応じて様々な技法が開発・適用されている。
無衝突プラズマの第一原理運動論シミュレーション手法として、分布関数を位相空間で直接離散化し、その時間発展をオイラー的に追跡するブラソフシミュレーションが提案されている。広く用いられている粒子法（Particle-In-Cell, PIC）と比較すると、粒子法固有の統計ノイズが発生しないため高解像度の計算が期待できる、並列計算による効率化が容易である、計算領域を広く取れる、といった利点がある一方で、より多くの計算機資源を必要とするため、超大規模並列計算機向けの技術である。
ブラソフシミュレーションでは、プラズマ分布関数を実空間と速度空間で更新するが、特に速度空間での数値分散や散逸は非物理的な波動励起やプラズマ加熱をもたらす。これは、磁化プラズマのブラソフシミュレーションにおいて特にシビアな問題となっている。実用的なブラソフシミュ
レーションのためには、数値誤差を極力抑えた高精度解法の開発が必須である。

1

そこで我々はこれまで、情報エントロピー保存の観点から、複数の高次モーメントを同時に更新することで物理量の波形を保持するマルチモーメント移流法を開発し、電磁ブラソフシミュレーションに適用してきた [3, 4]。さらに、安定性と使い勝手を改良した有限体積型マルチモーメント移流法を開発した（論文投稿中）。本プロジェクトでは、有限体積型マルチモーメント移流法を用いて、実空間 2 次元・速度空間 3 次元の 5 次元電磁ブラソフシミュレーションコードを開発した。
シミュレーションコードの性能評価のために高性能演算サーバ Fujitsu PRIMERGY CX400 を利用させていただいたので、その結果について報告する。

2 ブラソフ方程式系
無衝突プラズマの第一原理であるブラソフ方程式系は、ボルツマン方程式の衝突項を無視し、外力にローレンツ力を用いたブラソフ方程式 (1) と、電磁場の振る舞いを記述するマクスウェル方程式 (2)-(5) から成る。


∂ fs
∂ fs
qsv×B
∂ fs+v·+E+·= 0, (s = p, e)∂t∂xmsc∂v∂E= c∇ × B − 4πj,∂t∂B= −c∇ × E,∂t
∇ · E = 4πρ,∇ · B = 0,X Z
ρ=qs fs dv,
s=p,e

j=X

(1)
(2)
(3)
(4)
(5)
(6)

Zqs

(7)

vfs dv.

s=p,e

ここで、E(x), B(x) は電磁場、 j(x) は電流密度、ρ(x) は電荷密度, c は光速, qs は電荷, ms は質量, fs (v, x) は粒子種 s（陽子 p ないしは電子 e）の分布関数である。ブラソフ方程式 (1) を速度空間 v で積分する、もしくはマクスウェル方程式 (2) の発散を取ることで、電荷保存則が導かれる。

∂ρ+ ∇ · j = 0.
∂t(8)ブラソフ方程式 (1) をそのままの形で取り扱うことは数値的に複雑なので、実空間と速度空間の演算子を分離する方法 [1] が良く用いられる。

∂ fs∂ fs+v·= 0,∂t
∂x
qsv×B
∂ fs
∂ fs+E+·= 0.
∂tmsc∂v(9)(10)

このようにすると、式 (9) は実空間で速度一定の線形移流であり、さらに時間刻み ∆t の間で電磁場が一定と近似すると、式 (10) は速度空間での線形移流と剛体回転になるため、線形移流方程式2の高精度数値解法を適用することが出来る。ブラソフシミュレーションでは、これらとマクスウェル方程式を交互に解くことで、システムの時間更新を行う。

3 有限体積型マルチモーメント移流法
マルチモーメント移流法 [3, 4] は、分布関数の速度空間更新 (式 10）に特化した手法である。
式 (10) は磁場周りのジャイロ運動を含んでおり、これによる数値的安定性の低下や急速な数値拡散が問題である。そこで本手法では、情報エントロピー保存の観点から、物理量の高次モーメントを同時に更新することで、安定かつ正確な時間更新を行う。本プロジェクトでは、安定性と使い勝手を改良した有限体積型マルチモーメント移流法（論文投稿中）を用いる。
式 (10) の一般形

∂f∂∂∂+(Fx f ) +(Fy f ) +(Fz f ) = 0,∂t∂vx
∂vy
∂vz

(11)

を速度空間で積分することで、(l + m + n) 次モーメント M l,m,n の時間発展方程式を得る。

M

l,m,n

1
=
l!m!n!

ZZZ
vxl vym vzn f dv,(12) l ZZ ZZZvym vzn
H(l)Fx vxl−1 vym vzn
vx
Fxf dvy dvz +f dv
l!
m!n!
(l − 1)!m!n!
 m ZZ ZZZZvyH(m)Fy vxl vym−1 vzn∂vzn vxl− dvy
Fyf dvz dvx +f dv∂vy m!
n!l!
l!(m − 1)!n!
! ZZZZZZvxl vym
H(n)Fz vxl vym vzn−1vzn∂
Fzf dvx dvy +f dv. (13)− dvz∂vz n!
l!m!
l!m!(n − 1)!

∂ M l,m,n=−∂tZ

∂
dvx
∂vx

（H(l) = 1 for l ≥ 1, 0 for l = 0）. 本手法では、l, m, n ≥ 0, l + m + n ≤ 2 を満たす速度空間でl,m,n

区分化されたモーメント Mi,j,k を式 (13) に従って更新する。物理量を高次精度補間する際に、モーメント量を拘束条件とする。0 次モーメント（質量）の保存に加えて、1 次中心モーメント（中心値）や 2 次モーメント（分散）の保存も保証される点が特徴である。
図 1 は 3 次元移流回転問題（式 10 に相当）を有限体積型マルチモーメント移流法と従来の 3 次精度風上法で解いた結果である。有限体積型マルチモーメント移流法は複数の物理量を更新するので、同程度のメモリ使用量で比較するために、3 次精度風上法では多くのグリッドを用いている。にもかかわらず、3 次精度風上法では急速な数値拡散により分布がぼやけてしまっている。一方で有限体積型マルチモーメント移流法では、500 回転後でも分布が十分保持されている。このことか
ら、有限体積型マルチモーメント移流法は分布関数の速度空間更新に非常に適した手法といえる。

34 ブラソフシミュレーション
本プロジェクトでは、実空間 2 次元・速度空間 3 次元の 5 次元電磁ブラソフシミュレーション
コードを開発した。時間更新する方程式は、ブラソフ方程式の実空間移流 (9)、速度空間移流 (10)、マクスウェル方程式 (2), (3) である。また、マクスウェル方程式の更新では、拘束条件 (4), (5) を保証する必要がある。速度空間移流では有限体積型マルチモーメント移流法を、実空間移流では従来の 3 次精度保存型セミラグランジュ法を用いる。電流 (7) の計算では、電荷保存法 [5] を用いる。
この方法は、電荷保存則 (8) を満たすように電流を計算するので、式 (4) が自動的に満たされる。
マクスウェル方程式の更新では、3 次精度の Lax-Friedrichs 法を用いる。電磁場解析で広く用いられている Finite-Difference Time-Domain （FDTD） 法と比較すると、散逸性誤差のために解像度が低下するが、数値分散が抑えられる。磁場に関しては Constrained Transport （CT）法 [2] のスタガード配置を採用することで、式 (5) を離散化レベルで満たす。以上より、本シミュレーショ
ンコードは、拘束条件としてのポアソン方程式 (4), (5) を解く必要がないように設計されている。
システムの時間更新は次の 2 次精度リープフロッグ法で行う。（1）n ステップ時の分布関数を速度空間で 1 時間刻み分更新する。∗ fs ← n fs （2）次に分布関数を実空間で 1 時間刻み分更新し、(n + 1) ステップ時の分布関数とする。n+1 fs ← ∗ fs（3）(n + 1/2) ステップ時の電流を計算する。
（4）(n + 1/2) ステップ時の電流を用いて、n ステップ時の電磁場を (n + 1) ステップ時に更新する。(n+1 E, n+1 B) ← (n E, n B)
シミュレーションコードは、MPI と OpenMP のハイブリッド並列を施した。ブラソフシミュ
レーションでは電流を計算するために分布関数の速度空間でのリダクション演算が必要なため（式7）、実空間のみを MPI 領域分割し、各 MPI プロセス内で OpenMP スレッドを立ち上げる。これにより、電流計算の際に集団通信を要求しない。さらに、ポアソン方程式 (4), (5) を解く必要がないように設計されていること、用いているスキームが全て陽的解法であることから、本コードは時間更新の間一切の集団通信を要求しない。
その一方で、5 次元位相空間分布関数のうち領域分割されているのは実空間 2 次元だけなので、通信量が多くなる点が問題である。そこで、分布関数を速度空間で更新する際に、次の方法で MPI通信の隠蔽を試みる。（1）実空間で MPI 通信に関わる「のりしろ」に位置している分布関数の速度空間更新を先に行う。（2）実空間で MPI 通信に関わらない位置の分布関数を速度空間で更新しつつ、先に更新した「のりしろ」位置の分布関数の通信を同時に行う。

5 結果
シミュレーションコードの性能評価のため、高性能演算サーバ Fujitsu PRIMERGY CX400 で線形プラズマ波動の計算を行い、強いスケーリングと実行効率の測定を行った。グリッド数は速度空間を 323 とし、実空間を 2562 から 10242 と変化させ、20 ステップの計算に要する時間を測定した。総メモリ使用量は、グリッド数 ×10（モーメント数）×2（粒子種）×8（倍精度） Byte と4見積もられる。実空間と速度空間の境界条件はそれぞれ周期境界、開放境界とした。用いた言語はC/C++、強いスケーリング測定ではインテルコンパイラを、実行効率測定では富士通コンパイラを用いた。コンパイルオプションは-O2,-ipo,-xAVX,-openmp 相当を用いた。
図 2 はグリッド数 323 × 10242 での計算例であり、磁力線方向の電場 Ex の空間分布を示している。ランダウ減衰により 12 プラズマ振動後では短波長モードが少し減衰している。
表 1 及び図 3 は、強いスケーリング測定結果を示している。MPI プロセスあたりの OpenMP
スレッド数は 4 とした*1 。128 コアから 16384 コアまで完璧にスケーリングしていることがわかる。隣接する 2 点の並列数 m, n(> m) での実行時間 Tm , Tn から並列化率 α を求める。

α=1 − n1

Tm − T n .
1Tm − 1 − mTn(14)結果、並列化率は 99.998% 以上となった。システムの総コア数が 23616 コアなので、システムの大多数を使用してもスケーリングが保たれていることになる。この理由として、本コードが集団通信を要求しないことや、MPI 通信の隠蔽が効果的であることが考えられる*2 。
次に、富士通コンパイラを用いて実行効率の測定を行った。この測定ではグリッド数を 323 ×10242 とし、表 1 の上 4 つと同様の計算を行った。実行時間はそれぞれ、5767 秒（128 ノード）、2866 秒（256 ノード）、1537 秒（512 ノード）、739 秒（1024 ノード）となり、インテルコンパイラを用いた場合に比べて 1.5 倍ほど性能が劣化した。実行効率は約 9% であった。よって、インテルコンパイラを用いた場合の実行効率は約 14% 程度と推定される。プロファイラの診断から、シミュレーション全体のうち、有限体積型マルチモーメント移流法による速度空間更新が 8 割以上を占めていることがわかったので、今後は有限体積型マルチモーメント移流法の性能チューニングを行い、シミュレーションコードの高速化を図りたい。

謝辞
高性能演算サーバ Fujitsu PRIMERGY CX400 の利用機会を与えてくださった九州大学 情報基盤研究開発センター及び利用に際してサポートしてくださった深沢圭一郎助教に深く感謝したします。

参考文献[1] C. Z. Cheng and G. Knorr. The integration of the Vlasov equation in configuration space.
Journal of Computational Physics, 22:330–351, November 1976.

*1 並列数を固定して MPI プロセスあたりの OpenMP スレッド数を 8 に増やしたところ、約 1.5 倍の性能劣化が確認された。
*2 PRIMERGY CX400 では MPI 通信を隠蔽しない計算を行わなかったため実際の効果の程は不明だが、別のシステムでそのような小規模計算を行い、通信の隠蔽による実行時間の短縮を確認している。

5

[2] C. R. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flows - A constrainedtransport method. Astrophysical Journal , 332:659–677, September 1988.
[3] T. Minoshima, Y. Matsumoto, and T. Amano. Multi-moment advection scheme for Vlasov
simulations. Journal of Computational Physics, 230:6800–6823, July 2011.
[4] T. Minoshima, Y. Matsumoto, and T. Amano. Multi-moment advection scheme in threedimension for Vlasov simulations of magnetized plasma. Journal of Computational Physics,
236:81–95, March 2013.
[5] T. Umeda, K. Togano, and T. Ogino. Two-dimensional full-electromagnetic Vlasov codewith conservative scheme and its application to magnetic reconnection. Computer Physics
Communications, 180:365–374, March 2009.

6図1非対称ガウス分布の 3 次元移流回転問題。（上段）有限体積型マルチモーメント移流法を用いた結果。左から、初期、100 回転後、500 回転後の分布。グリッド数は 323 . （下段）3 次
精度風上法を用いた結果。左から、初期、100 回転後、500 回転後の分布。グリッド数は 683 .

図 2 電磁ブラソフシミュレーションによる磁力線方向の電場分布。（a）初期、（b）12 プラズマ振動後。

7

1•107

Runtime x cores

8•106

6•106

4•106

2•106

01001000
Number of cores10000

図 3 Fujitsu PRIMERGY CX400 を用いたブラソフシミュレーションコードの強いスケーリング測定。インテルコンパイラ使用。横軸はコア数、縦軸は実行時間 × コア数。

グリッド数

ノード数

MPI プロセス数

コア数

実行時間（秒）323 × 10242

128

512

2048

3891

256

1024

4096

1956

512

2048

8192

999

1024

4096

16384

492

32

128

512

3877

64

256

1024

1959

128

512

2048

955

256

1024

4096

479

8

32

128

3896

16

64

256

1946

32

128

512

963

64

256

1024

483

323 × 5122

323 × 2562

表 1 Fujitsu PRIMERGY CX400 を用いたブラソフシミュレーションコードの強いスケーリング測定。インテルコンパイラ使用。MPI プロセスあたりの OpenMP スレッド数は 4.

8