複数プロセスにまたがる座標データから密度プロファイルを作成

これはMPI Advent Calendar 2017の19日目の記事です。気がついたらMPI AdCの7本目の記事となりました。この辺でもう遠慮しようと思います。

はじめに

「MPIを使って非自明並列をする」というと、OpenMPによる共有メモリ並列や、PGASを用いた仮想的な共有メモリ並列に比べて難しい、面倒くさい、というイメージがあります。正直「面倒くさい」というところは否定しませんが、難しくはありません。むしろ「必要なデータは明示的に通信しなければならない」という縛りは、共有メモリモデルの並列化に比べて圧倒的にエンバグ確率を減らすと思います。また、バグった時のデバッグも、スレッド並列のタイミングに依存するバグに比べて、プロセス間通信のデバッグの方が通信のタイミングが明示的にわかる分だけ楽です。っていうか、マルチスレッドで「あるデータが特別な状態になっている時にスレッドが割り込んできた時にのみ起きるバグ」みたいな奴は、組む前から気をつけるのも難しいし、この手のバグを入れてしまった場合、デバッグは正直地獄です。

とは言え、(うまくいけば)ディレクティブ一発で並列化できちゃう他の手法に比べて、MPIで通信を明示的に書くのは面倒なことは確かでしょう。また、「MPIで並列プログラムを書く」といっても、自明並列(いわゆる馬鹿パラ)では簡単すぎ、かといっていきなりガチ並列というのも難しすぎます。

そんなわけで、少しだけ非自明で、でもなるべく簡単な、それでいて実践的なMPIプログラムの例として、領域分割によりプロセスごとに異なる座標情報を持った状態から、系全体の密度場を計算するプログラムを組んでみたいと思います。

ソースコードは以下においておきます。
https://github.com/kaityo256/density_reduce

シリアル版

分子動力学法(Molecular Dynamics simulation, MD)では、初期条件として原子の位置と速度を受け取り、時間発展をさせ、適当なインターバルで観測を行ってそれをファイルに保存する、ということを考えます。

これを模して、初期条件として原子の座標を受け取り、時間発展をさせずに密度場を保存するコードを書いてみましょう。その際、空間をグリッドに切って、それぞれのセル内の数密度を密度場として保存することにします。さらに、その密度場データをVTKファイルに変換して、ParaViewで可視化することを考えましょう。密度を保存する際にいきなりVTKファイルにしても良いですが、ファイルサイズが大きくなってしまうため、とりあえず生でグリッドデータを吐いておいて、後でVTKに保存するほうがいろいろ便利です。

というわけで、こんな感じのプログラム群を書きます¹。

最初にmakedump.rbで初期条件を作ります。原子配置を作って、doubleの(x,y,z)データをずらずら並べたtest.dumpというファイルを作るスクリプトです。

次にdump2dat.cppで、初期条件を受け取って、密度場のデータファイルtest.datを作ります。これがMDシミュレーションを模したもの(シミュレーションのシミュレーション？)となります。

それ以降はシミューレーションのポストプロセスとなります。MDを行って作成された大量の(今回は一つだけですが)datファイルを、可視化のためにVTKファイルに変換し、最終的にParaViewで可視化します。

まずは並列化をしない、シリアル版のコードを書いてみましょう。どんなに簡単なコードでも、まず等価な動作をするシリアル版を組むのは並列プログラミングの基本です。

それぞれ特に難しくないので、スクリプトやソースを載せます。

makedump.rb

L = 100
N = 100000
a = []
C = L*0.5
R = L*0.4
loop do
  theta = 2.0*Math::PI*rand
  r = rand**(1.0/3.0)
  z = rand*2.0 - 1.0
  x = C + R*r*Math::sqrt(1-z**2)*Math::cos(theta)
  y = C + R*r*Math::sqrt(1-z**2)*Math::sin(theta)
  z = C + R*r*z
  a.push x,y,z
  break if a.size >= N*3
