OpenMHD コード

Japanese | English

OpenMHD は Fortran 90 で書かれた 磁気流体 (MHD) シミュレーションコードです。 シンプルな実装で動作が軽く、 IDL/Python3 の解析・可視化ルーチンも付属しているため、 シミュレーション入門に最適です。 シリアル版と MPI/OpenMP 並列版が同梱されていて、 スーパーコンピューターでも PC でも 同じコードを走らせることができます。

OpenMHD コードは、作者の MHD リコネクション研究 [1,2] の並列コードをパッケージ化したものです。 そして、その後の改良を取り入れて、随時更新を続けています。 MPI 通信コードの一部は 東大地惑・横山先生の CANSコード に由来しています。

ファイル

tar.gz ファイルを直接ダウンロードしてください。

開発中のコードは GitHub で管理しています。以下のように git コマンドを使うと、リポジトリの中身を複製することができます。

$ git clone https://github.com/zenitani/OpenMHD.git
動作環境

動作を確認している環境は以下の通りです。 あまり凝ったことをしていないので、 ほとんどの環境でそのまま動作するはずです。 MPI 並列計算には、MPI-2 規格に準拠した MPI ライブラリが必要です。

  1. コンパイラ
    • Intel Fortran v2013, v2015, v2016, v2017
    • gfortran 4.8, 4.9, 5.4, 6.1, 7.1
    • Cray Fortran コンパイラ
    • 富士通コンパイラ
  2. MPI ライブラリ(オプション)
  3. 可視化環境(オプション)
基礎方程式と計算手法

次のような保存形の磁気流体方程式を採用しています。 \begin{align} \frac{\partial \rho}{\partial t} &+ \nabla \cdot ( \rho \vec{v} ) = 0, \\ \frac{\partial \rho \vec{v}}{\partial t} &+ \nabla \cdot ( \rho\vec{v}\vec{v} + p_T\overleftrightarrow{I} - \vec{B}\vec{B} ) = 0, \\ \frac{\partial e}{\partial t} &+ \nabla \cdot \Big( (e+p_T )\vec{v} - (\vec{v}\cdot\vec{B}) \vec{B} + \eta \vec{j} \times \vec{B} \Big) = 0, \\ \frac{\partial \vec{B}}{\partial t} &+ \nabla \cdot ( \vec{v}\vec{B} - \vec{B}\vec{v} ) + \nabla \times (\eta \vec{j}) + \nabla \psi = 0, \\ \frac{\partial \psi}{\partial t} &+ c_h^2 \nabla \cdot \vec{B} = - \Big(\frac{c_h^2}{c_p^2}\Big) \psi, \end{align} このうち、$p_T=p+B^2/2$ は磁気圧を含めた総圧、 $\overleftrightarrow{I}$ は単位テンソル、 $e=p/(\Gamma-1) + \rho v^2/2 + B^2/2$ はエネルギー密度、 $\Gamma=5/3$ は断熱定数です。 ここでは $4\pi$ や $\mu_0$ などの係数が消える単位系を採用していますから、 音速は $c_s \equiv \sqrt{\Gamma p / \rho }$、アルヴェン速度は $c_A \equiv B / \sqrt{\rho}$ になります。

$\psi$ は計算結果を補正するための仮想ポテンシャルです。 本来、磁場はソレノイダル条件 $\nabla \cdot \vec{B} = 0$ を満たしますが、 シミュレーションでは数値エラーで $\nabla \cdot \vec{B} \ne 0$ となることがあります。 そうすると、磁気トポロジーが変わるだけでなく 磁気単極子に対して非物理的な力が働いて計算結果を壊してしまいます。 そこで、仮想ポテンシャル $\psi$ を介して $\nabla \cdot \vec{B}$ の数値エラーを 減衰・拡散させます(Hyperbolic divergence cleaning 法 [3])。 $c_h,c_p$ は減衰・拡散度合いを調整するパラメーターです。

OpenMHD コードは、数値流束を評価する際に、 隣接するセル境界での一次元リーマン問題を考える近似リーマン解法を使っています。 その中でも5種類の不連続面に対応した HLLD 法 [4] を採用しています。 HLLD 法は、リーマン問題を適度な計算コストで精度良く近似するため、 MHD シミュレーションでは非常によく使われています。 本稿執筆(2017年)時点では業界標準と言って良いでしょう。

コードの計算精度は時間・空間ともに2次精度です。 2段の Runge=Kutta 法あるいは TVD Runge=Kutta 法で時間発展を解くとともに、 物理量(基本量)を MUSCL 補間して2次の空間精度を達成しています。 ψのソース項についてのみ、Runge=Kutta の前後で 0.5Δt ずつ、 解析解 $\psi = \psi_0 \exp [ - ({c_h^2}/{c_p^2}) t ]$ を使っています。 他の技術的な詳細は、文献 [1,2] 及びその参考文献を参照してください。 OpenMHD を取り巻く研究状況については、「生存圏研究」の解説記事 [5] も参考になると思います。

参考文献
  1. S. Zenitani & T. Miyoshi, Phys. Plasmas 18, 022105 (2011)
  2. S. Zenitani, Phys. Plasmas 22, 032114 (2015)
  3. A. Dedner et al., J. Comput. Phys. 175, 645 (2002)
  4. T. Miyoshi & K. Kusano, J. Comput. Phys. 208, 315 (2005)
  5. 銭谷誠司、生存圏研究 13, 27 (2017)
参考資料