


# Start-to-Finish Example: Applying Boundary Conditions in Curvilinear Coordinates for Rank-0, Rank-1, and Symmetric Rank-2 Tensors in Three Dimensions

## Author: Zach Etienne
### Formatting improvements courtesy Brandon Clark


## This module documents and validates basic boundary condition algorithms for curvilinear coordinate systems (e.g., Spherical, Cylindrical), based on prescription in the [SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658).

## Introduction:
Following the prescription in the [SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658), we will implement curvilinear boundary conditions for rank-0, rank-1, and symmetric rank-2 tensors in three dimensions; as this is the same dimension and highest rank needed for BSSN.

### Introduction to parity conditions:

Suppose we have a vector $v^\rho$ defined at ghostzone $(-\rho,\phi,z)$ ($\rho>0$) in cylindrical coordinates. This will map to an interior point at $(\rho,\phi+\pi,z)$. At this point, the direction of the $\hat{\rho}$ unit vector flips sign. Thus we cannot simply set the value of $v^\rho$ to the value it possesses at interior point $(\rho,\phi+\pi,z)$; that would result in a sign error. Instead we have
\begin{align}
v^\rho(-\rho,\phi,z)&=-v^\rho(\rho,\phi+\pi,z) \\
&= \mathbf{e}^\rho\left(-\rho,\phi,z\right) \cdot \mathbf{e}^\rho\left(\rho,\phi+\pi,z\right)v^\rho(\rho,\phi+\pi,z),
\end{align}
where $\mathbf{e}^\rho\left(\rho,\phi,z\right)$ is the $\rho$ unit vector evaluated at point $(\rho,\phi,z)$, and $\mathbf{e}^\rho\left(-\rho,\phi,z\right) \cdot \mathbf{e}^\rho\left(\rho,\phi+\pi,z\right)$ is the dot product of the two unit vectors, which must evaluate to $\pm 1$ (i.e., the **parity**). Contrast this with scalars, which do not possess a sense of direction/parity.

### Basic algorithm for applying boundary conditions in curvilinear coordinates:

At each ghost zone grid point $\mathbf{d}_{\rm gz}=(x_0,x_1,x_2)$, we will do the following:

1. Evaluate the Cartesian coordinate $\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\right)$, corresponding to this grid point. Then evaluate the inverse $\mathbf{d}_{\rm new}=\left(x_0(x,y,z),x_1(x,y,z),x_2(x,y,z)\right)$. 
    1. If $\mathbf{d}_{\rm new} \ne \mathbf{d}_{\rm gz}$, then the ghost zone grid point maps to a point in the grid interior, *which is exactly the case described in the above section*. To distinguish this case from an "outer boundary condition", we shall henceforth refer to it variously as an application of an "interior", "inner", or "parity" boundary condition.
    1. If $\mathbf{d}_{\rm new} \equiv \mathbf{d}_{\rm gz}$, then the ghost zone grid point is on the outer boundary of the grid, and standard outer boundary conditions should be applied.

### Applying parity conditions to arbitrary-rank tensors:

Above we presented the strategy for applying parity boundary conditions to a single component of a vector. Here we outline the generic algorithm for arbitrary-rank tensors.

Continuing the discussion from the previous section, we assume $\mathbf{d}_{\rm new} \ne \mathbf{d}_{\rm gz}$ (otherwise we would apply the *outer* boundary condition algorithm). Next suppose we are given a generic rank-$N$ tensor ($N>0$).

1. The first component of the rank-$N$ tensor corresponds to some direction with unit vector $\mathbf{e}^i$; e.g., $v^r$ corresponds to the $\mathbf{e}^r$ direction. Compute the dot product of the unit vector $\mathbf{e}^i$ evaluated at points $\mathbf{d}_{\rm gz}$ and $\mathbf{d}_{\rm new}$. Define this dot product as $P_1$ ("$P$" for "parity"):
$$
P_1 = \mathbf{e}^i\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^i\left(\mathbf{d}_{\rm new}\right).
$$
1. $P_1$ will take the value of $\pm 1$, depending on the unit-vector direction and the points $\mathbf{d}_{\rm gz}$ and $\mathbf{d}_{\rm new}$
1. Repeat the above for the remaining components of the rank-$N$ tensor $j\in \{2,3,...,N\}$, storing each $P_j$.
1. The tensor mapping from $\mathbf{d}_{\rm gz}$ to $\mathbf{d}_{\rm new}$ for this tensor $T^{ijk...}_{mnp...}$ will be given by
$$
T^{ijk...}_{lmn...}(x_0,x_1,x_2)_{\rm gz} = \prod_{\ell=1}^N P_\ell T^{ijk...}_{mnp...}(x_0,x_1,x_2)_{\rm new}.
$$

In this formulation of BSSN, we only need to deal with rank-0, rank-1, and *symmetric* rank-2 tensors. Further, our basis consists of 3 directions, so there are a total of 
+ 1 parity condition (the trivial +1) for scalars (rank-0 tensors)
+ 3 parity conditions for all rank-1 tensors (corresponding to each direction)
+ 6 parity conditions for all *symmetric* rank-2 tensors (corresponding to the number of elements in the lower or upper triangle of a $3\times3$ matrix, including the diagonal)

Thus we must keep track of the behavior of **10 separate parity conditions**, which can be evaluated once the numerical grid has been set up, for all time.

<a id='toc'></a>

# Table of Contents
$$\label{toc}$$

This notebook is organized as follows

