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 1D shock tests commonly used within the literature, i.e. the Balsara shock tube tests, found in Table 1 here. These tests are performed in flat spacetime.
This notebook is organized as follows
ShockTests_1D
NRPy+ ModuleHere, 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()
Each test defines the GRMHD primitive variables as piece-wise functions, with left and right sides separated by a discontinuity or shock. These primitive variables are the fluid density $\rho$, pressure $p$, velocity $v^i$ and magnetic field $B^i$. For ease in transcriptions we follow Balasra (2001) and define these quantities in Cartesian coordinates. As stated by Balasra (2001), "the initial discontinuity was placed at the center of the computational domain", so we choose to place them at $x=0$.
To aid in speedy calculations, instead of using if
statements to define these piece-wise functions, we use the Min_Max_and_Piecewise_Expressions
NRPy+ module, which replaces if
statements with absolute value calls. In C these are done using the fabs
function.
How do we use Min_Max_and_Piecewise_Expressions
???
import Min_Max_and_Piecewise_Expressions as noif
from sympy import Rational as rl
def balsara1(x, bound = 0.0):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho_left = rl(1.0)
rho_right = rl(1,8)
rho = noif.coord_less_bound(x, bound)*rho_left \
+ noif.coord_geq_bound(x,bound)*rho_right
press_left = rl(1.0)
press_right = rl(1,10)
press = noif.coord_less_bound(x, bound)*press_left \
+ noif.coord_geq_bound(x,bound)*press_right
BU[0] = rl(1,2)
By_left = rl(1.0)
By_right = rl(-1.0)
BU[1] = noif.coord_less_bound(x, bound)*By_left \
+ noif.coord_geq_bound(x,bound)*By_right
return rho, press, vU, BU
def balsara2(x, bound = 0.0):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho = rl(1.0)
press_left = rl(30.0)
press_right = rl(1.0)
press = noif.coord_less_bound(x, bound)*press_left \
+ noif.coord_geq_bound(x,bound)*press_right
BU[0] = rl(1,2)
By_left = rl(6.0)
By_right = rl(7,10)
BU[1] = BU[2] = noif.coord_less_bound(x, bound)*By_left \
+ noif.coord_geq_bound(x,bound)*By_right
return rho, press, vU, BU
def balsara3(x, bound = 0.0):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho = rl(1.0)
press_left = rl(1000.0)
press_right = rl(1,10)
press = noif.coord_less_bound(x, bound)*press_left \
+ noif.coord_geq_bound(x,bound)*press_right
BU[0] = rl(10.)
By_left = rl(7.0)
By_right = rl(7,10)
BU[1] = BU[2] = noif.coord_less_bound(x, bound)*By_left \
+ noif.coord_geq_bound(x,bound)*By_right
return rho, press, vU, BU
def balsara4(x, bound = 0.0):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho = rl(1.0)
press = rl(1,10)
vx_left = rl(999,1000)
vx_right = rl(-999,1000)
vU[0] = noif.coord_less_bound(x, bound)*vx_left \
+ noif.coord_geq_bound(x,bound)*vx_right
BU[0] = rl(10.)
By_left = rl(7.0)
By_right = rl(-7.0)
BU[1] = BU[2] = noif.coord_less_bound(x, bound)*By_left \
+ noif.coord_geq_bound(x,bound)*By_right
return rho, press, vU, BU
def balsara5(x, bound = 0.0):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho_left = rl(108,100)
rho_right = rl(1.0)
rho = noif.coord_less_bound(x, bound)*rho_left \
+ noif.coord_geq_bound(x,bound)*rho_right
press_left = rl(95, 100)
press_right = rl(1.0)
press = noif.coord_less_bound(x, bound)*press_left \
+ noif.coord_geq_bound(x,bound)*press_right
vx_left = rl(4,10)
vx_right = rl(-45,100)
vU[0] = noif.coord_less_bound(x, bound)*vx_left \
+ noif.coord_geq_bound(x,bound)*vx_right
vy_left = rl(3,10)
vy_right = rl(-2,10)
vU[1] = noif.coord_less_bound(x, bound)*vy_left \
+ noif.coord_geq_bound(x,bound)*vy_right
vU[2] = rl(2,10)
BU[0] = rl(2.)
By_left = rl(3,10)
By_right = rl(-7,10)
BU[1] = noif.coord_less_bound(x, bound)*By_left \
+ noif.coord_geq_bound(x,bound)*By_right
Bz_left = rl(3,10)
Bz_right = rl(5,10)
BU[2] = noif.coord_less_bound(x, bound)*Bz_left \
+ noif.coord_geq_bound(x,bound)*Bz_right
return rho, press, vU, BU
Some GRMHD codes elect to evolve the vector potential $A_i$ instead of the magnetic field $B^i$ directly, either to maintain the divergence-free nature of the magnetic fields across grid boundaries or throughout evolution. While we can convert from $A_i \rightarrow B^i$ using
$$ B^i = \epsilon^{ijk} \partial_j A_k, $$where $\epsilon^{ijk}$ is the Levi-Civita tensor, going from $B^i \rightarrow A_i$ is non-trivial for generic magnetic fields. However, given that the fields here are piece-wise constant, for Cartesian coordinates and flat spacetime we can therefore write
$$ A_i = \left( z B^y, x B^z, y B^x \right). $$Here we create a function to deal with implementing the above.
def BtoA_piecewise_constant_Cart_Flat(x, y, z, BU):
AD = ixp.zerorank1()
AD[0] = z*BU[1]
AD[1] = x*BU[2]
AD[2] = y*BU[0]
return AD
Here we define initial data for pure hydrodynamics, i.e. $B^i = 0$. In this section, we construct test problem 1 from Marti and Muller (2003), table 7. With $v^i = 0$ everywhere, we define the density and pressure profiles as
\begin{align} \rho &= \left \{ \begin{array}{lll} 10.0 & \mbox{if} & x < 0.5 \\ 1.0 & \mbox{if} & x \geq 0.5 \end{array} \right. \\ p &= \left \{ \begin{array}{lll} 40.0/3.0 & \mbox{if} & x < 0.5 \\ 0.0 & \mbox{if} & x \geq 0.5 \end{array} \right. \\ \end{align}def hydro1(x, bound = 0.5):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho_left = rl(10.)
rho_right = rl(1.0)
rho = noif.coord_less_bound(x, bound)*rho_left \
+ noif.coord_geq_bound(x,bound)*rho_right
press_left = rl(40, 3)
press_right = rl(0.0)
press = noif.coord_less_bound(x, bound)*press_left \
+ noif.coord_geq_bound(x,bound)*press_right
return rho, press, vU, BU
Similar to the previous section, here we construct test problem 2 from Marti and Muller (2003), table 7. Again with $v^i = 0$ everywhere, we define the density and pressure profiles as
\begin{align} \rho &= 1.0 \\ p &= \left \{ \begin{array}{lll} 1000.0 & \mbox{if} & x < 0.5 \\ 0.01 & \mbox{if} & x \geq 0.5 \end{array} \right. \\ \end{align}def hydro2(x, bound = 0.5):
vU = ixp.zerorank1()
BU = ixp.zerorank1()
rho = rl(1.0)
press_left = rl(1000.0)
press_right = rl(1, 100)
press = noif.coord_less_bound(x, bound)*press_left \
+ noif.coord_geq_bound(x,bound)*press_right
return rho, press, vU, BU
ShockTests_1D
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_1D.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_1D as st_1d
x = rfm.xx_to_Cart[0]
y = rfm.xx_to_Cart[1]
z = rfm.xx_to_Cart[2]
compare_results(balsara1, st_1d.balsara1, x)
compare_results(balsara2, st_1d.balsara2, z)
compare_results(balsara3, st_1d.balsara3, y)
compare_results(balsara4, st_1d.balsara4, z)
compare_results(balsara5, st_1d.balsara5, y)
compare_results(hydro1, st_1d.hydro1, z)
compare_results(hydro2, st_1d.hydro2, x)
_,_,_, BU = balsara5(x, bound=0.23)
AD_A = BtoA_piecewise_constant_Cart_Flat(x,y,z, BU)
AD_B = st_1d.BtoA_piecewise_constant_Cart_Flat(x,y,z, BU)
for i in range(DIM):
assert_zero(AD_A[i] - AD_B[i])
print('All assertions have passed!')
All assertions have passed!
sqrtgammaDET = sp.sympify(1.0)
LeviCivitaTensorUUU = ixp.LeviCivitaTensorUUU_dim3_rank3(sqrtgammaDET)
_, _, _, BU_test = balsara5(y, bound=3.0)
AD_test = BtoA_piecewise_constant_Cart_Flat(rfm.Cartx, rfm.Carty, rfm.Cartz, BU_test)
BU_compare = ixp.zerorank1()
for i in range(DIM):
for j in range(DIM):
for k in range(DIM):
BU_compare[i] += LeviCivitaTensorUUU[i][j][k]*sp.diff(AD_test[k], rfm.Cart[j])
for w in range(DIM):
assert_zero(sp.simplify(BU_test[w] - BU_compare[w]))
print('BtoA passes!')
BtoA passes!
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-1D_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-1D_ShockTests")
Created Tutorial-GRMHD-1D_ShockTests.tex, and compiled LaTeX file to PDF file Tutorial-GRMHD-1D_ShockTests.pdf