{
"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",
"# Equations of General Relativistic Hydrodynamics (GRHD)\n",
"\n",
"## Authors: Zach Etienne & Patrick Nelson\n",
"\n",
"## This notebook documents and constructs a number of quantities useful for building symbolic (SymPy) expressions for the equations of general relativistic hydrodynamics (GRHD), using the same (Valencia) formalism as `IllinoisGRMHD`.\n",
"\n",
"**Notebook Status:** <font color='orange'><b> Self-Validated </b></font>\n",
"\n",
"**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented. (TODO)**\n",
"\n",
"## Introduction\n",
"\n",
"We write the equations of general relativistic hydrodynamics in conservative form as follows (adapted from Eqs. 41-44 of [Duez et al](https://arxiv.org/pdf/astro-ph/0503420.pdf):\n",
"\n",
"\\begin{eqnarray}\n",
"\\ \\partial_t \\rho_* &+& \\partial_j \\left(\\rho_* v^j\\right) = 0 \\\\\n",
"\\partial_t \\tilde{\\tau} &+& \\partial_j \\left(\\alpha^2 \\sqrt{\\gamma} T^{0j} - \\rho_* v^j \\right) = s \\\\\n",
"\\partial_t \\tilde{S}_i &+& \\partial_j \\left(\\alpha \\sqrt{\\gamma} T^j{}_i \\right) = \\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i},\n",
"\\end{eqnarray}\n",
"where we assume $T^{\\mu\\nu}$ is the stress-energy tensor of a perfect fluid:\n",
"$$\n",
"T^{\\mu\\nu} = \\rho_0 h u^{\\mu} u^{\\nu} + P g^{\\mu\\nu},\n",
"$$\n",
"the $s$ source term is given in terms of ADM quantities via\n",
"\n",
"$$\n",
"s = \\alpha \\sqrt{\\gamma}\\left[\\left(T^{00}\\beta^i\\beta^j + 2 T^{0i}\\beta^j + T^{ij} \\right)K_{ij}\n",
"- \\left(T^{00}\\beta^i + T^{0i} \\right)\\partial_i\\alpha \\right],\n",
"$$\n",
"\n",
"and \n",
"\\begin{align}\n",
"v^j &= \\frac{u^j}{u^0} \\\\\n",
"\\rho_* &= \\alpha\\sqrt{\\gamma} \\rho_0 u^0 \\\\\n",
"h &= 1 + \\epsilon + \\frac{P}{\\rho_0}.\n",
"\\end{align}\n",
"\n",
"Also, we will write the 4-metric in terms of the ADM 3-metric, lapse, and shift using standard equations.\n",
"\n",
"Thus the full set of input variables includes:\n",
"* Spacetime quantities:\n",
" * ADM quantities $\\alpha$, $\\beta^i$, $\\gamma_{ij}$, $K_{ij}$\n",
"* Hydrodynamical quantities:\n",
" * Rest-mass density $\\rho_0$\n",
" * Pressure $P$\n",
" * Internal energy $\\epsilon$\n",
" * 4-velocity $u^\\mu$\n",
"\n",
"For completeness, the rest of the conservative variables are given by\n",
"\\begin{align}\n",
"\\tilde{\\tau} &= \\alpha^2\\sqrt{\\gamma} T^{00} - \\rho_* \\\\\n",
"\\tilde{S}_i &= \\alpha \\sqrt{\\gamma} T^0{}_i\n",
"\\end{align}\n",
"\n",
"### A Note on Notation\n",
"\n",
"As is standard in NRPy+, \n",
"\n",
"* Greek indices refer to four-dimensional quantities where the zeroth component indicates temporal (time) component.\n",
"* Latin indices refer to three-dimensional quantities. This is somewhat counterintuitive since Python always indexes its lists starting from 0. As a result, the zeroth component of three-dimensional quantities will necessarily indicate the first *spatial* direction.\n",
"\n",
"For instance, in calculating the first term of $b^2 u^\\mu u^\\nu$, we use Greek indices:\n",
"\n",
"```python\n",
"T4EMUU = ixp.zerorank2(DIM=4)\n",
"for mu in range(4):\n",
" for nu in range(4):\n",
" # Term 1: b^2 u^{\\mu} u^{\\nu}\n",
" T4EMUU[mu][nu] = smallb2*u4U[mu]*u4U[nu]\n",
"```\n",
"\n",
"When we calculate $\\beta_i = \\gamma_{ij} \\beta^j$, we use Latin indices:\n",
"```python\n",
"betaD = ixp.zerorank1(DIM=3)\n",
"for i in range(3):\n",
" for j in range(3):\n",
" betaD[i] += gammaDD[i][j] * betaU[j]\n",
"```\n",
"\n",
"As a corollary, any expressions involving mixed Greek and Latin indices will need to offset one set of indices by one: A Latin index in a four-vector will be incremented and a Greek index in a three-vector will be decremented (however, the latter case does not occur in this tutorial notebook). This can be seen when we handle $\\frac{1}{2} \\alpha \\sqrt{\\gamma} T^{\\mu \\nu}_{\\rm EM} \\partial_i g_{\\mu \\nu}$:\n",
"```python\n",
"# \\alpha \\sqrt{\\gamma} T^{\\mu \\nu}_{\\rm EM} \\partial_i g_{\\mu \\nu} / 2\n",
"for i in range(3):\n",
" for mu in range(4):\n",
" for nu in range(4):\n",
" S_tilde_rhsD[i] += alpsqrtgam * T4EMUU[mu][nu] * g4DD_zerotimederiv_dD[mu][nu][i+1] / 2\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"Each family of quantities is constructed within a given function (**boldfaced** below). This notebook is organized as follows\n",
"\n",
"\n",
"1. [Step 1](#importmodules): Import needed NRPy+ & Python modules\n",
"1. [Step 2](#stressenergy): Define the stress-energy tensor $T^{\\mu\\nu}$ and $T^\\mu{}_\\nu$:\n",
" * **compute_enthalpy()**, **compute_T4UU()**, **compute_T4UD()**: \n",
"1. [Step 3](#primtoconserv): Writing the conservative variables in terms of the primitive variables:\n",
" * **compute_sqrtgammaDET()**, **compute_rho_star()**, **compute_tau_tilde()**, **compute_S_tildeD()**\n",
"1. [Step 4](#grhdfluxes): Define the fluxes for the GRHD equations\n",
" 1. [Step 4.a](#rhostarfluxterm): Define $\\rho_*$ flux term for GRHD equations:\n",
" * **compute_vU_from_u4U__no_speed_limit()**, **compute_rho_star_fluxU()**:\n",
" 1. [Step 4.b](#taustildesourceterms) Define $\\tilde{\\tau}$ and $\\tilde{S}_i$ flux terms for GRHD equations:\n",
" * **compute_tau_tilde_fluxU()**, **compute_S_tilde_fluxUD()**\n",
"1. [Step 5](#grhdsourceterms): Define source terms on RHSs of GRHD equations\n",
" 1. [Step 5.a](#ssourceterm): Define $s$ source term on RHS of $\\tilde{\\tau}$ equation:\n",
" * **compute_s_source_term()**\n",
" 1. [Step 5.b](#stildeisourceterm): Define source term on RHS of $\\tilde{S}_i$ equation\n",
" 1. [Step 5.b.i](#fourmetricderivs): Compute $g_{\\mu\\nu,i}$ in terms of ADM quantities and their derivatives:\n",
" * **compute_g4DD_zerotimederiv_dD()**\n",
" 1. [Step 5.b.ii](#stildeisource): Compute source term of the $\\tilde{S}_i$ equation: $\\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i}$:\n",
" * **compute_S_tilde_source_termD()**\n",
"1. [Step 6](#convertvtou): Conversion of $v^i$ to $u^\\mu$ (Courtesy Patrick Nelson):\n",
" * **u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit()**, **u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit()**\n",
"1. [Step 7](#declarevarsconstructgrhdeqs): Declare ADM and hydrodynamical input variables, and construct GRHD equations\n",
"1. [Step 8](#code_validation): Code Validation against `GRHD.equations` NRPy+ module\n",
"1. [Step 9](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='importmodules'></a>\n",
"\n",
"# Step 1: Import needed NRPy+ & Python modules \\[Back to [top](#toc)\\]\n",
"$$\\label{importmodules}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:31.702854Z",
"iopub.status.busy": "2021-03-07T17:11:31.702078Z",
"iopub.status.idle": "2021-03-07T17:11:32.031485Z",
"shell.execute_reply": "2021-03-07T17:11:32.032077Z"
}
},
"outputs": [],
"source": [
"# Step 1: Import needed core NRPy+ modules\n",
"from outputC import nrpyAbs # NRPy+: Core C code output module\n",
"import NRPy_param_funcs as par # NRPy+: parameter interface\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='stressenergy'></a>\n",
"\n",
"# Step 2: Define the stress-energy tensor $T^{\\mu\\nu}$ and $T^\\mu{}_\\nu$ \\[Back to [top](#toc)\\]\n",
"$$\\label{stressenergy}$$\n",
"\n",
"Recall from above that\n",
"\n",
"$$\n",
"T^{\\mu\\nu} = \\rho_0 h u^{\\mu} u^{\\nu} + P g^{\\mu\\nu},\n",
"$$\n",
"where $h = 1 + \\epsilon + \\frac{P}{\\rho_0}$. Also \n",
"\n",
"$$\n",
"T^\\mu{}_{\\nu} = T^{\\mu\\delta} g_{\\delta \\nu}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.042090Z",
"iopub.status.busy": "2021-03-07T17:11:32.041286Z",
"iopub.status.idle": "2021-03-07T17:11:32.044195Z",
"shell.execute_reply": "2021-03-07T17:11:32.043650Z"
}
},
"outputs": [],
"source": [
"# Step 2.a: First define h, the enthalpy:\n",
"def compute_enthalpy(rho_b,P,epsilon):\n",
" global h\n",
" h = 1 + epsilon + P/rho_b\n",
"\n",
"# Step 2.b: Define T^{mu nu} (a 4-dimensional tensor)\n",
"def compute_T4UU(gammaDD,betaU,alpha, rho_b,P,epsilon,u4U):\n",
" global T4UU\n",
"\n",
" compute_enthalpy(rho_b,P,epsilon)\n",
" # Then define g^{mu nu} in terms of the ADM quantities:\n",
" import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
" AB4m.g4UU_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
"\n",
" # Finally compute T^{mu nu}\n",
" T4UU = ixp.zerorank2(DIM=4)\n",
" for mu in range(4):\n",
" for nu in range(4):\n",
" T4UU[mu][nu] = rho_b * h * u4U[mu]*u4U[nu] + P*AB4m.g4UU[mu][nu]\n",
"\n",
"# Step 2.c: Define T^{mu}_{nu} (a 4-dimensional tensor)\n",
"def compute_T4UD(gammaDD,betaU,alpha, T4UU):\n",
" global T4UD\n",
" # Next compute T^mu_nu = T^{mu delta} g_{delta nu}, needed for S_tilde flux.\n",
" # First we'll need g_{alpha nu} in terms of ADM quantities:\n",
" import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
" AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
" T4UD = ixp.zerorank2(DIM=4)\n",
" for mu in range(4):\n",
" for nu in range(4):\n",
" for delta in range(4):\n",
" T4UD[mu][nu] += T4UU[mu][delta]*AB4m.g4DD[delta][nu]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='primtoconserv'></a>\n",
"\n",
"# Step 3: Writing the conservative variables in terms of the primitive variables \\[Back to [top](#toc)\\]\n",
"$$\\label{primtoconserv}$$\n",
"\n",
"Recall from above that the conservative variables may be written as\n",
"\\begin{align}\n",
"\\rho_* &= \\alpha\\sqrt{\\gamma} \\rho_0 u^0 \\\\\n",
"\\tilde{\\tau} &= \\alpha^2\\sqrt{\\gamma} T^{00} - \\rho_* \\\\\n",
"\\tilde{S}_i &= \\alpha \\sqrt{\\gamma} T^0{}_i\n",
"\\end{align}\n",
"\n",
"$T^{\\mu\\nu}$ and $T^\\mu{}_\\nu$ have already been defined $-$ all in terms of primitive variables. Thus we'll just need $\\sqrt{\\gamma}=$`gammaDET`, and all conservatives can then be written in terms of other defined quantities, which themselves are written in terms of primitive variables and the ADM metric."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.051679Z",
"iopub.status.busy": "2021-03-07T17:11:32.050840Z",
"iopub.status.idle": "2021-03-07T17:11:32.054156Z",
"shell.execute_reply": "2021-03-07T17:11:32.053458Z"
}
},
"outputs": [],
"source": [
"# Step 3: Writing the conservative variables in terms of the primitive variables\n",
"def compute_sqrtgammaDET(gammaDD):\n",
" global sqrtgammaDET\n",
" _gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)\n",
" sqrtgammaDET = sp.sqrt(gammaDET)\n",
"\n",
"def compute_rho_star(alpha, sqrtgammaDET, rho_b,u4U):\n",
" global rho_star\n",
" # Compute rho_star:\n",
" rho_star = alpha*sqrtgammaDET*rho_b*u4U[0]\n",
"\n",
"def compute_tau_tilde(alpha, sqrtgammaDET, T4UU,rho_star):\n",
" global tau_tilde\n",
" tau_tilde = alpha**2*sqrtgammaDET*T4UU[0][0] - rho_star\n",
"\n",
"def compute_S_tildeD(alpha, sqrtgammaDET, T4UD):\n",
" global S_tildeD\n",
" S_tildeD = ixp.zerorank1(DIM=3)\n",
" for i in range(3):\n",
" S_tildeD[i] = alpha*sqrtgammaDET*T4UD[0][i+1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='grhdfluxes'></a>\n",
"\n",
"# Step 4: Define the fluxes for the GRHD equations \\[Back to [top](#toc)\\]\n",
"$$\\label{grhdfluxes}$$\n",
"\n",
"<a id='rhostarfluxterm'></a>\n",
"\n",
"## Step 4.a: Define $\\rho_*$ flux term for GRHD equations \\[Back to [top](#toc)\\]\n",
"$$\\label{rhostarfluxterm}$$\n",
"\n",
"Recall from above that\n",
"\\begin{array}\n",
"\\ \\partial_t \\rho_* &+ \\partial_j \\left(\\rho_* v^j\\right) = 0.\n",
"\\end{array}\n",
"\n",
"Here we will define the $\\rho_* v^j$ that constitutes the flux of $\\rho_*$, first defining $v^j=u^j/u^0$:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.060154Z",
"iopub.status.busy": "2021-03-07T17:11:32.059543Z",
"iopub.status.idle": "2021-03-07T17:11:32.061704Z",
"shell.execute_reply": "2021-03-07T17:11:32.062146Z"
}
},
"outputs": [],
"source": [
"# Step 4: Define the fluxes for the GRHD equations\n",
"# Step 4.a: vU from u4U may be needed for computing rho_star_flux from u4U\n",
"def compute_vU_from_u4U__no_speed_limit(u4U):\n",
" global vU\n",
" # Now compute v^i = u^i/u^0:\n",
" vU = ixp.zerorank1(DIM=3)\n",
" for j in range(3):\n",
" vU[j] = u4U[j+1]/u4U[0]\n",
"\n",
"# Step 4.b: rho_star flux\n",
"def compute_rho_star_fluxU(vU, rho_star):\n",
" global rho_star_fluxU\n",
" rho_star_fluxU = ixp.zerorank1(DIM=3)\n",
" for j in range(3):\n",
" rho_star_fluxU[j] = rho_star*vU[j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='taustildesourceterms'></a>\n",
"\n",
"## Step 4.b: Define $\\tilde{\\tau}$ and $\\tilde{S}_i$ flux terms for GRHD equations \\[Back to [top](#toc)\\]\n",
"$$\\label{taustildesourceterms}$$\n",
"\n",
"Recall from above that\n",
"\\begin{array}\n",
"\\ \\partial_t \\tilde{\\tau} &+ \\partial_j \\underbrace{\\left(\\alpha^2 \\sqrt{\\gamma} T^{0j} - \\rho_* v^j \\right)} &= s \\\\\n",
"\\partial_t \\tilde{S}_i &+ \\partial_j \\underbrace{\\left(\\alpha \\sqrt{\\gamma} T^j{}_i \\right)} &= \\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i}.\n",
"\\end{array}\n",
"\n",
"Here we will define all terms that go inside the $\\partial_j$'s on the left-hand side of the above equations (i.e., the underbraced expressions):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.069449Z",
"iopub.status.busy": "2021-03-07T17:11:32.068400Z",
"iopub.status.idle": "2021-03-07T17:11:32.071901Z",
"shell.execute_reply": "2021-03-07T17:11:32.071240Z"
}
},
"outputs": [],
"source": [
"# Step 4.c: tau_tilde flux\n",
"def compute_tau_tilde_fluxU(alpha, sqrtgammaDET, vU,T4UU, rho_star):\n",
" global tau_tilde_fluxU\n",
" tau_tilde_fluxU = ixp.zerorank1(DIM=3)\n",
" for j in range(3):\n",
" tau_tilde_fluxU[j] = alpha**2*sqrtgammaDET*T4UU[0][j+1] - rho_star*vU[j]\n",
"\n",
"# Step 4.d: S_tilde flux\n",
"def compute_S_tilde_fluxUD(alpha, sqrtgammaDET, T4UD):\n",
" global S_tilde_fluxUD\n",
" S_tilde_fluxUD = ixp.zerorank2(DIM=3)\n",
" for j in range(3):\n",
" for i in range(3):\n",
" S_tilde_fluxUD[j][i] = alpha*sqrtgammaDET*T4UD[j+1][i+1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='grhdsourceterms'></a>\n",
"\n",
"# Step 5: Define source terms on RHSs of GRHD equations \\[Back to [top](#toc)\\]\n",
"$$\\label{grhdsourceterms}$$\n",
"\n",
"<a id='ssourceterm'></a>\n",
"\n",
"## Step 5.a: Define $s$ source term on RHS of $\\tilde{\\tau}$ equation \\[Back to [top](#toc)\\]\n",
"$$\\label{ssourceterm}$$\n",
"\n",
"\n",
"Recall again from above the $s$ source term on the right-hand side of the $\\tilde{\\tau}$ evolution equation is given in terms of ADM quantities and the stress-energy tensor via\n",
"$$\n",
"s = \\underbrace{\\alpha \\sqrt{\\gamma}}_{\\text{Term 3}}\\left[\\underbrace{\\left(T^{00}\\beta^i\\beta^j + 2 T^{0i}\\beta^j + T^{ij} \\right)K_{ij}}_{\\text{Term 1}}\n",
"\\underbrace{- \\left(T^{00}\\beta^i + T^{0i} \\right)\\partial_i\\alpha}_{\\text{Term 2}} \\right],\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.080112Z",
"iopub.status.busy": "2021-03-07T17:11:32.079512Z",
"iopub.status.idle": "2021-03-07T17:11:32.081738Z",
"shell.execute_reply": "2021-03-07T17:11:32.082223Z"
}
},
"outputs": [],
"source": [
"def compute_s_source_term(KDD,betaU,alpha, sqrtgammaDET,alpha_dD, T4UU):\n",
" global s_source_term\n",
" s_source_term = sp.sympify(0)\n",
" # Term 1:\n",
" for i in range(3):\n",
" for j in range(3):\n",
" s_source_term += (T4UU[0][0]*betaU[i]*betaU[j] + 2*T4UU[0][i+1]*betaU[j] + T4UU[i+1][j+1])*KDD[i][j]\n",
" # Term 2:\n",
" for i in range(3):\n",
" s_source_term += -(T4UU[0][0]*betaU[i] + T4UU[0][i+1])*alpha_dD[i]\n",
" # Term 3:\n",
" s_source_term *= alpha*sqrtgammaDET"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='stildeisourceterm'></a>\n",
"\n",
"## Step 5.b: Define source term on RHS of $\\tilde{S}_i$ equation \\[Back to [top](#toc)\\]\n",
"$$\\label{stildeisourceterm}$$\n",
"\n",
"Recall from above\n",
"$$\n",
"\\partial_t \\tilde{S}_i + \\partial_j \\left(\\alpha \\sqrt{\\gamma} T^j{}_i \\right) = \\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i}.\n",
"$$\n",
"Our goal here will be to compute\n",
"$$\n",
"\\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i}.\n",
"$$\n",
"\n",
"<a id='fourmetricderivs'></a>\n",
"\n",
"### Step 5.b.i: Compute $g_{\\mu\\nu,i}$ in terms of ADM quantities and their derivatives \\[Back to [top](#toc)\\]\n",
"$$\\label{fourmetricderivs}$$\n",
"\n",
"\n",
"To compute $g_{\\mu\\nu,i}$ we need to evaluate the first derivative of $g_{\\mu\\nu}$ in terms of ADM variables.\n",
"\n",
"We are given $\\gamma_{ij}$, $\\alpha$, and $\\beta^i$, and the 4-metric is given in terms of these quantities via\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",
"$$\n",
"\n",
"Thus \n",
"$$\n",
"g_{\\mu\\nu,k} = \\begin{pmatrix} \n",
"-2 \\alpha\\alpha_{,i} + \\beta^j_{,k} \\beta_j + \\beta^j \\beta_{j,k} & \\beta_{i,k} \\\\\n",
"\\beta_{j,k} & \\gamma_{ij,k}\n",
"\\end{pmatrix},\n",
"$$\n",
"where $\\beta_{i} = \\gamma_{ij} \\beta^j$, so\n",
"$$\n",
"\\beta_{i,k} = \\gamma_{ij,k} \\beta^j + \\gamma_{ij} \\beta^j_{,k}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.093301Z",
"iopub.status.busy": "2021-03-07T17:11:32.092617Z",
"iopub.status.idle": "2021-03-07T17:11:32.094609Z",
"shell.execute_reply": "2021-03-07T17:11:32.095047Z"
}
},
"outputs": [],
"source": [
"def compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD):\n",
" global g4DD_zerotimederiv_dD\n",
" # Eq. 2.121 in B&S\n",
" betaD = ixp.zerorank1(DIM=3)\n",
" for i in range(3):\n",
" for j in range(3):\n",
" betaD[i] += gammaDD[i][j]*betaU[j]\n",
"\n",
" betaDdD = ixp.zerorank2(DIM=3)\n",
" for i in range(3):\n",
" for j in range(3):\n",
" for k in range(3):\n",
" # Recall that betaD[i] = gammaDD[i][j]*betaU[j] (Eq. 2.121 in B&S)\n",
" betaDdD[i][k] += gammaDD_dD[i][j][k]*betaU[j] + gammaDD[i][j]*betaU_dD[j][k]\n",
"\n",
" # Eq. 2.122 in B&S\n",
" g4DD_zerotimederiv_dD = ixp.zerorank3(DIM=4)\n",
" for k in range(3):\n",
" # Recall that g4DD[0][0] = -alpha^2 + betaU[j]*betaD[j]\n",
" g4DD_zerotimederiv_dD[0][0][k+1] += -2*alpha*alpha_dD[k]\n",
" for j in range(3):\n",
" g4DD_zerotimederiv_dD[0][0][k+1] += betaU_dD[j][k]*betaD[j] + betaU[j]*betaDdD[j][k]\n",
"\n",
" for i in range(3):\n",
" for k in range(3):\n",
" # Recall that g4DD[i][0] = g4DD[0][i] = betaD[i]\n",
" g4DD_zerotimederiv_dD[i+1][0][k+1] = g4DD_zerotimederiv_dD[0][i+1][k+1] = betaDdD[i][k]\n",
" for i in range(3):\n",
" for j in range(3):\n",
" for k in range(3):\n",
" # Recall that g4DD[i][j] = gammaDD[i][j]\n",
" g4DD_zerotimederiv_dD[i+1][j+1][k+1] = gammaDD_dD[i][j][k]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='stildeisource'></a>\n",
"\n",
"### Step 5.b.ii: Compute source term of the $\\tilde{S}_i$ equation: $\\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{stildeisource}$$\n",
"\n",
"Now that we've computed `g4DD_zerotimederiv_dD`$=g_{\\mu\\nu,i}$, the $\\tilde{S}_i$ evolution equation source term may be quickly constructed."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.101114Z",
"iopub.status.busy": "2021-03-07T17:11:32.100478Z",
"iopub.status.idle": "2021-03-07T17:11:32.102468Z",
"shell.execute_reply": "2021-03-07T17:11:32.102940Z"
}
},
"outputs": [],
"source": [
"# Step 5.b.ii: Compute S_tilde source term\n",
"def compute_S_tilde_source_termD(alpha, sqrtgammaDET,g4DD_zerotimederiv_dD, T4UU):\n",
" global S_tilde_source_termD\n",
" S_tilde_source_termD = ixp.zerorank1(DIM=3)\n",
" for i in range(3):\n",
" for mu in range(4):\n",
" for nu in range(4):\n",
" S_tilde_source_termD[i] += sp.Rational(1,2)*alpha*sqrtgammaDET*T4UU[mu][nu]*g4DD_zerotimederiv_dD[mu][nu][i+1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='convertvtou'></a>\n",
"\n",
"# Step 6: Conversion of $v^i$ to $u^\\mu$ (Courtesy Patrick Nelson) \\[Back to [top](#toc)\\]\n",
"$$\\label{convertvtou}$$\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 the largest allowed Lorentz factor as $\\Gamma_{\\rm max}$.\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)} \\to \\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",
"$$\n",
"and the remaining components $u^i$ via\n",
"$$\n",
"u^i = u^0 v^i.\n",
"$$\n",
"\n",
"In summary our algorithm for computing $u^{\\mu}$ from $v^i = \\frac{u^i}{u^0}$ is as follows:\n",
"\n",
"1. Choose a maximum Lorentz factor $\\Gamma_{\\rm max}$=`GAMMA_SPEED_LIMIT`, and define $v^i_{(n)} = \\frac{1}{\\alpha}\\left( \\frac{u^i}{u^0} + \\beta^i\\right)$.\n",
"1. Compute $R=\\gamma_{ij}v^i_{(n)}v^j_{(n)}=1 - \\frac{1}{\\Gamma^2}$\n",
"1. If $R \\le 1 - \\frac{1}{\\Gamma_{\\rm max}^2}$, then skip the next step.\n",
"1. Otherwise if $R > 1 - \\frac{1}{\\Gamma_{\\rm max}^2}$ then adjust $v^i_{(n)}\\to \\sqrt{\\frac{1 - \\frac{1}{\\Gamma_{\\rm max}^2}}{R}}v^i_{(n)}$, which will force $R=R_{\\rm max}$.\n",
"1. Given the $R$ computed in the above step, $u^0 = \\frac{1}{\\alpha \\sqrt{1-R}}$, and $u^i=u^0 v^i$.\n",
"\n",
"While the above algorithm is quite robust, its `if()` statement in the fourth step is not very friendly to NRPy+ or an optimizing C compiler, as it would require NRPy+ to generate separate C kernels for each branch of the `if()`. Let's instead try the following trick, which Roland Haas taught us. \n",
"\n",
"Define $R^*$ as\n",
"\n",
"$$\n",
"R^* = \\frac{1}{2} \\left(R_{\\rm max} + R - |R_{\\rm max} - R| \\right).\n",
"$$\n",
"\n",
"If $R>R_{\\rm max}$, then $|R_{\\rm max} - R|=R - R_{\\rm max}$, and we get:\n",
"\n",
"$$\n",
"R^* = \\frac{1}{2} \\left(R_{\\rm max} + R - (R - R_{\\rm max}) \\right) = \\frac{1}{2} \\left(2 R_{\\rm max}\\right) = R_{\\rm max}\n",
"$$\n",
"\n",
"If $R\\le R_{\\rm max}$, then $|R_{\\rm max} - R|=R_{\\rm max} - R$, and we get:\n",
"\n",
"$$\n",
"R^* = \\frac{1}{2} \\left(R_{\\rm max} + R - (R_{\\rm max} - R) \\right) = \\frac{1}{2} \\left(2 R\\right) = R\n",
"$$\n",
"\n",
"Then we can rescale *all* $v^i_{(n)}$ via\n",
"\n",
"$$\n",
"v^i_{(n)} \\to v^i_{(n)} \\sqrt{\\frac{R^*}{R}},\n",
"$$\n",
"\n",
"though we must be very careful to carefully handle the case in which $R=0$. To avoid any problems in this case, we simply adjust the above rescaling by adding a tiny number [`TINYDOUBLE`](https://en.wikipedia.org/wiki/Tiny_Bubbles) to $R$ in the denominator, typically `1e-100`:\n",
"\n",
"$$\n",
"v^i_{(n)} \\to v^i_{(n)} \\sqrt{\\frac{R^*}{R + {\\rm TINYDOUBLE}}}.\n",
"$$\n",
"\n",
"Finally, $u^0$ can be immediately and safely computed, via:\n",
"$$\n",
"u^0 = \\frac{1}{\\alpha \\sqrt{1-R^*}},\n",
"$$\n",
"and $u^i$ via \n",
"$$\n",
"u^i = u^0 v^i = u^0 \\left(\\alpha v^i_{(n)} - \\beta^i\\right).\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.115663Z",
"iopub.status.busy": "2021-03-07T17:11:32.114726Z",
"iopub.status.idle": "2021-03-07T17:11:32.116970Z",
"shell.execute_reply": "2021-03-07T17:11:32.117443Z"
}
},
"outputs": [],
"source": [
"# Step 6.a: Convert Valencia 3-velocity v_{(n)}^i into u^\\mu, and apply a speed limiter\n",
"# Speed-limited ValenciavU is output to rescaledValenciavU global.\n",
"def u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha,betaU,gammaDD, ValenciavU):\n",
" # Inputs: Metric lapse alpha, shift betaU, 3-metric gammaDD, Valencia 3-velocity ValenciavU\n",
" # Outputs (as globals): u4U_ito_ValenciavU, rescaledValenciavU\n",
"\n",
" # R = gamma_{ij} v^i v^j\n",
" R = sp.sympify(0)\n",
" for i in range(3):\n",
" for j in range(3):\n",
" R += gammaDD[i][j]*ValenciavU[i]*ValenciavU[j]\n",
"\n",
" thismodule = \"GRHD\"\n",
" # The default value isn't terribly important here, since we can overwrite in the main C code\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",
" # Now, we set Rstar = min(Rmax,R):\n",
" # If R < Rmax, then Rstar = 0.5*(Rmax+R-Rmax+R) = R\n",
" # If R >= Rmax, then Rstar = 0.5*(Rmax+R+Rmax-R) = Rmax\n",
" Rstar = sp.Rational(1, 2) * (Rmax + R - nrpyAbs(Rmax - R))\n",
"\n",
" # We add TINYDOUBLE to R below to avoid a 0/0, which occurs when\n",
" # ValenciavU == 0 for all Valencia 3-velocity components.\n",
" # \"Those tiny *doubles* make me warm all over\n",
" # with a feeling that I'm gonna love you till the end of time.\"\n",
" # - Adapted from Connie Francis' \"Tiny Bubbles\"\n",
" TINYDOUBLE = par.Cparameters(\"#define\",thismodule,\"TINYDOUBLE\",1e-100)\n",
"\n",
" # The rescaled (speed-limited) Valencia 3-velocity\n",
" # is given by, v_{(n)}^i = sqrt{Rstar/R} v^i\n",
" global rescaledValenciavU\n",
" rescaledValenciavU = ixp.zerorank1(DIM=3)\n",
" for i in range(3):\n",
" # If R == 0, then Rstar == 0, so sqrt( Rstar/(R+TINYDOUBLE) )=sqrt(0/1e-100) = 0\n",
" # If your velocities are of order 1e-100 and this is physically\n",
" # meaningful, there must be something wrong with your unit conversion.\n",
" rescaledValenciavU[i] = ValenciavU[i]*sp.sqrt(Rstar/(R + TINYDOUBLE))\n",
"\n",
" # Finally compute u^mu in terms of Valenciav^i\n",
" # u^0 = 1/(alpha-sqrt(1-R^*))\n",
" global u4U_ito_ValenciavU\n",
" u4U_ito_ValenciavU = ixp.zerorank1(DIM=4)\n",
" u4U_ito_ValenciavU[0] = 1/(alpha*sp.sqrt(1-Rstar))\n",
" # u^i = u^0 ( alpha v^i_{(n)} - beta^i ), where v^i_{(n)} is the Valencia 3-velocity\n",
" for i in range(3):\n",
" u4U_ito_ValenciavU[i+1] = u4U_ito_ValenciavU[0] * (alpha * rescaledValenciavU[i] - betaU[i])\n",
"\n",
"# Step 6.b: Convert v^i into u^\\mu, and apply a speed limiter.\n",
"# Speed-limited vU is output to rescaledvU global.\n",
"def u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit(alpha,betaU,gammaDD, vU):\n",
" ValenciavU = ixp.zerorank1(DIM=3)\n",
" for i in range(3):\n",
" ValenciavU[i] = (vU[i] + betaU[i])/alpha\n",
" u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha,betaU,gammaDD, ValenciavU)\n",
" # Since ValenciavU is written in terms of vU,\n",
" # u4U_ito_ValenciavU is actually u4U_ito_vU\n",
" global u4U_ito_vU\n",
" u4U_ito_vU = ixp.zerorank1(DIM=4)\n",
" for mu in range(4):\n",
" u4U_ito_vU[mu] = u4U_ito_ValenciavU[mu]\n",
" # Finally compute the rescaled (speed-limited) vU\n",
" global rescaledvU\n",
" rescaledvU = ixp.zerorank1(DIM=3)\n",
" for i in range(3):\n",
" rescaledvU[i] = alpha * rescaledValenciavU[i] - betaU[i]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='declarevarsconstructgrhdeqs'></a>\n",
"\n",
"# Step 7: Declare ADM and hydrodynamical input variables, and construct GRHD equations \\[Back to [top](#toc)\\]\n",
"$$\\label{declarevarsconstructgrhdeqs}$$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.129334Z",
"iopub.status.busy": "2021-03-07T17:11:32.128086Z",
"iopub.status.idle": "2021-03-07T17:11:32.446181Z",
"shell.execute_reply": "2021-03-07T17:11:32.447047Z"
}
},
"outputs": [],
"source": [
"# First define hydrodynamical quantities\n",
"u4U = ixp.declarerank1(\"u4U\", DIM=4)\n",
"rho_b,P,epsilon = sp.symbols('rho_b P epsilon',real=True)\n",
"\n",
"# Then ADM quantities\n",
"gammaDD = ixp.declarerank2(\"gammaDD\",\"sym01\",DIM=3)\n",
"KDD = ixp.declarerank2(\"KDD\" ,\"sym01\",DIM=3)\n",
"betaU = ixp.declarerank1(\"betaU\", DIM=3)\n",
"alpha = sp.symbols('alpha', real=True)\n",
"\n",
"# First compute stress-energy tensor T4UU and T4UD:\n",
"compute_T4UU(gammaDD,betaU,alpha, rho_b,P,epsilon,u4U)\n",
"compute_T4UD(gammaDD,betaU,alpha, T4UU)\n",
"\n",
"# Next sqrt(gamma)\n",
"compute_sqrtgammaDET(gammaDD)\n",
"\n",
"# Compute conservative variables in terms of primitive variables\n",
"compute_rho_star( alpha, sqrtgammaDET, rho_b,u4U)\n",
"compute_tau_tilde(alpha, sqrtgammaDET, T4UU,rho_star)\n",
"compute_S_tildeD( alpha, sqrtgammaDET, T4UD)\n",
"\n",
"# Then compute v^i from u^mu\n",
"compute_vU_from_u4U__no_speed_limit(u4U)\n",
"\n",
"# Next compute fluxes of conservative variables\n",
"compute_rho_star_fluxU( vU, rho_star)\n",
"compute_tau_tilde_fluxU(alpha, sqrtgammaDET, vU,T4UU, rho_star)\n",
"compute_S_tilde_fluxUD( alpha, sqrtgammaDET, T4UD)\n",
"\n",
"# Then declare derivatives & compute g4DD_zerotimederiv_dD\n",
"gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\",DIM=3)\n",
"betaU_dD = ixp.declarerank2(\"betaU_dD\" ,\"nosym\",DIM=3)\n",
"alpha_dD = ixp.declarerank1(\"alpha_dD\" ,DIM=3)\n",
"compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
"\n",
"# Then compute source terms on tau_tilde and S_tilde equations\n",
"compute_s_source_term(KDD,betaU,alpha, sqrtgammaDET,alpha_dD, T4UU)\n",
"compute_S_tilde_source_termD( alpha, sqrtgammaDET,g4DD_zerotimederiv_dD, T4UU)\n",
"\n",
"# Then compute the 4-velocities in terms of an input Valencia 3-velocity testValenciavU[i]\n",
"testValenciavU = ixp.declarerank1(\"testValenciavU\",DIM=3)\n",
"u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha,betaU,gammaDD, testValenciavU)\n",
"\n",
"# Finally compute the 4-velocities in terms of an input 3-velocity testvU[i] = u^i/u^0\n",
"testvU = ixp.declarerank1(\"testvU\",DIM=3)\n",
"u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit(alpha,betaU,gammaDD, testvU)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 8: Code Validation against `GRHD.equations` NRPy+ module \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"As a code validation check, we verify agreement in the SymPy expressions for the GRHD equations generated in\n",
"1. this tutorial versus\n",
"2. the NRPy+ [GRHD.equations](../edit/GRHD/equations.py) module."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.457940Z",
"iopub.status.busy": "2021-03-07T17:11:32.456810Z",
"iopub.status.idle": "2021-03-07T17:11:32.470829Z",
"shell.execute_reply": "2021-03-07T17:11:32.471356Z"
}
},
"outputs": [],
"source": [
"import GRHD.equations as Ge\n",
"\n",
"# First compute stress-energy tensor T4UU and T4UD:\n",
"Ge.compute_T4UU(gammaDD,betaU,alpha, rho_b,P,epsilon,u4U)\n",
"Ge.compute_T4UD(gammaDD,betaU,alpha, Ge.T4UU)\n",
"\n",
"# Next sqrt(gamma)\n",
"Ge.compute_sqrtgammaDET(gammaDD)\n",
"\n",
"# Compute conservative variables in terms of primitive variables\n",
"Ge.compute_rho_star( alpha, Ge.sqrtgammaDET, rho_b,u4U)\n",
"Ge.compute_tau_tilde(alpha, Ge.sqrtgammaDET, Ge.T4UU,Ge.rho_star)\n",
"Ge.compute_S_tildeD( alpha, Ge.sqrtgammaDET, Ge.T4UD)\n",
"\n",
"# Then compute v^i from u^mu\n",
"Ge.compute_vU_from_u4U__no_speed_limit(u4U)\n",
"\n",
"# Next compute fluxes of conservative variables\n",
"Ge.compute_rho_star_fluxU ( Ge.vU, Ge.rho_star)\n",
"Ge.compute_tau_tilde_fluxU(alpha, Ge.sqrtgammaDET, Ge.vU,Ge.T4UU, Ge.rho_star)\n",
"Ge.compute_S_tilde_fluxUD (alpha, Ge.sqrtgammaDET, Ge.T4UD)\n",
"\n",
"# Then declare derivatives & compute g4DD_zerotimederiv_dD\n",
"# gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\",DIM=3)\n",
"# betaU_dD = ixp.declarerank2(\"betaU_dD\" ,\"nosym\",DIM=3)\n",
"# alpha_dD = ixp.declarerank1(\"alpha_dD\" ,DIM=3)\n",
"Ge.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
"\n",
"# Finally compute source terms on tau_tilde and S_tilde equations\n",
"Ge.compute_s_source_term(KDD,betaU,alpha, Ge.sqrtgammaDET,alpha_dD, Ge.T4UU)\n",
"Ge.compute_S_tilde_source_termD( alpha, Ge.sqrtgammaDET,Ge.g4DD_zerotimederiv_dD,Ge.T4UU)\n",
"\n",
"GetestValenciavU = ixp.declarerank1(\"testValenciavU\")\n",
"Ge.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, GetestValenciavU)\n",
"\n",
"GetestvU = ixp.declarerank1(\"testvU\")\n",
"Ge.u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit( alpha, betaU, gammaDD, GetestvU)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.542274Z",
"iopub.status.busy": "2021-03-07T17:11:32.505403Z",
"iopub.status.idle": "2021-03-07T17:11:32.545869Z",
"shell.execute_reply": "2021-03-07T17:11:32.545317Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ALL TESTS PASSED!\n"
]
}
],
"source": [
"def comp_func(expr1,expr2,basename,prefixname2=\"Ge.\"):\n",
" if str(expr1-expr2)!=\"0\":\n",
" print(basename+\" - \"+prefixname2+basename+\" = \"+ str(expr1-expr2))\n",
" return 1\n",
" return 0\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",
"expr_list = []\n",
"exprcheck_list = []\n",
"namecheck_list = []\n",
"\n",
"namecheck_list.extend([\"sqrtgammaDET\",\"rho_star\",\"tau_tilde\",\"s_source_term\"])\n",
"exprcheck_list.extend([Ge.sqrtgammaDET,Ge.rho_star,Ge.tau_tilde,Ge.s_source_term])\n",
"expr_list.extend([sqrtgammaDET,rho_star,tau_tilde,s_source_term])\n",
"for mu in range(4):\n",
" namecheck_list.extend([gfnm(\"u4_ito_ValenciavU\",mu),gfnm(\"u4U_ito_vU\",mu)])\n",
" exprcheck_list.extend([ Ge.u4U_ito_ValenciavU[mu], Ge.u4U_ito_vU[mu]])\n",
" expr_list.extend( [ u4U_ito_ValenciavU[mu], u4U_ito_vU[mu]])\n",
" for nu in range(4):\n",
" namecheck_list.extend([gfnm(\"T4UU\",mu,nu),gfnm(\"T4UD\",mu,nu)])\n",
" exprcheck_list.extend([Ge.T4UU[mu][nu],Ge.T4UD[mu][nu]])\n",
" expr_list.extend([T4UU[mu][nu],T4UD[mu][nu]])\n",
" for delta in range(4):\n",
" namecheck_list.extend([gfnm(\"g4DD_zerotimederiv_dD\",mu,nu,delta)])\n",
" exprcheck_list.extend([Ge.g4DD_zerotimederiv_dD[mu][nu][delta]])\n",
" expr_list.extend([g4DD_zerotimederiv_dD[mu][nu][delta]])\n",
"\n",
"\n",
"for i in range(3):\n",
" namecheck_list.extend([gfnm(\"S_tildeD\",i),gfnm(\"vU\",i),gfnm(\"rho_star_fluxU\",i),\n",
" gfnm(\"tau_tilde_fluxU\",i),gfnm(\"S_tilde_source_termD\",i),\n",
" gfnm(\"rescaledValenciavU\",i), gfnm(\"rescaledvU\",i)])\n",
" exprcheck_list.extend([Ge.S_tildeD[i],Ge.vU[i],Ge.rho_star_fluxU[i],\n",
" Ge.tau_tilde_fluxU[i],Ge.S_tilde_source_termD[i],\n",
" Ge.rescaledValenciavU[i],Ge.rescaledvU[i]])\n",
" expr_list.extend([S_tildeD[i],vU[i],rho_star_fluxU[i],\n",
" tau_tilde_fluxU[i],S_tilde_source_termD[i],\n",
" rescaledValenciavU[i],rescaledvU[i]])\n",
" for j in range(3):\n",
" namecheck_list.extend([gfnm(\"S_tilde_fluxUD\",i,j)])\n",
" exprcheck_list.extend([Ge.S_tilde_fluxUD[i][j]])\n",
" expr_list.extend([S_tilde_fluxUD[i][j]])\n",
"\n",
"num_failures = 0\n",
"for i in range(len(expr_list)):\n",
" num_failures += comp_func(expr_list[i],exprcheck_list[i],namecheck_list[i])\n",
"\n",
"import sys\n",
"if num_failures == 0:\n",
" print(\"ALL TESTS PASSED!\")\n",
"else:\n",
" print(\"ERROR: \" + str(num_failures) + \" TESTS DID NOT PASS\")\n",
" sys.exit(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Step 9: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
"$$\\label{latex_pdf_output}$$\n",
"\n",
"The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
"[Tutorial-GRHD_Equations-Cartesian.pdf](Tutorial-GRHD_Equations-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": 13,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:32.551267Z",
"iopub.status.busy": "2021-03-07T17:11:32.550195Z",
"iopub.status.idle": "2021-03-07T17:11:36.834756Z",
"shell.execute_reply": "2021-03-07T17:11:36.835554Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-GRHD_Equations-Cartesian.tex, and compiled LaTeX file to\n",
" PDF file Tutorial-GRHD_Equations-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-GRHD_Equations-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
}