プラズマシミュレーションにおける粒子法
梅田 隆行
名古屋大学太陽地球環境研究所

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  vn
 dt

その他（境界処理、ロードバランスetc）

t –Dt/2tJ 0t1 E  B  0 J  2c t  E   Bt

t +Dt/2
t +Dt

t +3Dt/2x 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/cell4/cell16/cell
64/cell
256/cell
1024/cell

時刻 (wpet)

2nd order
1/cell4/cell16/cell
64/cell
256/cell
1024/cell

時刻 (wpet)PIC法の数値ノイズ② 原因• エイリアジング– 右に進んでたと思ってたら、なぜか左にも進んでたらしい…

• 数値チェレンコフ放射（電磁場の数値分散＋エイリアジング）
𝜔=0.8 ck

• 形状関数の次数をw/wpe上げる以外に、
今のところ解決策無し

𝑐𝑘22
𝜔 =sin+𝜔𝑝𝑒 22k 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 CVelocity

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コードの性能②
FX10Velocity

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コードの性能③FX10CX400+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 forparticles 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）