Athena++ ver.21.0
Orbital Advection と Shearing Box の使い方
小野智弘（東京工業大学）

2020年度CfCA 流体学校

1

高い品質の結果を低い計算コストで得たい計算コストを抑えたい

良い結果が欲しい
解像度を上げる相反

解像度を下げる

追加の物理を考慮

簡単な物理だけ考える

良い問題設定計算手法の工夫◼ 研究課題を遂行するための問題設定を練り上げる◼ 明らかにしたい事柄から解像度や考慮する物理を決定する◼ 数値計算で使う計算手法を吟味する➔ 実際に数値計算を開始する
2020年度CfCA 流体学校

2

Athena++ ver.21.0 で使える『計算手法の工夫』
解像度を上げるために◼ ノード数を増やしより多くの計算リソースを使用する◼ Mesh Refinementを使用する (SMRでOK) ? SMR : AMR

局所系の使用◼ Shearing Box (円盤局所系)

高速化モジュールの使用◼ Super-Time-Stepping の使用 (拡散項が入った計算を高速化)◼ Orbital Advection (ケプラー回転円盤の計算を高速化)
2020年度CfCA 流体学校

3

Shearing Box (SB) 系
円盤の局所系

SB系での基礎方程式
𝜕𝜌
+∇∙
𝜕𝑡

𝜌𝒗 = 0𝜕𝐸+∇∙
𝜕𝑡𝐸𝒗 − 𝑩 ⋅ 𝒗 𝑩 = 2𝑞Ω20 𝑥𝑣𝑥𝑑𝒗
1= − ∇𝑃 + 2𝑞Ω20 𝑥 𝒊Ƹ
𝑑𝑡
𝜌

潮汐力

෡×𝒗
− 2Ω0 𝒌コリオリ力

shear rate

メリット

デメリット

• グローバル計算に比べ低コスト• 境界条件の性質が良い• 曲率項が無視されている• 大スケールの現象は追えない
2020年度CfCA 流体学校

4

SB系特有の境界条件
SB系では 𝑥 方向にシアー周期境界を用いる場合がある
𝑦

密度、x2方向以外の速度場、磁場に対して
𝑤 𝑥in , 𝑦, 𝑧 = 𝑤 𝑥out , 𝑦 + 𝑡𝛿𝑣, 𝑧𝑥in

0𝑥 𝑥out

x2方向の速度場に対して
𝑣𝑦 𝑥in , 𝑦, 𝑧 = 𝑣𝑦 𝑥out , 𝑦 + 𝑡𝛿𝑣, 𝑧 + 𝛿𝑣

エネルギーに対して
𝐸 𝑥in , 𝑦, 𝑧 = 𝐸 𝑥out , 𝑦 + 𝑡𝛿𝑣, 𝑧 + 𝜌𝛿𝑣𝑣𝑦 + 𝜌𝛿𝑣 2 /2

シアー速度 𝑣𝑦0 = −𝑞Ω0 𝑥
2020年度CfCA 流体学校

5

Shearing Box の使い方
Input ファイル内で① x1方向の境界条件を “shear_periodic” にする<mesh>
ix1_bc= shear_periodic #ox1_bc= shear_periodic #② シアリングボックスのパラメータを設定<orbital_advection>qshear= 値 # シアー率𝑞 デフォルト: 0, ケプラー回転: -1.5
Omega0= 値 # 系の角速度Ω0 デフォルト: 0
shboxcoord = 値 # 2次元x-z 平面を扱う時は2, それ以外は1

ただし、(qshear !=0.0 && Omega0 != 0.0) が必要

2020年度CfCA 流体学校

6Shearing Box 使用上の制限◼ Cartesian座標のみ◼ シアー周期境界はx1方向の境界条件としてのみ有効◼ シアーの流れは x2方向 (shboxcoord == 1)
or x3方向 (shboxcoord == 2 && 二次元)
◼ Time Integratorは 3次精度以下 (rk1, vl2, rk2, rk3)
◼ SB系でReconstructionの精度を4次以上にするメリットは少ない◼ x2方向に一様グリッド (x2rat=1.0)◼ x1方向の境界に面するMesh Blockについて以下の条件を満たすことが必要① シアー周期境界に面するMesh Blockはx2方向に同じレベル② シアー周期境界で面する可能性のあるMesh Blockは同じRefinementレベル◼ 自己重力との併用はまだ未実装（近いうちに）◼ その他はだいたい使えるはず
2020年度CfCA 流体学校

7

Orbital Advection (OA) アルゴリズム
FARGO (Fast Advection in Rotating Gaseous Objects) アルゴリズム とも呼ばれる
点源重力天体の周りにある円盤の回転速はおよそ Keplerian 𝑣𝐾 =𝐺𝑀/𝑅

𝑅 ⟶ 0 の時 𝑣𝐾 ⟶ ∞ なので 𝑣𝐾 ≫ 𝑐𝑠 になることがある➔ 速い Kepler 回転速度のために (復習 CFL条件 )円盤計算ではタイムステップdtが小さくなりがち➔ Kepler回転成分が既知・定常であることを利用しOperator Splitting法でKepler回転成分だけを分離して解く手法が OA◼ OAを使用すると dt が大きくなり、円盤計算の計算コストを下げることができる◼ 円盤の局所系である Shearing Box 系でも使用可能
2020年度CfCA 流体学校

8

