九州大学先端的計算科学研究プロジェクト 平成 27 年度成果報告書 2016 年 4 月

月周辺プラズマ環境の高並列プラズマ粒子シミュレーション（課題名：月面近傍プラズマ・ダスト環境の高並列プラズマ粒子シミュレーション）
三宅洋平 1、西野真木 2、木倉佳祐 1、臼井英之 1
1. 神戸大学大学院システム情報学研究科／計算科学教育センター
2. 名古屋大学宇宙地球環境研究所
概要 九州大学 HA8000-tc/HT210 を用いた大規模高並列プラズマ粒子シミュレーションにより、月などの非磁化「小型天体」表面とプラズマ間の相互作用を解明する。具体的には、月面に相当する固体表面境界をプラズマ空間内部に設定できる衛星プラズマ環境シミュレーターEMSES を用いて、月縦孔地形や、月面磁気異常と太陽風プラズマの間の相互作用を解析する。これに加えて、独自の動的負荷分散手法 OhHelp を適用した
EMSES コードに関して、SIMD-aware なプラズマ粒子計算実装法を探求することで、メニーコアプロセッサ上での高効率計算が可能なシミュレーションコードの開発を目指す。これらの取り組みを通してプラズマ粒子シミュレーション計算のさらなる高度化に向けた技術的知見を得ることを目標とする。

1. 研究の目的と意義間的な粗密による負荷不均衡が常に問題となるが、本研究の目的は、九州大学 HITACHI HA8000-EMSES にはこれを解決するための動的負荷分散tc/HT210 を用いた大規模高並列プラズマ粒子シミ技法 OhHelp が実装されており、これまでに「京」ュレーションにより、人工衛星や月などの「小型天コンピューター、九州大学、京都大学、名古屋大学体」と太陽風プラズマ間の相互作用を解明するこの大規模並列計算機システム上で高効率の並列計とである。人工衛星や月面などの非プラズマ物体算を行ってきた実績がある。当該シミュレーショ周辺のプラズマは境界領域プラズマと呼ばれ、シンコードのさらなる性能向上のために、共有メモースや光電子層、電位バリア形成といった著しいリ並列化や SIMD-aware な計算アルゴリズムの開空間非一様性を有するため、その全容はいまだに発などアーキテクチャ階層を意識した最適化が必理解されていない。本研究では、近年「かぐや」衛要である。本レポートではプラズマ粒子計算のカ星観測で発見され、将来の着陸探査対象や月面基ーネル演算である電流計算ルーチンの SIMD 最適地建設候補地として注目が集まっている月面上の化の取り組みも紹介する。

縦孔に着目し、複雑な形状を有する月固体表面と2. 計算手法の概要

太陽風プラズマとの電気力学的相互作用を自己無動着に解き進め、縦孔周辺静電環境の 3 次元的描

衛星プラズマ環境シミュレーターEMSES は 3像を明らかにする。また月面の磁気異常と太陽風次元電磁粒子シミュレーション法を基本とし、こプラズマの相互作用によって、上空に形成されるれに固体表面の内部境界の扱いを導入したものと小型磁気圏境界面付近での電子のダイナミクスをなっている[1]。計算空間は格子によって離散化さ解き明かす。

れており、電磁場と電流／電荷密度場はオイラー

これに加えて、数値シミュレーション手法とし量として格子点上に定義される。一方、プラズマてのプラズマ粒子（Particle-In-Cell）計算の高度化

分布関数は連続的な位置座標を取りうる代表粒子にも取り組む。次節に紹介するように、本課題で用によって記述される。そしてこれらの粒子と場のいる衛星プラズマ環境シミュレーターEMSES は

量は Particle-In-Cell と呼ばれる手法により互い領域分割に基づく Flat-MPI 並列化が施されておに結び付けられる。基礎方程式としては、電磁場り、分散メモリ型並列計算機による大規模計算が発展を記述する Maxwell 式と荷電粒子ダイナミ可能である。また単純な領域分割法では粒子の空クスを記述する運動方程式が用いられ、これらを1

九州大学先端的計算科学研究プロジェクト 平成 27 年度成果報告書 2016 年 4 月

図 1:月縦孔近傍の帯電ダスト挙動のテスト粒子解析シミュレーション．
相互に解き進めることによりプラズマと電磁場間縦孔構造や地下空洞を直接探査する UZUME 計の相互作用を自己無動着に再現する。

画が立案されている中、将来の着陸探査への影響

