ここでは多次元シミュレーションとその可視化の方法を学ぶために、典型的な磁気流体シミュレーションのテスト問題であるOrszag-Tang Vortexを実行してみましょう。以下の例は講習会期間中にXC50で実行することを想定していますので、コンフィギュレーション等については自身の環境に合わせて適宜読み替えて下さい。

1. コードの設定とコンパイル
まずは1次元の場合と同様にコードのコンフィギュレーションとコンパイルを行います。多次元計算を行うのに特別な設定は必要ありませんが、Orszag-Tang Vortexの計算のためには"--prob orszag_tang"を指定します。また、磁気流体計算のために"-b"オプションを指定します。

> python configure.py --prob orszag_tang -b -mpi --mpiccmd CC
Your Athena++ distribution has now been configured with the following options:
  Problem generator:          orszag_tang
  Coordinate system:          cartesian
  Equation of state:          adiabatic
  Riemann solver:             hlld
  Magnetic fields:            ON
  Number of scalars:          0
  Special relativity:         OFF
  General relativity:         OFF
  Frame transformations:      OFF
  Self-Gravity:               OFF
  Super-Time-Stepping:        OFF
  Debug flags:                OFF
  Code coverage flags:        OFF
  Linker flags:
  Floating-point precision:   double
  Number of ghost cells:      2
  MPI parallelism:            ON
  OpenMP parallelism:         OFF
  FFT:                        OFF
  HDF5 output:                OFF
  Compiler:                   g++
  Compilation command:        CC  -O3 -std=c++11
> make clean
> make

2. 計算用ディレクトリの作成と実行ファイル・パラメータファイルのコピー
これも前回と同じです。

> cd /work/[username]
> mkdir orszag_tang
> cd orszag_tang
> cp ~/athena/bin/athena .
> cp ~/athena/inputs/mhd/athinput.orszag-tang .

3. 計算の実行
これも前回と同じです。以下のようなスクリプトファイルを用意しましょう。

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

cd ${PBS_O_WORKDIR}

aprun -n 1 ./athena -i athinput.orszag-tang > log

計算が終了すると大量のVTKファイルが生成されているはずです。これらを手元の計算機にダウンロードしてください。

VisItによる可視化

ではこのVTKファイルを可視化してみましょう（注：一部使用しているVisItのバージョンが古い画像がありますが、内容に相違はありません）。VTK（Visualization Tool Kit）フォーマットは様々なソフトウェアでサポートされているのでお好みのものがあればそれを用いてもらって構いません。ここではLawrence Livermore National Laboratoryによって開発されているVisItを利用する方法を紹介します。

1. VisItでファイルを開く
VisItを起動すると操作用と表示用の二つのウィンドウが開きます。まず操作用のウィンドウからOpenボタンを押してファイルを開きます。

VisItは自動的に連番のファイルをグループ化してくれます。"OrszagTang.block0.out2.*.vtk database"を選択して"OK"をクリックします。

2. Pseudocolorプロット
2次元の画像を可視化する一番簡単な方法は疑似カラーマップ（Pseudocolor）です。これを用いて密度の2次元画像を表示してみましょう。"Add"ボタンから"Pseudocolor"→"rho"の順に辿ってオブジェクトを追加します。

追加できたら"Draw"ボタンをクリックします。表示用ウィンドウに真っ青な四角が表示されますが焦る必要はありません。これは最初の出力時刻における密度分布、即ち初期条件が表示されています。 操作用ウィンドウのスライダーを操作するか、その下にある左右のボタンをクリックして時間を進めてみましょう。

初期条件から複雑な渦が生成されて行くのが見えれば成功です。

3. プロットを操作する
作成した画像をインタラクティブに操作することができます。表示用ウィンドウのギアのようなボタンをクリックしてから画像をドラッグすると左右に移動（3次元の場合は回転）できます。画像の一部分を拡大するには虫眼鏡のボタンをクリックしてから興味ある領域をドラッグするか、ホイールまたはピンチ操作を行います。 折れ線グラフのアイコンをクリックするとラインアウトモードになり、画像の上でドラッグするとその直線に沿った折れ線グラフを作成できます。

操作のし過ぎで訳が分からなったら、カメラに緑色の矢印が描かれているアイコン（下図の左から二番目）をクリックすると画像をリセットできます。

4. 表示する物理量の変更
表示する物理量を変更するには、操作用ウィンドウの"Variables"ボタン（ウィンドウサイズが小さいと"Draw"ボタンの右に隠れてしまうことがありますが、この場合は矢印をクリックするかウィンドウの幅を広げて下さい）をクリックするかまたは"Pseudocolor"のアイコンを右クリックしてメニューから"Variables"を選択し、表示したい変数を選択します。例えばガス圧の分布をみるために"press"変数を表示してみましょう。