end

IO.binwrite("test.dump", a.pack("d*"))
puts "Saved as test.dump"

100×100×100の空間に、球状に一様に10万個の原子をばらまくスクリプトです。

dump2dat.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <assert.h>

const double L = 100.0;
const int G = 51;

struct Atom {
  double x, y, z;
};

void
make_dat(const std::string filename, std::vector<Atom> &atoms) {
  std::vector<int> num(G * G * G, 0);
  const double s = L / G;
  for (auto a : atoms) {
    int ix = static_cast<int>(a.x / s);
    int iy = static_cast<int>(a.y / s);
    int iz = static_cast<int>(a.z / s);
    int index = ix + iy * G + iz * G * G;
    assert(index < G * G * G);
    num[index]++;
  }
  const double v_inv = 1.0 / s / s / s;
  std::vector<double> density(G * G * G);
  for (int i = 0; i < G * G * G; i++) {
    density[i] = static_cast<double>(num[i]) * v_inv;
  }
  std::ofstream ofs(filename);
  ofs.write((char*)density.data(), sizeof(double)*density.size());
  std::cout << "Saved as " << filename << std::endl;
}

void
load_dump(const std::string filename, std::vector<Atom> &atoms) {
  std::ifstream ifs(filename);
  ifs.seekg(0, ifs.end);
  int size = ifs.tellg() / sizeof(double) / 3;
  ifs.seekg(0, ifs.beg);
  atoms.resize(size);
  for (int i = 0; i < size; i++) {
    ifs.read((char*)(&atoms[i].x), sizeof(double));
    ifs.read((char*)(&atoms[i].y), sizeof(double));
    ifs.read((char*)(&atoms[i].z), sizeof(double));
  }
}

int
main(void) {
  std::vector<Atom> atoms;
  load_dump("test.dump", atoms);
  make_dat("test.dat", atoms);
}

粒子データを受け取って、密度場を作るプログラムです。51×51×51のグリッドで空間を区切って、それぞれのセルに入っている原子の数を数え、その数密度をtest.datとして保存します。グリッドの分割数がキリの悪い数字になっているのはわざとです。

dat2vtk.rb

require 'pathname'

if ARGV.size == 0
  puts "usage: ruby dat2vtk.rb datfile"
  exit
end
inputfile = ARGV[0]

data = IO.binread(inputfile)
density = data.unpack("d*")
dim = 51

outputfile = Pathname(inputfile).sub_ext(".vtk").to_s

f = open(outputfile,"w")
f.puts <<"EOS"
# vtk DataFile Version 1.0
density
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS #{dim} #{dim} #{dim}
ORIGIN 0.0 0.0 0.0
SPACING 1.0 1.0 1.0
POINT_DATA #{dim**3}
SCALARS intensity float
LOOKUP_TABLE default
EOS
density.each do |d|
  f.puts d
end

puts "#{inputfile} -> #{outputfile}"

密度場をVTKファイルに変換するスクリプトです。VTKレガシーフォマットはヘッダをつけて生データをずらずら並べるだけなのでとても楽です。

以上、一連の動作をやってみましょう。

$ ruby makedump.rb  # 初期条件ファイル作成
Saved as test.dump
$ g++ dump2dat.cpp; ./a.out # 密度場ファイルの作成
Saved as test.dat
$ ruby dat2vtk.rb test.dat # VTKファイル作成
test.dat -> test.vtk

できたtest.vtkをParaViewで読み込んで「Volume」でボリュームレンダリングするとこんな感じになります。

それっぽい図ができました。

MPI並列版

さて、先程のコードのMPI並列版を考えます。並列化手法としては空間分割を採用し、並列化するのは時間発展のところのみとしましょう。つまり、各プロセスが別々の空間を割り当てられており、それぞれに原子の情報を持っているとします。この状態から全体の密度場を計算します。並列化するのはdump2dat.cppのところだけです。これを並列化したdump2dat_mpi.cppを作ります。

