Athena++

Static/Adaptive Mesh Refinement

宇宙物理学で扱う流体力学の問題では多くの場合重力による収縮や高マッハ数の衝撃波による圧縮などで、しばしば非常に小さな構造が発生します。問題によりますが、興味のある領域が計算領域の一部である場合はその領域にだけ高解像度の配置することで、比較的少ない計算量で必要な解像度を実現することができます。特にそのような構造が固定されている場合はStatic Mesh Refinement(SMR)で十分ですが、構造が様々な場所で発生したり移動する場合はAdaptive Mesh Refinement(AMR)を用いる必要があります。ここではこれらの技術を使った計算例について説明します。

Static Mesh Refinement

SMRでは事前に指定した領域に高解像度の格子を配置します。これはAMRと比べると柔軟ではありませんが、高い性能を出しやすく解析も容易で、これで十分な問題も多数存在します。ここでは再び3次元デカルト座標系でのblast waveを例にとり、爆発の中心付近に高解像度の格子を配置してみましょう。(注：これはむしろ説明のための例で、あまり物理的に重要な設定ではありません。むしろ衝撃波を分解したいという場合が多いでしょうから、そのような場合にはAMRを使用することになります)。

SMRやAMRを用いるのに特別なコンフィギュレーションは必要ありませんが、SMRやAMRを用いると複雑なデータ構造を伴うのでHDF5を使うことを推奨します。以下のようにコードを設定してコンパイルしましょう。

> python configure.py --mpiccmd CC --prob blast -b -mpi -hdf5
> make clean ; make

では次にパラメータファイル（athinput.blast_smrとします）を編集して計算領域の中心部、爆発源の周囲に高解像度の格子を配置します。まず<mesh>ブロックの"refinement"パラメータを"static"に指定します。Athena++ではMeshBlockは並列化や計算の単位であると同時にリファインメントの単位でもあり、一つの格子を倍の解像度の格子8個（3次元の場合。2次元ならば4個、1次元ならば2個）に分割するようになっています。そのためSMR・AMRを利用するためにはまず計算領域を複数のMeshBlockに分割する必要があります。MeshBlockを小さくすればより柔軟な格子配置を実現できますが、一方で計算の効率は悪くなるため、性能と格子配置のバランスを取って格子サイズを適切に設定する必要があります。ここでは64³のルートグリッドを16³のMeshBlockに分割しています。

<mesh>
nx1        = 64         # Number of zones in X1-direction
x1min      = -1.0       # minimum value of X1
x1max      =  1.0       # maximum value of X1
ix1_bc     = periodic   # inner-X1 boundary flag
ox1_bc     = periodic   # outer-X1 boundary flag

nx2        = 64         # Number of zones in X2-direction
x2min      = -1.0       # minimum value of X2
x2max      =  1.0       # maximum value of X2
ix2_bc     = periodic   # inner-X2 boundary flag
ox2_bc     = periodic   # outer-X2 boundary flag

nx3        = 64         # Number of zones in X3-direction
x3min      = -1.0       # minimum value of X3
x3max      =  1.0       # maximum value of X3
ix3_bc     = periodic   # inner-X3 boundary flag
ox3_bc     = periodic   # outer-X3 boundary flag

num_threads = 1         # Number of OpenMP threads per process
refinement = static

<meshblock>
nx1    =  16
nx2    =  16
nx3    =  16

次に高解像度で分解する領域を指定します。以下では中心の半分の領域に二倍の分解能の格子を配置します。

<refinement1>
x1min=-0.5
x1max= 0.5
x2min=-0.5
x2max= 0.5
x3min=-0.5
x3max= 0.5
level=1

<refinement2>、<refinement3>というように複数の領域を指定することも可能です。これらの領域が重複していても構いませんし、必要な領域を精密に指定する必要はありません。コードは指定された領域を必要な解像度で覆うために最小限の格子を自動的に生成します。より高解像度の格子が必要な場合は"level"パラメータを大きくします。1レベル増やすごとに解像度が2倍になります。

