#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # Tutorial-IllinoisGRMHD: outer_boundaries.C
# 
# ## Authors: Leo Werneck & Zach Etienne
# 
# <font color='red'>**This module is currently under development**</font>
# 
# ## This notebook explains the outer boundary conditions imposed on the quantities evolved within `IllinoisGRMHD`, specifically observing vector potential, hydrodynamic, and conservative variables. 
# 
# ### Required and recommended citations:
# 
# * **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).
# * **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).
# * **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618)).

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$
# 
# This module is organized as follows
# 
# 0. [Step 0](#src_dir): **Source directory creation**
# 1. [Step 1](#introduction): **Introduction**
# 1. [Step 2](#outer_boundaries__c): **`outer_boundaries.C`**
#     1. [Step 2.a](#outer_boundaries__amu): *The vector potential variables*
#         1. [Step 2.a.i](#outer_boundaries__amu__linear_extrapolation): Defining the linear extrapolation operators
#         1. [Step 2.a.ii](#outer_boundaries__amu__applying_bcs): Applying outer boundary conditions to $A_{\mu}$
#     1. [Step 2.b](#outer_boundaries__hydro_vars): *The hydrodynamic variables*
#         1. [Step 2.b.i](#outer_boundaries__hydro_vars__zero_deriv_outflow): Defining the zero derivative, outflow operators
#         1. [Step 2.b.ii](#outer_boundaries__hydro_vars__applying_bcs): Applying boundary conditions to $\left\{P,\rho_{b},v^{i}\right\}$
#     1. [Step 2.c](#outer_boundaries__conservatives): *The conservative variables*
# 1. [Step 3](#code_validation): **Code validation**
# 1. [Step 4](#latex_pdf_output): **Output this notebook to $\LaTeX$-formatted PDF file**

# <a id='src_dir'></a>
# 
# # Step 0: Source directory creation \[Back to [top](#toc)\]
# $$\label{src_dir}$$
# 
# We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory if it does not exist yet.

# In[1]:


# 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")


# <a id='introduction'></a>
# 
# # Step 1: Introduction \[Back to [top](#toc)\]
# $$\label{introduction}$$

# <a id='outer_boundaries__c'></a>
# 
# # Step 2: `outer_boundaries.C` \[Back to [top](#toc)\]
# $$\label{outer_boundaries__c}$$
# 
# 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](https://arxiv.org/pdf/1501.07276.pdf). 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\}$.

# <a id='outer_boundaries__amu'></a>
# 
# ## Step 2.a: The vector potential variables \[Back to [top](#toc)\]
# $$\label{outer_boundaries__amu}$$
# 
# <a id='outer_boundaries__amu__linear_extrapolation'></a>
# 
# ### Step 2.a.i: Defining the linear extrapolation operators \[Back to [top](#toc)\]
# $$\label{outer_boundaries__amu__linear_extrapolation}$$
# 
# 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.

# In[2]:


get_ipython().run_cell_magic('writefile', '$outfile_path__outer_boundaries__C', '/*******************************************************\n * Outer boundaries are handled as follows:\n * (-1) Update RHS quantities, leave RHS quantities zero on all outer ghostzones (including outer AMR refinement, processor, and outer boundaries)\n * ( 0) Let MoL update all evolution variables\n * ( 1) Apply outer boundary conditions (BCs) on A_{\\mu}\n * ( 2) Compute B^i from A_i everywhere, synchronize B^i\n * ( 3) Call con2prim to get primitives on interior pts\n * ( 4) Apply outer BCs on {P,rho_b,vx,vy,vz}.\n * ( 5) (optional) set conservatives on outer boundary.\n *******************************************************/\n\n#include "cctk.h"\n#include <cstdio>\n#include <cstdlib>\n#include <cmath>\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n\n#include "IllinoisGRMHD_headers.h"\n#include "IllinoisGRMHD_EoS_lowlevel_functs.C"\n#include "inlined_functions.C"\n\n#define IDX(i,j,k) CCTK_GFINDEX3D(cctkGH,(i),(j),(k))\n\n#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)];\n#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)];\n#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)];\n\n#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)];\n#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)];\n#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)];\n')


