./jisshuu-2009-1.html

ToC

4. N 体シミュレーション実習

ここからが実習本番です。準備は大変だったかもしれませんが、実際に研究に使われているパッケージですので、役に立つこともあると思います。

この実習では、

いくつかの簡単に作成できる系の進化をシミュレーションしてみる。
その時に、グローバルな物理量、例えばビリアル比等の進化を調べる。
計算精度と、その結果への影響について調べる。

Cold Collapse とは、「冷たい」、つまり、粒子の速度分散等が系のサイズを支えるだけの大きさがない粒子分布を作って、そこから時間発展をさせるとどんなことが起こるか、という問題です。もちろん、宇宙でそのままのことが起こるわけではないのですが、現在のビッグバンモデルでは、例えば銀河といった現在ある構造は、一様な宇宙膨張から重力ゆらぎが成長していって、膨張から外れて収縮に移ってできるわけです。大雑把な理解としては、膨張から収縮に移る、ということはどこかで最大膨張になったところがあり、そこではその中の粒子は全部止まっているわけですから、 Cold Collapse というのはそういう状態のモデルであると考えることができます。

もっとも、現在の CDM 宇宙論では、小さいスケールの構造が先にできるので、銀河系のような大きなサイズのものが一気に Cold Collapse でできるわけではありませんが、それでも Cold Collapse で何が起きるか、ということを理解しておくことは宇宙における構造形成を理解する基本といえます。

4.1.2. 初期条件の作成とその検証

能書きはさておき、まず初期条件を作ってみます。

nemo の中に役に立ちそうなプログラムはあるかをみます。nemo のプログラムでモデルを作るものは大体 mk なんとかという名前なので、ls してみると、、、

 %ls $NEMOBIN/mk*
 /home/makino/work/nemo_cvs/bin/mkbaredisk
 /home/makino/work/nemo_cvs/bin/mkconfig
 /home/makino/work/nemo_cvs/bin/mkcube
 /home/makino/work/nemo_cvs/bin/mkdisk
 /home/makino/work/nemo_cvs/bin/mkexpdisk
 /home/makino/work/nemo_cvs/bin/mkexphot
 /home/makino/work/nemo_cvs/bin/mkflowdisk
 /home/makino/work/nemo_cvs/bin/mkhom
 /home/makino/work/nemo_cvs/bin/mkhomsph
 /home/makino/work/nemo_cvs/bin/mkkd95
 /home/makino/work/nemo_cvs/bin/mknemo
 /home/makino/work/nemo_cvs/bin/mkommod
 /home/makino/work/nemo_cvs/bin/mkop73
 /home/makino/work/nemo_cvs/bin/mkorbit
 /home/makino/work/nemo_cvs/bin/mkpdoc
 /home/makino/work/nemo_cvs/bin/mkplummer
 /home/makino/work/nemo_cvs/bin/mkpolytrope
 /home/makino/work/nemo_cvs/bin/mksphere
 /home/makino/work/nemo_cvs/bin/mkspiral
 /home/makino/work/nemo_cvs/bin/mktestdisk

やたら沢山あります。ディスクっぽいものが多いですが、ここでは mkhomsph というのを使ってみます。

 %mkhomsph help=
 mkhomsph out=??? nbody=2048 rmin=0 rmax=1.2 2t/w=0.0 vmax=0 power=0 seed=0 zerocm=t headline= VERSION=1.4b

 %mkhomsph test4k 4096
 %glnemo test4k

という感じで、どんなものができるかみてましょう。

Figure 19: mkhomsph test4k 4096 でできた初期モデルを glnemo で見る。

なんとなくそれらしい丸いものができている、というのはわかると思います。しかし、粒子分布は本当に一様でしょうか？これを確認するためには、中心からの距離によって粒子を並べかえて、累積質量分布を作ってみるのが一つの方法です。もちろん、半径方向の密度プロファイルでもよくて、これは単に累積質量分布を半径で微分したものですが、 N体で作った粒子分布は離散的なので微分は簡単ではありません。

並べかえには snapsort というプログラムが使えます。どんなオプションがあるかは、ヘルプオプションや man ページで確認して下さい。

 %snapsort test4k test4k.sort rank=r

で、半径方向に並べかえてくれます。次に、「半径方向の累積質量分布」を作って、これを半径をグラフにしてみたい、ということになります。

