Athena++

3次元磁気流体シミュレーションと並列計算

さてここからはより実用的な計算問題へと進みましょう。3次元の計算には多大な計算コストがかかるために並列化が必要になります。Athena++はOpenMPによる共有メモリ並列化とMPIによる分散メモリ並列化に対応しており、両方を同時に使用することもできます。詳しい説明はここでは省きますが、共有メモリ並列化とは同じ物理的なメモリにアクセスできるノード内での並列化、分散メモリ並列化とはネットワーク越しのノード間での並列化（ノード内でも利用できます）と理解すればよいでしょう。ここではXC50のような大規模並列計算機で主に使用することになるMPIによる並列化にしぼって説明します。

MPIによる分散メモリ並列化

MPIは分散メモリ並列化の事実上の標準です。ここではCray XC50でMPIを用いた並列計算により、超新星爆発などのモデルとして使われる3次元のblast waveの計算を行います。また3次元データの可視化方法やHDF5を用いた並列出力についても説明します。では早速進めていきましょう。

1. コードをMPIとHDF5を有効にしてコンパイルする
MPIによる並列化には-mpiオプションを、HDF5による出力を有効にするには-hdf5オプションを付けてコンパイルします。XC50ではモジュールを読み込むだけで利用できます（環境設定参照）が、一般にはHDF5ライブラリ（MPIで利用する場合は並列版）をインストールする必要があります。HDF5による出力は必須ではありませんが、この後AMRやSMRを利用する際にはVTKによる出力は大量のファイルを生成してしまうので大規模並列には不便です。

> cd ~/athena
> python configure.py --mpiccmd CC --prob blast -b -mpi -hdf5
Your Athena++ distribution has now been configured with the following options:
  Problem generator:          blast
  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. パラメータファイルを編集し領域分割を指定する
Athena++では計算領域を（論理的に）同じ形状のMeshBlockと呼ばれる単位に分割することで並列化を行います。このMeshBlockは並列化の単位であると同時に、後程説明するSMRやAMRを用いる際の分割の単位でもあります。ではまず計算用ディレクトリを作成してblast waveのパラメータファイルを開いてみましょう。

> cd /work/[username]
> mkdir blast
> cd blast
> cp ~/athena/bin/athena .
> cp ~/athena/inputs/mhd/athinput.blast .

提供されているパラメータファイルでは64³の計算領域が指定されています。

<mesh>
nx1        = 64         # Number of zones in X1-direction
nx2        = 64         # Number of zones in X2-direction
nx3        = 64         # Number of zones in X3-direction

このままでも実行できますが、ここではこれを32³のMeshBlockに分割します。 "<mesh>" ブロックの後（本当はどこでも構わないのですが）に"<meshblock>ブロックを作成し、サイズを指定します。

<meshblock>
nx1        = 32         # Number of zones per MeshBlock in X1-direction
nx2        = 32         # Number of zones per MeshBlock in X2-direction
nx3        = 32         # Number of zones per MeshBlock in X3-direction

この例は計算領域を2×2×2 = 8のMeshBlockを生成します。当然のことですが、ルートグリッド（Mesh全体）のサイズはMeshBlockで割り切れる数でなければなりません。これをMPIプロセスに分配して並列計算を行うわけですが、各プロセスは1個または複数個のMeshBlockを持つことができます。またその数は等分配であることが望ましいですが、そうである必要もありません。この例では1～8プロセスを使用した並列計算が可能です。より大きい計算を試したい場合、Meshのサイズを拡大するのが良いでしょう。例えば

<mesh>
nx1        = 128        # Number of zones in X1-direction
nx2        = 128        # Number of zones in X2-direction
nx3        = 128        # Number of zones in X3-direction

<meshblock>
nx1        = 32         # Number of zones per MeshBlock in X1-direction
nx2        = 32         # Number of zones per MeshBlock in X2-direction
nx3        = 32         # Number of zones per MeshBlock in X3-direction

とすれば128³の領域を64個の32³のMeshBlockで分割することができます。なお、ここでの例はあくまでチュートリアルのためで、結果を早く得るために小さめの計算領域を設定しています。1プロセス当たりの計算量が多い方が計算時間は長くなりますが計算効率は上がります。問題や環境（そして忍耐力）にもよりますが、Cray XC50での単純な磁気流体計算では、コア当たり64³程度をお勧めします。適切な計算リソースの決め方については後程説明しますが、共有資産であるスーパーコンピューターを適切に利用することはユーザーの責任であることを忘れないでください。

3. MPIによる並列計算の実行
MPIの並列計算の実行方法は環境によって微妙に異なりますが、ここではCfCAのCray XC50の場合に限って説明します。他の計算機システムについてはそれぞれのマニュアルを参照してください。まずジョブを投入するためのスクリプトを作成します。例えば"blast.sh"というファイルを以下の内容で作成して下さい。

#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 -t 00:08:00 > log

