三宅洋平
神戸大学計算科学教育センター

y-miyake@eagle.kobe-u.ac.jp
1

概要

宇宙空間中の小型天体・人工衛星は電離気体であるプラズマとの接触により帯電する．またそれらの固体表面が宇宙プラズマの境界となり，その周辺ではプラズマ環境そのものが大きく乱される．
将来衛星ミッションに向け、このような衛星プラズマ相互作用を事前に定量評価する必要がある．本講演では，スーパーコンピュータ向けに最適化された電磁粒子プラズマ粒子シミュレーションを用いた小型天体・プラズマ相互作用解析の最新の取り組みを紹介する．
2

目次

1. 宇宙プラズマ環境と小型天体・人工衛星の相互作用
2. 宇宙プラズマ現象の数値的アプローチ
3. 宇宙プラズマ粒子計算の並列化技法

4. 解析事例
5. 今後の展望とまとめ

3

対象とする宇宙空間
我々が対象としている「宇宙」－ 地球近傍の人類生存基盤領域－ 太陽から地球に至る領域－ 太陽系惑星領域

4太陽活動による宇宙電磁環境への影響

copyright@NASA5

太陽－地球環境：太陽風と固有の地球磁気圏

太陽

地球太陽風：
磁場を伴った高速プラズマ流

地球磁気圏• 領域毎に全く異なるプラズマ環境• フレアなどの太陽活動に伴うダイナミックな変動→宇宙天気

6

宇宙プラズマ• 宇宙の99.9%以上はプラズマ状態• プラズマ： 物質の第4の状態。荷電粒子（電子、イオン）の集合体• 宇宙環境利用には宇宙プラズマの理解が必要

ローレンツ力（電場の力、磁場の力）

7

宇宙プラズマ中での人工衛星帯電
e +++++ e +eee
e + ee+ +e ++ee
e+e+ e+ e
e e ++ e +e+ e
e +ee
e+++e+e ++ e +e
e +e + ee e

8帯電による人工衛星障害
2003年10月
みどり2号
全損事故

MLI
（多層断熱材）

ハーネス
オーロラ帯でのMLI内面の帯電

故障推定箇所（大電力ハーネス部）

トリガ
放電

MLI内面とハーネス間の放電

ハーネス間の持続放電

持続放電の波及

104本全ての焼損
9衛星帯電・放電による衛星故障

1973年～1997年の衛星障害326件の統計
米国Aero-space社と米空軍、アメリカ海洋大気局による調査結果

半分以上が衛星帯電および放電による障害
10

数値シミュレーションの必要性

宇宙プラズマ環境をどうやって評価するか個々のパーツを地上プラズマチャンバー内でテストする。
•プラズマ環境の実現が容易ではない•チャンバー壁の影響がある•プラズマ現象が複雑（異方性、分散性、非線形性）•大型システムの評価が困難

プラズマシミュレーションを数値チェンバーとして利用

11

宇宙プラズマ現象の数値的アプローチ12

プラズマダイナミクスの記述法②オイラー(Euler)的記述
時刻tにおけるプラズマの状態は？

位相空間 (x,y,z,vx,vy,vz)
における存在確率（密度）v
体積要素①ラグランジュ(Lagrange)的記述

粒子位置

f (x,v,t)dv

位相空間密度

個々の粒子の位置、速度を追う

X2(t), V2(t)…

xdx

X1(t) , V1(t)位相空間密度の記述
運動の記述

個々の粒子の運動方程式

位相空間密度fの輸送方程式（ブラゾフ方程式）

粒子（Particle-in-cell）モデル

流体モデル

13何を解き進めるのか。
宇宙プラズマの動き（Dynamics）とそれによって生じる電磁場変動（いわゆる電波）。互いに影響しあう。
電波（電界、磁界からなる）を解く方程式（Maxwell方程式）

B
E  t1 E
  B  oJ  2c t


E 


o

B  0プラズマダイナミクスを解く方程式（運動方程式）

dvi qs (E  vi  B)
dt msdxi
 vidt
2つの方程式を組み合わせて、電磁場（電波）とプラズマの動きを解き進める。
14

プラズマ粒子シミュレーション
超粒子