ここでは、 PGPLOT に簡易なユーザーインターフェースをつけた wip というプログラムでグラフを書き、そのためのデータは snapprint の出力を awkで加工してみます。

 %snapprint test4k.sort options=r,m | awk '{macc += $2; print macc, $0}' > tab.out

ruby や perl を使ってもよいですが、こういう使い方なら awk が簡単です。ここでは、ソートした test4k.sort から、粒子の位置の半径と質量が1行になった出力を作り、それを awk で処理します。 awk でやっていることは macc という変数に 1行の2つめの数 ($2) を加算し、加算した macc と読み込んだ行全体を出力する (print macc, $0) というものです。 {} でくくった中を1行毎に実行する、というのが awk の基本動作で、さらに '' でくくっているのはこれ全体をシェルに解釈させないで awk に渡すためです。 tab.out がどんなファイルかみてみましょう。

 % head tab.out
 0.000244141 0.0480947 0.000244141 
 0.000488282 0.0717042 0.000244141 
 0.000732423 0.120009 0.000244141 
 0.000976564 0.137944 0.000244141 
 0.00122071 0.154036 0.000244141 
 0.00146485 0.157324 0.000244141 
 0.00170899 0.163142 0.000244141 
 0.00195313 0.163527 0.000244141 
 0.00219727 0.164162 0.000244141 
 0.00244141 0.164364 0.000244141

こんな感じ(2列目の距離の数字は乱数の初期値によって違うはず) で、最初のほうでは距離の数字が小さく

 % tail tab.out
 0.997804 1.21369 0.000244141 
 0.998048 1.2137 0.000244141 
 0.998293 1.21397 0.000244141 
 0.998537 1.2148 0.000244141 
 0.998781 1.21576 0.000244141 
 0.999025 1.21585 0.000244141 
 0.999269 1.21661 0.000244141 
 0.999513 1.21795 0.000244141 
 0.999757 1.21842 0.000244141 
 1 1.2211 0.000244141

最後のほうでは大きくなっていればもっともらしいです。グラフは

 %wip
 Could not open user's WIP initialization file.
 Setting up default device [/XWINDOW]
 Welcome to WIP Version: 2.3 22jan98

 WIP> device /xs
 WIP> data tab.out
 WIP> xcol 2
 WIP> ycol 1
 WIP> limit
 WIP> box
 WIP> conn

こんな感じでコマンドをいれると図 20 のようなグラフがでるはずです。なお、 wip は end で終了します。

Figure 20: mkhomsph test4k 4096 でできた初期モデルの半径方向質量累積分布

これは、誤差の範囲で、

と一致しているはずです。折角なので、これも計算してみましょう。

 %snapprint test4k.sort options=r,m | awk '{macc += $2; rs=$1/1.2; print macc, $0, rs*rs*rs}' > tab.out

 wip
 Could not open user's WIP initialization file.
 Setting up default device [/XWINDOW]
 Welcome to WIP Version: 2.3 22jan98

 WIP> device /xw
 WIP> data tab.out
 WIP> xcol 2
 WIP> ycol 1
 WIP> era
 WIP> limit
 WIP> box
 WIP> conn
 WIP> ycol 4
 WIP> color 2
 WIP> conn
 WIP> end

が4列目のデータになっているはずなので、それを ycol 4 で y軸方向のデータとして読み、 color 2 で赤に変えてプロットしています。グラフは載せませんが、ちゃんとあっていることを各自確認して下さい。

人が作ったプログラムだから、自分が作ったものよりは信用できますが、だからといってチェックしないで使ってはいけません。これは、プログラムが間違っている可能性は(それまでに沢山の人が使っているものでも)必ずあるし、また、プログラムが正しくても自分の使いかたが間違っていることもある、というかこちらは良くある、からです。

4.1.3. 時間積分

さて、時間積分です。今回は、 hackcode1_qp という、 nemo に初めからはいっているプログラムを使ってみることにします。これは、 Barnes-Hut ツリーアルゴリズムで有名な Barnes 本人が1985年頃に書いたものです。この、nemo というパッケージ自体が、その頃にプリンストン高等研究所にいた Piet Hut, Josh Barnes と、現在も開発、メンテナンスを続けている Peter Teuben の3人が共同で始めたもので、構造化されたバイナリファイルや、引数で数式を与えるとコンパイルして自分自身にリンクする、といった仕掛けの最初の実装は Barnes によるものです。

