next up previous
Next: 4 時間積分 Up: 計算天文学 II 第10回 データ構造とアルゴリズム(2) Previous: 2 力の計算

3 各関数のテスト

というわけで、一応木構造をつくって重力計算をするためのデータ構造、アル ゴリズムとその実現をきわめて簡単に紹介した。これを実際に多体問題の数値 解を求めるプログラムにするためには、時間積分、初期条件の生成/読み込み、 結果の表示/解析等いろんなことをするプログラムを作る必要がある。

それらに入る前に、とりあえず今まで出来たところを動かしてみることにする。

割合いろんな関数を作っていろんなことをやらせるので、全体を動かして時に ちゃんと動いているかどうかをどうやって判断し、さらに動いていない時にど こがおかしいかを見つけ出すのはそんなに簡単ではなくなる。

ある程度大きな、沢山の関数からできたプログラムをテスト・デバッグするた めの基本的な考え方は以下の 2 つである。

  1. それぞれの関数を単体でテストする。
  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;
}
ここでは、以下のようなチェックを行なっている。

  1. ノードに粒子が1つだけなら、それがちゃんと中にあるかどうか。
  2. そうでなければ、自分の子供は正しい位置、大きさを持つかどうか。
で、さらに再帰的に子供についても同じチェックをする。

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 も変えて調べる必要 があるがこれは課題ということにする。



Jun Makino
平成15年1月27日