というわけで、一応木構造をつくって重力計算をするためのデータ構造、アル ゴリズムとその実現をきわめて簡単に紹介した。これを実際に多体問題の数値 解を求めるプログラムにするためには、時間積分、初期条件の生成／読み込み、 結果の表示／解析等いろんなことをするプログラムを作る必要がある。
それらに入る前に、とりあえず今まで出来たところを動かしてみることにする。
割合いろんな関数を作っていろんなことをやらせるので、全体を動かして時に ちゃんと動いているかどうかをどうやって判断し、さらに動いていない時にど こがおかしいかを見つけ出すのはそんなに簡単ではなくなる。
ある程度大きな、沢山の関数からできたプログラムをテスト・デバッグするた めの基本的な考え方は以下の 2 つである。
もちろん、プログラムをデバッグする時の大きな問題は、
ということである。特に後者がデバッグを難しいものにする。
というわけで、まず、ツリー構造を作り、それが正しいかどうか調べる。それから実際に重力を計算し、それもツリー構造を使わないで総当りで計算したものとくらべてみよう。メインプログラムは以下のような感じになる。
pt minus 1 pt
#ifdef TEST
int main()
{
particle * pp;
int n;
cerr << "Enter n:";
cin >> n ;
pp = new particle[n];
double rsize = 1.0;
create_uniform_sphere(pp, n, 0 , rsize);
for(int i =0;i <n;i++){
PRC(i); PRL(pp[i].pos);
}
bhnode * bn = NULL;
int nnodes = n*2;
bn = new bhnode[nnodes];
bn->assign_root(vector(0.0), rsize*2, pp, n);
int heap_remainder = nnodes-1;
bhnode * btmp = bn + 1;
bn->create_tree_recursive(btmp, heap_remainder);
PRL(bn->sanity_check());
bn->set_cm_quantities();
bn->dump();
double eps2 = 0.01;
calculate_uncorrected_gravity_direct(pp, n, eps2);
cerr << "Direct force \n";
particle * p = pp;
for(int i = 0; i<n; i++){
PR(i); PR(p->pos);PR(p->phi); PRL(p->acc);
p++;
}
cerr << "Tree force \n";
clear_acc_and_phi(pp, n);
p = pp;
for(int i = 0; i<n; i++){
calculate_gravity_using_tree(p, bn, eps2, 0.5);
PR(i); PR(p->pos);PR(p->phi); PRL(p->acc);
p++;
}
return 0;
}
#endif
プログラム全体は
http://grape.astron.s.u-tokyo.ac.jp/%7Emakino/kougi/keisan_tenmongakuII/programs/simpletree/BHtree.C
にあるので細かい関数については省略するが、関数 sanity_check は ちょっと中を見ておこう。
pt minus 1 pt
int bhnode::sanity_check()
{
int i;
int iret = 0;
if (nparticle == 1 ){
// this is the lowest level node. Things to check:
// particle is in the cell
particle * p = pfirst;
if(inbox(pos,p->pos,l)){
cerr << "Error, particle out of box ... \n";
dump();
return 1;
}
}else{
// This is the non-leaf node. Check the position and side
// length of the child cells and then check recursively..
for(i=0;i<8;i++){
if (child[i] != NULL){
int err = 0;
err = child[i]->sanity_check();
if (l*0.5 != child[i]->l) err += 2;
vector relpos = cpos-child[i]->cpos;
for (int k = 0 ; k<ndim;k++){
if (fabs(relpos[k]) !=l*0.25)err += 4;
}
if (err){
cerr << "Child " << i << " Error type = " << err << endl;
dump();
}
iret += err;
}
}
}
return iret;
}
ここでは、以下のようなチェックを行なっている。
inbox はこんなものである。
pt minus 1 pt
int inbox(vector & cpos, // center of the box
vector & pos, // position of the particle
real l) // length of one side of the box
{
for(int i = 0; i< ndim; i++){
if (fabs(pos[i]-cpos[i]) > l*0.5) return 1;
}
return 0;
}
ここではベクトルの各要素についてそれぞれ調べる必要がある。
さて、この関数が「エラー」を返して来たとしても、少なくとも以下の4通りの可能性が ある。
というわけで、最初は粒子 1 つとか 2 つで何ができるかみるといったところ から始めないと後が大変ではある。逆に数が少なければいろいろ書き出して人 間がチェックするのも不可能ではない。 以下の関数はツリー構造そのものを再帰的にプリントアウトする。 spc は単に引数の数だけ空白文字を印刷する関数である。
pt minus 1 pt
void bhnode::dump(int indent)
{
int i;
spc(indent); cerr << "node center pos " << cpos ;
cerr << endl;
spc(indent); cerr << "node cm " << pos << " m " << mass ;
if (nparticle == 1){
cerr << " IS LEAF" ;PRL(nparticle);
particle * p = pfirst;
for(i = 0; i < nparticle; i++){
for(int j=0;j<indent+2;j++)cerr << " ";
PRL(p->pos);
p= p->next;
}
}else{
cerr << " IS _not_ LEAF ";PRL(nparticle);
for(i=0;i<8;i++){
if (child[i] != NULL){
child[i]->dump(indent + 2);
}
}
}
}
以下は N=2 での実行結果である
pt minus 1 pt
BHtree
Enter n:2
nbody = 2, power_index = 0
i = 0, pp[i].pos = -0.20707 0.680971 -0.293328
i = 1, pp[i].pos = -0.106833 -0.362614 0.772857
bn->sanity_check() = 0
node center pos 0 0 0
node cm -0.156952 0.159178 0.239765 m 1 IS _not_ LEAF nparticle = 2
node center pos -0.5 -0.5 0.5
node cm -0.106833 -0.362614 0.772857 m 0.5 IS LEAFnparticle = 1
p->pos = -0.106833 -0.362614 0.772857
node center pos -0.5 0.5 -0.5
node cm -0.20707 0.680971 -0.293328 m 0.5 IS LEAFnparticle = 1
p->pos = -0.20707 0.680971 -0.293328
Direct force
i = 0 p->pos = -0.20707 0.680971 -0.293328 p->phi = -0.33364 p->acc = 0.014891 -0.155032 0.158389
i = 1 p->pos = -0.106833 -0.362614 0.772857 p->phi = -0.33364 p->acc = -0.014891 0.155032 -0.158389
Tree force
i = 0 p->pos = -0.20707 0.680971 -0.293328 p->phi = -0.33364 p->acc = 0.014891 -0.155032 0.158389
i = 1 p->pos = -0.106833 -0.362614 0.772857 p->phi = -0.33364 p->acc = -0.014891 0.155032 -0.158389
それらしいツリー構造ができていること、力の計算結果が総当りとツリーを使っ たもので一致していることがわかる。もちろん、粒子が2個の 時にはツリーを使っても相手は1個なので、結果は完全に一致しなければおか しい。まあ、確かに一致しているというわけである。
2粒子でそれらしく動いたからといって正しいということにはならない。もっ と大きな粒子数でいろいろ調べ、さらに opening angle も変えて調べる必要 があるがこれは課題ということにする。