#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # Tutorial-IllinoisGRMHD: apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
# 
# ## Authors: Leo Werneck & Zach Etienne
# 
# <font color='red'>**This module is currently under development**</font>
# 
# ## This tutorial explains two major functions within IllinoisGRMHD that are used to ensure that the results obtained throughout the simulation are physically sound. It covers methods for ensuring the positive-definiteness of the metric, imposing physical limits on the primitive variables, and re-computing conservative variables.
# 
# ### Required and recommended citations:
# 
# * **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).
# * **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).
# * **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618)).

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$
# 
# This module is organized as follows
# 
# 0. [Step 0](#src_dir): **Source directory creation**
# 1. [Step 1](#introduction): **Introduction**
# 1. [Step 2](#apply_tau_floor): **The `apply_tau_floor` function**
#     1. [Step 2.a](#positive_definiteness_of_the_metric): *Positive-definiteness of the metric*
#     1. [Step 2.b](#barbi_barb_i_barb2_and_barb): *Computing $\bar{B}^{i}$, $\bar{B}_{i}$, $\bar{B}^{2}$, and $\bar{B}$*
#     1. [Step 2.c](#barbdots_hatbarbdots_and_sdots): *Computing $\bar{B}\cdot\tilde{S}$, $\hat{\bar{B}}\cdot\tilde{S}$, and $\tilde{S}^{2}$*
#     1. [Step 2.d](#modifying_tau): *Modifying $\tilde{\tau}$*
#         1. [Step 2.d.i](#wm_sm_and_wmin): $W_{m}$, $\tilde{S}_{m}$, and $W_{\min}$
#         1. [Step 2.d.ii](#tau_min): $\tau_{\min}$
#     1. [Step 2.e](#modifying_tilde_s_i): *Modifying $\tilde{S}_{i}$*
# 1. [Step 3](#enforce_limits_on_primitives_and_recompute_conservs): **The `IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs` function**
#     1. [Step 3.a](#enforce_limits_on_primitives): *Enforcing physical limits on the primitive variables*
#     1. [Step 3.b](#bi_b2_and_u_i): *Computing $b^{i}$, $b^{2}$, and $u_{i}$*
#     1. [Step 3.c](#gupmunu): *The physical ADM inverse 4-metric, $g^{\mu\nu}$*
#     1. [Step 3.d](#tupmunu): *The stress-energy tensor, $T^{\mu\nu}$*
#     1. [Step 3.e](#gdnmunu): *The physical ADM 4-metric, $g_{\mu\nu}$*
#     1. [Step 3.f](#tdnmunu): *The stress-energy tensor, $T_{\mu\nu}$*
#     1. [Step 3.g](#recompute_conservs): *Recomputing the conservative variables*
# 1. [Step 4](#code_validation): **Code validation**
# 1. [Step 5](#latex_pdf_output): **Output this notebook to $\LaTeX$-formatted PDF file**

# <a id='src_dir'></a>
# 
# # Step 0: Source directory creation \[Back to [top](#toc)\]
# $$\label{src_dir}$$
# 
# We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory if it does not exist yet.

# In[1]:


# 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")


# <a id='introduction'></a>
# 
# # Step 1: Introduction \[Back to [top](#toc)\]
# $$\label{introduction}$$
# 
# 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)](https://arxiv.org/pdf/1112.0568.pdf) 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.

# <a id='apply_tau_floor'></a>
# 
# # Step 2: The `apply_tau_floor` function \[Back to [top](#toc)\]
# $$\label{apply_tau_floor}$$
# 
# 
# Here we will simply declare the `apply_tau_floor()` and store some of the input in useful variables.

# In[2]:


