Practical assignments for MHD simulations using Athena and HARM2D numerical schemes

These assignments are to be carried out a your own pace. At the end of the course, I will collect the answers to the items marked in red with asterisks (*). These will be used by the conference organizers to decide on the amount of financial support for the future schools.

We have secured the access to a computational cluster in Pushchino. You can use it in order to run the simulations instead of your laptop. You can run the problems on your laptop and/or on the cluster. The cluster has the advantage that running there does not drain your battery and the laptop does not get too hot.

Note that the cluster is not intended for carrying out visualization. For that, you will need to copy the simulation results from the cluster to your laptop and use python/VisIt/ParaView installed on your laptop to visualize/analyze those results.

How to login to the cluster:


Here, replace ## with the number assigned to you.

How to download the codes:

You can get Athena here or, if on the cluster, via
cp ~philippo/athena4.1.tar .

You can get HARM2D here or, if on the cluster, via
cp ~atchekho/harm2d.tar .

On the cluster, to run HARM2D you can use the following submit script (copy the following into a file named

#PBS -N test
#PBS -l nodes=1
#PBS -l walltime=01:00:00
#PBS -o test.out
#PBS -e test.err
cd /globaltmp/test##/run #put your home directory here
mpiexec -np 1 ./harm

To run the job, place the executable into the /globaltmp/test##/run directory and type:

qsub #to submit the job

For Athena, replace ./harm with ./athena

In order to check the status of the jobs, do either of the two:
showq #to show all running jobsqstat -u test### #a way to look at your jobs
To kill your job, do:

qdel <jobid>

where <jobid> is the job identifier given to you by the call to showq or qstat.

Please do not hesitate to ask me questions if you have any difficulties running on the cluster! My email is .

1D problems

Athena: hydro problems

1) Familiarize yourself with the Athena MHD code. See the tutorial at
2) If you have not already, download the latest Athena code from and compile the code using the above tutorial.
3) Following the tutorial, run Athena's 1D Sod shock tube test:

- make a plot of the initial and final states of the problem in density, pressure, and velocity
- run another 1D hydro problem: "noh" problem. In this problem, a 1D flow slams against a reflecting wall, shocks, and the shock propagates back into the oncoming flow. Run this test problem by appropriately modifying Athena's input file. Explain what physically is going on in the problem. That is, answer the following questions:

* How large is the density jump in this problem? (*)

* Show that in the limit of strong shock the post/preshock density ratio is rho_2/rho_1 = (gamma+1)/(gamma-1). (*)

* Compute analytically the velocity of the shock in this problem assuming the shock is strong (i.e, neglecting gas pressure in the pre-shock region). (*)

* Make a movie of the time-evolution of density, velocity, or thermal pressure. (*)

Athena: magnetized problems:

4) Athena: Run a 1D MHD problem, e.g., "linear_wave1d" or "cshock1d". Make plots/movies as in item 1.

HARM2D: hydro problems

To run problems with HARM2D and to analyze the results, please follow this tutorial.

5) HARM2D: Run a 1D GR hydro problem: Bondi accretion, i.e., set WHICHPROBLEM to BONDI_PROBLEM_1D in decs.h.Plot the profiles of density at a series of a few times in a simulation. Determine the position of the sonic surface. (*)

HARM: magnetized problems

6) HARM2D: Run 1D monopole problem (i.e., set WHICHPROBLEM to MONOPOLE_PROBLEM_1D in decs.h)
  • Measure the ratio OmegaF/OmegaH. How does it compare with the standard value, 0.5?
  • What is the initial magnetization (near the black hole) of the plasma? Hint: look at the quantity \sigma_0 = b^2/4\pi\rho c^2 = bsq/rho . Here b^2/4\pi = bsq is the square of the fluid frame magnetic field and \rho = rho is the fluid frame density.
  • Make a plot of Lorentz factor vs. radius, Gamma(r). Hint: Gamma = alpha * uu[0]. What value does the Lorentz factor saturate at? How does it compare to the near-black hole magnetization, \sigma_0, of the magnetosphere determined above?

2D problems

1) Athena: 2D Rayleigh-Taylor and/or KH-instabilities (*)
2) HARM2D: 2D BZ-Michel monopole problem, i.e., set WHICHPROBLEM to BZ_MONOPOLE_2D in decs.h . (*)
  • Plot the location of the surface in R-z plane at which the radial contravariant component of velocity, u^r, vanishes. Hint: plot the contour of uu[1] == 0.
  • Plot the dependence of the ratio OmegaF/OmegaH at the horizon as a function of the angle, \theta.
  • Compare the prediction of power by the standard Blandford & Znajek (1977) power formula to the simulation results.
3) HARM2D: 2D monopole problem, i.e., set WHICHPROBLEM to MONOPOLE_PROBLEM_2D in decs.h. Note: runtime is about 16 hours, so run this problem overnight on the cluster.
  • Compare Gamma(r) to that in the above 1D monopole problem. In which case is the Lorentz factor higher? Why?
  • Verify that the value of the Lorentz factor at the fast surface is Gamma = (bsq/rho)^0.5. What is the value of Gamma at the fast surface?

4) HARM2D: Torus problem. Note: runtime is about 16 hours, so run this problem overnight on the cluster.
  • Make a movie of the simulation: logarithm of density shown with color contours overlaid with magnetic field lines. Please come and talk to me for guidance on how to make the movie!

3D problems

1) Athena: choose any 3D problem of your choice, run it (remember to adjust resolution such that you complete it in a reasonable time given the number of processors you are using to run it), and make a 3D volume-rendering of it using VisIt or ParaView. For instance, run 3D magnetic implosion problem. Visualize with VisIt including tracing magnetic field lines. (*)