IllinoisGRMHD
¶Notebook Status: Self-Validated
Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below. Additional validation tests may have been performed, but are as yet, undocumented. (TODO)
We write the equations of general relativistic hydrodynamics in conservative form as follows (adapted from Eqs. 41-44 of Duez et al):
\begin{eqnarray} \partial_t \tilde{S}_i &+& \partial_j \left(\alpha \sqrt{\gamma} T^j_{{\rm EM}i} \right) = \frac{1}{2} \alpha\sqrt{\gamma} T^{\mu\nu}_{\rm EM} g_{\mu\nu,i}, \end{eqnarray}where we assume $T^{\mu\nu}_{\rm EM}$ is the electromagnetic stress-energy tensor: $$ T^{\mu\nu}_{\rm EM} = b^2 u^{\mu} u^{\nu} + \frac{b^2}{2} g^{\mu\nu} - b^\mu b^\nu, $$ and $$ v^j = \frac{u^j}{u^0} \\ $$
Some quantities can be defined in precisely the same way they are defined in the GRHD equations. *Therefore, we will not define special functions for generating these quantities, and instead refer the user to the appropriate functions in the GRHD module* Namely,
GRHD.compute_S_tildeD(alpha, sqrtgammaDET, T4UD)
GRHD.compute_S_tilde_fluxUD(alpha, sqrtgammaDET, T4UD)
GRHD.compute_S_tilde_source_termD(alpha, sqrtgammaDET,g4DD_zerotimederiv_dD, T4UU)
Also we will write the 4-metric in terms of the ADM 3-metric, lapse, and shift using standard equations.
Thus the full set of input variables include:
Additionally, we need to evolve the vector $A_i$ according to $$ \partial_t A_i = \epsilon_{ijk} v^j B^k - \partial_i (\alpha \Phi - \beta^j A_j) $$ and the scalar potential $\Phi$ according to $$ \partial_t [\sqrt{\gamma} \Phi] + \partial_j (\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi]) = - \xi \alpha [\sqrt{\gamma} \Phi], $$ where $\xi$ determines the strength of the damping term. This is typically set to $0.1$.
As is standard in NRPy+,
For instance, in calculating the first term of $b^2 u^\mu u^\nu$, we use Greek indices:
T4EMUU = ixp.zerorank2(DIM=4)
for mu in range(4):
for nu in range(4):
# Term 1: b^2 u^{\mu} u^{\nu}
T4EMUU[mu][nu] = smallb2*u4U[mu]*u4U[nu]
When we calculate $\beta_i = \gamma_{ij} \beta^j$, we use Latin indices:
betaD = ixp.zerorank1(DIM=3)
for i in range(3):
for j in range(3):
betaD[i] += gammaDD[i][j] * betaU[j]
As a corollary, any expressions involving mixed Greek and Latin indices will need to offset one set of indices by one: A Latin index in a four-vector will be incremented and a Greek index in a three-vector will be decremented (however, the latter case does not occur in this tutorial notebook). This can be seen when we handle $\frac{1}{2} \alpha \sqrt{\gamma} T^{\mu \nu}_{\rm EM} \partial_i g_{\mu \nu}$:
# \alpha \sqrt{\gamma} T^{\mu \nu}_{\rm EM} \partial_i g_{\mu \nu} / 2
for i in range(3):
for mu in range(4):
for nu in range(4):
S_tilde_rhsD[i] += alpsqrtgam * T4EMUU[mu][nu] * g4DD_zerotimederiv_dD[mu][nu][i+1] / 2
Each family of quantities is constructed within a given function (boldfaced below). This notebook is organized as follows
GRFFE.equations
NRPy+ module for GRMHDGRFFE.equations
NRPy+ module for GRFFE# Step 1: Import needed core NRPy+ modules
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
We are given $B^i$, the magnetic field as measured by a normal observer, yet $T^{\mu\nu}_{\rm EM}$ depends on $b^{\mu}$, the magnetic field as measured by an observer comoving with the plasma $B^{\mu}_{\rm (u)}$, divided by $\sqrt{4\pi}$.
In the ideal MHD limit, $B^{\mu}_{\rm (u)}$ is orthogonal to the plasma 4-velocity $u^\mu$, which sets the $\mu=0$ component.
$B^{\mu}_{\rm (u)}$ is related to the magnetic field as measured by a normal observer $B^i$ via a simple projection (Eq 21 in Duez et al (2005)), which results in the expressions (Eqs 23 and 24 in Duez et al (2005)):
\begin{align} \sqrt{4\pi} b^0 = B^0_{\rm (u)} &= \frac{u_j B^j}{\alpha} \\ \sqrt{4\pi} b^i = B^i_{\rm (u)} &= \frac{B^i + (u_j B^j) u^i}{\alpha u^0}\\ \end{align}Further, $B^i$ is related to the actual magnetic field evaluated in IllinoisGRMHD, $\tilde{B}^i$ via
$$B^i = \frac{\tilde{B}^i}{\gamma},$$where $\gamma$ is the determinant of the spatial 3-metric:
# Step 2.a: Define B^i = Btilde^i / sqrt(gamma)
def compute_B_notildeU(sqrtgammaDET, B_tildeU):
global B_notildeU
B_notildeU = ixp.zerorank1(DIM=3)
for i in range(3):
B_notildeU[i] = B_tildeU[i]/sqrtgammaDET
Next we compute Eqs 23 and 24 in Duez et al (2005):
\begin{align} \sqrt{4\pi} b^0 = B^0_{\rm (u)} &= \frac{u_j B^j}{\alpha} \\ \sqrt{4\pi} b^i = B^i_{\rm (u)} &= \frac{B^i + (u_j B^j) u^i}{\alpha u^0}. \end{align}In doing so, we will store the scalar $u_j B^j$ to u4_dot_B_notilde
:
# Step 2.b.i: Define b^mu.
def compute_smallb4U(gammaDD,betaU,alpha, u4U,B_notildeU, sqrt4pi):
global smallb4U
import BSSN.ADMBSSN_tofrom_4metric as AB4m
AB4m.g4DD_ito_BSSN_or_ADM("ADM",gammaDD,betaU,alpha)
u4D = ixp.zerorank1(DIM=4)
for mu in range(4):
for nu in range(4):
u4D[mu] += AB4m.g4DD[mu][nu]*u4U[nu]
smallb4U = ixp.zerorank1(DIM=4)
u4_dot_B_notilde = sp.sympify(0)
for i in range(3):
u4_dot_B_notilde += u4D[i+1]*B_notildeU[i]
# b^0 = (u_j B^j)/[alpha * sqrt(4 pi)]
smallb4U[0] = u4_dot_B_notilde / (alpha*sqrt4pi)
# b^i = [B^i + (u_j B^j) u^i]/[alpha * u^0 * sqrt(4 pi)]
for i in range(3):
smallb4U[i+1] = (B_notildeU[i] + u4_dot_B_notilde*u4U[i+1]) / (alpha*u4U[0]*sqrt4pi)
However, in the case of pure GRFFE (that is, if we assume the pressure and density are zero), we can make an additional simplifying assumption. When deriving the equations of GRFFE, one has a choice as to how to define the velocity. Following the example set by McKinney (2006), section 2.1, and Paschalidis (2013), Appendix A, we choose a referencing frame comoving with the plasma that fulfills the ideal MHD condition (i.e., the electric field is zero). However, in GRFFE, this choice yields a one-parameter family of timelike vectors. By choosing the drift velocity $v^i = u^i/u^0$ that minimizes the Lorentz factor, we find that the four velocity is orthogonal to the magnetic field, or $u_j B^j = 0$. With this assumption,
\begin{align} \sqrt{4\pi} b^0 &= 0 \\ \sqrt{4\pi} b^i &= \frac{B^i}{\alpha u^0}. \end{align}This simplification also gives the inversion from $\tilde{S}_i$ to the drift velocity $v^i$ a closed form, greatly simplifying the conservative-to-primitive solver and removing the need to
# Step 2.b.ii: Define b^mu when u4 and B are orthogonal
def compute_smallb4U_with_driftvU_for_FFE(gammaDD,betaU,alpha, u4U,B_notildeU, sqrt4pi):
global smallb4_with_driftv_for_FFE_U
import BSSN.ADMBSSN_tofrom_4metric as AB4m
AB4m.g4DD_ito_BSSN_or_ADM("ADM",gammaDD,betaU,alpha)
u4D = ixp.zerorank1(DIM=4)
for mu in range(4):
for nu in range(4):
u4D[mu] += AB4m.g4DD[mu][nu]*u4U[nu]
smallb4_with_driftv_for_FFE_U = ixp.zerorank1(DIM=4)
# b^0 = 0
smallb4_with_driftv_for_FFE_U[0] = 0
# b^i = B^i / [alpha * u^0 * sqrt(4 pi)]
for i in range(3):
smallb4_with_driftv_for_FFE_U[i+1] = B_notildeU[i] / (alpha*u4U[0]*sqrt4pi)
Finally we compute smallbsquared
=$b^2 = b_{\mu} b^{\mu} = g_{\mu \nu} b^{\nu}b^{\mu}$:
# Step 2.c: Define b^2.
def compute_smallbsquared(gammaDD,betaU,alpha, smallb4U):
global smallbsquared
import BSSN.ADMBSSN_tofrom_4metric as AB4m
AB4m.g4DD_ito_BSSN_or_ADM("ADM",gammaDD,betaU,alpha)
smallbsquared = sp.sympify(0)
for mu in range(4):
for nu in range(4):
smallbsquared += AB4m.g4DD[mu][nu]*smallb4U[mu]*smallb4U[nu]
Recall from above that
$$ T^{\mu\nu}_{\rm EM} = b^2 u^{\mu} u^{\nu} + \frac{1}{2} b^2 g^{\mu\nu} - b^\mu b^\nu. $$Also
$$ T^\mu_{{\rm EM}\nu} = T^{\mu\delta}_{\rm EM} g_{\delta \nu} $$# Step 3.a: Define T_{EM}^{mu nu} (a 4-dimensional tensor)
def compute_TEM4UU(gammaDD,betaU,alpha, smallb4U, smallbsquared,u4U):
global TEM4UU
# Then define g^{mu nu} in terms of the ADM quantities:
import BSSN.ADMBSSN_tofrom_4metric as AB4m
AB4m.g4UU_ito_BSSN_or_ADM("ADM",gammaDD,betaU,alpha)
# Finally compute T^{mu nu}
TEM4UU = ixp.zerorank2(DIM=4)
for mu in range(4):
for nu in range(4):
TEM4UU[mu][nu] = smallbsquared*u4U[mu]*u4U[nu] \
+ sp.Rational(1,2)*smallbsquared*AB4m.g4UU[mu][nu] \
- smallb4U[mu]*smallb4U[nu]
# Step 3.b: Define T^{mu}_{nu} (a 4-dimensional tensor)
def compute_TEM4UD(gammaDD,betaU,alpha, TEM4UU):
global TEM4UD
# Next compute T^mu_nu = T^{mu delta} g_{delta nu}, needed for S_tilde flux.
# First we'll need g_{alpha nu} in terms of ADM quantities:
import BSSN.ADMBSSN_tofrom_4metric as AB4m
AB4m.g4DD_ito_BSSN_or_ADM("ADM",gammaDD,betaU,alpha)
TEM4UD = ixp.zerorank2(DIM=4)
for mu in range(4):
for nu in range(4):
for delta in range(4):
TEM4UD[mu][nu] += TEM4UU[mu][delta]*AB4m.g4DD[delta][nu]
Now, we will turn our attention to the induction equation, which evolves $A_i$, as well as an additional electromagnetic gauge evolution equation, which evolves $\Phi$. For a cell-centered formulation, they are as follows:
$$
\partial_t A_i = \epsilon_{ijk} v^j B^k - \partial_i (\alpha \Phi - \beta^j A_j),
$$
from Eq. 17 of the original GiRaFFE
paper, and
$$
\partial_t [\sqrt{\gamma} \Phi] + \partial_j (\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi]) = - \xi \alpha [\sqrt{\gamma} \Phi],
$$
from Eq. 19. When it comes to taking the derivatives in these equations, we will write these functions with the intention of computing the operand first, storing it as a gridfunction, and then finite differencing that in a later step.
We'll now take a closer looks at the induction equation, the right-hand side of which has two terms: $$ \partial_t A_i = \underbrace{\epsilon_{ijk} v^j B^k}_{\rm "Flux"\ term} - \partial_i (\underbrace{\alpha \Phi - \beta^j A_j}_{\rm "Source"\ term}), $$ The flux term here is a simple cross-product between the drift velocity and the magnetic field. The source term is the gradient of a gauge-dependent combination of the lapse $\alpha$, the shift $\beta^i$, the vector potential $A_i$, and the scalar potential $\Phi$; we will write a function to compute that operand, and save the finite-difference derivative until later.
def compute_AD_flux_term(sqrtgammaDET,driftvU,BU):
# Levi-Civita tensor for cross products
LeviCivitaDDD = ixp.LeviCivitaTensorDDD_dim3_rank3(sqrtgammaDET)
global A_fluxD
A_fluxD = ixp.zerorank1()
for i in range(3):
for j in range(3):
for k in range(3):
# \epsilon_{ijk} v^j B^k
A_fluxD[i] += LeviCivitaDDD[i][j][k]*driftvU[j]*BU[k]
def compute_AD_source_term_operand_for_FD(sqrtgammaDET,betaU,alpha,psi6Phi,AD):
Phi = psi6Phi/sqrtgammaDET
global AevolParen
# \alpha \Phi - \beta^j A_j
AevolParen = alpha * Phi
for j in range(3):
AevolParen += -betaU[j] * AD[j]
Next, we will build the expressions for the RHS of the evolution equation $$ \partial_t [\sqrt{\gamma} \Phi] = -\partial_j (\underbrace{\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi]}_{\rm Gauge\ evolution\ term}) - \underbrace{\xi \alpha [\sqrt{\gamma} \Phi]}_{\rm Damping\ term}. $$ Once again, we will do this in two parts. First, we will compute the operand of the divergence in the flux term, leaving the finite-difference derivative for later; then, we will compute the damping term.
def compute_psi6Phi_rhs_flux_term_operand(gammaDD,sqrtgammaDET,betaU,alpha,AD,psi6Phi):
gammaUU,_gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
AU = ixp.zerorank1()
# Raise the index on A in the usual way:
for i in range(3):
for j in range(3):
AU[i] += gammaUU[i][j] * AD[j]
global PhievolParenU
PhievolParenU = ixp.zerorank1(DIM=3)
for j in range(3):
# \alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi]
PhievolParenU[j] += alpha*sqrtgammaDET*AU[j] - betaU[j]*psi6Phi
def compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping):
# - \xi \alpha [\sqrt{\gamma} \Phi]
# Combine the divergence and the damping term
global psi6Phi_damping
psi6Phi_damping = - xi_damping * alpha * psi6Phi
The GRFFE equations are given by the induction equation (handled later) and the evolution equation for $\tilde{S}_i$, the Poynting one-form:
\begin{eqnarray} \partial_t \tilde{S}_i &+& \partial_j \left(\alpha \sqrt{\gamma} T^j_{{\rm EM}i} \right) = \frac{1}{2} \alpha\sqrt{\gamma} T^{\mu\nu}_{\rm EM} g_{\mu\nu,i}. \end{eqnarray}Notice that all terms in this equation ($\tilde{S}_i$, $\left(\alpha \sqrt{\gamma} T^j_{{\rm EM}i} \right)$, and the source term) are identical to those in the evolution equation for $\tilde{S}_i$ in general relativistic hydrodynamics (GRHD); one need only replace $T^{\mu\nu}$ of GRHD with the $T^{\mu\nu}_{\rm EM}$ defined above.
Thus we will reuse expressions from the general relativistic hydrodynamics (GRHD) module:
# First define hydrodynamical quantities
u4U = ixp.declarerank1("u4U", DIM=4)
B_tildeU = ixp.declarerank1("B_tildeU", DIM=3)
AD = ixp.declarerank1("D", DIM=3)
B_tildeU = ixp.declarerank1("B_tildeU", DIM=3)
psi6Phi = sp.symbols('psi6Phi', real=True)
# Then ADM quantities
gammaDD = ixp.declarerank2("gammaDD","sym01",DIM=3)
betaU = ixp.declarerank1("betaU", DIM=3)
alpha = sp.symbols('alpha', real=True)
# Then numerical constant
sqrt4pi = sp.symbols('sqrt4pi', real=True)
xi_damping = sp.symbols('xi_damping', real=True)
# First compute stress-energy tensor T4UU and T4UD:
import GRHD.equations as GRHD
GRHD.compute_sqrtgammaDET(gammaDD)
compute_B_notildeU(GRHD.sqrtgammaDET, B_tildeU)
compute_smallb4U(gammaDD,betaU,alpha, u4U, B_notildeU, sqrt4pi)
compute_smallbsquared(gammaDD,betaU,alpha, smallb4U)
compute_TEM4UU(gammaDD,betaU,alpha, smallb4U, smallbsquared,u4U)
compute_TEM4UD(gammaDD,betaU,alpha, TEM4UU)
# Compute conservative variables in terms of primitive variables
GRHD.compute_S_tildeD( alpha, GRHD.sqrtgammaDET, TEM4UD)
S_tildeD = GRHD.S_tildeD
# Next compute fluxes of conservative variables
GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET, TEM4UD)
S_tilde_fluxUD = GRHD.S_tilde_fluxUD
# Then declare derivatives & compute g4DDdD
gammaDD_dD = ixp.declarerank3("gammaDD_dD","sym01",DIM=3)
betaU_dD = ixp.declarerank2("betaU_dD" ,"nosym",DIM=3)
alpha_dD = ixp.declarerank1("alpha_dD" ,DIM=3)
GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)
# Finally compute source terms on tau_tilde and S_tilde equations
GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET, GRHD.g4DD_zerotimederiv_dD, TEM4UU)
S_tilde_source_termD = GRHD.S_tilde_source_termD
# We must also compute the terms for the induction equation and gauge evolution equation
GRHD.compute_vU_from_u4U__no_speed_limit(u4U) # We need the drift velocity
compute_AD_flux_term(GRHD.sqrtgammaDET,GRHD.vU,B_notildeU)
compute_AD_source_term_operand_for_FD(GRHD.sqrtgammaDET,betaU,alpha,psi6Phi,AD)
compute_psi6Phi_rhs_flux_term_operand(gammaDD,GRHD.sqrtgammaDET,betaU,alpha,AD,psi6Phi)
compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping)
GRFFE.equations
NRPy+ module for GRMHD [Back to top]¶As a code validation check, we verify agreement in the SymPy expressions for the GRFFE equations generated in
import GRFFE.equations as GRFFE
# First compute B^i from Btilde^i:
GRFFE.compute_B_notildeU(GRHD.sqrtgammaDET, B_tildeU)
# Then compute b^mu and b^2:
GRFFE.compute_smallb4U(gammaDD, betaU, alpha, u4U, GRFFE.B_notildeU, sqrt4pi)
GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)
# Next construct stress-energy tensor T4UU and T4UD:
GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,u4U)
GRFFE.compute_TEM4UD(gammaDD,betaU,alpha, GRFFE.TEM4UU)
# Compute conservative variables in terms of primitive variables
GRHD.compute_S_tildeD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)
Ge_S_tildeD = GRHD.S_tildeD
# Next compute fluxes of conservative variables
GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)
Ge_S_tilde_fluxUD = GRHD.S_tilde_fluxUD
# Then declare derivatives & compute g4DDdD
# gammaDD_dD = ixp.declarerank3("gammaDD_dD","sym01",DIM=3)
# betaU_dD = ixp.declarerank2("betaU_dD" ,"nosym",DIM=3)
# alpha_dD = ixp.declarerank1("alpha_dD" ,DIM=3)
GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)
# Finally compute source terms on S_tilde equations
GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD,GRFFE.TEM4UU)
Ge_S_tilde_source_termD = GRHD.S_tilde_source_termD
# We must also compute the terms for the induction equation and gauge evolution equation
# GRHD.compute_vU_from_u4U__no_speed_limit(u4U) # We need the drift velocity
GRFFE.compute_AD_flux_term(GRHD.sqrtgammaDET,GRHD.vU,B_notildeU)
GRFFE.compute_AD_source_term_operand_for_FD(GRHD.sqrtgammaDET,betaU,alpha,psi6Phi,AD)
GRFFE.compute_psi6Phi_rhs_flux_term_operand(gammaDD,GRHD.sqrtgammaDET,betaU,alpha,AD,psi6Phi)
GRFFE.compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping)
def comp_func(expr1,expr2,basename,prefixname2="GRFFE."):
if str(expr1-expr2)!="0":
print(basename+" - "+prefixname2+basename+" = "+ str(expr1-expr2))
return 1
return 0
def gfnm(basename,idx1,idx2=None,idx3=None):
if idx2 is None:
return basename+"["+str(idx1)+"]"
if idx3 is None:
return basename+"["+str(idx1)+"]["+str(idx2)+"]"
return basename+"["+str(idx1)+"]["+str(idx2)+"]["+str(idx3)+"]"
expr_list = []
exprcheck_list = []
namecheck_list = []
# PhievolParenU
# psi6Phi_damping
# A_fluxD
# AevolParen
for mu in range(4):
for nu in range(4):
namecheck_list.extend([gfnm("TEM4UU",mu,nu),gfnm("TEM4UD",mu,nu)])
exprcheck_list.extend([GRFFE.TEM4UU[mu][nu],GRFFE.TEM4UD[mu][nu]])
expr_list.extend([TEM4UU[mu][nu],TEM4UD[mu][nu]])
for i in range(3):
namecheck_list.extend([gfnm("S_tildeD",i),gfnm("S_tilde_source_termD",i),gfnm("A_fluxD",i),gfnm("PhievolParenU",i)])
exprcheck_list.extend([Ge_S_tildeD[i],Ge_S_tilde_source_termD[i],GRFFE.A_fluxD[i],GRFFE.PhievolParenU[i]])
expr_list.extend([S_tildeD[i],S_tilde_source_termD[i],A_fluxD[i],PhievolParenU[i]])
for j in range(3):
namecheck_list.extend([gfnm("S_tilde_fluxUD",i,j)])
exprcheck_list.extend([Ge_S_tilde_fluxUD[i][j]])
expr_list.extend([S_tilde_fluxUD[i][j]])
namecheck_list.extend(["AevolParen","psi6Phi_damping"])
exprcheck_list.extend([GRFFE.AevolParen,GRFFE.psi6Phi_damping])
expr_list.extend([AevolParen,psi6Phi_damping])
num_failures = 0
for i in range(len(expr_list)):
num_failures += comp_func(expr_list[i],exprcheck_list[i],namecheck_list[i])
import sys
if num_failures == 0:
print("ALL TESTS PASSED!")
else:
print("ERROR: "+str(num_failures)+" TESTS DID NOT PASS")
sys.exit(1)
ALL TESTS PASSED!
GRFFE.equations
NRPy+ module for GRFFE [Back to top]¶Additionally, we verify agreement in the SymPy expressions for the GRFFE equations generated in
for the case $u_j B^j = 0$.
We will only recompute $b^\mu$ and the expressions that depend on it.
# Generate the expressions within the tutorial, starting with:
compute_smallb4U_with_driftvU_for_FFE(gammaDD,betaU,alpha, u4U, B_notildeU, sqrt4pi)
compute_smallbsquared(gammaDD,betaU,alpha, smallb4_with_driftv_for_FFE_U)
compute_TEM4UU(gammaDD,betaU,alpha, smallb4_with_driftv_for_FFE_U, smallbsquared,u4U)
compute_TEM4UD(gammaDD,betaU,alpha, TEM4UU)
# Compute conservative variables in terms of primitive variables
GRHD.compute_S_tildeD( alpha, GRHD.sqrtgammaDET, TEM4UD)
S_tildeD = GRHD.S_tildeD
# Next compute fluxes of conservative variables
GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET, TEM4UD)
S_tilde_fluxUD = GRHD.S_tilde_fluxUD
# Then declare derivatives & compute g4DDdD
gammaDD_dD = ixp.declarerank3("gammaDD_dD","sym01",DIM=3)
betaU_dD = ixp.declarerank2("betaU_dD" ,"nosym",DIM=3)
alpha_dD = ixp.declarerank1("alpha_dD" ,DIM=3)
GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)
# Finally compute source terms on tau_tilde and S_tilde equations
GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET, GRHD.g4DD_zerotimederiv_dD, TEM4UU)
S_tilde_source_termD = GRHD.S_tilde_source_termD
# Now compute the expressions from the module
# Then compute b^mu and b^2:
GRFFE.compute_smallb4U_with_driftvU_for_FFE(gammaDD, betaU, alpha, u4U, GRFFE.B_notildeU, sqrt4pi)
GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4_with_driftv_for_FFE_U)
# Next construct stress-energy tensor T4UU and T4UD:
GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4_with_driftv_for_FFE_U, GRFFE.smallbsquared,u4U)
GRFFE.compute_TEM4UD(gammaDD,betaU,alpha, GRFFE.TEM4UU)
# Compute conservative variables in terms of primitive variables
GRHD.compute_S_tildeD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)
Ge_S_tildeD = GRHD.S_tildeD
# Next compute fluxes of conservative variables
GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)
Ge_S_tilde_fluxUD = GRHD.S_tilde_fluxUD
# Then declare derivatives & compute g4DDdD
# gammaDD_dD = ixp.declarerank3("gammaDD_dD","sym01",DIM=3)
# betaU_dD = ixp.declarerank2("betaU_dD" ,"nosym",DIM=3)
# alpha_dD = ixp.declarerank1("alpha_dD" ,DIM=3)
GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)
# Finally compute source terms on S_tilde equations
GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD,GRFFE.TEM4UU)
Ge_S_tilde_source_termD = GRHD.S_tilde_source_termD
all_passed=True
expr_list = []
exprcheck_list = []
namecheck_list = []
for mu in range(4):
for nu in range(4):
namecheck_list.extend([gfnm("TEM4UU",mu,nu),gfnm("TEM4UD",mu,nu)])
exprcheck_list.extend([GRFFE.TEM4UU[mu][nu],GRFFE.TEM4UD[mu][nu]])
expr_list.extend([TEM4UU[mu][nu],TEM4UD[mu][nu]])
for i in range(3):
namecheck_list.extend([gfnm("S_tildeD",i),gfnm("S_tilde_source_termD",i)])
exprcheck_list.extend([Ge_S_tildeD[i],Ge_S_tilde_source_termD[i]])
expr_list.extend([S_tildeD[i],S_tilde_source_termD[i]])
for j in range(3):
namecheck_list.extend([gfnm("S_tilde_fluxUD",i,j)])
exprcheck_list.extend([Ge_S_tilde_fluxUD[i][j]])
expr_list.extend([S_tilde_fluxUD[i][j]])
for i in range(len(expr_list)):
comp_func(expr_list[i],exprcheck_list[i],namecheck_list[i])
import sys
if all_passed:
print("ALL TESTS PASSED!")
else:
print("ERROR: AT LEAST ONE TEST DID NOT PASS")
sys.exit(1)
ALL TESTS 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-GRFFE_Equations-Cartesian.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-GRFFE_Equations-Cartesian")
Created Tutorial-GRFFE_Equations-Cartesian.tex, and compiled LaTeX file to PDF file Tutorial-GRFFE_Equations-Cartesian.pdf