学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

jh180018-NAJ

超並列宇宙プラズマ粒子シミュレーションの研究
三宅洋平（神戸大学・課題代表者）、中島浩（京都大学）、加藤雄人（東北大学）
小路真史（名古屋大学）、臼井英之（神戸大学）、大村善治（京都大学）、岩下武史（北海道大学）
深沢圭一郎（京都大学）、佐伯拓哉（神戸大学）
概要 本研究の目的は、Particle-in-Cell 方式に基づく宇宙プラズマ粒子シミュレーションを最新の大規模並列計算機アーキテクチャへ効率的に実装し、宇宙プラズマ理工学の重要問題に適用することである。今年度は、メニーコア向け高効率プラズマ粒子計算に関して、集団通信アルゴリズムの再検討を通して、プロセス間通信の最適化を実施する。
また計算解法自体も陰解法の導入により、より大規模な計算対象に適用可能なモデルとする。これに加え、月や小惑星などの周辺で予想されるプラズマ電磁じょう乱現象、電子ハイブリッドシミュレーションによる、地球放射線帯の相対論的高エネルギー電子加速現象、の物理機構を大規模粒子シミュレーション解析により解明する。

1. 共同研究に関する情報

負荷不均衡について、PIC シミュレーショ

(1) 共同研究を実施した拠点名

ンは、膨大な数のプラズマ代表粒子と空間格

本課題は、神戸大学・京都大学・東北大学・子点上の電磁場の相互作用計算がシミュレー

名古屋大学・北海道大学の複数の研究者が参ションの大部分を占める。このとき、分割さ加する共同研究であり、京都大学の大規模計れた小領域間で粒子の粗密が生じると、負荷算資源を利用する。

不均衡が発生し、並列効率を著しく損なう。

(2) 共同研究分野：超大規模数値計算系応用分野

分散メモリ並列レベルでは、我々が独自に開(3) 参加研究者の役割分担発した動的負荷分散技法 OhHelp により、こ参加研究者の役割分担は過年度の体制に、の問題に対処している。

大学院生を加えた構成である。詳細は H26 年

不規則アクセスに関しては、格子点配列へ度最終報告書の表 1 を参照願いたい。宇宙プのアクセスパターンが粒子位置に依存して決ラズマ理工学分野の中でも粒子モデルシミュ

まる、いわゆる間接アクセスを有したままだ

レーションの専門家と流体計算の専門家、そと、メニーコアプロセッサが有する SIMD 演して計算機科学分野研究者・技術者が有機的算機構を有効活用することができない。このに協力して行う共同研究体制となっている。

問題に関して、我々は所属セルに応じた厳密な粒子ビニングを行うことで、複数粒子間で2. 研究の目的と意義間接アクセスを排することを提案した。

本研究の目的は、Particle-in-Cell (PIC) 方

上記 2 点の個別要素技術に加え、それらを式に基づく宇宙プラズマシミュレーションを組み合わせてハイブリッド並列実装を実現す

対象に、近年主流となりつつあるメニーコアる上で問題となる、粒子ビンのオーバーフロ型並列計算システムへの高効率実装を実現すーなどの例外的なイベントへの対策も検討しることにある。この目標に向け、不均一な粒た。これらの検討を通して開発したプロトタ

子分布に起因する負荷不均衡、粒子が所属すイプ版プラズマ粒子シミュレータによるるセルに依存して参照先が決定されることにXeon Phi KNL 上での性能評価では、粒子均よる、格子点配列への不規則アクセス、の問一なケースにおいて良好なスケーラビリティ

題を解決する必要がある。と、Xeon Phi と通常の Xeon プロセッサの理1

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

論的な性能差に見合った性能向上が得られた。

度な並列処理関連技術を開発するに至った。

一方で、粒子分布が極端な不均一なケースで実際に OhHelp では、木構造を用いた負荷均

は、プロセス間の通信性能が原因とみられる衡状況の再帰的なチェックなど、非数値的な性能劣化が顕著となった。ポスト「京」など発想を随所に活用しており、計算機科学から次世代のスーパーコンピュータシステムでは計算科学への貢献の代表例となっている。

