1次元流体・磁気流体シミュレーション
流体力学の数値シミュレーションは非常に複雑で様々な要素がありますが、まずは習うより慣れよという言葉に従ってコードを動かしてみましょう。ここでは数値計算コードのテストにも良く使われる衝撃波管問題を扱います。これは一言で言えば「ある仕切りを挟んで左側と右側に異なる流体を置いて、仕切りを取り払ったらどのように時間進化するか」という問題です。解析解が存在するのでテスト問題として有用なだけでなく、現代的な流体力学シミュレーションで使われるRiemann Solverの基礎となっている重要な問題です。ここではこの問題を題材にコードの動かし方を学びます。
この箱で表示されているものは計算機での実行例になります。 > ← この記号（プロンプト）で始まっている行はコマンドを表します。
0. コードのアップロードとディレクトリの移動
まず（もしまだであれば）コードのファイルを全てXC50にアップロードしましょう。初期設定が済んでいない人は環境設定のページを参照して準備をしてから次に進んでください。以下ではコードを/home/[username]/athena（以下このディレクトリをコードのルートディレクトリと呼びます。[username]の部分を各自のユーザー名に読み替えて下さい）に展開したとして説明します。ログインしたらコードのルートディレクトリに移動します。
> pwd /home/[username] > cd athena > pwd /home/[username]/athena
Sodの衝撃波管問題
最初に磁場のない流体力学での衝撃波管問題の典型例としてSodの衝撃波管を計算してみます。
1a. コードの設定（コンフィギュレーション） - XC50の場合
Athena++は様々な計算に対応しており、行う計算に応じて必要な機能（オプション）を選択します。これらにはコンフィギュレーション時に指定するものと実行時に指定するものがありますが、まずはコードのコンフィギュレーションを行います。これはルートディレクトリ直下のconfigure.pyというPythonスクリプトで行い、様々なオプションを指定することで使用する機能を選択します。まずどのようなオプションがあるのか、以下のようにして表示してみてください。
> python configure.py -h
usage: configure.py [-h]
[--prob {blast,cpaw,default_pgen,disk,dmr,eos_test,fft,field_loop,field_loop_poles,from_array,gr_blast,
gr_bondi,gr_geodesic_infall,gr_linear_wave,gr_mhd_inflow,gr_shock_tube,gr_torus,hb3,hgb,jeans,jet,jgg,kh,linear_
wave,lw_implode,magnoh,mignone_advection,noh,orszag_tang,poisson,resist,rotor,rt,scalar_diff,shk_cloud,shock_tube,
shu_osher,slotted_cylinder,ssheet,strat,turb,twoibw,visc}]
[--coord {cartesian,cylindrical,spherical_polar,minkowski,sinusoidal,tilted,schwarzschild,kerr-schild,gr_user}]
[--eos {adiabatic,isothermal,general/eos_table,general/hydrogen,general/ideal}]
[--flux {default,hlle,hllc,hlld,roe,llf}]
[--nghost NGHOST] [--nscalars NSCALARS] [-b] [-sts] [-s]
[-g] [-t] [-debug] [-coverage] [-float] [-mpi] [-omp]
[--grav {none,fft}] [-fft] [--fftw_path FFTW_PATH] [-hdf5]
[-h5double] [--hdf5_path HDF5_PATH]
[--cxx {g++,g++-simd,icpc,icpc-debug,icpc-phi,cray,bgxlc++,clang++,clang++-simd,clang++-apple}]
[--ccmd CCMD] [--mpiccmd MPICCMD] [--gcovcmd GCOVCMD]
[--cflag CFLAG] [--include INCLUDE] [--lib_path LIB_PATH]
[--lib LIB]
...（長いので省略）...
これらのうち最低限必要なのは初期条件を設定するproblem generatorを設定する`--prob'オプションで、それ以外は問題や環境に応じて設定します。詳細についてはマニュアルのコンフィギュレーションの項を参照して頂くことにして、ここではまず最も典型的なSodの衝撃波管問題を解いてみます。以下のようにコンフィギュレーションスクリプトを実行してください。
> python configure.py --prob shock_tube -mpi --mpiccmd CC Your Athena++ distribution has now been configured with the following options: Problem generator: shock_tube Coordinate system: cartesian Equation of state: adiabatic Riemann solver: hllc Magnetic fields: OFF 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
"--prob"オプションでsrc/pgen/ディレクトリから初期条件を設定するためのソースコードを選択し、"-mpi"でMPIによる並列化を有効にしています（ただしここでは実際には１コアで実行します）。性能的にはIntelコンパイラがおすすめです（`--cxx icpc'を指定します）が、講習会ではコンパイル時間の節約のためにg++を使用します（特別な指定は不要です）。Cray XC50上ではコンパイラに依らずコンパイルコマンドは`CC'ですので、`--mpiccmd CC'でコンパイラコマンドを指定します。
1b. コードの設定 - 手元のg++の場合
手元の計算機で非並列（シリアル）で実行したい場合、g++を利用する場合は次のように指定します。"--prob"オプションで実行する問題設定を切り替えます（src/pgen/以下の初期条件生成ファイルを切り替えています）。g++はデフォルトなのでコンパイラを指定するオプションは不要です。また、ここでは説明しませんが、Mac上で実行する場合は標準コンパイラになっているclang++を指定するのが良いでしょう。
> python configure.py --prob shock_tube Your Athena++ distribution has now been configured with the following options: Problem generator: shock_tube Coordinate system: cartesian Equation of state: adiabatic Riemann solver: hllc Magnetic fields: OFF 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: OFF OpenMP parallelism: OFF FFT: OFF HDF5 output: OFF Compiler: g++ Compilation command: g++ -O3 -std=c++11
2. コードのコンパイル
前項でコンパイルに必要なMakefileが生成されたので、これを用いてコードをコンパイルします。まず、コンフィギュレーションを変更した後は念のため一時ファイルを消去します。
> make clean
そしてmakeします。（-jオプションを付けると並列でコンパイルされるのでコンパイル時間を短縮できますが、それでもこのコンパイルには2-3分程度時間がかかります。講習会ではログインノードの負荷を下げるため-jはつけないでコンパイルしてください。)
> make
これで実行ファイル`athena'がbin/ディレクトリに生成されます。なおXC50でIntelコンパイラの場合に表示される"icpc: command line warning #10121: overriding '-xCORE-AVX512' with '-xhost'"という警告は無視して構いません。また、cray-c++のエラーが出る場合にはPrgEnvが正しく指定されていません。環境設定を確認して、PrgEnvを切り替えてください。
3a. 計算用のディレクトリ作成 - XC50上の場合
出力されるファイルを整理するために、１計算ごとに新しいディレクトリを作ることをお勧めします。CfCAのXC50では大容量のデータは/work以下に格納することになっていますので、ここに計算用ディレクトリを作成して移動しましょう。
> cd /work/[username] > mkdir sod > cd sod
3b. 計算用のディレクトリ作成 - 自分のPCの場合
出力されるファイルを整理するために、１計算ごとに新しいディレクトリを作ることをお勧めします。ホームディレクトリ以下にworkディレクトリを作成し、その下にそれぞれ問題ごとに計算用ディレクトリを作成することにします。以下の作業はこのディレクトリで行います。
> mkdir ~/work > cd ~/work > mkdir sod > cd sod
4. 実行ファイルと計算パラメータファイルのコピー
実行ファイル"athena"を計算用ディレクトリにコピーします。更に、Sodの衝撃波管問題のパラメータファイルのサンプルが/home/[username]/athena/inputs/hydro/athinput.sodにありますので、これも計算用ディレクトリにコピーします。
> cp ~/athena/bin/athena . > cp ~/athena/inputs/hydro/athinput.sod .
5. ジョブ投入用ファイルの作成 - XC50等ジョブ管理システムのある場合
XC50を初めとする大型計算機では、多数のユーザーに計算資源を効率的かつ公平に提供するために、ジョブ管理システムによってCPU資源の割り当てが行われます。計算を実行するためには「計算に必要なノード（計算機）の数」や「計算時間」等の必要な情報と、実行のためのコマンドを指定したスクリプトファイルが必要になります。（手元の計算機で実行する場合はこの項目は不要なので飛ばしてください。）XC50ではPBSというシステムが使われていますが、機関によっては別のシステムが使われていることがあり、その場合はそれぞれ異なる書式が必要になります。ここではjob.shというファイルで講習会用のキューに1CPUで計算を実行するスクリプトを作成します。一般の場合の詳細についてはCfCA提供のマニュアルを参照してください。job.shというファイルを作成して以下の内容を入力して下さい（手元の計算機で作成して転送しても構いません）。[queue]には環境設定のページを参照して、利用するキュー名で置き換えてください。
#PBS -q [queue]
#PBS -l nodes=1
#PBS -l walltime=0:10:00
#PBS -N sod
cd ${PBS_O_WORKDIR}
aprun -n 1 ./athena -i athinput.sod > log
例えば講習会期間初日では以下のようにします。
#PBS -q R6968428
#PBS -l nodes=1
#PBS -l walltime=0:10:00
#PBS -N sod
cd ${PBS_O_WORKDIR}
aprun -n 1 ./athena -i athinput.sod > log
まず最初の#PBSから始まる行でジョブ管理システムに対して要求するリソースに関する設定をします。"-q"はジョブを投入するキューを指定します。"-l nodes"では確保するノード数を指定し、ここでは1ノードでの実行を指定しています。XC50では1ノードは40コアです。"-l walltime"では最大の実行時間を指定しますが、ここでは短時間で終了するので10分としています。"-N"ではジョブの名前を指定できますが、自分がわかりやすいものであれば何でも構いません。
後半は実行するコマンドで、まず"cd ${PBS_O_WORKDIR}"で実行時ディレクトリ（この場合ジョブを投入したディレクトリ）に移動します。"aprun"はXC50上での計算実行するためのコマンド（他の並列計算機におけるmpiexec/mpirunに対応）で、"-n 1"で1プロセス（非並列）での実行することを示します。"athena -i [input file]"で計算に使用するパラメータファイルを指定します。"> log"は標準出力（通常画面に出力される内容）をlogというファイルに出力（リダイレクト）するという意味です。
6a. 計算の実行 - XC50の場合
では計算を実行しましょう。qsubコマンドを実行してジョブスクリプトをキューに投入します。
> qsub job.sh
計算の実行状況を確認するにはqstatコマンドを使います。
> qstat ...（略）... 1234567.sdb sod [username] 00:00:01 R debug 1234568.sdb sod [username] 0 Q debug
STATの欄が"Q"であれば実行待ち状態、"R"であれば実行中であることを意味します。計算が終了すればリストから消えます。長いリストから自分のジョブだけを表示させるには"qstat -u [username]"とします。計算の途中経過を見るにはtailコマンドでlogファイルを表示してみましょう。
> tail log
計算が終了したらlogファイルを確認してみましょう。
> less log RootGrid = 1 x 1 x 1 MeshBlock 0, rank = 0, lx1 = 0, lx2 = 0, lx3 = 0, level = 0 is=2 ie=257 x1min=-0.5 x1max=0.5 js=0 je=0 x2min=-0.5 x2max=0.5 ks=0 ke=0 x3min=-0.5 x3max=0.5 Setup complete, entering main loop... cycle=0 time=0.00000000000000e+00 dt=1.32055352301331e-03 cycle=1 time=1.32055352301331e-03 dt=8.96163836433927e-04 ... cycle=350 time=2.50000000000000e-01 dt=7.11985791492906e-04 Terminating on time limit time=2.50000000000000e-01 cycle=350 tlim=2.50000000000000e-01 nlim=-1 cpu time used = 5.00000007450581e-02 zone-cycles/cpu_second = 1.79200000000000e+06
6b. 計算の実行 - 手元の計算機で実行する場合
手元の計算機で、シリアル実行する場合はもっと簡単です。単純にパラメータファイルを指定して実行すればOKです。
> ./athena -i athinput.sod RootGrid = 1 x 1 x 1 MeshBlock 0, rank = 0, lx1 = 0, lx2 = 0, lx3 = 0, level = 0 is=2 ie=257 x1min=-0.5 x1max=0.5 js=0 je=0 x2min=-0.5 x2max=0.5 ks=0 ke=0 x3min=-0.5 x3max=0.5 Setup complete, entering main loop... cycle=0 time=0.00000000000000e+00 dt=1.32055352301331e-03 cycle=1 time=1.32055352301331e-03 dt=8.96163836433927e-04 ...
バックグラウンドで実行する場合は"> log &"等を付けて実行すると良いでしょう。
7. 計算結果の表示
計算終了後には次のようなファイルが計算用ディレクトリ内に出力されているはずです。
> ls athinput.sod Sod.block0.out1.00006.tab Sod.block0.out1.00013.tab Sod.block0.out1.00000.tab Sod.block0.out1.00007.tab Sod.block0.out1.00014.tab Sod.block0.out1.00001.tab Sod.block0.out1.00008.tab Sod.block0.out1.00015.tab Sod.block0.out1.00002.tab Sod.block0.out1.00009.tab Sod.block0.out1.00016.tab Sod.block0.out1.00003.tab Sod.block0.out1.00010.tab Sod.block0.out1.00017.tab Sod.block0.out1.00004.tab Sod.block0.out1.00011.tab Sod.block0.out1.00018.tab Sod.block0.out1.00005.tab Sod.block0.out1.00012.tab Sod.block0.out1.00019.tab
.tabファイルは単なるテキストでデータを列挙したいファイルでdt=0.01ごとに出力されています。.hstファイルは計算の履歴を格納するファイルですが、ここではとりあえず置いておいて、ここでは.tabファイルをgnuplotで表示してみましょう。データを手元の計算機に転送して（X Windowを転送すればXC50上でもできますが、ログインノードに負荷をかけるのは好ましくないので）、gnuplotを起動します。データの2列目がx座標、3列目が密度なので以下のようにプロットします。
> gnuplot gnuplot> plot "Sod.block0.out1.00025.tab" using 2:3 with lines
このような画像が表示できれば成功です。一番右の不連続が衝撃波、次の不連続面が接触不連続面、滑らかに密度が変化している領域が膨張波と呼ばれる構造です。後はファイルの番号を変えて時間進化を見たり、他の変数（tabファイルのヘッダ部に説明があります）を表示して、初めてのシミュレーション結果を堪能して下さい。
Brio-Wuの衝撃波管問題
次に、磁場のある場合の例としてBrio & Wuの衝撃波管問題を計算します。
8a. コードの再設定とコンパイル - XC50の場合
まずコードのルートディレクトリに戻り、磁場を有効にする"-b"オプションを付けて再度コンフィギュレーションを行います。以下のようにコマンドを実行してください。）
> cd ~/athena > python configure.py --prob shock_tube -mpi --mpiccmd CC -b Your Athena++ distribution has now been configured with the following options: Problem generator: shock_tube 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
コードのコンフィギュレーションを行った後は必ず"make clean"で一時ファイルを消去するのを忘れないでください。もし不安なら、"make"する前は必ず"make clean"としても構いません（少し余計に時間がかかるだけです）。
8b. コードの再設定とコンパイル - 手元の計算機の場合
手元の計算機でも同様です。"--prob shock_tube"に加えて磁場を有効にする"-b"オプションを付けて"configure.py"を実行します。
> cd ~/athena > python configure.py --prob shock_tube -b ... > make clean > make
コードのコンフィギュレーションを行った後は必ずmake cleanで一時ファイルを消去するのを忘れないでください。
9. 計算の準備
同様に、計算データを出力するためのディレクトリを作成して実行ファイルとBrio-Wu衝撃波管用のパラメータファイルをコピーし、（XC50上で実行する場合は）ジョブ投入用スクリプトファイルを作成します。
> cd /work/[username] > mkdir bw > cd bw > cp ~/athena/bin/athena . > cp ~/athena/inputs/mhd/athinput.bw .
#PBS -q [queue]
#PBS -l nodes=1
#PBS -l walltime=0:10:00
#PBS -N bw
cd ${PBS_O_WORKDIR}
aprun -n 1 ./athena -i athinput.bw > log
10. 計算の実行、結果の表示
では計算を実行して、終了したら結果を表示してみましょう。
gnuplot> plot "Brio-Wu.block0.out1.00040.tab" using 2:3 with lines
このような結果が得られたでしょうか。磁気流体計算では磁場無しの場合と比べ磁場と関係した様々な波が現れ衝撃波管問題の解もより複雑になりますが、ここではその詳細は教科書等に譲ります。
追加課題
これで１次元計算の基本は終わりですが、余裕のある人はコードに慣れるためにももう少し色々試してみましょう。
1. Riemann Solverの変更
Athena++では数値流束の計算にはHLLCが、磁気流体力学の計算ではHLLDが標準のRiemann Solverとして選択されます。これらはどちらも計算コストに対して比較的精度が高く安定なので標準となっていますが、問題によっては他のSolverを利用することが良い場合もあります。Athena++では流体力学ではHLLE、HLLCとRoe法が、磁気流体力学ではHLLE、HLLDとRoe法がサポートされているので、これらを試してみて下さい。Riemann Solverを変更するにはコンフィギュレーションで--fluxオプションを指定します。指定できるオプションについては"python configure.py -h"で表示されるヘルプメッセージを参照してください。
> python configure.py --prob shock_tube --flux roe > python configure.py --prob shock_tube -b --flux hlle
Riemann Solverを変更すると解はどのように振る舞うでしょうか？HLLE法は接触不連続面やAlfven波を分解しないために波の構造が鈍ってしまいます。Roe法は波の方程式を線形化して解くもので、全ての波を分解できるので精度は良いですが強い衝撃波や膨張波などで不安定になりやすいという欠点があります。HLLD法はMiyoshi & Kusano 2005で提案された高性能かつ高精度で安定性も高い、現代の宇宙物理学における磁気流体計算の事実上の標準ともいえる解法です。非常に丁寧に書かれた教育的な論文ですので、是非一度論文を読んでみることをお薦めします。特殊相対論や等温状態方程式への拡張も存在します。
2. その他の衝撃波管問題
衝撃波管問題は初期条件の左右の状態を変えることで多様な振る舞いをします。Athena++では流体力学のテスト問題としてinputs/hydro/athinput.einfeldt1125とathinput.einfeldt1203、磁気流体力学ではinputs/mhd/athinput.rj2aが用意してありますのでこれらを試してみましょう。更に余裕があれば（パラメータファイルの変更は後に扱いますが）パラメータファイルの最下部にある<problem>以下のパラメータが左右の初期状態を表していますので、これらを変更して遊んでみるのも良いでしょう。衝撃波管問題を系統的に調べた論文としてRyu and Jones 1995, ApJ, 442, 228がありますので参考にして下さい。
3. 高次精度・特性変数による補間
Athena++ではデフォルトの2次精度Piecewise Linear Methodの他に3次精度Picewise Parabolic Methodもサポートしています。これを用いると、高次精度の補間を行うため滑らかな解をより精度良く計算することができるようになります。これを有効にするにはまず"--nghost 4"としてコードをconfigureします。その後パラメータファイルの<time>ブロック内にある"xorder = 2"を"xorder = 3"に書き換えます。結果をプロットして比較してみましょう。
もう一つ精度の改善に寄与する（ことがある）のは補間変数の変更です。Athena++では標準ではprimitive variables（密度・速度・圧力）を用いて空間補間を行いますが、これを特性変数characteristic variablesに変更できます。これを有効にするには"xorder"を"2c"または"3c"に設定します。これにより変数変換に伴う計算コストは増えますが、特に強い衝撃波などでの解の振る舞いが改善します。