J. Plasma Fusion Res. Vol.89, No.4 (2013)２
４５‐２６０講座 核融合プラズマシミュレーションの技法−大規模並列計算環境の活用−

５．粒子シミュレーションのコーディング技法
5. Coding Techniques of Particle Simulations
内 藤 裕 志，佐 竹 真 介１）
NAITOU Hiroshi and SATAKE Shinsuke１）
山口大学大学院理工学研究科，１）核融合科学研究所核融合理論シミュレーション研究系（原稿受付：２
０
１
３年３月１
５日）

この章ではプラズマ中の MHD 現象や輸送現象の解析に用いられる PIC 法とモンテカルロ法を題材に，粒子シミュレーションのコーディングにおける並列化，高速化のポイントを紹介する．PIC 法では領域分割と粒子分割のハイブリッド並列化，モンテカルロ法では乱数発生法の並列化とメモリへのランダムアクセスの抑制について主に解説する．
Keywords:
Particle-In-Cell, Monte-Carlo, parallel computing, domain decomposition, random numbers

５．
１ 序論

の新古典輸送を計算するために，クーロン衝突項を含むドプラズマ・核融合分野に応用される代表的な粒子シミュ

リフト運動論方程式を粒子シミュレーションとして解くた

レーション法は Particle-In-Cell 法（PIC 法）とモンテカル

めにも応用される．いずれの手法においても，衝突散乱過ロ法が挙げられる．両者とも，プラズマが荷電粒子の集ま程を模擬するために乱数を使うことからモンテカルロ法とりであるという粒子的描像を反映させて作られたシミュ呼ばれる．

レーションモデルであり，３章の MHD シミュレーショ

PIC 法とモンテカルロ法は粒子コードとして計算の並列

ン，４章の Vlasov シミュレーションで取り扱った３次元配化・高速化に関する技法に共通する部分も多いが，それぞ位空間の電磁流体，あるいは５次元位相空間（３次元配位れの手法の中でネックとなる部分には違いがある．PIC 法

空間＋２次元速度空間）における荷電粒子の分布関数を連では荷電粒子が作る電磁場の計算にかかるコストがモンテ

続体として扱う手法とは異なる数値計算上の工夫が必要とカルロ法に比べると大きく，粒子軌道計算の並列化と場のなる．

計算の並列化をどのように組み合わせて実装するかが大規

２つの粒子シミュレーション手法がプラズマ・核融合分

模並列化の鍵となる．５．
２節では３次元 PIC コードを例に，

野のどのような研究に応用されているかをみてみると，並列化コードのプログラミング法を，順を追って解説す

PIC 法は磁気リコネクション現象のように外部電磁場に加る．ま た Gpic-MHD コ ー ド で の 高 並 列 化 の 例 を 紹 介 すえてプラズマ自身が作り出す電磁場が強く影響する現象をる．一方，モンテカルロコードを並列化する上でまず必要調べる際や，電磁波と荷電粒子の相互作用，荷電粒子の集となるのは乱数の並列化であり，５．
３節でいくつかの方法

団的振る舞い，強いビーム成分の存在するプラズマの研究を紹介する．また，モンテカルロ法ではシミュレーションなどに用いられる．荷電粒子間の近接相互作用，いわゆる粒子の分布からアンサンブル平均量やモーメント量を取る

クーロン衝突は無視され，デバイ長以上のスケールでの電ことが多いが，その際にメモリへのランダムアクセスによ

荷の不均一による長距離相関に着目したモデルである．

る計算速度の劣化をどう防ぐかについても５．
３節で解説すPIC 法ではシミュレーション粒子１つ１つが荷電粒子を多る．

数集めた「超粒子」とみなされ，実際の荷電粒子の運動を

５．
２ PIC シミュレーション

少ないシミュレーション粒子で代表させて現象を模擬す
る．

プラズマは自由に動き回る電子とイオンの集合体であ一方モンテカルロシミュレーションはプラズマ・核融合

る．電子とイオンは外部磁場や外部電場に加えて，自分自分野において主に輸送現象の解析に応用されている．例え

身が作る電磁場の影響下で運動するため，きわめて複雑か

ば核融合反応で発生する高速アルファ粒子や周辺部からコつ神秘的な振る舞いを示す．この超多体的性質を持つプラアプラズマに流入する不純物イオンの輸送，その他荷電粒ズマを，コンピュータでシミュレーションする手法の一つ子以外にも中性子やプラズマ周辺領域の中性粒子がどう拡として粒子的手法がある．粒子的手法では，個々の荷電粒散するかを調べる際に用いられる．あるいはコアプラズマ

子の運動を支配するニュートン・ローレンツの運動方程式

authors’ address: Graduate School of Science and Engineering,Yamaguchi University, Ube 755-8611 Japan, 1)National Institute for Fusion Science, Toki
509-5292 Japane-mail: naitou@yamaguchi-u.ac.jp, satake@nifs.ac.jp

２４５!2013 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013

を，電磁場の時間・空間変化を支配するマクスウェルの方

間隔を %で表す．格子点上での，電荷密度，静電ポテン

程式と連立して解くことにより追跡する．粒子的手法は，基礎方程式は き わめ て シ ン プ ル で あ る が，巨 大 なコン

シャル，電場を表す配列 )&
，(&
，!&
#! '
#! '
#! 'を用意す#る．ここで ! は各格子点の座標を示す．以下では，場の量

ピュータ資源を要求するため，第一原理シミュレーションの（格子点間）補間や粒子の電荷の格子点への分配は線形の一つとして分類される．

補間を使用することとする．

粒子的手法を用いたコードは普通 PIC（Particle-In-Cell）

時 間 を１ス テ ッ プ（%,）進 め る た め に は 以 下 に 示 す

コードと呼ばれる．PICコードの詳しい解説は文献［１］に示

SOURCE，FIELD，PUSH のループを繰り返せばよい．空されている．PIC では，個々の粒子は空間を自由に飛び回

間 は %で 規 格 化 さ れ て い る と す る．特 定 の 粒 子

るのに対して，電場や磁場の場の量の計算には（メッシュ

#+&&
,
'
"
("
)'で代表
.
/
'
#&
&
+
&"
+
&"
+
&'は，以下に示すように，やグリッドと呼ばれる）空間格子を用いる．空間格子の最されるセルに属している．

小単位をセルと呼ぶ．１個の荷電粒子は空間内のどこかの- （!$&
-#"）
"&
+
&#'

１セル内に存在することになる．現在または近未来のコン.
. （!$&
.#"）
+
&#("&

ピュータでは，プラズマ中のすべての電子とイオンの振る/ （!$&
/
/#"）
+
&#)"&

舞いをシミュレーションすることは不可能である．PIC
コードで追跡する電子（イオン）は多数の電子（イオン）の，
ループの入り口では，すべての粒子の位置と速度 #+&&,'

電荷と質量を一つにまとめた超粒子であることに注意す

"+&&
,
!
%,
%
!!
'が与えられている．またループの入り口ではる．超粒子はその近傍の格子点上（線形補間の場合は超粒

毎回 )&#! 'にすべてゼロを代入する．

子の属するセルの各頂点）の電磁場の影響のみを受けて運１．SOURCE：各粒子の電荷を各粒子の属するセルの８個動する．
本節の構成は以下のとおりである．５．
２．
１では，PIC コーの頂点に分配し加算する（charge assignment）．加算する
配列は )&#! 'である．以下で矢印は代入を示す．

ドの基本構成を簡単に紹介する．５．
２．
２では，読者が実際にPIC コードを並列化する際の参考になるように，簡単な並

)&
'
"
("
)'

列化からより複雑な並列化に向かってコードを変更してい

'
"
("
)'
( )&
" *+&
-'
.'
/
"!&
"!&
"!&
&
&
'

くステップを提案する．５．
２．
３では，ジャイロ運動論に基づいた PIC コードである Gpic-MHD（Gyrokinetic PIC code

)&
'
"
("
)'
""

'
"
("
)'
( )&
""

)&
'
"
(""
"
)'

'
"
(""
"
)'
( )&

)&
'
"
(""
"
)'
""

'
"
(""
"
)'
( )&
""

)&
'
"
("
)""'

'
"
("
)""'
( )&

)&
'
"
("
)""'
""

'
"
("
)""'
( )&
""

)&
'
"
(""
"
)""'

'
"
(""
"
)""'
( )&

" *+&
-&
.'
/
"!&
"!&
&
'

for MHD simulation）の核融合科学研究所の Plasma Simulator（ 日立 SR 16000 ）と IFERC-CSC の HELIOS （ Bull 社

" *+&
-'
&
.&
/
"!&
"!&
'

Bullx B５１０）の両スーパーコンピュータ上での並列化の経験と現状について解説する．５．
２．
４では粒子コードの並列

" *+&
&
.&
/
"!&
'

化のまとめを述べる．
５．
２．
１ PIC コードの基本構成

" *+&
-'
.'
&
/
"!&
"!&
&

簡単のため３次元の静電近似コードを例として PIC コードの基本構成を概説する．直角座標中の各辺が "-，".，"/

" *+&
-&
.'
&
/
"!&

の直方体のシステム中に，電子（+#$）と１種類のイオン（+#%）からなる超粒子プラズマが存在している．+種の粒子の電荷を *+，質量を (+，個数を #+とする．時間 ,での

" *+&
-'
&
.
&
/
"!&
)&
'
"
(""
"
)""'( )&
'
"
(""
"
)""'
""
""

，"+&&
（&#"
粒子の位置と速度を #+&&"#+）とする．PIC
,
,
'
'

" *+&
&
.
&
/

コードは荷電粒子の運動を初期値問題として解くため時間差分を用いる．時間差分には，PIC コードの標準的手法で，すべての粒子に対して上記の処理を行った後，適当な定数

粒子の位置と速度が !
!
%,だけ異なった時間で計算される%をかければ，時間 ,での格子点上での電荷密度が計算され

蛙飛び法（leap-frog-method）を用いることにする．ここでたことになる．

%,は時間差分のステップ幅である．蛙飛び法では時間微分
２．FIELD：格子点上の電荷密度 )&#! 'から時間 ,での電が中心差分で近似されるように差分化を行っている．中心差分を用いると，実際の運動方程式と同様に，差分化した場 !&#! 'を計算する．通常は，ポアッソン方程式により運動方程式も時間反転に関して対称（time-reversal）にな#! 'から静電ポテンシャル(&#! 'を求め，(&#! 'の勾配)&

