というわけで、一応木構造をつくるためのデータ構造、アルゴリズムとその実 現を紹介した。これを実際に多体問題の数値解を求めるプログ ラムにするためには、ノードの物理データの計算、重力計算、時間積分、初期 条件の生成／読み込み、結果の表示／解析等いろんなことをするプログラムを 作る必要がある。
それらに入る前に、とりあえず今まで出来たところを動かしてみることにする。
割合いろんな関数を作っていろんなことをやらせるので、全体を動かして時に ちゃんと動いているかどうかをどうやって判断し、さらに動いていない時にど こがおかしいかを見つけ出すのはそんなに簡単ではなくなる。
ある程度大きな、沢山の関数からできたプログラムをテスト・デバッグするた めの基本的な考え方は以下の 2 つである。
もちろん、プログラムをデバッグする時の大きな問題は、
ということである。特に後者がデバッグを難しいものにする。
というわけで、まず、ツリー構造を作り、それが正しいかどうか調べる。メインプログラムは以下のような感じになる。
pt minus 1 pt
#ifdef TREETEST
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(myvector(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->dumptree();
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;
myvector 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(myvector & cpos, // center of the box
myvector & 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::dumptree(int indent)
{
int i;
spc(indent); cerr << "node center pos " << cpos ;
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]->dumptree(indent + 2);
}
}
}
}
以下は 
での実行結果である
pt minus 1 pt
simpletree>BHtreetest
Enter n:4
nbody = 4, power_index = 0
i = 0, pp[i].pos = -0.815405 -0.0255656 0.0535006
i = 1, pp[i].pos = -0.0911332 -0.533643 0.662584
i = 2, pp[i].pos = 0.863463 0.136119 0.112189
i = 3, pp[i].pos = -0.49528 -0.403606 0.751962
bn->sanity_check() = 0
node center pos 0 0 0 IS _not_ LEAF nparticle = 4
node center pos -0.5 -0.5 0.5 IS _not_ LEAF nparticle = 3
node center pos -0.75 -0.25 0.25 IS LEAFnparticle = 1
p->pos = -0.815405 -0.0255656 0.0535006
node center pos -0.25 -0.75 0.75 IS LEAFnparticle = 1
p->pos = -0.0911332 -0.533643 0.662584
node center pos -0.25 -0.25 0.75 IS LEAFnparticle = 1
p->pos = -0.49528 -0.403606 0.751962
node center pos 0.5 0.5 0.5 IS LEAFnparticle = 1
p->pos = 0.863463 0.136119 0.112189
それらしいツリー構造ができていることがわかる。
来週は、力の計算とか、残りのところをプログラムする。