以前のAthenaと比べてAthena++の大きな変更点の一つは円筒座標や球面極座標などの直交曲線座標系に対応していることです（一般相対性理論もこの機能の延長上にあります）。この章では実際に円筒座標・球面極座標でのシミュレーションを行います。

円筒座標・球面極座標でのシミュレーション

デカルト座標以外での計算を行うには"--coord"オプションで座標系を指定してコンフィギュレーションを行います。対応している座標系のリストは"python configure.py -h"で表示することができます。例えば磁場無しのblast waveの計算を球面極座標系で実行するには以下のように指定します。

> python configure.py --mpiccmd CC --prob blast --coord spherical_polar -mpi -hdf5

"--coord spherical_polar"で球面極座標座標を指定しています。円筒座標にする場合には"--coord cylindrical"を指定してください。なお、blast.cppは円筒座標・球面座標にも対応していますが、Athena++が自動的に座標系を変換した初期条件を生成するわけではありません。実際の計算ではそれぞれの問題について使用する座標系に対応した初期条件を生成するコードが必要になります。

Athena++の円筒座標系ではx1=R, x2=ϕ（方位角）, x3=zになっており、球面極座標ではx1=r, x2=θ（極からの角度）, x3=ϕ（方位角）となっています。2次元計算を行う場合はx1, x2のみが使われますので、現時点ではAthena++では円筒座標系での2次元軸対称シミュレーションはサポートされていません。このような場合は球座標を用いるか、円筒座標でϕ方向に少しだけメッシュを取るか、あるいは独自にCoordinatesクラスを書く必要があります。

コードのコンパイルができたら球座標の場合は"inputs/hydro/athinput.blast_sph"、円筒座標の場合は"inputs/hydro/athinput.blast_cyl"を使って計算を実行します。この例では計算領域の中心から外れた点から球形のblast waveを発生します。このような座標系は等方ではないため、blast waveが球形を維持できるかどうかは自明ではありません。歪みが小さく保たれていることを確認しましょう。

データの解析

VisItで円筒座標や球面極座標のデータを可視化するのは簡単です。デカルト座標に変換するために"Transform"オペレータを利用します。以下では球座標の場合を例に取って説明します。

PseudocolorプロットにTransformオペレータを追加し、Pseudocolorプロットの左にある三角形をクリックしてオペレータのリストを展開します。Transformオペレータをダブルクリックして設定用ウィンドウを開き、まずCoordinateタブからInputに"Spherical"、Outputに"Cartesian"を選択します。いつも通り"Apply"してから"Draw"をクリックします。この時点では計算領域全体に対応する球（の一部）が表示されます。

次のThreeSliceオペレータを追加しましょう。そのままではスライスする平面がblast waveの中心に合っていないので、オペレータをダブルクリックしてスライスの位置を下図のように調整します。調整後は忘れずに"Apply"してから制御用ウィンドウの"Draw"をクリックします。


これで三次元球座標を用いたblast waveの計算の可視化ができました。blast waveが期待通り球対称性を保って進化することを確認しましょう。

追加課題

非一様格子間隔
多くのシミュレーションでは格子の間隔を一定に取りますが、特に円筒座標や極座標では格子間隔を一様にとると原点付近または外側で格子が大きく歪んでしまいます。このような状況は、計算精度に影響を及ぼすだけでなく短い格子間隔によって時間刻みが短くなるなどあまり好ましいものではありません。φ方向の格子間隔はrΔφに比例するので、動径方向の格子間隔も半径に比例するように取ることで各格子のアスペクト比を一定に近く保つことができます。Athena++では任意の格子間隔を設定することができますが、隣り合う格子の間隔が一定の比になるような格子であれば特に簡単に設定することができます。例えばパラメータファイルで次のように指定します。

<mesh>
nx1        = 64        # Number of zones in X1-direction
x1min      = 1.0       # minimum value of X1
x1max      = 3.0       # maximum value of X1
x1rat      = 1.017314  # ratio between neighboring cells

どのような比が適切かは計算領域のサイズとメッシュ数の両方によって決まります。動径方向の格子点数がN点の時、r^Nが計算領域の内側と外側の比に一致するように決める、つまりx1rat = (x1max/x1min)^(1/N)とするのが一つの目安ですが（上の例はそのようになっています）、問題に応じて適切に調節して下さい。また、任意の格子間隔を設定することも可能ですが、そのような複雑な格子を設定したい場合はマニュアルを参照してください。