と、これは余談ですが、 gadget とか、もっと速いコードを使ってもよいですが、今日の実習ではそこまではしないで、ユーザーインターフェースが簡単で安全に使える hackcode1_qp を使うことにします。例によって help をみてみます。

 %hackcode1_qp help=
 hackcode1_qp in= out= restart= continue= save= nbody=128 seed=123 cencon=false freq=32.0 eps=0.05 tol=1.0 fcells=1.0 options=mass,phase tstop=2.0 freqout=4.0 minor_freqout=32.0 debug=false VERSION=1.4

色々謎なパラメータがあるので、 man hackcode1_qp 、、、あれ、そんなの知らない、とでますね。 man hackcode1 でマニュアルを見ることができます。 qp がつくと4重極モーメントを使って、より正確に力を計算するものになっています。

試しに、入力、出力だけを指定して、あとはデフォルトでやってみましょう。

 %hackcode1_qp test4k test4k.outsnap

4.1.4. エネルギー誤差のチェック

画面に一杯出力がでますが、なんだかわからないかもしれません。最後の数行は

        tnow       T+U       T/U     nttot     nbavg     ncavg   cputime
       2.000   -0.4579   -0.7211    796020        39       155      0.13

                cm pos    0.0007   -0.0011   -0.0015
                cm vel    0.0026   -0.0014   -0.0037

        particle data written

という感じのはずです。(細かい数字は違います) tnow は時刻、 T+U は系のトータルエネルギー、 T/U は運動エネルギー T とポテンシャルエネルギー U の比、いわゆるビリアル比です。力学平衡の状態ならこれは -0.5 です。その次の3つの数字はツリーコードでの計算量を示していて、 1ステップで粒子への力を計算した回数、1粒子への、粒子からの力の数の平均、1粒子への、ツリーノードからの力の数の平均、となっています。その下の cm pos, cm vel は、系全体の重心の位置と速度です。大体0ですが少し動いているのがわかります。

ウインドウをスクロールして上のほうに戻ると、最初のほうの出力は

 Hack code

        nbody        freq         eps         tol
         4096       32.00      0.0500      1.0000

         options: mass,phase

         tnow       T+U       T/U     nttot     nbavg     ncavg   cputime
        0.000   -0.4983   -0.0000    425175        17        86      0.00

                 cm pos   -0.0000   -0.0000    0.0000
                 cm vel    0.0000    0.0000    0.0000

         particle data written

         tnow       T+U       T/U     nttot     nbavg     ncavg   cputime
        0.031   -0.4983   -0.0003    424385        17        86      0.01

                 cm pos   -0.0000    0.0000   -0.0000
                 cm vel   -0.0000    0.0000   -0.0000

というような感じだと思います。freq はタイムステップの逆数、 eps は重力ソフトニングの値、 tol はツリーコードの opening angle で、小さくすると精度は上がりますが計算量も増えます。ここでエネルギーである T+U の値を比べると、 -0.4983 だったのが -0.4579 になって、大きく変わってしまっています。これでは駄目そうな感じがするので、タイムステップを小さくしてみましょう。

 %hackcode1_qp test4k test4k.outsnap freq=128

 Hack code

        nbody        freq         eps         tol
         4096      128.00      0.0500      1.0000

         options: mass,phase
 ### Fatal error [hackcode1_qp]: stropen: file "test4k.outsnap" already exists

nemo のプログラムは、出力ファイルがすでにあったらアボートするようになっています。

 %rm test4k.outsnap ; hackcode1_qp test4k test4k.outsnap freq=128

        tnow       T+U       T/U     nttot     nbavg     ncavg   cputime
       2.000   -0.5006   -0.7081    780717        38       151      0.92

                cm pos    0.0012   -0.0003   -0.0008
                cm vel    0.0038    0.0004   -0.0020

freq=128 を指定して、タイムステップをデフォルトの 1/4 にしてみた結果です。劇的にエネルギー保存が改善されているのがわかります。最初のでは 10% 程度だったものが、 0.5% 程度になっています。

理論的には、hackcode1 で使っている leap frog 公式では、エネルギーの誤差は時間ステップの2乗に比例します。なので、これは大体正しい振る舞いといえます。

