初期条件・境界条件・ソース関数等の導入

さてここでは自分の研究に合わせて初期条件や境界条件の関数を自分で定義して使用する方法を説明します。これまでと異なりここからはプログラミングがメインになります。

Athena++の約束事

具体的な説明に入る前に、初期条件等を書くのに最低限必要なAthena++の機能や約束事について説明します。

1. AthenaArray
Athena++では多次元のデータを格納するのにC/C++の配列ではなくAthenaArrayと呼ばれるテンプレートを使用しています。多次元配列を定義するにはまずAthenaArrayを<>の中に使用する変数の型を入れて宣言します。次に配列のメモリを確保するためにNewAthenaArrayメソッドを使用します（注：data.NewAthenaArray()の"."（ピリオド）はdataというクラスのNewAthenaArrayという関数または変数を呼び出すという意味です。もしdataがクラスへのポインタであればdata->NewAthenaArrayとなることに注意して下さい）。NewAthenaArrayは1～5個の変数を取り、それぞれ1～5次元配列の各方向の要素数を渡します。配列の配置はC/C++風であり、一番最後の変数がメモリ上連続に配置されます。確保した配列にはdata(k,j,i)のように（Fortranの配列とよく似た文法で）読み書きすることができます。使用後はDeleteAthenaArray()を呼んで配列を解放します。

AthenaArray<Real> data; // 倍精度浮動小数点型
AthenaArray<int> nums; // 整数型
data.NewAthenaArray(64,64,64); // 64^3 3D array
nums.NewAthenaArray(256); // 1D array
for(int k=0; k<64; k++) {
  for(int j=0; j<64; j++) {
    for(int i=0; i<64; i++) {
      data(k,j,i)=(Real)(i+j+k);
      std::cout << "data = " << data(k,j,i) << std::endl;
    }
  }
}
data.DeleteAthenaArray();
nums.DleteAthenaArray();

2. 格子の構造とクラスの関係
Athena++では計算領域を"Mesh"と呼び、全てのプロセスがこの情報を保持するMeshクラスを持っています。更にその領域を（論理的に）同じ形をした"MeshBlock"という単位に分割し、各プロセスが一つまたは複数のMeshBlockを所有しています。このMeshBlockはMeshクラス以下に配列として格納されています。更にMeshBlock間の位置関係を計算するためのツリー構造がMeshBlockTreeクラスに構築されており、これもMeshクラスに格納されています。各MeshBlockはその下に座標を格納・計算するCoordinatesクラス、流体力学部の情報を格納・計算するHydroクラス、磁場の情報を格納・計算するFieldクラス、境界条件を処理するBoundaryValuesクラス、状態方程式を処理するEquationOfStateクラスなどを保持しています。AMR・SMRにおけるリファインメントのレベルもMeshBlockのパラメータに過ぎません。長々と書いてきましたが、この概念を図示すると以下のようになります。
Mesh and Class Structure

3. 有用な変数の格納場所とその意味
初期条件や境界条件その他を操作する際に様々な情報にアクセスする必要があります。以下に良く使う情報の格納場所と型を列挙します。

u,wの最初の次元（変数）は物理量を表します。各成分はuに対しては密度=IDN、運動量=IM1,IM2,IM3、（定義されている場合）全エネルギー=IEN、wに対しては密度=IDN、速度＝IVX,IVY,IVZ、ガス圧力=IPRを用いて指定します。同様にbccの最初の次元は方向を表し、IB1,IB2,IB3の順にx1,x2,x3方向に対応しています。

4. 入力パラメータの読み取り
Athena++ではファイルから読み取ったパラメータをInputParametersクラスの連想配列に格納します。この読み出しには以下のような関数を使います。<problem>ブロック以外の計算自体に関わるパラメータも読みだすことができます。 

ParameterInput *pin;
// ...
Real r0    = pin->GetReal("problem","radius"); // Real型の<problem>のradius変数を読み取る
Real amp   = pin->GetOrAddReal("problem","amp",0.0); // 同様、ただし未定義なら0.0をセット
Real mode  = pin->GetInteger("problem","mode"); // 整数型の<problem>のmode変数を読み取る
Real nw    = pin->GetOrAddInteger("problem","nwave",nx1); // 同様、ただし未定義ならnx1
std::string problemid, refinement; // #include <string>が必要
problemid  = pin->GetOrString("job","problem_id"); // 文字列型のproblem IDの取得
refinement = pin->GetOrAddString("mesh","refinement","none"); // 同様、未定義なら"none"

5. 座標の定義
Athena++でのセルと座標の定義は以下の図のようになっています。
Cell structure
セル中心のインデックスが(i,j,k)の時、その中心の座標は(x1v(i), x2v(j), x3v(k))、その左側の表面の座標はそれぞれx1f(i), x2f(j), x3f(k)、右側の表面の座標はx1f(i+1), x2f(j+1), x3f(k+1)となっています。その他計算に必要な体積や表面積、辺の長さは以下のCoordinatesクラスの関数で取得できます。 

