in_progress/Tutorial-GiRaFFE_NRPy-Characteristic_Speeds.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",
    "# The Characteristic GRFFE speeds\n",
    "## Author: Patrick Nelson\n",
    "\n",
    "This notebook documents the function from the original `GiRaFFE` that calculates the speeds of the characteristic GRFFE hydrodynamics waves.\n",
    "\n",
    "**Notebook Status:** <font color=green><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This code has been validated to round-off level agreement with the corresponding code in the original `GiRaFFE` by means of its implementation in [Tutorial-GiRaFFE_NRPy-Stilde-flux](Tutorial-GiRaFFE_NRPy-Stilde-flux.ipynb) and [Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten](Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten.ipynb)\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [GiRaFFE_NRPy/GiRaFFE_NRPy_Characteristic_Speeds.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Characteristic_Speeds.py)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "Here, we will find the speeds at which the hydrodynamics waves propagate. We start from the speed of light (since FFE deals with very diffuse plasmas), which is $c=1.0$ in our chosen units. We then find the speeds $c_+$ and $c_-$ on each face with the function `find_cp_cm`; then, we find minimum and maximum speeds possible from among those. A complete derivation of these speeds is included [below](#derive_speed) as an appendix; here, for brevity, we simply start with the implementation of the code in the original `GiRaFFE`. That code itself is adapted from the more general `IllinoisGRMHD` code; we will implement several more simplifications to help improve performance.\n",
    "\n",
    "Below is the source code for `find_cp_cm`, edited to work with the NRPy+ version of `GiRaFFE`. One edit we need to make in particular is to the term `psim4*gupii` in the definition of `c`; that was written assuming the use of the conformal metric $\\tilde{g}^{ii}$. Since we are not using that here, and are instead using the ADM metric, we should not multiply by $\\psi^{-4}$.\n",
    "\n",
    "```c\n",
    "static inline void find_cp_cm(REAL &cplus,REAL &cminus,const REAL v02,const REAL u0,\n",
    "                              const REAL vi,const REAL lapse,const REAL shifti,\n",
    "                              const REAL gammadet,const REAL gupii) {\n",
    "  const REAL u0_SQUARED=u0*u0;\n",
    "  const REAL ONE_OVER_LAPSE_SQUARED = 1.0/(lapse*lapse);\n",
    "  // sqrtgamma = psi6 -> psim4 = gammadet^(-1.0/3.0)\n",
    "  const REAL psim4 = pow(gammadet,-1.0/3.0);\n",
    "  //Find cplus, cminus:\n",
    "  const REAL a = u0_SQUARED * (1.0-v02) + v02*ONE_OVER_LAPSE_SQUARED;\n",
    "  const REAL b = 2.0* ( shifti*ONE_OVER_LAPSE_SQUARED * v02 - u0_SQUARED * vi * (1.0-v02) );\n",
    "  const REAL c = u0_SQUARED*vi*vi * (1.0-v02) - v02 * ( gupii -\n",
    "                                                               shifti*shifti*ONE_OVER_LAPSE_SQUARED);\n",
    "  REAL detm = b*b - 4.0*a*c;\n",
    "  //ORIGINAL LINE OF CODE:\n",
    "  //if(detm < 0.0) detm = 0.0;\n",
    "  //New line of code (without the if() statement) has the same effect:\n",
    "  detm = sqrt(0.5*(detm + fabs(detm))); /* Based on very nice suggestion from Roland Haas */\n",
    "  \n",
    "  cplus = 0.5*(detm-b)/a;\n",
    "  cminus = -0.5*(detm+b)/a;\n",
    "  if (cplus < cminus) {\n",
    "    const REAL cp = cminus;\n",
    "    cminus = cplus;\n",
    "    cplus = cp;\n",
    "  }\n",
    "}\n",
    "```\n",
    "Comments documenting this have been excised for brevity, but are reproduced in $\\LaTeX$ [below](#derive_speed).\n",
    "\n",
    "<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](#prelim): Preliminaries\n",
    "1. [Step 2](#simplifications): Apply GRFFE assumptions to simplify expressions\n",
    "1. [Step 3](#char_speeds): Compute the characteristic speeds\n",
    "    1. [Step 3.a](#cp_cm): Compute $c_+$ and $c_-$\n",
    "    1. [Step 3.b](#cmax_cmin): Compute $c_\\max$ and $c_\\min$\n",
    "1. [Step 4](#code_validation): Code Validation against `GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds` NRPy+ Module\n",
    "1. [Step 5](#derive_speed): Complete Derivation of the Wave Speeds\n",
    "1. [Step 6](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='prelim'></a>\n",
    "\n",
    "# Step 1: Preliminaries \\[Back to [top](#toc)\\]\n",
    "$$\\label{prelim}$$\n",
    "\n",
    "This first block of code imports the core NRPy+ functionality after first adding the main NRPy+ directory to the path. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "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 sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='simplifications'></a>\n",
    "\n",
    "# Step 2: Apply GRFFE assumptions to simplify expressions \\[Back to [top](#toc)\\]\n",
    "$$\\label{simplifications}$$\n",
    "\n",
    "We could use this code directly, but there's substantial improvement we can make by changing the code into a NRPyfied form. Note the `if` statement; NRPy+ does not know how to handle these, so we must eliminate it if we want to leverage NRPy+'s full power. (Calls to `fabs()` are also cheaper than `if` statements.) This can be done if we rewrite this, taking inspiration from the other eliminated `if` statement documented in the above code block:\n",
    "```c\n",
    "  cp = 0.5*(detm-b)/a;\n",
    "  cm = -0.5*(detm+b)/a;\n",
    "  cplus  = 0.5*(cp+cm+fabs(cp-cm));\n",
    "  cminus = 0.5*(cp+cm-fabs(cp-cm));\n",
    "```\n",
    "This can be simplified further, by substituting `cp` and `cm` into the below equations and eliminating terms as appropriate. First note that `cp+cm = -b/a` and that `cp-cm = detm/a`. Thus,\n",
    "```c\n",
    "  cplus  = 0.5*(-b/a + fabs(detm/a));\n",
    "  cminus = 0.5*(-b/a - fabs(detm/a));\n",
    "```\n",
    "This fulfills the original purpose of the `if` statement in the original code because we have guaranteed that $c_+ \\geq c_-$.\n",
    "\n",
    "This leaves us with an expression that can be much more easily NRPyfied. So, we will rewrite the following in NRPy+, making only minimal changes to be proper Python. However, it turns out that we can make this even simpler. In GRFFE, $v_0^2$ is guaranteed to be exactly one. In GRMHD, this speed was calculated as $$v_{0}^{2} = v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right),$$ where the Alfv&eacute;n speed $v_{\\rm A}^{2}$ $$v_{\\rm A}^{2} = \\frac{b^{2}}{\\rho_{b}h + b^{2}}.$$ So, we can see that when the density $\\rho_b$ goes to zero, $v_{0}^{2} = v_{\\rm A}^{2} = 1$. Then \n",
    "\\begin{align}\n",
    "a &= (u^0)^2 (1-v_0^2) + v_0^2/\\alpha^2 \\\\\n",
    "&= 1/\\alpha^2 \\\\\n",
    "b &= 2 \\left(\\beta^i v_0^2 / \\alpha^2 - (u^0)^2 v^i (1-v_0^2)\\right) \\\\\n",
    "&= 2 \\beta^i / \\alpha^2 \\\\\n",
    "c &= (u^0)^2 (v^i)^2 (1-v_0^2) - v_0^2 \\left(\\gamma^{ii} - (\\beta^i)^2/\\alpha^2\\right) \\\\\n",
    "&= -\\gamma^{ii} + (\\beta^i)^2/\\alpha^2,\n",
    "\\end{align}\n",
    "are simplifications that should save us some time; we can see that $a \\geq 0$ is guaranteed. Note that we also force `detm` to be positive. Thus, `detm/a` is guaranteed to be positive itself, rendering the calls to `nrpyAbs()` superfluous. Furthermore, we eliminate any dependence on the Valencia 3-velocity and the time compoenent of the four-velocity, $u^0$. This leaves us free to solve the quadratic in the familiar way: $$c_\\pm = \\frac{-b \\pm \\sqrt{b^2-4ac}}{2a}$$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='char_speeds'></a>\n",
    "\n",
    "# Step 3: Compute the characteristic speeds \\[Back to [top](#toc)\\]\n",
    "$$\\label{char_speeds}$$\n",
    "\n",
    "<a id='cp_cm'></a>\n",
    "\n",
    "## Step 3.a: Compute $c_+$ and $c_-$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{cp_cm}$$\n",
    "\n",
    "We now have the expressions we need for $c_+$ and $c_-$:\n",
    "$$c_\\pm = \\frac{-b \\pm \\sqrt{b^2-4ac}}{2a},$$\n",
    "where\n",
    "\\begin{align}\n",
    "a &= 1/\\alpha^2 \\\\\n",
    "b &= 2 \\beta^i / \\alpha^2 \\\\\n",
    "c &= -\\gamma^{ii} + (\\beta^i)^2/\\alpha^2.\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll write this as a function so that we can calculate the expressions on-demand for any choice of i\n",
    "def find_cp_cm(lapse,shifti,gammaUUii):\n",
    "    # Inputs:  u0,vi,lapse,shift,gammadet,gupii\n",
    "    # Outputs: cplus,cminus\n",
    "\n",
    "    # a = 1/(alpha^2)\n",
    "    a = sp.sympify(1)/(lapse*lapse)\n",
    "    # b = 2 beta^i / alpha^2\n",
    "    b = sp.sympify(2) * shifti /(lapse*lapse)\n",
    "    # c = -g^{ii} + (beta^i)^2 / alpha^2\n",
    "    c = - gammaUUii + shifti*shifti/(lapse*lapse)\n",
    "\n",
    "    # Now, we are free to solve the quadratic equation as usual. We take care to avoid passing a\n",
    "    # negative value to the sqrt function.\n",
    "    detm = b*b - sp.sympify(4)*a*c\n",
    "\n",
    "    import Min_Max_and_Piecewise_Expressions as noif\n",
    "    detm = sp.sqrt(noif.max_noif(sp.sympify(0),detm))\n",
    "    global cplus,cminus\n",
    "    cplus  = sp.Rational(1,2)*(-b/a + detm/a)\n",
    "    cminus = sp.Rational(1,2)*(-b/a - detm/a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cmax_cmin'></a>\n",
    "\n",
    "## Step 3.b: Compute $c_\\max$ and $c_\\min$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{cmax_cmin}$$\n",
    "\n",
    "In flat spacetime, where $\\alpha=1$, $\\beta^i=0$, and $\\gamma^{ij} = \\delta^{ij}$, $c_+ > 0$ and $c_- < 0$. For the HLLE solver, we will need both `cmax` and `cmin` to be positive; we also want to choose the speed that is larger in magnitude because overestimating the characteristic speeds will help damp unwanted oscillations. (However, in GRFFE, we only get one $c_+$ and one $c_-$, so we only need to fix the signs here.) Hence, the following function.  \n",
    "\n",
    "We will now write a function in NRPy+ similar to the one used in the old `GiRaFFE`, allowing us to generate the expressions with less need to copy-and-paste code; the key difference is that this one will be in Python, and generate optimized C code integrated into the rest of the operations. Notice that since we eliminated the dependence on velocities, none of the input quantities are different on either side of the face. So, this function won't really do much besides guarantee that `cmax` and `cmin` are positive, but we'll leave the machinery here since it is likely to be a useful guide to somebody who wants to something similar. We use the same technique as above to replace the `if` statements inherent to the `MAX()` and `MIN()` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll write this as a function, and call it within HLLE_solver, below.\n",
    "def find_cmax_cmin(flux_dirn,gamma_faceDD,beta_faceU,alpha_face):\n",
    "    # Inputs:  flux direction flux_dirn, Inverse metric gamma_faceUU, shift beta_faceU,\n",
    "    #          lapse alpha_face, metric determinant gammadet_face\n",
    "    # Outputs: maximum and minimum characteristic speeds cmax and cmin\n",
    "    # First, we need to find the characteristic speeds on each face\n",
    "    gamma_faceUU,unusedgammaDET = ixp.generic_matrix_inverter3x3(gamma_faceDD)\n",
    "    # Original needed for GRMHD\n",
    "#     find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])\n",
    "#     cpr = cplus\n",
    "#     cmr = cminus\n",
    "#     find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])\n",
    "#     cpl = cplus\n",
    "#     cml = cminus\n",
    "    find_cp_cm(alpha_face,beta_faceU[flux_dirn],gamma_faceUU[flux_dirn][flux_dirn])\n",
    "    cp = cplus\n",
    "    cm = cminus\n",
    "\n",
    "    # The following algorithms have been verified with random floats:\n",
    "\n",
    "    global cmax,cmin\n",
    "    # Now, we need to set cmax to the larger of cpr,cpl, and 0\n",
    "\n",
    "    import Min_Max_and_Piecewise_Expressions as noif\n",
    "    cmax = noif.max_noif(cp,sp.sympify(0))\n",
    "\n",
    "    # And then, set cmin to the smaller of cmr,cml, and 0\n",
    "    cmin = -noif.min_noif(cm,sp.sympify(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4:  Code Validation against `GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds` NRPy+ Module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for the `GiRaFFE` evolution equations and auxiliary quantities we intend to use between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Characteristic_Speeds.py) module.\n",
    "\n",
    "This first validation directly compares the sympy expressions. This is generally quicker and more reliable, but might overlook some complexities in implementing the C code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Checking the flux in direction 0\n",
      "Checking the flux in direction 1\n",
      "Checking the flux in direction 2\n",
      "ALL TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "all_passed=True\n",
    "def comp_func(expr1,expr2,basename,prefixname2=\"Sf.\"):\n",
    "    if str(expr1-expr2)!=\"0\":\n",
    "        print(basename+\" - \"+prefixname2+basename+\" = \"+ str(expr1-expr2))\n",
    "        all_passed=False\n",
    "\n",
    "def gfnm(basename,idx1,idx2=None,idx3=None):\n",
    "    if idx2 is None:\n",
    "        return basename+\"[\"+str(idx1)+\"]\"\n",
    "    if idx3 is None:\n",
    "        return basename+\"[\"+str(idx1)+\"][\"+str(idx2)+\"]\"\n",
    "    return basename+\"[\"+str(idx1)+\"][\"+str(idx2)+\"][\"+str(idx3)+\"]\"\n",
    "\n",
    "# We'll first declare the inputs to the function.\n",
    "gammaDD = ixp.declarerank2(\"gammaDD\",\"sym01\",DIM=3)\n",
    "betaU = ixp.declarerank1(\"betaU\",DIM=3)\n",
    "alpha = sp.sympify(\"alpha\")\n",
    "import GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds as chsp\n",
    "\n",
    "for flux_dirn in range(3):\n",
    "    expr_list = []\n",
    "    exprcheck_list = []\n",
    "    namecheck_list = []\n",
    "\n",
    "    print(\"Checking the flux in direction \"+str(flux_dirn))\n",
    "    find_cmax_cmin(flux_dirn,gammaDD,betaU,alpha)\n",
    "    chsp.find_cmax_cmin(flux_dirn,gammaDD,betaU,alpha)\n",
    "\n",
    "    namecheck_list.extend([\"cmax\",\"cmin\"])\n",
    "    exprcheck_list.extend([chsp.cmax,chsp.cmin])\n",
    "    expr_list.extend([cmax,cmin])\n",
    "\n",
    "    for i in range(len(expr_list)):\n",
    "        comp_func(expr_list[i],exprcheck_list[i],namecheck_list[i])\n",
    "\n",
    "if all_passed:\n",
    "    print(\"ALL TESTS PASSED!\")\n",
    "else:\n",
    "    print(\"ERROR: AT LEAST ONE TEST DID NOT PASS\")\n",
    "    sys.exit(1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='derive_speed'></a>\n",
    "\n",
    "# Step 5: Complete Derivation of the Wave Speeds \\[Back to [top](#toc)\\]\n",
    "$$\\label{derive_speed}$$\n",
    "\n",
    "This computes phase speeds in the direction given by `flux_dirn`. Note that we replace the full dispersion relation with a simpler one, which overestimates the maximum speeds by a factor of ~2. See full discussion around Eqs. 49 and 50 in [Duez, et al.](http://arxiv.org/pdf/astro-ph/0503420.pdf). In summary, we solve the dispersion relation (in, e.g., the $x$-direction) with a wave vector of $k_\\mu = (-\\omega,k_x,0,0)$. So, we solve the approximate dispersion relation $\\omega_{\\rm cm}^2 = [v_A^2 + c_s^2 (1-v_A^2)]k_{\\rm cm}^2$ for the wave speed $\\omega/k_x$, where the sound speed $c_s = \\sqrt{\\Gamma P/(h \\rho_0)}$, the Alfv&eacute;n speed $v_A = 1$ (in GRFFE), $\\omega_{\\rm cm} = -k_\\mu k^\\mu$ is the frequency in the comoving frame, $k_{\\rm cm}^2 = K_\\mu K^\\mu$ is the wavenumber squared in the comoving frame, and $K_\\mu = (g_{\\mu\\nu} + u_\\mu u_\\nu)k^\\nu$ is the part of the wave vector normal to the four-velocity $u^\\mu$. See below for a complete derivation.\n",
    "\n",
    "What follows is a complete derivation of the quadratic we solve. We start from the following relations:\n",
    "\\begin{align}\n",
    "w_{\\rm cm} &= (-k_0 u^0 - k_x u^x) \\\\\n",
    "k_{\\rm cm}^2 &= K_{\\mu} K^{\\mu}, \\\\\n",
    "K_{\\mu} K^{\\mu} &= (g_{\\mu a} + u_{\\mu} u_a) k^a  g^{\\mu b} [ (g_{c b} + u_c u_b) k^c ] \\\\\n",
    "\\end{align}\n",
    "The last term of the above can be written as follow:\n",
    "$$\n",
    "(g_{c b} + u_{c} u_{b}) k^c = (\\delta^{\\mu}_c + u_c u^{\\mu} ) k^c\n",
    "$$\n",
    "\n",
    "Then,\n",
    "\\begin{align}\n",
    "K_{\\mu} K^{\\mu} &= (g_{\\mu a} + u_{\\mu} u_a) k^a  g^{\\mu b} [ (g_{c b} + u_c u_b) k^c ] \\\\\n",
    "                 &= (g_{\\mu a} + u_{\\mu} u_a) k^a  (\\delta^{\\mu}_c + u_c u^{\\mu} ) k^c \\\\\n",
    "                 &=[(g_{\\mu a} + u_{\\mu} u_a) \\delta^{\\mu}_c + (g_{\\mu a} + u_{\\mu} u_a) u_c u^{\\mu} ] k^c k^a \\\\\n",
    "                 &=[(g_{c a} + u_c u_a) + (u_c u_a -  u_a u_c] k^c k^a \\\\\n",
    "                 &=(g_{c a} + u_c u_a) k^c k^a \\\\\n",
    "                 &= k_a k^a + u^c u^a k_c k_a \\\\\n",
    "k^a = g^{\\mu a} k_{\\mu} &= g^{0 a} k_0 + g^{x a} k_x \\\\\n",
    "k_a k^a &= k_0 g^{0 0} k_0 + k_x k_0 g^{0 x} + g^{x 0} k_0 k_x + g^{x x} k_x k_x \\\\\n",
    "         &= g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2 \\\\\n",
    "u^c u^a k_c k_a &= (u^0 k_0 + u^x k_x) (u^0 k_0 + u^x k_x) = (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 \\\\\n",
    "(k_0 u^0)^2  + 2 k_x u^x k_0 u^0 + (k_x u^x)^2 &= v_0^2 [ (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 + g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2] \\\\\n",
    "(1-v_0^2) (u^0 k_0 + u^x k_x)^2 &= v_0^2 (g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2) \\\\\n",
    "(1-v_0^2) (u^0 k_0/k_x + u^x)^2 &= v_0^2 (g^{00} (k_0/k_x)^2 + 2 g^{x0} k_0/k_x + g^{xx}) \\\\\n",
    "(1-v_0^2) (u^0 X + u^x)^2 &= v_0^2 (g^{00} X^2 + 2 g^{x0} X + g^{xx}) \\\\\n",
    "(1-v_0^2) ((u^0)^2 X^2 + 2 u^x (u^0) X + (u^x)^2) &= v_0^2 (g^{00} X^2 + 2 g^{x0} X + g^{xx}) \\\\\n",
    "0 &= X^2 ( (1-v_0^2) (u^0)^2 - v_0^2 g^{00}) + X (2 u^x u^0 (1-v_0^2) - 2 v_0^2 g^{x0}) + (1-v_0^2) (u^x)^2 - v_0^2 g^{xx} \\\\\n",
    "a &= (1-v_0^2) (u^0)^2 - v_0^2 g^{00} = (1-v_0^2) (u^0)^2 + v_0^2/\\alpha^2 \\leftarrow {\\rm VERIFIED} \\\\\n",
    "b &= 2 u^x u^0 (1-v_0^2) - 2 v_0^2 \\beta^x/\\alpha^2 \\leftarrow {\\rm VERIFIED,\\ } X\\rightarrow -X, {\\rm because\\ } X = -w/k_1, {\\rm \\ and\\ we\\ are\\ solving\\ for} -X. \\\\\n",
    "c &= (1-v_0^2) (u^x)^2 - v_0^2 (\\gamma^{xx}\\psi^{-4} - (\\beta^x/\\alpha)^2) \\leftarrow {\\rm VERIFIED} \\\\\n",
    "v_0^2 &= v_A^2 + c_s^2 (1 - v_A^2) \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: 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-GiRaFFE_NRPy-Characteristic_Speeds.pdf](Tutorial-GiRaFFE_NRPy-Characteristic_Speeds.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": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[NbConvertApp] ERROR | Error while converting 'Tutorial-GiRaFFE_NRPy-Characteristic_Speeds.ipynb'\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/nbconvertapp.py\", line 418, in export_single_notebook\n",
      "    output, resources = self.exporter.from_filename(notebook_filename, resources=resources)\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/exporters/exporter.py\", line 181, in from_filename\n",
      "    return self.from_file(f, resources=resources, **kw)\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/exporters/exporter.py\", line 199, in from_file\n",
      "    return self.from_notebook_node(nbformat.read(file_stream, as_version=4), resources=resources, **kw)\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/exporters/latex.py\", line 77, in from_notebook_node\n",
      "    return super().from_notebook_node(nb, resources, **kw)\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/exporters/templateexporter.py\", line 384, in from_notebook_node\n",
      "    output = self.template.render(nb=nb_copy, resources=resources)\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/jinja2/environment.py\", line 1090, in render\n",
      "    self.environment.handle_exception()\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/jinja2/environment.py\", line 832, in handle_exception\n",
      "    reraise(*rewrite_traceback_stack(source=source))\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/jinja2/_compat.py\", line 28, in reraise\n",
      "    raise value.with_traceback(tb)\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/index.tex.j2\", line 8, in top-level template code\n",
      "    ((* extends cell_style *))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/style_jupyter.tex.j2\", line 176, in top-level template code\n",
      "    \\prompt{(((prompt)))}{(((prompt_color)))}{(((execution_count)))}{(((extra_space)))}\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/base.tex.j2\", line 7, in top-level template code\n",
      "    ((*- extends 'document_contents.tex.j2' -*))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/document_contents.tex.j2\", line 51, in top-level template code\n",
      "    ((*- block figure scoped -*))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/display_priority.j2\", line 5, in top-level template code\n",
      "    ((*- extends 'null.j2' -*))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/null.j2\", line 30, in top-level template code\n",
      "    ((*- block body -*))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/base.tex.j2\", line 206, in block \"body\"\n",
      "    ((( super() )))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/null.j2\", line 32, in block \"body\"\n",
      "    ((*- block any_cell scoped -*))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/null.j2\", line 85, in block \"any_cell\"\n",
      "    ((*- block markdowncell scoped-*)) ((*- endblock markdowncell -*))\n",
      "  File \"/home/penelson/nrpyvirtualenv/share/jupyter/nbconvert/templates/latex/document_contents.tex.j2\", line 68, in block \"markdowncell\"\n",
      "    ((( cell.source | citation2latex | strip_files_prefix | convert_pandoc('markdown+tex_math_double_backslash', 'json',extra_args=[]) | resolve_references | convert_pandoc('json','latex'))))\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/filters/pandoc.py\", line 24, in convert_pandoc\n",
      "    return pandoc(source, from_format, to_format, extra_args=extra_args)\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/utils/pandoc.py\", line 52, in pandoc\n",
      "    check_pandoc_version()\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/utils/pandoc.py\", line 100, in check_pandoc_version\n",
      "    v = get_pandoc_version()\n",
      "  File \"/home/penelson/nrpyvirtualenv/lib/python3.8/site-packages/nbconvert/utils/pandoc.py\", line 77, in get_pandoc_version\n",
      "    raise PandocMissing()\n",
      "nbconvert.utils.pandoc.PandocMissing: Pandoc wasn't found.\n",
      "Please check that pandoc is installed:\n",
      "https://pandoc.org/installing.html\n"
     ]
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'pdflatex'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-2cc4dcfe5a63>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcmdline_helper\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mcmd\u001b[0m    \u001b[0;31m# NRPy+: Multi-platform Python command-line interface\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mcmd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moutput_Jupyter_notebook_to_LaTeXed_PDF\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Tutorial-GiRaFFE_NRPy-Characteristic_Speeds\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlocation_of_template_file\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"..\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/nrpyvirtualenv/nrpytutorial/cmdline_helper.py\u001b[0m in \u001b[0;36moutput_Jupyter_notebook_to_LaTeXed_PDF\u001b[0;34m(notebookname, location_of_template_file, verbose)\u001b[0m\n\u001b[1;32m    194\u001b[0m                          +r\" --log-level='WARN' \"+notebookname+\".ipynb\",verbose=False)\n\u001b[1;32m    195\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0m_i\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# _i is an unused variable.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 196\u001b[0;31m         \u001b[0mExecute_input_string\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr\"pdflatex -interaction=batchmode \"\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mnotebookname\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m\".tex\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    197\u001b[0m     \u001b[0mdelete_existing_files\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnotebookname\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m\".out \"\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mnotebookname\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m\".aux \"\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mnotebookname\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m\".log\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    198\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/nrpyvirtualenv/nrpytutorial/cmdline_helper.py\u001b[0m in \u001b[0;36mExecute_input_string\u001b[0;34m(input_string, file_to_redirect_stdout, verbose)\u001b[0m\n\u001b[1;32m    161\u001b[0m     \u001b[0mfilename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"tmp.txt\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    162\u001b[0m     \u001b[0;32mwith\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'w'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mwriter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'rb'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbuffering\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mreader\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile_to_redirect_stdout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'wb'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mrdirect\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 163\u001b[0;31m         \u001b[0mprocess\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msubprocess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstdout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mrdirect\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mwriter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    164\u001b[0m         \u001b[0;32mwhile\u001b[0m \u001b[0mprocess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpoll\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    165\u001b[0m             \u001b[0;31m# https://stackoverflow.com/questions/21689365/python-3-typeerror-must-be-str-not-bytes-with-sys-stdout-write/21689447\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.8/subprocess.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)\u001b[0m\n\u001b[1;32m    852\u001b[0m                             encoding=encoding, errors=errors)\n\u001b[1;32m    853\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 854\u001b[0;31m             self._execute_child(args, executable, preexec_fn, close_fds,\n\u001b[0m\u001b[1;32m    855\u001b[0m                                 \u001b[0mpass_fds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcwd\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0menv\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    856\u001b[0m                                 \u001b[0mstartupinfo\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcreationflags\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshell\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.8/subprocess.py\u001b[0m in \u001b[0;36m_execute_child\u001b[0;34m(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)\u001b[0m\n\u001b[1;32m   1700\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0merrno_num\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1701\u001b[0m                         \u001b[0merr_msg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstrerror\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merrno_num\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1702\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mchild_exception_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merrno_num\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merr_msg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merr_filename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1703\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0mchild_exception_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr_msg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1704\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'pdflatex'"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GiRaFFE_NRPy-Characteristic_Speeds\",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
}