Athena++では粗い格子を細かく分解する際に粗い格子を保持しないため、MeshBlockを1個リファインすると7個のMeshBlockが生成されます。そのため必ずしもMeshBlockの数がきれいに割り切れる数になるとは限りませんが、効率よい計算をするための設計なので心配は不要です。ただしこのために、実際に計算する前に生成されたMeshBlockの数を把握する必要があります。そこでまずコードを"-m [nprocess]"オプションを付けて実行します。[nprocess]はプロセス数で、指定したプロセス数でのロードバランスを表示しますが、MeshBlockの数を計算するには1としておけば問題ありません。以下ではコマンドラインで実行した例を示していますが、XC50ではデバッグキューを用いて1コアで実行すると良いでしょう。

>~/athena/bin/athena -i athinput.example -m 1
Root grid = 4 x 4 x 4 MeshBlocks
Total number of MeshBlocks = 120
Number of physical refinement levels = 1
Number of logical  refinement levels = 3
  Physical level = 0 (logical level = 2): 56 MeshBlocks, cost = 56
  Physical level = 1 (logical level = 3): 64 MeshBlocks, cost = 64
Number of parallel ranks = 1
  Rank = 0: 120 MeshBlocks, cost = 120
Load Balancing:
  Minimum cost = 1, Maximum cost = 1, Average cost = 1

See the 'mesh_structure.dat' file for a complete list of MeshBlocks.
Use 'python ../vis/python/plot_mesh.py' or gnuplot to visualize the mesh structure.

この例ではレベル0（ルートグリッド）に56個、レベル1（解像度2倍）に64個、合計120個のMeshBlockが生成されました。-mオプションの後の数値を変更することで指定したプロセス数でロードバランスがどうなっているかを確認できます。可能な限り等分配になるように計算を実行するようにしましょう。また生成された格子の配置が"mesh_structure.dat"ファイルに出力されているので、これを可視化することで計算格子の配置を確認することができます。これにはgnuplotを用います：

> gnuplot
> splot "mesh_structure.dat" w l

Mesh Structure
中心部分に高解像度の格子が生成されていることが分かると思います。では、これを8並列で計算してみましょう。全部で120MeshBlockあるので、1プロセスあたり15MeshBlockを割り当てることになります。

#PBS -q [queue]
#PBS -l nodes=1
#PBS -l walltime=0:10:00
#PBS -N blast

cd ${PBS_O_WORKDIR}

aprun -n 8 -N 8 -S 4 ./athena -i athinput.blast_smr > log

VisItを用いた計算の可視化は今までと同じように行うことができます。VisItでは"Mesh"プロットを用いることで格子の配置を表示することができます。
IndexSelect1
そのままでは全てのセルが表示されるために非常に込み合った画像になってしまうので、可視化に使うデータ点を間引く"IndexSelect"オペレータを用い、IncrをMeshBlockのサイズと同じにすることでMeshBlockの境界を抽出します。（"Apply operators to all plots"のチェックが外れていることに注意してください。）
IndexSelect2
Blast Wave with SMR

Adaptive Mesh Refinement

計算結果に適合して動的に最適な格子を配置する技術がAdaptive Mesh Refinementです。この機能を利用するには格子を生成または消去する条件を指定する必要があり、Athena++ではこれはユーザー定義関数を用いることで実現されます。ユーザー定義関数については後でより詳細に説明するので、ここではとりあえずサンプルを実装してAMRとはどのようなものかを見ることにします。

格子の生成（リファインメント）・消滅（デリファインメント）の条件は様々で、こうすれば十分というレシピが確立されているわけではなく、問題に合わせて検討する必要があります。よく使われる条件としては密度や圧力などの物理量の勾配（1階微分）または曲率（2階微分）が一定値を超えた時にリファインすることで衝撃波などの構造を捕捉するものがあります。また星形成等の自己重力系ではジーンズ波長を最低4メッシュで分解必要があることが知られており、これを条件として用いるジーンズ条件も良く使われます。

ここでは例として2次元のKelvin-Helmholtz不安定性を考えてみましょう。これは速度差のある二層の境界面が不安定となって大きな渦が生じる不安定で、稀に等間隔に並んだ渦状の雲として観測されることで有名ですが、ジェットと星間物質の間の混合など宇宙物理学でも重要な応用があります。境界層を高解像度で分解する条件としては、速度シアの大きさを使うのが自然です。"src/pgen/kh.cpp"にはこれが既に実装しされていますので、このファイルを開いて眺めてみましょう。コードの先頭付近にリファインメント条件関数のプロトタイプが宣言されています。