Coordinates *pcoord;
AthenaArray<Real> len, area, vol;
//...
//(k,j,il～iu)の値をx1方向の一次元配列として連続的に取得する関数群
pcoord->CellVolume(k, j, il, iu, vol);
pcoord->Face1Area(k, j, il, iu, area);
pcoord->Face2Area(k, j, il, iu, area);
pcoord->Face3Area(k, j, il, iu, area);
pcoord->Edge1Length(k, j, il, iu, len);
pcoord->Edge2Length(k, j, il, iu, len);
pcoord->Edge3Length(k, j, il, iu, len);
std::cout << "cell volume = " << vol(i) << std::endl;
//(k,j,i)一点の値を返す関数群
Real cellvol=pcoord->GetCellVolume(k, j, i);
Real x1face=pcoord->GetFace1Area(k, j, i);  // (i-1/2, j,     k    )
Real x2face=pcoord->GetFace2Area(k, j, i);  // (i,     j-1/2, k    )
Real x3face=pcoord->GetFace3Area(k, j, i);  // (i,     j,     k-1/2)
Real e1len=pcoord->GetEdge1Length(k, j, i); // (i,     j-1/2, k-1/2)
Real e2len=pcoord->GetEdge2Length(k, j, i); // (i-1/2, j,     k-1/2)
Real e3len=pcoord->GetEdge3Length(k, j, i); // (i-1/2, j-1/2, k    )

6. 計算設定のフラグ・グローバル変数
初期条件の計算に当たってコードがどのように設定されているかを知る必要があります。これらのフラグはsrc/defs.hppに格納されています。重要なものは以下の二つです。

またMPI並列計算を実行している時に各プロセスのランクや並列数が必要になることがあります。これらはMPIの関数を呼び出しても取得できますが、globals.hpp内のGlobals::my_rank, Globals::nranksを利用することもできます。

7. 密度・圧力フロア
特に高マッハ数、強磁場、強い冷却を含むような計算は密度や温度が負になり計算が破綻してしまうケースがあります。しかしこのような破綻は計算領域のごく一部だけで発生し、そこさえ乗り越えれば大域的な計算に（あまり）悪影響を与えないということも多々あります。このようなケースを避けるためにしばしば密度や圧力に最低値を設定することがあります。Athena++では以下のように密度または圧力フロアを指定することができます。これはバッドノウハウに類するテクニックなので使用する場合には十分注意して、実際に大きな影響がないことを確認して使用して下さい。 

<hydro>
dfloor = 0.0001
pfloor = 0.0001

Athena++ではデフォルトのフロアとして密度・圧力共に10-35程度の値が指定されています。通常これらは十分小さい値ですが、計算をcgs単位系で行う場合などはこのフロアでは大きすぎる場合があります。その場合は計算に現れるよりも十分小さい値をフロアに指定してください。ただし、次の項目で述べるように、計算には極端に大きいまたは小さい値が現れないような単位系を取ることをお勧めします。

8. 数値計算用の単位系
Athena++では計算に使用する単位系を指定していません。そのためcgsやSI等「普通の」単位系で計算をすることも可能ですが、宇宙物理学では通常の単位系を使うと値が非常に大きくなったり小さくなったりすることがあります。計算機で使用する浮動小数点変数で表現できる範囲には制限があるため、極端に大きいまたは小さい数が計算に現れるのは好ましくありませんし、計算結果を見る際にも直観的ではありません。そのため、通常は対象とする系に典型的なスケールが数値計算用の単位として用いること（これを規格化 normalizationと言います）をお勧めします。といってもピンとこない人もいるかと思いますので、具体的に簡単な例で計算してみましょう。まず磁場のない流体力学の場合を考えます。
Hydrodynamics Equations
ここに現れる物理量が全て（例えば）もともとcgs系で定義されていたとしましょう。これらの変数のうち流体力学方程式の場合は任意の三つ（エネルギーを解かない場合は二つ）の次元について、好みの単位を取ることができます。ここでは密度ρ、長さL、速度uについて単位を取り直すことにします。ここで、新しい単位とする値をそれぞれρ0、L0、u0とします。例えば密度としてはρ0=1/cc * mH、長さとしてはL0=1pc、速度の単位は典型的な温度での音速 u0=cs(T=100K)を与えるなどが典型的な例です。元の物理量から新しい単位系で測った物理量に変換するためには元の物理量をこれらの新しい単位で割ってやればよく、元の単位系で記述された物理量を以下のように書き直すことができます。
Scaling Three Primary Variables
ここでρ'、L'、u'などの'を付けた変数は新しい単位で測った物理量を表しています。この時に、他の物理量も同時に、新しい単位系に整合的になるように変換する必要があります。例えば質量M、時間t、圧力pなどの単位は密度・長さ・速度の組み合わせで記述できますから、新しい単位系では以下のような規則によって変換されます。 Scaling Other Variables
新しい単位系での流体力学方程式を得るためには、元の方程式を新しい単位で記述した変数に全て置き換えるだけです。ただしこの際注意が必要なのは、時間・位置での微分演算子もそれぞれ1/t、1/Lの次元を持つため、整合的に変換する必要があるということです。あとは単純な掛け算割り算の計算を行えば以下のような新しい単位系で記述された流体力学方程式が得られます。
Scaled Hydrodynamics Equations
この結果を見ると、[ ]の中の新しい単位系で記述された流体力学方程式はもとのcgs/SI系で記述された流体力学方程式と全く同じ（全体として定数倍されただけ）形をしていることがわかります。即ちどんな単位で物理量を記述しようとも単位の変換が整合的であれば方程式の形は変わらないため、同じシミュレーションコードをそのまま適用することができるのです。