EMSES の分散メモリ並列は、均等領域分割にを知る上でも縦孔周辺プラズマおよび帯電ダスト基づいており、各 MPI プロセスは分割された各部環境を事前評価することは重要である。そこで分領域とそこに含まれる粒子に関する処理を担当

OhHelp 適用済の全粒子プラズマシミュレーショする。しかしこの単純な領域分割法のみでは粒子

ンコード EMSES を、月面昼側の月縦孔周辺プラの空間的な粗密による負荷不均衡が問題となる。

ズマ環境の定量評価に適用した。

そこで OhHelp では、一つのプロセスを除く全て

本解析では、Marius の丘に発見されている縦孔のプロセスが本来の担当（1 次担当領域）とは別ののデータを参考に、そのサイズを直径 50 m、深さ部分領域を一つだけ担当し（2 次担当領域）、その45 m と定め、その上空も含めて、おおよそ 200領域に含まれる粒子の一部分について電磁場とのm×200 m×1 km の空間を計算領域とし、最大で相互作用計算を行うことで、この問題を解決してHA8000-tc の 128 ノードを利用した計算により、いる[2]。これに加えて、OhHelp 特有の粒子移送①縦孔近傍のプラズマ分布および電気環境、②そパターンを考慮して低コストな MPI 通信方式をの中での帯電ダストの挙動を明らかにした。

考案し、さらなる最適化を行った[3]。

一般的にプラズマに固体表面が接している状態
では、プラズマ荷電粒子が表面に捕捉されること

3. シミュレーション成果

により、帯電が生じる。縦孔内部では、日向と日(1) 月縦孔周辺プラズマ・ダスト環境の大規模全陰で帯電値が大きく異なり、大きな電位差（40 V

粒子シミュレーション程度）が生じている。これは、日照面では先述し「かぐや」衛星による地形カメラ観測により、た光電子放出により、正に帯電するのに対し、日月面に特徴的な縦孔構造が発見されている。この陰面からは光電子は放出されず、周辺を浮遊する
2

九州大学先端的計算科学研究プロジェクト 平成 27 年度成果報告書 2016 年 4 月

電子の一部が表面に付着することで負に帯電するな条件下では、太陽風プラズマと磁気異常磁場間ためである[4]。の相互作用により小型磁気圏が形成される。また

この帯電現象は月面のみならず、月表層に堆積磁気圏境界層における太陽風電子とイオンの磁場している微小な塵であるダストに対しても発生すに対する応答差に起因する静電界が磁気圏形成にる。帯電したダストは、月縦孔周辺の特異な電気重要な役割を果たす。

環境のもとで浮遊し、特徴的なダイナミクスを示

本シミュレーションでは、特に磁気圏境界層にすことが予想されるが、こうした帯電ダストが月おける電子ダイナミクスによって決定づけられる面探査用のローバーに付着することにより、機器電流構造に着目して、詳細解析を実施した。シミの故障などの障害を引き起こす恐れがある。そこュレーション結果から、境界層低緯度領域においでダストの挙動を明らかにすることを目的に、テては朝側から夕方側に向かう電子ドリフト運動がスト粒子シミュレーション解析を実施した。

顕著であるが、中高緯度領域においては逆に夕方

図 1 に縦孔周辺でのダスト分布のスナップショ側から朝側に向かう電子の流れが主に見えた。こットを示す。初期状態ではダストは月面上に堆積のような二層の電子ドリフト運動の構造により磁しているという仮定の下での計算であるが、静電気圏層間側、すなわち磁気異常上空の南北両半球気力の作用により、ダストが浮遊していることがにおいて電子電流の渦的構造が見られた。特に赤図から確認できる。ダストの持つ電荷は日向では道面における電子電流には、電子－イオン間の荷正、日陰では負であるが、一方で日向および日陰

電分離によって発生する電場と磁気異常磁場によの月面からはそれぞれ上向きおよび下向きの電場

る E×B ドリフトと、磁気圏境界層付近での電子が発生しているため、ダストが浮かび上がる向き

圧力の急速な減少に伴う反磁性ドリフト電流が大に静電気力が作用する。図からは判断できないが、きな役割を果たしていることが明らかとなった。

この静電気力と下向きの重力との釣り合いにより、 これらの結果は小型磁気圏の構造、ひいてはそのダストが上下に振動する運動が確認されている。

周辺でのイオンのダイナミクスにも大きく影響を現在はダストの直径と質量をある一定値に固定し及ぼすという点で、月プラズマ科学上、きわめてた計算であるが、今後はダストの直径に分布があ重要な成果ということができる[5]。

