J. Plasma Fusion Res. Vol.96, No.2 (2020)５
７‐６３

解説

超並列粒子法シミュレーションプログラム自動生成ツールの紹介〜並列プログラミングが初心者にもできる！〜
Introduction to Auto-Generator Tool of Massive Parallel Simulation Code野 村 昴 太 郎１，２），沼 田 龍 介３），八 柳 祐 一４），行 方 大 輔１，２），岩 澤 全 規１，２），牧 野 淳 一 郎１，２）
NOMURA Kentaro１，２）, NUMATA Ryusuke３）, YATSUYANAGI Yuichi４）, NAMEKATA Daisuke１，２）,
IWASAWA Masaki１，２）and MAKINO Junichiro１，２）１）神戸大学惑星科学研究センター，２）理化学研究所 計算科学研究センター，３）

兵庫県立大学大学院シミュレーション学研究科，４）静岡大学教育学部（原稿受付：２
０
１
９年１
０月１
８日）

本稿では，私達が開発・公開している多体シミュレーションプログラム開発フレームワーク「FDPS（Framework for Developing Particle Simulator）」を紹介する．FDPS は，粒子シミュレーションを研究に使っている多くの研究者が，並列化や計算機アーキテクチャ固有のチューニングに多大な時間を費やすことなく，自分の扱いたい問題向けのシミュレーションプログラムを容易に作成できるようになることを目標として開発したフレームワークである．
Keywords:
FDPS, particle simulation, plasma simulation

１．はじめに"複数の演算を並列に行う SIMD 演算ユニットを（コア

我々が FDPS（Framework for Developing Particle Simu-

内 SIMD）

lator）を開発した動機は，現在，大規模な粒子シミュレー"複数もつ演算コアを（スーパースカラー）

ションをする，特に，そのためのプログラムを開発するこ"複数もつプロセッサチップを（マルチコア）とが非常に大変になっている，ということである．FDPS"複数ネットワークで接続した（分散メモリ並列）

はこの開発を容易にすることをめざす．
３０年前であれば，研究に使える粒子シミュレーションの構成を持つ．つまり，４つの違うレベルで並列動作する複プログラムを研究者，あるいは学生が自分で書くことはそ

数のユニットを持つ．しかも，それぞれのレベルで，どうれほど大変ではなかった．普通に書いた上で，ベクトル化

いうふうにすれば上手く並列化できるか，が違う．この中等をすればよかったからである．

で，「スーパースカラー」というのは，CPU が，機械語プ

しかし，現在では，そう簡単ではない．まず，1CPU の計

ログラムの中から，並行して実行できそうなものを同時に算機でも，CPU コアが複数あり，コアの中には SIMD ユ

実行する機能である．したがって，これはある程度ハード

ニットがあるため，これらを有効に使えるかどうかで数十ウェアがやる．しかし，他の３つのレベルのために固有のから数百倍性能が変わる．また，キャッシュメモリも有効最適化が必要で，しかも分散メモリ並列では MPI を使ってに使わないと逆に性能が桁で落ちることもある．さらにプログラム全体を書き換え，例えば空間分割をして粒子をスーパーコンピュータを，と思うと MPI を使った並列化が移動させるといった処理も書く必要がでてくる．

必要である．これらをあわせるとそもそも個人では困難，

並列化の複数のレベルの扱いをさらに困難なものにしてというレベルの大変な作業になる．

いるのが，メモリ階層の存在である．現代の計算機では，大変になっている要因を整理してみると，次のようにま演算器のほうがメモリより速いので，演算器の性能を生かとめられる．

すためには容量は小さいけれど高速なキャッシュメモリを使う．

"並列化の階層の増加

しかし，この，キャッシュメモリは，並列計算機とは決"メモリボトルネックの発生とそれに対応したメモリ階して相性のいいものではない．マルチコアプロセッサで層の増加・複雑化

は，それぞれのコアが独立にキャッシュをもちたいわけだ

以下，それぞれについて簡単に述べる．

が，そうすると，あるコアがメモリのどこかに書いても，現代の典型的な大規模並列計算機は，その同じアドレスのデータを別のコアが自分のキャッシュ
corresponding author’s e-mail: kentaro.nomura@riken.jp

５７!2020 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.96, No.2 February 2020

２．FDPS の基本的な考え方にもっているとか，あるいはまだ誰かのキャッシュの中でしか更新されていなくて主記憶には古いデータがあるもの

大規模並列化が可能で実行性能も高い MPI 並列化粒子