というわけで、Athena++でシミュレーションを行う際にはその単位はユーザーが自身で決めることができますし、逆に言えばどのような単位系を設定するかは全てユーザーに任されています。上で述べたのは単位の決め方の一例であって、流体力学方程式であれば独立な変数三つについては自由に決めることができます（例えば質量・時間・長さなどでも可）。

磁気流体力学の場合でも基本は同じです。
MHD Equations
ここで注意点が二つあります。まず、天文学では電磁場の方程式の基本的な単位系としてcgs Gauss単位系を用いることが慣例ですが、Athena++をはじめ多くのシミュレーションコードでは方程式を更に簡単にするために磁場をB/sqrt(4π)→Bのように置き換えたものを使用しています。これは言い換えればsqrt(4π)Gaussを磁束密度の基本的な単位としていることに対応しており、Gauss等次元量で磁束密度を表現したいときは注意が必要です。次に磁場Bの二乗B2が圧力と同じ次元を持っているため磁場の単位は独立ではなく、磁場の変換は任意に決めることはできません。上の例でいえば磁束密度は
MHD Equations
のようにスケールされることに留意してください（1/2乗？と思うかもしれませんが、そのまま計算して問題ありません）。なお数値計算ではGaussなどの物理的な単位系を用いるよりも、より直観的な無次元量であるプラズマβ=pgas/pmagやAlfven速度と音速の比a=valf/cs等を用いて磁場の強度を指定・測定することが多々あります。

なお追加の物理過程を導入する場合は注意が必要です。例えば輻射輸送または放射冷却を導入するとダストや分子・原子のオパシティという特徴的なスケールを持つ量が出てきます。この場合でも方程式を別の単位系に変換することは可能ですが、オパシティなども整合的になるように変換する必要があることに注意してください。他によく使う計算用の単位系としては、点源重力場が存在する場合GM=1としたり、自己重力を考慮する場合には4πG=1とすることがあります。

初期条件ファイル

ここでは初期条件を生成するプログラムの書き方について説明します。初期条件を生成するソースコード（problem generatorと呼びます）はsrc/pgenディレクトリに格納されており、その中の使用するファイルをコンフィギュレーションスクリプトの--probオプションで指定します。

Problem generatorファイルには以下のような関数を含みます。

これらを実現するために以下のようなインターフェースが提供されており、これらのうち必要なものを記述します（定義されない場合はデフォルト＝なにもしないの物が使用されます）。具体例はsrc/pgen以下の各ファイルを参考にして下さい。

MeshBlock::ProblemGenerator

Athena++で計算を行うのに最も重要なのが初期条件を設定するこのProblemGenerator関数です。この関数ではアクティブ領域の保存量（密度、運動量、全エネルギー、phydro->u）とMHD計算の場合はセル表面の磁場（pfield->bの3成分）を全て指定する必要があります。セル中心の物理量についてはi=is～ie, j=js～je, k=ks～ke、磁場についてはこれに加えて各法線方向について1セル多く（x1方向ならi=is～ie+1）の値を計算して下さい。MHDの場合セル中心での磁場（pfield->bcc）は計算上必要ありませんが、セル中心での全エネルギーを計算する際に必要になる場合があります。primitive変数についてはコードによって自動的に計算されます。またMHD計算ではdivB=0を全域で満たさなければなりません。そのため、自明でない磁場構造を初期条件に与える場合にはベクトルポテンシャルを用いて計算する必要があります。またSMR等を利用する場合にレベル間で整合的な初期条件になるように指定しなければなりません（これについてはsrc/pgen/linear_wave.cppの例を参照してください）。以下は線形波動伝播テスト計算（を簡略化したもの）の例で、各格子の辺上でベクトルポテンシャル（A1,A2,A3は指定された座標でのベクトルポテンシャルを与える関数）を計算し、そこからStokesの定理を用いて各セル表面の磁場を計算しています。（注：この例では初期磁場は空間的に一定(bx0,by0,bz0)なので自明ですがあえてベクトルポテンシャルを使用しています。）

