Tests

has been tested extensively in linear, nonlinear, special and general relativistic regimes. The tests below are grouped according to the physical model being solved, with subsections describing individual tests.

Extended MHD

Linear modes

It is important to check that any numerical implementation of the EMHD model can reproduce the corresponding linear theory with an error that falls off at the expected order of spatio-temporal discretization. In order to perform this test, one needs to first derive the linear theory of the EMHD model.

The governing equations of the EMHD model are considerably more complex than those of ideal MHD. In particular, the inclusion of both anisotropic pressure and conduction, which are sourced by spatio-temporal derivatives projected along the magnetic field lines make it challenging to derive the linear theory of this model. This is especially so for arbitrary inclinations of the wave vectors with the magnetic field . Even in special cases where and , the derivation is prone to errors if done manually. To address this issue, we have written a general linear analysis package 1, built on top of the 2 computer algebra system, which takes as input the governing equations of any model, and generates the characteristic matrix of the corresponding linear theory. The eigenvectors of this matrix are then used as initial conditions in , and their numerical evolutions checked against the corresponding analytic solutions.

With , we are able to compute the linear modes of the EMHD model for a relativistically hot configuration where rest mass energy density is comparable to the internal energy density , including both anisotropic pressure and conduction, and for any inclinations of the wave vector with the magnetic field .

To thoroughly test the numerical implementation, we choose a mode which excites the variables , with a wave vector , that is misaligned with the background magnetic field . Evidently, both the wave propagation vector and the background magnetic field are misaligned with the numerical grid.

The initial conditions are set with the following eigenmode having eigenvalue = -0.5533585207638108 - 3.626257128688849.

Variable Background state Perturbed value
1. -0.518522524082246-0.1792647678001878
1. 0.5516170736393813
0. 0.008463122479547856+0.011862022608466367
0. -0.16175466371870734-0.034828080823603294
0. 0.
0.1 -0.05973794979640743-0.03351707506150924
0.3 0.02986897489820372+0.016758537530754618
0. 0.
0. 0.5233486841539429+0.04767672501939605
0. 0.2909106062057659+0.021594520553365606

All the variables are initialized as , where is the variable, is the background state, is the perturbation and is the amplitude.

To run this test
* Open src/CMakeLists.txt
* Set the problem with set(PROBLEM "linear_modes")
* Open src/problem/linear_modes/CMakeLists.txt
* Set this mode using set(PROBLEM "FULL_EMHD_2D")

The test has been run with the following options in src/problem/linear_modes/CMakeLists.txt

# Time stepping options: EXPLICIT, IMEX or IMPLICIT
set(TIME_STEPPING "IMEX")
set(DT "0.01")
set(DT_DUMP ".1")
set(START_TIME "0.")
set(FINAL_TIME ".5")
set(START_DUMP_COUNTER "0")
set(COURANT ".9")
set(MAX_DT_INCREMENT "1.3")

# Domain size. If the problem is 1D then N2 is ignored.
set(N1 "512")
set(N2 "512")

# Geometry
set(METRIC "MINKOWSKI")
set(EPS    "1e-5")

# Domain
set(X1_A  "0.")
set(X1_B  "1.")
set(X2_A  "0.")
set(X2_B  "1.")

# Boundary conditions
set(PHYSICAL_BOUNDARY_LEFT_EDGE   "PERIODIC")
set(PHYSICAL_BOUNDARY_RIGHT_EDGE  "PERIODIC")
set(PHYSICAL_BOUNDARY_TOP_EDGE    "PERIODIC")
set(PHYSICAL_BOUNDARY_BOTTOM_EDGE "PERIODIC")

# Reconstrution options
# MP5, MONOTONIZED_CENTRAL or MIN_MOD
set(RECONSTRUCTION "MONOTONIZED_CENTRAL")

# Initial condition parameters 
# Mode options: 
# 1) ENTROPY_WAVE_1D
# 2) HYDRO_SOUND_MODE_1D
# 3) CONDUCTION_STABLE_1D
# 4) CONDUCTION_STABLE_2D
# 5) VISCOSITY_2D
# 6) VISCOSITY_1D
# 7) ALFVEN_2D
# 8) FIREHOSE
# 9) FULL_EMHD_2D
set(AMPLITUDE "1e-8")
set(MODE      "FULL_EMHD_2D")

The plot below shows that all evolved variables in converge to their respective analytic solutions at the expected second order. linear_modes_convergence

EMHD Shock solutions

The presence of sufficient viscosity can smoothen a shock and connect the left and the right sides with a well-defined solution. The hyperbolic nature of the dissipation in the EMHD model leads to new features in the shock structure which have been qualitatively described in Chandra et. al. Here, we solve the shock structure in the EMHD model as a boundary value problem with the left and the right states fixed to their values given by the Rankine-Hugoniot jump conditions. We then use this as a reference solution to check the EMHD shock solutions obtained from , which solves the EMHD equations as an initial value problem.

The boundary value solutions are obtained using a global Newton root finder. We are looking for a steady state nonlinear solution of the EMHD equations and hence set the time derivatives . Since we are interested in the shock sub-structure, we approximate all spatial derivatives by central differences with a truncation error . Thus we have a set of coupled discrete nonlinear equations , where are the primitive variables at and is the chosen spatial resolution of the numerical grid. The system is iterated upon starting from a smooth initial guess using the Newton’s method combined with a numerical Jacobian assembled to machine precision. The iterations are continued till we achieve machine precision error .

