/nbody-with-center

Simple planetary ring code with FDPS

Primary LanguageC++MIT LicenseMIT

User's Guide for nbody-with-center

                   2019/11/26 牧野
                   2021/09/10 アップデート
                   2021/12/05 アップデート

はじめに

この文章では、 nbody-with-center の最低限の使い方について述べる。FDPS を使った惑星リング系計算のサンプルコードである。あんまりテストしてないのでちゃんとあってるかどうか知らないが、エネルギー保存は(ダンピング入ってるので厳密ではないがそこそこ)する。

コンパイル方法等

基本的に FDPS の sample/c++/nbody の下にある nbody.cpp と同じだが、

  • PS_PATH を適切に (fdps の src があるディレクトリをさすように)設定す る

  • PHANTOM_GRAPE はサポートされない

というあたりが違いである。PHANTOM_GRAPE がなくて遅いので、誰か改良して下さいな

 make nbody-with-center :            OpenMP サポート、重心近似
 make nbody-with-center-quad :       OpenMP サポート、4重極近似
 make nbody-with-center-mpi :   OpenMP+MPI サポート、重心近似
 make nbody-with-center-quad-mpi :   OpenMP+MPI サポート、4重極近似

の4種類の実行ファイルをつくれる。MPI で高精度にするなら最後のを使うこと。

実行

実行時オプションについて説明する

-b: 詳細な実行時間プロファイルを生成する。

-D: スナップショット出力の時間間隔である。デフォルトは 1

-d: ログ出力の時間間隔である。デフォルトは 1/8

-E: FDPS の ExchangeLETMode を指定した値にする。デフォルトは 2(EXCHANGE_LET_P2P_FAST)

-e: 衝突時のダンピング係数 eta である。デフォルトは 0。 K オプションを0にしなければそちらから設定される。

-g: 粒子分布の内側カットオフ半径。この半径よりも内側にはいった粒子は取 り除かれる。除かれた粒子は remove.dat に書かれる。MPI の時には removennnn.dat とランク毎のファイルになる

-I: 初期分布を内部生成する。角度方向の粒子数 nx, 半径方向 ny, 半径方 向の幅 dr, 粒子質量 m をコンマで区切り、スペースはいれないで与える。 半径は内側が1、外側が 1+dr となる。例: -I 1000,20,0.1,5e-10 (デフォルト値)

-i: 入力スナップショットの名前。デフォルトはなし。この時は -I オプションのデフォルト値でディスクが生成される。入力ファイルの形式は

 time
 n
 id1 mass1 x1 y1 z1 vx1 vy1 vz1
 ...
 idn massn xn yn zn vxn vyn vzn

である。サンプルファイル ringin を見ること。

-K: 粒子同士が衝突した時ににはねかえるまでの時間(振動周期の半分)をタ イムステップで表した値。デフォルトは 32 である。この値が正である と -k, -e の値に関わらずこれと S オプションから係数が計算される。

-k: 衝突時のバネ定数 kappa である。デフォルトは0。 K オプションを0にしなければそちらから設定される。

-l: ツリーの leaf limit である。デフォルトは 8

-M: 中心星の質量である。デフォルトは 1

-n: ツリーの n_group_limit である。デフォルトは 64

-o: 出力ファイルを置くディレクトリである。デフォルトは ./results である。 出力はスナップショット毎に別のファイルになる。形式は入力ファイル と同じ。また、 t-de.dat にエネルギー、remove.dat に取り除いた粒子が記録される

-p: 重力ソフトニングである。デフォルトは 0

-R: MPI 並列の時の半径方向の領域分割の数。デフォルトは1である。リング の幅とプロセス数を考慮して適切な値を設定すること。

-r: 粒子の物理半径である。デフォルトは0

-S: 反発係数。デフォルト 0.5。垂直方向の反発係数である。接線方向は0。

-s: 時間刻みである。デフォルトは 1.0//1024

-t: ツリーの見込み角 theta である。デフォルトは0.5

-T: シミュレーション終了時刻である。デフォルトは10

-h: ヘルプメッセージを出力して終了する。

初期条件生成

単に動いているかどうかの確認程度だが

ring.rb

を実行すると、標準出力に、半径 1 から 1.1 の間に、半径方向に20個、角度 方向に 1000 個の粒子を、中心質量が1としてケプラー回転の速度を与えた粒 子を生成する。質量は 5e-8 である。引数に4個パラメータを m n dr massの 順番で与えると (m, n は半径方向及び角度方向) これらのパラメータを変更 できる。

実行例

nbody-with-center-quad -i ringin_hot -T 0.5 -r 0.004 -K 128

エネルギーの正確な値はコンパイラやハードウェアによるが、出力の最後のほうが

time:  0.0000000 etot,k,p,l=-4.709488e-04 4.770623e-04 -9.480111e-04 0.000000e+00 energy error: -0.000000e+00
time:  0.1250000 etot,k,p,l=-4.709500e-04 4.806427e-04 -9.540743e-04 2.481589e-06 energy error: +2.479234e-06
time:  0.2500000 etot,k,p,l=-4.709494e-04 4.806416e-04 -9.541657e-04 2.574786e-06 energy error: +1.198108e-06
time:  0.3750000 etot,k,p,l=-4.709492e-04 4.806520e-04 -9.542009e-04 2.599748e-06 energy error: +7.774599e-07

で、エネルギー誤差が 1e-5 以下なら問題なく動作している。エネルギー保存の精度はあまり高くないが、これは粒子間の非弾性衝突によるエネルギー散逸の計算精度が低いためで、大きな問題ではない。

実行例その2

mpirun -np 16 nbody-with-center-quad-mpi -T 1000 -r 0.01 -K 128 -i moontest20k.in -o result -D 0.125 -g 0.95 -S 0.1

時刻500くらいまで回すと月っぽいものができてくる。

アニメーション(86MB).

更新履歴

2021/12/5

  • 粒子の非弾性衝突によるエネルギー散逸を考慮したエネルギー保存チェックにした
  • 主天体にぶつかった粒子を取り除くオプションを追加した
  • オプションの説明、実行例等を追加した。

2021/9/10

  • FDPS 7.0 に対応して、円筒座標での計算をサポートした。
  • オプションを増やし、デフォルトで反発係数が 0.5 になるようにした。
  • 富岳用の Makefile Makefile.a64fx を追加した。