{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\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:** Validated \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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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 }