Tutorial-u0_smallb_Poynting-Cartesian.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",
    "# Computing the 4-Velocity Time-Component $u^0$, the Magnetic Field Measured by a Comoving Observer $b^{\\mu}$, and the Poynting Vector $S^i$\n",
    "\n",
    "## Authors: Zach Etienne & Patrick Nelson\n",
    "\n",
    "## This notebook computes key parameters in relativistic MHD, utilizing the [BSSN.ADMBSSN_tofrom_4metric](../edit/BSSN/ADMBSSN_tofrom_4metric.py) NRPy+ module and the Valencia 3-velocity.  The computations are then validated against the [`u0_smallb_Poynting__Cartesian.py`](../edit/u0_smallb_Poynting__Cartesian/u0_smallb_Poynting__Cartesian.py) NRPy+ module.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated against a trusted code (the hand-written smallbPoynET in WVUThorns_diagnostics, which itself is based on expressions in IllinoisGRMHD... which was validated against the original GRMHD code of the Illinois NR group)\n",
    "\n",
    "### NRPy+ Source Code for this module: [u0_smallb_Poynting__Cartesian.py](../edit/u0_smallb_Poynting__Cartesian/u0_smallb_Poynting__Cartesian.py)\n",
    "\n",
    "[comment]: <> (Introduction: TODO)"
   ]
  },
  {
   "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](#u0bu): Computing $u^0$ and $b^{\\mu}$\n",
    "    1. [Step 1.a](#4metric): Compute the 4-metric $g_{\\mu\\nu}$ and its inverse $g^{\\mu\\nu}$ from the ADM 3+1 variables, using the [`BSSN.ADMBSSN_tofrom_4metric`](../edit/BSSN/ADMBSSN_tofrom_4metric.py) ([**tutorial**](Tutorial-ADMBSSN_tofrom_4metric.ipynb)) NRPy+ module\n",
    "    1. [Step 1.b](#u0): Compute $u^0$ from the Valencia 3-velocity\n",
    "    1. [Step 1.c](#uj): Compute $u_j$ from $u^0$, the Valencia 3-velocity, and $g_{\\mu\\nu}$\n",
    "    1. [Step 1.d](#gamma): Compute $\\gamma=$ `gammaDET` from the ADM 3+1 variables\n",
    "    1. [Step 1.e](#beta): Compute $b^\\mu$\n",
    "1. [Step 2](#poynting_flux): Defining the Poynting Flux Vector $S^{i}$\n",
    "    1. [Step 2.a](#g): Computing $g^{i\\nu}$\n",
    "    1. [Step 2.b](#s): Computing $S^{i}$\n",
    "1. [Step 3](#code_validation): Code Validation against `u0_smallb_Poynting__Cartesian` NRPy+ module\n",
    "1. [Step 4](#appendix): Appendix: Proving Eqs. 53 and 56 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf)\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='u0bu'></a>\n",
    "\n",
    "# Step 1: Computing $u^0$ and $b^{\\mu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{u0bu}$$\n",
    "\n",
    "First some definitions. The spatial components of $b^{\\mu}$ are simply the magnetic field as measured by an observer comoving with the plasma $B^{\\mu}_{\\rm (u)}$, divided by $\\sqrt{4\\pi}$. In addition, in the ideal MHD limit, $B^{\\mu}_{\\rm (u)}$ is orthogonal to the plasma 4-velocity $u^\\mu$, which sets the $\\mu=0$ component. \n",
    "\n",
    "Note also that $B^{\\mu}_{\\rm (u)}$ is related to the magnetic field as measured by a *normal* observer $B^i$ via a simple projection (Eq 21 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf)), which results in the expressions (Eqs 23 and 24 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf)):\n",
    "\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",
    "\n",
    "$B^i$ is related to the actual magnetic field evaluated in IllinoisGRMHD, $\\tilde{B}^i$ via\n",
    "\n",
    "$$B^i = \\frac{\\tilde{B}^i}{\\gamma},$$\n",
    "\n",
    "where $\\gamma$ is the determinant of the spatial 3-metric.\n",
    "\n",
    "The above expressions will require that we compute\n",
    "1. the 4-metric $g_{\\mu\\nu}$ from the ADM 3+1 variables\n",
    "1. $u^0$ from the Valencia 3-velocity\n",
    "1. $u_j$ from $u^0$, the Valencia 3-velocity, and $g_{\\mu\\nu}$\n",
    "1. $\\gamma$ from the ADM 3+1 variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='4metric'></a>\n",
    "\n",
    "## Step 1.a: Compute the 4-metric $g_{\\mu\\nu}$ and its inverse $g^{\\mu\\nu}$ from the ADM 3+1 variables, using the [`BSSN.ADMBSSN_tofrom_4metric`](../edit/BSSN/ADMBSSN_tofrom_4metric.py) ([**tutorial**](Tutorial-ADMBSSN_tofrom_4metric.ipynb)) NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{4metric}$$\n",
    "\n",
    "We are given $\\gamma_{ij}$, $\\alpha$, and $\\beta^i$ from ADMBase, so let's first compute \n",
    "\n",
    "$$\n",
    "g_{\\mu\\nu} = \\begin{pmatrix} \n",
    "-\\alpha^2 + \\beta^k \\beta_k & \\beta_i \\\\\n",
    "\\beta_j & \\gamma_{ij}\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:18:33.811044Z",
     "iopub.status.busy": "2021-03-07T17:18:33.809895Z",
     "iopub.status.idle": "2021-03-07T17:18:34.206377Z",
     "shell.execute_reply": "2021-03-07T17:18:34.206911Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Initialize needed Python/NRPy+ modules\n",
    "import sympy as sp                         # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par             # NRPy+: Parameter interface\n",
    "import indexedexp as ixp                   # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "from outputC import outputC                # NRPy+: Basic C code output functionality\n",
    "import BSSN.ADMBSSN_tofrom_4metric as AB4m # NRPy+: ADM/BSSN <-> 4-metric conversions\n",
    "\n",
    "# Set spatial dimension = 3\n",
    "DIM=3\n",
    "\n",
    "thismodule = \"smallbPoynET\"\n",
    "\n",
    "# Step 1.a: Compute the 4-metric $g_{\\mu\\nu}$ and its inverse\n",
    "#           $g^{\\mu\\nu}$ from the ADM 3+1 variables, using the\n",
    "#           BSSN.ADMBSSN_tofrom_4metric NRPy+ module\n",
    "import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "gammaDD,betaU,alpha = AB4m.setup_ADM_quantities(\"ADM\")\n",
    "AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "g4DD = AB4m.g4DD\n",
    "AB4m.g4UU_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "g4UU = AB4m.g4UU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='u0'></a>\n",
    "\n",
    "## Step 1.b: Compute $u^0$ from the Valencia 3-velocity \\[Back to [top](#toc)\\]\n",
    "$$\\label{u0}$$\n",
    "\n",
    "According to Eqs. 9-11 of [the IllinoisGRMHD paper](https://arxiv.org/pdf/1501.07276.pdf), the Valencia 3-velocity $v^i_{(n)}$ is related to the 4-velocity $u^\\mu$ via\n",
    "\n",
    "\\begin{align}\n",
    "\\alpha v^i_{(n)} &= \\frac{u^i}{u^0} + \\beta^i \\\\\n",
    "\\implies u^i &= u^0 \\left(\\alpha v^i_{(n)} - \\beta^i\\right)\n",
    "\\end{align}\n",
    "\n",
    "Defining $v^i = \\frac{u^i}{u^0}$, we get\n",
    "\n",
    "$$v^i = \\alpha v^i_{(n)} - \\beta^i,$$\n",
    "\n",
    "and in terms of this variable we get\n",
    "\n",
    "\\begin{align}\n",
    "g_{00} \\left(u^0\\right)^2 + 2 g_{0i} u^0 u^i + g_{ij} u^i u^j &= \\left(u^0\\right)^2 \\left(g_{00} + 2 g_{0i} v^i + g_{ij} v^i v^j\\right)\\\\\n",
    "\\implies u^0 &= \\pm \\sqrt{\\frac{-1}{g_{00} + 2 g_{0i} v^i + g_{ij} v^i v^j}} \\\\\n",
    "&= \\pm \\sqrt{\\frac{-1}{(-\\alpha^2 + \\beta^2) + 2 \\beta_i v^i + \\gamma_{ij} v^i v^j}} \\\\\n",
    "&= \\pm \\sqrt{\\frac{1}{\\alpha^2 - \\gamma_{ij}\\left(\\beta^i + v^i\\right)\\left(\\beta^j + v^j\\right)}}\\\\\n",
    "&= \\pm \\sqrt{\\frac{1}{\\alpha^2 - \\alpha^2 \\gamma_{ij}v^i_{(n)}v^j_{(n)}}}\\\\\n",
    "&= \\pm \\frac{1}{\\alpha}\\sqrt{\\frac{1}{1 - \\gamma_{ij}v^i_{(n)}v^j_{(n)}}}\n",
    "\\end{align}\n",
    "\n",
    "Generally speaking, numerical errors will occasionally drive expressions under the radical to either negative values or potentially enormous values (corresponding to enormous Lorentz factors). Thus a reliable approach for computing $u^0$ requires that we first rewrite the above expression in terms of the Lorentz factor squared: $\\Gamma^2=\\left(\\alpha u^0\\right)^2$:\n",
    "\\begin{align}\n",
    "u^0 &= \\pm \\frac{1}{\\alpha}\\sqrt{\\frac{1}{1 - \\gamma_{ij}v^i_{(n)}v^j_{(n)}}}\\\\\n",
    "\\implies \\left(\\alpha u^0\\right)^2 &= \\frac{1}{1 - \\gamma_{ij}v^i_{(n)}v^j_{(n)}} \\\\\n",
    "\\implies \\gamma_{ij}v^i_{(n)}v^j_{(n)} &= 1 - \\frac{1}{\\left(\\alpha u^0\\right)^2} \\\\\n",
    "&= 1 - \\frac{1}{\\Gamma^2}\n",
    "\\end{align}\n",
    "\n",
    "In order for the bottom expression to hold true, the left-hand side must be between 0 and 1. Again, this is not guaranteed due to the appearance of numerical errors. In fact, a robust algorithm will not allow $\\Gamma^2$ to become too large (which might contribute greatly to the stress energy of a given gridpoint), so let's define $\\Gamma_{\\rm max}$, the largest allowed Lorentz factor. \n",
    "\n",
    "Then our algorithm for computing $u^0$ is as follows:\n",
    "\n",
    "If\n",
    "$$R=\\gamma_{ij}v^i_{(n)}v^j_{(n)}>1 - \\frac{1}{\\Gamma_{\\rm max}^2},$$ \n",
    "then adjust the 3-velocity $v^i$ as follows:\n",
    "\n",
    "$$v^i_{(n)} = \\sqrt{\\frac{1 - \\frac{1}{\\Gamma_{\\rm max}^2}}{R}}v^i_{(n)}.$$\n",
    "\n",
    "After this rescaling, we are then guaranteed that if $R$ is recomputed, it will be set to its ceiling value $R=R_{\\rm max} = 1 - \\frac{1}{\\Gamma_{\\rm max}^2}$.\n",
    "\n",
    "Then, regardless of whether the ceiling on $R$ was applied, $u^0$ can be safely computed via\n",
    "\n",
    "$$\n",
    "u^0 = \\frac{1}{\\alpha \\sqrt{1-R}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:18:34.291188Z",
     "iopub.status.busy": "2021-03-07T17:18:34.275195Z",
     "iopub.status.idle": "2021-03-07T17:18:34.340673Z",
     "shell.execute_reply": "2021-03-07T17:18:34.341284Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "/* Function for computing u^0 from Valencia 3-velocity. */\n",
      "/* Inputs: ValenciavU[], alpha, gammaDD[][], GAMMA_SPEED_LIMIT (C parameter) */\n",
      "/* Output: u0=u^0 and velocity-limited ValenciavU[] */\n",
      "\n",
      "const double R = ((ValenciavU0)*(ValenciavU0))*gammaDD00 + 2*ValenciavU0*ValenciavU1*gammaDD01 + 2*ValenciavU0*ValenciavU2*gammaDD02 + ((ValenciavU1)*(ValenciavU1))*gammaDD11 + 2*ValenciavU1*ValenciavU2*gammaDD12 + ((ValenciavU2)*(ValenciavU2))*gammaDD22;\n",
      "const double Rmax = 1 - 1/((GAMMA_SPEED_LIMIT)*(GAMMA_SPEED_LIMIT));\n",
      "if(R <= Rmax) {\n",
      "  u0 = 1/(alpha*sqrt(-((ValenciavU0)*(ValenciavU0))*gammaDD00 - 2*ValenciavU0*ValenciavU1*gammaDD01 - 2*ValenciavU0*ValenciavU2*gammaDD02 - ((ValenciavU1)*(ValenciavU1))*gammaDD11 - 2*ValenciavU1*ValenciavU2*gammaDD12 - ((ValenciavU2)*(ValenciavU2))*gammaDD22 + 1));\n",
      "}\n",
      " else {\n",
      "  const double tmprescale_1 = sqrt((1 - 1/((GAMMA_SPEED_LIMIT)*(GAMMA_SPEED_LIMIT)))/(((ValenciavU0)*(ValenciavU0))*gammaDD00 + 2*ValenciavU0*ValenciavU1*gammaDD01 + 2*ValenciavU0*ValenciavU2*gammaDD02 + ((ValenciavU1)*(ValenciavU1))*gammaDD11 + 2*ValenciavU1*ValenciavU2*gammaDD12 + ((ValenciavU2)*(ValenciavU2))*gammaDD22));\n",
      "  ValenciavU0 = ValenciavU0*tmprescale_1;\n",
      "  ValenciavU1 = ValenciavU1*tmprescale_1;\n",
      "  ValenciavU2 = ValenciavU2*tmprescale_1;\n",
      "  u0 = fabs(GAMMA_SPEED_LIMIT)/alpha;\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"ValenciavU\",DIM=3)\n",
    "\n",
    "# Step 1: Compute R = 1 - 1/max(Gamma)\n",
    "R = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        R += gammaDD[i][j]*ValenciavU[i]*ValenciavU[j]\n",
    "\n",
    "GAMMA_SPEED_LIMIT = par.Cparameters(\"REAL\",thismodule,\"GAMMA_SPEED_LIMIT\",10.0) # Default value based on\n",
    "                                                                                # IllinoisGRMHD.\n",
    "                                                                                # GiRaFFE default = 2000.0\n",
    "Rmax = 1 - 1/(GAMMA_SPEED_LIMIT*GAMMA_SPEED_LIMIT)\n",
    "\n",
    "rescaledValenciavU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    rescaledValenciavU[i] = ValenciavU[i]*sp.sqrt(Rmax/R)\n",
    "\n",
    "rescaledu0 = 1/(alpha*sp.sqrt(1-Rmax))\n",
    "regularu0 =  1/(alpha*sp.sqrt(1-R))\n",
    "\n",
    "computeu0_Cfunction  = \"\"\"\n",
    "/* Function for computing u^0 from Valencia 3-velocity. */\n",
    "/* Inputs: ValenciavU[], alpha, gammaDD[][], GAMMA_SPEED_LIMIT (C parameter) */\n",
    "/* Output: u0=u^0 and velocity-limited ValenciavU[] */\\n\\n\"\"\"\n",
    "\n",
    "computeu0_Cfunction += outputC([R,Rmax],[\"const double R\",\"const double Rmax\"],\"returnstring\",\n",
    "                               params=\"includebraces=False,CSE_varprefix=tmpR,outCverbose=False\")\n",
    "\n",
    "computeu0_Cfunction += \"if(R <= Rmax) \"\n",
    "computeu0_Cfunction += outputC(regularu0,\"u0\",\"returnstring\",\n",
    "                               params=\"includebraces=True,CSE_varprefix=tmpnorescale,outCverbose=False\")\n",
    "computeu0_Cfunction += \" else \"\n",
    "computeu0_Cfunction += outputC([rescaledValenciavU[0],rescaledValenciavU[1],rescaledValenciavU[2],rescaledu0],\n",
    "                               [\"ValenciavU0\",\"ValenciavU1\",\"ValenciavU2\",\"u0\"],\"returnstring\",\n",
    "                               params=\"includebraces=True,CSE_varprefix=tmprescale,outCverbose=False\")\n",
    "\n",
    "print(computeu0_Cfunction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='uj'></a>\n",
    "\n",
    "## Step 1.c: Compute $u_j$ from $u^0$, the Valencia 3-velocity, and $g_{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{uj}$$\n",
    "\n",
    "The basic equation is\n",
    "\n",
    "\\begin{align}\n",
    "u_j &= g_{\\mu j} u^{\\mu} \\\\\n",
    "&= g_{0j} u^0 + g_{ij} u^i \\\\\n",
    "&= \\beta_j u^0 + \\gamma_{ij} u^i \\\\\n",
    "&= \\beta_j u^0 + \\gamma_{ij} u^0 \\left(\\alpha v^i_{(n)} - \\beta^i\\right) \\\\\n",
    "&= u^0 \\left(\\beta_j + \\gamma_{ij} \\left(\\alpha v^i_{(n)} - \\beta^i\\right) \\right)\\\\\n",
    "&= \\alpha u^0 \\gamma_{ij} v^i_{(n)} \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:18:34.352644Z",
     "iopub.status.busy": "2021-03-07T17:18:34.351433Z",
     "iopub.status.idle": "2021-03-07T17:18:34.354734Z",
     "shell.execute_reply": "2021-03-07T17:18:34.354155Z"
    }
   },
   "outputs": [],
   "source": [
    "u0 = par.Cparameters(\"REAL\",thismodule,\"u0\",1e300) # Will be overwritten in C code. Set to crazy value to ensure this.\n",
    "\n",
    "uD = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        uD[j] += alpha*u0*gammaDD[i][j]*ValenciavU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='beta'></a>\n",
    "\n",
    "## Step 1.d: Compute $b^\\mu$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{beta}$$\n",
    "\n",
    "We compute $b^\\mu$ from the above expressions:\n",
    "\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",
    "\n",
    "$B^i$ is exactly equal to the $B^i$ evaluated in IllinoisGRMHD/GiRaFFE.\n",
    "\n",
    "Pulling this together, we currently have available as input:\n",
    "+ $\\tilde{B}^i$\n",
    "+ $u_j$\n",
    "+ $u^0$,\n",
    "\n",
    "with the goal of outputting now $b^\\mu$ and $b^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:18:34.386008Z",
     "iopub.status.busy": "2021-03-07T17:18:34.385080Z",
     "iopub.status.idle": "2021-03-07T17:18:34.387571Z",
     "shell.execute_reply": "2021-03-07T17:18:34.388076Z"
    }
   },
   "outputs": [],
   "source": [
    "M_PI = par.Cparameters(\"#define\",thismodule,\"M_PI\",\"\")\n",
    "BU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"BU\",DIM=3)\n",
    "\n",
    "# uBcontraction = u_i B^i\n",
    "uBcontraction = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    uBcontraction += uD[i]*BU[i]\n",
    "\n",
    "# uU = 3-vector representing u^i = u^0 \\left(\\alpha v^i_{(n)} - \\beta^i\\right)\n",
    "uU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    uU[i] = u0*(alpha*ValenciavU[i] - betaU[i])\n",
    "\n",
    "smallb4U = ixp.zerorank1(DIM=4)\n",
    "smallb4U[0] = uBcontraction/(alpha*sp.sqrt(4*M_PI))\n",
    "for i in range(DIM):\n",
    "    smallb4U[1+i] = (BU[i] + uBcontraction*uU[i])/(alpha*u0*sp.sqrt(4*M_PI))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='poynting_flux'></a>\n",
    "\n",
    "# Step 2: Defining the Poynting Flux Vector $S^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{poynting_flux}$$\n",
    "\n",
    "The Poynting flux is defined in Eq. 11 of [Kelly *et al.*](https://arxiv.org/pdf/1710.02132.pdf) (note that we choose the minus sign convention so that the Poynting luminosity across a spherical shell is $L_{\\rm EM} = \\int (-\\alpha T^i_{\\rm EM\\ 0}) \\sqrt{\\gamma} d\\Omega = \\int S^r \\sqrt{\\gamma} d\\Omega$, as in [Farris *et al.*](https://arxiv.org/pdf/1207.3354.pdf):\n",
    "\n",
    "$$\n",
    "S^i = -\\alpha T^i_{\\rm EM\\ 0} = -\\alpha\\left(b^2 u^i u_0 + \\frac{1}{2} b^2 g^i{}_0 - b^i b_0\\right)\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='s'></a>\n",
    "\n",
    "## Step 2.a: Computing $S^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{s}$$\n",
    "\n",
    "Given $g^{\\mu\\nu}$ computed above, we focus first on the $g^i{}_{0}$ term by computing \n",
    "$$\n",
    "g^\\mu{}_\\delta = g^{\\mu\\nu} g_{\\nu \\delta},\n",
    "$$\n",
    "and then the rest of the Poynting flux vector can be immediately computed from quantities defined above:\n",
    "$$\n",
    "S^i = -\\alpha T^i_{\\rm EM\\ 0} = -\\alpha\\left(b^2 u^i u_0 + \\frac{1}{2} b^2 g^i{}_0 - b^i b_0\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:18:34.437504Z",
     "iopub.status.busy": "2021-03-07T17:18:34.427172Z",
     "iopub.status.idle": "2021-03-07T17:18:34.457798Z",
     "shell.execute_reply": "2021-03-07T17:18:34.457115Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a.i: compute g^\\mu_\\delta:\n",
    "g4UD = ixp.zerorank2(DIM=4)\n",
    "for mu in range(4):\n",
    "    for delta in range(4):\n",
    "        for nu in range(4):\n",
    "            g4UD[mu][delta] += g4UU[mu][nu]*g4DD[nu][delta]\n",
    "\n",
    "# Step 2.a.ii: compute b_{\\mu}\n",
    "smallb4D = ixp.zerorank1(DIM=4)\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        smallb4D[mu] += g4DD[mu][nu]*smallb4U[nu]\n",
    "\n",
    "# Step 2.a.iii: compute u_0 = g_{mu 0} u^{mu} = g4DD[0][0]*u0 + g4DD[i][0]*uU[i]\n",
    "u_0 = g4DD[0][0]*u0\n",
    "for i in range(DIM):\n",
    "    u_0 += g4DD[i+1][0]*uU[i]\n",
    "\n",
    "# Step 2.a.iv: compute b^2, setting b^2 = smallb2etk, as gridfunctions with base names ending in a digit\n",
    "#          are forbidden in NRPy+.\n",
    "smallb2etk = sp.sympify(0)\n",
    "for mu in range(4):\n",
    "    smallb2etk += smallb4U[mu]*smallb4D[mu]\n",
    "\n",
    "# Step 2.a.v: compute S^i\n",
    "PoynSU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    PoynSU[i] = -alpha * (smallb2etk*uU[i]*u_0 + sp.Rational(1,2)*smallb2etk*g4UD[i+1][0] - smallb4U[i+1]*smallb4D[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation against `u0_smallb_Poynting__Cartesian` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for u0, smallbU, smallb2etk, and PoynSU between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [u0_smallb_Poynting__Cartesian module](../edit/u0_smallb_Poynting__Cartesian/u0_smallb_Poynting__Cartesian.py)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:18:34.466263Z",
     "iopub.status.busy": "2021-03-07T17:18:34.465184Z",
     "iopub.status.idle": "2021-03-07T17:18:34.488107Z",
     "shell.execute_reply": "2021-03-07T17:18:34.488639Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PASSED: u0 C code matches!\n",
      "u0etc.smallb4U[0] - smallb4U[0] = 0\n",
      "u0etc.smallb4U[1] - smallb4U[1] = 0\n",
      "u0etc.smallb4U[2] - smallb4U[2] = 0\n",
      "u0etc.smallb4U[3] - smallb4U[3] = 0\n",
      "u0etc.smallb2etk - smallb2etk = 0\n",
      "u0etc.PoynSU[0] - PoynSU[0] = 0\n",
      "u0etc.PoynSU[1] - PoynSU[1] = 0\n",
      "u0etc.PoynSU[2] - PoynSU[2] = 0\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import u0_smallb_Poynting__Cartesian.u0_smallb_Poynting__Cartesian as u0etc\n",
    "u0etc.compute_u0_smallb_Poynting__Cartesian(gammaDD,betaU,alpha,ValenciavU,BU)\n",
    "\n",
    "if u0etc.computeu0_Cfunction != computeu0_Cfunction:\n",
    "    print(\"FAILURE: u0 C code has changed!\")\n",
    "    sys.exit(1)\n",
    "else:\n",
    "    print(\"PASSED: u0 C code matches!\")\n",
    "\n",
    "for i in range(4):\n",
    "    print(\"u0etc.smallb4U[\"+str(i)+\"] - smallb4U[\"+str(i)+\"] = \"\n",
    "              + str(u0etc.smallb4U[i]-smallb4U[i]))\n",
    "\n",
    "print(\"u0etc.smallb2etk - smallb2etk = \" + str(u0etc.smallb2etk-smallb2etk))\n",
    "\n",
    "for i in range(DIM):\n",
    "    print(\"u0etc.PoynSU[\"+str(i)+\"] - PoynSU[\"+str(i)+\"] = \"\n",
    "              + str(u0etc.PoynSU[i]-PoynSU[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='appendix'></a>\n",
    "\n",
    "# Step 4: Appendix: Proving Eqs. 53 and 56 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf)\n",
    "$$\\label{appendix}$$\n",
    "\n",
    "$u^\\mu u_\\mu = -1$ implies\n",
    "\n",
    "\\begin{align}\n",
    "g^{\\mu\\nu} u_\\mu u_\\nu &= g^{00} \\left(u_0\\right)^2 + 2 g^{0i} u_0 u_i + g^{ij} u_i u_j = -1 \\\\\n",
    "\\implies &g^{00} \\left(u_0\\right)^2 + 2 g^{0i} u_0 u_i + g^{ij} u_i u_j + 1 = 0\\\\\n",
    "& a x^2 + b x + c = 0\n",
    "\\end{align}\n",
    "\n",
    "Thus we have a quadratic equation for $u_0$, with solution given by\n",
    "\n",
    "\\begin{align}\n",
    "u_0 &= \\frac{-b \\pm \\sqrt{b^2 - 4 a c}}{2 a} \\\\\n",
    "&= \\frac{-2 g^{0i}u_i \\pm \\sqrt{\\left(2 g^{0i} u_i\\right)^2 - 4 g^{00} (g^{ij} u_i u_j + 1)}}{2 g^{00}}\\\\\n",
    "&= \\frac{-g^{0i}u_i \\pm \\sqrt{\\left(g^{0i} u_i\\right)^2 - g^{00} (g^{ij} u_i u_j + 1)}}{g^{00}}\\\\\n",
    "\\end{align}\n",
    "\n",
    "Notice that (Eq. 4.49 in [Gourgoulhon](https://arxiv.org/pdf/gr-qc/0703035.pdf))\n",
    "$$\n",
    "g^{\\mu\\nu} = \\begin{pmatrix} \n",
    "-\\frac{1}{\\alpha^2} & \\frac{\\beta^i}{\\alpha^2} \\\\\n",
    "\\frac{\\beta^i}{\\alpha^2} & \\gamma^{ij} - \\frac{\\beta^i\\beta^j}{\\alpha^2}\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "so we have\n",
    "\n",
    "\\begin{align}\n",
    "u_0 &= \\frac{-\\beta^i u_i/\\alpha^2 \\pm \\sqrt{\\left(\\beta^i u_i/\\alpha^2\\right)^2 + 1/\\alpha^2 (g^{ij} u_i u_j + 1)}}{1/\\alpha^2}\\\\\n",
    "&= -\\beta^i u_i \\pm \\sqrt{\\left(\\beta^i u_i\\right)^2 + \\alpha^2 (g^{ij} u_i u_j + 1)}\\\\\n",
    "&= -\\beta^i u_i \\pm \\sqrt{\\left(\\beta^i u_i\\right)^2 + \\alpha^2 \\left(\\left[\\gamma^{ij} - \\frac{\\beta^i\\beta^j}{\\alpha^2}\\right] u_i u_j + 1\\right)}\\\\\n",
    "&= -\\beta^i u_i \\pm \\sqrt{\\left(\\beta^i u_i\\right)^2 + \\alpha^2 \\left(\\gamma^{ij}u_i u_j + 1\\right) - \\beta^i\\beta^j u_i u_j}\\\\\n",
    "&= -\\beta^i u_i \\pm \\sqrt{\\alpha^2 \\left(\\gamma^{ij}u_i u_j + 1\\right)}\\\\\n",
    "\\end{align}\n",
    "\n",
    "Now, since \n",
    "\n",
    "$$\n",
    "u^0 = g^{\\alpha 0} u_\\alpha = -\\frac{1}{\\alpha^2} u_0 + \\frac{\\beta^i u_i}{\\alpha^2},\n",
    "$$\n",
    "\n",
    "we get\n",
    "\n",
    "\\begin{align}\n",
    "u^0 &= \\frac{1}{\\alpha^2} \\left(u_0 + \\beta^i u_i\\right) \\\\\n",
    "&= \\pm \\frac{1}{\\alpha^2} \\sqrt{\\alpha^2 \\left(\\gamma^{ij}u_i u_j + 1\\right)}\\\\\n",
    "&= \\pm \\frac{1}{\\alpha} \\sqrt{\\gamma^{ij}u_i u_j + 1}\\\\\n",
    "\\end{align}\n",
    "\n",
    "By convention, the relativistic Gamma factor is positive and given by $\\alpha u^0$, so we choose the positive root. Thus we have derived Eq. 53 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf):\n",
    "\n",
    "$$\n",
    "u^0 = \\frac{1}{\\alpha} \\sqrt{\\gamma^{ij}u_i u_j + 1}.\n",
    "$$\n",
    "\n",
    "Next we evaluate \n",
    "\n",
    "\\begin{align}\n",
    "u^i &= u_\\mu g^{\\mu i} \\\\\n",
    "&= u_0 g^{0 i} + u_j g^{i j}\\\\\n",
    "&= u_0 \\frac{\\beta^i}{\\alpha^2} + u_j \\left(\\gamma^{ij} - \\frac{\\beta^i\\beta^j}{\\alpha^2}\\right)\\\\\n",
    "&= \\gamma^{ij} u_j + u_0 \\frac{\\beta^i}{\\alpha^2} - u_j \\frac{\\beta^i\\beta^j}{\\alpha^2}\\\\\n",
    "&= \\gamma^{ij} u_j + \\frac{\\beta^i}{\\alpha^2} \\left(u_0  - u_j \\beta^j\\right)\\\\\n",
    "&= \\gamma^{ij} u_j - \\beta^i u^0,\\\\\n",
    "\\implies v^i &= \\frac{\\gamma^{ij} u_j}{u^0} - \\beta^i\n",
    "\\end{align}\n",
    "\n",
    "which is equivalent to Eq. 56 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf). Notice in the last step, we used the above definition of $u^0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: 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-u0_smallb_Poynting-Cartesian.pdf](Tutorial-u0_smallb_Poynting-Cartesian.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:18:34.493584Z",
     "iopub.status.busy": "2021-03-07T17:18:34.492496Z",
     "iopub.status.idle": "2021-03-07T17:18:38.174534Z",
     "shell.execute_reply": "2021-03-07T17:18:38.173401Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[pandoc warning] Duplicate link reference `[comment]' \"source\" (line 22, column 1)\n",
      "Created Tutorial-u0_smallb_Poynting-Cartesian.tex, and compiled LaTeX file\n",
      "    to PDF file Tutorial-u0_smallb_Poynting-Cartesian.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-u0_smallb_Poynting-Cartesian\")"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}