int RefinementCondition(MeshBlock *pmb);

次に、AMRの条件にこの関数を用いることを指定するため、Mesh::InitUserMeshData関数内でEnrollUserRefinementConditionを呼び出します。また、パラメータファイルからthrという値を読み取って、これをthresholdという変数に保存しています。この関数は計算領域全体に関わる関数や変数を初期化するのに用いるもので、他にも追加の物理過程をセットする時などに使います（実際にここでは解析用の変数を設定するコードがありますが、とりあえずは無視してください）。

void Mesh::InitUserMeshData(ParameterInput *pin) {
  vflow = pin->GetReal("problem","vflow");
  iprob = pin->GetInteger("problem","iprob");

  if (adaptive) {
    threshold = pin->GetReal("problem", "thr");
    EnrollUserRefinementCondition(RefinementCondition);
  }
  if (iprob == 4 && NSCALARS > 0) {
    AllocateUserHistoryOutput(1);
    EnrollUserHistoryOutput(0, PassiveDyeEntropy, "tot-S");
  }
  return;
}

RefinementCondition関数の本体はコードの下部にあります。この関数は速度シアを測定して、各MeshBlock内の一点でも規定値thresholdを上回った場合はそのMeshBlockをリファインするようフラグを立て（正の整数値を返す）、逆に全ての点で規定値の半分を下回った場合はそのMeshBlockをデリファインしても安全であるというフラグ（負の整数値を返す）を立て、どちらでもない場合はそのまま何もしない（0を返す）という判定をする関数です。なお実際の格子の生成・消滅はそのMeshBlock内だけではなく周辺のMeshBlockとの整合性（隣接するブロックは同一レベルか1レベル差しか許されないという条件があります）を保ちつつ指定した条件を満たすように計算されるので、意図したものよりも多くの格子が生成されることがあることに注意してください。

int RefinementCondition(MeshBlock *pmb) {
  AthenaArray<Real> &w = pmb->phydro->w;
  Real vgmax = 0.0;
  for (int k=pmb->ks; k<=pmb->ke; k++) {
    for (int j=pmb->js-1; j<=pmb->je+1; j++) {
      for (int i=pmb->is-1; i<=pmb->ie+1; i++) {
        Real vgy = std::abs(w(IVY,k,j,i+1) - w(IVY,k,j,i-1))*0.5;
        Real vgx = std::abs(w(IVX,k,j+1,i) - w(IVX,k,j-1,i))*0.5;
        Real vg  = std::sqrt(vgx*vgx + vgy*vgy);
        if (vg > vgmax) vgmax = vg;
      }
    }
  }
  if (vgmax > threshold) return 1;
  if (vgmax < 0.5*threshold) return -1;
  return 0;
}

よく似た例が二重マッハ反射テスト問題 ("src/pgen/dmr.cpp")でも使われていますので参考にして下さい。内容を確認したら、コードをconfigureしてmakeし直します。

> python configure.py --mpiccmd CC --prob kh -mpi -hdf5
> make clean ; make

その次にパラメータファイルを見てみましょう。"inputs/hydro/athinput.kh_amr16"を開いてみてください。

<mesh>
nx1        = 256       # Number of zones in X1-direction
x1min      = -0.5      # minimum value of X1
x1max      =  0.5      # maximum value of X1
ix1_bc     = periodic  # inner-X1 boundary flag
ox1_bc     = periodic  # inner-X1 boundary flag

nx2        = 256       # Number of zones in X2-direction
x2min      = -0.5      # minimum value of X2
x2max      =  0.5      # maximum value of X2
ix2_bc     = periodic  # inner-X2 boundary flag
ox2_bc     = periodic  # inner-X2 boundary flag

nx3        = 1         # Number of zones in X3-direction
x3min      = -0.5      # minimum value of X3
x3max      = 0.5       # maximum value of X3
ix3_bc     = periodic  # inner-X3 boundary flag
ox3_bc     = periodic  # inner-X3 boundary flag