void MeshBlock::ProblemGenerator(ParameterInput *pin)
{
  // Initialize the magnetic fields.  Note wavevector, eigenvectors
  // and other variables are set in InitUserMeshData

  if (MAGNETIC_FIELDS_ENABLED) {
    AthenaArray a1,a2,a3;
    int nx1 = (ie-is)+1 + 2*(NGHOST);
    int nx2 = (je-js)+1 + 2*(NGHOST);
    int nx3 = (ke-ks)+1 + 2*(NGHOST);
    a1.NewAthenaArray(nx3,nx2,nx1);
    a2.NewAthenaArray(nx3,nx2,nx1);
    a3.NewAthenaArray(nx3,nx2,nx1);

    // wave amplitudes
    dby = amp*rem[NWAVE-2][wave_flag];
    dbz = amp*rem[NWAVE-1][wave_flag];

    // calculate vector potential
    for (int k=ks; k<=ke+1; k++) {
      for (int j=js; j<=je+1; j++) {
        for (int i=is; i<=ie+1; i++) {
          a1(k,j,i) = A1(pcoord->x1v(i), pcoord->x2f(j), pcoord->x3f(k));
          a2(k,j,i) = A2(pcoord->x1f(i), pcoord->x2v(j), pcoord->x3f(k));
          a3(k,j,i) = A3(pcoord->x1f(i), pcoord->x2f(j), pcoord->x3v(k));
        }
      }
    }

    // Initialize interface fields
    for (int k=ks; k<=ke; k++) {
      for (int j=js; j<=je; j++) {
        for (int i=is; i<=ie+1; i++) {
          pfield->b.x1f(k,j,i) = (a3(k  ,j+1,i) - a3(k,j,i))/pcoord->dx2f(j) -
                                 (a2(k+1,j  ,i) - a2(k,j,i))/pcoord->dx3f(k);
        }
      }
    }

    for (int k=ks; k<=ke; k++) {
      for (int j=js; j<=je+1; j++) {
        for (int i=is; i<=ie; i++) {
          pfield->b.x2f(k,j,i) = (a1(k+1,j,i  ) - a1(k,j,i))/pcoord->dx3f(k) -
                                 (a3(k  ,j,i+1) - a3(k,j,i))/pcoord->dx1f(i);
        }
      }
    }

    for (int k=ks; k<=ke+1; k++) {
      for (int j=js; j<=je; j++) {
        for (int i=is; i<=ie; i++) {
         pfield->b.x3f(k,j,i) = (a2(k,j  ,i+1) - a2(k,j,i))/pcoord->dx1f(i) -
                                (a1(k,j+1,i  ) - a1(k,j,i))/pcoord->dx2f(j);
        }
      }
    }
    a1.DeleteAthenaArray();
    a2.DeleteAthenaArray();
    a3.DeleteAthenaArray();
  }

  // initialize conserved variables
  for (int k=ks; k<=ke; k++) {
    for (int j=js; j<=je; j++) {
      for (int i=is; i<=ie; i++) {
        Real x = cos_a2*(pcoord->x1v(i)*cos_a3 + pcoord->x2v(j)*sin_a3)
               + pcoord->x3v(k)*sin_a2;
        Real sn = sin(k_par*x);

        phydro->u(IDN,k,j,i) = d0 + amp*sn*rem[0][wave_flag];

        Real mx = d0*vflow + amp*sn*rem[1][wave_flag];
        Real my = amp*sn*rem[2][wave_flag];
        Real mz = amp*sn*rem[3][wave_flag];

        phydro->u(IM1,k,j,i) = mx*cos_a2*cos_a3 - my*sin_a3 - mz*sin_a2*cos_a3;
        phydro->u(IM2,k,j,i) = mx*cos_a2*sin_a3 + my*cos_a3 - mz*sin_a2*sin_a3;
        phydro->u(IM3,k,j,i) = mx*sin_a2                    + mz*cos_a2;

        if (NON_BAROTROPIC_EOS) {
          phydro->u(IEN,k,j,i) = p0/gm1 + 0.5*d0*u0*u0 + amp*sn*rem[4][wave_flag];
          if (MAGNETIC_FIELDS_ENABLED) {
            phydro->u(IEN,k,j,i) += 0.5*(bx0*bx0+by0*by0+bz0*bz0);
          }
        }
      }
    }
  }

  return;
}

Mesh::InitUserMeshData

この関数はMesh全体に関わるデータを確保・初期化したり、ユーザー定義関数を指定するためのものです。この関数は計算の開始時、またはリスタート時に実行されます。主に以下のようなユーザー定義関数・変数が指定できます。詳細はこの後の個別の説明を参照してください。

境界条件の指定