を主記憶から読んだ時にその誰かのキャッシュの中の新し法シミュレーションコードを「だれでも」開発できるよういデータが来ないといけないといった問題が起きる．こうにする，というのが，私たちが FDPS によって実現しよう

いったややこしい状況でもちゃんと整合性がある結果を保としていることである．ここで，「だれでも開発できる」

証するのが「コヒーレントキャッシュ」というもので，そというのは，具体的には，FDPS を使うとのためにコア間で複雑なやりとりをする大規模なハード!MPI を使った並列化は FDPS が勝手にやる

ウェアが必要になる．

!OpenMP を使った並列化も FDPS が勝手にやる

このようなシステムで性能を出すためにはいかにして!長距離相互作用も短距離相互作用も，FDPS を使う人

キャッシュを制御するかが重要になり，キャッシュの特性

は粒子間の相互作用を計算する関数さえ用意すればよを考慮したプログラムを書く必要がでてくる．最近のプロ

くてツリー法とかネイバーリストを使った高速化はセッサでは３レベルから４レベルのキャッシュを持つよう

FDPS が勝手にやるになっている．しかし，密行列乗算ならともかく，それ以外の計算アルゴリズムではこのような複数レベルのキャッということを意味する．

シュを有効に使えるとは限らないし，できたとしてもプロ我々が非常に単純な粒子系のプログラムを書く時にはグラムは極めて複雑なものにならざるを得ない．

１．MPI や OpenMP を使った並列化はしない

このため，単純なアルゴリズムでも，最新の高性能プロ

２．長距離力も短距離力もまずは全粒子からの力を単純にセッサで性能を出すのは容易なことではなくなっている．

計算する

実際，今この文章を読んでおられる読者の中でも，俺は３．相互作用計算ループも単純に書く．SIMD 化とか特にMPI で並列化してキャッシュも有効利用し SIMD 演算器も意識しない使ってプラズマシミュレーションプログラムを書いた，あるいは書ける，あるいは書く気がある，という人はあまりというふうにしたい．そういうふうに単純に作ったプログ

多くないのでは，と思う．

ラムを，あまり手をかけないで並列化・高速化できること

では，どうすればいいのか？というのがここでの問題でが目標である．

ある．多くの場合にとられているアプローチは，大規模な

そんなうまい話が本当にあるか，というわけだが，粒子

ソフトウェアを開発チームを作ることでなんとか開発しよ系のシミュレーションプログラムはどういう構造をしていう，というもので，実際，様々な応用分野で，粒子法の並るか，を考えてみると，基本的な構造は以下のようになる．

列化されたプログラムが公開され，利用可能になっている．

１．粒子の初期分布を作る（ファイルから読む・内部で生しかし，そういったプログラムは，あらかじめ開発グ成する）

ループが実装した機能を，開発グループがターゲットにし２．粒子間相互作用を計算して，各粒子の加速度を求める

たマシン・OS で使うことしかできないのが普通である．

３．速度をアップデートする（時間刻みの中央まで）

もちろん，オープンソースで公開されているものでは，原４．位置をアップデートする

理的にはソースコードを修正していろいろな機能を実装で５．ステップ２に戻る

きるわけであるが，巨大で複雑なプログラムで，さらに特定のアーキテクチャ向けの最適化されたコードを修正して

ここでは時間積分にはリープフロッグを使うとした．粒動くようにするのは容易なことではないのは，やってみよ

子法ではリープフロッグや，それと等価な方法を使うことうと思ったことがある人は良くご存じのことかと思う．

が多いと思う．ルンゲクッタ等を使うなら相互作用の計算なので，自分でプログラムを開発できるようにしたいわを中間結果を使って行う必要があるが，いずれにしても相けだが，それにはどうすればいいのか，というのが我々の互作用を計算する部分とそれ以外をする部分があるのは同問題意識である．我々が提案する方法は，特に MPI にかかじである．

わるような複雑な並列化とそのために必要なプログラム

これを，MPI で並列化された，領域分割し，領域毎に自と，実際に扱う系の記述や時間積分の方法の記述とを明確

分の担当の粒子を持つプログラムにするとすれば，に分離することである．

１．粒子の初期分布を作る（ファイルから読む・内部で生明確に分離，と書くのは簡単だが，実際にどのように実成する）現するか，十分に色々なことを表現するにはどうするか，

２．領域分割のしかたを決める

計算速度がでるようにするにはどうするのか，と，いろい

３．粒子を担当するプロセスが持つように転送するろな問題がある．本解説では，まず基本的な考え方を２節

４．粒子間相互作用を計算して，各粒子の加速度を求める

