This module is currently under development
This module is organized as follows
apply_tau_floor
functionIllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs
functionWe will now use the cmdline_helper.py NRPy+ module to create the source directory within the IllinoisGRMHD
NRPy+ directory if it does not exist yet.
# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: 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 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)
# Step 0c: Create the output file path
outfile_path__tau_and_prim_limits__C = os.path.join(IGM_src_dir_path,"apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C")
In this tutorial notebook, we will discuss how we adjust our conservative variables given that our primitive variables are in the physical range.
For a given set of primitive variables $\left\{\rho_{b},P,v^i,B^i\right\}$ in the physical range (i.e. $\rho_{b}\geq0$, $P\geq0$ and $\epsilon\geq0$), the corresponding conservative variables $\left\{\rho_{\star},\tilde{\tau},\tilde{S}_{i},\tilde{B}^{i}\right\}$ must satisfy certain inequalities (see appendix A of Etienne et al. (2012) for the full discussion). Here we provide a practical recipe to impose these inequalities approximately to reduce inversion failures, which occur mainly in regions with very low density in the artificial “atmosphere” or inside the BH horizon where high accuracy is difficult to maintain but not crucial.
%%writefile $outfile_path__tau_and_prim_limits__C
void eigenvalues_3by3_real_sym_matrix(CCTK_REAL & lam1, CCTK_REAL & lam2, CCTK_REAL & lam3,
CCTK_REAL M11, CCTK_REAL M12, CCTK_REAL M13, CCTK_REAL M22, CCTK_REAL M23, CCTK_REAL M33);
static inline int apply_tau_floor(const int index,const CCTK_REAL tau_atm,const CCTK_REAL rho_b_atm,const CCTK_REAL Psi6threshold,CCTK_REAL *PRIMS,CCTK_REAL *METRIC,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,output_stats &stats,eos_struct &eos, CCTK_REAL *CONSERVS) {
Writing ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
We start by verifying if the metric $\gamma_{ij}$ is positive definite. Notice that although we expect this to always be true, the metric may lose its positive-definiteness due to numerical error during the evolution, especially in the region deep inside the BH, near the “puncture”.
To verify whether or not the metric is positive definite, we analyse its eigenvectors. If the metrix is not positive definite, we reset $\gamma_{ij}\to\psi^{4}\tilde{\gamma}_{ij}$, where $\tilde{\gamma}_{ij}$ corresponds to the 3D flat metric tensor.
%%writefile -a $outfile_path__tau_and_prim_limits__C
//First apply the rho_star floor:
//rho_star = alpha u0 Psi6 rho_b, alpha u0 > 1, so if rho_star < Psi6 rho_b_atm, then we are GUARANTEED that we can reset to atmosphere.
//if(CONSERVS[RHOSTAR] < 1e4*METRIC_LAP_PSI4[PSI6]*rho_b_atm) {
//if(CONSERVS[RHOSTAR] < 2*METRIC_LAP_PSI4[PSI6]*rho_b_atm) {
CCTK_REAL lam1,lam2,lam3;
eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,METRIC[GXX], METRIC[GXY], METRIC[GXZ], METRIC[GYY], METRIC[GYZ], METRIC[GZZ]);
if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0.0) {
// Metric is not positive-defitive, reset the metric to be conformally-flat.
METRIC[GXX] = 1.0;
METRIC[GXY] = 0.0;
METRIC[GXZ] = 0.0;
METRIC[GYY] = 1.0;
METRIC[GYZ] = 0.0;
METRIC[GZZ] = 1.0;
METRIC_PHYS[GUPXX] = METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPXY] = 0.0;
METRIC_PHYS[GUPXZ] = 0.0;
METRIC_PHYS[GUPYY] = METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPYZ] = 0.0;
METRIC_PHYS[GUPZZ] = METRIC_LAP_PSI4[PSIM4];
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
We then set
$$ \boxed{\bar{B}^{i} = \frac{B^{i}}{\sqrt{4\pi}}}\ , $$and
$$ \bar{B}_{i} = \gamma_{ij}\bar{B}^{j} \implies \boxed{ \left\{ \begin{matrix} \bar{B}_{x} = \gamma_{xx}\bar{B}^{x} + \gamma_{xy}\bar{B}^{y} + \gamma_{xz}\bar{B}^{z}\\ \bar{B}_{y} = \gamma_{yx}\bar{B}^{x} + \gamma_{yy}\bar{B}^{y} + \gamma_{yz}\bar{B}^{z}\\ \bar{B}_{z} = \gamma_{zx}\bar{B}^{x} + \gamma_{zy}\bar{B}^{y} + \gamma_{zz}\bar{B}^{z} \end{matrix} \right. }\ , $$then
$$ \bar{B}^{2} \equiv B_{i}B^{i} \implies \boxed{\bar{B}^{2} = B_{x}B^{x} + B_{y}B^{y} + B_{z}B^{z}}\ , $$and finally
$$ \boxed{\bar{B} \equiv \sqrt{\bar{B}^{2}}}\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
//Next, prepare for the tau and stilde fixes:
CCTK_REAL Bxbar = PRIMS[BX_CENTER]*ONE_OVER_SQRT_4PI,Bybar = PRIMS[BY_CENTER]*ONE_OVER_SQRT_4PI,Bzbar = PRIMS[BZ_CENTER]*ONE_OVER_SQRT_4PI;
CCTK_REAL Bbar_x = METRIC_PHYS[GXX]*Bxbar + METRIC_PHYS[GXY]*Bybar + METRIC_PHYS[GXZ]*Bzbar;
CCTK_REAL Bbar_y = METRIC_PHYS[GXY]*Bxbar + METRIC_PHYS[GYY]*Bybar + METRIC_PHYS[GYZ]*Bzbar;
CCTK_REAL Bbar_z = METRIC_PHYS[GXZ]*Bxbar + METRIC_PHYS[GYZ]*Bybar + METRIC_PHYS[GZZ]*Bzbar;
CCTK_REAL Bbar2 = Bxbar*Bbar_x + Bybar*Bbar_y + Bzbar*Bbar_z;
CCTK_REAL Bbar = sqrt(Bbar2);
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
The next part of the code is written to prevent floating-point underflow. We compute $\bar{B}$ in a different way. We start by evaluating
$$ \bar{B}_{\rm check} = \left|\bar{B}^{x}\right| + \left|\bar{B}^{y}\right| + \left|\bar{B}^{z}\right|\ , $$and verifying whether that is a very small, positive number. Then, we determine the largest component of $\bar{B}_{\rm check}$:
$$ \bar{B}_{\max} = \max\left(\left|\bar{B}^{x}\right|,\left|\bar{B}^{y}\right|,\left|\bar{B}^{z}\right|\right)\ . $$Then, we rescale $\bar{B}_{i}$ and $\bar{B}^{i}$ using
$$ \left(\bar{B}^{i}\right)_{\rm tmp} \equiv \frac{\bar{B}^{i}}{\bar{B}_{\max}}\ ,\quad \left(\bar{B}_{i}\right)_{\rm tmp} \equiv \frac{\bar{B}_{i}}{\bar{B}_{\max}}\ , $$and finally recompute $\bar{B}$
$$ \bar{B} = \left[\left(\bar{B}_{i}\right)_{\rm tmp}\left(\bar{B}^{i}\right)_{\rm tmp}\right]\bar{B}_{\max}\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
CCTK_REAL check_B_small = fabs(Bxbar)+fabs(Bybar)+fabs(Bzbar);
if (check_B_small>0 && check_B_small<1.e-150) {
// need to compute Bbar specially to prevent floating-point underflow
CCTK_REAL Bmax = fabs(Bxbar);
if (Bmax < fabs(Bybar)) Bmax=fabs(Bybar);
if (Bmax < fabs(Bzbar)) Bmax=fabs(Bzbar);
CCTK_REAL Bxtmp=Bxbar/Bmax, Bytemp=Bybar/Bmax, Bztemp=Bzbar/Bmax;
CCTK_REAL B_xtemp=Bbar_x/Bmax, B_ytemp=Bbar_y/Bmax, B_ztemp=Bbar_z/Bmax;
Bbar = sqrt(Bxtmp*B_xtemp + Bytemp*B_ytemp + Bztemp*B_ztemp)*Bmax;
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Then we compute
$$ \bar{B} \cdot \tilde{S} = \bar{B}_{i}\tilde{S}^{i} = \bar{B}_{x}\tilde{S}^{x} + \bar{B}_{y}\tilde{S}^{y} + \bar{B}_{z}\tilde{S}^{z}\ . $$and
$$ \hat{\bar{B}}\cdot\tilde{S} \equiv \frac{\bar{B} \cdot \tilde{S}}{\bar{B}}\ . $$However, if $\bar{B} \ll 1$, we set $\hat{\bar{B}}\cdot\tilde{S}=0$.
%%writefile -a $outfile_path__tau_and_prim_limits__C
CCTK_REAL BbardotS = Bxbar*CONSERVS[STILDEX] + Bybar*CONSERVS[STILDEY] + Bzbar*CONSERVS[STILDEZ];
CCTK_REAL hatBbardotS = BbardotS/Bbar;
if (Bbar<1.e-300) hatBbardotS = 0.0;
// Limit hatBbardotS
//CCTK_REAL max_gammav = 100.0;
//CCTK_REAL rhob_max = CONSERVS[RHOSTAR]/METRIC_LAP_PSI4[PSI6];
//CCTK_REAL hmax = 1.0 + 2.0*rhob_max;
//CCTK_REAL abs_hatBbardotS_max = sqrt(SQR(max_gammav)-1.0)*CONSERVS[RHOSTAR]*hmax;
//if (fabs(hatBbardotS) > abs_hatBbardotS_max) {
// CCTK_REAL fac_reduce = abs_hatBbardotS_max/fabs(hatBbardotS);
// CCTK_REAL hatBbardotS_max = hatBbardotS*fac_reduce;
// CCTK_REAL Bbar_inv = 1.0/Bbar;
// CCTK_REAL hat_Bbar_x = Bbar_x*Bbar_inv;
// CCTK_REAL hat_Bbar_y = Bbar_y*Bbar_inv;
// CCTK_REAL hat_Bbar_z = Bbar_z*Bbar_inv;
// CCTK_REAL sub_fact = hatBbardotS_max - hatBbardotS;
// CONSERVS[STILDEX] += sub_fact*hat_Bbar_x;
// CONSERVS[STILDEY] += sub_fact*hat_Bbar_y;
// CONSERVS[STILDEZ] += sub_fact*hat_Bbar_z;
// hatBbardotS = hatBbardotS_max;
// BbardotS *= fac_reduce;
// CONSERVS[STILDEX] = CONSERVS[STILDEX]; CONSERVS[STILDEY] = CONSERVS[STILDEY]; CONSERVS[STILDEZ] = CONSERVS[STILDEZ];
//}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Next we compute
$$ \tilde{S}^{2} \equiv \tilde{S} \cdot \tilde{S} = \gamma^{ij}\tilde{S}_{i}\tilde{S}_{j}\ , $$i.e.
$$ \boxed{ \begin{align} \tilde{S}^{2} &= \gamma^{xx}\left(\tilde{S}_{x}\right)^{2} + \gamma^{yy}\left(\tilde{S}_{y}\right)^{2} + \gamma^{zz}\left(\tilde{S}_{z}\right)^{2}\\ &+2\left( \gamma^{xy}\tilde{S}_{x}\tilde{S}_{y} +\gamma^{xz}\tilde{S}_{x}\tilde{S}_{z} +\gamma^{yz}\tilde{S}_{y}\tilde{S}_{z} \right) \end{align} }\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
CCTK_REAL sdots= METRIC_PHYS[GUPXX]*SQR(CONSERVS[STILDEX])+METRIC_PHYS[GUPYY]*SQR(CONSERVS[STILDEY])+METRIC_PHYS[GUPZZ]*SQR(CONSERVS[STILDEZ])+2.0*
(METRIC_PHYS[GUPXY]*CONSERVS[STILDEX]*CONSERVS[STILDEY]+METRIC_PHYS[GUPXZ]*CONSERVS[STILDEX]*CONSERVS[STILDEZ]+METRIC_PHYS[GUPYZ]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]);
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Then we compute other useful quantities, which are eqs. (A52), (A53), and (A54) in appendix A of Etienne et al. (2012)
$$ \begin{align} W_{m} &= \psi^{-6}\left[\left(\hat{\bar{B}}\cdot\tilde{S}\right)^{2}+\rho_{\star}^{2}\right]^{1/2}\ ,\\ \tilde{S}_{m}^{2} &= \frac{W_{m}^{2}\tilde{S}^{2} + \left(\bar{B}\cdot\tilde{S}\right)^{2}\left(\bar{B}^{2}+2W_{m}\right)}{\left(W_{m}+\bar{B}^{2}\right)^{2}}\ ,\\ W_{\min} &= \psi^{-6}\left(S_{m}^{2}+\rho_{\star}^{2}\right)^{1/2}\ ,\\ \end{align} $$respectively (notice the slightly different notation between the equations above and the one used in the paper).
%%writefile -a $outfile_path__tau_and_prim_limits__C
CCTK_REAL Wm = sqrt(SQR(hatBbardotS)+ SQR(CONSERVS[RHOSTAR]))/METRIC_LAP_PSI4[PSI6];
CCTK_REAL Sm2 = (SQR(Wm)*sdots + SQR(BbardotS)*(Bbar2+2.0*Wm))/SQR(Wm+Bbar2);
CCTK_REAL Wmin = sqrt(Sm2 + SQR(CONSERVS[RHOSTAR]))/METRIC_LAP_PSI4[PSI6];
CCTK_REAL sdots_fluid_max = sdots;
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
%%writefile -a $outfile_path__tau_and_prim_limits__C
//tau fix, applicable when B==0 and B!=0:
if(CONSERVS[TAUENERGY] < 0.5*METRIC_LAP_PSI4[PSI6]*Bbar2) {
CONSERVS[TAUENERGY] = tau_atm+0.5*METRIC_LAP_PSI4[PSI6]*Bbar2;
stats.failure_checker+=1000000;
}
CCTK_REAL tau_fluid_min = CONSERVS[TAUENERGY] - 0.5*METRIC_LAP_PSI4[PSI6]*Bbar2 - (Bbar2*sdots - SQR(BbardotS))*0.5/(METRIC_LAP_PSI4[PSI6]*SQR(Wmin+Bbar2));
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Then we verify if $\tilde{\tau}_{\min} \geq \tilde{\tau}_{\rm atm}$. If $\tilde{\tau}_{\min} < \tilde{\tau}_{\rm atm}$, then reset $\tilde{\tau}$:
$$ \tilde{\tau} = \tilde{\tau}_{\min} + \frac{\psi^{6}}{2}\bar{B}^{2} + \frac{\bar{B}^{2}\tilde{S}^{2} - \left(\bar{B}\cdot\tilde{S}\right)^{2}}{2\psi^{6}\left(W_{\min}+\bar{B}^{2}\right)^{2}}\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
//Apply Stilde fix when B==0.
//if(PRIMS[BX_CENTER]==0 && PRIMS[BY_CENTER]==0 && PRIMS[BZ_CENTER]==0 && (METRIC_LAP_PSI4[PSI6]>30.0 || CONSERVS[RHOSTAR]/METRIC_LAP_PSI4[PSI6]<100*rho_b_atm)) {
//if(check_B_small < 1.e-300) {
/**********************************
* Piecewise Polytropic EOS Patch *
* Computing Patm *
**********************************/
/* This modification of the code trades the variable
* "gamma_equals2" by the already defined function
* pow().
*
* Also, assuming that Patm < rho_ppoly_tab[0], we skip
* the declaration of new variables to store the
* values of K_ppoly_tab[0] and Gamma_ppoly_tab[0]. Thus:
* -----------------------------------------
* | P_{atm} = K_{0} * rho_{atm}^{Gamma_{0}} |
* -----------------------------------------
*/
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_atm);
CCTK_REAL Patm = eos.K_ppoly_tab[polytropic_index]*pow(rho_b_atm,eos.Gamma_ppoly_tab[polytropic_index]);
if(check_B_small*check_B_small < Patm*1e-32) {
CCTK_REAL rhot=CONSERVS[TAUENERGY]*(CONSERVS[TAUENERGY]+2.0*CONSERVS[RHOSTAR]);
CCTK_REAL safetyfactor = 0.999999;
//if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) safetyfactor=0.99;
if(sdots > safetyfactor*rhot) {
CCTK_REAL rfactm1 = sqrt((safetyfactor*rhot)/sdots);
CONSERVS[STILDEX]*=rfactm1;
CONSERVS[STILDEY]*=rfactm1;
CONSERVS[STILDEZ]*=rfactm1;
stats.failure_checker+=10000000;
}
} else if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) {
//Apply new Stilde fix.
if (tau_fluid_min < tau_atm*1.001) {
tau_fluid_min = tau_atm*1.001;
CONSERVS[TAUENERGY] = tau_fluid_min + 0.5*METRIC_LAP_PSI4[PSI6]*Bbar2 + (Bbar2*sdots - SQR(BbardotS))*0.5/(METRIC_LAP_PSI4[PSI6]*SQR(Wmin+Bbar2));
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Then we check if $\tilde{S}^{2} \leq \tilde{\tau}_{\min}\left(\tilde{\tau}_{\min}+2\rho_{\star}\right)$. If not, we reset $\tilde{S}_{i}$
$$ \tilde{S}_{i}\to \tilde{S}_{i}\sqrt{\frac{\tilde{\tau}_{\min}\left(\tilde{\tau}_{\min}+2\rho_{\star}\right)}{\tilde{S}^{2}}} $$%%writefile -a $outfile_path__tau_and_prim_limits__C
CCTK_REAL LHS = tau_fluid_min*(tau_fluid_min+2.0*CONSERVS[RHOSTAR]);
CCTK_REAL RHS = sdots_fluid_max;
CCTK_REAL safetyfactor = 0.999999;
if(safetyfactor*LHS < RHS) {
CCTK_REAL rfactm1 = sqrt((safetyfactor*LHS)/RHS);
CONSERVS[STILDEX]*=rfactm1;
CONSERVS[STILDEY]*=rfactm1;
CONSERVS[STILDEZ]*=rfactm1;
stats.failure_checker+=100000000;
}
}
return 0;
}
/***********************************************************/
/***********************************************************/
/***********************************************************/
/***********************************************************/
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
The IllinoisGRMHD_ enforce_limits_on_primitives_and_recompute_conservs function
[Back to top]¶We start by imposing physical limits on the primitive variables $\left\{\rho_{b},P,v^{i}\right\}$, using:
enforce_pressure_floor_ceiling()
: documented in the inlined_functions.C
tutorial module.impose_speed_limit_output_u0()
: documented in the inlined_functions.C
tutorial module.%%writefile -a $outfile_path__tau_and_prim_limits__C
void IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs(const int already_computed_physical_metric_and_inverse,CCTK_REAL *U,struct output_stats &stats,eos_struct &eos,
CCTK_REAL *METRIC,CCTK_REAL g4dn[4][4],CCTK_REAL g4up[4][4], CCTK_REAL *TUPMUNU,CCTK_REAL *TDNMUNU,CCTK_REAL *CONSERVS) {
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX];
SET_LAPSE_PSI4(METRIC_LAP_PSI4,METRIC);
// Useful debugging tool, can be used to track fixes:
//CCTK_REAL rho_b_orig=U[RHOB],P_orig=U[PRESSURE],vx_orig=U[VX],vy_orig=U[VY],vz_orig=U[VZ];
/***********************************************************/
// Enforce limits on pressure, density, and v^i
/***********************************************************/
// Density floor:
U[RHOB] = MAX(U[RHOB],rho_b_atm);
// Density ceiling:
U[RHOB] = MIN(U[RHOB],rho_b_max);
// Next set h, the enthalpy:
CCTK_REAL h_enthalpy, P_cold,eps_cold,dPcold_drho,eps_th,Gamma_cold; /* <- Note that in setting h, we need to define several
* other variables. Though some will be unused later
* in this function, they may be useful in other
* functions */
compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(U,eos,Gamma_th,P_cold,eps_cold,dPcold_drho,eps_th,h_enthalpy,Gamma_cold);
// Pressure floor & ceiling:
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,U[RHOB]);
enforce_pressure_floor_ceiling(stats,eos.K_ppoly_tab[polytropic_index],P_cold,METRIC_LAP_PSI4[PSI6],Psi6threshold,U[RHOB],rho_b_atm, U[PRESSURE]);
// Possibly adjusted pressure, so recompute eps & h:
CCTK_REAL eps = eps_cold + (U[PRESSURE]-P_cold)/(Gamma_th-1.0)/U[RHOB];
h_enthalpy = 1.0 + eps + U[PRESSURE]/U[RHOB];
CCTK_REAL uUP[4];
impose_speed_limit_output_u0(METRIC,U,METRIC_LAP_PSI4[PSI4],METRIC_LAP_PSI4[LAPSEINV],stats, uUP[0]);
// Compute u^i. We've already set uUP[0] in the lines above.
for(int ii=0;ii<3;ii++) uUP[UX+ii] = uUP[0]*U[VX+ii];
// Useful debugging tool, part 2: can be used to track fixes:
//if(P_orig!=U[PRESSURE] || rho_b_orig!=U[RHOB] || vx_orig!=U[VX] || vy_orig!=U[VY] || vz_orig!=U[VZ]) {
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Next we compute quantities associated with $b^{i}$, using the compute_smallba_b2_and_u_i_over_u0_psi4()
function which is documented in the inlined_functions.C
tutorial module. We also compute $u_{i}$ from the output of the compute_smallba_b2_and_u_i_over_u0_psi4()
function.
%%writefile -a $outfile_path__tau_and_prim_limits__C
/***************************************************************/
// COMPUTE TUPMUNU, TDNMUNU, AND CONSERVATIVES FROM PRIMITIVES
/***************************************************************/
// Compute b^{\mu}, b^2, and u_i/(u^0 Psi4)
CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI = METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;
CCTK_REAL u_x_over_u0_psi4,u_y_over_u0_psi4,u_z_over_u0_psi4;
CCTK_REAL smallb[NUMVARS_SMALLB];
compute_smallba_b2_and_u_i_over_u0_psi4(METRIC,METRIC_LAP_PSI4,U,uUP[0],ONE_OVER_LAPSE_SQRT_4PI,
u_x_over_u0_psi4,u_y_over_u0_psi4,u_z_over_u0_psi4,smallb);
// Compute u_i; we compute u_0 below.
CCTK_REAL uDN[4] = { 1e200, u_x_over_u0_psi4*uUP[0]*METRIC_LAP_PSI4[PSI4],u_y_over_u0_psi4*uUP[0]*METRIC_LAP_PSI4[PSI4],u_z_over_u0_psi4*uUP[0]*METRIC_LAP_PSI4[PSI4] };
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Then we start evaluating the physical ADM inverse 4-metric, $g^{\mu\nu}$ (eq. 2.119 in Numerical Relativity, by Baumgarte & Shapiro)
$$ \boxed{ g^{\mu\nu} = \begin{pmatrix} -\alpha^{-2} & \alpha^{-2}\beta^{i}\\ \alpha^{-2}\beta^{j} & \gamma^{ij} - \alpha^{-2}\beta^{i}\beta^{j} \end{pmatrix} }\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
// Precompute some useful quantities, for later:
CCTK_REAL alpha_sqrt_gamma=METRIC_LAP_PSI4[LAPSE]*METRIC_LAP_PSI4[PSI6];
CCTK_REAL rho0_h_plus_b2 = (U[RHOB]*h_enthalpy + smallb[SMALLB2]);
CCTK_REAL P_plus_half_b2 = (U[PRESSURE]+0.5*smallb[SMALLB2]);
if(!already_computed_physical_metric_and_inverse) {
// If you like, see Eq 2.119 in Numerical Relativity, by Baumgarte & Shapiro:
CCTK_REAL ONE_OVER_LAPSE_SQUARED = SQR(METRIC_LAP_PSI4[LAPSEINV]);
// g^{\mu \nu} = upper four-metric.
//CCTK_REAL g4up[4][4];
g4up[0][0] = -ONE_OVER_LAPSE_SQUARED;
g4up[0][1] = ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX];
g4up[0][2] = ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTY];
g4up[0][3] = ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTZ];
g4up[1][1] = METRIC[GUPXX]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX]*METRIC[SHIFTX];
g4up[1][2] = METRIC[GUPXY]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX]*METRIC[SHIFTY];
g4up[1][3] = METRIC[GUPXZ]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX]*METRIC[SHIFTZ];
g4up[2][2] = METRIC[GUPYY]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTY]*METRIC[SHIFTY];
g4up[2][3] = METRIC[GUPYZ]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTY]*METRIC[SHIFTZ];
g4up[3][3] = METRIC[GUPZZ]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTZ]*METRIC[SHIFTZ];
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Then we compute the stress-energy tensor, $T^{\mu\nu}$, which is given by eq. (8) in Etienne et al. (2015):
$$ \boxed{ T^{\mu\nu} = \left(\rho_{0}h+b^{2}\right)u^{\mu}u^{\nu} + \left(P + \frac{b^{2}}{2}\right)g^{\mu\nu} - b^{\mu}b^{\nu} }\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
int count;
// Next compute T^{\mu \nu}. This is very useful when computing fluxes and source terms in the GRMHD evolution equations.
// (Eq. 33 in http://arxiv.org/pdf/astro-ph/0503420.pdf):
// T^{mn} = (rho_0 h + b^2) u^m u^n + (P + 0.5 b^2) g^{mn} - b^m b^n, where m and n both run from 0 to 3.
count=0; for(int ii=0;ii<4;ii++) for(int jj=ii;jj<4;jj++) { TUPMUNU[count] = rho0_h_plus_b2*uUP[ii]*uUP[jj] + P_plus_half_b2*g4up[ii][jj] - smallb[SMALLBT+ii]*smallb[SMALLBT+jj]; count++; }
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Next we wish to compute the stress-energy tensor, $T_{\mu\nu}$, which is given by
$$ T_{\mu\nu} = \left(\rho_{0}h+b^{2}\right)u_{\mu}u_{\nu} + \left(P + \frac{b^{2}}{2}\right)g_{\mu\nu} - b_{\mu}b_{\nu}\ . $$We start by setting up the physical ADM 4-metric, $g_{\mu\nu}$ (eq. 2.122 in Numerical Relativity, by Baumgarte & Shapiro)
$$ \boxed{g_{\mu\nu} = \begin{pmatrix} -\alpha^{2} + \beta_{\ell}\beta^{\ell} & \beta_{i}\\ \beta_{j} & \gamma_{ij} \end{pmatrix}}\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
// Next compute T_{\mu \nu}
// T_{mn} = (rho_0 h + b^2) u_m u_n + (P + 0.5 b^2) g_{mn} - b_m b_n, where m and n both run from 0 to 3.
if(!already_computed_physical_metric_and_inverse) {
CCTK_REAL LAPSE_SQUARED=SQR(METRIC_LAP_PSI4[LAPSE]);
CCTK_REAL BETADN[4],BETAUP[4] = { 0.0, METRIC[SHIFTX],METRIC[SHIFTY],METRIC[SHIFTZ] };
lower_4vector_output_spatial_part(METRIC_LAP_PSI4[PSI4],METRIC,BETAUP, BETADN);
//CCTK_REAL g4dn[4][4];
// g_{00} = - alpha^2 + gamma_{ij} beta^i beta^j = - alpha^2 beta_i beta^i
g4dn[0][0] = -LAPSE_SQUARED + (BETAUP[1]*BETADN[1] + BETAUP[2]*BETADN[2] + BETAUP[3]*BETADN[3]);
// g_{0i} = gamma_{ij} beta^j = beta_i
g4dn[0][1] = g4dn[1][0] = BETADN[1];
g4dn[0][2] = g4dn[2][0] = BETADN[2];
g4dn[0][3] = g4dn[3][0] = BETADN[3];
// g_{ij} = gamma_{ij} = Psi^4 \tilde{gamma_ij}
g4dn[1][1] = METRIC[GXX]*METRIC_LAP_PSI4[PSI4];
g4dn[1][2] = g4dn[2][1] = METRIC[GXY]*METRIC_LAP_PSI4[PSI4];
g4dn[1][3] = g4dn[3][1] = METRIC[GXZ]*METRIC_LAP_PSI4[PSI4];
g4dn[2][2] = METRIC[GYY]*METRIC_LAP_PSI4[PSI4];
g4dn[2][3] = g4dn[3][2] = METRIC[GYZ]*METRIC_LAP_PSI4[PSI4];
g4dn[3][3] = METRIC[GZZ]*METRIC_LAP_PSI4[PSI4];
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
%%writefile -a $outfile_path__tau_and_prim_limits__C
CCTK_REAL smallb_lower[NUMVARS_SMALLB];
// FIXME: This could be replaced by the function call
// lower_4vector_output_spatial_part(psi4,smallb,smallb_lower);
// b_a = b^c g_{ac}
for(int ii=0;ii<4;ii++) { smallb_lower[SMALLBT+ii]=0; for(int jj=0;jj<4;jj++) smallb_lower[SMALLBT+ii] += smallb[SMALLBT+jj]*g4dn[ii][jj]; }
// Compute u_0, as we've already computed u_i above.
uDN[0]=0.0; for(int jj=0;jj<4;jj++) uDN[0] += uUP[jj]*g4dn[0][jj];
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
Finally, compute
$$ \boxed{ T_{\mu\nu} = \left(\rho_{0}h+b^{2}\right)u_{\mu}u_{\nu} + \left(P + \frac{b^{2}}{2}\right)g_{\mu\nu} - b_{\mu}b_{\nu} }\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
// Compute T_{\mu \nu}
if(update_Tmunu) {
count=0; for(int ii=0;ii<4;ii++) for(int jj=ii;jj<4;jj++) { TDNMUNU[count] = rho0_h_plus_b2*uDN[ii]*uDN[jj] + P_plus_half_b2*g4dn[ii][jj] - smallb_lower[SMALLBT+ii]*smallb_lower[SMALLBT+jj]; count++; }
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
We now have all the ingredients to recompute the conservative variables
$$ \boxed{ \left\{ \begin{align} \rho_{\star} &= \alpha\sqrt{\gamma}\rho_{b}u^{0}\\ \tilde{S}_{i} &= \left(\rho_{\star}h + \alpha\sqrt{\gamma}u^{0}b^{2}\right)u_{i} - \alpha\sqrt{\gamma}b^{0}b_{i}\\ \tilde{\tau} &= \alpha^{2}\sqrt{\gamma}T^{00} - \rho_{\star} \end{align} \right. }\ . $$%%writefile -a $outfile_path__tau_and_prim_limits__C
//CCTK_VInfo(CCTK_THORNSTRING,"YAY %e",TDNMUNU[0]);
// Finally, compute conservatives:
CONSERVS[RHOSTAR] = alpha_sqrt_gamma * U[RHOB] * uUP[0];
CONSERVS[STILDEX] = CONSERVS[RHOSTAR]*h_enthalpy*uDN[1] + alpha_sqrt_gamma*(uUP[0]*smallb[SMALLB2]*uDN[1] - smallb[SMALLBT]*smallb_lower[SMALLBX]);
CONSERVS[STILDEY] = CONSERVS[RHOSTAR]*h_enthalpy*uDN[2] + alpha_sqrt_gamma*(uUP[0]*smallb[SMALLB2]*uDN[2] - smallb[SMALLBT]*smallb_lower[SMALLBY]);
CONSERVS[STILDEZ] = CONSERVS[RHOSTAR]*h_enthalpy*uDN[3] + alpha_sqrt_gamma*(uUP[0]*smallb[SMALLB2]*uDN[3] - smallb[SMALLBT]*smallb_lower[SMALLBZ]);
// tauL = alpha^2 sqrt(gamma) T^{00} - CONSERVS[RHOSTAR]
CONSERVS[TAUENERGY] = METRIC_LAP_PSI4[LAPSE]*alpha_sqrt_gamma*(rho0_h_plus_b2*SQR(uUP[0]) + P_plus_half_b2*(-SQR(METRIC_LAP_PSI4[LAPSEINV])) - SQR(smallb[SMALLBT])) - CONSERVS[RHOSTAR];
// FIXME: We have TUPMUNU already, should replace ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ by TUPMUNU[0].
}
Appending to ../src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code
# First download the original IllinoisGRMHD source code
import urllib
from os import path
original_IGM_file_url = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C"
original_IGM_file_name = "apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs-original.C"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)
# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
try:
original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
# If all else fails, hope wget does the job
!wget -O $original_IGM_file_path $original_IGM_file_url
# Perform validation
Validation__tau_and_prims_limits__C = !diff $original_IGM_file_path $outfile_path__tau_and_prim_limits__C
if Validation__tau_and_prims_limits__C == []:
# If the validation passes, we do not need to store the original IGM source code file
!rm $original_IGM_file_path
print("Validation test for apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C: FAILED!")
# We also print out the difference between the code generated
# in this tutorial module and the original IGM source code
print("Diff:")
for diff_line in Validation__tau_and_prims_limits__C:
print(diff_line)
Validation test for apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C: FAILED! Diff: 5,9c5 < < CCTK_REAL gamma = eos.gamma_tab[0]; < CCTK_REAL kpoly = eos.k_tab[0]; < int gamma_equals2 = 1; < if (fabs(gamma-2.0) > 1.e-10) gamma_equals2 = 0; --- > 34a31 > 43a41 > 53a52 > 78a78 > 81a82 > 86a88 > 94c96,97 < --- > > 98,105c101,119 < CCTK_REAL Patm; < // FIXME: the operation below could be replaced by < // Patm = kpoly*fasterpow_ppm_reconstruct(rho_b_atm,gamma) < if (gamma_equals2==1) { < Patm = kpoly*rho_b_atm*rho_b_atm; < } else { < Patm = kpoly*pow(rho_b_atm,gamma); < } --- > > /********************************** > * Piecewise Polytropic EOS Patch * > * Computing Patm * > **********************************/ > /* This modification of the code trades the variable > * "gamma_equals2" by the already defined function > * pow(). > * > * Also, assuming that Patm < rho_ppoly_tab[0], we skip > * the declaration of new variables to store the > * values of K_ppoly_tab[0] and Gamma_ppoly_tab[0]. Thus: > * ----------------------------------------- > * | P_{atm} = K_{0} * rho_{atm}^{Gamma_{0}} | > * ----------------------------------------- > */ > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_atm); > CCTK_REAL Patm = eos.K_ppoly_tab[polytropic_index]*pow(rho_b_atm,eos.Gamma_ppoly_tab[polytropic_index]); > 115c129 < CONSERVS[STILDEZ]*=rfactm1; --- > CONSERVS[STILDEZ]*=rfactm1; 120c134 < if (tau_fluid_min < tau_atm*1.001) { --- > if (tau_fluid_min < tau_atm*1.001) { 123a138 > 136c151 < --- > 147a163 > 149a166 > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 150a168 > #endif 152c170 < CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX]; --- > CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX]; 154c172 < --- > 167,171c185,189 < CCTK_REAL h_enthalpy=0, P_cold=0,eps_cold=0,dPcold_drho,eps_th,gamma_cold; /* <- Note that in setting h, we need to define several < * other variables. Though some will be unused later < * in this function, they may be useful in other < * functions */ < compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold(U,eos,P_cold,eps_cold,dPcold_drho,eps_th,h_enthalpy,gamma_cold); --- > CCTK_REAL h_enthalpy, P_cold,eps_cold,dPcold_drho,eps_th,Gamma_cold; /* <- Note that in setting h, we need to define several > * other variables. Though some will be unused later > * in this function, they may be useful in other > * functions */ > compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(U,eos,Gamma_th,P_cold,eps_cold,dPcold_drho,eps_th,h_enthalpy,Gamma_cold); 174c192,193 < enforce_pressure_floor_ceiling(stats,eos.K_poly,P_cold,METRIC_LAP_PSI4[PSI6],Psi6threshold,U[RHOB],rho_b_atm, U[PRESSURE]); --- > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,U[RHOB]); > enforce_pressure_floor_ceiling(stats,eos.K_ppoly_tab[polytropic_index],P_cold,METRIC_LAP_PSI4[PSI6],Psi6threshold,U[RHOB],rho_b_atm, U[PRESSURE]); 177c196 < CCTK_REAL eps = eps_cold + (U[PRESSURE]-P_cold)/(gamma_th-1.0)/U[RHOB]; --- > CCTK_REAL eps = eps_cold + (U[PRESSURE]-P_cold)/(Gamma_th-1.0)/U[RHOB]; 187a207 > 199c219,220 < --- > > 201c222 < CCTK_REAL alpha_sqrt_gamma=METRIC_LAP_PSI4[LAPSE]*METRIC_LAP_PSI4[PSI6]; --- > CCTK_REAL alpha_sqrt_gamma=METRIC_LAP_PSI4[LAPSE]*METRIC_LAP_PSI4[PSI6]; 208c229 < --- > 222a244 > 228a251 > 235c258 < --- > 251a275 > 260a285 > 265a291 > 276a303 >
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-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means).
latex_nrpy_style_path = os.path.join(nrpy_dir_path,"latex_nrpy_style.tplx")
#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.tex
!rm -f Tut*.out Tut*.aux Tut*.log