HARM2D Tutorial by Sasha Tchekhovskoy

How to download HARM2D, choose the problem, compile, and run the code:

You can download a version of HARM2D that I modified for this school (by adding several additional physical problems) by following this link:

Original HARM2D can be downloaded here .

To install the code, you do:

unzip -d harm2d harm2d.zip
cd harm2d
make clean

If you would like to use ICC instead of GCC, then edit makefile and replace "USEICC = 0" with "USEICC = 1".

Before running one of the problems with HARM2D, you choose the problem you want to run in the file decs.h.

You do so by modifying decs.h. There are 6 problems to choose from:

//which problem

#define BZ_MONOPOLE_2D 3
#define BONDI_PROBLEM_1D 5
#define BONDI_PROBLEM_2D 6

The default choice is


This is a magnetized torus accretion problem in 2D. For starters, however, it might be easier to start with a 1D problem without magnetic field, e.g., with BONDI_PROBLEM_1D.

To run the code, do the following:

make clean

• Compile:

make clean

• Run:


As the code runs, it produces sequential dump files in the 'dumps' subdirectory, i.e., dumps/dump###. It also produces dumps/gdump file that contains the information about the metric and the grid.

How to read in the output into Python

You can use any visualization program you are used to (Python, Gnuplot, IDL, Matlab, etc)
• Here are the scripts for reading in the data in Python: harm_script.py
• For a great Python tutorial, check out http://www.learnpython.org/

To use these Python scripts, do:

Start Interactive Python Shell:

ipython --pylab

You should do this from the location of the harm executable.

• Initialize the scripts:

%run -i harm_script.py

Read in grid information:


Read in dump file:


where you can replace "000" with any other number of a dump file (provided it was already written out).

Once you read in a dump file, you can print its current time via

print t

Understanding the grid

HARM2D is capable of using non-uniform grids. The examples below will automatically take care of converting from the grid coordinates, which are x1, x2, x3, to the physical coordinates, which are r, theta, phi.

How to visualize the output in Python

• Python Visualization Examples

For a 1D test: scroll down to the bottom of harm_script.py and replace

if False:
#1D plot example


if True:
#1D plot example

For a 2D test: scroll down to the bottom of harm_script.py and replace

if False:
#2D plot example


if True:
#2D plot example

• Plot the 1D radial dependence of density along the equatorial plane

plt.plot(r[:,ny/2,0], rho[:,ny/2,0])
#make the coordinates log-log

Plot 2D (R-z) contours of logarithm of density in color (with color bar):


Here xmax = 100 and ymax = 50 set the plotting range. You can omit these arguments, and a default range will be used.

Compute magnetic flux function:

psi = psicalc()

Overplot magnetic field lines in black:


Compute Lorentz factor:

alpha = (-guu[0,0])**(-0.5) #lapse
gamma = alpha*uu[0]

Compute energy output of a black hole:

aux() #compute OmegaF and T^\mu_\nu
dEr = -gdet*Tud[1,0]*_dx2*_dx3
Er = dEr.sum(2).sum(1) #sum in phi (axis=2) and theta (axis=1)

Plot total radial energy flux vs. radius:

x=r[:,0,0] #x is now the 1D radial coordinate
plt.plot(x,Er) #make a 1D plot of energy flux vs. x
plt.xscale("log") #logarithmic x-axis
plt.xlim(rhor,20) #limit x-range to 20r_g

• Plot the angular frequency of field line rotation vs. theta on the black hole horizon:

OmegaF = omegaf2 #the value of omegaf2 is computed inside aux()
OmegaH = a / (2*rhor)
rhor = 1+(1-a*a)**0.5
ihor = ti[r>=rhor][0] #compute the index of the cell at the location of the horizon

Show ergosphere with thick red line:

rergo = 1+np.sqrt(1-(a*np.cos(h))**2)

Show black hole with black circle:

#generate ellipse object
#import the Ellipse function
from matplotlib.patches import Ellipse
#create Ellipse "actor"
el = Ellipse(
(0,0), #origin of ellipse
2*rhor, 2*rhor, #major axes of ellipse
facecolor='k', #color of ellipse
alpha=1) #opacity of ellipse
ax=plt.gca() #get current axes
art=ax.add_artist(el) #show object in axes
art.set_zorder(20) #bring ellipse to front
plt.draw() #make sure it's drawn

Save figure file


Clear the current window


Open a new window