get_ipython().run_cell_magic('writefile', '$outfile_path__tau_and_prim_limits__C', 'void eigenvalues_3by3_real_sym_matrix(CCTK_REAL & lam1, CCTK_REAL & lam2, CCTK_REAL & lam3,\n                                      CCTK_REAL M11, CCTK_REAL M12, CCTK_REAL M13, CCTK_REAL M22, CCTK_REAL M23, CCTK_REAL M33);\n\nstatic 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) {\n')


# <a id='positive_definiteness_of_the_metric'></a>
# 
# ## Step 2.a: Positive-definiteness of the metric \[Back to [top](#toc)\]
# $$\label{positive_definiteness_of_the_metric}$$
# 
# 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](https://en.wikipedia.org/wiki/Definiteness_of_a_matrix#Eigenvalues). 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.

# In[3]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  //First apply the rho_star floor:\n\n  //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.\n  //if(CONSERVS[RHOSTAR] < 1e4*METRIC_LAP_PSI4[PSI6]*rho_b_atm) {\n  //if(CONSERVS[RHOSTAR] < 2*METRIC_LAP_PSI4[PSI6]*rho_b_atm) {\n\n  CCTK_REAL lam1,lam2,lam3;\n  eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,METRIC[GXX], METRIC[GXY], METRIC[GXZ], METRIC[GYY], METRIC[GYZ], METRIC[GZZ]);\n  if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0.0) {\n    // Metric is not positive-defitive, reset the metric to be conformally-flat.\n    METRIC[GXX] = 1.0;\n    METRIC[GXY] = 0.0;\n    METRIC[GXZ] = 0.0;\n    METRIC[GYY] = 1.0;\n    METRIC[GYZ] = 0.0;\n    METRIC[GZZ] = 1.0;\n    METRIC_PHYS[GUPXX] = METRIC_LAP_PSI4[PSIM4];\n    METRIC_PHYS[GUPXY] = 0.0;\n    METRIC_PHYS[GUPXZ] = 0.0;\n    METRIC_PHYS[GUPYY] = METRIC_LAP_PSI4[PSIM4];\n    METRIC_PHYS[GUPYZ] = 0.0;\n    METRIC_PHYS[GUPZZ] = METRIC_LAP_PSI4[PSIM4];\n  }\n')


# <a id='barbi_barb_i_barb2_and_barb'></a>
# 
# ## Step 2.b: Computing $\bar{B}^{i}$, $\bar{B}_{i}$, $\bar{B}^{2}$, and $\bar{B}$ \[Back to [top](#toc)\]
# $$\label{barbi_barb_i_barb2_and_barb}$$
# 
# 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}}}\ .
# $$

# In[4]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  //Next, prepare for the tau and stilde fixes:\n\n  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;\n\n  CCTK_REAL Bbar_x = METRIC_PHYS[GXX]*Bxbar + METRIC_PHYS[GXY]*Bybar + METRIC_PHYS[GXZ]*Bzbar;\n  CCTK_REAL Bbar_y = METRIC_PHYS[GXY]*Bxbar + METRIC_PHYS[GYY]*Bybar + METRIC_PHYS[GYZ]*Bzbar;\n  CCTK_REAL Bbar_z = METRIC_PHYS[GXZ]*Bxbar + METRIC_PHYS[GYZ]*Bybar + METRIC_PHYS[GZZ]*Bzbar;\n  CCTK_REAL Bbar2 = Bxbar*Bbar_x + Bybar*Bbar_y + Bzbar*Bbar_z;\n  CCTK_REAL Bbar = sqrt(Bbar2);\n')


# The next part of the code is written to prevent [floating-point underflow](https://en.wikipedia.org/wiki/Arithmetic_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}\ .
# $$