で，またいくつかの FDPS をつかったサンプルを３節で紹５．速度をアップデートする（時間刻みの中央まで）介する．４節はまとめである．

６．位置をアップデートする

１，２，４節は牧野他 FDPS チームが，３節はユーザ代表

７．ステップ２に戻るとして野村，沼田，八柳が担当した．

となる．主な違いは，ステップ２，３が入ること，ステップ
５８

Commentary

Introduction to Auto-Generator Tool of Massive Parallel Simulation CodeK. Nomura et al.

４では，各プロセスは自分の担当の粒子の加速度を計算す

ンピュータ「京」のような計算機でも高い性能を出すことるのに必要な情報を他のプロセスからもらってこなければができるし，その同じプログラムをノートパソコンで走らいけないことである．

せることもできる．

ステップ１は，大規模並列では入力ファイルを並列に読もちろん，実際にはコンパイルされたプログラムではむ必要があったりして若干面倒だが，難しいものではな

FDPS 経由で MPI が呼ばれている．ということは，ユーザ

い．ステップ５，６も，単純に全粒子に対して，計算されたが書いた同一のプログラムがコンパイルしなおすだけで複加速度を使って時間積分の公式を適用するだけで，特に難数の MPI プロセスで並列に実行されたり，１コアで実行さしいことはない．積分公式がなにか違うものでも，全粒子れたりする．MPI 実行の場合には，FDPS が勝手に領域分に適用する，ということには変わりない．

割をして粒子をプロセス間で交換する．

FDPS がすることは，上のステップ２，
３，４のそれぞれ

このため，ユーザプログラムから見ると知らないうちにについて，そのためのライブラリ関数を提供することであ勝手に粒子が入れ替わって，粒子の数も FDPS が変更する．ステップ２，３のためには，ライブラリ関数は粒子を表る．もちろん，MPI プロセス全部で見ると，どこかに粒子すデータがどういうものかを知る必要があるし，また，があるが，調べなければどこにあるかはわからない．これ

MPIプロセス間で計算時間の差がでないようにするためには領域分割する粒子系プログラムでは必ず起こることで，なんらかの計算負荷に関する情報をもらう必要もある．

調べられるようにしておくためには，粒子データ構造体のライブラリが空間分割を決めるためには，例えば粒子の中に粒子のインデックスをいれておくのが普通のやり方で座標の配列をもらうことができればいいが，粒子を転送すある．

るためにはその粒子がどういうデータから構成されるか，
メモリにどう配置されているかをライブラリが知る必要が

３．実例

ある．

３．
１ 単純な分子動力学計算

FDPS では，C++言語のクラスを使って，オブジェクト

本節では，希ガスの単原子分子を模した Lennard-Jones

として粒子を表現してもらうことで，ユーザプログラム側（LJ）粒子のシミュレーションコード例について説明する．

で定義した粒子データの操作をライブラリ側で行うことを

LJ 粒子 "と #の２体間相互作用は以下のポテンシャルで表可能にしている．具体的には，粒子が「オブジェクト」でされる．

表現されることで，粒子の代入（コピー）が可能になり，まけ取って相互作用を計算する関数をユーザ側で定義するこ