プロットの色やレンジなどの設定を変更するにはPseudocolorアイコンをダブルクリックします。下の例では表示を0.06から0.6の対数プロットに変更し、更に黒→オレンジ→白と変化するカラーマップを選択しています。設定を変更した後は"Apply"ボタンをクリックして変更を適用します。


これで基本的な2次元のプロットは作成できるようになりました。よりカッコいい画像を作成するにはVisItのマニュアルやチュートリアルを参照してください。
なお、ゼロや負の値を含むデータを対数プロットしようとするとエラーになります。この場合は設定を変更した上で再度"Draw"すればOKです。

ソースコードを見てみよう

さて、ここまではほぼブラックボックスでシミュレーションを行ってきましたが、コードの中でどのように初期条件を生成しているのかを見てみましょう。src/pgen/orszag_tang.cppを適当なエディタで開いてみてください。詳しいことは後ほど説明しますので、ここでは初期条件の設定の雰囲気だけでもつかんでもらえればOKです。このファイルの中ではMeshBlock::ProblemGeneratorという関数が定義されており、この関数内で初期条件を設定します。Athena++では初期条件には各セルの流体力学の保存量（phydro->u）と（MHDの場合）各セル表面の磁場（pfield->b.x?f）を指定する必要があります。初期条件の計算に必要な座標はpcoord->x?v（セル中心の座標）やx?f（セル表面の座標）、dx?f（セル表面間の距離＝セルの長さ）等を用います。Athena++では多次元配列はC++言語のテンプレートを用いて実装しているため通常のC言語の配列とは異なりますが、その意味は直観的に理解できると思います。phydro->uの第一引数は物理量を表しており、IDN=密度、IM1=x1方向の運動量、IEN=全エネルギー等を表しています。以下の部分では密度と運動量の初期値を設定しています。

  for (int k=ks; k<=ke; k++) {
    for (int j=js; j<=je; j++) {
      for (int i=is; i<=ie; i++) {
        phydro->u(IDN,k,j,i) = d0;
        phydro->u(IM1,k,j,i) =  d0*v0*std::sin(TWO_PI*pcoord->x2v(j));
        phydro->u(IM2,k,j,i) = -d0*v0*std::sin(TWO_PI*pcoord->x1v(i));
        phydro->u(IM3,k,j,i) = 0.0;
      }
    }
  }

また、次の部分ではこれより上で計算したベクトルポテンシャル（のz成分）azを用いて、各セル表面での磁場を計算しています。

  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) = (az(j+1,i) - az(j,i))/pcoord->dx2f(j);
      }
    }
  }
  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) = (az(j,i) - az(j,i+1))/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) = 0.0;
      }
    }
  }

Pythonによる解析

VisItによる可視化は重要ですが、データをより定量的に解析したい場合はPythonなどのスクリプト言語を用いるのが便利です。この講習会ではPythonを用いた解析は扱いませんが、Athena++ではPythonに計算結果を読み込むための機能（/vis/python/以下のスクリプト）を提供しています。詳しくは公式WikiのPythonによるデータ読み込みのページを参照してください。

画像・動画の作成方法

VisItで作成した画像を保存するには操作用ウィンドウのFileメニューからSet save options...で保存サイズや保存先ディレクトリを指定した後、同じくFileメニューのSave windowで保存することができます。この際、作成した図をそのまま保存するためにはoptionの"Screen capture"にチェックを付けて実行します。またプロットのラベルや軸を操作したい場合にはControls→Annotationメニューからそれぞれ編集することが可能です。詳細はマニュアルを参照してください。

動画を作成する場合には同様にFileメニューのSave movieをクリックします。ウィザードに従って設定を入力すれば、VisItが自動的に一連の画像から動画を作成してくれます。このウィザードでは直接MPEG動画を作成することもできますが、エンコードエンジンの性能に問題があるのかあまり画質がよろしくないことがあるので、PNGなどで連番画像を作成したのちにお好みの動画エンコーダー（mencoderやffmpeg、TMPGEncなどで調べてください）を用いてパラメータを細かく調整してエンコードすることを個人的にはお勧めします。

追加課題

1. 断面図
Orszag-Tang Vortexは多次元磁気流体数値計算コードのテストとしてしばしば使われます。過去の文献を探して、自分のシミュレーションの結果と比較してみましょう。特に、Lineoutモードを使って折れ線グラフを描き、適当な断面図での分布を例えばAthenaコードの論文の図と比べてみましょう。

2. Rotor、Kelvin Helmholz不安定性
同じく多次元磁気流体計算のテストとして使われる問題にRotorテストやKelvin Helmholz不安定性があります。余裕があればこれらの問題も実行してみましょう。--prob rotor (kh)でコンフィギュレーションを行い、inputs/mhd/athinput.rotor (inputs/hydro/athinput.kh-shearまたはkh-slip)のパラメータファイルを使って計算することができます。

3. 論文
Athena++で採用しているアルゴリズムは基本的にはStone & Gardiner 2009のvan-Leer+CT法です。興味のある方は論文を読んでみて下さい。