# <a id='outer_boundaries__amu__applying_bcs'></a>
# 
# ### Step 2.a.ii: Applying outer boundary conditions to $A_{\mu}$ \[Back to [top](#toc)\]
# $$\label{outer_boundaries__amu__applying_bcs}$$
# 
# 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](https://einsteintoolkit.org/usersguide/UsersGuidech9.html#x13-81000C1.6) for further details):
# 
# 1. `cctk_lsh[i]`: the number of *total* number of grid points along direction $x^{i}$, used *by each processor*.
# 2. `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:
#     1. `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}}$
#     1. `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}}$
#     1. `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}}$
#     1. `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}}$
#     1. `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}}$
#     1. `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}}$

# In[3]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__outer_boundaries__C', '\n\n/*********************************************\n * Apply outer boundary conditions on A_{\\mu}\n ********************************************/\nextern "C" void IllinoisGRMHD_outer_boundaries_on_A_mu(CCTK_ARGUMENTS) {\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n  if(CCTK_EQUALS(EM_BC,"frozen")) return;\n\n  bool Symmetry_none=false; if(CCTK_EQUALS(Symmetry,"none")) Symmetry_none=true;\n\n  int levelnumber = GetRefinementLevel(cctkGH);\n\n  IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh,  gxx,gxy,gxz,gyy,gyz,gzz,alp,\n                                                                gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,\n                                                                gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,\n                                                                phi_bssn,psi_bssn,lapm1);\n\n  // Don\'t apply approximate outer boundary conditions on initial data, which should be defined everywhere, or on levels != [coarsest level].\n  if(cctk_iteration==0 || levelnumber!=0) return;\n\n  if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2])\n    CCTK_VError(VERR_DEF_PARAMS,"ERROR: IllinoisGRMHD outer BC driver does not support unequal number of ghostzones in different directions!");\n  for(int which_bdry_pt=0;which_bdry_pt<cctk_nghostzones[0];which_bdry_pt++) {\n    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\n    int jmax=cctk_lsh[1]-cctk_nghostzones[1]+which_bdry_pt;\n    int kmax=cctk_lsh[2]-cctk_nghostzones[2]+which_bdry_pt;\n\n    int imin=cctk_nghostzones[0]-which_bdry_pt-1; // for cctk_nghostzones==3, this goes {2,1,0}\n    int jmin=cctk_nghostzones[1]-which_bdry_pt-1;\n    int kmin=cctk_nghostzones[2]-which_bdry_pt-1;\n\n    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); }\n    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); }\n    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); }\n\n    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); }\n    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); }\n    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); }\n  }\n}\n')


# <a id='outer_boundaries__hydro_vars'></a>
# 
# ## Step 2.b: The hydrodynamic variables \[Back to [top](#toc)\]
# $$\label{outer_boundaries__hydro_vars}$$
# 
# <a id='outer_boundaries__hydro_vars__zero_deriv_outflow'></a>
# 
# ### Step 2.b.i: Defining the zero derivative, outflow operators \[Back to [top](#toc)\]
# $$\label{outer_boundaries__hydro_vars__zero_deriv_outflow}$$
# 
# 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}
# }\ .
# $$

# In[4]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__outer_boundaries__C', '\n\n#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)];\n#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)];\n#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)];\n\n#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)];\n#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)];\n#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)];\n\n\n#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.;\n#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.;\n#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.;\n\n#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.;\n#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.;\n#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.;\n')


# <a id='outer_boundaries__hydro_vars__applying_bcs'></a>
# 
# ### Step 2.b.ii: Applying boundary conditions to $\left\{P,\rho_{b},v^{i}\right\}$ \[Back to [top](#toc)\]
# $$\label{outer_boundaries__hydro_vars__applying_bcs}$$
# 
# 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](#outer_boundaries__amu__applying_bcs), in case clarifications are needed.

