{
"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",
"# Tutorial-IllinoisGRMHD: harm_u2p_util.c\n",
"\n",
"## Authors: Leo Werneck & Zach Etienne\n",
"\n",
"<font color='red'>**This module is currently under development**</font>\n",
"\n",
"## In this tutorial module we explain the utility functions needed by the conservative-to-primitive algorithm used by `HARM`. This module will likely be absorbed by another one once we finish documenting the code.\n",
"\n",
"### Required and recommended citations:\n",
"\n",
"* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
"* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
"* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This module is organized as follows\n",
"\n",
"0. [Step 0](#src_dir): **Source directory creation**\n",
"1. [Step 1](#introduction): **Introduction**\n",
"1. [Step 2](#harm_utoprim_2d__c__eos_indep): **EOS independent routines**\n",
" 1. [Step 2.a](#raise_g): *The `raise_g()` function*\n",
" 1. [Step 2.b](#lower_g): *The `lower_g()` function*\n",
" 1. [Step 2.c](#ncov_calc): *The `ncov_calc()` function*\n",
"1. [Step 3](#harm_utoprim_2d__c__eos_dep): **EOS dependent routines**\n",
" 1. [Step 3.a](#pressure_rho0_u): *The `pressure_rho0_u()` function*\n",
" 1. [Step 3.b](#pressure_rho0_w): *The `pressure_rho0_w()` function*\n",
"1. [Step n-1](#code_validation): **Code validation**\n",
"1. [Step n](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='src_dir'></a>\n",
"\n",
"# Step 0: Source directory creation \\[Back to [top](#toc)\\]\n",
"$$\\label{src_dir}$$\n",
"\n",
"We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory, if it does not exist yet."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Step 0: Creation of the IllinoisGRMHD source directory\n",
"# Step 0a: Add NRPy's directory to the path\n",
"# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
"import os,sys\n",
"nrpy_dir_path = os.path.join(\"..\",\"..\")\n",
"if nrpy_dir_path not in sys.path:\n",
" sys.path.append(nrpy_dir_path)\n",
"\n",
"# Step 0b: Load up cmdline_helper and create the directory\n",
"import cmdline_helper as cmd\n",
"IGM_src_dir_path = os.path.join(\"..\",\"src\")\n",
"cmd.mkdir(IGM_src_dir_path)\n",
"\n",
"# Step 0c: Create the output file path\n",
"outfile_path__harm_u2p_util__c = os.path.join(IGM_src_dir_path,\"harm_u2p_util.c\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='introduction'></a>\n",
"\n",
"# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
"$$\\label{introduction}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='harm_utoprim_2d__c__eos_indep'></a>\n",
"\n",
"# Step 2: EOS independent routines` \\[Back to [top](#toc)\\]\n",
"$$\\label{harm_utoprim_2d__c__eos_indep}$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing ../src/harm_u2p_util.c\n"
]
}
],
"source": [
"%%writefile $outfile_path__harm_u2p_util__c\n",
"#ifndef __HARM_U2P_UTIL__C__\n",
"#define __HARM_U2P_UTIL__C__\n",
"/*\n",
" -------------------------------------------------------------------------------\n",
" Copyright 2005 Scott C. Noble, Charles F. Gammie,\n",
" Jonathan C. McKinney, and Luca Del Zanna\n",
"\n",
"\n",
" This file is part of PVS-GRMHD.\n",
"\n",
" PVS-GRMHD is free software; you can redistribute it and/or modify\n",
" it under the terms of the GNU General Public License as published by\n",
" the Free Software Foundation; either version 2 of the License, or\n",
" (at your option) any later version.\n",
"\n",
" PVS-GRMHD is distributed in the hope that it will be useful,\n",
" but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
" MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n",
" GNU General Public License for more details.\n",
"\n",
" You should have received a copy of the GNU General Public License\n",
" along with PVS-GRMHD; if not, write to the Free Software\n",
" Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA\n",
"\n",
" -------------------------------------------------------------------------------\n",
"*/\n",
"\n",
"// Function prototypes for this file:\n",
"static void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM]);\n",
"static void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM]);\n",
"static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]);\n",
"static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u);\n",
"static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w);\n",
"\n",
"// Inlined function used by this file\n",
"static inline void compute_P_cold__eps_cold(eos_struct eos,CCTK_REAL rho_in, CCTK_REAL &P_cold,CCTK_REAL &eps_cold);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='raise_g'></a>\n",
"\n",
"## Step 2.a: The `raise_g()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{raise_g}$$\n",
"\n",
"This is a simple function, used to *raise* the indices of a *covariant vector* using an inverse metric, $g^{\\mu\\nu}$. Usually the vector is 4-dimensional and $g^{\\mu\\nu}$ is the inverse physical ADM 4-metric, but the function can be used with arbitrary vectors and metrics. In other words, given a vector $v_{\\mu}$ and an inverse metric $g^{\\mu\\nu}$ the function outputs\n",
"\n",
"$$\n",
"\\boxed{v^{\\mu} = g^{\\mu\\nu}v_{\\nu}}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/harm_u2p_util.c\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__harm_u2p_util__c\n",
"\n",
"\n",
"/**********************************************************************\n",
" raise_g():\n",
"\n",
" -- calculates the contravariant form of a covariant tensor,\n",
" using the inverse of the metric;\n",
"***********************************************************************/\n",
"static void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM])\n",
"{\n",
" int i,j;\n",
"\n",
" for(i=0;i<NDIM;i++) {\n",
" vcon[i] = 0. ;\n",
" for(j=0;j<NDIM;j++)\n",
" vcon[i] += gcon[i][j]*vcov[j] ;\n",
" }\n",
"\n",
" return ;\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='lower_g'></a>\n",
"\n",
"## Step 2.b: The `lower_g()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{lower_g}$$\n",
"\n",
"This is a simple function, used to *lower* the indices of a *contravariant vector* using a metric, $g_{\\mu\\nu}$. Usually the vector is 4-dimensional and $g_{\\mu\\nu}$ is the physical ADM 4-metric, but the function can be used with arbitrary vectors and metrics. In other words, given a vector $v^{\\mu}$ and a metric $g_{\\mu\\nu}$ the function outputs\n",
"\n",
"$$\n",
"\\boxed{v_{\\mu} = g_{\\mu\\nu}v^{\\nu}}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/harm_u2p_util.c\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__harm_u2p_util__c\n",
"\n",
"\n",
"/**********************************************************************\n",
" lower_g():\n",
"\n",
" -- calculates the ocvariant form of a contravariant tensor\n",
" using the metric;\n",
"***********************************************************************/\n",
"static void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM])\n",
"{\n",
" int i,j;\n",
"\n",
" for(i=0;i<NDIM;i++) {\n",
" vcov[i] = 0. ;\n",
" for(j=0;j<NDIM;j++)\n",
" vcov[i] += gcov[i][j]*vcon[j] ;\n",
" }\n",
"\n",
" return ;\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='ncov_calc'></a>\n",
"\n",
"## Step 2.c: The `ncov_calc()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{ncov_calc}$$\n",
"\n",
"This simple function sets the covariant normal vector $n_{\\mu} = \\left(-\\alpha,0,0,0\\right)$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/harm_u2p_util.c\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__harm_u2p_util__c\n",
"\n",
"\n",
"/**********************************************************************\n",
" ncov_calc():\n",
"\n",
" -- calculates the covariant form of the normal vector to our\n",
" spacelike hypersurfaces ala the ADM formalism.\n",
"\n",
" -- requires the inverse metric;\n",
"***********************************************************************/\n",
"static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM])\n",
"{\n",
" CCTK_REAL lapse ;\n",
" int i;\n",
"\n",
" lapse = sqrt(-1./gcon[0][0]) ;\n",
"\n",
" ncov[0] = -lapse ;\n",
" for( i = 1; i < NDIM; i++) {\n",
" ncov[i] = 0. ;\n",
" }\n",
"\n",
" return ;\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='harm_utoprim_2d__c__eos_dep'></a>\n",
"\n",
"# Step 3: EOS dependent routines\n",
"$$\\label{harm_utoprim_2d__c__eos_dep}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='pressure_rho0_u'></a>\n",
"\n",
"## Step 3.a: The `pressure_rho0_u()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{pressure_rho0_u}$$\n",
"\n",
"The $\\Gamma$-law EOS implemented in `HARM` is\n",
"\n",
"$$\n",
"p_{\\Gamma}\\left(\\rho_{b},u\\right) = \\left(\\Gamma-1\\right)u\\ ,\n",
"$$\n",
"\n",
"where\n",
"\n",
"$$\n",
"u = \\rho_{b}\\epsilon\\ .\n",
"$$\n",
"\n",
"Thus, the pre-PPEOS Patch version of this function was\n",
"\n",
"```c\n",
"/**************************************************\n",
" The following functions assume a Gamma-law EOS:\n",
"***************************************************/\n",
"\n",
"/* \n",
" pressure as a function of rho0 and u \n",
" this is used by primtoU and Utoprim_?D\n",
"*/\n",
"static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u)\n",
"{\n",
" return((GAMMA - 1.)*u) ;\n",
"}\n",
"```\n",
"\n",
"In the case of a hybrid EOS, however, we have $p_{\\rm hybrid}=p_{\\rm hybrid}\\left(\\rho_{b},\\epsilon\\right)$. To obtain $p_{\\rm hybrid}\\left(\\rho_{b},u\\right)$ we use:\n",
"\n",
"$$\n",
"p_{\\rm hybrid} = P_{\\rm cold}\\left(\\rho_{b}\\right) + \\left(\\Gamma_{\\rm th}-1\\right)\\rho_{b}\\left[\\epsilon - \\epsilon_{\\rm cold}\\left(\\rho_{b}\\right)\\right]\n",
"\\implies\n",
"\\boxed{\n",
"p_{\\rm hybrid}\\left(\\rho_{b},u\\right) = P_{\\rm cold}\\left(\\rho_{b}\\right) + \\left(\\Gamma_{\\rm th}-1\\right)\\left[u - \\rho_{b}\\epsilon_{\\rm cold}\\left(\\rho_{b}\\right)\\right]\\ ,\n",
"}\n",
"$$\n",
"\n",
"where\n",
"\n",
"$$\n",
"\\left\\{\n",
"\\begin{align}\n",
"P_{\\rm cold}\\left(\\rho_{b}\\right) = K_{\\rm poly}\\rho_{b}^{\\Gamma_{\\rm poly}}\\ ,\\\\\n",
"\\epsilon_{\\rm cold}\\left(\\rho_{b}\\right) = C + \\frac{P_{\\rm cold}}{\\rho_{b}\\left(\\Gamma_{\\rm poly}-1\\right)}\\ ,\n",
"\\end{align}\n",
"\\right.\n",
"$$\n",
"\n",
"where $C$ is a precomputed integration constant which guarantees the continuity of $\\epsilon_{\\rm cold}\\left(\\rho_{b}\\right)$ and the subscript \"poly\" indicates the local polytropic EOS quantity."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/harm_u2p_util.c\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__harm_u2p_util__c\n",
"\n",
"/**************************************************\n",
" The following functions assume a Gamma-law EOS:\n",
"***************************************************/\n",
"\n",
"/*\n",
"pressure as a function of rho0 and u\n",
"this is used by primtoU and Utoprim_?D\n",
"*/\n",
"static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u)\n",
"{\n",
"\n",
" // Set up Gamma_th:\n",
"#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
" DECLARE_CCTK_PARAMETERS;\n",
"#endif\n",
"\n",
" // Compute P_cold, eps_cold\n",
" CCTK_REAL P_cold, eps_cold;\n",
" compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);\n",
"\n",
" /* Compute the pressure as a function of rho_b (rho0) and\n",
" * u = rho_b * eps, using our hybrid EOS:\n",
" * .-------------------------------------------------------------.\n",
" * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |\n",
" * .-------------------------------------------------------------.\n",
" */\n",
" return( P_cold + (Gamma_th - 1.0)*(u - rho0*eps_cold) );\n",
"\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='pressure_rho0_w'></a>\n",
"\n",
"## Step 3.2: The `pressure_rho0_w()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{pressure_rho0_w}$$\n",
"\n",
"The $\\Gamma$-law EOS implemented in `HARM` is\n",
"\n",
"$$\n",
"p_{\\Gamma} = \\left(\\Gamma-1\\right)u\\ .\n",
"$$\n",
"\n",
"We want now to obtain $p_{\\Gamma}\\left(\\rho_{b},w\\right)$, where\n",
"\n",
"$$\n",
"w = u + \\rho_{b} + p\\ .\n",
"$$\n",
"\n",
"Then\n",
"\n",
"$$\n",
"p_{\\Gamma} = \\left(\\Gamma-1\\right)\\left(w - \\rho_{b} - p_{\\Gamma}\\right)\n",
"\\implies\n",
"\\boxed{\n",
"p_{\\Gamma}\\left(\\rho_{b},w\\right) = \\frac{\\left(\\Gamma-1\\right)}{\\Gamma}\\left(w-\\rho_{b}\\right)\n",
"}\\ .\n",
"$$\n",
"\n",
"Thus, the pre-PPEOS Patch version of this function was\n",
"\n",
"```c\n",
"/* \n",
" pressure as a function of rho0 and w = rho0 + u + p \n",
" this is used by primtoU and Utoprim_1D\n",
"*/\n",
"static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w)\n",
"{\n",
" return((GAMMA-1.)*(w - rho0)/GAMMA) ;\n",
"}\n",
"```\n",
"\n",
"For our hybrid EOS, we have\n",
"\n",
"$$\n",
"\\begin{align}\n",
"p_{\\rm hybrid} &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\rho_{b}\\left[\\epsilon - \\epsilon_{\\rm cold}\\right]\\\\\n",
"&= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left[u - \\rho_{b}\\epsilon_{\\rm cold}\\right]\\\\\n",
"&= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left[w-\\rho_{b}-p_{\\rm hybrid} - \\rho_{b}\\epsilon_{\\rm cold}\\right]\\\\\n",
"&= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left[w - \\rho_{b}\\left(1+\\epsilon_{\\rm cold}\\right)\\right]- \\left(\\Gamma_{\\rm th}-1\\right)p_{\\rm hybrid}\\\\\n",
"\\implies\n",
"&\\boxed{\n",
"p_{\\rm hybrid}\\left(\\rho_{b},w\\right) = \\frac{P_{\\rm cold}}{\\Gamma_{\\rm th}} + \\frac{\\left(\\Gamma_{\\rm th}-1\\right)}{\\Gamma_{\\rm th}}\\left[w - \\rho_{b}\\left(1+\\epsilon_{\\rm cold}\\right)\\right]\n",
"}\n",
"\\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/harm_u2p_util.c\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__harm_u2p_util__c\n",
"\n",
"\n",
"/*\n",
" pressure as a function of rho0 and w = rho0 + u + p\n",
" this is used by primtoU and Utoprim_1D\n",
"*/\n",
"static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w)\n",
"{\n",
"\n",
" // Set up Gamma_th:\n",
"#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
" DECLARE_CCTK_PARAMETERS;\n",
"#endif\n",
"\n",
" // Compute P_cold, eps_cold\n",
" CCTK_REAL P_cold, eps_cold;\n",
" compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);\n",
"\n",
" /* Compute the pressure as a function of rho_b (rho0) and\n",
" * w = u + rho_b + p, using our hybrid EOS:\n",
" * ----------------------------------------------------------------------------\n",
" * | p(rho_b,w) = ( P_cold + (Gamma_th-1)*( w - rho_b*(1+eps_cold) ) )/Gamma_th |\n",
" * ----------------------------------------------------------------------------\n",
" */\n",
" return( (P_cold + (Gamma_th-1.0)*( w - rho0*(1.0+eps_cold) ) )/Gamma_th );\n",
"}\n",
"#endif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 4: Code validation \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Validation test for harm_u2p_util.c: FAILED!\n",
"Diff:\n",
"1c1,2\n",
"< \n",
"---\n",
"> #ifndef __HARM_U2P_UTIL__C__\n",
"> #define __HARM_U2P_UTIL__C__\n",
"4c5\n",
"< Copyright 2005 Scott C. Noble, Charles F. Gammie, \n",
"---\n",
"> Copyright 2005 Scott C. Noble, Charles F. Gammie,\n",
"28,144c29,33\n",
"< void primtoU_g( CCTK_REAL prim[], CCTK_REAL gcov[][NDIM], CCTK_REAL gcon[][NDIM], CCTK_REAL gdet, CCTK_REAL U[] );\n",
"< static void ucon_calc_g(CCTK_REAL prim[],CCTK_REAL gcov[][NDIM],CCTK_REAL gcon[][NDIM],CCTK_REAL ucon[]);\n",
"< static void raise_g(CCTK_REAL vcov[], CCTK_REAL gcon[][NDIM], CCTK_REAL vcon[]);\n",
"< static void lower_g(CCTK_REAL vcon[], CCTK_REAL gcov[][NDIM], CCTK_REAL vcov[]);\n",
"< static void ncov_calc(CCTK_REAL gcon[][NDIM],CCTK_REAL ncov[]) ;\n",
"< static void bcon_calc_g(CCTK_REAL prim[],CCTK_REAL ucon[],CCTK_REAL ucov[],CCTK_REAL ncov[],CCTK_REAL bcon[]); \n",
"< static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u);\n",
"< static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w);\n",
"< int gamma_calc_g(CCTK_REAL *pr, CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL *gamma);\n",
"< \n",
"< /********************************************************************** \n",
"< primtoU_g(): \n",
"< \n",
"< -- calculates the conserved variables from the primitive variables \n",
"< and the metric;\n",
"< -- assumes that the conserved and primitive variables are defined ala HARM:\n",
"< \n",
"< / rho u^t \\\n",
"< U = | T^t_\\mu + rho u^t | sqrt(-det(g_{\\mu\\nu}))\n",
"< \\ B^i /\n",
"< \n",
"< / rho \\\n",
"< P = | uu |\n",
"< | \\tilde{u}^i |\n",
"< \\ B^i /\n",
"< \n",
"< **************************************************************************/\n",
"< \n",
"< \n",
"< void primtoU_g(\n",
"< CCTK_REAL prim[NPR], /* primitive variables */\n",
"< CCTK_REAL gcov[NDIM][NDIM], /* covariant (index dn) form of metric */\n",
"< CCTK_REAL gcon[NDIM][NDIM], /* contravariant (index up) form of metric */\n",
"< CCTK_REAL gdet, /* sqrt of -1 times det(g_{\\mu \\nu}) */\n",
"< CCTK_REAL U[NPR] /* matrix of derivatives */\n",
"< ) {\n",
"< int i ;\n",
"< CCTK_REAL rho0 ;\n",
"< static CCTK_REAL ucon[NDIM],ucov[NDIM],bcon[NDIM],bcov[NDIM],ncov[NDIM] ;\n",
"< CCTK_REAL gamma,n_dot_b,bsq,u,p,w, alpha ;\n",
"< \n",
"< \n",
"< /* Calculate auxiliary quantities: */\n",
"< alpha = 1.0/sqrt(-gcon[0][0]);\n",
"< \n",
"< ucon_calc_g(prim,gcov,gcon,ucon) ;\n",
"< lower_g(ucon,gcov,ucov) ;\n",
"< ncov_calc(gcon,ncov) ;\n",
"< \n",
"< gamma = -ncov[0]*ucon[0] ;\n",
"< \n",
"< bcon_calc_g(prim,ucon,ucov,ncov,bcon) ;\n",
"< lower_g(bcon,gcov,bcov) ;\n",
"< \n",
"< n_dot_b = 0. ;\n",
"< for(i=0;i<4;i++) n_dot_b += ncov[i]*bcon[i] ;\n",
"< bsq = 0. ;\n",
"< for(i=0;i<4;i++) bsq += bcov[i]*bcon[i] ;\n",
"< \n",
"< rho0 = prim[RHO] ;\n",
"< u = prim[UU] ;\n",
"< p = pressure_rho0_u(rho0,u) ;\n",
"< w = rho0 + u + p ;\n",
"< \n",
"< // Now set the conserved variables themselves, using HARM's definition:\n",
"< U[RHO] = ucon[0]*rho0 ;\n",
"< \n",
"< for( i = 0; i < 4; i++) {\n",
"< U[QCOV0+i] = gamma*(w + bsq)*ucov[i] \n",
"< - (p + bsq/2.)*ncov[i] \n",
"< + n_dot_b*bcov[i] ;\n",
"< \n",
"< U[QCOV0+i] /= alpha;\n",
"< }\n",
"< \n",
"< U[QCOV0] = U[QCOV0] + U[RHO];\n",
"< U[BCON1] = prim[BCON1] ;\n",
"< U[BCON2] = prim[BCON2] ;\n",
"< U[BCON3] = prim[BCON3] ;\n",
"< \n",
"< for(i = 0; i < NPR; i++ ) {\n",
"< U[i] *= gdet;\n",
"< }\n",
"< \n",
"< return ;\n",
"< }\n",
"< \n",
"< /********************************************************************** \n",
"< ucon_calc_g(): \n",
"< \n",
"< -- calculates the contravariant (up) components of the four-velocity\n",
"< given the primitive variables, of which the velocity is \n",
"< \\tilde{u}^i = \\gamma v^j where v^j is the velocity of the \n",
"< flow w.r.t a normal observer to the coordinates;\n",
"< \n",
"< -- also requires the metric and inverse metric;\n",
"< \n",
"< -- assumes:\n",
"< \n",
"< / rho \\\n",
"< P = | uu |\n",
"< | \\tilde{u}^i |\n",
"< \\ B^i /\n",
"< \n",
"< ******************************************************************/\n",
"< static void ucon_calc_g(CCTK_REAL prim[NPR],CCTK_REAL gcov[NDIM][NDIM],CCTK_REAL gcon[NDIM][NDIM],\n",
"< CCTK_REAL ucon[NDIM])\n",
"< {\n",
"< CCTK_REAL u_tilde_con[4] ;\n",
"< CCTK_REAL u_tilde_sq ;\n",
"< CCTK_REAL gamma,lapse ;\n",
"< int i,j;\n",
"< \n",
"< u_tilde_con[0] = 0. ;\n",
"< u_tilde_con[1] = prim[UTCON1] ;\n",
"< u_tilde_con[2] = prim[UTCON2] ;\n",
"< u_tilde_con[3] = prim[UTCON3] ;\n",
"---\n",
"> static void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM]);\n",
"> static void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM]);\n",
"> static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]);\n",
"> static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u);\n",
"> static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w);\n",
"146,156c35,36\n",
"< u_tilde_sq = 0. ;\n",
"< for(i=0;i<NDIM;i++)\n",
"< for(j=0;j<NDIM;j++)\n",
"< u_tilde_sq += gcov[i][j]*u_tilde_con[i]*u_tilde_con[j] ;\n",
"< u_tilde_sq = fabs(u_tilde_sq) ;\n",
"< \n",
"< gamma = sqrt(1. + u_tilde_sq) ;\n",
"< \n",
"< lapse = sqrt(-1./gcon[0][0]) ;\n",
"< \n",
"< for(i=0;i<NDIM;i++) ucon[i] = u_tilde_con[i] - lapse*gamma*gcon[0][i] ;\n",
"---\n",
"> // Inlined function used by this file\n",
"> static inline void compute_P_cold__eps_cold(eos_struct eos,CCTK_REAL rho_in, CCTK_REAL &P_cold,CCTK_REAL &eps_cold);\n",
"158,159d37\n",
"< return ;\n",
"< }\n",
"161c39\n",
"< /********************************************************************** \n",
"---\n",
"> /**********************************************************************\n",
"163,164c41,42\n",
"< \n",
"< -- calculates the contravariant form of a covariant tensor, \n",
"---\n",
"> \n",
"> -- calculates the contravariant form of a covariant tensor,\n",
"166c44\n",
"< ******************************************************************/\n",
"---\n",
"> ***********************************************************************/\n",
"173c51\n",
"< for(j=0;j<NDIM;j++) \n",
"---\n",
"> for(j=0;j<NDIM;j++)\n",
"180c58,59\n",
"< /********************************************************************** \n",
"---\n",
"> \n",
"> /**********************************************************************\n",
"182,183c61,62\n",
"< \n",
"< -- calculates the ocvariant form of a contravariant tensor \n",
"---\n",
"> \n",
"> -- calculates the ocvariant form of a contravariant tensor\n",
"185c64\n",
"< ******************************************************************/\n",
"---\n",
"> ***********************************************************************/\n",
"192c71\n",
"< for(j=0;j<NDIM;j++) \n",
"---\n",
"> for(j=0;j<NDIM;j++)\n",
"199,200d77\n",
"< /********************************************************************** \n",
"< ncov_calc(): \n",
"202c79,82\n",
"< -- calculates the covariant form of the normal vector to our \n",
"---\n",
"> /**********************************************************************\n",
"> ncov_calc():\n",
"> \n",
"> -- calculates the covariant form of the normal vector to our\n",
"206,207c86,87\n",
"< ******************************************************************/\n",
"< static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]) \n",
"---\n",
"> ***********************************************************************/\n",
"> static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM])\n",
"215c95\n",
"< for( i = 1; i < NDIM; i++) { \n",
"---\n",
"> for( i = 1; i < NDIM; i++) {\n",
"222,292d101\n",
"< /********************************************************************** \n",
"< bcon_calc_g(): \n",
"< \n",
"< -- using the primitive variables, contra-/co-variant 4-vel., \n",
"< and covariant normal vector, calculate the contravariant \n",
"< form of the magnetic 4-vector b^\\mu (the small \"b\" in HARM);\n",
"< -- assumes:\n",
"< \n",
"< / rho \\\n",
"< P = | uu |\n",
"< | \\tilde{u}^i |\n",
"< \\ B^i /\n",
"< ******************************************************************/\n",
"< static void bcon_calc_g(CCTK_REAL prim[NPR],CCTK_REAL ucon[NDIM],CCTK_REAL ucov[NDIM],\n",
"< CCTK_REAL ncov[NDIM],CCTK_REAL bcon[NDIM]) \n",
"< {\n",
"< static CCTK_REAL Bcon[NDIM] ;\n",
"< CCTK_REAL u_dot_B ;\n",
"< CCTK_REAL gamma ;\n",
"< int i ;\n",
"< \n",
"< // Bcon = \\mathcal{B}^\\mu of the paper:\n",
"< Bcon[0] = 0. ;\n",
"< for(i=1;i<NDIM;i++) Bcon[i] = -ncov[0] * prim[BCON1+i-1] ;\n",
"< \n",
"< u_dot_B = 0. ;\n",
"< for(i=0;i<NDIM;i++) u_dot_B += ucov[i]*Bcon[i] ;\n",
"< \n",
"< gamma = -ucon[0]*ncov[0] ;\n",
"< for(i=0;i<NDIM;i++) bcon[i] = (Bcon[i] + ucon[i]*u_dot_B)/gamma ;\n",
"< }\n",
"< \n",
"< \n",
"< /********************************************************************** \n",
"< gamma_calc_g(): \n",
"< \n",
"< -- using the primitive variables, contra-/co-variant 4-vel., \n",
"< and covariant normal vector, calculate the contravariant \n",
"< form of the magnetic 4-vector b^\\mu (the small \"b\" in HARM);\n",
"< \n",
"< -- assumes:\n",
"< \n",
"< / rho \\\n",
"< P = | uu |\n",
"< | \\tilde{u}^i |\n",
"< \\ B^i /\n",
"< ******************************************************************/\n",
"< int gamma_calc_g(CCTK_REAL *pr, CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL *gamma)\n",
"< {\n",
"< CCTK_REAL utsq ;\n",
"< \n",
"< utsq = gcov[1][1]*pr[UTCON1]*pr[UTCON1]\n",
"< + gcov[2][2]*pr[UTCON2]*pr[UTCON2]\n",
"< + gcov[3][3]*pr[UTCON3]*pr[UTCON3]\n",
"< + 2.*( gcov[1][2]*pr[UTCON1]*pr[UTCON2]\n",
"< + gcov[1][3]*pr[UTCON1]*pr[UTCON3]\n",
"< + gcov[2][3]*pr[UTCON2]*pr[UTCON3]) ;\n",
"< \n",
"< if(utsq<0.0){\n",
"< if(fabs(utsq)>1E-10){ // then assume not just machine precision\n",
"< return (1);\n",
"< }\n",
"< else utsq=1E-10; // set floor\n",
"< }\n",
"< \n",
"< *gamma = sqrt(1. + utsq) ;\n",
"< \n",
"< return(0) ;\n",
"< }\n",
"< \n",
"< \n",
"297,299c106,108\n",
"< /* \n",
"< pressure as a function of rho0 and u \n",
"< this is used by primtoU and Utoprim_?D\n",
"---\n",
"> /*\n",
"> pressure as a function of rho0 and u\n",
"> this is used by primtoU and Utoprim_?D\n",
"301c110\n",
"< static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u)\n",
"---\n",
"> static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u)\n",
"302a112,114\n",
"> \n",
"> // Set up Gamma_th:\n",
"> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
"304c116,129\n",
"< return((gamma_th /* <- Should be local polytropic Gamma factor */ - 1.)*u) ;\n",
"---\n",
"> #endif\n",
"> \n",
"> // Compute P_cold, eps_cold\n",
"> CCTK_REAL P_cold, eps_cold;\n",
"> compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);\n",
"> \n",
"> /* Compute the pressure as a function of rho_b (rho0) and\n",
"> * u = rho_b * eps, using our hybrid EOS:\n",
"> * .-------------------------------------------------------------.\n",
"> * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |\n",
"> * .-------------------------------------------------------------.\n",
"> */\n",
"> return( P_cold + (Gamma_th - 1.0)*(u - rho0*eps_cold) );\n",
"> \n",
"308,310c133,134\n",
"< \n",
"< /* \n",
"< pressure as a function of rho0 and w = rho0 + u + p \n",
"---\n",
"> /*\n",
"> pressure as a function of rho0 and w = rho0 + u + p\n",
"313c137\n",
"< static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w)\n",
"---\n",
"> static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w)\n",
"315,317d138\n",
"< DECLARE_CCTK_PARAMETERS;\n",
"< return((gamma_th /* <- Should be local polytropic Gamma factor */ -1.)*(w - rho0)/gamma_th /* <- Should be local polytropic Gamma factor */ ) ;\n",
"< }\n",
"318a140,143\n",
"> // Set up Gamma_th:\n",
"> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
"> DECLARE_CCTK_PARAMETERS;\n",
"> #endif\n",
"319a145,157\n",
"> // Compute P_cold, eps_cold\n",
"> CCTK_REAL P_cold, eps_cold;\n",
"> compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);\n",
"> \n",
"> /* Compute the pressure as a function of rho_b (rho0) and\n",
"> * w = u + rho_b + p, using our hybrid EOS:\n",
"> * ----------------------------------------------------------------------------\n",
"> * | p(rho_b,w) = ( P_cold + (Gamma_th-1)*( w - rho_b*(1+eps_cold) ) )/Gamma_th |\n",
"> * ----------------------------------------------------------------------------\n",
"> */\n",
"> return( (P_cold + (Gamma_th-1.0)*( w - rho0*(1.0+eps_cold) ) )/Gamma_th );\n",
"> }\n",
"> #endif\n"
]
}
],
"source": [
"# Verify if the code generated by this tutorial module\n",
"# matches the original IllinoisGRMHD source code\n",
"\n",
"# First download the original IllinoisGRMHD source code\n",
"import urllib\n",
"from os import path\n",
"\n",
"original_IGM_file_url = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/harm_u2p_util.c\"\n",
"original_IGM_file_name = \"harm_u2p_util-original.c\"\n",
"original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
"\n",
"# Then download the original IllinoisGRMHD source code\n",
"# We try it here in a couple of ways in an attempt to keep\n",
"# the code more portable\n",
"try:\n",
" original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
"except:\n",
" try:\n",
" original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
" except:\n",
" # If all else fails, hope wget does the job\n",
" !wget -O $original_IGM_file_path $original_IGM_file_url\n",
"\n",
"# Perform validation\n",
"Validation__harm_u2p_util__c = !diff $original_IGM_file_path $outfile_path__harm_u2p_util__c\n",
"\n",
"if Validation__harm_u2p_util__c == []:\n",
" # If the validation passes, we do not need to store the original IGM source code file\n",
" !rm $original_IGM_file_path\n",
" print(\"Validation test for harm_u2p_util.c: PASSED!\")\n",
"else:\n",
" # If the validation fails, we keep the original IGM source code file\n",
" print(\"Validation test for harm_u2p_util.c: FAILED!\")\n",
" # We also print out the difference between the code generated\n",
" # in this tutorial module and the original IGM source code\n",
" print(\"Diff:\")\n",
" for diff_line in Validation__harm_u2p_util__c:\n",
" print(diff_line)"
]
},
{
"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-IllinoisGRMHD__harm_u2p_util.pdf](Tutorial-IllinoisGRMHD__harm_u2p_util.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": 9,
"metadata": {},
"outputs": [],
"source": [
"latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n",
"#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__harm_u2p_util.ipynb\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex\n",
"!rm -f Tut*.out Tut*.aux Tut*.log"
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}