原子配置の読み込み

まず、原子配置を各プロセスに読み込むことを考えましょう。原子の読み込みはシミュレーションの最初に一度だけ行われる処理であるため、実行時間は気にしないことにします。なので、

• 全てのプロセスが同じ初期条件ファイルを読み込みに行く
• 自分の担当領域にいる原子だけ、自分のvectorに追加する

という方針でいきます。そのために、3次元の領域を担当するクラスを作って置くと便利です。ここではMyrectという手抜きクラスを実装します。

struct Myrect {
  double sx, sy, sz;
  double ex, ey, ez;
  Myrect(int rank, int d[3], double L) {
    int ix = rank % d[0];
    int iy = (rank / d[0]) % d[1];
    int iz = (rank / d[0] / d[1]);
    double gx = L / d[0];
    double gy = L / d[1];
    double gz = L / d[2];
    sx = gx * ix;
    ex = sx + gx;
    sy = gx * iy;
    ey = sy + gy;
    sz = gz * iz;
    ez = sz + gz;
  }
  bool is_inside(double x, double y, double z) {
    if (x < sx)return false;
    if (x >= ex)return false;
    if (y < sy)return false;
    if (y >= ey)return false;
    if (z < sz)return false;
    if (z >= ez)return false;
    return true;
  }
};

クラスMyrectは、コンストラクタで自分のランク、分割数、シミュレーションボックスの一辺の長さ(立方体であることを仮定しています)を受け取り、自分の担当領域の直方体、(sx, sy, sz) - (ex, ey, ez)を計算します。その後、それを使って、ある座標が自分の内部にあるかどうか判定するis_insideを実装しています。

分割数は、例えばプロセス数が12の場合には3×2×2みたいに、どのようにプロセス数を分割するかを表します。ここではMPI_Dims_createの分割をそのまま採用しています。

  int rank, size;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
  int d[3] = {};
  MPI_Dims_create(size, 3, d);
  Myrect rect(rank, d, L);

密度場の計算

もし、以下の図の様にグリッドの分割ラインとMPIによる空間分割がちょうど一致していたら、それぞれのプロセスにおいてそれぞれのセルの密度場を計算し、後でうまくまとめれば良いでしょう。

しかし、MPIの空間分割ラインと局所物理量の計算のためのセルが一致していない時にはどうすればよいでしょうか。

プロセスが担当する領域が中途半端にまたがったセルが存在するため、面倒なことになります。

もちろん、「常にプロセスの分割数とグリッドの分割数が整合するように考える」というのも一つの解でしょう。しかし、後々の処理のためにグリッドの分割数は2のベキで行いたいが、MPIのプロセス数はどうしても2のベキで取れない、ということがよく発生します。

そこで、富豪的プログラミングをすることにしましょう。つまり、全てのプロセスが、全空間のセル情報を持っておき、そこに自分の情報だけ書き込み、後でそのセル情報をAllreduceすることにします。密度場、速度場、圧力場、温度場といった、通常興味のある物理量は、ほとんど「全部足してから何か(多くの場合原子数)で平均」という処理ができます。今回も、密度場の計算のため、各プロセスは各セルに何個の原子がいるかをstd::vector<int>に保存し、後でAllreduceしてしまうことにしましょう。

並列版のコード

以上の方針で実装した並列版のコードはこんな感じになります。

dump2dat_mpi.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <assert.h>
#include <mpi.h>

const double L = 100.0;
const int G = 51;

template<class T>
void
myallreduce(void *s, void *r, int count) {
  printf("Unsupported Type\n");
  MPI_Abort(MPI_COMM_WORLD, 1);
}