4.1.5. アニメーション

結果を snapplot や glnemo で見てみます。

 %snapplot  test4k.outsnap
 %glnemo  test4k.outsnap

どちらも、パラパラ漫画ふうですが、出力間隔が大き過ぎて気分がでません。間隔を変えるオプションは freqout なので、これも指定してもう一度実行してみます。

 %rm test4k.outsnap ; hackcode1_qp test4k test4k.outsnap freq=128 freqout=32

今度は snapplot 等でもう少し気分がでると思います。

もうちょっとがんばって、ムービーを作ってみましょう。glnemo で、アニメーションのアイコン(左のメニューの一番下)を押すと、アニメーションのメニューがでます。ここでまず、 Recording/Playing の「O」を押すとアニメーションの実行記録が始まります。そこで、左側の play ボタンを押し、表示が終わる数フレーム前に(ここが極めて重要) Recording/Playing の「X」を押します。

それが済んだら、入力スナップショットを巻き戻し(play の3角ボタンの下の巻き戻しっぽいアイコン)、アニメーションのメニューウインドウのほうの play の「>」ボタンを押して、

アニメーションが動くこと
最後に、Snapshot ... end of snapshot reached! の警告ウインドウがでないこと

を確認して下さい。警告ウインドウがでてしまったら、 Recording からやり直して下さい。

警告ウインドウがでたもので、以下の Rendering を実行すると、何故か X が固まります。大変危険ですのでここは十分に注意して下さい。

上手くいったら、 Rendaring のメニューに移動して、スナップショットを巻き戻した上で、「Start」をクリックします。そうすると、どこかに (デフォルトは /tmp/render/frame.xxxxx.png, xxxxx は連番)フレーム毎の png ファイルができます。これを、例えば animation gif に変換するなら、ファイルができているディレクトリで

 %convert *.png anim.gif

でOKです。gif アニメは、ブラウザで見ることもできますが、何が起こっているか、の物理を直観的に理解するためには、1コマずつ進めるとか、逆に戻るとかが簡単にできるツールが必須です。結構、なかなか適当なものがありません。 gif アニメの表示には xanim が便利ですが、開発が止まって 10年ほどたつこともあり Cygwin や Linux の最近のディストリビューションにははいっていません。 Gnome があれば、GImageView はあまり便利ではないですが一応それらしいことができます。Windows でも Explorer で沢山の画像ファイルを見る、というのが簡単な方法のような気がします。

4.1.6. エネルギー誤差の時間変化グラフを作る

さて、さっきは、エネルギーの最終の値をチェックしましたが、これだけでは普通は十分ではありません。最後の値は良くても、途中でデタラメな値になっていることもあるからです。これをチェックするためには、時間変化のグラフを作る必要があります。

hackcode1_qp は、出力のスナップショットファイルの中に、エネルギーや他の色々な情報を埋め込んでいます。これは

 %tsf  test4k.outsnap

で見ることができます。最後のほうの出力が

 set SnapShot
   set Parameters
     int Nobj 4096 
     double Time 2.00000 
   tes
   set Particles
     int CoordSystem 66306 
     double Mass[4096] 0.000244141 0.000244141 0.000244141 0.000244141 
       0.000244141 0.000244141 0.000244141 0.000244141 0.000244141 
       0.000244141 0.000244141 0.000244141 0.000244141 0.000244141 
       0.000244141 0.000244141 0.000244141 0.000244141 0.000244141 
       . . .
     double PhaseSpace[4096][2][3] 0.153245 0.0136813 -0.205303 1.79778 
       -0.760285 -0.378773 -0.454063 0.919896 -0.552425 -0.704771 
       1.34618 -1.00689 -0.983543 -0.662581 -0.412249 -1.66848 -1.23127 
       -0.691078 0.153962 0.0313744 -0.0690967 1.43265 -0.379112 
       . . .
   tes
   set Diagnostics
     double Energy[3] -0.500630 1.21438 -1.71501 
     double KETensor[3][3] 0.377060 0.0249498 -0.0361606 0.0249498 
       0.380638 -0.0253557 -0.0361606 -0.0253557 0.456682 
     double PETensor[3][3] -0.413097 0.0255456 0.0248519 0.0249459 
       -0.417775 0.0157845 0.0236208 0.0148032 -0.599691 
     double AMTensor[3][3] 0.00000 -0.000262578 0.000150161 0.000262578 
       0.00000 -0.000102650 -0.000150161 0.000102650 0.00000 
     double CMPhaseSpace[2][3] 0.00120219 -0.000280852 -0.000787409 
       0.00384288 0.000367961 -0.00195136 
     double cputime 0.778333 
   tes
 tes

