This module is currently under development
This module is organized as follows
Step 0: Source directory creation
Step 1: Introduction
Step 2: driver_conserv_to_prims.C
Step 2.a: Converting ADM quantities to BSSN quantities and enforcing $\bar\gamma=1$
Step 2.b: Applying equatorial symmetry
Step 2.c: Setting up the variables needed by HARM
Step 2.d: Determining the primitives variables from the conservatives variables
Step 2.e: Enforcing physical limits on primitives and recomputing the conservatives variables
Step 2.f: Updating conservative and primitive gridfunctions
Step 2.g: Diagnostics and debugging tools
Step 3: harm_primitives_lowlevel.C
HARM
HARM
conservative-to-primitive solverStep 4: font_fix_hybrid_EOS.C
Step 5: harm_primitives_headers.h
Step 6: Code validation
Step 7: Output this notebook to $\LaTeX$-formatted PDF file
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__driver_conserv_to_prims__C = os.path.join(IGM_src_dir_path,"driver_conserv_to_prims.C")
outfile_path__harm_primitives_lowlevel__C = os.path.join(IGM_src_dir_path,"harm_primitives_lowlevel.C")
outfile_path__font_fix_hybrid_EOS__C = os.path.join(IGM_src_dir_path,"font_fix_hybrid_EOS.C")
outfile_path__harm_primitives_headers__h = os.path.join(IGM_src_dir_path,"harm_primitives_headers.h")
driver_conserv_to_prims.C
[Back to top]¶We start here by creating the driver_conserv_to_prims.C
file and loading all files used by it. Note that of the files loaded, we have the following IllinoisGRMHD
files:
harm_primitives_headers.h
: we will discuss this file in step 4 of this tutorial module.inlined_functions.C
: this file is discussed in the inlined_functions NRPy tutorial moduleapply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C
: this file is discussed in the apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs NRPy tutorial module%%writefile $outfile_path__driver_conserv_to_prims__C
/* We evolve forward in time a set of functions called the
* "conservative variables", and any time the conserv's
* are updated, we must solve for the primitive variables
* (rho, pressure, velocities) using a Newton-Raphson
* technique, before reconstructing & evaluating the RHSs
* of the MHD equations again.
*
* This file contains the driver routine for this Newton-
* Raphson solver. Truncation errors in conservative
* variables can lead to no physical solutions in
* primitive variables. We correct for these errors here
* through a number of tricks described in the appendices
* of http://arxiv.org/pdf/1112.0568.pdf.
*
* This is a wrapper for the 2d solver of Noble et al. See
* harm_utoprim_2d.c for references and copyright notice
* for that solver. This wrapper was primarily written by
* Zachariah Etienne & Yuk Tung Liu, in 2011-2013.
*
* For optimal compatibility, this wrapper is licensed under
* the GPL v2 or any later version.
*
* Note that this code assumes a simple gamma law for the
* moment, though it would be easy to extend to a piecewise
* polytrope. */
// Standard #include's
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <ctime>
#include <cstdlib>
#ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER
#include "standalone_conserv_to_prims_main_function.h"
#else
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"
#include "IllinoisGRMHD_headers.h"
#include "harm_primitives_headers.h"
#include "harm_u2p_util.c"
#include "inlined_functions.C"
#include "apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C"
extern "C" void IllinoisGRMHD_conserv_to_prims(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
// We use proper C++ here, for file I/O later.
using namespace std;
#endif
/**********************************
* Piecewise Polytropic EOS Patch *
* Setting up the EOS struct *
**********************************/
/*
* The short piece of code below takes care
* of initializing the EOS parameters.
* Please refer to the "inlined_functions.C"
* source file for the documentation on the
* function.
*/
eos_struct eos;
initialize_EOS_struct_from_input(eos);
Writing ../src/driver_conserv_to_prims.C
Here we compute the conformal metric, $\bar\gamma_{ij}$, and its inverse, $\bar\gamma^{ij}$, from the physical metric $\gamma_{ij}$. We also compute $\phi$, the conformal factor, and $\psi\equiv e^{\phi}$. Finally, we enforce the constraint $\bar\gamma = \det\left(\bar\gamma_{ij}\right) = 1$. The entire procedure is explained in detail in the convert ADM to BSSN and enforce determinant constraint NRPy tutorial module.
%%writefile -a $outfile_path__driver_conserv_to_prims__C
// These BSSN-based variables are not evolved, and so are not defined anywhere that the grid has moved.
// Here we convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.
IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh, gxx,gxy,gxz,gyy,gyz,gzz,alp,
gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
phi_bssn,psi_bssn,lapm1);
Appending to ../src/driver_conserv_to_prims.C
We then use the CardGrid3D ETK thorn to apply equatorial symmetry to our problem.
%%writefile -a $outfile_path__driver_conserv_to_prims__C
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
if(CCTK_EQUALS(Symmetry,"equatorial")) {
// SET SYMMETRY GHOSTZONES ON ALL CONSERVATIVE VARIABLES!
int ierr=0;
ierr+=CartSymGN(cctkGH,"IllinoisGRMHD::grmhd_conservatives");
// FIXME: UGLY. Filling metric ghostzones is needed for, e.g., Cowling runs.
ierr+=CartSymGN(cctkGH,"lapse::lapse_vars");
ierr+=CartSymGN(cctkGH,"bssn::BSSN_vars");
ierr+=CartSymGN(cctkGH,"bssn::BSSN_AH");
ierr+=CartSymGN(cctkGH,"shift::shift_vars");
if(ierr!=0) CCTK_VError(VERR_DEF_PARAMS,"IllinoisGRMHD ERROR (grep for it, foo!) :(");
}
#endif
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
//Start the timer, so we can benchmark the primitives solver during evolution.
// Slower solver -> harder to find roots -> things may be going crazy!
//FIXME: Replace this timing benchmark with something more meaningful, like the avg # of Newton-Raphson iterations per gridpoint!
/*
struct timeval start, end;
long mtime, seconds, useconds;
gettimeofday(&start, NULL);
*/
int failures=0,font_fixes=0,vel_limited_ptcount=0;
int pointcount=0;
int failures_inhoriz=0;
int pointcount_inhoriz=0;
int pressure_cap_hit=0;
CCTK_REAL error_int_numer=0,error_int_denom=0;
int imin=0,jmin=0,kmin=0;
int imax=cctk_lsh[0],jmax=cctk_lsh[1],kmax=cctk_lsh[2];
int rho_star_fix_applied=0;
long n_iter=0;
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
#pragma omp parallel for reduction(+:failures,vel_limited_ptcount,font_fixes,pointcount,failures_inhoriz,pointcount_inhoriz,error_int_numer,error_int_denom,pressure_cap_hit,rho_star_fix_applied,n_iter) schedule(static)
for(int k=kmin;k<kmax;k++)
for(int j=jmin;j<jmax;j++)
for(int i=imin;i<imax;i++) {
int index = CCTK_GFINDEX3D(cctkGH,i,j,k);
int ww;
CCTK_REAL METRIC[NUMVARS_FOR_METRIC],dummy=0;
ww=0;
// FIXME: NECESSARY?
//psi_bssn[index] = exp(phi[index]);
METRIC[ww] = phi_bssn[index];ww++;
METRIC[ww] = dummy; ww++; // Don't need to set psi.
METRIC[ww] = gtxx[index]; ww++;
METRIC[ww] = gtxy[index]; ww++;
METRIC[ww] = gtxz[index]; ww++;
METRIC[ww] = gtyy[index]; ww++;
METRIC[ww] = gtyz[index]; ww++;
METRIC[ww] = gtzz[index]; ww++;
METRIC[ww] = lapm1[index]; ww++;
METRIC[ww] = betax[index]; ww++;
METRIC[ww] = betay[index]; ww++;
METRIC[ww] = betaz[index]; ww++;
METRIC[ww] = gtupxx[index]; ww++;
METRIC[ww] = gtupyy[index]; ww++;
METRIC[ww] = gtupzz[index]; ww++;
METRIC[ww] = gtupxy[index]; ww++;
METRIC[ww] = gtupxz[index]; ww++;
METRIC[ww] = gtupyz[index]; ww++;
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
CCTK_REAL PRIMS[MAXNUMVARS];
ww=0;
PRIMS[ww] = rho_b[index]; ww++;
PRIMS[ww] = P[index]; ww++;
PRIMS[ww] = vx[index]; ww++;
PRIMS[ww] = vy[index]; ww++;
PRIMS[ww] = vz[index]; ww++;
PRIMS[ww] = Bx[index]; ww++;
PRIMS[ww] = By[index]; ww++;
PRIMS[ww] = Bz[index]; ww++;
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
CCTK_REAL CONSERVS[NUM_CONSERVS] = {rho_star[index], mhd_st_x[index],mhd_st_y[index],mhd_st_z[index],tau[index]};
Appending to ../src/driver_conserv_to_prims.C
Then we load the lapse function, $\alpha$, and $\psi$ related variables into the $\rm METRIC\_LAP\_PSI4$ array. Notice that this is done using the ${\rm SET\_LAPSE\_PSI4}()$ "function", which is defined in the IllinoisGRMHD_headers.h
file from IllinoisGRMHD
. The "function" itself is quite simple:
#define SET_LAPSE_PSI4(array_name,METRIC) { \
array_name[LAPSE] = METRIC[LAPM1]+1.0; \
array_name[PSI2] = exp(2.0*METRIC[PHI]); \
array_name[PSI4] = SQR(array_name[PSI2]); \
array_name[PSI6] = array_name[PSI4]*array_name[PSI2]; \
array_name[PSIM4] = 1.0/array_name[PSI4]; \
array_name[LAPSEINV] = 1.0/array_name[LAPSE]; \
}
defining the quantities $\left\{\alpha, \psi^{2}, \psi^{4}, \psi^{6}, \psi^{-4},\alpha^{-1}\right\}$, where $\psi=e^{\phi}$.
%%writefile -a $outfile_path__driver_conserv_to_prims__C
CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX];
SET_LAPSE_PSI4(METRIC_LAP_PSI4,METRIC);
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
CCTK_REAL METRIC_PHYS[NUMVARS_FOR_METRIC];
METRIC_PHYS[GXX] = METRIC[GXX]*METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GXY] = METRIC[GXY]*METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GXZ] = METRIC[GXZ]*METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GYY] = METRIC[GYY]*METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GYZ] = METRIC[GYZ]*METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GZZ] = METRIC[GZZ]*METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GUPXX] = METRIC[GUPXX]*METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPXY] = METRIC[GUPXY]*METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPXZ] = METRIC[GUPXZ]*METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPYY] = METRIC[GUPYY]*METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPYZ] = METRIC[GUPYZ]*METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPZZ] = METRIC[GUPZZ]*METRIC_LAP_PSI4[PSIM4];
Appending to ../src/driver_conserv_to_prims.C
We then evaluate
$$ \beta_{i} = \gamma_{ij}\beta^{j} \implies \boxed{ \left\{ \begin{align} \beta_{x} &= \gamma_{xx}\beta^{x} + \gamma_{xy}\beta^{y} + \gamma_{xz}\beta^{z}\\ \beta_{y} &= \gamma_{yx}\beta^{x} + \gamma_{yy}\beta^{y} + \gamma_{yz}\beta^{z}\\ \beta_{z} &= \gamma_{zx}\beta^{x} + \gamma_{zy}\beta^{y} + \gamma_{zz}\beta^{z} \end{align} \right. }\ , $$and
$$ \boxed{\beta^{2} \equiv \beta_{i}\beta^{i} = \beta_{x}\beta^{x} + \beta_{y}\beta^{y} + \beta_{z}\beta^{z}}\ . $$%%writefile -a $outfile_path__driver_conserv_to_prims__C
CCTK_REAL TUPMUNU[10],TDNMUNU[10];
CCTK_REAL shift_xL = METRIC_PHYS[GXX]*METRIC[SHIFTX] + METRIC_PHYS[GXY]*METRIC[SHIFTY] + METRIC_PHYS[GXZ]*METRIC[SHIFTZ];
CCTK_REAL shift_yL = METRIC_PHYS[GXY]*METRIC[SHIFTX] + METRIC_PHYS[GYY]*METRIC[SHIFTY] + METRIC_PHYS[GYZ]*METRIC[SHIFTZ];
CCTK_REAL shift_zL = METRIC_PHYS[GXZ]*METRIC[SHIFTX] + METRIC_PHYS[GYZ]*METRIC[SHIFTY] + METRIC_PHYS[GZZ]*METRIC[SHIFTZ];
CCTK_REAL beta2L = shift_xL*METRIC[SHIFTX] + shift_yL*METRIC[SHIFTY] + shift_zL*METRIC[SHIFTZ];
Appending to ../src/driver_conserv_to_prims.C
We then setup the ADM 4-metric and its inverse. We refer the reader to eqs. (2.119) and (2.122) from Baumgarte & Shapiro's Numerical Relativity (2010), which are repeated here for the sake of the reader (in reverse order)
$$ \boxed{g_{\mu\nu} = \begin{pmatrix} -\alpha^{2} + \beta_{\ell}\beta^{\ell} & \beta_{i}\\ \beta_{j} & \gamma_{ij} \end{pmatrix}}\ . $$and
$$ \boxed{ g^{\mu\nu} = \begin{pmatrix} -\alpha^{-2} & \alpha^{-2}\beta^{i}\\ \alpha^{-2}\beta^{j} & \gamma^{ij} - \alpha^{-2}\beta^{i}\beta^{j} \end{pmatrix} }\ . $$%%writefile -a $outfile_path__driver_conserv_to_prims__C
// Compute 4-metric, both g_{\mu \nu} and g^{\mu \nu}.
// This is for computing T_{\mu \nu} and T^{\mu \nu}. Also the HARM con2prim lowlevel function requires them.
CCTK_REAL g4dn[4][4],g4up[4][4];
g4dn[0][0] = -SQR(METRIC_LAP_PSI4[LAPSE]) + beta2L;
g4dn[0][1] = g4dn[1][0] = shift_xL;
g4dn[0][2] = g4dn[2][0] = shift_yL;
g4dn[0][3] = g4dn[3][0] = shift_zL;
g4dn[1][1] = METRIC_PHYS[GXX];
g4dn[1][2] = g4dn[2][1] = METRIC_PHYS[GXY];
g4dn[1][3] = g4dn[3][1] = METRIC_PHYS[GXZ];
g4dn[2][2] = METRIC_PHYS[GYY];
g4dn[2][3] = g4dn[3][2] = METRIC_PHYS[GYZ];
g4dn[3][3] = METRIC_PHYS[GZZ];
CCTK_REAL alpha_inv_squared=SQR(METRIC_LAP_PSI4[LAPSEINV]);
g4up[0][0] = -1.0*alpha_inv_squared;
g4up[0][1] = g4up[1][0] = METRIC[SHIFTX]*alpha_inv_squared;
g4up[0][2] = g4up[2][0] = METRIC[SHIFTY]*alpha_inv_squared;
g4up[0][3] = g4up[3][0] = METRIC[SHIFTZ]*alpha_inv_squared;
g4up[1][1] = METRIC_PHYS[GUPXX] - METRIC[SHIFTX]*METRIC[SHIFTX]*alpha_inv_squared;
g4up[1][2] = g4up[2][1] = METRIC_PHYS[GUPXY] - METRIC[SHIFTX]*METRIC[SHIFTY]*alpha_inv_squared;
g4up[1][3] = g4up[3][1] = METRIC_PHYS[GUPXZ] - METRIC[SHIFTX]*METRIC[SHIFTZ]*alpha_inv_squared;
g4up[2][2] = METRIC_PHYS[GUPYY] - METRIC[SHIFTY]*METRIC[SHIFTY]*alpha_inv_squared;
g4up[2][3] = g4up[3][2] = METRIC_PHYS[GUPYZ] - METRIC[SHIFTY]*METRIC[SHIFTZ]*alpha_inv_squared;
g4up[3][3] = METRIC_PHYS[GUPZZ] - METRIC[SHIFTZ]*METRIC[SHIFTZ]*alpha_inv_squared;
Appending to ../src/driver_conserv_to_prims.C
Instead of declaring new variables to store the currently known values of the conservative variables $\left\{\rho_{\star}, \tilde{S}_{i}, \tilde{\tau}\right\}$, we will simply use the flux variables $\left\{\rho_{\star}^{\rm flux}, \tilde{S}_{i}^{\rm flux}, \tilde{\tau}^{\rm flux}\right\}$ which are used by IllinoisGRMHD
as temporary storage. This is done for debugging purposes.
We also store the original values in the variables $\left\{\rho_{\star}^{\rm orig}, \tilde{S}_{i}^{\rm orig}, \tilde{\tau}^{\rm orig}\right\}$.
%%writefile -a $outfile_path__driver_conserv_to_prims__C
//FIXME: might slow down the code.
if(isnan(CONSERVS[RHOSTAR]*CONSERVS[STILDEX]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]*CONSERVS[TAUENERGY]*PRIMS[BX_CENTER]*PRIMS[BY_CENTER]*PRIMS[BZ_CENTER])) {
CCTK_VInfo(CCTK_THORNSTRING,"NAN FOUND: i,j,k = %d %d %d, x,y,z = %e %e %e , index=%d st_i = %e %e %e, rhostar = %e, tau = %e, Bi = %e %e %e, gij = %e %e %e %e %e %e, Psi6 = %e",
i,j,k,x[index],y[index],z[index],index,
CONSERVS[STILDEX],CONSERVS[STILDEY],CONSERVS[STILDEZ],CONSERVS[RHOSTAR],CONSERVS[TAUENERGY],
PRIMS[BX_CENTER],PRIMS[BY_CENTER],PRIMS[BZ_CENTER],METRIC_PHYS[GXX],METRIC_PHYS[GXY],METRIC_PHYS[GXZ],METRIC_PHYS[GYY],METRIC_PHYS[GYZ],METRIC_PHYS[GZZ],METRIC_LAP_PSI4[PSI6]);
}
// Here we use _flux variables as temp storage for original values of conservative variables.. This is used for debugging purposes only.
rho_star_flux[index] = CONSERVS[RHOSTAR];
st_x_flux[index] = CONSERVS[STILDEX];
st_y_flux[index] = CONSERVS[STILDEY];
st_z_flux[index] = CONSERVS[STILDEZ];
tau_flux[index] = CONSERVS[TAUENERGY];
CCTK_REAL rho_star_orig = CONSERVS[RHOSTAR];
CCTK_REAL mhd_st_x_orig = CONSERVS[STILDEX];
CCTK_REAL mhd_st_y_orig = CONSERVS[STILDEY];
CCTK_REAL mhd_st_z_orig = CONSERVS[STILDEZ];
CCTK_REAL tau_orig = CONSERVS[TAUENERGY];
Appending to ../src/driver_conserv_to_prims.C
In this part of the code we determine the primitives variables from the conservatives variables. Please note that this is only the driver function. The algorithm is discussed on Step 3 below.
We start by calling the apply_tau_floor()
function to ensure that the value of $\tilde\tau$ is not out of physical range. This function is discussed in the apply $\tilde\tau$ floor and enforce limits on primitives NRPy+ tutorial module.
Notice that this is only performed when $\rho_{\star}>0$. If this is not the case, we fix the issue by setting $\rho_{b} = \rho_{b}^{\rm atm}$ and the conservative to primitive algorithm is skipped altogether.
When $\rho_{\star}>0$, we call upon the primitive to conservatives algorithm through the harm_primitives_gammalaw_lowlevel()
function, described below.
%%writefile -a $outfile_path__driver_conserv_to_prims__C
int check=0;
struct output_stats stats;
stats.n_iter=0;
stats.vel_limited=0;
stats.failure_checker=0;
if(CONSERVS[RHOSTAR]>0.0) {
// Apply the tau floor
apply_tau_floor(index,tau_atm,rho_b_atm,Psi6threshold,PRIMS,METRIC,METRIC_PHYS,METRIC_LAP_PSI4,stats,eos, CONSERVS);
stats.font_fixed=0;
for(int ii=0;ii<3;ii++) {
check = harm_primitives_gammalaw_lowlevel(index,i,j,k,x,y,z,METRIC,METRIC_PHYS,METRIC_LAP_PSI4,
CONSERVS,PRIMS, g4dn,g4up, stats,eos);
if(check==0) ii=4;
else stats.failure_checker+=100000;
}
} else {
stats.failure_checker+=1;
// Set to atmosphere if rho_star<0.
//FIXME: FOR GAMMA=2 ONLY:
PRIMS[RHOB] = rho_b_atm;
/* Set P = P_cold */
int polytropic_index = find_polytropic_K_and_Gamma_index(eos, rho_b_atm);
CCTK_REAL K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];
CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];
PRIMS[PRESSURE] = K_ppoly_tab*pow(rho_b_atm,Gamma_ppoly_tab);
PRIMS[VX] =-METRIC[SHIFTX];
PRIMS[VY] =-METRIC[SHIFTY];
PRIMS[VZ] =-METRIC[SHIFTZ];
rho_star_fix_applied++;
}
Appending to ../src/driver_conserv_to_prims.C
We now enforce that the values of the primitives variables are physically meaningful. This is done by calling the IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs()
function, which is described in the enforce physical limits on primitives and recomputive conservatives NRPy+ tutorial module.
%%writefile -a $outfile_path__driver_conserv_to_prims__C
// Enforce limits on primitive variables and recompute conservatives.
static const int already_computed_physical_metric_and_inverse=1;
IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs(already_computed_physical_metric_and_inverse,PRIMS,stats,eos,METRIC,g4dn,g4up, TUPMUNU,TDNMUNU,CONSERVS);
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
rho_star[index] = CONSERVS[RHOSTAR];
mhd_st_x[index] = CONSERVS[STILDEX];
mhd_st_y[index] = CONSERVS[STILDEY];
mhd_st_z[index] = CONSERVS[STILDEZ];
tau[index] = CONSERVS[TAUENERGY];
// Set primitives, and/or provide a better guess.
rho_b[index] = PRIMS[RHOB];
P[index] = PRIMS[PRESSURE];
vx[index] = PRIMS[VX];
vy[index] = PRIMS[VY];
vz[index] = PRIMS[VZ];
Appending to ../src/driver_conserv_to_prims.C
%%writefile -a $outfile_path__driver_conserv_to_prims__C
if(update_Tmunu) {
ww=0;
eTtt[index] = TDNMUNU[ww]; ww++;
eTtx[index] = TDNMUNU[ww]; ww++;
eTty[index] = TDNMUNU[ww]; ww++;
eTtz[index] = TDNMUNU[ww]; ww++;
eTxx[index] = TDNMUNU[ww]; ww++;
eTxy[index] = TDNMUNU[ww]; ww++;
eTxz[index] = TDNMUNU[ww]; ww++;
eTyy[index] = TDNMUNU[ww]; ww++;
eTyz[index] = TDNMUNU[ww]; ww++;
eTzz[index] = TDNMUNU[ww];
}
//Finally, we set h, the enthalpy:
//CCTK_REAL eps = P[index]/rho_b[index]/(GAMMA-1.0);
//h[index] = 1.0 + P[index]/rho_b[index] + eps;
/***************************************************************************************************************************/
// DIAGNOSTICS:
//Pressure cap hit?
/* FIXME
CCTK_REAL P_cold = rho_b[index]*rho_b[index];
if(P[index]/P_cold > 0.99*1e3 && rho_b[index]>100.0*rho_b_atm) {
if(exp(phi[index]*6.0) <= Psi6threshold) pressure_cap_hit++;
}
*/
//Now we compute the difference between original & new conservatives, for diagnostic purposes:
error_int_numer += fabs(tau[index] - tau_orig) + fabs(rho_star[index] - rho_star_orig) +
fabs(mhd_st_x[index] - mhd_st_x_orig) + fabs(mhd_st_y[index] - mhd_st_y_orig) + fabs(mhd_st_z[index] - mhd_st_z_orig);
error_int_denom += tau_orig + rho_star_orig + fabs(mhd_st_x_orig) + fabs(mhd_st_y_orig) + fabs(mhd_st_z_orig);
if(stats.font_fixed==1) font_fixes++;
vel_limited_ptcount+=stats.vel_limited;
if(check!=0) {
failures++;
if(exp(METRIC[PHI]*6.0)>Psi6threshold) {
failures_inhoriz++;
pointcount_inhoriz++;
}
}
pointcount++;
/***************************************************************************************************************************/
failure_checker[index] = stats.failure_checker;
n_iter += stats.n_iter;
}
/*
gettimeofday(&end, NULL);
seconds = end.tv_sec - start.tv_sec;
useconds = end.tv_usec - start.tv_usec;
mtime = ((seconds) * 1000 + useconds/1000.0) + 0.999; // We add 0.999 since mtime is a long int; this rounds up the result before setting the value. Here, rounding down is incorrect.
solutions per second: cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2] / ((CCTK_REAL)mtime/1000.0),
*/
if(CCTK_Equals(verbose, "essential") || CCTK_Equals(verbose, "essential+iteration output")) {
CCTK_VInfo(CCTK_THORNSTRING,"C2P: Lev: %d NumPts= %d | Fixes: Font= %d VL= %d rho*= %d | Failures: %d InHoriz= %d / %d | Error: %.3e, ErrDenom: %.3e | %.2f iters/gridpt",
(int)GetRefinementLevel(cctkGH),
pointcount,font_fixes,vel_limited_ptcount,rho_star_fix_applied,
failures,
failures_inhoriz,pointcount_inhoriz,
error_int_numer/error_int_denom,error_int_denom,
(double)n_iter/( (double)(cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]) ));
}
if(pressure_cap_hit!=0) {
//CCTK_VInfo(CCTK_THORNSTRING,"PRESSURE CAP HIT %d TIMES! Outputting debug file!",pressure_cap_hit);
}
// Very useful con2prim debugger. If the primitives (con2prim) solver fails, this will output all data needed to
// debug where and why the solver failed. Strongly suggested for experimenting with new fixes.
if(conserv_to_prims_debug==1 && error_int_numer/error_int_denom > 0.05) {
ofstream myfile;
char filename[100];
srand(time(NULL));
sprintf(filename,"primitives_debug-%e.dat",error_int_numer/error_int_denom);
//Alternative, for debugging purposes as well:
//srand(time(NULL));
//sprintf(filename,"primitives_debug-%d.dat",rand());
myfile.open (filename, ios::out | ios::binary);
//myfile.open ("data.bin", ios::out | ios::binary);
myfile.write((char*)cctk_lsh, 3*sizeof(int));
myfile.write((char*)&GAMMA_SPEED_LIMIT, 1*sizeof(CCTK_REAL));
myfile.write((char*)&rho_b_max, 1*sizeof(CCTK_REAL));
myfile.write((char*)&rho_b_atm, 1*sizeof(CCTK_REAL));
myfile.write((char*)&tau_atm, 1*sizeof(CCTK_REAL));
myfile.write((char*)&Psi6threshold, 1*sizeof(CCTK_REAL));
myfile.write((char*)&update_Tmunu, 1*sizeof(int));
myfile.write((char*)&neos, 1*sizeof(int));
myfile.write((char*)&Gamma_th, 1*sizeof(CCTK_REAL));
myfile.write((char*)&K_ppoly_tab0, 1*sizeof(CCTK_REAL));
myfile.write((char*)Gamma_ppoly_tab_in, neos*sizeof(CCTK_REAL));
myfile.write((char*)rho_ppoly_tab_in, (neos-1)*sizeof(CCTK_REAL));
int fullsize=cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];
myfile.write((char*)x, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)y, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)z, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char *)failure_checker, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTtt, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTtx, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTty, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTtz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTxx, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTxy, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTxz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTyy, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTyz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)eTzz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)alp, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)gxx, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)gxy, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)gxz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)gyy, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)gyz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)gzz, fullsize*sizeof(CCTK_REAL));
myfile.write((char *)psi_bssn, fullsize*sizeof(CCTK_REAL));
myfile.write((char*)phi_bssn, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtxx, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtxy, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtxz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtyy, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtyz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtzz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtupxx, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtupxy, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtupxz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtupyy, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtupyz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)gtupzz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)betax, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)betay, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)betaz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)lapm1, (fullsize)*sizeof(CCTK_REAL));
// HERE WE USE _flux variables as temp storage for original values of conservative variables.. This is used for debugging purposes only.
myfile.write((char*)tau_flux, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)st_x_flux, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)st_y_flux, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)st_z_flux, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)rho_star_flux, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)Bx, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)By, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)Bz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)vx, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)vy, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)vz, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)P, (fullsize)*sizeof(CCTK_REAL));
myfile.write((char*)rho_b,(fullsize)*sizeof(CCTK_REAL));
int checker=1063; myfile.write((char*)&checker,sizeof(int));
myfile.close();
CCTK_VInfo(CCTK_THORNSTRING,"Finished writing %s",filename);
}
#ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER
return 0; // int main() requires an integer be returned
#endif
}
#include "harm_primitives_lowlevel.C"
Appending to ../src/driver_conserv_to_prims.C
harm_primitives_lowlevel.C
[Back to top]¶We will now begin documenting the harm_primitives_lowlevel.C
code. We will start by declaring the harm_primitives_gammalaw_lowlevel()
function and loading up the EOS parameters:
where we are currently assuming the single polytrope EOS
$$ P = \kappa \rho_{b}^{\Gamma}\ . $$HARM
[Back to top]¶In this section we will describe the final set of manipulations that are required before calling the conservative to primitive algorithm. These manipulations are necessary because we make use of the HARM
software, which uses a different set of conservative/primitives variables than IllinoisGRMHD
.
Notice that IllinoisGRMHD
uses the set of conservative variables $\boldsymbol{C}_{\rm IGM} = \left\{\rho_{\star},\tilde{\tau},\tilde{S}_{i},\tilde{B}^{i}\right\}$ and the primitive variables $\boldsymbol{P}_{\rm IGM} = \left\{\rho_{b},P,v^{i},B^{i}\right\}$, while HARM
makes use of the conservative variables $\boldsymbol{C}_{\rm HARM} = \left\{\sqrt{-g}\rho_{b}u^{0},\sqrt{-g}\left(T^{0}_{\ 0}+\rho_{b}u^{0}\right),\sqrt{-g}T^{0}_{\ i},\sqrt{-g}B^{i}_{\rm HARM}\right\}$ and primitive variables $\boldsymbol{P}_{\rm HARM} = \left\{\rho_{b},u,\tilde{u}^{i},B^{i}_{\rm HARM}\right\}$. The key step here is then writing HARM
's variables in terms of IllinoisGRMHD
's variables.
The first variable we will need to compute is ${\rm detg} := \sqrt{-g}$, where $g=\det(g_{\mu\nu})$, with $g_{\mu\nu}$ the physical ADM 4-metric evaluated above. To relate $\sqrt{-g}$ with IllinoisGRMHD
's variables, we will make use of the well known relation (see eq. 2.124 in Baumgarte & Shapiro's Numerical Relativity)
%%writefile $outfile_path__harm_primitives_lowlevel__C
inline int harm_primitives_gammalaw_lowlevel(const int index,const int i,const int j,const int k,CCTK_REAL *X,CCTK_REAL *Y,CCTK_REAL *Z,
CCTK_REAL *METRIC,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,
CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,
CCTK_REAL g4dn[NDIM][NDIM],CCTK_REAL g4up[NDIM][NDIM],
struct output_stats &stats, eos_struct &eos) {
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
Writing ../src/harm_primitives_lowlevel.C
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
// declare some variables for HARM.
CCTK_REAL U[NPR];
CCTK_REAL prim[NPR];
CCTK_REAL detg = METRIC_LAP_PSI4[LAPSE]*METRIC_LAP_PSI4[PSI6]; // == alpha sqrt{gamma} = alpha Psi^6
// Check to see if the metric is positive-definite.
// Note that this will slow down the code, and if the metric doesn't obey this, the run is probably too far gone to save,
// though if it happens deep in the horizon, it might resurrect the run.
/*
CCTK_REAL lam1,lam2,lam3;
CCTK_REAL M11 = METRIC[GXX], M12=METRIC[GXY], M13=METRIC[GXZ], M22=METRIC[GYY], M23=METRIC[GYZ], M33=METRIC[GZZ];
eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,M11, M12, M13, M22, M23, M33);
if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0.0) {
// Metric is not positive-defitive, reset the physical metric to be conformally-flat.
METRIC_PHYS[GXX] = METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GXY] = 0.0;
METRIC_PHYS[GXZ] = 0.0;
METRIC_PHYS[GYY] = METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GYZ] = 0.0;
METRIC_PHYS[GZZ] = METRIC_LAP_PSI4[PSI4];
METRIC_PHYS[GUPXX] = METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPXY] = 0.0;
METRIC_PHYS[GUPXZ] = 0.0;
METRIC_PHYS[GUPYY] = METRIC_LAP_PSI4[PSIM4];
METRIC_PHYS[GUPYZ] = 0.0;
METRIC_PHYS[GUPZZ] = METRIC_LAP_PSI4[PSIM4];
}
*/
Appending to ../src/harm_primitives_lowlevel.C
Now we must relate HARM
's $B^{i}_{\rm HARM}$ with the variables used by IllinoisGRMHD
. We can start by looking at eqs. (23), (24), and (31) in Duez et al. (2005) to find the following useful relations
Now, if we look at eqs. (16) and (17) in Gammie & McKinney (2003) we find the following relations
$$ {\rm HARM:}\ \left\{ \begin{align} b^{0} &= u_{i}B^{i}_{\rm HARM}\ ,\\ b^{i} &= \frac{B^{i}_{\rm HARM} + b^{0}u^{i}}{u^{0}}\ , \end{align} \right. $$from which we can then find the relation
$$ \boxed{B^{i}_{\rm HARM} = \frac{B^{i}}{\alpha\sqrt{4\pi}}}\ . $$%%writefile -a $outfile_path__harm_primitives_lowlevel__C
// Note that ONE_OVER_SQRT_4PI gets us to the object
// referred to as B^i in the Noble et al paper (and
// apparently also in the comments to their code).
// This is NOT the \mathcal{B}^i, which differs by
// a factor of the lapse.
CCTK_REAL BxL_over_alpha_sqrt_fourpi = PRIMS[BX_CENTER]*METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;
CCTK_REAL ByL_over_alpha_sqrt_fourpi = PRIMS[BY_CENTER]*METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;
CCTK_REAL BzL_over_alpha_sqrt_fourpi = PRIMS[BZ_CENTER]*METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;
Appending to ../src/harm_primitives_lowlevel.C
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
CCTK_REAL rho_b_oldL = PRIMS[RHOB];
CCTK_REAL P_oldL = PRIMS[PRESSURE];
CCTK_REAL vxL = PRIMS[VX];
CCTK_REAL vyL = PRIMS[VY];
CCTK_REAL vzL = PRIMS[VZ];
Appending to ../src/harm_primitives_lowlevel.C
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
/*
-- 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.
// / 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 / //
(above equations have been fixed by Yuk Tung & Zach)
*/
// U[NPR] = conserved variables (current values on input/output);
// g4dn[NDIM][NDIM] = covariant form of the 4-metric ;
// g4up[NDIM][NDIM] = contravariant form of the 4-metric ;
// gdet = sqrt( - determinant of the 4-metric) ;
// prim[NPR] = primitive variables (guess on input, calculated values on
// output if there are no problems);
// U[1] =
// U[2-4] = stildei + rhostar
CCTK_REAL rho_star_orig = CONSERVS[RHOSTAR];
CCTK_REAL mhd_st_x_orig = CONSERVS[STILDEX];
CCTK_REAL mhd_st_y_orig = CONSERVS[STILDEY];
CCTK_REAL mhd_st_z_orig = CONSERVS[STILDEZ];
CCTK_REAL tau_orig = CONSERVS[TAUENERGY];
Appending to ../src/harm_primitives_lowlevel.C
We will now start preparing the primitives for the conservative-to-primitive algorithm, which employs the Newton-Raphson method.
We offer the algorithm the initial guess
$$ \boxed{ {\rm Guess\ \#1:}\ \left\{ \begin{align} \rho_{b}^{\rm guess} &= \frac{\rho_{\star}}{\psi^{6}}\\ P_{\rm guess} &= \kappa \left(\rho_{b}^{\rm guess}\right)^{\Gamma}\\ u0 &= \frac{1}{\alpha}\\ v^{i} &= -\beta^{i} \end{align} \right. }\ . $$If this initial guess causes the Newton-Raphson method to fail, we try a second guess
$$ \boxed{ {\rm Guess\ \#2:}\ \left\{ \begin{align} \rho_{b}^{\rm guess} &= 100\rho_{b}^{\rm atm}\\ P_{\rm guess} &= \kappa \left(\rho_{b}^{\rm guess}\right)^{\Gamma}\\ u0 &= \frac{1}{\alpha}\\ v^{i} &= -\beta^{i} \end{align} \right. }\ . $$%%writefile -a $outfile_path__harm_primitives_lowlevel__C
// Other ideas for setting the gamma speed limit
//CCTK_REAL GAMMA_SPEED_LIMIT = 100.0;
//if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) GAMMA_SPEED_LIMIT=500.0;
//if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) GAMMA_SPEED_LIMIT=100.0;
//FIXME: Only works if poisoning is turned on. Otherwise will access unknown memory. This trick alone speeds up the whole code (Cowling) by 2%.
//int startguess=0;
//if(std::isnan(PRIMS[VX])) startguess=1;
int startguess=1;
CCTK_REAL u0L=1.0;
CCTK_REAL K_ppoly_tab,Gamma_ppoly_tab;
for(int which_guess=startguess;which_guess<3;which_guess++) {
int check;
if(which_guess==1) {
//Use a different initial guess:
rho_b_oldL = CONSERVS[RHOSTAR]/METRIC_LAP_PSI4[PSI6];
/**********************************
* Piecewise Polytropic EOS Patch *
* Finding Gamma_ppoly_tab and K_ppoly_tab *
**********************************/
/* Here we use our newly implemented
* find_polytropic_K_and_Gamma() function
* to determine the relevant polytropic
* Gamma and K parameters to be used
* within this function.
*/
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL);
K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];
Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];
// After that, we compute P_cold
P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab);
u0L = METRIC_LAP_PSI4[LAPSEINV];
vxL = -METRIC[SHIFTX];
vyL = -METRIC[SHIFTY];
vzL = -METRIC[SHIFTZ];
}
if(which_guess==2) {
//Use atmosphere as initial guess:
rho_b_oldL = 100.0*rho_b_atm;
/**********************************
* Piecewise Polytropic EOS Patch *
* Finding Gamma_ppoly_tab and K_ppoly_tab *
**********************************/
/* Here we use our newly implemented
* find_polytropic_K_and_Gamma() function
* to determine the relevant polytropic
* Gamma and K parameters to be used
* within this function.
*/
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL);
K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];
Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];
// After that, we compute P_cold
P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab);
u0L = METRIC_LAP_PSI4[LAPSEINV];
vxL = -METRIC[SHIFTX];
vyL = -METRIC[SHIFTY];
vzL = -METRIC[SHIFTZ];
}
Appending to ../src/harm_primitives_lowlevel.C
We will now relate $\boldsymbol{C}_{\rm HARM}$ to $\boldsymbol{C}_{\rm IGM}$. As previously mentioned, we have
$$ \begin{align} \boldsymbol{C}_{\rm IGM} &= \left\{\rho_{\star},\tilde{\tau},\tilde{S}_{i},\tilde{B}^{i}\right\}\ ,\\ \boldsymbol{C}_{\rm HARM} &= \left\{\sqrt{-g}\rho_{b}u^{0},\sqrt{-g}\left(T^{0}_{\ 0}+\rho_{b}u^{0}\right),\sqrt{-g}T^{0}_{\ i},\sqrt{-g}B^{i}_{\rm HARM}\right\} \equiv \left\{\rho_{\rm HARM},\tilde{\tau}_{\rm HARM},\tilde{S}_{i}^{\rm HARM},\tilde{B}^{i}_{\rm HARM}\right\}\ . \end{align} $$The following relations immediately hold
$$ \boxed{ \begin{align} \rho^{\rm HARM}_{\star} &= \rho_{\star}\\ \tilde{S}_{i}^{\rm HARM} &= \tilde{S}_{i} \end{align} }\ . $$As previously derived in Step 3.a.ii, we then have
$$ \boxed{\tilde{B}^{i}_{\rm HARM} = \underbrace{\left(\alpha\sqrt{\gamma}\right)}_{\rm detg}\underbrace{\left(\frac{B^{i}}{\alpha\sqrt{4\pi}}\right)}_{\rm BiL\_over\_alpha\_sqrt\_fourpi}}\ . $$Finally, we must relate $\tilde{\tau}_{\rm HARM} = \sqrt{-g}\left(T^{0}_{\ 0}+\rho_{b}u^{0}\right)$ to the IllinoisGRMHD
variables. First, consider the identity
Consider, also, the identity
$$ \tilde{\tau} = \alpha^{2}\sqrt{\gamma}T^{00} - \rho_{\star} \implies T^{00} = \frac{\left(\tilde\tau + \rho_{\star}\right)}{\alpha^{2}\sqrt{\gamma}}\ . $$Then
$$ \begin{align} T^{0}_{0} &= g_{0\mu}T^{0\mu}\\ &= g_{00}T^{00} + g_{0i}T^{0i}\\ &= \left(-\alpha^{2}+\beta_{\ell}\beta^{\ell}\right)T^{00} + \beta_{i}T^{0i}\\ &= -\alpha^{2}T^{00} + \beta_{i}\left[\beta^{i}T^{00} + T^{i0}\right]\\ &= -\alpha^{2}\left[\frac{\left(\tilde\tau + \rho_{\star}\right)}{\alpha^{2}\sqrt{\gamma}}\right] + \frac{\beta^{i}S_{i}}{\alpha}\\ \implies T^{0}_{0} &= -\frac{\tilde\tau+\rho_{\star}}{\sqrt{\gamma}} + \frac{\beta^{i}S_{i}}{\alpha}\ . \end{align} $$Finally, we can compute
$$ \begin{align} \tilde\tau_{\rm HARM} &= \sqrt{-g}\left(T^{0}_{0} + \rho_{b}u^{0}\right)\\ &= \alpha\sqrt{\gamma}T^{0}_{0} + \alpha\sqrt{\gamma}\rho_{b}u^{0}\\ &= \alpha\sqrt{\gamma}\left[-\frac{\tilde\tau+\rho_{\star}}{\sqrt{\gamma}} + \frac{\beta^{i}S_{i}}{\alpha}\right] + \rho_{\star}\\ &= -\alpha\tilde{\tau} -\alpha\rho_{\star} + \beta^{i}\left(\sqrt{\gamma}S_{i}\right) + \rho_{\star}\\ \implies &\boxed{\tilde\tau_{\rm HARM} = -\alpha\tilde{\tau} - \left(\alpha-1\right)\rho_{\star} + \beta^{i}\tilde{S}_{i}}\ . \end{align} $$%%writefile -a $outfile_path__harm_primitives_lowlevel__C
// Fill the array of conserved variables according to the wishes of Utoprim_2d.
U[RHO] = CONSERVS[RHOSTAR];
U[UU] = -CONSERVS[TAUENERGY]*METRIC_LAP_PSI4[LAPSE] - (METRIC_LAP_PSI4[LAPSE]-1.0)*CONSERVS[RHOSTAR] +
METRIC[SHIFTX]*CONSERVS[STILDEX] + METRIC[SHIFTY]*CONSERVS[STILDEY] + METRIC[SHIFTZ]*CONSERVS[STILDEZ] ; // note the minus sign on tau
U[UTCON1] = CONSERVS[STILDEX];
U[UTCON2] = CONSERVS[STILDEY];
U[UTCON3] = CONSERVS[STILDEZ];
U[BCON1] = detg*BxL_over_alpha_sqrt_fourpi;
U[BCON2] = detg*ByL_over_alpha_sqrt_fourpi;
U[BCON3] = detg*BzL_over_alpha_sqrt_fourpi;
Appending to ../src/harm_primitives_lowlevel.C
Keep in mind that for the primitive variables, we are not intereted in providing exact relations. Instead, what we are interested in doing is providing an educated guess of what it should look like in hopes that the Newton-Raphson method then converges to the correct values. With that in mind, the primitive variables are then initialized to:
$$ \boxed{ \begin{align} \rho_{\rm HARM} &= \rho_{b}\\ u &= \frac{P}{\Gamma_{\rm poly} - 1}\\ \tilde{u}^{i} &= u^{0}\left(v^{i}+\beta^{i}\right)\\ B^{i}_{\rm HARM} &= \frac{B^{i}}{\alpha\sqrt{4\pi}} \end{align}\ , } $$where $\Gamma_{\rm poly}$ stands for the local polytropic $\Gamma$-factor.
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
CCTK_REAL uL = P_oldL/(Gamma_ppoly_tab - 1.0);
CCTK_REAL utxL = u0L*(vxL + METRIC[SHIFTX]);
CCTK_REAL utyL = u0L*(vyL + METRIC[SHIFTY]);
CCTK_REAL utzL = u0L*(vzL + METRIC[SHIFTZ]);
prim[RHO] = rho_b_oldL;
prim[UU] = uL;
prim[UTCON1] = utxL;
prim[UTCON2] = utyL;
prim[UTCON3] = utzL;
prim[BCON1] = BxL_over_alpha_sqrt_fourpi;
prim[BCON2] = ByL_over_alpha_sqrt_fourpi;
prim[BCON3] = BzL_over_alpha_sqrt_fourpi;
Appending to ../src/harm_primitives_lowlevel.C
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
/*************************************************************/
// CALL HARM PRIMITIVES SOLVER:
check = Utoprim_2d(eos, U, g4dn, g4up, detg, prim,stats.n_iter);
// Note that we have modified this solver, so that nearly 100%
// of the time it yields either a good root, or a root with
// negative epsilon (i.e., pressure).
/*************************************************************/
Appending to ../src/harm_primitives_lowlevel.C
If the conservative-to-primitive solver fails to converge, we apply the procedure suggested by Font et al. (1998). The algorithm can be summarized as the requirement that
$$ P=P_{\rm cold} = \kappa \rho^{\Gamma_{\rm cold}}_{b}\ , $$and then recomputing the velocities $u_{i}$. We will describe this procedure in detail in Step 4 below.
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
// Use the new Font fix subroutine
int font_fix_applied=0;
if(check!=0) {
font_fix_applied=1;
CCTK_REAL u_xl=1e100, u_yl=1e100, u_zl=1e100; // Set to insane values to ensure they are overwritten.
/************************
* New Font fix routine *
************************/
check = font_fix__hybrid_EOS(u_xl,u_yl,u_zl, CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4, eos);
Appending to ../src/harm_primitives_lowlevel.C
Now we evaluate
$$ \tilde{u}^{i} = \gamma^{ij}u_{i} \implies \boxed{ \left\{ \begin{align} \tilde{u}^{x} &= \gamma^{xx}u_{x} + \gamma^{xy}u_{y} + \gamma^{xz}u_{z}\\ \tilde{u}^{y} &= \gamma^{yx}u_{x} + \gamma^{yy}u_{y} + \gamma^{yz}u_{z}\\ \tilde{u}^{z} &= \gamma^{zx}u_{x} + \gamma^{zy}u_{y} + \gamma^{zz}u_{z} \end{align} \right. } $$%%writefile -a $outfile_path__harm_primitives_lowlevel__C
//Translate to HARM primitive now:
prim[UTCON1] = METRIC_PHYS[GUPXX]*u_xl + METRIC_PHYS[GUPXY]*u_yl + METRIC_PHYS[GUPXZ]*u_zl;
prim[UTCON2] = METRIC_PHYS[GUPXY]*u_xl + METRIC_PHYS[GUPYY]*u_yl + METRIC_PHYS[GUPYZ]*u_zl;
prim[UTCON3] = METRIC_PHYS[GUPXZ]*u_xl + METRIC_PHYS[GUPYZ]*u_yl + METRIC_PHYS[GUPZZ]*u_zl;
if (check==1) {
CCTK_VInfo(CCTK_THORNSTRING,"Font fix failed!");
CCTK_VInfo(CCTK_THORNSTRING,"i,j,k = %d %d %d, stats.failure_checker = %d x,y,z = %e %e %e , index=%d st_i = %e %e %e, rhostar = %e, Bi = %e %e %e, gij = %e %e %e %e %e %e, Psi6 = %e",i,j,k,stats.failure_checker,X[index],Y[index],Z[index],index,mhd_st_x_orig,mhd_st_y_orig,mhd_st_z_orig,rho_star_orig,PRIMS[BX_CENTER],PRIMS[BY_CENTER],PRIMS[BZ_CENTER],METRIC_PHYS[GXX],METRIC_PHYS[GXY],METRIC_PHYS[GXZ],METRIC_PHYS[GYY],METRIC_PHYS[GYZ],METRIC_PHYS[GZZ],METRIC_LAP_PSI4[PSI6]);
}
}
stats.failure_checker+=font_fix_applied*10000;
stats.font_fixed=font_fix_applied;
/*************************************************************/
Appending to ../src/harm_primitives_lowlevel.C
The procedure we follow here is similar to the one discussed in the inlined functions tutorial module.
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
if(check==0) {
//Now that we have found some solution, we first limit velocity:
//FIXME: Probably want to use exactly the same velocity limiter function here as in mhdflux.C
CCTK_REAL utx_new = prim[UTCON1];
CCTK_REAL uty_new = prim[UTCON2];
CCTK_REAL utz_new = prim[UTCON3];
//Velocity limiter:
CCTK_REAL gijuiuj = METRIC_PHYS[GXX]*SQR(utx_new ) +
2.0*METRIC_PHYS[GXY]*utx_new*uty_new + 2.0*METRIC_PHYS[GXZ]*utx_new*utz_new +
METRIC_PHYS[GYY]*SQR(uty_new) + 2.0*METRIC_PHYS[GYZ]*uty_new*utz_new +
METRIC_PHYS[GZZ]*SQR(utz_new);
CCTK_REAL au0m1 = gijuiuj/( 1.0+sqrt(1.0+gijuiuj) );
u0L = (au0m1+1.0)*METRIC_LAP_PSI4[LAPSEINV];
// *** Limit velocity
stats.vel_limited=0;
if (au0m1 > 0.9999999*(GAMMA_SPEED_LIMIT-1.0)) {
CCTK_REAL fac = sqrt((SQR(GAMMA_SPEED_LIMIT)-1.0)/(SQR(1.0+au0m1) - 1.0));
utx_new *= fac;
uty_new *= fac;
utz_new *= fac;
gijuiuj = gijuiuj * SQR(fac);
au0m1 = gijuiuj/( 1.0+sqrt(1.0+gijuiuj) );
// Reset rho_b and u0
u0L = (au0m1+1.0)*METRIC_LAP_PSI4[LAPSEINV];
prim[RHO] = rho_star_orig/(METRIC_LAP_PSI4[LAPSE]*u0L*METRIC_LAP_PSI4[PSI6]);
stats.vel_limited=1;
stats.failure_checker+=1000;
} //Finished limiting velocity
Appending to ../src/harm_primitives_lowlevel.C
%%writefile -a $outfile_path__harm_primitives_lowlevel__C
//The Font fix only sets the velocities. Here we set the pressure & density HARM primitives.
if(font_fix_applied==1) {
prim[RHO] = rho_star_orig/(METRIC_LAP_PSI4[LAPSE]*u0L*METRIC_LAP_PSI4[PSI6]);
//Next set P = P_cold:
CCTK_REAL P_cold;
/**********************************
* Piecewise Polytropic EOS Patch *
* Finding Gamma_ppoly_tab and K_ppoly_tab *
**********************************/
/* Here we use our newly implemented
* find_polytropic_K_and_Gamma() function
* to determine the relevant polytropic
* Gamma and K parameters to be used
* within this function.
*/
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,prim[RHO]);
K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];
Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];
// After that, we compute P_cold
P_cold = K_ppoly_tab*pow(prim[RHO],Gamma_ppoly_tab);
prim[UU] = P_cold/(Gamma_ppoly_tab-1.0);
} //Finished setting remaining primitives if there was a Font fix.
/* Set rho_b */
PRIMS[RHOB] = prim[RHO];
/***************
* PPEOS Patch *
* Hybrid EOS *
***************
*/
/* We now compute the pressure as a function
* of rhob, P_cold, eps_cold, and u = rhob*eps,
* using the function pressure_rho0_u(), which
* implements the equation:
* .-------------------------------------------------------------.
* | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |
* .-------------------------------------------------------------.
*/
PRIMS[PRESSURE] = pressure_rho0_u(eos, prim[RHO],prim[UU]);
/* Already set u0L. */
PRIMS[VX] = utx_new/u0L - METRIC[SHIFTX];
PRIMS[VY] = uty_new/u0L - METRIC[SHIFTY];
PRIMS[VZ] = utz_new/u0L - METRIC[SHIFTZ];
return 0;
} else {
//If we didn't find a root, then try again with a different guess.
}
}
CCTK_VInfo(CCTK_THORNSTRING,"Couldn't find root from: %e %e %e %e %e, rhob approx=%e, rho_b_atm=%e, Bx=%e, By=%e, Bz=%e, gij_phys=%e %e %e %e %e %e, alpha=%e",
tau_orig,rho_star_orig,mhd_st_x_orig,mhd_st_y_orig,mhd_st_z_orig,rho_star_orig/METRIC_LAP_PSI4[PSI6],rho_b_atm,PRIMS[BX_CENTER],PRIMS[BY_CENTER],PRIMS[BZ_CENTER],METRIC_PHYS[GXX],METRIC_PHYS[GXY],METRIC_PHYS[GXZ],METRIC_PHYS[GYY],METRIC_PHYS[GYZ],METRIC_PHYS[GZZ],METRIC_LAP_PSI4[LAPSE]);
return 1;
}
//#include "harm_u2p_util.c"
#include "harm_utoprim_2d.c"
#include "eigen.C"
#include "font_fix_hybrid_EOS.C"
Appending to ../src/harm_primitives_lowlevel.C
font_fix_hybrid_EOS.C
[Back to top]¶The Font et al. algorithm (henceforth Font Fix algorithm) can be summarized as follows. Font fixes occur at the atmospheric region, so we start by assuming that $P$ is given only by its cold part, i.e.
$$ P = P_{\rm cold} = K_{\rm atm} \rho_{b}^{\Gamma_{\rm atm}}\ , $$where $K_{\rm atm}$ and $\Gamma_{\rm atm}$ are the constants used in the atmosphere. Then, the specific internal energy is given by
$$ \begin{align} \epsilon &= \epsilon_{\rm cold}\\ &= \int d\rho \frac{P_{\rm cold}}{\rho^{2}}\\ &= K_{\rm atm}\int d\rho \rho^{\Gamma_{\rm atm}-2}\\ &= \frac{K_{\rm atm}\rho^{\Gamma_{\rm atm}-1}}{\Gamma_{\rm atm}-1}\ . \end{align} $$Having computed $P$ and $\epsilon$, we can compute the enthalpy, $h$, giving
$$ \begin{align} h &= 1 + \epsilon + \frac{P}{\rho}\\ &= 1 + \frac{K_{\rm atm}\rho^{\Gamma_{\rm atm}-1}}{\Gamma_{\rm atm}-1} + K_{\rm atm}\rho^{\Gamma_{\rm atm} - 1}\\ &= 1 + \left(\frac{1}{\Gamma_{\rm atm}-1}+1\right)K_{\rm atm}\rho^{\Gamma_{\rm atm} - 1}\\ \implies &\boxed{ h = 1 + \left(\frac{\Gamma_{\rm atm}}{\Gamma_{\rm atm}-1}\right)K_{\rm atm} \rho^{\Gamma_{\rm atm}-1} }\ . \end{align} $$We then run an iterative process that updates $\rho$ in a consistent way, based on the value of $h$. We now describe this process and its implementation.
We start by computing all basic quantities needed by the Font Fix algorithm:
$$ \boxed{ \begin{align} \bar{B}^{i} &= \frac{B^i}{\sqrt{4\pi}}\\ \bar{B}_{i} &= \gamma_{ij}\bar{B}^{j}\\ \bar{B}^{2} &= \bar{B}_{i}\bar{B}^{i}\\ \bar{B} &= \sqrt{\bar{B}^{2}}\\ \bar{B}\cdot\tilde{S} &= \bar{B}^{i}\tilde{S}_{i}\\ (\bar{B}&\cdot\tilde{S})^{2}\\ \hat{\bar{B}}\cdot\tilde{S} &= \hat{\bar{B}}^{i}\tilde{S}_{i} \equiv \left(\frac{\bar{B}^{i}}{\bar{B}}\right)\tilde{S}_{i}\\ \tilde{S}\cdot\tilde{S} &= \gamma^{ij}\tilde{S}_{i}\tilde{S}_{j} \end{align} }\ . $$%%writefile $outfile_path__font_fix_hybrid_EOS__C
/**********************************
* Piecewise Polytropic EOS Patch *
* Font fix: function call *
**********************************/
inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos) {
CCTK_REAL Bxbar = PRIMS[BX_CENTER]*ONE_OVER_SQRT_4PI;
CCTK_REAL Bybar = PRIMS[BY_CENTER]*ONE_OVER_SQRT_4PI;
CCTK_REAL Bzbar = PRIMS[BZ_CENTER]*ONE_OVER_SQRT_4PI;
CCTK_REAL Bbar_x = METRIC_PHYS[GXX]*Bxbar + METRIC_PHYS[GXY]*Bybar + METRIC_PHYS[GXZ]*Bzbar;
CCTK_REAL Bbar_y = METRIC_PHYS[GXY]*Bxbar + METRIC_PHYS[GYY]*Bybar + METRIC_PHYS[GYZ]*Bzbar;
CCTK_REAL Bbar_z = METRIC_PHYS[GXZ]*Bxbar + METRIC_PHYS[GYZ]*Bybar + METRIC_PHYS[GZZ]*Bzbar;
CCTK_REAL B2bar = Bxbar*Bbar_x + Bybar*Bbar_y + Bzbar*Bbar_z;
CCTK_REAL Bbar = sqrt(B2bar);
CCTK_REAL check_B_small = fabs(Bxbar)+fabs(Bybar)+fabs(Bzbar);
if (check_B_small>0 && check_B_small<1.e-150) {
// need to compute B2bar specially to prevent floating-point underflow
CCTK_REAL Bmax = fabs(Bxbar);
if (Bmax < fabs(Bybar)) Bmax=fabs(Bybar);
if (Bmax < fabs(Bzbar)) Bmax=fabs(Bzbar);
CCTK_REAL Bxtmp=Bxbar/Bmax, Bytemp=Bybar/Bmax, Bztemp=Bzbar/Bmax;
CCTK_REAL B_xtemp=Bbar_x/Bmax, B_ytemp=Bbar_y/Bmax, B_ztemp=Bbar_z/Bmax;
Bbar = sqrt(Bxtmp*B_xtemp + Bytemp*B_ytemp + Bztemp*B_ztemp)*Bmax;
}
CCTK_REAL BbardotS = Bxbar*CONSERVS[STILDEX] + Bybar*CONSERVS[STILDEY] + Bzbar*CONSERVS[STILDEZ];
CCTK_REAL BbardotS2 = BbardotS*BbardotS;
CCTK_REAL hatBbardotS = BbardotS/Bbar;
if (Bbar<1.e-300) hatBbardotS = 0.0;
CCTK_REAL Psim6 = 1.0/METRIC_LAP_PSI4[PSI6];
// Limit hatBbardotS
//CCTK_REAL max_gammav = 100.0;
//CCTK_REAL rhob_max = CONSERVS[RHOSTAR]*Psim6;
//CCTK_REAL hmax = 1.0 + gam_gamm1_kpoly*pow(rhob_max,gam1);
//CCTK_REAL abs_hatBbardotS_max = sqrt(SQR(max_gammav)-1.0)*CONSERVS[RHOSTAR]*hmax;
//if (fabs(hatBbardotS) > abs_hatBbardotS_max) {
// CCTK_REAL fac_reduce = abs_hatBbardotS_max/fabs(hatBbardotS);
// CCTK_REAL hatBbardotS_max = hatBbardotS*fac_reduce;
// CCTK_REAL Bbar_inv = 1.0/Bbar;
// CCTK_REAL hat_Bbar_x = Bbar_x*Bbar_inv;
// CCTK_REAL hat_Bbar_y = Bbar_y*Bbar_inv;
// CCTK_REAL hat_Bbar_z = Bbar_z*Bbar_inv;
// CCTK_REAL sub_fact = hatBbardotS_max - hatBbardotS;
// CONSERVS[STILDEX] += sub_fact*hat_Bbar_x;
// CONSERVS[STILDEY] += sub_fact*hat_Bbar_y;
// CONSERVS[STILDEZ] += sub_fact*hat_Bbar_z;
// hatBbardotS = hatBbardotS_max;
// BbardotS *= fac_reduce;
// BbardotS2 = BbardotS*BbardotS;
//}
CCTK_REAL sdots = METRIC_PHYS[GUPXX]*SQR(CONSERVS[STILDEX]) + METRIC_PHYS[GUPYY]*SQR(CONSERVS[STILDEY]) + METRIC_PHYS[GUPZZ]*SQR(CONSERVS[STILDEZ])
+ 2.0*( METRIC_PHYS[GUPXY]*CONSERVS[STILDEX]*CONSERVS[STILDEY] + METRIC_PHYS[GUPXZ]*CONSERVS[STILDEX]*CONSERVS[STILDEZ]
+ METRIC_PHYS[GUPYZ]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]);
Writing ../src/font_fix_hybrid_EOS.C
We start by looking at the dot product $\tilde{S}^{2}$. Recall that
$$ \tilde{S}_{i} = \left(\rho_{\star}h + \alpha\sqrt{\gamma}\, u^{0}b^{2}\right)u_{i}-\alpha\sqrt{\gamma}\, b^{0}b_{i}\ . $$If $\tilde{S}^{2} = 0$, then we must be in a region where $u_{i} = 0 = b_{i}$. In this case, we return
$$ \begin{align} \rho_{b} &= \psi^{-6}\rho_{\star}\ ,\\ u^{i} &= 0\ , \end{align} $$and terminate the function call.
%%writefile -a $outfile_path__font_fix_hybrid_EOS__C
CCTK_REAL rhob;
if (sdots<1.e-300) {
rhob = CONSERVS[RHOSTAR]*Psim6;
u_x=0.0; u_y=0.0; u_z=0.0;
return 0;
}
/* This test has some problem.
if (fabs(BbardotS2 - sdots*B2bar) > 1e-8) {
CCTK_VInfo(CCTK_THORNSTRING,"(Bbar dot S)^2, Bbar^2 * sdotS, %e %e",SQR(BbardotS),sdots*B2bar);
CCTK_VInfo(CCTK_THORNSTRING,"Cauchy-Schwartz inequality is violated!");
}
*/
Appending to ../src/font_fix_hybrid_EOS.C
If $\tilde{S}^{2} \neq 0$, then we move on to the iterative procedure previously mentioned. We start by setting the initial data based on eqs. (A52), (A53), and (A59) found in Appendix A of Zachariah et al. (2012):
$$ \boxed{ \begin{align} W_{0} &= \psi^{-6}\sqrt{\left(\hat{\bar{B}}\cdot\tilde{S}\right)^{2} + \rho_{\star}^{2}}\\ S_{{\rm fluid},0}^{2} &=\frac{W_{0}^{2}\left(\tilde{S}\cdot\tilde{S}\right)+\left(\bar{B}\cdot\tilde{S}\right)^{2}\left(\bar{B}^{2} + 2W_{0}\right)}{\left(W_{0} + \bar{B}^{2}\right)^{2}}\\ \rho_{0} &= \frac{\psi^{-6}\rho_{\star}}{\sqrt{1+\frac{S_{{\rm fluid},0}^{2}}{\rho_{\star}^{2}}}} \end{align} }\ . $$%%writefile -a $outfile_path__font_fix_hybrid_EOS__C
// Initial guess for W, S_fluid and rhob
CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6;
CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar);
CCTK_REAL rhob0 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]));
Appending to ../src/font_fix_hybrid_EOS.C
We now perform the following iterative process, which is again described in Appendix A of Zachariah et al. (2012). We refer the reader to eqs. (A60), (A61), and (A62).
This is done using the function font_fix__rhob_loop()
, which is documented in the inlined functions tutorial notebook.
%%writefile -a $outfile_path__font_fix_hybrid_EOS__C
//****************************************************************
// FONT FIX
// Impose Font fix when HARM primitives solver fails to find
// acceptable set of primitives.
//****************************************************************
/* Set the maximum number of iterations */
int maxits = 500;
/* Set the allowed tolerance */
CCTK_REAL tol = 1.e-15;
/* Declare basic variables */
int font_fix_status;
/**********************
* FONT FIX MAIN LOOP *
**********************
* Perform the font fix routine until convergence
* is obtained and the algorithm returns with no
* error. Every time the Font fix fails, increase
* the tolerance by a factor of 10.
*/
int font_fix_attempts = 5;
CCTK_REAL font_fix_tol_factor = 10.0;
for(int n=0; n<font_fix_attempts; n++) {
tol *= pow(font_fix_tol_factor,n);
font_fix_status = font_fix__rhob_loop(maxits,tol, W0,Sf20,Psim6,sdots,BbardotS2,B2bar, CONSERVS,eos, rhob0,rhob);
rhob0 = rhob;
if(font_fix_status==0) break;
}
Appending to ../src/font_fix_hybrid_EOS.C
Finally, we return $\rho_{b}$ and $u_{i}$. In the relations below, $N$ indicates the last computed value of $\rho$ obtained by our iterative process. The quantities evaluated here are
$$ \boxed{ \begin{align} \rho_{b} &= \rho_{N}\\ \gamma_{v} &= \frac{\psi^{-6}\rho_{\star}}{\rho_{b}}\\ f_{1} &= \frac{\psi^{6}\left(\bar{B}\cdot\tilde{S}\right)}{\gamma_{v}\rho_{\star} h}\\ f_{2} &= \left(\rho_{\star}h + \psi^{6}\frac{\bar{B}^{2}}{\gamma_{v}}\right)^{-1}\\ u_{i} &= f_{2}\left(\tilde{S}_{i} + f_{1}\bar{B}_{i}\right) \end{align} }\ . $$%%writefile -a $outfile_path__font_fix_hybrid_EOS__C
//**************************************************************************************************************
/* Font fix works! */
/* First compute P_cold, eps_cold, then h = h_cold */
CCTK_REAL P_cold, eps_cold;
compute_P_cold__eps_cold(eos,rhob, P_cold,eps_cold);
CCTK_REAL h = 1.0 + eps_cold + P_cold/rhob;
/* Then compute gamma_v using equation (A19) in
* Etienne et al. (2011) [https://arxiv.org/pdf/1112.0568.pdf]
* .-----------------------------------------.
* | gamma_v = psi^{-6} * (rho_star / rho_b) |
* .-----------------------------------------.
*/
CCTK_REAL gammav = CONSERVS[RHOSTAR]*Psim6/rhob;
/* Finally, compute u_{i} */
CCTK_REAL rhosh = CONSERVS[RHOSTAR]*h;
CCTK_REAL fac1 = METRIC_LAP_PSI4[PSI6]*BbardotS/(gammav*rhosh);
CCTK_REAL fac2 = 1.0/(rhosh + METRIC_LAP_PSI4[PSI6]*B2bar/gammav);
u_x = fac2*(CONSERVS[STILDEX] + fac1*Bbar_x);
u_y = fac2*(CONSERVS[STILDEY] + fac1*Bbar_y);
u_z = fac2*(CONSERVS[STILDEZ] + fac1*Bbar_z);
return 0;
}
Appending to ../src/font_fix_hybrid_EOS.C
%%writefile $outfile_path__harm_primitives_headers__h
/***********************************************************************************
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
***********************************************************************************/
#ifndef HARM_PRIMITIVES_HEADERS_H_
#define HARM_PRIMITIVES_HEADERS_H_
static const int NPR =8;
static const int NDIM=4;
/* Adiabatic index used for the state equation */
//#define GAMMA (2.0)
static const CCTK_REAL G_ISOTHERMAL = 1.0;
/* use K(s)=K(r)=const. (G_ATM = GAMMA) of time or T = T(r) = const. of time (G_ATM = 1.) */
/*
#define USE_ISENTROPIC 1
#if( USE_ISENTROPIC )
#define G_ATM GAMMA
#else
#define G_ATM G_ISOTHERMAL
#endif
*/
static const int MAX_NEWT_ITER=30; /* Max. # of Newton-Raphson iterations for find_root_2D(); */
//#define MAX_NEWT_ITER 300 /* Max. # of Newton-Raphson iterations for find_root_2D(); */
static const CCTK_REAL NEWT_TOL =1.0e-10; /* Min. of tolerance allowed for Newton-Raphson iterations */
static const CCTK_REAL MIN_NEWT_TOL=1.0e-10; /* Max. of tolerance allowed for Newton-Raphson iterations */
static const int EXTRA_NEWT_ITER=0; /* ZACH SAYS: Original value = 2. But I don't think this parameter > 0 is warranted. Just slows the code for no reason, since our tolerances are fine. */
static const CCTK_REAL NEWT_TOL2 =1.0e-15; /* TOL of new 1D^*_{v^2} gnr2 method */
static const CCTK_REAL MIN_NEWT_TOL2=1.0e-10; /* TOL of new 1D^*_{v^2} gnr2 method */
static const CCTK_REAL W_TOO_BIG =1.e20; /* \gamma^2 (\rho_0 + u + p) is assumed
to always be smaller than this. This
is used to detect solver failures */
static const CCTK_REAL UTSQ_TOO_BIG =1.e20; /* \tilde{u}^2 is assumed to be smaller
than this. Used to detect solver
failures */
static const CCTK_REAL FAIL_VAL =1.e30; /* Generic value to which we set variables when a problem arises */
static const CCTK_REAL NUMEPSILON=(2.2204460492503131e-16);
/* some mnemonics */
/* for primitive variables */
static const int RHO =0;
static const int UU =1;
static const int UTCON1 =2;
static const int UTCON2 =3;
static const int UTCON3 =4;
static const int BCON1 =5;
static const int BCON2 =6;
static const int BCON3 =7;
/* for conserved variables */
static const int QCOV0 =1;
static const int QCOV1 =2;
static const int QCOV2 =3;
static const int QCOV3 =4;
/********************************************************************************************/
// Function prototype declarations:
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);
inline int harm_primitives_gammalaw_lowlevel(const int index,const int i,const int j,const int k,CCTK_REAL *X,CCTK_REAL *Y,CCTK_REAL *Z,
CCTK_REAL *METRIC,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,
CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,
CCTK_REAL g4dn[NDIM][NDIM],CCTK_REAL g4up[NDIM][NDIM],
output_stats &stats,eos_struct &eos);
inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos);
void eigenvalues_3by3_real_sym_matrix(CCTK_REAL & lam1, CCTK_REAL & lam2, CCTK_REAL & lam3,
CCTK_REAL M11, CCTK_REAL M12, CCTK_REAL M13, CCTK_REAL M22, CCTK_REAL M23, CCTK_REAL M33);
/********************************************************************************************/
#endif /* HARM_PRIMITIVES_HEADERS_H_ */
Writing ../src/harm_primitives_headers.h
# 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/driver_conserv_to_prims.C"
original_IGM_file_name = "driver_conserv_to_prims-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__driver_conserv_to_prims__C = !diff $original_IGM_file_path $outfile_path__driver_conserv_to_prims__C
if Validation__driver_conserv_to_prims__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 driver_conserv_to_prims.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for driver_conserv_to_prims.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__driver_conserv_to_prims__C:
print(diff_line)
Validation test for driver_conserv_to_prims.C: FAILED! Diff: 1c1 < /* We evolve forward in time a set of functions called the --- > /* We evolve forward in time a set of functions called the 3,4c3,4 < * are updated, we must solve for the primitive variables < * (rho, pressure, velocities) using a Newton-Raphson --- > * are updated, we must solve for the primitive variables > * (rho, pressure, velocities) using a Newton-Raphson 6c6 < * of the MHD equations again. --- > * of the MHD equations again. 9,12c9,12 < * Raphson solver. Truncation errors in conservative < * variables can lead to no physical solutions in < * primitive variables. We correct for these errors here < * through a number of tricks described in the appendices --- > * Raphson solver. Truncation errors in conservative > * variables can lead to no physical solutions in > * primitive variables. We correct for these errors here > * through a number of tricks described in the appendices 15,16c15,16 < * This is a wrapper for the 2d solver of Noble et al. See < * harm_utoprim_2d.c for references and copyright notice --- > * This is a wrapper for the 2d solver of Noble et al. See > * harm_utoprim_2d.c for references and copyright notice 19,20c19,20 < * < * For optimal compatibility, this wrapper is licensed under --- > * > * For optimal compatibility, this wrapper is licensed under 23c23 < * Note that this code assumes a simple gamma law for the --- > * Note that this code assumes a simple gamma law for the 27,28c27 < < #include "cctk.h" --- > // Standard #include's 32d30 < #include <sys/time.h> 35a34,39 > > > #ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER > #include "standalone_conserv_to_prims_main_function.h" > #else > #include "cctk.h" 41a46 > #include "harm_u2p_util.c" 50a56 > #endif 52,54c58,68 < // FIXME: This is definitely one of the places we need to look at later. < // Will leave this comment to emphasize the "FIXME" below. < // FIXME: only for single gamma-law EOS. --- > /********************************** > * Piecewise Polytropic EOS Patch * > * Setting up the EOS struct * > **********************************/ > /* > * The short piece of code below takes care > * of initializing the EOS parameters. > * Please refer to the "inlined_functions.C" > * source file for the documentation on the > * function. > */ 56,63c70,71 < eos.neos=neos; < eos.K_poly=K_poly; < eos.rho_tab[0]=rho_tab[0]; < eos.P_tab[0]=P_tab[0]; < eos.gamma_th=gamma_th; < eos.eps_tab[0]=eps_tab[0]; < eos.k_tab[0]=k_tab[0]; eos.k_tab[1]=k_tab[1]; < eos.gamma_tab[0]=gamma_tab[0]; eos.gamma_tab[1]=gamma_tab[1]; --- > initialize_EOS_struct_from_input(eos); > 71a80,81 > > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 82a93,94 > #endif > 93,95d104 < int gamma_equals2 = 1; < if (fabs(gamma_th-2.0) > 1.e-10) gamma_equals2 = 0; < 109a119 > 139a150 > 150a162 > 152a165 > 155c168,169 < --- > > 169a184 > 176a192 > 190c206 < --- > 202a219 > 210c227 < --- > 223a241 > 245,250c263,270 < PRIMS[RHOB] =rho_b_atm; < if(gamma_equals2) { < PRIMS[PRESSURE]=K_poly*SQR(rho_b_atm); < } else { < PRIMS[PRESSURE]=K_poly*pow(rho_b_atm,gamma_th); < } --- > PRIMS[RHOB] = rho_b_atm; > > /* Set P = P_cold */ > int polytropic_index = find_polytropic_K_and_Gamma_index(eos, rho_b_atm); > CCTK_REAL K_ppoly_tab = eos.K_ppoly_tab[polytropic_index]; > CCTK_REAL Gamma_ppoly_tab = eos.K_ppoly_tab[polytropic_index]; > PRIMS[PRESSURE] = K_ppoly_tab*pow(rho_b_atm,Gamma_ppoly_tab); > 257a278 > 261a283 > 274a297 > 304c327 < error_int_numer += fabs(tau[index] - tau_orig) + fabs(rho_star[index] - rho_star_orig) + --- > error_int_numer += fabs(tau[index] - tau_orig) + fabs(rho_star[index] - rho_star_orig) + 341c364 < --- > 346c369 < // Very useful con2prim debugger. If the primitives (con2prim) solver fails, this will output all data needed to --- > // Very useful con2prim debugger. If the primitives (con2prim) solver fails, this will output all data needed to 348c371 < if(conserv_to_prims_debug==1) { --- > if(conserv_to_prims_debug==1 && error_int_numer/error_int_denom > 0.05) { 353,354c376,379 < sprintf(filename,"primitives_debug-%e.dat",rho_star[CCTK_GFINDEX3D(cctkGH,3,15,6)]); < //Alternative, for debugging purposes as well: sprintf(filename,"primitives_debug-%d.dat",rand()); --- > sprintf(filename,"primitives_debug-%e.dat",error_int_numer/error_int_denom); > //Alternative, for debugging purposes as well: > //srand(time(NULL)); > //sprintf(filename,"primitives_debug-%d.dat",rand()); 358a384,386 > myfile.write((char*)&GAMMA_SPEED_LIMIT, 1*sizeof(CCTK_REAL)); > > myfile.write((char*)&rho_b_max, 1*sizeof(CCTK_REAL)); 364,373c392,398 < CCTK_REAL gamma_th=2.0; < myfile.write((char*)&gamma_th, 1*sizeof(CCTK_REAL)); < myfile.write((char*)&neos, 1*sizeof(int)); < < myfile.write((char*)gamma_tab, (neos+1)*sizeof(CCTK_REAL)); < myfile.write((char*)k_tab, (neos+1)*sizeof(CCTK_REAL)); < < myfile.write((char*)eps_tab, neos*sizeof(CCTK_REAL)); < myfile.write((char*)rho_tab, neos*sizeof(CCTK_REAL)); < myfile.write((char*)P_tab, neos*sizeof(CCTK_REAL)); --- > myfile.write((char*)&update_Tmunu, 1*sizeof(int)); > > myfile.write((char*)&neos, 1*sizeof(int)); > myfile.write((char*)&Gamma_th, 1*sizeof(CCTK_REAL)); > myfile.write((char*)&K_ppoly_tab0, 1*sizeof(CCTK_REAL)); > myfile.write((char*)Gamma_ppoly_tab_in, neos*sizeof(CCTK_REAL)); > myfile.write((char*)rho_ppoly_tab_in, (neos-1)*sizeof(CCTK_REAL)); 378a404,424 > > myfile.write((char *)failure_checker, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTtt, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTtx, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTty, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTtz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTxx, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTxy, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTxz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTyy, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTyz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)eTzz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)alp, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)gxx, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)gxy, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)gxz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)gyy, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)gyz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)gzz, fullsize*sizeof(CCTK_REAL)); > myfile.write((char *)psi_bssn, fullsize*sizeof(CCTK_REAL)); > 399c445 < --- > 421c467 < CCTK_VInfo(CCTK_THORNSTRING,"Finished writing..."); --- > CCTK_VInfo(CCTK_THORNSTRING,"Finished writing %s",filename); 422a469,473 > > #ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER > return 0; // int main() requires an integer be returned > #endif > 425a477 >
# 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_primitives_lowlevel.C"
original_IGM_file_name = "harm_primitives_lowlevel-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_primitives_lowlevel__C = !diff $original_IGM_file_path $outfile_path__harm_primitives_lowlevel__C
if Validation__harm_primitives_lowlevel__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_primitives_lowlevel.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for harm_primitives_lowlevel.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_primitives_lowlevel__C:
print(diff_line)
Validation test for harm_primitives_lowlevel.C: FAILED! Diff: 6c6 < --- > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER 7a8 > #endif 9,12d9 < CCTK_REAL kpoly = eos.k_tab[0]; < CCTK_REAL gamma = eos.gamma_tab[0]; < int gamma_equals2 = 1; < if (fabs(gamma-2.0) > 1.e-10) gamma_equals2 = 0; 15c12 < CCTK_REAL U[NPR]; --- > CCTK_REAL U[NPR]; 42a40 > 46c44 < // This is NOT the \mathcal{B}^i, which differs by --- > // This is NOT the \mathcal{B}^i, which differs by 51a50 > 57a57 > 60,61c60,61 < between the two sets of definitions for U and P. The user may < wish to alter the translation as they see fit. --- > between the two sets of definitions for U and P. The user may > wish to alter the translation as they see fit. 64,72c64,72 < // / 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 / // 76c76 < --- > 81c81 < // prim[NPR] = primitive variables (guess on input, calculated values on --- > // prim[NPR] = primitive variables (guess on input, calculated values on 84c84 < // U[1] = --- > // U[1] = 92a93 > 103a105 > CCTK_REAL K_ppoly_tab,Gamma_ppoly_tab; 111,115c113,130 < if (gamma_equals2==1) { < P_oldL = kpoly*rho_b_oldL*rho_b_oldL; < } else { < P_oldL = kpoly*pow(rho_b_oldL,gamma); < } --- > > /********************************** > * Piecewise Polytropic EOS Patch * > * Finding Gamma_ppoly_tab and K_ppoly_tab * > **********************************/ > /* Here we use our newly implemented > * find_polytropic_K_and_Gamma() function > * to determine the relevant polytropic > * Gamma and K parameters to be used > * within this function. > */ > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL); > K_ppoly_tab = eos.K_ppoly_tab[polytropic_index]; > Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index]; > > // After that, we compute P_cold > P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab); > 125,130c140,157 < // GAMMA=2 ONLY: < if (gamma_equals2==1) { < P_oldL = kpoly*rho_b_oldL*rho_b_oldL; < } else { < P_oldL = kpoly*pow(rho_b_oldL,gamma); < } --- > > /********************************** > * Piecewise Polytropic EOS Patch * > * Finding Gamma_ppoly_tab and K_ppoly_tab * > **********************************/ > /* Here we use our newly implemented > * find_polytropic_K_and_Gamma() function > * to determine the relevant polytropic > * Gamma and K parameters to be used > * within this function. > */ > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL); > K_ppoly_tab = eos.K_ppoly_tab[polytropic_index]; > Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index]; > > // After that, we compute P_cold > P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab); > 136a164 > 139c167 < U[UU] = -CONSERVS[TAUENERGY]*METRIC_LAP_PSI4[LAPSE] - (METRIC_LAP_PSI4[LAPSE]-1.0)*CONSERVS[RHOSTAR] + --- > U[UU] = -CONSERVS[TAUENERGY]*METRIC_LAP_PSI4[LAPSE] - (METRIC_LAP_PSI4[LAPSE]-1.0)*CONSERVS[RHOSTAR] + 148c176,177 < CCTK_REAL uL = P_oldL/(gamma_th /* <- Should be local polytropic gamma factor */ - 1.0); --- > > CCTK_REAL uL = P_oldL/(Gamma_ppoly_tab - 1.0); 161a191 > 164,166c194,196 < check = Utoprim_2d(U, g4dn, g4up, detg, prim,stats.n_iter); < // Note that we have modified this solver, so that nearly 100% < // of the time it yields either a good root, or a root with --- > check = Utoprim_2d(eos, U, g4dn, g4up, detg, prim,stats.n_iter); > // Note that we have modified this solver, so that nearly 100% > // of the time it yields either a good root, or a root with 170c200,201 < // Use the new Font fix subroutine --- > > // Use the new Font fix subroutine 175,179c206,210 < if (gamma_equals2==1) { < check = font_fix_gamma_equals2(u_xl,u_yl,u_zl,CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4,eos); < } else { < check = font_fix_general_gamma(u_xl,u_yl,u_zl,CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4,eos); < } --- > /************************ > * New Font fix routine * > ************************/ > check = font_fix__hybrid_EOS(u_xl,u_yl,u_zl, CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4, eos); > 192a224 > 223a256 > 229,234c262,280 < if (gamma_equals2==1) { < P_cold = kpoly*prim[RHO]*prim[RHO]; // <- FIXME: GAMMA=2 ONLY < } else { < P_cold = kpoly*pow(prim[RHO],gamma); < } < prim[UU] = P_cold/(gamma_th /* <- Should be local polytropic gamma factor */-1.0); --- > > /********************************** > * Piecewise Polytropic EOS Patch * > * Finding Gamma_ppoly_tab and K_ppoly_tab * > **********************************/ > /* Here we use our newly implemented > * find_polytropic_K_and_Gamma() function > * to determine the relevant polytropic > * Gamma and K parameters to be used > * within this function. > */ > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,prim[RHO]); > K_ppoly_tab = eos.K_ppoly_tab[polytropic_index]; > Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index]; > > // After that, we compute P_cold > P_cold = K_ppoly_tab*pow(prim[RHO],Gamma_ppoly_tab); > > prim[UU] = P_cold/(Gamma_ppoly_tab-1.0); 237,239c283,301 < PRIMS[RHOB] = prim[RHO]; < PRIMS[PRESSURE] = (gamma_th /* <- Should be local polytropic gamma factor */-1.0)*prim[UU]; < //Already set u0L. --- > /* Set rho_b */ > PRIMS[RHOB] = prim[RHO]; > > /*************** > * PPEOS Patch * > * Hybrid EOS * > *************** > */ > /* We now compute the pressure as a function > * of rhob, P_cold, eps_cold, and u = rhob*eps, > * using the function pressure_rho0_u(), which > * implements the equation: > * .-------------------------------------------------------------. > * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) | > * .-------------------------------------------------------------. > */ > PRIMS[PRESSURE] = pressure_rho0_u(eos, prim[RHO],prim[UU]); > > /* Already set u0L. */ 254c316 < #include "harm_u2p_util.c" --- > //#include "harm_u2p_util.c" 257c319,320 < #include "font_fix_gamma_law.C" --- > #include "font_fix_hybrid_EOS.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/font_fix_gamma_law.C"
original_IGM_file_name = "font_fix_gamma_law-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__font_fix_gamma_law__C = !diff $original_IGM_file_path $outfile_path__font_fix_hybrid_EOS__C
if Validation__font_fix_gamma_law__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 font_fix_gamma_law.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for font_fix_gamma_law.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__font_fix_gamma_law__C:
print(diff_line)
Validation test for font_fix_gamma_law.C: FAILED! Diff: 1,10d0 < /*----------------------------------------------------------------------------- < * < * Apply "Font Fix", of Font et al., ONLY when the primitives (con2prim) solver < * fails. Note that this version is OPTIMIZED for gamma=2 polytrope EOS only. < * Application of this fix GUARANTEES (on an analytic level, at least) < * that con2prim will succeed. < * < * FIXME: There exist better fixes/strategies now. See McKinney et al's work. < * < -----------------------------------------------------------------------------*/ 12,39d1 < inline int font_fix_gamma_equals2(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos) { < CCTK_REAL rhob; < CCTK_REAL kpoly2 = 2.0*eos.K_poly; < CCTK_REAL tol = 1.e-15; < CCTK_REAL Bxbar = PRIMS[BX_CENTER]*ONE_OVER_SQRT_4PI; < CCTK_REAL Bybar = PRIMS[BY_CENTER]*ONE_OVER_SQRT_4PI; < CCTK_REAL Bzbar = PRIMS[BZ_CENTER]*ONE_OVER_SQRT_4PI; < CCTK_REAL Bbar_x = METRIC_PHYS[GXX]*Bxbar + METRIC_PHYS[GXY]*Bybar + METRIC_PHYS[GXZ]*Bzbar; < CCTK_REAL Bbar_y = METRIC_PHYS[GXY]*Bxbar + METRIC_PHYS[GYY]*Bybar + METRIC_PHYS[GYZ]*Bzbar; < CCTK_REAL Bbar_z = METRIC_PHYS[GXZ]*Bxbar + METRIC_PHYS[GYZ]*Bybar + METRIC_PHYS[GZZ]*Bzbar; < CCTK_REAL B2bar = Bxbar*Bbar_x + Bybar*Bbar_y + Bzbar*Bbar_z; < CCTK_REAL Bbar = sqrt(B2bar); < CCTK_REAL check_B_small = fabs(Bxbar)+fabs(Bybar)+fabs(Bzbar); < if (check_B_small>0 && check_B_small<1.e-150) { < // need to compute B2bar specially to prevent floating-point underflow < CCTK_REAL Bmax = fabs(Bxbar); < if (Bmax < fabs(Bybar)) Bmax=fabs(Bybar); < if (Bmax < fabs(Bzbar)) Bmax=fabs(Bzbar); < CCTK_REAL Bxtmp=Bxbar/Bmax, Bytemp=Bybar/Bmax, Bztemp=Bzbar/Bmax; < CCTK_REAL B_xtemp=Bbar_x/Bmax, B_ytemp=Bbar_y/Bmax, B_ztemp=Bbar_z/Bmax; < Bbar = sqrt(Bxtmp*B_xtemp + Bytemp*B_ytemp + Bztemp*B_ztemp)*Bmax; < } < < CCTK_REAL BbardotS = Bxbar*CONSERVS[STILDEX] + Bybar*CONSERVS[STILDEY] + Bzbar*CONSERVS[STILDEZ]; < CCTK_REAL BbardotS2 = BbardotS*BbardotS; < CCTK_REAL hatBbardotS = BbardotS/Bbar; < if (Bbar<1.e-300) hatBbardotS = 0.0; < CCTK_REAL Psim6 = 1.0/METRIC_LAP_PSI4[PSI6]; 41,165c3,7 < // Limit hatBbardotS < //CCTK_REAL max_gammav = 100.0; < //CCTK_REAL rhob_max = CONSERVS[RHOSTAR]*Psim6; < //CCTK_REAL hmax = 1.0 + kpoly2*rhob_max; < //CCTK_REAL abs_hatBbardotS_max = sqrt(SQR(max_gammav)-1.0)*CONSERVS[RHOSTAR]*hmax; < //if (fabs(hatBbardotS) > abs_hatBbardotS_max) { < // CCTK_REAL fac_reduce = abs_hatBbardotS_max/fabs(hatBbardotS); < // CCTK_REAL hatBbardotS_max = hatBbardotS*fac_reduce; < // CCTK_REAL Bbar_inv = 1.0/Bbar; < // CCTK_REAL hat_Bbar_x = Bbar_x*Bbar_inv; < // CCTK_REAL hat_Bbar_y = Bbar_y*Bbar_inv; < // CCTK_REAL hat_Bbar_z = Bbar_z*Bbar_inv; < // CCTK_REAL sub_fact = hatBbardotS_max - hatBbardotS; < // CONSERVS[STILDEX] += sub_fact*hat_Bbar_x; < // CONSERVS[STILDEY] += sub_fact*hat_Bbar_y; < // CONSERVS[STILDEZ] += sub_fact*hat_Bbar_z; < // hatBbardotS = hatBbardotS_max; < // BbardotS *= fac_reduce; < // BbardotS2 = BbardotS*BbardotS; < //} < < CCTK_REAL sdots = METRIC_PHYS[GUPXX]*SQR(CONSERVS[STILDEX]) + METRIC_PHYS[GUPYY]*SQR(CONSERVS[STILDEY]) + METRIC_PHYS[GUPZZ]*SQR(CONSERVS[STILDEZ]) < + 2.0*( METRIC_PHYS[GUPXY]*CONSERVS[STILDEX]*CONSERVS[STILDEY] + METRIC_PHYS[GUPXZ]*CONSERVS[STILDEX]*CONSERVS[STILDEZ] < + METRIC_PHYS[GUPYZ]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]); < < if (sdots<1.e-300) { < rhob = CONSERVS[RHOSTAR]*Psim6; < u_x=0.0; u_y=0.0; u_z=0.0; < return 0; < } < < /* < if (fabs(BbardotS2 - sdots*B2bar) > 1e-8) { < CCTK_VInfo(CCTK_THORNSTRING,"(Bbar dot S)^2, Bbar^2 * sdotS, %e %e",SQR(BbardotS),sdots*B2bar); < CCTK_VInfo(CCTK_THORNSTRING,"Cauchy-Schwartz inequality is violated!"); < } < */ < < // Initial guess for W, S_fluid and rhob < CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6; < CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar); < CCTK_REAL rhob0 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR])); < CCTK_REAL W=W0,Sf2=Sf20,rhob1=rhob0; < < //**************************************************************** < // FONT FIX < // Impose Font fix when HARM primitives solver fails to find < // acceptable set of primitives. < //**************************************************************** < bool fontcheck=true; < < int itcount = 0, maxits=300; < while(fontcheck && itcount < maxits) { < itcount++; < W0 = W; < Sf20 = Sf2; < rhob0 = rhob1; < < // first, find rhob for the given S_fluid^2 < CCTK_REAL h = 1.0 + kpoly2*rhob0; < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < while( fabs(rhob1-rhob0) > rhob1*tol) { < rhob0 = rhob1; < h = 1.0 + kpoly2*rhob0; < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < } < < h = 1.0 + kpoly2*rhob1; < W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6; < Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar); < if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false; < } < < if (itcount>=maxits) { < // Increase tol and try again < maxits*=100; < tol *=10.0; < itcount = 0; < fontcheck=true; < while(fontcheck && itcount < maxits) { < itcount++; < W0 = W; < Sf20 = Sf2; < rhob0 = rhob1; < < // first, find rhob for the given S_fluid^2 < CCTK_REAL h = 1.0 + kpoly2*rhob0; < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < while( fabs(rhob1-rhob0) > rhob1*tol) { < rhob0 = rhob1; < h = 1.0 + kpoly2*rhob0; < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < } < < h = 1.0 + kpoly2*rhob1; < W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6; < Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar); < if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false; < } < } < //************************************************************************************************************** < < if(fontcheck==true) { < return 1; < } < < // Font fix works, now compute u_i < rhob = rhob1; < CCTK_REAL h = 1.0+kpoly2*rhob; < CCTK_REAL gammav = CONSERVS[RHOSTAR]*Psim6/rhob; < CCTK_REAL rhosh = CONSERVS[RHOSTAR]*h; < CCTK_REAL fac1 = METRIC_LAP_PSI4[PSI6]*BbardotS/(gammav*rhosh); < CCTK_REAL fac2 = 1.0/(rhosh + METRIC_LAP_PSI4[PSI6]*B2bar/gammav); < u_x = fac2*(CONSERVS[STILDEX] + fac1*Bbar_x); < u_y = fac2*(CONSERVS[STILDEY] + fac1*Bbar_y); < u_z = fac2*(CONSERVS[STILDEZ] + fac1*Bbar_z); < return 0; < } < < inline int font_fix_general_gamma(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos) { < CCTK_REAL gamma=eos.gamma_tab[0]; < CCTK_REAL rhob; < CCTK_REAL tol = 1.e-15; < CCTK_REAL gam1 = gamma-1.0; < CCTK_REAL gam_gamm1_kpoly = gamma/gam1*eos.K_poly; --- > /********************************** > * Piecewise Polytropic EOS Patch * > * Font fix: function call * > **********************************/ > inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos) { 216a59,60 > > CCTK_REAL rhob; 227a72,73 > > 229,230c75,76 < CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6; < CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar); --- > CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6; > CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar); 232c78 < CCTK_REAL W=W0,Sf2=Sf20,rhob1=rhob0; --- > 240c86,87 < bool fontcheck=true; --- > /* Set the maximum number of iterations */ > int maxits = 500; 242,262c89,90 < int itcount = 0, maxits=500; < while(fontcheck && itcount < maxits) { < itcount++; < W0 = W; < Sf20 = Sf2; < rhob0 = rhob1; < < // first, find rhob for the given S_fluid^2 < CCTK_REAL h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1); < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < while( fabs(rhob1-rhob0) > rhob1*tol) { < rhob0 = rhob1; < h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1); < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < } < < h = 1.0 + gam_gamm1_kpoly*pow(rhob1,gam1); < W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6; < Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar); < if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false; < } --- > /* Set the allowed tolerance */ > CCTK_REAL tol = 1.e-15; 264,290c92,110 < if (itcount>=maxits) { < // Increase tol and try again < fontcheck=true; < tol *=100.0; < itcount = 0; < while(fontcheck && itcount < maxits) { < itcount++; < W0 = W; < Sf20 = Sf2; < rhob0 = rhob1; < < // first, find rhob for the given S_fluid^2 < CCTK_REAL h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1); < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < while( fabs(rhob1-rhob0) > rhob1*tol) { < rhob0 = rhob1; < h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1); < rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h)); < } < < h = 1.0 + gam_gamm1_kpoly*pow(rhob1,gam1); < W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6; < Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar); < if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false; < } < } < //************************************************************************************************************** --- > /* Declare basic variables */ > int font_fix_status; > > /********************** > * FONT FIX MAIN LOOP * > ********************** > * Perform the font fix routine until convergence > * is obtained and the algorithm returns with no > * error. Every time the Font fix fails, increase > * the tolerance by a factor of 10. > */ > int font_fix_attempts = 5; > CCTK_REAL font_fix_tol_factor = 10.0; > for(int n=0; n<font_fix_attempts; n++) { > > tol *= pow(font_fix_tol_factor,n); > font_fix_status = font_fix__rhob_loop(maxits,tol, W0,Sf20,Psim6,sdots,BbardotS2,B2bar, CONSERVS,eos, rhob0,rhob); > rhob0 = rhob; > if(font_fix_status==0) break; 292,293d111 < if(fontcheck==true) { < return 1; 296,298c114,128 < // Font fix works, now compute u_i < rhob = rhob1; < CCTK_REAL h = 1.0+gam_gamm1_kpoly*pow(rhob,gam1); --- > > //************************************************************************************************************** > > /* Font fix works! */ > /* First compute P_cold, eps_cold, then h = h_cold */ > CCTK_REAL P_cold, eps_cold; > compute_P_cold__eps_cold(eos,rhob, P_cold,eps_cold); > CCTK_REAL h = 1.0 + eps_cold + P_cold/rhob; > > /* Then compute gamma_v using equation (A19) in > * Etienne et al. (2011) [https://arxiv.org/pdf/1112.0568.pdf] > * .-----------------------------------------. > * | gamma_v = psi^{-6} * (rho_star / rho_b) | > * .-----------------------------------------. > */ 299a130,131 > > /* Finally, compute u_{i} */
# 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_primitives_headers.h"
original_IGM_file_name = "harm_primitives_headers-original.h"
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_primitives_headers__h = !diff $original_IGM_file_path $outfile_path__harm_primitives_headers__h
if Validation__harm_primitives_headers__h == []:
# 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_primitives_headers.h: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for harm_primitives_headers.h: 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_primitives_headers__h:
print(diff_line)
Validation test for harm_primitives_headers.h: FAILED! Diff: 2c2 < Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble, --- > Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble, 7c7 < 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,10c9,10 < 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 12c12 < an accretion disk model. --- > an accretion disk model. 14c14 < 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 17c17 < [1] Gammie, C. F., McKinney, J. C., \& Toth, G.\ 2003, --- > [1] Gammie, C. F., McKinney, J. C., \& Toth, G.\ 2003, 20c20 < [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,24c23,24 < < Further, we strongly encourage you to obtain the latest version of --- > > Further, we strongly encourage you to obtain the latest version of 59c59 < #if( USE_ISENTROPIC ) --- > #if( USE_ISENTROPIC ) 89c89 < static const int RHO =0; --- > static const int RHO =0; 106c106 < 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], 115,117c115 < inline int font_fix_general_gamma(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos); < inline int font_fix_gamma_equals2(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos); < --- > inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos); 123a122 >
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__driver_conserv_to_prims.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__the_conservative_to_primitive_algorithm.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.tex
!rm -f Tut*.out Tut*.aux Tut*.log