# In[5]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n  CCTK_REAL check_B_small = fabs(Bxbar)+fabs(Bybar)+fabs(Bzbar);\n  if (check_B_small>0 && check_B_small<1.e-150) {\n    // need to compute Bbar specially to prevent floating-point underflow\n    CCTK_REAL Bmax = fabs(Bxbar);\n    if (Bmax < fabs(Bybar)) Bmax=fabs(Bybar);\n    if (Bmax < fabs(Bzbar)) Bmax=fabs(Bzbar);\n    CCTK_REAL Bxtmp=Bxbar/Bmax, Bytemp=Bybar/Bmax, Bztemp=Bzbar/Bmax;\n    CCTK_REAL B_xtemp=Bbar_x/Bmax, B_ytemp=Bbar_y/Bmax, B_ztemp=Bbar_z/Bmax;\n    Bbar = sqrt(Bxtmp*B_xtemp + Bytemp*B_ytemp + Bztemp*B_ztemp)*Bmax;\n  }\n')


# <a id='barbdots_hatbarbdots_and_sdots'></a>
# 
# ## Step 2.c: Computing $\bar{B}\cdot\tilde{S}$, $\hat{\bar{B}}\cdot\tilde{S}$, and $\tilde{S}^{2}$ \[Back to [top](#toc)\]
# $$\label{barbdots_hatbarbdots_and_sdots}$$
# 
# 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$.

# In[6]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n  CCTK_REAL BbardotS = Bxbar*CONSERVS[STILDEX] + Bybar*CONSERVS[STILDEY] + Bzbar*CONSERVS[STILDEZ];\n  CCTK_REAL hatBbardotS = BbardotS/Bbar;\n  if (Bbar<1.e-300) hatBbardotS = 0.0;\n\n  // Limit hatBbardotS\n  //CCTK_REAL max_gammav = 100.0;\n  //CCTK_REAL rhob_max = CONSERVS[RHOSTAR]/METRIC_LAP_PSI4[PSI6];\n  //CCTK_REAL hmax = 1.0 + 2.0*rhob_max;\n  //CCTK_REAL abs_hatBbardotS_max = sqrt(SQR(max_gammav)-1.0)*CONSERVS[RHOSTAR]*hmax;\n  //if (fabs(hatBbardotS) > abs_hatBbardotS_max) {\n  //   CCTK_REAL fac_reduce = abs_hatBbardotS_max/fabs(hatBbardotS);\n  //   CCTK_REAL hatBbardotS_max = hatBbardotS*fac_reduce;\n  //   CCTK_REAL Bbar_inv = 1.0/Bbar;\n  //   CCTK_REAL hat_Bbar_x = Bbar_x*Bbar_inv;\n  //   CCTK_REAL hat_Bbar_y = Bbar_y*Bbar_inv;\n  //   CCTK_REAL hat_Bbar_z = Bbar_z*Bbar_inv;\n  //   CCTK_REAL sub_fact = hatBbardotS_max - hatBbardotS;\n  //   CONSERVS[STILDEX] += sub_fact*hat_Bbar_x;\n  //   CONSERVS[STILDEY] += sub_fact*hat_Bbar_y;\n  //   CONSERVS[STILDEZ] += sub_fact*hat_Bbar_z;\n  //   hatBbardotS = hatBbardotS_max;\n  //   BbardotS *= fac_reduce;\n  //   CONSERVS[STILDEX] = CONSERVS[STILDEX]; CONSERVS[STILDEY] = CONSERVS[STILDEY]; CONSERVS[STILDEZ] = CONSERVS[STILDEZ];\n  //}\n')


# 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}
# }\ .
# $$

# In[7]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  CCTK_REAL sdots= METRIC_PHYS[GUPXX]*SQR(CONSERVS[STILDEX])+METRIC_PHYS[GUPYY]*SQR(CONSERVS[STILDEY])+METRIC_PHYS[GUPZZ]*SQR(CONSERVS[STILDEZ])+2.0*\n    (METRIC_PHYS[GUPXY]*CONSERVS[STILDEX]*CONSERVS[STILDEY]+METRIC_PHYS[GUPXZ]*CONSERVS[STILDEX]*CONSERVS[STILDEZ]+METRIC_PHYS[GUPYZ]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]);\n')