We find that the major contribution to the shock structure comes from the pressure anisotropy with the role of the heat conduction in determining the values of the thermodynamic variables inside the shock being marginal. The EMHD equations have higher order corrections that we expect to contribute in strong nonlinear regimes, and indeed we see that the shock structure differs as we turn on and turn off these terms.

There exists an upper limit to the strength of the shock that can be solved for using the EMHD model. Higher mach number shocks require a larger viscosity (or ) to smoothly connect the left and the right states. Beyond a certain critical value of , the theory loses hyperbolicity, and eventually causality and stability. The root of this problem lies in the fact that ultimately, the theory is a second order perturbation , about an equilibrium and as the dissipative effects become stronger, the validity of the expansion breaks down. We obtain the thresholds of this breakdown by expanding about an non-equilibrium using a background heat flux and assembling the characteristic matrix using . As is increased, hyperbolicity is first broken upon the appearance of imaginary values in the eigenvalues. As is further increased, the real values of the eigenvalues exceed , indicating a breakdown of causality. Finally, there is a loss of stability with the linear modes growing exponentially with the fastest growth at .

rho_viscosity rho_viscosity dP_viscosity

Hydrostatic Conducting Atmosphere

Relativisitic conduction in curved space-time contains qualitatively new features when compared to non-relativistic conduction because the heat flux is driven by red-shifted temperature gradients where is the four-acceleration. We test this effect with a hydrostatic fluid configuration in a Schwarzschild metric.

atmosphere_heat_flux convergence_atmosphere

Bondi Inflow without backreaction

bondi_heat_flux bondi_pressure_anisotropy

Anisotropic Conduction Snake Test

Magneto Thermal Instability

radial_mag_field_growth

Heat Flux Driven buoyancy Instability

theta_mag_field_growth

Firehose Instability

firehose_growth

Ideal MHD tests

Sound modes

Komissarov shock tests

Relativistic MHD shock test suite from Komissarov (1999)3 and reproduced with corrections by Gammie, McKinney and Toth4. The tests involve starting with left and right states for all thermodynamic variables and letting the system evolve to a desired time. The solution is then compared to the known solutions in the literature. The initial states for each test are tabulated below. To run this test suite

The test have all been run with the following options in src/problem/shock_tests/CMakeLists.txt

# Time stepping options: EXPLICIT, IMEX or IMPLICIT
set(TIME_STEPPING "IMEX")
set(DT "1e-5")
set(START_TIME "0.")
set(FINAL_TIME "100.")
set(COURANT "0.5")
set(MAX_DT_INCREMENT "1.3")

# Domain size. If the problem is 1D then N2 is ignored.
set(COMPUTE_DIM "1")
set(N1 "512")

# Physics variables
set(ADIABATIC_INDEX "4./3")
set(CONDUCTION "OFF")
set(VISCOSITY  "ON")

# Geometry
set(METRIC "MINKOWSKI")
set(EPS    "1e-5")

# Domain
set(X1_A  "-4.")
set(X1_B  "4.")

# Boundary conditions
set(PHYSICAL_BOUNDARY_LEFT_EDGE   "OUTFLOW")
set(PHYSICAL_BOUNDARY_RIGHT_EDGE  "OUTFLOW")

# Reconstrution options
# MP5, MONOTONIZED_CENTRAL or MIN_MOD
set(RECONSTRUCTION "MONOTONIZED_CENTRAL")
Fast shock
Variable Left state Right state
1. 25.48
1. 367.5
25. 1.091
0. 0.3923
0. 0.
20. 20.
25.02 49.
0. 0.

fast_shock_rho fast_shock_ux

Slow shock
Variable Left state Right state
1. 3.323
10. 55.36
1.53 0.9571
0. -0.6822
0. 0.
10. 10.
18.28 14.49
0. 0.

slow_shock_rho slow_shock_ux

Switch on slow
Variable Left state Right state
1.78e-3 0.01
0.1 1.
-0.765 0.
-1.386 0.
0. 0.
1. 1.
1.022 0.
0. 0.

switch_on_rho switch_on_ux

Switch off fast
Variable Left state Right state
0.1 0.562
1. 10.
-2. -0.212
0. -0.590
0. 0.
2. 2.
0. 4.71
0. 0.

switch_off_rho switch_off_ux

Alfven wave
Variable Left state Right state
1. 1.
1. 1.
0. 3.7
0. 5.76
0. 0.
3. 3.
3. -6.857
0. 0.
Shock tube 1
Variable Left state Right state
1. 0.1
1000. 1.
0. 0.
0. 0.
0. 0.
1. 1.
0. 0.
0. 0.

shock_tube_1_rho shock_tube_1_ux

Shock tube 2
Variable Left state Right state
1. 0.1
30. 1.
0. 0.
0. 0.
0. 0.
0. 0.
20. 0.
0. 0.

shock_tube_2_rho shock_tube_2_ux

Collision
Variable Left state Right state
1. 1.
1. 1.
5. -5.
0. 0.
0. 0.
10. 10.
10. -10.
0. 0.

collision_rho collision_ux

Orzag-Tang

Komissarov Magnetized Cylindrical Explosion

  1. : A framework for automated linear analysis. Hosted on SageMathCloud

  2. Sage Mathematics Software Version 6.7

  3. A Godunov-type scheme for relativistic magnetohydrodynamics, Komissarov (1999)

  4. HARM: A Numerical Scheme for General Relativistic Magnetohydrodynamics, Gammie, McKinney, and Toth (2003)