FDPSを使ってみる その2

はじめに

粒子系の並列シミュレーションフレームワークFDPSを使ってみる。
その1では、とりあえずMPIを使ってビルド、実行するところまでやった。次は系を定義し、領域分割するところまでやってみる。

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

系の定義

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

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

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

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

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

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

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

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++;
      }
    }
  }
}

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

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

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

ドメインの定義

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

  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としよう。その後、それぞれのプロセスの粒子数を表示させる。

  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()で取得できる。

実行するとこうなる。

$ 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のままである。

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

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

実行結果はこうなる。

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

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;
  }
}

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

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

実行するとこうなる。

$ 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でプロットしてみる。

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

まとめ

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

その3へ続く。

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