というものに変えるのが普通である。これをどのような場合にやってよいか、ま
た、 にどの程度の大きさを取るべきかはここでは議論しないが、
ブラックホールの影響を考える場合には、これらについてはこれらの回りの星
の軌道が大きく変わってしまうような値にしてはいけない。
ちなみに、 この のことを通常ソフトニング長という。面
倒なので単にソフトニングということも多い。
これは小さいほど計算は正しいわけだが、何も考えないで小さくすると計算
が大変なので、計算がおかしくならない範囲で大きくしたい。 SMBH
と IMBH の質量は非常に大きく違うので、このことは、以下の4種類の粒子の
組合せでソフトニングを変える必要があるということを意味する。
-
普通の星-普通の星
-
SMBH-普通の星
-
IMBH-普通の星
-
SMBH-IMBH
重力ポテンシャルは作用・反作用の法則を満たす必要があるから、 普通
の星から例えば SMBH への重力を計算する時と、その逆を計算する時のソフト
ニングは同じでないといけない。つまり、上の4つの場合わけが発生する。
というようなことを考えて某君(以下A君)が書いたのが上のプログラムである。
これは get_force_host、 get_force_host_bh と get_force_initial の3つの
関数がある、重力計算だけのファイルの中身を示している。
プログラムを見ると、なんとなく繰り返しが多いなあという気がする。
r2 = eps2;
xdotv = 0.0;
for (k = 0; k < 3; k++) {
dx[k] = x1[j][k] - xi[ii][k] ;
dv[k] = v1[j][k] - vi[ii][k] ;
r2 += dx[k] * dx[k];
xdotv += dx[k] * dv[k];
}
このあとはちょっと違うことがある
というのが粒子間の重力計算の本体だが、これが都合5回でてきていること
がわかる。 get_force_host は時間積分の間に重力を計算する関数で、
これはある時刻では一部の粒子への力しか計算しない。このため、重力計算
自体では対称性をつかわないで1つの粒子への力を求めるのに他の全粒子から
の力をそのまま計算する。動かす粒子の数が ni 個、全粒子の数は n
個となっている。 get_force_initial のほうは全粒子から全粒子へ
の力を計算するので、対称性を使って一度距離を計算したらそれを使って両方
の粒子の加速度を計算することで、計算量が半分とはいかないがある程度
減らしている。 get_force_host_bh のほうはブラックホール粒子からの力
を計算する。
まず、 get_force_host とその _bh がついた版のほうをみてみよう。同じコー
ドが2回でてくるならそれを関数にできるので、そうすることを考えてみる。
void pairwise_force(double xj[3],
double vj[3],
double mj,
double eps2,
double x[3],
double v[3],
double a[3],
double j[3],
double *p)
{
double r2, xdotv, r2inv, rinv, r3inv, r5inv, xdotvr5inv;
double dx[3], dv[3];
int k;
r2 = eps2;
xdotv = 0.0;
for (k = 0; k < 3; k++) {
dx[k] = xj[k] - x[k] ;
dv[k] = vj[k] - v[k] ;
r2 += dx[k] * dx[k];
xdotv += dx[k] * dv[k];
}
r2inv = 1.0 / r2 ;
rinv = sqrt(r2inv);
r3inv = rinv * r2inv;
r5inv = r3inv * r2inv;
xdotvr5inv = 3.0 * xdotv * r5inv;
for (k = 0; k < 3; k++) {
a[k] += mj * r3inv * dx[k];
j[k] += mj * (r3inv * dv[k] - xdotvr5inv * dx[k]);
}
*p -= mj * rinv;
}
/* ----------------------------------- */
void get_force_host(int ni, int nj,
double *m, double (*x1)[3], double (*v1)[3],
double (*xi)[3], double (*vi)[3],
double (*ai)[3], double (*ji)[3],
double *eps2i, int *indexi, double *pi)
{
int ii, j, k, idi;
for(ii=0;ii<ni;ii++){
for (k = 0; k < 3; k++) {
ai[ii][k] = 0.0;
ji[ii][k] = 0.0;
}
pi[ii] = 0.0;
}
// ALL <== FS
for(ii=0;ii<ni;ii++){
idi = indexi[ii];
for(j=0;j<nj;j++){
if(idi!=j){
pairwise_force(x1[j],v1[j],m[j],eps2i[ii],
xi[ii],vi[ii],ai[ii],ji[ii],pi+ii);
}
}
}
}
void get_force_from_bh(int ni, int nj,
double *m,
double (*x1)[3], double (*v1)[3],
double (*xi)[3], double (*vi)[3],
double (*ai)[3], double (*ji)[3],
double *eps2, int *indexi,double *pi,
int BHnum, int *BHi, double eps_bb2)
{
int ii, j, k,idi;
double eps2;
for(ii=0;ii<ni;ii++){
idi = indexi[ii];
// --- ALL <== BH ---
for (j = nj; j < (nj+BHnum) ; j++) {
if(idi < nj ) eps2 = eps2[j];
if(idi >= nj ) eps2 = eps_bb2;
if(idi!=j){
pairwise_force(x1[j],v1[j],m[j],eps2,
xi[ii],vi[ii],ai[ii],ji[ii],pi+ii);
}
}
}
}
こんなの。元の関数2つが 110行くらいだったのが、関数3つ合わせても
80行になってだいぶ短くなった。しかし、書き直したものを良くみると
(良くみなくても元からそうなのだが)新しく作った pairwise_force を持つ
2つの2重ループ自体が同じ構造をしています。違うのは eps2 の与えかただけ
である。ですから、この eps2 の与えかたのロジックを整理すればもうちょっと
簡単になるはずである。
与えかたのロジックを調べてみると、なんだか良くわからない。_bh がつ
かないほうでは eps2 という配列を使っている。これに対してつくほうでは、
自分のインデックスである idi が nj より大きい時は eps_bb2 というものを
使い、それ以外の場合には配列を使っている。
ということは、このプログラムでは、ブラックホール粒子は nj より先のとこ
ろに入っているということである。 このことを素直に表現するなら、
#include <stdio.h>
#include <math.h>
#include "bh.h"
/* ----------------------------------- */
void pairwise_force(double xj[3],
double vj[3],
double mj,
double eps2,
double x[3],
double v[3],
double a[3],
double j[3],
double *p)
{
double r2, xdotv, r2inv, rinv, r3inv, r5inv, xdotvr5inv;
double dx[3], dv[3];
int k;
r2 = eps2;
xdotv = 0.0;
for (k = 0; k < 3; k++) {
dx[k] = xj[k] - x[k] ;
dv[k] = vj[k] - v[k] ;
r2 += dx[k] * dx[k];
xdotv += dx[k] * dv[k];
}
r2inv = 1.0 / r2 ;
rinv = sqrt(r2inv);
r3inv = rinv * r2inv;
r5inv = r3inv * r2inv;
xdotvr5inv = 3.0 * xdotv * r5inv;
for (k = 0; k < 3; k++) {
a[k] += mj * r3inv * dx[k];
j[k] += mj * (r3inv * dv[k] - xdotvr5inv * dx[k]);
}
*p -= mj * rinv;
}
/* ----------------------------------- */
void get_force_host(int ni, int nj,
double *m, double (*x1)[3], double (*v1)[3],
double (*xi)[3], double (*vi)[3],
double (*ai)[3], double (*ji)[3],
double int *indexi, double *pi)
{
int ii, j, k, idi;
eps2[3];
for(ii=0;ii<ni;ii++){
for (k = 0; k < 3; k++) {
ai[ii][k] = 0.0;
ji[ii][k] = 0.0;
}
pi[ii] = 0.0;
}
// ALL <== FS
for(ii=0;ii<ni;ii++){
idi = indexi[ii];
if(idi!=j){
for(j=0;j<nj;j++){
if (idi<nj){
if (j<nj){
eps2=eps2_ff;
}else if (j==nj){
eps2=eps2_bh1;
}else{
eps2=eps2_bh2;
}
}else{
if (j>nj){
eps2=eps2_bh;
}else if (idi==nj){
eps2=eps2_bh1;
}else{
eps2=eps2_bh2;
}
}
pairwise_force(x1[j],v1[j],m[j],eps2,
xi[ii],vi[ii],ai[ii],ji[ii],pi+ii);
}
}
}
}
のように、多少複雑だが条件分岐で与えればいい(このコードはコンパイルで
きない)。もっとも、もうちょっとましな方法もある。以下をみてみよう。
#include <stdio.h>
#include <math.h>
#include "bh.h"
/* ----------------------------------- */
void pairwise_force(double xj[3],
double vj[3],
double mj,
double eps2,
double x[3],
double v[3],
double a[3],
double j[3],
double *p)
{
double r2, xdotv, r2inv, rinv, r3inv, r5inv, xdotvr5inv;
double dx[3], dv[3];
int k;
r2 = eps2;
xdotv = 0.0;
for (k = 0; k < 3; k++) {
dx[k] = xj[k] - x[k] ;
dv[k] = vj[k] - v[k] ;
r2 += dx[k] * dx[k];
xdotv += dx[k] * dv[k];
}
r2inv = 1.0 / r2 ;
rinv = sqrt(r2inv);
r3inv = rinv * r2inv;
r5inv = r3inv * r2inv;
xdotvr5inv = 3.0 * xdotv * r5inv;
for (k = 0; k < 3; k++) {
a[k] += mj * r3inv * dx[k];
j[k] += mj * (r3inv * dv[k] - xdotvr5inv * dx[k]);
}
*p -= mj * rinv;
}
void force_from_array(int index,
int nj;
double xj[][3],
double vj[][3],
double mj[],
double eps2,
double x[3],
double v[3],
double a[3],
double jerk[3],
double *p)
{
int j;
for(j=0;j<nj;j++){
if (j!=index){
pairwise_force(xj[j],vj[j],mj[j],eps2,xi,v,a,jerk,p);
}
}
}
void get_force_host(int ni, int nj,
double *m, double (*x1)[3], double (*v1)[3],
double (*xi)[3], double (*vi)[3],
double (*ai)[3], double (*ji)[3],
double int *indexi, double *pi)
{
int ii, j, k, idi;
for(ii=0;ii<ni;ii++){
for (k = 0; k < 3; k++) {
ai[ii][k] = 0.0;
ji[ii][k] = 0.0;
}
pi[ii] = 0.0;
}
for(ii=0;ii<ni;ii++){
double eps2 = eps2_ff;
int i = indexi[ii];
if (i<nj){
pairwise_force(x1[nj],v1[nj],m[nj],eps2_bh1,
xi[i],vi[i],ai[i],ji[i],pi+i);
pairwise_force(x1[nj+1],v1[nj+1],m[nj+1],eps2_bh2,
xi[i],vi[i],ai[i],ji[i],pi+i);
}else{
eps2 = eps2_bh1;
jbh=nj+1;
if (i>nj) {
eps2 = eps2_bh2;
jbh=nj;
}
pairwise_force(x1[jbh],v1[jbh],m[jbh],eps2_bb,
xi[i],vi[i],ai[i],ji[i],pi+i);
}
force_from_array(nj,idi,x1,v1,m,eps2,
xi[i],vi[i],ai[i],ji[i],pi+i);
}
}
ここでは、「多数の星からの近くを同じソフトニングで計算する」という機能
を関数にまとめることで、コードとして見通しがよいものにしている。これに
はおまけもあって、このような形で整理されていれば、この部分だけを SIMD
命令とか GRAPE とかいったもので高速化する、というのが容易になる。
初期化の部分では、対称性を利用して計算量を半分強にしているが、問題は
それによって得るもの、失うものがどれくらいあるかである。2倍弱計算量が
増えることを気にしなければ、ここはもちろん普通に重力を計算する関数を
全粒子について呼ぶだけですむ。実際のシミュレーションのことを考えると
初期化に必要な時間は誤差の範囲であり、この部分を2倍弱高速化することに
あまり意味はない。また、SIMD 等での高速化を図るなら、対称性の利用のた
めに別ルーチンを書いているのは余計な手間を増やすだけになる。
もしも、常に初期化のほうに時間がかかる(時間積分プログラムに場合には原
理的にありえないが)というような場合には、ここでの2倍に意味があるかもし
れないが、そういう状況であるとわかっている場合でなければここでわざわざ
手間を増やすようなことはしてはいけない。
さて、このように順番に見ていくと、何故最初の例のようなプログラムを書く
のか?ということが疑問に感じられるかもしれない。というか、疑問に感じな
いようではちょっと問題だが、でも、実際に自分が書いたプログラムのことを
ふりかえってみた時に最初の例のようなことをしていない、と言い切ることが
できる人はあまりいないのではないかと思う。
最初の例のようなプログラムができるのは、ブラックホールとかないプログラ
ムから、場当り的にプログラムを変更していってできたものをそのまま放置す
るからである。つまり、元々は
#include <stdio.h>
#include <math.h>
#include "bh.h"
#include <stdio.h>
#include <math.h>
#include "bh.h"
void get_force(int ni, int nj,
double *m, double (*x1)[3], double (*v1)[3],
double (*xi)[3], double (*vi)[3],
double (*ai)[3], double (*ji)[3],
double eps2, int *indexi, double *pi)
{
int ii, j, k, idi;
double rinv, r2, r2inv, r3inv, r5inv;
double xdotv, xdotvr5inv,r3invdx;
double r3invdvetc;
double dx[3], dv[3];
for(ii=0;ii<ni;ii++){
for (k = 0; k < 3; k++) {
ai[ii][k] = 0.0;
ji[ii][k] = 0.0;
}
pi[ii] = 0.0;
}
for(ii=0;ii<ni;ii++){
idi = indexi[ii];
for(j=0;j<nj;j++){
if(idi!=j){
r2 = eps2;
xdotv = 0.0;
for (k = 0; k < 3; k++) {
dx[k] = x1[j][k] - xi[ii][k] ;
dv[k] = v1[j][k] - vi[ii][k] ;
r2 += dx[k] * dx[k];
xdotv += dx[k] * dv[k];
}
r2inv = 1.0 / r2 ;
rinv = sqrt(r2inv);
r3inv = rinv * r2inv;
r5inv = r3inv * r2inv;
xdotvr5inv = 3.0 * xdotv * r5inv;
for (k = 0; k < 3; k++) {
r3invdx = r3inv * dx[k];
r3invdvetc = r3inv * dv[k] - xdotvr5inv * dx[k];
ai[ii][k] += m[j] * r3invdx;
ji[ii][k] += m[j] * r3invdvetc;
}
pi[ii] -= m[j] * rinv;
}
}
}
}
void get_force_initial(int n, double eps2, double *m,
double (*x0)[3], double (*v0)[3],
double (*a0)[3], double (*j0)[3], double *p)
{
int i, j, k, b;
double rinv, r2, r2inv, r3inv, r5inv;
double xdotvr5inv,r3invdx, r3invdvetc;
double xdotv, dx[3], dv[3];
for (i = 0; i < n; i++) {
for (k = 0; k < 3; k++) {
a0[i][k] = 0.0;
j0[i][k] = 0.0;
}
p[i] = 0.0;
}
for (i = 0; i < (n -BHnum); i++) {
for (j = i+1; j < (n-BHnum); j++) {
r2 = eps2[i];
xdotv = 0.0;
for (k = 0; k < 3; k++) {
dx[k] = x0[j][k] - x0[i][k] ;
dv[k] = v0[j][k] - v0[i][k] ;
r2 += dx[k] * dx[k];
xdotv += dx[k] * dv[k];
}
r2inv = 1.0 / r2 ;
rinv = sqrt(r2inv);
r3inv = rinv * r2inv;
r5inv = r3inv * r2inv;
xdotvr5inv = 3.0 * xdotv * r5inv;
for (k = 0; k < 3; k++) {
r3invdx = r3inv * dx[k];
r3invdvetc = r3inv * dv[k] - xdotvr5inv * dx[k];
a0[i][k] += m[j] * r3invdx;
a0[j][k] -= m[i] * r3invdx;
j0[i][k] += m[j] * r3invdvetc;
j0[j][k] += - m[i] * r3invdvetc;
}
p[i] -= m[j] * rinv;
p[j] -= m[i] * rinv;
}
}
}
というような感じの、まあ、力の計算が重複しているけど2つくらいで初期化
と時間積分で違うからまあいいか、というプログラムだったものに、特殊な
場合としてブラックホール1つ、もう一つ、と付け加えていった時に、さらに
ソフトニングを配列にすればよいかと思って変更したけどよく考えるとそれで
は駄目だったのでもうちょっといじって動くようにした、というのが最初の
コードである。
問題は、どのようにすれば最初のコードのようなものを書かないですむか、と
いうことだが、個人のレベルとしては、これは、
書いてしまったと
気が付いたら、あるいは
指摘されたら
直す
というだけである。最初のコードが駄目なものであるのはみればわかるし、ど
うなっているべきかの判断もそれほど難しいことではない。問題は、あくまで
も、実際に直すかどうか、である。直したほうが後で楽だが、「直す」という
作業は常に大変なことにように思え、また現在なんとか動いているものに
手をつけたくない、とか、「動いているものを変更するな」と偉い人がいって
いたとか、そういう、直さないための理由や口実を見つけてくることは簡単だ
からである。
これは、色々な人がいっているように、意志、あるいは習慣の問題で、本人が
努力して身につけるしかない。が、ここでのポイントは、最初のコードのよう
なものを長年にわたって書き続ける、ということは、それだけの間研究者なり
プログラマなりとしての人生をドブに捨てている、ということである。よりま
しなやり方をとっていれば短い時間ででき、より進歩できるはずなのに、
それができない、ということになるからである。