Athena++ではデフォルトの境界条件として周期境界(periodic)、流出境界（outflow、より正確にはゼロ勾配）、反射境界(reflecting)の三つが提供されており、これらはパラメータファイル中で指定するだけで使用することができます。更に特殊な境界条件として3次元球面極座標の場合、x2(θ)方向の境界で計算領域が極軸に接している（θ=0またはπ=3.141592653589793）の場合に限り極境界（polar）を利用することができます。これは極軸を自然に扱うためのもので、幾つか条件（φ方向が0～2πの周期境界でなければならない、φ方向の格子点数が偶数でなければならない、φ方向のMeshBlock数が1または偶数でなければならない）があります。

これらで対応できない一般の境界条件を設定するためにはユーザー定義境界条件を用います。このためにはまずパラメータファイルで境界条件にuserを指定します。その上で以下の型に従った境界条件関数を定義し、Mesh::InitUserMeshData内で設定（enroll）します。 

void MyBnd_ix1(MeshBlock *pmb, Coordinates *pco, AthenaArray<Real> &prim, FaceField &b,
               Real time, Real dt, int is, int ie, int js, int je, int ks, int ke);
//...
void Mesh::InitUserMeshData(ParameterInput *pin)
{
  //...
  EnrollUserBoundaryFunction(BoundaryFace::inner_x1, MyBnd_ix1);
  //...
  return;
}

ユーザー定義境界条件関数内では各方向毎に指定された範囲内のprimitive variablesを指定します。以下はx1の左端から細いジェットを入射させるような境界条件の例です。

void MyBnd_ix1(MeshBlock *pmb, Coordinates *pco, AthenaArray<Real> &prim, FaceField &b,
               Real time, Real dt, int is, int ie, int js, int je, int ks, int ke)
{
  // set primitive variables in inlet ghost zones
  for(int k=ks; k<=ke; ++k){
  for(int j=js; j<=je; ++j){
    for(int i=1; i<=(NGHOST); ++i){
      Real rad = sqrt(SQR(pco->x2v(j)-x2_0) + SQR(pco->x3v(k)-x3_0));
      if(rad <= r_jet){
        prim(IDN,k,j,is-i) = d_jet;
        prim(IVX,k,j,is-i) = vx_jet;
        prim(IVY,k,j,is-i) = vy_jet;
        prim(IVZ,k,j,is-i) = vz_jet;
        prim(IPR,k,j,is-i) = p_jet;
      } else{
        prim(IDN,k,j,is-i) = prim(IDN,k,j,is);
        prim(IVX,k,j,is-i) = prim(IVX,k,j,is);
        prim(IVY,k,j,is-i) = prim(IVY,k,j,is);
        prim(IVZ,k,j,is-i) = prim(IVZ,k,j,is);
        prim(IPR,k,j,is-i) = prim(IPR,k,j,is);
      }
    }
  }}

  // set magnetic field in inlet ghost zones
  if (MAGNETIC_FIELDS_ENABLED) {
    for(int k=ks; k<=ke; ++k){
      for(int j=js; j<=je; ++j){
        for(int i=1; i<=(NGHOST); ++i){
          Real rad = sqrt(SQR(pco->x2v(j)-x2_0) + SQR(pco->x3v(k)-x3_0));
          if(rad <= r_jet){
            b.x1f(k,j,is-i) = bx_jet;
          } else{
            b.x1f(k,j,is-i) = b.x1f(k,j,is);
          }
        }
      }
    }

    for (int k=ks; k<=ke; ++k) {
      for (int j=js; j<=je+1; ++j) {
        for (int i=1; i<=(NGHOST); ++i) {
          Real rad = sqrt(SQR(pco->x2v(j)-x2_0) + SQR(pco->x3v(k)-x3_0));
          if(rad <= r_jet){
            b.x2f(k,j,is-i) = by_jet;
          } else{
            b.x2f(k,j,is-i) = b.x2f(k,j,is);
          }
        }
      }
    }

    for (int k=ks; k<=ke+1; ++k) {
      for (int j=js; j<=je; ++j) {
        for (int i=1; i<=(NGHOST); ++i) {
          Real rad = sqrt(SQR(pco->x2v(j)-x2_0) + SQR(pco->x3v(k)-x3_0));
          if(rad <= r_jet){
            b.x3f(k,j,is-i) = bz_jet;
          } else{
            b.x3f(k,j,is-i) = b.x3f(k,j,is);
          }
        }
      }
    }
  }
  return;
}

境界条件関数で指定すべき領域のループ範囲は以下のテーブルの通りです。ここではNGHOST=2を仮定していますので、より大きいNGHOSTを使用している場合にはinner側は下端を、outer側は上端をそれぞれ必要な数だけ増やす必要があります。

