{
"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",
"[comment]: <> (Abstract: TODO)\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",
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}