#!/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>
# 
# # IllinoisGRMHD: Basic equations and modules
# 
# ## Authors: Leo Werneck & Zach Etienne
# 
# <font color='red'>**This module is currently under development**</font>
# 
# ## This module introduces the basic equations solved by IllinoisGRMHD and provides an introduction to the modules inside the code.
# 
# ## Introduction:
# 
# [`IllinoisGRMHD`](http://arxiv.org/abs/1501.07276) solves the equations of General Relativistic MagnetoHydroDynamics (GRMHD) using a high-resolution shock capturing scheme. It is a rewrite of the Illinois Numerical Relativity (ILNR) group's GRMHD code, and generates results that agree to roundoff error with that original code. Its feature set coincides with the features of the ILNR group's recent code (ca. 2009-2014), which was used in their modeling of the following systems:
# 
# 1. Magnetized circumbinary disk accretion onto binary black holes
# 2. Magnetized black hole-neutron star mergers
# 3. Magnetized Bondi flow, Bondi-Hoyle-Littleton accretion
# 4. White dwarf-neutron star mergers
# 
# `IllinoisGRMHD` is particularly good at modeling GRMHD flows into black holes without the need for excision. Its [HARM-based conservative-to-primitive solver](https://arxiv.org/abs/astro-ph/0512420) has also been modified to check the physicality of conservative variables prior to primitive inversion and moves them into the physical range if they become unphysical.
# 
# Currently, IllinoisGRMHD consists of
# 
# 1. the Piecewise Parabolic Method (PPM) for reconstruction, 
# 2. the Harten, Lax, van Leer (HLL/HLLE) approximate Riemann solver, and
# 3. a modified HARM Conservative-to-Primitive solver. 
# 
# `IllinoisGRMHD` evolves the vector potential $A_{\mu}$ (on staggered grids) instead of the magnetic fields ($B^i$) directly, to guarantee that the magnetic fields will remain divergenceless even at AMR boundaries. On uniform resolution grids, this vector potential formulation produces results equivalent to those generated using the standard, staggered flux-CT scheme. This scheme is based on that of [Del Zanna *et al.* (2003)](https://arxiv.org/abs/astro-ph/0210618).
# 
# ### 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
# 
# 1. [Step 1](#basic_equations): **The equations of GRMHD**
# 1. [Step 2](#suitable_equations): **Recasting the equations of GRMHD into a more useful form**
#     1. [Step 2.1](#einstein_equations): Einstein's equations
#     1. [Step 2.2](#conserv_of_baryon_number): Conservation of baryon number
#     1. [Step 2.3](#conserv_of_em_tensor): Conservation of energy-momentum
#     1. [Step 2.4](#em_evo_equations): The evolution equations for the electromagnetic field
#     1. [Step 2.5](#equation_of_state): The equation of state (EOS)
# 1. [Step 3](#latex_pdf_output): **Output this notebook to $\LaTeX$-formatted PDF file**