ることを想定した計算を実施し、平均的なダスト4. 計算コードの高度化および性能評価密度分布を算出することを目指す。

OhHelp 適用済プラズマ粒子コードのさらなる(2) 月磁気異常‐太陽風プラズマ相互作用の大規高度化に向けた検討として、計算負荷の大部分を模全粒子シミュレーション

占める電流計算に着目し、SIMD ベクトル化を促

上述の解析は、月面上空およそ 1 km までの描進する実装法を探求した。Xeon Phi など近年、普像であるが、より上空では太陽風プラズマと月に及が進むメニーコアプロセッサにおける SIMD 演

局所的に存在する磁気異常の間の電磁的相互作用算機構は、NEC SX シリーズに代表されるようなが電気環境の決定に大きな役割を果たす。そこで、ベクトル計算機におけるパイプライン処理に概念

Reiner Gamma と呼ばれる月磁気異常をモデルケ上は類似しているが、リストベクトルのような強ースとして採用し、プラズマ粒子シミュレーショ力な機構を持たないことや、ベクトルレジスタのンによる解析を実施した。ダイポール中心から磁数が少ないなどの理由により、複雑なデータアク気圧と太陽風動圧が釣り合う点までの距離を磁気セスパターンを有するアプリケーションへの適用

異常の代表長 L とすると、磁気異常ではその代表性は必ずしも高くない。また SIMD 命令発行の可長 L が太陽風電子のジャイロ半径よりも十分大き否は最終的にはコンパイラの判断に委ねるしかな

く、イオンのジャイロ半径より小さい。このよういため、カーネル計算部のプログラム構造やデー3九州大学先端的計算科学研究プロジェクト 平成 27 年度成果報告書 2016 年 4 月

タ構造を変更しながら、SIMD 命令が発行可能なE[][][]/B[][][]実装を見つけ出す、といった試行錯誤を繰り返す

for(p in a block)

こととなる。

4848

Update vp
Update xp
Scatter J((xp-vp))
Scatter J((xp))

電流計算では、各粒子がつくる電流要素を、それらの位置座標に隣接する格子点上の電流格子点

Update vp
Update xp
Scatter J((xp-vp))
Scatter J((xp))

48

配列に配分する。この粒子位置は実行時まで未知

48

J[][][]数であるため、格子点配列へのアクセスはしばし

図 2: SIMD 最適化前の PIC カーネルループ。

ばランダムになる。
Xeon Phi は 512 bit 幅の SIMD ベクトルを持つfor(each cell)

ため、原理上は 8 つの粒子を同時に処理すること

E[][][]/B[][][]

が当面の目標であるが、従来のプログラムにおい48

Scalarized e/bては、ループ内のプログラム構造が複雑すぎるとfor(part.in cell)

して、Intel コンパイルによる SIMD 命令発行がUpdate vp

Update vp

for(part.in cell)
Scatter J((xp))
Update xp

なされていなかった。この主な要因は、電磁場や電流配列に対する間接的なアクセスにある。粒子

Scatter J((xp))
Update xp
24

Scalarized j

ループの各繰り返しにおいて、速度更新の際には粒子に隣接する 8 格子点分の電磁場ベクトル（6J[][][]

Sort particles成分）、すなわち 48 成分の E((xp))や B((xp))を参for(each cell)for(part.in cell)

照する。ここで xp は浮動小数点であらわされた粒

Scatter J((xp))

子位置、関数はそれを周辺の 8 格子点に整数化す

Scatter J((xp))
24

Scalarized jる処理を表す。同様に電流計算においては、粒子J[][][]

図 3: SIMD 最適化後の PIC カーネルループ。

が存在するセルの辺上に定義される電流ベクトル
12 成分に参照書き込みが行われるが、それを粒子で共通化する．

移動前後それぞれについて行う必要があるため、

B) A がコンパイラにとって明らかになるよう、

計 24 成分の J((xpvpt))および J((xp))にアクセ電磁場や電流の参照書き込みにスカラー変数スする必要がある。ここで vp とt はそれぞれ粒子を用いる．

速度と時間ステップ幅を表し、(xpvpt)により粒A の対策は、粒子の所属セルをキーとしたソーテ子移動前の位置を示している。ここでコンパイラ