スクリプトの内容は単一コアで実行する場合とほとんど変わりません。ここでは１ノード・８コアを用いて並列計算を行います。"#PBS -l nodes=1"で1ノード使用することを指定し、aprunのオプション"-n 8"で全体で８プロセス使用すること、"-N 8"でノードあたり８プロセス使用することを指定します。更に"-S 4"は、XC50は１ノード当たり２ソケット（CPUが2機）あるので、１ソケットあたり４プロセス使用することを意味します。基本的には"-n"だけ指定すれば十分ですが、各ノードの全コア（40コア）を使用しない場合はこれらの指定によりできる限り均等にプロセスを割り当てた方が性能が向上します。そのような場合は"-N"は全プロセス数/ノード数、"-S"は"-N"の半分になるように指定すれば良いでしょう。またathenaの"-t 00:08:00"オプションを用いて計算時間が8分間を超えた場合自動的に停止するよう設定しています。このファイルを用いて計算を実行します。

> qsub blast.sh

4. VTKファイルを結合する（機能として紹介しますが、飛ばして5に進んでも構いません。）
上記の設定のままで実行するとMeshBlockごとにVTKファイルが生成されるため大量のVTKファイルが生成されます。このVTKファイルを一つ一つ読み込むこともできますが、これらを結合するプログラムが"vis/vtk"以下の"join_vtk++"です。これを利用するには以下のようにコンパイルして、結合したいファイルを引数に実行します。（注：XC50でもログインノード上で実行する非並列コードのコンパイルにはg++が利用できます。）ただし、AMR・SMRを使用する場合はこの方法でファイルを結合することはできませんので、次項で説明するHDF5を利用してください。

> g++ -o join_vtk++ join_vtk++.cpp -lm

> ./join_vtk++ -o Blast.00000.vtk Blast.block0.out1.00000.vtk Blast.block1.out1.00000.vtk ...

5. HDF5による出力
VTKファイルでも結合すれば単一のファイルとして扱うことができますが、逐一面倒ですし結合スクリプトはAMRやSMRには対応していません。HDF5を利用すれば出力ごとにデータとVisItで読み込むためのXDMFファイルの合わせて２ファイルだけになりますし、AMRやSMRのデータをそのままVisItで扱うことができます。これを有効にするにはパラメータファイルを以下のように変更します。

<output1>
file_type  = hdf5       # HDF5 data dump
variable   = prim       # variables to be output
dt         = 0.1        # time increment between outputs

ghost_zones = true # enables ghost zone output

これを用いた可視化については以下で説明します。なお"ghost_zones = true"は境界条件に使用されているゴーストゾーンの情報も出力することを指定しています。これは必ずしも必要ではありませんが、VisItで磁力線を描く際にこれをonにしておかないとMeshBlockの切れ目で磁力線が途切れてしまうので有効にしておくことをお勧めします。

３次元シミュレーションの可視化

引き続きVisItを用いた可視化手法について説明します。

1. データを読み込む
前項で説明した通り、VisItを起動してOpenボタンを押してファイルを読み込みます。VTKファイルを読み込むには"*.vtk" databaseを、HDF5ファイルを読み込むには"*.athdf"ではなく"*.athdf.xdmf" databaseを読み込みます。このxdmfファイルは、VisItにどのようにHDF5を解釈するかを教えるためのものです。

2. Expressionsファイルを開く (HDF5のみ)
HDF5は全てのデータをスカラーとして保存しています。VisItのExpressionsという機能を用いることで、既存の変数から新しい変数を組み立てることができます。解析によく使われる幾つかの標準的な変数を"athena/vis/visit/athdf_MHD_primitive.xml"に用意されていますので、VisItの操作用ウィンドウのメニューから"Controls → Expressions"とクリックしてExpressionsウィンドウを開きます。"Load"ボタンを押して上記のファイルを読み込み、"Apply"してから"Dismiss"してウィンドウを閉じます。これで読み込んだ変数が"Variables"のリストに表示されます。

3. Pseudocolorプロットとスライスオペレータ
まずガス圧の分布を表示してみましょう。"Add → Pseudocolor → gas_pressure"の順にクリックしたのちに"Draw"ボタンを押します。おそらく真っ青な四角が表示されることと思いますが、これは計算領域の表面だけが表示されているからです。中身を見るためには色々方法がありますが、まずは断面図を使ってみましょう。VisItではプロットに対してデータを編集するオペレータを追加することができます。ここではx,y,zの断面図を表示する"ThreeSlice"オペレータを使います。"Operators → Slicing → ThreeSlice"と辿って、"Draw"をクリックします。表示用ウィンドウ内でドラッグすることで回転したり、操作用ウィンドウのスライダーを動かして時間進化を見てみましょう。爆発波が伝搬していく様子が見えるはずです。

各オペレータにもパラメータがあり、変更することでプロットを調整できます。"Pseudocolor"アイコンの左にある三角形をクリックするとオペレータのリストが展開されます。"ThreeSlice"アイコンをダブルクリックするとオペレータのパラメータを変更するウィンドウが開きます。この場合はスライスする平面を変更することができます。オペレータを解除するには×印を、またオペレータが複数適用されている場合にはアイコンの右にある上下の三角印で順序を変更することができます。

