3 各関数のテスト

前回のテストプログラムに機能を追加して、 実際に重力を計算し、ツリー構造を使わないで総当りで計算したものとくらべてみよう。メインプログラムは以下のような感じになる。

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(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->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

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