通信コストが現在以上に切実なオーバーヘッ

4. 前年度までに得られた研究成果の概要

ド要因になると考えられ、上述の問題に正面から取り組む必要がある。

領域分割方式の超並列計算が困難とされて実問題への応用に関しては、OhHelp 適用

きたプラズマ粒子 PIC シミュレーションに動済の全粒子シミュレーションコード EMSES

的負荷分散技法 OhHelp を適用した結果、マにより、月惑星プラズマ／ダスト環境の大規ルチコア CPU 上では粒子局在化による性能模 3 次元解析を推進している。数 1000 コア劣化を大幅に低減することに成功した。またオーダーの大規模演算により、プラズマ荷電所属セルに応じた厳密な粒子ビニング法を提粒子による天体表面の帯電と、静電気力によ案し、その結果として従来の PIC コードで問って天体上空を浮遊するダスト微粒子のダイ題となっていたコンパイラによる SIMD 命令

ナミクスを解き明かす。

発行の阻害要因を排することに成功した。こ上述の課題を遂行することは OhHelp 法のれらの要素技術を実装して開発されたプロト有効性の実証という数値アルゴリズム的観点、タイプシミュレータは、近年普及が進むメニおよび超大規模モデルによる宇宙プラズマ粒ーコア型のクラスタシステムにおいて、粒子子シミュレーションの実現と最先端研究成果分布が均一の場合に良好なスケーラビリティの創出という宇宙科学的観点の双方においてを有することが確認された。一方で、粒子が大きな意義がある。

著しく局在している状況においては、事前の課題申請時に記載した本年度の技術開発計予想を上回る性能劣化が確認された。この性画を以下に再掲する。

能劣化は、動的負荷分散に伴って必要となる①現状のコード開発において主な性能劣化要プラズマ電流データの All-reduce 処理に起因因として特定されている、プラズマ電流縮約することが判明した。現在適切な All-reduce

通信の実装最適化を行う。

通信アルゴリズムの選択により、この性能劣②計算時間削減のための、計算アルゴリズム化を低減する方法を模索している。

上の改良を行う。具体的には電磁界ソルバー

応用分野における成果としては、本共同研に陰解法を導入する。

究プロジェクトで開発を推進してきた高性能

これらの点についての達成状況を 6 節に報告

計算技法により、時間変動する宇宙プラズマする。

環境下での衛星帯電現象の数値モデリングを開始した。この取組みにより、より広範な宇

3. 当拠点公募型共同研究として実施した意義

宙天気現象に対する衛星電位応答の理解が可

本研究は、アルゴリズム開発、コード開発、能になると強く期待される。

大規模シミュレーション実行の全ての段階で、理学的応用としては、地球放射線帯での相複数の大学の計算科学者と同センターの計算対論的電子生成機構において重要な役割を果機科学者との緊密な連携によって実施されてたすと考えられているプラズマ波動・コーラいる。このような研究体制により、動的負荷ス放射の発生過程を再現する計算機実験を実分散技法 OhHelp や粒子ビニング法などの高施し、その波動特性が背景磁場強度の空間構
2

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図1

粒子シミュレーションによって得られた月縦孔・地下空洞の表面電位分布。

造の変化に応じて大きく変化することを明らであり、月面電位の平衡状態は、荷電粒子種かにした。

間の電流平衡が満たされるように決定される。
図 1 は、粒子シミュレーションで求めた、5. 今年度の研究成果の詳細

月昼側に位置する縦孔・地下空洞周辺の表面

近年、月に複数発見されている縦孔地形に電位分布を、カラーマップで表示したもので近年注目が集まっている。2017 年には、月探ある。太陽天頂角は 30である。まず縦孔外の査機「かぐや」の観測により、それらの一つ月面電位は+5 V と光電子温度 2.2 eV の数倍が溶岩チューブにつながっていることが裏付の正電位となっており、一般に知られているけられた。縦孔や地下空洞は、電磁気環境、月昼側電位に一致している。縦孔・地下空洞

