はじめに

すごーい！きみは粒子系の並列シミュレーションフレームワークFDPSを使ってみるフレンズなんだね！

その5までで、Lennard-Jones系の並列計算ができた。これでもう実用的なプログラムを走らせることができるかもしれない。そんなわけで液滴衝突をやってみる。

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

初期条件作成

(2L ,L, L)のサイズを持つ直方体の系を用意し、x座標の左右にボールをふたつ配置してぶつけることにする。ドメインの設定はこんな感じ。

  const PS::F64 L = 40;
  PS::DomainInfo dinfo;
  dinfo.initialize();
  dinfo.setBoundaryCondition(PS::BOUNDARY_CONDITION_PERIODIC_XYZ);
  dinfo.setPosRootDomain(PS::F64vec(-L, -L*0.5, -L*0.5), PS::F64vec(L, L*0.5, L*0.5));

初期条件の作成ルーチンをどうするのが正しいのかわからないが、とりあえず一度粒子情報を持つstd::vectorに突っ込んでからシステムにコピーすることにしよう。

struct Atom {
  PS::F64vec p;
  PS::F64vec q;
};

こんなクラスを作って、

std::vector<Atom> atoms;

にひたすら情報を突っ込んで、あとでPS::ParticleSystem<FPLJ> lj_system;につっこむ。

そのために、こんな関数を作った。

template<class C>
void
put_ball(PS::F64vec c, PS::F64 r, PS::F64 xv, double density, C &atoms) {
  const double s = 1.0 / pow(density * 0.25, 1.0 / 3.0);
  const int is = static_cast<int>(2 * r / s);
  const double hs = s * 0.5;
  PS::F64vec dx(hs, 0, 0);
  PS::F64vec dy(0, hs, 0);
  PS::F64vec dz(0, 0, hs);
  for (int iz = 0; iz < is; iz++) {
    for (int iy = 0; iy < is; iy++) {
      for (int ix = 0; ix < is; ix++) {
        PS::F64vec b(ix * s - r, iy * s - r, iz * s - r);
        add_atom(b + c, c, r, xv, atoms);
        add_atom(b + c + dx, c, r, xv, atoms);
        add_atom(b + c + dy, c, r, xv, atoms);
        add_atom(b + c + dz, c, r, xv, atoms);
      }
    }
  }
}

中心座標cに、半径r、密度densityのFCCに組んだ球を作り、x軸方向にxvの初速を与える。これを二回よんでやってからコピーすれば初期条件完成。

template<class C>
int
make_conf(double density, PS::F64 L, int rank, C &lj_system) {
  std::vector<Atom> atoms;
  put_ball(PS::F64vec(-L * 0.5, 0, 0), L * 0.25, 2, density, atoms);
  put_ball(PS::F64vec(L * 0.5, 0, 0), L * 0.25, -2, density, atoms);
  const int N = atoms.size();
  if (rank == 0) {
    lj_system.setNumberOfParticleLocal(N);
    for (int i = 0; i < atoms.size() ; i++) {
      lj_system[i].id = i;
      lj_system[i].p = atoms[i].p;
      lj_system[i].q = atoms[i].q;
    }
  } else {
    lj_system.setNumberOfParticleLocal(0);
  }
  return N;
}

結果の保存

あとで可視化するために、cdview形式のファイルを吐いておく。複数のプロセスからまとめて一つのファイルに吐く必要がある。一度ルートプロセスに情報を持ってきても良いが、面倒なので、同じファイルに追記していくことにしよう。

template<class C>
void
export_cdview(int rank, int procs, C &system) {
  static int index = 0;
  char filename[256];
  sprintf(filename, "conf%03d.cdv", index);
  if (0 == rank) {
    std::ofstream ofs(filename);
    std::cout << filename << std::endl;
    ofs.close();
  }
  for (int r = 0; r < procs; r++) {
    PS::Comm::barrier();
    if (rank == r) {
      std::ofstream ofs(filename, std::ios::app);
      const PS::S32 n = system.getNumberOfParticleLocal();
      for (PS::S32 i = 0; i < n; i++) {
        ofs << i << " ";
        ofs << rank << " ";
        ofs << system[i].q.x << " ";
        ofs << system[i].q.y << " ";
        ofs << system[i].q.z << " ";
        ofs << std::endl;
      }
    }
  }
  index++;
}

プロセスの数だけループを回し、自分の番が来た時にファイルを出力する。MPIプログラムで手抜きでファイルを出力したい時によく使う。cdviewには粒子の種類を指定するカラムがあるが、そこにプロセス番号を指定している。これで、異なるプロセスが担当する粒子が異なる色で表示される。

時間発展

初期条件を作ってしまえば、時間発展はステップ5までと同じ。メインループはこんな感じになろう。

  const PS::F64 dt = 0.002;
  PS::F64 s_time = 0.0;
  const int LOOP = 10000;
  const int INTERVAL = 100;
  dinfo.decomposeDomainAll(lj_system);
  lj_system.exchangeParticle(dinfo);
  export_cdview(rank, procs, lj_system);

  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);
    lj_system.adjustPositionIntoRootDomain(dinfo);
    dinfo.decomposeDomainAll(lj_system);
    lj_system.exchangeParticle(dinfo);
    if (i % INTERVAL == 0) {
      export_cdview(rank, procs, lj_system);
    }
    s_time += dt;
  }

二次のシンプレクティック積分をして、たまにcdview形式で保存している。

結果

以上を4プロセスで実行してみた。初期条件はこんな感じ。2つのプロセスで1つのボールを担当するのが理想だが、もう一つのボールの粒子も少しだけ含まれてしまっている。

ただし、この状態は数ステップ後に解消され、ちゃんと2つのプロセスで一つのボールを担当するようになる。

衝突したところ

ぶつかって扁平な形になる。

アニメーションGIFにしたもの。

なお、可視化はcdviewコンパチなコード、Avogadro Challenge Viewer, acvを用いた¹。

まとめ

FDPSを使って相互作用のある粒子系を定義し、　液滴衝突の並列計算をしてみた。とりあえず1プロセスより2プロセス、2プロセスより4プロセスの方が早いことは確認したが、それ以上は速度がサチるようだ。

上記のコードは以下のリポジトリのstep6においてある。実行するとcdview形式のファイルを連番で吐くので、興味のある人は可視化を試みられたい。

というわけで「FDPSを使ってみる」はおしまい。お疲れさまでした＞僕

Macでのビルドがacvの方が楽だったため。 ↩