OpenMHD code
![[screenshot]](images/openmhd-banner13.png)
Japanese | English
OpenMHD is a multidimensional finite-volume code for magnetohydrodynamics (MHD). The code is written in modern Fortran. It is parallelized by using MPI-3.1, OpenMP, and CUDA.
OpenMHD was originally developed for our studies on magnetic reconnection in space physics [1,2]. Over the past 10 years, many improvements have been made to the code. The code has been made publicly available in the hope that others may find it useful.
The following version is available in .tar.gz format.
The source code is hosted on GitHub.
OpenMHD solves the following equations of magnetohydrodynamics. \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} where $p_T=p+B^2/2$ is the total pressure, $\overleftrightarrow{I}$ is the unit tensor, $e=p/(\Gamma-1) + \rho v^2/2 + B^2/2$ is the energy density, $\Gamma=5/3$ is the adiabatic index, and $\psi$ is a virtual potential for hyperbolic divergence cleaning [3]. The second-order Runge=Kutta methods, the second-order MUSCL scheme, and the HLLD flux solver [4] are employed. The source term for $\psi$ is handled by an operator splitting method and an analytic solution $\psi = \psi_0 \exp [ - ({c_h^2}/{c_p^2}) t ]$. Other numerical techniques are documented in Refs. [1,2] and references therein.
This page explains how to use OpenMHD.
Users with a Google Colaboratory account can try running calculations and analyzing data in the cloud.
Download the tar.gz file from the OpenMHD homepage.
Extracting the downloaded file with a tar command will create many directories full of related files.
(The following operations require a Linux/macOS terminal environment.
Windows does not handle symbolic links within the UNIX directory structure.)
$ gunzip openmhd-20240130.tar.gz $ tar xvf openmhd-20240130.tar
Move the files to an appropriate directory.
To get started, navigate to the 2D_basic directory.
$ cd openmhd-20240130/2D_basic/
OpenMHD compiles source code by using the make command.
First, use a text editor to open the file Makefile and change the variable F90 to set the compiler.
For MPI parallel computing,
specify an MPI-specific Fortran compiler.
For serial computing, simply specify an ordinary Fortran compiler.
F90 = gfortran
Then use the make command to compile the program.
$ make
MPI parallel code is included with OpenMHD 2D problems. The default setting is to compile both serial and parallel versions. If you only want to compile the serial version, just using
$ make run
is fine as well.
Please refer to the contents of Makefile for the rules for this process.
After successful compilation, you should have a serial executable (a.out) and
an MPI parallel executable (ap.out).
Next, run the program. For the serial version, specify the executable directly on the command line.
$ ./a.out
You can also use
$ ./a.out | tee output.log
The directory you are in (current directory) is represented as "." in UNIX, so instead of writing "a.out", write "./a.out" to indicate that the a.out file in the current directory is to be executed. The tee command displays the program's (standard) output on the screen and records it in a file (output.log).
When using the MPI parallel version, specify the number of processes
by editing the variable mpi_nums,
found around the top of the main file (e.g. mainp.f90).
Use (/4, 1/) to divide the computation area into 4 x 1, and
(/2, 2/) to divide it into 2 x 2. The (2) in mpi_nums(2) means 2-dimensional; there is no need to change it.
integer, parameter :: mpi_nums(2) = (/4, 1/) ! MPI numbers
After using the make command to compile the program,
use the mpirun command to run the executable.
Here, the -np option can be used to set the number of parallel processes to 4.
Make sure the number of processes here matches the number of processes specified in the mpi_nums variable.
$ mpirun -np 4 ./ap.out $ mpirun -np 4 ./ap.out | tee output.log
When the calculation is complete, the results are output into the data/ directory.
$ ls data/ field-00000.dat field-00007.dat field-00014.dat field-00021.dat field-00028.dat field-00035.dat field-00001.dat field-00008.dat field-00015.dat field-00022.dat field-00029.dat field-00036.dat field-00002.dat field-00009.dat field-00016.dat field-00023.dat field-00030.dat field-00037.dat field-00003.dat field-00010.dat field-00017.dat field-00024.dat field-00031.dat field-00038.dat field-00004.dat field-00011.dat field-00018.dat field-00025.dat field-00032.dat field-00039.dat field-00005.dat field-00012.dat field-00019.dat field-00026.dat field-00033.dat field-00040.dat field-00006.dat field-00013.dat field-00020.dat field-00027.dat field-00034.dat
On supercomputers,
commands such as qsub, qstat, etc., are often used to
submit jobs to compute nodes.
For guidance on using the job submission system,
please refer to the user manual for your supercomputer system.
When the calculation is successfully completed, plot the results and check them. OpenMHD provides Python (Matplotlib) and IDL visualization routines for every physics problem. Python supports two methods of analysis and visualization: IPython and Jupyter.
While there are a variety of methods for analysis and visualization in Python, we recommend using IPython.
To begin, open a working IPython shell with the ipython3 command.
A prompt in the form of In [number] will then be displayed.
$ ipython3 --pylab
Python 3.12.7 (v3.12.7:0b05ead877f, Sep 30 2024, 23:18:00) [Clang 13.0.0 (clang-1300.0.29.30)]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.29.0 -- An enhanced Interactive Python. Type '?' for help.
Using matplotlib backend: MacOSX
In [1]:
Next, try executing a sample script (plot.py).
%run -i is an IPython command; enter it starting from "%".
The -i option means files will be executed in interactive mode (see: IPython tutorial).
In [2]: %run -i plot.py
Once you execute this, you should see a plot similar to the one below.
If the plot does not appear, use plt.show() to display it explicitly.
In [3]: plt.show()
Those interested in the method for making plots can find it described in plot.py in the Python language.
Here, at the beginning of the file,
import openmhd
it loads the library for OpenMHD data files.
This library is actually located via a symbolic link
in visualization/python/openmhd.py and
it provides the following two subroutines.
data_readdata_read_from_bigendian
Then the data file (field/data-00020.dat) is specified and read in at the following location.
x, y are arrays of position coordinates,
t is time, and data is an array of physical quantities.
and
data is an array of physical quantities.
# reading the data ...
x,y,t,data = openmhd.data_read("data/field-00020.dat")
The range of data to be read can be refined by optional arguments
such as ix1, ix2, jx1, jx2 (array index numbers) and xrange, yrange (range of coordinates).
# reading the data (partial domain: [ix1,ix2] x [jx1,jx2])
x,y,t,data = openmhd.data_read("data/field-00020.dat",ix1=0,ix2=100,jx1=11)
x,y,t,data = openmhd.data_read("data/field-00020.dat",xrange=[0.0,3.14159],yrange=[0.0,3.14159])
If endianness differs between
the execution environment (supercomputer) and the analysis environment (PC, etc.),
it can be converted in Python.
When reading Big endian data in a Little endian environment,
use the data_read_from_bigendian routine.
Besides endian conversion, its operation is the same as data_read.
The data (data[ix,jx,9]) is read into a 3-dimensional array of size ix, jx, 9,
where the last dimension represents physical quantities.
For example, plasma pressure is in data[:,:,3] but
you can use a user-friendly form of data[:,:,pr] using the dummy variable at the beginning of the file.
The final dummy variable ps is the virtual potential $\psi$ for
the $\nabla \cdot \vec{B}$ problem (See: "Introduction > Basic Equations")
You can continue to work with these data and plots in the IPython shell. For example, the following will change the color bar from red to blue.
In [4]: myimg.set_cmap('RdBu_r')
Useful commands are left in the comments of the sample files; please try them out.
You can output multiple images and combine them into a single movie.
For 2D vortex problems (2D_basic),
a sample file (plot_movie.py) is provided for movie creation.
This sample creates sequentially numbered files.
Note that a directory named "movie" will be automatically created.
In [5]: %run -i plot_movie.py
For other physics problems, please prepare scripts, refering this sample.
# ---- loop for movies ---------------
for ii in range(0,41):
# reading the data ...
x,y,t,data = openmhd.data_read(ii)
Note that range (0,41) in Python will indicate 0,1,...40.
Check the movie/ directory once execution is complete.
$ ls movie/ output-00000.png output-00011.png output-00022.png output-00033.png output-00001.png output-00012.png output-00023.png output-00034.png output-00002.png output-00013.png output-00024.png output-00035.png output-00003.png output-00014.png output-00025.png output-00036.png output-00004.png output-00015.png output-00026.png output-00037.png output-00005.png output-00016.png output-00027.png output-00038.png output-00006.png output-00017.png output-00028.png output-00039.png output-00007.png output-00018.png output-00029.png output-00040.png output-00008.png output-00019.png output-00030.png output-00009.png output-00020.png output-00031.png output-00010.png output-00021.png output-00032.png
Then combine the separate image files into a single movie file using one of the following four methods.
It is possible to perform analysis in Python using Jupyter notebooks. If Jupyter is not installed, you can install jupyter and ipywidgets by executing the following command on the command line. (The ipywidgets is optional; it is only required for Jupyter's interactive features)
$ pip3 install jupyter ipywidgets [ $ pip3 install -U jupyter ipywidgets ]
Running the following command on the command line will allow you to see a Jupyter notebook in your browser screen.
$ jupyter-notebook plot.ipynb
Then by selecting a cell and pressing Shift + Enter,
you can execute the Python program in that cell.
The image plot will be displayed inline in the notebook.
The Python program is essentially the same as those introduced in Analysis and Visualization (Python/IPython). The 2D vortex problem (2D_basic) notebook also contains sample scripts for creating movies.
Last, press Ctrl + C in the terminal where the jupyter-notebook command was executed to terminate Jupyter.
After that, restart Jupyter, select the appropriate cell in your notebook, and execute the contents with Shift + Enter. For some physics problems, plot contents and color maps can be interactively changed.
We can also use Interactive Data Language (IDL) for analysis and visualization.
Start up IDL with the idl command.
IDL> is the IDL prompt.
$ idl IDL 8.7.2 (darwin x86_64 m64). (c) 2019, Harris Geospatial Solutions, Inc. Licensed for use by: Kobe University - Zenitani License: ***-******* IDL>
Execute the batch file (batch.pro) using the .r command as follows.
IDL> .r batch
You should see a plot like the following.
Methods for creating plots are described in batch.pro in the IDL language.
Anyone interested should take a look.
The OpenMHD-specific routine is loaded in a part below the top of the file.
The subroutine file is located at visualization/idl/data_read.pro via symbolic link.
;; find and compile data_read routine resolve_routine, "data_read"
Then load an appropriate data file (field/data-00020.dat).
x, y are arrays of position coordinates,
t is time, and data is an array of physical quantities.
;; reading the data ... data_read,data,x,y,t,'data/field-00020.dat'
The range of data to be read can be contolled by optional arguments such as ix1, ix2, jx1, jx2 (array index numbers) and xrange, yrange (range of coordinates).
;; reading the data (partial domain: [ix1,ix2] x [jx1,jx2]) data_read,data,x,y,t,'data/field-00020.dat',ix1=0,ix2=100,jx1=11 data_read,data,x,y,t,'data/field-00020.dat',xrange=[0.0,3.14159],yrange=[0.0,3.14159])
If the endianness is different between the execution environment (supercomputer) and the analysis environment (PC, etc.),
add the /swap_endian keyword to the openr function in the data_read.pro file to convert endianness when reading data.
; no record marker
openr,/get_lun,unit,filename,/swap_endian
The data is read into a three dimensional array of size ix, jx, 9,
in which the last dimension represents physical quantities.
For example, plasma pressure is in data[*,*,3] but
you can use a user-friendly form of data[*,*,pr] using the dummy variable at the beginning of the file.
When making a movie, refer to the sample file (batch_movie.pro) to output a large number of images.
IDL> .r batch_movie
When you are finished executing this file, the image files can be found in the movie/ directory, just like in Python.
$ ls movie/ output-00000.png output-00011.png output-00022.png output-00033.png output-00001.png output-00012.png output-00023.png output-00034.png output-00002.png output-00013.png output-00024.png output-00035.png output-00003.png output-00014.png output-00025.png output-00036.png output-00004.png output-00015.png output-00026.png output-00037.png output-00005.png output-00016.png output-00027.png output-00038.png output-00006.png output-00017.png output-00028.png output-00039.png output-00007.png output-00018.png output-00029.png output-00040.png output-00008.png output-00019.png output-00030.png output-00009.png output-00020.png output-00031.png output-00010.png output-00021.png output-00032.png
Then combine the separate image files into a single movie using one of the following methods.
OpenMHD results should be able to be visualized in languages and tools other than Python/IDL. Here we provide information that may prove useful in the development of visualization routines.
visualization/gnuplot in the OpenMHD distribution file.
The Python image processing library Pillow converts the image files
into a single animated GIF file.
To begin, install Pillow with pip.
$ pip3 install pillow [or $ pip install pillow]
Then run the following script in IPython, etc. Note that the name of the package to import is not Pillow but PIL.
import glob
from PIL import Image
images = []
for png in sorted(glob.glob("movie/output-*.png")):
im = Image.open(png)
images.append(im)
images[0].save('animation_pillow.gif', save_all=True, append_images=images[1:], optimize=False, duration=100, loop=0)
The duration option of the last save() method specifies the display time (in milliseconds) of each frame. By changing this value, you can adjust the playback speed of the video.
Reference
Using the Python interface to the computer vision library OpenCV allows you to convert separate images into an MP4 movie. There are a variety of installation methods. For example, you can install opencv-python with pip.
$ pip3 install opencv-python [or $ pip install opencv-python]
Then execute the following script in IPython, etc.
This sample should be included in the notebook file (plot.ipynb).
The package name of OpenCV is cv2.
import glob
import cv2
images = []
for png in sorted(glob.glob("movie/output-*.png")):
im = cv2.imread(png)
height, width, layers = im.shape
images.append(im)
out = cv2.VideoWriter("animation_cv2.mp4", cv2.VideoWriter_fourcc(*'MP4V'), 10.0, (width, height))
for i in range(len(images)):
out.write(images[i])
out.release()
The second argument of cv2.VideoWriter is the codec, the third is the video playback rate (fps), and the fourth is the video size (width and height).
Changing the fps value will alter the video playback speed.
Reference
ImageMagick can combine multiple image files into one.
The following command outputs the image files in the movie/ directory generated by the above 2D problems to an animated GIF.
$ convert movie/output-000*.png result.gif
The format of the output movie is automatically determined by the file name extension. For example,
$ convert movie/output-000*.png result.mp4
will generate an MPEG4 movie. MPEG4 conversion function calls FFmpeg, described below, so it also requires FFmpeg.
These tools can be installed on Ubuntu Linux as shown below.
$ sudo apt install imagemagick $ sudo apt install imagemagick ffmpeg
Reference
FFmpeg can also combine multiple image files into one.
The following command outputs the image files in the movie/ directory to an MPEG4 movie.
Compared to a GIF animation, the file size is considerably smaller. (However, the image quality is somewhat lower.)
$ ffmpeg -i movie/output-%05d.png result.mp4
The playback speed of the movie can be adjusted by the -r option.
$ ffmpeg -r 10 -i movie/output-%05d.png -r 10 result.mp4
FFmpeg can be installed on Ubuntu Linux in the following way.
$ sudo apt install ffmpeg
References
Sample Python 3 (matplotlib) and IDL routines are provided. For Python 3, please run the 'plot.py' script in the following way.
$ ipython3 --pylab
In [1]: %run -i plot.py
For IDL, please run the 'batch' script in the following way.
$ idl IDL> .r batch
The user can configure
simulation options in the main file (main.f90, mainp.f90) and
the initial configuration in the model file (model.f90, modelp.f90 ...).
The electric resistivity for magnetic reconnection can be configured
in the set_eta.f90 file.
Please edit the mpi_nums variable
in a main file (mainp.f90 or mainp.cuf).
integer, parameter :: mpi_nums(2) = (/4, 1/) ! MPI numbers
The user needs to set up boundary conditions
for both conserved variables (U) and MUSCL-interporlated primitive variables (VL and VR).
One might want to modify the boundary conditions
in the bc.f90 file for a serial version or
in the mpibc.f90 file for a parallel version.
For the parallel version, please look at the bc_periodicity parameters in the main file (mainp.f90).
In the .true. direction, the parallel communication module automatically transfers the data.
The user does not need to write a periodic BC code in the mpibc.f90 file.
In the .false. direction,
the user needs to prepare a boundary procedure in the mpibc.f90 file.
logical, parameter :: bc_periodicity(2) = (/.true., .false./)
Probably the allreduce communication fails when updating the timestep Δt.
We empirically see this error in the combination of macOS + Intel fortran + MPI library.
There is a workaround for this problem. Please contact us.
Set the n_start parameter to non-zero.
No, it is fixed to 1:1 (Δx = Δy).
We use Lorentz-Heaviside-style units, because they are convenient for calculations.
According to our benchmark, the HLLD flux solver is the most expensive part.
We encourage you to cite Phys. Plasmas 22, 032114 (2015) paper, because much of the code and documentation was written as a byproduct of this work.
For URL, please refer the github URL (https://github.com/zenitani/OpenMHD) or the entry (https://ascl.net/1604.001) for Astrophysics Source Code Library (ASCL). This URL (https://sci.nao.ac.jp/MEMBER/zenitani/openmhd-e.html) may not be stable.
Here are snippets for your convenience.
S. Zenitani, Phys. Plasmas, 22, 032114 (2015). S. Zenitani, Astrophysics Source Code Library, record ascl:1604.001 (2016).
Zenitani, S. 2015, PhPl, 22, 032114 Zenitani, S. 2016, OpenMHD: Godunov-type code for ideal/resistive magnetohydrodynamics, Astrophysics Source Code Library, ascl:1604.001
Please check whether the program uses OpenMP or not. Some compilers automatically activate OpenMP, even though the user does not specify so.
The following problems are pre-configured.
(1D_basic)(1D_relativisticJet)(2D_basic)(2D_KH)(2D_reconnection)(3D_basic)A Riemann problem is an initial value problem, in which one sets up two different states across a boundary (x=0.0) to see the time evolution of the system. The system evolves in a self-similar manner, accompanied by a complex structure. Since OpenMHD employs a Riemann solver that considers the Riemann problem at the cell boundaries, the Riemann problem is the most important problem to check the accuracy of the code. The initial configuration file (model.f90) contains parameters for typical test problems. The compound wave problem discussed in Brio & Wu (1988) is shown below.
Some problems are calculated with specific heat ratio $\Gamma=2$. If necessary, change the variable gamma in the parameter file (param.h). Also, make sure to change it back to $\Gamma=5/3$ later, because all the projects share param.h.
In active galactic nuclei (AGN) environments, jets are believed to emanate from the central object at relativistic velocities. Relativistic jets are accompanied by exotic phenomena. For example, Aloy & Rezzolla (2006) reported a new type of fluid acceleration in the boundary layer around the jet. This may sound counter-intuitive, but it is an adiabatic-type acceleration under the relativistic perfect fluid approximation. Zenitani et al. (2010) confirmed this phenomenon with 1-D relativistic MHD simulations and derived simple scaling laws.
The Orszag & Tang (1979) vortex problem is one of the most popular problems in 2-D. This problem is fun to watch, because many shocks and complex patterns are generated.
This is also an ideal problem to see the importance of the solenoidal condition ($\nabla\cdot\vec{B}=0$). To demonstrate this, let us comment out line 102 of the main.f90 file. If we set the parameter "ch" to zero, the solenoidal condition is no longer corrected. Your program will crash at some point, once unphysical errors are accumulated.
The Kelvin-Helmoholtz (KH) instability is a hydrodynamic instability that occurs in the velocity shear. OpenMHD has a KH vortex problem in a magnetohydrodynamic fluid, discussed in Matsumoto & Hoshino (2004). In the initial condition file (model.f90), we use the following settings: We consider a velocity shear in the vertical (y) direction. Fluids travel in the left on the upper side and to the right on the lower side. We also consider a density gradient by controlling a density in the lower side α. When α=1, the density is uniform, and when α=0.2, the upper and lower densities are five dimes different. The magnetic field is perpendicular to the 2-D plane.
\begin{align} \vec{V}(y) &= -0.5 V_0 \tanh(y) \vec{e}_x\\ N(y) &= \frac{N_0}{2} ( (1+\alpha) + (1-\alpha) \tanh(y) ) \\ \vec{B}(y) &= B_0 \vec{e}_z \end{align}The below figure shows the results for a resolution of 480 [X] x 600 [Y], and α=0.2. Secondary structures (K-H or R-T instability) grow at the boundary inside the KH vortex.
The magnetohydrodynamic KH instability is sensitive to the direction of the magnetic field. For example, let us change the initial conditions as follows. You will see that the in-plane magnetic field slows down the instability.
U(i,j,bx) = 0.d0
U(i,j,bz) = 0.d0
U(i,j,bx) = 0.2d0
U(i,j,bz) = 0.d0
Magnetic reconnection occurs in plasmas, when the oppositely-directed magnetic field lines are met each other. Magnetic reconnection requires the resistive MHD approximation, which incorporates electric resistivity in the Ohm's law. OpenMHD includes an improved version of Zenitani & Miyoshi (2011, 2015)'s code. The initial condition is as follows.
\begin{align} \vec{B}(y) &= B_0 \tanh(y/L) {\hat{x}},\\ \vec{v} &= 0,\\ \rho(y) &= \rho_0 \big[ 1 + \cosh^{-2}(y/L)/ \beta_{\rm in} \big], \\ p(y) &= p_0 \big[ \beta_{\rm in} + \cosh^{-2}(y/L) \big], \end{align}This $\beta_{\rm in}$ is the plasma-to-magnetic pressure ratio in the inflow region. For simplicity, the temperature ($T=p/\rho$) is set uniform. We further impose (1) an enhanced electric resistivity and (2) a weak magnetic perturbation around the origin, in order to trigger reconnection there.
\begin{align} \eta(x,y) &= \eta_0 + (\eta_1-\eta_0)\cosh^{-2} ( \sqrt{x^2+y^2} ) \\ \delta A_z &= 0.03 \times 2LB_0 \exp[-(x^2+y^2)/(2L)^2 ] \end{align}
The above figure shows the horizontal plasma velocity ($V_x$), calculated with default settings. One can see a fast plasma jet from the reconnection point. The jet velocity is approximated by the Alfven velocity in the inflow region. When the jet speed exceeds the sound speed, the reconnection jet generates various shocks and discontinuities in the downstream. These structures are schematically illustrated in the "plasmoid diagram" (a review article [in Japanese]).
This test problem advects a tilted magnetic loop. The initial configuration is as follows. Here, A is the vector potential and the subscripts 1,2,3 indicate the components in tilted coordinates. The magnetic field B is calculated from $\vec{B} = \nabla \times \vec{A}$. For more detail, see the Gardiner & Stone (2008) paper.
\begin{align} \rho &= 1, \\ p &= 1, \\ \vec{v} &= (1,1,2) \\ A_3(r) &= 0.001 ~ {\rm max}( 0.3 - r, 0 ) \\ r &\equiv \sqrt{x^2_1+x_2^2} \end{align}The following panel shows the magnetic energy density, after a certain timestep. The loops are essentially advected by the ambient plasma flow.
ifort, mpiifort to ifx, mpiifort -fc=ifc.flux_resistive_uni).flux_solver) now contains the hyperbolic-divergence-cleaning solver (flux_glm).xrange, yrange are added to the data_read routine.flux_resistive solver in the 2020-08-30 CPU code. (Courtesy of Mr. Nashita)glm_ss2, flux_resistive subroutines.ipywidgets for all physics problems. (Courtesy of Dr. Wada)data_read) to load the data. Now one can use a file path as well as a file number.n_start parameter is negative, the program searches an appropriate restart file.As of May 2024, the following papers have been published by OpenMHD users.