Tutorial-GRHD_Equations-Cartesian-c-code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# C code generation of GRHD Equations\n",
    "\n",
    "## Authors: Phil Chang & Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This notebook demonstrates the C code generation of General Relativistic HydroDynamics (GRHD) equations in conservative form using a specific state vector ${\\boldsymbol{\\mathcal{U}}}=(\\rho_*,\\tilde{S},\\tilde{\\tau})$. The process of transitioning between primitive and conservative variables, computing flux and source terms, as well as performing Lorentz boosts and other transformations is detailed.\n",
    "[comment]: <> (omit white space somehow?)\n",
    "$\\newcommand{\\be}{\\begin{equation}}$\n",
    "$\\newcommand{\\ee}{\\end{equation}}$\n",
    "$\\newcommand{\\grad}{{\\boldsymbol{\\nabla}}}$\n",
    "$\\newcommand{\\vel}{{\\boldsymbol{v}}}$\n",
    "$\\newcommand{\\mom}{{\\boldsymbol{p}}}$\n",
    "$\\newcommand{\\ddt}[1]{{\\frac{\\partial #1}{\\partial t}}}$\n",
    "$\\newcommand{\\ddx}[1]{{\\frac{\\partial #1}{\\partial x}}}$\n",
    "$\\newcommand{\\state}{{\\boldsymbol{\\mathcal{U}}}}$\n",
    "$\\newcommand{\\charge}{{\\boldsymbol{U}}}$\n",
    "$\\newcommand{\\psicharge}{{\\boldsymbol{\\psi}}}$\n",
    "$\\newcommand{\\lapse}{\\alpha}$\n",
    "$\\newcommand{\\shift}{\\boldsymbol{\\beta}}$\n",
    "$\\newcommand{\\rhostar}{{\\rho_*}}$\n",
    "$\\newcommand{\\tautilde}{{\\tilde{\\tau}}}$\n",
    "$\\newcommand{\\Svectilde}{{\\tilde{\\boldsymbol{S}}}}$\n",
    "$\\newcommand{\\rtgamma}{{\\sqrt{\\gamma}}}$\n",
    "$\\newcommand{\\T}[2]{{T^{#1 #2}}}$\n",
    "$\\newcommand{\\uvec}{{\\boldsymbol{u}}}$\n",
    "$\\newcommand{\\Vvec}{{\\boldsymbol{\\mathcal{V}}}}$\n",
    "$\\newcommand{\\vfluid}{{\\boldsymbol{v}_{\\rm f}}}$\n",
    "$\\newcommand{\\vVal}{{\\tilde{\\boldsymbol{v}}}}$\n",
    "\n",
    "$\\newcommand{\\flux}{{\\boldsymbol{\\mathcal{F}}}}$\n",
    "$\\newcommand{\\fluxV}{{\\boldsymbol{F}}}$\n",
    "$\\newcommand{\\source}{{\\boldsymbol{\\mathcal{S}}}}$\n",
    "$\\newcommand{\\sourceV}{{\\boldsymbol{S}}}$\n",
    "\n",
    "$\\newcommand{\\area}{{\\boldsymbol{A}}}$\n",
    "$\\newcommand{\\normal}{{\\hat{\\boldsymbol{n}}}}$\n",
    "$\\newcommand{\\pt}{{\\boldsymbol{p}}}$\n",
    "$\\newcommand{\\nb}{{\\boldsymbol{n}}}$\n",
    "$\\newcommand{\\meshv}{{\\boldsymbol{w}}}$\n",
    "$\\newcommand{\\facev}{{\\boldsymbol{\\tilde{w}}_{ij}}}$\n",
    "$\\newcommand{\\facer}{{\\boldsymbol{\\tilde{r}}_{ij}}}$\n",
    "$\\newcommand{\\meshr}{{\\boldsymbol{r}}}$\n",
    "$\\newcommand{\\cmr}{{\\boldsymbol{c}}}$\n",
    "\n",
    "## Introduction: \n",
    "We start out with the ** GRHD ** equations in conservative form with the state vector $\\state=(\\rhostar, \\Svectilde, \\tautilde)$:\n",
    "\\begin{equation}\n",
    "\\ddt{\\state} + \\grad\\cdot\\flux = \\source,\n",
    "\\end{equation}\n",
    "where $\\rhostar = \\lapse\\rho\\rtgamma u^0$, $\\Svectilde = \\rhostar h \\uvec$, $\\tautilde = \\lapse^2\\rtgamma \\T00 - \\rhostar$. The associated set of primitive variables are $(\\rho, \\vel, \\epsilon)$, which are the rest mass density, fluid 3-velocity, and internal energy (measured in the rest frame).  \n",
    "\n",
    "The flux, $\\flux$ is given by\n",
    "\\begin{equation}\n",
    " \\flux=(\\rhostar \\vel, \\lapse\\rtgamma\\T{j}{\\beta}g_{\\beta i}, \\lapse^2\\rtgamma\\T0j - \\rhostar\\vel\n",
    "\\end{equation}\n",
    "where $\\vel$ is the 3-velocity, and $\\source = (0, \\frac 1 2 \\lapse\\rtgamma \\T{\\lapse}{\\beta}g_{\\lapse\\beta,i}, s)$ is the source function, and\n",
    "\\begin{equation}\n",
    "s = \\lapse\\rtgamma\\left[\\left(\\T00\\beta^i\\beta^j + 2\\T0i\\beta^j\\right)K_{ij} - \\left(\\T00\\beta^i + \\T0i\\right)\\partial_i\\lapse\\right]\n",
    "\\end{equation}\n",
    "The stress-energy tensor for a perfect fluid is written as \n",
    "\\begin{equation}\n",
    "\\T{\\mu}{\\nu} = \\rho h u^{\\mu} u^{\\nu} + P g^{\\mu\\nu},\n",
    "\\end{equation}\n",
    "where $h = 1 + \\epsilon + P/\\rho$ is the specific enthalpy and $u^{\\mu}$ are the respective components of the four velocity.  \n",
    "\n",
    "Noting that the mass $\\flux$ is defined in terms of $\\rhostar$ and $\\vel$, we need to first find a mapping between $\\vel$ and $u$.  \n",
    "\n",
    "### Alternative formulation\n",
    "\n",
    "The Athena++ folks have an alternative formulation that might be superior.  \n",
    "Begin with the continuity equation\n",
    "\\begin{equation}\n",
    "\\grad_{\\mu}\\rho u^{\\mu} = 0,\n",
    "\\end{equation}\n",
    "where $\\grad$ is the covariant derivative.  This can be mapped directly to \n",
    "\\begin{equation}\n",
    "\\partial_{0} \\sqrt{-g}\\rho u^0 + \\partial_i\\sqrt{-g} \\rho u^0 v^i = 0 \n",
    "\\end{equation}\n",
    "which we can identify with $\\rhostar = \\alpha\\rtgamma \\rho u^0$ because $\\sqrt{-g} = \\alpha\\rtgamma$.\n",
    "\n",
    "Now the second equation is the conservation of energy-momentum which we write as\n",
    "\\begin{equation}\n",
    "\\grad_{\\nu}T^{\\nu}_{\\mu} = 0 \n",
    "\\end{equation}\n",
    "writing this out we have \n",
    "\\begin{equation}\n",
    "\\partial_0 g_{\\mu\\alpha}T^{\\alpha 0} + \\partial_i g_{\\mu\\alpha}T^{\\alpha i} - \\Gamma_{\\mu\\alpha}^{\\gamma} g_{\\gamma\\beta}T^{\\alpha\\beta} = 0 \n",
    "\\end{equation}\n",
    "Noting that\n",
    "\\begin{equation}\n",
    "\\Gamma^{\\alpha}_{\\beta\\gamma} = \\frac 1 2 g^{\\alpha\\delta}\\left(\\partial_{\\gamma}g_{\\beta\\delta} + \\partial_{\\beta}g_{\\gamma\\delta} - \\partial_{\\delta}g_{\\beta\\gamma}\\right)\n",
    "\\end{equation}\n",
    "Writing this all out, we note the last term is\n",
    "\\begin{equation}\n",
    "\\Gamma_{\\mu\\alpha}^{\\gamma} g_{\\gamma\\beta}T^{\\alpha\\beta} =\n",
    "\\frac 1 2 g^{\\gamma\\delta}\\left(\\partial_{\\alpha}g_{\\mu\\delta} + \\partial_{\\mu}g_{\\alpha \\delta} - \\partial_{\\delta}g_{\\mu \\alpha}\\right) T_{\\gamma}^{\\alpha} = \n",
    "\\frac 1 2 \\left(\\partial_{\\alpha}g_{\\mu\\delta} + \\partial_{\\mu}g_{\\alpha \\delta} - \\partial_{\\delta}g_{\\mu \\alpha}\\right)\n",
    "T^{\\alpha\\delta}\n",
    "\\end{equation}\n",
    "We sum over $\\alpha$ and $\\delta$, but note that we are antisymmetric in first and last terms in $\\alpha$ and $\\delta$ in the () but symmetric in $T_{\\alpha\\delta}$ so we have\n",
    "\\begin{equation}\n",
    "\\Gamma_{\\mu\\alpha}^{\\gamma} g_{\\gamma\\beta}T^{\\alpha\\beta} = \\frac 1 2 \\partial_{\\mu}g_{\\alpha \\delta} T^{\\alpha\\delta}\n",
    "\\end{equation}\n",
    "\n",
    "Thus we have \n",
    "\\begin{equation}\n",
    "\\partial_0 T^{0}_{\\mu} + \\partial_i T^{i}_{\\mu} = \\frac 1 2 \\partial_{\\mu}g_{\\alpha \\delta} T^{\\alpha\\delta}\n",
    "\\end{equation}\n",
    "The $\\mu = (1,2,3)$, we almost get back the equations in the standard formulation\n",
    "\\begin{equation}\n",
    "\\partial_0 \\rho h u^0 u_i + \\partial_j T^j_i = \\frac 1 2 \\partial_{i}g_{\\alpha \\delta} T^{\\alpha\\delta},\n",
    "\\end{equation}\n",
    "which modulo factors of $\\lapse\\rtgamma$ in front is the same as the \"standard\" equations.\n",
    "\n",
    "The $T^0_0$ term is more interesting. Here we have\n",
    "\\begin{equation}\n",
    "\\partial_0 (\\rho h u^0 u_0 + + \\partial_j T^j_i = \\frac 1 2 \\partial_{0}g_{\\alpha \\delta} T^{\\alpha\\delta},\n",
    "\\end{equation}\n",
    "\n",
    "However, the disadvantage is that we need the time derivative of the metric."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#mapping): Primitive to Conservative Mapping\n",
    "1. [Step 2](#zach): Compute $u^0$ from the Valencia 3-velocity (Zach step)\n",
    "1. [Step 3](#flux): Compute the flux\n",
    "1. [Step 4](#source): Source Terms\n",
    "1. [Step 5](#rotation): Rotation\n",
    "1. [Step 6](#solver): Conservative to Primitive Solver\n",
    "1. [Step 7](#lorentz): Lorentz Boosts\n",
    "1. [Step 8](#TmunuSph): Compute $T^{\\mu\\nu}$ in Cartesian Coordinates\n",
    "1. [Step 9](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mapping'></a>\n",
    "\n",
    "# Step 1: Primitive to Conservative Mapping\n",
    "$$\\label{mapping}$$\n",
    "\n",
    "We want to make a mapping from the primitives to conserved variables following Zach notebook:\n",
    "\\begin{equation}\n",
    "(\\rho, \\vel, \\epsilon) \\rightarrow (\\rhostar = \\lapse\\rho\\rtgamma u^0, \\Svectilde = \\rhostar h \\uvec, \\tautilde = \\lapse^2\\rtgamma \\T00 - \\rhostar).\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:22.838667Z",
     "iopub.status.busy": "2021-03-07T17:11:22.824256Z",
     "iopub.status.idle": "2021-03-07T17:11:23.933919Z",
     "shell.execute_reply": "2021-03-07T17:11:23.932608Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+prim2Con.h\"\n",
      "{\n",
      "  const double tmp_0 = (1.0/((GAMMA_SPEED_LIMIT)*(GAMMA_SPEED_LIMIT)));\n",
      "  const double tmp_3 = (1.0/((mi.alpha)*(mi.alpha)));\n",
      "  const double tmp_4 = mi.betaX + vx;\n",
      "  const double tmp_5 = mi.gamDDxx*tmp_3*((tmp_4)*(tmp_4));\n",
      "  const double tmp_7 = mi.betaY + vy;\n",
      "  const double tmp_8 = mi.gamDDyy*tmp_3*((tmp_7)*(tmp_7));\n",
      "  const double tmp_10 = mi.betaZ + vz;\n",
      "  const double tmp_11 = mi.gamDDzz*((tmp_10)*(tmp_10))*tmp_3;\n",
      "  const double tmp_14 = mi.gamDDxy*tmp_3*tmp_4*tmp_7;\n",
      "  const double tmp_15 = mi.gamDDxz*tmp_10*tmp_3*tmp_4;\n",
      "  const double tmp_16 = mi.gamDDyz*tmp_10*tmp_3*tmp_7;\n",
      "  const double tmp_20 = (1.0/2.0)*fabs(-tmp_0 - tmp_11 - 2*tmp_14 - 2*tmp_15 - 2*tmp_16 - tmp_5 - tmp_8 + 1);\n",
      "  const double tmp_21 = (1.0/2.0)*tmp_0 - 1.0/2.0*tmp_11 - tmp_14 - tmp_15 - tmp_16 + tmp_20 - 1.0/2.0*tmp_5 - 1.0/2.0*tmp_8 + 1.0/2.0;\n",
      "  const double tmp_22 = (1.0/sqrt(tmp_21));\n",
      "  const double tmp_23 = sqrt((-1.0/2.0*tmp_0 + (1.0/2.0)*tmp_11 + tmp_14 + tmp_15 + tmp_16 - tmp_20 + (1.0/2.0)*tmp_5 + (1.0/2.0)*tmp_8 + 1.0/2.0)/(TINYDOUBLE + tmp_11 + 2*tmp_14 + 2*tmp_15 + 2*tmp_16 + tmp_5 + tmp_8));\n",
      "  const double tmp_24 = -mi.betaX + tmp_23*tmp_4;\n",
      "  const double tmp_25 = sqrt(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy);\n",
      "  const double tmp_26 = rho*tmp_22*tmp_25;\n",
      "  const double tmp_27 = p*tmp_3;\n",
      "  const double tmp_28 = rho*tmp_3*(ie + p/rho + 1)/tmp_21;\n",
      "  const double tmp_29 = -tmp_27 + tmp_28;\n",
      "  const double tmp_30 = mi.betaX*tmp_27 + tmp_24*tmp_28;\n",
      "  const double tmp_31 = mi.betaY*tmp_27 + tmp_28*(-mi.betaY + tmp_23*tmp_7);\n",
      "  const double tmp_32 = mi.betaZ*tmp_27 + tmp_28*(-mi.betaZ + tmp_10*tmp_23);\n",
      "  const double tmp_33 = mi.alpha*tmp_25;\n",
      "  u4U1 = tmp_22*tmp_24/mi.alpha;\n",
      "  con[iRhoStar] = tmp_26;\n",
      "  con[iSx] = tmp_33*(mi.gamDDxx*tmp_30 + mi.gamDDxy*tmp_31 + mi.gamDDxz*tmp_32 + tmp_29*(mi.betaX*mi.gamDDxx + mi.betaY*mi.gamDDxy + mi.betaZ*mi.gamDDxz));\n",
      "  con[iSy] = tmp_33*(mi.gamDDxy*tmp_30 + mi.gamDDyy*tmp_31 + mi.gamDDyz*tmp_32 + tmp_29*(mi.betaX*mi.gamDDxy + mi.betaY*mi.gamDDyy + mi.betaZ*mi.gamDDyz));\n",
      "  con[iSz] = tmp_33*(mi.gamDDxz*tmp_30 + mi.gamDDyz*tmp_31 + mi.gamDDzz*tmp_32 + tmp_29*(mi.betaX*mi.gamDDxz + mi.betaY*mi.gamDDyz + mi.betaZ*mi.gamDDzz));\n",
      "  con[iTau] = ((mi.alpha)*(mi.alpha))*tmp_25*tmp_29 - tmp_26;\n",
      "}\n",
      "Wrote to file \"NRPY+detg.h\"\n",
      "/*\n",
      " *  Original SymPy expression:\n",
      " *  \"detg = mi.alpha*sqrt(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy)\"\n",
      " */\n",
      "{\n",
      "  detg = mi.alpha*sqrt(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy);\n",
      "}\n",
      "Wrote to file \"NRPY+gamUU.h\"\n",
      "/*\n",
      " *  Original SymPy expressions:\n",
      " *  \"[gamUUxx = (mi.gamDDyy*mi.gamDDzz - mi.gamDDyz**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy),\n",
      " *    gamUUxy = (-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy),\n",
      " *    gamUUxz = (mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy),\n",
      " *    gamUUyy = (mi.gamDDxx*mi.gamDDzz - mi.gamDDxz**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy),\n",
      " *    gamUUyz = (-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy),\n",
      " *    gamUUzz = (mi.gamDDxx*mi.gamDDyy - mi.gamDDxy**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy)]\"\n",
      " */\n",
      "{\n",
      "  const double tmp_5 = (1.0/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy));\n",
      "  gamUUxx = tmp_5*(mi.gamDDyy*mi.gamDDzz - ((mi.gamDDyz)*(mi.gamDDyz)));\n",
      "  gamUUxy = tmp_5*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz);\n",
      "  gamUUxz = tmp_5*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy);\n",
      "  gamUUyy = tmp_5*(mi.gamDDxx*mi.gamDDzz - ((mi.gamDDxz)*(mi.gamDDxz)));\n",
      "  gamUUyz = tmp_5*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz);\n",
      "  gamUUzz = tmp_5*(mi.gamDDxx*mi.gamDDyy - ((mi.gamDDxy)*(mi.gamDDxy)));\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "import GRHD.equations as Ge  # NRPy: Implementation of GRHD equations in Cartesian coordinates\n",
    "import indexedexp as ixp     # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import sympy as sp           # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "from outputC import outputC  # NRPy+: Basic C code output functionality\n",
    "\n",
    "# declare gammaDD\n",
    "gammaDD = ixp.zerorank2()\n",
    "components = [\"xx\", \"xy\", \"xz\", \"yy\", \"yz\", \"zz\"]\n",
    "names = \"\"\n",
    "for comp in components :\n",
    "    names = names + \"mi.gamDD{0} \".format(comp)\n",
    "\n",
    "gxx, gxy, gxz, gyy, gyz, gzz = sp.symbols( names)\n",
    "\n",
    "gammaDD[0][0] = gxx\n",
    "gammaDD[0][1] = gxy\n",
    "gammaDD[0][2] = gxz\n",
    "gammaDD[1][0] = gxy\n",
    "gammaDD[1][1] = gyy\n",
    "gammaDD[1][2] = gyz\n",
    "gammaDD[2][0] = gxz\n",
    "gammaDD[2][1] = gyz\n",
    "gammaDD[2][2] = gzz\n",
    "\n",
    "#declare alpha\n",
    "alpha = sp.symbols( \"mi.alpha\")\n",
    "\n",
    "#declare beta\n",
    "betaU = ixp.zerorank1()\n",
    "for i, comp in enumerate([\"X\", \"Y\", \"Z\"]) :\n",
    "    betaU[i] = sp.symbols( \"mi.beta{0}\".format(comp), real=True)\n",
    "\n",
    "#now get the primitives\n",
    "rho_b, epsilon, P = sp.symbols(\"rho ie p\")\n",
    "\n",
    "#get the 3-velocities\n",
    "vU = ixp.zerorank1()\n",
    "for i, comp in enumerate( [\"vx\", \"vy\", \"vz\"]) :\n",
    "    vU[i] = sp.symbols(\"{0}\".format(comp))\n",
    "\n",
    "Ge.u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit(alpha,betaU,gammaDD, vU)\n",
    "\n",
    "u4U = Ge.u4U_ito_vU\n",
    "# Zach says: Probably want to adopt speed-limited vU[i], Ge.rescaledvU[i], here, a la:\n",
    "# for i in range(3):\n",
    "# ... vU[i] = Ge.rescaledvU[i]\n",
    "\n",
    "# First compute stress-energy tensor T4UU and T4UD:\n",
    "Ge.compute_T4UU(gammaDD,betaU,alpha, rho_b,P,epsilon,u4U)\n",
    "Ge.compute_T4UD(gammaDD,betaU,alpha, Ge.T4UU)\n",
    "\n",
    "# Next sqrt(gamma)\n",
    "Ge.compute_sqrtgammaDET(gammaDD)\n",
    "\n",
    "# Compute conservative variables in terms of primitive variables\n",
    "Ge.compute_rho_star( alpha, Ge.sqrtgammaDET, rho_b, u4U)\n",
    "Ge.compute_tau_tilde(alpha, Ge.sqrtgammaDET, Ge.T4UU,Ge.rho_star)\n",
    "Ge.compute_S_tildeD( alpha, Ge.sqrtgammaDET, Ge.T4UD)\n",
    "\n",
    "# Zach says: Why only output u^x? Debugging reasons?\n",
    "outputC([u4U[1],     Ge.rho_star, Ge.S_tildeD[0], Ge.S_tildeD[1], Ge.S_tildeD[2], Ge.tau_tilde],\n",
    "        [\"u4U1\", \"con[iRhoStar]\",     \"con[iSx]\",     \"con[iSy]\",     \"con[iSz]\",  \"con[iTau]\"],\n",
    "        filename=\"NRPY+prim2Con.h\", params=\"outCverbose=False\")\n",
    "!cat NRPY+prim2Con.h\n",
    "\n",
    "outputC([Ge.sqrtgammaDET*alpha], [\"detg\"], filename=\"NRPY+detg.h\")\n",
    "!cat NRPY+detg.h\n",
    "\n",
    "\n",
    "gammaUU, gammabarDet = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "outputC([gammaUU[0][0],gammaUU[0][1],gammaUU[0][2],gammaUU[1][1],gammaUU[1][2],gammaUU[2][2]],\n",
    "        [    \"gamUUxx\",    \"gamUUxy\",    \"gamUUxz\",    \"gamUUyy\",    \"gamUUyz\",    \"gamUUzz\"],\n",
    "        filename=\"NRPY+gamUU.h\")\n",
    "!cat NRPY+gamUU.h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='flux'></a>\n",
    "\n",
    "# Step 3: Compute the flux\n",
    "$$\\label{flux}$$\n",
    "\n",
    "The fluxes are as follows\n",
    "\\begin{equation}\n",
    "\\frac{\\partial}{\\partial t} \n",
    "\\begin{pmatrix}\n",
    "\\rhostar\\\\\n",
    "\\Svectilde\\\\\n",
    "\\tautilde\n",
    "\\end{pmatrix} + \\frac{\\partial}{\\partial x^j}\\begin{pmatrix} \\rhostar v^j\\\\\n",
    "\\lapse\\rtgamma T^j_i\\\\ \\lapse^2\\rtgamma T^{0j} - \\rhostar v^j\n",
    "\\end{pmatrix}  = \\begin{pmatrix} 0 \\\\ \\frac 1 2 \\lapse\\rtgamma T^{\\alpha\\beta}g_{\\alpha\\beta,i} \\\\ s \\end{pmatrix}\n",
    "\\end{equation}\n",
    "so the flux is \n",
    "\\begin{equation}\n",
    "\\mathcal{F} = \\begin{pmatrix} \\rhostar v^i \\\\ \\lapse\\rtgamma T^i_k \\\\ \\lapse^2\\rtgamma T^{0i} - \\rhostar v^i\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "In the moving-mesh formalism, the flux is just taken along the x directions so we have\n",
    "\\begin{equation}\n",
    "\\mathcal{F} = \\begin{pmatrix} \\rhostar v^1 \\\\ \\lapse\\rtgamma T^1_k \\\\ \\lapse^2\\rtgamma T^{01} - \\rhostar v^1\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "Note that we will need to rotate $T^{\\mu\\nu}$ and $g_{\\mu\\nu}$ to get the right orientation.\n",
    "In order to do this, we must first compute the stress energy tensor:\n",
    "\\begin{equation}\n",
    "T^{\\mu\\nu} = \\rho h u^{\\mu}u^{\\nu} + Pg^{\\mu\\nu} = \\rho h (u^0)^2v^iv^j + P g^{\\mu\\nu}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:23.994875Z",
     "iopub.status.busy": "2021-03-07T17:11:23.958983Z",
     "iopub.status.idle": "2021-03-07T17:11:24.625036Z",
     "shell.execute_reply": "2021-03-07T17:11:24.623498Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+calFlux.h\"\n",
      "{\n",
      "  const double tmp_6 = mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy;\n",
      "  const double tmp_7 = sqrt(tmp_6);\n",
      "  const double tmp_8 = (1.0/((GAMMA_SPEED_LIMIT)*(GAMMA_SPEED_LIMIT)));\n",
      "  const double tmp_11 = (1.0/((mi.alpha)*(mi.alpha)));\n",
      "  const double tmp_12 = mi.betaX + vx;\n",
      "  const double tmp_13 = mi.gamDDxx*tmp_11*((tmp_12)*(tmp_12));\n",
      "  const double tmp_15 = mi.betaY + vy;\n",
      "  const double tmp_16 = mi.gamDDyy*tmp_11*((tmp_15)*(tmp_15));\n",
      "  const double tmp_18 = mi.betaZ + vz;\n",
      "  const double tmp_19 = mi.gamDDzz*tmp_11*((tmp_18)*(tmp_18));\n",
      "  const double tmp_22 = tmp_11*tmp_12*tmp_15;\n",
      "  const double tmp_24 = mi.gamDDxz*tmp_11*tmp_12*tmp_18;\n",
      "  const double tmp_25 = mi.gamDDyz*tmp_11*tmp_15*tmp_18;\n",
      "  const double tmp_26 = 2*mi.gamDDxy*tmp_22;\n",
      "  const double tmp_29 = (1.0/2.0)*fabs(-tmp_13 - tmp_16 - tmp_19 - 2*tmp_24 - 2*tmp_25 - tmp_26 - tmp_8 + 1);\n",
      "  const double tmp_30 = -mi.gamDDxy*tmp_22 - 1.0/2.0*tmp_13 - 1.0/2.0*tmp_16 - 1.0/2.0*tmp_19 - tmp_24 - tmp_25 + tmp_29 + (1.0/2.0)*tmp_8 + 1.0/2.0;\n",
      "  const double tmp_31 = rho*tmp_7/sqrt(tmp_30);\n",
      "  const double tmp_32 = p*tmp_11;\n",
      "  const double tmp_33 = sqrt((mi.gamDDxy*tmp_22 + (1.0/2.0)*tmp_13 + (1.0/2.0)*tmp_16 + (1.0/2.0)*tmp_19 + tmp_24 + tmp_25 - tmp_29 - 1.0/2.0*tmp_8 + 1.0/2.0)/(TINYDOUBLE + tmp_13 + tmp_16 + tmp_19 + 2*tmp_24 + 2*tmp_25 + tmp_26));\n",
      "  const double tmp_34 = -mi.betaX + tmp_12*tmp_33;\n",
      "  const double tmp_35 = rho*tmp_11*(ie + p/rho + 1)/tmp_30;\n",
      "  const double tmp_36 = tmp_34*tmp_35;\n",
      "  const double tmp_37 = mi.betaX*tmp_32 + tmp_36;\n",
      "  const double tmp_41 = faceVelocity[0]*norm[0] + faceVelocity[1]*norm[1] + faceVelocity[2]*norm[2];\n",
      "  const double tmp_42 = mi.betaX*mi.gamDDxx + mi.betaY*mi.gamDDxy + mi.betaZ*mi.gamDDxz;\n",
      "  const double tmp_43 = -tmp_32 + tmp_35;\n",
      "  const double tmp_44 = -mi.betaY + tmp_15*tmp_33;\n",
      "  const double tmp_46 = mi.betaY*tmp_32 + tmp_35*tmp_44;\n",
      "  const double tmp_47 = -mi.betaZ + tmp_18*tmp_33;\n",
      "  const double tmp_48 = mi.betaZ*tmp_32 + tmp_35*tmp_47;\n",
      "  const double tmp_49 = mi.alpha*tmp_7;\n",
      "  const double tmp_50 = tmp_41*tmp_49;\n",
      "  const double tmp_51 = (1.0/(tmp_6));\n",
      "  const double tmp_52 = p*(-((mi.betaX)*(mi.betaX))*tmp_11 + tmp_51*(mi.gamDDyy*mi.gamDDzz - ((mi.gamDDyz)*(mi.gamDDyz)))) + ((tmp_34)*(tmp_34))*tmp_35;\n",
      "  const double tmp_54 = p*(-mi.betaX*mi.betaY*tmp_11 + tmp_51*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz)) + tmp_36*tmp_44;\n",
      "  const double tmp_56 = p*(-mi.betaX*mi.betaZ*tmp_11 + tmp_51*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy)) + tmp_36*tmp_47;\n",
      "  const double tmp_58 = norm[0]*tmp_49;\n",
      "  const double tmp_59 = p*(-((mi.betaY)*(mi.betaY))*tmp_11 + tmp_51*(mi.gamDDxx*mi.gamDDzz - ((mi.gamDDxz)*(mi.gamDDxz)))) + tmp_35*((tmp_44)*(tmp_44));\n",
      "  const double tmp_60 = p*(-mi.betaY*mi.betaZ*tmp_11 + tmp_51*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz)) + tmp_35*tmp_44*tmp_47;\n",
      "  const double tmp_61 = norm[1]*tmp_49;\n",
      "  const double tmp_62 = p*(-((mi.betaZ)*(mi.betaZ))*tmp_11 + tmp_51*(mi.gamDDxx*mi.gamDDyy - ((mi.gamDDxy)*(mi.gamDDxy)))) + tmp_35*((tmp_47)*(tmp_47));\n",
      "  const double tmp_63 = norm[2]*tmp_49;\n",
      "  const double tmp_64 = mi.betaX*mi.gamDDxy + mi.betaY*mi.gamDDyy + mi.betaZ*mi.gamDDyz;\n",
      "  const double tmp_66 = mi.betaX*mi.gamDDxz + mi.betaY*mi.gamDDyz + mi.betaZ*mi.gamDDzz;\n",
      "  const double tmp_67 = ((mi.alpha)*(mi.alpha))*tmp_7;\n",
      "  temp_rho_star = tmp_31;\n",
      "  temp_T4UU01 = tmp_37;\n",
      "  flux[iRhoStar] = norm[0]*tmp_31*vx + norm[1]*tmp_31*vy + norm[2]*tmp_31*vz - tmp_31*tmp_41;\n",
      "  flux[iSx] = -tmp_50*(mi.gamDDxx*tmp_37 + mi.gamDDxy*tmp_46 + mi.gamDDxz*tmp_48 + tmp_42*tmp_43) + tmp_58*(mi.gamDDxx*tmp_52 + mi.gamDDxy*tmp_54 + mi.gamDDxz*tmp_56 + tmp_37*tmp_42) + tmp_61*(mi.gamDDxx*tmp_54 + mi.gamDDxy*tmp_59 + mi.gamDDxz*tmp_60 + tmp_42*tmp_46) + tmp_63*(mi.gamDDxx*tmp_56 + mi.gamDDxy*tmp_60 + mi.gamDDxz*tmp_62 + tmp_42*tmp_48);\n",
      "  flux[iSy] = -tmp_50*(mi.gamDDxy*tmp_37 + mi.gamDDyy*tmp_46 + mi.gamDDyz*tmp_48 + tmp_43*tmp_64) + tmp_58*(mi.gamDDxy*tmp_52 + mi.gamDDyy*tmp_54 + mi.gamDDyz*tmp_56 + tmp_37*tmp_64) + tmp_61*(mi.gamDDxy*tmp_54 + mi.gamDDyy*tmp_59 + mi.gamDDyz*tmp_60 + tmp_46*tmp_64) + tmp_63*(mi.gamDDxy*tmp_56 + mi.gamDDyy*tmp_60 + mi.gamDDyz*tmp_62 + tmp_48*tmp_64);\n",
      "  flux[iSz] = -tmp_50*(mi.gamDDxz*tmp_37 + mi.gamDDyz*tmp_46 + mi.gamDDzz*tmp_48 + tmp_43*tmp_66) + tmp_58*(mi.gamDDxz*tmp_52 + mi.gamDDyz*tmp_54 + mi.gamDDzz*tmp_56 + tmp_37*tmp_66) + tmp_61*(mi.gamDDxz*tmp_54 + mi.gamDDyz*tmp_59 + mi.gamDDzz*tmp_60 + tmp_46*tmp_66) + tmp_63*(mi.gamDDxz*tmp_56 + mi.gamDDyz*tmp_60 + mi.gamDDzz*tmp_62 + tmp_48*tmp_66);\n",
      "  flux[iTau] = norm[0]*(-tmp_31*vx + tmp_37*tmp_67) + norm[1]*(-tmp_31*vy + tmp_46*tmp_67) + norm[2]*(-tmp_31*vz + tmp_48*tmp_67) - tmp_41*(-tmp_31 + tmp_43*tmp_67);\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# Next compute fluxes of conservative variables\n",
    "Ge.compute_rho_star_fluxU(                            vU,        Ge.rho_star)\n",
    "Ge.compute_tau_tilde_fluxU(alpha, Ge.sqrtgammaDET,    vU,Ge.T4UU,Ge.rho_star)\n",
    "Ge.compute_S_tilde_fluxUD( alpha, Ge.sqrtgammaDET,       Ge.T4UD)\n",
    "\n",
    "\n",
    "normD = ixp.zerorank1()\n",
    "normD[0], normD[1], normD[2] = sp.symbols(\"norm[0] norm[1] norm[2]\", real=True)\n",
    "\n",
    "faceVelU = ixp.zerorank1()\n",
    "faceVelU[0], faceVelU[1], faceVelU[2] = sp.symbols(\"faceVelocity[0] faceVelocity[1] faceVelocity[2]\", real=True)\n",
    "# Zach says: don't forget to limit the velocities after they are computed!\n",
    "\n",
    "faceVelNorm = sp.sympify(0)\n",
    "for i in range(3) :\n",
    "    faceVelNorm += normD[i]*faceVelU[i]\n",
    "\n",
    "exprArray = []\n",
    "nameArray = []\n",
    "exprArray.append( Ge.rho_star)\n",
    "nameArray.append( \"temp_rho_star\")\n",
    "exprArray.append( Ge.T4UU[0][1])\n",
    "nameArray.append( \"temp_T4UU01\")\n",
    "\n",
    "rho_star_flux = sp.sympify(0)\n",
    "for i in range(3) :\n",
    "    rho_star_flux += Ge.rho_star_fluxU[i]*normD[i]\n",
    "rho_star_flux -= Ge.rho_star*faceVelNorm\n",
    "\n",
    "exprArray.append( rho_star_flux)\n",
    "nameArray.append( \"flux[iRhoStar]\")\n",
    "\n",
    "\n",
    "tau_tilde_flux = sp.sympify(0)\n",
    "for i in range(3) :\n",
    "    tau_tilde_flux += Ge.tau_tilde_fluxU[i]*normD[i]\n",
    "\n",
    "tau_tilde_flux -= Ge.tau_tilde*faceVelNorm\n",
    "\n",
    "S_tilde_fluxD = ixp.zerorank1()\n",
    "\n",
    "for i in range(3) :\n",
    "    S_tilde_fluxD[i] -= Ge.S_tildeD[i]*faceVelNorm\n",
    "    for j in range(3) :\n",
    "        S_tilde_fluxD[i] += Ge.S_tilde_fluxUD[j][i]*normD[j]\n",
    "\n",
    "for j, comp in enumerate([\"x\",\"y\", \"z\"]) :\n",
    "    exprArray.append( S_tilde_fluxD[j])\n",
    "    nameArray.append( \"flux[iS{0}]\".format(comp))\n",
    "\n",
    "exprArray.append( tau_tilde_flux)\n",
    "nameArray.append( \"flux[iTau]\")\n",
    "\n",
    "#for expr, name in zip( exprArray, nameArray) :\n",
    "#    print( name)\n",
    "outputC(exprArray, nameArray, filename=\"NRPY+calFlux.h\", params=\"outCverbose=False\")\n",
    "\n",
    "!cat NRPY+calFlux.h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='source'></a>\n",
    "\n",
    "# Step 4: Source Terms\n",
    "$$\\label{source}$$\n",
    "\n",
    "The sources terms are for mass, momentum and energy are: \n",
    "\\begin{equation}\n",
    "\\source = (0, \\frac 1 2 \\lapse\\rtgamma \\T{\\alpha}{\\beta}g_{\\alpha\\beta,i}, s),\n",
    "\\end{equation}\n",
    "For a time stationary metric $s\\neq 0$, so we will ignore this until the next section.  As for the rest, we need to define derivatives of the metric.  Suppose I have done this already.  Then the code for the source terms is:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:24.656686Z",
     "iopub.status.busy": "2021-03-07T17:11:24.635149Z",
     "iopub.status.idle": "2021-03-07T17:11:25.132213Z",
     "shell.execute_reply": "2021-03-07T17:11:25.130743Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+calSources.h\"\n",
      "{\n",
      "  const double tmp_0 = (1.0/(h));\n",
      "  const double tmp_1 = (1.0/2.0)*tmp_0;\n",
      "  const double tmp_2 = tmp_1*(-mi_minus[0].alpha + mi_plus[0].alpha);\n",
      "  const double tmp_9 = mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy;\n",
      "  const double tmp_10 = sqrt(tmp_9);\n",
      "  const double tmp_12 = (1.0/((mi.alpha)*(mi.alpha)));\n",
      "  const double tmp_13 = p*tmp_12;\n",
      "  const double tmp_14 = (1.0/((GAMMA_SPEED_LIMIT)*(GAMMA_SPEED_LIMIT)));\n",
      "  const double tmp_16 = mi.betaX + vx;\n",
      "  const double tmp_17 = mi.gamDDxx*tmp_12*((tmp_16)*(tmp_16));\n",
      "  const double tmp_19 = mi.betaY + vy;\n",
      "  const double tmp_20 = mi.gamDDyy*tmp_12*((tmp_19)*(tmp_19));\n",
      "  const double tmp_22 = mi.betaZ + vz;\n",
      "  const double tmp_23 = mi.gamDDzz*tmp_12*((tmp_22)*(tmp_22));\n",
      "  const double tmp_26 = tmp_12*tmp_16*tmp_19;\n",
      "  const double tmp_28 = mi.gamDDxz*tmp_12*tmp_16*tmp_22;\n",
      "  const double tmp_29 = mi.gamDDyz*tmp_12*tmp_19*tmp_22;\n",
      "  const double tmp_30 = 2*mi.gamDDxy*tmp_26;\n",
      "  const double tmp_33 = (1.0/2.0)*fabs(-tmp_14 - tmp_17 - tmp_20 - tmp_23 - 2*tmp_28 - 2*tmp_29 - tmp_30 + 1);\n",
      "  const double tmp_34 = rho*tmp_12*(ie + p/rho + 1)/(-mi.gamDDxy*tmp_26 + (1.0/2.0)*tmp_14 - 1.0/2.0*tmp_17 - 1.0/2.0*tmp_20 - 1.0/2.0*tmp_23 - tmp_28 - tmp_29 + tmp_33 + 1.0/2.0);\n",
      "  const double tmp_35 = ((mi.alpha)*(mi.alpha))*tmp_10*(-tmp_13 + tmp_34);\n",
      "  const double tmp_36 = (1.0/(tmp_9));\n",
      "  const double tmp_37 = sqrt((mi.gamDDxy*tmp_26 - 1.0/2.0*tmp_14 + (1.0/2.0)*tmp_17 + (1.0/2.0)*tmp_20 + (1.0/2.0)*tmp_23 + tmp_28 + tmp_29 - tmp_33 + 1.0/2.0)/(TINYDOUBLE + tmp_17 + tmp_20 + tmp_23 + 2*tmp_28 + 2*tmp_29 + tmp_30));\n",
      "  const double tmp_38 = -mi.betaX + tmp_16*tmp_37;\n",
      "  const double tmp_39 = mi.alpha*tmp_10;\n",
      "  const double tmp_40 = (1.0/4.0)*tmp_0*tmp_39;\n",
      "  const double tmp_41 = tmp_40*(p*(-((mi.betaX)*(mi.betaX))*tmp_12 + tmp_36*(mi.gamDDyy*mi.gamDDzz - ((mi.gamDDyz)*(mi.gamDDyz)))) + tmp_34*((tmp_38)*(tmp_38)));\n",
      "  const double tmp_42 = -mi.betaY + tmp_19*tmp_37;\n",
      "  const double tmp_43 = tmp_40*(p*(-((mi.betaY)*(mi.betaY))*tmp_12 + tmp_36*(mi.gamDDxx*mi.gamDDzz - ((mi.gamDDxz)*(mi.gamDDxz)))) + tmp_34*((tmp_42)*(tmp_42)));\n",
      "  const double tmp_44 = -mi.betaZ + tmp_22*tmp_37;\n",
      "  const double tmp_45 = tmp_40*(p*(-((mi.betaZ)*(mi.betaZ))*tmp_12 + tmp_36*(mi.gamDDxx*mi.gamDDyy - ((mi.gamDDxy)*(mi.gamDDxy)))) + tmp_34*((tmp_44)*(tmp_44)));\n",
      "  const double tmp_47 = tmp_34*tmp_38;\n",
      "  const double tmp_48 = tmp_1*tmp_39;\n",
      "  const double tmp_49 = tmp_48*(p*(-mi.betaX*mi.betaY*tmp_12 + tmp_36*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz)) + tmp_42*tmp_47);\n",
      "  const double tmp_50 = tmp_48*(p*(-mi.betaX*mi.betaZ*tmp_12 + tmp_36*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy)) + tmp_44*tmp_47);\n",
      "  const double tmp_52 = tmp_48*(p*(-mi.betaY*mi.betaZ*tmp_12 + tmp_36*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz)) + tmp_34*tmp_42*tmp_44);\n",
      "  const double tmp_53 = tmp_1*(-mi_minus[1].alpha + mi_plus[1].alpha);\n",
      "  const double tmp_54 = tmp_1*(-mi_minus[2].alpha + mi_plus[2].alpha);\n",
      "  vSource[0] = -tmp_2*tmp_35 + tmp_41*(-mi_minus[0].gamDDxx + mi_plus[0].gamDDxx) + tmp_43*(-mi_minus[0].gamDDyy + mi_plus[0].gamDDyy) + tmp_45*(-mi_minus[0].gamDDzz + mi_plus[0].gamDDzz) + tmp_49*(-mi_minus[0].gamDDxy + mi_plus[0].gamDDxy) + tmp_50*(-mi_minus[0].gamDDxz + mi_plus[0].gamDDxz) + tmp_52*(-mi_minus[0].gamDDyz + mi_plus[0].gamDDyz);\n",
      "  vSource[1] = -tmp_35*tmp_53 + tmp_41*(-mi_minus[1].gamDDxx + mi_plus[1].gamDDxx) + tmp_43*(-mi_minus[1].gamDDyy + mi_plus[1].gamDDyy) + tmp_45*(-mi_minus[1].gamDDzz + mi_plus[1].gamDDzz) + tmp_49*(-mi_minus[1].gamDDxy + mi_plus[1].gamDDxy) + tmp_50*(-mi_minus[1].gamDDxz + mi_plus[1].gamDDxz) + tmp_52*(-mi_minus[1].gamDDyz + mi_plus[1].gamDDyz);\n",
      "  vSource[2] = -tmp_35*tmp_54 + tmp_41*(-mi_minus[2].gamDDxx + mi_plus[2].gamDDxx) + tmp_43*(-mi_minus[2].gamDDyy + mi_plus[2].gamDDyy) + tmp_45*(-mi_minus[2].gamDDzz + mi_plus[2].gamDDzz) + tmp_49*(-mi_minus[2].gamDDxy + mi_plus[2].gamDDxy) + tmp_50*(-mi_minus[2].gamDDxz + mi_plus[2].gamDDxz) + tmp_52*(-mi_minus[2].gamDDyz + mi_plus[2].gamDDyz);\n",
      "  eSource = tmp_39*(tmp_2*(-mi.betaX*tmp_13 - tmp_47) + tmp_53*(-mi.betaY*tmp_13 - tmp_34*tmp_42) + tmp_54*(-mi.betaZ*tmp_13 - tmp_34*tmp_44));\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# FIXME: Assume static spacetime with KDD = betaU = betaU_dD = 0\n",
    "KDD = ixp.zerorank2()\n",
    "betaU = ixp.zerorank1()\n",
    "betaU_dD = ixp.zerorank2()\n",
    "\n",
    "# Set+evaluate derivatives of alpha, performing 2nd-order finite difference\n",
    "alpha_dD = ixp.zerorank1()\n",
    "h = sp.symbols(\"h\")\n",
    "for i in range(3) :\n",
    "    alpha_plus, alpha_minus = sp.symbols(\"mi_plus[{0}].alpha mi_minus[{0}].alpha\".format(i))\n",
    "    alpha_dD[i] = (alpha_plus - alpha_minus)/(2*h)\n",
    "\n",
    "# Set+evaluate derivatives of gamma_{ij}, performing 2nd-order finite difference\n",
    "gammaDD_dD = ixp.zerorank3()\n",
    "components = [\"xx\", \"xy\", \"xz\", \"yy\", \"yz\", \"zz\"]\n",
    "for i in range(3) :\n",
    "    names_plus = \"\"\n",
    "    names_minus = \"\"\n",
    "    for comp in components :\n",
    "        names_plus  = names_plus  +  \"mi_plus[{0}].gamDD{1} \".format(i, comp)\n",
    "        names_minus = names_minus + \"mi_minus[{0}].gamDD{1} \".format(i, comp)\n",
    "\n",
    "    gxx_plus, gxy_plus, gxz_plus, gyy_plus, gyz_plus, gzz_plus = sp.symbols( names_plus)\n",
    "    gxx_minus, gxy_minus, gxz_minus, gyy_minus, gyz_minus, gzz_minus = sp.symbols( names_minus)\n",
    "\n",
    "    gammaDD_dD[0][0][i] = (gxx_plus - gxx_minus)/(2*h)\n",
    "    gammaDD_dD[0][1][i] = (gxy_plus - gxy_minus)/(2*h)\n",
    "    gammaDD_dD[0][2][i] = (gxz_plus - gxz_minus)/(2*h)\n",
    "    gammaDD_dD[1][0][i] = (gxy_plus - gxy_minus)/(2*h)\n",
    "    gammaDD_dD[1][1][i] = (gyy_plus - gyy_minus)/(2*h)\n",
    "    gammaDD_dD[1][2][i] = (gyz_plus - gyz_minus)/(2*h)\n",
    "    gammaDD_dD[2][0][i] = (gxz_plus - gxz_minus)/(2*h)\n",
    "    gammaDD_dD[2][1][i] = (gyz_plus - gyz_minus)/(2*h)\n",
    "    gammaDD_dD[2][2][i] = (gzz_plus - gzz_minus)/(2*h)\n",
    "\n",
    "# Compute g_{mu nu, i} based on ADM quantities & derivatives defined above\n",
    "Ge.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
    "\n",
    "# Compute source terms for tau tilde & S tilde:\n",
    "Ge.compute_s_source_term(KDD,betaU,alpha, Ge.sqrtgammaDET,alpha_dD,  Ge.T4UU)\n",
    "Ge.compute_S_tilde_source_termD(   alpha, Ge.sqrtgammaDET,Ge.g4DD_zerotimederiv_dD, Ge.T4UU)\n",
    "\n",
    "exprArray = []\n",
    "nameArray = []\n",
    "\n",
    "#momentum terms\n",
    "for i in range(3) :\n",
    "    exprArray.append( Ge.S_tilde_source_termD[i])\n",
    "    nameArray.append( \"vSource[{0}]\".format(i))\n",
    "\n",
    "#tau term\n",
    "exprArray.append( Ge.s_source_term)\n",
    "nameArray.append( \"eSource\")\n",
    "outputC( exprArray, nameArray, filename=\"NRPY+calSources.h\", params=\"outCverbose=False\")\n",
    "\n",
    "!cat NRPY+calSources.h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='solver'></a>\n",
    "\n",
    "# Step 6: Conservative to Primitive Solver\n",
    "$$\\label{solver}$$\n",
    "\n",
    "We now discuss reverse mapping from conservative to primitive variables.\n",
    "Given the lapse, shift vector, and $\\rtgamma$, the mapping between primitive and conserved variables is straightforward.  However, the reverse is not as simple.  In GRMHD, the conservative to primitive solver is amplified by the inclusion of the magnetic field, leading to rather sophisticated root-finding strategies.  The failure rates of these algorithms are low (??), but since this algorithm may be executed several times per timestep for every gridpoint, even a low failure can give unacceptable collective failure rates.  However, for purely polytropic equations of state, e.g., $P\\propto\\rho^{\\Gamma_1}$, the conservative to primitive variable solver is greatly simplified.  \n",
    "\n",
    "To construct the conservative-to-primitive variable solver, we restrict ourselves to polytropic equations of states\n",
    "\\begin{equation}\n",
    "P = P_0\\left(\\frac{\\rho}{\\rho_0}\\right)^{\\Gamma_1} \\quad\\textrm{and}\\quad \\epsilon = \\epsilon_0\\left(\\frac{\\rho}{\\rho_0}\\right)^{\\Gamma_1-1},\n",
    "\\end{equation}\n",
    "where $P_0$, $\\rho_0$, and $\\epsilon_0$ are the fiducial pressure, density, and internal energy, and we have used the relation $P = (\\Gamma_1 - 1)\\rho\\epsilon$.  \n",
    "\n",
    "For such a polytropic equation of state, the energy equation is redundant and effectively we are only concerned with the continuity and momentum equations. The conservative variables of concern are $\\rhostar$ and $\\Svectilde$.  Noting that the shift, $\\alpha$, and $\\rtgamma$ are provided by the Einsteins field equation solver, we can write\n",
    "\\begin{equation}\n",
    "u^0 = \\frac{\\rhostar}{\\alpha\\rtgamma\\rho} = u^0(\\rho)  \\quad\\textrm{and}\\quad \\uvec = \\frac{\\Svectilde}{\\alpha\\rtgamma\\rho h} = \\uvec(\\rho).\n",
    "\\end{equation}\n",
    "Noting that the four-velocity $u^2 = g_{\\mu\\nu}u^{\\mu}u^{\\nu} = g^{00}u^0u^0 + 2g^{0i}u^0\\uvec^i + g_{ij}\\uvec^i\\uvec^j = -1$, we have\n",
    "\\begin{equation}\n",
    " 0 = f(\\rho)\\equiv \\alpha^2\\gamma\\rho^2h^2 + \\left(-\\lapse^2 + \\shift\\cdot\\shift\\right)\\rhostar^2h^2 + 2h\\rhostar\\shift\\cdot\\Svectilde + \\Svectilde\\cdot\\Svectilde,\n",
    "\\end{equation}\n",
    "which is an implicit equation of either $\\rho$ or $u^0$, where $h(\\rho = \\rhostar/(\\alpha\\rtgamma u^0)) = 1 + \\gamma_1 \\epsilon$ which can be inverted by standard nonlinear root-finding algorithms, e.g., Newton-Raphson. \n",
    "\n",
    "We put this all together to define a function, $f(\\rho)$, whose root is zero that we will find via Newton-Raphson.  \n",
    "\n",
    "Several checks must be performed:\n",
    "\n",
    "1. $\\rhostar > 0$ : This check is performed at the very beginning\n",
    "\n",
    "2. $\\rho > \\rho_{\\rm min}$ : This check is performed after the fact\n",
    "\n",
    "3. $u_0 < \\alpha^{-1}\\Gamma_{\\rm max}$ : This check is performed after the fact as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:25.208787Z",
     "iopub.status.busy": "2021-03-07T17:11:25.173004Z",
     "iopub.status.idle": "2021-03-07T17:11:25.628665Z",
     "shell.execute_reply": "2021-03-07T17:11:25.626995Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+rhoRoot.h\"\n",
      "Wrote to file \"NRPY+Stilde2.h\"\n",
      "/*\n",
      " *  Original SymPy expression:\n",
      " *  \"rootRho = con[iSx]**2*(mi.gamDDyy*mi.gamDDzz - mi.gamDDyz**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + 2*con[iSx]*con[iSy]*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + 2*con[iSx]*con[iSz]*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + con[iSy]**2*(mi.gamDDxx*mi.gamDDzz - mi.gamDDxz**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + 2*con[iSy]*con[iSz]*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + con[iSz]**2*(mi.gamDDxx*mi.gamDDyy - mi.gamDDxy**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + rhostar**2*(gamma*p_0*(rho/rho_0)**(gamma - 1)/(rho_0*(gamma - 1)) + 1)**2 + rhostar**2*(2.0*gamma*p_0*(rho/rho_0)**(gamma - 1)/(rho_0*(gamma - 1)) + 2.0)*(con[iSx]*mi.betaX + con[iSy]*mi.betaY + con[iSz]*mi.betaZ)/(mi.alpha*rho*sqrt(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy)) + rhostar**4.0*(gamma*p_0*(rho/rho_0)**(gamma - 1)/(rho_0*(gamma - 1)) + 1)**2.0*(-mi.alpha**2 + mi.betaX**2*mi.gamDDxx + 2*mi.betaX*mi.betaY*mi.gamDDxy + 2*mi.betaX*mi.betaZ*mi.gamDDxz + mi.betaY**2*mi.gamDDyy + 2*mi.betaY*mi.betaZ*mi.gamDDyz + mi.betaZ**2*mi.gamDDzz)/(mi.alpha**2*rho**2*(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy))\"\n",
      " */\n",
      "{\n",
      "  const double tmp_1 = (1.0/(rho_0));\n",
      "  const double tmp_3 = gamma*p_0*tmp_1*pow(rho*tmp_1, gamma - 1)/(gamma - 1);\n",
      "  const double tmp_11 = mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy;\n",
      "  const double tmp_12 = (1.0/(tmp_11));\n",
      "  const double tmp_13 = 2*con[iSx]*tmp_12;\n",
      "  rootRho = ((con[iSx])*(con[iSx]))*tmp_12*(mi.gamDDyy*mi.gamDDzz - ((mi.gamDDyz)*(mi.gamDDyz))) + ((con[iSy])*(con[iSy]))*tmp_12*(mi.gamDDxx*mi.gamDDzz - ((mi.gamDDxz)*(mi.gamDDxz))) + 2*con[iSy]*con[iSz]*tmp_12*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz) + con[iSy]*tmp_13*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz) + ((con[iSz])*(con[iSz]))*tmp_12*(mi.gamDDxx*mi.gamDDyy - ((mi.gamDDxy)*(mi.gamDDxy))) + con[iSz]*tmp_13*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy) + ((rhostar)*(rhostar))*((tmp_3 + 1)*(tmp_3 + 1)) + ((rhostar)*(rhostar))*(2.0*tmp_3 + 2.0)*(con[iSx]*mi.betaX + con[iSy]*mi.betaY + con[iSz]*mi.betaZ)/(mi.alpha*rho*sqrt(tmp_11)) + ((rhostar)*(rhostar)*(rhostar)*(rhostar))*tmp_12*((tmp_3 + 1)*(tmp_3 + 1))*(-((mi.alpha)*(mi.alpha)) + ((mi.betaX)*(mi.betaX))*mi.gamDDxx + 2*mi.betaX*mi.betaY*mi.gamDDxy + 2*mi.betaX*mi.betaZ*mi.gamDDxz + ((mi.betaY)*(mi.betaY))*mi.gamDDyy + 2*mi.betaY*mi.betaZ*mi.gamDDyz + ((mi.betaZ)*(mi.betaZ))*mi.gamDDzz)/(((mi.alpha)*(mi.alpha))*((rho)*(rho)));\n",
      "}\n",
      "/*\n",
      " *  Original SymPy expression:\n",
      " *  \"Stilde2 = con[iSx]**2*(mi.gamDDyy*mi.gamDDzz - mi.gamDDyz**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + 2*con[iSx]*con[iSy]*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + 2*con[iSx]*con[iSz]*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + con[iSy]**2*(mi.gamDDxx*mi.gamDDzz - mi.gamDDxz**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + 2*con[iSy]*con[iSz]*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy) + con[iSz]**2*(mi.gamDDxx*mi.gamDDyy - mi.gamDDxy**2)/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*mi.gamDDyz**2 - mi.gamDDxy**2*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - mi.gamDDxz**2*mi.gamDDyy)\"\n",
      " */\n",
      "{\n",
      "  const double tmp_5 = (1.0/(mi.gamDDxx*mi.gamDDyy*mi.gamDDzz - mi.gamDDxx*((mi.gamDDyz)*(mi.gamDDyz)) - ((mi.gamDDxy)*(mi.gamDDxy))*mi.gamDDzz + 2*mi.gamDDxy*mi.gamDDxz*mi.gamDDyz - ((mi.gamDDxz)*(mi.gamDDxz))*mi.gamDDyy));\n",
      "  const double tmp_6 = 2*con[iSx]*tmp_5;\n",
      "  Stilde2 = ((con[iSx])*(con[iSx]))*tmp_5*(mi.gamDDyy*mi.gamDDzz - ((mi.gamDDyz)*(mi.gamDDyz))) + ((con[iSy])*(con[iSy]))*tmp_5*(mi.gamDDxx*mi.gamDDzz - ((mi.gamDDxz)*(mi.gamDDxz))) + 2*con[iSy]*con[iSz]*tmp_5*(-mi.gamDDxx*mi.gamDDyz + mi.gamDDxy*mi.gamDDxz) + con[iSy]*tmp_6*(-mi.gamDDxy*mi.gamDDzz + mi.gamDDxz*mi.gamDDyz) + ((con[iSz])*(con[iSz]))*tmp_5*(mi.gamDDxx*mi.gamDDyy - ((mi.gamDDxy)*(mi.gamDDxy))) + con[iSz]*tmp_6*(mi.gamDDxy*mi.gamDDyz - mi.gamDDxz*mi.gamDDyy);\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "DIM = 3\n",
    "# Declare rank-1 contravariant (\"v\") vector\n",
    "vU = ixp.declarerank1(\"vU\")\n",
    "shiftU = ixp.zerorank1()\n",
    "rho, gamma1 = sp.symbols(\"rho gamma\")\n",
    "Sx, Sy, Sz = sp.symbols(\"con[iSx] con[iSy] con[iSz]\")\n",
    "p0, rho0, rhostar = sp.symbols(\"p_0 rho_0 rhostar\")\n",
    "# Declare rank-2 covariant gmunu\n",
    "#gammaDD = ixp.declarerank2(\"gammaDD\",\"sym01\")\n",
    "StildeD = ixp.declarerank1(\"StildeD\")\n",
    "lapse, beta_x, beta_y, beta_z = sp.symbols( \"mi.alpha mi.betaX mi.betaY mi.betaZ\")\n",
    "rtgamma = Ge.sqrtgammaDET\n",
    "\n",
    "shiftU[0] = beta_x\n",
    "shiftU[1] = beta_y\n",
    "shiftU[2] = beta_z\n",
    "\n",
    "StildeD[0] = Sx\n",
    "StildeD[1] = Sy\n",
    "StildeD[2] = Sz\n",
    "\n",
    "# gamma = rtgamma*rtgamma <- unused\n",
    "lapse2 = lapse*lapse\n",
    "uU0 = rhostar/(lapse*rtgamma*rho)\n",
    "epsilon = p0/rho0*(rho/rho0)**(gamma1 - 1)/(gamma1 - 1)\n",
    "h = 1 + gamma1*epsilon\n",
    "\n",
    "beta2 = 0.\n",
    "\n",
    "for i in range(DIM) :\n",
    "    for j in range(DIM) :\n",
    "        beta2 += gammaDD[i][j] * shiftU[i]*shiftU[j]\n",
    "\n",
    "betaDotStilde = 0\n",
    "for i in range(DIM) :\n",
    "    betaDotStilde += shiftU[i]*StildeD[i]\n",
    "\n",
    "# Note that this is |Stilde|^2, where the absolute value denotes\n",
    "#    that this is not a proper contraction of Stilde_i, as\n",
    "#    Stilde^i is NOT equal to gamma^{ij} Stilde_j (to understand\n",
    "#    why this is, notice that Stilde_i is proportional to the\n",
    "#    *4D* stress-energy tensor.)\n",
    "Stilde2 = 0\n",
    "for i in range(DIM) :\n",
    "    for j in range(DIM) :\n",
    "        Stilde2 += gammaUU[i][j] * StildeD[i]*StildeD[j]\n",
    "\n",
    "f = rhostar**2*h**2 + (-lapse2 + beta2)*rhostar**2.*h**2.*uU0**2 + 2.*h*rhostar*betaDotStilde*uU0 + Stilde2\n",
    "\n",
    "outputC(f,\"rootRho\",filename=\"NRPY+rhoRoot.h\")\n",
    "outputC(Stilde2, \"Stilde2\", filename=\"NRPY+Stilde2.h\")\n",
    "!cat NRPY+rhoRoot.h\n",
    "!cat NRPY+Stilde2.h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The root solve above finds $\\rho$, which then allows us to get \n",
    "\\begin{equation}\n",
    "u^0 = \\frac{\\rhostar}{\\alpha\\rtgamma\\rho}\\quad\\textrm{and}\\quad \\vel = \\frac{\\uvec}{u^0} = \\frac{\\Svectilde}{\\rhostar h(\\rho)}.\n",
    "\\end{equation}\n",
    "and thus we can find the rest of the primitives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:25.704569Z",
     "iopub.status.busy": "2021-03-07T17:11:25.668942Z",
     "iopub.status.idle": "2021-03-07T17:11:25.721658Z",
     "shell.execute_reply": "2021-03-07T17:11:25.720998Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+getv.h\"\n"
     ]
    }
   ],
   "source": [
    "#rhostar = sp.symbols(\"rhostar\")\n",
    "#StildeU = ixp.declarerank1(\"StildeU\")\n",
    "velU = ixp.zerorank1()\n",
    "\n",
    "#lapse, rtgamma, rho, gamma1, c = sp.symbols(\"lapse rtgamma rho gamma1 c\")\n",
    "rho, rhostar = sp.symbols(\"testPrim[iRho] con[iRhoStar]\")\n",
    "\n",
    "u0 = rhostar/(lapse*rtgamma*rho)\n",
    "epsilon = p0/rho0*(rho/rho0)**(gamma1 - 1)/(gamma1 - 1)\n",
    "\n",
    "h = 1. + gamma1*epsilon\n",
    "\n",
    "for i in range(DIM) :\n",
    "    for j in range(DIM) :\n",
    "        velU[i] += gammaUU[i][j]*StildeD[j]/(rhostar * h)/u0\n",
    "\n",
    "outputC([h,u0,velU[0],velU[1],velU[2]], [\"h\", \"u0\",\"testPrim[ivx]\", \"testPrim[ivy]\", \"testPrim[ivz]\"],filename=\"NRPY+getv.h\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "<a id='lorentz'></a>\n",
    "\n",
    "# Step 7: Lorentz Boosts\n",
    "$$\\label{lorentz}$$\n",
    "\n",
    "We need to boost to the frame of the moving face.  The boost is\n",
    "\\begin{equation}\n",
    "B(\\beta) =\\begin{pmatrix}\n",
    "\\gamma & -\\beta\\gamma n_x & -\\beta\\gamma n_y & -\\beta\\gamma n_z  \\\\\n",
    "-\\beta\\gamma n_x  & 1 + (\\gamma-1)n_x^2 & (\\gamma-1)n_x n_y & (\\gamma-1)n_x n_z\\\\\n",
    "-\\beta\\gamma n_x  & (\\gamma-1)n_y n_x & 1 + (\\gamma-1)n_y^2 & (\\gamma-1)n_y n_z\\\\\n",
    "-\\beta\\gamma n_x & (\\gamma-1) n_z n_x & (\\gamma-1)n_z n_x & 1 + (\\gamma-1)n_z^2 \n",
    "\\end{pmatrix} \n",
    "\\end{equation}\n",
    "And the boost is $X' = B(\\beta) X$, where $X'$ and $X$ are four vectors.\n",
    "\n",
    "So the rest of this is straightforward.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='TmunuSph'></a>\n",
    "\n",
    "# Step 8: Compute $T^{\\mu\\nu}$ in Cartesian Coordinates\n",
    "$$\\label{TmunuSph}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:25.770971Z",
     "iopub.status.busy": "2021-03-07T17:11:25.735052Z",
     "iopub.status.idle": "2021-03-07T17:11:25.931834Z",
     "shell.execute_reply": "2021-03-07T17:11:25.931078Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+getT4UU.h\"\n"
     ]
    }
   ],
   "source": [
    "# declare gammaDD\n",
    "gammaDD = ixp.zerorank2()\n",
    "components = [\"xx\", \"xy\", \"xz\", \"yy\", \"yz\", \"zz\"]\n",
    "names = \"\"\n",
    "for comp in components :\n",
    "    names = names + \"mi.gamDD{0} \".format(comp)\n",
    "\n",
    "g11, g12, g13, g22, g23, g33 = sp.symbols( names)\n",
    "\n",
    "gammaDD[0][0] = g11\n",
    "gammaDD[0][1] = g12\n",
    "gammaDD[0][2] = g13\n",
    "gammaDD[1][0] = g12\n",
    "gammaDD[1][1] = g22\n",
    "gammaDD[1][2] = g23\n",
    "gammaDD[2][0] = g13\n",
    "gammaDD[2][1] = g23\n",
    "gammaDD[2][2] = g33\n",
    "\n",
    "#declare alpha\n",
    "alpha = sp.symbols( \"mi.alpha\")\n",
    "\n",
    "#declare beta\n",
    "betaU = ixp.zerorank1()\n",
    "for i, comp in enumerate([\"X\", \"Y\", \"Z\"]) :\n",
    "    betaU[i] = 0 #NEED A BETTER WAY sp.symbols( \"mi.beta{0}\".format(comp), real=True)\n",
    "\n",
    "#now get the primitives\n",
    "rho_b, epsilon, P = sp.symbols(\"rho ie press\")\n",
    "\n",
    "#get the 3-velocities\n",
    "vU = ixp.zerorank1()\n",
    "for i, comp in enumerate( [\"vx\", \"vy\", \"vz\"]) :\n",
    "    vU[i] = sp.symbols(\"{0}\".format(comp))\n",
    "\n",
    "Ge.u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit(alpha,betaU,gammaDD, vU)\n",
    "\n",
    "u4U = Ge.u4U_ito_vU\n",
    "\n",
    "# First compute stress-energy tensor T4UU and T4UD in Spherical Coordinates:\n",
    "Ge.compute_T4UU(gammaDD,betaU,alpha, rho_b,P,epsilon,u4U)\n",
    "Ge.compute_T4UD(gammaDD,betaU,alpha, Ge.T4UU)\n",
    "\n",
    "outputC([Ge.T4UU[0][0],Ge.T4UU[1][1],Ge.T4UU[2][2],Ge.T4UU[3][3] ], [\"T4UU_diag[0]\", \"T4UU_diag[1]\",\"T4UU_diag[2]\", \"T4UU_diag[3]\"],filename=\"NRPY+getT4UU.h\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 9: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "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\n",
    "[Tutorial-GRHD_Equations-Cartesian-c-code.pdf](Tutorial-GRHD_Equations-Cartesian-c-code.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:25.937823Z",
     "iopub.status.busy": "2021-03-07T17:11:25.936725Z",
     "iopub.status.idle": "2021-03-07T17:11:29.427877Z",
     "shell.execute_reply": "2021-03-07T17:11:29.426808Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-GRHD_Equations-Cartesian-c-code.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-GRHD_Equations-Cartesian-c-code.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GRHD_Equations-Cartesian-c-code\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  },
  "nteract": {
   "version": "nteract-on-jupyter@1.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}