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, as documented below. The initial data has validated against the original GiRaFFE
, as documented here.
We need to "feed" our giraffe with initial data to evolve. There are several different choices of initial data we can use here; here, we will only be implementing the "Split Monopole" initial data, given by Table 3 in the original paper. This solution is based on the Blandford-Znajek force-free monopole; it is an approximation for the case of small spin with the solution inverted in the lower hemisphere. The vector potential is
\begin{align}
A_r &= -\frac{aC}{8}\cos\theta \left( 1 + \frac{4M}{r} \right), \\
A_\theta &= 0, \\
A_\phi &= M^2 C [1-\cos \theta + a^2 f(r) \cos \theta \sin^2 \theta],
\end{align}
and the electric field is
\begin{align}
E_r &= -\frac{C a^3}{8\alpha M^3} f'(r) \cos \theta \sin^2 \theta \\
E_\theta &= -\frac{Ca}{8\alpha}[\sin \theta + a^2 f(r) \sin \theta (2 \cos^2 \theta-\sin^2 \theta) ] - \beta^r \sqrt{\gamma} \frac{a C}{8 r^2}\left( 1+\frac{4M}{r}\right) \\
E_\phi &= \frac{\beta^r}{\alpha M} Ca^2 f'(r) \cos \theta \sin^2 \theta,
\end{align}
where
\begin{align}
f(r) =& \ \frac{r^2(2r-3M)}{8M^3} L \left(\frac{2M}{r}\right) \\
&+ \frac{M^2+3Mr-6r^2}{12M^2} \ln \frac{r}{2M} \\
&+ \frac{11}{72} + \frac{M}{3r} + \frac{r}{2M} - \frac{r^2}{M^2}, \\
L(x) =& \ {\rm Li}_2(x) + \frac{1}{2} \ln x \ln (1-x).
\end{align}
The function ${\rm Li}_2(x)$ is known as the dilogarithm function, defined as
$$ {\rm Li}_2(x) = -\int_{0}^{1} \frac{\ln(1-tx)}{t} dt = \sum_{k=1}^{\infty} \frac{x^k}{k^2}. $$
Now, to use this initial data scheme, we need to transform the above into the quantities actually tracked by GiRaFFE
and HydroBase: $A_i$, $B^i$, $\tilde{S}_i$, $v^i$, and $\Phi$. Of these quantities, GiRaFFEfood
will only set $A_i$, $v^i$, and $\Phi=0$, then call a separate function to calculate $\tilde{S}_i$; GiRaFFE
itself will call a function to set $B^i$ before the time-evolution begins. This can be done with eqs. 16 and 18, here given in that same order:
\begin{align}
v_{(n)}^i &= \frac{\epsilon^{ijk} E_j B_k}{B^2} \\
B^i &= \frac{[ijk]}{\sqrt{\gamma}} \partial_j A_k \\
\end{align}
In the simulations, $B^i$ will be calculated numerically from $A_i$; however, it will be useful to analytically calculate $B^i$ to use calculating the initial $v^i$.
This notebook is organized as follows
GiRaFFEfood_NRPy.GiRaFFEfood_NRPy
NRPy+ Module# 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
from outputC import nrpyAbs
import NRPy_param_funcs as par # NRPy+: Parameter interface
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_Common_Functions as gfcf # Some useful functions for GiRaFFE initial data.
import reference_metric as rfm # NRPy+: Reference metric support
import Min_Max_and_Piecewise_Expressions as noif
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
# Step 1a: Set commonly used parameters.
thismodule = "GiRaFFEfood_NRPy_Split_Monopole"
# The solution depends on a constant C
C_SM = par.Cparameters("REAL",thismodule,["C_SM"], 1.0)
import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gf # Import this now to set parameters
We begin by coding the function $f(r)$ with the inputs $r$ and mass $M$. We avoid calling the function simply f
out of an abundance of caution, as we do not want to risk overwriting an identically named function elsewhere.
\begin{align}
f(r) =& \ \frac{r^2(2r-3M)}{8M^3} L \left(\frac{2M}{r}\right) \\
&+ \frac{M^2+3Mr-6r^2}{12M^2} \ln \frac{r}{2M} \\
&+ \frac{11}{72} + \frac{M}{3r} + \frac{r}{2M} - \frac{r^2}{2M^2}, \\
\end{align}
where $L(x) = {\rm Li}_2 (x) + \frac{1}{2} \ln x \ln (1-x)$ and ${\rm Li}_2 (x)$ is known as the dilogarithm function. We will use the C library gsl
to handle this special function. In order to do so, we must tell NRPy+ that nrpyDilog
will be our code word that means "use the gsl
dilogarithm function, gsl_sf_dilog
". This is done by simply creating a new sympy function using Function()
and then adding the name-value pair to the dictionary custom_functions_for_SymPy_ccode
.
nrpyDilog = sp.Function('nrpyDilog')
from outputC import custom_functions_for_SymPy_ccode
custom_functions_for_SymPy_ccode["nrpyDilog"] = "gsl_sf_dilog"
def f_of_r(r,M):
if par.parval_from_str("drop_fr"):
return sp.sympify(0)
x = sp.sympify(2)*M/r
L = sp.sympify(0) + \
noif.coord_greater_bound(x,sp.sympify(0))*noif.coord_less_bound(x,sp.sympify(1))*nrpyDilog(x)\
+sp.Rational(1,2)*sp.log(noif.coord_greater_bound(x,sp.sympify(0))*x + noif.coord_leq_bound(x,sp.sympify(1)))\
*sp.log(noif.coord_less_bound(x,sp.sympify(1))*(sp.sympify(1)-x) + noif.coord_geq_bound(x,sp.sympify(1)))
f = r*r*(sp.sympify(2)*r-sp.sympify(3)*M)*sp.Rational(1,8)/(M**3)*L\
+(M*M+sp.sympify(3)*M*r-sp.sympify(6)*r*r)*sp.Rational(1,12)/(M*M)*sp.log(r*sp.Rational(1,2)/M)\
+sp.Rational(11,72) + M*sp.Rational(1,3)/r + r*sp.Rational(1,2)/M - r*r*sp.Rational(1,2)/(M*M)
return f
We will also need the derivative of $f(r)$: \begin{align} f'(r) =& \ \frac{2r(2r-3M)+2r^2}{8M^3} L \left(\frac{2M}{r}\right) + \frac{r^2(2r-3M)}{8M^3} L' \left(\frac{2M}{r}\right) \left( -\frac{2M}{r^2} \right) \\ &+ \frac{3M-12r}{12M^2} \ln \frac{r}{2M} + \frac{M^2+3Mr-6r^2}{12M^2} \frac{2M}{r} \frac{1}{2M} \\ &- \frac{M}{3r^2} + \frac{1}{2M} - \frac{2r}{M^2}. \\ \end{align}
Because $$\frac{\partial {\rm Li}_2 (x)}{\partial x} = \frac{{\rm Li}_1 (x)}{x} = \frac{-\ln (1-x)}{x},$$ we know that \begin{align} L'(x) &= \frac{-\ln (1-x)}{x} + \frac{\ln (1-x)}{2x} - \frac{\ln (x)}{2-2x} \\ &= -\frac{1}{2} \left( \frac{\ln (1-x)}{x} + \frac{\ln(x)}{1-x} \right). \end{align} We simplify this some. \begin{align} f'(r) =& \ \frac{3r^2-3Mr}{4M^3} L \left(\frac{2M}{r}\right) - \frac{2r-3M}{4M^2} L' \left(\frac{2M}{r}\right)\\ &+ \frac{3M-12r}{12M^2} \ln \frac{r}{2M} + \frac{M^2+3Mr-6r^2}{12rM^2} \\ &- \frac{M}{3r^2} + \frac{1}{2M} - \frac{r}{M^2}. \\ \end{align}
def fp_of_r(r,M):
if par.parval_from_str("drop_fr"):
return sp.sympify(0)
x = sp.sympify(2)*M/r
L = sp.sympify(0) + \
noif.coord_greater_bound(x,sp.sympify(0))*noif.coord_less_bound(x,sp.sympify(1))*nrpyDilog(x)\
+sp.Rational(1,2)*sp.log(noif.coord_greater_bound(x,sp.sympify(0))*x + noif.coord_leq_bound(x,sp.sympify(1)))\
*sp.log(noif.coord_less_bound(x,sp.sympify(1))*(sp.sympify(1)-x) + noif.coord_geq_bound(x,sp.sympify(1)))
Lp = sp.sympify(0) + noif.coord_greater_bound(x,sp.sympify(0))*noif.coord_less_bound(x,sp.sympify(1)) * -sp.Rational(1,2) *\
(sp.log(noif.coord_less_bound(x,sp.sympify(1))*(sp.sympify(1)-x) + noif.coord_geq_bound(x,sp.sympify(1)))/(x+sp.sympify(1e-100))\
+sp.log(noif.coord_greater_bound(x,sp.sympify(0))*x + noif.coord_leq_bound(x,sp.sympify(1)))/(sp.sympify(1)-x+sp.sympify(1e-100)))
fp = sp.sympify(3)*r*(r-M)*sp.Rational(1,4)/(M**3)*L + (sp.sympify(2)*r-sp.sympify(3)*M)*sp.Rational(1,4)/(M*M)*Lp\
+(sp.sympify(3)*M-12*r)*sp.Rational(1,12)/(M*M)*sp.log(r*sp.Rational(1,2)/M) + (M*M+sp.sympify(3)*M*r-sp.sympify(6)*r*r)*sp.Rational(1,12)/(r*M*M)\
-M*sp.Rational(1,3)/(r*r) + sp.Rational(1,2)/M - r/(M*M)
return fp
Now, we will code the components of the vector potential $A_i$ in spherical coordinates and a spherical basis. The outputs from these functions can then be easily converted to other coordinate systems by giving the spherical coordinates as inputs in terms of the desired coordinates (e.g., if we want to use Cartesian coordinates, we pass $r = \sqrt{x^2+y^2+z^2}$ and so on). They can also by transformed into any other basis using the appropriate Jacobian matrix.
We will code each component as its own function to more easily apply the appropriate staggering. \begin{align} A_r &= -\frac{aC}{8} \left| \cos \theta \right| \left( 1 + \frac{4M}{r} \right) \sqrt{1 + \frac{2M}{r}}, \\ \end{align}
def Ar_SM(r,theta,phi, **params):
M = params["M"]
a = params["a"]
# A_r = -aC/8 * cos \theta ( 1 + 4M/r ) \sqrt{1 + 2M/r}
return -a*C_SM*sp.Rational(1,8)*nrpyAbs(sp.cos(theta))*(sp.sympify(1)+sp.sympify(4)*M/r)*sp.sqrt(sp.sympify(1)+sp.sympify(2)*M/r)
def Ath_SM(r,theta,phi, **params):
# A_\theta = 0
return sp.sympify(0)
def Aph_SM(r,theta,phi, **params):
M = params["M"]
a = params["a"]
# A_\phi = M^2 C [1-\cos \theta + a^2 f(r) cos \theta sin^2 \theta]
return M*M*C_SM*(sp.sympify(1)-nrpyAbs(sp.cos(theta))+a*a*f_of_r(r,M)*sp.cos(theta)*sp.sin(theta)**2)
Next, we will code up the Valencia velocity, which will require us to first code the electric and magnetic fields. The magnetic field is simply $$B^i = \frac{[ijk]}{\sqrt{\gamma}} \partial_j A_k,$$ which gives \begin{align} B^r &= \frac{C \alpha M^2}{r^2} + \frac{C \alpha a^2 M^2}{2r^4} \left[ -2\cos \theta + \left(\frac{r}{M}\right)^2 (1+3 \cos 2\theta) f(r) \right], \\ B^\theta &= - \frac{C \alpha a^2}{r^2} \sin \theta \cos \theta f'(r), \\ B^\phi &= -\frac{C \alpha a M}{8r^2} \left( 1 + \frac{4M}{r}\right) . \end{align}
The electric field is \begin{align} E_r &= -\frac{C a^3}{8\alpha M^3} f'(r) \cos \theta \sin^2 \theta \\ E_\theta &= -\frac{Ca}{8\alpha}[\sin \theta + a^2 f(r) \sin \theta (2 \cos^2 \theta-\sin^2 \theta) ] - \beta^r \sqrt{\gamma} \frac{a C}{8 r^2}\left( 1+\frac{4M}{r}\right) \\ E_\phi &= \frac{\beta^r}{\alpha M} Ca^2 f'(r) \cos \theta \sin^2 \theta, \end{align}
We can then calculate the the velocity as $$v_{(n)}^i = \frac{\epsilon^{ijk} E_j B_k}{B^2}.$$
def ValenciavU_func_SM(**params):
M = params["M"]
a = params["a"]
alpha = params["alpha"]
betaU = params["betaU"] # Note that this must use a spherical basis!
gammaDD = params["gammaDD"] # Note that this must use a Cartesian basis!
sqrtgammaDET = params["sqrtgammaDET"] # This must be spherical
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
theta = rfm.xxSph[1]
phi = rfm.xxSph[2]
z = rfm.xx_to_Cart[2]
split_C_SM = noif.coord_geq_bound(z,sp.sympify(0))*C_SM - noif.coord_less_bound(z,sp.sympify(0))*C_SM
BsphU = ixp.zerorank1()
BsphU[0] = split_C_SM*alpha*M*M/(r*r) + \
split_C_SM*alpha*a*a*M*M*sp.Rational(1,2)/(r**4)*(-sp.sympify(2)*sp.cos(theta) + (r/M)**2*(sp.sympify(1)+sp.sympify(3)*sp.cos(sp.sympify(2)*theta))*f_of_r(r,M))
BsphU[1] = -split_C_SM*alpha*a*a/(r*r) * sp.sin(theta)*sp.cos(theta)*fp_of_r(r,M)
BsphU[2] = -split_C_SM*alpha*a*M*sp.Rational(1,8)/(r*r)*(sp.sympify(1)+sp.sympify(4)*M/r)
EsphD = ixp.zerorank1()
EsphD[0] = -split_C_SM*a**3/(sp.sympify(8)*alpha*M**3)*fp_of_r(r,M)*sp.cos(theta)*sp.sin(theta)**2
EsphD[1] = -split_C_SM*a*sp.Rational(1,8)/alpha*(sp.sin(theta) + a*a*f_of_r(r,M)*sp.sin(theta)*(sp.sympify(2)*sp.cos(theta)**2-sp.sin(theta)**2)) - \
betaU[0]*sqrtgammaDET*a*split_C_SM*sp.Rational(1,8)/(r*r)*(sp.sympify(1)+sp.sympify(4)*M/r)
EsphD[2] = betaU[0]/(alpha*M)*split_C_SM*a*a*fp_of_r(r,M)*sp.cos(theta)*sp.sin(theta)**2
ED = rfm.basis_transform_vectorD_from_rfmbasis_to_Cartesian(gfcf.Jac_dUrfm_dDCartUD, EsphD)
BU = rfm.basis_transform_vectorU_from_rfmbasis_to_Cartesian(gfcf.Jac_dUCart_dDrfmUD, BsphU)
return gfcf.compute_ValenciavU_from_ED_and_BU(ED, BU, gammaDD)
GiRaFFEfood_NRPy.GiRaFFEfood_NRPy
NRPy+ module [Back to top]¶Here, as a code validation check, we verify agreement in the SymPy expressions for the GiRaFFE
Exact Wald initial data equations we intend to use between
import BSSN.ShiftedKerrSchild as sks
sks.ShiftedKerrSchild(True)
import reference_metric as rfm # NRPy+: Reference metric support
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
gammaSphDD = ixp.zerorank2()
for i in range(3):
for j in range(3):
gammaSphDD[i][j] += sks.gammaSphDD[i][j].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
gammaDD = rfm.basis_transform_tensorDD_from_rfmbasis_to_Cartesian(gfcf.Jac_dUrfm_dDCartUD,gammaSphDD)
unused_gammaUU,gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
sqrtgammaDET = sp.sqrt(gammaDET)
betaU = ixp.zerorank1()
for i in range(3):
betaU[i] += sks.betaSphU[i].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
alpha = sks.alphaSph.subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
A_smD = gfcf.Axyz_func_spherical(Ar_SM,Ath_SM,Aph_SM,stagger_enable = True,M=sks.M,a=sks.a,KerrSchild_radial_shift=sks.r0)
Valenciav_smD = ValenciavU_func_SM(M=sks.M,a=sks.a,KerrSchild_radial_shift=sks.r0,alpha=alpha,betaU=betaU,gammaDD=gammaDD,sqrtgammaDET=sqrtgammaDET)
gf.GiRaFFEfood_NRPy_generate_initial_data(ID_type = "SplitMonopole", stagger_enable = True,M=sks.M,a=sks.a,KerrSchild_radial_shift=sks.r0,alpha=alpha,betaU=betaU,gammaDD=gammaDD,sqrtgammaDET=sqrtgammaDET)
def consistency_check(quantity1,quantity2,string):
if quantity1-quantity2==0:
print(string+" is in agreement!")
else:
print(string+" does not agree!")
sys.exit(1)
for i in range(3):
consistency_check(Valenciav_smD[i],gf.ValenciavU[i],"ValenciavU"+str(i))
consistency_check(A_smD[i],gf.AD[i],"AD"+str(i))
ValenciavU0 is in agreement! AD0 is in agreement! ValenciavU1 is in agreement! AD1 is in agreement! ValenciavU2 is in agreement! AD2 is in agreement!
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-GiRaFFEfood_NRPy-Split_Monopole.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-GiRaFFEfood_NRPy-Split_Monopole")
Created Tutorial-GiRaFFEfood_NRPy-Split_Monopole.tex, and compiled LaTeX file to PDF file Tutorial-GiRaFFEfood_NRPy-Split_Monopole.pdf