This module is currently under development
HARM
. This module will likely be absorbed by another one once we finish documenting the code.¶We will now use the cmdline_helper.py NRPy+ module to create the source directory within the IllinoisGRMHD
NRPy+ directory, if it does not exist yet.
# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
# Step 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)
# Step 0c: Create the output file path
outfile_path__harm_u2p_util__c = os.path.join(IGM_src_dir_path,"harm_u2p_util.c")
%%writefile $outfile_path__harm_u2p_util__c
#ifndef __HARM_U2P_UTIL__C__
#define __HARM_U2P_UTIL__C__
/*
-------------------------------------------------------------------------------
Copyright 2005 Scott C. Noble, Charles F. Gammie,
Jonathan C. McKinney, and Luca Del Zanna
This file is part of PVS-GRMHD.
PVS-GRMHD is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
PVS-GRMHD is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with PVS-GRMHD; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-------------------------------------------------------------------------------
*/
// Function prototypes for this file:
static void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM]);
static void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM]);
static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]);
static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u);
static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w);
// Inlined function used by this file
static inline void compute_P_cold__eps_cold(eos_struct eos,CCTK_REAL rho_in, CCTK_REAL &P_cold,CCTK_REAL &eps_cold);
Writing ../src/harm_u2p_util.c
raise_g()
function [Back to top]¶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
$$ \boxed{v^{\mu} = g^{\mu\nu}v_{\nu}}\ . $$%%writefile -a $outfile_path__harm_u2p_util__c
/**********************************************************************
raise_g():
-- calculates the contravariant form of a covariant tensor,
using the inverse of the metric;
***********************************************************************/
static void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM])
{
int i,j;
for(i=0;i<NDIM;i++) {
vcon[i] = 0. ;
for(j=0;j<NDIM;j++)
vcon[i] += gcon[i][j]*vcov[j] ;
}
return ;
}
Appending to ../src/harm_u2p_util.c
lower_g()
function [Back to top]¶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
$$ \boxed{v_{\mu} = g_{\mu\nu}v^{\nu}}\ . $$%%writefile -a $outfile_path__harm_u2p_util__c
/**********************************************************************
lower_g():
-- calculates the ocvariant form of a contravariant tensor
using the metric;
***********************************************************************/
static void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM])
{
int i,j;
for(i=0;i<NDIM;i++) {
vcov[i] = 0. ;
for(j=0;j<NDIM;j++)
vcov[i] += gcov[i][j]*vcon[j] ;
}
return ;
}
Appending to ../src/harm_u2p_util.c
%%writefile -a $outfile_path__harm_u2p_util__c
/**********************************************************************
ncov_calc():
-- calculates the covariant form of the normal vector to our
spacelike hypersurfaces ala the ADM formalism.
-- requires the inverse metric;
***********************************************************************/
static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM])
{
CCTK_REAL lapse ;
int i;
lapse = sqrt(-1./gcon[0][0]) ;
ncov[0] = -lapse ;
for( i = 1; i < NDIM; i++) {
ncov[i] = 0. ;
}
return ;
}
Appending to ../src/harm_u2p_util.c
pressure_rho0_u()
function [Back to top]¶The $\Gamma$-law EOS implemented in HARM
is
where
$$ u = \rho_{b}\epsilon\ . $$Thus, the pre-PPEOS Patch version of this function was
/**************************************************
The following functions assume a Gamma-law EOS:
***************************************************/
/*
pressure as a function of rho0 and u
this is used by primtoU and Utoprim_?D
*/
static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u)
{
return((GAMMA - 1.)*u) ;
}
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:
$$ 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] \implies \boxed{ 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]\ , } $$where
$$ \left\{ \begin{align} P_{\rm cold}\left(\rho_{b}\right) = K_{\rm poly}\rho_{b}^{\Gamma_{\rm poly}}\ ,\\ \epsilon_{\rm cold}\left(\rho_{b}\right) = C + \frac{P_{\rm cold}}{\rho_{b}\left(\Gamma_{\rm poly}-1\right)}\ , \end{align} \right. $$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.
%%writefile -a $outfile_path__harm_u2p_util__c
/**************************************************
The following functions assume a Gamma-law EOS:
***************************************************/
/*
pressure as a function of rho0 and u
this is used by primtoU and Utoprim_?D
*/
static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u)
{
// Set up Gamma_th:
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
// Compute P_cold, eps_cold
CCTK_REAL P_cold, eps_cold;
compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);
/* Compute the pressure as a function of rho_b (rho0) and
* u = rho_b * eps, using our hybrid EOS:
* .-------------------------------------------------------------.
* | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |
* .-------------------------------------------------------------.
*/
return( P_cold + (Gamma_th - 1.0)*(u - rho0*eps_cold) );
}
Appending to ../src/harm_u2p_util.c
pressure_rho0_w()
function [Back to top]¶The $\Gamma$-law EOS implemented in HARM
is
We want now to obtain $p_{\Gamma}\left(\rho_{b},w\right)$, where
$$ w = u + \rho_{b} + p\ . $$Then
$$ p_{\Gamma} = \left(\Gamma-1\right)\left(w - \rho_{b} - p_{\Gamma}\right) \implies \boxed{ p_{\Gamma}\left(\rho_{b},w\right) = \frac{\left(\Gamma-1\right)}{\Gamma}\left(w-\rho_{b}\right) }\ . $$Thus, the pre-PPEOS Patch version of this function was
/*
pressure as a function of rho0 and w = rho0 + u + p
this is used by primtoU and Utoprim_1D
*/
static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w)
{
return((GAMMA-1.)*(w - rho0)/GAMMA) ;
}
For our hybrid EOS, we have
$$ \begin{align} p_{\rm hybrid} &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\rho_{b}\left[\epsilon - \epsilon_{\rm cold}\right]\\ &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left[u - \rho_{b}\epsilon_{\rm cold}\right]\\ &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left[w-\rho_{b}-p_{\rm hybrid} - \rho_{b}\epsilon_{\rm cold}\right]\\ &= 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}\\ \implies &\boxed{ 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] } \end{align} $$%%writefile -a $outfile_path__harm_u2p_util__c
/*
pressure as a function of rho0 and w = rho0 + u + p
this is used by primtoU and Utoprim_1D
*/
static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w)
{
// Set up Gamma_th:
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
// Compute P_cold, eps_cold
CCTK_REAL P_cold, eps_cold;
compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);
/* Compute the pressure as a function of rho_b (rho0) and
* w = u + rho_b + p, using our hybrid EOS:
* ----------------------------------------------------------------------------
* | p(rho_b,w) = ( P_cold + (Gamma_th-1)*( w - rho_b*(1+eps_cold) ) )/Gamma_th |
* ----------------------------------------------------------------------------
*/
return( (P_cold + (Gamma_th-1.0)*( w - rho0*(1.0+eps_cold) ) )/Gamma_th );
}
#endif
Appending to ../src/harm_u2p_util.c
# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code
# First download the original IllinoisGRMHD source code
import urllib
from os import path
original_IGM_file_url = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/harm_u2p_util.c"
original_IGM_file_name = "harm_u2p_util-original.c"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)
# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
try:
original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
# If all else fails, hope wget does the job
!wget -O $original_IGM_file_path $original_IGM_file_url
# Perform validation
Validation__harm_u2p_util__c = !diff $original_IGM_file_path $outfile_path__harm_u2p_util__c
if Validation__harm_u2p_util__c == []:
# If the validation passes, we do not need to store the original IGM source code file
!rm $original_IGM_file_path
print("Validation test for harm_u2p_util.c: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for harm_u2p_util.c: FAILED!")
# We also print out the difference between the code generated
# in this tutorial module and the original IGM source code
print("Diff:")
for diff_line in Validation__harm_u2p_util__c:
print(diff_line)
Validation test for harm_u2p_util.c: FAILED! Diff: 1c1,2 < --- > #ifndef __HARM_U2P_UTIL__C__ > #define __HARM_U2P_UTIL__C__ 4c5 < Copyright 2005 Scott C. Noble, Charles F. Gammie, --- > Copyright 2005 Scott C. Noble, Charles F. Gammie, 28,144c29,33 < void primtoU_g( CCTK_REAL prim[], CCTK_REAL gcov[][NDIM], CCTK_REAL gcon[][NDIM], CCTK_REAL gdet, CCTK_REAL U[] ); < static void ucon_calc_g(CCTK_REAL prim[],CCTK_REAL gcov[][NDIM],CCTK_REAL gcon[][NDIM],CCTK_REAL ucon[]); < static void raise_g(CCTK_REAL vcov[], CCTK_REAL gcon[][NDIM], CCTK_REAL vcon[]); < static void lower_g(CCTK_REAL vcon[], CCTK_REAL gcov[][NDIM], CCTK_REAL vcov[]); < static void ncov_calc(CCTK_REAL gcon[][NDIM],CCTK_REAL ncov[]) ; < static void bcon_calc_g(CCTK_REAL prim[],CCTK_REAL ucon[],CCTK_REAL ucov[],CCTK_REAL ncov[],CCTK_REAL bcon[]); < static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u); < static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w); < int gamma_calc_g(CCTK_REAL *pr, CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL *gamma); < < /********************************************************************** < primtoU_g(): < < -- calculates the conserved variables from the primitive variables < and the metric; < -- assumes that the conserved and primitive variables are defined ala HARM: < < / rho u^t \ < U = | T^t_\mu + rho u^t | sqrt(-det(g_{\mu\nu})) < \ B^i / < < / rho \ < P = | uu | < | \tilde{u}^i | < \ B^i / < < **************************************************************************/ < < < void primtoU_g( < CCTK_REAL prim[NPR], /* primitive variables */ < CCTK_REAL gcov[NDIM][NDIM], /* covariant (index dn) form of metric */ < CCTK_REAL gcon[NDIM][NDIM], /* contravariant (index up) form of metric */ < CCTK_REAL gdet, /* sqrt of -1 times det(g_{\mu \nu}) */ < CCTK_REAL U[NPR] /* matrix of derivatives */ < ) { < int i ; < CCTK_REAL rho0 ; < static CCTK_REAL ucon[NDIM],ucov[NDIM],bcon[NDIM],bcov[NDIM],ncov[NDIM] ; < CCTK_REAL gamma,n_dot_b,bsq,u,p,w, alpha ; < < < /* Calculate auxiliary quantities: */ < alpha = 1.0/sqrt(-gcon[0][0]); < < ucon_calc_g(prim,gcov,gcon,ucon) ; < lower_g(ucon,gcov,ucov) ; < ncov_calc(gcon,ncov) ; < < gamma = -ncov[0]*ucon[0] ; < < bcon_calc_g(prim,ucon,ucov,ncov,bcon) ; < lower_g(bcon,gcov,bcov) ; < < n_dot_b = 0. ; < for(i=0;i<4;i++) n_dot_b += ncov[i]*bcon[i] ; < bsq = 0. ; < for(i=0;i<4;i++) bsq += bcov[i]*bcon[i] ; < < rho0 = prim[RHO] ; < u = prim[UU] ; < p = pressure_rho0_u(rho0,u) ; < w = rho0 + u + p ; < < // Now set the conserved variables themselves, using HARM's definition: < U[RHO] = ucon[0]*rho0 ; < < for( i = 0; i < 4; i++) { < U[QCOV0+i] = gamma*(w + bsq)*ucov[i] < - (p + bsq/2.)*ncov[i] < + n_dot_b*bcov[i] ; < < U[QCOV0+i] /= alpha; < } < < U[QCOV0] = U[QCOV0] + U[RHO]; < U[BCON1] = prim[BCON1] ; < U[BCON2] = prim[BCON2] ; < U[BCON3] = prim[BCON3] ; < < for(i = 0; i < NPR; i++ ) { < U[i] *= gdet; < } < < return ; < } < < /********************************************************************** < ucon_calc_g(): < < -- calculates the contravariant (up) components of the four-velocity < given the primitive variables, of which the velocity is < \tilde{u}^i = \gamma v^j where v^j is the velocity of the < flow w.r.t a normal observer to the coordinates; < < -- also requires the metric and inverse metric; < < -- assumes: < < / rho \ < P = | uu | < | \tilde{u}^i | < \ B^i / < < ******************************************************************/ < static void ucon_calc_g(CCTK_REAL prim[NPR],CCTK_REAL gcov[NDIM][NDIM],CCTK_REAL gcon[NDIM][NDIM], < CCTK_REAL ucon[NDIM]) < { < CCTK_REAL u_tilde_con[4] ; < CCTK_REAL u_tilde_sq ; < CCTK_REAL gamma,lapse ; < int i,j; < < u_tilde_con[0] = 0. ; < u_tilde_con[1] = prim[UTCON1] ; < u_tilde_con[2] = prim[UTCON2] ; < u_tilde_con[3] = prim[UTCON3] ; --- > static void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM]); > static void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM]); > static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]); > static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u); > static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w); 146,156c35,36 < u_tilde_sq = 0. ; < for(i=0;i<NDIM;i++) < for(j=0;j<NDIM;j++) < u_tilde_sq += gcov[i][j]*u_tilde_con[i]*u_tilde_con[j] ; < u_tilde_sq = fabs(u_tilde_sq) ; < < gamma = sqrt(1. + u_tilde_sq) ; < < lapse = sqrt(-1./gcon[0][0]) ; < < for(i=0;i<NDIM;i++) ucon[i] = u_tilde_con[i] - lapse*gamma*gcon[0][i] ; --- > // Inlined function used by this file > static inline void compute_P_cold__eps_cold(eos_struct eos,CCTK_REAL rho_in, CCTK_REAL &P_cold,CCTK_REAL &eps_cold); 158,159d37 < return ; < } 161c39 < /********************************************************************** --- > /********************************************************************** 163,164c41,42 < < -- calculates the contravariant form of a covariant tensor, --- > > -- calculates the contravariant form of a covariant tensor, 166c44 < ******************************************************************/ --- > ***********************************************************************/ 173c51 < for(j=0;j<NDIM;j++) --- > for(j=0;j<NDIM;j++) 180c58,59 < /********************************************************************** --- > > /********************************************************************** 182,183c61,62 < < -- calculates the ocvariant form of a contravariant tensor --- > > -- calculates the ocvariant form of a contravariant tensor 185c64 < ******************************************************************/ --- > ***********************************************************************/ 192c71 < for(j=0;j<NDIM;j++) --- > for(j=0;j<NDIM;j++) 199,200d77 < /********************************************************************** < ncov_calc(): 202c79,82 < -- calculates the covariant form of the normal vector to our --- > /********************************************************************** > ncov_calc(): > > -- calculates the covariant form of the normal vector to our 206,207c86,87 < ******************************************************************/ < static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]) --- > ***********************************************************************/ > static void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]) 215c95 < for( i = 1; i < NDIM; i++) { --- > for( i = 1; i < NDIM; i++) { 222,292d101 < /********************************************************************** < bcon_calc_g(): < < -- using the primitive variables, contra-/co-variant 4-vel., < and covariant normal vector, calculate the contravariant < form of the magnetic 4-vector b^\mu (the small "b" in HARM); < -- assumes: < < / rho \ < P = | uu | < | \tilde{u}^i | < \ B^i / < ******************************************************************/ < static void bcon_calc_g(CCTK_REAL prim[NPR],CCTK_REAL ucon[NDIM],CCTK_REAL ucov[NDIM], < CCTK_REAL ncov[NDIM],CCTK_REAL bcon[NDIM]) < { < static CCTK_REAL Bcon[NDIM] ; < CCTK_REAL u_dot_B ; < CCTK_REAL gamma ; < int i ; < < // Bcon = \mathcal{B}^\mu of the paper: < Bcon[0] = 0. ; < for(i=1;i<NDIM;i++) Bcon[i] = -ncov[0] * prim[BCON1+i-1] ; < < u_dot_B = 0. ; < for(i=0;i<NDIM;i++) u_dot_B += ucov[i]*Bcon[i] ; < < gamma = -ucon[0]*ncov[0] ; < for(i=0;i<NDIM;i++) bcon[i] = (Bcon[i] + ucon[i]*u_dot_B)/gamma ; < } < < < /********************************************************************** < gamma_calc_g(): < < -- using the primitive variables, contra-/co-variant 4-vel., < and covariant normal vector, calculate the contravariant < form of the magnetic 4-vector b^\mu (the small "b" in HARM); < < -- assumes: < < / rho \ < P = | uu | < | \tilde{u}^i | < \ B^i / < ******************************************************************/ < int gamma_calc_g(CCTK_REAL *pr, CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL *gamma) < { < CCTK_REAL utsq ; < < utsq = gcov[1][1]*pr[UTCON1]*pr[UTCON1] < + gcov[2][2]*pr[UTCON2]*pr[UTCON2] < + gcov[3][3]*pr[UTCON3]*pr[UTCON3] < + 2.*( gcov[1][2]*pr[UTCON1]*pr[UTCON2] < + gcov[1][3]*pr[UTCON1]*pr[UTCON3] < + gcov[2][3]*pr[UTCON2]*pr[UTCON3]) ; < < if(utsq<0.0){ < if(fabs(utsq)>1E-10){ // then assume not just machine precision < return (1); < } < else utsq=1E-10; // set floor < } < < *gamma = sqrt(1. + utsq) ; < < return(0) ; < } < < 297,299c106,108 < /* < pressure as a function of rho0 and u < this is used by primtoU and Utoprim_?D --- > /* > pressure as a function of rho0 and u > this is used by primtoU and Utoprim_?D 301c110 < static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u) --- > static CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u) 302a112,114 > > // Set up Gamma_th: > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 304c116,129 < return((gamma_th /* <- Should be local polytropic Gamma factor */ - 1.)*u) ; --- > #endif > > // Compute P_cold, eps_cold > CCTK_REAL P_cold, eps_cold; > compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold); > > /* Compute the pressure as a function of rho_b (rho0) and > * u = rho_b * eps, using our hybrid EOS: > * .-------------------------------------------------------------. > * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) | > * .-------------------------------------------------------------. > */ > return( P_cold + (Gamma_th - 1.0)*(u - rho0*eps_cold) ); > 308,310c133,134 < < /* < pressure as a function of rho0 and w = rho0 + u + p --- > /* > pressure as a function of rho0 and w = rho0 + u + p 313c137 < static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w) --- > static CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w) 315,317d138 < DECLARE_CCTK_PARAMETERS; < return((gamma_th /* <- Should be local polytropic Gamma factor */ -1.)*(w - rho0)/gamma_th /* <- Should be local polytropic Gamma factor */ ) ; < } 318a140,143 > // Set up Gamma_th: > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER > DECLARE_CCTK_PARAMETERS; > #endif 319a145,157 > // Compute P_cold, eps_cold > CCTK_REAL P_cold, eps_cold; > compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold); > > /* Compute the pressure as a function of rho_b (rho0) and > * w = u + rho_b + p, using our hybrid EOS: > * ---------------------------------------------------------------------------- > * | p(rho_b,w) = ( P_cold + (Gamma_th-1)*( w - rho_b*(1+eps_cold) ) )/Gamma_th | > * ---------------------------------------------------------------------------- > */ > return( (P_cold + (Gamma_th-1.0)*( w - rho0*(1.0+eps_cold) ) )/Gamma_th ); > } > #endif
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 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).
latex_nrpy_style_path = os.path.join(nrpy_dir_path,"latex_nrpy_style.tplx")
#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__harm_u2p_util.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex
!rm -f Tut*.out Tut*.aux Tut*.log