空間格子
（電磁界、電流）

粒子の速度値を近傍格子点に配分し、格子点で電流計算
各格子点において電流から電磁界更新(Maxwell式)

各粒子の速度、位置更新（運動方程式）15

プラズマ粒子計算モデルの妥当性評価
1. 基本物理量の確認エネルギー保存、プラズマ波動分散関係…
2. 基礎現象の再現（理論との比較）
ラングミュアプローブ、OML理論、シース理論…
3. 科学衛星観測との比較
例）プラズマ航跡電場計測との比較
衛星観測
magnetotail lobes

just outside magnetopause

CLUSTER衛星（Credit: ESA）

シミュレーションfSC = + 40 V

fSC = + 20 V
16

計算例：衛星搭載イオンエンジンプラズマ挙動

ムービー再生17

スーパーコンピュータの利用

http://jp.fujitsu.com/about/tech/k/

 3次元かつ数値ノイズの少ないシミュレーション解析‒ 膨大な数の粒子（～1010個＝必要主記憶～TB）を解き進める必要‒ スーパーコンピュータの利用を前提
 近年のスーパーコンピュータは超並列システム‒ シミュレータの階層的並列化は必須

18プラズマ粒子計算の階層的並列化「プロセス並列化×スレッド並列化×SIMDベクトル化」

Process2C) SIMDベクトル化（複数粒子同時処理）
Via
粒子処理ループ構造の変更
粒子位置に基づくソーティング

Process3

Thrd2 Thrd3
Thrd0 Thrd1A) 高プロセス並列化
via
動的負荷分散OhHelp

MPI通信B) 高スレッド並列化
via
マルチカラースケジューリング

小計算領域
19

Process0

Process1

プラズマ粒子計算の分散メモリ並列技法

20

粒子分割法（例えば、[Martino et al., Parallel Comput., 27(3), 2001]• 全プロセスが全空間の電磁場情報を保有／計算• 粒子数を均等分割• プロセス毎に求められた電流を総和（all-reduce）○ 実装とても簡単
× 領域大メモリ不足
× 領域大＆並列度大電流計算コスト大大規模 all-reduce• 1セルあたりの粒子数大（≒領域小）21• 低並列数（102並列） • 低空間次元数（1次元）…ならば、これでもオッケー。

粒子分割PICコードの並列性能
3次元PICコード

別の方法
（後述）

1次元PICコード
別の方法
（後述）

粒子分割法

粒子分割法

⇒数10プロセスで飽和⇒数1000プロセスで飽和
22

（静的）領域分割法• 空間を小領域に分割（1～3次元）• 小領域を各プロセスに割り当て• 割り当てられた小領域内で粒子－格子の相互作用計算○ 実装まずまず簡単
× 粒子集中メモリ不足
× 粒子数のばらつき負荷不均衡並列効率低下• 粒子の粗密が少ない
23（一様粒子分布、非一様性が1or2次元的）…ならば、これでもオッケー。

（静的）領域分割法のMPI通信

通信

小領域1 小領域2+足し込み

J

通信

+足し込み

J

通信

24

より技巧的な分散メモリ並列技法

25動的領域分割法（例えば、[Thackera et al., Proc. HPCSA, 2003]• 空間を小領域に分割• 小領域を各プロセスに割り当て• 小領域内粒子数が均衡するように境界を動的に変更△遠隔作用の近似計算向き△負荷均衡複数領域担当（完全均衡は困難）△境界処理やや面倒
×低密度領域が大×負荷変動への対応が面倒

Octree○粒子数完全均衡
×低密度領域が大×負荷変動への対応が面倒
×境界処理かなり面倒

Orthogonal Recursive
Bisection (ORB)

26

（いろいろ登場したので）手法間比較
N: ノード数, P: 全空間中の総粒子数, D: 全空間の格子点数
For efficiency, particle & field solving costs should scale as O (P/N) and O (D/N).

Approach

粒子負荷(Ideal : Worst)

場の負荷(Ideal : Worst)

1. 粒子分割

P/N : P/N

D:D

2. 静的領域分割

P/N : P

D/N : D/N

3. 動的領域分割

P/N : P/N(?)