ですから、

Time がある行の数字
Energy がある行の最初の数字

をプロットすればよさそうです。また awk を使うことにすると

 %tsf test4k.outsnap | awk '/Time/{t=$3}/Energy/{print t, $3}' > tab.out

です。 /Time/{t=$3}は、入力行が正規表現 /Time/ と一致したら、括弧内を実行、で、 Time という文字列w含む行の3個めの文字列(区切りはスペース)を t という変数に代入、/Energy/{print t, $3} は同様に Energy を含む行にきたら、 t の値と 3個めの文字列をプリント、となります。

wip でグラフを書いた例が図 21です。結構駄目ですね。

Figure 21: エネルギーの時間変化

T=1.5 くらいでずれが非常に大きくなっています。これは、系全体が非常に小さく収縮し、速度も大きくなるためですが、10% はいかにも大きいので、この場合にはタイムステップをもっと小さくしたほうが安心できます。

4.1.7. 初期条件を変えてみる

エネルギー誤差が大きいのは、タイムステップが大き過ぎるせいではありますが、初期条件が「一様」という特別なものであるためでもあります。

密度一様な球のコラプスは、密度一様なまま進む、という特徴があります。実際の計算では粒子数が有限なので1点にはあつまらないですが、原理的には1点に集まってしまうわけで、そうするとどんなにタイムステップを短くしても上手くいきません。そういう、破綻する場合の積分方法の話も講義では少し触れたかもしれませんが、実習ではそこまで高度なことはしません。

逆にいうと、密度が一様でなければ、そこまでひどいことは起きないわけです。密度一様でなくて、半径の-1乗にしてみます。

 %mkhomsph test4kb 4096 power=1

これも、累積質量分布がちゃんと出来ていることを確認しましょう。

 % hackcode1_qp test4kb test4kb.outsnap freq=128 freqout=32 tstop=10

で少し長い計算をして glnemo でアニメーションを作ったりしてみましょう。

少しありえないほど細長いですが、楕円銀河のように見えないこともないものが出来ているはずです。エネルギーエラーも、ずっと小さくなっていることが確認できます。

エネルギーエラーを計算する awk のスクリプトを少し変えると、ビリアル比 T/U が計算できます。この場合にビリアル比のグラフを出す、というのが最初の課題、ということにします。 Energy の行で、 4 個めが T, 5個めが U なので、 $3 の代わりに $4/$5 をプリントさせればよいわけです。

さて、タイムステップをどんどん小さくしていくと、エネルギーエラーは小さくなりますが、ある程度以下にはならないことがわかります。それはなぜか、また、さらに小さくするにはどうすればよいか、を調べるのは発展課題とします。

4.1.8. 密度分布を書く

前に、初期条件をチェックするには累積質量分布を使いましたが、「どのような銀河ができたのか」といった物理的な性質を見るのには累積分布はあまり便利ではありません。なので、半径方向の空間密度分布を書いてみます。

nemo には radprof というプログラムがあるのですが、動作が気に喰わないというか今一つ使えないので、snapprint の出力から書いてみることにします。前に書いたように、どうしてもN体では粒子数によるゆらぎがあります。ここでは粒子 128 個毎に、その2つの粒子にはさまれた球殼の平均密度を、2つの粒子の半径の平均の半径での密度とします。

式で書いたほうがわかりやすいですね。粒子が半径順に並んでいて、半径がで質量がとすると、粒子にはさまれた範囲の密度を

とします。これを計算するプログラムですが、C でも Fortran でも好きな言語でよいですが、以下は awk の例です。

 {
     m+= $2
     if ( NR % 128 == 0){
        r1 = $1
        print (r1+r0)/2, m/(r1*r1*r1-r0*r0*r0)*0.2387
        m=0
        r0=r1
     }
 }

awk には、

未定義の変数は値が 0
変数の値は文字列または実数、どっちになるかは適当
暗黙に、1行づつ読んで処理、というループになっている
1行の、複数のフィールドの数字があらかじめ $1, $2, ... という名前で使える