# <a id='modifying_tau'></a>
# 
# ## Step 2.d: Modifying $\tilde{\tau}$ \[Back to [top](#toc)\]
# $$\label{modifying_tau}$$
# 
# <a id='wm_sm_and_wmin'></a>
# 
# ### Step 2.d.i: $W_{m}$, $\tilde{S}_{m}$, and $W_{\min}$ \[Back to [top](#toc)\]
# $$\label{wm_sm_and_wmin}$$
# 
# Then we compute other useful quantities, which are eqs. (A52), (A53), and (A54) in appendix A of [Etienne *et al.* (2012)](https://arxiv.org/pdf/1112.0568.pdf)
# 
# $$
# \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).

# In[8]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  CCTK_REAL Wm = sqrt(SQR(hatBbardotS)+ SQR(CONSERVS[RHOSTAR]))/METRIC_LAP_PSI4[PSI6];\n  CCTK_REAL Sm2 = (SQR(Wm)*sdots + SQR(BbardotS)*(Bbar2+2.0*Wm))/SQR(Wm+Bbar2);\n  CCTK_REAL Wmin = sqrt(Sm2 + SQR(CONSERVS[RHOSTAR]))/METRIC_LAP_PSI4[PSI6];\n  CCTK_REAL sdots_fluid_max = sdots;\n')


# <a id='tau_min'></a>
# 
# ### Step 2.d.ii: $\tilde{\tau}_{\min}$ \[Back to [top](#toc)\]
# $$\label{tau_min}$$
# 
# Next we evaluate
# 
# $$
# \tilde{\tau}_{\min} = \tilde{\tau} - \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}}
# $$

# In[9]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  //tau fix, applicable when B==0 and B!=0:\n  if(CONSERVS[TAUENERGY] < 0.5*METRIC_LAP_PSI4[PSI6]*Bbar2) {\n    CONSERVS[TAUENERGY] = tau_atm+0.5*METRIC_LAP_PSI4[PSI6]*Bbar2;\n    stats.failure_checker+=1000000;\n  }\n\n  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));\n')


# 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}}\ .
# $$

# In[10]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  //Apply Stilde fix when B==0.\n  //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)) {\n  //if(check_B_small < 1.e-300) {\n\n  /**********************************\n   * Piecewise Polytropic EOS Patch *\n   *         Computing Patm         *\n   **********************************/\n  /* This modification of the code trades the variable\n   * "gamma_equals2" by the already defined function\n   * pow().\n   *\n   * Also, assuming that Patm < rho_ppoly_tab[0], we skip\n   * the declaration of new variables to store the\n   * values of K_ppoly_tab[0] and Gamma_ppoly_tab[0]. Thus:\n   *  -----------------------------------------\n   * | P_{atm} = K_{0} * rho_{atm}^{Gamma_{0}} |\n   *  -----------------------------------------\n   */\n  int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_atm);\n  CCTK_REAL Patm = eos.K_ppoly_tab[polytropic_index]*pow(rho_b_atm,eos.Gamma_ppoly_tab[polytropic_index]);\n\n  if(check_B_small*check_B_small < Patm*1e-32) {\n    CCTK_REAL rhot=CONSERVS[TAUENERGY]*(CONSERVS[TAUENERGY]+2.0*CONSERVS[RHOSTAR]);\n    CCTK_REAL safetyfactor = 0.999999;\n    //if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) safetyfactor=0.99;\n\n    if(sdots > safetyfactor*rhot) {\n      CCTK_REAL rfactm1 = sqrt((safetyfactor*rhot)/sdots);\n      CONSERVS[STILDEX]*=rfactm1;\n      CONSERVS[STILDEY]*=rfactm1;\n      CONSERVS[STILDEZ]*=rfactm1;\n      stats.failure_checker+=10000000;\n    }\n  } else if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) {\n    //Apply new Stilde fix.\n    if (tau_fluid_min < tau_atm*1.001) {\n      tau_fluid_min = tau_atm*1.001;\n      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));\n    }\n')


