はじめに
粒子系の並列シミュレーションフレームワークFDPSを使ってみる。
その2では粒子系の定義、領域分割までやった。まずは相互作用の無い系で時間発展させてみる(基本)。
コードはここにおいてある。
https://github.com/kaityo256/fdps_sample
変数名
-
PS::ParticleSystem<FPLJ> lj_system粒子系 -
PS::DomainInfo dinfo領域情報
なお、前回はテンプレートがどのような実体を持つか明示するために
dump< PS::ParticleSystem<FPLJ> >(rank, lj_system);
などとしていたが、今回は型を省略して
dump(rank, lj_system);
となどと呼ぶことにする。
時間発展
時間発展を記述するには、運動量と座標を更新すれば良い。相互作用がなければ、単に座標を運動量で更新するだけである。記述すべき関数driftはこうなる。
template<class Tpsys>
void
drift(Tpsys & system, const PS::F64 dt) {
const PS::S32 n = system.getNumberOfParticleLocal();
for(PS::S32 i = 0; i < n; i++) {
system[i].q += system[i].p * dt;
}
}
PS::ParticleSystem::getNumberOfParticleLocal()で、「自分の領域が持っている粒子数」が取得できるので、それでループを回せば良い。引数として時間刻みdtを受け取っている。
これを単に呼べばよいが、問題は周期境界条件の補正である。全体の領域(root domainと呼ばれている)をはみ出した粒子がいる状態でドメイン間の粒子交換を行おうとするとFDPSに怒られる。
PS_ERROR: A particle is out of root domain
function: exchangeParticle, line: 773, file:
/path_to/FDPS/src/particle_system.hpp
FDPSには境界条件補正のための関数、PS::ParticleSystem::adjustPositionIntoRootDomainが用意されているので、それを呼び出してやれば良い。この関数を実行するためには、フレームワーク側に座標をいじらせる方法を提供する必要がある(粒子クラスはユーザが自由に記述できるため、フレームワーク側はどれが座標かわからないから)。これは、その2で作ったFPLJクラスにsetPosというsetterを設定すれば良い。
class FPLJ {
public:
PS::F64vec p;
PS::F64vec q;
PS::F64vec getPos() const {
return q;
}
void setPos(PS::F64vec nq) {
q = nq;
}
};
この状態で、
lj_system.adjustPositionIntoRootDomain(dinfo);
を呼べば、情報が更新される。
以上から、時間発展の1ステップは以下のようにかける。
drift(lj_system, dt);
lj_system.adjustPositionIntoRootDomain(dinfo);
dinfo.decomposeDomainAll(lj_system);
lj_system.exchangeParticle(dinfo);
領域分割は毎ステップ呼ぶ必要は無いが、まずはサンプルということで気にしないことにする。
ついでに、各プロセスごとの粒子座標をダンプする関数も修正して、あとでアニメを作れるようにしておこう。
template<class System>
void
dump(const int rank, System &sys) {
static int step = 0;
char filename[256];
sprintf(filename,"dump%02d_%03d.dat",rank,step);
step++;
std::ofstream ofs(filename);
const PS::S32 n = sys.getNumberOfParticleLocal();
for (PS::S32 i = 0; i < n ; i++) {
ofs << sys[i].q.x << " ";
ofs << sys[i].q.y << " ";
ofs << sys[i].q.z << std::endl;
}
}
std::stringstream使うのが面倒になったのでsprintfで。
これを適当なインターバルで呼び出すと、最終的に時間発展ループは以下のようになる。
const PS::F64 dt = 0.01;
const int LOOP = 1000;
const int INTERVAL = 10;
dinfo.decomposeDomainAll(lj_system);
lj_system.exchangeParticle(dinfo);
for(int i=0;i<LOOP;i++){
if(i%INTERVAL == 0){
dump(rank, lj_system);
printf("(%d/%d) n= %d\n", rank, procs, lj_system.getNumberOfParticleLocal());
}
drift(lj_system, dt);
lj_system.adjustPositionIntoRootDomain(dinfo);
dinfo.decomposeDomainAll(lj_system);
lj_system.exchangeParticle(dinfo);
}
結果
実行すると、各ステップで各プロセスの持つ粒子数を出力しながら結果をダンプする。システムサイズは前回より小さくしてL=10とし、1000粒子系とした。
$ make
$ mpirun -np 4 ./a.out
******** FDPS has successfully begun. ********
(snip)
(0/4) n= 200
(2/4) n= 200
(3/4) n= 300
(1/4) n= 300
(snip)
(0/4) n= 254
(3/4) n= 254
******** FDPS has successfully finished. ********
するとdump??_???.datが大量にできるので、それをこんなgnuplotスクリプトで処理する。
set pointsize 2
set size square
do for [i=0:99] {
set term png
set out sprintf("plot%03d.png",i)
plot sprintf("dump00_%03d.dat",i)\
,sprintf("dump01_%03d.dat",i)\
,sprintf("dump02_%03d.dat",i)\
,sprintf("dump03_%03d.dat",i)
}
これで連番pngを作ってconvertすればアニメーションGIFを作ることができる。
$ gnuplot dump.plt
$ convert plot*.png plot.gif
こうしてできたのが以下のGIF。異なるプロセスに所属する粒子が異なる色で表現されている。
ただしく計算できているっぽい。
まとめ
FDPSで相互作用無しの粒子系(要するに理想気体)を計算させてみた。シリアル処理なら簡単だが、理想気体でもMPIを使って並列化しようとすると領域間の粒子のやりとりを書かなければならず面倒である。FDPSはその粒子のやりとりを隠蔽してくれる。ここまで来たら相互作用を書いて時間発展させるのもすぐできそうな気がする。
以下のリポジトリのstep3でmake testすれば、コンパイルからGIF作成までやってくれるはず。
その4に続く。