といった特性があるので、割合簡潔に処理を表現できているのがわかります。なお、 awk なんていう大昔の言語は嫌だ、もうちょっと現代的なものがよい、という向きには、例えば Ruby を使うなら

 r0=m=0
 while s=gets()
   a=s.chomp.split
     m+= a[1].to_f
     if ( $. % 128 == 0)
        r1 = a[0].to_f
        print (r1+r0)/2, " ",m/(r1*r1*r1-r0*r0*r0)*0.2387, "\n"
        m=0
        r0=r1
     end
 end

という感じですが、これでは awk に比べて面倒なだけでメリットはありません。どちらの場合でも

  %snapprint test4k.sort options=r,m | awk -f density.awk > tab.out

とか

  %snapprint test4k.sort options=r,m | awk -f density.rb > tab.out

とかすれば、それらしいデータがでているはずです。(もちろん、上のコードをそういう名前でファイルにセーブする必要があります) 前と同様にグラフを書いてみましょう。折角 ruby を使うのでもう少し気の効いたことをするなら、プログラムを2行ばかり変更して

 r0=m=0
 open("|snapprint #{ARGV[0]} options=r,m"){|io|
   while s= io.gets()
     a=s.chomp.split
     m+= a[1].to_f
     if ( $. % 128 == 0)
       r1 = a[0].to_f
       print (r1+r0)/2, " ",m/(r1*r1*r1-r0*r0*r0)*0.2387, "\n"
       m=0
       r0=r1
     end
   end
 }

として見る、といったこともできます。この場合は、 ruby プログラムの中で snapprint #{ARGV[0]} options=r,m を実行して、その出力を io.gets() で読む、というふうにしています。こうしておくと、

  %ruby density2.rb test4k.sort

という形で実行できるので、沢山のファイルについて密度分布を書く、といった場合にやりやすいし、後になってこのプログラムの使いかたをすっかり忘れた時にも、 snapprint をどう使ったかは書いてあるので安全です。まあ、ヘルプとか、色々つけておくとあとで助かるのですが、世間で普通に使われている方法は面倒であまり使いやすくないように思います。プリンストン高等研究所の Piet Hut (Barnes-Hut ツリーコードの Hut) と牧野が少し前に「コマンドラインオプションのパーザはどうあるべきか」だけで本一冊分くらいの原稿を書いたことがあるので、そういうのに興味がある人は読んで見ると得ることがあるかもしれません。

グラフは、この場合には

 %wip
 Could not open user's WIP initialization file.
 Setting up default device [/XWINDOW]
 Welcome to WIP Version: 2.3 22jan98
 
 WIP> device /xw
 WIP> data tab.out
 WIP> xcol 1
 WIP> ycol 2
 WIP> limit
 WIP> era
 WIP> box
 WIP> conn

という感じで書けると思います。これを、シミュレーションの最後のスナップショットに適用してみましょう。最後のスナップショットに snapsort とかを適用するには、最後のスナップショットを切り出すコマンドを使うのが便利です。

 %snaptrim test4kb.outsnap - times=10 | snapsort - - rank=r |\
 ? snapprint - options=r,m | ruby density.rb > tab.out
 ### nemo Debug Info: r m 
 ### nemo Debug Info: time = 10  npart =   1     ndiag =   1     outputing particles

4.1.9. 提出課題

test4k.outsnap と test4kb.outsnap のそれぞれについて、最終状態での密度分布を計算し、両対数グラフにして提出せよ。
test4kb.outsnap のグラフはスナップショットの絵と比べて「もっともらしい」、つまり、構造の特徴を捉えていると考えられるか？もしも捉えていないならそれはなぜで、どう処理を変更すればよいか検討し、「もっともらしい」ものになった結果を提出せよ。
test4k からの数値積分について、時間刻みを2桁以上の範囲で変化させて、 t=0 から t=2 まで積分した時の最大誤差と最終誤差を時間刻みの関数としてプロットしたグラフを作り、提出せよ。

なお、提出用のファイルを作るには、画面のスクリーンショットとかでもいけなくはないですが、例えば TeX の場合 ps ファイルを作るのが便利です。 WIP の場合、

 %wip
 Could not open user's WIP initialization file.
 Setting up default device [/XWINDOW]
 Welcome to WIP Version: 2.3 22jan98

 WIP> device xxx.ps/vcps

