プラズマシミュレーション
における粒子法
梅田隆行
名古屋大学太陽地球環境研究所
Particle-In-Cell (PIC) 法
• 基本的な手法は1970年頃に完成
• (宇宙)プラズマ科学&電波科学からスタート
• 古典物理的な第一原理基礎方程式
荷
電
粒
子
の
運
動
Sir Isaac Newton Hendrik Antoon Lorentz
James Clerk Maxwell
電
磁
界
Newton’s law (1678?) Lorentz force (1892?)
Maxwell’s equation (1864?)
基礎方程式
 dx n
 vn

 dt
 dv
 n  qn E  v n  B 

mn
 dt
x  ( x, y , z )
v  (v x , v y , v z )
1 E

  B  0 J  2


c t

  E   B

t

FDTD法で数値的に満たす
B  0
(Yee格子と呼ばれる staggered格子)
K. S. Yee, IEEE Trans. Antenn. Propagat.,
AP-14, 302 (1966).

J 
0
t
電荷保存法


E 


0

J   qn v n

n

O. Buneman, J. Comput. Phys., 1, 517 (1967).
J. P. Boris, Proc. Fourth Conf. Num. Sim. Plasmas, ed. J. P. Boris and R. A. Shanny,
p.3 (Naval Research Laboratory, Washington D. C., 1970).
プラズマPICシミュレーション
• フル（全）PIC：第一原理、電子・イオン共に
すべての粒子の運動を解く
◎「ある制限内」では厳密。現在では「厳密性」は
薄れつつある…
– 無衝突性
– 格子内に十分な（～数十個以上の）粒子が存在
• （初期値が微妙に違う）同じ計算を何回もやって
平均を取る、ということはしない
• 応用分野では、「PICの計算は正しい」という
誤解が蔓延…
デバイ遮蔽
※イメージです(笑)
・電子
○イオン
• 無衝突な荷電粒子の集団中に存在する、ある点電荷が作る静電
ポテンシャルは、真空中の値のexp(-r/lD)になる。
• 点電荷からの距離がlDになると、点電荷の電場の影響はほぼ
無い。（exp(-1)=0.368）
応用分野では
⇒格子幅はlD若しくはそれ未満である必要がある。
たぶん無視
⇒幅lDの範囲に内に十分な粒子数がある。
されている…
プラズマPICコードのカーネル
• Leap-Frog (時間2次精度)
qn
 dv n
E  v n  B 


mn
 dt

 dx n  v
n

 dt
その他（境界処理、
ロードバランスetc）
t –Dt/2

J 
0
t
1 E

  B  0 J  2


c t

  E   B

t