軌道移流法でやってることഥ + 𝒗𝑲 という形に分ける（ただし、𝒗𝑲 は軌道速度場）
速度場を 𝒗 = 𝒗
ഥが速度場だと考えた時の保存変数を 𝒖ഥ とすると𝒗
基礎方程式
𝜕ഥ𝒖 𝜕𝑭𝟏 (ഥ𝒖) 𝜕𝑭𝟐 (ഥ𝒖) 𝜕𝑭𝟑 (ഥ𝒖)𝜕ഥ𝒖ഥ++++ 𝒗𝑲=𝜳 𝒖
𝜕𝑡
𝜕𝑥
𝜕𝑦𝜕𝑧
𝜕𝑦

赤は普通のMHD基礎方程式と同一、青が軌道速度による移流項
この式を Operator Splitting 法を使って𝜕ഥ𝒖
𝜕𝑭𝟏 (ഥ𝒖)𝜕𝑭𝟐 (ഥ𝒖)𝜕𝑭𝟑 (ഥ𝒖)ഥ+++=𝜳 𝒖
𝜕𝑡𝜕𝑥
𝜕𝑦𝜕𝑧
𝜕ഥ𝒖
𝜕ഥ
𝒖+ 𝒗𝑲 = 0𝜕𝑡𝜕𝑦(operator A) 流体パート(operator B) 移流パートという2つの処理 (A, B) に分けて順次解いていく
2020年度CfCA 流体学校

9それぞれのオペレータ処理をどうやって解くの？
• Operator A: 流体パート
𝜕ഥ𝒖 𝜕𝑭𝟏 (ഥ𝒖) 𝜕𝑭𝟐 (ഥ𝒖) 𝜕𝑭𝟑 (ഥ𝒖)ഥ+++=𝜳 𝒖𝜕𝑡𝜕𝑥𝜕𝑦𝜕𝑧
普通のMHD基礎方程式と同じなので、Riemann solver を用いて解く• Operator B: 軌道移流パート
𝜕ഥ𝒖𝜕ഥ
𝒖+ 𝒗𝑲=0𝜕𝑡
𝜕𝑦
軌道速度場 𝒗𝑲 は既知・定常➔ 解析解が存在するのでCFL条件を気にしなくて良い
Athena++ の OA では時間1次・2次精度の operator splitting に対応(STS の 時間2次精度と同様)
2020年度CfCA 流体学校

10

OAにおける時間精度の違い
OA 1次精度のメリット➢ OA使用による 1 cycle辺りの計算コスト増大が小さい (20%程度)➢ 移流を解く際に発生する数値拡散を軽減できる

OA 2次精度のメリット➢ 時間積分由来のエラーが卓越する場合でも、結果が2次精度になる➢ 他の operator splitting を使ったモジュール (ex. STS) と併用しても高い性能
使い分けの基準[上級]
流体パートと移流パートの非可換性が高い時は OA 2次精度

[初級]

STSと併用する時は必ずOA2次精度
まずOA2次精度を使い、OA1次精度でも大丈夫そうなら使う
2020年度CfCA 流体学校

11

OAの有効化
inputファイル内で<orbital_advection>
OAorder = 値 # 0:OAなし (デフォルト)、1: OA1次精度、2: OA2次精度[Cartesian座標]<orbital_advection> qshear, Omega0 の両方が0でないor
ユーザー定義軌道速度の設定✓ <orbital_advection> shboxcoord = 1✓ 軌道速度の方向はx2方向で2次元 or 3次元問題[cylindrical座標 & spherical_polar座標]<problem> gm が0でない or
ユーザー定義軌道速度の設定✓ 軌道移流の方向はx2方向で2次元 or 3次元問題 (--coord=cylindrical)✓ 軌道移流の方向はx3方向、3次元問題 (--coord=spherical_polar)
2020年度CfCA 流体学校

12

OAを使う上での注意点✓ 計算内で使う速度場は軌道速度場を差っ引いたもの
初期条件、境界条件を与える際に注意が必要✓ Time Integratorは 3次以下 (rk1, vl2, rk2, rk3)
✓ hydro・fields・passive scalar の全てに対応✓ Mesh Refinementとの併用が可能✓ 拡散物理を考慮した際も使用可能、STSとの併用も可✓ STSとの併用時、時間2次精度OAの使用を強く推奨✓ 相対論と一緒には使用不可✓ シアリングボックスで使う時は軌道移流法の使用を推奨

2020年度CfCA 流体学校

13

OA使用時の出力
Inputファイル<output#>
orbital_system = 値
値 = false の時
速度場 𝒗 が出力データ （デフォルト）
値 = true の時
ഥ = 𝒗𝑲 − 𝒗 が出力データ
速度場 𝒗
このオプションは file_type = rst 以外のoutputで利用可能

2020年度CfCA 流体学校

14

回転座標系
cylindrical or spherical_polar 座標使用時にinputファイル内で<orbital_advection>
Omega0=値 #➢ Omega0 の回転角速度を持った回転座標系 (コリオリ力 + 潮汐力)
➢ 回転方向はOAの軌道速度と同じ。

2020年度CfCA 流体学校

15

参考となるpgen•ssheet.cpp
shwaveテスト(e.g., Johnson & Gammie 2005)
2次元等温流体における Shearing Box テスト•jgg.cpp
MHD shwaveテスト(Johnson 2007)
3次元等温磁気流体におけるShearing Box テスト•hb3.cpp
2次元MRI問題 (Hawly & Balbus 1992)
2次元磁気回転不安定性の Shearing Box テスト•hgb.cpp
3次元MRI問題 (Hawly & Balbus 1992)
3次元磁気回転不安定性の Shearing Box テスト•disk.cpp

回転円盤のテスト

2020年度CfCA 流体学校

16