Primitive variables i-loweri-upperj-lowerj-upperk-lowerk-upper
ix1is-2is-1jsjekske
ox1ie+1ie+2jsjekske
ix2isiejs-2js-1kske
ox2isieje+1je+2kske
ix3isiejsjeks-2ks-1
ox3isiejsjeke+1ke+2
b.x1f i-loweri-upperj-lowerj-upperk-lowerk-upper
ix1is-2is-1jsjekske
ox1ie+2ie+3jsjekske
ix2isie+1js-2js-1kske
ox2isie+1je+1je+2kske
ix3isie+1jsjeks-2ks-1
ox3isie+1jsjeke+1ke+2
b.x2f i-loweri-upperj-lowerj-upperk-lowerk-upper
ix1is-2is-1jsje+1kske
ox1ie+1ie+2jsje+1kske
ix2isiejs-2js-1kske
ox2isieje+2je+3kske
ix3isiejsje+1ks-2ks-1
ox3isiejsje+1ke+1ke+2
b.x3f i-loweri-upperj-lowerj-upperk-lowerk-upper
ix1is-2is-1jsjekske+1
ox1ie+1ie+2jsjekske+1
ix2isiejs-2js-1kske+1
ox2isieje+1je+2kske+1
ix3isiejsjeks-2ks-1
ox3isiejsjeke+2ke+3

境界条件は些細なようで数値計算を支配する極めて重要なものです。特に磁場を含む場合境界条件によっては計算が不安定になることがあり、安定な計算を行うために試行錯誤が必要になる場合があります。このあたりは多少経験と物理的な勘が必要になりますが、まずは試してみましょう。

ソース関数と時間刻み

次にユーザー定義のソース関数、例えば局所的な加熱冷却関数などを実装する方法を紹介します。ただしAthena++が提供するインターフェースは陽解法で計算することができる源泉項に限ります（サブサイクリングを無理やり実装することはできますが、時間精度が一次精度に落ちてしまいます。また陰解法が必要な場合には独自の関数を記述する必要があります）。

ユーザー定義ソース関数は以下のプロトタイプを持ち、境界条件と同様にInitUserMeshData関数内でEnrollUserExplicitSourceFunctionを用いて指定します。

void MySource(MeshBlock *pmb, const Real time, const Real dt,
              const AthenaArray<Real> &prim, const AthenaArray<Real> &bcc,
              AthenaArray<Real> &cons);
//...
void Mesh::InitUserMeshData(ParameterInput *pin)
{
  //...
  EnrollUserExplicitSourceFunction(MySource);
  //...
  return;
}

ソース関数内では引数として渡されたprimitive変数(prim)とセル中心の磁場(bcc)を用いて保存量(cons)をアップデートします。例えばある一定のタイムスケールtauで目標温度t0に戻るような冷却関数は以下のように記述することができます。

void cooling(MeshBlock *pmb, const Real time, const Real dt,
             const AthenaArray<Real> &prim, const AthenaArray<Real> &bcc,
             AthenaArray<Real> &cons)
{
  Real g = pmb->peos->GetGamma();
  Real tau = 0.01;
  Real t0 = 10.0;
  for (int k=pmb->ks; k<=pmb->ke; ++k) {
    for (int j=pmb->js; j<=pmb->je; ++j) {
      for (int i=pmb->is; i<=pmb->ie; ++i) {
        Real temp = (g-1.0)*prim(IPR,k,j,i)/prim(IDN,k,j,i);
        cons(IEN,k,j,i) -= dt*prim(IDN,k,j,i)*(temp - t0)/tau/(g-1.0);
      }
    }
  }
  return;
}

前述の通りこの関数は陽的に積分していますので、時間刻みは冷却時間tauを十分分解する必要があります（あるいはサブサイクルを用いる）。ユーザー定義の時間刻みを指定する場合は以下のようにします。この関数ではあるMeshNBlock内で満たすべき最小の時間刻みを計算して返します。なおこの数値にはCFL数はかからないことに注意して下さい。

Real MyTimeStep(MeshBlock *pmb)
{
  Real min_dt=FLT_MAX;
  for (int k=pmb->ks; k<=pmb->ke; ++k) {
    for (int j=pmb->js; j<=pmb->je; ++j) {
      for (int i=pmb->is; i<=pmb->ie; ++i) {
        Real dt;
        dt = ... // calculate your own time step here
        min_dt = std::min(min_dt, dt);
      }
    }
  }
  return min_dt;
}

//...
void Mesh::InitUserMeshData(ParameterInput *pin)
{
  //...
  EnrollUserTimeStepFunction(MyTimeStep);
  //...
  return;
}

なお球座標または2次元の円筒座標系における中心点源重力と、任意の座標系での一方向への一様重力加速度はコードに内蔵されており、前者は<problem>のGMパラメータに有限の値を指定することで、後者は<hydro>ブロックにgrav_acc1, grav_acc2, grav_acc3パラメータを指定すると自動的に有効になります。

MeshBlock::InitUserMeshBlockData