# <a id='modifying_tilde_s_i'></a>
# 
# ## Step 2.e: Modifying $\tilde{S}_{i}$ \[Back to [top](#toc)\]
# $$\label{modifying_tilde_s_i}$$
# 
# 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}}}
# $$

# In[11]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n    CCTK_REAL LHS = tau_fluid_min*(tau_fluid_min+2.0*CONSERVS[RHOSTAR]);\n    CCTK_REAL RHS = sdots_fluid_max;\n\n    CCTK_REAL safetyfactor = 0.999999;\n    if(safetyfactor*LHS < RHS) {\n      CCTK_REAL rfactm1 = sqrt((safetyfactor*LHS)/RHS);\n      CONSERVS[STILDEX]*=rfactm1;\n      CONSERVS[STILDEY]*=rfactm1;\n      CONSERVS[STILDEZ]*=rfactm1;\n      stats.failure_checker+=100000000;\n    }\n  }\n\n\n\n  return 0;\n}\n\n/***********************************************************/\n/***********************************************************/\n/***********************************************************/\n/***********************************************************/\n\n')


# <a id='enforce_limits_on_primitives_and_recompute_conservs'></a>
# 
# # Step 3: `The IllinoisGRMHD_ enforce_limits_on_primitives_and_recompute_conservs function` \[Back to [top](#toc)\]
# $$\label{enforce_limits_on_primitives_and_recompute_conservs}$$
# 
# <a id='enforce_limits_on_primitives'></a>
# 
# ## Step 3.a: Enforcing physical limits on the primitive variables \[Back to [top](#toc)\]
# $$\label{enforce_limits_on_primitives}$$
# 
# We start by imposing physical limits on the primitive variables $\left\{\rho_{b},P,v^{i}\right\}$, using:
# 
# 1. $\rho_{b} \to \min\left(\rho_{b},\rho_{b,{\rm atm}}\right)$
# 1. `enforce_pressure_floor_ceiling()`: documented in the [`inlined_functions.C`](Tutorial-IllinoisGRMHD__inlined_functions.ipynb) tutorial module.
# 1. `impose_speed_limit_output_u0()`: documented in the [`inlined_functions.C`](Tutorial-IllinoisGRMHD__inlined_functions.ipynb) tutorial module.

# In[12]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', "\n\nvoid 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,\n                                                                       CCTK_REAL *METRIC,CCTK_REAL g4dn[4][4],CCTK_REAL g4up[4][4], CCTK_REAL *TUPMUNU,CCTK_REAL *TDNMUNU,CCTK_REAL *CONSERVS) {\n#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n  DECLARE_CCTK_PARAMETERS;\n#endif\n\n  CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX];\n  SET_LAPSE_PSI4(METRIC_LAP_PSI4,METRIC);\n\n  // Useful debugging tool, can be used to track fixes:\n  //CCTK_REAL rho_b_orig=U[RHOB],P_orig=U[PRESSURE],vx_orig=U[VX],vy_orig=U[VY],vz_orig=U[VZ];\n\n  /***********************************************************/\n  // Enforce limits on pressure, density, and v^i\n  /***********************************************************/\n  // Density floor:\n  U[RHOB] = MAX(U[RHOB],rho_b_atm);\n  // Density ceiling:\n  U[RHOB] = MIN(U[RHOB],rho_b_max);\n\n  //   Next set h, the enthalpy:\n  CCTK_REAL h_enthalpy,  P_cold,eps_cold,dPcold_drho,eps_th,Gamma_cold; /* <- Note that in setting h, we need to define several\n                                                                         *    other variables. Though some will be unused later\n                                                                         *    in this function, they may be useful in other\n                                                                         *    functions */\n  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);\n\n  // Pressure floor & ceiling:\n  int polytropic_index = find_polytropic_K_and_Gamma_index(eos,U[RHOB]);\n  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]);\n\n  // Possibly adjusted pressure, so recompute eps & h:\n  CCTK_REAL eps = eps_cold + (U[PRESSURE]-P_cold)/(Gamma_th-1.0)/U[RHOB];\n  h_enthalpy = 1.0 + eps + U[PRESSURE]/U[RHOB];\n\n  CCTK_REAL uUP[4];\n  impose_speed_limit_output_u0(METRIC,U,METRIC_LAP_PSI4[PSI4],METRIC_LAP_PSI4[LAPSEINV],stats, uUP[0]);\n  // Compute u^i. We've already set uUP[0] in the lines above.\n  for(int ii=0;ii<3;ii++) uUP[UX+ii] = uUP[0]*U[VX+ii];\n\n  // Useful debugging tool, part 2: can be used to track fixes:\n  //if(P_orig!=U[PRESSURE] || rho_b_orig!=U[RHOB] || vx_orig!=U[VX] || vy_orig!=U[VY] || vz_orig!=U[VZ]) {\n")


