This module is currently under development
HARM
.¶This module is organized as follows
Utoprim_2d()
functionUtoprim_new_body()
functionUtoprim_new_body()
functionvsq_calc()
functionx1_of_x0()
functionvalidate_x()
functiongeneral_newton_raphson()
functionfunc_vsq()
functionpressure_W_vsq()
functiondpdW_calc_vsq()
functiondpdvsq_calc()
functionWe 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_utoprim_2d__c = os.path.join(IGM_src_dir_path,"harm_utoprim_2d.c")
Let us now start documenting the harm_utoprim_2d.c
, which is a part of the Harm
code. Our main reference throughout this discussion will be the required citation Noble et al. (2006).
We will start with the code's required preamble.
%%writefile $outfile_path__harm_utoprim_2d__c
#ifndef __HARM_UTOPRIM_2D__C__
#define __HARM_UTOPRIM_2D__C__
/***********************************************************************************
Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble,
Gabor Toth, and Luca Del Zanna
HARM version 1.0 (released May 1, 2006)
This file is part of HARM. HARM is a program that solves hyperbolic
partial differential equations in conservative form using high-resolution
shock-capturing techniques. This version of HARM has been configured to
solve the relativistic magnetohydrodynamic equations of motion on a
stationary black hole spacetime in Kerr-Schild coordinates to evolve
an accretion disk model.
You are morally obligated to cite the following two papers in his/her
scientific literature that results from use of any part of HARM:
[1] Gammie, C. F., McKinney, J. C., \& Toth, G.\ 2003,
Astrophysical Journal, 589, 444.
[2] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006,
Astrophysical Journal, 641, 626.
Further, we strongly encourage you to obtain the latest version of
HARM directly from our distribution website:
http://rainman.astro.uiuc.edu/codelib/
HARM 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.
HARM 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 HARM; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
***********************************************************************************/
/*************************************************************************************/
/*************************************************************************************/
/*************************************************************************************
utoprim_2d.c:
---------------
Uses the 2D method:
-- solves for two independent variables (W,v^2) via a 2D
Newton-Raphson method
-- can be used (in principle) with a general equation of state.
-- Currently returns with an error state (>0) if a negative rest-mass
density or internal energy density is calculated. You may want
to change this aspect of the code so that it still calculates the
velocity and so that you can floor the densities. If you want to
change this aspect of the code please comment out the "return(retval)"
statement after "retval = 5;" statement in Utoprim_new_body();
******************************************************************************/
static const int NEWT_DIM=2;
// Declarations:
static CCTK_REAL vsq_calc(CCTK_REAL W,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);
static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[],long &n_iter);
static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter, void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *, CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);
static void func_vsq( eos_struct eos, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);
static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D ) ;
static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) ;
static CCTK_REAL dpdW_calc_vsq(CCTK_REAL W, CCTK_REAL vsq);
static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D);
/**********************************************************************/
/******************************************************************
Utoprim_2d():
-- Driver for new prim. var. solver. The driver just translates
between the two sets of definitions for U and P. The user may
wish to alter the translation as they see fit. Note that Greek
indices run 0,1,2,3 and Latin indices run 1,2,3 (spatial only).
/ rho u^t \
U = | T^t_t + rho u^t | sqrt(-det(g_{\mu\nu}))
| T^t_i |
\ B^i /
/ rho \
P = | uu |
| \tilde{u}^i |
\ B^i /
Arguments:
U[NPR] = conserved variables (current values on input/output);
gcov[NDIM][NDIM] = covariant form of the metric ;
gcon[NDIM][NDIM] = contravariant form of the metric ;
gdet = sqrt( - determinant of the metric) ;
prim[NPR] = primitive variables (guess on input, calculated values on
output if there are no problems);
-- NOTE: for those using this routine for special relativistic MHD and are
unfamiliar with metrics, merely set
gcov = gcon = diag(-1,1,1,1) and gdet = 1. ;
******************************************************************/
Writing ../src/harm_utoprim_2d.c
Utoprim_2d()
function [Back to top]¶The Utoprim_2d()
function is the driver function of the HARM
conservative-to-primitive algorithm. We remind you from the definitions of primitive and conservative variables used in the code:
Let
$$ \tilde{B}^{i}_{\rm HARM} \equiv \sqrt{-g}B^{i}_{\rm HARM}\ . $$The code starts by relating
$$ \boxed{B^{i}_{\rm HARM} = \frac{\tilde{B}^{i}_{\rm HARM}}{\sqrt{-g}}}\ , $$and setting
$$ \boxed{\alpha = \frac{1}{\sqrt{-g^{00}}}} \ . $$%%writefile -a $outfile_path__harm_utoprim_2d__c
int Utoprim_2d(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM],
CCTK_REAL gdet, CCTK_REAL prim[NPR], long &n_iter)
{
CCTK_REAL U_tmp[NPR], prim_tmp[NPR];
int i, ret;
CCTK_REAL alpha;
if( U[0] <= 0. ) {
return(-100);
}
/* First update the primitive B-fields */
for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] / gdet ;
/* Set the geometry variables: */
alpha = 1.0/sqrt(-gcon[0][0]);
Appending to ../src/harm_utoprim_2d.c
Utoprim_new_body()
function [Back to top]¶The conservative-to-primitive algorithm uses the Utoprim_new_body()
function. However, this function assumes a different set of primitive/conservative variables. Thus, we must perform the proper conversion. First, let us ease on the notation a bit by defining:
Below we list the main differences in the conservative variables:
Utoprim_2d() |
Utoprim_new_body() |
---|---|
$\color{blue}{\textbf{Conservatives}}$ | $\color{red}{\textbf{Conservatives}}$ |
$\color{blue}{\rho_{\star}}$ | $\color{red}{\frac{\alpha}{\sqrt{-g}}\rho_{\star}}$ |
$\color{blue}{u_{\star}}$ | $\color{red}{\frac{\alpha}{\sqrt{-g}}\left(u_{\star}-\rho_{\star}\right)}$ |
$\color{blue}{\tilde{S}_{i}}$ | $\color{red}{\frac{\alpha}{\sqrt{-g}}\tilde{S}_{i}}$ |
$\color{blue}{\tilde{B}^{i}_{\rm HARM}}$ | $\color{red}{\frac{\alpha}{\sqrt{-g}}\tilde{B}^{i}_{\rm HARM}}$ |
These are necessary conversions because while Utoprim_2d()
assumes the set of conservatives above, Utoprim_new_body()
assumes
Let us first pause to understand the table above. From definition (15) in Noble et al. (2006) and the discussion just below it, we know that $\gamma = \alpha u^{0}$. Thus
$$ \rho_{\star} = \sqrt{-g}\rho_{b}u^{0} = \sqrt{-g}\left(\frac{\gamma}{\alpha}\rho_{b}\right)\implies\boxed{\gamma \rho_{b} = \frac{\alpha}{\sqrt{-g}}\rho_{\star}}\ . $$Then we have
$$ u_{\star} = \sqrt{-g}\left(T^{0}_{\ \ 0} + \rho_{b}u^{0}\right)= \sqrt{-g}\left(T^{0}_{\ \ 0} + \frac{\rho_{\star}}{\sqrt{-g}}\right) = \sqrt{-g}T^{0}_{\ \ 0} + \rho_{\star} \implies \boxed{\alpha T^{0}_{\ \ 0} = \frac{\alpha}{\sqrt{-g}}\left(u_{\star}-\rho_{\star}\right)}\ . $$The other two relations are more straightforward. We have
$$ \tilde{S}_{i} = \sqrt{-g}T^{0}_{\ \ i} \implies \boxed{\alpha T^{0}_{\ \ i} = \frac{\alpha}{\sqrt{-g}}\tilde{S}_{i}}\ , $$and
$$ \tilde{B}^{i}_{\rm HARM} = \sqrt{-g}B^{i}_{\rm HARM}\implies \boxed{\alpha B^{i}_{\rm HARM} = \frac{\alpha}{\sqrt{-g}}\tilde{B}^{i}_{\rm HARM}}\ . $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/* Transform the CONSERVED variables into the new system */
U_tmp[RHO] = alpha * U[RHO] / gdet;
U_tmp[UU] = alpha * (U[UU] - U[RHO]) / gdet ;
for( i = UTCON1; i <= UTCON3; i++ ) {
U_tmp[i] = alpha * U[i] / gdet ;
}
for( i = BCON1; i <= BCON3; i++ ) {
U_tmp[i] = alpha * U[i] / gdet ;
}
Appending to ../src/harm_utoprim_2d.c
Below we list the necessary transformations on the primitive variables:
Utoprim_2d() |
Utoprim_new_body() |
---|---|
$\color{blue}{\textbf{Primitives}}$ | $\color{red}{\textbf{Primitives}}$ |
$\color{blue}{\rho_{b}}$ | $\color{red}{\rho_{b}}$ |
$\color{blue}{u}$ | $\color{red}{u}$ |
$\color{blue}{\tilde{u}^{i}}$ | $\color{red}{\tilde{u}^{i}}$ |
$\color{blue}{B^{i}_{\rm HARM}}$ | $\color{red}{\alpha B^{i}_{\rm HARM}}$ |
After this slight modification, we call the Utoprim_new_body()
function. If it returns without errors, then the variables ${\rm prim\_tmp}$ will now contain the values of the primitives. We then update the ${\rm prim}$ variables with these newly computed values.
%%writefile -a $outfile_path__harm_utoprim_2d__c
/* Transform the PRIMITIVE variables into the new system */
for( i = 0; i < BCON1; i++ ) {
prim_tmp[i] = prim[i];
}
for( i = BCON1; i <= BCON3; i++ ) {
prim_tmp[i] = alpha*prim[i];
}
ret = Utoprim_new_body(eos, U_tmp, gcov, gcon, gdet, prim_tmp,n_iter);
/* Transform new primitive variables back if there was no problem : */
if( ret == 0 || ret == 5 || ret==101 ) {
for( i = 0; i < BCON1; i++ ) {
prim[i] = prim_tmp[i];
}
}
return( ret ) ;
}
Appending to ../src/harm_utoprim_2d.c
%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************/
/**********************************************************************************
Utoprim_new_body():
-- Attempt an inversion from U to prim using the initial guess prim.
-- This is the main routine that calculates auxiliary quantities for the
Newton-Raphson routine.
-- assumes that
/ rho gamma \
U = | alpha T^t_\mu |
\ alpha B^i /
/ rho \
prim = | uu |
| \tilde{u}^i |
\ alpha B^i /
return: (i*100 + j) where
i = 0 -> Newton-Raphson solver either was not called (yet or not used)
or returned successfully;
1 -> Newton-Raphson solver did not converge to a solution with the
given tolerances;
2 -> Newton-Raphson procedure encountered a numerical divergence
(occurrence of "nan" or "+/-inf" ;
j = 0 -> success
1 -> failure: some sort of failure in Newton-Raphson;
2 -> failure: utsq<0 w/ initial p[] guess;
3 -> failure: W<0 or W>W_TOO_BIG
4 -> failure: v^2 > 1
5 -> failure: rho,uu <= 0 ;
**********************************************************************************/
static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM],
CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[NPR], long &n_iter)
{
CCTK_REAL x_2d[NEWT_DIM];
CCTK_REAL QdotB,Bcon[NDIM],Bcov[NDIM],Qcov[NDIM],Qcon[NDIM],ncov[NDIM],ncon[NDIM],Qsq,Qtcon[NDIM];
CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,utsq,vsq;
int i,j, n, retval, i_increase;
n = NEWT_DIM ;
// Assume ok initially:
retval = 0;
Appending to ../src/harm_utoprim_2d.c
We start by computing basic quantities from the input variables. Notice that this conservative-to-primitive algorithm does not need to update the magnetic field, thus
$$ \boxed{B_{\rm prim}^{i} = B_{\rm conserv}^{i}}\ . $$Since they are both equal, we will not distinguish between prim and conserve in what follows. We also set $B^{0} = 0$. Then we define
$$ \boxed{Q_{\mu} \equiv \alpha T^{0}_{\ \ \mu}}\ . $$From these, the following quantities are then computed:
$$ \boxed{ \begin{align} B_{i} &= g_{i\mu}B^{\mu}\\ Q^{\mu} &= g^{\mu\nu}Q_{\nu}\\ B^{2} &= B_{i}B^{i}\\ Q\cdot B &= Q_{\mu}B^{\mu}\\ \left(Q\cdot B\right)^{2} &= \left(Q\cdot B\right)\left(Q\cdot B\right)\\ n_{\mu} &= \left(-\alpha,0,0,0\right)\\ n^{\mu} &= g^{\mu\nu}n_{\nu}\\ \left(Q\cdot n\right) &= Q^{\mu}n_{\mu}\\ Q^{2} &= Q_{\mu}Q^{\mu}\\ \tilde{Q}^{2} &= Q^{2} + \left(Q\cdot n\right)\left(Q\cdot n\right)\\ D &\equiv \gamma \rho_{b} \end{align} }\ . $$%%writefile -a $outfile_path__harm_utoprim_2d__c
for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;
// Calculate various scalars (Q.B, Q^2, etc) from the conserved variables:
Bcon[0] = 0. ;
for(i=1;i<4;i++) Bcon[i] = U[BCON1+i-1] ;
lower_g(Bcon,gcov,Bcov) ;
for(i=0;i<4;i++) Qcov[i] = U[QCOV0+i] ;
raise_g(Qcov,gcon,Qcon) ;
CCTK_REAL Bsq = 0. ;
for(i=1;i<4;i++) Bsq += Bcon[i]*Bcov[i] ;
QdotB = 0. ;
for(i=0;i<4;i++) QdotB += Qcov[i]*Bcon[i] ;
CCTK_REAL QdotBsq = QdotB*QdotB ;
ncov_calc(gcon,ncov) ;
// FIXME: The exact form of n^{\mu} can be found
// in eq. (2.116) and implementing it
// directly is a lot more efficient than
// performing n^{\mu} = g^{\mu\nu}n_{nu}
raise_g(ncov,gcon,ncon);
CCTK_REAL Qdotn = Qcon[0]*ncov[0] ;
Qsq = 0. ;
for(i=0;i<4;i++) Qsq += Qcov[i]*Qcon[i] ;
CCTK_REAL Qtsq = Qsq + Qdotn*Qdotn ;
CCTK_REAL D = U[RHO] ;
Appending to ../src/harm_utoprim_2d.c
The quantity $W$ is defined as
$$ W \equiv w\gamma^{2}\ , $$where
$$ \begin{align} w &= \rho_{b} + u + p\ ,\\ \gamma^{2} &= 1 + g_{ij}\tilde{u}^{i}\tilde{u}^{j}\ . \end{align} $$Thus the quantities $g_{ij}\tilde{u}^{i}\tilde{u}^{j}$ and then $\gamma^{2}$ and $\gamma$. Thus, by computing $\rho_{b}$ and $p$ from the input variables, i.e. $D$, one can determine $w$ and then compute the value of $W$ from the input values (previous iteration), which we denote by $W_{\rm last}$.
Dependecy note: Note that this function depends on the pressure_rho0_u()
function, which is not EOS independent.
%%writefile -a $outfile_path__harm_utoprim_2d__c
/* calculate W from last timestep and use for guess */
utsq = 0. ;
for(i=1;i<4;i++)
for(j=1;j<4;j++) utsq += gcov[i][j]*prim[UTCON1+i-1]*prim[UTCON1+j-1] ;
if( (utsq < 0.) && (fabs(utsq) < 1.0e-13) ) {
utsq = fabs(utsq);
}
if(utsq < 0. || utsq > UTSQ_TOO_BIG) {
retval = 2;
return(retval) ;
}
gammasq = 1. + utsq ;
gamma = sqrt(gammasq);
// Always calculate rho from D and gamma so that using D in EOS remains consistent
// i.e. you don't get positive values for dP/d(vsq) .
rho0 = D / gamma ;
u = prim[UU] ;
p = pressure_rho0_u(eos, rho0,u) ;
w = rho0 + u + p ;
W_last = w*gammasq ;
// Make sure that W is large enough so that v^2 < 1 :
i_increase = 0;
while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq )
- QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq))
&& (i_increase < 10) ) {
W_last *= 10.;
i_increase++;
}
Appending to ../src/harm_utoprim_2d.c
Then we use equation (28) in Noble et al. (2006) to determine $v^{2}$:
$$ \boxed{v^{2} = \frac{\tilde{Q}^{2}W^{2} + \left(Q\cdot B\right)^{2}\left(B^{2}+2W\right)}{\left(B^{2}+W\right)^{2}W^{2}}}\ . $$This is done by calling the x1_of_x0()
function, where $x_{0} = W$ and $x_{1} = v^{2}$, which itself calls the vsq_calc()
function which implements the boxed equation above.
After we have $\left\{W_{\rm last},v^{2}_{\rm last}\right\}$ we use them as the initial guess for the general_newton_raphson()
, which returns the updated values $\left\{W,v^{2}\right\}$.
All functions mentioned above are documented in this tutorial notebook, so look at the Table of Contents for more information.
%%writefile -a $outfile_path__harm_utoprim_2d__c
// Calculate W and vsq:
x_2d[0] = fabs( W_last );
x_2d[1] = x1_of_x0( W_last , Bsq,QdotBsq,Qtsq,Qdotn,D) ;
retval = general_newton_raphson( eos, x_2d, n, n_iter, func_vsq, Bsq,QdotBsq,Qtsq,Qdotn,D) ;
W = x_2d[0];
vsq = x_2d[1];
/* Problem with solver, so return denoting error before doing anything further */
if( (retval != 0) || (W == FAIL_VAL) ) {
retval = retval*100+1;
return(retval);
}
else{
if(W <= 0. || W > W_TOO_BIG) {
retval = 3;
return(retval) ;
}
}
// Calculate v^2:
if( vsq >= 1. ) {
vsq = 1.-2.e-16;
//retval = 4;
//return(retval) ;
}
Appending to ../src/harm_utoprim_2d.c
Now that we have $\left\{W,v^{2}\right\}$, we recompute the primitive variables. We start with
$$ \left\{ \begin{align} \tilde{g} &\equiv \sqrt{1-v^{2}}\\ \gamma &= \frac{1}{\tilde{g}} \end{align} \right. \implies \boxed{\rho_{b} = D\tilde{g}}\ . $$Then, we determine the pressure $p$ using the pressure_rho0_w()
function and
Dependecy note: Note that this function depends on the pressure_rho0_w()
function, which is not EOS independent.
Finally, we can obtain $\tilde{u}^{i}$ using eq. 31 in Noble et al. (2006)
$$ \boxed{ \tilde{u}^{i} = \frac{\gamma}{\left(W+B^{2}\right)}\left[\tilde{Q}^{i} + \frac{\left(Q\cdot B\right)}{W}B^{i}\right] }\ , $$where
$$ \tilde{Q}^{i} = Q^{i} + \left(Q\cdot n\right)n^{i}\ . $$%%writefile -a $outfile_path__harm_utoprim_2d__c
// Recover the primitive variables from the scalars and conserved variables:
gtmp = sqrt(1. - vsq);
gamma = 1./gtmp ;
rho0 = D * gtmp;
w = W * (1. - vsq) ;
p = pressure_rho0_w(eos, rho0,w) ;
u = w - (rho0 + p) ; // u = rho0 eps, w = rho0 h
if( (rho0 <= 0.) || (u <= 0.) ) {
// User may want to handle this case differently, e.g. do NOT return upon
// a negative rho/u, calculate v^i so that rho/u can be floored by other routine:
retval = 5;
//return(retval) ;
}
/*
if(retval==5 && fabs(u)<1e-16) {
u = fabs(u);
CCTK_VInfo(CCTK_THORNSTRING,"%e\t%e\t%e",1.0-w/(rho0 + p),rho0,p);
retval=0;
}
*/
prim[RHO] = rho0 ;
prim[UU] = u ;
for(i=1;i<4;i++) Qtcon[i] = Qcon[i] + ncon[i] * Qdotn;
for(i=1;i<4;i++) prim[UTCON1+i-1] = gamma/(W+Bsq) * ( Qtcon[i] + QdotB*Bcon[i]/W ) ;
/* set field components */
for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;
/* done! */
return(retval) ;
}
Appending to ../src/harm_utoprim_2d.c
vsq_calc()
function [Back to top]¶This function implements eq. (28) in Noble et al. (2006) to determine $v^{2}$:
$$ \boxed{v^{2} = \frac{\tilde{Q}^{2}W^{2} + \left(Q\cdot B\right)^{2}\left(B^{2}+2W\right)}{\left(B^{2}+W\right)^{2}W^{2}}}\ . $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************/
/****************************************************************************
vsq_calc():
-- evaluate v^2 (spatial, normalized velocity) from
W = \gamma^2 w
****************************************************************************/
static CCTK_REAL vsq_calc(CCTK_REAL W,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D)
{
CCTK_REAL Wsq,Xsq;
Wsq = W*W ;
Xsq = (Bsq + W) * (Bsq + W);
return( ( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq) );
}
Appending to ../src/harm_utoprim_2d.c
%%writefile -a $outfile_path__harm_utoprim_2d__c
/********************************************************************
x1_of_x0():
-- calculates v^2 from W with some physical bounds checking;
-- asumes x0 is already physical
-- makes v^2 physical if not;
*********************************************************************/
static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D )
{
CCTK_REAL vsq;
CCTK_REAL dv = 1.e-15;
vsq = fabs(vsq_calc(x0,Bsq,QdotBsq,Qtsq,Qdotn,D)) ; // guaranteed to be positive
return( ( vsq > 1. ) ? (1.0 - dv) : vsq );
}
Appending to ../src/harm_utoprim_2d.c
%%writefile -a $outfile_path__harm_utoprim_2d__c
/********************************************************************
validate_x():
-- makes sure that x[0,1] have physical values, based upon
their definitions:
*********************************************************************/
static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] )
{
CCTK_REAL dv = 1.e-15;
/* Always take the absolute value of x[0] and check to see if it's too big: */
x[0] = fabs(x[0]);
x[0] = (x[0] > W_TOO_BIG) ? x0[0] : x[0];
x[1] = (x[1] < 0.) ? 0. : x[1]; /* if it's too small */
x[1] = (x[1] > 1.) ? (1. - dv) : x[1]; /* if it's too big */
return;
}
Appending to ../src/harm_utoprim_2d.c
general_newton_raphson()
function [Back to top]¶This function implements a multidimensional Newton-Raphson method. We will not make the effort of explaining the algorithm exhaustively since it is pretty standard, so we will settle for a summary of the method.
Given a system of $N$ non-linear of equations and $N$ variables, $\left\{\vec{F}\!\left(\vec{x}\right),\vec{x}\right\}$, the Newton-Raphson method attempts to determine the root vector, $\vec{x}_{\star}$, iteratively through
$$ \begin{align} \vec{x}_{n+1} = \vec{x}_{n} - J^{-1}_{F}\!\left(\vec{x}_{n}\right)\vec{F}\!\left(\vec{x}\right)\ , \end{align} $$where $J^{-1}_{F}$ is the Jacobian matrix
$$ \left(J_{F}\right)^{i}_{\ \ j} = \frac{\partial F^{i}}{\partial x^{j}}\ . $$The index $n$ above is an iteration index and $\vec{x}_{n+1}$ represents an improved approximation to $\vec{x}_{\star}$ when compared to $\vec{x}_{n}$.
%%writefile -a $outfile_path__harm_utoprim_2d__c
/************************************************************
general_newton_raphson():
-- performs Newton-Rapshon method on an arbitrary system.
-- inspired in part by Num. Rec.'s routine newt();
*****************************************************************/
static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter,
void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [],
CCTK_REAL [][NEWT_DIM], CCTK_REAL *,
CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D)
{
CCTK_REAL f, df, dx[NEWT_DIM], x_old[NEWT_DIM];
CCTK_REAL resid[NEWT_DIM], jac[NEWT_DIM][NEWT_DIM];
CCTK_REAL errx, x_orig[NEWT_DIM];
int id, i_extra, doing_extra;
int keep_iterating;
// Initialize various parameters and variables:
errx = 1. ;
df = f = 1.;
i_extra = doing_extra = 0;
for( id = 0; id < n ; id++) x_old[id] = x_orig[id] = x[id] ;
n_iter = 0;
/* Start the Newton-Raphson iterations : */
keep_iterating = 1;
while( keep_iterating ) {
(*funcd) (eos, x, dx, resid, jac, &f, &df, n, Bsq,QdotBsq,Qtsq,Qdotn,D); /* returns with new dx, f, df */
/* Save old values before calculating the new: */
errx = 0.;
for( id = 0; id < n ; id++) {
x_old[id] = x[id] ;
}
/* Make the newton step: */
for( id = 0; id < n ; id++) {
x[id] += dx[id] ;
}
/****************************************/
/* Calculate the convergence criterion */
/****************************************/
errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]);
/****************************************/
/* Make sure that the new x[] is physical : */
/****************************************/
validate_x( x, x_old ) ;
/*****************************************************************************/
/* If we've reached the tolerance level, then just do a few extra iterations */
/* before stopping */
/*****************************************************************************/
if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) {
doing_extra = 1;
}
if( doing_extra == 1 ) i_extra++ ;
if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0))
|| (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) {
keep_iterating = 0;
}
n_iter++;
} // END of while(keep_iterating)
/* Check for bad untrapped divergences : */
if( (finite(f)==0) || (finite(df)==0) ) {
return(2);
}
if( fabs(errx) > MIN_NEWT_TOL){
//CCTK_VInfo(CCTK_THORNSTRING,"%d %e %e %e %e",n_iter,f,df,errx,MIN_NEWT_TOL);
return(1);
}
if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){
return(0);
}
if( fabs(errx) <= NEWT_TOL ){
return(0);
}
return(0);
}
Appending to ../src/harm_utoprim_2d.c
%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************/
/*********************************************************************************
func_vsq():
-- calculates the residuals, and Newton step for general_newton_raphson();
-- for this method, x=W,vsq here;
Arguments:
x = current value of independent var's (on input & output);
dx = Newton-Raphson step (on output);
resid = residuals based on x (on output);
jac = Jacobian matrix based on x (on output);
f = resid.resid/2 (on output)
df = -2*f; (on output)
n = dimension of x[];
*********************************************************************************/
static void func_vsq(eos_struct eos, CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
CCTK_REAL jac[][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,
CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D)
{
CCTK_REAL W, vsq, Wsq, p_tmp, dPdvsq, dPdW;
CCTK_REAL t11;
CCTK_REAL t16;
CCTK_REAL t18;
CCTK_REAL t2;
CCTK_REAL t21;
CCTK_REAL t23;
CCTK_REAL t24;
CCTK_REAL t25;
CCTK_REAL t3;
CCTK_REAL t35;
CCTK_REAL t36;
CCTK_REAL t4;
CCTK_REAL t40;
CCTK_REAL t9;
// vv TESTING vv
// CCTK_REAL D,gtmp,gamma,rho0,w,p,u;
// ^^ TESTING ^^
W = x[0];
vsq = x[1];
Wsq = W*W;
// vv TESTING vv
/*
D = U[RHO] ;
gtmp = sqrt(1. - vsq);
gamma = 1./gtmp ;
rho0 = D * gtmp;
w = W * (1. - vsq) ;
p = pressure_rho0_w(eos, rho0,w) ;
u = w - (rho0 + p) ;
if(u<=0 && 1==1) {
vsq = 0.9999999 * (1.0-(rho0+p)/W);
w = W * (1. - vsq) ;
p = pressure_rho0_w(eos, rho0,w) ;
u = w - (rho0 + p) ;
//CCTK_VInfo(CCTK_THORNSTRING,"%e check",u);
}
*/
// ^^ TESTING ^^
p_tmp = pressure_W_vsq( eos, W, vsq , D);
dPdW = dpdW_calc_vsq( W, vsq );
dPdvsq = dpdvsq_calc( eos, W, vsq, D );
// These expressions were calculated using Mathematica, but made into efficient
// code using Maple. Since we know the analytic form of the equations, we can
// explicitly calculate the Newton-Raphson step:
t2 = -0.5*Bsq+dPdvsq;
t3 = Bsq+W;
t4 = t3*t3;
t9 = 1/Wsq;
t11 = Qtsq-vsq*t4+QdotBsq*(Bsq+2.0*W)*t9;
t16 = QdotBsq*t9;
t18 = -Qdotn-0.5*Bsq*(1.0+vsq)+0.5*t16-W+p_tmp;
t21 = 1/t3;
t23 = 1/W;
t24 = t16*t23;
t25 = -1.0+dPdW-t24;
t35 = t25*t3+(Bsq-2.0*dPdvsq)*(QdotBsq+vsq*Wsq*W)*t9*t23;
t36 = 1/t35;
dx[0] = -(t2*t11+t4*t18)*t21*t36;
t40 = (vsq+t24)*t3;
dx[1] = -(-t25*t11-2.0*t40*t18)*t21*t36;
//detJ = t3*t35; // <- set but not used...
jac[0][0] = -2.0*t40;
jac[0][1] = -t4;
jac[1][0] = t25;
jac[1][1] = t2;
resid[0] = t11;
resid[1] = t18;
*df = -resid[0]*resid[0] - resid[1]*resid[1];
*f = -0.5 * ( *df );
}
Appending to ../src/harm_utoprim_2d.c
%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************
**********************************************************************
The following routines specify the equation of state. All routines
above here should be indpendent of EOS. If the user wishes
to use another equation of state, the below functions must be replaced
by equivalent routines based upon the new EOS.
**********************************************************************
**********************************************************************/
Appending to ../src/harm_utoprim_2d.c
pressure_W_vsq()
function [Back to top]¶This function computes $p\left(W,v^{2}\right)$. For a $\Gamma$-law equation of state,
$$ p_{\Gamma} = \left(\Gamma-1\right)u\ , $$and with the definitions
$$ \begin{align} \gamma^{2} &= \frac{1}{1-v^{2}}\ ,\\ W &= \gamma^{2}w\ ,\\ D &= \gamma\rho_{b}\ ,\\ w &= \rho_{b} + u + p\ , \end{align} $$we have
$$ \begin{align} p_{\Gamma} &= \left(\Gamma-1\right)u\\ &= \left(\Gamma-1\right)\left(w - \rho_{b} - p_{\Gamma}\right)\\ &= \left(\Gamma-1\right)\left(\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\right) - \left(\Gamma-1\right)p_{\Gamma}\\ \implies &\boxed{ p_{\Gamma} = \frac{\left(\Gamma-1\right)}{\Gamma}\left(\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\right) }\ . \end{align} $$Thus, the pre-PPEOS Patch version of this function was
/**********************************************************************/
/**********************************************************************
pressure_W_vsq():
-- Gamma-law equation of state;
-- pressure as a function of W, vsq, and D:
**********************************************************************/
static CCTK_REAL pressure_W_vsq(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)
{
CCTK_REAL gtmp;
gtmp = 1. - vsq;
return( (GAMMA - 1.) * ( W * gtmp - D * sqrt(gtmp) ) / GAMMA );
}
We are now, however, interested in the hybrid EOS of the form
$$ p_{\rm hybrid} = P_{\rm cold} + P_{\rm th}\ , $$where $P_{\rm cold}$ is given by a single or piecewise polytrope EOS,
$$ P_{\rm cold} = K_{i}\rho_{b}^{\Gamma_{i}}\ , $$$P_{\rm th}$ accounts for thermal effects and is given by
$$ P_{\rm th} = \left(\Gamma_{\rm th} - 1\right)\epsilon_{\rm th}\ , $$and
$$ \begin{align} \epsilon \equiv \frac{u}{\rho_{b}} &= \epsilon_{\rm th}+\epsilon_{\rm cold}\ ,\\ \epsilon_{\rm cold} &= \int d\rho \frac{P_{\rm cold}(\rho)}{\rho^{2}}\ . \end{align} $$We then have
$$ \begin{align} p_{\rm hybrid} &= P_{\rm cold} + P_{\rm th}\\ &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\rho_{b}\epsilon_{\rm th}\\ &= 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 - \frac{D}{\gamma}\epsilon_{\rm cold}\right)\\ &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left(w - \rho_{b} - p_{\rm hybrid} - \frac{D}{\gamma}\epsilon_{\rm cold}\right)\\ &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left(\frac{W}{\gamma^{2}} - \frac{D}{\gamma} - \frac{D}{\gamma}\epsilon_{\rm cold}\right)-\left(\Gamma_{\rm th}-1\right)p_{\rm hybrid}\\ &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left[\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\left(1+\epsilon_{\rm cold}\right)\right]-\left(\Gamma_{\rm th}-1\right)p_{\rm hybrid}\\ \implies &\boxed{ p_{\rm hybrid} = \frac{P_{\rm cold}}{\Gamma_{\rm th}} + \frac{\left(\Gamma_{\rm th}-1\right)}{\Gamma_{\rm th}}\left[\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\left(1+\epsilon_{\rm cold}\right)\right] } \end{align} $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************/
/**********************************************************************
pressure_W_vsq():
-- Hybrid single and piecewise polytropic equation of state;
-- pressure as a function of P_cold, eps_cold, W, vsq, and D:
**********************************************************************/
static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)
{
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
// Compute gamma^{-2} = 1 - v^{2} and gamma^{-1}
CCTK_REAL inv_gammasq = 1.0 - vsq;
CCTK_REAL inv_gamma = sqrt(inv_gammasq);
// Compute rho_b = D / gamma
CCTK_REAL rho_b = D*inv_gamma;
// Compute P_cold and eps_cold
CCTK_REAL P_cold, eps_cold;
compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold);
// Compute p = P_{cold} + P_{th}
return( ( P_cold + (Gamma_th - 1.0)*( W*inv_gammasq - D*inv_gamma*( 1.0 + eps_cold ) ) )/Gamma_th );
}
Appending to ../src/harm_utoprim_2d.c
dpdW_calc_vsq()
function [Back to top]¶This function computes $\frac{\partial p\left(W,v^{2}\right)}{\partial W}$. For a $\Gamma$-law equation of state, remember that
$$ p_{\Gamma} = \frac{\left(\Gamma-1\right)}{\Gamma}\left(\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\right)\ , $$which then implies
$$ \boxed{\frac{\partial p_{\Gamma}}{\partial W} = \frac{\Gamma-1}{\Gamma \gamma^{2}} = \frac{\left(\Gamma-1\right)\left(1-v^{2}\right)}{\Gamma}}\ . $$Thus, the pre-PPEOS Patch version of this function was
/**********************************************************************/
/**********************************************************************
dpdW_calc_vsq():
-- partial derivative of pressure with respect to W;
**********************************************************************/
static CCTK_REAL dpdW_calc_vsq(CCTK_REAL W, CCTK_REAL vsq)
{
return( (GAMMA - 1.) * (1. - vsq) / GAMMA ) ;
}
For the case of a hybrid, single or piecewise polytropic EOS, we have
$$ p_{\rm hybrid} = \frac{P_{\rm cold}}{\Gamma_{\rm th}} + \frac{\left(\Gamma_{\rm th}-1\right)}{\Gamma_{\rm th}}\left[\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\left(1+\epsilon_{\rm cold}\right)\right]\ . $$It is important to notice that the cold components of $p_{\rm hybrid}$ are not functions of $W$, but instead functions of $D$: $P_{\rm cold} = P_{\rm cold}(\rho_{b}) = P_{\rm cold}(D)$ and $\epsilon_{\rm cold} = \epsilon_{\rm cold}(\rho_{b}) = \epsilon_{\rm cold}(D)$. Thus
$$ \boxed{\frac{\partial p_{\rm hybrid}}{\partial W} = \frac{\Gamma_{\rm th}-1}{\Gamma_{\rm th} \gamma^{2}} = \frac{\left(\Gamma_{\rm th}-1\right)\left(1-v^{2}\right)}{\Gamma_{\rm th}}}\ . $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************/
/**********************************************************************
dpdW_calc_vsq():
-- partial derivative of pressure with respect to W;
**********************************************************************/
static CCTK_REAL dpdW_calc_vsq(CCTK_REAL W, CCTK_REAL vsq)
{
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
return( (Gamma_th - 1.0) * (1.0 - vsq) / Gamma_th ) ;
}
Appending to ../src/harm_utoprim_2d.c
dpdvsq_calc()
function [Back to top]¶This function computes $\frac{\partial p\left(W,v^{2}\right)}{\partial W}$. For a $\Gamma$-law equation of state, remember that
$$ p_{\Gamma} = \frac{\left(\Gamma-1\right)}{\Gamma}\left(\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\right) = \frac{\left(\Gamma-1\right)}{\Gamma}\left[W\left(1-v^{2}\right) - D\sqrt{1-v^{2}}\right]\ , $$which then implies
$$ \boxed{\frac{\partial p_{\Gamma}}{\partial\left(v^{2}\right)} = \frac{\Gamma-1}{\Gamma}\left(\frac{D}{2\sqrt{1-v^{2}}}-W\right)} \ . $$Thus, the pre-PPEOS Patch version of this function was
/**********************************************************************/
/**********************************************************************
dpdvsq_calc():
-- partial derivative of pressure with respect to vsq
**********************************************************************/
static CCTK_REAL dpdvsq_calc(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)
{
return( (GAMMA - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / GAMMA ) ;
}
For the case of a hybrid, single or piecewise polytropic EOS, we have
$$ p_{\rm hybrid} = \frac{P_{\rm cold}}{\Gamma_{\rm th}} + \frac{\left(\Gamma_{\rm th}-1\right)}{\Gamma_{\rm th}}\left[\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\left(1+\epsilon_{\rm cold}\right)\right]\ . $$Let us thus begin by setting the necessary parameters from the hybrid EOS.
%%writefile -a $outfile_path__harm_utoprim_2d__c
/**********************************************************************/
/**********************************************************************
dpdvsq_calc():
-- partial derivative of pressure with respect to vsq
**********************************************************************/
static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)
{
// This sets Gamma_th
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
// Set gamma and rho
CCTK_REAL gamma = 1.0/sqrt(1.0 - vsq);
CCTK_REAL rho_b = D/gamma;
// Compute P_cold and eps_cold
CCTK_REAL P_cold, eps_cold;
compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold);
// Set basic polytropic quantities
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b);
CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];
Appending to ../src/harm_utoprim_2d.c
Next, remember that $P_{\rm cold} = P_{\rm cold}(\rho_{b}) = P_{\rm cold}(D,v^{2})$ and also $\epsilon_{\rm cold} = \epsilon_{\rm cold}(D,v^{2})$. Therefore, we must start by finding the derivatives of $P_{\rm cold}$ and $\epsilon_{\rm cold}$ with respect to $v^{2}$.
Let us first notice that
$$ \frac{\partial\gamma}{\partial\left(v^{2}\right)} = \frac{\partial}{\partial\left(v^{2}\right)}\left[\frac{1}{\sqrt{1-v^{2}}}\right] = \frac{1}{2}\left(1-v^{2}\right)^{-3/2} = \frac{\gamma^{3}}{2}\ . $$Thus, for a general power
$$ \frac{\partial\gamma^{a}}{\partial\left(v^{2}\right)} = a\gamma^{a-1}\frac{\partial\gamma}{\partial\left(v^{2}\right)} = a\gamma^{a-1}\left(\frac{\gamma^{3}}{2}\right) = \frac{a}{2}\gamma^{a+2} $$Thus we have
$$ \begin{align} \frac{\partial P_{\rm cold}}{\partial \left(v^{2}\right)} &= \frac{\partial}{\partial\left(v^{2}\right)}\left(K_{\rm poly}\rho_{b}^{\Gamma_{\rm poly}}\right)\\ &= \frac{\partial}{\partial\left(v^{2}\right)}\left[K_{\rm poly}\left(\frac{D}{\gamma}\right)^{\Gamma_{\rm poly}}\right]\\ &= K_{\rm poly}D^{\Gamma_{\rm poly}}\frac{\partial}{\partial\left(v^{2}\right)}\left[\gamma^{-\Gamma_{\rm poly}/2}\right]\\ &=K_{\rm poly}D^{\Gamma_{\rm poly}}\left[\frac{-\Gamma_{\rm poly}/2}{2}\gamma^{-\Gamma_{\rm poly}/2 + 2}\right]\\ &=K_{\rm poly}\left(\frac{D}{\gamma}\right)^{\Gamma_{\rm poly}}\gamma^{-\frac{\Gamma_{\rm poly}}{2} + 2 + \Gamma_{\rm poly}}\\ \implies &\boxed{ \frac{\partial P_{\rm cold}}{\partial \left(v^{2}\right)} = \gamma^{2+\frac{\Gamma_{\rm poly}}{2}}P_{\rm cold}}\ . \end{align} $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/* Now we implement the derivative of P_cold with respect
* to v^{2}, given by
* ----------------------------------------------------
* | dP_cold/dvsq = gamma^{2 + Gamma_{poly}/2} P_{cold} |
* ----------------------------------------------------
*/
CCTK_REAL dPcold_dvsq = P_cold * pow(gamma,2.0 + 0.5*Gamma_ppoly_tab);
Appending to ../src/harm_utoprim_2d.c
Now, obtaining $\epsilon_{\rm cold}$ from $P_{\rm cold}$ requires an integration and, therefore, generates an integration constant. Since we are interested in a derivative of $\epsilon_{\rm cold}$, however, we will simply drop the constant altogether. Remember that:
$$ \epsilon_{\rm cold} = K_{\rm poly}\int d\rho_{b} \rho_{b}^{\Gamma_{\rm poly}-2} = \frac{K_{\rm poly}\rho_{b}^{\Gamma_{\rm poly}-1}}{\Gamma_{\rm poly}-1} = \frac{P_{\rm cold}}{\rho_{b}\left(\Gamma_{\rm poly}-1\right)} = \frac{\gamma P_{\rm cold}}{D\left(\Gamma_{\rm poly}-1\right)}\ . $$Thus
$$ \begin{align} \frac{\partial \epsilon_{\rm cold}}{\partial \left(v^{2}\right)} &= \frac{1}{D\left(\Gamma_{\rm poly}-1\right)}\left[\gamma\frac{\partial P_{\rm cold}}{\partial \left(v^{2}\right)} + P_{\rm cold}\frac{\partial\gamma}{\partial \left(v^{2}\right)}\right]\\ &=\frac{1}{D\left(\Gamma_{\rm poly}-1\right)}\left[\gamma\frac{\partial P_{\rm cold}}{\partial \left(v^{2}\right)} + P_{\rm cold}\left(\frac{\gamma^{3}}{2}\right)\right]\\ \implies &\boxed{ \frac{\partial \epsilon_{\rm cold}}{\partial \left(v^{2}\right)} = \frac{\gamma}{D\left(\Gamma_{\rm poly}-1\right)}\left[\frac{\partial P_{\rm cold}}{\partial \left(v^{2}\right)} + \frac{\gamma^{2} P_{\rm cold}}{2}\right]\ . } \end{align} $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/* Now we implement the derivative of eps_cold with respect
* to v^{2}, given by
* -----------------------------------------------------------------------------------
* | deps_cold/dvsq = gamma/(D*(Gamma_ppoly_tab-1)) * (dP_cold/dvsq + gamma^{2} P_cold / 2) |
* -----------------------------------------------------------------------------------
*/
CCTK_REAL depscold_dvsq = ( gamma/(D*(Gamma_ppoly_tab-1.0)) ) * ( dPcold_dvsq + 0.5*gamma*gamma*P_cold );
Appending to ../src/harm_utoprim_2d.c
Finally, remembering that
$$ \begin{align} p_{\rm hybrid} &= \frac{P_{\rm cold}}{\Gamma_{\rm th}} + \frac{\left(\Gamma_{\rm th}-1\right)}{\Gamma_{\rm th}}\left[\frac{W}{\gamma^{2}} - \frac{D}{\gamma}\left(1+\epsilon_{\rm cold}\right)\right]\ ,\\ \frac{\partial\gamma^{a}}{\partial\left(v^{2}\right)} &= \frac{a}{2}\gamma^{a+2}\ , \end{align} $$we have
$$ \boxed{ \frac{\partial p_{\rm hybrid}}{\partial\left(v^{2}\right)} = \frac{1}{\Gamma_{\rm th}}\left\{\frac{\partial P_{\rm cold}}{\partial\left(v^{2}\right)} + \left(\Gamma_{\rm th}-1\right)\left[-W + \frac{D\gamma}{2}\left(1+\epsilon_{\rm cold}\right) - \frac{D}{\gamma}\frac{\partial \epsilon_{\rm cold}}{\partial\left(v^{2}\right)}\right]\right\}\ . } $$%%writefile -a $outfile_path__harm_utoprim_2d__c
/* Now we implement the derivative of p_hybrid with respect
* to v^{2}, given by
* -----------------------------------------------------------------------------
* | dp/dvsq = Gamma_th^{-1}( dP_cold/dvsq |
* | + (Gamma_{th}-1)*(-W |
* | + D gamma (1 + eps_cold)/2 |
* | - (D/gamma) * deps_cold/dvsq) ) |
* -----------------------------------------------------------------------------
*/
return( ( dPcold_dvsq + (Gamma_th-1.0)*( -W + D*gamma*(1+eps_cold)/2.0 - D*depscold_dvsq/gamma ) )/Gamma_th );
}
/******************************************************************************
END OF UTOPRIM_2D.C
******************************************************************************/
#endif
Appending to ../src/harm_utoprim_2d.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_utoprim_2d.c"
original_IGM_file_name = "harm_utoprim_2d-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_utoprim_2d__c = !diff $original_IGM_file_path $outfile_path__harm_utoprim_2d__c
if Validation__harm_utoprim_2d__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_utoprim_2d.c: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for harm_utoprim_2d.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_utoprim_2d__c:
print(diff_line)
Validation test for harm_utoprim_2d.c: FAILED! Diff: 0a1,2 > #ifndef __HARM_UTOPRIM_2D__C__ > #define __HARM_UTOPRIM_2D__C__ 2c4 < Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble, --- > Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble, 7c9 < This file is part of HARM. HARM is a program that solves hyperbolic --- > This file is part of HARM. HARM is a program that solves hyperbolic 9,10c11,12 < shock-capturing techniques. This version of HARM has been configured to < solve the relativistic magnetohydrodynamic equations of motion on a --- > shock-capturing techniques. This version of HARM has been configured to > solve the relativistic magnetohydrodynamic equations of motion on a 12c14 < an accretion disk model. --- > an accretion disk model. 14c16 < You are morally obligated to cite the following two papers in his/her --- > You are morally obligated to cite the following two papers in his/her 17c19 < [1] Gammie, C. F., McKinney, J. C., \& Toth, G.\ 2003, --- > [1] Gammie, C. F., McKinney, J. C., \& Toth, G.\ 2003, 20c22 < [2] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, --- > [2] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, 23,24c25,26 < < Further, we strongly encourage you to obtain the latest version of --- > > Further, we strongly encourage you to obtain the latest version of 49c51 < utoprim_2d.c: --- > utoprim_2d.c: 52c54 < Uses the 2D method: --- > Uses the 2D method: 54,55c56,57 < Newton-Raphson method < -- can be used (in principle) with a general equation of state. --- > Newton-Raphson method > -- can be used (in principle) with a general equation of state. 58,60c60,62 < density or internal energy density is calculated. You may want < to change this aspect of the code so that it still calculates the < velocity and so that you can floor the densities. If you want to --- > density or internal energy density is calculated. You may want > to change this aspect of the code so that it still calculates the > velocity and so that you can floor the densities. If you want to 68c70 < // Declarations: --- > // Declarations: 70,72c72,74 < static int Utoprim_new_body(CCTK_REAL U[], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[],long &n_iter); < static int general_newton_raphson( CCTK_REAL x[], int n, long &n_iter, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *, CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D); < static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D); --- > static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[],long &n_iter); > static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter, void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *, CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D); > static void func_vsq( eos_struct eos, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D); 74c76 < static CCTK_REAL pressure_W_vsq(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) ; --- > static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) ; 76c78 < static CCTK_REAL dpdvsq_calc(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D); --- > static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D); 82c84 < --- > 84c86 < between the two sets of definitions for U and P. The user may --- > between the two sets of definitions for U and P. The user may 89,97c91,99 < / rho u^t \ < U = | T^t_t + rho u^t | sqrt(-det(g_{\mu\nu})) < | T^t_i | < \ B^i / < < / rho \ < P = | uu | < | \tilde{u}^i | < \ B^i / --- > / rho u^t \ > U = | T^t_t + rho u^t | sqrt(-det(g_{\mu\nu})) > | T^t_i | > \ B^i / > > / rho \ > P = | uu | > | \tilde{u}^i | > \ B^i / 101c103 < U[NPR] = conserved variables (current values on input/output); --- > U[NPR] = conserved variables (current values on input/output); 105,107c107,109 < prim[NPR] = primitive variables (guess on input, calculated values on < output if there are no problems); < --- > prim[NPR] = primitive variables (guess on input, calculated values on > output if there are no problems); > 109c111 < unfamiliar with metrics, merely set --- > unfamiliar with metrics, merely set 114c116,117 < int Utoprim_2d(CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], --- > > int Utoprim_2d(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], 119c122 < int i, ret; --- > int i, ret; 122c125 < if( U[0] <= 0. ) { --- > if( U[0] <= 0. ) { 131c134,135 < --- > > 141a146 > 150c155 < ret = Utoprim_new_body(U_tmp, gcov, gcon, gdet, prim_tmp,n_iter); --- > ret = Utoprim_new_body(eos, U_tmp, gcov, gcon, gdet, prim_tmp,n_iter); 152c157 < /* Transform new primitive variables back if there was no problem : */ --- > /* Transform new primitive variables back if there was no problem : */ 163a169 > 171,172c177,178 < -- This is the main routine that calculates auxiliary quantities for the < Newton-Raphson routine. --- > -- This is the main routine that calculates auxiliary quantities for the > Newton-Raphson routine. 174,177c180,183 < -- assumes that < / rho gamma \ < U = | alpha T^t_\mu | < \ alpha B^i / --- > -- assumes that > / rho gamma \ > U = | alpha T^t_\mu | > \ alpha B^i / 181,184c187,190 < / rho \ < prim = | uu | < | \tilde{u}^i | < \ alpha B^i / --- > / rho \ > prim = | uu | > | \tilde{u}^i | > \ alpha B^i / 187,188c193,194 < return: (i*100 + j) where < i = 0 -> Newton-Raphson solver either was not called (yet or not used) --- > return: (i*100 + j) where > i = 0 -> Newton-Raphson solver either was not called (yet or not used) 190c196 < 1 -> Newton-Raphson solver did not converge to a solution with the --- > 1 -> Newton-Raphson solver did not converge to a solution with the 192c198 < 2 -> Newton-Raphson procedure encountered a numerical divergence --- > 2 -> Newton-Raphson procedure encountered a numerical divergence 194,196c200,202 < < j = 0 -> success < 1 -> failure: some sort of failure in Newton-Raphson; --- > > j = 0 -> success > 1 -> failure: some sort of failure in Newton-Raphson; 198,199c204,205 < 3 -> failure: W<0 or W>W_TOO_BIG < 4 -> failure: v^2 > 1 --- > 3 -> failure: W<0 or W>W_TOO_BIG > 4 -> failure: v^2 > 1 204c210 < static int Utoprim_new_body(CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], --- > static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], 217a224 > 237a245,248 > // FIXME: The exact form of n^{\mu} can be found > // in eq. (2.116) and implementing it > // directly is a lot more efficient than > // performing n^{\mu} = g^{\mu\nu}n_{nu} 248a260 > 255c267 < if( (utsq < 0.) && (fabs(utsq) < 1.0e-13) ) { --- > if( (utsq < 0.) && (fabs(utsq) < 1.0e-13) ) { 265c277 < --- > 267c279 < // i.e. you don't get positive values for dP/d(vsq) . --- > // i.e. you don't get positive values for dP/d(vsq) . 270c282 < p = pressure_rho0_u(rho0,u) ; --- > p = pressure_rho0_u(eos, rho0,u) ; 276c288 < // Make sure that W is large enough so that v^2 < 1 : --- > // Make sure that W is large enough so that v^2 < 1 : 278c290 < while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) --- > while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) 284,285c296,298 < < // Calculate W and vsq: --- > > > // Calculate W and vsq: 288c301 < retval = general_newton_raphson( x_2d, n, n_iter, func_vsq, Bsq,QdotBsq,Qtsq,Qdotn,D) ; --- > retval = general_newton_raphson( eos, x_2d, n, n_iter, func_vsq, Bsq,QdotBsq,Qtsq,Qdotn,D) ; 292c305 < --- > 311a325 > 318c332 < p = pressure_rho0_w(rho0,w) ; --- > p = pressure_rho0_w(eos, rho0,w) ; 322c336 < // User may want to handle this case differently, e.g. do NOT return upon --- > // User may want to handle this case differently, e.g. do NOT return upon 342c356 < --- > 353c367,368 < /**********************************************************************/ --- > > /**********************************************************************/ 355,358c370,373 < vsq_calc(): < < -- evaluate v^2 (spatial, normalized velocity) from < W = \gamma^2 w --- > vsq_calc(): > > -- evaluate v^2 (spatial, normalized velocity) from > W = \gamma^2 w 364c379 < --- > 371a387 > 374,375c390,391 < x1_of_x0(): < --- > x1_of_x0(): > 382c398 < static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D ) --- > static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D ) 386,387d401 < < vsq = fabs(vsq_calc(x0,Bsq,QdotBsq,Qtsq,Qdotn,D)) ; // guaranteed to be positive 388a403 > vsq = fabs(vsq_calc(x0,Bsq,QdotBsq,Qtsq,Qdotn,D)) ; // guaranteed to be positive 390c405,406 < return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); --- > > return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); 393a410 > 396,398c413,415 < validate_x(): < < -- makes sure that x[0,1] have physical values, based upon --- > validate_x(): > > -- makes sure that x[0,1] have physical values, based upon 400c417 < --- > 403c420 < static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) --- > static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) 405c422 < --- > 408c425 < /* Always take the absolute value of x[0] and check to see if it's too big: */ --- > /* Always take the absolute value of x[0] and check to see if it's too big: */ 411c428 < --- > 419a437 > 422c440 < general_newton_raphson(): --- > general_newton_raphson(): 429,431c447,449 < static int general_newton_raphson( CCTK_REAL x[], int n, long &n_iter, < void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], < CCTK_REAL [][NEWT_DIM], CCTK_REAL *, --- > static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter, > void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], > CCTK_REAL [][NEWT_DIM], CCTK_REAL *, 443c461 < errx = 1. ; --- > errx = 1. ; 452c470,472 < while( keep_iterating ) { --- > while( keep_iterating ) { > > (*funcd) (eos, x, dx, resid, jac, &f, &df, n, Bsq,QdotBsq,Qtsq,Qdotn,D); /* returns with new dx, f, df */ 454,455d473 < (*funcd) (x, dx, resid, jac, &f, &df, n, Bsq,QdotBsq,Qtsq,Qdotn,D); /* returns with new dx, f, df */ < 484c502 < --- > 491c509 < if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) --- > if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) 509c527 < } --- > } 522a541 > 525c544 < func_vsq(): --- > func_vsq(): 540c559 < static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], --- > static void func_vsq(eos_struct eos, CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 545c564 < --- > 568c587 < --- > 579c598 < p = pressure_rho0_w(rho0,w) ; --- > p = pressure_rho0_w(eos, rho0,w) ; 586c605 < p = pressure_rho0_w(rho0,w) ; --- > p = pressure_rho0_w(eos, rho0,w) ; 594,595c613,614 < < p_tmp = pressure_W_vsq( W, vsq , D); --- > > p_tmp = pressure_W_vsq( eos, W, vsq , D); 597c616 < dPdvsq = dpdvsq_calc( W, vsq, D ); --- > dPdvsq = dpdvsq_calc( eos, W, vsq, D ); 599,601c618,620 < // These expressions were calculated using Mathematica, but made into efficient < // code using Maple. Since we know the analytic form of the equations, we can < // explicitly calculate the Newton-Raphson step: --- > // These expressions were calculated using Mathematica, but made into efficient > // code using Maple. Since we know the analytic form of the equations, we can > // explicitly calculate the Newton-Raphson step: 627c646 < --- > 636,642c655,662 < /********************************************************************** < ********************************************************************** < < The following routines specify the equation of state. All routines < above here should be indpendent of EOS. If the user wishes < to use another equation of state, the below functions must be replaced < by equivalent routines based upon the new EOS. --- > > /********************************************************************** > ********************************************************************** > > The following routines specify the equation of state. All routines > above here should be indpendent of EOS. If the user wishes > to use another equation of state, the below functions must be replaced > by equivalent routines based upon the new EOS. 646a667 > 648,652c669,673 < /********************************************************************** < pressure_W_vsq(): < < -- Gamma-law equation of state; < -- pressure as a function of W, vsq, and D: --- > /********************************************************************** > pressure_W_vsq(): > > -- Hybrid single and piecewise polytropic equation of state; > -- pressure as a function of P_cold, eps_cold, W, vsq, and D: 654c675 < static CCTK_REAL pressure_W_vsq(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) --- > static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) 655a677,678 > > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 656a680,691 > #endif > > // Compute gamma^{-2} = 1 - v^{2} and gamma^{-1} > CCTK_REAL inv_gammasq = 1.0 - vsq; > CCTK_REAL inv_gamma = sqrt(inv_gammasq); > > // Compute rho_b = D / gamma > CCTK_REAL rho_b = D*inv_gamma; > > // Compute P_cold and eps_cold > CCTK_REAL P_cold, eps_cold; > compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold); 658,661c693,694 < CCTK_REAL gtmp; < gtmp = 1. - vsq; < < return( (gamma_th /* <- Should be local polytropic Gamma factor */ - 1.) * ( W * gtmp - D * sqrt(gtmp) ) / gamma_th /* <- Should be local polytropic Gamma factor */ ); --- > // Compute p = P_{cold} + P_{th} > return( ( P_cold + (Gamma_th - 1.0)*( W*inv_gammasq - D*inv_gamma*( 1.0 + eps_cold ) ) )/Gamma_th ); 665a699 > 667,669c701,703 < /********************************************************************** < dpdW_calc_vsq(): < --- > /********************************************************************** > dpdW_calc_vsq(): > 673a708,709 > > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 675c711,713 < return( (gamma_th /* <- Should be local polytropic Gamma factor */ - 1.) * (1. - vsq) / gamma_th /* <- Should be local polytropic Gamma factor */ ) ; --- > #endif > > return( (Gamma_th - 1.0) * (1.0 - vsq) / Gamma_th ) ; 678a717 > 680,682c719,721 < /********************************************************************** < dpdvsq_calc(): < --- > /********************************************************************** > dpdvsq_calc(): > 685c724 < static CCTK_REAL dpdvsq_calc(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) --- > static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) 686a726,728 > > // This sets Gamma_th > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 688c730,772 < return( (gamma_th /* <- Should be local polytropic Gamma factor */ - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gamma_th /* <- Should be local polytropic Gamma factor */ ) ; --- > #endif > > > // Set gamma and rho > CCTK_REAL gamma = 1.0/sqrt(1.0 - vsq); > CCTK_REAL rho_b = D/gamma; > > // Compute P_cold and eps_cold > CCTK_REAL P_cold, eps_cold; > compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold); > > // Set basic polytropic quantities > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b); > CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index]; > > > /* Now we implement the derivative of P_cold with respect > * to v^{2}, given by > * ---------------------------------------------------- > * | dP_cold/dvsq = gamma^{2 + Gamma_{poly}/2} P_{cold} | > * ---------------------------------------------------- > */ > CCTK_REAL dPcold_dvsq = P_cold * pow(gamma,2.0 + 0.5*Gamma_ppoly_tab); > > > /* Now we implement the derivative of eps_cold with respect > * to v^{2}, given by > * ----------------------------------------------------------------------------------- > * | deps_cold/dvsq = gamma/(D*(Gamma_ppoly_tab-1)) * (dP_cold/dvsq + gamma^{2} P_cold / 2) | > * ----------------------------------------------------------------------------------- > */ > CCTK_REAL depscold_dvsq = ( gamma/(D*(Gamma_ppoly_tab-1.0)) ) * ( dPcold_dvsq + 0.5*gamma*gamma*P_cold ); > > /* Now we implement the derivative of p_hybrid with respect > * to v^{2}, given by > * ----------------------------------------------------------------------------- > * | dp/dvsq = Gamma_th^{-1}( dP_cold/dvsq | > * | + (Gamma_{th}-1)*(-W | > * | + D gamma (1 + eps_cold)/2 | > * | - (D/gamma) * deps_cold/dvsq) ) | > * ----------------------------------------------------------------------------- > */ > return( ( dPcold_dvsq + (Gamma_th-1.0)*( -W + D*gamma*(1+eps_cold)/2.0 - D*depscold_dvsq/gamma ) )/Gamma_th ); 692c776 < /****************************************************************************** --- > /****************************************************************************** 694a779 > #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_utoprim_2d.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_utoprim_2d.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_utoprim_2d.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_utoprim_2d.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_utoprim_2d.tex
!rm -f Tut*.out Tut*.aux Tut*.log