およびその影響を受けるダスト環境において内部では、部分的に太陽光が遮断され、日陰特異な様相を示す。本研究では、プラズマ物の領域が生じる。日向と日陰でそれぞれ正と理分野や人工衛星帯電解析に用いられてきた負の帯電が生じており、日向ではおよそ+20粒子シミュレーション技術を、月縦孔・地下V、日陰では−40 V 程度の電位が形成される。

空洞周辺の電磁気・ダスト環境の解析に応用

特に日向の縦孔底面からは、縦孔外部と同様し、縦孔周辺の帯電微粒子環境が縦孔探査にの条件で光電子が放出されるが、最終的な電及ぼす影響を評価する。

位は縦孔内部の方が大きい結果となっている。

本計算では、Marius の丘に発見されている

上述のようなプロセスにより月表面のレゴ縦孔のデータを参考に、そのサイズを直径 50リス層が帯電し、その一部はレゴリス層直上m、深さ 45 m とする。また縦孔の下部が、天の特異な電場環境のもとで浮遊し、特徴的な井までの高さが 20 m の地下空洞につながっダイナミクスを示す。これらの浮遊ダストがている状況を想定し解析を実施した。月面で月面探査機に付着することにより、機体表面の荷電粒子の吸収・放出過程は、月レゴリスの摩耗や異常帯電、熱特性の予期しない変化層に蓄積される電荷の符号と大きさを規定し、が引き起こされる恐れがある。そこでダスト表面電位を決定づける。一般に固体表面におの挙動を明らかにすることを念頭に、テストけるプラズマ電流は、物体の表面電位の関数粒子解析を実施した。今回は、静電気力が月として表される。これは月面においても同様面重力に打ち勝ち、浮遊する条件を満たすダ3

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図2

月縦孔・地下空洞周辺の代表的なダストの軌跡。それぞれ初期に (a) 縦孔外、(b) 縦孔内日向領域、(c) 縦孔内日陰領域、に配置されたダストの運動を示す。

ストとして、数十 nm の微小ダストを計算のな運動の非対称性は日向・日陰空間中の荷電対象とする。

粒子密度の違いに起因しているが判明した。

図 2 に縦孔周辺でのダスト挙動パターンを

上記の結果に基づくと、数十 nm スケール数例示す。まず図 2(a)に示す縦孔外の正帯電の微小ダストは，縦孔や地下空洞の日陰領域

ダストは月面上を浮遊し、最終的に縦孔日陰に引き寄せられやすいことが，強く示唆され領域に到達していることが図から確認できる。

る。今後，ダストサイズに対する依存性を調

ダストの持つ電荷は日向では正であるが、同査し，日向‐日陰を横切るダストフラックス

時に日向の月面からは上向きの電場が発生しの見積もりを実施する予定である。ているため、ダストが浮かび上がる向きに静電気力が作用する。図 1 で見たように縦孔の6. 今年度の進捗状況と今後の展望

日陰領域は負に帯電しているため、浮遊した(1) 上半期の活動：プラズマ電流密度縮約処理にまま縦孔直上に到達したダストは、当該領域かかるプロセス間通信方式の検討に引き込まれる結果となる。また図 2(b)およまず 2 節の最終段落で言及した①の課題にび(c)に示すあらかじめ縦孔内部の日向と日陰対する進捗を報告する。これまでに開発したに配置されたダストは、両者とも初め日向－要素技術を組み込んだプロトタイプシミュレ

日陰境界を交差するような運動を示す。このータによる、京都大学 Cray XC40 の Intel

時点では、日向・日陰内ダストの挙動は対称

Xeon Phi KNL 7250 上での性能評価では、粒的であるが、その後、経過は異なる。日陰か子均一なケースで良好なスケーラビリティがら日向に移動したダスト（図 2(b)）は数 10 秒得られたものの、粒子分布が不均一なケース～数分の時間後に再び、日陰領域に帰っていでは性能劣化が顕著であった（前年度最終報

