GiRaFFEfood Initial Data for GiRaFFE¶Notebook Status: In Progress
Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module through the main initial data modules that depend on it.
We will need to "feed" our giraffe with initial data to evolve. There are several different choices of initial data we can use here; while each represents different physical systems, they all have some steps in common with each other. To avoid code duplication, we will first write several functions that we will use for all of them.
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: Import the NRPy+ core modules and set the reference metric to Cartesian
import NRPy_param_funcs as par
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import reference_metric as rfm
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
# Step 1a: Set commonly used parameters.
thismodule = __name__
First, we will write a function to generate the vector potential from input functions for each component. This function will also apply the correct coordinate staggering if the input is set as such. That is, in the staggered prescription, $A_x$ is sampled at $(i,j+1/2,k+1/2)$, $A_y$ at $(i+1/2,j,k+1/2)$, and $A_z$ at $(i+1/2,j+1/2,k)$.
We will first do this for initial data that are given with Cartesian vector components.
# Generic function for all 1D tests: Compute Ax,Ay,Az
def Axyz_func_Cartesian(Ax_func,Ay_func,Az_func, stagger_enable, **params):
x = rfm.xx_to_Cart[0]
y = rfm.xx_to_Cart[1]
z = rfm.xx_to_Cart[2]
AD = ixp.zerorank1()
# First Ax
if stagger_enable:
y += sp.Rational(1,2)*gri.dxx[1]
z += sp.Rational(1,2)*gri.dxx[2]
AD[0] = Ax_func(x,y,z, **params)
# Then Ay
if stagger_enable:
x += sp.Rational(1,2)*gri.dxx[0]
y -= sp.Rational(1,2)*gri.dxx[1]
z += sp.Rational(1,2)*gri.dxx[2]
AD[1] = Ay_func(x,y,z, **params)
# Finally Az
if stagger_enable:
x += sp.Rational(1,2)*gri.dxx[0]
y += sp.Rational(1,2)*gri.dxx[1]
z -= sp.Rational(1,2)*gri.dxx[2]
AD[2] = Az_func(x,y,z, **params)
return AD
However, some initial data are given with spherical vector components. These will need to be handled slightly differently. We will need to perform a change of basis; fortunately, sympy provides us with powerful tools to do so. We will use these to first write a generic function to transform vectors to a Cartesian basis.
# Generic function to convert contravariant vectors from a spherical to Cartesian basis.
def change_basis_spherical_to_Cartesian_D(somevector_sphD):
# Use the Jacobian matrix to transform the vectors to Cartesian coordinates.
drrefmetric__dx_0UDmatrix = sp.Matrix([[sp.diff(rfm.xxSph[0],rfm.xx[0]), sp.diff(rfm.xxSph[0],rfm.xx[1]), sp.diff(rfm.xxSph[0],rfm.xx[2])],
[sp.diff(rfm.xxSph[1],rfm.xx[0]), sp.diff(rfm.xxSph[1],rfm.xx[1]), sp.diff(rfm.xxSph[1],rfm.xx[2])],
[sp.diff(rfm.xxSph[2],rfm.xx[0]), sp.diff(rfm.xxSph[2],rfm.xx[1]), sp.diff(rfm.xxSph[2],rfm.xx[2])]])
somevectorD = ixp.zerorank1()
for i in range(3):
for j in range(3):
somevectorD[i] += drrefmetric__dx_0UDmatrix[(j,i)]*somevector_sphD[j]
return somevectorD
# Generic function to convert covariant vectors from a spherical to Cartesian basis.
def change_basis_spherical_to_Cartesian_U(somevector_sphU):
# Use the Jacobian matrix to transform the vectors to Cartesian coordinates.
drrefmetric__dx_0UDmatrix = sp.Matrix([[sp.diff(rfm.xxSph[0],rfm.xx[0]), sp.diff(rfm.xxSph[0],rfm.xx[1]), sp.diff(rfm.xxSph[0],rfm.xx[2])],
[sp.diff(rfm.xxSph[1],rfm.xx[0]), sp.diff(rfm.xxSph[1],rfm.xx[1]), sp.diff(rfm.xxSph[1],rfm.xx[2])],
[sp.diff(rfm.xxSph[2],rfm.xx[0]), sp.diff(rfm.xxSph[2],rfm.xx[1]), sp.diff(rfm.xxSph[2],rfm.xx[2])]])
dx__drrefmetric_0UDmatrix = drrefmetric__dx_0UDmatrix.inv()
somevectorU = ixp.zerorank1()
for i in range(3):
for j in range(3):
somevectorU[i] += dx__drrefmetric_0UDmatrix[(j,i)]*somevector_sphU[j]
return somevectorU
# Generic function for all 1D tests: Compute Ax,Ay,Az
def Axyz_func_spherical(Ar_func,At_func,Ap_func, stagger_enable, **params):
if "KerrSchild_radial_shift" in params:
KerrSchild_radial_shift = params["KerrSchild_radial_shift"]
r = rfm.xxSph[0] + KerrSchild_radial_shift # We are setting the data up in Shifted Kerr-Schild coordinates
else:
r = rfm.xxSph[0] # Some other coordinate system
theta = rfm.xxSph[1]
phi = rfm.xxSph[2]
AsphD = ixp.zerorank1()
# First Ax
if stagger_enable:
y += sp.Rational(1,2)*gri.dxx[1]
z += sp.Rational(1,2)*gri.dxx[2]
AsphD[0] = Ar_func(r,theta,phi, **params)
# Then Ay
if stagger_enable:
x += sp.Rational(1,2)*gri.dxx[0]
y -= sp.Rational(1,2)*gri.dxx[1]
z += sp.Rational(1,2)*gri.dxx[2]
AsphD[1] = At_func(r,theta,phi, **params)
# Finally Az
if stagger_enable:
x += sp.Rational(1,2)*gri.dxx[0]
y += sp.Rational(1,2)*gri.dxx[1]
z -= sp.Rational(1,2)*gri.dxx[2]
AsphD[2] = Ap_func(r,theta,phi, **params)
# Use the Jacobian matrix to transform the vectors to Cartesian coordinates.
AD = change_basis_spherical_to_Cartesian(AsphD)
return AD
This function computes the Valenciea 3-velocity from input electric and magnetic fields. It can also take the three-metric $\gamma_{ij}$ as an optional input; if this is not set, the function defaults to flat spacetime.
# Generic function for all 1D tests: Valencia 3-velocity from ED and BU
def compute_ValenciavU_from_ED_and_BU(ED, BU, gammaDD=None):
# Now, we calculate v^i = ([ijk] E_j B_k) / B^2,
# where [ijk] is the Levi-Civita symbol and B^2 = \gamma_{ij} B^i B^j$ is a trivial dot product in flat space.
LeviCivitaSymbolDDD = ixp.LeviCivitaSymbol_dim3_rank3()
B2 = sp.sympify(0)
# In flat spacetime, use the Minkowski metric; otherwise, use the input metric.
if gammaDD is None:
gammaDD = ixp.zerorank2()
for i in range(3):
gammaDD[i][i] = sp.sympify(1)
for i in range(3):
for j in range(3):
B2 += gammaDD[i][j] * BU[i] * BU[j]
BD = ixp.zerorank1()
for i in range(3):
for j in range(3):
BD[i] = gammaDD[i][j]*BU[j]
ValenciavU = ixp.zerorank1()
for i in range(3):
for j in range(3):
for k in range(3):
ValenciavU[i] += LeviCivitaSymbolDDD[i][j][k] * ED[j] * BD[k] / B2
return ValenciavU
def GiRaFFEfood_NRPy_generate_initial_data(ID_type = "DegenAlfvenWave", stagger_enable = False,**params):
global AD, ValenciavU
if ID_type == "ExactWald":
AD = gfcf.Axyz_func_spherical(gfew.Ar_EW,gfew.Ath_EW,gfew.Aph_EW,stagger_enable,**params)
ValenciavU = gfew.ValenciavU_func_EW(**params)
elif ID_type == "MagnetosphericWald":
AD = gfcf.Axyz_func_spherical(gfmw.Ar_MW,gfmw.Ath_MW,gfmw.Aph_MW,stagger_enable,**params)
ValenciavU = gfmw.ValenciavU_func_MW(**params)
elif ID_type == "SplitMonopole":
AD = gfcf.Axyz_func_spherical(gfsm.Ar_SM,gfsm.Ath_SM,gfsm.Aph_SM,stagger_enable,**params)
ValenciavU = gfsm.ValenciavU_func_SM(**params)
elif ID_type == "AlfvenWave":
AD = gfcf.Axyz_func_Cartesian(gfaw.Ax_AW,gfaw.Ay_AW,gfaw.Az_AW, stagger_enable, **params)
ValenciavU = gfaw.ValenciavU_func_AW(**params)
elif ID_type == "FastWave":
AD = gfcf.Axyz_func_Cartesian(gffw.Ax_FW,gffw.Ay_FW,gffw.Az_FW, stagger_enable, **params)
ValenciavU = gffw.ValenciavU_func_FW(**params)
elif ID_type == "DegenAlfvenWave":
AD = gfcf.Axyz_func_Cartesian(gfdaw.Ax_DAW,gfdaw.Ay_DAW,gfdaw.Az_DAW, stagger_enable, **params)
ValenciavU = gfdaw.ValenciavU_func_DAW(**params)
elif ID_type == "ThreeWaves":
AD = gfcf.Axyz_func_Cartesian(gftw.Ax_TW,gftw.Ay_TW,gftw.Az_TW, stagger_enable, **params)
ValenciavU = gftw.ValenciavU_func_TW(**params)
elif ID_type == "FFE_Breakdown":
AD = gfcf.Axyz_func_Cartesian(gffb.Ax_FB,gffb.Ay_FB,gffb.Az_FB, stagger_enable, **params)
ValenciavU = gffb.ValenciavU_func_FB(**params)
elif ID_type == "AlignedRotator":
AD = gfcf.Axyz_func_spherical(gfar.Ar_AR,gfar.Ath_AR,gfar.Aph_AR, stagger_enable, **params)
ValenciavU = gfar.ValenciavU_func_AR(**params)