This module is currently under development
IllinoisGRMHD
, specifically observing vector potential, hydrodynamic, and conservative variables.¶This module is organized as follows
outer_boundaries.C
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__outer_boundaries__C = os.path.join(IGM_src_dir_path,"outer_boundaries.C")
outer_boundaries.C
[Back to top]¶The strategy used to set the outer boundary for the primitives, $\left\{P,\rho_{b},v^{i}\right\}$, and for the scalar and vector potentials, $\left\{\left[\sqrt{\gamma}\Phi\right],A_{i}\right\}$, follows Eqs. (39) and (40) of the original release paper of IllinoisGRMHD. For example, if we are trying to apply boundary conditions along the $x$-direction, we would have
$$ \boxed{ E_{i+1} = \left\{ \begin{align} E_{i}\ , &{\rm\ if\ } E\in\left\{P,\rho_{b},v^{y},v^{z}\right\},{\rm\ or\ } E=v^{x}\ {\rm and\ } v^{x}\geq0\\ 0\ , &{\rm\ if\ } E=v^{x}\ {\rm and\ } v^{x}<0\\ 2E_{i} - E_{i-1}\ , &{\rm\ if\ } E\in\left\{\left[\sqrt{\gamma}\Phi\right],A_{x},A_{y},A_{z}\right\} \end{align} \right. }\ , $$for the ghostzone points along the positive $x$-direction, and
$$ \boxed{ E_{i-1} = \left\{ \begin{align} E_{i}\ , &{\rm\ if\ } E\in\left\{P,\rho_{b},v^{y},v^{z}\right\},{\rm\ or\ } E=v^{x}\ {\rm and\ } v^{x}\geq0\\ 0\ , &{\rm\ if\ } E=v^{x}\ {\rm and\ } v^{x}<0\\ 2E_{i} - E_{i+1}\ , &{\rm\ if\ } E\in\left\{\left[\sqrt{\gamma}\Phi\right],A_{x},A_{y},A_{z}\right\} \end{align} \right. }\ , $$for the ghostzone points along the negative $x$-direction.In this way, linear extrapolation outer boundary conditions are applied to the vector potential variables $\left\{\left[\sqrt{\gamma}\Phi\right],A_{i}\right\}$ and zero-derivative, outflow outer boundary conditions are applied to the hydrodynamic variables $\left\{P,\rho_{b},v^{i}\right\}$.
We start by applying outer boundary conditions to $\left\{\left[\sqrt{\gamma}\Phi\right],A_{i}\right\}$. We follow the prescription described above:
$$ \boxed{ \begin{align} \text{Positive direction: }E_{i+1} = 2E_{i} - E_{i-1}\ , &{\rm\ if\ } E\in\left\{\left[\sqrt{\gamma}\Phi\right],A_{x},A_{y},A_{z}\right\}\\ \text{Negative direction: }E_{i-1} = 2E_{i} - E_{i+1}\ , &{\rm\ if\ } E\in\left\{\left[\sqrt{\gamma}\Phi\right],A_{x},A_{y},A_{z}\right\} \end{align} }\ , $$which uses a linear extrapolation outer boundary condition.
%%writefile $outfile_path__outer_boundaries__C
/*******************************************************
* Outer boundaries are handled as follows:
* (-1) Update RHS quantities, leave RHS quantities zero on all outer ghostzones (including outer AMR refinement, processor, and outer boundaries)
* ( 0) Let MoL update all evolution variables
* ( 1) Apply outer boundary conditions (BCs) on A_{\mu}
* ( 2) Compute B^i from A_i everywhere, synchronize B^i
* ( 3) Call con2prim to get primitives on interior pts
* ( 4) Apply outer BCs on {P,rho_b,vx,vy,vz}.
* ( 5) (optional) set conservatives on outer boundary.
*******************************************************/
#include "cctk.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "IllinoisGRMHD_headers.h"
#include "IllinoisGRMHD_EoS_lowlevel_functs.C"
#include "inlined_functions.C"
#define IDX(i,j,k) CCTK_GFINDEX3D(cctkGH,(i),(j),(k))
#define XMAX_OB_LINEAR_EXTRAP(FUNC,imax) for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) FUNC[IDX(imax,j,k)] = 2.0 * FUNC[IDX(imax-1,j,k)] - FUNC[IDX(imax-2,j,k)];
#define YMAX_OB_LINEAR_EXTRAP(FUNC,jmax) for(int k=0;k<cctk_lsh[2];k++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,jmax,k)] = 2.0 * FUNC[IDX(i,jmax-1,k)] - FUNC[IDX(i,jmax-2,k)];
#define ZMAX_OB_LINEAR_EXTRAP(FUNC,kmax) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,j,kmax)] = 2.0 * FUNC[IDX(i,j,kmax-1)] - FUNC[IDX(i,j,kmax-2)];
#define XMIN_OB_LINEAR_EXTRAP(FUNC,imin) for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) FUNC[IDX(imin,j,k)] = 2.0 * FUNC[IDX(imin+1,j,k)] - FUNC[IDX(imin+2,j,k)];
#define YMIN_OB_LINEAR_EXTRAP(FUNC,jmin) for(int k=0;k<cctk_lsh[2];k++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,jmin,k)] = 2.0 * FUNC[IDX(i,jmin+1,k)] - FUNC[IDX(i,jmin+2,k)];
#define ZMIN_OB_LINEAR_EXTRAP(FUNC,kmin) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,j,kmin)] = 2.0 * FUNC[IDX(i,j,kmin+1)] - FUNC[IDX(i,j,kmin+2)];
Writing ../src/outer_boundaries.C
Now we apply boundary conditions to $A_{\mu}$. The code below is pretty straightforward, but it is useful to understand the following cctk
variables (refer to the Cactus Variables paragraph of section C1.6.2 of the Einstein Toolkit UserGuide for further details):
cctk_lsh[i]
: the number of total number of grid points along direction $x^{i}$, used by each processor.cctk_bbox[i]
: an array of integers that tell if the boundary gridpoints used by each processor are internal (i.e. artificial) or physical (i.e. actual boundary points). The variable follows the pattern:cctk_bbox[0]
: Direction: $x$ | Orientation: $+$ | Returns $\color{red}{0}$ if the boundary is $\color{red}{\text{artificial}}$ and $\color{blue}{1}$ if it is $\color{blue}{\text{physical}}$cctk_bbox[1]
: Direction: $x$ | Orientation: $-$ | Returns $\color{red}{0}$ if the boundary is $\color{red}{\text{artificial}}$ and $\color{blue}{1}$ if it is $\color{blue}{\text{physical}}$cctk_bbox[2]
: Direction: $y$ | Orientation: $+$ | Returns $\color{red}{0}$ if the boundary is $\color{red}{\text{artificial}}$ and $\color{blue}{1}$ if it is $\color{blue}{\text{physical}}$cctk_bbox[3]
: Direction: $y$ | Orientation: $-$ | Returns $\color{red}{0}$ if the boundary is $\color{red}{\text{artificial}}$ and $\color{blue}{1}$ if it is $\color{blue}{\text{physical}}$cctk_bbox[4]
: Direction: $z$ | Orientation: $+$ | Returns $\color{red}{0}$ if the boundary is $\color{red}{\text{artificial}}$ and $\color{blue}{1}$ if it is $\color{blue}{\text{physical}}$cctk_bbox[5]
: Direction: $z$ | Orientation: $-$ | Returns $\color{red}{0}$ if the boundary is $\color{red}{\text{artificial}}$ and $\color{blue}{1}$ if it is $\color{blue}{\text{physical}}$%%writefile -a $outfile_path__outer_boundaries__C
/*********************************************
* Apply outer boundary conditions on A_{\mu}
********************************************/
extern "C" void IllinoisGRMHD_outer_boundaries_on_A_mu(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
if(CCTK_EQUALS(EM_BC,"frozen")) return;
bool Symmetry_none=false; if(CCTK_EQUALS(Symmetry,"none")) Symmetry_none=true;
int levelnumber = GetRefinementLevel(cctkGH);
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);
// Don't apply approximate outer boundary conditions on initial data, which should be defined everywhere, or on levels != [coarsest level].
if(cctk_iteration==0 || levelnumber!=0) return;
if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2])
CCTK_VError(VERR_DEF_PARAMS,"ERROR: IllinoisGRMHD outer BC driver does not support unequal number of ghostzones in different directions!");
for(int which_bdry_pt=0;which_bdry_pt<cctk_nghostzones[0];which_bdry_pt++) {
int imax=cctk_lsh[0]-cctk_nghostzones[0]+which_bdry_pt; // for cctk_nghostzones==3, this goes {cctk_lsh-3,cctk_lsh-2,cctk_lsh-1}; outer bdry pt is at cctk_lsh-1
int jmax=cctk_lsh[1]-cctk_nghostzones[1]+which_bdry_pt;
int kmax=cctk_lsh[2]-cctk_nghostzones[2]+which_bdry_pt;
int imin=cctk_nghostzones[0]-which_bdry_pt-1; // for cctk_nghostzones==3, this goes {2,1,0}
int jmin=cctk_nghostzones[1]-which_bdry_pt-1;
int kmin=cctk_nghostzones[2]-which_bdry_pt-1;
if(cctk_bbox[1]) { XMAX_OB_LINEAR_EXTRAP(Ax,imax); XMAX_OB_LINEAR_EXTRAP(Ay,imax); XMAX_OB_LINEAR_EXTRAP(Az,imax); XMAX_OB_LINEAR_EXTRAP(psi6phi,imax); }
if(cctk_bbox[3]) { YMAX_OB_LINEAR_EXTRAP(Ax,jmax); YMAX_OB_LINEAR_EXTRAP(Ay,jmax); YMAX_OB_LINEAR_EXTRAP(Az,jmax); YMAX_OB_LINEAR_EXTRAP(psi6phi,jmax); }
if(cctk_bbox[5]) { ZMAX_OB_LINEAR_EXTRAP(Ax,kmax); ZMAX_OB_LINEAR_EXTRAP(Ay,kmax); ZMAX_OB_LINEAR_EXTRAP(Az,kmax); ZMAX_OB_LINEAR_EXTRAP(psi6phi,kmax); }
if(cctk_bbox[0]) { XMIN_OB_LINEAR_EXTRAP(Ax,imin); XMIN_OB_LINEAR_EXTRAP(Ay,imin); XMIN_OB_LINEAR_EXTRAP(Az,imin); XMIN_OB_LINEAR_EXTRAP(psi6phi,imin); }
if(cctk_bbox[2]) { YMIN_OB_LINEAR_EXTRAP(Ax,jmin); YMIN_OB_LINEAR_EXTRAP(Ay,jmin); YMIN_OB_LINEAR_EXTRAP(Az,jmin); YMIN_OB_LINEAR_EXTRAP(psi6phi,jmin); }
if((cctk_bbox[4]) && Symmetry_none) { ZMIN_OB_LINEAR_EXTRAP(Ax,kmin); ZMIN_OB_LINEAR_EXTRAP(Ay,kmin); ZMIN_OB_LINEAR_EXTRAP(Az,kmin); ZMIN_OB_LINEAR_EXTRAP(psi6phi,kmin); }
}
}
Appending to ../src/outer_boundaries.C
We now apply outer boundary conditions to $\left\{P,\rho_{b},v^{i}\right\}$, imposing zero derivative, outflow boundary conditions. We follow the prescription described above:
$$ \boxed{ \begin{matrix} \text{Positive direction: }E_{i+1} = \left\{ \begin{matrix} E_{i}\ , &{\rm\ if\ } E\in\left\{P,\rho_{b},v^{y},v^{z}\right\},{\rm\ or\ } E=v^{x}\ {\rm and\ } v^{x}\geq0\\ 0\ , &{\rm\ if\ } E=v^{x}\ {\rm and\ } v^{x}<0 \end{matrix} \right.\\ \text{Negative direction: }E_{i-1} = \left\{ \begin{matrix} E_{i}\ , &{\rm\ if\ } E\in\left\{P,\rho_{b},v^{y},v^{z}\right\},{\rm\ or\ } E=v^{x}\ {\rm and\ } v^{x}\geq0\\ 0\ , &{\rm\ if\ } E=v^{x}\ {\rm and\ } v^{x}<0 \end{matrix} \right. \end{matrix} }\ . $$%%writefile -a $outfile_path__outer_boundaries__C
#define XMAX_OB_SIMPLE_COPY(FUNC,imax) for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) FUNC[IDX(imax,j,k)] = FUNC[IDX(imax-1,j,k)];
#define YMAX_OB_SIMPLE_COPY(FUNC,jmax) for(int k=0;k<cctk_lsh[2];k++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,jmax,k)] = FUNC[IDX(i,jmax-1,k)];
#define ZMAX_OB_SIMPLE_COPY(FUNC,kmax) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,j,kmax)] = FUNC[IDX(i,j,kmax-1)];
#define XMIN_OB_SIMPLE_COPY(FUNC,imin) for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) FUNC[IDX(imin,j,k)] = FUNC[IDX(imin+1,j,k)];
#define YMIN_OB_SIMPLE_COPY(FUNC,jmin) for(int k=0;k<cctk_lsh[2];k++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,jmin,k)] = FUNC[IDX(i,jmin+1,k)];
#define ZMIN_OB_SIMPLE_COPY(FUNC,kmin) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) FUNC[IDX(i,j,kmin)] = FUNC[IDX(i,j,kmin+1)];
#define XMAX_INFLOW_CHECK(vx,imax) for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) if(vx[IDX(imax,j,k)]<0.) vx[IDX(imax,j,k)]=0.;
#define YMAX_INFLOW_CHECK(vy,jmax) for(int k=0;k<cctk_lsh[2];k++) for(int i=0;i<cctk_lsh[0];i++) if(vy[IDX(i,jmax,k)]<0.) vy[IDX(i,jmax,k)]=0.;
#define ZMAX_INFLOW_CHECK(vz,kmax) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) if(vz[IDX(i,j,kmax)]<0.) vz[IDX(i,j,kmax)]=0.;
#define XMIN_INFLOW_CHECK(vx,imin) for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) if(vx[IDX(imin,j,k)]>0.) vx[IDX(imin,j,k)]=0.;
#define YMIN_INFLOW_CHECK(vy,jmin) for(int k=0;k<cctk_lsh[2];k++) for(int i=0;i<cctk_lsh[0];i++) if(vy[IDX(i,jmin,k)]>0.) vy[IDX(i,jmin,k)]=0.;
#define ZMIN_INFLOW_CHECK(vz,kmin) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) if(vz[IDX(i,j,kmin)]>0.) vz[IDX(i,j,kmin)]=0.;
Appending to ../src/outer_boundaries.C
As with the previous case, applying the boundary conditions is a straightforward procedure. We refer the reader to the cctk
quantities discussed in Step 2.a.ii, in case clarifications are needed.
%%writefile -a $outfile_path__outer_boundaries__C
/*******************************************************
* Apply outer boundary conditions on {P,rho_b,vx,vy,vz}
* It is better to apply BCs on primitives than conservs,
* because small errors in conservs can be greatly
* amplified in con2prim, sometimes leading to unphysical
* primitives & unnecessary fixes.
*******************************************************/
extern "C" void IllinoisGRMHD_outer_boundaries_on_P_rho_b_vx_vy_vz(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
if(CCTK_EQUALS(Matter_BC,"frozen")) return;
bool Symmetry_none=false; if(CCTK_EQUALS(Symmetry,"none")) Symmetry_none=true;
int levelnumber = GetRefinementLevel(cctkGH);
// Don't apply approximate outer boundary conditions on initial data, which should be defined everywhere, or on levels != [coarsest level].
if(cctk_iteration==0 || levelnumber!=0) return;
int ENABLE=1;
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);
//if(levelnumber<=11110) {
if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2])
CCTK_VError(VERR_DEF_PARAMS,"ERROR: IllinoisGRMHD outer BC driver does not support unequal number of ghostzones in different directions!");
for(int which_bdry_pt=0;which_bdry_pt<cctk_nghostzones[0];which_bdry_pt++) {
int imax=cctk_lsh[0]-cctk_nghostzones[0]+which_bdry_pt; // for cctk_nghostzones==3, this goes {cctk_lsh-3,cctk_lsh-2,cctk_lsh-1}; outer bdry pt is at cctk_lsh-1
int jmax=cctk_lsh[1]-cctk_nghostzones[1]+which_bdry_pt;
int kmax=cctk_lsh[2]-cctk_nghostzones[2]+which_bdry_pt;
int imin=cctk_nghostzones[0]-which_bdry_pt-1; // for cctk_nghostzones==3, this goes {2,1,0}
int jmin=cctk_nghostzones[1]-which_bdry_pt-1;
int kmin=cctk_nghostzones[2]-which_bdry_pt-1;
// Order here is for compatibility with old version of this code.
/* XMIN & XMAX */
// i=imax=outer boundary
if(cctk_bbox[1]) { XMAX_OB_SIMPLE_COPY(P,imax); XMAX_OB_SIMPLE_COPY(rho_b,imax); XMAX_OB_SIMPLE_COPY(vx,imax); XMAX_OB_SIMPLE_COPY(vy,imax); XMAX_OB_SIMPLE_COPY(vz,imax); if(ENABLE) XMAX_INFLOW_CHECK(vx,imax); }
// i=imin=outer boundary
if(cctk_bbox[0]) {
XMIN_OB_SIMPLE_COPY(P,imin); XMIN_OB_SIMPLE_COPY(rho_b,imin); XMIN_OB_SIMPLE_COPY(vx,imin); XMIN_OB_SIMPLE_COPY(vy,imin); XMIN_OB_SIMPLE_COPY(vz,imin); if(ENABLE) XMIN_INFLOW_CHECK(vx,imin); }
/* YMIN & YMAX */
// j=jmax=outer boundary
if(cctk_bbox[3]) { YMAX_OB_SIMPLE_COPY(P,jmax); YMAX_OB_SIMPLE_COPY(rho_b,jmax); YMAX_OB_SIMPLE_COPY(vx,jmax); YMAX_OB_SIMPLE_COPY(vy,jmax); YMAX_OB_SIMPLE_COPY(vz,jmax); if(ENABLE) YMAX_INFLOW_CHECK(vy,jmax); }
// j=jmin=outer boundary
if(cctk_bbox[2]) {
YMIN_OB_SIMPLE_COPY(P,jmin); YMIN_OB_SIMPLE_COPY(rho_b,jmin); YMIN_OB_SIMPLE_COPY(vx,jmin); YMIN_OB_SIMPLE_COPY(vy,jmin); YMIN_OB_SIMPLE_COPY(vz,jmin); if(ENABLE) YMIN_INFLOW_CHECK(vy,jmin); }
/* ZMIN & ZMAX */
// k=kmax=outer boundary
if(cctk_bbox[5]) { ZMAX_OB_SIMPLE_COPY(P,kmax); ZMAX_OB_SIMPLE_COPY(rho_b,kmax); ZMAX_OB_SIMPLE_COPY(vx,kmax); ZMAX_OB_SIMPLE_COPY(vy,kmax); ZMAX_OB_SIMPLE_COPY(vz,kmax); if(ENABLE) ZMAX_INFLOW_CHECK(vz,kmax); }
// k=kmin=outer boundary
if((cctk_bbox[4]) && Symmetry_none) {
ZMIN_OB_SIMPLE_COPY(P,kmin); ZMIN_OB_SIMPLE_COPY(rho_b,kmin); ZMIN_OB_SIMPLE_COPY(vx,kmin); ZMIN_OB_SIMPLE_COPY(vy,kmin); ZMIN_OB_SIMPLE_COPY(vz,kmin); if(ENABLE) ZMIN_INFLOW_CHECK(vz,kmin); }
}
Appending to ../src/outer_boundaries.C
After we have applied boundary conditions to our primitives (i.e. hydrodynamics) variables, we make sure their values lie within the physical range and then recompute the conservatives. Notice that the boundary conditions are then not applied directly to the conservative variables. The reason why the code is structured in this way is because small variations in the values of the conservative variables can cause the conservative-to-primitive algorithm to fail.
%%writefile -a $outfile_path__outer_boundaries__C
/**********************************
* 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);
#pragma omp parallel for
for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
if(((cctk_bbox[0]) && i<cctk_nghostzones[0]) ||
((cctk_bbox[1]) && i>=cctk_lsh[0]-cctk_nghostzones[0]) ||
((cctk_bbox[2]) && j<cctk_nghostzones[1]) ||
((cctk_bbox[3]) && j>=cctk_lsh[1]-cctk_nghostzones[1]) ||
((cctk_bbox[4]) && k<cctk_nghostzones[2] && CCTK_EQUALS(Symmetry,"none")) ||
((cctk_bbox[5]) && k>=cctk_lsh[2]-cctk_nghostzones[2])) {
int index = CCTK_GFINDEX3D(cctkGH,i,j,k);
int ww;
CCTK_REAL METRIC[NUMVARS_FOR_METRIC],dummy=-1e100; // Set dummy to insane value, to ensure it isn't being used.
ww=0;
//psi[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++;
CCTK_REAL U[MAXNUMVARS];
ww=0;
U[ww] = rho_b[index]; ww++;
U[ww] = P[index]; ww++;
U[ww] = vx[index]; ww++;
U[ww] = vy[index]; ww++;
U[ww] = vz[index]; ww++;
U[ww] = Bx[index]; ww++;
U[ww] = By[index]; ww++;
U[ww] = Bz[index]; ww++;
struct output_stats stats;
CCTK_REAL CONSERVS[NUM_CONSERVS],TUPMUNU[10],TDNMUNU[10];
const int already_computed_physical_metric_and_inverse=0;
CCTK_REAL g4dn[4][4],g4up[4][4];
IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs(already_computed_physical_metric_and_inverse,U,stats,eos,METRIC,g4dn,g4up, TUPMUNU,TDNMUNU,CONSERVS);
rho_b[index] = U[RHOB];
P[index] = U[PRESSURE];
vx[index] = U[VX];
vy[index] = U[VY];
vz[index] = U[VZ];
rho_star[index]=CONSERVS[RHOSTAR];
tau[index] =CONSERVS[TAUENERGY];
mhd_st_x[index]=CONSERVS[STILDEX];
mhd_st_y[index]=CONSERVS[STILDEY];
mhd_st_z[index]=CONSERVS[STILDEZ];
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];
}
//if(i==5 && j==5 && k==5) CCTK_VInfo(CCTK_THORNSTRING,"%e %e %e %e",eTtt[index],eTtx[index],eTty[index],eTxy[index]);
//CCTK_VInfo(CCTK_THORNSTRING,"YAY: "); for(ww=0;ww<10;ww++) CCTK_VInfo(CCTK_THORNSTRING,"%e ",TDNMUNU[ww]); CCTK_VInfo(CCTK_THORNSTRING,"");
}
}
}
Appending to ../src/outer_boundaries.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/outer_boundaries.C"
original_IGM_file_name = "outer_boundaries-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__outer_boundaries__C = !diff $original_IGM_file_path $outfile_path__outer_boundaries__C
if Validation__outer_boundaries__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 outer_boundaries.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for outer_boundaries.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__outer_boundaries__C:
print(diff_line)
Validation test for outer_boundaries.C: FAILED! Diff: 19a20 > #include "IllinoisGRMHD_EoS_lowlevel_functs.C" 31a33 > 53c55 < if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2]) --- > if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2]) 67c69 < --- > 73a76 > 91a95 > 95c99 < * because small errors in conservs can be greatly --- > * because small errors in conservs can be greatly 120c124 < if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2]) --- > if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2]) 148c152 < if(cctk_bbox[5]) { ZMAX_OB_SIMPLE_COPY(P,kmax); ZMAX_OB_SIMPLE_COPY(rho_b,kmax); ZMAX_OB_SIMPLE_COPY(vx,kmax); ZMAX_OB_SIMPLE_COPY(vy,kmax); ZMAX_OB_SIMPLE_COPY(vz,kmax); if(ENABLE) ZMAX_INFLOW_CHECK(vz,kmax); } --- > if(cctk_bbox[5]) { ZMAX_OB_SIMPLE_COPY(P,kmax); ZMAX_OB_SIMPLE_COPY(rho_b,kmax); ZMAX_OB_SIMPLE_COPY(vx,kmax); ZMAX_OB_SIMPLE_COPY(vy,kmax); ZMAX_OB_SIMPLE_COPY(vz,kmax); if(ENABLE) ZMAX_INFLOW_CHECK(vz,kmax); } 155c159,170 < // 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. > */ 157,164c172 < 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); 212c220 < --- > 246a255 >
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__outer_boundaries.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__outer_boundaries.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__outer_boundaries.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__outer_boundaries.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__outer_boundaries.tex
!rm -f Tut*.out Tut*.aux Tut*.log