t
t +Dt/2
t +Dt
t +3Dt/2
x v J E B
粒子形状
（shape/form factor/function）
• 零次（Nearest Grid Point）
Q(i)=Q(i)+1.0
ほぼ使われない
• 一次（Square-Shaped Cloud）
Q(i )=Q(i )+0.5+ax
Q(i+1)=Q(i+1)+0.5-ax
• 二次（Triangular-Shaped Could）
Q(i-1)=Q(i-1)+(0.5 -ax)**2*0.5
Q(i )=Q(i )+ 0.75-ax **2
Q(i+1)=Q(i+1)+(0.5 +ax)**2*0.5
１次精度
Q(i ,j )=Q(i ,j )+(0.5+ax)*(0.5+ay)
Q(i+1,j )=Q(i+1,j )+(0.5-ax)*(0.5+ay)
Q(i ,j+1)=Q(i ,j+1)+(0.5+ax)*(0.5-ay)
Q(i+1,j+1)=Q(i+1,j+1)+(0.5-ax)*(0.5-ay)
ロードが2個(Q(i,j),Q(i,j+1))、ストアが2個
加減算が8個、掛け算が4個
ロードが3個、ストアが3個
⇒要求B/F:16/12=1.33
加減算が18個、掛け算が33個
実際：加減算が4個、掛け算が4個
⇒要求B/F:24/51=0.47
⇒要求B/F:16/8=2
実際：加減算が6個、掛け算が19個
２次精度
⇒要求B/F:24/25=0.96
Q(i-1,j-1)=Q(i-1,j-1)+(0.5 -ax)**2*(0.5 -ay)**2*0.25
Q(i ,j-1)=Q(i ,j-1)+(0.75-ax*ax)*(0.5 -ay)**2*0.5
Q(i+1,j-1)=Q(i+1,j-1)+(0.5 +ax)**2*(0.5 -ay)**2*0.25
Q(i-1,j )=Q(i-1,j )+(0.5 -ax)**2*(0.75-ay*ay)
Q(i ,j )=Q(i ,j )+(0.75-ax*ax)*(0.75-ay*ay)
Q(i+1,j )=Q(i+1,j )+(0.5 +ax)**2*(0.75-ay*ay)
Q(i-1,j+1)=Q(i-1,j+1)+(0.5 -ax)**2*(0.5 +ay)**2*0.25
Q(i ,j+1)=Q(i ,j+1)+(0.75-ax*ax)*(0.5 +ay)**2*0.5
Q(i+1,j+1)=Q(i+1,j+1)+(0.5 -ax)**2*(0.5 +ay)**2*0.25
PIC法の数値ノイズ① 性質
• セルあたりの粒子数が多いほどノイズが低い
• 形状関数の次数が高いほど、ノイズの蓄積
が少ない
電磁場エネルギー
1st order
1/cell
4/cell
16/cell
64/cell
256/cell
1024/cell
時刻 (wpet)
2nd order
1/cell
4/cell
16/cell
64/cell
256/cell
1024/cell
時刻 (wpet)
PIC法の数値ノイズ② 原因
• エイリアジング
– 右に進んでたと思ってたら、
なぜか左にも進んでたらしい…
• 数値チェレンコフ放射（電磁場の数値分散
＋エイリアジング）
𝜔=0.8 ck
• 形状関数の次数を
上げる以外に、
今のところ解決策無し
w/wpe
𝑐𝑘
2
2
𝜔 =sin
+𝜔𝑝𝑒 2
2
k c/wpe
計算科学的観点からのPIC①
Why PIC ??
• 計算量を減らす：N体：N^2 ⇒ N+M (M:格子数)
• 電磁波の存在：粒子間のN体相互作用だけでは
電磁波の伝搬が解けない
– 電磁波の放射性（電磁場が広がる）より、格子で解い
たほうがよい
• 現在のPCでは数億粒子数（数GB）
• スパコンでは数百～千億個オーダー
（現在の世界最高は数十兆個オーダー）
Progress in particle simulations
(measured in terms of number of particles)
（中村琢磨氏
@LANL提供）
ゴードンベル
2008 Finalist
並列化手法は
門外不出・・・
- Present-day petascale supercomputers allow us to treat ~trillion particles.
 corresponding to ~50-100di3 scale 3D fully kinetic (PIC) simulations
