というわけで、一応木構造をつくるためのデータ構造、アルゴリズムとその実 現を紹介した。これを実際に多体問題の数値解を求めるプログ ラムにするためには、ノードの物理データの計算、重力計算、時間積分、初期 条件の生成/読み込み、結果の表示/解析等いろんなことをするプログラムを 作る必要がある。
それらに入る前に、とりあえず今まで出来たところを動かしてみることにする。
割合いろんな関数を作っていろんなことをやらせるので、全体を動かして時に ちゃんと動いているかどうかをどうやって判断し、さらに動いていない時にど こがおかしいかを見つけ出すのはそんなに簡単ではなくなる。
ある程度大きな、沢山の関数からできたプログラムをテスト・デバッグするた めの基本的な考え方は以下の 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
それらしいツリー構造ができていることがわかる。
来週は、力の計算とか、残りのところをプログラムする。