く。しかし当初、日向から日陰に移動したダ告書を参照）。この性能劣化については、MPIストについては、数 10 分の数値計算では、こ通信、特にプラズマ電流の All-reduce 処理にうした帰還運動は見られなかった。このよう伴うオーバーヘッドが主要因であることが既4

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図3

リングシフト通信のマルチスレッド化の効果．に判明している。

で、導入が可能であり、またプロセス内の演

上記の課題を踏まえ、前年度下半期より、算のみが対象のため、マルチスレッド化の効

当該 MPI 通信部分の最適化の取り組みを開用は比較的予測しやすい。一方で、前者につ始した。OhHelp 適用シミュレータで出現す

いては、先述した特殊な All-reduce 処理に適る、データ量：大、かつプロセス数：小の多数用する前段階として、比較的簡単なデータ通同時並行 All-Reduce 処理という特殊な通信信モデルにマルチスレッド実行を適用するこパターンにおいては、Ring 方式による All-とで、その効果を事前評価することが必要でreduce 通信アルゴリズムが比較的有効であある。
今回は京都大学 XC40 を 16 ノード使用し、ることを見出した。
また上述の性能分析においては、KNL プロ

様々なノード内プロセス数 P とスレッド数 T

セッサ内の全コアをスレッド並列に用いるこ（PT=64 コア）の組み合わせに対して、リとが有効でないことも判明している。これはングシフトのデータ通信を行う実験を行った。

現状の実装では、各プロセス内で単一のスレこれは All-reduce の通信アルゴリズムとしてッドのみが MPI 通信を担当していることとRing 方式が有効であるという、先述の検討結関係している。すなわち KNL プロセッサ内果に基づく実験内容である。通信するデータで 1 プロセスのみ起動し、前コアをスレッドサイズはスレッドあたり 80 KB とする。この並列に用いた場合、ノード間通信に関与する

ときノード当たりではおよそ 5 MB 程度と計コアはノード内で 1 つと極端に少ないため、算され、実際のプラズマ粒子シミュレーショノードが本来持つ通信バンド幅を充足していンで用いるノード当たりの電流データ量に対ないと考えられる。そこで今年度からは電流応したサイズとなる。複数スレッドによる通All-reduce 処理のマルチスレッド化も検討す

信実装には MPI_Thread_Multiple を用いた。
図 3 にノード当たりのプロセス数（プロセる。
実際のマルチスレッド化の対象としては、ス当たりのスレッド数）をパラメータとした、All-reduce に必要な通信自体と、通信された通信性能の測定結果を示す。図中の青線は従データの縮約計算の部分に分けられる。後者

来どおり単一スレッドのみが通信を実施する

は、All-reduce 通信を自作の実装とすること

場合、緑および赤線は複数スレッドで通信を5

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

実施した場合である（緑と赤の違いについて縮することが可能である。

は後述）。図からマルチスレッド化は通信性能

本課題では、Hoshino [1986]により提案さの改善においては逆効果であることが判明しれ、Matsumoto et al. [2012]等に採用されてた。特に、ノード内のプロセス数が少ない（＝いる定式化を適用することとした。詳細は省1 プロセス内のスレッド数が多い）ケースに略するが本定式化を用いると各格子点上の磁おいて、この傾向が顕著である。

界 B の時間発展を表現する連立一次方程式が得られる。Hoshino [1986]らの手法では、こ複数スレッドによる通信が逆効果になった
原因としては、以下が考えらえる。

れを共役勾配法に基づく反復法を用いて解い・単一スレッド通信の性能が最大化されるている。しかし今回我々は上式を離散フーリ

48 procs/node では、既にノード間のリンクエ空間に写像した上で、直説法により解く方性能がほぼ充当されているため、複数スレッ法（スペクトル法）を採用した。
構築した電磁界ソルバーを用いた 2 次元電ド化による恩恵が少ない。
・1 procs/node（64 threads/proc）では、複数磁波伝搬のテスト計算結果を図 4 に示す。図スレッド間の排他制御によるオーバーヘッド

中の(a)～(d)は 4 種類の異なる時間刻み幅をが大きくなり、性能劣化につながっている。