# In[5]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__outer_boundaries__C', '\n\n\n/*******************************************************\n * Apply outer boundary conditions on {P,rho_b,vx,vy,vz}\n * It is better to apply BCs on primitives than conservs,\n * because small errors in conservs can be greatly\n * amplified in con2prim, sometimes leading to unphysical\n * primitives & unnecessary fixes.\n *******************************************************/\nextern "C" void IllinoisGRMHD_outer_boundaries_on_P_rho_b_vx_vy_vz(CCTK_ARGUMENTS) {\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n  if(CCTK_EQUALS(Matter_BC,"frozen")) return;\n\n  bool Symmetry_none=false; if(CCTK_EQUALS(Symmetry,"none")) Symmetry_none=true;\n\n  int levelnumber = GetRefinementLevel(cctkGH);\n\n  // Don\'t apply approximate outer boundary conditions on initial data, which should be defined everywhere, or on levels != [coarsest level].\n  if(cctk_iteration==0 || levelnumber!=0) return;\n\n  int ENABLE=1;\n\n  IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh,  gxx,gxy,gxz,gyy,gyz,gzz,alp,\n                                                                gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,\n                                                                gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,\n                                                                phi_bssn,psi_bssn,lapm1);\n\n  //if(levelnumber<=11110) {\n  if(cctk_nghostzones[0]!=cctk_nghostzones[1] || cctk_nghostzones[0]!=cctk_nghostzones[2])\n    CCTK_VError(VERR_DEF_PARAMS,"ERROR: IllinoisGRMHD outer BC driver does not support unequal number of ghostzones in different directions!");\n  for(int which_bdry_pt=0;which_bdry_pt<cctk_nghostzones[0];which_bdry_pt++) {\n    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\n    int jmax=cctk_lsh[1]-cctk_nghostzones[1]+which_bdry_pt;\n    int kmax=cctk_lsh[2]-cctk_nghostzones[2]+which_bdry_pt;\n\n    int imin=cctk_nghostzones[0]-which_bdry_pt-1; // for cctk_nghostzones==3, this goes {2,1,0}\n    int jmin=cctk_nghostzones[1]-which_bdry_pt-1;\n    int kmin=cctk_nghostzones[2]-which_bdry_pt-1;\n\n    // Order here is for compatibility with old version of this code.\n    /* XMIN & XMAX */\n    // i=imax=outer boundary\n    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); }\n    // i=imin=outer boundary\n    if(cctk_bbox[0]) {\n      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); }\n\n    /* YMIN & YMAX */\n    // j=jmax=outer boundary\n    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); }\n    // j=jmin=outer boundary\n    if(cctk_bbox[2]) {\n      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); }\n\n    /* ZMIN & ZMAX */\n    // k=kmax=outer boundary\n    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); }\n    // k=kmin=outer boundary\n    if((cctk_bbox[4]) && Symmetry_none) {\n      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); }\n  }\n\n')


# <a id='outer_boundaries__conservatives'></a>
# 
# ## Step 2.c: The conservative variables \[Back to [top](#toc)\]
# $$\label{outer_boundaries__conservatives}$$
# 
# 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](Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.ipynb). 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.