4. 速度ベクトルの表示
VisItでは様々なプロットを重ねて表示することができます。ここでは速度ベクトルを重ねて表示してみましょう。まず、標準の設定ではVisItは全てのプロットに対して同じオペレータを適用するようになっているので、"Apply operators to all plots"のチェックを外してこれを解除します。

次に"Add → Vectors → gas_velocity"としてVectorsプロットを追加します。このアイコンをダブルクリックしてパラメータ設定ウィンドウを開きます。そのままでは表示される矢印が少ないので、まずは数を増やしてみましょう。また、矢印が小さい場合には"Glyphs"タブ中のスケールを変更することでサイズを調整することができます。また、ベクトルの起点（vector origin）は"middle"にした方が自然なプロットになります。

設定をしたら"Apply"をクリックした後操作用ウィンドウの"Draw"をクリックします。パラメータを調整しながら、見やすいプロットを作ります。

5. 磁力線の描画
磁力線（あるいは速度場の流線）の構造は磁気流体シミュレーションにおいて重要な情報です。磁力線の描画には"Pseudocolor"に"IntegralCurve"オペレータを用います。"Add → Pseudocolor → operators → IntegralCurve → bcc（VTK）またはbfield（HDF5）"と辿ってプロットを追加したら、IntegralCurveオペレータをダブルクリックして設定ウィンドウを開きます。流線の描画には多数の設定が必要です。

流線はベクトル場を基準点から積分していくことで計算されます。そこでまず基準となる点の集まりを指定する必要があります。ここでは簡単に"Source"タブの"Source type"を"Plane"（平面）に指定します。更に"Normal"（法線ベクトル）を"1 0 0"（x軸方向）に指定します。この設定は問題によって調整して下さい。他のSource typeも試してみると良いでしょう。次に"Samples"をX,Yそれぞれ10に指定します。これでこの平面から10x10=100本の磁力線が計算されます。最後に、X,Yの"distance"を2.0にして計算領域全体が含まれるようにしましょう。

次に"Integration"タブの"Integration direction"（基準面から積分する方向）を"Both"にします。

次に"Appearance"タブに移動します。"Color by"を"Vector magnitude"に指定すると磁場の強度が色で表されます。あるいは"Solid"を指定すると単色で磁力線が表示されます。最後に"Apply"をクリックして設定を反映させます。流線の色については通常のPseudocolorの指定と同様に、Pseudocolorプロットをダブルクリックして出てくるウィンドウから指定できます。

操作用ウィンドウに戻って"Draw"をクリックしましょう。計算領域を斜めに横切る磁が表示されたと思います。この計算設定では比較的強い磁場を入れているために爆発によっても磁力線はあまり大きくは歪みません。磁力線に対してどちらの方向により早く情報（衝撃波）が伝播するか観察してみましょう。

このようにして表示された磁力線の時間進化の解釈には注意が必要です。なぜなら、IntegralCurveでプロットされた流線はその瞬間の基準点を通る磁力線を表示しているだけで、基準点は時間的に固定されているので流線の運動は実際の磁力線の運動とは対応していないからです。もし個々の磁力線の運動を追跡したいならば、基準点も流体の時間発展と整合的に動かしてやらなければなりません。

リスタートファイル

CfCAをはじめとして多くのスーパーコンピューターでは一回当たりの実行時間が制限されています。しかし実際の科学的研究には長時間の計算が必要なことが多く、何回もジョブを乗り継ぐ必要があります。Athena++には計算の途中過程をファイルに保存し、そのファイルから計算を継続する機能が用意されています。リスタートファイルを出力するためにはパラメータファイルに以下のようなoutputブロックを追加します。

<output2>
file_type   = rst      # restart dump
dt          = 1.0      # time increment between outputs

これによりシミュレーション内時間で1.0ごとに拡張子.rstのリスタートファイルを生成します。また、計算が制限時間などで正常に終了した場合、*.final.rstというファイルを生成します。このファイルを使用して計算を再開するには-iオプションの代わりに-rオプションでリスタートファイルを指定します

aprun -n 8 ./athena -r blast.final.rst -t 00:08:00 > log

追加課題

1. 計算パラメータの変更
この衝撃波の振る舞いを見るために、初期のパラメータを変更してみましょう。パラメータファイルの"<problem>"以下は個々の問題に固有のパラメータです。例えば爆発波の圧力や半径、磁力線の強度や方向を変更することができます。また磁場無しの計算もやってみましょう。"-b"なしでconfigureしなおして、"inputs/hydro/athinput.blast"を使って実行してみてください。

2. Rayleigh-Taylor不安定性
比較的計算時間がかかる問題なので全員の課題とはしませんが、Rayleigh-Taylor不安定性（高温低密度なガスの上に低温高密度なガスが成層している系）も流体・磁気流体シミュレーションの典型的な問題です。"--prob rt"でコンフィギュレーションを行いinputs/hydro/athinput.rt3dを用いて計算できますのでこちらも試してみて下さい。パラメータに応じて不安定の波長がどのように変わるでしょうか。またこの問題は<problem>b0パラメータを指定することで（磁場有のコードを使えば）磁場を入れることができますが、磁場により不安定性はどのように変わるでしょうか？