用いた計算結果のスナップショットを示して

また XC40 では、複数スレッド間の排他制おり、各ケースの時間刻み幅をクーラン数御を最小化するオプションが提供されており、

𝑐∆𝑡/∆𝑟（𝑐、∆𝑡、∆𝑟はそれぞれ光速、および時

図 3 の赤線はそれを有効にした場合の結果で間と空間の刻み幅）によって表示している。

あるが、その効果は限定的であった。

陽解法を用いた従来コードでは、クーラン数が1/√𝑁d （𝑁d は空間次元数．テスト計算で(2) 下半期の活動：陰解法電磁界ソルバーの導入

は𝑁d = 2）を超えると計算が不安定化し、数による時間更新刻み幅制約条件の緩和値発散が生じる。陰解法の適用により、クー

次に 2 節最終段落の課題②の進捗を報告すラン数に対する制約が取り払われ、全てのケる。大規模 PIC シミュレーションの高速化をースにおいて安定的に電磁波伝搬が再現され実現するためには、計算機アーキテクチャにていることがわかる。この場合、解くべき振合わせた最適化に加え、計算アルゴリズムそ動数（周波数）を十分に解像できる時間刻みのものの高度化も必要である。特に改善が見幅で計算を進めることができるため、従来に込まれるのは、電磁界ソルバーの部分である。

比べ各段な計算時間短縮が可能となった。

これまで用いていた陽解法では、時間刻み幅

現在、開発された陰解法電磁界ソルバーをは光速に対するクーラン条件によって拘束さ利用したプラズマ粒子シミュレーションの動れる。一方でプラズマ荷電粒子の興味あるダ作テストを実施しており、2019 年度中には実イナミクスは、光が 1 セルを伝搬する時間よ用化できる見込みである。なお上述のスペクりはるかにゆっくりとしているため、上述のトル法を用いた陰解法には、3 次元のフーリ

時間刻み幅は本来の目的に照らして過剰に小エ変換が必要である。フーリエ空間の並列化

さく、計算時間が増大する結果となっている。

においては、All-to-all 処理によるデータ転置

この問題を解決するため、今年度下半期でに基づく標準的な実装を用いている。このは電磁界ソルバーに陰解法を導入した。これ

All-to-all 通信のため、電磁界ソルバー単体でにより、時間刻み幅は上述のクーラン条件では、大規模並列計算時のスケーラビリティがはなく、興味ある現象を十分に解像できる値劣化するものの、プラズマ粒子計算部分に比に設定することができ、大幅に計算時間を短べて演算時間が小さいため、PIC シミュレー
6

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図4

陰解法適用後の電波伝搬テスト結果．(b)～(d)は陽解法における不安定領域に対応する．

ション全体のパフォーマンスに与える影響はべて大きな時間刻み幅を使用することが可能軽微であると予想される。

である。本年度中に電磁界ソルバー単体の実装作業が完了し、今後 PIC 計算に組み込んで(3) 課題遂行の総括と今後の展望テスト計算と性能評価を実施する予定である。

本年度は、大規模並列計算システムに適合課題申請時の研究計画に準拠した工程を推した PIC シミュレーションのプロセス間通信進し、一部期待した成果が得られなかった部

高効率化の検討と、時間刻み幅に対する拘束分はあるものの、着実に研究は進展している。

条件の緩和を目指した改良計算アルゴリズム

なお、上述したように、特定の通信処理のコの実装に取り組んだ。

ストが全体の計算時間の中で大きな割合を占開発中のメニーコア向け PIC コードにおけめるようになってきているが、これは数年にる性能劣化要因は既に特定されており、該当わたるコード最適化検討で他の部分の計算コするプラズマ電流縮約通信について検討を実ストを大幅に削減してきた結果でもある。こ施した。結果として、当該 MPI 通信そのもののような状況でさらなる高効率化を進めるにのマルチスレッド化は通信性能の向上には有は、アルゴリズムそのものの改良など、多方効ではなく、All-reduce 処理全体の性能を向面からの検討が必要になると考えられる。

