Relativity

Special Relativity (SR)

Special relativistic dynamics can be enabled with the -s configure option. For example, to compile the code for a relativistic MHD shock tube, we can use the following commands:

cd ~/athena
./configure.py -s -b --prob gr_shock_tube --coord cartesian
make clean
make -j

Note that the plain "shock_tube" problem generator does not support relativity; the "gr_shock_tube" does however. We specify the default Cartesian coordinates for clarity. All coordinate systems that work for Newtonian dynamics work for SR as well. In fact, these are the only coordinates that can be used with SR. In particular, "Minkowski" coordinates are reserved for exercising the full machinery of general relativity (below) with a trivial metric. The -s flag changes the equation of state and Riemann solver to be compatible with relativity, but the rest of the code (time integrator, reconstruction algorithm, etc.) remains unchanged.

Now we can run the code. Several SR-specific input files are included in the repository. Here we will run one of the shock tubes from Mignone, Ugliano, & Bodo 2009.

cd ~/work
~/athena/bin/athena -i ~/athena/inputs/mhd_sr/athinput.mub_1

As with the Newtonian shock tube example, this input file will produce .tab files. We can use the Python line-plotting script to read these files and plot the data:

~/athena/vis/python/plot_lines.py mhd_shock_rel_1.block0.out1.00001.tab x1v \
    Etot,dens show -c r,b -l '$E$,$D$' --x_min=-0.5 --x_max=0.5 --x_label '$x$'

For more details, consult the full User Guide entry on Special Relativity.

General Relativity (GR)

General relativity can be enabled with the -g configure flag. This and -s are mutually exclusive. Like with SR, GR requires the code to use particular methods for the equation of state and Riemann problem. In addition, it requires specification of a metric (i.e. a coordinate system). The metric must be stationary. There are several options available, including Minkowski, Schwarzschild, and Kerr-Schild.

For a demonstration, we will set up a hydrodynamic equilibrium torus around a spinning black hole as per Fishbone & Moncrief 1976, seeded with a weak magnetic field. Because general-relativistic problems are very often global problems in spherical-like coordinates, mesh refinement is critical to performance if one wishes to resolve the equatorial plane. Here we will employ static mesh refinement, and therefore HDF5 output. (See also the tutorial on mesh refinement.)

First we configure and compile the code.

cd ~/athena
./configure.py -g -b --prob gr_torus --coord=kerr-schild --flux hlle \
    --nghost 4 -hdf5
make clean
make -j

Note the GR (-g) and MHD (-b) flags. One can also include -t, which performs a frame transformation at each Riemann problem, allowing the use of advanced SR-enabled Riemann solvers like HLLD. The HLLE solver does not need a frame transformation, so we may as well not incur the expense. Here we specify the spherical-like Kerr-Schild coordinates, which describe Kerr spacetime. MPI is omitted for simplicity, but in practice one will need many processes to evolve such a global simulation.

Running the code with the default input file will produce a simple torus setup.

cd ~/work
~/athena/bin/athena -i ~/athena/inputs/mhd_gr/athinput.fm_torus time/nlim=0

We can visualize this with the Python plotting script for spherical data.

~/athena/vis/python/plot_spherical.py fm_torus.prim.00000.athdf rho show \
    --vmin=1.0e-6 --vmax=1.0e0 --logc -s Bcc --stream_density 2.0

This plots a poloidal slice of density, together with streamlines for the magnetic field.

The default configuration is somewhat coarse. It consists of a single MeshBlock of 64³ cells. In practice, simulations in spherical coordinates often make heavy use of static mesh refinement. The input file has a refinement region already specified, covering at least the region 5π/16 < θ < 11π/16. If we divide the theta-direction into 8 MeshBlocks of 8 cells each, this refinement criterion will refine the 4 closest to the midplane, π/4 < θ < 3π/4.

~/athena/bin/athena -i ~/athena/inputs/mhd_gr/athinput.fm_torus time/nlim=0 \
    job/problem_id=fm_torus_refined mesh/refinement=static meshblock/nx2=8

We can do a quick check of the refinement with the same plotting script.

~/athena/vis/python/plot_spherical.py fm_torus_refined.prim.00000.athdf Levels \
    show -c cool

Going back to density and magnetic field, the torus is much smoother.

~/athena/vis/python/plot_spherical.py fm_torus_refined.prim.00000.athdf rho \
    show --vmin=1.0e-6 --vmax=1.0e0 --logc -s Bcc --stream_density 2.0

Note that MeshBlocks that are thin in the theta-direction are often necessary. This is because the root grid should not have too many cells in the theta-direction, or else the timestep will be slow even without refinement, while at the same time there must be sufficiently many MeshBlocks in the root grid to refine some of them without inducing refinement in the ones touching the poles.

Torus initial conditions are ubiquitous workhorses in GRMHD simulations, and there are several features in the associated problem generator and input file that may prove useful:

• Retrograde and tilted disks are fully supported.
• In addition to the density isocontour magnetic field illustrated above, one can use a vector potential that creates multiple small alternating loops of poloidal field. This is scaled with gas pressure, so approximately constant plasma beta can be attained.
• The history file will contain fluxes of mass, energy, angular momentum, and magnetic field integrated on spheres at any desired radii. (Note that many of the other history variables included in all history dumps do not properly account for the metric.)
• Density and pressure floors can be set as power laws in x1-position.
• Sinusoidal perturbations to the initial velocity can be activated.
• Surfaces of constant theta can be nonuniform, concentrating resolution near the midplane. This is controlled by making mesh/x2rat negative and setting coord/h to some value between 0 and 1. The surfaces will be uniform in the modified polar coordinate introduced in the HARM method paper.
• The user-defined radial boundary conditions limit the primitive (normal observer) radial velocity to 0 if material would otherwise flow into the active region. This is much easier to impose than limiting coordinate u¹.
• The boundary conditions do their best to naturally extend the magnetic field in a way that exactly keeps the ghost cells divergence-free.

For more details on the implementation of GR in Athena++, consult the full User Guide entry on General Relativity.