# <a id='basic_equations'></a>
# 
# # Step 1: The equations of GRMHD \[Back to [top](#toc)\]
# $$\label{basic_equations}$$
# 
# The basic equations solved by `IllinoisGRMHD` are
# 
# 1. Einstein's equations: $R^{\mu\nu} + \frac{1}{2}g^{\mu\nu}R = 8\pi T^{\mu\nu} = 8\pi \left(T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\right)$ ;
# 1. Conservation of baryon number: $\nabla_{\mu}\left(\rho_{0} u^{\mu}\right) = 0$ ;
# 1. Conservation of energy-momentum: $\nabla_{\mu}T^{\mu\nu} = \nabla_{\mu}\left(T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\right) = 0$ ;
# 1. Dual of homogeneous Maxwell's equations: $\nabla_{\nu} F^{*\mu\nu}=\frac{1}{\sqrt{-g}}\partial_{\nu}\left(\sqrt{-g}F^{*\mu\nu}\right)=0$ ,
# 
# where $g_{\mu\nu}$ is the ADM 4-metric, $R_{\mu\nu}$ and $R$ are the Ricci tensor and scalar, respectively, constructed from $g_{\mu\nu}$, $\nabla_{\mu}$ is the covariant derivative compatible with $g_{\mu\nu}$, $g\equiv\det\left(g_{\mu\nu}\right)$, $\rho_0$ is the rest-mass density of the fluid, $u^{\mu}$ is the four-velocity of the fluid, $F^{\mu\nu}$ is the Faraday tensor and $F^{*\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}$ its dual ($\epsilon^{\mu\nu\rho\sigma}$ is the Levi-Civita symbol).
# 
# The final equation that must be solved is the equation of state (EOS) for the matter. Currently `IllinoisGRMHD` implements a [hybrid EOS of the form](https://ui.adsabs.harvard.edu/abs/1993A&A...268..360J/abstract)
# 
# 5. Equation of State (EOS): $P(\rho_{0},\epsilon) = P_{\rm cold}(\rho_{0}) + \left(\Gamma_{\rm th}-1\right) \rho_{0} \left[\epsilon-\epsilon_{\rm cold}(\rho_{0})\right]$ ,
# 
# where P is the pressure, $\epsilon$ the specific internal energy (the subscript cold indicate the cold components of these quantities), and $\Gamma_{\rm th}$ is a constant parameter which determines the conversion efficiency of kinetic to thermal energy at shocks.
# 
# <a id='suitable_equations'></a>
# 
# # Step 2: Recasting the equations of GRMHD into a more useful form \[Back to [top](#toc)\]
# $$\label{suitable_equations}$$
# 
# In this step, we will write down the equations used by `IllinoisGRMHD` in the form they are implemented. To give an example, the GRMHD equations are written in the *conservative form*
# 
# $$
# \partial_{t}\boldsymbol{C} + \boldsymbol{\nabla}\cdot\boldsymbol{F} = \boldsymbol{S}\ ,
# $$
# 
# where $\boldsymbol{C} = \left\{\rho_{\star},\tilde\tau,\tilde{S}_{i},\tilde{B}^{i}\right\}$ is the vector of conservative variables, $\boldsymbol{F}$ is the flux vector, and $\boldsymbol{S}$ the vector of source terms (explicit expressions for the components of $\boldsymbol{C}$, $\boldsymbol{F}$, and $\boldsymbol{S}$ can be found below).
# 
# <a id='einstein_equations'></a>
# 
# ## Step 2.1: Einstein's equations \[Back to [top](#toc)\]
# $$\label{einstein_equations}$$
# 
# `IllinoisGRMHD` solves Einstein's field equations (with $G_{\rm N}=1=c$) in the presence of matter sources,
# 
# $$
# G^{\mu\nu} = 8\pi T^{\mu\nu}\ ,
# $$
# 
# where $G^{\mu\nu} \equiv R^{\mu\nu} + \frac{1}{2}g^{\mu\nu}R$ is the Einstein tensor and the total energy-momentum tensor, $T^{\mu\nu}$, is the sum of the matter and electromagnetic energy-momentum tensors:
# 
# $$
# T^{\mu\nu} = T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\ .
# $$
# 
# Einstein's field equations are solved using the BSSN formalism (see [this tutorial module](../../Tutorial-BSSN_formulation.ipynb) for an overview).
# 
# <a id='conserv_of_baryon_number'></a>
# 
# ## Step 2.2: Conservation of baryon number \[Back to [top](#toc)\]
# $$\label{conserv_of_baryon_number}$$
# 
# The conservation of the baryon number equation is written as (cf. eqs. (6), (7), and (18) in [Etienne *et al.*](https://arxiv.org/pdf/1501.07276.pdf))
# 
# $$
# \boxed{\partial_{t}\rho_{\star} + \partial_{j}\left(\rho_{\star}v^{j}\right) = 0}\ ,
# $$
# 
# where $\rho_{\star}\equiv \alpha \sqrt{\gamma}\rho_{0}u^{0}$ and $v^{i}\equiv u^{i}/u^{0}$.
# 
# <a id='conserv_of_em_tensor'></a>
# 
# ## Step 2.3: Conservation of energy-momentum \[Back to [top](#toc)\]
# $$\label{conserv_of_em_tensor}$$
# 
# In the ideal MHD limit, we can write down the total energy-momentum tensor as (cf. eq. (8) in [Etienne *et al.*](https://arxiv.org/pdf/1501.07276.pdf) and the discussion before it)
# 
# $$
# T^{\mu \nu} = (\rho_0 h +b^2) u^{\mu} u^{\nu} + \left( P + \frac{b^2}{2}\right) g^{\mu \nu} - b^{\mu} b^{\nu}
# $$
# 
# The spatial components of the energy-momentum conservation equation give (cf. eq. (18) of [Etienne *et al.*](https://arxiv.org/pdf/1501.07276.pdf) and Eqs. (35) and (36) of [Duez et al.](https://arxiv.org/pdf/astro-ph/0503420.pdf))
# 
# $$
# \boxed{\partial_{t}\tilde{S}_{i} + \partial_{j}\left(\alpha\sqrt{\gamma} T^{j}_{\ i}\right) = \frac{1}{2}\sqrt{\gamma} T^{\alpha\beta}\partial_{i}g_{\alpha\beta}}\ ,
# $$
# 
# where
# 
# $$
# \tilde{S}_{i} = \sqrt{\gamma}S_{i} = \alpha \sqrt{\gamma} T^{0}_{\ i} = \left(\rho_{\star}h + \alpha u^{0}\sqrt{\gamma}b^{2}\right)u_{i} - \alpha \sqrt{\gamma} b^{0} b^{i}\ .
# $$
# 
# The time component of the energy-momentum conservation equation gives (cf. eq. (18) of [Etienne *et al.*](https://arxiv.org/pdf/1501.07276.pdf))
# 
# $$
# \boxed{\partial_{t}\tilde{\tau} + \partial_{i}\left(\alpha^{2}\sqrt{\gamma}T^{0i} - \rho_{\star}v^{i}\right) = s}\ ,
# $$
# 
# where
# 
# \begin{align}
# \tilde{\tau} &= \sqrt{\gamma}n_{\mu}n_{\nu} - \rho_{\star} = \alpha^{2}\sqrt{\gamma}T^{00} - \rho_{\star}\ ,\\
# s &= -\alpha \sqrt{\gamma} T^{\mu\nu} \nabla_{\nu} n_{\mu} = \alpha\sqrt{\gamma}\left[\left(T^{00}\beta^{i}\beta^{j} + 2T^{0i}\beta^{j} + T^{ij}\right)K_{ij} - \left(T^{00}\beta^{i} + T^{0i}\right)\partial_{i}\alpha\right]\ ,
# \end{align}
# 
# with $n_{\mu} = \left(\alpha,0,0,0\right)$ being the normal vector and $K_{ij}$ the extrinsic curvature.
# 
# <a id='em_evo_equations'></a>
# 
# ## Step 2.4: The evolution equations for the electromagnetic field \[Back to [top](#toc)\]
# $$\label{em_evo_equations}$$
# 
# From the spatial components of the dual of Maxwell's equations,
# 
# $$
# \nabla_{\nu} F^{*\mu\nu}=\frac{1}{\sqrt{-g}}\partial_{\nu}\left(\sqrt{-g}F^{*\mu\nu}\right)=0\ ,
# $$
# 
# we get the magnetic induction equation, which in conservative form may be written as (cf. eq. (12) [Etienne *et al.*](https://arxiv.org/pdf/1501.07276.pdf))
# 
# $$
# \partial_{t}\tilde{B}^{i} + \partial_{j}\left(v^{j}\tilde{B}^{i} - v^{i}\tilde{B}^{j}\right) = 0\ ,
# $$
# 
# where $\tilde{B}^{i} = \sqrt{\gamma}B^{i}$. We must also guarantee that no magnetic monopoles form, via the constraint
# 
# $$
# \partial_{i}\tilde{B}^{i} = 0\ ,
# $$
# 
# which is the time component of the dual of Maxwell's equations. Satisfying this constraint equation while evolving the magnetic field forward in time via the evolution equation above turns out to be a nontrivial endeavor, *particularly* on AMR grids. Instead, we choose to evolve the magnetic 4-vector potential $\mathcal{A}_{\mu}$ instead of the magnetic fields directly, so that
# 
# \begin{align}
# \mathcal{A}_{\mu} &= \Phi n_{\mu} + A_{\mu}\ ,\\
# \tilde{B}^{i} &= \tilde{\epsilon}^{ijk}\partial_{j}A_{k}\ ,
# \end{align}
# 
# where $A_{\mu}$ is purely spatial ($A_{\mu}n^{\mu} = 0$) and $\Phi$ is the EM scalar potential.
# 
# Special finite difference operators for the vector potential are defined in IllinoisGRMHD so that the divergence of a curl is zero to roundoff error, which implies that the divergence of $\tilde{B}^{i}$ (as defined above) is zero and the condition $\partial_{i}\tilde{B}^{i}=0$ is satisfied automatically, *even on AMR grids*.
# 
# In terms of $A_{i}$, the induction equation becomes
# 
# $$
# \boxed{\partial_t A_i = \tilde\epsilon_{ijk} v^j \tilde{B}^k - \partial_i (\alpha \Phi - \beta^j A_j)}\ .
# $$
# 
# The final equation comes from choosing the covariant version of the "[generalized Lorenz gauge condition](https://arxiv.org/pdf/1207.3354.pdf)" that was introduced by the Illinois Relativity group,
# 
# $$
# \nabla_{\mu}\mathcal{A}^{\mu} = \xi n_{\mu} \mathcal{A}^{\mu}\ ,
# $$
# 
# where $\xi$ is a parameter with dimensions 1/Length, chosen so that the CFL condition remains satisfied. This gauge choice yields the additional evolution equation
# 
# $$
# \boxed{\partial_t [\sqrt{\gamma} \Phi] + \partial_j (\alpha \sqrt{\gamma} A^j - \beta^j [\sqrt{\gamma}\Phi]) = -\xi\alpha \sqrt{\gamma}\Phi}\ .
# $$
# 
# Note that in `IllinoisGRMHD` the evolved variable is $\sqrt{\gamma}\Phi$, not $\Phi$.
# 
# <a id='equation_of_state'></a>
# 
# ## Step 2.5: The equation of state (EOS) \[Back to [top](#toc)\]
# $$\label{equation_of_state}$$
# 
# `IllinoisGRMHD` currently implements a hybrid EOS of the form
# 
# $$
# \boxed{P(\rho_{0},\epsilon) = P_{\rm cold}(\rho_{0}) + \left(\Gamma_{\rm th}-1\right) \rho_{0} \left[\epsilon-\epsilon_{\rm cold}(\rho_{0})\right]}\ .
# $$
# 
# The function $\epsilon_{\rm cold}$ is related to $P_{\rm cold}$ by the first law of thermodynamics,
# 
# $$
# \epsilon_{\rm cold}(\rho_{0}) = \int \frac{P_{\rm cold}(\rho_{0})}{\rho_{0}^{2}}d\rho_{0}\ .
# $$
# 
# Currently, all functions within `IllinoisGRMHD` support piecewise-defined $P_{\rm cold}(\rho_{0})$ (the so-called "piecewise polytrope" EOS) with up to nine different polytropic indices.
# 
# In typical runs of the code, particularly on the ones presented in the [release paper](https://arxiv.org/pdf/1501.07276.pdf), the $\Gamma$-law EOS, $P = (\Gamma - 1)\rho_{0}\epsilon$, is adopted. This corresponds to setting $P_{\rm cold} = (\Gamma - 1)\rho_{0}\epsilon_{\rm cold}$ in the boxed equation above, which is equivalent to $P_{\rm cold} = \kappa\rho_{0}^{\Gamma}$ (with constant $\kappa$), and $\Gamma_{\rm th} = \Gamma$. In the absence of shocks, $\epsilon = \epsilon_{\rm cold}$ so that $P = P_{\rm cold}$.
# 
# `IllinoisGRMHD` has also performed successful binary neutron star mergers using the widely adopted piecewise polytropic EOSs of [Read *et al.*](https://arxiv.org/abs/0812.2163). In the [EOS low-level functions tutorial](Tutorial-IllinoisGRMHD__EoS_lowlevel_functs.ipynb), you will find special `Python` functions designed to help the user who wishes to use [Read *et al.*](https://arxiv.org/abs/0812.2163) EOSs or any other piecewise polytropic EOS.

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

# In[1]:


#!jupyter nbconvert --to latex --template ../../latex_nrpy_style.tplx Tutorial-IllinoisGRMHD__Overview.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Overview.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Overview.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Overview.tex
get_ipython().system('rm -f Tut*.out Tut*.aux Tut*.log')