上する上では、縮約の演算に限ったマルチスレッド化が有効であることが強く示唆された。 7.
PIC シミュレーションに含まれる電磁界ソ研究成果リスト(1) 学術論文
1.

ルバーに陰解法を導入することで、従来に比7Paulsson, J.J.P., Y. Miyake, W. J. Miloch,

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月and H. Usui, Effects of Booms of Sounding

2.

8.

Rockets in Flowing Plasmas, Phys. Plasmas,

the Nonlinear Processes of Whistler-Mode26, 32902, doi:10.1063/1.5051414, 2019.

Chorus Generation in the Earth's Inner

Ahmad, N., D. Herdiwijaya, T. Djamaluddin,

Magnetosphere, 13th ISSS, Los Angeles,

H. Usui, and Y. Miyake, Diagnosing Low

2018.

Earth Orbit Satellite Anomalies using

3.

9.

Usui, Numerical Simulation of Spacecraft

Geomagnetic

Charging

Processes

Planets, and Space, 70, doi:10.1186/s40623-

PlasmaEnvironment,018-0852-2, 2018.

Charging Technology Conference, Kobe,

Katoh, Y., Y. Omura, Y. Miyake, H. Usui,

2018.

H.

Perturbations",

Nakashima,

GenerationofEarth,

Dependence

WhistlerMode

15thSpacecraft

ChorusUsui, Numerical Simulation of Spacecraft

Emissions on the Temperature Anisotropy

Charging

and Density of Energetic Electrons in thePlasma Environment, 2nd URSI AT-RASC,

Earth’s Inner Magnetosphere. J. Geophys.

Gran Canaria, 2018.

Processes

inTime-varying(4) 国内会議発表
11. 三宅洋平, 西野真木, 月面近傍プラズマ・ダス

Miyake, Y.,, Y. Funaki, M.N. Nishino, and H.

ト環境の粒子モデルシミュレーション, 第Usui, Particle Simulations of Electric and

144 回 SGEPSS 総会および講演会, 名古屋,Dust Environment near the Lunar Vertical2018.

Hole, Proc. 8th International Conference on

12. 三宅洋平, 粒子モデル プラズマ運動論シミ

the Physics of Dusty Plasmas, 1925, 020001,ュレーションの高効率計算技法, 第二回ポス2018.

ト「京」時代の天体形成シミュレーション研究会, 神戸, 2018.

(3) 国際会議発表

13. 三宅洋平, 寸村良樹, 佐伯拓哉, 中島浩, メニ

Katoh, Y., Rapid Acceleration of Relativistic
Wave-Particleーコア型クラスタ向け高効率プラズマ粒子シInteractions in Planetary Magnetospheres,ミュレータの開発, 日本地球惑星科学連合大Electronsthrough会, 幕張, 2018.

International Symposium on "Physics at the

14. 三宅洋平, 月縦孔および地下空洞におけるプ

Cosmic Frontier", Toyama, 2019.

ラズマ・ダスト輸送の数値シミュレーション,Nakashima, H., An Implementation of PIC

日本地球惑星科学連合大会, 幕張, 2018.

Method with OhHelp Load Balancing

15. 三宅洋平, 粒子モデルプラズマ運動論シミュ

Library for Manycore Processors, 24th

7.

Time-varying

10. Miyake, Y., T. Kiriyama, Y. Katoh, and H.

(2) 国際会議プロシーディングス

6.

inofRes., 123, 1165–1177, 2018.

5.

Miyake, Y., T. Kiriyama, Y. Katoh, and H.

NOAA-15 Electron Data Associated withand4.

Katoh, Y., Electron Hybrid Simulations of

Numerical EXperiment of Tokamak (NEXT)

レーションの高効率計算技法, 第二回ポスト

Workshop, Rokkasho, 2018.

「京」時代の天体形成シミュレーション研究会, 神戸, 2018.

Miyake, Y., Numerical Simulations of Dust
Plasma Environment near Lunar Surface,(5) その他（特許，プレス発表，著書等）

13th International Symposium of Space

該当なし

Simulations (ISSS), Los Angeles, 2018.
8