と、device コマンドで出力を ps ファイル(この場合、名前が xxx.ps) にすることができます。この後で

 WIP> viewport 0.2 0.9 0.3 0.8

を実行すると大体正方形のグラフになるはずです。

4.2. 円盤銀河

実際の実習ではこの辺で既に力つきるかもしれませんが、折角なのでもう少し気の効いた銀河モデルを作ってみましょう。nemo には Kuijken と Dubinski の作った、 GalactICS という銀河モデル作成パッケージがはいっています。

 %mkkd95 testkd

で何かができます。 glnemo でぐるぐる回してみると、ハロー、ディスク、バルジからできているらしいことがわかります。

4.2.1. 単位系の修正

N体計算をする時の割合大事な話として、どういう単位系で考えるか、ということがあります。例えば、銀河系の計算をするとして、 kpc, Myr, solar mass といった単位系を使っていたのでは、物理的な理解が必ずしも簡単ではありません。

従って、多くの場合に、いわゆる standard unit あるいは N-body unit といわれる単位系を使います。これは、系の質量が1、重力定数も1、系の全エネルギーが -1/4 であるような単位系です。

何故 -1/4 で -1 でないのか、というと、この場合ポテンシャルエネルギーが -1/2 であり、全粒子ペアについての平均の逆数、つまり粒子間距離の調和平均をとると、それが1になっている、ということで -1/4 を選んでいます。 nemo の多くのプログラムはこれに従っているか、大体従っています。例えば mkplummer の場合にはそうなっています。

しかし、GalactICS は作った人にそういう考え方はないので、全く違う単位系 (但し、重力定数は1)ででてきます。なので、ここではまず、GalactICS の出力を standard unit に変換してみます。

元の系の質量が、ポテンシャルエネルギーが、運動エネルギーがであったとして、質量を倍、位置座標を倍、速度を倍すれば、新しい質量とエネルギーは

となります。この時に、ビリアル比 T/V を変化させないように速度をスケールするには、

となるわけです。というわけで、まずできたモデルの質量とエネルギーを計算してみます。

 % hackforce_qp testkd - tol=0.5 eps=0.025 |\
 ? snapprint - options=m,"m*(phi+vx*vx+vy*vy+vz*vz)*0.5" |\
 ? awk '{mtot += $1; etot += $2}END{ print mtot, etot}'
 ### nemo Debug Info: m m*(phi+vx*vx+vy*vy+vz*vz)*0.5 
 ### nemo Debug Info: initial rsize: 4.000000    rmin: -2.000000  -2.000000  -2.000000
 ### nemo Debug Info:   final rsize: 128.000000    rmin: -38.000000  -54.000000  -86.000000
 6.21204 -2.48885

ちょっと意味不明ですね。 "|" はパイプで、左のコマンドの出力を右のコマンドの入力にします。 nemo のプログラムは、ファイル名に "-" を指定すると出力なら標準出力、入力なら標準入力になるので、上の例のようにパイプでつなぐことができます。このために、中間ファイルとかを作る必要が(少なくとも見かけ上は、、、)なく、複雑な処理をパイプで実現できます。

この例では、最初の hackforce_qp は、ツリーコードでポテンシャルと加速度を計算して、それをスナップショットに付け加えて出力するものです。次の snapprint では、各粒子の質量と、「運動エネルギーとポテンシャルエネルギーの半分の合計」を書かせています。ポテンシャルエネルギーを半分にするのは、ペアを2重に数えてしまっているのを補正するためです。最後の awk では、質量、エネルギーをそれぞれ合計していって、最後に (END {} のところで)出力しています。

さて、これで質量とエネルギーがわかったので、あとは電卓でスケールを計算して、

 %snapscale testkd testkd.scaled mscale=0.1609777 \
 ?   rscale=0.25798248826425 vscale=0.789928431392108