# In[6]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__outer_boundaries__C', '\n\n  /**********************************\n   * Piecewise Polytropic EOS Patch *\n   *   Setting up the EOS struct    *\n   **********************************/\n  /*\n   * The short piece of code below takes care\n   * of initializing the EOS parameters.\n   * Please refer to the "inlined_functions.C"\n   * source file for the documentation on the\n   * function.\n   */\n  eos_struct eos;\n  initialize_EOS_struct_from_input(eos);\n\n#pragma omp parallel for\n  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++) {\n        if(((cctk_bbox[0]) && i<cctk_nghostzones[0]) ||\n           ((cctk_bbox[1]) && i>=cctk_lsh[0]-cctk_nghostzones[0]) ||\n           ((cctk_bbox[2]) && j<cctk_nghostzones[1]) ||\n           ((cctk_bbox[3]) && j>=cctk_lsh[1]-cctk_nghostzones[1]) ||\n           ((cctk_bbox[4]) && k<cctk_nghostzones[2] && CCTK_EQUALS(Symmetry,"none")) ||\n           ((cctk_bbox[5]) && k>=cctk_lsh[2]-cctk_nghostzones[2])) {\n          int index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n          int ww;\n\n          CCTK_REAL METRIC[NUMVARS_FOR_METRIC],dummy=-1e100; // Set dummy to insane value, to ensure it isn\'t being used.\n          ww=0;\n          //psi[index] = exp(phi[index]);\n          METRIC[ww] = phi_bssn[index];ww++;\n          METRIC[ww] = dummy;          ww++; // Don\'t need to set psi.\n          METRIC[ww] = gtxx[index];    ww++;\n          METRIC[ww] = gtxy[index];    ww++;\n          METRIC[ww] = gtxz[index];    ww++;\n          METRIC[ww] = gtyy[index];    ww++;\n          METRIC[ww] = gtyz[index];    ww++;\n          METRIC[ww] = gtzz[index];    ww++;\n          METRIC[ww] = lapm1[index];   ww++;\n          METRIC[ww] = betax[index];   ww++;\n          METRIC[ww] = betay[index];   ww++;\n          METRIC[ww] = betaz[index];   ww++;\n          METRIC[ww] = gtupxx[index];  ww++;\n          METRIC[ww] = gtupyy[index];  ww++;\n          METRIC[ww] = gtupzz[index];  ww++;\n          METRIC[ww] = gtupxy[index];  ww++;\n          METRIC[ww] = gtupxz[index];  ww++;\n          METRIC[ww] = gtupyz[index];  ww++;\n\n          CCTK_REAL U[MAXNUMVARS];\n          ww=0;\n          U[ww] = rho_b[index]; ww++;\n          U[ww] = P[index];     ww++;\n          U[ww] = vx[index];    ww++;\n          U[ww] = vy[index];    ww++;\n          U[ww] = vz[index];    ww++;\n          U[ww] = Bx[index];    ww++;\n          U[ww] = By[index];    ww++;\n          U[ww] = Bz[index];    ww++;\n\n          struct output_stats stats;\n          CCTK_REAL CONSERVS[NUM_CONSERVS],TUPMUNU[10],TDNMUNU[10];\n\n          const int already_computed_physical_metric_and_inverse=0;\n          CCTK_REAL g4dn[4][4],g4up[4][4];\n          IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs(already_computed_physical_metric_and_inverse,U,stats,eos,METRIC,g4dn,g4up, TUPMUNU,TDNMUNU,CONSERVS);\n\n          rho_b[index] = U[RHOB];\n          P[index]     = U[PRESSURE];\n          vx[index]    = U[VX];\n          vy[index]    = U[VY];\n          vz[index]    = U[VZ];\n\n          rho_star[index]=CONSERVS[RHOSTAR];\n          tau[index]     =CONSERVS[TAUENERGY];\n          mhd_st_x[index]=CONSERVS[STILDEX];\n          mhd_st_y[index]=CONSERVS[STILDEY];\n          mhd_st_z[index]=CONSERVS[STILDEZ];\n\n          if(update_Tmunu) {\n            ww=0;\n            eTtt[index] = TDNMUNU[ww]; ww++;\n            eTtx[index] = TDNMUNU[ww]; ww++;\n            eTty[index] = TDNMUNU[ww]; ww++;\n            eTtz[index] = TDNMUNU[ww]; ww++;\n            eTxx[index] = TDNMUNU[ww]; ww++;\n            eTxy[index] = TDNMUNU[ww]; ww++;\n            eTxz[index] = TDNMUNU[ww]; ww++;\n            eTyy[index] = TDNMUNU[ww]; ww++;\n            eTyz[index] = TDNMUNU[ww]; ww++;\n            eTzz[index] = TDNMUNU[ww];\n          }\n          //if(i==5 && j==5 && k==5) CCTK_VInfo(CCTK_THORNSTRING,"%e %e %e %e",eTtt[index],eTtx[index],eTty[index],eTxy[index]);\n          //CCTK_VInfo(CCTK_THORNSTRING,"YAY: "); for(ww=0;ww<10;ww++) CCTK_VInfo(CCTK_THORNSTRING,"%e ",TDNMUNU[ww]); CCTK_VInfo(CCTK_THORNSTRING,"");\n        }\n      }\n}\n\n')


# <a id='code_validation'></a>
# 
# # Step 3: Code validation \[Back to [top](#toc)\]
# $$\label{code_validation}$$
# 
# First, we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook.

# In[7]:


# 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
        get_ipython().system('wget -O $original_IGM_file_path $original_IGM_file_url')

# Perform validation
Validation__outer_boundaries__C  = get_ipython().getoutput('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
    get_ipython().system('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)


# <a id='latex_pdf_output'></a>
# 
# # Step 4: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
# $$\label{latex_pdf_output}$$
# 
# 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](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).

# In[8]:


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
get_ipython().system('rm -f Tut*.out Tut*.aux Tut*.log')