この関数はMesh::InitUserMeshDataとよく似ていますが、各MeshBlockについて主にユーザー定義のデータ領域を確保するために使用します。MeshBlock::ProblemGeneratorとの違いはProblemGeneratorが計算の最初に一度だけ呼ばれるのに対しこの関数はリスタートの際も実行されるため、計算に必要な一時データを記憶するためのメモリを確保したり必要なファイルを読み込むために利用することができます。

ユーザー定義データ領域

ユーザーは任意の場所で計算に必要なデータ領域を確保することができますが、Athena++ではシミュレーションをリスタートする際に自動的に保存・復元されるデータ領域を確保する機能を提供しています。このデータ領域はRealまたは整数型のAthenaArrayとして定義され、MeshクラスまたはMeshBlockクラス内に任意のサイズを確保することができます。この機能を利用するにはMesh::InitUserMeshDataまたはMeshBlock::InitUserMeshBlockData内で以下の関数で必要なデータ領域の数を宣言します。

すると以下のように宣言した数のAthenaArrayが確保されるので、NewAthenaArrayを用いて必要なデータ領域を確保します。

ちょっとわかりにくいと思うので、例を示します。以下では2つのReal型の領域を確保し、一つは１変数だけ、もう一つは16x16の二次元の領域を確保しています。

void Mesh::InitUserMeshData(ParameterInput *pin)
{
  AllocateRealUserMeshDataField(2);
  ruser_mesh_data[0].NewAthenaArray(1);
  ruser_mesh_data[1].NewAthenaArray(16,16);
  ruser_mesh_data[0](0)=0.0;
  for(int j=0; j<16; j++) {
    for(int i=0; i<16; i++) {
      // initialization
      ruser_mesh_data[1](j,i)=0.0;
    }
  }
  return;
}

Meshクラス内のユーザー定義データ領域は全プロセスで共通の大域的な情報を格納することを意図しています。このデータはリスタートファイル内に自動的に保存されますが、その際にはマスターノード（rank 0）のデータが利用され、プロセス間の整合性を取るのはユーザーの責任になります。

同様に、三次元の整数型の領域を各MeshBlock内に確保するためには

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin)
{
  AllocateIntUserMeshBlockDataField(1);
  iuser_meshblock_data[0].NewAthenaArray(16,16,16);
  for(int k=0; k<16; k++)
    for(int j=0; j<16; j++) {
      for(int i=0; i<16; i++) {
        // initialization
        iuser_meshblock_data[0](k,j,i)=0;
      }
    }
  }
  return;
}

とすればOKです。この領域は個々のMeshBlock内で独立であり、共有はされません。

これらの領域は計算終了時に自動的に解放されるため、手動で削除する必要はありません。

MeshBlock::UserWorkInLoop

この関数は計算中毎ステップの最後に呼び出され（注：サブステップの途中では実行されません）、主に計算途中のデータを解析するために使用することを目的としています。計算データを操作することは一般には推奨されません - そのような場合にはユーザー定義のソース関数または境界条件を利用して下さい。

MeshBlock::UserWorkBeforeOutput

この関数はUserWorkInLoopとよく似ていますが、ファイル出力の直前にだけ実行されます。そのため、次項で説明するユーザー定義出力変数を計算するのに便利です。

Mesh::UserWorkAfterLoop

この関数は計算終了時に一度だけ呼ばれる関数で、確保したメモリの後処理や計算結果の集計作業に用いることを意図しています。

ユーザー定義出力変数

ユーザー定義出力変数(コード内ではuser_out_var)はデータ出力に特化したもう一つのデータ領域で、VisIt等で追加の変数を他の変数と同様に可視化するのに利用することができます。ユーザー定義出力変数は各セルごとに定義され、MeshBlockと同じ形のAthenaArrayとして確保されます。この機能を有効するには、まず必要な変数の数をMeshBlock::InitUserMeshBlockData関数内でAllocateUserOutputVariables関数で宣言します。 

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin)
{
  AllocateUserOutputVariables(2);
  return;
}

この関数によりMeshBlock::nuser_out_varが2に設定され4次元のAthenaArray<Real> user_out_var配列が確保されます。次にこの配列にデータを格納するためのコードをMeshBlock::UserWorkInLoopに記述します。次の例ではガス温度とプラズマβを格納しています（注：これらはVisItやPythonでも簡単に計算できるので自明な例です）:

void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin)
{
  for(int k=ks; k<=ke; k++) {
    for(int j=js; j<=je; j++) {
      for(int i=is; i<=ie; i++) {
        Real pmag = 0.5*(SQR(pfield->bcc(IB1,k,j,i))
                        +SQR(pfield->bcc(IB2,k,j,i))
                        +SQR(pfield->bcc(IB3,k,j,i)));
        user_out_var(0,k,j,i) = phydro->w(IPR,k,j,i)/phydro->w(IDN,k,j,i);
        user_out_var(1,k,j,i) = phydro->w(IPR,k,j,i)/pmag;
      }
    }
  }
}

