OpenMHD コード

Japanese | English

OpenMHD は Fortran 90 で書かれた 磁気流体 (MHD) シミュレーションコードです。 シンプルな実装で動作が軽く、 Python/IDL の解析・可視化ルーチンを同梱しているため、 計算実行から可視化までを手軽に体験できるようになっています。 シリアル版、MPI/OpenMP並列版、GPU版、MPI/GPU版を同梱していますので、 PCでもスーパーコンピューターでも、CPU/GPUを問わず 同じコードを走らせることができます。

OpenMHD コードは、作者の MHD リコネクション研究 [1,2] の並列コードをパッケージ化したものです。 そして、その後の改良を取り入れて、随時更新を続けています。

ファイル

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

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

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

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

  1. コンパイラ
  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 \mathcal{E}}{\partial t} &+ \nabla \cdot \Big( (\mathcal{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}$ は単位テンソル、 $\mathcal{E}=p/(\Gamma-1) + \rho v^2/2 + B^2/2$ はエネルギー密度、 $\Gamma=5/3$ は断熱定数です。 これらの式は $4\pi$ や $\mu_0$ などの係数が消える Lorentz-Heaviside 風の表式を使っています。 音速は $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 シミュレーションでは非常によく使われています。 本稿執筆(2024年)時点では業界標準と言って良いでしょう。

コードの計算精度は時間・空間ともに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)
参考資料