refinement  = adaptive # AMR
derefine_count = 10    # allow derefinement after 5 steps
numlevel    = 4        # number of AMR levels

<meshblock>
nx1        = 16        # Number of zones in X1-direction
nx2        = 16        # Number of zones in X2-direction
nx3        = 1         # Number of zones in X3-direction

<hydro>
iso_sound_speed = 1.0 
gamma           = 1.4       # gamma = C_p/C_v

<problem>
iprob = 5
amp   = 0.01
thr   = 0.005
drat  = 2.0
vflow = 0.5
b0    = 0.1
sigma = 0.2
a     = 0.01

この例ではリファインメントのレベル数をルートレベルを含め最大4階層の格子を生成するよう指定しています（numlevels=4）。言い換えれば、ルートレベルに対して最大2³ = 8倍高い解像度の格子までが生成されます。"derefine_count"パラメータはMeshBlockが不必要に生成・消滅を繰り返すのを防ぐために、生成されたMeshBlockがすぐに消去されないように保持するステップ数（この例では10）を指定します。最後の<problem>ブロックでこの問題に関するパラメータを幾つか設定しています。このうちthrが上で見た、AMRの格子を生成・破壊するための基準値になっています。MeshBlockのサイズとして、この例では16²セルを指定しています。繰り返しますがMeshBlockのサイズは性能と柔軟性のバランスで決めるべきものですので、本計算を始める前に自分の問題に適した分解能や分解条件を事前に試しておく必要があります。この計算はかなり時間のかかる計算なのでMPIを利用して並列実行しましょう。ジョブスクリプトは以下のようにすると良いでしょう。この例では2ノード80コアを使用しています。

#PBS -q [queue]
#PBS -l nodes=2
#PBS -l walltime=0:30:00
#PBS -N kh

cd ${PBS_O_WORKDIR}
aprun -n 80 ./athena -i athinput.kh_amr16 > log

計算を実行すると、コードはまず初期条件を生成した後にリファインメントの条件をチェックし、初期条件がAMR条件を満たすように高解像度の格子を収束するまで繰り返し生成します。これにより初期の設定よりも大幅に多いMeshBlockが生成され、次のような警告が表示されることがあります。これが意図したものである場合は無視しても構いません。

### Warning in Mesh::Initialize
The number of MeshBlocks increased more than twice during initialization.
More computing power than you expected may be required.

AMRデータの可視化

AMRのデータもSMRの場合と同様にVisItで可視化することができます。
Blast Wave with AMR

なおこれは磁場なしの例ですが、"inputs/mhd/athinput.kh_amr16"に磁場有版のパラメータファイルも用意してありますので、余裕があれば試してみてください。

AMR計算の２段階初期化

この項は講習会では扱わないので必要に応じて試してください。AMRでは、生成された格子が多すぎて計算コストがかかり過ぎる場合があります。一度は初期条件を計算しないと生成されるMeshBlockの数を把握することはできないので、このような場合は「2段階初期化」を行います。まず適当なコア数で初期条件を生成して一度リスタートファイルに保存し（リスタートファイルの生成手順は３次元計算のページを参照してください）、そこから必要なコア数を指定して計算を再開します。以下のようにtime/nlim=0をコマンドラインに付加することで初期条件だけを生成します。

> aprun -n 16 ~/athena/bin/athena -i athinput.blast time/nlim=0 > log

次に-mオプションを利用して生成されたMeshBlock数を確認し、使用するコア数を決定します。ただし、AMRではMeshBlockの数が減ることもあり、MeshBlockの数が並列数を下回ると実行できなくなるため、1コア当たり数個のMeshBlockを割り当てるよう並列数を設定します。

> aprun -n 1 ~/athena/bin/athena -r amrexample.out1.00000.rst -m 64

そしてこのリスタートファイルから計算を再開します。この時time/nlimに本来の計算ステップのリミットを（最初の一回だけでいいので）指定して下さい。

> aprun -n 64 ~/athena/bin/athena -r amrexample.out1.00000.rst time/nlim=100000 > log