Notebook Status: Validated
Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below.
Define shocks.
Shocks can occur commonly in astrophysical contexts. Among other methods, shocks may be modeled using codes that solve the equations of magnetohydrodynamics. Here we will develop 1-3D shock tests for MHD codes. Within NRPy+, we will specifically use these tests to validate the NRPy+ version of IllinoisGRMHD. We will borrow the same tests used by the Spritz code for its validation, described here.
In this notebook, we document and implement 2D shock tests commonly used within the literature. Specifically, we will implement the cylindrical explosion, magnetic rotor, and magnetic loop advection tests.
This notebook is organized as follows
Here, we will import the NRPy+ core modules, set the reference metric to Cartesian, and set commonly used NRPy+ parameters. We will also set up a parameter to determine what initial data is set up, although it won't do much yet.
# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
# Step 0.a: Import the NRPy+ core modules and set the reference metric to Cartesian
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par # NRPy+: Parameter interface
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm # NRPy+: Reference metric support
par.set_parval_from_str("grid::DIM", 3)
DIM = par.parval_from_str("grid::DIM")
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
This test is described in section 4.2.1 of Cipolletta et al. Briefly, the test consists of a dense, over-pressurized central region surrounded by an ambient medium, with some smoothing function connecting the two regions. The entire domain has a uniform magnetic field
Similar to our implementation of the 1D Balsara tests, we will make use of Min_Max_and_Piecewise_Expressions
NRPy+ module, to replace if
statements with fabs
calls. The inner region goes out to radius $r_{in} = 0.8$, and the outer region begins at $r_{out} = 1.0$. For the inner region we have $\rho_{in} = 10^{-2}, \ p_{in} = 1.0$. The outer region has $\rho_{out} = 10^{-4}, \ p_{out} = 3 \times 10^{-5}$.
The density and pressure profiles are then given by
\begin{align} \rho\left( r \right) &= \left \{ \begin{array}{lll} \rho_{in} & \mbox{if} & r \leq r_{in} \\ \exp \left[ \frac{\left( r_{out} - r \right) \ln\rho_{in} + \left( r - r_{in} \right) \ln\rho_{out}}{r_{out} - r_{in}} \right] & \mbox{if} & r_{in} < r < r_{out} \\ \rho_{out} & \mbox{if} & r \geq r_{out} \end{array} \right.\\ p\left( r \right) &= \left \{ \begin{array}{lll} p_{in} & \mbox{if} & r \leq r_{in} \\ \exp \left[ \frac{\left( r_{out} - r \right) \ln p_{in} + \left( r - r_{in} \right) \ln p_{out}}{r_{out} - r_{in}} \right] & \mbox{if} & r_{in} < r < r_{out} \\ p_{out} & \mbox{if} & r \geq r_{out} \end{array} \right. \\ \end{align}We also have $v^i = 0$ and $B^x = 0.1, \ B^y = B_z = 0$.
import Min_Max_and_Piecewise_Expressions as noif
from sympy import Rational as rl
def cylindrical_explosion(r, r_in=0.8, r_out=1.0):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
r_out_minus_r_in = r_out - r_in
r_out_minus_r = r_out - r
r_minus_r_in = r - r_in
rho_in = rl(1,100)
rho_out = rl(1,10000)
rho_mid = sp.exp( (r_out_minus_r*sp.ln(rho_in) + r_minus_r_in*sp.ln(rho_out) ) / (r_out_minus_r_in) )
rho = noif.coord_leq_bound(r, r_in)*rho_in\
+noif.coord_greater_bound(r, r_in)*noif.coord_less_bound(r, r_out)*rho_mid\
+noif.coord_geq_bound(r, r_out)*rho_out
press_in = rl(1.0)
press_out = rl(3,100000)
press_mid = sp.exp( (r_out_minus_r*sp.ln(press_in) + r_minus_r_in*sp.ln(press_out) ) / (r_out_minus_r_in) )
press = noif.coord_leq_bound(r, r_in)*press_in\
+noif.coord_greater_bound(r, r_in)*noif.coord_less_bound(r, r_out)*press_mid\
+noif.coord_geq_bound(r, r_out)*press_out
BU[0] = rl(1,10)
return rho, press, vU, BU
This test, described in section 4.2.2 of Cipolletta et al, consists of a rapidly spinning dense, central region surrounded by a static ambient medium. Pressure and magnetic field are uniform throughout the entire domain. For the fluid density we have $r_{in} = 0.1, \ \rho_{in} = 10.0,\ \rho_{out} = 1.0$, giving
\begin{align} \rho\left( r \right) &= \left \{ \begin{array}{lll} \rho_{in} & \mbox{if} & r \leq r_{in} \\ \rho_{out} & \mbox{if} & r > r_{in} \end{array} \right.\\ \end{align}The inner region is rotating with uniform angular velocity $\Omega = 9.95$, i.e.
\begin{align} v^\phi\left( r \right) &= \left \{ \begin{array}{lll} \Omega & \mbox{if} & r \leq r_{in} \\ \rho_{out} & \mbox{if} & r > r_{in} \end{array} \right.\\ \end{align}Below we basis transform $v^i = \left(0, \Omega, 0 \right)$ from cylindrical coordinates to Cartesian coordinates. The entire domain has $B^i = \left( 1.0, 0, 0 \right), \ p = 1.0$.
def magnetic_rotor(r, r_in=0.1, Omega=9.95, cartx=rfm.Cartx, carty=rfm.Carty):
cart_list = [cartx, carty]
vU_cyl = ixp.zerorank1()
BU = ixp.zerorank1()
rho_in = rl(10.0)
rho_out = rl(1.0)
rho = noif.coord_leq_bound(r, r_in)*rho_in\
+noif.coord_greater_bound(r, r_in)*rho_out
press = rl(1.0)
par.set_parval_from_str("reference_metric::CoordSystem","Cylindrical")
rfm.reference_metric()
vU_Cyl = ixp.zerorank1()
vU_Cyl[1] = noif.coord_leq_bound(r, r_in)*Omega
Jac_dUCart_dDrfmUD, Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()
vU = rfm.basis_transform_vectorU_from_rfmbasis_to_Cartesian(Jac_dUCart_dDrfmUD, vU_Cyl)
for i in range(DIM):
for j in range(DIM):
vU[i] = vU[i].subs(rfm.xx[j], rfm.Cart_to_xx[j])
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
for i in range(DIM):
for j in range(2):
vU[i] = vU[i].subs(rfm.Cart[j], cart_list[j])
BU[0] = rl(1.0)
return rho, press, vU, BU
This test, described in section 4.2.3 of Cipolletta et al, consists of a magnetized circular loop propagating through an ambient, unmagnetized medium at constant velocity.
The magnetic field is given by
\begin{align} B^x,B^y &= \left \{ \begin{array}{lll} A_{loop}y/r, A_{loop}x/r & \mbox{if} & r < R_{loop} \\ 0 & \mbox{if} & r \geq R_{loop} \end{array} \right.\\ B^z &= 0 \end{align}The corresponding vector potential may be written as $A_i = \left(0,0, \max\left[0, A_{loop} \left(R_{loop} - r\right)\right] \right)$. Cipolletta et al chooses $\rho=1.0, p=3.0, A_{loop}=0.001, R_{loop}=0.3$
The velocity field is given by
\begin{align} v^x,v^y &= \left \{ \begin{array}{lll} 1/12, 1/24 & \mbox{if} & r < R_{loop} \\ 0 & \mbox{if} & r \geq R_{loop} \end{array} \right.\\ \end{align}For $v^z$, we have the option to set it as $0$ or $1/24$. Cipolletta et al show results for the latter.
def loop_advection(r, R_loop=rl(3,10), A_loop=rl(1/1000), vz_equals_zero=False):
vU = ixp.zerorank1()
AD = ixp.zerorank1()
rho = rl(1.0)
press = rl(3.0)
vx_in = rl(1,12)
vy_in = rl(1,24)
vU[0] = noif.coord_less_bound(r, R_loop)*vx_in
vU[1] = noif.coord_less_bound(r, R_loop)*vy_in
if not vz_equals_zero:
vU[2] = noif.coord_less_bound(r, R_loop)*vy_in
AD[2] = noif.max_noif(0, A_loop*(R_loop - r))
return rho, press, vU, AD
ShockTests_2D
NRPy+ Module [Back to top]¶Here, as a code validation check, we verify agreement in the SymPy expressions for the 1D shock tests we intend to use between
ShockTests_2D.py
module.# here we modify the check_zero function
from UnitTesting.assert_equal import check_zero
def assert_zero(exp):
assert check_zero(exp)==True
def compare_results(functionA, functionB, x, bound=0.0):
rho_A, press_A, vU_A, BU_A = functionA(x, bound)
rho_B, press_B, vU_B, BU_B = functionB(x, bound)
assert_zero(rho_A - rho_B)
assert_zero(press_A - press_B)
for i in range(DIM):
assert_zero(vU_A[i] - vU_B[i])
assert_zero(BU_A[i] - BU_B[i])
import ShockTests_2D as st_2d
r = rfm.xxSph[0]
compare_results(cylindrical_explosion, st_2d.cylindrical_explosion, r)
compare_results(magnetic_rotor, st_2d.magnetic_rotor, r)
compare_results(loop_advection, st_2d.loop_advection, r)
print('All assertions have passed!')
All assertions have passed!
The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-GRMHD-2D_ShockTests.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GRMHD-2D_ShockTests")
Created Tutorial-GRMHD-2D_ShockTests.tex, and compiled LaTeX file to PDF file Tutorial-GRMHD-2D_ShockTests.pdf