# <a id='bi_b2_and_u_i'></a>
# 
# ## Step 3.b: Computing $b^{i}$, $b^{2}$, and $u_{i}$ \[Back to [top](#toc)\]
# $$\label{bi_b2_and_u_i}$$
# 
# 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-IllinoisGRMHD__inlined_functions.ipynb) tutorial module. We also compute $u_{i}$ from the output of the `compute_smallba_b2_and_u_i_over_u0_psi4()` function.

# In[13]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  /***************************************************************/\n  // COMPUTE TUPMUNU, TDNMUNU, AND  CONSERVATIVES FROM PRIMITIVES\n  /***************************************************************/\n  // Compute b^{\\mu}, b^2, and u_i/(u^0 Psi4)\n  CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI = METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;\n  CCTK_REAL u_x_over_u0_psi4,u_y_over_u0_psi4,u_z_over_u0_psi4;\n  CCTK_REAL smallb[NUMVARS_SMALLB];\n  compute_smallba_b2_and_u_i_over_u0_psi4(METRIC,METRIC_LAP_PSI4,U,uUP[0],ONE_OVER_LAPSE_SQRT_4PI,\n                                          u_x_over_u0_psi4,u_y_over_u0_psi4,u_z_over_u0_psi4,smallb);\n  // Compute u_i; we compute u_0 below.\n  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] };\n')


# <a id='gupmunu'></a>
# 
# ## Step 3.c: The physical ADM inverse 4-metric, $g^{\mu\nu}$ \[Back to [top](#toc)\]
# $$\label{gupmunu}$$
# 
# Then we start evaluating the physical ADM inverse 4-metric, $g^{\mu\nu}$ (eq. 2.119 in [Numerical Relativity, by Baumgarte & Shapiro](https://books.google.com/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC))
# 
# $$
# \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}
# }\ .
# $$