D/N : ~D

4. 他に方法は？

？

？

27

領域分割方式PICにおける粒子負荷不均衡の問題⇒粒子・格子連成計算に共通に現れる課題
粒子数

粒子が
局在

03

13233302

12223201

11213100

102030

ブロック領域分割

各小領域に所属する粒子数33
計算ノード

負荷不均衡！
28動的負荷分散アルゴリズムの必要性

OhHelp: 負荷分散戦略（[Nakashima et al., Proc. ICS, 2009]など）
0313233302 03 12 21 22 23 3200 3132103301
11 20 21 01 31 02
1130
13220010203033

One-handed help
(OhHelp)

粒子数

粒子の「一部」と場の「コピー」を委譲

平均粒子数

ノード番号 33 00 32 01 30 10 13 03 23 20 31 02 11 21 12 22「粒子数」と「小領域サイズ」双方の観点で負荷均衡化

29

OhHelp: 手助けノード（２次）割当のアルゴリズム

低負荷ノードが平均粒子数を持つように高負荷ノードの粒子の一部を担当

高負荷ノード側は一時的に自分が平均以下になっても粒子を供出
後で誰かからもらえるav. #p

33 00 32 01 30 10 13 03 23 20 31 02 11 21 12 22 Node ID•各ノードは他のノードを一つだけ「助ける」．例：“33” helps only “22”.
30•各ノードは複数のノードに「助けられる」場合がある．例：“22” ishelped by “23”, “32”, and “33”.

（OhHelpも含めてあらためて）手法間比較
N: ノード数, P: 全空間中の総粒子数, D: 全空間の格子点数
For efficiency, particle & field solving costs should scale as O (P/N) and O (D/N).

Approach

粒子負荷(Ideal : Worst)

場の負荷(Ideal : Worst)

1. 粒子分割

P/N : P/N

D:D

2. 静的領域分割

P/N : P

D/N : D/N

3. 動的領域分割

P/N : P/N(?)

D/N : ~DP/N : (1+a) (P/N)

D/N : 2D/N

4. OhHelp31

OhHelp実装PICコードの並列性能
並列手法間の比較（3次元PIC）

粒子分布の影響（OhHelp法）

OhHelp
均一粒子分布

粒子分割法

局在粒子分布

プラズマ粒子計算の超並列化の最大の障壁である負荷不均衡を（許容可能な追加コストで）解決32

共有メモリ並列・SIMDベクトル化

33

シミュレーションモデル
共有メモリ（スレッド）並列• 粒子分割 or 領域分割？ ⇒ 領域分割
（プロセス毎の小領域を内部でさらに分割）• 色付けによるアクセス衝突回避アルゴリズム

jp
jp

SIMDベクトル化• 基本的にはコンパイラの最適化機能に依拠• PICカーネル計算（プラズマ電流）に含まれるJ( d(xp) ) += jp

jp

J：電流配列、d：粒子隣接格子点を求める関数、jp：プラズマ粒子が作る電流の計算が自動ベクトル化を阻害• 粒子所属セルに基づく粒子並び替えによる解決

34

Elapsed time[sec]共有メモリ並列・SIMD最適化の効果

good

x17.6 (6.6)

Xeon Phi
5120D

x3.1 (2.5)

Xeon E5-2600

x3.9 (2.4)

SPARC64 (FX10)

メニーコアプロセッサ（Intel Xeon Phi）など新世代のプロセッサでは、このレベルの最適化が（も）特に重要
35

計算科学・計算機科学の連携
JHPCN(学際大規模情報基盤共同利用・共同研究拠点) collaboration research「超並列宇宙プラズマ粒子シミュレーションの研究」

神大：三宅、臼井
東北大：加藤名大：小路etc.

Applications•Spacecraft-plasma interactions•Plasma wave generation & particle acceleration
processes in Earth’s magnetosphere

京大：中島、深沢
北大：岩下etc.

HPC•Parallel & distributed algorithms•High-performance computational& programming models

最新のコンピュータアーキテクチャに適合した実用コードをコデザイン
互いの分野の研究内容・文化をよく理解する必要
36

解析事例

37

太陽探査衛星と太陽コロナプラズマの相互作用