る．なお蛙飛び法は，１段法でありながら，２次の Rungeから !&#! 'を計算する．

$のオーダー）が保証さ-Kutta 法と同等の精度（誤差は&%,''%#(#)$!

れるアルゴリズムになっている．まず，すべての粒子の位

!#!%(

，"+&&
置と速度の初期値 #+&&
!
%,
!'
%
!!
'を与える．例えば位

３．PUSH：時間 ,での各粒子の属するセルの８個の頂点置は空間的に一様に，速度はガウス乱数を用いてマクス
ウェル分布に従うように設定する．場の量は空間的に離散の 電 場 !&#! 'か ら 個 々 の 荷 電 粒 子 の 位 置 で の 電 場化された等間隔の格子の格子点上でのみ計算される．格子

!&
#+&&
,
'
'を求め，
２４６

Lecture Note

5. Coding Techniques of Particle Simulations

H. Naitou and S. Satake

る場合は，Runge-Kutta 法や予測子・修正子法などの多段

"%
$'"%
"%
(
)&
*&
+
#
"
$"
%&
"!%
"!%
"!%
&
&
$%
%
%
&

法を用いる必要がある．５．
２．
３で紹介する Gpic-MHD コー

"%
"%
)%
*&
+
#
!
$!
%&
"!%
"!%
%
&
""

ドの場合は２次の Runge-Kutta 法（２段法）を用いている．

"%
"%
)&
%
*%
+
#
"
$""
"
%&
"!%
"!%
&

時間積分に多段法を使用する場合や，内挿公式に高次のス

"%
"%
%&
)
%
*%
+
#
"
$""
"
"!%
&
""

プラインを用いる場合でも，ここで述べた PIC コードの基

"%
)&
*&
%
+"%
#
"
$"
%""&
"!%
"!%
%

本構造は本質的に変わらない．なお PIC コードの基礎をも

"%
)%
*&
%
+"%
#
"
$"
%""&
"!%
""

う少し詳しく知りたい読者は，本学会誌の講座「プラズマ

"%
)&
%
*
%
+"%
#
"
$""
"
%""&
"!%

計算機シミュレーション入門」の中の「２．粒子シミュレー

"%
#
"
$""
"
%""&
)
%
*
%
+"%
""

ションの基礎」［２］を参照されたい．
５．
２．
２ PIC コードの並列化

この電場を使って，１ステップだけ荷電粒子の位置と速度

読者がこれから並列化 PIC コードを作成することを仮定を進める（particle pushing または particle acceleration）．

して，最短の並列化の順序を提案したい．以下に示すよう

