{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Tutorial-IllinoisGRMHD: harm_u2p_util.c\n", "\n", "## Authors: Leo Werneck & Zach Etienne\n", "\n", "**This module is currently under development**\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": [ "\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": [ "\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": [ "\n", "\n", "# Step 1: Introduction \\[Back to [top](#toc)\\]\n", "$$\\label{introduction}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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": [ "\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\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\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": [ "\n", "\n", "# Step 3: EOS dependent routines\n", "$$\\label{harm_utoprim_2d__c__eos_dep}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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": [ "\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": [ "\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 // 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 for(j=0;j \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 for(j=0;j /**********************************************************************\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;i1E-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": [ "\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 }