This module is currently under development
This module is organized as follows
reconstruct_set_of_prims_PPM()
function
slope_limit()
functionsteepen_rho()
functionmonotonize()
functioncompute_P_cold__Gamma_cold()
functionftilde_gf_compute()
functionftilde_compute()
functionWe will now use the cmdline_helper.py NRPy+ module to create the source directory within the IllinoisGRMHD
NRPy+ directory, if it does not exist yet.
# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
# Step 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)
# Step 0c: Create the output file path
outfile_path__reconstruct_set_of_prims_PPM__C = os.path.join(IGM_src_dir_path,"reconstruct_set_of_prims_PPM.C")
outfile_path__loop_defines_reconstruction__h = os.path.join(IGM_src_dir_path,"loop_defines_reconstruction.h")
In this tutorial module, we will go through the implementation of the piecewise parabolic method (PPM), introduced by Colella & Woodward (1984) (which shall henceforth be our main reference), used by IllinoisGRMHD
.
The piecewise parabolic method (PPM) is an algorithm used to construct the values of primitive variables, $U$, at cell interfaces. The interpolation procedure alone can lead to unstable evolutions. To remedy this, we introduce three different techniques:
These algorithms are intended to also produce narrower profiles near the vicinity of a shock. These steps tend to reduce the third-order accuracy of the interpolation code, but only in cases where a third-order interpolation algorithm would produce worse results (e.g. at local extrema).
The algorithmic flow of the code is as follows:
where $\delta U^{\rm slope-lim}$ is referred to as the slope-limited gradient of $U$ and
\begin{align} \delta U_{i} &\equiv U_{i} - U_{i-1}\ ,\\ \delta_{m} U_{i} &\equiv \frac{\delta U_{i} + \delta U_{i+1}}{2} = \frac{U_{i+1} - U_{i-1}}{2}\ . \end{align}with $K_{0}$ a problem dependent constant.
where
\begin{align} \rho^{\rm MC}_{r,i+1} &= \rho_{i+1} - \frac{1}{2}\delta\rho^{\rm slope-lim}_{i+1}\ ,\\ \rho^{\rm MC}_{l,i+0} &= \rho_{i-1} + \frac{1}{2}\delta\rho^{\rm slope-lim}_{i-1}\ ,\\ \eta_{i} &= \max\left\{0,\min\left[\eta_{1}\left(\tilde\eta_{i}-\eta_{2}\right),1\right]\right\}\ , \end{align}with $\eta_{1}$ and $\eta_{2}$ constants and
$$ \tilde\eta_{i} = \left\{ \begin{matrix} 0&, \ {\rm if}\ \delta\rho_{i} = 0\ ,\\ -\frac{1}{6}\left(\frac{\delta^{2}\rho_{i+1} - \delta^{2}\rho_{i-1}}{2\delta\rho_{i}}\right)&, \ {\rm otherwise}\ , \end{matrix} \right. $$and finally
\begin{align} \delta\rho_{i+0} &= \frac{\rho_{i+1} - \rho_{i-1}}{2}\ ,\\ \delta^{2}\rho_{i-1} &= \rho_{i+0} - 2\rho_{i-1} + \rho_{i-2}\ ,\\ \delta^{2}\rho_{i+1} &= \rho_{i+2} - 2\rho_{i+1} + \rho_{i+0}\ . \end{align}$^{\color{red}\ddagger}$: note that the common notation is $\rho_{0}$, but we will use $\rho_{b}$ to represent the baryonic matter density.
where
\begin{align} \tilde{f} &= \min\left[1,w\max\left(0,q_{1}\right)\right]\ ,\\ w &= \left\{ \begin{matrix} 1\ , &\ {\rm if}\ q_{2} > \epsilon_{2}\ {\rm and}\ q_{2}\left(v^{\rm flux\ dirn}_{i-1}-v^{\rm flux\ dirn}_{i+1}\right)>0\ \left({\rm inside\ shock}\right)\ ,\\ 0\ , &\ {\rm otherwise}\ \left({\rm outside\ shock}\right)\ , \end{matrix} \right. \end{align}and
\begin{align} q_{1} &= \left(\frac{\delta P_{1}}{\delta P_{2}}-\omega_{1}\right)\omega_{2}\ ,\\ q_{2} &= \frac{\left|\delta P_{1}\right|}{\min\left(P_{i+1},P_{i-1}\right)}\ , \end{align}with $\omega_{1}$ and $\omega_{2}$ constants and $\delta P_{n} \equiv P_{i+n} - P_{i-n}$.
where
\begin{align} \delta U &\equiv U_{r} - U_{l}\ ,\\ \delta_{m}U &\equiv \frac{U_{r} + U_{l}}{2}\ . \end{align}We then perform the following shift
$$ \boxed{ \begin{matrix} U_{i-1/2+\epsilon} = U_{l,i}^{\rm old} = U_{r,i}^{\rm new}\ ,\\ U_{i-1/2-\epsilon} = U_{r,i-1}^{\rm old} = U_{l,i}^{\rm new}\ .\\ \end{matrix} } $$loop_defines_reconstruction.h
header file [Back to top]¶This header file defines useful quantities to be used throughout the reconstruct_set_of_prims_PPM.C
code. They are:
%%writefile $outfile_path__loop_defines_reconstruction__h
#ifndef LOOP_DEFINES_RECONSTRUCTION_H_
#define LOOP_DEFINES_RECONSTRUCTION_H_
#define LOOP_DEFINE(gz_shift_lo,gz_shift_hi, ext,flux_dirn, ijkgz_lo_hi,gz_lo,gz_hi) \
for(int rr=1;rr<=3;rr++) { \
ijkgz_lo_hi[rr][0]= gz_lo[rr]; \
ijkgz_lo_hi[rr][1]=ext[rr-1]-gz_hi[rr]; \
} \
ijkgz_lo_hi[flux_dirn][0] += gz_shift_lo; \
ijkgz_lo_hi[flux_dirn][1] -= gz_shift_hi; \
/* The following line is valid C99 */ \
_Pragma("omp parallel for private(U,dU,slope_lim_dU,Ur,Ul)") \
for(int k=ijkgz_lo_hi[3][0];k<ijkgz_lo_hi[3][1];k++) \
for(int j=ijkgz_lo_hi[2][0];j<ijkgz_lo_hi[2][1];j++) \
for(int i=ijkgz_lo_hi[1][0];i<ijkgz_lo_hi[1][1];i++)
// This define only sets indices.
// FIXME: benchmark with and without the if() statement.
// FIXME: try without index_arr being defined in all directions.
#define SET_INDEX_ARRAYS(IMIN,IMAX,flux_dirn) \
int max_shift=(MAXNUMINDICES/2); \
/* DEBUGGING ONLY: if(IMIN<-max_shift || IMAX>max_shift) CCTK_VError(VERR_DEF_PARAMS,"FIX MAXNUMINDICES!"); */ \
int index_arr[4][MAXNUMINDICES]; \
for(int idx=IMIN;idx<=IMAX;idx++) { \
index_arr[flux_dirn][idx+max_shift]= \
CCTK_GFINDEX3D(cctkGH, \
i+idx*kronecker_delta[flux_dirn][0], \
j+idx*kronecker_delta[flux_dirn][1], \
k+idx*kronecker_delta[flux_dirn][2]); \
}
#define SET_INDEX_ARRAYS_3DBLOCK(IJKLOHI) \
int max_shift=(MAXNUMINDICES/2); \
int index_arr_3DB[MAXNUMINDICES][MAXNUMINDICES][MAXNUMINDICES]; \
for(int idx_k=IJKLOHI[4];idx_k<=IJKLOHI[5];idx_k++) for(int idx_j=IJKLOHI[2];idx_j<=IJKLOHI[3];idx_j++) for(int idx_i=IJKLOHI[0];idx_i<=IJKLOHI[1];idx_i++) { \
index_arr_3DB[idx_k+max_shift][idx_j+max_shift][idx_i+max_shift]=CCTK_GFINDEX3D(cctkGH,i+idx_i,j+idx_j,k+idx_k); \
}
#endif /* LOOP_DEFINES_RECONSTRUCTION_H_ */
Writing ../src/loop_defines_reconstruction.h
reconstruct_set_of_prims_PPM.C
code file [Back to top]¶We then initialize the reconstruct_set_of_prims_PPM.C
code file with a basic preamble, some simple definitions, and function headers. Notice that these functions are defined in the other Steps of this tutorial notebook. See the table of contents above for more information.
%%writefile $outfile_path__reconstruct_set_of_prims_PPM__C
/*****************************************
* PPM Reconstruction Interface.
* Zachariah B. Etienne (2013)
*
* This version of PPM implements the standard
* Colella & Woodward PPM, though modified as in GRHydro
* to have 3 ghostzones instead of 4.
*****************************************/
#define MINUS2 0
#define MINUS1 1
#define PLUS0 2
#define PLUS1 3
#define PLUS2 4
#define MAXNUMINDICES 5
// ^^^^^^^^^^^^^ Be _sure_ to define MAXNUMINDICES appropriately!
// You'll find the #define's for LOOP_DEFINE and SET_INDEX_ARRAYS inside:
#include "loop_defines_reconstruction.h"
static inline CCTK_REAL ftilde_compute(const int flux_dirn,CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES]);
static inline CCTK_REAL slope_limit(CCTK_REAL dU,CCTK_REAL dUp1);
static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],
CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold,
CCTK_REAL *rho_br_ppm,CCTK_REAL *rho_bl_ppm);
static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold);
static inline void monotonize(CCTK_REAL U,CCTK_REAL &Ur,CCTK_REAL &Ul);
Writing ../src/reconstruct_set_of_prims_PPM.C
reconstruct_set_of_prims_PPM()
function [Back to top]¶The reconstruct_set_of_prims_PPM()
function receives as input the direction in which to perform the reconstruction ($\rm flux\_dirn$), the number of primitives which are to be reconstructed ($\rm num\_prims\_to\_reconstruct$), which primitives are to be reconstructed ($\rm which\_prims\_to\_reconstruct$), the equation of state ($\rm eos$), and the array which stores the primitives ($\rm in\_prims$). The reconstructed primitive values are then stored in the output arrays $\rm out\_prims\_r$ and $\rm out\_prims\_l$.
Note: you can find more information on the ETK-specific parameters (such as ${\rm cctkGH}$ and $\rm cctk\_lsh$) by looking at the Cactus Variables paragraph of section C1.6.2 of the Einstein Toolkit UserGuide.
Notice that we will start a loop that runs from 0 to $\rm num\_prims\_to\_reconstruct$, meaning that the discussion here will apply to each of the primitives that we wish to reconstruct.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,eos_struct &eos,
gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,CCTK_REAL *ftilde_gf, CCTK_REAL *temporary) {
DECLARE_CCTK_PARAMETERS;
CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],dU[MAXNUMVARS][MAXNUMINDICES],slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],
Ur[MAXNUMVARS][MAXNUMINDICES],Ul[MAXNUMVARS][MAXNUMINDICES];
int ijkgz_lo_hi[4][2];
for(int ww=0;ww<num_prims_to_reconstruct;ww++) {
int whichvar=which_prims_to_reconstruct[ww];
if(in_prims[whichvar].gz_lo[flux_dirn]!=0 || in_prims[whichvar].gz_hi[flux_dirn]!=0) {
CCTK_VError(VERR_DEF_PARAMS,"TOO MANY GZ'S! WHICHVAR=%d: %d %d %d : %d %d %d DIRECTION %d",whichvar,
in_prims[whichvar].gz_lo[1],in_prims[whichvar].gz_lo[2],in_prims[whichvar].gz_lo[3],
in_prims[whichvar].gz_hi[1],in_prims[whichvar].gz_hi[2],in_prims[whichvar].gz_hi[3],flux_dirn);
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
We start by reading in the input primitive gridfunction and storing them to a variable $U$. Notice that for a given direction and a given point $i$, we will know: $\left\{U_{i-2},U_{i-1},U_{i},U_{i+1},U_{i+2}\right\}$.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// *** LOOP 1: Interpolate to Ur and Ul, which are face values ***
// You will find that Ur depends on U at MINUS1,PLUS0, PLUS1,PLUS2, and
// Ul depends on U at MINUS2,MINUS1,PLUS0,PLUS1.
// However, we define the below loop from MINUS2 to PLUS2. Why not split
// this up and get additional points? The reason is that later on,
// Ur and Ul depend on ftilde, which is defined from MINUS2 to PLUS2,
// so we would lose those points anyway.
LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
SET_INDEX_ARRAYS(-2,2,flux_dirn);
/* *** LOOP 1a: READ INPUT *** */
// Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U.
for(int ii=MINUS2;ii<=PLUS2;ii++) U[whichvar][ii] = in_prims[whichvar].gf[index_arr[flux_dirn][ii]];
Appending to ../src/reconstruct_set_of_prims_PPM.C
We will need $\delta U_{i} \equiv U_{i} - U_{i-1}$ in order to compute $U_{r}$ and $U_{l}$. We will then compute (notice the notation change $i\to i+0$ so that it is easier to understand the C code that follows):
\begin{align} \delta U_{i-1} &= U_{i-1} - U_{i-2}\ ,\\ \delta U_{i+0} &= U_{i+0} - U_{i-1}\ ,\\ \delta U_{i+1} &= U_{i+1} - U_{i+0} \ ,\\ \delta U_{i+2} &= U_{i+2} - U_{i+1}\ . \end{align}After evaluating the differences $\delta U$, we compute the slope-limited $\delta U$, $\delta U^{\rm slope-lim}$, using the slope_limit()
function.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
/* *** LOOP 1b: DO COMPUTATION *** */
/* First, compute simple dU = U(i) - U(i-1), where direction of i
* is given by flux_dirn, and U is a primitive variable:
* {rho_b,P,vx,vy,vz,Bx,By,Bz}. */
// Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1),
// and dU(i+2)
dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];
dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1];
dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0];
dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1];
// Then, compute slope-limited dU, using MC slope limiter:
slope_lim_dU[whichvar][MINUS1]=slope_limit(dU[whichvar][MINUS1],dU[whichvar][PLUS0]);
slope_lim_dU[whichvar][PLUS0] =slope_limit(dU[whichvar][PLUS0], dU[whichvar][PLUS1]);
slope_lim_dU[whichvar][PLUS1] =slope_limit(dU[whichvar][PLUS1], dU[whichvar][PLUS2]);
Appending to ../src/reconstruct_set_of_prims_PPM.C
We now compute $U_{r}$ and $U_{l}$. Keep in mind that $U_{r,i} = U_{i+1/2}$, while $U_{l,i} = U_{i-1/2}$. The implemented equation follows eq. A1 in Duez et al. (2005), but with the standardd PPM coefficient of $\frac{1}{6}$ (i.e. eq. A1 with $\frac{1}{8}\to\frac{1}{6}$). Keep in mind that we simplify the equation slightly before implementing it:
\begin{align} U_{r,i+0} &= U_{i+0} + \frac{1}{2}\left(U_{i+1} - U_{i+0}\right) + \frac{1}{6}\left(\delta U^{\rm slope-lim}_{i+0} - \delta U^{\rm slope-lim}_{i+1}\right) \Rightarrow \boxed{U_{r,i+0} = \frac{1}{2}\left(U_{i+1} + U_{i+0}\right) + \frac{1}{6}\left(\delta U^{\rm slope-lim}_{i+0} - \delta U^{\rm slope-lim}_{i+1}\right)}\ ,\\ U_{l,i+0} &= U_{i-1} + \frac{1}{2}\left(U_{i+0} - U_{i-1}\right) + \frac{1}{6}\left(\delta U^{\rm slope-lim}_{i-1} - \delta U^{\rm slope-lim}_{i+0}\right) \Rightarrow \boxed{U_{l,i+0} = \frac{1}{2}\left(U_{i+0} + U_{i-1}\right) + \frac{1}{6}\left(\delta U^{\rm slope-lim}_{i-1} - \delta U^{\rm slope-lim}_{i+0}\right)}\ . \end{align}After this step, the values of $U_{r,l,i+0}$ are stored as outputs.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// Finally, compute face values Ur and Ul based on the PPM prescription
// (Eq. A1 in http://arxiv.org/pdf/astro-ph/0503420.pdf, but using standard 1/6=(1.0/6.0) coefficient)
// Ur[PLUS0] represents U(i+1/2)
// We applied a simplification to the following line: Ur=U+0.5*(U(i+1)-U) + ... = 0.5*(U(i+1)+U) + ...
Ur[whichvar][PLUS0] = 0.5*(U[whichvar][PLUS1] + U[whichvar][PLUS0] ) + (1.0/6.0)*(slope_lim_dU[whichvar][PLUS0] - slope_lim_dU[whichvar][PLUS1]);
// Ul[PLUS0] represents U(i-1/2)
// We applied a simplification to the following line: Ul=U(i-1)+0.5*(U-U(i-1)) + ... = 0.5*(U+U(i-1)) + ...
Ul[whichvar][PLUS0] = 0.5*(U[whichvar][PLUS0] + U[whichvar][MINUS1]) + (1.0/6.0)*(slope_lim_dU[whichvar][MINUS1] - slope_lim_dU[whichvar][PLUS0]);
/* *** LOOP 1c: WRITE OUTPUT *** */
// Store right face values to {rho_br,Pr,vxr,vyr,vzr,Bxr,Byr,Bzr},
// and left face values to {rho_bl,Pl,vxl,vyl,vzl,Bxl,Byl,Bzl}
out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ur[whichvar][PLUS0];
out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0];
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
Following the procedure described in Step 4, we will now steepen $\rho_{b}$ using the steepen_rho()
function. Keep in mind that although we will loop over all primitives which are set for reconstruction, the steepening procedure is applied to $\rho_{b}$ only.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// *** LOOP 2: STEEPEN RHOB ***
// Note that this loop applies ONLY to RHOB.
if(whichvar==RHOB) {
LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
SET_INDEX_ARRAYS(-2,2,flux_dirn);
// Set rho and P separately, since within this loop,
// 1) steepen_rho() depends on RHOB(MINUS2,MINUS1,PLUS0,PLUS1,PLUS2)
// Read in all primitives between MINUS2 & PLUS2. Store to U.
for(int ii=MINUS2;ii<=PLUS2;ii++) U[RHOB][ii] = in_prims[RHOB ].gf[index_arr[flux_dirn][ii]];
for(int ii=MINUS1;ii<=PLUS1;ii++) U[PRESSURE][ii] = in_prims[PRESSURE].gf[index_arr[flux_dirn][ii]];
Ur[RHOB][PLUS0] = out_prims_r[RHOB].gf[index_arr[flux_dirn][PLUS0]];
Ul[RHOB][PLUS0] = out_prims_l[RHOB].gf[index_arr[flux_dirn][PLUS0]];
dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];
dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1];
dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0];
dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1];
slope_lim_dU[whichvar][MINUS1]=slope_limit(dU[whichvar][MINUS1],dU[whichvar][PLUS0]);
//slope_lim_dU[whichvar][PLUS0] =slope_limit(dU[whichvar][PLUS0], dU[whichvar][PLUS1]);
slope_lim_dU[whichvar][PLUS1] =slope_limit(dU[whichvar][PLUS1], dU[whichvar][PLUS2]);
// Steepen rho
// DEPENDENCIES: RHOB face values, RHOB(MINUS2,MINUS1,PLUS0,PLUS1,PLUS2), P(MINUS1,PLUS0,PLUS1), and slope_lim_dU[RHOB](MINUS1,PLUS1)
CCTK_REAL P_cold,Gamma_cold;
compute_P_cold__Gamma_cold(U[RHOB][PLUS0],eos, P_cold,Gamma_cold);
steepen_rho(U,slope_lim_dU, Gamma_th,P_cold,Gamma_cold, Ur[RHOB],Ul[RHOB]);
// Output rho
out_prims_r[RHOB].gf[index_arr[flux_dirn][PLUS0]] = Ur[RHOB][PLUS0];
out_prims_l[RHOB].gf[index_arr[flux_dirn][PLUS0]] = Ul[RHOB][PLUS0];
}
}
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
The flattening procedure modifies $U_{r}$ and $U_{l}$ via
$$ \boxed{ \begin{matrix} U_{r,i+0} = U_{i+0}\tilde{f} + U_{r,i+0}\left(1-\tilde{f}\right)\\ U_{l,i+0} = U_{i+0}\tilde{f} + U_{l,i+0}\left(1-\tilde{f}\right) \end{matrix} }\ , $$where $\tilde{f}$ is computed by the ftilde_compute()
function, described in Step 8.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
/* ORIGINAL PPM REQUIRES AT LEAST 4 GHOSTZONES, which can add
* significantly to the size of AMR ref. boundaries.
* To reduce to 3 ghostzones, we comment the following lines out:
* if ((P[indexp1] - P[indexm1]) <= 0.0) {
* f = MAX(ftilde,ftilde_p1);
* } else {
* f = MAX(ftilde,ftilde_m1);
* }
*/
// *** LOOP 3: FLATTEN BASED ON FTILDE AND MONOTONIZE ***
for(int ww=0;ww<num_prims_to_reconstruct;ww++) {
int whichvar=which_prims_to_reconstruct[ww];
// ftilde() depends on P(MINUS2,MINUS1,PLUS1,PLUS2)
LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
SET_INDEX_ARRAYS(0,0,flux_dirn);
U[whichvar][PLUS0] = in_prims[whichvar].gf[index_arr[flux_dirn][PLUS0]];
Ur[whichvar][PLUS0] = out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]];
Ul[whichvar][PLUS0] = out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]];
// ftilde_gf was computed in the function compute_ftilde_gf(), called before this routine
CCTK_REAL ftilde = ftilde_gf[index_arr[flux_dirn][PLUS0]];
// ...and then flatten (local operation)
Ur[whichvar][PLUS0] = U[whichvar][PLUS0]*ftilde + Ur[whichvar][PLUS0]*(1.0-ftilde);
Ul[whichvar][PLUS0] = U[whichvar][PLUS0]*ftilde + Ul[whichvar][PLUS0]*(1.0-ftilde);
Appending to ../src/reconstruct_set_of_prims_PPM.C
Furthermore, we apply the monotonize()
function, as described in Step 5.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// Then monotonize
monotonize(U[whichvar][PLUS0],Ur[whichvar][PLUS0],Ul[whichvar][PLUS0]);
out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ur[whichvar][PLUS0];
out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0];
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
At this point, we have
\begin{align} U_{r,i} &= U_{i+1/2}\ ,\\ U_{l,i} &= U_{i-1/2}\ . \end{align}To keep things consistent, we shift indices to get
$$ \boxed{ \begin{matrix} U_{i-1/2+\epsilon} = U_{l,i}^{\rm old} = U_{r,i}^{\rm new}\ ,\\ U_{i-1/2-\epsilon} = U_{r,i-1}^{\rm old} = U_{l,i}^{\rm new}\ ,\\ \end{matrix} } $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// Ur depends on ftilde, which depends on points of U between MINUS2 and PLUS2
out_prims_r[whichvar].gz_lo[flux_dirn]+=2;
out_prims_r[whichvar].gz_hi[flux_dirn]+=2;
// Ul depends on ftilde, which depends on points of U between MINUS2 and PLUS2
out_prims_l[whichvar].gz_lo[flux_dirn]+=2;
out_prims_l[whichvar].gz_hi[flux_dirn]+=2;
}
// *** LOOP 4: SHIFT Ur AND Ul ***
/* Currently face values are set so that
* a) Ur(i) represents U(i+1/2), and
* b) Ul(i) represents U(i-1/2)
* Here, we shift so that the indices are consistent:
* a) U(i-1/2+epsilon) = oldUl(i) = newUr(i)
* b) U(i-1/2-epsilon) = oldUr(i-1) = newUl(i)
* Note that this step is not strictly necessary if you keep
* track of indices when computing the flux. */
for(int ww=0;ww<num_prims_to_reconstruct;ww++) {
int whichvar=which_prims_to_reconstruct[ww];
LOOP_DEFINE(3,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
SET_INDEX_ARRAYS(-1,0,flux_dirn);
temporary[index_arr[flux_dirn][PLUS0]] = out_prims_r[whichvar].gf[index_arr[flux_dirn][MINUS1]];
}
LOOP_DEFINE(3,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
SET_INDEX_ARRAYS(0,0,flux_dirn);
// Then shift so that Ur represents the gridpoint at i-1/2+epsilon,
// and Ul represents the gridpoint at i-1/2-epsilon.
// Ur(i-1/2) = Ul(i-1/2) = U(i-1/2+epsilon)
// Ul(i-1/2) = Ur(i+1/2 - 1) = U(i-1/2-epsilon)
out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]];
out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = temporary[index_arr[flux_dirn][PLUS0]];
}
// Ul was just shifted, so we lost another ghostzone.
out_prims_l[whichvar].gz_lo[flux_dirn]+=1;
out_prims_l[whichvar].gz_hi[flux_dirn]+=0;
// As for Ur, we didn't need to get rid of another ghostzone,
// but we did ... seems wasteful!
out_prims_r[whichvar].gz_lo[flux_dirn]+=1;
out_prims_r[whichvar].gz_hi[flux_dirn]+=0;
}
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
slope_limit()
function [Back to top]¶We will now show the definition of $\delta U^{\rm slope-lim}$. The reason why we introduce this slope-limited procedure is twofold. First, it leads to steeper representations of discontinuities, and second it guarantees that $U_{i+1/2}$ lies inside the range $\left[U_{i},U_{i+1}\right]$.
We start by defining
\begin{align} \delta U_{i} &\equiv U_{i} - U_{i-1}\ ,\\ \delta_{m} U_{i} &\equiv \frac{\delta U_{i} + \delta U_{i+1}}{2} = \frac{U_{i+1} - U_{i-1}}{2}\ . \end{align}%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// Set SLOPE_LIMITER_COEFF = 2.0 for MC, 1 for minmod
#define SLOPE_LIMITER_COEFF 2.0
//Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996)
// [note the factor of 2 missing in the |a_{j+1} - a_{j}| term].
// Recall that dU = U_{i} - U_{i-1}.
static inline CCTK_REAL slope_limit(CCTK_REAL dU,CCTK_REAL dUp1) {
if(dU*dUp1 > 0.0) {
//delta_m_U=0.5 * [ (u_(i+1)-u_i) + (u_i-u_(i-1)) ] = (u_(i+1) - u_(i-1))/2 <-- first derivative, second-order; this should happen most of the time (smooth flows)
CCTK_REAL delta_m_U = 0.5*(dU + dUp1);
Appending to ../src/reconstruct_set_of_prims_PPM.C
Next we implement the slope-limited $\delta U$ as
$$ \boxed{\delta U^{\rm slope-lim} \equiv \left\{ \begin{matrix} {\rm sign}\left(\delta_{m} U_{i}\right)\min\left(\left|\delta_{m} U_{i}\right|,c\left|\delta U_{i}\right|,c\left|\delta U_{i+1}\right|\right) & ,\ {\rm if}\ dU_{i}dU_{i+1} > 0\\ 0 &,\ {\rm otherwise} \end{matrix} \right.}\ . $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// EXPLANATION OF BELOW LINE OF CODE.
// In short, sign_delta_a_j = sign(delta_m_U) = (0.0 < delta_m_U) - (delta_m_U < 0.0).
// If delta_m_U>0, then (0.0 < delta_m_U)==1, and (delta_m_U < 0.0)==0, so sign_delta_a_j=+1
// If delta_m_U<0, then (0.0 < delta_m_U)==0, and (delta_m_U < 0.0)==1, so sign_delta_a_j=-1
// If delta_m_U==0,then (0.0 < delta_m_U)==0, and (delta_m_U < 0.0)==0, so sign_delta_a_j=0
int sign_delta_m_U = (0.0 < delta_m_U) - (delta_m_U < 0.0);
//Decide whether to use 2nd order derivative or first-order derivative, limiting slope.
return sign_delta_m_U*MIN(fabs(delta_m_U),MIN(SLOPE_LIMITER_COEFF*fabs(dUp1),SLOPE_LIMITER_COEFF*fabs(dU)));
}
return 0.0;
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
steepen_rho()
function [Back to top]¶The steepening procedure within the PPM algorithm is applied only to $\rho_{b}$. The idea here is to produce narrower profiles near the vicinity of a contact discontinuity.
A NOTE ON NOTATION: in the discussion below we will refer to $\rho$ as $\rho$ to keep the notation a bit lighter. No confusion should arise from this since there is no other quantity $\rho$ involved.
We start the algorithm by computing
\begin{align} \delta\rho_{i+0} &= \frac{\rho_{i+1} - \rho_{i-1}}{2}\ ,\\ \delta^{2}\rho_{i-1} &= \rho_{i+0} - 2\rho_{i-1} + \rho_{i-2}\ ,\\ \delta^{2}\rho_{i+1} &= \rho_{i+2} - 2\rho_{i+1} + \rho_{i+0}\ . \end{align}%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// standard Colella-Woodward parameters:
// K0 = 0.1d0, eta1 = 20.0, eta2 = 0.05, epsilon = 0.01d0
#define K0 0.1
#define ETA1 20.0
#define ETA2 0.05
#define EPSILON 0.01
static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold,
CCTK_REAL *rho_br_ppm,CCTK_REAL *rho_bl_ppm) {
// Next compute centered differences d RHOB and d^2 RHOB
CCTK_REAL d1rho_b = 0.5*(U[RHOB][PLUS1] - U[RHOB][MINUS1]);
CCTK_REAL d2rho_b_m1 = U[RHOB][PLUS0] - 2.0*U[RHOB][MINUS1] + U[RHOB][MINUS2];
CCTK_REAL d2rho_b_p1 = U[RHOB][PLUS2] - 2.0*U[RHOB][PLUS1] + U[RHOB][PLUS0];
Appending to ../src/reconstruct_set_of_prims_PPM.C
Then we evaluate
$$ \Gamma = \left.\left(\frac{\partial P}{\partial\rho}\right)\middle/\left(\frac{P}{\rho}\right)\right. = \Gamma_{\rm th} + \left(\Gamma_{\rm cold} - \Gamma_{\rm th}\right)\frac{P_{\rm cold}}{P}\ . $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// Compute effective Gamma = (partial P / partial rho0)_s /(P/rho0)
CCTK_REAL Gamma = Gamma_th + (Gamma_cold-Gamma_th)*P_cold/U[PRESSURE][PLUS0];
Appending to ../src/reconstruct_set_of_prims_PPM.C
Next the contact discontinuity condition, eq. (3.2) of Colella & Woodward (1984), is checked:
$$ \Gamma K_{0}\frac{\left|\rho_{i+1}-\rho_{i-1}\right|}{\min\left(\rho_{i+1},\rho_{i-1}\right)} \geq \frac{\left|P_{i+1}-P_{i-1}\right|}{\min\left(P_{i+1},P_{i-1}\right)}\ , $$where $K_{0}$ is a problem dependent constant. Keep in mind that we implement the quantity
$$ \boxed{{\rm contact\_discontinuity\_check} \equiv \Gamma K_{0}\left|\rho_{i+1}-\rho_{i-1}\right|\min\left(P_{i+1},P_{i-1}\right) - \left|P_{i+1}-P_{i-1}\right|\min\left(\rho_{i+1},\rho_{i-1}\right)}\ , $$and verify whether ${\rm contact\_discontinuity\_check} \geq 0$ to verify the discontinuity condition. We also define the quantities
$$ \boxed{{\rm second\_deriv\_check} \equiv - \delta^{2}\rho_{i-1}\delta^{2}\rho_{i+1}}\ , $$and
$$ \boxed{{\rm relative\_change\_check} \equiv 2\left|\delta\rho\right| - \epsilon\min\left(\rho_{i+1},\rho_{i-1}\right)}\ , $$where again $\epsilon$ is a constant. The contact discontinuity condition is then satisfied when all three quantities inside the boxes above are non-negative. When that is the case, we evaluate
$$ \boxed{\eta_{i} = \max\left\{0,\min\left[\eta_{1}\left(\tilde\eta_{i}-\eta_{2}\right),1\right]\right\}}\ , $$where $\eta_{1}$ and $\eta_{2}$ are constants and
$$ \boxed{\tilde\eta_{i} = \left\{ \begin{matrix} 0&, \ {\rm if}\ \delta\rho_{i} = 0\ ,\\ -\frac{1}{6}\left(\frac{\delta^{2}\rho_{i+1} - \delta^{2}\rho_{i-1}}{2\delta\rho_{i}}\right)&, \ {\rm otherwise}\ . \end{matrix} \right.} $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
CCTK_REAL contact_discontinuity_check = Gamma*K0*fabs(U[RHOB][PLUS1]-U[RHOB][MINUS1])*
MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1])
-fabs(U[PRESSURE][PLUS1]-U[PRESSURE][MINUS1])*MIN(U[RHOB][PLUS1],U[RHOB][MINUS1]);
CCTK_REAL second_deriv_check = -d2rho_b_p1*d2rho_b_m1;
CCTK_REAL relative_change_check = fabs(2.0*d1rho_b) - EPSILON*MIN(U[RHOB][PLUS1],U[RHOB][MINUS1]);
if(contact_discontinuity_check >= 0.0 && second_deriv_check >= 0.0
&& relative_change_check >= 0.0) {
CCTK_REAL eta_tilde=0.0;
if (fabs(d1rho_b) > 0.0) {
eta_tilde = -(1.0/6.0)*(d2rho_b_p1-d2rho_b_m1)/(2.0*d1rho_b);
}
CCTK_REAL eta = MAX(0.0,MIN(ETA1*(eta_tilde - ETA2),1.0));
Appending to ../src/reconstruct_set_of_prims_PPM.C
We then apply the monotonized central (MC) scheme of van Leer (see also Step 3 for a discussion on the quantities $\delta\rho^{\rm slope-lim}_{i}$ below),
\begin{align} \rho^{\rm MC}_{r,i+1} &= \rho_{i+1} - \frac{1}{2}\delta\rho^{\rm slope-lim}_{i+1}\ ,\\ \rho^{\rm MC}_{l,i+0} &= \rho_{i-1} + \frac{1}{2}\delta\rho^{\rm slope-lim}_{i-1}\ , \end{align}so that, finally, the steepening algorithm sets
$$ \begin{matrix} \rho_{r}\to \rho_{r}(1-\eta) + \rho^{\rm MC}_{r}\eta\ ,\\ \rho_{l}\to \rho_{l}(1-\eta) + \rho^{\rm MC}_{l}\eta\ , \end{matrix} $$or, as implemented below:
$$ \boxed{\begin{matrix} \rho_{r,i+0}\rightarrow \rho_{r,i+0}\left(1-\eta_{i+0}\right) + \rho^{\rm MC}_{r,i+1}\eta_{i+0}\ ,\\ \rho_{l,i+0}\rightarrow \rho_{r,i+0}\left(1-\eta_{i+0}\right) + \rho^{\rm MC}_{l,i+0}\eta_{i+0}\ . \end{matrix}} $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// Next compute Urp1 and Ul for RHOB, using the MC prescription:
// Ur_p1 = U_p1 - 0.5*slope_lim_dU_p1
CCTK_REAL rho_br_mc_p1 = U[RHOB][PLUS1] - 0.5*slope_lim_dU[RHOB][PLUS1];
// Ul = U_m1 + 0.5*slope_lim_dU_m1
// Based on this line of code, Ur[index] = a_j - \delta_m a_j / 2. (cf. Eq. 65 in Marti & Muller's "PPM Method for 1D Relativistic Hydro." paper)
// So: Ur[indexp1] = a_{j+1} - \delta_m a_{j+1} / 2. This is why we have rho_br_mc[indexp1]
CCTK_REAL rho_bl_mc = U[RHOB][MINUS1] + 0.5*slope_lim_dU[RHOB][MINUS1];
rho_bl_ppm[PLUS0] = rho_bl_ppm[PLUS0]*(1.0-eta) + rho_bl_mc*eta;
rho_br_ppm[PLUS0] = rho_br_ppm[PLUS0]*(1.0-eta) + rho_br_mc_p1*eta;
}
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
monotonize()
function [Back to top]¶The value $U_{i+1/2}$ will be assigned to $U_{l,i}$ and $U_{r,i-1}$ for most values of $i$, but in some cases this would lead to incorrect interpolation results. Near discontinuities, the value of either $U_{l}$, $U_{r}$, or both need to be adjusted.
Consider, then, the following quantities:
\begin{align} \delta U &\equiv U_{r} - U_{l}\ ,\\ \delta_{m}U &\equiv \frac{U_{r} + U_{l}}{2}\ . \end{align}%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
static inline void monotonize(CCTK_REAL U,CCTK_REAL &Ur,CCTK_REAL &Ul) {
CCTK_REAL dU = Ur - Ul;
CCTK_REAL mU = 0.5*(Ur+Ul);
Appending to ../src/reconstruct_set_of_prims_PPM.C
Then, following Eq. (1.10) of Colella & Woodward (1984), we will check the following three cases:
$$ \boxed{\text{Case 1: if}\ \left(U_{r} - U\right)\left(U - U_{l}\right)\leq 0 \ \text{then:}\ \left\{ \begin{matrix} U_{r} \to U\ ,\\ \ U_{l}\to U\ . \end{matrix} \right.} $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
if ( (Ur-U)*(U-Ul) <= 0.0) {
Ur = U;
Ul = U;
return;
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
if ( dU*(U-mU) > (1.0/6.0)*SQR(dU)) {
Ul = 3.0*U - 2.0*Ur;
return;
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
if ( dU*(U-mU) < -(1.0/6.0)*SQR(dU)) {
Ur = 3.0*U - 2.0*Ul;
return;
}
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
compute_P_cold__Gamma_cold()
function [Back to top]¶This part of the code evaluates $P_{\rm cold}$ and $\Gamma_{\rm cold}$ for the equations of state (EOS) presented in Eqs. 13-16 of Stephens et al. (2008).
First, if $\rho_{b} = 0$, then $P_{\rm cold} = 0$ and $\Gamma_{\rm cold}$ simply receives its tabulated value.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold) {
// This code handles equations of state of the form defined
// in Eqs 13-16 in http://arxiv.org/pdf/0802.0200.pdf
// Default in case rho_b == 0.0
if(rho_b==0.0) { P_cold = 0.0; Gamma_cold = eos.Gamma_ppoly_tab[0]; return; }
Appending to ../src/reconstruct_set_of_prims_PPM.C
Next we consider the case where the EOS is given by a single-polytrope
$$ \boxed{P_{\rm cold} = \kappa \rho_{b}^{\Gamma_{\rm cold}}}\ , $$and also the piecewise polytrope EOS
$$ \boxed{ P_{\rm cold} = \left\{ \begin{matrix} K_{0}\rho^{\Gamma_{0}} & , & \rho \leq \rho_{0}\\ K_{1}\rho^{\Gamma_{1}} & , & \rho_{0} \leq \rho \leq \rho_{1}\\ \vdots & & \vdots\\ K_{j}\rho^{\Gamma_{j}} & , & \rho_{j-1} \leq \rho \leq \rho_{j}\\ \vdots & & \vdots\\ K_{N-1}\rho^{\Gamma_{N-1}} & , & \rho_{N-2} \leq \rho \leq \rho_{N-1}\\ K_{N}\rho^{\Gamma_{N}} & , & \rho \geq \rho_{N-1} \end{matrix} \right. }\ . $$Notice that we left the fact that $\Gamma_{i} \equiv \Gamma_{{\rm cold},i}$ implicit above.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
/***********************************
* Piecewise Polytropic EOS Patch *
* Computing P_cold and Gamma_cold *
***********************************/
int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b);
Gamma_cold = eos.Gamma_ppoly_tab[polytropic_index];
P_cold = eos.K_ppoly_tab[polytropic_index]*pow(rho_b,Gamma_cold);
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
#define OMEGA1 0.75
#define OMEGA2 10.0
#define EPSILON2 0.33
static void ftilde_gf_compute(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,gf_and_gz_struct *in_prims,CCTK_REAL *ftilde_gf) {
int ijkgz_lo_hi[4][2];
CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES];
/*Remove gcc unused variable warning/error Re: Pragma statement in loop define:*/
CCTK_REAL dU,slope_lim_dU,Ur,Ul; dU=slope_lim_dU=Ur=Ul=0.0; dU*=0;
// Compute ftilde, which is used for flattening left and right face values
LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[VX+(flux_dirn-1)].gz_lo,in_prims[VX+(flux_dirn-1)].gz_hi) {
SET_INDEX_ARRAYS(-2,2,flux_dirn);
for(int ii=MINUS2;ii<=PLUS2;ii++) U[PRESSURE][ii] = in_prims[PRESSURE].gf[index_arr[flux_dirn][ii]];
U[VX+(flux_dirn-1)][MINUS1] = in_prims[VX+(flux_dirn-1)].gf[index_arr[flux_dirn][MINUS1]];
U[VX+(flux_dirn-1)][PLUS1] = in_prims[VX+(flux_dirn-1)].gf[index_arr[flux_dirn][PLUS1]];
// Compute ftilde, which is used for flattening left and right face values
// DEPENDENCIES: P(MINUS2,MINUS1,PLUS1,PLUS2) and v^m(MINUS1,PLUS1), where m=flux_dirn={1,2,3}={x,y,z}.
ftilde_gf[index_arr[flux_dirn][PLUS0]] = ftilde_compute(flux_dirn,U);
}
}
Appending to ../src/reconstruct_set_of_prims_PPM.C
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
static inline CCTK_REAL ftilde_compute(const int flux_dirn,CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES]) {
CCTK_REAL dP1 = U[PRESSURE][PLUS1] - U[PRESSURE][MINUS1];
CCTK_REAL dP2 = U[PRESSURE][PLUS2] - U[PRESSURE][MINUS2];
Appending to ../src/reconstruct_set_of_prims_PPM.C
Then we modify the standard PPM algorithm slightly by introducing the following conditions:
\begin{align} {\rm if}\ \left|\frac{\delta P_{1}}{\delta_{m}P_{1}}\right| = 0\ {\rm then\ set}\ \delta P_{1}=0\ ,\\ {\rm if}\ \left|\frac{\delta P_{2}}{\delta_{m}P_{2}}\right| = 0\ {\rm then\ set}\ \delta P_{2}=0\ , \end{align}where
\begin{align} \delta_{m} P_{1} &\frac{\equiv P_{i+1} + P_{i-1}}{2}\ ,\\ \delta_{m} P_{2} &\frac{\equiv P_{i+2} + P_{i-2}}{2}\ . \end{align}Note that if the first condition above is satisfied then we are not inside a shock, while if the second condition is triggered alone there may be a shock.
%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// MODIFICATION TO STANDARD PPM:
// Cure roundoff error issues when dP1==0 or dP2==0 to 15 or more significant digits.
CCTK_REAL avg1=0.5*(U[PRESSURE][PLUS1] + U[PRESSURE][MINUS1]);
CCTK_REAL avg2=0.5*(U[PRESSURE][PLUS2] + U[PRESSURE][MINUS2]);
if(fabs(dP1)/avg1<1e-15) dP1=0.0; /* If this is triggered, there is NO shock */
if(fabs(dP2)/avg2<1e-15) dP2=0.0; /* If this is triggered alone, there may be a shock. Otherwise if triggered with above, NO shock. */
Appending to ../src/reconstruct_set_of_prims_PPM.C
Next we set
$$ {\rm dP1\_over\_dP2} = \left\{ \begin{matrix} \frac{\delta P_{1}}{\delta P_{2}} &,\ {\rm if}\ \delta P_{2} \neq 0\\ 1 &,\ {\rm otherwise} \end{matrix} \right. $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
CCTK_REAL dP1_over_dP2=1.0;
if (dP2 != 0.0) dP1_over_dP2 = dP1/dP2;
Appending to ../src/reconstruct_set_of_prims_PPM.C
We then construct
\begin{align} q_{1} &= \left({\rm dP1\_over\_dP2}-\omega_{1}\right)\omega_{2}\ ,\\ q_{2} &= \frac{\left|\delta P_{1}\right|}{\min\left(P_{i+1},P_{i-1}\right)}\ . \end{align}%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
CCTK_REAL q1 = (dP1_over_dP2-OMEGA1)*OMEGA2;
CCTK_REAL q2 = fabs(dP1)/MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1]);
Appending to ../src/reconstruct_set_of_prims_PPM.C
We then initialize a new variable, $w$, to
$$ w = \left\{ \begin{matrix} 1\ , &\ {\rm if}\ q_{2} > \epsilon_{2}\ {\rm and}\ q_{2}\left(v^{\rm flux\ dirn}_{i-1}-v^{\rm flux\ dirn}_{i+1}\right)>0\ \left({\rm inside\ shock}\right)\ ,\\ 0\ , &\ {\rm otherwise}\ \left({\rm outside\ shock}\right)\ , \end{matrix} \right. $$where $v^{\rm flux\ dirn}$ represents either $v^{x}$, $v^{y}$, or $v^{z}$, dependding on the flux direction. Finally, we get
$$ \boxed{\tilde{f} = \min\left[1,w\max\left(0,q_{1}\right)\right]}\ . $$%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C
// w==0 -> NOT inside a shock
CCTK_REAL w=0.0;
// w==1 -> inside a shock
if (q2 > EPSILON2 && q2*( (U[VX+(flux_dirn-1)][MINUS1]) - (U[VX+(flux_dirn-1)][PLUS1]) ) > 0.0) w = 1.0;
return MIN(1.0, w*MAX(0.0,q1));
}
Appending to ../src/reconstruct_set_of_prims_PPM.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/loop_defines_reconstruction.h"
original_IGM_file_name = "loop_defines_reconstruction-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__loop_defines_reconstruction__h = !diff $original_IGM_file_path $outfile_path__loop_defines_reconstruction__h
if Validation__loop_defines_reconstruction__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 loop_defines_reconstruction.h: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for loop_defines_reconstruction.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__loop_defines_reconstruction__h:
print(diff_line)
Validation test for loop_defines_reconstruction.h: FAILED! Diff: 17c17 < // This define only sets indices. --- > // This define only sets indices. 39a40 >
# 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/reconstruct_set_of_prims_PPM.C"
original_IGM_file_name = "reconstruct_set_of_prims_PPM-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__reconstruct_set_of_prims_PPM__C = !diff $original_IGM_file_path $outfile_path__reconstruct_set_of_prims_PPM__C
if Validation__reconstruct_set_of_prims_PPM__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 reconstruct_set_of_prims_PPM.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for reconstruct_set_of_prims_PPM.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__reconstruct_set_of_prims_PPM__C:
print(diff_line)
Validation test for reconstruct_set_of_prims_PPM.C: FAILED! Diff: 5c5 < * This version of PPM implements the standard --- > * This version of PPM implements the standard 7c7 < * to have 3 ghostzones instead of 4. --- > * to have 3 ghostzones instead of 4. 24c24 < CCTK_REAL gamma_th,CCTK_REAL P_cold,CCTK_REAL gamma_cold, --- > CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold, 26c26 < static inline void compute_P_cold__gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &gamma_cold); --- > static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold); 28a29 > 31a33,34 > DECLARE_CCTK_PARAMETERS; > 44c47,48 < --- > > 55c59 < // Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U. --- > // Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U. 57c61,62 < --- > > 59,60c64,65 < /* First, compute simple dU = U(i) - U(i-1), where direction of i < * is given by flux_dirn, and U is a primitive variable: --- > /* First, compute simple dU = U(i) - U(i-1), where direction of i > * is given by flux_dirn, and U is a primitive variable: 62c67 < // Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1), --- > // Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1), 64,67c69,72 < dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2]; < dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1]; < dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0]; < dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1]; --- > dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2]; > dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1]; > dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0]; > dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1]; 74c79,80 < // Finally, compute face values Ur and Ul based on the PPM prescription --- > > // Finally, compute face values Ur and Ul based on the PPM prescription 87c93 < out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0]; --- > out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0]; 89a96 > 104,107c111,114 < dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2]; < dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1]; < dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0]; < dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1]; --- > dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2]; > dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1]; > dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0]; > dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1]; 115,117c122,124 < CCTK_REAL P_cold=0,gamma_cold=0; < compute_P_cold__gamma_cold(U[RHOB][PLUS0],eos, P_cold,gamma_cold); < steepen_rho(U,slope_lim_dU, eos.gamma_th,P_cold,gamma_cold, Ur[RHOB],Ul[RHOB]); --- > CCTK_REAL P_cold,Gamma_cold; > compute_P_cold__Gamma_cold(U[RHOB][PLUS0],eos, P_cold,Gamma_cold); > steepen_rho(U,slope_lim_dU, Gamma_th,P_cold,Gamma_cold, Ur[RHOB],Ul[RHOB]); 125a133 > 142c150 < --- > 152a161 > 158a168 > 160c170 < out_prims_r[whichvar].gz_lo[flux_dirn]+=2; --- > out_prims_r[whichvar].gz_lo[flux_dirn]+=2; 185c195 < // Then shift so that Ur represents the gridpoint at i-1/2+epsilon, --- > // Then shift so that Ur represents the gridpoint at i-1/2+epsilon, 195c205 < // As for Ur, we didn't need to get rid of another ghostzone, --- > // As for Ur, we didn't need to get rid of another ghostzone, 202a213 > 206,207c217,218 < //Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996) < // [note the factor of 2 missing in the |a_{j+1} - a_{j}| term]. --- > //Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996) > // [note the factor of 2 missing in the |a_{j+1} - a_{j}| term]. 212a224 > 224a237 > 231c244 < static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL gamma_th,CCTK_REAL P_cold,CCTK_REAL gamma_cold, --- > static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold, 238a252 > 240,242c254,257 < CCTK_REAL Gamma = gamma_th + (gamma_cold-gamma_th)*P_cold/U[PRESSURE][PLUS0]; < CCTK_REAL contact_discontinuity_check = Gamma*K0*fabs(U[RHOB][PLUS1]-U[RHOB][MINUS1])* < MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1]) --- > CCTK_REAL Gamma = Gamma_th + (Gamma_cold-Gamma_th)*P_cold/U[PRESSURE][PLUS0]; > > CCTK_REAL contact_discontinuity_check = Gamma*K0*fabs(U[RHOB][PLUS1]-U[RHOB][MINUS1])* > MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1]) 247c262 < if(contact_discontinuity_check >= 0.0 && second_deriv_check >= 0.0 --- > if(contact_discontinuity_check >= 0.0 && second_deriv_check >= 0.0 254a270 > 263c279 < rho_bl_ppm[PLUS0] = rho_bl_ppm[PLUS0]*(1.0-eta) + rho_bl_mc*eta; --- > rho_bl_ppm[PLUS0] = rho_bl_ppm[PLUS0]*(1.0-eta) + rho_bl_mc*eta; 268a285 > 272,273c289,291 < < if ( (Ur-U)*(U-Ul) <= 0.0) { --- > > > if ( (Ur-U)*(U-Ul) <= 0.0) { 278c296,297 < if ( dU*(U-mU) > (1.0/6.0)*SQR(dU)) { --- > > if ( dU*(U-mU) > (1.0/6.0)*SQR(dU)) { 281a301 > 288c308,309 < static inline void compute_P_cold__gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &gamma_cold) { --- > > static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold) { 293c314,322 < if(rho_b==0.0) { P_cold = 0.0; gamma_cold = eos.gamma_tab[0]; return; } --- > if(rho_b==0.0) { P_cold = 0.0; Gamma_cold = eos.Gamma_ppoly_tab[0]; return; } > > /*********************************** > * Piecewise Polytropic EOS Patch * > * Computing P_cold and Gamma_cold * > ***********************************/ > int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b); > Gamma_cold = eos.Gamma_ppoly_tab[polytropic_index]; > P_cold = eos.K_ppoly_tab[polytropic_index]*pow(rho_b,Gamma_cold); 295,311d323 < if(eos.neos==1) { < // Eq. 14 of http://arxiv.org/pdf/0802.0200.pdf : < // P_{cold} = K_i rho_i^{\Gamma_i} < P_cold = eos.k_tab[0]*fasterpow_ppm_reconstruct(rho_b,eos.gamma_tab[0]); < gamma_cold = eos.gamma_tab[0]; < return; < } < for(int nn=1;nn<eos.neos;nn++) { < if (rho_b <= eos.rho_tab[nn] && rho_b > eos.rho_tab[nn-1]) { < P_cold = eos.k_tab[nn]*fasterpow_ppm_reconstruct(rho_b,eos.gamma_tab[nn]); < gamma_cold = eos.gamma_tab[nn]; < } < } < if (rho_b > eos.rho_tab[eos.neos-1]) { < P_cold = eos.k_tab[eos.neos]*fasterpow_ppm_reconstruct(rho_b,eos.gamma_tab[eos.neos]); < gamma_cold = eos.gamma_tab[eos.neos]; < } 313a326 > 334a348 > 338a353 > 345a361 > 348a365 > 351a369 > 359a378 >
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_reconstruct_set_of_prims_PPM.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__reconstruct_set_of_prims_PPM.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD_reconstruct__set_of_prims_PPM.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD_reconstruct__set_of_prims_PPM.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD_reconstruct__set_of_prims_PPM.tex
!rm -f Tut*.out Tut*.aux Tut*.log