Solar Probe Plus (Credit: NASA & JHU/APL)•高温高密度の太陽コロナプラズマ•衛星から超高密度の光電子・二次電子放出 衛星帯電特性は地球周回衛星と同じか？
 科学衛星観測への影響は？

38

シミュレーションモデル
z

イオン密度：ni

太陽風

プラズマ流

yx
衛星進行方向Solar Probe Plus(SPP)
ウェイク

x太陽コロナプラズマ（電子、イオン）、光電子、二次電子を考慮
39

シミュレーションによる衛星帯電予測
Case

PIC simulation

Theory*太陽近傍(0.04 AU)
水星周辺

20.5 V

+2.9 V+9.9 V+7.4 V地球周辺+12.8 V+9.5 V*OML理論[Mott-Smith andLangmuir, 1926]による値

光電子＋

太陽最近傍では、理論と大きく異なり衛星が負に帯電（地球周辺では理論値と大体一致）

40

衛星周辺の空間電位分布
負の電位バリア(-35 V)e−(V)負電位領域（35 V）＝電子にとっては障壁
e−ee−

e−
e−

−e−

e−

e−e−

衛星(20 V)

e−e−

e−

e−

衛星表面（ 20V）

光電子電流を大きく制限(背景電子流入による負電荷蓄積)＞(光電子放出による正電荷蓄積)→衛星は負に帯電太陽近傍環境特有の現象

41

自然の小型天体（月、小惑星）への応用
太陽風プラズマ－月固体表面間の電磁気的相互作用（局所的には月の磁気異常と太陽風プラズマの相互作用）Sun
42[Credit: Halekas、Delory, Farrell and Stubbs]

月縦孔地形の電気環境のシミュレーション解析
Marius hills

かぐや衛星により、クレーターとは「異なる」縦孔地形の発見⇒月・火星の縦孔・地下空洞探査計画UZUME

マリウスの丘の縦孔

将来の直接探査・月面基地利用を見据えて…
月縦孔の電気環境および帯電ダスト環境の数値予測

43

計算モデル・条件
q

孔の深さ
45 ｍzyx

月面： 1. 絶縁性
2. 荷電粒子吸着
背景プラズマ：
太陽方向から到来
光電子放出：
太陽光到来角を考慮
格子点幅：50 cm
解析空間：
200  200  1000 m3
使用計算資源：2048コア、3.5 TB
44

計算結果：月縦孔表面の帯電分布
q = 30

日照面は正に帯電
電界：50 V/m

電界：20 V/m

縦孔底面がより高電位45

電位差（> 40 V）

計算結果：月縦孔周辺の帯電ダスト挙動
縦孔周辺の電気環境（月面電位）

帯電ダストの挙動予測“Particle-in-Cell計算”

出力データ利用“ダスト挙動(Test-particle)計算”

今後帯電ダストの月面着陸探査機に付着した場合の、影響を評価し、月縦孔探査ミッション設計に役立てる。

46

今後の展望：ポスト「京」に向けた宇宙環境シミュレータの高度化
ポスト「京」で重点的に取り組むべき社会的・科学的課題（萌芽的課題）┗太陽系外惑星の誕生と太陽系内惑星環境変動の解明┗太陽活動による地球環境変動の解明
グローバル磁気圏スケール
シミュレーション（流体モデル）…に参画人工衛星近傍スケールシミュレーション（本モデル）

47

各時空間スケールシミュレーションコードの次世代アーキテクチャ向け最適化、およびスケール間の連成研究に向けた技術的検討

まとめ• 宇宙空間を満たすプラズマと小型天体や人工衛星の相互作用は宇宙探査ミッションに多大な影響を与えうる。
• 人類未踏領域における人工衛星プラズマ相互作用や月や小惑星などの小型天体周辺のプラズマ電気環境に関して、研究の余地が十分に残されている。
• 仮想的に宇宙環境を再現可能なシミュレーションは、上記の問題を解決する上で非常に有用なツールである。
• 本質的に3次元、非均一な系を対象とする本問題ではスパコンの利用は必須である。計算機科学分野との緊密な連携により大きな成果を上げることができる。
48

ご清聴ありがとうございました

49