FDPSを使ってみる その5

# はじめに

粒子系の並列シミュレーションフレームワーク[FDPS](https://github.com/FDPS/FDPS)を使ってみる。
[その4](http://qiita.com/kaityo256/items/2b35dac3c35d49c2eaad)ではカットオフのあるLennard-Jones粒子系の二体衝突のシリアル計算を行った。次は全く同じ計算の並列版を実行し、軌道、エネルギーが一致するか確認する(基本)。

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

# 並列化にあたっての修正

まずはmakefileに、FDPSの並列化マクロなどを追加指定する。

```mf
CPPFLAGS += -DPARTICLE_SIMULATOR_MPI_PARALLEL
LDFLAGS=-lmpi -lmpi_cxx
```

フレームワークの初期化直後にMPIのrankとプロセス数も取得しておこう。

```cpp
  PS::Initialize(argc, argv);
  const int rank = PS::Comm::getRank();
  const int procs = PS::Comm::getNumberOfProc();
```

粒子を追加するにあたって、とりあえずプロセス0に突っ込んで、粒子の分配はフレームワークに任せることにする。

```cpp
  const int N = 2;
  if (rank == 0){
    lj_system.setNumberOfParticleLocal(N);
    lj_system[0].id = 0;
    lj_system[0].p = PS::F64vec(1.0, 0.0, 0.0);
    lj_system[0].q = PS::F64vec(2.5, 5.0, 5.0);
    lj_system[1].id = 1;
    lj_system[1].p = PS::F64vec(-1.0, 0.0, 0.0);
    lj_system[1].q = PS::F64vec(7.5, 5.0, 5.0);
  }else{
    lj_system.setNumberOfParticleLocal(0);
  }
  dinfo.decomposeDomainAll(lj_system);  // 領域分割
  lj_system.exchangeParticle(dinfo); // 粒子の分配
```

エネルギーの計算だが、既に二体関数であるポテンシャルエネルギーを一体関数として分配してあるため、AllReduceするだけで良い。そのラッパーである`PS::Comm::getSum`を呼ぶと、エネルギーを返す関数は

```cpp
template<class Tpsys>
PS::F64
energy(const Tpsys &system){
  PS::F64 e = 0.0;
  const PS::S32 n = system.getNumberOfParticleLocal();
  for(PS::S32 i = 0; i < n; i++) {
    FPLJ a = system[i];
    e += (a.p*a.p)*0.5;
    e += a.potential;
  }
  return PS::Comm::getSum(e); //←ここを修正
}
```

と書ける。

#　結果の保存

各プロセスごとに、自分が持っている粒子を個別に吐く関数を作っておく。

```cpp
template<class Tpsys>
void
dump(std::ostream &os, const PS::F64 s_time, Tpsys &system){
  const PS::S32 n = system.getNumberOfParticleLocal();
  if (n == 0)return;
  os << s_time << " ";
  for (PS::S32 i = 0; i < n; i++) {
    os << system[i].q.x << " ";
  }
  os << std::endl;
}
```

# メインループ

計算のメインループはこんな感じになる。

```cpp
  char filename[256];
  sprintf(filename,"trac%02d.dat",rank);
  std::ofstream ofs(filename);

  for(int i=0;i<LOOP;i++){
    drift(lj_system, dt*0.5);
    tree_lj.calcForceAllAndWriteBack(CalcForceLJ(), lj_system, dinfo);
    kick(lj_system, dt);
    drift(lj_system, dt*0.5);
    dump(ofs, s_time, lj_system);
    tree_lj.calcForceAllAndWriteBack(CalcForceLJ(), lj_system, dinfo);
    PS::F64 e = energy(lj_system);
    if (rank==0){
      std::cout << s_time << " " << e << std::endl;
    }
    s_time += dt;
  }
```

結果の出力まわり以外、シリアル版とほとんどかわらない。なお、エネルギーの計算ルーチン`energy`の内部でMPI_Allreduceを呼んでいるため、

```cpp
    if (rank==0){
      std::cout << s_time << " " << energy(lj_system) << std::endl;
    }
```

みたいなことをするとデッドロックする[^2]。なので一度外で変数に受けてから表示するようにしている。

# 結果

1プロセスでシリアル実行した場合と、2プロセスで並列実行した場合の粒子の軌道は以下の通り。

![image](https://qiita-image-store.s3.amazonaws.com/0/79744/ef3acc0c-c9f3-df99-5b4c-c59abb935481.png)

それぞれのプロセスが1粒子ずつ担当し、かつシリアル実行した場合と軌道が一致していることがわかる。

エネルギーの時間発展はこんな感じ。

![image](https://qiita-image-store.s3.amazonaws.com/0/79744/4acd2be6-1985-11e6-bd75-ab6e7b040eb5.png)

当然のことながら、完全に一致している。

# まとめ

FDPSを使って相互作用のある粒子系を定義し、二粒子衝突の計算を並列化してみた。シリアル版がかけてしまえば、並列版に修正するのは容易だと思う。

以下のリポジトリのstep5でmake testすれば、コンパイルして、シリアル実行、2並列実行した結果の比較グラフまで作ってくれるはず。

https://github.com/kaityo256/fdps_sample

[その6](http://qiita.com/kaityo256/items/95e7a8076622ad27aac3)に続く。

[^2]: ・・・というようなことを意識しなくて良いフレームワークが設計できれば、より美しいと思う。