とすれば一件落着なのですが、これだと、GalactICS で銀河モデルを作る度に電卓で面倒な計算をして、数字を人間が打ち込む(まあ、 xcalc とか使っていればコピペもできますが)ことになります。式に従って計算する、というのは計算機のほうが人間よりも得意に決まっているので、これも計算機にやらせてみましょう。

 #!/bin/csh -f
 #
 # scaletoheggie.csh
 #
 # scale a snapshot to the standarr (aka Heggie) units
 #
 # Usage: scaletoheggie.csh snapshot_file
 # output file: snapshot_file.scaled
 snapscale $1 $1.scaled \
 `hackforce_qp $1 - tol=0.5 eps=0.025 |\
   snapprint - options=m,"m*(phi+vx*vx+vy*vy+vz*vz)*0.5" |\
   awk '{mtot += $1; etot += $2}END{a=1/mtot;b=-a*a*etot/0.25; print "ms=" a, "rs=" b, "vs=" sqrt(a/b)}'`

こんな感じのスクリプトを作って、実行すると

 %csh -f scaletoheggie.csh testkd
 ### nemo Debug Info: initial rsize: 4.000000    rmin: -2.000000  -2.000000  -2.000000
 ### nemo Debug Info: m m*(phi+vx*vx+vy*vy+vz*vz)*0.5 
 ### nemo Debug Info:   final rsize: 128.000000    rmin: -38.000000  -54.000000  -86.000000
 ### Warning [snapscale]: Resolving partially matched keyword ms= into mscale=
 ### Warning [snapscale]: Resolving partially matched keyword rs= into rscale=
 ### Warning [snapscale]: Resolving partially matched keyword vs= into vscale=

 %tsf testkd.scaled
 char Headline[12] "-1250698778"
 char History[80] "snapscale testkd testkd.scaled ms=0.160978 rs=0.25798
   3 vs=0.789928 VERSION=3.2b"
 char History[82] "tabtos galaxy ../testkd nbody,time mass,pos,vel headl
   ine=-1250698778 VERSION=1.5a"
 set SnapShot
   set Parameters
     int Nobj 18000 
     double Time 0.00000 
   tes
   set Particles
     int CoordSystem 66306 
     double Mass[18000] 1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 
       1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 
       1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 
       1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 1.75180e-05 
       . . .
     double PhaseSpace[18000][2][3] -0.0729338 0.257252 0.00723250 
       -0.607239 -0.289207 0.212093 -0.261455 -0.103321 0.0670730 
       0.159549 -0.373926 0.129851 -0.849309 -0.927010 -0.00670115 
       0.560895 -0.433799 0.0387852 0.0114802 0.336401 -0.0339088 
       . . .
   tes
 tes

で、 tsf の出力から、 snapscale に渡っているオプションはそれらしいことがわかります。本当にエネルギー等が正しいかどうかはもう一度全エネルギーと質量を計算して確認しましょう。上のスクリプトでは `command`、つまり、バッククォートで囲ったコマンドは、そのコマンドの実行結果に展開される、という shell の割合便利な機能を使って、awk の出力を snapscale のオプションの格好で出したものを直接 snapscale に渡しています。 ms= とか書くと、 nemo のコマンドライン解釈ルーチンが、一致するオプションを見つけてきます。この場合は、 in out mscale rscale vscale の順でオプションを与えているので ms= とかは書かなくてもよいのですが、ここでは書いてみました。

では、適当な初期条件で円盤銀河2つをぶつけてみましょう。

 %snaprotate testkd.scaled testkd.rotated theta=45 spinvector=1,0,0
 %snapstack  testkd.rotated testkd.rotated  merger.in deltar=5,0,0\
 ? deltav=-0.7,0.35,0 zerocm=true
 %hackcode1_qp merger.in merger.out freq=128 eps=0.02 tstop=20\
 ? freqout=16 minor_freqout=16 &

これはとても長い時間がかかるので、時々 glnemo で様子をみながら休憩、というところで、今日の実習はおしまいにします。

 %glnemo merger.out select=0:7999,8000:17999,18000:25999,26000:35999 blending=t

とすると、ディスク、ハローで色を変えることができます。また、オプションメニューからさらに色々変更することもできます。

ここで作った銀河モデルは粒子数 18000 で、円盤銀河の色々な性質を調べるには到底十分とはいえません。従って、実際の研究でシミュレーションをするには、 hackcode1_qp よりももっと速くシミュレーションできる必要があります。

これには、

プログラムを改良して速くする
速い計算機を使う。使えるようにプログラムを改良する
速い計算機を作る

といった、様々なアプローチが行われています。今日やった計算は 36000 粒子ですが、現在の宇宙論研究の最先端では100億を超える粒子数を使った計算も行われています。

ToC