template<> void myallreduce<int>(void *s, void *r, int count) {
  MPI_Allreduce(s, r, count, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
}

template <class T>
void
AllreduceVector(std::vector<T> &sendbuf, std::vector<T> &recvbuf) {
  const int count = sendbuf.size();
  recvbuf.resize(count);
  std::fill(recvbuf.begin(), recvbuf.end(), 0);
  myallreduce<T>(sendbuf.data(), recvbuf.data(), count);
}

struct Atom {
  double x, y, z;
  Atom(double _x, double _y, double _z) {
    x = _x;
    y = _y;
    z = _z;
  }
};

struct Myrect {
  double sx, sy, sz;
  double ex, ey, ez;
  Myrect(int rank, int d[3], double L) {
    int ix = rank % d[0];
    int iy = (rank / d[0]) % d[1];
    int iz = (rank / d[0] / d[1]);
    double gx = L / d[0];
    double gy = L / d[1];
    double gz = L / d[2];
    sx = gx * ix;
    ex = sx + gx;
    sy = gx * iy;
    ey = sy + gy;
    sz = gz * iz;
    ez = sz + gz;
  }
  bool is_inside(double x, double y, double z) {
    if (x < sx)return false;
    if (x >= ex)return false;
    if (y < sy)return false;
    if (y >= ey)return false;
    if (z < sz)return false;
    if (z >= ez)return false;
    return true;
  }
};

void
make_dat(const std::string filename, std::vector<Atom> &atoms, int rank) {
  std::vector<int> num(G * G * G, 0);
  const double s = L / G;
  for (auto a : atoms) {
    int ix = static_cast<int>(a.x / s);
    int iy = static_cast<int>(a.y / s);
    int iz = static_cast<int>(a.z / s);
    int index = ix + iy * G + iz * G * G;
    assert(index < G * G * G);
    num[index]++;
  }
  std::vector<int> num_recv;
  AllreduceVector(num, num_recv);
  const double v_inv = 1.0 / s / s / s;
  std::vector<double> density(G * G * G);
  for (int i = 0; i < G * G * G; i++) {
    density[i] = static_cast<double>(num_recv[i]) * v_inv;
  }
  if (rank == 0) {
    std::ofstream ofs(filename);
    ofs.write((char*)density.data(), sizeof(double)*density.size());
    std::cout << "Saved as " << filename << std::endl;
  }
}

void
load_dump(const std::string filename, std::vector<Atom> &atoms, Myrect &rect) {
  std::ifstream ifs(filename);
  ifs.seekg(0, ifs.end);
  int size = ifs.tellg() / sizeof(double) / 3;
  ifs.seekg(0, ifs.beg);
  for (int i = 0; i < size; i++) {
    double x, y, z;
    ifs.read((char*)(&x), sizeof(double));
    ifs.read((char*)(&y), sizeof(double));
    ifs.read((char*)(&z), sizeof(double));
    if (rect.is_inside(x, y, z)) {
      atoms.push_back(Atom(x, y, z));
    }
  }
}

void
save_txt(int rank, std::vector<Atom> &atoms) {
  char filename[256];
  sprintf(filename, "rank%02d.txt", rank);
  std::ofstream ofs(filename);
  for (auto &a : atoms) {
    ofs << a.x << " " << a.y << " " << a.z << std::endl;
  }
}

int
main(int argc, char **argv) {
  MPI_Init(&argc, &argv);
  int rank, size;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
  int d[3] = {};
  MPI_Dims_create(size, 3, d);
  Myrect rect(rank, d, L);
  std::vector<Atom> atoms;
  load_dump("test.dump", atoms, rect);
  save_txt(rank, atoms);
  make_dat("test_mpi.dat", atoms, rank);
  MPI_Finalize();
}

処理の流れはこんな感じです。

• 自分の領域の確認(Myrect rect(rank, d, L);)
• load_dumpで、自分の領域にいる原子だけ読み込む
• make_datで、シリアル版と同様にセル内の原子数を数えてnumsに保存
• nums_recvにnumsをAllreduceした結果を受け取り、その情報を使って密度場を計算

結局、通信はAllreduce一回だけですね。

実行結果

まず、読み込みがうまくいったかどうかの確認のため、各プロセスが担当する原子をrank%02d.txtというファイル名で保存しています。例えば4分割だとこんな感じです。

$ g++ dump2dat_mpi.cpp -lmpi -lmpi_cxx -o b.out
$ mpirun --oversubscribe -np 4 ./b.out
Saved as test_mpi.dat
$ gnuplot
gnuplot>  sp "rank00.txt","rank01.txt","rank02.txt","rank03.txt"

ちゃんと四分割されているようですね。

また、並列版で作成された密度場のファイルが、シリアル版と一致しているか確認しましょう。

$ diff test.dat test_mpi.dat

ちゃんと一致しています。バイナリが同じなので、これをVTKに変換して可視化しても、全く同じ画像を得ることができます。

まとめ

並列MDシミュレーションにおける、密度場の保存＆可視化をやってみました。並列MDコードを持っている人は、上記のコードを少し改造すれば簡単にVTK用のファイルを作成することができるでしょう。

MPIによる領域分割コードは確かに面倒です。ですが、短距離古典MDに関して言うならば、領域の「のりしろ」の通信と、今回のような全体通信くらいしか発生しないので、さほど難しいことは発生しません。ただ面倒なだけです。もちろん、並列化によって入ってくるバグもあると思いますが、今回のように「同じ動作をするシリアルコードを書く」という作業を徹底すれば、並列化に起因するバグはほとんど入りません²。

個人的な印象として、MPIによる並列プログラムの難しさは、最初からなるべく効率的な並列化コードを目指して破綻することに起因するのではないかと思っています。僕も、今回のようなコードを以前書いた時には、「各プロセスが全空間の状態を持っていたらメモリ的に破綻するのではないか」と考えて、グリッドセルと領域分割を必ず整合させて、自分の担当部分だけを通信して・・・とかやっていました。でも、よく考えたら、結局ルートプロセスは全空間の状態を保持することになるわけで、ルートプロセスがメモリ的に保持できるなら、他のプロセスだって保持できるのですよね(多くの場合、ルートプロセスが最もメモリを食うから)。なので、メモリ使用量がカツカツで、flat-MPIでノード内に多数のプロセスを立ち上げてる、とかいう状況でなければ、今回みたいに「全プロセスで全空間の情報を持って、後でallreduceしちゃう」方式の方がいろいろ楽です³。

MPIのAdCの記事を何本も書いておいてなんですが、僕はMPIという仕組みがあまり好きではありません。好きではありませんが、PGASやらOpenMPやらで(仮想)共有メモリで凝ったことをやるくらいならMPIで明示的に通信を書いてしまった方がいろいろ楽だと思ってます。でもまぁ、そのうちBetter MPI、もしくはMPI的なことをネイティブに仕様に取り込んだプログラム言語が出てこないかなぁ、と夢想中です⁴。っていうかディレクティブ嫌い。

参考

• MPIでいい感じに因数分解 MPI_Dims_createの説明
• Boost MPI Libraryでstd::vectorをall_reduceする Boost MPIを使ったり使わなかったりしてstd::vectorをAllreduce

---

1. 今気がついたけど、パワーポイントをそのままキャプチャしたから、校正の赤線入っちゃってる・・・ ↩

2. ただし、規模によって顕在化するバグは防げません。たとえば数十億分の一の確率でで発生するバグとか、ローカルでテストしている時にはスルーしても、大規模並列時にはほぼ確実に踏んでしまうとか、あまりメモリを気にせず書いたコードが、プロセス数に比例してメモリ消費が増え、数万プロセスで破綻するとか・・・。 ↩

3. まぁ、「メモリ使用量がカツカツで、かつflat-MPIでノード内にプロセスを立ち上げている」状況って、わりとあるわけですが・・・。 ↩

4. まぁ、本当に欲しければ自分で作れよ、という話なんだけどさ・・・。それはわかってるんだけどさ・・・ ↩