計算科学的観点からのPIC②
性能を出しにくい…
⑭
④
③
⑤
⑨
⑪
⑫
⑱
②
⑦
⑮
⑬
①
⑰
⑥
⑲
⑧
⑯
⑳
⑩
電磁場
ex(i,j), ex(i+1,j), …
ex(i,j+1), ex(i+1,j+1), …
…
by(i,j), by(i+1,j), …
粒子
x(n), x(n+1), x(n+2), …
y(n), y(n+1), y(n+2), …
…
vz(n), vz(n+1), vz(n+2), …
• 格子データはメモリ上で連続
• 粒子データの物理的な座標はメモリ上の配置に無関係
• 粒子の座標から格子データを参照すると、ランダム
アクセスになる
PICコードの性能①
• 2次元、144,000,000粒子/ノード、2次精度
• カーネルごとに分けて測定
– Velocity：Buneman-Boris（ランダムロード
pex=(0.5-ax)**2*0.5*ex(i-1)+(0.75-ax*ax)*ex(i)
+(0.5+ax)**2*0.5*ex(i+1)
– Current：電荷保存法（ランダムロード＆ストア）
jz(i-1)=jz(i-1)+(0.5-ax)**2*0.5*qvz(n)
jz(i )=jz(i )+(0.75-ax*ax)
*qvz(n)
jz(i+1)=jz(i+1)+(0.5+ax)**2*0.5*qvz(n)
– Maxwell：FDTD法（2次中心差分、メモリバンド依存）
bz(i)=bz(i)+dt/dx*(ey(i)-ey(i-1))
– Etc：境界処理など
PICコードの性能①
• カーネルごとに分けて測定（10ステップ分）
FX10 – 16t
CX400 – 24t
+ Fujitsu C
CX400 – 24t
+ Intel C
Velocity
37.89 sec
9.517 sec
14.29 sec
Current
99.69 sec
40.93 sec
80.09 sec
Maxwell
3.02E-02 sec
1.99E-02 sec
3.02E-0.2 sec
Etc
0.8783 sec
0.4273 sec
0.7811 sec
– ノードの性能差： FX10:CX400~1:2
• 70%以上がCurrent, 20%以上がVelocity
以降、この２つのルーチンのみに注目
PICコードの性能②
FX10
Velocity
CX400+ifort
CX400+frt
Velocity
Current
FX10
CX400+ifort
CX400+frt
Current
FX10
CX400+frt
FX10
CX400+ifort
CX400+ifort
CX400+frt
PICコードの性能③
•
•
•
•
•
•
コア性能： FX10:CX400 = 1:1.48
Velocity 実効性能： FX10:CX400 = 1:2
Current 実効性能： FX10:CX400 = 1:2.5
スレッド性能：Fujitsu >> Intel
Velocity は割とスケール
Current のスケーラビリティは低い
– ランダムアクセスにより、隣のスレッドのアクセス
がキャッシュのデータを追い出す
粒子データの並べ替え
⑭
④
③
⑤
⑨
⑪
⑫
⑱
②
⑦
⑬
①
⑰
⑥
⑲
⑧
⑯
⑳
⑭
⑫
⑥
⑮
⑲
⑳
⑩
①
⑬
⑦
②
⑧
⑮
⑯
⑨
⑩
⑰
⑱
⑪
③
④
⑤
• 格子データのメモリ上の配置の順番に、粒子データを
並べ替える
ソート版PICコードの性能①
Velocity
ソート無
Velocity
ソート有
Current
ソート無
Current
ソート有
Etc(境界
＋ソート)
FX10 – 16t
CX400 – 24t
+ frt
CX400 – 24t
+ ifort
37.89 sec
9.517 sec
14.29 sec
10.91 sec
(x3.7)
99.69 sec
6.823 sec
(x1.4)
40.93 sec
9.369 sec
(x1.5)
80.09 sec
13.29 sec
(x7.7)
5.193 sec
(x7.9)
7.985 sec
(x10)
70.29 sec
8.233 sec
16.87 sec
ソート版PICコードの性能②
Velocity
FX10
CX400+ifort
CX400+frt
FX10
Velocity
CX400+frt
FX10
Current
CX400+frt
CX400+ifort
FX10
Current
CX400+frt
CX400+ifort
CX400+ifort
ソート版PICコードの性能③
FX10
CX400+frt
CX400+ifort
CX400+frt
FX10
CX400+ifort
• 自作のスレッド化バケツソートを使用。
• 計算アルゴリズムの一部（カウント作業）は、
currentのゼロ次版と同じ ⇒ ランダムアクセス
• ソートは基本的に整数計算 ⇒ FX10は苦手
並列版PICコードの性能④
•
•
•
•
コア性能： FX10:CX400 = 1:1.48
Velocity 実効性能： FX10:CX400 = 1:1.1
Current 実効性能： FX10:CX400 = 1:2
スレッド性能：Fujitsuはほぼスケール
Intelは高スレッド数で性能低下
• ソーティング（バケツソート）のスレッド性能が
低い・・・ ⇒コーディングに工夫が必要・・・
Particle-In-Cell計算の並列化の問題点

粒子分割（場は冗長に計算）：cannot
sustain large space domain and global
operations on it.
－スレッド並列はこれ。大容量メモリが必要な理由
－Reductionの負荷が無視できるぐらい粒子数が
多いとスケール



領域分割：cannot work when particles
are distributed non-uniformly.
動的領域分割： also has pitfalls for
particles populating a small subdomain
too densely.
OhHelp：省略(笑)
（三宅洋平氏@神戸大提供）
ノード間並列化の最近の話題
（昔はよかった編①）
• プラズマ粒子は磁力線に沿って動く
• 単位時間ステップで単位格子以上は動かない
一様プラズマ
磁力線がほぼ固定
⇒磁力線に垂直方向
に分割
B0
ノード間並列化の最近の話題
（昔はよかった編②）
y/li
• 問題に依存した分割
x/li
衝撃波：磁場と密度がほぼ比例
高密度領域が左へ伝搬
⇒衝撃波面方向に分割
電流層不安定性（リコネクション）
⇒電流層方向に分割
ノード間並列化の最近の話題
y/li
• 計算領域が増えたが故に・・・
y/li
x/ li
x /l i
最大で密度8倍差
ミニ磁気圏：密度100倍差以上
動的負荷分散①（藤本氏@天文台）
AMRだからうまくいっている？
28
動的負荷分散②（臼井@神戸大）
実際の仕事は沼波氏@核融合研
オーダリングは
毎ステップやってない
29
プラズマPICシミュレーションの
課題（まとめ）
• ランダムアクセス⇒粒子データの並べ替え
– （課題）並列ソーティングの性能
– （課題）粒子軌道の追跡
• 数値ノイズの軽減⇒高次形状関数
– エイリアジング・数値チェレンコフ放射の軽減
– B/F値の軽減
– 他の方法は？？
• ノード間並列時のロードバランス
– OhHelp
– 動的領域分割（D-1次元ソート + Morton order）

Download Report