1. [Step 1](#initializenrpy): Set core NRPy+ parameters for numerical grids and reference metric
1. [Step 2](#ccode_bcs): Implement C code for implementation of curvilinear coordinate boundary conditions
    1. [Step 2.a](#register_gfs): Register gridfunctions of all 10 parity types; output gridfunction aliases to `CurviBoundaryConditions/gridfunction\_defines.h`
    1. [Step 2.b](#assign_parity): Assign the correct basic parity type to each test gridfunction
        1. [Step 2.b.i](#parity_digits): Set parity type for each gridfunction, based on the digits at the end of its name
        1. [Step 2.b.ii](#dot): Set up unit-vector dot products (=parity) for each of the 10 parity condition types
    1. [Step 2.c](#gzmap): Implement a modified version of the scalar wave in curvilinear coordinates boundary condition ghost zone mapping routine
    1. [Step 2.d](#bpcs): Construct the C loops that apply boundary and parity conditions
1. [Step 3](#validate): Set up test data for Curvilinear Boundary Conditions code validation
1. [Step 4](#mainc): `CurviBC_Playground.c` : The Main C Code
1. [Step 5](#senr_compare): Validation: Compare with original SENR results
1. [Step 6](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF file
        
    

<a id='initializenrpy'></a>

# Step 1: Set core NRPy+ parameters for numerical grids and reference metric \[Back to [top](#toc)\]
$$\label{initializenrpy}$$

We import needed NRPy+ modules and set reference_metric::CoordSystem=Spherical.


```python
# First we import needed core NRPy+ modules
from outputC import *
import NRPy_param_funcs as par
import grid as gri
import loop as lp
import indexedexp as ixp
import finite_difference as fin
import reference_metric as rfm

# Set spatial dimension (must be 3 for BSSN)
DIM = 3
par.set_parval_from_str("grid::DIM",DIM)

# Set the finite differencing order to 4; although this module doesn't compute 
#   finite difference derivatives, it does define NGHOSTS, which depends
#   on finite_difference::FD_CENTDERIVS_ORDER being set.
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER",4)

# Then we set the coordinate system for the numerical grid
par.set_parval_from_str("reference_metric::CoordSystem","Spherical")
rfm.reference_metric()
```

<a id='ccode_bcs'></a>

# Step 2: Implement C code for implementation of curvilinear coordinate boundary conditions \[Back to [top](#toc)\]
$$\label{ccode_bcs}$$


<a id='register_gfs'></a>

## Step 2.a: Register gridfunctions of all 10 parity types; output gridfunction aliases to `CurviBoundaryConditions/gridfunction_defines.h` \[Back to [top](#toc)\]
$$\label{register_gfs}$$ 

Here we 

1. Register within NRPy+ one gridfunction per each of the 10 parity conditions, and then 
1. output to file the corresponding gridfunction aliases, so the C code can access each gridfunction by its human-friendly alias (e.g., `test_gfs[RANKONEU0GF][idx]` instead of `test_gfs[6][idx]`).


```python
# Step 2.A.1: Register gridfunctions of all 10 parity types

# 6 gridfunctions, corresponding to all unique rank-2 tensor components:
ranktwosymmDD = ixp.register_gridfunctions_for_single_rank2("AUX","ranktwosymmDD", "sym01")
# 3 gridfunctions, corresponding to all unique rank-1 tensor components:
rankoneU = ixp.register_gridfunctions_for_single_rank1("AUX","rankoneU")
# 1 rank-0 (scalar) gridfunction
rankzero = ixp.gri.register_gridfunctions("AUX","rankzero")

# Step 2.A.2. Output gridfunction #define aliases to file:
!mkdir CurviBoundaryConditions 2>/dev/null # 2>/dev/null: Don't throw an error if the directory already exists.
evolved_variables_list, auxiliary_variables_list, auxevol_variables_list = \
        gri.output__gridfunction_defines_h__return_gf_lists("CurviBoundaryConditions")
```

<a id='assign_parity'></a>

## Step 2.b: Assign the correct basic parity type to each test gridfunction \[Back to [top](#toc)\]
$$\label{assign_parity}$$

The appropriate symbolic dot products determining parity condition are assigned to each gridfunction based on the following numbering:

Tensor type | Parity type (Step 2.B.1) | Dot product(s) determining parity condition (Step 2.B.2)
--- | --- | ---
Scalar (Rank-0 tensor) | 0 | (*none*)
Rank-1 tensor in **i0** direction | 1 | $\mathbf{e}^0\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^0\left(\mathbf{d}_{\rm new}\right)$
Rank-1 tensor in **i1** direction | 2 | $\mathbf{e}^1\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^1\left(\mathbf{d}_{\rm new}\right)$
Rank-1 tensor in **i2** direction | 3 | $\mathbf{e}^2\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^2\left(\mathbf{d}_{\rm new}\right)$
Rank-2 tensor in **i0-i0** direction | 4 | $\left[\mathbf{e}^0\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^0\left(\mathbf{d}_{\rm new}\right)\right]^2 = 1$
Rank-2 tensor in **i0-i1** direction | 5 | $\left[\mathbf{e}^0\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^0\left(\mathbf{d}_{\rm new}\right)\right]\left[\mathbf{e}^1\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^1\left(\mathbf{d}_{\rm new}\right)\right]$
Rank-2 tensor in **i0-i2** direction | 6 | $\left[\mathbf{e}^0\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^0\left(\mathbf{d}_{\rm new}\right)\right]\left[\mathbf{e}^2\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^2\left(\mathbf{d}_{\rm new}\right)\right]$
Rank-2 tensor in **i1-i1** direction | 7 | $\left[\mathbf{e}^1\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^1\left(\mathbf{d}_{\rm new}\right)\right]^2 = 1$
Rank-2 tensor in **i1-i2** direction | 8 | $\left[\mathbf{e}^1\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^1\left(\mathbf{d}_{\rm new}\right)\right]\left[\mathbf{e}^2\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^2\left(\mathbf{d}_{\rm new}\right)\right]$
Rank-2 tensor in **i2-i2** direction | 9 | $\left[\mathbf{e}^2\left(\mathbf{d}_{\rm gz}\right) \cdot \mathbf{e}^2\left(\mathbf{d}_{\rm new}\right)\right]^2 = 1$

In [Step 2.b.i](#parity_digits) below, we set the middle column of the table for each gridfunction, and [Step 2.b.ii](#dot) sets the rightmost column.

<a id='parity_digits'></a>

### Step 2.b.i: Set parity type for each gridfunction, based on the digits at the end of its name \[Back to [top](#toc)\]
$$\label{parity_digits}$$

For example, if the gridfunction name ends with "01", then (based on the table above) the `set_parity_types()` function below will set the parity_type of that gridfunction to 5. We can be assured this is a robust algorithm, because `gri.register_gridfunctions()` in [grid.py](../edit/grid.py) will throw an error if a gridfunction base name ends in an integer. 

After each parity type is found, we store the parity type of each gridfunction to `const int8_t arrays` `evol_gf_parity` and `aux_gf_parity`, appended to the end of `CurviBoundaryConditions/gridfunction\_defines.h`.


```python
# Step 2.B.1: set the parity conditions on all gridfunctions in gf_list,
#       based on how many digits are at the end of their names
def set_parity_types(gf_list):
    parity_type = []
    for i in range(len(gf_list)):
        varname = gf_list[i]
        parity_type__orig_len = len(parity_type)
        if  len(varname)>2:
            if   varname[-2] == "0" and varname[-1] == "0": # In Python, a[-1] points to the last
                                                            # element of a list; a[-2] the
                                                            # second-to-last element, etc.
                parity_type.append(4)
            elif varname[-2] == "0" and varname[-1] == "1":
                parity_type.append(5)
            elif varname[-2] == "0" and varname[-1] == "2":
                parity_type.append(6)
            elif varname[-2] == "1" and varname[-1] == "1":
                parity_type.append(7)
            elif varname[-2] == "1" and varname[-1] == "2":
                parity_type.append(8)
            elif varname[-2] == "2" and varname[-1] == "2":
                parity_type.append(9)
        if len(varname)>1 and len(parity_type) == parity_type__orig_len:
            if   varname[-1] == "0":
                parity_type.append(1)
            elif varname[-1] == "1":
                parity_type.append(2)
            elif varname[-1] == "2":
                parity_type.append(3)
        if varname[len(varname)-1].isdigit() == False:
            parity_type.append(0)

        if len(parity_type) == parity_type__orig_len:
            print("Error: Could not figure out parity type for evolved variable: "+varname)
            exit(1)
    return parity_type

evol_parity_type = set_parity_types(evolved_variables_list)
aux_parity_type = set_parity_types(auxiliary_variables_list)
auxevol_parity_type = set_parity_types(auxevol_variables_list)

with open("CurviBoundaryConditions/gridfunction_defines.h", "a") as file:
    file.write("\n\n/* PARITY TYPES FOR ALL GRIDFUNCTIONS.\n")
    file.write("   SEE \"Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide.ipynb\" FOR DEFINITIONS. */\n")

    if len(evolved_variables_list) > 0:
        file.write("const int8_t evol_gf_parity[" + str(len(evolved_variables_list)) + "] = { ")
        for i in range(len(evolved_variables_list) - 1):
            file.write(str(evol_parity_type[i]) + ", ")
        file.write(str(evol_parity_type[len(evolved_variables_list) - 1]) + " };\n")

    if len(auxiliary_variables_list) > 0:
        file.write("const int8_t aux_gf_parity[" + str(len(auxiliary_variables_list)) + "] = { ")
        for i in range(len(auxiliary_variables_list) - 1):
            file.write(str(aux_parity_type[i]) + ", ")
        file.write(str(aux_parity_type[len(auxiliary_variables_list) - 1]) + " };\n")

    if len(auxevol_variables_list) > 0:
        file.write("const int8_t auxevol_gf_parity[" + str(len(auxevol_variables_list)) + "] = { ")
        for i in range(len(auxevol_variables_list) - 1):
            file.write(str(auxevol_parity_type[i]) + ", ")
        file.write(str(auxevol_parity_type[len(auxevol_variables_list) - 1]) + " };\n")


for i in range(len(evolved_variables_list)):
    print("Evolved gridfunction \""+evolved_variables_list[i]+"\" has parity type "+str(evol_parity_type[i])+".")
for i in range(len(auxiliary_variables_list)):
    print("Auxiliary gridfunction \""+auxiliary_variables_list[i]+"\" has parity type "+str(aux_parity_type[i])+".")
for i in range(len(auxevol_variables_list)):
    print("AuxEvol gridfunction \""+auxevol_variables_list[i]+"\" has parity type "+str(auxevol_parity_type[i])+".")

print("Output from CurviBoundaryConditions/gridfunction_defines.h (notice the parity types appended to the bottom):")
!cat CurviBoundaryConditions/gridfunction_defines.h
```

    Auxiliary gridfunction "rankoneU0" has parity type 1.
    Auxiliary gridfunction "rankoneU1" has parity type 2.
    Auxiliary gridfunction "rankoneU2" has parity type 3.
    Auxiliary gridfunction "ranktwosymmDD00" has parity type 4.
    Auxiliary gridfunction "ranktwosymmDD01" has parity type 5.
    Auxiliary gridfunction "ranktwosymmDD02" has parity type 6.
    Auxiliary gridfunction "ranktwosymmDD11" has parity type 7.
    Auxiliary gridfunction "ranktwosymmDD12" has parity type 8.
    Auxiliary gridfunction "ranktwosymmDD22" has parity type 9.
    Auxiliary gridfunction "rankzero" has parity type 0.
    Output from CurviBoundaryConditions/gridfunction_defines.h (notice the parity types appended to the bottom):
    /* This file is automatically generated by NRPy+. Do not edit. */
    
    /* EVOLVED VARIABLES: */
    #define NUM_EVOL_GFS 0
    
    
     /* AUXILIARY VARIABLES: */
    #define NUM_AUX_GFS 10
    #define RANKONEU0GF	0
    #define RANKONEU1GF	1
    #define RANKONEU2GF	2
    #define RANKTWOSYMMDD00GF	3
    #define RANKTWOSYMMDD01GF	4
    #define RANKTWOSYMMDD02GF	5
    #define RANKTWOSYMMDD11GF	6
    #define RANKTWOSYMMDD12GF	7
    #define RANKTWOSYMMDD22GF	8
    #define RANKZEROGF	9
    
    
     /* AUXEVOL VARIABLES: */
    #define NUM_AUXEVOL_GFS 0
    
    
    /* PARITY TYPES FOR ALL GRIDFUNCTIONS.
       SEE "Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide.ipynb" FOR DEFINITIONS. */
    const int8_t aux_gf_parity[10] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 0 };


<a id='dot'></a>

###  Step 2.b.ii: Set up unit-vector dot products (=parity) for each of the 10 parity condition types \[Back to [top](#toc)\]
$$\label{dot}$$

First we fill in the parity condition arrays. These take as input $(x_0,x_1,x_2)_{\rm in}$ and $(x_0,x_1,x_2)_{\rm IB}$, and output the necessary dot product(s) for each parity type. To wit (as described above), there are 10 parity types for BSSN evolved variables, which include tensors up to and including rank-2:


```python
# Step 2.B.2: Set up unit-vector dot products (=parity) for each of the 10 parity condition types
parity = ixp.zerorank1(DIM=10)
UnitVectors_inner = ixp.zerorank2()
xx0_inbounds,xx1_inbounds,xx2_inbounds = sp.symbols("xx0_inbounds xx1_inbounds xx2_inbounds", real=True)
for i in range(3):
    for j in range(3):
        UnitVectors_inner[i][j] = rfm.UnitVectors[i][j].subs(rfm.xx[0],xx0_inbounds).subs(rfm.xx[1],xx1_inbounds).subs(rfm.xx[2],xx2_inbounds)
# Type 0: scalar
parity[0] = sp.sympify(1)
# Type 1: i0-direction vector or one-form
# Type 2: i1-direction vector or one-form
# Type 3: i2-direction vector or one-form
for i in range(3):
    for Type in range(1,4):
        parity[Type] += rfm.UnitVectors[Type-1][i]*UnitVectors_inner[Type-1][i]
# Type 4: i0i0-direction rank-2 tensor
# parity[4] = parity[1]*parity[1]
# Type 5: i0i1-direction rank-2 tensor
# Type 6: i0i2-direction rank-2 tensor
# Type 7: i1i1-direction rank-2 tensor
# Type 8: i1i2-direction rank-2 tensor
# Type 9: i2i2-direction rank-2 tensor
count = 4
for i in range(3):
    for j in range(i,3):
        parity[count] = parity[i+1]*parity[j+1]
        count = count + 1

lhs_strings = []
for i in range(10):
    lhs_strings.append("parity["+str(i)+"]")
outputC(parity,lhs_strings, "CurviBoundaryConditions/set_parity_conditions.h")

print("\n\nExample: parity type 1's dot product is given by: \n"+lhs_strings[1]+" = "+str(parity[1]))
```

    Wrote to file "CurviBoundaryConditions/set_parity_conditions.h"
    
    
    Example: parity type 1's dot product is given by: 
    parity[1] = sin(xx1)*sin(xx1_inbounds)*sin(xx2)*sin(xx2_inbounds) + sin(xx1)*sin(xx1_inbounds)*cos(xx2)*cos(xx2_inbounds) + cos(xx1)*cos(xx1_inbounds)


<a id='gzmap'></a>

##  Step 2.c: Implement a modified version of the scalar wave in curvilinear coordinates boundary condition ghost zone mapping routine \[Back to [top](#toc)\]
$$\label{gzmap}$$

We implement a modified version of the scalar wave in curvilinear coordinates boundary condition ghost zone mapping routine, so that it also defines the parity conditions

`set_up_bc_gz_map_and_parity_conditions()` below implements a two-step algorithm:

1. Find locations to where outer ghost zone gridpoints map (this is identical to the algorithm used in the [scalar wave in curvilinear coordinates tutorial](Tutorial-Start_to_Finish-ScalarWaveCurvilinear.ipynb), which is described above in terms of $\mathbf{d}_{\rm gz}$ and $\mathbf{d}_{\rm new}$)
    1. As described in the [scalar wave in curvilinear coordinates tutorial](Tutorial-Start_to_Finish-ScalarWaveCurvilinear.ipynb), this requires first mapping from the curvilinear coordinate gridpoints $(x_0,x_1,x_2)$ in the outer ghost zones to the corresponding Cartesian grid points $(x,y,z)$ in the grid interior or outer boundary (handled by xxminmax, xxCart, and Cart_to_xx defined below).
    1. The interior gridpoint to which each ghost zone maps will be stored in the *ghostzone_map* data structure:
```C
typedef struct ghostzone_map {
   short i0,i1,i2;
} gz_map;
```
1. At each ghost zone gridpoint, find and store the correct parity condition type for each gridfunction up to rank 2 (the highest rank in the BSSN RHSs). As described above, there are a total of 10 possible parity conditions, which each have a unique behavior at a given ghost zone. Thus:
    1. at each ghost zone, we store all 10 parity conditions, in the data structure:
```C
typedef struct parity_conditions {
  int8_t parity[10];
} parity_condition;
```
    1. In the above data structure, `parity[i]` can only take a value of $\pm 1$, which is why we store it in the smallest C integer data type, `int8_t` (a one-byte, signed integer).
    
In the below code blocks:

1. we first output basic curvilinear$\leftrightarrow$Cartesian grid mapping C code, which is required by Step 1 of the above algorithm, and
1. then we implement Steps 1 and 2 of the above algorithm.


```python
# Step 2.C: Modified version of the scalar wave in curvilinear coordinates boundary 
#           condition ghost zone mapping routine, so that it also defines the parity 
#           conditions.

# First output code needed for mapping from any given curvilinear coordinate gridpoint 
#  to the Cartesian coordinate in the grid interior (xxCart), and then find the 
#  corresponding gridpoint index in the grid interior (Cart_to_xx; xxminmax).
# Generic coordinate NRPy+ file output, Part 1: output the conversion from (x0,x1,x2) to Cartesian (x,y,z)
outputC([rfm.xxCart[0],rfm.xxCart[1],rfm.xxCart[2]],["xCart[0]","xCart[1]","xCart[2]"],
        "CurviBoundaryConditions/xxCart.h")
# Generic coordinate NRPy+ file output, Part 2: output the coordinate bounds xxmin[] and xxmax[]:
with open("CurviBoundaryConditions/xxminmax.h", "w") as file:
    file.write("const REAL xxmin[3] = {"+str(rfm.xxmin[0])+","+str(rfm.xxmin[1])+","+str(rfm.xxmin[2])+"};\n")
    file.write("const REAL xxmax[3] = {"+str(rfm.xxmax[0])+","+str(rfm.xxmax[1])+","+str(rfm.xxmax[2])+"};\n")
# Generic coordinate NRPy+ file output, Part 3: output the conversion from Cartesian (x,y,z) to interior/OB (x0,x1,x2)
outputC([rfm.Cart_to_xx[0],rfm.Cart_to_xx[1],rfm.Cart_to_xx[2]],
        ["Cart_to_xx0_inbounds","Cart_to_xx1_inbounds","Cart_to_xx2_inbounds"],
        "CurviBoundaryConditions/Cart_to_xx.h")
```

    Wrote to file "CurviBoundaryConditions/xxCart.h"
    Wrote to file "CurviBoundaryConditions/Cart_to_xx.h"



```python
%%writefile CurviBoundaryConditions/curvilinear_parity_and_outer_boundary_conditions.h

// First we define the struct that will be used to store the 10 parity conditions at all gridpoints:
// We store the 10 parity conditions in a struct consisting of 10 integers, one for each condition.
// Note that these conditions can only take one of two values: +1 or -1.
typedef struct parity_conditions {
  int8_t parity[10];
} parity_condition;

typedef struct ghostzone_map {
  short i0,i1,i2;
} gz_map;

void set_bc_parity_conditions(REAL parity[10], const REAL xx0,const REAL xx1,const REAL xx2, 
                              const REAL xx0_inbounds,const REAL xx1_inbounds,const REAL xx2_inbounds) {
    #include "set_parity_conditions.h"
}

void set_up_bc_gz_map_and_parity_conditions(const int Nxx_plus_2NGHOSTS[3], REAL *xx[3], 
                                            const REAL dxx[3], const REAL xxmin[3], const REAL xxmax[3], 
                                            gz_map *bc_gz_map, parity_condition *bc_parity_conditions) {
  LOOP_REGION(0,Nxx_plus_2NGHOSTS[0],0,Nxx_plus_2NGHOSTS[1],0,Nxx_plus_2NGHOSTS[2]) {
    // First find Cartesian coordinate corresponding to (x_0,x_1,x_2)=(xx[0][i0],xx[1][i1],xx[2][i2]):
    REAL xCart[3];
    xxCart(xx, i0,i1,i2, xCart);
    REAL Cartx = xCart[0];
    REAL Carty = xCart[1];
    REAL Cartz = xCart[2];
    
    // Next find the (i0_inbounds,i1_inbounds,i2_inbounds) corresponding to the above Cartesian coordinate.
    //   If (i0_inbounds,i1_inbounds,i2_inbounds) is in a ghost zone, then it must equal (i0,i1,i2), and
    //      the point is an outer boundary point.
    //   Otherwise (i0_inbounds,i1_inbounds,i2_inbounds) is in the grid interior, and data at (i0,i1,i2)
    //      must be replaced with data at (i0_inbounds,i1_inbounds,i2_inbounds), but multiplied by the
    //      appropriate parity condition (+/- 1).
    REAL Cart_to_xx0_inbounds,Cart_to_xx1_inbounds,Cart_to_xx2_inbounds;
#include "Cart_to_xx.h"
    int i0_inbounds = (int)( (Cart_to_xx0_inbounds - xxmin[0] - (1.0/2.0)*dxx[0] + ((REAL)NGHOSTS)*dxx[0])/dxx[0] + 0.5 ); 
    int i1_inbounds = (int)( (Cart_to_xx1_inbounds - xxmin[1] - (1.0/2.0)*dxx[1] + ((REAL)NGHOSTS)*dxx[1])/dxx[1] + 0.5 );
    int i2_inbounds = (int)( (Cart_to_xx2_inbounds - xxmin[2] - (1.0/2.0)*dxx[2] + ((REAL)NGHOSTS)*dxx[2])/dxx[2] + 0.5 );

    REAL xCart_orig[3]; for(int ii=0;ii<3;ii++) xCart_orig[ii] = xCart[ii];
    xxCart(xx, i0_inbounds,i1_inbounds,i2_inbounds, xCart);

#define EPS_ABS 1e-8
    if(fabs( (double)(xCart_orig[0] - xCart[0]) ) > EPS_ABS ||
       fabs( (double)(xCart_orig[1] - xCart[1]) ) > EPS_ABS ||
       fabs( (double)(xCart_orig[2] - xCart[2]) ) > EPS_ABS) {
      printf("Error. Cartesian disagreement: ( %.15e %.15e %.15e ) != ( %.15e %.15e %.15e )\n",
             (double)xCart_orig[0],(double)xCart_orig[1],(double)xCart_orig[2],
             (double)xCart[0],(double)xCart[1],(double)xCart[2]);
      exit(1);
    }

    if(i0_inbounds-i0 == 0 && i1_inbounds-i1 == 0 && i2_inbounds-i2 == 0) {
      bc_gz_map[IDX3(i0,i1,i2)].i0=-1;
      bc_gz_map[IDX3(i0,i1,i2)].i1=-1;
      bc_gz_map[IDX3(i0,i1,i2)].i2=-1;
      for(int which_parity=0; which_parity<10; which_parity++) {
        bc_parity_conditions[IDX3(i0,i1,i2)].parity[which_parity] = 1;
      }
    } else {
      bc_gz_map[IDX3(i0,i1,i2)].i0=i0_inbounds;
      bc_gz_map[IDX3(i0,i1,i2)].i1=i1_inbounds;
      bc_gz_map[IDX3(i0,i1,i2)].i2=i2_inbounds;
      const REAL xx0 = xx[0][i0];
      const REAL xx1 = xx[1][i1];
      const REAL xx2 = xx[2][i2];
      const REAL xx0_inbounds = xx[0][i0_inbounds];
      const REAL xx1_inbounds = xx[1][i1_inbounds];
      const REAL xx2_inbounds = xx[2][i2_inbounds];
      REAL REAL_parity_array[10];
      set_bc_parity_conditions(REAL_parity_array,  xx0,xx1,xx2, xx0_inbounds,xx1_inbounds,xx2_inbounds);
      for(int whichparity=0;whichparity<10;whichparity++) {
          //printf("Good? Parity %d evaluated to %e\n",whichparity,(double)REAL_parity_array[whichparity]);
          // Perform sanity check on parity array output: should be +1 or -1 to within 8 significant digits:
          if( (REAL_parity_array[whichparity]  > 0 && fabs(REAL_parity_array[whichparity] - (+1)) > 1e-8) ||
              (REAL_parity_array[whichparity] <= 0 && fabs(REAL_parity_array[whichparity] - (-1)) > 1e-8) ) {
              printf("Error. Parity evaluated to %e , which is not within 8 significant digits of +1 or -1.",REAL_parity_array[whichparity]);
              exit(1);
          }
          if(REAL_parity_array[whichparity] < 0.0) bc_parity_conditions[IDX3(i0,i1,i2)].parity[whichparity] = -1;
          if(REAL_parity_array[whichparity] > 0.0) bc_parity_conditions[IDX3(i0,i1,i2)].parity[whichparity] = +1;
      }
    }
  }
}
```

    Overwriting CurviBoundaryConditions/curvilinear_parity_and_outer_boundary_conditions.h


<a id='bpcs'></a>

## Step 2.d: Construct the C loops that apply boundary and parity conditions \[Back to [top](#toc)\]
$$\label{bpcs}$$

Now that `bc_gz_map` and `bc_parity_conditions` have been set at all gridpoints for all gridfunctions, construct the C loops that apply these boundary and parity conditions.

The algorithm below loops over all gridfunctions at all boundary gridpoints, starting from the innermost ghost zone in all 6 faces of the cube and working its way out (to avoid data dependency issues).

1. At the innermost ghost zone, only gridpoints corresponding to outer boundaries are updated first, followed by another pass in which only gridpoints corresponding to inner (parity) boundaries are updated. This ensures that points on an outer boundary that actually maps to another point on the outer boundary are properly updated.
1. The above process is then repeated at the next-to-innermost ghost zone, and then the next-to-next, etc., until all ghost zones have been filled. This ordering must be obeyed, or data will be filled in with unknown values.


```python
%%writefile -a CurviBoundaryConditions/curvilinear_parity_and_outer_boundary_conditions.h

// Part P6: Declare boundary condition OB_UPDATE macro,
//          which updates a single face of the 3D grid cube
//          with
//          1. quadratic polynomial extrapolation, if the face
//             corresponds to an outer boundary, or
//          2. parity condition, if the face maps to a point
//             in the grid interior.
const int MAXFACE = -1;
const int NUL     = +0;
const int MINFACE = +1;

#define OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \
  LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) {                  \
    const int idx3 = IDX3(i0,i1,i2);                                    \
    if(bc_gz_map[idx3].i0 == -1) {                                      \
      gfs[IDX4(which_gf,i0,i1,i2)] =                                    \
        +3.0*gfs[IDX4(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)]    \
        -3.0*gfs[IDX4(which_gf,i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)]    \
        +1.0*gfs[IDX4(which_gf,i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)];   \
    }                                                                   \
  }

#define OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \
  LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) {                  \
    const int idx3 = IDX3(i0,i1,i2);                                    \
    if(bc_gz_map[idx3].i0 != -1) {                                      \
      gfs[IDX4(which_gf,i0,i1,i2)] =                                    \
        ( (REAL)bc_parity_conditions[idx3].parity[gfs_parity[which_gf]] )* \
        gfs[IDX4(which_gf,                                              \
                 bc_gz_map[idx3].i0,                                    \
                 bc_gz_map[idx3].i1,                                    \
                 bc_gz_map[idx3].i2)];                                  \
    }                                                                   \
  }

// Part P7: Boundary condition driver routine: Apply BCs to all six
//          boundary faces of the cube, filling in the innermost
//          ghost zone first, and moving outward.
void apply_bcs(const int Nxx[3],const int Nxx_plus_2NGHOSTS[3],
               gz_map *bc_gz_map,parity_condition *bc_parity_conditions,int num_gfs,const int8_t *gfs_parity, REAL *gfs) {
#pragma omp parallel for
  for(int which_gf=0;which_gf<num_gfs;which_gf++) {
    int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };
    int imax[3] = { Nxx_plus_2NGHOSTS[0]-NGHOSTS, Nxx_plus_2NGHOSTS[1]-NGHOSTS, Nxx_plus_2NGHOSTS[2]-NGHOSTS };
    for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {

      // First apply OUTER boundary conditions,
      //   in case an INNER (parity) boundary point
      //   needs data at the outer boundary:
      // After updating each face, adjust imin[] and imax[] 
      //   to reflect the newly-updated face extents.
      OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;
      OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;

      OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;
      OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;

      OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE); imin[2]--;
      OB_UPDATE_OUTER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE); imax[2]++;

      // Then apply INNER (parity) boundary conditions:
      for(int ii=0;ii<3;ii++) {imin[ii]++; imax[ii]--;}
      // After updating each face, adjust imin[] and imax[] 
      //   to reflect the newly-updated face extents.
      OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;
      OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;

      OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;
      OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;

      OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE); imin[2]--;
      OB_UPDATE_INNER(which_gf, bc_gz_map,bc_parity_conditions, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE); imax[2]++;

    }
  }
}
```

    Appending to CurviBoundaryConditions/curvilinear_parity_and_outer_boundary_conditions.h


<a id='validate'></a>

# Step 3: Set up test data for Curvilinear Boundary Conditions code validation \[Back to [top](#toc)\]
$$\label{validate}$$

We will validate this curvilinear boundary condition module by comparing its results with the original (trusted) SENR code, as follows:

* **Discrete data test**:
    1. Fill all 10 gridfunctions at each gridpoint with the unique gridpoint integer index `IDX3(i0,i1,i2)`
    1. Apply curvilinear boundary conditions
    1. Compare output data at all gridpoints with those from the original SENR code. Agreement should be perfect.
    
Another (future, to-be-implemented) test, which will enable us to validate coordinate systems that do not exist within the original SENR code, is described below:

* **Smooth data test** (TODO):
    1. Fill all 10 gridfunctions with data that are smooth in the Cartesian basis.
    1. Apply Jacobian transformation to all data points, to convert to curvilinear basis
    1. Apply curvilinear boundary conditions
    1. Apply Jacobian transformation to all data points, to convert back to Cartesian basis
    1. Compute difference between original Cartesian data and transformed data. Difference should be zero (to within roundoff) at all points except those that are influenced by outer boundary conditions.


```python
%%writefile CurviBoundaryConditions/CurviBC_Discrete_Test.h

void Discrete_initial_data(const int Nxx_plus_2NGHOSTS[3],REAL *in_gfs) {
    #pragma omp parallel for
    for(int i=0;i<Nxx_plus_2NGHOSTS[0]*Nxx_plus_2NGHOSTS[1]*Nxx_plus_2NGHOSTS[2]*NUM_AUX_GFS;i++) {
        in_gfs[i] = (REAL)i;
    }
}
```

    Overwriting CurviBoundaryConditions/CurviBC_Discrete_Test.h


<a id='mainc'></a>

# Step 4: `CurviBC_Playground.c`: The Main C Code \[Back to [top](#toc)\]
$$\label{mainc}$$



```python
# Part P0: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER
with open("CurviBoundaryConditions/NGHOSTS.h", "w") as file:
    file.write("// Part P0: Set the number of ghost zones, from NRPy+'s FD_CENTDERIVS_ORDER\n")
    # Upwinding in BSSN requires that NGHOSTS = FD_CENTDERIVS_ORDER/2 + 1 <- Notice the +1.
    file.write("#define NGHOSTS "+str(int(par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")/2)+1)+"\n")
```


```python
%%writefile CurviBoundaryConditions/CurviBC_Playground.c

// Step P1: Import needed header files
#include "NGHOSTS.h" // A NRPy+-generated file, which is set based on FD_CENTDERIVS_ORDER.
#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "time.h"
#include "stdint.h" // Needed for Windows GCC 6.x compatibility

// Step P2: Add needed #define's to set data type, the IDX4() macro, and the gridfunctions
// Step P2a: set REAL=double, so that all floating point numbers are stored to at least ~16 significant digits.
#define REAL double

// Step P3: Set free parameters
// Step P3a: Free parameters for the numerical grid

// Spherical coordinates parameter:
// Set RMAX, the scale of the numerical domain, 
//  and the max radius in Spherical coordinates.
#define RMAX 7.5
// SinhSpherical coordinates parameters:
const REAL AMPL    = RMAX;
const REAL SINHW   = 0.2;
// Cylindrical coordinates parameters:
const REAL ZMIN   = -RMAX;
const REAL ZMAX   =  RMAX;
const REAL RHOMAX =  RMAX;
// Cartesian coordinates parameters:
const REAL xmin = -RMAX, xmax = RMAX;
const REAL ymin = -RMAX, ymax = RMAX;
const REAL zmin = -RMAX, zmax = RMAX;
// SymTP coordinates parameters:
const REAL bScale =   0.5;
const REAL AMAX   =  RMAX;

// Step P4: Declare the IDX4(gf,i,j,k) macro, which enables us to store 4-dimensions of
//          data in a 1D array. In this case, consecutive values of "i" 
//          (all other indices held to a fixed value) are consecutive in memory, where 
//          consecutive values of "j" (fixing all other indices) are separated by 
//          Nxx_plus_2NGHOSTS[0] elements in memory. Similarly, consecutive values of
//          "k" are separated by Nxx_plus_2NGHOSTS[0]*Nxx_plus_2NGHOSTS[1] in memory, etc.
#define IDX4(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS[0] * ( (j) + Nxx_plus_2NGHOSTS[1] * ( (k) + Nxx_plus_2NGHOSTS[2] * (g) ) ) )
#define IDX3(i,j,k) ( (i) + Nxx_plus_2NGHOSTS[0] * ( (j) + Nxx_plus_2NGHOSTS[1] * (k) ) )
// Assuming idx = IDX3(i,j,k). Much faster if idx can be reused over and over:
#define IDX4pt(g,idx)   ( (idx) + (Nxx_plus_2NGHOSTS[0]*Nxx_plus_2NGHOSTS[1]*Nxx_plus_2NGHOSTS[2]) * (g) )

// Step P5: Set #define's for BSSN gridfunctions. C code generated above
#include "gridfunction_defines.h"

#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \
  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)

void xxCart(REAL *xx[3],const int i0,const int i1,const int i2, REAL xCart[3]) {
    REAL xx0 = xx[0][i0];
    REAL xx1 = xx[1][i1];
    REAL xx2 = xx[2][i2];
#include "xxCart.h"
}

// Step P6: Include basic functions needed to impose curvilinear
//          parity and boundary conditions.
#include "curvilinear_parity_and_outer_boundary_conditions.h"

// Step P7: Declare the function for the exact solution. time==0 corresponds to the initial data.
#include "CurviBC_Discrete_Test.h"

// main() function:
// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates
// Step 1: Write test data to gridfunctions
// Step 2: Overwrite all data in ghost zones with NaNs
// Step 3: Output relative error between numerical and exact solution.
// Step 4: Print gridfunction data after curvilinear boundary conditions have been applied:
// Step 5: Free all allocated memory
int main(int argc, const char *argv[]) {

    // Step 0a: Read command-line input, error out if nonconformant
    if(argc != 5 || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < NGHOSTS) {
        fprintf(stderr,"Error: Expected one command-line argument: ./CurviBC_Playground Nx0 Nx1 Nx2 [test type: Smooth or Discrete],\n");
        fprintf(stderr,"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\n");
        fprintf(stderr,"Nx[] MUST BE larger than NGHOSTS (= %d)\n",NGHOSTS);
        exit(1);
    }
    // Step 0b: Set test type to Smooth or Discrete
    char test_type[100];
    snprintf(test_type,100,"%s",argv[4]);
    if(strncmp("Smooth",test_type,100)!=0 && strncmp("Discrete",test_type,100)!=0) {
        fprintf(stderr,"Error: test type = %s not supported. Choose Smooth or Discrete.\n",test_type);
        exit(1);
    }
    // Step 0b: Set up numerical grid structure, first in space...
    if(atoi(argv[1])%2 != 0 || atoi(argv[2])%2 != 0 || atoi(argv[3])%2 != 0) {
        fprintf(stderr,"Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\n");
        fprintf(stderr,"       For example, in case of angular directions, proper symmetry zones will not exist.\n");
        exit(1);
    }
    const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };
    const int Nxx_plus_2NGHOSTS[3] = { Nxx[0]+2*NGHOSTS, Nxx[1]+2*NGHOSTS, Nxx[2]+2*NGHOSTS };
    const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS[0]*Nxx_plus_2NGHOSTS[1]*Nxx_plus_2NGHOSTS[2];
    #include "xxminmax.h"

    // Step 0c: Allocate memory for gridfunctions
    REAL *test_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUX_GFS * Nxx_plus_2NGHOSTS_tot);

    // Step 0d: Set up space and time coordinates
    // Step 0d.i: Set \Delta x^i on uniform grids.
    REAL dxx[3];
    for(int i=0;i<3;i++) dxx[i] = (xxmax[i] - xxmin[i]) / ((REAL)Nxx[i]);

    // Step 0d.ii: Set up uniform coordinate grids
    REAL *xx[3];
    for(int i=0;i<3;i++) {
        xx[i] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS[i]);
        for(int j=0;j<Nxx_plus_2NGHOSTS[i];j++) {
            xx[i][j] = xxmin[i] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*dxx[i]; // Cell-centered grid.
        }
    }

    // Step 0e: Find ghostzone mappings and parities:
    gz_map *bc_gz_map = (gz_map *)malloc(sizeof(gz_map)*Nxx_plus_2NGHOSTS_tot);
    parity_condition *bc_parity_conditions = (parity_condition *)malloc(sizeof(parity_condition)*Nxx_plus_2NGHOSTS_tot);
    set_up_bc_gz_map_and_parity_conditions(Nxx_plus_2NGHOSTS,xx,dxx,xxmin,xxmax,  bc_gz_map, bc_parity_conditions);

    // Step 1: Write test data to gridfunctions
    if(strncmp("Discrete",test_type,100)==0) {
        Discrete_initial_data(Nxx_plus_2NGHOSTS, test_gfs);
    } else {
        fprintf(stderr,"Sorry, curvilinear boundary conditions test = %s not yet supported. Feel free to contribute!\n",test_type);
    }
    
    // Step 2: Overwrite all data in ghost zones with NaNs
    LOOP_REGION(0,Nxx_plus_2NGHOSTS[0], 0,Nxx_plus_2NGHOSTS[1], 0,Nxx_plus_2NGHOSTS[2]) {
        for(int gf=0;gf<NUM_AUX_GFS;gf++) {
            const int idx4 = IDX4(gf,i0,i1,i2);
            if(i0 < NGHOSTS || i0 >= Nxx_plus_2NGHOSTS[0]-NGHOSTS) test_gfs[idx4] = +(0.0 / 0.0);
            if(i1 < NGHOSTS || i1 >= Nxx_plus_2NGHOSTS[1]-NGHOSTS) test_gfs[idx4] = +(0.0 / 0.0);
            if(i2 < NGHOSTS || i2 >= Nxx_plus_2NGHOSTS[2]-NGHOSTS) test_gfs[idx4] = +(0.0 / 0.0);
        }
    }

    // Step 3: Apply curvilinear boundary conditions
    apply_bcs(Nxx, Nxx_plus_2NGHOSTS, bc_gz_map,bc_parity_conditions, NUM_AUX_GFS,aux_gf_parity, test_gfs);
    
    // Step 4: Print gridfunction data after curvilinear boundary conditions have been applied:
    LOOP_REGION(0,Nxx_plus_2NGHOSTS[2], 0,Nxx_plus_2NGHOSTS[1], 0,Nxx_plus_2NGHOSTS[0]) {
        printf("%d %d %d | ",i0,i1,i2);
        for(int gf=0;gf<NUM_AUX_GFS;gf++) {
            const int idx4 = IDX4(gf,i0,i1,i2);
            if(!isnan(test_gfs[idx4])) {
                printf("%d ",(int)test_gfs[idx4]);
            } else {
                fprintf(stderr,"ERROR: found NaN %d %d %d %d %d\n",gf,i0,i1,i2,NUM_AUX_GFS);
                exit(1);
            }
        }
        printf("\n");
    }

    // Step 5: Free all allocated memory
    free(bc_parity_conditions);
    free(bc_gz_map);
    free(test_gfs);
    for(int i=0;i<3;i++) free(xx[i]);
    return 0;
}
```

    Overwriting CurviBoundaryConditions/CurviBC_Playground.c



```python
import cmdline_helper as cmd
cmd.C_compile("CurviBoundaryConditions/CurviBC_Playground.c", "CurviBC_Playground")
cmd.Execute("CurviBC_Playground", "4 4 4 Discrete", "CurviBoundaryConditions/out4x4x4-Spherical-NGHOSTS4oFD.txt")
print("Finished this code cell.")
```

    Compiling executable...
    Executing `gcc -Ofast -fopenmp -march=native -funroll-loops CurviBoundaryConditions/CurviBC_Playground.c -o CurviBC_Playground -lm`...
    Finished executing in 0.821406126022 seconds.
    Finished compilation.
    Executing `taskset -c 0,1,2,3 ./CurviBC_Playground 4 4 4 Discrete`...
    Finished executing in 0.213351964951 seconds.
    Finished this code cell.


<a id='senr_compare'></a>

# Step 5: Validation: Compare with original SENR results \[Back to [top](#toc)\]
$$\label{senr_compare}$$



```python
import sys

import filecmp
if filecmp.cmp('CurviBoundaryConditions/out4x4x4-Spherical-NGHOSTS4oFD.txt',
               'CurviBoundaryConditions/SENRout4x4x4-Spherical_NGHOSTS4oFD.txt') == False:
    print("ERROR: Spherical boundary conditions do not agree!")
    sys.exit(1)
else:
    print("Spherical boundary condition comparison test between this tutorial notebook & trusted original SENR code: PASSED")
```

    Spherical boundary condition comparison test between this tutorial notebook & trusted original SENR code: PASSED



```python
import filecmp
import CurviBoundaryConditions.CurviBoundaryConditions as cbcs
cbcs.Set_up_CurviBoundaryConditions()

print("Running discrete CurviBC_Playground test on Spherical data, using NRPy+ CurviBoundaryConditions.CurviBoundaryConditions module.")
cmd.C_compile("CurviBoundaryConditions/CurviBC_Playground.c", "CurviBC_Playground")
cmd.Execute("CurviBC_Playground", "4 4 4 Discrete", "CurviBoundaryConditions/CBCmodule-out4x4x4-Spherical-NGHOSTS4oFD.txt")

# Then we set the coordinate system for the numerical grid
par.set_parval_from_str("reference_metric::CoordSystem","Cylindrical")
rfm.reference_metric()
cbcs.Set_up_CurviBoundaryConditions()
print("Running discrete CurviBC_Playground test on Cylindrical data, using NRPy+ CurviBoundaryConditions.CurviBoundaryConditions module.")
cmd.C_compile("CurviBoundaryConditions/CurviBC_Playground.c", "CurviBC_Playground")
cmd.Execute("CurviBC_Playground", "4 4 4 Discrete", "CurviBoundaryConditions/CBCmodule-out4x4x4-Cylindrical-NGHOSTS4oFD.txt")
print("Finished this code cell.")

print("\n")

if filecmp.cmp('CurviBoundaryConditions/CBCmodule-out4x4x4-Spherical-NGHOSTS4oFD.txt',
               'CurviBoundaryConditions/SENRout4x4x4-Spherical_NGHOSTS4oFD.txt') == False:
    print("ERROR: Spherical boundary conditions do not agree!")
    sys.exit(1)
else:
    print("Spherical boundary condition comparison test between CBC & trusted original SENR code: PASSED")

if filecmp.cmp('CurviBoundaryConditions/CBCmodule-out4x4x4-Cylindrical-NGHOSTS4oFD.txt',
               'CurviBoundaryConditions/SENRout4x4x4-Cylindrical_NGHOSTS4oFD.txt') == False:
    print("ERROR: Cylindrical boundary conditions do not agree!")
    sys.exit(1)
else:
    print("Cylindrical boundary condition comparison test between CBC & trusted original SENR code: PASSED")
```

    Wrote to file "CurviBoundaryConditions/gridfunction_defines.h"
    Wrote to file "CurviBoundaryConditions/set_parity_conditions.h"
    Wrote to file "CurviBoundaryConditions/xxCart.h"
    Wrote to file "CurviBoundaryConditions/xxminmax.h"
    Wrote to file "CurviBoundaryConditions/Cart_to_xx.h"
    Running discrete CurviBC_Playground test on Spherical data, using NRPy+ CurviBoundaryConditions.CurviBoundaryConditions module.
    Compiling executable...
    Executing `gcc -Ofast -fopenmp -march=native -funroll-loops CurviBoundaryConditions/CurviBC_Playground.c -o CurviBC_Playground -lm`...
    Finished executing in 0.821196079254 seconds.
    Finished compilation.
    Executing `taskset -c 0,1,2,3 ./CurviBC_Playground 4 4 4 Discrete`...
    Finished executing in 0.214704990387 seconds.
    Wrote to file "CurviBoundaryConditions/gridfunction_defines.h"
    Wrote to file "CurviBoundaryConditions/set_parity_conditions.h"
    Wrote to file "CurviBoundaryConditions/xxCart.h"
    Wrote to file "CurviBoundaryConditions/xxminmax.h"
    Wrote to file "CurviBoundaryConditions/Cart_to_xx.h"
    Running discrete CurviBC_Playground test on Cylindrical data, using NRPy+ CurviBoundaryConditions.CurviBoundaryConditions module.
    Compiling executable...
    Executing `gcc -Ofast -fopenmp -march=native -funroll-loops CurviBoundaryConditions/CurviBC_Playground.c -o CurviBC_Playground -lm`...
    Finished executing in 0.818950891495 seconds.
    Finished compilation.
    Executing `taskset -c 0,1,2,3 ./CurviBC_Playground 4 4 4 Discrete`...
    Finished executing in 0.215003013611 seconds.
    Finished this code cell.
    
    
    Spherical boundary condition comparison test between CBC & trusted original SENR code: PASSED
    Cylindrical boundary condition comparison test between CBC & trusted original SENR code: PASSED


<a id='latex_pdf_output'></a>

# Step 6: 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-Start_to_Finish-Curvilinear_BCs.pdf](Tutorial-Start_to_Finish-Curvilinear_BCs.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)


```python
!jupyter nbconvert --to latex --template latex_nrpy_style.tplx Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb
!pdflatex -interaction=batchmode Tutorial-Start_to_Finish-Curvilinear_BCs.tex
!pdflatex -interaction=batchmode Tutorial-Start_to_Finish-Curvilinear_BCs.tex
!pdflatex -interaction=batchmode Tutorial-Start_to_Finish-Curvilinear_BCs.tex
!rm -f Tut*.out Tut*.aux Tut*.log
```

    [NbConvertApp] Converting notebook Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb to latex
    [NbConvertApp] Writing 110822 bytes to Tutorial-Start_to_Finish-Curvilinear_BCs.tex
    This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)
     restricted \write18 enabled.
    entering extended mode
    This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)
     restricted \write18 enabled.
    entering extended mode
    This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)
     restricted \write18 enabled.
    entering extended mode


b42e310f0d463ed5ff6940292c5ab501cdcb9f3d
ipynb
Jupyter Notebook
notebook/Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb
kazewong/nrpytutorial
cc511325f37f01284b2b83584beb2a452556b3fb
BSD-2-Clause
notebook/Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb
kazewong/nrpytutorial
cc511325f37f01284b2b83584beb2a452556b3fb
BSD-2-Clause
notebook/Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb
kazewong/nrpytutorial
cc511325f37f01284b2b83584beb2a452556b3fb
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES