FDPSを使ってみる その3

# はじめに

粒子系の並列シミュレーションフレームワーク[FDPS](https://github.com/FDPS/FDPS)を使ってみる。
[その2](http://qiita.com/kaityo256/items/8da0ac002a6d8d0fc62d)では粒子系の定義、領域分割までやった。まずは相互作用の無い系で時間発展させてみる(基本)。

コードはここにおいてある。
https://github.com/kaityo256/fdps_sample

# 変数名

* `PS::ParticleSystem<FPLJ> lj_system` 粒子系
* `PS::DomainInfo dinfo` 領域情報

なお、前回はテンプレートがどのような実体を持つか明示するために

```cpp
dump< PS::ParticleSystem<FPLJ> >(rank, lj_system);
```

などとしていたが、今回は型を省略して

```cpp
dump(rank, lj_system);
```

となどと呼ぶことにする。

# 時間発展

時間発展を記述するには、運動量と座標を更新すれば良い。相互作用がなければ、単に座標を運動量で更新するだけである。記述すべき関数`drift`はこうなる。

```cpp
template<class Tpsys>
void
drift(Tpsys & system, const PS::F64 dt) {
  const PS::S32 n = system.getNumberOfParticleLocal();
  for(PS::S32 i = 0; i < n; i++) {
    system[i].q  += system[i].p * dt;
  }
}
```

`PS::ParticleSystem::getNumberOfParticleLocal()`で、「自分の領域が持っている粒子数」が取得できるので、それでループを回せば良い。引数として時間刻み`dt`を受け取っている。

これを単に呼べばよいが、問題は周期境界条件の補正である。全体の領域(root domainと呼ばれている)をはみ出した粒子がいる状態でドメイン間の粒子交換を行おうとするとFDPSに怒られる。

```shell-session
PS_ERROR: A particle is out of root domain 
function: exchangeParticle, line: 773, file: 
/path_to/FDPS/src/particle_system.hpp
```

FDPSには境界条件補正のための関数、`PS::ParticleSystem::adjustPositionIntoRootDomain`が用意されているので、それを呼び出してやれば良い。この関数を実行するためには、フレームワーク側に座標をいじらせる方法を提供する必要がある(粒子クラスはユーザが自由に記述できるため、フレームワーク側はどれが座標かわからないから)。これは、[その2](http://qiita.com/kaityo256/items/8da0ac002a6d8d0fc62d)で作った`FPLJ`クラスに`setPos`というsetterを設定すれば良い。

```cpp
class FPLJ {
public:
  PS::F64vec p;
  PS::F64vec q;
  PS::F64vec getPos() const {
    return q;
  }
  void setPos(PS::F64vec nq) {
    q = nq;
  }
};
```

この状態で、

```
lj_system.adjustPositionIntoRootDomain(dinfo);
```

を呼べば、情報が更新される。

以上から、時間発展の1ステップは以下のようにかける。

```cpp
    drift(lj_system, dt);
    lj_system.adjustPositionIntoRootDomain(dinfo);
    dinfo.decomposeDomainAll(lj_system);
    lj_system.exchangeParticle(dinfo);
```

領域分割は毎ステップ呼ぶ必要は無いが、まずはサンプルということで気にしないことにする。

ついでに、各プロセスごとの粒子座標をダンプする関数も修正して、あとでアニメを作れるようにしておこう。

```cpp
template<class System>
void
dump(const int rank, System &sys) {
  static int step = 0;
  char filename[256];
  sprintf(filename,"dump%02d_%03d.dat",rank,step);
  step++;
  std::ofstream ofs(filename);
  const PS::S32 n = sys.getNumberOfParticleLocal();
  for (PS::S32 i = 0; i < n ; i++) {
    ofs << sys[i].q.x << " ";
    ofs << sys[i].q.y << " ";
    ofs << sys[i].q.z << std::endl;
  }
}
```

`std::stringstream`使うのが面倒になったので`sprintf`で。

これを適当なインターバルで呼び出すと、最終的に時間発展ループは以下のようになる。

```cpp
  const PS::F64 dt = 0.01;
  const int LOOP = 1000;
  const int INTERVAL = 10;
  dinfo.decomposeDomainAll(lj_system);
  lj_system.exchangeParticle(dinfo);
  for(int i=0;i<LOOP;i++){
    if(i%INTERVAL == 0){
      dump(rank, lj_system);
      printf("(%d/%d) n= %d\n", rank, procs, lj_system.getNumberOfParticleLocal());
    }
    drift(lj_system, dt);
    lj_system.adjustPositionIntoRootDomain(dinfo);
    dinfo.decomposeDomainAll(lj_system);
    lj_system.exchangeParticle(dinfo);
  }
```

# 結果

実行すると、各ステップで各プロセスの持つ粒子数を出力しながら結果をダンプする。システムサイズは前回より小さくして`L=10`とし、1000粒子系とした。

```shell-session
$ make
$ mpirun -np 4 ./a.out
******** FDPS has successfully begun. ********
(snip)
(0/4) n= 200
(2/4) n= 200
(3/4) n= 300
(1/4) n= 300
(snip)
(0/4) n= 254
(3/4) n= 254
******** FDPS has successfully finished. ********
```

するとdump??_???.datが大量にできるので、それをこんなgnuplotスクリプトで処理する。

```gnuplot:dump.plt
set pointsize 2
set size square
do for [i=0:99] {
  set term png
  set out sprintf("plot%03d.png",i)
  plot sprintf("dump00_%03d.dat",i)\
  ,sprintf("dump01_%03d.dat",i)\
  ,sprintf("dump02_%03d.dat",i)\
  ,sprintf("dump03_%03d.dat",i)
}
```

これで連番pngを作ってconvertすればアニメーションGIFを作ることができる。

```shell-session
$ gnuplot dump.plt
$ convert plot*.png plot.gif
```

こうしてできたのが以下のGIF。異なるプロセスに所属する粒子が異なる色で表現されている。

![plot.gif](https://qiita-image-store.s3.amazonaws.com/0/79744/3fd2595e-f2ed-60be-adb1-3d8d9ce2caec.gif)

ただしく計算できているっぽい。

# まとめ

FDPSで相互作用無しの粒子系(要するに理想気体)を計算させてみた。シリアル処理なら簡単だが、理想気体でもMPIを使って並列化しようとすると領域間の粒子のやりとりを書かなければならず面倒である。FDPSはその粒子のやりとりを隠蔽してくれる。ここまで来たら相互作用を書いて時間発展させるのもすぐできそうな気がする。

以下のリポジトリのstep3で`make test`すれば、コンパイルからGIF作成までやってくれるはず。

https://github.com/kaityo256/fdps_sample

[その4](http://qiita.com/kaityo256/items/2b35dac3c35d49c2eaad)に続く。