FDPSを使ってみる その2

# はじめに

粒子系の並列シミュレーションフレームワーク[FDPS](https://github.com/FDPS/FDPS)を使ってみる。
[その1](http://qiita.com/kaityo256/items/c4955adaad718fd399bc)では、とりあえずMPIを使ってビルド、実行するところまでやった。次は系を定義し、領域分割するところまでやってみる。

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

# 系の定義

FDPSでは、粒子の集まりを`PS::ParticleSystem`で管理する。これに粒子を定義したクラスを入れて実体化させる。ここではLennard-Jones粒子系を想定し、こんな定義を書いてみる。

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

`PS::F64vec`はデフォルトで3次元ベクトルを表現する[^1]。とりあえず運動量`p`と座標`q`だけ定義しておく。あとで領域分割をする際に、ユーザーが定義した変数のうちどれが座標を表すかを教えるため、`getPos()`という関数を定義する必要がある。粒子の最低限の定義はこれで良い。

さて、このクラスを`PS::ParticleSystem`につっこんで実体化させる。

```cpp
  PS::ParticleSystem<FPLJ> lj_system;
  lj_system.initialize();
```

これで粒子数0の粒子系が実体化された。ユーザ側からはとりあえず`std::vector<FPLJ>`みたいなものと思っておけば良い。

これに粒子を追加するには、まず`setNumberOfParticleLocal`で粒子数を教えて、その後`[]`演算子で値を入力すれば良い。`std::vector::resize`だと思えば良いと思う。例えば、一辺の長さ`L`の系に、数密度1で立方格子状に粒子を配置するコードはこんな感じになる。

```cpp
template<class System>
void
make_conf(System &sys, const PS::F64 L) {
  std::mt19937 mt(1);
  std::uniform_real_distribution<double> ud(0.0, 1.0);
  const double V0 = 1.0;
  const int il = static_cast<int>(L);
  PS::S32 n_total = L * L * L;
  sys.setNumberOfParticleLocal(n_total);
  int n = 0;
  for (int ix = 0; ix < L; ix++) {
    for (int iy = 0; iy < L; iy++) {
      for (int iz = 0; iz < L; iz++) {
        double z = ud(mt) * 2.0 - 1.0;
        double phi = ud(mt) * M_PI;
        sys[n].p.x = V0 * sqrt(1 - z * z) * cos(phi);
        sys[n].p.y = V0 * sqrt(1 - z * z) * sin(phi);
        sys[n].p.z = V0 * z;
        sys[n].q.x = ix + 0.5;
        sys[n].q.y = iy + 0.5;
        sys[n].q.z = iz + 0.5;
        n++;
      }
    }
  }
}
```

ついでに運動量もランダムに与えている。この関数を

```cpp
make_conf< PS::ParticleSystem<FPLJ> >(lj_system, L);
```

として呼んでやれば、系に粒子が追加される。

# ドメインの定義

FDPSにはドメインを管理するクラス`PS::DomeinInfo`がある。まずはシステムサイズ`L`、周期的境界条件のドメインを定義する。

```cpp
  const PS::F64 L = 20;
  PS::DomainInfo dinfo;
  dinfo.initialize();
  dinfo.setBoundaryCondition(PS::BOUNDARY_CONDITION_PERIODIC_XYZ);
  dinfo.setPosRootDomain(PS::F64vec(0, 0, 0), PS::F64vec(L, L, L));
```

さいしょに初期化し、境界条件を与えるところまでは良いと思う。最後の`PS::DomainInfo::setPosRootDomain`は系の直方体を与えるところで、二つのベクトルにより定義する。

# 領域分割によるロードバランス

まずは何もせず、プロセス0にのみ粒子を配置し、残りは0としよう。その後、それぞれのプロセスの粒子数を表示させる。

```cpp
  PS::Initialize(argc, argv);
  const int rank = PS::Comm::getRank();
  const int procs = PS::Comm::getNumberOfProc();
  PS::ParticleSystem<FPLJ> lj_system;
  lj_system.initialize();

  const PS::F64 L = 20;
  PS::DomainInfo dinfo;
  dinfo.initialize();
  dinfo.setBoundaryCondition(PS::BOUNDARY_CONDITION_PERIODIC_XYZ);
  dinfo.setPosRootDomain(PS::F64vec(0, 0, 0), PS::F64vec(L, L, L));
  if (rank == 0) {
    make_conf< PS::ParticleSystem<FPLJ> >(lj_system, L);
  } else {
    lj_system.setNumberOfParticleLocal(0);
  }
  //あとでここに領域分割と粒子交換を追加する
  PS::Comm::barrier();
  printf("(%d/%d) n= %d\n", rank, procs, lj_system.getNumberOfParticleLocal());
  PS::Finalize();
```

各プロセスの粒子数は`PS::ParticleSystem<FPLJ>::getNumberOfParticleLocal()`で取得できる。

実行するとこうなる。

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

当然のことながらプロせス0番が全ての粒子をもっており、残りは0のままである。

さて、この状態で領域分割と粒子の再配置をさせる。粒子追加処理の直後に以下の二行を追加すれば良い。

```cpp
  dinfo.decomposeDomainAll(lj_system);
  lj_system.exchangeParticle(dinfo);
```

実行結果はこうなる。

```shell-session
$ mpirun -np 4 ./a.out
(snip)
******** FDPS has successfully begun. ********
(0/4) n= 2200
(1/4) n= 1800
(2/4) n= 1200
(3/4) n= 2800
******** FDPS has successfully finished. ********
```

粒子が4プロセスに再配置された。どのように分配されたか見るため、各プロセスの持つ粒子をファイルに吐いてみる。そういうことをするためのフレームワークも用意されているっぽいが、ここでは自前で書いてしまう。

```cpp
template<class System>
void
dump(const int rank, System &sys) {
  std::stringstream ss;
  ss << "dump" << rank << ".dat";
  std::ofstream ofs(ss.str());
  const int n = sys.getNumberOfParticleLocal();
  for (int i = 0; i < n ; i++) {
    ofs << sys[i].q.x << " ";
    ofs << sys[i].q.y << " ";
    ofs << sys[i].q.z << std::endl;
  }
}
```

これを、粒子交換直後に以下のように呼ぶ。

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

実行するとこうなる。

```shell-session
$ mpirun -np 4 ./a.out    
(snip)
******** FDPS has successfully begun. ********
(0/4) n= 2200
(2/4) n= 1200
(3/4) n= 2800
(1/4) n= 1800
******** FDPS has successfully finished. ********
```

粒子が分割された。dump0.dat、dump1.dat・・・が作られているのでgnuplotでプロットしてみる。

![image](https://qiita-image-store.s3.amazonaws.com/0/79744/c2631570-e72b-4e29-7aeb-60800ac6fdf4.png)

系は三次元だが、4プロセスだとxy平面で分割される。何も設定せずにつっこんだのでロードバランスがやや悪いが、適切に設定すればちゃんと分割されると思われる。

# まとめ

粒子系を定義し、ドメイン分割をさせてみた。きちんとロードバランスを取るにはいろいろノウハウが必要なのだろうが、面倒な粒子の再配置やドメインの定義が数行で実現できるのは良い感じ。

[その3](http://qiita.com/kaityo256/items/36d267be911ac45e848e)へ続く。

[^1]: コンパイルオプションにより2次元にすることもできる。