"
"#
%
"
$$
(）
##!"" !! "$ （if $
"
#!$
$
$
$
"
#
"
#
（１）
$
!'&#
"#
"
#$
$
!
（else）%とで，ライブラリ側は粒子データの中身を知らないままで#と "はそれぞれポテンシャルの深さと粒子の大きさを特相互作用を計算する関数を呼び出せる．

徴づけるパラメータで，$(はカットオフ距離である．本解た粒子の位置や電荷・質量を返す関数を提供してもらうことで，長距離力のためのツリー構造を作ることもできる．
さらに，相互作用を計算する関数は，粒子データ自体を受説では簡単のために #""""としている．

FDPS はこのように C++言語の機能を使って実装されているが，現在は Fortran 言語（２００３以降），C 言語にも対応

紙面の都合上，最低限必要な部分のみを解説していく．

している．近代的な Fortran では，C++のクラスに相当す

ある程度，C もしくは C++の知識があることを前提としてる「構造型」があり，また，C 言語で作った構造体や関数をいる．Fortran インターフェースを用いる場合でも，同様

Fortran 側からアクセスする方法も，言語仕様として定義の情報を書く必要があるので参考にはなると思われる．

された．これらの機能を利用して，C 言語や Fortranで書か

コード全体については https://github.com/FDPS/JSPF̲

れたユーザプログラムからも FDPS を利用可能にしていSample を参照してほしい．
まず，コードの本体である main.cpp である．ここでは，る．
FDPS を使って粒子系シミュレーションプログラムを書

１．FDPS のヘッダーファイルのインクルード

くには

２．粒子データクラス（クラス名 FP，EP，Force）の定義!粒子を「構造体」で表現する

３．粒子間相互作用関数（関数名 Kernel）の定義!粒子間相互作用をその粒子構造を引数にとる関数で書く

４．main 関数の定義が行われている．FDPS のヘッダーファイルのインクルーという FDPS 側の要請に従う必要があるが，それによっドはファイルの先頭で行われている．

て，ユーザプログラムは MPI による並列化をほとんど意識することなく，MPI 並列でないプログラムと同程度，ない

1 #include <particle_simulator.hpp>

しはより少ない手間で書くことができる．また，コンパイ

2 using namespace PS;

ル時のフラグを変えるだけで，MPI を使わないこともでき

3 using namespace std;

るし，OpenMP での並列化も，また MPI と OpenMP を組み合わせたハイブリッド並列もできる．したがって，ハイ

また，ここでは記述を短くするために名前空間 PS（Parti-

ブリッド並列でないと大規模な実行ができないスーパーコ

cleSimulator の略）と std（C++標準ライブラリ）で定義さ５９

Journal of Plasma and Fusion Research Vol.96, No.2 February 2020

れている変数や関数，クラスを使うことを宣言している

7（PS::が省略されている）．

8

粒子データクラス（構造体）の定義では，粒子が持たな

9くてはならない情報と FDPS が必要とするデータを操作す

10

F64 pi = force[i].p;for(S32 j=0; j<nj; j++){F64vec rij = ri − epj[j].r;
const F64 r2 = rij * rij;if(r2==0.0 || r2>rc2) continue;
const F64 r2i = 1.0/r2;
const F64 r6i = r2i * r2i * r2i;fi += r6i*(48.0*r6i−24.0)*r2i * rij;pi += 4.0*r6i*(r6i−1.0);}
force[i].f = fi;force[i].p = 0.5*pi;

11る関数を定義する．FP クラスについて見てみよう．FP

12は FullParticle クラスの略で，時間発展や相互作用を計算

13

するにあたって必要な情報全てを持つ粒子データクラスで

14

ある．

15
16

1 struct FP{

17

F64vec r,v,f; // 座標，速度，力F64p;// ポテンシャル
void copyFromForce(const Force& _f){f = _f.f; p = _f.p;}
F64vec getPos() const {return r;};void setPos(const F64vec& _r){r = _r;}
void kick(const F64 dt){v += f*dt;}
void drift(const F64 dt){r += v*dt;}

2
3
4
5
6
7
8
9

18

}

19
20 }

main関数の定義では，シミュレーションを行うのに必要なインスタンスを作成し，実際にシミュレーションを行う．main 関 数 で は，ま ず 最 初 に PS::Initialize を 呼 び 出
す．この関数は FDPS に必要な初期化を行う（主に MPI

10 };プロセスの初期化など）．
粒子データクラスが持つメンバ変数として座標，速度，力，1 Initialize(argc,argv);

ポテンシャルエネルギーと，FDPS がこのデータを操作する の に 必 要 な メ ン バ 関 数（copyFromForce，getPos，

次に ParticleSystem や DomainInfo，TreeForForce を定

setPos）が定義されている．また，kick と drift という粒子

義した粒子データクラスを使用してインスタンスにして，の時間発展を行うメンバ関数が定義されている．このよう各インスタンスにおいて Initialize 関数を呼び出す．本コーに，FDPS が必要としていないメンバ関数もユーザの必要

ドでは，初期条件として粒子を格子上に配置しているが，に応じて追加することができる．粒子データクラスの定義

ここでは解説を省略する．

では，更に Force と EP（EssentialParticle の略）クラスを
1 ParticleSystem<FP> ps;

定義している．それぞれ，相互作用計算を行う際の結果

2 ps.initialize();（力とポテンシャル）と相互作用計算に必要（Essential）な
データ（ここでは座標のみ）を持ったクラスである．それ

1 DomainInfo di;

ぞれについて，FDPS が必要とする関数が少しずつ異なる．

2 di.initialize(0.3);

このように相互作用計算のみに使われる粒子データクラス

3 di.decomposeDomainAll(ps);を定義して相互作用計算のときに利用することにより，メ4 ps.exchangeParticle(di);モリ使用量の削減や最適化を行いやすくしている．EP ク

5 TreeForForceShort<Force,EP,EP>::Scatter t;

ラスは，!粒子と "粒子（相互作用を受ける粒子とおよぼす

