はじめに

すごーい！きみは粒子系の並列シミュレーションフレームワーク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;
  }