#'"%
"%
$'"%
(
!
$(
(
!
$(
&'#
(
$(
$'&
$
$
"!
&
$#'"%
!!
&
"%
&
&

に小規模の並列化から大規模の並列化へとプログラムを変

$'"%
(
(
(
!
$(
$(
$
"$(
&
$$'"%
&
"#'"%
"!
&

更していく順序に従い解説する．

なお，磁場がある場合は次式を用いる．次式も時間に対す（１）並列化前の最適化

る中心差分を採用していて，時間反転に関して対称であ（２）スレッド並列化

る．

（３）プロセス並列化１（粒子分割）!

# $#'"%
"%
$'"%
(
!
$(
!
&'#
(
$(
$'&
$
$%
!!
&
"!
&
&

（４）プロセス並列化２（１次元領域分割＋粒子分割）!%$'"%#,$#!"#!#!
!
&'#
(
$(
$'&
$%
&
&

（５）プロセス並列化３（多次元領域分割＋粒子分割）

"

!

,

# $# "##%
!%
$'"%
&'#
(
$(
#
$'&
&
&

シミュレーションしたい物理現象によっては，小規模ま

#
$'"%
!
&'#
!%
(
$(
$'&
"")
!
$%
'
&
&
*
(

たは中規模の並列化で終了することも可である．また（２）

#'"%
"%
$'"%
(
!
$(
!
&'#
(
$(
$'&
$
$%
"!
&
$#""!
&
&
$'"%
(
(
(
!
$(
$(
$
"$(
&
$$'"%
&
"#'"%
"!
&

は場合によっては省略可能である．
（１）並列化前の最適化

すべての粒子に対してこの処理をした後，シミュレー

並列化前の PIC コードは存在して い る と す る．ま ず

ション領域（直方体）から出た荷電粒子の処理をする（粒は，この PIC コードを最適化しておくことが重要で あ子に関する境界条件の適用）．

る．基本的にはSOURCEとPUSHのサブルーチンを最適化すれば良い．実際には，最近のコンパイラは優れていて，，#'"%ループの出口では全粒子の $'"%
(
(
!
$(
$
"$(
&
"!
&が

よく知られている最適化手法を適用してもそれほどの差が求まっている．(+ ("$(としてこのループを所定のスでないことを認識される結果になる場合が多いと思われテップ数だけ繰り返せば PIC シミュレーションは終了する．ただし，グリッドサイズが大きい場合にはキャッシュる．

の関係で差がでてくるので注意する．このステップは PC

SOURCE では荷電粒子の電荷を空間的にランダムな格でも可能である．一見無駄なステップであるが，コードが子点に対応する配列の成分に加える scatter の計算を実行

並列化され複雑化する前に一度試しておくことは重要だとしている．また PUSH では，粒子の感じる電場は空間的に思われる．

ランダムな格子点に対応する配列の成分を引用するgather（２）スレッド並列化の計算を実行している．なお，荷電粒子が自分自身の作る

最初にスレッド並列化をする．OpenMP または自動並列

電磁場で自分自身を加速しないために，SOURCEでの電荷化を用いる．自動並列化は機種に依存するため可搬性を考の分配公式と PUSH での電場の内挿公式は同一にしておく

えるなら OpenMP で書いておいたほうがよいかもしれな

必要がある．PIC コードでは粒子に関する統計的精度を維い．OpenMP で 書 い て お い て も，コ ン パ イ ル の 際 に持するため，１セル内に数十から数百（場合によってはそ

OpenMP を使用しないで自動並列化を有効にすることもれ以上）の粒子を含めてシミュレーションを実行する．こ可能である．最近の PC では１チップ上に４コアが搭載さのため一般に，場の量の計算と比較して，粒子に関する計れているものがあり，intel 系のチップでは，１物理コアを算が支配的である．したがってSOURCEとPUSHの並列化

２論理コアとして動作させる hyper-threading 技術が利用

が PIC コードの並列化に特有の問題である．FIELD でのポできるため，８並列までのスレッド並列化が可能である．

アッソン方程式の解法では高速フーリエ変換か差分法を用

パブリックドメインのコンパイラ gfortran 等も OpenMP

いる．場の量の計算は一般の数値計算の教科書にあるものに対応している．この段階でFFTを使用している場合は並と本質的に変わらない．

列化 FFT を利用するように変更しておく．

電磁コードの場合には場の量の計算は少し複雑になる．

PUSH の計算は簡単にスレッド並列化できる．これは，また，ここでの解説で用いた蛙飛び法は，ニュートン・

同時に異なるスレッドから同じ配列の成分にアクセスがローレンツの運動法定式を時間積分するのに適した１段法あった場合も，配列の値の引用のみ で あ る た め，コ ンである．ドリフト近似を用いた運動方程式等を時間積分すピュータが適切に処理をするような設計になっていること
２４７

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013による．ただし，メモリの取り合いがスレッド並列化の効

レッド数が多くなるにつれて高速化が飽和しやすい傾向が果を限定的にすることに注意する．

ある．プロセス並列の場合は大きな粒度で並列化できるた

SOURCE では配列の成分へのランダムな加算計算を実め並列化性能を引き出しやすい．

行する必要があり，そのままでは並列化が困難である．こ（３）プロセス並列化１（粒子分割）

れは異なるスレッドから同時に同じ配列の成分に加算する

プロセス並列化を開始する．
（２）でスレッド並列化を済ま

アクセスがあった場合の動作が保証されていないことによしていると，ここからはスレッド並列とプロセス並列のハる．実 際，OpenMP で は，そ の ま ま で は 並 列 化 で き な

イブリッド並列になる．まず，第一段階として粒子分割

い．このため並列化のためのアルゴリズムを工夫する必要（particle decomposition）の手法により並列化する．粒子分がある．筆者はまだ試みていないが，場の量の複数のコ割の基本概念を図１に示す．MPI を用いる．場の量（電場ピー配列を用いてスレッド毎に計算した後，スレッド間のと磁場）にレプリカ（コピー）を用いて，各プロセスに対総和計算を行う［（３）の粒子分割に似た手法］とか，スレッ応させる．各プロセスは，レプリカとそのレプリカ内の粒

ド毎に場の量の異なる領域を担当する［（４）の領域分割に子の計算を行う．粒子の位置と速度の初期値は，プロセス似た手法］等のア ル ゴ リ ズ ム が 適 用 で き る と 思 われる

毎 に 変 え て お く 必 要 が あ る．レ プ リ カ 数 を !$&'*［３］．SOURCE 部に現れる演算はモンテカルロ法でもボト（"!'(&$%))）とすると，各プロセスは全体の粒子の !
!
!$&'*ルネックとなるものであり，スレッド並列化に関するいく

個の粒子を取り扱うため，粒子分割と呼ばれている．粒子

つかのアイデアが５．
３．
２節で解説されているので参考にさに関する演算は各プロセスで完全に独立である．プロセス

れたい．幸い Plasma Simulator では，HITACHI の自動並間の粒子の移動はない．PIC コードは基本的には粒子の計列化コンパイラが（実際どのような手法で並列化している

算が，場の量の計算より圧倒的に多い．このため，粒子分のが不明であるが）加算に対応しているので，DO ループ

割は非常に有効な並列化アルゴリズムである．
SOURCEでは各レプリカの電荷密度と電流密度（静電近の直前に似コードの場合は電荷密度のみ）を計算する．SOURCE*poption parallel, sum（加算する配列名）の直後に，全レプリカ間の総和を計算し，各プロセスに分という並列化制御行を入れておくだけでよい．なお自動並

配する総和計算（SUM）を実行すれば並列化は完了す

列化の場合は，かなり複雑なループも無理に並列化してしる．総和計算のアルゴリズムで一般によく知られているもまう危険性があるため，並列化で逆に遅くなっているルーのを紹介する．各プロセスは，それぞれのデータ（配列）を

プを探し出して，保持している．プロセス数を "! とする．プロセス数が８（!"#）の例を図２に示す．"番目のステップでは，全*poption noparallel

データを""個ずつのデータの組に分け，それぞれの組の中という並列化制御行を入れて，並列化を禁止することが重で ""!!離れたデータのペアを作る．ペア間でデータの送受要になる．

信を行い，２個のデータの和を計算する．! ステップで総本節の本筋からは離れるが，加算のスレッド並列化は

和計算は完了する．総和計算は，MPI の関数である MPIGPU（Graphic Processing Unit）を用いた並列化の場合に_ALLREDUCE を用いても計算できるが，読者には MPIでプも問題になる．GPU に含まれる多数の演算コアを，図形処

ロセス並列化する場合の練習問題として一度作成して実行

理 で は な く，一 般 の 計 算 に 利 用 す る こ と は，GPGPU

してみることを勧める．また，並列化したコードは将来総（General-purpose

和計算をカスタマイズしたい場合の基本コードとしても利computingongraphic

processing

用できる．

units）と呼ばれる．GPU を使用する場合は，各粒子を異な各プロセス（レプリカ）は，総和計算の結果の電荷密度るスレッドで計算するように並列化することができないた
め，セル毎にスレッド並列化するアルゴリズムが存在する［４］．各スレッドは各セルに対応してそのセル内の粒子の処理をする．この場合には PUSH の後で，粒子を各セルに再配置する処理が必要になる．領域分割の手法を１セルの極限にまで拡張し適用したような例になっている．
並列化が終了した段階で，スレッド数を変更して高速化
性能（計算時間の逆数）のスケーリングを確認しておくこと．ただし，スレッドあたりの粒子数は十分大きくなるよ
うにする．粒子数が十分大きくない場合はすぐにスケーリングが飽和してしまうので注意する．小規模の並列化で十分な場合はここで終了である．ただし，小規模の並列化の場合はこのステップを省略して次のステップでのプロセス
並列（粒子分割）のみを採用する選択もある．通常スレッド並列はDOループ単位の並列化であり，DOループが呼ば

図１

れる毎に，処理を異なるコアに分割して計算するため，ス
２４８粒子分割によるプロセス並列化の基本概念．

Lecture Note

5. Coding Techniques of Particle Simulations

H. Naitou and S. Satake

分割のプロセス並列化の概念を図３に示す．場の量に関する領域分割は３．
２節で述べられているものと基本的に同じである．隣り合う分割領域（subdomain）に共有される袖領域［ガードセル（guard cell）］の取り扱いに注意すればよい．
以下，すべての方向に高速フーリエ変換（FFT）を利用するとして解説する．３．
３節で述べられた転置（transpose）の手法を用いる．まず "，#方向について FFT（順 FFT）を実行する．$方向の FFT を実行する前に，分割の方向を例
えば "方向に転置しておく．転置後は，$方向の FFT の計算が異なるプロセスにまたがることはない．フーリエ空間
で場の量を求めたら，$方向の逆 FFT を実行し，その後で分割の方向を "方向から $方向になるように逆転置を行い，
図２

最後に "
，#方向の逆 FFT を実行すれば，実空間での場の総和計算のプロセス並列化のアルゴリズム．

量が求まる．
と電流密度を用いて，それぞれ電場と磁場を計算する．こ粒子コードの場合は，分割領域がそれぞれ多数の粒子をの部分は全プロセスでまったく同じ計算をすることになる含んでいる．このため，PUSH で１ステップ粒子の位置とので冗長計算になっている．

速度を進める毎に，異なる分割領域に出ていく粒子を行き先の分割領域に送り出し，異なる分割領域から入ってくる

粒子分割の１ループの手順をまとめると以下のようにな粒子を受け取る必要がある．このプロセスを MOVE で表す

る．ここで追加した部分を網掛で示している．
１．SOURCE

ことにする．いかなる粒子も隣接する分割領域にしか到達

２．SUM（総和計算）

できないことを仮定してプログラムするのが簡単である．

３．FIELD

このため，時間ステップ幅は最も高速の粒子でも，隣の分

４．PUSH

割領域にしか到達しないように選択しておく必要がある．

粒子並列の並列化性能を確認するために，各プロセスのプロセス間のデータ通信はパケット毎に行われるため，ま粒子数を固定して，プロセス数（＝レプリカの数）を変化

とまった量を同時に送受信する必要がある．このため，各させて計算時間のスケーリングを計測する．総和計算にか分割領域で，送り出す粒子と受け取る粒子を保存するバッ

かる時間は，レプリカ数を !#()+とすると，!#()+'(%"!#()+

ファー用配列を準備する必要がある．分割領域から出ていに比例する．総和計算を除いた部分は各プロセスで完全にく全粒子のデータを送信用バッファー配列（行き先は２方独立であるため，計算時間に変化がない．このため総和計向あるので各２配列を用意する必要がある）に格納し，隣算を除いた計算時間は，理論的にも実測でもプロセス数に接する分割領域に送信する．このとき受信も同時に行う．

関係なくまったく同一になる．総和計算にかかる計算時間

MPI _Isend/Irecv を用いた非同期送受信を行う方法はがその他の部分の計算時間と同程度になった時点で，粒子２．
４節で述べられている．受信用バッファー用配列（これも並列による性能向上は飽和することになる［粒子シミュ

２方向に対応して２配列を用意する必要がある）で受けレーションでは超粒子のノイズを低減して粒子数無限大取った粒子の情報を，各分割領域の粒子用配列に格納す（無衝突の極限）へのスケーリングを調べる必要がある場

る．各分割領域の粒子数は平均して一様になるが，分割領合がある．上記のスケーリングはこの場合のスケーリング

域間の粒子の移動があるため，各分割領域の粒子数は時間の導出の際にも必要になる］．ここでの並列化は各プロセスに対する負荷分布が一様である．粒子並列のプログラム
は基本的には総和のプロセスを導入するだけであるので簡単に短時間で並列化可能である．ここで並列化を終了しても，シミュレーションしたい物理現象によっては十分である．場の量の配列が大きく（場合によっては各コアのメモ
リに入りきれない場合を含む），場の量の計算時間が粒子を１ステップ進める時間と比較して無視できない場合は領域分割による並列化を採用する必要がある．
（４）プロセス並列化２（１次元領域分割＋粒子分割）
１次元の領域分割（domain decomposition）をコードに組み込む．ここでは場の量を $方向に，!$&*,個に等分割するとして説明する．簡単のため，"，#，$方向のメッシュ数!"，!#，!$は同じとする．まず，!#()+!!にして，純粋図３

な１次元領域分割コードを完成させる．純粋な１次元領域
２４９

領域分割によるプロセス並列化の基本概念．

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013

うになる．

的に変化する．このため，粒子に関する情報を格納する配列の大きさは十分余裕をもって設定しておく必要がある．

!$)+&"!-'+&#%!!*+)$&,,"!-'+&#%!!%(.0!!$)*/

この段階で，分割数に関するスケーリングをよく検証しておく．各分割領域の粒子数を固定して分割数 !%(.0を増やこの等間隔の領域分割は不安定性等が生じても分割領域

していくと，理想的には場の量の計算は !
!
!%(.0に比例して間で粒子密度の非一様性が大きくない物理現象の場合に有小さくなる．実際には袖領域や MOVE の計算に関するプロ効であることに注意する．ビーム不安定性のように粒子のセス間通信が加わる．FFT を使用する場合はフーリエ変換

バンチングによる大きな疎密が生じる場合には，分割領域

されたデータがモード毎に実部，虚部の順に入るため，境界の動的変更などの動的領域分割手法を取り入れる必要

フーリエ空間での計算が異なる分割領域（厳密にはフーリがある．

エ空間での分割領域になっている）に属しないためには，

領域分割と粒子分割を併用した場合の各ステップに対応

分割可能な分割数の最大値は !%(.0"!%!
"になる場合があするループの構造は以下のように表される．

ることに注意する．通常，分割数は最大値に固定してよい

１．SOURCEが，念のため計算時間の，分割数に対するスケーリングを

２．SUM（分割領域毎の総和計算＋ガードセルの処理）実測し，最速の分割数を確認しておくことが重要である．

３．FIELD（ ＋TRANSPOSE ）コードが正しく動いていることが検証できたら，総和計

４．PUSH

算の部分を分割領域に対応するようにコードを書き換えれ

５．MOVE（分割領域間の粒子の移動）

ばこの段階でのプログラムの変更は終了する．この場合は通常，１次元の領域分割数は使用可能な計算機のコア数

各分割領域が，それぞれのレプリカをもつことになる．領よりはるかに少ない．このため，領域分割と粒子分割を併域分割と粒子分割によるプロセス並列化の基本概念を図４用したハイブリッドプロセス並列を使用する必要がある．

に示す．ここでは，グループという概念を導入した方が説

このハイブリッドプロセス並列は，最近では domain clon-

明しやすい．グループ数はレプリカ数 !$)*/に等しい．それ

ing と呼ばれることがある．

ぞれのグループ内の粒子の運動はグループ内で閉じていこの段階で数千コアの高並列コンピュータの性能を十分る．このため MOVE の計算はグループ内の通信のみで達成

引きだせるコードになっていると思われる．この程度の並される．また TRANSPOSE と袖領域の計算もグリープ内で列化で十分予定していたシミュレーションが可能ならここ閉じている．分割領域毎の総和計算 SUM のみがグループで並列化を終了してよい．さらに大規模のシミュレーショ

をまたがった通信を必要とする．

ンを，数万コアから数十万コアの高並列コンピュータで実

ここで，!%(.0を最適値に固定し，プロセス毎の粒子数を行したい場合には，次に示す多次元の領域分割に進む必要

固定して，計算時間の !$)*/依存性を調べ，!$)*/をどの程がある．

度まで問題なく増やせるか確認しておく．!$)*/が大きくな（５）プロセス並列化３（多次元領域分割＋粒子分割）

ると，総和計算の時間が無視できなくなるので，!-'+&#%を

２次元の領域分割に対応したコードを作成する．１次元

増やした方が，計算時間が短縮される場合がある．このたの領域分割を２方向に繰り返すことで２次元領域分割版のめ，使用コア数を固定し，コピー数を増やした場合と，ス

コードを作成することができる．さらに必要ならもう１方

レッド数を増やした場合の最適値がどのようになるかを検向に領域分割し３次元の領域分割コードを作成する．ここ証しておく必要がある．ここで，全使用コア数は以下のよでそれぞれの方向の分割数を，!%(.0，!%(./，!"#$%とする

図４

領域分割と粒子分割によるプロセス並列化の基本概念．

２５０Lecture Note

5. Coding Techniques of Particle Simulationsと，３次元領域分割版での使用コア数は，H. Naitou and S. Satake

ル方向を示し，この方向には周期境界条件が適用されている（なお，Gpic-MHD コードの最終目標はトロイダル版で!%*,'"!-(,'$&!!&)./!!&).0!!&).1!!%*+0

ある．）基礎方程式系は文献［５］に示されている．荷電粒になるため，数万から数十万の高並列コンピュータの計算子（マーカー）の運動方程式はジャイロ運動論で表される．

機資源を十分利用可能になる．

これはジャイロ運動論 Vlasov 方程式の特性曲線に対応しここでの例では，直角座標系を採用しているので特に問ている．またデルタエフ（"）法を使用している．"""法では題はない．ただし複雑な座標系を使用している場合は，多の和とし分布関数 "
を熱平衡成分 ""!とそこからの変化分 "次元の領域分割は容易ではない．例えば円筒座標の場合をて表し，マーカー法を用いて ""の変化を計算する．各マー

考える．これは５．
２．
３で述べる Gpic-MHD コードの場合にカーは重み（weight）)'##($を持つ．マーカーの電荷と質対応している．軸方向は問題なく領域分割できる．方位角量は，&')'##，$')'##"法では，運動方程式に(($$になる．"

方向は軸上の特異点の存在のため領域分割の方向に採用す

加えて，重みの時間変化を表す式を同時に解く．

るのは難しい（軸上では粒子が任意の方位角方向分割領域

GYR3D の基礎方程式と，トカマクの鋸歯状振動の崩壊に移動する可能性があるため MOVE の計算が複雑にな過程に関連する無衝突内部キンクモードの非線形シミュる）．このため３次元の領域分割は実際的でない．２次元

レーションは文献［５，６］に示されている．これは PIC コー

目の領域分割の方向として半径方向を採用する．半径方向

ドで，無衝突磁気再結合を含む巨視的・運動論的 MHD 現の分割は，各分割領域に含まれる粒子数を一定にすると，

象（トカマクの内部崩壊現象に対応）を世界で初めてシ

非等間隔の分割になる．この場合，場の量の計算負荷は半ミュレーションしたものである．GYR3D の並列化は領域径方向の分割領域毎に一様でなくなる．このため，半径方分割と粒子分割のプロセス並列化に対応していた［６］．日

向の領域分割の効果は軸方向の領域分割の効果と比べると本語の解説としては，本学会誌の小特集「超並列計算機のすこし限定的である．軸方向と半径方向の２次元領域分割

シミュレーションとプラズマ」中の「１．プラズマの粒子の例を図５に示す．

シミュレーションコードの並列化」［７］を参照されたい．

５．
２．
３ Gpic-MHD コードの並列化

GYR3D は，直角座標で書かれていたため，シミュレーショ

ジャイロ運動論シミュレーションは第４章で紹介された

ンしたい物理に不必要な短波長で高周波数のアルベン波をように PIC 法と，Vlasov シミュレーション法がある．ここ含むため時間ステップ幅をきわめて小さくとる必要があ

では，PIC 法に基づくジャイロ運動論コードの一つである

り，計算機に対する負荷はきわめて大きいものであった．

Gpic-MHD コード（Gyrokinetic PIC code for MHD simula-

このため，円筒座標に対応した Gpic-MHD コードを開発しtion）の並列化について解説する．Gpic-MHD は運動論的効

た．円筒座標には，半径方向のみに依存する平衡解が存在

果を 含 み な が ら，流 体 系 の コ ー ド と比 較 し て「ク ロ ーすること，トカマク中の固有モードの表示に適した軸方向

ジャー」の問題がないため，流体コードの結果の正当性・と方位角方向のフーリエモード展開が利用できること等の健全性をチェックするためにも有用である．

利点ある．Gpic-MHD コードには，単一ヘリシティを仮定

Gpic-MHD は直角座標に対応した GYR３D コードを基礎した２次元版 Gpic-MHD と，マルチヘリシティに対応したとしている．どちらも，トカマクの第０次近似である直線

３次 元 版 Gpic-MHD が 存 在 す る．２次 元 版 Gpic-MHD

トカマクをシミュレーションしている．*方向がトロイダ

は，単一ヘリシティのモードのみを含むため，時間ステッ
プ幅を大きくとれること，基本的に２次元コードであることで，３次元コードと比較してコンピュータに対する負荷%""が小さく，特に $""!
（$，%はポロイダルとトロイ
ダルのモード数を表す）の内部キンクモードのシミュレー
ションに最適である．また３次元版 Gpic-MHD も不必要な
高周波モードをフーリエ空間で除去できるため，適切な時
間ステップ幅を選択することができる．
Gpic-MHD を Plasma Simulator にインストールし，並列化性能を検証した．ただしここで示す結果は，増強前のPlasma Simulator（Phase-1）を 使 用 し た 結 果 で あ る．
Plasma Simulator では１物理コアが２論理コアに対応している．２次元版２D-Gpic-MHD は，スレッド並列（自動並列）とプロセス並列（MPI）によるハイブリッド並列により並列化した．プロセス並列は粒子分割による．８１９２論理コアまで良好な並列化性能が得られることを実証した［８］．
３次元版３D-Gpic-MHD は，場の量の計算が増大するため，並列化のため領域分割を利用している．トロイダル方

図５

軸方向と半径方向の２次元領域分割．
（文献［１０］Fig.1 より転載）．

向のみの１次元領域分割と粒子分割を利用した場合の例［８，９］を以下に示す．なお，並列化はスレッド並列（%-(,'$&"#）
２５１

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013とプロセス並列（!./-)+00#!)-.2"!*,13）のハイブリッド並る．Plasma Simulator では１
６３８４コアまで，HELIOS では列になっている．メッシュは !""!!"!##"#(""#'""#'

３２７６８コアまでのスケーリングを検証しているが，特に顕とし，コア当たりの平均粒子数を百万個に固定し，１０００ス

著な飽和は観測されていない．

テップの計算を実行した．!*,13#&$とし，!)-.2を１から標準的なジャイロ運動論的 PIC コードのアルゴリズムに

６４（!)-/+#"#'!'"(#）に変更して FLOPS 値の論理コア数

加えて先進的アルゴリズムの開発研究も行った．標準的な依存性を調べたものを図６に示す．比較的良好なスケーリ

ジャイロ運動論的 PIC コードは，大規模・高ベータのプラ

ングがPlasma Simulator の最大コア数まで得られることズマを取り扱う場合，電流密度から磁場を求める際に大きを実 証 し た．な お，８１９２論 理 コ ア を 使 用 し て 得 ら れ たな誤差が生じることが知られている（「キャンセレーショ

2.4 TFLOPS は Plasma Simulator の理論最大値 77 TFLOPSン」の問題）．この問題の解決法として拡張 split-weight-の３．
１％である．

scheme があるが，我々は場の量の計算に渦方程式と磁力

このスケーリングを主なループの構成要素に分けて分析

線方向の一般化オームの法則を用いる方法を提唱し，このし た．プ ロ セ ス 間 通 信 を 必 要 と し な い PUSH，方法による新しい 2D-Gpic-MHD を作成した．また，新しい

SOURCE，FFT（FIELD中のFFTだけを取り出したもの）に

2D-Gpic-MHD を用いて大規模・高ベータの領域でも精度関する計算時間は論理コア数に依存せずほとんど一定である こ と が 実 証 さ れ た（図７参 照）．通 信 に 関 す る 部 分
MOVE，TRANSPOSE，SUM の計算時間の論理コア数依存性を図８に示す．SUM の計算時間が !)-.2の増加とともに増大するのは自然である．TRANPOSE と MOVE の計算はグループ内にとどまるため，もし各グループで通信回線が独立であれば，計算時間は一定のはずである．このため，
TRANPOSE と MOVE の計算時間の，!)-.2の増加による増大は，通信回線が共通であることによる通信容量の制限に起因すると思われる．
トカマクの計算では有理面近傍の微細構造を解像するためにも，半径方向の解像度を増やす必要がある．半径方向の メ ッ シ ュ 数 を 増 や し，!""!!"!##"!#%""#'""#'
として，計算時間を測定した．トロイダル方向と半径方向の２次元領域分割の場合１次元の領域分割より高速化されていることを実測した．これは，場の量の計算時間が半径方向の領域分割により短縮されていることによる．２次元の領域分割についても論理コア数に対する良好なスケーリ

図７

通信に無関係な部分の計算時間の論理コア数依存性．
（文献［９］
Fig.3 より転載）

図８

通信に関系する部分の計算時間の論理コア数依存性．
（文献［９］
Fig.4 より転載）ングを得た［９，１０］．増強後（Phase-2）の Plasma Simulatorと HELIOS でのスケーリングについては現在検証中であ図６

FLOPS 値の論理コア数依存性．
（文献［９］
Fig.2 より転載）．

２５２

Lecture Note

5. Coding Techniques of Particle Simulations

H. Naitou and S. Satake""!!の運動論的内部キンクモードがシミュ
よく !!!

もアルゴリズムを変える必要がなく，収束速度の問題を別

レーションできることを実証した［１１］．

にすれば決定論的なアルゴリズムに比べて多次元化が容易

今後のスケジュールとしては，新しいアルゴリズムを採であるという利点をもつ．また，扱う問 題 に 系 か ら の用 し た３D-Gpic-MHD の 作 成 と ト ロ イ ダ ル 版３D-Gpic-

フィードバックのない場合（例えば時間的に固定された場

MHD の作成がある．

の中でのテスト粒子の拡散問題や，モンテカルロ積分法な

５．
２．
４ PIC コードの並列化のまとめど），解いている問題を独立事象の和とみなすことができ，非常に並列化効率の高い計算が実行可能である．

５．
２節では PIC コードの並列化の手順について段階的に解説するとともに，Gpic-MHD コードの高並列コンピュー

本節で取り上げるのは５．
１節で紹介したモンテカルロ法

タでの並列化の現況について報告した．Gpic-MHD の並列のプラズマ・核融合分野への応用例のうち，特にプラズマ化に関しては，（１）プロセス並列における２次元の領域分

中の輸送現象をドリフト運動論あるいはジャイロ運動論モ

割の場合の最適化と（２）スレッド並列のさらなる最適化がデルで，かつ粒子シミュレーションとして扱う事例を念頭課題とし て 残 さ れた．Gpic-MHD コ ー ド の 並 列 化では，においた並列化コーディング技法であるが，乱数の並列化

Plasma Simulator で自動並列化コンパイラが利用できたたに関して（５．
３．
１）はテスト粒子的モンテカルロ法や熱浴モめ，特別なアルゴリズムは採用しなかった．高度のチュー

ンテカルロ法，モンテカルロ積分法など，統計的サンプリ

ニングが可能な OpenMP を利用し，プロセス内のスレッド

ングを多ノード数で並列実行する場合にも役に立つであろ間での粒子分割や領域分割を取り入れた並列化は今後の課

う．またランダムアクセスの抑制（５．
３．
２）や MPI_REDUCE_SCATTER を用いた並列化の応用例（５．
３．
３）については，

題である．

PIC 法など他の粒子コードを記述する上でもヒントになる

超並列計算になればなるほど，個々のコンピュータに対する最適化が必要になる．異なるコンピュータの個性もあと思われる．

るので，高度の並列化はきわめて専門色の強いものにな

５．
３．
１ 乱数の並列化

る．計算物理の研究者は機種依存性の少ない，物理モデル

粒子シミュレーションモデルのドリフト運動論やジャイと並列アルゴリズムの開発に専念して，機種依存性の強いロ運動論におけるモンテカルロ法とは，フォッカー・プラ

最適化は，専門家にまかせるのが賢明であると思われる．

ンク方程式の形で記述されたクーロン相互作用による２体

なお，次世代の高並列コンピュータでは，GPU のような衝突項を乱数によって模擬することであり，具体的には

計算アクセラレータが搭載される可能性が高い．マルチコ

マーカー粒子の速度 %を，ある確率密度分布 #"
#%
!
&
!
%
!
$
#

アやメニーコアに対するスレッド並列化に加えて，計算アに従うランダムな #%だけ変位させることである．ここで，クセラレータによるスレッド並列化のための準備も必要と#は一般的に各マーカーの位置 &と速度 %，および時刻 $感じる．

におけるプラズマの背景パラメータ（密度，温度など）に

PIC コードの並列化については，プラズマ・核融合学会

依存する．クーロン衝突項をモンテカルロ法で模擬する具主催の第２２回専門講習会「プラズマ・核融合分野での計算体的な方法については［１４‐１６］を参照するとよい．

機シミュレーション技法とその応用」（２００９年１２月１６日，東

さて，このようなモンテカルロシミュレーションを並列

京大学山上会館）で「PC クラスタの原理と基礎」という題

計算機で実行する場合，乱数をどのように生成すればよい目で講演した．特に２次元粒子コードについて解説し，プか考えよう．今，MPI 並列数が nmpi で，1MPI プロセスあ

ログラムの一部も掲載しているので参考にされたい［１２］．

たり np 個のマーカーを用い，５．
２節の PIC 法と同様に場のまた，２次元の静電近似 PIC コードの非並列版と粒子並列量のレプリカを nmpi 個使った並列化を行うとする（総化版のコードは学会 HP［１３］からダウンロードできるのでマーカー数は nmpi×np 個）．乱数は衝突項の計算で #%参考にされたい．

を与える場合だけでなく，マーカーの初期分布をある種の統計分布（マクスウェル分布など）に従って発生させる場
合にも必要となる．任意の確率密度分布 #に従う #%は，

５．
３ モンテカルロシミュレーション
一口にモンテカルロシミュレーションと言っても，そこ［０，１］区間の一様乱数をもとに直接法，棄却法，合成法には広範囲な研究対象と実装方法がある．何らかの物理現［１７］などを用いて生成することができるから，各 MPI プ

象（あるいは社会現象や生態系などでもよい）のシミュ

ロセスで必要なのは［０，１］区間で発生させた np 個の乱数でレーションをする際に，その素過程が確率論的な現象であある．ただし，半導体の熱雑音などを利用した物理乱数

るとみなせる場合，あるいは素過程を方程式として記述す［１８］ではない，数値計算ライブラリなどに用意されているることができないか困難である場合に，ある確率分布を

乱数は正確には疑似乱数と呼ばれるものであり，一般的に持って振る舞うモデルとして近似できる場合に用いられる

ある種の漸化式を解くことによって統計的に十分乱数のよのがモンテカルロ法である．さらに，ある型の偏微分方程

うに振る舞う＊１数列を決定論的プロセスから生成するもの

式の境界値問題や，複雑な境界形状をもつ積分区間での数である．数値ライブラリ毎に形は異なるが，非並列コード

値積分の技法として用いられることもある．一般的に，モの場合以下のような手順で呼び出される．

ンテカルロ法を用いた数値解法は扱う問題が多次元化して＊１ 疑似乱数がどれだけ「乱数っぽい」のか，またそれをどう検定するのという議論は，それだけでかなりの研究がなされている
奥の深い研究テーマである．興味のある読者は
［１
７，
１
９，
２
０］
やそこに紹介されている文献を参照することをお勧めする．

２５３

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013

（サンプルコードも間違っている．なぜなら，一般的に乱数の初期化ルーチ

5.A）

ンは引数 id を基に内部で別の乱数ルーチンを使って iwk

integer :: i,iwk(nwk),id,npの初期値を決めており，idの値に大きな違いがあっても乱real :: rnd(np),dv(np),v(np)数列の始点位置が大きく離れることは保証されないからで

id=123456

ある．実際上は，このような間違ったやり方で乱数発生を

call init_rnd(iwk,id)

並列化しても目立った問題が起こらないかも知れない．し! 乱数発生ルーチンの初期化

かしそれはたまたまあってはならない相関関係が目立たず

call gen_rand01(rnd,iwk,np)にシミュレーションが実行できたというだけで，並列数! [0,1]区間の乱数列を配列 rnd に入れる

nmpi を増やした時に常に正しく動く保証はない．

call rnd_to_dv(rnd,dv,np)

それでは，並列コードの中で乱数はどのように呼び出せ! rnd を用いてある確率分布に従う dv を求める

ばよいかというと，次のような方法が挙げられる．

do i=1,np

１．単一の乱数列を，MPI プロセス毎に重複しない区間__v(i)=v(i)+dv(i)を選んで利用する．

end do２．互 い に 独 立 で あ る と い う 保 証 の あ る 乱 数 列 を，

ここで，iwk は乱数発生のための作業配列であり，その大nmpi 個用意して利用する．

きさ nwk は利用する疑似乱数ルーチンに依存する．iwk

３．物理乱数を利用する．

は疑似乱数を計算するのに必要な漸化式の履歴が保存されている「状態行列」であり，漸化式の続きを計算するため

まず，１．については前提として周期長の非常に長い疑似にその値を保持しておく必要がある．init_rnd は iwk

乱数ルーチンが必要であり，その実装の仕方として２通りの初期値を与えるもので，計算の最初に１回呼び出せばよある．１つ目は，各 MPI

い．また，途中で止めた計算を再開する際には，init_

いに重ならないように飛び飛びに読んでいくもので，以下プロセスが同一の乱数列から互rnd を呼ぶ代わりに前回終了時にファイルに記録しておいのようになる．
（サンプルでは mpi_init，mpi_rank な

た iwk を読み込めばよい．init_rnd の引数 id の値は乱ど MPI コードの実行に必要な部分は省略してある．）

数の初期値を決める．ここで注意すべき点は，疑似乱数列（サンプルコードはそれぞれ固有の周期長をもっており，idはその数列のど

5.B）

integer :: i,iwk(nwk),id,np,myrank,nmpi

こから乱数を読み始めるかを指定しているに過ぎないとい! myrank=[0,1,2,…,nmpi-1]

うことである．衝突による速度変化 !!を乱数で模擬する

real :: rnd0(nmpi*np),rnd(np),dv(np),v(np)ためには，疑似乱数列 rnd が相関を持たない白色ノイズのid=123456

ように振る舞う必要がある．実際には決定論的に生成され! 注：この値は各 rank で共通のものを使う.

ているので完全に白色ノイズというわけではないが，メ

call init_rnd(iwk,id)

ジャーな数値ライブラリで使われている周期長の長い（現call gen_rand01(rnd0,iwk,nmpi*np)

在では #"!! 以上のものが一般的である）に関しては，疑似

do i=1,np

乱数列全体に渡って統計的によく乱数のように振る舞うと__rnd(i)= rnd(i)=rnd0(np*myrank+i) !

考えてよい．ただし，だからといって次のようにモンテカ

end doルロコードを MPI 並列化することは間違いである．

call rnd_to_dv(rnd,dv,np)…

１．乱数列初期化の id を，MPI プロセスごとに適当にずらす．例：id=myrank*100+1

注意点として，id と np の値は各 rank で同一にしなけれ

２．他の乱数ルーチンを用いて，各 MPI プロセス毎にばいけない．図９にこの例の模式図を示す．ただしこの方id を乱数で決める．

法は各 rank で必要な乱数より余分な乱数列を発生させて３．１．や２．だけでは心配なので，そのように各 MPIいることにな り，乱数生成にかかる時間は非並列時のプロセスでバラバラの id の値を使って初期化した

nmpi 倍になる．並列数が多くなった場合，どんなに乱数生

２つの乱数ルーチンの平均値を採用する．

成ルーチンが速くても計算コストが無視できなくなる可能

なぜならば，単に id の値を変えるということは乱数列の読み取りの始点を動かすことに過ぎず，ある MPI プロセスで j 番目に読まれた乱数は他のプロセスにおいて j+k 番目に読まれる値と同一であるからである．もし k の値が偶然小さかった場合，この２つの MPI プロセスには強い相関関係があるということになり，MPI 並列化されたコードの各プロセスで計算される衝突項が独立事象であることが保障されないのである．「周期長が長い乱数列だから id をプ

図９

ロセス毎に大きく離せば重ならないのでは」という考え方
２５４

サンプルコード 5.B における MPI 並列コードにおける乱数の並列化例．nmpi=4 の場合を示している．

Lecture Note

5. Coding Techniques of Particle Simulations

H. Naitou and S. Satake性がある．逆にいえば，マーカー初期分布の生成を並列化とまったく同じでよい．現在の計算機の性能であれば npする場合など開始時に１回使うだけなら，この方法で十分*nmpi*nmax の値がかなり大きくても，この事前準備にかかもしれない．

かる時間は大したものではないであろう．メインのシミュ

重複しないように乱数列を読むもう一つの方法は，あらレーションにおいては乱数生成にかかる時間は非並列の場

かじめ非常に遠く離れた位置の状態行列 iwk の値を nmpi

合と変わらず，理想的な並列実行性能向上が見込める．も通り用意しておき，それをシミュレーションの開始時に読し同じパラメータのシミュレーションを（統計性を見るた

み込むことである．図１０に模式図を示す．今，各 MPI プロ

めに）乱数だけ変えて複数回試行したい場合は，シミュセスで np 個のマーカーを用いたモンテカルロ計算で，充レーションの最初に十分長い乱数発生の「から回し」を，分多く見積もって nmax 回衝突計算などに乱数を呼び出す

試行ごとに長さを変えて行えばよい．

可能性があるとする．この場合，各 rank の乱数列の呼び

疑似乱数ルーチンは漸化式を解いているため，!!!回

出しの始点を np*nmax ずつずらしておけば，シミュレー

先の状態行列を求めるには基本的には 5.C のように N 回実

ションの中でお互いに重複する乱数列を使う事は起こら

際にから回ししてやる必要があるが，最近の疑似乱数生成

ず，互いに独立な計算をしていることになる．そのような法の中には１つずつ逐次実行せずともある整数 n 回分先のiwk の初期値は，次のような非 MPI コードを事前に回して状態行列にジャンプできるものがある．そのような疑似乱

求めておく．

数ルーチンを利用できる環境であれば，5.C のから回し（サンプルコード

ループを置き換えることができるし，またその実行速度が5.C）

十分速いようであれば 5.C に相当する部分をシミュレー

integer, parameter :: nmpi=***,np=***,&&

ションコードの先頭に組み込んで，計算開始時に各ランク

nmax=***,nwk=***毎に十分離れた乱数列の始点を決めればよい．

integer :: iwk0(nwk,nmpi), iwk(nwk),&

次に，並列化の方法２．で挙げた，独立性の保証された並& i,id,k,m

列 乱 数 列 の 利 用 と し て，疑 似 乱 数 ル ー チ ン Mersenne

real :: rnd(nmpi*np)
id=123456
call init_rnd(iwk,id)

Twister（MT）［１９］とその並列化のための Dynamic Ceartor（DC）を紹介する．これらは松本眞氏（東大数理科学研究科）によって開発されたもので，ソースコードがオープン__do m=1,nmpi____do k=1,nmax ! 乱数発生のから回しループ______call gen_rand01(rnd,iwk,np)になっているので利用したい方はぜひ参照してほしい［２１］．Mersenne Twister は単体の疑似乱数としても良好な統計的性質と非常に長い周期長などの特長を持ってい____end do

る．Dynamic Creator は大まかに 説 明 す る と，こ の MT____iwk0(:,m)=iwk(:)

を規定する漸化式の係数を，互いの乱数列が相関をもたな__end do

いように複数組生成してくれるものである．すなわち__write(10) iwk0

nmpi 並列でモンテカルロシミュレーションする際に，

end

nmpi 個の「異なる乱数列表」を用意してくれるというわけ

そして，メインプログラム実行時に init_rnd を呼ぶ代わである．公開されている DC のソースコードは C 言語で書りにかれているため FORTRAN ユーザーにはわかりにくいかもしれないが，事前準備として get_mt_parameter_id

read(10) iwk0に異なる整数値 id を入力として与えて得られた構造体 mt

iwk(:)=iwk0(:,myrank+1)_struct の３つ組（mts->aaa, mts->maskB, mts->

として iwk の初期値を各 rank に与えればよい．その後のmaskC）の値をファイルに記録しておく．MPI コード内で乱数ルーチンの呼び出し方は非 MPI のサンプルコード 5.A

MT を並列に利用する際にはそのファイルを読み込み，
MPI ランク毎に異なる（aaa, maskB, maskC）の組を MT
サブルーチン内のパラメータとして与えればよい（ネット上 に 公 開 さ れ て い る 作 者 の MT の ソ ー ス コ ー ド
genmtrand.c を参照せよ＊２）．ただしこの DC は実行に時間がかかるので，数千並列を超える計算の準備にはそれなりの日数が必要となるだろう．なお，最近 MT の状態行列
を""回分ジャンプさせる手法が発表されたので，これを用いて１．の方法で MT の並列化を実現するのもよい．それ

図１０ サンプルコード 5.C における MPI 並列コードにおける乱数の並列化例．

に関する情報やサンプルコードが公開されている［２２］．
最後に，３．の手法については物理乱数ボードが組み込

＊２ ［２
１］
のホームページ中に，有志の方が作られた Fortran９０の MT ソースコードも公開されている．ただし，DC で生成した３つ組変数を利用する際，Fortran９
０サンプルコード中にある TSHFTU(y)=ishft(y,-11)は DC のソースコードで使われている shift_00(=12)の値に合わせ，TSHFTU(y)=ishft(y,-12)に修正する必要がある．

２５５

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013

まれた計算機が利用できる環境がないと使えないので，詳て，いくつかの方法が考えられる．なお，ここでは単純化

しい説明は省く．物理乱数の利点は，どんな疑似乱数よりのため１次元空間の分割を考えているが，５．
３．
３節のコー

も乱数らしい乱数である，という点に尽きる．その一方で，ド 5.H のように多次元空間における空間セルインデックス

再現性がないためコードのデバッグ用途には使えないこの１次元化を行えば，本節で示す考え方は５．
２．
１節で解説

と，以前に比べると高速化されたとはいえ，乱数列の生成

したような PIC 法における SOURCE プロセスのスレッド並

速度はまだ疑似乱数に及ばない点が挙げられる．疑似乱数列にも応用可能である．

の品質に不安や関心のあるユーザーは利用を試みるのもよ

まず一つ目は，各区間 k に入るマーカーのインデックス

いだろう．MPI 並列コードから利用する際は，各 MPI プロと総数を予め調べて記録しておく方法である．

セスが走っているノードに乱数発生器が付いているか否か（サンプルコードの対応関係が問題になり，例えば発生器を持つノードで実5.E）__nb(:)=0

行されているランクが他のランクに乱数列を MPI_SCAT-__do i=1,np

TER するなどの工夫が必要になる．なお，統計数理研究所____k=int(x(i)/dx)+1

ではインターネットを通じていくつかの物理乱数器や上記____nb(k)=nb(k)+1の MT で発生させた乱数（３０分毎に更新）をダウンロード____ib(nb(k),k)=iできるサービスを提供している［２３］．

! 区間 k に入る nb(k)番目のマーカーのインデックス

５．
３．
２ 粒子コードにおけるランダムアクセスの抑制__end do

この節では一旦 MPI による分散並列の話は置いて，１!$OMP PARALLEL DO

ノード内の並列メモリ上でのスレッド並列化に目を向けよ__do k=1,nx

う．粒子コードの場合，２．
３章で説明したように OpenMP____u(k)=0

指示行を用いてマーカーのインデックスに関するループ____do j=1, nb(k)（例：サンプルコード 5.A 中の do i＝1,np）をスレッド並______i= ib(j,k)

列化するのが最も容易だ．しかし，マーカーは系の中をバ______u(k)=u(k)+v(i)

ラバラに走り回るので，その位置はインデックス i に関し____end doて不連続になる．モンテカルロ法や PIC では系を微小区間____u(k)=u(k)/nb(k)に切り，各区間にいるマーカーのアンサンブル平均量（密__end do

度や平均速度など）を計算することが多い．例えば，マー!$OMP END PARALLEL DO

カーの位置 x が !!"!"!の範囲にある場合，これを nx
個の区間に分割して平均速度 u を求める計算は次のように

このようにすると２番目のループは区間 k に対するスレッなる．

ド並列になり，コード 5.D より効率がよくなると期待でき（サンプルコード

る．ただし，各区間に入るマーカー数 nb(k)がなるべく均

5.D）

等であることが望ましい．ところで，１番目のループは順__dx=xb/nx; u(:)=0; nb(:)=0

序に依存するので，このままではスレッド並列化できな!$OMP PARALLEL DO reduction(+:u,nb)

い．このようなループも並列化するには，以下のような明__do i=1,np

示的な並列化が必要である．

____k=int(x(i)/dx)+1____u(k)=u(k)+v(i)（サンプルコード____nb(k)=nb(k)+1__integer, allocatable :: kst(:),ked(:)__end do__nomp=OMP_GET_MAX_THREADS()!$OMP END PARALLEL DO

5.F）!スレッド並列数を取得する関数! ここで各スレッドの部分和が u,nb に集められる.

__allocate(kst(0:nomp-1),ked(0:nomp-1))__u(:)=u(:)/nb(:)__do ith=0,nomp-1____kst(ith)=ith*nx/nomp+1

これで u と nb に関する総和演算がマーカーインデックス____ked(ith)=(ith+1)*nx/nomp

i に関してスレッド並列化される．各スレッドはそれぞれ__end doが一時的に u と nb の部分和を記憶するための作業配列を__nb(:)=0

用意し，doループの終了後に全スレッド分の総和を集計す!$OMP PARALLEL DO private(ith,i,k) &る．し か し 多 次 元 領 域 分 割 の 場 合 や 分 割 数 が 多 い と，&__shared (nb,ib)

キャッシュに作業配列が乗り切らなくなったり，最後に各__do ith=0,nomp-1

スレッドの部分和を u や nb にかき集める作業に時間がか____do i=1,np

かったりして並列実行性能が出ない場合がある．また，各______k=int(x(i)/dx)+1

マーカーがどの区間に入るかを示す k の値が，i に対して______if((k<kst(ith)).or.(k>ked(ith))&

不連続に変化することも，このような区間ごとの統計量の______& cycle

計算を遅くする要因である．このような場合の対処法とし
２５６

Lecture Note

5. Coding Techniques of Particle Simulations

H. Naitou and S. Satake______nb(k)=nb(k)+1

こす恐れがある．このようなケースはどうチューニングす______ib(nb(k),k)=iればよいであろうか．トーラスプラズマの新古典輸送コー____end do

ドを一例として取り上げる．
新古典輸送コードでは，磁気座標系 '
+
!
%
!
&
(で与えられ__end do

た三次元磁場中の案内中心軌道を解くために，i 番目の!$OMP END PARALLEL DO

マーカー位置'+!
%!
&(*における磁場の値 "*やその微分値がこの例では，ith のループに対してスレッド並列化され必要である．ところでトーラスプラズマ中の荷電粒子は磁る．第 ith スレッドは kst(ith)<=k<=ked(ith)の区間

力線に沿って速く動くが，磁気面を横切る +方向のドリフに入るマーカーの数とインデックスだけをみて，その情報

ト運動は比較的遅い．したがって，磁場データを '
+
!
%
!
&
(

を共有メモリ上の配列 nb，ib の各スレッドが担当していの３次元メッシュデータではなく，+方向にスプライン展

る区間に直接書き込む（shared 指定のため）．担当区間に

開，'
%
!
&
(方向にはフーリエ級数展開した次のような形で重複がないので，配列の同じ要素に複数のスレッドが同時与える（注：sin 成分をもつ場合もある）．

に書き込むといったエラーは起きない．

#(

上のようにローカル区間ごとのマーカーのインデックス

"'
+
!
%
!
&
/
(
.0)
(
,'
,(
-'
,(
(
#%*!(
#&
*% ! ", '
*!/
**
+(

を調べる方法は，そのインデックスを何回も使いまわせる(%"

場合にはメリットが大きいが，もし知りたい統計量がコー

ただし，+方向のメッシュ点数を #+，+の最大値を a とし

ド5.Dのように１つだけなら，サンプル5.Fの$OMP PARAL-て，)/%'，+%*，/#&/#)/とする．したがっ
#
-1
/
)/
'
(
*
+%+

LEL DO 以下を次のように書いてもよい．
（サンプルコードて !&%,+であり，"( '
),
*%+
*!+
*(はスプライン補間+")係数によって5.F'）__nb(:)=0; u(:)=0

$
#
", '
),
#'!( $),$'!( $),%$!'!( "$"'!( $),
*(
*"$
* "$
&

!$OMP PARALLEL DO private(ith,i,k) &___& shared (u,nb)の形で与えられる．この記述の利点の一つは，場の値だけ__do ith=0,nomp-1

でなくその '
‐方向微分も容易に計算できることであ
+
!
%
!
&
(

____do i=1,np

る．ところで，軌道を解く毎に各マーカーがその瞬間に滞______k=int(x(i)/dx)+1

在している区間 l に対応したスプラ イ ン 補 間 テ ー ブ ル______if((k<kst(ith)).or.(k>ked(ith))&#)*個のデータ
$'
!!$(
'
!
( を参照していたのでは，#(#%______& cycle

を毎回読み込むことになる．一般的に磁気面が複雑な形状______nb(k)=nb(k)+1をしていたり，プラズマ $値が高くシャフラノフシフトが______u(k)=u(k)+v(i)大きいほど磁場の'‐依存性を表現するのに必要なモー%!
&(____end do

ド数 Nm を大きくとる必要がある． データ参照を必要最__end do

低限に抑えるために，１ステップ前の各マーカー位置の l!$OMP END PARALLEL DOの値とそこでのスプライン係数を配列 li，ci に保存して__u(:)=u(:)/nb(:)

おき，以下のようにして磁場データを参照する．

この他にも総和演算に対する明示的なスレッド並列化の方（サンプルコード法として，複数の統計量（u の他に v の n 乗モーメントのinteger :: li(np)5.G）

和など）を求める必要がある場合，スレッドごとに担当す

real:: r(np),dsi(np),ci(0:3,Nm,np),&

る統計量を明示的に割り振ってしまうという方法もある___& clmn(0:3,Nm,Nr)（例：ith=0 なら v の和を u に，ith=1 なら v**2 の和を

do i=1,np

q に，…）．これはスレッド並列数をある固定値で実行する__l=r(i)/dr

ことがプログラム作成時にわかっている場合にとられる戦__rl=dr*l略である．どのようにスレッド並列した場合に最も速いか__dsi(i)=r(i)-rl

は，np や nx の大きさ，スレッド数 ith，またコンピュー__if(l==li(i)) cycle

タのキャッシュサイズやアーキテクチャにも依存するの__ci(:,:,i)= clmn(:,:,l)! スプラインテーブルを参照

で，各自様々なパターンで試してみてほしい．

__li(i)=l ! li の値を更新

ここまでは配列への書き込みがあるランダムアクセスのループに関しての話であったが，モンテカルロコードでは

end do各マーカー位置における場の量（磁場，電場，温度など）を!

軌道追跡しながら読み込む必要がある．バラバラに動いて

do i=1,np

いるマーカーのインデックス順にマーカー位置の場のデー__Bi(i)=0

タにアクセスするとベクトル計算機の場合はバンク衝突，__ds=dsi(i)キャッシュの小さいスカラ計算機だとキャッシュミスを起__do m=1,Nm

２５７

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013

つけたのは，内側ループが終わるたびに並列スレッドを閉____bmn=ci(0,m,i)+(ci(1,m,i)+(ci(2,m,i)&

じていてはオーバーヘッドが大きくなるので，その抑制の_______& +ci(3,m,i)*ds)*ds)*ds

ためである．

____cs=cos(cm(m)*th(i)-cn(m)*zt(i))

５．
３．
３ MPI̲REDUCE̲SCATTER を使った多次元領域分____Bi(i)=Bi(i)+bmn*cs

割の並列化__end do

粒子コードにおける領域分割については５．
２．
２節でも取

end doり扱ったが，ここでは少し変わった方法として，MPI_RE-

このようにすれば，ciの値はマーカーが１ステップ前にいDUCE_SCATTER を使った多次元領域分割における並列化

た r 方向メッシュ li から他のメッシュに移動した場合だの例を紹介しよう．まず問題の背景として，新古典輸送やけ磁場データ clmn を参照することになる．また，Bi の計

微 視 的 乱 流 シ ミ ュ レ ー シ ョ ン は５次 元 位 相 空 間

算ループにおいても ci はメモリに保存された順番に連続

!
!
"
!
#
!
"$!
"!#におけるプラズマ分布関数を扱う必要があ"

アクセスされる．なお実際にどの程度データ参照数が減るることを指摘しておく．４章の Vlasov シミュレーションでかは解いている問題設定に強く依存する．一般的にイオンは，この位相空間を領域分割し，差分方程式の形で分布関に比べ電子の方が !方向のドリフト幅が小さく，その一方数の時間発展を解いた．一方 PIC 法やモンテカルロ法なで磁力線方向への運動が速いため，より細かい時間刻みでど，粒子コードを使ったプラズマシミュレーションでは，軌道を計算する必要があるので，ここで紹介した方法は電並列化はマーカーに関して行うのがまず基本で，５次元位子軌道計算においてより有益であるといえる．

相空間における場の量に関してはすべての MPI プロセス

余談であるが，計算機によっては cos，sin 計算がベクが同じレプリカをもっていることになる．今，５次元位相トル型命令に置き換えられるため，三角関数だけをまとめ空間を単純化して(x , v)と書く．それぞれの次元を Nx，たループとして回した方が速い場合と，そうでない場合とNv個に分割したセルに関して，そのセルに入るマーカーががある．前者の場 合，cs を 配 列 cs(np)と し て 宣言し，持つ情報から決まる場の量を求めるには，まず各 MPI プロ

cs(i)=cos(cm(m)*th(i)- cn(m)*zt(i))の計算だ け

セス内で集計をとり，次に全 MPI プロセスで総 和を取を独立したループとして先に回した方が速い．また，スる，という過程が必要になる．

レッド並列を i のループに関して掛ければよい並列実効性（サンプルコード

能を期待できるが，粒子数が非常に多くて np*Nm*4 の大きさのデータ ci がキャッシュに乗り切らない場合，上のdx=(xmax-xmin)/Nx ; dv=(vmax-vmin)/Nv

書き方では性能を引き出せない場合がある．その場合，ci

do i=1,npの i と m の次元を入れ替えて宣言しておき，後半のループ__j=int(x(i)/dx)+1を次のように並列実行させる．
（サンプルコード

5.H）

a_tmp(:,:)=0__k=int(v(i)/dv)+1

5.G'）__jk=(k-1)*Nx+j! (j,k)2 次元を 1 次元にする__Bi(:)=0__a_tmp(1,jk)=a_tmp(1,jk)+w(i)!$OMP PARALLEL, private(i,m,ccm,ccn,&__a_tmp(2,jk)=a_tmp(2,jk)+w(i)*x(i)___& ds,bmn,cs) shared(Bi)__a_tmp(3,jk)=a_tmp(3,jk)+w(i)*v(i)……

__do m=1,Nm____ccm=cm(m)__a_tmp(n,jk)=a_tmp(n,jk)+w(i)*……

____ccn=cn(m)end do!$OMP DO

icnt=n*Nx*Nv____do i=1,np

call MPI_ALLREDUCE(a_tmp,a,icnt, &______ds=dsi(i)_________&_________mpi_real8,mpi_sum, &______bmn=ci(0,i,m)+(ci(1,i,m)&_________&_________mpi_comm_world,ierr)& +(ci(2,i,m)+ci(3,i,m)*ds)*ds)*ds

ここで，w(i)は一般にマーカーのウェイト（重み）と呼ば______cs=cos(ccm*th(i)-ccn*zt(i))

れるものであり，上の例では各セルに入ったマーカーの x，______Bi(i)=Bi(i)+bmn*csv，等に対する w の重み付きモーメントを計算している．必____end do

要に応じて５．
３．
２節で説明したスレッド並列における総和!$OMP END DO

演算のチューニングを施すものとする．また，２次元位相__end do

空間のセル位置(j,k)を，変数 jk を用いて１次元化して!$OMP END PARALLELいる（より多次元の場合も同様にして次元を減らせる）．

この例でもスレッド並列は i のループに対して掛けること

さて，このようにして求められた n 種類のモーメント量に注意．内側ループでは，各スレッドは配列 ci のうち長さa(1:n,jk)から何かしらの演算を行って m 種類の新しい

4*np/nomp の部分区間を参照するだけなの で，キ ャ ッ

場の量 b(1:m,jk)を位相空間の各セル(j,k）に対して求シュに乗りやすい．なお!$OMP PARALLEL を外側ループにめる必要があるとする（例えば各セルの中で離散データ
２５８

Lecture Note

5. Coding Techniques of Particle Simulations

w(i)を(x,v)に関して滑らかに補完した場 %"
&!
$"を求
!
#!

H. Naitou and S. Satakecall MPI_REDUCE_SCATTER(a_tmp,a_loc, &めるなど）．この時，a から b を求める計算量がそれなりに______&_____icnta,mpi_real8, &大きいなら，b の計算を領域分割して MPI 並列実行したい______&_____mpi_sum,mpi_comm_world,ierr)と当然思うだろう．サンプルコード 5.H ではすべての MPIプロセスに a の全区間分の集計結果のコピーが配信される

これで a_tmp の全 MPI プロセスの総和のうち，先頭からが，b を計算するには a のうち各プロセスが担当する部分

icnta(0)個の要素が rank=0 に，次の incta(1)個の要区間の情報だけあればよい．このようなケースに使えるの素が rank=1 に，…という形に振り分けられる．位相空間

が MPI_REDUCE_SCATTER である．

についてインデックス jk を用いて１次元化してあるた

まず，MPI プロセス毎に担当する位相空間の区間の始め，どう振り分けされたかイメージが付きやすいだろう

点，終点を決める．ここで総 MPI プロセス数を nproc，各（図１１）．次に各プロセスの担当区間で a_loc から b_loc

MPI プロセスの rank を myrank とする．

を求めてから，すべてのプロセスに b_loc の結果を送信す（サンプルコード

る．それは以下のように行われる．

5.I）

integer, dimension(0:nproc-1) :: idispb,&

call solve_a_to_b(a_loc,l,b_loc,m,mjk)__&_jkst,jked,icnta,icntb

call MPI_ALLGATHERV(b_loc,icntb(myrank),&

Nxv=Nx*Nv______&_____mpi_real8, b,icntb,idispb, &do mp=0,nproc-1______&_____mpi_real8,mpi_comm_world,ierr)__if(nproc<Nxv) then

これで各プロセスの b_loc の icntb（myrank）個の計算____jkst(mp)=mp*Nxv/nproc+1

結果が，mpi_comm_world に含まれる全てのプロセス上____jked(mp)=(mp+1)*Nxv/nprocの配列 b の，先頭から idispb(myrank)+1 番目〜idispb____icnta(mp)=(jked(mp)-jkst(mp)+1)*n(myrank)+ icntb(myrank)に 収 め ら れ る（図１２）．な____icntb(mp)=(jked(mp)-jkst(mp)+1)*m

お，icnt[a,b](myrank)=0 のプロセスも MPI による送__else受信に参加しなければならないことに注意．こうして a____if(mp<Nxv) thenから b の値を領域分割して分散計算し，結果を全プロセス______jkst(mp)=mp+1 ; jked(mp)=mp+1で共有するという目的は達せられた．

______icnta(mp)=n ;icntb(mp)=m

ここで取り上げた２つの集団通信には，事前にデータの____else

分配方法が記述された配列を与えておく必要がある（サン______jkst(mp)=0; jked(mp)=0; &プル中の icnta，icntb，idispb）．３章，４章で取り上_________& icnta(mp)=0; icntb(mp)=0

げた新しい communicator を定義して通信先を指定する____end if

ものに比べるとやや煩雑であるが，icnta などの計算方法____end if

を理解さえすれば様々な応用ができるだろう．なお，筆者__end if!
__idispb(0)=0__do mp=1,nproc-1____idispb(mp)=idispb(mp-1)+icntb(mp-1)__end do
（このサンプルは Nxv<nproc で担当する区間がない rankが生じる場合にも対応している）
また，a のうち担当区間の総和を受け取る作業配列 a_loc と，担当区間における b の値を格納する b_loc を各
rank 毎に次のように宣言しておく．

図１
１ MPI_REDUCE_SCATTER による総和演算データの振り分け．

real, allocatable :: a_loc(:,:),b_loc(:,:)
if(icnta(myrank).ne.0) then__mjk=jked(myrank)- jkst(myrank)+1__allocate(a_loc(n,mjk),b_loc(m,mjk))
end if
準備は以上である．各 MPI プロセスで求めた a_tmp から
a_loc へ振り分けながら総和通信を行うには，コード 5.Hの MPI_ALLREDUCE を次のように置き換える．
図１
２ MPI_ALLGATHERV による各担当区間の計算結果の通信．

２５９

Journal of Plasma and Fusion Research Vol.89, No.4 April 2013

がこの並列化法を取り入れている新古典輸送コード

参考文献

FORTEC-3D［２４］では，５次元位相空間の領域分割数は［１］C.K. Birdsall and C.K. Langdon, Plasma Physics via Computer Simulation (MacGraw-Hill, New York, 1985).
［２］内藤裕志：プラズマ・核融合学会誌 74, 470 (1998).
［３］S. Briguglio et al., in Recent Advances in Parallel Virtual
Machine and Message Passing Interface, Lecture Notes
in Computer Science 2840, 180 (Springer, 2003).
［４］V.K. Decyk and T.V. Singh, Computer Physics Communications 182, 641 (2011).
［５］H. Naitou et al., Phys. Plasmas 2, 4257 (1995).
［６］H. Naitou et al., J. Plasma Fusion Res. 72, 259 (1996).
［７］内藤裕志，徳田伸二：プラズマ・核融合学会誌 72, 737(1996).
［８］H. Naitou et al., J. Plasma Fusion Res. SERIES 8, 1158
(2009).
［９］H. Naitou et al., Progress in Nuclear Science and Technology 2, 657 (2011).
［１０］H. Naitou et al., Plasma Fusion Res. 6, 2401084 (2011).
［１１］H. Naitou et al., Plasma Sci. Technol. 13, 528 (2011).
［１２］http://www.jspf.or.jp/conference/specialist/22prog.
html［１３］http://www.jspf.or.jp/journal/8904koza.html［１４］X.Q. Xu and M.N. Rosenbluth, Phys. Fluids B 3, 627 (1991).
［１５］Z. Lin et al., Phys. Plasmas 2, 2975 (1995).
［１６］A.H. Boozer and G. Kuo-Petravic, Phys. Fluids 24, 851(1981).
［１
７］津田孝夫：モンテカルロ法とシミュレーション（三訂版）（培風館，１９９５）．
［１８］田村義保 他：日本統計学会誌 35, 201 (2006).
［１９］松本 眞：日本統計学会誌 35, 165 (2006).
［２０］D.E. Knuth，有沢 誠 他訳：The Art of Computer Programming (2) Seminumerical algorithms 日本語版（ア
スキー，２
０
０
４）
．
［２
１］http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/mt.
html［２２］http://theo.phys.sci.hiroshima-u.ac.jp/ishikawa/PRNG/
mt̲stream.html［２３］http://random.ism.ac.jp/random/index.php［２４］S. Satake et al., Plasma Fusion Res. 3, S1062 (2008).

１０６程度のオーダーであり，それに対して nproc＝５１２程度
までならこの通信を使った MPI 並列での計算時間の短縮が確認できた（Plasma Simulator および HELIOS 上）．それ以上の並列数にしてもほとんど実行性能は向上しなかったが，多ノード間での MPI 集団通信に付随するオーバーヘッ
ドが増える事と，一方で各プロセス上の b_loc の担当区間
が小さくなり，通信時間に比べ無視できるほど a_loc から
b_loc を求める計算時間が短くなることの両方がその原因となっている．nproc>>1000 のような次世代スパコンの環境でもスケーラビリティを得るためには，５．
２．
２節の「（４）プロセス並列化２」で紹介された領域分割とレプリカによる粒子分割を組み合わせたより高度な並列化が必要となるかもしれない．

５．
４ おわりに４回に渡って連載してきました本講座は今回で終了です．計算機の進化は日進月歩であり，大規模なシミュレー
ションを効率よく実行し，今までに扱えなかった現象をシ
ミュレーションするためには計算機の進化に合わせたプロ
グラミングの知識，技法が必要とされますが，そのような知識は市販されている教科書やシミュレーション結果に関する論文を読んだだけでは得ることができず，そのことが新しくシミュレーション研究を始めようという初心者に対する障壁になっているのではないかという危惧がありました．本講座の一番の目標はそのような方々のために，実際の大規模並列化コードがどのようなことに配慮してどう書かれているのかという実例を紹介しながら，並列化コードの最適化・高速化を図るための基礎知識を学んでいただくことにありました．本講座の内容が今後ますます盛んになっていくと予想されるスーパーコンピュータを活用した
プラズマ・核融合分野のシミュレーション研究を志す方々の一助になれば幸いです．

し内 藤 裕 志山口大学大学院理工学研究科・教授．気がつけば，還暦を過ぎてしまいました．「少年老い易く学なりがたし」．年齢に抵抗して頑張っています．趣味はハングルの勉強と韓国ドラマの鑑賞です．ジャイロ運動論的粒子シミュレー
ションとジャイロ簡約 MHD シミュレーションの両方からプラズマの MHD 現象を解明するのが研究テーマです．

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

ひろ

さ

たけ

しん

すけ

佐 竹 真 介核融合科学研究所核融合理論シミュレー
ション研究系准教授．専門は磁場閉じ込め
プラズマ中の新古典輸送や両極性径電場，
新古典粘性等のシミュレーション研究．趣
味はドライブとサッカー観戦，ギター演奏．最近の海外での日本人サッカー選手の活躍を見て，研究者も世界を相手に挑戦し，渡り合わなければと励まされる気持ちです．

!!!!!!!!!!!!!!!!!!!!!!!!!!

とう!!!!!!!!!!!!!!

ない!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

２６０!!!!!!!!!!!!!!!!!!!!!!!!!!