ィングで実現される。ここで、粒子ループ中ではは関数によって計算されるインデックスに規則粒子位置更新が行われるため、その際に所属セル性を見つけることができないため、8 つの粒子そが変化する可能性があり、その場合は A が満たされぞれに対する E((xp))、B((xp))、J((xpvpt))、れなくなる点に注意する必要がある。そこで従来そして J((xp))から構成される SIMD ベクトルの長大な粒子ループを、速度更新、位置更新前の48+24=72 本をレジスタに登録することで SIMD電流計算、位置更新後の電流計算という 3 つのル命令発行を試みる。しかし、現実には Xeon Phi でープに分割することで、各ループの中で常に A を登録可能な最大ベクトル数は 32 であるため最終満たすことが保証される。

的には SIMD 命令発行が不可であると判断されてB については、これまで 1 重の粒子ループであしまうのである。

ったものを、セルに関する外側ループとそのセル

この問題を解決するためには、次の 2 つの対策に所属する粒子の内側ループから成る、2 重ルーをとればよいと考えられる[6]。

プ構造に変更する。これによりセルループの内側A) 参照書き込みを行う格子点を 8 つの粒子の間では、(xp)は共通であるから、E((xp))、B((xp))
4九州大学先端的計算科学研究プロジェクト 平成 27 年度成果報告書 2016 年 4 月

を e や b、J((xp))を j といったスカラー変数にコタアクセスパターンを考慮することにより、カーピーすることが可能である。内側での粒子ループ

ネル計算部分の SIMD ベクトル化を促進できるこ

では、元の配列ではなく、e、b、そして j を対象にとが明らかとなった。今後は、これまでに蓄積し計算を実施する。これらを踏まえ、変更前後のルたノード間並列、ノード内並列、SIMD ベクトループ構造を表したのが、図 2 および 3 である。こ化の各技術要素を有機的に統合することで、エクのような変更を施すことにより、PIC 計算のカーサスケール時代のメニーコア型スーパーコンピュ

ネル部分である粒子速度位置更新、電流計算につータに対応可能なプラズマ粒子シミュレーション

いて、コンパイルによる SIMD 命令発行が可能と手法を確立していきたい。

なることを確認した。
現在、上述の最適化法の鍵となる、厳密な粒子ソ

6. 文献ーティング部分について、実装を進めている。こ(1) 参考文献れが完了し、正しい動作が確認され次第、性能測[1] Miyake and Usui, Phys. Plasmas, 16, 062904,

定によって SIMD 最適化の効果を評価する予定で2009.

ある。実際には、ソーティングはその原理上、完[2] Nakashima et al., Proc. 23rd ICS, 90-99, 2009.

全な SIMD 最適化実装が極めて困難な処理であり、 [3] Miyake and Nakashima, Proc. 11th IEEE
性能のボトルネックになり得る点を十分に認識しISPA, 2013.ている。このため、ソーティングを SIMD 処理可[6]

能な部分（粒子へのラベリング）と不可能な部分

Engineering, 46, 81-94, 2015.

Nakashima,

Computers&Electrical（メモリ上でのデータの並べ替え）を仕分けし、SIMD 実行可能部分の割合を最大化するような実(2) 当該プロジェクトによる研究成果装を検討中である[7]。

[4] Miyake, Y. and M. N. Nishino, Electrostatic
environment near lunar vertical hole: 3D
plasma particle simulations, Icarus, Vol. 260, pp.

5. まとめと今後の展望

301-307, 2015.

先端的計算科学研究プロジェクトによる九州大学 HA8000-tc/HT210 計算機資源を活用して、月[5] Usui, H., T. Matsubara, Y. Miyake, M. N.

周辺のプラズマ環境、より具体的には①月縦孔周Nishino, and J. Wang, Electron dynamics in the

辺の特徴的な電気環境、②縦孔電気環境下での帯

boundary current of a mini-magnetosphere

電ダストの挙動、③月磁気異常と太陽風プラズマ

above aの相互作用によって形成される小型磁気圏境界にsubmitted.

おける電子ダイナミクスの役割、に関して、重要[7] Nakashima, H., K. Kikura, Y. Summura, andな知見を得られた。今後、月上空 1 km までの物

Y. Miyake, Large Scale Manycore-Aware PIC

理描像である①および②の知見と、上空数 10～

Simulation with Efficient Particle Binning,

100 km で生起する③の物理描像の間の関連性をsubmitted.

解明することにより、月プラズマ科学の包括的な理解が得られることが期待される。
また計算科学的観点からは、OhHelp 技法により既に分散メモリ並列に対応した EMSES シミュレーターについて、さらなる計算コードの高度化を目指して SIMD 最適化の検討を行った。プラズ
マ粒子シミュレーション特有の演算フローとデー
5

lunar magnetic anomaly, to be