{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# `GiRaFFE_NRPy`: Basic Equations\n",
    "\n",
    "## Authors: Patrick Nelson &\n",
    "\n",
    "In this tutorial, we introduce the basic GRFFE equations that we wish to solve. \n",
    "\n",
    "Throughout this section and beyond, we will work in geometrized units such that $G=c=1$. As is standard in NRPy+, \n",
    "* Greek indices refer to four-dimensional quantities where the zeroth component indicates temporal (time) component,\n",
    "* and Latin indices refer to three-dimensional quantities. This is somewhat counterintuitive since Python always indexes its lists starting from 0. As a result, the zeroth component of three-dimensional quantities will necessarily indicate the first *spatial* direction.\n",
    "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. For further detail and examples, see [this notebook](../Tutorial-GRFFE_Equations-Cartesian.ipynb). \n",
    "\n",
    "To describe spacetime, we will use a standard 3+1 decomposition in which the line element is given as \n",
    "$$\n",
    "ds^2 = -\\alpha^2 dt^2 + \\gamma_{ij} \\left( dx^i + \\beta^i dt \\right) \\left( dx^j + \\beta^j dt \\right).\n",
    "$$\n",
    "where $\\alpha$ is the lapse function, $\\beta^i$ is the shift vector, $\\gamma_{ij}$ is the spatial three metric. Given two adjacent spatial hypersurfaces with coordinate times $t_0$ and $t_0+dt$, $\\alpha dt$ is the proper time interval between the two hypersurfaces, $\\beta^i$ represents the magnitude of the spatial coordinate shift between them, and $\\gamma_{ij}$ is the spatial three metric within a hypersurface.\n",
    "\n",
    "We will use the pure electromagnetic stress energy tensor, \n",
    "$$\n",
    "T^{\\mu\\nu}_{\\rm EM} = b^2 u^{\\mu} u^{\\nu} + \\frac{1}{2} b^2 g^{\\mu\\nu} - b^\\mu b^\\nu,\n",
    "$$\n",
    "where \n",
    "\\begin{align}\n",
    "\\sqrt{4\\pi} b^0 = B^0_{\\rm (u)} &= \\frac{u_j B^j}{\\alpha} \\\\\n",
    "\\sqrt{4\\pi} b^i = B^i_{\\rm (u)} &= \\frac{B^i + (u_j B^j) u^i}{\\alpha u^0},\\\\\n",
    "\\end{align}\n",
    "and $u^\\mu$ is the four-velocity. We will not couple this stress-energy tensor to the Einstein Tensor because the plasma we consider is extremely diffuse. \n",
    "\n",
    "**Derivation from $F^{\\mu\\nu}$, including the constraints of GRFFE**\n",
    "\n",
    "In place of the magnetic field $B^i$, we choose to evolve the vector potential $A_i$. By doing so, we guarantee that the magnetic field is divergenceless to round-off level. We also evolve the Poynting flux $\\tilde{S}_i$, which in GRFFE, represents the momentum flux in the plasma, and $\\psi^6 \\Phi$, where $\\psi^6=\\sqrt{\\gamma}$ and $\\gamma$ is the determinant of the three metric.\n",
    "\n",
    "We thus arrive at the following set of equations to evolve:\n",
    "\\begin{align}\n",
    "\\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} \\partial_i g_{\\mu \\nu} \\\\\n",
    "\\partial_t [\\sqrt{\\gamma} \\Phi] &= -\\partial_j (\\alpha\\sqrt{\\gamma}A^j - \\beta^j [\\sqrt{\\gamma} \\Phi]) - \\xi \\alpha [\\sqrt{\\gamma} \\Phi] \\\\\n",
    "\\partial_t A_i &= \\epsilon_{ijk} v^j B^k - \\partial_i (\\alpha \\Phi - \\beta^j A_j) \\\\\n",
    "\\end{align}\n",
    "\n",
    "Here, the four velocity $u^\\mu$ is the velocity of the plasma as measured by a normal observer. There are several possible choices of three-velocity that can be made; the ones which we will use here are the drift velocity $v^i$ and the Valencia three-velocity $\\bar{v}^i$. While we have chosen to use the Valencia three-velocity in this version of `GiRaFFE`, we have also frequently made use of a relationship expressing this in terms of the drift velocity, which was used in the original `GiRaFFE`. The usefulness of this relationship to drift velocity extends beyond merely translating the original code. As discussed in [Paschalidis, et al.](https://arxiv.org/pdf/1310.3274.pdf), Sec. III.A (just above Eq. 45, with a proof in Appendix A), there is a one-parameter family of velocity definitions that fulfill the GRFFE conditions. The drift velocity sets this parameter to 0, which minimizes the Lorentz factor and *guarantees* that the four-velocity and magnetic fields are orthogonal to each other. This simplifies the form of $b^\\mu$ and quantities that depend on it. \n",
    "\n",
    "This must be taken into account in developing unit tests, because NRPy+'s GRFFE module defaults to using a definition of $b^\\mu$ that does not assume that this criterion is met, while the original `GiRaFFE` code assumes this in its C2P and P2C solvers. So, if we do not guarantee that our test data fulfills this criterion, these two different routines will produce different results. We will now go through the derivation of the equation used by `GiRaFFE` from first principles to show where this extra term appears.\n",
    "\n",
    "<font color='yellow'><b>Terrence says: What extra term?</b></font>\n",
    "\n",
    "This is the equation used by `GiRaFFE`, pulled from Eqs. 47 and 85 of [this paper](https://arxiv.org/abs/1310.3274):\n",
    "$$\\tilde{S}_i = \\gamma_{ij} \\frac{\\bar{v}^j \\sqrt{\\gamma}B^2}{4 \\pi},$$\n",
    "or $$\\tilde{S}_i = \\gamma_{ij} \\frac{(v^j + \\beta^j) \\sqrt{\\gamma}B^2}{4 \\pi \\alpha},$$\n",
    "where $\\bar{v}^j$ is the Valencia 3-velocity and $v^j$ is the drift velocity.\n",
    "\n",
    "In IllinoisGRMHD (IGM), the expression used is $\\tilde{S}_i = \\alpha \\sqrt{\\gamma} T^0_{{\\rm EM}i},$\n",
    "where \n",
    "\\begin{align}\n",
    "T^{\\mu\\nu}_{\\rm EM} &= b^2 u^{\\mu} u^{\\nu} + \\frac{1}{2} b^2 g^{\\mu\\nu} - b^\\mu b^\\nu \\\\\n",
    "b^0 &= \\frac{u_j B^j}{\\sqrt{4\\pi} \\alpha} \\\\\n",
    "b^i &= \\frac{B^i + (u_j B^j) u^i}{\\sqrt{4\\pi} \\alpha u^0} \\\\\n",
    "u^i &= u^0 v^i.\n",
    "\\end{align}\n",
    "**The above has been pulled from the C2P unit test. This is probably a better permanent home for this discussion.**\n",
    "\n",
    "<font color='yellow'><b>Terrence says: Can we remove the above?</b></font>\n",
    "\n",
    "\n",
    "We also use the Valencia three-velocity $\\bar{v}^i$ and the magnetic field $B^i$ as primitives, which can be defined in terms of the conservative variables:\n",
    "\\begin{align}\n",
    "B^i &= \\epsilon^{ijk} \\partial_j A_k \\\\\n",
    "\\bar{v}^i &= 4 \\pi \\frac{\\gamma^{ij} {\\tilde S}_j}{\\sqrt{\\gamma} B^2} \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os,sys\n",
    "nrpy_dir_path = os.path.join(\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "\n",
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"GiRaFFE_NRPy_Tutorial\",location_of_template_file=os.path.join(\"..\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}