6 t.initialize(3*n*n*n,0.0,64,256);
7 t.calcForceAllAndWriteBack(Kernel,ps,di);粒子）それぞれについて定義できるが，本コードでは両方に同じ EP クラスを用いている．最も簡単にコードを書く

以下のループでは実際にシミュレーションを進めている．

場合，Force クラス，EP クラスを定義せずに FP クラスを

時間発展（kick，
drift）を行い，空間の分割（DomainInfo::
decomposeDomainAll），粒 子 の 分 配（ParticleSystem::

全ての用途で使い回すこともできる．

exchangeParticle），力 の 計 算（TreeForForce::calcForce

相互作用関数の定義では，上記のポテンシャルから力の計算を行っている．この関数は，!と "粒子のリストとそのAllAndWriteBack）を行っている．

長さ（epi，epj と ni，nj）を受 け 取 っ て，結 果 の リ スト
1 S64 nl = ps.getNumberOfParticleLocal();（force）に計算結果を書き込んでいる．FDPS は FP のリス

2 for(int s=0;s<10000;s++){

トから !と "粒子のリストを作り，相互作用関数を呼び出

3

し，結果を FP に書き戻す．この際，分散メモリ上の並列計

4算については，復数のプロセス間での通信などが必要に応5

じて行われる．

6
7

1 void Kernel(const EP *epi,const S32 ni,
2
3
4
5
6

8

const EP *epj,const S32 nj,
Force *force){
const F64 rc2 = RCUT*RCUT;
for(S32 i=0; i<ni; i++){F64vec ri = epi[i].r, fi = force[i].f;

9
10
11
12
13

６０

for(int i=0;i<nl;i++){ps[i].kick(dth);
ps[i].drift(dt);}ps.adjustPositionIntoRootDomain(di);di.decomposeDomainAll(ps);ps.exchangeParticle(di);nl = ps.getNumberOfParticleLocal();t.calcForceAllAndWriteBack(Kernel,ps,di);for(int i=0;i<nl;i++) ps[i].kick(dth);
const F64 pot = AccumulatePotential(ps);Commentary

Introduction to Auto-Generator Tool of Massive Parallel Simulation Code

下付き添え字 $""""は粒子番号，!
*"$，"$は粒子の位置，

const F64 kin = CalcKineticEnergy(ps);if(Comm::getRank()==0)
cout << scientific << pot << ” ” << kin<< ” ” << pot+kin << endl;

14
15
16
17

K. Nomura et al.

速度，($，&$は粒子の電荷，質量である．ハット記号は無次元量をあらわし，各変数は以下のように無次元化されて
いる．

18 }

3
!
$"

最後に PS::Finalize を行い，終了処理を行う．
1 Finalize();

!
$
，
%)!

3
"
$"

&$
，
&!

(
($" $
(!

3$"
&

次 に Makefile を 見 て み よ う．FDPS̲PATH は githubからダウンロードしてきたファイルのパスである．FDPS

"$
，
*1!##
)

3")
)
'/!

下付き添え字０は，参照粒子をあらわす．

を 用 い た コ ー ド を コ ン パ イ ル す る 際 に は＄
（FDPS̲

#
'!(!
，
'/!#
"(
!&!

PATH）/src/particle̲simulator.hpp を include する必要があるためヘッダーファイルの検索パスを追加している．

(
!#
!
%)!#
#
" '!(!

（３）

MPI を用いて並列化する場合は，MPI 用コンパイラを用いは，て，PARTICLE̲SIMULATOR̲MPI̲PARALLEL マ ク ロ

長であり，これらを用いると熱速度は *1!")#%)!'/!で与を 有 効 に す る 必 要 が あ る．OpenMP を 利 用 す る 際 は，$えられる．$!#%
は参照プラズマのプラズマパラ&'!%)!

OpenMP 用 の オ プ シ ョ ン を 有 効 化 し て PARTICLE̲

メータであるが，Coulomb 相互作用の強さを規定するパラ

SIMULATOR̲THREAD̲PARALLEL マクロを有効化す

メータである．なお，密度と長さの単位は独立ではないたることでスレッド並列化を行うことができる．このとき，め系全体の（平均）密度 '*2-は重み + を用いて以下のよう
3$'．ここに無次元化される：'*2-"+"#
!$"&
+#
%$ '
"#
!
&

時間発展などのユーザー記述部分はスレッド並列化されな

それぞれ参照プラズマのプラズマ周波数，Debye)!

$
で，'!"+#となるので，$!"%&+ に他ならない．
%)!

いので注意されたい．また，両方を用いるハイブリッド並列化も可能である．

コードは，FDPS に含まれるサンプルコード nbody の重力相互作用を Coulomb 相互作用に置き換えればほぼ完成

1 INC+= ‐I$(FDPS_PATH)/src

である．周期境界条件を実現するためには，Coulomb 力を

2 CXX=g++粒子―粒子と粒子―メッシュの相互作用の和で表現する

3 CXXFLAGS= ‐O2 $(INC)
4 #FL+= ‐DPARTICLE_SIMULATOR_THREAD_PARALLEL

Particle-Mesh 法（FDPS の拡張機能として提供されてい5 #FL+= ‐DPARTICLE_SIMULATOR_MPI_PARALLELる）を用いるが，この部分もサンプルコード p3m を参考に

6 SRC=main.cpp

作成した．
サイズ !"%)!の領域の中心に固定された電荷 (.のイオ78 all:

argon.out
9 argon.out:$ (SRC) Makefile
10$(CXX) $(FL) $(SRC) -o argon.out

ンを置き，","'""(#個の電子によるポテンシャルの遮蔽のシミュレーションを行う．$!""!とした．イオンの電荷と質量を (."",(,，&.""!"!&,とする（立体の添え字

LJ 粒子のサンプルコードの解説は以上である．シリアル，,はそれぞれイオン，電子をあらわす）．電子の初期位置.

なコードをコンパイル時にマクロを変更することで簡単にはランダムに配置し，初速はゼロとする．系の全エネル

並列化が行えることがわかってもらえただろうか．

ギーが保存するミクロカノニカルアンサンブルを考えてい３．
２ プラズマの例１：Coulomb 相互作用するプラズマのるため，電子はある程度の時間が経過した後で温度一定のシミュレーション

熱平衡状態に達する（図１）．平衡状態における系全体の分子動力学モデルを用いたプラズマの研究として，強結

平均電子温度は #,0*1##!$%%!となった．このとき電子のプ合プラズマ（下記で定義するプラズマパラメータ $が１よ

ラ ズ マ パ ラ メ ー タ は $,$"!$，Debye 長 は %),#
!$!
!
#$

り小さい）の物性を取り扱ったものがある［１］．ここでは，となる．図２に -方向 の 静 電 ポ テ ン シ ャ ル 分 布 を 示 す

結合の強さを限定せず，主に弱結合の核融合・宇宙プラズ

（,#
!"!
!
!"!
!
&とした）．電子によるポテンシャルの&，.#マのダイナミクスを議論する目的で粒子モデルに基づくシ遮蔽がよく再現できることが確認された（境界付近でのずミュレーションを行う．同様のアプローチの研究が［２，３］

れは，サンプル粒子数が十分でないためであると思われにある．また，Escande らは粒子モデルを用いて理論的にる）．ここでは，系の全エネルギーが一定のアンサンブル

Debye 遮 蔽 や Landau 減 衰 を 議 論 し て い る［４］．本 節 でを考えているが，温度一定のカノニカルアンサンブルを考は，周期境界領域内（領域は一辺 !の立方体とする）でえるほうが実際のプラズマとの比較がしやすい．温度を一Coulomb 相互作用するプラズマの集団現象として Debye

定に保つ熱浴の実装を検討している．また，ここでは粒子

遮蔽のシミュレーション例を示す．Coulomb 相互作用する

数が１０４程度に留まり，FDPS の威力は十分発揮されている" 個の粒子の運動方程式は以下のように書ける．

とは言えないが，今後大規模計算を実施する際にその効果

3
+!
$
3
""
$，
3
+)

"
3
3
3
3
!
+"
" ! (
$(
%
$
$
"
$.
3 $!%"""%%
3$ (
3
3
"$ &
!
+)
$!!
%(

を体感できることを期待している．
（２）

６１

Journal of Plasma and Fusion Research Vol.96, No.2 February 2020

ラズマ実験では電場と磁場により流体内部が制御可能な流体実験を行うことができ，渦結晶に代表される自己組織化現象を観測してきた［８，９］．非粘性渦度方程式は，点渦法で数値的に解くことが可能なため，筆者（八柳）は，純電子
プラズマ実験で見られる自己組織化現象とリンクした理論／シミュレーション研究を展開してきた．
２次元点渦系では，渦度方程式
#
'
&'(
!
&
)
#
$(
#
&'(
!
&
!
&
"'$
(
)
)
)
%!
'
&

（４）

の「形式的」解を Dirac のデルタ関数 $(#)で表現された点の集合体とする．点の総数は ! で表す．
!

図１

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

電子温度の時間変化．温度は x，y，z 各方向でほぼ等しい
ため，３方向の平均温度を取った．

$%"

（５）

ここで &'(#!
&)は渦度，#$，#$は $番目の点渦の強さを表
す循環，２次元位置ベクトルである．点渦の時間発展は，離散化された Biot-Savart 積分となる．
!
&
#
#%
%#
" ! #(
$!#
%)
$
%!
$
#
%%%"(&
%&
#
#
*
%$
)
$!#
%*
!
&
#!# )
(
#%
"
" ! #$ $ % #
（６）
%%%"
#
#
*
$!#
%*

右 辺 第２項 は 半 径 "の 円 形 境 界 を 表 現 す る
図２##$%"###*$"$*に置かれた鏡像渦からの影響を表す．
数値シミュレーションでは，適当な初期分布で #$を初期

静電ポテンシャルの空間分布．

化したあと，相互作用（６）を計算する．今回は，わかりや３．
３ プラズマの例２：純電子プラズマ実験に対応した２すい例として，Kelvin-Helmholtz（Diocotron）不安定性の次元点渦シミュレーション計算を行った．結果を図３に示す．初期設定から求めた線

本節では，FDPS の開発チームより，２次元系の計算が形不安定モードは，不安定な順に２と３である［１０］．その

正しく行えるか，チェックを行ってほしい旨の依頼を受け中から，数値揺らぎにより，多くのエネルギーが蓄積されたので，点渦系シミュレーションを FDPS を用いて実装した モ ー ド３が 立 ち 上 が っ て い る と 考 え ら れ る．ま た，検証を行った．なお，点渦系は，モデルに内在する運動論#%"#!
!
"$!のプロットから，境界をまたいで越境してし的効果により［５］，現象を再現できる下限の粒子数で計算まう粒子がいないことから，円形境界も鏡像渦によって十を行った結果が，一番理にかなったものになるケースが多分な精度で計算されていることがわかる．

いため，FDPS が得意とする多粒子，多並列でのシミュ

紙面の残りを使って，今回，著者がつまずいた点，特にレーションを行う意味が薄い．よって，本解説は，FDPS

MPIに慣れていないことに由来する困難について述べてみで２次元シミュレーションも作成できることを示すことをたい．FDPS では領域分割技法である Tree 法や，疎結合並第一目的とし，ベンチマーク結果は掲載しない．

列計算のためのAPIであるMPIに関する知識はまったく必

２次元点渦シミュレーションでターゲットとする物理系

要ない．しかし，FDPS 側でどのような内部処理を行っては，純電子プラズマ実験である．一般的なプラズマと異ないるのか，その概要を理解していないと，どうにも解決でり，電気的中性条件が破れたプラズマを非中性プラズマときないバグを埋め込む結果となる＊２．粒子系におけるクー

よび［６］，特に電子のみから構成される場合，純電子プラ

ロン力の計算は，力を計算する場所の粒子（$粒子）それぞ

ズマとよぶ．軸方向に強磁場を掛けた真空容器内におけるれについて，力の源泉となる粒子（%粒子）からの相互作用

磁場に垂直な断面内の電子の運動は，Larmor 半径＝０のを計算する．このとき，非 MPI 環境下であれば，粒子に関極限（案内中心近似）で，"#!ドリフト運動となり，その連付けられる添字の値が変化することはないので，同一粒時間発展方程式は２次元非圧縮流体方程式である Euler 方子間で発散する相互作用計算を回避する目的で，ループ内

程式（正確には，回転微分をとった渦度方程式）で記述さに お い て $&%%と い う 条 件 が 使 用 可 能 で あ る．し か し，れる［７］．このとき，流れ関数が静電ポテンシャル，渦度

FDPS の場合，各粒子は時々刻々変化する位置に応じて領が電子の数密度に比例する＊１．この関係により，純電子プ

域ごとに配分しなおされ，各 MPI ノードへ渡される $粒子＊１ この関係の導出過程はベクトル解析のよい練習問題なので，学生の皆さんはぜひチャレンジしてみてください．
＊２ 開発チームのみなさま，問題解決に協力していただき，大変感謝しております．

６２

Commentary

Introduction to Auto-Generator Tool of Massive Parallel Simulation CodeK. Nomura et al.

る．指摘されてみれば当たり前の結論であるが，MPI でのプログラミング経験がないと，障壁になりやすいと感じた
ので，敢えてここに tips として掲載した．

４．おわりに本稿では，粒子系大規模並列化のためのフレームワーク
FDPS を紹介した．現代の大規模計算では MPI 等での並列化が必須になっているが，FDPS を使うことでその部分は自分で開発しなくても，高性能な並列プログラムを実現できる．本稿が読者の皆様の研究のためのプログラム開発の一助になれば幸いである．

参考文献

図３

［１］S. Hamaguchi, J. Plasma Fusion Res. 78, 313 (2002).
［２］A. Panarese et al., J. Plasma Phys. 84, 905840308 (2018).
［３］C.-P. Wang and Y. Nishimura, IEEE Trans. Plasma Sci.
47,1196 (2018).
［４］D.F. Escande et al., Rev. Mod. Plasma Phys. 2, 9 (2018).
［５］Y. Yatsuyanagi and T. Hatori, Fluid Dyn. Res. 47, 065506(2015).
［６］毛利明博：プラズマ・核融合学会誌 77, 213 (2001).
［７］際本泰士：プラズマ・核融合学会誌 77, 329 (2001).
［８］D.Z. Jin and D.H.E. Dubin, Phys. Rev. Lett. 84, 1443 (2000).
［９］A. Sanpei, et al., Phys. Rev. E 68, 016404 (2003).
［１０］R.C. Davidson, Physics of nonneutral plasmas (AddisonWesley, Reading, 1990) §6.

円形境界半径を１として，ドーナツの外側半径 ＝ 0.8，内側
半径 ＝ 0.46，総粒子数 ＝ 4177，ドーナツの内部も点渦で満たされた円盤の剛体回転タイムスケール ＝ 10.8．

の個数をほぼ同一にするような仕組みをもっている．このため，粒子の添字が時々刻々変化するため，自己相互作用計算を排除する目的で !"!"という条件は使用できず，!
粒子と "粒子の距離が0.0超か否かにより判定することにな

６３

Journal of Plasma and Fusion Research Vol.96, No.2 February 2020

たろう

野 村 昴 太 郎
神戸大学大学院理学研究科惑星科学研究センター．分子動力学（MD）シミュレーショ
ンを用いたナノチューブ閉じ込め系の水に関する相挙動研究を行ってきた．あわせて，Xeon Phi や GPU，PEZY-SC など様々なプロセッサを用いた MD シミュレーションコード開発を行った．現在は，N体，MD，SPH 向けに各種プロセッサに最適化された相互作
用計算カーネル開発および，プレファードネットワークスと開発している深層学習向けプロセッサのシミュレーション向け応用を進めている．

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

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

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

けん

やつ やなぎ ゆういち

八 柳 祐 一静岡大学教育学部 准教授．長距離相関に支配された N 体系の自己組織化に関する理論／シミュレーション研究．最近，講談社
青い鳥文庫版のシャーロック・ホームズ・シリーズを再読中．既読の新潮文庫版と比べて，１
９世紀末のロンドンの時代背景描写が小中学生向けに豊富なので，再発見が多数（冷や汗）．アイリーンの印象が変わりました．

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

むら!!!!!!!!!!!!!!

の!!!!!!!!!!!!!!!!!!

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

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

沼 田 龍 介兵庫県立大学大学院シミュレーション学研究科 准教授．２
００
４年東京大学大学院新領域創成科学研究科修了・博士（科学）．磁気リコネクションや乱流などの非線形複雑現象のシミュレーション研究を行っています．粒子モデルでプラズマの乱流は再現できるでしょうか．昨年，KOBE 六甲
全山縦走大会に参加，公称 55 km を約１
５時間かけて完走しま
した．人生で最長距離行でしたがなんとかなるものです．

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

!!!!!

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

りゅう すけ

なめ

かた

だい

すけ

行 方 大 輔!!!!!!!!!!!!!!!!!!!!!!!!!!

た

!!!!!

ぬま

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

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

!!!!!!!!!!!!!!!!!!!!!!!!!!
まさき

神戸大学理学研究科 特命准教授 ２００
９年東京大学総合文化研究科広域科学専攻修了
博士（学術）．専門は数値天文学．またそのためのソフトウェアやアルゴリズムの研究開発も行っている．

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

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

さわ

まき の じゅんいちろう

牧野 淳 一郎!!!!!!!!!!!!!!

いわ岩 澤 全 規!!!!!!!!!!!!!!!!!!!!!!!!!!

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

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

神戸大学理学研究科惑星学専攻／理研計算科学研究センター．球状星団の多体シミュ
レーションから研究をスタートし，多体問題向け専用計算機 GRAPE の開発から，「京」，富岳の開発にもかかわってきました．今はプレファー
ドネットワークスと進めている深層学習向けプロセッサの様々なシミュレーションへの応用を進めています．

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

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

６４