# In[14]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  // Precompute some useful quantities, for later:\n  CCTK_REAL alpha_sqrt_gamma=METRIC_LAP_PSI4[LAPSE]*METRIC_LAP_PSI4[PSI6];\n  CCTK_REAL rho0_h_plus_b2 = (U[RHOB]*h_enthalpy + smallb[SMALLB2]);\n  CCTK_REAL P_plus_half_b2 = (U[PRESSURE]+0.5*smallb[SMALLB2]);\n\n  if(!already_computed_physical_metric_and_inverse) {\n    // If you like, see Eq 2.119 in Numerical Relativity, by Baumgarte & Shapiro:\n    CCTK_REAL ONE_OVER_LAPSE_SQUARED = SQR(METRIC_LAP_PSI4[LAPSEINV]);\n\n    // g^{\\mu \\nu} = upper four-metric.\n    //CCTK_REAL g4up[4][4];\n    g4up[0][0] = -ONE_OVER_LAPSE_SQUARED;\n    g4up[0][1] = ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX];\n    g4up[0][2] = ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTY];\n    g4up[0][3] = ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTZ];\n    g4up[1][1] = METRIC[GUPXX]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX]*METRIC[SHIFTX];\n    g4up[1][2] = METRIC[GUPXY]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX]*METRIC[SHIFTY];\n    g4up[1][3] = METRIC[GUPXZ]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTX]*METRIC[SHIFTZ];\n    g4up[2][2] = METRIC[GUPYY]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTY]*METRIC[SHIFTY];\n    g4up[2][3] = METRIC[GUPYZ]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTY]*METRIC[SHIFTZ];\n    g4up[3][3] = METRIC[GUPZZ]*METRIC_LAP_PSI4[PSIM4] - ONE_OVER_LAPSE_SQUARED*METRIC[SHIFTZ]*METRIC[SHIFTZ];\n  }\n')


# <a id='tupmunu'></a>
# 
# ## Step 3.d: The stress-energy tensor, $T^{\mu\nu}$ \[Back to [top](#toc)\]
# $$\label{tupmunu}$$
# 
# Then we compute the stress-energy tensor, $T^{\mu\nu}$, which is given by eq. (8) in [Etienne *et al.* (2015)](https://arxiv.org/pdf/1501.07276.pdf):
# 
# $$
# \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}
# }\ .
# $$

# In[15]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  int count;\n  // Next compute T^{\\mu \\nu}. This is very useful when computing fluxes and source terms in the GRMHD evolution equations.\n  // (Eq. 33 in http://arxiv.org/pdf/astro-ph/0503420.pdf):\n  // 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.\n  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++; }\n')


# <a id='gdnmunu'></a>
# 
# ## Step 3.e: The physical ADM 4-metric, $g_{\mu\nu}$ \[Back to [top](#toc)\]
# $$\label{gdnmunu}$$
# 
# 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](https://books.google.com/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC))
# 
# $$
# \boxed{g_{\mu\nu} =
# \begin{pmatrix}
# -\alpha^{2} + \beta_{\ell}\beta^{\ell} & \beta_{i}\\
# \beta_{j} & \gamma_{ij}
# \end{pmatrix}}\ .
# $$

# In[16]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  // Next compute T_{\\mu \\nu}\n  // 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.\n  if(!already_computed_physical_metric_and_inverse) {\n    CCTK_REAL LAPSE_SQUARED=SQR(METRIC_LAP_PSI4[LAPSE]);\n    CCTK_REAL BETADN[4],BETAUP[4] = { 0.0, METRIC[SHIFTX],METRIC[SHIFTY],METRIC[SHIFTZ] };\n    lower_4vector_output_spatial_part(METRIC_LAP_PSI4[PSI4],METRIC,BETAUP, BETADN);\n\n    //CCTK_REAL g4dn[4][4];\n    // g_{00} = - alpha^2 + gamma_{ij} beta^i beta^j = - alpha^2 beta_i beta^i\n    g4dn[0][0] = -LAPSE_SQUARED + (BETAUP[1]*BETADN[1] + BETAUP[2]*BETADN[2] + BETAUP[3]*BETADN[3]);\n    // g_{0i} =  gamma_{ij} beta^j = beta_i\n    g4dn[0][1] = g4dn[1][0] = BETADN[1];\n    g4dn[0][2] = g4dn[2][0] = BETADN[2];\n    g4dn[0][3] = g4dn[3][0] = BETADN[3];\n    // g_{ij} =  gamma_{ij} = Psi^4 \\tilde{gamma_ij}\n    g4dn[1][1] =              METRIC[GXX]*METRIC_LAP_PSI4[PSI4];\n    g4dn[1][2] = g4dn[2][1] = METRIC[GXY]*METRIC_LAP_PSI4[PSI4];\n    g4dn[1][3] = g4dn[3][1] = METRIC[GXZ]*METRIC_LAP_PSI4[PSI4];\n    g4dn[2][2] =              METRIC[GYY]*METRIC_LAP_PSI4[PSI4];\n    g4dn[2][3] = g4dn[3][2] = METRIC[GYZ]*METRIC_LAP_PSI4[PSI4];\n    g4dn[3][3] =              METRIC[GZZ]*METRIC_LAP_PSI4[PSI4];\n  }\n')