なおゴーストゾーンの値も出力する場合はこれらのループの範囲はis-NGHOSTからie+NGHOST等に変更してください。user_out_varは計算終了時に自動的に破棄されますので、ユーザーがこれらを解放する必要はありません。

この出力を有効にするにはインプットファイルで以下のように指定します。

<output2>
file_type   = hdf5
dt          = 0.1
variable    = uov  # (or user_out_var)

課題

さて、この講習会の目標はAthena++を用いて自身の科学的計算を行うことです。講習会の残りの時間を使って、自身の興味ある問題を実装してみて下さい。こういうことがしてみたいのだけれどどう実現すればいいか、など気軽に相談して下さい。

今すぐ思いつかない、またはアイディアはあるが現状のAthena++の機能では実現できない等の場合は、こちらから以下のような例題を提案しますので、お好みで選んで実装してみて下さい。これらは現在ではそのまま論文になるほどではなりませんが、良いアイディアがあれば論文となるような研究の基礎に十分なり得るテーマです。あるいは、提供されているpgen内のソースコードのうち興味ある問題をまず解読し、改造を試みるのも良いでしょう。

Papaloizou & Pringle不安定性
磁場のない比角運動量のトーラスが非軸対称に歪む不安定性で、現在は他のより強力な不安定性が知られているためややアカデミックな問題だと考えられていますが、かつては角運動量輸送やQPOとの関係でよく研究されていました。Papaloizou & Pringle 1984が原典で、Stone et al. 1999は２次元の流体計算ですが、この初期条件を３次元化して非軸対称の摂動を加えることでそのままこの不安定性の初期条件に使えます。球座標で計算してみましょう。トーラス部にSMRを利用して高解像度の格子を配置するのもお勧めです。例を配布したコード（2021年3月12日更新）内のpgen/sphtorus.cpp及びinputs/hydro/athinput.torusに用意しましたが、できればまずは自分で書いてみて下さい。なおこれらのサンプルは比較的短時間で終わるよう分解能を下げていますが、乱流を十分正確に分解するには高い解像度が必要です。

磁気回転不安定性
ケプラー回転等、差動回転している磁場を持つ円盤は多くの場合磁気回転不安定により乱流状態になります。この不安定性は一般的かつ強力で、降着円盤内の角運動量・質量輸送を担う重要な仮定であると考えられています。非常に幅広い研究がこれまでになされていますが、上のPP不安定性の初期条件にそのまま磁場を入れた例がStone & Pringle 2001で調べられています。まずは３次元でトロイダル（φ方向）磁場のみ入れて計算するのが簡単ですが、ベクトルポテンシャルを用いてポロイダル磁場を初期条件とする計算も是非試してみて下さい。この問題も例を配布したコード内のpgen/mhdtorus.cpp及びinputs/mhd/athinput.torusに用意しましたが、できればまずは自分で書いてみて下さい。

円盤-惑星相互作用
原始惑星系円盤中に惑星があると、その重力摂動により円盤に密度波が立ちます。この密度波は角運動量を輸送するため惑星や円盤の進化に影響を及ぼします。2次元円筒座標で円盤を表現し、ソース関数で惑星の重力を導入して密度波が生じる様子を確認しましょう。de Val-Borro et al. 2006に初期条件・境界条件の説明や多数の計算例があります。

衝撃波と熱的不安定性
熱的不安定による乱流の駆動は星間乱流の起源として有望視されています。衝突する衝撃波の熱的不安定性は超新星爆発と星間物質の相互作用のモデルとして非常によく研究されています。例としてInoue & Inutsuka 2008などを参考に加熱冷却関数を実装して衝撃波を衝突させ、自然に乱流を駆動させてみましょう。

Parker不安定性
成層大気の中で磁力線が浮上するParker不安定性は太陽表面での磁場の活動や銀河のダイナモなどで重要になる基礎的なプロセスです。やや古典的ですがMatsumoto et al. 1993は太陽表面でのParker不安定性のモデルで、今日の計算機であれば同等以上の結果が得られるはずです。

惑星風と恒星風の相互作用
中心星からの輻射などにより短周期のガス惑星の表面からは質量放出が起こりますが、その分布は輻射の当たり方等により非等方になります。Stone and Proga 2009は惑星からのウィンドと恒星からのウィンドが相互作用する様子を2次元軸対称シミュレーションしたものです。例えばウィンドに角運動量を与えたり磁場の影響を考慮する等様々な拡張が考えられます。

Jetと星間物質の相互作用
原始星やコンパクトオブジェクトからのジェットと星間物質の相互作用もしばしば観測され良く研究されているテーマです。単純に星間空間にジェットを導入するだけでも衝撃波やKH不安定性等の構造が現れますが、そこに加熱冷却関数や星間物質の非一様性を取り入れればより複雑な研究になり得ます。Stone & Hardee 2000等が参考になります。