# <a id='tdnmunu'></a>
# 
# ## Step 3.f: The stress-energy tensor, $T_{\mu\nu}$ \[Back to [top](#toc)\]
# $$\label{tdnmunu}$$
# 
# Then we compute $b_{\mu} = g_{\mu\nu}b^{\nu}$ and $u_{0} = g_{0\mu}u^{\mu}$.

# In[17]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', "\n\n  CCTK_REAL smallb_lower[NUMVARS_SMALLB];\n  // FIXME: This could be replaced by the function call\n  //           lower_4vector_output_spatial_part(psi4,smallb,smallb_lower);\n  // b_a = b^c g_{ac}\n  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]; }\n\n  // Compute u_0, as we've already computed u_i above.\n  uDN[0]=0.0; for(int jj=0;jj<4;jj++) uDN[0] += uUP[jj]*g4dn[0][jj];\n")


# 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}
# }\ .
# $$

# In[18]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  // Compute T_{\\mu \\nu}\n  if(update_Tmunu) {\n    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++; }\n  }\n')


# <a id='recompute_conservs'></a>
# 
# ## Step 3.g: Recomputing the conservative variables \[Back to [top](#toc)\]
# $$\label{recompute_conservs}$$
# 
# 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.
# }\ .
# $$

# In[19]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__tau_and_prim_limits__C', '\n\n  //CCTK_VInfo(CCTK_THORNSTRING,"YAY %e",TDNMUNU[0]);\n\n  // Finally, compute conservatives:\n  CONSERVS[RHOSTAR] = alpha_sqrt_gamma * U[RHOB] * uUP[0];\n  CONSERVS[STILDEX] = CONSERVS[RHOSTAR]*h_enthalpy*uDN[1] + alpha_sqrt_gamma*(uUP[0]*smallb[SMALLB2]*uDN[1] - smallb[SMALLBT]*smallb_lower[SMALLBX]);\n  CONSERVS[STILDEY] = CONSERVS[RHOSTAR]*h_enthalpy*uDN[2] + alpha_sqrt_gamma*(uUP[0]*smallb[SMALLB2]*uDN[2] - smallb[SMALLBT]*smallb_lower[SMALLBY]);\n  CONSERVS[STILDEZ] = CONSERVS[RHOSTAR]*h_enthalpy*uDN[3] + alpha_sqrt_gamma*(uUP[0]*smallb[SMALLB2]*uDN[3] - smallb[SMALLBT]*smallb_lower[SMALLBZ]);\n  // tauL = alpha^2 sqrt(gamma) T^{00} - CONSERVS[RHOSTAR]\n  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];\n  // FIXME: We have TUPMUNU already, should replace           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ by TUPMUNU[0].\n}\n\n\n')


# <a id='code_validation'></a>
# 
# # Step 4: Code validation \[Back to [top](#toc)\]
# $$\label{code_validation}$$
# 
# First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook.

# In[20]:


# 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
        get_ipython().system('wget -O $original_IGM_file_path $original_IGM_file_url')

# Perform validation
Validation__tau_and_prims_limits__C  = get_ipython().getoutput('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
    get_ipython().system('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)


# <a id='latex_pdf_output'></a>
# 
# # Step 5: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
# $$\label{latex_pdf_output}$$
# 
# 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](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).

# In[21]:


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
get_ipython().system('rm -f Tut*.out Tut*.aux Tut*.log')